text:
module CurlConformingFESpacesTests

using Test
using LinearAlgebra
using Gridap.FESpaces
using Gridap.Geometry
using Gridap.TensorValues
using Gridap.CellData
using Gridap.Fields
using Gridap.ReferenceFEs

domain =(0,1,0,1,0,1)
partition = (3,3,3)
model = CartesianDiscreteModel(domain,partition)

order = 2

u(x) = VectorValue(x[1]*x[1],x[1]*x[1]*x[1],0.0)
# u(x) = x

reffe = ReferenceFE(nedelec,order)

V = TestFESpace(model,reffe,dirichlet_tags = [21,22])
# dirichlet_tags = "boundary")
test_single_field_fe_space(V)

U = TrialFESpace(V,u)

uh = interpolate(u,U)

e = u - uh

Ω = Triangulation(model)
dΩ = Measure(Ω,order)

el2 = sqrt(sum( ∫( e⋅e )*dΩ ))
@test el2 < 1.0e-10

# using Gridap.Visualization

# writevtk(trian,"trian",nsubcells=10,cellfields=["uh"=>uh])

end # module


text:
module PipsNLP

using JuMP, Plasmo
using Libdl
using MPI
using Distributed
using DataStructures

export pipsnlp_solve, distribute_optigraph, retrieve_solution

include("pips_nlp_data.jl")
include("distribute.jl")
include("pipsnlp_interface.jl")

using .PipsNlpInterface

end # module

ext: jl


text:
using ArcadeLearningEnvironment
using CartesianGeneticProgramming
using Dates
using IICGP
using Images
using Random
using YAML

function fill_width_with_zeros(x::Matrix{UInt8}, width::Int64)
    if size(x)[2] < width
        x = cat(x, zeros(UInt8, size(x)[1], width-size(x)[2]), dims=2)
    end
    x
end

function width_equalize(v::Matrix{UInt8}, w::Matrix{UInt8})
    width = max(size(v)[2], size(w)[2])
    fill_width_with_zeros(v, width), fill_width_with_zeros(w, width)
end

function cat_along_dim(v::Vector{Matrix{UInt8}}, dim::Int64=2)
    v_cat = v[1]
    if length(v) > 1
        for i in 2:length(v)
            v_cat = cat(v_cat, v[i], dims=dim)
        end
    end
    v_cat
end

function rotate!(v::Vector{Matrix{UInt8}})
    for i in eachindex(v)
        v[i] = transpose(v[i])
    end
end

function buffer_snapshot(enco::CGPInd, active::Vector{Bool})
    s = enco.buffer[1:enco.n_in]
    b = enco.buffer[enco.n_in+1:end][active[enco.n_in+1:end]]
    rotate!(s)
    rotate!(b)
    n_col = max(length(s), length(b))
    s_cat = cat_along_dim(s)
    b_cat = cat_along_dim(b)
    s_cat, b_cat = width_equalize(s_cat, b_cat)
    cat(s_cat, b_cat, dims=1)
end

function save_rgb(game::Game, dir::String, frame::Int64)
    fname = joinpath(dir, string(frame, "_rgb.png"))
    save_screen_png(game, fname)
end

function save_state(s::Vector{Matrix{UInt8}}, dir::String, frame::Int64)
    for i in eachindex(s)
        fname = joinpath(dir, string(frame, "_s", i, ".png"))
        save(fname, transpose(s[i]))
    end
end

function save_enco_buffer(enco::CGPInd, dir::String, frame::Int64)
    for i in eachindex(enco.buffer)
        if i < enco.n_in + 1 || enco.nodes[i].active
            fname = joinpath(dir, string(frame, "_e", i, ".png"))
            save(fname, transpose(enco.buffer[i]))
        end
    end
end

function save_features(f::Vector{Matrix{Float64}}, dir::String, frame::Int64,
                       outputs::Vector{Int16})
    for i in eachindex(f)
        fname = joinpath(dir, string(frame, "_f", outputs[i], ".png"))
        save(fname, transpose(f[i]))
    end
end

function save_cont_buffer(cont::CGPInd, dir::String, frame::Int64)
    data = Dict()
    for i in eachindex(cont.buffer)
        if i < cont.n_in + 1 || cont.nodes[i].active
            data[string(i)] = cont.buffer[i]
        end
    end
    fname = joinpath(dir, string(frame, "_c.yaml"))
    YAML.write_file(fname, data)
end

function save_metadata(metadata::Dict, dir::String, frame::Int64)
    fname = joinpath(dir, string(frame, "_m.yaml"))
    YAML.write_file(fname, metadata)
end

function get_metadata(
    action::T,
    is_sticky::Bool,
    score::Float64,
    e_activated::Vector{Int16},
    e_output::Vector{Int16},
    c_activated::Vector{Int16},
    c_output::Int16
) where {T <: Union{Int32, Int64}}
    metadata = Dict(
        "action"=>action,
        "is_sticky"=>is_sticky,
        "score"=>score,
        "encoder"=>Dict("activated"=>e_activated, "outputs"=>e_output),
        "controller"=>Dict("activated"=>c_activated, "outputs"=>[c_output])
    )
    metadata
end

function recur_activated!(ind, n, activated)
    if !activated[n] # first time seeing this node
        activated[n] = true
        is_input = ind.nodes[n].f == CartesianGeneticProgramming.f_null
        if !is_input
            recur_activated!(ind, ind.nodes[n].x, activated)
            fname = String(Symbol(ind.nodes[n].f))
            arity = IICGP.CGPFunctions.arity[fname]
            if arity > 1
                recur_activated!(ind, ind.nodes[n].y, activated)
            end
        end
    end
end

function find_activated_nodes(ind::CGPInd, out_node::Int16)
    activated = falses(length(ind.nodes))
    recur_activated!(ind, out_node, activated)
    activated_nodes = Vector{Int16}()
    for i in eachindex(activated)
        if activated[i]
            push!(activated_nodes, convert(Int16, i))
        end
    end
    activated_nodes
end

function find_activated_encoder_outputs(
    enco::CGPInd,
    redu::Reducer,
    cont::CGPInd,
    c_activated::Vector{Int16}
)
    @assert redu.parameters["type"] == "pooling"
    fsize = redu.parameters["size"]^2
    e_activated = Vector{Int16}()
    for node in c_activated
        if node < cont.n_in+1 # activated node is a controller input
            activated_enco_output_index = convert(Int64, ceil(node/fsize))
            push!(e_activated, enco.outputs[activated_enco_output_index])
        end
    end
    e_activated
end

function visu_dualcgp_ingame(
    enco::CGPInd,
    redu::Reducer,
    cont::CGPInd,
    game::String,
    seed::Int64,
    max_frames::Int64,
    grayscale::Bool,
    downscale::Bool,
    stickiness::Float64;
    do_save::Bool,
    do_display::Bool,
    buffer_path::String
)
    Random.seed!(seed)
    mt = MersenneTwister(seed)
    g = Game(game, seed)
    img_size = size(get_state(g, grayscale, downscale)[1])
    IICGP.reset!(redu) # zero the buffers
    reward = 0.0
    frames = 1
    prev_action = Int32(0)
    prev_chosen_output = 1
    features = Vector{Matrix{Float64}}()
    active = [enco.nodes[i].active for i in eachindex(enco.nodes)]

    # Init rendering buffer
    if do_display
        visu = []
    end

    while ~game_over(g.ale)
        s = get_state(g, grayscale, downscale)
        rgb = get_rgb(g)
        is_sticky = rand(mt) < stickiness
        if frames == 1
            is_sticky = false
        end
        if !is_sticky
            features, output = IICGP.process_f(enco, redu, cont, s)
            chosen_output = argmax(output)
            action = g.actions[chosen_output]
        else
            chosen_output = prev_chosen_output
            action = prev_action
        end

        # Scan which nodes were activated
        c_output = cont.outputs[chosen_output]
        c_activated = find_activated_nodes(cont, c_output)
        e_output = find_activated_encoder_outputs(enco, redu, cont, c_activated)
        e_activated = Vector{Int16}()
        for node in e_output
            push!(e_activated, find_activated_nodes(enco, node)...)
        end

        # Saving
        if do_save
            metadata = get_metadata(action, is_sticky, reward, e_activated,
                                    e_output, c_activated, c_output)
            save_metadata(metadata, buffer_path, frames)
            save_rgb(g, buffer_path, frames)
            save_enco_buffer(enco, buffer_path, frames)
            save_features(features, buffer_path, frames, enco.outputs)
            save_cont_buffer(cont, buffer_path, frames)
        end

        # Rendering
        if do_display
            snap = buffer_snapshot(enco, active)
            if frames == 1
                visu = snap
            else
                visu = cat(visu, snap, dims=3)
            end
        end
        reward += act(g.ale, action)
        frames += 1
        prev_action = action
        prev_chosen_output = chosen_output
        if frames > max_frames
            break
        end
    end
    close!(g)
    if do_display
        plot_pipeline(visu)
    end

    println()
    println("-"^20)
    println("Game             : ", game)
    println("Seed             : ", seed)
    println("Total return     : ", reward)
    println("Number of frames : ", frames)
    println("-"^20)
    reward
end

function save_graph_struct(
    inds::Array{CGPInd,1},
    saving_dir::String,
    actions::Vector{Int32}
)
    for ind in inds
        data = Dict()
        data["n_in"] = ind.n_in
        data["n_out"] = ind.n_out
        data["nodes"] = []
        data["fs"] = []
        data["edges"] = []
        data["outputs"] = ind.outputs
        for i in eachindex(ind.nodes)
            is_active = ind.nodes[i].active
            if is_active
                x = ind.nodes[i].x
                fname = String(Symbol(ind.nodes[i].f))
                arity = IICGP.CGPFunctions.arity[fname]
                push!(data["nodes"], i)
                push!(data["fs"], fname)
                push!(data["edges"], (x, i))
                if arity == 2
                    y = ind.nodes[i].y
                    push!(data["edges"], (y, i))
                end
            end
        end
        is_controller = typeof(ind.buffer[1]) == Float64
        if is_controller
            data["actions"] = actions
        end
        graph_name = is_controller ? "controller.yaml" : "encoder.yaml"
        graph_path = joinpath(saving_dir, graph_name)
        YAML.write_file(graph_path, data)
    end
end

function get_minimal_actions_set(game::String)
    g = Game(game, 0)
    actions = g.actions
    close!(g)
    actions
end

function visu_ingame(
    exp_dir::String,
    game::String,
    max_frames::Int64;
    do_save::Bool,
    do_display::Bool,
    seed::Int64
)
    cfg = cfg_from_exp_dir(exp_dir)
    # seed = cfg["seed"]
    stickiness = cfg["stickiness"]
    grayscale = cfg["grayscale"]
    downscale = cfg["downscale"]
    is_dualcgp = haskey(cfg, "encoder")
    @assert is_dualcgp
    enco, redu, cont = get_last_dualcgp(exp_dir, game, cfg)
    mini_actions = get_minimal_actions_set(game)
    if do_save
        graph_path = joinpath(exp_dir, "graphs")
        mkpath(graph_path)
        buffer_path = joinpath(exp_dir, "buffers")
        rm(buffer_path, recursive=true, force=true)
        mkpath(buffer_path)
        save_graph_struct([enco, cont], graph_path, mini_actions)
    end
    reward = visu_dualcgp_ingame(
        enco, redu, cont, game, seed, max_frames, grayscale,
        downscale, stickiness, do_save=do_save,
        do_display=do_display, buffer_path=buffer_path
    )
    reward
end

function test_manyvis(
    exp_dir::String,
    game::String,
    max_frames::Int64;
    do_save::Bool,
    do_display::Bool
)
    max_reward = -Inf
    for t = 1:1000
        seed = t
        reward = visu_ingame(exp_dir, game, max_frames, do_save=do_save,
                             do_display=do_display, seed=seed)
        max_reward = max(max_reward, reward)
        if reward > 21
            break
        end
        println("Max total return: ", max_reward)
    end
end

rootdir = joinpath(homedir(), "Documents/git/ICGP-results/")
resdir = joinpath(rootdir, "results/")
min_date = DateTime(2021, 09, 01)
max_date = DateTime(2021, 10, 03)
games = ["boxing", "asteroids", "breakout", "freeway", "gravitar", "riverraid", "space_invaders"]
games = ["freeway"]
reducers = ["pooling"] # Array{String,1}() # ["pooling"]
exp_dirs, games = get_exp_dir(resdir, min_date=min_date, max_date=max_date,
                              games=games, reducers=reducers)
max_frames = 10000
render_graph = false
seed = 0

for i in eachindex(exp_dirs)
    # Generate images (may display / save)
    visu_ingame(exp_dirs[i], games[i], max_frames,
                do_save=true, do_display=false, seed)
    #test_manyvis(exp_dirs[i], games[i], max_frames, do_save=true, do_display=false)

    # Launch python script
    if render_graph
        exp_dir = exp_dirs[i]
        run(`python3.8 pytexgraph.py $exp_dir`)
    end
end


text:
SubWorker{F,T,A,S} = RemoteChannel{Channel{Vector{SubProblem{F,T,A,S}}}}
ScenarioProblemChannel{S} = RemoteChannel{Channel{StochasticPrograms.ScenarioProblems{S}}}
Work = RemoteChannel{Channel{Int}}

function load_subproblems!(subworkers::Vector{SubWorker{F,T,A,S}},
                           scenarioproblems::AbstractScenarioProblems,
                           x::AbstractVector,
                           subsolver::SubSolver) where {F <: AbstractFeasibility,
                                                        T <: AbstractFloat,
                                                        A <: AbstractVector,
                                                        S <: LQSolver}
    # Create subproblems on worker processes
    @sync begin
        for w in workers()
            subworkers[w-1] = RemoteChannel(() -> Channel{Vector{SubProblem{F,T,A,S}}}(1), w)
            @async load_worker!(scenarioproblems, w, subworkers[w-1], x, subsolver)
        end
    end
end

function load_worker!(sp::ScenarioProblems,
                      w::Integer,
                      worker::SubWorker,
                      x::AbstractVector,
                      subsolver::SubSolver)
    n = StochasticPrograms.nscenarios(sp)
    (nscen, extra) = divrem(n, nworkers())
    prev = [nscen + (extra + 2 - p > 0) for p in 2:(w-1)]
    start = isempty(prev) ? 1 : sum(prev) + 1
    stop = min(start + nscen + (extra + 2 - w > 0) - 1, n)
    prev = [begin
            jobsize = nscen + (extra + 2 - p > 0)
            ceil(Int, jobsize)
            end for p in 2:(w-1)]
    start_id = isempty(prev) ? 0 : sum(prev)
    πs = [probability(sp.scenarios[i]) for i = start:stop]
    return remotecall_fetch(init_subworker!,
                            w,
                            worker,
                            sp.parent,
                            sp.problems[start:stop],
                            πs,
                            x,
                            subsolver,
                            start_id)
end

function load_worker!(sp::DScenarioProblems,
                      w::Integer,
                      worker::SubWorker,
                      x::AbstractVector,
                      subsolver::SubSolver)
    prev = [sp.scenario_distribution[p-1] for p in 2:(w-1)]
    start_id = isempty(prev) ? 0 : sum(prev)
    return remotecall_fetch(init_subworker!,
                            w,
                            worker,
                            sp[w-1],
                            x,
                            subsolver,
                            start_id)
end

function init_subworker!(subworker::SubWorker{F,T,A,S},
                         parent::JuMP.Model,
                         submodels::Vector{JuMP.Model},
                         πs::A,
                         x::A,
                         subsolver::SubSolver,
                         start_id::Integer) where {F, T <: AbstractFloat, A <: AbstractArray, S <: LQSolver}
    subproblems = Vector{SubProblem{F,T,A,S}}(undef, length(submodels))
    for (i,submodel) = enumerate(submodels)
        y₀ = convert(A, rand(submodel.numCols))
        subproblems[i] = SubProblem(submodel, parent, start_id + i, πs[i], x, y₀, get_solver(subsolver), F)
    end
    put!(subworker, subproblems)
    return nothing
end

function init_subworker!(subworker::SubWorker{F,T,A,S},
                         scenarioproblems::ScenarioProblemChannel,
                         x::A,
                         subsolver::SubSolver,
                         start_id::Integer) where {F, T <: AbstractFloat, A <: AbstractArray, S <: LQSolver}
    sp = fetch(scenarioproblems)
    subproblems = Vector{SubProblem{F,T,A,S}}(undef, StochasticPrograms.nsubproblems(sp))
    for (i,submodel) = enumerate(sp.problems)
        y₀ = convert(A, rand(sp.problems[i].numCols))
        subproblems[i] = SubProblem(submodel, sp.parent, start_id + i, probability(sp.scenarios[i]), x, y₀, get_solver(subsolver), F)
    end
    put!(subworker, subproblems)
    return nothing
end

function resolve_subproblems!(subworker::SubWorker{F,T,A,S}, x::AbstractVector, cutqueue::CutQueue{T}, aggregator::AbstractAggregator, t::Integer, metadata::MetaData) where {F <: AbstractFeasibility, T <: AbstractFloat, A <: AbstractArray, S <: LQSolver}
    # Fetch all subproblems stored in worker
    subproblems::Vector{SubProblem{F,T,A,S}} = fetch(subworker)
    if isempty(subproblems)
        # Workers has nothing do to, return.
        return nothing
    end
    # Aggregation policy
    aggregation::AbstractAggregation = aggregator(length(subproblems), T)
    # Solve subproblems
    for subproblem ∈ subproblems
        update_subproblem!(subproblem, x)
        cut = subproblem()
        aggregate_cut!(cutqueue, aggregation, metadata, t, cut, x)
    end
    flush!(cutqueue, aggregation, metadata, t, x)
    return nothing
end

function work_on_subproblems!(subworker::SubWorker{F,T,A,S},
                              work::Work,
                              finalize::Work,
                              cutqueue::CutQueue{T},
                              decisions::Decisions{A},
                              metadata::MetaData,
                              aggregator::AbstractAggregator) where {F, T <: AbstractFloat, A <: AbstractArray, S <: LQSolver}
    subproblems::Vector{SubProblem{F,T,A,S}} = fetch(subworker)
    if isempty(subproblems)
       # Workers has nothing do to, return.
       return nothing
    end
    aggregation::AbstractAggregation = aggregator(length(subproblems), T)
    quit = false
    while true
        t::Int = try
            if isready(finalize)
                quit = true
                take!(finalize)
            else
                wait(work)
                take!(work)
            end
        catch err
            if err isa InvalidStateException
                # Master closed the work/finalize channel. Worker finished
                return nothing
            end
        end
        t == -1 && continue
        x::A = fetch(decisions,t)
        for subproblem in subproblems
            update_subproblem!(subproblem, x)
            cut = subproblem()
            !quit && aggregate_cut!(cutqueue, aggregation, metadata, t, cut, x)
        end
        !quit && flush!(cutqueue, aggregation, metadata, t, x)
        if quit
            # Worker finished
            return nothing
        end
    end
end

function eval_second_stage(subworkers::Vector{<:SubWorker}, x::AbstractVector)
    partial_objectives = Vector{Float64}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_objectives[i] = remotecall_fetch(calculate_subobjective, w, subworkers[w-1], x)
        end
    end
    return sum(partial_objectives)
end

function calculate_subobjective(subworker::SubWorker{F,T,A,S}, x::A) where {F, T <: AbstractFloat, A <: AbstractArray, S <: LQSolver}
    subproblems::Vector{SubProblem{F,T,A,S}} = fetch(subworker)
    if length(subproblems) > 0
        return sum([subproblem.π*subproblem(x) for subproblem in subproblems])
    else
        return zero(T)
    end
end

function fill_submodels!(subworkers::Vector{<:SubWorker}, x::AbstractVector, scenarioproblems::ScenarioProblems)
    j = 0
    @sync begin
        for w in workers()
            n = remotecall_fetch((sw)->length(fetch(sw)), w, subworkers[w-1])
            for i = 1:n
                k = i+j
                @async fill_submodel!(scenarioproblems.problems[k],remotecall_fetch((sw,i,x)->begin
                    sp = fetch(sw)[i]
                    sp(x)
                    get_solution(sp)
                end,
                w,
                subworkers[w-1],
                i,
                x)...)
            end
            j += n
        end
    end
    return nothing
end

function fill_submodels!(subworkers::Vector{<:SubWorker}, x::AbstractVector, scenarioproblems::DScenarioProblems)
    @sync begin
        for w in workers()
            @async remotecall_fetch(fill_submodels!,
                                    w,
                                    subworkers[w-1],
                                    x,
                                    scenarioproblems[w-1])
        end
    end
    return nothing
end

function fill_submodels!(subworker::SubWorker{F,T,A,S},
                         x::A,
                         scenarioproblems::ScenarioProblemChannel) where {F <: AbstractFeasibility, T <: AbstractFloat, A <: AbstractArray, S <: LQSolver}
    sp = fetch(scenarioproblems)
    subproblems::Vector{SubProblem{F,T,A,S}} = fetch(subworker)
    for (i, submodel) in enumerate(sp.problems)
        subproblems[i](x)
        fill_submodel!(submodel, subproblems[i])
    end
    return nothing
end


text:
include("Include.jl")

# extra:
using Flux
using Flux: @epochs
using BSON: @save

# load training set -
full_training_data_frame = load_tga_training_data()

# filter -
has_TM_flag = 0
experimental_data_table = filter([:visitid, :TM] => (x, y) -> ((x == 2 || x == 3) && y == has_TM_flag), full_training_data_frame)

# initialize storage for the training data -
training_data = Vector{Tuple{Vector{Float32},Vector{Float32}}}()

# build a model architecture -
deep_coag_model = Chain(Dense(12, 16, σ), Dense(16, 5));

# setup a loss function -
loss(x, y) = Flux.Losses.mae(deep_coag_model(x), y; agg = mean)

# pointer to params -
ps = Flux.params(deep_coag_model)

# # use old school gradient descent -
opt = Momentum(0.1, 0.95)

# main training loop -
(P, D) = size(experimental_data_table)
for i ∈ 1:P

    for j = 1:P

        # leave one out -
        if (i != j)

            # get input and output data -
            input_data = convert.(Float32, Vector(experimental_data_table[j, 3:14]))
            output_data = convert.(Float32, Vector(experimental_data_table[j, 15:19]))
            data_example = (input_data, output_data)

            # capture -
            push!(training_data, data_example)
        end
    end

    # ok, so have the training data for this case -
    # train -
    @epochs 12000 Flux.train!(loss, ps, training_data, opt)

    # save -
    model_name = "deep_coag_model-L$(i)O-TF-TM-$(has_TM_flag)-V2-V3.bson"
    model_file_path = joinpath(_PATH_TO_MODELS, model_name)
    @save model_file_path deep_coag_model

    # need to empty the training data -
    empty!(training_data)
end

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

mutable struct APPStruct

    sourcepath::String;  # path to source code
    currentdir::String;  # current directory
    codename::String; # MDP
    version::String;  # Set version
    appname::String;  # application name
    platform::String; # CPU or GPU
    cpucompiler::String; # Path to CPU compiler
    mpicompiler::String; # Path to MPI compiler
    gpucompiler::String; # Path to GPU compiler
    mpirun::String;      # Path to MPI run command and MPI run options
    cpuflags::String;    # options for CPU compiler
    gpuflags::String;    # options for GGU compiler
    cpumacros::String;    # options for CPU compiler
    gpumacros::String;    # options for GGU compiler
    cuda_arch::String;
    enzyme::String;
    unitstyle::String;
    ensemblemode::String;
    descriptor::String;
    globaloutputs::String;
    peratomoutputs::String;
    potentialfile::String;# APP model file name
    configfile::String;# APP model file name
    weightfile::String;# APP model file name
    model::String;# used to indicate APP model
    buildexec::Int64;
    modelfile::String;# APP model file name
    modelnumber::Int64;
    preprocessmode::Int64; # preprocessing mode
    configmode::Int64;
    weightmode::Int64;
    mpiprocs::Int64; # number of MPI ranks
    globalfreq::Int64;
    peratomfreq::Int64;
    unitstylenum::Int64;
    ensemblemodenum::Int64;
    outputflag::Int64

    datapath::String;
    datafile::String;
    dataformat::String;
    datamode::Int64;

    traininglist::Array{Int64,2}; # a list of configurations for training the potential
    validatelist::Array{Int64,2}; # a list of configurations for validating the potential

    dim::Int64; # physical dimension
    nconfigs::Int64; # number of configurations
    maxnumneighbors::Int64; # maximum number of neighbors allowed
    ntimesteps::Int64;# number of time steps
    nab::Int64;# number of atoms per block
    backend::Int64;# backend computing platform
    training::Int64;  # 0 -> no training, 1 -> Linear regression, 2 -> Gaussian process, 3 -> Neural net
    dftdata::Int64;   # 0 -> no data, 1 -> energies only, 2 -> forces only, 3 -> energies and forces
    runMD::Int64;        # 0 no MD simulation, 1 -> run MD simulation
    potentialform::Int64;# 0 -> empirical potential, 1 -> ML potential, 2 -> combined potential
    #neighpair::Int64;    # 0 -> full neighbor list, 1 -> half neighbor list for pair potentials
    #neighcell::Int64;    # 0 -> O(N^2) algorithm, 1 -> Cell-linked list algorithm to form neighbor list
    decomposition::Int64;# 0 -> force decomposition, 1 -> atom decomposition
    energycal::Int64;    # turns energy calculation on or off
    forcecal::Int64;     # turns force calculation on or off
    stresscal::Int64;    # turns stress calculation on or off

    neighskin::Float64;
    neighcell::Int64;       # 0 -> O(N^2) algorithm, 1 -> Cell-linked list algorithm to form neighbor list
    neighmax::Int64;      # maximum number of neighbors allowed
    neighevery::Int64;      # perform neighbor rebuild list for "every" iterations
    neighdelay::Int64;      # delay neighbor rebuild
    neighcheck::Int64;      # 0 -> neighbor rebuild by delay and every settings,
                             # 1 -> neighbor rebuild occurs if some atom has moved more than half the skin distance
    neighpair::Int64;       # 0 -> full neighbor list, 1 -> half neighbor list for pair potentials

    time::Float64;       # initial time
    dt::Float64;         # time step

    coeff::Array{Float64,2};   # coefficients for the potential
    rcutsqmax::Float64;  # square of maximum cutoff radius
    boxoffset::Array{Float64,2}; # offset for simulation box to account for periodic boundary conditions
    bcs::Array{Int64,2}; # boundary conditions
    pbc::Array{Int64,2};  # periodic boundary conditions

    we::Float64;  # energy weight
    wf::Float64;  # force weight

    # atom types
    natomtype::Int64;  # number of atom types
    atomspecies::Array{String,1};
    atommasses::Array{Float64,2};
    atomcharges::Array{Float64,2};
    atomnumbers::Array{Int64,2};

    # molecue types
    nmoletype::Int64;  # number of molecule types
    moleinfo::Array{Int64,2}; # information about each molecule type

    # simulation and solver parameters
    ndims::Array{Int64,2};     # a list of integers indicating dimensions
    flag::Array{Int64,2};      # a list of flags
    simparam::Array{Float64,2}; # simulation parameters
    solparam::Array{Float64,2}; # solver parameters
    nveparam::Array{Float64,2}; # NVE parameters
    nvtparam::Array{Float64,2}; # NVT parameters
    createvelocity::Array{Float64,2}; # NVT parameters
    setatomtypefraction::Array{Float64,2};

    snaprcutfac::Float64;
    snaprfac0::Float64;
    snaprmin0::Float64;
    snaptwojmax::Int64;
    snapbzeroflag::Int64;
    snapquadraticflag::Int64;
    snapswitchflag::Int64;
    snapchemflag::Int64;
    snapbnormflag::Int64;
    snapwselfallflag::Int64;
    snapnelem::Int64;
    snapncoeff::Int64;
    snaparam::Array{Float64,2};
    snapelemradius::Array{Float64,2};
    snapelemweight::Array{Float64,2};
    snapcoeff::Array{Float64,2};
    snapcoefffile::String;

    # machine learning potentials
    muml::Array{Float64,2};   # coefficients of ML potential
    K::Int64;          # degree of radial basis functions
    L::Int64;          # degree of spherical harmonics
    rcutml::Float64;    # cut-off radius for machine learning potential
    descriptornum::Int64; # descriptor flag: 0 -> Spherical Harmonics Bessel
    spectrum::Int64;   # spectrum flag: 0-> power spectrum, 1-> bispectrum, 2-> both power and bispectrum
    chemtype::Int64;   # 0 -> single atom-type basis functions, 1 -> double atom-type basis functions
    ncmuml::Int64;     # length of mu0

    # empirical potential parameters
    eta::Array{Float64,2};   # hyperparameters for emprical potentials
    kappa::Array{Int64,2}; # flag parameters for emprical potentials

    # nonbonded single potentials
    np1a::Int64;           # number of nonbonded single potentials
    pot1a::Array{Int64,2};     # list of nonbonded single potentials
    mu1a::Array{Float64,2};  # parameters for nonbonded single potentials
    ncmu1a::Int64;         # length of mu1a

    # bonded single potentials
    np1b::Int64;           # number of bonded single potentials
    pot1b::Array{Int64,2};      # list of bonded single potentials
    mu1b::Array{Float64,2};   # parameters for bonded single potentials
    ncmu1b::Int64;         # length of mu1a
    atom1b::Array{Int64,2};     # list of bonded atom types

    # nonbonded pair potentials
    np2a::Int64;           # number of nonbonded pair potentials
    pot2a::Array{Int64,2};     # list of nonbonded pair potentials
    mu2a::Array{Float64,2};  # parameters for all nonbonded pair potentials
    ncmu2a::Int64;         # length of mu2a
    rcut2a::Array{Float64,2}; # cut-off radius for each nonbonded pair potential

    # bonded pair potentials
    np2b::Int64;               # number of bonded pair potentials
    pot2b::Array{Int64,2};       # list of bonded pair potentials
    mu2b::Array{Float64,2};  # parameters for all bonded pair potentials
    ncmu2b::Int64;             # length of mu2b
    rcut2b::Array{Float64,2};   # cut-off radius for each nonbonded pair potential
    atom2b::Array{Int64,2}; # list of bonded atom pair types

    # two-body bond order potentials
    np2c::Int64;               # number of two-body bond order potentials
    pot2c::Array{Int64,2};         # list of two-body bond order potentials
    mu2c::Array{Float64,2};  # parameters for all two-body bond order potentials
    ncmu2c::Int64;             # length of mu2c
    rcut2c::Array{Float64,2};    # cut-off radius for each two-body bond order potential
    atom2c::Array{Int64,2};        # list of bond porder atom types

    # nonbonded triplet potentials
    np3a::Int64;           # number of nonbonded triplet potentials
    pot3a::Array{Int64,2};     # list of nonbonded triplet potentials
    mu3a::Array{Float64,2};  # parameters for all nonbonded triplet potentials
    ncmu3a::Int64;         # length of mu3a
    rcut3a::Array{Float64,2}; # cut-off radius for each nonbonded triplet potential

    # bonded triplet potentials
    np3b::Int64;               # number of bonded triplet potentials
    pot3b::Array{Int64,2};         # list of bonded triplet potentials
    mu3b::Array{Float64,2};  # parameters for all bonded triplet potentials
    ncmu3b::Int64;             # length of mu3b
    rcut3b::Array{Float64,2};    # cut-off radius for each nonbonded triplet potential
    atom3b::Array{Int64,2}; # list of bonded atom triplet types

    # three-body bond order potentials
    np3c::Int64;               # number of three-body bond order potentials
    pot3c::Array{Int64,2};         # list of three-body bond order potentials
    mu3c::Array{Float64,2};  # parameters for all three-body bond order potentials
    ncmu3c::Int64;             # length of mu3c
    rcut3c::Array{Float64,2};    # cut-off radius for each three-body bond order potential
    atom3c::Array{Int64,2};   # list of three-body bond order atom types

    # nonbonded quadruplet potentials
    np4a::Int64;           # number of nonbonded quadruplet potentials
    pot4a::Array{Int64,2};     # list of nonbonded quadruplet potentials
    mu4a::Array{Float64,2};  # parameters for all nonbonded quadruplet potentials
    ncmu4a::Int64;          # length of mu4a
    rcut4a::Array{Float64,2}; # cut-off radius for each nonbonded quadruplet potential

    # bonded quadruplet potentials
    np4b::Int64;               # number of bonded quadruplet potentials
    pot4b::Array{Int64,2};         # list of bonded quadruplet potentials
    mu4b::Array{Float64,2};  # parameters for all bonded quadruplet potentials
    ncmu4b::Int64;             # length of mu3b
    rcut4b::Array{Float64,2};    # cut-off radius for each nonbonded quadruplet potential
    atom4b::Array{Int64,2}; # list of bonded atom quadruplet types

    ncx::Int64; # number of compoments of x
    ncv::Int64; # number of compoments of v
    nce::Int64; # number of compoments of e
    ncf::Int64; # number of compoments of f
    ncq::Int64; # number of compoments of q
    ncmu::Int64; # number of compoments of mu
    nceta::Int64; # number of compoments of eta
    nckappa::Int64; # number of compoments of kappa

    lattice
    domain
    region
    setforce
    lineforce
    planeforce
    wallharmonic
    walllj126
    wallmorse
    wallreflect
    setvelocity
    velocityrescaling

    APPStruct() = new();
end

function initializeapp(cdir,sourcepath,version="version0")
    app = APPStruct();

    app.sourcepath = sourcepath;
    app.currentdir = cdir;
    app.codename = "MDP";
    app.version = version;
    app.appname = "app";
    app.platform = "cpu";
    app.cpucompiler = "clang++";
    app.mpicompiler = "mpicxx";
    app.gpucompiler = "clang++";
    app.cuda_arch = "sm_60";
    app.enzyme = "";
    app.mpirun = "mpirun";
    app.cpuflags = "-O2 -ldl -lm -lblas -llapack";
    app.gpuflags = "-lcudart -lcublas";
    app.cpumacros = "";
    app.gpumacros = "";
    app.unitstyle = "";
    app.descriptor = "";
    app.ensemblemode = "";
    app.potentialfile = "";
    app.configfile = "";
    app.weightfile = "";
    app.buildexec = 0;
    app.configmode = 0;
    app.weightmode = 0;
    app.preprocessmode = 1;
    app.mpiprocs = 1;

    app.datapath = "";
    app.datafile = "";
    app.dataformat = "";
    app.datamode = 1;
    app.outputflag = 0

    app.traininglist = reshape([], 0, 2);
    app.validatelist = reshape([], 0, 2);

    app.dim = 3;         # physical dimension
    app.nconfigs = 1;    # number of configurations
    app.maxnumneighbors = 10; # maximum number of neighbors allowed
    app.ntimesteps = 0;  # number of time steps
    app.nab = 4096;      # number of atoms per block
    app.backend = 0;     # backend computing platform
    app.training = 0;     # 0 -> no training, 1 -> Linear regression, 2 -> Gaussian process, 3 -> Neural net
    app.dftdata = 0;      # 0 -> no data, 1 -> energies only, 2 -> forces only, 3 -> energies and forces
    app.runMD = 0;        # 0 no MD simulation, 1 -> run MD simulation
    app.potentialform = 0;# 0 -> empirical potential, 1 -> ML potential, 2 -> combined potential
    app.decomposition = 0;# 0 -> force decomposition, 1 -> atom decomposition
    app.energycal = 0;    # turns energy calculation on or off
    app.forcecal = 0;     # turns force calculation on or off
    app.stresscal = 0;    # turns stress calculation on or off
    app.time = 0.0;       # initial time
    app.dt = 0.0;         # time step
    app.rcutsqmax = 0.0;  # square of maximum cutoff radius
    app.boxoffset = [0.0 0.0 0.0]; # offset for simulation box to account for periodic boundary conditions

    app.globalfreq = 1;       # frequency to print and save default global quantities (pe, ke, ce, te, temp, press)
    app.peratomfreq = 1;     # frequency to save default peratom quantities (id, t, x, f)
    app.globaloutputs = "";    # additional global outputs (stresses, com, vcm)
    app.peratomoutputs = "";   # additional peratom outputs (velocity, mass, virial, image)

    app.neighskin = 0.0;
    app.neighcell = 0;       # 0 -> O(N^2) algorithm, 1 -> Cell-linked list algorithm to form neighbor list
    app.neighmax = 100;      # maximum number of neighbors allowed
    app.neighevery = 1;      # perform neighbor rebuild list for "every" iterations
    app.neighdelay = 0;      # delay neighbor rebuild
    app.neighcheck = 1;      # 0 -> neighbor rebuild by delay and every settings,
                             # 1 -> neighbor rebuild occurs if some atom has moved more than half the skin distance
    app.neighpair = 0;       # 0 -> full neighbor list, 1 -> half neighbor list for pair potentials

    app.coeff = reshape([], 0, 2);
    app.bcs = [0 0 0 0 0 0]; # boundary conditions
    app.pbc = [1 1 1];       # periodic boundary conditions

    app.we = 1.0; # energy weight
    app.wf = 1.0; # force weight

    # atom types
    app.natomtype = 0;  # number of atom types
    app.atommasses = reshape([], 0, 2);
    app.atomcharges = reshape([], 0, 2);
    app.atomnumbers = reshape([], 0, 2);

    # molecue types
    app.nmoletype = 0;  # number of molecule types
    app.moleinfo = reshape([], 0, 2); # information about each molecule type

    # simulation and solver parameters
    app.ndims = reshape([], 0, 2);     # a list of integers indicating dimensions
    app.flag = reshape([], 0, 2);      # a list of flags
    app.simparam = reshape([], 0, 2); # simulation parameters
    app.solparam = reshape([], 0, 2); # solver parameters
    app.nveparam = reshape([], 0, 2); # NVE parameters
    app.nvtparam = reshape([], 0, 2); # NVT parameters
    app.createvelocity = reshape([], 0, 2); # NVT parameters
    app.setatomtypefraction = reshape([], 0, 2);

    # machine learning potentials
    app.muml = reshape([], 0, 2);   # coefficients of ML potential
    app.K = 0;          # degree of radial basis functions
    app.L = 0;          # degree of spherical harmonics
    app.rcutml = 0.0;    # cut-off radius for machine learning potential
    app.spectrum = 2;   # spectrum flag: 0-> power spectrum, 1-> bispectrum, 2-> both power and bispectrum
    app.chemtype = 0;   # 0 -> single atom-type basis functions, 1 -> double atom-type basis functions
    app.ncmuml = 0;     # length of mu0

    app.snaprcutfac = 0;
    app.snaptwojmax = 0;
    app.snaprfac0 = 0;
    app.snaprmin0 = 0;
    app.snapbzeroflag = 1;
    app.snapquadraticflag = 0;
    app.snapswitchflag = 1;
    app.snapchemflag = 0;
    app.snapbnormflag = 0;
    app.snapwselfallflag = 0;
    app.snapnelem = 0;    # # of elements
    app.snapncoeff = 0;  # # of coeffients per element
    app.snaparam = reshape([], 0, 2); # snap parameters
    app.snapelemradius = reshape([], 0, 2); # radius per element
    app.snapelemweight = reshape([], 0, 2);   # weight per element
    app.snapcoeff = reshape([], 0, 2);
    app.snapcoefffile = "";

    # empirical potential parameters
    app.eta = reshape([], 0, 2);   # hyperparameters for emprical potentials
    app.kappa = reshape([], 0, 2); # flag parameters for emprical potentials

    # nonbonded single potentials
    app.np1a = 0;           # number of nonbonded single potentials
    app.pot1a = reshape([], 0, 2);     # list of nonbonded single potentials
    app.mu1a = reshape([], 0, 2);  # parameters for nonbonded single potentials
    app.ncmu1a = 0;         # length of mu1a

    # bonded single potentials
    app.np1b = 0;           # number of bonded single potentials
    app.pot1b = reshape([], 0, 2);      # list of bonded single potentials
    app.mu1b = reshape([], 0, 2);   # parameters for bonded single potentials
    app.ncmu1b = 0;         # length of mu1a
    app.atom1b = reshape([], 0, 2);     # list of bonded atom types

    # nonbonded pair potentials
    app.np2a = 0;           # number of nonbonded pair potentials
    app.pot2a = reshape([], 0, 2);     # list of nonbonded pair potentials
    app.mu2a = reshape([], 0, 2);  # parameters for all nonbonded pair potentials
    app.ncmu2a = 0;         # length of mu2a
    app.rcut2a = reshape([],0,2); # cut-off radius for each nonbonded pair potential

    # bonded pair potentials
    app.np2b = 0;               # number of bonded pair potentials
    app.pot2b = reshape([], 0, 2);       # list of bonded pair potentials
    app.mu2b = reshape([], 0, 2);  # parameters for all bonded pair potentials
    app.ncmu2b = 0;             # length of mu2b
    app.rcut2b = reshape([],0,2);   # cut-off radius for each nonbonded pair potential
    app.atom2b = reshape([], 0, 2); # list of bonded atom pair types

    # two-body bond order potentials
    app.np2c = 0;               # number of two-body bond order potentials
    app.pot2c = reshape([], 0, 2);         # list of two-body bond order potentials
    app.mu2c = reshape([], 0, 2);  # parameters for all two-body bond order potentials
    app.ncmu2c = 0;             # length of mu2c
    app.rcut2c = reshape([],0,2);    # cut-off radius for each two-body bond order potential
    app.atom2c = reshape([], 0, 2);        # list of bond porder atom types

    # nonbonded triplet potentials
    app.np3a = 0;           # number of nonbonded triplet potentials
    app.pot3a = reshape([], 0, 2);     # list of nonbonded triplet potentials
    app.mu3a = reshape([], 0, 2);  # parameters for all nonbonded triplet potentials
    app.ncmu3a = 0;         # length of mu3a
    app.rcut3a = reshape([],0,2); # cut-off radius for each nonbonded triplet potential

    # bonded triplet potentials
    app.np3b = 0;               # number of bonded triplet potentials
    app.pot3b = reshape([], 0, 2);         # list of bonded triplet potentials
    app.mu3b = reshape([], 0, 2);  # parameters for all bonded triplet potentials
    app.ncmu3b = 0;             # length of mu3b
    app.rcut3b = reshape([],0,2);    # cut-off radius for each nonbonded triplet potential
    app.atom3b = reshape([], 0, 2); # list of bonded atom triplet types

    # three-body bond order potentials
    app.np3c = 0;               # number of three-body bond order potentials
    app.pot3c = reshape([], 0, 2);         # list of three-body bond order potentials
    app.mu3c = reshape([], 0, 2);  # parameters for all three-body bond order potentials
    app.ncmu3c = 0;             # length of mu3c
    app.rcut3c = reshape([],0,2);    # cut-off radius for each three-body bond order potential
    app.atom3c = reshape([], 0, 2);   # list of three-body bond order atom types

    # nonbonded quadruplet potentials
    app.np4a = 0;           # number of nonbonded quadruplet potentials
    app.pot4a = reshape([], 0, 2);     # list of nonbonded quadruplet potentials
    app.mu4a = reshape([], 0, 2);  # parameters for all nonbonded quadruplet potentials
    app.ncmu4a = 0;          # length of mu4a
    app.rcut4a = reshape([],0,2); # cut-off radius for each nonbonded quadruplet potential

    # bonded quadruplet potentials
    app.np4b = 0;               # number of bonded quadruplet potentials
    app.pot4b = reshape([], 0, 2);         # list of bonded quadruplet potentials
    app.mu4b = reshape([], 0, 2);  # parameters for all bonded quadruplet potentials
    app.ncmu4b = 0;             # length of mu3b
    app.rcut4b = reshape([],0,2);    # cut-off radius for each nonbonded quadruplet potential
    app.atom4b = reshape([], 0, 2); # list of bonded atom quadruplet types

    app.ncx = 0; # number of compoments of x
    app.ncv = 0; # number of compoments of v
    app.nce = 0; # number of compoments of e
    app.ncf = 0; # number of compoments of f
    app.ncq = 0; # number of compoments of q
    app.ncmu = 0; # number of compoments of mu
    app.nceta = 0; # number of compoments of eta
    app.nckappa = 0; # number of compoments of kappa

    app.lattice = [];
    app.region = [];
    app.domain = [];
    app.setforce = nothing
    app.lineforce = nothing
    app.planeforce = nothing
    app.wallharmonic = nothing
    app.walllj126 = nothing
    app.wallmorse = nothing
    app.wallreflect = nothing
    app.setvelocity = nothing
    app.velocityrescaling = nothing

    return app;

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
#=
    Other statistics quoted in text.
=#

##############################################################
## Load model and results
##############################################################

include("../Results/MainResult_Setup.jl")
@load "Results/MainResult.jld2"

##############################################################
## Lower bound on σ_λ
##############################################################

σλ = Matrix(C)[8,:]
quantile(σλ,[0.025,0.975])

##############################################################
## ρ_Rβ < 0
##############################################################

ρRβ = Matrix(C)[14,:]
count(ρRβ .< 0.0) / length(ρRβ)

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
export NeuralNetworkQ, update!

using Flux
using CUDAapi

import Flux:params
import Zygote:gradient
import .Utils:to_device, to_host

# TODO: Add static size check for legal input data type.

"""
    NeuralNetworkQ(;kwargs...) -> NeuralNetworkQ{D, M, O, P}

Use neural networks to generate estimations of state-action values.

# Keywords

- `model::M`: describes the network structure.
- `optimizer::O`: defines how to update parameters given grads.
- `parameters::P=params(model)`: the parameters of `model`.
- `device::Symbol=:cpu`: the param `D` of `NeuralNetworkQ`, specify where to run the model. Supported keywords are:
    - `:cpu`
    - `:gpu`, if the specified device is `:gpu`, then we will automatically change the device to one of the following device types according to the `backend(model)`:
        - `:Knet_gpu`, means the model is `Knet` that runs on gpu.
        - `:Zygote_gpu`, means the model is `Flux` with `Zygote` that runs on gpu.

# Fields

- `model::M`
- `optimizer::O`
- `params::P`

"""
struct NeuralNetworkQ{D, M, O, P} <: AbstractQApproximator
    model::M
    optimizer::O
    params::P
end

function NeuralNetworkQ(;model::M, optimizer::O, parameters::P=params(model), device::Symbol=:cpu) where {M, O, P}
    if device == :gpu
        has_cuda() || throw(ArgumentError("the specified device is gpu but can not find cuda!"))
        device = Symbol(backend(model), "_gpu")
    end
    m = to_device(Val(device), model)
    ps = params(m)
    NeuralNetworkQ{device, typeof(m), O, typeof(ps)}(m, optimizer, ps)
end

params(Q::NeuralNetworkQ) = Q.params

to_device(Q::NeuralNetworkQ{D}, x) where D = to_device(Val(D), x)

"get Q value of the batch"
batch_estimate(Q::NeuralNetworkQ, states) = Q.model(states)

"get Q value of all actions"
(Q::NeuralNetworkQ)(s) = Q.model(to_device(Q, reshape(s, size(s)..., 1))) |> to_host |> drop_last_dim

"get Q value of some specific action"
(Q::NeuralNetworkQ)(s, a) = Q(s)[a]

function update!(Q::NeuralNetworkQ, gs)
    Flux.Optimise.update!(Q.optimizer, Q.params, gs)
end

Base.copyto!(dest::NeuralNetworkQ, src::NeuralNetworkQ) = Flux.loadparams!(dest.model, src.params)

gradient(f, Q::NeuralNetworkQ) = gradient(f, Val(backend(Q.model)), params(Q))


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

using JuliaFEM, Test

abstract type PlaneStressElasticityProblem <: AbstractProblem end

function PlaneStressElasticityProblem(dim::Int=2, elements=[])
    return Problem{PlaneStressElasticityProblem}(dim, elements)
end

function get_unknown_field_name(::Type{P}) where P<:PlaneStressElasticityProblem
    return "displacement"
end

function get_unknown_field_type(::Type{P}) where P<:PlaneStressElasticityProblem
    return Vector{Float64}
end

function get_residual_vector(problem::Problem{PlaneStressElasticityProblem}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)

    basis = element(ip, time)
    dbasis = element(ip, time, Val{:grad})

    # material parameters
    E = element("youngs modulus", ip, time)
    nu = element("poissons ratio", ip, time)
    mu = E/(2*(1+nu))
    la = E*nu/((1+nu)*(1-2*nu))
    la = 2*la*mu/(la + 2*mu)  # <- correction for 2d

    # elasticity formulation
    u = element("displacement", ip, time, variation)
    gradu = element("displacement", ip, time, Val{:grad}, variation)
    F = I + gradu

    E = 1/2*(F'*F - I)
    S = la*trace(E)*I + 2*mu*E
    r = F*S*dbasis

    b = element("displacement volume load", ip, time)
    r -= b*basis

    return vec(r)
end

function test_residual_form()
    # create model -- start
    element = Quad4([1, 2, 3, 4])
    element["geometry"] = Vector[[0.0,0.0], [10.0,0.0], [10.0,1.0], [0.0,1.0]]
    element["youngs modulus"] = 500.0
    element["poissons ratio"] = 0.3
    element["displacement volume load"] = Vector[[0.0,-10.0], [0.0,-10.0], [0.0,-10.0], [0.0,-10.0]]
    element["displacement"] = (0.0 => Vector{Float64}[zeros(2) for i=1:length(element)])
    problem = PlaneStressElasticityProblem()
    push!(problem, element)
    # create model -- end

    free_dofs = [3, 4, 5, 6]
    solve!(problem, free_dofs, 0.0)  # launch a newton solver for single element
    disp = element("displacement", [1.0, 1.0], 0.0)
    @info("displacement at tip: $disp")

    # verified using Code Aster.
    @test isapprox(disp[2], -8.77303119819776E+00)
end


text:
immutable NGramModel
    ml::Dict{UTF8String, Float64}
    n::Integer
end

getindex(m::NGramModel, gram::String) = get(m.ml, utf8(gram), 0.0)

function NGramModel(sentences::Vector{String}, n)

    # Tokenize string
    tokenize(s) = TextAnalysis.tokenize(TextAnalysis.EnglishLanguage, s)

    # NGramize tokens
    function ngramize(words)
        # Use "*" symbol start and "STOP" for end of sentence
        words = vcat("*", "*", words, "STOP")
        n_words = length(words)
        tokens = Dict{UTF8String, Int}()

        for m in 1:n
            for index in 1:(n_words - m + 1)
                token = join(words[index:(index + m - 1)], " ")
                tokens[token] = get(tokens, token, 0) + 1
            end
        end

        tokens
    end

    # Create a NGramDocument from a gram dict
    document(gram) = NGramDocument(gram, n)
    togeneric(docs) = convert(Array{GenericDocument, 1}, docs)

    # Create corpus
    corpus= Corpus(map(x -> x |> tokenize |> ngramize |> document, sentences) |> togeneric)

    # Create lexicon
    update_lexicon!(corpus)

    # Create model with maximum likelihood estimation
    # => Linear approximation
    ml = Dict{UTF8String, Float64}()
    words = collect(UTF8String, keys(lexicon(corpus)))
    n_words = length(words)
    gramcount = lexicon(corpus)

    function likelihood(gram)
        println(join(gram, " "))
        res = 0

        for i = 1:(n-1)
            println("$(i) => $(join(gram[i:end], " ")) / $(join(gram[i:(end - 1)], " "))")
            res += (1 / n) * gramcount[join(gram[i:end], " ")] / gramcount[join(gram[i:(end - 1)], " ")]
        end
        println("$(n) => $(gram[end]) / $(n_words)")
        res += (1 / n) * gramcount[gram[end]] / n_words

        res
    end

    for w in words
        trigram = tokenize(w)

        if length(trigram) == n
            ml[w] = likelihood(trigram)
        end
    end

    NGramModel(ml, n)
end


text:
using Optim
using FANN
using Base.Test

net = ANN([10, 1], [:linear], b=0.1, errorfunc=:linear)

X = randn(10, 100)
y = sum(X, 1)
dset = DataSet(X, y)

# initial condition
w0 = weights(net)

# objective fucntion
f(w) = (weights!(net, w); mse(net, dset))
optimize(f, w0, ftol=1e-6) |> println

weights!(net, w0)
train!(net, dset, max_epochs=1000, epochs_between_reports=0, desired_error=1e-6)


text:
using Distributions
import Distributions.logpdf

# TODO: iid is currently a bit weird. We'd like it to allow both a specific multiplicity and "on demand" (unspecified but flexible). It's clear how to do either, but not yet if this will allow a single representation.

export iid
struct iid{N,D}
    size :: N
    dist :: D
end

iid(t::Int...) = dist -> iid(t,dist)

# iid(dist) = iid(Nothing, dist)

function Base.iterate(d::iid)
    n = prod(d.size) * length(d.dist)
    return (d.dist, (n, 1))
end

function Base.iterate(d::iid, state)
    (n, count) = state
    count < n && return (d.dist, (n, count + 1))
    return nothing
end

# Distributions.cdf(iid, x) = cdf.(iid.dist, x)
# Distributions.quantile(iid, x) = quantile.(iid.dist, x)

import Base.length
Base.length(d::iid) = prod(d.size)

import Base.eltype
Base.eltype(d::iid) = typeof(d.dist)

function Base.rand(d::iid)
    T = typeof(rand(d.dist))
    x = Array{T}(undef, d.size)
    for cartix in CartesianIndices(x)
        ix = Tuple(cartix)
        @inbounds setindex!(x,rand(d.dist), ix...)
    end
    x
end

function Distributions.logpdf(d::iid,x)
    s = zero(Float64)
    Δs(xj) = logpdf(d.dist, xj)

    @inbounds @simd for j in eachindex(x)
        s += Δs(x[j])
    end
    s
end

# function Base.rand(d::Union{iid, HalfNormal}, n::Int)
#     x1 = rand(d)
#     x = Array{typeof(x1),1}(undef, n)
#     for j in eachindex(x)
#         @inbounds x[j] = rand(d)
#     end
#     x
# end


text:
Graphs.flip_direction(edge::Edge{Int32}) = Edge(-edge.data)

@testset "graphs" begin
    @testset "disconnected" begin
        Random.seed!(8)
        graph = DirectedGraph{Vertex{Int64}, Edge{Float64}}()
        verts = [Vertex(rand(Int64)) for i = 1 : 10]
        for v in verts
            add_vertex!(graph, v)
        end
        @test num_vertices(graph) == length(verts)
        @test num_edges(graph) == 0
        for v in verts
            @test length(in_edges(v, graph)) == 0
            @test length(out_edges(v, graph)) == 0
            @test length(in_neighbors(v, graph)) == 0
            @test length(out_neighbors(v, graph)) == 0
        end
        @test isempty(setdiff(vertices(graph), verts))
        @test isempty(edges(graph))
        show(devnull, graph)
    end

    @testset "tree graph" begin
        Random.seed!(9)
        graph = DirectedGraph{Vertex{Int64}, Edge{Float64}}()
        root = Vertex(rand(Int64))
        add_vertex!(graph, root)
        nedges = 15
        for i = 1 : nedges
            parent = rand(vertices(graph))
            child = Vertex(rand(Int64))
            edge = Edge(rand())
            add_edge!(graph, parent, child, edge)
        end
        @test num_vertices(graph) == nedges + 1
        @test num_edges(graph) == nedges
        for v in vertices(graph)
            if v == root
                @test length(in_edges(v, graph)) == 0
                @test length(out_edges(v, graph)) > 0
                @test length(in_neighbors(v, graph)) == 0
                @test length(out_neighbors(v, graph)) > 0
            else
                @test length(in_edges(v, graph)) == 1
                @test length(in_neighbors(v, graph)) == 1
            end
        end
        show(devnull, graph)
    end

    @testset "remove_vertex!" begin
        Random.seed!(10)
        graph = DirectedGraph{Vertex{Int64}, Edge{Float64}}()
        edge = Edge(rand())
        add_edge!(graph, Vertex(rand(Int64)), Vertex(rand(Int64)), edge)
        for v in vertices(graph)
            @test_throws ErrorException remove_vertex!(graph, v)
        end

        graph = DirectedGraph{Vertex{Int64}, Edge{Float64}}()
        for i = 1 : 100
            add_vertex!(graph, Vertex(i))
        end
        for i = 1 : num_vertices(graph) - 1
            add_edge!(graph, rand(vertices(graph)), rand(vertices(graph)), Edge(Float64(i)))
        end
        original = deepcopy(graph)

        vertex = rand(collect(filter(v -> isempty(in_edges(v, graph)) && isempty(out_edges(v, graph)), vertices(graph))))
        remove_vertex!(graph, vertex)
        @test vertex ∉ vertices(graph)
        for v in vertices(graph)
            v_orig = vertices(original)[findfirst(v_orig -> v.data == v_orig.data, vertices(original))]
            for (e, e_orig) in zip(in_edges(v, graph), in_edges(v_orig, original))
                @test e.data == e_orig.data
            end
            for (e, e_orig) in zip(out_edges(v, graph), out_edges(v_orig, original))
                @test e.data == e_orig.data
            end
        end
    end

    @testset "remove_edge!" begin
        Random.seed!(11)
        graph = DirectedGraph{Vertex{Int64}, Edge{Float64}}()
        for i = 1 : 100
            add_vertex!(graph, Vertex(i))
        end
        for i = 1 : num_vertices(graph) - 1
            add_edge!(graph, rand(vertices(graph)), rand(vertices(graph)), Edge(Float64(i)))
        end
        original = deepcopy(graph)

        edge = rand(edges(graph))
        remove_edge!(graph, edge)
        @test edge ∉ edges(graph)
        @test num_edges(graph) == num_edges(original) - 1
        for v in vertices(graph)
            @test edge ∉ in_edges(v, graph)
            @test edge ∉ out_edges(v, graph)
        end

        for e in edges(graph)
            e_orig = edges(original)[findfirst(e_orig -> e.data == e_orig.data, edges(original))]
            @test source(e, graph).data == source(e_orig, original).data
            @test target(e, graph).data == target(e_orig, original).data
        end
    end

    @testset "rewire!" begin
        Random.seed!(12)
        graph = DirectedGraph{Vertex{Int64}, Edge{Float64}}()
        for i = 1 : 100
            add_vertex!(graph, Vertex(i))
        end
        for i = 1 : num_vertices(graph) - 1
            add_edge!(graph, rand(vertices(graph)), rand(vertices(graph)), Edge(Float64(i)))
        end
        original = deepcopy(graph)

        edge = rand(edges(graph))
        oldsource = source(edge, graph)
        oldtarget = target(edge, graph)
        non_source_vertices = delete!(Set(vertices(graph)), oldsource)
        non_target_vertices = delete!(Set(vertices(graph)), oldtarget)
        newsource = rand(collect(non_source_vertices))
        newtarget = rand(collect(non_target_vertices))

        rewire!(graph, edge, newsource, newtarget)
        @test source(edge, graph) == newsource
        @test target(edge, graph) == newtarget
        @test edge ∈ out_edges(newsource, graph)
        @test edge ∈ in_edges(newtarget, graph)
        @test edge ∉ out_edges(oldsource, graph)
        @test edge ∉ in_edges(oldtarget, graph)

        @test map(x -> x.data, vertices(original)) == map(x -> x.data, vertices(graph))
        for e in filter(e -> e != edge, edges(graph))
            e_orig = edges(original)[findfirst(e_orig -> e.data == e_orig.data, edges(original))]
            @test source(e, graph).data == source(e_orig, original).data
            @test target(e, graph).data == target(e_orig, original).data
        end
    end

    @testset "replace_edge!" begin
        Random.seed!(13)
        graph = DirectedGraph{Vertex{Int64}, Edge{Float64}}()
        for i = 1 : 100
            add_vertex!(graph, Vertex(i))
        end
        for i = 1 : num_vertices(graph) - 1
            add_edge!(graph, rand(vertices(graph)), rand(vertices(graph)), Edge(Float64(i)))
        end
        original = deepcopy(graph)

        for i = 1 : 10
            old_edge = rand(edges(graph))
            src = source(old_edge, graph)
            dest = target(old_edge, graph)
            new_edge = Edge(NaN)
            replace_edge!(graph, old_edge, new_edge)

            @test Graphs.edge_index(old_edge) == -1
            @test all(map(x -> x.data, vertices(graph)) .== map(x -> x.data, vertices(original)))
            @test new_edge ∈ in_edges(dest, graph)
            @test old_edge ∉ in_edges(dest, graph)
            @test new_edge ∈ out_edges(src, graph)
            @test old_edge ∉ out_edges(src, graph)
            @test source(new_edge, graph) == src
            @test target(new_edge, graph) == dest
            @test isnan(new_edge.data)
        end
    end

    @testset "SpanningTree" begin
        Random.seed!(14)
        rootdata = 0

        # graph1: tree grown incrementally
        graph1 = DirectedGraph{Vertex{Int64}, Edge{Int32}}()
        root1 = Vertex(rootdata)
        add_vertex!(graph1, root1)
        tree1 = SpanningTree(graph1, root1)

        # graph2: tree constructed after graph is built
        graph2 = DirectedGraph{Vertex{Int64}, Edge{Int32}}()
        root2 = Vertex(rootdata)
        add_vertex!(graph2, root2)

        nedges = 15
        for i = 1 : nedges
            parentind = rand(1 : num_vertices(graph1))
            childdata = i
            edgedata = Int32(i + 3)
            add_edge!(tree1, vertices(tree1)[parentind], Vertex(childdata), Edge(edgedata))
            add_edge!(graph2, vertices(graph2)[parentind], Vertex(childdata), Edge(edgedata))
        end

        tree2 = SpanningTree(graph2, root2)

        @test all(map(x -> x.data, vertices(tree1)) == map(x -> x.data, vertices(tree2)))
        for (v1, v2) in zip(vertices(tree1), vertices(tree2))
            if v1 == root(tree1)
                @test v2 == root(tree2)
            else
                @test edge_to_parent(v1, tree1).data == edge_to_parent(v2, tree2).data

                outedgedata1 = map(x -> x.data, collect(edges_to_children(v1, tree1)))
                outedgedata2 = map(x -> x.data, collect(edges_to_children(v2, tree2)))
                @test isempty(setdiff(outedgedata1, outedgedata2))

                @test isempty(setdiff(out_edges(v1, graph1), edges_to_children(v1, tree1)))
                @test isempty(setdiff(out_edges(v2, graph2), edges_to_children(v2, tree2)))
            end
        end

        tree = tree1
        show(devnull, tree)

        @test_throws AssertionError add_edge!(tree, rand(vertices(tree)), rand(vertices(tree)), Edge(rand(Int32)))

        for src in vertices(tree)
            for dest in vertices(tree)
                src_ancestors = ancestors(src, tree)
                dest_ancestors = ancestors(dest, tree)
                lca = lowest_common_ancestor(src, dest, tree)
                p = TreePath(src, dest, tree)

                show(devnull, p)
                @inferred collect(p)

                @test source(p) == src
                @test target(p) == dest

                source_to_lca = collect(edge for edge in p if direction(edge, p) == PathDirections.up)
                target_to_lca = reverse!(collect(edge for edge in p if direction(edge, p) == PathDirections.down))

                for (v, v_ancestors, pathsegment) in [(src, src_ancestors, source_to_lca); (dest, dest_ancestors, target_to_lca)]
                    if v == root(tree)
                        @test tree_index(v, tree) == 1
                        @test lca == v
                        @test length(v_ancestors) == 1
                        @test isempty(pathsegment)
                    end
                    for v_ancestor in v_ancestors
                        @test tree_index(v_ancestor, tree) <= tree_index(v, tree)
                    end
                    @test lca ∈ v_ancestors
                    if v != lca
                        @test source(last(pathsegment), tree) == lca
                    end
                end

                for v in intersect(src_ancestors, dest_ancestors)
                    @test tree_index(v, tree) <= tree_index(lca, tree)
                    if v != lca
                        @test v ∉ (v -> source(v, tree)).(source_to_lca)
                        @test v ∉ (v -> source(v, tree)).(target_to_lca)
                    end
                end
                for v in setdiff(src_ancestors, dest_ancestors)
                    @test tree_index(v, tree) > tree_index(lca, tree)
                end
            end
        end

        for i = 1 : 10
            old_edge = rand(edges(tree))

            src = source(old_edge, tree)
            dest = target(old_edge, tree)
            old_id = Graphs.edge_index(old_edge)

            # make sure that replacing edge with itself doesn't mess with anything
            replace_edge!(tree, old_edge, old_edge)
            @test source(old_edge, tree) == src
            @test target(old_edge, tree) == dest
            @test Graphs.edge_index(old_edge) == old_id

            # replace with a new edge
            d = Int32(-10 * i)
            new_edge = Edge(d)
            replace_edge!(tree, old_edge, new_edge)

            @test Graphs.edge_index(old_edge) == -1
            @test new_edge ∈ in_edges(dest, tree)
            @test old_edge ∉ in_edges(dest, tree)
            @test new_edge ∈ out_edges(src, tree)
            @test old_edge ∉ out_edges(src, tree)
            @test source(new_edge, tree) == src
            @test target(new_edge, tree) == dest
            @test new_edge.data == d
            @test edge_to_parent(dest, tree) == new_edge
            @test new_edge ∈ edges_to_children(src, tree)
        end

        original = deepcopy(graph2)
        edgemap = Dict(zip(edges(original), edges(graph2)))
        newroot = rand(setdiff(vertices(graph2), [root2]))
        flipped_edge_map = Dict{Edge{Int32}, Edge{Int32}}()
        newtree = SpanningTree(graph2, newroot, flipped_edge_map)

        @test !isempty(flipped_edge_map)
        for (oldedge, newedge) in edgemap
            flipped = haskey(flipped_edge_map, newedge)
            if flipped
                newedge = flipped_edge_map[newedge]
            end

            old_source_ind = Graphs.vertex_index(source(oldedge, original))
            old_target_ind = Graphs.vertex_index(target(oldedge, original))
            new_source_ind = Graphs.vertex_index(source(newedge, graph2))
            new_target_ind = Graphs.vertex_index(target(newedge, graph2))

            if flipped
                @test oldedge.data == -newedge.data
                @test new_source_ind == old_target_ind
                @test new_target_ind == old_source_ind
            else
                @test oldedge.data == newedge.data
                @test new_source_ind == old_source_ind
                @test new_target_ind == old_target_ind
            end
        end
    end

    @testset "reindex!" begin
        Random.seed!(15)
        graph = DirectedGraph{Vertex{Int64}, Edge{Float64}}()
        for i = 1 : 100
            add_vertex!(graph, Vertex(i))
        end
        for i = 1 : num_vertices(graph) - 1
            add_edge!(graph, rand(vertices(graph)), rand(vertices(graph)), Edge(Float64(i)))
        end
        newvertices = shuffle(vertices(graph))
        newedges = shuffle(edges(graph))
        reindex!(graph, newvertices, newedges)
        @test all(vertices(graph) .== newvertices)
        @test all(edges(graph) .== newedges)
        @test all(Graphs.vertex_index.(vertices(graph)) .== 1 : num_vertices(graph))
        @test all(Graphs.edge_index.(edges(graph)) .== 1 : num_edges(graph))
    end

    @testset "map-like constructor" begin
        Random.seed!(16)
        graph = DirectedGraph{Vertex{Int32}, Edge{Float32}}()
        for i = Int32(1) : Int32(100)
            add_vertex!(graph, Vertex(i))
        end
        for i = 1 : num_vertices(graph) - 1
            add_edge!(graph, rand(vertices(graph)), rand(vertices(graph)), Edge(Float32(i)))
        end
        mappedgraph = DirectedGraph(x -> Vertex(Int64(x.data), x.id), x -> Edge(Float64(x.data), x.id), graph)

        @test vertextype(mappedgraph) == Vertex{Int64}
        @test edgetype(mappedgraph) == Edge{Float64}
        @test all(v1.data == v2.data for (v1, v2) in zip(vertices(graph), vertices(mappedgraph)))
        @test all(e1.data == e2.data for (e1, e2) in zip(edges(graph), edges(mappedgraph)))
        @test all(source(e1, graph).data == source(e2, mappedgraph).data for (e1, e2) in zip(edges(graph), edges(mappedgraph)))
        @test all(target(e1, graph).data == target(e2, mappedgraph).data for (e1, e2) in zip(edges(graph), edges(mappedgraph)))

        inedgesmatch(v1, v2) = all(e1.data == e2.data for (e1, e2) in zip(in_edges(v1, graph), in_edges(v2, mappedgraph)))
        @test all(inedgesmatch(v1, v2) for (v1, v2) in zip(vertices(graph), vertices(mappedgraph)))

        outedgesmatch(v1, v2) = all(e1.data == e2.data for (e1, e2) in zip(out_edges(v1, graph), out_edges(v2, mappedgraph)))
        @test all(outedgesmatch(v1, v2) for (v1, v2) in zip(vertices(graph), vertices(mappedgraph)))
    end
end


text:
function extract_set!(
    sets::Dict{String,TV}, setname::AbstractString, file, offset=0
) where {TI,TV<:AbstractVector{TI}}
    sets[setname] = Int[]
    vector = sets[setname]

    pattern_single = r"^(\d+)"
    pattern_subset = r"^([^,]+)"
    line = readline(file)
    m = match(stopping_pattern, line)
    while m isa Nothing
        m = match(pattern_single, line)
        if m != nothing
            push!(vector, parse(TI, m[1]) - offset)
        else
            m = match(pattern_subset, line)
            if m != nothing
                subsetname = String(m[1])
                if haskey(sets, subsetname)
                    append!(vector, sets[subsetname])
                end
            end
        end
        line = readline(file)
        m = match(stopping_pattern, line)
    end
    return line
end


text:
"""
    SEPP

Abstract type for the self-exciting point process model.
"""
abstract type SEPP end

"""
    baseline(sepp::SEPP)::Real

Get the baseline of an SEPP.
"""
baseline(sepp::SEPP)::Real = sepp.μ

"""
    exp_decay(sepp::SEPP)::Real

Get the exponential decay of an SEPP.

For now all SEPP models in that package have an exponential kernel, hence a method for any SEPP.
"""
exp_decay(sepp::SEPP)::Real = sepp.γ

"""
    self_excitement_factor(sepp::SEPP)::Real

Get the factor by which the volatility function gets multiplied before it is added to the baseline in the rate.
"""
self_excitement_factor(sepp::SEPP)::Real = sepp.ϕ

"""
    params(sepp::SEPP)::AbstractDict

Get the paramaters of the model.
"""
function params end

include("SEPP/SEPPExpKern.jl")
include("SEPP/SEMPPExpKern.jl")


text:
function test_solmat(model;
    ϵ=1.0e-6,
    tsim=0.1,
    ctrl=(m, k)->nothing,
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    verbose=false,
    T=Float64,
    kwargs...)

    @testset "$(string(model))" begin
        # mechanism
        mechanism = get_mechanism(model,
            timestep=timestep,
            gravity=gravity;
            kwargs...)
        initialize!(mechanism, model)

        # simulate
        storage = simulate!(mechanism, tsim, ctrl,
            record=true,
            verbose=verbose,
            opts=SolverOptions(rtol=ϵ, btol=ϵ))

        # Set data
        Nb = length(mechanism.bodies)
        data = Dojo.get_data(mechanism)
        Dojo.set_data!(mechanism, data)
        sol = Dojo.get_solution(mechanism)
        attjac = Dojo.attitude_jacobian(data, Nb)

        # IFT
        solmat = Dojo.full_matrix(mechanism.system)
        # finite diff
        fd_solmat = finite_difference_solution_matrix(mechanism, data, sol,
            δ=1.0e-5,
            verbose=verbose)
        @test norm(fd_solmat + solmat, Inf) < ϵ
    end
    return nothing
end

function finite_difference_solution_matrix(mechanism::Mechanism, data::AbstractVector, sol::AbstractVector;
    δ=1.0e-8,
    verbose=false)

    nsol = length(sol)
    jac = zeros(nsol, nsol)

    Dojo.set_data!(mechanism, data)
    Dojo.set_solution!(mechanism, sol)

    for i = 1:nsol
        verbose && println("$i / $nsol")
        solp = deepcopy(sol)
        solm = deepcopy(sol)
        solp[i] += δ
        solm[i] -= δ
        rp = Dojo.evaluate_residual!(deepcopy(mechanism), data, solp)
        rm = Dojo.evaluate_residual!(deepcopy(mechanism), data, solm)
        jac[:, i] = (rp - rm) / (2δ)
    end
    return jac
end

function control!(mechanism, k;
    u=0.1)
    for joint in mechanism.joints
        nu = Dojo.input_dimension(joint,
            ignore_floating_base=false)
        su = mechanism.timestep * u * sones(nu)
        Dojo.set_input!(joint, su)
    end
end

#TODO: decrease ϵ tolerance once we replace finite-difference methods
# Flying after 0.1 sec simulation
@testset "Flying" begin
    tsim = 0.1
    test_solmat(:atlas,      tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7)
    test_solmat(:atlas,      tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7)
    test_solmat(:atlas,      tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1e3, damper=5e2)
    test_solmat(:quadruped,  tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:block,        tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7)
    test_solmat(:snake,      tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:slider,     tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:pendulum,   tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:npendulum,  tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:nslider,    tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:twister,    tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:sphere,     tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, contact_type=:nonlinear)
    test_solmat(:sphere,     tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, contact_type=:linear)
    test_solmat(:sphere,     tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, contact_type=:impact)
    test_solmat(:ant,        tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7)
    test_solmat(:halfcheetah,tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7)
end

# In contact with the ground after 0.4 sec simulation
@testset "In contact" begin
    tsim = 0.4
    test_solmat(:atlas,      tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7)
    test_solmat(:atlas,      tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1e3, damper=5e2)
    test_solmat(:quadruped,  tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:block,        tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7)
    test_solmat(:snake,      tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:slider,     tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:pendulum,   tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:npendulum,  tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:nslider,    tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:twister,    tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, spring=1.0, damper=0.2)
    test_solmat(:sphere,     tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7)
    test_solmat(:sphere,     tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, contact_type=:nonlinear)
    test_solmat(:sphere,     tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, contact_type=:linear)
    test_solmat(:sphere,     tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7, contact_type=:impact)
    test_solmat(:ant,        tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7)
    test_solmat(:halfcheetah,tsim=tsim, ctrl=(m,k)->control!(m,k,u=0.1), ϵ=1.0e-7)
end


text: # IO functions for saving/reading native data


text:
using MOCNeutronTransport
using Plots
pyplot()
T = Float64
Nsoln = 20
Nset = [ 1, 2, 3, 4, 5, 10, 15, 20]
f = area

p₁ = Point_3D(T, 0)
p₂ = Point_3D(T, 2)
p₃ = Point_3D(T, 2, 3)
p₄ = Point_3D(T, 0, 3)
p₅ = Point_3D(T, 3//2, 1//2)
p₆ = Point_3D(T, 5//2, 3//2)
p₇ = Point_3D(T, 3//2, 5//2)
p₈ = Point_3D(T, 0,    3//2)
p₈ = Point_3D(T, 0, 3//2, 2)
quad8 = Quadrilateral8_3D((p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈))

soln = f(quad8; N=Nsoln)
err = [abs(soln - f(quad8; N = n)) for n = Nset]
plot(Nset[1:length(Nset)-1] .* Nset[1:length(Nset)-1],
     err[1:length(Nset)-1],
     title = "Quadrature Points vs Quadrilateral8_3D area error (Nsoln = $Nsoln)",
     yaxis=:log
    )
xlabel!("Quadrature Points")
ylabel!("area error")

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport