text:
module DomTerm

include("display.jl")
include("utilities.jl")

isDomterm = false

export
    # status info
    dt_status, dt_list,
    # stylesheets
    dt_listst, dt_printst, dt_addst,
    # screen setting
    dt_reversevideo, dt_clear, dt_winsizex,
    # general setting
    dt_settings,
    # file display
    hput, htmlcat, svgcat, jpegcat, pngcat, gifcat,
    # helper functions
    displayinline,
    # macros
    @dt

function __init__()
    check_domterm()
    pushdisplay(DomTermInline())
end # __init__

end # DomTerm


text:
import Base.clamp
export clamp
function clamp(x::AFArray, lo, hi)
  max(min(x,hi),lo)
end

ext: jl


text:
function erarawregion(
    emod::Dict, epar::Dict, ereg::Dict,
    start::TimeType, finish::TimeType;
    path::AbstractString="",
    region::AbstractString="GLB"
)

    if path == ""; eroot = eraroot(actionID); else; eroot = eraroot(path,actionID); end

    @info "$(Dates.now()) - Extracting $(uppercase(emod["dataset"])) $(epar["name"]) data for the entire $(gregionfullname(region)) region ..."

    if !isdir(eraregfolder(ereg,eroot))
        error("$(Dates.now()) - No data has been downloaded from $(uppercase(emod["dataset"])) $(epar["name"]) in the $(gregionfullname(region))")
    end

    dvec,dys,dyf,ndy = extractdate(start,finish); ndates = length(dvec);
    if emod["datasetID"] == 1; nt = 24; else; nt = 4; end

    lon = ereg["lon"]; lat = ereg["lat"]; rlon,rlat,rinfo = gregiongridvec(region,lon,lat);
    nlon = length(rlon); nlat = length(rlat); data = zeros(Float32,nlon,nlat,ndy*nt);

    for ii = 1 : ndates; dateii = dvec[ii];

        fol = erarawfolder(epar,ereg,eroot);
        if !isdir(fol)
            error("$(Dates.now()) - There is no data for $(uppercase(emod["dataset"])) $(epar["name"]) for $(yrmo2dir(dateii)).")
        end

        fnc = joinpath(fol,erarawname(emod,epar,ereg,dateii));
        ds  = Dataset(fnc,"r"); vds = ds[varname];

        if     ii == 1 && ii != ndates
            moday = daysinmonth(dateii);
            ibeg = (dys-1)*nt + 1; iend = (moday+1-dys)*nt;
            data[:,:,1:iend] = vds[:,:,ibeg:end];
        elseif ii != 1 && ii == ndates
            ibeg = iend+1; iend = dyf*nt;
            data[:,:,ibeg:end] = vds[:,:,1:iend];
        elseif ii == 1 && ii == ndates
            ibeg = (dys-1)*nt + 1; iend = dyf*nt;
            data = vds[:,:,ibeg:iend];
        else
            moday = daysinmonth(dateii);
            ibeg = iend+1; iend = ibeg-1 + moday*nt;
            data[:,:,ibeg:iend] = vds[:];
        end

        close(ds);

    end

    @info "$(Dates.now()) - $(uppercase(emod["dataset"])) $(epar["name"]) data for the entire $(gregionfullname(region)) region has been extracted."

    return datavec,info,[rlon,rlat]

end

function erarawpoint(
    emod::Dict, epar::Dict, ereg::Dict,
    start::TimeType, finish::TimeType;
    coord::Array{<:Real,1},
    path::AbstractString="",
    region::AbstractString="GLB"
)

    if path == ""; eroot = eraroot(actionID); else; eroot = eraroot(path,actionID); end
    if length(coord) != 2
        error("$(Dates.now()) - Coordinate vector must be in the form [lon,lat]")
    end

    @info "$(Dates.now()) - Extracting $(uppercase(emod["dataset"])) $(epar["name"]) data at coordinates $(coord) ..."

    if !isdir(eraregfolder(ereg,eroot))
        error("$(Dates.now()) - No data has been downloaded from $(uppercase(emod["dataset"])) $(epar["name"]) in the $(gregionfullname(region))")
    end

    dvec,dys,dyf,ndy = extractdate(start,finish); ndates = length(dvec);
    if emod["datasetID"] == 1; nt = 24; else; nt = 4; end

    lon = ereg["lon"]; lat = ereg["lat"]; rlon,rlat,rinfo = gregiongridvec(region,lon,lat);
    plon,plat = coord; ilon,ilat = regionpoint(plon,plat,rlon,rlat);
    data = zeros(Float32,ndy*nt);

    for ii = 1 : ndates; dateii = dvec[ii];

        fol = erarawfolder(epar,ereg,eroot);
        if !isdir(fol)
            error("$(Dates.now()) - There is no data for $(uppercase(emod["dataset"])) $(epar["name"]) for $(yrmo2dir(dateii)).")
        end

        fnc = joinpath(fol,erarawname(emod,epar,ereg,dateii));
        ds  = Dataset(fnc,"r"); vds = ds[varname];

        if     ii == 1 && ii != ndates
            moday = daysinmonth(dateii);
            ibeg = (dys-1)*nt + 1; iend = (moday+1-dys)*nt;
            data[1:iend] = vds[ilon,ilat,ibeg:end];
        elseif ii != 1 && ii == ndates
            ibeg = iend+1; iend = dyf*nt;
            data[ibeg:end] = vds[ilon,ilat,1:iend];
        elseif ii == 1 && ii == ndates
            ibeg = (dys-1)*nt + 1; iend = dyf*nt;
            data = vds[ilon,ilat,ibeg:iend];
        else
            moday = daysinmonth(dateii);
            ibeg = iend+1; iend = ibeg-1 + moday*nt;
            data[ibeg:iend] = vds[ilon,ilat,:];
        end

        close(ds);

    end

    @info "$(Dates.now()) - $(uppercase(emod["dataset"])) $(epar["name"]) data for the coordinates $(coord) has been extracted."

    return datavec,info

end

function erarawgrid(
    emod::Dict, epar::Dict, ereg::Dict,
    start::TimeType, finish::TimeType;
    grid::Array{<:Real,1},
    path::AbstractString="",
    region::AbstractString="GLB"
)

    if path == ""; eroot = eraroot(actionID); else; eroot = eraroot(path,actionID); end
    if length(grid) != 4
        error("$(Dates.now()) - Grid vector must be in the form [N,S,E,W]")
    end

    @info "$(Dates.now()) - Extracting $(uppercase(emod["dataset"])) $(epar["name"]) data for the [N,S,E,W] bounds $(grid) ..."

    if !isdir(eraregfolder(ereg,eroot))
        error("$(Dates.now()) - No data has been downloaded from $(uppercase(emod["dataset"])) $(epar["name"]) in the $(gregionfullname(region))")
    end

    dvec,dys,dyf,ndy = extractdate(start,finish); ndates = length(dvec);
    if emod["datasetID"] == 1; nt = 24; else; nt = 4; end

    lon = ereg["lon"]; lat = ereg["lat"]; isgridinregion(grid,region);
    rlon,rlat,rinfo = gregiongridvec(region,lon,lat);
    glon,glat,ginfo = regiongridvec(grid,rlon,rlat); iWE,iNS = ginfo["IDvec"];
    nlon = length(glon); nlat = length(glat);
    data = zeros(Float32,nlon,nlat,ndy*nt);

    for ii = 1 : ndates; dateii = dvec[ii];

        fol = erarawfolder(epar,ereg,eroot);
        if !isdir(fol)
            error("$(Dates.now()) - There is no data for $(uppercase(emod["dataset"])) $(epar["name"]) for $(yrmo2dir(dateii)).")
        end

        fnc = joinpath(fol,erarawname(emod,epar,ereg,dateii));
        ds  = Dataset(fnc,"r"); vds = ds[varname];

        if     ii == 1 && ii != ndates
            moday = daysinmonth(dateii);
            ibeg = (dys-1)*nt + 1; iend = (moday+1-dys)*nt;
            data[:,:,1:iend] = vds[iWE,iNS,ibeg:end];
        elseif ii != 1 && ii == ndates
            ibeg = iend+1; iend = dyf*nt;
            data[:,:,ibeg:end] = vds[iWE,iNS,1:iend];
        elseif ii == 1 && ii == ndates
            ibeg = (dys-1)*nt + 1; iend = dyf*nt;
            data = vds[iWE,iNS,ibeg:iend];
        else
            moday = daysinmonth(dateii);
            ibeg = iend+1; iend = ibeg-1 + moday*nt;
            data[:,:,ibeg:iend] = vds[iWE,iNS,:];
        end

        close(ds);

    end

    @info "$(Dates.now()) - $(uppercase(emod["dataset"])) $(epar["name"]) data within the [N,S,E,W] bounds $(grid) has been extracted."

    return datavec,info,[glon,glat]

end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct LinkedModelFile <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function LinkedModelFile(; kwargs...)
        obj = new(meta(LinkedModelFile), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct LinkedModelFile
const __meta_LinkedModelFile = Ref{ProtoMeta}()
function meta(::Type{LinkedModelFile})
    ProtoBuf.metalock() do
        if !isassigned(__meta_LinkedModelFile)
            __meta_LinkedModelFile[] = target = ProtoMeta(LinkedModelFile)
            allflds = Pair{Symbol,Union{Type,String}}[:linkedModelFileName => StringParameter, :linkedModelSearchPath => StringParameter]
            meta(target, LinkedModelFile, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_LinkedModelFile[]
    end
end
function Base.getproperty(obj::LinkedModelFile, name::Symbol)
    if name === :linkedModelFileName
        return (obj.__protobuf_jl_internal_values[name])::StringParameter
    elseif name === :linkedModelSearchPath
        return (obj.__protobuf_jl_internal_values[name])::StringParameter
    else
        getfield(obj, name)
    end
end

mutable struct LinkedModel <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function LinkedModel(; kwargs...)
        obj = new(meta(LinkedModel), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct LinkedModel
const __meta_LinkedModel = Ref{ProtoMeta}()
function meta(::Type{LinkedModel})
    ProtoBuf.metalock() do
        if !isassigned(__meta_LinkedModel)
            __meta_LinkedModel[] = target = ProtoMeta(LinkedModel)
            allflds = Pair{Symbol,Union{Type,String}}[:linkedModelFile => LinkedModelFile]
            oneofs = Int[1]
            oneof_names = Symbol[Symbol("LinkType")]
            meta(target, LinkedModel, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_LinkedModel[]
    end
end
function Base.getproperty(obj::LinkedModel, name::Symbol)
    if name === :linkedModelFile
        return (obj.__protobuf_jl_internal_values[name])::LinkedModelFile
    else
        getfield(obj, name)
    end
end

export LinkedModel, LinkedModelFile


text:
export EthHdr, IpFlags, IpHdr,
       UdpHdr, TcpFlags, TcpHdr,
       IcmpHdr, DecPkt, decode_pkt

type EthHdr
    dest_mac::AbstractString
    src_mac::AbstractString
    ptype::UInt16
    EthHdr() = new("","",0)
end # type EthHdr

type IpFlags
    reserved::Bool
    dont_frag::Bool
    more_frags::Bool
    IpFlags() = new(false,false,false)
end # type IpFlags

type IpHdr
    version::UInt8
    length::UInt8
    services::UInt8
    totlen::UInt16
    id::UInt16
    flags::IpFlags
    frag_offset::UInt16
    ttl::UInt8
    protocol::UInt8
    checksum::Bool
    src_ip::AbstractString
    dest_ip::AbstractString
    IpHdr() = new(0,0,0,0,0,IpFlags(),0,0,0,false,"","")
end # type IpHdr

type TcpFlags
    reserved::Bool
    nonce::Bool
    cwr::Bool
    ecn::Bool
    urgent::Bool
    ack::Bool
    push::Bool
    reset::Bool
    syn::Bool
    fin::Bool
    TcpFlags() = new(false,false,false,false,false,
                     false,false,false,false,false)
end # type TcpFlags

type TcpHdr
    src_port::UInt16
    dest_port::UInt16
    seq::UInt32
    ack::UInt32
    offset::UInt8
    flags::TcpFlags
    window::UInt16
    checksum::UInt16
    uptr::UInt16
    data::Vector{UInt8}
    TcpHdr() = new(0,0,0,0,0,TcpFlags(),0,0,0, Vector{UInt8}(0))
end # type TcpHdr

type UdpHdr
    src_port::UInt16
    dest_port::UInt16
    length::UInt16
    checksum::UInt16
    data::Vector{UInt8}
    UdpHdr() = new(0,0,0,0,Vector{UInt8}(0))
end # type UdpHdr

type IcmpHdr
    ptype::UInt8
    code::UInt8
    checksum::UInt16
    identifier::UInt16
    seqno::UInt16
    IcmpHdr() = new(0,0,0,0,0)
end # type IcmpHdr

type DecPkt
    datalink::EthHdr
    network::IpHdr
    protocol::Any
    DecPkt() = new(EthHdr(), IpHdr(), nothing)
end # type DecPkt

@inline function getindex_he{T}(::Type{T}, b::Vector{UInt8}, i)
    # When 0.4 support is dropped, add @boundscheck
    checkbounds(b, i + sizeof(T) - 1)
    return unsafe_load(Ptr{T}(pointer(b, i)))
end

@inline getindex_be{T}(::Type{T}, b::Vector{UInt8}, i) = hton(getindex_he(T, b, i))

#----------
# decode ethernet header
#----------
function decode_eth_hdr(d::Array{UInt8})
    eh = EthHdr()
    eh.dest_mac = string(hex(d[1], 2), ":", hex(d[2], 2), ":", hex(d[3], 2), ":",
                         hex(d[4], 2), ":", hex(d[5], 2), ":", hex(d[6], 2))
    eh.src_mac  = string(hex(d[7], 2), ":", hex(d[8], 2), ":", hex(d[9], 2), ":",
                         hex(d[10], 2), ":", hex(d[11], 2), ":", hex(d[12], 2))
    eh.ptype    = getindex_be(UInt16, d, 13)
    eh
end # function decode_eth_hdr

#----------
# calculate IP checksum
#----------
function ip_checksum(buf::Array{UInt8})
    sum::UInt64 = 0
    for i in 1:2:length(buf)
        pair = getindex_he(UInt16, buf, i)
        sum += pair
        if (sum & 0x80000000) != 0
            sum = (sum & 0xFFFF) + (sum >> 16)
        end
    end

    while ((sum >> 16) != 0)
        sum = (sum & 0xFFFF) + (sum >> 16)
    end
    ~sum
end # function ip_checksum

#----------
# decode IP header
#----------
function decode_ip_hdr(d::Array{UInt8})
    iph = IpHdr()
    iph.version     = (d[1] & 0xf0) >> 4
    iph.length      = (d[1] & 0x0f) * 4
    if ip_checksum(d[1:iph.length]) == 0xFFFFFFFFFFFF0000
        iph.checksum = true
    end
    iph.services    = d[2]
    iph.totlen      = getindex_be(UInt16, d, 3)
    iph.id          = getindex_be(UInt16, d, 5)

    # set flags
    flags = IpFlags()
    flags.reserved   = (d[7] & (1 << 7)) > 0
    flags.dont_frag  = (d[7] & (1 << 6)) > 0
    flags.more_frags = (d[7] & (1 << 5)) > 0
    iph.flags        = flags

    iph.frag_offset = getindex_be(UInt16, d, 7) & 0x7ff
    iph.ttl         = d[9]
    iph.protocol    = d[10]
    iph.src_ip      = string(Int(d[13]), ".", Int(d[14]), ".", Int(d[15]), ".", Int(d[16]))
    iph.dest_ip     = string(Int(d[17]), ".", Int(d[18]), ".", Int(d[19]), ".", Int(d[20]))
    iph
end # function decode_ip_hdr

#----------
# decode TCP header
#----------
function decode_tcp_hdr(d::Array{UInt8})
    tcph = TcpHdr()
    tcph.src_port  = getindex_be(UInt16, d, 1)
    tcph.dest_port = getindex_be(UInt16, d, 3)
    tcph.seq       = getindex_be(UInt32, d, 5)
    tcph.ack       = getindex_be(UInt32, d, 9)
    tcph.offset    = (d[13] & 0xf0) >> 4

    # set flags
    flags = TcpFlags()
    flags.reserved = ((d[13] & 0x0e) >> 1) > 0
    flags.nonce    = (d[13] & 1) > 0
    flags.cwr      = (d[14] & (1 << 7)) > 0
    flags.ecn      = (d[14] & (1 << 6)) > 0
    flags.urgent   = (d[14] & (1 << 5)) > 0
    flags.ack      = (d[14] & (1 << 4)) > 0
    flags.push     = (d[14] & (1 << 3)) > 0
    flags.reset    = (d[14] & (1 << 2)) > 0
    flags.syn      = (d[14] & (1 << 1)) > 0
    flags.fin      = (d[14] & 1) > 0
    tcph.flags     = flags

    tcph.window    = getindex_be(UInt16, d, 15)
    tcph.checksum  = getindex_be(UInt16, d, 17)
    tcph.uptr      = getindex_be(UInt16, d, 19)
    tcph.data      = d[tcph.offset * 4 + 1:end]
    tcph
end # function decode_tcp_hdr

#----------
# decode UDP header
#----------
function decode_udp_hdr(d::Array{UInt8})
    udph = UdpHdr()
    udph.src_port  = getindex_be(UInt16, d, 1)
    udph.dest_port = getindex_be(UInt16, d, 3)
    udph.length    = getindex_be(UInt16, d, 5)
    udph.checksum  = getindex_be(UInt16, d, 7)
    udph.data      = d[9:end]
    udph
end # function decode_udp_hdr

#----------
# decode ICMP header
#----------
function decode_icmp_hdr(d::Array{UInt8})
    icmph = IcmpHdr()
    icmph.ptype      = d[1]
    icmph.code       = d[2]
    icmph.checksum   = getindex_be(UInt16, d, 3)
    icmph.identifier = getindex_be(UInt16, d, 5)
    icmph.seqno      = getindex_be(UInt16, d, 7)
    icmph
end # function decode_icmp_hdr

#----------
# decode ethernet packet
#----------
function decode_pkt(pkt::Array{UInt8})
    decoded           = DecPkt()
    decoded.datalink  = decode_eth_hdr(pkt)
    iphdr             = decode_ip_hdr(pkt[15:end])
    decoded.network   = iphdr

    proto = nothing
    if (iphdr.protocol == 1)
        proto = decode_icmp_hdr(pkt[15 + iphdr.length:end])
    elseif (iphdr.protocol == 6)
        proto = decode_tcp_hdr(pkt[15 + iphdr.length:end])
    elseif (iphdr.protocol == 17)
        proto = decode_udp_hdr(pkt[15 + iphdr.length:end])
    end

    decoded.protocol = proto
    decoded
end # function decode_pkt


text:
## Methods related to the SEXP (pointer to SEXPREC type) in R

## length methods for the vector types of SEXP's.  To handle long vectors the length method
## should be extended to check for the sentinal value (typemin(Cint)) followed by extraction
## of the Int64 length.  The size methods check for a dim attribute and return a tuple.
for N in [LGLSXP,INTSXP,REALSXP,CPLXSXP,STRSXP,VECSXP,EXPRSXP]
    @eval begin
        Base.length(s::SEXP{$N}) = unsafe_load(convert(Ptr{Cint},s.p+loffset),1)
        function Base.size(s::SEXP{$N})
            dd = reval(lang2(dimSymbol,s))
            isa(dd,SEXP{INTSXP}) || return (int64(length(s)),)
            vv = vec(dd)
            ntuple(length(vv),i->int64(vv[i]))
        end
    end
end

## Element extraction and iterators
for (N,elt) in ((STRSXP,:STRING_ELT),(VECSXP,:VECTOR_ELT),(EXPRSXP,:VECTOR_ELT))
    @eval begin
        function Base.getindex(s::SEXP{$N},I::Number)  # extract a single element
            0 < I ≤ length(s) || throw(BoundsError())
            sexp(ccall(($(string(elt)),libR),Ptr{Void},(Ptr{Void},Cint),s,I-1))
        end
        Base.start(s::SEXP{$N}) = 0  # start,next,done and eltype provide an iterator
        Base.next(s::SEXP{$N},state) = (state += 1;(s[state],state))
        Base.done(s::SEXP{$N},state) = state ≥ length(s)
        Base.eltype(s::SEXP{$N}) = SEXP
    end
end
Base.eltype(s::SEXP{STRSXP}) = SEXP{CHARSXP} # be more specific for STRSXP

## ToDo: write an isNA method for Complex128 - not sure how it is defined though.
isNA(x::Cdouble) = x == R_NaReal
isNA(x::Cint) = x == R_NaInt
isNA(x::Union(ASCIIString,UTF8String)) = x == bytestring(R_NaString)
isNA(a::Array) = reshape(bitpack([isNA(aa) for aa in a]),size(a))

## bytestring copies the contents of the 0-terminated string at the Ptr{Uint8} address
Base.bytestring(s::SEXP{CHARSXP}) = bytestring(ccall((:R_CHAR,libR),Ptr{Uint8},(Ptr{Void},),s))

Base.vec(s::SEXP{STRSXP}) = [bytestring(ss) for ss in s]

for (N,typ) in ((LGLSXP,:Int32),(INTSXP,:Int32),(REALSXP,:Float64),(CPLXSXP,:Complex128))
    @eval begin
        function Base.vec(s::SEXP{$N})
            ccall((:R_PreserveObject,libR),Void,(Ptr{Void},),s)
            rv = pointer_to_array(convert(Ptr{$typ},s.p+voffset),length(s))
            finalizer(rv,x->ccall((:R_ReleaseObject,libR),Void,(Ptr{Void},),s))
            rv
        end
    end
end

## Not sure what to do about R's Logical vectors (SEXP{10}) They are stored as Int32 and can
## have missing values.  The DataArray method must copy the values if it is to produce Bool's.
## For the time being, I will leave them as Int32's.
for N in [LGLSXP,INTSXP,REALSXP,CPLXSXP,STRSXP]
    @eval begin
        DataArrays.DataArray(s::SEXP{$N}) = (rv = reshape(vec(s),size(s));DataArray(rv,isNA(rv)))
        dataset(s::SEXP{$N}) = DataArray(s)
    end
end

## ToDo: expand these definitions to check for names and produce NamedArrays?

@doc "dataset method for SEXP{INTSXP} must check if the argument is a factor"->
function dataset(s::SEXP{INTSXP})
    isFactor(s) || return DataArray(s)
    ## refs array uses a zero index where R has a missing value, R_NaInt
    refs = DataArrays.RefArray(map!(x -> x == R_NaInt ? zero(Int32) : x,vec(s)))
    compact(PooledDataArray(refs,R.levels(s)))
end

Base.names(s::SEXP) = vec(sexp(ccall((:Rf_getAttrib,libR),Ptr{Void},
                                       (Ptr{Void},Ptr{Void}),s,namesSymbol)))

function dataset(s::SEXP{VECSXP})
    val = [dataset(v) for v in s]
    R.inherits(s,"data.frame") ? DataFrame(val,Symbol[symbol(nm) for nm in names(s)]) : val
end

dataset(st::ASCIIString) = dataset(reval(rparse(st)))

@doc "Evaluate Symbol s as an R dataset"->
dataset(s::Symbol) = dataset(reval(s))

@doc "extract the value of symbol s in the environment e"->
Base.getindex(e::SEXP{ENVSXP},s::Symbol) =
    sexp(ccall((:Rf_findVarInFrame,libR),Ptr{Void},(Ptr{Void},Ptr{Void}),e,sexp(s)))

@doc "assign value v to symbol s in the environment e"->
Base.setindex!(e::SEXP{ENVSXP},v::SEXP,s::Symbol) =
    ccall((:Rf_setVar,libR),Void,(Ptr{Void},Ptr{Void},Ptr{Void}),sexp(s),v,e)
Base.setindex!{T<:Number}(e::SEXP{ENVSXP},v::Array{T},s::Symbol) = setindex!(e,sexp(v),s)
Base.setindex!(e::SEXP{ENVSXP},v::Number,s::Symbol) = setindex!(e,sexp(v),s)
Base.setindex!(e::SEXP{ENVSXP},v::Real,s::Symbol) = setindex!(e,sexp(v),s)

function preserve(s::SEXP)
    ccall((:R_PreserveObject,libR),Void,(Ptr{Void},),s)
    finalizer(s,x -> ccall((:R_ReleaseObject,libR),Void,(Ptr{Void},),x))
    s
end

function sexp(v::BitVector)             # handled separately
    l = length(v)
    vv = sexp(ccall((:Rf_allocVector,libR),Ptr{Void},(Cint,Cptrdiff_t),LGLSXP,l))
    copy!(pointer_to_array(convert(Ptr{Int32},vv.p+voffset),l),v)
    preserve(vv)
end

for (typ,rnm,tag,rtyp) in ((:Bool,:Logical,LGLSXP,:Int32),
                           (:Complex128,:Complex,CPLXSXP,:Complex128),
                           (:Integer,:Integer,INTSXP,:Int32),
                           (:Real,:Real,REALSXP,:Float64))
    @eval begin
        function sexp(v::$typ)
            preserve(sexp(ccall(($(string("Rf_Scalar",rnm)),libR),Ptr{Void},($rtyp,),v)))
        end
        function sexp{T<:$typ}(v::Vector{T})
            l = length(v)
            vv = sexp(ccall((:Rf_allocVector,libR),Ptr{Void},(Cint,Cptrdiff_t),$tag,l))
            copy!(pointer_to_array(convert(Ptr{$rtyp},vv.p+voffset),l),v)
            preserve(vv)
        end
        function sexp{T<:$typ}(m::Matrix{T})
            p,q = size(m)
            vv = sexp(ccall((:Rf_allocMatrix,libR),Ptr{Void},(Cint,Cint,Cint),$tag,p,q))
            copy!(pointer_to_array(convert(Ptr{$rtyp},vv.p+voffset),l),m)
            preserve(vv)
        end
        function sexp{T<:$typ}(a::Array{T,3})
            p,q,r = size(a)
            vv = sexp(ccall((:Rf_alloc3DArray,libR),Ptr{Void},(Cint,Cint,Cint,Cint),$tag,p,q,r))
            copy!(pointer_to_array(convert(Ptr{$rtyp},vv.p+voffset),l),a)
            preserve(vv)
        end
        function sexp{T<:$typ}(a::Array{T})
            rdims = sexp([size(a)...])
            vv = sexp(ccall((:Rf_allocArray,libR),Ptr{Void},(Cint,Ptr{Void}),$tag,rdims))
            copy!(pointer_to_array(convert(Ptr{$rtyp},vv.p+voffset),l),a)
            preserve(vv)
        end
    end
end

sexp(s::Symbol) = sexp(ccall((:Rf_install,libR),Ptr{Void},(Ptr{Uint8},),string(s)))
sexp(st::Union(ASCIIString,UTF8String)) = sexp(ccall((:Rf_mkString,libR),Ptr{Void},(Ptr{Uint8},),st))


text:
using Plots;gr()

T = 1000
dt = 0.001
y_0 = 1

# オイラー法
function euler_method(y_0, T, dt)
    y = zeros(T)
    y[1] = y_0

    for t = 1:T-1
        y[t+1] = (1+5*dt)y[t] + 4*dt
    end

    return y
end

function exact_solution(t)
    (y_0+4//5)*exp.(5*t) - 4//5
end

t = 0:dt:(T-1)dt
y = euler_method(y_0, T, dt)

plot(t, exact_solution, color=:cornflowerblue, label="Exact",xlabel="t",ylabel="y")
plot!(t, y, color=:coral, linestyle=:dash, label="Numerical",xlabel="t",ylabel="y")
savefig("1st_order_euler.png")


text:
module LieAlgebraRankCondition

	using LinearAlgebra

	export lieprod,isli,rankcondition,ad, OutputSet

	import Base.show

	"""
		OutputSet

	It is an internal type used to print the output of `rankcondition` function. Once the `rankcondition` function has been used, just access the fields according to the example below.

	# Examples
	```
	julia-repl
	julia> s = rankcondition(A,B)
	julia> s.LI # returns true or false
	julia> s.elements # returns a sequence of arrays (L.I arrays)
	julia> s.path # show the sequence of Lie products used to get the solution `s`.
	```
	"""
	struct OutputSet
		path :: Vector{String}
		LI :: Bool
		elements :: Vector
	end

	"""
		lieprod

	An auxiliary function to compute the lie product of two matrix: AB-BA

	# Examples
	```
	julia-repl
	julia> A = rand(4,4)
	julia> B = rand(4,4)
	julia> C = lieprod(A,b)
	```
	returns the lie product of A and B and storage in C.
	"""
	function lieprod(A::Array{Float64},B::Array{Float64})
		return A*B-B*A
	end

	"""
		isli

	An auxiliary function to test if a set of matrix is an L. I set.

	# Examples
	```
	julia-repl
	julia> C=[A,B]
	julia> t = isli(C)
	```
	returns true or false and storage the answer in t.
	"""
	function isli(C::Vector)
		D = zeros(length(C),16)
		for k=1:length(C)
			D[k,:]=C[k]
		end
		if rank(D)==length(C)
			return true
		else
			return false
		end
	end
	"""
		rankcondition

	It is the main function. See the example to run.

	# Example
	```
	julia-repl
	julia> A = [1 1 0 0; -1.0 1 0 0;0 0 -1 0.5;0 0 -0.5 -1]
	julia> B = [2.0 0.0 0.0 0.0; 0.0 -1.5 -0.1 0.0; 0.0 0.1 -1.5 0.0; 0.0 0.0 0.0 1.0]
	julia> s = rankcondition(A,B)
	```
	returns s like an OutputSet with three fields. This OutputSet type is used to store L.I. matrices, bool values and a convenient path to reproduce the same solution. There is an optional argument used to define the number of L.I. elements.
	# Example
	```
	julia-repl
	julia> A = [1 1 0 0; -1.0 1 0 0;0 0 -1 0.5;0 0 -0.5 -1]
	julia> B = [ 2.0   0.0   0.0  0.0; 0.0  -1.5  -0.1  0.0; 0.0   0.1  -1.5  0.0; 0.0   0.0   0.0  1.0]
	julia> s = rankcondition(A,B,dim=15)
	```
	By default `dim` parameter can be omitted. You can change this parameter if it is convenient for you.
	"""
	function rankcondition(A,B,dim=15)
		partial = [A,B]
		path = ["A","B"]
		search_index = ones(Int,dim)
		if isli(partial)==true
			push!(partial,lieprod(partial[1],partial[2]))
			push!(path,"[$(path[1]),$(path[2])]")
			k = 3
			while 3<=k<=dim
				next_step = false
				search_partial_solution = true
				while search_partial_solution == true
					if isli(partial)==true
						push!(partial,lieprod(partial[search_index[k]],partial[k]))
						push!(path,"[$(path[search_index[k]]),$(path[k])]")
						next_step = true
						search_partial_solution = false
					elseif search_index[k]<=k-1
						search_index[k] = search_index[k]+1
						partial[k] = lieprod(partial[search_index[k]],partial[k-1])
						path[k] = "[$(path[search_index[k]]),$(path[k-1])]"
						search_partial_solution = false
					end
				end
				if next_step == false
					deleteat!(partial,k)
					deleteat!(path,k)
					search_index[k] = 1
					k = k-1
					search_index[k]=search_index[k]+1
					partial[k] = lieprod(partial[search_index[k]],partial[k-1])
					path[k] = "[$(path[search_index[k]]),$(path[k-1])]"
				else
					k = k+1
				end
			end
			if k<=3
				println(" 🚑 No L.I sets were found!")
				return OutputSet([""],false,[])
			end
			if isli(partial[1:dim]) == true
				println(" 🎈 A set of L.I arrays with $(dim) elements was found.")
				return OutputSet(path[1:dim],true,partial[1:dim])
			end
		else
			println(" 🚑 No L.I sets with $(dim) elements were found!")
			return OutputSet([""],false,[])

		end
	end

	"""
		ad(n,A,B)

		This functins is used to compute [A[A,...,[A,B]]

	# Examples
	```
	julia-repl
	julia> ad(2,A,B)
	```
	returns lieprod(A,lieprod(A,B))
	"""
	function ad(n,M1,M2)
		M = lieprod(M1,M2)
		for k=2:n
			M = lieprod(M1,M)
		end
		return M
	end

	#####################
	# print output	    #			                                                                    #
	#####################

	function Base.show(io::IO,ots::OutputSet)
		if ots.LI == true
			print(io,"\n 👉 Here are the set and the path:  \n")
			for k = 1:length(ots.elements)
				print(io, "\n ↪  Element $(k):")
				print(io, " $(ots.path[k]) = \n")
				display(ots.elements[k])
			end
		else
			print(io, "\n ")
		end
	end
end # module


text:
@recipe function f{CT <: AbstractPredictiveAsymmetryTest, DT, RT}(analysis::CausalAnalysis{CT,DT,RT})
    pas = analysis.result
    ymax = maximum(pas) * 1.1
    ηs = get_ηs(analysis.test)
    ηs = ηs[ηs .> 0]

    if CT isa NormalisedPredictiveAsymmetryTest
        ylabel --> L"$\mathbb{A}(\eta)$"
    elseif CT isa PredictiveAsymmetryTest
        ylabel --> L"$\mathcal{A}(\eta)$"
    end
    seriestype --> :line
    xlabel --> L"$\eta$"
    ylabel --> L"$\mathbb{A}(\eta)$"
    label --> ""
    size --> (382, 300)
    fg_legend --> :transparent
    bg_legend --> :transparent

    (ηs, pas)
end


text:
abstract type MomentumBC end
abstract type MomentumDragBC end

"""
    Impenetrable(drag::MomentumDragBC) :: MomentumBC

Defines an impenetrable wall model for momentum. This implies:
  - no flow in the direction normal to the boundary, and
  - flow parallel to the boundary is subject to the `drag` condition.
"""
struct Impenetrable{D <: MomentumDragBC} <: MomentumBC
    drag::D
end

"""
    FreeSlip() :: MomentumDragBC

No surface drag on momentum parallel to the boundary.
"""
struct FreeSlip <: MomentumDragBC end

function atmos_momentum_boundary_state!(
    nf::NumericalFluxFirstOrder,
    bc_momentum::Impenetrable{FreeSlip},
    atmos,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    bctype,
    t,
    args...,
)
    state⁺.ρu -= 2 * dot(state⁻.ρu, n) .* SVector(n)
end
function atmos_momentum_boundary_state!(
    nf::NumericalFluxGradient,
    bc_momentum::Impenetrable{FreeSlip},
    atmos,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    bctype,
    t,
    args...,
)
    state⁺.ρu -= dot(state⁻.ρu, n) .* SVector(n)
end
function atmos_momentum_normal_boundary_flux_second_order!(
    nf,
    bc_momentum::Impenetrable{FreeSlip},
    atmos,
    args...,
) end

"""
    NoSlip() :: MomentumDragBC

Zero momentum at the boundary.
"""
struct NoSlip <: MomentumDragBC end

function atmos_momentum_boundary_state!(
    nf::NumericalFluxFirstOrder,
    bc_momentum::Impenetrable{NoSlip},
    atmos,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    bctype,
    t,
    args...,
)
    state⁺.ρu = -state⁻.ρu
end
function atmos_momentum_boundary_state!(
    nf::NumericalFluxGradient,
    bc_momentum::Impenetrable{NoSlip},
    atmos,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    bctype,
    t,
    args...,
)
    state⁺.ρu = zero(state⁺.ρu)
end
function atmos_momentum_normal_boundary_flux_second_order!(
    nf,
    bc_momentum::Impenetrable{NoSlip},
    atmos,
    args...,
) end

"""
    DragLaw(fn) :: MomentumDragBC

Drag law for momentum parallel to the boundary. The drag coefficient is
`C = fn(state, aux, t, normPu_int)`, where `normPu_int⁻` is the internal speed
parallel to the boundary.
"""
struct DragLaw{FN} <: MomentumDragBC
    fn::FN
end
function atmos_momentum_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    bc_momentum::Impenetrable{DL},
    atmos,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    bctype,
    t,
    args...,
) where {DL <: DragLaw}
    atmos_momentum_boundary_state!(
        nf,
        Impenetrable(FreeSlip()),
        atmos,
        state⁺,
        aux⁺,
        n,
        state⁻,
        aux⁻,
        bctype,
        t,
        args...,
    )
end
function atmos_momentum_normal_boundary_flux_second_order!(
    nf,
    bc_momentum::Impenetrable{DL},
    atmos,
    fluxᵀn,
    n,
    state⁻,
    diff⁻,
    hyperdiff⁻,
    aux⁻,
    state⁺,
    diff⁺,
    hyperdiff⁺,
    aux⁺,
    bctype,
    t,
    state1⁻,
    diff1⁻,
    aux1⁻,
) where {DL <: DragLaw}

    u1⁻ = state1⁻.ρu / state1⁻.ρ
    Pu1⁻ = u1⁻ - dot(u1⁻, n) .* SVector(n)
    normPu1⁻ = norm(Pu1⁻)
    # NOTE: difference from design docs since normal points outwards
    C = bc_momentum.drag.fn(state⁻, aux⁻, t, normPu1⁻)
    τn = C * normPu1⁻ * Pu1⁻
    # both sides involve projections of normals, so signs are consistent
    fluxᵀn.ρu += state⁻.ρ * τn
    fluxᵀn.ρe += state⁻.ρu' * τn
end


text:
using Luxor, Test, Random

# test empty drawings

Drawing()
[Drawing() for x in 1:10]
[Drawing(x, x) for x in 1:10]
Drawing(NaN, NaN, :rec)
[Drawing(x, x, :rec) for x in 1:10]
[Drawing(x, x, :image) for x in 1:10]

p = Point(1, 2)
x, y = p

@test x == 1.0
@test y == 2.0

bb = BoundingBox(p, 2p)

rp = rand(bb)

@test rp in bb

@testset "center3pts" begin
    @test center3pts(Point(1,1), Point(3,1), Point(2,0)) == (Point(2.0, 1.0), 1.0)
    Random.seed!(42)
    p1 = Point(rand(), rand())
    p2 = Point(rand(), rand())
    p3 = Point(rand(), rand())
    c, r = center3pts(p1, p2, p3)
    @test distance(p1, c) ≈ distance(p2, c) ≈ distance(p2, c)
end

println("...finished basic-test")


text:
module ODEConvergenceTester

import LinearAlgebra
import Logging
import OrdinaryDiffEq

compute_err_norm_default(u::FT, v::FT) where {FT <: Real} = LinearAlgebra.norm(u - v)
compute_err_norm_default(u, v) = LinearAlgebra.norm(u .- v)

Base.@kwdef struct ConvergenceTestResults{FT}
    p::Vector{FT}
    err::Vector{FT}
    dts::Vector{FT}
    p_expected::FT
end

"""
    refinement_study(
        integrator₀,
        dts::Vector{<:Real};
        compute_err_norm::Function = compute_err_norm_default,
    )

    refinement_study(
        integrator₀;
        compute_err_norm::Function = compute_err_norm_default,
        refinement_range::UnitRange = 1:3,
    )

Estimates and reports the convergence rates.

## Arguments
 - `integrator₀` a OrdinaryDiffEq integrator
 - `dts` a vector of timesteps
 - `compute_err_norm = (x,y) -> norm(x,y)` a function for computing
    the norm of the solution
 - `refinement_range = 1:3` a `UnitRange` of refinements.
   `2:4` is finer than `1:3`

## Theory /  analysis

Consider the error for a time-stepping scheme

```
    err = C Δt^p+ H.O.T
    err_k ≈ C Δt_k^p
    err_k/err_m ≈ Δt_k^p/Δt_m^p
    log(err_k/err_m) ≈ log((Δt_k/Δt_m)^p)
    log(err_k/err_m) ≈ p*log(Δt_k/Δt_m)
    log(err_k/err_m)/log(Δt_k/Δt_m) ≈ p
```

Or, put differently, for `refinement_range = 1:3`:

```
            |f3 - f2|   /
    p = log ---------  / log (rf)
            |f2 - f1| /
```

Where
 - `f1` is the finest resolution
 - `f3` is the coarsest resolution
 - `rf` (where `rf` > 1) is the refinement factor

## References

    Roache, P. J. Quantification of Uncertainty in Computational
    Fluid Dynamics. Annu. Rev. Fluid Mech. 29, 123–160 (1997).
    De Vahl Davis, G. Natural convection of air in a square cavity: a
    benchmark solution. Int. J. Num. Methods Fluids 3, 249–264 (1983).

    http://www.grc.nasa.gov/WWW/wind/valid/tutorial/spatconv.html
"""
function refinement_study(
        integrator₀,
        dts::Vector{dtType};
        compute_err_norm::Function = compute_err_norm_default,
        print_report::Bool = true,
    ) where {dtType <: Real}

    n_refinements = length(dts)
    # This can be pretty expensive, depending on the problem
    # so let's limit the number of refinements
    @assert 3 ≤ n_refinements "Need at least 3 runs to compute convergence order"

    # Create deepcopy of integrators to prevent accidental mutation
    integrators = map(1:n_refinements) do n
        deepcopy(integrator₀)
    end

    # Refinement factors
    refinement_factor = map(1:(n_refinements - 1)) do i
        dts[i] / dts[i + 1]
    end

    # Compute new t_final such that n_steps[i]*dt[i] = t_final
    t_final₀ = last(integrator₀.sol.prob.tspan)
    n_steps₀ = round(t_final₀ / dts[1])
    n_steps = map(i -> n_steps₀ * 2^(i - 1), 1:n_refinements)
    t_final = map(enumerate(dts)) do (i, dt)
        n_steps[i] * dt
    end

    if print_report
        @info "------ Convergence parameters ------"
        @info "nsteps                    : $n_steps"
        @info "refinement factors        : $refinement_factor"
        @info "dts                       : $dts"
        @info "tfinal (rounded)          : $t_final"
        @info "--- Running convergence study... ---"
    end

    for i in 1:n_refinements
        dt = dts[i]
        OrdinaryDiffEq.set_proposed_dt!(integrators[i], dt)
        print("@timing iteration $i:")
        @time Logging.with_logger(Logging.NullLogger()) do
            for n in 1:n_steps[i]
                OrdinaryDiffEq.step!(integrators[i], dt)
            end
        end
        u = integrators[i].u
    end
    err = map(1:(n_refinements - 1)) do i
        u = deepcopy(integrators[i].u)
        v = deepcopy(integrators[i + 1].u)
        compute_err_norm(u, v)
    end
    p = map(1:(length(err) - 1)) do i
        log(err[i] / err[i + 1]) / log(max(refinement_factor...))
    end

    FT = eltype(p)
    p_expected = FT(OrdinaryDiffEq.alg_order(integrator₀.sol.alg))
    if print_report
        @info "----- Convergence study output -----"
        @info "Errors                    : $(err)"
        @info "Convergence orders        : $(p)"
        @info "Expected convergence order: $(p_expected)"
        @info "------------------------------------"
    end
    return ConvergenceTestResults(; p, err, dts, p_expected)
end

function refinement_study(
        integrator₀;
        refinement_range::UnitRange = 1:3,
        compute_err_norm::Function = compute_err_norm_default,
        print_report::Bool = true,
    )
    t_final₀ = last(integrator₀.sol.prob.tspan)
    # For this analysis, Δt must be ~ 2^n, so we must
    # find a scale that will result in a reasonable timestep:
    # 2^(scale) = t_final
    # i.e., dt_max should scale linearly with t_final, while
    # still being of the form 2^n
    scale = round(log(t_final₀) / log(2))
    buffer = 3 # use a buffer to avoid extremely coarse timesteps
    dt_max = 2^(scale - buffer)
    refinement_factor = 2
    dts = map(refinement_range) do i
        dt_max / (refinement_factor^(i - 1))
    end
    return refinement_study(
        integrator₀,
        dts;
        print_report=print_report,
        compute_err_norm=compute_err_norm,
    )
end

end # module


text:
@testset "Refinement" begin
  @testset "TriRefinement" begin
    grid = CartesianGrid{T}(3, 3)
    ref1 = refine(grid, TriRefinement())
    ref2 = refine(ref1, TriRefinement())

    if visualtests
      p1 = plot(grid, fillcolor=false)
      p2 = plot(ref1, fillcolor=false)
      p3 = plot(ref2, fillcolor=false)
      p = plot(p1, p2, p3, layout=(1,3), size=(900,300))
      @test_reference "data/trirefine-$T.png" p
    end
  end

  @testset "QuadRefinement" begin
    points = P2[(0,0), (1,0), (0,1), (1,1), (0.25,0.25), (0.75,0.25), (0.5,0.75)]
    connec = connect.([(1,2,6,5),(1,5,7,3),(2,4,7,6),(3,7,4)])
    mesh   = SimpleMesh(points, connec)
    ref1   = refine(mesh, QuadRefinement())
    ref2   = refine(ref1, QuadRefinement())
    ref3   = refine(ref2, QuadRefinement())

    if visualtests
      p1 = plot(ref1, fillcolor=false)
      p2 = plot(ref2, fillcolor=false)
      p3 = plot(ref3, fillcolor=false)
      p = plot(p1, p2, p3, layout=(1,3), size=(900,300))
      @test_reference "data/quadrefine-$T.png" p
    end
  end

  @testset "CatmullClark" begin
    points = P2[(0,0), (1,0), (0,1), (1,1), (0.5,0.5)]
    connec = connect.([(1,2,5),(2,4,5),(4,3,5),(3,1,5)])
    mesh   = SimpleMesh(points, connec)
    ref1   = refine(mesh, CatmullClark())
    ref2   = refine(ref1, CatmullClark())
    ref3   = refine(ref2, CatmullClark())

    if visualtests
      p1 = plot(ref1, fillcolor=false)
      p2 = plot(ref2, fillcolor=false)
      p3 = plot(ref3, fillcolor=false)
      p = plot(p1, p2, p3, layout=(1,3), size=(900,300))
      @test_reference "data/catmullclark-1-$T.png" p
    end

    points = P2[(0,0), (1,0), (0,1), (1,1), (0.25,0.25), (0.75,0.25), (0.5,0.75)]
    connec = connect.([(1,2,6,5),(1,5,7,3),(2,4,7,6),(3,7,4)])
    mesh   = SimpleMesh(points, connec)
    ref1   = refine(mesh, CatmullClark())
    ref2   = refine(ref1, CatmullClark())
    ref3   = refine(ref2, CatmullClark())

    if visualtests
      p1 = plot(ref1, fillcolor=false)
      p2 = plot(ref2, fillcolor=false)
      p3 = plot(ref3, fillcolor=false)
      p = plot(p1, p2, p3, layout=(1,3), size=(900,300))
      @test_reference "data/catmullclark-2-$T.png" p
    end

    points = P3[(0,0,0),(1,0,0),(1,1,0),(0,1,0),(0,0,1),(1,0,1),(1,1,1),(0,1,1)]
    connec = connect.([(1,4,3,2),(5,6,7,8),(1,2,6,5),(3,4,8,7),(1,5,8,4),(2,3,7,6)])
    mesh   = SimpleMesh(points, connec)
    ref1   = refine(mesh, CatmullClark())
    ref2   = refine(ref1, CatmullClark())
    ref3   = refine(ref2, CatmullClark())

    if visualtests
      p1 = plot(ref1, fillcolor=false)
      p2 = plot(ref2, fillcolor=false)
      p3 = plot(ref3, fillcolor=false)
      p = plot(p1, p2, p3, layout=(1,3), size=(900,300))
      @test_reference "data/catmullclark-3-$T.png" p
    end
  end
end


text:
export NaiveNPU, RealNaiveNPU, RealNPU, NPU

"""
  NPU(in::Int, out::Int; initRe=glorot_uniform, initIm=Flux.zeros, initg=init05)

Neural Power Unit that can learn arbitrary power functions by using a complex
weights. Uses gating on inputs to simplify learning. In 1D the layer looks
like:

    g = min(max(g, 0), 1)
    r = abs(x) + eps(T)
    r = g*r + (1-g)*T(1)
    k = r < 0 ? pi : 0.0
    exp(Re*log(r) - Im*k) * cos(Re*k + Im*log(r))
"""
struct NPU{M<:AbstractMatrix,V<:AbstractVector}
    Re::M
    Im::M
    g::V
end

function NPU(in::Int, out::Int;
    initRe=glorot_uniform, initIm=Flux.zeros, initg=init05)
    Re = initRe(out, in)
    Im = initIm(out, in)
    g  = initg(in)
    NPU(Re,Im,g)
end

Flux.@functor NPU

function mult(Re::AbstractMatrix{T}, Im::AbstractMatrix{T}, g::AbstractVector{T}, x::AbstractArray{T}) where T
    g = gateclip(g)

    r = abs.(x) .+ eps(T)
    r = g .* r .+ (1 .- g)

    k = signk(x)
    k = g .* k

    exp.(Re*log.(r) - Im*k) .* cos.(Re*k + Im*log.(r))
end
(l::NPU)(x) = mult(l.Re, l.Im, l.g, x)

"""
    RealNPU(in::Int, out::Int; initRe=glorot_uniform, initg=init05)

NPU without imaginary weights.
"""
struct RealNPU{Tw<:AbstractMatrix,Tg<:AbstractVector}
    Re::Tw
    g::Tg
end

RealNPU(in::Int, out::Int; initRe=Flux.glorot_uniform, initg=init05) =
    RealNPU(initRe(out,in), initg(in))

Flux.@functor RealNPU

function mult(Re::AbstractMatrix{T}, g::AbstractVector{T}, x::AbstractArray{T}) where T
    g = gateclip(g)

    r = abs.(x) .+ eps(T)
    r = g .* r .+ (1 .- g)

    k = signk(x)
    k = g .* k

    z = exp.(Re * log.(r)) .* cos.(Re*k)
end

(l::RealNPU)(x) = mult(l.Re, l.g, x)

"""
    NaiveNPU(in::Int, out::Int; initRe=glorot_uniform, initIm=zeros)

`NPU` without relevance gating mechanism.
"""
struct NaiveNPU{T<:AbstractMatrix}
    Re::T
    Im::T
end

function NaiveNPU(in::Int, out::Int; initRe=glorot_uniform, initIm=Flux.zeros)
    Re = initRe(out, in)
    Im = initIm(out, in)
    NaiveNPU(Re, Im)
end

Flux.@functor NaiveNPU

function mult(Re::AbstractMatrix{T}, Im::AbstractMatrix{T}, x::AbstractArray{T}) where T
    r = abs.(x) .+ eps(T)
    k = signk(x)
    exp.(Re*log.(r) - Im*k) .* cos.(Re*k + Im*log.(r))
end

(l::NaiveNPU)(x) = mult(l.Re, l.Im, x)

"""
  RealNaiveNPU(in::Int, out::Int; initRe=glorot_uniform)

NaiveNPU without imaginary weights.
"""
struct RealNaiveNPU{T<:AbstractMatrix}
    Re::T
end

RealNaiveNPU(in::Int, out::Int; initRe=glorot_uniform) = RealNaiveNPU(initRe(out,in))
Flux.@functor RealNaiveNPU

function mult(Re::AbstractMatrix{T}, x::AbstractArray{T}) where T
    r = abs.(x) .+ eps(T)
    k = signk(x)
    exp.(Re * log.(r)) .* cos.(Re*k)
end

(l::RealNaiveNPU)(x) = mult(l.Re, x)

init05(s...) = Flux.ones(s...) ./ 2
signk(x::AbstractArray{T}) where T = max.(-sign.(x), 0) * T(π)
gateclip(g::AbstractVector) = min.(max.(g, 0), 1)

function ChainRulesCore.rrule(::typeof(gateclip), g)
    ĝ = gateclip(g)
    gateclip_pullback(ΔΩ) = (NoTangent(), ΔΩ)
    return ĝ, gateclip_pullback
end

function _repr(l)
    T = string(nameof(typeof(l)))
    "$T(in=$(size(l.Re,2)), out=$(size(l.Re,1)))"
end
Base.show(io::IO, l::NPU) = print(io, _repr(l))
Base.show(io::IO, l::NaiveNPU) = print(io, _repr(l))
Base.show(io::IO, l::RealNPU) = print(io, _repr(l))
Base.show(io::IO, l::RealNaiveNPU) = print(io, _repr(l))


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

code = create_basis(
    :Seg2,
    "2 node linear segment/line element",
    ( (-1.0,), ( 1.0,) ),
    :(1 + u))
eval(code)

code = create_basis(
    :Seg3,
    "3 node quadratic segment/line element",
    ( (-1.0,), (1.0,), (0.0,)),
    :(1 + u + u^2))
eval(code)


text:
using IndianTrafficSim
using Test

@testset "IndianTrafficSim.jl" begin
    p = plot_environment!()
    @test true == true
end


text:
@deprecate seq2ref(i::Integer, aln::Alignment) seq2ref(aln::Alignment, i::Integer)
@deprecate ref2seq(i::Integer, aln::Alignment) ref2seq(aln::Alignment, i::Integer)
@deprecate seq2ref(i::Integer, alnseq::AlignedSequence) seq2ref(aln::AlignedSequence, i::Integer)
@deprecate ref2seq(i::Integer, alnseq::AlignedSequence) ref2seq(aln::AlignedSequence, i::Integer)
@deprecate seq2ref(i::Integer, aln::PairwiseAlignment) seq2ref(aln::PairwiseAlignment, i::Integer)
@deprecate ref2seq(i::Integer, aln::PairwiseAlignment) ref2seq(aln::PairwiseAlignment, i::Integer)


text:
using HDF5
#using LegendHDF5IO: readdata, writedata
using SolidStateDetectors
using Unitful

using Plots;
pyplot(fmt = :png)

T = Float32

function SaveSimulation(filename::String, simulation::Simulation{T})::Nothing
    @info "Saving simulation"
    if !isfile(filename) mkpath(dirname(filename)) end
    HDF5.h5open(filename, "w") do h5f
        writedata(h5f, "Simulation", NamedTuple(sim))
    end
end

function ReadSimulation(filename)::Union{Nothing, Simulation{T}}
    if isfile(filename)
        @info "Reading simulation"
        simulation = HDF5.h5open(filename, "r") do h5f
            Simulation(readdata(h5f, "Simulation"))
        end
        set_charge_drift_model!(simulation, ADLChargeDriftModel())
        return simulation
    end
    return nothing
end

function CalculateDefaultDetectorFields!(sim::Simulation{T},
    driftConfigFile::Union{Missing, String} = missing,
    max_ref = 4)

    CalculateElectricPotential!(sim, max_ref)
    CalculateElectricField!(sim)
    CalculateWeightingPotential!(sim, max_ref)
    SetChargeDriftModel!(sim, driftConfigFile)
end

function SetupDefaultSimulation(geomConfigFile::String, CCDName::String="")::Simulation{T}
    simulation = Simulation{T}(geomConfigFile)

    if CCDName != ""
        include(CCDName)
        simulation.detector.semiconductors[1].charge_density_model = ccdm
    end

    apply_initial_state!(simulation, ElectricPotential)
    for c in simulation.detector.contacts
        apply_initial_state!(simulation, WeightingPotential, c.id)
    end
    return simulation
end

function SetBiasVoltage!(sim, contact, voltage)
    sim.detector.contacts[contact].potential = voltage
end

function CalculateElectricPotential!(simulation::Simulation{T}, max_ref = 4, verbose = false)
    calculate_electric_potential!(simulation, depletion_handling = true, max_refinements = max_ref)
    if verbose
        println("Active volume: ", get_active_volume(simulation.point_types))
    end
end

function CalculateWeightingPotential!(simulation::Simulation{T}, max_ref = 2)
    # TODO Add options
    for contact in simulation.detector.contacts
        calculate_weighting_potential!(simulation, contact.id, max_refinements = max_ref, n_points_in_φ = 2, verbose = false)
    end
end

function CalculateElectricField!(simulation::Simulation{T}, points_in_φ = 72, verbose = true)
    calculate_electric_field!(simulation, n_points_in_φ = points_in_φ)
    if verbose
        println("Capacitance: ", calculate_capacitance(simulation))
    end
end

function SetChargeDriftModel!(simulation::Simulation{T}, driftConfigFile::Union{Missing, String})
    charge_drift_model = ADLChargeDriftModel(driftConfigFile, T=T)
    set_charge_drift_model!(simulation, charge_drift_model)
    calculate_drift_fields!(simulation)
end

#-------------------------------------------------------------------------------
# Plotting
#

function PlotDetector(simulation::Simulation{T})
    p = plot(simulation.detector, size = (700, 700), title = "Detector geometry")
    return p
end

function PlotElectricPotential(simulation::Simulation{T})
    p = plot(
    plot(simulation.electric_potential), # initial electric potential (boundary conditions)
    plot(simulation.point_types), # map of different point types: fixed point / inside or outside detector volume / depleted/undepleted
    plot(simulation.ρ), # charge density distribution
    plot(simulation.ϵ), # dielectric distribution
    layout = (1, 4), size = (1400, 700)
    )
    return p
end

function PlotWeightingPotential(simulation::Simulation{T}, contact::Int64)
    p = plot(
    plot(simulation.weighting_potentials[contact]),
    plot(simulation.point_types),
    plot(simulation.electric_potential),
    layout = (1, 3), size = (1000, 700)
    )
    return p
end

max_stars_repo_name: rjtayl/NDP

max_issues_repo_name: rjtayl/NDP

max_forks_repo_name: rjtayl/NDP


text:
import DataFrames

"""
"""
function predictionsassoctodataframe(
        probabilitiesassoc::AbstractDict,
        label_names::AbstractVector = [],
        )
    if length(label_names) == 0
        label_names = sort(unique(collect(keys(probabilitiesassoc))))
    end
    result = DataFrames.DataFrame()
    for j = 1:length(label_names)
        result[label_names[j]] = probabilitiesassoc[label_names[j]]
    end
    return result
end


text:
#### Cauchy matrix

export Cauchy

struct Cauchy{T} <: AbstractMatrix{T}
    x::Vector{T} #
    y::Vector{T} #
end # immutable

function Cauchy(k::Number)
         Cauchy(collect(1:k),collect(1:k))
end

function Cauchy(x::Vector)
         Cauchy(x,x)
end

# Define its size

size(A::Cauchy, dim::Integer) = length(A.x)
size(A::Cauchy)= size(A,1), size(A,1)

# Index into a Cauchy
function getindex(A::Cauchy,i::Integer,j::Integer)
    return 1.0/(A.x[i]+A.y[j])
end # getindex

# Dense version of Cauchy
Matrix(A::Cauchy) = [A[i,j] for i=1:size(A,1), j=1:size(A,2)]


text:
type MockResource
  resource :: Vector{Int}
  rc       :: RefCounter
end

function alloc_resource(n::Int)
  MockResource(zeros(Int, n), RefCounter(1))
end
function share_resource(r :: MockResource)
  MockResource(r.resource, ref(r.rc))
end
function free_resource(r :: MockResource)
  if dec(r.rc) == 0
    fill!(r.resource, -1)
  end
end

function test_ref_count()
  println("-- Testing simple reference counting...")

  res = alloc_resource(1)
  res2 = share_resource(res)

  free_resource(res)
  @test res2.resource[1] != -1

  free_resource(res2)
  @test res2.resource[1] == -1
end

test_ref_count()


text:
import Base: copy!
import Base: splice!
import Base: append!
import Base: push!
import Base: resize!
import Base: setindex!
import Base: getindex
import Base: map
import Base: length
import Base: eltype
import Base: ndims
import Base: size

abstract type GPUArray{T, NDim} <: AbstractArray{T, NDim} end

length(A::GPUArray) = prod(size(A))
eltype(b::GPUArray{T, NDim}) where {T, NDim} = T
endof(A::GPUArray) = length(A)
ndims(A::GPUArray{T, NDim}) where {T, NDim} = NDim
size(A::GPUArray) = (A.len,)
size(A::GPUArray, i::Integer) = i <= ndims(A) ? A.size[i] : 1

function checkdimensions(value::Array, ranges::Union{Integer, UnitRange}...)
    array_size   = size(value)
    indexes_size = map(length, ranges)
    (array_size != indexes_size) && throw(DimensionMismatch("asigning a $array_size to a $(indexes_size) location"))
    true
end
function to_range(index)
    map(index) do val
        isa(val, Integer) && return val:val
        isa(val, Range) && return val
        @error "Indexing only defined for integers or ranges. Found: $val"
    end
end
setindex!(A::GPUArray{T, N}, value::Union{T, Array{T, N}}) where {T, N} = (A[1] = value)

function setindex!(A::GPUArray{T, N}, value, indexes...) where {T, N}
    ranges = to_range(indexes)
    v = isa(value, T) ? [value] : convert(Array{T,N}, value)
    setindex!(A, v, ranges...)
    nothing
end

setindex!(A::GPUArray{T, 2}, value::Vector{T}, i::Integer, range::UnitRange) where {T} =
   (A[i, range] = reshape(value, (length(value),1)))

function setindex!(A::GPUArray{T, N}, value::Array{T, N}, ranges::UnitRange...) where {T, N}
    checkbounds(A, ranges...)
    checkdimensions(value, ranges...)
    gpu_setindex!(A, value, ranges...)
    nothing
end

function update!(A::GPUArray{T, N}, value::Array{T, N}) where {T, N}
    if length(A) != length(value)
        if isa(A, Buffer)
            resize!(A, length(value))
        elseif isa(A, Texture) && ndims(A) == 2
            resize_nocopy!(A, size(value))
        else
            @error "Dynamic resizing not implemented for $(typeof(A))"
        end
    end
    dims = map(x->1:x, size(A))
    A[dims...] = value
    nothing
end

function getindex(A::GPUArray{T, N}, i::Int) where {T, N}
    checkbounds(A, i)
    gpu_getindex(A, i:i)[1] # not as bad as its looks, as so far gpu data must be loaded into an array anyways
end
function getindex(A::GPUArray{T, N}, ranges::UnitRange...) where {T, N}
    checkbounds(A, ranges...)
    gpu_getindex(A, ranges...)
end

resize!(A::GPUArray{T, NDim}, dims::Int...) where {T, NDim} = resize!(A, dims)
function resize!(A::GPUArray{T, NDim}, newdims::NTuple{NDim, Int}) where {T, NDim}
    newdims == size(A) && return A
    gpu_resize!(A, newdims)
    A
end

copy!(a::GPUArray, a_offset::Int, b::Vector,   b_offset::Int, amount::Int) = _copy!(a, a_offset, b, b_offset, amount)
copy!(a::Vector,   a_offset::Int, b::GPUArray, b_offset::Int, amount::Int) = _copy!(a, a_offset, b, b_offset, amount)
copy!(a::GPUArray, a_offset::Int, b::GPUArray, b_offset::Int, amount::Int) = _copy!(a, a_offset, b, b_offset, amount)

#don't overwrite Base.copy! with a::Vector, b::Vector
function _copy!(a::Union{Vector, GPUArray}, a_offset::Int, b::Union{Vector, GPUArray}, b_offset::Int, amount::Int)
    (amount <= 0) && return nothing
    @assert a_offset > 0 && (a_offset-1) + amount <= length(a) "a_offset $a_offset, amount $amount, lengtha $(length(a))"
    @assert b_offset > 0 && (b_offset-1) + amount <= length(b) "b_offset $b_offset, amount $amount, lengthb $(length(b))"
    unsafe_copy!(a, a_offset, b, b_offset, amount)
    return nothing
end

# Interface:
gpu_data(t) = @error "gpu_data not implemented for: $(typeof(t)). This happens, when you call data on an array, without implementing the GPUArray interface"
gpu_resize!(t) = @error "gpu_resize! not implemented for: $(typeof(t)). This happens, when you call resize! on an array, without implementing the GPUArray interface"
gpu_getindex(t) = @error "gpu_getindex not implemented for: $(typeof(t)). This happens, when you call getindex on an array, without implementing the GPUArray interface"
gpu_setindex!(t) = @error "gpu_setindex! not implemented for: $(typeof(t)). This happens, when you call setindex! on an array, without implementing the GPUArray interface"
max_dim(t) = @error "max_dim not implemented for: $(typeof(t)). This happens, when you call setindex! on an array, without implementing the GPUArray interface"

# const BaseSerializer = if isdefined(Base, :AbstractSerializer)
#     Base.AbstractSerializer
# elseif isdefined(Base, :SerializationState)
#     Base.SerializationState
# else
#     error("No Serialization type found. Probably unsupported Julia version")
# end

# function Base.serialize(s::BaseSerializer, t::T) where T<:GPUArray
#     Base.serialize_type(s, T)
#     serialize(s, Array(t))
# end
# function Base.deserialize(s::BaseSerializer, ::Type{T}) where T<:GPUArray
#     A = deserialize(s)
#     T(A)
# end

export data
export resize
export GPUArray

export update!

export gpu_data
export gpu_resize!
export gpu_getindex
export gpu_setindex!
export max_dim


text:
###
# Error

immutable ArchiveRetry <: Exception
    msg
end
immutable ArchiveFailed <: Exception
    msg
end
immutable ArchiveFatal <: Exception
    msg
end

_la_error_msg(archive::Archive) = error_string(archive)
_la_error_msg(other) = ""

@noinline function _la_error(err::Cint, obj=nothing)
    err == Status.EOF && throw(EOFError())
    err == Status.RETRY && throw(ArchiveRetry(_la_error_msg(obj)))
    if err == Status.WARN
        warn("LibArchive: $(_la_error_msg(obj))")
        return
    end
    err == Status.FAILED && throw(ArchiveFailed(_la_error_msg(obj)))
    err == Status.FATAL && throw(ArchiveFatal(_la_error_msg(obj)))
    error("Unknown error $err")
end

macro _la_call(name, types, args...)
    call_ex = esc(:(ccall(($(QuoteNode(name)), $libarchive),
                          Cint, $types, $(args...))))
    if length(args) >= 1
        error_expr = :(_la_error(status, $(esc(args[1]))))
    else
        error_expr = :(_la_error(status))
    end
    quote
        status = $call_ex
        if status != Status.OK
            $error_expr
        end
    end
end


text:
using Phylo
using Compat: undef
using .RCall
using .RCall: protect, unprotect, rcall_p, RClass, isObject, isS4
import .RCall: rcopy

function rcopy(::Type{T}, rt::Ptr{VecSxp}) where T <: AbstractTree
    if !isObject(rt) || isS4(rt) || rcopy(rcall_p(:class, rt)) != "phylo"
        error("Object is not of S3 phylo class, aborting")
    end

    if !rcopy(rcall_p(Symbol("is.rooted"), rt))
        error("Cannot currently translate unrooted trees")
    end

    dict = rcopy(Dict{Symbol, Any}, rt)
    nodes = dict[:tip_label]
    tree = NamedTree(nodes)
    edges = dict[:edge]
    nnode = dict[:Nnode]
    lengths = dict[:edge_length]
    nontips = nnode
    append!(nodes, addnodes!(tree, nontips))

    for edge in 1:size(edges, 1)
        addbranch!(tree,
                   nodes[edges[edge, 1]], nodes[edges[edge, 2]],
                   lengths[edge])
    end

    validate(tree) || warn("Tree does not internally validate")
    return tree
end

import .RCall.rcopytype

rcopytype(::Type{RClass{:phylo}}, s::Ptr{VecSxp}) = NamedTree

import .RCall.sexp

function sexp(tree::T) where T <: AbstractTree
    validate(tree) || warn("Tree does not internally validate")

    tipnames = collect(nodenamefilter(isleaf, tree))
    root = collect(nodenamefilter(isroot, tree))
    if (length(root) != 1)
        error("Can't currently translate tree with > 1 roots")
    end
    nontips = collect(nodenamefilter(isinternal, tree))
    tor = Dict{Symbol, Any}()
    tor[:Nnode] = length(nontips) + length(root)
    tor[Symbol("tip.label")] = tipnames
    nodes = copy(tipnames)
    push!(nodes, root[1])
    append!(nodes, nontips)
    bi = branchiter(tree)
    lengths = Vector{Float64}(undef, length(bi))
    edges = Matrix{Int32}(undef, length(lengths), 2)
    index = 1
    for branch in bi
        lengths[index] = getlength(branch)
        edges[index, :] = indexin([src(branch), dst(branch)], nodes)
        index += 1
    end
    tor[:edge] = edges
    tor[Symbol("edge.length")] = lengths
    sobj = protect(sexp(tor))
    #setattrib!(sobj, :order, sexp("cladewise"))
    setclass!(sobj, sexp("phylo"))
    unprotect(1)
    return sobj
end


text:
# This file includes testset for DDESystem

import DifferentialEquations: MethodOfSteps, Vern9

@testset "DDESystemTestSet" begin
    @info "Running DDESystemTestSet ..."

    # DDESystem construction
    out = zeros(1)
    tau = 1
    constlags = [tau]
    histfunc(out, u, t) = (out .= 1.)
    function statefunc(dx, x, h, u, t)
        h(out, u, t - tau) # Update `out`.
        dx[1] = out[1] + x[1]
    end
    outputfunc(x, u, t) = x
    ds = DDESystem(righthandside=statefunc, history=histfunc, readout=outputfunc, state=[1.], input=nothing, output=Outport(1), alg=MethodOfSteps(Vern9()), constlags=constlags, depslags=nothing)
    # ds = DDESystem((statefunc, histfunc), outputfunc, [1.], 0., nothing, Outport(1), alg=MethodOfSteps(Tsit5()), modelkwargs=(constant_lags=constlags,))
    @test typeof(ds.trigger) == Inpin{Float64}
    @test typeof(ds.handshake) == Outpin{Bool}
    @test ds.input === nothing
    @test isa(ds.output, Outport)
    @test length(ds.output) == 1
    @test ds.integrator.sol.prob.constant_lags == constlags
    @test ds.integrator.sol.prob.dependent_lags === nothing
    @test ds.integrator.sol.prob.neutral == false
    @test ds.integrator.alg == Vern9()

    # Driving DDESystem
    iport = Inport(1)
    trg  = Outpin()
    hnd = Inpin{Bool}()
    connect!(ds.output, iport)
    connect!(trg, ds.trigger)
    connect!(ds.handshake, hnd)
    tsk = launch(ds)
    tsk2 = @async while true
        all(take!(iport) .=== NaN) && break
    end
    for t in 1 : 10
        put!(trg, t)
        take!(hnd)
        @test ds.t == t
        @test [read(pin.link.buffer) for pin in iport] == ds.state
    end
    put!(trg, NaN)
    sleep(0.1)
    @test istaskdone(tsk)
    put!(ds.output, NaN * ones(length(ds.output)))
    sleep(0.1)
    @test istaskdone(tsk2)

    # DDESystem with input
    # hist2 = History(histfunc, constlags, ())
    function statefunc2(dx, x, h, u, t)
        h(out, u, t - tau) # Update `out`.
        dx[1] = out[1] + x[1] + sin(u[1](t)) + cos(u[2](t))
    end
    outputfunc2(x, u, t) = x
    ds = DDESystem(righthandside=statefunc2, history=histfunc, readout=outputfunc2, state=[1.], input=Inport(2),
        output=Outport(1), constlags=constlags, depslags=nothing)
    @test isa(ds.input, Inport)
    @test isa(ds.output, Outport)
    @test length(ds.input) == 2
    @test length(ds.output) == 1
    @test typeof(ds.integrator.sol.prob.p) <: Interpolant
    @test size(ds.integrator.sol.prob.p.timebuf) == (3,)
    @test size(ds.integrator.sol.prob.p.databuf) == (2, 3)
    oport = Outport(2)
    iport = Inport(1)
    trg = Outpin()
    hnd = Inpin{Bool}()
    connect!(oport, ds.input)
    connect!(ds.output, iport)
    connect!(trg, ds.trigger)
    connect!(ds.handshake, hnd)
    tsk = launch(ds)
    tsk2 = @async while true
        all(take!(iport) .=== NaN) && break
    end
    for t in 1 : 10
        put!(trg, t)
        put!(oport, [t, 2t])
        take!(hnd)
        @test ds.t == t
        @test [read(pin.link.buffer) for pin in iport] == ds.state
    end
    put!(trg, NaN)
    sleep(0.1)
    @test istaskdone(tsk)
    put!(ds.output, NaN * ones(length(ds.output)))
    sleep(0.1)
    @test istaskdone(tsk2)

    # Test defining new DDESystem
    # The type must be mutable
    @test_throws Exception @eval @def_dde_system struct DDESystem{CL, DL, RH, HST, RO, ST, IP, OP}
        constlags::CL
        depslags::DL
        righthandside::RH
        history::HST
        readout::RO
        state::ST
        input::IP
        output::OP
    end

    # The type must be a subtype of AbstractDAESystem
    @test_throws Exception @eval @def_dde_system mutable struct DDESystem{CL, DL, RH, HST, RO, ST, IP, OP}
        constlags::CL
        depslags::DL
        righthandside::RH
        history::HST
        readout::RO
        state::ST
        input::IP
        output::OP
    end

    # The type must be a subtype of AbstractDAESystem
    @test_throws Exception @eval @def_dde_system mutable struct DDESystem{CL, DL, RH, HST, RO, ST, IP, OP} <: MyDummyAbstractDDESystem
        constlags::CL
        depslags::DL
        righthandside::RH
        history::HST
        readout::RO
        state::ST
        input::IP
        output::OP
    end

    @info "Done DDESystemTestSet ..."
end  # testset


text:
using CSVFiles
using DataFrames
using Distributions
using FileIO
using MacroTools
using StatsBase
using IterTools

include("mcs_types.jl")
include("EmpiricalDistribution.jl")
include("montecarlo.jl")
include("lhs.jl")
include("sobol.jl")
include("defmcs.jl")

export
    @defsim, generate_trials!, run, save_trial_inputs, _save_trial_results, set_models!,
    EmpiricalDistribution, ReshapedDistribution, RandomVariable, TransformSpec, CorrelationSpec,
    SimulationDef, SimulationInstance, AbstractSimulationData,
    LHSData, LatinHypercubeSimulationDef, MCSData, MonteCarloSimulationDef, SobolData, SobolSimulationDef,
    INNER, OUTER, sample!, analyze, MonteCarloSimulationInstance, LatinHypercubeSimulationInstance,
    SobolSimulationInstance


text:
with_repo_access("test lookup raises error if object type does not match") do test_repo, path
    # blob
    @test_throws LibGitError{:Invalid,:NotFound} lookup_tree(test_repo, Oid("fa49b077972391ad58037050f2a75f74e3671e92"))
    # commit
    @test_throws LibGitError{:Invalid,:NotFound} lookup_tree(test_repo, Oid("8496071c1b46c854b31185ea97743be6a8774479"))
    # tag
    @test_throws LibGitError{:Invalid,:NotFound} lookup_tree(test_repo, Oid("0c37a5391bbff43c37f0d0371823a5509eed5b1d"))
end

with_repo_access() do test_repo, path
    id = Oid("c4dc1555e4d4fa0e0c9c3fc46734c7c35b3ce90b")
    test_tree = lookup_tree(test_repo, id)

    context("test read tree data") do
        @test id == Oid(test_tree)
        @test isa(test_tree, GitTree)
        @test length(test_tree) == 3
        @test Oid(test_tree[1]) == Oid("1385f264afb75a56a5bec74243be9b367ba4ca08")
        @test Oid(test_tree[2]) == Oid("fa49b077972391ad58037050f2a75f74e3671e92")
    end

    context("test read tree entry data") do
        bent = test_tree[1]
        tent = test_tree[3]

        @test name(bent) == "README"
        @test isa(bent, GitTreeEntry{GitBlob})

        @test name(tent) == "subdir"
        @test isa(tent, GitTreeEntry{GitTree})
        @test Oid(tent) == Oid("619f9935957e010c419cb9d15621916ddfcc0b96")
        @test isa(lookup(test_repo, Oid(tent)), GitTree)
    end

    context("test get entry by id") do
        bentid = Oid("1385f264afb75a56a5bec74243be9b367ba4ca08")
        bent = test_tree[bentid]
        @test name(bent) == "README"
        @test isa(bent, GitTreeEntry{GitBlob})
    end

    context("test get entry by id returns nothing if no id") do
        nadaid = Oid("1385f264afb75a56a5bec74243be9b367ba4ca07")
        nada = test_tree[nadaid]
        @test nada == nothing
    end

    context("test tree iteration") do
        @test length(collect(test_tree)) == length(test_tree)
        tes = sort(collect(test_tree))
        str = join([name(te) for te in tes], ":")
        @test str == "README:subdir:new.txt"
    end

    context("test tree walk only trees") do
        walk(GitTree, test_tree, :postorder) do res
            root, entry = res
            @assert isa(entry, GitTreeEntry{GitTree})
        end
    end

    context("test tree walk only blob") do
        walk(GitBlob, test_tree, :postorder) do res
            root, entry = res
            @test isa(entry, GitTreeEntry{GitBlob})
        end
    end

    context("test iterate subtrees") do
        for t in foreach(GitTree, test_tree)
            @test isa(t, GitTreeEntry{GitTree})
        end
    end

    context("test iterate subtree blobs") do
        for b in foreach(GitBlob, test_tree)
            @test isa(b, GitTreeEntry{GitBlob})
        end
    end
end

#TODO: more tree builder tests
with_tmp_repo_access("test tree builder") do test_repo, path
  builder = GitTreeBuilder(test_repo)
  insert!(builder, "README.txt",
                   Oid("1385f264afb75a56a5bec74243be9b367ba4ca08"),
                   33188)
  id = write!(builder)
  obj = test_repo[id]
  @test isa(obj, GitTree)
  @test length(raw(obj)) == 38
end


text:
using Interpolations

x  = [NaN, NaN, 2.0, 1.0, 2.0, NaN, 2.5, 2.3, NaN, 5.6, 8.7, NaN]
g = 1:length(x)

intp = interpolate_nans(g, x, Linear(), NaN)

@test isnan(intp[1])
@test isnan(intp[end])
@test any(isnan.(intp[3:end-1])) == false

intp = interpolate_nans(g, x, Linear(), Flat(OnGrid()))

@test intp[1] == intp[3]
@test intp[end] == intp[end-1]
@test any(isnan.(intp[3:end-1])) == false

x1 = interpolate_nans(1:20, x, Linear(), NaN)
x2 = fill_nans(x, Linear())

all(x1[.!(isnan.(x1))] .≈ x2[.!(isnan.(x2))])


text:
using TuringGLM
using Documenter
using Pluto: Configuration.CompilerOptions
using PlutoStaticHTML

tutorials = [
    "Linear Regression",
    "Logistic Regression",
    "Poisson Regression",
    "Negative Binomial Regression",
    "Robust Regression",
    "Hierarchical Models",
    "Custom Priors",
]

include("build.jl")

build()
md_files = markdown_files()
T = [t => f for (t, f) in zip(tutorials, md_files)]

DocMeta.setdocmeta!(TuringGLM, :DocTestSetup, :(using TuringGLM); recursive=true)

makedocs(;
    modules=[TuringGLM],
    authors="Jose Storopoli <thestoropoli@gmail.com>, Rik Huijzer <t.h.huijzer@rug.nl>, and contributors",
    repo="https://github.com/TuringLang/TuringGLM.jl/blob/{commit}{path}#{line}",
    sitename="TuringGLM.jl",
    format=Documenter.HTML(;
        assets=["assets/favicon.ico"],
        canonical="https://TuringLang.github.io/TuringGLM.jl",
        # Using MathJax3 since Pluto uses that engine too.
        mathengine=Documenter.MathJax3(),
        prettyurls=get(ENV, "CI", "false") == "true",
    ),
    pages=["Home" => "index.md", "Tutorials" => T, "API reference" => "api.md"],
)

deploydocs(; repo="github.com/TuringLang/TuringGLM.jl", devbranch="main")

# Useful for local development.
cd(pkgdir(TuringGLM))


text:
const AutocompleteOptions = Vector{Tuple{String, Vector{String}}}

to_autocomplete_options(options::Union{AbstractArray, Tuple}) = vecmap(Tuple, options)
to_autocomplete_options(options::AbstractDict) = [Tuple(p) for p in pairs(options)]
to_autocomplete_options(options::AbstractObservable) = lift(to_autocomplete_options, options)

struct Autocomplete
    value::Observable{String}
    options::Observable{AutocompleteOptions}
    list::List
end

function Autocomplete(value::Observable, options::Observable{AutocompleteOptions})
    return Autocomplete(value, options, List(Observable(String[]), value))
end

function Autocomplete(value::Observable, options′)
    options::Observable{AutocompleteOptions} = to_autocomplete_options(options′)
    return Autocomplete(value, options)
end

function jsrender(session::Session, wdg::Autocomplete)

    hidden, keydown = wdg.list.hidden,wdg.list.keydown

    list = jsrender(session, wdg.list)

    input = DOM.input(
        onfocus=js"JSServe.update_obs($(hidden), false)",
        onblur=js"""
            const tgt = event.relatedTarget;
            if (tgt && $(list).contains(tgt)) {
                this.focus();
            } else {
                JSServe.update_obs($(hidden), true);
            }
        """,
        class="w-full",
        type="text",
        autocomplete="off",
        spellcheck=false,
        value=wdg.value,
        oninput=js"JSServe.update_obs($(wdg.value), this.value)",
        onkeydown=js"event.key == 'Escape' ? this.blur() : JSServe.update_obs($(keydown), event.key)"
    )

    ui = DOM.div(input, list)

    onjs(session, wdg.value, js"""
        function (value) {
            const options = JSServe.get_observable($(wdg.options));
            const pre = options.map(e => e[0]);
            const post = options.map(e => e[1]);
            const lastSpace = value.lastIndexOf(' ');
            const lastColon = value.lastIndexOf(':');
            const idx = Math.max(lastSpace, lastColon);
            let list
                if (lastSpace < lastColon) {
                const key = value.slice(lastSpace + 1, lastColon);
                const extendedKey = value.slice(value.lastIndexOf(' ', lastSpace - 1) + 1, lastColon)
                const target = post[pre.indexOf(extendedKey)] || post[pre.indexOf(key)]
                list = target.map(str => str + ' ');
            } else {
                list = pre.map((str, i) => str + (post[i].length ? ':' : ' '));
            }
            const slice = value.slice(idx + 1, value.length)
            const keys = list.filter(text => text && text.toLowerCase().startsWith(slice.toLowerCase()));
            const entries = keys.map(function (key) {
                const completion = value.slice(0, idx + 1) + key;
                return {key, value: completion};
            });
            JSServe.update_obs($(wdg.list.entries), entries);
        }
    """)

    notify!(wdg.value)

    return jsrender(session, ui)
end

struct RichTextField
    name::String
    widget::Autocomplete
    default::String
    parsed::Vector{Call}
    confirmedvalue::Observable{String}
end

function RichTextField(name, widget::Autocomplete, default)
    rtf = RichTextField(name, widget, default, Call[], Observable(widget.value[]))
    parse!(rtf)
    return rtf
end

function RichTextField(name, options, default)
    wdg = Autocomplete(Observable(default), options)
    return RichTextField(name, wdg, default)
end

function parse!(rtf::RichTextField)
    value = rtf.widget.value[]
    empty!(rtf.parsed)
    append!(rtf.parsed, compute_calls(value))
    rtf.confirmedvalue[] = value
    return
end

function reset!(rtf::RichTextField)
    rtf.widget.value[] = rtf.default
    parse!(rtf)
end

maybethrow(e::Exception; strict) = strict ? throw(e) : nothing

function extract_call(rtf::RichTextField; strict=true)
    calls = rtf.parsed
    length(calls) == 1 || maybethrow(ArgumentError("Field $(rtf.name) must have a unique addend"); strict)
    return first(calls)
end

function extract_positional_argument(rtf::RichTextField; strict=true)
    positionals = extract_call(rtf; strict).positional
    length(positionals) == 1 || maybethrow(ArgumentError("Field $(rtf.name) must have a unique positional argument"); strict)
    return first(positionals)
end

extract_positional_arguments(rtf::RichTextField; strict=true) = extract_call(rtf; strict).positional

function extract_positional_arguments(rtf::RichTextField, n::Int; strict=true)
    positionals = extract_positional_arguments(rtf; strict)
    length(positionals) == n || maybethrow(ArgumentError("Field $(rtf.name) must have exactly $n positional arguments"); strict)
    return positionals
end

extract_named_arguments(rtf::RichTextField; strict=true) = extract_call(rtf; strict).named

function extract_all_arguments(rtf::RichTextField; strict=true)
    positional = extract_positional_arguments(rtf; strict)
    named = extract_named_arguments(rtf; strict)
    return positional, named
end

function extract_function(rtf::RichTextField; strict=true)
    fs = extract_call(rtf; strict).fs
    length(fs) == 1 || maybethrow(ArgumentError("Field $(rtf.name) must have a unique function"); strict)
    return first(fs)
end

function jsrender(session::Session, rtf::RichTextField)
    label = DOM.p(class="text-blue-800 text-xl font-semibold py-4 w-full text-left", rtf.name)
    ui = DOM.div(class="mb-4", label, rtf.widget)
    return jsrender(session, ui)
end


text:
"""
    Chebyshev(n)

Basis of Chebyshev polynomials up to degree `n`.
"""
struct Chebyshev <: AbstractBasis
    n::Int
end

Base.length(B::Chebyshev) = B.n

approxpoints(B::Chebyshev) = ChebyshevPoints(B.n)

fftwtype(::Type{T}) where {T <: FFTW.fftwNumber} = T
fftwtype(::Type{T}) where {T <: Real} = Float64
fftwtype(::Type{T}) where {T <: Complex} = ComplexF64
function approxtransform(B::Chebyshev, f)
    n = length(B)
    T = fftwtype(eltype(f))
    n == 1 && return convert(Array{T},f)
    c = FFTW.r2r(f,FFTW.REDFT00,1)
    d = inv(real(T)(n-1)).*(i->isodd(i) ? 1 : -1).(1:n)
    d[1] /= 2; d[end] /= 2
    return d .* c
end

function iterate_basis(B::Chebyshev, x)
    length(B) == 0 && return nothing
    T0 = one(x)
    return T0,(2,T0,T0)
end
function iterate_basis(B::Chebyshev, x, (i,T0,T1))
    i > length(B) && return nothing
    if i == 2
        T0,T1 = T1, val(x)
    else
        T0,T1 = T1, 2 * (x*T1) - T0
    end
    return T1,(i+1,T0,T1)
end

# Use FFT to evaluate in Chebyshev points
function evaltransform(B::Chebyshev, x::ChebyshevPoints, c)
    # Warning: this code is not tested for mixed precision / BigFloat computations
    @assert length(B) == size(c,1)
    T = fftwtype(promote_type(eltype(x), eltype(c)))

    nB = length(B)-1
    nx = length(x)-1

    if nB == 0
        return one.(x) .* c
    elseif nx == 0
        return Matrix(B,x)*c
    else
        k = (nB+nx-1)÷nx
        pnB = k*nx
        O = zeros(eltype(c), (pnB - nB, Base.tail(size(c))...))
        d = inv(real(T)(2)).*(i->isodd(i) ? 1 : -1).(1:pnB+1); d[1] = 1; d[end] *= 2
        idx = ( 1:k:k*(nx+1), ntuple(i->1:size(c,i+1), ndims(c)-1)... )
        return FFTW.r2r(d.*[c;O], FFTW.REDFT00,1)[idx...]
    end
end


text:
import PlotlyBase
import JSON

function DashBase.to_dash(p::PlotlyBase.Plot)
    data = JSON.lower(p)
    pop!(data, :config, nothing)
    return data
end


text:
using HardSphereDynamics, StaticArrays

table = HardSphereDynamics.RectangularBox(SA[-0.5, -0.5, -1.0],
                SA[+0.5, +0.5, +3.0])

fluid = HardSphereFluid{3,Float64}(table, N, r)

initial_condition!(fluid, lower=table.lower, upper=-table.lower)
collision_type = ElasticCollision()

flow_type = ExternalFieldFlow(SA[0.0, 0.0, -10.0])
event_handler = AllToAll(fluid, flow_type)

simulation =  HardSphereSimulation(
    fluid, event_handler, flow_type, collision_type);

states, times = evolve!(simulation, 0.01, 100);

using Makie

visualize_3d(states, sleep_step=0.005, lower=table.lower, upper=-table.lower)

max_stars_repo_name: dpsanders/HardSphereDynamics

max_forks_repo_name: dpsanders/HardSphereDynamics


text:
module Samplers

export GibbsSampler, SequentialMonteCarlo, SequentialImportanceSampler, MetropolisHastingsSampler
export MetropolisHastingsSampler, MetropolisWithinGibbsSampler

abstract type Sampler end

abstract type MetropolisWithinGibbsSampler <: Sampler  end

struct GibbsSampler <: MetropolisWithinGibbsSampler
    num_iterations::Int64
    num_burn_in::Int64
    skip::Int64
    assignment_types::Vector{String}
    function GibbsSampler(;num_iterations::Int64=1000, num_burn_in::Int64=0, skip::Int64=1,
                          assignment_types::Vector{String}=["random"])
        if num_iterations < 1
            error("Number of iterations should be positive.")
        end
        if num_burn_in < 0
            error("Number of burn in iterations should be positive.")
        end
        if !(all(map(x -> x in ["random", "all same cluster", "all different clusters"], assignment_types)))
            error("Initial assignment types can only be one of 'random', 'all same cluster' or 'all different clusters'")
        end
        return new(num_iterations, num_burn_in, skip, assignment_types)
    end
end

mutable struct MetropolisHastingsSampler <: MetropolisWithinGibbsSampler
    num_iterations::Int64
    num_burn_in::Int64
    proposal_radius::Int64
    skip::Int64
    adaptive::Bool
    assignment_types::Vector{String}
    function MetropolisHastingsSampler(;num_iterations=1000, num_burn_in=1000, proposal_radius=5, skip=1,
                                       adaptive=false, assignment_types::Vector{String}=["random"])
        if num_iterations < 1
            error("Number of iterations should be positive.")
        end
        if num_burn_in < 0
            error("Number of burn in iterations should be positive.")
        end
        if !(all(map(x -> x in ["random", "all same cluster", "all different clusters"], assignment_types)))
            error("Initial assignment types can only be one of 'random', 'all same cluster' or 'all different clusters'")
        end
        return new(num_iterations, num_burn_in, proposal_radius, skip, adaptive, assignment_types)
    end
end

struct SequentialMonteCarlo <: Sampler
    num_particles::Int64
    ess_threshold::Float64
    function SequentialMonteCarlo(;num_particles::Int64=1000, ess_threshold::T=0.5) where
                                  {T <: Real}
        if num_particles < 1
            error("Number of particles should be positive.")
        end
        if ess_threshold < 0
            error("ESS Threshold should be non-negative.")
        end
        if typeof(ess_threshold) <: AbstractFloat && ess_threshold > 1
            error("ESS Threshold should be either a non-negative integer or a fraction")
        end
        if ess_threshold < 1
            ess_threshold = num_particles*ess_threshold
        end
        return new(num_particles, ess_threshold)
    end
end

struct SequentialImportanceSampler <: Sampler
    num_particles::Int64
    function SequentialImportanceSampler(;num_particles=1000)
        return new(num_particles)
    end
end

end

max_stars_repo_name: realseanla/ntl-stick-breaking-julia

max_issues_repo_name: realseanla/ntl-stick-breaking-julia

max_forks_repo_name: realseanla/ntl-stick-breaking-julia


text:
@sde_model BlackScholes dS = r*S*dt + σ*S*dW

function _bs_transition_distr(m, s, x0)
  @unpack r, σ = m
  μ = log(x0) + (r-0.5*σ^2)*s.Δt
  d = sqrt(s.Δt)*σ
  LogNormal(μ, d)
end

function sample(model::BlackScholes, scheme::Exact, t0, x0)
  d = _bs_transition_distr(model, scheme, x0)
  rand(d)
end

function transition(model::BlackScholes, scheme::Exact, t0, x0, x1)
  d = _bs_transition_distr(model, scheme, x0)
  pdf(d, x1)
end

function logtransition(model::BlackScholes, scheme::Exact, t0, x0, x1)
  d = _bs_transition_distr(model, scheme, x0)
  logpdf(d, x1)
end


text:
tt = 0.:dt:T

Random.seed!(2)
# generate one long path
T_long = 10.0
tt_long = 0.:dt:T_long
W_long = sample(tt_long, Wiener())
X_long = solve(Euler(), x0, W_long, P)
# write long forward to csv file
f = open(outdir*"longforward_fh.csv","w")
head = "time, component, value \n"
write(f, head)
longpath = [Any[tt_long[j], d, X_long.yy[j][d]] for d in 1:2, j in 1:10:length(X_long) ][:]
writedlm(f, longpath, ',')
close(f)

# simulate forwards, on the shorter interval
Random.seed!(3)
W = sample(tt, Wiener())
X = solve(Euler(), x0, W, P)
XX = [X]
samples = 100
# draw more paths
for j in 2:samples
    W = sample(tt, Wiener())
    X = solve(Euler(), x0, W, P)
    push!(XX, X)
end

# write forwards to csv file
f = open(outdir*"forwards_fh.csv","w")
head = "samplenr, time, component, value \n"
write(f, head)
iterates = [Any[i, tt[j], d, XX[i].yy[j][d]] for d in 1:2, j in 1:10:length(X), i in 1:samples ][:]
writedlm(f, iterates, ',')
close(f)

# simulate bridges by forward simulating and waiting, in case of a extreme conditioning
#v = 0.6   # conditioning in extreme case
v = 1.1

Random.seed!(3)
W = sample(tt, Wiener())
X = solve(Euler(), x0, W, P)
XX = []
samples = 30
s = 0
# draw more paths
while true
    sample!(W, Wiener())
    solve!(Euler(), X, x0, W, P)
    if norm(L*X.yy[end] .- v) < 0.01
        push!(XX, copy(X))
        s += 1
        s % 10 == 0 && println(".")
        s >= samples && break
    end
end

# write 'forward-bridges' to csv file

f = open(outdir*"bridges_extreme_blindforward.csv","w")
head = "samplenr, time, component, value \n"
write(f, head)
iterates = [Any[i, tt[j], d, XX[i].yy[j][d]] for d in 1:2, j in 1:10:length(X), i in 1:samples ][:]
writedlm(f, iterates, ',')
close(f)

# simulate bridges by forward simulating and waiting, in case of a extreme conditioning
#v = 0.6   # conditioning in extreme case
v = -1

Random.seed!(3)
W = sample(tt, Wiener())
X = solve(Euler(), x0, W, P)
XX = []
samples = 30
s = 0
# draw more paths
while true
    sample!(W, Wiener())
    solve!(Euler(), X, x0, W, P)
    if norm(L*X.yy[end] .- v) < 0.01
        push!(XX, copy(X))
        s += 1
        s % 10 == 0 && println(".")
        s >= samples && break
    end
end

# write 'forward-bridges' to csv file
f = open(outdir*"bridges_first_blindforward.csv","w")
head = "samplenr, time, component, value \n"
write(f, head)
iterates = [Any[i, tt[j], d, XX[i].yy[j][d]] for d in 1:2, j in 1:10:length(X), i in 1:samples ][:]
writedlm(f, iterates, ',')
close(f)

tt = τ(T).(0.:dt:T)


text:
using Pkg
# First, activate the project: this should load all dependencies, etc
Pkg.activate(".")
# `instantiate` actually ensures that the dependency are exactly at the same state than on my computer
Pkg.instantiate()

# Load the module
using VinbergsAlgorithmNF

# Load Hecke, which deals with everything number theoretic
using Hecke

# And LinearAlgebra
using LinearAlgebra

# And CoxeterDiagrams, which contains the algo to check co{compactness,finiteness}
using CoxeterDiagrams

# **Important**: if you want the code to run fast enough, always call `toggle(false)` first.
# I peppered the code with asserts at some performance critical places, and calling `toggle(false)` disables those asserts.
using ToggleableAsserts
toggle(false)

# This is just a shortcut to not have to write `VinbergsAlgorithmNF` in full everytime
VA = VinbergsAlgorithmNF

# This is "orthogonal sum"
⊕ = VA.Lat.:(⊕)

# These are the usual lattices (defined in `src/some_lattices.jl`)
H = VA.Lat.U()
block_diag = VA.Lat.block_diag

# Define the field
ℚ, a = Hecke.rationals_as_number_field()

# Quadratic forms from Turkalj in dimension 5 (up to rational isometry)

Turkalj_dim_5 = [
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1],
    [0, -1, 0, 0, 1, 0, -1, -2, 0, 0, 1, -1, 0, 0, 2, -1, 1, -1, 0, 0, -1, 2, -1, 1, 1, 1, 1, -1, 2, 0, 0, -1, -1, 1, 0, 2],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0, 1, 2],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 3],
    [0, -1, 0, 0, -1, -2, -1, -2, -1, 1, 0, 0, 0, -1, 2, 0, 0, 0, 0, 1, 0, 2, 1, 1, -1, 0, 0, 1, 10, 2, -2, 0, 0, 1, 2, 12],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -2, 1, 0, 0, 0, 0, 1, 2],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -2],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 6],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -6],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 7],
    [0, -1, 0, 0, 0, -2, -1, 0, 0, 0, 1, -2, 0, 0, 2, -1, 1, 1, 0, 0, -1, 4, 1, 1, 0, 1, 1, 1, 4, 1, -2, -2, 1, 1, 1, -4],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 15],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -2, 1, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 3],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -2, 1, 0, 0, 0, 0, 1, 7],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 10, 5, 0, 0, 0, 0, 5, 10],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, -6],
    [0, -1, -1, -3, 2, -6, -1, 0, 0, 1, -1, -3, -1, 0, 2, -1, 0, 0, -3, 1, -1, 12, -3, -4, 2, -1, 0, -3, 12, -4, -6, -3, 0, -4, -4, -20],
    [0, 0, 0, 0, 0, 1, 0, 2, 0, 0, -1, -1, 0, 0, 2, 1, 0, 1, 0, 0, 1, 2, 0, 0, 0, -1, 0, 0, 4, 1, 1, -1, 1, 0, 1, -6],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 14],
    [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -10],
    [2, 1, -1, 1, 0, 0, 1, 2, -1, 0, 0, 0, -1, -1, 2, -1, 0, 0, 1, 0, -1, 2, 0, 0, 0, 0, 0, 0, 2, -1, 0, 0, 0, 0, -1, -10],
    [2, 0, 0, 0, 0, 1, 0, 2, 1, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 2, 0, 1, 0, 0, 0, 0, 2, -1, 1, 0, 0, 1, -1, -6]
]

function is_good(das)
    m, n = size(das.D, 1), size(das.D, 2)
    @assert m==n "!!! Something wrong with the diagram matrix !!! \n"
    for i in 1:n
        face_neighbours = das.D[i,:]
        if all(x->(x<=2), face_neighbours)
            return true
        end
    end
    return false
end

function data_and_diagrams_dim_5()
    ind = 0
    for L in Turkalj_dim_5
        ind = ind + 1
        print("Starting Vinberg algorithm for no. ", ind, "\n")
        vd = VA.VinbergData(ℚ, reshape(L, 6, 6))
        (status, (roots, dict, das)) = VA.next_n_roots!(vd, n=128)
        print("Number of roots: ", length(roots), "\n")
        print("All roots found: \n")
        print(roots, "\n")
        print("\n")
        @assert status==true "!!! Error at no. "*string(ind)*" !!!\n"
        if is_good(das)
            print("Has an isolated root (good) \n")
            CoxeterDiagrams.save_png("/Users/sasha/VinbergsAlgorithmNF/scratchpad/Turkalj-dim-5/Good/tur5no"*string(ind)*".png", das)
        else
            print("Does *not* have an isolated root (bad) \n")
            CoxeterDiagrams.save_png("/Users/sasha/VinbergsAlgorithmNF/scratchpad/Turkalj-dim-5/Bad/tur5no"*string(ind)*".png", das)
        end
    end
end

# Drawing the diagrams and saving them each in a separate PNG file

data_and_diagrams_dim_5()

# b as in "bad"

print("\n")
print("Commensurable lattices for bad lattices ... \n")
print("\n")

for b in [1, 2, 3, 5, 11]
    L = Turkalj_dim_5[b]
    print("Starting Vinberg algorithm for no. ", b, " commensurable \n")
    c = [1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 2]
    q = c*reshape(L, 6, 6)*transpose(c)
    vd = VA.VinbergData(ℚ, q)
    (status, (roots, dict, das)) = VA.next_n_roots!(vd, n=128)
    print("Number of roots: ", length(roots), "\n")
    print("All roots found: \n")
    print(roots, "\n")
    print("\n")
    @assert status==true "!!! Error at no. "*string(b)*" !!! \n"
    @assert is_good(das) "!!! Is not good (yet) !!! \n"
    CoxeterDiagrams.save_png("/Users/sasha/VinbergsAlgorithmNF/scratchpad/Turkalj-dim-5/Bad/tur5no"*string(b)*"comm.png", das)
end

max_stars_repo_name: sashakolpakov/crystallographic-packings

max_issues_repo_name: sashakolpakov/crystallographic-packings

max_forks_repo_name: sashakolpakov/crystallographic-packings


text:
"""

The exported names are
$(EXPORTS)
"""
module Integration

using DocStringExtensions
using Test
using QuadGK: gauss
using FastGaussQuadrature: gaussjacobi
using FastGaussQuadrature: gausslegendre
using FillArrays

using Gridap.Helpers
using Gridap.Fields
using Gridap.TensorValues
using Gridap.Arrays

export Quadrature
export GenericQuadrature
export TensorProductQuadrature
export DuffyQuadrature
export num_points
export get_coordinates
export get_weights
export num_dims
export num_point_dims
export test_quadrature

include("Quadratures.jl")

include("TensorProductQuadratures.jl")

include("DuffyQuadratures.jl")

end # module


text:
if !haskey(ENV, "BUILDKITE")
    import Pkg
    Pkg.develop(Pkg.PackageSpec(; path = dirname(@__DIR__)))
end
using Test

using KernelAbstractions
const KA = KernelAbstractions

import CUDAKernels
const CK = CUDAKernels

import Thermodynamics
const TD = Thermodynamics

import UnPack

using Random
using LinearAlgebra

import RootSolvers
const RS = RootSolvers

import CLIMAParameters
const CP = CLIMAParameters

struct EarthParameterSet <: CP.AbstractEarthParameterSet end
const param_set = EarthParameterSet()

if get(ARGS, 1, "Array") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
    CUDA.allowscalar(false)
    device(::Type{T}) where {T <: CUDA.CuArray} = CK.CUDADevice()
else
    ArrayType = Array
    device(::Type{T}) where {T <: Array} = CK.CPU()
end

@show ArrayType

@kernel function test_thermo_kernel!(
    param_set,
    dst::AbstractArray{FT},
    e_int,
    ρ,
    p,
    q_tot,
) where {FT}
    i = @index(Group, Linear)
    @inbounds begin

        ts = TD.PhaseEquil_ρeq(param_set, FT(ρ[i]), FT(e_int[i]), FT(q_tot[i]))
        dst[1, i] = TD.air_temperature(param_set, ts)

        ts_ρpq = TD.PhaseEquil_ρpq(
            param_set,
            FT(ρ[i]),
            FT(p[i]),
            FT(q_tot[i]),
            true,
            100,
            RS.RegulaFalsiMethod,
        )
        dst[2, i] = TD.air_temperature(param_set, ts_ρpq)
    end
end

# Since we use `rand` to generate the ProfileSet,
# just initialize on the CPU, and provide convert
# function to move arrays to the GPU.
convert_profile_set(ps::TD.TestedProfiles.ProfileSet, ArrayType, slice) =
    TD.TestedProfiles.ProfileSet(
        ArrayType(ps.z[slice]),
        ArrayType(ps.T[slice]),
        ArrayType(ps.p[slice]),
        ArrayType(ps.RS[slice]),
        ArrayType(ps.e_int[slice]),
        ArrayType(ps.h[slice]),
        ArrayType(ps.ρ[slice]),
        ArrayType(ps.θ_liq_ice[slice]),
        ArrayType(ps.q_tot[slice]),
        ArrayType(ps.q_liq[slice]),
        ArrayType(ps.q_ice[slice]),
        TD.PhasePartition.(ps.q_tot[slice], ps.q_liq[slice], ps.q_ice[slice]),
        ArrayType(ps.RH[slice]),
        ArrayType(ps.e_pot[slice]),
        ArrayType(ps.u[slice]),
        ArrayType(ps.v[slice]),
        ArrayType(ps.w[slice]),
        ArrayType(ps.e_kin[slice]),
        ps.phase_type,
    )

@testset "Thermodynamics - kernels" begin
    FT = Float32
    dev = device(ArrayType)
    profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, Array)
    slice = Colon()
    profiles = convert_profile_set(profiles, ArrayType, slice)

    n_profiles = length(profiles.z)
    n_vars = length(propertynames(profiles))
    d_dst = ArrayType(Array{FT}(undef, 2, n_profiles))
    fill!(d_dst, 0)

    UnPack.@unpack e_int, ρ, p, q_tot = profiles

    work_groups = (1,)
    ndrange = (n_profiles,)
    kernel! = test_thermo_kernel!(dev, work_groups)
    event = kernel!(param_set, d_dst, e_int, ρ, p, q_tot, ndrange = ndrange)
    wait(dev, event)

    ts_correct =
        TD.PhaseEquil_ρeq.(param_set, Array(ρ), Array(e_int), Array(q_tot))
    @test all(Array(d_dst)[1, :] .≈ TD.air_temperature.(param_set, ts_correct))

    ts_correct =
        TD.PhaseEquil_ρpq.(
            param_set,
            Array(ρ),
            Array(p),
            Array(q_tot),
            true,
            100,
            RS.RegulaFalsiMethod,
        )
    @test all(Array(d_dst)[2, :] .≈ TD.air_temperature.(param_set, ts_correct))

end


text:
@inline function abs_fp_dd(x::T) where {T<:AbstractFloat}
    return abs(x), zero(T)
end

@inline function negabs_fp_dd(x::T) where {T<:AbstractFloat}
    return -abs(x), zero(T)
end

@inline function neg_fp_dd(x::T) where {T<:AbstractFloat}
    return -x, zero(T)
end

@inline function inv_fp_dd(x::T) where {T<:AbstractFloat}
    hi, lo = two_inv(x)
    return hi, lo
end

@inline function square_fp_dd(a::T) where {T<:AbstractFloat}
    zhi = a * a
    zlo = fma(a, a, -zhi)
    return zhi, zlo
end

@inline function cube_fp_dd(a::T) where {T<:AbstractFloat}
    zhi, zlo = square_fp_dd(a)
    zhi, zlo = mul_dddd_dd((zhi, zlo), (a, zero(T)))
    return zhi, zlo
end

@inline function powr2_fp_dd(x::T) where {T<:AbstractFloat}
    return powr2_(x)
end

@inline function powr3_fp_dd(x::T) where {T<:AbstractFloat}
    return powr3_(x)
end

@inline function powr4_fp_dd(x::T) where {T<:AbstractFloat}
    hi, lo3 = powr3_fp_dd(x)
    hi, lo = mul_(x, hi)
    lo = fma(lo3, x, lo)
    return hi, lo
end

@inline function powr5_fp_dd(x::T) where {T<:AbstractFloat}
    hi, lo4 = powr4_fp_dd(x)
    hi, lo = mul_(x, hi)
    lo = fma(lo4, x, lo)
    return hi, lo
end

@inline function powr6_fp_dd(x::T) where {T<:AbstractFloat}
    hi, lo5 = powr5_fp_dd(x)
    hi, lo = mul_(x, hi)
    lo = fma(lo5, x, lo)
    return hi, lo
end

@inline function root2_fp_dd(x::T) where {T<:AbstractFloat}
    return two_sqrt(x)
end

@inline function root3_fp_dd(x::T) where {T<:AbstractFloat}
    return two_cbrt(x)
end

@inline function root4_fp_dd(x::T) where {T<:AbstractFloat}
     return two_sqrt(sqrt(x))
end

@inline function root5_fp_dd(x::T) where {T<:AbstractFloat}
    throw(ErrorException("not implemented"))
end

@inline function root6_fp_dd(x::T) where {T<:AbstractFloat}
    return two_cbrt(sqrt(x))
end


text:
using Reproducible, Test

import Reproducible: CodeBlock

module TestMod
end

@testset "CodeBlock" begin
    cb = CodeBlock("x = 1\ny = 2", TestMod)
    @test strip(cb[1].input) == "x = 1"
    @test strip(cb[2].input) == "y = 2"
    @test cb[1].output == 1
    @test cb[2].output == 2
end


text:
export Gates

module Gates

export GateOperator, QuantumCircuit
export x, y, z, h, Id, cx, cy, rk, swap

"Struct to represent operators acting on specific qubits"
mutable struct GateOperator
    n::Integer
    qubits::Array{Integer, 1}
    gate::Array{T, 2} where T <: Number
    GateOperator(qubits::Array{<:Integer, 1}, gate::Array{<:Number, 2}) =
        new(length(qubits), qubits, gate)
end

"Struct to represent a quantum circuit"
mutable struct QuantumCircuit
    n::Integer
    ops::Array{GateOperator, 1}
    QuantumCircuit(N) = new(N, Array{GateOperator, 1}())
end

"""
Define some simple gates
"""
x = [[0., 1.] [1., 0.]]
z = [[1., 0.] [0., -1.]]
y = .- z * x .* im
s = [[1., 0.] [0., 1im]]
t = [[1., 0.] [0., exp(π*1im/4)]]
h = 1/sqrt(2)*[[1., 1.] [1., -1.]]
Id = [[1., 0.] [0., 1.]]
cx = [[1., 0., 0., 0.] [0., 1., 0., 0.] [0., 0., 0., 1.] [0., 0., 1., 0.]]
cy = [[1., 0., 0., 0.] [0., 1., 0., 0.] [0., 0., 0., -1im] [0., 0., 1im, 0.]]
rk = k -> [[1., 0., 0., 0.] [0., 1., 0., 0.] [0., 0., 1., 0.] [0., 0., 0., exp(2π*1im/2^k)]]
swap = [[1. ,0., 0., 0.] [0., 0., 1., 0.] [0., 1., 0., 0] [0., 0., 0., 1.]]

# IBM physical gates
u3(θ, ϕ, λ) = [[cos(θ/2), exp(1im*λ)*sin(θ/2)] [-exp(1im*ϕ)*sin(θ/2), exp(1im*(λ + ϕ))*cos(θ/2)]]
u2(ϕ, λ) = u3(0, ϕ, λ)
u1(λ) = u3(0, 0, λ)
end # module