text:
module utilitiesOCIM

# package code goes here
"""
    Iabove, Ibelow = buildIndexMatrices()

Build the matrices of indices above and below.

These are shifted-diagonal matrices that allow for
easy construction of the particle flux divergence `PFD` operator.
"""
function buildIaboveIbelow()
  n = nlon * nlat * (ndepth+1)
  In = speye(n)
  idx = zeros(nlat,nlon,ndepth+1)
  idx[:] = 1:n
  idxabove = idx[:,:,[ndepth+1; 1:ndepth]] # use a periodic upward shift
  idxbelow = idx[:,:,[2:ndepth+1; 1]]      # use a periodic downward shift
  Iabove = In[idxabove[:],:]
  Ibelow = In[idxbelow[:],:]
  # keep only indices of wet boxes
  Iabove = Iabove[iwet,iwet]
  Ibelow = Ibelow[iwet,iwet]
  return Iabove, Ibelow
end
const Iabove, Ibelow = buildIaboveIbelow()

"""
    Iabove = buildIabove(wet3d)

Build the shifted-diagonal sparse matrix of the indices of above neighbours.

Iabove[i,j] = 1 if the box represented by the linear index i
lies directly above the box represented by the linear index j.
"""
function buildIabove(wet3d)
  nlat, nlon, ndepth = size(wet3d)
  iwet = find(wet3d)
  n = nlon * nlat * (ndepth + 1)
  In = speye(n)
  idx = zeros(nlat, nlon, ndepth + 1)
  idx[:] = 1:n
  idx .= idx[:, :, [ndepth + 1 ; 1:ndepth]] # shift
  return In[idx[:], :][iwet, iwet]
end

function buildIabove(wet3d)

end

end # module


text:
# roessler chaotic phase diff
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random, OrdinaryDiffEq, LinearAlgebra

a = 0.2
b = 0.2
c = 5.7
amu = 0.02
ω1 = 1.0 + amu
ω2 = 1.0 - amu
k1 = 0.01
k2 = k1

y_init = [0.1, 0.2, 0., 0.11, 0.19, 0.1]

croesslers = Systems.coupled_roessler(y_init; ω1, ω2, a, b, c, k1, k2)
diffeq = (alg = Vern9(), reltol = 1e-9, abstol = 1e-9)
Δt = 0.01
T = 1000
t = 0:Δt:T
L = length(t)÷10

ks = [0.03, 0.04] # estimated from the Lyapunov exponents plot
bbox2 = Dict(:boxstyle => "round,pad=0.3", :facecolor=>"white", :alpha => 0.8)

fig, axs = subplots(1,2; figsize = (0.66figx, figy))

function phase(x, y)
    f = atan.(y, x)
    for i in 2:length(f)
        if f[i] - f[i-1] < 0
            f[i:end] .+= 2π
        end
    end
    return f
end

for (i, k) in enumerate(ks)
    set_parameter!(croesslers, 6, k)
    set_parameter!(croesslers, 7, k)
    tr = trajectory(croesslers, T; Ttr = 200, Δt, diffeq...)
    x1 = tr[:, 1]
    x2 = tr[:, 2]
    y1 = tr[:, 4]
    y2 = tr[:, 5]
    f1 = phase(x1, x2)
    f2 = phase(y1, y2)
    axs[2].plot(t, f1 .- f2, label = "\$k=$(k)\$")

    if i == 1
        c = COLORS[1]
        axs[1].plot(x1[1:L], x2[1:L]; color = "0.5", linewidth = 2)
        j = 250
        axs[1].quiver(zeros(1), zeros(1), [x1[j+2]], [x2[j+2]],
            color = c,
            scale = 1., scale_units = "xy", zorder = 9,
            width = 0.01
        )
        axs[1].plot([x1[j+1]], [x2[j+1]], marker = "o", zorder = 99, ms = 8,
        ls = "None", mec = c, mfc = :white,
            # scale = 1., scale_units = "xy", zorder = 9,
            # width = 0.02
        )
        θ = atan(x2[j], x1[j])
        r = norm([x1[j], x2[j]])
        a = matplotlib.patches.Wedge((0,0), 0.3r;
        theta1 = 0, theta2 = 180θ/π, color = c, alpha = 0.5)
        axs[1].add_artist(a)
        axs[1].text(0.3x1[j-10], 0.3x2[j-10], "\$\\alpha\$"; color = c)
    end
end

axs[1].set_xlabel("\$x_1\$"; labelpad = -20)
axs[1].set_ylabel("\$x_2\$"; labelpad = -20)
axs[1].set_xticks([-10, 0, 10])
axs[1].set_xticklabels([-10, "", 10])
axs[2].set_xlabel("\$t\$"; labelpad = -20)
axs[2].set_ylabel("\$\\Delta \\alpha\$")
axs[2].set_xticks([0,t[end]])
axs[2].set_xlim(0,t[end])
axs[2].legend()
add_identifiers!(fig)
fig.tight_layout(;pad=0.33)
wsave(plotsdir("9", "coupled_roessler_phasediff"), fig)

ext: jl

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
@def Ax_min begin #(U)

	# minimum longitudinal acceleration for given speed
	Ax_min = Array(Float64,(length(U),1))
	for i in eachindex(U)
		Ax_min[i,1] = AXC[5]*U[i]^3 + AXC[6]*U[i]^2 + AXC[7]*U[i] + AXC[8]
	end
	Ax_min
end


text:
#ubounds-resolve.jl - contains functions which operate on two unum values that
#are meant to go into ubounds.

if options[:devmode]
  macro check_resolve()
    esc(quote
      options[:devmode] && (lower >= upper) && throw(ArgumentError("unum constructors must be well-ordered"))
    end)
  end
else
  macro check_resolve(); :(); end
end

@universal function set_zero_exp!(value::Unum, to_match::Unum)
  if (decode_exp(to_match) >= 0)
    value.esize = z16
    value.fsize = z16
  else
    test_exponent = decode_exp(to_match.esize - o16, z64)
    if (test_exponent - max_fsize(FSS)) > decode_exp(to_match)
      value.esize = to_match.esize
      value.fsize = to_match.fsize
    else
      value.esize = to_match.esize - o16
      value.fsize = test_exponent - decode_exp(to_match)
    end
  end
end

doc"""
  `Unums.resolve_as_utype!(::Unum, ::Unum)`
  checks to see if the two unums, which must both be ulps, can be resolved to a
  single unum.  If that is possible, it returns the new unum.  Otherwise, it
  returns a ubound with the same values.

  Importantly, this function must be allowed to take possession of the two
  numbers.
"""
@universal function resolve_as_utype!(lower::Unum, upper::Unum)
  @ensure_ulp(lower)
  @ensure_ulp(upper)

  (lower == upper) && return lower

  #check to make sure the signs are the same.
  (@signof(lower) == @signof(upper)) || return B(lower, upper)

  #check to see if they're both zero ulps
  if is_zero_ulp(lower) && is_zero_ulp(upper)
    _lexp = decode_exp(lower) - lower.fsize
    _uexp = decode_exp(upper) - upper.fsize
    (_lexp < _uexp) && return upper
    return lower
  end

  #if only one is a zero_ulp.
  if is_zero_ulp(lower) && (decode_exp(upper) < 0)
    #check to see if we can subsume.
    lubu = lub(upper)
    test_exponent = decode_exp(lubu.esize - o16, z64)
    if (lubu.fraction == z64)
      if (test_exponent - max_fsize(FSS) <= decode_exp(lubu))
        lower.esize = lubu.esize - o16
        lower.fsize = test_exponent - decode_exp(lubu)
        return lower
      end
    end
  end

  if is_zero_ulp(upper) && (decode_exp(lower) < 0)
    glbl = glb(lower)
    test_exponent = decode_exp(glbl.esize - o16, z64)
    if (glbl.fraction == z64) && (test_exponent - max_fsize(FSS) <= decode_exp(glbl))
      upper.esize = glbl.esize - o16
      upper.fsize = test_exponent - decode_exp(glbl)
      return upper
    end
  end

  #check to make sure the exponents are the same.
  (decode_exp(lower) == decode_exp(upper)) || return B(lower, upper)

  (lower.fraction == upper.fraction) && return lower

  #discriminate between negative and positive numbers, which have their bounds
  #differently.
  local cfsize::UInt16
  if is_positive(lower)
    #firstly, the inner ubound is privileged and must be longer than the upper.
    cfsize = contract_outer_fsize(upper.fraction, upper.fsize)
    cfsize = max(cfsize, contract_inner_fsize(lower.fraction, lower.fsize))
  else
    cfsize = contract_outer_fsize(lower.fraction, lower.fsize)
    cfsize = max(cfsize, contract_inner_fsize(upper.fraction, upper.fsize))
  end

  same_till_fsize(upper.fraction, lower.fraction, cfsize) || return B(lower, upper)

  #special case where we can't coalesce them into two unums.
  if cfsize == 0
    lower.fsize = 0
    upper.fsize = 0

    if is_negative(upper)
      frac_top!(lower)
    else
      frac_top!(upper)
    end

    return B(lower, upper)
  end
  lower.fsize = cfsize - o16
  frac_trim!(lower, lower.fsize)
  return lower
end

doc"""
  `Unums.resolve_lower(::Unum)`
  `Unums.resolve_lower(::Ubound)`
  Returns a Unum that represents the lower end of a Ubound, or just the Unum itself.
"""
@universal resolve_lower(value::Unum) = value
@universal resolve_lower(value::Ubound) = value.lower

doc"""
  `Unums.resolve_upper(::Unum)`
  `Unums.resolve_upper(::Ubound)`
  Returns a Unum that represents the upper end of a Ubound, or just the Unum itself.
"""
@universal resolve_upper(value::Unum) = value
@universal resolve_upper(value::Ubound) = value.upper

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
    load_network_data(setup::Dict, path::AbstractString, inputs_nw::Dict)

Function for reading input parameters related to the electricity transmission network
"""
#DEV NOTE:  add DC power flow related parameter inputs in a subsequent commit
function load_network_data(setup::Dict, path::AbstractString, inputs_nw::Dict)

    # Network zones inputs and Network topology inputs
    network_var = DataFrame(CSV.File(joinpath(path,"Network.csv"), header=true), copycols=true)

    # Number of zones in the network
    inputs_nw["Z"] = size(findall(s -> (startswith(s, "z")) & (tryparse(Float64, s[2:end]) !== nothing), names(network_var)),1)
    Z = inputs_nw["Z"]
    # Number of lines in the network
    inputs_nw["L"]=size(collect(skipmissing(network_var[!,:Network_Lines])),1)

    # Topology of the network source-sink matrix
    start = findall(s -> s == "z1", names(network_var))[1]
    inputs_nw["pNet_Map"] = Matrix{Float64}(network_var[1:inputs_nw["L"],start:start+inputs_nw["Z"]-1])

    # Transmission capacity of the network (in MW)
    if setup["ParameterScale"] ==1  # Parameter scaling turned on - adjust values of subset of parameter values to GW
        inputs_nw["pTrans_Max"] = convert(Array{Float64}, collect(skipmissing(network_var[!,:Line_Max_Flow_MW])))/ModelScalingFactor  # convert to GW
    else # no scaling
        inputs_nw["pTrans_Max"] = convert(Array{Float64}, collect(skipmissing(network_var[!,:Line_Max_Flow_MW])))
    end

    if setup["Trans_Loss_Segments"] == 1 ##Aaron Schwartz Please check
        # Line percentage Loss - valid for case when modeling losses as a fixed percent of absolute value of power flows
        inputs_nw["pPercent_Loss"] = convert(Array{Float64}, collect(skipmissing(network_var[!,:Line_Loss_Percentage])))
    elseif setup["Trans_Loss_Segments"] >= 2
        # Transmission line voltage (in kV)
        inputs_nw["kV"] = convert(Array{Float64}, collect(skipmissing(network_var[!,:Line_Voltage_kV])))
        # Transmission line resistance (in Ohms) - Used when modeling quadratic transmission losses
        inputs_nw["Ohms"] = convert(Array{Float64}, collect(skipmissing(network_var[!,:Line_Resistance_ohms])))
    end

    # Maximum possible flow after reinforcement for use in linear segments of piecewise approximation
    inputs_nw["pTrans_Max_Possible"] = zeros(Float64, inputs_nw["L"])

    if setup["NetworkExpansion"]==1
        if setup["ParameterScale"] ==1  # Parameter scaling turned on - adjust values of subset of parameter values
            # Read between zone network reinforcement costs per peak MW of capacity added
            inputs_nw["pC_Line_Reinforcement"] = convert(Array{Float64}, collect(skipmissing(network_var[!,:Line_Reinforcement_Cost_per_MWyr])))/ModelScalingFactor # convert to million $/GW/yr with objective function in millions
            # Maximum reinforcement allowed in MW
            #NOTE: values <0 indicate no expansion possible
            inputs_nw["pMax_Line_Reinforcement"] = convert(Array{Float64}, collect(skipmissing(network_var[!,:Line_Max_Reinforcement_MW])))/ModelScalingFactor # convert to GW
        else
            # Read between zone network reinforcement costs per peak MW of capacity added
            inputs_nw["pC_Line_Reinforcement"] = convert(Array{Float64}, collect(skipmissing(network_var[!,:Line_Reinforcement_Cost_per_MWyr])))
            # Maximum reinforcement allowed in MW
            #NOTE: values <0 indicate no expansion possible
            inputs_nw["pMax_Line_Reinforcement"] = convert(Array{Float64}, collect(skipmissing(network_var[!,:Line_Max_Reinforcement_MW])))
        end
        for l in 1:inputs_nw["L"]
            if inputs_nw["pMax_Line_Reinforcement"][l] > 0
                inputs_nw["pTrans_Max_Possible"][l] = inputs_nw["pTrans_Max"][l] + inputs_nw["pMax_Line_Reinforcement"][l]
            else
                inputs_nw["pTrans_Max_Possible"][l] = inputs_nw["pTrans_Max"][l]
            end
        end
    else
        inputs_nw["pTrans_Max_Possible"] = inputs_nw["pTrans_Max"]
    end

    # Transmission line (between zone) loss coefficient (resistance/voltage^2)
    inputs_nw["pTrans_Loss_Coef"] = zeros(Float64, inputs_nw["L"])
    for l in 1:inputs_nw["L"]
        # For cases with only one segment
        if setup["Trans_Loss_Segments"] == 1
            inputs_nw["pTrans_Loss_Coef"][l] = inputs_nw["pPercent_Loss"][l]
        elseif setup["Trans_Loss_Segments"] >= 2
            # If zones are connected, loss coefficient is R/V^2 where R is resistance in Ohms and V is voltage in Volts
            if setup["ParameterScale"] ==1  # Parameter scaling turned on - adjust values of subset of parameter values
                inputs_nw["pTrans_Loss_Coef"][l] = (inputs_nw["Ohms"][l]/10^6)/(inputs_nw["kV"][l]/10^3)^2 *ModelScalingFactor # 1/GW ***
            else
                inputs_nw["pTrans_Loss_Coef"][l] = (inputs_nw["Ohms"][l]/10^6)/(inputs_nw["kV"][l]/10^3)^2 # 1/MW
            end
        end
    end

    ## Sets and indices for transmission losses and expansion
    inputs_nw["TRANS_LOSS_SEGS"] = setup["Trans_Loss_Segments"] # Number of segments used in piecewise linear approximations quadratic loss functions
    inputs_nw["LOSS_LINES"] = findall(inputs_nw["pTrans_Loss_Coef"].!=0) # Lines for which loss coefficients apply (are non-zero);

    if setup["NetworkExpansion"] == 1
        # Network lines and zones that are expandable have non-negative maximum reinforcement inputs
        inputs_nw["EXPANSION_LINES"] = findall(inputs_nw["pMax_Line_Reinforcement"].>=0)
        inputs_nw["NO_EXPANSION_LINES"] = findall(inputs_nw["pMax_Line_Reinforcement"].<0)
    end

    println("Network.csv Successfully Read!")

    return inputs_nw, network_var
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
coordvalue(s, onechar) = reduce((x, c) -> (x << 1) ⊻ (c == onechar), s; init=0)

seatid(s) = 8coordvalue(s[1:7], 'B') + coordvalue(s[8:10], 'R')

part1(seats) = maximum(seats)

part2(seats) = seats[findfirst((seats[2:end] - seats[1:end-1]) .== 2)] + 1

let seats = sort!(map(seatid, readlines("input.txt")))
    @info "Solutions" part1(seats) part2(seats)
end

max_stars_repo_name: dglmoore/advent-of-code

max_issues_repo_name: dglmoore/advent-of-code

max_forks_repo_name: dglmoore/advent-of-code


text:
using FeedbackParticleFilters, Distributions, Random, PDMats

println("Testing LinearDiffusionStateModel.jl:")

@testset "LinearDiffusionStateModel.jl" begin

    A      = [0. 1.; -0.7 -0.3]
    B      = zeros(2,1)
    B[2,1] = sqrt(2.)
    init   = MvNormal(LinearAlgebra.Diagonal(FillArrays.Fill(1., 2)))

    print("  inner constructor")
    model = LinearDiffusionStateModel(A, B, init)
    B2 = zeros(Int, 2, 1)
    B2[2,1] = 2
    print("-")
    @test_throws MethodError LinearDiffusionStateModel(A, B2, init) # wrong element type of B2
    B3 = zeros(3, 1)
    print("-")
    @test_throws DimensionMismatch LinearDiffusionStateModel(A, B3, init) # wrong size of B3
    println("DONE")

    print("  method find_stationary_variance")
    Sigma = find_stationary_variance(A, B*B')
    print("-")
    @test Sigma == [4.746796851882334 -5.543824348282954e-5; -5.543824348282954e-5 3.3192089048776765]
    println("DONE")

    print("  method drift")
    print("-")
    @test drift(model) == A
    println("DONE")

    print("  method diffusion")
    print("-")
    @test diffusion(model) == B
    println("DONE")

    print("  method initial_condition")
    print("-")
    @test initial_condition(model) == model.init
    println("DONE")

    print("  method initialize")
    print("-")
    x0 = initialize(model)
    @test length(x0) == 2
    println("DONE")

    print("  method state_dim")
    print("-")
    @test state_dim(model) == 2
    println("DONE")

    print("  method noise_dim")
    print("-")
    @test noise_dim(model) == 1
    println("DONE")

    print("  method drift_function")
    print("-")
    @test drift_function(model)([1., 3.]) == [3.0, -1.5999999999999999]
    println("DONE")

    print("  method diffusion_function")
    print("-")
    @test diffusion_function(model)([1., 3.]) == hcat([0.0, 1.4142135623730951]...)'
    println("DONE")

    print("  calling instance of DiffusionStateModel")
    x0 = ones(2)
    print("-")
    @test model(x0, 0.01) != x0
    xx0 = hcat(ones(2), zeros(2))
    print("-")
    @test model(xx0, 0.01) != xx0
    println("DONE")

    print("  method propagate!")
    x1 = copy(x0)
    propagate!(x1, model, 0.01)
    @test x1 != x0
    print("-")
    xx1 = copy(xx0)
    propagate!(xx1, model, 0.01)
    print("-")
    @test xx1 != xx0
    println("DONE")

end; #LinearDiffusionStateModel.jl

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
#=
P26 (**) Generate the combinations of K distinct objects chosen from the N elements of a list
In how many ways can a committee of 3 be chosen from a group of 12 people? We all know that there are C(12,3) = 220 possibilities (C(N,K) denotes the well-known binomial coefficients). For pure mathematicians, this result may be great. But we want to really generate all the possibilities in a list.

Example:
* (combination 3 '(a b c d e f))
((A B C) (A B D) (A B E) ... )
=#

function filter_colls(colls, len)
    filter(coll -> length(coll) == len, colls)
end

function append_indices(indices, last_index)
    next_index = indices[end] + 1
    map(i -> vcat(indices, i), next_index:last_index)
end

function combination(n, coll)
    function _combination(index_list, count)
        if count == n
            index_list
        else
            new_index_list = map(indices -> append_indices(indices, length(coll)), index_list)
            flattened = reduce(vcat, new_index_list)
            _combination(filter_colls(flattened, count + 1), count + 1)
        end
    end

    index_list = _combination(map(idx -> [idx], 1:length(coll)), 1)
    map(indices -> map(i -> coll[i], indices), index_list)
end

coll = [:a, :b, :c, :d, :e, :f]
combos = combination(3, coll)
@assert length(combos) == 20
@assert combos[1] == [:a, :b, :c]
@assert combos[2] == [:a, :b, :d]
@assert combos[3] == [:a, :b, :e]
@assert combos[4] == [:a, :b, :f]
@assert combos[5] == [:a, :c, :d]
@assert combos[6] == [:a, :c, :e]
@assert combos[7] == [:a, :c, :f]
@assert combos[8] == [:a, :d, :e]
@assert combos[9] == [:a, :d, :f]
@assert combos[10] == [:a, :e, :f]
@assert combos[11] == [:b, :c, :d]
@assert combos[12] == [:b, :c, :e]
@assert combos[13] == [:b, :c, :f]
@assert combos[14] == [:b, :d, :e]
@assert combos[15] == [:b, :d, :f]
@assert combos[16] == [:b, :e, :f]
@assert combos[17] == [:c, :d, :e]
@assert combos[18] == [:c, :d, :f]
@assert combos[19] == [:c, :e, :f]
@assert combos[20] == [:d, :e, :f]

println("Tests passed: JL-26.jl")


text:
"Abstract type that deals with global boundary-boundary matrices"
abstract BBSolver

type BBMatrix{T<:Number} <: BBSolver
    nb::Int
    nbslv::Int
    A::Array{T,2}
    fact::Base.LinAlg.LU{T,Array{T,2}}
    function BBMatrix(dof::DofMap)
        nb = nbmodes(dof)
        nbslv = nbslvmodes(dof)
        new(nb, nbslv, zeros(T, nbslv, nbslv))
    end
end

function assemble!{T<:Number}(Ag::BBMatrix{T}, Ae, m)

    nm = length(m)
    n1 = Ag.nbslv+1
    for i = 1:nm
        ig = m[i]
        if ig < n1
            for k = 1:nm
                kg = m[k]
                if kg < n1
                    Ag.A[kg,ig] += Ae[k,i]
                end
            end
        end
    end
end

function trf!(Ag::BBMatrix)
    Ag.fact = lufact!(Ag.A)
end
trs!(Ag::BBMatrix, x) = A_ldiv_B!(Ag.fact, x)

"Global boundary-boundary Symmetric Positive-Definite matrix"
abstract BBSym <: BBSolver

type BBSymMatrix{T<:Number} <: BBSym
    nb::Int
    nbslv::Int
    A::Array{T,2}
    fact::Base.LinAlg.Cholesky{T,Array{T,2}}
    function BBSymMatrix(dof::DofMap)
        nb = nbmodes(dof)
        nbslv = nbslvmodes(dof)
        new(nb, nbslv, zeros(T, nbslv, nbslv))
    end

end

function assemble!{T<:Number}(Ag::BBSymMatrix{T}, Ae, m)

    nm = length(m)
    n1 = Ag.nbslv+1
    for i = 1:nm
        ig = m[i]
        if ig < n1
            for k = 1:nm
                kg = m[k]
                if kg < n1
                    Ag.A[kg,ig] += Ae[k,i]
                end
            end
        end
    end
end

function trf!(Ag::BBSymMatrix)
    Ag.fact = cholfact!(Ag.A)
end
trs!(Ag::BBSymMatrix, x) = A_ldiv_B!(Ag.fact, x)

"Global boundary-boundary symmetric tridiagonal matrices"
type BBSymTri{T<:Number} <: BBSym
    "Number of boundary modes"
    nb::Int
    "Number of boundary modes that should be solved"
    nbslv::Int
    "Diagonal elements of the matrix"
    D::Array{T,1}
    "Sub-diagonal elements"
    E::Array{T,1}

    "Tridiagonal Matrix"
    tri::SymTridiagonal{T}
    "Factorization of symmetric tridiagonal matrix"
    fact:: LDLt{T,SymTridiagonal{T}}

    function BBSymTri(dof::DofMap)
        nb = nbmodes(dof)
        nbslv = nbslvmodes(dof)
        D = zeros(T,nbslv)
        E = zeros(T, max(nbslv-1,0))
        tri = SymTridiagonal(D, E)
        new(nb, nbslv, D, E, tri)
    end
end

"LU (or Choleksy) decomposition of boundary-boundary system"
function trf!(Ag::BBSymTri)
    Ag.fact = ldltfact!(Ag.tri)
end

"Solving linear system using LU (or Cholesky) decomposition"
trs!(Ag::BBSymTri, x) = A_ldiv_B!(Ag.fact, x)

"""
Assembles the global boundary-boundary matrix.
"""
function assemble!{T}(Ag::BBSymTri{T}, Ae, m)
    np1 = Ag.nbslv + 1
    D = Ag.D
    E = Ag.E

    ig = m[1]
    kg = m[2]

    if ig < np1
        D[ig] += Ae[1,1]
        if kg < np1
            D[kg] += Ae[2,2]
            E[min(ig,kg)] += Ae[1,2]
        end
    elseif kg < np1
        D[kg] += Ae[2,2]
    end

end

"Global boundary-boundary tridiagonal matrices"
type BBTri{T<:Number} <: BBSolver
    "Number of boundary modes"
    nb::Int
    "Number of boundary modes that should be solved"
    nbslv::Int
    "Diagonal elements of the matrix"
    D::Array{T,1}
    "Lower sub-diagonal elements"
    Dl::Array{T,1}
    "Upper sub-diagonal elements"
    Du::Array{T,1}
    "Tridiagonal object"
    tri::Tridiagonal{T}
    fact::Base.LinAlg.LU{T,Tridiagonal{T}}
    """
    Inner constructor. First the matrix must be assembled.
    """
    function BBTri(dof::DofMap)
        nb = nbmodes(dof)
        nbslv = nbslvmodes(dof)
        D = zeros(T,nbslv)
        Dl = zeros(T, max(nbslv-1, 0))
        Du = zeros(T, max(nbslv-1, 0))
        tri = Tridiagonal(Dl, D, Du)
        new(nb, nbslv, D, Dl, Du, tri)
    end
end

function assemble!{T}(Ag::BBTri{T}, Ae, m)
    np1 = Ag.nbslv + 1
    D = Ag.D
    Dl = Ag.Dl
    Du = Ag.Du

    ig = m[1]
    kg = m[2]

    if ig < np1
        D[ig] += Ae[1,1]
        if kg < np1
            D[kg] += Ae[2,2]
            if ig < kg
                Du[ig] += Ae[1,2]
                Dl[ig] += Ae[2,1]
            elseif kg < ig
                Du[kg] = Ae[2,1]
                Dl[kg] = Ae[1,2]
            end
        end
    elseif kg < np1
        D[kg] += Ae[2,2]
    end

end

function trf!(Ag::BBTri)
    Ag.fact = lufact!(Ag.tri)
    return
end
trs!(Ag::BBTri, x) = A_ldiv_B!(Ag.fact, x)

type BBTriP{T<:Number} <: BBSolver
    "Number of independent boundary modes"
    nb::Int
    "Number of boundary modes that should be solved"
    nbslv::Int
    "Lower diagonal"
    Dl::Vector{T}
    "Main Diagonal"
    D::Vector{T}
    "Upper diagonal"
    Du::Vector{T}
    "Auxiliary memory"
    x2::Vector{T}
    an1::T
    b1::T
    bn1::T
    cn::T
    cn1::T
    tri::Tridiagonal{T}
    fact::LinAlg.LU{T,Tridiagonal{T}}

    function BBTriP(dof::DofMap)
        nb = nbmodes(dof)
        nbslv = nbslvmodes(dof)

        nbslv = nb-1
        D = zeros(T,nbslv)
        Dl = zeros(T,nbslv-1)
        Du = zeros(T,nbslv-1)
        x2 = zeros(T,nbslv)
        tri = Tridiagonal(Dl, D, Du)
        new(nb, nbslv, Dl, D, Du, x2, zero(T), zero(T), zero(T), zero(T), zero(T), tri)

    end

end

function assemble!{T}(Ag::BBTriP{T}, Ae, m)

    m1 = m[1]
    m2 = m[2]

    nb = Ag.nb
    nbslv = Ag.nb-1

    ig = m1
    kg = m2

    D = Ag.D
    Dl = Ag.Dl
    Du = Ag.Du

    if m2 < nb && m2 != 1 #Not the last element
        np1 = Ag.nbslv + 1

        if ig < np1
            D[ig] += Ae[1,1]
            if kg < np1
                D[kg] += Ae[2,2]
                if ig < kg
                    Du[ig] += Ae[1,2]
                    Dl[ig] += Ae[2,1]
                elseif kg < ig
                    Du[kg] = Ae[2,1]
                    Dl[kg] = Ae[1,2]
                end
            end
        elseif kg < np1
            D[kg] += Ae[2,2]
        end
    elseif m2 == nb # It is the next to last
        Ag.an1   += Ae[2,2]
        Ag.cn    += Ae[1,2]
        Ag.bn1   += Ae[2,1]
        D[nbslv] += Ae[1,1]
    elseif m2 == 1
        Ag.b1    += Ae[2,1]
        Ag.cn1   += Ae[1,2]
        D[1]     += Ae[2,2]
        Ag.an1   += Ae[1,1]
    end

end

function trf!(Ag::BBTriP)
    Ag.fact = lufact!(Ag.tri)
    return
end
function trs!(Ag::BBTriP, x)
    n = Ag.nb
    n1 = n-1
    qn1 = view(x, 1:n1, :)
    nrhs = size(x,2)
    A_ldiv_B!(Ag.fact, qn1)
    x2 = Ag.x2
    for i = 1:nrhs
        fill!(x2, 0)
        x2[1] = -Ag.b1
        x2[n1] = -Ag.cn
        A_ldiv_B!(Ag.fact, x2)
        x[n,i] = (x[n,i] - Ag.cn1*qn1[1,i] - Ag.bn1*qn1[n1,i]) /
            (Ag.an1 + Ag.cn1*x2[1] + Ag.bn1*x2[n1])
        for k = 1:n1
            x[k,i] += x[n,i]*x2[k]
        end
    end
    x

end

type BBSymTriP{T<:Number} <: BBSym
    "Number of independent boundary modes"
    nb::Int
    "Number of boundary modes that should be solved"
    nbslv::Int
    "Main Diagonal"
    D::Vector{T}
    "Upper diagonal"
    Du::Vector{T}
    "Auxiliary memory"
    x2::Vector{T}
    an1::T
    cn::T
    cn1::T
    "Tridiagonal Matrix"
    tri::SymTridiagonal{T}
    "Factorization of symmetric tridiagonal matrix"
    fact:: LDLt{T,SymTridiagonal{T}}

    function BBSymTriP(dof::DofMap)
        nb = nbmodes(dof)
        nbslv = nbslvmodes(dof)

        nbslv = nb-1
        D = zeros(T,nb-1)
        Du = zeros(T,nb-2)
        x2 = zeros(T,nb-1)
        tri = SymTridiagonal(D, Du)
        new(nb, nbslv, D, Du, x2, zero(T), zero(T), zero(T), tri)
    end
end

function assemble!{T<:Number}(Ag::BBSymTriP{T}, Ae, m)

    m1 = m[1]
    m2 = m[2]

    nb = Ag.nb
    nbslv = Ag.nb-1

    ig = m1
    kg = m2

    D = Ag.D
    Du = Ag.Du

    if m2 < nb && m2 != 1 # Not the last element (if periodic)

        np1 = Ag.nbslv + 1

        if ig < np1
            D[ig] += Ae[1,1]
            if kg < np1
                D[kg] += Ae[2,2]
                if ig < kg
                    Du[ig] += Ae[1,2]
                elseif kg < ig
                    Du[kg] = Ae[2,1]
                end
            end
        elseif kg < np1
            D[kg] += Ae[2,2]
        end
    elseif m2 == nb # It is  the next to last
        Ag.an1   += Ae[2,2]
        Ag.cn    += Ae[1,2]
        D[nbslv] += Ae[1,1]
    elseif m2 == 1
        Ag.cn1   += Ae[1,2]
        D[1]     += Ae[2,2]
        Ag.an1   += Ae[1,1]
    end

end

function trf!(Ag::BBSymTriP)
    Ag.fact = ldltfact!(Ag.tri)
    return
end
function trs!(Ag::BBSymTriP, x)
    n = Ag.nb
    n1 = n-1
    qn1 = view(x, 1:n1, :)
    nrhs = size(x,2)
    A_ldiv_B!(Ag.fact, qn1)
    x2 = Ag.x2
    for i = 1:nrhs
        fill!(x2, 0)
        x2[1] = -Ag.cn1
        x2[n1] = -Ag.cn
        A_ldiv_B!(Ag.fact, x2)
        x[n,i] = (x[n,i] - Ag.cn1*qn1[1,i] - Ag.cn*qn1[n1,i]) /
            (Ag.an1 + Ag.cn1*x2[1] + Ag.cn*x2[n1])
        for k = 1:n1
            x[k,i] += x[n,i]*x2[k]
        end
    end
    x
end

@enum MatType FULLMAT TRIMAT TRIPMAT

"""
Generic matrix interface for 1d problems.

If the problem is very small, a full matrix will be used.
If it is larger and periodic, a BBTriP will be used.
Else a BBTri is used.
"""
type BBMatrix1d{T<:Number} <: BBSolver
    nb::Int
    nbslv::Int
    mtype::MatType
    periodic::Bool
    mat::BBMatrix{T}
    tri::BBTri{T}
    trip::BBTriP{T}
    #BBMatrix1d(nb, nbslv, mtype, periodic) = new(nb, nbslv, mtype, periodic)
    function BBMatrix1d(dof::DofMap)
        nb = nbmodes(dof)
        nbslv = nbslvmodes(dof)
        periodic = false
        if dof.bmap[1,1] == dof.bmap[2,end]
            periodic = true
        end

        if nbslv <= 3
            mtype = FULLMAT
        elseif periodic == false
            mtype = TRIMAT
        else
            mtype = TRIPMAT
        end

        bb = new(nb, nbslv, mtype, periodic)

        if mtype == FULLMAT
            bb.mat = BBMatrix{T}(dof)
        elseif mtype == TRIMAT
            bb.tri = BBTri{T}(dof)
        else
            bb.trip = BBTriP{T}(dof)
        end

        return bb
    end
end

"""
Generic symmetric matrix interface for 1d problems.

If the problem is very small, a full matrix will be used.
If it is larger and periodic, a BBSymTriP will be used.
Else a BBSymTri is used.
"""
type BBSymMatrix1d{T<:Number} <: BBSolver
    nb::Int
    nbslv::Int
    mtype::MatType
    periodic::Bool
    mat::BBSymMatrix{T}
    tri::BBSymTri{T}
    trip::BBSymTriP{T}
    function BBSymMatrix1d(dof::DofMap)
        nb = nbmodes(dof)
        nbslv = nbslvmodes(dof)
        periodic = false
        if dof.bmap[1,1] == dof.bmap[2,end]
            periodic = true
        end

        if nbslv <= 3
            mtype = FULLMAT
        elseif periodic == false
            mtype = TRIMAT
        else
            mtype = TRIPMAT
        end

        bb = new(nb, nbslv, mtype, periodic)

        if mtype == FULLMAT
            bb.mat = BBSymMatrix{T}(dof)
        elseif mtype == TRIMAT
            bb.tri = BBSymTri{T}(dof)
        else
            bb.trip = BBSymTriP{T}(dof)
        end

        return bb
    end
end

function assemble!{T<:Number}(Ag::Union{BBMatrix1d{T},BBSymMatrix1d{T}}, Ae, m)
    mt = Ag.mtype
    if mt == FULLMAT
        assemble!(Ag.mat, Ae, m)
    elseif mt == TRIMAT
        assemble!(Ag.tri, Ae, m)
    else
        assemble!(Ag.trip, Ae, m)
    end

end

function trf!(Ag::Union{BBMatrix1d,BBSymMatrix1d})
    mt = Ag.mtype
    if mt == FULLMAT
        trf!(Ag.mat)
    elseif mt == TRIMAT
        trf!(Ag.tri)
    else
        trf!(Ag.trip)
    end
end

function trs!(Ag::Union{BBMatrix1d,BBSymMatrix1d}, x)
    mt = Ag.mtype
    if mt == FULLMAT
        trs!(Ag.mat,x)
    elseif mt == TRIMAT
        trs!(Ag.tri,x)
    else
        trs!(Ag.trip,x)
    end

end


text:
struct QUInt{N}
  bits::Bits{N}
end

QUInt{N}(v::Integer = 0) where N = QUInt{N}(Bits{N}(v))

const uint_types = [UInt8,UInt16,UInt32,UInt64]
@generated uint_type(::Type{QUInt{N}}) where N = uint_types[ceil(Int,N/8)]

bits(x::QUInt) = x.bits
repr(::Type{QUInt{N}}, bits) where N = uint_type(QUInt{N})(int(bits))
Base.show(io::IO, x::QUInt) = print(io, "QUInt{$(nbits(x.bits))}()")

for N in [4,8,16,32,64]
  QUIntN = Symbol(:QUInt,N)
  @eval export $QUIntN
  @eval const $QUIntN = QUInt{$N}
end


text:
using Test, Documenter, Parameters, StructMapping

@with_kw struct A
    a::Float64
    b::String
end

dict_a = Dict("a"=>1.0, "b"=>"b")

@testset "basic" begin
    @test convertdict(A, dict_a) == A(1.0, "b")
end

@with_kw struct A2
    a::Float64
    b::Vector{String}
end

dict_a2 = Dict("a"=>1.0, "b"=>["b"])

@testset "basic2" begin
    a2 = convertdict(A2, dict_a2)
    @test a2.a == 1.0
    @test a2.b == ["b"]
end

@with_kw struct B
    a::A
    b::Int64 = 0
end

dict_b = Dict("a"=>dict_a, "b"=>4)

@testset "nested" begin
    @test convertdict(B, dict_b) == B(A(1.0, "b"), 4)
end

dict_b2 = Dict("a"=>dict_a)

@testset "default" begin
    @test convertdict(B, dict_b2) == B(A(1.0, "b"), 0)
end

@with_kw struct C
    a::Vector{A}
end

dict_a2 = Dict("a"=>2.0, "b"=>"b2")
dict_c = Dict("a"=>[dict_a, dict_a2])

@testset "vector" begin
    @test convertdict(C, dict_c).a == [A(1.0, "b"), A(2.0, "b2")]
end

@with_kw struct D
    a::Union{A, Nothing} = nothing
end

@testset "union" begin
    @test convertdict(D, Dict()) == D()
    @test convertdict(D, dict_b2) == D(A(1.0, "b"))
end

@with_kw struct E
    a::Union{Vector{A}, Nothing} = nothing
end

@testset "union_vector" begin
    @test convertdict(E, Dict()) == E()
    @test convertdict(E, dict_c).a == [A(1.0, "b"), A(2.0, "b2")]
end

@with_kw struct F
    b::B
    d::D
end

@testset "deeply_nested" begin
    f = convertdict(F, Dict("b"=>dict_b, "d"=>dict_b2))
    @test f.b == B(A(1.0, "b"), 4)
    @test f.d == D(A(1.0, "b"))
end

@with_kw struct G
    a::A
    g::Dict{String, Int64}
end

dict_g = Dict("a"=>dict_a, "g"=>Dict("g1"=>0, "g2"=>1))

@testset "dict" begin
    g = convertdict(G, dict_g)
    @test g.a == A(1.0, "b")
    @test g.g == Dict("g1"=>0, "g2"=>1)
end

@testset "doctest" begin
    DocMeta.setdocmeta!(StructMapping, :DocTestSetup, :(using StructMapping); recursive=true)
    doctest(StructMapping)
end


text:
@testset "$TEST $G" begin

g = readgraph(:lesmis, G)
@test typeof(g) == G
@test nv(g) == 77
@test ne(g) == 254
@test has_self_loops(g) == false

g = readgraph(:serengetifoodweb, G)
@test typeof(g) == DG
@test ne(g) == 592
@test nv(g) == 161
@test has_self_loops(g) == true

g = readgraph(:adjnoun, G)
@test nv(g) ==  112
@test ne(g) ==   425
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:as22july06, G)
@test nv(g) == 22963
@test ne(g) ==  48436
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:astroph, G)
@test nv(g) == 16706
@test ne(g) ==  121251
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:celegansneural, G)
@test nv(g) ==  297
# @test ne(g) == 2359
@test ne(g) == 2345 # multiedges?
@test typeof(g) == DG
@test has_self_loops(g) == false

g = readgraph(:condmat, G)
@test nv(g) == 16726
@test ne(g) ==  47594
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:condmat2003, G)
@test nv(g) ==31163
@test ne(g) == 120029
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:condmat2005, G)
@test nv(g) ==40421
@test ne(g) == 175693
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:dolphins, G)
@test nv(g) ==  62
@test ne(g) ==   159
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:emailenron, G)
@test nv(g) == 36692
@test ne(g) ==  183831
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:football, G)
@test nv(g) ==  115
# @test ne(g) == 615 # multiedges?
@test ne(g) ==  613
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:hepth, G)
@test nv(g) == 8361
@test ne(g) ==  15751
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:karate, G)
@test nv(g) ==  34
@test ne(g) ==   78
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:lesmis, G)
@test nv(g) ==  77
@test ne(g) ==   254
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:netscience, G)
@test nv(g) ==  1589
@test ne(g) ==   2742
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:pgpstrong2009, G)
@test nv(g) ==39796
@test ne(g) == 301498
@test typeof(g) == DG
@test has_self_loops(g) == false

g = readgraph(:polblogs, G)
@test nv(g) ==  1490
# @test ne(g) ==  19090
@test ne(g) ==  19025 #multiedges
@test typeof(g) == DG
@test has_self_loops(g) == true

g = readgraph(:polbooks, G)
@test nv(g) ==  105
@test ne(g) ==  441
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:power, G)
@test nv(g) ==  4941
@test ne(g) ==   6594
@test typeof(g) == G
@test has_self_loops(g) == false

g = readgraph(:serengetifoodweb, G)
@test nv(g) ==  161
@test ne(g) ==   592
@test typeof(g) == DG
@test has_self_loops(g) == true

end #testset


text:
function TriSquareMesh(
    m::Int,
    n::Int = m,
)::Tuple{Array{Float64,2},Array{Int,2}}

    p = [[x; y] for x = 0.0:1.0/Float64(m):1.0, y = 0.0:1.0/Float64(n):1.0][:]
    t = [
        [
            [i + (j - 1) * (m + 1) i + (j - 1) * (m + 1) + 1 i + j * (m + 1)],
            [i + (j - 1) * (m + 1) + 1 i + j * (m + 1) + 1 i + j * (m + 1)],
        ] for i = 1:m, j = 1:n
    ]

    return hcat(p...), hcat(hcat(t'...)...)
end

function QuadSquareMesh(
    m::Int,
    n::Int = m,
)::Tuple{Array{Float64,2},Array{Int,2}}

    p = [[x; y] for x = 0.0:1.0/Float64(m):1.0, y = 0.0:1.0/Float64(n):1.0][:]
    t = [ [i + (j - 1) * (m + 1)    i + (j - 1) * (m + 1) + 1  i + j * (m + 1) + 1   i + j * (m + 1)]
          for j = 1:m, i = 1:n ]

    return hcat(p...), hcat(t'...)
end

function SquareMesh(
    m::Int,
    n::Int,
    elemtype::Int,
)::Tuple{Array{Float64,2},Array{Int,2}}

    if (elemtype==0)
        p, t = TriSquareMesh(m,n)
    else
        p, t = QuadSquareMesh(m,n)
    end

    return p, t
end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
# Tests for io.jl

test_pdb_filepath_1 = testfile("1AKE.pdb")
test_dssp_filepath_1 = testfile("4AKE.dssp")
test_pdb_filepath_2 = testfile("1CTR_H.pdb")
test_dssp_filepath_2 = testfile("1CTR.dssp")
test_pocket_points = testfile("1AKE_pocket_points.pdb")
test_ligsite_filepath = testfile("1DVM_pocket_all.pdb")

@testset "IO" begin
    dssps = readdssp(test_dssp_filepath_1)
    @test dssps["21A"] == 'H'
    @test dssps["43A"] == ' '

    atoms = readpdb(test_pdb_filepath_2)
    dssps = readdssp(test_dssp_filepath_2, atoms)
    @test length(dssps) == 142
    @test dssps["75A"] == '-'

    atoms = readpdb(test_pdb_filepath_1)
    @test isa(atoms, Array{Atom,1})
    @test length(atoms) == 3312
    atom = atoms[10]
    @test atom.atom_name == "CA"
    @test atom.res_n == 2
    @test atom.coords[2] == 49.969
    atoms = readpdb(test_pdb_filepath_1, hetatm=true)
    @test length(atoms) == 3804

    @test fixstring("Test string", 4) == "Test"
    @test fixstring('T', 4) == "   T"
    @test fixstring(0.123456, 5) == "0.123"
    @test fixstring(123, 5) == "  123"

    @test spaceatomname("CA", "C") == " CA "
    @test spaceatomname("CAB", "C") == " CAB"
    @test spaceatomname("1HB1", "H") == "1HB1"
    @test spaceatomname("CABA", "C") == "CABA"

    atoms = readpdb(test_pdb_filepath_1)
    # Temp file which is removed at the end
    temp_filepath = tempname()
    writepdb(temp_filepath, atoms)
    test_atoms = readpdb(temp_filepath)
    @test length(test_atoms) == 3312
    atom = test_atoms[10]
    @test atom.atom_name == "CA"
    @test atom.res_n == 2
    @test atom.coords[2] == 49.969
    rm(temp_filepath)

    pock_points = readpocketpoints(test_pocket_points)
    @test length(pock_points) == 18
    @test size(pock_points[2]) == (3,54)
    @test pock_points[6] == [
        37.08   37.08   37.08 ;
        42.827  42.827  43.827;
        31.49   32.49   31.49 ;
    ]

    coords, vols = readligsite(test_ligsite_filepath)
    @test size(coords) == (3,43)
    @test length(vols) == 43
    @test coords[:,5] == [8.863, 28.126, -27.957]
    @test vols[5] == 21

    lines = readpdblines(test_pdb_filepath_1)
    @test length(lines) == 3317
    @test lines[10] == "ATOM     10  CA  ARG A   2      27.437  49.969  37.786  1.00 25.76           C  "
end


text:
## support routines

function clopts(opts...)
    args = ["", opts...]
    API.LLVMParseCommandLineOptions(length(args), args, C_NULL)
end

"""
    add_symbol(name, ptr)

Permanently add the symbol `name` with the value `ptr`. These symbols are searched
before any libraries.
"""
add_symbol(name, ptr) = LLVM.API.LLVMAddSymbol(name, ptr)

"""
    load_library_permantly(path)

This function permanently loads the dynamic library at the given path.
It is safe to call this function multiple times for the same library.
"""
load_library_permantly(path) = LLVM.API.LLVMLoadLibraryPermanently(path)

"""
    find_symbol(name)

Search the global symbols for `name` and return the pointer to it.
"""
find_symbol(name) = LLVM.API.LLVMSearchForAddressOfSymbol(name)


text:
module LawTests

using Test
using FillArrays
using Gridap.Arrays
using Gridap.ReferenceFEs
using Gridap.Fields
using Gridap.CellData
using Gridap.Fields: OtherMockBasis

ndofs = num_dofs(QUAD4)
b = get_shapefuns(QUAD4)
c = fill(1,ndofs)
f = OtherMockBasis{2}(ndofs)

ncells = 10
cl = fill(c,ncells)
fl = Fill(f,ncells)
ϕl = lincomb(fl,cl)

degree = 3
quad = CellQuadrature(degree,[QUAD,],Fill(1,ncells))
q = get_coordinates(quad)
dv = GenericCellField(Fill(b,ncells),ϕl,Val(true),Fill((Base.OneTo(ndofs),),ncells),Val((:,)))
u = convert_to_cell_field(x->x[2]^2,ϕl)

@law g(u) = u^2
@law h(u,dv) = (1+u)*dv

cf = g(u)
@test isa(cf,CellField)
@test get_metasize(cf) == ()
test_array(evaluate(cf,q),evaluate(operate(g,u),q),≈)

cf = h(u,dv)
@test isa(cf,CellField)
@test get_metasize(cf) == (:,)
test_array(evaluate(cf,q),evaluate(operate(h,u,dv),q),≈)

end # module


text:
"""
Lightweight type to handle CUDA texture objects inside CUDAnative.jl kernels.
Textures are fetched through indexing operations on `CuTexture`/`CuDeviceTexture` objects, e.g., `cutexture2d[0.2f0, 0.2f0]`.
"""
struct CuDeviceTexture{T,N}
    handle::CUtexObject
    dims::Dims{N}
end

Adapt.adapt_storage(::CUDAnative.Adaptor, t::CuTexture{T,N}) where {T,N} = CuDeviceTexture{T,N}(t.handle, size(t.mem))

@inline function tex1d(texObject::Int64, x::Float32)::Tuple{Int32,Int32,Int32,Int32}
    Base.llvmcall(("declare [4 x i32] @llvm.nvvm.tex.unified.1d.v4s32.f32(i64, float)",
        "%3 =  call [4 x i32] @llvm.nvvm.tex.unified.1d.v4s32.f32(i64 %0, float %1)\nret [4 x i32] %3"),
        Tuple{Int32,Int32,Int32,Int32},
        Tuple{Int64,Float32}, texObject, x)
end
@inline function tex2d(texObject::Int64, x::Float32, y::Float32)::Tuple{Int32,Int32,Int32,Int32}
    Base.llvmcall(("declare [4 x i32] @llvm.nvvm.tex.unified.2d.v4s32.f32(i64, float, float)",
        "%4 =  call [4 x i32] @llvm.nvvm.tex.unified.2d.v4s32.f32(i64 %0, float %1, float %2)\nret [4 x i32] %4"),
        Tuple{Int32,Int32,Int32,Int32},
        Tuple{Int64,Float32,Float32}, texObject, x, y)
end
@inline function tex3d(texObject::Int64, x::Float32, y::Float32, z::Float32)::Tuple{Int32,Int32,Int32,Int32}
    Base.llvmcall(("declare [4 x i32] @llvm.nvvm.tex.unified.3d.v4s32.f32(i64, float, float, float)",
        "%5 =  call [4 x i32] @llvm.nvvm.tex.unified.3d.v4s32.f32(i64 %0, float %1, float %2, float %3)\nret [4 x i32] %5"),
        Tuple{Int32,Int32,Int32,Int32},
        Tuple{Int64,Float32,Float32,Float32}, texObject, x, y, z)
end
@inline texXD(t::CuDeviceTexture{<:Any,1}, x::Real)::Tuple{Int32,Int32,Int32,Int32} = tex1d(reinterpret(Int64, t.handle), convert(Float32, x))
@inline texXD(t::CuDeviceTexture{<:Any,2}, x::Real, y::Real)::Tuple{Int32,Int32,Int32,Int32} = tex2d(reinterpret(Int64, t.handle), convert(Float32, x), convert(Float32, y))
@inline texXD(t::CuDeviceTexture{<:Any,3}, x::Real, y::Real, z::Real)::Tuple{Int32,Int32,Int32,Int32} = tex3d(reinterpret(Int64, t.handle), convert(Float32, x), convert(Float32, y), convert(Float32, z))

@inline reconstruct(::Type{T}, x::Int32) where {T <: Union{Int32,UInt32,Int16,UInt16,Int8,UInt8}} = unsafe_trunc(T, x)
@inline reconstruct(::Type{Float32}, x::Int32) = reinterpret(Float32, x)
@inline reconstruct(::Type{Float16}, x::Int32) = convert(Float16, reinterpret(Float32, x))

@inline reconstruct(::Type{T}, i32_x4::NTuple{4,Int32}) where T = reconstruct(T, i32_x4[1])
@inline reconstruct(::Type{NTuple{2,T}}, i32_x4::NTuple{4,Int32}) where T = (reconstruct(T, i32_x4[1]),
                                                                             reconstruct(T, i32_x4[2]))
@inline reconstruct(::Type{NTuple{4,T}}, i32_x4::NTuple{4,Int32}) where T = (reconstruct(T, i32_x4[1]),
                                                                             reconstruct(T, i32_x4[2]),
                                                                             reconstruct(T, i32_x4[3]),
                                                                             reconstruct(T, i32_x4[4]))

@inline function cast(::Type{T}, x) where {T}
    @assert sizeof(T) == sizeof(x)
    r = Ref(x)
    GC.@preserve r begin
       unsafe_load(convert(Ptr{T}, Base.unsafe_convert(Ptr{Cvoid}, r)))
    end
end

@inline function _fetch(t::CuDeviceTexture{T}, idcs::NTuple{<:Any,Real}) where {T}
    i32_x4 = texXD(t, idcs...)
    Ta = cuda_texture_alias_type(T)
    cast(T, reconstruct(Ta, i32_x4))
end

@inline _expand(x, n) = x * n
@inline (t::CuDeviceTexture{T,1})(x::R) where {T,R <: Real} = _fetch(t, _expand.((x,), t.dims))
@inline (t::CuDeviceTexture{T,2})(x::R, y::R) where {T,R <: Real} = _fetch(t, _expand.((x, y), t.dims))
@inline (t::CuDeviceTexture{T,3})(x::R, y::R, z::R) where {T,R <: Real} = _fetch(t, _expand.((x, y, z), t.dims))

@inline _offset(x) = x - 0.5f0
@inline Base.getindex(t::CuDeviceTexture{T,1}, x::R) where {T,R <: Real} = _fetch(t, _offset.((x,)))
@inline Base.getindex(t::CuDeviceTexture{T,2}, x::R, y::R) where {T,R <: Real} = _fetch(t, _offset.((x, y)))
@inline Base.getindex(t::CuDeviceTexture{T,3}, x::R, y::R, z::R) where {T,R <: Real} = _fetch(t, _offset.((x, y, z)))


text:
"""
    NegativeBinomial

* Parametrization
parametrized in r, p

* Score

* Fisher Information

* `time_varying_params` map.

* Default link
"""
NegativeBinomial

function score!(score_til::Matrix{T}, y::Int, ::Type{NegativeBinomial}, param::Matrix{T}, t::Int) where T
    score_til[t, 1] = digamma(y + param[t, 1]) - digamma(param[t, 1]) + log(param[t, 2])
    score_til[t, 2] = param[t, 1] / param[t, 2] - y / (one(T) - param[t, 2])
    return
end

function fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{NegativeBinomial}, param::Matrix{T}, t::Int) where T
    return error("Fisher information not implemented for NegativeBinomial distribution.")
end

function log_likelihood(::Type{NegativeBinomial}, y::Vector{Int}, param::Matrix{T}, n::Int) where T
    loglik = zero(T)
    for t in 1:n
        loglik += loggamma(y[t] + param[t, 1]) - logfactorial(y[t]) - loggamma(param[t, 1]) +
                param[t, 1] * log(param[t, 2]) + y[t] * log(1 - param[t, 2])
    end
    return -loglik
end

# Links
function link!(param_tilde::Matrix{T}, ::Type{NegativeBinomial}, param::Matrix{T}, t::Int) where T
    param_tilde[t, 1] = link(LogLink, param[t, 1], zero(T))
    param_tilde[t, 2] = link(LogitLink, param[t, 2], zero(T), one(T))
    return
end
function unlink!(param::Matrix{T}, ::Type{NegativeBinomial}, param_tilde::Matrix{T}, t::Int) where T
    param[t, 1] = unlink(LogLink, param_tilde[t, 1], zero(T))
    param[t, 2] = unlink(LogitLink, param_tilde[t, 2], zero(T), one(T))
    return
end
function jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{NegativeBinomial}, param::Matrix{T}, t::Int) where T
    aux.jac[1] = jacobian_link(LogLink, param[t, 1], zero(T))
    aux.jac[2] = jacobian_link(LogitLink, param[t, 2], zero(T), one(T))
    return
end

# utils
function update_dist(::Type{NegativeBinomial}, param::Matrix{T}, t::Int) where T
    return NegativeBinomial(param[t, 1], param[t, 2])
end

function params_sdm(d::NegativeBinomial)
    return Distributions.params(d)
end

function num_params(::Type{NegativeBinomial})
    return 2
end


text:
using Test
using TestSetExtensions
using FluxTraining
using Animations
using Colors
using FluxTraining: EpochEnd, LearningRate, getdataiter, protect, Protected,
    ProtectedException, Read, Write, Events, Phases
using Flux: trainable

include("./testdata.jl");


text:
function maximize(
    domain::AbstractDomain,
    u::AbstractEigenfunction,
    λ::arb,
    n::Integer;
    store_trace = false,
    show_progress = false,
)
    N = length(coefficients(u))
    f = t -> u(boundary_parameterization(t, domain, n), λ, boundary = n)

    if show_progress
        enclose_progress(state) = begin
            sticky =
                radius(state.maxenclosure) / abs(state.maxenclosure) <= 1e-2 ? :done : true
            @info "Computing enclosure of maximum" state sticky = sticky
            false
        end
    else
        enclose_progress = nothing
    end

    enclosemaximum(
        f,
        zero(λ),
        one(λ),
        absmax = true,
        evaltype = :taylor,
        n = N,
        atol = 0,
        rtol = 1e-2,
        store_trace = store_trace,
        extended_trace = store_trace,
        callback = enclose_progress,
    )
end

function maximize(
    domain::AbstractDomain,
    u::AbstractEigenfunction,
    λ::arb;
    store_trace = false,
    show_progress = false,
)
    boundaries = active_boundaries(domain, u)
    m = zero(λ)

    traces = Dict()
    for boundary in boundaries
        m2 = maximize(
            domain,
            u,
            λ,
            boundary,
            store_trace = store_trace,
            show_progress = show_progress,
        )
        if store_trace
            m2, trace = m2
            traces[boundary] = trace
        end
        m = max(m, m2)
        if !isfinite(m)
            return m
        end
    end

    if store_trace
        return m, traces
    else
        return m
    end
end

function enclose_eigenvalue(
    domain::AbstractDomain,
    u::AbstractEigenfunction,
    λ::arb;
    rigorous_norm = true,
    show_progress = false,
    store_trace = false,
    extended_trace = false,
)
    @timeit_debug "maximize" m =
        maximize(domain, u, λ, store_trace = extended_trace, show_progress = show_progress)
    if extended_trace
        m, trace = m
    end

    @timeit_debug "norm" begin
        if rigorous_norm
            n = norm(domain, u, λ)
        else
            n = norm(domain, u, λ, numpoints = 4length(coefficients(u)), warn = false)
        end
    end

    ϵ = sqrt(area(domain)) * m / n
    lower = λ / (1 + getinterval(ϵ)[2])
    upper = λ / (1 - getinterval(ϵ)[2])
    if lower <= upper
        enclosure = setinterval(lower, upper)
    else
        enclosure = ball(λ, domain.parent(Inf))
    end

    if extended_trace
        return enclosure, n, m, trace
    elseif store_trace
        return enclosure, n, m
    else
        return enclosure
    end
end

function enclose_eigenvalue_approx(
    domain::AbstractDomain,
    u::AbstractEigenfunction,
    λ::arb;
    max_numpoints::Integer = 16length(coefficients(u)),
    norm_numpoints::Integer = 8length(coefficients(u)),
    store_trace::Bool = false,
    extended_trace::Bool = false,
)
    ## Approximate maximum on boundary
    # Points to evaluate u on
    pts, bds = boundary_points(domain, u, length(coefficients(u)), max_numpoints)
    values = similar(pts, arb)
    Threads.@threads for i in eachindex(pts)
        values[i] = u(pts[i], λ, boundary = bds[i])
    end
    m = zero(λ)
    for v in values
        m = max(m, abs(v))
    end

    ## Approximate norm
    n = norm(domain, u, λ, numpoints = norm_numpoints, warn = false)

    ## Compute enclosure
    ϵ = sqrt(area(domain)) * m / n
    lower = λ / (1 + getinterval(ϵ)[2])
    upper = λ / (1 - getinterval(ϵ)[2])
    if lower <= upper
        enclosure = setinterval(lower, upper)
    else
        enclosure = ball(λ, domain.parent(Inf))
    end

    if extended_trace
        return enclosure, n, m, Dict("values" => values, "boundaries" => bds)
    elseif store_trace
        return enclosure, n, m
    else
        return enclosure
    end
end


text:
@testset "Translation" begin
    # crummy string translation to test against
    standard_genetic_code_dict = Dict{String,Char}(
        "AAA" => 'K', "AAC" => 'N', "AAG" => 'K', "AAU" => 'N',
        "ACA" => 'T', "ACC" => 'T', "ACG" => 'T', "ACU" => 'T',
        "AGA" => 'R', "AGC" => 'S', "AGG" => 'R', "AGU" => 'S',
        "AUA" => 'I', "AUC" => 'I', "AUG" => 'M', "AUU" => 'I',
        "CAA" => 'Q', "CAC" => 'H', "CAG" => 'Q', "CAU" => 'H',
        "CCA" => 'P', "CCC" => 'P', "CCG" => 'P', "CCU" => 'P',
        "CGA" => 'R', "CGC" => 'R', "CGG" => 'R', "CGU" => 'R',
        "CUA" => 'L', "CUC" => 'L', "CUG" => 'L', "CUU" => 'L',
        "GAA" => 'E', "GAC" => 'D', "GAG" => 'E', "GAU" => 'D',
        "GCA" => 'A', "GCC" => 'A', "GCG" => 'A', "GCU" => 'A',
        "GGA" => 'G', "GGC" => 'G', "GGG" => 'G', "GGU" => 'G',
        "GUA" => 'V', "GUC" => 'V', "GUG" => 'V', "GUU" => 'V',
        "UAA" => '*', "UAC" => 'Y', "UAG" => '*', "UAU" => 'Y',
        "UCA" => 'S', "UCC" => 'S', "UCG" => 'S', "UCU" => 'S',
        "UGA" => '*', "UGC" => 'C', "UGG" => 'W', "UGU" => 'C',
        "UUA" => 'L', "UUC" => 'F', "UUG" => 'L', "UUU" => 'F',

        # translatable ambiguities in the standard code
        "CUN" => 'L', "CCN" => 'P', "CGN" => 'R', "ACN" => 'T',
        "GUN" => 'V', "GCN" => 'A', "GGN" => 'G', "UCN" => 'S'
    )

    function string_translate(seq::AbstractString)
        @assert length(seq) % 3 == 0
        aaseq = Vector{Char}(undef, div(length(seq), 3))
        for i in 1:3:length(seq) - 3 + 1
            aaseq[div(i, 3) + 1] = standard_genetic_code_dict[seq[i:i+2]]
        end
        return String(aaseq)
    end

    function check_translate(seq::AbstractString)
        return string_translate(seq) == String(translate(RNASequence(seq)))
    end

    global reps = 10
    for len in [1, 10, 32, 1000, 10000, 100000]
        @test all(Bool[check_translate(random_translatable_rna(len)) for _ in 1:reps])
    end

    # ambiguous codons
    @test translate(rna"YUGMGG") == aa"LR"
    @test translate(rna"GAYGARGAM") == aa"DEX"

    # BioSequences{RNAAlphabet{2}}
    @test translate(BioSequence{RNAAlphabet{2}}("AAAUUUGGGCCC")) == translate(rna"AAAUUUGGGCCC")
    @test translate(BioSequence{DNAAlphabet{2}}("AAATTTGGGCCC")) == translate(dna"AAATTTGGGCCC")

    # DNASequence
    @test translate(dna"ATGTAA") == aa"M*"

    # Alternative start codons
    @test translate(rna"GUGUAA", alternative_start = true) == aa"M*"

    @test_throws Exception translate(rna"ACGUACGU")  # can't translate non-multiples of three
    # can't translate N
    @test_throws Exception translate(rna"ACGUACGNU", allow_ambiguous_codons=false)

    # issue #133
    @test translate(rna"GAN") == aa"X"
end


text:
function DiscreteND1D(μ::Float64, δ::Float64, x::Vector{Float64}, xstep::Float64) # x should be sorted array with constant step
    if μ < x[1] || μ > x[end]
        error("The mean value is out of input array")
    end
    #xstep = x[2] - x[1]
    CDF = Vector{Float64}()
    D = Normal(μ, δ)
    #for i=1:length(x)
    for xi in x
        #cdfi = cdf(D, x[i]+xstep/2) - cdf(D, x[i]-xstep/2)
        cdfi = cdf(D, xi+xstep/2) - cdf(D, xi-xstep/2)
        if cdfi < 5.0e-3
            cdfi = 0.0
        end
        push!(CDF, cdfi)
    end
    return normalize!(CDF, 1)
    #return map(x-> round(x/Sum; digits=4), CDF)
end

function ContinuousND2D(Means::NTuple{2, Float64}, StdDeviations::NTuple{2, Float64})
    CovarianceMat = [StdDeviations[1]^2 0.0; 0.0 StdDeviations[2]^2]
    D = MvNormal([Means...], CovarianceMat)
    return D
end

function Sample2D(rng::AbstractRNG, Means::NTuple{2, Float64}, StdDeviations::NTuple{2, Float64})
    D = ContinuousND2D(Means, StdDeviations)
    pair = [0.0, 0.0]
    rand!(rng, D, pair)
    return pair
end


text:
# This code is not used yet
# It was moved from the work-in-progress code in proxy.jl

struct GenInfo
    code
    typeCode
    deps
    classList
    methodDicts
    fielddicts
end

struct GenArgInfo
    name::Symbol
    javaType::Type
    juliaType
    spec
end

const genned = Set()

hasClass(name::AbstractString) = hasClass(Symbol(name))
hasClass(name::Symbol) = name in genned
hasClass(gen, name::AbstractString) = hasClass(gen, Symbol(name))
hasClass(gen, name::Symbol) = name in genned || haskey(gen.methodDicts, string(name))

function genTypeDecl(name::AbstractString, supername::Symbol, gen)
    if string(name) != "String" && !haskey(types, Symbol(name)) && !haskey(gen.methodDicts, name)
        typeName = typeNameFor(name)
        push!(gen.typeCode, :(abstract type $typeName <: $supername end))
    end
end

function registerclass(name::AbstractString, classType::Type)
    registerclass(Symbol(name), classType)
end
function registerclass(name::Symbol, classType::Type)
    if !(classType <: Union{Array, String}) && !haskey(types, name)
        types[name] = classType
    end
    infoFor(classforname(string(name)))
end

gen(name::Symbol; genmode=:none, print=false, eval=true) = _gen(classforname(string(name)), genmode, print, eval)
gen(name::AbstractString; genmode=:none, print=false, eval=true) = _gen(classforname(name), genmode, print, eval)
gen(pxy::JProxy{T, C}) where {T, C} = gen(C)
gen(class::JClass; genmode=:none, print=false, eval=true) = _gen(class, genmode, eval)
function _gen(class::JClass, genmode, print, evalResult)
    n = legalClassName(class)
    gen = GenInfo()
    genClass(class, gen)
    if genmode == :deep
        while !isempty(gen.deps)
            cls = pop!(gen.deps)
            !hasClass(gen, cls) && genClass(classforname(string(cls)), gen)
        end
    else
        while !isempty(gen.deps)
            cls = pop!(gen.deps)
            !hasClass(gen, cls) && genType(classforname(string(cls)), gen)
        end
    end
    expr = :(begin $(gen.typeCode...); $(gen.code...); $(genClasses(getname.(gen.classList))...); end)
    if print
        for e in expr.args
            println(e)
        end
    end
    evalResult && eval(expr)
end

function genType(class, gen::GenInfo)
    name = getname(class)
    sc = superclass(class)
    push!(genned, Symbol(legalClassName(class)))
    if !isNull(sc)
        if !(Symbol(legalClassName(sc)) in genned)
            genType(getcomponentclass(sc), gen)
        end
        supertype = typeNameFor(sc)
        cType = componentType(supertype)
        genTypeDecl(name, cType, gen)
    else
        genTypeDecl(name, :java_lang, gen)
    end
end

genClass(class::JClass, gen::GenInfo) = genClass(class, gen, infoFor(class))
function genClass(class::JClass, gen::GenInfo, info::JClassInfo)
    name = getname(class)
    if !(Symbol(name) in genned)
        gen.fielddicts[legalClassName(class)] = fielddict(class)
        push!(gen.classList, class)
        sc = superclass(class)
        #@verbose("SUPERCLASS OF $name is $(isNull(sc) ? "" : "not ")null")
        push!(genned, Symbol(legalClassName(class)))
        if !isNull(sc)
            supertype = typeNameFor(sc)
            cType = componentType(supertype)
            !hasClass(gen, cType) && genClass(sc, gen)
            genTypeDecl(name, cType, gen)
        else
            genTypeDecl(name, :java_lang, gen)
        end
        genMethods(class, gen, info)
    end
end

GenInfo() = GenInfo([], [], Set(), [], Dict(), Dict())

function GenArgInfo(index, info::JMethodInfo, gen::GenInfo)
    javaType = info.argTypes[index]
    GenArgInfo(Symbol("a" * string(index)), javaType, argType(javaType, gen), argSpec(javaType, gen))
end

argType(t, gen) = t
argType(::Type{JavaObject{Symbol("java.lang.String")}}, gen) = String
argType(::Type{JavaObject{Symbol("java.lang.Object")}}, gen) = :JLegalArg
argType(::Type{<: Number}, gen) = Number
argType(typ::Type{JavaObject{T}}, gen) where T = :(JProxy{<:$(typeNameFor(T, gen)), T})

argSpec(t, gen) = t
argSpec(::Type{JavaObject{Symbol("java.lang.String")}}, gen) = String
argSpec(::Type{JavaObject{Symbol("java.lang.Object")}}, gen) = :JObject
argSpec(::Type{<: Number}, gen) = Number
argSpec(typ::Type{JavaObject{T}}, gen) where T = :(JProxy{<:$(typeNameFor(T, gen)), T})
argSpec(arg::GenArgInfo) = arg.spec

typeNameFor(T::Symbol, gen::GenInfo) = typeNameFor(string(T), gen)
function typeNameFor(T::AbstractString, gen::GenInfo)
    aType, dims = arrayinfo(T)
    c = dims != 0 ? aType : T
    csym = Symbol(c)
    if (dims == 0 || length(c) > 1) && !(csym in gen.deps) && !hasClass(gen, csym) && !get(typeInfo, c, genericFieldInfo).primitive
        push!(gen.deps, csym)
    end
    typeNameFor(T)
end

function argCode(arg::GenArgInfo)
    argname = arg.name
    if arg.juliaType == String
        argname
    elseif arg.juliaType == JLegalArg
        :(box($argname))
    elseif arg.juliaType == Number
        :($(arg.javaType)($argname))
    else
        argname
    end
end

function fieldEntry((name, fld))
    fieldType = JavaObject{Symbol(legalClassName(fld.owner))}
    name => :(jfield($(fld.typeInfo.class), $(string(name)), $fieldType))
end

function genMethods(class, gen, info)
    methodList = listmethods(class)
    classname = legalClassName(class)
    gen.methodDicts[classname] = methods = Dict()
    typeName = typeNameFor(classname, gen)
    classVar = Symbol("class_" * string(typeName))
    fieldsVar = Symbol("fields_" * string(typeName))
    methodsVar = Symbol("staticMethods_" * string(typeName))
    push!(gen.code, :($classVar = classforname($classname)))
    push!(gen.code, :($fieldsVar = Dict([$([fieldEntry(f) for f in gen.fielddicts[classname]]...)])))
    push!(gen.code, :($methodsVar = Set($([string(n) for (n, m) in info.methods if any(x->x.static, m)]))))
    push!(gen.code, :(function Base.getproperty(p::JProxy{T, C}, name::Symbol) where {T <: $typeName, C}
                      if (f = get($fieldsVar, name, nothing)) != nothing
                              getField(p, name, f)
                          else
                              JMethodProxy(name, $typeName, p, emptyset)
                          end
                      end))
    for nameSym in sort(collect(keys(info.methods)))
        name = string(nameSym)
        multiple = length(info.methods[nameSym]) > 1
        symId = 0
        for minfo in info.methods[nameSym]
            owner = javaType(minfo.owner)
            if isSame(class.ptr, minfo.owner.ptr)
                symId += 1
                args = (GenArgInfo(i, minfo, gen) for i in 1:length(minfo.argTypes))
                argDecs = (:($(arg.name)::$(arg.juliaType)) for arg in args)
                methodIdName = Symbol("method_" * string(typeName) * "__" * name * (multiple ? string(symId) : ""))
                callinfo = jnicalls[minfo.typeInfo.classname]
                push!(gen.code, :($methodIdName = getmethodid($(minfo.static), $classVar, $name, $(legalClassName(minfo.returnClass)), $(legalClassName.(minfo.argClasses)))))
                push!(gen.code, :(function (pxy::JMethodProxy{Symbol($name), <: $typeName})($(argDecs...))::$(genReturnType(minfo, gen))
                                      @verbose($("Generated method $name$(multiple ? "(" * string(symId) * ")" : "")"))
                                      $(genConvertResult(minfo.typeInfo.convertType, minfo, :(call(pxy.obj, $methodIdName, $(static ? callinfo.static : callinfo.instance), $(minfo.typeInfo.juliaType), ($(argSpec.(args)...),), $((argCode(arg) for arg in args)...)))))
                                  end))
            end
        end
    end
    push!(gen.code, :(push!(genned, Symbol($(legalClassName(class))))))
end

function genReturnType(methodInfo, gen)
    t = methodInfo.typeInfo.convertType
    if methodInfo.typeInfo.primitive || t <: String || t == Nothing
        t
    else
        :(JProxy{<:$(typeNameFor(methodInfo.returnType, gen))})
    end
end

genConvertResult(toType::Type{Bool}, info, expr) = :($expr != 0)
genConvertResult(toType::Type{String}, info, expr) = :(unsafe_string($expr))
genConvertResult(toType::Type{<:JBoxTypes}, info, expr) = :(unbox($(toType.parameters[1]), $expr))
function genConvertResult(toType, info, expr)
    if isVoid(info) || info.typeInfo.primitive
        expr
    else
        :(asJulia($toType, $expr))
    end
end

genClasses(classNames) = (:(registerclass($name, $(Symbol(typeNameFor(name))))) for name in reverse(classNames))


text:
using WebAPI
using HTTP
using JSON3
using Sockets
using Test

@testset "WebAPI.jl" begin
    include("params_query.jl")
    include("req_body.jl")
    include("res_constructors.jl")
    include("add_route.jl")
    include("utility_functions.jl")
    include("prints.jl")
end


text:
#=
# Active
normal at size except for active beam: normal at true dist

# Passive
normal at detect radius, except if in detect radius: normal in active beam

# beams
angles follow RHR from x,y axis (they start on x and go CCW)

Beam | covers (deg)
-------------------
1    | (0,45]
2    | (45,90]
etc.
=#

struct SubObsDist
    abeam::Int
    an::Normal{Float64}
    n::Normal{Float64}
end

function rand(rng::AbstractRNG, d::SubObsDist)
    o = MVector{8, Float64}()
    for i in 1:length(o)
        if i == d.abeam
            o[i] = rand(rng, d.an)
        else
            o[i] = rand(rng, d.n)
        end
    end
    return SVector(o)
end

function pdf(d::SubObsDist, o::Vec8)
    p = 1.0
    for i in 1:length(o)
        if i == d.abeam
            p *= pdf(d.an, o[i])
        else
            p *= pdf(d.n, o[i])
        end
    end
    return p
end

function active_beam(rel_pos::Pos)
    angle = atan2(rel_pos[2], rel_pos[1])
    while angle <= 0.0
        angle += 2*pi
    end
    bm = ceil(Int, 8*angle/(2*pi))
    return clamp(bm, 1, 8)
end

function observation(p::SubHuntPOMDP, a::Int, sp::SubState)
    rel_pos = sp.target - sp.own
    dist = norm(rel_pos)
    if a == PING
        abeam = active_beam(rel_pos)
        an = Normal(dist, p.active_std)
        n = Normal(p.size, p.active_std)
    else
        if dist <= p.passive_detect_radius
            abeam = active_beam(rel_pos)
            an = Normal(dist, p.passive_detected_std)
        else
            abeam = -1
            an = Normal(p.passive_detect_radius, p.passive_std)
        end
        n = Normal(p.passive_detect_radius, p.passive_std)
    end
    SubObsDist(abeam, an, n)
end


text:
#######################################################################
#
# Example of how to use Python and the XlsxWriter module to write
# simple array formulas.
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()
	wb = Workbook("array_formula.xlsx")

	ws = add_worksheet!(wb)

	# Write some test data.
	write!(ws, "B1", 500)
	write!(ws, "B2", 10)
	write!(ws, "B5", 1)
	write!(ws, "B6", 2)
	write!(ws, "B7", 3)
	write!(ws, "C1", 300)
	write!(ws, "C2", 15)
	write!(ws, "C5", 20234)
	write!(ws, "C6", 21003)
	write!(ws, "C7", 10000)

	# Write an array formula that returns a single value
	write_formula!(ws, "A1", "{=SUM(B1:C1*B2:C2)}")
	# 9500 in A1

	# Same as above but more verbose.
	write_array_formula!(ws, "A2", "{=SUM(B1:C1*B2:C2)}")
	# 9500 in A2

	# Write an array formula that returns a range of values
	write_array_formula!(ws, "A5:A7", "{=TREND(C5:C7,B5:B7)}")

	close(wb)

	isfile("array_formula.xlsx")
end

test()


text:
# This file is part of Kpax3. License is MIT.

function split!(ij::Vector{Int},
                neighbours::Vector{Int},
                S::Int,
                data::Matrix{UInt8},
                priorR::PriorRowPartition,
                priorC::PriorColPartition,
                settings::KSettings,
                support::MCMCSupport,
                state::AminoAcidState)
  # number of clusters after the split
  k = state.k + 1

  initsupportsplit!(ij, k, data, priorC, settings, support)

  hi = state.R[ij[1]]

  # logarithm of the product of sequential probabilities
  lq = 0.0

  # temporary / support variables
  u = 0
  lcp = zeros(Float64, 2)
  lw = zeros(Float64, 2)
  z = 0.0
  p = 0.0

  # allocate the neighbours of i and j
  for l in 1:S
    u = neighbours[l]
    lcp[1] = lcp[2] = 0.0

    # compute p(x_{u} | x_{hi,1:(u-1)}) and p(x_{u} | x_{hj,1:(u-1)})
    for b in 1:support.m
      lcp[1] += computeclusteriseqprobs!(data[b, u], b, priorC, support)
      lcp[2] += computeclusterjseqprobs!(data[b, u], b, priorC, support)
    end

    lw[1] = clusterweight(support.vi, priorR)
    lw[2] = clusterweight(support.vj, priorR)

    # (w1 * p1) / (w1 * p1 + w2 * p2) = 1 / (1 + (w2 * p2) / (w1 * p1))
    # => e^(-log(1 + e^(log(w2) + log(p2) - log(w1) - log(p1))))
    z = -log1p(exp(lw[2] + lcp[2] - lw[1] - lcp[1]))
    p = exp(z)

    if rand() <= p
      updateclusteri!(u, data, support)
      lq += z
    else
      updateclusterj!(u, data, support)
      lq += log1p(-p)
    end
  end

  updatelogmargliki!(priorC, support)
  updatelogmarglikj!(priorC, support)

  logratiopriorrowsplit!(k, priorR, support)
  logmargliksplit!(state.cl, state.k, hi, priorC, support)

  ratio = exp(support.lograR + support.logmlikcandidate - support.logmlik - lq)

  if ratio >= 1 || ((ratio > 0) && (rand() <= ratio))
    performsplit!(hi, k, priorC, settings, support, state)
  end

  nothing
end

function initsupportsplit!(ij::Vector{Int},
                           k::Int,
                           data::Matrix{UInt8},
                           priorC::AminoAcidPriorCol,
                           settings::KSettings,
                           support::MCMCSupport)
  resizesupport!(support, k)

  support.vi = 1
  support.vj = 1

  support.ui[1] = ij[1]
  support.uj[1] = ij[2]

  for b in 1:support.m
    initclusteriweights!(data[b, ij[1]], b, k, priorC, support)
    initclusterjweights!(data[b, ij[2]], b, k, priorC, support)
  end

  nothing
end

function performsplit!(hi::Int,
                       k::Int,
                       priorC::PriorColPartition,
                       settings::KSettings,
                       support::MCMCSupport,
                       state::AminoAcidState)
  resizestate!(state, k, settings)

  hj = findfirst(state.emptycluster)

  for a in 1:support.vj
    state.R[support.uj[a]] = hj
  end

  state.emptycluster[hj] = false

  h = 0
  for a in 1:length(state.emptycluster)
    if !state.emptycluster[a]
      h += 1
      state.cl[h] = a
    end
  end

  state.k = h

  state.v[hi] = support.vi
  state.v[hj] = support.vj

  for b in 1:support.m
    state.n1s[hi, b] = support.ni[b]
    state.n1s[hj, b] = support.nj[b]

    support.lp[1, hi, b] = support.lpi[1, b]
    support.lp[2, hi, b] = support.lpi[2, b]
    support.lp[3, hi, b] = support.lpi[3, b]
    support.lp[4, hi, b] = support.lpi[4, b]

    support.lp[1, hj, b] = support.lpj[1, b]
    support.lp[2, hj, b] = support.lpj[2, b]
    support.lp[3, hj, b] = support.lpj[3, b]
    support.lp[4, hj, b] = support.lpj[4, b]
  end

  copyto!(state.unit[hi], 1, support.ui, 1, support.vi)

  resize!(state.unit[hj], support.vj)
  copyto!(state.unit[hj], 1, support.uj, 1, support.vj)

  state.logpR += support.lograR

  support.logmlik = support.logmlikcandidate

  nothing
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(cell_centroid),PolytopalMesh{2,Float64,jFEMTools.Cell{2,4,4,1}},Int64})
    precompile(Tuple{typeof(jFEMTools.gradient_value),jFEMTools.Monomials{2,1,Float64},Int64,Tensors.Tensor{1,2,Float64,2}})
    precompile(Tuple{Type{jFEMTools.Monomials},Int64,Int64,Tensors.Tensor{1,2,Float64,2},Float64})
    precompile(Tuple{typeof(jFEMTools.value),jFEMTools.Monomials{2,1,Float64},Int64,Tensors.Tensor{1,2,Float64,2}})
    precompile(Tuple{typeof(jFEMTools.getverticesindices),PolytopalMesh{2,Float64,jFEMTools.Cell{2,4,4,1}},jFEMTools.EdgeIndex})
    precompile(Tuple{Type{DofHandler},PolytopalMesh{2,Float64,jFEMTools.Cell{2,3,3,1}},TrialFunction{PoissonVirtualElement{2}}})
    precompile(Tuple{typeof(rectangle_mesh),Type{jFEMTools.Cell{2,3,3,1}},Tuple{Int64,Int64},Tensors.Tensor{1,2,Float64,2},Tensors.Tensor{1,2,Float64,2}})
    precompile(Tuple{typeof(cell_diameter),PolytopalMesh{2,Float64,jFEMTools.Cell{2,3,3,1}},Int64})
    precompile(Tuple{typeof(jFEMTools.gradient_value),jFEMTools.Monomials{2,2,Float64},Int64,Tensors.Tensor{1,2,Float64,2}})
    precompile(Tuple{typeof(rectangle_mesh),Type{jFEMTools.Cell{2,4,4,1}},Tuple{Int64,Int64},Tensors.Tensor{1,2,Float64,2},Tensors.Tensor{1,2,Float64,2}})
    precompile(Tuple{typeof(apply!),SparseMatrixCSC{Bool,Int64},Array{Float64,1},Dirichlet{Float64}})
    precompile(Tuple{typeof(jFEMTools.value),jFEMTools.Monomials{2,2,Float64},Int64,Tensors.Tensor{1,2,Float64,2}})
    precompile(Tuple{typeof(get_vertices_matrix),PolytopalMesh{2,Float64,jFEMTools.Cell{2,3,3,1}}})
    precompile(Tuple{typeof(jFEMTools.create_sparsity_pattern),DofHandler{2,Float64}})
    precompile(Tuple{typeof(cell_volume),PolytopalMesh{2,Float64,jFEMTools.Cell{2,3,3,1}},Int64})
    precompile(Tuple{typeof(jFEMTools.get_2Dtesselation),Array{Tensors.Tensor{1,2,Float64,2},1}})
    precompile(Tuple{typeof(assemble!),jFEMTools.Assembler{Float64},Array{Int64,1},Array{Float64,2}})
    precompile(Tuple{typeof(cell_volume),PolytopalMesh{2,Float64,jFEMTools.Cell{2,4,4,1}},Int64})
    precompile(Tuple{typeof(jFEMTools.value),jFEMTools.Monomials{2,0,Float64},Int64,Tensors.Tensor{1,2,Float64,2}})
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
using Base.Test
using Meganet

for TYPE=[Float64,Float32]
    K      = getDenseKernel(TYPE,[32,18])
    nex    = 8
    Bin    = randn(TYPE,nFeatOut(K),4)
    Bout   = randn(TYPE,nFeatIn(K),3)
    nLayer = getTVNormLayer(TYPE,[8,4])
    L      = getDoubleSymLayer(TYPE,K,nLayer,Bin=Bin,Bout=Bout)
    @testset "doubleSymLayer (dense/TV) $TYPE" begin
     testAbstractMeganetElement(L)
    end

    K      = getDenseKernel(TYPE,[32,18])
    nex    = 8
    Bin    = randn(TYPE,nFeatOut(K),4)
    Bout   = randn(TYPE,nFeatIn(K),3)
    nLayer = getBatchNormLayer(TYPE,[8,4])
    L      = getDoubleSymLayer(TYPE,K,nLayer,Bin=Bin,Bout=Bout)
    @testset "doubleSymLayer (dense/Batch) $TYPE" begin
     testAbstractMeganetElement(L)
    end

    nImg   = [32 32]
    nc     = 16
    nex    = 50
    K      = getSparseConvKernel2D(TYPE,nImg,[3,3,1,nc])
    Bin    = randn(TYPE,nFeatOut(K),4)
    Bout   = randn(TYPE,nFeatIn(K),3)
    nLayer = getBatchNormLayer(TYPE,[prod(nImg),nc],isTrainable=false)
    L      = getDoubleSymLayer(TYPE,K,nLayer,Bin=Bin,Bout=Bout)
    @testset "doubleSymLayer (conv/Batch/not trainable) $TYPE" begin
     testAbstractMeganetElement(L,nex=nex)
    end

    nImg   = [8 4]
    nc     = 3
    nex    = 4
    K      = getSparseConvKernel2D(TYPE,nImg,[3,3,1,nc])
    Bin    = randn(TYPE,nFeatOut(K),4)
    Bout   = randn(TYPE,nFeatIn(K),3)
    nLayer = getBatchNormLayer(TYPE,[prod(nImg),nc])
    L      = getDoubleSymLayer(TYPE,K,nLayer,Bin=Bin,Bout=Bout)
    @testset "doubleSymLayer (conv/Batch) $TYPE" begin
     testAbstractMeganetElement(L,nex=nex)
    end

    nImg   = [16 8]
    nc     = 6
    nex    = 8
    K      = getSparseConvKernel2D(TYPE,nImg,[3,3,1,nc])
    Bin    = randn(TYPE,nFeatOut(K),4)
    Bout   = randn(TYPE,nFeatIn(K),3)
    nLayer = getTVNormLayer(TYPE,[prod(nImg),nc])
    L      = getDoubleSymLayer(TYPE,K,nLayer,Bin=Bin,Bout=Bout)
    @testset "doubleSymLayer (conv/TV) $TYPE" begin
     testAbstractMeganetElement(L)
    end
end


text:
import Base: +, -, ==, <, <=, hash

struct Priority<:Real
    time::Float64
    priority::UInt
end

"""
    PriorityTime(time)
Construct `PriorityTime` with `priority` randomly generated using `rand()`.
"""

Priority(time::Float64) = Priority(time, rand(UInt))

function +(x::Priority, y::Priority)
    Priority(x.time+y.time, x.priority+y.priority)
end

function -(x::Priority, y::Priority)
    Priority(x.time-y.time, x.priority-y.priority)
end

function ==(x::Priority, y::Priority)
    x.time==y.time && x.priority==y.priority
end

function <(x::Priority, y::Priority)
    x.time < y.time && return true
    x.time == y.time && x.priority < y.priority && return true
    return false
end

function <=(x::Priority, y::Priority)
    x.time < y.time && return true
    x.time == y.time && x.priority <= y.priority && return true
    return false
end

function hash(x::Priority, h::UInt)
    hash((x.time, x.priority), h)
end


text:
"""
    crosscorrelate(series, template, [element_type=Float64], normalize_template=true)

Return the normalized cross-correlation between `series` and `template`.

If `series` is a vector of ``n`` elements and template a vector of ``N`` elements,
cross-correlation will be a vector of ``n - N + 1`` elements of type `element_type`.
The element of cross-correlation ``χ_k`` is Pearson correlation coefficient between
``x^{(k)}_i`` denoting the elements of `series` ``x_{i + k}`` and
the elements of `template` ``y_i``, with ``i = 1, 2, \\dots, N``

``
χ_k = \\frac{\\sum{\\left(x^{(k)}_i - μ_{x^{(k)}}\\right)
\\left(y_i - μ_y\\right)}}{N σ_{x^{(k)}} σ_y} \\, ,
``

where ``μ`` and ``σ`` denotes the mean and variance respectively.

If `normalize_template` is set to false then `template` is assumed to have mean and std
respectively equal to 0.0 and 1.0.

See the [Wikipedia page](https://en.wikipedia.org/wiki/Cross-correlation#Normalization)
on cross-correlation for more details.

# Examples

```jldoctest
julia> crosscorrelate(sin.(0:0.25pi:2pi), [1, 1+√2, 1])
7-element Vector{Float64}:
  0.23258781949447394
  1.000000000000001
  0.23258781949447402
  7.401486830834377e-17
 -0.23258781949447394
 -1.0000000000000007
 -0.23258781949447394
```
"""
function crosscorrelate(series::AbstractVector{T1}, template::AbstractVector{T2}, element_type::Type{T3}=Float64;
    normalize_template=true) where {T1 <: Number, T2 <: Number, T3 <: AbstractFloat}

    if isempty(series)
        throw(ArgumentError("Series must be a non-empty vector."))
    elseif isempty(template)
        throw(ArgumentError("Template must be a non-empty vector."))
    elseif size(series) < size(template)
        throw(DimensionMismatch("Template is longer than series."))
    end

    if normalize_template
        template_mean, template_std = mean_and_std(template, corrected=false)
        y = @. (template - template_mean) / template_std
    else
        y = template
    end

    N = length(template)
    cc = Vector{element_type}(undef, length(series) - N + 1)
    x = similar(y)
    for n in eachindex(cc)
        series_view = view(series, n:n + N - 1)
        series_view_std = std(series_view, corrected=false)
        @. x = series_view / series_view_std
        cc[n] = dot(x, y) / N
    end
    cc
end

"""
    maxfilter(x, tolerance)

Return a vector of the same size of `x` whose element at index `n` is the maximum
of the elements of `x` which are at most `tolerance` apart from `x[n]`
(i.e. in a window of at most `2 * tolerance + 1` elements).

# Examples

```jldoctest
julia> maxfilter(sin.(0:0.25pi:2pi), 1)
9-element Vector{Float64}:
  0.7071067811865475
  1.0
  1.0
  1.0
  0.7071067811865476
  1.2246467991473532e-16
 -0.7071067811865475
 -2.4492935982947064e-16
 -2.4492935982947064e-16
```
"""
function maxfilter(x::AbstractVector, l)
    y = similar(x)
    for n in eachindex(x)
        lower = max(n - l, firstindex(x))
        upper = min(n + l, lastindex(x))
        y[n] = maximum(view(x, lower:upper))
    end
    y
end

"""
    stack(correlations, offsets)

Return the average cross-correlation after aligning `correlations`.
Each cross-correlation `correlations[n]` is shifted to the right by `offsets[n]`.

# Examples

```jldoctest
julia> stack([[0, 1.0, 0, 0], [0, 0, 1.0, 0]], [1, 2])
3-element OffsetArray(::Vector{Float64}, 0:2) with eltype Float64 with indices 0:2:
 0.0
 1.0
 0.0
```
"""
function stack(correlations::AbstractVector{T1}, offsets::AbstractVector{T2}) where {T1 <: AbstractVector, T2 <: Integer}
    if isempty(correlations)
        throw(ArgumentError("Cross-correlations vector must be non-empty."))
    elseif length(correlations) != length(offsets)
        throw(DimensionMismatch("Cross-correlations and offsets vectors must have the same length."))
    end
    stackedcorrelations = OffsetVector.(correlations, -offsets)
    start = maximum(series -> firstindex(series), stackedcorrelations)
    stop = minimum(series -> lastindex(series), stackedcorrelations)
    OffsetVector(mean(series -> view(series, start:stop), stackedcorrelations), start:stop)
end

"""
    correlatetemplate(data, template, offsets, tolerance, [element_type=Float64])

Return the cross-correlation between `data` and `template`. The cross-correlations for each series
are aligned using `offsets` and averaged. If `tolerance`` is not zero, then the average accounts for possible
misplacement of each series by `tolerance` sample, and return the average of the maximum cross-correlation
compatible with that misplacement.
"""
function correlatetemplate(data, template, offsets, tolerance, element_type=Float64)
    if isempty(data)
        throw(ArgumentError("Data must be non-empty."))
    elseif !(length(data) == length(template) == length(offsets))
        throw(DimensionMismatch("Data, template and shifts must have the same length."))
    end
    correlations = similar(data, Vector{element_type})
    Threads.@threads for n = eachindex(data)
        correlations[n] = maxfilter(crosscorrelate(data[n], template[n], element_type), tolerance)
    end
    stack(correlations, offsets)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdtransfer

"""
    'feffery'_antdtransfer(;kwargs...)

An AntdTransfer component.

Keyword arguments:
- `id` (String; optional)
- `className` (String; optional)
- `dataSource` (optional): . dataSource has the following type: Array of lists containing elements 'key', 'title'.
Those elements have the following types:
  - `key` (String | Real; optional)
  - `title` (String; optional)s
- `disabled` (Bool; optional)
- `height` (String; optional)
- `moveDirection` (String; optional)
- `moveKeys` (Array; optional)
- `operations` (Array; optional)
- `pagination` (optional): . pagination has the following type: Bool | lists containing elements 'pageSize'.
Those elements have the following types:
  - `pageSize` (Real; optional)
- `showSearch` (Bool; optional)
- `showSelectAll` (Bool; optional)
- `style` (Dict; optional)
- `targetKeys` (Array; optional)
- `titles` (Array; optional)
"""
function 'feffery'_antdtransfer(; kwargs...)
        available_props = Symbol[:id, :className, :dataSource, :disabled, :height, :moveDirection, :moveKeys, :operations, :pagination, :showSearch, :showSelectAll, :style, :targetKeys, :titles]
        wild_props = Symbol[]
        return Component("'feffery'_antdtransfer", "AntdTransfer", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: Beddow/feffery-antd-components

max_issues_repo_name: limmyforget/feffery-antd-components

max_forks_repo_name: limmyforget/feffery-antd-components


text:
using BinaryBuilder

# Collection of sources required to build FastTransforms
name = "FastTransforms"
version = v"0.3.2"
sources = [
    ArchiveSource("https://github.com/MikaelSlevinsky/FastTransforms/archive/v$(version).tar.gz",
                  "ef401079751e570d3dcfd4ac3c40e7d7d5858a44ea66c62808cf319a68c223c8"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/FastTransforms-*
if [[ ${target} == x86_64-w64-mingw32 ]]; then
    export CFLAGS="-O3 -mavx -fno-asynchronous-unwind-tables "
else
    export CFLAGS="-O3 -mavx "
fi
if [[ ${nbits} == 64 ]]; then
    SYMBOL_DEFS=()
    SYMBOLS=(dgemm dtrmm dtrmv dtrsm sgemm strmm strsm ztrmm)
    for sym in ${SYMBOLS[@]}; do
        SYMBOL_DEFS+=("-Dcblas_${sym}=cblas_${sym}64_")
    done
    CFLAGS+=${SYMBOL_DEFS[@]}
    BLAS=openblas64_
else
    BLAS=openblas
fi
make assembly CC=gcc
make lib CC=gcc FT_PREFIX=${prefix} FT_BLAS=${BLAS} FT_FFTW_WITH_COMBINED_THREADS=1
mv -f libfasttransforms.${dlext} ${libdir}
"""

platforms = expand_gfortran_versions(supported_platforms())
platforms = [p for p in platforms if BinaryBuilder.proc_family(p) == :intel]

# The products that we will ensure are always built
products = [
    LibraryProduct("libfasttransforms", :libfasttransforms),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("FFTW_jll"),
    Dependency("MPFR_jll"),
    Dependency("OpenBLAS_jll"),
]

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


text:
module TryMultiDispatch
import Base.>=, Base.push!
using PatternDispatch.Patterns, PatternDispatch.Recode, PatternDispatch.Dispatch

abstract DNode

type Decision <: DNode
    intent::Intension
    pass
    fail
    seq::Vector

    Decision(intent::Intension, pass, fail) = new(intent, pass, fail)
end

type Method
    sig::Pattern
    body
end

type MethodCall <: DNode
    m::Method
    bind_seq::Vector
    bindings::Dict{Symbol,Node}

    MethodCall(m::Method) = new(m)
end

type MethodNode
    m
    gt::Set{MethodNode}
end
MethodNode(m) = MethodNode(m, Set{MethodNode}())

>=(x::MethodNode, y::MethodNode) = intentof(x) >= intentof(y)

type NoMethod <: DNode; end
const nomethod = NoMethod()

intentof(m::Method)     = m.sig.intent
intentof(m::MethodNode) = intentof(m.m)
intentof(::NoMethod)    = anything

type MethodTable
    name::Symbol
    top::MethodNode
    f::Function

    function MethodTable(name::Symbol)
        mt = new(name, MethodNode(nomethod))
        mt.f = rebuilder(mt)
        mt
    end
end

rebuilder(mt::MethodTable) = (args...)->begin
        mt.f = create_dispatch(mt)
        mt.f(args...)
    end

dispatch(mt::MethodTable, args) = mt.f(args...)

function push!(mt::MethodTable, m::Method)
    insert!(mt.top, MethodNode(m))
    mt.f = rebuilder(mt)
end

function code_dispatch(mt::MethodTable)
    dtree = build_dtree(mt.top, subtreeof(mt.top))

    seq_dispatch!(ResultsDict(), dtree)
    code = code_dispatch(dtree)
    fdef = :(($argsym...)->$code)
end
create_dispatch(mt::MethodTable) = eval(code_dispatch(mt)) # todo: which eval?

# ---- update method DAG ------------------------------------------------------

insert!(at::MethodNode, m::MethodNode) = insert!(Set{MethodNode}(), at, m)
function insert!(seen::Set{MethodNode}, at::MethodNode, m::MethodNode)
    if has(seen, at); return; end
    push!(seen, at)
    if m >= at
        if at >= m
            at.m = m.m  # at == m
            return true
        end
        # m > at
        push!(m.gt, at)
        at_above_m = false
    else
        at_above_m = any([insert!(seen, below, m) for below in at.gt])
    end
    if !at_above_m
        if at >= m
            del_each(at.gt, m.gt)
            push!(at.gt, m)
            at_above_m = true
        end
    end
    at_above_m
end

# ---- create decision tree ---------------------------------------------------

firstitem(iter) = next(iter, start(iter))[1]

subtreeof(m::MethodNode) = (sub = Set{MethodNode}(); addsubtree!(sub, m); sub)
function addsubtree!(seen::Set{MethodNode}, m::MethodNode)
    if has(seen, m); return; end
    push!(seen, m)
    for below in m.gt; addsubtree!(seen, below); end
end

function build_dtree(top::MethodNode, ms::Set{MethodNode})
    if isempty(top.gt) || length(ms) == 1
        isa(top.m, NoMethod) ? nomethod : MethodCall(top.m)
    else
        pivot = firstitem(top.gt & ms)
        below = subtreeof(pivot)
        pass = build_dtree(pivot, ms & below)
        fail = build_dtree(top,   ms - below)
        Decision(intentof(pivot), pass, fail)
    end
end

# ---- seq_dispatch! ----------------------------------------------------------

seq_dispatch!(results::ResultsDict, d::DNode) = nothing
function seq_dispatch!(results::ResultsDict, m::MethodCall)
    s = Sequence(m.m.sig.intent, results) # shouldn't need an Intension...
    for node in values(m.m.sig.bindings);  sequence!(s, node)  end
    m.bind_seq = s.seq
    m.bindings = (Symbol=>Node)[name => results[node]
                                for (name,node) in m.m.sig.bindings]
end
function seq_dispatch!(results::ResultsDict, d::Decision)
    results_fail = copy(results)
    s = Sequence(d.intent, results)
    for g in guardsof(d.intent);  sequence!(s, Guard(g))  end
    d.seq = s.seq

    seq_dispatch!(results, d.pass)
    seq_dispatch!(results_fail, d.fail)
end

function code_dispatch(::NoMethod)
    :( error("No matching pattern method found") )
end
function code_dispatch(m::MethodCall)
    prebind = encoded(m.bind_seq)
    binds = { :( $name = $(resultof(node))) for (name,node) in m.bindings }
    quote
        $(prebind...)
        $(binds...)
        $(m.m.body)
    end
end
function code_dispatch(d::Decision)
    pred = code_predicate(d.seq)
    pass = code_dispatch(d.pass)
    fail = code_dispatch(d.fail)
    code = :( if $pred; $pass; else; $fail; end )
    #length(pre) == 0 ? code : (quote; $(pre...); $code; end)
    code
end

# ---- test code ----

m1 = Method((@qpat (x::Int,)),    :1)
m2 = Method((@qpat (x::String,)), :2)
m3 = Method((@qpat (x,)),         :3)

# node1 = MethodNode(m1, Set{MethodNode}())
# node2 = MethodNode(m2, Set{MethodNode}())
# node3 = MethodNode(m3, Set{MethodNode}(node1, node2))
# top   = MethodNode(nomethod, Set{MethodNode}(node3))
# mnodes = Set{MethodNode}(top, node1, node2, node3)

# dc = Decision(m2.sig.intent, MethodCall(m2), MethodCall(m3))
# db = Decision(m1.sig.intent, MethodCall(m1), dc)
# da = Decision(m3.sig.intent, db, nomethod)

# top = MethodNode(nomethod)
# for m in [m3,m2,m1];  insert!(top, MethodNode(m))  end
# mnodes = subtreeof(top)

# dtree = build_dtree(top, mnodes)

# seq_dispatch!(ResultsDict(), dtree)
# code = code_dispatch(dtree)
# println(code)
# fdef = :(($argsym...)->$code)

# f = eval(fdef)

mt = MethodTable(:f)
for m in [m3,m2,m1];  push!(mt, m)  end

@assert mt.f(5)     == 1
@assert mt.f("foo") == 2
@assert mt.f(5.0)   == 3

end # module


text:
###
### This test is to check against the reference zarr implementation in Python
### We save some data in Julia and python and test if it is still the same
### when read from teh other language
###
@testset "Python zarr implementation" begin

using PyCall
import PyCall: @py_str
#If we are on conda, import zarr
pyimport_conda("zarr","zarr")

#Create some directories
proot = tempname()
mkpath(proot)
pjulia = joinpath(proot,"julia")
ppython = joinpath(proot,"python")

#First create an array in Julia and read with python zarr
groupattrs = Dict("String attribute"=>"One", "Int attribute"=>5, "Float attribute"=>10.5)
g = zgroup(pjulia,attrs=groupattrs)

# Test all supported data types and compressors
import Zarr: NoCompressor, BloscCompressor, ZlibCompressor, MaxLengthString
using Random: randstring
dtypes = (UInt8, UInt16, UInt32, UInt64,
    Int8, Int16, Int32, Int64,
    Float16, Float32, Float64,
    Complex{Float32}, Complex{Float64},
    Bool,MaxLengthString{10,UInt8},MaxLengthString{10,UInt32})
compressors = ("no"=>NoCompressor(), "blosc"=>BloscCompressor(cname="zstd"),"zlib"=>ZlibCompressor())
testarrays = Dict(t=>(t<:AbstractString) ? [randstring(maximum(i.I)) for i in CartesianIndices((1:10,1:6,1:2))] : rand(t,10,6,2) for t in dtypes)

for t in dtypes, co in compressors
    compstr, comp = co
    att = Dict("This is a nested attribute"=>Dict("a"=>5))
    a = zcreate(t, g,string("a",t,compstr),10,6,2,attrs=att, chunks = (5,2,2),compressor=comp)
    a[:,:,:] = testarrays[t]
end
# Test reading in python
py"""
import zarr
g = zarr.open_group($pjulia)
gatts = g.attrs
"""

#Test group attributes
@test py"gatts['String attribute']" == "One"
@test py"gatts['Int attribute']" == 5
@test py"gatts['Float attribute']" == 10.5

dtypesp = ("uint8","uint16","uint32","uint64",
    "int8","int16","int32","int64",
    "float16","float32","float64",
    "complex64", "complex128","bool","S10","U10")

#Test accessing arrays from python and reading data
for i=1:length(dtypes), co in compressors
    compstr,comp = co
    t = dtypes[i]
    tp = dtypesp[i]
    arname = string("a",t,compstr)
    py"""
    ar=g[$arname]
    """

    @test py"ar.attrs['This is a nested attribute']" == Dict("a"=>5)
    @test py"ar.dtype==$tp"
    @test py"ar.shape" == (2,6,10)
    if t<:MaxLengthString
      pyar = py"ar[:,:,:]"
      jar = [get(get(get(pyar,k-1),j-1),i-1) for i in 1:10, j in 1:6, k in 1:2]
      @test jar == testarrays[t]
    else
      @test py"ar[:,:,:]" == permutedims(testarrays[t],(3,2,1))
    end
end

## Now the other way around, we create a zarr array using the python lib and read back into julia
data = rand(Int32,2,6,10)
py"""
import numcodecs
g = zarr.group($ppython)
g.attrs["groupatt"] = "Hi"
z1 = g.create_dataset("a1", shape=(2,6,10),chunks=(1,2,3), dtype='i4')
z1[:,:,:]=$data
z1.attrs["test"]={"b": 6}
z2 = g.create_dataset("a2", shape=(5,),chunks=(5,), dtype='S1', compressor=numcodecs.Zlib())
z2[:]=[k for k in 'hallo']
zarr.consolidate_metadata($ppython)
"""

#Open in Julia
g = zopen(ppython)
@test g isa Zarr.ZGroup
@test g.attrs["groupatt"] == "Hi"
a1 = g["a1"]
@test a1 isa ZArray
@test a1[:,:,:]==permutedims(data,(3,2,1))
@test a1.attrs["test"]==Dict("b"=>6)
# Test reading the string array
@test String(g["a2"][:])=="hallo"

# And test for consolidated metadata
# Delete files so we make sure they are not accessed
rm(joinpath(ppython,".zattrs"))
rm(joinpath(ppython,"a1",".zattrs"))
rm(joinpath(ppython,"a1",".zarray"))
rm(joinpath(ppython,"a2",".zarray"))
g = zopen(ppython, "w", consolidated=true)
@test g isa Zarr.ZGroup
@test g.attrs["groupatt"] == "Hi"
a1 = g["a1"]
@test a1 isa ZArray
@test a1[:,:,:]==permutedims(data,(3,2,1))
@test a1.attrs["test"]==Dict("b"=>6)
@test storagesize(a1) == 960
@test sort(Zarr.subkeys(a1.storage,"a1"))[1:5] == ["0.0.0","0.0.1","0.0.2","0.0.3","0.1.0"]
a1[:,1,1] = 1:10
@test a1[:,1,1] == 1:10
# Test reading the string array
@test String(g["a2"][:])=="hallo"
end

@testset "Python datetime types" begin
using Dates, Test, Zarr
vd = Date(1970,1,1):Day(1):Date(1970,6,30) |> collect
vt = DateTime(1970,1,1):Second(1):DateTime(1970,1,1,2,0,0)|> collect
ad = ZArray(vd)
at = ZArray(vt)
@test eltype(ad)==Zarr.DateTime64{Day}
@test eltype(at)==Zarr.DateTime64{Millisecond}
@test DateTime.(at[:]) == vt[:]
@test Date.(ad[:]) == vd[:]

p = tempname()
g = zgroup(p)
for pt in [Week, Day, Hour, Minute, Second,
        Millisecond]

    if pt <: DatePeriod
        vd = range(Date(1970,1,1),step = pt(1), length=100)
        a = zcreate(Zarr.DateTime64{pt},g,string(pt),100)
        a[:] = vd
    else
        vd = range(DateTime(1970,1,1),step = pt(1), length=100)
        a = zcreate(Zarr.DateTime64{pt},g,string(pt),100)
        a[:] = vd
    end
end

zarr = pyimport("zarr")
np = pyimport("numpy")

g_julia = zopen(p)
g_python = zarr.open(p)

for unit in ["Week", "Day", "Hour", "Minute", "Second",
        "Millisecond"]
    @test_py np.datetime64(g_julia[unit][1] |> DateTime |> string) == get(getproperty(g_python,unit),0)
    @test_py np.datetime64(g_julia[unit][10] |> DateTime |> string) ==  get(getproperty(g_python,unit),9)
    @test_py np.datetime64(g_julia[unit][100] |> DateTime |> string) == get(getproperty(g_python,unit),99)
end

end


text:
#
# Calculates requests per second from client to server without response, but
# using only naked sockets.
# Results: rps: 13160.318995561785, time: 75.98600006103516
#
type T
  x::Int
end
server = listen(2000)
sock   = accept(server)
function run()
  local t::T = T(rand(1:100))
  for i=1:1000000 serialize(sock, t) end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo