text:
module Azure

using Swagger
using HTTP
using JSON

# API versions
const _module_versions = Dict{Module,String}()
const _api_versions = Dict{DataType,String}()

# include Azure REST protocol (not Swagger)
include("rest_api_protocol.jl")

# inlcude sub modules for individual API groups
include("DataLakeStore/DataLakeStoreAccountManagementClient/DataLakeStoreAccountManagementClient.jl")
_module_versions[DataLakeStoreAccountManagementClient] = "2016-11-01"
_api_versions[DataLakeStoreAccountManagementClient.AccountApi] = "2016-11-01"
_api_versions[DataLakeStoreAccountManagementClient.FirewallRulesApi] = "2016-11-01"
_api_versions[DataLakeStoreAccountManagementClient.LocationsApi] = "2016-11-01"
_api_versions[DataLakeStoreAccountManagementClient.OperationsApi] = "2016-11-01"
_api_versions[DataLakeStoreAccountManagementClient.TrustedIdProvidersApi] = "2016-11-01"

include("DataLakeStore/DataLakeStoreFileSystemManagementClient/DataLakeStoreFileSystemManagementClient.jl")
_module_versions[DataLakeStoreFileSystemManagementClient] = "2016-11-01"
_api_versions[DataLakeStoreFileSystemManagementClient.FileSystemApi] = "2016-11-01"

include("Storage/StorageManagementClient/StorageManagementClient.jl")
_module_versions[StorageManagementClient] = "2017-06-01"
_api_versions[StorageManagementClient.OperationsApi] = "2017-06-01"
_api_versions[StorageManagementClient.SkusApi] = "2017-06-01"
_api_versions[StorageManagementClient.StorageAccountsApi] = "2017-06-01"
_api_versions[StorageManagementClient.UsageApi] = "2017-06-01"

include("Compute/ComputeManagementClient/ComputeManagementClient.jl")
_module_versions[ComputeManagementClient] = "2017-03-30"
_api_versions[ComputeManagementClient.AvailabilitySetsApi] = "2017-03-30"
_api_versions[ComputeManagementClient.ImagesApi] = "2017-03-30"
_api_versions[ComputeManagementClient.SkusApi] = "2017-03-30"
_api_versions[ComputeManagementClient.UsageApi] = "2017-03-30"
_api_versions[ComputeManagementClient.VirtualMachineExtensionImagesApi] = "2017-03-30"
_api_versions[ComputeManagementClient.VirtualMachineExtensionsApi] = "2017-03-30"
_api_versions[ComputeManagementClient.VirtualMachineImagesApi] = "2017-03-30"
_api_versions[ComputeManagementClient.VirtualMachineScaleSetExtensionsApi] = "2017-03-30"
_api_versions[ComputeManagementClient.VirtualMachineScaleSetRollingUpgradesApi] = "2017-03-30"
_api_versions[ComputeManagementClient.VirtualMachineScaleSetVMsApi] = "2017-03-30"
_api_versions[ComputeManagementClient.VirtualMachineScaleSetsApi] = "2017-03-30"
_api_versions[ComputeManagementClient.VirtualMachineSizesApi] = "2017-03-30"
_api_versions[ComputeManagementClient.VirtualMachinesApi] = "2017-03-30"

include("Compute/DiskResourceProviderClient/DiskResourceProviderClient.jl")
_module_versions[DiskResourceProviderClient] = "2017-03-30"
_api_versions[DiskResourceProviderClient.DisksApi] = "2017-03-30"
_api_versions[DiskResourceProviderClient.SnapshotsApi] = "2017-03-30"

include("Compute/RunCommandsClient/RunCommandsClient.jl")
_module_versions[RunCommandsClient] = "2017-03-30"
_api_versions[RunCommandsClient.VirtualMachineRunCommandsApi] = "2017-03-30"
_api_versions[RunCommandsClient.VirtualMachinesApi] = "2017-03-30"

include("Compute/ContainerServiceClient/ContainerServiceClient.jl")
_module_versions[ContainerServiceClient] = "2017-01-31"
_api_versions[ContainerServiceClient.ContainerServiceApi] = "2017-01-31"
_api_versions[ContainerServiceClient.ContainerServicesApi] = "2017-01-31"

include("Network/ApplicationGatewayClient/ApplicationGatewayClient.jl")
_module_versions[ApplicationGatewayClient] = "2017-03-01"
_api_versions[ApplicationGatewayClient.ApplicationGatewaysApi] = "2017-03-01"

include("Network/CheckDnsAvailabilityClient/CheckDnsAvailabilityClient.jl")
_module_versions[CheckDnsAvailabilityClient] = "2017-03-01"
_api_versions[CheckDnsAvailabilityClient.DefaultApi] = "2017-03-01"

include("Network/ExpressRouteCircuitClient/ExpressRouteCircuitClient.jl")
_module_versions[ExpressRouteCircuitClient] = "2017-03-01"
_api_versions[ExpressRouteCircuitClient.ExpressRouteCircuitArpTableApi] = "2017-03-01"
_api_versions[ExpressRouteCircuitClient.ExpressRouteCircuitAuthorizationsApi] = "2017-03-01"
_api_versions[ExpressRouteCircuitClient.ExpressRouteCircuitPeeringsApi] = "2017-03-01"
_api_versions[ExpressRouteCircuitClient.ExpressRouteCircuitRoutesTableApi] = "2017-03-01"
_api_versions[ExpressRouteCircuitClient.ExpressRouteCircuitRoutesTableSummaryApi] = "2017-03-01"
_api_versions[ExpressRouteCircuitClient.ExpressRouteCircuitStatsApi] = "2017-03-01"
_api_versions[ExpressRouteCircuitClient.ExpressRouteCircuitsApi] = "2017-03-01"
_api_versions[ExpressRouteCircuitClient.ExpressRouteServiceProvidersApi] = "2017-03-01"

include("Network/LoadBalancerClient/LoadBalancerClient.jl")
_module_versions[LoadBalancerClient] = "2017-03-01"
_api_versions[LoadBalancerClient.LoadBalancersApi] = "2017-03-01"

include("Network/NetworkManagementClient/NetworkManagementClient.jl")
_module_versions[NetworkManagementClient] = "2017-03-01"

include("Network/NetworkInterfaceClient/NetworkInterfaceClient.jl")
_module_versions[NetworkInterfaceClient] = "2017-03-01"
_api_versions[NetworkInterfaceClient.NetworkInterfacesApi] = "2017-03-01"

include("Network/NetworkSecurityGroupClient/NetworkSecurityGroupClient.jl")
_module_versions[NetworkSecurityGroupClient] = "2017-03-01"
_api_versions[NetworkSecurityGroupClient.NetworkSecurityGroupsApi] = "2017-03-01"
_api_versions[NetworkSecurityGroupClient.SecurityRulesApi] = "2017-03-01"

include("Network/NetworkWatcherClient/NetworkWatcherClient.jl")
_module_versions[NetworkWatcherClient] = "2017-03-01"
_api_versions[NetworkWatcherClient.NetworkWatchersApi] = "2017-03-01"
_api_versions[NetworkWatcherClient.PacketCapturesApi] = "2017-03-01"

include("Network/PublicIpAddressClient/PublicIpAddressClient.jl")
_module_versions[PublicIpAddressClient] = "2017-03-01"
_api_versions[PublicIpAddressClient.PublicIPAddressesApi] = "2017-03-01"

include("Network/RouteFilterClient/RouteFilterClient.jl")
_module_versions[RouteFilterClient] = "2017-03-01"
_api_versions[RouteFilterClient.RouteFilterRulesApi] = "2017-03-01"
_api_versions[RouteFilterClient.RouteFiltersApi] = "2017-03-01"

include("Network/RouteTableClient/RouteTableClient.jl")
_module_versions[RouteTableClient] = "2017-03-01"
_api_versions[RouteTableClient.RouteTablesApi] = "2017-03-01"
_api_versions[RouteTableClient.RoutesApi] = "2017-03-01"

include("Network/ServiceCommunityClient/ServiceCommunityClient.jl")
_module_versions[ServiceCommunityClient] = "2017-03-01"
_api_versions[ServiceCommunityClient.BgpServiceCommunitiesApi] = "2017-03-01"

include("Network/UsageClient/UsageClient.jl")
_module_versions[UsageClient] = "2017-03-01"
_api_versions[UsageClient.UsagesApi] = "2017-03-01"

include("Network/VirtualNetworkClient/VirtualNetworkClient.jl")
_module_versions[VirtualNetworkClient] = "2017-03-01"
_api_versions[VirtualNetworkClient.DefaultApi] = "2017-03-01"
_api_versions[VirtualNetworkClient.SubnetsApi] = "2017-03-01"
_api_versions[VirtualNetworkClient.VirtualNetworkPeeringsApi] = "2017-03-01"
_api_versions[VirtualNetworkClient.VirtualNetworksApi] = "2017-03-01"

include("Network/VirtualNetworkGatewayClient/VirtualNetworkGatewayClient.jl")
_module_versions[VirtualNetworkGatewayClient] = "2017-03-01"
_api_versions[VirtualNetworkGatewayClient.LocalNetworkGatewaysApi] = "2017-03-01"
_api_versions[VirtualNetworkGatewayClient.VirtualNetworkGatewayConnectionsApi] = "2017-03-01"
_api_versions[VirtualNetworkGatewayClient.VirtualNetworkGatewaysApi] = "2017-03-01"

include("Network/VmssNetworkInterfaceClient/VmssNetworkInterfaceClient.jl")
_module_versions[VmssNetworkInterfaceClient] = "2017-03-01"
_api_versions[VmssNetworkInterfaceClient.NetworkInterfacesApi] = "2017-03-01"

include("Resource/ResourceManagementClient/ResourceManagementClient.jl")
_module_versions[ResourceManagementClient] = "2017-05-10"
_api_versions[ResourceManagementClient.DeploymentOperationsApi] = "2017-05-10"
_api_versions[ResourceManagementClient.DeploymentsApi] = "2017-05-10"
_api_versions[ResourceManagementClient.ProvidersApi] = "2017-05-10"
_api_versions[ResourceManagementClient.ResourceGroupsApi] = "2017-05-10"
_api_versions[ResourceManagementClient.ResourcesApi] = "2017-05-10"
_api_versions[ResourceManagementClient.TagsApi] = "2017-05-10"

include("Resource/SubscriptionClient/SubscriptionClient.jl")
_module_versions[SubscriptionClient] = "2016-06-01"
_api_versions[SubscriptionClient.SubscriptionsApi] = "2016-06-01"
_api_versions[SubscriptionClient.TenantsApi] = "2016-06-01"

include("Resource/PolicyClient/PolicyClient.jl")
_module_versions[PolicyClient] = "2016-04-01"
_api_versions[PolicyClient.PolicyAssignmentsApi] = "2016-04-01"
_api_versions[PolicyClient.PolicyDefinitionsApi] = "2016-04-01"

include("Commerce/UsageManagementClient/UsageManagementClient.jl")
_module_versions[UsageManagementClient] = "2015-06-01-preview"
_api_versions[UsageManagementClient.RateCardApi] = "2015-06-01-preview"
_api_versions[UsageManagementClient.UsageAggregatesApi] = "2015-06-01-preview"

# Storage services
include("Storage/StorageServices/StorageServices.jl")

# helper methods to assist in authentication, logging and such
include("helper.jl")

end # module Azure


text:
@testset "Test mesh reader" begin

using HDiscontinuousGalerkin
using Tensors

# Test read external mesh in Triangle format
root_file = "mesh/figure2.1"
mesh = parse_mesh_triangle(root_file)

# General checks
@test getncells(mesh) == 4
@test getnnodes(mesh) == 5
@test n_faces_per_cell(mesh) == 3
for cell in mesh.cells
    @test volume(get_coordinates(cell,  mesh)) ≈ 1/4
end
@test getcells_matrix(mesh) == [2 3 5; 4 1 5; 5 3 4; 1 2 5]
@test get_vertices_matrix(mesh) == [0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0; 0.5 0.5]
# Check expected data for cell 1
@test mesh.cells[1].nodes == (2,3,5)
@test mesh.cells[1].faces == (1,2,3)
@test [face_orientation(mesh,1,i) for i in 1:3] == [true,false,true]
@test cell_diameter(mesh,1) == 1.0
@test getfaceset(mesh, "boundary") == Set([3,6,7,8])
@test get_coordinates(1, mesh) == [Vec{2}([1.0, 1.0]), Vec{2}([0.5, 0.5])]
@test getnfaces(mesh.cells[1]) == 3

#Test generated mesh
mesh = rectangle_mesh(TriangleCell, (2,2), Vec{2}((0.0,0.0)), Vec{2}((1.0,1.0)))
@test getncells(mesh) == 8
@test getnnodes(mesh) == 9
@test n_faces_per_cell(mesh) == 3
for cell in mesh.cells
    @test volume(get_coordinates(cell,  mesh)) ≈ 1/8
end
@test getcells_matrix(mesh) == [1 2 4; 2 5 4; 2 3 5; 3 6 5; 4 5 7; 5 8 7; 5 6 8; 6 9 8]
@test get_vertices_matrix(mesh) == [0.0 0.0; 0.5 0.0; 1.0 0.0; 0.0 0.5; 0.5 0.5; 1.0 0.5; 0.0 1.0; 0.5 1.0; 1.0 1.0]
# Check expected data for cell 1
@test mesh.cells[1].nodes == (1,2,4)
@test mesh.cells[1].faces == (1,2,3)
@test [face_orientation(mesh,1,i) for i in 1:3] == [true,false,true]
@test cell_diameter(mesh,1) == sqrt(2)/2
@test getfaceset(mesh, "boundary") == Set([3,7,9,16,2,11,12,15])
@test getnfaces(mesh.cells[1]) == 3

end

ext: jl


text:
module Utilities

export log2i, log2dim, group_shift, lmove, bsubspace, bcomspace,
    bmask, anyone, literal_pow, ismatch, Hoistable, VarDefs,
    hoist!, hoisted_expr, def!, block, @def, @defnew, eval_interp, var_exprs,
    flatten_block

using OrderedCollections
using YaoLocations: plain, AbstractLocations

function block(f, it)
    ret = Expr(:block)
    for xs in it
        push!(ret.args, f(xs...))
    end
    return ret
end

function flatten_block(ex)
    ex isa Expr || return ex
    ex.head === :block && return _flatten_block(ex)
    return Expr(ex.head, map(flatten_block, ex.args)...)
end

function _flatten_block(ex)
    ex isa Expr || return ex
    ex.head === :block || return ex
    ret = Expr(:block)
    for stmt in ex.args
        if stmt isa Expr && stmt.head === :block
            append!(ret.args, flatten_block(stmt).args)
        else
            push!(ret.args, flatten_block(stmt))
        end
    end
    return ret
end

function eval_interp(m::Module, ex)
    ex isa Expr || return ex
    ex.head === :$ && return Base.eval(m, ex.args[1])
    return Expr(ex.head, map(x->eval_interp(m, x), ex.args)...)
end

struct Hoistable
    origin::OrderedDict{Any, Any}
    hoisted::OrderedDict{Any, Any}
end

function Hoistable(X::Matrix)
    origin = OrderedDict()
    for i in axes(X, 1), j in axes(X, 2)
        origin[(i, j)] = X[i, j]
    end
    return Hoistable(origin, OrderedDict())
end

function hoist!(h::Hoistable, key, others...)
    key = isempty(others) ? key : (key, others...)
    haskey(h.hoisted, key) && return h.hoisted[key]
    @gensym hoisted
    h.hoisted[key] = hoisted
    return hoisted
end

function Base.getindex(h::Hoistable, key, others...)
    key = isempty(others) ? key : (key, others...)
    haskey(h.hoisted, key) && return h.hoisted[key]
    return h.origin[key]
end

function hoisted_expr(h::Hoistable)
    ret = Expr(:block)
    for (key, sym) in h.hoisted
        push!(ret.args, :($sym = $(h.origin[key])))
    end
    return ret
end

struct VarDefs
    mangled::OrderedDict{Symbol, Symbol}
    symbols::OrderedDict{Symbol, Any}
end

VarDefs() = VarDefs(OrderedDict{Symbol, Symbol}(), OrderedDict{Symbol, Any}())

function def!(defs::VarDefs, sym, expr; force=false)
    mangled = getfield(defs, :mangled)
    symbols = getfield(defs, :symbols)
    force || haskey(symbols, sym) && return defs
    symbols[sym] = expr
    mangled[sym] = gensym(sym)
    return defs
end

macro def(defs, ex::Expr)
    ex.head === :(=) || error("expect an assignment expr")
    quote
        $def!($defs, $(QuoteNode(ex.args[1])), $(Expr(:quote, ex.args[2])))
        $(ex.args[1]) = $defs.$(ex.args[1])
    end |> esc
end

macro defnew(defs, ex::Expr)
    ex.head === :(=) || error("expect an assignment expr")
    quote
        $def!($defs, $(QuoteNode(ex.args[1])), $(Expr(:quote, ex.args[2])); force=true)
        $(ex.args[1]) = $defs.$(ex.args[1])
    end |> esc
end

function var_exprs(defs::VarDefs)
    mangled = getfield(defs, :mangled)
    symbols = getfield(defs, :symbols)
    isempty(symbols) && return
    ret = Expr(:block)
    for (sym, expr) in symbols
        push!(ret.args, Expr(:(=), mangled[sym], expr))
    end
    return ret
end

function Base.getproperty(defs::VarDefs, name::Symbol)
    mangled = getfield(defs, :mangled)
    haskey(mangled, name) && return mangled[name]
    error("$name is not defined")
end

function Base.show(io::IO, defs::VarDefs)
    mangled = getfield(defs, :mangled)
    symbols = getfield(defs, :symbols)
    println(io, "global variables:")
    for (k, (sym, expr)) in enumerate(symbols)
        print(io, "  ", mangled[sym], ": ", sym, " = ", expr)
        if k != length(symbols)
            println(io)
        end
    end
end

@inline function literal_pow(::Val{:im}, p::Int)
    m = mod(p, 4)
    return if m == 0
        Complex{Int}(1, 0)
    elseif m == 1
        Complex{Int}(0, 1)
    elseif m == 2
        Complex{Int}(-1, 0)
    else
        Complex{Int}(0, -1)
    end
end
# fast pow of ``exp(im * pi / d)^x``
@generated function literal_pow(::Val{:exp}, ::Val{d}, x::Int) where d
    root = Expr(:if)
    curr = root
    for k in 0:(2d-2)
        push!(curr.args, :(m == $k))
        push!(curr.args, exp(im * π * k / d))
        next = Expr(:elseif)
        push!(curr.args, next)
        curr = next
    end

    push!(curr.args, :(m == $(2d-2)))
    push!(curr.args, exp(im * π * (2d-2) / d))
    push!(curr.args, exp(im * π * (2d-1) / d))

    return quote
        m = mod(x, 2d)
        $root
    end
end

"""
    log2i(x::Integer) -> Integer

Return log2(x), this integer version of `log2` is fast but only valid for number equal to 2^n.
"""
function log2i end

for N in [8, 16, 32, 64, 128]
    T = Symbol(:Int, N)
    UT = Symbol(:UInt, N)
    @eval begin
        log2i(x::$T) =
            !signbit(x) ? ($(N - 1) - leading_zeros(x)) :
            throw(ErrorException("nonnegative expected ($x)"))
        log2i(x::$UT) = $(N - 1) - leading_zeros(x)
    end
end

log2dim(x::AbstractVector) = log2i(size(x, 1))
log2dim(x::AbstractMatrix) = log2i(size(x, 2))

@inline lmove(b::Int, mask::Int, k::Int)::Int = (b & ~mask) << k + (b & mask)

"""
    bmask(::Type{T}) where T <: Integer -> zero(T)
    bmask([T::Type], positions::Int...) -> T
    bmask([T::Type], range::UnitRange{Int}) -> T
Return an integer mask of type `T` where `1` is the position masked according to
`positions` or `range`. Directly use `T` will return an empty mask `0`.
"""
function bmask end

bmask(args...) = bmask(Int, args...)
bmask(::Type{T}) where {T<:Integer} = zero(T)
bmask(::Type{T}, positions::Int...) where {T<:Integer} = bmask(T, positions)

function bmask(::Type{T}, itr) where {T<:Integer}
    isempty(itr) && return 0
    ret = zero(T)
    for b in itr
        ret += one(T) << (b - 1)
    end
    return ret
end

@inline function bmask(::Type{T}, range::UnitRange{Int})::T where {T<:Integer}
    ((one(T) << (range.stop - range.start + 1)) - one(T)) << (range.start - 1)
end

"""
    anyone(index::Integer, mask::Integer) -> Bool
Return `true` if any masked position of index is 1.
# Example
`true` if any masked positions is 1.
```jldoctest
julia> anyone(0b1011, 0b1001)
true
julia> anyone(0b1011, 0b1100)
true
julia> anyone(0b1011, 0b0100)
false
```
"""
anyone(index::T, mask::T) where {T<:Integer} = (index & mask) != zero(T)

"""
    ismatch(index::Integer, mask::Integer, target::Integer) -> Bool

Return `true` if bits at positions masked by `mask` equal to `1` are equal to `target`.

## Example

```julia
julia> n = 0b11001; mask = 0b10100; target = 0b10000;

julia> ismatch(n, mask, target)
true
```
"""
ismatch(index::T, mask::T, target::T) where {T<:Integer} = (index & mask) == target

@inline function _group_shift(masks::Vector{Int}, shift_len::Vector{Int}, k::Int, k_prv::Int)
    # if current position in the contiguous region
    # since these bits will be moved together with
    # the first one, we don't need to generate a
    # new mask
    if k == k_prv + 1
        shift_len[end] += 1
    else
        # we generate a bit mask where the 1st to k-th bits are 1
        push!(masks, bmask(0:k-1))
        push!(shift_len, 1)
    end
end

group_shift(locations::AbstractLocations) = group_shift(plain(locations))
complement_group_shift(n::Int, locations::AbstractLocations) = complement_group_shift(n, plain(locations))

@inline function group_shift(locations)
    masks = Int[]
    shift_len = Int[]
    k_prv = -1
    for k in locations
        _group_shift(masks, shift_len, k, k_prv)
        k_prv = k
    end
    return masks, shift_len
end

@inline function complement_group_shift(n::Int, locations)
    masks = Int[]
    shift_len = Int[]
    k_prv = -1
    for k in 1:n
        k in locations && continue
        _group_shift(masks, shift_len, k, k_prv)
        k_prv = k
    end
    return masks, shift_len
end

struct BitSubspace
    n::Int # number of bits in fullspace
    sz_subspace::Int # size of the subspace
    n_shifts::Int # number of shifts
    masks::Vector{Int} # shift masks
    shift_len::Vector{Int} # length of each shift
end

function Base.getindex(s::BitSubspace, i::Int)
    index = i - 1
    @inbounds for k in 1:s.n_shifts
        index = lmove(index, s.masks[k], s.shift_len[k])
    end
    return index
end

Base.firstindex(s::BitSubspace) = 1
Base.lastindex(s::BitSubspace) = s.sz_subspace
Base.length(s::BitSubspace) = s.sz_subspace
Base.eltype(::BitSubspace) = Int

function Base.iterate(s::BitSubspace, st::Int = 1)
    st <= length(s) || return
    return s[st], st + 1
end

function bsubspace(n::Int, locs)
    @assert issorted(locs)
    masks, shift_len = group_shift(locs)
    BitSubspace(n, 1 << (n - length(locs)), length(masks), masks, shift_len)
end

function bcomspace(n::Int, locs)
    @assert issorted(locs)
    masks, shift_len = complement_group_shift(n, locs)
    BitSubspace(n, 1 << length(locs), length(masks), masks, shift_len)
end

function Base.show(io::IO, ::MIME"text/plain", s::BitSubspace)
    indent = get(io, :indent, 0)
    println(io, " "^indent, s.sz_subspace, "-element BitSubspace:")
    if s.sz_subspace < 5
        for k in 1:s.sz_subspace
            print(io, " "^(indent+1), string(s[k]; base=2, pad=s.n))
            if k != s.sz_subspace
                println(io)
            end
        end
    else # never print more than 4 elements
        println(io, " "^(indent+1), string(s[1]; base=2, pad=s.n))
        println(io, " "^(indent+1), string(s[2]; base=2, pad=s.n))
        println(io, " "^(indent+1), "⋮")
        println(io, " "^(indent+1), string(s[end-1]; base=2, pad=s.n))
        print(io, " "^(indent+1), string(s[end]; base=2, pad=s.n))
    end
end

end # Utilities


text:
abstract type PaddleStatelessModule end
abstract type PaddleStatelessLayer end

# a rough solution for General Net
struct PaddleStatelessGeneralNet<:PaddleStatelessModule
    NN::PyObject
end

function (stateless_module::PaddleStatelessGeneralNet)(params::Vector, inputs; kwinputs...)
    map((p,p_new)->p.set_value(p_new), stateless_module.NN.parameters(), params)
    out = stateless_module.NN(inputs)
    return out
end

struct PaddleStatelessFCNet<:PaddleStatelessModule
    layers::Vector{PaddleStatelessLayer}
end

function (stateless_module::PaddleStatelessFCNet)(params::Vector, inputs; kwinputs...)
    out = PyNULL()
    copy!(out, inputs)
    state = 1
    for layer in stateless_module.layers
        state = PaddleLayerForward!(out, params, state, layer)
    end
    return out
end

function PaddleStatelessModule(paddle_module)
    if pybuiltin("isinstance")(paddle_module, paddle.nn.Sequential)
        layers = PaddleStatelessLayer[]
        for layer in paddle_module
            if pybuiltin("isinstance")(layer, paddle.nn.Linear)
                push!(layers, PaddleLinear(layer.weight.shape...))
            elseif pybuiltin("isinstance")(layer, (paddle.nn.Sigmoid, paddle.nn.ReLU, paddle.nn.Tanh))
                push!(layers, PaddleActivation(layer))
            else
                throw(error("Unsupported layer type"))
            end
        end
        paddel_fc_net = PaddleStatelessFCNet(layers)
        return paddel_fc_net
    else
        paddle_net = PaddleStatelessGeneralNet(paddle_module)
        return paddle_net
    end
end

# wrap paddle's layer
struct PaddleLinear<:PaddleStatelessLayer
    features_ins::Int
    features_outs::Int
end

Base.show(io::IO, f::PaddleLinear) = print(io, "Linear($(f.features_ins), $(f.features_outs))")

function PaddleLayerForward!(out::PyObject, params::Vector, state::Int, L::PaddleLinear)
    weight, state = iterate(params, state)
    bias, state = iterate(params, state)
    copy!(out, paddle.matmul(out, weight))
    copy!(out, paddle.add(out, bias))
    return state
end

struct PaddleActivation<:PaddleStatelessLayer
    act::PyObject
end

Base.show(io::IO, f::PaddleActivation) = print(io, "$(f.act)")

function PaddleLayerForward!(out::PyObject, params::Vector, state::Int, L::PaddleActivation)
    copy!(out, L.act(out))
    return state
end

# easy constructor for full connected network
function PaddleFCNet(dim_ins, dim_outs, num_layers, hidden_size; dtype="float32", activation="sigmoid")
    act = PyNULL()
    if activation == "sigmoid"
        copy!(act, paddle.nn.Sigmoid())
    elseif activation == "tanh"
        copy!(act, paddle.nn.Tanh())
    elseif activation == "relu"
        copy!(act, paddle.nn.ReLU())
    else
        throw(error("Unsupported activation type"))
    end

    if !(dtype in ["float16", "float32", "float64"])
        throw(error("Unsupported data type, use \"float16\", \"float32\" or \"float64\""))
    end

    pyparams = PyObject[]
    layers = PaddleStatelessLayer[]
    lsize, rsize = 0, 0
    for i in 1:num_layers
        if i == 1
            lsize = dim_ins
            rsize = hidden_size
        elseif i == num_layers
            lsize = hidden_size
            rsize = dim_outs
        else
            lsize = hidden_size
            rsize = hidden_size
        end
        push!(pyparams, paddle.static.create_parameter(shape=PyVector([lsize, rsize]), dtype=dtype))
        push!(pyparams, paddle.static.create_parameter(shape=PyVector([rsize]), dtype=dtype))
        push!(layers, PaddleLinear(lsize, rsize))
        if i != num_layers
            push!(layers, PaddleActivation(act))
        end
    end
    jlparams = map(x->DLPack.wrap(x, pyto_dlpack), pyparams)
    dtype = eltype(jlparams[1])
    paddel_fc_net = PaddleStatelessFCNet(layers)
    return PaddleModuleWrapper(paddel_fc_net, dtype, jlparams)
end

max_stars_repo_name: PangHua/InfiniteScience

max_issues_repo_name: PangHua/InfiniteScience

max_forks_repo_name: PangHua/InfiniteScience


text:
function read_svm_data(pathname::String)
	mat_x_i = Array(Int, 0)
	mat_x_j = Array(Int, 0)
	mat_x_v = Array(Float64, 0)
	vec_y = Array(Float64, 0)
	io = open(pathname, "r")
	j = 1
	for line in eachline(io)
		fields = split(chomp(line), " ")
		push!(vec_y, parsefloat(Float64, fields[1]))
		for k in 2:(length(fields) - 1)
			feature_index, feature_value = split(fields[k], ":")
			push!(mat_x_i, parseint(Int, feature_index))
			push!(mat_x_j, j)
			push!(mat_x_v, parsefloat(Float64, feature_value))
		end
		j += 1
	end
	close(io)
	return vec_y, sparse(mat_x_i, mat_x_j, mat_x_v)
end


text:
export newtonpolygon, vertices, vertexindices, lattices, facets, subdivision

"""
    newtonpolygon(p::AbstractPolynomial; lowerhull=false)

Construct a [`NewtonPolygon`](@ref) from the support of the polynomial `p`. `lowerhull` indicates
whether the lower or upper convex hull should be used to construct the regular subdivision
of the Newton polygon.
"""
function newtonpolygon(p::MP.AbstractPolynomial; kwargs...)
    @assert MP.nvariables(p) == 2 "Expected a bivariate polynomial."
    lattices = map(t -> SVector{2}(MP.exponents(t)), MP.terms(p))
    coeffs = realcoefficients(p)
    return newtonpolygon(lattices, coeffs; kwargs...)
end
function newtonpolygon(p::MP.AbstractPolynomial{<:Complex})
    @assert MP.nvariables(p) == 2 "Expected a bivariate polynomial."
    lattices = map(t -> SVector{2}(MP.exponents(t)), MP.terms(p))
    return newtonpolygon(lattices)
end

function newtonpolygon(lattice_points::Vector{SVector{2, Int}}, lift::Union{Nothing,AbstractVector{<:Real}}=nothing; lowerhull=false)
    P = Polyhedra.polyhedron(Polyhedra.vrep(lattice_points), CDDLib.Library())
    Polyhedra.removevredundancy!(P)
    if Polyhedra.dim(P) != 2
        throw(AssertionError("Newton polygon is not fulldimensional"))
    end

    vpts = sort_counterclockwise(collect(Polyhedra.points(P)))
    vertices = map(vpts) do v
        for (i,p) in enumerate(lattice_points)
            p == v && return i
        end
    end

    # facets
    idx = Int[]
    for p in vpts
        for vid in vertices
            if p == lattice_points[vid]
                push!(idx, vid)
                break
            end
        end
    end
    # iterate reverse to obtain counter-clockwise order
    facets = [SVector(idx[i], idx[i+1]) for i in 1:length(idx)-1]
    push!(facets, SVector(idx[end], idx[1]))

    # subdivison
    subdivision = Vector{Int}[]
    if lift !== nothing
        lifted_pts = map((p, w) -> [p;w], lattice_points, lift)
        lifted_P = Polyhedra.polyhedron(Polyhedra.vrep(lifted_pts), CDDLib.Library())
        Polyhedra.removevredundancy!(lifted_P)
        if Polyhedra.dim(lifted_P) == 2
            push!(subdivision, copy(vertices))
        else
            hspaces = Polyhedra.halfspaces(lifted_P)
            for (idx, H) in zip(eachindex(hspaces), hspaces)
                if (lowerhull ? H.a[end] < 0 : H.a[end] > 0)
                    ipts = Polyhedra.incidentpoints(lifted_P, idx)
                    pts = sort_counterclockwise(map(p -> SVector(p[1], p[2]), ipts))
                    push!(subdivision, map(pts) do p
                        for (i, v) in enumerate(lattice_points)
                            p[1] == v[1] && p[2] == v[2] && return i
                        end
                    end)
                end
            end
        end
    end

    return NewtonPolygon(lattice_points, vertices, facets, subdivision)
end

function getsemihull(ps::Vector{PT}, sign_sense, counterclockwise, yray = nothing) where PT
    hull = PT[]
    if length(ps) == 0
        return hull
    end
    prev = sign_sense == 1 ? first(ps) : last(ps)
    cur = prev
    for j in (sign_sense == 1 ? (2:length(ps)) : ((length(ps)-1):-1:1))
        while prev != cur && counterclockwise(cur - prev, ps[j] - prev) >= 0
            cur = prev
            pop!(hull)
            if !isempty(hull)
                prev = last(hull)
            end
        end
        if yray !== nothing && counterclockwise(ps[j] - cur, yray) >= 0
            break
        else
            push!(hull, cur)
            prev = cur
            cur = ps[j]
        end
    end
    push!(hull, cur)
    hull
end

function sort_counterclockwise(ps)
    sort!(ps, by = first)
    counterclockwise(p1, p2) = dot(cross([p1; 0], [p2; 0]), [0, 0, 1])
    top = getsemihull(ps,  1, counterclockwise)
    bot = getsemihull(ps, -1, counterclockwise)
    if !isempty(top) && !isempty(bot)
        @assert top[end] == bot[1]
        pop!(top)
    end
    if !isempty(bot) && !isempty(top)
        @assert bot[end] == top[1]
        pop!(bot)
    end
    sort!(top, rev=true)
    sort!(bot)
    [top; bot]
end

function newtonpolygon(exponents::Matrix)
    @assert size(exponents, 1) == 2 "Expected a bivariate polynomial"
    lattices = [SVector(exponents[1, j], exponents[2, j]) for j=1:size(exponents, 2)]
    newtonpolygon(lattices)
end

"""
    lattices(newton_polygon)

The lattice points of the Newton polygon
"""
lattices(np::NewtonPolygon) = np.lattices

"""
    vertexindices(newton_polygon)

The indices of the vertices of the Newton polygon.
"""
vertexindices(np::NewtonPolygon) = np.vertices

"""
    vertices(newton_polygon)

The coordinates of the vertices of the Newton polygon.
"""
vertices(np::NewtonPolygon) = np.lattices[np.vertices]

"""
    facets(newton_polygon)

The facets of the Newton polygon. Each facet is a `SVector{2}` with the coordinates
of the start and end vertex of the facet.
"""
facets(np::NewtonPolygon) = map(f -> SVector(np.lattices[f[1]], np.lattices[f[2]]), np.facets)

iscontained(np::NewtonPolygon, p) = iscontained(vertices(np), p)
on_boundary(np::NewtonPolygon, p) = on_boundary(vertices(np), p)

"""
    subdivision(newton_polygon)

The regular subdivision of the Newton polygon. Each simplex is a vector with the coordinates
of the vertices of the simplex.
"""
function subdivision(np::NewtonPolygon)
    map(np.subdivision) do s
        np.lattices[s]
    end
end

@recipe function plot(newt::NewtonPolygon; showsubdivision=false)
    xs, ys = unzip(vertices(newt))
    xlims --> (minimum(xs) - 1, maximum(xs) + 1)
    ylims --> (minimum(ys) - 1, maximum(ys) + 1)
    border --> nothing
    legend --> nothing
    aspect_ratio --> :equal
    @series begin
        seriestype := :shape
        fill --> :transparent
        linewidth --> 1
        [xs; first(xs)], [ys; first(ys)]
    end

    if showsubdivision
        for simplex in subdivision(newt)
            @series begin
                seriestype := :shape
                fill --> :transparent
                linewidth --> 1
                unzip([simplex..., simplex[1]])
            end
        end
    end

    @series begin
        seriestype := :scatter
        markercolor --> :black
        markerstrokewidth --> 0
        markersize --> 6
        unzip(lattices(newt))
    end
end

"""
    convexhull_vertices(line::Vector{Tuple{Int, Int}})

Returns a list of indices of the vertices of the convex hull.
"""
function convexhull_vertices(points::Vector{NTuple{2, T}}) where {T<:Real}
    convexhull_vertices(SVector.(points))
end
function convexhull_vertices(points::Vector{SVector{2, T}}) where {T<:Real}
    P = Polyhedra.polyhedron(Polyhedra.vrep(points), CDDLib.Library())
    Polyhedra.removevredundancy!(P)
    if Polyhedra.dim(P) != 2
        throw(ArgumentError("points are not full dimensional"))
    end
    vpts = sort_counterclockwise(collect(Polyhedra.points(P)))
    vertices = map(vpts) do v
        for (i,p) in enumerate(points)
            p == v && return i
        end
    end
    vertices
end


text:
module BoundaryDescriptors

using Gridap

export BoundaryDescriptor

struct BoundaryDescriptor
  facet_to_cell::IndexCellNumber{<:Integer}
  facet_to_lfacet::IndexCellNumber{<:Integer}
  cell_to_polytope::IndexCellValue{<:Polytope}
  cell_phi::CellGeomap
  #TODO
  #facet_to_perm::IndexCellNumber{<:Integer}
end

end # module


text:
# This file constains sink tools for the objects of Causal.

"""
    @def_sink ex

where `ex` is the expression to define to define a new AbstractSink component type. The usage is as follows:
```julia
@def_sink struct MySink{T1,T2,T3,...,TN, A} <: AbstractSink
    param1::T1 = param1_default     # optional field
    param2::T2 = param2_default     # optional field
    param3::T3 = param3_default     # optional field
        ⋮
    paramN::TN = paramN_default     # optional field
    action::A = action_function     # mandatory field
end
```
Here, `MySink` has `N` parameters and `action` function

!!! warning
    `action` function must have a method `action(sink::MySink, t, u)` where `t` is the time data and `u` is the data flowing into the sink.

!!! warning
    New static system must be a subtype of `AbstractSink` to function properly.

# Example
```jldoctest
julia> @def_sink struct MySink{A} <: AbstractSink
       action::A = actionfunc
       end

julia> actionfunc(sink::MySink, t, u) = println(t, u)
actionfunc (generic function with 1 method)

julia> sink = MySink();

julia> sink.action(sink, ones(2), ones(2) * 2)
[1.0, 1.0][2.0, 2.0]
```
"""
macro def_sink(ex)
    ex.args[2].head == :(<:) && ex.args[2].args[2] == :AbstractSink ||
        error("Invalid usage. The type should be a subtype of AbstractSink.\n$ex")
    fields = quote
        trigger::$(TRIGGER_TYPE_SYMBOL) = Inpin()
        handshake::$(HANDSHAKE_TYPE_SYMBOL) = Outpin{Bool}()
        callbacks::$(CALLBACKS_TYPE_SYMBOL) = nothing
        name::Symbol = Symbol()
        id::$(ID_TYPE_SYMBOL) = Causal.uuid4()
        input::$(INPUT_TYPE_SYMBOL) = Inport()
        buflen::Int = 64
        plugin::$(PLUGIN_TYPE_SYMBOL) = nothing
        timebuf::$(TIMEBUF_TYPE_SYMBOL) = Buffer(buflen)
        # TIP: Since column-wise reading and writing data is more performant, we construct buffers of the form
        # Buffer(d, l) where d is the dimension and l is the length of buffer.
        databuf::$(DATABUF_TYPE_SYMBOL) = length(input) == 1 ? Buffer(buflen) :  Buffer(length(input), buflen)
        sinkcallback::$(SINK_CALLBACK_TYPE_SYMBOL) = plugin === nothing ?
            Callback(sink->ishit(databuf), sink->action(sink, outbuf(timebuf), outbuf(databuf)), true, id) :
            Callback(sink->ishit(databuf), sink->action(sink, outbuf(timebuf), plugin.process(outbuf(databuf))), true, id)
    end, [TRIGGER_TYPE_SYMBOL, HANDSHAKE_TYPE_SYMBOL, CALLBACKS_TYPE_SYMBOL, ID_TYPE_SYMBOL, INPUT_TYPE_SYMBOL, PLUGIN_TYPE_SYMBOL, TIMEBUF_TYPE_SYMBOL, DATABUF_TYPE_SYMBOL, SINK_CALLBACK_TYPE_SYMBOL]
    _append_common_fields!(ex, fields...)
    deftype(ex)
end

##### Define sink library

# ----------------------------- Writer --------------------------------
"""
    Writer(input=Inport(); buflen=64, plugin=nothing, callbacks=nothing, name=Symbol(uuid4()),
        path=joinpath(tempdir(), string(name)))
Constructs a `Writer` whose input bus is `input`. `buflen` is the length of the internal buffer of `Writer`. If not nothing, `plugin` is used to processes the incomming data. `path` determines the path of the file of `Writer`.

!!! note
    The type of `file` of `Writer` is [`JLD2`](https://github.com/JuliaIO/JLD2.jl).

!!! warning
    When initialized, the `file` of `Writer` is closed. See [`open(writer::Writer)`](@ref) and [`close(writer::Writer)`](@ref).
"""
@def_sink mutable struct Writer{A, FL} <: AbstractSink
    action::A = write!
    path::String = joinpath(tempdir(), string(uuid4()))
    file::FL = (f = jldopen(path, "w"); close(f); f)
end

"""
    write!(writer, td, xd)

Writes `xd` corresponding to `xd` to the file of `writer`.

# Example
```julia
julia> w = Writer(Inport(1))
Writer(path:/tmp/e907d6ad-8db2-4c4a-9959-5b8d33d32156.jld2, nin:1)

julia> open(w)
Writer(path:/tmp/e907d6ad-8db2-4c4a-9959-5b8d33d32156.jld2, nin:1)

julia> write!(w, 0., 10.)
10.0

julia> write!(w, 1., 20.)
20.0

julia> w.file
JLDFile /tmp/e907d6ad-8db2-4c4a-9959-5b8d33d32156.jld2 (read/write)
 ├─🔢 0.0
 └─🔢 1.0

julia> w.file[string(0.)]
10.0
```
"""
write!(writer::Writer, td, xd) = fwrite!(writer.file, td, xd)
fwrite!(file, td, xd) = file[string(td)] = xd

"""
    read(writer::Writer, flatten=false)

Read the contents of the file of `writer` and returns the sorted content of the file. If `flatten` is `true`, the content is also flattened.
"""
read(writer::Writer; flatten=true) = fread(writer.file.path, flatten=flatten)

"""
    fread(path::String)

Reads the content of `jld2` file and returns the sorted file content.
"""
function fread(path::String; flatten=false)
    content = load(path)
    data = SortedDict([(eval(Meta.parse(key)), val) for (key, val) in zip(keys(content), values(content))])
    if flatten
        t = vcat(reverse.(keys(data), dims=1)...)
        if typeof(data) <: SortedDict{T1, T2, T3} where {T1, T2<:AbstractVector, T3}
            x = vcat(reverse.(values(data), dims=1)...)
        elseif typeof(data) <: SortedDict{T1, T2, T3} where {T1, T2<:AbstractMatrix, T3}
            x = collect(hcat(reverse.(values(data), dims=2)...)')
        end
        return t, x
    else
        return data
    end
end

flatten(content) = (collect(vcat(keys(content)...)), collect(vcat(values(content)...)))

"""
    mv(writer::Writer, dst; force::Bool=false)

Moves the file of `writer` to `dst`. If `force` is `true`, the if `dst` is not a valid path, it is forced to be constructed.

# Example
```julia
julia> mkdir(joinpath(tempdir(), "testdir1"))
"/tmp/testdir1"

julia> mkdir(joinpath(tempdir(), "testdir2"))
"/tmp/testdir2"

julia> w = Writer(Inport(), path="/tmp/testdir1/myfile.jld2")
Writer(path:/tmp/testdir1/myfile.jld2, nin:1)

julia> mv(w, "/tmp/testdir2")
Writer(path:/tmp/testdir2/myfile.jld2, nin:1)
```
"""
function mv(writer::Writer, dst; force::Bool=false)
    # id = writer.id
    id = basename(writer.file.path)
    dstpath = joinpath(dst, string(id))
    srcpath = writer.file.path
    mv(srcpath, dstpath, force=force)
    writer.file.path = dstpath  # Update the file path
    writer
end

"""
    cp(writer::Writer, dst; force=false, follow_symlinks=false)

Copies the file of `writer` to `dst`. If `force` is `true`, the if `dst` is not a valid path, it is forced to be constructed. If `follow_symlinks` is `true`, symbolinks are followed.

# Example
```julia
julia> mkdir(joinpath(tempdir(), "testdir1"))
"/tmp/testdir1"

julia> mkdir(joinpath(tempdir(), "testdir2"))
"/tmp/testdir2"

julia> w = Writer(Inport(), path="/tmp/testdir1")
Writer(path:/tmp/testdir1.jld2, nin:1)

julia> cp(w, "/tmp/testdir2")
Writer(path:/tmp/testdir2/1e72bad1-9800-4ca0-bccd-702afe75e555, nin:1)
```
"""
function cp(writer::Writer, dst; force=false, follow_symlinks=false)
    # id = writer.id
    id = basename(writer.file.path)
    dstpath = joinpath(dst, string(id))
    cp(writer.file.path, dstpath, force=force, follow_symlinks=follow_symlinks)
    writer
end

"""
    open(writer::Writer)

Opens `writer` by opening the its `file` in  `read/write` mode. When `writer` is not openned, it is not possible to write data in `writer`. See also [`close(writer::Writer)`](@ref)
"""
open(writer::Writer) = (writer.file = jldopen(writer.file.path, "a"); writer)

"""
    close(writer::Writer)

Closes `writer` by closing its `file`. When `writer` is closed, it is not possible to write data in `writer`. See also [`open(writer::Writer)`](@ref)
"""
close(writer::Writer) =  (close(writer.file); writer)

# ----------------------------- Printer --------------------------------
"""
  Printer(input=Inport(); buflen=64, plugin=nothing, callbacks=nothing, name=Symbol()) where T
Constructs a `Printer` with input bus `input`. `buflen` is the length of its internal `buflen`. `plugin` is data proccessing tool.
"""
@def_sink mutable struct Printer{A} <: AbstractSink
    action::A = print
end

import Base.print

"""
    print(printer::Printer, td, xd)

Prints `xd` corresponding to `xd` to the console.
"""
print(printer::Printer, td, xd) = print("For time", "[", td[1], " ... ", td[end], "]", " => ", xd, "\n")

"""
    open(printer::Printer)

Does nothing. Just a common interface function ot `AbstractSink` interface.
"""
open(printer::Printer) = printer

"""
    close(printer::Printer)

Does nothing. Just a common interface function ot `AbstractSink` interface.
"""
close(printer::Printer) =  printer

# ----------------------------- Scope --------------------------------
"""
    Scope(input=Inport(), args...; buflen::Int=64, plugin=nothing, callbacks=nothing, name=Symbol(), kwargs...)
Constructs a `Scope` with input bus `input`. `buflen` is the length of the internal buffer of `Scope`. `plugin` is the additional data processing tool. `args`,`kwargs` are passed into `plots(args...; kwargs...))`. See (https://github.com/JuliaPlots/Plots.jl) for more information.

!!! warning
    When initialized, the `plot` of `Scope` is closed. See [`open(sink::Scope)`](@ref) and [`close(sink::Scope)`](@ref).
"""
@def_sink mutable struct Scope{A, PA, PK, PLT} <: AbstractSink
    action::A = update!
    pltargs::PA = ()
    pltkwargs::PK = NamedTuple()
    plt::PLT = plot(pltargs...; pltkwargs...)
end

"""
    update!(s::Scope, x, yi)

Updates the series of the plot windows of `s` with `x` and `yi`.
"""
function update!(s::Scope, x, yi)
    y = collect(hcat(yi...)')
    plt = s.plt
    subplots = plt.subplots
    clear.(subplots)
    plot!(plt, x, y, xlim=(x[1], x[end]), label="")  # Plot the new series
    gui()
end

clear(sp::Plots.Subplot) = popfirst!(sp.series_list)  # Delete the old series

"""
    close(sink::Scope)

Closes the plot window of the plot of `sink`.
"""
close(sink::Scope) = closeall()

"""
    open(sink::Scope)

Opens the plot window for the plots of `sink`.
"""
open(sink::Scope) = Plots.isplotnull() ? (@warn "No current plots") : gui()

##### Pretty printing

show(io::IO, writer::Writer) = print(io, "Writer(path:$(writer.file.path), nin:$(length(writer.input)))")
show(io::IO, printer::Printer) = print(io, "Printer(nin:$(length(printer.input)))")
show(io::IO, scp::Scope) = print(io, "Scope(nin:$(length(scp.input)))")

##### Deprecated

# function unfasten!(sink::AbstractSink)
#     callbacks = sink.callbacks
#     sid = sink.id
#     typeof(callbacks) <: AbstractVector && disable!(callbacks[callback.id == sid for callback in callbacks])
#     typeof(callbacks) <: Callback && disable!(callbacks)
#     sink
# end


text:
using Weave

Weave.set_chunk_defaults!(:echo => false, :results => "raw")

include("../src/lotka-volterra-2d-singular.jl")

weave("weave/lotka-volterra-2d-singular-vprk-gauss.jmd",
         out_path = "build/lotka-volterra-2d-singular",
         doctype = "github",
         mod = LotkaVolterra2dSingularSRK)

weave("weave/lotka-volterra-2d-singular-vprk-lobatto-ab.jmd",
         out_path = "build/lotka-volterra-2d-singular",
         doctype = "github",
         mod = LotkaVolterra2dSingularSRK)

weave("weave/lotka-volterra-2d-singular-vprk-lobatto-ba.jmd",
         out_path = "build/lotka-volterra-2d-singular",
         doctype = "github",
         mod = LotkaVolterra2dSingularSRK)

weave("weave/lotka-volterra-2d-singular-vprk-radau.jmd",
         out_path = "build/lotka-volterra-2d-singular",
         doctype = "github",
         mod = LotkaVolterra2dSingularSRK)

weave("weave/lotka-volterra-2d-singular-firk.jmd",
         out_path = "build/lotka-volterra-2d-singular",
         doctype = "github",
         mod = LotkaVolterra2dSingularSRK)

weave("weave/lotka-volterra-2d-singular-srk.jmd",
         out_path = "build/lotka-volterra-2d-singular",
         doctype = "github",
         mod = LotkaVolterra2dSingularSRK)


text:
using Distributions

struct InfiniteSimulation{F<:Number,W<:World, P<:Plotter} <: Simulation
    time::Base.RefValue{F}
    t_sheep_boredom::Base.RefValue{F}
    dt::F
    world::W
    plotter::P
end

InfiniteSimulation(;
    time = 0.0,
    dt::F,
    world::W,
    t_sheep_boredom = rand(Exponential(world.meanSheepDiffusionTime / world.sheeps.current_sheep[])),
    plotter::P = VoidPlotter()
    ) where {F<:Number,W<:World,P<:Plotter} = InfiniteSimulation{F,W,P}( Ref(time), Ref(t_sheep_boredom), dt, world, plotter )

function runsim( sim::InfiniteSimulation )
    p = plot( sim.plotter, sim.world, sim.time[] )
    io = IOBuffer()
    print(io, "\e[25l")
    @sync begin
    is_running = Base.RefValue(true)
    @async begin
        readavailable(stdin)
        is_running[] = false
        return nothing
    end # async
    @async begin
        while is_running[]
            iterate!( sim )

            plot!( io, p, sim.plotter, sim.world, sim.time[] )
        end # while
    end # async
    end # sync
    print(io, "\e[25h")
    endtime = sim.time[]
    sim.time[] = zero(sim.time[])
    return endtime
end # function


text:
# ADMM tracking example
# Aggregator node in Julia
#
# This file is part of the openBuildNet simulation framework developed at EPFL.
#
# Authors: Truong X. Nghiem (xuan.nghiem@epfl.ch)

function print_help()
  helpstring ="""aggregator.jl <N>

  where <N> is the number of local nodes.
  The environment variable OBNMQTT can be set to the address of the MQTT server.
  If it is not defined, the localhost address tcp://localhost:1883 will be used."""
  println(helpstring)
end

if length(ARGS) < 1
  println("Not enough input arguments.")
  print_help()
  exit(1)
end

Nnode = tryparse(Int64, ARGS[1])
if isnull(Nnode) || Nnode.value < 1
  println("Invalid number of local nodes.")
  print_help()
  exit(2)
end
Nnode = Nnode.value

using openbuildnet

const WORKSPACE = "disttracking"
const MAIN_UPDATE = 0
const maxIter = 50
const epsPrimal = 1e-5 # primal convergence criterion
const epsDual = 1e-2 # dual convergence criterion

# parameters and internal variables
Nhorizon = 5
rho = 10000

totalSteps = 50
reference = (rand(totalSteps) - 0.5) * Nnode / 2

# The following variables should be initialized inside the INIT callback
lambda = zeros(Nhorizon,1)
err = zeros(Nhorizon,1)

refHorizon = zeros(Nhorizon)  # reference for the horizon

Gammas = fill(rho*eye(Nhorizon), Nnode)
deltas = fill(zeros(Nhorizon), Nnode)

Σ = eye(Nnode) + ones(Nnode, Nnode)   # Constant matrix to explicitly solve the primal update

step = 1
admm_events = falses(Nnode,1)
in_admm_iters = false     # true if and only if ADMM iterations are ongoing
iter = 0

# define and reset ADMM variables
function reset_admm_vars()
  global U, Z, Z_prev, MU
  U = zeros(Nhorizon,Nnode)
  Z = zeros(Nhorizon,Nnode)
  Z_prev = zeros(Nhorizon,Nnode)
  MU = zeros(Nhorizon,Nnode)
end

# Vector of input and output ports
InputsU = Vector{OBNInputAbstract}()
OutputsGamma = Vector{OBNOutputAbstract}()
OutputsDelta = Vector{OBNOutputAbstract}()

# Computation of one iteration at the master node
function iter_computation()
  global MU += (U - Z)

  # local reference update
  global Z_prev = Z
  global Z = (U + MU .+ (refHorizon - lambda)) / Σ
  global err = sum(Z,2) - refHorizon

  # global dual update
  global lambda += err

  # Write to the output port(s) to send values to the local nodes immediately.
  global Gammas, deltas
  for k = 1:Nnode
    # Gammas[k] = rho*eye(Nhorizon)
    deltas[k] = rho * (-Z[:,k] + MU[:,k])
    sendsync(OutputsGamma[k], Gammas[k])
    sendsync(OutputsDelta[k], deltas[k])
  end
end

# ADMM iteration function for an input port corresponding to local node #k
# It is called every time a local node sends U<i> values to the aggregator.
# This callback should perform the ADMM iteration whenever all U<i>'s have been received.
# It should also check the termination condition.
# The computation should only be performed if and only if in_admm_iters = true; and when it finishes, it must reset this variable to stop the ADMM iterations.
function admmiter(port, k)
  global in_admm_iters, admm_events, U, iter

  if !in_admm_iters
    return;
  end

  # Mark the port event
  if admm_events[k]
    warn("During ADMM iteration, values from port U$k have been received multiple times.")
  end
  admm_events[k] = true
  U[:,k] = get(port)    # save the values

  # If all have been received, perform the iteration computation
  if all(admm_events)
    iter += 1

    fill!(admm_events, false) # reset tracking

    # Compute residual
    residualRef = sum(Z,2) - refHorizon
    convPrimal = sqrt(sum((U-Z).^2) + sum(residualRef.^2))
    convDual = sqrt(rho^2 * sum((Z - Z_prev).^2))

    # check convergence
    if ((convPrimal <= epsPrimal && convDual <= epsDual) || iter >= maxIter)
      println("Stop ADMM at iteration $iter: primal = $convPrimal, dual = $convDual")
      iter = 0
      in_admm_iters = false   # Stop iterations
      reset_admm_vars()
    else
      iter_computation()
    end
  end
end

# Create the node
node = OBNNode("Aggregator", WORKSPACE, get(ENV, "OBNMQTT", "tcp://localhost:1883"))

# Create individual input(s) for each local node to receive values from that node, for example x1 to xN as vectors.
for k = 1:Nnode
  push!(InputsU, create_input(node, "U$k", Vector{Float64}))
  push!(OutputsGamma, create_output(node, "Gamma$k", Matrix{Float64}))
  push!(OutputsDelta, create_output(node, "delta$k", Vector{Float64}))

  # Add callback on the U port
  on_receive(admmiter, InputsU[k], InputsU[k], k)
end

on_init(node) do
  println("INIT")
  reset_admm_vars()
end

on_term(node) do
  println("TERMINATE at $(sim_time(node,:s))")
end

on_block_output(node, MAIN_UPDATE) do
  # the current time, 0 if this is the initial step
  println("Start ADMM at time $(sim_time(node,:s)) (s)")

  global refHorizon = reference[step:(step-1+Nhorizon)]
  global step += 1

  # Start ADMM
  global in_admm_iters = true

  # Tell local nodes to start iterations
  iter_computation()

  # Loop to process all port events until it finishes, or a long timeout of 60s.
  while in_admm_iters
    if !process_port_events(node, 60)
      in_admm_iters = false
      stop(node, true)
      error("Timeout error while waiting for local nodes to update.")
    end
  end
end

println("Ready to start simulation...")
status = runsafe(node, 30)
println("Final status: ", status)
writecsv("reference.txt", reference)  # write the reference to a CSV file


text:
module TestOperators
using Test
using Operators
using LinearAlgebra

@testset "basic operators" begin
    # evaluation functional
    f = x->sin(sum(x))
    x = randn()
    δ = Operators.Evaluation(x)
    @test δ(f) ≈ f(x)

    # 2d evaluation functional test
    x = randn(2)
    δ = Operators.Evaluation(x)
    @test δ(f) ≈ f(x)
    @test Operators.islinear(δ)

    # sum
    δ = Operators.Evaluation.(x)
    δ = Operators.Sum(δ...)
    @test δ(f) ≈ sum(f.(x))

    # scaling
    δ = 2δ
    @test Operators.isfunctional(δ)
    @test Operators.islinear(δ)
    @test δ(f) ≈ 2sum(f.(x))

    # dot

    # shift
    x = randn()
    y = randn()
    S = Operators.Shift(x)
    δ = Operators.Evaluation(y)
    @test δ(S(f)) ≈ f(y-x)

# composition
# println(Operators.Composition(δ, S)(f))

end

@testset "scalar differential operators" begin
    f(x) = sin(sum(x))
    # del
    ∇ = Operators.del
    # derivative
    x = randn()
    δ = Operators.δ(x)
    @test δ(∇(f)) ≈ cos(x)
    # gradient
    x = randn(2)
    δ = Operators.Evaluation(x)
    @test  δ(∇(f)) ≈ cos(sum(x)) * ones(2)

    # laplacian
    # Δ = Operators.Laplacian()
    Δ = Operators.Δ
    x = randn()
    δ = Operators.Evaluation(x)
    @test δ(Δ(f)) ≈ -sin(x)

    x = randn(2)
    δ = Operators.Evaluation(x)
    @test  δ(Δ(f)) ≈ sum(-sin(sum(x)) * ones(2))
    # ODE operator
end

@testset "vector differential operators" begin

    g(x) = [sin(sum(x)), cos(sum(x))]
    Operators.outputvalue(::typeof(g)) = Operators.VectorValued()

    # vector del, scalar input
    ∇ = Operators.∇
    @test ∇(g)(0.) ≈ [1, 0]

    # vector del, vector input
    x = randn(2)
    @test ∇(g)(x) ≈ [[cos(sum(x)) cos(sum(x))]; [-sin(sum(x)) -sin(sum(x))]]

    # vector laplacian, scalar input
    Δ = Operators.Δ
    println(Δ(g)(1.))

    # vector laplacian, vector input
    # println(Δ(g)(randn(2)))

    # curl
    f(x, y, z) = [y, -x, zero(z)]
    f(x::AbstractVector) = f(x...)
    x = randn(3)
    @test (∇ × f)(x) ≈ [0, 0, -2]

    f(x, y, z) = [zero(x), -x^2, zero(z)]
    f(x::AbstractVector) = f(x...)
    x = randn(3)
    @test (∇ × f)(x) ≈ [0, 0, -2x[1]]
end

end # TestOperators


text:
workspace()
using HoTRG_lemon.SimulatorModule
using HoTRG_lemon.SpinModule
using HoTRG_lemon.LatticeModule
using TensorOperations
using TensorMatrices_lemon

dimM = 5
inititer = 0
trotterparameter = 0.01
trotteriteration = 40

isinginfo = SpinInfo("quantum_ising_2",1.0e-13,0.0)
trotterinfo = TrotterInfo(trotterparameter, trotteriteration)
ising = QuantumIsingModel(isinginfo, trotterinfo)
lattice = buildLattice("quantum_2d_fractal", ising)
simulator = buildSimulator(lattice,dimM, inititer)

fieldrange = linspace(0.1, 3.5, 20)
simulatorQuantum(fieldrange, simulator; filename = "data_dim" * string(dimM)  * "_init" *  string(inititer) * ".txt")


text:
abstract RegERM
abstract RegressionSolver

# FIX: could not use ``method`` as keyword directly due to ``invoke` in RidgeReg, see https://github.com/JuliaLang/julia/issues/7045
optimize(method::RegERM; optimizer::Symbol=:l_bfgs) = optimize(method, optimizer)

function optimize(method::RegERM, optimizer::Symbol=:l_bfgs)
    check_hyperparameters(method.params)

    # init model
    model = Model(method.X, method.y, method.regression_type, method.kernel)

    if optimizer == :sgd
        model.theta = solve(model, method, SGDSolver(), method.X, method.y)
    elseif optimizer == :l1_rda && typeof(method.params) == L1RDAParameters
        model.theta = solve(model, method, L1RDASolver(), method.X, method.y)
    elseif optimizer == :l_bfgs
        model.theta = solve(model, method, LBFGSSolver(), method.X, method.y)
    else
        throw(ArgumentError("Unknown optimizer=$(optimizer) or mismatched hyperparameters"))
    end
    model
end

function check_arguments(X::Matrix, y::Vector, regression_type::Symbol)
    (n, m) = size(X)
    if (n != length(y))
        throw(DimensionMismatch("Dimensions of X and y mismatch."))
    end
    if (regression_type == :binomial)
        if (sort(unique(y)) != [-1,1])
            throw(ArgumentError("Class labels have to be either -1 or 1"))
        end
    elseif (regression_type == :multinomial)
        if (typeof(y) != Array{Int,1} || any(y.<1))
            throw(ArgumentError("Classes have to be positive integer values"))
        end
    elseif (regression_type == :ordinal)
        # something additional to check for regression?
    else
        throw(ArgumentError("Unknown regression_type=$(regression_type)"))
    end
end

# Pretty-print
function Base.show(io::IO, model::RegERM)
    println(io, "$(methodname(model))")
    println(io, repeat("-", length(methodname(model))))
    println(io, "number of examples:       $(model.n)")
    println(io, "number of features:       $(model.m)")
    println(io, "kernel function:          $(model.kernel)")
end


text:
"""
Create a System by parsing power-system parameters from a Matpower, PTI, or JSON file and do
some data checks.
"""
function parse_standard_files(file::String; kwargs...)::System

    # function `parse_file` is in pm_io/common.jl
    data = parse_file(file)

    # Check for at least one bus in input file
    if (length(data["bus"]) < 1)
        @error "There are no buses in this file"
    end

    # in pm2ps_parser.jl
    sys = pm2ps_dict(data; kwargs...)

    return sys

end


text:
@testset "interface" begin
    @test typeof(SurrogateModelOptim.Options()) == SurrogateModelOptim.Options
    @test_throws MethodError SurrogateModelOptim.Options(not_a_valid_arugment=true)
    @test_throws AssertionError SurrogateModelOptim.Options(smooth=true)
    @test typeof(SurrogateModelOptim.Options(smooth=false)) == SurrogateModelOptim.Options
    @test typeof(SurrogateModelOptim.Options(smooth=:variable)) == SurrogateModelOptim.Options
    @test typeof(SurrogateModelOptim.Options(smooth=:single)) == SurrogateModelOptim.Options
    @test typeof(SurrogateModelOptim.Options(smooth=:single_user)) == SurrogateModelOptim.Options
end


text:
using ArgParse, PackageCompiler

function parse_commandline()
    s = ArgParseSettings(description="Builds system image for diffeqtorch")

    @add_arg_table! s begin
        "--replace"
            help = "Build system image, replacing default"
            action = :store_true
        "--restore"
            help = "Restore default system image"
            action = :store_true
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()

    packages = [:DifferentialEquations, :DiffEqSensitivity]

    if parsed_args["restore"]
        println("Restore default system image")
        PackageCompiler.restore_default_sysimage()
    elseif parsed_args["replace"]
        println("Build system image, replace default")
        create_sysimage(packages, precompile_execution_file=string(@__DIR__, "/precompile.jl"), replace_default=true)
    else
        println("Build system image")
        create_sysimage(packages, sysimage_path=ENV["JULIA_SYSIMAGE_DIFFEQTORCH"], precompile_execution_file=string(@__DIR__, "/precompile.jl"))
    end

end

main()

max_stars_repo_name: michaeldeistler/diffeqtorch

max_issues_repo_name: michaeldeistler/diffeqtorch

max_forks_repo_name: michaeldeistler/diffeqtorch


text:
module MapApplyTests

using Test
using Gridap
using Gridap.CachedArrays
using TensorValues

using ..MapsMocks

a = VectorValue(10,10)
b = VectorValue(15,20)
p1 = VectorValue(1,1)
p2 = VectorValue(2,2)
p3 = VectorValue(3,3)
p = [p1,p2,p3]
m = MockMap(a)

m2 = apply(-,m,broadcast=true)

r = evaluate(m,p)
r2 = -r
test_map(m2,p,r2)

m3 = apply(-,m,m2,broadcast=true)
r3 = r .- r2
test_map(m3,p,r3)

m3 = apply(-,m,p,broadcast=true)
r3 = r .- p
test_map(m3,p,r3)

a = [1.0,2.0,3.0]
m3 = apply(outer,m,a,broadcast=true)
r3 = broadcast(outer,r,a)
test_map(m3,p,r3)

a = 2.0
m3 = apply(outer,m,a,broadcast=true)
r3 = broadcast(outer,r,a)
test_map(m3,p,r3)

a = VectorValue(10.1,10.3)
m = MockMap(a)
r = evaluate(m,p)
m2 = apply(-,m,broadcast=true)
r2 = -r
test_map(m2,p,r2)

end # module


text:
@muladd function perform_step!(integrator,cache::TauLeapingConstantCache)
  @unpack t,dt,uprev,u,W,p,P,c = integrator
  tmp = c(uprev, p, t, P.dW, nothing)
  integrator.u = uprev .+ tmp

  if integrator.opts.adaptive
    if typeof(integrator.alg) <: TauLeaping
      oldrate = P.cache.currate
      newrate = P.cache.rate(integrator.u,p,t+dt)
      EEstcache = @. abs(newrate - oldrate) / max(50integrator.opts.reltol*oldrate,integrator.rate_constants/integrator.dt)
      integrator.EEst = maximum(EEstcache)
      if integrator.EEst <= 1
        P.cache.currate = newrate
      end
    elseif typeof(integrator.alg) <: CaoTauLeaping
      # Calculate τ as EEst
    end
  end
end

@muladd function perform_step!(integrator,cache::TauLeapingCache)
  @unpack t,dt,uprev,u,W,p,P,c = integrator
  @unpack tmp, newrate, EEstcache = cache
  c(tmp, uprev, p, t, P.dW, nothing)
  @.. u = uprev + tmp

  if integrator.opts.adaptive
    if typeof(integrator.alg) <: TauLeaping
      oldrate = P.cache.currate
      P.cache.rate(newrate,u,p,t+dt)
      @.. EEstcache =  abs(newrate - oldrate) / max(50integrator.opts.reltol*oldrate,integrator.rate_constants/integrator.dt)
      integrator.EEst = maximum(EEstcache)
      if integrator.EEst <= 1
        P.cache.currate .= newrate
      end
    elseif typeof(integrator.alg) <: CaoTauLeaping
      # Calculate τ as EEst
    end
  end
end

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq


text:
using CSV, DataFrames, Statistics

sa(str::AbstractString) = let (c, time) = readlines(`../bin/sa $str`)
    parse(Int64, c), parse(Float64, time)
end

df = CSV.File("data/str_sa.csv") |> DataFrame
select!(df, :string, :nLetters, :nSeqs, :SA, :nPaths, :runtime)
rename!(df, Dict(:runtime => :runtime_sa))

@time asa = sa.(df.string)

df.ASA = first.(asa)
df.runtime_asa = last.(asa)

δ = df.SA - df.ASA
@show extrema(δ)
@show mean(δ) std(δ)

CSV.write("data/sa_asa.csv", df)

max_stars_repo_name: dglmoore/pathways

max_issues_repo_name: dglmoore/pathways

max_forks_repo_name: dglmoore/pathways


text:
abstract type LoadCategory end
struct ConstantCurrent <: LoadCategory end
struct ConstantImpedance <: LoadCategory end
struct ConstantPower <: LoadCategory end

const LOAD_MAP = Dict(
    PSY.LoadModels.ConstantCurrent => ConstantCurrent,
    PSY.LoadModels.ConstantImpedance => ConstantImpedance,
    PSY.LoadModels.ConstantPower => ConstantPower,
)


text:
include("triangle.jl");
include("genmesh.jl");
include("assemble_BEM.jl");
import .NEPCore.compute_Mder
import .NEPCore.size
import .NEPCore.compute_Mlincomb

using LinearAlgebra
"""
    struct BEM_NEP <: NEP

Represents a boundary element discretization. Create with a `mesh::Vector{Triangle}` using `nep=BEM_NEP(mesh)`.

"""
struct BEM_NEP <: NEP
    mesh::Vector{Triangle}
    n::Int
    gauss_order::Int
end

function BEM_NEP(mesh::Vector,gauss_order=3)
    return BEM_NEP(mesh,size(mesh,1),gauss_order);
end

function size(nep::BEM_NEP)
    return (nep.n,nep.n);
end
function size(nep::BEM_NEP,dim)
    return nep.n;
end
#
function compute_Mder(nep::BEM_NEP,λ::Number,der::Int=0)
    return assemble_BEM(λ, nep.mesh, nep.gauss_order,der)[:,:,1];
end
# Delegate the compute Mlincomb functions. This is quite inefficient.
compute_Mlincomb(nep::BEM_NEP,λ::Number,V::Union{AbstractMatrix,AbstractVector}, a::Vector) = compute_Mlincomb_from_Mder(nep,λ,V,a)
compute_Mlincomb(nep::BEM_NEP,λ::Number,V::Union{AbstractMatrix,AbstractVector}) = compute_Mlincomb(nep,λ,V, ones(eltype(V),size(V,2)))

function bem_fichera(N::Int=3)
    mesh=gen_ficheramesh(N)
    nep=BEM_NEP(mesh);
    precompute_quad!(nep.mesh,nep.gauss_order)
    return nep;
end


text:
module TestReshape

using Test, DataFrames, Random, Logging, PooledArrays
const ≅ = isequal

@testset "the output of unstack" begin
    df = DataFrame(Fish = CategoricalArray{Union{String, Missing}}(["Bob", "Bob", "Batman", "Batman"]),
                   Key = CategoricalArray{Union{String, Missing}}(["Mass", "Color", "Mass", "Color"]),
                   Value = Union{String, Missing}["12 g", "Red", "18 g", "Grey"])
    # Check that reordering levels does not confuse unstack
    levels!(df[!, 1], ["XXX", "Bob", "Batman"])
    levels!(df[!, 2], ["YYY", "Color", "Mass"])
    #Unstack specifying a row column
    df2 = unstack(df, :Fish, :Key, :Value)
    @test levels(df[!, 1]) == ["XXX", "Bob", "Batman"] # make sure we did not mess df[!, 1] levels
    @test levels(df[!, 2]) == ["YYY", "Color", "Mass"] # make sure we did not mess df[!, 2] levels
    #Unstack without specifying a row column
    df3 = unstack(df, :Key, :Value)
    #The expected output, XXX level should be dropped as it has no rows with this key
    df4 = DataFrame(Fish = Union{String, Missing}["Bob", "Batman"],
                    Color = Union{String, Missing}["Red", "Grey"],
                    Mass = Union{String, Missing}["12 g", "18 g"])
    @test df2 ≅ df4
    @test typeof(df2[!, :Fish]) <: CategoricalVector{Union{String, Missing}}
    # first column stays as CategoricalArray in df3
    @test df3 == df4
    #Make sure unstack works with missing values at the start of the value column
    df[1,:Value] = missing
    df2 = unstack(df, :Fish, :Key, :Value)
    #This changes the expected result
    df4[1,:Mass] = missing
    @test df2 ≅ df4

    df = DataFrame(Fish = CategoricalArray{Union{String, Missing}}(["Bob", "Bob", "Batman", "Batman"]),
                   Key = CategoricalArray{Union{String, Missing}}(["Mass", "Color", "Mass", "Color"]),
                   Value = Union{String, Missing}["12 g", "Red", "18 g", "Grey"])
    levels!(df[!, 1], ["XXX", "Bob", "Batman"])
    levels!(df[!, 2], ["YYY", "Color", "Mass"])
    df2 = unstack(df, :Fish, :Key, :Value, renamecols=x->string("_", uppercase(x), "_"))
    df3 = unstack(df, :Key, :Value, renamecols=x->string("_", uppercase(x), "_"))
    df4 = DataFrame(Fish = Union{String, Missing}["Bob", "Batman"],
                    _COLOR_ = Union{String, Missing}["Red", "Grey"],
                    _MASS_ = Union{String, Missing}["12 g", "18 g"])
    @test df2 == df4
    @test df3 == df4

    #The same as above but without CategoricalArray
    df = DataFrame(Fish = ["Bob", "Bob", "Batman", "Batman"],
                   Key = ["Mass", "Color", "Mass", "Color"],
                   Value = ["12 g", "Red", "18 g", "Grey"])
    #Unstack specifying a row column
    df2 = unstack(df, :Fish, :Key, :Value)
    #Unstack without specifying a row column
    df3 = unstack(df, :Key, :Value)
    #The expected output, XXX level should be dropped as it has no rows with this key
    df4 = DataFrame(Fish = ["Batman", "Bob"],
                    Color = ["Grey", "Red"],
                    Mass = ["18 g", "12 g"])
    @test df2 ≅ df4
    @test typeof(df2[!, :Fish]) <: Vector{String}
    # first column stays as CategoricalArray in df3
    @test df3 == df4
    #Make sure unstack works with missing values at the start of the value column
    allowmissing!(df, :Value)
    df[1,:Value] = missing
    df2 = unstack(df, :Fish, :Key, :Value)
    #This changes the expected result
    allowmissing!(df4, :Mass)
    df4[2,:Mass] = missing
    @test df2 ≅ df4

    df = DataFrame(Fish = ["Bob", "Bob", "Batman", "Batman"],
                   Key = ["Mass", "Color", "Mass", "Color"],
                   Value = ["12 g", "Red", "18 g", "Grey"])
    df2 = unstack(df, :Fish, :Key, :Value, renamecols=x->string("_", uppercase(x), "_"))
    df3 = unstack(df, :Key, :Value, renamecols=x->string("_", uppercase(x), "_"))
    df4 = DataFrame(Fish = ["Batman", "Bob"],
                    _COLOR_ = ["Grey", "Red"],
                    _MASS_ = ["18 g", "12 g"])
    @test df2 == df4
    @test df3 == df4

    # test empty set of grouping variables
    @test_throws ArgumentError unstack(df, Int[], :Key, :Value)
    @test_throws ArgumentError unstack(df, r"xxxxx", :Key, :Value)
    @test_throws ArgumentError unstack(df, Symbol[], :Key, :Value)
    @test_throws ArgumentError unstack(stack(DataFrame(rand(10, 10))),
                                  :id, :variable, :value)
    @test_throws TypeError unstack(df, :Key, :Value, renamecols=Symbol)

    # test missing value in grouping variable
    mdf = DataFrame(id=[missing,1,2,3], a=1:4, b=1:4)
    @test unstack(stack(mdf, Not(:id)), :id, :variable, :value)[1:3,:] == sort(mdf)[1:3,:]
    @test unstack(stack(mdf, Not(1)), :id, :variable, :value)[1:3,:] == sort(mdf)[1:3,:]
    @test unstack(stack(mdf, Not(:id)), :id, :variable, :value)[:, 2:3] == sort(mdf)[:, 2:3]
    @test unstack(stack(mdf, Not(1)), :id, :variable, :value)[:, 2:3] == sort(mdf)[:, 2:3]

    # test more than one grouping column
    wide = DataFrame(id = 1:12,
                     a  = repeat([1:3;], inner = [4]),
                     b  = repeat([1:4;], inner = [3]),
                     c  = randn(12),
                     d  = randn(12))
    w2 = wide[:, [1, 2, 4, 5]]
    rename!(w2, [:id, :a, :_C_, :_D_])
    long = stack(wide)
    wide3 = unstack(long, [:id, :a], :variable, :value)
    @test wide3 == wide[:, [1, 2, 4, 5]]
    wide3 = unstack(long, [:id, :a], :variable, :value,
                    renamecols=x->string("_", uppercase(string(x)), "_"))
    @test wide3 == w2

    wide3 = unstack(long, r"^[ia]", :variable, :value)
    @test wide3 == wide[:, [1, 2, 4, 5]]
    wide3 = unstack(long, r"^[ia]", :variable, :value,
                    renamecols=x->string("_", uppercase(string(x)), "_"))
    @test wide3 == w2
end

@testset "unstack promotion to support missing values" begin
    df = DataFrame([repeat(1:2, inner=4), repeat('a':'d', outer=2), collect(1:8)],
                   [:id, :variable, :value])
    udf = unstack(df, :variable, :value)
    @test udf == unstack(df, :variable, :value) == unstack(df, :id, :variable, :value)
    @test udf == DataFrame([Union{Int, Missing}[1, 2], Union{Int, Missing}[1, 5],
                            Union{Int, Missing}[2, 6], Union{Int, Missing}[3, 7],
                            Union{Int, Missing}[4, 8]], [:id, :a, :b, :c, :d])
    @test isa(udf[!, 1], Vector{Int})
    @test all(isa.(eachcol(udf)[2:end], Vector{Union{Int, Missing}}))
    df = DataFrame([categorical(repeat(1:2, inner=4)),
                       categorical(repeat('a':'d', outer=2)), categorical(1:8)],
                   [:id, :variable, :value])
    udf = unstack(df, :variable, :value)
    @test udf == unstack(df, :variable, :value) == unstack(df, :id, :variable, :value)
    @test udf == DataFrame([Union{Int, Missing}[1, 2], Union{Int, Missing}[1, 5],
                            Union{Int, Missing}[2, 6], Union{Int, Missing}[3, 7],
                            Union{Int, Missing}[4, 8]], [:id, :a, :b, :c, :d])
    @test isa(udf[!, 1], CategoricalVector{Int})
    @test all(isa.(eachcol(udf)[2:end], CategoricalVector{Union{Int, Missing}}))
end

@testset "duplicate entries in unstack warnings" begin
    df = DataFrame(id=Union{Int, Missing}[1, 2, 1, 2],
                   id2=Union{Int, Missing}[1, 2, 1, 2],
                   variable=["a", "b", "a", "b"], value=[3, 4, 5, 6])
    @test_logs (:warn, "Duplicate entries in unstack at row 3 for key 1 and variable a.") unstack(df, :id, :variable, :value)
    @test_logs (:warn, "Duplicate entries in unstack at row 3 for key (1, 1) and variable a.") unstack(df, :variable, :value)
    a, b = with_logger(NullLogger()) do
        unstack(df, :id, :variable, :value), unstack(df, :variable, :value)
    end
    @test a ≅ DataFrame(id = [1, 2], a = [5, missing], b = [missing, 6])
    @test b ≅ DataFrame(id = [1, 2], id2 = [1, 2], a = [5, missing], b = [missing, 6])

    df = DataFrame(id=1:2, variable=["a", "b"], value=3:4)
    @test_nowarn unstack(df, :id, :variable, :value)
    @test_nowarn unstack(df, :variable, :value)
    a = unstack(df, :id, :variable, :value)
    b = unstack(df, :variable, :value)
    @test a ≅ b ≅ DataFrame(id = [1, 2], a = [3, missing], b = [missing, 4])

    df = DataFrame(variable=["x", "x"], value=[missing, missing], id=[1,1])
    @test_logs (:warn, "Duplicate entries in unstack at row 2 for key 1 and variable x.") unstack(df, :variable, :value)
    @test_logs (:warn, "Duplicate entries in unstack at row 2 for key 1 and variable x.") unstack(df, :id, :variable, :value)
end

@testset "missing values in colkey" begin
    df = DataFrame(id=[1, 1, 1, missing, missing, missing, 2, 2, 2],
                   variable=["a", "b", missing, "a", "b", "missing", "a", "b", "missing"],
                   value=[missing, 2.0, 3.0, 4.0, 5.0, missing, 7.0, missing, 9.0])
    @test_logs (:warn, "Missing value in variable variable at row 3. Skipping.") unstack(df, :variable, :value)
    udf = with_logger(NullLogger()) do
        unstack(df, :variable, :value)
    end
    @test names(udf) == [:id, :a, :b, :missing]
    @test udf[!, :missing] ≅ [missing, 9.0, missing]
    df = DataFrame(id=[1, 1, 1, missing, missing, missing, 2, 2, 2],
                   id2=[1, 1, 1, missing, missing, missing, 2, 2, 2],
                   variable=["a", "b", missing, "a", "b", "missing", "a", "b", "missing"],
                   value=[missing, 2.0, 3.0, 4.0, 5.0, missing, 7.0, missing, 9.0])
    @test_logs (:warn, "Missing value in variable variable at row 3. Skipping.") unstack(df, 3, 4)
    udf = with_logger(NullLogger()) do
        unstack(df, 3, 4)
    end
    @test names(udf) == [:id, :id2, :a, :b, :missing]
    @test udf[!, :missing] ≅ [missing, 9.0, missing]
end

@testset "stack-unstack correctness" begin
    Random.seed!(1234)
    x = DataFrame(rand(100, 50))
    x[!, :id] = [1:99; missing]
    x[!, :id2] = string.("a", x[!, :id])
    x[!, :s] = [i % 2 == 0 ? randstring() : missing for i in 1:100]
    allowmissing!(x, :x1)
    x[1, :x1] = missing
    y = stack(x, Not([:id, :id2]))
    @test y ≅ stack(x, Not(r"id"))
    @test y ≅ stack(x, Not(Not(Not(r"id"))))
    z = unstack(y, :id, :variable, :value)
    @test all(isequal(z[!, n], x[!, n]) for n in names(z))
    z = unstack(y, :variable, :value)
    @test all(isequal(z[!, n], x[!, n]) for n in names(x))
end

@testset "reshape" begin
    d1 = DataFrame(a = Array{Union{Int, Missing}}(repeat([1:3;], inner = [4])),
                b = Array{Union{Int, Missing}}(repeat([1:4;], inner = [3])),
                c = Array{Union{Float64, Missing}}(randn(12)),
                d = Array{Union{Float64, Missing}}(randn(12)),
                e = Array{Union{String, Missing}}(map(string, 'a':'l')))

    @test names(stack(d1, :a)) == [:b, :c, :d, :e, :variable, :value]
    d1s = stack(d1, [:a, :b])
    @test d1s == stack(d1, r"[ab]")
    @test d1s == stack(d1, Not(r"[cde]"))
    @test d1s == stack(d1, Not(Not(r"[ab]")))
    d1s2 = stack(d1, [:c, :d])
    @test d1s2 == stack(d1, r"[cd]")
    @test d1s2 == stack(d1, Not([1, 2, 5]))
    d1s3 = stack(d1)
    d1m = stack(d1, Not([:c, :d, :e]))
    @test d1m == stack(d1, Not(r"[cde]"))
    @test d1s[1:12, :c] == d1[!, :c]
    @test d1s[13:24, :c] == d1[!, :c]
    @test d1s2 == d1s3
    @test names(d1s) == [:c, :d, :e, :variable, :value]
    @test d1s == d1m
    d1m = stack(d1[:, [1,3,4]], Not(:a))
    @test names(d1m) == [:a, :variable, :value]

    # Test naming of measure/value columns
    d1s_named = stack(d1, [:a, :b], variable_name=:letter, value_name=:someval)
    @test d1s_named == stack(d1, r"[ab]", variable_name=:letter, value_name=:someval)
    @test names(d1s_named) == [:c, :d, :e, :letter, :someval]
    d1m_named = stack(d1[:, [1,3,4]], Not(:a), variable_name=:letter, value_name=:someval)
    @test names(d1m_named) == [:a, :letter, :someval]

    # test empty measures or ids
    dx = stack(d1, [], [:a])
    @test dx == stack(d1, r"xxx", r"a")
    @test size(dx) == (0, 3)
    @test names(dx) == [:a, :variable, :value]
    dx = stack(d1, :a, [])
    @test dx == stack(d1, r"a", r"xxx")
    @test size(dx) == (12, 2)
    @test names(dx) == [:variable, :value]
    dx = stack(d1, [:a], [])
    @test dx == stack(d1, r"a", r"xxx")
    @test size(dx) == (12, 2)
    @test names(dx) == [:variable, :value]
    dx = stack(d1, [], :a)
    @test dx == stack(d1, r"xxx", r"a")
    @test size(dx) == (0, 3)
    @test names(dx) == [:a, :variable, :value]

    @test stack(d1, :a, view=true) == stack(d1, [:a], view=true)
    @test all(isa.(eachcol(stack(d1, :a, view=true)),
                   [fill(DataFrames.RepeatedVector, 5);
                    DataFrames.StackedVector]))
    @test all(isa.(eachcol(stack(d1, Not([:b, :c, :d, :e]), view=true)),
                   [fill(DataFrames.RepeatedVector, 5);
                    DataFrames.StackedVector]))

    # Tests of RepeatedVector and StackedVector indexing
    d1s = stack(d1, [:a, :b], view=true)
    @test d1s == stack(d1, r"[ab]", view=true)
    @test d1s[!, 4] isa DataFrames.RepeatedVector
    @test ndims(d1s[!, 4]) == 1
    @test ndims(typeof(d1s[!, 4])) == 1
    @test d1s[!, 5] isa DataFrames.StackedVector
    @test ndims(d1s[!, 5]) == 1
    @test ndims(typeof(d1s[!, 2])) == 1
    @test d1s[!, 4][[1,24]] == ["a", "b"]
    @test d1s[!, 5][[1,24]] == [1, 4]
    @test_throws ArgumentError d1s[!, 4][true]
    @test_throws ArgumentError d1s[!, 5][true]
    @test_throws ArgumentError d1s[!, 4][1.0]
    @test_throws ArgumentError d1s[!, 5][1.0]

    d1ss = stack(d1, [:a, :b], view=true)
    @test d1ss[!, 4][[1,24]] == ["a", "b"]
    @test d1ss[!, 4] isa DataFrames.RepeatedVector
    d1ss = stack(d1, [:a, :b], view=true, variable_eltype=String)
    @test d1ss[!, 4][[1,24]] == ["a", "b"]
    @test d1ss[!, 4] isa DataFrames.RepeatedVector
    d1ss = stack(d1, [:a, :b], view=true, variable_eltype=Symbol)
    @test d1ss[!, 4][[1,24]] == [:a, :b]
    @test d1ss[!, 4] isa DataFrames.RepeatedVector

    # Those tests check indexing RepeatedVector/StackedVector by a vector
    @test d1s[!, 4][trues(24)] == d1s[!, 4]
    @test d1s[!, 5][trues(24)] == d1s[!, 5]
    @test d1s[!, 4][:] == d1s[!, 4]
    @test d1s[!, 5][:] == d1s[!, 5]
    @test d1s[!, 4][1:24] == d1s[!, 4]
    @test d1s[!, 5][1:24] == d1s[!, 5]
    @test [d1s[!, 4][1:12]; d1s[!, 4][13:24]] == d1s[!, 4]
    @test [d1s[!, 5][1:12]; d1s[!, 5][13:24]] == d1s[!, 5]

    d1s2 = stack(d1, [:c, :d], view=true)
    @test d1s2 == stack(d1, r"[cd]", view=true)
    d1s3 = stack(d1, view=true)
    d1m = stack(d1, Not([:c, :d, :e]), view=true)
    @test d1m == stack(d1, Not(r"[cde]"), view=true)
    @test d1s[1:12, :c] == d1[!, :c]
    @test d1s[13:24, :c] == d1[!, :c]
    @test d1s2 == d1s3
    @test names(d1s) == [:c, :d, :e, :variable, :value]
    @test d1s == d1m
    d1m = stack(d1[:, [1,3,4]], Not(:a), view=true)
    @test names(d1m) == [:a, :variable, :value]

    d1s_named = stack(d1, [:a, :b], variable_name=:letter, value_name=:someval, view=true)
    @test d1s_named == stack(d1, r"[ab]", variable_name=:letter, value_name=:someval, view=true)
    @test names(d1s_named) == [:c, :d, :e, :letter, :someval]
    d1m_named = stack(d1, Not([:c, :d, :e]), variable_name=:letter, value_name=:someval, view=true)
    @test d1m_named == stack(d1, Not(r"[cde]"), variable_name=:letter, value_name=:someval, view=true)
    @test names(d1m_named) == [:c, :d, :e, :letter, :someval]

    d1s[!, :id] = Union{Int, Missing}[1:12; 1:12]
    d1s2[!, :id] =  Union{Int, Missing}[1:12; 1:12]
    d1us = unstack(d1s, :id, :variable, :value)
    d1us2 = unstack(d1s2, :id, :variable, :value)
    d1us3 = unstack(d1s2, :variable, :value)
    @test d1us[!, :a] == d1[!, :a]
    @test d1us2[!, :d] == d1[!, :d]
    @test d1us2[!, 3] == d1[!, :d]
    @test d1us3[!, :d] == d1[!, :d]
    @test d1us3 == unstack(d1s2)

    # test unstack with exactly one key column that is not passed
    df1 = stack(DataFrame(rand(10,10)))
    df1[!, :id] = 1:100
    @test size(unstack(df1, :variable, :value)) == (100, 11)
    @test unstack(df1, :variable, :value) ≅ unstack(df1)

    # test empty keycol
    @test_throws ArgumentError unstack(stack(DataFrame(rand(3,2))), :variable, :value)
end

@testset "column names duplicates" begin
    wide = DataFrame(id = 1:12,
                     a  = repeat([1:3;], inner = [4]),
                     b  = repeat([1:4;], inner = [3]),
                     c  = randn(12),
                     d  = randn(12))
    w1 = wide[:, [1, 2, 4, 5]]
    rename!(w1, [:id, :a, :D, :d])
    w2 = wide[:, [1, 2, 4, 5]]
    rename!(w2, [:id, :a, :_D_, :_d_])
    long = stack(wide)
    long.variable = replace.(string.(long.variable), Ref("c" => "D"))
    wide3 = unstack(long, [:id, :a], :variable, :value)
    @test wide3 == w1
    wide3 = unstack(long, :id, :variable, :value)
    @test wide3 == w1[!, Not(2)]
    @test_throws ArgumentError unstack(long, [:id, :a], :variable, :value,
                                       renamecols=uppercase)
    wide3 = unstack(long, [:id, :a], :variable, :value,
                    renamecols=x->string("_", x, "_"))
    @test wide3 == w2
    wide3 = unstack(long, :id, :variable, :value,
                    renamecols=x->string("_", x, "_"))
    @test wide3 == w2[!, Not(2)]
end

@testset "flatten single column" begin
    df_vec = DataFrame(a = [1, 2], b = [[1, 2], [3, 4]])
    df_tup = DataFrame(a = [1, 2], b = [(1, 2), (3, 4)])
    ref = DataFrame(a = [1, 1, 2, 2], b = [1, 2, 3, 4])
    @test flatten(df_vec, :b) == flatten(df_tup, :b) == ref
    df_mixed_types = DataFrame(a = [1, 2], b = [[1, 2], ["x", "y"]])
    ref_mixed_types = DataFrame(a = [1, 1, 2, 2], b = [1, 2, "x", "y"])
    @test flatten(df_mixed_types, :b) == ref_mixed_types
    df_three = DataFrame(a = [1, 2, 3], b = [[1, 2], [10, 20], [100, 200, 300]])
    ref_three = DataFrame(a = [1, 1, 2, 2, 3, 3, 3], b = [1, 2, 10, 20, 100, 200, 300])
    @test flatten(df_three, :b) == ref_three
    df_gen = DataFrame(a = [1, 2], b = [(i for i in 1:5), (i for i in 6:10)])
    ref_gen = DataFrame(a = [fill(1, 5); fill(2, 5)], b = collect(1:10))
    @test flatten(df_gen, :b) == ref_gen
    df_miss = DataFrame(a = [1, 2], b = [Union{Missing, Int}[1, 2], Union{Missing, Int}[3, 4]])
    ref = DataFrame(a = [1, 1, 2, 2], b = [1, 2, 3, 4])
    @test flatten(df_miss, :b) == ref
    v1 = [[1, 2], [3, 4]]
    v2 = [[5, 6], [7, 8]]
    v = [v1, v2]
    df_vec_vec = DataFrame(a = [1, 2], b = v)
    ref_vec_vec = DataFrame(a = [1, 1, 2, 2], b = [v1 ; v2])
    @test flatten(df_vec_vec, :b) == ref_vec_vec
    df_cat = DataFrame(a = [1, 2], b = [CategoricalArray([1, 2]), CategoricalArray([1, 2])])
    df_flat_cat = flatten(df_cat, :b)
    ref_cat = DataFrame(a = [1, 1, 2, 2], b = [1, 2, 1, 2])
    @test df_flat_cat == ref_cat
    @test df_flat_cat.b isa CategoricalArray
end

@testset "flatten multiple columns" begin
    df = DataFrame(a = [1, 2], b = [[1, 2], [3, 4]], c = [[5, 6], [7, 8]])
    @test flatten(df, []) == df
    ref = DataFrame(a = [1, 1, 2, 2], b = [1, 2, 3, 4], c = [5, 6, 7, 8])
    @test flatten(df, [:b, :c]) == ref
    @test flatten(df, [:c, :b]) == ref
    @test flatten(df, 2:3) == ref
    @test flatten(df, r"[bc]") == ref
    @test flatten(df, Not(:a)) == ref
    @test flatten(df, Between(:b, :c)) == ref
    df_allcols = DataFrame(b = [[1, 2], [3, 4]], c = [[5, 6], [7, 8]])
    ref_allcols = DataFrame(b = [1, 2, 3, 4], c = [5, 6, 7, 8])
    @test flatten(df_allcols, All()) == ref_allcols
    @test flatten(df_allcols, :) == ref_allcols
    df_bad = DataFrame(a = [1, 2], b = [[1, 2], [3, 4]], c = [[5, 6], [7]])
    @test_throws ArgumentError flatten(df_bad, [:b, :c])
end

@testset "test RepeatedVector for categorical" begin
    v = categorical(["a", "b", "c"], ordered=true)
    levels!(v, ["b", "c", "a"])
    rv = DataFrames.RepeatedVector(v, 1, 1)
    @test isordered(v)
    # uncomment after CategoricalArrays.jl is fixed
    # @test isordered(categorical(v))
    @test levels(v) == ["b", "c", "a"]
    @test levels(categorical(v)) == ["b", "c", "a"]

    v = categorical(["a", "b", "c"])
    levels!(v, ["b", "c", "a"])
    rv = DataFrames.RepeatedVector(v, 1, 1)
    @test !isordered(v)
    # uncomment after CategoricalArrays.jl is fixed
    # @test !isordered(categorical(v))
    @test levels(v) == ["b", "c", "a"]
    @test levels(categorical(v)) == ["b", "c", "a"]
end

@testset "stack categorical test" begin
    Random.seed!(1234)
    d1 = DataFrame(a = repeat([1:3;], inner = [4]),
                   b = repeat([1:4;], inner = [3]),
                   c = randn(12),
                   d = randn(12),
                   e = map(string, 'a':'l'))
    d1s = stack(d1, [:d, :c])
    @test d1s.variable isa CategoricalVector{String}
    @test levels(d1s.variable) == ["d", "c"]
    d1s = stack(d1, [:d, :c], view=true)
    @test d1s.variable isa DataFrames.RepeatedVector{<:CategoricalString}
    @test levels(d1s.variable) == ["d", "c"]
    @test d1s[:, 4] isa CategoricalVector{String}
    @test levels(d1s[:, 4]) == ["d", "c"]

    d1s = stack(d1, [:d, :c], variable_eltype=String)
    @test d1s.variable isa PooledVector{String}
    @test levels(d1s.variable) == ["c", "d"]
    d1s = stack(d1, [:d, :c], view=true, variable_eltype=String)
    @test d1s.variable isa DataFrames.RepeatedVector{String}
    @test levels(d1s.variable) == ["c", "d"]
    @test d1s[:, 4] isa Vector{String}
    @test levels(d1s[:, 4]) == ["c", "d"]

    d1s = stack(d1, [:d, :c], variable_eltype=Symbol)
    @test d1s.variable isa Vector{Symbol}
    @test levels(d1s.variable) == [:c, :d]
    d1s = stack(d1, [:d, :c], view=true, variable_eltype=Symbol)
    @test d1s.variable isa DataFrames.RepeatedVector{Symbol}
    @test levels(d1s.variable) == [:c, :d]
    @test d1s[:, 4] isa Vector{Symbol}
    @test levels(d1s[:, 4]) == [:c, :d]

    d2 = categorical(d1, :)
    levels!(d2.a, [2, 1, 3])
    ordered!(d2.a, true)
    ref_levels = shuffle!(unique([levels(d2.c); levels(d2.d)]))
    levels!(d2.c, ref_levels)
    ordered!(d2.c, true)
    levels!(d2.d, ref_levels)
    ordered!(d2.d, true)
    d2s = stack(d2, [:d, :c])
    for col in eachcol(d2s)
        @test col isa CategoricalVector
    end
    @test levels(d2s.value) == ref_levels
    @test isordered(d2s.value)
    @test levels(d2.a) == levels(d2s.a)
    @test levels(d2.b) == levels(d2s.b)
    @test levels(d2.e) == levels(d2s.e)
    @test isordered(d2.a) == isordered(d2s.a)
    @test isordered(d2.b) == isordered(d2s.b)
    @test isordered(d2.e) == isordered(d2s.e)
end

@testset "test stack eltype" begin
    df = DataFrame(rand(4,5))
    sdf = stack(df)
    @test eltype(sdf.variable) <: CategoricalString
    @test eltype(typeof(sdf.variable)) <: CategoricalString
    @test eltype(sdf.value) <: Float64
    @test eltype(typeof(sdf.value)) <: Float64
    sdf2 = first(sdf, 3)
    @test eltype(sdf2.variable) <: CategoricalString
    @test eltype(typeof(sdf2.variable)) <: CategoricalString
    @test eltype(sdf2.value) <: Float64
    @test eltype(typeof(sdf2.value)) <: Float64
end

end # module


text:
# Read in from JLD
data, TTT, RRR, CCC, ZZ, DD, F_ϵ, F_u, s_init = jldopen("reference/tpf_main_inputs.jld", "r") do file
    read(file, "data"),
    read(file, "TTT"),
    read(file, "RRR"),
    read(file, "CCC"),
    read(file, "ZZ"),
    read(file, "DD"),
    read(file, "F_epsilon"),
    read(file, "F_u"),
    read(file, "s_init")
end

# Tune algorithm
tuning = Dict(:r_star => 2., :c_init => 0.3, :target_accept_rate => 0.4,
              :resampling_method => :systematic, :n_mh_steps => 1,
              :n_particles => 1000, :n_presample_periods => 0,
              :allout => true)

# Define Φ and Ψ (can't be saved to JLD)
Φ(s_t::AbstractVector{Float64}, ϵ_t::AbstractVector{Float64}) = TTT*s_t + RRR*ϵ_t + CCC
Ψ(s_t::AbstractVector{Float64}) = ZZ*s_t + DD

# Load in test inputs and outputs
test_file_inputs = load("reference/tpf_aux_inputs.jld")
test_file_outputs = load("reference/tpf_aux_outputs.jld")

## Correction and Associated Auxiliary Function Tests
φ_old = test_file_inputs["phi_old"]
norm_weights = test_file_inputs["norm_weights"]
coeff_terms = test_file_inputs["coeff_terms"]
log_e_1_terms = test_file_inputs["log_e_1_terms"]
log_e_2_terms = test_file_inputs["log_e_2_terms"]
inc_weights = test_file_inputs["inc_weights"]
HH = cov(F_u)
s_t_nontemp = test_file_inputs["s_t_nontemp"]

weight_kernel!(coeff_terms, log_e_1_terms, log_e_2_terms, φ_old, Ψ, data[:, 47], s_t_nontemp, det(HH), inv(HH);
               initialize = false, parallel = false)
φ_new = next_φ(φ_old, coeff_terms, log_e_1_terms, log_e_2_terms, length(data[:,47]), tuning[:r_star], 2)
correction!(inc_weights, norm_weights, φ_new, coeff_terms, log_e_1_terms, log_e_2_terms, length(data[:,47]))

@testset "Corection and Auxiliary Tests" begin
    @test coeff_terms[1] ≈ test_file_outputs["coeff_terms"][1]
    @test log_e_1_terms[1] ≈ test_file_outputs["log_e_1_terms"][1]
    @test log_e_2_terms[1] ≈ test_file_outputs["log_e_2_terms"][1]
    @test φ_new ≈ test_file_outputs["phi_new"]
    @test inc_weights[1] ≈ test_file_outputs["inc_weights"][1]
end

## Selection Tests
s_t1_temp = test_file_inputs["s_t1_temp"]
ϵ_t = test_file_inputs["eps_t"]

srand(47)
selection!(norm_weights, s_t1_temp, s_t_nontemp,ϵ_t, resampling_method = tuning[:resampling_method])
@testset "Selection Tests" begin
    @test s_t1_temp[1] ≈ test_file_outputs["s_t1_temp"][1]
    @test s_t_nontemp[1] ≈ test_file_outputs["s_t_nontemp"][1]
    @test ϵ_t[1] ≈ test_file_outputs["eps_t"][1]
end

## Mutation Tests
QQ = cov(F_ϵ)
accept_rate = test_file_inputs["accept_rate"]
c = test_file_inputs["c"]

c = update_c(c, accept_rate, tuning[:target_accept_rate])
srand(47)
mutation!(Φ, Ψ, QQ, det(HH), inv(HH), φ_new, data[:,47], s_t_nontemp, s_t1_temp, ϵ_t, c, tuning[:n_mh_steps])

@testset "Mutation Tests" begin
    @test s_t_nontemp[1] ≈ test_file_outputs["s_t_nontemp_mutation"][1]
    @test ϵ_t[1] ≈ test_file_outputs["eps_t_mutation"][1]
end

## Whole TPF Tests
srand(47)
out_no_parallel = tempered_particle_filter(data, Φ, Ψ, F_ϵ, F_u, s_init; tuning..., verbose = :none, parallel = false)
srand(47)
out_parallel_one_worker = tempered_particle_filter(data, Φ, Ψ, F_ϵ, F_u, s_init; tuning..., verbose = :none, parallel = true)
@testset "TPF tests" begin
    @test out_no_parallel[1] ≈ -302.99967306704133
    @test out_parallel_one_worker[1] ≈ -306.8211172094595
end


text:
#!/usr/bin/env julia

export pe2

function pe2(n::Integer=4*10^6)
    function fib(n::Integer)
        if n==1||n==2
            n
        else
            fib(n-1)+fib(n-2);
        end
    end
    i=1;
    while fib(i)<n
        i+=1;
    end
    filter(iseven,fib.(1:i-1))|>sum
end

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe2()|>print
end

max_stars_repo_path: 2. Even Fibonacci numbers/2.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 2. Even Fibonacci numbers/2.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 2. Even Fibonacci numbers/2.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
using JOLI
using Revise
using WAVEFORM
using LinearAlgebra

n = 100*[1;1;1];
d = 10.0*[1;1;1];
o = 0.0*[1;1;1];
v0 = 2000;
freq = 4.0;
λ = v0/freq;

t0 = 0.0;
f0 = 10.0;
unit = "m/s";
xsrc = [(n.*d)[1]/2];
ysrc = [(n.*d)[2]/2];
zsrc = [100.0];
xrec = range(0.0,stop=(n.*d)[1]::Float64,length=n[1]);
yrec = range(0.0,stop=(n.*d)[2]::Float64,length=n[2]);
zrec = [100.0];
freqs = [freq];
nsrc = length(xsrc)*length(ysrc)*length(zsrc);
nfreq = length(freqs)
model = Model{Int64,Float64}(n,d,o,t0,f0,unit,freqs,xsrc,ysrc,zsrc,xrec,yrec,zrec);

comp_n = n;
comp_d = d;
comp_o = o;
npml = convert(Int,λ/minimum(comp_d))*ones(Int64,(2,3));
npml = 30*ones(Int64,(2,3));
scheme = WAVEFORM.helm3d_operto27;
cut_pml = true;
implicit_matrix = true;
srcfreqmask = trues(nsrc,nfreq);
misfit = WAVEFORM.least_squares;
lsopts = LinSolveOpts(solver=:fgmres,maxinnerit=5,tol=1e-6);
lsopts.precond = :mlgmres;
opts = PDEopts{Int64,Float64}(scheme,comp_n,comp_d,comp_o,cut_pml,implicit_matrix,npml,misfit,srcfreqmask,lsopts);

nλ = maximum(n.*d/λ);
ppλ = λ/maximum(d);
v = v0*ones(Float64,n...);
v = vec(v);
opts.implicit_matrix = true;

(H,comp_grid,T,DT_adj,P) = helmholtz_system(v,model,freq,opts);

# ensure linear system is actually solved
nt = comp_grid.comp_n
q = zeros(eltype(H),tuple(nt...));
q[div.(nt,2)...] = 1.0
u = H\vec(q)
@test norm(H*u-vec(q)) < 1e-6*norm(q)
ut = H'\vec(q);
@test norm(H'*ut-vec(q)) < 1e-6*norm(q)

# test computed green's function ≈ analytical green's function in the non-pml region
(xt,yt,zt) = odn_to_grid(comp_grid.comp_o,comp_grid.comp_d,comp_grid.comp_n)
idx = findmax(abs.(q))[2]
qCartesian = CartesianIndices(q)
(ix,iy,iz) = qCartesian[idx].I
ω = (2*pi*freq)/v0;
R = [sqrt((x-xt[ix])^2 + (y-yt[ix])^2 + (z-zt[ix])^2) for x in xt, y in yt, z in zt];
G = [prod(model.d)*exp(im*ω*r)/(4*π*r) for r in R];
r = reshape(u,size(G))-G;
r[ix,iy,iz] = 0;
@test norm(comp_grid.comp_to_phys_grid*vec(r))/norm(comp_grid.comp_to_phys_grid*vec(u)) < 0.01

# test derivative behaviour
n = 50*[1;1;1];

xsrc = [(n.*d)[1]/2];
ysrc = [(n.*d)[2]/2];
zsrc = [100.0];
xrec = range(0.0,stop=(n.*d)[1]::Float64,length=n[1]);
yrec = range(0.0,stop=(n.*d)[2]::Float64,length=n[2]);
zrec = [100.0];
freqs = [freq];
nsrc = length(xsrc)*length(ysrc)*length(zsrc);
nfreq = length(freqs)
model = Model{Int64,Float64}(n,d,o,t0,f0,unit,freqs,xsrc,ysrc,zsrc,xrec,yrec,zrec);

comp_n = n;
comp_d = d;
comp_o = o;
npml = convert(Int,λ/minimum(comp_d))*ones(Int64,(2,3));
npml = 10*ones(Int64,(2,3));
scheme = WAVEFORM.helm3d_operto27;
cut_pml = true;
implicit_matrix = true;
srcfreqmask = trues(nsrc,nfreq);
misfit = WAVEFORM.least_squares;
lsopts = LinSolveOpts(solver=:fgmres,maxinnerit=5,tol=1e-10);
lsopts.precond = :mlgmres;
opts = PDEopts{Int64,Float64}(scheme,comp_n,comp_d,comp_o,cut_pml,implicit_matrix,npml,misfit,srcfreqmask,lsopts);

nλ = maximum(n.*d/λ);
ppλ = λ/maximum(d);
v = v0*ones(Float64,n...);
v = vec(v);
opts.implicit_matrix = true;

(H,comp_grid,T,DT_adj,P) = helmholtz_system(v,model,freq,opts);
nt = comp_grid.comp_n
q = zeros(eltype(H),tuple(nt...));
q[div.(nt,2)...] = 1.0;
q = vec(q);
u0 = H\q;
δv = randn(size(v));
δvx = comp_grid.phys_to_comp_grid*δv;
δu = H\(-T(u)*δvx);
step = 10.^(-5.0:1.0);
err0 = zeros(length(step))
err1 = zeros(length(step))
for i=1:length(step)
    h = step[i]
    H1 = helmholtz_system(v+δv*h,model,freq,opts)[1];
    u1 = H1\q;
    err0[i] = norm(u1-u0);
    err1[i] = norm(u1-u0-h*δu)
end

@test abs(median(diff(log10.(err0)))-1) < 0.2
@test abs(median(diff(log10.(err1)))-2) < 0.2