text:
struct Automatic; end
const automatic = Automatic()

function _js_array(x::AbstractDict; process=string, placeholder=nothing)
    v = OrderedDict[OrderedDict("key" => key, "val" => i, "id" => "id"*randstring()) for (i, (key, val)) in enumerate(x)]
    placeholder !== nothing && pushfirst!(v, OrderedDict("key" => placeholder, "val" => 0, "id" => "id"*randstring()))
    return v
end

function _js_array(x::AbstractArray; process=string, placeholder=nothing)
    v = OrderedDict[OrderedDict("key" => process(val), "val" => i, "id" => "id"*randstring()) for (i, val) in enumerate(x)]
    placeholder !== nothing && pushfirst!(v, OrderedDict("key" => placeholder, "val" => 0, "id" => "id"*randstring()))
    return v
end

function _js_array(o::AbstractObservable; process=string, placeholder=nothing)
    map(t -> _js_array(t; process=process, placeholder=placeholder), o)
end

struct Vals2Idxs{T} <: AbstractVector{T}
    vals::Vector{T}
    vals2idxs::Dict{T, Int}
    function Vals2Idxs(v::AbstractArray{T}) where {T}
        vals = convert(Vector{T}, v)
        idxs = 1:length(vals)
        vals2idxs = Dict{T, Int}(zip(vals, idxs))
        new{T}(vals, vals2idxs)
    end
end

Vals2Idxs(v::AbstractDict) = Vals2Idxs(collect(values(v)))

Base.parent(d::Vals2Idxs) = d.vals

Base.get(d::Vals2Idxs, key, default = 0) = get(d.vals2idxs, key, default)
Base.get(d::Vals2Idxs, key::Integer, default = 0) = get(d.vals2idxs, key, default)
getmany(d::Vals2Idxs{T}, key::AbstractArray{<:T}, default = 0) where {T} =
    filter(t -> t!= 0, map(x -> get(d, x), key))

Base.getindex(d::Vals2Idxs, x::Int) = get(d.vals, x, nothing)

Base.size(d::Vals2Idxs) = size(d.vals)

function valueindexpair(value, vals2idxs, args...; multiple = false, rev = false)
    _get = multiple ? getmany : get
    f = x -> _get(vals2idxs[], x)
    g = x -> getindex(vals2idxs[], x)
    p = ObservablePair(value, args..., f=f, g=g)
    on(vals2idxs) do x
        p.excluded[rev+1](p[rev+1][])
    end
    p
end

function initvalueindex(value, index, vals2idxs;
    multiple = false, default = multiple ? eltype(vals2idxs[])[] : first(vals2idxs[]), rev = false)

    if value === automatic
        value = (index === nothing) ? default : vals2idxs[][Observables.to_value(index)]
    end
    (value isa AbstractObservable) || (value = Observable{Any}(value))
    if index === nothing
        p = valueindexpair(value, vals2idxs; multiple = multiple, rev = rev)
        index = p.second
    else
        (index isa AbstractObservable) || (index = Observable{Any}(index))
        p = valueindexpair(value, vals2idxs, index; multiple = multiple, rev = rev)
    end
    return p
end

"""
```
dropdown(options::AbstractDict;
         value = first(values(options)),
         label = nothing,
         multiple = false)
```

A dropdown menu whose item labels are the keys of options.
If `multiple=true` the observable will hold an array containing the values
of all selected items
e.g. `dropdown(OrderedDict("good"=>1, "better"=>2, "amazing"=>9001))`

`dropdown(values::AbstractArray; kwargs...)`

`dropdown` with labels `string.(values)`
see `dropdown(options::AbstractDict; ...)` for more details
"""
dropdown(T::WidgetTheme, options; kwargs...) =
    dropdown(T::WidgetTheme, Observable{Any}(options); kwargs...)

"""
```
dropdown(options::AbstractObservable;
         value = first(values(options[])),
         label = nothing,
         multiple = false)
```

A dropdown menu whose `options` are a given `Observable`. Set the `Observable` to some other
value to update the options in real time.

## Examples

```julia
options = Observable(["a", "b", "c"])
wdg = dropdown(options)
options[] = ["c", "d", "e"]
```

Note that the `options` can be modified from the widget directly:

```julia
wdg[:options][] = ["c", "d", "e"]
```
"""
function dropdown(::WidgetTheme, options::AbstractObservable;
    attributes=PropDict(),
    placeholder = nothing,
    label = nothing,
    multiple = false,
    value = automatic,
    index = nothing,
    className = "",
    style = PropDict(),
    div_select = nothing,
    kwargs...)

    div_select !== nothing && warn("`div_select` keyword is deprecated", once=true)
    multiple && (attributes[:multiple] = true)
    vals2idxs = map(Vals2Idxs, options)
    p = initvalueindex(value, index, vals2idxs, multiple = multiple)
    value, index = p.first, p.second

    bind = multiple ? "selectedOptions" : "value"
    option_array = _js_array(options, placeholder=placeholder)
    disablePlaceholder =
        js"""
        function(option, item) {
            ko.applyBindingsToNode(option, {disable: item.val == 0}, item);
        }
        """

    attrDict = merge(
        Dict(Symbol("data-bind") => "options : options_js, $bind : index, optionsText: 'key', optionsValue: 'val', valueAllowUnset: true, optionsAfterRender: disablePlaceholder"),
        attributes
    )

    className = mergeclasses(getclass(:dropdown, multiple), className)
    div_select === nothing && (div_select = node(:div, className = className))
    template = node(:select; attributes = attrDict, kwargs...)() |> div_select
    label != nothing && (template = vbox(label, template))
    ui = knockout(template, ["index" => index, "options_js" => option_array];
        methods = ["disablePlaceholder" => disablePlaceholder])
    slap_design!(ui)
    Widget{:dropdown}(["options"=>options, "index" => ui["index"]], scope = ui, output = value, layout = node(:div, className = "field interact-widget")∘Widgets.scope)
end

multiselect(T::WidgetTheme, options; kwargs...) =
    multiselect(T, Observable{Any}(options); kwargs...)

function multiselect(T::WidgetTheme, options::AbstractObservable; container=node(:div, className=:field), wrap=identity,
    label = nothing, typ="radio", wdgtyp=typ, stack=true, skip=1em, hskip=skip, vskip=skip,
    value = automatic, index = nothing, kwargs...)
    vals2idxs = map(Vals2Idxs, options)
    p = initvalueindex(value, index, vals2idxs, multiple = (typ != "radio"))
    value, index = p.first, p.second

    s = gensym()
    option_array = _js_array(options)
    entry = wrap(InteractBase.entry(s; typ=typ, wdgtyp=wdgtyp, stack=stack, kwargs...))
    (entry isa Tuple )|| (entry = (entry,))
    template = container(attributes = Dict("data-bind" => "foreach : options_js"))(
        entry...
    )
    ui = knockout(template, ["index" => index, "options_js" => option_array])
    if (label != nothing)
        ui.dom = stack ? vbox(label, CSSUtil.vskip(vskip), ui.dom) : hbox(label, CSSUtil.hskip(hskip), ui.dom)
    end
    slap_design!(ui)
    Widget{:radiobuttons}(["options"=>options, "index" => ui["index"]], scope = ui, output = value, layout = node(:div, className = "field interact-widget")∘Widgets.scope)
end

function entry(T::WidgetTheme, s; className="", typ="radio", wdgtyp=typ, stack=(typ!="radio"), kwargs...)
    className = mergeclasses(getclass(:input, wdgtyp), className)
    f = stack ? node(:div, className="field") : tuple
    f(
        node(:input, className = className, attributes = Dict("name" => s, "type" => typ, "data-bind" => "checked : \$root.index, checkedValue: val, attr : {id : id}"))(),
        node(:label, attributes = Dict("data-bind" => "text : key, attr : {for : id}"))
    )
end

"""
```
radiobuttons(options::AbstractDict;
             value::Union{T, Observable} = first(values(options)))
```

e.g. `radiobuttons(OrderedDict("good"=>1, "better"=>2, "amazing"=>9001))`

`radiobuttons(values::AbstractArray; kwargs...)`

`radiobuttons` with labels `string.(values)`
see `radiobuttons(options::AbstractDict; ...)` for more details

```
radiobuttons(options::AbstractObservable; kwargs...)
```

Radio buttons whose `options` are a given `Observable`. Set the `Observable` to some other
value to update the options in real time.

## Examples

```julia
options = Observable(["a", "b", "c"])
wdg = radiobuttons(options)
options[] = ["c", "d", "e"]
```

Note that the `options` can be modified from the widget directly:

```julia
wdg[:options][] = ["c", "d", "e"]
```
"""
radiobuttons(T::WidgetTheme, vals; kwargs...) =
    multiselect(T, vals; kwargs...)

"""
```
checkboxes(options::AbstractDict;
         value = valtype(options)[]
```

A list of checkboxes whose item labels are the keys of options.
The observable will hold an array containing the values
of all selected items,
e.g. `checkboxes(OrderedDict("good"=>1, "better"=>2, "amazing"=>9001))`

`checkboxes(values::AbstractArray; kwargs...)`

`checkboxes` with labels `string.(values)`
see `checkboxes(options::AbstractDict; ...)` for more details

```
checkboxes(options::AbstractObservable; kwargs...)
```

Checkboxes whose `options` are a given `Observable`. Set the `Observable` to some other
value to update the options in real time.

## Examples

```julia
options = Observable(["a", "b", "c"])
wdg = checkboxes(options)
options[] = ["c", "d", "e"]
```

Note that the `options` can be modified from the widget directly:

```julia
wdg[:options][] = ["c", "d", "e"]
```

To create checkboxes that are already checked off when shown, one can use the `value` keyword:

```julia
options = ["a", "b", "c"]
value = ["a", "c"]
checkboxes(options, value = value)
```

The boxes "a" and "c" will be checked off when the checkboxes widget is shown.
Both `options` and `value` can be Observables.
"""
checkboxes(T::WidgetTheme, options; kwargs...) =
    Widget{:checkboxes}(multiselect(T, options; typ="checkbox", kwargs...))

"""
```
toggles(options::AbstractDict;
         value = first(values(options)))
```

A list of toggle switches whose item labels are the keys of options.
Tthe observable will hold an array containing the values
of all selected items,
e.g. `toggles(OrderedDict("good"=>1, "better"=>2, "amazing"=>9001))`

`toggles(values::AbstractArray; kwargs...)`

`toggles` with labels `string.(values)`
see `toggles(options::AbstractDict; ...)` for more details

```
toggles(options::AbstractObservable; kwargs...)
```

Toggles whose `options` are a given `Observable`. Set the `Observable` to some other
value to update the options in real time.

## Examples

```julia
options = Observable(["a", "b", "c"])
wdg = toggles(options)
options[] = ["c", "d", "e"]
```

Note that the `options` can be modified from the widget directly:

```julia
wdg[:options][] = ["c", "d", "e"]
```
"""
toggles(T::WidgetTheme, options; kwargs...) =
    Widget{:toggles}(multiselect(T, options; typ="checkbox", wdgtyp="toggle", kwargs...))

for wdg in [:togglebuttons, :tabs]
    @eval $wdg(T::WidgetTheme, options; kwargs...) = $wdg(T::WidgetTheme, Observable(options); kwargs...)
end

"""
`togglebuttons(options::AbstractDict; value::Union{T, Observable}, multiple=false)`

Creates a set of toggle buttons whose labels are the keys of options. Set `multiple=true`
to allow multiple (or zero) active buttons at the same time.

`togglebuttons(values::AbstractArray; kwargs...)`

`togglebuttons` with labels `string.(values)`
see `togglebuttons(options::AbstractDict; ...)` for more details

```
togglebuttons(options::AbstractObservable; kwargs...)
```

Togglebuttons whose `options` are a given `Observable`. Set the `Observable` to some other
value to update the options in real time.

## Examples

```julia
options = Observable(["a", "b", "c"])
wdg = togglebuttons(options)
options[] = ["c", "d", "e"]
```

Note that the `options` can be modified from the widget directly:

```julia
wdg[:options][] = ["c", "d", "e"]
```
"""
function togglebuttons(T::WidgetTheme, options::AbstractObservable;
    className = "",
    activeclass = getclass(:button, "active"),
    multiple = false,
    index = nothing, value = automatic,
    container = node(:div, className = getclass(:togglebuttons)), wrap=identity,
    label = nothing, readout = false, vskip = 1em, kwargs...)

    vals2idxs = map(Vals2Idxs, options)
    p = initvalueindex(value, index, vals2idxs; multiple = multiple)
    value, index = p.first, p.second

    className = mergeclasses("interact-widget", getclass(:button), className)

    updateMethod = multiple ? js"""
    function (val) {
        var id = this.index.indexOf(val);
        id > -1 ? this.index.splice(id, 1) : this.index.push(val);
    }
    """ : js"function (val) {this.index(val)}"
    active = multiple ? "\$root.index().includes(val)" : "\$root.index() == val"

    updateSelected = js_lambda("\$root.update(val)")

    btn = node(:span,
        node(:label, attributes = Dict("data-bind" => "text : key")),
        attributes=Dict("data-bind"=>
        "click: $updateSelected, css: {'$activeclass' : $active, '$className' : true}"),
    )
    option_array = _js_array(options)
    template = container(attributes = "data-bind" => "foreach : options_js")(wrap(btn))

    label != nothing && (template = flex_row(wdglabel(label), template))
    ui = knockout(template, ["index" => index, "options_js" => option_array], methods = Dict("update" => updateMethod))
    slap_design!(ui)

    w = Widget{:togglebuttons}(["options"=>options, "index" => ui["index"], "vals2idxs" => vals2idxs];
                               scope = ui,
                               output = value,
                               layout = t -> div(Widgets.scope(t), className = "interact-widget"))
    if readout
        w[:display] = mask(map(parent, vals2idxs); index = index)
        w.layout = t -> div(Widgets.scope(t), CSSUtil.vskip(vskip), t[:display], className = "interact-widget")
    end
    w
end

"""
`tabs(options::AbstractDict; value::Union{T, Observable})`

Creates a set of tabs whose labels are the keys of options. The label can be a link.

`tabs(values::AbstractArray; kwargs...)`

`tabs` with labels `values`
see `tabs(options::AbstractDict; ...)` for more details

```
tabs(options::AbstractObservable; kwargs...)
```

Tabs whose `options` are a given `Observable`. Set the `Observable` to some other
value to update the options in real time.

## Examples

```julia
options = Observable(["a", "b", "c"])
wdg = tabs(options)
options[] = ["c", "d", "e"]
```

Note that the `options` can be modified from the widget directly:

```julia
wdg[:options][] = ["c", "d", "e"]
```
"""
function tabs(T::WidgetTheme, options::AbstractObservable;
    className = "",
    activeclass = getclass(:tab, "active"),
    index = nothing, value = automatic,
    container=node(:div, className = getclass(:tabs)), wrap=identity,
    label = nothing, readout = false, vskip = 1em, kwargs...)

    vals2idxs = map(Vals2Idxs, options)
    p = initvalueindex(value, index, vals2idxs; default = first(vals2idxs[]))
    value, index = p.first, p.second

    className = mergeclasses("interact-widget", getclass(:tab), className)
    updateSelected = js_lambda("\$root.index(val)")
    tab = node(:li,
        wrap(node(:a, attributes = Dict("data-bind" => "text: key"))),
        attributes=Dict("data-bind"=>
        "click: $updateSelected, css: {'$activeclass' : \$root.index() == val, '$className' : true}"),
    )
    option_array = _js_array(options)
    template = container(
        node(:ul, attributes = "data-bind" => "foreach : options_js")(tab)
    )

    label != nothing && (template = flex_row(wdglabel(label), template))
    ui = knockout(template, ["index" => index, "options_js" => option_array])
    slap_design!(ui)

    w = Widget{:tabs}(["options"=>options, "index" => ui["index"], "vals2idxs" => vals2idxs];
        scope = ui, output = value, layout = t -> div(Widgets.scope(t), className = "interact-widget"))
    if readout
        w[:display] = mask(map(parent, vals2idxs); index = index)
        w.layout = t -> div(Widgets.scope(t), CSSUtil.vskip(vskip), t[:display], className = "interact-widget")
    end
    w
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct OSProfile <: SwaggerModel
    computerName::Any # spec type: Union{ Nothing, String } # spec name: computerName
    adminUsername::Any # spec type: Union{ Nothing, String } # spec name: adminUsername
    adminPassword::Any # spec type: Union{ Nothing, String } # spec name: adminPassword
    customData::Any # spec type: Union{ Nothing, String } # spec name: customData
    windowsConfiguration::Any # spec type: Union{ Nothing, WindowsConfiguration } # spec name: windowsConfiguration
    linuxConfiguration::Any # spec type: Union{ Nothing, LinuxConfiguration } # spec name: linuxConfiguration
    secrets::Any # spec type: Union{ Nothing, Vector{VaultSecretGroup} } # spec name: secrets

    function OSProfile(;computerName=nothing, adminUsername=nothing, adminPassword=nothing, customData=nothing, windowsConfiguration=nothing, linuxConfiguration=nothing, secrets=nothing)
        o = new()
        validate_property(OSProfile, Symbol("computerName"), computerName)
        setfield!(o, Symbol("computerName"), computerName)
        validate_property(OSProfile, Symbol("adminUsername"), adminUsername)
        setfield!(o, Symbol("adminUsername"), adminUsername)
        validate_property(OSProfile, Symbol("adminPassword"), adminPassword)
        setfield!(o, Symbol("adminPassword"), adminPassword)
        validate_property(OSProfile, Symbol("customData"), customData)
        setfield!(o, Symbol("customData"), customData)
        validate_property(OSProfile, Symbol("windowsConfiguration"), windowsConfiguration)
        setfield!(o, Symbol("windowsConfiguration"), windowsConfiguration)
        validate_property(OSProfile, Symbol("linuxConfiguration"), linuxConfiguration)
        setfield!(o, Symbol("linuxConfiguration"), linuxConfiguration)
        validate_property(OSProfile, Symbol("secrets"), secrets)
        setfield!(o, Symbol("secrets"), secrets)
        o
    end
end # type OSProfile

const _property_map_OSProfile = Dict{Symbol,Symbol}(Symbol("computerName")=>Symbol("computerName"), Symbol("adminUsername")=>Symbol("adminUsername"), Symbol("adminPassword")=>Symbol("adminPassword"), Symbol("customData")=>Symbol("customData"), Symbol("windowsConfiguration")=>Symbol("windowsConfiguration"), Symbol("linuxConfiguration")=>Symbol("linuxConfiguration"), Symbol("secrets")=>Symbol("secrets"))
const _property_types_OSProfile = Dict{Symbol,String}(Symbol("computerName")=>"String", Symbol("adminUsername")=>"String", Symbol("adminPassword")=>"String", Symbol("customData")=>"String", Symbol("windowsConfiguration")=>"WindowsConfiguration", Symbol("linuxConfiguration")=>"LinuxConfiguration", Symbol("secrets")=>"Vector{VaultSecretGroup}")
Base.propertynames(::Type{ OSProfile }) = collect(keys(_property_map_OSProfile))
Swagger.property_type(::Type{ OSProfile }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_OSProfile[name]))}
Swagger.field_name(::Type{ OSProfile }, property_name::Symbol) =  _property_map_OSProfile[property_name]

function check_required(o::OSProfile)
    true
end

function validate_property(::Type{ OSProfile }, name::Symbol, val)
end

ext: jl


text:
function addSavingRAM_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

    #List of nonempty arrays
    l = []

    if length(abm.declaredSymbols["Local"])>0
        push!(l,:(Core.Array(localV)[1:N,:]))
    else
        push!(l,:(Core.Array{Float64,2}(undef,0,2)))
    end

    if length(abm.declaredSymbols["Identity"])>0
        push!(l,:(Core.Array(identityV)[1:N,:]))
    else
        push!(l,:(Core.Array{Int,2}(undef,0,2)))
    end

    if length(abm.declaredSymbols["Global"])>0
        push!(l,:(Core.Array(globalV)))
    else
        push!(l,:(Core.Array{Float64,1}()))
    end

    if length(abm.declaredSymbols["GlobalArray"])>0
        list = string("[copy(",abm.declaredSymbols["GlobalArray"][1],")")
        for i in abm.declaredSymbols["GlobalArray"][2:end]
            list = string(list,",copy(",i,")")
        end
        list = string(list,"]")
        push!(l,:(Core.Array($list)))
    elseif length(abm.declaredSymbols["GlobalArray"]) > 0
        list = string("[copy(",abm.declaredSymbols["GlobalArray"][1],"")
        list = string(list,")]")
        push!(l,Meta.parse(string("Core.Array(",list,")")))
    else
        push!(l,:(Core.Array{Core.Array{Float64},1}()))
    end

    if length(abm.declaredSymbols["Medium"]) > 0
        g = []
        if abm.dims >= 1
            if space.medium[1].minBoundaryType == "Periodic"
                push!(g,:(2:Nx_-1))
            else
                push!(g,:(1:Nx_))
            end
        end
        if abm.dims >= 2
            if space.medium[2].minBoundaryType == "Periodic"
                push!(g,:(2:Ny_-1))
            else
                push!(g,:(1:Ny_))
            end
        end
        if abm.dims >= 3
            if space.medium[3].minBoundaryType == "Periodic"
                push!(g,:(2:Nz_-1))
            else
                push!(g,:(1:Nz_))
            end
        end

        push!(l,:(Core.Array(mediumV)[$(g...),:]))
    else
        push!(l,:(Core.Array{Float64,1}()))
    end

    push!(p.declareVar.args,:(commRAM_ = CommunityInTime()))

    push!(p.execInit.args,
        :(begin
            ob = Community($(abm.dims),t,N,com.declaredSymbols_,$(l...))
            push!(commRAM_,ob)
        end)
    )
    push!(p.execInloop.args,
        :(begin
            if t >= tSave
                tSave += dtSave
                ob = Community($(abm.dims),t+dt,N,com.declaredSymbols_,$(l...))
                push!(commRAM_,ob)
            end
        end)
    )

    # push!(p.execAfter.args,
    # :(begin
    #     ob = Community(t+dt,N,com.declaredSymbols_,$(l...))
    #     push!(commRAM_,ob)
    #     end)
    # )

    push!(p.returning.args,:(commRAM_))

    return
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
module PeriodicDarcyTests

using Gridap
using Test
using LinearAlgebra

u(x) = VectorValue(x[1]*(x[1]-1)*(2x[2]-1.0),-x[2]*(x[2]-1.0)*(2x[1]-1.0))
p(x) = x[2]-0.5
f(x) = VectorValue(x[1]*(x[1]-1)*(2x[2]-1.0),-x[2]*(x[2]-1.0)*(2x[1]-1.0)+1.0)
g(x) = 0.0

domain = (0,1,0,1,0)
order = 1
partition = (4,4)
model = CartesianDiscreteModel(domain,partition; isperiodic=(true,false))

V = FESpace(model,ReferenceFE(raviart_thomas,Float64,order),conformity=:Hdiv,
      dirichlet_tags=collect(1:6))

Q = FESpace(model,ReferenceFE(lagrangian,Float64,order);
      conformity=:L2, constraint=:zeromean)

U = TrialFESpace(V,u)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V, Q])
X = MultiFieldFESpace([U, P])

trian = Triangulation(model)
degree = 2*(order)
dΩ = Measure(trian,degree)
x = get_physical_coordinate(trian)

a((u, p),(v, q)) = ∫( u⋅v - (∇⋅v)*p + q*(∇⋅u) )*dΩ

b(( v, q)) = ∫( v⋅f + q*g)*dΩ

op = AffineFEOperator(a,b,X,Y)
xh = solve(op)
uh, ph = xh

eu = u - uh
ep = p - ph

l2(v) = sqrt(sum(∫(v⋅v)*dΩ))
eu_l2 = l2(eu)
ep_l2 = l2(ep)

tol = 1.0e-9
@test eu_l2 < tol
@test ep_l2 < tol

end # module


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/MarketDamages.jl")

    marketdamages = addmarketdamages(m)

    set_param!(m, :MarketDamages, :rtl_realizedtemperature, readpagedata(m, "test/validationdata/$valdir/rtl_realizedtemperature.csv"))
    set_param!(m, :MarketDamages, :rcons_per_cap_SLRRemainConsumption, readpagedata(m,"test/validationdata/$valdir/rcons_per_cap_SLRRemainConsumption.csv"))
    set_param!(m, :MarketDamages, :rgdp_per_cap_SLRRemainGDP, readpagedata(m,"test/validationdata/$valdir/rgdp_per_cap_SLRRemainGDP.csv"))
    set_param!(m, :MarketDamages, :atl_adjustedtolerableleveloftemprise, readpagedata(m,"test/validationdata/atl_adjustedtolerableleveloftemprise_market.csv"))
    set_param!(m, :MarketDamages, :imp_actualreduction, readpagedata(m,"test/validationdata/$valdir/imp_actualreduction_market.csv"))
    set_param!(m, :MarketDamages, :isatg_impactfxnsaturation, 28.333333333333336)

    p = load_parameters(m)
    p["y_year_0"] = 2015.
    p["y_year"] = Mimi.dim_keys(m.md, :time)
    set_leftover_params!(m, p)

    run(m)

    ## Not a very good test, but damages are based on Burke et al. now anyway
    iref = m[:MarketDamages, :iref_ImpactatReferenceGDPperCap]
    iref_compare = readpagedata(m, "test/validationdata/$valdir/iref_ImpactatReferenceGDPperCap_econ_page.csv")
    @test iref ≈ iref_compare rtol=1e-1

    # rcons_per_cap = m[:MarketDamages, :rcons_per_cap_MarketRemainConsumption]
    # rcons_per_cap_compare = readpagedata(m, "test/validationdata/rcons_per_cap_MarketRemainConsumption.csv")
    # @test rcons_per_cap ≈ rcons_per_cap_compare rtol=1e-1

    # rgdp_per_cap = m[:MarketDamages, :rgdp_per_cap_MarketRemainGDP]
    # rgdp_per_cap_compare = readpagedata(m, "test/validationdata/rgdp_per_cap_MarketRemainGDP.csv")
    # @test rgdp_per_cap ≈ rgdp_per_cap_compare rtol=1e-2
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
using JuMP

function test_view_subarray()
  model = Model()
  @variable(model, x[1:3])
  @NLobjective(model, Min, sum(x[i]^4 for i = 1:3))
  @NLconstraint(model, sum(x[i]^2 for i = 1:3) == 4.0)
  @NLconstraint(model, sum(x[i] for i = 1:3) == 1.0)
  nlp = MathProgNLPModel(model)

  x = rand(5)
  functions = [obj, grad, hess, cons, jac]
  for foo in functions, I = [1:3, 1:2:5, [2;4;5], [4;1;3]]
    @test foo(nlp, x[I]) == foo(nlp, @view x[I])
  end
end

test_view_subarray()


text:
# Copy of RadiativeTransfer.Scattering, had problems including it otherwise

"""
    type AbstractPolarizationType
Abstract Polarization type
"""
abstract type AbstractPolarizationType  end

"""
    struct Stokes_IQUV{FT<:AbstractFloat}
A struct which defines full Stokes Vector ([I,Q,U,V]) RT code
# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct Stokes_IQUV{FT<:AbstractFloat} <: AbstractPolarizationType
    "Number of Stokes components (int)"
    n::Int = 4
    "Vector of length `n` for ... (see eq in Sanghavi )"
    D::Array{FT}  = [1.0, 1.0, -1.0, -1.0]
    "Incoming Stokes vector for scalar only"
    I₀::Array{FT} = [1.0, 0.0, 0.0, 0.0] #assuming completely unpolarized incident stellar radiation
end

"""
    struct Stokes_IQU{FT<:AbstractFloat}
A struct which defines Stokes Vector ([I,Q,U]) RT code
# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct Stokes_IQU{FT<:AbstractFloat} <: AbstractPolarizationType
    "Number of Stokes components (int)"
    n::Int = 3
    "Vector of length `n` for ... (see eq in Sanghavi )"
    D::Array{FT}  = [1.0, 1.0, -1.0]
    "Incoming Stokes vector for scalar only"
    I₀::Array{FT} = [1.0, 0.0, 0.0] #assuming linearly unpolarized incident stellar radiation
end

"""
    struct Stokes_I{FT<:AbstractFloat}
A struct which define scalar I only RT code
# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct Stokes_I{FT<:AbstractFloat} <: AbstractPolarizationType
    "Number of Stokes components (int)"
    n::Int = 1
    "Vector of length `n` for ... (see eq in Sanghavi )"
    D::Array{FT} = [1.0]
    "Incoming Stokes vector for scalar only"
    I₀::Array{FT} = [1.0]
end


text:
function draw_global_panel(gui_data::GuiData, app_data::Model.AppData, sock::Comm.SocClient)
    if CImGui.CollapsingHeader("Global")
        CImGui.PushItemWidth(150)

        gui_data.buffer = Model.prep_field(Model.simulation(app_data.model), 50)
        returned = CImGui.InputText("Active Simulation", gui_data.buffer, 50, CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        # println(gui_data.buffer)
        if returned
            # println("enter: [", gui_data.buffer16, "]")
            # range = findfirst("\0", gui_data.buffer16)
            # value = String(SubString(gui_data.buffer16, 1:(range[1] - 1)))
            # Model.set_duration!(app.model, value)
            Model.set_simulation!(app_data.model, gui_data.buffer)
        end

        CImGui.PopItemWidth()

        CImGui.PushItemWidth(80)

        gui_data.buffer = Model.prep_field(Model.duration(app_data.model), 20)
        returned = CImGui.InputText("Duration", gui_data.buffer, 20, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
        # println("enter: [", gui_data.buffer16, "]")
        # range = findfirst("\0", gui_data.buffer16)
        # value = String(SubString(gui_data.buffer16, 1:(range[1] - 1)))
        # Model.set_duration!(app.model, value)
            Model.set_duration!(app_data.model, gui_data.buffer)
        end
        CImGui.SameLine(200)

        gui_data.buffer = Model.prep_field(Model.time_scale(app_data.model), 10)
        returned = CImGui.InputText("Time Scale", gui_data.buffer, 10, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_time_scale!(app_data.model, gui_data.buffer)
        end
        CImGui.SameLine(400)

        gui_data.buffer = Model.prep_field(Model.spans(app_data.model), 10)
        returned = CImGui.InputText("Spans", gui_data.buffer, 10, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_spans!(app_data.model, gui_data.buffer)
        end

        CImGui.PopItemWidth()
    end
end


text:
using Test
using ClimateMachine
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies
using ClimateMachine.MPIStateArrays
using Impero, Printf, MPI, LinearAlgebra, Statistics, GaussQuadrature
using ImperoHooks
import ImperoHooks: launch_volume_divergence!, launch_interface_divergence!

ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const mpicomm = MPI.COMM_WORLD
const FT = Float64
Ω = Circle(-1,1) × Circle(-1,1) × Circle(-1,1)
dims = ndims(Ω)

ClimateMachine.gpu_allowscalar(true)

# Define Grid: might want to loop over element sizes and polynomial orders
grid = DiscontinuousSpectralElementGrid(Ω, elements = (10,10,10), polynomialorder = (4,4,4), array = ArrayType)

x, y, z = coordinates(grid)
device = array_device(x)
dim = ndims(Ω)
N = round(Int, size(x)[1]^(1/dim)) - 1
dependencies = nothing

# Initialize State s
F1 = MPIStateArray{FT}(mpicomm, ArrayType, size(x)[1], size(x)[2], 1)
F2 = MPIStateArray{FT}(mpicomm, ArrayType, size(x)[1], size(x)[2], 1)
F3 = MPIStateArray{FT}(mpicomm, ArrayType, size(x)[1], size(x)[2], 1)
flux = MPIStateArray{FT}(mpicomm, ArrayType, size(x)[1], size(x)[2] , 3)

v_flux_divergence = F1 .* 0
s_flux_divergence = F1 .* 0

allindices = collect(1:length(s_flux_divergence))
interior   = setdiff(1:length(s_flux_divergence), grid.vmap⁻)

## Test Block 1: Incompressible Flow Field
# Impero should make this nicer
@. F1 = z
@. F2 = x
@. F3 = y + x

flux[:,:,1:1] .= F1
flux[:,:,2:2] .= F2
flux[:,:,3:3] .= F3
analytic_flux_divergence = @. 0 * x

event = launch_volume_divergence!(v_flux_divergence, flux, grid)
wait(event)
event = launch_interface_divergence!(s_flux_divergence, flux, grid)
wait(event)

tol = eps(1e6)
L∞(x) = maximum(abs.(x))
# The divergence operation is : -v_flux_divergence + s_flux_divergence (note the sign convention)
@testset "Incompressible Flow Field" begin
    computed_divergence = -v_flux_divergence + s_flux_divergence
    @test L∞(computed_divergence - analytic_flux_divergence ) < tol
    @test L∞(v_flux_divergence - s_flux_divergence) < tol
    # test to see if only the exterior terms were affected
    @test L∞(s_flux_divergence[interior]) < eps(1.0)
end

## Test Block 2: Compressible Flow Field
@. F1 =  sin(π*x)
@. F2 =  sin(π*z)
@. F3 =  cos(π*y)
ClimateMachine.gpu_allowscalar(true)
flux[:,:,1:1] .= F1
flux[:,:,2:2] .= F2
flux[:,:,3:3] .= F3
analytic_flux_divergence = @. π*cos(π*x)

event = launch_volume_divergence!(v_flux_divergence, flux, grid)
wait(event)
event = launch_interface_divergence!(s_flux_divergence, flux, grid)
wait(event)

tol = 1e-3 # for this test should be a function of polynomial order / element size
# The divergence operation is : -v_flux_divergence + s_flux_divergence (note the sign convention)
@testset "Compressible Flow Field" begin
    computed_divergence = -v_flux_divergence + s_flux_divergence
    @test L∞(computed_divergence - analytic_flux_divergence ) < tol
    # test to see if only the exterior terms were affected
    @test L∞(s_flux_divergence[interior]) < eps(1.0)
end

## Test Block 3: Compressible Flow Field
a = 1
b = 1
c = 1
@. F1 =  a * sin(π*x)
@. F2 =  b * sin(π*y)
@. F3 =  c * cos(π*z)
ClimateMachine.gpu_allowscalar(true)
flux[:,:,1:1] .= F1
flux[:,:,2:2] .= F2
flux[:,:,3:3] .= F3
analytic_flux_divergence = @. a * π*cos(π*x) + b*π*cos(π*y) - c*π*sin(π*z)

event = launch_volume_divergence!(v_flux_divergence, flux, grid)
wait(event)
event = launch_interface_divergence!(s_flux_divergence, flux, grid)
wait(event)

tol = 1e-3 # for this test should be a function of polynomial order / element size
# The divergence operation is : -v_flux_divergence + s_flux_divergence (note the sign convention)
@testset "Compressible Flow Field" begin
    computed_divergence = -v_flux_divergence + s_flux_divergence
    @test L∞(computed_divergence - analytic_flux_divergence ) < tol
    # test to see if only the exterior terms were affected
    @test L∞(s_flux_divergence[interior]) < eps(1.0)
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

abstract type AbstractLShapedExecution end
# Execution API #
# ------------------------------------------------------------
num_thetas(lshaped::AbstractLShaped) = num_thetas(lshaped, lshaped.execution)
initialize_subproblems!(lshaped::AbstractLShaped, scenarioproblems::AbstractScenarioProblems) = initialize_subproblems!(lshaped.execution, scenarioproblems)
mutate_subproblems!(mutator::Function, lshaped::AbstractLShaped) = mutate_subproblems!(mutator, lshaped.execution)
finish_initilization!(lshaped::AbstractLShaped) = finish_initilization!(lshaped, lshaped.execution)
restore_subproblems!(lshaped::AbstractLShaped) = restore_subproblems!(lshaped, lshaped.execution)
solve_master!(lshaped::AbstractLShaped) = solve_master!(lshaped, lshaped.execution)
iterate!(lshaped::AbstractLShaped) = iterate!(lshaped, lshaped.execution)
start_workers!(lshaped::AbstractLShaped) = start_workers!(lshaped, lshaped.execution)
close_workers!(lshaped::AbstractLShaped) = close_workers!(lshaped, lshaped.execution)
resolve_subproblems!(lshaped::AbstractLShaped) = resolve_subproblems!(lshaped, lshaped.execution)
current_decision(lshaped::AbstractLShaped) = current_decision(lshaped, lshaped.execution)
calculate_objective_value(lshaped::AbstractLShaped) = calculate_objective_value(lshaped, lshaped.execution)
timestamp(lshaped::AbstractLShaped) = timestamp(lshaped, lshaped.execution)
incumbent_decision(lshaped::AbstractLShaped, t::Integer, regularizer::AbstractRegularization) = incumbent_decision(lshaped, t, regularizer, lshaped.execution)
incumbent_objective(lshaped::AbstractLShaped, t::Integer, regularizer::AbstractRegularization) = incumbent_objective(lshaped, t, regularizer, lshaped.execution)
incumbent_trustregion(lshaped::AbstractLShaped, t::Integer, regularizer::AbstractRegularization) = incumbent_trustregion(lshaped, t, regularizer, lshaped.execution)
readd_cuts!(lshaped::AbstractLShaped, consolidation::Consolidation) = readd_cuts!(lshaped, consolidation, lshaped.execution)
subobjectives(lshaped::AbstractLShaped) = subobjectives(lshaped, lshaped.execution)
set_subobjectives(lshaped::AbstractLShaped, Qs::AbstractVector) = set_subobjectives(lshaped, Qs, lshaped.execution)
model_objectives(lshaped::AbstractLShaped) = model_objectives(lshaped, lshaped.execution)
set_model_objectives(lshaped::AbstractLShaped, Qs::AbstractVector) = set_model_objectives(lshaped, Qs, lshaped.execution)
fill_submodels!(lshaped::AbstractLShaped, scenarioproblems::AbstractScenarioProblems) = fill_submodels!(lshaped, scenarioproblems, lshaped.execution)

# ------------------------------------------------------------
include("common.jl")
include("serial.jl")
include("distributed.jl")
include("synchronous.jl")
include("asynchronous.jl")


text:
module Signals

using Statistics
include("DSP.jl")

"""

Output xN after adding noise to x

SNR:: in decibel scale
"""
function addnoise!(xN, x, SNR)

	# store zero mean version of x in xN
	xavg=0.0
	for i in eachindex(x)
		xavg += x[i]
	end
	xavg /= length(x)

	for i in eachindex(x)
		xN[i] = x[i] - xavg
	end

	σx=var(xN)

	σxN=sqrt(σx^2*inv(10^(SNR/10.)))

	# factor to be multiplied to each scalar
	α=sqrt(σxN)
	for i in eachindex(x)
		xN[i] = x[i] + α*randn()
	end
end

"""
Construct Toy Green's functions
Decaying peaks, control number of events, and their positions, depending on bfrac and efrac.
* `afrac` : control amplitude of peaks
"""
function toy_reflec_green!(x; bfrac=0.0, c=1.0, afrac=1.0)
	nt=size(x,1)
	nr=size(x,2)
	c=c*sqrt((nr^2)/(nt-1))

	it0=1+round(Int,bfrac*nt)
	for ir in 1:nr
		it=round(Int, it0+(ir*ir)*inv(c)*inv(c))
		if(it<=nt)
			x[it, ir]=1.0*afrac
		end
	end
	return x
end

function toy_direct_green!(x; bfrac=0.0, c=1.0, afrac=1.0)
	nt=size(x,1)
	nr=size(x,2)
	c=c*((nr)/(nt-1))

	it0=1+round(Int,bfrac*nt)
	for ir in 1:nr
		it=round(Int, it0+(ir)*inv(c))
		if(it<=nt)
			x[it, ir]=1.0*afrac
		end
	end
	return x

end

#=

	a=1+round(Int,bfrac*nt)
	b=nt-round(Int,efrac*nt)
	(b-a<nevents+2) && error("not enough samples")

	mult_shift=zeros(Int64,nevents+1)
	itvec=zeros(Int64,nevents+1)
	for ir in 1:nr
		itvec[:]=0
		for ie in 1:nevents
			if(ie==1)
				# first event, direct arrival
				itvec[ie]=rand(DiscreteUniform(a,b-nevents-1))
			else
				# all other events after first event
				itvec[ie]=rand(DiscreteUniform(itvec[1]+1,b-1))
				if(mult_shift[ie]==0)
					mult_shift[ie]=rand(DiscreteUniform(1,b-1-itvec[1]))
				end

			end
			sgn=1.
			# randomly choose sign for all other events
			if(ie≠1)
				sgn=(rand(Bernoulli())==1) ? 1. : -1.
			end
			x[itvec[ie],ir]=sgn*afrac[ie]
		end
	end
	return x
end

=#

end # module


text:
using Symbolics
using RoboDojo
using LinearAlgebra
using SparseArrays

using MeshCat
using GeometryBasics
using Colors
using CoordinateTransformations
using Rotations

# visualize
function set_background!(vis::Visualizer; top_color=RGBA(1,1,1.0), bottom_color=RGBA(1,1,1.0))
    setprop!(vis["/Background"], "top_color", top_color)
    setprop!(vis["/Background"], "bottom_color", bottom_color)
end
vis = Visualizer()
open(vis)
set_background!(vis)

function hat(x)
    return [0 -x[3] x[2];
            x[3] 0 -x[1];
           -x[2] x[1] 0]
end

function L_mult(x)
    [x[1] -transpose(x[2:4]);
     x[2:4] x[1] * I(3) + hat(x[2:4])]
end

# right quaternion multiply as matrix
function R_mult(x)
    [x[1] -transpose(x[2:4]); x[2:4] x[1] * I(3) - hat(x[2:4])]
end

# rotation matrix
function rotation_matrix(q)
    H = [zeros(1, 3); I(3)]
    transpose(H) * L_mult(q) * transpose(R_mult(q)) * H
end

nz = 10
ny = 12
nw = nz + ny
nθ = 12
@variables z[1:nz] y[1:ny] θ[1:nθ] κ[1:1]

function unpack_primals(z)
    pa = z[0 .+ (1:3)]
    pb = z[3 .+ (1:3)]
    ta = z[6 .+ (1:1)]
    tb = z[7 .+ (1:1)]
    sa = z[8 .+ (1:1)]
    sb = z[9 .+ (1:1)]
    return pa, pb, ta, tb, sa, sb
end

function unpack_duals(y)
   λa = y[0 .+ (1:4)]
   λb = y[4 .+ (1:4)]

   νa = y[8 .+ (1:2)]
   νb = y[10 .+ (1:2)]

   return λa, λb, νa, νb
end

function unpack_data(θ)
    pa1 = θ[0 .+ (1:3)]
    pa2 = θ[3 .+ (1:3)]

    pb1 = θ[6 .+ (1:3)]
    pb2 = θ[9 .+ (1:3)]

    return pa1, pa2, pb1, pb2
end

function objective(z, θ)
    pa, pb, ta, tb, sa, sb = unpack_primals(z)
    d = pa - pb
    return dot(d, d) + 1.0e-3 * (ta[1] - 0.5)^2 + 1.0e-3 * (tb[1] - 0.5)^2
end

function constraints(z, θ)
    pa, pb, ta, tb, sa, sb = unpack_primals(z)
    pa1, pa2, pb1, pb2 = unpack_data(θ)

    [
        pa - (pa2 + ta[1] * (pa1 - pa2));
        sa[1] - (1.0 - ta[1]);
        pb - (pb2 + tb[1] * (pb1 - pb2));
        sb[1] - (1.0 - tb[1]);
    ]
end

function lagrangian(z, y, θ)
    # initialize
    L = 0.0

    # primals
    pa, pb, ta, tb, sa, sb = unpack_primals(z)

    # duals
    λa, λb, νa, νb = unpack_duals(y)

    # objective
    J = objective(z, θ)
    L += J

    # constraints
    c = constraints(z, θ)
    L += sum([λa[1] λa[2] λa[3] λa[4] λb[1] λb[2] λb[3] λb[4]] * c)

    # inequalities
    L -= ta[1] * νa[1]
    L -= sa[1] * νa[2]

    L -= tb[1] * νb[1]
    L -= sb[1] * νb[2]

    return L
end

L = lagrangian(z, y, θ)
Lz = Symbolics.gradient(L, z)

function residual(w, θ, κ)
    # primals
    z = w[1:10]
    pa, pb, ta, tb, sa, sb = unpack_primals(z)

    # duals
    y = w[10 .+ (1:12)]
    λa, λb, νa, νb = unpack_duals(y)

    # Lagrangian
    lag = lagrangian(z, y, θ) #+ 1.0e-2 * dot(λa, λa) - 1.0e-2 * dot(νa, νa) - 1.0e-5 * dot(νb, νb)
    lagz = Symbolics.gradient(lag, z)
    con = constraints(z, θ)
    comp = [
            ta[1] * νa[1];
            sa[1] * νa[2];
            tb[1] * νb[1];
            sb[1] * νb[2];
           ]

    res = [
            lagz;
            con;
            comp .- κ;
          ]

    return res
end

w = [z; y]
r = residual([z; y], θ, κ)
rw = Symbolics.jacobian(r, [z; y])
rθ = Symbolics.jacobian(r, θ)

r_func = eval(Symbolics.build_function(r, w, θ, κ)[2])
rw_func = eval(Symbolics.build_function(rw, w, θ)[2])
rθ_func = eval(Symbolics.build_function(rθ, w, θ)[2])

# pre-allocate
r0 = zeros(nw)
rw0 = zeros(nw, nw)
rθ0 = zeros(nw, nθ)

# random
# w0 = randn(nw)
# θ0 = randn(nθ)
# κ0 = randn(1)

# r_func(r0, w0, θ0, κ0)
# rw_func(rw0, w0, θ0)
# rθ_func(rθ0, w0, θ0)

function rw_func_reg(rw, w, θ)
    rw_func(rw, w, θ)
    rw .+= Diagonal([1.0e-5 * ones(10); -1.0e-5 * ones(8); -1.0e-5 * ones(4)])
    return
end

## setup
xa = [-1.0; 0.0; 0.0]
qa = Quaternion(RotY(0.5 * π) * RotX(0.0 * π))
qa = [qa.w; qa.x; qa.y; qa.z]
ra = 0.1
ha = 0.2

pa1 = xa + rotation_matrix(qa) * [0.0; 0.0;  0.5 * ha]
pa2 = xa + rotation_matrix(qa) * [0.0; 0.0; -0.5 * ha]

xb = [1.0; 0.0; 0.0]
qb = Quaternion(RotZ(0.0 * π) * RotY(0.0 * π) * RotX(0.0))
qb = [qb.w; qb.x; qb.y; qb.z]
rb = 0.1
hb = 0.2
pb1 = xb + rotation_matrix(qb) * [0.0; 0.0;  0.5 * hb]
pb2 = xb + rotation_matrix(qb) * [0.0; 0.0; -0.5 * hb]

## initialization
w0 = [xa; xb; 0.5; 0.5; 1.0; 1.0; zeros(8); ones(4)]
θ0 = [pa1; pa2; pb1; pb2]
κ0 = [1.0]

r_func(r0, w0, θ0, κ0)
rw_func_reg(rw0, w0, θ0)
# rw_func(rw0, w0, θ0)
rθ_func(rθ0, w0, θ0)
cond(rw0)
rank(rw0)
rw0 \ r0
rw0 \ rθ0

## solver
idx = RoboDojo.IndicesOptimization(
    nw, nw,
    [collect(7:10), collect(19:22)], [collect(7:10), collect(19:22)],
    Vector{Vector{Vector{Int}}}(), Vector{Vector{Vector{Int}}}(),
    collect(1:18),
    collect(19:22),
    Vector{Int}(),
    Vector{Vector{Int}}(),
    collect(19:22),
)

ip = RoboDojo.interior_point(w0, θ0;
    s = RoboDojo.Euclidean(length(w0)),
    idx = idx,
    r! = r_func, rz! = rw_func_reg, rθ! = rθ_func,
    r  = zeros(idx.nΔ),
    rz = zeros(idx.nΔ, idx.nΔ),
    rθ = zeros(idx.nΔ, length(θ0)),
    opts = RoboDojo.InteriorPointOptions(
            undercut=5.0,
            γ_reg=0.1,
            r_tol=1e-5,
            κ_tol=1e-5,
            max_ls=50,
            ϵ_min=0.05,
            diff_sol=true,
            verbose=true))

RoboDojo.interior_point_solve!(ip)

pa = ip.z[0 .+ (1:3)]
pb = ip.z[3 .+ (1:3)]
dir = (pa - pb) ./ norm(pa - pb)

ip.z[6 .+ (1:1)]
ip.z[7 .+ (1:1)]
ip.z[8 .+ (1:1)]
ip.z[9 .+ (1:1)]

ip.z[18 .+ (1:4)]

# ∂pa∂xa = ip.δz[0 .+ (1:3), 0  .+ (1:3)]
# ∂pa∂qa = ip.δz[0 .+ (1:3), 3  .+ (1:4)]
# ∂pa∂xb = ip.δz[0 .+ (1:3), 7  .+ (1:3)]
# ∂pa∂qb = ip.δz[0 .+ (1:3), 10 .+ (1:4)]

# ∂pb∂xa = ip.δz[3 .+ (1:3), 0  .+ (1:3)]
# ∂pb∂qa = ip.δz[3 .+ (1:3), 3  .+ (1:4)]
# ∂pb∂xb = ip.δz[3 .+ (1:3), 7  .+ (1:3)]
# ∂pb∂qb = ip.δz[3 .+ (1:3), 10 .+ (1:4)]

# visualize
vis = Visualizer()
open(vis)

# capsule a
cyl1 = GeometryBasics.Cylinder(Point(0.0, 0.0, -0.5 * ha), Point(0.0, 0.0, 0.5 * ha), ra)
sph1 = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), ra)
color1 = Colors.RGBA(0.7, 0.7, 0.7, 0.5);
setobject!(vis[:cap1][:cyl], cyl1, MeshPhongMaterial(color=color1))
setobject!(vis[:cap1][:base], sph1, MeshPhongMaterial(color=color1))
settransform!(vis[:cap1][:base], Translation(0.0, 0.0, -0.5 * ha))
setobject!(vis[:cap1][:tip], sph1, MeshPhongMaterial(color=color1))
settransform!(vis[:cap1][:tip], Translation(0.0, 0.0, 0.5 * ha))

# capsule b
cyl2 = GeometryBasics.Cylinder(Point(0.0, 0.0, -0.5 * hb), Point(0.0, 0.0, 0.5 * hb), rb)
sph2 = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), rb)
color2 = Colors.RGBA(0.7, 0.7, 0.7, 0.5);
setobject!(vis[:cap2][:cyl], cyl2, MeshPhongMaterial(color=color2))
setobject!(vis[:cap2][:base], sph2, MeshPhongMaterial(color=color2))
settransform!(vis[:cap2][:base], Translation(0.0, 0.0, -0.5 * hb))
setobject!(vis[:cap2][:tip], sph2, MeshPhongMaterial(color=color2))
settransform!(vis[:cap2][:tip], Translation(0.0, 0.0, 0.5 * hb))

# set configuration
settransform!(vis[:cap1], compose(Translation(xa), LinearMap(Quaternion(qa..., false))))
settransform!(vis[:cap2], compose(Translation(xb), LinearMap(Quaternion(qb..., false))))

# closest points
color_cp = Colors.RGBA(0.0, 0.0, 0.0, 1.0);
cs1 = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), 0.025)
cs2 = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), 0.025)
setobject!(vis[:cs1], cs1, MeshPhongMaterial(color=color_cp))
setobject!(vis[:cs2], cs2, MeshPhongMaterial(color=color_cp))
settransform!(vis[:cs1], Translation(pa - dir * ra))
settransform!(vis[:cs2], Translation(pb + dir * rb))


text:
@testset "pointer (on device)" begin

@testset "unsafe_load & unsafe_store!" begin

@eval struct LoadableStruct
    a::Int64
    b::UInt8
end
Base.one(::Type{LoadableStruct}) = LoadableStruct(1,1)
Base.zero(::Type{LoadableStruct}) = LoadableStruct(0,0)

@testset for T in (Int8, UInt16, Int32, UInt32, Int64, UInt64, Int128,
                   Float32,Float64,
                   LoadableStruct),
             cached in (false, true)
    d_a = Mem.upload(ones(T))
    d_b = Mem.upload(zeros(T))

    ptr_a = CUDAnative.DevicePtr{T,AS.Global}(Base.unsafe_convert(Ptr{T}, d_a))
    ptr_b = CUDAnative.DevicePtr{T,AS.Global}(Base.unsafe_convert(Ptr{T}, d_b))

    @test Mem.download(T, d_a) != Mem.download(T, d_b)
    if cached
        @on_device Base.unsafe_store!($ptr_b, CUDAnative.unsafe_cached_load($ptr_a))
    else
        @on_device Base.unsafe_store!($ptr_b, CUDAnative.unsafe_load($ptr_a))
    end
    @test Mem.download(T, d_a) == Mem.download(T, d_b)
end

end

end


text:
module Baseball

using HTTP
using DataFrames
using CSV
using EzXML
using Dates
using TimeZones
using JSON

include("utils.jl")
include("data/mlb_game_files/mlb_game_files.jl")
include("data/mlb_savant/mlb_savant.jl")
include("data/weather/darksky.jl")
include("atmosphere.jl")

end # module


text:
@recipe function f(p::Passive{T}) where {T}
    seriestype --> :csg
    linecolor --> :grey
    fillcolor --> :grey
    fillalpha --> 0.2
    l = p.name != "" ? p.name : "Passive $(p.id)"
    @series begin
        #add empty line so that label is shown with alpha = 1
        seriestype := :path
        label --> l
        linewidth := 1
        T[], T[]
    end
    label := ""
    p.geometry
end

@recipe function f(sc::Semiconductor{T}) where {T}
    seriestype --> :csg
    linecolor --> :black
    @series begin
        #add empty line so that label is shown with alpha = 1
        seriestype := :path
        label --> "Semiconductor"
        linewidth := 1
        T[], T[]
    end
    label := ""
    sc.geometry
end

@recipe function f(contact::Contact{T}) where {T}
    seriestype --> :csg
    linecolor --> contact.id
    fillcolor --> contact.id
    fillalpha --> 0.2
    l = contact.name != "" ? "$(contact.name) (id: $(contact.id))" : "Contact - id: $(contact.id)"
    @series begin
        #add empty line so that label is shown with alpha = 1
        seriestype := :path
        label --> l
        linewidth := 1
        T[], T[]
    end
    label := ""
    contact.geometry
end

@recipe function f(det::SolidStateDetector; show_semiconductor = false, show_passives = true)

    show_normal --> false

    if show_semiconductor
        @series begin
            det.semiconductor
        end
    end
    for c in det.contacts
        @series begin
            c
        end
    end
    if show_passives && !ismissing(det.passives)
        for p in det.passives
            @series begin
                p
            end
        end
    end
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module Config_ACASX_GM

export defineSimParams, defineSim

using SISLES
using SISLES.GenerativeModel
using AdaptiveStressTesting

function defineSimParams(;encounter_number::Int64 = 1,
                         encounter_seed::UInt64 = UInt64(0),
                         nmac_r::Float64 = 500.0,
                         nmac_h::Float64 = 100.0,
                         max_steps::Int64 = 51,
                         num_aircraft::Int64 = 2,
                         encounter_model::Symbol = :LLCEMDBN, #:LLCEMDBN, :StarDBN
                         encounter_equipage::Symbol = :EvE, #:EvE, :EvU
                         response_model::Symbol = :ICAO, #:ICAO
                         cas_model::Symbol = :CCAS, #:CCAS, :ADD
                         dynamics_model::Symbol = :LLADM, #:LLADM
                         end_on_nmac::Bool = true,
                         encounter_file::AbstractString = Pkg.dir("SISLES/src/Encounter/CorrAEMImpl/params/cor.txt"),
                         command_method::Symbol = :DBN,
                         initial_sample_file::AbstractString = Pkg.dir("RLESCAS/encounters/initial.txt"),
                         transition_sample_file::AbstractString = Pkg.dir("RLESCAS/encounters/transition.txt"),
                         #libcas::AbstractString = Pkg.dir("CCAS/libcas0.8.6/lib/libcas"), #empty if using :ADD
                         #libcas_config::AbstractString = Pkg.dir("CCAS/libcas0.8.6/parameters/0.8.5.standard.r13.xa.config.txt")
                         #libcas::AbstractString = Pkg.dir("CCAS/libcas0.9.0/lib/libcas"),
                         #libcas_config::AbstractString = Pkg.dir("CCAS/libcas0.9.0/parameters/0.9.0.r14.rev2_3_4candidate07_active.config.txt")
                         #libcas::AbstractString = Pkg.dir("CCAS/libcas0.9.2/lib/libcas"),
                         #libcas_config::AbstractString = Pkg.dir("CCAS/libcas0.9.2/parameters/0.9.2.r14.rev3_7candidate08_active.config.txt")
                         #libcas::AbstractString = Pkg.dir("CCAS/libcas0.9.3/lib/libcas"),
                         #libcas_config::AbstractString = Pkg.dir("CCAS/libcas0.9.3/parameters/0.9.3.standard.r14.xa.config.txt")
                         #libcas::AbstractString = Pkg.dir("CCAS/libcas0.10.0/lib/libcas"),
                         #libcas_config::AbstractString = Pkg.dir("CCAS/libcas0.10.0/parameters/0.10.0.standard.r15_pre25iter93.xa.tcas.config.txt")
                         libcas::AbstractString = Pkg.dir("CCAS/libcas0.10.1/lib/libcas"),
                         libcas_config::AbstractString = Pkg.dir("CCAS/libcas0.10.1/parameters/0.10.1.standard.r15mtf.tcas.xa.config.txt"),
                         libcas2::Union{Void,String}=nothing,
                         libcas2_config::Union{Void,String}=nothing
                         )
  p = ACASX_GM_params()

  p.encounter_number = encounter_number
  p.encounter_seed = encounter_seed
  p.nmac_r = nmac_r
  p.nmac_h = nmac_h
  p.max_steps = max_steps
  p.num_aircraft = num_aircraft
  p.response_model = response_model
  p.cas_model = cas_model
  p.encounter_equipage = encounter_equipage
  p.dynamics_model = dynamics_model
  p.end_on_nmac = end_on_nmac
  p.command_method = command_method
  p.encounter_file = encounter_file
  p.initial_sample_file = initial_sample_file
  p.transition_sample_file = transition_sample_file
  p.libcas = libcas
  p.libcas_config = libcas_config

  if libcas2 == nothing || libcas2_config == nothing
      return p #single sim run
  end

  p2 = deepcopy(p)
  p2.libcas = libcas2
  p2.libcas_config = libcas2_config

  return (p, p2) #dual sim run
end

#single sim run
defineSim(p::ACASX_GM_params) = ACASX_GM(p)

#dual sim run
defineSim(p::Tuple{ACASX_GM_params,ACASX_GM_params}) = DualSim(ACASX_GM(p[1]), ACASX_GM(p[2]))

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
export Jackknife, JackknifeSet
export jackknife

type Jackknife <: ScalarObservable
  xs :: Vector{Float64}
end

Jackknife(jk::Jackknife, f::Function) = Jackknife(map(f,jk.xs))

function jk_helper(xs)
  s = sum(xs)
  n = length(xs)-1
  ret = similar(xs)
  for i in 1:n+1
    ret[i] = (s-xs[i])/n
  end
  return ret
end

function Jackknife(o::ScalarObservable)
  return count(o) > 0 ? Jackknife(jk_helper(o.bins)) : Jackknife(zeros(0))
end

count(jk::Jackknife) = length(jk.xs)

function mean(jk::Jackknife)
  if isempty(jk)
    return NaN
  else
    return mean(jk.xs)
  end
end
function stderror(jk::Jackknife)
  n = count(jk)
  if n == 0
    return NaN
  elseif n == 1
    return Inf
  else
    m2 = sumabs2(jk.xs)
    m2 /= n
    m = mean(jk)
    sigma2 = m2 - m*m
    sigma2 *= n-1
    sigma2 = maxzero(sigma2)
    return sqrt(sigma2)
  end
end

unary_functions = (
  :-,
  :sin, :cos, :tan,
  :sind, :cosd, :tand,
  :sinpi, :cospi,
  :sinh, :cosh, :tanh,
  :asin, :acos, :atan,
  :asind, :acosd, :atand,
  :sec, :csc, :cot,
  :secd, :cscd, :cotd,
  :asec, :acsc, :acot,
  :asecd, :acscd, :acotd,
  :sech, :csch, :coth,
  :asinh, :acosh, :atanh,
  :asech, :acsch, :acoth,
  :sinc, :cosc,
  :log, :log2, :log10, :log1p,
  :exp, :exp2, :exp10, :expm1,
  :abs, :abs2,
  :sqrt, :cbrt,
  :erf, :erfc, :erfcx,
  :erfinv, :erfcinv,
  :gamma, :lgamma, :lfact,
  :digamma, :invdigamma, :trigamma,
  :airyai, :airyprime, :airyaiprime,
  :airybi, :airybiprime,
  :besselj0, :besselj1,
  :bessely0, :bessely1,
  :eta, :zeta
)

for op in unary_functions
  eval( Expr(:import, :Base, op) )
  eval( Expr(:export, op) )
  @eval ($op)(jk::Jackknife) = Jackknife(jk, $op)
end

binary_functions = (
  :+, :-, :*, :/, :\
)

for op in binary_functions
  eval( Expr(:import, :Base, op) )
  eval( Expr(:export, op) )
  @eval ($op)(jk::Jackknife, rhs::Real) = Jackknife(jk, lhs->($op)(lhs,rhs))
  @eval ($op)(lhs::Real, jk::Jackknife) = Jackknife(jk, rhs->($op)(lhs,rhs))
  op_bw = symbol("."*string(op))
  @eval ($op)(lhs::Jackknife, rhs::Jackknife) = Jackknife( ($op_bw)(lhs.xs, rhs.xs))
end

typealias JackknifeSet MCObservableSet{Jackknife}

jackknife(obs::ScalarObservable) = Jackknife(obs)
function jackknife(obsset :: MCObservableSet)
  JK = JackknifeSet()
  for (k,v) in obsset
    JK[k] = Jackknife(v)
  end
  return JK
end


text:
int_rules_5_5_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 5.5.1*u*(a+b*arcsec(c*x))^n =#
    @apply_utils Antiderivative(asec(~(c') * ~x), ~x) => ~x * asec(~c * ~x) - (1 / ~c) * Antiderivative(1 / (~x * sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x) <-- FreeQ(~c, ~x)
    @apply_utils Antiderivative(acsc(~(c') * ~x), ~x) => ~x * acsc(~c * ~x) + (1 / ~c) * Antiderivative(1 / (~x * sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x) <-- FreeQ(~c, ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * asec(~(c') * ~x)) ^ ~n, ~x) => (1 / ~c) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * sec(~x) * tan(~x), ~x), ~x, asec(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * acsc(~(c') * ~x)) ^ ~n, ~x) => (-1 / ~c) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * csc(~x) * cot(~x), ~x), ~x, acsc(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * asec(~(c') * ~x)) / ~x, ~x) => -(Subst(Antiderivative((~a + ~b * acos(~x / ~c)) / ~x, ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * acsc(~(c') * ~x)) / ~x, ~x) => -(Subst(Antiderivative((~a + ~b * asin(~x / ~c)) / ~x, ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * asec(~(c') * ~x)), ~x) => ((~d * ~x) ^ (~m + 1) * (~a + ~b * asec(~c * ~x))) / (~d * (~m + 1)) - ((~b * ~d) / (~c * (~m + 1))) * Antiderivative((~d * ~x) ^ (~m - 1) / sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * acsc(~(c') * ~x)), ~x) => ((~d * ~x) ^ (~m + 1) * (~a + ~b * acsc(~c * ~x))) / (~d * (~m + 1)) + ((~b * ~d) / (~c * (~m + 1))) * Antiderivative((~d * ~x) ^ (~m - 1) / sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * asec(~(c') * ~x)) ^ ~n, ~x) => (1 / (~c) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * sec(~x) ^ (~m + 1) * tan(~x), ~x), ~x, asec(~c * ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ(~n) && (IntegerQ(~m) && (GtQ(~n, 0) || LtQ(~m, -1))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * acsc(~(c') * ~x)) ^ ~n, ~x) => (-1 / (~c) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * csc(~x) ^ (~m + 1) * cot(~x), ~x), ~x, acsc(~c * ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ(~n) && (IntegerQ(~m) && (GtQ(~n, 0) || LtQ(~m, -1))))
    @apply_utils Antiderivative((~(a') + ~(b') * asec(~(c') * ~x)) / (~(d') + ~(e') * ~x), ~x) => ((((((~a + ~b * asec(~c * ~x)) * log(1 + ((~e - sqrt(-((~c) ^ 2) * (~d) ^ 2 + (~e) ^ 2)) * ℯ ^ (im * asec(~c * ~x))) / (~c * ~d))) / ~e + ((~a + ~b * asec(~c * ~x)) * log(1 + ((~e + sqrt(-((~c) ^ 2) * (~d) ^ 2 + (~e) ^ 2)) * ℯ ^ (im * asec(~c * ~x))) / (~c * ~d))) / ~e) - ((~a + ~b * asec(~c * ~x)) * log(1 + ℯ ^ (2 * im * asec(~c * ~x)))) / ~e) - (~b / (~c * ~e)) * Antiderivative(log(1 + ((~e - sqrt(-((~c) ^ 2) * (~d) ^ 2 + (~e) ^ 2)) * ℯ ^ (im * asec(~c * ~x))) / (~c * ~d)) / ((~x) ^ 2 * sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x)) - (~b / (~c * ~e)) * Antiderivative(log(1 + ((~e + sqrt(-((~c) ^ 2) * (~d) ^ 2 + (~e) ^ 2)) * ℯ ^ (im * asec(~c * ~x))) / (~c * ~d)) / ((~x) ^ 2 * sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x)) + (~b / (~c * ~e)) * Antiderivative(log(1 + ℯ ^ (2 * im * asec(~c * ~x))) / ((~x) ^ 2 * sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * acsc(~(c') * ~x)) / (~(d') + ~(e') * ~x), ~x) => (((((~a + ~b * acsc(~c * ~x)) * log(1 - (im * (~e - sqrt(-((~c) ^ 2) * (~d) ^ 2 + (~e) ^ 2)) * ℯ ^ (im * acsc(~c * ~x))) / (~c * ~d))) / ~e + ((~a + ~b * acsc(~c * ~x)) * log(1 - (im * (~e + sqrt(-((~c) ^ 2) * (~d) ^ 2 + (~e) ^ 2)) * ℯ ^ (im * acsc(~c * ~x))) / (~c * ~d))) / ~e) - ((~a + ~b * acsc(~c * ~x)) * log(1 - ℯ ^ (2 * im * acsc(~c * ~x)))) / ~e) + (~b / (~c * ~e)) * Antiderivative(log(1 - (im * (~e - sqrt(-((~c) ^ 2) * (~d) ^ 2 + (~e) ^ 2)) * ℯ ^ (im * acsc(~c * ~x))) / (~c * ~d)) / ((~x) ^ 2 * sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x) + (~b / (~c * ~e)) * Antiderivative(log(1 - (im * (~e + sqrt(-((~c) ^ 2) * (~d) ^ 2 + (~e) ^ 2)) * ℯ ^ (im * acsc(~c * ~x))) / (~c * ~d)) / ((~x) ^ 2 * sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x)) - (~b / (~c * ~e)) * Antiderivative(log(1 - ℯ ^ (2 * im * acsc(~c * ~x))) / ((~x) ^ 2 * sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * asec(~(c') * ~x)), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * asec(~c * ~x))) / (~e * (~m + 1)) - (~b / (~c * ~e * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 1) / ((~x) ^ 2 * sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acsc(~(c') * ~x)), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * acsc(~c * ~x))) / (~e * (~m + 1)) + (~b / (~c * ~e * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 1) / ((~x) ^ 2 * sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asec(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * asec(~c * ~x), u, ~x) - ((~b * ~c * ~x) / sqrt((~c) ^ 2 * (~x) ^ 2)) * Antiderivative(SimplifyIntegrand(u / (~x * sqrt((~c) ^ 2 * (~x) ^ 2 - 1)), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~p, 0) || ILtQ(~p + 1 / 2, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acsc(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * acsc(~c * ~x), u, ~x) + ((~b * ~c * ~x) / sqrt((~c) ^ 2 * (~x) ^ 2)) * Antiderivative(SimplifyIntegrand(u / (~x * sqrt((~c) ^ 2 * (~x) ^ 2 - 1)), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~p, 0) || ILtQ(~p + 1 / 2, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asec(~(c') * ~x)) ^ ~(n'), ~x) => -(Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~x / ~c)) ^ ~n) / (~x) ^ (2 * (~p + 1)), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acsc(~(c') * ~x)) ^ ~(n'), ~x) => -(Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~x / ~c)) ^ ~n) / (~x) ^ (2 * (~p + 1)), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asec(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt((~x) ^ 2)) / ~x) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~x / ~c)) ^ ~n) / (~x) ^ (2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~p + 1 / 2) && (GtQ(~e, 0) && LtQ(~d, 0)))))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acsc(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt((~x) ^ 2)) / ~x) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~x / ~c)) ^ ~n) / (~x) ^ (2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~p + 1 / 2) && (GtQ(~e, 0) && LtQ(~d, 0)))))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asec(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt(~d + ~e * (~x) ^ 2)) / (~x * sqrt(~e + ~d / (~x) ^ 2))) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~x / ~c)) ^ ~n) / (~x) ^ (2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~p + 1 / 2) && Not(GtQ(~e, 0) && LtQ(~d, 0)))))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acsc(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt(~d + ~e * (~x) ^ 2)) / (~x * sqrt(~e + ~d / (~x) ^ 2))) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~x / ~c)) ^ ~n) / (~x) ^ (2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~p + 1 / 2) && Not(GtQ(~e, 0) && LtQ(~d, 0)))))
    @apply_utils Antiderivative(~x * (~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asec(~(c') * ~x)), ~x) => ((~d + ~e * (~x) ^ 2) ^ (~p + 1) * (~a + ~b * asec(~c * ~x))) / (2 * ~e * (~p + 1)) - ((~b * ~c * ~x) / (2 * ~e * (~p + 1) * sqrt((~c) ^ 2 * (~x) ^ 2))) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p + 1) / (~x * sqrt((~c) ^ 2 * (~x) ^ 2 - 1)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && NeQ(~p, -1)
    @apply_utils Antiderivative(~x * (~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acsc(~(c') * ~x)), ~x) => ((~d + ~e * (~x) ^ 2) ^ (~p + 1) * (~a + ~b * acsc(~c * ~x))) / (2 * ~e * (~p + 1)) + ((~b * ~c * ~x) / (2 * ~e * (~p + 1) * sqrt((~c) ^ 2 * (~x) ^ 2))) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p + 1) / (~x * sqrt((~c) ^ 2 * (~x) ^ 2 - 1)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && NeQ(~p, -1)
    @apply_utils Antiderivative((~(f') * ~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asec(~(c') * ~x)), ~x) => With([u = IntHide((~f * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * asec(~c * ~x), u, ~x) - ((~b * ~c * ~x) / sqrt((~c) ^ 2 * (~x) ^ 2)) * Antiderivative(SimplifyIntegrand(u / (~x * sqrt((~c) ^ 2 * (~x) ^ 2 - 1)), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~p], ~x) && (IGtQ(~p, 0) && Not(ILtQ((~m - 1) / 2, 0) && GtQ(~m + 2 * ~p + 3, 0)) || (IGtQ((~m + 1) / 2, 0) && Not(ILtQ(~p, 0) && GtQ(~m + 2 * ~p + 3, 0)) || ILtQ((~m + 2 * ~p + 1) / 2, 0) && Not(ILtQ((~m - 1) / 2, 0))))
    @apply_utils Antiderivative((~(f') * ~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acsc(~(c') * ~x)), ~x) => With([u = IntHide((~f * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * acsc(~c * ~x), u, ~x) + ((~b * ~c * ~x) / sqrt((~c) ^ 2 * (~x) ^ 2)) * Antiderivative(SimplifyIntegrand(u / (~x * sqrt((~c) ^ 2 * (~x) ^ 2 - 1)), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~p], ~x) && (IGtQ(~p, 0) && Not(ILtQ((~m - 1) / 2, 0) && GtQ(~m + 2 * ~p + 3, 0)) || (IGtQ((~m + 1) / 2, 0) && Not(ILtQ(~p, 0) && GtQ(~m + 2 * ~p + 3, 0)) || ILtQ((~m + 2 * ~p + 1) / 2, 0) && Not(ILtQ((~m - 1) / 2, 0))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asec(~(c') * ~x)) ^ ~(n'), ~x) => -(Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~x / ~c)) ^ ~n) / (~x) ^ (~m + 2 * (~p + 1)), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (IntegerQ(~m) && IntegerQ(~p)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acsc(~(c') * ~x)) ^ ~(n'), ~x) => -(Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~x / ~c)) ^ ~n) / (~x) ^ (~m + 2 * (~p + 1)), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (IntegerQ(~m) && IntegerQ(~p)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asec(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt((~x) ^ 2)) / ~x) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~x / ~c)) ^ ~n) / (~x) ^ (~m + 2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (IntegerQ(~p + 1 / 2) && (GtQ(~e, 0) && LtQ(~d, 0))))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acsc(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt((~x) ^ 2)) / ~x) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~x / ~c)) ^ ~n) / (~x) ^ (~m + 2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (IntegerQ(~p + 1 / 2) && (GtQ(~e, 0) && LtQ(~d, 0))))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asec(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt(~d + ~e * (~x) ^ 2)) / (~x * sqrt(~e + ~d / (~x) ^ 2))) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~x / ~c)) ^ ~n) / (~x) ^ (~m + 2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (IntegerQ(~p + 1 / 2) && Not(GtQ(~e, 0) && LtQ(~d, 0))))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acsc(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt(~d + ~e * (~x) ^ 2)) / (~x * sqrt(~e + ~d / (~x) ^ 2))) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~x / ~c)) ^ ~n) / (~x) ^ (~m + 2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (IntegerQ(~p + 1 / 2) && Not(GtQ(~e, 0) && LtQ(~d, 0))))))
    @apply_utils Antiderivative(~u * (~(a') + ~(b') * asec(~(c') * ~x)), ~x) => With([v = IntHide(~u, ~x)], Dist(~a + ~b * asec(~c * ~x), v, ~x) - (~b / ~c) * Antiderivative(SimplifyIntegrand(v / ((~x) ^ 2 * sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x), ~x) <-- InverseFunctionFreeQ(v, ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~u * (~(a') + ~(b') * acsc(~(c') * ~x)), ~x) => With([v = IntHide(~u, ~x)], Dist(~a + ~b * acsc(~c * ~x), v, ~x) + (~b / ~c) * Antiderivative(SimplifyIntegrand(v / ((~x) ^ 2 * sqrt(1 - 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x), ~x) <-- InverseFunctionFreeQ(v, ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * asec(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable(~u * (~a + ~b * asec(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * acsc(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable(~u * (~a + ~b * acsc(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
end

max_stars_repo_path: src/rules/5 Inverse trig functions/5.5 Inverse secant/5.5.1 u (a+b arcsec(c x))^n.jl


text:
"""
    add_uniaxial_metal!(structure,id,E,ν,ρ,α,E₂,f,fᵤ)
向structure加入单轴金属材料
# Arguments
- `structure`: Structure类型实例
- `id`: 材料id
- `E::Integer`: 杨氏模量
- `ν::Integer`: 泊松比
- `ρ`: 质量密度
- `α`: 热膨胀系数
- `E₂`: 屈服后强化模量
- `f`: 屈服点
- `fᵤ`: 极限强度
"""
function add_uniaxial_metal!(structure::Structure,id,E=2e11,ν=0.3,ρ=7850,α=1e-7,E₂=1e11,f=0.345,fᵤ=0.420)::Int
    id=string(id)
    if id in keys(structure.materials)
        throw("material with id "*id*" already exists!")
    end
    hid=length(structure.materials)+1
    material=UniAxialMetal(id,hid,E,ν,ρ,α,E₂,f,fᵤ)
    structure.materials[id]=material
    return 0
end

"""
    add_beam_section!(structure,id,sec_type,size1 [,size2 [,size3 [,size4 [,size5 [,size6 [,size7 [,size8]]]]]]])
向structure加入常用beam截面
# Arguments
- `structure`: Structure类型实例
- `id`: 截面id
- `sec_type`: 截面类型，1-I形截面，2-H形截面，3-箱形截面，4-圆管，5-实心圆形，6-实心矩形
- `size1 ～ size8`: 截面尺寸
"""
function add_beam_section!(structure::Structure,id,sec_type,size1=0,size2=0,size3=0,size4=0,size5=0,size6=0,size7=0,size8=0)::Int
    id=string(id)
    if id in keys(structure.sections)
        throw("section with id "*id*" already exists!")
    end
    hid=length(structure.sections)+1
    if SectionType(sec_type)==Enums.GENERAL_SECTION
        h=float(size1)
        b=float(size2)
        structure.sections[id]=RectangleSection(id,hid,h,b)
    elseif SectionType(sec_type)==Enums.ISECTION
        h=float(size1)
        b1=float(size2)
        b2=float(size3)
        tw=float(size4)
        tf1=float(size5)
        tf2=float(size6)
        structure.sections[id]=ISection(id,hid,h,b1,b2,tw,tf1,tf2)
    elseif SectionType(sec_type)==Enums.HSECTION
        h=float(size1)
        b=float(size2)
        tw=float(size3)
        tf=float(size4)
        structure.sections[id]=HSection(id,hid,h,b,tw,tf)
    elseif SectionType(sec_type)==Enums.BOX
        h=float(size1)
        b=float(size2)
        tw=float(size3)
        tf=float(size4)
        structure.sections[id]=BoxSection(id,hid,h,b,tw,tf)
    elseif SectionType(sec_type)==Enums.PIPE
        d=float(size1)
        t=float(size2)
        structure.sections[id]=PipeSection(id,hid,d,t)
    elseif SectionType(sec_type)==Enums.CIRCLE
        d=float(size1)
        structure.sections[id]=CircleSection(id,hid,d)
    else
        throw("sec_type error!")
    end
    return 0
end

"""
    add_general_section!(structure,id,A,I₂,I₃,J,As₂,As₃,W₂,W₃)
向structure实例加入一般beam截面
# Arguments
- `structure`: Structure类型实例
- `id`: 截面id
- `A`: 截面面积
- `I₂`: 绕2轴惯性矩
- `I₃`: 绕3轴惯性矩
- `J`: 抗扭常数
- `As₂`: 2轴抗剪截面
- `As₃`: 3轴抗剪截面
- `W₂`: 2轴抗弯模量
- `W₃`: 3轴抗弯模量
"""
function add_general_section!(structure::Structure,id,A,I₂,I₃,J,As₂,As₃,W₂,W₃)::Int
    id=string(id)
    if id in keys(structure.sections)
        throw("section with id "*id*" already exists!")
    end
    hid=length(structure.sections)+1
    structure.sections[id]=BeamSection(id,hid,A,I₂,I₃,J,As₂,As₃,W₂,W₃)
    return 0
end

"""
    add_node!(structure,id,x,y,z)
向structure实例加入节点
# Arguments
- `structure`: Structure类型实例
- `id`: 节点id
- `x y z`: 节点坐标
"""
function add_node!(structure::Structure,id,x,y,z)::Int
    id=string(id)
    if id in keys(structure.nodes)
        throw("node id "*string(id)*" already existed")
    end
    hid=length(structure.nodes)+1
    node=Node(id,hid,Float64(x),Float64(y),Float64(z))
    structure.nodes[id]=node
    return 0
end

"""
    set_nodal_restraint!(structure,id,u1,u2,u3,r1,r2,r3)
设置节点约束
# Arguments
- `structure`: Structure类型实例
- `id`: 节点id
- `u1,u2,u3,r1,r2,r3::Bool`: 自由度约束，true为约束，false不约束
"""
function set_nodal_restraint!(structure::Structure,id,u1::Bool,u2::Bool,u3::Bool,r1::Bool,r2::Bool,r3::Bool)::Int
    id=string(id)
    if !(id in keys(structure.nodes))
        throw("node id "*string(id)*" does't existed")
    end
    structure.nodes[id].restraints=[u1,u2,u3,r1,r2,r3]
    return 0
end

"""
    set_nodal_spring!(structure,id,u1,u2,u3,r1,r2,r3)
设置节点弹簧
# Arguments
- `structure`: Structure类型实例
- `id`: 节点id
- `u1,u2,u3,r1,r2,r3`::Number`: 弹簧刚度
"""
function set_nodal_spring!(structure::Structure,id,u1,u2,u3,r1,r2,r3)::Int
    id=string(id)
    if !(id in keys(structure.nodes))
        throw("node id "*string(id)*" does't existed")
    end
    structure.nodes[id].spring=Float64.([u1,u2,u3,r1,r2,r3])
    return 0
end

"""
    set_nodal_mass!(structure,id,u1,u2,u3,r1,r2,r3)
设置节点附加质量
# Arguments
- `structure`: Structure类型实例
- `id`: 节点id
- `u1,u2,u3,r1,r2,r3`::Number`: 各自由度集中质量
"""
function set_nodal_mass!(structure::Structure,id,u1,u2,u3,r1,r2,r3)::Int
    id=string(id)
    if !(id in keys(structure.nodes))
        throw("node id "*string(id)*" does't existed")
    end
    structure.nodes[id].mass=Float64.([u1,u2,u3,r1,r2,r3])
    return 0
end

function add_link!(structure::Structure,id,i,j,E,A,ρ)
    hid=length(structure.links)+1
    node1=structure.nodes[findfirst(x->x.id==string(i),structure.nodes)]
    node2=structure.nodes[findfirst(x->x.id==string(j),structure.nodes)]
    link=Link(id,hid,node1,node2,Float64(E),Float64(A),Float64(ρ))
    push!(structure.links,link)
end

function add_cable!(structure::Structure,id,i,j,E,A,ρ)
    hid=length(structure.cables)+1
    node1=structure.nodes[findfirst(x->x.id==string(i),structure.nodes)]
    node2=structure.nodes[findfirst(x->x.id==string(j),structure.nodes)]
    cable=Cable(id,hid,node1,node2,E,A,ρ)
    push!(structure.cables,cable)
end

"""
    add_beam!(structure,id,i,j,mat_id,sec_id)
向structure实例加入beam单元
# Arguments
- `structure`: Structure类型实例
- `id`: 单元id
- `i,j`: 节点id
- `mat_id`: 材料id
- `sec_id`: 截面id
"""
function add_beam!(structure::Structure,id,i,j,mat_id,sec_id)::Int
    id=string(id)
    hid=length(structure.beams)+1
    if id in keys(structure.beams)
        throw("beam id "*string(id)*" already exists!")
    end
    i,j=string(i),string(j)
    if !(i in keys(structure.nodes))
        throw("node with id "*string(i)*" doesn't exist!")
    end
    if !(j in keys(structure.nodes))
        throw("node with id "*string(j)*" doesn't exist!")
    end
    if !(mat_id in keys(structure.materials))
        throw("material with id "*string(mat_id)*" doesn't exist!")
    end
    if !(sec_id in keys(structure.sections))
        throw("section with id "*string(sec_id)*" doesn't exist!")
    end
    node1=structure.nodes[i]
    node2=structure.nodes[j]
    material=structure.materials[mat_id]
    section=structure.sections[sec_id]
    beam=Beam(id,hid,node1,node2,material,section)
    structure.beams[id]=beam
    return 0
end

"""
    set_beam_orient!(structure,id,degree)
设置beam单元旋转指向
# Arguments
- `structure`: Structure类型实例
- `id`: 单元id
- `degree`: 绕单元1轴旋转角度
"""
function set_beam_orient!(structure::Structure,id,degree)::Int
    id=string(id)
    if !(id in keys(structure.beams))
        throw("beam id "*string(id)*" doesn't exists!")
    end
    rad=π/180*degree
    Rₓ=[1 0 0;
        0 cos(rad) -sin(rad);
        0 sin(rad) cos(rad)]
    structure.beams[id].T[1:3,1:3]*=Rₓ
    structure.beams[id].T[4:6,4:6]*=Rₓ
    structure.beams[id].T[7:9,7:9]*=Rₓ
    structure.beams[id].T[10:12,10:12]*=Rₓ
    return 0
end

"""
    set_beam_release!(structure,id,fᵢ₁,fᵢ₂,fᵢ₃,m₁₁,mᵢ₂,mᵢ₃,fⱼ₁,fⱼ₂,fⱼ₃,mⱼ₁,mⱼ₂,mⱼ₃)
设置beam单元自由度释放
# Arguments
- `structure`: Structure类型实例
- `id`: 单元id
- `fᵢ₁,fᵢ₂,fᵢ₃,m₁₁,mᵢ₂,mᵢ₃,fⱼ₁,fⱼ₂,fⱼ₃,mⱼ₁,mⱼ₂,mⱼ₃::Bool`: 释放自由度，true表示释放，false不释放
"""
function set_beam_release!(structure::Structure,id,fᵢ₁::Bool,fᵢ₂::Bool,fᵢ₃::Bool,m₁₁::Bool,mᵢ₂::Bool,mᵢ₃::Bool,fⱼ₁::Bool,fⱼ₂::Bool,fⱼ₃::Bool,mⱼ₁::Bool,mⱼ₂::Bool,mⱼ₃::Bool)::Int
    id=string(id)
    if !(id in keys(structure.beams))
        throw("beam id "*string(id)*" doesn't exists!")
    end
    structure.beams[id].release=[fᵢ₁,fᵢ₂,fᵢ₃,m₁₁,mᵢ₂,mᵢ₃,fⱼ₁,fⱼ₂,fⱼ₃,mⱼ₁,mⱼ₂,mⱼ₃]
    return 0
end

"""
    add_quad!(structure,id,i,j,k,l,mat_id,t)
向structure实例加入quad单元
# Arguments
- `structure`: Structure类型实例
- `id`: 单元id
- `i,j,k,l`: 节点id
- `mat_id`: 材料id
- `t`: 截面厚度
"""
function add_quad!(structure::Structure,id,i,j,k,l,mat_id,t;t2=0,t3=0,t4=0,elm_type="DKGQ",mass_type="concentrate")::Int
    id=string(id)
    hid=length(structure.quads)+1
    if id in keys(structure.quads)
        throw("quad id "*string(id)*" already exists!")
    end
    i,j,k,l=string(i),string(j),string(k),string(l)
    for n_i in (i,j,k,l)
        if !(n_i in keys(structure.nodes))
            throw("node with id "*n_i*" doesn't exist!")
        end
    end
    if !(string(mat_id) in keys(structure.materials))
        throw("material with id "*string(mat_id)*" doesn't exist!")
    end
    node1=structure.nodes[i]
    node2=structure.nodes[j]
    node3=structure.nodes[k]
    node4=structure.nodes[l]
    material=structure.materials[mat_id]
    quad=Quad(id,hid,node1,node2,node3,node4,material,t,t2,t3,t4,elm_type,mass_type)
    structure.quads[id]=quad
    return 0
end

"""
    set_quad_rotation!(structure,id,degree)
设置单元局部坐标系转角
# Arguments
- `structure`: Structure类型实例
- `id`: 单元id
- `degree`: 旋转角度（ᵒ）
"""
function set_quad_rotation!(structure::Structure,id,degree)::Int
    id=string(id)
    if !(id in keys(structure.quads))
        throw("quad id "*string(id)*" doesn't exists!")
    end
    rad=π/180*degree
    Rₓ=[cos(rad) -sin(rad) 0;
        sin(rad) cos(rad) 0;
        0 0 1]
    for i in 1:4
        structure.quads[id].T[6i-5:6i-3,6i-5:6i-3]*=Rₓ
        structure.quads[id].T[6i-2:6i,6i-2:6i]*=Rₓ
    end
    return 0
end

"""
    add_tria!(structure,id,i,j,k,mat_id,t)
向structure实例加入Tria单元
# Arguments
- `structure`: Structure类型实例
- `id`: 单元id
- `i,j,k`: 节点id
- `mat_id`: 材料id
- `t`: 截面厚度
"""
function add_tria!(structure::Structure,id,i,j,k,mat_id,t;t2=0,t3=0,t4=0,elm_type="DKGQ",mass_type="concentrate")::Int
    id=string(id)
    hid=length(structure.trias)+1
    if id in keys(structure.trias)
        throw("tria id "*string(id)*" already exists!")
    end
    i,j,k=string(i),string(j),string(k)
    for n_i in (i,j,k)
        if !(n_i in keys(structure.nodes))
            throw("node with id "*n_i*" doesn't exist!")
        end
    end
    if !(mat_id in keys(structure.materials))
        throw("material with id "*string(j)*" doesn't exist!")
    end
    node1=structure.nodes[i]
    node2=structure.nodes[j]
    node3=structure.nodes[k]
    material=structure.materials[mat_id]
    tria=Tria(id,hid,node1,node2,node3,material,t,t2,t3,elm_type,mass_type)
    structure.trias[id]=tria
    return 0
end

"""
    set_tria_rotation!(structure,id,degree)
设置三角单元局部坐标系转角
# Arguments
- `structure`: Structure类型实例
- `id`: 单元id
- `degree`: 旋转角度（ᵒ）
"""
function set_tria_rotation!(structure::Structure,id,degree)::Int
    id=string(id)
    if !(id in keys(structure.trias))
        throw("tria id "*string(id)*" doesn't exists!")
    end
    rad=π/180*degree
    Rₓ=[cos(rad) -sin(rad) 0;
        sin(rad) cos(rad) 0;
        0 0 1]
    for i in 1:3
        structure.quads[id].T[6i-5:6i-3,6i-5:6i-3]*=Rₓ
        structure.quads[id].T[6i-2:6i,6i-2:6i]*=Rₓ
    end
    return 0
end

"""
    set_damp_constant!(structure,ζ)
设置structure实例阻尼矩阵为常数阻尼
# Arguments
- `structure`: Structure类型实例
- `ζ`: 阻尼比
"""
function set_damp_constant!(structure::Structure,ζ)::Int
    structure.damp="constant"
    structure.ζ₁=ζ
    return 0
end

"""
    set_damp_Rayleigh!(structure,α,β)
设置structure实例的阻尼矩阵为Rayleigh阻尼
# Arguments
- `structure`: Structure类型实例
- `α`: 刚度矩阵系数
- `β`: 质量矩阵系数
"""
function set_damp_Rayleigh!(structure::Structure,α,β)::Int
    structure.damp="rayleigh"
    structure.ζ₁=α
    structure.ζ₂=β
    return 0
end


text:
convert(::Type{EventTraceData}, C::EventChannel) = EventTraceData(C)

function convert(::Type{EventTraceData}, S::T) where T <:  GphysData
  (T == EventTraceData) && (return deepcopy(S))
  if T != SeisData
    S = convert(SeisData, S)
  end
  TD = EventTraceData(getfield(S, :n))
  for f in datafields
    if (f in unindexed_fields) == false
      setfield!(TD, f, deepcopy(getfield(S, f)))
    end
  end
  return TD
end
EventTraceData(S::T) where T <: GphysData = convert(EventTraceData, S)

function convert(::Type{SeisData}, TD::EventTraceData)
  S = SeisData(getfield(TD, :n))
  for f in datafields
    if (f in unindexed_fields) == false
      setfield!(S, f, deepcopy(getfield(TD, f)))
    end
  end
  return S
end

function convert(::Type{SeisChannel}, D::EventChannel)
  C = SeisChannel()
  for f in datafields
    setfield!(C, f, deepcopy(getfield(D, f)))
  end
  return C
end

function convert(::Type{EventChannel}, C::T) where T <: GphysChannel
  (T == EventChannel) && (return deepcopy(C))
  if T != SeisChannel
    S = convert(SeisChannel, S)
  end
  D = EventChannel()
  for f in datafields
    setfield!(D, f, deepcopy(getfield(C, f)))
  end
  return D
end
EventChannel(C::T) where T <: GphysChannel = convert(EventChannel, C)

function unsafe_convert(::Type{SeisData}, TD::EventTraceData)
  S = SeisData(getfield(TD, :n))
  for f in datafields
    if (f in unindexed_fields) == false
      setfield!(S, f, getfield(TD, f))
    end
  end
  return S
end

function unsafe_convert(::Type{EventTraceData}, S::SeisData)
  TD = EventTraceData(getfield(S, :n))
  for f in datafields
    if (f in unindexed_fields) == false
      setfield!(TD, f, getfield(S, f))
    end
  end
  return TD
end

push!(TD::EventTraceData, C::SeisChannel) = push!(TD, convert(EventChannel, C))


text:
# Template to define new Loss Functions

############################################################################
### Type
############################################################################

struct SemRidge{FT, GT, HT} <: SemLossFunction
    α

    F::FT
    G::GT
    H::HT
end

############################################################################
### Constructor
############################################################################

function SemRidge(α, n_par; parameter_type = Float64)

    return SemRidge(
        α,

        zeros(parameter_type, 1),
        zeros(parameter_type, n_par),
        zeros(parameter_type, n_par, n_par))
end

############################################################################
### functors
############################################################################

function (ridge::SemRidge)(parameters, F, G, H, model)

    if F
        F = ridge.α*norm(parameters)^2
        mytype.F[1] = F
    end

    if G
        G = 2*ridge.α*parameters
        mytype.G .= G
    end

    if H
        stop("hessian for MyType is not implemented")
    end

end


text:
@recipe(Violin, x, y) do scene
    Theme(;
        default_theme(scene, Poly)...,
        side = :both,
        width = 0.8,
        mediancolor = :white,
        show_median = false
    )
end

conversion_trait(x::Type{<:Violin}) = SampleBased()

function plot!(plot::Violin)
    width, side = plot[:width], plot[:side]

    signals = lift(plot[1], plot[2], width, side) do x, y, bw, vside
        meshes = GeometryTypes.GLPlainMesh[]
        lines = Pair{Point2f0, Point2f0}[]
        for (key, idxs) in finduniquesorted(x)
            v = view(y, idxs)
            spec = (x = key, kde = kde(v), median = median(v))
            min, max = extrema_nan(spec.kde.density)
            xl = reverse(spec.x .- spec.kde.density .* (0.5*bw/max))
            xr = spec.x .+ spec.kde.density .* (0.5*bw/max)
            yl = reverse(spec.kde.x)
            yr = spec.kde.x

            x_coord = vside == :left ? xl : vside == :right ? xr : vcat(xr, xl)
            y_coord = vside == :left ? yl : vside == :right ? yr : vcat(yr, yl)
            mesh = GeometryTypes.GLPlainMesh(Point2f0.(x_coord, y_coord))
            push!(meshes, mesh)
            median_left = Point2f0(vside == :right ? spec.x : spec.x-(0.5*bw), spec.median)
            median_right = Point2f0(vside == :left ? spec.x : spec.x+(0.5*bw), spec.median)
            push!(lines, median_left => median_right)
        end
        return meshes, lines
    end
    mesh!(plot, lift(first, signals), color = plot[:color], visible = plot[:visible])
    linesegments!(plot, lift(last, signals), color = plot[:mediancolor], visible = plot[:show_median])
end


text:
module Profiler

# profiledata.jl
export  FunctionData,
        XNode,
        XTensor,
        enums, times, gettime, bytes,
        can_select_algo,
        is_persistent,
        isarg,
        nodes, tensors,
        unx,
        producer,
        consumer,
        users,
        locations,
        get_configs,
        getconfig,
        _setup!,
        configs_for,
        live_tensors

# cache.jl
export CPUKernelCache, GPUKernelCache

# kernels.jl
export profile, read_timing_data

# compare.jl
export compare

# stdlib deps
using Serialization, Dates

using ..Utils
using ProgressMeter
using JuMP
using PersistentArrays
using Statistics
using DocStringExtensions
import nGraph
import ..AutoTM._env_context

include("functiondata.jl")
include("cache.jl")
include("kernels.jl")
include("compare.jl")
include("inspect.jl")

end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
const FaceMesh{Dim,T,Element} = Mesh{Dim,T,Element,<:FaceView{Element}}

coordinates(mesh::FaceMesh) = coordinates(getfield(mesh, :simplices))
faces(mesh::FaceMesh) = faces(getfield(mesh, :simplices))

function texturecoordinates(mesh::AbstractMesh)
    hasproperty(mesh, :uv) && return mesh.uv
    hasproperty(mesh, :uvw) && return mesh.uvw
    return nothing
end

function normals(mesh::AbstractMesh)
    hasproperty(mesh, :normals) && return mesh.normals
    return nothing
end

const GLTriangleElement = Triangle{3,Float32}
const GLTriangleFace = TriangleFace{GLIndex}
const PointWithUV{Dim,T} = PointMeta{Dim,T,Point{Dim,T},(:uv,),Tuple{Vec{2,T}}}
const PointWithNormal{Dim,T} = PointMeta{Dim,T,Point{Dim,T},(:normals,),Tuple{Vec{3,T}}}
const PointWithUVNormal{Dim,T} = PointMeta{Dim,T,Point{Dim,T},(:normals, :uv),
                                           Tuple{Vec{3,T},Vec{2,T}}}
const PointWithUVWNormal{Dim,T} = PointMeta{Dim,T,Point{Dim,T},(:normals, :uvw),
                                            Tuple{Vec{3,T},Vec{3,T}}}

"""
    TriangleMesh{Dim, T, PointType}

Abstract Mesh with triangle elements of eltype `T`.
"""
const TriangleMesh{Dim,T,PointType} = AbstractMesh{TriangleP{Dim,T,PointType}}

"""
    PlainMesh{Dim, T}

Triangle mesh with no meta information (just points + triangle faces)
"""
const PlainMesh{Dim,T} = TriangleMesh{Dim,T,Point{Dim,T}}
const GLPlainMesh{Dim} = PlainMesh{Dim,Float32}
const GLPlainMesh2D = GLPlainMesh{2}
const GLPlainMesh3D = GLPlainMesh{3}

"""
    UVMesh{Dim, T}

PlainMesh with texture coordinates meta at each point.
`uvmesh.uv isa AbstractVector{Vec2f}`
"""
const UVMesh{Dim,T} = TriangleMesh{Dim,T,PointWithUV{Dim,T}}
const GLUVMesh{Dim} = UVMesh{Dim,Float32}
const GLUVMesh2D = UVMesh{2}
const GLUVMesh3D = UVMesh{3}

"""
    NormalMesh{Dim, T}

PlainMesh with normals meta at each point.
`normalmesh.normals isa AbstractVector{Vec3f}`
"""
const NormalMesh{Dim,T} = TriangleMesh{Dim,T,PointWithNormal{Dim,T}}
const GLNormalMesh{Dim} = NormalMesh{Dim,Float32}
const GLNormalMesh2D = GLNormalMesh{2}
const GLNormalMesh3D = GLNormalMesh{3}

"""
    NormalUVMesh{Dim, T}

PlainMesh with normals and uv meta at each point.
`normalmesh.normals isa AbstractVector{Vec3f}`
`normalmesh.uv isa AbstractVector{Vec2f}`
"""
const NormalUVMesh{Dim,T} = TriangleMesh{Dim,T,PointWithUVNormal{Dim,T}}
const GLNormalUVMesh{Dim} = NormalUVMesh{Dim,Float32}
const GLNormalUVMesh2D = GLNormalUVMesh{2}
const GLNormalUVMesh3D = GLNormalUVMesh{3}

"""
    NormalUVWMesh{Dim, T}

PlainMesh with normals and uvw (texture coordinates in 3D) meta at each point.
`normalmesh.normals isa AbstractVector{Vec3f}`
`normalmesh.uvw isa AbstractVector{Vec3f}`
"""
const NormalUVWMesh{Dim,T} = TriangleMesh{Dim,T,PointWithUVWNormal{Dim,T}}
const GLNormalUVWMesh{Dim} = NormalUVWMesh{Dim,Float32}
const GLNormalUVWMesh2D = GLNormalUVWMesh{2}
const GLNormalUVWMesh3D = GLNormalUVWMesh{3}

function decompose_triangulate_fallback(primitive::Meshable; pointtype, facetype)
    positions = decompose(pointtype, primitive)
    faces = decompose(facetype, primitive)
    # If faces returns nothing for primitive, we try to triangulate!
    if faces === nothing
        # triangulation.jl
        faces = decompose(facetype, positions)
    end
    return positions, faces
end

"""
    mesh(primitive::GeometryPrimitive;
         pointtype=Point, facetype=GLTriangle,
         uv=nothing, normaltype=nothing)

Creates a mesh from `primitive`.

Uses the element types from the keyword arguments to create the attributes.
The attributes that have their type set to nothing are not added to the mesh.
Note, that this can be an `Int` or `Tuple{Int, Int}``, when the primitive is grid based.
It also only losely correlates to the number of vertices, depending on the algorithm used.
#TODO: find a better number here!
"""
function mesh(primitive::Meshable; pointtype=Point, facetype=GLTriangleFace, uv=nothing,
              normaltype=nothing)

    positions, faces = decompose_triangulate_fallback(primitive; pointtype=pointtype, facetype=facetype)

    # We want to preserve any existing attributes!
    attrs = attributes(primitive)
    # Make sure this doesn't contain position, we'll add position explicitely via meta!
    delete!(attrs, :position)

    if uv !== nothing
        # this may overwrite an existing :uv, but will only create a copy
        # if it has a different eltype, otherwise it should replace it
        # with exactly the same instance - which is what we want here
        attrs[:uv] = decompose(UV(uv), primitive)
    end

    if normaltype !== nothing
        primitive_normals = normals(primitive)
        if primitive_normals !== nothing
            attrs[:normals] = decompose(normaltype, primitive_normals)
        else
            # Normals not implemented for primitive, so we calculate them!
            n = normals(positions, faces; normaltype=normaltype)
            if n !== nothing # ok jeez, this is a 2d mesh which cant have normals
                attrs[:normals] = n
            end
        end
    end
    return Mesh(meta(positions; attrs...), faces)
end

"""
    mesh(polygon::AbstractVector{P}; pointtype=P, facetype=GLTriangleFace,
         normaltype=nothing)

Create a mesh from a polygon given as a vector of points, using triangulation.
"""
function mesh(polygon::AbstractVector{P}; pointtype=P, facetype=GLTriangleFace,
              normaltype=nothing) where {P<:AbstractPoint{2}}

    return mesh(Polygon(polygon); pointtype=pointtype, facetype=facetype,
                normaltype=normaltype)
end

function mesh(polygon::AbstractPolygon{Dim,T}; pointtype=Point{Dim,T},
              facetype=GLTriangleFace, normaltype=nothing) where {Dim,T}

    positions, faces = decompose_triangulate_fallback(polygon; pointtype=pointtype, facetype=facetype)

    if normaltype !== nothing
        n = normals(positions, faces; normaltype=normaltype)
        positions = meta(positions; normals=n)
    end
    return Mesh(positions, faces)
end

function triangle_mesh(primitive::Meshable{N}; nvertices=nothing) where {N}
    if nvertices !== nothing
        @warn("nvertices argument deprecated. Wrap primitive in `Tesselation(primitive, nvertices)`")
        primitive = Tesselation(primitive, nvertices)
    end
    return mesh(primitive; pointtype=Point{N,Float32}, facetype=GLTriangleFace)
end

function uv_mesh(primitive::Meshable{N,T}) where {N,T}
    return mesh(primitive; pointtype=Point{N,Float32}, uv=Vec2f, facetype=GLTriangleFace)
end

function uv_normal_mesh(primitive::Meshable{N}) where {N}
    return mesh(primitive; pointtype=Point{N,Float32}, uv=Vec2f, normaltype=Vec3f,
                facetype=GLTriangleFace)
end

function normal_mesh(points::AbstractVector{<:AbstractPoint},
                     faces::AbstractVector{<:AbstractFace})
    _points = decompose(Point3f, points)
    _faces = decompose(GLTriangleFace, faces)
    return Mesh(meta(_points; normals=normals(_points, _faces)), _faces)
end

function normal_mesh(primitive::Meshable{N}; nvertices=nothing) where {N}
    if nvertices !== nothing
        @warn("nvertices argument deprecated. Wrap primitive in `Tesselation(primitive, nvertices)`")
        primitive = Tesselation(primitive, nvertices)
    end
    return mesh(primitive; pointtype=Point{N,Float32}, normaltype=Vec3f,
                facetype=GLTriangleFace)
end

"""
    volume(triangle)

Calculate the signed volume of one tetrahedron. Be sure the orientation of your
surface is right.
"""
function volume(triangle::Triangle) where {VT,FT}
    v1, v2, v3 = triangle
    sig = sign(orthogonal_vector(v1, v2, v3) ⋅ v1)
    return sig * abs(v1 ⋅ (v2 × v3)) / 6
end

"""
    volume(mesh)

Calculate the signed volume of all tetrahedra. Be sure the orientation of your
surface is right.
"""
function volume(mesh::Mesh) where {VT,FT}
    return sum(volume, mesh)
end

function Base.merge(meshes::AbstractVector{<:Mesh})
    return if isempty(meshes)
        error("No meshes to merge")
    elseif length(meshes) == 1
        return meshes[1]
    else
        m1 = meshes[1]
        ps = copy(coordinates(m1))
        fs = copy(faces(m1))
        for mesh in Iterators.drop(meshes, 1)
            append!(fs, map(f -> f .+ length(ps), faces(mesh)))
            append!(ps, coordinates(mesh))
        end
        return Mesh(ps, fs)
    end
end

"""
    pointmeta(mesh::Mesh; meta_data...)

Attaches metadata to the coordinates of a mesh
"""
function pointmeta(mesh::Mesh; meta_data...)
    points = coordinates(mesh)
    attr = attributes(points)
    delete!(attr, :position) # position == metafree(points)
    # delete overlapping attributes so we can replace with `meta_data`
    foreach(k -> delete!(attr, k), keys(meta_data))
    return Mesh(meta(metafree(points); attr..., meta_data...), faces(mesh))
end

function pointmeta(mesh::Mesh, uv::UV)
    return pointmeta(mesh; uv=decompose(uv, mesh))
end

function pointmeta(mesh::Mesh, normal::Normal)
    return pointmeta(mesh; normals=decompose(normal, mesh))
end

"""
    pop_pointmeta(mesh::Mesh, property::Symbol)
Remove `property` from point metadata.
Returns the new mesh, and the property!
"""
function pop_pointmeta(mesh::Mesh, property::Symbol)
    points = coordinates(mesh)
    attr = attributes(points)
    delete!(attr, :position) # position == metafree(points)
    # delete overlapping attributes so we can replace with `meta_data`
    m = pop!(attr, property)
    return Mesh(meta(metafree(points); attr...), faces(mesh)), m
end

"""
    facemeta(mesh::Mesh; meta_data...)

Attaches metadata to the faces of a mesh
"""
function facemeta(mesh::Mesh; meta_data...)
    return Mesh(coordinates(mesh), meta(faces(mesh); meta_data...))
end

function attributes(hasmeta::Mesh)
    return Dict{Symbol,Any}((name => getproperty(hasmeta, name)
                             for name in propertynames(hasmeta)))
end