text:
using ArgParse,CairoMakie

CairoMakie.activate!(type="pdf")

function parse_argument()
end

function normal_dist_prob_desity(μ,σ,x)
    (ℯ^(-(x-μ)^2/(2σ)))/sqrt(2*π*σ)
end

function main()
    graph_name = "normal_cumulative_dist"

    xs::Array{Float64} = -4.0:0.1:4.0

    graph_params = [Dict("params"=>Dict("μ"=>0.0,"σ"=>1.0),   "linestyle"=>nothing, "label"=>"N(0.0,1.0)"),
                    Dict("params"=>Dict("μ"=>0.0,"σ"=>0.5),   "linestyle"=>:dash,   "label"=>"N(0.0,0.5)"),
                    Dict("params"=>Dict("μ"=>-1.0,"σ"=>0.75), "linestyle"=>:dot,    "label"=>"N(-1.0,0.75)")]

    fig = Figure(resolution = (800, 600))

    axis = Axis(fig[1,1], title = "normal cumulative distribution function")

    for params in graph_params
        ys::Array{Float64} = []
        accum::Float64 = 0.0
        μ = params["params"]["μ"]
        σ = params["params"]["σ"]
        for x in xs
            accum += normal_dist_prob_desity(μ,σ,x)
            push!(ys,accum)
        end
        lines!(fig[1,1], xs, ys, linestyle=params["linestyle"], label=params["label"])
    end
    axislegend(axis, position = :rb)

    save("$graph_name.pdf", fig)

end

main()

max_stars_repo_name: seisyuu-hantatsushi/julia

max_issues_repo_name: seisyuu-hantatsushi/julia

max_forks_repo_name: seisyuu-hantatsushi/julia


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

module REPLTests
import ..Pkg # ensure we are using the correct Pkg

using Pkg
using Pkg.Types: manifest_info, EnvCache, Context
import Pkg.Types.PkgError
using UUIDs
using Test
import LibGit2

using ..Utils

@testset "help" begin
    pkg"?"
    pkg"?  "
    pkg"?add"
    pkg"? add"
    pkg"?    add"
    pkg"help add"
    @test_throws PkgError pkg"helpadd"
end

temp_pkg_dir() do project_path
    with_pkg_env(project_path; change_dir=true) do;
        pkg"generate HelloWorld"
        LibGit2.close((LibGit2.init(".")))
        cd("HelloWorld")

        @test_throws PkgError pkg"dev Example#blergh"

        @test_throws PkgError pkg"generate 2019Julia"
        pkg"generate Foo"
        pkg"dev Foo"
        mv(joinpath("Foo", "src", "Foo.jl"), joinpath("Foo", "src", "Foo2.jl"))
        @test_throws PkgError pkg"dev Foo"
        ###
        mv(joinpath("Foo", "src", "Foo2.jl"), joinpath("Foo", "src", "Foo.jl"))
        write(joinpath("Foo", "Project.toml"), """
            name = "Foo"
        """
        )
        @test_throws PkgError pkg"dev Foo"
        write(joinpath("Foo", "Project.toml"), """
            uuid = "b7b78b08-812d-11e8-33cd-11188e330cbe"
        """
        )
        @test_throws PkgError pkg"dev Foo"
    end
end

temp_pkg_dir(;rm=false) do project_path; cd(project_path) do;
    tmp_pkg_path = mktempdir()

    pkg"activate ."
    pkg"add Example@0.5"
    @test isinstalled(TEST_PKG)
    v = Pkg.dependencies()[TEST_PKG.uuid].version
    pkg"rm Example"
    pkg"add Example, Random"
    pkg"rm Example Random"
    pkg"add Example,Random"
    pkg"rm Example,Random"
    pkg"add Example#master"

    # Test upgrade --fixed doesn't change the tracking (https://github.com/JuliaLang/Pkg.jl/issues/434)
    entry = Pkg.Types.manifest_info(Context(), TEST_PKG.uuid)
    @test entry.repo.rev == "master"
    pkg"up --fixed"
    entry = Pkg.Types.manifest_info(Context(), TEST_PKG.uuid)
    @test entry.repo.rev == "master"

    pkg"test Example"
    @test isinstalled(TEST_PKG)
    @test Pkg.dependencies()[TEST_PKG.uuid].version > v

    pkg2 = "UnregisteredWithProject"
    pkg2_uuid = UUID("58262bb0-2073-11e8-3727-4fe182c12249")
    p2 = git_init_package(tmp_pkg_path, joinpath(@__DIR__, "test_packages/$pkg2"))
    Pkg.REPLMode.pkgstr("add $p2")
    Pkg.REPLMode.pkgstr("pin $pkg2")
    # FIXME: this confuses the precompile logic to know what is going on with the user
    # FIXME: why isn't this testing the Pkg after importing, rather than after freeing it
    #@eval import Example
    #@eval import $(Symbol(pkg2))
    @test Pkg.dependencies()[pkg2_uuid].version == v"0.1.0"
    Pkg.REPLMode.pkgstr("free $pkg2")
    @test_throws PkgError Pkg.REPLMode.pkgstr("free $pkg2")
    Pkg.test("UnregisteredWithProject")

    write(joinpath(p2, "Project.toml"), """
        name = "UnregisteredWithProject"
        uuid = "58262bb0-2073-11e8-3727-4fe182c12249"
        version = "0.2.0"
        """
    )
    LibGit2.with(LibGit2.GitRepo, p2) do repo
        LibGit2.add!(repo, "*")
        LibGit2.commit(repo, "bump version"; author = TEST_SIG, committer=TEST_SIG)
        pkg"update"
        @test Pkg.dependencies()[pkg2_uuid].version == v"0.2.0"
        Pkg.REPLMode.pkgstr("rm $pkg2")

        c = LibGit2.commit(repo, "empty commit"; author = TEST_SIG, committer=TEST_SIG)
        c_hash = LibGit2.GitHash(c)
        Pkg.REPLMode.pkgstr("add $p2#$c")
    end

    mktempdir() do tmp_dev_dir
    withenv("JULIA_PKG_DEVDIR" => tmp_dev_dir) do
        pkg"develop Example"

        # Copy the manifest + project and see that we can resolve it in a new environment
        # and get all the packages installed
        proj = read("Project.toml", String)
        manifest = read("Manifest.toml", String)
        cd_tempdir() do tmp
            old_depot = copy(DEPOT_PATH)
            try
                empty!(DEPOT_PATH)
                write("Project.toml", proj)
                write("Manifest.toml", manifest)
                mktempdir() do depot_dir
                    pushfirst!(DEPOT_PATH, depot_dir)
                    pkg"instantiate"
                    @test Pkg.dependencies()[pkg2_uuid].version == v"0.2.0"
                end
            finally
                empty!(DEPOT_PATH)
                append!(DEPOT_PATH, old_depot)
            end
        end # cd_tempdir
    end # withenv
    end # mktempdir
end # cd
end # temp_pkg_dir

# issue #904: Pkg.status within a git repo
temp_pkg_dir() do path
    pkg2 = "UnregisteredWithProject"
    p2 = git_init_package(path, joinpath(@__DIR__, "test_packages/$pkg2"))
    Pkg.activate(p2)
    Pkg.status() # should not throw
    Pkg.REPLMode.pkgstr("status") # should not throw
end

temp_pkg_dir() do project_path; cd(project_path) do
    mktempdir() do tmp
        mktempdir() do depot_dir
            old_depot = copy(DEPOT_PATH)
            try
                empty!(DEPOT_PATH)
                pushfirst!(DEPOT_PATH, depot_dir)
                withenv("JULIA_PKG_DEVDIR" => tmp) do
                    # Test an unregistered package
                    p1_path = joinpath(@__DIR__, "test_packages", "UnregisteredWithProject")
                    p1_new_path = joinpath(tmp, "UnregisteredWithProject")
                    cp(p1_path, p1_new_path)
                    Pkg.REPLMode.pkgstr("develop $(p1_new_path)")
                    Pkg.REPLMode.pkgstr("build; precompile")
                    @test realpath(Base.find_package("UnregisteredWithProject")) == realpath(joinpath(p1_new_path, "src", "UnregisteredWithProject.jl"))
                    @test Pkg.dependencies()[UUID("58262bb0-2073-11e8-3727-4fe182c12249")].version == v"0.1.0"
                    Pkg.test("UnregisteredWithProject")
                end
            finally
                empty!(DEPOT_PATH)
                append!(DEPOT_PATH, old_depot)
            end
        end # withenv
    end # mktempdir
    # nested
    mktempdir() do other_dir
        mktempdir() do tmp;
            cd(tmp)
            pkg"generate HelloWorld"
            cd("HelloWorld") do
                with_current_env() do
                    uuid1 = Pkg.generate("SubModule1")["SubModule1"]
                    uuid2 = Pkg.generate("SubModule2")["SubModule2"]
                    pkg"develop SubModule1"
                    mkdir("tests")
                    cd("tests")
                    pkg"develop ../SubModule2"
                    @test Pkg.dependencies()[uuid1].version == v"0.1.0"
                    @test Pkg.dependencies()[uuid2].version == v"0.1.0"
                    # make sure paths to SubModule1 and SubModule2 are relative
                    manifest = Pkg.Types.Context().env.manifest
                    @test manifest[uuid1].path == "SubModule1"
                    @test manifest[uuid2].path == "SubModule2"
                end
            end
            cp("HelloWorld", joinpath(other_dir, "HelloWorld"))
            cd(joinpath(other_dir, "HelloWorld"))
            with_current_env() do
                # Check that these didn't generate absolute paths in the Manifest by copying
                # to another directory
                @test Base.find_package("SubModule1") == joinpath(pwd(), "SubModule1", "src", "SubModule1.jl")
                @test Base.find_package("SubModule2") == joinpath(pwd(), "SubModule2", "src", "SubModule2.jl")
            end
        end
    end
end # cd
end # temp_pkg_dir

# activate
temp_pkg_dir() do project_path
    cd_tempdir() do tmp
        path = pwd()
        pkg"activate ."
        @test Base.active_project() == joinpath(path, "Project.toml")
        # tests illegal names for shared environments
        @test_throws Pkg.Types.PkgError pkg"activate --shared ."
        @test_throws Pkg.Types.PkgError pkg"activate --shared ./Foo"
        @test_throws Pkg.Types.PkgError pkg"activate --shared Foo/Bar"
        @test_throws Pkg.Types.PkgError pkg"activate --shared ../Bar"
        # check that those didn't change te enviroment
        @test Base.active_project() == joinpath(path, "Project.toml")
        mkdir("Foo")
        cd(mkdir("modules")) do
            pkg"generate Foo"
        end
        pkg"develop modules/Foo"
        pkg"activate Foo" # activate path Foo over deps Foo
        @test Base.active_project() == joinpath(path, "Foo", "Project.toml")
        pkg"activate ."
        #=@test_logs (:info, r"activating new environment at ")))=# pkg"activate --shared Foo" # activate shared Foo
        @test Base.active_project() == joinpath(Pkg.envdir(), "Foo", "Project.toml")
        pkg"activate ."
        rm("Foo"; force=true, recursive=true)
        pkg"activate Foo" # activate path from developed Foo
        @test Base.active_project() == joinpath(path, "modules", "Foo", "Project.toml")
        pkg"activate ."
        #=@test_logs (:info, r"activating new environment at ")=# pkg"activate ./Foo" # activate empty directory Foo (sidestep the developed Foo)
        @test Base.active_project() == joinpath(path, "Foo", "Project.toml")
        pkg"activate ."
        #=@test_logs (:info, r"activating new environment at ")=# pkg"activate Bar" # activate empty directory Bar
        @test Base.active_project() == joinpath(path, "Bar", "Project.toml")
        pkg"activate ."
        pkg"add Example" # non-deved deps should not be activated
        #=@test_logs (:info, r"activating new environment at ")=# pkg"activate Example"
        @test Base.active_project() == joinpath(path, "Example", "Project.toml")
        pkg"activate ."
        cd(mkdir("tests"))
        pkg"activate Foo" # activate developed Foo from another directory
        @test Base.active_project() == joinpath(path, "modules", "Foo", "Project.toml")
        tmpdepot = mktempdir(tmp)
        tmpdir = mkpath(joinpath(tmpdepot, "environments", "Foo"))
        push!(Base.DEPOT_PATH, tmpdepot)
        pkg"activate --shared Foo" # activate existing shared Foo
        @test Base.active_project() == joinpath(tmpdir, "Project.toml")
        pop!(Base.DEPOT_PATH)
        pkg"activate" # activate home project
        @test Base.ACTIVE_PROJECT[] === nothing
        # expansion of ~
        if !Sys.iswindows()
            pkg"activate ~/Foo"
            @test Base.active_project() == joinpath(homedir(), "Foo", "Project.toml")
        end
    end
end

# path should not be relative when devdir() happens to be in project
# unless user used dev --local.
temp_pkg_dir() do depot
    cd_tempdir() do tmp
        uuid = UUID("7876af07-990d-54b4-ab0e-23690620f79a") # Example
        pkg"activate ."
        withenv("JULIA_PKG_DEVDIR" => joinpath(pwd(), "dev")) do
            pkg"dev Example"
            @test manifest_info(Context(), uuid).path == joinpath(pwd(), "dev", "Example")
            pkg"dev --shared Example"
            @test manifest_info(Context(), uuid).path == joinpath(pwd(), "dev", "Example")
            pkg"dev --local Example"
            @test manifest_info(Context(), uuid).path == joinpath("dev", "Example")
        end
    end
end

test_complete(s) = Pkg.REPLMode.completions(s, lastindex(s))
apply_completion(str) = begin
    c, r, s = test_complete(str)
    str[1:prevind(str, first(r))]*first(c)
end

# Autocompletions
temp_pkg_dir() do project_path; cd(project_path) do
    @testset "tab completion" begin
        pkg"registry add General" # instantiate the `General` registry to complete remote package names
        pkg"activate ."
        c, r = test_complete("add Exam")
        @test "Example" in c
        c, r = test_complete("rm Exam")
        @test isempty(c)

        Pkg.REPLMode.pkgstr("develop $(joinpath(@__DIR__, "test_packages", "PackageWithDependency"))")

        c, r = test_complete("rm PackageWithDep")
        @test "PackageWithDependency" in c
        c, r = test_complete("rm -p PackageWithDep")
        @test "PackageWithDependency" in c
        c, r = test_complete("rm --project PackageWithDep")
        @test "PackageWithDependency" in c
        c, r = test_complete("rm Exam")
        @test isempty(c)
        c, r = test_complete("rm -p Exam")
        @test isempty(c)
        c, r = test_complete("rm --project Exam")
        @test isempty(c)

        c, r = test_complete("rm -m PackageWithDep")
        @test "PackageWithDependency" in c
        c, r = test_complete("rm --manifest PackageWithDep")
        @test "PackageWithDependency" in c
        c, r = test_complete("rm -m Exam")
        @test "Example" in c
        c, r = test_complete("rm --manifest Exam")
        @test "Example" in c

        c, r = test_complete("rm PackageWithDep")
        @test "PackageWithDependency" in c
        c, r = test_complete("rm Exam")
        @test isempty(c)
        c, r = test_complete("rm -m Exam")
        c, r = test_complete("rm -m Exam")
        @test "Example" in c

        pkg"add Example"
        c, r = test_complete("rm Exam")
        @test "Example" in c
        c, r = test_complete("up --man")
        @test "--manifest" in c
        c, r = test_complete("rem")
        @test "remove" in c
        @test apply_completion("rm E") == "rm Example"
        @test apply_completion("add Exampl") == "add Example"

        # stdlibs
        c, r = test_complete("add Stat")
        @test "Statistics" in c
        c, r = test_complete("add Lib")
        @test "LibGit2" in c
        c, r = test_complete("add REPL")
        @test "REPL" in c

        # upper bounded
        c, r = test_complete("add Chu")
        @test !("Chunks" in c)

        # local paths
        mkpath("testdir/foo/bar")
        c, r = test_complete("add ")
        @test Sys.iswindows() ? ("testdir\\\\" in c) : ("testdir/" in c)
        @test "Example" in c
        @test apply_completion("add tes") == (Sys.iswindows() ? "add testdir\\\\" : "add testdir/")
        @test apply_completion("add ./tes") == (Sys.iswindows() ? "add ./testdir\\\\" : "add ./testdir/")
        c, r = test_complete("dev ./")
        @test (Sys.iswindows() ? ("testdir\\\\" in c) : ("testdir/" in c))

        # complete subdirs
        c, r = test_complete("add testdir/f")
        @test Sys.iswindows() ? ("foo\\\\" in c) : ("foo/" in c)
        @test apply_completion("add testdir/f") == (Sys.iswindows() ? "add testdir/foo\\\\" : "add testdir/foo/")
        # dont complete files
        touch("README.md")
        c, r = test_complete("add RE")
        @test !("README.md" in c)

        # Expand homedir and
        if !Sys.iswindows()
            dirname = "JuliaPkgTest744a757c-d313-11e9-1cac-118368d5977a"
            tildepath = "~/$dirname"
            try
                mkdir(expanduser(tildepath))
                c, r = test_complete("dev ~/JuliaPkgTest744a75")
                @test joinpath(homedir(), dirname, "") in c
            finally
                rm(expanduser(tildepath); force = true)
            end
            c, r = test_complete("dev ~")
            @test joinpath(homedir(), "") in c
        end

        # activate
        pkg"activate --shared FooBar"
        pkg"add Example"
        pkg"activate ."
        c, r = test_complete("activate --shared ")
        @test "FooBar" in c

        # invalid options
        c, r = test_complete("rm -rf ")
        @test isempty(c)

        # parse errors should not throw
        _ = test_complete("add \"Foo")
        # invalid option should not throw
        _ = test_complete("add -z Foo")
        _ = test_complete("add --dontexist Foo")
    end # testset
end end

temp_pkg_dir() do project_path; cd(project_path) do
    mktempdir() do tmp
        cp(joinpath(@__DIR__, "test_packages", "BigProject"), joinpath(tmp, "BigProject"))
        cd(joinpath(tmp, "BigProject"))
        with_current_env() do
            # the command below also tests multiline input
            pkg"""
                dev RecursiveDep2
                dev RecursiveDep
                dev SubModule
                dev SubModule2
                add Random
                add Example
                add JSON
                build
            """
            @eval using BigProject
            pkg"build BigProject"
            @test_throws PkgError pkg"add BigProject"
            # the command below also tests multiline input
            Pkg.REPLMode.pkgstr("""
                test SubModule
                test SubModule2
                test BigProject
                test
                """)
            json_uuid = Pkg.project().dependencies["JSON"]
            current_json = Pkg.dependencies()[json_uuid].version
            old_project = read("Project.toml", String)
            open("Project.toml"; append=true) do io
                print(io, """

                [compat]
                JSON = "0.18.0"
                """
                )
            end
            pkg"up"
            @test Pkg.dependencies()[json_uuid].version.minor == 18
            write("Project.toml", old_project)
            pkg"up"
            @test Pkg.dependencies()[json_uuid].version == current_json
        end
    end
end; end

temp_pkg_dir() do project_path
    cd(project_path) do
        @testset "add/remove using quoted local path" begin
            # utils
            setup_package(parent_dir, pkg_name) = begin
                mkdir(parent_dir)
                cd(parent_dir) do
                    withenv("USER" => "Test User") do
                        Pkg.generate(pkg_name)
                    end
                    cd(pkg_name) do
                        git_init_and_commit(joinpath(project_path, parent_dir, pkg_name))
                    end #cd pkg_name
                end # cd parent_dir
            end

            # extract uuid from a Project.toml file
            extract_uuid(toml_path) = begin
                uuid = ""
                for line in eachline(toml_path)
                    m = match(r"uuid = \"(.+)\"", line)
                    if m !== nothing
                        uuid = m.captures[1]
                        break
                    end
                end
                return uuid
            end

            # testing local dir with space in name
            dir_name = "space dir"
            pkg_name = "WeirdName77"
            setup_package(dir_name, pkg_name)
            uuid = extract_uuid("$dir_name/$pkg_name/Project.toml")
            Pkg.REPLMode.pkgstr("add \"$dir_name/$pkg_name\"")
            @test isinstalled((name=pkg_name, uuid = UUID(uuid)))
            Pkg.REPLMode.pkgstr("remove \"$pkg_name\"")
            @test !isinstalled((name=pkg_name, uuid = UUID(uuid)))

            # testing dir name with significant characters
            dir_name = "some@d;ir#"
            pkg_name = "WeirdName77"
            setup_package(dir_name, pkg_name)
            uuid = extract_uuid("$dir_name/$pkg_name/Project.toml")
            Pkg.REPLMode.pkgstr("add \"$dir_name/$pkg_name\"")
            @test isinstalled((name=pkg_name, uuid = UUID(uuid)))
            Pkg.REPLMode.pkgstr("remove '$pkg_name'")
            @test !isinstalled((name=pkg_name, uuid = UUID(uuid)))

            # more complicated input
            ## pkg1
            dir1 = "two space dir"
            pkg_name1 = "name1"
            setup_package(dir1, pkg_name1)
            uuid1 = extract_uuid("$dir1/$pkg_name1/Project.toml")

            ## pkg2
            dir2 = "two'quote'dir"
            pkg_name2 = "name2"
            setup_package(dir2, pkg_name2)
            uuid2 = extract_uuid("$dir2/$pkg_name2/Project.toml")

            Pkg.REPLMode.pkgstr("add '$dir1/$pkg_name1' \"$dir2/$pkg_name2\"")
            @test isinstalled((name=pkg_name1, uuid = UUID(uuid1)))
            @test isinstalled((name=pkg_name2, uuid = UUID(uuid2)))
            Pkg.REPLMode.pkgstr("remove '$pkg_name1' $pkg_name2")
            @test !isinstalled((name=pkg_name1, uuid = UUID(uuid1)))
            @test !isinstalled((name=pkg_name2, uuid = UUID(uuid2)))

            Pkg.REPLMode.pkgstr("add '$dir1/$pkg_name1' \"$dir2/$pkg_name2\"")
            @test isinstalled((name=pkg_name1, uuid = UUID(uuid1)))
            @test isinstalled((name=pkg_name2, uuid = UUID(uuid2)))
            Pkg.REPLMode.pkgstr("remove '$pkg_name1' \"$pkg_name2\"")
            @test !isinstalled((name=pkg_name1, uuid = UUID(uuid1)))
            @test !isinstalled((name=pkg_name2, uuid = UUID(uuid2)))
        end
    end
end

@testset "parse package url win" begin
    @test typeof(Pkg.REPLMode.parse_package_identifier("https://github.com/abc/ABC.jl";
                                                       add_or_develop=true)) == Pkg.Types.PackageSpec
end

@testset "unit test for REPLMode.promptf" begin
    function set_name(projfile_path, newname)
        sleep(1.1)
        project = Pkg.TOML.parsefile(projfile_path)
        project["name"] = newname
        open(projfile_path, "w") do io
            Pkg.TOML.print(io, project)
        end
    end

    with_temp_env("SomeEnv") do
        @test Pkg.REPLMode.promptf() == "(SomeEnv) pkg> "
    end

    env_name = "Test2"
    with_temp_env(env_name) do env_path
        projfile_path = joinpath(env_path, "Project.toml")
        @test Pkg.REPLMode.promptf() == "($env_name) pkg> "

        newname = "NewName"
        set_name(projfile_path, newname)
        @test Pkg.REPLMode.promptf() == "($newname) pkg> "
        cd(env_path) do
            @test Pkg.REPLMode.promptf() == "($newname) pkg> "
        end
        @test Pkg.REPLMode.promptf() == "($newname) pkg> "

        newname = "NewNameII"
        set_name(projfile_path, newname)
        cd(env_path) do
            @test Pkg.REPLMode.promptf() == "($newname) pkg> "
        end
        @test Pkg.REPLMode.promptf() == "($newname) pkg> "
    end
end

@testset "test" begin
    temp_pkg_dir() do project_path; cd_tempdir() do tmpdir; with_temp_env() do;
        Pkg.add("Example")
        @test_throws PkgError Pkg.REPLMode.pkgstr("test --project Example")
        Pkg.REPLMode.pkgstr("test --coverage Example")
        Pkg.REPLMode.pkgstr("test Example")
    end
    end
    end
end

@testset "activate" begin
    temp_pkg_dir() do project_path; cd_tempdir() do tmpdir; with_temp_env() do;
        mkdir("Foo")
        pkg"activate"
        default = Base.active_project()
        pkg"activate Foo"
        @test Base.active_project() == joinpath(pwd(), "Foo", "Project.toml")
        pkg"activate"
        @test Base.active_project() == default
    end end end
end

@testset "status" begin
    temp_pkg_dir() do project_path
        pkg"""
        add Example Random
        status
        status -m
        status Example
        status Example=7876af07-990d-54b4-ab0e-23690620f79a
        status 7876af07-990d-54b4-ab0e-23690620f79a
        status Example Random
        status -m Example
        """
        # --diff option
        @test_logs (:warn, r"diff option only available") pkg"status --diff"
        @test_logs (:warn, r"diff option only available") pkg"status -d"
        git_init_and_commit(project_path)
        @test_logs () pkg"status --diff"
        @test_logs () pkg"status -d"
    end
end

@testset "subcommands" begin
    temp_pkg_dir() do project_path; cd_tempdir() do tmpdir; with_temp_env() do
        Pkg.REPLMode.pkg"package add Example"
        @test isinstalled(TEST_PKG)
        Pkg.REPLMode.pkg"package rm Example"
        @test !isinstalled(TEST_PKG)
    end end end
end

@testset "REPL API `up`" begin
    # errors
    temp_pkg_dir() do project_path; with_temp_env() do;
        @test_throws PkgError Pkg.REPLMode.pkgstr("up --major --minor")
    end end
end

end # module

ext: jl


text:
function postprocess!(doc::Org)
    # Radio Links
    radios = filter(c -> c isa RadioTarget, doc.components)
    if length(radios) > 0
        paragraphs = filter(c -> c isa Paragraph, doc.components)
        footnoterefs = filter(c -> c isa FootnoteReference{<:Any, Vector{Object}}, doc.components)
        markups = filter(c -> c isa TextMarkup{Vector{Object}}, doc.components)
        for radio in radios
            for heading in doc.headings
                heading.title = radio_replace(heading.title, radio)
            end
            for paragraph in paragraphs
                paragraph.contents = radio_replace(paragraph.contents, radio)
            end
            for markup in markups
                markup.contents = radio_replace(markup.contents, radio)
            end
            for fnref in footnoterefs
                fnref.definition = radio_replace(fnref.definition, radio)
            end
        end
    end
end

function radio_replace(objs::Vector{Object}, radio::RadioTarget)
    if length(radio.contents) == 1 && radio.contents[1] isa TextPlain
        radio_replace_singular(objs, radio)
    else
        needles = radio.contents
        function radioat(i)
            (objs[i] == needles[1] || # 1 matches, or
             (objs[i] isa TextPlain && needles[1] isa TextPlain &&
              endswith(objs[i].text, needles[1].text))) && # 1 ends with needles[1]
                objs[i+1:i+length(needles)-2] == needles[2:end-1] && # 2-(end-1) matches
                (objs[i+length(needles)-1] == needles[end] || # end matches, or
                 (objs[i+length(needles)-1] isa TextPlain && needles[end] isa TextPlain &&
                  startswith(objs[i+length(needles)-1].text, needles[end].text))) # end start with radio[end]
        end
        newobjs = Object[]
        i = 1
        while i <= length(objs)
            if radioat(i)
                if needles[1] isa TextPlain &&
                    objs[i] != needles[1]
                    # objs[1] must be a TextPlain that ends in needles[1]
                    push!(newobjs, TextPlain(@inbounds @view objs[i].text[1:end-ncodeunits(needles[1].text)]))
                end
                push!(newobjs, RadioLink(radio))
                i += length(needles)-1
                if needles[end] isa TextPlain &&
                    objs[i-1] != needles[end]
                    push!(newobjs, TextPlain(@inbounds @view objs[i-1].text[ncodeunits(needles[end].text):end]))
                end
            else
                push!(newobjs, objs[i])
            end
            i += 1
        end
        newobjs
    end
end

function radio_replace_singular(objs::Vector, radio::RadioTarget)
    @assert length(radio.contents) == 1 && radio.contents[1] isa TextPlain
    needle = radio.contents[1].text
    alphnum(c::Char) = c in 'A':'Z' || c in 'a':'z' || c in '0':'9'
    newobjs = Object[]
    for obj in objs
        if obj isa TextPlain && occursin(needle, obj.text)
            occurances = findall(needle, obj.text)
            point = 1
            for occ in occurances
                if (occ.start == 1 || !alphnum(obj.text[prevind(obj.text, occ.start)])) &&
                   (occ.stop == length(obj.text) || !alphnum(obj.text[nextind(obj.text, occ.stop)]))
                    push!(newobjs, TextPlain(@inbounds @view obj.text[point:occ.start-1]))
                    push!(newobjs, RadioLink(radio))
                    point = occ.stop + 1
                end
            end
            if point <= lastindex(obj.text)
                push!(newobjs, TextPlain(@inbounds @view obj.text[point:end]))
            end
        else
            push!(newobjs, obj)
        end
    end
    newobjs
end


text:
#MAXIMIN or MINIMAX neighbourhood
function analyse_NH(
    NH_start::Neighbourhood,
    ata_model::MinimaxModel;
    fF = true,
    n_fill = 1,
    opt_feas = 0.9,
    max_conv = 1,
    start_temp = 1000.0,
    geom_temp = 0.1,
    n_item_sample = 1,
    n_test_sample = 1,
    verbosity = 1,
    start_time = 0,
    max_time = 1000,
)
    if fF == true
        NH_start.obj = zeros(Float64, ata_model.settings.T)
    end
    NH₁ = Neighbourhood()
    NH₀ = Neighbourhood()
    NH⁺ = Neighbourhood()
    NH₁ = _mycopy(NH_start, NH₁)
    NH₀ = _mycopy(NH_start, NH₀)
    NH⁺ = _mycopy(NH_start, NH⁺)
    f_star = ones(2) .* Inf
    f_evals = 0
    t = copy(start_temp)
    T = ata_model.settings.T
    n_items = ata_model.settings.n_items
    n_fs = ata_model.settings.n_fs
    fs_counts = ata_model.settings.fs.counts * ones(Float64, T)'

    round = 1

    if n_fill > 0
        #Fill up
        println("Fill-up starting...")
        for round = 1:n_fill
            NH₁ = _mycopy(NH_start, NH₁)
            warmup = fill(true, T)
            while any(warmup) #filling forms
                v = findfirst(warmup .== true)
                constraints = ata_model.constraints[v]
                #println("ol = ", NH₁.ol)
                fₜ = (1 - opt_feas) * (NH₁.infeas + NH₁.ol) - (opt_feas * NH₁.obj)
                mm = fₜ[v]
                for v2 in findall(warmup .== true)
                    if fₜ[v2] > mm
                        mm = copy(fₜ[v2])
                        v = copy(v2)
                    end
                end

                #filling
                n_t = LinearAlgebra.dot(NH₁.x[:, v], ata_model.settings.fs.counts)
                #println("test ", v, " chosen, length was: ", n_t)
                #try to add other items, the first time it goes over n_max it stops
                if n_t < constraints.length_max
                    if opt_feas == 0 || fF == true
                        NH_add = find_best_itemᵥ(
                            NH₁,
                            v,
                            ata_model.settings.iu,
                            ata_model.settings.fs,
                            constraints,
                            ata_model.settings.forced0[v],
                            n_items,
                            n_fs,
                            ata_model.settings.ol_max[:, v],
                            ata_model.settings.to_apply,
                        )
                        NH_add.f = opt_feas * NH_add.f
                    else
                        NH_add = find_best_itemᵥ(
                            NH₁,
                            ata_model.obj.cores[v],
                            opt_feas,
                            v,
                            ata_model.settings.iu,
                            ata_model.settings.fs,
                            constraints,
                            ata_model.settings.forced0[v],
                            n_items,
                            n_fs,
                            ata_model.settings.ol_max[:, v],
                            ata_model.settings.to_apply,
                        )
                    end
                    if ata_model.settings.to_apply[3]
                        NH_add.ol = eval_overlap(
                            NH_add.x,
                            fs_counts,
                            ata_model.settings.ol_max,
                            T,
                            NH_add.ol,
                        )
                    end
                    Printf.@printf "."
                    n_t = LinearAlgebra.dot(NH_add.x[:, v], ata_model.settings.fs.counts)
                    #println("length for test ", v, ": ", n_t)
                    if n_t <= constraints.length_max
                        NH₁ = _mycopy(NH_add, NH₁)
                    else
                        warmup[v] = false
                        #NH₁ = _mycopy(NH_add, NH₁)
                        Printf.@printf("\n f₁:  %16.3f", NH₁.f)
                        println("-Test ", v, " filled up with ", n_t, " items, ")
                    end
                else
                    warmup[v] = false
                    println("-Test ", v, " filled up with ", n_t, " items, ")
                end
            end
        end #end of round, round = nRound
        println("End of fill-up")
        NH₀ = _mycopy(NH₁, NH₀)
        NH₀.f = _comp_f(NH₀, opt_feas)
        NH⁺ = _mycopy(NH₀, NH⁺)
        if sum(NH₀.infeas + NH₀.ol) + NH₀.iu == 0
            println("Feasible solution found in fill-up")
            fF = false
        else
            println("Feasible solution not found in fill-up")
        end
    end
    if verbosity > 1
        print_neighbourhood(NH⁺)
    end
    coverage_ok = 0
    if n_test_sample > T
        n_test_sample = T
    end
    convergence = 0

    while coverage_ok == 0
        NH₁ = _mycopy(NH₀, NH₁)
        weights = (1 - opt_feas) .* (NH₀.infeas + NH₀.ol) - opt_feas .* (NH₀.obj)
        weights = (weights .- minimum(weights)) .+ 1
        weights = weights ./ sum(weights)
        weights = StatsBase.ProbabilityWeights(weights)
        #determine test order
        #testOrder = StatsBase.sample(collect(1:T), weights, n_test_sample, replace = false)
        #Random.shuffle!(testOrder)
        testOrder = sortperm(weights, rev = true)
        #println("test order = ", Int.(testOrder))
        v₂ = 0
        exit = 0
        #iteratorTestItem = vec(collect(Iterators.product(collect(1:n_item_sample), testOrder[1:n_test_sample])))
        #println(iteratorTestItem[1:nItemoStatsBase.sample])
        #it = 0
        #xnew = copy(NH₀.x)
        while exit == 0 && v₂ < n_test_sample #it<size(iteratorTestItem, 1) #
            #it+= 1
            v₂ += 1
            exit = 0
            v = testOrder[v₂]
            x_forced0ₜ = ata_model.settings.forced0[v]
            Constraintsₜ = ata_model.constraints[v]
            coreₜ = ata_model.obj.cores[v]
            ol_maxₜ = ata_model.settings.ol_max[:, v]
            #it<size(iteratorTestItem, 1) #
            #v = iteratorTestItem[it][2]
            #NH₀.x = copy(xnew)
            taken_items = findall(isone.(NH₀.x[:, v]))
            if n_item_sample > size(taken_items, 1)
                nI = Int(size(taken_items, 1))
            else
                nI = n_item_sample
            end
            taken_items = Random.shuffle!(taken_items) # !removed
            # exit2 = 0
            # if iteratorTestItem[it][1]>size(taken_items, 1)
            # 	exit2 = 1
            # end
            #if exit2 == 0
            # if v!= v2
            # 	v = copy(v2)
            #
            # end
            #h = taken_items[iteratorTestItem[it][1]]
            #println("test ", v₂, " of ", size(testOrder, 1))
            #fix test features
            h₂ = 0
            while exit == 0 && h₂ < nI
                NH₁ = _mycopy(NH₀, NH₁)
                h₂ += 1
                #println("item ", h₂, " of ", size(taken_items, 1))
                #try to remove h
                h = taken_items[h₂]
                #iu = copy(iu₀)
                #if rand()>0.0
                NH₁.x[h, v] = zero(Float64)
                #iu[h]-= 1
                taken = 0
                #else
                #	taken = 1
                #end
                if sum(NH₁.x[:, v] .* ata_model.settings.fs.counts) >=
                   Constraintsₜ.length_min
                    NH₁.infeas[v], x_Iₜ = check_feas(
                        ata_model.settings.fs,
                        Constraintsₜ,
                        NH₁.x[:, v],
                        n_fs,
                        n_items,
                    )
                    if ata_model.settings.to_apply[1]
                        iu_max = sum(NH₁.x, dims = 2)[:, 1] - ata_model.settings.iu.max
                        iu_max = sum_pos(iu_max)
                    else
                        iu_max = 0
                    end
                    if ata_model.settings.to_apply[2]
                        iu_min = -sum(NH₁.x, dims = 2)[:, 1] + ata_model.settings.iu.min
                        iu_min = sum_pos(iu_min)
                    else
                        iu_min = 0
                    end
                    NH₁.iu = iu_max + iu_min
                    if ata_model.settings.to_apply[3]
                        NH₁.ol = eval_overlap(
                            NH₁.x,
                            fs_counts,
                            ata_model.settings.ol_max,
                            T,
                            NH₁.ol,
                        )
                    end
                    #NH₁.ol[v] = eval_overlapᵥ(NH₁.x[:, v], NH₁.x, ata_model.settings.fs.counts, ol_maxₜ, v)
                    if fF == false
                        NH₁.obj[v] = eval_TIFₜ(x_Iₜ, coreₜ)
                    end
                    NH₁.f = _comp_f(NH₁, opt_feas)
                    f_evals += 1
                    if (NH₁.f <= NH₀.f)
                        #switch item
                        #NH₀.ol = eval_overlap(NH₀.x, fs_counts, ata_model.settings.ol_max, T, NH₀.ol)
                        #NH₀.f = _comp_f(NH₀, opt_feas)
                        NH₀ = _mycopy(NH₁, NH₀)
                        # println("better to remove, new f₀ = ", NH₀.f)
                        if (NH₁.f < NH⁺.f)
                            exit = 1
                            convergence = 0
                            NH⁺ = _mycopy(NH₁, NH⁺)
                            if verbosity == 2
                                Printf.@printf("\n f⁺ : %5.3f", NH⁺.f)
                                Printf.@printf(" in test : %5d", v)
                            else
                                Printf.@printf "+"
                            end
                        end
                    else
                        p = _sig_c(-(NH₁.f - NH₀.f) / t) # -(NH₁.f - NH₀.f) / t always negative
                        if p < 1e-10
                            p = zeros(Float64)[1]
                        end
                        if (rand() < p)
                            #remove
                            #exit = 1
                            NH₀ = _mycopy(NH₁, NH₀)
                            #NH₀.ol = eval_overlap(NH₀.x, fs_counts, ata_model.settings.ol_max, T, NH₀.ol)
                            #NH₀.f = _comp_f(NH₀, opt_feas)
                            if verbosity == 2
                                Printf.@printf("\n f₀ : %5.3f", NH₀.f)
                                Printf.@printf(" in test : %5d", v)
                            else
                                Printf.@printf "_"
                            end
                        end
                    end
                end
                #try to switch
                idxₜ₂ = findall(iszero.(NH₁.x[:, v]))
                idxₜ₂ = Random.shuffle!(idxₜ₂)#i₂ = 1, ..., I !removed
                betterFound = 0
                i₃ = 0
                x₋₁ = copy(NH₁.x)
                while betterFound == 0 && i₃ < size(idxₜ₂, 1)
                    i₃ += 1
                    i₂ = idxₜ₂[i₃]
                    if x_forced0ₜ[i₂] && i₂ != h
                        NH₁ = _mycopy(NH₀, NH₁)
                        NH₁.x = copy(x₋₁)
                        #NH₁.x[h, v] = zero(Float64)
                        NH₁.x[i₂, v] = one(Float64)
                        # iu = copy(iu₀)
                        # iu[i₂]+= 1
                        # iu[h]-= 1
                        if ata_model.settings.to_apply[1]
                            iu_max = sum(NH₁.x, dims = 2)[:, 1] - ata_model.settings.iu.max
                            iu_max = sum_pos(iu_max)
                        else
                            iu_max = 0
                        end
                        if ata_model.settings.to_apply[2]
                            iu_min = -sum(NH₁.x, dims = 2)[:, 1] + ata_model.settings.iu.min
                            iu_min = sum_pos(iu_min)
                        else
                            iu_min = 0
                        end
                        NH₁.iu = iu_max + iu_min
                        NH₁.infeas[v], x_Iₜ = check_feas(
                            ata_model.settings.fs,
                            Constraintsₜ,
                            NH₁.x[:, v],
                            n_fs,
                            n_items,
                        )
                        if fF == false
                            NH₁.obj[v] = eval_TIFₜ(x_Iₜ, coreₜ)
                        end
                        if ata_model.settings.to_apply[3]
                            NH₁.ol = eval_overlap(
                                NH₁.x,
                                fs_counts,
                                ata_model.settings.ol_max,
                                T,
                                NH₁.ol,
                            )
                        end
                        #NH₁.ol[v] = eval_overlapᵥ(NH₁.x[:, v], NH₁.x, ata_model.settings.fs.counts, ol_maxₜ, v)
                        NH₁.f = _comp_f(NH₁, opt_feas)
                        if (NH₁.f <= NH₀.f)
                            #switch item
                            #NH₀ = _mycopy(NH₁, NH₀)
                            # NH₀.ol = eval_overlap(NH₀.x, fs_counts, ata_model.settings.ol_max, T, NH₀.ol)
                            # NH₀.f = _comp_f(NH₀, opt_feas)
                            NH₀ = _mycopy(NH₁, NH₀)
                            # println("better to switch, new f₀ = ", NH₀.f)
                            if (NH₁.f < NH⁺.f)
                                exit = 1
                                betterFound = 1
                                convergence = 0
                                NH⁺ = _mycopy(NH₁, NH⁺)
                                if verbosity == 2
                                    Printf.@printf("\n f⁺ : %5.3f", NH⁺.f)
                                    Printf.@printf(" in test : %5d", v)
                                else
                                    Printf.@printf "+"
                                end
                            end
                        else
                            p = _sig_c(-(NH₁.f - NH₀.f) / t)
                            if p < 1e-10
                                p = zeros(Float64)[1]
                            end
                            if (rand() < p)
                                #remove
                                NH₀ = _mycopy(NH₁, NH₀)
                                # NH₀.ol = eval_overlap(NH₀.x, fs_counts, ata_model.settings.ol_max, T, NH₀.ol)
                                # NH₀.f = _comp_f(NH₀, opt_feas)
                                if verbosity == 2
                                    Printf.@printf("\n f₀ : %5.3f", NH₀.f)
                                    Printf.@printf(" in test : %5d", v)
                                else
                                    Printf.@printf "_"
                                end
                            end
                        end
                    end
                end #end of betterFound (betterFound = 1)
            end #end of itemorder h₂ (exit = 1)

        end #end of testorder v₂ (exit = 1)
        if sum(NH₀.infeas + NH₀.ol) + NH₀.iu <= 0 && fF == true
            fF = false
            for v = 1:T
                x_Iₜ = fs_to_items(
                    NH₀.x[:, v],
                    ata_model.settings.n_items,
                    ata_model.settings.fs.items,
                )
                NH₀.obj[v] = eval_TIFₜ(x_Iₜ, ata_model.obj.cores[v])
            end
            NH₀.f = _comp_f(NH₀, opt_feas)
            NH⁺ = _mycopy(NH₀, NH⁺)
        end
        f_star[1] = copy(NH₀.f)
        #if exit == 0
        if f_star[2] == f_star[1]
            convergence += 1
            Printf.@printf(" %3d", convergence)
        end
        #println("convergence is ", convergence)
        #how many equal f₀ in the last iterations?
        if (f_star[2] == f_star[1] && convergence == max_conv) ||
           time() - start_time >= max_time
            coverage_ok = 1
            #t += 1
        else
            t *= geom_temp
            coverage_ok = 0
            if f_star[1] <= f_star[2]
                f_star[2] = copy(f_star[1])
            end
        end
        if verbosity == 2
            println("\n")
            Printf.@printf("\n  f⁺:	%5.3f", NH⁺.f)
            Printf.@printf("\n  f₀:	%5.3f", NH₀.f)
            Printf.@printf("\n  Convergence:	%5d", convergence)
            println("\n")
        end
    end #end of NH coverage (coverage_ok = 1)
    print_neighbourhood(NH⁺)
    return NH⁺::Neighbourhood
end


text:
using JuliaMBD
using Test
using DifferentialEquations
using Plots

#=@testset "XmlToJulia" begin
    f = open("DCMtest.xml", "r")
    data = read(f, String)
    close(f)
    println(toJulia(data))
end

@testset "JuliaSimulate" begin
    @model MSD begin
        @parameter M g k D

        @blk In = InBlock()
        @blk Consta = ConstantBlock(M*g)
        @blk Gain1 = GainBlock(D)
        @blk Gain2 = GainBlock(k)
        @blk Gain3 = GainBlock(1/M)
        @blk Inte = IntegratorBlock()
        @blk Inte1 = IntegratorBlock(initialcondition=M*g/k)
        @blk Out = OutBlock()

        @connect In + Consta - Gain1 - Gain2 => Gain3
        @connect Gain3 => Inte
        @connect Inte => Inte1
        @connect Inte => Gain1
        @connect Inte1 => Out
        @connect Inte1 => Gain2
    end

    @model MSDtest begin
        @parameter M g k D ff

        @blk MSDsys = MSD(M=M, g=g, k=k, D=D) outport:x
        @blk Pulse = PulseGeneratorBlock(amplitude=ff, period=20, pulsewidth=50) outport:f

        @connect Pulse => MSDsys
    end

    m = MSDtest(M=10, g=9.8, k=100, D=10, ff=10)
    sol = @simulate(m, tspan=(0.0, 60.0), scope=(x, f))
    sol.graph
end

@testset "JuliaMBD.jl" begin
    f = open("DCMtest.xml", "r")
    data = read(f, String)
    close(f)
    eval(Meta.parse(toJulia(data)))

    @model DCM_test begin
        @parameter R K_e K_tau J_M J_I D L v_M

        @blk Step = StepBlock(steptime=1, finalvalue=v_M) outport:v_M
        @blk DCMmodel = DCM(R=R, K_e=K_e, K_tau=K_tau, J_M=J_M, J_I=J_I, D=D, L=L) outport[1]:omega outport[2]:i_M

        @connect Step => DCMmodel
    end

    m = DCM_test(R=5.7, K_e=7.16e-2, K_tau=7.2e-2, J_M=0.11e-3, J_I=1.3e-3, D=6.0e-5, L=0.2, v_M=24)
    sol = @simulate(m, tspan=(0.0, 10.0), scope=(omega, i_M, v_M))
    sol.graph

end

@testset "JuliaMBD_macro.jl" begin
    @xmlmodel "DCMtest.xml"

    @model DCM_test begin
        @parameter R K_e K_tau J_M J_I D L v_M

        @blk Step = StepBlock(steptime=1, finalvalue=v_M) outport:v_M
        @blk DCMmodel = DCM(R=R, K_e=K_e, K_tau=K_tau, J_M=J_M, J_I=J_I, D=D, L=L) outport[1]:omega outport[2]:i_M

        @connect Step => DCMmodel
    end

    m = DCM_test(R=5.7, K_e=7.16e-2, K_tau=7.2e-2, J_M=0.11e-3, J_I=1.3e-3, D=6.0e-5, L=0.2, v_M=24)
    sol = @simulate(m, tspan=(0.0, 10.0), scope=(omega, i_M, v_M))
    sol.graph

end

@testset "JuliaMBD_mxfile.jl" begin
    @xmlmodel "DCMexample.xml"

    @model DCM_test begin
        @parameter R K_e K_tau J_M J_I D L v_M

        @blk Step = StepBlock(steptime=1, finalvalue=v_M) outport:v_M
        @blk DCMmodel = DCM(R=R, K_e=K_e, K_tau=K_tau, J_M=J_M, J_I=J_I, D=D, L=L) outport[1]:omega outport[2]:i_M

        @connect Step => DCMmodel
    end

    m = DCM_test(R=5.7, K_e=7.16e-2, K_tau=7.2e-2, J_M=0.11e-3, J_I=1.3e-3, D=6.0e-5, L=0.2, v_M=24)
    sol = @simulate(m, tspan=(0.0, 10.0), scope=(omega, i_M, v_M))
    sol.graph

end

@testset "JuliaMBD_mxfile_intId.jl" begin
    @xmlmodel "DCMexample2.xml"

    @model DCM_test begin
        @parameter R K_e K_tau J_M J_I D L v_M

        @blk Step = StepBlock(steptime=1, finalvalue=v_M) outport:v_M
        @blk DCMmodel = DCM(R=R, K_e=K_e, K_tau=K_tau, J_M=J_M, J_I=J_I, D=D, L=L) outport[1]:omega outport[2]:i_M

        @connect Step => DCMmodel
    end

    m = DCM_test(R=5.7, K_e=7.16e-2, K_tau=7.2e-2, J_M=0.11e-3, J_I=1.3e-3, D=6.0e-5, L=0.2, v_M=24)
    sol = @simulate(m, tspan=(0.0, 10.0), scope=(omega, i_M, v_M))
    sol.graph

end

@testset "JuliaMBD_mxfile_3.jl" begin
    @xmlmodel "DCMexample3.xml"

    @model DCM_test begin
        @parameter R K_e K_tau J_M J_I D L v_M

        @blk Step = StepBlock(steptime=1, finalvalue=v_M) outport:v_M
        @blk DCMmodel = DCM(R=R, K_e=K_e, K_tau=K_tau, J_M=J_M, J_I=J_I, D=D, L=L) outport[1]:omega outport[2]:i_M

        @connect Step => DCMmodel
    end

    m = DCM_test(R=5.7, K_e=7.16e-2, K_tau=7.2e-2, J_M=0.11e-3, J_I=1.3e-3, D=6.0e-5, L=0.2, v_M=24)
    sol = @simulate(m, tspan=(0.0, 10.0), scope=(omega, i_M, v_M))
    sol.graph

end

@testset "JuliaMBD_compress.jl" begin
    @xmlmodel "DCMexample4.xml"

    @model DCM_test begin
        @parameter R K_e K_tau J_M J_I D L v_M

        @blk Step = StepBlock(steptime=1, finalvalue=v_M) outport:v_M
        @blk DCMmodel = DCM(R=R, K_e=K_e, K_tau=K_tau, J_M=J_M, J_I=J_I, D=D, L=L) outport[1]:omega outport[2]:i_M

        @connect Step => DCMmodel
    end

    m = DCM_test(R=5.7, K_e=7.16e-2, K_tau=7.2e-2, J_M=0.11e-3, J_I=1.3e-3, D=6.0e-5, L=0.2, v_M=24)
    sol = @simulate(m, tspan=(0.0, 10.0), scope=(omega, i_M, v_M))
    sol.graph

end

@testset "suion.jl" begin
    @xmlmodel "suion.xml"
end

@testset "ProductBlockTest" begin
    @model MotorDriver begin
        @parameter alpha_i

        @blk In1 = InBlock()
        @blk In2 = InBlock()
        @blk In3 = InBlock()
        @blk Gain = GainBlock(1/100)
        @blk Gain1 = GainBlock(alpha_i)
        #@blk Add = AddBlock([:+, :*]) inport[1]:in1 inport[2]:in2
        @blk Product = ProductBlock() inport[1]:in1 inport[2]:in2
        @blk Out = OutBlock()
        @blk Out1 = OutBlock()

        @connect In1 => in1
        @connect In2 => Gain
        @connect Gain => in2
        @connect Product => Out
        @connect In3 => Gain1
        @connect Gain1 => Out1
    end

    @model MDTest begin
        @parameter alpha_i u_M_d i_M_d Vs

        @blk Con = ConstantBlock(Vs)
        @blk Ramp = RampBlock(slope=u_M_d, starttime=0, initialoutput=0)
        @blk Ramp1 = RampBlock(slope=i_M_d, starttime=0, initialoutput=0)
        @blk MD = MotorDriver(alpha_i=alpha_i) inport[1]:Vs inport[2]:u_M inport[3]:i_M outport[1]:v_A outport[2]:v_i

        @connect Con => Vs
        @connect Ramp => u_M
        @connect Ramp1 => i_M
    end
    m = MDTest(alpha_i=5/37.5, u_M_d=10, i_M_d=3, Vs=24)
    sol = @simulate(m, tspan=(0.0, 10.0), scope=(v_A, Ramp, v_i, Ramp1))
    sol.graph
end

@testset "MotorDriverProductBlock.jl" begin
    @xmlmodel "MotorDriverP.xml"
end

@testset "SaturationTest.jl" begin
    @xmlmodel "SaturationTest.xml"
end

@testset "QuantizerTest.jl" begin
    f = open("QuantizerTest.xml", "r")
    data = read(f, String)
    close(f)
    println(toJulia(data))
    @xmlmodel "QuantizerTest.xml"

    @model ADConverterTest begin
        @parameter V_A_min V_A_max alpha_A

        @blk Ramp = RampBlock(slope=1, starttime=0, initialoutput=5)
        @blk ADC = ADConverter(V_A_min=V_A_min, V_A_max=V_A_max, alpha_A=alpha_A)

        @connect Ramp => ADC
    end
    m = ADConverterTest(V_A_min=0, V_A_max=0, alpha_A=(2^10 - 1)/5)
    sol = @simulate(m, tspan=(0.0, 10.0), scope=(Ramp, ADC))
    sol.graph
end

@testset "ModTest.jl" begin
    @xmlmodel "ModTest.xml"

    @model PulseGeneratorTest begin
        @parameter m alpha_P d_P_d

        @blk Ramp = RampBlock(slope=d_P_d, starttime=0, initialoutput=-50)
        @blk PG = PulseGenerator(m=m, alpha_P=alpha_P)

        @connect Ramp => PG
    end
    m = PulseGeneratorTest(m=8, alpha_P=100/(2^8-1), d_P_d=5)
    sol = @simulate(m, tspan=(0.0, 70.0), scope=(Ramp, PG))
    sol.graph
end

@testset "testADCModel" begin
    f = open("testADC.xml", "r")
    data = read(f, String)
    close(f)
    println(toJulia(data))
end

@testset "testADC.xml" begin
    @xmlmodel "testADC.xml"

    @model ADConverterTest begin
        @parameter V_A_max V_A_min alpha_A

        @blk Ramp = RampBlock(slope=1, starttime=0, initialoutput=-5)
        @blk ADC = newADC(V_A_max=V_A_max, V_A_min=V_A_min, alpha_A=alpha_A)

        @connect Ramp => ADC
    end
    m = ADConverterTest(V_A_min=0, V_A_max=5, alpha_A=(2^10 - 1)/5)
    sol = @simulate(m, tspan=(0.0, 15.0), scope=(Ramp, ADC))
    sol.graph
end

@testset "MDtest" begin
    f = open("MDtest.xml", "r")
    data = read(f, String)
    close(f)
    println(toJulia(data))

    @xmlmodel "MDtest.xml"

    @model MDTest begin
        @parameter alpha_i u_M_d i_M_d Vs

        @blk Con = ConstantBlock(Vs)
        @blk Ramp = RampBlock(slope=u_M_d, starttime=0, initialoutput=0)
        @blk Ramp1 = RampBlock(slope=i_M_d, starttime=0, initialoutput=0)
        @blk MD = MotorDriver(alpha_i=alpha_i) inport[1]:Vs inport[2]:u_M inport[3]:i_M outport[1]:v_A outport[2]:v_i

        @connect Con => Vs
        @connect Ramp => u_M
        @connect Ramp1 => i_M
    end
    m = MDTest(alpha_i=5/37.5, u_M_d=10, i_M_d=3, Vs=24)
    sol = @simulate(m, tspan=(0.0, 10.0), scope=(v_A, Ramp, v_i, Ramp1))
    sol.graph
end

@testset "MDtestNumber" begin
    f = open("MDtestNumber.xml", "r")
    data = read(f, String)
    close(f)
    println(toJulia(data))

    @xmlmodel "MDtestNumber.xml"
end=#

#=@testset "MDtest2" begin
    #=f = open("MDtest2.xml", "r")
    data = read(f, String)
    close(f)
    println(toJulia(data))
    =#
    @modelcontents "MDtest2.xml"

    @xmlmodel "MDtest2.xml"

    @model MDTest begin
        @parameter alpha_i u_M_d i_M_d Vs

        @blk Con = ConstantBlock(Vs)
        @blk Ramp = RampBlock(slope=u_M_d, starttime=0, initialoutput=0)
        @blk Ramp1 = RampBlock(slope=i_M_d, starttime=0, initialoutput=0)
        @blk MD = MotorDriver(alpha_i=alpha_i) inport[1]:Vs inport[2]:u_M inport[3]:i_M outport[1]:v_A outport[2]:v_i

        @connect Con => Vs
        @connect Ramp => u_M
        @connect Ramp1 => i_M
    end
    m = MDTest(alpha_i=5/37.5, u_M_d=10, i_M_d=3, Vs=24)
    sol = @simulate(m, tspan=(0.0, 10.0), scope=(v_A, Ramp, v_i, Ramp1))
    sol.graph
end

@testset "DCMD" begin
    f = open("DCMotorDisk.xml", "r")
    data = read(f, String)
    close(f)
    println(toJulia(data))

    @xmlmodel "DCMotorDisk.xml"
end=#

#=@testset "MSDModel" begin
    @modelcontents "MSDModel.xml"

    @xmlmodel "MSDModel.xml"

    @model MSDTest begin
        @parameter M g k D p_cycle f p_width

        @blk Pulse = PulseGeneratorBlock(amplitude=f, period=p_cycle, pulsewidth=p_width, phasedelay=10)
        @blk MSDModel = MSD(M=M, g=g, k=k, D=D)

        @connect Pulse => MSDModel
    end

    m = MSDTest(M=10, g=9.8, k=100, D=10, p_cycle=20, f=10, p_width=50)
    sol = @simulate(m, tspan=(0.0, 60.0), scope=(MSDModel))
    sol.graph
end=#

#=@testset "TempModel" begin
    @modelcontents "TempModel.xml"

    @xmlmodel "TempModel.xml"

    @model TempTest begin
        @parameter C R d0 q1

        @blk Step = StepBlock(steptime=100, initialvalue=0, finalvalue=q1)
        @blk TempModel = Temp(C=C, R=R, d0=d0)

        @connect Step => TempModel
    end

    m = TempTest(C=4187, R=0.1, d0=20, q1=300)
    sol = @simulate(m, tspan=(0.0, 1800), scope=(Step, TempModel))
    sol.graph
end

@testset "systemtest" begin
    @modelcontents "systemtest.xml"
end

@testset "systemtest2" begin
    @modelcontents "systemtest2.xml"
end=#

#=@testset "MDTest" begin
    @xmlmodel "MDtest2.xml"

    @xmlmodel "MDTest.xml"

    @model MDTest2 begin
        @parameter alpha_i u_M_d i_M_d Vs

        @blk Ramp = RampBlock(slope=u_M_d, starttime=0, initialoutput=0)
        @blk Ramp1 = RampBlock(slope=i_M_d, starttime=0, initialoutput=0)
        @blk MD = MDTest(alpha_i=alpha_i, Vs=Vs) inport[1]:u_M inport[2]:i_M outport[1]:v_A outport[2]:v_i

        @connect Ramp => u_M
        @connect Ramp1 => i_M
    end
    m = MDTest2(alpha_i=5/37.5, u_M_d=10, i_M_d=3, Vs=24)
    sol = @simulate(m, tspan=(0.0, 10.0), scope=(v_A, Ramp, v_i, Ramp1))
    sol.graph
end=#

@testset "inport, outport" begin
    @modelcontents "InOutTest.xml"
    @xmlmodel "InOutTest.xml"
    @model PulseGeneratorTest begin
        @parameter m alpha_P
        @parameter d_P_d

        @blk Ramp = RampBlock(slope=d_P_d, initialoutput=-50)
        @blk PG = PulseGenerator(m=m, alpha_P=alpha_P)

        @connect Ramp => PG
    end
    m = PulseGeneratorTest(m=8, alpha_P=100/(2^8-1), d_P_d=5)
    sol = @simulate(m, tspan=(0.0, 70.0), scope=(Ramp, PG))
    sol.graph
end

@testset "MotorDriver" begin
    @modelcontents "MotorDriver.xml"
    @xmlmodel "MotorDriver.xml"

    @modelcontents "MDriverTest.xml"
    @xmlmodel "MDriverTest.xml"

    m = MDriverTest(alpha_i=5/37.5, u_M_d=10, i_M_d=3, Vs=24)
    sol = @simulate(m, tspan=(0.0, 10.0), scope=(azRXO18pTX6HvRUZAOTaL3, Ramp, Ramp1))
    sol.graph
end

@testset "MotorDisk" begin
    @modelcontents "MotorDisk.xml"
    @xmlmodel "MotorDisk.xml"
end

@testset "Tachogenerator" begin
    @modelcontents "Tachogenerator.xml"
    @xmlmodel "Tachogenerator.xml"
end

@testset "Plant" begin
    @modelcontents "PlantModel.xml"
    @xmlmodel "PlantModel.xml"
end

@testset "Plant_test" begin
    @modelcontents "PlantModel_test.xml"
end


text:
# execution: julia -p 4 para_twindragon.jl

using Images
using PyPlot
using ProgressMeter
using LatexPrint

############
# settings #
############

# type
@everywhere const T = Float64

# number system setting
@everywhere const mat  = T[-1 -1 ; 1 -1]  #[-1 -1; 1 -1]
@everywhere const digs = T[[-1,-1] [1,-1]]

#magic constants
@everywhere const array_length = 500 #yet to be documented
@everywhere const radius = T(20)

#bounds
@everywhere const xmin= T(-1.7)
const xmax = T(1.4)

@everywhere const ymin= T(-1.3)
const ymax = T(2.2)

#sampling
@everywhere const dx = T(0.01)
@everywhere const dy = T(0.01)

#picture
const picture_title = "Twindragon"
const picture_dpi = 1000

#file
const filename = "twindragon_para" #warning, files will be overwritten without mercy (or questions)

###################
# end of settings #
###################

@everywhere include("Twindragon.jl")

const dimx = floor(Int, (xmax - xmin) / dx)
const dimy = floor(Int, (ymax - ymin) / dy)

image = SharedArray(T, dimx, dimy)

@everywhere oldset = SimpleSet()
@everywhere newset = SimpleSet()

@everywhere const mdigits = -1 * digs
@everywhere const number_of_digits = size(digs, 2)

@everywhere function compute_chunk!(image, jx, jy)
    @inbounds image[jx, jy] = twindragon([xmin + (jx - T(0.5))*dx, ymin + (jy - T(0.5))*dy], oldset, newset, ball)
end

jx = 1
jy = 1

# function to produce the next work item from the queue
function nextidx()
  global jx; global jy; global dimy; idx=(jx, jy); jy+=1;
  if jy > dimy
    jy = 1; jx+=1;
  end
  return idx
end

function run!(img)
  np = nprocs()  # determine the number of processes available
  progress_meter = Progress(dimx * dimy, 1)

  @sync begin
    for p in 1:np
      if p != myid() || np == 1
        @async begin
          while true
            idx = nextidx()
            if idx[1] > dimx
              break
            end
            remotecall_fetch(compute_chunk!, p, image, idx[1], idx[2])
            next!(progress_meter)
          end
        end
      end
    end
  end
end

println("Starting...")
tic()
run!(image)
toc()

# raw image
save("$(filename).png", grayim(transpose(flipdim(image,2))))

# pyplot insanity
PyPlot.matplotlib[:rc]("text", usetex=true)
PyPlot.matplotlib[:rc]("text.latex",preamble="\\usepackage{amsmath}") #to print matrices in latex
PyPlot.matplotlib[:rcParams]["text.latex.unicode"] = true #probably not needed
set_delims("(", ")") #set array delimiters for latex printing

imshow(transpose(flipdim(image,2)), cmap="Reds", interpolation="lanczos", extent=[xmin,xmax,ymin,ymax])

title(picture_title)

s = string("M = ",latex_form(mat),",D = \\left \\{ ",latex_form(digs[:,1]))
for d in 2:number_of_digits
    s = string(s,",",latex_form(digs[:,d]))
end
s = string(s," \\right \\} ")
text(xmin,ymin - (ymax-ymin)/10,latexstring(replace(s,"\n","")),fontsize=8)

savefig("$(filename)_pyplot.png", dpi = picture_dpi)

max_stars_repo_name: kalvotom/julia-examples

max_issues_repo_name: kalvotom/julia-examples

max_forks_repo_name: kalvotom/julia-examples


text:
#=
type jlDerivatives
   dinvdpos::Vector{mjtNum} #Matrix{mjtNum}
   dinvdvel::Vector{mjtNum} #Matrix{mjtNum}
   dinvdacc::Vector{mjtNum} #Matrix{mjtNum}
   daccdpos::Vector{mjtNum} #Matrix{mjtNum}
   daccdvel::Vector{mjtNum} #Matrix{mjtNum}
   daccdfrc::Vector{mjtNum} #Matrix{mjtNum}
end

type dinv
   dpos::Vector{mjtNum} #Matrix{mjtNum}
   dvel::Vector{mjtNum} #Matrix{mjtNum}
   dacc::Vector{mjtNum} #Matrix{mjtNum}
   dinv(nv::Integer) = new(Vector{mjtNum}(nv*nv),
                           Vector{mjtNum}(nv*nv),
                           Vector{mjtNum}(nv*nv))
end

type dacc
   dpos::Vector{mjtNum} #Matrix{mjtNum}
   dvel::Vector{mjtNum} #Matrix{mjtNum}
   dfrc::Vector{mjtNum} #Matrix{mjtNum}
   dacc(nv::Integer) = new(Vector{mjtNum}(nv*nv),
                           Vector{mjtNum}(nv*nv),
                           Vector{mjtNum}(nv*nv))
end
=#

function fwdworker(m::jlModel, dmain::jlData, d::jlData,
                   tid::Integer,
                   block1::Matrix{mjtNum}, block2::Matrix{mjtNum}, block3::Matrix{mjtNum})
   nv = length(d.qvel)
   eps = 1e-6
   nwarmup = 3
   nthread = Threads.nthreads() # TODO here?

   mark = mj.MARKSTACK(d)

   center = unsafe_wrap(Array, mj.stackAlloc(d.d, nv), nv)
   warmstart = unsafe_wrap(Array, mj.stackAlloc(d.d, nv), nv)
   angvel = unsafe_wrap(Array, mj.stackAlloc(d.d, 3), 3)
   quat = unsafe_wrap(Array, mj.stackAlloc(d.d, 4), 4)

   # prepare static schedule: range of derivative columns to be computed by this thread
   chunk = Integer(floor( (nv + nthread-1) / nthread ))
   istart = (tid-1) * chunk + 1
   iend = min(istart + chunk - 1, nv)
   println(istart, " ", chunk, " ", iend)

   # copy state and control from dmain to thread-specific d
   set(d, :time, get(dmain, :time)) #d.time = dmain.time
   copyto!(d.qpos, dmain.qpos)
   copyto!(d.qvel, dmain.qvel)
   copyto!(d.qacc, dmain.qacc)
   copyto!(d.qacc_warmstart, dmain.qacc_warmstart)
   copyto!(d.qfrc_applied, dmain.qfrc_applied)
   copyto!(d.xfrc_applied, dmain.xfrc_applied)
   copyto!(d.ctrl, dmain.ctrl)

   # run full computation at center point (usually faster than copying dmain)
   mj.forward(m, d)
   for rep=1:nwarmup
      mj.forwardSkip(m, d, Int(mj.STAGE_VEL), 1)
   end

   # save output for center point and warmstart (needed in forward only)
   copyto!(center, d.qacc)
   copyto!(warmstart, d.qacc_warmstart)

   # finite-difference over force or acceleration: skip = mjSTAGE_VEL
   for i=istart:iend
      # perturb selected target
      d.qfrc_applied[i] += eps

      # evaluate dynamics, with center warmstart
      copyto!(d.qacc_warmstart, warmstart)
      mj.forwardSkip(m, d, Int(mj.STAGE_VEL), 1)

      # undo perturbation
      d.qfrc_applied[i] = dmain.qfrc_applied[i]

      # compute column i of derivative 2
      block3[:,i] = (d.qacc - center)/eps
   end

   # finite-difference over velocity: skip = mjSTAGE_POS
   for i=istart:iend
      # perturb velocity
      d.qvel[i] += eps

      # evaluate dynamics, with center warmstart
      copyto!(d.qacc_warmstart, warmstart)
      mj.forwardSkip(m, d, Int(mj.STAGE_POS), 1)

      # undo perturbation
      d.qvel[i] = dmain.qvel[i]

      # compute column i of derivative 1
      block2[:,i] = (d.qacc - center)/eps
   end

   # finite-difference over position: skip = mjSTAGE_NONE
   for i=istart:iend
      # get joint id for this dof
      jid = m.dof_jntid[i] + 1

      # get quaternion address and dof position within quaternion (-1: not in quaternion)
      quatadr = -1
      dofpos = 0
      if m.jnt_type[jid] == Int(mj.JNT_BALL)
         quatadr = m.jnt_qposadr[jid]
         dofpos = (i-1) - m.jnt_dofadr[jid]
      elseif m.jnt_type[jid] == Int(mj.JNT_FREE) && i>=m.jnt_dofadr[jid]+4
         quatadr = m.jnt_qposadr[jid] + 3
         dofpos = (i-1) - m.jnt_dofadr[jid] - 3
      end

      # apply quaternion or simple perturbation
      if quatadr>=0
         angvel .= 0.0
         angvel[dofpos+1] = eps # already +1 from i
         quat .= d.qpos[(quatadr+1):(quatadr+4)]
         mju_quatIntegrate(quat, angvel, 1.0)
         d.qpos[(quatadr+1):(quatadr+4)] .= quat
      else
         d.qpos[i] += eps
      end

      # evaluate dynamics, with center warmstart
      copyto!(d.qacc_warmstart, warmstart)
      mj.forwardSkip(m, d, Int(mj.STAGE_NONE), 1)

      # undo perturbation
      copyto!(d.qpos, dmain.qpos)

      # compute column i of derivative 0
      block1[:,i] = (d.qacc - center)/eps
   end

   mj.FREESTACK(d, mark)
end

function invworker(m::jlModel, dmain::jlData, d::jlData,
                   tid::Integer,
                   block1::Matrix{mjtNum}, block2::Matrix{mjtNum}, block3::Matrix{mjtNum})

   nv = length(d.qvel)
   eps = 1e-6
   nwarmup = 3
   nthread = Threads.nthreads() # TODO here?

   mark = mj.MARKSTACK(d)

   center = unsafe_wrap(Array, mj.stackAlloc(d.d, nv), nv)
   warmstart = unsafe_wrap(Array, mj.stackAlloc(d.d, nv), nv)
   angvel = unsafe_wrap(Array, mj.stackAlloc(d.d, 3), 3)
   quat = unsafe_wrap(Array, mj.stackAlloc(d.d, 4), 4)

   # prepare static schedule: range of derivative columns to be computed by this thread
   chunk = Integer(floor( (nv + nthread-1) / nthread ))
   istart = (tid-1) * chunk + 1
   iend = min(istart + chunk - 1, nv)
   println(istart, " ", chunk, " ", iend)

   # copy state and control from dmain to thread-specific d
   set(d, :time, get(dmain, :time)) #d.time = dmain.time
   copyto!(d.qpos, dmain.qpos)
   copyto!(d.qvel, dmain.qvel)
   copyto!(d.qacc, dmain.qacc)
   copyto!(d.qacc_warmstart, dmain.qacc_warmstart)
   copyto!(d.qfrc_applied, dmain.qfrc_applied)
   copyto!(d.xfrc_applied, dmain.xfrc_applied)
   copyto!(d.ctrl, dmain.ctrl)

   # run full computation at center point (usually faster than copying dmain)
   mj.inverse(m, d)

   # save output for center point and warmstart (needed in forward only)
   copyto!(center, d.qfrc_inverse)
   copyto!(warmstart, d.qacc_warmstart)

   # finite-difference over force or acceleration: skip = mjSTAGE_VEL
   for i=istart:iend
      # perturb selected target
      d.qacc[i] += eps

      # evaluate dynamics, with center warmstart
      mj.inverseSkip(m, d, Int(mj.STAGE_VEL), 1)

      # undo perturbation
      d.qacc[i] = dmain.qacc[i]

      # compute column i of derivative 2
      block3[:,i] = (d.qfrc_inverse - center)/eps
   end

   # finite-difference over velocity: skip = mjSTAGE_POS
   for i=istart:iend
      # perturb velocity
      d.qvel[i] += eps

      # evaluate dynamics, with center warmstart
      mj.inverseSkip(m, d, Int(mj.STAGE_POS), 1)

      # undo perturbation
      d.qvel[i] = dmain.qvel[i]

      # compute column i of derivative 1
      block2[:,i] = (d.qfrc_inverse - center)/eps
   end

   # finite-difference over position: skip = mjSTAGE_NONE
   for i=istart:iend
      # get joint id for this dof
      jid = m.dof_jntid[i] + 1

      # get quaternion address and dof position within quaternion (-1: not in quaternion)
      quatadr = -1
      if m.jnt_type[jid] == Int(mj.JNT_BALL)
         quatadr = m.jnt_qposadr[jid]
         dofpos = (i-1) - m.jnt_dofadr[jid]
      elseif m.jnt_type[jid] == Int(mj.JNT_FREE) && i>=m.jnt_dofadr[jid]+4
         quatadr = m.jnt_qposadr[jid] + 3
         dofpos = (i-1) - m.jnt_dofadr[jid] - 3
      end

      # apply quaternion or simple perturbation
      if quatadr>=0
         angvel .= 0.0
         angvel[dofpos+1] = eps # already +1 from i
         quat .= d.qpos[(quatadr+1):(quatadr+4)]
         mju_quatIntegrate(quat, angvel, 1.0)
         d.qpos[(quatadr+1):(quatadr+4)] .= quat
      else
         d.qpos[i] += eps
      end

      # evaluate dynamics, with center warmstart
      mj.inverseSkip(m, d, Int(mj.STAGE_NONE), 1)

      # undo perturbation
      copyto!(d.qpos, dmain.qpos)

      # compute column i of derivative 0
      block1[:,i] = (d.qfrc_inverse - center)/eps
   end

   mj.FREESTACK(d, mark)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct LoadBalancerPropertiesFormat <: SwaggerModel
    frontendIPConfigurations::Any # spec type: Union{ Nothing, Vector{FrontendIPConfiguration} } # spec name: frontendIPConfigurations
    backendAddressPools::Any # spec type: Union{ Nothing, Vector{BackendAddressPool} } # spec name: backendAddressPools
    loadBalancingRules::Any # spec type: Union{ Nothing, Vector{LoadBalancingRule} } # spec name: loadBalancingRules
    probes::Any # spec type: Union{ Nothing, Vector{Probe} } # spec name: probes
    inboundNatRules::Any # spec type: Union{ Nothing, Vector{InboundNatRule} } # spec name: inboundNatRules
    inboundNatPools::Any # spec type: Union{ Nothing, Vector{InboundNatPool} } # spec name: inboundNatPools
    outboundRules::Any # spec type: Union{ Nothing, Vector{OutboundRule} } # spec name: outboundRules
    resourceGuid::Any # spec type: Union{ Nothing, String } # spec name: resourceGuid
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function LoadBalancerPropertiesFormat(;frontendIPConfigurations=nothing, backendAddressPools=nothing, loadBalancingRules=nothing, probes=nothing, inboundNatRules=nothing, inboundNatPools=nothing, outboundRules=nothing, resourceGuid=nothing, provisioningState=nothing)
        o = new()
        validate_property(LoadBalancerPropertiesFormat, Symbol("frontendIPConfigurations"), frontendIPConfigurations)
        setfield!(o, Symbol("frontendIPConfigurations"), frontendIPConfigurations)
        validate_property(LoadBalancerPropertiesFormat, Symbol("backendAddressPools"), backendAddressPools)
        setfield!(o, Symbol("backendAddressPools"), backendAddressPools)
        validate_property(LoadBalancerPropertiesFormat, Symbol("loadBalancingRules"), loadBalancingRules)
        setfield!(o, Symbol("loadBalancingRules"), loadBalancingRules)
        validate_property(LoadBalancerPropertiesFormat, Symbol("probes"), probes)
        setfield!(o, Symbol("probes"), probes)
        validate_property(LoadBalancerPropertiesFormat, Symbol("inboundNatRules"), inboundNatRules)
        setfield!(o, Symbol("inboundNatRules"), inboundNatRules)
        validate_property(LoadBalancerPropertiesFormat, Symbol("inboundNatPools"), inboundNatPools)
        setfield!(o, Symbol("inboundNatPools"), inboundNatPools)
        validate_property(LoadBalancerPropertiesFormat, Symbol("outboundRules"), outboundRules)
        setfield!(o, Symbol("outboundRules"), outboundRules)
        validate_property(LoadBalancerPropertiesFormat, Symbol("resourceGuid"), resourceGuid)
        setfield!(o, Symbol("resourceGuid"), resourceGuid)
        validate_property(LoadBalancerPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type LoadBalancerPropertiesFormat

const _property_map_LoadBalancerPropertiesFormat = Dict{Symbol,Symbol}(Symbol("frontendIPConfigurations")=>Symbol("frontendIPConfigurations"), Symbol("backendAddressPools")=>Symbol("backendAddressPools"), Symbol("loadBalancingRules")=>Symbol("loadBalancingRules"), Symbol("probes")=>Symbol("probes"), Symbol("inboundNatRules")=>Symbol("inboundNatRules"), Symbol("inboundNatPools")=>Symbol("inboundNatPools"), Symbol("outboundRules")=>Symbol("outboundRules"), Symbol("resourceGuid")=>Symbol("resourceGuid"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_LoadBalancerPropertiesFormat = Dict{Symbol,String}(Symbol("frontendIPConfigurations")=>"Vector{FrontendIPConfiguration}", Symbol("backendAddressPools")=>"Vector{BackendAddressPool}", Symbol("loadBalancingRules")=>"Vector{LoadBalancingRule}", Symbol("probes")=>"Vector{Probe}", Symbol("inboundNatRules")=>"Vector{InboundNatRule}", Symbol("inboundNatPools")=>"Vector{InboundNatPool}", Symbol("outboundRules")=>"Vector{OutboundRule}", Symbol("resourceGuid")=>"String", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ LoadBalancerPropertiesFormat }) = collect(keys(_property_map_LoadBalancerPropertiesFormat))
Swagger.property_type(::Type{ LoadBalancerPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_LoadBalancerPropertiesFormat[name]))}
Swagger.field_name(::Type{ LoadBalancerPropertiesFormat }, property_name::Symbol) =  _property_map_LoadBalancerPropertiesFormat[property_name]

function check_required(o::LoadBalancerPropertiesFormat)
    true
end

function validate_property(::Type{ LoadBalancerPropertiesFormat }, name::Symbol, val)
end


text:
## Fit counties with significant sero data

## Do inference in parallel

using Distributed

## add number of processors desired
# addprocs(2)
params = (exename = `nice -19 /Applications/Julia-1.7.app/Contents/Resources/julia/bin/julia`,
    dir = "./Github/KenyaCoVaccinesPrivate")
addprocs([("samandfi@10.0.0.6", 6)]; params...)
addprocs(6)
# params = (exename = `nice -19 /Applications/Julia-1.6.app/Contents/Resources/julia/bin/julia`,
#     dir = "./Github/KenyaCoVaccinesPrivate")
# addprocs([("samandfi@10.0.0.6", :auto)]; params...)
# addprocs(2)

@everywhere begin
    using Pkg
    Pkg.activate(".")
end

## Load relevant code into scope for each worker

@everywhere using JLD2, LinearAlgebra, TransformVariables, LogDensityProblems, MCMCChains, Dates
@everywhere import KenyaCoVaccines

## Load data
# @everywhere begin
using DataFrames
@load("data/linelist_data_with_pos_neg_20feb_to_27sept_c_age__pcr.jld2")
@load("data/cleaned_serology_data_c_age__combined.jld2")
@load("data/deaths_20210919.jld2")
@load("data/p_IH.jld2")
@load("data/p_ID.jld2")

#Load contact structure and population size data
M_Kenya_ho, M_Kenya_other, M_Kenya_school, M_Kenya_work = KenyaCoVaccines.get_rescaled_contact_matrices("data/agemixingmatrix_Kenya_all_types.jld2")
N_kenya = KenyaCoVaccines.get_population_size_matrix("data/2019_census_age_pyramids_counties.csv")

# Parameter transformations and initial KE guess for HMC

trans = as((
    β₀ = as(Real, 0.0, 5.0),  #β₀= as(Real, 0.0, 7.0),
    β_home = as(Real, 0.0, 5.0),
    β_school = as(Real, 0.0, 5.0),
    β_other = as(Real, 0.0, 5.0),
    β_work = as(Real, 0.0, 5.0),
    ϵ = as(Real, 0.0, 0.999),
    χ = as(Real, 0.1, 30.0),
    p_test = as(Real, 0.01, 5),
    E₀ = as(Real, 0.0, 10e3),
    inc_R_αvar = as(Real, 0.0, 1.0),
    time_scale_αvar = as(Real, 0.01, 0.5),
    mid_point_αvar = as(Real, 0.0, 4.0),
    inc_R_δvar = as(Real, 0.001, 1.5),
    time_scale_δvar = as(Real, 0.01, 0.5),
    mid_point_δvar = as(Real, 0.0, 4.0),
    init_scale = as(Real, 0.0, 2.0)
))

D = Diagonal(0.1 * ones(TransformVariables.dimension(trans)))

# end

## Name for each county with substantial serological data, and matching priors

semiurbancounties_lowserodata = ["Kajiado", "Kiambu", "Machakos", "Isiolo", "Kirinyaga", "Murang'a",
    "Trans Nzoia", "Nyandarua", "Meru", "Bungoma", "Nyamira", "Kakamega",
    "Laikipia", "Taita Taveta", "Makueni", "Elgeyo-Marakwet",
    "Tharaka-Nithi", "Kericho", "Bomet", "Homa Bay", "Busia", "Migori",
    "Nandi", "Vihiga"]

ruralcounties_lowserodata = ["Garissa", "Kitui",
    "Lamu", "Marsabit", "Narok", "Tana River", "Turkana", "Wajir",
    "Baringo", "Mandera", "Samburu", "West Pokot"]
counties_with_informed_priors = [semiurbancounties_lowserodata; ruralcounties_lowserodata]

# counties_with_informed_priors = ["Kiambu"
# # , "Machakos", "Murang'a", "Meru",
#     # "Kakamega", "Makueni", "Kericho", "Busia", "Vihiga"]

#Baseline seropos for low density coverage counties at December 1st
#Based on the unweighted seroprevalence from Ngere et al in Nairobi mid-Sept scaled down 33% -> 25% unweighted seroprevalence
init_sero_baseline = [(34 + 74) / (179 + 244), (100 + 95 + 51) / (265 + 241 + 134), 21 / 61, 9 / 40, 9 / 40, 9 / 40] .* (25 / 33)

## Parallelised loop for fitting each county

@distributed for i = 1:36
    name = counties_with_informed_priors[i]
    model = KenyaCoVaccines.CoVAreaModel(name, KenyaCoVaccines.ll_onegroup_newvariant_infboost, KenyaCoVaccines.priors_onegroup_alpha_delta_variant_noncities_fitted;
        linelist_data_with_pos_neg = linelist_data_with_pos_neg,
        serology_data = serology_data,
        average_sero_init = init_sero_baseline,
        deaths = deaths_data,
        pop_data = N_kenya,
        M_county_ho = M_Kenya_ho,
        M_county_other = M_Kenya_other,
        M_county_school = M_Kenya_school,
        M_county_work = M_Kenya_work,
        rapid = false,
        scenario = 0,
        startdate = Date(2020, 12, 1),
        enddate = Date(2021, 9, 24))

    model.init_sero = init_sero_baseline #Enforcing this because maybe more reliable with low data

    KenyaCoVaccines.inferparameters!(model, 2000, trans, 0.05, D; serowaningrate = 0.0)
    @save("modelfits/$(name)_model.jld2", model)
    println("Finished inference for county $(name)")

end
##
rmprocs(workers())


text:
# See https://github.com/spglib/spglib/blob/444e061/python/spglib/spglib.py#L415-L463 and https://github.com/unkcpz/LibSymspg.jl/blob/f342e72/src/cell-reduce-api.jl#L3-L35
"""
    standardize_cell(cell::Cell; to_primitive=false, no_idealize=false, symprec=1e-5)

Return standardized cell.

The standardized unit cell is generated from an input unit cell structure and
its symmetry found by the symmetry search. The choice of the setting for each
space group type is as explained for [`get_dataset`](@ref).
"""
function standardize_cell(
    cell::Cell;
    to_primitive = false,
    no_idealize = false,
    symprec = 1e-5,
)
    @unpack lattice, positions, types = _expand_cell(cell)
    to_primitive = Base.cconvert(Cint, to_primitive)
    no_idealize = Base.cconvert(Cint, no_idealize)
    num_atom = Base.cconvert(Cint, length(types))
    allocations = 4
    _positions = Matrix{Cdouble}(undef, 3, num_atom * allocations)
    _types = Vector{Cint}(undef, num_atom * allocations)
    _positions[:, 1:num_atom] = positions
    _types[1:num_atom] = types
    num_atom_std = ccall(
        (:spg_standardize_cell, libsymspg),
        Cint,
        (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint, Cint, Cint, Cdouble),
        lattice,
        _positions,
        _types,
        num_atom,
        to_primitive,
        no_idealize,
        symprec,
    )  # Note: not `num_atom`!
    if num_atom_std <= 0
        throw(SpglibError("Cell standardization failed!"))
    end
    # We have to `transpose` back because of `_expand_cell`!
    return Cell(transpose(lattice), _positions[:, 1:num_atom_std], _types[1:num_atom_std])
end

"""
    find_primitive(cell::Cell, symprec=1e-5)

Find the primitive cell of an input unit cell.

This function is now a shortcut of `standardize_cell` with `to_primitive = true`
and `no_idealize = false`.
"""
find_primitive(cell::Cell, symprec = 1e-5) =
    standardize_cell(cell; to_primitive = true, no_idealize = false, symprec = symprec)

"""
    refine_cell(cell::Cell, symprec=1e-5)

Return refined cell.

The standardized crystal structure is obtained from a non-standard crystal
structure which may be slightly distorted within a symmetry recognition
tolerance, or whose primitive vectors are differently chosen, etc.
This function is now a shortcut of `standardize_cell` with `to_primitive = false`
and `no_idealize = false`.
"""
refine_cell(cell::Cell, symprec = 1e-5) =
    standardize_cell(cell; to_primitive = false, no_idealize = false, symprec = symprec)


text:
"""
SubProblem

* `np`: a TransitNetwork object that contains the network information
* `model`: a JuMP model
* `src`: src[u] indicates whether the source is at stop u
* `snk`: snk[u] indicates whether the sink is at stop v
* `edg`: edg[u,v] indicates whether edge (u,v) is used
* `ingraph`: ingraph[u] is an expression indicating whether u is used in the line
* `srv`: srv[u,v] indicates whether commute (u,v) can be served by the line directly (one leg)
* `srv2`: srv2[u,v] indicates whether commute (u,v) can be served by the line with one transfer (two legs)
* `dists`: dists[u,v] is the distance of edge (u,v)
* `outneighbors`: outneighbors[u] is a vector containing other stations v for which (u, v) is a valid edge
* `inneighbors`: inneighbors[u] is a vector containing other stations v for which (v, u) is a valid edge
* `nlegs`: passed by the MasterProblem; 1 (direct-route model) or 2 (single-transfer model)
* `auxinfo`: stores auxiliary solver information
"""
mutable struct SubProblem
    np::TransitNetwork
    model::JuMP.Model
    src::Vector{JuMP.Variable}
    snk::Vector{JuMP.Variable}
    edg::JuMP.JuMPDict{JuMP.Variable}
    ingraph
    srv::JuMP.JuMPArray
    srv2
    dists::Dict{Tuple{Int,Int},Float64}
    outneighbors::Vector{Vector{Int}}
    inneighbors::Vector{Vector{Int}}
    nlegs::Int
    auxinfo::Dict{Symbol,Any}
end

"""
SubProblem with Edge Processing

Constructs an acyclic graph that only has edges within some `delta` tolerance of
`direction`.

### Keyword Arguments
* `nlegs`: the (maximum) number of legs of each commute (should be same as in the MasterProblem).
* `solver`: The solver being used to solve the problem.
* `maxdist`: the threshold distance for two stations to be considered.
    neighbors. The distance is measured based on the `gridtype` in
    `rmp::MasterProblem`.
* `direction`: The direction for which we run the subproblem with preprocessed
    edge set E(direction, delta).
* `delta`: tolerance within `direction` for edges in the graph to obey.
* `maxlength`: the maximum number of edges in a path.
* `traveltimes`: whether to constrain travel times on the network
* `detail`: whether to use a more detailed heuristic to compute shortest paths for travel times
"""
function SubProblem(
        rmp::MasterProblem;
        nlegs::Int = rmp.options.nlegs,
        solver = Gurobi.GurobiSolver(OutputFlag = 0),
        maxdist::Float64 = 0.5,
        direction::Vector{Float64} = [0.0,1.0],
        delta::Float64 = 1.0, # 1 - (delta in paper) -- 1.0 in the code corresponds to 0.0 in the paper
        maxlength::Int = 30, # maximum number of edges in a path
        traveltimes::Bool = false,
        detail::Bool = false
    )
    np = rmp.np
    nstns = np.nstations

    # construct graph and ensure there are no cycles
    if rmp.np.dists == nothing
        dists = Dict{Tuple{Int,Int},Float64}()
    else
        dists = rmp.np.dists
    end
    outneighbors = [Int[] for u in 1:nstns]
    inneighbors  = [Int[] for u in 1:nstns]
    graph = LightGraphs.DiGraph(nstns)
    for u in 1:nstns, v in (u+1):nstns
        d = edgecost(np, u, v)
        b = dir(np, u, v)
        sim = dot(direction, b) / norm(direction) / norm(b)
        if (d < maxdist) && (1 - abs(sim) <= delta)
            if sim >= 0
                if rmp.np.dists == nothing
                    dists[u,v] = d
                end
                push!(outneighbors[u], v)
                push!(inneighbors[v], u)
                LightGraphs.add_edge!(graph, (u, v))
            else
                if rmp.np.dists == nothing
                    dists[v,u] = d
                end
                push!(outneighbors[v], u)
                push!(inneighbors[u], v)
                LightGraphs.add_edge!(graph, (v, u))
            end
        end
    end
    @assert !LightGraphs.is_cyclic(graph)

    # construct base JuMP model
    sp, src, snk, edg, srv, ingraph = basemodel(
        rmp, inneighbors, outneighbors, maxlength, rmp.linelist, solver
    )
    # add variables and constraints for single-transfer model if relevant
    if nlegs == 2
        srv2 = transfermodel(rmp, sp, srv, ingraph)
    else
        srv2 = nothing
    end

    # initializes empty dictionary for auxiliary solve info to be stored
    auxinfo = initializedict()

    # adds cuts for limiting travel time on network
    function cuttraveltimes(cb)
        # get line from the variable values
        visited = falses(np.nstations) # whether a node has been visited
        visited[setdiff(1:np.nstations, findall(JuMP.getvalue(ingraph) .> 0))] .= true
        source_val = findfirst(round.(JuMP.getvalue.(src)) .> 0.1)
        sink_val = findfirst(round.(JuMP.getvalue.(snk)) .> 0.1)
        @assert source_val > 0
        @assert sink_val > 0
        @assert source_val != sink_val
        simplepath = getpath( # bus line
            source_val, source_val, sink_val, edg, visited, outneighbors
        )
        # add the cut
        if traveltimes
            addlazytraveltimes(cb,
                np, ingraph, edg, outneighbors, rmp.options.distparam,
                detail,
                auxinfo,
                simplepath)
        end
    end
    if traveltimes
        JuMP.addlazycallback(sp, cuttraveltimes)
    end

    SubProblem(
        np, sp, src, snk, edg, ingraph, srv, srv2,
        dists, outneighbors, inneighbors, nlegs, auxinfo
    )
end

"""
SubProblem with Cutting Planes

Constructs a graph (not necessarily acyclic) to be optimized over by using
cutting planes to remove subtours.

### Keyword Arguments
* `nodeset`: The set of stations under consideration (default all).
* `nlegs`: the (maximum) number of legs of each commute (should be same as in the MasterProblem).
* `solver`: The solver being used to solve the problem.
* `maxdist`: the threshold distance for two stations to be considered.
    neighbors. The distance is measured based on the `gridtype` in
    `rmp::MasterProblem`.
* `maxlength`: the maximum number of edges in a path.
* `traveltimes`: whether to constrain travel times on the network
* `detail`: whether to use a more detailed heuristic to compute shortest paths for travel times
"""
function SubProblemCP(
        rmp::MasterProblem;
        nodeset::Vector{Int} = Vector(1:rmp.np.nstations),
        nlegs::Int = rmp.options.nlegs,
        solver = Gurobi.GurobiSolver(OutputFlag = 0),
        maxdist::Float64 = 0.5,
        maxlength::Int = 30,
        traveltimes::Bool = false,
        detail::Bool = false
    )
    np = rmp.np
    nstns = np.nstations

    # construct graph and ensure there are no cycles
    if rmp.np.dists == nothing
        dists = Dict{Tuple{Int,Int},Float64}()
    else
        dists = rmp.np.dists
    end
    outneighbors = [Int[] for u in 1:nstns]
    inneighbors  = [Int[] for u in 1:nstns]
    for u in 1:nstns
        for v in (u+1):nstns
            !in(u, nodeset) && !in(v, nodeset) && continue
            d = edgecost(np, u, v)
            if d < maxdist
                if rmp.np.dists == nothing
                    dists[u,v] = d
                    dists[v,u] = d
                end
                push!(outneighbors[u], v)
                push!(inneighbors[v], u)
                push!(outneighbors[v], u)
                push!(inneighbors[u], v)
            end
        end
    end

    # construct base JuMP model
    sp, src, snk, edg, srv, ingraph = basemodel(
        rmp, inneighbors, outneighbors, maxlength, rmp.linelist, solver
    )
    # add variables and constraints for single-transfer model if relevant
    if nlegs == 2
        srv2 = transfermodel(rmp, sp, srv, ingraph)
    else
        srv2 = nothing
    end
    # some easy cuts to start
    JuMP.@constraint(sp,
        [u=1:np.nstations,v=outneighbors[u]],
        edg[u,v] + edg[v,u] <= 1)

    # initializes empty dictionary for auxiliary solve info to be stored
    auxinfo = initializedict()

    # subtour elimination
    function removecycles(cb)
        visited = falses(np.nstations) # whether a node has been visited
        visited[setdiff(1:np.nstations, findall(JuMP.getvalue(ingraph) .> 0))] .= true
        source_val = findfirst(round.(JuMP.getvalue.(src)) .> 0.1)
        sink_val = findfirst(round.(JuMP.getvalue.(snk)) .> 0.1)
        @assert source_val > 0
        @assert sink_val > 0
        @assert source_val != sink_val
        simplepath = getpath( # bus line
            source_val, source_val, sink_val, edg, visited, outneighbors
        )
        while length(findall(.!visited)) > 0 # search for subtours
            cyclenodes = getpath(
                findfirst(.!visited), source_val, sink_val, edg, visited,
                outneighbors
            )
            if length(cyclenodes) > 1
                if traveltimes && !checknodes(np, rmp.options.distparam, cyclenodes)
                    expr = sum(ingraph[u] for u in cyclenodes)
                else
                    expr = sum(sum(edg[u,v]
                                   for v in intersect(outneighbors[u], cyclenodes))
                               for u in cyclenodes)
                end
                JuMP.@lazyconstraint(cb, expr <= length(cyclenodes) - 1)
                auxinfo[:cycle] += 1
                auxinfo[:nlazy] += 1
                break
            end
        end
        if traveltimes
            addlazytraveltimes(cb,
                np, ingraph, edg, outneighbors, rmp.options.distparam,
                detail,
                auxinfo,
                simplepath)
        end
    end
    JuMP.addlazycallback(sp, removecycles)

    SubProblem(
        np, sp, src, snk, edg, ingraph, srv, srv2,
        dists, outneighbors, inneighbors, nlegs, auxinfo
    )
end

"""
Creates an empty Dict saving solve information
"""
function initializedict()
    auxinfo = Dict{Symbol,Any}()
    auxinfo[:nlazy] = 0
    auxinfo[:travelshort] = 0
    auxinfo[:traveltourn] = 0
    auxinfo[:cycle] = 0
    auxinfo[:checknodes] = Dict{Vector{Int},Bool}()
    auxinfo
end