text:
#=
# Steady frictional flow through a constant-area duct
This notebook demonstrates the use of tools for computing flow through a constant-area
 duct with frictional effects, also known as **Fanno flow**
=#

# ### Set up the module
using Gasdynamics1D
#-
using Plots
using LaTeXStrings

#=
### Some general aspects of Fanno flow
Most calculations with frictional flow involve the relationship between Mach
number and the reference length $L^*$, representing the length required from that
point to bring the flow to Mach number equal to 1. The relationship is based on
the dimensionless parameter

$$\dfrac{f L^*}{D}$$

where $f$ is the Darcy friction factor and $D$ is the duct diameter. Let's plot
 this relationship for air:
=#
Mrange = range(0.1,8,length=801)
fLDarray = []
for M in Mrange
    push!(fLDarray,value(FLStarOverD(MachNumber(M),FannoFlow,gas=Air)))
end
#-
plot(Mrange,fLDarray,xlim=(0,8),ylim=(0,10),xticks=0:1:8,xlabel="Mach number",ylabel=L"fL^*/D",legend=false)

#=
Notice that the reference length goes to zero at Mach number 1, as it should, by
its definition. It should also be noted that there are a few values of $fL^*/D$ for
which there are two possible Mach numbers---one subsonic and one supersonic. This
is analogous to the situation in isentropic flow with area changes.
=#
#=
### Example 1
Flow enters a duct of diameter 2 cm at Mach number equal to 0.1 and leaves at
Mach number equal to 0.5. Assume that the friction factor is equal to 0.024 throughout.

(a) What is the length of the duct?

(b) How much longer would the duct have to be to reach Mach number equal to 1 at the exit?
=#
#=
First, we use $M_1 = 0.1$ to calculate $f L_1^*/D$. We will do the same with $M_2$
 to calculate $f L_2^*/D$. The actual length $L$ is given by the difference between
 $L_1^*$ and $L_2^*$.
=#
f = FrictionFactor(0.024)
D = Diameter(2u"cm")
#-
M1 = MachNumber(0.1)
fL1star_over_D = FLStarOverD(M1,FannoFlow)
#-
M2 = MachNumber(0.5)
fL2star_over_D = FLStarOverD(M2,FannoFlow)
# Now calculate $L_1^*$ and $L_2^*$, and $L$ from their difference:
L1star = Length(fL1star_over_D*D/f)
L2star = Length(fL2star_over_D*D/f)
L = Length(L1star-L2star)
# The additional length needed to bring the flow to Mach number 1 is, by definition, $L_2^*$:
L2star
# So we only need 0.89 m more to bring the flow to sonic (choked) conditions.

#=
###  Example 2
Consider a flow of air through a circular pipe of diameter 3 cm, starting at
stagnation pressure 200 kPa, stagnation temperature 500 K, and velocity 100 m/s.
The friction factor in the pipe can be assumed to be 0.02 throughout.

(a) What length of pipe would be needed to drive the flow to Mach number 1 at the exit?

(b) What is the mass flow rate through the pipe at those conditions?

(c) Suppose the pipe were 30 m long. How would the mass flow rate change?
=#

# For (a), we are seeking $L_1^*$. First set the known conditions:
D = Diameter(3u"cm")
p01 = StagnationPressure(200u"kPa")
T01 = StagnationTemperature(500)
u1 = Velocity(100)
f = FrictionFactor(0.02)

#=
Now calculate the Mach number at the entrance. For this, we need $c_1$, the speed
of sound. We get this by using the known stagnation temperature $T_{01} = 500$ K
and velocity $u_1 = 100$ m/s. We will calculate the temperature $T_1$ from these
(using $h_1 = h_{01} - u_1^2/2$), and $c_1$ from that:
=#
h01 = StagnationEnthalpy(T01)
h1 = Enthalpy(h01,u1) # this computes h1 = h01 - u1^2/2
T1 = Temperature(h1)
c1 = SoundSpeed(T1)
M1 = MachNumber(u1/c1)
# Now we can compute $fL_1^*/D$ and then $L_1^*$. This is what we seek.
fL1starD = FLStarOverD(M1,FannoFlow)
L1star = Length(fL1starD*D/f)

#=
so we need a pipe of 16.6 m to reach Mach number 1.

Now let's compute the mass flow rate. We need the area and density. For the density,
we will first get the stagnation density (from $p_{01}$ and $T_{01}$) and then get
 density from the isentropic relation. (We are not using isentropic relations
 between two different points in the duct; we are using it to relate *local*
 stagnation conditions to *local* conditions.)
=#
A = Area(D)
ρ01 = StagnationDensity(p01,T01)
ρ1 = Density(ρ01,M1,Isentropic)
# and the mass flow rate is:
ṁ = MassFlowRate(ρ1*u1*A)

#=
(c) Now suppose the pipe is 30 m long. This becomes our new $L_1^*$, and all of
the conditions at the entrance must adjust accordingly. First find the Mach number
=#
Lstar = L = Length(30u"m")
fLstar_over_D = FLOverD(f*Lstar/D)
M1 = SubsonicMachNumber(fLstar_over_D,FannoFlow)

#=
Note that the entrance Mach number is lower. The flow slows down due to the longer pipe.
Now find the new values of $\rho_1$ and $u_1$ (by finding $c_1$):
=#
T1 = Temperature(T01,M1,Isentropic)
ρ1 = Density(ρ01,M1,Isentropic)
c1 = SoundSpeed(T1)
u1 = Velocity(c1,M1)
ṁ = MassFlowRate(ρ1*u1*A)

#=
So the mass flow rate is smaller than it was with the shorter pipe! For the same
reservoir conditions, less mass flow rate is delivered to a longer (choked) pipe.
=#


text:
module RustModels

import NLsolve
import Distributions: wsample
import StatsBase.StatisticalModel
importall DynamicDiscreteModels

# importall Markov
import HiddenMarkovModels: rsm, nsm, z2q, q2z

export sa2y, y2sa, xs2k, k2xs,
		coef!, rand, loglikelihood, dim,mle,
		rustmodel, checkdp

include("rustmodelcore.jl")
include("indexmanipulation.jl")
include("dynamicprogram.jl")
include("classicalrustmodel.jl")

end

ext: jl


text:
#use NTuple because of simpler syntax with CartesianIndex
function idx2pos(i::NTuple{N,T}, dx::NTuple{N,T2}, xmin) where {N,T<:Int,T2<:Real}
    x = SVector(@. (i-0.5) * dx + xmin)
end

function pos2idx_right(x::SVector{N,T}, dx::NTuple{N,T}, xmin) where {N,T<:Real}
    i = convert(NTuple{N,Int}, (@. floor((x-xmin)/dx+0.5)).data )
end
function pos2idx_left(x::SVector{N,T}, dx::NTuple{N,T}, xmin) where {N,T<:Real}
    i = convert(NTuple{N,Int}, (@. floor((x-xmin)/dx+1.5)).data ) #an extra +1 to include the left-most cell
end

function pos2idxNGP(x::SVector{N,T}, dx::NTuple{N,T}, xmin) where {N,T<:Real}
    i = convert(NTuple{N,Int}, (@. ceil((x-xmin)/dx)).data )
end

#we only need to do periodic wrapping for idx (when they fall out of bound of the array)
#no periodic wrapping needed for pos
function periodic_idx(i::T, Ngrid::T, pbc::Bool) where {T<:Int}
    if pbc == false
        return i
    end
    i = i <= 0 ? i + Ngrid : i
    i = i > Ngrid ? i - Ngrid : i
    return i
end

#loop over particles; much faster and doesn't require a tree
function map_particle_to_2Dgrid_loopP(
    field::Vector{T}, volume::Vector{T}, X::Array{SVector{N,T},1}, hsml::Vector{T},
    Xmin::NTuple{N,T}, Xmax::NTuple{N,T}; xaxis::Int=1,yaxis::Int=2, column::Bool=true,
    ngrids::NTuple{N,Int}=(100,100,100), pbc::NTuple{N,Bool}=(true,true,true)) where {N,T}

    map = zeros(T, ngrids[xaxis], ngrids[yaxis]); #each thread has its own slice
    ΔX = @. (Xmax - Xmin) / ngrids
    @show ΔX

    #in the special case of a slice (Nx/y/z = 1), the following hack can speed up a lot without touching the code
    xmin,xmax,Δx=trick(Xmin,Xmax,ΔX,ngrids)

    @time @inbounds for p in eachindex(X)
        Istart = CartesianIndex(pos2idx_left( X[p] .- hsml[p], Δx, xmin))
        Iend   = CartesianIndex(pos2idx_right(X[p] .+ hsml[p], Δx, xmin))

        hsml_inv = 1.0 / hsml[p]
        f_v = field[p] * volume[p] * hsml_inv^3

        for Idx in Istart:Iend
            #@show Idx
            #only for the periodic dimensions
            iw = periodic_idx.(Idx.I, ngrids, pbc)

            #drop the out-of-bound idx (can only happen along the non-periodic dimensions after wrapping)
            if any(iw .<= 0) || any(iw .> ngrids) continue end

            #calculate distance using unwrapped idx
            r = norm(idx2pos(Idx.I,Δx, xmin) - X[p]) * hsml_inv

            map[iw[xaxis],iw[yaxis]] += kernel_cubic(r) * f_v
        end
    end
    res = map
    losdim = findfirst(((1,2,3).!=xaxis) .& ((1,2,3).!=yaxis)) #line-of-sight dimension

    #If this is a slice, we can't do column density, so use column=false instead.
    column = ngrids[losdim] == 1 ? false : column

    fac = column ? Δx[losdim] : 1.0 / ngrids[losdim]
    return res.*=fac
end

#loop over particles; much faster and doesn't require a tree
function map_particle_to_2Dgrid_loopP_thread(
    field::Vector{T}, volume::Vector{T}, X::Array{SVector{N,T},1}, hsml::Vector{T},
    Xmin::NTuple{N,T}, Xmax::NTuple{N,T}; xaxis::Int=1,yaxis::Int=2, column::Bool=true,
    ngrids::NTuple{N,Int}=(100,100,100), pbc::NTuple{N,Bool}=(true,true,true)) where {N,T}

    map = zeros(T, ngrids[xaxis], ngrids[yaxis], nthreads()); #each thread has its own slice
    ΔX = @. (Xmax - Xmin) / ngrids

    #in the special case of a slice (Nx/y/z = 1), the following hack can speed up a lot without touching the code
    xmin,xmax,Δx=trick(Xmin,Xmax,ΔX,ngrids)

    @time @inbounds @threads for p in eachindex(X)
        Istart = CartesianIndex(pos2idx_left( X[p] .- hsml[p], Δx, xmin))
        Iend   = CartesianIndex(pos2idx_right(X[p] .+ hsml[p], Δx, xmin))

        hsml_inv = 1.0 / hsml[p]
        f_v = field[p] * volume[p] * hsml_inv^3

        for Idx in Istart:Iend
            #only for the periodic dimensions
            iw = periodic_idx.(Idx.I, ngrids, pbc)

            #drop the out-of-bound idx (can only happen along the non-periodic dimensions after wrapping)
            if any(iw .<= 0) || any(iw .> ngrids) continue end

            #calculate distance using unwrapped idx
            r = norm(idx2pos(Idx.I,Δx, xmin) - X[p]) * hsml_inv

            map[iw[xaxis],iw[yaxis],threadid()] += kernel_cubic(r) * f_v
        end
    end
    res = dropdims(sum(map, dims=3),dims=3) #sum over threads
    losdim = findfirst(((1,2,3).!=xaxis) .& ((1,2,3).!=yaxis)) #line-of-sight dimension

    #If this is a slice, we can't do column density, so use column=false instead.
    column = ngrids[losdim] == 1 ? false : column

    fac = column ? Δx[losdim] : 1.0 / ngrids[losdim]
    return res.*=fac
end

#loop over particles; much faster and doesn't require a tree
function map_particle_to_3Dgrid_loopP_thread(
    field::Vector{T}, volume::Vector{T}, X::Array{SVector{N,T},1}, hsml::Vector{T},
    xmin::NTuple{N,T}, xmax::NTuple{N,T};
    ngrids::NTuple{N,Int}=(100,100,100),
    pbc::NTuple{N,Bool}=(true,true,true)) where {N,T}

    map = zeros(T, ngrids[1], ngrids[2], ngrids[3], nthreads()); #each thread has its own slice
    Δx = @. (xmax - xmin) / ngrids

    if findfirst(ngrids.==1) != nothing
        println("ngrids=", ngrids)
        error("this is a slice! use map_particle_to_2Dgrid_loopP_thread() instead...")
    end

    @time @inbounds @threads for p in eachindex(X)
        Istart = CartesianIndex(pos2idx_left( X[p] .- hsml[p], Δx, xmin))
        Iend   = CartesianIndex(pos2idx_right(X[p] .+ hsml[p], Δx, xmin))

        hsml_inv = 1.0 / hsml[p]
        f_v = field[p] * volume[p] * hsml_inv^3

        for Idx in Istart:Iend
            #only for the periodic dimensions
            iw = periodic_idx.(Idx.I, ngrids, pbc)

            #drop the out-of-bound idx (can only happen along the non-periodic dimensions after wrapping)
            if any(iw .<= 0) || any(iw .> ngrids) continue end

            #calculate distance using unwrapped idx
            r = norm(idx2pos(Idx.I,Δx, xmin) - X[p]) * hsml_inv

            idx = CartesianIndex(iw)
            map[idx,threadid()] += kernel_cubic(r) * f_v
        end
    end
    return dropdims(sum(map, dims=4),dims=4) #sum over threads
end

function trick(xmin::NTuple{N,T}, xmax::NTuple{N,T}, Δx::NTuple{N,T}, ngrids::NTuple{N,Int}) where {N,T}
    islice = findfirst(ngrids.==1)  #assume there is only one element = 1, which is the slice dimension
    if islice != nothing #we're doing a slice
        mask = ((1,0,0),(0,1,0),(0,0,1))[islice]
        xmin = xmin .- mask.*maximum(Δx).*1e3
        xmax = xmax .+ mask.*maximum(Δx).*1e3
        Δx = @. (xmax - xmin) / ngrids
        @show Δx, (xmax .+ xmin) ./ 2
    end
    #strangely this doesn't work in the thread version (which got even slower)
    #okay, in order to be fast, xmin,xmax,Δx need to have explicit type declaration (but why?)
    #turns out that we just need a diferent name Xmin,Xmax,ΔX (but why?)
    return xmin,xmax,Δx
end

#loop over particles; much faster and doesn't require a tree
function map_particle_to_2Dgrid_loopP_noCar(
    field::Vector{T}, volume::Vector{T}, X::Array{SVector{N,T},1}, hsml::Vector{T},
    Xmin::NTuple{N,T}, Xmax::NTuple{N,T}; xaxis::Int=1,yaxis::Int=2, column::Bool=true,
    ngrids::NTuple{N,Int}=(100,100,100), pbc::NTuple{N,Bool}=(true,true,true)) where {N,T}

    map = zeros(T, ngrids[xaxis], ngrids[yaxis]); #each thread has its own slice
    ΔX = @. (Xmax - Xmin) / ngrids
    @show ΔX

    #in the special case of a slice (Nx/y/z = 1), the following hack can speed up a lot without touching the code
    xmin,xmax,Δx=trick(Xmin,Xmax,ΔX,ngrids)

    @time @inbounds for p in eachindex(X)
    #@time @inbounds for i in eachindex(X)
        istart = pos2idx_left( X[p] .- hsml[p], Δx, xmin)
        iend   = pos2idx_right(X[p] .+ hsml[p], Δx, xmin)
        hsml_inv = 1.0 / hsml[p]
        f_v = field[p] * volume[p] * hsml_inv^3
        for k in istart[3]:iend[3], j in istart[2]:iend[2], i in istart[1]:iend[1]
            #@show (i,j,k)
            #only for the periodic dimensions
            iw = periodic_idx.((i,j,k), ngrids, pbc)

            #drop the out-of-bound idx (can only happen along the non-periodic dimensions after wrapping)
            if any(iw .<= 0) || any(iw .> ngrids) continue end

            #calculate distance using unwrapped idx
            r = norm(idx2pos((i,j,k),Δx, xmin) - X[p]) * hsml_inv

            map[iw[xaxis],iw[yaxis]] += kernel_cubic(r) * f_v
        end
    end
    res = map
    losdim = findfirst(((1,2,3).!=xaxis) .& ((1,2,3).!=yaxis)) #line-of-sight dimension

    #If this is a slice, we can't do column density, so use column=false instead.
    column = ngrids[losdim] == 1 ? false : column

    fac = column ? Δx[losdim] : 1.0 / ngrids[losdim]
    return res.*=fac
end

#nearest-grid-point (NGP): super fast but super noisy too
function map_particle_to_3Dgrid_NGP(field::Vector{T}, mass::Vector{T},
        X::Array{SVector{N,T},1}, ngrids::NTuple{N,Int}, xmin::NTuple{N,T}, xmax::NTuple{N,T}) where {N,T}

    map = zeros(T, ngrids[1], ngrids[2], ngrids[3]); #each thread has its own slice
    counts = zeros(T, ngrids[1], ngrids[2], ngrids[3]); #each thread has its own slice
    Δx = @. (xmax - xmin) / ngrids

    @time @inbounds for i in eachindex(X)
        #idx = periodic_idx(pos2idx(X[i], Δx))
        idx = (pos2idxNGP(X[i], Δx, xmin))
        if any(idx .<= 0) || any(idx .> ngrids) continue end
        #Idx = CartesianIndex(idx.data)
        Idx = CartesianIndex(idx)
        #map[Idx] += mass[i] / Δx^3
        map[Idx] += field[i]
        counts[Idx] += 1.0
    end
    map[counts.>0] ./= counts[counts.>0]
    return map
end

#nearest-grid-point (NGP): super fast but super noisy too
function map_particle_to_3Dgrid_NGP_thread(field::Vector{T}, mass::Vector{T},
        X::Array{SVector{N,T},1}, ngrids::NTuple{N,Int}, xmin::NTuple{N,T}, xmax::NTuple{N,T}) where {N,T}

    map = zeros(T, ngrids[1], ngrids[2], ngrids[3], nthreads()); #each thread has its own slice
    counts = zeros(T, ngrids[1], ngrids[2], ngrids[3], nthreads()); #each thread has its own slice
    Δx = @. (xmax - xmin) / ngrids

    @time @inbounds @threads for i in eachindex(X)
        #idx = periodic_idx(pos2idx(X[i], Δx))
        idx = (pos2idxNGP(X[i], Δx, xmin))
        if any(idx .<= 0) || any(idx .> ngrids) continue end
        #Idx = CartesianIndex(idx.data)
        Idx = CartesianIndex(idx)
        #map[Idx] += mass[i] / Δx^3
        map[Idx,threadid()] += field[i]
        counts[Idx,threadid()] += 1.0
    end
    #map[counts.>0] ./= counts[counts.>0]

    map_all = dropdims(sum(map, dims=4),dims=4) #sum over threads
    counts_all = dropdims(sum(counts, dims=4),dims=4) #sum over threads
    map_all[counts_all.>0] ./= counts_all[counts_all.>0]
    return map_all
end

#=
#loop over pixels, very slow and need a tree
function map_particle_to_grid(field::Vector{T}, volume::Vector{T},
        X::Array{SVector{N,T},1}, hsml::Vector{T}, boxsizes::SVector{N,T}, tree::Node{N,T}) where {N,T}
    map = zeros(T, Ngrid_x, Ngrid_y, Ngrid_z);
    Δx = BOXSIZE_X / Ngrid_x
    Δy = BOXSIZE_Y / Ngrid_y
    Δz = BOXSIZE_Z / Ngrid_z

    @time for i in CartesianIndices(map)
        #@show i.I
        Xgrid = SVector(i.I .* (Δx,Δy,Δz))
        #@show Xgrid
        idx_ngbs = get_scatter_ngb_tree(Xgrid, tree, boxsizes)
        for k in eachindex(idx_ngbs)
            j = idx_ngbs[k]
            dx = nearest.(X[j] - Xgrid, boxsizes)
            dr = norm(dx)
            Wij = kernel_cubic(dr/hsml[j]) / hsml[j]^N #divide by hsml[i]^N after the loop
            #if i == 1 @show typeof.((mass[j] , rho[j] , Wij, dr, dx)) end
            map[i] += field[j] * (volume[j] * Wij)
            #aaa = @SMatrix zeros(T,N,N)
            #aaa += 1.0
        end
    end
    return map
end
=#


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

#####################
# lattice utilities #
#####################

function rewrap(@nospecialize(t), @nospecialize(u))
    if isa(t, TypeVar) || isa(t, Type)
        return rewrap_unionall(t, u)
    end
    return t
end

isType(@nospecialize t) = isa(t, DataType) && t.name === _TYPE_NAME

# true if Type{T} is inlineable as constant T
# requires that T is a singleton, s.t. T == S implies T === S
isconstType(@nospecialize t) = isType(t) && hasuniquerep(t.parameters[1])

# test whether type T has a unique representation, s.t. T == S implies T === S
function hasuniquerep(@nospecialize t)
    # typeof(Bottom) is special since even though it is a leaftype,
    # at runtime, it might be Type{Union{}} instead, so don't attempt inference of it
    t === typeof(Union{}) && return false
    t === Union{} && return true
    isa(t, TypeVar) && return false # TypeVars are identified by address, not equality
    iskindtype(typeof(t)) || return true # non-types are always compared by egal in the type system
    isconcretetype(t) && return true # these are also interned and pointer comparable
    if isa(t, DataType) && t.name !== Tuple.name && !isvarargtype(t) # invariant DataTypes
        return _all(hasuniquerep, t.parameters)
    end
    return false
end

function has_nontrivial_const_info(@nospecialize t)
    isa(t, PartialStruct) && return true
    return isa(t, Const) && !isdefined(typeof(t.val), :instance) && !(isa(t.val, Type) && hasuniquerep(t.val))
end

# Subtyping currently intentionally answers certain queries incorrectly for kind types. For
# some of these queries, this check can be used to somewhat protect against making incorrect
# decisions based on incorrect subtyping. Note that this check, itself, is broken for
# certain combinations of `a` and `b` where one/both isa/are `Union`/`UnionAll` type(s)s.
isnotbrokensubtype(@nospecialize(a), @nospecialize(b)) = (!iskindtype(b) || !isType(a) || hasuniquerep(a.parameters[1]))

argtypes_to_type(argtypes::Array{Any,1}) = Tuple{anymap(widenconst, argtypes)...}

function isknownlength(t::DataType)
    isvatuple(t) || return true
    return length(t.parameters) > 0 && isa(unwrap_unionall(t.parameters[end]).parameters[2], Int)
end

# test if non-Type, non-TypeVar `x` can be used to parameterize a type
function valid_tparam(@nospecialize(x))
    if isa(x, Tuple)
        for t in x
            isa(t, Symbol) || isbits(t) || return false
        end
        return true
    end
    return isa(x, Symbol) || isbits(x)
end

# return an upper-bound on type `a` with type `b` removed
# such that `return <: a` && `Union{return, b} == Union{a, b}`
function typesubtract(@nospecialize(a), @nospecialize(b))
    if a <: b && isnotbrokensubtype(a, b)
        return Bottom
    end
    if isa(a, Union)
        return Union{typesubtract(a.a, b),
                     typesubtract(a.b, b)}
    end
    return a # TODO: improve this bound?
end

function tvar_extent(@nospecialize t)
    while t isa TypeVar
        t = t.ub
    end
    return t
end

_typename(@nospecialize a) = Union{}
_typename(a::TypeVar) = Core.TypeName
function _typename(a::Union)
    ta = _typename(a.a)
    tb = _typename(a.b)
    ta === tb && return ta # same type-name
    (ta === Union{} || tb === Union{}) && return Union{} # threw an error
    (ta isa Const && tb isa Const) && return Union{} # will throw an error (different type-names)
    return Core.TypeName # uncertain result
end
_typename(union::UnionAll) = _typename(union.body)
_typename(a::DataType) = Const(a.name)

function tuple_tail_elem(@nospecialize(init), ct::Vector{Any})
    t = init
    for x in ct
        # FIXME: this is broken: it violates subtyping relations and creates invalid types with free typevars
        t = tmerge(t, tvar_extent(unwrapva(x)))
    end
    return Vararg{widenconst(t)}
end

function countunionsplit(atypes::Union{SimpleVector,Vector{Any}})
    nu = 1
    for ti in atypes
        if isa(ti, Union)
            nu, ovf = Core.Intrinsics.checked_smul_int(nu, unionlen(ti::Union))
            if ovf
                return typemax(Int)
            end
        end
    end
    return nu
end

# take a Tuple where one or more parameters are Unions
# and return an array such that those Unions are removed
# and `Union{return...} == ty`
function switchtupleunion(@nospecialize(ty))
    tparams = (unwrap_unionall(ty)::DataType).parameters
    return _switchtupleunion(Any[tparams...], length(tparams), [], ty)
end

function _switchtupleunion(t::Vector{Any}, i::Int, tunion::Vector{Any}, @nospecialize(origt))
    if i == 0
        tpl = rewrap_unionall(Tuple{t...}, origt)
        push!(tunion, tpl)
    else
        ti = t[i]
        if isa(ti, Union)
            for ty in uniontypes(ti::Union)
                t[i] = ty
                _switchtupleunion(t, i - 1, tunion, origt)
            end
            t[i] = ti
        else
            _switchtupleunion(t, i - 1, tunion, origt)
        end
    end
    return tunion
end

# unioncomplexity estimates the number of calls to `tmerge` to obtain the given type by
# counting the Union instances, taking also into account those hidden in a Tuple or UnionAll
function unioncomplexity(u::Union)
    inner = max(unioncomplexity(u.a), unioncomplexity(u.b))
    return inner == 0 ? 0 : 1 + inner
end
function unioncomplexity(t::DataType)
    t.name === Tuple.name || return 0
    c = 1
    for ti in t.parameters
        ci = unioncomplexity(ti)
        if ci > c
            c = ci
        end
    end
    return c
end
unioncomplexity(u::UnionAll) = max(unioncomplexity(u.body), unioncomplexity(u.var.ub))
unioncomplexity(@nospecialize(x)) = 0

function improvable_via_constant_propagation(@nospecialize(t))
    if isconcretetype(t) && t <: Tuple
        for p in t.parameters
            p === DataType && return true
        end
    end
    return false
end

max_stars_repo_name: Dictino/julia

max_issues_repo_name: Dictino/julia

max_forks_repo_name: Dictino/julia


text:
# Implement some extensions to Euclidean and Squared-Euclidean distances.
for (d, D) in [(:sqeuclidean, :SqEuclidean), (:euclidean, :Euclidean)]
    @eval begin
        function ew(d::$D, x::AV{<:Real}, x′::AV{<:Real})
            return colwise($D(dtol), reshape(x, 1, :), reshape(x′, 1, :))
        end
        function pw(d::$D, x::AV{<:Real}, x′::AV{<:Real})
            return pw($D(dtol), reshape(x, 1, :), reshape(x′, 1, :); dims=2)
        end

        ew(::$D, x::AV{T}) where {T<:Real} = zeros(T, length(x))
        pw(::$D, x::AV{<:Real}) = pairwise($D(dtol), x, x)

        ew(::$D, x::ColVecs{<:Real}, x′::ColVecs{<:Real}) = colwise($D(dtol), x.X, x′.X)
        pw(::$D, x::ColVecs{<:Real}, x′::ColVecs{<:Real}) = pw($D(dtol), x.X, x′.X; dims=2)

        ew(::$D, x::ColVecs{T}) where {T<:Real} = zeros(T, length(x))
        pw(::$D, x::ColVecs{<:Real}) = pairwise($D(dtol), x.X; dims=2)
    end
end

function rrule(::typeof(pairwise), ::Euclidean, x::AbstractVector{<:Real})
    D, pb = Zygote.pullback(
        x->sqrt.(max.(pairwise(SqEuclidean(dtol), x), eps(eltype(x))^2)), x,
    )

    function pairwise_Euclidean_pullback(Δ)
        return (NO_FIELDS, DoesNotExist(), first(pb(Δ)))
    end

    return D, pairwise_Euclidean_pullback
end


text:
using SparseTensors
using Test
using InteractiveUtils: subtypes

TESTLIST = [
    "AbstractSparseTensor" => "abstract.jl",
    "SparseTensorCOO" => "coo.jl",
    "SparseTensorDOK" => "dok.jl",
]

@testset verbose = true "$name" for (name, path) in TESTLIST
    include(path)
end


text:
struct Eulerian
end

function calcEdgeVelocities(vm::Eulerian,
                            state::HydroSnapshot)::Array{Float64, 1}
    return [0.0 for x in 1:length(state.grid)]
end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
#sample (x,y) from a matrix of joint probabilities.
function wsample2(mu)
	dx,dy=size(mu)
	ind2sub((dx,dy),wsample(1:dx*dy,vec(mu)))
end

function rand(model::DynamicDiscreteModel,T::Int)
	#throw error if not calibrated
	data=Array(Int,T)
	dx,dy=size(model.mu)
	x,data[1]=wsample2(model.mu)
	for t=2:T
		x,data[t]=wsample2(reshape(model.m[x,data[t-1],:,:],(dx,dy)))
	end
	data
end

#simulate iid individuals with heterogeneous number of periods of observation
function rand(model::DynamicDiscreteModel,T::Array{Int,1})
	n=length(T)
	data=Array(Array,n)
	for i=1:n
		data[i]=rand(model,T[i])
	end
	data
end

# (T,n) -> n iid individuals with T periods of observation each
rand(model::DynamicDiscreteModel,T,n)=rand(model,fill(T,n))
# rand(model::StatisticalModel,Tn::Tuple{Int,Int})=rand(model,Tn[1],Tn[2])


text:
# based on the tcp part of DonHo's bot
# https://github.com/DonHonerbrink/bruhbot

using Base: StatusUninit, StatusInit, StatusConnecting, StatusOpen,
    StatusClosing, StatusClosed, StatusActive, StatusEOF, StatusPaused

@enum(Status, uninit=StatusUninit, init=StatusInit, connecting=StatusConnecting,
      opened=StatusOpen, act=StatusActive, closing=StatusClosing,
      closed=StatusClosed, Eof=StatusEOF, paused=StatusPaused)

function tcp_create_sock(tcp_sock::TCPSocket)
    tcp_sock = TCPSocket()
end

function tcp_close_sock(tcp_sock::TCPSocket)
    tcp_sock.status = StatusClosed
    return Status(tcp_sock.status)
end

function tcp_connect(tcp_sock::TCPSocket, hostname::String, port::Integer)
    Sockets.connect!(tcp_sock, hostname, port)
    return Status(tcp_sock.status)
end

function tcp_send(tcp_sock::TCPSocket, buffer::String)
    write(tcp_sock, buffer)
    return Status(tcp_sock.status)
end

function tcp_recv(tcp_sock::TCPSocket)
    line = ""
    line = seek_sock(tcp_sock)
    if line != ""
        return Status(tcp_sock.status), line
    else
        return Status(tcp_sock.status), nothing
    end
end

function seek_sock(tcp_sock::TCPSocket)
    line = ""
    @async begin
        while !eof(tcp_sock)
            line = readline(tcp_sock)
            println(stdout, line)
        end
    end
    return line
end

function read_msg(msg::Vector{UInt8})::Union{String, Vector{UInt8}}
    result = convert(String, msg)
    if !isvalid(String, result)
        return msg
    end
    return result
end

max_stars_repo_name: Iaggelis/bald-bot

max_issues_repo_name: Iaggelis/bald-bot

max_forks_repo_name: Iaggelis/bald-bot


text:
function ChainRulesCore.rrule(
    ::typeof(*), A::AbstractSparseMatrix{T}, B::DenseInputVecOrMat
) where T
    C = A * B
    function mul_pullback(C̄)
        _, i_A, i_B, _, _ = (~imul!)(AD.GVar(C, unthunk(C̄)), AD.GVar(A), AD.GVar(B), AD.GVar(1.), AD.GVar(1.))
        return ChainRulesCore.NoTangent(), AD.grad(i_A), AD.grad(i_B)
    end
    return C, mul_pullback
end

function ChainRulesCore.rrule(
    ::typeof(*), xA::Adjoint{T, <:AbstractSparseMatrix}, B::DenseInputVecOrMat
) where T
    C = xA * B
    function mul_pullback(C̄)
        _, i_xA, i_B, _, _ = (~imul!)(AD.GVar(C, C̄), AD.GVar(xA), AD.GVar(B), AD.GVar(1.), AD.GVar(1.))
        return ChainRulesCore.NoTangent(), AD.grad(i_xA), AD.grad(i_B)
    end
    return C, mul_pullback
end

function ChainRulesCore.rrule(
    ::typeof(*), X::DenseMatrixUnion,A::AbstractSparseMatrix{T}
) where T
    C = X * A
    function mul_pullback(C̄)
        _, i_X, i_A, _, _ = (~imul!)(AD.GVar(C, C̄), AD.GVar(X), AD.GVar(A), AD.GVar(1.), AD.GVar(1.))
        return ChainRulesCore.NoTangent(), AD.grad(i_X), AD.grad(i_A)
    end
    return C, mul_pullback
end

function ChainRulesCore.rrule(
    ::typeof(*), X::Adjoint{T1, <:DenseMatrixUnion}, A::AbstractSparseMatrix{T2}
) where {T1, T2}
    C = X * A
    function mul_pullback(C̄)
        _, i_X, i_A, _, _ = (~imul!)(AD.GVar(C, C̄), AD.GVar(X), AD.GVar(A), AD.GVar(1.), AD.GVar(1.))
        return ChainRulesCore.NoTangent(), AD.grad(i_X), AD.grad(i_A)
    end
    return C, mul_pullback
end

function ChainRulesCore.rrule(
    ::typeof(dot), A::AbstractSparseMatrix{T},B::AbstractSparseMatrix{T}
) where T
    r = dot(A, B)
    function dot_pullback(r̄)
        _, i_A, i_B = (~idot)(AD.GVar(r, r̄), AD.GVar(A), AD.GVar(B))
        return ChainRulesCore.NoTangent(), AD.grad(i_A), AD.grad(i_B)
    end
    return r, dot_pullback
end

function ChainRulesCore.rrule(
    ::typeof(dot), x::AbstractVector, A::AbstractSparseMatrix{T1}, y::AbstractVector{T2}
) where {T1, T2}
    r = dot(x, A, y)
    function dot_pullback(r̄)
        _, i_x, i_A, i_y = (~idot)(AD.GVar(r, r̄), AD.GVar(x), AD.GVar(A), AD.GVar(y))
        return ChainRulesCore.NoTangent(), AD.grad(i_x), AD.grad(i_A), AD.grad(i_y)
    end
    return r, dot_pullback
end

function ChainRulesCore.rrule(
    ::typeof(dot), x::SparseVector, A::AbstractSparseMatrix{T1}, y::SparseVector{T2}
) where {T1, T2}
    r = dot(x, A, y)
    function dot_pullback(r̄)
        _, i_x, i_A, i_y = (~idot)(AD.GVar(r, r̄), AD.GVar(x), AD.GVar(A), AD.GVar(y))
        return ChainRulesCore.NoTangent(), AD.grad(i_x), AD.grad(i_A), AD.grad(i_y)
    end
    return r, dot_pullback
end

function _qr_pullback(Ȳ::Tangent, F, A)
    ∂X = qr_rev(F.Q, F.R, Ȳ.Q, Ȳ.R, A)
    return (NoTangent(), ∂X)
end
_qr_pullback(Ȳ::AbstractThunk, F, A) = _qr_pullback(unthunk(Ȳ), F, A)
function ChainRulesCore.rrule(::typeof(private_qr), X::AbstractArray)
    F = private_qr(X)
    qr_pullback(ȳ) = _qr_pullback(ȳ, F, X)
    return F, qr_pullback
end

function ChainRulesCore.rrule(::typeof(getproperty), F::T, x::Symbol) where T <: Normal_QR
    function getproperty_qr_pullback(Ȳ)
        C = Tangent{T}
        ∂F = if x === :Q
            C(Q=Ȳ,)
        elseif x === :R
            C(R=Ȳ,)
        end
        return NoTangent(), ∂F, NoTangent()
    end
    return getproperty(F, x), getproperty_qr_pullback
end

function qr_rev_fullrank(q, r, dq, dr)
    dqnot0 = !(dq isa NoTangent)
    drnot0 = !(dr isa NoTangent)
    if (!dqnot0 && !drnot0)
        return NoTangent()
    end
    ex = drnot0 && dqnot0 ? r*dr' - dq'*q : (dqnot0 ? -dq'*q  : r*dr')
    b = dqnot0 ? q*copyltu!(ex)+dq : q*copyltu!(ex)
    return Matrix(trtrs!('U', 'N', 'N', r, do_adjoint(b))')
end

function qr_rev(q, r, dq, dr, A)
    dqnot0 = !(dq isa NoTangent)
    drnot0 = !(dr isa NoTangent)
    (!dqnot0 && !drnot0) && return NoTangent()
    size(r, 1) == size(r, 2) && return qr_rev_fullrank(q, r, dq ,dr)
    M, N = size(r)
    B = view(A,:,M+1:N)
    U = view(r,:,1:M)
    D = view(r,:,M+1:N)
    if drnot0
        dD = view(dr,:,M+1:N);
        da = qr_rev_fullrank(q, U, dqnot0 ? dq+B*dD' : B*dD', view(dr,:,1:M))
        db = q*dD
    else
        da = qr_rev_fullrank(q, U, dq, nothing);
        db = zero(B)
    end
    return hcat(da, db)
end


text:
"""
    get_author(pdf::String) -> String

Extract the author from the PDF.

# Arguments
- `pdf::String`: $(PDF_ARGUMENT)

# Returns
- `String`: the author

# Throws
- Exceptions from: [`get_author_title`](@ref)

# Example
```jldoctest; output = false
using PDFHighlights

path_to_pdf = joinpath(
    pathof(PDFHighlights) |> dirname |> dirname,
    "test",
    "pdf",
    "TestPDF.pdf"
)

get_author(path_to_pdf) == "Pavel Sobolev"

# output

true
```
"""
function get_author(pdf::String)::String
    author, _ = get_author_title(pdf)
    return author
end


text:
module FranklinTemplates

export newsite, addtemplate, modify

const FRANKLIN_PATH     = dirname(pathof(FranklinTemplates))
const TEMPL_PATH        = joinpath(FRANKLIN_PATH, "templates")
const LIST_OF_TEMPLATES = (
    "sandbox",
    "sandbox-extended",
    "basic",
    "jemdoc",
    "just-the-docs",
    "hyde",
    "hypertext",
    "lanyon",
    "minimal-mistakes",
    "pure-sm",
    "tufte",
    "vela")

include("utils.jl")

end # module


text:
# Author: Ashwani Rathee
# Program to Find Largest Number Among Three Numbers
println("Enter 3 Numbers to find the largest : ")
first=readline()
first=parse(Int64,first)

second=readline()
second=parse(Int64,second)

third=readline()
third=parse(Int64,third)

if first >second && first>third
    println("First Number is the largest: ",first)
elseif second>first && second>third
    println("Second number is the largest: ",second)
else
    println("Third number is the largest: ",third)
end

max_stars_repo_name: ashwani-rathee/Learning--Julia

max_issues_repo_name: ashwani-rathee/Learning--Julia

max_forks_repo_name: ashwani-rathee/Learning--Julia


text:
using PyDSTool, DataStructures

name = "Calcium channel model"
pars = Dict{String,Any}(
          "vl"=> -60,
          "vca"=> 120,
          "i" => 0,
          "gl" => 2,
          "gca" => 4,
          "c" => 20,
          "v1" => -1.2,
          "v2" => 18)
vars = Dict{String,Any}(
          "v" => "( i + gl * (vl - v) - gca * 0.5 * (1 + tanh( (v-v1)/v2 )) * (v-vca) )/c",
          "w" => "v-w")
ics = Dict{String,Any}(
          "v" => 0,
          "w" => 0)
tdomain = [0;30]

dsargs = build_ode(name,ics,pars,vars,tdomain)

#Solve the ODE
#d = solve_ode(dsargs)
#using Plots
#plot(d[:t],d[:v])

#Bifurcation Plots
ode = ds[:Generator][:Vode_ODEsystem](dsargs)
ode[:set](pars = Dict("i"=>-220))
ode[:set](ics  = Dict("v"=>-170))
PC = ds[:ContClass](ode)

bif = bifurcation_curve(PC,"EP-C",["i"],
                          max_num_points=450,
                          max_stepsize=2,min_stepsize=1e-5,
                          stepsize=2e-2,loc_bif_points="all",
                          save_eigen=true,name="EQ1",
                          print_info=true,calc_stab=true)

@test length(bif.changes) == 2
#=
using Plots
plot(bif,(:i,:v))
=#

bif = bifurcation_curve(PC,"LP-C",["i","gca"],
                        max_num_points=200,initpoint="EQ1:LP2",
                        max_stepsize=2,min_stepsize=1e-5,
                        stepsize=2e-2,loc_bif_points="CP",
                        save_eigen=true,name="SN1",
                        print_info=true,calc_stab=true,
                        solver_sequence=[:forward,:backward])

@test length(bif.changes) == 0

#=
using Plots
plot(bif,(:i,:gca))
=#


text:
export real_matrix, complex_matrix

"""
`real_matrix(x::SimpleQuaternion)` returns a 4-by-4 real matrix that represents
`x`.
"""
function real_matrix(x::SimpleQuaternion)
    a, b, c, d = get_parts(x)
    M = [
        a -b -c -d
        b a -d c
        c d a -b
        d -c b a
    ]
    return M
end

"""
`complex_matrix(x::SimpleQuaternion)` returns a 2-by-2 complex matrix
that represents `x`.
"""
function complex_matrix(x::SimpleQuaternion)
    a, b, c, d = get_parts(x)
    M = [
         a+b*im     c+d*im
        -c+d*im     a-b*im
    ]
    return M
end

function SimpleQuaternion(A::AbstractMatrix{T}) where T
    r,c = size(A)
    if r!=c || !in(r,[2,4])
        error("Matrix must be 2-by-2 or 4-by-4; got $r-by-$c.")
    end

    if r==4
        v = A[:,1]
        x = SimpleQuaternion(v...)
        if real_matrix(x) != A
            error("Matrix not derived from a quaternion")
        end
        return x
    end

    a,b = reim(A[1,1])
    c,d = reim(A[1,2])
    x = SimpleQuaternion(a,b,c,d)
    if complex_matrix(x) != A
        error("Matrix not derived from a quaternion")
    end
    return x
end


text:
"Rotates vector counter-clockwise with respecto to the x-axis."
rotx(θ::Real) = [1		0		0
				 0	cos(θ)	-sin(θ);
				 0	sin(θ)	cos(θ)]
"Rotates vector counter-clockwise with respecto to the y-axis."
roty(θ::Real) = [cos(θ) 0	sin(θ);
				 0		1		0;
				-sin(θ)	0 	cos(θ)]
"Rotates vector counter-clockwise with respecto to the z-axis."
rotz(θ::Real) = [cos(θ) -sin(θ)	0;
				 sin(θ) cos(θ)	0;
				 0 		0 		1]
#####################
## Gradient OBJECT ##
#####################
# @everywhere begin
"""
    Grad(A,T)

Square gradient object.

# Arguments
- `A::Real` := Gradient amplitude [T].
- `T::Real` := Gradient duration [s].

# Examples
```julia-repl
julia> Grad(1,1)
Grad(1, 1)
```
"""
mutable struct Grad
	A::Real #Amplitud [T]
	T::Real #Duration of sequence [s]
    function Grad(A,T)
		T < 0 ? error("Gradient duration must be positive.") : new(A, T)
    end
end
#Gradient operations
*(x::Grad,α::Real) = Grad(α*x.A,x.T)
*(α::Real,x::Grad) = Grad(α*x.A,x.T)
*(A::Matrix{Float64},GR::Matrix{Grad}) = begin
	N, M = size(GR)
	[sum(A[i,1:N] .* GR[:,j]) for i=1:N, j=1:M]
end
import Base.zero
zero(::Grad) = Grad(0,0)
/(x::Grad,α::Real) = Grad(x.A/α,x.T)
+(x::Grad,y::Grad) = Grad(x.A+y.A,max(x.T,y.T))
+(x::Array{Grad,1}, y::Array{Grad,1}) = [x[i]+y[i] for i=1:length(x)]
-(x::Grad) = -1*x
-(x::Grad,y::Grad) = Grad(x.A-y.A,max(x.T,y.T))
#Gradient functions
vcat(x::Array{Grad,1},y::Array{Grad,1}) = [i==1 ? x[j] : y[j] for i=1:2,j=1:length(x)]
vcat(x::Array{Grad,1},y::Array{Grad,1},z::Array{Grad,1}) = [i==1 ? x[j] : i==2 ? y[j] : z[j] for i=1:3,j=1:length(x)]
getproperty(x::Vector{Grad}, f::Symbol) = getproperty.(x,f)
getproperty(x::Matrix{Grad}, f::Symbol) = begin
	if f == :x
		x[1,:].A
	elseif f == :y
		x[2,:].A
	elseif f == :z && size(x,1) == 3
		x[3,:].A
	elseif f == :T
		x[1,:].T
	elseif f == :A
		getproperty.(x,:A)
	end
end

"""
	delay(T)

Auxiliary function, creates Grad object of duration `T` [s] and strength 0 [T].
"""
delay(T::Real) = Grad(0,T)
"""
	Grad_fun(f,T,N)

Generates **arbitrary gradient waveform** defined by function `f` in the interval t∈[0,`T`].
It uses `N` square gradients uniformly spaced in the interval.

# Arguments
 - `f::Function` := Gradient waveform.
 - `T::Real`     := Gradient duration.
 - `N::Integer=1`:= Number of sample points.

# Examples
```julia-repl
julia> Grad_fun(x-> sin(π*x),1,4)
1×4 Array{Grad,2}:
 Grad(0.0, 0.333333)  Grad(0.866025, 0.333333)  Grad(0.866025, 0.333333)  Grad(1.22465e-16, 0.333333)
```
"""
Grad(f::Function,T::Real,N::Int64=300) = begin
	Grads = [Grad(f(t),T/N) for t = range(0,stop=T,length=N)]
	reshape(Grads,(1,N))
end

Grad(A::Vector, T::Real) = begin
	N = length(A)
	Δt = T/N
	Grads = [Grad(A[i],Δt) for i = 1:N]
	reshape(Grads,(1,N))
end


text:
const BASE = "https://statsapi.mlb.com/api/"

"""
Endpoints' information that `MLBStatsAPI.jl` supports.
Users can reach this information by `meta()` or `note()` function.
"""
const ENDPOINTS = Dict(
    # game
    "game" => Dict(
        "url" => BASE * "version/game/gamePk/feed/live",
        "version" => "v1.1",
        "params" => ["timecode", "hydrate", "fields"],
        "required" => ["gamePk"],
        "pathparams" => ["version", "gamePk"],
    ),
    # schedule
    "schedule" => Dict(
        "url" => BASE * "version/schedule",
        "version" => "v1",
        "params" => [
            "scheduleType",
            "eventTypes",
            "hydrate",
            "teamId",
            "leagueId",
            "sportId",
            "gamePk",
            "gamePks",
            "venueIds",
            "gameTypes",
            "date",
            "startDate",
            "endDate",
            "opponentId",
            "fields",
        ],
        "required" => ["sportId", "gamePk", "gamePks"],
        "pathparams" => ["version"],
    ),
    # attendance
    "attendance" => Dict(
        "url" => BASE * "version/attendance",
        "version" => "v1",
        "params" => [
            "teamId",
            "leagueId",
            "season",
            "date",
            "leagueListId",
            "gameType",
            "fields",
        ],
        "required" => ["teamId", "leagueId", "leagueListId"],
        "pathparams" => ["version"],
    ),
    # awards
    "awards" => Dict(
        "url" => BASE * "version/awards/awardId/recipients",
        "version" => "v1",
        "params" => ["sportId", "leagueId", "season", "hydrate", "fields"],
        "required" => ["awardId"],
        "pathparams" => ["version", "awardId"],
    ),
    # draft
    "draft" => Dict(
        "url" => BASE * "version/draft/prospects/year",
        "version" => "v1",
        "params" => [
            "limit",
            "fields",
            "round",
            "name",
            "school",
            "state",
            "country",
            "position",
            "teamId",
            "playerId",
            "bisPlayerId",
        ],
        "required" => ["year"],
        "pathparams" => ["version", "year"],
    ),
    # people
    "people" => Dict(
        "url" => BASE * "version/people",
        "version" => "v1",
        "params" => ["personIds", "hydrate", "fields"],
        "required" => ["personIds"],
        "pathparams" => ["version"],
    ),
    # standings
    "standings" => Dict(
        "url" => BASE * "version/standings",
        "version" => "v1",
        "params" => [
            "leagueId",
            "season",
            "standingsTypes",
            "date",
            "hydrate",
            "fields",
        ],
        "required" => ["leagueId"],
        "pathparams" => ["version"],
    ),
    # teams
    "teams" => Dict(
        "url" => BASE * "version/teams",
        "version" => "v1",
        "params" => [
            "season",
            "activeStatus",
            "leagueIds",
            "sportIds",
            "gameType",
            "hydrate",
            "fields",
        ],
        "required" => [],
        "pathparams" => ["version"],
    ),
    # meta
    "meta" => Dict(
        "url" => BASE * "version/type",
        "version" => "v1",
        "params" => [],
        "required" => ["type"],
        "pathparams" => ["version", "type"],
    )
)


text:
module SeasonalStreamflowForecasts

import Base.show
import HTTP
using CSV: File
using JSON
using DataFrames: DataFrame, select!
using Dates: Date, year, month
using Printf: @sprintf

const SSF_URL = "http://www.bom.gov.au/water/ssf/"
const SITES_URL = SSF_URL * "content/images/forecast_site_geojson.json"

const SITE_PROPERTIES = ["name", "ID", "AWRC", "description",
				         "drainage", "basin", "area", "areaUnits"]

const HEADER_DELIM = "#"

export SSF, get_forecasts, close!

"""
    ssf = SSF()

Open the connection and download the information of SSF service sites e.g. name, ID, AWRC ID and description.

# Fields
* `sites`: Site information table

# Examples
```julia
julia> ssf = SSF();
julia> ssf.sites
215×8 DataFrame. Omitted printing of 5 columns
│ Row │ name                   │ ID              │ AWRC     │
│     │ String                 │ String          │ String   │
├─────┼────────────────────────┼─────────────────┼──────────┤
│ 1   │ upstreamofbaileysgrave │ G9070142        │ G9070142 │
│ 2   │ coenracecourse         │ 922101B         │ 922101B  │
│ 3   │ monument               │ 927001B         │ 927001B  │
...
```
"""
struct SSF
	sites::DataFrame

	function SSF()
		sites = get_sites()
		new(sites)
	end
end

"""
    data, header = get_forecasts(ssf::SSF, site_id::AbstractString, fc_date::Date)

Return seasonal forecasts of a site.

# Arguments
* `ssf` : SSF object
* `site_id`: AWRC ID of the site. The ID can found in the table from `get_sites()`
* `fc_date`: Forecast date

# Examples
```julia
julia> using Dates
julia> data, header = get_forecasts(ssf,"410730", Date(2020,8,1));
julia> data
5000×7 DataFrame. Omitted printing of 4 columns
│ Row  │ Member No. │  Streamflow Forecast (GL) Aug │  Streamflow Forecast (GL) Aug - Sep │
│      │ Int64      │ Float64                       │ Float64                             │
├──────┼────────────┼───────────────────────────────┼─────────────────────────────────────┤
│ 1    │ 1          │ 6.801                         │ 10.187                              │
│ 2    │ 2          │ 5.094                         │ 11.032                              │
│ 3    │ 3          │ 7.776                         │ 20.112                              │
...
julia> println(header)

Australian Bureau of Meteorology
Seasonal Streamflow Forecasts

Forecast data
...
```
"""
function get_forecasts(ssf::SSF,
                       site_id::AbstractString,
                       fc_date::Date)::Tuple{DataFrame,String}
    row = ssf.sites[ssf.sites.ID .== site_id, :]
    isempty(row) && return (DataFrame(), "")

    data_url = get_url(site_id, fc_date, row.drainage[1], row.basin[1])

	r = HTTP.get(data_url)

	body_buf = IOBuffer(String(r.body))

	header = extract_header!(body_buf, HEADER_DELIM)
	new_header = prune_header(header, HEADER_DELIM)

	body_buf = seek(body_buf, 0)
	data = DataFrame(File(body_buf, comment=HEADER_DELIM))

	return data, new_header
end

function show(io::IO, ssf::SSF)
    isempty(ssf.sites) && return
    println("The SSF Sites:")
    show(io, first(ssf.sites, 6))
    println(io)
    println(io, "...")
    return
end

"""
    close!(ssf::SSF)

Close the connection and reset the SSF object.
"""
function close!(ssf::SSF)
    select!(ssf.sites, Int[])
    return
end

"""
    sites, header = get_sites()

Download site information e.g. ID, AWRC ID and descriptions.
"""
function get_sites()::DataFrame
	r = HTTP.get(SITES_URL)

	features = JSON.parse(String(r.body))["stations"]["features"]

	sites = DataFrame()
	for key in SITE_PROPERTIES
		sites[!, key] = [x["properties"][key] for x in features]
	end

	return sites
end

"""
    url = get_url(site_id::AbstractString, fc_date::Date,
                  drainage::AbstractString, basin::AbstractString) => AbstractString

Return the URL to download data of a site.
"""
function get_url(site_id::AbstractString, fc_date::Date,
                 drainage::AbstractString, basin::AbstractString)::AbstractString
	fc_year = year(fc_date)
	fc_month = @sprintf("%02d", month(fc_date))

	url = SSF_URL * "content/$drainage/$basin/fc/$fc_year/$fc_month/$(site_id)_FC_10_$(fc_year)_$(fc_month)_table.csv"

    return url
end

"""
    header = extract_header!(body_buf::Base.GenericIOBuffer{Array{UInt8,1}}, delim::AbstractString)

Extract the header document of the data file. Note that it moves the position of body_buf.
"""
function extract_header!(body_buf::Base.GenericIOBuffer{Array{UInt8,1}},
                         delim::AbstractString)::Array{String,1}
    header = String[]
    for line in eachline(body_buf)
        startswith(line, delim) ? push!(header, line) : break
    end
    return header
end

"""
    new_header = prune_header(header::AbstractString, delim::AbstractString)

Prune the header document. It drop prefixes from the raw document.
"""
function prune_header(header::Array{String,1}, delim::AbstractString)::String
    new_header = String[]
    for line in header
        #TODO: Give more information about the format error.
        startswith(line, delim) || throw(IOError("A header line does not start with $(delim)."))
        the_line = rstrip(line[3:end])
        push!(new_header, the_line)
    end
    return join(new_header, "\n")
end

end # module


text:
using .NamedDims: NamedDimsArray
dimname(a::NamedDimsArray{N},i) where N = N[i]
dimnames(a::NamedDimsArray{N}) where N = N
getdata(a::NamedDimsArray) = parent(a)
function yaxcreate(::Type{<:NamedDimsArray},data, dnames, dvals, atts)
  n = ntuple(i->dnames[i],ndims(data))
  NamedDimsArray(data,n)
end


text:
module Imogen

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 3
end

using DataFrames, Distances, NearestNeighbors, SpecialFunctions, Statistics, Reexport

@reexport using Random

export InfoDist, observe!, clear!, estimate
export Entropy, entropy!, entropy
export MutualInfo, mutualinfo!, mutualinfo
export ActiveInfo, activeinfo!, activeinfo
export SpecificInfo, specificinfo!, specificinfo
export TransferEntropy, transferentropy!, transferentropy
export AbstractVertex, AbstractUnnamedVertex, AbstractNamedVertex, id, name, payload, above, below
export UnnamedVertex, Vertex, clone
export Hasse, top, bottom, vertices, zero!, prune, graphviz
export WilliamsBeer, pid, pid!
export Significance, EmpiricalSig, AnalyticSig, sig, @sig
export histories, box, encodehistories

include("core.jl")
include("methods.jl")

include("util.jl")

include("entropy.jl")

include("mi.jl")

include("ai.jl")

include("si.jl")

include("te.jl")

include("lattice.jl")
include("pid.jl")
include("williamsbeer.jl")

include("sig.jl")

end


text:
using BenchmarkTools
using DataFrames
using CSV
using Statistics
using StatsBase: sample, quantile
using Plots

using Revise
using EvoTrees
using EvoTrees: sigmoid, logit
using Profile
using ProfileView

# prepare a dataset
features = rand(10_000) .* 5
X = reshape(features, (size(features)[1], 1))
Y = sin.(features) .* 0.5 .+ 0.5
Y = logit(Y) + randn(size(Y))
Y = sigmoid(Y)
𝑖 = collect(1:size(X,1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace = false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1))+1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

# linear
params1 = EvoTreeRegressor(
    loss=:linear, metric=:mse,
    nrounds=100, nbins = 100,
    λ = 0.5, γ=0.1, η=0.05,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)

@time model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 25)
# 67.159 ms (77252 allocations: 28.06 MiB)
@btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval)
# Profile.clear()  # in case we have any previous profiling data
# @profile grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 25)
# ProfileView.view()

# @btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval, print_every_n = 25, metric=:mae)
@time pred_train_linear = predict(model, X_train)
@time pred_eval_linear = predict(model, X_eval)
mean(abs.(pred_train_linear .- Y_train))
sqrt(mean((pred_train_linear .- Y_train) .^ 2))

# logistic / cross-entropy
params1 = EvoTreeRegressor(
    loss=:logistic, metric = :logloss,
    nrounds=100, nbins = 100,
    λ = 0.5, γ=0.1, η=0.05,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)

@time model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 25)
# 218.040 ms (123372 allocations: 34.71 MiB)
# @btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval)
@time pred_train_logistic = predict(model, X_train)
@time pred_eval_logistic = predict(model, X_eval)
sqrt(mean((pred_train_logistic .- Y_train) .^ 2))

# Poisson
params1 = EvoTreeRegressor(
    loss=:poisson, metric = :poisson,
    nrounds=100, nbins = 100,
    λ = 0.5, γ=0.1, η=0.05,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)
@time model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 25)
# @btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval)
@time pred_train_poisson = predict(model, X_train)
@time pred_eval_poisson = predict(model, X_eval)
sqrt(mean((pred_train_poisson .- Y_train) .^ 2))

# L1
params1 = EvoTreeRegressor(
    loss=:L1, α=0.5, metric = :mae,
    nrounds=100, nbins=100,
    λ = 0.5, γ=0.1, η=0.05,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)
@time model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 25)
@time pred_train_L1 = predict(model, X_train)
@time pred_eval_L1 = predict(model, X_eval)
sqrt(mean((pred_train_L1 .- Y_train) .^ 2))

x_perm = sortperm(X_train[:,1])
plot(X_train, Y_train, ms = 1, mcolor = "gray", mscolor = "lightgray", background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("feature"), yaxis = ("target"), legend = true, label = "")
plot!(X_train[:,1][x_perm], pred_train_linear[x_perm], color = "navy", linewidth = 1.5, label = "Linear")
plot!(X_train[:,1][x_perm], pred_train_logistic[x_perm], color = "darkred", linewidth = 1.5, label = "Logistic")
plot!(X_train[:,1][x_perm], pred_train_poisson[x_perm], color = "green", linewidth = 1.5, label = "Poisson")
plot!(X_train[:,1][x_perm], pred_train_L1[x_perm], color = "gold", linewidth = 1.5, label = "L1")
savefig("regression_sinus.png")

###############################
## Quantiles
###############################
# q50
params1 = EvoTreeRegressor(
    loss=:quantile, α=0.5, metric=:quantile,
    nrounds=100, nbins = 100,
    λ = 0.1, γ=1e-3, η=0.05,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)

@time model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 25)
# 116.822 ms (74496 allocations: 36.41 MiB) for 100 iterations
# @btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval)
@time pred_train_q50 = predict(model, X_train)
sum(pred_train_q50 .< Y_train) / length(Y_train)

# q20
params1 = EvoTreeRegressor(
    loss=:quantile, α=0.2, metric=:quantile,
    nrounds=200, nbins = 100,
    λ = 0.1, γ=1e-3, η=0.05,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)
@time model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 25)
@time pred_train_q20 = predict(model, X_train)
sum(pred_train_q20 .< Y_train) / length(Y_train)

# q80
params1 = EvoTreeRegressor(
    loss=:quantile, α=0.8, metric=:quantile,
    nrounds=200, nbins = 100,
    λ = 0.1, γ=1e-3, η=0.05,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)
@time model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 25)
@time pred_train_q80 = predict(model, X_train)
sum(pred_train_q80 .< Y_train) / length(Y_train)

x_perm = sortperm(X_train[:,1])
plot(X_train, Y_train, ms = 1, mcolor = "gray", mscolor = "lightgray", background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("feature"), yaxis = ("target"), legend = true, label = "")
plot!(X_train[:,1][x_perm], pred_train_q50[x_perm], color = "navy", linewidth = 1.5, label = "Median")
plot!(X_train[:,1][x_perm], pred_train_q20[x_perm], color = "darkred", linewidth = 1.5, label = "Q20")
plot!(X_train[:,1][x_perm], pred_train_q80[x_perm], color = "green", linewidth = 1.5, label = "Q80")
savefig("quantiles_sinus.png")


text:
using Reactive
using Interact

export trim_viz, function_viz, fzeros_viz,
       bisection_viz, limit_viz,
       derivative_viz, critical_pts_viz,
       newton_viz, riemann_viz

"""

Some visualizations in the `MTH229` package:

- `trim_viz(f, a, b). Shows a graph of `f` with controls to trim large negative and positive `y` values.

- `functions_viz(f, a, b)`. Shows where a function is positive, increasing, or concave up over `[a,b]`.

- `fzeros_viz(f, a, b)`. Shows zeros of a function over an interval using `fzeros`

- `bisection_viz(f, a , b)`. For a bracketing interval, shows first few bisection steps.

- `limit_viz(f, c, dir="+")`. Shows limiting value as `x` approaches `c`.

- `derivative_viz(f, c)`. Shows the secant line and its approach to the tangent line for different values of "`h`."

- `critical_pts_viz(f, a, b)`. Shows zeros of the derivative found with `fzeros`. (So *may* show some places where derivative is undefined too.)

- `newton_viz(f, x0). A simple function to visualize some `steps` of Newton's method. The values of `a` and `b` are optional, but, if not set, the x-viewing window will be determined by the points of the sequence of approximations

- `riemann_viz(f, a, b)`. Show the Riemann sum partition in action.

"""
visualizations = nothing
export visualizations

"""

`trim_viz(f::Function, a, b)`: plot `f` over `[a,b]` with trimming of large and small values

"""
function trim_viz(f::Function, a=-5, b=5)
    trim(f, lo, hi) = x -> f(x) < lo ? NaN : (f(x) > hi ? NaN : f(x))
    @manipulate for lo in -[2^i for i in 1:10], hi in [2^i for i in 1:10]
        plot(trim(f, lo, hi), a, b)
    end
end

"""

Simple illustration of how to find positive, increasing or concavity graphically. Uses `plotif` helper.

```
function_viz(sin, a, b)
```
"""
function function_viz(f::Function, a=-5, b=5)
    @manipulate for u in ["positive", "increasing", "concave up"]
        if u == "positive"
            p = plotif(f, f, a, b, title="Highlighting where f > 0")
        elseif u == "increasing"
            p = plotif(f, D(f), a, b, title="Showng where f' > 0")
        else
            p = plotif(f, D(f,2), a, b, title="Showing where f'' > 0")
        end
        p
    end
end

"""

fzeros_viz: show zeros of a function over `[a,b]`. Uses `fzeros`, so may miss some.

```
fzeros_viz(f, a, b)
```
"""
function fzeros_viz(f::Function, a, b)
    zs = fzeros(f, a, b)
    plot(f, a, b, legend=false)
    scatter!(zs, 0*zs)
end

"""

`bisection_viz(f, a, b)`: for a bracketing interval `[a,b]`, show first few steps of the bisection method

```
bisection_viz(sin, 3, 4)
```
"""
function bisection_viz(f, a, b)
    lo,hi = a, b
    rt = fzero(f, lo, hi)

    @manipulate for n=slider(1:8, value=1, label="n")
        cvged = false
        a,b = lo, hi
        for i in 1:n
            c = (a + b)/2
            for j in 2:i
                if f(c) == 0
                    cvged = true
                    break
                elseif f(a)*f(c) < 0
                    b = c
                else
                    a = c
                end
            end
        end

        p = plot(f, lo, hi, linewidth=3, legend=false)
        scatter!(p, [rt], [0], markersize=3)
        if cvged
            scatter!(p, [rt], [0], markersize=6)
        else
            plot!(p, [a, (a + b)/2], [0,0], linewidth=3)
            plot!(p, [(a + b)/2, b], [0,0], linewidth=3)
        end
        p
    end
end

"""
`limit_viz(f, c)`. Show graph of `f` and how `x` approaches `c`

```
f(x) = sin(x) / x
limit_viz(f, 0)
```
"""
function limit_viz(f::Function, c, dir="+")
    L = N(limit(f, c))
    @manipulate for i in slider(0:5, value=0, label="n")
        h = (dir == "+" ? (1/2)^i : -(1/2)^i)
        s = (dir == "+" ? "+" : "-")
        val = f(c + h)
        plot(f, c-1, c + 1, linewidth=3, title="At c $s (1/2)^$i f is $(round(val,3))", legend=false)
        plot!([c-1, c+1], L * [1,1])
        plot!([c-1, c+1], val * [1,1])
        scatter!([c+h], [f(c+h)], markersize=3)
    end
end

"""
Simple illustration of how secant line approaches tangent line as `h` goes to `0`. This plots `f` from `[c-1/2, c + 1/2]` and draws a secant line for different values of `h`.

Example:
```
derivative_viz(sin, pi/4)
```
"""
function derivative_viz(f::Function, c, a=c-1/2, b=c+1/2)

    @manipulate for j in slider(1:8, value=1, label="h=(1/2)^i")
	h = (1/2)^j
        m = (f(c + h) - f(c)) / h
        err = abs(f'(c) - m)
        e = Int(round(log2(1 / err), 0))

	p = plot(f, a, b, title="Plotting secant line with h=(1/2)^$j, error is about (1/2)^(-$e)", linewidth=3, legend=false)
	m = (f(c + h) - f(c)) / h
	plot!(p, x -> f(c) + m * (x-c), a, b, linewidth=3)
        scatter!(p, [c, c+h], [f(c), f(c+h)], markersize=3)
#	plot!(p, x -> f(c) + f'(c)*(x-c))
        p
    end

end

"""

critical_pts_viz: show critical points of a function over `[a,b]`. Uses `fzeros`, so may miss some.

```
critical_pts_viz(f, a, b)
```
"""
function critical_pts_viz(f::Function, a, b)
    zs = fzeros(f', a, b)
    plot(f, a, b, legend=false)
    scatter!(zs, f.(zs))
end

"""
newton_viz

Simple visualization of Newton's method

```
f(x) = log(x) - 1/10
newton_viz(f, 1) # newton_viz(f, 1, steps=20)
"""
function newton_viz(f::Function, x0, a=Inf, b = -Inf; steps=10)
    @manipulate for nsteps = slider(1:steps, value=1, label="no. steps")
        newton_vis(f, x0, a, b; steps=nsteps)
    end
end

"""

Riemann visualization.

Show right-riemann approximation for different values of `n`.

```
riemann_viz(sin, 0, pi)
```

"""
function riemann_viz(f::Function, a, b)
    @manipulate for j in slider(1:8, value=1, label="no. rectangles")
        n = 2^j
        val = riemann(f, a, b, n)
        act_val = quadgk(f, a, b)[1]
        err = abs(val - act_val)
        e = Int(round(log2(1 / abs(val - act_val)), 0))
        p = plot(f, a, b, linewidth=3, legend=false, title="Riemann sum with n=2^$j: $(round(val,3)). Error is about 2^(-$e)")
        delta = (b-a)/n
        for i in 1:n
            plot!(p, a + delta * [i-1, i, i, i-1, i-1], f(a + i * delta) * [0,0,1,1,0], color=:blue)
        end
        p
    end
end


text:
using LinearAlgebra

include("./kernel.jl")
include("./transform.jl")
include("./incremental.jl")
include("./kernel_system.jl")

@doc """
"""
struct BTG{G<:AbstractTransform,K<:AbstractCorrelation,C}
    horizon::Int
    n::Int
    g::G
    k::K
    comp::ComputeData
    comp_k::Vector{C}
end
function BTG()
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
include("../src/pressurecorrelations.jl")

@testset "Friction factors" begin
    @test ChenFrictionFactor(35700, 2.259, 0.001)/4 ≈ 0.0063 atol = 0.001

    @test ChenFrictionFactor(5.42e5, 2.295, 0.0006)/4 ≈ 0.0046 atol=0.001 # per economides
    @test SerghideFrictionFactor(5.42e5, 2.295, 0.0006)/4 ≈ 0.0046 atol=0.001

    @test ChenFrictionFactor(71620, 2.441, 0.0002) ≈ 0.021 atol = 0.005 #per Takacs
    @test SerghideFrictionFactor(71620, 2.441, 0.0002) ≈ 0.021 atol = 0.005

    @test ChenFrictionFactor(101000, 2.295, 0.0002)/4 ≈ 0.005 atol = 0.005 #per economides
    @test SerghideFrictionFactor(101000, 2.295, 0.0002)/4 ≈ 0.005 atol = 0.005
end

@testset "Superficial velocities" begin

    @test liquidvelocity_superficial(375, 0, 2.441, 1.05, 0) ≈ 0.79 atol = 0.01
    @test gasvelocity_superficial(375, 0, 480, 109.7, 2.441, 0.0389) ≈ 1.93 atol = 0.01

end

#%% Beggs and Brill flowmap

#= B&B flowmap graphical test, to make sure the correct mapping is generated
# Re-run after any modifications to B&B flowmap function.
# Note that this does not establish stability at the region boundaries.

λ_l = 10 .^ collect(-4:0.1:0);
N_Fr = 10 .^ collect(-1:0.1:3);

function crossjoin(x, y)

    output = Array{Float64, 2}(undef, length(x) * length(y), 2)

    index = 1
    @inbounds for i in x
        @inbounds for j in y
            output[index, 1] = i
            output[index, 2] = j
            index += 1
        end
    end

    return output
end

testdata = crossjoin(λ_l, N_Fr);
flowpattern = map(BeggsAndBrillFlowMap, testdata[:,1], testdata[:,2]);

using Gadfly

plot(x = testdata[:,1], y = testdata[:,2], color = flowpattern,
Geom.rectbin,
Scale.x_log10, Scale.y_log10, Coord.Cartesian(xmin=-4, xmax=0, ymin=-1, ymax=3))
=#

@testset "Beggs and Brill vertical - single segment" begin

    inclination = 0
    v_sl = 0.79
    v_sg = 1.93
    v_m = 2.72
    ρ_l = 49.9
    ρ_g = 1.79
    id = 2.441
    σ_l = 8
    μ_l = 3
    μ_g = 0.02
    roughness = 0.0006
    pressure_est = 346.6

    @test BeggsAndBrill(1, 1, inclination, id, v_sl, v_sg, ρ_l, ρ_g, σ_l, μ_l, μ_g, roughness, pressure_est, SerghideFrictionFactor, true, false) ≈ 0.170 atol = 0.01

end


text:
#Linear ARD Covariance Function

@doc """
# Description
Constructor for the ARD linear kernel (covariance)

k(x,x') = xᵀL⁻²x', where L = diag(ℓ₁,ℓ₂,...)
# Arguments:
* `ll::Vector{Float64}`: Log of the length scale ℓ
""" ->
type LinArd <: Kernel
    ll::Vector{Float64}      # Log of Length scale
    dim::Int                 # Number of hyperparameters
    LinArd(ll::Vector{Float64}) = new(ll,size(ll,1))
end

function cov(lin::LinArd, x::Vector{Float64}, y::Vector{Float64})
    ell = exp(lin.ll)
    K = dot(x./ell,y./ell)
    return K
end

function cov(lin::LinArd, X::Matrix{Float64}, data::EmptyData)
    ell = exp(lin.ll)
    return (X./ell)' * (X./ell)
end

get_params(lin::LinArd) = lin.ll
get_param_names(lin::LinArd) = get_param_names(lin.ll, :ll)
num_params(lin::LinArd) = lin.dim

function set_params!(lin::LinArd, hyp::Vector{Float64})
    length(hyp) == lin.dim || throw(ArgumentError("Linear ARD kernel has $(lin.dim) parameters"))
    lin.ll = hyp
end

function grad_kern(lin::LinArd, x::Vector{Float64}, y::Vector{Float64})
    ell = exp(lin.ll)
    dK_ell = -2.0*(x./ell).*(y./ell)
    return dK_ell
end

function grad_stack!(stack::AbstractArray, lin::LinArd, X::Matrix{Float64}, data::EmptyData)
    ell = exp(lin.ll)
    d, nobsv = size(X)
    for j in 1:d
        stack[:,:,j] = ((-2.0/ell[j]^2).*vec(X[j,:])) * vec(X[j,:])'
    end
end


text:
using JSON

# replace get_odata for testing
function CBSOData3.get_odata(url)
    if url == "http://opendata.cbs.nl/ODataFeed/odata/80137ned?\$format=json"
        content = JSON.parse(String(open("meta_1.json") do io read(io) end))
        return content
    elseif url == "https://opendata.cbs.nl/ODataFeed/odata/80137ned/TableInfos?\$format=json"
        content = JSON.parse(String(open("meta_2.json") do io read(io) end))
        return content
    elseif url == "https://opendata.cbs.nl/ODataFeed/odata/80137ned/CategoryGroups?\$format=json"
        content = JSON.parse(String(open("meta_3.json") do io read(io) end))
        return content
    elseif url == "https://opendata.cbs.nl/ODataFeed/odata/80137ned/Persoonskenmerken?\$format=json"
        content = JSON.parse(String(open("meta_4.json") do io read(io) end))
        return content
    elseif url == "https://opendata.cbs.nl/ODataFeed/odata/80137ned/BedrijfstakkenSBI93?\$format=json"
        content = JSON.parse(String(open("meta_5.json") do io read(io) end))
        return content
    elseif url == "https://opendata.cbs.nl/ODataFeed/odata/80137ned/Perioden?\$format=json"
        content = JSON.parse(String(open("meta_6.json") do io read(io) end))
        return content
    elseif url == "https://opendata.cbs.nl/ODataFeed/odata/80137ned/DataProperties?\$format=json"
        content = JSON.parse(String(open("meta_7.json") do io read(io) end))
        return content
    elseif url == "http://opendata.cbs.nl/ODataFeed/odata/80137ned/TypedDataSet/?\$format=json"
        content = JSON.parse(String(open("meta_8.json") do io read(io) end))
        return content
    elseif url == "http://opendata.cbs.nl/ODataCatalog/Tables?\$format=json"
        content = JSON.parse(String(open("tables_1.json") do io read(io) end))
        return content
    elseif url == "https://opendata.cbs.nl/ODataFeed/OData/83599NED/UntypedDataSet/?\$format=json"
        content = JSON.parse(String(open("tables_2.json") do io read(io) end))
        return content
    elseif url == "https://opendata.cbs.nl/ODataFeed/OData/83599NED/UntypedDataSet/?\$format=json&\$skip=10000"
        content = JSON.parse(String(open("tables_3.json") do io read(io) end))
        return content
    else
        println("Fout:" * url)
        return Dict("value"=>nothing)
    end
end


text:
# -----------------------------------------------------------------------------
# Implementation of Atomistic MolecularDynamicsResult API
# -----------------------------------------------------------------------------

"""
    MollyResult <: MolecularDynamicsResult

The result generated from running an `MollySimulator`.

**Field descriptions**
- `system::System` the `System` resulting from the simulation
- `simulator::MollySimulator` the `MollySimulator` that produced the result
"""
struct MollyResult <: MolecularDynamicsResult
    system::System
    simulator::MollySimulator
end

function get_time_range(result::MollyResult)
    s = result.simulator
    Δt = s.Δt * s.stride
    start = s.t₀ + Δt - s.Δt
    steps = s.steps ÷ s.stride
    start:Δt:start+Δt*steps
end
get_num_bodies(result::MollyResult) = length(result.system)
get_bounding_box(result::MollyResult) = bounding_box(result.system)
get_boundary_conditions(result::MollyResult) = boundary_conditions(result.system)

reference_temperature(result::MollyResult) = coupling_reference_temperature(result.simulator.coupling)
coupling_reference_temperature(coupling) = coupling.temperature
coupling_reference_temperature(::NoCoupling) = missing

get_positions(result::MollyResult) = position(result.system)
get_velocities(result::MollyResult) = velocity(result.system)
get_particles(result::MollyResult) = AtomsBase.Atom.(result.system.atoms_data, get_positions(result), get_velocities(result))
get_system(result::MollyResult) = result.system

get_positions(result::MollyResult, t::Integer) = result.system.loggers["c"].coords[t]
get_velocities(result::MollyResult, t::Integer) = result.system.loggers["v"].velocities[t]
get_particles(result::MollyResult, t::Integer) = AtomsBase.Atom.(result.system.atoms_data, get_positions(result, t), get_velocities(result, t))

temperature(result::MollyResult, t::Integer) = result.system.loggers["t"].temperatures[t]
kinetic_energy(result::MollyResult, t::Integer) = result.system.loggers["k"].energies[t]
potential_energy(result::MollyResult, t::Integer) = result.system.loggers["p"].energies[t]

"""
    visualize(result::MollyResult, filename::String; kwargs...)

Animate a `MollyResult` using Molly's `visualize` function and store the result in a file.
Any kwargs will be passed directly to the `visualize` function.

This function is only available when GLMakie is imported.
"""
Molly.visualize(result::MollyResult, filename::String; kwargs...) = visualize(result.system.loggers["c"], result.system.box_size, filename; kwargs...)


text:
module ModV5Server
export V5Server,run

using Sockets
using ProtoBuf
using ..V5RPC
import ..V5RPC:on_event,get_team_info
using ..PB
using ..ModV5Packet

struct V5Server
    strategy::Strategy
    endpoint::Tuple{IPAddr,Integer}
end # struct

V5Server(strategy::Strategy, port::Integer) = V5Server(strategy, (IPv4(0), port))

function recvpkt(sock::UDPSocket)
    (endpoint, data) = recvfrom(sock)
    io = IOBuffer(data)
    return (endpoint, read(io, V5Packet))
end

function sendpkt(sock::UDPSocket, pkt::V5Packet, endpoint::Sockets.InetAddr)
    io = IOBuffer()
    write(io, pkt)
    bytes = take!(io)
    send(sock, endpoint.host, endpoint.port, bytes)
end

function Base.run(srv::V5Server)
    sock = UDPSocket()
    (ip, port) = srv.endpoint
    bind(sock, ip, port) || error("Failed to bind socket on $ip:$port")
    lastid = nothing
    lastret = nothing
    try
        while true
            (ep, pkt) = recvpkt(sock)
            if pkt.requestid == lastid
                ret = lastret
            else
                ret = dispatch(srv, pkt)
                lastid = pkt.requestid
                lastret = ret
            end
            sendpkt(sock, ret, ep)
        end
    catch ex
        ex isa InterruptException || rethrow(ex)
    finally
        close(sock)
    end
end

function dispatch(srv::V5Server, pkt::V5Packet)::V5Packet
    io = PipeBuffer()
    write(io, pkt.payload)
    call = readproto(io, RPCCall())
    callsym = which_oneof(call, :method)
    callobj = get_field(call, callsym)
    params = fieldnames(typeof(callobj)) .|> n->getfield(callobj, n)
    fn = getfield(V5RPC, callsym)
    ret = fn(srv.strategy, params...)
    packed = pack(Val(callsym), ret)
    V5Packet(packed != nothing ? packed : UInt8[], pkt.requestid)
end

pack(data::ProtoType)::Vector{UInt8} = begin
    io = IOBuffer()
    writeproto(io, data)
    take!(io)
end

pack(::Val{:on_event}, data) = nothing

pack(::Val{:get_team_info}, data::TeamInfo) =
    pack(GetTeamInfoResult(team_info = data))
pack(m::Val{:get_team_info}, data::String) =
    pack(m, TeamInfo(team_name = data, version = Version.V1_1))

pack(::Val{:get_instruction}, data::Vector{Wheel}, info::ControlInfo) =
    pack(GetInstructionResult(wheels = data, command = info))
function pack(m::Val{:get_instruction}, data::Vector)
    pack(m, (data, 0))
end
function pack(m::Val{:get_instruction}, data::Tuple{Vector, Integer})
    wdata, idata = data
    @assert length(wdata) == 5
    wheels = map(wdata) do (lspd, rspd)
        Wheel(left_speed = lspd, right_speed = rspd)
    end
    info = ControlInfo(command = idata)
    pack(m, wheels, info)
end

pack(::Val{:get_placement}, data::Placement) = pack(GetPlacementResult(placement = data))
function pack(m::Val{:get_placement}, data::Vector)
    @assert length(data) == 5
    robots = map(data) do (x, y, r)
        Robot(position = Vector2(x = x, y = y), rotation = r, wheel = Wheel(left_speed = 0, right_speed = 0))
    end
    ball = Ball(position = Vector2(x = 0, y = 0))
    pack(m, Placement(robots = robots, ball = ball))
end

end # module


text:
module WebSockets

# This module implements the server side of the WebSockets protocol. Some
# things would need to be added to implement a WebSockets client, such as
# masking of sent frames.
#
# WebSockets expects to be used with HttpServer to provide the HttpServer
# for accepting the HTTP request that begins the opening handshake. WebSockets
# implements a subtype of the WebSocketInterface from HttpServer; this means
# that you can create a WebSocketsHandler and pass it into the constructor for
# an http server.
#
# Future improvements:
# 1. Logging of refused requests and closures due to bad behavior of client.
# 2. Better error handling (should we always be using "error"?)
# 3. Unit tests with an actual client -- to automatically test the examples.
# 4. Send close messages with status codes.
# 5. Allow users to receive control messages if they want to.

using HttpCommon
using HttpServer
using Codecs
using MbedTLS
using Compat; import Compat.String

export WebSocket,
       WebSocketHandler,
       write,
       read,
       close,
       send_ping,
       send_pong

const TCPSock = VERSION < v"0.4.0-dev" ? Base.TcpSocket : Base.TCPSocket
if VERSION < v"0.4.0-dev"
  init_socket(sock) = nothing
else
  init_socket(sock) = Base.buffer_writes(sock) # Buffer writes to socket till flush(sock)
end

# A WebSocket is a wrapper over a TcpSocket. It takes care of wrapping outgoing
# data in a frame and unwrapping (and concatenating) incoming data.
type WebSocket
  id::Int
  socket::TCPSock
  is_closed::Bool
  sent_close::Bool

  function WebSocket(id::Int,socket::TCPSock)
    init_socket(socket)
    new(id,socket, !isopen(socket), false)
  end
end

# WebSocket Frames
#
#      0                   1                   2                   3
#      0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1
#     +-+-+-+-+-------+-+-------------+-------------------------------+
#     |F|R|R|R| opcode|M| Payload len |    Extended payload length    |
#     |I|S|S|S|  (4)  |A|     (7)     |             (16/64)           |
#     |N|V|V|V|       |S|             |   (if payload len==126/127)   |
#     | |1|2|3|       |K|             |                               |
#     +-+-+-+-+-------+-+-------------+ - - - - - - - - - - - - - - - +
#     |     Extended payload length continued, if payload len == 127  |
#     + - - - - - - - - - - - - - - - +-------------------------------+
#     |                               |Masking-key, if MASK set to 1  |
#     +-------------------------------+-------------------------------+
#     | Masking-key (continued)       |          Payload Data         |
#     +-------------------------------- - - - - - - - - - - - - - - - +
#     :                     Payload Data continued ...                :
#     + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +
#     |                     Payload Data continued ...                |
#     +---------------------------------------------------------------+
#

# Opcode values
#  *  %x0 denotes a continuation frame
const OPCODE_CONTINUATION = 0x00
#  *  %x1 denotes a text frame
const OPCODE_TEXT = 0x1
#  *  %x2 denotes a binary frame
const OPCODE_BINARY = 0x2
#  *  %x3-7 are reserved for further non-control frames
#
#  *  %x8 denotes a connection close
const OPCODE_CLOSE = 0x8
#  *  %x9 denotes a ping
const OPCODE_PING = 0x9
#  *  %xA denotes a pong
const OPCODE_PONG = 0xA
#  *  %xB-F are reserved for further control frames

# Write the raw frame to a bufffer
function write_fragment(io::IO, islast::Bool, data::Array{UInt8}, opcode)
  l = length(data)
  b1::UInt8 = (islast ? 0b1000_0000 : 0b0000_0000) | opcode

  # TODO: Do the mask xor thing??
  # 1. set bit 8 to 1,
  # 2. set a mask
  # 3. xor data with mask

  if l <= 125
    write(io, b1)
    write(io, @compat UInt8(l))
    write(io, data)
  elseif l <= typemax(UInt16)
    write(io, b1)
    write(io, @compat UInt8(126))
    write(io, hton(@compat UInt16(l)))
    write(io, data)
  elseif l <= typemax(UInt64)
    write(io, b1)
    write(io, @compat UInt8(127))
    write(io, hton(@compat UInt64(l)))
    write(io, data)
  else
    error("Attempted to send too much data for one websocket fragment\n")
  end
end

# A version of send_fragment for text data.
function write_fragment(io::IO, islast::Bool, data::String, opcode)
  write_fragment(io, islast, data.data, opcode)
end

if VERSION < v"0.4.0-dev"
  function locked_write(io::IO, islast::Bool, data, opcode)
    buf = IOBuffer()
    write_fragment(buf, islast, data, opcode)
    write(io, takebuf_array(buf))
  end
else
  function locked_write(io::IO, islast::Bool, data, opcode)
    isa(io, TCPSock) && lock(io.lock)
    try
      write_fragment(io, islast, data, opcode)
    finally
      if isa(io, TCPSock)
        flush(io)
        unlock(io.lock)
      end
    end
  end
end

# Write text data; will be sent as one frame.
function Base.write(ws::WebSocket,data::String)
  if ws.is_closed
    @show ws
    error("Attempted write to closed WebSocket\n")
  end
  locked_write(ws.socket, true, data, OPCODE_TEXT)
end

# Write binary data; will be sent as one frmae.
function Base.write(ws::WebSocket, data::Array{UInt8})
  if ws.is_closed
    @show ws
    error("attempt to write to closed WebSocket\n")
  end
  locked_write(ws.socket, true, data, OPCODE_BINARY)
end

# Send a ping message, optionally with data.
function write_ping(io::IO, data = "")
  locked_write(io, true, data, OPCODE_PING)
end

send_ping(ws, data...) = write_ping(ws.socket, data...)

# Send a pong message, optionally with data.
function write_pong(io::IO, data = "")
  locked_write(io, true, data, OPCODE_PONG)
end

send_pong(ws, data...) = write_pong(ws.socket, data...)

# Send a close message.
function Base.close(ws::WebSocket)
    # Tell client to close connection
    locked_write(ws.socket, true, "", OPCODE_CLOSE)
    ws.is_closed = true

    # Wait till client responds with an OPCODE_CLOSE
    while true
      wsf = read_frame(ws.socket)
      # ALERT: stuff might get lost in ether here
      is_control_frame(wsf) || continue
      wsf.opcode == OPCODE_CLOSE || continue
      break
    end
    close(ws.socket)
end

# A WebSocket is closed if the underlying TCP socket closes, or if we send or
# receive a close message.
Base.isopen(ws::WebSocket) = !ws.is_closed && isopen(ws.socket)

# Represents one (received) message frame.
type WebSocketFragment
  is_last::Bool
  rsv1::Bool
  rsv2::Bool
  rsv3::Bool
  opcode::UInt8  # This is actually a UInt4 value.
  is_masked::Bool
  payload_len::UInt64
  maskkey::Vector{UInt8}  # This will be 4 bytes on frames from the client.
  data::Vector{UInt8}  # For text messages, this is a String.
end

# This constructor handles conversions from bytes to bools.
function WebSocketFragment(
   fin::UInt8
  ,rsv1::UInt8
  ,rsv2::UInt8
  ,rsv3::UInt8
  ,opcode::UInt8
  ,masked::UInt8
  ,payload_len::UInt64
  ,maskkey::Vector{UInt8}
  ,data::Vector{UInt8})

  WebSocketFragment(
      fin != 0
    , rsv1 != 0
    , rsv2 != 0
    , rsv3 != 0
    , opcode
    , masked != 0
    , payload_len
    , maskkey
    , data)
end

# Control frames have opcodes with the highest bit = 1.
is_control_frame(msg::WebSocketFragment) = (msg.opcode & 0b0000_1000) > 0

# Respond to pings, ignore pongs, respond to close.
function handle_control_frame(ws::WebSocket,wsf::WebSocketFragment)
  if wsf.opcode == OPCODE_CLOSE
    # Reply with an empty CLOSE frame
    locked_write(ws.socket, true, "", OPCODE_CLOSE)
    ws.is_closed = true
    wait(ws.socket.closenotify)
  elseif wsf.opcode == OPCODE_PING
    write_pong(ws.socket,wsf.data)
  elseif wsf.opcode == OPCODE_PONG
    # Nothing to do here; no reply is needed for a pong message.
  else  # %xB-F are reserved for further control frames
    error("Unknown opcode $(wsf.opcode)")
  end
end

# Read a frame: turn bytes from the websocket into a WebSocketFragment.
function read_frame(io::IO)
  a = read(io,UInt8)
  fin    = a & 0b1000_0000 >>> 7  # If fin, then is final fragment
  rsv1   = a & 0b0100_0000  # If not 0, fail.
  rsv2   = a & 0b0010_0000  # If not 0, fail.
  rsv3   = a & 0b0001_0000  # If not 0, fail.
  opcode = a & 0b0000_1111  # If not known code, fail.
  # TODO: add validation somewhere to ensure rsv, opcode, mask, etc are valid.

  b = read(io,UInt8)
  mask = b & 0b1000_0000 >>> 7  # If not 1, fail.

  if mask != 1
      error("WebSocket reader cannot handle incoming messages without mask. " *
            "See http://tools.ietf.org/html/rfc6455#section-5.3")
  end

  payload_len::UInt64 = b & 0b0111_1111
  if payload_len == 126
    payload_len = ntoh(read(io,UInt16))  # 2 bytes
  elseif payload_len == 127
    payload_len = ntoh(read(io,UInt64))  # 8 bytes
  end

  maskkey = Array(UInt8,4)
  for i in 1:4
    maskkey[i] = read(io,UInt8)
  end

  data = Array(UInt8, payload_len)
  for i in 1:payload_len
    d = read(io, UInt8)
    d = d $ maskkey[mod(i - 1, 4) + 1]
    data[i] = d
  end

  return WebSocketFragment(fin,rsv1,rsv2,rsv3,opcode,mask,payload_len,maskkey,data)
end

# Read one non-control message from a WebSocket. Any control messages that are
# read will be handled by the handle_control_frame function. This function will
# not return until a full non-control message has been read. If the other side
# doesn't ever complete it's message, this function will never return. Only the
# data (contents/body/payload) of the message will be returned from this
# function.
function Base.read(ws::WebSocket)
  if ws.is_closed
    error("Attempt to read from closed WebSocket")
  end
  frame = read_frame(ws.socket)

  # Handle control (non-data) messages.
  if is_control_frame(frame)
    # Don't return control frames; they're not interesting to users.
    handle_control_frame(ws,frame)

    # Recurse to return the next data frame.
    return read(ws)
  end

  # Handle data message that uses multiple fragments.
  if !frame.is_last
    return vcat(frame.data, read(ws))
  end

  return frame.data
end

#
# WebSocket Handshake
#

# This function transforms a websocket client key into the server's accept
# value. This is done in three steps:
#   1. Concatenate key with magic string from RFC.
#   2. SHA1 hash the resulting base64 string.
#   3. Encode the resulting number in base64.
# This function then returns the string of the base64-encoded value.
function generate_websocket_key(key)
    hashed_key = digest(MD_SHA1, key*"258EAFA5-E914-47DA-95CA-C5AB0DC85B11")
    String(encode(Base64, hashed_key))
end

# Responds to a WebSocket handshake request.
# Checks for required headers; sends Response(400) if they're missing or bad.
# Otherwise, transforms client key into accept value, and sends Reponse(101).
function websocket_handshake(request,client)
  if !haskey(request.headers, "Sec-WebSocket-Key")
    Base.write(client.sock, Response(400))
    return
  end
  if get(request.headers, "Sec-WebSocket-Version", "13") != "13"
    response = Response(400)
    response.headers["Sec-WebSocket-Version"] = "13"
    Base.write(client.sock, response)
    return
  end

  key = request.headers["Sec-WebSocket-Key"]
  if length(decode(Base64,key)) != 16 # Key must be 16 bytes
    Base.write(client.sock, Response(400))
    return
  end
  resp_key = generate_websocket_key(key)

  response = Response(101)
  response.headers["Upgrade"] = "websocket"
  response.headers["Connection"] = "Upgrade"
  response.headers["Sec-WebSocket-Accept"] = resp_key
  Base.write(client.sock, response)
end

# Implement the WebSocketInterface, for compatilibility with HttpServer.
immutable WebSocketHandler <: HttpServer.WebSocketInterface
    handle::Function
end

import HttpServer: handle, is_websocket_handshake
function handle(handler::WebSocketHandler, req::Request, client::HttpServer.Client)
    websocket_handshake(req, client)
    sock = WebSocket(client.id, client.sock)
    handler.handle(req, sock)
    isopen(sock) && close(sock)
end
function is_websocket_handshake(handler::WebSocketHandler, req::Request)
    is_get = req.method == "GET"
    # "upgrade" for Chrome and "keep-alive, upgrade" for Firefox.
    is_upgrade = contains(lowercase(get(req.headers, "Connection", "")),"upgrade")
    is_websockets = lowercase(get(req.headers, "Upgrade", "")) == "websocket"
    return is_get && is_upgrade && is_websockets
end

end # module WebSockets


text:
using Revise
using AdFem
using SparseArrays
using LinearAlgebra
using PyPlot
αm = 0.0
αf = 0.0
β2 = 0.5*(1 - αm + αf)^2
γ = 0.5 - αm + αf

m = 40
n = 20
h = 0.01
NT = 100
Δt = 1/NT
bdedge = []
for j = 1:n
    push!(bdedge, [(j-1)*(m+1)+m+1 j*(m+1)+m+1])
end
bdedge = vcat(bdedge...)

bdnode = Int64[]
for j = 1:n+1
    push!(bdnode, (j-1)*(m+1)+1)
end

M = compute_fem_mass_matrix1(m, n, h)
S = spzeros((m+1)*(n+1), (m+1)*(n+1))
M = [M S;S M]
H = diagm(0=>[1,1,0.5])
K = 0.1
σY = 0.03
# σY = 1000.

state = zeros(2(m+1)*(n+1),NT+1)
velo = zeros(2(m+1)*(n+1),NT+1)
acce = zeros(2(m+1)*(n+1),NT+1)
stress = zeros(NT+1, 4*m*n, 3)
internal_variable = zeros(NT+1, 4*m*n)

t = 0.0
for i = 1:NT
    @info i
        ##### STEP 1: Computes the external force ####

    T = eval_f_on_boundary_edge((x,y)->0.02*sin(2π*i*Δt), bdedge, m, n, h)
    # T = eval_f_on_boundary_edge((x,y)->0.0, bdedge, m, n, h)
    T = [zeros(length(T)) -T]
    T = compute_fem_traction_term(T, bdedge, m, n, h)
    f1 = eval_f_on_gauss_pts((x,y)->0., m, n, h)
    f2 = eval_f_on_gauss_pts((x,y)->0., m, n, h)
    # f2 = eval_f_on_gauss_pts((x,y)->0.1, m, n, h)
    F = compute_fem_source_term(f1, f2, m, n, h)
    fext = F+T

        ##### STEP 2: Extract Variables ####

    u = state[:,i]
    ∂∂u = acce[:,i]
    ∂u = velo[:,i]

    ε0 = eval_strain_on_gauss_pts(u, m, n, h)
    σ0 = stress[i,:,:]
    α0 = internal_variable[i,:]

        ##### STEP 3: Newton Iteration ####

    global t += (1 - αf)*Δt
    ∂∂up = ∂∂u[:]
    iter = 0
    while true
        iter += 1
        # @info iter
        up = (1 - αf)*(u + Δt*∂u + 0.5 * Δt^2 * ((1 - β2)*∂∂u + β2*∂∂up)) + αf*u
        global fint, stiff, α, σ = compute_planestressplasticity_stress_and_stiffness_matrix(
            up, ε0, σ0, α0, K, σY, H, m, n, h
        )
        res = M * (∂∂up *(1 - αm) + αm*∂∂u)  + fint - fext
        A = M*(1 - αm) + (1 - αf) * 0.5 * β2 * Δt^2 * stiff
        A, _ = fem_impose_Dirichlet_boundary_condition(A, bdnode, m, n, h)
        res[[bdnode; bdnode .+ (m+1)*(n+1)]] .= 0.0

        Δ∂∂u = A\res
        ∂∂up -= Δ∂∂u
        err = norm(res)
        # @info err
        if err<1e-8
            break
        end
    end
    global t += αf*Δt

        ##### STEP 3: Update State Variables ####

    u += Δt * ∂u + Δt^2/2 * ((1 - β2) * ∂∂u + β2 * ∂∂up)
    ∂u += Δt * ((1 - γ) * ∂∂u + γ * ∂∂up)
    stress[i+1,:,:] = σ
    internal_variable[i+1,:] = α

    state[:,i+1] = u
    acce[:,i+1] = ∂∂up
    velo[:,i+1] = ∂u
end

x = []
y = []
for j= 1:n+1
    for i = 1:m+1
        push!(x, (i-1)*h)
        push!(y, (j-1)*h)
    end
end

for i = 1:5:NT+1
    close("all")
    scatter(x+state[1:(m+1)*(n+1), i], y+state[(m+1)*(n+1)+1:end, i])
    scatter(x[m+1]+state[m+1, i],
            y[m+1]+state[(m+1)*(n+1)+m+1, i], color="red")
    xlabel("x")
    ylabel("y")
    k = string(i)
    k = repeat("0", 3-length(k))*k
    title("t = $k")
    ylim(-0.05,0.25)
    xlim(-0.01, 0.45)
    gca().invert_yaxis()
    savefig("u$k.png")

    close("all");
    plot(1:i, -state[(m+1)*(n+1)+m+1, 1:i])
    xlim(0, NT+2)
    ylim(0, 0.03)
    scatter(i, -state[(m+1)*(n+1)+m+1, i], color="red")
    savefig("du$k.png")
end

run(`convert -delay 10 -loop 0 u*.png plasticity_u.gif`)
run(`convert -delay 10 -loop 0 du*.png plasticity_du.gif`)


text:
function contact_normal_τ_direct!(τ,sim_data::SimData,h,Hi,envj::EnvironmentJacobian,dyn_bias,u0,v0,x_upper::AbstractArray{U},n::Int) where U
    num_contacts = length(sim_data.env.contacts)
    env = sim_data.env
    lower_vs = sim_data.normal_vs[n]
    lower_cs = sim_data.normal_cs[n]
    lower_options = sim_data.normal_options[n]

    ϕAs = []
    ϕbs = []
    for i = 1:num_contacts
        J = envj.contact_jacobians[i].J
        ϕ = envj.contact_jacobians[i].ϕ
        N = envj.contact_jacobians[i].N

        ϕA = h^2*N*Hi*J
        ϕb = N*(h^2*Hi*(dyn_bias - u0) .- h*v0) .- ϕ

        push!(ϕAs,ϕA)
        push!(ϕbs,ϕb)
    end

    function eval_obj_(x::AbstractArray{L}) where L
        obj = 0.

        for i = 1:num_contacts
            c_n = lower_vs(x, Symbol("c_n", i))
            obj += c_n'*c_n
        end

        obj
    end

    function eval_cons_(x::AbstractArray{L}) where L
        # TODO use the cs
        g = []

        for i = 1:num_contacts
            c_n = lower_vs(x, Symbol("c_n", i))
            g = vcat(g, -c_n)
            g = vcat(g, ϕAs[i]*vcat(c_n, zeros(4)) + ϕbs[i])
        end

        g
    end

    fres = DiffResults.HessianResult(zeros(U, lower_vs.num_vars))
    gres = DiffResults.JacobianResult(zeros(U, lower_cs.num_cons), zeros(U, lower_vs.num_vars))
    solver_fn_ = generate_autodiff_solver_fn(eval_obj_,fres,eval_cons_,gres,lower_cs.eqs,lower_cs.ineqs)

    x0 = zeros(lower_vs.num_vars)

    xopt, info = auglag(solver_fn_, lower_cs.num_eqs, lower_cs.num_ineqs, x0, lower_options)

    # TODO include the total weight here, not in J
    τ .= mapreduce(+, enumerate(envj.contact_jacobians)) do (i,cj)
        contact_τ(cj, lower_vs(xopt, Symbol("c_n", i)), zeros(4))
    end

    xopt
end

function contact_normal_τ_direct_osqp!(τ,sim_data::SimData,h,Hi,envj::EnvironmentJacobian,dyn_bias,u0,v0,x_upper::AbstractArray{U},n::Int) where U
    num_contacts = length(sim_data.env.contacts)
    env = sim_data.env
    lower_vs = sim_data.normal_vs[n]
    lower_cs = sim_data.normal_cs[n]
    lower_options = sim_data.normal_options[n]

    ϕAs = []
    ϕbs = []
    for i = 1:num_contacts
        J = envj.contact_jacobians[i].J
        ϕ = envj.contact_jacobians[i].ϕ
        N = envj.contact_jacobians[i].N

        ϕA = h^2*N*Hi*J
        ϕb = N*(h^2*Hi*(dyn_bias - u0) .- h*v0) .- ϕ

        push!(ϕAs,ϕA)
        push!(ϕbs,ϕb)
    end

    xopt = Array{U,1}(undef, lower_vs.num_vars)
    P = ones(1,1)
    q = [0.]
    l = [0., -1e19]
    for i = 1:num_contacts
        A = Array{eltype(ϕAs[i]),2}(undef, 2, 1)
        A[1,1] = 1.
        A[2,1] = ϕAs[i][1,1]
        u = Array{eltype(ϕbs[i]),1}(undef, 2)
        u[1] = 1e19
        u[2] = -ϕbs[i][1]
        xopt_i = osqp(P,q,A,l,u,sim_data.normal_osqp_models[n])
        xopt[lower_vs(Symbol("c_n", i))] = xopt_i
    end

    # TODO include the total weight here, not in J
    τ .= mapreduce(+, enumerate(envj.contact_jacobians)) do (i,cj)
        contact_τ(cj, lower_vs(xopt, Symbol("c_n", i)), zeros(4))
    end

    xopt
end

function contact_friction_τ_direct!(τ,sim_data::SimData,h,Hi,envj::EnvironmentJacobian,dyn_bias,u0,v0,x_upper::AbstractArray{U},x_normal::AbstractArray{K},n::Int) where {U,K}
    num_contacts = length(sim_data.env.contacts)
    env = sim_data.env
    normal_vs = sim_data.normal_vs[n]
    lower_vs = sim_data.fric_vs[n]
    lower_cs = sim_data.fric_cs[n]
    lower_options = sim_data.fric_options[n]

    Qds = []
    rds = []
    for i = 1:num_contacts
        J = envj.contact_jacobians[i].J
        ϕ = envj.contact_jacobians[i].ϕ
        N = envj.contact_jacobians[i].N

        Qd = h^2*J'*Hi*J
        rd = J'*(h^2*Hi*(dyn_bias - u0) .- h*v0)

        push!(Qds,Qd)
        push!(rds,rd)
    end

    function eval_obj_(x::AbstractArray{L}) where L
        obj = 0.

        for i = 1:num_contacts
            c_n = normal_vs(x_normal, Symbol("c_n", i))
            β = lower_vs(x, Symbol("β", i))
            z = vcat(c_n, β)
            obj += .5*z'*Qds[i]*z + z'*rds[i]
        end

        obj
    end

    function eval_cons_(x::AbstractArray{L}) where L
        # TODO use the cs
        g = []

        for i = 1:num_contacts
            c_n = normal_vs(x_normal, Symbol("c_n", i))
            β = lower_vs(x, Symbol("β", i))
            g = vcat(g, -β)
            g = vcat(g, sum(β) .- env.contacts[i].obstacle.μ * c_n)
        end

        g
    end

    fres = DiffResults.HessianResult(zeros(U, lower_vs.num_vars))
    gres = DiffResults.JacobianResult(zeros(U, lower_cs.num_cons), zeros(U, lower_vs.num_vars))
    solver_fn_ = generate_autodiff_solver_fn(eval_obj_,fres,eval_cons_,gres,lower_cs.eqs,lower_cs.ineqs)

    x0 = zeros(lower_vs.num_vars)

    xopt, info = auglag(solver_fn_, lower_cs.num_eqs, lower_cs.num_ineqs, x0, lower_options)

    # TODO include the total weight here, not in J
    τ .= mapreduce(+, enumerate(envj.contact_jacobians)) do (i,cj)
        contact_τ(cj, normal_vs(x_normal, Symbol("c_n", i)), lower_vs(xopt, Symbol("β", i)))
    end

    xopt
end

function contact_friction_τ_direct!(τ,sim_data::SimData,h,Hi,envj::EnvironmentJacobian,dyn_bias,u0,v0,x_upper::AbstractArray{U},n::Int) where U
    num_contacts = length(sim_data.env.contacts)
    env = sim_data.env
    upper_vs = sim_data.vs
    lower_vs = sim_data.fric_vs[n]
    lower_cs = sim_data.fric_cs[n]
    lower_options = sim_data.fric_options[n]

    Qds = []
    rds = []
    for i = 1:num_contacts
        J = envj.contact_jacobians[i].J
        ϕ = envj.contact_jacobians[i].ϕ
        N = envj.contact_jacobians[i].N

        Qd = h^2*J'*Hi*J
        rd = J'*(h^2*Hi*(dyn_bias - u0) .- h*v0)

        push!(Qds,Qd)
        push!(rds,rd)
    end

    function eval_obj_(x::AbstractArray{L}) where L
        obj = 0.

        for i = 1:num_contacts
            c_n = upper_vs(x_upper, Symbol("c_n", i, "_", n))
            β = lower_vs(x, Symbol("β", i))
            z = vcat(c_n, β)
            obj += .5*z'*Qds[i]*z + z'*rds[i]
        end

        obj
    end

    function eval_cons_(x::AbstractArray{L}) where L
        # TODO use the cs
        g = []

        for i = 1:num_contacts
            c_n = upper_vs(x_upper, Symbol("c_n", i, "_", n))
            β = lower_vs(x, Symbol("β", i))
            g = vcat(g, -β)
            g = vcat(g, sum(β) .- env.contacts[i].obstacle.μ * c_n)
        end

        g
    end

    fres = DiffResults.HessianResult(zeros(U, lower_vs.num_vars))
    gres = DiffResults.JacobianResult(zeros(U, lower_cs.num_cons), zeros(U, lower_vs.num_vars))
    solver_fn_ = generate_autodiff_solver_fn(eval_obj_,fres,eval_cons_,gres,lower_cs.eqs,lower_cs.ineqs)

    x0 = zeros(lower_vs.num_vars)

    xopt, info = auglag(solver_fn_, lower_cs.num_eqs, lower_cs.num_ineqs, x0, lower_options)

    # TODO include the total weight here, not in J
    τ .= mapreduce(+, enumerate(envj.contact_jacobians)) do (i,cj)
        contact_τ(cj, upper_vs(x_upper, Symbol("c_n", i, "_", n)), lower_vs(xopt, Symbol("β", i)))
    end

    xopt
end

function contact_friction_τ_direct_osqp!(τ,sim_data::SimData,h,Hi,envj::EnvironmentJacobian,dyn_bias,u0,v0,x_upper::AbstractArray{U},x_normal::AbstractArray{K},n::Int) where {U,K}
    num_contacts = length(sim_data.env.contacts)
    env = sim_data.env
    normal_vs = sim_data.normal_vs[n]
    lower_vs = sim_data.fric_vs[n]
    lower_cs = sim_data.fric_cs[n]
    lower_options = sim_data.fric_options[n]

    Qds = []
    rds = []
    for i = 1:num_contacts
        J = envj.contact_jacobians[i].J
        ϕ = envj.contact_jacobians[i].ϕ
        N = envj.contact_jacobians[i].N

        Qd = h^2*J'*Hi*J
        rd = J'*(h^2*Hi*(dyn_bias - u0) .- h*v0)

        push!(Qds,Qd)
        push!(rds,rd)
    end

    xopt = Array{U,1}(undef, lower_vs.num_vars)
    # TODO use the β selector
    β_max = 1e8
    A = [1. 0. 0. 0.; 0. 1. 0. 0.; 0. 0. 1. 0.; 0. 0. 0. 1.; 1. 1. 1. 1.]
    l = [0., 0., 0., 0., -4. * β_max]
    for i = 1:num_contacts
        c_n = normal_vs(x_normal, Symbol("c_n", i))[1]
        P = Qds[i][2:end,2:end]
        P = .5 * P' * P
        P[abs.(P) .< 1e-6] .= 0.
        P += I*1e-8
        q = (rds[i][2:end] + .5*Qds[i][1,2:end]*c_n + .5*Qds[i][2:end,1]*c_n)
        u = Array{eltype(x_normal), 1}([β_max, β_max, β_max, β_max, env.contacts[i].obstacle.μ*c_n])
        xopt_i = osqp(P,q,A,l,u,sim_data.fric_osqp_models[n])
        xopt[lower_vs(Symbol("β", i))] = xopt_i
    end

    # TODO include the total weight here, not in J
    τ .= mapreduce(+, enumerate(envj.contact_jacobians)) do (i,cj)
        contact_τ(cj, normal_vs(x_normal, Symbol("c_n", i)), lower_vs(xopt, Symbol("β", i)))
    end

    xopt
end

function contact_friction_τ_direct_osqp!(τ,sim_data::SimData,h,Hi,envj::EnvironmentJacobian,dyn_bias,u0,v0,x_upper::AbstractArray{U},n::Int) where U
    num_contacts = length(sim_data.env.contacts)
    env = sim_data.env
    upper_vs = sim_data.vs
    lower_vs = sim_data.fric_vs[n]
    lower_cs = sim_data.fric_cs[n]
    lower_options = sim_data.fric_options[n]

    Qds = []
    rds = []
    for i = 1:num_contacts
        J = envj.contact_jacobians[i].J
        ϕ = envj.contact_jacobians[i].ϕ
        N = envj.contact_jacobians[i].N

        Qd = h^2*J'*Hi*J
        rd = J'*(h^2*Hi*(dyn_bias - u0) .- h*v0)

        push!(Qds,Qd)
        push!(rds,rd)
    end

    xopt = Array{U,1}(undef, lower_vs.num_vars)
    # TODO use the β selector
    β_max = 1e8
    A = [1. 0. 0. 0.; 0. 1. 0. 0.; 0. 0. 1. 0.; 0. 0. 0. 1.; 1. 1. 1. 1.]
    l = [0., 0., 0., 0., -4. * β_max]
    for i = 1:num_contacts
        c_n = upper_vs(x_upper, Symbol("c_n", i, "_", n))[1]
        P = Qds[i][2:end,2:end]
        P = .5 * P' * P
        P[abs.(P) .< 1e-6] .= 0.
        P += I*1e-8
        q = (rds[i][2:end] + .5*Qds[i][1,2:end]*c_n + .5*Qds[i][2:end,1]*c_n)
        u = Array{eltype(x_upper), 1}([β_max, β_max, β_max, β_max, env.contacts[i].obstacle.μ*c_n])
        xopt_i = osqp(P,q,A,l,u,sim_data.fric_osqp_models[n])
        xopt[lower_vs(Symbol("β", i))] = xopt_i
    end

    # TODO include the total weight here, not in J
    τ .= mapreduce(+, enumerate(envj.contact_jacobians)) do (i,cj)
        contact_τ(cj, upper_vs(x_upper, Symbol("c_n", i, "_", n)), lower_vs(xopt, Symbol("β", i)))
    end

    xopt
end

function contact_τ_indirect!(τ::AbstractArray{T},sim_data::SimData,envj::EnvironmentJacobian{T},x::AbstractArray{T}) where T
    # TODO: parallel
    τ .= mapreduce(+, enumerate(envj.contact_jacobians)) do (i,cj)
        contact_τ(cj, sim_data.vs(x, Symbol("c_n", i)), sim_data.vs(x, Symbol("β", i)))
    end
end

function contact_τ_indirect!(τ::AbstractArray{T},sim_data::SimData,envj::EnvironmentJacobian{T},x::AbstractArray{T},n::Int) where T
    # TODO: parallel
    τ .= mapreduce(+, enumerate(envj.contact_jacobians)) do (i,cj)
        contact_τ(cj, sim_data.vs(x, Symbol("c_n", i, "_", n)), sim_data.vs(x, Symbol("β", i, "_", n)))
    end
end