text:
function abs(obj::DerivativeStructure)
    return jcall(obj, "abs", DerivativeStructure, ())
end

function acos(obj::DerivativeStructure)
    return jcall(obj, "acos", DerivativeStructure, ())
end

function acosh(obj::DerivativeStructure)
    return jcall(obj, "acosh", DerivativeStructure, ())
end

function add(obj::DerivativeStructure, arg0::DerivativeStructure)
    return jcall(obj, "add", DerivativeStructure, (DerivativeStructure,), arg0)
end

function add(obj::DerivativeStructure, arg0::jdouble)
    return jcall(obj, "add", DerivativeStructure, (jdouble,), arg0)
end

function asin(obj::DerivativeStructure)
    return jcall(obj, "asin", DerivativeStructure, ())
end

function asinh(obj::DerivativeStructure)
    return jcall(obj, "asinh", DerivativeStructure, ())
end

function atan(obj::DerivativeStructure)
    return jcall(obj, "atan", DerivativeStructure, ())
end

function atan2(::Type{DerivativeStructure}, arg0::DerivativeStructure, arg1::DerivativeStructure)
    return jcall(DerivativeStructure, "atan2", DerivativeStructure, (DerivativeStructure, DerivativeStructure), arg0, arg1)
end

function atan2(obj::DerivativeStructure, arg0::DerivativeStructure)
    return jcall(obj, "atan2", DerivativeStructure, (DerivativeStructure,), arg0)
end

function atanh(obj::DerivativeStructure)
    return jcall(obj, "atanh", DerivativeStructure, ())
end

function cbrt(obj::DerivativeStructure)
    return jcall(obj, "cbrt", DerivativeStructure, ())
end

function ceil(obj::DerivativeStructure)
    return jcall(obj, "ceil", DerivativeStructure, ())
end

function compose(obj::DerivativeStructure, arg0::Vector{jdouble})
    return jcall(obj, "compose", DerivativeStructure, (Vector{jdouble},), arg0)
end

function copy_sign(obj::DerivativeStructure, arg0::DerivativeStructure)
    return jcall(obj, "copySign", DerivativeStructure, (DerivativeStructure,), arg0)
end

function copy_sign(obj::DerivativeStructure, arg0::jdouble)
    return jcall(obj, "copySign", DerivativeStructure, (jdouble,), arg0)
end

function cos(obj::DerivativeStructure)
    return jcall(obj, "cos", DerivativeStructure, ())
end

function cosh(obj::DerivativeStructure)
    return jcall(obj, "cosh", DerivativeStructure, ())
end

function divide(obj::DerivativeStructure, arg0::DerivativeStructure)
    return jcall(obj, "divide", DerivativeStructure, (DerivativeStructure,), arg0)
end

function divide(obj::DerivativeStructure, arg0::jdouble)
    return jcall(obj, "divide", DerivativeStructure, (jdouble,), arg0)
end

function equals(obj::DerivativeStructure, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function exp(obj::DerivativeStructure)
    return jcall(obj, "exp", DerivativeStructure, ())
end

function expm1(obj::DerivativeStructure)
    return jcall(obj, "expm1", DerivativeStructure, ())
end

function floor(obj::DerivativeStructure)
    return jcall(obj, "floor", DerivativeStructure, ())
end

function get_all_derivatives(obj::DerivativeStructure)
    return jcall(obj, "getAllDerivatives", Vector{jdouble}, ())
end

function get_exponent(obj::DerivativeStructure)
    return jcall(obj, "getExponent", jint, ())
end

function get_factory(obj::DerivativeStructure)
    return jcall(obj, "getFactory", DSFactory, ())
end

function get_field(obj::DerivativeStructure)
    return jcall(obj, "getField", Field, ())
end

function get_free_parameters(obj::DerivativeStructure)
    return jcall(obj, "getFreeParameters", jint, ())
end

function get_order(obj::DerivativeStructure)
    return jcall(obj, "getOrder", jint, ())
end

function get_partial_derivative(obj::DerivativeStructure, arg0::Vector{jint})
    return jcall(obj, "getPartialDerivative", jdouble, (Vector{jint},), arg0)
end

function get_real(obj::DerivativeStructure)
    return jcall(obj, "getReal", jdouble, ())
end

function get_value(obj::DerivativeStructure)
    return jcall(obj, "getValue", jdouble, ())
end

function hash_code(obj::DerivativeStructure)
    return jcall(obj, "hashCode", jint, ())
end

function hypot(::Type{DerivativeStructure}, arg0::DerivativeStructure, arg1::DerivativeStructure)
    return jcall(DerivativeStructure, "hypot", DerivativeStructure, (DerivativeStructure, DerivativeStructure), arg0, arg1)
end

function hypot(obj::DerivativeStructure, arg0::DerivativeStructure)
    return jcall(obj, "hypot", DerivativeStructure, (DerivativeStructure,), arg0)
end

function linear_combination(obj::DerivativeStructure, arg0::DerivativeStructure, arg1::DerivativeStructure, arg2::DerivativeStructure, arg3::DerivativeStructure)
    return jcall(obj, "linearCombination", DerivativeStructure, (DerivativeStructure, DerivativeStructure, DerivativeStructure, DerivativeStructure), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::DerivativeStructure, arg0::DerivativeStructure, arg1::DerivativeStructure, arg2::DerivativeStructure, arg3::DerivativeStructure, arg4::DerivativeStructure, arg5::DerivativeStructure)
    return jcall(obj, "linearCombination", DerivativeStructure, (DerivativeStructure, DerivativeStructure, DerivativeStructure, DerivativeStructure, DerivativeStructure, DerivativeStructure), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::DerivativeStructure, arg0::DerivativeStructure, arg1::DerivativeStructure, arg2::DerivativeStructure, arg3::DerivativeStructure, arg4::DerivativeStructure, arg5::DerivativeStructure, arg6::DerivativeStructure, arg7::DerivativeStructure)
    return jcall(obj, "linearCombination", DerivativeStructure, (DerivativeStructure, DerivativeStructure, DerivativeStructure, DerivativeStructure, DerivativeStructure, DerivativeStructure, DerivativeStructure, DerivativeStructure), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function linear_combination(obj::DerivativeStructure, arg0::Vector{DerivativeStructure}, arg1::Vector{DerivativeStructure})
    return jcall(obj, "linearCombination", DerivativeStructure, (Vector{DerivativeStructure}, Vector{DerivativeStructure}), arg0, arg1)
end

function linear_combination(obj::DerivativeStructure, arg0::Vector{jdouble}, arg1::Vector{DerivativeStructure})
    return jcall(obj, "linearCombination", DerivativeStructure, (Vector{jdouble}, Vector{DerivativeStructure}), arg0, arg1)
end

function linear_combination(obj::DerivativeStructure, arg0::jdouble, arg1::DerivativeStructure, arg2::jdouble, arg3::DerivativeStructure)
    return jcall(obj, "linearCombination", DerivativeStructure, (jdouble, DerivativeStructure, jdouble, DerivativeStructure), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::DerivativeStructure, arg0::jdouble, arg1::DerivativeStructure, arg2::jdouble, arg3::DerivativeStructure, arg4::jdouble, arg5::DerivativeStructure)
    return jcall(obj, "linearCombination", DerivativeStructure, (jdouble, DerivativeStructure, jdouble, DerivativeStructure, jdouble, DerivativeStructure), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::DerivativeStructure, arg0::jdouble, arg1::DerivativeStructure, arg2::jdouble, arg3::DerivativeStructure, arg4::jdouble, arg5::DerivativeStructure, arg6::jdouble, arg7::DerivativeStructure)
    return jcall(obj, "linearCombination", DerivativeStructure, (jdouble, DerivativeStructure, jdouble, DerivativeStructure, jdouble, DerivativeStructure, jdouble, DerivativeStructure), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function log(obj::DerivativeStructure)
    return jcall(obj, "log", DerivativeStructure, ())
end

function log10(obj::DerivativeStructure)
    return jcall(obj, "log10", DerivativeStructure, ())
end

function log1p(obj::DerivativeStructure)
    return jcall(obj, "log1p", DerivativeStructure, ())
end

function multiply(obj::DerivativeStructure, arg0::DerivativeStructure)
    return jcall(obj, "multiply", DerivativeStructure, (DerivativeStructure,), arg0)
end

function multiply(obj::DerivativeStructure, arg0::jdouble)
    return jcall(obj, "multiply", DerivativeStructure, (jdouble,), arg0)
end

function multiply(obj::DerivativeStructure, arg0::jint)
    return jcall(obj, "multiply", DerivativeStructure, (jint,), arg0)
end

function negate(obj::DerivativeStructure)
    return jcall(obj, "negate", DerivativeStructure, ())
end

function new_instance(obj::DerivativeStructure, arg0::jdouble)
    return jcall(obj, "newInstance", DerivativeStructure, (jdouble,), arg0)
end

function pow(::Type{DerivativeStructure}, arg0::jdouble, arg1::DerivativeStructure)
    return jcall(DerivativeStructure, "pow", DerivativeStructure, (jdouble, DerivativeStructure), arg0, arg1)
end

function pow(obj::DerivativeStructure, arg0::DerivativeStructure)
    return jcall(obj, "pow", DerivativeStructure, (DerivativeStructure,), arg0)
end

function pow(obj::DerivativeStructure, arg0::jdouble)
    return jcall(obj, "pow", DerivativeStructure, (jdouble,), arg0)
end

function pow(obj::DerivativeStructure, arg0::jint)
    return jcall(obj, "pow", DerivativeStructure, (jint,), arg0)
end

function reciprocal(obj::DerivativeStructure)
    return jcall(obj, "reciprocal", DerivativeStructure, ())
end

function remainder(obj::DerivativeStructure, arg0::DerivativeStructure)
    return jcall(obj, "remainder", DerivativeStructure, (DerivativeStructure,), arg0)
end

function remainder(obj::DerivativeStructure, arg0::jdouble)
    return jcall(obj, "remainder", DerivativeStructure, (jdouble,), arg0)
end

function rint(obj::DerivativeStructure)
    return jcall(obj, "rint", DerivativeStructure, ())
end

function root_n(obj::DerivativeStructure, arg0::jint)
    return jcall(obj, "rootN", DerivativeStructure, (jint,), arg0)
end

function scalb(obj::DerivativeStructure, arg0::jint)
    return jcall(obj, "scalb", DerivativeStructure, (jint,), arg0)
end

function signum(obj::DerivativeStructure)
    return jcall(obj, "signum", DerivativeStructure, ())
end

function sin(obj::DerivativeStructure)
    return jcall(obj, "sin", DerivativeStructure, ())
end

function sin_cos(obj::DerivativeStructure)
    return jcall(obj, "sinCos", FieldSinCos, ())
end

function sinh(obj::DerivativeStructure)
    return jcall(obj, "sinh", DerivativeStructure, ())
end

function sqrt(obj::DerivativeStructure)
    return jcall(obj, "sqrt", DerivativeStructure, ())
end

function subtract(obj::DerivativeStructure, arg0::DerivativeStructure)
    return jcall(obj, "subtract", DerivativeStructure, (DerivativeStructure,), arg0)
end

function subtract(obj::DerivativeStructure, arg0::jdouble)
    return jcall(obj, "subtract", DerivativeStructure, (jdouble,), arg0)
end

function tan(obj::DerivativeStructure)
    return jcall(obj, "tan", DerivativeStructure, ())
end

function tanh(obj::DerivativeStructure)
    return jcall(obj, "tanh", DerivativeStructure, ())
end

function taylor(obj::DerivativeStructure, arg0::Vector{jdouble})
    return jcall(obj, "taylor", jdouble, (Vector{jdouble},), arg0)
end

function to_degrees(obj::DerivativeStructure)
    return jcall(obj, "toDegrees", DerivativeStructure, ())
end

function to_radians(obj::DerivativeStructure)
    return jcall(obj, "toRadians", DerivativeStructure, ())
end


text:
"""

Plotting of symbolic objects.

The `Plots` package provide a uniform interface to many of `Julia`'s
plotting packages. `SymPy` plugs into `Plots`' "recipes."

The basic goal is that when `Plots` provides an interface for function
objects, this package extends the interface to symbolic expressions.

In particular:

* `plot(ex::Sym, a, b; kwargs...)` will plot a function evaluating `ex` over [a,b]

Example. Here we use the default backend for `Plots` to make a plot:

```
using Plots
@syms x
plot(x^2 - 2x, 0, 4)
```

* `plot(ex1, ex2, a, b; kwargs...)` will plot the two expressions in a parametric plot over the interval `[a,b]`.

Example:

```
@syms x
plot(sin(2x), cos(3x), 0, 4pi) ## also
```

For a few backends (those that support `:path3d`) a third symbolic
expression may be added to have a 3d parametric plot rendered:

```
plot(sin(x), cos(x), x, 0, 4pi) # helix in 3d
```

* `plot(xs, ys, expression)` will make a contour plot (for many backends).

```
@syms x y
plot(range(0,stop=5, length=50), range(0,stop=5, length=50), x*y)
```

* To plot the surface  `z=ex(x,y)` over a region we have `Plots.surface`. For example,

```
@syms x y
surface(-5:5, -5:5, 25 - x^2 - y^2)
```

* a vectorfield plot can (inefficiently but directly) be produced following this example:

```
function vfieldplot(fx, fy; xlim=(-5,5), ylim=(-5,5), n=8)
    xs = range(xlim[1], stop=xlim[2], length=n)
    ys = range(ylim[1], stop=ylim[2], length=n)

    us = vec([x for x in xs, y in ys])
    vs = vec([y for x in xs, y in ys])
    fxs = vec([fx(x,y) for x in xs, y in ys])
    fys = vec([fy(x,y) for x in xs, y in ys])

    mxs = maximum(abs.(filter(!isinf, filter(!isnan, fxs))))
    mys = maximum(abs.(filter(!isinf, filter(!isnan, fys))))
    d = 1/2 * max((xlim[2]-xlim[1])/mxs, (ylim[2]-ylim[1])/mys) / n

    quiver(us, vs, quiver=(fxs.*d, fys.*d))

end
fx = (x + y) / sqrt(x^2 + y^2)
fy = (x - y) / sqrt(x^2 + y^2)
vfieldplot(fx, fy)

```

* To plot two or more functions at once, the style `plot([ex1, ex2], a, b)` does not work. Rather, use
    `plot(ex1, a, b); plot!(ex2)`, as in:

```
@syms x
plot(sin(x), 0, 2pi)
plot!(cos(x))
```
----

Some graphics provided by `SymPy` are available if `PyPlot` is installed, such as:

* `sympy.plotting.plot3d_parametric_surface`
* `sympy.plotting.plot_implicit`

Plot the parametrically defined surface `[exs[1](u,v), exs[2](u,v), exs[3](u,v)]` over `[a0,a1] x
[b0,b1]`. The specification of the variables uses a tuple of the form
`(Sym, Real, Real)` following the style of SymPy in `integrate`, say,
where disambiguation of variable names is needed.

```
@syms theta, phi
r = 1
sympy.plotting.plot3d_parametric_surface((r*sin(theta)*sin(phi), r*sin(theta)*cos(phi), r*cos(theta)),
                        (theta, 0, pi), (phi, 0, pi/2))
```

* `sympy.plotting.plot_implicit(equation, (xvar, x0, x1), (yvar, y0, y1))` will plot implicitly the equation.

```
@syms x y
sympy.plotting.plot_implicit(Eq(x^2+ y^2,3), (x, -2, 2), (y, -2, 2))  # draw a circle
```

"""
sympy_plotting = nothing
export sympy_plotting

## Recipes for hooking into Plots

##
@recipe f(::Type{T}, v::T) where {T<:Sym} = lambdify(v)

## for vectors of expressions
## This does not work. See: https://github.com/JuliaPlots/RecipesBase.jl/issues/19
#@recipe f(ss::AbstractVector{Sym}) = lambdify.(ss)
#@recipe  function f{T<:Array{Sym,1}}(::Type{T}, ss::T)  Function[lambdify(s) for s in ss]  end

## A vector field plot can be visualized as an n × n collection of arrows
## over the region xlims × ylims
## These arrows are defined by:
## * fx, fy giving the components of each. These are callbable objects, such as
##   (x,y) -> sin(y)
## * Fyx  A function F(y,x), useful for visualizing first-order ODE y'=F(y(x),x).
##   note reverse order of y and x.
## The vectors are scaled so as not to overlap.

"""
`VectorField(fx, fy`): create an object that can be `plot`ted as a vector field.

A vectorfield plot draws arrows at grid points proportional to `[fx(x_i,y_i), fy(x_i,y_i)]` to visualize the field generated by `[fx, fy]`.

The plot command: `plot(VectorField(fx, fy), xlims=(-5,5), ylims=(-5,5), n=8)` will draw the vectorfield. This uses the default values, so the same graph would be rendered by `plot(VectorField(fx,fy))`.

To faciliate the visualization of solution to the ODE y' = F(x, y(x)), the call `plot(VectorField(F))` will work. (The order is x then y, though often this is written as F(y(x),x).)

`SymPy` objects can be passed to `VectorField`, but this is a bit
fragile, as they must each have two variables so that they can be
called with two variables.  (E.g., `y(1,2)` will be `1` not `2`, as
might be intended.)

Examples:
```
using Plots

fx(x,y) = sin(y); fy(x,y) = cos(y)
plot(VectorField(fx, fy), xlims=(-2pi, 2pi), ylims=(-2pi,2pi))

# plot field of y' = 3y*x over (-5,5) x (-5,5)
F(x,y) = 3*y*x
plot(VectorField(F))

# plot field and solution u' = u*(1-u)
@syms x u()
F(x,y) = y*(1-y)
out = dsolve(u'(x) - F(x, u(x)), x, (u, 0, 1))
plot(VectorField(F), xlims=(0,5), ylims=(0,2))
plot!(rhs(out))
```

"""
struct VectorField
    fx
    fy
end
VectorField(f) = VectorField((x,y) -> 1.0, f)
export VectorField

@recipe function f(F::VectorField; n=8)

    xlims = get(plotattributes,:xlims, (-5,5))
    ylims = get(plotattributes, :ylims, (-5,5))

    xs = repeat(range(xlims[1], stop=xlims[2], length=n), inner=(n,))
    ys = repeat(range(ylims[1], stop=ylims[2], length=n), outer=(n,))

    us, vs = broadcast(F.fx, xs, ys), broadcast(F.fy, xs, ys)

    delta = min((xlims[2]-xlims[1])/n, (ylims[2]-ylims[1])/n)
    m = maximum([norm([u,v]) for (u,v) in zip(us, vs)])

    lambda = delta / m

    x := xs
    y := ys
    quiver := (lambda*us, lambda*vs)
    seriestype := :quiver
    ()
end

## ---------------------

## These functions give acces to SymPy's plotting module. They will work if PyPlot is installed, but may otherwise cause an error

## surface plot xvar = Tuple(Sym, Real, Real)
##
"""

Render a parametrically defined surface plot.

Example:
```
@syms u, v
plot_parametric_surface((u*v,u-v,u+v), (u,0,1), (v,0,1))
```

This uses `PyPlot`, not `Plots` for now.
"""
function plot_parametric_surface(exs::Tuple{Sym,Sym,Sym},
                                 xvar=(-5.0, 5.0),
                                 yvar=(-5.0, 5.0),
                                 args...;
                                 kwargs...)

    sympy.plotting.plot3d_parametric_surface(exs..., args...; kwargs...)

end
export plot_parametric_surface

"""
Plot an implicit equation

```
@syms x y
plot_implicit(Eq(x^2+ y^2,3), (x, -2, 2), (y, -2, 2))
```

"""
plot_implicit(ex, args...; kwargs...) = sympy.plotting.plot_implicit(ex, args...; kwargs...)
export plot_implicit

ext: jl

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
module GLText
using ModernGL
using GLAbstraction
using GLFW
using Compat

export getfont, build_line_indexes, findline

function getfont()
	standardFont
end

function inittext()
	global standardFont = GLFont(Pkg.dir("GLText", "src", "VeraMono"))
end
init_after_context_creation(inittext)

include("types.jl")
#include("textField.jl")
#include("render.jl")

end # module


text:
# Modified from Shreyas Kowshik's implementation.

@with_kw mutable struct PPOSolver
    episode_length::Int64 = 100 # Length of each episode
    terminate_horizon::Bool = false # Zero reward at end of horizon
    resume::Bool = false # Load saved policy and resume.

    # Optimization parameters
    η::Float64 = 1e-3 # learning rate η
    optimizer = ADAM(η) # Optimizer for the value function and policy neural network

    # Generalized advantage estimate (GAE) parameters
    gae_γ::Float64 = 0.99 # GAE-Gamma
    gae_λ::Float64 = 0.95 # GAE-Lambda

    # Training parameters
    num_episodes::Int64 = 10 # Number of epochs of interaction (equal to number of policy updates)
    ppo_epochs::Int64 = 10 # Number of steps of interaction (state-action pairs) for the agent and the environment in each episode
    batch_size::Int64 = episode_length # Size of the batches to perfom the update on
    cₚ::Float64 = 1.0 # policy loss coefficient
    cᵥ::Float64 = 1.0 # value loss coefficient
    cₑ::Float64 = 0.01 # entropy loss coefficient

    # PPO parameters
    ϵ::Float64 = 0.1 # PPO gradient clipping

    # Policy parameters
    hidden_layer_size::Int64 = 30 # Policy and value hidden layer size
    policy_type::Symbol = :continuous # :discrete or :continuous (TODO. infer)
    log_std::Float64 = log(1) # Initial log standard deviation used in the log-STD network

    # Environment parameters
    state_size::Int64 = 1 # Size of the state/observation space.
    action_size::Int64 = 1 # Size of the action space.

    # Frequencies and verbosity
    save::Bool = false
    save_frequency::Int64 = 10
    verbose_frequency::Int64 = 5
    verbose::Bool = false
    show_progress::Bool = true
end

mutable struct PPOPlanner{P<:Union{MDP,POMDP}}
    solver::PPOSolver
    mdp::P
    env::Union{MDPEnvironment,POMDPEnvironment}
    policy::Union{CategoricalPolicy,DiagonalGaussianPolicy,Nothing}
end

# Uses RLInterface environment wrapper
POMDPs.solve(solver::PPOSolver, mdp::MDP) = PPOPlanner(set_action_size!(solver, mdp), mdp, MDPEnvironment(mdp), nothing)
POMDPs.solve(solver::PPOSolver, mdp::POMDP) = PPOPlanner(set_action_size!(solver, mdp), mdp, POMDPEnvironment(mdp), nothing)

# NOTE: action and playback shared between TRPO and PPO (see PolicyOptimization.jl)

function loss(solver::PPOSolver, policy, states::Array, actions::Array, advantages::Array, returns::Array, old_log_probs::Array)
    new_log_probs = log_prob(policy, states, actions)

    # Surrogate loss computations
    ratio = exp.(new_log_probs .- old_log_probs)
    surr1 = ratio .* advantages
    surr2 = clamp.(ratio, (1.0 - solver.ϵ), (1.0 + solver.ϵ)) .* advantages
    policy_loss = mean(min.(surr1, surr2))

    value_predicted = policy.value_net(states)
    value_loss = mean((value_predicted .- returns).^2)

    entropy_loss = mean(entropy(policy, states))

    return -solver.cₚ*policy_loss + solver.cᵥ*value_loss - solver.cₑ*entropy_loss
end

function ppo_update!(solver::PPOSolver, policy, states::Array, actions::Array, advantages::Array, returns::Array, old_log_probs::Array, kl_vars)
    model_params = Flux.params(get_policy_params(policy)..., get_value_params(policy)...)

    # Calculate gradients
    gs = Tracker.gradient(() -> loss(solver, policy, states, actions, advantages, returns, old_log_probs), model_params)

    # Take a step of optimization
    update!(solver.optimizer, model_params, gs)
end


text:
function getcentroids(problem::AbstractTopOptProblem)
    dh = problem.ch.dh
    return map(CellIterator(dh)) do cell
        return Vector(mean(cell.coords))
    end
end

abstract type AbstractMLModel end

# in -- params --> out -> filter -> compliance === loss
# params -- in --> out -> filter -> compliance === loss

@params struct NeuralNetwork <: AbstractMLModel
    model::Any
    init_params::Any
    params_to_out::Any
    in_to_out::Any
end
function NeuralNetwork(nn_model, input_coords::AbstractVector)
    f = x -> nn_model(x)[1]
    @assert all(0 .<= f.(input_coords) .<= 1)
    p, re = Flux.destructure(nn_model)
    return NeuralNetwork(
        nn_model, Float64.(p), p -> getindex.(re(p).(input_coords), 1), nn_model
    )
end
function NeuralNetwork(nn_model, problem::AbstractTopOptProblem)
    centroids = getcentroids(problem)
    m, s = mean(centroids), std(centroids)
    scentroids = map(centroids) do c
        (c .- m) ./ s
    end
    return NeuralNetwork(nn_model, scentroids)
end

@params struct PredictFunction <: Function
    model::AbstractMLModel
end
function (pf::PredictFunction)(in)
    return pf.model.in_to_out(in)
end

@params struct TrainFunction <: Function
    model::AbstractMLModel
end
function (tf::TrainFunction)(p)
    return tf.model.params_to_out(p)
end


text:
module TestParseWDI

using Test
using WorldBankData

us_gnp_data = Any[  Dict( "total"=>23,"per_page"=>"25000","pages"=>1,"page"=>1 ),
                 Any[ Dict("date"=>"2012","value"=>"52340","indicator"=>Dict("id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"),"country"=>Dict("id"=>"US","value"=>"United States"),"decimal"=>"0"),
                   Dict("date"=>"2011","value"=>"50650","indicator"=>Dict("id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"),"country"=>Dict("id"=>"US","value"=>"United States"),"decimal"=>"0"),
                   Dict("date"=>"2010","value"=>"48960","indicator"=>Dict("id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"),"country"=>Dict("id"=>"US","value"=>"United States"),"decimal"=>"0"),
                   Dict("date"=>"2009","value"=>"48040","indicator"=>Dict("id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"),"country"=>Dict("id"=>"US","value"=>"United States"),"decimal"=>"0"),
                   Dict("date"=>"2008","value"=>"49350","indicator"=>Dict("id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"),"country"=>Dict("id"=>"US","value"=>"United States"),"decimal"=>"0"),
                   Dict("date"=>"2007","value"=>"48640","indicator"=>Dict("id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"),"country"=>Dict("id"=>"US","value"=>"United States"),"decimal"=>"0"),
                   Dict("date"=>"2006","value"=>"48080","indicator"=>Dict("id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"),"country"=>Dict("id"=>"US","value"=>"United States"),"decimal"=>"0")
                 ]
              ]

us_gnp = WorldBankData.parse_wdi("NY.GNP.PCAP.CD",us_gnp_data[2],2006,2012)

@test us_gnp[!, :year] == Float64[2012, 2011, 2010, 2009, 2008, 2007, 2006]
@test us_gnp[!, :NY_GNP_PCAP_CD] == Float64[52340, 50650, 48960, 48040, 49350, 48640, 48080]

end


text:
# using Compose
# stroke = Compose.stroke;
# rectangle = Compose.rectangle;
# circle = Compose.circle;
# UB = (0,0,1,1);
# contextdim = (0.5,0.5,0.5,0.5);
# (x1, y1) = (0.5,0.0);
# c0 = context();
# c1 = context(contextdim...;units=UnitBox(UB...));
# (x0,y0) = ctparent((x1,y1), contextdim;UnitBox=UB)
# compose(c0,
#     (c0, circle(x0,y0,0.09), stroke("cyan")),
#     (c1, circle(x1,y1,0.1), stroke("red"))
# 	, fill(nothing)
#     )

@testset "Check if `ctparent` works fine:" begin

	UB = (0,0,2,2);
	contextdim = (0.0,0.5,0.5,0.5);
	(x1, y1) = (1,1);
	@test ctparent((x1,y1), contextdim;UnitBox=UB) == (0.25, 0.75)

	UB = (0,0,3,2);
	contextdim = (0.0,0.2,0.4,0.1);
	(x1, y1) = (1.5,1.0);
	@test ctparent((x1,y1), contextdim;UnitBox=UB) == (0.2, 0.25)

	UB = (0,0,1,1);
	contextdim = (0.5,0.5,0.5,0.5);
	(x1, y1) = (0.5,0.0);
	@test ctparent((x1,y1), contextdim;UnitBox=UB) == (0.75, 0.5)
	# @test ctparent((x1,y1), contextdim;UnitBox=UB) ==
end


text:
import Base: convert, haskey, setindex!, getindex
export slider, togglebutton, button,
       checkbox, textbox, textarea,
       radiobuttons, dropdown, select,
       togglebuttons, html, latex,
       progress, widget

const Empty = VERSION < v"0.4.0-dev" ? Nothing : Void

### Input widgets

########################## Slider ############################

type Slider{T<:Number} <: InputWidget{T}
    signal::Signal{T}
    label::AbstractString
    value::T
    range::Range{T}
    readout_format::AbstractString
    continuous_update::Bool
end

# differs from median(r) in that it always returns an element of the range
medianelement(r::Range) = r[(1+length(r))>>1]

slider(args...) = Slider(args...)
"""
    slider(range; value, signal, label="", continuous_update=true)

Create a slider widget with the specified `range`. Optionally specify
the starting `value` (defaults to the median of `range`), provide the
(Reactive.jl) `signal` coupled to this slider, and/or specify a string
`label` for the widget.
"""
slider{T}(range::Range{T};
          value=medianelement(range),
          signal::Signal{T}=Signal(value),
          label="",
          readout_format=T <: Integer ? "d" : ".3f",
          continuous_update=true) =
              Slider(signal, label, value, range, readout_format, continuous_update)

######################### Checkbox ###########################

type Checkbox <: InputWidget{Bool}
    signal::Signal{Bool}
    label::AbstractString
    value::Bool
end

checkbox(args...) = Checkbox(args...)

"""
    checkbox(value=false; label="", signal)

Provide a checkbox with the specified starting (boolean)
`value`. Optional provide a `label` for this widget and/or the
(Reactive.jl) `signal` coupled to this widget.
"""
checkbox(value::Bool; signal=Signal(value), label="") =
    Checkbox(signal, label, value)
checkbox(; label="", value=false, signal=Signal(value)) =
    Checkbox(signal, label, value)

###################### ToggleButton ########################

type ToggleButton <: InputWidget{Bool}
    signal::Signal{Bool}
    label::AbstractString
    value::Bool
end

togglebutton(args...) = ToggleButton(args...)

togglebutton(; label="", value=false, signal=Signal(value)) =
    ToggleButton(signal, label, value)

"""
    togglebutton(label=""; value=false, signal)

Create a toggle button. Optionally specify the `label`, the initial
state (`value=false` is off, `value=true` is on), and/or provide the
(Reactive.jl) `signal` coupled to this button.
"""
togglebutton(label; kwargs...) =
    togglebutton(label=label; kwargs...)

######################### Button ###########################

type Button{T} <: InputWidget{T}
    signal::Signal{T}
    label::AbstractString
    value::T
end

button(; value=nothing, label="", signal=Signal(value)) =
    Button(signal, label, value)

"""
    button(label; value=nothing, signal)

Create a push button. Optionally specify the `label`, the `value`
emitted when then button is clicked, and/or the (Reactive.jl) `signal`
coupled to this button.
"""
button(label; kwargs...) =
    button(label=label; kwargs...)

######################## Textbox ###########################

type Textbox{T} <: InputWidget{T}
    signal::Signal{T}
    label::AbstractString
    @compat range::Union{Empty, Range}
    value::T
end

function empty(t::Type)
    if is(t, Number) zero(t)
    elseif is(t, AbstractString) ""
    end
end

function Textbox(; label="",
                 value=utf8(""),
                 # Allow unicode characters even if initiated with ASCII
                 typ=typeof(value),
                 range=nothing,
                 signal=Signal(typ, value))
    if isa(value, AbstractString) && range != nothing
        throw(ArgumentError(
               "You cannot set a range on a string textbox"
             ))
    end
    Textbox{typ}(signal, label, range, value)
end

textbox(;kwargs...) = Textbox(;kwargs...)

"""
    textbox(value=""; label="", typ=typeof(value), range=nothing, signal)

Create a box for entering text. `value` is the starting value; if you
don't want to provide an initial value, you can constrain the type
with `typ`. Optionally provide a `label`, specify the allowed range
(e.g., `-10.0:10.0`) for numeric entries, and/or provide the
(Reactive.jl) `signal` coupled to this text box.
"""
textbox(val; kwargs...) =
    Textbox(value=val; kwargs...)
textbox(val::AbstractString; kwargs...) =
    Textbox(value=utf8(val); kwargs...)

parse_msg{T<:Number}(w::Textbox{T}, val::AbstractString) = parse_msg(w, parse(T, val))
function parse_msg{T<:Number}(w::Textbox{T}, val::Number)
    v = convert(T, val)
    if isa(w.range, Range)
        # force value to stay in range
        v = max(first(w.range),
                min(last(w.range), v))
    end
    v
end

######################### Textarea ###########################

type Textarea{AbstractString} <: InputWidget{AbstractString}
    signal::Signal{AbstractString}
    label::AbstractString
    value::AbstractString
end

textarea(args...) = Textarea(args...)

textarea(; label="",
         value="",
         signal=Signal(value)) =
    Textarea(signal, label, value)

"""
    textarea(value=""; label="", signal)

Creates an extended text-entry area. Optionally provide a `label`
and/or the (Reactive.jl) `signal` associated with this widget. The
`signal` updates when you type.
"""
textarea(val; kwargs...) =
    textarea(value=val; kwargs...)

##################### SelectionWidgets ######################

immutable OptionDict
    keys::Vector
    dict::Dict
end

Base.getindex(x::OptionDict, y) = getindex(x.dict, y)
Base.haskey(x::OptionDict, y) = haskey(x.dict, y)
Base.keys(x::OptionDict) = x.keys
Base.values(x::OptionDict) = [x.dict[k] for k in keys(x)]
function Base.setindex!(x::OptionDict, v, k)
    if !haskey(x.dict, k)
        push!(x.keys, k)
    end
    x.dict[k] = v
    v
end
type Options{view, T} <: InputWidget{T}
    signal::Signal
    label::AbstractString
    value::T
    value_label::AbstractString
    options::OptionDict
    icons::AbstractArray
    tooltips::AbstractArray
end

Options(view::Symbol, options::OptionDict;
        label = "",
        value_label=first(options.keys),
        value=options[value_label],
        icons=[],
        tooltips=[],
        typ=typeof(value),
        signal=Signal(value)) =
            Options{view, typ}(signal, label, value, value_label, options, icons, tooltips)

addoption(opts, v::NTuple{2}) = opts[string(v[1])] = v[2]
addoption(opts, v) = opts[string(v)] = v
function Options(view::Symbol,
                    options::AbstractArray;
                    kwargs...)
    opts = OptionDict(Any[], Dict())
    for v in options
        addoption(opts, v)
    end
    Options(view, opts; kwargs...)
end

function Options(view::Symbol,
                    options::Associative;
                    kwargs...)
    opts = OptionDict(Any[], Dict())
    for (k, v) in options
        opts[string(k)] = v
    end
    Options(view, opts; kwargs...)
end

"""
    dropdown(choices; label="", value, typ, icons, tooltips, signal)

Create a "dropdown" widget. `choices` can be a vector of
options. Optionally specify the starting `value` (defaults to the
first choice), the `typ` of elements in `choices`, supply custom
`icons`, provide `tooltips`, and/or specify the (Reactive.jl) `signal`
coupled to this widget.

# Examples

    a = dropdown(["one", "two", "three"])

To link a callback to the dropdown, use

    f = dropdown(["turn red"=>colorize_red, "turn green"=>colorize_green])
    map(g->g(image), signal(f))
"""
dropdown(opts; kwargs...) =
    Options(:Dropdown, opts; kwargs...)

"""
radiobuttons: see the help for `dropdown`
"""
radiobuttons(opts; kwargs...) =
    Options(:RadioButtons, opts; kwargs...)

select(opts; kwargs...) =
    Options(:Select, opts; kwargs...)

"""
togglebuttons: see the help for `dropdown`
"""
togglebuttons(opts; kwargs...) =
    Options(:ToggleButtons, opts; kwargs...)

### Output Widgets

export Latex, Progress

Base.@deprecate html(value; label="")  HTML(value)

type Latex <: Widget
    label::AbstractString
    value::AbstractString
end
latex(label, value::AbstractString) = Latex(label, value)
latex(value::AbstractString; label="") = Latex(label, value)
latex(value; label="") = Latex(label, mimewritable("application/x-latex", value) ? stringmime("application/x-latex", value) : stringmime("text/latex", value))

## # assume we already have Latex
## writemime(io::IO, m::MIME{symbol("application/x-latex")}, l::Latex) =
##     write(io, l.value)

type Progress <: Widget
    label::AbstractString
    value::Int
    range::Range
end

progress(args...) = Progress(args...)
progress(;label="", value=0, range=0:100) =
    Progress(label, value, range)

# Make a widget out of a domain
widget(x::Signal, label="") = x
widget(x::Widget, label="") = x
widget(x::Range, label="") = slider(x, label=label)
widget(x::AbstractVector, label="") = togglebuttons(x, label=label)
widget(x::Associative, label="") = togglebuttons(x, label=label)
widget(x::Bool, label="") = checkbox(x, label=label)
widget(x::AbstractString, label="") = textbox(x, label=label, typ=AbstractString)
widget{T <: Number}(x::T, label="") = textbox(typ=T, value=x, label=label)


text:
export Domain,GlobalData,updateStates!,updateTimeDependentEssentialBoundaryCondition!,
    setConstantNodalForces!, setGeometryPoints!, setConstantDirichletBoundary!, getExternalForce,
    commitHistory, getBodyForce, getEdgeForce, getNGauss, getDofs, getDStrain, getCoords, getState,
    getGaussPoints, getElems, getStressHistory, getStrainHistory, getStateHistory

@doc raw"""
    GlobalData

Store data for finite element updates. Assume the problem has n freedoms,

- `state`: a vector of length $n$. Displacement array at the **current** time, only for **active** freedoms.
   The ordering is based on the equation number.
- `Dstate`: a vector of length $n$. Displacement array at the **previous** time.
- `velo`: a vector of length $n$. Velocity array at the **current** time.
- `acce`: a vector of length $n$. Acceleration array at the **current** time.
- `time`: float, current time.
- `M`: a matrix of size $n\times n$ spares mass matrix
- `Mlumped`: a vector of length $n$ lumped mass array
- `MID`: Float64[n, nd1] off-diagonal part of the mass matrix, between the active freedoms and the time-dependent Dirichlet freedoms, assume there are nd time-dependent Dirichlet freedoms
- `EBC_func`: displacement $d$, velocity $v$, and acceleration $a$ from time-dependent Dirichlet boundary conditions

$$d, v, a = \text{EBC\_func}(\text{time})$$

The length of each output is the same as number of "-2" in `EBC` array. The ordering is direction major, i.e., $u_1, u_3, \ldots, v_1, v_3, \ldots$

- `FBC_func`: time-dependent load boundary condition.

$$f = \text{FBC\_func}(\text{time})$$

Here $f$ is a vector. Its length is the same as number of "-2" in `FBC` array. The ordering is direction major, i.e., $u_1, u_3, \ldots, v_1, v_3, \ldots$

- `Body_func`: time-dependent/independent body force function.

$$f = \text{Body\_func}(x_{\text{array}}, y_{\text{array}}, \text{time})$$

Here $f$ is a vector or a matrix (second dimension is 2) depending on the dimension of state variables.
The output is a $N\times n_{\text{dim}}$ matrix, where $N$ is the length of $x_{\text{array}}$ or $y_{\text{array}}$, and $n_{\text{dim}}$ is the dimension of the problem (1 or 2).

- `Edge_func`: time-dependent/independent traction load.

$$f = \text{Edge\_func}(x_{\text{array}}, y_{\text{array}}, \text{time}, \text{id})$$

Here $f$ is a vector. Its length is the same as the length of $x_{\text{array}}$ or $y_{\text{array}}$.

"""
mutable struct GlobalData
    state::Array{Float64}    #u
    Dstate::Array{Float64}   #uk
    velo::Array{Float64}     #∂u
    acce::Array{Float64}     #∂∂u
    time::Float64
    M::Union{SparseMatrixCSC{Float64,Int64},Array{Float64}}
    Mlumped::Array{Float64}
    MID::Array{Float64}

    EBC_func::Union{Function,Nothing}  #time dependent Dirichlet boundary condition f(t)
    FBC_func::Union{Function,Nothing}  #time dependent nodal force  boundary condition f(t)
    Body_func::Union{Function,Nothing} #body force function f(x_array, y_array, t)
    Edge_func::Union{Function,Nothing} #edge traction function f(x_array, y_array, t, id)

end

function Base.show(io::IO, z::GlobalData)
    yes = "✔"
    no = "✘"
print(io, """GlobalData with $(length(z.state)) active DOFs, time=$(z.time)
Mass matrix ... $(length(z.M)==0 ? no : yes)
EBC_func    ... $(isnothing(z.EBC_func) ? no : yes)
FBC_func    ... $(isnothing(z.FBC_func) ? no : yes)
Body_func   ... $(isnothing(z.Body_func) ? no : yes)
Edge_func   ... $(isnothing(z.Edge_func) ? no : yes)
""")
end

@doc raw"""
    GlobalData(state::Union{Array{Float64,1},Missing},
            Dstate::Union{Array{Float64,1},Missing},
            velo::Union{Array{Float64,1},Missing},
            acce::Union{Array{Float64,1},Missing},
            neqs::Int64,
            EBC_func::Union{Function, Nothing}=nothing, FBC_func::Union{Function, Nothing}=nothing,
            Body_func::Union{Function,Nothing}=nothing, Edge_func::Union{Function,Nothing}=nothing)

The size of `state`, `Dstate`, `velo`, `acce` must be `neqs`, i.e., the active DOFs.
If they are missing, they are treated as zeros.

`state`, `Dstate`, `velo`, `acce` can be missing, in which case they are interpreted as zeros.
"""
function GlobalData(state::Union{Array{Float64,1},Missing},Dstate::Union{Array{Float64,1},Missing},
        velo::Union{Array{Float64,1},Missing},acce::Union{Array{Float64,1},Missing},
        neqs::Int64,
        EBC_func::Union{Function, Nothing}=nothing, FBC_func::Union{Function, Nothing}=nothing,
        Body_func::Union{Function,Nothing}=nothing, Edge_func::Union{Function,Nothing}=nothing)
    time = 0.0
    M = Float64[]
    Mlumped = Float64[]
    MID = Float64[]
    state = coalesce(state, zeros(neqs))
    Dstate = coalesce(Dstate, zeros(neqs))
    velo = coalesce(velo, zeros(neqs))
    acce = coalesce(acce, zeros(neqs))
    GlobalData(state, Dstate, velo, acce, time, M, Mlumped, MID, EBC_func, FBC_func, Body_func, Edge_func)
end

@doc raw"""
    Domain

Date structure for the computatational domain.

- `nnodes`: Int64, number of nodes (each quadratical quad element has 9 nodes)
- `nodes`: Float64[nnodes, ndims], coordinate array of all nodes
- `neles`: number of elements
- `elements`: a list of `neles` element arrays, each element is a struct
- `ndims`: Int64, dimension of the problem space
- `state`: a matrix of size `nnodes×ndims`. **Current** displacement of all nodal freedoms, `state[1:nnodes]` are for the first direction.
- `Dstate`: `nnodes×ndims`. **Previous** displacement of all nodal freedoms, `Dstate[1:nnodes]` are for the first direction.
- `LM`:  `Int64[neles][ndims]`, LM(e,d) is the global equation number (active freedom number) of element e's d th freedom,

         ∘ -1 means fixed (time-independent) Dirichlet

         ∘ -2 means time-dependent Dirichlet

         ∘ >0 means the global equation number

- `DOF`: a matrix of size `neles×ndims`, DOF(e,d) is the global freedom (including both active and inactive DOFs) number of element e's d th freedom.
- `ID`:  a matrix of size `nnodes×ndims`. `ID(n,d)` is the equation number (active freedom number) of node n's $d$-th freedom,

         ∘ -1 means fixed (time-independent) Dirichlet

         ∘ -2 means time-dependent Dirichlet

         ∘ >0 means the global equation number

- `neqs`:  Int64,  number of equations, a.k.a., active freedoms
- `eq_to_dof`:  an integer vector of length `neqs`, map from to equation number (active freedom number) to the freedom number (Int64[1:nnodes] are for the first direction)
- `dof_to_eq`:  a bolean array of size `nnodes×ndims`, map from freedom number(Int64[1:nnodes] are for the first direction) to booleans (active freedoms(equation number) are true)
- `EBC`:  Int64[nnodes, ndims], EBC[n,d] is the displacement boundary condition of node n's dth freedom,
           -1 means fixed(time-independent) Dirichlet boundary nodes
           -2 means time-dependent Dirichlet boundary nodes
- `g`:  Float64[nnodes, ndims], values for fixed(time-independent) Dirichlet boundary conditions of node n's dth freedom,
- `FBC`: Int64[nnodes, ndims], `FBC[n,d]`` is the force load boundary condition of node n's dth freedom,
           -1 means constant(time-independent) force load boundary nodes
           -2 means time-dependent force load boundary nodes
- `fext`:  Float64[neqs], constant (time-independent) nodal forces on these freedoms
- `Edge_Traction_Data`: `n × 3` integer matrix for natural boundary conditions.

1. `Edge_Traction_Data[i,1]` is the element id,

2. `Edge_Traction_Data[i,2]` is the local edge id in the element, where the force is exterted (should be on the boundary, but not required)

3. `Edge_Traction_Data[i,3]` is the force id, which should be consistent with the last component of the `Edge_func` in the `Globdat`

- `time`: Float64, current time
- `npoints`: Int64, number of points (each quadratical quad element has 4 points, npoints==nnodes, when porder==1)
- `node_to_point`: Int64[nnodes]:map from node number to point point, -1 means the node is not a geometry point
"""
mutable struct Domain
    nnodes::Int64
    nodes::Array{Float64}
    neles::Int64
    elements::Array
    ndims::Int64
    state::Array{Float64}
    Dstate::Array{Float64}
    LM::Array{Array{Int64}}
    DOF::Array{Array{Int64}}
    ID::Array{Int64}
    neqs::Int64
    eq_to_dof::Array{Int64}
    dof_to_eq::Array{Bool}
    EBC::Array{Int64}  # Dirichlet boundary condition
    g::Array{Float64}  # Value for Dirichlet boundary condition
    FBC::Array{Int64}  # Nodal force boundary condition
    fext::Array{Float64}  # Value for Nodal force boundary condition
    edge_traction_data::Array{Int64} # traction force location, [element id, local edge id, force id]
    time::Float64

    npoints::Int64     # number of mesh points(the same as nodes, when porder==1)
    node_to_point::Array{Int64} # map from node to point, -1 means the node is not a geometry point

    ii_stiff::Array{Int64}
    jj_stiff::Array{Int64}
    vv_stiff_ele_indptr::Array{Int64}
    vv_stiff::Array{Float64}

    ii_dfint_dstress::Array{Int64}
    jj_dfint_dstress::Array{Int64}
    vv_dfint_dstress_ele_indptr::Array{Int64}
    vv_dfint_dstress::Array{Float64}

    ii_dstrain_dstate::Array{Int64}
    jj_dstrain_dstate::Array{Int64}
    vv_dstrain_dstate_ele_indptr::Array{Int64}
    vv_dstrain_dstate::Array{Float64}

    history::Dict{String, Array{Array{Float64}}}

end

function Base.show(io::IO, z::Domain)
    yes = "✔"
    no = "✘"
    print(io, """Domain with $(length(z.elements)) elements, $(z.nnodes) nodes and $(z.neqs) active equations
    edge_traction_data ... $(size(z.edge_traction_data, 1)==0 ? no : yes)
    """)
end

function Base.:copy(g::GlobalData)
    names = propertynames(g)
    args = [copy(getproperty(g, n)) for n in names]
    GlobalData(args...)
end

function Base.:copy(g::Domain)
    names = propertynames(g)
    args = [copy(getproperty(g, n)) for n in names]
    Domain(args...)
end

@doc raw"""
    Domain(nodes::Array{Float64}, elements::Array, ndims::Int64 = 2,
    EBC::Union{Missing, Array{Int64}} = missing, g::Union{Missing, Array{Float64}} = missing, FBC::Union{Missing, Array{Int64}} = missing,
    f::Union{Missing, Array{Float64}} = missing, edge_traction_data::Array{Int64,2}=zeros(Int64,0,3))

Creating a finite element domain.

- `nodes`: coordinate array of all nodes, a `nnodes × 2` matrix
- `elements`: element array. Each element is a material struct, e.g., [`PlaneStrain`](@ref).
- `ndims`: dimension of the problem space. For 2D problems, ndims = 2.
- `EBC`:  `nnodes × ndims` integer matrix for essential boundary conditions
    `EBC[n,d]`` is the displacement boundary condition of node `n`'s $d$-th freedom,

    ∘ -1: fixed (time-independent) Dirichlet boundary nodes

    ∘ -2: time-dependent Dirichlet boundary nodes

- `g`:  `nnodes × ndims` double matrix, values for fixed (time-independent) Dirichlet boundary conditions of node `n`'s $d$-th freedom,
- `FBC`: `nnodes × ndims` integer matrix for nodal force boundary conditions.
FBC[n,d] is the force load boundary condition of node n's dth freedom,

∘ -1 means constant(time-independent) force load boundary nodes

∘ -2 means time-dependent force load boundary nodes

- `f`:  `nnodes × ndims` double matrix, values for constant (time-independent) force load boundary conditions of node n's $d$-th freedom,

- `Edge_Traction_Data`: `n × 3` integer matrix for natural boundary conditions.
`Edge_Traction_Data[i,1]` is the element id,
`Edge_Traction_Data[i,2]` is the local edge id in the element, where the force is exterted (should be on the boundary, but not required)
`Edge_Traction_Data[i,3]` is the force id, which should be consistent with the last component of the Edge_func in the Globdat

For time-dependent boundary conditions (`EBC` or `FBC` entries are -2), the corresponding `f` or `g` entries are not used.
"""
function Domain(nodes::Array{Float64}, elements::Array, ndims::Int64 = 2,
    EBC::Union{Missing, Array{Int64}} = missing, g::Union{Missing, Array{Float64}} = missing, FBC::Union{Missing, Array{Int64}} = missing,
    f::Union{Missing, Array{Float64}} = missing, edge_traction_data::Array{Int64,2}=zeros(Int64,0,3))

    nnodes = size(nodes,1)
    neles = size(elements,1)
    state = zeros(nnodes * ndims)
    Dstate = zeros(nnodes * ndims)
    LM = Array{Int64}[]
    DOF = Array{Int64}[]
    ID = Int64[]
    neqs = 0
    eq_to_dof = Int64[]
    dof_to_eq = zeros(Bool, nnodes * ndims)
    fext = Float64[]

    npoints = -1
    node_to_point = Int64[]

    history = Dict("state"=>Array{Float64}[], "acc"=>Array{Float64}[], "fint"=>Array{Float64}[],
                "fext"=>Array{Float64}[], "strain"=>[], "stress"=>[], "time"=>Array{Float64}[])

    EBC = coalesce(EBC, zeros(Int64, nnodes, ndims))
    FBC = coalesce(FBC, zeros(Int64, nnodes, ndims))
    g = coalesce(g, zeros(Float64, nnodes, ndims))
    f = coalesce(f, zeros(Float64, nnodes, ndims))

    domain = Domain(nnodes, nodes, neles, elements, ndims, state, Dstate,
    LM, DOF, ID, neqs, eq_to_dof, dof_to_eq,
    EBC, g, FBC, fext, edge_traction_data, 0.0, npoints, node_to_point,
    Int64[], Int64[], Int64[], Float64[],
    Int64[], Int64[], Int64[], Float64[],
    Int64[], Int64[], Int64[], Float64[], history)

    #set fixed(time-independent) Dirichlet boundary conditions
    setConstantDirichletBoundary!(domain, EBC, g)
    #set constant(time-independent) force load boundary conditions
    setConstantNodalForces!(domain, FBC, f)
    assembleSparseMatrixPattern!(domain)

    domain
end

@doc """
    commitHistory(domain::Domain)

Update current step strain and stress in the history map of the `domain`.
This is essential for visualization and time dependent constitutive relations.
"""
function commitHistory(domain::Domain)
    for e in domain.elements
        commitHistory(e)
    end

    if options.save_history>=1
        # 1D, nstrain=1; 2D, nstrain=3
        eledim = domain.elements[1].eledim
        nstrain = div((eledim + 1)*eledim, 2)
        ngp = domain.neles * length(domain.elements[1].weights)
        if nstrain==1
            strain = zeros(ngp)
            stress = zeros(ngp)
            k = 1
            for e in domain.elements
                for igp in e.mat
                    strain[k] = igp.ε0
                    stress[k] = igp.σ0
                    k += 1
                end
            end
        else
            strain = zeros(ngp, nstrain)
            stress = zeros(ngp, nstrain)
            k = 1
            for e in domain.elements
                for igp in e.mat
                    strain[k,:] = igp.ε0
                    stress[k,:] = igp.σ0
                    k += 1
                end
            end
        end
        push!(domain.history["strain"], strain)
        push!(domain.history["stress"], stress)
    end
end

@doc """
    setConstantDirichletBoundary!(self::Domain, EBC::Array{Int64}, g::Array{Float64})

Bookkeepings for time-independent Dirichlet boundary conditions. Only called once in the constructor of `domain`.
It updates the fixed (time-independent Dirichlet boundary) state entries and builds both LM and DOF arrays.

- `self`: Domain
- `EBC`:  Int64[nnodes, ndims], EBC[n,d] is the displacement boundary condition of node n's dth freedom,

  ∘ -1 means fixed(time-independent) Dirichlet boundary nodes

  ∘ -2 means time-dependent Dirichlet boundary nodes

- `g`:  Float64[nnodes, ndims], values for fixed (time-independent) Dirichlet boundary conditions of node n's dth freedom,
""" ->
function setConstantDirichletBoundary!(domain::Domain, EBC::Array{Int64}, g::Array{Float64})

    # ID(n,d) is the global equation number of node n's dth freedom,
    # -1 means fixed (time-independent) Dirichlet
    # -2 means time-dependent Dirichlet

    nnodes, ndims = domain.nnodes, domain.ndims
    neles, elements = domain.neles, domain.elements
    #ID = zeros(Int64, nnodes, ndims) .- 1

    ID = copy(EBC)

    eq_to_dof, dof_to_eq = Int64[], zeros(Bool, nnodes * ndims)
    neqs = 0
    for idof = 1:ndims
      for inode = 1:nnodes
          if (EBC[inode, idof] == 0)
              neqs += 1
              ID[inode, idof] = neqs
              push!(eq_to_dof,inode + (idof-1)*nnodes)
              dof_to_eq[(idof - 1)*nnodes + inode] = true
          elseif (EBC[inode, idof] == -1)
              #update state fixed (time-independent) Dirichlet boundary conditions
              domain.state[inode + (idof-1)*nnodes] = g[inode, idof]
          end
      end
    end

    domain.ID, domain.neqs, domain.eq_to_dof, domain.dof_to_eq = ID, neqs, eq_to_dof, dof_to_eq

    # LM(e,d) is the global equation number of element e's d th freedom
    LM = Array{Array{Int64}}(undef, neles)
    for iele = 1:neles
      el_nodes = getNodes(elements[iele])
      ieqns = ID[el_nodes, :][:]
      LM[iele] = ieqns
    end
    domain.LM = LM

    # DOF(e,d) is the global dof number of element e's d th freedom

    DOF = Array{Array{Int64}}(undef, neles)
    for iele = 1:neles
      el_nodes = getNodes(elements[iele])
      if domain.ndims==1
        DOF[iele] = el_nodes
      elseif domain.ndims==2
        DOF[iele] = [el_nodes;[idof + nnodes for idof in el_nodes]]
      else
        error("NOT IMPLEMENTED YET")
      end

    end
    domain.DOF = DOF

end

@doc """

Bookkeepings for time-independent Nodal force boundary conditions. Only called once in the constructor of `domain`.
It updates the fixed (time-independent Nodal forces) state entries and builds both LM and DOF arrays.

- `self`: Domain
- `FBC`:  Int64[nnodes, ndims], FBC[n,d] is the displacement boundary condition of node n's dth freedom,

    ∘ -1 means fixed (time-independent) Nodal force freedoms

    ∘ -2 means time-dependent Nodal force freedoms

- `f`:  Float64[nnodes, ndims], values for fixed (time-independent) Neumann boundary conditions of node n's dth freedom,

#The name is misleading
"""
function setConstantNodalForces!(self::Domain, FBC::Array{Int64}, f::Array{Float64})

    fext = zeros(Float64, self.neqs)
    # ID(n,d) is the global equation number of node n's dth freedom, -1 means no freedom

    nnodes, ndims, ID = self.nnodes, self.ndims, self.ID
    for idof = 1:ndims
      for inode = 1:nnodes
          if (FBC[inode, idof] == -1)
              if  ID[inode, idof] <= 0
                error("Node $inode is both a Dirichlet node and a force node.")
              end
              fext[ID[inode, idof]] += f[inode, idof]
          end
        end
    end
    self.fext = fext
end

@doc """
    updateStates!(domain::Domain, globaldat::GlobalData)
    update time-dependent Dirichlet boundary condition to globaldat.time

Update `state` and `Dstate` in `domain`. This includes

- Copy state variable values for active DOFs from `globaldat`

- Set time-dependent essential boundary conditions using `globaldat.EBC`

The time-independent boundary conditions are inherented from last time step.
"""
function updateStates!(domain::Domain, globaldat::GlobalData)

    domain.state[domain.eq_to_dof] = globaldat.state[:]

    domain.time = globaldat.time
    if options.save_history>=1
        push!(domain.history["state"], copy(domain.state))
        push!(domain.history["acc"], copy(globaldat.acce))
    end

    updateTimeDependentEssentialBoundaryCondition!(domain, globaldat)

    domain.Dstate[:] = domain.state
end

@doc """
    updateTimeDependentEssentialBoundaryCondition!(domain::Domain, globaldat::GlobalData)

If there exists time-dependent Dirichlet boundary conditions, `updateTimeDependentEssentialBoundaryCondition!` must be called to update
the boundaries in `domain`. This function is called by [`updateStates!`](@ref)

This function updates `state` data in `domain`.
"""
function updateTimeDependentEssentialBoundaryCondition!(domain::Domain, globaldat::GlobalData)
    if globaldat.EBC_func ≠ nothing
        disp, _, _ = globaldat.EBC_func(globaldat.time) # user defined time-dependent boundary
        dof_id = 0

        #update state of all nodes
        for idof = 1:domain.ndims
            for inode = 1:domain.nnodes
                if (domain.EBC[inode, idof] == -2)
                    dof_id += 1
                    domain.state[inode + (idof-1)*domain.nnodes] = disp[dof_id]
                end
            end
        end
    end

end

@doc """
    getExternalForce(self::Domain, globaldat::GlobalData, fext::Union{Missing,Array{Float64}}=missing)

Computes external force vector at `globaldat.time`,
This includes all the body force, external load, and internal force caused by acceleration.
"""
function getExternalForce(domain::Domain, globaldat::GlobalData)
    #Update time-independent nodal force
    fext = copy(domain.fext)

    #Update time-dependent nodal force
    if !isnothing(globaldat.FBC_func)
        ID = domain.ID
        nodal_force = globaldat.FBC_func(globaldat.time) # user defined time-dependent boundary
        # @info nodal_force
        dof_id = 0
        #update fext for active nodes (length of neqs)
        for idof = 1:domain.ndims
            for inode = 1:domain.nnodes
                if (domain.FBC[inode, idof] == -2)
                    dof_id += 1
                    @assert ID[inode, idof] > 0
                    fext[ID[inode, idof]] += nodal_force[dof_id]
                end
            end
        end
    end

    #Update the acceleration effect from the time-dependent Dirichlet boundary condition
    if !isnothing(globaldat.EBC_func)
        MID = globaldat.MID
        _, _, acce = globaldat.EBC_func(globaldat.time)
        fext -= MID * acce
    end

    #Update time-dependent and time-independent body force
    fbody = getBodyForce(domain, globaldat, globaldat.time)

    #Update time-dependent and time-independent edge traction force
    fedge = getEdgeForce(domain, globaldat, globaldat.time)

    fext + fbody + fedge
end

@doc raw"""
    getBodyForce(domain::Domain, globdat::GlobalData, time::Float64)

Computes the body force vector $F_\mathrm{body}$ of length `neqs`
- `globdat`: GlobalData
- `domain`: Domain, finite element domain, for data structure
- `Δt`:  Float64, current time step size
"""
function getBodyForce(domain::Domain, globdat::GlobalData, time::Float64)

    Fbody = zeros(Float64, domain.neqs)
    neles = domain.neles

    if isnothing(globdat.Body_func)
        return Fbody
    end

    # Loop over the elements in the elementGroup
    for iele  = 1:neles
        element = domain.elements[iele]

        gauss_pts = getGaussPoints(element)
        fvalue = globdat.Body_func(gauss_pts[:,1], gauss_pts[:,2], time)

        fbody = getBodyForce(element, fvalue)

      # Assemble in the global array
        el_eqns = getEqns(domain, iele)
        el_eqns_active = (el_eqns .>= 1)
        Fbody[el_eqns[el_eqns_active]] += fbody[el_eqns_active]
    end

    return Fbody
end

@doc raw"""
    getEdgeForce(domain::Domain, globdat::GlobalData, time::Float64)

Computes the edge force vector $F_\mathrm{edge}$ defined in `domain.edge_traction_data`
- `globdat`: GlobalData
- `domain`: Domain, finite element domain, for data structure
- `time`:  Float64, current time step size
"""
function getEdgeForce(domain::Domain, globdat::GlobalData, time::Float64)
    Fedge = zeros(Float64, domain.neqs)
    neles = domain.neles

    if (isnothing(globdat.Edge_func) && size(domain.edge_traction_data, 1)!=0) ||
        (!isnothing(globdat.Edge_func) && size(domain.edge_traction_data, 1)==0)
        @warn("`GlobalData` does not have (has) `Edge_func` but `Domain` has (does not have) `edge_traction_data`")
    end

    if isnothing(globdat.Edge_func)
        return Fedge
    end

    # Loop over the elements in the elementGroup

    for i = 1:size(domain.edge_traction_data)[1]

        iele, iedge, ifunc = domain.edge_traction_data[i, :]

        element = domain.elements[iele]

        gauss_pts = getEdgeGaussPoints(element, iedge)

        fvalue = globdat.Edge_func(gauss_pts[:,1], gauss_pts[:,2], time, ifunc)

        fedge = getEdgeForce(element, iedge, fvalue)

      # Assemble in the global array
        el_eqns = getEqns(domain, iele)
        el_eqns_active = (el_eqns .>= 1)
        Fedge[el_eqns[el_eqns_active]] += fedge[el_eqns_active]
    end

    return Fedge
end

@doc """
    getCoords(domain::Domain, el_nodes::Array{Int64})

Get the coordinates of several nodes (possibly in one element)
- `domain`: Domain
- `el_nodes`: Int64[n], node array

Return: Float64[n, ndims], the coordinates of these nodes
"""
function getCoords(domain::Domain, el_nodes::Array{Int64})
    return domain.nodes[el_nodes, :]
end

@doc """
    getDofs(domain::Domain, iele::Int64)

Get the global freedom numbers of the element
- `domain`: Domain
- `iele`: Int64, element number

Return: Int64[], the global freedom numbers of the element (ordering in local element ordering)
"""
function getDofs(domain::Domain, iele::Int64)
    return domain.DOF[iele]
end

@doc """
    getNGauss(domain::Domain)

Gets the total number of Gauss quadrature points.
"""
function getNGauss(domain::Domain)
    ng = 0
    for e in domain.elements
        ng += length(e.weights)
    end
    ng
end

@doc """
    getEqns(domain::Domain, iele::Int64)

Gets the equation numbers (active freedom numbers) of the element.
This excludes both the time-dependent and time-independent Dirichlet boundary conditions.
"""
function getEqns(domain::Domain, iele::Int64)
    return domain.LM[iele]
end

@doc """
    getState(domain::Domain, el_dofs::Array{Int64})

Get the displacements of several nodes (possibly in one element)
- `domain`: Domain
- `el_nodes`: Int64[n], node array

Return: Float64[n, ndims], the displacements of these nodes
"""
function getState(domain::Domain, el_dofs::Array{Int64})
    return domain.state[el_dofs]
end

@doc """
    getDstate(domain::Domain, el_dofs::Array{Int64})

Get the displacements of several nodes (possibly in one element) at the previous time step
- `domain`: Domain
- `el_nodes`: Int64[n], node array

Return: Float64[n, ndims], the displacements of these nodes at the previous time step
""" ->

function getDstate(domain::Domain, el_dofs::Array{Int64})
    return domain.Dstate[el_dofs]
end

@doc raw"""
    getGaussPoints(domain::Domain)

Returns all Gauss points as a $n_g\times 2$ matrix, where $n_g$ is the total number of Gauss points.
"""
function getGaussPoints(domain::Domain)
    v = []
    for e in domain.elements
        vg = getGaussPoints(e)
        push!(v, vg)
    end
    vcat(v...)
end

@doc raw"""
    getStrain(domain::Domain)

Computes the strain from the `domain` data. The output is $n_g\times 3$ matrix, where $n_g$ is the total number of Gauss points.
Each row is the strain tensor $\begin{bmatrix} \epsilon_{xx} & \epsilon_{yy} & \gamma_{xy} \end{bmatrix}$
"""
function getStrain(domain::Domain)
    E = Array{Array{Float64}}(undef, domain.neles)
    for (k,e) in enumerate(domain.elements)
        el_dof = getDofs(domain, k)
        state = getState(domain, el_dof)
        E[k], _ = getStrain(e, state)
    end
    vcat(E...)
end

@doc raw"""
    getDStrain(domain::Domain)

Computes the strain **at last time step** from the `domain` data. The output is $n_g\times 3$ matrix, where $n_g$ is the total number of Gauss points.
Each row is the strain tensor $\begin{bmatrix} \epsilon_{xx} & \epsilon_{yy} & \gamma_{xy} \end{bmatrix}$
"""
function getDStrain(domain::Domain)
    E = Array{Array{Float64}}(undef, domain.neles)
    for (k,e) in enumerate(domain.elements)
        el_dof = getDofs(domain, k)
        Dstate = getDstate(domain, el_dof)
        E[k], _ = getDStrain(e, Dstate)
    end
    vcat(E...)
end

@doc raw"""
    getStress(domain::Domain, Δt::Float64 = 0.0; save_trace::Bool = false)

Returns the stress based on `domain.state` and `domain.Dstate`. If `save_trace` is true,
the stress is also saved to `domain.stress`, which is useful for visualization.

The output is $n_g\times 3$ matrix, where $n_g$ is the total number of Gauss points.
Each row is the strain tensor $\begin{bmatrix} \sigma_{xx} & \sigma_{yy} & \sigma_{xy} \end{bmatrix}$
"""
function getStress(domain::Domain, Δt::Float64 = 0.0; save_trace::Bool = false)
    Ret = zeros(getNGauss(domain), 3)
    local ∂E∂u, E, DE
    l = 1
    for (k,elem) in enumerate(domain.elements)
        nnodes = length(elem.elnodes)
        el_dof = getDofs(domain,k)
        state = getState(domain, el_dof)
        Dstate = getDstate(domain, el_dof)
        u = state[1:nnodes]; v = state[nnodes+1:2*nnodes]
        Du = Dstate[1:nnodes]; Dv = Dstate[nnodes+1:2*nnodes]

        for k = 1:length(elem.weights)
            # #@show "Gaussian point ", k
            g1 = elem.dhdx[k][:,1]; g2 = elem.dhdx[k][:,2]

            ux = u'*g1; uy = u'*g2; vx = v'*g1; vy = v'*g2
            Dux = Du'*g1; Duy = Du'*g2; Dvx = Dv'*g1; Dvy = Dv'*g2

            if isa(elem, SmallStrainContinuum)
                ∂E∂u = [g1   zeros(nnodes)    g2;
                        zeros(nnodes)    g2   g1;]
                E = [ux; vy; uy+vx]
                DE = [Dux; Dvy; Duy+Dvx]
            elseif isa(elem, FiniteStrainContinuum)
                ∂E∂u = [g1+ux*g1 uy*g2    g2 + g2*ux+g1*uy;
                vx*g1    g2+vy*g2 g1 + g1*vy+g2*vx;]
                E = [ux+0.5*(ux*ux+vx*vx); vy+0.5*(uy*uy+vy*vy); uy+vx+ux*uy+vx*vy]
                DE = [Dux+0.5*(Dux*Dux+Dvx*Dvx); Dvy+0.5*(Duy*Duy+Dvy*Dvy); Duy+Dvx+Dux*Duy+Dvx*Dvy]
            else
                error("NOT IMPLEMENTED")
            end

            S, dS_dE = getStress(elem.mat[k], E, DE, Δt)
            Ret[l,:] = S
            l += 1

            if save_trace
                elem.stress[k] = S
            end
        end
    end
    Ret
end

@doc raw"""
    getElems(domain::Domain)

Returns the element connectivity matrix $n_e \times 4$. This function implicitly assumes that all elements are quadrilateral.
"""
function getElems(domain::Domain)
    elem = zeros(Int64, domain.neles, 4)
    for (k,e) in enumerate(domain.elements)
        elem[k,:] = e.elnodes
    end
    return elem
end

"""
    getStressHistory(domain::Domain)

Returns the stress history.
"""
function getStressHistory(domain::Domain)
    vcat([reshape(x, 1, size(x,1),size(x,2)) for x in domain.history["stress"]]...)
end

"""
    getStrainHistory(domain::Domain)

Returns the strain history.
"""
function getStrainHistory(domain::Domain)
    vcat([reshape(x, 1, size(x,1),size(x,2)) for x in domain.history["strain"]]...)
end

"""
    getStateHistory(domain::Domain)

Returns the state history.
"""
function getStateHistory(domain::Domain)
    vcat([reshape(x, 1, size(x,1)) for x in domain.history["state"]]...)
end

#------------------------------------------------------------------------------------
# Maybe useful in the future

@doc """
    Compute constant stiff, dfint_dstress, dstrain_dstate matrix patterns
    - 'self': Domain
"""
function assembleSparseMatrixPattern!(self::Domain)

    neles = self.neles
    eledim = self.elements[1].eledim
    nstrain = div((eledim + 1)*eledim, 2)
    ngps_per_elem = length(self.elements[1].weights)
    neqs = self.neqs

    ii_stiff = Int64[]; jj_stiff = Int64[]; vv_stiff_ele_indptr = ones(Int64, neles+1);
    ii_dfint_dstress = Int64[]; jj_dfint_dstress = Int64[]; vv_dfint_dstress_ele_indptr = ones(Int64, neles+1);
    ii_dstrain_dstate = Int64[]; jj_dstrain_dstate = Int64[]; vv_dstrain_dstate_ele_indptr = ones(Int64, neles+1);

    neles = self.neles

    # Loop over the elements in the elementGroup
    for iele  = 1:neles
      element = self.elements[iele]

      el_eqns = getEqns(self,iele)

      el_dofs = getDofs(self,iele)

      el_state  = getState(self, el_dofs)

      gp_ids = (iele-1)*ngps_per_elem+1 : iele*ngps_per_elem

      # Assemble in the global array
      el_eqns_active = el_eqns .>= 1
      el_eqns_active_idx = el_eqns[el_eqns_active]
      # K[el_eqns[el_eqns_active], el_eqns[el_eqns_active]] += stiff[el_eqns_active,el_eqns_active]

      el_eqns_active_idx = el_eqns[el_eqns_active]

      for j = 1:length(el_eqns_active_idx)
        for i = 1:length(el_eqns_active_idx)
          push!(ii_stiff, el_eqns_active_idx[i])
          push!(jj_stiff, el_eqns_active_idx[j])
          #push!(vv_stiff, stiff_active[i,j])
        end
      end
      vv_stiff_ele_indptr[iele+1] = vv_stiff_ele_indptr[iele] + length(el_eqns_active_idx)*length(el_eqns_active_idx)

      for j = 1:ngps_per_elem*nstrain
        for i = 1:length(el_eqns_active_idx)
          push!(ii_dfint_dstress, el_eqns_active_idx[i])
          push!(jj_dfint_dstress, (iele-1)*ngps_per_elem*nstrain+j)
          #push!(vv_dfint_dstress, dfint_dstress_active[i,j])
        end
      end
      vv_dfint_dstress_ele_indptr[iele+1] = vv_dfint_dstress_ele_indptr[iele] + ngps_per_elem*nstrain*length(el_eqns_active_idx)

      for j = 1:length(el_eqns_active_idx)
        for i = 1:ngps_per_elem*nstrain

          push!(ii_dstrain_dstate, (iele-1)*ngps_per_elem*nstrain+i)
          push!(jj_dstrain_dstate, el_eqns_active_idx[j])
          #push!(vv_dstrain_dstate, dstrain_dstate_active[i,j])
        end
      end
      vv_dstrain_dstate_ele_indptr[iele+1] = vv_dstrain_dstate_ele_indptr[iele] + ngps_per_elem*nstrain*length(el_eqns_active_idx)

    end

    self.ii_stiff = ii_stiff; self.jj_stiff = jj_stiff;
    self.vv_stiff_ele_indptr = vv_stiff_ele_indptr; self.vv_stiff = similar(ii_stiff)

    self.ii_dfint_dstress = ii_dfint_dstress; self.jj_dfint_dstress = jj_dfint_dstress;
    self.vv_dfint_dstress_ele_indptr = vv_dfint_dstress_ele_indptr; self.vv_dfint_dstress= similar(ii_dfint_dstress)

    self.ii_dstrain_dstate = ii_dstrain_dstate; self.jj_dstrain_dstate = jj_dstrain_dstate;
    self.vv_dstrain_dstate_ele_indptr = vv_dstrain_dstate_ele_indptr; self.vv_dstrain_dstate = similar(ii_dstrain_dstate)

  end

@doc """
    In the constructor
    Update the node_to_point map

    - `self`: Domain, finit element domain
    - 'npoints': Int64, number of points (each quadratical quad element has 4 points, npoints==nnodes, when porder==1)
    - 'node_to_point': Int64[nnodes]:map from node number to point point, -1 means the node is not a geometry point

"""
function setGeometryPoints!(self::Domain, npoints::Int64, node_to_point::Array{Int64})
    self.npoints = npoints
    self.node_to_point = node_to_point
end


text:
# can consolidate the remaining simulate method if we have default meanspecs for univariate/multivariate
# proper multivariate meanspec, include return prediction in predict
# implement correlations, covariances, residuals in terms of update!, and move them from DCC to multivariate
"""
	DOW29
Stock returns, in procent, from 03/19/2008 through 04/11/2019, for tickers
AAPL, IBM, XOM, KO, MSFT, INTC, MRK, PG, VZ, WBA, V, JNJ, PFE, CSCO,
TRV, WMT, MMM, UTX, UNH, NKE, HD, BA, AXP, MCD, CAT, GS, JPM, CVX, DIS.
"""
const DOW29 = readdlm(joinpath(dirname(pathof(ARCHModels)), "data", "dow29.csv"), ',')

"""
    MultivariateStandardizedDistribution{T, d} <: Distribution{Multivariate, Continuous}

Abstract supertype that multivariate standardized distributions inherit from.
"""
abstract type MultivariateStandardizedDistribution{T, d} <: Distribution{Multivariate, Continuous} end

"""
    MultivariateVolatilitySpec{T, d} <: VolatilitySpec{T}

Abstract supertype that multivariate volatility specifications inherit from.
"""
abstract type MultivariateVolatilitySpec{T, d} <: VolatilitySpec{T} end

"""
	MultivariateARCHModel{T<:AbstractFloat,
	     				  d,
						  VS<:MultivariateVolatilitySpec{T, d},
						  SD<:MultivariateStandardizedDistribution{T, d},
						  MS<:MeanSpec{T}
						 } <: ARCHModel
"""
mutable struct MultivariateARCHModel{T<:AbstractFloat,
									 d,
                 				     VS<:MultivariateVolatilitySpec{T, d},
                 		  	  	     SD<:MultivariateStandardizedDistribution{T, d},
                 				     MS<:MeanSpec{T}
                 				   	} <: ARCHModel
    spec::VS
    data::Matrix{T}
    dist::SD
    meanspec::Vector{MS}
	fitted::Bool
    function MultivariateARCHModel{T, d, VS, SD, MS}(spec, data, dist, meanspec, fitted) where {T, d, VS, SD, MS}
        new(spec, data, dist, meanspec, fitted)
    end
end

dof(am::MultivariateARCHModel{T, d}) where {T, d} = nparams(typeof(am.spec)) + nparams(typeof(am.dist)) + d * nparams(eltype(am.meanspec))

function loglikelihood(am::MultivariateARCHModel)
	sigs = covariances(am)
	z = residuals(am; standardized=true, decorrelated=true)
	n, d = size(am.data)
	return -.5 * (n * d * log(2π) + sum(logdet.(cholesky.(sigs))) + sum(z.^2))
end

"""
	MultivariateARCHModel(spec::MultivariateVolatilitySpec, data::Matrix;
          			  	  dist=MultivariateStdNormal,
					  	  meanspec::[NoIntercept{T}() for _ in 1:d]
		  			  	  fitted::Bool=false
					  	  )
Create a MultivariateARCHModel.
"""
function MultivariateARCHModel(spec::VS,
							   data::Matrix{T};
          					   dist::SD=MultivariateStdNormal{T, d}(),
          				 	   meanspec::Vector{MS}=[NoIntercept{T}() for _ in 1:d], # should come up with a proper multivariate version
		  			 		   fitted::Bool=false
          					  ) where {T<:AbstractFloat,
							   		   d,
                    			 	   VS<:MultivariateVolatilitySpec{T, d},
                   					   SD<:MultivariateStandardizedDistribution,
                   					   MS<:MeanSpec
                   			 		  }
    MultivariateARCHModel{T, d, VS, SD, MS}(spec, data, dist, meanspec, fitted)
end

"""
    predict(am::MultivariateARCHModel, what=:covariance)
Form a 1-step ahead prediction from `am`. `what` controls which object is predicted.
The choices are `:covariance` (the default) or `:correlation`.
"""
function predict(am::MultivariateARCHModel; what=:covariance)
    Ht = covariances(am)
    Rt = correlations(am)
    H = uncond(am.spec)
    R = to_corr(H)
    zt = residuals(am; decorrelated=false)
    at = residuals(am; standardized=false, decorrelated=false)
	T = size(am.data)[1]
	zt = [zt[t, :] for t in 1:T]
	at = [at[t, :] for t in 1:T]
    update!(Ht, Rt, H, R, zt, at, typeof(am.spec), coef(am.spec))
    if what == :covariance
        return Ht[end]
    elseif what == :correlation
        return Rt[end]
    else
        error("Prediction target $what unknown.")
    end
end

# documented in general
fit(am::MultivariateARCHModel; algorithm=BFGS(), autodiff=:forward, kwargs...) = fit(typeof(am.spec), am.data; dist=typeof(am.dist), meanspec=am.meanspec[1], algorithm=algorithm, autodiff=autodiff, kwargs...) # hacky. need multivariate version

# documented in general
function fit!(am::MultivariateARCHModel; algorithm=BFGS(), autodiff=:forward, kwargs...)
    am2 = fit(typeof(am.spec), am.data; meanspec=am.meanspec[1], method=am.spec.method, dist=typeof(am.dist), algorithm=algorithm, autodiff=autodiff, kwargs...)
    am.spec = am2.spec
    am.dist = am2.dist
    am.meanspec = am2.meanspec
    am.fitted = true
    am
end

# documented in general
function simulate(spec::MultivariateVolatilitySpec{T2, d}, nobs;
                  warmup=100,
                  dist::MultivariateStandardizedDistribution{T2}=MultivariateStdNormal{T2, d}(),
                  meanspec::Vector{<:MeanSpec{T2}}=[NoIntercept{T2}() for i = 1:d]
                  ) where {T2<:AbstractFloat, d}
    data = zeros(T2, nobs, d)
	_simulate!(data, spec; warmup=warmup, dist=dist, meanspec=meanspec)
	return MultivariateARCHModel(spec, data; dist=dist, meanspec=meanspec, fitted=false)
end

function _simulate!(data::Matrix{T2}, spec::MultivariateVolatilitySpec{T2, d};
                  warmup=100,
                  dist::MultivariateStandardizedDistribution{T2}=MultivariateStdNormal{T2, d}(),
                  meanspec::Vector{<:MeanSpec{T2}}=[NoIntercept{T2}() for i = 1:d]
                  ) where {T2<:AbstractFloat, d}
	@assert warmup >= 0

	T, d2 = size(data)
	@assert d == d2
	simdata = zeros(T2, T + warmup, d)
	r1 = presample(typeof(spec))
	r2 = maximum(presample.(meanspec))
	r = max(r1, r2)
	r = max(r, 1) # make sure this works for, e.g., ARCH{0}; CircularBuffer requires at least a length of 1
	Ht = CircularBuffer{Matrix{T2}}(r)
	Rt = CircularBuffer{Matrix{T2}}(r)
	zt = CircularBuffer{Vector{T2}}(r)
	at = CircularBuffer{Vector{T2}}(r)

	@inbounds begin
        H = uncond(spec)
		R = to_corr(H)
		all(eigvals(H) .> 0) || error("Model is nonstationary.")
		themean = zeros(T2, d)
        for t = 1: warmup + T
			for i = 1:d
				if t > r2
					ht = getindex.(Ht, i, i)
					lht = log.(ht)
					themean[i] = mean(getindex.(at, i), ht, lht, view(simdata, :, i), meanspec[i], meanspec[i].coefs, t)
				else
					themean[i] = uncond(meanspec[i])
				end
			end

			if t>r1
                update!(Ht, Rt, H, R, zt, at, typeof(spec), coef(spec))
            else
				push!(Ht, H)
				push!(Rt, R)
            end

			z = rand(dist)
			push!(zt, cholesky(Rt[end], check=false).L * z)
			push!(at, sqrt.(diag(Ht[end])) .* zt[end])
			simdata[t, :] .= themean + at[end]
        end
    end
    data .= simdata[warmup + 1 : end, :]
end


text:
function configPotential(y, nodepot, edgepot, es::EdgeStruct)
  nnodes = size(nodepot, 1)
  nedges = size(es.edgeends,1)
  pot = 1
  for n = 1:nnodes
    pot = pot*nodepot[n, y[n]]
  end

  for e = 1:nedges
    i = es.edgeends[e,1]
    j = es.edgeends[e,2]
    pot = pot*edgepot[y[i],y[j],e];
  end
  return pot
end


text:
using Test
using DataValues

@testset "DataValueArray: Constructor" begin

# test Inner Constructor
    @test_throws ArgumentError DataValueArray([1, 2, 3, 4], [true, false, true])

# test (A::AbstractArray, m::Array{Bool}) constructor
    v = [1, 2, 3, 4]
    dv = DataValueArray(v, fill(false, size(v)))

    m = [1 2; 3 4]
    dm = DataValueArray(m, fill(false, size(m)))

    t = Array{Int}(undef, 2, 2, 2)
    t[1:2, 1:2, 1:2] .= 1

    dt = DataValueArray(t, fill(false, size(t)))
    dv = DataValueArray(v, [false, false, false, false])

    y2 = DataValueArray([1, 2, 3, 4, 5, 6],
                    [true, false, false, false, false ,false])
    @test isa(y2, DataValueVector{Int})
    @test y2.isna[1]

# test (::AbstractArray) constructor
    dv = DataValueArray(v)
    @test isa(dv, DataValueVector{Int})

    y = DataValueArray([1, 2, 3, 4, 5, 6])
    @test isa(y, DataValueVector{Int})

    z = DataValueArray(1.:6.)
    @test isa(z, DataValueVector{Float64})

# test (::Type{T}, dims::Dims) constructor
    u1 = DataValueArray{Int}((5,))
    u2 = DataValueArray{Int}((2, 2))
    u3 = DataValueArray{Int}((2, 2, 2))
    @test isa(u1, DataValueVector{Int})
    @test isa(u2, DataValueMatrix{Int})
    @test isa(u3, DataValueArray{Int,3})

# test (::Type{T}, dims::Int...) constructor
    x1 = DataValueArray{Int}(2)
    x2 = DataValueArray{Int}(2, 2)
    x3 = DataValueArray{Int}(2, 2, 2)
    @test isa(x1, DataValueVector{Int})
    @test isa(x2, DataValueMatrix{Int})
    @test isa(x3, DataValueArray{Int,3})

# test DataValueArray{T}(dims::Dims)
    d1, d2 = rand(1:100), rand(1:100)
    X1 = DataValueArray{Int}((d1,))
    X2 = DataValueArray{Int}((d1, d2))
    X3 = DataValueArray{Int}(undef, (d1, d2))
    @test isequal(X1, DataValueArray(Array{Int}(undef, d1), fill(true, (d1,))))
    @test isequal(X2, DataValueArray(Array{Int}(undef, d1, d2), fill(true, (d1, d2))))
    @test isequal(X2, X3)
    for i in 1:5
        m = rand(3:5)
        dims = tuple([ rand(1:5) for i in 1:m ]...)
        X3 = DataValueArray{Int}(dims)
        @test isequal(X3, DataValueArray(Array{Int}(undef, dims), fill(true, dims)))
    end

# test DataValueArray{T,N}(dims::Dims)
    d1, d2 = rand(1:100), rand(1:100)
    X1 = DataValueArray{Int,1}((d1,))
    X2 = DataValueArray{Int,2}((d1, d2))
    X3 = DataValueArray{Int,2}(undef, (d1, d2))
    @test isequal(X1, DataValueArray(Array{Int}(undef, d1), fill(true, (d1,))))
    @test isequal(X2, DataValueArray(Array{Int}(undef, d1, d2), fill(true, (d1, d2))))
    @test isequal(X2, X3)
    for i in 1:5
        m = rand(3:5)
        dims = tuple([ rand(1:5) for i in 1:m ]...)
        X3 = DataValueArray{Int,m}(dims)
        @test isequal(X3, DataValueArray(Array{Int}(undef, dims), fill(true, dims)))
    end

# test DataValueArray{T}(dims::Int...)
    d1, d2 = rand(1:100), rand(1:100)
    X1 = DataValueArray{Int}(d1)
    X2 = DataValueArray{Int}(d1, d2)
    X3 = DataValueArray{Int}(undef, d1, d2)
    @test isequal(X1, DataValueArray(Array{Int}(undef, d1), fill(true, d1)))
    @test isequal(X2, DataValueArray(Array{Int}(undef, d1, d2), fill(true, d1, d2)))
    @test isequal(X2, X3)
    for i in 1:5
        m = rand(3:5)
        dims = [ rand(1:5) for i in 1:m ]
        X3 = DataValueArray{Int}(dims...)
        @test isequal(X3, DataValueArray(Array{Int}(undef, dims...), fill(true, dims...)))
    end

# test DataValueArray{T}(dims::Int...)
    d1, d2 = rand(1:100), rand(1:100)
    X1 = DataValueArray{Int,1}(d1)
    X2 = DataValueArray{Int,2}(d1, d2)
    X3 = DataValueArray{Int,2}(undef, d1, d2)
    @test isequal(X1, DataValueArray(Array{Int,1}(undef, d1), fill(true, d1)))
    @test isequal(X2, DataValueArray(Array{Int,2}(undef, d1, d2), fill(true, d1, d2)))
    @test isequal(X2, X3)
    for i in 1:5
        m = rand(3:5)
        dims = [ rand(1:5) for i in 1:m ]
        X3 = DataValueArray{Int,length(dims)}(dims...)
        @test isequal(X3, DataValueArray(Array{Int}(undef, dims...), fill(true, dims...)))
    end

# test (A::AbstractArray{DataValue}) constructor
    z = DataValueArray([1, NA, 2, NA, 3])
    @test isa(z, DataValueVector{Int})
    @test z.isna[2]
    @test z.isna[4]

# test DataValueArray{T}()
    X = DataValueArray{Int}()
    @test isequal(size(X), ())

# test conversion from arrays, arrays of DataValues and DataValueArrays
    miss1 = [false,false,false,false]
    miss2 = [false,false,false,true]
    for (a, miss) in zip(([1, 2, 3, 4],
                        # 1:4, # Currently does not work on 0.4, cf. JuliaLang/julia#16265
                        DataValue{Int}[DataValue(1), DataValue(2), DataValue(3), DataValue()],
                        DataValueArray(DataValue{Int}[DataValue(1), DataValue(2), DataValue(3), DataValue()])),
                        (miss1, miss2, miss2))
        @test isa(DataValueArray(a), DataValueArray{Int,1})
        @test isequal(DataValueArray(a), DataValueArray{Int,1}([1,2,3,4], miss))
        @test isa(DataValueArray{Int}(a), DataValueArray{Int,1})
        @test isequal(DataValueArray{Int}(a), DataValueArray{Int,1}([1,2,3,4], miss))
        @test isa(DataValueArray{Float64}(a), DataValueArray{Float64,1})
        @test isequal(DataValueArray{Float64}(a), DataValueArray{Float64,1}([1.0,2.0,3.0,4.0], miss))
        @test isa(DataValueArray{Int,1}(a), DataValueArray{Int,1})
        @test isequal(DataValueArray{Int,1}(a), DataValueArray{Int,1}([1,2,3,4], miss))
        @test isa(DataValueArray{Float64,1}(a), DataValueArray{Float64,1})
        @test isequal(DataValueArray{Float64,1}(a), DataValueArray{Float64,1}([1.0,2.0,3.0,4.0], miss))

        @test isa(convert(DataValueArray, a), DataValueArray{Int,1})
        @test isequal(convert(DataValueArray, a), DataValueArray{Int,1}([1,2,3,4], miss))
        @test isa(convert(DataValueArray{Int}, a), DataValueArray{Int,1})
        @test isequal(convert(DataValueArray{Int}, a), DataValueArray{Int,1}([1,2,3,4], miss))
        @test isa(convert(DataValueArray{Float64}, a), DataValueArray{Float64,1})
        @test isequal(convert(DataValueArray{Float64}, a), DataValueArray{Float64,1}([1.0,2.0,3.0,4.0], miss))
        @test isa(convert(DataValueArray{Int,1}, a), DataValueArray{Int,1})
        @test isequal(convert(DataValueArray{Int,1}, a), DataValueArray{Int,1}([1,2,3,4], miss))
        @test isa(convert(DataValueArray{Float64,1}, a), DataValueArray{Float64,1})
        @test isequal(convert(DataValueArray{Float64,1}, a), DataValueArray{Float64,1}([1.0,2.0,3.0,4.0], miss))
    end

# test conversion from Array{DataValue} without element type (issue #145)
    @test isa(DataValueArray(DataValue[DataValue(1), DataValue(2), DataValue(3), DataValue()]), DataValueArray{Any})
    @test isa(DataValueArray{Int}(DataValue[DataValue(1), DataValue(2), DataValue(3), DataValue()]), DataValueArray{Int})
    @test isa(DataValueArray{Float64}(DataValue[DataValue(1), DataValue(2), DataValue(3), DataValue()]), DataValueArray{Float64})

    @test isa(convert(DataValueArray, DataValue[DataValue(1), DataValue(2), DataValue(3), DataValue()]), DataValueArray{Any})
    @test isa(convert(DataValueArray{Int}, DataValue[DataValue(1), DataValue(2), DataValue(3), DataValue()]), DataValueArray{Int})
    @test isa(convert(DataValueArray{Float64}, DataValue[DataValue(1), DataValue(2), DataValue(3), DataValue()]), DataValueArray{Float64})

# test conversion from Any array
    @test convert(DataValueArray{Float64}, Any[1., 2., 3.]) isa DataValueVector{Float64}
    @test convert(DataValueArray{Float64}, Any[1., NA, 3.]) isa DataValueVector{Float64}
    @test convert(DataValueArray{Float64}, Any[DataValue(1.), DataValue{Float64}(), 3.]) isa DataValueVector{Float64}

# converting a DataValueArray to unqualified type DataValueArray should be no-op
    m = rand(10:100)
    A = rand(m)
    M = rand(Bool, m)
    X = DataValueArray(A, M)
    @test X === convert(DataValueArray, X)

end