text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

using Core.Intrinsics: llvmcall

import Base: setindex!, getindex, unsafe_convert
import Base.Sys: ARCH, WORD_SIZE

export
    Atomic,
    atomic_cas!,
    atomic_xchg!,
    atomic_add!, atomic_sub!,
    atomic_and!, atomic_nand!, atomic_or!, atomic_xor!,
    atomic_max!, atomic_min!,
    atomic_fence

# Disable 128-bit types on 32-bit Intel sytems due to LLVM problems;
# see <https://github.com/JuliaLang/julia/issues/14818> (fixed on LLVM 3.9)
# 128-bit atomics do not exist on AArch32.
if (VersionNumber(Base.libllvm_version) < v"3.9-" && ARCH === :i686) ||
        startswith(string(ARCH), "arm")
    const inttypes = (Int8, Int16, Int32, Int64,
                      UInt8, UInt16, UInt32, UInt64)
else
    const inttypes = (Int8, Int16, Int32, Int64, Int128,
                      UInt8, UInt16, UInt32, UInt64, UInt128)
end
const floattypes = (Float16, Float32, Float64)
# TODO: Support Bool, Ptr
const atomictypes = (inttypes..., floattypes...)
typealias IntTypes Union{inttypes...}
typealias FloatTypes Union{floattypes...}
typealias AtomicTypes Union{atomictypes...}

"""
    Threads.Atomic{T}

Holds a reference to an object of type `T`, ensuring that it is only
accessed atomically, i.e. in a thread-safe manner.

Only certain "simple" types can be used atomically, namely the
bitstypes integer and float-point types. These are `Int8`...`Int128`,
`UInt8`...`UInt128`, and `Float16`...`Float64`.

New atomic objects can be created from a non-atomic values; if none is
specified, the atomic object is initialized with zero.

Atomic objects can be accessed using the `[]` notation:

```jldoctest
julia> x = Threads.Atomic{Int}(3)
Base.Threads.Atomic{Int64}(3)

julia> x[] = 1
1

julia> x[]
1
```

Atomic operations use an `atomic_` prefix, such as `atomic_add!`,
`atomic_xchg!`, etc.
"""
type Atomic{T<:AtomicTypes}
    value::T
    Atomic{T}() where T<:AtomicTypes = new(zero(T))
    Atomic{T}(value) where T<:AtomicTypes = new(value)
end

Atomic() = Atomic{Int}()

"""
    Threads.atomic_cas!{T}(x::Atomic{T}, cmp::T, newval::T)

Atomically compare-and-set `x`

Atomically compares the value in `x` with `cmp`. If equal, write
`newval` to `x`. Otherwise, leaves `x` unmodified. Returns the old
value in `x`. By comparing the returned value to `cmp` (via `===`) one
knows whether `x` was modified and now holds the new value `newval`.

For further details, see LLVM's `cmpxchg` instruction.

This function can be used to implement transactional semantics. Before
the transaction, one records the value in `x`. After the transaction,
the new value is stored only if `x` has not been modified in the mean
time.

```jldoctest
julia> x = Threads.Atomic{Int}(3)
Base.Threads.Atomic{Int64}(3)

julia> Threads.atomic_cas!(x, 4, 2);

julia> x
Base.Threads.Atomic{Int64}(3)

julia> Threads.atomic_cas!(x, 3, 2);

julia> x
Base.Threads.Atomic{Int64}(2)
```
"""
function atomic_cas! end

"""
    Threads.atomic_xchg!{T}(x::Atomic{T}, newval::T)

Atomically exchange the value in `x`

Atomically exchanges the value in `x` with `newval`. Returns the **old**
value.

For further details, see LLVM's `atomicrmw xchg` instruction.

```jldoctest
julia> x = Threads.Atomic{Int}(3)
Base.Threads.Atomic{Int64}(3)

julia> Threads.atomic_xchg!(x, 2)
3

julia> x[]
2
```
"""
function atomic_xchg! end

"""
    Threads.atomic_add!{T}(x::Atomic{T}, val::T)

Atomically add `val` to `x`

Performs `x[] += val` atomically. Returns the **old** value.

For further details, see LLVM's `atomicrmw add` instruction.

```jldoctest
julia> x = Threads.Atomic{Int}(3)
Base.Threads.Atomic{Int64}(3)

julia> Threads.atomic_add!(x, 2)
3

julia> x[]
5
```
"""
function atomic_add! end

"""
    Threads.atomic_sub!{T}(x::Atomic{T}, val::T)

Atomically subtract `val` from `x`

Performs `x[] -= val` atomically. Returns the **old** value.

For further details, see LLVM's `atomicrmw sub` instruction.

```jldoctest
julia> x = Threads.Atomic{Int}(3)
Base.Threads.Atomic{Int64}(3)

julia> Threads.atomic_sub!(x, 2)
3

julia> x[]
1
```
"""
function atomic_sub! end

"""
    Threads.atomic_and!{T}(x::Atomic{T}, val::T)

Atomically bitwise-and `x` with `val`

Performs `x[] &= val` atomically. Returns the **old** value.

For further details, see LLVM's `atomicrmw and` instruction.

```jldoctest
julia> x = Threads.Atomic{Int}(3)
Base.Threads.Atomic{Int64}(3)

julia> Threads.atomic_and!(x, 2)
3

julia> x[]
2
```
"""
function atomic_and! end

"""
    Threads.atomic_nand!{T}(x::Atomic{T}, val::T)

Atomically bitwise-nand (not-and) `x` with `val`

Performs `x[] = ~(x[] & val)` atomically. Returns the **old** value.

For further details, see LLVM's `atomicrmw nand` instruction.

```jldoctest
julia> x = Threads.Atomic{Int}(3)
Base.Threads.Atomic{Int64}(3)

julia> Threads.atomic_nand!(x, 2)
3

julia> x[]
-3
```
"""
function atomic_nand! end

"""
    Threads.atomic_or!{T}(x::Atomic{T}, val::T)

Atomically bitwise-or `x` with `val`

Performs `x[] |= val` atomically. Returns the **old** value.

For further details, see LLVM's `atomicrmw or` instruction.

```jldoctest
julia> x = Threads.Atomic{Int}(5)
Base.Threads.Atomic{Int64}(5)

julia> Threads.atomic_or!(x, 7)
5

julia> x[]
7
```
"""
function atomic_or! end

"""
    Threads.atomic_xor!{T}(x::Atomic{T}, val::T)

Atomically bitwise-xor (exclusive-or) `x` with `val`

Performs `x[] \$= val` atomically. Returns the **old** value.

For further details, see LLVM's `atomicrmw xor` instruction.

```jldoctest
julia> x = Threads.Atomic{Int}(5)
Base.Threads.Atomic{Int64}(5)

julia> Threads.atomic_xor!(x, 7)
5

julia> x[]
2
```
"""
function atomic_xor! end

"""
    Threads.atomic_max!{T}(x::Atomic{T}, val::T)

Atomically store the maximum of `x` and `val` in `x`

Performs `x[] = max(x[], val)` atomically. Returns the **old** value.

For further details, see LLVM's `atomicrmw max` instruction.

```jldoctest
julia> x = Threads.Atomic{Int}(5)
Base.Threads.Atomic{Int64}(5)

julia> Threads.atomic_max!(x, 7)
5

julia> x[]
7
```
"""
function atomic_max! end

"""
    Threads.atomic_min!{T}(x::Atomic{T}, val::T)

Atomically store the minimum of `x` and `val` in `x`

Performs `x[] = min(x[], val)` atomically. Returns the **old** value.

For further details, see LLVM's `atomicrmw min` instruction.

```jldoctest
julia> x = Threads.Atomic{Int}(7)
Base.Threads.Atomic{Int64}(7)

julia> Threads.atomic_min!(x, 5)
7

julia> x[]
5
```
"""
function atomic_min! end

unsafe_convert{T}(::Type{Ptr{T}}, x::Atomic{T}) = convert(Ptr{T}, pointer_from_objref(x))
setindex!{T}(x::Atomic{T}, v) = setindex!(x, convert(T, v))

const llvmtypes = Dict(
    Bool => "i1",
    Int8 => "i8", UInt8 => "i8",
    Int16 => "i16", UInt16 => "i16",
    Int32 => "i32", UInt32 => "i32",
    Int64 => "i64", UInt64 => "i64",
    Int128 => "i128", UInt128 => "i128",
    Float16 => "i16", # half
    Float32 => "float",
    Float64 => "double",
)
inttype{T<:Integer}(::Type{T}) = T
inttype(::Type{Float16}) = Int16
inttype(::Type{Float32}) = Int32
inttype(::Type{Float64}) = Int64

alignment{T}(::Type{T}) = ccall(:jl_alignment, Cint, (Csize_t,), sizeof(T))

# All atomic operations have acquire and/or release semantics, depending on
# whether the load or store values. Most of the time, this is what one wants
# anyway, and it's only moderately expensive on most hardware.
for typ in atomictypes
    lt = llvmtypes[typ]
    ilt = llvmtypes[inttype(typ)]
    rt = VersionNumber(Base.libllvm_version) >= v"3.6" ? "$lt, $lt*" : "$lt*"
    irt = VersionNumber(Base.libllvm_version) >= v"3.6" ? "$ilt, $ilt*" : "$ilt*"
    if VersionNumber(Base.libllvm_version) >= v"3.8"
        @eval getindex(x::Atomic{$typ}) =
            llvmcall($"""
                     %rv = load atomic $rt %0 acquire, align $(alignment(typ))
                     ret $lt %rv
                     """, $typ, Tuple{Ptr{$typ}}, unsafe_convert(Ptr{$typ}, x))
        @eval setindex!(x::Atomic{$typ}, v::$typ) =
            llvmcall($"""
                     store atomic $lt %1, $lt* %0 release, align $(alignment(typ))
                     ret void
                     """, Void, Tuple{Ptr{$typ},$typ}, unsafe_convert(Ptr{$typ}, x), v)
    else
        if typ <: Integer
            @eval getindex(x::Atomic{$typ}) =
                llvmcall($"""
                         %rv = load atomic $rt %0 acquire, align $(alignment(typ))
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}}, unsafe_convert(Ptr{$typ}, x))
            @eval setindex!(x::Atomic{$typ}, v::$typ) =
                llvmcall($"""
                         store atomic $lt %1, $lt* %0 release, align $(alignment(typ))
                         ret void
                         """, Void, Tuple{Ptr{$typ},$typ}, unsafe_convert(Ptr{$typ}, x), v)
        else
            @eval getindex(x::Atomic{$typ}) =
                llvmcall($"""
                         %iptr = bitcast $lt* %0 to $ilt*
                         %irv = load atomic $irt %iptr acquire, align $(alignment(typ))
                         %rv = bitcast $ilt %irv to $lt
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}}, unsafe_convert(Ptr{$typ}, x))
            @eval setindex!(x::Atomic{$typ}, v::$typ) =
                llvmcall($"""
                         %iptr = bitcast $lt* %0 to $ilt*
                         %ival = bitcast $lt %1 to $ilt
                         store atomic $ilt %ival, $ilt* %iptr release, align $(alignment(typ))
                         ret void
                         """, Void, Tuple{Ptr{$typ},$typ}, unsafe_convert(Ptr{$typ}, x), v)
        end
    end
    # Note: atomic_cas! succeeded (i.e. it stored "new") if and only if the result is "cmp"
    if VersionNumber(Base.libllvm_version) >= v"3.5"
        if typ <: Integer
            @eval atomic_cas!(x::Atomic{$typ}, cmp::$typ, new::$typ) =
                llvmcall($"""
                         %rs = cmpxchg $lt* %0, $lt %1, $lt %2 acq_rel acquire
                         %rv = extractvalue { $lt, i1 } %rs, 0
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                         unsafe_convert(Ptr{$typ}, x), cmp, new)
        else
            @eval atomic_cas!(x::Atomic{$typ}, cmp::$typ, new::$typ) =
                llvmcall($"""
                         %iptr = bitcast $lt* %0 to $ilt*
                         %icmp = bitcast $lt %1 to $ilt
                         %inew = bitcast $lt %2 to $ilt
                         %irs = cmpxchg $ilt* %iptr, $ilt %icmp, $ilt %inew acq_rel acquire
                         %irv = extractvalue { $ilt, i1 } %irs, 0
                         %rv = bitcast $ilt %irv to $lt
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                         unsafe_convert(Ptr{$typ}, x), cmp, new)
        end
    else
        if typ <: Integer
            @eval atomic_cas!(x::Atomic{$typ}, cmp::$typ, new::$typ) =
                llvmcall($"""
                         %rv = cmpxchg $lt* %0, $lt %1, $lt %2 acq_rel
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                         unsafe_convert(Ptr{$typ}, x), cmp, new)
        else
            @eval atomic_cas!(x::Atomic{$typ}, cmp::$typ, new::$typ) =
                llvmcall($"""
                         %iptr = bitcast $lt* %0 to $ilt*
                         %icmp = bitcast $lt %1 to $ilt
                         %inew = bitcast $lt %2 to $ilt
                         %irv = cmpxchg $ilt* %iptr, $ilt %icmp, $ilt %inew acq_rel
                         %rv = bitcast $ilt %irv to $lt
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                         unsafe_convert(Ptr{$typ}, x), cmp, new)
        end
    end
    for rmwop in [:xchg, :add, :sub, :and, :nand, :or, :xor, :max, :min]
        rmw = string(rmwop)
        fn = Symbol("atomic_", rmw, "!")
        if (rmw == "max" || rmw == "min") && typ <: Unsigned
            # LLVM distinguishes signedness in the operation, not the integer type.
            rmw = "u" * rmw
        end
        if typ <: Integer
            @eval $fn(x::Atomic{$typ}, v::$typ) =
                llvmcall($"""
                         %rv = atomicrmw $rmw $lt* %0, $lt %1 acq_rel
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}, $typ}, unsafe_convert(Ptr{$typ}, x), v)
        else
            rmwop == :xchg || continue
            @eval $fn(x::Atomic{$typ}, v::$typ) =
                llvmcall($"""
                         %iptr = bitcast $lt* %0 to $ilt*
                         %ival = bitcast $lt %1 to $ilt
                         %irv = atomicrmw $rmw $ilt* %iptr, $ilt %ival acq_rel
                         %rv = bitcast $ilt %irv to $lt
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}, $typ}, unsafe_convert(Ptr{$typ}, x), v)
        end
    end
end

# Provide atomic floating-point operations via atomic_cas!
const opnames = Dict{Symbol, Symbol}(:+ => :add, :- => :sub)
for op in [:+, :-, :max, :min]
    opname = get(opnames, op, op)
    @eval function $(Symbol("atomic_", opname, "!")){T<:FloatTypes}(var::Atomic{T}, val::T)
        IT = inttype(T)
        old = var[]
        while true
            new = $op(old, val)
            cmp = old
            old = atomic_cas!(var, cmp, new)
            reinterpret(IT, old) == reinterpret(IT, cmp) && return new
            # Temporary solution before we have gc transition support in codegen.
            ccall(:jl_gc_safepoint, Void, ())
        end
    end
end

"""
    Threads.atomic_fence()

Insert a sequential-consistency memory fence

Inserts a memory fence with sequentially-consistent ordering
semantics. There are algorithms where this is needed, i.e. where an
acquire/release ordering is insufficient.

This is likely a very expensive operation. Given that all other atomic
operations in Julia already have acquire/release semantics, explicit
fences should not be necessary in most cases.

For further details, see LLVM's `fence` instruction.
"""
atomic_fence() = llvmcall("""
                          fence seq_cst
                          ret void
                          """, Void, Tuple{})


text:
# LAxis and LScene currently have to be AbstractScenes because they need to work with
# the plotting pipeline, so they are not LObjects and manually get added here
const Layoutable = Union{LAxis, LScene, LObject}

# almost like in AbstractPlotting
# make fields type inferrable
# just access attributes directly instead of via indexing detour

@generated Base.hasfield(x::T, ::Val{key}) where {T<:Layoutable, key} = :($(key in fieldnames(T)))

@inline function Base.getproperty(x::T, key::Symbol) where T <: Layoutable
    if hasfield(x, Val(key))
        getfield(x, key)
    else
        x.attributes[key]
    end
end

@inline function Base.setproperty!(x::T, key::Symbol, value) where T <: Layoutable
    if hasfield(x, Val(key))
        setfield!(x, key, value)
    else
        x.attributes[key][] = value
    end
end

# propertynames should list fields and attributes
function Base.propertynames(layoutable::T) where T <: Layoutable
    [fieldnames(T)..., keys(layoutable.attributes)...]
end

# treat all layoutables as scalars when broadcasting
Base.Broadcast.broadcastable(l::Layoutable) = Ref(l)

function Base.show(io::IO, ::T) where T <: LObject
    print(io, "$T()")
end

ext: jl


text:
@testset "select" begin
    m = GBMatrix([[1,2,3] [4,5,6] [7,8,9]])
    s = select(tril, m)
    @test s[1,2] === nothing && s[3,1] == 3
    s = select(<, m, 6)
    @test s[2,2] == 5 && s[3,3] === nothing
end


text:
# Trying out derivatives of cel_bulirsch with respect to
# k_c, p, a, b:
include("cel_bulirsch.jl")

using ForwardDiff
using DiffRules

dcel_dn = m -> ForwardDiff.derivative(cel_bulirsch,m)

function cel_bad_grad(r::T,b::T) where {T <: Real}
x=[r,b]
function diff_cel(x::Array{T,1}) where {T <: Real}
r,b=x
onembpr2 = (1-b-r)*(1+b+r); onembmr2=(r-b+1)*(b-r+1); fourbr = 4b*r
k2 = onembpr2/fourbr+1
k2c = -onembpr2/fourbr; kc = sqrt(b+r-1)*sqrt(b+r+1)/sqrt(4*b*r)
#p=(b-r)^2*k2c; a=0.0; b_cel=3*k2c*(b-r)*(b+r)
p=(b-r)^2*k2c; a=0.0; b_cel=3*(b+r)*(b-r)*k2c
return cel_bulirsch(k2,kc,p,zero(b),b_cel)
end
out = DiffResults.GradientResult(x)
# Compute the Jacobian (and value):
out = ForwardDiff.gradient!(out,diff_cel,x)
# Place the value in the cel_b vector:
cel_b = DiffResults.value(out)
# And, place the Jacobian in an array:
cel_gradient= DiffResults.gradient(out)
return cel_b,cel_gradient
end

function cel_grad(kc::T,p::T,a::T,b::T) where {T <: Real}
  # Computes the derivative of cel(k_c,p,a,b) with respect to (k_c,p,a,b)
  # Create a vector for use with ForwardDiff
  x=[kc,p,a,b]
  # Now, define a wrapper of cel_bulirsch for use with ForwardDiff:
  function diff_cel(x::Array{T,1}) where {T <: Real}
  # x should be a four-element vector with values [kc,p,a,b]
  kc,p,a,b = x
  return cel_bulirsch(1.0-kc*kc,kc,p,a,b)
  end

  # Set up a type to store s_n and it's Jacobian with respect to x:
  out = DiffResults.GradientResult(x)
  # Compute the Jacobian (and value):
  out = ForwardDiff.gradient!(out,diff_cel,x)
  # Place the value in the cel_b vector:
  cel_b = DiffResults.value(out)
  # And, place the Jacobian in an array:
  cel_gradient= DiffResults.gradient(out)
return cel_b,cel_gradient
end

function cel_grad_anal(kc::T,p::T,a::T,b::T) where {T <: Real}
# Now, compute the gradients analytically:
kc2 = kc*kc
dcel_dkc = -kc/(p-kc2)*(cel_bulirsch(1-kc*kc,kc,kc2,a,b)-cel_b)
#dcel_dp = (2*p*b-p*a-b)/(2*p*(1-p))*cel_bulirsch(1-kc*kc,kc,p,zero(p),one(p))+
#  (b-a*p)/(2*p*(1-p)*(p-kc2))*cel_bulirsch(1-kc*kc,kc,one(p),one(p),kc2) -
#  (b-a*p)/(2*(1-p)*(p-kc2))*cel_bulirsch(1-kc*kc,kc,p,one(p),one(p))
dcel_dp = ((kc2*(b+a*p-2*b*p)+p*(3*b*p-a*p^2-2*b))*cel_bulirsch(1-kc*kc,kc,p,zero(p),one(p))+
  (b-a*p)*cel_bulirsch(1-kc*kc,kc,one(p),1-p,kc2-p))/(2*p*(1-p)*(p-kc2))
dcel_da = cel_bulirsch(1-kc*kc,kc,p,one(p),zero(p))
dcel_db = cel_bulirsch(1-kc*kc,kc,p,zero(p),one(p))
return [dcel_dkc,dcel_dp,dcel_da,dcel_db]
end

function cel_grad_num(kc::T,p::T,a::T,b::T) where {T <: Real}
# Now, compute the gradients numerically:
dq = 1e-24
dcel_dkc = convert(Float64,(cel_bulirsch(big(1)-(big(kc)+big(dq))^2,big(kc)+big(dq),big(p),big(a),big(b))-
  cel_bulirsch(big(1)-(big(kc)-big(dq))^2,big(kc)-big(dq),big(p),big(a),big(b)))/(2*big(dq)))
dcel_dp = convert(Float64,(cel_bulirsch(big(1)-big(kc)^2,big(kc),big(p)+big(dq),big(a),big(b))-
  cel_bulirsch(big(1)-big(kc)^2,big(kc),big(p)-big(dq),big(a),big(b)))/(2*big(dq)))
dcel_da = convert(Float64,(cel_bulirsch(big(1)-big(kc)^2,big(kc),big(p),big(a)+big(dq),big(b))-
  cel_bulirsch(big(1)-big(kc)^2,big(kc),big(p),big(a)-big(dq),big(b)))/(2*big(dq)))
dcel_db = convert(Float64,(cel_bulirsch(big(1)-big(kc)^2,big(kc),big(p),big(a),big(b)+big(dq))-
  cel_bulirsch(big(1)-big(kc)^2,big(kc),big(p),big(a),big(b)-big(dq)))/(2*big(dq)))
return [dcel_dkc,dcel_dp,dcel_da,dcel_db]
end

r = 100.0; b = 100.0-1e-12
onembpr2 = (1-b-r)*(1+b+r); onembmr2=(r-b+1)*(b-r+1); fourbr = 4b*r
k2 = onembpr2/fourbr+1
k2c = -onembpr2/fourbr; kc = sqrt(k2c)
p=(b-r)^2*k2c; a=0.0; b_cel=3*k2c*(b-r)*(b+r)

#kc = rand(); p = rand(); a=rand(); b= rand()
cel_b,cel_gradient= cel_grad(kc,p,a,b_cel)
dcel_anal = cel_grad_anal(kc,p,a,b_cel)
dcel_num  = cel_grad_num(kc,p,a,b_cel)
println("d cel/d kc: ",dcel_anal[1]," ",cel_gradient[1]," ",dcel_num[1]," diff 1-2: ",dcel_anal[1]-cel_gradient[1]," ",dcel_num[1]-cel_gradient[1])
println("d cel/d p: ",dcel_anal[2]," ",cel_gradient[2]," ",dcel_num[2]," diff 1-2: ",dcel_anal[2]-cel_gradient[2]," ",dcel_num[2]-cel_gradient[2])
println("d cel/d a: ",dcel_anal[3]," ",cel_gradient[3]," ",dcel_num[3]," diff 1-2: ",dcel_anal[3]-cel_gradient[3]," ",dcel_num[3]-cel_gradient[3])
println("d cel/d b: ",dcel_anal[4]," ",cel_gradient[4]," ",dcel_num[4]," diff 1-2: ",dcel_anal[4]-cel_gradient[4]," ",dcel_num[4]-cel_gradient[4])

function cel_bad_anal(r::T,b::T) where {T <: Real}
onembpr2 = (1-b-r)*(1+b+r); onembmr2=(r-b+1)*(b-r+1); fourbr = 4b*r
k2 = onembpr2/fourbr+1
k2c = -onembpr2/fourbr; kc = sqrt(k2c)
p=(b-r)^2*k2c; a=zero(r); b_cel=3*k2c*(b-r)*(b+r)
dcel_anal = cel_grad_anal(kc,p,a,b_cel)
dceldr = (dcel_anal[1]*(1-(b-r)*(b+r))/(2*kc)+dcel_anal[2]*(1-3r^2-b^2)*(b-r)*(b+r)-dcel_anal[4]*3*(3*r^4+b^4-b^2-r^2+4*b*r^3))/(4*b*r^2)
dceldb = (dcel_anal[1]*(1+(b-r)*(b+r))/(2*kc)+dcel_anal[2]*(r^2+3b^2-1)*(b-r)*(b+r)+dcel_anal[4]*3*(3*b^4+r^4-b^2-r^2+4*b^3*r))/(4*b^2*r)
return [dceldr,dceldb]
end

dcel_bad_anal = cel_bad_anal(big(r),big(b))
println("dcel/dr: ",dcel_bad_anal[1]," dcel/db: ",dcel_bad_anal[2])
dcel_bad_anal = cel_bad_anal(r,b)
println("dcel/dr: ",dcel_bad_anal[1]," dcel/db: ",dcel_bad_anal[2])
cel_bad,cel_bad_gradient=cel_bad_grad(big(r),big(b))
println("dcel/dr: ",cel_bad_gradient[1]," dcel/db: ",cel_bad_gradient[2])
cel_bad,cel_bad_gradient=cel_bad_grad(r,b)
println("dcel/dr: ",cel_bad_gradient[1]," dcel/db: ",cel_bad_gradient[2])

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
using StructTypes, JSON3

sanitize(t::Tuple) = map(sanitize, t)
sanitize(t::NamedTuple) = map(sanitize, t)
sanitize(x::Function) = repr(x)
sanitize(d::Dict) = Dict([sanitize(k)=>sanitize(d[k]) for k in keys(d)])
sanitize(a::Array) = sanitize.(a)
sanitize(x) = x

StructTypes.StructType(::Type{<:Dagger.Processor}) = StructTypes.CustomStruct()
StructTypes.lower(proc::Dagger.Processor) = repr(proc)
StructTypes.StructType(::Type{<:Dagger.Chunk}) = StructTypes.CustomStruct()
StructTypes.lower(chunk::Dagger.Chunk) = repr(chunk)
StructTypes.StructType(::Type{<:Dagger.Event}) = StructTypes.CustomStruct()
function StructTypes.lower(ev::Dagger.Event{T}) where T
    kind = T
    category = ev.category
    timestamp = ev.timestamp
    id = ev.id
    timeline = if kind == :move
        sanitize((;f=ev.timeline.f,id=ev.timeline.id))
    elseif kind == :evict
        (;)
    else
        sanitize(ev.timeline)
    end
    return (;kind, category, timestamp, id, timeline)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GoogleSparseHash_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GoogleSparseHash")
JLLWrappers.@generate_main_file("GoogleSparseHash", UUID("6345bf58-7ab3-5e61-945d-7ed95101af62"))
end  # module GoogleSparseHash_jll


text:
using CSV, DataFrames, Pipe, Dates

df = DataFrame(CSV.File("data/zhvi.csv"))
select!(df, :RegionName, 294:297)
println(df)
# Select by column
select(df, :RegionName, 4:5)

# Filter (note: the data frame is in the 2nd argument)
filter(:RegionName => ==("Abilene, TX"), df)
filter("2020-01-31" => >(400_000), df)

# Sort
sort(df, :RegionName)
sort(df, :RegionName, rev = true)

# Transform (add new columns)
transform(df, :RegionName => ByRow(length) => :RegionNameLength)

# Not tidy
#=
915×5 DataFrame
│ Row │ RegionName                         │ 2020-01-31 │ 2020-02-29 │ 2020-03-31 │ 2020-04-30 │
│     │ String                             │ Float64    │ Float64    │ Float64    │ Float64    │
├─────┼────────────────────────────────────┼────────────┼────────────┼────────────┼────────────┤
│ 1   │ United States                      │ 247060.0   │ 248046.0   │ 249140.0   │ 250271.0   │
│ 2   │ New York, NY                       │ 485111.0   │ 486070.0   │ 486979.0   │ 488002.0   │
│ 3   │ Los Angeles-Long Beach-Anaheim, CA │ 675561.0   │ 680734.0   │ 685911.0   │ 690062.0   │
│ 4   │ Chicago, IL                        │ 242376.0   │ 242743.0   │ 243262.0   │ 243627.0   │
⋮
│ 911 │ Lamesa, TX                         │ 73791.0    │ 73892.0    │ 73978.0    │ 74134.0    │
│ 912 │ Craig, CO                          │ 189339.0   │ 189215.0   │ 189224.0   │ 189312.0   │
│ 913 │ Pecos, TX                          │ 106328.0   │ 107507.0   │ 108640.0   │ 109261.0   │
│ 914 │ Vernon, TX                         │ 69788.0    │ 69656.0    │ 69520.0    │ 69430.0    │
│ 915 │ Ketchikan, AK                      │ 309378.0   │ 309427.0   │ 309458.0   │ 309240.0   │
=#

# Turn to long format
# These are all the same
stack(df, Not(:RegionName); variable_name = :Date)
stack(df, Not(1); variable_name = :Date)
stack(df, 2:5; variable_name = :Date)

# Let's assign to a variable
sdf = stack(df, Not(:RegionName); variable_name = :Date)

describe(sdf, :eltype)
# Interesting type for the Date column
#=
│ Row │ variable   │ eltype                          │
│     │ Symbol     │ DataType                        │
├─────┼────────────┼─────────────────────────────────┤
│ 1   │ RegionName │ String                          │
│ 2   │ Date       │ CategoricalValue{String,UInt32} │
│ 3   │ value      │ Float64                         │
=#

typeof(sdf.Date)
#=
CategoricalArray{String,1,UInt32,String,CategoricalValue{String,UInt32},Union{}}
=#

# Let's fix it.
using Dates
sdf.Date = [Date(get(x)) for x in sdf.Date]

# Convert to wide format again?
unstack(sdf, :Date, :value)

# Group the rows by certain column(s)
groupby(sdf, :Date)

# Summarize the grouped data
combine(groupby(sdf, :Date), :value => mean => :avg)
#=
│ Row │ Date       │ avg      │
│     │ Date       │ Float64  │
├─────┼────────────┼──────────┤
│ 1   │ 2020-01-31 │ 188880.0 │
│ 2   │ 2020-02-29 │ 189514.0 │
│ 3   │ 2020-03-31 │ 190150.0 │
│ 4   │ 2020-04-30 │ 190746.0 │
=#

# Using Pipe.jl to build a transformation pipeline
@pipe sdf |>
    groupby(_, :Date) |>
    combine(_, :value => mean => :avg)

# Row-wise computation
select(df,
    :RegionName,
    ["2020-01-31", "2020-02-29"] => ByRow(-) => :JanFebDiff)

# Joins

# Let’s say we have a reference data frame
country = DataFrame(Name = ["United States"])

# Joining data is easy
innerjoin(df, country, on = [:RegionName => :Name])
leftjoin(df, country, on = [:RegionName => :Name])
rightjoin(df, country, on = [:RegionName => :Name])
outerjoin(df, country, on = [:RegionName => :Name])

semijoin(df, country, on = [:RegionName => :Name])
antijoin(df, country, on = [:RegionName => :Name])


text:
# plotting the results of the experiments

import JLD2
import FileIO: load
import BenchmarkTools: Trial, Parameters
using PGFPlotsX, LinearAlgebra

push!(PGFPlotsX.CUSTOM_PREAMBLE, "\\usepgfplotslibrary{colorbrewer}")

fig_theme = @pgf {
    "cycle list/Dark2-8",
    mark_options = "solid",
    line_width = "1.5pt",
    grid = "major"}

# %% plot potentials and solution

N_unit = 128
N_v = 4
N_c = 5
N_k = 256

r_super, V_grid, W_grid = load("1d/example_$(N_unit)_$(N_v)_$(N_c)_$(N_k).jld2", "r_super", "V_grid", "W_grid")

plot_indices = findall(-7.5 .<= r_super .<= 7.5)

fig = @pgf TikzPicture(Axis(
    {fig_theme...,
    xlabel = "\$r\$", ylabel = "",
    xmin = -7.5, xmax = 7.5,
    ymin = 0, ymax = 12,
    legend_pos = "north west",
    "width = 0.4\\textwidth",
    "height = 0.3\\textwidth",
    scale_only_axis = true},
    PlotInc({"no markers"},
        Table(; x = r_super[plot_indices], y = V_grid[plot_indices])),
    LegendEntry("\$V\$"),
    PlotInc({"no markers"},
        Table(; x = r_super[plot_indices], y = W_grid[plot_indices])),
    LegendEntry("\$W\$"),
    PlotInc({"dashed", "no markers"},
        Table(; x = [0.0, 0.0], y = [0.0, 12.0])),
    LegendEntry("\$\\Omega\$"),
    "\\pgfplotsset{cycle list shift=-1}",
    PlotInc({"dashed", "no markers"},
        Table(; x = [1.5, 1.5], y = [0.0, 12.0]))
))

pgfsave("1d/1d_potentials.tex", fig, include_preamble = false)
pgfsave("1d/1d_potentials.pdf", fig, include_preamble = false)

# %% plot band structure

E_v, E_c, k_bz = load("1d/example_$(N_unit)_$(N_v)_$(N_c)_$(N_k).jld2", "E_v", "E_c", "k_bz")
# ev, ef = load("1d/H_exact_$(N_unit)_$(N_v)_$(N_c)_$(N_k).jld2", "ev", "ef")
ev, ef = load("1d/eigs_$(N_unit)_$(N_v)_$(N_c)_$(N_k)_$(1e-5).jld2", "ev", "ef")

k_1bz = mod.(k_bz .+ 0.5, 1.0) .- 0.5
p = sortperm(k_1bz)

i_gs = findmin(real.(ev))[2]
ef_vck = reshape(abs2.(ef[:, i_gs]), N_v, N_c, N_k)
marker_v = 10 .* sqrt.([sum(ef_vck, dims=2)[iv, 1, ik] for ik in p, iv in 1:N_v])
marker_c = 10 .* sqrt.([sum(ef_vck, dims=1)[1, ic, ik] for ik in p, ic in 1:N_c])

fig = @pgf TikzPicture(Axis(
    {fig_theme...,
    title = "", xlabel = "\$k\$", ylabel = "\$\\epsilon_k\$",
    ylabel_shift = "-3mm",
    xmin = -0.5, xmax = 0.5,
    ymin = -10, ymax = 110,
    legend_pos = "north west",
    "width = 0.4\\textwidth",
    "height = 0.3\\textwidth",
    scale_only_axis = true,
    "colormap/Dark2-8"},
    "\\pgfplotsset{cycle list shift=1}",
    PlotInc({mark = "*", "scatter",
        mark_options = {draw_opacity = 0.0, fill_opacity = 0.3},
        point_meta = "{\\thisrow{markercolor}}",
        "visualization depends on = \\thisrow{markersize} \\as \\perpointmarksize",
        "scatter/@pre marker code/.append style={/tikz/mark size=\\perpointmarksize}"},
        Table(; x = k_1bz[p], y = E_c[1, p], markercolor = [1; 8; 2 * ones(254)], markersize = marker_c[:, 1])),
    LegendEntry("conduction bands"),
    "\\pgfplotsset{cycle list shift=-1}",
    PlotInc({mark = "*", "scatter",
        mark_options = {draw_opacity = 0.0, fill_opacity = 0.3},
        point_meta = "{\\thisrow{markercolor}}",
        "visualization depends on = \\thisrow{markersize} \\as \\perpointmarksize",
        "scatter/@pre marker code/.append style={/tikz/mark size=\\perpointmarksize}"},
        Table(; x = k_1bz[p], y = E_v[4, p], markercolor = [1; 8; 1 * ones(254)], markersize = marker_v[:, 4])),
    LegendEntry("valence bands"),
    "\\pgfplotsset{cycle list shift=-2}",
    PlotInc({"no markers"},
        Table(; x = k_1bz[p], y = E_v[1, p])),
    "\\pgfplotsset{cycle list shift=-3}",
    PlotInc({"no markers"},
        Table(; x = k_1bz[p], y = E_v[2, p])),
    "\\pgfplotsset{cycle list shift=-4}",
    PlotInc({"no markers"},
        Table(; x = k_1bz[p], y = E_v[3, p])),
    "\\pgfplotsset{cycle list shift=-4}",
    PlotInc({"no markers"},
        Table(; x = k_1bz[p], y = E_c[2, p])),
    "\\pgfplotsset{cycle list shift=-5}",
    PlotInc({"no markers"},
        Table(; x = k_1bz[p], y = E_c[3, p])),
    "\\pgfplotsset{cycle list shift=-6}",
    PlotInc({"no markers"},
        Table(; x = k_1bz[p], y = E_c[4, p])),
    "\\pgfplotsset{cycle list shift=-7}",
    PlotInc({"no markers"},
        Table(; x = k_1bz[p], y = E_c[5, p]))
))

pgfsave("1d/1d_band_structure.tex", fig, include_preamble = false)
pgfsave("1d/1d_band_structure.pdf", fig, include_preamble = false)

# %% plot error

N_μ_vec, errors_Z_vv, errors_Z_cc, errors_Z_vc =  load("1d/errors_Z.jld2", "N_μ_vec", "errors_Z_vv", "errors_Z_cc", "errors_Z_vc")

fig = @pgf Axis(
    {fig_theme...,
    title = "Error in ISDF", xlabel = "\$N_\\mu^{ij}\$",
    legend_entries = {"\$Z^{cc}\$", "\$Z^{vv}\$", "\$Z^{vc}\$"},
    xmin = N_μ_vec[1], xmax = N_μ_vec[end], ymin = 1e-10, ymax = 1e0,
    ymode = "log",
    "width = 0.4\\textwidth",
    "height = 0.3\\textwidth",
    scale_only_axis = true},
    PlotInc({"solid", mark = "*"},
        Table(; x = N_μ_vec, y = errors_Z_cc)),
    PlotInc({"dashed", mark = "cube*"},
        Table(; x = N_μ_vec, y = errors_Z_vv)),
    PlotInc({"dotted", mark = "triangle*"},
        Table(; x = N_μ_vec, y = errors_Z_vc))
)

pgfsave("1d/1d_errors_isdf.tex", fig, include_preamble = false)
pgfsave("1d/1d_errors_isdf.pdf", fig, include_preamble = false)

# %% error in the first eigenvalue and spectral function for different N_μ

Z_tol_vec, errors_optical_absorption, errors_ground_state_energy = load("1d/errors_spectrum.jld2", "Z_tol_vec", "errors_optical_absorption", "errors_ground_state_energy")

fig = @pgf LogLogAxis(
    {fig_theme...,
    title = "Error in Spectrum of \$H\$", xlabel = "\$Z_{tol}\$",
    legend_entries = {"error in first eigenvalue", "error in spectral function"},
    xmin = Z_tol_vec[end], xmax = Z_tol_vec[1], ymin = 1e-10, ymax = 5e4,
    # ymode = "log",
    legend_pos = "north west",
    "width = 0.4\\textwidth",
    "height = 0.3\\textwidth",
    scale_only_axis = true},
    PlotInc({"solid", mark = "*"},
        Table(; x = Z_tol_vec, y = errors_optical_absorption)),
    PlotInc({"dashed", mark = "cube*"},
        Table(; x = Z_tol_vec, y = errors_ground_state_energy))
)

pgfsave("1d/1d_errors_spectrum.tex", fig, include_preamble = false)
pgfsave("1d/1d_errors_spectrum.pdf", fig, include_preamble = false)

# %% plot benchmark results

N_k_vec, setup_times, evaluation_times =  load("1d/benchmark.jld2", "N_k_vec", "setup_times", "evaluation_times")

fig = @pgf LogLogAxis(
    {fig_theme...,
    title = "Run Time Scaling", xlabel = "\$N_k\$", ylabel = "run time [s]",
    legend_style = "{at={(0.6,0.2)},anchor=west}",
    ymin = 1e-4,
    "width = 0.4\\textwidth",
    "height = 0.3\\textwidth",
    scale_only_axis = true},
    PlotInc({mark = "*"},
        Table(; x = N_k_vec, y = setup_times)),
    LegendEntry("initial setup"),
    PlotInc({mark = "square*"},
        Table(; x = N_k_vec, y = evaluation_times)),
    LegendEntry("matrix vector product"),
    PlotInc({"dashed"},
        Table(; x = N_k_vec, y = 1.5 * 1e-4 * N_k_vec)),
        LegendEntry("\$O(N_k)\$")
)

pgfsave("1d/1d_timings_k.tex", fig, include_preamble = false)
pgfsave("1d/1d_timings_k.pdf", fig, include_preamble = false)


text:
include("../utils.jl")

@testset "tex_str" begin
    i = 5
    s = t"\$x+##i\$"
    @test s // "\\\$x+5\\\$"
end

@testset "c_str" begin
    @test c"red" == colorant"red"
end

@testset "col2str" begin
    c = c"red"
    s = G.col2str(c)
    @test s // "rgba(1.0,0.0,0.0,1.0)"
    s = G.col2str(c, str=true)
    @test s // "rgba_1_0_0_0_0_0_1_0_"
end

@testset "fl2str" begin
    s = G.fl2str(0.5555555)
    @test s // "0.5556"
    s = G.fl2str(Rational(1, 2))
    @test s // "0.5"
end

@testset "vec2str, dlist" begin
    s = G.vec2str(["ab", "cd"])
    @test s // "\"ab\" \"cd\""
    s = G.vec2str(["ab", "cd"]; sep=",")
    @test s // "\"ab\",\"cd\""

    @test G.vec2str([pi,pi,pi]) // "3.1416 3.1416 3.1416"

    @test G.dlist(1:3) // "d1 d2 d3"
end


text:
function subcircuits(circuit,terminals = "RCLP")
    subcircs = Circuit[]; karva = ""; parameters = [];
    sublength = 0
    coding_length = length(circuit.karva)
    tree = karva_to_tree(circuit.karva,circuit.parameters)
    saplings = subtrees(tree)
    for subtree in saplings
        karva = tree_to_karva(subtree,div((coding_length-1),2),terminals)
        parameters = get_tree_parameters(subtree)
        sublength = length(parameters)
        push!(subcircs,Circuit(karva,vcat(parameters,karva_parameters(karva[sublength+1:end])),nothing))
    end
    return subcircs
end

function subtrees(tree)
    trees = []
    components = listcomponents(tree)
    for component in components
        push!(trees,removecomponent(tree,component))
    end
    return trees
end

function redundacy_testing(circuit,measurements,frequencies,terminals = "RCLP",cutoffratio = 0.80)
    if count(isoperation,circuit.karva[1:3]) < 2
        return circuit
    end
    subcircs = subcircuits(circuit,terminals)
    evaluate_fitness!(subcircs,measurements,frequencies)
    candidate = minimum(subcircs)
    fitnessratio = circuit.fitness/candidate.fitness
    if fitnessratio >= cutoffratio
        return candidate
    else
        return circuit
    end
end

function get_subcircuits(circuit,measurements,frequencies,terminals = "RCLP")
    subcircs = subcircuits(circuit,terminals)
    evaluate_fitness!(subcircs,measurements,frequencies)
    return subcircs
end

function removeredundancy(circuit,measurements,frequencies,terminals = "RCPL",cutoff = 0.80)
    if count(isoperation,circuit.karva[1:3]) == 1
        return circuit
    end
    subcircs = get_subcircuits(circuit,measurements,frequencies,terminals)
    candidate = minimum(subcircs)
    redundancy = true
    while redundancy && count(isoperation,candidate.karva[1:3])>1
        if circuit.fitness/candidate.fitness > cutoff
            circuit = candidate
        else
            redundancy = false
        end
        subcircs = get_subcircuits(circuit,measurements,frequencies,terminals)
        candidate = minimum(subcircs)
    end
    return circuit
end


text:
using HDF5, JLD

module JLDTemp
using HDF5, JLD
include("JLDTest.jl")

function create()
    x = JLDTest(int16(5))  # int16 makes this work on 0.2
    jldopen("require.jld", "w") do file
        addrequire(file, joinpath(Pkg.dir(), "HDF5", "test", "JLDTest.jl"))
        write(file, "x", x, rootmodule="JLDTemp")
    end
end
end

JLDTemp.create()

x = jldopen("require.jld") do file
    read(file, "x")
end
@assert typeof(x) == JLDTest
@assert x.data == 5
rm("require.jld")


text:
using PdfArchiver
using Test

@testset "PdfArchiver.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: rryi/PdfArchiver

max_issues_repo_name: rryi/PdfArchiver

max_forks_repo_name: rryi/PdfArchiver


text:
struct Morlet <: ContWaveClass
    σ::Float64 # \sigma is the time/space trade-off. as sigma->0, the spacial resolution increases; below 5, there is a danger of being non-analytic. Default is 5.8
    κσ::Float64
    cσ::Float64
end
"""
    morl = Morlet(σ::T) where T<: Real

    return the Morlet wavelet with first central frequency parameter σ, which controls the time-frequency trade-off. As σ goes to zero, all of the information becomes spatial. Default is `` 2π``, and it is recommended that you stay within 3-12. Above this range, the overlap between wavelets becomes impractically small. Below it, the mean subtraction term approaches the magnitude of the wavelet, so they become a sum of Gaussians rather than Gaussians.
"""
function Morlet(σ::T) where T <: Real
    κσ = exp(-σ^2 / 2)
    cσ = 1. / sqrt(1 + κσ^2 - 2 * exp(-3 * σ^2 / 4))
    Morlet(σ, κσ, cσ)
end
Morlet() = Morlet(2π)
class(::Morlet) = "Morlet"; name(::Morlet) = "morl"; vanishingmoments(::Morlet) = 0; isAnalytic(::Morlet) = true
const morl = Morlet()
Base.show(io::IO, x::Morlet) = print(io, "Morlet mean $(x.σ)")

struct Morse <: ContinuousWavelets.ContWaveClass
    ga::Float64
    be::Float64
    cf::Float64
end
"""
    morse = Morse(ga::T,be::T,cf::T) where T<: Float64

    return the Morse wavelet with the central frequency parameter cf, gamma parameter ga and beta parameter be.
"""
function Morse_convert(ga::Real, be::Real, cf::Real)
    ga, be, cf = Float64.(ga), Float64.(be), Float64.(cf)
    Morse(ga,be,cf)
end

Morse() = Morse_convert(4, 2, 1)
class(::Morse) = "Morse"; name(::Morse) = "morse";
vanishingmoments(::Morse) = 0; isAnalytic(::Morse) = true
const morse = Morse()
Base.show(io::IO, x::Morse) = print(io, "Morse gamma = $(x.ga), Morse beta = $(x.be), center frequency = $(x.cf)")

# TODO: include a sombrero wavelet, which is dog2.

# Parameterized classes

# abstract type Dog <: ContWaveClass end
# abstract type Paul <: ContWaveClass end

# continuous parameterized
for (TYPE, NAMEBASE, MOMENTS, RANGE, ISAN) in ((:Paul, "paul", -1, 1:20, true), # moments? TODO: is this a good range of parameters?
        (:Dog, "dog",  -1, 0:6, false),)
    @eval begin
        struct $TYPE{N} <: ContWaveClass end# $TYPE end
        class(::$TYPE) = $(string(TYPE))
        name(::$TYPE{N}) where N = string($NAMEBASE, N)
        vanishingmoments(::$TYPE{N}) where N = -1
        order(::$TYPE{N}) where N = N # either order for Paul wavelets, or number of derivatives for DOGs
        Base.show(io::IO, x::$TYPE{N}) where N = print(io, $NAMEBASE * " order $(N)")
        isAnalytic(::$TYPE) = $(ISAN)
    end
    for NUM in RANGE
        CONSTNAME = Symbol(NAMEBASE, NUM)
        @eval begin
            const $CONSTNAME = $TYPE{$NUM}()                  # type shortcut
            export $CONSTNAME
    end
end
end

# this is a bit broken, but could unify these
# abstract type SupaType end
# macro makeType(FAMILY, FORPRINTING, PARAMS, RANGE)
#     parsedParams = mapreduce(x->"$(x[1])::$(x[2])", (a,b) ->a*"\n"*b, @eval($PARAMS))
#     println(parsedParams)
#     @eval begin
#         localParams =
#         struct $FAMILY{N} <: SupaType
#             :($(parsedParams))
#         end
#         family(::$FAMILY) = $(string(FAMILY))
#         order(::$FAMILY{N}) where N = N
#     end
#     # constants for convienence
#     for NUM in RANGE
#         CONSTNAME = Symbol(FAMILY, NUM)
#         @eval begin
#             const $CONSTNAME = $FAMILY{$NUM}()
#             export $CONSTNAME
#         end
#     end
# end
# adapting the orthogonal wavelet classes to the continuous case
struct ContOrtho{OWT} <: ContWaveClass where OWT <: OrthoFilter
    o::OWT
end
ContOrtho(o::W) where W <: WT.OrthoWaveletClass = ContOrtho{typeof(wavelet(o))}(wavelet(o))
class(a::ContOrtho{OWT}) where OWT = "Continuous $(class(a.o))"
name(a::ContOrtho{OWT}) where OWT = "c$(name(a.o))"
vanishingmoments(a::ContOrtho{OWT}) where OWT = vanishingmoments(a.o)
isAnalytic(::ContOrtho) = false
qmf(w::ContOrtho) = w.o.qmf
Base.show(io::IO, x::ContOrtho{W}) where W = print(io, "Continuous $(x.o.name)")

# Single classes
const cHaar = ContOrtho(WT.haar)
const cBeyl = ContOrtho(WT.beyl)
const cVaid = ContOrtho(WT.vaid)
# parametric orthogonal
for (TYPE, NAMEBASE, RANGE) in ((:Daubechies, "Db", 1:10),
        (:Coiflet, "Coif", 2:2:8),
        (:Symlet, "Sym", 4:10),
        (:Battle, "Batt", 2:2:6),)
    for NUM in RANGE
        CONSTNAME = Symbol(string("c", NAMEBASE, NUM))
        @eval begin
            const $CONSTNAME = ContOrtho(WT.$TYPE{$NUM}())        # type shortcut
            export $CONSTNAME
    end
end
end

### averaging types (for now this is already exported by wavelets, so just use that)

abstract type Average end
struct Dirac <: Average end
struct Father <: Average end
struct NoAve <: Average end

function Base.show(io::IO, cf::Dirac)
    print(io, "Dirac")
end
function Base.show(io::IO, cf::Father)
    print(io, "Father Wavelet")
end
function Base.show(io::IO, cf::NoAve)
    print(io, "No Averaging")
end


text:
function measure_feature(property::BasicTopology, t::IndexedTable, labels::AbstractArray, N::Int = 0)
    N = N == 0 ? maximum(labels) : N
    fill_properties(property, t)
end

function fill_properties(property::BasicTopology, t::IndexedTable)
    t = property.holes ? determine_holes(t) : t
    t = property.euler_number ? compute_euler_number(t) : t
end

function compute_euler_number(t)
    @transform t {euler₄ = (1/4)*(:Q₁ - :Q₃ + 2*:Qₓ), euler₈ = (1/4)*(:Q₁ - :Q₃ - 2*:Qₓ)}
end

function determine_holes(t)
    # Number of holes equals the number of connected components (i.e. 1) minus
    # the Euler number.
    @transform t {holes = 1 - (1/4)*(:Q₁ - :Q₃ - 2*:Qₓ)}
end


text:
"""
    function compact_ped(ped)
---
This function is to simplify inbreeding calculation.
With *A* matrix, full-sibs have the same inbreeding value.
This function returns a dataframe, taking the first ID of a full sibship.
There is also a column counting the sibship size.
"""
function compact_ped(ped)
    t = DataFrame(g8n=Int16[], id=Int32[], ic=Float64[])
    pa = ma = count = 0
    for x in eachrow(ped)
        if x.sir != pa && x.dam != ma
            pa, ma = x.sir, x.dam
            push!(t, (x.g8n, x.id, 0.))
        end
    end
    v = combine(groupby(ped, [:g8n, :sir, :dam]), nrow => :count)
    t.count = v.count[2:end]
    t
end

"""
    function summarize(ped, snp, qtl)
---
# Objective
Summarize the results from function `breeding`.
Return results in a DataFrame

The following are calculated
1. mean TBV of each trait
2. genetic variance of each trait
3. mean inbreeding coefficient of each generation based on A matrix

# Arguments
Here,
- `ped` is the pedigree of the production/1st trait, where selection is
also taking place.
- `snp` is the SNP genotype matrix of the production trait.
- `qtl` is the `QTL` vector for the binary (2nd) traits.

# Inbreeding using A matrix
- Recursive method is used, as here we have large full sibship,
    only one calculation is needed.
- Inbreeding calculation is also multiple threaded

# Returns
- mean inbreeding coefficients by generation
- mean TBV by generation
"""
function summarize(ped, snp, qtl)
    bgt = alleles2gt(view(snp, qtl.pos, :))
    ped.clg = bgt'qtl.effect
    rst = combine(groupby(ped, :g8n),
                  :tbv => mean => :mpbv,
                  :tbv => var => :pvg,
                  :clg => mean => :mcbv,
                  :clg => var => :cvg)
    mpd = Matrix(select(ped, :sir, :dam)) # my pedigree, for func `kinship`

    # Inbreeding using A matrix, create a compact pedigree below
    cpd = compact_ped(ped)
    @Threads.threads for i in 1:nrow(cpd)
        cpd[i, :ic] = kinship(mpd, cpd[i, :id], cpd[i,:id]) - 1
    end
    # mean inbreeding coefficients with A matrix
    mica = combine(groupby(cpd, :g8n),
                   [:ic, :count] => ((x, y) -> x'y/sum(y)) => :mica).mica
    rst.mica = [0; mica]
    rst
end

"""
    function sumsum(df)
---
This `sumsum` is to summarize some summarized results.
The results, in a `DataFrame`, are grouped by repeats.
Then, the `combine` function is used.
"""
function sumsum(df)
    gp = groupby(df, :rpt)
    m = copy(gp[1])
    for j in 2:ncol(m)-1
        for i in 2:length(gp)
            m[:, j] += gp[i][:, j]
        end
        m[:, j] ./= last(df).rpt
    end
    m
end

"""
    function sum_qtl_frq(prd, snp, qtl)
---
QTL frequency changes in the production
"""
function sum_qtl_frq(prd, snp, qtl)
    gg = groupby(prd, :g8n)     # group on generation
    fq = Float64[]
    for grp in gg
        fra = 2first(grp).id - 1
        til = 2last(grp).id
        qa = view(snp, qtl.pos, fra:til)
        append!(fq, mean(qa, dims=2))
    end
    reshape(fq, length(gg), :)
end

"""
    function allele_frq_flux(prd, snp, loci)
---
Return the frequency change of `loci` over generations.
"""
function allele_frq_flux(prd, snp, loci)
    gp = groupby(prd, :g8n)
    ng = length(gp)
    df = DataFrame()
    ft = []                     # [[from, to]]
    for df in gp
        fra, til = 2first(df.id) - 1, 2last(df.id)
        push!(ft, [fra, til])
    end
    for l in loci
        frq = Float64[]
        for (fra, til) in ft
            push!(frq, mean(view(snp, l, fra:til)))
        end
        df[!, string(l)] = frq
    end
    df
end

"""
    function bp_summary(prd, snp, qtl, dir)
---
Summarize simulation results of the `b`reeding `p`rogram.
The figures and tables will be saved in `dir`.
"""
function bp_summary(prd, snp, qtl; dir=".")
    isdir(dir) || mkpath(dir)
    # construct pedigree
    n₀ = length(
        unique(
            Matrix(
                select(
                    filter(row ->row.g8n ==1, prd), :sir, :dam))))
    ped = begin                 # this was the way I recorded the pedigree
        n = n₀ + nrow(prd)      # total ID number, including base
        ped = zeros(Int32, n, 2) # the first `n₀` are automatically 0s
        # parents of generation 1 refer to base
        pm₁ = Matrix(select(filter(row -> row.g8n == 1, prd), :sir, :dam))
        n₁ = size(pm₁)[1]
        ped[n₀+1:n₀+n₁, :] = pm₁
        # parents of generation 2+ are about descendants
        pm₂ = Matrix(select(filter(row -> row.g8n > 1, prd), :sir, :dam))
        pm₂ .+= n₀
        ped[n₀+n₁+1:end, :] = pm₂
        ped
    end

    @info "Inbreeding using A matrix"
    begin                 # generation average inbreeding coefficients
        # only need to record generation, and one ID number of its many full sibs.
        n = size(ped)[1]
        g8n = zeros(Int, n)
        g8n[n₀+1:end] = prd.g8n
        df = DataFrame(g8n=Int[], id=Int[], cnt=Int[], ic=Float64[])
        a = b = 0
        for (g, id, pa, ma) in eachrow([g8n 1:n ped])
            g == 0 && continue
            if pa == a && ma == b
                df.cnt[end] += 1
                continue
            end
            push!(df, [g, id, 1, 0])
            a, b = pa, ma
        end
        Threads.@threads for i in 1:nrow(df)
            df.ic[i] = kinship(ped, df.id[i], df.id[i]) - 1
        end
        boxplot(df.g8n, df.ic, label="Inbreeding measured by A")
        savefig("$dir/inbreeding-changes-with-A.pdf")
    end

    @info "Inbreeding using homozygous loci"
    begin
        m = size(snp)[1]
        df = DataFrame(g8n = prd.g8n, ic = zeros(nrow(prd)))
        Threads.@threads for k in 1:nrow(prd)
            df.ic[k] = m - sum(snp[:, 2k-1] .⊻ snp[:, 2k])
        end
        df.ic ./= m
        t = mean(filter(row -> row.g8n == 1, df).ic)
        df.ic ./= t
        df.ic .-= 1
        boxplot(df.g8n, df.ic, label="Inbreeding measured with homozygotes ratio")
        savefig("$dir/inbreeding-changes-with-homozygous-loci.pdf")
    end

    # Something not right below.  Will visit this later.
    #@info "Inbreeding using G matrix"
    #begin
    #    p = mean(snp, dims = 2)
    #    tp = 2p
    #    s2pq = sum((1 .- p) .* tp)
    #    r2pq = 1 / s2pq
    #    df = DataFrame(g8n = prd.g8n, ic = zeros(nrow(prd)))
    #    for i in 1:nrow(prd)
    #        z = (snp[:, 2i-1] + snp[:, 2i]) - tp
    #        df.ic[i] = (z'z)[1]*r2pq - 1
    #    end
    #    boxplot(df.g8n, df.ic, label="Inbreeding measured with G")
    #    savefig("doc/fig/inbreeding-changes-with-G.pdf")
    #end

    gp = groupby(prd, :g8n)
    @info "Genetic variance and inbreeding over generations"
    begin                                 # plot vₐ² trends of 2 traits
        p = combine(gp, :tbv => var => :vprd)
        b = combine(gp, :t2c => var => :vbin)
        plot(p.vprd, label="Production", dpi=300)
        plot!(b.vbin, label="Binary")
        savefig("$dir/Va-changes.pdf")
    end
    @info join(["",
                "QTL frequency changes",
                "  - if negtive a, switch p to 1-p",
                "  - bubble plot, sizes are with respect to initial frq"], '\n')
    begin
        qg = GEAS.alleles2gt(snp[qtl[2].pos, :])
        n1 = nrow(filter(row -> row.g8n == 1, prd))
        p1 = mean(qg[:, 1:n1], dims=2)./2
        n2 = nrow(filter(row -> row.g8n == prd.g8n[end], prd))
        p2 = mean(qg[:, end-n2:end], dims=2)./2
        scatter(p1, p2 .- p1, ms=p1.*10,
                label="Frequency changes",
                xlabel="MAF",
                ylabel="frequency change",
                dpi=300)
        savefig("$dir/frequency-changes.pdf")
    end
end

function sum_improve(prd)
    gp = groupby(prd, :g8n)
    @info "Genetic variance and inbreeding over generations"
    p = combine(gp, :tbv => mean => :mprd)
    b = combine(gp, :t2c => mean => :mbin)
    p, b
end

"""
    function plot_results(df, fig)
---
Summarize results in `DataFrame` `df`.
Results will be plotted in `fig`.
Before I used `JLD2` for data storage.
It was not stable.
So I turned to `CSV`, which needs to be read into a DataFrame `df`.
```
"""
function plot_results(df, fig)
    ng8n = maximum(df.g8n)
    rst = DataFrame()
    for gp in groupby(df, :method)
        for gn in groupby(gp, :g8n)
            tmp = combine(gn,
                          :mpbv => mean => :mpbv,
                          :pvg  => mean => :pvg,
                          :mcbv => mean => :mcbv,
                          :cvg  => mean => :cvg,
                          :mica => mean => :mica)
            tmp.method .= gp.method[1]
            tmp.g8n .= gn.g8n[1]
            append!(rst, tmp)
        end
    end
    bymethod = groupby(rst, :method)
    scalefontsizes(0.65)
    p1 = plot(bymethod[1].g8n, bymethod[1].mpbv,
              ylabel = "Mean TBV, production",
              legend = :topleft,
              label = bymethod[1].method[1],
              dpi = 300,
              xticks = 0:ng8n);
    plot!(p1, bymethod[2].g8n, bymethod[2].mpbv, label = bymethod[2].method[1]);
    plot!(p1, bymethod[3].g8n, bymethod[3].mpbv, label = bymethod[3].method[1]);

    p2 = plot(bymethod[1].g8n, bymethod[1].mcbv,
              ylabel = "Mean TBV, challenge",
              legend = :topleft,
              label = bymethod[1].method[1],
              dpi = 300,
              xticks = 0:ng8n);
    plot!(p2, bymethod[2].g8n, bymethod[2].mcbv, label = bymethod[2].method[1]);
    plot!(p2, bymethod[3].g8n, bymethod[3].mcbv, label = bymethod[3].method[1]);

    p3 = plot(bymethod[1].g8n, bymethod[1].pvg,
              ylabel = "Var(TBV), production",
              legend = :topright,
              label = bymethod[1].method[1],
              xticks = 0:ng8n);
    plot!(p3, bymethod[2].g8n, bymethod[2].pvg, label = bymethod[2].method[1]);
    plot!(p3, bymethod[3].g8n, bymethod[3].pvg, label = bymethod[3].method[1]);

    p4 = plot(bymethod[1].g8n, bymethod[1].cvg,
              xlabel = "Generation",
              ylabel = "Var(TBV), challenge",
              legend = :topright,
              label = bymethod[1].method[1],
              xticks = 0:ng8n);
    plot!(p4, bymethod[2].g8n, bymethod[2].cvg, label = bymethod[2].method[1]);
    plot!(p4, bymethod[3].g8n, bymethod[3].cvg, label = bymethod[3].method[1]);

    p5 = plot(bymethod[1].g8n, bymethod[1].mica,
              ylabel = "Mean inbreeding, pedigree",
              legend = :topleft,
              label = bymethod[1].method[1],
              xticks = 0:ng8n);
    plot!(p5, bymethod[2].g8n, bymethod[2].mica, label = bymethod[2].method[1]);
    plot!(p5, bymethod[3].g8n, bymethod[3].mica, label = bymethod[3].method[1]);

    plot(p1, p3, p5, p2, p4);
    savefig(fig)
    scalefontsizes()            # restore the sizes, or other plots are affected
end

"""
    function table_results(df)
---
This function will table the simulation results.
Return a DataFrame only.
"""
function table_results(df)
    ng8n = maximum(df.g8n)
    rst = DataFrame()
    for gp in groupby(df, :method)
        for gn in groupby(gp, :g8n)
            tmp = combine(gn,
                          :mpbv => mean => :mpbv,
                          :pvg  => mean => :pvg,
                          :mcbv => mean => :mcbv,
                          :cvg  => mean => :cvg,
                          :mica => mean => :mica)
            tmp.method .= gp.method[1]
            tmp.g8n .= gn.g8n[1]
            append!(rst, tmp)
        end
    end
    return rst
end

"""
    function table_binary(df, scenario)
---
Select mean breeding values, and its variance, plus the method used
from `df`.  This function transpose these 3 columns
"""
function table_binary(df)
    out = DataFrame()
    for d in groupby(df, :method)
        dic = Dict("x1" => "g_0")
        t = DataFrame(Matrix(select(d, :mcbv, :cvg))', :auto)
        for i in 2:ncol(t)
            dic["x$i"] = "g_$(i-1)"
        end
        rename!(t, dic)
        t.stats = ["EBV", "Var"]
        t.method = repeat([d.method[1]], 2)
        append!(out, select(t, :method, :stats, r"g"))
    end
    out
end

function table_scenario(df)
    out = DataFrame()
    for r in groupby(df, :scenario)
        t = table_results(r)
        b = table_binary(t)
        b.scenario = repeat([r.scenario[1]], nrow(b))
        append!(out, select(b, :scenario, :method, :stats, r"g"))
    end
    out
end


text:
function +(A::LRMat,B::LRMat)
    assert(A.height == B.height);
    assert(A.width == B.width);

    C = HMat{T}();
    C.height = A.height;
    C.width = A.width;
    C.trg = A.trg;
    C.src = A.src;
    C.level = A.level;

    C.EPS = A.EPS;
    C.MAXRANK = A.MAXRANK;
    C.MINN = A.MINN;

    if A.blockType == LOWRANK
        C.blockType = LOWRANK;
        if B.blockType == LOWRANK
            C.UMat = [A.UMat B.UMat];
            C.VMat = [A.VMat B.VMat];
        else
            (Utmp,Vtmp) = hh2l(B);
            C.UMat = [A.UMat Utmp];
            C.VMat = [A.VMat Vtmp];
        end
        hcompress(C);
    elseif A.blockType == DENSE
        C.blockType = DENSE;
        C.DMat = A.DMat + hh2d(B);
    elseif A.blockType == HMAT
        C.blockType = HMAT;
        C.childHMat = Array(HMat{T},4,4);
        hoffset = 0;
        for i = 1:4
            woffset = 0;
            for j = 1:4
                if B.blockType == LOWRANK
                    BH = HMat();
                    BH.height = A.childHMat[i,j].height;
                    BH.width = A.childHMat[i,j].width;
                    BH.blockType = LOWRANK;
                    BH.UMat = B.UMat[hoffset+(1:BH.height),:];
                    BH.VMat = B.VMat[woffset+(1:BH.width),:];
                    woffset += BH.width;
                    if j == 4
                        hoffset += BH.height;
                    end
                    C.childHMat[i,j] = hadd(A.childHMat[i,j],BH);
                elseif B.blockType == DENSE
                    BH = HMat();
                    BH.height = A.childHMat[i,j].height;
                    BH.width = A.childHMat[i,j].width;
                    BH.blockType = DENSE;
                    BH.EPS = B.EPS;
                    BH.MAXRANK = B.MAXRANK;
                    BH.MINN = B.MINN;
                    BH.DMat = B.DMat[hoffset+(1:BH.height),woffset+(1:BH.width)];
                    woffset += BH.width;
                    if j == 4
                        hoffset += BH.height;
                    end
                    C.childHMat[i,j] = hadd(A.childHMat[i,j],BH);
                else
                    C.childHMat[i,j] = hadd(A.childHMat[i,j],B.childHMat[i,j]);
                end
            end
        end
    end
    return C;
end

function hadd{T<:Number}(A::HMat{T},UMat::Matrix{T},VMat::Matrix{T})
    assert(A.height == size(UMat,1));
    assert(A.width == size(VMat,1));

    C = HMat{T}();
    C.height = A.height;
    C.width = A.width;
    C.trg = A.trg;
    C.src = A.src;
    C.level = A.level;

    C.EPS = A.EPS;
    C.MAXRANK = A.MAXRANK;
    C.MINN = A.MINN;

    if A.blockType == LOWRANK
        C.blockType = LOWRANK;
        C.UMat = [A.UMat UMat];
        C.VMat = [A.VMat VMat];
        hcompress(C);
    elseif A.blockType == DENSE
        C.blockType = DENSE;
        C.DMat = A.DMat + UMat*VMat';
    elseif A.blockType == HMAT
        C.blockType = HMAT;
        C.childHMat = Array(HMat{T},4,4);
        hoffset = 0;
        for i = 1:4
            woffset = 0;
            for j = 1:4
                C.childHMat[i,j] = hadd(A.childHMat[i,j],UMat[hoffset+(1:A.childHMat[i,j].height),:],VMat[woffset+(1:A.childHMat[i,j].width),:]);
                woffset += A.childHMat[i,j].width;
                if j == 4
                    hoffset += A.childHMat[i,j].height;
                end
            end
        end
    end
    return C;
end

function hadd!(A::HMat,UMat::AbstractMatrix,VMat::AbstractMatrix)
    assert(A.height == size(UMat,1));
    assert(A.width == size(VMat,1));

    if A.blockType == LOWRANK
        A.UMat = [A.UMat UMat];
        A.VMat = [A.VMat VMat];
        hcompress(A);
    elseif A.blockType == DENSE
        A.DMat += UMat*VMat';
    elseif A.blockType == HMAT
        hoffset = 0;
        for i = 1:4
            woffset = 0;
            for j = 1:4
                urange = hoffset .+ (1:A.childHMat[i,j].height)
                vrange = woffset .+ (1:A.childHMat[i,j].width)
                hadd!(A.childHMat[i,j],
                      sub(UMat,urange,:),
                      sub(VMat,vrange,:))
                woffset += A.childHMat[i,j].width;
                if j == 4
                    hoffset += A.childHMat[i,j].height;
                end
            end
        end
    end
end

function hadddiag{T<:Number}(A::HMat{T},alpha)
    C = HMat{T}();
    C.height = A.height;
    C.width = A.width;
    C.trg = A.trg;
    C.src = A.src;
    C.level = A.level;

    C.EPS = A.EPS;
    C.MAXRANK = A.MAXRANK;
    C.MINN = A.MINN;

    if A.blockType == LOWRANK
        C.blockType = LOWRANK;
        C.UMat = A.UMat;
        C.VMat = A.VMat;
    elseif A.blockType == DENSE
        C.blockType = DENSE;
        C.DMat = A.DMat;
        if C.trg == C.src
            C.DMat += alpha*eye(C.height);
        end
    elseif A.blockType == HMAT
        C.blockType = HMAT;
        C.childHMat = Array(HMat{T},4,4);
        for i = 1:4, j = 1:4
            C.childHMat[i,j] = hadddiag(A.childHMat[i,j],alpha);
        end
    end
    return C;
end

function hadddiag!(A::HMat,alpha)
    if A.blockType == DENSE
        if A.trg == A.src
            A.DMat += alpha*eye(A.height);
        end
    elseif A.blockType == HMAT
        for i = 1:4
            hadddiag!(A.childHMat[i,i],alpha);
        end
    end
end

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat


text:
#!/usr/bin/env julia
array1 = [(2n + 1)^2 for n in 1:5]
array2 = [sqrt(i) for i in array1]
println(typeof(1:5), " ", typeof(array1), " ", typeof(array2))
1:5, array1, array2
    for j in 1:(n - i)
      if a[j] > a[j + 1]
        a[j], a[j + 1] = a[j + 1], a[j]
      end
    end
  end
  return a
end

data = [65, 51, 32, 12, 23, 84, 68, 1]
bubbleSort!(data)
for i in 1:length(data)
  print(i," ")
  global s
  #This usage of the ‘global‘ keyword is not needed in Jupyter
  #But elsewhere without it:
  ##ERROR: LoadError: UndefVarError: s not defined
  s += beta * data[i]
  data[i] *= -1
end
# print(i) #Would cause ERROR: LoadError: UndefVarError: i not defined

println("\nSum of data in external scope: ", s)

function sumData(beta)
  s = 0
  #try adding the prefix global
  for i in 1:length(data)
    s += data[i] + gamma
  end
  return s
end

println("Sum of data in a function: ", sumData(beta/2))
@show s


text:
function with_midi(f::Function; port::Int, driver_class)
    driver = driver_class()
    try
        println("Opening port $port")
        driver.open_port(port)
        f(driver)
    catch e
        println(e)
    finally
        println("Closing port $port")
        driver.close_port()
    end
end

function transmit_midi_input_to(observable)
    function(midi_in)
        while true
            msg = midi_in.get_message()
            if (msg != nothing)
                ((signal, pitch, velocity), delta_time) = msg
                if signal == SUSTAIN_ON_SIGNAL
                    break
                end
                if signal == NOTE_ON_SIGNAL || signal == NOTE_OFF_SIGNAL
                    R.next!(observable, msg)
                end
            else
                sleep(0.001)
            end
        end
        R.complete!(observable)
    end
end

max_stars_repo_name: miksu/FluidIntonation

max_issues_repo_name: miksu/FluidIntonation

max_forks_repo_name: miksu/FluidIntonation


text:
function execute(problem::Problem{:elasticFlagAggFem}; kwargs...)

   # Problem setting (Default FSI-2)
   println("Setting Elastic flag problem parameters")
   Um = _get_kwarg(:Um,kwargs,1.0)
   H = _get_kwarg(:H,kwargs,0.41)
   ⌀ = _get_kwarg(:D,kwargs,0.1)
   # Solid properties
   E_s = _get_kwarg(:E_s,kwargs,1.4e6)
   ν_s = _get_kwarg(:nu_s,kwargs,0.4)
   ρ_s = _get_kwarg(:rho_s,kwargs,1.0e4)
   # Fluid properties
   ρ_f = _get_kwarg(:rho_f,kwargs,1.0e3)
   Re = _get_kwarg(:Re,kwargs, 100.0)
   μ_f = ρ_f * Um * ⌀ / Re
   γ_f = _get_kwarg(:gamma_f,kwargs,1.0)
   # Mesh properties
   E_m = _get_kwarg(:E_m,kwargs,1.0)
   ν_m = _get_kwarg(:nu_m,kwargs,-0.1)
   α_m = _get_kwarg(:alpha_m,kwargs,1.0e-5)
   weight_strategy = _get_kwarg(:alpha_m_weight,kwargs,"constant")
   # Time stepping
   t0 = _get_kwarg(:t0,kwargs,0.0)
   tf = _get_kwarg(:tf,kwargs,0.1)
   dt = _get_kwarg(:dt,kwargs,0.1)
   θ  = _get_kwarg(:theta,kwargs,0.5)
   # Post-process
   is_vtk = _get_kwarg(:is_vtk,kwargs,false)

  ## Build geometry
  # Cilinder
  R = 0.05
  C = Point(0.2,0.2)
  cylinder = disk(R,x0=C)
  # Flag
  C1 = Point(0.2,0.19)
  d1 = VectorValue(0.4,0.0)
  d2 = VectorValue(0.0,0.02)
  C2 = Point(0.6,0.21)
  n1 = VectorValue(1.0,0.0)
  n2 = VectorValue(0.0,1.0)
  p1 = plane(x0=C1,v=-n1)
  p2 = plane(x0=C1,v=-n2)
  p3 = plane(x0=C2,v=n1)
  p4 = plane(x0=C2,v=n2)
  flag = intersect(intersect(p1,p2),intersect(p3,p4))
  #flag = !(quadrilateral(x0=C1,d1=d1,d2=d2))
  # Fluid and solid domains
  solid_geo = setdiff(flag,cylinder,name="solid")
  not_fluid = union(cylinder,flag)
  fluid_geo = !(not_fluid,name="fluid")

  # Background model
  pmin = Point(0.0,0.0)
  pmax = Point(2.2,0.41)
  partition = (80,40)
  bgmodel = CartesianDiscreteModel(pmin,pmax,partition)
  dp = pmax - pmin
  h = dp[1]/80

  # Fluid and solid models
  #cutgeo = cut(bgmodel,union(fluid_geo,solid_geo))
  cutgeo = cut(bgmodel,fluid_geo)
  model_fluid = DiscreteModel(cutgeo,"fluid")
  #model_solid = DiscreteModel(cutgeo,"solid")

  # Boundary labels
  labels = get_face_labeling(bgmodel)
  add_tag_from_tags!(labels,"wall",[2,4,5,6])
  add_tag_from_tags!(labels,"inlet",[1,3,7])
  add_tag_from_tags!(labels,"outlet",[8])

  # Define BC functions
  println("Defining Boundary conditions")
  u_in(x, t) = VectorValue(1.5 * Um * x[2] * (H - x[2]) / ((H / 2)^2), 0.0)
  u_noSlip(x, t) = VectorValue(0.0, 0.0)
  u_in(t::Real) = x -> u_in(x, t)
  u_noSlip(t::Real) = x -> u_noSlip(x, t)

  # Triangulations
  println("Defining triangulations")
  trian = Triangulation(bgmodel)
  trian_Ω = Triangulation(cutgeo)
  #trian_solid = Triangulation(model_solid)
  trian_fluid = Triangulation(model_fluid)
  #trian_Γi = InterfaceTriangulation(model_fluid,model_solid)
  trian_Γi = EmbeddedBoundary(cutgeo)
  #trian_boundary_Γi = get_left_boundary(trian_Γi)
  n_Γi = get_normal_vector(trian_Γi)
  hΓᵢ = reindex(cell_measure(trian_Γi,trian),trian_Γi)

  # Quadratures
  println("Defining quadratures")
  order = _get_kwarg(:order,kwargs,2)
  degree = 2*order
  bdegree = 2*order
  quad_Ω = CellQuadrature(trian_Ω,2*order)
  #quad_solid = CellQuadrature(trian_solid,degree)
  quad_fluid = CellQuadrature(trian_fluid,degree)
  quad_Γi = CellQuadrature(trian_Γi,bdegree)

  # Cell aggregation
  strategy = AggregateAllCutCells()
  aggregates = aggregate(strategy,cutgeo)

  # Test FE Spaces
  println("Defining FE spaces")
  Vf = TestFESpace(
  model=model_fluid,
  valuetype=VectorValue{2,Float64},
  reffe=:QLagrangian,
  order=order,
  conformity =:H1,
  dirichlet_tags=["inlet", "wall"],
  dof_space=:physical
  )
  # Vs = TestFESpace(
  # model=model_solid,
  # valuetype=VectorValue{2,Float64},
  # reffe=:Lagrangian,
  # order=order,
  # conformity =:H1,
  # dirichlet_tags=["fixed"]
  # )
  Vfser = TestFESpace(
  reffe=:SLagrangian,
  conformity=:L2, # we don't neet to impose continuity since we only use the cell dof basis / shapefuns
  valuetype=VectorValue{2,Float64},
  model=model_fluid,
  order=order,
  dof_space=:physical
  )
  Qf = TestFESpace(
  reffe=:PLagrangian,
  conformity=:L2,
  valuetype=Float64,
  model=model_fluid,
  order=order-1,
  dof_space=:physical
  )
  V = AgFEMSpace(Vf,aggregates,Vfser)
  Q = AgFEMSpace(Qf,aggregates)
  U = TrialFESpace(V,[u_in(0.0),u_noSlip(0.0)])
  P = TrialFESpace(Q)
  Y = MultiFieldFESpace([V,Q])
  X = MultiFieldFESpace([U,P])

  # Stokes operator
  println("Defining Stokes operator")
  function A_Ω(x,y)
    u, p = x
    v, q = y
    ∇(v)⊙∇(u) - q*(∇⋅u) - (∇⋅v)*p
  end
  function B_Ω(y)
    v, q = y
    0.0
  end
  γ = order*(order+1)
  function A_Γd(x,y)
    u, p = x
    v, q = y
    (γ/h)*v⋅u - v⋅(n_Γi⋅∇(u)) - (n_Γi⋅∇(v))⋅u + (p*n_Γi)⋅v + (q*n_Γi)⋅u
  end
  function B_Γd(y)
    v, q = y
    (γ/h)*v⋅u_noSlip(0.0) - (n_Γi⋅∇(v))⋅u_noSlip(0.0) + (q*n_Γi)⋅u_noSlip(0.0)
  end
  t_Ω = LinearFETerm(A_Ω,trian_Ω,quad_Ω)
  t_Γd = AffineFETerm(A_Γd,B_Γd,trian_Γi,quad_Γi)
  # res_ST(x,y) = WeakForms.stokes_residual(x,y,μ_f,VectorValue(0.0,0.0))
  # res_ST_Γd(x,y) = WeakForms.stokes_residual_Γd(x,y,n_Γi,μ_f,γ_f,hΓᵢ,u_noSlip(0.0))
  # jac_ST_Γd(x,dx,y) = WeakForms.stokes_residual_Γd(dx,y,n_Γi,μ_f,γ_f,hΓᵢ,u_noSlip(0.0))
  # t_ST_Ωf = FETerm(res_ST,trian_fluid,quad_fluid)
  # t_ST_Γd = LinearFETerm(res_ST_Γd,trian_Γi,quad_Γi)
  op_ST = FEOperator(X,Y,t_Ω,t_Γd)#,t_ST_Γd)

  # Solve Stokes problem
  xh = solve(op_ST)
  writevtk(trian_fluid,"stokes.vtu",cellfields=["vh"=>restrict(xh[1],trian_fluid),"ph"=>restrict(xh[2],trian_fluid)])
end


text:
# Autogenerated wrapper script for yaml_cpp_jll for armv6l-linux-musleabihf-cxx11
export libyaml_cpp

JLLWrappers.@generate_wrapper_header("yaml_cpp")
JLLWrappers.@declare_library_product(libyaml_cpp, "libyaml-cpp.so.0.7")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libyaml_cpp,
        "lib/libyaml-cpp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This example demonstrates a job scheduling through adding the
# number 100 to every component of the vector data. The root
# assigns one element to each worker to compute the operation.
# When the worker is finished, the root sends another element
# until each element is added 100
# Inspired on
# https://www.hpc.ntnu.no/ntnu-hpc-group/vilje/user-guide/software/mpi-and-mpi-io-training-tutorial/basic-mpi/job-queue

using MPI

function job_queue(data,f)
    MPI.Init()

    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    world_size = MPI.Comm_size(comm)
    nworkers = world_size - 1

    root = 0

    MPI.Barrier(comm)
    T = eltype(data)
    N = size(data)[1]
    send_mesg = Array{T}(undef, 1)
    recv_mesg = Array{T}(undef, 1)

    if rank == root # I am root

        idx_recv = 0
        idx_sent = 1

        new_data = Array{T}(undef, N)
        # Array of workers requests
        sreqs_workers = Array{MPI.Request}(undef,nworkers)
        # -1 = start, 0 = channel not available, 1 = channel available
        status_workers = ones(nworkers).*-1

        # Send message to workers
        for dst in 1:nworkers
            if idx_sent > N
                break
            end
            send_mesg[1] = data[idx_sent]
            sreq = MPI.Isend(send_mesg, dst, dst+32, comm)
            idx_sent += 1
            sreqs_workers[dst] = sreq
            status_workers[dst] = 0
            print("Root: Sent number $(send_mesg[1]) to Worker $dst\n")
        end

        # Send and receive messages until all elements are added
        while idx_recv != N
            # Check to see if there is an available message to receive
            for dst in 1:nworkers
                if status_workers[dst] == 0
                    (flag, status) = MPI.Test!(sreqs_workers[dst])
                    if flag
                        status_workers[dst] = 1
                    end
                end
            end
            for dst in 1:nworkers
                if status_workers[dst] == 1
                    ismessage, status = MPI.Iprobe(dst,dst+32, comm)
                    if ismessage
                        # Receives message
                        MPI.Recv!(recv_mesg, dst, dst+32, comm)
                        idx_recv += 1
                        new_data[idx_recv] = recv_mesg[1]
                        print("Root: Received number $(recv_mesg[1]) from Worker $dst\n")
                        if idx_sent <= N
                            send_mesg[1] = data[idx_sent]
                            # Sends new message
                            sreq = MPI.Isend(send_mesg, dst, dst+32, comm)
                            idx_sent += 1
                            sreqs_workers[dst] = sreq
                            status_workers[dst] = 1
                            print("Root: Sent number $(send_mesg[1]) to Worker $dst\n")
                        end
                    end
                end
            end
        end

        for dst in 1:nworkers
            # Termination message to worker
            send_mesg[1] = -1
            sreq = MPI.Isend(send_mesg, dst, dst+32, comm)
            sreqs_workers[dst] = sreq
            status_workers[dst] = 0
            print("Root: Finish Worker $dst\n")
        end

        MPI.Waitall!(sreqs_workers)
        print("Root: New data = $new_data\n")
    else # If rank == worker
        # -1 = start, 0 = channel not available, 1 = channel available
        status_worker = -1
        while true
            sreqs_workers = Array{MPI.Request}(undef,1)
            ismessage, status = MPI.Iprobe(root, rank+32, comm)

            if ismessage
                # Receives message
                MPI.Recv!(recv_mesg, root, rank+32, comm)
                # Termination message from root
                if recv_mesg[1] == -1
                    print("Worker $rank: Finish\n")
                    break
                end
                print("Worker $rank: Received number $(recv_mesg[1]) from root\n")
                # Apply function (add number 100) to array
                send_mesg = f(recv_mesg)
                sreq = MPI.Isend(send_mesg, root, rank+32, comm)
                sreqs_workers[1] = sreq
                status_worker = 0
            end
            # Check to see if there is an available message to receive
            if status_worker == 0
                (flag, status) = MPI.Test!(sreqs_workers[1])
                if flag
                    status_worker = 1
                end
            end
        end
    end
    MPI.Barrier(comm)
    MPI.Finalize()
end

f = x -> x.+100
data = collect(1:10)
job_queue(data,f)


text:
input = joinpath(@__DIR__, "input")
raw = readlines(input)

struct Rule
    outer
    inner
    quantity
end

rcontainer = r"^([a-z]+ [a-z]+) bags contain "
rbags = r"(\d+) ([a-z]+ [a-z]+)"

function parserule(rule)
    color = match(rcontainer, rule).captures[1]
    return map(m -> Rule(color, m.captures[2], parse(Int, m.captures[1])), collect(eachmatch(rbags, rule)))
end

rules = reduce(append!, map(parserule, raw))

function cancontain(rules, color)
    return map(r -> Pair(r.inner, r.quantity), filter(r -> r.outer == color, rules))
end

function containedby(rules, color)
    return unique(map(r -> r.outer, filter(r -> r.inner == color, rules)))
end

function countcontainers(rules, colors)
    queue = Set(colors)
    tried = Set{String}()
    while !isempty(queue)
        color = pop!(queue)
        more = containedby(rules, color)
        for m in more
            push!(queue, m)
        end
        push!(tried, color)
    end
    return length(tried)
end

function countcontained(rules, colors)
    queue = Set(colors)
    count = 0
    while !isempty(queue)
        next = pop!(queue)
        col = next.first
        qty = next.second
        count += qty
        for (k, v) in cancontain(rules, col)
            push!(queue, Pair(k, qty * v))
        end
    end
    return count
end

p1 = countcontainers(rules, containedby(rules, "shiny gold"))
p2 = countcontained(rules, cancontain(rules, "shiny gold"))

@assert(p1 == 226)
@assert(p2 == 9569)

println("-----------------------------------------------------------------------")
println("handy haversacks -- part one :: $p1")
println("handy haversacks -- part two :: $p2")
println("-----------------------------------------------------------------------")


text:
@testset "Check physical conversion between different quantities" begin

    @test 12 * 20 / 2 / 100 ≈ calc_NA(100, 20, 12)
    @test 1.33 ≈ calc_NA(100e-3, 200e-3, 1.33)

    λ =  500e-9
    n = 1.123
    k = 2 * π / λ * n
    @test calc_k(λ, n) ≈ k
    @test calc_κ(λ, n) ≈ k / 2 / π
    @test calc_λ(k, n) ≈ λ

    n = 1
    λ =  500e-9
    k = 2 * π / λ * n
    @test calc_k(λ) ≈ k
    @test calc_κ(λ) ≈ k / 2 / π
    @test calc_λ(k) ≈ λ

end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
"""
    struct Scope <: Any
Hold a pointer to a `clang::Scope` object.
"""
struct Scope
    ptr::CXScope
end

Base.unsafe_convert(::Type{CXScope}, x::Scope) = x.ptr
Base.cconvert(::Type{CXScope}, x::Scope) = x