text:
let
    #= line 1 =#
    var"|Prim|b_op_Getter"::Function = function op_Getter(x)
                #= none:1 =#
                #= none:4 =#
                function (a,)
                    #= none:4 =#
                    #= none:7 =#
                    x(a)
                end
            end::Function
    var"|Prim|c_unit"::Symbol = Symbol("{}")::Symbol
    var"|Prim|d_Vec"::Symbol = Symbol("Vec@10_|Prim|e")::Symbol
    var"|Prim|f_make"::Function = function (var"|Prim|g_t",)
                #= prim.mlfs:14 =#
                var"|Prim|g_t"::Function
            end::Function
    var"|Prim|h_mkVec"::Function = (_->begin
                    #= none:1 =#
                    #= none:3 =#
                    []
                end)::Function
    var"|Prim|i_pushVec"::Any = (x->begin
                    #= none:1 =#
                    #= none:3 =#
                    function (a,)
                        #= none:3 =#
                        #= none:6 =#
                        push!(x, a)
                        #= none:8 =#
                        ()
                    end
                end)::Any
    var"|Prim|j_println"::Function = function (x,)
                #= none:1 =#
                #= none:4 =#
                println(x)
                #= none:6 =#
                ()
            end::Function
    var"|Prim|k"::Function = function (_,)
            #= prim.mlfsf:1 =#
            var"|Prim|j_println"
        end
    var"|Prim|l"::Function = function (_,)
            #= prim.mlfsf:1 =#
            var"|Prim|i_pushVec"
        end
    var"|Prim|m"::Function = function (_,)
            #= prim.mlfsf:1 =#
            var"|Prim|h_mkVec"
        end
    var"|Prim|n"::Function = function (_,)
            #= prim.mlfsf:1 =#
            var"|Prim|f_make"
        end
    var"|Prim|o"::Function = function (_,)
            #= prim.mlfsf:1 =#
            Symbol("Vec@10_|Prim|e")
        end
    var"|Prim|p"::Function = function (_,)
            #= prim.mlfsf:1 =#
            Symbol("{}")
        end
    var"|Prim|q"::Function = function (_,)
            #= prim.mlfsf:1 =#
            var"|Prim|b_op_Getter"
        end
    var"|Prim|"::Any = Symbol("Module Prim")
    Prim::Any = Symbol("namespace {Module Prim, {{make, |Prim|f_make, forall t@16162419886138674179 a@163180927404684273.implicit[(Type t@16162419886138674179)->t@16162419886138674179 a@163180927404684273]->(Type t@16162419886138674179)->t@16162419886138674179 a@163180927404684273}, {unit, |Prim|c_unit, Type {}}, {mkVec, |Prim|h_mkVec, implicit[forall a@13733476492383310116.(Type Vec@10_|Prim|e)->Vec@10_|Prim|e a@13733476492383310116]}, {println, |Prim|j_println, forall a@13197780277297027335.a@13197780277297027335->{}}, {op_Getter, |Prim|b_op_Getter, forall a@3636708016426571584 b@17888154148990276253 c@17426973612936557375.(field {a@3636708016426571584, b@17888154148990276253, c@17426973612936557375})->a@3636708016426571584->c@17426973612936557375}, {Vec, |Prim|d_Vec, Type Vec@10_|Prim|e}, {pushVec, |Prim|i_pushVec, implicit[forall a@4941567182984758855.field {Vec@10_|Prim|e a@4941567182984758855, push, a@4941567182984758855->{}}]}}}")
    (var"|Prim|j_println"::Function)(1::Int64)::Tuple
    (var"|Prim|j_println"::Function)(:i64::Symbol)::Tuple
    var"|main|r_a"::Any = (var"|Prim|f_make"(var"|Prim|h_mkVec"::Function)::Function)(Symbol("Vec@10_|Prim|e")::Symbol)::Any
    var"|main|s_c"::Tuple = ((op_Getter(var"|Prim|i_pushVec"::Any))(var"|main|r_a"::Any)::Function)(1::Int64)::Tuple
    var"|main|t_mk2"::Function = function (var"|main|u_x",)
                #= main.mlfs:17 =#
                (var"|main|u_x"::Any, var"|main|u_x"::Any)::Tuple
            end::Function
    (var"|Prim|j_println"::Function)((var"|main|t_mk2"::Function)(var"|main|r_a"::Any)::Tuple)::Tuple
    var"|main|v_index"::Function = (x::Vector->begin
                    #= none:1 =#
                    #= none:3 =#
                    i::Int64->begin
                            #= none:3 =#
                            #= none:5 =#
                            x[i]
                        end
                end)::Function
    var"|main|w_CanAdd1"::Symbol = Symbol("CanAdd1@24_|main|x")::Symbol
    var"|main|y_int_add1"::Function = (_->begin
                    #= none:1 =#
                    #= none:3 =#
                    x::Int64->begin
                            #= none:3 =#
                            #= none:5 =#
                            x + 1
                        end
                end)::Function
    var"|main|z_float_add1"::Function = (_->begin
                    #= none:1 =#
                    #= none:3 =#
                    x::Float64->begin
                            #= none:3 =#
                            #= none:5 =#
                            x + 1.0
                        end
                end)::Function
    var"|main|ab_vec_head_add1"::Function = function (var"|main|bb_instance",)
                #= main.mlfs:35 =#
                function (var"|main|cb_x",)
                        #= main.mlfs:35 =#
                        ((var"|main|bb_instance"::Function)(Symbol("CanAdd1@24_|main|x")::Symbol)::Function)(((var"|main|v_index"::Function)(var"|main|cb_x"::Any)::Function)(1::Int64)::Any)::Any
                    end::Function
            end::Function
    (var"|Prim|j_println"::Function)((var"|main|ab_vec_head_add1"(var"|main|y_int_add1"::Function)::Function)(var"|main|r_a"::Any)::Int64)::Tuple
    var"|main|db"::Function = function (_,)
            #= main.mlfsf:1 =#
            var"|main|ab_vec_head_add1"
        end
    var"|main|eb"::Function = function (_,)
            #= main.mlfsf:1 =#
            var"|main|z_float_add1"
        end
    var"|main|fb"::Function = function (_,)
            #= main.mlfsf:1 =#
            var"|main|y_int_add1"
        end
    var"|main|gb"::Function = function (_,)
            #= main.mlfsf:1 =#
            Symbol("CanAdd1@24_|main|x")
        end
    var"|main|hb"::Function = function (_,)
            #= main.mlfsf:1 =#
            var"|main|v_index"
        end
    var"|main|ib"::Function = function (_,)
            #= main.mlfsf:1 =#
            var"|main|t_mk2"
        end
    var"|main|jb"::Function = function (_,)
            #= main.mlfsf:1 =#
            var"|main|s_c"
        end
    var"|main|kb"::Function = function (_,)
            #= main.mlfsf:1 =#
            var"|main|r_a"
        end
    var"|main|lb"::Function = function (_,)
            #= main.mlfsf:1 =#
            var"|Prim|i_pushVec"
        end
    var"|main|mb"::Function = function (_,)
            #= main.mlfsf:1 =#
            Symbol("Vec@10_|Prim|e")
        end
    var"|main|nb"::Function = function (_,)
            #= main.mlfsf:1 =#
            var"|Prim|b_op_Getter"
        end
    var"|main|ob"::Function = function (_,)
            #= main.mlfsf:1 =#
            var"|Prim|j_println"
        end
    var"|main|pb"::Function = function (_,)
            #= main.mlfsf:1 =#
            var"|Prim|h_mkVec"
        end
    var"|main|qb"::Function = function (_,)
            #= main.mlfsf:1 =#
            Symbol("{}")
        end
    var"|main|rb"::Function = function (_,)
            #= main.mlfsf:1 =#
            var"|Prim|f_make"
        end
    var"|main|sb"::Function = function (_,)
            #= main.mlfsf:1 =#
            var"|Prim|"
        end
    var"|main|"::Any = Symbol("Module main")
    main::Any = Symbol("namespace {Module main, {{Prim, |Prim|, Module Prim}, {make, |Prim|f_make, forall t@16162419886138674179 a@163180927404684273.implicit[(Type t@16162419886138674179)->t@16162419886138674179 a@163180927404684273]->(Type t@16162419886138674179)->t@16162419886138674179 a@163180927404684273}, {c, |main|s_c, {}}, {mkVec, |Prim|h_mkVec, implicit[forall a@13733476492383310116.(Type Vec@10_|Prim|e)->Vec@10_|Prim|e a@13733476492383310116]}, {int_add1, |main|y_int_add1, implicit[(Type CanAdd1@24_|main|x)->i64->i64]}, {op_Getter, |Prim|b_op_Getter, forall a@3636708016426571584 b@17888154148990276253 c@17426973612936557375.(field {a@3636708016426571584, b@17888154148990276253, c@17426973612936557375})->a@3636708016426571584->c@17426973612936557375}, {float_add1, |main|z_float_add1, implicit[(Type CanAdd1@24_|main|x)->f64->f64]}, {CanAdd1, |main|w_CanAdd1, Type CanAdd1@24_|main|x}, {index, |main|v_index, forall a@5978395343969346419.(Vec@10_|Prim|e a@5978395343969346419)->i64->a@5978395343969346419}, {unit, |Prim|c_unit, Type {}}, {println, |Prim|j_println, forall a@13197780277297027335.a@13197780277297027335->{}}, {mk2, |main|t_mk2, forall a@5152008298992179848.a@5152008298992179848->{a@5152008298992179848, a@5152008298992179848}}, {a, |main|r_a, Vec@10_|Prim|e i64}, {vec_head_add1, |main|ab_vec_head_add1, forall a@7907858595278931228.implicit[(Type CanAdd1@24_|main|x)->a@7907858595278931228->a@7907858595278931228]->(Vec@10_|Prim|e a@7907858595278931228)->a@7907858595278931228}, {Vec, |Prim|d_Vec, Type Vec@10_|Prim|e}, {pushVec, |Prim|i_pushVec, implicit[forall a@4941567182984758855.field {Vec@10_|Prim|e a@4941567182984758855, push, a@4941567182984758855->{}}]}}}")
    0
end


text:
immutable RectType
    x::Cint
    y::Cint
    width::Cint
    height::Cint
end

immutable ClassificationType
    id::Cint
    confidence::Cfloat
end

immutable CompType
    rect::RectType
    neighbors::Cint
    classification::ClassificationType
end

immutable CCVSize
    width::Cint
    height::Cint
end

ext: jl


text:
import Base: convert

for (V, W) in [(DVertex, IVertex), (IVertex, DVertex)]
  @eval function convert(::Type{$W{T}}, v::$V, cache = ODict()) where T
    haskey(cache, v) && return cache[v]
    w = cache[v] = $W{T}(value(v))
    thread!(w, [convert($W{T}, v′, cache) for v′ in inputs(v)]...)
  end
  @eval convert(::Type{$W}, v::Vertex) = convert($W{eltype(v)}, v)
end


text:
"""
    Trace(stat, b=500, f=value)

Wrapper around an OnlineStat that stores `b` "snapshots" (a fixed copy of the OnlineStat).  The snapshots
are taken at approximately equally-spaced intervals between 0 and the current `nobs`.  The main use
case is visualizing state changes as observations are added.

# Example

    using OnlineStats, Plots

    o = Trace(Mean(), 10)

    fit!(o, 1:100)

    OnlineStats.snapshots(o)

    plot(o)
"""
mutable struct Trace{T, O <: OnlineStat{T}} <: OnlineStat{T}
    parts::Vector{Pair{Tuple{Int,Int}, O}}
    b::Int
    n::Int
end
Trace(o::OnlineStat, b::Int=100) = Trace([(1,1) => o], b, 0)

name(o::Trace, args...) = "Trace($(name(o.parts[end][2], args...)))"
value(o::Trace) = value(o.parts[end][2])

snapshots(o::Trace) = last.(o.parts)

function _fit!(o::Trace, y)
    n = o.n += 1
    parts = o.parts
    (a, b), stat = parts[end]
    a ≤ n ≤ b ?
        _fit!(stat, y) :
        push!(parts, (n, n + ceil(Int, n / o.b) - 1) => fit!(copy(stat), y))
    if length(parts) > o.b && nobs(parts[end][2]) == parts[end][1][2]
        trace_merge_next!(parts)
    end
end

function trace_merge_next!(parts)
    n = parts[1][1][2] - parts[1][1][1]
    i = 1
    for (j, p) in enumerate(parts)
        a, b = p[1]
        b - a < n && (i = j; break;)
    end
    a, b = parts[i], parts[i+1]
    parts[i] = (min(a[1][1], b[1][1]), max(a[1][2], b[1][2])) => b[2]
    deleteat!(parts, i + 1)
end


text:
module DiffEqUncertainty

using DiffEqBase, Statistics, Distributions, Reexport
@reexport using Quadrature
using KernelDensity

include("probints.jl")
include("koopman.jl")

# Type Piracy, should upstream
Base.eltype(K::UnivariateKDE)  = eltype(K.density)
Base.minimum(K::UnivariateKDE) = minimum(K.x)
Base.maximum(K::UnivariateKDE) = maximum(K.x)

export ProbIntsUncertainty,AdaptiveProbIntsUncertainty
export expectation, centralmoment, Koopman, MonteCarlo

end


text:
import Base: parse

"""
Abstract type for Command Line Interface Patterns
"""
abstract type CLIPattern end

"""
Command Line Interface Flags, e.g --flag, -f
"""
abstract type CLIFlag <: CLIPattern end

==(lhs::CLIFlag, rhs::CLIFlag) = (lhs.name == rhs.name) && (lhs.val == rhs.val)

"""
    CLILongFlag(name::String, val::String)

Command line long flags, startswith `--`. e.g `--flag`
"""
struct CLILongFlag <: CLIFlag
    name::String
    val::String
end

"""
    CLIShortFlag(name::Char, val::String)

Command line short flags, startswith `-`. e.g `-f`
"""
struct CLIShortFlag <: CLIFlag
    name::Char
    val::String
end

pattern(::Type{CLILongFlag}) = r"^--([a-zA-Z\-_]+[0-9]*)"
pattern(::Type{CLIShortFlag}) = r"^-([a-zA-Z]{1})"

islongflag(s::AbstractString) = ismatch(pattern(CLILongFlag), s)
isshortflag(s::AbstractString) = ismatch(pattern(CLIShortFlag), s)
iscmdarg(s::AbstractString) = !islongflag(s) && !isshortflag(s)

function parse(::Type{CLILongFlag}, s::AbstractString)
    m = match(pattern(CLILongFlag), s)
    m === nothing && return nothing

    name = m.captures[1]
    s = s[length(m.match)+1:end]
    m = match(r"^=(.+)", s)

    if m === nothing
        return CLILongFlag(name, "")
    else
        return CLILongFlag(name, m.captures[1])
    end
end

function parse(::Type{CLIShortFlag}, s::AbstractString)
    m = match(pattern(CLIShortFlag), s)
    m === nothing && return nothing

    name = m.captures[1]
    s = s[length(m.match)+1:end]

    m = match(r"^=(.+)", s)
    if m != nothing
        return CLIShortFlag(name[1], m.captures[1])
    end

    m = match(r"^([^=]+)", s)

    if m === nothing
        return CLIShortFlag(name[1], "")
    else
        return CLIShortFlag(name[1], m.captures[1])
    end
end

"""
get flags from an command line argument stream.

if greedy is `false`, it stops at the first non-flag
argument, or it reads all the flags and put ptr in
stream to `eof`.
"""
function getflags(stream::ARGStream; greedy=false)
    flags = []
    while !eof(stream)
        pos = position(stream)
        arg = readarg(stream)
        flag = parse(CLILongFlag, arg)
        if flag === nothing
            flag = parse(CLIShortFlag, arg)
        end

        if flag === nothing
            if !greedy
                seek(stream, pos)
                break
            end
        else
            push!(flags, flag)
        end
    end
    return flags
end

"""
get command arguments from stream. if greedy is `true`
it reads all possible arguments.
"""
function getcmds(stream::ARGStream; greedy=false)
    args = String[]
    while !eof(stream)
        pos = position(stream)
        arg = readarg(stream)

        if iscmdarg(arg)
            push!(args, arg)
        else
            greedy || return args
        end
    end
    return args
end

"""
find corresponding option from the list according to command line flags
"""
function findoption(clioption::CLIFlag, options::Vector{T}) where {T <: AbstractOption}
    for each in options
        if clioption.name == each
            return each
        end
    end
    return nothing
end

"""
eval and parse all arguments from command line according to signature
"""
function getargs(cmds::Vector{String}, sigs::Vector{DataType})
    if length(cmds) != length(sigs)
        throw(CLIError("Expect $(length(sigs)) positional arguments (got $(length(cmds)))"))
    end

    args = []
    for (cmd, sig) in zip(cmds, sigs)
        val = eval(parse(sig, cmd))
        isa(val, sig) || throw(CLIArgError(sig, typeof(val)))
        push!(args, val)
    end
    return args
end

"""
eval and parse all keywords according to command line flags
"""
function getkeywords(flags, options)
    keywords = []
    for each in flags
        option = findoption(each, options)
        if option === nothing
            throw(OptionError(each.name))
        else
            push!(keywords, parse(each.val, option))
        end
    end
    return keywords
end

"""
flush the error infomation to IO if the function throws. rethrow the error
if `raise` is true (default).
"""
function flusherror(f::Function, io::IO, cmd::AbstractCommand; raise=true)
    result = try
        f()
    catch e
        if isa(e, CLIException)
            push!(e, name(cmd))
        end
        raise ? rethrow(e) : errorinfo(io, e)
        !(raise) && exit(-1)
    end
    return result
end

export resolve

"""
resolve command line inputs according to registered command. throw
all the errors if `debug` is `true`.
"""
function resolve(stream::ARGStream, cmd::MainCommand; debug=false)
    flags = getflags(stream)

    flusherror(STDERR, cmd; raise=debug) do
        keywords = getkeywords(flags, cmd.options)
        isempty(keywords) ? nothing : cmd(;keywords...)
    end

    flusherror(STDERR, cmd; raise=debug) do
        resolve(stream, cmd.cmds)
    end
end

function resolve(stream::ARGStream, cmd::NodeCommand)
    flags = getflags(stream)

    flusherror(STDERR, cmd) do
        keywords = getkeywords(flags, [help])
        cmd(;keywords...)
        resolve(stream, cmd.cmds)
    end
end

function resolve(stream::ARGStream, cmd::LeafCommand)
    start = position(stream)
    flags = getflags(stream;greedy=true)
    seek(stream, start)
    cmds = getcmds(stream;greedy=true)

    flusherror(STDERR, cmd) do
        keywords = getkeywords(flags, cmd.sig.options)
        args = getargs(cmds, cmd.sig.args)
        result = cmd(args...;keywords...)
        result === nothing ? nothing: println(STDOUT, result)
    end
end

function resolve(stream::ARGStream, cmds::Vector{AbstractCommand})
    eof(stream) && return 0

    arg = readarg(stream)
    subcmd = findcmd(arg, cmds)
    if subcmd === nothing
        throw(CmdError(arg))
    else
        resolve(stream, subcmd)
    end
    return 0
end


text:
module VDPTag2MCTSVisualize

using Plots, ContinuumWorld, POMDPToolbox, MCTS, VDPTag2

export VDPTag2MCTSVis

struct VDPTag2MCTSVis
    mdp::VDPTagMDP
    s::TagState
    actions::Vector{Vector{Float64}}
    next_states::Vector{Vector{TagState}}
    ns::Vector{Int}
    radii::Vector{Vector{Float64}}
end
VDPTag2MCTSVis(mdp::VDPTagMDP, s::TagState=TagState()) = VDPTag2MCTSVis(mdp, s, Vector{Float64}[],
    Vector{TagState}[], Int[], Vector{Float64}[])

function circle(x,y,r)
    c = Shape(Plots.partialcircle(0,2π,20,r))
    translate!(c, x, y)
    c
end

@recipe function f(vis::VDPTag2MCTSVis, h::SimHistory, index::Int)
    mdp = vis.mdp
    aspect_ratio := :equal
    xlim --> (-5, 5)
    ylim --> (-5, 5)
    @series begin
        label := "actions"
        alpha := 0.25
        rng=MersenneTwister(0)
        [translate!(Shape(Plots.partialcircle(ang-r, ang+r, 20, 1.0+0.4rand(rng)-0.2)),vis.s.agent...) for (ang,r) in zip(vis.actions[index],vis.radii[index]) if r < Inf]
    end
    @series begin
        label := "path"
        x = [s.agent[1] for s in state_hist(h)[1:end-1]]
        y = [s.agent[2] for s in state_hist(h)[1:end-1]]
        x, y
    end
    @series begin
        a = action_hist(h)[end]
        if a isa TagAction && a.look
            color := :blue
        else
            color := :red
        end
        s = state_hist(h)[end-1]
        label := "current agent position"
        pts = Plots.partialcircle(0, 2*pi, 100, mdp.tag_radius)
        x, y = Plots.unzip(pts)
        x+s.agent[1], y+s.agent[2]
    end
    @series begin
        seriestype := :scatter
        label := "current target"
        pos = state_hist(h)[end-1].target
        color --> :orange
        [pos[1]], [pos[2]]
    end
end
function MCTS.notify_listener(vis::VDPTag2MCTSVis,dsb::DSBPlanner,s,a,sp,r,snode,sanode,spnode)
    if s == vis.s
        tree = get(dsb.tree)
        sol = dsb.solver
        actions = [tree.a_labels[c] for c in tree.children[snode]]
        if isempty(vis.actions) || (actions != vis.actions[end])  #only push the deltas
            push!(vis.actions, actions)
            #push!(vis.next_states, [s+action for action in actions])
            push!(vis.ns, tree.total_n[snode])
            push!(vis.radii, [sol.r0_action/tree.total_n[snode]^sol.lambda_action for action in actions])
        end
    end
end
function MCTS.notify_listener(vis::VDPTag2MCTSVis,asb::ASBPlanner,s,a,sp,r,snode,sanode,spnode)
    if s == vis.s
        tree = get(asb.tree)
        sol = asb.solver
        actions = [tree.a_labels[c] for c in tree.children[snode]]
        if isempty(vis.actions) || (actions != vis.actions[end])  #only push the deltas
            push!(vis.actions, actions)
            #push!(vis.next_states, [s+action for action in actions])
            push!(vis.ns, tree.total_n[snode])
            push!(vis.radii, [tree.a_radius[x] for x in tree.children[snode]])
        end
    end
end
end # module


text:
st=Structure()
lcset=LoadCaseSet()
add_uniaxial_metal!(st,"steel",2e11,0.2,7849.0474)
add_general_section!(st,"frame",4.26e-3,3.301e-6,6.572e-5,9.651e-8,1e-3,1e-3,0,0)

N=600
for i in 0:N
    add_node!(st,i,3i/100,2i/100,2i/100)
end

for i in 0:N-1
    add_beam!(st,i,i,i+1,"steel","frame")
end

set_nodal_restraint!(st,0,true,true,true,true,true,true)

add_static_case!(lcset,"DL",1.1)
@time begin
assembly=assemble!(st,lcset,path=PATH)
end

@time begin
solve(assembly)
end

r=result_nodal_reaction(assembly,"DL",0)
@test r≈[-2.83076e-10,-9.04954e-10,8926.15,53556.9,-80335.3,-3.29473e-8] rtol=1e-2
r=result_nodal_displacement(assembly,"DL",N)
@test r≈[0.454025, 0.302683, -0.98385, -0.0336002, 0.0504003, -5.77725e-13] rtol=1e-2


text:
using MLDatasets
using Flux
using BSON
using Random
using Statistics
using Base.Iterators: partition
using Flux: crossentropy, onehotbatch, onecold

accuracy(x, y) = mean(onecold(cpu(m(x))) .== onecold(cpu(y)))

function reshape_data(X::AbstractArray{T, 3}, y::AbstractVector) where T
    s = size(X)
    return reshape(X, s[1], s[2], 1, s[3]), reshape(y, 1, :)
end

function train_model!(m, X, y;
    opt=ADAM(0.001),
    batch_size=128,
    n_epochs=10,
    file_name="")

    loss(x, y) = crossentropy(m(x), y)

    batches_train = map(partition(randperm(size(y, 2)), batch_size)) do inds
        return (gpu(X[:, :, :, inds]), gpu(y[:, inds]))
    end

    for i in 1:n_epochs
        println("Iteration " * string(i))
        Flux.train!(loss, params(m), batches_train, opt)
    end

    !isempty(file_name) && BSON.bson(file_name, m=m|>cpu)

    return
end

function load_data(dataset; onehot=false, T=Float32)
    classes = 0:9
    X_train, y_train = reshape_data(dataset.traindata(T)...)
    X_test, y_test = reshape_data(dataset.testdata(T)...)
    y_train = T.(y_train)
    y_test = T.(y_test)

    if onehot
        y_train = onehotbatch(y_train[:], classes)
        y_test = onehotbatch(y_test[:], classes)
    end

    return X_train, y_train, X_test, y_test
end

dataset = MLDatasets.MNIST
T = Float32
X_train, y_train, X_test, y_test = load_data(dataset; T=T, onehot=true)

m = Chain(
    Conv((2,2), 1=>16, sigmoid),
    MaxPool((2,2)),
    Conv((2,2), 16=>8, sigmoid),
    MaxPool((2,2)),
    flatten,
    Dense(288, size(y_train,1)), softmax) |> gpu

file_name = joinpath("data", "mnist_sigmoid.bson")
train_model!(m, X_train, y_train; file_name=file_name, n_epochs=100)

accuracy(X_test |> gpu, y_test |> gpu)

max_stars_repo_name: rzehumat/JuliaCourse

max_issues_repo_name: rzehumat/JuliaCourse

max_forks_repo_name: rzehumat/JuliaCourse


text:
Mathy = Union{ModelingToolkit.Operation, ModelingToolkit.Variable, Number}
#---------------------------------------------------------------
#---------------------------------------------------------------
# Addition of Differentials

function Base.:+(x::Differential, y::Mathy)
    x + Differential([DTag()], [y])
end

function Base.:+(x::Mathy, y::Differential)
    Differential([DTag()], [x]) + y
end

function Base.:+(x::Differential, y::Differential)
    tags = tagUnion(x,y)
    values = [((try x[tag] catch 0 end) + (try y[tag] catch 0 end)) for tag in tags]
    Differential(tags, values)
end

#---------------------------------------------------------------
#---------------------------------------------------------------
# Multiplication of Differentials

function Base.:*(x::Differential, y::Mathy)
    Differential(tag => value * y for (tag, value) in x.terms)
end

function Base.:*(y::Mathy, x::Differential)
    Differential(tag => y * value for (tag, value) in x.terms)
end

function Base.:*(x::Differential, y::Differential)
    out = Differential([],[0])
    for (k1,v1) in x.terms
        for (k2,v2) in y.terms
            out += Differential([k1*k2], [v1*v2])
        end
    end
    out
end

function Base.:*(t1::DTag,t2::DTag)
    if vcat(t1.tag,t2.tag) |> hasDuplicates
        DTag(-1)
    else
        DTag(vcat(t1.tag,t2.tag) |> sort)
    end
end

#---------------------------------------------------------------
#---------------------------------------------------------------
# constructors for new unary and binary operations

function unaryOp(f::Function, Df::Function)
    function (Dx::Differential)
        f(Dx[1:end-1]) + Df(Dx[1:end-1])*Dx[end]
    end
end

function binaryOp(f::Function, Dfx::Function, Dfy::Function)
    function (Dx::Differential, Dy::Differential)
        f(Dx[1:end-1],Dy[1:end-1]) + Dfx(Dx[1:end-1], Dy[1:end-1])*Dx[end] + Dfy(Dx[1:end-1], Dy[1:end-1])*Dy[end]
    end
end

Base.:-(Dx::Differential, y::Mathy) = unaryOp(x -> x - y, x -> 1)(Dx)
Base.:-(x::Mathy, Dy::Differential) = unaryOp(y -> x - y, y -> -1)(Dy)
Base.:-(Dx::Differential, Dy::Differential) = binaryOp((x,y) -> x - y,
                                                       (x,y) -> 1,
                                                       (x,y) -> -1)(Dx,Dy)

Base.:/(Dx::Differential,y::Mathy) = unaryOp(Dx -> Dx/y, Dx -> 1/y)(Dx)
Base.:/(x::Mathy,Dy::Differential) = unaryOp(Dy -> x/Dy, Dy -> -1/(Dy)^2)(Dy)
Base.:/(Dx::Differential,Dy::Differential) = binaryOp((x,y) -> x/y,
                                                      (x,y) -> 1/y,
                                                      (x,y) -> -x/y^2)(Dx,Dy)

Base.:^(Dx::Differential,y::Mathy) = unaryOp(Dx -> Dx^y, Dx -> y*Dx^(y-1))(Dx)
Base.:^(Dx::Differential,y::Int) = unaryOp(Dx -> Dx^y, Dx -> y*Dx^(y-1))(Dx)
Base.:^(x::Mathy,Dy::Differential) = unaryOp(Dy -> x^Dy, Dy -> log(x)*x^Dy)(Dy)

Base.log(Dx::Differential) = unaryOp(log, x -> 1/x)(Dx)

#---------------------------------------------------------------
#---------------------------------------------------------------
# Derivatives

function makeDiff()
    global tagCount += 1
    Differential([DTag(tagCount)],[1])
end

function extractDiff(Dx::Differential, ϵ::Differential)
    tag = (ϵ |> getTagList)[1]
    out = Differential(tagRemove(t, tag) => v for (t, v) in Dx.terms)
    if getTagList(out) == [DTag()]
        out[1]
    else
        out
    end
end
extractDiff(ex::Operation, ϵ::Differential) = 0

function D(f::Function)
    ϵ = makeDiff()
    x -> extractDiff(f(x + ϵ), ϵ)
end

Operation(*, [1,2])

max_stars_repo_name: MasonProtter/ModellingToolkitDeriv

max_issues_repo_name: MasonProtter/ModellingToolkitDeriv

max_forks_repo_name: MasonProtter/ModellingToolkitDeriv


text:
# SPDX-License-Identifier: MIT

"""
	abstract type ExponentialFamilyDistribution end

Abstract base type for distributions member of the exponential family.

| Methods            | Brief description |
|:-------------------|:------------------|
| `η(p)`             | Returns the natural parameters of the distribution |
| `ξ(p[, η])`        | Parametrization of the distribution as a function of the natural parameters |
| `μ(p)`             | Returns the expectation of the sufficient statistics |
| `A(p[, η])`        | Log-normalizer of the distribution as a function of the natural parameters |
| `kldiv(p, q[, μ])` | Returns the Kullback-Leibler divergence ``D[p || q].`` |
| `sample(p)`        | Draw a sample from `p` |
| `unpack(p, μ)      | Unpack the vector of sufficient statistics |

"""
abstract type ExponentialFamilyDistribution end

μ(p::ExponentialFamilyDistribution) = gradient(η -> A(p, η), η(p))[1]

function kldiv(p₁::ExponentialFamilyDistribution,
               p₂::ExponentialFamilyDistribution,
               μ₁::AbstractVector = μ(p₁))
	η₁, η₂ = η(p₁), η(p₂)
    A(p₂, η₂) - A(p₁, η₁) - dot(η₂ - η₁, μ₁)
end

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
function excess_gibbs_free_energy(model::ActivityModel,p,T,z)
    γ = activity_coefficient(model,p,T,z)
    return sum(z[i]*R̄*T*log(γ[i]) for i ∈ @comps)
end

function eos(model::ActivityModel,V,T,z)
    Σz = sum(z)
    lnΣz = log(Σz)
    g_E = excess_gibbs_free_energy(model,V,T,z)
    pures = model.puremodel
    p = sum(z[i]*pressure(pures[i],V,T) for i ∈ @comps)/Σz
    g_ideal = sum(z[i]*R̄*T*(log(z[i])-lnΣz) for i ∈ @comps)
    g_pure = sum(z[i]*VT_gibbs_free_energy(pures[i],V,T) for i ∈ @comps)
    return g_E+g_ideal+g_pure-p*V
end

function eos_res(model::ActivityModel,V,T,z)
    Σz = sum(z)
    g_E = excess_gibbs_free_energy(model,V,T,z)
    pures = model.puremodel
    g_pure_res = sum(VT_chemical_potential_res(pures[i],V,T)[1]*z[i] for i ∈ @comps)
    p = sum(z[i]*pressure(pures[i],V,T) for i ∈ @comps)/Σz
    p_res = p - Σz*R̄*T/V
    return g_E+g_pure_res-p_res*V
end

function bubble_pressure(model::ActivityModel,T,x)
    sat = saturation_pressure.(model.puremodel,T)
    p_sat = [tup[1] for tup in sat]
    γ     = activity_coefficient(model,1e-4,T,x)
    p     = sum(x.*γ.*p_sat)
    y     = x.*γ.*p_sat ./ p
    return (p,y)
end

function bubble_temperature(model::ActivityModel,p,x;T0=nothing)
    f(z) = Obj_bubble_temperature(model,z,p,x)
    if T0===nothing
        pure = model.puremodel
        sat = saturation_temperature.(pure,p)
        Ti   = zero(x)
        for i ∈ 1:length(x)
            if isnan(sat[i][1])
                Tc,pc,vc = crit_pure(pure[i])
                g(x) = p-pressure(pure[i],vc,x,[1.])
                Ti[i] = Roots.find_zero(g,(Tc))
            else
                Ti[i] = sat[i][1]
            end
        end
        T = Roots.find_zero(f,(minimum(Ti)*0.9,maximum(Ti)*1.1))
    else
        T = Roots.find_zero(f,T0)
    end
    p,y = bubble_pressure(model,T,x)
    return (T,y)
end

function Obj_bubble_temperature(model::ActivityModel,T,p,x)
    sat = saturation_pressure.(model.puremodel,T)
    p_sat = [tup[1] for tup in sat]
    γ     = activity_coefficient(model,1e-4,T,x)
    y     = x.*γ.*p_sat ./ p
    return sum(y)-1
end


text:
using Documenter, DiffEqBayes

include("pages.jl")

makedocs(
    sitename="DiffEqBayes.jl",
    authors="Chris Rackauckas, Vaibhav Kumar Dixit et al.",
    clean=true,
    doctest=false,
    modules=[DiffEqBayes],

    format=Documenter.HTML(assets=["assets/favicon.ico"],
                           canonical="https://diffeqbayes.sciml.ai/stable/"),

    pages=pages
)

deploydocs(
    repo="github.com/SciML/DiffEqBayes.jl.git";
    push_preview=true
)


text:
function Process(sim::Simulation, func::Function, args...)
  task = Task(()->func(sim, args...))
  target = timeout(sim)
  return Process(task, target)
end

function execute(sim::Simulation, ev::Event, proc::Process)
  try
    sim.active_proc = Nullable(proc)
    value = consume(proc.task, ev.value)
    sim.active_proc = Nullable{Process}()
    if istaskdone(proc.task)
      schedule(sim, proc.ev, value=value)
    end
  catch exc
    if !isempty(proc.ev.callbacks)
      schedule(sim, proc.ev, value=exc)
    else
      rethrow(exc)
    end
  end
end

function yield(sim::Simulation, target::Event) :: Any
  if target.state == processed
    return target.value
  end
  proc = get(sim.active_proc)
  proc.target = target
  proc.resume = append_callback(proc.target, execute, proc)
  value = produce(nothing)
  if isa(value, Exception)
    throw(value)
  end
  return value
end

function yield(sim::Simulation, target::Process) :: Any
  yield(sim, target.ev)
end

function interrupt(sim::Simulation, proc::Process, cause::Any=nothing) :: Event
  if !istaskdone(proc.task)
    remove_callback(proc.target, proc.resume)
    proc.target = timeout(sim, priority=true, value=InterruptException(cause))
    proc.resume = append_callback(proc.target, execute, proc)
  end
  return timeout(sim, priority=true)
end


text:
# An MTRBasis for d = 0 or 1 has functions:
#   \{a_{j}(z)\}_{j=1}^{J}
#   \{b_{k}(u)\}_{k=1}^{K}
#   \{\int_{u}^{v} b_{k}(u')\, du\}_{k=1}^{K}
# The basis is then
#   \{a_{j}(z)b_{k}(u) : j = 1,...,J, k = 1...K\}
#
# The rationale for this formulation is:
#   - Allows me to integrate with respect to u fairly easily
#   which would be difficult with say just a single set of functions b_{k}(u,z)
#   - Allows me to keep track of the coefficients by (d,j,k), which is
#   important for imposing certain types of constraints.
@with_kw struct MTRBasis
    a::Vector{Function}
    b::Vector{Function}
    ib::Vector{Function}

    function MTRBasis(a,b,ib)
        @assert length(b) == length(ib)
        new(a, b, ib)
    end
end
export MTRBasis

# Take multiple MTRBasis and concatenate them
function MTRBasis(bases::Vector{MTRBasis})
    anew = vcat([basis.a for basis in bases]...)
    bnew = vcat([basis.b for basis in bases]...)
    ibnew = vcat([basis.ib for basis in bases]...)
    MTRBasis(anew, bnew, ibnew)
end

function bernstein_basis(K)
    a = [(z) -> 1]
    b = [u -> bernstein_polynomial(u, k, K) for k in 0:K]
    ib = [(u, v) -> integrate_bernstein_polynomial(u, v, k, K) for k in 0:K]
    MTRBasis(a, b, ib)
end
export bernstein_basis

function interacted_bernstein_basis(K; notℓ = 1)
    @assert notℓ in [1, 2]

    basis = bernstein_basis(K)
    a = [z -> 1, z -> z[notℓ]]
    b = vcat(basis.b)
    ib = vcat(basis.ib)

    MTRBasis(a, b, ib)
end
export interacted_bernstein_basis

################################################################################
# The kth Bernstein polynomial of degree K, evaluated at u
#   for k = 0,1,...,K
# The formula is:
#   B_{k,K}(u) = nchoosek(K, k) u^(k)*(u^(K-k))
################################################################################
function bernstein_polynomial(u, k::Integer, K::Integer)
    if ((u < 0) | (u > 1))
        throw(DomainError("Should have u in [0,1]."))
    end
    bernstein_errorcheck(k, K)
    return binomial(K,k) * (u.^k) .* ((1 .- u) .^ (K-k));
end
export bernstein_polynomial

function bernstein_errorcheck(k::Integer, K::Integer)
    if ((k < 0) | (k > K))
        throw(DomainError("Should have k in 0,1,...,K."))
    end
    if (K < 0)
        throw(DomainError("Should have K >= 0."))
    end
end

################################################################################
# Integrate B_{k,K}(u) over u
#   where B_{k,K} is a Bernstein polynomial
#
# Input:
#   u: An N-vector of lower integration points
#   v: An N-vector of upper integration points
#   k, K: kth Bernstein polynomial of degree K
#       for k = 0,1,...,K
# Output:
#   An N x 1 column vector with ith component:
#       \int_{u(i)}^{v(i)} b_{k,K}(u)\, du
################################################################################
function integrate_bernstein_polynomial(u::Array{<:Real, 1},
                                        v::Array{<:Real, 1},
                                        k::Integer, K::Integer)
    @assert size(u) == size(v)
    @assert all(u .>= 0) & all(u .<= 1)
    @assert all(v .>= 0) * all(v .<= 1)
    @assert all(u .<= v)
    bernstein_errorcheck(k, K)

    BI = zeros(size(u));
    for i = k:1:K
        BI =  BI +
              binomial(K, i) *
              binomial(i, k) *
              (-1)^(i - k) .*
              (v.^(i+1) - u.^(i+1)) ./ (i+1);
    end
    return BI
end

function integrate_bernstein_polynomial(u::Real, v::Real,
                                        k::Integer, K::Integer)
    integrate_bernstein_polynomial([u], [v], k, K)[1]
end
export integrate_bernstein_polynomial

function constantspline_basis(knots::Vector{<:Real})
    a = [(z) -> 1]
    @assert all(0 .<= knots .<= 1)
    @assert (0 in knots) & (1 in knots)
    unique!(knots)
    sort!(knots)
    b = [u -> Int.(knots[k-1] .<= u .<= knots[k]) for k in 2:length(knots)]
    ib = [(u, v) -> max(0, (min(v, knots[k]) - max(u, knots[k-1])))
          for k in 2:length(knots)]
    MTRBasis(a, b, ib)
end
export constantspline_basis

function interacted_constantspline_basis(knots::Vector{<:Real}; notℓ = 1)
    @assert notℓ in [1, 2]

    basis = constantspline_basis(knots)
    a = [z -> 1, z -> z[notℓ]]
    b = vcat(basis.b)
    ib = vcat(basis.ib)

    MTRBasis(a, b, ib)
end
export interacted_constantspline_basis

function switching_constantspline_basis(knots::Vector{<:Real}; notℓ = 1)
    @assert notℓ in [1, 2]

    basis = interacted_constantspline_basis(knots)
    a = [z -> (1 - z[notℓ]), z -> z[notℓ]]
    MTRBasis(a, basis.b, basis.ib)
end
export switching_constantspline_basis


text:
export setexperiments!, setexperiment!, setsampledata!, setsamplemap!, setmetadata!
import DataStructures
import DataFrames
import SummarizedExperiments

"""
    setexperiments!(x, value)

Set the experiments in the `MultiAssayExperiment` `x` to the `OrderedDict` `value`.
This returns a reference to the modified `x`.

# Examples
```jldoctest
julia> using MultiAssayExperiments

julia> x = exampleobject();

julia> y = copy(experiments(x));

julia> delete!(y, "foo");

julia> setexperiments!(x, y);

julia> collect(keys(experiments(x)))
1-element Vector{String}:
 "bar"
```
"""
function setexperiments!(x::MultiAssayExperiment, value::DataStructures.OrderedDict{String,SummarizedExperiments.SummarizedExperiment})
    x.experiments = value;
    return x
end

"""
    setexperiment!(x[, i], value)

Set experiment `i` in `MultiAssayExperiment` `x` to the `SummarizedExperiment` `value`.
This returns a reference to the modified `x`.

`i` may be a positive integer, in which case it should be no greater than the length of `experiments(x)`.
It may also be a string specifying a new or existing experiment in `x`.
If omitted, we set the first experiment by default.

# Examples
```jldoctest
julia> using MultiAssayExperiments;

julia> x = exampleobject();

julia> size(experiment(x, 2))
(50, 8)

julia> val = experiment(x);

julia> setexperiment!(x, 2, val);

julia> size(experiment(x, 2))
(100, 10)
```
"""
function setexperiment!(x::MultiAssayExperiment, value::SummarizedExperiments.SummarizedExperiment)
    return setexperiment!(x, 1, value)
end

function setexperiment!(x::MultiAssayExperiment, i::Int, value::SummarizedExperiments.SummarizedExperiment)
    counter = 1
    for key in keys(experiments(x))
        if i == counter
            x.experiments[key] = value
            return x
        end
        counter += 1
    end
    throw(BoundsError("experiment " * string(i) * " is out of range (" * length(experiments(x)) * " experiments available)"))
end

function setexperiment!(x::MultiAssayExperiment, i::String, value::SummarizedExperiments.SummarizedExperiment)
    x.experiments[i] = value
    return x
end

"""
    setsampledata!(x, value)

Set the sample data in the `MultiAssayExperiment` `x` to the `DataFrame` `value`.

The returned object should contain `name` as the first column, containing a vector of unique strings.
If `check = true`, the function will check the validity of the sample data before returning it.

# Examples
```jldoctest
julia> using MultiAssayExperiments

julia> x = exampleobject();

julia> sd = copy(sampledata(x));

julia> sd[!,"stuff"] = [rand() for i in 1:size(sd)[1]];

julia> setsampledata!(x, sd);

julia> names(sampledata(x))
3-element Vector{String}:
 "name"
 "disease"
 "stuff"
```
"""
function setsampledata!(x::MultiAssayExperiment, value::DataFrames.DataFrame)
    check_sampledata(value)
    x.sampledata = value
    return x
end

"""
    setsamplemap!(x, value)

Set the sample mapping in the `MultiAssayExperiment` `x` to a `DataFrame` `value`.
This returns a reference to the modified `x`.

`value` should contain the `sample`, `experiment` and `colname` columns in that order.
Each column should contain a vector of strings:

- Values of `sample` may (but are not required to) correspond to the names of samples in `sampledata(x)`.
- Values of `experiment` may (but are not required to) correspond to the keys of `experiments(x)`.
- Values of `colname` should (but are not required to) correspond to the columns of the corresponding `SummarizedExperiment` in the `experiment` of the same row.

This correspondence is used for convenient subsetting and extraction, e.g., [`expandsampledata`](@ref), [`filtersamplemap`](@ref).
However, values in the sample mapping columns need not have a 1:1 match to their corresponding target;
any values unique to one or the other will be ignored in the relevant methods.
This allows users to flexibly manipulate the object without constantly hitting validity checks.

It is legal (but highly unusual) for a given combination of `experiment` and `colname` to occur more than once.
This may incur warnings in methods like [`expandsampledata`](@ref).

# Examples
```jldoctest
julia> using MultiAssayExperiments

julia> x = exampleobject();

julia> y = samplemap(x)[1:10,:];

julia> setsamplemap!(x, y);

julia> size(samplemap(x))[1]
10
```
"""
function setsamplemap!(x::MultiAssayExperiment, value::DataFrames.DataFrame)
    check_samplemap(value)
    x.samplemap = value
    return x
end

"""
    setmetadata!(x, value)

Set the metadata of a `MultiAssayExperiment` `x` to a dictionary `value`.
This returns a reference to the modified `x`.

# Examples
```jldoctest
julia> using MultiAssayExperiments

julia> x = exampleobject();

julia> meta = copy(metadata(x));

julia> meta["version"] = "0.2.0";

julia> setmetadata!(x, meta);

julia> metadata(x)["version"]
"0.2.0"
```
"""
function setmetadata!(x::MultiAssayExperiment, value::Dict{String,Any})
    x.metadata = value
    return x
end


text:
# Contains all functions handling geometry

"""
    Calculate spherical coordinates arrays representing the surface of the particle
This can be either 1D arrays for axi-symmetric particles or 2D arrays for general 3D particles
"""
function get_r_θ_ϕ_arrays(geometry; angular_resolution=0.5, θ_array=nothing, ϕ_array=nothing)
    # TODO: use a package that takes a mesh and do integrals
    # check: https://github.com/jareeger/Smooth_Closed_Surface_Quadrature_RBF-julia

    if geometry isa String
        if geometry in ["wire", "cylinder"]
        elseif geometry in ["spheroid", "ellipsoid"]
        end
    elseif geometry isa Function
        # TODO: check if the function has one or two arguments

    else
        throw(DomainError("Please define geometry as 'String', 'Function', or 'Filepath'"))
    end
end

function create_array_from_elements(x1, x2, x3)
    # TODO: @Alok, should I use SVector? what are the drawbacks?
    # return [x1,x2,x3]
    return SVector(x1, x2, x3)
end

"""
    Claculate radial coordinate of ellipsoidal particle surface, given θ points
returns r and n̂ coordinate as a function of θ_array
"""
function ellipsoid(rx, rz, θ_array)
    r = rx .* rz ./ sqrt.((rx .* cos.(θ_array)).^2 + (rz .* sin.(θ_array)).^2)
    ∂r_by_∂θ = (rx.^2 - rz.^2) / (rx.^2 * rz.^2) .* r.^3 .* sin.(θ_array) .* cos.(θ_array)

    n̂_r_comp = r ./ sqrt.(r.^2 + ∂r_by_∂θ.^2)
    n̂_θ_comp = -∂r_by_∂θ ./ sqrt.(r.^2 + ∂r_by_∂θ.^2)
    n̂_ϕ_comp = zeros(Real, size(r)) # TODO: I need to make the type the same as the type of r

    # n̂ = create_array_from_elements.(n̂_r_comp, n̂_θ_comp, n̂_ϕ_comp) # TODO: verify that you can delete this line
    # n̂ = SVector.(n̂_r_comp, n̂_θ_comp, n̂_ϕ_comp)
    # n̂ = Vector.(n̂_r_comp, n̂_θ_comp, n̂_ϕ_comp)
    n̂ = reshape([Vector([n̂_r_comp[id], n̂_θ_comp[id], n̂_ϕ_comp[id]]) for id in eachindex(n̂_r_comp)], size(n̂_r_comp))
    return r, n̂
end

"""
    Create 1D mesh, then get r, θ and n̂_θ_comp for points between each two consecutive points defined by r_edges_array, θ_edges_array
It returns a Matrix, first and second columns are for r,θ, respectively. The third, fourth and fifth columns are for the three compoonents of the normal direction.
"""
function create_1D_mesh_and_get_r_n̂(r_edges_array::AbstractVector{R}, θ_edges_array::AbstractVector{R}) where R <: Real
    z = find_point_splitting_line_between_two_points_and_normal_vector_to_this_line.(r_edges_array[1:end - 1], θ_edges_array[1:end - 1], r_edges_array[2:end], θ_edges_array[2:end]);
    """
    r_array = [z[id][1] for id = 1 : length(z)]
    θ_array = [z[id][2] for id = 1 : length(z)]
    n̂_array = [z[id][3] for id = 1 : length(z)]
    """
    return vcat(z...)
end

function create_1D_mesh_and_get_r_n̂_with_plotting(r_edges_array::AbstractVector{R}, θ_edges_array::AbstractVector{R}) where R <: Real
    z = create_1D_mesh_and_get_r_n̂(r_edges_array, θ_edges_array);
end

"""
    Find a point (r,θ) that split a line connecting two points (r1, θ1), and (r2, θ2) into two equal splits, and find the normal unit vector is spherical coordinates
"""
function find_point_splitting_line_between_two_points_and_normal_vector_to_this_line(r1::R, θ1::R, r2::R, θ2::R) where R <: Real
    S = calculate_distance_between_two_points(r1, θ1, r2, θ2)
    angle_between_r1_and_S = asin(r2 * sin(θ2 - θ1) / S)
    # to fix the inverse sin of obtuse angle
    if r2 > r1; angle_between_r1_and_S = pi - angle_between_r1_and_S; end
    r = sqrt(r1^2 + (S / 2)^2 - 2 * r1 * (S / 2) * cos(angle_between_r1_and_S))
    angle_between_r1_and_r = asin(S / 2 * sin(angle_between_r1_and_S) / r) # TODO: do the check for the "obtuse" angle
    θ = angle_between_r1_and_r + θ1

    γ = pi - angle_between_r1_and_S - angle_between_r1_and_r

    n̂_θ_comp = -cos(γ)
    n̂_r_comp = sin(γ)

    # println("r1=$r1, θ1=$θ1, r2=$r2, θ2=$θ2, r=$r, S=$S, θ=$θ, angle_between_r1_and_S/pi=$(angle_between_r1_and_S/pi), asin($(r2 * sin(θ2-θ1) / S)), θ2-θ1=$(θ2-θ1)")
    return hcat(r, θ, n̂_r_comp, n̂_θ_comp, zero(R))
    # return hcat(r, θ, n̂_θ_comp)
end

"""
    Find distance between two points (r1, θ1), and (r2, θ2)
"""
function calculate_distance_between_two_points(r1::R, θ1::R, r2::R, θ2::R) where R <: Real
    return sqrt(r1^2 + r2^2 - 2 * r1 * r2 * cos(θ2 - θ1))
end

"""
    Calculate the total length of a 1D mesh, defined as r and θ arrays.
"""
function calculate_total_length_of_1D_mesh(r_array::AbstractVector{R}, θ_array::AbstractVector{R}) where R <: Real
    return sum(
        calculate_distance_between_two_points.(
            r_array[1:end - 1], θ_array[1:end - 1], r_array[2:end], θ_array[2:end]
        )
    )
end

"""
    Calculate the surface area of an axisymmetric particle, defined by a 1D mesh.
"""
function calculate_surface_area_of_axisymmetric_particle(r_array::AbstractVector{R}, θ_array::AbstractVector{R}) where R <: Real
    return sum(
        calculate_lateral_surface_area_of_truncated_cone.(
            r_array[1:end - 1], θ_array[1:end - 1], r_array[2:end], θ_array[2:end]
        )
    )
end

"""
    Calculate surface area of a truncated cone defined by two points revolving about the Z axis.
"""
function calculate_lateral_surface_area_of_truncated_cone(r1::R, θ1::R, r2::R, θ2::R) where R <: Real
    r = r1 * sin(θ1)
    RR = r2 * sin(θ2)
    s = abs(r1 * cos(θ1) - r2 * cos(θ2))
    return pi * (r + RR) * sqrt((r - RR)^2 + s^2)
end

"""
    Calculate the volume of an axisymmetric particle, defined by a 1D mesh.
"""
function calculate_volume_of_axisymmetric_particle(r_array::AbstractVector{R}, θ_array::AbstractVector{R}) where R <: Real
    return sum(
        calculate_volume_of_truncated_cone.(
            r_array[1:end - 1], θ_array[1:end - 1], r_array[2:end], θ_array[2:end]
        )
    )
end

"""
    Calculate volume of a truncated cone defined by two points revolving about the Z axis.
"""
function calculate_volume_of_truncated_cone(r1::R, θ1::R, r2::R, θ2::R) where R <: Real
    r = r1 * sin(θ1)
    RR = r2 * sin(θ2)
    s = abs(r1 * cos(θ1) - r2 * cos(θ2))
    return 1 / 3 * pi * (r^2 + r * RR + RR^2) * s
end

"""
    Calculate surface area for an ellipsoid
    following equations in SMARTIES paper: 1. Somerville, W. R. C., Auguié, B. & Le Ru, E. C. Smarties: User-friendly codes for fast and accurate calculations of light scattering by spheroids. J. Quant. Spectrosc. Radiat. Transf. 174, 39–55 (2016).
    verified against results from this area and volume calculator https://planetcalc.com/149/
"""
function get_volume_area_for_ellipsoid(rx, rz)
    h = rx / rz
    if h < 1; h = 1 / h; end  # h is the aspect ratio that is larger than 1
    e = sqrt(h^2 - 1) / h

    volume = 4 / 3 * pi * rx^2 * rz  # eq. 9
    if rx > rz  # oblate
        area_surface = (
            2 * pi * rx^2 * (1 + ((1 - e^2) / e) * atanh(e))
        )  # eq. 11
    elseif rx == rz  # sphere
        area_surface = 4 * pi * rx^2
    else  # prolate
        area_surface = 2 * pi * rx^2 * (1 + rz / (rx * e) * asin(e))  # eq. 11
    end
    return volume, area_surface
end

function get_volume_area_for_geometry(geometry_name, geometry_parameters...)
    if geometry_name in ["spheroid", "ellipsoid"]
        return get_volume_area_for_ellipsoid(geometry_parameters...)
    end
end

function convert_coordinates_Cart2Sph(x, y, z)
    r = sqrt((x^2) + (y^2) + (z^2))
    theta = atan(sqrt((x^2) + (y^2)), z)
    phi = atan(y, x)
    return hcat(r, theta, phi)
end

function convert_coordinates_Sph2Cart(r, theta, phi)

    x = r * sin(theta) * cos(phi)
    y = r * sin(theta) * sin(phi)
    z = r * cos(theta)

    return hcat(x, y, z)
end

"""
    Given r and θ arrays, this function returns another r and θ arrays covering the same range of θ, with normal direction vector
Note that the first and last elements of θ_array are assumed to be 0 and pi, respectively (or numbers very close to these). We assume that the n̂ is pointing at +ve and -ve z directions at these points.
"""
function get_r_θ_n̂_arrays_from_r_θ_arrays_axisymmetrix(r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}) where R <: Real
    # TODO: I need to make sure that the first and last elements of θ_array are 0 and pi, respectively (or numbers very close to these). Note that we assume that the n̂ is pointing at +ve and -ve z directions at these points.
    r_theta_n̂ = create_1D_mesh_and_get_r_n̂(r_array, θ_array)
    return vcat(
        [r_array[1], θ_array[1], one(R), zero(R), zero(R)]',
        r_theta_n̂,
        [r_array[end], θ_array[end], one(R), zero(R), zero(R)]',
    )
end

"""
Increase mesh density by 4x, by dividing every line element into two equal parts.
Works only for 1D mesh.
"""
function double_mesh_density(r_array, θ_array)
    z = Tmatrix.create_1D_mesh_and_get_r_n̂(r_array, θ_array)
    r_middles = z[:,1]
    θ_middles = z[:,2]
    # return hcat(r_middles, θ_middles)

    """
    r_all = vcat(r_array,r_middles)
    r_all[1:2:end] = r_array
    r_all[2:2:end] = r_middles

    θ_all = vcat(θ_array,θ_middles)
    θ_all[1:2:end] = θ_array
    θ_all[2:2:end] = θ_middles
    return hcat(r_all, θ_all)
    """

    return hcat(
        add_elements_from_two_arrays_in_turn(r_array, r_middles),
        add_elements_from_two_arrays_in_turn(θ_array, θ_middles),
        )
end

"""
Increase mesh density by 4x, by dividing every line element into four equal parts.
Works only for 1D mesh.
"""
function quadruple_mesh_density(r_array, θ_array)
    r_theta_array_double = double_mesh_density(r_array, θ_array)
    return double_mesh_density(r_theta_array_double[:,1], r_theta_array_double[:,2])
end

"""
Increase mesh density by 8x, by dividing every line element into eight equal parts.
Works only for 1D mesh.
"""
function octuple_mesh_density(r_array, θ_array)
    r_theta_array_quadruple = quadruple_mesh_density(r_array, θ_array)
    return double_mesh_density(r_theta_array_quadruple[:,1], r_theta_array_quadruple[:,2])
end

"""
Increase mesh density by 16x, by dividing every line element into sixteen equal parts.
Works only for 1D mesh.
"""
function sexdecuple_mesh_density(r_array, θ_array)
    r_theta_array_octuple = octuple_mesh_density(r_array, θ_array)
    return double_mesh_density(r_theta_array_octuple[:,1], r_theta_array_octuple[:,2])
end

"""
julia> a = [1,2,3]; b = [10,20];
julia> add_elements_from_two_arrays_in_turn(a,b)
julia> 5-element Vector{Int64}:
        1
        10
        2
        20
        3
"""
function add_elements_from_two_arrays_in_turn(array_1, array_2)
    return [isodd(i) ? array_1[Int((i + 1) / 2)] : array_2[Int(i / 2)] for i = 1:(length(array_1) + length(array_2))]
end


text:
using GLM, DataFrames, Gadfly, Distributions
include("misc.jl")

function avPlot(dfrm::RegressionModel, variable::Symbol)
    varNames = names(dfrm.mf.df)
    # get e_yx
    lhs = dfrm.mf.terms.eterms[1]
    rhs = extractrhs(dfrm)
    rhsnew = replace(rhs.rhsstring, "+"*string(variable), "")
    rhsnew = rhsnew[2:end]
    rhsnew = parse(rhsnew)
    fnew = Formula(lhs, rhsnew)
    newfit = fit(LinearModel, fnew, dfrm.mf.df)
    e_yx = residuals(newfit)
    # get e_zx
    X = dfrm.mf.df[2:end]
    Z = X[variable]
    X = X[~[(name in [variable]) for name in names(X)]]
    intercept = ones(Int, nrow(X), 1)
    X = convert(Array, X)
    X = hcat(intercept, X)
    H = *(*(X, inv(*(transpose(X), X))), transpose(X))
    e_zx = *((eye(size(H)[1])-H), Z)
    df = DataFrame(e_yx=e_yx, e_zx=e_zx)
    # get alpha
    avlm = fit(LinearModel, @formula(e_yx~0+e_zx), df)
    alpha = coef(avlm)
    # plot
    xs_range = abs(maximum(df[:e_zx]) - minimum(df[:e_zx]))
    xs = linspace(minimum(df[:e_zx])-xs_range/10, maximum(df[:e_zx])+xs_range/10)
    ys = alpha .* xs
    xy = DataFrame(xs=xs, ys=ys)
    plot(layer(df, x="e_zx", y="e_yx", Geom.point), Guide.xlabel(string(variable)*" | others"),
         Guide.ylabel(string(varNames[1])*" | others"),
         Guide.title("Added-Variable Plot"),
         layer(xy, x="xs", y="ys", Theme(default_color=color("red")), Geom.line))
end

function vif(dfrm::RegressionModel)
    X = dfrm.mf.df[2:end]
    rhs = extractrhs(dfrm)
    result = DataFrame(variable=rhs.rhsarray, vif=0.0)
    i = 1
    for var in rhs.rhsarray
        lhs = parse(var)
        rhsnew = replace(rhs.rhsstring, "+"*var, "")
        rhsnew = rhsnew[2:end]
        rhsnew = parse(rhsnew)
        fnew = Formula(lhs, rhsnew)
        newfit = fit(LinearModel, fnew, X)
        r2 = rsquared(newfit)
        result[:vif][i] = 1 / (1-r2)
        i = i + 1
    end
    result
end

function rsquared(dfrm::RegressionModel)
    SStot = sum((dfrm.model.rr.y - mean(dfrm.model.rr.y)).^2)
    SSres = sum((dfrm.model.rr.y - dfrm.model.rr.mu).^2)
    return (1-(SSres/SStot))
end

function adjrsquared(dfrm::RegressionModel)
    SStot = sum((dfrm.model.rr.y - mean(dfrm.model.rr.y)).^2)
    SSres = sum((dfrm.model.rr.y - dfrm.model.rr.mu).^2)
    n = size(dfrm.model.rr.y, 1)  #number of samples
    if dfrm.mf.terms.intercept
        p = size(dfrm.mm.m, 2) - 1
    else
        p = size(dfrm.mm.m, 2)
    end
    return 1- ( (SSres/(n-p-1)) / (SStot/(n-1)) )
end

function rstudent(dfrm::RegressionModel)
    SSres = sum((dfrm.model.rr.y - dfrm.model.rr.mu).^2)
    n = size(dfrm.model.rr.y, 1)  #number of samples
    if dfrm.mf.terms.intercept
        p = size(dfrm.mm.m, 2) - 1
    else
        p = size(dfrm.mm.m, 2)
    end
    sigma2 = SSres / (n-p)
    X = dfrm.mm.m
    H = *(*(X, inv(*(transpose(X), X))), transpose(X))
    h = diag(H)
    r = residuals(dfrm) ./ (sqrt(sigma2 .* (1 - h)))
    return r
end

function jackknife(dfrm::RegressionModel)
    r = rstudent(dfrm)
    n = size(dfrm.model.rr.y, 1)  #number of samples
    if dfrm.mf.terms.intercept
        p = size(dfrm.mm.m, 2) - 1
    else
        p = size(dfrm.mm.m, 2)
    end
    t = r .* sqrt((n-p-1) ./ (n-p-r.^2))
    return t
end

function halfnorm(dfrm::RegressionModel)
    N = size(dfrm.model.rr.y, 1)
    n = 1:N
    X = dfrm.mm.m
    H = *(*(X, inv(*(transpose(X), X))), transpose(X))
    h = diag(H)
    labels = sortperm(h)
    h = h[labels]
    U = (N+n) / (2*N+1)
    d = Normal()
    u = quantile(d, U)
    # prepare labels for potential outliers
    df = DataFrame(u=u, h=h, label="")
    outlierlabel = labels[u.>2]
    no_outliers = length(outlierlabel)
    for i in (length(labels)-no_outliers+1):length(labels)
        df[:label][i] = string(labels[i])
    end
    plot(df, x="u", y="h", label="label", Geom.point, Geom.label, Guide.xlabel("Half-normal quantiles"), Guide.ylabel("Leverages"))
end

function cooksdistance(dfrm::RegressionModel; plotit::Bool=false)
    r = rstudent(dfrm)
    X = dfrm.mm.m
    H = *(*(X, inv(*(transpose(X), X))), transpose(X))
    h = diag(H)
    p = size(dfrm.mm.m, 2)
    d = (r.^2 .* h) ./ (p*(1-h))
    if !plotit
        return d
    else
        # identify outliers
        n = length(d)
        cutoff = 4/(n-p)
        labels = 1:n
        outlierlabel = labels[d.>cutoff]
        df = DataFrame(n=labels, d=d, label="")
        for i in outlierlabel
            df[:label][i] = string(i)
        end
        plot(df, x="n", y="d", label="label", Guide.xlabel("Case number"), Guide.ylabel("Cook's distance"), Geom.bar, Geom.label)
    end
end

max_stars_repo_name: atantos/RegTools

max_issues_repo_name: atantos/RegTools

max_forks_repo_name: atantos/RegTools


text:
"""
Describes a single parameter to a model function.
"""
@compat abstract type ModelParameter end
@compat const ModelParameters = Dict{Symbol, ModelParameter}

"""
Describes a free (non-constant) parameter to a model function.
"""
type FreeParameter <: ModelParameter

    "Index of the parameter in the function's argument list."
    position::Integer

    "The initial guess used for fitting."
    value::AbstractFloat

    "Best-fit value of the parameter, if it exists."
    fit_value::Nullable{AbstractFloat}

    "Uncertainty on the best-fit value of the parameter, if it exists."
    fit_uncertainty::Nullable{AbstractFloat}

    """
        FreeParameter(position, guess)

    Creates a new FreeParameter object with indexed at `position` with
    initial guess `guess`.
    """
    function FreeParameter(position, guess)
        instance = new()

        instance.position = position
        instance.value    = guess

        instance.fit_value       = Nullable{AbstractFloat}()
        instance.fit_uncertainty = Nullable{AbstractFloat}()

        instance
    end

end

"""
Describes a fixed (constant) parameter to a model function.
"""
type FixedParameter <: ModelParameter

    "Index of the parameter in the function's argument list."
    position::Integer

    "The value to which the parameter is fixed."
    value::AbstractFloat

    """
        FixedParameter(position, value)

    Creates a new FixedParameter object with indexed at `position` with
    constant value `value`.
    """
    function FixedParameter(position, value)
        instance = new()

        instance.position = position
        instance.value    = value

        instance
    end

end


text:
using QuarkDroplets
using Test

@testset "QuarkDroplets.jl" begin

    R = π
    filename = string("results/", R, ".jld")
    if isfile(filename)
        @info "Simulation already ran"
    else
        xroot, lroot, ϵroot, κ, norm, coef, energ, j, mass = njlradius(R)
        saveradius(R, xroot, lroot, ϵroot, κ, norm, coef, energ, j, mass)
    end
    rm(filename)

end

max_stars_repo_name: avramescudana/QuarkDroplets

max_issues_repo_name: avramescudana/QuarkDroplets

max_forks_repo_name: avramescudana/QuarkDroplets


text:
using Random, LinearAlgebra, BlockArrays
using BlockArrays: _BlockArray

function dense_BlockMatrix_BlockVector_mul_tests(rng, X, y)
    Ps = blocksizes(X, 1)
    z = _BlockArray([randn(rng, P) for P in Ps], Ps)

    @test mul!(z, X, y) isa AbstractBlockVector
    @test Vector(mul!(z, X, y)) ≈ Matrix(X) * Vector(y)

    @test X * y isa AbstractBlockVector
    @test Vector(X * y) ≈ Matrix(X) * Vector(y)

    z̄ = _BlockArray([randn(rng, P) for P in Ps], Ps)
    z, back = Zygote.pullback(*, X, y)
    X̄, ȳ = back(z̄)
    @test X̄ isa AbstractBlockMatrix
    @test ȳ isa AbstractBlockVector

    z_dense, back_dense = Zygote.pullback(*, Matrix(X), Vector(y))
    X̄_dense, ȳ_dense = back_dense(Vector(z̄))
    @test Vector(z) ≈ z_dense
    @test Matrix(X̄) ≈ X̄_dense
    @test Vector(ȳ) ≈ ȳ_dense

    X̄, ȳ = back(Vector(z̄))
    @test Matrix(X̄) ≈ X̄_dense
    @test Vector(ȳ) ≈ ȳ_dense
end

function dense_BlockMatrix_BlockMatrix_mul_tests(rng, X, Y)
    Ps, Qs = blocksizes(X, 1), blocksizes(Y, 2)
    Z = _BlockArray([randn(rng, P, Q) for P in Ps, Q in Qs], Ps, Qs)

    @test mul!(Z, X, Y) isa AbstractBlockMatrix
    @test Matrix(mul!(Z, X, Y)) ≈ Matrix(X) * Matrix(Y)

    @test X * Y isa AbstractBlockMatrix
    @test Matrix(X * Y) ≈ Matrix(X) * Matrix(Y)

    Z̄ = _BlockArray([randn(rng, P, Q) for P in Ps, Q in Qs], Ps, Qs)
    Z, back = Zygote.pullback(*, X, Y)
    X̄, Ȳ = back(Z̄)
    @test X̄ isa AbstractBlockMatrix
    @test Ȳ isa AbstractBlockMatrix

    Z_dense, back_dense = Zygote.pullback(*, Matrix(X), Matrix(Y))
    X̄_dense, Ȳ_dense = back_dense(Matrix(Z̄))
    @test Matrix(Z) ≈ Z_dense
    @test Matrix(X̄) ≈ X̄_dense
    @test Matrix(Ȳ) ≈ Ȳ_dense

    X̄, Ȳ = back(Matrix(Z̄))
    @test Matrix(X̄) ≈ X̄_dense
    @test Matrix(Ȳ) ≈ Ȳ_dense
end

@testset "fdm stuff" begin
    rng, Ps, Qs = MersenneTwister(123456), [5, 4], [3, 2, 1]
    X = _BlockArray([randn(rng, P, Q) for P in Ps, Q in Qs], Ps, Qs)
    vec_X, from_vec = FiniteDifferences.to_vec(X)
    @test vec_X isa Vector
    @test from_vec(vec_X) == X
end


text:
"""
    trajectory(trajs)

Convert an array of datum to a concatenated NamedTuple data.

# Examples
```jldoctest
julia> using LazyFym

julia> trajectory([(; t=1, x=[1, 2, 3]), (; t=2, x=[2, 3, 4])])
(t = [1, 2], x = [[1, 2, 3], [2, 3, 4]])
```
"""
function trajectory(trajs)
    _trajs = trajs |> collect
    all_keys = union([keys(traj) for traj in _trajs]...)
    get_values(key) = [get(traj, key, missing) for traj in _trajs]
    all_values = all_keys |> Map(get_values)
    return (; zip(all_keys, all_values)...)  # NamedTuple
end

"""
    catTrajectory(multiple_trajs)

Concatenate a collection of evaluated trajectories.
It is probably useful to concatenate the result of parallel simulation.
Applying `catTrajectory` to a single trajectory or
`evaluate` to multiple trajectories would yield undesirable results.

# Examples
```jldoctest
julia> using LazyFym

julia> catTrajectory([(t = [1, 2], x = [[1, 2, 3], [2, 3, 4]]), (t = [3, 4], x = [[3, 4, 5], [4, 5, 6]])])
(t = [1, 2, 3, 4], x = [[1, 2, 3], [2, 3, 4], [3, 4, 5], [4, 5, 6]])

julia> catTrajectory([(; t=1, x=[1, 2, 3]), (; t=2, x=[2, 3, 4])])
(t = [1, 2], x = [1, 2, 3, 2, 3, 4])

julia> evaluate([(t = [1, 2], x = [[1, 2, 3], [2, 3, 4]]), (t = [3, 4], x = [[3, 4, 5], [4, 5, 6]])])
(t = [[1, 2], [3, 4]], x = [[[1, 2, 3], [2, 3, 4]], [[3, 4, 5], [4, 5, 6]]])
```
"""
function catTrajectory(multiple_trajs)
    _multiple_trajs = multiple_trajs |> collect
    all_keys = union([keys(trajs) for trajs in _multiple_trajs]...)
    get_values(key) = vcat([get(trajs, key, missing) for trajs in _multiple_trajs]...)
    all_values = all_keys |> Map(get_values)
    return (; zip(all_keys, all_values)...)  # NamedTuple
end


text:
using SPIRV
using Graphs
using Test

resource(filename) = joinpath(@__DIR__, "resources", filename)

modules = [
    "vert.spv",
    "frag.spv",
    "comp.spv",
    "decorations.spv",
]

@testset "SPIRV.jl" begin
    @testset "Testing SPIR-V module $file" for file in modules
        r = resource(file)
        pmod = PhysicalModule(r)
        mod = SPIRV.Module(pmod)

        @testset "Assembly/disassembly isomorphisms" begin
            pmod_reconstructed = PhysicalModule(mod)
            @test pmod == pmod_reconstructed

            @test sizeof(assemble(pmod)) == stat(r).size
            @test assemble(pmod) == assemble(mod)

            tmp = IOBuffer()
            write(tmp, pmod_reconstructed)
            seekstart(tmp)
            @test PhysicalModule(tmp) == pmod_reconstructed
        end
    end

    @testset "Intermediate Representation" begin
        mod = SPIRV.Module(resource("vert.spv"))
        ir = IR(mod)
        f = ir.fdefs[4]
        @test length(f.cfg.blocks) == nv(f.cfg.graph) == count(==(SPIRV.OpLabel), map(x -> x.opcode, SPIRV.body(f)))
        mod2 = SPIRV.Module(ir)
        pmod2 = PhysicalModule(mod2)
        @test SPIRV.Module(pmod2) == mod2

        mod = SPIRV.Module(resource("comp.spv"))
        ir = IR(mod)
        f = ir.fdefs[52]
        @test length(f.cfg.blocks) == nv(f.cfg.graph) == count(==(SPIRV.OpLabel), map(x -> x.opcode, SPIRV.body(f)))
        mod2 = SPIRV.Module(ir)
        pmod2 = PhysicalModule(mod2)
        @test SPIRV.Module(pmod2) == mod2
    end
end


text:
"Parameters for right-hand side function of fully implicit Runge-Kutta methods."
type NonlinearFunctionParametersFIRK{T,FT} <: NonlinearFunctionParameters{T}
    f::FT
    Δt::T

    d::Int
    s::Int

    a::Matrix{T}
    c::Vector{T}

    t::T
    x::Vector{T}
    y::Vector{T}

    X::Matrix{T}
    Y::Matrix{T}
    F::Matrix{T}

    tX::Vector{T}
    tF::Vector{T}

    function NonlinearFunctionParametersFIRK(f, Δt, d, s, a, c)
        # create solution vectors
        x = zeros(T,d)
        y = zeros(T,d)

        # create internal stage vectors
        X = zeros(T,d,s)
        Y = zeros(T,d,s)
        F = zeros(T,d,s)

        # create temporary vectors
        tX = zeros(T,d)
        tF = zeros(T,d)

        new(f, Δt, d, s, a, c, 0, x, y, X, Y, F, tX, tF)
    end
end

"Compute stages of fully implicit Runge-Kutta methods."
function function_stages!{T,FT}(y::Vector{T}, b::Vector{T}, params::NonlinearFunctionParametersFIRK{T,FT})
    local tᵢ::T

    for i in 1:params.s
        tᵢ = params.t + params.Δt * params.c[i]

        for k in 1:params.d
            # copy y to Y
            params.Y[k,i] = y[params.d*(i-1)+k]

            # compute X
            params.X[k,i] = params.x[k] + params.Δt * params.Y[k,i]
        end

        # compute f(X)
        simd_copy_xy_first!(params.tX, params.X, i)
        params.f(tᵢ, params.tX, params.tF)
        simd_copy_yx_first!(params.tF, params.F, i)
    end

    # compute b = - (Y-AF)
    for i in 1:params.s
        for k in 1:params.d
            b[params.d*(i-1)+k] = - params.Y[k,i]
            for j in 1:params.s
                b[params.d*(i-1)+k] += params.a[i,j] * params.F[k,j]
            end
        end
    end
end

"Fully implicit Runge-Kutta integrator."
immutable IntegratorFIRK{T, FT, ST, IT} <: Integrator{T}
    equation::ODE{T,FT}
    tableau::TableauFIRK{T}
    Δt::T

    solver::ST
    iguess::InitialGuess{T, FT, IT}

    x::Array{T,1}
    y::Array{T,1}
    X::Array{T,2}
    Y::Array{T,2}
    F::Array{T,2}
end

function IntegratorFIRK{T,FT}(equation::ODE{T,FT}, tableau::TableauFIRK{T}, Δt::T;
                              nonlinear_solver=QuasiNewtonSolver,
                              interpolation=HermiteInterpolation{T})
    D = equation.d
    S = tableau.s

    # create solution vector for internal stages / nonlinear solver
    z = zeros(T, D*S)

    # create params
    params = NonlinearFunctionParametersFIRK{T,FT}(equation.f, Δt, D, S, tableau.a, tableau.c)

    # create solver
    solver = nonlinear_solver(z, params)

    # create initial guess
    iguess = InitialGuess(interpolation, equation, Δt)

    # create integrator
    IntegratorFIRK{T, FT, typeof(solver), typeof(iguess.int)}(equation, tableau, Δt, solver, iguess, params.x, params.y, params.X, params.Y, params.F)
end

"Integrate ODE with fully implicit Runge-Kutta integrator."
function integrate!{T, FT, ST, IT}(int::IntegratorFIRK{T, FT, ST, IT}, sol::SolutionODE{T})
    # loop over initial conditions
    for m in 1:sol.n0
        # copy initial conditions from solution
        simd_copy_xy_first!(int.x, sol, 0, m)

        # initialise initial guess
        initialize!(int.iguess, sol.t[0], int.x)

        for n in 1:sol.ntime
            # set time for nonlinear solver
            int.solver.Fparams.t = sol.t[n]

            # copy previous solution to initial guess
            update!(int.iguess, sol.t[n], int.x)

            # compute initial guess for internal stages
            for i in 1:int.tableau.s
                evaluate(int.iguess, int.y, int.tableau.c[i])
                for k in 1:int.equation.d
                    int.solver.x[int.equation.d*(i-1)+k] = int.y[k]
                end
            end

            # call nonlinear solver
            solve!(int.solver)

            if !solverStatusOK(int.solver.status, int.solver.params)
                println(int.solver.status)
            end

            # compute final update
            simd_mult!(int.y, int.F, int.tableau.b)
            simd_axpy!(int.Δt, int.y, int.x)

            # copy to solution
            if mod(n, sol.nsave) == 0
                simd_copy_yx_first!(int.x, sol, div(n, sol.nsave), m)
            end
        end
    end
    nothing
end

"Integrate partitioned ODE with fully implicit Runge-Kutta integrator."
function integrate!(int::IntegratorFIRK, s::SolutionPODE)
    # TODO
end


text:
function test3()
  # copy test3 files to temp so we can test writing to them
  cp("test3.asc","temp\\test3.asc",force = true)
  cp("test3.log","temp\\test3.log",force = true)

  filename = "temp\\test3.asc"
  # exectuablepath = null string will not run LTspice.exe.  Test parsing only.
  sim = LTspiceSimulation(filename,executablepath="")
  @test LTspice.does_circuitfilearray_file_match(sim)
  show(IOBuffer(),sim)

  parameternamesverify = ["a","b","c","d","e","f","g","h","i","j","k","l"]
  for (verify,parameter,value) in zip(parameternamesverify,parameternames(sim),parametervalues(sim))
    @test parameter == verify
    @test value == sim[parameter]
  end

  # parameters
  @test sim["a"]≈10.0
  @test sim["b"]≈8.0
  @test sim["c"]≈2.0
  @test sim["d"]≈100.0
  @test sim["e"]≈2.6e-12
  @test sim["f"]≈1.0e7
  @test sim["g"]≈1.0e-9
  @test sim["h"]≈1.0e-15
  @test sim["i"]≈0.005
  @test sim["j"]≈1.276e6
  @test sim["k"]≈-3.45e6
  @test sim["l"]≈454.5

  # measurements
  @test sim["x"]≈1.00394
  @test sim["z"]≈0.019685
  @test sim["y"]≈0.984252

  for key in parameternames(sim)
    sim[key] = 1.0
  end
  dummyread = sim["x"]
  sim_b = LTspiceSimulation(filename, executablepath="")
  for key in parameternames(sim_b)
    @test sim_b[key]≈1.0
  end

  for key in parameternames(sim)
    sim[key] = 2.0
  end
  dummyread = measurementvalues(sim)
  sim_b = LTspiceSimulation(filename, executablepath="")
  for key in parameternames(sim_b)
    @test sim_b[key]≈2.0
  end
  show(IOBuffer(),sim_b)
end
test3()


text:
function ZipJarCrawler(arg0::ClassLoader, arg1::JString)
    return ZipJarCrawler((ClassLoader, JString), arg0, arg1)
end

function ZipJarCrawler(arg0::File)
    return ZipJarCrawler((File,), arg0)
end

function ZipJarCrawler(arg0::JString)
    return ZipJarCrawler((JString,), arg0)
end

function ZipJarCrawler(arg0::URL)
    return ZipJarCrawler((URL,), arg0)
end

function feed(obj::ZipJarCrawler, arg0::Pattern, arg1::DataLoader, arg2::DataProvidersManager)
    return jcall(obj, "feed", jboolean, (Pattern, DataLoader, DataProvidersManager), arg0, arg1, arg2)
end