text:
export draw_xaxis, draw_yaxis, draw_axes, draw_xtick, draw_ytick

"""
`draw_xaxis(x)` draws a thin arrow from the origin to `(x,0)`.

`draw_xaxis(x1,x2)` is equivalent to `draw_xaxis(x1); draw_xaxis(x2)`.

If no x-value is given, draw a pair of arrows from the origin
to values 10% larger than those returned by `xlims()`.
"""
function draw_xaxis(x::Real; opts...)
    draw_vector(x, 0; color = :black, linewidth = 0.5, legend = false, opts...)
end

function draw_xaxis(x1::Real, x2::Real; opts...)
    draw_xaxis(x1; opts...)
    draw_xaxis(x2; opts...)
end

function draw_xaxis(; opts...)
    x1, x2 = xlims()
    draw_xaxis(x1; opts...)
    draw_xaxis(x2; opts...)
end

"""
`draw_yaxis(y)` draws a thin arrow from the origin to `(0,y)`.

`draw_yaxis(y1,y2)` is equivalent to `draw_yaxis(x1); draw_yaxis(x2)`.

If no y-value is given, draw a pair of arrows from the origin
to values 10% larger than those returned by `ylims()`.
"""
function draw_yaxis(y::Real; opts...)
    draw_vector(0, y; color = :black, linewidth = 0.5, opts...)
end

function draw_yaxis(y1::Real, y2::Real; opts...)
    draw_yaxis(y1; opts...)
    draw_yaxis(y2; opts...)
end

function draw_yaxis(; opts...)
    y1, y2 = ylims()
    draw_yaxis(y1; opts...)
    draw_yaxis(y2; opts...)
end

"""
`draw_axes()` invokes `draw_xaxis()` and `draw_yaxis()`.
"""
function draw_axes(; opts...)
    draw_xaxis(; opts...)
    draw_yaxis(; opts...)
end

const _DEFAULT_TICK_LEN = 0.2

"""
`draw_xtick(x::Real,len)` draws a tick mark on the x-axis with total length `len`.
If `len` is omitted, use `SimpleDrawing._DEFAULT_TICK_LEN`.

`draw_xtick(xlist,len)` calls `draw_xtick` for each value in `xlist`.
"""
function draw_xtick(x::Real, len = _DEFAULT_TICK_LEN; opts...)
    draw_segment(x, -len / 2, x, len / 2, color = :black, linewidth = 0.5, opts...)
end

function draw_xtick(xlist, len = _DEFAULT_TICK_LEN; opts...)
    for x in xlist
        draw_xtick(x, len, opts...)
    end
    plot!()
end

"""
`draw_ytick(y::Real,len)` draws a tick mark on the y-axis with total length `len`.
If `len` is omitted, use `SimpleDrawing._DEFAULT_TICK_LEN`.

`draw_xtick(ylist,len)` calls `draw_ytick` for each value in `xyist`.
"""
function draw_ytick(y::Real, len = _DEFAULT_TICK_LEN; opts...)
    draw_segment(-len / 2, y, len / 2, y, color = :black, linewidth = 0.5, opts...)
end

function draw_ytick(ylist, len = _DEFAULT_TICK_LEN; opts...)
    for y in ylist
        draw_ytick(y, len, opts...)
    end
    plot!()
end


text:
function extrapolate_left!(f, r, rₘᵢₙ = 1e-7)
  m = 2
  for i in eachindex(r)
      if r[i] ≥ rₘᵢₙ
          m = i
          break
      end
  end
  tg = ((f[m-1] - f[m])*(r[m+1] - r[m])^2 - (f[m+1] - f[m])*(r[m-1] - r[m])^2) /
       ((r[m+1] - r[m])*(r[m-1] - r[m])*(r[m+1] - r[m-1]))
  f0 = f[m] - tg*r[m]
  for i=1:m-1
      f[i] = tg*r[i] + f0
  end
  return m, tg, f0
end

"""
    radial_shr_eq(V::Function = r -> -1/r, x::AbstractRange{T} = -30.0:0.1:5.0;
    conf = 1, μ = 1.0, α = 1e5) where {T}

    radial_shr_eq(V::Array{S}, x::AbstractRange{T};
    conf = 1, μ = 1.0, α = 1e5) where {S, T}

Solve Schrödinger equation ([atomic units](https://en.wikipedia.org/wiki/Hartree_atomic_units) are assumed)
```math
-\\frac{1}{2\\mu}\\left[\\frac{1}{r^{2}}\\frac{\\partial}{\\partial r}\\left(r^{2}\\frac{\\partial R_{nl}}{\\partial r}\\right)-\\frac{l(l+1)}{r^{2}}R_{nl}(r)\\right]+V(r)R_{nl}(r)=E_{nl}R_{nl}(r)
```
for [the particle in a spherically symmetric potential.](https://en.wikipedia.org/wiki/Particle_in_a_spherically_symmetric_potential)

## On input:

* `V` - potential, provided as an explicit function ``V(r)`` (by default `r -> -1/r`), or an array of values on a radial grid.
* `x` - uniform grid (`AbstractRange` object), such as ``r_i = \\exp(x_i)``. Default grid is `-30.0:0.1:5.0`
* `conf` - electronic configuration of interest, by default: ground state `conf = ((1.0),)`
* `μ` - [reduced mass](https://en.wikipedia.org/wiki/Reduced_mass), by default `μ = 1.0` (electron mass in a.u.)
* `α` - parameter ``α = 10^5`` is used in the matrix pencil ``\\mathrm{H} + α\\mathrm{S}`` for the generalized eigenproblem.

## On output function returns:

* sum of one-particle energies: ``E = \\sum_{i=1}^{N} n_i \\epsilon_i``
* particles density: ``\\rho(x) = \\frac{1}{4\\pi} \\sum_{i=1}^{N} n_i y_{i}^{2}(x)``
* orbitals ``y_i(x)``, corresponding eigenvalues ``\\epsilon_i`` (energy levels), azimuthal and radial quantum numbers ``l``, ``n_r``, and level populations ``n_i`` (as listed in `conf`).

## Examples
[Hydrogen atom](https://en.wikipedia.org/wiki/Hydrogen_atom)

```jldoctest
julia> radial_shr_eq(r -> -1/r, conf = conf_enc("1s1")).energy ≈ -1/2
true

julia> radial_shr_eq(r -> -1/r, conf = conf_enc("2s1")).energy ≈ -1/8
true

julia> radial_shr_eq(r -> -1/r, conf = conf_enc("2p1")).energy ≈ -1/8
true

julia> radial_shr_eq(r -> -1/r, conf = conf_enc("3s1")).energy ≈ -1/18
true

julia> radial_shr_eq(r -> -1/r, conf = conf_enc("3p1")).energy ≈ -1/18
true

julia> radial_shr_eq(r -> -1/r, conf = conf_enc("3d1")).energy ≈ -1/18
true
```
"""
function radial_shr_eq(V::AbstractArray,
                       x::AbstractRange = -35.0:0.1:5.0;
                    conf::Array{Array{T,1},1} = [[1]],
                       μ::Real = 1,
                       Α::Real = 1e5) where T <: Real

  if length(V) != length(x)
    throw(DimensionMismatch("V and r have different length"))
  end

  r, r², n, dx = exp.(x), exp.(2x), length(x), step(x)

  D = 1:n+1:n*n

  H = -1/2μ*laplacian(x); HD = H[D] .+ V .* r²
  S = copy(H)

  ε = similar(r); ρ = zero(r); ∑nᵢεᵢ = 0.0

  ψ = Dict()
  # the equation is solved separately for each subshell: s, p, d, f
  for l=0:length(conf)-1
    @. H[D] = HD + 1/2μ * (l + 1/2)^2
    @. @views S[D] = H[D] + Α * r²
    θ, y = eigen(Symmetric(H), Symmetric(S))
    @. ε = Α * θ / (1 - θ)
    for (nᵣ, nᵢ) in enumerate(conf[l+1])
      @views y[:, nᵣ] ./= sqrt(∫(dx, y[:, nᵣ] .^ 2 .* r²))
      @views ρ .+= nᵢ / 4π * y[:, nᵣ] .^ 2
      ∑nᵢεᵢ += nᵢ * ε[nᵣ]
      ψ[(nᵣ = nᵣ - 1, l = l)] = (ϵᵢ = ε[nᵣ], nᵢ = nᵢ, ψᵢ = y[:, nᵣ])
    end
  end
  @. ρ /= r
  extrapolate_left!(ρ, r)
  @. ρ *= r
  return (energy = ∑nᵢεᵢ, density = ρ, orbitals = ψ)
end
radial_shr_eq(V::Function = r -> -1/r, x::AbstractRange = -35.0:0.1:5.0; conf = 1, μ = 1, Α = 1e5) =
              radial_shr_eq(V.(exp.(x)), x, conf = conf, μ = μ, Α = Α)

ext: jl


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

"""
Utilities for capturing, filtering and presenting streams of log events.
Normally you don't need to import `Logging` to create log events; for this
the standard logging macros such as `@info` are already exported by `Base`
and available by default.
"""
module Logging

# Import the CoreLogging implementation into Logging as new const bindings.
# Doing it this way (rather than with import) makes these symbols accessible to
# tab completion.
for sym in [
    :LogLevel, :BelowMinLevel, :Debug, :Info, :Warn, :Error, :AboveMaxLevel,
    :AbstractLogger,
    :NullLogger,
    :handle_message, :shouldlog, :min_enabled_level, :catch_exceptions,
    Symbol("@debug"),
    Symbol("@info"),
    Symbol("@warn"),
    Symbol("@error"),
    Symbol("@logmsg"),
    :with_logger,
    :current_logger,
    :global_logger,
    :disable_logging,
    :SimpleLogger]
    @eval const $sym = Base.CoreLogging.$sym
end

export
    AbstractLogger,
    LogLevel,
    NullLogger,
    @debug,
    @info,
    @warn,
    @error,
    @logmsg,
    with_logger,
    current_logger,
    global_logger,
    disable_logging,
    SimpleLogger,
    ConsoleLogger

include("ConsoleLogger.jl")

# The following are also part of the public API, but not exported:
#
# 1. Log levels:
#   BelowMinLevel, Debug, Info, Warn, Error, AboveMaxLevel,
#
# 2. AbstractLogger message related functions:
#  handle_message, shouldlog, min_enabled_level, catch_exceptions,

function __init__()
    global_logger(ConsoleLogger(stderr))
end

end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
#
# This file is part of the Simulate.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

# @static if VERSION >= v"1.3"
#     using Base.Threads: @spawn
# else
#     @eval const $(Symbol("@spawn")) = $(Symbol("@async"))
# end

"""
    register!(clk::Clock, p::Prc)

Register a Prc to a clock. Check its id and change it apropriately.
"""
function register!(clk::Clock, p::Prc)
    id = p.id
    while haskey(clk.processes, id)
        if isa(id, Float64)
            id = nextfloat(id)
        elseif isa(id, Int)
            id += 1
        elseif isa(id, String)
            s = split(id, "#")
            id = length(s) > 1 ? chop(id)*string(s[end][end]+1) : id*"#1"
        else
            throw(ArgumentError("process id $id is duplicate, cannot convert!"))
        end
    end
    clk.processes[id] = p
    p.id = id
end

"""
    loop(p::Prc, start::Channel, cycles::Number)

Put a [`Prc`](@ref) in a loop which can be broken by a `ClockException`.

# Arguments
- `p::Prc`:
- `start::Channel`: a channel to ensure that a process starts,
- `cycles=Inf`: determine, how often the loop should be run.
"""
function loop(p::Prc, start::Channel, cycles::Number)
    take!(start)
    if threadid() > 1
        p.clk = pclock(p.clk, threadid())
    end

    while cycles > 0
        try
            p.f(p.clk, p.arg...; p.kw...)
        catch exc
            if isa(exc, ClockException)
                exc.ev == Stop() && break
            end
            rethrow(exc)
        end
        cycles -= 1
    end
    p.clk.processes = delete!(p.clk.processes, p.id)
end

"""
    startup!(p::Prc, cycles::Number, spawn::Bool)

Start a `Prc` as a task in a loop.
"""
function startup!(c::AbstractClock, p::Prc, cycles::Number, spawn::Bool)

    function startit()
        start = Channel{Int}(0)
        p.task = @async loop(p, start, cycles)
        put!(start, 1)  # block until the process has started
    end

    if spawn
        if threadid() == 1
            @assert !isempty(c.ac) "no parallel clocks available!"
            sid = spawnid(c)
            if sid == 1
                startit()
            else
                talk(c, sid, Register((p=p, cycles=cycles)) )
            end
        else
            @threads for i in 1:nthreads()
                if threadid() == 1
                    process!(c.master[], p, cycles, spawn=spawn)
                end
            end
        end
    else
        startit()
    end
    register!(p.clk, p)
end

"""
```
process!([clk::Clock], p::Prc, cycles=Inf; spawn::Bool=false)
```
Register a [`Prc`](@ref) to a clock, start it as an asynchronous process and
return the `id` it was registered with. It can then be found under `clk.processes[id]`.

# Arguments
- `c::AbstractClock`: `Clock` or `ActiveClock`, if not provided, the process runs
    under `𝐶`,
- `p::Prc`: it contains a function and its arguments,
- `cycles::Number=Inf`: number of cycles the process should run,
- `spawn::Bool=false`: if true, the process may be scheduled on another thread
    in parallel and registered to the thread specific clock.

!!! note
    `spawn`ing a process is possible only with parallel clocks setup with
    [`PClock`](@ref) or [`fork!`](@ref).

"""
function process!(c::AbstractClock, p::Prc, cycles::Number=Inf; spawn::Bool=false)
    p.clk = c
    startup!(c, p, cycles, spawn)
    p.id
end
process!(p::Prc, cycles=Inf) = process!(𝐶, p, cycles)

"wakeup a process waiting for a `Condition`"
wakeup(c::Condition) = (notify(c); yield())

"""
```
delay!(clk::Clock, t::Number)
```
Delay a process for a time interval `t` on the clock `clk`. Suspend the calling
process until being reactivated by the clock at the appropriate time.

# Arguments
- `clk::Clock`: if not provided, the delay goes to `𝐶`.
- `t::Number`: the time interval for the delay.
"""
function delay!(clk::Clock, t::Number)
    c = Condition()
    event!(clk, Fun(wakeup, c), after, t)
    wait(c)
    # c = Channel{Int}()
    # event!(clk, Fun(wakeup, c), after, t)
    # take!(c)
end

"""
```
delay!(clk::Clock, T::Timing, t::Number)
```

Used for delaying a process *until* a given time t.

# Arguments
- `clk::Clock`: if no clock is given, the delay goes to 𝐶,
- `T::Timing`: only `until` is accepted,
- `t::Number`: delay until time t if t > clk.time, else give a warning.
"""
function delay!(clk::Clock, T::Timing, t::Number)
    @assert T == until "bad Timing $T for delay!"
    if t > clk.time
        c = Condition()
        event!(clk, Fun(wakeup, c), t)
        wait(c)
    else
        now!(clk, Fun(println, stderr, "warning: delay until $t ≤ τ=$(tau(clk))"))
    end
end

"""
```
wait!(clk::Clock, cond::Action; scope::Module=Main)
```
Wait on a clock for a condition to become true. Suspend the calling process
until the given condition is true.

# Arguments
- `clk::Clock`: if no clock is supplied, the delay goes to `𝐶`,
- `cond::Action`: a condition is an expression or function
    or an array or tuple of them. It is true only if all expressions or functions
    therein return true,
- `scope::Module=Main`: evaluation scope for given expressions.
"""
function wait!(clk::Clock, cond::Action; scope::Module=Main)
    if all(evExec(cond))   # all conditions met
        return         # return immediately
    else
        c = Condition()
        event!(clk, Fun(wakeup, c), cond, scope=scope)
        wait(c)
    end
end

"""
    interrupt!(p::Prc, ev::ClockEvent, value=nothing)

Interrupt a `Prc` by throwing a `ClockException` to it.
"""
function interrupt!(p::Prc, ev::ClockEvent, value=nothing)
    schedule(p.task, ClockException(ev, value), error=true)
    yield()
end

"Stop a Prc"
stop!(p::Prc, value=nothing) = interrupt!(p, Stop(), value)

"""
```
now!(clk::Clock, op::Action)
```
Tell the clock to execute an operation. Thus it cannot proceed before the op is finished.

!!! note
    This is needed for IO-operations of tasks. IO-operations yield the task to
    the scheduler and the scheduler may invoke the clock before giving control
    back to the task. In that case the clock will proceed and the task has gone
    out of sync with the clock. Use `now!` to avoid this situation!

# Arguments
- `clk::Clock`,
- `op::Action`: operation to execute.
"""
now!(clk::Clock, ex::Action) = event!(clk, ex, clk.time)


text:
export ExplicitSolver, ExplicitSolverTime, GeneralizedAlphaSolver,
GeneralizedAlphaSolverTime, ExplicitStaticSolver, ImplicitStaticSolver,
compute_boundary_info, compute_external_force

@doc raw"""
    ExplicitSolverTime(Δt::Float64, NT::Int64)

Returns the times for explicit solver. Boundary conditions and external forces should be given at these times.
"""
function ExplicitSolverTime(Δt::Float64, NT::Int64)
    U = zeros(NT)
    for i = 1:NT
        U[i] = i*Δt
    end
    U
end

@doc raw"""
    ExplicitSolver(globdat::GlobalData, domain::Domain,
        d0::Union{Array{Float64, 1}, PyObject},
        v0::Union{Array{Float64, 1}, PyObject},
        a0::Union{Array{Float64, 1}, PyObject},
        Δt::Float64, NT::Int64,
        H::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
        Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
        ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
        abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; strain::String = "small")

Differentiable Explicit Solver.

- `d0`, `v0`, `a0`: initial **full** displacement, velocity, and acceleration.

- `Δt`: time step

- `Hs`: linear elasticity matrix at each Gauss point

- `Fext`: external force, $\mathrm{NT}\times n$, where $n$ is the active dof.
   The external force includes all **body forces**, **external load forces** (also called **edge forces** in NNFEM) and **boundary acceleration-induced forces**.

- `ubd`, `abd`: boundary displacementt and acceleration, $\mathrm{NT}\times m$, where $m$ is **time-dependent** boundary DOF.
   Time-independent boundary conditions are extracted from `domain`.

- `strain_type` (default = "small"): small strain or finite strain
"""
function ExplicitSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject},
    v0::Union{Array{Float64, 1}, PyObject},
    a0::Union{Array{Float64, 1}, PyObject},
    Δt::Float64, NT::Int64,
    H::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
    Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    abd::Union{Array{Float64, 2}, PyObject, Missing}=missing;
    strain_type::String = "small")
    assembleMassMatrix!(globdat, domain)
    if !(strain_type in ["small", "finite"])
        error("Only small strain or finite strain are supported. Unknown strain type $strain_type")
    end
    init_nnfem(domain)
    M = factorize(constant(globdat.M))
    bddof = findall(domain.EBC[:] .== -2)
    fixed_bddof = findall(domain.EBC[:] .== -1)

    Fext, ubd, abd, H = convert_to_tensor([Fext, ubd, abd, H], [Float64, Float64, Float64, Float64])

    function condition(i, tas...)
        i<=NT
    end
    function body(i, tas...)
        d_arr, v_arr, a_arr = tas
        u, ∂u, ∂∂u = read(d_arr, i), read(v_arr, i), read(a_arr, i)

        u +=  Δt*∂u + 0.5*Δt*Δt*∂∂u
        ∂u += 0.5*Δt * ∂∂u

        if !ismissing(abd)
            u = scatter_update(u, bddof, ubd[i])
        end
        if length(sum(fixed_bddof))>0
            u = scatter_update(u, fixed_bddof, d0[fixed_bddof])
        end

        if strain_type=="small"
            ε = s_eval_strain_on_gauss_points(u, domain)
        else
            ε = f_eval_strain_on_gauss_points(u, domain)
        end
        if length(size(H))==2
            σ = tf.matmul(ε, H)
        else
            σ = batch_matmul(H, ε)
        end
        if strain_type=="small"
            fint  = s_compute_internal_force_term(σ, domain)
        else
            fint  = f_compute_internal_force_term(σ, u, domain)
        end
        if ismissing(Fext)
            fext = zeros(length(fint))
        else
            fext = Fext[i]
        end
        ∂∂up = vector(findall(domain.dof_to_eq), M\(fext - fint), domain.nnodes*2)

        if !ismissing(abd)
            ∂∂up = scatter_update(∂∂up, bddof, abd[i])
        end

        ∂u += 0.5 * Δt * ∂∂up

        i+1, write(d_arr, i+1, u), write(v_arr, i+1, ∂u), write(a_arr, i+1, ∂∂up)
    end

    arr_d = TensorArray(NT+1); arr_d = write(arr_d, 1, d0)
    arr_v = TensorArray(NT+1); arr_v = write(arr_v, 1, v0)
    arr_a = TensorArray(NT+1); arr_a = write(arr_a, 1, a0)
    i = constant(1, dtype=Int32)
    tas = [arr_d, arr_v, arr_a]
    _, d, v, a = while_loop(condition, body, [i, tas...])
    d, v, a = stack(d), stack(v), stack(a)
    sp = (NT+1, 2domain.nnodes)
    set_shape(d, sp), set_shape(v, sp), set_shape(a, sp)
end

"""
    compute_boundary_info(domain::Domain, globdat::GlobalData, ts::Array{Float64})

Computes the boundary information `ubd` and `abd`
"""
function compute_boundary_info(domain::Domain, globdat::GlobalData, ts::Array{Float64})
    bd = domain.EBC[:] .== -2
    if sum(bd)==0
        return missing, missing
    end

    @assert globdat.EBC_func!=nothing

    ubd = zeros(length(ts), sum(bd))
    abd = zeros(length(ts), sum(bd))
    for i = 1:length(ts)
        time = ts[i]
        ubd[i,:], _, abd[i,:] = globdat.EBC_func(time) # user defined time-dependent boundary
    end

    return ubd, abd
end

"""
    compute_external_force(domain::Domain, globdat::GlobalData, ts::Union{Array{Float64}, Missing} = missing)

Computes the external force (body force, edge force and force due to boundary acceleration).
"""
function compute_external_force(domain::Domain, globdat::GlobalData, ts::Union{Array{Float64}, Missing} = missing)
    if ismissing(ts)
        globdat.time = 0.0
        return getExternalForce(domain, globdat)
    end

    NT = length(ts)
    fext = zeros(NT, domain.neqs)
    for i = 1:length(ts)
        time = ts[i]
        globdat.time = time
        fext[i, :] = getExternalForce(domain, globdat)
    end
    globdat.time = 0.0
    fext
end
compute_external_force(globdat::GlobalData, domain::Domain, ts::Union{Array{Float64}, Missing} = missing) = compute_external_force(domain, globdat, ts)

@doc raw"""
    GeneralizedAlphaSolverTime(Δt::Float64, NT::Int64;ρ::Float64 = 0.0)

Returns the times for the generalized $\alpha$ solver. Boundary conditions and external forces should be given at these times.
"""
function GeneralizedAlphaSolverTime(Δt::Float64, NT::Int64;ρ::Float64 = 0.0)
    U = zeros(NT)
    @assert 0<=ρ<=1
    αm = (2ρ-1)/(1+ρ)
    αf = ρ/(1+ρ)
    β2 = 0.5*(1 - αm + αf)^2
    γ = 0.5 - αm + αf
    t = 0
    for i = 1:NT
        t += (1 - αf)*Δt
        U[i] = t
    end
    U
end

@doc raw"""
    GeneralizedAlphaSolver(globdat::GlobalData, domain::Domain,
        d0::Union{Array{Float64, 1}, PyObject},
        v0::Union{Array{Float64, 1}, PyObject},
        a0::Union{Array{Float64, 1}, PyObject},
        Δt::Float64, NT::Int64,
        Hs::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
        Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
        ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
        abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; ρ::Float64 = 0.0)

Differentiable Generalized $\alpha$ scheme. This is an extension of [`αscheme`](https://kailaix.github.io/ADCME.jl/dev/alphascheme/)
provided in ADCME. This function does not support damping and variable time step (for efficiency).

- `d0`, `v0`, `a0`: initial **full** displacement, velocity, and acceleration.

- `Δt`: time step

- `Hs`: linear elasticity matrix at each Gauss point

- `Fext`: external force, $\mathrm{NT}\times n$, where $n$ is the active dof.
  The external force includes all **body forces**, **external load forces** (also called **edge forces** in NNFEM) and **boundary acceleration-induced forces**.

- `ubd`, `abd`: boundary displacementt and acceleration, $\mathrm{NT}\times m$, where $m$ is boundary DOF.
  Time-independent boundary conditions are extracted from `domain`.

`GeneralizedAlphaSolver` does not support finite-strain continuum yet.
"""
function GeneralizedAlphaSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject},
    v0::Union{Array{Float64, 1}, PyObject},
    a0::Union{Array{Float64, 1}, PyObject},
    Δt::Float64, NT::Int64,
    Hs::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
    Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; ρ::Float64 = 0.0)
    assembleMassMatrix!(globdat, domain)
    @assert 0<=ρ<=1
    init_nnfem(domain)
    αm = (2ρ-1)/(1+ρ)
    αf = ρ/(1+ρ)
    β2 = 0.5*(1 - αm + αf)^2
    γ = 0.5 - αm + αf

    Fext, ubd, abd, Hs = convert_to_tensor([Fext, ubd, abd, Hs], [Float64, Float64, Float64, Float64])
    M = constant(globdat.M)
    stiff = s_compute_stiffness_matrix(Hs, domain)
    A = M*(1 - αm) + (1 - αf) * 0.5 * β2 * Δt^2 * stiff
    A = factorize(A)
    bddof = findall(domain.EBC[:] .== -2)
    fixed_bddof = findall(domain.EBC[:] .== -1)
    nbddof = findall(domain.dof_to_eq)

    function condition(i, tas...)
        i<=NT
    end
    function body(i, tas...)
        d_arr, v_arr, a_arr = tas
        u, ∂u, ∂∂u = read(d_arr, i), read(v_arr, i), read(a_arr, i)
        if ismissing(Fext)
            fext = zeros(domain.neqs)
        else
            fext = Fext[i]
        end
        ∂∂up = ∂∂u
        up =  (1 - αf)*(u + Δt*∂u + 0.5 * Δt * Δt * ((1 - β2)*∂∂u + β2*∂∂up)) + αf*u
        if !ismissing(abd)
            up = scatter_update(up, bddof, ubd[i])
        end
        if length(sum(fixed_bddof))>0
            up = scatter_update(up, fixed_bddof, d0[fixed_bddof])
        end

        ε = s_eval_strain_on_gauss_points(up, domain)
        if length(size(Hs))==2
            σ = tf.matmul(ε, Hs)
        else
            # @info Hs, ε
            σ = batch_matmul(Hs, ε)
        end
        fint  = s_compute_internal_force_term(σ, domain)
        res = M * (∂∂up[nbddof] *(1 - αm) + αm*∂∂u[nbddof])  + fint - fext
        Δ = -(A\res)
        ∂∂up= scatter_add(∂∂up, nbddof, Δ)
        if !ismissing(abd)
            ∂∂up = scatter_update(∂∂up, bddof, abd[i])
        end

        # updaet
        u += Δt * ∂u + Δt^2/2 * ((1 - β2) * ∂∂u + β2 * ∂∂up)
        ∂u += Δt * ((1 - γ) * ∂∂u + γ * ∂∂up)

        if length(sum(fixed_bddof))>0
            u = scatter_update(u, fixed_bddof, d0[fixed_bddof])
        end

        i+1, write(d_arr, i+1, u), write(v_arr, i+1, ∂u), write(a_arr, i+1, ∂∂up)
    end

    arr_d = TensorArray(NT+1); arr_d = write(arr_d, 1, d0)
    arr_v = TensorArray(NT+1); arr_v = write(arr_v, 1, v0)
    arr_a = TensorArray(NT+1); arr_a = write(arr_a, 1, a0)
    i = constant(1, dtype=Int32)
    tas = [arr_d, arr_v, arr_a]
    _, d, v, a = while_loop(condition, body, [i, tas...])
    d, v, a = stack(d), stack(v), stack(a)
    sp = (NT+1, 2domain.nnodes)
    set_shape(d, sp), set_shape(v, sp), set_shape(a, sp)
end

@doc raw"""
    GeneralizedAlphaSolver(globdat::GlobalData, domain::Domain,
        d0::Union{Array{Float64, 1}, PyObject},
        v0::Union{Array{Float64, 1}, PyObject},
        a0::Union{Array{Float64, 1}, PyObject},
        Δt::Float64, NT::Int64,
        Cs::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
        Hs::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
        Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
        ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
        abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; ρ::Float64 = 0.0)

Solve linear dynamical structural problem using a generalized alpha solver.
`Cs` and `Hs` are the corresponding linear viscosity and linear elasticity matrix.
"""
function GeneralizedAlphaSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject},
    v0::Union{Array{Float64, 1}, PyObject},
    a0::Union{Array{Float64, 1}, PyObject},
    Δt::Float64, NT::Int64,
    Cs::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
    Hs::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
    Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; ρ::Float64 = 0.0)
    assembleMassMatrix!(globdat, domain)
    @assert 0<=ρ<=1
    init_nnfem(domain)
    αm = (2ρ-1)/(1+ρ)
    αf = ρ/(1+ρ)
    β2 = 0.5*(1 - αm + αf)^2
    γ = 0.5 - αm + αf

    Cs, Fext, ubd, abd, Hs = convert_to_tensor([Cs, Fext, ubd, abd, Hs], [Float64, Float64, Float64, Float64, Float64])
    M = constant(globdat.M)
    if length(size(Hs))==2
        Hs = reshape(repeat(Hs, getNGauss(domain)), (getNGauss(domain), 3, 3))
    end
    if length(size(Cs))==2
        Cs = reshape(repeat(Cs, getNGauss(domain)), (getNGauss(domain), 3, 3))
    end

    stiff = s_compute_stiffness_matrix(Hs, domain)
    damp = s_compute_stiffness_matrix(Cs, domain)
    A = M*(1 - αm) + (1-αf)*Δt*γ*damp +  (1 - αf) * 0.5 * β2 * Δt^2 * stiff
    A = factorize(A)
    bddof = findall(domain.EBC[:] .== -2)
    fixed_bddof = findall(domain.EBC[:] .== -1)
    nbddof = findall(domain.dof_to_eq)

    function condition(i, tas...)
        i<=NT
    end
    function body(i, tas...)
        d_arr, v_arr, a_arr = tas
        u, ∂u, ∂∂u = read(d_arr, i), read(v_arr, i), read(a_arr, i)
        if ismissing(Fext)
            fext = zeros(domain.neqs)
        else
            fext = Fext[i]
        end
        ∂∂up = ∂∂u
        up =  (1 - αf)*(u + Δt*∂u + 0.5 * Δt * Δt * ((1 - β2)*∂∂u + β2*∂∂up)) + αf*u
        if !ismissing(abd)
            up = scatter_update(up, bddof, ubd[i])
        end
        if length(sum(fixed_bddof))>0
            up = scatter_update(up, fixed_bddof, d0[fixed_bddof])
        end

        ε = s_eval_strain_on_gauss_points(up, domain)
        if length(size(Hs))==2
            σ = tf.matmul(ε, Hs)
        else
            # @info Hs, ε
            σ = batch_matmul(Hs, ε)
        end
        fint  = s_compute_internal_force_term(σ, domain)

        vn = ∂u + Δt * (1- γ)*∂∂u
        vf = (1-αf)*vn + αf*∂u
        dotε = s_eval_strain_on_gauss_points(vf, domain)
        if length(size(Cs))==2
            visc = tf.matmul(dotε, Cs)
        else
            visc = batch_matmul(Cs, dotε)
        end
        fv = s_compute_internal_force_term(visc, domain)

        res = M * (∂∂up[nbddof] *(1 - αm) + αm*∂∂u[nbddof])  + fv + fint - fext
        Δ = -(A\res)
        ∂∂up= scatter_add(∂∂up, nbddof, Δ)
        if !ismissing(abd)
            ∂∂up = scatter_update(∂∂up, bddof, abd[i])
        end

        # updaet
        u += Δt * ∂u + Δt^2/2 * ((1 - β2) * ∂∂u + β2 * ∂∂up)
        ∂u += Δt * ((1 - γ) * ∂∂u + γ * ∂∂up)

        if length(sum(fixed_bddof))>0
            u = scatter_update(u, fixed_bddof, d0[fixed_bddof])
        end

        i+1, write(d_arr, i+1, u), write(v_arr, i+1, ∂u), write(a_arr, i+1, ∂∂up)
    end

    arr_d = TensorArray(NT+1); arr_d = write(arr_d, 1, d0)
    arr_v = TensorArray(NT+1); arr_v = write(arr_v, 1, v0)
    arr_a = TensorArray(NT+1); arr_a = write(arr_a, 1, a0)
    i = constant(1, dtype=Int32)
    tas = [arr_d, arr_v, arr_a]
    _, d, v, a = while_loop(condition, body, [i, tas...])
    d, v, a = stack(d), stack(v), stack(a)
    sp = (NT+1, 2domain.nnodes)
    set_shape(d, sp), set_shape(v, sp), set_shape(a, sp)
end

############################# NN based constitutive models #############################

@doc raw"""
    ExplicitSolver(globdat::GlobalData, domain::Domain,
        d0::Union{Array{Float64, 1}, PyObject},
        v0::Union{Array{Float64, 1}, PyObject},
        a0::Union{Array{Float64, 1}, PyObject},
        Δt::Float64, NT::Int64,
        nn::Function,
        Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
        ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
        abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; strain_type::String = "small"))

Similar to [`ExplicitSolver`](@ref); however, the constituve relation from $\epsilon$ to $\sigma$ must be provided by
the function `nn`.
"""
function ExplicitSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject},
    v0::Union{Array{Float64, 1}, PyObject},
    a0::Union{Array{Float64, 1}, PyObject},
    Δt::Float64, NT::Int64,
    nn::Function,
    Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; strain_type::String = "small")
    assembleMassMatrix!(globdat, domain)
    if !(strain_type in ["small", "finite"])
        error("Only small strain or finite strain are supported. Unknown strain type $strain_type")
    end

    init_nnfem(domain)
    M = factorize(constant(globdat.M))
    bddof = findall(domain.EBC[:] .== -2)
    fixed_bddof = findall(domain.EBC[:] .== -1)

    Fext, ubd, abd = convert_to_tensor([Fext, ubd, abd], [Float64, Float64, Float64])

    function condition(i, tas...)
        i<=NT
    end
    function body(i, tas...)
        d_arr, v_arr, a_arr = tas
        u, ∂u, ∂∂u = read(d_arr, i), read(v_arr, i), read(a_arr, i)

        u +=  Δt*∂u + 0.5*Δt*Δt*∂∂u
        ∂u += 0.5*Δt * ∂∂u

        if !ismissing(abd)
            u = scatter_update(u, bddof, ubd[i])
        end

        if length(sum(fixed_bddof))>0
            u = scatter_update(u, fixed_bddof, d0[fixed_bddof])
        end

        if strain_type=="small"
            ε = s_eval_strain_on_gauss_points(u, domain)
        else
            ε = f_eval_strain_on_gauss_points(u, domain)
        end

        σ = nn(ε)

        if strain_type=="small"
            fint  = s_compute_internal_force_term(σ, domain)
        else
            fint  = f_compute_internal_force_term(σ, u, domain)
        end
        if ismissing(Fext)
            fext = zeros(length(fint))
        else
            fext = Fext[i]
        end
        ∂∂up = vector(findall(domain.dof_to_eq), M\(fext - fint), domain.nnodes*2)

        if !ismissing(abd)
            ∂∂up = scatter_update(∂∂up, bddof, abd[i])
        end
        if length(sum(fixed_bddof))>0
            ∂∂up = scatter_update(∂∂up, fixed_bddof, a0[fixed_bddof])
        end

        ∂u += 0.5 * Δt * ∂∂up

        i+1, write(d_arr, i+1, u), write(v_arr, i+1, ∂u), write(a_arr, i+1, ∂∂up)
    end

    arr_d = TensorArray(NT+1); arr_d = write(arr_d, 1, d0)
    arr_v = TensorArray(NT+1); arr_v = write(arr_v, 1, v0)
    arr_a = TensorArray(NT+1); arr_a = write(arr_a, 1, a0)
    i = constant(1, dtype=Int32)
    tas = [arr_d, arr_v, arr_a]
    _, d, v, a = while_loop(condition, body, [i, tas...])
    d, v, a = stack(d), stack(v), stack(a)
    sp = (NT+1, 2domain.nnodes)
    set_shape(d, sp), set_shape(v, sp), set_shape(a, sp)
end

@doc raw"""
    ExplicitSolver(globdat::GlobalData, domain::Domain,
        d0::Union{Array{Float64, 1}, PyObject},
        v0::Union{Array{Float64, 1}, PyObject},
        a0::Union{Array{Float64, 1}, PyObject},
        σ0::Union{Array{Float64, 1}, PyObject},
        ε0::Union{Array{Float64, 1}, PyObject},
        Δt::Float64, NT::Int64,
        nn::Function,
        Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
        ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
        abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; strain_type::String = "small")

Similar to [`ExplicitSolver`](@ref); however, the constitutive relation has the form

$$\sigma^{n+1} = \mathrm{nn}(\epsilon^{n+1}, \epsilon^n, \sigma^n)$$

Here the strain and stress are $n \times 3$ tensors. $n$ is the total number of Gaussian points and can be
obtained via `getNGauss(domain)`.
"""
function ExplicitSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject},
    v0::Union{Array{Float64, 1}, PyObject},
    a0::Union{Array{Float64, 1}, PyObject},
    σ0::Union{Array{Float64, 2}, PyObject},
    ε0::Union{Array{Float64, 2}, PyObject},
    Δt::Float64, NT::Int64,
    nn::Function,
    Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; strain_type::String = "small")
    assembleMassMatrix!(globdat, domain)
    if !(strain_type in ["small", "finite"])
        error("Only small strain or finite strain are supported. Unknown strain type $strain_type")
    end

    init_nnfem(domain)
    M = factorize(constant(globdat.M))
    bddof = findall(domain.EBC[:] .== -2)
    fixed_bddof = findall(domain.EBC[:] .== -1)

    Fext, ubd, abd = convert_to_tensor([Fext, ubd, abd], [Float64, Float64, Float64])

    function condition(i, tas...)
        i<=NT
    end
    function body(i, tas...)
        d_arr, v_arr, a_arr, σ_arr, ε_arr = tas
        u, ∂u, ∂∂u = read(d_arr, i), read(v_arr, i), read(a_arr, i)
        σc, εc = read(σ_arr, i), read(ε_arr, i)

        u +=  Δt*∂u + 0.5*Δt*Δt*∂∂u
        ∂u += 0.5*Δt * ∂∂u

        if !ismissing(abd)
            u = scatter_update(u, bddof, ubd[i])
        end
        if length(sum(fixed_bddof))>0
            u = scatter_update(u, fixed_bddof, d0[fixed_bddof])
        end

        if strain_type=="small"
            ε = s_eval_strain_on_gauss_points(u, domain)
        else
            ε = f_eval_strain_on_gauss_points(u, domain)
        end

        σ = nn(ε, εc, σc)
        if strain_type=="small"
            fint  = s_compute_internal_force_term(σ, domain)
        else
            fint  = f_compute_internal_force_term(σ, u, domain)
        end

        if ismissing(Fext)
            fext = zeros(length(fint))
        else
            fext = Fext[i]
        end
        ∂∂up = vector(findall(domain.dof_to_eq), M\(fext - fint), domain.nnodes*2)

        if !ismissing(abd)
            ∂∂up = scatter_update(∂∂up, bddof, abd[i])
        end
        if length(sum(fixed_bddof))>0
            ∂∂up = scatter_update(∂∂up, fixed_bddof, a0[fixed_bddof])
        end

        ∂u += 0.5 * Δt * ∂∂up

        i+1, write(d_arr, i+1, u), write(v_arr, i+1, ∂u), write(a_arr, i+1, ∂∂up), write(σ_arr, i+1, σ), write(ε_arr, i+1, ε)
    end

    arr_d = TensorArray(NT+1); arr_d = write(arr_d, 1, d0)
    arr_v = TensorArray(NT+1); arr_v = write(arr_v, 1, v0)
    arr_a = TensorArray(NT+1); arr_a = write(arr_a, 1, a0)
    arr_σ = TensorArray(NT+1); arr_σ = write(arr_σ, 1, σ0)
    arr_ε = TensorArray(NT+1); arr_ε = write(arr_ε, 1, ε0)
    i = constant(1, dtype=Int32)
    tas = [arr_d, arr_v, arr_a, arr_σ, arr_ε]
    _, d, v, a, σ, ε = while_loop(condition, body, [i, tas...])
    d, v, a, σ, ε = stack(d), stack(v), stack(a), stack(σ), stack(ε)
    sp = (NT+1, 2domain.nnodes)
    set_shape(d, sp), set_shape(v, sp), set_shape(a, sp), set_shape(σ, (NT+1, getNGauss(domain),3)), set_shape(ε, (NT+1, getNGauss(domain),3))
end

function ExplicitStaticSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject},
    v0::Union{Array{Float64, 1}, PyObject},
    Δt::Float64, NT::Int64,
    nn::Function,
    Fext::Union{Array{Float64, 1}, PyObject, Missing}=missing; strain_type::String = "small")
    assembleMassMatrix!(globdat, domain)
    @assert length(findall(domain.EBC[:] .== -2))==0
    @assert length(findall(domain.FBC[:] .== -2))==0
    a0 = tf.zeros_like(d0)
    Fext = repeat(constant(reshape(Fext, 1, :)), NT, 1)
    load_vec = reshape(Array((1:NT)/NT), NT, 1)
    Fext = Fext .* load_vec
    d, _, _ = ExplicitSolver(globdat, domain, d0, v0, a0, Δt, NT, nn, Fext, missing, missing; strain_type = strain_type)
    return d
end

@doc raw"""
    LinearStaticSolver(globdat::GlobalData, domain::Domain,
        d0::Union{Array{Float64, 1}, PyObject},
        Hs::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
        Fext::Union{Array{Float64, 1}, PyObject, Missing}=missing)

An AD-capable static linear elasticity solver.

- `d0` has length `2domain.nnodes` (full DOFs)

- `Hs` is either a $3\times 3$ matrix or $N\times 3\times 3$ tensor, representing elasticity matrix.
"""
function LinearStaticSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject},
    Hs::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
    Fext::Union{Array{Float64, 1}, PyObject, Missing}=missing)
    @assert length(findall(domain.EBC[:] .== -2))==0
    @assert length(findall(domain.FBC[:] .== -2))==0
    init_nnfem(domain)
    Hs = constant(Hs)
    if length(size(Hs))==2
        Hs = repeat(reshape(Hs, (-1,)), getNGauss(domain))
        Hs = reshape(Hs, (getNGauss(domain), 3, 3) )
    end
    Fext = coalesce(Fext, zeros(domain.neqs))
    d0, Fext = convert_to_tensor([d0, Fext], [Float64, Float64])
    ε = s_eval_strain_on_gauss_points(d0, domain)
    σ = batch_matmul(Hs, ε)
    fint = s_compute_internal_force_term(σ, domain)
    stiff = s_compute_stiffness_matrix(Hs, domain)
    u0 = stiff\(Fext - fint)
    nbdnode = findall(domain.EBC[:].!=-1)
    d = scatter_update(d0, nbdnode, u0)
    return d
end

function ImplicitStaticSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject},
    Δt::Float64, NT::Int64,
    Hs::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
    Fext::Union{Array{Float64, 1}, PyObject, Missing}=missing;
    strain_type::String = "small", method::String = "newton")
    if !(strain_type in ["small", "finite"])
        error("Only small strain or finite strain are supported. Unknown strain type $strain_type")
    end
    assembleMassMatrix!(globdat, domain)
    @assert length(findall(domain.EBC[:] .== -2))==0
    @assert length(findall(domain.FBC[:] .== -2))==0

    Hs = constant(Hs)
    if length(size(Hs))==2
        Hs = repeat(reshape(Hs, (-1,)), getNGauss(domain))
        Hs = reshape(Hs, (getNGauss(domain), 3, 3) )
    end
    if method=="incremental"
        a0 = tf.zeros_like(d0)
        v0 = tf.zeros_like(d0)
        Fext = repeat(constant(reshape(Fext, 1, :)), NT, 1)
        d, _, _ = GeneralizedAlphaSolver(globdat, domain, d0, v0, a0, Δt, NT, Hs, Fext, missing, missing)
    elseif method=="newton"
        if !(strain_type in ["small"])
            error("Only small strain or finite strain are supported. Unknown strain type $strain_type")
        end
        d0, Fext = convert_to_tensor([d0, Fext], [Float64, Float64])
        stiff = s_compute_stiffness_matrix(Hs, domain)
        u0 = stiff\Fext
        nbdnode = findall(domain.EBC[:].!=-1)
        d = scatter_add(d0, nbdnode, u0)
    end
    return d
end

@doc raw"""
    ImplicitStaticSolver(globdat::GlobalData, domain::Domain,
        d0::Union{Array{Float64, 1}, PyObject},
        nn::Function, θ::Union{Array{Float64, 1}, PyObject},
        Fext::Union{Array{Float64, 1}, PyObject, Missing}=missing)

Solves the static problem

$$K(u) = F$$

using Newton's method. Users provide `nn`, which is a function that outputs stress and stress sensitivity given the strain tensor.

$$\sigma, \frac{\partial \sigma}{\partial \epsilon} = \mathrm{nn}(\epsilon, \theta)$$

- `d0`: a vector of length `2domain.nnodes`, the fixed Dirichlet DOF should be populated with boundary values.

- `Fext`: external force, a vector of length `domain.neqs`
"""
function ImplicitStaticSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject},
    nn::Function, θ::Union{Array{Float64, 1}, PyObject},
    Fext::Union{Array{Float64, 1}, PyObject, Missing}=missing)
    init_nnfem(domain)
    assembleMassMatrix!(globdat, domain)
    @assert length(findall(domain.EBC[:] .== -2))==0
    @assert length(findall(domain.FBC[:] .== -2))==0
    nbdnode = findall(domain.EBC[:].!=-1)
    K0 = spdiag(domain.nnodes*2)
    function residual_and_jac(θ, u)
        ε = s_eval_strain_on_gauss_points(u, domain)
        σ, Hs = nn(ε, θ)
        stiff = s_compute_stiffness_matrix(Hs, domain)
        fint = s_compute_internal_force_term(σ, domain)
        # @info fint, σ, Hs, stiff
        δ = scatter_update(d0, nbdnode, fint - Fext)
        K = scatter_update(K0, nbdnode, nbdnode, stiff)
        # @info δ, K
        δ, K
    end
    d0, Fext, θ = convert_to_tensor([d0, Fext, θ], [Float64, Float64, Float64])
    u0 = newton_raphson_with_grad(residual_and_jac, d0, θ)
    return u0
end


text:
module GameObjectManager

export IGameObject
export GameObject

using CoreExtended
using JLGEngine
using ..Management

""" TODO """
IComponent = JLGEngine.IComponent

""" TODO """
type GameObject <: IComponent
	id::Symbol
	enabled::Bool
	components::SortedDict{Symbol, IComponent}

	#transform::Transform

	GameObject(id=:_) = new(id,true,SortedDict{Symbol, IComponent}(Forward))
end

""" TODO """
setComponent(this::GameObject, c::IComponent) = (this.components[k] = c; c.gameObject=this)

""" TODO """
getComponent(this::GameObject, T::DataType) = getComponents(this, T).first

""" TODO """
function getComponents(this::GameObject, T::DataType)
	components = SortedDict{Symbol, IComponent}(Forward)
	for (k,v) in this.components
		if typeof(v) == T
			components[k] = v
		end
	end
	components
end

#addComponent(this::GameObject, c::Transform) = (setComponent(this,c); this.transform=c)

# preset for manager: includes code here
#presetManager(GameObject)
m=Management.presetManager(GameObject)

end #GameObjectManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
# # Dataset Previews

# First, initialize by adding tools and declare floating type
using GriddingMachine
using PkgUtility
using Plots
using Plots.PlotMeasures

ENV["GKSwstype"]="100";
FT = Float32;
#------------------------------------------------------------------------------

# Then, define a function to plot the dataset
function preview_data(ds::GriddedDataset{FT}, ind::Int)
    ## preview data
    return heatmap(view(ds.data,:,:,ind)',
                   origin="lower",
                   aspect_ratio=1,
                   xticks=[],
                   yticks=[],
                   c=:viridis,
                   size=(700,300),
                   framestyle=:none)
end

function preview_data(ds::GriddedDataset{FT}, ind::Int, clim::Tuple)
    ## preview data
    return heatmap(view(ds.data,:,:,ind)',
                   origin="lower",
                   aspect_ratio=1,
                   xticks=[],
                   yticks=[],
                   c=:viridis,
                   clim=clim,
                   size=(700,300),
                   framestyle=:none)
end
#------------------------------------------------------------------------------

# ## Leaf level datasets
# ### Leaf nitrogen content
LNC_LUT = load_LUT(LeafNitrogenButler{FT}(), 1);
preview_data(LNC_LUT, 1)
#------------------------------------------------------------------------------

LNC_LUT = load_LUT(LeafNitrogenBoonman{FT}(), 1);
preview_data(LNC_LUT, 1)
#------------------------------------------------------------------------------

# ### Leaf phosphorus content
LPC_LUT = load_LUT(LeafPhosphorus{FT}(), 1);
preview_data(LPC_LUT, 1)
#------------------------------------------------------------------------------

# ### Specific leaf area
SLA_LUT = load_LUT(LeafSLAButler{FT}(), 1);
preview_data(SLA_LUT, 1)
#------------------------------------------------------------------------------

SLA_LUT = load_LUT(LeafSLABoonman{FT}(), 1);
preview_data(SLA_LUT, 1)
#------------------------------------------------------------------------------

# ### Vcmax
VCM_LUT = load_LUT(VcmaxOptimalCiCa{FT}(), 1);
preview_data(VCM_LUT, 1)
#------------------------------------------------------------------------------

# ## Stand level datasets
# ### Canopy height
CHT_LUT = load_LUT(CanopyHeightGLAS{FT}(), 1);
preview_data(CHT_LUT, 1)
#------------------------------------------------------------------------------

CHT_LUT = load_LUT(CanopyHeightBoonman{FT}(), 1);
preview_data(CHT_LUT, 1)
#------------------------------------------------------------------------------

# ### Clumping index
## global clumping index
CLI_LUT = load_LUT(ClumpingIndexMODIS{FT}(), "12X", "1Y", 1);
preview_data(CLI_LUT, 1, (0.4,1))
#------------------------------------------------------------------------------

## global clumping index per PFT
CLI_LUT = load_LUT(ClumpingIndexPFT{FT}(), 1);
anim = @animate for i ∈ 1:size(CLI_LUT.data,3)
    preview_data(CLI_LUT, i, (0.4,1));
end
gif(anim, fps=1)
#------------------------------------------------------------------------------

# ### Gross primary productivity
## GPP MPI
GPP_LUT = load_LUT(GPPMPIv006{FT}(), 2005, "2X", "8D", 1);
anim = @animate for i ∈ 1:46
    preview_data(GPP_LUT, i, (0,10));
end
gif(anim, fps=5)
#------------------------------------------------------------------------------

## GPP VPM
GPP_LUT = load_LUT(GPPVPMv20{FT}(), 2005, "5X", "8D", 1);
anim = @animate for i ∈ 1:46
    preview_data(GPP_LUT, i, (0,10));
end
gif(anim, fps=5)
#------------------------------------------------------------------------------

# ### Leaf area index
## Annual data
LAI_LUT = load_LUT(LAIMODISv006{FT}(), 2005, "2X", "8D", 1);
anim = @animate for i ∈ 1:46
    preview_data(LAI_LUT, i, (0,6));
end
gif(anim, fps=5)
#------------------------------------------------------------------------------

## monthly mean of multiple years
LAI_LUT = load_LUT(LAIMonthlyMean{FT}(), 1);
anim = @animate for i ∈ 1:size(LAI_LUT.data,3)
    preview_data(LAI_LUT, i, (0,6));
end
gif(anim, fps=1)
#------------------------------------------------------------------------------

# ### Normalized difference vegetation index
NDV_LUT = load_LUT(NDVIAvhrr{FT}(), 2018, "20X", "1M", 1);
anim = @animate for i ∈ 1:size(NDV_LUT.data,3)
    preview_data(NDV_LUT, i, (0,1));
end
gif(anim, fps=1)
#------------------------------------------------------------------------------

# ### Near infrared reflectance of vegetation
NIV_LUT = load_LUT(NIRvAvhrr{FT}(), 2018, "20X", "1M", 1);
anim = @animate for i ∈ 1:size(NIV_LUT.data,3)
    preview_data(NIV_LUT, i, (0,1));
end
gif(anim, fps=1)
#------------------------------------------------------------------------------

# ### Near infrared reflectance of vegetation with offset
NIO_LUT = load_LUT(NIRoAvhrr{FT}(), 2018, "20X", "1M", 1);
anim = @animate for i ∈ 1:size(NIO_LUT.data,3)
    preview_data(NIO_LUT, i, (0,1));
end
gif(anim, fps=1)
#------------------------------------------------------------------------------

# ### Net primary productivity
NPP_LUT = load_LUT(NPPModis{FT}(), 1);
NPP_LUT.data .*= 1e9;
preview_data(NPP_LUT, 1)
#------------------------------------------------------------------------------

# ### PFT maps
PFT_LUT = load_LUT(PFTPercentCLM{FT}(), 1);
anim = @animate for i ∈ 1:17
    preview_data(PFT_LUT, i);
end
gif(anim, fps=1)
#------------------------------------------------------------------------------

# ### Soil color class
SCC_LUT = load_LUT(SoilColor{FT}(), 1);
preview_data(SCC_LUT, 1)
#------------------------------------------------------------------------------

# ### Soil van Genuchten alpha
VGA_LUT = load_LUT(VGMAlphaJules{FT}(), "12X", "1Y", 1);
anim = @animate for i ∈ 1:4
    preview_data(VGA_LUT, i);
end
gif(anim, fps=1)
#------------------------------------------------------------------------------

# ### Soil van Genuchten log(n)
VGN_LUT = load_LUT(VGMLogNJules{FT}(), "12X", "1Y", 1);
anim = @animate for i ∈ 1:4
    preview_data(VGN_LUT, i);
end
gif(anim, fps=1)
#------------------------------------------------------------------------------

# ### Soil van Genuchten residual SWC
VGT_LUT = load_LUT(VGMThetaRJules{FT}(), "12X", "1Y", 1);
anim = @animate for i ∈ 1:4
    preview_data(VGT_LUT, i);
end
gif(anim, fps=1)
#------------------------------------------------------------------------------

# ### Soil van Genuchten saturated SWC
VGT_LUT = load_LUT(VGMThetaSJules{FT}(), "12X", "1Y", 1);
anim = @animate for i ∈ 1:4
    preview_data(VGT_LUT, i);
end
gif(anim, fps=1)
#------------------------------------------------------------------------------

# ### Sun induced fluorescence
SIF_LUT = load_LUT(SIFTropomi740{FT}(), 2018, "1X", "1M", 1);
anim = @animate for i ∈ 1:12
    preview_data(SIF_LUT, i, (0,3.5));
end
gif(anim, fps=3)
#------------------------------------------------------------------------------

# ### Surface area
SUA_LUT = load_LUT(SurfaceAreaCLM{FT}(), 1);
preview_data(SUA_LUT, 1)
#------------------------------------------------------------------------------

# ### Tree density
TDT_LUT = load_LUT(TreeDensity{FT}(), "12X", "1Y", 1);
preview_data(TDT_LUT, 1, (0, 150000))
#------------------------------------------------------------------------------

# ### Wood density
TDT_LUT = load_LUT(WoodDensity{FT}(), 1);
preview_data(TDT_LUT, 1)
#------------------------------------------------------------------------------

# ## Land surface
# ### Land mask
LMK_LUT = load_LUT(LandMaskERA5{FT}(), 1);
preview_data(LMK_LUT, 1)
#------------------------------------------------------------------------------


text:
"""
Generate figure 2 which the difference between anti-lensing and inverse lensing
===================================================== """

"""
Load packages. Set grid etc.
----------------------------------------------------- """
srcpath  = "/Users/ethananderes/Dropbox/BayesLense/src/"
savepath = "/Users/ethananderes/Dropbox/BayesLense/paper/"

include(srcpath*"Interp.jl")
include(srcpath*"cmb.jl")
include(srcpath*"fft.jl")

using PyPlot, PyCall, Interp
@pyimport scipy.interpolate as scii

seed = 1
srand(seed)

parlr = setpar(2.0, 2^9,  0, 0, 3000, 600, srcpath);
parhr = setpar(0.5, 2^11, 0, 0, 3000, 600, srcpath);

""" Simulate the lensing potential, define functions
-----------------------------------------------"""
phik   = fft2(randn(size(parhr.grd.x))./parhr.grd.deltx, parhr) .* sqrt(parhr.cPP)
phix   = ifft2r(phik, parhr.grd.deltk)
phidx1 = ifft2r(im .* parhr.grd.k1 .* phik, parhr)
phidx2 = ifft2r(im .* parhr.grd.k2 .* phik, parhr)

function griddata(x::Matrix, y::Matrix, z::Matrix, xi::Matrix,yi::Matrix)
	xpd, ypd = Interp.perodic_padd_xy(x, y, 0.1)
	zpd = Interp.perodic_padd_z(z, 0.1)
	points = [xpd[:] ypd[:]]
	# grid = (vec(xi[1,:]), vec(yi[:,1]))
	grid = (xi, yi)
	zi = scii.griddata(points, zpd[:], grid, method = "cubic")
end

lensex = parhr.grd.x + phidx1
lensey = parhr.grd.y + phidx2
antix_dsp = griddata(
	lensex,
	lensey,
	-phidx1,
	parlr.grd.x,
	parlr.grd.y
)
antiy_dsp = griddata(
	lensex,
	lensey,
	-phidx2,
	parlr.grd.x,
	parlr.grd.y
)

function helmholtz(ax::Matrix, bx::Matrix, par::SpectrumGrids)
	# (ax,bx) is the vector field defined in pixel space
	ak, bk = fft2(ax, par), fft2(bx, par)
	k₁, k₂ = par.grd.k1, par.grd.k2
	adk = (k₁./k₂) .* ((k₁./k₂) .* ak + bk) ./ ((k₁./k₂).^2 .+ 1) # curl free pat
	ack = ak - adk
	divk = adk ./ ( im * k₁)
	crlk = ack ./ (-im * k₂)
	divk[par.grd.k1 .== 0.0] = bk[par.grd.k1 .== 0.0] ./ ( im * k₂[par.grd.k1 .== 0.0])
	divk[par.grd.k2 .== 0.0] = ak[par.grd.k2 .== 0.0] ./ ( im * k₁[par.grd.k2 .== 0.0])
	crlk[par.grd.k1 .== 0.0] = ak[par.grd.k1 .== 0.0] ./ (-im * k₂[par.grd.k1 .== 0.0])
	crlk[par.grd.k2 .== 0.0] = bk[par.grd.k2 .== 0.0] ./ ( im * k₁[par.grd.k2 .== 0.0])
	divk[par.grd.r .<= 0.0] = 0.0
	crlk[par.grd.r .<= 0.0] = 0.0
	divx, crlx = ifft2r(divk, par), ifft2r(crlk, par)
	divx, divk, crlx, crlk
end
divx, divk, crlx, crlk = helmholtz(antix_dsp, antiy_dsp, parlr)

"""
figure2a.pdf
-----------------------------------"""
imshow(
  -phix[1:4:end, 1:4:end],
  interpolation = "nearest",
  origin="lower",
  extent=(180/pi)*[minimum(parlr.grd.x), maximum(parlr.grd.x),minimum(parlr.grd.x), maximum(parlr.grd.x)]
)
title(L"$-\phi$")
xlabel("degrees")
ylabel("degrees")
colorbar(format="%.e")
savefig(joinpath(savepath, "figure2a.pdf"), dpi=300, bbox_inches="tight", pad_inches=0)

"""
figure2b.pdf
-----------------------------------"""
imshow(
  divx,
  interpolation = "nearest",
  origin="lower",
  extent=(180/pi)*[minimum(parlr.grd.x), maximum(parlr.grd.x),minimum(parlr.grd.x), maximum(parlr.grd.x)]
)
title(L"$-\phi^{inv}$")
xlabel("degrees")
ylabel("degrees")
colorbar(format="%.e")
savefig(joinpath(savepath, "figure2b.pdf"), dpi=300, bbox_inches="tight", pad_inches=0)

"""
figure2c.pdf
-----------------------------------"""
imshow(
  divx + phix[1:4:end, 1:4:end],
  interpolation = "nearest",
  origin="lower",
  extent=(180/pi)*[minimum(parlr.grd.x), maximum(parlr.grd.x),minimum(parlr.grd.x), maximum(parlr.grd.x)]
)
title(L"$\phi-\phi^{inv}$")
xlabel("degrees")
ylabel("degrees")
colorbar(format="%.e")
savefig(joinpath(savepath, "figure2c.pdf"), dpi=300, bbox_inches="tight", pad_inches=0)

"""
figure2d.pdf
-----------------------------------"""
imshow(
  crlx,
  vmin = -1.2e-8, vmax = 1.2e-8,
  interpolation = "nearest",
  origin="lower",
  extent=(180/pi)*[minimum(parlr.grd.x), maximum(parlr.grd.x),minimum(parlr.grd.x), maximum(parlr.grd.x)]
)
title(L"$-\psi^{inv}$")
xlabel("degrees")
ylabel("degrees")
colorbar(format="%.e")
savefig(joinpath(savepath, "figure2d.pdf"), dpi=300, bbox_inches="tight", pad_inches=0)

max_stars_repo_name: EthanAnderes/BayesianCmbLensing

max_issues_repo_name: EthanAnderes/BayesianCmbLensing

max_forks_repo_name: EthanAnderes/BayesianCmbLensing


text:
pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

using Plots
using Printf
using Oceananigans
using Oceananigans.Units
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity

grid = RectilinearGrid(size=16, z=(-64, 0), topology=(Flat, Flat, Bounded))

closure = CATKEVerticalDiffusivity()

Qᵇ = 1e-8
Qᵘ = 0.0
Qᵛ = 0.0

u★ = (Qᵘ^2 + Qᵛ^2)^(1/4)
w★ = Qᵇ * grid.Δzᵃᵃᶜ

u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))
v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ))
b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ))

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    tracers = (:b, :e),
                                    buoyancy = BuoyancyTracer(),
                                    boundary_conditions = (; b=b_bcs),
                                    closure = closure)

N² = 1e-5
bᵢ(x, y, z) = N² * z
set!(model, b = bᵢ)

z = znodes(model.tracers.b)

b = view(interior(model.tracers.b), 1, 1, :)
e = view(interior(model.tracers.e), 1, 1, :)
Kc = view(interior(model.diffusivity_fields.Kᶜ), 1, 1, :)
Ke = view(interior(model.diffusivity_fields.Kᵉ), 1, 1, :)

b_plot = plot(b, z, linewidth = 2, label = "t = 0", xlabel = "Buoyancy", ylabel = "z", legend=:bottomright)
e_plot = plot(e, z, linewidth = 2, label = "t = 0", xlabel = "TKE", ylabel = "z", legend=:bottomright)

simulation = Simulation(model, Δt = 20.0, stop_time = 48hours)

run!(simulation)

plot!(b_plot, b, z, linewidth = 2, label = @sprintf("t = %s", prettytime(model.clock.time)))
plot!(e_plot, e, z, linewidth = 2, label = @sprintf("t = %s", prettytime(model.clock.time)))
K_plot = plot(Kc, z, linewidth = 2, linestyle=:dash, label = @sprintf("Kᶜ, t = %s", prettytime(model.clock.time)), legend=:bottomright, xlabel="Diffusivities")
plot!(K_plot, Ke, z, linewidth = 3, alpha=0.6, label = @sprintf("Kᵉ, t = %s", prettytime(model.clock.time)))

eb_plot = plot(b_plot, e_plot, K_plot, layout=(1, 3), size=(1200, 600))

display(eb_plot)


text:
using ArgParse

function make_arg_parser()
  s = ArgParseSettings()
  @add_arg_table s begin
    "--data"
      default = "demo"
    "--n-epoch"
        arg_type = Int
        default = 100
    "--regu-coef"
      arg_type = Float64
      default = 0.01
  end
  return s
end


text:
export get_discrete_Grad

"""
2D version
 input:
# Arguments
 - n1 : number of grid points in the first dimension
 - n2 : number of grid points in the second dimension
 - h1 : distance between grid points in the first dimension
 - h2 : distance between grid points in the second dimension
 - TD_type : type of derivative operator as a string

 output :
 - TD_OP : transform domain operator as a sparse matrix
"""
function get_discrete_Grad(n1,n2,h1::TF,h2::TF,TD_type::String) where {TF<:Real}
    TI=Int64

    #define difference matrix D acting on vectorized model using Kronecker products
    Ix = SparseMatrixCSC{TF}(LinearAlgebra.I,n1,n1) #x
    Iz = SparseMatrixCSC{TF}(LinearAlgebra.I,n2,n2) #z
    Dx = spdiagm(0 => ones(TF,n1-1)*-1, 1 => ones(TF,n1-1)*1); Dx = Dx[1:end-1,:] ./ h1
    Dz = spdiagm(0 => ones(TF,n2-1)*-1, 1 => ones(TF,n2-1)*1); Dz = Dz[1:end-1,:] ./ h2

    if TD_type=="D_z"
      D_OP = kron(Dz,Ix) #D2z
    elseif TD_type=="D_x"
      D_OP = kron(Iz,Dx) #D2x
    elseif  TD_type=="TV" || TD_type=="D2D"
      D2z  = kron(Dz,Ix)
      D2x  = kron(Iz,Dx)
      D_OP = vcat(D2z,D2x) #D2D
    end

    return D_OP
end

"""
  3D version
  input: n1 : number of grid points in the first dimension
         n2 : number of grid points in the second dimension
         n3 : number of grid points in the third dimension
         h1 : distance between grid points in the first dimension
         h2 : distance between grid points in the second dimension
         h3 : distance between grid points in the third dimension
         TD_type : type of derivative operator as a string
   output : TD_OP : transform domain operator as a sparse matrix
"""
function get_discrete_Grad(n1,n2,n3,h1::TF,h2::TF,h3::TF,TD_type::String) where {TF<:Real}
    TI = Int64

    #define difference matrix D acting on vectorized model using Kronecker products
    Ix = SparseMatrixCSC{TF}(LinearAlgebra.I,n1,n1) #x
    Iy = SparseMatrixCSC{TF}(LinearAlgebra.I,n2,n2) #x
    Iz = SparseMatrixCSC{TF}(LinearAlgebra.I,n3,n3) #z
    Dx = spdiagm(0 => ones(TF,n1-1)*-1, 1 => ones(TF,n1-1)*1); Dx = Dx[1:end-1,:] ./ h1
    Dy = spdiagm(0 => ones(TF,n2-1)*-1, 1 => ones(TF,n2-1)*1); Dy = Dy[1:end-1,:] ./ h2
    Dz = spdiagm(0 => ones(TF,n3-1)*-1, 1 => ones(TF,n3-1)*1); Dz = Dz[1:end-1,:] ./ h3

    if TD_type=="D_z"
      D_OP = kron(Dz,Iy,Ix)
    elseif TD_type=="D_y"
      D_OP = kron(Iz,Dy,Ix)
    elseif TD_type=="D_x"
      D_OP = kron(Iz,Iy,Dx)
    elseif TD_type=="TV" || TD_type=="D3D"
      D3z = kron(Dz,Iy,Ix)
      D3y = kron(Iz,Dy,Ix)
      D3x = kron(Iz,Iy,Dx)
      D_OP = vcat(D3z,D3y,D3x)
    end

    return D_OP
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
# a variety of rankings
#
# Please refer to http://en.wikipedia.org/wiki/Ranking#Strategies_for_assigning_rankings
# to see the definitions of a variety of ranking strategies
#
# The implementations here follow this wikipedia page.
#

function _check_randparams(rks, x, p)
    n = length(rks)
    length(x) == length(p) == n || raise_dimerror()
    return n
end

# Ordinal ranking ("1234 ranking") -- use the literal order resulted from sort
function ordinalrank!(rks::RealArray, x::RealArray, p::IntegerArray)
    n = _check_randparams(rks, x, p)

    if n > 0
        i = 1
        while i <= n
            rks[p[i]] = i
            i += 1
        end
    end

    return rks
end

"""
    ordinalrank(x)

Return the ordinal ranking ("1234" ranking) of a real-valued array.
All items in `x` are given distinct, successive ranks based on their
position in `sort(x)`.
"""
ordinalrank(x::RealArray) = ordinalrank!(Array(Int, size(x)), x, sortperm(x))

# Competition ranking ("1224" ranking) -- resolve tied ranks using min
function competerank!(rks::RealArray, x::RealArray, p::IntegerArray)
    n = _check_randparams(rks, x, p)

    if n > 0
        p1 = p[1]
        v = x[p1]
        rks[p1] = k = 1

        i = 2
        while i <= n
            pi = p[i]
            xi = x[pi]
            if xi == v
                rks[pi] = k
            else
                rks[pi] = k = i
                v = xi
            end
            i += 1
        end
    end

    return rks
end

"""
    competerank(x)

Return the standard competition ranking ("1224" ranking) of a real-valued
array. Items that compare equal are given the same rank, then a gap is left
in the rankings the size of the number of tied items - 1.
"""
competerank(x::RealArray) = competerank!(Array(Int, size(x)), x, sortperm(x))

# Dense ranking ("1223" ranking) -- resolve tied ranks using min
function denserank!(rks::RealArray, x::RealArray, p::IntegerArray)
    n = _check_randparams(rks, x, p)

    if n > 0
        p1 = p[1]
        v = x[p1]
        rks[p1] = k = 1

        i = 2
        while i <= n
            pi = p[i]
            xi = x[pi]
            if xi == v
                rks[pi] = k
            else
                rks[pi] = (k += 1)
                v = xi
            end
            i += 1
        end
    end

    return rks
end

"""
    denserank(x)

Return the dense ranking ("1223" ranking) of a real-valued array. Items that
compare equal receive the same ranking, and the next subsequent rank is
assigned with no gap.
"""
denserank(x::RealArray) = denserank!(Array(Int, size(x)), x, sortperm(x))

# Tied ranking ("1 2.5 2.5 4" ranking) -- resolve tied ranks using average
function tiedrank!(rks::RealArray, x::RealArray, p::IntegerArray)
    n = _check_randparams(rks, x, p)

    if n > 0
        v = x[p[1]]

        s = 1  # starting index of current range
        e = 2  # pass-by-end index of current range
        while e <= n
            cx = x[p[e]]
            if cx != v
                # fill average rank to s : e-1
                ar = (s + e - 1) / 2
                for i = s : e-1
                    rks[p[i]] = ar
                end
                # switch to next range
                s = e
                v = cx
            end
            e += 1
        end

        # the last range (e == n+1)
        ar = (s + n) / 2
        for i = s : n
            rks[p[i]] = ar
        end
    end

    return rks
end

# order (aka. rank), resolving ties using the mean rank
"""
    tiedrank(x)

Return the tied ranking, also called fractional or "1 2.5 2.5 4" ranking,
of a real-valued array. Items that compare equal receive the mean of the
rankings they would have been assigned under ordinal ranking.
"""
tiedrank(x::RealArray) = tiedrank!(Array(Float64, size(x)), x, sortperm(x))


text:
# Global component registry: @defcomp stores component definitions here
global const _compdefs = Dict{ComponentId, ComponentDef}()

compdefs() = collect(values(_compdefs))

compdef(comp_id::ComponentId) = _compdefs[comp_id]

function compdef(comp_name::Symbol)
    matches = collect(Iterators.filter(obj -> name(obj) == comp_name, values(_compdefs)))
    count = length(matches)

    if count == 1
        return matches[1]
    elseif count == 0
        error("Component $comp_name was not found in the global registry")
    else
        error("Multiple components named $comp_name were found in the global registry")
    end
end

compdefs(md::ModelDef) = values(md.comp_defs)

compkeys(md::ModelDef) = keys(md.comp_defs)

hascomp(md::ModelDef, comp_name::Symbol) = haskey(md.comp_defs, comp_name)

compdef(md::ModelDef, comp_name::Symbol) = md.comp_defs[comp_name]

function reset_compdefs(reload_builtins=true)
    empty!(_compdefs)

    if reload_builtins
        compdir = joinpath(dirname(@__FILE__), "..", "components")
        load_comps(compdir)
    end
end

first_period(comp_def::ComponentDef) = comp_def.first
first_period(md::ModelDef, comp_def::ComponentDef) = first_period(comp_def) === nothing ? time_labels(md)[1] : first_period(comp_def)

last_period(comp_def::ComponentDef) = comp_def.last
last_period(md::ModelDef, comp_def::ComponentDef) = last_period(comp_def) === nothing ? time_labels(md)[end] : last_period(comp_def)

# Return the module object for the component was defined in
compmodule(comp_id::ComponentId) = comp_id.module_name

compname(comp_id::ComponentId) = comp_id.comp_name

function Base.show(io::IO, comp_id::ComponentId)
    print(io, "$(comp_id.module_name).$(comp_id.comp_name)")
end

"""
    name(def::NamedDef) = def.name

Return the name of `def`.  Possible `NamedDef`s include `DatumDef`, and `ComponentDef`.
"""
name(def::NamedDef) = def.name

number_type(md::ModelDef) = md.number_type

numcomponents(md::ModelDef) = length(md.comp_defs)

function dump_components()
    for comp in compdefs()
        println("\n$(name(comp))")
        for (tag, objs) in ((:Variables, variables(comp)), (:Parameters, parameters(comp)), (:Dimensions, dimensions(comp)))
            println("  $tag")
            for obj in objs
                println("    $(obj.name) = $obj")
            end
        end
    end
end

"""
    new_comp(comp_id::ComponentId, verbose::Bool=true)

Add an empty `ComponentDef` to the global component registry with the given
`comp_id`. The empty `ComponentDef` must be populated with calls to `addvariable`,
`addparameter`, etc.
"""
function new_comp(comp_id::ComponentId, verbose::Bool=true)
    if verbose
        if haskey(_compdefs, comp_id)
            @warn "Redefining component $comp_id"
        else
            @info "new component $comp_id"
        end
    end

    comp_def = ComponentDef(comp_id)
    _compdefs[comp_id] = comp_def
    return comp_def
end

"""
    delete!(m::ModelDef, component::Symbol

Delete a `component` by name from a model definition `m`.
"""
function Base.delete!(md::ModelDef, comp_name::Symbol)
    if ! haskey(md.comp_defs, comp_name)
        error("Cannot delete '$comp_name' from model; component does not exist.")
    end

    delete!(md.comp_defs, comp_name)

    ipc_filter = x -> x.src_comp_name != comp_name && x.dst_comp_name != comp_name
    filter!(ipc_filter, md.internal_param_conns)

    epc_filter = x -> x.comp_name != comp_name
    filter!(epc_filter, md.external_param_conns)
end

#
# Dimensions
#
"""
    add_dimension!(comp::ComponentDef, name)

Add a dimension name to a `ComponentDef`, with an optional Dimension definition.
The definition is included only for the case of anonymous dimensions, which are
indicated in `@defcomp` with an Int rather than a symbol name. In this case, the
Int (e.g., 4) is converted to a dimension of the same length, e.g., `Dimension(4)`,
and assigned a new name `Symbol(4)``
"""
function add_dimension!(comp::ComponentDef, name)
    # create the Dimension and store it in the ComponentDef until build time
    dim = (name isa Int) ? Dimension(name) : nothing
    comp.dimensions[Symbol(name)] = dim
end

add_dimension!(comp_id::ComponentId, name) = add_dimension!(compdef(comp_id), name)

dimensions(comp_def::ComponentDef) = keys(comp_def.dimensions)

dimensions(def::DatumDef) = def.dimensions

dimensions(comp_def::ComponentDef, datum_name::Symbol) = dimensions(datumdef(comp_def, datum_name))

dim_count(def::DatumDef) = length(def.dimensions)

datatype(def::DatumDef) = def.datatype

description(def::DatumDef) = def.description

unit(def::DatumDef) = def.unit

function first_and_step(md::ModelDef)
    keys::Vector{Int} = time_labels(md) # labels are the first times of the model runs
    return first_and_step(keys)
end

function first_and_step(values::Vector{Int})
     return values[1], (length(values) > 1 ? values[2] - values[1] : 1)
end

function time_labels(md::ModelDef)
    keys::Vector{Int} = dim_keys(md, :time)
    return keys
end

function check_parameter_dimensions(md::ModelDef, value::AbstractArray, dims::Vector, name::Symbol)
    for dim in dims
        if haskey(md, dim)
            if isa(value, NamedArray)
                labels = names(value, findnext(isequal(dim), dims, 1))
                dim_vals = dim_keys(md, dim)
                for i in 1:length(labels)
                    if labels[i] != dim_vals[i]
                        error("Labels for dimension $dim in parameter $name do not match model's index values")
                    end
                end
            end
        else
            error("Dimension $dim in parameter $name not found in model's dimensions")
        end
    end
end

function datum_size(md::ModelDef, comp_def::ComponentDef, datum_name::Symbol)
    dims = dimensions(comp_def, datum_name)
    if dims[1] == :time
        time_length = getspan(md, comp_def)[1]
        rest_dims = filter(x->x!=:time, dims)
        datum_size = (time_length, dim_counts(md, rest_dims)...,)
    else
        datum_size = (dim_counts(md, dims)...,)
    end
    return datum_size
end

dimensions(md::ModelDef) = md.dimensions
dimensions(md::ModelDef, dims::Vector{Symbol}) = [dimension(md, dim) for dim in dims]
dimension(md::ModelDef, name::Symbol) = md.dimensions[name]

dim_count_dict(md::ModelDef) = Dict([name => length(value) for (name, value) in dimensions(md)])
dim_counts(md::ModelDef, dims::Vector{Symbol}) = [length(dim) for dim in dimensions(md, dims)]
dim_count(md::ModelDef, name::Symbol) = length(dimension(md, name))

dim_key_dict(md::ModelDef) = Dict([name => collect(keys(dim)) for (name, dim) in dimensions(md)])
dim_keys(md::ModelDef, name::Symbol) = collect(keys(dimension(md, name)))

dim_values(md::ModelDef, name::Symbol) = collect(values(dimension(md, name)))
dim_value_dict(md::ModelDef) = Dict([name => collect(values(dim)) for (name, dim) in dimensions(md)])

Base.haskey(md::ModelDef, name::Symbol) = haskey(md.dimensions, name)

isuniform(md::ModelDef) = md.is_uniform

# Helper function invoked when the user resets the time dimension with set_dimension!
# This function calls set_run_period! on each component definition to reset the first and last values.
function reset_run_periods!(md, first, last)
    for comp_def in compdefs(md)
        changed = false
        first_per = first_period(comp_def)
        last_per  = last_period(comp_def)

        if first_per !== nothing && first_per < first
            @warn "Resetting $(comp_def.name) component's first timestep to $first"
            changed = true
        else
            first = first_per
        end

        if last_per !== nothing && last_per > last
            @warn "Resetting $(comp_def.name) component's last timestep to $last"
            changed = true
        else
            last = last_per
        end

        if changed
            set_run_period!(comp_def, first, last)
        end
    end
    nothing
end

"""
    set_dimension!(md::ModelDef, name::Symbol, keys::Union{Int, Vector, Tuple, Range})

Set the values of `md` dimension `name` to integers 1 through `count`, if `keys` is
an integer; or to the values in the vector or range if `keys` is either of those types.
"""
set_dimension!(md::ModelDef, name::Symbol, keys::Union{Int, Vector, Tuple, AbstractRange}) = set_dimension!(md, name, Dimension(keys))

function set_dimension!(md::ModelDef, name::Symbol, dim::Dimension)
    redefined = haskey(md, name)
    if redefined
        @warn "Redefining dimension :$name"
    end

    if name == :time
        k = [keys(dim)...]
        md.is_uniform = isuniform(k)
        if redefined
            reset_run_periods!(md, k[1], k[end])
        end
    end

    md.dimensions[name] = dim
    return dim
end

# helper functions used to determine if the provided time values are
# a uniform range.

function all_equal(values)
    return all(map(val -> val == values[1], values[2:end]))
end

function isuniform(values)
   if length(values) == 0
        return false
   else
        return all_equal(diff(collect(values)))
   end
end

#needed when time dimension is defined using a single integer
function isuniform(values::Int)
    return true
end

# function isuniform(values::AbstractRange{Int})
#     return isuniform(collect(values))
# end

#
# Parameters
#

external_params(md::ModelDef) = md.external_params

function addparameter(comp_def::ComponentDef, name, datatype, dimensions, description, unit, default)
    p = DatumDef(name, datatype, dimensions, description, unit, :parameter, default)
    comp_def.parameters[name] = p
    return p
end

function addparameter(comp_id::ComponentId, name, datatype, dimensions, description, unit, default)
    addparameter(compdef(comp_id), name, datatype, dimensions, description, unit, default)
end

"""
    parameters(comp_def::ComponentDef)

Return a list of the parameter definitions for `comp_def`.
"""
parameters(comp_def::ComponentDef) = values(comp_def.parameters)

"""
    parameters(comp_id::ComponentDef)

Return a list of the parameter definitions for `comp_id`.
"""
parameters(comp_id::ComponentId) = parameters(compdef(comp_id))

"""
    parameter_names(md::ModelDef, comp_name::Symbol)

Return a list of all parameter names for a given component `comp_name` in a model def `md`.
"""
parameter_names(md::ModelDef, comp_name::Symbol) = parameter_names(compdef(md, comp_name))

parameter_names(comp_def::ComponentDef) = [name(param) for param in parameters(comp_def)]

parameter(md::ModelDef, comp_name::Symbol, param_name::Symbol) = parameter(compdef(md, comp_name), param_name)

function parameter(comp_def::ComponentDef, name::Symbol)
    try
        return comp_def.parameters[name]
    catch
        error("Parameter $name was not found in component $(comp_def.name)")
    end
end

function parameter_unit(md::ModelDef, comp_name::Symbol, param_name::Symbol)
    param = parameter(md, comp_name, param_name)
    return param.unit
end

function parameter_dimensions(md::ModelDef, comp_name::Symbol, param_name::Symbol)
    param = parameter(md, comp_name, param_name)
    return param.dimensions
end

"""
    set_param!(m::ModelDef, comp_name::Symbol, name::Symbol, value, dims=nothing)

Set the parameter `name` of a component `comp_name` in a model `m` to a given `value`. The
`value` can by a scalar, an array, or a NamedAray. Optional argument 'dims' is a
list of the dimension names ofthe provided data, and will be used to check that
they match the model's index labels.
"""
function set_param!(md::ModelDef, comp_name::Symbol, param_name::Symbol, value, dims=nothing)
    # perform possible dimension and labels checks
    if isa(value, NamedArray)
        dims = dimnames(value)
    end

    if dims !== nothing
        check_parameter_dimensions(md, value, dims, param_name)
    end

    comp_param_dims = parameter_dimensions(md, comp_name, param_name)
    num_dims = length(comp_param_dims)

    comp_def = compdef(md, comp_name)
    data_type = datatype(parameter(comp_def, param_name))
    dtype = data_type == Number ? number_type(md) : data_type

    if length(comp_param_dims) > 0
        # convert the number type and, if NamedArray, convert to Array
        if dtype <: AbstractArray
            value = convert(dtype, value)
        else
            value = convert(Array{dtype, num_dims}, value)
        end

        if comp_param_dims[1] == :time
            T = eltype(value)

            if num_dims == 0
                values = value
            else
                # Want to use the first from the comp_def if it has it, if not use ModelDef
                first = first_period(md, comp_def)

                if isuniform(md)
                    _, stepsize = first_and_step(md)
                    values = TimestepArray{FixedTimestep{first, stepsize}, T, num_dims}(value)
                else
                    times = time_labels(md)
                    #use the first from the comp_def
                    first_index = findfirst(isequal(first), times)
                    values = TimestepArray{VariableTimestep{(times[first_index:end]...,)}, T, num_dims}(value)
                end
            end
        else
            values = value
        end

        set_external_array_param!(md, param_name, values, comp_param_dims)

    else # scalar parameter case
        value = convert(dtype, value)
        set_external_scalar_param!(md, param_name, value)
    end

    connect_param!(md, comp_name, param_name, param_name)
    nothing
end

#
# Variables
#
variables(comp_def::ComponentDef) = values(comp_def.variables)

variables(comp_id::ComponentId) = variables(compdef(comp_id))

function variable(comp_def::ComponentDef, var_name::Symbol)
    try
        return comp_def.variables[var_name]
    catch
        error("Variable $var_name was not found in component $(comp_def.comp_id)")
    end
end

variable(comp_id::ComponentId, var_name::Symbol) = variable(compdef(comp_id), var_name)

variable(md::ModelDef, comp_name::Symbol, var_name::Symbol) = variable(compdef(md, comp_name), var_name)

"""
    variable_names(md::ModelDef, comp_name::Symbol)

Return a list of all variable names for a given component `comp_name` in a model def `md`.
"""
variable_names(md::ModelDef, comp_name::Symbol) = variable_names(compdef(md, comp_name))

variable_names(comp_def::ComponentDef) = [name(var) for var in variables(comp_def)]

function variable_unit(md::ModelDef, comp_name::Symbol, var_name::Symbol)
    var = variable(md, comp_name, var_name)
    return var.unit
end

function variable_dimensions(md::ModelDef, comp_name::Symbol, var_name::Symbol)
    var = variable(md, comp_name, var_name)
    return var.dimensions
end

# Add a variable to a ComponentDef
function addvariable(comp_def::ComponentDef, name, datatype, dimensions, description, unit)
    v = DatumDef(name, datatype, dimensions, description, unit, :variable)
    comp_def.variables[name] = v
    return v
end

# Add a variable to a ComponentDef referenced by ComponentId
function addvariable(comp_id::ComponentId, name, datatype, dimensions, description, unit)
    addvariable(compdef(comp_id), name, datatype, dimensions, description, unit)
end

#
# Other
#

# Return the number of timesteps a given component in a model will run for.
function getspan(md::ModelDef, comp_name::Symbol)
    comp_def = compdef(md, comp_name)
    return getspan(md, comp_def)
end

function getspan(md::ModelDef, comp_def::ComponentDef)
    first = first_period(md, comp_def)
    last  = last_period(md, comp_def)
    times = time_labels(md)
    first_index = findfirst(isequal(first), times)
    last_index  = findfirst(isequal(last), times)
    return size(times[first_index:last_index])
end

function set_run_period!(comp_def::ComponentDef, first, last)
    comp_def.first = first
    comp_def.last = last
    return nothing
end

#
# Model
#
const NothingInt    = Union{Nothing, Int}
const NothingSymbol = Union{Nothing, Symbol}

function _add_anonymous_dims!(md::ModelDef, comp_def::ComponentDef)
    for (name, dim) in filter(pair -> pair[2] !== nothing, comp_def.dimensions)
        @info "Setting dimension $name to $dim"
        set_dimension!(md, name, dim)
    end
end

"""
    add_comp!(md::ModelDef, comp_def::ComponentDef; first=nothing, last=nothing, before=nothing, after=nothing)

Add the component indicated by `comp_def` to the model indcated by `md`. The component is added at the
end of the list unless one of the keywords, `first`, `last`, `before`, `after`. If the `comp_name`
differs from that in the `comp_def`, a copy of `comp_def` is made and assigned the new name.
"""
function add_comp!(md::ModelDef, comp_def::ComponentDef, comp_name::Symbol;
                   first::NothingInt=nothing, last::NothingInt=nothing,
                   before::NothingSymbol=nothing, after::NothingSymbol=nothing)

    # check that a time dimension has been set
    if !haskey(dimensions(md), :time)
        error("Cannot add component to model without first setting time dimension.")
    end

    # check that first and last are within the model's time index range
    time_index = dim_keys(md, :time)

    if first !== nothing && first < time_index[1]
        error("Cannot add component $name with first time before first of model's time index range.")
    end

    if last !== nothing && last > time_index[end]
        error("Cannot add component $name with last time after end of model's time index range.")
    end

    if before !== nothing && after !== nothing
        error("Cannot specify both 'before' and 'after' parameters")
    end

    # Check if component being added already exists
    if hascomp(md, comp_name)
        error("Cannot add two components of the same name ($comp_name)")
    end

    # Create a shallow copy of the original but with the new name
    # TBD: Why do we need to make a copy here? Sort this out.
    if compname(comp_def.comp_id) != comp_name
        comp_def = copy_comp_def(comp_def, comp_name)
    end

    set_run_period!(comp_def, first, last)

    _add_anonymous_dims!(md, comp_def)

    if before === nothing && after === nothing
        md.comp_defs[comp_name] = comp_def   # just add it to the end
    else
        new_comps = OrderedDict{Symbol, ComponentDef}()

        if before !== nothing
            if ! hascomp(md, before)
                error("Component to add before ($before) does not exist")
            end

            for i in compkeys(md)
                if i == before
                    new_comps[comp_name] = comp_def
                end
                new_comps[i] = md.comp_defs[i]
            end

        else    # after !== nothing, since we've handled all other possibilities above
            if ! hascomp(md, after)
                error("Component to add before ($before) does not exist")
            end

            for i in compkeys(md)
                new_comps[i] = md.comp_defs[i]
                if i == after
                    new_comps[comp_name] = comp_def
                end
            end
        end

        md.comp_defs = new_comps
        # println("md.comp_defs: $(md.comp_defs)")
    end

    # Set parameters to any specified defaults
    for param in parameters(comp_def)
        if param.default !== nothing
            set_param!(md, comp_name, name(param), param.default)
        end
    end

    return nothing
end

"""
    add_comp!(md::ModelDef, comp_id::ComponentId; comp_name::Symbol=comp_id.comp_name,
        first=nothing, last=nothing, before=nothing, after=nothing)

Add the component indicated by `comp_id` to the model indicated by `md`. The component is added at the end of
the list unless one of the keywords, `first`, `last`, `before`, `after`. If the `comp_name`
differs from that in the `comp_def`, a copy of `comp_def` is made and assigned the new name.
"""
function add_comp!(md::ModelDef, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name;
                   first::NothingInt=nothing, last::NothingInt=nothing,
                   before::NothingSymbol=nothing, after::NothingSymbol=nothing)
    # println("Adding component $comp_id as :$comp_name")
    add_comp!(md, compdef(comp_id), comp_name, first=first, last=last, before=before, after=after)
end

"""
    replace_comp!(md::ModelDef, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name;
        first::NothingInt=nothing, last::NothingInt=nothing,
        before::NothingSymbol=nothing, after::NothingSymbol=nothing,
        reconnect::Bool=true)

Replace the component with name `comp_name` in model definition `md` with the
component `comp_id` using the same name. The component is added in the same
position as the old component, unless one of the keywords `before` or `after`
is specified. The component is added with the same first and last values,
unless the keywords `first` or `last` are specified. Optional boolean argument
`reconnect` with default value `true` indicates whether the existing parameter
connections should be maintained in the new component.
"""
function replace_comp!(md::ModelDef, comp_id::ComponentId, comp_name::Symbol=comp_id.comp_name;
                           first::NothingInt=nothing, last::NothingInt=nothing,
                           before::NothingSymbol=nothing, after::NothingSymbol=nothing,
                           reconnect::Bool=true)

    if ! haskey(md.comp_defs, comp_name)
        error("Cannot replace '$comp_name'; component not found in model.")
    end

    # Get original position if new before or after not specified
    if before === nothing && after === nothing
        comps = collect(keys(md.comp_defs))
        n = length(comps)
        if n > 1
            idx = findfirst(isequal(comp_name), comps)
            if idx == n
                after = comps[idx - 1]
            else
                before = comps[idx + 1]
            end
        end
    end

    # Get original first and last if new run period not specified
    old_comp = md.comp_defs[comp_name]
    first = first === nothing ? old_comp.first : first
    last = last === nothing ? old_comp.last : last

    if reconnect
        # Assert that new component definition has same parameters and variables needed for the connections

        new_comp = compdef(comp_id)

        function _compare_datum(dict1, dict2)
            set1 = Set([(k, v.datatype, v.dimensions) for (k, v) in dict1])
            set2 = Set([(k, v.datatype, v.dimensions) for (k, v) in dict2])
            return set1 >= set2
        end

        # Check incoming parameters
        incoming_params = map(ipc -> ipc.dst_par_name, internal_param_conns(md, comp_name))
        old_params = filter(pair -> pair.first in incoming_params, old_comp.parameters)
        new_params = new_comp.parameters
        if !_compare_datum(new_params, old_params)
            error("Cannot replace and reconnect; new component does not contain the same definitions of necessary parameters.")
        end

        # Check outgoing variables
        outgoing_vars = map(ipc -> ipc.src_var_name, filter(ipc -> ipc.src_comp_name == comp_name, md.internal_param_conns))
        old_vars = filter(pair -> pair.first in outgoing_vars, old_comp.variables)
        new_vars = new_comp.variables
        if !_compare_datum(new_vars, old_vars)
            error("Cannot replace and reconnect; new component does not contain the same definitions of necessary variables.")
        end

        # Check external parameter connections
        remove = []
        for epc in external_param_conns(md, comp_name)
            param_name = epc.param_name
            if ! haskey(new_params, param_name)  # TODO: is this the behavior we want? don't error in this case? just (warn)?
                @warn "Removing external parameter connection from component $comp_name; parameter $param_name no longer exists in component."
                push!(remove, epc)
            else
                old_p = old_comp.parameters[param_name]
                new_p = new_params[param_name]
                if new_p.dimensions != old_p.dimensions
                    error("Cannot replace and reconnect; parameter $param_name in new component has different dimensions.")
                end
                if new_p.datatype != old_p.datatype
                    error("Cannot replace and reconnect; parameter $param_name in new component has different datatype.")
                end
            end
        end
        filter!(epc -> !(epc in remove), md.external_param_conns)

        # Delete the old component from comp_defs, leaving the existing parameter connections
        delete!(md.comp_defs, comp_name)
    else
        # Delete the old component and all its internal and external parameter connections
        delete!(md, comp_name)
    end

    # Re-add
    add_comp!(md, comp_id, comp_name; first=first, last=last, before=before, after=after)

end

"""
    copy_comp_def(comp_def::ComponentDef, comp_name::Symbol)

Create a mostly-shallow copy of `comp_def` (named `comp_name`), but make a deep copy of its
ComponentId so we can rename the copy without affecting the original.
"""
function copy_comp_def(comp_def::ComponentDef, comp_name::Symbol)
    comp_id = comp_def.comp_id
    obj     = ComponentDef(comp_id)

    # Use the comp_id as is, since this identifies the run_timestep function, but
    # use an alternate name to reference it in the model's component list.
    obj.name = comp_name

    obj.variables  = comp_def.variables
    obj.parameters = comp_def.parameters
    obj.dimensions = comp_def.dimensions
    obj.first      = comp_def.first
    obj.last       = comp_def.last

    return obj
end

"""
    copy_external_params(md::ModelDef)

Make copies of ModelParameter subtypes representing external parameters of model `md`.
This is used both in the copy() function below, and in the MCS subsystem
to restore values between trials.

"""
function copy_external_params(md::ModelDef)
    external_params = Dict{Symbol, ModelParameter}(key => copy(obj) for (key, obj) in md.external_params)
    return external_params
end

Base.copy(obj::ScalarModelParameter{T}) where T = ScalarModelParameter{T}(copy(obj.value))

Base.copy(obj::ArrayModelParameter{T})  where T = ArrayModelParameter{T}(copy(obj.values), obj.dimensions)

function Base.copy(obj::TimestepVector{T_ts, T}) where {T_ts, T}
    return TimestepVector{T_ts, T}(copy(obj.data))
end

function Base.copy(obj::TimestepMatrix{T_ts, T}) where {T_ts, T}
    return TimestepMatrix{T_ts, T}(copy(obj.data))
end

function Base.copy(obj::TimestepArray{T_ts, T, N}) where {T_ts, T, N}
    return TimestepArray{T_ts, T, N}(copy(obj.data))
end

"""
    copy(md::ModelDef)

Create a copy of a ModelDef `md` object that is not entirely shallow, nor completely deep.
The aim is to copy the full structure, reusing references to immutable elements.
"""
function Base.copy(md::ModelDef)
    mdcopy = ModelDef(md.number_type)
    mdcopy.module_name = md.module_name

    merge!(mdcopy.comp_defs, md.comp_defs)

    mdcopy.dimensions = deepcopy(md.dimensions)

    # These are vectors of immutable structs, so we can (shallow) copy them safely
    mdcopy.internal_param_conns = copy(md.internal_param_conns)
    mdcopy.external_param_conns = copy(md.external_param_conns)

    # Names of external params that the ConnectorComps will use as their :input2 parameters.
    mdcopy.backups = copy(md.backups)
    mdcopy.external_params = copy_external_params(md)

    mdcopy.sorted_comps = md.sorted_comps === nothing ? nothing : copy(md.sorted_comps)

    mdcopy.is_uniform = md.is_uniform

    return mdcopy
end


text:
const MP = MultivariatePolynomials
@testset "Division" begin
    @testset "GCD and LCM" begin
        Mod.@polyvar x y z
        @test gcd(x*y, y) == y
        @test lcm(x*y, y) == x*y
        @test gcd(x*y^2*z, y*z^3) == y*z
        @test lcm(x*y^2*z, y*z^3) == x*y^2*z^3
        @test gcd(x^2*y^7*z^3, x^4*y^5*z^2) == x^2*y^5*z^2
        @test lcm(x^2*y^7*z^3, x^4*y^5*z^2) == x^4*y^7*z^3
    end
    # Taken from
    # Ideals, Varieties, and Algorithms
    # Cox, Little and O'Shea, Fourth edition
    # They have been adapted to the grlex ordering
    @testset "Divides" begin
        Mod.@polyvar x y z
        @test divides(leadingmonomial(x+y), x) # leadingmonomial(x+y) will be x^1*y^0 -> tricky test !
        @test !divides(leadingmonomial(x^2+y), x)
        @test divides(x*y, x^2*y)
        @test divides(x*y, x*y^2)
        @test divides(y*z, x*y*z)
        @test !divides(y*z, x*z)
        @test !divides(x^2*y, x*y^2)
    end
    @testset "Leading function" begin
        Mod.@polyvar x y z
        @test @inferred(monic(2x^2 + 4y + 2)) == x^2 + 2y + 1
        # See page 60
        p = 4x*y^2*z + 4z^2 + 7x^2*z^2 - 5x^3
        @test leadingcoefficient(p) == 7
        @test leadingmonomial(p) == x^2*z^2
        @test leadingterm(p) == 7x^2*z^2
    end
    @testset "Division examples" begin
        Mod.@polyvar x y
        @test iszero(@inferred MP.proddiff(2x*y, 3y^2*x))
        @test (@inferred div(x*y^2 + 1, x*y + 1)) == y
        @test (@inferred rem(x*y^2 + 1, x*y + 1)) == -y + 1
        @test (@inferred div(x*y^2 + x, y)) == x*y
        @test (@inferred rem(x*y^2 + x, y)) == x
        @test (@inferred rem(x^4 + x^3 + (1+1e-10)*x^2 + 1, x^2 + x + 1)) == 1
        @test (@inferred rem(x^4 + x^3 + (1+1e-10)*x^2 + 1, x^2 + x + 1; ztol=1e-11)) ≈ -((1+1e-10)-1)x + 1
    end
    @testset "Division by multiple polynomials examples" begin
        function testdiv(p, ps)
            q, r = @inferred divrem(p, ps)
            @test p == dot(q, ps) + r
            q, r
        end
        Mod.@polyvar x y
        # Example 1
        q, r = testdiv(x*y^2 + 1, [x*y + 1, y + 1])
        @test q == [y, -1]
        @test r == 2
        # Example 2
        q, r = testdiv(x^2*y + x*y^2 + y^2, [x*y - 1, y^2 - 1])
        @test q == [x + y, 1]
        @test r == x + y + 1
        # Example 4
        q, r = testdiv(x^2*y + x*y^2 + y^2, [y^2 - 1, x*y - 1])
        @test q == [x + 1, x]
        @test r == 2x + 1
        # Example 5
        q, r = testdiv(x*y^2 - x, [x*y - 1, y^2 - 1])
        @test q == [y, 0]
        @test r == -x + y
        q, r = testdiv(x*y^2 - x, [y^2 - 1, x*y - 1])
        @test q == [x, 0]
        @test r == 0

        @test (@inferred rem(x^2*y + (1+1e-10)*x*y + 1, [x^2 + x, y + 1])) == 1
        @test (@inferred rem(x^2*y + (1+1e-10)*x*y + 1, [x^2 + x, y + 1]; ztol=1e-11)) == -((1+1e-10)-1)x + 1
    end
    @testset "gcd" begin
        Mod.@polyvar x
        @test -2(x + 1) == @inferred gcd(x^2 - 1, x^2 + 2x + 1)
        @test  2(x + 1) == @inferred gcd(x^2 + 2x + 1, x^2 - 1)
        @test    x + 1  == @inferred gcd(x^2 - 1, x + 1)
        @test    x + 1  == @inferred gcd(x + 1, x^2 - 1)
        @test    x + 1  == @inferred gcd(x - x, x + 1)
        @test    x + 1  == @inferred gcd(x + 1, x - x)
        @test        0  == @inferred gcd(x - x, x^2 - x^2)
        @test        0  == @inferred gcd(x^2 - x^2, x - x)
    end
    @testset "lcm" begin
        Mod.@polyvar x
        @test -(x-1) * (x+1)^2 / 2 == @inferred lcm(x^2 - 1, x^2 + 2x + 1)
        @test  (x-1) * (x+1)^2 / 2 == @inferred lcm(x^2 + 2x + 1, x^2 - 1)
        @test x^2 - 1 == @inferred lcm(x^2 - 1, x - 1)
        @test x^2 - 1 == @inferred lcm(x - 1, x^2 - 1)
        @test       0 == @inferred lcm(x - x, x + 1)
        @test       0 == @inferred lcm(x + 1, x - x)
        @test       0 == @inferred lcm(x - x, x^2 - x^2)
        @test       0 == @inferred lcm(x^2 - x^2, x - x)
    end
end


text:
export SineWave

"""
    SineWave(amplitude, period, phase)
    SineWave(:from_phase, amplitude, period, phase)
    SineWave(:from_peak, amplitude, period, peak_time)

Creates a sine wave, i.e., something of the form

    ϕ(t) = A sin( 2π t / T + θ₀),

where `A` is the `amplitude`, `T` is the period, and `θ₀` is the phase.
"""
struct SineWave <: Sequence
  amplitude
  period
  phase
end

function SineWave(type::Symbol, amplitude, period, third)
  if type == :from_phase
    return SineWave(amplitude, period, third)
  elseif type == :from_peak
    phase = π / 2 - 2π * third / period
    return SineWave(amplitude, period, phase)
  end
end

function sample(seq::SineWave, t::AbstractVector)
  return seq.amplitude * sin.(2π * t / seq.period .+ seq.phase)
end