text:
# Test items in source.jl

file = joinpath(dir,"testfile.txt")

@testset "Source Generation" begin
    b = """
       abcd    10102017
       efgh123410102018"""

    @test_throws ArgumentError FWF.Source("", [1])
    @test_throws ArgumentError FWF.Source(file, Vector{Int}())
    @test_throws ArgumentError FWF.Source(file, [4,4,8], header=["A", "B"])
    @test_throws ArgumentError FWF.Source(file, [4,4,8], types=["A", "B"])
    @test_throws ArgumentError FWF.Source(file, [4,4,8], types=[Int64, 1234, String])
    @test_throws ArgumentError FWF.Source(file, [4,4,8], types=[String, Float32, DateFormat("mmddyyyy")])

    tmp = FWF.Source(file, [4,4,8])
    show(IOBuffer(), tmp) # Will retult in exception if problem unfortunately there is no @test_notthrows
    @test Data.header(tmp.schema)[1] == "Column1"
    @test Data.header(tmp.schema)[2] == "Column2"
    @test Data.header(tmp.schema)[3] == "Column3"

    tmp = FWF.Source(open(file, "r"), [4,4,8], header=true)
    @test Data.header(tmp.schema)[1] == "abcd"
    @test Data.header(tmp.schema)[2] == "1234"
    @test Data.header(tmp.schema)[3] == "10102017"

    tmp = FWF.Source(IOBuffer(b), [4,4,8], header=true)
    @test Data.header(tmp.schema)[1] == "abcd"
    @test Data.header(tmp.schema)[2] == "Column2"
    @test Data.header(tmp.schema)[3] == "10102017"
    @test Data.types(tmp.schema)[1] == Union{Missing, String}
    @test Data.types(tmp.schema)[2] == Union{Missing, String}
    @test Data.types(tmp.schema)[3] == Union{Missing, String}
    tmp = FWF.Source(file, [4,4,8], types=[String, Int, Float64])
    @test Data.types(tmp.schema)[1] == Union{Missing, String}
    @test Data.types(tmp.schema)[2] == Union{Missing, Int64}
    @test Data.types(tmp.schema)[3] == Union{Missing, Float64}
    tmp = FWF.Source(file, [4,4,8], types=[String, Int, DateFormat("mmddyyyy")], missings=["NA","***"])
    @test Data.types(tmp.schema)[1] == Union{Missing, String}
    @test Data.types(tmp.schema)[2] == Union{Missing, Int64}
    @test Data.types(tmp.schema)[3] == Union{Missing, Date}
    @test tmp.options.dateformats[3] == DateFormat("mmddyyyy")
    @test "***" in tmp.options.missingvals
    @test "NA" in tmp.options.missingvals
end

@testset "Functions" begin
    b = """
    abc
    def"""

    # Range building
    @test_throws ArgumentError FWF.calculate_ranges([-1,2,3])
    @test_throws ArgumentError FWF.calculate_ranges([-1:2])
    @test_throws ArgumentError FWF.calculate_ranges([1:2, 1:5])
    @test_throws ArgumentError FWF.calculate_ranges([1:5, 6:5])
    tmp = FWF.calculate_ranges([1,4,6])
    @test tmp[1] == 1:1
    @test tmp[2] == 2:5
    @test tmp[3] == 6:11
    tmp = FWF.calculate_ranges([1:1, 2:5, 6:11])
    @test tmp[1] == 1:1
    @test tmp[2] == 2:5
    @test tmp[3] == 6:11

end


text:
Pkg.init()
Pkg.add("Clustering")
Pkg.add("Color")
Pkg.add("DataFrames")
Pkg.add("DataStructures")
Pkg.add("Distributions")
Pkg.add("FixedPointNumbers")
Pkg.add("Gadfly")
Pkg.add("HDF5")
Pkg.add("IJulia")
Pkg.add("Images")
Pkg.add("JSON")
Pkg.add("MCMC")
Pkg.add("MeCab")
Pkg.add("ODBC")
Pkg.add("PyPlot")
Pkg.add("RDatasets")

ext: jl

max_stars_repo_name: mrkn/omni-jupyter

max_issues_repo_name: mrkn/omni-jupyter

max_forks_repo_name: mrkn/omni-jupyter


text:
# This file was generated, do not modify it. # hide
xgb  = XGBoostClassifier()
xgbm = machine(xgb, X, y)


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

function intersection(f, g::Geometry, p::Polygon)
  for t in triangulate(p)
    I = intersection(g, t)
    if type(I) != NoIntersection
      return I
    end
  end
  return @IT NoIntersection nothing f
end


text:
using ResultTypes
using Compat.Test
using Nullables

@testset "ResultTypes" begin

@testset "Result" begin
    @testset "Basic result" begin
        x = Result(2)
        @test unwrap(x) === 2
        @test_throws ErrorException unwrap_error(x)
        @test iserror(x) === false
    end

    @testset "Result with error type" begin
        x = Result(2, DivideError)
        @test unwrap(x) === 2
        @test_throws ErrorException unwrap_error(x)
        @test iserror(x) === false
    end

    @testset "Malformed result" begin
        x = Result(Nullable{Int}(), Nullable{DivideError}())
        @test_throws ErrorException unwrap(x)
        @test_throws ErrorException unwrap_error(x)
    end
end

@testset "ErrorResult" begin
    @testset "Basic error" begin
        x = ErrorResult(Int, "Basic Error")
        @test_throws ErrorException unwrap(x)
        e = unwrap_error(x)
        @test isa(e, ErrorException)
        @test e.msg == "Basic Error"
        @test iserror(x) === true
    end

    @testset "Empty error" begin
        x = ErrorResult(Int)
        @test_throws ErrorException unwrap(x)
        e = unwrap_error(x)
        @test isa(e, ErrorException)
        @test iserror(x) === true
    end

    @testset "Typed Error" begin
        x = ErrorResult(Int, DivideError())
        @test_throws DivideError unwrap(x)
        @test isa(unwrap_error(x), DivideError)
        @test iserror(x) === true
    end
end

@testset "Convert" begin
    @testset "From Result" begin
        x = Result(2)
        @test convert(Int, x) === 2
        @test convert(Float64, x) === 2.0
        @test_throws MethodError convert(String, x)
    end

    @testset "From ErrorResult" begin
        x = ErrorResult(Int, "Foo")
        @test_throws ErrorException convert(Int, x)
    end

    @testset "To Result" begin
        x = convert(Result{Int, ErrorException}, 2.0)
        @test unwrap(x) === 2
        @test isa(x, Result{Int, ErrorException})
    end

    @testset "To ErrorResult" begin
        x = convert(Result{Int, DivideError}, DivideError())
        @test unwrap_error(x) == DivideError()
        @test isa(x, Result{Int, DivideError})
    end
end

@testset "Return" begin
    let
        function foo(x::Int, y::Int)::Result{Int, DivideError}
            try
                return div(x, y)
            catch e
                return e
            end
        end

        x = foo(3, 2)
        @test isa(x, Result{Int, DivideError})
        @test unwrap(x) === 1

        y = foo(1, 0)
        @test isa(y, Result{Int, DivideError})
        @test unwrap_error(y) == DivideError()
    end
end

@testset "Show" begin
    @testset "Result" begin
        x = Result(2)

        @test isa(string(x), String)
    end

    @testset "ErrorResult" begin
        x = ErrorResult(Int, "Basic Error")
        @test isa(string(x), String)
    end
end

end


text:
using VisualCircuits, Compose

width, height = 7.5, 4
nand1 = Gate(5, 1, :NAND)
nand2 = Gate(5, 3, :NAND)
nand3 = Gate(3, nand1[:in1][2], :NAND, [:in1, :in2, :in3])
nand4 = Gate(3, nand2[:in2][2], :NAND, [:in1, :in2, :in3])

circuit = begin
    compose(
    # Range of coordinate system
    context(units=UnitBox(0, 0, width, height)),

    # Gates
    render(nand1), render(nand2), render(nand3), render(nand4),

    # Connect to edges (point, x_value)
    connect(nand1[:out], width-1),
    connect(nand2[:out], width-1),
    connect(nand3[:in1], 1),
    connect(nand4[:in3], 1),
    connect((1.5, 2), 1),

    # Connect Gates
    connect(nand1[:in1], nand3[:out]),
    connect(nand2[:in2], nand4[:out]),

    # crossing lines
    connect(
        nand1[:in2],
        (4, nand1[:in2][2]),
        (4, 1.75), (6, 2.25),
        (6, nand2[:out][2]),
        nand2[:out]
    ),
    connect(
        nand2[:in1],
        (4, nand2[:in1][2]),
        (4, 2.25), (6, 1.75),
        (6, nand1[:out][2]),
        nand1[:out]
    ),
    # Clock line
    connect(
        nand3[:in2],
        (1.5, nand3[:in2][2]),
        (1.5, nand4[:in2][2]),
        nand4[:in2]
    ),
    # looping lines
    connect(nand3[:in3], (2, nand3[:in3][2]), (2, 1.5), (4, 1.5)),
    connect(nand4[:in1], (2, nand4[:in1][2]), (2, 2.5), (4, 2.5)),

    # connections/dots
    dots(
        (6, nand1[:out][2]),
        (6, nand2[:out][2]),
        (1.5, 2),
        (4, 1.5),
        (4, 2.5)
    ),

    # render text
    render(
        (0.5, nand3[:in1][2], "J"),
        (0.5, nand4[:in3][2], "K"),
        (0.5, 2, "C"),
        (width-0.5, nand1[:out][2], "Q"),
        (width-0.5, nand2[:out][2], "~Q")
    ),

    # compose options
    stroke("black"), linewidth(2mm)
    )
end

# draw(PDF("JK_FlipFlop.pdf", autosize(width, height)...), circuit)
draw(SVG("JK_FlipFlop.svg", autosize(width, height)...), circuit)


text:
immutable HammingDistance <: Distance
end

"""
The hamming distance between two sequenced collections, or strings, of equal
length is measured as the number of differences between them at each index.
"""
hamming() = HammingDistance()

# Shouldn't be repeating code here.
# Need to find a way to unify these two method headers.
function dist(::HammingDistance, x::DirectIndexString, y::DirectIndexString)
  dst = 0
  for i in 1 : min(length(x), length(y))
    @inbounds if x[i] != y[i]
      dst += 1
    end
  end
  return dst
end

function dist{T}(::HammingDistance, x::Vector{T}, y::Vector{T})
  dst = 0
  for i in 1 : min(length(x), length(y))
    @inbounds if x[i] != y[i]
      dst += 1
    end
  end
  return dst
end


text:
"""
    crc32(data)

Compute the CRC-32 checksum over the `data` input. `data` can be
`BufferedInputStream` or `Vector{UInt8}`.
"""
function crc32 end

function crc32(stream::BufferedInputStream)
    crc = crc32()
    while !eof(stream)
        crc = crc32(crc, pointer(stream.buffer), stream.available)
        stream.position = 1
        stream.available = 0
        BufferedStreams.fillbuffer!(stream)
    end
    return crc::UInt32
end

function crc32(data::Vector{UInt8})
    return crc32(crc32(), pointer(data), length(data))::UInt32
end

"""
    adler32(data)

Compute the Adler-32 checksum over the `data` input. `data` can be
`BufferedInputStream` or `Vector{UInt8}`.
"""
function adler32 end

function adler32(stream::BufferedInputStream)
    adler = adler32()
    while !eof(stream)
        adler = adler32(adler, pointer(stream.buffer), stream.available)
        stream.position = 1
        stream.available = 0
        BufferedStreams.fillbuffer!(stream)
    end
    return adler::UInt32
end

function adler32(data::Vector{UInt8})
    return adler32(adler32(), pointer(data), length(data))::UInt32
end


text:
using StaticLint, SymbolServer
using CSTParser, Test
using StaticLint: scopeof, bindingof, refof, errorof, check_all

server = StaticLint.FileServer();

function get_ids(x, ids = [])
    if StaticLint.typof(x) == CSTParser.IDENTIFIER
        push!(ids, x)
    else
        for a in x
            get_ids(a, ids)
        end
    end
    ids
end

function parse_and_pass(s)
    empty!(server.files)
    f = StaticLint.File("", s, CSTParser.parse(s, true), nothing, server)
    StaticLint.setroot(f, f)
    StaticLint.setfile(server, "", f)
    StaticLint.scopepass(f)
    return f.cst
end

function check_resolved(s)
    cst = parse_and_pass(s)
    IDs = get_ids(cst)
    [(refof(i) !== nothing) for i in IDs]
end

@testset "StaticLint" begin

@testset "Basic bindings" begin

@test check_resolved("""
x
x = 1
x
""")  == [false, true, true]

@test check_resolved("""
x, y
x = y = 1
x, y
""")  == [false, false, true, true, true, true]

@test check_resolved("""
x, y
x, y = 1, 1
x, y
""")  == [false, false, true, true, true, true]

@test check_resolved("""
M
module M end
M
""")  == [false, true, true]

@test check_resolved("""
f
f() = 0
f
""")  == [false, true, true]

@test check_resolved("""
f
function f end
f
""")  == [false, true, true]

@test check_resolved("""
f
function f() end
f
""")  == [false, true, true]

@test check_resolved("""
function f(a)
end
""")  == [true, true]

@test check_resolved("""
f, a
function f(a)
    a
end
f, a
""")  == [false, false, true, true, true, true, false]

@test check_resolved("""
x
let x = 1
    x
end
x
""")  == [false, true, true, false]

@test check_resolved("""
x,y
let x = 1, y = 1
    x, y
end
x, y
""")  == [false, false, true, true, true, true, false, false]

@test check_resolved("""
function f(a...)
    f(a)
end
""")  == [true, true, true, true]

@test check_resolved("""
for i = 1:1
end
""")  == [true]

@test check_resolved("""
[i for i in 1:1]
""")  == [true, true]

@test check_resolved("""
[i for i in 1:1 if i]
""")  == [true, true, true]

# @test check_resolved("""
# @deprecate f(a) sin(a)
# f
# """)  == [true, true, true, true, true, true]

@test check_resolved("""
@deprecate f sin
f
""")  == [true, true, true, true]

@test check_resolved("""
module Mod
f = 1
end
using .Mod: f
f
""") == [true, true, true, true, true]

@test check_resolved("""
module Mod
module SubMod
    f() = 1
end
using .SubMod: f
f
end
""") == [true, true, true, true, true, true]

@test check_resolved("""
struct T
    field
end
function f(arg::T)
    arg.field
end
""") == [true, true, true, true, true, true, true]

@test check_resolved("""
f(arg) = arg
""") == [1, 1, 1]

@test check_resolved("-(r::T) where T = r") == [1, 1, 1, 1]
@test check_resolved("[k * j for j = 1:10 for k = 1:10]") == [1, 1, 1, 1]
@test check_resolved("[k * j for j in 1:10 for k in 1:10]") == [1, 1, 1, 1]

@testset "inference" begin
    @test bindingof(parse_and_pass("f(arg) = arg")[1]).type == StaticLint.CoreTypes.Function
    @test bindingof(parse_and_pass("function f end")[1]).type == StaticLint.CoreTypes.Function
    @test bindingof(parse_and_pass("struct T end")[1]).type == StaticLint.CoreTypes.DataType
    @test bindingof(parse_and_pass("mutable struct T end")[1]).type == StaticLint.CoreTypes.DataType
    @test bindingof(parse_and_pass("abstract type T end")[1]).type == StaticLint.CoreTypes.DataType
    @test bindingof(parse_and_pass("primitive type T 8 end")[1]).type == StaticLint.CoreTypes.DataType
    @test bindingof(parse_and_pass("x = 1")[1][1]).type == StaticLint.CoreTypes.Int
    @test bindingof(parse_and_pass("x = 1.0")[1][1]).type == StaticLint.CoreTypes.Float64
    @test bindingof(parse_and_pass("x = \"text\"")[1][1]).type == StaticLint.CoreTypes.String
    @test bindingof(parse_and_pass("module A end")[1]).type == StaticLint.CoreTypes.Module
    @test bindingof(parse_and_pass("baremodule A end")[1]).type == StaticLint.CoreTypes.Module

    # @test parse_and_pass("function f(x::Int) x end")[1][2][3].binding.t == StaticLint.getsymbolserver(server)["Core"].vals["Function"]
    let cst = parse_and_pass("""
        struct T end
        function f(x::T) x end""")
        @test bindingof(cst[1]).type == StaticLint.CoreTypes.DataType
        @test bindingof(cst[2]).type == StaticLint.CoreTypes.Function
        @test bindingof(cst[2][2][3]).type == bindingof(cst[1])
        @test refof(cst[2][3][1]) == bindingof(cst[2][2][3])
    end
    let cst = parse_and_pass("""
        struct T end
        T() = 1
        function f(x::T) x end""")
        @test bindingof(cst[1]).type == StaticLint.CoreTypes.DataType
        @test bindingof(cst[3]).type == StaticLint.CoreTypes.Function
        @test bindingof(cst[3][2][3]).type == bindingof(cst[1])
        @test refof(cst[3][3][1]) == bindingof(cst[3][2][3])
    end

    let cst = parse_and_pass("""
        struct T end
        t = T()""")
        @test bindingof(cst[1]).type == StaticLint.CoreTypes.DataType
        @test bindingof(cst[2][1]).type == bindingof(cst[1])
    end

    let cst = parse_and_pass("""
        module A
        module B
        x = 1
        end
        module C
        import ..B
        B.x
        end
        end""")
        @test refof(cst[1][3][2][3][2][3][1]) == bindingof(cst[1][3][1][3][1][1])
    end

    let cst = parse_and_pass("""
        struct T0
            x
        end
        struct T1
            field::T0
        end
        function f(arg::T1)
            arg.field.x
        end""");
        @test refof(cst[3][3][1][1][1]) == bindingof(cst[3][2][3])
        @test refof(cst[3][3][1][1][3][1]) == bindingof(cst[2][3][1])
        @test refof(cst[3][3][1][3][1]) == bindingof(cst[1][3][1])
    end

    let cst = parse_and_pass("""
        raw"whatever"
        """)
        @test refof(cst[1][1]) !== nothing
    end
    let cst = parse_and_pass("""
        macro mac_str() end
        mac"whatever"
        """)
        @test refof(cst[2][1]) == bindingof(cst[1])
    end

    let cst = parse_and_pass("""
        [i * j for i = 1:10 for j = i:10]
        """)
        @test refof(cst[1][2][1][3][3][1]) == bindingof(cst[1][2][1][1][3][1])
    end
    let cst = parse_and_pass("""
        [i * j for i = 1:10, j = 1:10 for k = i:10]
        """)
        @test refof(cst[1][2][1][3][3][1]) == bindingof(cst[1][2][1][1][3][1])
    end

    let cst = parse_and_pass("""
        module Reparse
        end
        using .Reparse, CSTParser
        """)
        @test refof(cst[2][3]).val == bindingof(cst[1])
    end

    let cst = parse_and_pass("""
        module A
        A
        end
        """)
        @test scopeof(cst).names["A"] == scopeof(cst[1]).names["A"]
        @test refof(cst[1][2]) == bindingof(cst[1])
        @test refof(cst[1][3][1]) == bindingof(cst[1])
    end
    # let cst = parse_and_pass("""
    #     using Test: @test
    #     """)
    #     @test bindingof(cst[1][4]) !== nothing
    # end
    let cst = parse_and_pass("""
        sin(1,2,3)
        """)
        check_all(cst, StaticLint.LintOptions(),server)
        @test errorof(cst[1]) === StaticLint.IncorrectCallArgs
    end
    let cst = parse_and_pass("""
        for i in length(1) end
        for i in 1.1 end
        for i in 1 end
        for i in 1:1 end
        """)
        check_all(cst, StaticLint.LintOptions(),server)
        @test errorof(cst[1][2]) === StaticLint.IncorrectIterSpec
        @test errorof(cst[2][2]) === StaticLint.IncorrectIterSpec
        @test errorof(cst[3][2]) === StaticLint.IncorrectIterSpec
        @test errorof(cst[4][2]) === nothing
    end

    let cst = parse_and_pass("""
        [i for i in length(1) end]
        [i for i in 1.1 end]
        [i for i in 1 end]
        [i for i in 1:1 end]
        """)
        check_all(cst, StaticLint.LintOptions(),server)
        @test errorof(cst[1][2][3]) === StaticLint.IncorrectIterSpec
        @test errorof(cst[2][2][3]) === StaticLint.IncorrectIterSpec
        @test errorof(cst[3][2][3]) === StaticLint.IncorrectIterSpec
        @test errorof(cst[4][2][3]) === nothing
    end

    let cst = parse_and_pass("a == nothing")
        check_all(cst, StaticLint.LintOptions(),server)
        @test errorof(cst[1][2]) === StaticLint.NothingEquality
    end

    let cst = parse_and_pass("""
        struct Graph
            children:: T
        end

        function test()
            g = Graph()
            f = g.children
        end""")
        @test cst[2][3][2][3][3][1] in bindingof(cst[1][3][1]).refs
    end

    let cst = parse_and_pass("""
        __source__
        __module__
        macro m()
            __source__
            __module__
        end""")
        @test refof(cst[1]) === nothing
        @test refof(cst[2]) === nothing
        @test refof(cst[3][3][1]) !== nothing
        @test refof(cst[3][3][2]) !== nothing
    end
end

@testset "macros" begin
    @test check_resolved("""
    @enum(E,a,b)
    E
    a
    b
    """)  == [true, true, true, true, true, true, true]
    end

    @test check_resolved("""
    @enum E a b
    E
    a
    b
    """)  == [true, true, true, true, true, true, true]

    @test check_resolved("""
    @enum E begin
        a
        b
    end
    E
    a
    b
    """)  == [true, true, true, true, true, true, true]
end

@testset "tuple args" begin
    let cst = parse_and_pass("""
        function f((arg1, arg2))
            arg1, arg2
        end""")
        @test StaticLint.hasref(cst[1][3][1][1])
        @test StaticLint.hasref(cst[1][3][1][3])
    end

    let cst = parse_and_pass("""
        function f((arg1, arg2) = (1,2))
            arg1, arg2
        end""")
        @test StaticLint.hasref(cst[1][3][1][1])
        @test StaticLint.hasref(cst[1][3][1][3])
    end

    let cst = parse_and_pass("""
        function f((arg1, arg2)::Tuple{Int,Int})
            arg1, arg2
        end""")
        @test StaticLint.hasref(cst[1][3][1][1])
        @test StaticLint.hasref(cst[1][3][1][3])
    end
end

@testset "type params check" begin
    let cst = parse_and_pass("""
        f() where T
        f() where {T,S}
        f() where {T<:Any}
        """)
        StaticLint.check_all(cst, StaticLint.LintOptions(), server)
        @test StaticLint.errorof(cst[1][3]) == StaticLint.UnusedTypeParameter
        @test StaticLint.errorof(cst[2][4]) == StaticLint.UnusedTypeParameter
        @test StaticLint.errorof(cst[2][6]) == StaticLint.UnusedTypeParameter
        @test StaticLint.errorof(cst[3][4]) == StaticLint.UnusedTypeParameter
    end
    let cst = parse_and_pass("""
        f(x::T) where T
        f(x::T,y::S) where {T,S}
        f(x::T) where {T<:Any}
        """)
        StaticLint.check_all(cst, StaticLint.LintOptions(), server)
        @test !StaticLint.haserror(cst[1][3])
        @test !StaticLint.haserror(cst[2][4])
        @test !StaticLint.haserror(cst[2][6])
        @test !StaticLint.haserror(cst[3][4])
    end
end

@testset "overwrites_imported_function" begin
    let cst = parse_and_pass("""
        import Base:sin
        using Base:cos
        sin(x) = 1
        cos(x) = 1
        Base.tan(x) = 1
        """)

        @test StaticLint.overwrites_imported_function(bindingof(cst[3]))
        @test !StaticLint.overwrites_imported_function(bindingof(cst[4]))
        @test StaticLint.overwrites_imported_function(bindingof(cst[5]))
    end
end

@testset "pirates" begin
    let cst = parse_and_pass("""
        import Base:sin
        struct T end
        sin(x::Int) = 1
        sin(x::T) = 1
        sin(x::Array{T}) = 1
        """)
        StaticLint.check_for_pirates(cst[3])
        StaticLint.check_for_pirates(cst[4])
        @test errorof(cst[3]) === StaticLint.TypePiracy
        @test errorof(cst[4]) === nothing
    end
    let cst = parse_and_pass("""
        struct AreaIterator{T}
            array::AbstractMatrix{T}
            radius::Int
        end
        Base.eltype(::Type{AreaIterator{T}}) where T = Tuple{T, AbstractVector{T}}
        """)
        StaticLint.check_for_pirates(cst[2])
        @test errorof(cst[2]) === nothing
    end
    let cst = parse_and_pass("""
        import Base:sin
        abstract type T end
        sin(x::Array{T}) = 1
        sin(x::Array{<:T}) = 1
        sin(x::Array{Number}) = 1
        sin(x::Array{<:Number}) = 1
        """)
        StaticLint.check_for_pirates(cst[3])
        StaticLint.check_for_pirates(cst[4])
        StaticLint.check_for_pirates(cst[5])
        StaticLint.check_for_pirates(cst[6])
        @test errorof(cst[3]) === nothing
        @test errorof(cst[4]) === nothing
        @test errorof(cst[5]) === StaticLint.TypePiracy
        @test errorof(cst[6]) === StaticLint.TypePiracy
    end
end

@testset "docs for undescribed variables" begin
let cst = parse_and_pass("""
    \"\"\"
        somefunc() = true
    \"\"\"
    somefunc
    somefunc() = true
    """)
    @test StaticLint.hasref(cst[1][3])
    @test StaticLint.hasbinding(cst[1][3])
    @test refof(cst[1][3]) == bindingof(cst[1][3])
end
end

@testset "check_call" begin
    let cst = parse_and_pass("""
        sin(1)
        sin(1,2)
        """)
        StaticLint.check_call(cst[1], server)
        StaticLint.check_call(cst[2], server)
        @test StaticLint.errorof(cst[1]) === nothing
        @test StaticLint.errorof(cst[2]) == StaticLint.IncorrectCallArgs
    end

    let cst = parse_and_pass("""
        Base.sin(a,b) = 1
        function Base.sin(a,b)
            1
        end
        """)
        StaticLint.check_call(cst[1][1], server)
        @test StaticLint.errorof(cst[1][1]) === nothing
        StaticLint.check_call(cst[2][2], server)
        @test StaticLint.errorof(cst[2][2]) === nothing
    end

    let cst = parse_and_pass("""
        f(x) = 1
        f(1, 2)
        """)
        StaticLint.check_call(cst[2], server)
        @test StaticLint.errorof(cst[2]) === StaticLint.IncorrectCallArgs
    end

    let cst = parse_and_pass("""
        view([1], 1, 2, 3)
        """)
        StaticLint.check_call(cst[1], server)
        @test StaticLint.errorof(cst[1]) === nothing
    end

    let cst = parse_and_pass("""
        f(a...) = 1
        f(1)
        f(1, 2)
        """)
        StaticLint.check_call(cst[2], server)
        StaticLint.check_call(cst[3], server)
        @test StaticLint.errorof(cst[2]) === nothing
        @test StaticLint.errorof(cst[3]) === nothing
    end
    let cst = parse_and_pass("""
        function func(a, b)
            func(a...)
        end
        """)
        StaticLint.check_call(cst[1][3][1], server)
        m_counts = StaticLint.func_nargs(cst[1])
        call_counts = StaticLint.call_nargs(cst[1][3][1])
        @test StaticLint.errorof(cst[1][3][1]) === nothing
    end
    let cst = parse_and_pass("""
        function func(@nospecialize args...) end
        func(1, 2)
        """)
        @test StaticLint.func_nargs(cst[1]) == (0, typemax(Int), String[], false)
        StaticLint.check_call(cst[2], server)
        @test StaticLint.errorof(cst[2]) === nothing
    end
    let cst = parse_and_pass("""
        argtail(x, rest...) = 1
        tail(x::Tuple) = argtail(x...)
        """)
        @test StaticLint.func_nargs(cst[1]) == (1, typemax(Int), String[], false)
        StaticLint.check_call(cst[2], server)
        @test StaticLint.errorof(cst[2]) === nothing
    end
    let cst = parse_and_pass("""
        func(arg::Vararg{T,N}) where N = arg
        func(a,b)
        """)

        @test StaticLint.func_nargs(cst[1]) == (0, typemax(Int), String[], false)
        StaticLint.check_call(cst[2], server)
        @test StaticLint.errorof(cst[2]) === nothing
    end
    let cst = parse_and_pass("""
        function f(a, b; kw = kw) end
        f(1,2, kw = 1)
        """)
        StaticLint.check_call(cst[2], server)
        @test StaticLint.errorof(cst[2]) === nothing
    end
    let cst = parse_and_pass("""
        func(a,b,c,d) = 1
        func(a..., 2)
        """)
        StaticLint.call_nargs(cst[2])
        StaticLint.check_call(cst[2], server)
        @test StaticLint.errorof(cst[2]) === nothing
    end
    let cst = parse_and_pass("""
        @kwdef struct A
            x::Float64
        end
        A(x = 5.0)
        """)
        StaticLint.check_call(cst[2], server)
        @test StaticLint.errorof(cst[2]) === nothing
    end
end

@testset "check_modulename" begin
    let cst = parse_and_pass("""
        module Mod1
        module Mod11
        end
        end
        module Mod2
        module Mod2
        end
        end
        """)
        StaticLint.check_modulename(cst[1])
        StaticLint.check_modulename(cst[1][3][1])
        StaticLint.check_modulename(cst[2])
        StaticLint.check_modulename(cst[2][3][1])

        @test StaticLint.errorof(cst[1][2]) === nothing
        @test StaticLint.errorof(cst[1][3][1][2]) === nothing
        @test StaticLint.errorof(cst[2][2]) === nothing
        @test StaticLint.errorof(cst[2][3][1][2]) === StaticLint.InvalidModuleName
    end
end

if !(VERSION < v"1.3")
@testset "non-std var syntax" begin
    let cst = parse_and_pass("""
        var"name" = 1
        var"func"(arg) = arg
        function var"func1"() end
        name
        func
        func1
        """)
        StaticLint.collect_hints(cst)
        @test all(n in keys(cst.meta.scope.names) for n in ("name", "func"))
        @test StaticLint.hasref(cst[4])
        @test StaticLint.hasref(cst[5])
        @test StaticLint.hasref(cst[6])
    end
end
end

if false # Not to be run, requires JuMP
    @testset "JuMP macros" begin
let cst = parse_and_pass("""
    using JuMP
    model = Model()
    some_bound = 1
    @variable(model, x0)
    @variable(model, x1, somekw=1)
    @variable(model, x2 <= 1)
    @variable(model, x3 >= 1)
    @variable(model, 1 <= x4)
    @variable(model, 1 >= x5)
    @variable(model, x6 >= some_bound)
    # @variable(model, some_bound >= x7)
    """)
    @test isempty(StaticLint.collect_hints(cst))
end

let cst = parse_and_pass("""
    using JuMP
    model = Model()
    some_bound = 1
    @variable model x0
    @variable model x1 somekw=1
    @variable model x2 <= 1
    @variable model x3 >= 1
    @variable model 1 <= x4
    @variable model 1 >= x5
    @variable model x6 >= some_bound
    # @variable(model, some_bound >= x7)
    """)
    @test isempty(StaticLint.collect_hints(cst))
end

let cst = parse_and_pass("""
    using JuMP
    model = Model()
    some_bound = 1
    @variable(model, some_bound >= x7)
    """)
    @test !StaticLint.hasref(cst[4][5][3])
end

let cst = parse_and_pass("""
    using JuMP
    model = Model()
    some_bound = 1
    @expression(model, ex, some_bound >= 1)
    """)
    @test isempty(StaticLint.collect_hints(cst))
end

let cst = parse_and_pass("""
    using JuMP
    model = Model()
    @expression(model, expr, 1 == 1)
    @constraint(model, con1, expr)
    @constraint model con2 expr
    """)
    @test isempty(StaticLint.collect_hints(cst))
end
end
end

@testset "stdcall" begin
    let cst = parse_and_pass("""
        ccall(:GetCurrentProcess, stdcall, Ptr{Cvoid}, ())""")
        StaticLint.collect_hints(cst)
        @test isempty(StaticLint.collect_hints(cst))
    end
    let cst = parse_and_pass("""
        stdcall
        """)
        @test !StaticLint.hasref(cst[1])
    end
end

@testset "check_if_conds" begin
    let cst = parse_and_pass("""
        if true end
        """)
        StaticLint.check_if_conds(cst[1])
        @test cst[1][2].meta.error == StaticLint.ConstIfCondition
    end
    let cst = parse_and_pass("""
        if x = 1 end
        """)
        StaticLint.check_if_conds(cst[1])
        @test cst[1][2].meta.error == StaticLint.EqInIfConditional
    end
    let cst = parse_and_pass("""
        if a || x = 1 end
        """)
        StaticLint.check_if_conds(cst[1])
        @test cst[1][2].meta.error == StaticLint.EqInIfConditional
    end
    let cst = parse_and_pass("""
        if x = 1 && b end
        """)
        StaticLint.check_if_conds(cst[1])
        @test cst[1][2].meta.error == StaticLint.EqInIfConditional
    end
end

@testset "check redefinition of const" begin
    let cst = parse_and_pass("""
        T = 1
        struct T end
        """)
        StaticLint.check_const_decl(cst[2])
        @test cst[2].meta.error == StaticLint.CannotDeclareConst
    end
    let cst = parse_and_pass("""
        struct T end
        T = 1
        """)
        StaticLint.check_const_redef(cst[2][1])
        @test cst[2][1].meta.error == StaticLint.InvalidRedefofConst
    end
    let cst = parse_and_pass("""
        struct T end
        T() = 1
        """)
        StaticLint.check_const_redef(cst[2])
        @test cst[2].meta.error == nothing
    end
end

@testset "hoisting of inner constructors" begin
    let cst = parse_and_pass("""
        struct ASDF
            x::Int
            y::Int
            ASDF(x::Int) = new(x, 1)
        end
        ASDF() = something
        """)
        @test bindingof(cst[1]) === bindingof(cst[1][3][3]).prev
        @test bindingof(cst[1][3][3]) === bindingof(cst[2]).prev
    end
end

@testset "using of self" begin # e.g. `using StaticLint: StaticLint`
    let cst = parse_and_pass("""
        using Base.Filesystem: Filesystem
        """)
        @test StaticLint.hasref(cst[1][6])
    end
end
end


text:
module MFactor
#by Daniel

using Nemo
using Hecke

mutable struct disolveinfo{E}
    zero                ::E                         # the zero polynomial
    xalphas             ::Vector{E}                 # [x_i - alpha_i for i = 1..r]
    betas               ::Vector{Vector{E}}         # matrix of beta evaluations
    inv_prod_betas1     ::Vector{E}                 # info for inivariate solver
    prod_betas_coeffs   ::Vector{Vector{Vector{E}}} # matrix of taylor coeffs
    deltas              ::Vector{Vector{E}}         # answers
    delta_coeffs        ::Vector{Vector{Vector{E}}} # taylor coeffs of answer
end

############### subsets of 1, ..., n of size m ################################

mutable struct msubsets
    n::Int
    m::Int
end

function Base.iterate(a::msubsets)
    state = [i <= a.m ? true : false for i in 1:a.n]
    return (copy(state), state)
end

function Base.iterate(a::msubsets, state)
    i = 0
    while i < a.n && !state[1+i]
        i += 1
    end
    t1 = i
    while i < a.n && state[1+i]
        i += 1
    end
    t2 = i
    if t2 >= a.n
        return nothing
    end
    newstate = [t2 <= i < a.n ? state[1+i] : false for i in 0:a.n-1]
    newstate[1+t2] = true
    for i in 0:t2-t1-2
        newstate[1+i] = true
    end
    (copy(newstate), newstate)
end

# interface for enumerating evaluation points
# probably over-engineered and a simple random will do

function tuple_init(n::Int)
    return zeros(fmpz, n)
end

function tuple_saturate!(alpha::Vector{fmpz}, m::Int)
    n = length(alpha)
    @assert n > 0
    for i in m+1:n-1
        alpha[1 + m] += alpha[1 + i]
        alpha[1 + i] = 0
    end
    if m < n && alpha[1 + m] == 0
        for i in 1:m
            if alpha[i] != 0
                return
            end
        end
        alpha[1 + m] = 1
    end
end

function tuple_next!(alpha::Vector{fmpz})
    n = length(alpha)
    @assert n > 0
    sum = fmpz(0)
    for i in 1:n
        sum += alpha[i]
    end
    i = n - 1
    while i >= 0 && alpha[1 + i] == 0
        i -= 1
    end
    t1 = i
    while i >= 0 && alpha[1 + i] != sum
        i -= 1
    end
    t2 = i
    while i >= 0 && alpha[1 + i] == sum
        i -= 1
    end
    t3 = i
    if t1 > 0 && t1 != t2
        (alpha[1 + t1], alpha[1 + n - 1]) = (alpha[1 + n - 1], alpha[1 + t1])
        alpha[1 + n - 1] -= 1
        alpha[1 + t1 - 1] += 1
    elseif t1 > 0 && t1 == t2 && t3 >= 0
        alpha[1 + t3] += 1
        alpha[1 + t3 + 1] = 0
        alpha[1 + n - 1] = sum - 1
    else
        alpha[1 + n - 1] = alpha[1] + 1
        if n > 1
            alpha[1] = 0
        end
    end
end

########### various helpers for expansions in terms of (x - alpha) ############

# the evaluation of p at x = alpha. Don't use evaluate for this!!!
function taylor_rem(p::E, xalpha::E) where E
    q, r = divrem(p, xalpha)
    return r
end

function taylor_divrem(p::E, xalpha::E) where E
    return divrem(p, xalpha)
end

# set t to the coefficients of q when expanded in powers of (x - alpha)
function taylor_get_coeffs!(t::Vector{E}, q::E, xalpha::E) where E
    R = parent(q)
    empty!(t)
    while !iszero(q)
        q, r = divrem(q, xalpha)
        push!(t, r)
    end
end

function taylor_get_coeffs(q::E, xalpha::E) where E
    t = E[]
    taylor_get_coeffs(t, q, xalpha)
    return t
end

# opposite of taylor_get_coeffs
function taylor_from_coeffs(B::Vector{E}, xalpha::E, zero::E) where E
    a = zero
    for i in length(B):-1:1
        a = a*xalpha + B[i]
    end
    return a
end

function taylor_set_coeff!(t::Vector{E}, i::Int, c::E, zero::E) where E
    while 1 + i > length(t)
        push!(t, zero)
    end
    t[1 + i] = c
    while length(t) > 0 && iszero(t[end])
        popback!(t)
    end
end

#################### basic manipulation of Fac ################################

# a *= b^e
function mulpow!(a::Fac{T}, b::T, e::Int) where T
    @assert e >= 0
    if e < 1
        return
    end
    if isconstant(b)
        a.unit *= b^e
    elseif haskey(a.fac, b)
        a.fac[b] += e
    else
        a.fac[b] = e
    end
end

function mulpow!(a::Fac{T}, b::Fac{T}, e::Int) where T
    @assert !(a === b)
    @assert e >= 0
    if e < 1
        return
    end
    a.unit *= b.unit^e
    if isempty(a.fac)
        a.fac = b.fac
    elseif !isempty(b.fac)
        a.fac = merge((A, B) -> A + B*e, a.fac, b.fac)
    end
end

################ basic manipulation of polys #################################

function gcdcofactors(a, b)
    g = gcd(a, b)
    if iszero(g)
        @assert iszero(a)
        @assert iszero(b)
        return (g, a, b)
    else
        return (g, divexact(a, g), divexact(b, g))
    end
end

# remove the content with repect to variable v
function primitive_part(a::E, v::Int)::E where E
    R = parent(a)
    d = degree(a, v)
    g = R(0)
    for i in 0:d
        ci = coeff(a, Int[v], Int[i])
        if !iszero(ci)
            g = gcd(g, ci)::E
        end
    end
    if iszero(g)
        return a
    elseif isone(g)
        return a
    else
        return divexact(a, g)
    end
end

function get_lc(a::E, v::Int)::E where E
    d = degree(a, v)
    @assert d >= 0
    return coeff(a, Int[v], Int[d])
end

function set_lc(a::E, v::Int, b::E) where E
    R = parent(a)
    d = degree(a, v)
    @assert d >= 0
    diff = (b - coeff(a, Int[v], Int[d]))
    if iszero(diff)
        return a
    else
        return a + diff*gen(R, v)^d
    end
end

function make_monic(a::E)::E where E
    if length(a) < 1
        return a
    end
    lc = coeff(a, 1)
    if isone(a)
        return a
    else
        return inv(lc) * a
    end
end

function eval_one(a::E, v, alpha)::E where E
    R = parent(a)
    return divrem(a, gen(R, v) - alpha)[2]
end

################ generic diophantine solver ###################################

#=
    R is a multivariate ring K[X, x_1, ..., x_n]
    We will be dealing with polynomials in the first l + 1 variables, i.e.
    only the variables X, x_1, ..., x_l appear.
    l is length of alpha.

    For fixed B[1], ..., B[r], precompute info for solving

        t/(B[1]*...*B[r]) = delta[1]/B[1] + ... + delta[r]/B[r]

    for delta[1], ..., delta[r] given t in X, x_1, ..., x_l.

    alpha is an array of evaluation points
        x_1 = alpha[1]
        ...
        x_l = alpha[l]

    If, after evaluation at all x_j = alpha[j], the B[i] did not drop degree
    in X and are pairwise relatively prime in K[X], then the function returns
    true and a suitable disolveinfo. Otherwise, it returns false and junk

    Upon success, we are ready to solve for the delta[i] using disolve.
    The only divisions that disolve does in K are by the leading coefficients
    of the inv_prod_betas1 member of disolveinfo
=#

function disolveinit(B::Vector{E}, alpha::Vector, R) where E
    K = base_ring(R)
    n = nvars(R) - 1
    l = length(alpha)   # number of evaluated variables
    r = length(B)       # number of factors

    @assert n > 0
    @assert r > 1

    # betas is a 1+l by r matrix with the original B[i] in the last row
    # and successive evaluations in the previous rows
    betas = [E[R(0) for i in 1:r] for j in 0:l]

    # prod_betas is like betas, but has B[1]*...*B[r]/B[i] in the i^th column
    # we don't need prod_betas, but rather
    #   prod_betas_coeffs = taylor of prod_betas in x_(1+j) - alpha[j]
    prod_betas_coeffs = [[E[] for i in 1:r] for j in 0:l]

    # inv_prod_betas1[i] will hold invmod(B[1]*...*B[r]/B[i], B[i]) in R[X]
    # after all x_j = alpha[j] have been evaluated away
    inv_prod_betas1 = E[R(0) for i in 1:r]

    # working space that need not be reallocated on each recursive call
    deltas = [E[R(0) for i in 1:r] for j in 0:l]
    delta_coeffs = [[E[] for i in 1:r] for j in 0:l]

    xalphas = E[gen(R, 1 + j) - R(alpha[j]) for j in 1:l]

    I = disolveinfo{E}(R(0), xalphas,
                       betas, inv_prod_betas1, prod_betas_coeffs,
                       deltas, delta_coeffs)

    for i in 1:r
        I.betas[1+l][i] = deepcopy(B[i])
    end
    for j in l:-1:1
        for i in 1:r
            I.betas[j][i] = taylor_rem(I.betas[j + 1][i], xalphas[j])
        end
    end

    for i in 1:r
        if degree(I.betas[1][i], 1) != degree(I.betas[l + 1][i], 1)
            # degree drop in invariates in X
            return false, I
        end
    end

    # univariate ring for gcdx
    S, x = PolynomialRing(K, "x")
    sub = [k == 1 ? x : S(1) for k in 1:n+1]

    for j in 0:l
        for i in 1:r
            p = R(1)
            for k in 1:r
                if k != i
                    p *= I.betas[1 + j][k]
                end
            end
            if j > 0
                taylor_get_coeffs!(I.prod_betas_coeffs[1 + j][i], p, I.xalphas[j])
            else
                s = evaluate(betas[1][i], sub)
                t = evaluate(p, sub)
                g, s1, t1 = gcdx(s, t)
                if !isconstant(g)
                    # univariates are not pairwise prime
                    return false, I
                end
                I.inv_prod_betas1[i] = evaluate(divexact(t1, g), gen(R,1))
                # leave prod_betas_coeffs[1][i] as zero, it doesn't matter
            end
        end
    end

    return true, I
end

# The second return of disolve (the answer) is owned by I.
# So, at least its length should not be mutated.
function disolve(I::disolveinfo{E}, t::E, degs::Vector{Int}) where E
    return disolve_recursive(I, t, degs, length(I.xalphas))
end

function disolve_recursive(I::disolveinfo{E}, t::E, degs::Vector{Int},
                                                              lev::Int) where E
    @assert 0 <= lev <= length(I.xalphas)
    @assert lev <= length(degs)

    r = length(I.inv_prod_betas1)

    deltas = I.deltas[1 + lev]

    if lev == 0
        for i in 1:r
            Q, deltas[i] = divrem(t*I.inv_prod_betas1[i], I.betas[1][i])
            # TODO check the denominator of deltas[i] and possibly return false
        end
        return true, deltas
    end

    newdeltas = I.deltas[lev]
    delta_coeffs = I.delta_coeffs[1 + lev]

    for i in 1:r
        empty!(delta_coeffs[i])
    end

    for l in 0:degs[lev]
        t, newt = divrem(t, I.xalphas[lev])
        for s in 0:l-1
            for i in 1:r
                if s < length(delta_coeffs[i]) &&
                   l - s < length(I.prod_betas_coeffs[1 + lev][i])
                    newt -= delta_coeffs[i][1 + s] *
                            I.prod_betas_coeffs[1 + lev][i][1 + l - s]
                end
            end
        end

        if iszero(newt)
            continue
        end

        ok, newdeltas = disolve_recursive(I, newt, degs, lev - 1)
        if !ok
            return false, deltas
        end

        for i in 1:r
            if iszero(newdeltas[i])
                continue
            end
            if false && l + length(I.prod_betas_coeffs[1 + lev][i]) - 1 > degs[lev]
                return false, deltas
            end
            taylor_set_coeff!(delta_coeffs[i], l, newdeltas[i], I.zero)
        end
    end

    for i in 1:r
        deltas[i] = taylor_from_coeffs(delta_coeffs[i], I.xalphas[lev], I.zero)
    end

    return true, deltas
end

########### extremely thorough tests #####################

function map_down(Rp, a, mKp :: Map, pr::Int)
    M = MPolyBuildCtx(Rp)
    pk = prime(base_ring(a))^pr
    for (c, v) in zip(coeffs(a), exponent_vectors(a))
        @assert valuation(c) >= pr
        q = divexact(c, pk)  #should be a shift
        push_term!(M, mKp(q), v)
    end
    return finish(M)
end

function map_up(R, a, mKp :: Map, pr::Int)
    Rp = parent(a)
    Kp = base_ring(Rp)
    M = MPolyBuildCtx(R)
    pk = prime(base_ring(R))^pr

    for (c, v) in zip(coeffs(a), exponent_vectors(a))
        d = preimage(mKp, c)
        push_term!(M, d*pk, v)
    end
    return finish(M)
end

function lift_prime_power(a::fmpq_mpoly, fac::Vector{Generic.MPoly{qadic}},
                           alphas::Vector, degs::Vector{Int}, k::Int)

    r = length(fac)
    R = parent(fac[1])
    n = nvars(R)
    ZZ = base_ring(R)
    Kp, mKp = ResidueField(ZZ)
    Rp, x = PolynomialRing(Kp, n)

    pr = 0

    ok, I = disolveinit([map_down(Rp, f, mKp, pr) for f in fac], alphas, Rp)
    @assert ok

    a = map_coeffs(ZZ, a, parent = parent(fac[1]))

    for l in 1:k
        pr += 1
        pro = prod(fac)
        error = a - pro
#println("error: ", error)
        if iszero(error)
          @show "exact"
            break
        end
        t = map_down(Rp, error, mKp, pr)
        ok, deltas = disolve(I, t, degs)
        @assert ok
        for i in 1:r
            fac[i] += map_up(R, deltas[i], mKp, pr)
        end
#println("newfac: ", fac)
    end
    return fac
end

end

#=
println("------------- lifing to 5^10  --------------------")
Kxyz, (x, y, z) = PolynomialRing(ZZ, ["x", "y", "z"])
lift_prime_power(
    (111111*x+y+z+555555)*(x-222222*y+z)*(x+y+333333*z),
    [(111111*x+y+z), (x-2*y+z), (x+y+3*z)], # all the lc_x must be correct :-(
    [1, 2], # evaluation points for y, z
    [3, 3], # degree bounds for y, z
    5, 10)

println("----------- bad leading coefficient --------------")
Kxyz, (x, y, z) = PolynomialRing(QQ, ["x", "y", "z"])
@time F = mfactor_char_zero(x^24-y^24*z^12)
println("yay! ", F)

println("-------------- extraneous factors ----------------")
Kxyz, (x, y, z) = PolynomialRing(QQ, ["x", "y", "z"])
@time F = mfactor_char_zero(((x^2+y^2+z^2+2)*(x+1)*(y+2)*(z+3) + x*y*z)*(x^2+y^2+z^2+3))
println("yay! ", F)

println("--------------- number field ---------------------")
QA, A = PolynomialRing(QQ, "A")
K, a = number_field(A^3 - 2, "a")
Kxyz, (x, y, z) = PolynomialRing(K, ["x", "y", "z"])
println("yay! ", mfactor_char_zero(Kxyz(0)))
println("yay! ", mfactor_char_zero(Kxyz(1)))
println("yay! ", mfactor_char_zero(Kxyz(2)))
println("yay! ", mfactor_char_zero(2*x^2*y*(1+y)*(x^3-2)))
println("yay! ", mfactor_char_zero(2*x^2*y*(1+y)*(2*x^3-1)*(a*x+a^2*y+1)*(a^2*x+a*y+1)^2))
println("yay! ", mfactor_char_zero(x^6-2*y^3))
println("yay! ", mfactor_char_zero((x^2+z-1)*(x+a*z^3)))
println("yay! ", mfactor_char_zero((x+y+a)*(x+a*y+1)*(x+y+2*a)))
println("yay! ", mfactor_char_zero((x^2+y^2+a)*(x^3+a*y^3+1)*(x^2+y+2*a)))
println("yay! ", mfactor_char_zero((1+x+a*y)^3*(1+x+a^2*y)^5))
println("yay! ", mfactor_char_zero((x+y+z)*(a*x+y+z)*(x+a*y+z)*(x+y+a*z)))

println("------------ rational benchmarketing -------------")
function number_of_odd_divisors(n::Int)
    count = 0
    for i in 1:2:n
        if n%i == 0
            count += 1
        end
    end
    return count
end
Kxyzt, (x, y, z, t) = PolynomialRing(QQ, ["x", "y", "z", "t"])
for i in 1:15
    @time P = ((1+x+y+z+t)^i+1)*((1+x+y+z+t)^i+2)
    @time F = mfactor_char_zero(P)
    println("power: ", i, ", number of factors: ", length(F.fac), ", expected: ", 1 + number_of_odd_divisors(i))
end

=#


text:
module Model

using Parameters

"""Hill function without basal activation"""
@with_kw struct hill
    k
    n
end

function (par::hill)(s)
    @unpack k,n = par
    return s^n / (k^n + s^n )
end

"""Hill function with basal activation"""
@with_kw struct hill_basal
    k
    n
    x_basal
    x_max
end
function (par::hill_basal)(s)
    @unpack k, n, x_base, x_max = par

    return x_base + (x_max - x_base)*(s^n/(k^n+s^n))
end

"""Model: Positive Feedback"""
@with_kw mutable struct positive_feedback
    x= 0.0
    s=0.1
    kK=1.
    kA=1.3
    m=3.
    x_50=0.5
    x_tot=1.0
    kp=1.0

    #@assert  sum(x) <= x_tot
end

function (p::positive_feedback)(x,p_,t)
    @unpack s, kK, kA, m, x_50, x_tot, kp = p
    dx = (kK * s + kA * (x^m/(x_50^m + x^m))) * (x_tot-x) - kp * x

    return dx
end

end

#=

@parameters t s kK kA m x_50 x_tot kp
@variables x(t)
@derivatives D'~t
=#


text:
import Pkg; Pkg.activate("env")

using StringEncodings

include("momo.jl")
include("h.jl")
include("control.jl")

USAGE = """
unzip dela-fr-public.zip # il y a des personnes louches dedans

import Pkg; Pkg.activate("env")
Pkg.add("StringEncodings")

include("haikusaga.jl")
haikusaga("aimer fait mal") #HS0201
haikusaga("plus que hair") #HS0202
s = "Liure au pieu ?
Quelque ruse salique ?
Perles à la paupière ?
Sus, rue, résilie, Lili !"
control("plus que hair", s)
"""

"""
Un Momo a est un structure qui contient tous les mots débutant par a.c
Ces mots sont une liste a.t de momos

BIGMOMO est le Momo de tous les Momo's
Il part de rien

julia BIGMOMO.c == 0x00

julia> length(BIGMOMO.t) # 26
parce qu'il y a 26 lettres possibles pour commence un mot

julia> a = BIGMOMO.t[1];
julia> a.c # 0x61
julia> Char(a.c) # 'a'
Le premier mot commence par la lettre 'a' (c'est du hasard, pas de tri)

julia> length(a.t) # 27
Tiens, une lettre de trop ?
julia> sort(map(x -> Char(x.c), a.t)) # '\0' 'a' ...
Non, le caractère 0x00 est accepté, il signifie que le mot est terminé.
C'est à dire que "a" est un mot.
"""
BIGMOMO = buildmomo()

"""
BIGMOMO est le Momo qui commence par n'importe quoi.
Il n'y a pas de mot vide
julia> nm = findc(BIGMOMO, 0x00);
julia> @assert isnothing(nm)
"""

"""
Mais il y a le mot "l"
julia> nm = findc(BIGMOMO, UInt8('l'));
julia> nm0 = findc(nm, 0x00)
julia> @assert !isnothing(nm0)

julia> z = findc(BIGMOMO, 'z');
julia> @assert z.c == UInt8('z')
julia> length(z.t) # 10
julia> zi = findc(z, 'i');
julia> length(zi.t) # 11
julia> zip = findc(zi, 'p');
julia> sort(map(x -> Char(x.c), zip.t)) # '\0' 'p' 's'
Donc "zip" existe puisqu'il peut se terminer par '\0'
Il y a aussi "zipp..." et "zips..."
julia> zips = findc(zip, 's') # Momo(0x73, Momo[Momo(0x00, Momo[])])
En fait "zips" est le seul "zips..."

La fonction explore retoure la chaine donnée complétée avec le caractère donné
et le Momo des mots qui poursuivent le Momo donné avec le caractère donné
julia> zips = explore("devrait_etre_zip_", zip, UInt8('s'))

"""


text:
module DigitSets

    import Base: <, <=, ⊊, ⊈

    # Bitmask is used as an indirection mechanism to allow us to
    # construct e.g. a DigitSet containing the single digit "3"
    # as DigitSet(3), while also providing a way to construct
    # a DigitSet directly from the binary representation of an
    # integer.
    struct Bitmask{T<:Integer}
        data::T
    end

    struct DigitSet
        d::UInt16
        DigitSet(b::Bitmask) = new(b.data)
    end

    # Not strictly necessary, but I'd like to keep the logic for single
    # element set construction as simple as possible.
    function DigitSet(n::Integer)
        d = UInt16(0) | UInt16(1 << (n - 1))
        DigitSet(Bitmask(d))
    end

    DigitSet() = DigitSet(Bitmask(0))

    function DigitSet(itr)
        d = UInt16(0)
        # For each digit in a, set the corresponding
        # bit in d to 1.
        for n in itr
            d |= UInt16(1 << (n - 1))
        end
        DigitSet(Bitmask(d))
    end

    function DigitSet(ns::Integer...)
        DigitSet(ns)
    end

    function Base.iterate(ds::DigitSet, state=nothing)
        if state == nothing
            state = (ds.d, 0)
        end

        if state[1] <= 0 return nothing end

        (d, n) = state
        shift = 1 + trailing_zeros(d)
        n += shift
        d >>= shift
        return (n, (d, n))
    end

    Base.length(ds::DigitSet) = count_ones(ds.d)
    Base.in(n, ds::DigitSet) = (ds.d & (1 << (n - 1))) != 0
    Base.isempty(ds::DigitSet) = ds.d == 0
    Base.eltype(::Type{DigitSet}) = Int

    function Base.show(io::IO, ds::DigitSet)
      print(io, "DigitSet")
      print(io, "(")
      join(io, ds, ",")
      print(io, ")")
    end

    # Set operations
    Base.union(a::DigitSet, b::DigitSet) = DigitSet(Bitmask(a.d | b.d))
    Base.intersect(a::DigitSet, b::DigitSet) = DigitSet(Bitmask(a.d & b.d))
    Base.symdiff(a::DigitSet, b::DigitSet) = DigitSet(Bitmask(xor(a.d, b.d)))
    Base.setdiff(a::DigitSet, b::DigitSet) = DigitSet(Bitmask(a.d & (~b.d)))

    # Variadic versions of union, intersect, and symdiff. All of these
    # are associative.
    function Base.union(a::DigitSet, bs::DigitSet...)
        for b in bs
            a = union(a, b)
        end
        a
    end

    function Base.intersect(a::DigitSet, bs::DigitSet...)
        for b in bs
            a = intersect(a, b)
        end
        a
    end

    function Base.symdiff(a::DigitSet, bs::DigitSet...)
        for b in bs
            a = symdiff(a, b)
        end
    end

    # Ordering and subset relationships
    Base.issubset(a::DigitSet, b::DigitSet) = isempty(setdiff(a, b))
    <(a::DigitSet, b::DigitSet) = (length(a) < length(b)) && (a <= b)
    <=(a::DigitSet, b::DigitSet) = issubset(a, b)
    ⊊(a::DigitSet, b::DigitSet) = <(a, b)
    ⊈(a::DigitSet, b::DigitSet) = !⊆(a, b)

    export DigitSet

end


text:
"""
Copyright Marcus Greiff 2020

Tests of the continuous attitude FOF controller on SO(3)
"""

# Include modules and recompile if running test standalone
if !isdefined(Main, :runAllTests)
    using Test, AerialVehicleControl, LinearAlgebra
    recompile();
end

# Testing fixture
function fixture()
    # Randomize a set of reference states
    qref = rand(4,1); qref = qref./norm(qref); wref = rand(3,1); aref = rand(3,1); fref = 0.5;
    R = ref_state_qw_t(hcat(qref), hcat(wref), hcat(aref), fref)
    # Randomize a set of states
    q = rand(4,1); q = q./norm(q); w = rand(3,1);
    S = dyn_state_qw_t(hcat(q), hcat(w))
    # Randomize the controller struct
    J  = rand_PSD(3);
    Kw = rand_PSD(3);
    ki = rand(1,3);
    Cw = rand_PSD(3);
    cw, cR = 1.0, 1.0;
    a, b, c, d = 1.0, 2.0, 3.0, 4.0;
    Y0 = rand(3,1);
    Yi = rand(3,3);
    Vi = rand(3,3);
    # Form controller object
    C  = con_state_qw_fof_t(J, Kw, ki, Cw, cw, cR, a, b, c, d, Y0, Yi, Vi);

    return R, S, C, qref, wref, aref, q, w, J, Kw, ki, Cw, cw, cR, Y0, Yi, Vi;
end

@verbose(1, "Testing continuous FOF on SO(3)...")
@testset "Test the continuous FOF controller on SO(3)" begin
    @testset "Test the crosproduct sum function" begin
        # Fixture
        N = 5;
        YA, YAc, YAm = get_matrix_and_copy(randn(3,N));
        YB, YBc, YBm = get_matrix_and_copy(randn(3,N));
        g,  gc,  gm  = get_matrix_and_copy(rand(1,N).+0.1);
        o,  oc,  om  = get_matrix_and_copy(rand(3,1));

        out_J = zeros(3,1);
        for i = 1:N
            out_J[:,1] = out_J[:,1] + g[1,i] * cross(YA[:,i], YB[:,i]);
        end

        # Test
        status = ccall((:attitude_FOF_SO3_continuous_cross_terms, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{matrix_double_t}, Ref{matrix_double_t}, Ref{matrix_double_t}, Ref{matrix_double_t},),
            YAm,
            YBm,
            gm,
            om)

        @test isapprox(status, 1);

        # Assert that none of the input arguments have changed
        @test isapprox(YA, YAc);
        @test isapprox(YB, YBc);
        @test isapprox(g,  gc );
        # Assert that the output has changed
        @test !isapprox(o,  oc );
        # Assert that the output is computed correctly
        @test isapprox(o,  out_J);
    end
    @testset "Test computation of the control signals" begin
        # Fixture
        R, S, C, qr, wr, ar, q, w, J, Kw, ki, Cw, cw, cR, Y0, Yi, Vi = fixture();
        Rc = deepcopy(R);
        Sc = deepcopy(S);
        Cc = deepcopy(C);
        dt = Cdouble(0.1);

        # Compute the controls in Julia based on the current memory
        qhat  = zeros(4,1);
        qhat[1,1]  = S.quaternion[1];
        qhat[2,1]  = S.quaternion[2];
        qhat[3,1]  = S.quaternion[3];
        qhat[4,1]  = S.quaternion[4];
        what  = zeros(3,1);
        what[1,1]  = S.attrates[1];
        what[2,1]  = S.attrates[2];
        what[3,1]  = S.attrates[3];
        k1, k2, k3 = ki[1], ki[2], ki[3];

        taur       = J*ar - skew(J*wr)*wr;
        Rr         = quat_2_SO3(qr);
        Rhat       = quat_2_SO3(qhat);
        Rmat       = quat_2_SO3(q);
        v1, v2, v3 = Vi[:,1], Vi[:,2], Vi[:,3];

        whate  = wr - what;
        we     = wr - Y0;
        wtilde = what - Y0;

        vecA   = taur;
        vecB   = skew(J*whate)*wr;
        vecC   = Kw*whate;
        vecD   = (+k1*skew(Rr'*v1)*(Rhat'*v1)
                  +k2*skew(Rr'*v2)*(Rhat'*v2)
                  +k3*skew(Rr'*v3)*(Rhat'*v3));
        tau_J  = vecA + vecB + vecC + vecD;

        deltaR = -cR * (+k1 * skew(Rhat'*v1)*(Rr'*v1 + Yi[:,1])
                        +k2 * skew(Rhat'*v2)*(Rr'*v2 + Yi[:,2])
                        +k3 * skew(Rhat'*v3)*(Rr'*v3 + Yi[:,3]));

        deltaW = (-cw * J * skew(wr)*we
                  -cw * Kw * we
                  -Cw * wtilde);

        Rhatp1_J = Rmat * exp(skew( dt * (Y0 + deltaR)));
        whatp1_J = w + dt * inv(J) * (skew(J*Y0)*Y0 + tau_J + deltaW);

        # Call
        status = ccall((:update_attitude_FOF_SO3_continuous, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{ref_state_qw_t}, Ref{dyn_state_qw_t}, Ref{con_state_qw_fof_t}, Cdouble,),
            R,
            S,
            C,
            dt)

        # test that the function is called correctly
        @test isapprox(status, 1);

        # Test that the memory that should be unchanged remains unchanged
        @test isapprox(C.status,  Cc.status );
        @test isapprox(C.gain_cw, Cc.gain_cw);
        @test isapprox(C.gain_cR, Cc.gain_cR);
        @test isapprox(C.param_a, Cc.param_a);
        @test isapprox(C.param_b, Cc.param_b);
        @test isapprox(C.param_c, Cc.param_c);
        @test isapprox(C.param_d, Cc.param_d);

        for i = 1:3
            # The torque should have changed in the controller, but
            # everything else should be the same
            @test !isapprox(C.torques[i], Cc.torques[i]);
            @test isapprox(C.gain_ki[i], Cc.gain_ki[i]);
            @test isapprox(C.measuredGyrorates[i], Cc.measuredGyrorates[i]);
            for j = 1:3
                @test isapprox(C.inertia[i][j], Cc.inertia[i][j]);
                @test isapprox(C.invinertia[i][j], Cc.invinertia[i][j]);
                @test isapprox(C.gain_Kw[i][j], Cc.gain_Kw[i][j]);
                @test isapprox(C.gain_Cw[i][j], Cc.gain_Cw[i][j]);
                @test isapprox(C.measuredDirections[i][j], Cc.measuredDirections[i][j]);
                @test isapprox(C.globalDirections[i][j], Cc.globalDirections[i][j]);
            end
        end

        # The state quaternion and attitude rates should have changed, but
        # everything in the reference struct should be the same
        for i = 1:4
            @test !isapprox(S.quaternion[i], Sc.quaternion[i]);
            @test  isapprox(R.quaternion[i], Rc.quaternion[i]);
        end
        for i = 1:3
            @test !isapprox(S.attrates[i], Sc.attrates[i]);
            @test  isapprox(R.attrates[i], Rc.attrates[i]);
            @test  isapprox(R.attaccel[i], Rc.attaccel[i]);
        end

        # Verify that the moments have been computed correctly
        for i = 1:3
            @test isapprox(tau_J[i,1], C.torques[i]);
        end

        # Verify that the estimator has performed a correct update
        qhatp1 = zeros(4,1);
        for i = 1:4
            qhatp1[i,1] = S.quaternion[i];
        end
        Rhatp1_C = quat_2_SO3(qhatp1);

        for i = 1:3
            for j = 1:3
                @test isapprox(Rhatp1_J[i,j], Rhatp1_C[i,j]);
            end
        end

        for i = 1:3
            @test isapprox(whatp1_J[i], S.attrates[i]);
        end

    end
end


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

using HDF5
using JuliaFEM

function aster_parse_nodes(section; strip_characters=true)
    nodes = Dict{Any, Vector{Float64}}()
    has_started = false
    for line in split(section, '\n')
        m = matchall(r"[\w.-]+", line)
        if (length(m) != 1) && (!has_started)
            continue
        end
        if length(m) == 1
            if (m[1] == "COOR_2D") || (m[1] == "COOR_3D")
                has_started = true
                continue
            end
            if m[1] == "FINSF"
                break
            end
        end
        if length(m) == 4
            nid = m[1]
            if strip_characters
                nid = matchall(r"\d", nid)
                nid = parse(Int, nid[1])
            end
            nodes[nid] = float(m[2:end])
        end
    end
    return nodes
end

function parse(mesh, ::Type{Val{:CODE_ASTER_MAIL}})
    model = Dict()
    header = nothing
    data = []
    for line in split(mesh, '\n')
        length(line) != 0 || continue
        info("line: $line")
        if is_aster_mail_keyword(strip(line))
            header = parse_aster_header(line)
            empty!(data)
            continue
        end
        if line == "FINSF"
            info(data)
            header = nothing
            process_aster_section!(model, join(data, ""), header, Val{header[1]})
        end
    end
    return model
end

""" Code Aster binary file (.med). """
type MEDFile
    data :: Dict
end

function MEDFile(fn)
    return MEDFile(h5read(fn, "/"))
end

function get_mesh_names(med::MEDFile)
    return sort(collect(keys(med.data["FAS"])))
end

function get_nodes(med::MEDFile, nsets, mesh_name)
    increments = keys(med.data["ENS_MAA"][mesh_name])
    @assert length(increments) == 1
    increment = first(increments)
    nodes = med.data["ENS_MAA"][mesh_name][increment]["NOE"]
    node_ids = nodes["NUM"]
    nset_ids = nodes["FAM"]
    nnodes = length(node_ids)
    node_coords = nodes["COO"]
    dim = round(Int, length(node_coords)/nnodes)
    node_coords = reshape(node_coords, nnodes, dim)'
    d = Dict{Int64}{Tuple{Symbol, Vector{Float64}}}()
    for i=1:nnodes
        nset = Symbol(nsets[nset_ids[i]])
        d[node_ids[i]] = (nset, node_coords[:, i])
    end
    return d
end

function get_node_sets(med::MEDFile, mesh_name)
    ns = Dict{Int64, Symbol}(0 => :NALL)
    if !haskey(med.data["FAS"], mesh_name)
        warn("Mesh $mesh_name not found from med file.")
        meshes = get_mesh_names(med)
        all_meshes = join(meshes, ", ")
        warn("Available meshes: $all_meshes")
        error("Mesh $mesh_name not found.")
    end
    haskey(med.data["FAS"][mesh_name], "NOEUD") || return ns
    nsets = med.data["FAS"][mesh_name]["NOEUD"]
    for nset in keys(nsets)
        k = split(nset, "_")
        nset_id = parse(Int, k[2])
        nset_name = ascii(pointer(convert(Vector{UInt8}, nsets[nset]["GRO"]["NOM"][1])))
        ns[nset_id] = Symbol(nset_name)
    end
    return ns
end

function get_element_sets(med::MEDFile, mesh_name)
    es = Dict{Int64, Symbol}()
    if !haskey(med.data["FAS"][mesh_name], "ELEME")
        return es
    end
    elsets = med.data["FAS"][mesh_name]["ELEME"]
    for elset in keys(elsets)
        k = split(elset, '_')
        elset_id = parse(Int, k[2])
        elset_name = ascii(pointer(convert(Vector{UInt8}, elsets[elset]["GRO"]["NOM"][1])))
        es[elset_id] = Symbol(elset_name)
    end
    return es
end

function get_connectivity(med::MEDFile, elsets, mesh_name)
    elsets[0] = :OTHER
    increments = keys(med.data["ENS_MAA"][mesh_name])
    @assert length(increments) == 1
    increment = first(increments)
    all_elements = med.data["ENS_MAA"][mesh_name][increment]["MAI"]
    d = Dict{Int64, Tuple{Symbol, Symbol, Vector{Int64}}}()
    for eltype in keys(all_elements)
        elements = all_elements[eltype]
        elset_ids = elements["FAM"]
        element_ids = elements["NUM"]
        nelements = length(element_ids)
        element_connectivity = elements["NOD"]
        element_dim = round(Int, length(element_connectivity)/nelements)
        element_connectivity = reshape(element_connectivity, nelements, element_dim)'
        for i=1:nelements
            eltype = Symbol(eltype)
            elco = element_connectivity[:, i]
            elset = Symbol(elsets[elset_ids[i]])
            d[element_ids[i]] = (eltype, elset, elco)
        end
    end
    return d
end

""" Parse code aster .med file.

Paramters
---------
fn
    file name to parse
mesh_name :: optional
    mesh name, if several meshes in one file

Returns
-------
Dict containing fields "nodes" and "connectivity".

"""
function parse_aster_med_file(fn, mesh_name=nothing; debug=false)
    med = MEDFile(fn)
    mesh_names = get_mesh_names(med::MEDFile)
    all_meshes = join(mesh_names, ", ")
    if mesh_name == nothing
        length(mesh_names) == 1 || error("several meshes found from med, pick one: $all_meshes")
        mesh_name = mesh_names[1]
    else
        mesh_name in mesh_names || error("Mesh $mesh_name not found from mesh file $fn. Available meshes: $all_meshes")
    end
    nsets = get_node_sets(med, mesh_name)
    elsets = get_element_sets(med, mesh_name)
    if debug
        elset_names = join(values(elsets), ", ")
        info("Code Aster .med reader: found $(length(elsets)) element sets: $elset_names")
        nset_names = join(values(nsets), ", ")
        info("Code ASter .med reader: found $(length(nsets)) node sets: $nset_names")
    end
    nodes = get_nodes(med, nsets, mesh_name)
    conn = get_connectivity(med, elsets, mesh_name)
    result = Dict("nodes" => nodes, "connectivity" => conn)
    result["nodes"] = nodes
    result["connectivity"] = conn
    return result
end

# some glues about ordering, this is still a mystery..
# http://onelab.info/pipermail/gmsh/2008/003850.html
# http://caelinux.org/wiki/index.php/Proj:UNVConvert

#global const med_connectivity = Dict{Symbol, Vector{Int}}(
#    :Tet4 => [3, 2, 1, 4],
#    :Hex8 => [4, 8, 7, 3, 1, 5, 6, 2],  # ..?
#    :Tet10 => [3, 2, 1, 4, 6, 5, 7, 10, 9, 8])

global const med_connectivity = Dict{Symbol, Vector{Int}}(
    :Tet4   => [4,3,1,2],
    :Tet10  => [4,3,1,2,10,7,8,9,6,5],
    :Wedge6 => [4,5,6,1,2,3],
    :Hex8   => [4,8,7,3,1,5,6,2],
    :Hex20  => [4,8,7,3,1,5,6,2,20,15,19,11,12,16,14,10,17,13,18,9],
    :Hex27  => [4,8,7,3,1,5,6,2,20,15,19,11,12,16,14,10,17,13,18,9,24,25,26,23,21,22,27])

# element names in CA -> element names in JuliaFEM
global const mapping = Dict(

    :PO1 => :Poi1,

    :SE2 => :Seg2,
    :SE3 => :Seg3,
    :SE4 => :Seg4,

    :TR3 => :Tri3,
    :TR6 => :Tri6,
    :TR7 => :Tri7,

    :QU4 => :Quad4,
    :QU8 => :Quad8,
    :QU9 => :Quad9,

    :TE4 => :Tet4,
    :T10 => :Tet10,

    :PE6 => :Wedge6,
    :P15 => :Wedge15,
    :P18 => :Wedge18,

    :HE8 => :Hex8,
    :H20 => :Hex20,
    :H27 => :Hex27,

    :PY5 => :Pyramid5,
    :P13 => :Pyramid13,

    )

function aster_read_mesh(fn, mesh_name=nothing; reorder_element_connectivity=true)
    result = parse_aster_med_file(fn, mesh_name)
    mesh = Mesh()
    for (nid, (nset, ncoords)) in result["nodes"]
        add_node!(mesh, nid, ncoords)
        add_node_to_node_set!(mesh, nset, nid)
    end
    for (elid, (eltype, elset, elcon)) in result["connectivity"]
        haskey(mapping, eltype) || error("Code Aster .med reader: element type $eltype not found from mapping")
        add_element!(mesh, elid, mapping[eltype], elcon)
        add_element_to_element_set!(mesh, elset, elid)
    end
    if reorder_element_connectivity
        reorder_element_connectivity!(mesh, med_connectivity)
    end
    return mesh
end

""" Code Aster result file (.rmed). """
type RMEDFile
    data :: Dict
end

function RMEDFile(fn)
    return RMEDFile(h5read(fn, "/"))
end

""" Return nodes from result med file. """
function aster_read_nodes(rmed::RMEDFile)
    increments = keys(rmed.data["ENS_MAA"]["MAIL"])
    @assert length(increments) == 1
    increment = first(increments)
    nodes = rmed.data["ENS_MAA"]["MAIL"][increment]["NOE"]
    node_names = nodes["NOM"]
    node_coords = nodes["COO"]
    nnodes = length(node_names)
    dim = round(Int, length(node_coords)/nnodes)
    node_coords = reshape(node_coords, nnodes, dim)'
    stripper(node_name) = strip(ascii(pointer(convert(Vector{UInt8}, node_name))))
    node_names = map(stripper, node_names)
    # INFO: quite safe assumption is that id is in node name, i.e. N1 => 1, N123 => 123
    node_id(node_name) = parse(matchall(r"\d+", node_name)[1])
    node_ids = map(node_id, node_names)
    nodes = Dict([j => node_coords[:,j] for j in node_ids])
    return nodes
end

""" Read nodal field from rmed file. """
function aster_read_data(rmed::RMEDFile, field_name; field_type=:NODE,
                         info_fields=true, node_ids=nothing)

    if contains(field_name, "ELGA")
        field_type = :GAUSS
    end

    if node_ids == nothing
        nodes = aster_read_nodes(rmed)
        node_ids = sort(collect(keys(nodes)))
    end

    if info_fields
        field_names = keys(rmed.data["CHA"])
        all_fields = join(field_names, ", ")
        info("results: $all_fields")
    end

    chdata = rmed.data["CHA"]["RESU____$field_name"]
    @assert length(chdata) == 1
    increment = chdata[first(keys(chdata))]
    if field_type == :NODE
        data = increment["NOE"]["MED_NO_PROFILE_INTERNAL"]["CO"]
        results = Dict([j => data[j] for j in node_ids])
    else
        error("Unable to read result of type $field_type: not implemented")
    end
    return results
end


text:
function flow(t,t0,lat,long,x,vphi,d)

# computes the velocity and acceleration field on a lat,long
# grid at time t using a mode from macWave()
#
# input
#
# t    - time (years)
# t0   - origin time (years)
# colat - latitude (degree)
# long - longitude (degree)
# x    - meridional coordinate (x = cos(theta))
# vphi - solution for vphi at CMB (vs. x)
# d    - eigenvalue ( i x frequency )

#
# output
#
# vx   - phi component of velocity at CMB
# ax   - phi components of acceleration

# some constants
sectoyear = 365.25 * 24 * 60 * 60;
tsec = (t-t0) * sectoyear;
omega = 0.7292e-4;

# add solution at poles x = -1, 1
n = length(x)
vphi_ends = Array{ComplexF64}(undef,n+2);
x_ends = Array{Float64}(undef,n+2);
x_ends[1] = -1.0;
x_ends[2:n+1] = x;
x_ends[n+2] = 1.0;
vphi_ends[1] = 0.0;
vphi_ends[2:n+1] = vphi;
vphi_ends[n+2] = 0.0;

# allocate memory for velocity and acceleration on colat grid
vr = Array{Float64,1}(undef,length(lat));
vi = Array{Float64,1}(undef,length(lat));
v = Array{Float64,2}(undef,length(lat),length(long));
a = Array{Float64,2}(undef,length(lat),length(long));

# dimensional eigenvalue
dd = d * omega;

# evaluate latitude of wave solution
theta = -acos.(x_ends).+pi/2;

#  interpolate wave on to lat,long grid
nodes = (theta,);
itp = Interpolations.interpolate(nodes,real(vphi_ends),Gridded(Linear()));
vr = itp(lat);
itp = Interpolations.interpolate(nodes,imag(vphi_ends),Gridded(Linear()));
vi = itp(lat);

# copy real part onto grid
for j = 1 : length(lat)
    for k = 1 : length(long)

       phase = exp(dd * tsec);
       v[j,k] = real( (vr[j] + vi[j]*1im)*phase);
       a[j,k] = real( (vr[j] + vi[j]*1im)*dd*phase)*sectoyear;

    end
end

# normalize
vmax = maximum(v);
v = v / vmax;
a = a / vmax;

return v,a

end

max_stars_repo_name: bbuffett/MacWaves

max_issues_repo_name: bbuffett/MacWaves

max_forks_repo_name: bbuffett/MacWaves


text:
"""
    GeneralizedInverseGaussian(a, b, p)
The *Generalized Inverse Gaussian distribution* has probability density function
```math
f(x; a,b,p) = \\frac{(a/b)^{p/2}}{2K_p(\\sqrt{ab})}x^{(p-1)}e^{-(ax + b/x)/2}
x > 0
```

```julia
GeneralizedInverseGaussian(a,b,p)   # GeneralizedInverseGaussian distribution with parametrs a, b, and p

params(d)                           # Get the parameters, i.e. (a, b, p)
```

External Links

* [Generalized inverse Gaussian distribution on Wikipedia](https://en.wikipedia.org/wiki/Generalized_inverse_Gaussian_distribution)
* [Sampling implementation paper (Hörmann & Leydold)](https://doi.org/10.1007/s11222-013-9387-3)
* [Ratio-of-Uniforms with Mode Shift sampling (Dagpunar)](https://doi.org/10.1080/03610918908812785)
"""
struct GeneralizedInverseGaussian{T<:Real} <: ContinuousUnivariateDistribution
    a::T
    b::T
    p::T
    GeneralizedInverseGaussian{T}(a::T, b::T, p::T) where T = new{T}(a,b,p)
end

function GeneralizedInverseGaussian(a::T,b::T,p::T,check_args=true) where {T<:Real}
    check_args && @check_args(GeneralizedInverseGaussian, a > zero(a) && b > zero(b))
    return GeneralizedInverseGaussian{T}(a,b,p)
end

GeneralizedInverseGaussian(a::Real, b::Real, p::Real) = GeneralizedInverseGaussian(promote(a,b,p)...)
GeneralizedInverseGaussian(a::Integer, b::Integer, p::Integer) = GeneralizedInverseGaussian(float(a), float(b), float(p))

@distr_support GeneralizedInverseGaussian 0.0 Inf

#### Conversions
function convert(::Type{GeneralizedInverseGaussian{T}}, a::Real, b::Real, p::Real) where T<:Real
    GeneralizedInverseGaussian(T(a),T(b),T(p))
end

function convert(::Type{GeneralizedInverseGaussian{T}}, d::GeneralizedInverseGaussian{S}) where {T <: Real, S<: Real}
    GeneralizedInverseGaussian(T(d.a), T(d.b), T(d.p))
end

#### Parameters

params(d::GeneralizedInverseGaussian) = (d.a, d.b, d.p)
@inline partype(d::GeneralizedInverseGaussian{T}) where {T<:Real} = T

#### Statistics

mean(d::GeneralizedInverseGaussian) = ((a,b,p) = params(d); (sqrt(b) * besselk(p+1,sqrt(a*b))) / (sqrt(a) * besselk(p,sqrt(a*b))))

mode(d::GeneralizedInverseGaussian) = ((a,b,p) = params(d); ((p - 1) + sqrt((p - 1)^2 + a*b) / a))

function var(d::GeneralizedInverseGaussian)
    (a,b,p) = params(d)
    left = besselk(p+2,sqrt(a*b)) / besselk(p, sqrt(a*b))
    right = besselk(p+1,sqrt(a*b)) / besselk(p, sqrt(a*b))
    return (b/a) * (left - right^2)
end

#### Evaluation

function pdf(d::GeneralizedInverseGaussian{T}, x::Real) where {T<:Real}
    if x > 0
        (a,b,p) = params(d)
        top = (a/b)^(p/2)
        bot = 2*besselk(p,sqrt(a*b))
        right = x^(p-1) * exp(-(a*x + (b/x))/2)
        return (top/bot)*right
    end
    return zero(T)
end

cdf(d::GeneralizedInverseGaussian{T}, x::Real) where {T<:Real} = throw(MethodError(cdf, (d, x)))

function logpdf(d::GeneralizedInverseGaussian{T}, x::Real) where {T<:Real}
    if x > 0
        (a,b,p) = params(d)
        top = (p/2)*(log(a) - log(b))
        bot = log(2*besselk(p,sqrt(a*b)))
        left = (p - 1)*log(x)
        right = -(a*x + (b/x))/2
        return top - bot + left + right
    end
    return -T(Inf)
end

#### Sampling

#Extract a sample from the GeneralizedInverseGaussian distribution 'd'.
#The sampling procedure is implemented from [1].
#Different algorithms are used for the following scenarios:

#1. β > 1 or p > 1
#2. min(1/2, (2/3)*sqrt(1-p)) ≤ a*b ≤ 1 and abs(p) < 1
#3. 0 < a*b < min(1/2, (2/3)*sqrt(1-p)) and abs(p) < 1

#[1] Hörmann, W., Leydold, J. (2014).
#Generating generalized inverse Gaussian random variates.
#Stat Comput 24, 547–557.
#https://doi.org/10.1007/s11222-013-9387-3
function rand(rng::Random.AbstractRNG,d::GeneralizedInverseGaussian)
    (a,b,p) = params(d)
    α = sqrt(a/b)
    β = sqrt(a*b)
    λ = abs(p)
    if β > 1 || λ > 1
        x = sample_unif_mode_shift(λ,β)
    else
        β_bound = min(1/2, (2/3)*sqrt(1 - p))
        if β < 1 && β >= β_bound
            x = sample_unif_no_mode_shift(λ,β)
        elseif β < β_bound && β > 0
            x = concave_sample(λ,β)
        else
            throw(ArgumentError("None of the required conditions on the parameters are satisfied"))
        end
    end
    if p >= 0
        return x/α
    else
        return 1 / (α*x)
    end
end

#Sample the 2-parameter GeneralizedInverseGaussian distribution with parameters p and β using the Rejection method
#as described by Hörmann & Leydold (2014).
function concave_sample(p::Real,β::Real)
    m = β/((1 - p) + sqrt((1 - p)^2 + β^2))
    x_naut = β/(1 - p)
    x_star = max(x_naut,2/β)
    k1 = g(m,p,β)
    A1 = k1 * x_naut
    k2 = 0
    A2 = 0
    if x_naut < 2/β
        k2 = exp(-β)
        if p > 0
            A2 = k2 * ((2/β)^p - x_naut^p) / p
        else
            A2 = k2 * log(2/β^2)
        end
    end
    k3 = x_star^(p-1)
    A3 = 2 * k3 * exp(-x_star * β / 2) / β
    A = A1 + A2 + A3

    while true
        u = rand(Uniform(0,1))
        v = rand(Uniform(0,1))*A
        h = Inf
        if v <= A1
            x = x_naut * v / A1
            h = k1
        elseif v <= A1 + A2
            v = v - A1
            if p > 0
                x = (x_naut^p + (v*p/k2))^(1/p)
            else
                x = β*exp(v * exp(β))
            end
            h = k2 * x^(p-1)
        else
            v = v - (A1 + A2)
            x = -2 * log(exp(-x_star * β / 2) - (v * β) / (2 * k3)) / β
            h = k3 * exp(-x * β / 2)
        end
        if (u * h) <= g(x,p,β)
            return x
        end
    end
end

#Sample the 2-parameter GeneralizedInverseGaussian distribution with parameters p and β using the
#Ratio-of-Uniforms without mode shift as described by Hörmann & Leydold (2014).
function sample_unif_no_mode_shift(p::Real,β::Real)
    m = β/((1 - p) + sqrt((1 - p)^2 + β^2))
    x⁺= ((1 + p) + sqrt((1 + p)^2 + β^2))/β
    v⁺= sqrt(g(m,p,β))
    u⁺= x⁺ * sqrt(g(x⁺,p,β))
    while true
        u = rand(Uniform(0,1))*u⁺
        v = rand(Uniform(0,1))*v⁺
        x = u/v
        if v^2 <= g(x,p,β)
            return x
        end
    end
end

#Sample the 2-parameter GeneralizedInverseGaussian distribution with parameters p and β using the
#Ratio-of-Uniforms with mode shift as described by Hörmann & Leydold (2014) (originally from Dagpunar (1989)).
function sample_unif_mode_shift(p::Real,β::Real)
    m = (sqrt((p - 1)^2 + β^2) + (p-1)) / β
    a = -(2*(p+1)/β) - m
    b = (2*(p-1)/β) * m - 1
    p2 = b - (a^2)/3
    q = (2*a^3)/27 - (a*b/3) + m
    ϕ = acos(-(q/2) * sqrt(-27 / (p2^3)))
    x⁻ = sqrt((-4/3)*p2)*cos(ϕ/3 + (4/3)*π) - (a / 3)
    x⁺ = sqrt((-4/3)*p2)*cos(ϕ/3) - (a/3)
    v⁺ = sqrt(g(m,p,β))
    u⁻ = (x⁻ - m)*sqrt(g(x⁻,p,β))
    u⁺ = (x⁺ - m)*sqrt(g(x⁺,p,β))

    while true
        u = rand(Uniform(0,1))*(u⁺ - u⁻) + u⁻
        v = rand(Uniform(0,1))*v⁺
        x = (u / v) + m
        if x > 0 && v^2 <= g(x,p,β)
            return x
        end
    end
end

function g(x::Real,p::Real,β::Real)
    x^(p-1)*exp(-(β/2)*(x + (1/x)))
end


text:
"""
    solve_tsp(g, w; cutoff=Inf)

Given a graph `g` and an edgemap `w` containing the weight associated to each edge,
solves the Travelling Salesman Problem, returning the tour among all the cities with
minimal total weight.

Edges `e` with no associated value in `w` or with `w[e] > cutoff`
will not be considered for the optimal tour.

The algorithm uses a generic Integer Programming solver combined with a lazy
constraint augmentation procedure.

The package JuMP.jl and one of its supported solvers are required.

Returns a tuple `(status, W, tour)` containing:
- a solve `status` (indicating whether the problem was solved to optimality)
- the tototal weight `W` of the tour
- a vector `tour` containing the sequence of vertices in the optimal tour.

**Example**
```julia
pos = rand(30, 2)
g = CompleteGraph(30)
w = EdgeMap(g, e -> norm(pos[src(e)] - pos[dst(e)]))
status, W, tour = solve_tsp(g, w)
```
"""
function solve_tsp(g::G, w::AEdgeMap; cutoff=Inf, verb=false) where {G<:AGraph}
    h = G(nv(g))
    for e in edges(g)
        haskey(w, e) && w[e] <= cutoff && add_edge!(h, e)
    end
    return _solve_tsp(h, w, verb)
end

function _solve_tsp(g::AGraph, w::AEdgeMap, verb::Bool)
    elist = collect(edges(g))
    model = Model(with_optimizer(MIP_OPTIMIZER))
    @variable(model, y[elist], Bin)
    @objective(model, Min, sum(y[e]*w[e] for e in edges(g)))
    @constraint(model, c1[i=1:nv(g)], sum(y[sort(e)] for e in edges(g, i)) == 2)

    ## TODO add lazy callback when avalaible again in JuMP
    ## (see https://github.com/JuliaOpt/JuMP.jl/pull/1849)
    # addlazycallback(model,
    #     cb -> begin
    #         sol = getvalue(y)
    #         alltours = gettours(g, sol)
    #         verb && println("Found $(length(alltours)) tours of lengths $(length.(alltours)). Adding lazy constraints.")
    #         for tour in alltours
    #             length(tour) == nv(g) && break
    #             add_tour_constraint(cb, y, g, tour)
    #         end
    #     end)
    has_short_tours = true
    while has_short_tours
        optimize!(model)
        sol = Dict(e => value(y[e]) for e in elist)
        alltours = gettours(g, sol)
        if length(alltours) == 1
            @assert length(alltours[1]) == nv(g)
            has_short_tours = false
        else
            verb && println("@ Found $(length(alltours)) tours.")
            for tour in alltours
                verb && println("@@ Add constraint on tour of length $(length(tour)).")
                add_tour_constraint(model, y, g, tour)
            end
        end
    end
    status = termination_status(model)
    sol = Dict(e => value(y[e]) for e in elist)
    alltours = gettours(g, sol)
    cost = objective_value(model)
    return status, cost, alltours[1]
end

function gettours(g, sol; atol=1e-8)
    alltours = Vector{Vector{Int}}()
    notintour = [1:nv(g);]
    while !isempty(notintour)
        root = first(notintour)
        tour = Int[root]
        u  = root
        prev = -1
        while u != root || length(tour) == 1
            for e in out_edges(g, u)
                abs(sol[sort(e)]) < atol && continue
                v = dst(e)
                v == prev && continue
                push!(tour, v)
                u, prev = v, u
                break
            end
        end
        resize!(tour, length(tour)-1)
        @assert length(tour) > 2
        notintour = setdiff(notintour, tour)
        @assert length(notintour) > 2 || length(notintour) == 0 "notintour=$notintour"
        push!(alltours, tour)
    end
    @assert sum(length.(alltours)) == nv(g)
    return alltours
end

function add_tour_constraint(model, y, g, tour)
    aff = AffExpr()
    for (k, v) in enumerate(tour)
        vprev = k == 1 ? tour[end] : tour[k-1]
        vnext = k == length(tour) ? tour[1] : tour[k+1]
        for e in edges(g, v)
            if dst(e) ∉ [vprev,vnext]
                aff += y[sort(e)]
            end
        end
    end
    @constraint(model, aff >= 2)
end


text:
module SeriesIterators

export iterm, icumsum
export bernoulli

using IterTools
using Base.Iterators

import Base: iterate, IteratorSize, length, size

"""
    iterm(g, ind)

Assume a sequence of terms a[k] can be generated by

    a[i0], b = g(i0)
    a[i], b =  g(i, b) for subsequent indices

The `iterm(g, ind)` iterates the sequence, when `ind` is the
iterator of the indices; first index is `i0`.
"""
iterm(gen, it) = TermIterator(gen, it)
struct TermIterator{F,I}
    generate::F
    indices::I
end
function iterate(t::TermIterator)
    vs = iterate(t.indices)
    k, s = @IterTools.ifsomething vs
    v = t.generate(k)
    v[1], (s, v)
end
function iterate(t::TermIterator, s)
    sx, p = s
    vs = iterate(t.indices, sx)
    k, s = @IterTools.ifsomething vs
    v = t.generate(k, p)
    v[1], (s, v)
end
IteratorSize(::Type{<:TermIterator{<:Any,I}}) where I = IteratorSize(I)
length(t::TermIterator) = length(t.indices)
size(t::TermIterator) = size(t.indices)

function gen_bernoulli(::Type{T}) where T<:Integer
    B = Vector{Rational{T}}(undef, 2)
    B[1] = 1
    B[2] = -1 // 2
    Z = zero(Rational{T})
    function _bernoulli(n::Integer)
        b = B
        n <= 1 && return b[n+1]
        isodd(n) && return Z
        l = length(b)
        if l * 2 - 2 <= n
            ext = min(1_000_000 ÷ l^2, 40)
            m = max(n ÷ 2 + 2, l + ext)
            n0 = (length(b) - 1 ) * 2
            sizehint!(b, ((m + 80) >> 6) << 6)
            resize!(b, m)
            _bernoulli!(b, n0, (m - 1) * 2 - 1, n)
        end
        b[n÷2 + 2]
    end
    _bernoulli
end

function _bernoulli!(B::Vector{Rational{T}}, n0::Integer, n::Integer, nn::Integer) where T
    lastx = 0
    try
        for m = n0:2:n
            sum = -B[1] // T(m + 1)
            sum -= B[2]
            mok = T((m * (m - 1)) ÷ 2)
            for k = 2:2:m-1
                sum -= B[k÷2 + 2] // T(m - k + 1) * mok
                mok *= T(m - k) // T(k + 1)
                mok *= T(m - k - 1) // T(k + 2)
            end
            lastx = m ÷ 2 + 2
            B[lastx] = sum
        end
    catch
        if lastx > 0
            resize!(B, lastx)
        end
        if lastx < nn ÷ 2 + 2
            rethrow()
        end
    end
    nothing
end
"""
    bernoulli(n)

Return Bernoulli number `B_n`. `B_1 = -1/2`.
"""
const bernoulli = gen_bernoulli(BigInt)

struct CumSum{T}
    itr::T
end

"""
    icumsum(itr)

Return iterator, which lazily delivers the same results as `cumsum(itr)`.
"""
icumsum(itr) = CumSum(itr)
function iterate(c::CumSum)
    v, s = @IterTools.ifsomething iterate(c.itr)
    v, (s, v)
end
function iterate(c::CumSum, (s, p))
    v, s = @IterTools.ifsomething iterate(c.itr, s)
    v = p + v
    v, (s, v)
end
IteratorSize(::Type{<:CumSum{T}}) where T = IteratorSize(T)
length(c::CumSum) = length(c.itr)
size(c::CumSum) = size(c.itr)

end


text:
using BSONLite
using DataStructures
using JSON
using Test

test_names = [
              "array",
              "binary",
              "boolean",
              "code",
              "code_w_scope",
              "datetime",
              "dbpointer",
              "dbref",
              "decimal128-1",
              "decimal128-2",
              "decimal128-3",
              "decimal128-4",
              "decimal128-5",
              "decimal128-6",
              "decimal128-7",
              "document",
              "double",
              "int32",
              "int64",
              "maxkey",
              "minkey",
              "multi-type",
              "multi-type-deprecated",
              "null",
              "oid",
              "regex",
              "string",
              "symbol",
              "timestamp",
              "top",
              "undefined",
             ]

function test_corpus(name)
    path = joinpath(@__DIR__, "bson-corpus/tests", "$name.json")
    case = JSON.parse(String(read(path)))
    @testset "$(case["description"])" begin

        if haskey(case, "valid")
            @testset "Read/Write" begin
                for v in case["valid"]
                    @testset "$(v["description"])" begin
                        bson_str = v["canonical_bson"]
                        bson_bin = hex2bytes(bson_str)

                        # test minimal round trip read/write

                        result = bytes2hex(write_bson(read_bson(bson_bin, codec = :minimal)))
                        @test uppercase(result) == uppercase(bson_str)

                        codec = BSONCodec(OrderedDict)
                        result = bytes2hex(write_bson(read_bson(bson_bin, codec = codec)))
                        @test uppercase(result) == uppercase(bson_str)

                        if match(r"^decimal-", name) != nothing
                            json_str = v["canonical_extjson"]
                            json_dat = JSON.parse(json_str)
                            result = read_bson(bson_bin, codec = :json)
                            @test result == json_dat
                        end
                    end
                end
            end
        end

        if haskey(case, "decodeErrors")
            @testset "Decode errors" begin
                for v in case["decodeErrors"]
                    err = v["description"]
                    @testset "$err" begin
                        @test_throws Exception read_bson(hex2bytes(v["bson"]))
                    end
                end
            end
        end
    end
end

@testset "bson-corpus" begin
    test_corpus.(test_names)
end


text:
using Revise
using Plots
using DynamicalSystems
using Basins

# Multistability, phase diagrams, and intransitivity in the Lorenz-84 low-order atmospheric circulation model
# Chaos 18, 033121 (2008); https://doi.org/10.1063/1.2953589
@inline @inbounds function lorenz84(u, p, t)
    F = p[1]; G = p[2]; a = p[3]; b = p[4];
	x = u[1]; y = u[2]; z = u[3];
    dx = -y^2 -z^2 -a*x + a*F
    dy = x*y - y - b*x*z +G
	dz = b*x*y + x*z -z
    return SVector{3}(dx, dy, dz)
end

F=6.846; G=1.287; a=0.25; b=4.;
p= [F, G, a, b]
ds = ContinuousDynamicalSystem(lorenz84, rand(3), p)
xg=range(-1.,1.,length=200)
yg=range(-1.5,1.6,length=200)
pmap = poincaremap(ds, (3, 0.), Tmax=1e6; idxs = 1:2, rootkw = (xrtol = 1e-8, atol = 1e-8), reltol=1e-9)
@time bsn = basins_map2D(xg, yg, pmap)
plot(xg,yg,bsn.basin',seriestype=:heatmap)


text:
"""
    compute_fem_laplace_term1(u::Union{PyObject, Array{Float64, 1}},
        nu::Union{PyObject, Array{Float64, 1}},
        mesh::Mesh3)
"""
function compute_fem_laplace_term1(u::Union{PyObject, Array{Float64, 1}},
    nu::Union{PyObject, Array{Float64, 1}},
    mesh::Mesh3)
    @assert length(u) == mesh.ndof
    @assert length(nu) == get_ngauss(mesh)
    compute_laplace_term_mfem_ = load_op_and_grad(AdFem.libmfem3,"compute_laplace_term_mfem_t")
    u,nu = convert_to_tensor(Any[u,nu], [Float64,Float64])
    out = compute_laplace_term_mfem_(u,nu)
    set_shape(out, (mesh.ndof,))
end

##########################################

"""
    compute_fem_laplace_matrix1(kappa::PyObject, mesh::Mesh3)
"""
function compute_fem_laplace_matrix1(kappa::PyObject, mesh::Mesh3)
    @assert length(kappa) == get_ngauss(mesh)
    fem_laplace_scalar_ = load_op_and_grad(AdFem.libmfem3,"fem_laplace_scalar_t", multiple=true)
    kappa = convert_to_tensor(Any[kappa], [Float64]); kappa = kappa[1]
    indices, vv = fem_laplace_scalar_(kappa)
    n = mesh.ndof
    RawSparseTensor(indices, vv, n, n)
end

"""
    compute_fem_laplace_matrix1(kappa::Array{Float64,1}, mesh::Mesh3)
"""
function compute_fem_laplace_matrix1(kappa::Array{Float64,1}, mesh::Mesh3)
    @assert length(kappa) == get_ngauss(mesh)
    N = get_ngauss(mesh) * size(mesh.conn, 2)^2
    indices = zeros(Int64, 2N)
    vv = zeros(N)
    @eval ccall((:FemLaplaceScalarT_forward_Julia, $LIBMFEM3), Cvoid, (Ptr{Int64}, Ptr{Cdouble}, Ptr{Cdouble}), $indices, $vv, $kappa)
    indices = reshape(indices, 2, N)'|>Array
    sparse(indices[:,1] .+ 1, indices[:,2] .+ 1, vv, mesh.ndof, mesh.ndof)
end

"""
    compute_fem_laplace_matrix(kappa::Union{PyObject, Array{Float64, 1}}, mesh::Mesh3)
"""
function compute_fem_laplace_matrix(kappa::Union{PyObject, Array{Float64, 1}}, mesh::Mesh3)
    Z = compute_fem_laplace_matrix1(kappa, mesh)
    if isa(Z, SparseMatrixCSC)
        [Z spzeros(mesh.ndof, mesh.ndof)
        spzeros(mesh.ndof, mesh.ndof) Z]
    else
        [Z spzero(mesh.ndof)
        spzero(mesh.ndof) Z]
    end
end

compute_fem_laplace_matrix1(mesh::Mesh3) = compute_fem_laplace_matrix1(ones(get_ngauss(mesh)), mesh)
compute_fem_laplace_matrix(mesh::Mesh3) = compute_fem_laplace_matrix(ones(get_ngauss(mesh)), mesh)
##########################################

"""
    compute_fem_source_term1(f::PyObject, mesh::Mesh3)
"""
function compute_fem_source_term1(f::PyObject, mesh::Mesh3)
    @assert length(f)==get_ngauss(mesh)
    fem_source_scalar_ = load_op_and_grad(libmfem3,"fem_source_scalar_t")
    f = convert_to_tensor(Any[f], [Float64]); f = f[1]
    out = fem_source_scalar_(f)
    n = mesh.ndof
    set_shape(out, (n,))
end

"""
    compute_fem_source_term1(f::Array{Float64,1}, mesh::Mesh3)
"""
function compute_fem_source_term1(f::Array{Float64,1}, mesh::Mesh3)
    @assert length(f)==get_ngauss(mesh)
    out = zeros(mesh.ndof)
    @eval ccall((:FemSourceScalarT_forward_Julia, $LIBMFEM3), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}), $out, $f)
    out
end

"""
    compute_fem_source_term(f1::Union{PyObject,Array{Float64,2}},
        f2::Union{PyObject,Array{Float64,2}}, mesh::Mesh3)
"""
function compute_fem_source_term(f1::Union{PyObject,Array{Float64,1}},
    f2::Union{PyObject,Array{Float64,1}},
    f3::Union{PyObject,Array{Float64,1}},
    mesh::Mesh3)
    @assert length(f1)==length(f2)==length(f3)==get_ngauss(mesh)
    [compute_fem_source_term1(f1, mesh); compute_fem_source_term1(f2, mesh); compute_fem_source_term1(f3, mesh)]
end

function compute_fem_source_term(f::Union{PyObject,Array{Float64,1}},  mesh::Mesh3)
    @assert length(f)==3get_ngauss(mesh)
    N = get_ngauss(mesh)
    f1 = f[1:N]
    f2 = f[N+1:2N]
    f3 = f[2N+1:end]
    compute_fem_source_term(f1, f2, f3, mesh)
end

##########################################

"""
    compute_fem_mass_matrix1(ρ::Union{PyObject, Array{Float64, 1}},
            mmesh::Mesh3)
"""
function compute_fem_mass_matrix1(ρ::Union{PyObject, Array{Float64, 1}},
        mmesh::Mesh3)
        compute_fem_mass_matrix_mfem_t_ = load_op_and_grad(libmfem3,"compute_fem_mass_matrix_mfem_t", multiple=true)
        rho = convert_to_tensor(Any[ρ], [Float64]); rho = rho[1]
        indices, vv = compute_fem_mass_matrix_mfem_t_(rho)
        RawSparseTensor(indices, vv, mmesh.ndof, mmesh.ndof)
end

"""
    compute_fem_mass_matrix1(mmesh::Mesh3)
"""
compute_fem_mass_matrix1(mmesh::Mesh3) = compute_fem_mass_matrix1(ones(get_ngauss(mmesh)), mmesh)


text:
# ---
# title: 508. Most Frequent Subtree Sum
# id: problem508
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, Tree
# link: <https://leetcode.com/problems/most-frequent-subtree-sum/description/>
# hidden: true
# ---
#
# Given the root of a tree, you are asked to find the most frequent subtree sum.
# The subtree sum of a node is defined as the sum of all the node values formed
# by the subtree rooted at that node (including the node itself). So what is the
# most frequent subtree sum value? If there is a tie, return all the values with
# the highest frequency in any order.
#
# **Examples 1**
# Input:
#
#
#
#       5
#      /  \
#     2   -3
#
#
# return [2, -3, 4], since all the values happen only once, return all of them
# in any order.
#
# **Examples 2**
# Input:
#
#
#
#       5
#      /  \
#     2   -5
#
#
# return [2], since 2 happens twice, however -5 only occur once.
#
# **Note:** You may assume the sum of values in any subtree is in the range of
# 32-bit signed integer.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Archibald
using Test
using Distributions
using DataFrames
using LinearAlgebra

@testset "Over Sampling MCMC" begin

    function generate_data(;n = 10,
                            T = 100,
                            β = 0.25,
                            family = Bernoulli,
                            τ = 1,
                            link = x -> cdf(Normal(), x))
        i_index = repeat(1:n, inner = T)
        t_index = repeat(1:T, n)
        μ_x = rand(Normal(0, 1), n)
        α_i = rand(Normal(0, τ), n)
        x = vcat(rand.(Normal.(μ_x, 1), T )...)
        μ_vec = α_i[i_index] .+ x * β

        y = rand.(family.(link.(μ_vec)))

        return y, x, i_index, t_index, α_i
    end

    function linear_target_dist_oversample(y,
                                           x,
                                           i_index,
                                           single_parameters,
                                           over_parameters)
        β = single_parameters[1]
        τ = exp(single_parameters[2])
        α_matrix = over_parameters
        S = size(α_matrix, 2)
        ll = Vector{Float64}(undef, S)
        i_index = Int.(i_index)
        for s in 1:S
            resid_vec = y .- α_matrix[i_index, s] .- x*β
            ll[s] = sum(log.(pdf.(Normal.(0, 1), resid_vec)))
        end
        τ_prior = TruncatedNormal(0, 3, 0, Inf)
        τ_prior_dens = pdf(τ_prior, τ)
        param_probability = sum(ll) + log(τ_prior_dens)
        return param_probability
    end

    sim_n = 5
    sim_t = 100
    sim_τ = 0.5
    linear_data = generate_data(n = sim_n, T = sim_t, τ=sim_τ, link = x -> Normal(x, 1), family = identity)

    oversample_parameters = randn(sim_n, 3)
    chain_over = Archibald.MarkovChainOverSample(
        [0.1, 0.2],
        3,
        oversample_parameters,
        [linear_data[1], linear_data[2], linear_data[3]],
        (x_single, x_over, data) -> linear_target_dist_oversample(data[1],
                                                   data[2],
                                                   data[3],
                                                   x_single, x_over)
    )
    MetropolisHastings!(
        chain_over,
        (x) -> MvNormal(x, 0.1),
        (x) -> MvNormal(sim_n, x[2]),
        logs = true
        )

    mvn_draws = MetropolisHastings(
        1_000,
        chain_over,
        x -> MvNormal(x, 0.1), # jeez what an elegant design right? right?
        x -> MvNormal(sim_n, x[2]),
        logs = true,
        burn = 1_000)
    # now "for real" with lower proposal variance
    mvn_df = tidybayes(mvn_draws...)
    mvn_df
    @test isa(mvn_df, DataFrame)
    @test isa(mvn_draws[1][1], Vector)
    @test isa(mvn_draws[1][1][1], Float64)
    @test !isnan(mvn_draws[1][1][2])

end

@testset "Basic Metropolis Hastings" begin
    function generate_data(; N = 100)
        sigma_b = log(rand(InverseGamma(2, 2)))
        sigma_e = log(rand(InverseGamma(1, 1)))
        β = [1, 2, 3, 4]

        X = randn(2, 4, N)
        W = randn(2, 1, N)
        b = randn(1, N).*sigma_b
        e = randn(2, N).*sigma_e

        Y = Matrix{Float64}(undef, 2, N)
        for i in 1:N
            Y[:, i] = X[:, :, i] * β + W[:, :, i]*b[:, i] + e[:, i]
        end
        return Y, X, W, b, β, e, sigma_b, sigma_e
    end

    """
    Create target distribution for regression problem
    """
    function mvn_lik(y, X, W, parameters)
        N = size(X, 3)
        β = parameters[1:4]
        # constrain sigmas to be positive using exp(.)
        σ_b = exp(parameters[5])
        σ_e = exp(parameters[6])

        ll_i = Vector{Float64}(undef, N)
        for i in 1:N
            # Account for variance from Wb_i
            Σ_sqrt = σ_b^2 * W[:, :, i]*W[:, :, i]' + σ_e^2 .* Matrix(I,2,2)
            Σ = (Σ_sqrt + Σ_sqrt')./2 # make it Hermitian
            resid = y[:, i] - X[:, :, i]*β
            ll_i[i] = log(pdf(MvNormal([0, 0], Σ), resid))
        end

        # sum log lik contribution
        ll = sum(ll_i)
        # priors
        σ_b_prior = pdf(InverseGamma(2, 2), σ_b)
        σ_e_prior = pdf(InverseGamma(1, 1), σ_e)

        P = ll + log(sum(pdf.(Normal(0, 10), β))) + log(σ_b_prior) + log(σ_e_prior)
        return P
    end

    test_Y, test_X, test_W, test_b, test_β, test_e, test_σ_b, test_σ_e = generate_data(N = 1_00)

    # wrapper function to match arguments with previous P()
    mvn_target_pdf(x, data) =  mvn_lik(data[1], data[2], data[3], x)

    # Save true values for comparison later
    true_values = vcat(test_β, exp(test_σ_b), exp(test_σ_e), test_b[:])[:]
    true_params = DataFrame(
        [
            :true_value => true_values,
            :term => "parameter_" .* string.(1:size(true_values, 1)),
            :term_name => [
                "β_1",
                "β_2",
                "β_3",
                "β_4",
                "σ_b",
                "σ_e",
                "b_" .* string.(1:(size(true_values, 1) - 6))...
            ]
        ]
    )

    # Initialise chain
    mvn_markov_state = MarkovChain(
        fill(1.0, size(true_params, 1)),
        [test_Y, test_X, test_W],
        mvn_target_pdf
    )
    # burn in ad hoc for 20k
    mvn_draws = MetropolisHastings(
        1_000,
        mvn_markov_state,
        x -> MvNormal(x, 0.1), # jeez what an elegant design right? right?
        logs = true,
        burn = 1_000)
    # now "for real" with lower proposal variance
    mvn_df = tidybayes(mvn_draws...)
    @test isa(mvn_df, DataFrame)
    @test size(mvn_df, 1) == 106000
    @test isa(mvn_draws[1][1], Vector)
    @test isa(mvn_draws[1][1][1], Float64)
    @test !isnan(mvn_draws[1][1][2])
end

max_stars_repo_name: EdJeeOnGitHub/Archibald

max_issues_repo_name: EdJeeOnGitHub/Archibald

max_forks_repo_name: EdJeeOnGitHub/Archibald