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

using Test, Distributed, Random

Foo_module = :Foo4b3a94a1a081a8cb
Foo2_module = :F2oo4b3a94a1a081a8cb
FooBase_module = :FooBase4b3a94a1a081a8cb
@eval module ConflictingBindings
    export $Foo_module, $FooBase_module
    $Foo_module = 232
    $FooBase_module = 9134
end
using .ConflictingBindings

function precompile_test_harness(@nospecialize(f), testset::String)
    @testset "$testset" begin
        precompile_test_harness(f, true)
    end
end
function precompile_test_harness(@nospecialize(f), separate::Bool)
    load_path = mktempdir()
    load_cache_path = separate ? mktempdir() : load_path
    try
        pushfirst!(LOAD_PATH, load_path)
        pushfirst!(DEPOT_PATH, load_cache_path)
        f(load_path)
    finally
        rm(load_path, recursive=true, force=true)
        separate && rm(load_cache_path, recursive=true, force=true)
        filter!((≠)(load_path), LOAD_PATH)
        separate && filter!((≠)(load_cache_path), DEPOT_PATH)
    end
    nothing
end

# method root provenance

rootid(m::Module) = ccall(:jl_module_build_id, UInt64, (Any,), Base.parentmodule(m))
rootid(m::Method) = rootid(m.module)

function root_provenance(m::Method, i::Int)
    mid = rootid(m)
    isdefined(m, :root_blocks) || return mid
    idxs = view(m.root_blocks, 2:2:length(m.root_blocks))
    j = searchsortedfirst(idxs, i) - 1   # RLE roots are 0-indexed
    j == 0 && return mid
    return m.root_blocks[2*j-1]
end

struct RLEIterator{T}   # for method roots, T = UInt64 (even on 32-bit)
    items::Vector{Any}
    blocks::Vector{T}
    defaultid::T
end
function RLEIterator(roots, blocks, defaultid)
    T = promote_type(eltype(blocks), typeof(defaultid))
    return RLEIterator{T}(convert(Vector{Any}, roots), blocks, defaultid)
end
RLEIterator(m::Method) = RLEIterator(m.roots, m.root_blocks, rootid(m))
Base.iterate(iter::RLEIterator) = iterate(iter, (0, 0, iter.defaultid))
function Base.iterate(iter::RLEIterator, (i, j, cid))
    i += 1
    i > length(iter.items) && return nothing
    r = iter.items[i]
    while (j + 1 < length(iter.blocks) && i > iter.blocks[j+2])
        cid = iter.blocks[j+1]
        j += 2
    end
    return cid => r, (i, j, cid)
end

function group_roots(m::Method)
    mid = rootid(m)
    isdefined(m, :root_blocks) || return Dict(mid => m.roots)
    group_roots(RLEIterator(m.roots, m.root_blocks, mid))
end
function group_roots(iter::RLEIterator)
    rootsby = Dict{typeof(iter.defaultid),Vector{Any}}()
    for (id, r) in iter
        list = get!(valtype(rootsby), rootsby, id)
        push!(list, r)
    end
    return rootsby
end

precompile_test_harness("basic precompile functionality") do dir2
precompile_test_harness(false) do dir
    Foo_file = joinpath(dir, "$Foo_module.jl")
    Foo2_file = joinpath(dir, "$Foo2_module.jl")
    FooBase_file = joinpath(dir, "$FooBase_module.jl")

    write(FooBase_file,
          """
          false && __precompile__(false)
          module $FooBase_module
              import Base: hash, >
              struct fmpz end
              struct typeA end
              >(x::fmpz, y::Int) = Base.cmp(x, y) > 0
              function hash(a::typeA, h::UInt)
                  d = den(a)
                  return h
              end
          end
          """)
    write(Foo2_file,
          """
          module $Foo2_module
              export override
              override(x::Integer) = 2
              override(x::AbstractFloat) = Float64(override(1))
          end
          """)
    write(Foo_file,
          """
          module $Foo_module
              import $FooBase_module, $FooBase_module.typeA
              import $Foo2_module: $Foo2_module, override
              import $FooBase_module.hash
              import Test
              module Inner
                  import $FooBase_module.hash
                  using ..$Foo_module
                  import ..$Foo2_module
              end

              struct typeB
                  y::typeA
              end
              hash(x::typeB) = hash(x.y)

              # test that docs get reconnected
              @doc "foo function" foo(x) = x + 1
              include_dependency("foo.jl")
              include_dependency("foo.jl")
              module Bar
                  include_dependency("bar.jl")
              end
              @doc "Bar module" Bar # this needs to define the META dictionary via eval
              @eval Bar @doc "bar function" bar(x) = x + 2

              # test for creation of some reasonably complicated type
              struct MyType{T} end
              const t17809s = Any[
                    Tuple{
                        Type{Ptr{MyType{i}}},
                        Ptr{Type{MyType{i}}},
                        Array{Ptr{MyType{MyType{:sym}()}}(0), 0},
                        Val{Complex{Int}(1, 2)},
                        Val{3},
                        Val{nothing}}
                    for i = 0:25]

              # test that types and methods get reconnected correctly
              # issue 16529 (adding a method to a type with no instances)
              (::Task)(::UInt8, ::UInt16, ::UInt32) = 2

              # issue 16471 (capturing references to a kwfunc)
              Test.@test !isdefined(typeof(sin).name.mt, :kwsorter)
              Base.sin(::UInt8, ::UInt16, ::UInt32; x = 52) = x
              const sinkw = Core.kwfunc(Base.sin)

              # issue 16908 (some complicated types and external method definitions)
              abstract type CategoricalPool{T, R <: Integer, V} end
              abstract type CategoricalValue{T, R <: Integer} end
              struct NominalPool{T, R <: Integer, V} <: CategoricalPool{T, R, V}
                  index::Vector{T}
                  invindex::Dict{T, R}
                  order::Vector{R}
                  ordered::Vector{T}
                  valindex::Vector{V}
              end
              struct NominalValue{T, R <: Integer} <: CategoricalValue{T, R}
                  level::R
                  pool::NominalPool{T, R, NominalValue{T, R}}
              end
              struct OrdinalValue{T, R <: Integer} <: CategoricalValue{T, R}
                  level::R
                  pool::NominalPool{T, R, NominalValue{T, R}}
              end
              (::Union{Type{NominalValue}, Type{OrdinalValue}})() = 1
              (::Union{Type{NominalValue{T}}, Type{OrdinalValue{T}}})() where {T} = 2
              (::Type{Vector{NominalValue{T, R}}})() where {T, R} = 3
              (::Type{Vector{NominalValue{T, T}}})() where {T} = 4
              (::Type{Vector{NominalValue{Int, Int}}})() = 5

              # more tests for method signature involving a complicated type
              # issue 18343
              struct Pool18343{R, V}
                  valindex::Vector{V}
              end
              struct Value18343{T, R}
                  pool::Pool18343{R, Value18343{T, R}}
              end
              Base.convert(::Type{Some{S}}, ::Value18343{Some}) where {S} = 2
              Base.convert(::Type{Some{Value18343}}, ::Value18343{Some}) = 2
              Base.convert(::Type{Ref}, ::Value18343{T}) where {T} = 3

              # issue #28297
              mutable struct Result
                  result::Union{Int,Missing}
              end

              const x28297 = Result(missing)

              const d29936a = UnionAll(Dict.var, UnionAll(Dict.body.var, Dict.body.body))
              const d29936b = UnionAll(Dict.body.var, UnionAll(Dict.var, Dict.body.body))

              # issue #28998
              const x28998 = [missing, 2, missing, 6, missing,
                              missing, missing, missing,
                              missing, missing, missing,
                              missing, missing, 6]

              let some_method = which(Base.include, (Module, String,))
                    # global const some_method // FIXME: support for serializing a direct reference to an external Method not implemented
                  global const some_linfo = Core.Compiler.specialize_method(some_method,
                      Tuple{typeof(Base.include), Module, String}, Core.svec())
              end

              g() = override(1.0)
              Test.@test g() === 2.0 # compile this

              const abigfloat_f() = big"12.34"
              const abigfloat_x = big"43.21"
              const abigint_f() = big"123"
              const abigint_x = big"124"

              # issue #31488
              _v31488 = Base.StringVector(2)
              resize!(_v31488, 0)
              const a31488 = fill(String(_v31488), 100)

              const ptr1 = Ptr{UInt8}(1)
              ptr2 = Ptr{UInt8}(1)
              const ptr3 = Ptr{UInt8}(-1)
              const layout1 = Ptr{Int8}[Ptr{Int8}(0), Ptr{Int8}(1), Ptr{Int8}(-1)]
              const layout2 = Any[Ptr{Int8}(0), Ptr{Int16}(1), Ptr{Int32}(-1)]
              const layout3 = collect(x.match for x in eachmatch(r"..", "abcdefghijk"))::Vector{SubString{String}}

              # create a backedge that includes Type{Union{}}, to ensure lookup can handle that
              call_bottom() = show(stdout::IO, Union{})
              Core.Compiler.return_type(call_bottom, ())

              # check that @ccallable works from precompiled modules
              Base.@ccallable Cint f35014(x::Cint) = x+Cint(1)
          end
          """)
    # make sure `sin` didn't have a kwfunc (which would invalidate the attempted test)
    @test !isdefined(typeof(sin).name.mt, :kwsorter)

    # Issue #12623
    @test __precompile__(false) === nothing

    # Issue #21307
    Foo2 = Base.require(Main, Foo2_module)
    @eval $Foo2.override(::Int) = 'a'
    @eval $Foo2.override(::Float32) = 'b'

    Foo = Base.require(Main, Foo_module)
    Base.invokelatest() do # use invokelatest to see the results of loading the compile
        @test Foo.foo(17) == 18
        @test Foo.Bar.bar(17) == 19

        # Issue #21307
        @test Foo.g() === 97.0
        @test Foo.override(1.0e0) == Float64('a')
        @test Foo.override(1.0f0) == 'b'
        @test Foo.override(UInt(1)) == 2

        # Issue #15722
        @test Foo.abigfloat_f()::BigFloat == big"12.34"
        @test (Foo.abigfloat_x::BigFloat + 21) == big"64.21"
        @test Foo.abigint_f()::BigInt == big"123"
        @test Foo.abigint_x::BigInt + 1 == big"125"

        @test Foo.x28297.result === missing

        @test Foo.d29936a === Dict
        @test Foo.d29936b === Dict{K,V} where {V,K}

        @test Foo.x28998[end] == 6

        @test Foo.a31488 == fill("", 100)

        @test Foo.ptr1 === Ptr{UInt8}(1)
        @test Foo.ptr2 === Ptr{UInt8}(0)
        @test Foo.ptr3 === Ptr{UInt8}(-1)
        @test Foo.layout1::Vector{Ptr{Int8}} == Ptr{Int8}[Ptr{Int8}(0), Ptr{Int8}(0), Ptr{Int8}(-1)]
        @test Foo.layout2 == Any[Ptr{Int8}(0), Ptr{Int16}(0), Ptr{Int32}(-1)]
        @test typeof.(Foo.layout2) == [Ptr{Int8}, Ptr{Int16}, Ptr{Int32}]
        @test Foo.layout3 == ["ab", "cd", "ef", "gh", "ij"]
    end

    @eval begin function ccallable_test()
        Base.llvmcall(
        ("""declare i32 @f35014(i32)
            define i32 @entry() {
            0:
                %1 = call i32 @f35014(i32 3)
                ret i32 %1
            }""", "entry"
        ), Cint, Tuple{})
    end
    @test ccallable_test() == 4
    end

    cachedir = joinpath(dir, "compiled", "v$(VERSION.major).$(VERSION.minor)")
    cachedir2 = joinpath(dir2, "compiled", "v$(VERSION.major).$(VERSION.minor)")
    cachefile = joinpath(cachedir, "$Foo_module.ji")
    # use _require_from_serialized to ensure that the test fails if
    # the module doesn't reload from the image:
    @test_warn "@ccallable was already defined for this method name" begin
        @test_logs (:warn, "Replacing module `$Foo_module`") begin
            ms = Base._require_from_serialized(cachefile)
            @test isa(ms, Array{Any,1})
        end
    end

    @test_throws MethodError Foo.foo(17) # world shouldn't be visible yet
    Base.invokelatest() do # use invokelatest to see the results of loading the compile
        @test Foo.foo(17) == 18
        @test Foo.Bar.bar(17) == 19

        # Issue #21307
        @test Foo.g() === 97.0
        @test Foo.override(1.0e0) == Float64('a')
        @test Foo.override(1.0f0) == 'b'
        @test Foo.override(UInt(1)) == 2

        # issue #12284:
        @test string(Base.Docs.doc(Foo.foo)) == "foo function\n"
        @test string(Base.Docs.doc(Foo.Bar.bar)) == "bar function\n"
        @test string(Base.Docs.doc(Foo.Bar)) == "Bar module\n"

        modules, (deps, requires), required_modules = Base.parse_cache_header(cachefile)
        discard_module = mod_fl_mt -> (mod_fl_mt.filename, mod_fl_mt.mtime)
        @test modules == [ Base.PkgId(Foo) => Base.module_build_id(Foo) ]
        @test map(x -> x.filename, deps) == [ Foo_file, joinpath(dir, "foo.jl"), joinpath(dir, "bar.jl") ]
        @test requires == [ Base.PkgId(Foo) => Base.PkgId(string(FooBase_module)),
                            Base.PkgId(Foo) => Base.PkgId(Foo2),
                            Base.PkgId(Foo) => Base.PkgId(Test),
                            Base.PkgId(Foo) => Base.PkgId(string(FooBase_module)) ]
        srctxt = Base.read_dependency_src(cachefile, Foo_file)
        @test !isempty(srctxt) && srctxt == read(Foo_file, String)
        @test_throws ErrorException Base.read_dependency_src(cachefile, "/tmp/nonexistent.txt")
        # dependencies declared with `include_dependency` should not be stored
        @test_throws ErrorException Base.read_dependency_src(cachefile, joinpath(dir, "foo.jl"))

        modules, deps1 = Base.cache_dependencies(cachefile)
        @test Dict(modules) == merge(
            Dict(let m = Base.PkgId(s)
                    m => Base.module_build_id(Base.root_module(m))
                 end for s in
                 [ "Base", "Core", "Main",
                   string(Foo2_module), string(FooBase_module) ]),
            # plus modules included in the system image
            Dict(let m = Base.root_module(Base, s)
                     Base.PkgId(m) => Base.module_build_id(m)
                 end for s in
                [:ArgTools, :Artifacts, :Base64, :CRC32c, :Dates, :DelimitedFiles,
                 :Distributed, :Downloads, :FileWatching, :Future, :InteractiveUtils,
                 :LazyArtifacts, :LibCURL, :LibCURL_jll, :LibGit2, :Libdl, :LinearAlgebra,
                 :Logging, :Markdown, :Mmap, :MozillaCACerts_jll, :NetworkOptions, :Pkg, :Printf,
                 :Profile, :p7zip_jll, :REPL, :Random, :SHA, :Serialization, :SharedArrays, :Sockets,
                 :SparseArrays, :Statistics, :SuiteSparse, :TOML, :Tar, :Test, :UUIDs, :Unicode,
                 :nghttp2_jll]
            ),
        )
        @test discard_module.(deps) == deps1
        modules, (deps, requires), required_modules = Base.parse_cache_header(cachefile; srcfiles_only=true)
        @test map(x -> x.filename, deps) == [Foo_file]

        @test current_task()(0x01, 0x4000, 0x30031234) == 2
        @test sin(0x01, 0x4000, 0x30031234) == 52
        @test sin(0x01, 0x4000, 0x30031234; x = 9142) == 9142
        @test Foo.sinkw === Core.kwfunc(Base.sin)

        @test Foo.NominalValue() == 1
        @test Foo.OrdinalValue() == 1
        @test Foo.NominalValue{Int}() == 2
        @test Foo.OrdinalValue{Int}() == 2
        let T = Vector{Foo.NominalValue{Int}}
            @test isa(T(), T)
        end
        @test Vector{Foo.NominalValue{Int32, Int64}}() == 3
        @test Vector{Foo.NominalValue{UInt, UInt}}() == 4
        @test Vector{Foo.NominalValue{Int, Int}}() == 5
        @test all(i -> Foo.t17809s[i + 1] ===
            Tuple{
                Type{Ptr{Foo.MyType{i}}},
                Ptr{Type{Foo.MyType{i}}},
                Array{Ptr{Foo.MyType{Foo.MyType{:sym}()}}(0), 0},
                Val{Complex{Int}(1, 2)},
                Val{3},
                Val{nothing}},
            0:25)
        some_method = which(Base.include, (Module, String,))
        some_linfo = Core.Compiler.specialize_method(some_method, Tuple{typeof(Base.include), Module, String}, Core.svec())
        @test Foo.some_linfo::Core.MethodInstance === some_linfo

        ft = Base.datatype_fieldtypes
        PV = ft(Foo.Value18343{Some}.body)[1]
        VR = ft(PV)[1].parameters[1]
        @test ft(PV)[1] === Array{VR,1}
        @test pointer_from_objref(ft(PV)[1]) ===
              pointer_from_objref(ft(ft(ft(PV)[1].parameters[1])[1])[1])
        @test PV === ft(ft(PV)[1].parameters[1])[1]
        @test pointer_from_objref(PV) === pointer_from_objref(ft(ft(PV)[1].parameters[1])[1])
    end

    Nest_module = :Nest4b3a94a1a081a8cb
    Nest_file = joinpath(dir, "$Nest_module.jl")
    NestInner_file = joinpath(dir, "$(Nest_module)Inner.jl")
    NestInner2_file = joinpath(dir, "$(Nest_module)Inner2.jl")
    write(Nest_file,
        """
        module $Nest_module
        include("$(escape_string(NestInner_file))")
        end
        """)
    write(NestInner_file,
        """
        module NestInner
        include("$(escape_string(NestInner2_file))")
        end
        """)
    write(NestInner2_file,
        """
        f() = 22
        """)
    Nest = Base.require(Main, Nest_module)
    cachefile = joinpath(cachedir, "$Nest_module.ji")
    modules, (deps, requires), required_modules = Base.parse_cache_header(cachefile)
    @test last(deps).modpath == ["NestInner"]

    UsesB_module = :UsesB4b3a94a1a081a8cb
    B_module     = :UsesB4b3a94a1a081a8cb_B
    UsesB_file = joinpath(dir, "$UsesB_module.jl")
    B_file = joinpath(dir, "$(B_module).jl")
    write(UsesB_file,
        """
        module $UsesB_module
        using $B_module
        end
        """)
    write(B_file,
        """
        module $B_module
        export bfunc
        bfunc() = 33
        end
        """)
    UsesB = Base.require(Main, UsesB_module)
    cachefile = joinpath(cachedir, "$UsesB_module.ji")
    modules, (deps, requires), required_modules = Base.parse_cache_header(cachefile)
    id1, id2 = only(requires)
    @test Base.pkgorigins[id1].cachepath == cachefile
    @test Base.pkgorigins[id2].cachepath == joinpath(cachedir, "$B_module.ji")

    Baz_file = joinpath(dir, "Baz.jl")
    write(Baz_file,
          """
          true && __precompile__(false)
          module Baz
          baz() = 1
          end
          """)

    @test Base.compilecache(Base.PkgId("Baz")) == Base.PrecompilableError() # due to __precompile__(false)
    @eval using Baz
    @test Base.invokelatest(Baz.baz) == 1

    # Issue #12720
    FooBar1_file = joinpath(dir, "FooBar1.jl")
    write(FooBar1_file,
          """
          module FooBar1
              using FooBar
          end
          """)
    sleep(2) # give FooBar and FooBar1 different timestamps, in reverse order too
    FooBar_file = joinpath(dir, "FooBar.jl")
    write(FooBar_file,
          """
          module FooBar
          end
          """)

    cachefile = Base.compilecache(Base.PkgId("FooBar"))
    empty_prefs_hash = Base.get_preferences_hash(nothing, String[])
    @test cachefile == Base.compilecache_path(Base.PkgId("FooBar"), empty_prefs_hash)
    @test isfile(joinpath(cachedir, "FooBar.ji"))
    @test Base.stale_cachefile(FooBar_file, joinpath(cachedir, "FooBar.ji")) isa Vector
    @test !isdefined(Main, :FooBar)
    @test !isdefined(Main, :FooBar1)

    relFooBar_file = joinpath(dir, "subfolder", "..", "FooBar.jl")
    @test Base.stale_cachefile(relFooBar_file, joinpath(cachedir, "FooBar.ji")) isa (Sys.iswindows() ? Vector : Bool) # `..` is not a symlink on Windows
    mkdir(joinpath(dir, "subfolder"))
    @test Base.stale_cachefile(relFooBar_file, joinpath(cachedir, "FooBar.ji")) isa Vector

    @eval using FooBar
    fb_uuid = Base.module_build_id(FooBar)
    sleep(2); touch(FooBar_file)
    insert!(DEPOT_PATH, 1, dir2)
    @test Base.stale_cachefile(FooBar_file, joinpath(cachedir, "FooBar.ji")) === true
    @eval using FooBar1
    @test !isfile(joinpath(cachedir2, "FooBar.ji"))
    @test !isfile(joinpath(cachedir, "FooBar1.ji"))
    @test isfile(joinpath(cachedir2, "FooBar1.ji"))
    @test Base.stale_cachefile(FooBar_file, joinpath(cachedir, "FooBar.ji")) === true
    @test Base.stale_cachefile(FooBar1_file, joinpath(cachedir2, "FooBar1.ji")) isa Vector
    @test fb_uuid == Base.module_build_id(FooBar)
    fb_uuid1 = Base.module_build_id(FooBar1)
    @test fb_uuid != fb_uuid1

    # test checksum
    open(joinpath(cachedir2, "FooBar1.ji"), "a") do f
        write(f, 0x076cac96) # append 4 random bytes
    end
    @test Base.stale_cachefile(FooBar1_file, joinpath(cachedir2, "FooBar1.ji")) === true

    # test behavior of precompile modules that throw errors
    FooBar2_file = joinpath(dir, "FooBar2.jl")
    write(FooBar2_file,
          """
          module FooBar2
          error("break me")
          end
          """)
    @test_warn r"LoadError: break me\nStacktrace:\n \[1\] [\e01m\[]*error" try
            Base.require(Main, :FooBar2)
            error("the \"break me\" test failed")
        catch exc
            isa(exc, ErrorException) || rethrow()
            occursin("ERROR: LoadError: break me", exc.msg) && rethrow()
        end

    # Test that trying to eval into closed modules during precompilation is an error
    FooBar3_file = joinpath(dir, "FooBar3.jl")
    FooBar3_inc = joinpath(dir, "FooBar3_inc.jl")
    write(FooBar3_inc, "x=1\n")
    for code in ["Core.eval(Base, :(x=1))", "Base.include(Base, \"FooBar3_inc.jl\")"]
        write(FooBar3_file, code)
        @test_warn "Evaluation into the closed module `Base` breaks incremental compilation" try
                Base.require(Main, :FooBar3)
            catch exc
                isa(exc, ErrorException) || rethrow()
            end
    end

    # Test transitive dependency for #21266
    FooBarT_file = joinpath(dir, "FooBarT.jl")
    write(FooBarT_file,
          """
          module FooBarT
          end
          """)
    FooBarT1_file = joinpath(dir, "FooBarT1.jl")
    write(FooBarT1_file,
          """
          module FooBarT1
              using FooBarT
          end
          """)
    FooBarT2_file = joinpath(dir, "FooBarT2.jl")
    write(FooBarT2_file,
          """
          module FooBarT2
              using FooBarT1
          end
          """)
    Base.compilecache(Base.PkgId("FooBarT2"))
    write(FooBarT1_file,
          """
          module FooBarT1
          end
          """)
    rm(FooBarT_file)
    @test Base.stale_cachefile(FooBarT2_file, joinpath(cachedir2, "FooBarT2.ji")) === true
    @test Base.require(Main, :FooBarT2) isa Module
end
end

# method root provenance
# setindex!(::Dict{K,V}, ::Any, ::K) adds both compression and codegen roots
precompile_test_harness("code caching") do dir
    Bid = rootid(Base)
    Cache_module = :Cacheb8321416e8a3e2f1
    write(joinpath(dir, "$Cache_module.jl"),
          """
          module $Cache_module
              struct X end
              struct X2 end
              @noinline function f(d)
                  @noinline
                  d[X()] = nothing
              end
              @noinline fpush(dest) = push!(dest, X())
              function callboth()
                  f(Dict{X,Any}())
                  fpush(X[])
                  nothing
              end
              precompile(callboth, ())
          end
          """)
    Base.compilecache(Base.PkgId(string(Cache_module)))
    @eval using $Cache_module
    M = getfield(@__MODULE__, Cache_module)
    Mid = rootid(M)
    for name in (:f, :fpush, :callboth)
        func = getfield(M, name)
        m = only(collect(methods(func)))
        @test all(i -> root_provenance(m, i) == Mid, 1:length(m.roots))
    end
    m = which(setindex!, (Dict{M.X,Any}, Any, M.X))
    @test_broken M.X ∈ m.roots               # requires caching external compilation results
    Base.invokelatest() do
        Dict{M.X2,Any}()[M.X2()] = nothing
    end
    @test M.X2 ∈ m.roots
    groups = group_roots(m)
    @test_broken M.X ∈ groups[Mid]           # requires caching external compilation results
    @test M.X2 ∈ groups[rootid(@__MODULE__)]
    @test !isempty(groups[Bid])
    # PkgA loads PkgB, and both add roots to the same method (both before and after loading B)
    Cache_module2 = :Cachea1544c83560f0c99
    write(joinpath(dir, "$Cache_module2.jl"),
          """
          module $Cache_module2
              struct Y end
              @noinline f(dest) = push!(dest, Y())
              callf() = f(Y[])
              callf()
              using $(Cache_module)
              struct Z end
              @noinline g(dest) = push!(dest, Z())
              callg() = g(Z[])
              callg()
          end
          """)
    Base.compilecache(Base.PkgId(string(Cache_module2)))
    @eval using $Cache_module2
    M2 = getfield(@__MODULE__, Cache_module2)
    M2id = rootid(M2)
    dest = []
    Base.invokelatest() do  # use invokelatest to see the results of loading the compile
        M2.f(dest)
        M.fpush(dest)
        M2.g(dest)
        @test dest == [M2.Y(), M.X(), M2.Z()]
        @test M2.callf() == [M2.Y()]
        @test M2.callg() == [M2.Z()]
        @test M.fpush(M.X[]) == [M.X()]
    end
    mT = which(push!, (Vector{T} where T, Any))
    groups = group_roots(mT)
    # all below require caching external CodeInstances
    @test_broken M2.Y ∈ groups[M2id]
    @test_broken M2.Z ∈ groups[M2id]
    @test_broken M.X ∈ groups[Mid]
    @test_broken M.X ∉ groups[M2id]
end

# test --compiled-modules=no command line option
precompile_test_harness("--compiled-modules=no") do dir
    Time_module = :Time4b3a94a1a081a8cb
    write(joinpath(dir, "$Time_module.jl"),
          """
          module $Time_module
              time = Base.time()
          end
          """)
    Base.compilecache(Base.PkgId("Time4b3a94a1a081a8cb"))
    exename = `$(Base.julia_cmd()) --compiled-modules=yes --startup-file=no`
    testcode = """
        insert!(LOAD_PATH, 1, $(repr(dir)))
        insert!(DEPOT_PATH, 1, $(repr(dir)))
        using $Time_module
        getfield($Time_module, :time)
    """

    t1_yes = readchomp(`$exename --compiled-modules=yes -E $(testcode)`)
    t2_yes = readchomp(`$exename --compiled-modules=yes -E $(testcode)`)
    @test t1_yes == t2_yes

    t1_no = readchomp(`$exename --compiled-modules=no -E $(testcode)`)
    t2_no = readchomp(`$exename --compiled-modules=no -E $(testcode)`)
    @test t1_no != t2_no
    @test parse(Float64, t1_no) < parse(Float64, t2_no)
end

# test loading a package with conflicting namespace
precompile_test_harness("conflicting namespaces") do dir
    Test_module = :Test6c92f26
    write(joinpath(dir, "Iterators.jl"),
          """
          module Iterators
          end
          """)
    write(joinpath(dir, "$Test_module.jl"),
          """
          module $Test_module
               import Iterators # FIXME: use `using`
          end
          """)
    testcode = """
        insert!(LOAD_PATH, 1, $(repr(dir)))
        insert!(DEPOT_PATH, 1, $(repr(dir)))
        using $Test_module
        println(stderr, $Test_module.Iterators)
    """

    exename = `$(Base.julia_cmd()) --startup-file=no`
    let fname = tempname()
        try
            for i = 1:2
                @test readchomp(pipeline(`$exename -E $(testcode)`, stderr=fname)) == "nothing"
                @test read(fname, String) == "Iterators\n"
            end
        finally
            rm(fname, force=true)
        end
    end
end

precompile_test_harness("package_callbacks") do dir
    loaded_modules = Channel{Symbol}(32)
    callback = (mod::Base.PkgId) -> put!(loaded_modules, Symbol(mod.name))
    push!(Base.package_callbacks, callback)
    try
        Test1_module = :Teste4095a81
        Test2_module = :Teste4095a82
        Test3_module = :Teste4095a83

        write(joinpath(dir, "$(Test1_module).jl"),
              """
              module $(Test1_module)
              end
              """)
        Base.compilecache(Base.PkgId("$(Test1_module)"))

        write(joinpath(dir, "$(Test2_module).jl"),
              """
              module $(Test2_module)
                  using $(Test1_module)
              end
              """)
        Base.compilecache(Base.PkgId("$(Test2_module)"))

        @test !Base.isbindingresolved(Main, Test2_module)
        Base.require(Main, Test2_module)
        @test take!(loaded_modules) == Test1_module
        @test take!(loaded_modules) == Test2_module
        write(joinpath(dir, "$(Test3_module).jl"),
              """
              module $(Test3_module)
                  using $(Test3_module)
              end
              """)
        Base.require(Main, Test3_module)
        @test take!(loaded_modules) == Test3_module
    finally
        pop!(Base.package_callbacks)
    end
end

# Issue #19960
(f -> f())() do # wrap in function scope, so we can test world errors
    test_workers = addprocs(1)
    push!(test_workers, myid())
    save_cwd = pwd()
    temp_path = mktempdir()
    try
        cd(temp_path)
        load_path = mktempdir(temp_path)
        load_cache_path = mktempdir(temp_path)

        ModuleA = :Issue19960A
        ModuleB = :Issue19960B

        write(joinpath(load_path, "$ModuleA.jl"),
            """
            module $ModuleA
                import Distributed: myid
                export f
                f() = myid()
            end
            """)

        write(joinpath(load_path, "$ModuleB.jl"),
            """
            module $ModuleB
                using $ModuleA
                export g
                g() = f()
            end
            """)

        @everywhere test_workers begin
            pushfirst!(LOAD_PATH, $load_path)
            pushfirst!(DEPOT_PATH, $load_cache_path)
        end
        try
            @eval using $ModuleB
            uuid = Base.module_build_id(Base.root_module(Main, ModuleB))
            for wid in test_workers
                @test Distributed.remotecall_eval(Main, wid, quote
                        Base.module_build_id(Base.root_module(Main, $(QuoteNode(ModuleB))))
                    end) == uuid
                if wid != myid() # avoid world-age errors on the local proc
                    @test remotecall_fetch(g, wid) == wid
                end
            end
        finally
            @everywhere test_workers begin
                popfirst!(LOAD_PATH)
                popfirst!(DEPOT_PATH)
            end
        end
    finally
        cd(save_cwd)
        rm(temp_path, recursive=true)
        pop!(test_workers) # remove myid
        rmprocs(test_workers)
    end
end

# Ensure that module-loading plays nicely with Base.delete_method
# wrapped in function scope, so we can test world errors
precompile_test_harness("delete_method") do dir
    A_module = :Aedb164bd3a126418
    B_module = :Bedb164bd3a126418
    A_file = joinpath(dir, "$A_module.jl")
    B_file = joinpath(dir, "$B_module.jl")
    write(A_file,
          """
          module $A_module

          export apc, anopc

          apc(::Int, ::Int) = 1
          apc(::Any, ::Any) = 2

          anopc(::Int, ::Int) = 1
          anopc(::Any, ::Any) = 2

          end
          """)
    write(B_file,
          """
          module $B_module

          using $A_module

          bpc(x) = apc(x, x)
          bnopc(x) = anopc(x, x)

          precompile(bpc, (Int,))
          precompile(bpc, (Float64,))

          end
          """)
    A = Base.require(Main, A_module)
    for mths in (collect(methods(A.apc)), collect(methods(A.anopc)))
        Base.delete_method(mths[1])
    end
    B = Base.require(Main, B_module)
    @test Base.invokelatest(B.bpc, 1) == Base.invokelatest(B.bpc, 1.0) == 2
    @test Base.invokelatest(B.bnopc, 1) == Base.invokelatest(B.bnopc, 1.0) == 2
end

precompile_test_harness("Issues #19030 and #25279") do load_path
    ModuleA = :Issue19030
    write(joinpath(load_path, "$ModuleA.jl"),
        """
        module $ModuleA
            __init__() = push!(Base.package_callbacks, sym->nothing)
        end
        """)
    l0 = length(Base.package_callbacks)
    @eval using $ModuleA
    @test length(Base.package_callbacks) == l0 + 1
end

precompile_test_harness("Issue #25604") do load_path
    write(joinpath(load_path, "A25604.jl"),
        """
        module A25604
        using B25604
        using C25604
        end
        """)
    write(joinpath(load_path, "B25604.jl"),
        """
        module B25604
        end
        """)
    write(joinpath(load_path, "C25604.jl"),
        """
        module C25604
        using B25604
        end
        """)
    Base.compilecache(Base.PkgId("A25604"))
    @test_nowarn @eval using A25604
end

precompile_test_harness("Issue #26028") do load_path
    write(joinpath(load_path, "Foo26028.jl"),
        """
        module Foo26028
        module Bar26028
            x = 0
        end
        function __init__()
            include(joinpath(@__DIR__, "Baz26028.jl"))
        end
        end
        """)
    write(joinpath(load_path, "Baz26028.jl"),
        """
        module Baz26028
        import Foo26028.Bar26028.x
        end
        """)
    Base.compilecache(Base.PkgId("Foo26028"))
    @test_nowarn @eval using Foo26028
end

precompile_test_harness("Issue #29936") do load_path
    write(joinpath(load_path, "Foo29936.jl"),
          """
          module Foo29936
          const global m = Val{nothing}()
          const global h = Val{:hey}()
          wab = [("a", m), ("b", h),]
          end
          """)
    @eval using Foo29936
    @test [("Plan", Foo29936.m), ("Plan", Foo29936.h),] isa Vector{Tuple{String,Val}}
end

precompile_test_harness("Issue #25971") do load_path
    sourcefile = joinpath(load_path, "Foo25971.jl")
    write(sourcefile, "module Foo25971 end")
    chmod(sourcefile, 0o666)
    cachefile = Base.compilecache(Base.PkgId("Foo25971"))
    @test filemode(sourcefile) == filemode(cachefile)
    chmod(sourcefile, 0o600)
    cachefile = Base.compilecache(Base.PkgId("Foo25971"))
    @test filemode(sourcefile) == filemode(cachefile)
    chmod(sourcefile, 0o444)
    cachefile = Base.compilecache(Base.PkgId("Foo25971"))
    # Check writable
    @test touch(cachefile) == cachefile
end

precompile_test_harness("Issue #38312") do load_path
    TheType = """Array{Ref{Val{1}}, 1}"""
    write(joinpath(load_path, "Foo38312.jl"),
        """
        module Foo38312
        const TheType = $TheType
        end
        """)
    write(joinpath(load_path, "Bar38312.jl"),
        """
        module Bar38312
        const TheType = $TheType
        end
        """)
    Base.compilecache(Base.PkgId("Foo38312"))
    Base.compilecache(Base.PkgId("Bar38312"))
    @test pointer_from_objref((@eval (using Foo38312; Foo38312)).TheType) ===
          pointer_from_objref(eval(Meta.parse(TheType))) ===
          pointer_from_objref((@eval (using Bar38312; Bar38312)).TheType)
end

precompile_test_harness("Opaque Closure") do load_path
    write(joinpath(load_path, "OCPrecompile.jl"),
        """
        module OCPrecompile
        using Base.Experimental: @opaque
        f(x) = @opaque y->x+y
        end
        """)
    Base.compilecache(Base.PkgId("OCPrecompile"))
    f = (@eval (using OCPrecompile; OCPrecompile)).f
    @test Base.invokelatest(f, 1)(2) == 3
end

# issue #39405
precompile_test_harness("Renamed Imports") do load_path
    write(joinpath(load_path, "RenameImports.jl"),
          """
          module RenameImports
          import Base.Experimental as ex
          test() = ex
          end
          """)
    Base.compilecache(Base.PkgId("RenameImports"))
    @test (@eval (using RenameImports; RenameImports.test())) isa Module
end

# issue #41872 (example from #38983)
precompile_test_harness("No external edges") do load_path
    write(joinpath(load_path, "NoExternalEdges.jl"),
          """
          module NoExternalEdges
          bar(x::Int) = hcat(rand())
          @inline bar() = hcat(rand())
          bar(x::Float64) = bar()
          foo1() = bar(1)
          foo2() = bar(1.0)
          foo3() = bar()
          foo4() = hcat(rand())
          precompile(foo1, ())
          precompile(foo2, ())
          precompile(foo3, ())
          precompile(foo4, ())
          end
          """)
    Base.compilecache(Base.PkgId("NoExternalEdges"))
    @eval begin
        using NoExternalEdges
        @test only(methods(NoExternalEdges.foo1)).specializations[1].cache.max_world != 0
        @test only(methods(NoExternalEdges.foo2)).specializations[1].cache.max_world != 0
        @test only(methods(NoExternalEdges.foo3)).specializations[1].cache.max_world != 0
        @test only(methods(NoExternalEdges.foo4)).specializations[1].cache.max_world != 0
    end
end

@testset "issue 38149" begin
    M = Module()
    @eval M begin
        @nospecialize
        f(x, y) = x + y
        f(x::Int, y) = 2x + y
    end
    precompile(M.f, (Int, Any))
    precompile(M.f, (AbstractFloat, Any))
    mis = map(methods(M.f)) do m
        m.specializations[1]
    end
    @test any(mi -> mi.specTypes.parameters[2] === Any, mis)
    @test all(mi -> isa(mi.cache, Core.CodeInstance), mis)
end

max_stars_repo_name: adkabo/julia

max_issues_repo_name: adkabo/julia

max_forks_repo_name: adkabo/julia


text:
######
##FEM Heat dx Convergence Tests
######
using FiniteElementDiffEq, DiffEqDevTools

#Travis CI Test Setting
#Not good plots, but quick for unit tests
cs = cs_femheat_moving_dx

alg=FEMDiffEqHeatEuler(); println(alg) #Unstable due to μ
sim = test_convergence(cs,alg)

alg=FEMDiffEqHeatImplicitEuler(); println(alg)
sim2 = test_convergence(cs,alg)

alg=FEMDiffEqHeatCrankNicholson(); println(alg) #Bound by spatial discretization error at low dt, decrease dx for full convergence
cs = cs_femheat_moving_faster_dt
sim3 = test_convergence(cs,alg)

#plot(plot(sim),plot(sim2),plot(sim3),layout=@layout([a b c]),size=(1200,400))

#Returns true if all converge approximately dx^2
minimum([sim.𝒪est[:L2],sim2.𝒪est[:L2],sim3.𝒪est[:L2]] - 2 .<.1)

ext: jl


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(HDF5.__init__)})
    precompile(Tuple{typeof(HDF5.blosc_filter), UInt32, UInt64, Ptr{UInt32}, UInt64, Ptr{UInt64}, Ptr{Ptr{Nothing}}})
    precompile(Tuple{typeof(HDF5.blosc_set_local), Int64, Int64, Int64})
end


text:
# Copied from https://github.com/gridap/GridapGmsh.jl/blob/master/deps/build.jl
# to use the new officially registered Gmsh JLL

gmsh_root = nothing

if haskey(ENV,"GMSHROOT") && !isempty(ENV["GMSHROOT"])
  gmsh_root = ENV["GMSHROOT"]
  @info """
  Using the gmsh installation found via the GMSHROOT environment variable.
  GMSHROOT=$gmsh_root"""
else
  using gmsh_jll
  GMSH_FOUND = gmsh_jll.is_available()
  if GMSH_FOUND
    gmsh_root = gmsh_jll.artifact_dir
    @info """
    Using the gmsh installation installed via BinaryBuilder."""
  end
end

GMSH_FOUND = false
if gmsh_root !== nothing
  gmsh_bin = (Sys.iswindows() ? joinpath(gmsh_root,"bin","gmsh.exe") :
              joinpath(gmsh_root,"bin","gmsh"))
  if isfile(gmsh_bin)
    GMSH_FOUND = true
    @info "gmsh binary found in $gmsh_bin"
  else
    @warn """
    gmsh binary not found in $gmsh_bin
    Make sure that GMSHROOT points to a correct gmsh installation"""
  end
  gmsh_jl = joinpath(gmsh_root,"lib","gmsh.jl")
  if isfile(gmsh_jl)
    GMSH_FOUND = true
    @info "Using the gmsh Julia API found in $gmsh_jl"
  else
    @warn """
    gmsh Julia API not found in $gmsh_jl
    Make sure that you have installed the full gmsh SDK"""
  end
end

open("deps.jl","w") do f
  println(f, "# This file is automatically generated")
  println(f, "# Do not edit")
  println(f)
  println(f, :(const GMSH_FOUND = $GMSH_FOUND))
  println(f, :(const gmsh_jl    = $gmsh_jl))
end


text:
##############################################################################
##
## Solution method: minimize by sparse least square optimization
##
##############################################################################

##############################################################################
##
## Factor Model (no regressors)
##
##############################################################################

function StatsBase.fit!(t::Union{Type{Val{:levenberg_marquardt}}, Type{Val{:dogleg}}},
                         fp::FactorModel,
                         fs::FactorSolution;
                         maxiter::Integer = 10_000,
                         tol::Real = 1e-9,
                         lambda::Number = 0.0)
    # initialize
    iter = 0
    converged = true
    fullrank = rank(fp)
    fp = FactorModel(copy(fp.y), fp.sqrtw, fp.idrefs, fp.timerefs, 1)
    N = size(fs.idpool, 1)
    T = size(fs.timepool, 1)
    fg = FactorGradient(fp,
                        FactorSolution(Array{Float64}(undef,N), Array{Float64}(undef, T)),
                        FactorSolution(Array{Float64}(undef,N), Array{Float64}(undef, T))
                       )
    nls = LeastSquaresOptim.LeastSquaresProblem(
                view(fs, :, 1),
                similar(fp.y),
                (y, x) -> f!(y, x, fp),
                fg,
                g!)
    if t == Val{:levenberg_marquardt}
        optimizer = LeastSquaresOptim.LevenbergMarquardt(LeastSquaresOptim.LSMR())
    else
        optimizer = LeastSquaresOptim.Dogleg(LeastSquaresOptim.LSMR())
    end

    full = LeastSquaresOptim.LeastSquaresProblemAllocated(nls, optimizer)
    for r in 1:fullrank
        fsr = view(fs, :, r)
        full.x = fsr
        result = LeastSquaresOptim.optimize!(full,
            x_tol = 1e-32, g_tol = 1e-32, f_tol = tol,  iterations = maxiter)
        iter += result.mul_calls
        converged = result.converged && converged
        if r < fullrank
            rescale!(fsr)
            subtract_factor!(fp, fsr)
        end
    end
    # rescale factors and loadings so that factors' * factors = Id
    return rescale(fs), iter, converged
end

##############################################################################
##
## Interactive FixedEffectModel
##
##############################################################################

function StatsBase.fit!(t::Union{Type{Val{:levenberg_marquardt}}, Type{Val{:dogleg}}},
                         fp::InteractiveFixedEffectsModel{Rank},
                         fs::InteractiveFixedEffectsSolution{Rank};
                         maxiter::Integer = 10_000,
                         tol::Real = 1e-9,
                         lambda::Number = 0.0) where {Rank}
    scaleb = vec(sum(abs2, fp.X, dims = 1))
    fsT = InteractiveFixedEffectsSolutionT(fs.b, fs.idpool', fs.timepool')
    fg = InteractiveFixedEffectsGradientT(fp,
                similar(fsT),
                InteractiveFixedEffectsSolutionT(scaleb, similar(fsT.idpool), similar(fsT.timepool)))
    nls = LeastSquaresOptim.LeastSquaresProblem(fsT, similar(fp.y), (y, x) -> f!(y, x, fp), fg, g!)
    if t == Val{:levenberg_marquardt}
        optimizer = LeastSquaresOptim.LevenbergMarquardt(LeastSquaresOptim.LSMR())
    else
        optimizer = LeastSquaresOptim.Dogleg(LeastSquaresOptim.LSMR())
    end
    full = LeastSquaresOptim.LeastSquaresProblemAllocated(nls, optimizer)
    temp = similar(fp.y)
    result = LeastSquaresOptim.optimize!(full;
            x_tol = 1e-32, g_tol = 1e-32, f_tol = tol,  iterations = maxiter)
    # rescale factors and loadings so that factors' * factors = Id
    fs = InteractiveFixedEffectsSolution(fsT.b, fsT.idpool', fsT.timepool')
    return rescale(fs), result.mul_calls, result.converged
end

##############################################################################
##
## Methods used in optimize! in LeastSquares
##
##############################################################################

function LinearAlgebra.dot(fs1::AbstractArray{T}, fs2::AbstractArray{T})   where {T}
    out = zero(typeof(one(T) * one(T)))
    @inbounds @simd for i in eachindex(fs1)
        out += fs1[i] * fs2[i]
    end
    return out
end

function LinearAlgebra.dot(x::AbstractArray{T}, y::AbstractArray{T}, w::AbstractArray{T}) where {T}
    out = zero(typeof(one(T) * one(T)))
    @inbounds @simd for i in eachindex(x)
        out += w[i] * x[i] * y[i]
    end
    return out
end

for t in (FactorSolution, InteractiveFixedEffectsSolution, InteractiveFixedEffectsSolutionT, HalfInteractiveFixedEffectsSolution)
    vars = [:(fill!(x.$field, α)) for field in fieldnames(t)]
    expr = Expr(:block, vars...)
    @eval begin
        function Base.fill!(x::$t, α::Number)
             $expr
            return x
        end
    end

    vars = [:(rmul!(x.$field, α)) for field in fieldnames(t)]
    expr = Expr(:block, vars...)
    @eval begin
        function LinearAlgebra.rmul!(x::$t, α::Number)
            $expr
            return x
        end
    end

    vars = [:(clamp!(x.$field, α, β)) for field in fieldnames(t)]
    expr = Expr(:block, vars...)
    @eval begin
        function Base.clamp!(x::$t, α::Number, β::Number)
            $expr
            return x
        end
    end

    vars = [:(copyto!(x1.$field, x2.$field)) for field in fieldnames(t)]
    expr = Expr(:block, vars...)
    @eval begin
        function Base.copyto!(x1::$t, x2::$t)
            $expr
            return x1
        end
    end

    vars = [:(axpy!(α, x1.$field, x2.$field)) for field in fieldnames(t)]
    expr = Expr(:block, vars...)
    @eval begin
        function LinearAlgebra.axpy!(α::Number, x1::$t, x2::$t)
            $expr
            return x2
        end
    end

    vars = [:(map!(f, x1.$field, map(x -> x.$field, x2)...)) for field in fieldnames(t)]
    expr = Expr(:block, vars...)
    @eval begin
        function Base.map!(f, x1::$t,  x2::$t...)
            $expr
            return x1
        end
    end
    vars = [:(dot(x1.$field, x2.$field)) for field in fieldnames(t)]
    expr = Expr(:call, :+, vars...)
    @eval begin
        function LinearAlgebra.dot(x1::$t, x2::$t)
            $expr
        end
    end

    vars = [:(dot(x1.$field, x2.$field, w.$field)) for field in fieldnames(t)]
    expr = Expr(:call, :+, vars...)
    @eval begin
        function LinearAlgebra.dot(x1::$t, x2::$t, w::$t)
            $expr
        end
    end

    vars = [:(length(x.$field)) for field in fieldnames(t)]
    expr = Expr(:call, :+, vars...)
    @eval begin
        function Base.length(x::$t)
            $expr
        end
    end

    vars = [:(similar(x.$field)) for field in fieldnames(t)]
    expr = Expr(:call, t, vars...)
    @eval begin
        function Base.similar(x::$t)
            $expr
        end
    end

    vars = [:(x.$field) for field in fieldnames(t)]
    expr = Expr(:call, Iterators.flatten, Expr(:call, tuple, vars...))
    @eval begin
        Base.iterate(x::$t) = iterate($expr)
        Base.iterate(x::$t, state) = iterate($expr, state)
    end
end
LinearAlgebra.norm(x::AbstractFactorSolution) = sqrt(sum(abs2, x))

Base.eltype(fg::AbstractFactorSolution) = Float64

function safe_rmul!(x, β)
    if β != 1
        β == 0 ? fill!(x, zero(eltype(x))) : rmul!(x, β)
    end
end

##############################################################################
##
## Factor Gradient
##
##############################################################################

abstract type AbstractFactorGradient{Rank} end

struct FactorGradient{Rank, W, Rid, Rtime, Tid, Ttime, sTid, sTtime} <: AbstractFactorGradient{Rank}
    fp::FactorModel{Rank, W, Rid, Rtime}
    fs::FactorSolution{Rank, Tid, Ttime}
    scalefs::FactorSolution{Rank, sTid, sTtime}
end

function Base.size(fg::FactorGradient, i::Integer)
    if i == 1
        length(fg.fp.y)
    elseif i == 2
      length(fg.fs.idpool) + length(fg.fs.timepool)
    end
end

struct InteractiveFixedEffectsGradientT{Rank, W, Rid, Rtime, Tb, Tid, Ttime, sTb, sTid, sTtime} <: AbstractFactorGradient{Rank}
    fp::InteractiveFixedEffectsModel{Rank, W, Rid, Rtime}
    fs::InteractiveFixedEffectsSolutionT{Rank, Tb, Tid, Ttime}
    scalefs::InteractiveFixedEffectsSolutionT{Rank, sTb, sTid, sTtime}
end
LinearAlgebra.rank(f::AbstractFactorGradient{Rank}) where {Rank} = Rank

function Base.size(fg::InteractiveFixedEffectsGradientT, i::Integer)
    if i == 1
        length(fg.fp.y)
    elseif i == 2
      length(fg.fs.b) + length(fg.fs.idpool) + length(fg.fs.timepool)
    end
end

Base.size(fg::AbstractFactorGradient) = (size(fg, 1), size(fg, 2))
Base.eltype(fg::AbstractFactorGradient) = Float64
LinearAlgebra.adjoint(fg::AbstractFactorGradient) = Adjoint(fg)
LeastSquaresOptim.colsumabs2!(fs::AbstractFactorSolution, fg::AbstractFactorGradient) = copyto!(fs, fg.scalefs)

@generated function LinearAlgebra.mul!(y::AbstractVector{Float64}, fg::AbstractFactorGradient{Rank}, fs::AbstractFactorSolution, α::Number, β::Number) where {Rank}
    if Rank == 1
        ex = quote
            out += (fg.fs.idpool[idi] * fs.timepool[timei]
                + fg.fs.timepool[timei] * fs.idpool[idi])
        end
    else
        ex = quote
            Base.Cartesian.@nexprs $Rank r -> begin
                out += (fg.fs.idpool[r, idi] * fs.timepool[r, timei]
                         + fg.fs.timepool[r, timei] * fs.idpool[r, idi])
            end
        end
    end
    quote
        mα = convert(Float64, -α)
        mul!_X(y, fg.fp, fs, mα, β)
        @fastmath @inbounds @simd for i in 1:length(y)
            timei = fg.fp.timerefs[i]
            idi = fg.fp.idrefs[i]
            out = 0.0
            $ex
            y[i] += mα * fg.fp.sqrtw[i] * out
        end
        return y
    end
end

function mul!_X(y::AbstractVector{Float64}, fp::InteractiveFixedEffectsModel, fs::InteractiveFixedEffectsSolutionT, α::Number, β::Number)
    LinearAlgebra.BLAS.gemm!('N', 'N', α, fp.X, fs.b, convert(Float64, β), y)
end

function mul!_X(y::AbstractVector{Float64}, ::FactorModel, ::FactorSolution, ::Number, β::Number)
    safe_rmul!(y, β)
end

@generated function LinearAlgebra.mul!(fs::AbstractFactorSolution{Rank}, Cfg::Adjoint{T, U}, y::AbstractVector{Float64}, α::Number, β::Number) where {Rank, T, U <: AbstractFactorGradient}
    if Rank == 1
        ex = quote
           fs.idpool[idi] += sqrtwi * fg.fs.timepool[timei]
           fs.timepool[timei] += sqrtwi * fg.fs.idpool[idi]
        end
    else
        ex = quote
            Base.Cartesian.@nexprs $Rank r -> begin
                fs.idpool[r, idi] += sqrtwi * fg.fs.timepool[r, timei]
                fs.timepool[r, timei] += sqrtwi * fg.fs.idpool[r, idi]
            end
        end
    end
    quote
        fg = adjoint(Cfg)
        mα = convert(Float64, -α)
        mul!_X(fs, fg.fp, y, mα, β)
        @inbounds @simd for i in 1:length(y)
            sqrtwi = mα * y[i] * fg.fp.sqrtw[i]
            idi = fg.fp.idrefs[i]
            timei = fg.fp.timerefs[i]
            $ex
        end
        return fs
    end
end

function mul!_X(fs::InteractiveFixedEffectsSolutionT, fp::InteractiveFixedEffectsModel, y::AbstractVector{Float64}, α::Number, β::Number)
    safe_rmul!(fs, β)
    @inbounds @simd for k in 1:length(fs.b)
        out = zero(Float64)
        for i in 1:length(y)
           out += y[i] * fp.X[i, k]
        end
        fs.b[k] += α * out
    end
end

function mul!_X(fs::FactorSolution, ::FactorModel, ::AbstractVector{Float64}, ::Number, β::Number, ) safe_rmul!(fs, β)
end

##############################################################################
##
## f! for LeastSquaresOptim
##
##############################################################################

@generated function f!(out::Vector{Float64}, fs::AbstractFactorSolution{Rank}, fp::AbstractFactorModel{Rank}) where {Rank}
    if Rank == 1
        ex = quote
            out[i] -= sqrtwi * fs.idpool[idi] * fs.timepool[timei]
        end
    else
        ex = quote
           Base.Cartesian.@nexprs $Rank r -> begin
               out[i] -= sqrtwi * fs.idpool[r, idi] * fs.timepool[r, timei]
           end
        end
    end
    quote
        copyto!(out, fp.y)
        mul!_X(out, fp, fs, -1.0, 1.0)
        @fastmath @inbounds @simd for i in 1:length(out)
            sqrtwi = fp.sqrtw[i]
            idi = fp.idrefs[i]
            timei = fp.timerefs[i]
            $ex
        end
        return out
    end
end

##############################################################################
##
## g! for LeastSquaresOptim
##
##############################################################################

@generated function g!(fg::AbstractFactorGradient{Rank}, fs::AbstractFactorSolution{Rank}) where {Rank}
    if Rank == 1
        ex = quote
            fg.scalefs.idpool[idi] += abs2(sqrtwi * fg.fs.timepool[timei])
            fg.scalefs.timepool[timei] += abs2(sqrtwi * fg.fs.idpool[idi])
        end
    else
        ex = quote
            Base.Cartesian.@nexprs $Rank r -> begin
                fg.scalefs.idpool[r, idi] += abs2(sqrtwi * fg.fs.timepool[r, timei])
                fg.scalefs.timepool[r, timei] += abs2(sqrtwi * fg.fs.idpool[r, idi])
            end
        end
    end
    quote
        copyto!(fg.fs, fs)
        fill!(fg.scalefs.idpool, zero(Float64))
        fill!(fg.scalefs.timepool, zero(Float64))
        @inbounds @simd for i in 1:length(fg.fp.y)
            sqrtwi = fg.fp.sqrtw[i]
            idi = fg.fp.idrefs[i]
            timei = fg.fp.timerefs[i]
            $ex
        end
    end
end


text:
module MarkerFile

using ..LIKWID: LibLikwid

"""
Reads in the result file of an application run instrumented by the LIKWID Marker API.

*Note:* julia must have been started under `likwid-perfctr ... -m`.
"""
function read(fp::AbstractString)
    ret = LibLikwid.perfmon_readMarkerFile(fp)
    return ret ≥ 0
end

"""
Return the number of regions in an application run.
"""
numregions() = Int(LibLikwid.perfmon_getNumberOfRegions())

"""
Return the region tag for the region identified by `ridx`.
"""
regiontag(ridx) = unsafe_string(LibLikwid.perfmon_getTagOfRegion(ridx))

"""
Return the group id for the region identified by `ridx`.
"""
regiongroup(ridx) = Int(LibLikwid.perfmon_getGroupOfRegion(ridx))

"""
Return the number of events of the region identified by `ridx`.
"""
regionevents(ridx) = Int(LibLikwid.perfmon_getEventsOfRegion(ridx))

"""
Return the number of metrics of the region identified by `ridx`.
"""
regionmetrics(ridx) = Int(LibLikwid.perfmon_getMetricsOfRegion(ridx))

"""
Return the number of threads of the region identified by `ridx`.
"""
regionthreads(ridx) = Int(LibLikwid.perfmon_getThreadsOfRegion(ridx))

"""
Return the accumulated measurement time for the region identified by `rid` and the thread index `tidx`.
"""
regiontime(ridx, tidx) = LibLikwid.perfmon_getTimeOfRegion(ridx, tidx)

"""
Return the call count for the region identified by `rid` and the thread index `tidx`.
"""
regioncount(ridx, tidx) = LibLikwid.perfmon_getCountOfRegion(ridx, tidx)

"""
Return the call count for the region identified by `ridx`, the event index `eidx` and the thread index `tidx`.
"""
regionresult(ridx, eidx, tidx) = LibLikwid.perfmon_getResultOfRegionThread(ridx, eidx, tidx)

"""
Return the call count for the region identified by `ridx`, the metric index `midx` and the thread index `tidx`.
"""
regionmetric(ridx, midx, tidx) = LibLikwid.perfmon_getMetricOfRegionThread(ridx, midx, tidx)

end # module


text:
## Try to reproduce lexless and lexcmp from version v0.6.3.
## They were removed in v0.7.
## Note that these are in the Symata name space, not Base
## In principle we can name these lexcmp and lexless.
## But, we get deprecation warnings, as if these were the Base function.
## So, we rename them by prefixing with "sym"

"""
    lexcmp(x, y)

Compare `x` and `y` lexicographically and return -1, 0, or 1 depending on whether `x` is
less than, equal to, or greater than `y`, respectively. This function should be defined for
lexicographically comparable types, and `lexless` will call `lexcmp` by default.

```jldoctest
julia> lexcmp("abc", "abd")
-1

julia> lexcmp("abc", "abc")
0
```
"""
symlexcmp(x, y) = cmp(x, y)

"""
    lexless(x, y)

Determine whether `x` is lexicographically less than `y`.

```jldoctest
julia> lexless("abc", "abd")
true
```
"""
symlexless(x, y) = symlexcmp(x,y) < 0

# This has been removed from v0.7
# https://github.com/JuliaLang/julia/issues/5290#issuecomment-324703996
# The argument for removing lexcmp for Complex given in the comment is:
#   Remove lexless and lexcmp for complex numbers. Sorting is the main
#   issue here, which can be easily expressed with sort!(complexes, by=z->(real(z), real(im))).
# But, we are sorting arrays containing more or less anything, not just complex numbers,
function symlexcmp(a::Complex, b::Complex)
    c = cmp(real(a), real(b))
    c == 0 || return c
    cmp(imag(a), imag(b))
end


text:
mutable struct WasmtimeStore
    wasmtime_store_ptr::Ptr{wasmtime_store_t}
    wasmtime_context_ptr::Ptr{wasmtime_context_t}

    WasmtimeStore(store_ptr::Ptr{wasmtime_store_t},
		  ctx_ptr::Ptr{wasmtime_context_t}) = finalizer(wasmtime_store_delete, new(store_ptr, ctx_ptr))
end

function WasmtimeStore(engine)
    store = wasmtime_store_new(engine, C_NULL, C_NULL)
    context = wasmtime_store_context(store)
    WasmtimeStore(store, context)
end

Base.unsafe_convert(::Type{Ptr{wasmtime_store_t}}, store::WasmtimeStore) = store.wasmtime_store_ptr
Base.unsafe_convert(::Type{Ptr{wasmtime_context_t}}, store::WasmtimeStore) = store.wasmtime_context_ptr
Base.show(io::IO, ::WasmtimeStore) = write(io, "WasmtimeStore()")

add_fuel!(store, ⛽) = @wt_check wasmtime_context_add_fuel(store, ⛽)
function consume_fuel!(store, 🛻)
    remaining = Ref(0)
    @wt_check GC.@preserve remaining wasmtime_context_consume_fuel(store, 🛻, Base.pointer_from_objref(remaining))
    remaining[]
end

function fuel_consumed(store)
    remaining = Ref(0)
    fuel_enabled = GC.@preserve remaining wasmtime_context_fuel_consumed(store, Base.pointer_from_objref(remaining))
    if !fuel_enabled
        error("Fuel consumption is not enabled on the engine for this store")
    end
    remaining[]
end

function set_wasi!(store, wasi_config)
    @assert !wasi_config.owned "This WASI configuration is already used by anotherstore"
    @wt_check wasmtime_context_set_wasi(store, wasi_config)
    wasi_config.owned = true
    nothing
end


text:
#!/usr/bin/env julia

using FromFile

@from "Helper.jl" import Helper

function main()
    # fname = "example.txt"
    fname = "input.txt"

    data = Helper.read_lines_as_ints(fname)
    sonar = zip(data, data[2:end], data[3:end])

    numbers = [sum(t) for t in sonar]
    # println(numbers)

    cnt = 0
    for i in 2:length(numbers)
        prev = numbers[i-1]
        curr = numbers[i]
        if curr > prev
            cnt += 1
        end
    end
    println(cnt)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
for input in (1, 0x1, "1")
    s = D.Snowflake(input)
    @test s isa D.Snowflake
    @test s.n == 1
end

s = D.Snowflake(1)
@test s == 1
@test 1 == s
@test sprint(show, s) == "1"
@test JSON3.read("\"1\"", D.Snowflake) == D.Snowflake(1)
@test JSON3.write(s) == "\"1\""


text:
"""
    blockmodel(inblock, outblock[, backbone = blockbackbone(inblock)])

From a `backbone` model, construct a model suitable for learning
a mapping from `inblock` to `outblock`.
"""
function blockmodel end

blockmodel(inblock, outblock) = blockmodel(inblock, outblock, blockbackbone(inblock))

"""
    blockbackbone(inblock)

Create a default backbone that takes in block `inblock`.
"""
function blockbackbone end


text:
@testset "PadicLift" begin
  G = matrix(ResidueRing(ZZ, 8), 3, 3, [4,0,0, 0,2,1, 0,1,2])
  @test (2, 0, 2) == Hecke._last_block_index(G, 2)

  # We want to test for small and big residue rings
  for RR in [Int, ZZ]
    R = ResidueRing(ZZ, RR(3)^6)
    F1 = matrix(R, 3, 3, [2, 0, 0, 0, 2, 0, 0, 0, 1])
    F2 = matrix(R, 3, 3, [0, 1, 0, 1, 2, 2, 0, 2, 1])
    Flist = [F1,F2]
    G = matrix(R, 3, 3, [0, 1, 0, 1, 0, 0, 0, 0, 2])
    for F in Flist
      Flift = Hecke._hensel_qf_modular_odd(G, G, F, 1, 4)
      error = Flift*G*transpose(Flift) - G
      @test 4<=Hecke._min_val(error, 3)
    end

    k = FiniteField(2)[1]
    Y = matrix(k, 3, 3, [0,0,1, 0,0,1, 1,1,0])
    b = [k(i) for i in [1, 0, 0]]
    g = [k(i) for i in [0, 1, 0]]
    X = Hecke._solve_X(Y, b, g)
    @test Y == X + transpose(X)
    for i in 1:3
      @test b[i] == X[i,i] + sum([X[i,j]*g[j] for j in 1:3])
    end

    R = ResidueRing(ZZ, ZZ(3)^5)
    G = diagonal_matrix([R(i) for i in [3^2,1,1]])
    Z = G + matrix(R, 3, 3, [0,3^2,0, 3^2,0,0, 0,0,3])
    F = matrix(R, 3, 3, [1,0,0, 0,0,1, 0,1,0])
    Fn = @inferred Hecke._hensel_qf(Z, G, F, 1, 4, 3)
    @test(4<=Hecke._min_val(Z - F*G*transpose(F),3))

    p = ZZ(3)
    R = ResidueRing(ZZ,p^7)
    G = matrix(R, 6, 6, [0, 243, 0, 0, 0, 0, 243, 0, 0, 0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0])
    F = matrix(R, 6, 6, [0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 2, 1, 0, 1, 0, 0, 1, 2, 1, 0, 0, 0, 2, -1, 6, 3, 0, 1, 1, 9, 3, 6, 1, 0])
    Flift = Hecke.hensel_qf(G, F, 1, 6, p)
    @test Hecke._min_val(Flift*G*transpose(Flift)-G, p ) >= 6

    p = ZZ(2)
    R = ResidueRing(ZZ,p^10)
    U = matrix(R, 2, 2,[0, 1, 1 ,0])
    V = matrix(R, 2, 2,[2, 1, 1 ,2])
    G = diagonal_matrix([2*U, 2*U, V])
    F = matrix(R, 6, 6,[1, 0, 0, 0, 0, 0,
                        1, 1, 1, 1, 0, 0,
                        1, 0, 1, 0, 0, 0,
                        1, 0, 0, 1, 0, 0,
                        1, 0, 0, 0, 1, 1,
                        0, 0, 0, 1, 0, 1])
    Fl = @inferred Hecke.hensel_qf(G, F, 1, 6, p )
    @test 10<=Hecke._min_val(Fl*G*transpose(Fl)-G, p)

    p = ZZ(2)
    R = ResidueRing(ZZ,p^10)
    G = matrix(R, 4, 4,[0,1,0,0,
                        1,0,0,0,
                        0,0,1,0,
                        0,0,0,5])
    F = matrix(R, 4, 4, [1, 0, 0, 0,
                        1, 1, 1, 1,
                        1, 0, 0, 1,
                        1, 0, 1, 0])
    Fl = @inferred Hecke. _hensel_qf_modular_even(G, G, F, 1, 4)
    Er = Fl*G*transpose(Fl) - G
    @test 4<=Hecke._min_val(Er,2)
    @test 5<=Hecke._min_val(diagonal(Er),p)
    F = matrix(R, 4, 4, [0, 1, 0, 0,
                         1, 1, 1, 1,
                         0, 1, 1, 0,
                         0, 1, 0, 1])
    Fl = @inferred Hecke. _hensel_qf_modular_even(G, G, F, 1, 4)
    Er = Fl*G*transpose(Fl) - G
    @test 4<=Hecke._min_val(Er,2)
    @test 5<=Hecke._min_val(diagonal(Er),p)

    G = matrix(R, 4, 4, [2, 1, 0, 0,
                         1, 2, 0, 0,
                         0, 0, 3, 0,
                         0, 0, 0, 7])
    F = matrix(R, 4, 4, [1, 0, 0, 0,
                         1, 1, 0, 0,
                         0, 0, 0, 1,
                         0, 0, 1, 0])
    Fl = @inferred Hecke._hensel_qf_modular_even(G, G, F, 1, 5)
    Er = Fl*G*transpose(Fl) - G
    @test 5<=Hecke._min_val(Er, 2)
    @test 6<=Hecke._min_val(diagonal(Er), p)

    Z = G + identity_matrix(R, 4)*2^2
    Fl = @inferred Hecke. _hensel_qf_modular_even(Z, G, F, 1, 5)
    Er = Fl*G*transpose(Fl) - Z
    @test 5<=Hecke._min_val(Er,2)
    @test 6<=Hecke._min_val(diagonal(Er),p)

    R = ResidueRing(ZZ, ZZ(2)^5)
    G = matrix(R, 3, 3, [2, 1, 0,
                         1, 2, 0,
                         0, 0, 7])
    F = matrix(R, 3, 3, [0, 1, 0,
                         1, 1, 0,
                         0, 0, 1])
    Fl = @inferred Hecke. _hensel_qf_modular_even(G, G, F, 1, 8)
    Er = Fl*G*transpose(Fl) - G
    @test 8<=Hecke._min_val(Er, p)
    @test 9<=Hecke._min_val(diagonal(Er), p)

    G = matrix(R, 3, 3, [4,0,0,
                         0,2,1,
                         0,1,2])
    @test ([1, 2], [2, 0]) == @inferred Hecke._block_indices_vals(G, 2)
  end
end


text:
using Documenter
using MixedModelsSim

makedocs(
    sitename = "MixedModelsSim.jl",
    doctest = true,
    authors = "Phillip Alday, Douglas Bates, Lisa DeBruine, Reinhold Kliegl",
    pages = [
        "index.md",
        "Rapid Start" => "simulation.md",
        "Beginner Friendly Tutorial" => "simulation_tutorial.md",
    ],
)

deploydocs(; repo = "github.com/RePsychLing/MixedModelsSim.jl", push_preview = true)


text:
using sem, Arrow, ModelingToolkit, LinearAlgebra,
    SparseArrays, DataFrames, Optim, LineSearches,
    Statistics

cd("test")

## Observed Data
dat = DataFrame(Arrow.Table("comparisons/data_dem.arrow"))
par_ml = DataFrame(Arrow.Table("comparisons/par_dem_ml.arrow"))
par_ls = DataFrame(Arrow.Table("comparisons/par_dem_ls.arrow"))

dat =
    select(
        dat,
        [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8])
# observed
semobserved = SemObsCommon(data = Matrix{Float64}(dat))

## Model definition
@ModelingToolkit.variables x[1:31]
@Symbolics.variables x[1:31]

#x = rand(31)

S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
    0  0  0  0  0  0  0  0  0  0  0     0     1     0
    0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
    0  0  0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
    0  0  0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

par_order = [collect(21:34); collect(15:20); 2;3; 5;6;7; collect(9:14)]

start_val_ml = Vector{Float64}(par_ml.start[par_order])
start_val_ls = Vector{Float64}(par_ls.start[par_order])
start_val_snlls = Vector{Float64}(par_ls.start[par_order][21:31])

# loss
loss_ml = Loss([SemML(semobserved, [0.0], similar(start_val_ml))])
loss_ls = Loss([sem.SemWLS(semobserved, [0.0], similar(start_val_ml))])
loss_snlls = Loss([sem.SemSWLS(semobserved, [0.0], similar(start_val_ml))])

# imply
imply_ml = ImplySymbolic(A, S, F, x, start_val_ml)
imply_ls = sem.ImplySymbolicWLS(A, S, F, x, start_val_ls)
imply_snlls = sem.ImplySymbolicSWLS(A, S, F, x[21:31], start_val_snlls)

grad_ls = sem.∇SemWLS(loss_ls.functions[1], size(F, 1))

diff_ana =
    SemAnalyticDiff(
        BFGS(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()),# m = 100),
        #P = 0.5*inv(H0),
        #precondprep = (P, x) -> 0.5*inv(FiniteDiff.finite_difference_hessian(model_ls_ana, x))),
        Optim.Options(
            ;f_tol = 1e-10,
            x_tol = 1.5e-8),
            (grad_ls,))

model_ls_ana = Sem(semobserved, imply_ls, loss_ls, diff_ana)

grad_ml = sem.∇SemML_2()

diff_ana_ml =
    SemAnalyticDiff(
        LBFGS(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()),# P = one(rand(31, 31))),
        Optim.Options(
            ;f_tol = 1e-10,
            x_tol = 1.5e-8),
            (grad_ml,))

model_ml_ana = Sem(semobserved, imply_ml, loss_ml, diff_ana_ml)

# fit
solution_ana = sem_fit(model_ls_ana)
solution_ana_ml = sem_fit(model_ml_ana)

all(#
            abs.(solution_ana.minimizer .- par_ls.est[par_order]
                ) .< 0.05*abs.(par_ls.est[par_order]))

all(#
            abs.(solution_ana_ml.minimizer .- par_ml.est[par_order]
                ) .< 0.05*abs.(par_ml.est[par_order]))

imply_ls = sem.ImplySymbolicWLS(A, S, F, x, start_val_ls; hessian = true)

# J = rand(66)
# @benchmark imply_ls.hessian_fun(imply_ls.∇²Σ, J, start_val_ls)

hessian_ls = sem.∇²SemWLS(loss_ls.functions[1])

diff_ana_hes =
    SemAnalyticDiff(
        Newton(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()),# m = 100),
        #P = 0.5*inv(H0),
        #precondprep = (P, x) -> 0.5*inv(FiniteDiff.finite_difference_hessian(model_ls_ana, x))),
        Optim.Options(
            ;f_tol = 1e-10,
            x_tol = 1.5e-8),
            (grad_ls,),
            (hessian_ls,))

model_ls_hes = Sem(semobserved, imply_ls, loss_ls, diff_ana_hes)

#= randpar = rand(31)
grad = similar(randpar)
H = zeros(31, 31)

model_ls_hes(randpar, grad)
model_ls_hes(randpar, H)

maximum(abs.(grad' - FiniteDiff.finite_difference_jacobian(model_ls_ana, randpar)))
maximum(abs.(H - FiniteDiff.finite_difference_hessian(model_ls_ana, randpar))) =#

solution_hes = sem_fit(model_ls_hes)

all(#
            abs.(solution_hes.minimizer .- par_ls.est[par_order]
                ) .< 0.05*abs.(par_ls.est[par_order]))
################################# hessian #####################################
using BenchmarkTools, FiniteDiff

randpar = rand(31)
H = zeros(31,31)
J = rand(31)
grad = rand(31)

@benchmark model_ls_hes.imply.hessian_fun($H, $J, $randpar)
@benchmark model_ls_hes(randpar, grad)

grad - FiniteDiff.finite_difference_gradient(model_ls_ana, randpar)

@benchmark sem_fit(model_ls_hes)
@benchmark sem_fit(model_ls_ana)

@benchmark model_ls_hes(randpar)
@benchmark model_ls_ana(randpar)

grad = similar(randpar)

@benchmark model_ls_hes(randpar, grad)
@benchmark model_ls_ana(randpar, grad)

H = rand(31,31)

@benchmark model_ls_hes(randpar, H)
@benchmark FiniteDiff.finite_difference_hessian(model_ls_ana, randpar)

H0 = FiniteDiff.finite_difference_hessian(model_ls_ana, randpar)
model_ls_hes(randpar, H)
H - H0
#H0 = convert(Matrix{Float64}, H0)
#H0 = inv(H0)
@benchmark model_ls_ana.imply.∇Σ'*diff_ana.functions[1].V*model_ls_ana.imply.∇Σ
@benchmark model_ls_ana.imply.∇Σ'*diff_ana.functions[1].V*model_ls_ana.imply.∇Σ

#function mypre_H0(x)
#    return H0
#end

model_ls_ana(randpar)
gr = similar(randpar)
model_ls_ana(randpar, gr)
initial_P = 2*model_ls_ana.imply.∇Σ'*model_ls_ana.diff.functions[1].V*model_ls_ana.imply.∇Σ

sparse(abs.(initial_P-H0) .> 0.0001)

par = [x[i] for i = 1:31]
invia = sem.neumann_series(A)
imp_cov_sym = F*invia*S*permutedims(invia)*permutedims(F)
imp_cov_sym = Array(imp_cov_sym)
imp_cov_sym = ModelingToolkit.simplify.(imp_cov_sym)
#imp_cov_sym = Symbolics.simplify.(imp_cov_sym)
imp_cov_sym = imp_cov_sym[tril(trues(size(F, 1), size(F, 1)))]
∇Σ_sym = ModelingToolkit.sparsejacobian(imp_cov_sym, par)
∇²Σ_sym = ModelingToolkit.jacobian(vec(permutedims(∇Σ_sym)), x)

@time ∇Σ_sym = Symbolics.sparsejacobian(imp_cov_sym, par)
@time ∇²Σ_sym = Symbolics.jacobian(vec(permutedims(∇Σ_sym)), x)
∇Σ_sym = Array(∇Σ_sym)
H_array = [ModelingToolkit.sparsejacobian(∇Σ_sym[i, :], par) for i = 1:66]
@time H_array = [Symbolics.sparsejacobian(∇Σ_sym[i, :], par) for i = 1:66]

H_array = [ModelingToolkit.sparsehessian(imp_cov_sym[i], par) for i = 1:66]

H_array[2].nzval

function similar_sparse_float(S)
    S_tuple = findnz(S)
    S_tuple_new = (S_tuple[1], S_tuple[2], ones(size(S_tuple[3], 1)))
    new = sparse(S_tuple_new..., size(S)...)
end

H_pre = similar_sparse_float.(H_array)

using StatsBase
nobs = 80
npar = 100
nnd = Int64(nobs*(nobs+1)/2)

testsym = StatsBase.sample(imp_cov_sym, nnd)
testpar = StatsBase.sample(x, npar)

∇testsym = Symbolics.sparsejacobian(testsym, testpar)
∇testsym = Array(∇testsym)
#∇²Σ_sym = ModelingToolkit.jacobian(vec(permutedims(∇Σ_sym)), x)
@time H_array = [Symbolics.sparsejacobian(∇testsym[i, :], testpar) for i = 1:nnd]
# H_array = [sparse(∇²Σ_sym[((i-1)*31+1):i*31, :]) for i = 1:66]

#= jacobian_fun =
    eval(ModelingToolkit.build_function(
        ∇Σ_sym,
        x
    )[2]) =#

hessian_fun =
    eval(ModelingToolkit.build_function(
        ∇²Σ_sym,
        x
    )[2])

H = zeros(2046, 31)
@benchmark hessian_fun(H, randpar)
#H_array = [sparse(H[((i-1)*31+1):i*31, :]) for i = 1:66]

#model_ls_ana(randpar)

J = (-2*(grad_ls.s-model_ls_ana.imply.imp_cov)'*grad_ls.V)'

@variables J[1:nnd]

HT = zeros(Num, npar, npar)
for i in 1:nnd
    HT += J[i]*H_array[i]
end

@variables Jsym[1:66]

HT = zeros(Num, 31,31)
for i in 1:66
    HT += Jsym[i]*H_array[i]
end

HT = simplify.(HT)

hessian_fun =
    eval(ModelingToolkit.build_function(
        HT,
        Jsym,
        par
    )[2])

T = zeros(31,31)
randpar = rand(31)
J = rand(66)

@benchmark hessian_fun(T, J, randpar)

T

maximum(abs.(H0 - (initial_P + T)))

f(Hsum, x, J)

maximum(abs.(H0 - (initial_P+HT)))

vectorized_hessian_fun =
    eval(ModelingToolkit.build_function(
        H_array,
        x
    )[2])

vectorized_hessian_fun_alloc =
    eval(ModelingToolkit.build_function(
        H_array,
        x
    )[1])

function mypre(x)
    return initial_P
end

function precondprep_WLS(P, x, model)
    model.imply.gradient_fun(model.imply.∇Σ, x)
    H = 2*model.imply.∇Σ'*model.diff.functions[1].V*model.imply.∇Σ
    return H
end

A_ldiv_B!(pgr, P, gr) = copyto!(pgr, P \ gr)
dot(x, P, y) = dot(x, P*y)

D = duplication_matrix(observed.n_man)
S = inv(observed.obs_cov)
S = kron(S,S)
V = 0.5*(D'*S*D)

LinearAlgebra.ldiv!(P::Matrix{Float64}, b::Vector{Float64}) = P \ b
LinearAlgebra.ldiv!(x::Vector{Float64}, P::Matrix{Float64}, b::Vector{Float64}) = copyto!(x, P \ b)


text:
using Dates, Git

include("github_tools.jl")

@info "Uploading results to crstnbr/dqmc-benchmarks."
if isfile("result.json")
    dt = string(now())
    dt = replace(dt, ":"=>"--")
    branch = Git.branch()
    push("result.json"; filename="result_$(branch)_$(dt).json")
    if isfile("result.md")
        push("result.md"; filename="result_$(branch)_$(dt).md")
    else
        @warn "uploadresult.jl: No result.md. Pushed only result.json."
    end
else
    @warn "uploadresult.jl: Couldn't find result.json."
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
"""
streamplot(f::function, xinterval, yinterval; kwargs...)

f must either accept `f(::Point)` or `f(x::Number, y::Number)`.
f must return a Point2.

Example:
```julia
v(x::Point2{T}) where T = Point2f(x[2], 4*x[1])
streamplot(v, -2..2, -2..2)
```
## Attributes
$(ATTRIBUTES)

## Implementation
See the function `Makie.streamplot_impl` for implementation details.
"""
@recipe(StreamPlot, f, limits) do scene
    merge(
        Attributes(
            stepsize = 0.01,
            gridsize = (32, 32, 32),
            maxsteps = 500,
            colormap = theme(scene, :colormap),
            colorrange = Makie.automatic,
            arrow_size = 10,
            arrow_head = automatic,
            density = 1.0,
            quality = 16
        ),
        default_theme(scene, Lines) # so that we can theme the lines as needed.
    )
end

function convert_arguments(::Type{<: StreamPlot}, f::Function, xrange, yrange)
    xmin, xmax = extrema(xrange)
    ymin, ymax = extrema(yrange)
    return (f, Rect(xmin, ymin, xmax - xmin, ymax - ymin))
end
function convert_arguments(::Type{<: StreamPlot}, f::Function, xrange, yrange, zrange)
    xmin, xmax = extrema(xrange)
    ymin, ymax = extrema(yrange)
    zmin, zmax = extrema(zrange)
    mini = Vec3f(xmin, ymin, zmin)
    maxi = Vec3f(xmax, ymax, zmax)
    return (f, Rect(mini, maxi .- mini))
end

function convert_arguments(::Type{<: StreamPlot}, f::Function, limits::Rect)
    return (f, limits)
end

scatterfun(N) = N == 2 ? scatter! : meshscatter!

"""
streamplot_impl(CallType, f, limits::Rect{N, T}, resolutionND, stepsize)

Code adapted from an example implementation by Moritz Schauer (@mschauer)
from https://github.com/JuliaPlots/Makie.jl/issues/355#issuecomment-504449775

Background: The algorithm puts an arrow somewhere and extends the
streamline in both directions from there. Then, it chooses a new
position (from the remaining ones), repeating the the exercise until the
streamline gets blocked, from which on a new starting point, the process
repeats.

So, ideally, the new starting points for streamlines are not too close to
current streamlines.

Links:

[Quasirandom sequences](http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/)
"""
function streamplot_impl(CallType, f, limits::Rect{N, T}, resolutionND, stepsize, maxsteps=500, dens=1.0) where {N, T}
    resolution = to_ndim(Vec{N, Int}, resolutionND, last(resolutionND))
    mask = trues(resolution...) # unvisited squares
    arrow_pos = Point{N, Float32}[]
    arrow_dir = Vec{N, Float32}[]
    line_points = Point{N, Float32}[]
    colors = Float64[]
    line_colors = Float64[]
    dt = Point{N, Float32}(stepsize)
    mini, maxi = minimum(limits), maximum(limits)
    r = ntuple(N) do i
        LinRange(mini[i], maxi[i], resolution[i] + 1)
    end
    apply_f(x0, P) = if P <: Point
        f(x0)
    else
        f(x0...)
    end
    # see http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/
    ϕ = (MathConstants.φ, 1.324717957244746, 1.2207440846057596)[N]
    acoeff = ϕ.^(-(1:N))
    n_points = 0 # count visited squares
    ind = 0 # index of low discrepancy sequence
    while n_points < prod(resolution)*min(one(dens), dens) # fill up to 100*dens% of mask
        # next index from low discrepancy sequence
        c = CartesianIndex(ntuple(N) do i
            j = ceil(Int, ((0.5 + acoeff[i]*ind) % 1)*resolution[i])
            clamp(j, 1, size(mask, i))
        end)
        ind += 1
        if mask[c]
            x0 = Point(ntuple(N) do i
                first(r[i]) + (c[i] - 0.5) * step(r[i])
            end)
            point = apply_f(x0, CallType)
            if !(point isa Point2 || point isa Point3)
                error("Function passed to streamplot must return Point2 or Point3")
            end
            pnorm = norm(point)
            push!(arrow_pos, x0)
            push!(arrow_dir, point ./ pnorm)
            push!(colors, pnorm)
            mask[c] = false
            n_points += 1
            for d in (-1, 1)
                n_linepoints = 1
                x = x0
                ccur = c
                push!(line_points, Point{N, Float32}(NaN), x)
                push!(line_colors, 0.0, pnorm)
                while x in limits && n_linepoints < maxsteps
                    point = apply_f(x, CallType)
                    pnorm = norm(point)
                    x = x .+ d .* dt .* point ./ pnorm
                    if !(x in limits)
                        break
                    end
                    # WHAT? Why does point behave different from tuple in this
                    # broadcast
                    idx = CartesianIndex(searchsortedlast.(r, Tuple(x)))
                    if idx != ccur
                        if !mask[idx]
                            break
                        end
                        mask[idx] = false
                        n_points += 1
                        ccur = idx
                    end
                    push!(line_points, x)
                    push!(line_colors, pnorm)
                    n_linepoints += 1
                end
            end
        end
    end

    return (
        arrow_pos,
        arrow_dir,
        line_points,
        colors,
        line_colors,
    )
end

function plot!(p::StreamPlot)
    data = lift(p.f, p.limits, p.gridsize, p.stepsize, p.maxsteps, p.density) do f, limits, resolution, stepsize, maxsteps, density
        P = if applicable(f, Point2f(0)) || applicable(f, Point3f(0))
            Point
        else
            Number
        end
        streamplot_impl(P, f, limits, resolution, stepsize, maxsteps, density)
    end
    lines!(
        p,
        lift(x->x[3], data), color = lift(last, data), colormap = p.colormap, colorrange = p.colorrange,
        linestyle = p.linestyle,
        linewidth = p.linewidth,
        inspectable = p.inspectable,
        transparency = p.transparency
    )
    N = ndims(p.limits[])

    if N == 2 # && scatterplot.markerspace[] == Pixel (default)
        # Calculate arrow head rotations as angles. To avoid distortions from
        # (extreme) aspect ratios we need to project to pixel space and renormalize.
        scene = parent_scene(p)
        rotations = lift(scene.camera.projectionview, scene.px_area, data) do pv, pxa, data
            angles = map(data[1], data[2]) do pos, dir
                pstart = project(scene, pos)
                pstop = project(scene, pos .+ dir)
                pdir = pstop - pstart
                n = norm(pdir)
                if n == 0
                    zero(n)
                else
                    angle = acos(pdir[2] / n)
                    angle = ifelse(pdir[1] > 0, 2pi - angle, angle)
                end
            end
            Billboard(angles)
        end
    else
        rotations = map(x -> x[2], data)
    end

    scatterfun(N)(
        p,
        lift(first, data), markersize = p.arrow_size,
        marker = @lift(arrow_head(N, $(p.arrow_head), $(p.quality))),
        color = lift(x-> x[4], data), rotations = rotations,
        colormap = p.colormap, colorrange = p.colorrange,
        inspectable = p.inspectable, transparency = p.transparency
    )
end


text:
module ComoniconTestUtils

export @test_args,
    @test_kwargs,
    rand_argument,
    rand_option,
    rand_flag,
    rand_node_command,
    rand_leaf_command,
    rand_command,
    rand_input

using Faker
using Random
using Configurations
using Comonicon.AST
using Configurations: Maybe

const args = Ref{Vector{Any}}([])
const kwargs = Ref{Vector{Any}}([])

function test_function(xs...; kw...)
    args[] = [xs...]
    kwargs[] = collect(Any, kw)
    return 1
end

macro test_args(ex)
    quote
        @test $ComoniconTestUtils.args[] == $(ex)
    end |> esc
end

macro test_kwargs(ex)
    quote
        @test $ComoniconTestUtils.kwargs[] == $ex
    end |> esc
end

function rand_name(len = 6)
    a = rand('a':'z') # make sure first 1 are not -
    b = randstring(['a':'z'..., '-'], len - 2)
    c = rand('a':'z') # make sure last 3 are not -
    return a * b * c
end

function rand_description()
    Faker.paragraph(number = 4)
end

@option struct OptionParams
    name::Maybe{String} = nothing
    type = nothing
    short::Maybe{Bool} = nothing
    hint::Maybe{String} = nothing # if the value of hint is nothing, use Some(nothing)
    description::Maybe{String} = nothing
end

maybe(x, or) = isnothing(x) ? or : x

rand_option(; kw...) = rand_option(from_kwargs(OptionParams; kw...))

function rand_option(params::OptionParams)
    name = maybe(params.name, rand_name())
    type = maybe(params.type, rand([Int, Float32, Float64, String]))
    short = maybe(params.short, rand(Bool))
    hint = maybe(params.hint, rand([Faker.word(), nothing]))
    description = maybe(params.description, rand_description())

    # unpack Some(nothing)
    hint = hint isa Some ? hint.value : hint
    Option(; sym = Symbol(name), type, short, hint, description)
end

@option struct FlagParams
    name::Maybe{String} = nothing
    short::Maybe{Bool} = nothing
    description::Maybe{String} = nothing
end

rand_flag(; kw...) = rand_option(from_kwargs(FlagParams; kw...))

function rand_flag(params::FlagParams)
    name = maybe(params.name, rand_name())
    short = maybe(params.short, rand(Bool))
    description = maybe(params.description, rand_description())

    Flag(; sym = Symbol(name), short, description)
end

@option struct ArgParams
    name::Maybe{String} = nothing
    type = nothing
    require::Maybe{Bool} = nothing
    vararg::Maybe{Bool} = nothing
    default::Maybe{String} = nothing
    description::Maybe{String} = nothing
end

rand_argument(; kw...) = rand_argument(from_kwargs(ArgParams; kw...))

function rand_argument(params::ArgParams)
    name = maybe(params.name, rand_name())
    type = maybe(params.type, rand([Int, Float32, Float64, String]))
    require = maybe(params.require, rand(Bool))
    vararg = maybe(params.vararg, rand(Bool))
    default = params.require ? nothing : isnothing(params.default) ? Faker.word() : params.default
    description = maybe(params.description, rand_description())

    Argument(; name, type, vararg, require, default, description)
end

@option struct LeafParams
    name::Maybe{String} = nothing
    nrequire::Int = 2
    noptional::Int = 1
    vararg::Maybe{Bool} = nothing
    noptions::Int = 3
    nflags::Int = 3
    args::ArgParams = ArgParams() # don't generate vararg by default
    options::OptionParams = OptionParams()
    flags::FlagParams = FlagParams()
    description::Maybe{String} = nothing
end

rand_leaf_command(; kw...) = rand_leaf_command(from_kwargs(LeafParams; kw...))

function rand_leaf_command(params::LeafParams)
    name = maybe(params.name, rand_name())
    description = maybe(params.description, rand_description())
    options = Dict{String,Option}()
    flags = Dict{String,Flag}()
    for _ in 1:params.noptions
        opt = rand_option(params.options)
        options[opt.name] = opt
    end

    for _ in 1:params.nflags
        flg = rand_flag(params.flags)
        flags[flg.name] = flg
    end

    args = Argument[]
    for _ in 1:params.nrequire
        push!(
            args,
            rand_argument(;
                name = params.args.name,
                type = params.args.type,
                vararg = false,
                require = true,
                default = params.args.default,
                description = params.args.description,
            ),
        )
    end

    for _ in 1:params.noptional
        push!(
            args,
            rand_argument(;
                name = params.args.name,
                type = params.args.type,
                vararg = false,
                require = false,
                default = nothing,
                description = params.args.description,
            ),
        )
    end

    has_vararg = maybe(params.vararg, rand(Bool))

    if has_vararg
        vararg = rand_argument(;
            name = params.args.name,
            type = params.args.type,
            vararg = true,
            require = false,
            default = nothing,
            description = params.args.description,
        )
    else
        vararg = nothing
    end

    LeafCommand(; fn = test_function, name, args, vararg, options, flags, description)
end

@option struct NodeParams
    name::Maybe{String} = nothing
    nsubcmd::Int = 3
    max_depth::Int = 4 # maximum hierachy of node commands
    subcmd_contains_node::Bool = rand(Bool)
    description::Maybe{String} = nothing
    leaf::LeafParams = LeafParams()
end

rand_node_command(; kw...) = rand_node_command(from_kwargs(NodeParams; kw...))

function rand_node_command(params::NodeParams, depth::Int = 0)
    name = maybe(params.name, rand_name())
    description = maybe(params.description, rand_description())
    subcmds = Dict{String,Any}()
    for _ in 1:params.nsubcmd
        # limit max depth
        if params.subcmd_contains_node && rand() < 0.5 && depth < params.max_depth
            cmd = rand_node_command(params, depth + 1)
            subcmds[cmd.name] = cmd
        else
            cmd = rand_leaf_command(params.leaf)
            subcmds[cmd.name] = cmd
        end
    end

    NodeCommand(; name, subcmds, description)
end

@option struct CommandParams
    version::VersionNumber = v"1.2.3"
    has_node_command::Maybe{Bool} = nothing
    leaf::LeafParams = LeafParams()
    node::NodeParams = NodeParams()
end

rand_command(; kw...) = rand_command(from_kwargs(CommandParams; kw...))

function rand_command(params::CommandParams)
    has_node_command = maybe(params.has_node_command, rand(Bool))
    if has_node_command
        root = rand_node_command(params.node)
    else
        root = rand_leaf_command(params.leaf)
    end

    return Entry(; root, version = params.version)
end

@option struct InputParams
    max_nvararg::Int = 5
    dash::Bool = false
    shuffle::Bool = true # shuffle inputs
    variable_nflags::Bool = true
    variable_noptions::Bool = true
end

rand_input(x; kw...) = rand_input(x, from_kwargs(InputParams; kw...))

rand_input(::Type{Any}, ::InputParams) = Faker.word()
rand_input(::Type{String}, ::InputParams) = randstring('a':'z', 5)
rand_input(::Type{Int}, ::InputParams) = Faker.random_int()
rand_input(::Type{T}, ::InputParams) where {T} = string(rand(T))
rand_input(entry::Entry, params::InputParams) = rand_input(entry.root, params)

function rand_input(entry::LeafCommand, params::InputParams)
    args = String[]
    for each in entry.args
        if each.require
            push!(args, rand_input(each.type, params))
        elseif rand() < 0.5
            push!(args, rand_input(each.type, params))
        end
    end

    if !isnothing(entry.vararg)
        for _ in 1:rand(1:params.max_nvararg)
            push!(args, rand_input(entry.vararg.type, params))
        end
    end

    nflags = params.variable_nflags ? rand(1:length(entry.flags)) : length(entry.flags)
    noptions = params.variable_noptions ? rand(1:length(entry.options)) : length(entry.options)

    flags = String[]
    for (_, flag) in shuffle(collect(entry.flags))[1:nflags]
        if flag.short && rand(Bool)
            push!(flags, "-" * first(flag.name))
        else
            push!(flags, "--" * flag.name)
        end
    end

    options = []
    for (_, option) in shuffle(collect(entry.options))[1:noptions]
        short = option.short && rand(Bool)
        name = short ? "-" * first(option.name) : "--" * option.name
        value = rand_input(option.type, params)

        if rand(Bool) # use_assign
            push!(options, (name * "=" * value,))
        elseif short && rand(Bool) # use -o<value>
            push!(options, (name * value,))
        else # use arg
            push!(options, (name, value))
        end
    end

    inputs = String[]

    while !isempty(args) || !isempty(flags) || !isempty(options)
        if rand(Bool) && !isempty(args)
            push!(inputs, popfirst!(args))
        elseif rand(Bool) && !isempty(flags)
            push!(inputs, pop!(flags))

            # at least push one of them
        elseif !isempty(options)
            push!(inputs, pop!(options)...)
        elseif !isempty(args)
            push!(inputs, popfirst!(args))
        else
            push!(inputs, pop!(flags))
        end
    end
    return inputs
end

function rand_input(cmd::NodeCommand, params::InputParams)
    _, next = rand(cmd.subcmds)
    inputs = rand_input(next, params)
    pushfirst!(inputs, next.name)
    return inputs
end

end


text:
module AbstractControlProcesses

export AbstractProcess, PhysicalProcess, SimulatedProcess
export  num_outputs,
        num_inputs,
        outputrange,
        inputrange,
        isstable,
        isasstable,
        sampletime,
        periodic_wait,
        bias,
        control,
        measure,
        initialize,
        finalize

# Interface specification ===================================================================
"""
    AbstractProcess

Base abstract type for all lab processes. This should not be inherited from directly, see [`PhysicalProcess`](@ref), [`SimulatedProcess`](@ref)
"""
abstract type AbstractProcess end

"""
    PhysicalProcess <: AbstractProcess

Pysical processes should inherit from this abstract type.
"""
abstract type PhysicalProcess  <: AbstractProcess end

"""
    SimulatedProcess <: AbstractProcess

Simulated processes should inherit from this abstract type.
"""
abstract type SimulatedProcess <: AbstractProcess end

## Function definitions =====================================================================
"""
    ny = num_outputs(P::AbstractProcess)

Return the number of outputs (measurement signals) of the process.
"""
num_outputs(p::AbstractProcess) = error("Function not implemented for $(typeof(p))")

"""
    nu = num_inputs(P::AbstractProcess)

Return the number of inputs (control signals) of the process.
"""
num_inputs(p::AbstractProcess) = error("Function not implemented for $(typeof(p))")

"""
    range = outputrange(P::AbstractProcess)

Return the range of outputs (measurement signals) of the process. `range` is a vector of
tuples,  `length(range) = num_outputs(P), eltype(range) = Tuple(Real, Real)`
"""
outputrange(p::AbstractProcess) = error("Function not implemented for $(typeof(p))")

"""
    range = inputrange(P::AbstractProcess)

Return the range of inputs (control signals) of the process. `range` is a vector of
tuples, `length(range) = num_inputs(P), eltype(range) = Tuple(Real, Real)`
"""
inputrange(p::AbstractProcess) = error("Function not implemented for $(typeof(p))")

"""
    isstable(P::AbstractProcess)

Return true/false indicating whether or not the process is stable
"""
isstable(p::AbstractProcess) = error("Function not implemented for $(typeof(p))")

"""
    isasstable(P::AbstractProcess)

Return true/false indicating whether or not the process is asymptotically stable
"""
isasstable(p::AbstractProcess) = error("Function not implemented for $(typeof(p))")

"""
    h = sampletime(P::AbstractProcess)

Return the sample time of the process in seconds.
"""
sampletime(p::AbstractProcess) = error("Function not implemented for $(typeof(p))")

"""
    last_time = periodic_wait(P::AbstractProcess, last_time, dt)

Should continue execution at start of next period, i.e. when `dt` seconds has passed since `last_time`.

For a `PhysicalProcess` this has a default implementation which sleeps until `last_time + dt` using
`Libc.systemsleep`. Though this is more precise than using the standard `sleep` in Julia, it also
blocks the entire Julia thread it is running on which could affect execution in some cases.

For a `SimulatedProcess` it can either be running in realtime in a background process, in which case
it should be implemented similarly to the physical process, or it runs at full simulation speed, in
which case the environment should simulate a step corresponding to `dt` time when `periodic_wait` is called.

A good way to initialize `last_time` to the correct value is to call `last_time = periodic_wait(P, 0.0, 0.0)`.
"""
periodic_wait(p::AbstractProcess, last_time, dt) = error("Function not implemented for $(typeof(p))")
function periodic_wait(p::PhysicalProcess, last_time, dt)
    target_time = last_time + dt
    Libc.systemsleep(max(0, target_time - time()))
    return time()
end

"""
    b = bias(P::AbstractProcess)

Return an input bias for the process. This could be, i.e., the constant input u₀ around which
a nonlinear system is linearized, or whatever other bias might exist on the input.
`length(b) = num_inputs(P)`
"""
bias(p::AbstractProcess) = error("Function not implemented for $(typeof(p))")

"""
    control(P::AbstractProcess, u)

Send a control signal `u` to the process where `u` should have dimension `num_inputs(P)`.
"""
control(p::AbstractProcess, u::AbstractVector) = error("Function not implemented for $(typeof(p))")
control(p::AbstractProcess, u::Number) = control(p, [u])

"""
    y = measure(P::AbstractProcess)

Return a measurement from the process, `y` has length `num_outputs(P)`.
"""
measure(p::AbstractProcess) = error("Function not implemented for $(typeof(p))")

"""
    initialize(P::AbstractProcess)

This function is called before any control or measurement operations are performed. During a call to `initialize`, one might set up external communications etc. After control is done,
the function [`finalize`](@ref) is called.
"""
initialize(p::AbstractProcess) = nothing

"""
    finalize(P::AbstractProcess)

This function is called after any control or measurement operations are performed. During a call to `finalize`, one might finalize external communications etc. Before control is done,
the function [`initialize`](@ref) is called.
"""
Base.finalize(p::AbstractProcess) = nothing

"""
    show(io::IO, p::AbstractProcess)

To make the types show a bit nicer by default and not include all the
struct information.
"""
Base.show(io::IO, p::AbstractProcess) = print(io, "$(typeof(p))()")

end # module