text:
# Obtain AMR indicators from solver, then adapt mesh, and finally adapt solver
#
# If `only_refine` is true, no coarsening will be performed, independent of the indicator values.
# If `only_coarsen` is true, no refinement will be performed, independent of the indicator values.
# Passively adapted solvers in `passive_solvers` will be adapated according to the primary solver.
#
# Return true if anything was changed, false if no cells where coarsened/refined
function adapt!(mesh::TreeMesh, solver::AbstractSolver, time;
                only_refine=false, only_coarsen=false, passive_solvers=[])
  # Debug output
  globals[:verbose] && print("Begin adaptation...")

  # Alias for convenience
  tree = mesh.tree

  # Determine indicator value
  lambda = @timeit timer() "indicator" calc_amr_indicator(solver, mesh, time)

  # Get list of current leaf cells
  leaf_cell_ids = leaf_cells(tree)
  @assert length(lambda) == length(leaf_cell_ids) ("Indicator and leaf cell arrays have " *
                                                   "different length")

  # Set thresholds for refinement and coarsening
  indicator_threshold_refinement = parameter("indicator_threshold_refinement",  0.5)
  indicator_threshold_coarsening = parameter("indicator_threshold_coarsening", -0.5)

  # Determine list of cells to refine or coarsen
  to_refine = leaf_cell_ids[lambda .> indicator_threshold_refinement]
  to_coarsen = leaf_cell_ids[lambda .< indicator_threshold_coarsening]

  # Start by refining cells
  @timeit timer() "refine" if !only_coarsen && !isempty(to_refine)
    # Refine mesh
    refined_original_cells = @timeit timer() "mesh" refine!(tree, to_refine)

    # Refine solver
    @timeit timer() "solver" refine!(solver, mesh, refined_original_cells)
    if !isempty(passive_solvers)
      @timeit timer() "passive solvers" for ps in passive_solvers
        refine!(ps, mesh, refined_original_cells)
      end
    end
  else
    # If there is nothing to refine, create empty array for later use
    refined_original_cells = Int[]
  end

  # Then, coarsen cells
  @timeit timer() "coarsen" if !only_refine && !isempty(to_coarsen)
    # Since the cells may have been shifted due to refinement, first we need to
    # translate the old cell ids to the new cell ids
    if !isempty(to_coarsen)
      to_coarsen = original2refined(to_coarsen, refined_original_cells, mesh)
    end

    # Next, determine the parent cells from which the fine cells are to be
    # removed, since these are needed for the coarsen! function. However, since
    # we only want to coarsen if *all* child cells are marked for coarsening,
    # we count the coarsening indicators for each parent cell and only coarsen
    # if all children are marked as such (i.e., where the count is 2^ndims). At
    # the same time, check if a cell is marked for coarsening even though it is
    # *not* a leaf cell -> this can only happen if it was refined due to 2:1
    # smoothing during the preceding refinement operation.
    parents_to_coarsen = zeros(Int, length(tree))
    for cell_id in to_coarsen
      # If cell has no parent, it cannot be coarsened
      if !has_parent(tree, cell_id)
        continue
      end

      # If cell is not leaf (anymore), it cannot be coarsened
      if !is_leaf(tree, cell_id)
        continue
      end

      # Increase count for parent cell
      parent_id = tree.parent_ids[cell_id]
      parents_to_coarsen[parent_id] += 1
    end

    # Extract only those parent cells for which all children should be coarsened
    to_coarsen = collect(1:length(parents_to_coarsen))[parents_to_coarsen .== 2^ndims(mesh)]

    # Finally, coarsen mesh
    coarsened_original_cells = @timeit timer() "mesh" coarsen!(tree, to_coarsen)

    # Convert coarsened parent cell ids to the list of child cell ids that have
    # been removed, since this is the information that is expected by the solver
    removed_child_cells = zeros(Int, n_children_per_cell(tree) * length(coarsened_original_cells))
    for (index, coarse_cell_id) in enumerate(coarsened_original_cells)
      for child in 1:n_children_per_cell(tree)
        removed_child_cells[n_children_per_cell(tree) * (index-1) + child] = coarse_cell_id + child
      end
    end

    # Coarsen solver
    @timeit timer() "solver" coarsen!(solver, mesh, removed_child_cells)
    if !isempty(passive_solvers)
      @timeit timer() "passive solvers" for ps in passive_solvers
        coarsen!(ps, mesh, removed_child_cells)
      end
    end
  else
    # If there is nothing to coarsen, create empty array for later use
    coarsened_original_cells = Int[]
  end

  # Debug output
  globals[:verbose] && println("done (refined: $(length(refined_original_cells)), " *
                               "coarsened: $(length(coarsened_original_cells)), " *
                               "new number of cells/elements: " *
                               "$(length(tree))/$(solver.n_elements))")

  # Return true if there were any cells coarsened or refined, otherwise false
  return !isempty(refined_original_cells) || !isempty(coarsened_original_cells)
end

# After refining cells, shift original cell ids to match new locations
# Note: Assumes sorted lists of original and refined cell ids!
# Note: `mesh` is only required to extract ndims
function original2refined(original_cell_ids, refined_original_cells, mesh)
  # Sanity check
  @assert issorted(original_cell_ids) "`original_cell_ids` not sorted"
  @assert issorted(refined_original_cells) "`refined_cell_ids` not sorted"

  # Create array with original cell ids (not yet shifted)
  shifted_cell_ids = collect(1:original_cell_ids[end])

  # Loop over refined original cells and apply shift for all following cells
  for cell_id in refined_original_cells
    # Only calculate shifts for cell ids that are relevant
    if cell_id > length(shifted_cell_ids)
      break
    end

    # Shift all subsequent cells by 2^ndims ids
    shifted_cell_ids[(cell_id + 1):end] .+= 2^ndims(mesh)
  end

  # Convert original cell ids to their shifted values
  return shifted_cell_ids[original_cell_ids]
end


text:
using Test
using MultivariatePolynomials

# Allocating size for allocating a `BigInt`.
# Half size on 32-bit.
const BIGINT_ALLOC = Sys.WORD_SIZE == 64 ? 48 : 24

function alloc_test(f, n)
    f() # compile
    @test n == @allocated f()
end

struct CustomPoly{T, P<:AbstractPolynomial{T}} <: AbstractPolynomialLike{T}
    p::P
end
CustomPoly(p::AbstractPolynomial{T}) where T = CustomPoly{T, typeof(p)}(p)
MultivariatePolynomials.polynomial(p::CustomPoly) = p.p
MultivariatePolynomials.polynomial(p::CustomPoly, T::Type) = polynomial(p.p, T)
MultivariatePolynomials.variables(p::CustomPoly) = variables(p.p)
MultivariatePolynomials.monomialtype(::Type{<:CustomPoly{T, P}}) where {T, P} = monomialtype(P)
MultivariatePolynomials.constantmonomial(p::CustomPoly) = constantmonomial(p.p)

struct CustomTerms{T, P<:AbstractPolynomial{T}} <: AbstractPolynomialLike{T}
    p::P
end
CustomTerms(p::AbstractPolynomial{T}) where T = CustomTerms{T, typeof(p)}(p)
MultivariatePolynomials.terms(p::CustomTerms) = terms(p.p)
MultivariatePolynomials.variables(p::CustomTerms) = variables(p.p)
MultivariatePolynomials.monomialtype(::Type{<:CustomTerms{T, P}}) where {T, P} = monomialtype(P)
MultivariatePolynomials.constantmonomial(p::CustomPoly) = constantmonomial(p.p)

function _typetests(x, ::Type{T}) where T
    @test (@inferred coefficienttype(x)) == Int

    @test (@inferred monomialtype(x)) <: AbstractMonomial

    @test (@inferred termtype(x)) <: AbstractTerm{Int}
    @test (@inferred termtype(x, Float64)) <: AbstractTerm{Float64}

    @test (@inferred polynomialtype(x)) <: AbstractPolynomial{Int}
    @test (@inferred polynomialtype(x, Float64)) <: AbstractPolynomial{Float64}

    @test (@inferred monovectype(x)) <: AbstractArray{<:AbstractMonomial}
end

function typetests(x::Union{AbstractPolynomialLike{T}, Vector{<:AbstractPolynomialLike{T}}}) where T
    _typetests(x, T)
    _typetests(typeof(x), T)
end

ext: jl


text:
using Test, Sandbox

all_executors = Sandbox.all_executors

# Can we run `sudo` without a password?  If not, don't attempt to test the privileged runner
if !success(`sudo -k -n true`)
    all_executors = filter(exe -> exe != PrivilegedUserNamespacesExecutor, all_executors)
end

rootfs_dir = artifact"AlpineRootfs"
for executor in all_executors
    if !executor_available(executor)
        @error("Skipping $(executor) tests, as it does not seem to be available")
        continue
    end

    @testset "$(executor) Sandboxing" begin
        @testset "capturing stdout/stderr" begin
            stdout = IOBuffer()
            stderr = IOBuffer()
            config = SandboxConfig(
                Dict("/" => rootfs_dir);
                stdout,
                stderr,
            )
            with_executor(executor) do exe
                @test run(exe, config, `/bin/sh -c "echo stdout; echo stderr >&2"`)
                @test String(take!(stdout)) == "stdout\n";
                @test String(take!(stderr)) == "stderr\n";
            end
        end

        @testset "ignorestatus()" begin
            config = SandboxConfig(Dict("/" => rootfs_dir))
            with_executor(executor) do exe
                @test_throws ProcessFailedException run(exe, config, `/bin/sh -c "false"`)
                @test !run(exe, config, ignorestatus(`/bin/sh -c "false"`))
            end
        end

        @testset "environment passing" begin
            # Ensure all those pesky "special" variables make it through
            env = Dict(
                "PATH" => "for",
                "LD_LIBRARY_PATH" => "science",
                "DYLD_LIBRARY_PATH" => "you",
                "SHELL" => "monster",
            )
            stdout = IOBuffer()
            config = SandboxConfig(
                Dict("/" => rootfs_dir),
                Dict{String,String}(),
                env;
                stdout,
            )
            user_cmd = `/bin/sh -c "echo \$PATH \$LD_LIBRARY_PATH \$DYLD_LIBRARY_PATH \$SHELL"`
            with_executor(executor) do exe
                @test run(exe, config, user_cmd)
                @test String(take!(stdout)) == "for science you monster\n";
            end

            # Test that setting some environment onto `user_cmd` can override the `config` env:
            user_cmd = setenv(user_cmd, "DYLD_LIBRARY_PATH" => "my", "SHELL" => "friend")
            with_executor(executor) do exe
                @test run(exe, config, user_cmd)
                @test String(take!(stdout)) == "for science my friend\n";
            end
        end

        @testset "reading from maps" begin
            mktempdir() do dir
                open(joinpath(dir, "note.txt"), write=true) do io
                    write(io, "great success")
                end
                stdout = IOBuffer()
                config = SandboxConfig(
                    Dict("/" => rootfs_dir, "/glados" => dir);
                    stdout,
                )
                with_executor(executor) do exe
                    @test run(exe, config, `/bin/sh -c "cat /glados/note.txt"`)
                    @test String(take!(stdout)) == "great success";
                end
            end
        end

        @testset "writing to workspaces" begin
            mktempdir() do dir
                stdout = IOBuffer()
                config = SandboxConfig(
                    Dict("/" => rootfs_dir),
                    Dict("/glados" => dir);
                )
                @test run(executor(), config, `/bin/sh -c "echo aperture > /glados/science.txt"`)
                @test isfile(joinpath(dir, "science.txt"))
                @test String(read(joinpath(dir, "science.txt"))) == "aperture\n"
            end
        end

        @testset "pipelining" begin
            pipe = PipeBuffer()
            stdout = IOBuffer()
            first_config = SandboxConfig(
                Dict("/" => rootfs_dir),
                stdout = pipe,
            )
            second_config = SandboxConfig(
                Dict("/" => rootfs_dir),
                stdin = pipe,
                stdout = stdout,
            )
            @test run(executor(), first_config, `/bin/sh -c "echo 'ignore me'; echo 'pick this up foo'; echo 'ignore me as well'"`)
            @test run(executor(), second_config, `/bin/sh -c "grep foo"`)
            @test String(take!(stdout)) == "pick this up foo\n";
        end

        @testset "read-only mounts are really read-only" begin
            mktempdir() do dir
                read_only_dir = joinpath(dir, "read_only")
                read_write_dir = joinpath(dir, "read_write")
                mkdir(read_only_dir)
                mkdir(read_write_dir)
                stdout = IOBuffer()
                stderr = IOBuffer()
                config = SandboxConfig(
                    Dict("/" => rootfs_dir, "/read_only" => read_only_dir),
                    Dict("/read_write" => read_write_dir),
                    stdout = stdout,
                    stderr = stderr,
                    persist = false,
                )
                # Modifying the rootfs works, and is temporary; for docker containers this is modifying
                # the rootfs image, for userns this is all mounted within an overlay backed by a tmpfs,
                # because we have `persist` set to `false`.
                with_executor(executor) do exe
                    @test run(exe, config, `/bin/sh -c "echo aperture >> /bin/science && cat /bin/science"`)
                    @test String(take!(stdout)) == "aperture\n";
                    @test isempty(take!(stderr))
                    @test run(exe, config, `/bin/sh -c "echo aperture >> /bin/science && cat /bin/science"`)
                    @test String(take!(stdout)) == "aperture\n";
                    @test isempty(take!(stderr))

                    # An actual read-only mount will not allow writing, because it's truly read-only
                    @test !run(exe, config, ignorestatus(`/bin/sh -c "echo aperture >> /read_only/science && cat /read_only/science"`))
                    @test occursin("Read-only file system", String(take!(stderr)))

                    # A read-write mount, on the other hand, will be permanent
                    @test run(exe, config, `/bin/sh -c "echo aperture >> /read_write/science && cat /read_write/science"`)
                    @test String(take!(stdout)) == "aperture\n";
                    @test isempty(take!(stderr))
                    @test run(exe, config, `/bin/sh -c "echo aperture >> /read_write/science && cat /read_write/science"`)
                    @test String(take!(stdout)) == "aperture\naperture\n";
                    @test isempty(take!(stderr))
                end
            end
        end

        @testset "entrypoint" begin
            mktempdir() do dir
                read_only_dir = joinpath(dir, "read_only")
                mkdir(read_only_dir)
                stdout = IOBuffer()
                stderr = IOBuffer()
                config = SandboxConfig(
                    Dict("/" => rootfs_dir, "/read_only" => read_only_dir),
                    entrypoint = "/read_only/entrypoint",
                    stdout = stdout,
                    stderr = stderr,
                )

                # Generate an `entrypoint` script that mounts a tmpfs-backed overlayfs over our read-only mounts
                # Allowing us to write to those read-only mounts, but the changes are temporary
                open(joinpath(read_only_dir, "entrypoint"), write=true) do io
                    write(io, """
                    #!/bin/sh

                    echo entrypoint activated

                    mkdir /overlay_workdir
                    mount -t tmpfs -osize=1G tmpfs /overlay_workdir
                    mkdir -p /overlay_workdir/upper
                    mkdir -p /overlay_workdir/work
                    mount -t overlay overlay -olowerdir=/read_only -oupperdir=/overlay_workdir/upper -oworkdir=/overlay_workdir/work /read_only

                    exec "\$@"
                    """)
                end
                chmod(joinpath(read_only_dir, "entrypoint"), 0o755)

                # Modifying the read-only files now works, and is temporary
                with_executor(executor) do exe
                    @test run(exe, config, `/bin/sh -c "echo aperture >> /read_only/science && cat /read_only/science"`)
                    @test String(take!(stdout)) == "entrypoint activated\naperture\n";
                    @test isempty(take!(stderr))
                    @test run(exe, config, `/bin/sh -c "echo aperture >> /read_only/science && cat /read_only/science"`)
                    @test String(take!(stdout)) == "entrypoint activated\naperture\n";
                    @test isempty(take!(stderr))
                end
            end
        end

        @testset "persistence" begin
            mktempdir() do dir
                stdout = IOBuffer()
                stderr = IOBuffer()
                config = SandboxConfig(
                    Dict("/" => rootfs_dir),
                    stdout = stdout,
                    stderr = stderr,
                    persist = true,
                )

                # Modifying the read-only files is persistent within a single executor
                cmd = `/bin/sh -c "echo aperture >> /bin/science && cat /bin/science"`
                with_executor(executor) do exe
                    @test run(exe, config, cmd)
                    @test String(take!(stdout)) == "aperture\n";
                    @test isempty(take!(stderr))
                    @test run(exe, config, cmd)
                    @test String(take!(stdout)) == "aperture\naperture\n";
                    @test isempty(take!(stderr))
                end

                with_executor(executor) do exe
                    @test run(exe, config, cmd)
                    @test String(take!(stdout)) == "aperture\n";
                    @test isempty(take!(stderr))
                end
            end
        end
    end
end

@testset "default executor" begin
    stdout = IOBuffer()
    stderr = IOBuffer()
    config = SandboxConfig(
        Dict("/" => rootfs_dir);
        stdout,
        stderr,
    )
    with_executor() do exe
        @test run(exe, config, `/bin/sh -c "echo stdout; echo stderr >&2"`)
        @test String(take!(stdout)) == "stdout\n";
        @test String(take!(stderr)) == "stderr\n";
    end
end


text:
"""
    cKeplerFlow!(x_bek::Array{Float64,1}, v_bek::Array{Float64,1}, x0_bek::Array{Float64,1}, v0_bek::Array{Float64,1}, dt::Any, mu::Float64, traza::Int32)

Defines the position and speeds of a particle -after a time-step- orbiting another body given its initial position and speeds. It does the same as `keplerFlow` function, but it is faster (because it is written in C).

# Args

* `x_bek`: the position vector of the body after dt time from the initial values will be written in this array of D elements, where D is the number of used dimensions.
* `v_bek`: the velocity vector of the body after dt time from the initial values will be written in this array of D elements.
* `x0_bek`: position vector of the secondary body orbiting the primary body.
* `v0_bek`: speed vector of the secondary body orbiting the primary body.
* `dt`: time-step, the function returns the new position and speed vectors after dt time has passed.
* `mu`: standard gravitational parameter of the two bodies.
* `r0`: initial distance between the two bodies.
* `Gz`: array of the first four g-functions. Instead of defining it on each execution of the function, the same array is used in the whole simulation. That is why it's passed as an argument.
* `trace`: boolean value that shows information of the execution if true.
"""
cKeplerFlow!(x_bek::Array{Float64,1}, v_bek::Array{Float64,1}, x0_bek::Array{Float64,1}, v0_bek::Array{Float64,1}, dt::Any, mu::Float64, traza::Int32) = open(string(Base.source_path()[1:(end-8)],"kepler.c")) do f

    code = read(f, String)

    # liburutegiaren PATH-a hemen gordeko da
    Clib = mktemp()[1]

    # konpilatzeko, C-ko math liburutegia behar da
    # liburutegia konpilatuko dugu, -lm
    open(`gcc -fPIC -O3 -xc -shared -lm -o $(Clib * "." * Libdl.dlext) -`, "w") do f
         print(f, code)
    end

    # definitu Julia funtzio bat, C funtzioa deitzen duena
    cKeplerFlow!(x_bek::Array{Float64,1}, v_bek::Array{Float64,1}, x0_bek::Array{Float64,1}, v0_bek::Array{Float64,1}, dt::Any, mu::Float64, traza::Int32) =
        ccall(("kepler", Clib), Nothing,(Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},Float64,Float64,Int32),
        x_bek, v_bek, x0_bek, v0_bek, convert(Float64,dt), mu, traza)
end

max_stars_repo_name: salanueva/KeplerFlow

max_issues_repo_name: salanueva/KeplerFlow

max_forks_repo_name: salanueva/KeplerFlow


text:
@testset "560.subarray-sum-equals-k.jl" begin
    @test subarray_sum([1, 1, 1], 2) == 2
    @test subarray_sum([8, 17, -26, 7, 19, -7, 0, -15, -25, -30, 22, 9, 9, 11, 6, 5, 23, -28, -8, 5, 8, 11, -11, -20, -10, -23, 17, -24, 0, -30, -8, -1, -20, -1, 2, 21, -6, -2, 22, -24, 15, 7, -26, 23, 19, 8, -20, -8, 28, 16, -24, -26, 9, -14, -2, -25, -15, -27, -14, -26, -6, -27, 27, 3, 16, 8, -1, 17, 12, 24, -15, 22, -7, -27, -13, 20, -26, -30, 24, 15, 27, 30, -28, -5, 24, 12, -24, 17, 9, -7], 3) == 17
end


text:
# Interacting with the Petsc Options Database
export PetscOptionsSetValue, PetscOptionsClearValue, PetscOptionsView, PetscSetOptions, PetscClearOptions

typealias PetscViewer Ptr{Void}

function PetscOptionsSetValue(arg1::AbstractString,arg2::AbstractString)
    ccall((:PetscOptionsSetValue,petsc),PetscErrorCode,(Cstring, Cstring),arg1,arg2)
end

function PetscOptionsView(arg1::PetscViewer=C_NULL)
    ccall((:PetscOptionsView,petsc),PetscErrorCode,(PetscViewer,),arg1)
end

function PetscOptionsClearValue(arg1::AbstractString)
    ccall((:PetscOptionsClearValue,petsc),PetscErrorCode,(Cstring,),arg1)
end

"""
  Convenience wrapper for using a dictionary to set options
"""
function PetscSetOptions(opts::Dict)

  for (key, value) in opts
    PetscOptionsSetValue(key, value)
  end

end

function PetscClearOptions(opts::Dict)

  for key in keys(opts)
    PetscOptionsClearValue(key)
  end

end


text:
function RLCore.Experiment(
    ::Val{:JuliaRL},
    ::Val{:A2CGAE},
    ::Val{:CartPole},
    ::Nothing;
    save_dir = nothing,
    seed = 123,
)
    if isnothing(save_dir)
        t = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
        save_dir = joinpath(pwd(), "checkpoints", "JuliaRL_A2CGAE_CartPole_$(t)")
    end

    lg = TBLogger(joinpath(save_dir, "tb_log"), min_level = Logging.Info)
    rng = StableRNG(seed)
    N_ENV = 16
    UPDATE_FREQ = 10
    env = MultiThreadEnv([
        CartPoleEnv(; T = Float32, rng = StableRNG(hash(seed + i))) for i in 1:N_ENV
    ])
    ns, na = length(get_state(env[1])), length(get_actions(env[1]))
    RLBase.reset!(env, is_force = true)
    agent = Agent(
        policy = QBasedPolicy(
            learner = A2CGAELearner(
                approximator = ActorCritic(
                    actor = NeuralNetworkApproximator(
                        model = Chain(
                            Dense(ns, 256, relu; initW = glorot_uniform(rng)),
                            Dense(256, na; initW = glorot_uniform(rng)),
                        ),
                        optimizer = ADAM(1e-3),
                    ),
                    critic = NeuralNetworkApproximator(
                        model = Chain(
                            Dense(ns, 256, relu; initW = glorot_uniform(rng)),
                            Dense(256, 1; initW = glorot_uniform(rng)),
                        ),
                        optimizer = ADAM(1e-3),
                    ),
                ) |> cpu,
                γ = 0.99f0,
                λ = 0.97f0,
                actor_loss_weight = 1.0f0,
                critic_loss_weight = 0.5f0,
                entropy_loss_weight = 0.001f0,
                update_freq = UPDATE_FREQ
            ),
            explorer = BatchExplorer(GumbelSoftmaxExplorer(;)),
        ),
        trajectory = CircularArraySARTTrajectory(;
            capacity = UPDATE_FREQ,
            state = Matrix{Float32} => (ns, N_ENV),
            action = Vector{Int} => (N_ENV,),
            reward = Vector{Float32} => (N_ENV,),
            terminal = Vector{Bool} => (N_ENV,),
        ),
    )
    stop_condition = StopAfterStep(50_000)
    total_reward_per_episode = TotalBatchRewardPerEpisode(N_ENV)
    time_per_step = TimePerStep()
    hook = ComposedHook(
        total_reward_per_episode,
        time_per_step,
        DoEveryNStep() do t, agent, env
            with_logger(lg) do
                @info(
                    "training",
                    actor_loss = agent.policy.learner.actor_loss,
                    critic_loss = agent.policy.learner.critic_loss,
                    entropy_loss = agent.policy.learner.entropy_loss,
                    loss = agent.policy.learner.loss,
                )
                for i in 1:length(env)
                    if get_terminal(env[i])
                        @info "training" reward = total_reward_per_episode.rewards[i][end] log_step_increment =
                            0
                        break
                    end
                end
            end
        end,
    )
    Experiment(
        agent,
        env,
        stop_condition,
        hook,
        "# A2CGAE with CartPole",
    )
end


text:
function makeenv(;H = 0.8, dims=(100,100))
	sdm = rand(MidpointDisplacement(H), dims...)
end

function makesdm(;H = 0.8, dims=(100,100), cutoff= 0.75)
	sdm = rand(MidpointDisplacement(H), dims...)
	sdm = sdm ./ sum(sdm)
	ind = findall(x->x==1, sdm .<= quantile(vec(sdm), cutoff))
	sdm[ind] .= 0
	sdm = sdm ./ sum(sdm)
	sdm
end

function makeoccurrence(sdm; npts=100)
    occmat = zeros(size(sdm))
    for _ in 1:npts
        ind = rand(Categorical(vec(sdm)))
        occmat[ind] += 1
    end
    occmat
    return findall(i->i>0, occmat)
end

function makebon(sdm, n=50, padding=0.15)
	x,y = size(sdm)
	pts = rand((n,2)...)

	pts[:,1] = floor.((pts[:,1] .* (1-padding)*x) .+ (0.5padding*x))
	pts[:,2] = floor.((pts[:,2] .* (1-padding)*y) .+ (0.5padding*y))
	convert.(Int32,pts)
end


text:
module Sequent
using MLStyle
export @sequent, @semantics
@nospecialize
@data Judgement begin
    Assertion(
        line::LineNumberNode,
        rule_name::Symbol,
        vectorized::Bool,
        judge_op::Symbol,
        context,
        pattern,
        ret)

    EqMatch(line::LineNumberNode, pattern, value)
end

struct SequentDef
    conditions :: Vector{Judgement}
    case :: Assertion
end

macro sequent(name, formula)
    name isa Symbol || error("formula name must be a symbol!")
    @match formula Expr(:block, stmts...) =>
    begin
        conditions = Judgement[]
        line = __source__
        for stmt in stmts
            @switch stmt begin
            @case line::LineNumberNode
                continue
            @case :($a = $b)
                push!(conditions, EqMatch(line, a, b))
                continue
            @case Expr(:call,
                    :(=>) && let is_vec = false end ||
                    :(.=>) && let is_vec = true end,
                    Expr(:call, judge_op::Symbol, context, pattern),
                    result) &&
                    if startswith(string(judge_op), "⊢") end
                push!(conditions,
                      Assertion(
                        line, name, is_vec,
                        judge_op, context, pattern, result))
                continue
            end
        end
        length(conditions) >= 1 && conditions[end] isa Assertion ||
        error("malformed sequent definition $name at $__source__")

        case = pop!(conditions)
        SequentDef(conditions, case)
    end
end

function semantics(seqs::Vector{SequentDef}, __source__::LineNumberNode)
    groups = Dict{Symbol, Vector{SequentDef}}()
    for seq in seqs
        group = get!(groups, seq.case.judge_op) do
            SequentDef[]
        end
        push!(group, seq)
    end
    ctx_sym = gensym("ctx")
    term_sym = gensym("term")
    code = Expr(:block)
    for (funcname, group) in groups
        when_cases = Tuple{Expr, Any}[]
        for seq in group
            case_clauses = Expr(:block)
            main_case = Expr(:tuple, seq.case.context, seq.case.pattern)
            push!(
                case_clauses.args,
                :($main_case = ($ctx_sym, $term_sym)))
            for cond in seq.conditions
                @match cond begin
                    EqMatch(line, pattern, value) =>
                        begin
                            value = Expr(:block, line, value)
                            push!(case_clauses.args, :($pattern = $value))
                        end
                    Assertion(
                        line, name, is_vec,
                        judge_op, context, pattern, result) =>
                        begin
                            context === :_ && (context = :nothing)
                            fcall = is_vec ?
                                :(($term_sym -> $judge_op($context, $term_sym)).($pattern)) :
                                :($judge_op($context, $pattern))
                            fcall = Expr(:block, line, fcall)
                            push!(case_clauses.args, :($result = $fcall))
                        end
                end
            end
            push!(
                when_cases,
                (case_clauses, seq.case.ret)
            )
        end
        hd, tl = when_cases[1], @view when_cases[2:end]
        when_block = Expr(:block)
        push!(when_block.args, hd[2])
        for (binding, result) in tl
            push!(when_block.args, :(@when $binding))
            push!(when_block.args, result)
        end
        push!(when_block.args, :(@otherwise))
        push!(when_block.args, :($error("unrecognised input" * $string($term_sym))))
        match_expression =
            Expr(:let,
                hd[1],
                when_block)

        push!(
            code.args,
            Expr(
                :function,
                Expr(:call, funcname, ctx_sym, term_sym),
                Expr(:block,
                    Expr(
                        :macrocall,
                        GlobalRef(MLStyle, Symbol("@when")),
                        __source__,
                        match_expression))))
    end
    code
end

macro semantics(seqs)
    esc(semantics(__module__.eval(seqs), __source__))
end
@specialize
end


text:
mutable struct D4COption
    threshold::Float64

    # allow uninitialized creation
    D4COption() = new()
end

function d4c(x::Vector{Float64}, fs::Int, temporal_positions::Vector{Float64}, f0::Vector{Float64}, fft_size::Int)

    d4c_option = Ref(D4COption())
    @ccall libworld.InitializeD4COption(d4c_option::Ptr{D4COption})::Cvoid

    nbins = length(temporal_positions)
    nfreqs = (fft_size ÷ 2) + 1
    aperiodicity = Matrix{Float64}(undef, nfreqs, nbins) # transposed for C

    aperiodicity_pointer_vector = [Ref(aperiodicity, offset)
        for offset in 1:nfreqs:length(aperiodicity)]

    @ccall libworld.D4C(
        x::Ptr{Float64},
        length(x)::Cint,
        fs::Cint,
        temporal_positions::Ptr{Float64},
        f0::Ptr{Float64},
        length(f0)::Cint,
        fft_size::Cint,
        d4c_option::Ptr{D4COption},
        aperiodicity_pointer_vector::Ptr{Ptr{Cdouble}})::Cvoid

    aperiodicity
end


text:
abstract type AbstractFlow end

# The vector field `F(x, p)` associated to a Cauchy problem ẋ = F(x, p). Used for the differential of the shooting problem. The vector field is used like `vf(flow, x, p)`` and must return `F(x, p)``
function vf(::AbstractFlow, x, par; k...) end

# these functions are used in the Standard Shooting method
# the function implements the flow (or semigroup) `(x, p, t) -> flow(x, p, t)` associated to an autonomous Cauchy problem. Only the last time point must be returned in the form Named Tuple `(u = ..., t = t)`. In the case of Poincaré Shooting, one must be able to call the flow like `evolve(fl, x, par, Inf)`.
function evolve(::AbstractFlow, x, par, δt; k...) end

# The differential `dflow` of the flow *w.r.t.* `x`, `(x, p, dx, t) -> dflow(x, p, dx, t)`. One important thing is that we require `dflow(x, dx, t)` to return a Named Tuple: `(t = t, u = flow(x, p, t), du = dflow(x, p, dx, t))`, the last component being the value of the derivative of the flow.
function evolve(::AbstractFlow, x, par, dx, δt; k...) end

# [Optional] The function implements the flow (or semigroup) associated to an autonomous Cauchy problem `(x, p, t) -> flow(x, p, t)`. The whole solution on the time interval [0,t] must be returned. It is not strictly necessary to provide this, it is mainly used for plotting on the user side. In the case of Poincaré Shooting, one must be able to call the flow like `evolve(fl, Val(:Full), x, par, Inf)`.
function evolve(::AbstractFlow, ::Val{:Full}, x, par, δt; k...) end

# [Optional / Internal] Serial version of the flow. Used for Matrix based jacobian (Shooting and Poincaré Shooting) and diffPoincareMap. Must return a Named Tuple `(u = ..., t = t)`
function evolve(fl::AbstractFlow, ::Val{:SerialTimeSol}, x, par, δt; k...) end

# [Optional] Flow which returns the tuple `(t, u(t))`. Optional, mainly used for plotting on the user side.
function evolve(::AbstractFlow, ::Val{:TimeSol}, x, par, δt = Inf; k...) end

# [Optional] Serial version of `dflow`. Used internally when using parallel multiple shooting. Named Tuple `(u = ..., du = ..., t = t)`
function evolve(::AbstractFlow, ::Val{:SerialdFlow}, x, par, dx, tΣ; kwargs...) end

####################################################################################################
# Structures related to computing ODE/PDE Flows
"""
$(TYPEDEF)
$(TYPEDFIELDS)

# Simplified constructor(s)
We provide a simple constructor where you only pass the vector field `F`, the flow `ϕ` and its differential `dϕ`:

	fl = Flow(F, ϕ, dϕ)

# Simplified constructors for DifferentialEquations.jl

These are some simple constructors for which you only have to pass a `prob::ODEProblem` or `prob::EnsembleProblem` (for parallel computation) from `DifferentialEquations.jl` and an ODE time stepper like `Tsit5()`. Hence, you can do for example

	fl = Flow(prob, Tsit5(); kwargs...)

where `kwargs` is passed to `SciMLBase::solve`. If your vector field depends on parameters `p`, you can define a `Flow` using

	fl = Flow(prob, Tsit5(); kwargs...)

Finally, you can pass two `ODEProblem` where the second one is used to compute the variational equation:

	fl = Flow(prob1::ODEProblem, alg1, prob2::ODEProblem, alg2; kwargs...)

"""
@with_kw struct Flow{TF, Tf, Tts, Tff, Td, Tse, Tprob, TprobMono, Tfs, Tcb} <: AbstractFlow
	"The vector field `(x, p) -> F(x, p)` associated to a Cauchy problem. Used for the differential of the shooting problem."
	F::TF = nothing

	"The flow (or semigroup) `(x, p, t) -> flow(x, p, t)` associated to the Cauchy problem. Only the last time point must be returned in the form (u = ...)"
	flow::Tf = nothing

	"Flow which returns the tuple (t, u(t)). Optional, mainly used for plotting on the user side."
	flowTimeSol::Tts = nothing

	"[Optional] The flow (or semigroup) associated to the Cauchy problem `(x, p, t) -> flow(x, p, t)`. The whole solution on the time interval [0,t] must be returned. It is not strictly necessary to provide this, it is mainly used for plotting on the user side. Please use `nothing` as default."
	flowFull::Tff = nothing

	"The differential `dflow` of the flow *w.r.t.* `x`, `(x, p, dx, t) -> dflow(x, p, dx, t)`. One important thing is that we require `dflow(x, dx, t)` to return a Named Tuple: `(t = t, u = flow(x, p, t), du = dflow(x, p, dx, t))`, the last component being the value of the derivative of the flow."
	dflow::Td = nothing

	"[Optional] Serial version of dflow. Used internally when using parallel multiple shooting. Please use `nothing` as default."
	dfSerial::Tse = nothing

	"[Internal] store the ODEProblem associated to the flow of the Cauchy problem"
	prob::Tprob = nothing

	"[Internal] store the ODEProblem associated to the flow of the variational problem"
	probMono::TprobMono = nothing

	"[Internal] Serial version of the flow"
	flowSerial::Tfs = nothing

	"[Internal] Store possible callback"
	callback::Tcb = nothing
end

# constructors
Flow(F, fl, df = nothing) = Flow(F = F, flow = fl, dflow = df, dfSerial = df)

vf(fl::Flow, x, p) = fl.F(x, p)

evolve(fl::Flow, x, p, t; k...)                        = fl.flow(x, p, t; k...)
evolve(fl::Flow, x, p, dx, t; k...)                    = fl.dflow(x, p, dx, t; k...)
evolve(fl::Flow, ::Val{:Full}, x, p, t; k...)          = fl.flowFull(x, p, t; k...)

# for Poincaré Shooting
evolve(fl::Flow, ::Val{:SerialTimeSol}, x, p, t; k...)   = fl.flowSerial(x, p, t; k...)
evolve(fl::Flow, ::Val{:SerialdFlow}, x, p, dx, t; k...) = fl.dfSerial(x, p, dx, t; k...)


text:
entries = parse.(Int, readlines("input.txt"))

for (i, x) in enumerate(entries)
    for j in 1:i-1
        y = entries[j]
        if x + y == 2020
            @info "Part I" x y x*y
        end
        for k in 1:j-1
            z = entries[k]
            if x + y + z == 2020
                @info "Part II" x y z x*y*z
            end
        end
    end
end

max_stars_repo_name: dglmoore/advent-of-code

max_issues_repo_name: dglmoore/advent-of-code

max_forks_repo_name: dglmoore/advent-of-code


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.test.behaviour.concept.thing.relation;
#
# import typedb.client.api.concept.thing.Attribute;
# import typedb.client.api.concept.thing.Relation;
# import typedb.client.api.concept.thing.Thing;
# import typedb.client.common.Label;
# import io.cucumber.java.en.Then;
# import io.cucumber.java.en.When;
#
# import java.time.LocalDateTime;
# import java.util.List;
# import java.util.Map;
#
# import static typedb.client.test.behaviour.concept.thing.ThingSteps.get;
# import static typedb.client.test.behaviour.concept.thing.ThingSteps.put;
# import static typedb.client.test.behaviour.connection.ConnectionStepsBase.tx;
# import static typedb.client.test.behaviour.util.Util.assertThrows;
# import static org.junit.Assert.assertEquals;
# import static org.junit.Assert.assertFalse;
# import static org.junit.Assert.assertTrue;
#
# @SuppressWarnings("CheckReturnValue")
# public class RelationSteps {
#
#     @When("{var} = relation\\( ?{type_label} ?) create new instance")
#     public void relation_type_create_new_instance(String var, String typeLabel) {
#         put(var, tx().concepts().getRelationType(typeLabel).asRemote(tx()).create());
#     }
#
#     @Then("relation\\( ?{type_label} ?) create new instance; throws exception")
#     public void relation_type_create_new_instance_throws_exception(String typeLabel) {
#         assertThrows(() -> tx().concepts().getRelationType(typeLabel).asRemote(tx()).create());
#     }
#
#     @When("{var} = relation\\( ?{type_label} ?) create new instance with key\\( ?{type_label} ?): {int}")
#     public void relation_type_create_new_instance_with_key(String var, String type, String keyType, int keyValue) {
#         Attribute.Long key = tx().concepts().getAttributeType(keyType).asLong().asRemote(tx()).put(keyValue);
#         Relation relation = tx().concepts().getRelationType(type).asRemote(tx()).create();
#         relation.asRemote(tx()).setHas(key);
#         put(var, relation);
#     }
#
#     @When("{var} = relation\\( ?{type_label} ?) create new instance with key\\( ?{type_label} ?): {word}")
#     public void relation_type_create_new_instance_with_key(String var, String type, String keyType, String keyValue) {
#         Attribute.String key = tx().concepts().getAttributeType(keyType).asString().asRemote(tx()).put(keyValue);
#         Relation relation = tx().concepts().getRelationType(type).asRemote(tx()).create();
#         relation.asRemote(tx()).setHas(key);
#         put(var, relation);
#     }
#
#     @When("{var} = relation\\( ?{type_label} ?) create new instance with key\\( ?{type_label} ?): {datetime}")
#     public void relation_type_create_new_instance_with_key(String var, String type, String keyType, LocalDateTime keyValue) {
#         Attribute.DateTime key = tx().concepts().getAttributeType(keyType).asDateTime().asRemote(tx()).put(keyValue);
#         Relation relation = tx().concepts().getRelationType(type).asRemote(tx()).create();
#         relation.asRemote(tx()).setHas(key);
#         put(var, relation);
#     }
#
#     @When("{var} = relation\\( ?{type_label} ?) get instance with key\\( ?{type_label} ?): {long}")
#     public void relation_type_get_instance_with_key(String var1, String type, String keyType, long keyValue) {
#         put(var1, tx().concepts().getAttributeType(keyType).asLong().asRemote(tx()).get(keyValue).asRemote(tx()).getOwners()
#                 .filter(owner -> owner.getType().getLabel().equals(Label.of(type))).findFirst().orElse(null));
#     }
#
#     @When("{var} = relation\\( ?{type_label} ?) get instance with key\\( ?{type_label} ?): {word}")
#     public void relation_type_get_instance_with_key(String var1, String type, String keyType, String keyValue) {
#         put(var1, tx().concepts().getAttributeType(keyType).asString().asRemote(tx()).get(keyValue).asRemote(tx()).getOwners()
#                 .filter(owner -> owner.getType().getLabel().equals(Label.of(type))).findFirst().orElse(null));
#     }
#
#     @When("{var} = relation\\( ?{type_label} ?) get instance with key\\( ?{type_label} ?): {datetime}")
#     public void relation_type_get_instance_with_key(String var1, String type, String keyType, LocalDateTime keyValue) {
#         put(var1, tx().concepts().getAttributeType(keyType).asDateTime().asRemote(tx()).get(keyValue).asRemote(tx()).getOwners()
#                 .filter(owner -> owner.getType().getLabel().equals(Label.of(type))).findFirst().orElse(null));
#     }
#
#     @Then("relation\\( ?{type_label} ?) get instances contain: {var}")
#     public void relation_type_get_instances_contain(String typeLabel, String var) {
#         assertTrue(tx().concepts().getRelationType(typeLabel).asRemote(tx()).getInstances().anyMatch(i -> i.equals(get(var))));
#     }
#
#     @Then("relation\\( ?{type_label} ?) get instances do not contain: {var}")
#     public void relation_type_get_instances_do_not_contain(String typeLabel, String var) {
#         assertTrue(tx().concepts().getRelationType(typeLabel).asRemote(tx()).getInstances().noneMatch(i -> i.equals(get(var))));
#     }
#
#     @Then("relation\\( ?{type_label} ?) get instances is empty")
#     public void relation_type_get_instances_is_empty(String typeLabel) {
#         assertEquals(0, tx().concepts().getRelationType(typeLabel).asRemote(tx()).getInstances().count());
#     }
#
#     @When("relation {var} add player for role\\( ?{type_label} ?): {var}")
#     public void relation_add_player_for_role(String var1, String roleTypeLabel, String var2) {
#         get(var1).asRelation().asRemote(tx()).addPlayer(get(var1).asRelation().getType().asRemote(tx()).getRelates(roleTypeLabel), get(var2));
#     }
#
#     @When("relation {var} add player for role\\( ?{type_label} ?): {var}; throws exception")
#     public void relation_add_player_for_role_throws_exception(String var1, String roleTypeLabel, String var2) {
#         assertThrows(() -> get(var1).asRelation().asRemote(tx()).addPlayer(get(var1).asRelation().getType().asRemote(tx()).getRelates(roleTypeLabel), get(var2)));
#     }
#
#     @When("relation {var} remove player for role\\( ?{type_label} ?): {var}")
#     public void relation_remove_player_for_role(String var1, String roleTypeLabel, String var2) {
#         get(var1).asRelation().asRemote(tx()).removePlayer(get(var1).asRelation().getType().asRemote(tx()).getRelates(roleTypeLabel), get(var2));
#     }
#
#     @Then("relation {var} get players contain:")
#     public void relation_get_players_contain(String var, Map<String, String> players) {
#         Relation relation = get(var).asRelation();
#         players.forEach((rt, var2) -> assertTrue(relation.asRemote(tx()).getPlayersByRoleType().get(relation.getType().asRemote(tx()).getRelates(rt)).contains(get(var2.substring(1)))));
#     }
#
#     @Then("relation {var} get players do not contain:")
#     public void relation_get_players_do_not_contain(String var, Map<String, String> players) {
#         Relation relation = get(var).asRelation();
#         players.forEach((rt, var2) -> {
#             List<? extends Thing> p;
#             if ((p = relation.asRemote(tx()).getPlayersByRoleType().get(relation.getType().asRemote(tx()).getRelates(rt))) != null) {
#                 assertFalse(p.contains(get(var2.substring(1))));
#             }
#         });
#     }
#
#     @Then("relation {var} get players contain: {var}")
#     public void relation_get_players_contain(String var1, String var2) {
#         assertTrue(get(var1).asRelation().asRemote(tx()).getPlayers().anyMatch(p -> p.equals(get(var2))));
#     }
#
#     @Then("relation {var} get players do not contain: {var}")
#     public void relation_get_players_do_not_contain(String var1, String var2) {
#         assertTrue(get(var1).asRelation().asRemote(tx()).getPlayers().noneMatch(p -> p.equals(get(var2))));
#     }
#
#     @Then("relation {var} get players for role\\( ?{type_label} ?) contain: {var}")
#     public void relation_get_player_for_role_contain(String var1, String roleTypeLabel, String var2) {
#         assertTrue(get(var1).asRelation().asRemote(tx())
#                            .getPlayers(get(var1).asRelation().getType().asRemote(tx()).getRelates(roleTypeLabel))
#                            .anyMatch(p -> p.equals(get(var2))));
#     }
#
#     @Then("relation {var} get players for role\\( ?{type_label} ?) do not contain: {var}")
#     public void relation_get_player_for_role_do_not_contain(String var1, String roleTypeLabel, String var2) {
#         assertTrue(get(var1).asRelation().asRemote(tx())
#                            .getPlayers(get(var1).asRelation().getType().asRemote(tx()).getRelates(roleTypeLabel))
#                            .noneMatch(p -> p.equals(get(var2))));
#     }
# }


text:
"""
    combinationrank(xs...)

Returns the index of the combination `xs` in the combinatorial number system. The elements
of `xs` all have to be unique, greater than zero and sorted in ascending order.
"""
function combinationrank(xs...)
    @assert length(unique(xs)) == length(xs)
    @assert sort([xs...]) == [xs...]
    @assert all(xs .> 0)
    k = length(xs)
    1 + sum(binomial.(xs .- 1, 1:k))
end

"""
    combinationunrank(k, i)

Returns the `i`-th `k`-element combination in the combinatorial number system.
"""
function combinationunrank(k, idx)
    @assert k > 0
    @assert idx > 0
    v = zeros(Int, k)
    for i = k:-1:1
        ck = _max_ck(i, idx - 1)
        v[i] = ck + 1
        idx -= binomial(ck, i)
    end
    tuple(v...)
end

function _max_ck(k, idx)
    ck = 0
    while true
        binomial(ck + 1, k) > idx && return ck
        ck += 1
    end
    error("Unreachable reached.")
end

"""
    findcombinationdiff(c1::NTuple{N,Int}, c2::NTuple{N,Int}) -> (idx1, idx2) or nothing

If the two ordered combinations of the same length have all elements except one in common,
return a pair of indices of those elements. Otherwise, return `nothing`.
"""
function findcombinationdiff(c1::NTuple{N,Int}, c2::NTuple{N,Int}) where N
    idx1::Union{Nothing,Int}, idx2::Union{Nothing,Int} = nothing, nothing
    i = 1
    while i < N
        i1 = isnothing(idx1) ? i : i + 1
        i2 = isnothing(idx2) ? i : i + 1
        if c1[i1] == c2[i2]
            i += 1
        elseif (c1[i1] < c2[i2]) && isnothing(idx1)
            idx1 = i1
        elseif (c1[i1] > c2[i2]) && isnothing(idx2)
            idx2 = i2
        else
            return nothing
        end
    end
    isnothing(idx1) && (idx1 = N)
    isnothing(idx2) && (idx2 = N)
    c1[idx1] == c2[idx2] ? nothing : (idx1, idx2)
end


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

"""
!STATEV(1)=f:Porosity
!STATEV(2)=ep:equivalent plastic strain
!STATEV(3)=DEp: Hydrostatic inelastic strain increment
!STATEV(4)=DEp: Deviatoric inelastic strain increment
!STATEV(5)...STATEV(10)=EELAS
!STATEV(11)...STATEV(16)=EPLAS
!STATEV(17)...STATEV(22)=X

!Properties for porosity
f0=PROPS(1)!f0,initial porosity[GJS-700:0.02,GJV-450:0.05,GJL-250:0.125]
fc=PROPS(2)!fc,critical porosity[0.0325,0.08,1.0]
ff=PROPS(3)!ff,failure porosity[0.25,0.1,1.0]
fn=PROPS(4)!fn,volume fraction of graphite inclusions[0.04,0.1,0.25]
sn=PROPS(5)!sn,standard deviation of graphite inclusions[0.001,0.001,-]
en=PROPS(6)!en,mean value of graphite inclusions[0.0,0.0,-]
qinf=PROPS(7)!Isotropic hardening amplitude
b=PROPS(8)!Isotropic hardening exponent sflow=syield+qinf*(1-EXP(-b*ep))
C=PROPS(9)!Cyclic hardening parameter
ga=PROPS(10)!Cyclic saturation parameter. In von Mises plasticity C/ga=qinf,ga=b
EE=PROPS(11)!Young's modulus
nu=PROPS(12)!Poisson's ratio
syield=PROPS(13)!Yield stress
"""
function GursonMaterial(initial_porosity=0.05, critical_porosity=0.08, failure_porosity=0.1,
                        volume_fraction_of_graphite=0.1, sdev_of_graphite=0.001,
                        mean_of_graphite=0.0, isotropic_hardening_amplitude=100.0,
                        isotropic_hardening_exponent=1.0, cyclic_hardening_parameter=1.0,
                        cyclic_saturation_parameter=1.0, Youngs_modulus=169000.0,
                        Poissons_ratio=0.3, yield_stress=150.0)

    umat_params= [initial_porosity,critical_porosity,failure_porosity,volume_fraction_of_graphite,
                  sdev_of_graphite, mean_of_graphite, isotropic_hardening_amplitude,
                  isotropic_hardening_exponent, cyclic_hardening_parameter, cyclic_saturation_parameter,
                  Youngs_modulus, Poissons_ratio, yield_stress]

    return UmatMaterial(NTENS=6, NSTATV=22, NPROPS=13,
                        parameters=UmatParameterState(NPROPS=13, PROPS=umat_params),
                        lib_path=joinpath(lib_dir,"libgurson_porous_plasticity." * dlext))
end


text:
using QuadraticToBinary, Test, MathOptInterface

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOIT = MOI.DeprecatedTest
const MOIB = MOI.Bridges

optimizers = []
# include("solvers/gurobi.jl")
include("solvers/cbc.jl")

include("moi.jl")

const CONFIG_LOW_TOL = MOIT.Config{Float64}(atol = 1e-3, rtol = 1e-2, duals = false, infeas_certificates = false)
const CONFIG_VERY_LOW_TOL = MOIT.Config{Float64}(atol = 5e-3, rtol = 5e-2, duals = false, infeas_certificates = false)
const CONFIG = MOIT.Config{Float64}(duals = false, infeas_certificates = false)

function test_runtests()
    optimizer = MOI.instantiate(Cbc.Optimizer, with_bridge_type = Float64)
    model = QuadraticToBinary.Optimizer{Float64}(optimizer)
    MOI.set(model, MOI.Silent(), true)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            atol = 1e-3, rtol = 1e-2,
            exclude = Any[
                # MOI.DualStatus,
                # Cbc limitations
                MOI.ConstraintDual,
                MOI.DualObjectiveValue,
                MOI.ConstraintBasisStatus,
                MOI.VariableBasisStatus,
            ]
            ),
        exclude = String[
            # require bounds on variables:
            "test_quadratic_nonhomogeneous",
            "test_quadratic_nonconvex_constraint_integration",
            "test_quadratic_nonconvex_constraint_basic",
            "test_quadratic_integration",
            "test_quadratic_duplicate_terms",
            "test_quadratic_constraint_minimize",
            "test_quadratic_constraint_integration",
            "test_quadratic_constraint_basic",
            "test_quadratic_constraint_LessThan",
            "test_quadratic_constraint_GreaterThan",
            "test_quadratic_SecondOrderCone_basic",
            "test_objective_qp_ObjectiveFunction_zero_ofdiag",
            "test_objective_qp_ObjectiveFunction_edge_cases",
            "test_constraint_qcp_duplicate_off_diagonal",
            "test_constraint_qcp_duplicate_diagonal",
            # TODO: broken in Cbc:
            # https://github.com/jump-dev/Cbc.jl/blob/master/test/MOI_wrapper.jl
            # TODO(odow): bug in Cbc.jl
            "test_model_copy_to_UnsupportedAttribute",
            "test_model_ModelFilter_AbstractConstraintAttribute",
            # TODO(odow): bug in MOI
            "test_model_LowerBoundAlreadySet",
            "test_model_UpperBoundAlreadySet",
            # TODO(odow): upstream bug in Cbc
            "_Indicator_",
            "test_linear_SOS1_integration",
            "test_linear_SOS2_integration",
            "test_solve_SOS2_add_and_delete",
            # Can't prove infeasible.
            "test_conic_NormInfinityCone_INFEASIBLE",
            "test_conic_NormOneCone_INFEASIBLE",
            "test_solve_TerminationStatus_DUAL_INFEASIBLE",
        ],
    )
    #=
        Hard
    =#
    optimizer = MOI.instantiate(Cbc.Optimizer, with_bridge_type = Float64)
    model = QuadraticToBinary.Optimizer{Float64}(optimizer)
    MOI.set(model, QuadraticToBinary.FallbackUpperBound(), +10.0)
    MOI.set(model, QuadraticToBinary.FallbackLowerBound(), -10.0)
    MOI.set(model, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.TimeLimitSec(), 20.0)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            atol = 5e-3,
            rtol = 5e-2,
            exclude = Any[
                MOI.ConstraintDual,
                MOI.DualObjectiveValue
            ],
        ),
        include = String[
            "test_quadratic_nonhomogeneous",
            "test_quadratic_nonconvex_constraint_integration",
            "test_quadratic_nonconvex_constraint_basic",
            "test_quadratic_constraint_minimize",
            "test_quadratic_constraint_integration",
            "test_quadratic_constraint_basic",
            "test_quadratic_constraint_LessThan",
            "test_quadratic_constraint_GreaterThan",
            "test_objective_qp_ObjectiveFunction_zero_ofdiag",
            "test_objective_qp_ObjectiveFunction_edge_cases",
            "test_constraint_qcp_duplicate_off_diagonal",
            "test_constraint_qcp_duplicate_diagonal",
        ]
    )
    #=
        more precise bounds
    =#
    MOI.set(model, QuadraticToBinary.FallbackUpperBound(), 1.0)
    MOI.set(model, QuadraticToBinary.FallbackLowerBound(), 0.0)
    MOI.set(model, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.TimeLimitSec(), 80.0)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            atol = 5e-3,
            rtol = 5e-2,
            optimal_status = MOI.OPTIMAL,
            exclude = Any[
                MOI.ConstraintDual,
                MOI.DualObjectiveValue
            ],
        ),
        include = String[
            "test_quadratic_SecondOrderCone_basic",
            "test_quadratic_integration",
            "test_quadratic_duplicate_terms",
        ]
    )

    return
end
@testset "MOI Unit" begin
    @time test_runtests()
end

@testset "basic_constraint_tests" begin
    for opt in optimizers
        MODEL = QuadraticToBinary.Optimizer{Float64}(opt)
        MOIT.basic_constraint_tests(MODEL, CONFIG)
    end
end

@testset "readme" begin
    for opt in optimizers

        MOI.empty!(opt)
        @test MOI.is_empty(opt)

        model = QuadraticToBinary.Optimizer{Float64}(opt)

        x = MOI.add_variable(model)
        y = MOI.add_variable(model)

        MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
        MOI.add_constraint(model, y, MOI.GreaterThan(1.0))

        MOI.add_constraint(model, x, MOI.LessThan(10.0))
        MOI.add_constraint(model, y, MOI.LessThan(10.0))

        c = MOI.add_constraint(model, 1.0 * x * y, MOI.LessThan(4.0))

        MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            2.0 * x + y)
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

        MOI.optimize!(model)

        @assert MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL

        @assert MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @assert MOI.get(model, MOI.ObjectiveValue()) ≈ 9.0

        @assert MOI.get(model, MOI.VariablePrimal(), x) ≈ 4.0
        @assert MOI.get(model, MOI.VariablePrimal(), y) ≈ 1.0

        @assert MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 4.0
    end
end

@testset "MOI Deprecated Unit" begin
    for opt in optimizers
        MODEL = QuadraticToBinary.Optimizer{Float64}(opt)
        MOIT.unittest(MODEL, CONFIG, [
            "number_threads", # FIXME implement `MOI.NumberOfThreads`
            # INFEASIBLE_OR_UNBOUNDED instead of DUAL_INFEASIBLEÇ
            "solve_unbounded_model",
            # No quadratics (because of no bounds):
            "delete_soc_variables",
            "solve_qcp_edge_cases",
            "solve_qp_edge_cases",
            "solve_qp_zero_offdiag",
        ])
    end
end

@testset "ModelLike" begin
    # @test MOI.get(CACHED, MOI.SolverName()) == "COIN Branch-and-Cut (Cbc)"
    for opt in optimizers
        MODEL = QuadraticToBinary.Optimizer{Float64}(opt)
        @testset "default_objective_test" begin
            MOIT.default_objective_test(MODEL)
        end
        @testset "default_status_test" begin
            MOIT.default_status_test(MODEL)
        end
        @testset "nametest" begin
            MOIT.nametest(MODEL)
        end
        @testset "validtest" begin
            MOIT.validtest(MODEL)
        end
        @testset "emptytest" begin
            MOIT.emptytest(MODEL)
        end
        @testset "orderedindicestest" begin
            MOIT.orderedindicestest(MODEL)
        end
        @testset "copytest" begin
            # Requires VectorOfVariables
            # MOIT.copytest(MODEL, MOIU.CachingOptimizer(
            #     ModelForCachingOptimizer{Float64}(),
            #     Cbc.Optimizer()
            # ))
        end
    end
end

@testset "Continuous Linear" begin
    for opt in optimizers
        MODEL = QuadraticToBinary.Optimizer{Float64}(opt)
        MOIT.contlineartest(MODEL, CONFIG)
    end
end

@testset "Integer Linear" begin
    for opt in optimizers
        MODEL = QuadraticToBinary.Optimizer{Float64}(opt)
        MOIT.intlineartest(MODEL, CONFIG, [
            "indicator1", "indicator2", "indicator3", "indicator4",
            "int2",
        ])
    end
end

@testset "From MOI ncqcp" begin
    for opt in optimizers
        MODEL = QuadraticToBinary.Optimizer{Float64}(opt)
        ncqcp1test_mod(MODEL, CONFIG_LOW_TOL)
        ncqcp1test_mod2(MODEL, CONFIG_LOW_TOL)
        ncqcp2test_mod(MODEL, CONFIG_LOW_TOL)
    end
end

@testset "From MOI qcp" begin
    for opt in optimizers
        MODEL = QuadraticToBinary.Optimizer{Float64}(opt)
        MOI.set(MODEL, QuadraticToBinary.GlobalVariablePrecision(), 1e-4)
        qcp1test_mod(MODEL, CONFIG_LOW_TOL)
        qcp2test_mod(MODEL, CONFIG_LOW_TOL)
        qcp3test_mod(MODEL, CONFIG_LOW_TOL)
        qcp4test_mod(MODEL, CONFIG_LOW_TOL)
        qcp5test_mod(MODEL, CONFIG_LOW_TOL)
        #
        MOI.set(MODEL, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)
        socp1test_mod(MODEL, CONFIG_LOW_TOL)
    end
end

@testset "From MOI qp" begin
    for opt in optimizers
        MODEL = QuadraticToBinary.Optimizer{Float64}(opt)
        MOI.set(MODEL, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)
        qp1test_mod(MODEL, CONFIG_LOW_TOL)
        MOI.set(MODEL, QuadraticToBinary.GlobalVariablePrecision(), 1e-4)
        qp2test_mod(MODEL, CONFIG_VERY_LOW_TOL)
        qp3test_mod(MODEL, CONFIG_LOW_TOL)
    end
end


text:
"""
utils.jl defines utility methods.

------------------------------------------------------------------------------
COPYRIGHT/LICENSE. This file is part of the XYZ package. It is subject to
the license terms in the LICENSE file found in the top-level directory of
this distribution. No part of the XYZ package, including this file, may be
copied, modified, propagated, or distributed except according to the terms
contained in the LICENSE file.
------------------------------------------------------------------------------
"""
# --- Exports

export construct_ipc_url

# --- Method definitions

"""
    construct_ipc_url(ipc_dir::String, ipc_name::String)

Construct IPC URL from `ipc_dir` and `ipc_name`.
"""
function construct_ipc_url(ipc_dir::String, ipc_name::String)
    if occursin("ipc://", ipc_dir)
        if !startswith(ipc_dir, "ipc://")
            message = """Invalid `ipc_dir` "$ipc_dir"."""
            throw(ArgumentError(message))
        end
    elseif !startswith(ipc_dir, "ipc://")
        ipc_dir = "ipc://" * ipc_dir
    end

    return joinpath(ipc_dir, ipc_name)
end


text:
using DelaunayMaps
using Test

@testset "DelaunayMaps.jl" begin
    @testset "Tests for triangulation" begin
    @testset "Test function is_closer" begin
        points = [
            0  0
            1  0
            1  1
        ]
        @test DelaunayMaps.is_closer(points, 1, 2, 3) == true
        @test DelaunayMaps.is_closer(points, 1, 3, 2) == false
    end
    @testset "Test function is_aligned" begin
        points = [
            0  0
            1  0
            1  1
            2  0
        ]
        @test DelaunayMaps.is_aligned(points, 1, 2, 3) == false
        @test DelaunayMaps.is_aligned(points, 1, 2, 4) == true
    end
    @testset "Test function is_left_of" begin
        points = [
             0  0
             1  0
             1  1
             1 -1
            -1  1
            -1 -1
             2  0
        ]
        @test DelaunayMaps.is_left_of(points, 1, 2, 3) == true
        @test DelaunayMaps.is_left_of(points, 1, 2, 4) == false
        @test DelaunayMaps.is_left_of(points, 1, 2, 5) == true
        @test DelaunayMaps.is_left_of(points, 1, 2, 6) == false
        @test DelaunayMaps.is_left_of(points, 1, 2, 7) == true
    end
    @testset "Test function is_strictly_left_of" begin
        points = [
             0  0
             1  0
             1  1
             1 -1
            -1  1
            -1 -1
             2  0
        ]
        @test DelaunayMaps.is_strictly_left_of(points, 1, 2, 3) == true
        @test DelaunayMaps.is_strictly_left_of(points, 1, 2, 4) == false
        @test DelaunayMaps.is_strictly_left_of(points, 1, 2, 5) == true
        @test DelaunayMaps.is_strictly_left_of(points, 1, 2, 6) == false
        @test DelaunayMaps.is_strictly_left_of(points, 1, 2, 7) == false
    end
    @testset "Test function is_right_of" begin
        points = [
             0  0
             1  0
             1  1
             1 -1
            -1  1
            -1 -1
             2  0
        ]
        @test DelaunayMaps.is_right_of(points, 1, 2, 3) == false
        @test DelaunayMaps.is_right_of(points, 1, 2, 4) == true
        @test DelaunayMaps.is_right_of(points, 1, 2, 5) == false
        @test DelaunayMaps.is_right_of(points, 1, 2, 6) == true
        @test DelaunayMaps.is_right_of(points, 1, 2, 7) == true
    end
    @testset "Test function is_strictly_right_of" begin
        points = [
             0  0
             1  0
             1  1
             1 -1
            -1  1
            -1 -1
             2  0
        ]
        @test DelaunayMaps.is_strictly_right_of(points, 1, 2, 3) == false
        @test DelaunayMaps.is_strictly_right_of(points, 1, 2, 4) == true
        @test DelaunayMaps.is_strictly_right_of(points, 1, 2, 5) == false
        @test DelaunayMaps.is_strictly_right_of(points, 1, 2, 6) == true
        @test DelaunayMaps.is_strictly_right_of(points, 1, 2, 7) == false
    end

    @testset "Test of Triangle constructors" begin
        points = [
             0  0
             1  0
             1  1
        ]
        @testset "For two vertices" begin
            t = DelaunayMaps.Triangle(1, 2)
            @test t.v3 == 0
            @test t.t1 == t.t2 == t.t3
            next = t.t3
            @test next.v3 == 0
            @test next.v1 == t.v2
            @test next.v2 == t.v1
            @test next.t1 == next.t2 == next.t3 == t
        end
        @testset "For three vertices" begin
            t = DelaunayMaps.Triangle(points, 1, 2, 3)
            # v1 -> v2 -> v3 should be clockwise
            @test t.v1 == 1
            @test t.v2 == 3
            @test t.v3 == 2
            @test DelaunayMaps.is_left_of(points, t.v1, t.v3, t.v2)
            @test t.v1 == 1

            @test t.t1.t3 == t.t2.t3 == t.t3.t3 == t
            @test t.t1.t1 == t.t3.t2 == t.t2
            @test t.t1.t2 == t.t2.t1 == t.t3
            @test t.t2.t2 == t.t3.t1 == t.t1

            @test t.t2.v2 == t.t3.v1 == t.v1
            @test t.t1.v1 == t.t3.v2 == t.v2
            @test t.t1.v2 == t.t2.v1 == t.v3
        end
    end
    @testset "Test is_hull" begin
        points = [
             0  0
             1  0
             1  1
        ]
        t = DelaunayMaps.Triangle(points, 1, 2, 3)
        @test DelaunayMaps.is_hull(t.t1)
        @test !DelaunayMaps.is_hull(t)
    end
    @testset "Test hull navigation function" begin
        points = [
             0  0
             1  0
             1  1
        ]
        t = DelaunayMaps.Triangle(points, 1, 2, 3)
        @test DelaunayMaps.next(t.t1) == t.t2
        @test DelaunayMaps.prev(t.t1) == t.t3
    end
    @testset "Test find_base_tangente" begin
        points = [
             0  0
             1  0
             1  1
             2  0
             3  0
             3  1
        ]
        t1 = DelaunayMaps.Triangle(points, 1, 2, 3)
        t2 = DelaunayMaps.Triangle(points, 4, 5, 6)
        t3 = DelaunayMaps.Triangle(2, 3)
        t4 = DelaunayMaps.Triangle(4, 5)

        t_left, t_right = DelaunayMaps.find_base_tangente(points, t1.t1, t2.t2)
        @test (t_left.v1, t_right.v1) == (2, 4)
    end
    @testset "Test _is_above_base" begin
        @testset "Test _is_above_base_left" begin
            @testset "For two edges triangles" begin
                points = [
                    0 0
                    0 1
                    1 0
                ]
                t = DelaunayMaps.Triangle(1, 2)
                @test DelaunayMaps._is_above_base_left(t, points, 1, 3) == true
                points = [
                    0 -1
                    0 0
                    1 0
                ]
                t = DelaunayMaps.Triangle(2, 1)
                @test DelaunayMaps._is_above_base_left(t, points, 2, 3) == false
            end
            @testset "For classic triangles" begin
                points = [
                    0 0
                    1 2
                    2 1
                    3 0
                ]
                t = DelaunayMaps.Triangle(points, 1, 2, 3)
                @test DelaunayMaps._is_above_base_left(t, points, 1, 4) == true
                points = [
                    0 0
                    1 -1
                    1 -2
                    3 0
                ]
                t = DelaunayMaps.Triangle(points, 1, 2, 3)
                @test DelaunayMaps._is_above_base_left(t, points, 1, 4) == false
                points = [
                    0 1
                    1 0
                    2 0
                    3 0
                ]
                t = DelaunayMaps.Triangle(points, 1, 2, 3)
                @test DelaunayMaps._is_above_base_left(t, points, 3, 4) == true
                points = [
                    0 -1
                    1 0
                    2 0
                    3 0
                ]
                t = DelaunayMaps.Triangle(points, 1, 2, 3)
                @test DelaunayMaps._is_above_base_left(t, points, 3, 4) == false
                points = [
                    0 1
                    1 1
                    2 0
                    3 0
                ]
                t = DelaunayMaps.Triangle(points, 1, 2, 3)
                @test DelaunayMaps._is_above_base_left(t, points, 3, 4) == true
                points = [
                    0 -1
                    1 1
                    2 0
                    3 0
                ]
                t = DelaunayMaps.Triangle(points, 1, 2, 3)
                @test DelaunayMaps._is_above_base_left(t, points, 3, 4) == false
            end
        end
        @testset "Test _is_above_base_right" begin
            @testset "For two edges triangles" begin
                points = [
                    0 0
                    1 0
                    1 1
                ]
                t = DelaunayMaps.Triangle(2, 3)
                @test DelaunayMaps._is_above_base_right(t, points, 1, 2) == true
                points = [
                    0 0
                    1 0
                    1 -1
                ]
                t = DelaunayMaps.Triangle(3, 2)
                @test DelaunayMaps._is_above_base_right(t, points, 1, 2) == false
            end
            @testset "For classic triangles" begin
                points = [
                    0 0
                    1 3
                    2 1
                    3 0
                ]
                t = DelaunayMaps.Triangle(points, 2, 3, 4)
                @test DelaunayMaps._is_above_base_right(t, points, 1, 4) == true
                points = [
                    0 0
                    1 -3
                    2 -2
                    3 0
                ]
                t = DelaunayMaps.Triangle(points, 2, 3, 4)
                @test DelaunayMaps._is_above_base_right(t, points, 1, 4) == false
                points = [
                    0 0
                    1 1
                    2 0
                    3 1
                ]
                t = DelaunayMaps.Triangle(points, 2, 3, 4)
                @test DelaunayMaps._is_above_base_right(t, points, 1, 3) == true
                points = [
                    0 0
                    1 -1
                    2 0
                    3 -1
                ]
                t = DelaunayMaps.Triangle(points, 2, 3, 4)
                @test DelaunayMaps._is_above_base_right(t, points, 1, 3) == false
                points = [
                    0 0
                    1 0
                    2 1
                    3 1
                ]
                t = DelaunayMaps.Triangle(points, 2, 3, 4)
                @test DelaunayMaps._is_above_base_right(t, points, 1, 2) == true
                points = [
                    0 0
                    1 0
                    2 -1
                    3 -1
                ]
                t = DelaunayMaps.Triangle(points, 2, 3, 4)
                @test DelaunayMaps._is_above_base_right(t, points, 1, 2) == false
            end
        end
    end
    @testset "Test _next" begin
        @testset "For two edges triangles" begin
            t = DelaunayMaps.Triangle(1, 2)
            @test DelaunayMaps._next_left(t, 1) == t.t3
            @test DelaunayMaps._next_left(DelaunayMaps._next_left(t, 1), 1) == t
            @test DelaunayMaps._next_right(t, 1) == t.t3
            @test DelaunayMaps._next_right(DelaunayMaps._next_right(t, 1), 1) == t
        end
        @testset "For classic triangles" begin
            points = [
                0 0
                1 0
                1 1
            ]
            t = DelaunayMaps.Triangle(points, 1, 2, 3)
            @test DelaunayMaps._next_left(t, 1) == t.t3
            @test DelaunayMaps._next_left(t, 3) == t.t1
            @test DelaunayMaps._next_left(t, 2) == t.t2
            @test DelaunayMaps._next_right(t, 1) == t.t2
            @test DelaunayMaps._next_right(t, 3) == t.t3
            @test DelaunayMaps._next_right(t, 2) == t.t1
        end
    end

    @testset "I tested every corner case I could think of, but just in case, let's generate 100 graphs and test if they are delaunay 🙈" begin
        for i in 1:100
            m = Map(20, 20, 10)
            delaunay = DelaunayMaps.is_delaunay(m)
            @test delaunay == true
            if !delaunay
                break
            end
        end
    end
end
end


text:
### ============== ### ============== ###
##    Simple Vicksek Model  2D         ##
##    Martin Zumaya Hernandez          ##
##    EXAMPLE SIMULATION SCRIPT        ##
### ============== ### ============== ###

### ============ INCLUDE PACKAGES ============ ###

using CollectiveDynamics.SVM2D
using Plots

### ============ SYSTEM'S PARAMETERS ============ ###

N   = parse(Int, ARGS[1]) # number of particles
ρ   = parse(Float64, ARGS[2]) # density
η   = parse(Float64, ARGS[3]) # noise intensity
ω   = parse(Float64, ARGS[4]) # maximal angular velocity
φ   = parse(Float64, ARGS[5]) # particle field of view
T   = parse(Int, ARGS[6]) # integration time steps
rep = parse(Int, ARGS[7]) # ensemble index
viz = true # visualize steps

### =============== ### =============== ###

v0 = 0.2 # 1.0 # particle's speed
dt = 1.0 # integration timestep

### =============== ### =============== ###

M = 10 # number of boxes per dimension

# l = 0.5 # interaction range is double the distance each particle moves in one time step
#l = 0.1 # interaction range is ten times the distance each particle moves in one time step

# r0 = (v0 * dt) / l # interaction range
ratio = 1.03
# r0 = v0 / (ratio * ω) # interaction range
r0 = 1.0

L = sqrt(N / ρ) # size of box

### =============== ### =============== ###

cell_size = step(range(0., stop=L , length=M))

### ============ SYSTEM'S INITIALIZATION ============ ###

box = Box(L, M) # suqare box of size L and M partitions per dimension
flock = Flock(N, L, dt, v0, r0, η, ω, φ) # flock initialization

### ============== ### ============== ###

output_path = set_output_data_structure_vsk("SVM_2D", N, ρ)

pos_file = open(output_path * "/pos_$(rep).dat", "w+")
vel_file = open(output_path * "/vel_$(rep).dat", "w+")

times = [convert(Int, exp10(i)) for i in 0:T]

plt = !viz ? nothing : quiver(
    [ flock.pos[i][1] for i in 1:N ],
    [ flock.pos[i][2] for i in 1:N ],
    quiver=(
        [ flock.vel[i][1] for i in 1:N ],
        [ flock.vel[i][2] for i in 1:N ]
    ),
    legend=false
)
savefig(plt, output_path * "/step0.png")

### ============ TIME EVOLUTION ============ ###

for i in 1:(length(times) - 1)

    if i > 1

        for t in (times[i]+1):times[i+1]

            evolve_system(flock, box, cell_size)

            if t % times[i] == 0 || t % times[i-1] == 0
                println("//////// ", t)
                if viz
                    plt = quiver(
                        [ flock.pos[i][1] for i in 1:N ],
                        [ flock.pos[i][2] for i in 1:N ],
                        quiver=(
                            [ flock.vel[i][1] for i in 1:N ],
                            [ flock.vel[i][2] for i in 1:N ]
                        ),
                        legend=false
                    )
                    savefig(plt, output_path * "/step$(t).png")
                end
                write(pos_file, vcat(flock.pos...))
                write(vel_file, vcat(flock.vel...))
            end
        end

    else

        for t in (times[i]+1):times[i+1]

            evolve_system(flock, box, cell_size)

            if t % times[i] == 0
                println("//////// ", t)
                if viz
                    plt = quiver(
                        [ flock.pos[i][1] for i in 1:N ],
                        [ flock.pos[i][2] for i in 1:N ],
                        quiver=(
                            [ flock.vel[i][1] for i in 1:N ],
                            [ flock.vel[i][2] for i in 1:N ]
                        ),
                        legend=false
                    )
                    savefig(plt, output_path * "/step$(t).png")
                end
                write(pos_file, vcat(flock.pos...))
                write(vel_file, vcat(flock.vel...))
            end
        end

    end

end

close(pos_file)
close(vel_file)

println("Done all")


text:
module NWLLPriorityQueues
const DEBUG = false

### Node weight tracking data structure ###
#
# We track the weight (minimum degree, minimum fill, approximate minimum degree, or
# what have you) of each live supervariable (quotient graph node) in the following
# data structure. Specifically, we need be able to
# 1 --- Efficiently identify and remove a minimum-weight node from the data structure
#       in each gross iteration of the algorithm.
# 2 --- For each such gross iteration of the algorithm, we need efficiently update
#       the node weights of (typically several) nodes after eliminating the chosen
#       minimum-weight node from the quotient graph.
# 3 --- For each such gross iteration of the algorithm, we need efficiently remove
#       nodes (often several) that are merged into other nodes (superveriables)
#       due to indistinguishability.
# Additionally, once we've constructed the data structure, we never insert new nodes
# (outside of updating, where updating involves insertion).
#
# This seems to be the traditional data structure in the local preorderings literature

immutable NWLLPriorityQueue{T <: Signed}
    # Node weight linked-list priority queue
    order::T
    minweight::Vector{T}
    weights::Vector{T}
    fwdlinks::Vector{T}
    bwdlinks::Vector{T}
    headlinks::Vector{T}
    function NWLLPriorityQueue(order::T)
        minweight = T[order+1]
        fwdlinks = zeros(T, order) # FUTURE necessary to zero-initialize? no
        bwdlinks = zeros(T, order) # FUTURE necessary to zero-initialize? no
        headlinks = zeros(T, order) # FUTURE necessary to zero-initialize? yes
        weights = Vector{T}(order)
        new(order, minweight, weights, fwdlinks, bwdlinks, headlinks)
    end
end

@inline getnodeweight(pq::NWLLPriorityQueue, node) = pq.weights[node]
@inline setnodeweight!(pq::NWLLPriorityQueue, node, weight) = pq.weights[node] = weight

function enqueue!{T}(pq::NWLLPriorityQueue{T}, node::T, weight::T)
    fwdlink = pq.headlinks[weight+1]
    pq.headlinks[weight+1] = node
    if fwdlink != 0
        pq.bwdlinks[fwdlink] = node
    elseif pq.minweight[1] > weight
        pq.minweight[1] = weight
    end
    pq.fwdlinks[node] = fwdlink
    pq.bwdlinks[node] = 0
    pq.weights[node] = weight
end

function dequeue!{T}(pq::NWLLPriorityQueue{T}, node::T)
    weight = pq.weights[node]
    fwdlink = pq.fwdlinks[node]
    bwdlink = pq.bwdlinks[node]
    fwdlink != 0 && (pq.bwdlinks[fwdlink] = bwdlink)
    bwdlink != 0 ? (pq.fwdlinks[bwdlink] = fwdlink) : (pq.headlinks[weight+1] = fwdlink)
    if bwdlink == 0 && fwdlink == 0 && weight == pq.minweight[1]
        # FUTURE Get rid of unnnecessary branch in preceding comparison
        while true
            pq.minweight[1] += 1
            pq.minweight[1] == pq.order && break
            pq.headlinks[pq.minweight[1]+1] != 0 && break
        end
    end
    return node, weight
end
function dequeue!{T}(pq::NWLLPriorityQueue{T})
    node = pq.headlinks[pq.minweight[1]+1]
    dequeue!(pq, node)
end

@inline function update!{T}(pq::NWLLPriorityQueue{T}, node::T, weight::T)
    dequeue!(pq, node)
    enqueue!(pq, node, weight)
end
# FUTURE Write specialized update routine.

end


text:
include(joinpath(@__DIR__, "pedagogy_graph.jl"))
using Random

dims = [5, 5]
g = LG.grid(dims)

sfg = PRS.SinkFreeGraph(g)

rng = Random.MersenneTwister(123)
@time sample = PRS.generate_sample_prs(rng, sfg)

path = joinpath(@__DIR__, "output" "sink_free_graph", "sink_free_graph.pdf")
p = plot(sfg, sample, dims, path)


text:
function plot_results(sp, pv, w, d; s=1, stage_1=[:gci, :gco], stage_2=[:gci2, :gco2])
    plt_sto = plot()
    plt_invest = plot()
    plt = plot() # create separate plot for storage state of charge
    plot!(plt_invest, pv .* value(sp[1, :u_pv]), label="pv")
    plot!(plt_invest, w .* value(sp[1, :u_wind]), label="wind")
    plot!(plt_invest, d, label="demand")
    stor_flow_i = value.(sp[1, :sto_in]).data
    stor_flow_o = value.(sp[1, :sto_out]).data

    stor_charge = [-sum(stor_flow_i[1:t])+sum(stor_flow_o[1:t]) for t in value.(sp[1, :sto_in]).axes[1]] .+ 0.5*value(sp[1, :u_storage])
    plot!(plt_sto, value.(sp[1, :sto_in]).axes, stor_charge, label="global storage charge")

    for var in stage_1
        plot!(plt, value.(sp[1, var]).axes, value.(sp[1, var]).data, label=string(var))

    end
    for var in stage_2

        plot!(plt, value.(sp[2, var], s).axes, value.(sp[2, var], s).data, label=string(var)*string(s), linestyle=:dash, linewidth=2)
    end
    stor_flow_i = value.(sp[2, :sto_in2],s).data
    stor_flow_o = value.(sp[2, :sto_out2],s).data
    stor_charge = [-sum(stor_flow_i[1:t])+sum(stor_flow_o[1:t]) for t in value.(sp[2, :sto_in2], s).axes[1]] .+ 0.5*value(sp[1, :u_storage])
    plot!(plt_sto, value.(sp[2, :sto_in2], s).axes, stor_charge, label="stochastic storage charge in scenario $s", linestyle=:dash, linewidth=2)

    display(plot(plt_invest, plt, plt_sto, layout = (3,1)))
    #print(s)
end

function plot_flexibility(scen_results, ov)
    T = size(scen_results)[1]
    positive_flexibility = zeros(T ÷ 2)
    negative_flexibility = zeros(T ÷ 2)
    positive_potential = zeros(T ÷ 2)
    negative_potential = zeros(T ÷ 2)
    for i in 1:T
        if scen_results[i, 3] > 0
            if scen_results[i, 1] != Inf && scen_results[i, 1] != -Inf
                positive_flexibility[Int(scen_results[i, 2])] = scen_results[i, 1] - ov
                if scen_results != Inf
                    positive_potential[Int(scen_results[i, 2])] = scen_results[i, 3]
                end
            end
        else
            if scen_results[i, 1] != Inf && scen_results[i, 1] != -Inf
                negative_flexibility[Int(scen_results[i, 2])] = scen_results[i, 1] - ov
                if scen_results != -Inf
                    negative_potential[Int(scen_results[i, 2])] = scen_results[i, 3]
                end
            end
        end
    end
    plt_cost = plot()
    plt_pot = plot()
    plot!(plt_cost, positive_flexibility ./ positive_potential, label = "price of positive flexibility")
    plot!(plt_cost, negative_flexibility ./ negative_potential, label = "price of negative flexibility")
    plot!(plt_pot, positive_potential, fillrange = 0, fillalpha = 0.35, label = "positive flexibility potential")
    plot!(plt_pot, negative_potential, fillrange = 0, fillalpha = 0.35, label = "negative flexibility potential")
    display(plot(plt_cost, plt_pot, layout = (2, 1)))
    return positive_flexibility, negative_flexibility, positive_potential, negative_potential
en

function plot_difference(sp; s = 1)
    plt_gc = plot()
    plt_sto = plot()
    t_xi = scenarios(sp)[s].data.t_xi
    recovery_time = sp.stages[2].parameters[:recovery_time]

    plot!(plt_gc, -value.(sp[1, :gci])[t_xi+1:t_xi+recovery_time]+value.(sp[2, :gci2], s), xlabel = "time after the event, h", label = "gci2-gci")
    plot!(plt_gc, -value.(sp[1, :gco])[t_xi+1:t_xi+recovery_time]+value.(sp[2, :gco2], s), xlabel = "time after the event, h", label = "gco2-gco")

    stor_flow_i = value.(sp[1, :sto_in])[t_xi:t_xi+recovery_time]
    stor_flow_o = value.(sp[1, :sto_out])[t_xi:t_xi+recovery_time]
    stor_charge = [-sum(stor_flow_i[1:t]) + sum(stor_flow_o[1:t]) for t in 1:recovery_time+1] .+ 0.5 * value(sp[1, :u_storage])
    stor_flow_i2 = value.(sp[2, :sto_in2], s)
    stor_flow_o2 = value.(sp[2, :sto_out2], s)
    stor_charge2 = [-sum(stor_flow_i2[1:t]) + sum(stor_flow_o2[1:t]) for t in 1:(recovery_time+1)] .+ stor_charge[t_xi-1]

    plot!(plt_sto, stor_charge2-stor_charge, label = "soc2-soc")

    display(plot(plt_gc, plt_sto, layout = (2, 1)))
end


text:
# BitIndex
# --------
#
# Utils for indexing bits in a vector of unsigned integers (internal use only).
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md
#
#         index(i)-1        index(i)        index(i)+1
# ....|................|..X.............|................|....
#                          |<-offset(i)-|
#                      |<--- 64 bits -->|

"""
A type trait of any bitvector-like type, that specifies how many bits are
required per element.

In this case out bitvector-like types are BioSequences, as they store sequences
of biological symbols in a succinct, packed form.

This trait should be defined and overloaded for your own bitvector-like types
or own concrete BioSequence subtypes.

!!! tip
    In BioSequences.jl, this trait is defined as a trait of the sequence type.

    But as a generic fallback, a default `BitsPerElem` trait for a sequence type
    is determined by the sequence type's Alphabet type trait.

    For example if you made a new sequence type, with an `Alphabet` trait of
    `DNAAlphabet{4}` and did not overload a `BitsPerElem` method specifically
    for the type, calling `BitsPerElem` on a sequence of that type results in
    a fallback result of `BitsPerElem{4}`, as dictated by the `DNAAlphabet{4}`.

    This needs repeating: BitsPerElem is a trait of sequence types, not of
    alphabet types. The alphabet types determine default fallbacks which are
    typically used because they are reasonable, but they may be overruled by a
    developer writing a `BitsPerElem` overload specifically for their sequence
    type.
"""
struct BitsPerElem{N} end

BitsPerElem(::A) where {A<:NucleicAcidAlphabet{2}} = BitsPerElem{2}()
BitsPerElem(::A) where {A<:NucleicAcidAlphabet{4}} = BitsPerElem{4}()
BitsPerElem(s::BioSequence) = BitsPerElem(Alphabet(s))
BitsPerElem(k::Kmer) = BitsPerElem(Alphabet(k))

"""
A type trait of any bitvector-like type, that specifies whether packed elements
are stored in an unsigned integer in order from left to right, or from right to
left.

```
Left to right element order:

             1----------->N
MSB of UInt |.............X..| LSB of UInt

Right to left element order:

                  N<--------1
MSB of UInt |.....X..........| LSB of UInt
```
"""
abstract type ElementOrder end
struct LeftToRight <: ElementOrder end
struct RightToLeft <: ElementOrder end

ElementOrder(::Kmer{A,K,N}) where {A,K,N} = LeftToRight()

function registertype end

registertype(::Kmer{A,K,N}) where {A,K,N} = UInt64

struct BitIndex{N,W<:Unsigned,O<:ElementOrder}
    val::Int64
end

@inline function bitindex(::BitsPerElem{N}, ::Type{W}, ::O, i::Integer) where {N,W,O<:ElementOrder}
    return BitIndex{N,W,O}((i - 1) << trailing_zeros(N))
end

@inline function bitindex(x, i::Integer)
    return bitindex(BitsPerElem(x), registertype(x), ElementOrder(x), i)
end

"Compute the bitindex for the `i'th` base in `seq`'s bitvector-like datastore."
@inline function bitindex(seq::Kmer, i::Integer)
    return bitindex(BitsPerElem(seq), registertype(seq), ElementOrder(seq), i + n_unused(seq))
end

# TODO: Decide whether to define these in terms of firstindex and lastindex,
# instead of 1 and length(seq).
@inline firstbitindex(seq::Kmer) = bitindex(seq, 1 + n_unused(seq))
@inline lastbitindex(seq::Kmer)  = bitindex(seq, length(seq) + n_unused(seq))

@inline function eachbitindex(from::BitIndex{N,W,O}, to::BitIndex{N,W,O}) where {N,W,O}
    return from:N:to
end

@inline function eachbitindex(seq::Kmer, from = 1, to = length(seq))
    return eachbitindex(from, to)
end

index_shift(i::BitIndex{N,UInt64,O}) where {N,O} = 6
index_shift(i::BitIndex{N,UInt32,O}) where {N,O} = 5
index_shift(i::BitIndex{N,UInt16,O}) where {N,O} = 4
index_shift(i::BitIndex{N,UInt8,O})  where {N,O} = 3

offset_mask(i::BitIndex{N,W,O}) where {N,W,O} = UInt8(8 * sizeof(W)) - 0x01

"Get the index of the register the element pointed to by BitIndex `i`"
register(i::BitIndex) = (i.val >> index_shift(i)) + 1

"""
    offset(i::BitIndex{N,W,RightToLeft}) where {N,W}

Get the size of the right-shift of the register required to bring the element
pointed to by BitIndex `i`, to the `N` least significant digits.

```
               Before:
MSB of UInt |.....X..........| LSB of UInt
                  |--------->| Size of the right-shift required = offset
MSB of UInt |...............X| LSB of UInt
               After:
```
"""
@inline offset(i::BitIndex{N,W,RightToLeft}) where {N,W} = i.val & offset_mask(i)

"""
    offset(i::BitIndex{N,W,LeftToRight}) where {N,W}

Get the size of the right-shift of the register required to bring the element
pointed to by BitIndex `i`, to the `N` least significant digits.

```
               Before:
MSB of UInt |.....X..........| LSB of UInt
                  |--------->| Size of the right-shift required = offset
MSB of UInt |...............X| LSB of UInt
               After:
```
"""
@inline offset(i::BitIndex{N,W,LeftToRight}) where {N,W} = (8 * sizeof(W) - N) - (i.val & offset_mask(i))

"Create a bit mask covering the least significant `n` bits."
bitmask(bidx::BitIndex{N,W,O}) where {N,W,O} = (one(W) << N) - one(W)

Base.:+(i::BitIndex, n::Int) = typeof(i)(i.val + n)
Base.:-(i::BitIndex, n::Int) = typeof(i)(i.val - n)
Base.:-(i1::T, i2::T)     where {T<:BitIndex} = i1.val - i2.val
Base.:(==)(i1::T, i2::T)  where {T<:BitIndex} = i1.val == i2.val
Base.isless(i1::T, i2::T) where {T<:BitIndex} = isless(i1.val, i2.val)
Base.cmp(i1::T, i2::T)    where {T<:BitIndex} = cmp(i1.val, i2.val)

@inline function nextposition(i::BitIndex{N,W,O}) where {N,W,O}
    return i + N
end

@inline function prevposition(i::BitIndex{N,W,O}) where {N,W,O}
    return i - N
end

function Base.iterate(i::BitIndex, s = 1)
    if s == 1
        return (register(i), 2)
    elseif s == 2
        return (offset(i), 3)
    else
        return nothing
    end
end

Base.show(io::IO, i::BitIndex) = print(io, '(', register(i), ", ", offset(i), ')')

"Extract the element stored in a packed bitarray referred to by bidx."
@inline function extract_encoded_element(bidx::B, data) where {B<:BitIndex}
    @inbounds chunk = data[register(bidx)]
    offchunk = chunk >> offset(bidx)
    return offchunk & bitmask(bidx)
end


text:
struct SolverState{T}
    iteration::Int
    fnorm::T
    stepnorm::T
    metadata::Dict
end

function SolverState(i, fnorm)
    SolverState(Int(i), fnorm, oftype(fnorm, NaN), Dict())
end

function SolverState(i, fnorm, stepnorm)
    SolverState(Int(i), fnorm, stepnorm, Dict())
end

struct SolverTrace
    states::Vector{SolverState}
end

SolverTrace() = SolverTrace(Array{SolverState}(undef, 0))

function Base.show(io::IO, t::SolverState)
    @printf io "%6d   %14e   %14e\n" t.iteration t.fnorm t.stepnorm
    if !isempty(t.metadata)
        for (key, value) in t.metadata
            @printf io " * %s: %s\n" key value
        end
    end
    return
end

Base.push!(t::SolverTrace, s::SolverState) = push!(t.states, s)

Base.getindex(t::SolverTrace, i::Integer) = getindex(t.states, i)

function Base.setindex!(t::SolverTrace,
                        s::SolverState,
                        i::Integer)
    setindex!(t.states, s, i)
end

function Base.show(io::IO, t::SolverTrace)
    @printf io "Iter     f(x) inf-norm    Step 2-norm \n"
    @printf io "------   --------------   --------------\n"
    for state in t.states
        show(io, state)
    end
    return
end

function update!(tr::SolverTrace,
                 iteration::Integer,
                 fnorm::Real,
                 stepnorm::Real,
                 dt::Dict,
                 store_trace::Bool,
                 show_trace::Bool)
    ss = SolverState(iteration, fnorm, stepnorm, dt)
    if store_trace
        push!(tr, ss)
    end
    if show_trace
        show(ss)
    end
    return
end

mutable struct SolverResults{T,I<:AbstractArray{T},Z<:AbstractArray{T}}
    method::String
    initial_x::I
    zero::Z
    residual_norm::T
    iterations::Int
    x_converged::Bool
    xtol::T
    f_converged::Bool
    ftol::T
    trace::SolverTrace
    f_calls::Int
    g_calls::Int
end

function converged(r::SolverResults)
    return r.x_converged || r.f_converged
end

function Base.show(io::IO, r::SolverResults)
    @printf io "Results of Nonlinear Solver Algorithm\n"
    @printf io " * Algorithm: %s\n" r.method
    @printf io " * Starting Point: %s\n" string(r.initial_x)
    @printf io " * Zero: %s\n" string(r.zero)
    @printf io " * Inf-norm of residuals: %f\n" r.residual_norm
    @printf io " * Iterations: %d\n" r.iterations
    @printf io " * Convergence: %s\n" converged(r)
    @printf io "   * |x - x'| < %.1e: %s\n" r.xtol r.x_converged
    @printf io "   * |f(x)| < %.1e: %s\n" r.ftol r.f_converged
    @printf io " * Function Calls (f): %d\n" r.f_calls
    @printf io " * Jacobian Calls (df/dx): %d" r.g_calls
    return
end


text:
module SeqAdaptiveIS

using StatsBase, Random

# Struct for holding results of inference.
struct MCPosterior
    samples::AbstractMatrix
    logW::Vector
end

StatsBase.weights(P::MCPosterior) = softmax(P.logW)
# resample(P::MCPosterior, N::Int) = P.samples[rand(Categorical(weights(P)), N),:]
resample(rng::AbstractRNG, P::MCPosterior, N::Int) = P.samples[multicategorical_kitagawa(rng, weights(P), N),:]
resample(P::MCPosterior, N::Int) = resample(Random.GLOBAL_RNG, P, N)
Base.length(P::MCPosterior) = length(P.logW)
ess(P::MCPosterior) = eff_ss(weights(P))

# useful class for holding components of GMM (unlike native impl, cholesky is not calculated)
# Not especially useful outside of this package, so it is not exported by default.
struct GMMComp{T <: Real}
    pis::Vector{T}
    mus::AbstractMatrix{T}
    covs::AbstractArray{T}
end

# AMIS procedure (and GMM / wEM fitting)
include("inference.jl")

# Application to sequential problems
include("sequential.jl")

export amis, resample, ess, seq_amis

end # module

max_stars_repo_name: ornithos/SeqAdaptiveIS

max_issues_repo_name: ornithos/SeqAdaptiveIS

max_forks_repo_name: ornithos/SeqAdaptiveIS


text:
using Random
Random.seed!(12)

include("wiener.jl")
include("diffusion.jl")
include("euler.jl")
include("misc.jl")
include("VHK.jl")
include("lyap.jl")
include("guip.jl")
include("partialbridge.jl")
include("partialbridgenuH.jl")

include("linpro.jl")
include("linprobridge.jl")
include("timechange.jl")
include("uniformscaling.jl")
include("gamma.jl")
include("gaussian.jl")
include("bessel.jl")

include("onlinestat.jl")

include("expint.jl")

include("with_srand.jl") # run last


text:
using Plots
using Plots.PlotMeasures
using StatsPlots
using KernelDensity

function widen(x::Union{Tuple, Vector}, s::Number=0.1)
    #μ = (x[1] + x[2])/2
    w = x[2] - x[1]
    y = (x[1] - s*w, x[2] + s*w)
end
export widen
# ------------------------------------------------------------------------------------------------ #
#                          Plot recipe for Process types (e.g. label axes)                         #
# ------------------------------------------------------------------------------------------------ #

@recipe function f(P::Process; vars=1:size(P.X0)[1], transient=false, downsample=1, colormode=nothing)
    cmode = colormode
    linecolor --> :black
    markercolor --> :black

    if length(vars) == 1
        t = times(P, transient=transient)
        x = (t, timeseries(P, vars, transient=transient))
        seriestype --> :path
        xguide --> "t"
        yguide --> "x"
        label --> nothing
        title --> String(Symbol(P.process))
    else
        if length(vars) > 3
            vars = vars[1:3] # Can only plot in 3d
        end
        x = deepcopy(timeseries(P, vars, transient=transient))
        x = Tuple([x[:, i] for i in 1:size(x)[2]])
        if cmode != :velocity
            seriestype --> :scatter
            markersize --> 1
            markerstrokewidth --> 0
        end
        label --> nothing
    end
    if downsample != 1
        x = tuple([xx[1:downsample:end] for xx in x]...)
    end

    if cmode == :velocity
        vx = deepcopy(timeseries(P, transient=transient))
        vx = Tuple([vx[:, i] for i in 1:size(vx)[2]])
        velocity = sqrt.(sum([(r[2:end] .- collect(r[1:end-1])).^2 for r ∈ [vx[i] for i in 1:length(vx)]], dims=1)[1]) # So colours by velocity over all vars
        seriestype := :path
        line_z := velocity
        linealpha --> 0.1
        linewidth --> 1.0
        colorbar --> nothing
    elseif cmode == :fourth
        if length(vars) < 4
            @error "Cannot color by the fourth variable if there are not four variables specified"
        end
        thefourth = vec(x[4])
        vars = vars[1:3]
        seriestype := :path
        line_z := thefourth
        linealpha --> 0.1
        linewidth --> 1.0
        colorbar --> nothing
    elseif !isnothing(cmode) && cmode != :none
        @error "Not a supported cmode"
    end

    return x
end

include("Animations.jl")
include("Trajectories.jl")
include("Spectra.jl")
include("AMI.jl")
include("fourseas.jl")
include("SystemCharacterisation.jl")