text:
module SIMD

#=

# Various boolean types

# Idea (from <Gaunard-simd.pdf>): Use Mask{N,T} instead of booleans
# with different sizes

abstract Boolean <: Integer

for sz in (8, 16, 32, 64, 128)
    Intsz = Symbol(:Int, sz)
    UIntsz = Symbol(:UInt, sz)
    Boolsz = Symbol(:Bool, sz)
    @eval begin
        immutable $Boolsz <: Boolean
            int::$UIntsz
            $Boolsz(b::Bool) =
                new(ifelse(b, typemax($UIntsz), typemin($UIntsz)))
        end
        booltype(::Type{Val{$sz}}) = $Boolsz
        inttype(::Type{Val{$sz}}) = $Intsz
        uinttype(::Type{Val{$sz}}) = $UIntsz

        Base.convert(::Type{Bool}, b::$Boolsz) = b.int != 0

        Base. ~(b::$Boolsz) = $Boolsz(~b.int)
        Base. !(b::$Boolsz) = ~b
        Base. &(b1::$Boolsz, b2::$Boolsz) = $Boolsz(b1.int & b2.int)
        Base. |(b1::$Boolsz, b2::$Boolsz) = $Boolsz(b1.int | b2.int)
        Base.$(:$)(b1::$Boolsz, b2::$Boolsz) = $Boolsz(b1.int $ b2.int)

        Base. ==(b1::$Boolsz, b2::$Boolsz) = $Boolsz(b1.int == b2.int)
        Base. !=(b1::$Boolsz, b2::$Boolsz) = $Boolsz(b1.int != b2.int)
        Base. <(b1::$Boolsz, b2::$Boolsz) = $Boolsz(b1.int < b2.int)
        Base. <=(b1::$Boolsz, b2::$Boolsz) = $Boolsz(b1.int <= b2.int)
        Base. >(b1::$Boolsz, b2::$Boolsz) = $Boolsz(b1.int > b2.int)
        Base. >=(b1::$Boolsz, b2::$Boolsz) = $Boolsz(b1.int >= b2.int)
    end
end
Base.convert(::Type{Bool}, b::Boolean) = error("impossible")
Base.convert{I<:Integer}(::Type{I}, b::Boolean) = I(Bool(b))
Base.convert{B<:Boolean}(::Type{B}, b::Boolean) = B(Bool(b))
Base.convert{B<:Boolean}(::Type{B}, i::Integer) = B(i!=0)

booltype{T}(::Type{T}) = booltype(Val{8*sizeof(T)})
inttype{T}(::Type{T}) = inttype(Val{8*sizeof(T)})
uinttype{T}(::Type{T}) = uinttype(Val{8*sizeof(T)})

=#

# The Julia SIMD vector type

const BoolTypes = Union{Bool}
const IntTypes = Union{Int8, Int16, Int32, Int64, Int128}
const UIntTypes = Union{UInt8, UInt16, UInt32, UInt64, UInt128}
const IntegerTypes = Union{BoolTypes, IntTypes, UIntTypes}
const FloatingTypes = Union{Float16, Float32, Float64}
const ScalarTypes = Union{IntegerTypes, FloatingTypes}

const VE = Base.VecElement

export Vec
struct Vec{N,T<:ScalarTypes} # <: Number
    elts::NTuple{N,VE{T}}
    @inline Vec{N,T}(elts::NTuple{N, VE{T}}) where {N,T} = new{N,T}(elts)
end

function Base.show(io::IO, v::Vec{N,T}) where {N,T}
    print(io, "<$N x $T>[")
    join(io, [x.value for x in v.elts], ", ")
    print(io, "]")
end

# Type properties
Base.eltype(::Type{Vec{N,T}}) where {N,T} = T
Base.ndims( ::Type{Vec{N,T}}) where {N,T} = 1
Base.length(::Type{Vec{N,T}}) where {N,T} = N
Base.size(  ::Type{Vec{N,T}}) where {N,T} = (N,)
# TODO: This doesn't follow Base, e.g. `size([], 3) == 1`
Base.size(::Type{Vec{N,T}}, n::Integer) where {N,T} = (N,)[n]

Base.eltype(V::Vec) = eltype(typeof(V))
Base.ndims( V::Vec) = ndims(typeof(V))
Base.length(V::Vec) = length(typeof(V))
Base.size(  V::Vec) = size(typeof(V))
Base.size(  V::Vec, n::Integer) = size(typeof(V), n)

# Type conversion

# Create vectors from scalars or tuples
@generated function (::Type{Vec{N,T}})(x::S) where {N,T,S<:ScalarTypes}
    quote
        $(Expr(:meta, :inline))
        Vec{N,T}(tuple($([:(VE{T}(T(x))) for i in 1:N]...)))
    end
end
Vec{N,T}(xs::Tuple{}) where {N,T<:ScalarTypes} = error("illegal argument")
@generated function (::Type{Vec{N,T}})(xs::NTuple{N,S}) where {N,T,S<:ScalarTypes}
    quote
        $(Expr(:meta, :inline))
        Vec{N,T}(tuple($([:(VE{T}(T(xs[$i]))) for i in 1:N]...)))
    end
end
Vec(xs::NTuple{N,T}) where {N,T<:ScalarTypes} = Vec{N,T}(xs)

# Convert between vectors
@inline Base.convert(::Type{Vec{N,T}}, v::Vec{N,T}) where {N,T} = v

@inline Base.convert(::Type{Vec{N,R}}, v::Vec{N}) where {N,R} =
    Vec{N,R}(NTuple{N, R}(v))

@inline Tuple(v::Vec{N}) where {N} = ntuple(i -> v.elts[i].value, Val(N))
@inline NTuple{N, T}(v::Vec{N}) where{N, T} = ntuple(i -> convert(T, v.elts[i].value), Val(N))

@generated function Base. %(v::Vec{N,T}, ::Type{Vec{N,R}}) where {N,R,T}
    quote
        $(Expr(:meta, :inline))
        Vec{N,R}(tuple($([:(v.elts[$i].value % R) for i in 1:N]...)))
    end
end

# Convert vectors to tuples
@generated function Base.convert(::Type{NTuple{N,R}}, v::Vec{N,T}) where {N,R,T}
    quote
        $(Expr(:meta, :inline))
        tuple($([:(R(v.elts[$i].value)) for i in 1:N]...))
    end
end
@inline Base.convert(::Type{Tuple}, v::Vec{N,T}) where {N,T} =
    Base.convert(NTuple{N,T}, v)

# Promotion rules

# Note: Type promotion only works for subtypes of Number
# Base.promote_rule{N,T<:ScalarTypes}(::Type{Vec{N,T}}, ::Type{T}) = Vec{N,T}

Base.zero(::Type{Vec{N,T}}) where {N,T} = Vec{N,T}(zero(T))
Base.one(::Type{Vec{N,T}}) where {N,T} = Vec{N,T}(one(T))

# Floating point formats

int_type(::Type{Float16}) = Int16
int_type(::Type{Float32}) = Int32
int_type(::Type{Float64}) = Int64
# int_type(::Type{Float128}) = Int128
# int_type(::Type{Float256}) = Int256

uint_type(::Type{Float16}) = UInt16
uint_type(::Type{Float32}) = UInt32
uint_type(::Type{Float64}) = UInt64
# uint_type(::Type{Float128}) = UInt128
# uint_type(::Type{Float256}) = UInt256

significand_bits(::Type{Float16}) = 10
significand_bits(::Type{Float32}) = 23
significand_bits(::Type{Float64}) = 52
# significand_bits(::Type{Float128}) = 112
# significand_bits(::Type{Float256}) = 136

exponent_bits(::Type{T}) where {T<:FloatingTypes} =
    8*sizeof(T) - 1 - significand_bits(T)
sign_bits(::Type{T}) where {T<:FloatingTypes} = 1

significand_mask(::Type{T}) where {T<:FloatingTypes} =
    uint_type(T)(uint_type(T)(1) << significand_bits(T) - 1)
exponent_mask(::Type{T}) where {T<:FloatingTypes} =
    uint_type(T)(uint_type(T)(1) << exponent_bits(T) - 1) << significand_bits(T)
sign_mask(::Type{T}) where {T<:FloatingTypes} =
    uint_type(T)(1) << (significand_bits(T) + exponent_bits(T))

for T in (Float16, Float32, Float64)
    @assert sizeof(int_type(T)) == sizeof(T)
    @assert sizeof(uint_type(T)) == sizeof(T)
    @assert significand_bits(T) + exponent_bits(T) + sign_bits(T) == 8*sizeof(T)
    @assert significand_mask(T) | exponent_mask(T) | sign_mask(T) ==
        typemax(uint_type(T))
    @assert significand_mask(T) ⊻ exponent_mask(T) ⊻ sign_mask(T) ==
        typemax(uint_type(T))
end

# Convert Julia types to LLVM types

llvmtype(::Type{Bool}) = "i8"   # Julia represents Tuple{Bool} as [1 x i8]

# llvmtype(::Type{Bool8}) = "i8"
# llvmtype(::Type{Bool16}) = "i16"
# llvmtype(::Type{Bool32}) = "i32"
# llvmtype(::Type{Bool64}) = "i64"
# llvmtype(::Type{Bool128}) = "i128"

llvmtype(::Type{Int8}) = "i8"
llvmtype(::Type{Int16}) = "i16"
llvmtype(::Type{Int32}) = "i32"
llvmtype(::Type{Int64}) = "i64"
llvmtype(::Type{Int128}) = "i128"

llvmtype(::Type{UInt8}) = "i8"
llvmtype(::Type{UInt16}) = "i16"
llvmtype(::Type{UInt32}) = "i32"
llvmtype(::Type{UInt64}) = "i64"
llvmtype(::Type{UInt128}) = "i128"

llvmtype(::Type{Float16}) = "half"
llvmtype(::Type{Float32}) = "float"
llvmtype(::Type{Float64}) = "double"

# Type-dependent optimization flags
# fastflags{T<:IntTypes}(::Type{T}) = "nsw"
# fastflags{T<:UIntTypes}(::Type{T}) = "nuw"
# fastflags{T<:FloatingTypes}(::Type{T}) = "fast"

suffix(N::Integer, ::Type{T}) where {T<:IntegerTypes} = "v$(N)i$(8*sizeof(T))"
suffix(N::Integer, ::Type{T}) where {T<:FloatingTypes} = "v$(N)f$(8*sizeof(T))"

# Type-dependent LLVM constants
function llvmconst(::Type{T}, val) where T
    T(val) === T(0) && return "zeroinitializer"
    typ = llvmtype(T)
    "$typ $val"
end
function llvmconst(::Type{Bool}, val)
    Bool(val) === false && return "zeroinitializer"
    typ = "i1"
    "$typ $(Int(val))"
end
function llvmconst(N::Integer, ::Type{T}, val) where T
    T(val) === T(0) && return "zeroinitializer"
    typ = llvmtype(T)
    "<" * join(["$typ $val" for i in 1:N], ", ") * ">"
end
function llvmconst(N::Integer, ::Type{Bool}, val)
    Bool(val) === false && return "zeroinitializer"
    typ = "i1"
    "<" * join(["$typ $(Int(val))" for i in 1:N], ", ") * ">"
end
function llvmtypedconst(::Type{T}, val) where T
    typ = llvmtype(T)
    T(val) === T(0) && return "$typ zeroinitializer"
    "$typ $val"
end
function llvmtypedconst(::Type{Bool}, val)
    typ = "i1"
    Bool(val) === false && return "$typ zeroinitializer"
    "$typ $(Int(val))"
end

# Type-dependent LLVM intrinsics
llvmins(::Type{Val{:+}}, N, ::Type{T}) where {T <: IntegerTypes} = "add"
llvmins(::Type{Val{:-}}, N, ::Type{T}) where {T <: IntegerTypes} = "sub"
llvmins(::Type{Val{:*}}, N, ::Type{T}) where {T <: IntegerTypes} = "mul"
llvmins(::Type{Val{:div}}, N, ::Type{T}) where {T <: IntTypes} = "sdiv"
llvmins(::Type{Val{:rem}}, N, ::Type{T}) where {T <: IntTypes} = "srem"
llvmins(::Type{Val{:div}}, N, ::Type{T}) where {T <: UIntTypes} = "udiv"
llvmins(::Type{Val{:rem}}, N, ::Type{T}) where {T <: UIntTypes} = "urem"

llvmins(::Type{Val{:~}}, N, ::Type{T}) where {T <: IntegerTypes} = "xor"
llvmins(::Type{Val{:&}}, N, ::Type{T}) where {T <: IntegerTypes} = "and"
llvmins(::Type{Val{:|}}, N, ::Type{T}) where {T <: IntegerTypes} = "or"
llvmins(::Type{Val{:⊻}}, N, ::Type{T}) where {T <: IntegerTypes} = "xor"

llvmins(::Type{Val{:<<}}, N, ::Type{T}) where {T <: IntegerTypes} = "shl"
llvmins(::Type{Val{:>>>}}, N, ::Type{T}) where {T <: IntegerTypes} = "lshr"
llvmins(::Type{Val{:>>}}, N, ::Type{T}) where {T <: UIntTypes} = "lshr"
llvmins(::Type{Val{:>>}}, N, ::Type{T}) where {T <: IntTypes} = "ashr"

llvmins(::Type{Val{:(==)}}, N, ::Type{T}) where {T <: IntegerTypes} = "icmp eq"
llvmins(::Type{Val{:(!=)}}, N, ::Type{T}) where {T <: IntegerTypes} = "icmp ne"
llvmins(::Type{Val{:(>)}}, N, ::Type{T}) where {T <: IntTypes} = "icmp sgt"
llvmins(::Type{Val{:(>=)}}, N, ::Type{T}) where {T <: IntTypes} = "icmp sge"
llvmins(::Type{Val{:(<)}}, N, ::Type{T}) where {T <: IntTypes} = "icmp slt"
llvmins(::Type{Val{:(<=)}}, N, ::Type{T}) where {T <: IntTypes} = "icmp sle"
llvmins(::Type{Val{:(>)}}, N, ::Type{T}) where {T <: UIntTypes} = "icmp ugt"
llvmins(::Type{Val{:(>=)}}, N, ::Type{T}) where {T <: UIntTypes} = "icmp uge"
llvmins(::Type{Val{:(<)}}, N, ::Type{T}) where {T <: UIntTypes} = "icmp ult"
llvmins(::Type{Val{:(<=)}}, N, ::Type{T}) where {T <: UIntTypes} = "icmp ule"

llvmins(::Type{Val{:vifelse}}, N, ::Type{T}) where {T} = "select"

llvmins(::Type{Val{:+}}, N, ::Type{T}) where {T <: FloatingTypes} = "fadd"
llvmins(::Type{Val{:-}}, N, ::Type{T}) where {T <: FloatingTypes} = "fsub"
llvmins(::Type{Val{:*}}, N, ::Type{T}) where {T <: FloatingTypes} = "fmul"
llvmins(::Type{Val{:/}}, N, ::Type{T}) where {T <: FloatingTypes} = "fdiv"
llvmins(::Type{Val{:inv}}, N, ::Type{T}) where {T <: FloatingTypes} = "fdiv"
llvmins(::Type{Val{:rem}}, N, ::Type{T}) where {T <: FloatingTypes} = "frem"

llvmins(::Type{Val{:(==)}}, N, ::Type{T}) where {T <: FloatingTypes} = "fcmp oeq"
llvmins(::Type{Val{:(!=)}}, N, ::Type{T}) where {T <: FloatingTypes} = "fcmp une"
llvmins(::Type{Val{:(>)}}, N, ::Type{T}) where {T <: FloatingTypes} = "fcmp ogt"
llvmins(::Type{Val{:(>=)}}, N, ::Type{T}) where {T <: FloatingTypes} = "fcmp oge"
llvmins(::Type{Val{:(<)}}, N, ::Type{T}) where {T <: FloatingTypes} = "fcmp olt"
llvmins(::Type{Val{:(<=)}}, N, ::Type{T}) where {T <: FloatingTypes} = "fcmp ole"

llvmins(::Type{Val{:^}}, N, ::Type{T}) where {T <: FloatingTypes} =
    "@llvm.pow.$(suffix(N,T))"
llvmins(::Type{Val{:abs}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.fabs.$(suffix(N,T))"
llvmins(::Type{Val{:ceil}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.ceil.$(suffix(N,T))"
llvmins(::Type{Val{:copysign}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.copysign.$(suffix(N,T))"
llvmins(::Type{Val{:cos}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.cos.$(suffix(N,T))"
llvmins(::Type{Val{:exp}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.exp.$(suffix(N,T))"
llvmins(::Type{Val{:exp2}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.exp2.$(suffix(N,T))"
llvmins(::Type{Val{:floor}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.floor.$(suffix(N,T))"
llvmins(::Type{Val{:fma}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.fma.$(suffix(N,T))"
llvmins(::Type{Val{:log}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.log.$(suffix(N,T))"
llvmins(::Type{Val{:log10}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.log10.$(suffix(N,T))"
llvmins(::Type{Val{:log2}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.log2.$(suffix(N,T))"
llvmins(::Type{Val{:max}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.maxnum.$(suffix(N,T))"
llvmins(::Type{Val{:min}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.minnum.$(suffix(N,T))"
llvmins(::Type{Val{:muladd}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.fmuladd.$(suffix(N,T))"
llvmins(::Type{Val{:powi}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.powi.$(suffix(N,T))"
llvmins(::Type{Val{:round}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.rint.$(suffix(N,T))"
llvmins(::Type{Val{:sin}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.sin.$(suffix(N,T))"
llvmins(::Type{Val{:sqrt}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.sqrt.$(suffix(N,T))"
llvmins(::Type{Val{:trunc}}, N, ::Type{T}) where {T<:FloatingTypes} =
    "@llvm.trunc.$(suffix(N,T))"

# Convert between LLVM scalars, vectors, and arrays

function scalar2vector(vec, siz, typ, sca)
    instrs = []
    accum(nam, i) = i<0 ? "undef" : i==siz-1 ? nam : "$(nam)_iter$i"
    for i in 0:siz-1
        push!(instrs,
            "$(accum(vec,i)) = " *
                "insertelement <$siz x $typ> $(accum(vec,i-1)), " *
                "$typ $sca, i32 $i")
    end
    instrs
end

function array2vector(vec, siz, typ, arr, tmp="$(arr)_av")
    instrs = []
    accum(nam, i) = i<0 ? "undef" : i==siz-1 ? nam : "$(nam)_iter$i"
    for i in 0:siz-1
        push!(instrs, "$(tmp)_elem$i = extractvalue [$siz x $typ] $arr, $i")
        push!(instrs,
            "$(accum(vec,i)) = " *
                "insertelement <$siz x $typ> $(accum(vec,i-1)), " *
                "$typ $(tmp)_elem$i, i32 $i")
    end
    instrs
end

function vector2array(arr, siz, typ, vec, tmp="$(vec)_va")
    instrs = []
    accum(nam, i) = i<0 ? "undef" : i==siz-1 ? nam : "$(nam)_iter$i"
    for i in 0:siz-1
        push!(instrs,
            "$(tmp)_elem$i = extractelement <$siz x $typ> $vec, i32 $i")
        push!(instrs,
            "$(accum(arr,i)) = "*
                "insertvalue [$siz x $typ] $(accum(arr,i-1)), " *
                "$typ $(tmp)_elem$i, $i")
    end
    instrs
end

# TODO: change argument order
function subvector(vec, siz, typ, rvec, rsiz, roff, tmp="$(rvec)_sv")
    instrs = []
    accum(nam, i) = i<0 ? "undef" : i==rsiz-1 ? nam : "$(nam)_iter$i"
    @assert 0 <= roff
    @assert roff + rsiz <= siz
    for i in 0:rsiz-1
        push!(instrs,
            "$(tmp)_elem$i = extractelement <$siz x $typ> $vec, i32 $(roff+i)")
        push!(instrs,
            "$(accum(rvec,i)) = " *
                "insertelement <$rsiz x $typ> $(accum(rvec,i-1)), " *
                "$typ $(tmp)_elem$i, i32 $i")
    end
    instrs
end

function extendvector(vec, siz, typ, voff, vsiz, val, rvec, tmp="$(rvec)_ev")
    instrs = []
    accum(nam, i) = i<0 ? "undef" : i==siz+vsiz-1 ? nam : "$(nam)_iter$i"
    rsiz = siz + vsiz
    for i in 0:siz-1
        push!(instrs,
            "$(tmp)_elem$i = extractelement <$siz x $typ> $vec, i32 $i")
        push!(instrs,
            "$(accum(rvec,i)) = " *
                "insertelement <$rsiz x $typ> $(accum(rvec,i-1)), " *
                "$typ $(tmp)_elem$i, i32 $i")
    end
    for i in siz:siz+vsiz-1
        push!(instrs,
            "$(accum(rvec,i)) = " *
                "insertelement <$rsiz x $typ> $(accum(rvec,i-1)), $val, i32 $i")
    end
    instrs
end

# Element-wise access

export setindex
@generated function setindex(v::Vec{N,T}, x::Number, ::Type{Val{I}}) where {N,T,I}
    @assert isa(I, Integer)
    1 <= I <= N || throw(BoundsError())
    typ = llvmtype(T)
    ityp = llvmtype(Int)
    vtyp = "<$N x $typ>"
    decls = []
    instrs = []
    push!(instrs, "%res = insertelement $vtyp %0, $typ %1, $ityp $(I-1)")
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Vec{N,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{T}}, Tuple{NTuple{N,VE{T}}, T}, v.elts, T(x)))
    end
end

@generated function setindex(v::Vec{N,T}, x::Number, i::Int) where {N,T}
    typ = llvmtype(T)
    ityp = llvmtype(Int)
    vtyp = "<$N x $typ>"
    decls = []
    instrs = []
    push!(instrs, "%res = insertelement $vtyp %0, $typ %2, $ityp %1")
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        @boundscheck 1 <= i <= N || throw(BoundsError())
        Vec{N,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{T}}, Tuple{NTuple{N,VE{T}}, Int, T},
            v.elts, i-1, T(x)))
    end
end
setindex(v::Vec{N,T}, x::Number, i) where {N,T} = setindex(v, Int(i), x)

Base.getindex(v::Vec{N,T}, ::Type{Val{I}}) where {N,T,I} = v.elts[I].value
Base.getindex(v::Vec{N,T}, i) where {N,T} = v.elts[i].value

# Type conversion

@generated function Base.reinterpret(::Type{Vec{N,R}},
        v1::Vec{N1,T1}) where {N,R,N1,T1}
    @assert N*sizeof(R) == N1*sizeof(T1)
    typ1 = llvmtype(T1)
    vtyp1 = "<$N1 x $typ1>"
    typr = llvmtype(R)
    vtypr = "<$N x $typr>"
    decls = []
    instrs = []
    push!(instrs, "%res = bitcast $vtyp1 %0 to $vtypr")
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        Vec{N,R}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{R}}, Tuple{NTuple{N1,VE{T1}}}, v1.elts))
    end
end

# Generic function wrappers

# Functions taking one argument
@generated function llvmwrap(::Type{Val{Op}}, v1::Vec{N,T1},
        ::Type{R} = T1) where {Op,N,T1,R}
    @assert isa(Op, Symbol)
    typ1 = llvmtype(T1)
    vtyp1 = "<$N x $typ1>"
    typr = llvmtype(R)
    vtypr = "<$N x $typr>"
    ins = llvmins(Val{Op}, N, T1)
    decls = []
    instrs = []
    if ins[1] == '@'
        push!(decls, "declare $vtypr $ins($vtyp1)")
        push!(instrs, "%res = call $vtypr $ins($vtyp1 %0)")
    else
        if Op === :~
            @assert T1 <: IntegerTypes
            otherval = -1
        elseif Op === :inv
            @assert T1 <: FloatingTypes
            otherval = 1.0
        else
            otherval = 0
        end
        otherarg = llvmconst(N, T1, otherval)
        push!(instrs, "%res = $ins $vtyp1 $otherarg, %0")
    end
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        Vec{N,R}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{R}}, Tuple{NTuple{N,VE{T1}}}, v1.elts))
    end
end

# Functions taking one Bool argument
@generated function llvmwrap(::Type{Val{Op}}, v1::Vec{N,Bool},
        ::Type{Bool} = Bool) where {Op,N}
    @assert isa(Op, Symbol)
    btyp = llvmtype(Bool)
    vbtyp = "<$N x $btyp>"
    ins = llvmins(Val{Op}, N, Bool)
    decls = []
    instrs = []
    push!(instrs, "%arg1 = trunc $vbtyp %0 to <$N x i1>")
    otherarg = llvmconst(N, Bool, true)
    push!(instrs, "%res = $ins <$N x i1> $otherarg, %arg1")
    push!(instrs, "%resb = zext <$N x i1> %res to $vbtyp")
    push!(instrs, "ret $vbtyp %resb")
    quote
        $(Expr(:meta, :inline))
        Vec{N,Bool}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{Bool}}, Tuple{NTuple{N,VE{Bool}}}, v1.elts))
    end
end

# Functions taking two arguments
@generated function llvmwrap(::Type{Val{Op}}, v1::Vec{N,T1},
        v2::Vec{N,T2}, ::Type{R} = T1) where {Op,N,T1,T2,R}
    @assert isa(Op, Symbol)
    typ1 = llvmtype(T1)
    vtyp1 = "<$N x $typ1>"
    typ2 = llvmtype(T2)
    vtyp2 = "<$N x $typ2>"
    typr = llvmtype(R)
    vtypr = "<$N x $typr>"
    ins = llvmins(Val{Op}, N, T1)
    decls = []
    instrs = []
    if ins[1] == '@'
        push!(decls, "declare $vtypr $ins($vtyp1, $vtyp2)")
        push!(instrs, "%res = call $vtypr $ins($vtyp1 %0, $vtyp2 %1)")
    else
        push!(instrs, "%res = $ins $vtyp1 %0, %1")
    end
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        Vec{N,R}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{R}}, Tuple{NTuple{N,VE{T1}}, NTuple{N,VE{T2}}},
            v1.elts, v2.elts))
    end
end

# Functions taking two arguments, returning Bool
@generated function llvmwrap(::Type{Val{Op}}, v1::Vec{N,T1},
        v2::Vec{N,T2}, ::Type{Bool}) where {Op,N,T1,T2}
    @assert isa(Op, Symbol)
    btyp = llvmtype(Bool)
    vbtyp = "<$N x $btyp>"
    abtyp = "[$N x $btyp]"
    typ1 = llvmtype(T1)
    vtyp1 = "<$N x $typ1>"
    atyp1 = "[$N x $typ1]"
    typ2 = llvmtype(T2)
    vtyp2 = "<$N x $typ2>"
    atyp2 = "[$N x $typ2]"
    ins = llvmins(Val{Op}, N, T1)
    decls = []
    instrs = []
    if false && N == 1
        append!(instrs, array2vector("%arg1", N, typ1, "%0", "%arg1arr"))
        append!(instrs, array2vector("%arg2", N, typ2, "%1", "%arg2arr"))
        push!(instrs, "%cond = $ins $vtyp1 %arg1, %arg2")
        push!(instrs, "%res = zext <$N x i1> %cond to $vbtyp")
        append!(instrs, vector2array("%resarr", N, btyp, "%res"))
        push!(instrs, "ret $abtyp %resarr")
    else
        push!(instrs, "%res = $ins $vtyp1 %0, %1")
        push!(instrs, "%resb = zext <$N x i1> %res to $vbtyp")
        push!(instrs, "ret $vbtyp %resb")
    end
    quote
        $(Expr(:meta, :inline))
        Vec{N,Bool}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{Bool}}, Tuple{NTuple{N,VE{T1}}, NTuple{N,VE{T2}}},
            v1.elts, v2.elts))
    end
end

# Functions taking a vector and a scalar argument
# @generated function llvmwrap{Op,N,T1,T2,R}(::Type{Val{Op}}, v1::Vec{N,T1},
#         x2::T2, ::Type{R} = T1)
#     @assert isa(Op, Symbol)
#     typ1 = llvmtype(T1)
#     atyp1 = "[$N x $typ1]"
#     vtyp1 = "<$N x $typ1>"
#     typ2 = llvmtype(T2)
#     typr = llvmtype(R)
#     atypr = "[$N x $typr]"
#     vtypr = "<$N x $typr>"
#     ins = llvmins(Val{Op}, N, T1)
#     decls = []
#     instrs = []
#     append!(instrs, array2vector("%arg1", N, typ1, "%0", "%arg1arr"))
#     if ins[1] == '@'
#         push!(decls, "declare $vtypr $ins($vtyp1, $typ2)")
#         push!(instrs, "%res = call $vtypr $ins($vtyp1 %arg1, $typ2 %1)")
#     else
#         push!(instrs, "%res = $ins $vtyp1 %arg1, %1")
#     end
#     append!(instrs, vector2array("%resarr", N, typr, "%res"))
#     push!(instrs, "ret $atypr %resarr")
#     quote
#         $(Expr(:meta, :inline))
#         Vec{N,R}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
#             NTuple{N,R}, Tuple{NTuple{N,T1}, T2}, v1.elts, x2))
#     end
# end

# Functions taking two Bool arguments, returning Bool
@generated function llvmwrap(::Type{Val{Op}}, v1::Vec{N,Bool},
        v2::Vec{N,Bool}, ::Type{Bool} = Bool) where {Op,N}
    @assert isa(Op, Symbol)
    btyp = llvmtype(Bool)
    vbtyp = "<$N x $btyp>"
    ins = llvmins(Val{Op}, N, Bool)
    decls = []
    instrs = []
    push!(instrs, "%arg1 = trunc $vbtyp %0 to <$N x i1>")
    push!(instrs, "%arg2 = trunc $vbtyp %1 to <$N x i1>")
    push!(instrs, "%res = $ins <$N x i1> %arg1, %arg2")
    push!(instrs, "%resb = zext <$N x i1> %res to $vbtyp")
    push!(instrs, "ret $vbtyp %resb")
    quote
        $(Expr(:meta, :inline))
        Vec{N,Bool}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{Bool}}, Tuple{NTuple{N,VE{Bool}}, NTuple{N,VE{Bool}}},
            v1.elts, v2.elts))
    end
end

# Functions taking three arguments
@generated function llvmwrap(::Type{Val{Op}}, v1::Vec{N,T1},
        v2::Vec{N,T2}, v3::Vec{N,T3}, ::Type{R} = T1) where {Op,N,T1,T2,T3,R}
    @assert isa(Op, Symbol)
    typ1 = llvmtype(T1)
    vtyp1 = "<$N x $typ1>"
    typ2 = llvmtype(T2)
    vtyp2 = "<$N x $typ2>"
    typ3 = llvmtype(T3)
    vtyp3 = "<$N x $typ3>"
    typr = llvmtype(R)
    vtypr = "<$N x $typr>"
    ins = llvmins(Val{Op}, N, T1)
    decls = []
    instrs = []
    if ins[1] == '@'
        push!(decls, "declare $vtypr $ins($vtyp1, $vtyp2, $vtyp3)")
        push!(instrs,
            "%res = call $vtypr $ins($vtyp1 %0, $vtyp2 %1, $vtyp3 %2)")
    else
        push!(instrs, "%res = $ins $vtyp1 %0, %1, %2")
    end
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        Vec{N,R}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{R}},
            Tuple{NTuple{N,VE{T1}}, NTuple{N,VE{T2}}, NTuple{N,VE{T3}}},
            v1.elts, v2.elts, v3.elts))
    end
end

@generated function llvmwrapshift(::Type{Val{Op}}, v1::Vec{N,T},
                                  ::Type{Val{I}}) where {Op,N,T,I}
    @assert isa(Op, Symbol)
    if I >= 0
        op = Op
        i = I
    else
        if Op === :>> || Op === :>>>
            op = :<<
        else
            @assert Op === :<<
            if T <: Unsigned
                op = :>>>
            else
                op = :>>
            end
        end
        i = -I
    end
    @assert op in (:<<, :>>, :>>>)
    @assert i >= 0
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    ins = llvmins(Val{op}, N, T)
    decls = []
    instrs = []
    nbits = 8*sizeof(T)
    if (op === :>> && T <: IntTypes) || i < nbits
        count = llvmconst(N, T, min(nbits-1, i))
        push!(instrs, "%res = $ins $vtyp %0, $count")
        push!(instrs, "ret $vtyp %res")
    else
        zero = llvmconst(N, T, 0)
        push!(instrs, "return $vtyp $zero")
    end
    quote
        $(Expr(:meta, :inline))
        Vec{N,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{T}}, Tuple{NTuple{N,VE{T}}}, v1.elts))
    end
end

@generated function llvmwrapshift(::Type{Val{Op}}, v1::Vec{N,T},
                                  x2::Unsigned) where {Op,N,T}
    @assert isa(Op, Symbol)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    ins = llvmins(Val{Op}, N, T)
    decls = []
    instrs = []
    append!(instrs, scalar2vector("%count", N, typ, "%1"))
    nbits = 8*sizeof(T)
    push!(instrs, "%tmp = $ins $vtyp %0, %count")
    push!(instrs, "%inbounds = icmp ult $typ %1, $nbits")
    if Op === :>> && T <: IntTypes
        nbits1 = llvmconst(N, T, 8*sizeof(T)-1)
        push!(instrs, "%limit = $ins $vtyp %0, $nbits1")
        push!(instrs, "%res = select i1 %inbounds, $vtyp %tmp, $vtyp %limit")
    else
        zero = llvmconst(N, T, 0)
        push!(instrs, "%res = select i1 %inbounds, $vtyp %tmp, $vtyp $zero")
    end
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        # Note that this function might be called with out-of-bounds
        # values for x2, assuming that the results are then ignored
        Vec{N,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{T}}, Tuple{NTuple{N,VE{T}}, T}, v1.elts, x2 % T))
    end
end

@generated function llvmwrapshift(::Type{Val{Op}}, v1::Vec{N,T},
                                  x2::Integer) where {Op,N,T}
    if Op === :>> || Op === :>>>
        NegOp = :<<
    else
        @assert Op === :<<
        if T <: Unsigned
            NegOp = :>>>
        else
            NegOp = :>>
        end
    end
    ValOp = Val{Op}
    ValNegOp = Val{NegOp}
    quote
        $(Expr(:meta, :inline))
        ifelse(x2 >= 0,
               llvmwrapshift($ValOp, v1, unsigned(x2)),
               llvmwrapshift($ValNegOp, v1, unsigned(-x2)))
    end
end

@generated function llvmwrapshift(::Type{Val{Op}},
                                  v1::Vec{N,T},
                                  v2::Vec{N,U}) where {Op,N,T,U<:UIntTypes}
    @assert isa(Op, Symbol)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    ins = llvmins(Val{Op}, N, T)
    decls = []
    instrs = []
    push!(instrs, "%tmp = $ins $vtyp %0, %1")
    nbits = llvmconst(N, T, 8*sizeof(T))
    push!(instrs, "%inbounds = icmp ult $vtyp %1, $nbits")
    if Op === :>> && T <: IntTypes
        nbits1 = llvmconst(N, T, 8*sizeof(T)-1)
        push!(instrs, "%limit = $ins $vtyp %0, $nbits1")
        push!(instrs,
            "%res = select <$N x i1> %inbounds, $vtyp %tmp, $vtyp %limit")
    else
        zero = llvmconst(N, T, 0)
        push!(instrs,
            "%res = select <$N x i1> %inbounds, $vtyp %tmp, $vtyp $zero")
    end
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Vec{N,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{T}}, Tuple{NTuple{N,VE{T}}, NTuple{N,VE{T}}},
            v1.elts, (v2 % Vec{N,T}).elts))
    end
end

@generated function llvmwrapshift(::Type{Val{Op}},
                                  v1::Vec{N,T},
                                  v2::Vec{N,U}) where {Op,N,T,U<:IntegerTypes}
    if Op === :>> || Op === :>>>
        NegOp = :<<
    else
        @assert Op === :<<
        if T <: Unsigned
            NegOp = :>>>
        else
            NegOp = :>>
        end
    end
    ValOp = Val{Op}
    ValNegOp = Val{NegOp}
    quote
        $(Expr(:meta, :inline))
        vifelse(v2 >= 0,
                llvmwrapshift($ValOp, v1, v2 % Vec{N,unsigned(U)}),
                llvmwrapshift($ValNegOp, v1, -v2 % Vec{N,unsigned(U)}))
    end
end

# Conditionals

for op in (:(==), :(!=), :(<), :(<=), :(>), :(>=))
    @eval begin
        @inline Base.$op(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T} =
            llvmwrap(Val{$(QuoteNode(op))}, v1, v2, Bool)
    end
end
@inline function Base.isfinite(v1::Vec{N,T}) where {N,T<:FloatingTypes}
    U = uint_type(T)
    em = Vec{N,U}(exponent_mask(T))
    iv = reinterpret(Vec{N,U}, v1)
    iv & em != em
end
@inline Base.isinf(v1::Vec{N,T}) where {N,T<:FloatingTypes} = abs(v1) == Vec{N,T}(Inf)
@inline Base.isnan(v1::Vec{N,T}) where {N,T<:FloatingTypes} = v1 != v1
@inline function Base.issubnormal(v1::Vec{N,T}) where {N,T<:FloatingTypes}
    U = uint_type(T)
    em = Vec{N,U}(exponent_mask(T))
    sm = Vec{N,U}(significand_mask(T))
    iv = reinterpret(Vec{N,U}, v1)
    (iv & em == Vec{N,U}(0)) & (iv & sm != Vec{N,U}(0))
end
@inline function Base.signbit(v1::Vec{N,T}) where {N,T<:FloatingTypes}
    U = uint_type(T)
    sm = Vec{N,U}(sign_mask(T))
    iv = reinterpret(Vec{N,U}, v1)
    iv & sm != Vec{N,U}(0)
end

export vifelse
vifelse(c::Bool, x, y) = ifelse(c, x, y)
@generated function vifelse(v1::Vec{N,Bool}, v2::Vec{N,T},
        v3::Vec{N,T}) where {N,T}
    btyp = llvmtype(Bool)
    vbtyp = "<$N x $btyp>"
    abtyp = "[$N x $btyp]"
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    atyp = "[$N x $typ]"
    decls = []
    instrs = []
    if false && N == 1
        append!(instrs, array2vector("%arg1", N, btyp, "%0", "%arg1arr"))
        append!(instrs, array2vector("%arg2", N, typ, "%1", "%arg2arr"))
        append!(instrs, array2vector("%arg3", N, typ, "%2", "%arg3arr"))
        push!(instrs, "%cond = trunc $vbtyp %arg1 to <$N x i1>")
        push!(instrs, "%res = select <$N x i1> %cond, $vtyp %arg2, $vtyp %arg3")
        append!(instrs, vector2array("%resarr", N, typ, "%res"))
        push!(instrs, "ret $atyp %resarr")
    else
        push!(instrs, "%cond = trunc $vbtyp %0 to <$N x i1>")
        push!(instrs, "%res = select <$N x i1> %cond, $vtyp %1, $vtyp %2")
        push!(instrs, "ret $vtyp %res")
    end
    quote
        $(Expr(:meta, :inline))
        Vec{N,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{T}},
            Tuple{NTuple{N,VE{Bool}}, NTuple{N,VE{T}}, NTuple{N,VE{T}}},
            v1.elts, v2.elts, v3.elts))
    end
end

# Integer arithmetic functions

for op in (:~, :+, :-)
    @eval begin
        @inline Base.$op(v1::Vec{N,T}) where {N,T<:IntegerTypes} =
            llvmwrap(Val{$(QuoteNode(op))}, v1)
    end
end
@inline Base. !(v1::Vec{N,Bool}) where {N} = ~v1
@inline function Base.abs(v1::Vec{N,T}) where {N,T<:IntTypes}
    # s = -Vec{N,T}(signbit(v1))
    s = v1 >> Val{8*sizeof(T)}
    # Note: -v1 == ~v1 + 1
    (s ⊻ v1) - s
end
@inline Base.abs(v1::Vec{N,T}) where {N,T<:UIntTypes} = v1
# TODO: Try T(v1>0) - T(v1<0)
#       use a shift for v1<0
#       evaluate v1>0 as -v1<0 ?
@inline Base.sign(v1::Vec{N,T}) where {N,T<:IntTypes} =
    vifelse(v1 == Vec{N,T}(0), Vec{N,T}(0),
        vifelse(v1 < Vec{N,T}(0), Vec{N,T}(-1), Vec{N,T}(1)))
@inline Base.sign(v1::Vec{N,T}) where {N,T<:UIntTypes} =
    vifelse(v1 == Vec{N,T}(0), Vec{N,T}(0), Vec{N,T}(1))
@inline Base.signbit(v1::Vec{N,T}) where {N,T<:IntTypes} = v1 < Vec{N,T}(0)
@inline Base.signbit(v1::Vec{N,T}) where {N,T<:UIntTypes} = Vec{N,Bool}(false)

for op in (:&, :|, :⊻, :+, :-, :*, :div, :rem)
    @eval begin
        @inline Base.$op(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T<:IntegerTypes} =
            llvmwrap(Val{$(QuoteNode(op))}, v1, v2)
    end
end
@inline Base.copysign(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T<:IntTypes} =
    vifelse(signbit(v2), -abs(v1), abs(v1))
@inline Base.copysign(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T<:UIntTypes} = v1
@inline Base.flipsign(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T<:IntTypes} =
    vifelse(signbit(v2), -v1, v1)
@inline Base.flipsign(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T<:UIntTypes} = v1
@inline Base.max(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T<:IntegerTypes} =
    vifelse(v1>=v2, v1, v2)
@inline Base.min(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T<:IntegerTypes} =
    vifelse(v1>=v2, v2, v1)

@inline function Base.muladd(v1::Vec{N,T}, v2::Vec{N,T},
        v3::Vec{N,T}) where {N,T<:IntegerTypes}
    v1*v2+v3
end

# TODO: Handle negative shift counts
#       use vifelse
#       ensure vifelse is efficient
for op in (:<<, :>>, :>>>)
    @eval begin
        @inline Base.$op(v1::Vec{N,T}, ::Type{Val{I}}) where {N,T<:IntegerTypes,I} =
            llvmwrapshift(Val{$(QuoteNode(op))}, v1, Val{I})
        @inline Base.$op(v1::Vec{N,T}, x2::Unsigned) where {N,T<:IntegerTypes} =
            llvmwrapshift(Val{$(QuoteNode(op))}, v1, x2)
        @inline Base.$op(v1::Vec{N,T}, x2::Int) where {N,T<:IntegerTypes} =
            llvmwrapshift(Val{$(QuoteNode(op))}, v1, x2)
        @inline Base.$op(v1::Vec{N,T}, x2::Integer) where {N,T<:IntegerTypes} =
            llvmwrapshift(Val{$(QuoteNode(op))}, v1, x2)
        @inline Base.$op(v1::Vec{N,T},
                         v2::Vec{N,U}) where {N,T<:IntegerTypes,U<:UIntTypes} =
            llvmwrapshift(Val{$(QuoteNode(op))}, v1, v2)
        @inline Base.$op(v1::Vec{N,T},
                         v2::Vec{N,U}) where {N,T<:IntegerTypes,U<:IntegerTypes} =
            llvmwrapshift(Val{$(QuoteNode(op))}, v1, v2)
        @inline Base.$op(x1::T, v2::Vec{N,T}) where {N,T<:IntegerTypes} =
            $op(Vec{N,T}(x1), v2)
    end
end

# Floating point arithmetic functions

for op in (
        :+, :-,
        :abs, :ceil, :cos, :exp, :exp2, :floor, :inv, :log, :log10, :log2,
        :round, :sin, :sqrt, :trunc)
    @eval begin
        @inline Base.$op(v1::Vec{N,T}) where {N,T<:FloatingTypes} =
            llvmwrap(Val{$(QuoteNode(op))}, v1)
    end
end
@inline Base.exp10(v1::Vec{N,T}) where {N,T<:FloatingTypes} = Vec{N,T}(10)^v1
@inline Base.sign(v1::Vec{N,T}) where {N,T<:FloatingTypes} =
    vifelse(v1 == Vec{N,T}(0.0), Vec{N,T}(0.0), copysign(Vec{N,T}(1.0), v1))

for op in (:+, :-, :*, :/, :^, :copysign, :max, :min, :rem)
    @eval begin
        @inline Base.$op(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T<:FloatingTypes} =
            llvmwrap(Val{$(QuoteNode(op))}, v1, v2)
    end
end
@inline Base. ^(v1::Vec{N,T}, x2::Integer) where {N,T<:FloatingTypes} =
    llvmwrap(Val{:powi}, v1, Int(x2))
@inline Base.flipsign(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T<:FloatingTypes} =
    vifelse(signbit(v2), -v1, v1)

for op in (:fma, :muladd)
    @eval begin
        @inline function Base.$op(v1::Vec{N,T},
                v2::Vec{N,T}, v3::Vec{N,T}) where {N,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}, v1, v2, v3)
        end
    end
end

# Type promotion

# Promote scalars of all IntegerTypes to vectors of IntegerTypes, leaving the
# vector type unchanged

for op in (
        :(==), :(!=), :(<), :(<=), :(>), :(>=),
        :&, :|, :⊻, :+, :-, :*, :copysign, :div, :flipsign, :max, :min, :rem)
    @eval begin
        @inline Base.$op(s1::Bool, v2::Vec{N,Bool}) where {N} =
            $op(Vec{N,Bool}(s1), v2)
        @inline Base.$op(s1::IntegerTypes, v2::Vec{N,T}) where {N,T<:IntegerTypes} =
            $op(Vec{N,T}(s1), v2)
        @inline Base.$op(v1::Vec{N,T}, s2::IntegerTypes) where {N,T<:IntegerTypes} =
            $op(v1, Vec{N,T}(s2))
    end
end
@inline vifelse(c::Vec{N,Bool}, s1::IntegerTypes,
        v2::Vec{N,T}) where {N,T<:IntegerTypes} =
    vifelse(c, Vec{N,T}(s1), v2)
@inline vifelse(c::Vec{N,Bool}, v1::Vec{N,T},
        s2::IntegerTypes) where {N,T<:IntegerTypes} =
    vifelse(c, v1, Vec{N,T}(s2))

for op in (:muladd,)
    @eval begin
        @inline Base.$op(s1::IntegerTypes, v2::Vec{N,T},
                v3::Vec{N,T}) where {N,T<:IntegerTypes} =
            $op(Vec{N,T}(s1), v2, v3)
        @inline Base.$op(v1::Vec{N,T}, s2::IntegerTypes,
                v3::Vec{N,T}) where {N,T<:IntegerTypes} =
            $op(v1, Vec{N,T}(s2), v3)
        @inline Base.$op(s1::IntegerTypes, s2::IntegerTypes,
                v3::Vec{N,T}) where {N,T<:IntegerTypes} =
            $op(Vec{N,T}(s1), Vec{N,T}(s2), v3)
        @inline Base.$op(v1::Vec{N,T}, v2::Vec{N,T},
                s3::IntegerTypes) where {N,T<:IntegerTypes} =
            $op(v1, v2, Vec{N,T}(s3))
        @inline Base.$op(s1::IntegerTypes, v2::Vec{N,T},
                s3::IntegerTypes) where {N,T<:IntegerTypes} =
            $op(Vec{N,T}(s1), v2, Vec{N,T}(s3))
        @inline Base.$op(v1::Vec{N,T}, s2::IntegerTypes,
                s3::IntegerTypes) where {N,T<:IntegerTypes} =
            $op(v1, Vec{N,T}(s2), Vec{N,T}(s3))
    end
end

# Promote scalars of all ScalarTypes to vectors of FloatingTypes, leaving the
# vector type unchanged

for op in (
        :(==), :(!=), :(<), :(<=), :(>), :(>=),
        :+, :-, :*, :/, :^, :copysign, :flipsign, :max, :min, :rem)
    @eval begin
        @inline Base.$op(s1::ScalarTypes, v2::Vec{N,T}) where {N,T<:FloatingTypes} =
            $op(Vec{N,T}(s1), v2)
        @inline Base.$op(v1::Vec{N,T}, s2::ScalarTypes) where {N,T<:FloatingTypes} =
            $op(v1, Vec{N,T}(s2))
    end
end
@inline vifelse(c::Vec{N,Bool}, s1::ScalarTypes,
        v2::Vec{N,T}) where {N,T<:FloatingTypes} =
    vifelse(c, Vec{N,T}(s1), v2)
@inline vifelse(c::Vec{N,Bool}, v1::Vec{N,T},
        s2::ScalarTypes) where {N,T<:FloatingTypes} =
    vifelse(c, v1, Vec{N,T}(s2))

for op in (:fma, :muladd)
    @eval begin
        @inline Base.$op(s1::ScalarTypes, v2::Vec{N,T},
                v3::Vec{N,T}) where {N,T<:FloatingTypes} =
            $op(Vec{N,T}(s1), v2, v3)
        @inline Base.$op(v1::Vec{N,T}, s2::ScalarTypes,
                v3::Vec{N,T}) where {N,T<:FloatingTypes} =
            $op(v1, Vec{N,T}(s2), v3)
        @inline Base.$op(s1::ScalarTypes, s2::ScalarTypes,
                v3::Vec{N,T}) where {N,T<:FloatingTypes} =
            $op(Vec{N,T}(s1), Vec{N,T}(s2), v3)
        @inline Base.$op(v1::Vec{N,T}, v2::Vec{N,T},
                s3::ScalarTypes) where {N,T<:FloatingTypes} =
            $op(v1, v2, Vec{N,T}(s3))
        @inline Base.$op(s1::ScalarTypes, v2::Vec{N,T},
                s3::ScalarTypes) where {N,T<:FloatingTypes} =
            $op(Vec{N,T}(s1), v2, Vec{N,T}(s3))
        @inline Base.$op(v1::Vec{N,T}, s2::ScalarTypes,
                s3::ScalarTypes) where {N,T<:FloatingTypes} =
            $op(v1, Vec{N,T}(s2), Vec{N,T}(s3))
    end
end

# Reduction operations

# TODO: map, mapreduce

function getneutral(op::Symbol, ::Type{T}) where T
    zs = Dict{Symbol,T}()
    if T <: IntegerTypes
        zs[:&] = ~T(0)
        zs[:|] = T(0)
    end
    zs[:max] = typemin(T)
    zs[:min] = typemax(T)
    zs[:+] = T(0)
    zs[:*] = T(1)
    zs[op]
end

if VERSION >= v"0.7.0-beta2.195"
    nextpow2(n) = nextpow(2, n)
end

# We cannot pass in the neutral element via Val{}; if we try, Julia refuses to
# inline this function, which is then disastrous for performance
@generated function llvmwrapreduce(::Type{Val{Op}}, v::Vec{N,T}) where {Op,N,T}
    @assert isa(Op, Symbol)
    z = getneutral(Op, T)
    typ = llvmtype(T)
    decls = []
    instrs = []
    n = N
    nam = "%0"
    nold,n = n,nextpow2(n)
    if n > nold
        namold,nam = nam,"%vec_$n"
        append!(instrs,
            extendvector(namold, nold, typ, n, n-nold,
                llvmtypedconst(T,z), nam))
    end
    while n > 1
        nold,n = n, div(n, 2)
        namold,nam = nam,"%vec_$n"
        vtyp = "<$n x $typ>"
        ins = llvmins(Val{Op}, n, T)
        append!(instrs, subvector(namold, nold, typ, "$(nam)_1", n, 0))
        append!(instrs, subvector(namold, nold, typ, "$(nam)_2", n, n))
        if ins[1] == '@'
            push!(decls, "declare $vtyp $ins($vtyp, $vtyp)")
            push!(instrs,
                "$nam = call $vtyp $ins($vtyp $(nam)_1, $vtyp $(nam)_2)")
        else
            push!(instrs, "$nam = $ins $vtyp $(nam)_1, $(nam)_2")
        end
    end
    push!(instrs, "%res = extractelement <$n x $typ> $nam, i32 0")
    push!(instrs, "ret $typ %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            T, Tuple{NTuple{N,VE{T}}}, v.elts)
    end
end

@inline Base.all(v::Vec{N,T}) where {N,T<:IntegerTypes} = llvmwrapreduce(Val{:&}, v)
@inline Base.any(v::Vec{N,T}) where {N,T<:IntegerTypes} = llvmwrapreduce(Val{:|}, v)
@inline Base.maximum(v::Vec{N,T}) where {N,T<:FloatingTypes} =
    llvmwrapreduce(Val{:max}, v)
@inline Base.minimum(v::Vec{N,T}) where {N,T<:FloatingTypes} =
    llvmwrapreduce(Val{:min}, v)
@inline Base.prod(v::Vec{N,T}) where {N,T} = llvmwrapreduce(Val{:*}, v)
@inline Base.sum(v::Vec{N,T}) where {N,T} = llvmwrapreduce(Val{:+}, v)

@generated function Base.reduce(::Type{Val{Op}}, v::Vec{N,T}) where {Op,N,T}
    @assert isa(Op, Symbol)
    z = getneutral(Op, T)
    stmts = []
    n = N
    push!(stmts, :($(Symbol(:v,n)) = v))
    nold,n = n,nextpow2(n)
    if n > nold
        push!(stmts,
            :($(Symbol(:v,n)) = Vec{$n,T}($(Expr(:tuple,
                [:($(Symbol(:v,nold)).elts[$i]) for i in 1:nold]...,
                [z for i in nold+1:n]...)))))
    end
    while n > 1
        nold,n = n, div(n, 2)
        push!(stmts,
            :($(Symbol(:v,n,"lo")) = Vec{$n,T}($(Expr(:tuple,
                [:($(Symbol(:v,nold)).elts[$i]) for i in 1:n]...,)))))
        push!(stmts,
            :($(Symbol(:v,n,"hi")) = Vec{$n,T}($(Expr(:tuple,
                [:($(Symbol(:v,nold)).elts[$i]) for i in n+1:nold]...)))))
        push!(stmts,
            :($(Symbol(:v,n)) =
                $Op($(Symbol(:v,n,"lo")), $(Symbol(:v,n,"hi")))))
    end
    push!(stmts, :(v1[1]))
    Expr(:block, Expr(:meta, :inline), stmts...)
end

@inline Base.maximum(v::Vec{N,T}) where {N,T<:IntegerTypes} = reduce(Val{:max}, v)
@inline Base.minimum(v::Vec{N,T}) where {N,T<:IntegerTypes} = reduce(Val{:min}, v)

# Load and store functions

export valloc
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

export vload, vloada
@generated function vload(::Type{Vec{N,T}}, ptr::Ptr{T},
                          ::Type{Val{Aligned}} = Val{false}) where {N,T,Aligned}
    @assert isa(Aligned, Bool)
    ptyp = llvmtype(Int)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    decls = []
    instrs = []
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    flags = [""]
    if align > 0
        push!(flags, "align $align")
    end
    if VERSION < v"v0.7.0-DEV"
        push!(instrs, "%ptr = bitcast $typ* %0 to $vtyp*")
    else
        push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
    end
    push!(instrs, "%res = load $vtyp, $vtyp* %ptr" * join(flags, ", "))
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Vec{N,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{T}}, Tuple{Ptr{T}}, ptr))
    end
end

@inline vloada(::Type{Vec{N,T}}, ptr::Ptr{T}) where {N,T} =
    vload(Vec{N,T}, ptr, Val{true})

@inline function vload(::Type{Vec{N,T}},
                       arr::Union{Array{T,1},SubArray{T,1}},
                       i::Integer,
                       ::Type{Val{Aligned}} = Val{false}) where {N,T,Aligned}
    #TODO @boundscheck 1 <= i <= length(arr) - (N-1) || throw(BoundsError())
    vload(Vec{N,T}, pointer(arr, i), Val{Aligned})
end
@inline function vloada(::Type{Vec{N,T}},
                        arr::Union{Array{T,1},SubArray{T,1}},
                        i::Integer) where {N,T}
    vload(Vec{N,T}, arr, i, Val{true})
end

@generated function vload(::Type{Vec{N,T}}, ptr::Ptr{T},
                          mask::Vec{N,Bool},
                          ::Type{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 = []
    instrs = []
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end

    if VERSION < v"v0.7.0-DEV"
        push!(instrs, "%ptr = bitcast $typ* %0 to $vtyp*")
    else
        push!(instrs, "%ptr = inttoptr $ptyp %0 to $vtyp*")
    end
    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))
        Vec{N,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{T}}, Tuple{Ptr{T}, NTuple{N,VE{Bool}}}, ptr, mask.elts))
    end
end

@inline vloada(::Type{Vec{N,T}}, ptr::Ptr{T}, mask::Vec{N,Bool}) where {N,T} =
    vload(Vec{N,T}, ptr, mask, Val{true})

@inline function vload(::Type{Vec{N,T}},
                       arr::Union{Array{T,1},SubArray{T,1}},
                       i::Integer, mask::Vec{N,Bool},
                       ::Type{Val{Aligned}} = Val{false}) where {N,T,Aligned}
    #TODO @boundscheck 1 <= i <= length(arr) - (N-1) || throw(BoundsError())
    vload(Vec{N,T}, pointer(arr, i), mask, Val{Aligned})
end
@inline function vloada(::Type{Vec{N,T}},
                        arr::Union{Array{T,1},SubArray{T,1}}, i::Integer,
                        mask::Vec{N,Bool}) where {N,T}
    vload(Vec{N,T}, arr, i, mask, Val{true})
end

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

@inline vstorea(v::Vec{N,T}, ptr::Ptr{T}) where {N,T} = vstore(v, ptr, Val{true})

@inline function vstore(v::Vec{N,T},
                        arr::Union{Array{T,1},SubArray{T,1}},
                        i::Integer,
                        ::Type{Val{Aligned}} = Val{false}) where {N,T,Aligned}
    @boundscheck 1 <= i <= length(arr) - (N-1) || throw(BoundsError())
    vstore(v, pointer(arr, i), Val{Aligned})
end
@inline function vstorea(v::Vec{N,T}, arr::Union{Array{T,1},SubArray{T,1}},
                         i::Integer) where {N,T}
    vstore(v, arr, i, Val{true})
end

@generated function vstore(v::Vec{N,T}, ptr::Ptr{T},
                           mask::Vec{N,Bool},
                           ::Type{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 = []
    instrs = []
    if Aligned
        align = N * sizeof(T)
    else
        align = sizeof(T)   # This is overly optimistic
    end
    if VERSION < v"v0.7.0-DEV"
        push!(instrs, "%ptr = bitcast $typ* %1 to $vtyp*")
    else
        push!(instrs, "%ptr = inttoptr $ptyp %1 to $vtyp*")
    end
    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 %0, $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{NTuple{N,VE{T}}, Ptr{T}, NTuple{N,VE{Bool}}},
            v.elts, ptr, mask.elts)
    end
end

@inline vstorea(v::Vec{N,T}, ptr::Ptr{T}, mask::Vec{N,Bool}) where {N,T} =
    vstore(v, ptr, mask, Val{true})

@inline function vstore(v::Vec{N,T},
                        arr::Union{Array{T,1},SubArray{T,1}},
                        i::Integer,
                        mask::Vec{N,Bool},
                        ::Type{Val{Aligned}} = Val{false}) where {N,T,Aligned}
    #TODO @boundscheck 1 <= i <= length(arr) - (N-1) || throw(BoundsError())
    vstore(v, pointer(arr, i), mask, Val{Aligned})
end
@inline function vstorea(v::Vec{N,T},
                         arr::Union{Array{T,1},SubArray{T,1}},
                         i::Integer, mask::Vec{N,Bool}) where {N,T}
    vstore(v, arr, i, mask, Val{true})
end

# Vector shuffles

function shufflevector_instrs(N, T, I, two_operands)
    typ = llvmtype(T)
    vtyp2 = vtyp1 = "<$N x $typ>"
    M = length(I)
    vtyp3 = "<$M x i32>"
    vtypr = "<$M x $typ>"
    mask = "<" * join(map(x->string("i32 ", x), I), ", ") * ">"
    instrs = []
    v2 = two_operands ? "%1" : "undef"
    push!(instrs, "%res = shufflevector $vtyp1 %0, $vtyp2 $v2, $vtyp3 $mask")
    push!(instrs, "ret $vtypr %res")
    return M, [], instrs
end

export shufflevector
@generated function shufflevector(v1::Vec{N,T}, v2::Vec{N,T},
                                  ::Type{Val{I}}) where {N,T,I}
    M, decls, instrs = shufflevector_instrs(N, T, I, true)
    quote
        $(Expr(:meta, :inline))
        Vec{$M,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{$M,VE{T}},
            Tuple{NTuple{N,VE{T}}, NTuple{N,VE{T}}},
            v1.elts, v2.elts))
    end
end

@generated function shufflevector(v1::Vec{N,T}, ::Type{Val{I}}) where {N,T,I}
    M, decls, instrs = shufflevector_instrs(N, T, I, false)
    quote
        $(Expr(:meta, :inline))
        Vec{$M,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{$M,VE{T}},
            Tuple{NTuple{N,VE{T}}},
            v1.elts))
    end
end

end


text:
push!(LOAD_PATH,"../src/")
using Documenter, FourierAnalysis

makedocs(
   sitename="FourierAnalysis",
   authors="Marco Congedo, CNRS, France",
   modules=[FourierAnalysis],
   pages =  [
      "index.md",
      "Main Module" => "MainModule.md",
      "Tapering Window" => "tapers.md",
      "frequency domain" => Any[
                           "Spectral Estimations" => "spectra.md",
                           "Cross-Spectral Matrices" => "crossspectra.md",
                           "Coherence Matrices" => "coherence.md",
                           "Goertzel's Algorithms" => "goertzel.md",
      ],
      "time-frequency(TF) domain" => Any[
                           "TF Representations" => "timefrequency.md",
                           "TF Univariate Measures" => "timefrequencyuni.md",
                           "TF Bivariate Measures " => "timefrequencybi.md",
      ],
      "utilities" => Any[
                        "Plots" => "plots.md",
                        "Tools" => "tools.md",
                        "FFTW planners"  => "fftw.md",
                        "Filter Banks" => "filters.md",
                        "Hilbert Transform" => "hilbert.md",
      ]
   ]
)

# fix

deploydocs(
   # root
   # target = "build", # add this folder to .gitignore!
   repo = "github.com/Marco-Congedo/FourierAnalysis.jl.git",
   # branch = "gh-pages",
   # osname = "linux",
   # deps = Deps.pip("pygments", "mkdocs"),
   # devbranch = "dev",
   # devurl = "dev",
   # versions = ["stable" => "v^", "v#.#", devurl => devurl],
)

ext: jl


text:
# Script which produces the overview map
#
# Requires command line tools:
# - unzip
# - gdal_merge.py
# used coordinate system: CH1903+/LV95

using Rhone2020
const R20 = Rhone2020

import PyPlot
const Plt = PyPlot
#Plt.pygui(true) # figure appears in window
#Plt.pygui(false) # figure does not appear in window

using GeoData, Shapefile, ArchGDAL, DelimitedFiles, LaTeXStrings

# ------------------------------------------ #
#             Download data directly from SwissTopo & unzip          #
#             not necessary anymore as data is in data submission
# ------------------------------------------ #
#
# # surface DEM swisstopo SwissAlti3D (2m posting)
# files_surf = ["https://data.geo.admin.ch/ch.swisstopo.swissalti3d/swissalti3d_2019_267$(i)-116$(j)/swissalti3d_2019_267$(i)-116$(j)_2_2056_5728.tif" for i=1:3, j=0:2][:]
#
## Ortho photos
# files_ortho_highres = ["https://data.geo.admin.ch/ch.swisstopo.swissimage-dop10/swissimage-dop10_2020_267$(i)-116$(j)/swissimage-dop10_2020_267$(i)-116$(j)_0.1_2056.tif" for i=2:2, j=0:1][:]
# files_ortho = ["https://data.geo.admin.ch/ch.swisstopo.swissimage-dop10/swissimage-dop10_2020_267$(i)-116$(j)/swissimage-dop10_2020_267$(i)-116$(j)_2_2056.tif" for i=2:2, j=0:1][:]
#
# # outline of Switzerland
# # https://www.swisstopo.admin.ch/de/geodata/landscape/boundaries3d.html
# zipfile_CH = "https://cms.geo.admin.ch/ogd/topography/swissBOUNDARIES3D.zip"
# folder_CH = "BOUNDARIES_2021/DATEN/swissBOUNDARIES3D/SHAPEFILE_LV95_LN02/swissBOUNDARIES3D_1_3_TLM_LANDESGEBIET.*"
#
# # Bed DEM and outline from SwissGlacierThickness-R2020
# # https://www.research-collection.ethz.ch/handle/20.500.11850/434697
# #
# # Rhone is SGI B43/03
# zipfile_outline = "https://www.research-collection.ethz.ch/bitstream/handle/20.500.11850/434697/01_Summary.zip"
# zipfile_bed = "https://www.research-collection.ethz.ch/bitstream/handle/20.500.11850/434697/07_GlacierBed_SGI.zip"
# file_bed = "07_GlacierBed_SGI/B43-03_GlacierBed.tif"
#
# # download
# destdir = mkpath(joinpath(R20.datadir, "geodata/swisstopo/tmp-surf"))
# # zipfile_bed = R20.download_file_to_dir(zipfile_bed, destdir)
# # zipfile_outline = R20.download_file_to_dir(zipfile_outline, destdir)
# # # merge tifs into one
# files_surf = [R20.download_file_to_dir(f, destdir) for f in files_ortho]
# file_surf = joinpath(destdir, "surface_DEM.tif")
# run(`gdal_merge.py -of GTiff -o $file_surf $files_surf`)

# # # # files_surf = [R20.download_file_to_dir(f, destdir) for f in files_surf]
# destdir = mkpath(joinpath(R20.datadir, "geodata/swisstopo/tmp-ortho"))
# files_ortho = [R20.download_file_to_dir(f, destdir) for f in files_ortho]
# # # merge tifs into one
# file_ortho = joinpath(destdir, "ortho.tif")
# run(`gdal_merge.py -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -of GTiff -o $file_ortho $files_ortho`)

# zipfile_CH = R20.download_file_to_dir(zipfile_CH, destdir)
#
# # unzip (requires unzip)
# run(`unzip -o $zipfile_outline -d $destdir`)
# shpfile_outline = joinpath(destdir, "01_Summary/Summary.shp")
#
# run(`unzip -o $zipfile_bed $file_bed -d $destdir`)
# file_bed = joinpath(destdir, file_bed)
#
# run(`unzip -o $zipfile_CH $folder_CH -d $destdir`) # why is this still in the folders?
# shpfile_CH = joinpath(destdir, "BOUNDARIES_2021/DATEN/swissBOUNDARIES3D/SHAPEFILE_LV95_LN02/swissBOUNDARIES3D_1_3_TLM_LANDESGEBIET.shp")

# ------------------------------------------ #
#          Load data from data folder        #
# ------------------------------------------ #

fielddir = joinpath(@__DIR__, "../../data/geodata/field_data/")
ethdir = joinpath(@__DIR__, "../../data/geodata/research_collection_ethz/")
swisstopodir = joinpath(@__DIR__, "../../data/geodata/swisstopo/")

# shapefile glacier outline
shpfile_outline = ethdir * "outline_rhonegletscher/Summary.shp"
table = Shapefile.Table(shpfile_outline)
outline = nothing
for row in table
    if row.pk_sgi=="B43/03"
        global outline = Shapefile.shape(row)
    end
end

# shapefile CH boundary
shpfile_CH = swisstopodir * "outline_CH/swissBOUNDARIES3D_1_3_TLM_LANDESGEBIET.shp"
table = Shapefile.Table(shpfile_CH)
boundary = nothing
for row in table
    if row.NAME=="Schweiz"
        global boundary = Shapefile.shape(row)
    end
end

# shapefile streams
shpfile_stream1 = fielddir * "streams/stream_AM13_15.shp"
table = Shapefile.Table(shpfile_stream1)
stream1 = nothing
for row in table
    global stream1 = Shapefile.shape(row)
end

shpfile_stream2 = fielddir * "streams/stream_AM14.shp"
table = Shapefile.Table(shpfile_stream2)
stream2 = nothing
for row in table
    global stream2 = Shapefile.shape(row)
end

# bed dem
file_bed = ethdir * "B43-03_GlacierBed.tif"
bed = GDALarray(file_bed)[:,:,1] # the indexing brings it into memory

# surface dem
file_surf = swisstopodir * "surface_DEM.tif"
surf = GDALarray(file_surf)[1:5:end,1:5:end,1]
# thick = surf.-bed # does not work
# smooth it
function boxcar(A, window)
    out = similar(A)
    R = CartesianIndices(size(A))
    I1, Iend = first(R), last(R)
    for I in R
        n, s = 0, zero(eltype(out))
        I_ul = CartesianIndex(I1.I .* window)
        for J in CartesianIndices(UnitRange.(max(I1, I-I_ul).I , min(Iend, I+I_ul).I) )
            s += A[J]
            n += 1
        end
        out[I] = s/n
    end
    out
end
surf = boxcar(surf, 3)
#thick = boxcar(thick, 10)

# orthophoto (load as Array{Uint8, 3} for pyplot)
file_ortho = swisstopodir * "ortho.tif"
#ortho = Plt.imread(file_ortho)
ortho = ArchGDAL.readraster(file_ortho)
thin = 1
ortho = ortho[1:thin:end, 1:thin:end,:];

# borehole coordinates
BH_file = fielddir * "BH_coordinates.csv"
BH_data = readdlm(BH_file, ','; skipstart=1)
BH_easting = BH_data[:,2] .*1e-3 # in km
BH_northing = BH_data[:,3] .*1e-3 # in km

#######
## Plotting
#
# https://discourse.julialang.org/t/labels-for-levels-in-contour-plot-plots-jl

# cl = :grey
# p = contour(bed, levels=2000:50:3000,  ls=:dash,
#             contour_labels=true,
#             #color = cgrad([:black,:black,:black]),
#             color = cgrad([cl, cl, cl]),
#             aspect_ratio = :equal,
#             colorbar=false
#             )
# contour!(surf, levels=2000:50:3000,
#          contour_labels=true,
#          color = cgrad([:black,:black,:black]),
#          colorbar=false,
#          )
# plot!(outline, fillalpha=0,
#       lw=4,
#       #color=:blue,
#       xlim=(2671454, 2673415),
#       ylim=(1160739, 1162721),
#       ylabel = "Northing (m)",
#       xlabel = "Easting (m)",
#       title="",
#       axes=:box,
#       )

# using PyPlot instead
# https://discourse.julialang.org/t/labels-for-levels-in-contour-plot-plots-jl/8266/3

cl = :grey
bed[bed.==-9999] .= NaN
levels_surface = 2000.0:20:3000
level_labels_surf = map(x -> "$(round(Int,x))", levels_surface)
level_ind_surf = 1:1:length(level_labels_surf)
levels_bed = 2000.0:50:3000
level_labels_bed = map(x -> "$(round(Int,x))", levels_bed)
level_ind_bed = 1:1:length(level_labels_bed)

Plt.figure(figsize=(14, 6), dpi=300)

# font sizes
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["font.size"] = 14 # general
fs_panellabels = 19 # for panel labels a-c

# positions and sizes of figures
xmin_CH, xmax_CH, ymin_CH, ymax_CH = 2468.0, 2851.0, 1064.0, 1307.0
xmin_Rhone, xmax_Rhone, ymin_Rhone, ymax_Rhone = 2671.2, 2675.6, 1158.8, 1167.8
ratio_CH = (ymax_CH-ymin_CH)/(xmax_CH-xmin_CH)
ratio_Rhone = (ymax_Rhone-ymin_Rhone)/(xmax_Rhone-xmin_Rhone)

left = 0.05
lpad = 0.082
bottom = 0.1
width = 0.12
heightCH = width*ratio_CH*14/6
heightRhone = width*ratio_Rhone*14/6

# surface topography
ax = Plt.axes([left, bottom, 0.9, 0.85])
cs = Plt.contour(surf.dims[1].val./1e3, surf.dims[2].val./1e3, surf',
                levels=levels_surface, colors="k", linewidths=1)
Plt.xlim((2671.5, 2673.5))
Plt.ylim((1160.8, 1162.8))

Plt.clabel(cs, levels_surface[level_ind_surf], fmt = Dict(zip(levels_surface[level_ind_surf], level_labels_surf[level_ind_surf])),
           inline =true,
           inline_spacing=1)

# bed topography
cs = Plt.contour(bed.dims[1].val./1e3, bed.dims[2].val./1e3, bed',
                levels=levels_bed, colors="grey", linewidths=1)
Plt.xlim((2671.5, 2673.5))
Plt.ylim((1160.8, 1162.8))
Plt.clabel(cs, levels_bed[level_ind_bed], fmt = Dict(zip(levels_bed[level_ind_bed], level_labels_bed[level_ind_bed])),
           inline = 1,
           inline_spacing=1)

# glacier outline
outl = hcat([[p.x,p.y] for p in outline.points]...);
Plt.plot(outl[1,:]/1e3, outl[2,:]/1e3, lw=3)

# BH/AM coordinates
p1 = Plt.scatter(BH_easting[[1,2,4]], BH_northing[[1,2,4]], 30, "k", "o") # unactivated BH
p2 = Plt.scatter(BH_easting[[3,5]], BH_northing[[3,5]], 50, "k", "x") # AM

Plt.annotate("BH11", (BH_easting[1]-0.085, BH_northing[1]-0.005))
Plt.annotate("BH12", (BH_easting[2]-0.085, BH_northing[2]-0.01))
Plt.annotate("AM13/BH13", (BH_easting[3]-0.15, BH_northing[3]-0.005))
Plt.annotate("BH14", (BH_easting[4]+0.02, BH_northing[4]-0.015))
Plt.annotate("AM15/BH15", (BH_easting[5]-0.15, BH_northing[5]))

# stream
str1 = hcat([[p.x,p.y] for p in stream1.points]...);
str2 = hcat([[p.x,p.y] for p in stream2.points]...);
p4 = Plt.plot(str1[1,:]/1e3, str1[2,:]/1e3, "k--", lw=1)
Plt.plot(str2[1,:]/1e3, str2[2,:]/1e3, "k--", lw=1)
p3 = Plt.scatter(str1[1,1]/1e3, str1[2,1]/1e3, 30, "k", "<")
Plt.scatter(str2[1,1]/1e3, str2[2,1]/1e3, 30, "k", "<")

# workaround to make legend for 2D lines
L2D = Plt.matplotlib.lines.Line2D
custom_lines = [L2D([0], [0], color="k", linewidth=1, linestyle="--"),
                L2D([0], [0], color="k", linewidth=1),
                L2D([0], [0], color="grey", linewidth=1),
                L2D([0], [0], linewidth=3)]
Plt.legend((custom_lines[1], custom_lines[2], custom_lines[3], custom_lines[4],
            p1, p2, p3),
            ("Surface streams", L"\mathrm{Surface\,\,contour\,\,lines,\,\,20\,m\,\,spacing}", L"\mathrm{Bed\,\,contour\,\,lines,\,\,50\,m\,\,spacing}", "Glacier outline",
            "Not activated boreholes", "Artificial moulins", "Natural moulins"),
            framealpha = 1.0,
            loc = "lower right",
            ncol = 2,
            fontsize = 15,
            )

ymin, ymax = 1160.63, 1161.23
xmin, xmax = 2671.55, 2672.8
Plt.axis("equal")
Plt.axis("square")
Plt.xlim((xmin, xmax))
Plt.ylim((ymin, ymax))
Plt.xlabel("Easting (km)")
Plt.ylabel("Northing (km)")
# Plt.tick_params(labelleft=false, left=false, labelbottom=false, bottom=false)
Plt.text(0.88, 0.98, L"\bf{c}", transform=ax.transAxes, ha="left", va="top", fontsize=fs_panellabels+2) #
# Plt.hlines(ymax-0.25*height1/(height1+height2)*(ymax-ymin), xmin+0.895, xmin+0.995,
#           "k", linewidth=4)
# Plt.annotate(L"\mathrm{100\,m}", (xmin+0.92, ymax-0.22*height1/(height1+height2)*(ymax-ymin)), fontsize=12)
# Plt.arrow(xmin+0.945, ymin+0.38, 0, 0.055,
#          facecolor="black",
#          width=0.002,
#          head_width=0.01,
#          head_length=0.015
#          )
# Plt.annotate("N", (xmin+0.925, ymin+0.39), fontsize=28)
Plt.minorticks_off()

# glacier outline alone
ax = Plt.axes([left+lpad, bottom+(0.85-heightCH-heightRhone)*1/3, width, heightRhone]) #figsize=(12, 7))
Plt.plot(outl[1,1:1364]/1e3, outl[2,1:1364]/1e3, lw=2, zorder=1) # exclude some points at the end which are within the glacier and make outline ugly
Plt.hlines([ymin, ymax], xmin, xmax, "k", lw=2, zorder=2)
Plt.vlines([xmin, xmax], ymin, ymax, "k", lw=2, zorder=2)
Plt.axis("equal")
Plt.axis("square")
Plt.tick_params(labelleft=false, left=false, labelbottom=false, bottom=false)
Plt.xlim([xmin_Rhone, xmax_Rhone])
Plt.ylim([ymin_Rhone, ymax_Rhone])
Plt.text(xmin+0.5*(xmax-xmin), ymax+0.3*(ymax-ymin), "c", fontsize=fs_panellabels, ha="center")
Plt.text(0.06, 0.98, L"$\bf{b}$", transform=ax.transAxes, ha="left", va="top", fontsize=fs_panellabels)
Plt.hlines(ymin_Rhone+0.15*(ymax_Rhone-ymin_Rhone), xmin_Rhone+0.6*(xmax_Rhone-xmin_Rhone), xmin_Rhone+0.6*(xmax_Rhone-xmin_Rhone)+1.0,
           "k", linewidth=4)
Plt.annotate(L"\mathrm{1\,km}", (xmin_Rhone+0.6*(xmax_Rhone-xmin_Rhone), ymin_Rhone+0.18*(ymax_Rhone-ymin_Rhone)), fontsize=12)
Plt.arrow(0.55, 0.8, 0., -0.1, width=0.01, head_width=0.05, head_length=0.03, transform=ax.transAxes, facecolor="black")
Plt.text(0.65, 0.75, "ice flow", transform=ax.transAxes, ha="left", va="center", rotation="vertical")

# CH outline alone
ax = Plt.axes([left+lpad, 0.95-(0.85-heightCH-heightRhone)*2/3-heightCH, width, heightCH])
bCH = hcat([[p.x,p.y] for p in boundary.points]...);
Plt.plot(bCH[1,1:end-414]/1e3, bCH[2,1:end-414]/1e3, "k") # exclude some points, same reason as in glacier outline
Plt.scatter(BH_easting[1], BH_northing[1])
Plt.axis("equal")
Plt.tick_params(labelleft=false, left=false, labelbottom=false, bottom=false)
Plt.xlim([xmin_CH, xmax_CH])
Plt.ylim([ymin_CH, ymax_CH])
ymin_CH, ymax_CH = Plt.ylim()
xmin_CH, xmax_CH = Plt.xlim()
Plt.text(0.06, 1-0.02*heightRhone/heightCH, L"$\bf{a}$", transform=ax.transAxes, ha="left", va="top", fontsize=fs_panellabels)
Plt.text(BH_easting[1]-15, BH_northing[1]+30.0, "b", fontsize=fs_panellabels)

Plt.gcf()
Plt.savefig("../../products/paper_figures/figure1.png")

###
# Plot orthoimage and associated data
###
Plt.figure(figsize=(6*9/6,8*9/6), dpi=400)
Plt.imshow(permutedims(ortho, [2,1,3]), aspect="equal",
           extent= [2672, 2673, 1160, 1162])
# glacier outline
outl = hcat([[p.x,p.y] for p in outline.points]...);
Plt.plot(outl[1,:]/1e3, outl[2,:]/1e3, lw=3)

# BH/AM coordinates
p1 = Plt.scatter(BH_easting[[1,2,4]], BH_northing[[1,2,4]], 30, "k", "o") # unactivated BH
p2 = Plt.scatter(BH_easting[[3,5]], BH_northing[[3,5]], 50, "k", "x") # AM

Plt.annotate("BH11", (BH_easting[1]-0.0335, BH_northing[1]-0.005))
Plt.annotate("BH12", (BH_easting[2]-0.035, BH_northing[2]-0.01))
Plt.annotate("AM13/BH13", (BH_easting[3]-0.07, BH_northing[3]-0.005))
Plt.annotate("BH14", (BH_easting[4]+0.01, BH_northing[4]))
Plt.annotate("AM15/BH15", (BH_easting[5]-0.07, BH_northing[5]))

# stream
str1 = hcat([[p.x,p.y] for p in stream1.points]...);
str2 = hcat([[p.x,p.y] for p in stream2.points]...);
p4 = Plt.plot(str1[1,:]/1e3, str1[2,:]/1e3, "k--", lw=1)
Plt.plot(str2[1,:]/1e3, str2[2,:]/1e3, "k--", lw=1)
p3 = Plt.scatter(str1[1,1]/1e3, str1[2,1]/1e3, 30, "k", "<")
Plt.scatter(str2[1,1]/1e3, str2[2,1]/1e3, 30, "k", "<")

# full extent
# Plt.xlim([2672, 2673])
# Plt.ylim([1160, 1161.6])

Plt.xlim([2672.2, 2672.5])
Plt.ylim([1160.75, 1161.25])
Plt.xticks(2672.2:0.1:2672.5,2672.2:0.1:2672.5)
Plt.xlabel("Easting (km)")
Plt.ylabel("Northing (km)")

Plt.savefig("../../products/paper_figures/figure_s1.png", bbox_inches="tight")

Plt.close("all")


text:
include("../test_interfaces.jl")
nparts = 4
prun(test_interfaces,mpi,nparts)


text:
@testset "Accuracy test for Laffine on NonhomogeneousAbsorbing" begin
    regular_grids = [0:4, -1:6]
    irregular_grids = [[-1.0; 0.0; 2.0; 5.0; 9.0], [0.2; 0.3; 0.5; 0.9; 1.6; 2.5]]
    for x̄ in [regular_grids; irregular_grids],
        (L_bc_generator, L_generator, L_affine_generator) in [(L₁₋bc, L₁₋, L₁₋affine), (L₁₊bc, L₁₊, L₁₊affine), (L₂bc, L₂, L₂affine)],
        (S_lower, S_upper) in [(0.0, 1.0), (1.0, 0.0), (1.5, -1.5), (1.2, 1.3), (1.0, 1.0)]
        # Setup
        # RHS function
        f(x) = x^2

        # grid setup
        # x̄ = 0:4
        x = interiornodes(x̄)
        M = length(x)
        Δ_0p = x̄[2] - x̄[1]
        Δ_Mp = x̄[end] - x̄[end-1]

        # boundary condition
        bc = (NonhomogeneousAbsorbing(S_lower), NonhomogeneousAbsorbing(S_upper))

        # corresponding boundary condition matrix
        B = transpose([[1; zeros(M+1)] [zeros(M+1); 1]])
        b = [S_lower; S_upper]

        # operator on the interior, with bc applied
        L_bc = L_bc_generator(x̄, bc)
        # operator on the exterior
        L = L_generator(x̄)

        # solution for interior v from L_bc
        v_bc = L_bc \ (f.(x) + L_affine_generator(x̄, f.(x), bc))
        # solution for interior v from L_bc, using Laffine
        v_bc_Laffine = L_bc \ (f.(x) + Laffine(L, f.(x), bc))
        # solution for extended v from L by stacking up boundary condition matrices
        v̄ = [L; B] \ [f.(x); b]

        # test if extrapolated solutions are identical
        @test extrapolatetoboundary(v_bc, x̄, bc) ≈ v̄
        @test extrapolatetoboundary(v_bc_Laffine, x̄, bc) ≈ v̄

    end
end

@testset "Accuracy test for Laffine on Absorbing on interior nodes" begin
    regular_grids = [0:4, -1:6]
    irregular_grids = [[-1.0; 0.0; 2.0; 5.0; 9.0], [0.2; 0.3; 0.5; 0.9; 1.6; 2.5]]
    for x̄ in [regular_grids; irregular_grids],
        (L_bc_generator, L_generator, L_affine_generator) in [(L₁₋bc, L₁₋, L₁₋affine), (L₁₊bc, L₁₊, L₁₊affine), (L₂bc, L₂, L₂affine)],
        loc in 0:3,
        S_upper in [1.0; 0.0; -1.5; 1.3; 1.0]
        # Setup
        # RHS function
        f(x) = x^2

        # grid setup
        # x̄ = 0:4
        x = interiornodes(x̄)
        M = length(x)
        Δ_0p = x̄[2] - x̄[1]
        Δ_Mp = x̄[end] - x̄[end-1]

        # boundary condition
        bc = (Absorbing(loc), NonhomogeneousAbsorbing(S_upper))

        # corresponding boundary condition matrix
        B = transpose([[1; zeros(M+1)] [zeros(M+1); 1]])
        b = [0; S_upper]

        # operator on the interior, with bc applied
        L_bc = L_bc_generator(x̄, bc)
        # operator on the exterior
        L = L_generator(x̄)
        # apply interior boundary conditions (if needed)
        RHS = f.(x)
        if (bc[1].loc > 0)
            L[1:bc[1].loc,:] = zeros(bc[1].loc,M+2)
            L[1:bc[1].loc,2:(1+bc[1].loc)] = Diagonal(ones(bc[1].loc))
            RHS[1:bc[1].loc] .= zero(eltype(L))
        end

        # solution for interior v from L_bc
        v_bc = L_bc \ (f.(x) + L_affine_generator(x̄, f.(x), bc))
        # solution for interior v from L_bc, using Laffine
        v_bc_Laffine = L_bc \ (f.(x) + Laffine(L, f.(x), bc))
        # solution for extended v from L by stacking up boundary condition matrices
        v̄ = [L; B] \ [RHS; b]

        # test if extrapolated solutions are identical
        @test extrapolatetoboundary(v_bc, x̄, bc) ≈ v̄
        @test extrapolatetoboundary(v_bc_Laffine, x̄, bc) ≈ v̄

    end
end


text:
for i in 1:length(MOLPARAM)
    if length(MOLPARAM[i].I) > 1
        bcheb = MOLPARAM[i].hascheb
        ncheb = MOLPARAM[i].ncheb
        rlerr = MOLPARAM[i].maxrelerr
        coefs = MOLPARAM[i].cheb
        @test all(rlerr .<= 0.01)
        for j = 1:length(ncheb)
            @test ncheb[j] == length(coefs[j])
            if bcheb[j]
                @test !any(isnan.(coefs[j]))
            else
                @test length(coefs[j]) == 0
            end
        end
        @test sum(MOLPARAM[i].A) <= 1.001
    end
end


text:
min_double = 5e-324
status,len,point,buffer = Grisu.fastdtoa(min_double, Grisu.SHORTEST, 0,buffer)
@test status
@test "5" == trimrep(buffer)
@test -323 == point
fill!(buffer,0);

max_double = 1.7976931348623157e308
status,len,point,buffer = Grisu.fastdtoa(max_double, Grisu.SHORTEST, 0,buffer)
@test status
@test "17976931348623157" == trimrep(buffer)
@test 309 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(4294967272.0, Grisu.SHORTEST, 0,buffer)
@test status
@test "4294967272" == trimrep(buffer)
@test 10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(4.1855804968213567e298, Grisu.SHORTEST, 0,buffer)
@test status
@test "4185580496821357" == trimrep(buffer)
@test 299 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(5.5626846462680035e-309, Grisu.SHORTEST, 0,buffer)
@test status
@test "5562684646268003" == trimrep(buffer)
@test -308 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(2147483648.0, Grisu.SHORTEST, 0,buffer)
@test status
@test "2147483648" == trimrep(buffer)
@test 10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(3.5844466002796428e+298, Grisu.SHORTEST, 0,buffer)
@test !status  # Not all Grisu.fastdtoa variants manage to compute this number.
if status
  @test "35844466002796428" == trimrep(buffer)
  @test 299 == point
  fill!(buffer,0);
end

smallest_normal64 = 0x0010000000000000
v = reinterpret(Float64,smallest_normal64)
status,len,point,buffer = Grisu.fastdtoa(v, Grisu.SHORTEST, 0, buffer)
if status
  @test "22250738585072014" == trimrep(buffer)
  @test -307 == point
  fill!(buffer,0);
end

largest_denormal64 = 0x000FFFFFFFFFFFFF
v = reinterpret(Float64,largest_denormal64)
status,len,point,buffer = Grisu.fastdtoa(v, Grisu.SHORTEST, 0, buffer)
if status
  @test "2225073858507201" == trimrep(buffer)
  @test -307 == point
  fill!(buffer,0);
end

min_float = float32(1e-45)
status,len,point,buffer = Grisu.fastdtoa(min_float, Grisu.SHORTEST, 0,buffer)
@test status
@test "1" == trimrep(buffer)
@test -44 == point
fill!(buffer,0);

max_float = 3.4028234f38 #float32(3.4028234e38)
status,len,point,buffer = Grisu.fastdtoa(max_float, Grisu.SHORTEST, 0,buffer)
@test status
@test "34028235" == trimrep(buffer)
@test 39 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(float32(4294967272.0), Grisu.SHORTEST, 0,buffer)
@test status
@test "42949673" == trimrep(buffer)
@test 10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(float32(3.32306998946228968226e+35), Grisu.SHORTEST, 0,buffer)
@test status
@test "332307" == trimrep(buffer)
@test 36 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(float32(1.2341e-41), Grisu.SHORTEST, 0,buffer)
@test status
@test "12341" == trimrep(buffer)
@test -40 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(float32(3.3554432e7), Grisu.SHORTEST, 0,buffer)
@test status
@test "33554432" == trimrep(buffer)
@test 8 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(float32(3.26494756798464e14), Grisu.SHORTEST, 0,buffer)
@test status
@test "32649476" == trimrep(buffer)
@test 15 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(float32(3.91132223637771935344e37), Grisu.SHORTEST, 0,buffer)
if status  # Not all Grisu.fastdtoa variants manage to compute this number.
  @test "39113222" == trimrep(buffer)
  @test 38 == point
  fill!(buffer,0);
end

smallest_normal32 = 0x00800000
v = reinterpret(Float32,smallest_normal32)
status,len,point,buffer = Grisu.fastdtoa(v, Grisu.SHORTEST, 0, buffer)
if status
  @test "11754944" == trimrep(buffer)
  @test -37 == point
  fill!(buffer,0);
end

largest_denormal32 = 0x007FFFFF
v = reinterpret(Float32,largest_denormal32)
status,len,point,buffer = Grisu.fastdtoa(v, Grisu.SHORTEST, 0, buffer)
@test status
@test "11754942" == trimrep(buffer)
@test -37 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(1.0, Grisu.PRECISION, 3, buffer)
@test status
@test 3 >= len-1
@test "1" == trimrep(buffer)
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(1.5, Grisu.PRECISION, 10, buffer)
if status
  @test 10 >= len-1
  @test "15" == trimrep(buffer)
  @test 1 == point
  fill!(buffer,0);
end

min_double = 5e-324
status,len,point,buffer = Grisu.fastdtoa(min_double, Grisu.PRECISION, 5,buffer)
@test status
@test "49407" == trimrep(buffer)
@test -323 == point
fill!(buffer,0);

max_double = 1.7976931348623157e308
status,len,point,buffer = Grisu.fastdtoa(max_double, Grisu.PRECISION, 7,buffer)
@test status
@test "1797693" == trimrep(buffer)
@test 309 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(4294967272.0, Grisu.PRECISION, 14,buffer)
if status
  @test 14 >= len-1
  @test "4294967272" == trimrep(buffer)
  @test 10 == point
  fill!(buffer,0);
end

status,len,point,buffer = Grisu.fastdtoa(4.1855804968213567e298, Grisu.PRECISION, 17,buffer)
@test status
@test "41855804968213567" == trimrep(buffer)
@test 299 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(5.5626846462680035e-309, Grisu.PRECISION, 1,buffer)
@test status
@test "6" == trimrep(buffer)
@test -308 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(2147483648.0, Grisu.PRECISION, 5,buffer)
@test status
@test "21475" == trimrep(buffer)
@test 10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.fastdtoa(3.5844466002796428e+298, Grisu.PRECISION, 10,buffer)
@test status
@test 10 >= len-1
@test "35844466" == trimrep(buffer)
@test 299 == point
fill!(buffer,0);

smallest_normal64 = 0x0010000000000000
v = reinterpret(Float64,smallest_normal64)
status,len,point,buffer = Grisu.fastdtoa(v, Grisu.PRECISION, 17, buffer)
@test status
@test "22250738585072014" == trimrep(buffer)
@test -307 == point
fill!(buffer,0);

largest_denormal64 = 0x000FFFFFFFFFFFFF
v = reinterpret(Float64,largest_denormal64)
status,len,point,buffer = Grisu.fastdtoa(v, Grisu.PRECISION, 17, buffer)
@test status
@test 20 >= len-1
@test "22250738585072009" == trimrep(buffer)
@test -307 == point
fill!(buffer,0);

v = 3.3161339052167390562200598e-237
status,len,point,buffer = Grisu.fastdtoa(v, Grisu.PRECISION, 18, buffer)
@test status
@test "331613390521673906" == trimrep(buffer)
@test -236 == point
fill!(buffer,0);

v = 7.9885183916008099497815232e+191
status,len,point,buffer = Grisu.fastdtoa(v, Grisu.PRECISION, 4, buffer)
@test status
@test "7989" == trimrep(buffer)
@test 192 == point
fill!(buffer,0);


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

"""
    PointSet(coords)

A set of points with coordinates `coords`. Each point is represented by
a static vector or tuple. Alternatively, `coords` can be a matrix where
the number of rows equals the number of dimensions.

## Examples

Create a 2D point set with 100 points:

```julia
julia> PointSet(rand(2,100))
```

Or equivalently, using a vector of tuples:

```julia
julia> PointSet([(rand(),rand()) for i in 1:100])
```
"""
struct PointSet{T,N} <: SpatialDomain{T,N}
  coords::Vector{SVector{N,T}}
end

PointSet(coords::AbstractVector{<:NTuple}) =
  PointSet(SVector.(coords))

PointSet(coords::AbstractMatrix) =
  PointSet(SVector{size(coords,1)}.(eachcol(coords)))

nelms(ps::PointSet) = length(ps.coords)

coordinates!(buff, ps::PointSet, ind::Int) =
  @inbounds buff .= ps.coords[ind]

==(ps₁::PointSet, ps₂::PointSet) = ps₁.coords == ps₂.coords

# ------------
# IO methods
# ------------
function Base.show(io::IO, ps::PointSet{T,N}) where {N,T}
  npts = length(ps.coords)
  print(io, "$npts PointSet{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", ps::PointSet{T,N}) where {N,T}
  println(io, ps)
  Base.print_array(io, reduce(hcat, ps.coords))
end


text:
# Avoid ambiguity
function update_batch_vertex(g::GATConv, M::AbstractMatrix, X::CuMatrix)
    M = convert(typeof(X), M)
    update_batch_vertex(g, M, X)
end

function update_batch_vertex(g::GATConv, M::CuMatrix, X::CuMatrix)
    g.concat || (M = mean(M, dims=2))
    return M .+ g.bias
end


text:
using Query
using DataFrames

df = DataFrame(name=["John", "Sally", "Kirk"], age=[23., 42., 59.], children=[3,2,2])

x = @from i in df begin
    @group i by i.children into g
    @select {Key=key(g),Count=length(g)}
    @collect DataFrame
end

println(x)


text:
using Random
using LinearAlgebra
using Statistics
using Optim
using DataFrames
using CSV
using HTTP
using GLM
using FreqTables
using SMM
using DataFramesMeta
using Pipe,Distributions

function wrapper()
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS1-julia-intro/nlsw88.csv"
df = CSV.read(HTTP.get(url).body, DataFrame)
X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.married.==1

# 1)
bols = inv(X'*X)*X'*y
println(bols)

function ols_gmm(beta, X, y)
    g = y .- X*beta
    # I defined in LinearAlgebra
    J = g'*I*g
    return J
end
α̂_optim = optimize(a -> ols_gmm(a, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-8, iterations=100_000))
println(α̂_optim.minimizer)

# the answers are very close
# α̂_optim.minimizer.-bols |> println

#2)
#a) changed alpha_start to alpha_true to get the result quickly
include("ps2question5.jl")

#b)
alpha_true = [.1910213,-.0335262,.5963968,.4165052,-.1698368,-.0359784,1.30684,-.430997,.6894727,-.0104578,.5231634,-1.492475,
            -2.26748,-.0053001,1.391402,-.9849661,-1.398468,-.0142969,-.0176531,-1.495123,.2454891,-.0067267,-.5382892,-3.78975]

function gmm_logit(alpha, X, y)

    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==j
    end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

    num = zeros(N,J)
    dem = zeros(N)
    for j=1:J
        num[:,j] = exp.(X*bigAlpha[:,j])
        dem .+= num[:,j]
    end

    P = num./repeat(dem,1,J)

    loglike = -sum( bigY.*log.(P) )
    g = @pipe (bigY .- P) |> reshape(_,(J*N, 1))
    J = g'*I*g
    return J[1,1]
end

alpha_hat_optim = optimize(a -> gmm_logit(a, X, y), alpha_true, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
alpha_hat_mle = alpha_hat_optim.minimizer
println(alpha_hat_mle)

#c)
alpha_rand = rand(6*size(X,2))

alpha_hat_optim = optimize(a -> gmm_logit(a, X, y), alpha_rand, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
alpha_hat_mle = alpha_hat_optim.minimizer
println(alpha_hat_mle)

# I got different etimates which means that  the function isn't globally concave

#4)

N=size(X,1)
K=size(X,2)
J = length(unique(y))
Xrand = rand(N,K)
beta = rand(J)
# f)
#https://github.com/JuliaStats/Distributions.jl/blob/master/src/univariate/continuous/gumbel.jl
epsilon = rand(Gumbel(),N,J)
YY = zeros(N)
for i in 1:N
YY[i] = argmax(X[i]*beta[:] + epsilon[i,:])
end

# Q4

MA = SMM.parallelNormal() # Note: this line may take up to 5 minutes to execute
dc = SMM.history(MA.chains[1])
dc = dc[dc[:accepted].==true, :]
println(describe(dc))

# Q5
function ols_smm(θ, X, y, D)
	N=size(X,1)
	K=size(X,2)
	J = length(unique(y))
    β = θ[1:end-1]
    σ = θ[end]
    if length(β)==1
        β = β[1]
    end
    bigY = zeros(N,J)
	for j=1:J
            bigY[:,j] = y.==j
    end
	bigY2 = zeros(N,J)

    # N+1 moments in both model and data
    gmodel = zeros(N+1,D)
    # data moments are just the y vector itself
    # and the variance of the y vector
    gdata  = vcat(y,var(y))
    #### !!!!!!!!!!!!!!!!!!!!!!!!!!!!! ####
    # This is critical!                   #
    Random.seed!(1234)                    #
    # You must always use the same ε draw #
    # for every guess of θ!               #
    #### !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ###
    # simulated model moments
    for d=1:D
        ε = σ*randn(N)
        ỹ = X*β .+ ε
        gmodel[1:end-1,d] = ỹ
        gmodel[  end  ,d] = var(ỹ)
    end
    # criterion function
    err = vec(gdata .- mean(gmodel; dims=2))
    # weighting matrix is the identity matrix
    # minimize weighted difference between data and moments
    J = err'*I*err
    return J
end
return nothing
end
wrapper()


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

export
    avutil_version,
    avutil_configuration,
    avutil_license,
    av_get_media_type_string,
    av_get_picture_type_char,
    av_x_if_null,
    av_int_list_length_for_size,
    av_fopen_utf8,
    av_get_time_base_q

function avutil_version()
    ccall((:avutil_version,libavutil),UInt32,())
end

function avutil_configuration()
    ccall((:avutil_configuration,libavutil),Ptr{UInt8},())
end

function avutil_license()
    ccall((:avutil_license,libavutil),Ptr{UInt8},())
end

function av_get_media_type_string(media_type::AVMediaType)
    ccall((:av_get_media_type_string,libavutil),Ptr{UInt8},(AVMediaType,),media_type)
end

function av_get_picture_type_char(pict_type::AVPictureType)
    ccall((:av_get_picture_type_char,libavutil),UInt8,(AVPictureType,),pict_type)
end

function av_x_if_null(p,x)
    ccall((:av_x_if_null,libavutil),Ptr{Void},(Ptr{Void},Ptr{Void}),p,x)
end

function av_int_list_length_for_size(elsize::Integer,list,term::UInt64)
    ccall((:av_int_list_length_for_size,libavutil),UInt32,(UInt32,Ptr{Void},UInt64),elsize,list,term)
end

function av_fopen_utf8(path,mode)
    ccall((:av_fopen_utf8,libavutil),Ptr{FILE},(Ptr{UInt8},Ptr{UInt8}),path,mode)
end

function av_get_time_base_q()
    ccall((:av_get_time_base_q,libavutil),AVRational,())
end


text:
import Tensors
using jFEMTools
import WriteVTK

nel = (2,2)
LL = Tensors.Vec{2}((0.0,0.0))
UR = Tensors.Vec{2}((1.0,3.0))
mesh = rectangle_mesh(HexagonCell,(1,3),LL,UR);
vtk_file = vtk_grid("hexagon_mesh", mesh)
outfiles = WriteVTK.vtk_save(vtk_file)

#Plot mesh
using Makie
import AbstractPlotting
include("src/plot_recipes.jl")
#popdisplay(AbstractPlotting.PlotDisplay())
#AbstractPlotting.inline!(true)
scene = Scene(resolution = (400, 200))
plot!(scene, mesh)

nel = (2,2,2)
mesh3D_1 = hyper_rectagle_mesh2(HexahedronCell,nel)
vtk_file = vtk_grid("hexahedron_mesh", mesh3D_1, compress=false)
outfiles = WriteVTK.vtk_save(vtk_file)

mesh3D_2 = hyper_rectagle_mesh2(TetrahedronCell,nel)
vtk_file = vtk_grid("tetrahedron_mesh", mesh3D_2, compress=false)
outfiles = WriteVTK.vtk_save(vtk_file)

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
function [lam,r]=cart2Radial(x,y)
r=sqrt(x^2+y^2);
lam=atan2(y,x);
if lam<0.0
  lam=lam+2*pi;
end
end