text:
# nanmin.jl
#
#   Minimum, ignoring NaNs
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.6.0
# 02.04.2020
# Last Edit: 08.04.2021

nanmin(x::Array{T}) where T <: Number = minimum(filter(!isnan,x))

nanmin(x::Array{Union{T,Missing}}) where T <: Number = minimum(filter(!isnan,filter(!ismissing,x)))

function nanmin(x::Matrix{T},y::Integer) where T Union{S, Union{S, Missing}} where S <: Number
  return mapslices(nanmin,x,dims=y)
end


text:
import Base.+

function +(buffer::Array{String,1}, line::String)
    push!(buffer, line)
end

function read_reaction_file(path_to_file::String)::Array{String,1}

    # initialize -
    vff_file_buffer = String[]
    vff_reaction_array = Array{String,1}()

    # Read in the file -
    open("$(path_to_file)", "r") do file
        for line in eachline(file)
            +(vff_file_buffer,line)
        end
    end

    # process -
    for reaction_line ∈ vff_file_buffer

        # skip comments and empty lines -
        if (occursin("//", reaction_line) == false &&
            isempty(reaction_line) == false)

            # grab -
            push!(vff_reaction_array,reaction_line)
        end
    end

    # return -
    return vff_reaction_array
end

ext: jl


text:
function solve()
    n = parse(Int, readline())
    m = floor(Int, √n)

    ds = Set{Int}()

    for x=1:m
        if n % x == 0
            push!(ds, x, n / x)
        end
    end
    ds = sort(collect(ds))
    for d in ds
        println(d)
    end
end

solve()

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
module C2Julia

using Clang
using Clang.LibClang
import Clang: TokenList

using CSyntax

function translate(cursor::CLCursor)
    expr = Expr(:macrocall, Symbol("@error"), nothing, "not implemented yet or unknown bugs")
    return MetaExpr(expr, cursor)
end

include("types.jl")
export MetaExpr

include("literal.jl")
include("operator.jl")
include("expr.jl")
include("stmt.jl")
include("decl.jl")

export translate

end # module


text:
"""
    sinewave1D!(height, h₀, n, ϵ)

Creates a sine wave like height field with `n` full waves, average height `h₀` and displacement magnitude `ϵ`.
"""
function sinewave1D!(height, h₀, n::Int, ϵ, dims::Int)
    # Some helping dummy array
    hdummy = ones(size(height))
    # Used to perform the accumulate operation
    height .= 1.0
    # This we get n nice sine waves
    hdummy .= 2π .* n .* accumulate(+, height, dims=dims) ./ (size(height,dims)-1)
    # Still need to use the correct undulation ϵ and initial height
    height .= h₀ .* (1.0 .+ ϵ .* (map(sin, hdummy) .- 1/size(height,dims)))
    return nothing
end

"""
    randinterface!(height, h₀, ϵ)

Creates a random height field with average height `h₀` and displacement magnitude `ϵ`.
"""
function randinterface!(height, h₀, ϵ)
    # Some helping dummy array
    hdummy = ones(size(height))
    # Used to perform the accumulate operation
    height .= 1.0
    # This we get n nice sine waves
    randn!(hdummy)
    # Still need to use the correct undulation ϵ and initial height
    height .= h₀ .* (1.0 .+ ϵ .* hdummy)
    return nothing
end

"""
    singledroplet(T, size(θ,1), size(θ,2), radius, θ, center, device=false)

Generates a fluid configuration of a single droplet in the shape of spherical cap with contact angle `θ`, sphere radius `radius` and centered at `center`.

# Arguments

- `height::Array{undef, 2}`: numerical formate, either `Float64` or `Float32`
- `radius::AbstractFloat`: radius of the underlying sphere from which the spherical cap is cut off
- `θ::AbstractFloat`: contact angle in multiples of `π`
- `center::Tuple{Int,Int}`: x and y coordinates of the center of the droplet

# Examples

```jldoctest
julia> using Swalbe, Test

julia> rad = 50; θ = 1/3;

julia> height = Swalbe.singledroplet(ones(100,100), rad, θ, (50,50));

julia> @test maximum(height) == rad * (1 - cospi(θ)) # Simple geometry
Test Passed
  Expression: maximum(height) == rad * (1 - cospi(θ))
   Evaluated: 25.0 == 25.0

julia> argmax(height) # Which is constistent with the center!
CartesianIndex(50, 50)

```

# References

See also:
"""
function singledroplet(height, radius, θ, center)
    lx, ly = size(height)
    # area = 2π * radius^2 * (1- cospi(θ))
    @inbounds for i in 1:lx
        for j in 1:ly
            circ = sqrt((i-center[1])^2 + (j-center[2])^2)
            if circ <= radius
                height[i,j] = (cos(asin(circ/radius)) - cospi(θ)) * radius
            else
                height[i,j] = 0.05
            end
        end
    end
    @inbounds for i in 1:lx
        for j in 1:ly
            if height[i,j] < 0
                height[i,j] = 0.05
            end
        end
    end
    return height
end
# For a one dimensional droplet
function singledroplet(height::Vector, radius, θ, center)
    L = size(height)[1]
    # area = 2π * radius^2 * (1- cospi(θ))
    @inbounds for i in 1:L
        circ = sqrt((i-center)^2)
        if circ <= radius
            height[i] = (cos(asin(circ/radius)) - cospi(θ)) * radius
        else
            height[i] = 0.05
        end
    end
    @inbounds for i in 1:L
        if height[i] < 0
            height[i] = 0.05
        end
    end
    return height
end

"""
    two_droplets(sys)

Generates a fluid configuration of a two droplets in the shape of spherical cap with contact angles `θ₁`, `θ₂`, sphere radius `r₁`, `r₂` and centers at `center`.

This is work in progress, therefore so far it is only available for the lower dimension model.

# Arguments

- `r₁::Float64`: radius of the first sphere from which the cap is cut of
- `r₂::Float64`: radius of the second sphere from which the cap is cut of
- `θ₁::Float64`: contact angle of the first cap in multiples of `π`
- `θ₂::Float64`: contact angle of the second cap in multiples of `π`
- `center::Tuple{Int,Int}`: x coordinates of the centers of the droplets

# Examples

```jldoctest
julia> using Swalbe, Test

julia> rad = 45; θ = 1/4; sys = Swalbe.SysConst_1D(L=200);

julia> height = Swalbe.two_droplets(sys, r₁=rad, r₂=rad, θ₁=θ, θ₂=θ);

julia> @test maximum(height) ≈ rad * (1 - cospi(θ)) atol=0.01 # Simple geometry
Test Passed
  Expression: ≈(maximum(height), rad * (1 - cospi(θ)), atol = 0.01)
   Evaluated: ≈(13.178960261768571, 13.180194846605358; atol = 0.01)

```

# References

See also:
"""
function two_droplets(sys::SysConst_1D; r₁=230, r₂=230, θ₁=1/9, θ₂=1/9, center=(sys.L/3,2*sys.L/3))
    dum = zeros(sys.L)
    dum2 = zeros(sys.L)
    height = zeros(sys.L)
    # area = 2π * radius^2 * (1- cospi(θ))
    @inbounds for i in 1:sys.L
        circ = sqrt((i-center[1])^2)
        if circ <= r₁
            dum[i] = (cos(asin(circ/r₁)) - cospi(θ₁)) * r₁
        else
            dum[i] = 0.05
        end
    end

    @inbounds for i in 1:sys.L
        circ2 = sqrt((i-center[2])^2)
        if circ2 <= r₂
            dum2[i] = (cos(asin(circ2/r₂)) - cospi(θ₂)) * r₂
        else
            dum2[i] = 0.05
        end
    end

    @inbounds for i in 1:sys.L
        if dum[i] < 0
            dum[i] = 0.05
        end
        if dum2[i] < 0
            dum2[i] = 0.05
        end
    end

    height .= dum .+ dum2 .- 0.05
    @inbounds for i in 1:sys.L
        if height[i] < 0
            height[i] = 0.05
        end
    end
    return height
end

"""
    restart_from_height(data)

Restarts a simulation from already generated height data.

# Arguments

- `data::file`: Some file with computed data
- `kind::String`: File format, `.bson` or `.jld2` are valid options
- `timestep::Int`: The time step the height is read from file
- `size::Tuple{Int,Int}`: x and y limits of the computational domain

```julia
julia> using Swalbe, FileIO, Test

julia> h1 = rand(10,10); h2 = rand(10,10);

julia> f = Dict() # For storage
Dict{Any, Any}()

julia> f["h_1"] = vec(h1); f["h_2"] = vec(h2);

julia> save("file.jld2", f)

julia> h = Swalbe.restart_from_height("file.jld2", timestep=1, size=(10,10));

julia> @test all(h .== h1)
Test Passed

julia> rm("file.jld2")

```
"""
function restart_from_height(data; kind="jld2", timestep=0, size=(512,512))
    # Allocate the height
    height = zeros(size)
    # Use the correct file format
    if kind == "jld2"
        # Pipe the data to a dataframe
        df = load(data) |> DataFrame
        if timestep == 0
            height .= reshape(df[:,end], size[1], size[2])
        else
            height .= reshape(df[!, Symbol("h_$(timestep)")], size[1], size[2])
        end
    elseif kind == "bson"
        df = DataFrame(BSON.load(data))
        # Now choose only the timestep that is needed
        if timestep == 0
            height .= reshape(df[:,end], size[1], size[2])
        else
            height .= reshape(df[!, Symbol("h_$(timestep)")], size[1], size[2])
        end
    end

    return height

end

# ---------------------------------------------------------------- #
# ---------------------- Substrate patterns ---------------------- #
# ---------------------------------------------------------------- #
"""
    boxpattern(θ, θ₀; center=(size(θ,1)÷2, size(θ,2)÷2), δₐ=1/36, side=20)

Defines a quadratic box around a `center` with side length `side`.

# Arguments

- `θ::Array{undef,2}`: Contact angle map
- `θ₀`: Default contact angle `θ`
- `center::Tuple{Int, Int}`: center of the box pattern, default value is `(size(θ,1)÷2, size(θ,2)÷2)`
- `δₐ::Float64`: Contact angle contrast with the substrate, default is set to `1/36` ≈ 5 degrees difference
- `side::Int`: Length of the sides, default is set to `20`

# Examples

```jldoctest
julia> using Swalbe

julia> θ₀ = 1/9;

julia> pattern, polygon = Swalbe.boxpattern(ones(100,100), θ₀);

julia> polygon # Some cool thing we use to create the posize(θ,2)gones, a LazySet
LazySets.VPolygon{Float64, Vector{Float64}}([[40.0, 40.0], [60.0, 40.0], [60.0, 60.0], [40.0, 60.0]])

julia> pattern[50,50] == pattern[1,1] # In the center there is a different contact angle!
false

```

# References

Really not much to say here, check out [LazySets.jl](https://github.com/JuliaReach/LazySets.jl).

"""
function boxpattern(θ, θ₀; center=(size(θ,1)÷2, size(θ,2)÷2), δₐ=1/36, side=20)
    # Vertex vectors
    xcoords = [center[1]-side/2 center[1]-side/2 center[1]+side/2 center[1]+side/2]
    ycoords = [center[2]-side/2 center[2]+side/2 center[2]-side/2 center[2]+side/2]
    # Check if the vertices are on the grid
    for i in 1:length(xcoords)
        if xcoords[i] < 1 || xcoords[i] > size(θ,1)
            xcoords[i] = 1 # If not push them to the origin
        elseif ycoords[i] < 1 || ycoords[i] > size(θ,2)
            ycoords[i] = 1
        end
    end
    # Vertex coordinates
    vertex₁ = [xcoords[1], ycoords[1]]
    vertex₂ = [xcoords[2], ycoords[2]]
    vertex₃ = [xcoords[3], ycoords[3]]
    vertex₄ = [xcoords[4], ycoords[4]]
    # Build a posize(θ,2)gon
    P = VPolygon([vertex₁, vertex₂, vertex₃, vertex₄])
    for i in 1:size(θ,1), j in 1:size(θ,2)
        if [Float64(i), Float64(j)] ∈ P
            θ[j,i] = θ₀ + δₐ
        else
            θ[j,i] = θ₀
        end
    end
    return θ, P
end

"""
    ellipsepattern(θ, θ₀; center=(size(θ,1)÷2, size(θ,2)÷2), δₐ=1/36, a=10, b=5)

Creates an ellipse shaped contact angle defect with contact angle mismatch `δₐ`.

# Arguments

- `θ::Array{undef,2}`: Contact angle map
- `θ₀`: Default contact angle `θ`
- `center::Tuple{Int, Int}`: center of the created pattern, default values is `center = (size(θ,1)÷2, size(θ,2)÷2)`
- `δₐ::Float64`: contact angle mismatch between patch and rest of substrate, default is `δₐ = 1/36` or 5 degrees
- `a::Int`:: semimajor half ax of the ellipse, default value is `a=10`
- `b::Int`:: semiminor half ax of the ellipse, default value is `b=5`

# Examples

```jldoctest
julia> using Swalbe, Test

julia> θ₀ = 1/9;

julia> θ, P = Swalbe.ellipsepattern(ones(100,100), θ₀); # per default the center is in the middle!

julia> @test θ[1,1] == θ₀
Test Passed
  Expression: θ[1, 1] == θ₀
   Evaluated: 0.1111111111111111 == 0.1111111111111111

julia> @test θ[50,50] == θ₀ + 1/36 # The default increment, is about 5 degrees.
Test Passed
  Expression: θ[50, 50] == θ₀ + 1 / 36
   Evaluated: 0.1388888888888889 == 0.1388888888888889
```

# References
Nothing interesting here.

"""
function ellipsepattern(θ, θ₀; center=(size(θ,1)÷2, size(θ,2)÷2), δₐ=1/36, a=10, b=5)
    # Don't know why LazySets needs this, can not just parse center as argument of Ellipsoid constructor
    mid = ones(2)
    mid[1] = center[1]
    mid[2] = center[2]
    # The eigenvalues of the shape matrix, Q in the LazySets docs, needs to be spezified not the axis length
    shape₁ = a^(2)
    shape₂ = b^(2)
    P = Ellipsoid(mid, [shape₁ 0.0; 0.0 shape₂])
    for i=1:size(θ,1), j=1:size(θ,2)
        if [Float64(i), Float64(j)] ∈ P
            θ[j,i] = θ₀ + δₐ
        else
            θ[j,i] = θ₀
        end
    end
    return θ, P
end

"""
    trianglepattern(θ, θ₀; center=(size(θ,1)÷2, size(θ,2)÷2), δₐ=1/36, side=60)

Generates an equilateral triangle centered around `center` with contact angle contrast `δₐ` and side length `side`.

# Arguments
- `θ::Array{undef,2}`: Contact angle map
- `θ₀`: Default contact angle `θ`
- `center::Tuple(Int, Int)`: position of the center of the triangle
- `δₐ::Float64`: contact angle contrast with the rest of the substrate
- `side::Int`: length of the sides of the equilateral triangle

# Examples
```jldoctest
julia> using Swalbe

julia> θ, P = Swalbe.trianglepattern(ones(50,50), 1/9, side=20) # Returns a polygon and the contact angle field
([0.1111111111111111 0.1111111111111111 … 0.1111111111111111 0.1111111111111111; 0.1111111111111111 0.1111111111111111 … 0.1111111111111111 0.1111111111111111; … ; 0.1111111111111111 0.1111111111111111 … 0.1111111111111111 0.1111111111111111; 0.1111111111111111 0.1111111111111111 … 0.1111111111111111 0.1111111111111111], LazySets.VPolygon{Float64, Vector{Float64}}([[15.0, 19.226497308103742], [35.0, 19.226497308103742], [25.0, 36.547005383792516]]))

julia> P
LazySets.VPolygon{Float64, Vector{Float64}}([[15.0, 19.226497308103742], [35.0, 19.226497308103742], [25.0, 36.547005383792516]])

julia> θ[25,25]
0.1388888888888889
```

# References

See also:
"""
function trianglepattern(θ, θ₀; center=(size(θ,1)÷2, size(θ,2)÷2), δₐ=1/36, side=60)
    height = sqrt(3)/2*side
    # Vertices vectors
    xcoords = [center[1]-side/2 center[1]+side/2 center[1]]
    ycoords = [center[2]-height/3 center[2]-height/3 center[2]+2*height/3]
    # Check if the vertices are on the grid
    for i in 1:length(xcoords)
        if xcoords[i] < 1 || xcoords[i] > size(θ,1)
            xcoords[i] = 1 # If not push them to the origin
        elseif ycoords[i] < 1 || ycoords[i] > size(θ,2)
            ycoords[i] = 1
        end
    end
    vertex₁ = [xcoords[1], ycoords[1]]
    vertex₂ = [xcoords[2], ycoords[2]]
    vertex₃ = [xcoords[3], ycoords[3]]
    P = VPolygon([vertex₁, vertex₂, vertex₃])
    for i in 1:size(θ,1), j in 1:size(θ,2)
        if [Float64(i), Float64(j)] ∈ P
            θ[j,i] = θ₀ + δₐ
        else
            θ[j,i] = θ₀
        end
    end
    return θ, P
end


text:
# The double call to using GDAL and other GIS packages here is necessary
# to push the high-level API functions provided by those packages into Main
# See DataFrames for precedent
using GDAL

module GIS

using Images, ImageView
using GDAL

include("display.jl")

export display

end


text:
const KeyType = Union{Symbol, Int}

const Arguments = Vector{Any}
const NamedArguments = Dictionary{Symbol, Any}

arguments(x) = collect(Any, x)
namedarguments(x) = NamedArguments(keys(x), values(x))

function separate(f, d::AbstractDictionary)
    d1, d2 = empty(d), empty(d)
    for (k, v) in pairs(d)
        target = ifelse(f(v), d1, d2)
        insert!(target, k, v)
    end
    return d1, d2
end

# `f` takes a pair and returns a unique value
function map_pairs(f, s)
    ks, vs = keys(s), values(s)
    res = collect(Any, Iterators.map(f∘Pair,  ks, vs))
    return eltype(ks) <: Symbol ? NamedArguments(collect(Symbol, ks), res) : res
end


text:
"""
I'm not sure on the details, but I think this function can only allocate up to

524288 doubles

that is it may not be able to allocate arrays with more than about half a million elements, or 4 megabytes.
"""
@generated function alloca(::Val{N}, ::Type{T} = Float64, ::Val{Align} = Val{64}()) where {N, T, Align}
    typ = llvmtype(T)
    ptyp = llvmtype(Int)
    decls = String[]
    instrs = String[]
    push!(instrs, "%ptr = alloca $typ, i32 $N, align $Align")
    push!(instrs, "%iptr = ptrtoint $typ* %ptr to $ptyp")
    # push!(instrs, "%ptr = alloca i8, i32 $(N*sizeof(T)), align $Align")
    # push!(instrs, "%iptr = ptrtoint i8* %ptr to $ptyp")
    push!(instrs, "ret $ptyp %iptr")
    quote
        $(Expr(:meta,:inline))
        Base.llvmcall(
            $((join(decls, "\n"), join(instrs, "\n"))),
            Ptr{$T}, Tuple{}
        )
    end
end

function valloc(::Type{T}, N::Int, sz::Int) where {T}
    @assert N > 0
    @assert sz >= 0
    # We use padding to align the address of the first element, and
    # also to ensure that we can access past the last element up to
    # the next full vector width
    padding = N-1 + mod(-sz, N)
    mem = Vector{T}(undef, sz + padding)
    addr = Int(pointer(mem))
    off = mod(-addr, N * sizeof(T))
    @assert mod(off, sizeof(T)) == 0
    off = fld(off, sizeof(T))
    @assert 0 <= off <= padding
    res = view(mem, off+1 : off+sz)
    addr2 = Int(pointer(res))
    @assert mod(addr2, N * sizeof(T)) == 0
    res
end
function valloc(f, ::Type{T}, N::Int, sz::Int) where {T}
    mem = valloc(T, N, sz)
    @inbounds for i in 1:sz
        mem[i] = f(i)
    end
    mem
end

# @generated function sub2ind(dims::NTuple{N}, I::NTuple{N}) where {N}
#     ex = :(I[$N] - 1)
#     for i = (N - 1):-1:1
#         ex = :(I[$i] - 1 + dims[$i] * $ex)
#     end
#     quote
#         $(Expr(:meta, :inline))
#         $ex + 1
#     end
# end

@generated function vload(
    ::Union{Type{Vec{N,T}},Val{N}}, ptr::Ptr{T},
    ::Val{Aligned}, ::Val{Nontemporal}
    # ::Val{Aligned} = Val{false}(), ::Val{Nontemporal} = Val{false}()
) where {N,T,Aligned, Nontemporal}
    @assert isa(Aligned, Bool)
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    decls = String[]
    instrs = String[]
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    flags = [""]
    align > 0 && push!(flags, "align $align")
    Nontemporal && push!(flags, "!nontemporal !{i32 1}")
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
    push!(instrs, "%res = load $vtyp, $vtyp* %ptr" * join(flags, ", "))
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$N,$T}, Tuple{Ptr{$T}}, ptr)
    end
end
@generated function vload(
    ::Union{Type{Vec{N,T}},Val{N}}, ptr::Ptr{T}, mask::Vec{N,Bool}, ::Val{Aligned}# = Val{false}()
) where {N,T,Aligned}
    @assert isa(Aligned, Bool)
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    btyp = llvmtype(Bool)
    vbtyp = "<$N x $btyp>"
    decls = String[]
    instrs = String[]
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end

    push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
    push!(instrs, "%mask = trunc $vbtyp %1 to <$N x i1>")
    push!(decls,
        "declare $vtyp @llvm.masked.load.$(suffix(N,T))($vtyp*, i32, <$N x i1>, $vtyp)"
    )
    push!(instrs,
        "%res = call $vtyp @llvm.masked.load.$(suffix(N,T))($vtyp* %ptr, i32 $align, <$N x i1> %mask, $vtyp $(llvmconst(N, T, 0)))"
    )
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$N,$T}, Tuple{Ptr{$T}, Vec{$N,Bool}}, ptr, mask)
    end
end
@generated function vload(
    ::Union{Type{Vec{N,T}},Val{N}}, ptr::Ptr{T}, mask::U, ::Val{Aligned}# = Val{false}()
) where {N,T,Aligned,U<:Unsigned}
    @assert isa(Aligned, Bool)
    @assert 8sizeof(U) >= N
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    mtyp_input = llvmtype(U)
    mtyp_trunc = "i$N"
    decls = String[]
    instrs = String[]
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
    if mtyp_input == mtyp_trunc
        push!(instrs, "%mask = bitcast $mtyp_input %1 to <$N x i1>")
    else
        push!(instrs, "%masktrunc = trunc $mtyp_input %1 to $mtyp_trunc")
        push!(instrs, "%mask = bitcast $mtyp_trunc %masktrunc to <$N x i1>")
    end
    push!(decls,
        "declare $vtyp @llvm.masked.load.$(suffix(N,T))($vtyp*, i32, <$N x i1>, $vtyp)"
    )
    push!(instrs,
        "%res = call $vtyp @llvm.masked.load.$(suffix(N,T))($vtyp* %ptr, i32 $align, <$N x i1> %mask, $vtyp $(llvmconst(N, T, 0)))"
    )
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$N,$T}, Tuple{Ptr{$T}, $U}, ptr, mask)
    end
end

for v ∈ (:Vec, :SVec, :Val)
    vargs = [v === :Val ? :(::Val{W}) : :(::Type{$v{W,T}})]
    for ptr ∈ (:Ptr, :Pointer, :ZeroInitializedPointer)
        pargs = push!(copy(vargs), :(ptr::$ptr{T}))
        for index ∈ (true,false)
            if index
                icall = Union{Symbol,Expr}[:(Vec{W,T}), ptr == :Ptr ? :(ptr + i) : :(ptr.ptr + sizeof(T)*i)]
                iargs = push!(copy(pargs), :(i::Int))
            else
                icall = Union{Symbol,Expr}[:(Vec{W,T}), ptr == :Ptr ? :ptr : :(ptr.ptr)]
                iargs = pargs
            end
            for mask ∈ (true,false)
                if mask
                    margs = push!(copy(iargs), :(mask::Union{Vec{W,Bool},<:Unsigned}))
                    mcall = push!(copy(icall), :mask)
                    fopts = (:vload,:vloada)
                else
                    margs = iargs
                    mcall = icall
                    fopts = (:vload,:vloada,:vloadnt)
                end
                for f ∈ fopts
                    fcall = copy(mcall)
                    if f === :vload # aligned arg
                        push!(fcall, :(Val{false}()))
                    else
                        push!(fcall, :(Val{true}()))
                    end
                    if !mask
                        if f === :vloadnt # nontemporal arg
                            push!(fcall, :(Val{true}()))
                        else
                            push!(fcall, :(Val{false}()))
                        end
                    end
                    if ptr === :ZeroInitializedPointer
                        body = Expr(:tuple, :vbroadcast, :(Vec{W,T}), :(zero(T)))
                    else
                        body = Expr(:call, :vload, fcall...)
                        if v !== :Vec
                            body = :(SVec($body))
                        end
                    end
                    @eval @inline function $f($(margs...)) where {W,T}
                        $body
                    end
                end
            end
        end
    end
end

@generated function vstore!(
    ptr::Ptr{T}, v::Vec{N,T},
    ::Val{Aligned}, ::Val{Nontemporal}
    # ::Val{Aligned} = Val{false}(), ::Val{Nontemporal} = Val{false}()
) where {N,T,Aligned, Nontemporal}
    @assert isa(Aligned, Bool)
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    decls = String[]
    instrs = String[]
    if Aligned# || Nontemporal
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    flags = [""]
    align > 0 && push!(flags, "align $align")
    Nontemporal && push!(flags, "!nontemporal !{i32 1}")
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
    push!(instrs, "store $vtyp %1, $vtyp* %ptr" * join(flags, ", "))
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall(
            $((join(decls, "\n"), join(instrs, "\n"))),
            Cvoid, Tuple{Ptr{$T}, Vec{$N,$T}}, ptr, v
        )
    end
end

@generated function vstore!(
    ptr::Ptr{T}, v::Vec{N,T}, mask::Vec{N,Bool}, ::Val{Aligned}# = Val{false}()
) where {N,T,Aligned}
    @assert isa(Aligned, Bool)
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    btyp = llvmtype(Bool)
    vbtyp = "<$N x $btyp>"
    decls = String[]
    instrs = String[]
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
    push!(instrs, "%mask = trunc $vbtyp %2 to <$N x i1>")
    push!(decls,
        "declare void @llvm.masked.store.$(suffix(N,T))($vtyp, $vtyp*, i32, <$N x i1>)")
    push!(instrs,
        "call void @llvm.masked.store.$(suffix(N,T))($vtyp %1, $vtyp* %ptr, i32 $align, <$N x i1> %mask)")
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall(
            $((join(decls, "\n"), join(instrs, "\n"))),
            Cvoid, Tuple{Ptr{$T}, Vec{$N,$T}, Vec{$N,Bool}},
            ptr, v, mask
        )
    end
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{N,T}, mask::U, ::Val{Aligned}# = Val{false}()
) where {N,T,Aligned,U<:Unsigned}
    @assert isa(Aligned, Bool)
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    mtyp_input = llvmtype(U)
    mtyp_trunc = "i$N"
    decls = String[]
    instrs = String[]
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
    if mtyp_input == mtyp_trunc
        push!(instrs, "%mask = bitcast $mtyp_input %2 to <$N x i1>")
    else
        push!(instrs, "%masktrunc = trunc $mtyp_input %2 to $mtyp_trunc")
        push!(instrs, "%mask = bitcast $mtyp_trunc %masktrunc to <$N x i1>")
    end
    push!(decls,
        "declare void @llvm.masked.store.$(suffix(N,T))($vtyp, $vtyp*, i32, <$N x i1>)"
    )
    push!(instrs,
        "call void @llvm.masked.store.$(suffix(N,T))($vtyp %1, $vtyp* %ptr, i32 $align, <$N x i1> %mask)"
    )
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Cvoid, Tuple{Ptr{$T}, Vec{$N,$T}, $U},
            ptr, v, mask)
    end
end

for ptr ∈ (:Ptr, :Pointer, :ZeroInitializedPointer)
    pargs = [:(ptr::$ptr{T}), :(v::AbstractSIMDVector{W,T})]
    for index ∈ (true,false)
        if index
            icall = Union{Symbol,Expr}[ptr == :Ptr ? :(ptr + i) : :(ptr.ptr + sizeof(T)*i), :(extract_data(v))]
            iargs = push!(copy(pargs), :(i::Int))
        else
            icall = Union{Symbol,Expr}[ptr == :Ptr ? :ptr : :(ptr.ptr), :(extract_data(v))]
            iargs = pargs
        end
        for mask ∈ (true,false)
            if mask
                margs = push!(copy(iargs), :(mask::Union{Vec{W,Bool},<:Unsigned}))
                mcall = push!(copy(icall), :mask)
                fopts = (:vstore!,:vstorea!)
            else
                margs = iargs
                mcall = icall
                fopts = (:vstore!,:vstorea!,:vstorent!)
            end
            for f ∈ fopts
                fcall = copy(mcall)
                if f == :vstore! # aligned arg
                    push!(fcall, :(Val{false}()))
                else
                    push!(fcall, :(Val{true}()))
                end
                if !mask
                    if f == :vstorent! # nontemporal arg
                        push!(fcall, :(Val{true}()))
                    else
                        push!(fcall, :(Val{false}()))
                    end
                end
                body = Expr(:call, :vstore!, fcall...)
                @eval @inline function $f($(margs...)) where {W,T}
                    $body
                end
            end
        end
    end
end

@generated function vloadscope(
    ::Type{Vec{N,T}}, ptr::Ptr{T}, ::Val{Scope},
    ::Val{Aligned}, ::Val{Nontemporal}
) where {N, T, Scope, Aligned, Nontemporal}
    @assert isa(Aligned, Bool)
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    domain,scope,list = Scope::NTuple{3,Int}
    # domain, scope, list = "\"domain\"", "\"scope\"", "\"list\""
    decls = String["!$domain = !{!$domain}","!$scope = !{!$scope, !$domain}", "!$list = !{!$scope}"]
    # decls = String[]
    instrs = String[]
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    flags = [""]
    align > 0 && push!(flags, "align $align")
    push!(flags, "!alias.scope !$list")
    Nontemporal && push!(flags, "!nontemporal !{i32 1}")
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
    push!(instrs, "%res = load $vtyp, $vtyp* %ptr" * join(flags, ", "))
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$N,$T}, Tuple{Ptr{$T}}, ptr)
    end
end
@generated function vstorescope!(
    ptr::Ptr{T}, v::Vec{N,T}, ::Val{Scope},
    ::Val{Aligned}, ::Val{Nontemporal}
    # ::Val{Aligned} = Val{false}(), ::Val{Nontemporal} = Val{false}()
) where {N,T,Scope,Aligned, Nontemporal}
    @assert isa(Aligned, Bool)
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    domain,scope,list = Scope::NTuple{3,Int}
    decls = String["!$domain = !{!$domain}","!$scope = !{!$scope, !$domain}", "!$list = !{!$scope}"]
    # decls = String[]
    instrs = String[]
    if Aligned# || Nontemporal
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    flags = [""]
    align > 0 && push!(flags, "align $align")
    push!(flags, "!noalias !$list")
    Nontemporal && push!(flags, "!nontemporal !{i32 1}")
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
    push!(instrs, "store $vtyp %1, $vtyp* %ptr" * join(flags, ", "))
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall(
            $((join(decls, "\n"), join(instrs, "\n"))),
            Cvoid, Tuple{Ptr{$T}, Vec{$N,$T}}, ptr, v
        )
    end
end
function nonaliased_store_and_load!(storeptr::Ptr{T}, v::Vec{W,T}, loadptr::Ptr{T}) where {W,T}
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    decls = """!0 = !{!0}
!1 = !{!1, !0}
!2 = !{!1}
"""
    instrs = """
%spt = inttoptr $ptyp %0 to $vtyp
%lpt = inttoptr $ptyp %2 to $vtyp
store $vtyp %1, $vtyp* %spt align $(sizeof(T)), !noalias !2
%res = load $vtyp, $vtyp* %lpt align $(sizeof(T)), !noalias !2
ret $vtyp %res
"""
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall(
            ($decls, $instrs),
            Vec{$W,$T},
            Tuple{Ptr{$T},Vec{$W,$T},Ptr{$T}},
            storeptr, v, loadptr
        )
    end
end

@generated function vloadconstant(
    ::Type{Vec{N,T}}, ptr::Ptr{T},
    ::Val{Aligned}, ::Val{Nontemporal}
) where {N, T, Aligned, Nontemporal}
    @assert isa(Aligned, Bool)
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    decls = String[]
    instrs = String[]
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    flags = [""]
    align > 0 && push!(flags, "align $align")
    push!(flags, "!tbaa !13")
    Nontemporal && push!(flags, "!nontemporal !{i32 1}")
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
    push!(instrs, "%res = load $vtyp, $vtyp* %ptr" * join(flags, ", "))
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$N,$T}, Tuple{Ptr{$T}}, ptr)
    end
end
# @generated function vloadscope(
#     ::Type{Vec{N,T}}, ptr::Ptr{T}, ::Val{Scope},
#     ::Val{Aligned}, ::Val{Nontemporal}
# ) where {N, T, Scope, Aligned, Nontemporal}
#     @assert isa(Aligned, Bool)
#     ptyp = llvmtype(Int)
#     typ = llvmtype(T)
#     vtyp = "<$N x $typ>"
#     domain,scope,list = Scope::NTuple{3,Int}
#     decls = String["!$domain = !{!$domain}","!$scope = !{!$scope, !$domain}", "!$list = !{!$scope}"]
#     # decls = String[]
#     instrs = String[]
#     if Aligned
#         align = N * sizeof(T)
#     else
#         align = sizeof(T)   # This is overly optimistic
#     end
#     flags = [""]
#     align > 0 && push!(flags, "align $align")
#     push!(flags, "!alias.scope !$list")
#     Nontemporal && push!(flags, "!nontemporal !{i32 1}")
#     push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
#     push!(instrs, "%res = load $vtyp, $vtyp* %ptr" * join(flags, ", "))
#     push!(instrs, "ret $vtyp %res")
#     quote
#         $(Expr(:meta, :inline))
#         Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
#             Vec{$N,$T}, Tuple{Ptr{$T}}, ptr)
#     end
# end
# @generated function vstorescope!(
#     ptr::Ptr{T}, v::Vec{N,T}, ::Val{Scope},
#     ::Val{Aligned}, ::Val{Nontemporal}
#     # ::Val{Aligned} = Val{false}(), ::Val{Nontemporal} = Val{false}()
# ) where {N,T,Scope,Aligned, Nontemporal}
#     @assert isa(Aligned, Bool)
#     ptyp = llvmtype(Int)
#     typ = llvmtype(T)
#     vtyp = "<$N x $typ>"
#     domain,scope,list = Scope::NTuple{3,Int}
#     decls = String["!$domain = !{!$domain}","!$scope = !{!$scope, !$domain}", "!$list = !{!$scope}"]
#     # decls = String[]
#     instrs = String[]
#     if Aligned# || Nontemporal
#         align = N * sizeof(T)
#     else
#         align = sizeof(T)   # This is overly optimistic
#     end
#     flags = [""]
#     align > 0 && push!(flags, "align $align")
#     push!(flags, "!noalias !$list")
#     Nontemporal && push!(flags, "!nontemporal !{i32 1}")
#     push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
#     push!(instrs, "store $vtyp %1, $vtyp* %ptr" * join(flags, ", "))
#     push!(instrs, "ret void")
#     quote
#         $(Expr(:meta, :inline))
#         Base.llvmcall(
#             $((join(decls, "\n"), join(instrs, "\n"))),
#             Cvoid, Tuple{Ptr{$T}, Vec{$N,$T}}, ptr, v
#         )
#     end
# end

@generated function gather(
   ptr::Vec{N,Ptr{T}}, ::Val{Aligned} = Val{false}()
) where {N,T,Aligned}
    @assert isa(Aligned, Bool)
    ptyp = llvmtype(Int)
    vptyp = "<$N x $ptyp>"
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    vptrtyp = "<$N x $typ*>"
    decls = String[]
    instrs = String[]
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    mask = join((", i1 true" for i ∈ 2:N))
    push!(instrs, "%ptr = inttoptr $vptyp %0 to $vptrtyp")
    push!(decls,
        "declare $vtyp @llvm.masked.gather.$(suffix(N,T)).$(suffix(N,Ptr{T}))($vptrtyp, i32, <$N x i1>, $vtyp)")
    push!(instrs,
        "%res = call $vtyp @llvm.masked.gather.$(suffix(N,T)).$(suffix(N,Ptr{T}))($vptrtyp %ptr, i32 $align, <$N x i1> <i1 true$(mask)>, $vtyp undef)")#$(llvmconst(N, T, 0)))")
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall(
            $((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$N,$T}, Tuple{Vec{$N,Ptr{$T}}}, ptr
        )
    end
end

@generated function gather(
   ptr::Vec{N,Ptr{T}}, mask::U, ::Val{Aligned} = Val{false}()
) where {N,T,Aligned,U<:Unsigned}
    @assert isa(Aligned, Bool)
    @assert 8sizeof(U) >= N
    ptyp = llvmtype(Int)
    vptyp = "<$N x $ptyp>"
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    vptrtyp = "<$N x $typ*>"
    mtyp_input = llvmtype(U)
    mtyp_trunc = "i$N"
    decls = String[]
    instrs = String[]
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    push!(instrs, "%ptr = inttoptr $vptyp %0 to $vptrtyp")
    if mtyp_input == mtyp_trunc
        push!(instrs, "%mask = bitcast $mtyp_input %1 to <$N x i1>")
    else
        push!(instrs, "%masktrunc = trunc $mtyp_input %1 to $mtyp_trunc")
        push!(instrs, "%mask = bitcast $mtyp_trunc %masktrunc to <$N x i1>")
    end
    push!(decls,
        "declare $vtyp @llvm.masked.gather.$(suffix(N,T)).$(suffix(N,Ptr{T}))($vptrtyp, i32, <$N x i1>, $vtyp)")
    push!(instrs,
        "%res = call $vtyp @llvm.masked.gather.$(suffix(N,T)).$(suffix(N,Ptr{T}))($vptrtyp %ptr, i32 $align, <$N x i1> %mask, $vtyp $(llvmconst(N, T, 0)))")
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall(
            $((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$N,$T}, Tuple{Vec{$N,Ptr{$T}}, $U}, ptr, mask
        )
    end
end
@generated function scatter!(
    ptr::Vec{N,Ptr{T}}, v::Vec{N,T}, ::Val{Aligned} = Val{false}()
) where {N,T,Aligned}
    @assert isa(Aligned, Bool)
    ptyp = llvmtype(Int)
    vptyp = "<$N x $ptyp>"
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    vptrtyp = "<$N x $typ*>"
    mtyp_input = llvmtype(U)
    mtyp_trunc = "i$N"
    decls = String[]
    instrs = String[]
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    push!(instrs, "%ptr = inttoptr $vptyp %1 to $vptrtyp")
    mask = join((", i1 true" for i ∈ 2:N))
    # push!(decls, "declare void @llvm.masked.scatter.$(suffix(N,T)).$(suffix(N,Ptr{T}))($vtyp, $vptrtyp, i32, <$N x i1>)")
    # push!(instrs, "call void @llvm.masked.scatter.$(suffix(N,T)).$(suffix(N,Ptr{T}))($vtyp %0, $vptrtyp %ptr, i32 $align, <$N x i1> <i1 true$(mask)>)")
    push!(decls, "declare void @llvm.masked.scatter.$(suffix(N,T))($vtyp, $vptrtyp, i32, <$N x i1>)")
    push!(instrs, "call void @llvm.masked.scatter.$(suffix(N,T))($vtyp %0, $vptrtyp %ptr, i32 $align, <$N x i1> <i1 true$(mask)>)")
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Cvoid, Tuple{Vec{$N,$T}, Vec{$N,Ptr{$T}}, $U}, v, ptr, mask)
    end
end
@generated function scatter!(
    ptr::Vec{N,Ptr{T}}, v::Vec{N,T}, mask::U, ::Val{Aligned} = Val{false}()
) where {N,T,Aligned,U<:Unsigned}
    @assert isa(Aligned, Bool)
    @assert 8sizeof(U) >= N
    ptyp = llvmtype(Int)
    vptyp = "<$N x $ptyp>"
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    vptrtyp = "<$N x $typ*>"
    mtyp_input = llvmtype(U)
    mtyp_trunc = "i$N"
    decls = String[]
    instrs = String[]
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    push!(instrs, "%ptr = inttoptr $vptyp %1 to $vptrtyp")
    if mtyp_input == mtyp_trunc
        push!(instrs, "%mask = bitcast $mtyp_input %2 to <$N x i1>")
    else
        push!(instrs, "%masktrunc = trunc $mtyp_input %2 to $mtyp_trunc")
        push!(instrs, "%mask = bitcast $mtyp_trunc %masktrunc to <$N x i1>")
    end
    # push!(decls,
        # "declare void @llvm.masked.scatter.$(suffix(N,T)).$(suffix(N,Ptr{T}))($vtyp, $vptrtyp, i32, <$N x i1>)")
    # push!(instrs,
        # "call void @llvm.masked.scatter.$(suffix(N,T)).$(suffix(N,Ptr{T}))($vtyp %0, $vptrtyp %ptr, i32 $align, <$N x i1> %mask)")
    push!(decls, "declare void @llvm.masked.scatter.$(suffix(N,T))($vtyp, $vptrtyp, i32, <$N x i1>)")
    push!(instrs, "call void @llvm.masked.scatter.$(suffix(N,T))($vtyp %0, $vptrtyp %ptr, i32 $align, <$N x i1> %mask)")
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall(
            $((join(decls, "\n"), join(instrs, "\n"))),
            Cvoid, Tuple{Vec{$N,$T}, Vec{$N,Ptr{$T}}, $U}, v, ptr, mask
        )
    end
end
@inline function vadd(ptr::Ptr{T}, inds::Vec{W,I}) where {W,T,I<:Union{UInt,Int}}
    pirate_reinterpret(Vec{W,Ptr{T}}, vadd(vbroadcast(Vec{W,I}, reinterpret(I, ptr)), inds))
end
@inline function vmuladd(s::I, inds::Vec{W,I}, ptr::Ptr{T}) where {W,T,I<:Union{UInt,Int}}
    pirate_reinterpret(Vec{W,Ptr{T}}, vmuladd(vbroadcast(Vec{W,I}, s), inds, vbroadcast(Vec{W,I}, reinterpret(I, ptr))))
end
@inline scatter!(ptr::Ptr{T}, inds::Vec{N,I}, v::Vec{N,T}, mask::U) where {N,T,I<:IntegerTypes,U<:Unsigned} = scatter!(vmuladd(sizeof(T), inds, ptr), v, mask)
@inline scatter!(ptr::Ptr{T}, inds::Vec{N,I}, v::Vec{N,T}) where {N,T,I<:IntegerTypes} = scatter!(vmuladd(sizeof(T), inds, ptr), v)
@inline gather(ptr::Ptr{T}, inds::Vec{N,I}, mask::U) where {N,T,I<:IntegerTypes,U<:Unsigned} = gather(vmuladd(sizeof(T),inds,ptr), mask)
@inline gather(ptr::Ptr{T}, inds::Vec{N,I}) where {N,T,I<:IntegerTypes} = gather(vmuladd(sizeof(T),inds,ptr))

@inline vload(::Type{Vec{W,T}}, A::AbstractArray, args...) where {W,T} = vload(Vec{W,T}, VectorizationBase.vectorizable(A), args...)
@inline gather(A::AbstractArray, args...) = gather(VectorizationBase.vectorizable(A), args...)
@inline vstore!(A::AbstractArray, args...) = vstore!(VectorizationBase.vectorizable(A), args...)
@inline scatter!(A::AbstractArray, args...) = scatter!(VectorizationBase.vectorizable(A), args...)

@generated function lifetime_start!(ptr::Ptr{T}, ::Val{L}) where {L,T}
    ptyp = llvmtype(Int)
    decls = "declare void @llvm.lifetime.start(i64, i8* nocapture)"
    instrs = [
        "%ptr = inttoptr $ptyp %0 to i8*",
        "call void @llvm.lifetime.start(i64 $(L*sizeof(T)), i8* %ptr)",
        "ret void"
    ]
    quote
        $(Expr(:meta,:inline))
        Base.llvmcall(
            $((decls, join(instrs, "\n"))),
            Cvoid, Tuple{Ptr{$T}}, ptr
        )
    end
end
@generated function lifetime_end!(ptr::Ptr{T}, ::Val{L}) where {L,T}
    ptyp = llvmtype(Int)
    decls = "declare void @llvm.lifetime.end(i64, i8* nocapture)"
    instrs = [
        "%ptr = inttoptr $ptyp %0 to i8*",
        "call void @llvm.lifetime.end(i64 $(L*sizeof(T)), i8* %ptr)",
        "ret void"
    ]
    quote
        $(Expr(:meta,:inline))
        Base.llvmcall(
            $((decls, join(instrs, "\n"))),
            Cvoid, Tuple{Ptr{$T}}, ptr
        )
    end
end

@inline function lifetime_start!(ptr::Pointer{T}) where {T}
    SIMDPirates.lifetime_start!(pointer(ptr), Val{1}())
end
@inline function lifetime_start!(ptr::Ptr{T}) where {T}
    SIMDPirates.lifetime_start!(ptr, Val{1}())
end

@inline function lifetime_end!(ptr::Pointer{T}) where {T}
    SIMDPirates.lifetime_end!(pointer(ptr), Val{1}())
end
@inline function lifetime_end!(ptr::Ptr{T}) where {T}
    SIMDPirates.lifetime_end!(ptr, Val{1}())
end
# Fallback is to do nothing
@inline lifetime_start!(::Any) = nothing
@inline lifetime_end!(::Any) = nothing

@generated function noalias!(ptr::Ptr{T}) where {T}
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    decls = "define noalias $typ* @noalias($typ *%a) noinline { ret $typ* %a }"
    instrs = [
        "%ptr = inttoptr $ptyp %0 to $typ*",
        "%naptr = call $typ* @noalias($typ* %ptr)",
        "%jptr = ptrtoint $typ* %naptr to $ptyp",
        "ret $ptyp %jptr"
    ]
    quote
        $(Expr(:meta,:inline))
        Base.llvmcall(
            $((decls, join(instrs, "\n"))),
            Ptr{$T}, Tuple{Ptr{$T}}, ptr
        )
    end
end

@generated function compressstore!(
    ptr::Ptr{T}, v::Vec{N,T}, mask::U
) where {N,T,U<:Unsigned}
    @assert 8sizeof(U) >= N
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    mtyp_input = llvmtype(U)
    mtyp_trunc = "i$N"
    decls = String[]
    instrs = String[]
    push!(instrs, "%ptr = inttoptr $ptyp %1 to $typ*")
    if mtyp_input == mtyp_trunc
        push!(instrs, "%mask = bitcast $mtyp_input %2 to <$N x i1>")
    else
        push!(instrs, "%masktrunc = trunc $mtyp_input %2 to $mtyp_trunc")
        push!(instrs, "%mask = bitcast $mtyp_trunc %masktrunc to <$N x i1>")
    end
    push!(decls, "declare void @llvm.masked.compressstore.$(suffix(N,T))($vtyp, $typ*, <$N x i1>)")
    push!(instrs, "call void @llvm.masked.compressstore.$(suffix(N,T))($vtyp %0, $typ* %ptr, <$N x i1> %mask)")
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall(
            $((join(decls, "\n"), join(instrs, "\n"))),
            Cvoid, Tuple{Vec{$N,$T}, Ptr{$T}, $U}, v, ptr, mask
        )
    end
end

@generated function expandload!(
    ::Type{Vec{N,T}}, ptr::Ptr{T}, mask::U
) where {N,T,U<:Unsigned}
    @assert 8sizeof(U) >= N
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    vptrtyp = "<$N x $typ*>"
    mtyp_input = llvmtype(U)
    mtyp_trunc = "i$N"
    decls = String[]
    instrs = String[]
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $typ*")
    if mtyp_input == mtyp_trunc
        push!(instrs, "%mask = bitcast $mtyp_input %1 to <$N x i1>")
    else
        push!(instrs, "%masktrunc = trunc $mtyp_input %1 to $mtyp_trunc")
        push!(instrs, "%mask = bitcast $mtyp_trunc %masktrunc to <$N x i1>")
    end
    push!(decls, "declare $vtyp @llvm.masked.expandload.$(suffix(N,T))($typ*, <$N x i1>, $vtyp)")
    push!(instrs, "%res = call $vtyp @llvm.masked.expandload.$(suffix(N,T))($typ* %ptr, <$N x i1> %mask, $vtyp $(llvmconst(N, T, 0)))")
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall(
            $((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$N,$T}, Tuple{Ptr{$T}, $U}, ptr, mask
        )
    end
end

### LLVM syntax errors on invariants
# struct Invariant{L,T}
#     ivp::Ptr{Cvoid}
#     ptr::Ptr{T}
# end
# @generated function invariant_start!(ptr::Ptr{T}, ::Val{L}) where {L,T}
#     ptyp = llvmtype(Int)
#     decls = "declare {}* @llvm.invariant.start.p0i8(i64, i8* nocapture)"
#     instrs = [
#         "%ptr = inttoptr $ptyp %0 to i8*",
#         "%ivt = call {}* @llvm.invariant.start.p0i8(i64 $(L*sizeof(T)), i8* %ptr)",
#         "%ivp = ptrtoint {}* %ivt to $ptyp",
#         "ret %ivp"
#     ]
#     quote
#         $(Expr(:meta,:inline))
#         ivp = Base.llvmcall(
#             $((decls, join(instrs, "\n"))),
#             Ptr{Cvoid}, Tuple{Ptr{$T}}, ptr
#         )
#         Invariant{$(L*sizeof(T)),T}(ivp, ptr)
#     end
# end
# @generated function invariant_end!(ivp::Invariant{L}) where {L}
#     ptyp = llvmtype(Int)
#     decls = "declare void @llvm.invariant.end.p0i8({}*, i64, i8* nocapture)"
#     instrs = [
#         "%ivp = inttoptr $ptyp %0 to {}*",
#         "%ptr = inttoptr $ptyp %1 to i8*",
#         "call void @llvm.lifetime.end.p0i8({}* %ivp, i64 $(L), i8* %ptr)",
#         "ret void"
#     ]
#     quote
#         $(Expr(:meta,:inline))
#         Base.llvmcall(
#             $((decls, join(instrs, "\n"))),
#             Cvoid, Tuple{Ptr{$T}}, ivp.ivp, ivp.ptr
#         )
#     end
# end

@generated function vrange(::Val{W}) where {W}
    quote
        $(Expr(:meta,:inline))
        $(Expr(:tuple, [Core.VecElement(w) for w ∈ 0:W-1]...))
    end
end

@inline function vload(::Type{Vec{W,T}}, ptr::AbstractStridedPointer{T}) where {W,T}
    gather(vmuladd(stride(ptr), vrange(Val{W}()), ptr.ptr), Val{false}())
end
@inline function vload(::Type{Vec{W,T}}, ptr::AbstractStridedPointer{T}, U::Unsigned) where {W,T}
    gather(vmuladd(stride(ptr), vrange(Val{W}()), ptr.ptr), U, Val{false}())
end
@inline function vstore!(ptr::AbstractStridedPointer{T}, v::Vec{W,T}) where {W,T}
    scatter!(vmuladd(stride(ptr), vrange(Val{W}()), ptr.ptr), v, Val{false}())
end
@inline function vstore!(v::Vec{W,T}, ptr::AbstractStridedPointer{T}, U::Unsigned) where {W,T}
    scatter!(vmuladd(stride(ptr), vrange(Val{W}()), ptr.ptr), v, U, Val{false}())
end


text:
using Documenter, ModiaMedia

makedocs(
  modules  = [ModiaMedia],
  format   = :html,
  sitename = "ModiaMedia",
  authors  = "Martin Otter (DLR-SR), Hilding Elmqivst (Mogram), Chris Laughman (NERL)",
  pages    = [
     "Home"   => "index.md",
     "Library" => [
        "lib/Types.md",
        "lib/Functions.md"
        ]
  ]
)


text:
# License for this file: MIT (expat)
# Copyright 2017-2019, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Solids(Modia3D/Solids/_module.jl)
#

"""
    material = ElasticContactPairMaterial(;coefficientOfRestitution=0.0, slidingFrictionCoefficient=0.0,
                                           rotationalResistanceCoefficient=0.0, vsmall=0.01, wsmall=0.01)

Generates an `ElasticContactPairMaterial` object by providing the material properties of
two objects that are in contact to each other.

# Keyword Arguments
- coefficientOfRestitution: Coefficient of restitution between two objects (=0: inelastic ... =1: fully elastic).
- slidingFrictionCoefficient: Kinetic/sliding friction force coefficent between two objects (`>= 0.0`).
- rotationalResistanceCoefficient: Rotational resistance torque coefficient between two objects (`>= 0.0`).
  Its effect is that the contact torque is computed to reduce the relative angular velocity between two objects.
  For a ball, `rotationalResistanceCoefficient` is the (standard) rolling resistance coefficient.
- `vsmall` in [m/s]: Used for regularization when computing the unit vector in direction of
                     the relative tangential velocity to avoid a division by zero.
- `wsmall` in [rad/s]: Used for regularization when computing the unit vector in direction
                       of the relative angular velocity to avoid a division by zero.

# Example
```julia
import Modia3D
mat = Modia3D.ElasticContactPairMaterial(coefficientOfRestitution=0.5)
```
"""
mutable struct ElasticContactPairMaterial <: Modia3D.AbstractContactPairMaterial
    coefficientOfRestitution::Float64
    slidingFrictionCoefficient::Float64
    rotationalResistanceCoefficient::Float64
    vsmall::Float64
    wsmall::Float64

    function ElasticContactPairMaterial(; coefficientOfRestitution=0.0,
                                          slidingFrictionCoefficient=0.0,
                                          rotationalResistanceCoefficient=0.0,
                                          vsmall=0.01,
                                          wsmall=0.01)
        @assert(coefficientOfRestitution >= 0.0)
        @assert(coefficientOfRestitution <= 1.0)
        @assert(slidingFrictionCoefficient >= 0.0)
        @assert(rotationalResistanceCoefficient >= 0.0)
        new(coefficientOfRestitution, slidingFrictionCoefficient, rotationalResistanceCoefficient,
            vsmall, wsmall)
    end
end

"""
    material = combineContactPairMaterials(material1, material2)

Return a new `ElasticContactPairMaterial` object from the mean values of the ElasticContactPairMaterial
objects `material1` and `material2`.
"""
function combineContactPairMaterials(mat1::ElasticContactPairMaterial, mat2::ElasticContactPairMaterial)::ElasticContactPairMaterial
    mat = ElasticContactPairMaterial()
    mat.coefficientOfRestitution        = (mat1.coefficientOfRestitution        + mat2.coefficientOfRestitution)/2
    mat.slidingFrictionCoefficient      = (mat1.slidingFrictionCoefficient      + mat2.slidingFrictionCoefficient)/2
    mat.rotationalResistanceCoefficient = (mat1.rotationalResistanceCoefficient + mat2.rotationalResistanceCoefficient)/2
    mat.vsmall = (mat1.vsmall + mat2.vsmall)/2
    mat.wsmall = (mat1.wsmall + mat2.wsmall)/2
    return mat
end

struct TwoNamesKey
    name1::String
    name2::String
    function TwoNamesKey(name1::String, name2::String)
        if name1 < name2
            return new(name1, name2)
        else
            return new(name2, name1)
        end
    end
end

Base.:isequal(keyA::TwoNamesKey, keyB::TwoNamesKey) = keyA.name1 == keyB.name1 && keyA.name2 == keyB.name2
Base.:isless(keyA::TwoNamesKey, keyB::TwoNamesKey) = keyA.name1 < keyB.name1 || ( keyA.name1 == keyB.name1 && keyA.name2 < keyB.name2 )

"""
    readContactPairMaterialFromJSON(fileName)

Read a JSON file consisting of a dictionary of ContactPairMaterial instances from `fileName` and
return a `Dict{String, ContactPairMaterial}` dictionary.
"""
function readContactPairMaterialFromJSON(fileName::AbstractString)
    dict1 = JSON.parsefile(fileName)
    palette = Dict{TwoNamesKey, Modia3D.AbstractContactPairMaterial}()
    for (key1,value1) in dict1
        # Split a key "name1,name2" in two parts
        i = findfirst(",", key1)
        if typeof(i) == Nothing || i.start < 2 || i.start >= length(key1)
            println("... Wrong key \"", key1,"\" not included in contactPairMaterialPalette.")
        else
            key2 = TwoNamesKey( key1[1:i.start-1], key1[i.start+1:end] )
            obj  = ElasticContactPairMaterial()
            for (key3,value3) in value1
                setfield!(obj, Symbol(key3), value3)
            end
            palette[key2] = obj
        end
    end
    return palette
end

"""
    contactPairMaterialPalette

Dictionary of contact pair material data, see [`Modia3D.ElasticContactPairMaterial`](@ref)
"""
contactPairMaterialPalette = readContactPairMaterialFromJSON( joinpath(Modia3D.path, "palettes", "contactPairMaterials.json") )

struct NoContactPairMaterial <: Modia3D.AbstractContactPairMaterial
end

"""
    material = getContactPairMaterial(name1, name2)

Return a contact pair material object of the materials with the
names `name1` and `name2` from dictionary `Modia3D.contactPairMaterialPalette`.
"""
function getContactPairMaterial(name1::AbstractString, name2::AbstractString)::Modia3D.AbstractContactPairMaterial
    value = get(contactPairMaterialPalette, TwoNamesKey(name1,name2), NoContactPairMaterial())
    if typeof(value) == NoContactPairMaterial
        # Combination name1,name2 is not present
        value1 = get(contactPairMaterialPalette, TwoNamesKey(name1,name1), NoContactPairMaterial())
        value2 = get(contactPairMaterialPalette, TwoNamesKey(name2,name2), NoContactPairMaterial())
        if typeof(value1) == NoContactPairMaterial || typeof(value2) == NoContactPairMaterial
            error("No contact pair material (\"$name1\",\"$name2\") or\n",
                  "(\"$name1\",\"$name1\") and (\"$name2\",\"$name2\")\n",
                  "in Modia3D.contactPairMaterialPalette dictionary.")
        end
        return combineContactPairMaterials(value1, value2)
    end
    return value
end


text:
##############################################################################
#
# Basic printing of TextAnalysis types
#
##############################################################################

show(io::IO, d::AbstractDocument) = print(io, "A $(typeof(d))")
show(io::IO, crps::Corpus) = print(io, "A Corpus")
show(io::IO, dtm::DocumentTermMatrix) = print(io, "A DocumentTermMatrix")

##############################################################################
#
# Pretty printing of TextAnalysis types
#
##############################################################################

function summary(d::AbstractDocument)
    o = ""
    o *= "A $(typeof(d))\n"
    o *= " * Language: $(language(d))\n"
    o *= " * Title: $(title(d))\n"
    o *= " * Author: $(author(d))\n"
    o *= " * Timestamp: $(timestamp(d))\n"
    if contains(Any[TokenDocument, NGramDocument], typeof(d))
        o *= " * Snippet: ***SAMPLE TEXT NOT AVAILABLE***"
    else
        sample_text = replace(text(d)[1:50], r"\s+", " ")
        o *= " * Snippet: $(sample_text)"
    end
    return o
end

function summary(crps::Corpus)
    n = length(crps.documents)
    n_s = sum(map(d -> typeof(d) == StringDocument, crps.documents))
    n_f = sum(map(d -> typeof(d) == FileDocument, crps.documents))
    n_t = sum(map(d -> typeof(d) == TokenDocument, crps.documents))
    n_ng = sum(map(d -> typeof(d) == NGramDocument, crps.documents))
    o = ""
    o *= "A Corpus with $n documents:\n"
    o *= " * $n_s StringDocument's\n"
    o *= " * $n_f FileDocument's\n"
    o *= " * $n_t TokenDocument's\n"
    o *= " * $n_ng NGramDocument's\n\n"
    o *= "Corpus's lexicon contains $(lexicon_size(crps)) tokens\n"
    o *= "Corpus's index contains $(index_size(crps)) tokens"
    return o
end

function summary(dtm::DocumentTermMatrix)
    n, p = size(dtm.dtm)
    @sprintf "A %dx%d DocumentTermMatrix" n p
end

##############################################################################
#
# In the REPL, show the summary by default
#
##############################################################################

repl_show(io::IO, d::AbstractDocument) = print(io, summary(d))
repl_show(io::IO, crps::Corpus) = print(io, summary(crps))
repl_show(io::IO, dtm::DocumentTermMatrix) = print(io, summary(dtm))


text:
using PicoQuant
using LaTeXStrings
import Plots
#=
This script produces graphs showing the output distribution of a rqc converges
to the Porter-Thomas distribution as the depth of the rqc increases.
=#

# Set parameters of the rqc & number of samples to take for each value of depth.
depth = 30; n = 4; m = 4; samples = 100

# Allocate memory for entropy and inverse participation ratios.
average_entropy = zeros(depth); average_sqrd_entropy = zeros(depth)
IPR = zeros(5, depth); IPR_sqrd = zeros(5, depth)

# For each sample at each circuit depth, create a RQC, convert it to a
# TensorNetworkCircuit, contract the network to get the output state and
# calculate the entropy and inverse participation ratios.
for (i, d) in enumerate(1:depth)
    for _ = 1:samples
        rqc = create_RQC(n, m, d)

        # An InteractiveBackend is used by default
        rqc_tn = convert_qiskit_circ_to_network(rqc)
        add_input!(rqc_tn, "0"^length(rqc_tn.input_qubits))

        output_state = full_wavefunction_contraction!(rqc_tn, "vector")
        output_distribution = abs.(output_state).^2

        for k = 1:5
            ipr = 1/sum(output_distribution.^(2*k))
            IPR[k, d] += ipr
            IPR_sqrd[k, d] += ipr^2
        end

        entropy = -sum(output_distribution.*log.(output_distribution))
        average_entropy[i] += entropy
        average_sqrd_entropy[i] += entropy^2
    end
end

average_entropy /= samples
average_sqrd_entropy /= samples
IPR /= samples; IPR_sqrd /= samples

# Get the error in the average entropy.
error = sqrt.(abs.(average_sqrd_entropy - average_entropy.^2))./sqrt(samples)

# Plot the average entropy of the output distribution as a function of depth.
plt = Plots.plot(1:depth, average_entropy, xlabel = "Depth", ylabel = "Entropy",
                 label = "Average entropy of output distribution",
                 yerror = error,
                 title = "Convergence to Porter-Thomas Entropy")

N = 2^(n*m)
porter_thomas_entropy = log(N) - 1 + 0.5772156649

plt = Plots.plot!(plt, 1:depth, porter_thomas_entropy*ones(depth),
                  label = "Porter-Thomas entropy")
display(plt)
Plots.savefig("entropy_convergence.png")

# Divide by Porter-Thomas IPR and get the IPR error
N = BigFloat(N)
for k = 1:5
    IPR[k, :] *= factorial(2*k)/N^(2*k - 1)
    IPR_sqrd[k, :] *= (factorial(2*k)/(N^(2*k - 1)))^2
end
IPR_error = sqrt.(abs.(IPR_sqrd - IPR.^2))./sqrt(samples)

# Plot the average inverse participation ratio of the Porter-Thomas distribution
# divided by that of the RQC output distribution as a function of depth.
plt2 = Plots.plot(xlabel = "Depth", ylabel = L"N^{-k+1}k!/IPR^k",
                  title = "Convergence to Porter-Thomas IPR",
                  include_mathjax = "cdn")
for k = 1:5
    global plt2
    plt2 = Plots.plot!(plt2, 3:depth, IPR[k, 3:end],
                       label = "k = $(2*k)", yaxis = :log, color=k)
    plt2 = Plots.plot!(plt2, 3:depth, IPR[k, 3:end]+IPR_error[k, 3:end],
                       label = :none, yaxis = :log,
                       seriestype = :scatter, markershape = :hline, color=k)
end
display(plt2)
Plots.savefig("IPR_convergence.png")


text:
include("timearray.jl")
include("dataframe.jl")
include("auth_token.jl")


text:
# Julia wrapper for header: /usr/local/include/libavutil/fifo.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_fifo_alloc,
    av_fifo_free,
    av_fifo_reset,
    av_fifo_size,
    av_fifo_space,
    av_fifo_generic_read,
    av_fifo_generic_write,
    av_fifo_realloc2,
    av_fifo_drain,
    av_fifo_peek2,
    av_fifo_peek

function av_fifo_alloc(size::Integer)
    ccall((:av_fifo_alloc,libavutil),Ptr{AVFifoBuffer},(Uint32,),size)
end

function av_fifo_free(f)
    ccall((:av_fifo_free,libavutil),Void,(Ptr{AVFifoBuffer},),f)
end

function av_fifo_reset(f)
    ccall((:av_fifo_reset,libavutil),Void,(Ptr{AVFifoBuffer},),f)
end

function av_fifo_size(f)
    ccall((:av_fifo_size,libavutil),Cint,(Ptr{AVFifoBuffer},),f)
end

function av_fifo_space(f)
    ccall((:av_fifo_space,libavutil),Cint,(Ptr{AVFifoBuffer},),f)
end

function av_fifo_generic_read(f,dest,buf_size::Integer,func)
    ccall((:av_fifo_generic_read,libavutil),Cint,(Ptr{AVFifoBuffer},Ptr{Void},Cint,Ptr{Void}),f,dest,buf_size,func)
end

function av_fifo_generic_write(f,src,size::Integer,func)
    ccall((:av_fifo_generic_write,libavutil),Cint,(Ptr{AVFifoBuffer},Ptr{Void},Cint,Ptr{Void}),f,src,size,func)
end

function av_fifo_realloc2(f,size::Integer)
    ccall((:av_fifo_realloc2,libavutil),Cint,(Ptr{AVFifoBuffer},Uint32),f,size)
end

function av_fifo_drain(f,size::Integer)
    ccall((:av_fifo_drain,libavutil),Void,(Ptr{AVFifoBuffer},Cint),f,size)
end

function av_fifo_peek2(f,offs::Integer)
    ccall((:av_fifo_peek2,libavutil),Ptr{Uint8},(Ptr{AVFifoBuffer},Cint),f,offs)
end

function av_fifo_peek(f,offs::Integer)
    ccall((:av_fifo_peek,libavutil),Uint8,(Ptr{AVFifoBuffer},Cint),f,offs)
end


text:
using Random
using ProgressMeter
using PyMNE # not in Unfold.jl currently
using MixedModelsPermutations
#using BlockDiagonals # for olsranef
using StatsModels  # for olsranef

function cluster_permutation_test(mres::UnfoldLinearMixedModel,
                    dat::Array,
                    times::StepRangeLen,
                    coeffOfInterest = 2;
                    clusterFormingThreshold = 2,
                    adjacency = nothing,
                    nPerm = 500,
                    test_times = (minimum(times),maximum(times))

                    )

# from where to where to actually test (can be subrange)
@assert(test_times[1]<test_times[2])
tRange = findfirst(times.>=test_times[1]):findfirst(times.>=test_times[2])

# extract observed z-values
z_obs = [m.z for m in coefpvalues(mres.modelinfo) if String(m.coefname)==coefnames(mres.X.formulas)[2][coeffOfInterest]]

# cluster observed
obs_cluster = pymne_cluster(z_obs,clusterFormingThreshold;tRange = tRange,adjacency=adjacency)

# if there is no cluster in the real data, return empty p_df
if isempty(obs_cluster[2])
    p_df = DataFrame(:from=>[],:to=>[],:pval=>[],:coefname=>[])
    return p_df
end

# this is the meat of the function.
permDat = cluster_permutation(mres,dat,tRange,coeffOfInterest,nPerm)

# cluster permuted
perm_cluster = [pymne_cluster(permDat[:,p],clusterFormingThreshold;adjacency=adjacency) for p = 1:nPerm]

# get maximal cluster value
perm_H0 = maximum.([abs.(p[2]) for p in perm_cluster if ~isempty(p[2])])
perm_H0 = vcat(perm_H0,zeros(nPerm-length(perm_H0)))

# get pvalues
p_vals = PyMNE.stats.cluster_level._pval_from_histogram(obs_cluster[2], perm_H0, 0)

fs = 1 ./Float64(times.step)
# get p-values, we have to shift by the test_times[1] though
if clusterFormingThreshold == "tfce"
    # not sure how to handle this nicer...
    fromList = [o.start ./ fs - test_times[1] for o in obs_cluster[1]]
    toList = [o.stop ./fs - test_times[1] for o in obs_cluster[1]]
else
    fromList = [o[1].start ./ fs - test_times[1] for o in obs_cluster[1]]
toList = [o[1].stop ./fs - test_times[1] for o in obs_cluster[1]]
end
p_df = DataFrame(:from=>fromList,:to=>toList,:pval=>p_vals)

p_df[!,:coefname] .= coefnames(mres.X.formulas)[2][coeffOfInterest] # yeah this might break :S
return p_df
#println(p_vals)
#h = hist(perm_H0,bins=100)
#vlines!(h.axis,abs.(obs_cluster[2]),color="red")

end

cluster_permutation(args...;kwargs...) = cluster_permutation(MersenneTwister(1),args...;kwargs...)

function cluster_permutation(rng::AbstractRNG,mres,dat,tRange,coeffOfInterest,nPerm)
    permDat = Matrix{Float64}(undef,length(tRange),nPerm)
    mm_outer = Unfold.LinearMixedModel_wrapper(mres.X.formulas,dat[1,1,:],mres.X.Xs)

    chIx = 1 # for now we only support 1 channel anyway
    #
    #p = Progress(length(tRange))
    #Threads.@threads for tIx =1:length(tRange)
    #@showprogress "Processing Timepoints"
    for tIx =1:length(tRange)
        println(string(tIx))
        # splice in the correct data for residual calculation
        mm = deepcopy(mm_outer)
        mm.y .= dat[chIx,tRange[tIx],:]

        # set the previous calculated model-fit
        updateL!(setθ!(mm,Vector(mres.modelinfo.θ[tRange[tIx]])))

        # get the coefficient
        H0 = coef(mm)
        # set the one of interest to 0
        H0[coeffOfInterest] = 0
        # run the permutation
        # important here is to set the same seed to keep flip all time-points the same
        perm = permutation(deepcopy(rng),nPerm,mm;β=H0,blup_method=olsranef,use_threads=false); # constant rng to keep autocorr & olsranef for singular models

        # extract the z-value
        perm_z = [m.z for m in perm.coefpvalues if String(m.coefname)==coefnames(mm)[coeffOfInterest]]

        # save it
        permDat[tIx,:] = perm_z
        #next!(p)
    end
    return permDat
    end

    # function to call pymne -> _find_cluster function
    function pymne_cluster(data,clusterFormingThreshold;tRange = 1:size(data,1), adjacency=nothing)

        return PyMNE.stats.cluster_level._find_clusters(data,clusterFormingThreshold,adjacency=adjacency,include=tRange)
    end

    pymne_cluster(data,clusterFormingThreshold::String;kwargs...) = pymne_cluster(data,Dict(:start=>0,:step=>0.2);kwargs...)

function olsranefjf(model::LinearMixedModel{T}) where {T}

    l = size(model.reterms)[1]

    mat = Array{Any}(undef, l);
    code = Array{Any}(undef, l);
    ### I get the contrasts
    for i in 1:l
        trm = model.reterms[i];
        dim = size(trm.z)[1];
        cd = StatsModels.ContrastsMatrix(EffectsCoding(), trm.levels).matrix;
        cd = kron(cd,I(dim));
        code[i] = cd;
        mat[i] = trm*cd;
    end
    mat
    X = hcat(mat...)
    X1 = hcat(ones(size(X)[1]), X)
    fixef_res = response(model) - model.X*model.β;
    flatblups = X1'X1 \ X1'fixef_res;
    flatblups = deleteat!(flatblups,1)

    code_all = BlockDiagonal([x for x in code])

    flatblups = code_all*flatblups

    blups = Vector{Matrix{T}}()

    offset = 1
    for trm in model.reterms
        chunksize = size(trm, 2)
        ngrps = length(trm.levels)
        npreds = length(trm.cnames)
        re = Matrix{T}(reshape(view(flatblups, offset:(offset+chunksize-1)),
                               npreds, ngrps))
        offset += chunksize
        push!(blups, re)
    end
    return blups,dummy_scalings(model.reterms)
end


text:
using BDF, Compat.Test

#test splitting BDF at single time point
origFilePath = joinpath(dirname(@__FILE__), "Newtest17-256.bdf")
tSplitFilePath1 = joinpath(dirname(@__FILE__), "Newtest17-256_1.bdf")
tSplitFilePath2 = joinpath(dirname(@__FILE__), "Newtest17-256_2.bdf")
tSplitFilePathComb = joinpath(dirname(@__FILE__), "Newtest17-256_comb.bdf")

bdfHeader = readBDFHeader(origFilePath)
dats, evtTab, trigs, statusChan = readBDF(origFilePath)

splitBDFAtTime(origFilePath, 30)
dats1, evtTab1, trigs1, statusChan1 = readBDF(tSplitFilePath1)
dats2, evtTab2, trigs2, statusChan2 = readBDF(tSplitFilePath2)

writeBDF(tSplitFilePathComb, [dats1 dats2], [trigs1; trigs2], [statusChan1; statusChan2], bdfHeader["sampRate"][1])
datsComb, evtTabComb, trigsComb, statusChanComb = readBDF(tSplitFilePathComb)

rm(tSplitFilePath1)
rm(tSplitFilePath2)
rm(tSplitFilePathComb)

@test isequal(dats, datsComb)
@test isequal(evtTab, evtTabComb)
@test isequal(trigs, trigsComb)
@test isequal(statusChan, statusChanComb)


text:
using LabVIEW0

try
    # change the path accordingly
    include(get_script_path(raw"/users/eben60/Desktop/Julia/LabVIEW0.jl/src/LabVIEW/jl-scripts/09_combined_example.jl"))
    LabVIEW0.setglobals(isOK=true)
catch excep
    LabVIEW0.setglobals(isOK=false, excpn=(;excep, stack_trace=backtrace()))
end

server_0mq4lv(fns)


text:
@fastmath function ltc_gj!(du,u,p,t, n_in, n_neurons, w_sens, w_syns, w_gj)
  v = u
  stim, C, G_L, E_L, G_sens, s_sens, h_sens, E_sens, G_syns, s_syns, h_syns, E_syns, G_gj = ltc_gj_unpack_p(n_in,n_neurons,p)

  #I_senss = reshape([G_sens[s,n] * NNlib.σ((stim[s] - h_sens[s,n]) * s_sens[s,n]) * (v[n] - E_sens[s,n]) for s in 1:size(G_sens,1), n in 1:size(G_sens,2)], size(G_sens,1), :)
  #I_synss = reshape([G_syns[s,n] * NNlib.σ((v[s] - h_syns[s,n]) * s_syns[s,n]) * (v[n] - E_syns[s,n]) for s in 1:size(G_syns,1), n in 1:size(G_syns,2)], size(G_syns,1), :)

  #I_senss = Array{eltype(stim)}(undef, n_in, n_neurons)
  #I_synss = Array{eltype(stim)}(undef, n_neurons, n_neurons)
  I_senss = zeros(eltype(stim), n_in, n_neurons)
  I_synss = zeros(eltype(stim), n_neurons, n_neurons)
  I_gjs = zeros(eltype(stim), n_neurons, n_neurons)
  @inbounds @simd for n in 1:n_neurons
    for s in 1:n_in
      I_senss[s,n] = w_sens[s,n] * G_sens[s,n] * mysigm((stim[s] - h_sens[s,n]) * s_sens[s,n]) * (v[n] - E_sens[s,n])
    end
    for s in 1:n_neurons
      I_synss[s,n] = w_syns[s,n] * G_syns[s,n] * mysigm((v[s] - h_syns[s,n]) * s_syns[s,n]) * (v[n] - E_syns[s,n])
      I_gjs[s,n] = w_gj[s,n] * G_gj[s,n] * (v[n] - v[s])
    end
  end
  I_sens = reshape(sum(I_senss, dims=1), :)
  I_syns = reshape(sum(I_synss, dims=1), :)
  I_gj = reshape(sum(I_gjs, dims=1), :)

  # I_sens = zeros(eltype(stim), n_neurons)
  # I_syns = zeros(eltype(stim), n_neurons)
  # @inbounds for n in 1:n_neurons
  #   for s in 1:n_in
  #     #w_sens[s,n] == 0 && continue
  #     I_sens[n] += w_sens[s,n] * G_sens[s,n] * mysigm((stim[s] - h_sens[s,n]) * s_sens[s,n]) * (v[n] - E_sens[s,n])
  #   end
  #   # I_sens[n] += sum(G_sens[:,n] .* mysigm.((stim .- h_sens[:,n]) .* s_sens[:,n]) .* (v[n] .- E_sens[:,n]))
  #   for s in 1:n_neurons
  #     #w_sens[s,n] == 0 && continue
  #     I_syns[n] += w_syns[s,n] * G_syns[s,n] * mysigm((v[s] - h_syns[s,n]) * s_syns[s,n]) * (v[n] - E_syns[s,n])
  #   end
  #   # I_syns[n] += sum(G_syns[:,n] .* mysigm.((v .- h_syns[:,n]) .* s_syns[:,n]) .* (v[n] .- E_syns[:,n]))
  # end

  #du .= (-((G_L .* (v .- E_L)) .+ I_sens .+ I_syns) ./ C) #.* dZ
  @.. du = (-((G_L * (v - E_L)) + I_sens + I_syns + I_gj) / C) #.* dZ

  nothing
end

function ltc_gj_init_u0p(n_in, n_neurons; T=Float32)
  v = rand_uniform(T, 0.01, 0.1, n_neurons) #.± 0.1f0

  stim   = fill(T(0), n_in)
  C      = rand_uniform(T,  1.0,     1.002,   n_neurons)
  G_L    = rand_uniform(T,  0.0001,  0.1,     n_neurons)
  E_L    = rand_uniform(T, -0.3,     0.3,     n_neurons) #.± 0.1f0
  G_sens = rand_uniform(T,  0.0001,  0.1,     n_in, n_neurons)
  s_sens = rand_uniform(T,  3,       8,       n_in, n_neurons)
  h_sens = rand_uniform(T,  0.3,     0.8,     n_in, n_neurons)
  E_sens = rand_uniform(T, -0.9,     0.9,     n_in, n_neurons)
  G_syns = rand_uniform(T,  0.0001,  0.1,     n_neurons, n_neurons)
  s_syns = rand_uniform(T,  3,       8,       n_neurons, n_neurons)
  h_syns = rand_uniform(T,  0.3,     0.8,     n_neurons, n_neurons)
  E_syns = rand_uniform(T, -0.9,     0.9,     n_neurons, n_neurons)

  G_gj   = rand_uniform(T,  0.0001,  0.1,     n_neurons, n_neurons)

  u0 = v
  p = vcat( stim,
    C, G_L, E_L,
    vec(G_sens), vec(s_sens), vec(h_sens), vec(E_sens),
    vec(G_syns), vec(s_syns), vec(h_syns), vec(E_syns),
    vec(G_gj)
  )
  return u0, p
end

function ltc_gj_bounds(n_in, n_neurons; T=Float32)
  lb = vcat(
    [1.0 for _ in 1:n_neurons],               # C
    [0 for _ in 1:n_neurons],                 # G_L
    [-1 for _ in 1:n_neurons],                # E_L
    [0 for _ in 1:n_in*n_neurons],            # G_sens
    [2 for _ in 1:n_in*n_neurons],            # s_sens
    [-0.3 for _ in 1:n_in*n_neurons],          # h_sens
    [-1 for _ in 1:n_in*n_neurons],           # E_sens
    [0 for _ in 1:n_neurons*n_neurons],       # G_syns
    [2 for _ in 1:n_neurons*n_neurons],       # s_syns
    [-0.3 for _ in 1:n_neurons*n_neurons],     # h_syns
    [-1 for _ in 1:n_neurons*n_neurons],      # E_syns
    [0 for _ in 1:n_neurons*n_neurons],       # G_gj

    [-0.4 for _ in 1:n_neurons])             # v

  ub = vcat(
    [2.0 for _ in 1:n_neurons],
    [0.4 for _ in 1:n_neurons],
    [1 for _ in 1:n_neurons],
    [1 for _ in 1:n_in*n_neurons],
    [9 for _ in 1:n_in*n_neurons],
    [0.9 for _ in 1:n_in*n_neurons],
    [1 for _ in 1:n_in*n_neurons],
    [1 for _ in 1:n_neurons*n_neurons],
    [9 for _ in 1:n_neurons*n_neurons],
    [0.9 for _ in 1:n_neurons*n_neurons],
    [1 for _ in 1:n_neurons*n_neurons],
    [1 for _ in 1:n_neurons*n_neurons],       # G_gj

    [0.8 for _ in 1:n_neurons]) # v)

  T.(lb), T.(ub)
end

function ltc_gj_unpack_p(n_in,n_neurons,p)
  stiml = n_in
  Cl = n_neurons
  G_Ll = n_neurons
  E_Ll = n_neurons
  sensl = n_in*n_neurons
  synsl = n_neurons*n_neurons
  G_sensl = sensl
  s_sensl = sensl
  h_sensl = sensl
  E_sensl = sensl
  G_synsl = synsl
  s_synsl = synsl
  h_synsl = synsl
  E_synsl = synsl

  G_gjl = synsl

  s = 1
  stim = @view p[s:s+stiml-1]
  s += stiml
  C = @view p[s:s+Cl-1]
  s += Cl
  G_L = @view p[s:s+G_Ll-1]
  s += G_Ll
  E_L = @view p[s:s+E_Ll-1]
  s += E_Ll
  G_sens = @view p[reshape(s:s+G_sensl-1, n_in, n_neurons)]
  s += G_sensl
  s_sens = @view p[reshape(s:s+s_sensl-1, n_in, n_neurons)]
  s += s_sensl
  h_sens = @view p[reshape(s:s+h_sensl-1, n_in, n_neurons)]
  s += h_sensl
  E_sens = @view p[reshape(s:s+E_sensl-1, n_in, n_neurons)]
  s += E_sensl

  G_syns = @view p[reshape(s:s+G_synsl-1, n_neurons, n_neurons)]
  s += G_synsl
  s_syns = @view p[reshape(s:s+s_synsl-1, n_neurons, n_neurons)]
  s += s_synsl
  h_syns = @view p[reshape(s:s+h_synsl-1, n_neurons, n_neurons)]
  s += h_synsl
  E_syns = @view p[reshape(s:s+E_synsl-1, n_neurons, n_neurons)]
  s += E_synsl

  G_gj = @view p[reshape(s:s+G_gjl-1, n_neurons, n_neurons)]
  s += G_gjl

  return stim, C, G_L, E_L, G_sens, s_sens, h_sens, E_sens, G_syns, s_syns, h_syns, E_syns, G_gj
end

function LTCGJ(wiring, solver, sensealg;
  T=Float32, tspan=T.((0, 1)),# n_sens=n_neurons, n_out=n_neurons,
  #w_sens=ones(T, n_in, n_neurons),
  #w_syns=ones(T, n_neurons, n_neurons),
  #w_gj   = ones(T, n_neurons, n_neurons),
  mtkize=false, gen_jac=false, kwargs...)

  s_in = wiring.s_in
  n_neurons = wiring.n_total

  w_gj = ones(T, n_neurons, n_neurons)
  new_matrices = Dict(
    :w_gj => w_gj
  )
  w_sens = wiring.matrices[:w_sens]
  w_syns = wiring.matrices[:w_syns]
  wiring.matrices = merge(wiring.matrices, new_matrices)

  lb, ub = ltc_gj_bounds(s_in, n_neurons; T)
  u0, p = ltc_gj_init_u0p(s_in, n_neurons; T)
  dudt!(du,u,p,t) = ltc_gj!(du,u,p,t, s_in, n_neurons, w_sens, w_syns, w_gj)

  rnncell = BCTRNNCell(wiring, solver, sensealg, dudt!, u0, tspan, p, lb, ub; mtkize, gen_jac, kwargs...)
  MyRecur(rnncell)
end


text:
module UNSflow

using Dierckx
export Spline1D, derivative, evaluate

using ForwardDiff
export derivative

using PyPlot
export plot, scatter, figure, xlabel, ylabel, xlim, ylim, xticks, yticks, subplot, axes, legend, markers, axis

using NLsolve
export nlsolve, not_in_place

export camber_calc, update_boundpos, update_kinem, update_indbound,
update_downwash, update_a0anda1, place_tev, place_lev, update_a2toan,
mutual_ind, trapz, update_a2a3adot, update_bv, ind_vel, view_vorts,
wakeroll, lautat, lautat_wakeroll, ldvm, calc_forces, design_solve,
lesp_design_max, transfer_cm, theodorsen, anim_flow,
drone_trajectory_problem, simple_LLT, find_tstep, ldvm_E, update_adot,
ldvm_more, ldvm_E_more, calc_forces_more, calc_forces_E_more, interp,
QScorrect_lautat, QScorrect_ldvm, QSLLT_lautat, QSLLT_ldvm,
lautat_wakeroll_more, LLT_ldvm, update_externalvel

export KinemPar, KinemParwFlap, KinemPar2DOF, KinemPar2DFree,
MotionDef, KinemDef, KinemDefwFlap, EldUpDef, EldUptstartDef,
ConstDef, TwoDOFPar, TwoDFreePar, TwoDSurf, TwoDSurfwFlap,
TwoDSurf_2DOF, TwoDFreeSurf, TwoDVort, TwoDFlowField, KelvinCondition,
KelvinKutta, EldRampReturnDef, EldUpIntDef, EldUpInttstartDef, SinDef,
CosDef, TheoDef, TheoDefwFlap, TwoDFlowData, DelVortDef, patch,
ThreeDSurfSimple, KinemDef3D, ThreeDVort, ThreeDFieldSimple,
KelvinConditionLLTldvm, KelvinKuttaLLTldvm

include("types.jl")
include("calculations.jl")
include("solvers.jl")
include("postprocess.jl")

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -25.420040452488465 <= q <= 21.65452636072683)

                     add_NL_constraint(m, :(log(1 + exp(-0.9974376164994152 + 0.7314791518637902*log(1 + exp(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4]))) + 0.4470938857218645*log(1 + exp(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4]))) + 0.10010754470216687*log(1 + exp(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4]))) + -0.8594790184444285*log(1 + exp(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4]))) + -0.4658863904521331*log(1 + exp(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))))) + log(1 + exp(0.15396139316976853 + -0.4228510094486002*log(1 + exp(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4]))) + -0.8035341866400141*log(1 + exp(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4]))) + -0.07617073018159282*log(1 + exp(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4]))) + 0.3940887684836425*log(1 + exp(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4]))) + 0.6667343490308526*log(1 + exp(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))))) + log(1 + exp(0.9665409577623398 + -0.449175249339572*log(1 + exp(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4]))) + -0.2640107596849468*log(1 + exp(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4]))) + 0.7284419447850534*log(1 + exp(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4]))) + 0.5662248826250997*log(1 + exp(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4]))) + -0.21454884110410477*log(1 + exp(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))))) + log(1 + exp(-0.1582635851091898 + 0.5103556986884796*log(1 + exp(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4]))) + -0.6017569427714062*log(1 + exp(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4]))) + -0.23749573380464106*log(1 + exp(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4]))) + -0.2344836196764124*log(1 + exp(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4]))) + 0.07066758651932892*log(1 + exp(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))))) + log(1 + exp(-0.47638200661148833 + 0.9126631115783219*log(1 + exp(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4]))) + 0.35414400688245484*log(1 + exp(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4]))) + -0.36037763135774403*log(1 + exp(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4]))) + -0.18657915153735738*log(1 + exp(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4]))) + 0.7872605277677782*log(1 + exp(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"Market state in luxor strategy"
function luxormktstate(mafast::Float64, maslow::Float64)
  if isnan(mafast) || isnan(maslow)
    return :undefined
  end
  if mafast > maslow
    return :trendup
  end
  if mafast < maslow
    return :trenddown
  end
  # this would be rare
  return :undefined
end
# function luxormktstate(mafast::Float64, maslow::Float64)
#   # @match is overkill here but keeping it as a template
#   @match (mafast >= maslow) begin
#     true => :trendup
#     false => :trenddown
#     # reminds to always provide an unmatched option (no Nothing)
#     _ => :undefined
#   end
# end

"""
Target position and stop, limit prices (if any) for luxor strategy.
...
Returns `(poschg::Int64, Vector[limitprice, stopprice]`.
"""
function luxorposlogic(mktstate::Symbol,
                       mktchgh::Float64, mktchgl::Float64, pthresh::Float64,
                       targetqty::Int64,
                       position_actual_mut::Vector{Int64})
  posact = position_actual_mut[1]

  @match (posact == 0, mktstate) begin
    # enter long position
    (true, :trendup) => (targetqty, [mktchgh + pthresh, mktchgh])
    # exit from long position
    (false, :trenddown), if posact > 0 end => (-posact, Array(Float64, 0))
    # enter short position
    (true, :trenddown) => (-targetqty, [mktchgl - pthresh, mktchgl])
    # exit from short position
    (false, :trendup), if posact < 0 end => (-posact, Array(Float64, 0))

    # partial order fill: continue advancing position (buy/sell more)
    (false, :trendup), if 0 < posact < targetqty end => (targetqty-posact, [mktchgh + pthresh, mktchgh])
    (false, :trenddown), if 0 > posact > -targetqty end => (-abs(targetqty-posact), [mktchgl - pthresh, mktchgl])

    # position over target: move back to targetqty
    (false, :trendup), if posact > targetqty end => (-(posact-targetqty), Array(Float64, 0))
    (false, :trenddown), if posact < -targetqty end => (abs(posact-targetqty), Array(Float64, 0))

    # no match among action-triggering contidions above
    othercase => (0, Array(Float64, 0))
  end
end

"Target signal for luxor strategy."
function luxortarget(s_ohlc::Signal{OHLC},
                     ohlc_inds::Dict{Symbol,Int64},
                     position_actual_mut::Vector{Int64},
                     nsma_fast::Int64, nsma_slow::Int64,
                     pthreshold::Float64, targetqty::Int64)
  posact = position_actual_mut[1]
  targetqty > 0 || error("target quantity must be positive")

  # signals updating at each timestep (aka OHLC bar)
  buffpclose!(buff::Vector{Float64}, tohlc::OHLC) =
    sighistbuffer!(buff, tohlc[2][ohlc_inds[:close]])
  pcloseinit = s_ohlc.value[2][ohlc_inds[:close]]
  s_sma_fast = map(mean,
                   foldp(buffpclose!, initbuff(nsma_fast, pcloseinit), s_ohlc),
                   typ=Float64)
  s_sma_slow = map(mean,
                   foldp(buffpclose!, initbuff(nsma_slow, pcloseinit), s_ohlc),
                   typ=Float64)
  s_high = map(s -> s[2][ohlc_inds[:high]], s_ohlc, typ=Float64)
  s_low = map(s -> s[2][ohlc_inds[:low]], s_ohlc, typ=Float64)

  # market state signal
  s_mktstate = map(luxormktstate, s_sma_fast, s_sma_slow, typ=Symbol)

  # market state change Bool-signal
  s_mktchg = schange(s_mktstate)

  # signals updating only when market state changes
  s_mktchg_high = filterwhen(s_mktchg, s_high.value, s_high)
  s_mktchg_low = filterwhen(s_mktchg, s_low.value, s_low)

  # target position updates only when market state input changes
  s_target = map((mks, h, l) -> luxorposlogic(mks, h, l,
                                              pthreshold, targetqty,
                                              position_actual_mut),
                 s_mktstate, s_mktchg_high, s_mktchg_low)
  # targeting tuple to pass to the order processing function
  return s_target
end


text:
using ConstrainedDynamics

# Parameters
ex = [1.;0.;0.]
ey = [0.;1.;0.]
ez = [0.;0.;1.]

l1 = 2.0
x, y = .1, .1
r = l1 / 4
h = r / 10
box = Box(x, l1, y, l1, color = RGBA(1., 1., 0.))
cyl = Cylinder(r, h, 2 * r, color = RGBA(1., 0., 0.))
box2 = Box(h, 2r, 2r, 2r, color = RGBA(1., 0., 1.))

p0 = [0;l1 / 2;0]
p1 = -p0
p2 = [0;r;0]
p3 = -p2

# Initial orientation
q1 = Quaternion(RotY(pi / 2))

# Links
origin = Origin{Float64}()
link1 = Body(box)
link2 = Body(cyl)
link3 = Body(box2)

# Constraints
joint0to23 = EqualityConstraint(Revolute(origin, link2, p3 + p1, zeros(3), ez), Revolute(origin, link3, p2, zeros(3), ex))
joint1to23 = EqualityConstraint(Spherical(link1, link2, p1, p3), Spherical(link1, link3, p0, p3))

links = [link1;link2;link3]
constraints = [joint0to23;joint1to23]
shapes = [box,cyl,box2]

mech = Mechanism(origin, links, constraints, g = 0., tend = 20., shapes = shapes)
setPosition!(mech,link1,x = -p0)
setPosition!(mech,link2,x = p3 + p1)
setPosition!(mech,link3,x = -p2)
setVelocity!(mech,link2,ω = [0;0;1.])
setVelocity!(mech,link3,ω = [1.;0;0])
setForce!(mech,link2,τ = [0;0;1.])

simulate!(mech,save = true)
visualize!(mech)


text:
function linespmax(nodes::Nodes, lines::Lines)
    lines_pmax = JuMP.Containers.DenseAxisArray(zeros(length(nodes.id), length(nodes.id)), nodes.id, nodes.id)

    for id in lines.id
        j = lines.from[id]
        k = lines.to[id]

        lines_pmax[j, k] = lines.pmax[id]
        lines_pmax[k, j] = lines.pmax[id]
    end

    return lines_pmax
end

function bprime(nodes::Nodes, lines::Lines)
    # Create the B_prime matrix in one step
    B_prime = JuMP.Containers.DenseAxisArray(zeros(length(nodes.id), length(nodes.id)), nodes.id, nodes.id)

    for id in lines.id
        j = lines.from[id]
        k = lines.to[id]

        X = lines.reactance[id]				# reactance

        B_prime[j, k] = 1/X					# Fill YBUS with susceptance
        B_prime[k, j] = 1/X					# Symmetric matrix
    end

    for j in nodes.id
        B_prime[j, j] = -sum(B_prime[j, :])
    end

    return B_prime
end

max_issues_repo_name: bobbyxiong/redispatch-ptg

max_forks_repo_name: bobbyxiong/redispatch-ptg


text:
#=
Helper function for creating and manipulating qubit paulis.

Some code from https://github.com/blakejohnson/Qlab.jl

=#

export ## Methods
		pauli_mats,
		pauli_strs,
		pauli_decomp

const pX = complex128([[0.0 1.0;1.0 0.0]])
const pY = complex128([[0.0 -1.0im;1.0im 0.0]])
const pZ = complex128([[1.0 0.0;0.0 -1.0]])
const pI = eye(Complex128, 2)

function base_paulis(dim::Int)
	#Pauli projectors onto the bottom 2 levels of a dim-level system
	map({pI, pX, pY, pZ}) do p
		fullPauli = zeros(Complex128, dim, dim)
		fullPauli[1:2, 1:2] = p
		fullPauli
	end
end

function pauli_mats(n::Int, dim=2)
	#Multi-qubit Pauli operators for n qubits
	@assert n > 0 "You need at least 1 qubit!"
	n == 1 && return base_paulis(dim)
	map(x->kron(x...), product(base_paulis(dim), pauli_mats(n-1, dim)))
end

function pauli_mats(dims::Vector{Int})
	#Multi-qubit pauli operators for n multi-level systems
	length(dims) == 1 && return base_paulis(dims[1])
	map(x->kron(x...), product(base_paulis(dims[1]), pauli_mats(dims[2:end])))
end

function pauli_strs(n::Int; sortByWeight=false)
	#Return the strings of pauli operators, e.g.: II, IX, IY, IZ, ....
	if n == 0
		return [""]
	end
	ops = ["I" "X" "Y" "Z"]
	paulis = vec(ops .* pauli_strs(n-1))

	if (sortByWeight)
		sort!(paulis, by=hamming)
	end

	paulis
end

function hamming(pauliStr)
	# compute the Hamming weight of a Pauli operator (number of non-identity components)
	weight = 0
	for pauli in pauliStr
		weight += (pauli == 'I') ? 0 : 1
	end
	weight
end

function pauli_decomp(rho::Array{Complex128,2}, dims::Array{Int, 1}, cutoff=1e-3)
	#Provide a decomposition of a density matrix into qubit pauli operators
	dimScale = 2^length(dims)
	for (str, mat) in zip(pauli_strs(length(dims)), pauli_mats(dims))
		val = real(trace(rho*mat)) / dimScale
		if val > cutoff
			println(str, ": ", val)
		end
	end

end


text:
import Slurm: api

# show_flag may be ored
# SHOW_ALL
# SHOW_DETAIL
# SHOW_DETAIL2
# SHOW_MIXED

type JobTable
  jobs :: Ptr{api.job_info_msg_t}

  function JobTable(jobs :: Ptr{api.job_info_msg_t})
    jtable = new(r_jobs[])
    finalizer(jtable, delete)
    return jtable
  end
end

function delete(jt :: JobTable)
  if jt.jobs != C_NULL
    api.slurm_free_job_info_msg(jt.jobs)
    jt.jobs = C_NULL
  end
end

function JobTable(show_flag = api.SHOW_ALL)
  r_jobs = Ref{Ptr{api.job_info_msg_t}}()
  timestamp = round(Int, Dates.datetime2unix(now()))
  err = api.slurm_load_jobs(timestamp, r_jobs, show_flag)
  # TODO err

  JobTable(r_jobs[])
end

# slurm_load_job
# slurm_load_job_user
# slurm_load_job_yser

length(jt :: JobTable) = unsafe_load(jt.jobs).record_count
function getindex(jt :: JobTable, i)
  jobs = unsafe_load(jt.jobs).job_array
  unsafe_load(jobs, i)
end

typealias JobInfo api.slurm_job_info_t

# slurm_get_rem_time
# slurm_get_end_time

# slurm_job_cpus_allocated_on_node
# slurm_job_cpus_allocated_on_node_id

# slurm_job_cpus_allocated_str_on_node
# slurm_job_cpus_allocated_str_on_node_id

# slurm_notify_job

# slurm_pid2jobid

# slurm_complete_job
# slurm_kill_job


text:
using Test
using Sockets
using Base.Threads

include("../../src/JuliaBolt.jl")

using .JuliaBolt

@testset "Connection Tests" begin

    @testset "test_conn_timedout" begin
        address = (ip"127.0.0.1", UInt16(7687))
        connection = Connection(1, address, FakeSocket(address), max_connection_lifetime=0)
        @test timedout(connection)
    end

    @testset "test_conn_not_timedout_if_not_enabled" begin
        address = (ip"127.0.0.1", UInt16(7687))
        connection = Connection(1, address, FakeSocket(address), max_connection_lifetime=-1)
        @test !timedout(connection)
    end

    @testset "test_conn_not_timedout" begin
        address = (ip"127.0.0.1", UInt16(7687))
        connection = Connection(1, address, FakeSocket(address), max_connection_lifetime=999999999)
        @test !timedout(connection)
    end

end

mutable struct TestRig
    pool::Union{ConnectionPool, Nothing}
    TestRig() = new(nothing)
end

@testset "ConnectionPool Tests" begin
    function connector(address; kwargs...)
        return QuickConnection(FakeSocket(address))
    end

    function setUp(t::TestRig)
        t.pool = ConnectionPool(connector, (ip"127.0.0.1", UInt16(7687)))
    end

    function tearDown(t::TestRig)
        JuliaBolt.close(t.pool)
    end

    testRig = TestRig()

    function assert_pool_size(address, expected_active, expected_inactive, _pool=nothing)
        if _pool == nothing
            _pool = testRig.pool
        end

        connections = nothing
        try
            connections = _pool.connections[address]
        catch e
            @test expected_active == 0
            @test expected_inactive == 0
            return
        end

        @test expected_active == length([cx for cx in connections if cx.in_use])
        @test expected_inactive == length([cx for cx in connections if !cx.in_use])
    end

    @testset "test_can_acquire" begin
        setUp(testRig)
        address = (ip"127.0.0.1", UInt16(7687))
        connection = acquire_direct(testRig.pool, address)
        @test connection.address == address
        assert_pool_size(address, 1, 0)
        tearDown(testRig)
    end

    @testset "test_can_acquire_twice" begin
        setUp(testRig)
        address = (ip"127.0.0.1", UInt16(7687))
        connection_1 = acquire_direct(testRig.pool, address)
        connection_2 = acquire_direct(testRig.pool, address)
        @test connection_1.address == address
        @test connection_2.address == address
        @test connection_1 != connection_2
        assert_pool_size(address, 2, 0)
        tearDown(testRig)
    end

    @testset "test_can_acquire_two_addresses" begin
        setUp(testRig)
        address_1 = (ip"127.0.0.1", UInt16(7687))
        address_2 = (ip"127.0.0.1", UInt16(7474))
        connection_1 = acquire_direct(testRig.pool, address_1)
        connection_2 = acquire_direct(testRig.pool, address_2)
        @test connection_1.address == address_1
        @test connection_2.address == address_2
        assert_pool_size(address_1, 1, 0)
        assert_pool_size(address_2, 1, 0)
        tearDown(testRig)
    end

    @testset "test_can_acquire_and_release" begin
        setUp(testRig)
        address = (ip"127.0.0.1", UInt16(7687))
        connection = acquire_direct(testRig.pool, address)
        assert_pool_size(address, 1, 0)
        release(testRig.pool, connection)
        assert_pool_size(address, 0, 1)
    end

    @testset "test_releasing_twice" begin
        setUp(testRig)
        address = (ip"127.0.0.1", UInt16(7687))
        connection = acquire_direct(testRig.pool, address)
        release(testRig.pool, connection)
        assert_pool_size(address, 0, 1)
        release(testRig.pool, connection)
        assert_pool_size(address, 0, 1)
    end

    @testset "test_cannot_acquire_after_close" begin
        setUp(testRig)
        address = (ip"127.0.0.1", UInt16(7687))
        pool = ConnectionPool((a)->QuickConnection(FakeSocket(a)), address)
        JuliaBolt.close(pool)
        @test_throws ErrorException acquire_direct(pool, address)
    end

    @testset "test_in_use_count" begin
        setUp(testRig)
        address = (ip"127.0.0.1", UInt16(7687))
        @test in_use_connection_count(testRig.pool, address) == 0

        connection = acquire_direct(testRig.pool, address)
        @test in_use_connection_count(testRig.pool, address) == 1

        release(testRig.pool, connection)
        @test in_use_connection_count(testRig.pool, address) == 0
    end

    @testset "test_max_conn_pool_size" begin
        address = (ip"127.0.0.1", UInt16(7687))
        pool = ConnectionPool(connector, address, max_connection_pool_size=1, connection_acquisition_timeout=0)
        connection = acquire_direct(pool, address)
        @test in_use_connection_count(pool, address) == 1
        @test_throws ErrorException acquire_direct(pool, address)
        @test in_use_connection_count(pool, address) == 1
    end

    # @testset "test_multithread" begin
    #     address = (ip"127.0.0.1", UInt16(7687))
    #     pool = ConnectionPool(connector, address, max_connection_pool_size=5, connection_acquisition_timeout=10)

    #     connections = Vector(undef, 10)

    #     @threads for i in 1:10
    #         connections[i] = acquire_direct(pool, address)
    #     end

    #     assert_pool_size(address, 5, 0, pool)

    #     @threads for i in 1:10
    #         release(pool, connections[i])
    #     end

    #     assert_pool_size(address, 0, 5, pool)
    # end
end

max_stars_repo_name: virtualgraham/JuliaBolt

max_issues_repo_name: virtualgraham/JuliaBolt

max_forks_repo_name: virtualgraham/JuliaBolt


text:
# # unfished
# add_theme(:ggplot2_base,
#     bglegend = _invisible,
#     fg       = :white,
#     fglegend = _invisible,
#     fgguide  = :black)
#
# add_theme(:ggplot2,
#           base = :ggplot2_base,
#           bginside = :lightgray,
#           fg       = :lightgray,
#           fgtext   = :gray,
#           fglegend = :gray,
#           fgguide  = :black)
#
# add_theme(:ggplot2_grey, base = :ggplot2)
#
# add_theme(:ggplot2_bw,
#           base = :ggplot2_base,
#           bginside = :white,
#           fg       = :black,
#           fgtext   = :lightgray,
#           fglegend = :lightgray,
#           fgguide  = :black)

_themes[:ggplot2] = PlotTheme(
    bg = :white,
    bginside = :lightgray,
    bglegend = plot_color(:lightgray, 0.8),
    framestyle = :grid,
    gridcolor = :white,
    gridalpha = 1,
    fgtext = :grey,
    fglegend = :white,
    fgguide = :black,
)


text:
mutable struct MyStruct a end
x = MyStruct(2)
x.b = 3

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia