text:
@testset "Commands" begin
    delete_handler!(c, MessageCreate)

    # Adding handlers adds to MessageCreate.
    add_command!(c, :foo, (c, m) -> nothing)
    @test haskey(c.handlers[MessageCreate], :foo)

    # We can delete commands with a shortcut.
    delete_command!(c, :foo)
    @test !haskey(c.handlers[MessageCreate], :foo)

    # TODO: Finish these tests.
end


text:
using Printf

function foo()
    # TODO: Support POST
    # TODO: Foobar print
    printf("foobar")
    # TODO: End function
end

ext: jl

max_stars_repo_name: StephaneRob/leasot

max_issues_repo_name: StephaneRob/leasot

max_forks_repo_name: StephaneRob/leasot


text:
# This file is a part of ForwardDiffPullbacks.jl, licensed under the MIT License (MIT).

function dual_number end
function dual_value end
function dual_partials end
function dual_tagtype end

@inline dual_number(::Type{TagType}, x::Real, p::NTuple{N,Real}) where {TagType,N} = ForwardDiff.Dual{TagType}( x, p...)
@inline dual_value(x::Real) = ForwardDiff.value(x)
@inline dual_partials(x::Real) = ForwardDiff.partials(x)
@inline dual_tagtype(f::Any, ::Type{T}) where T = typeof(ForwardDiff.Tag(f, T))


text:
#
# Functions pertaining to instantiated models and their components
#

"""
    modeldef(mi)

Return the `ModelDef` contained by ModelInstance `mi`.
"""
modeldef(mi::ModelInstance) = mi.md

compmodule(obj::AbstractComponentInstance) = compmodule(obj.comp_id)
compname(obj::AbstractComponentInstance)   = compname(obj.comp_id)

"""
    add_comp!(obj::AbstractCompositeComponentInstance, ci::AbstractComponentInstance)

Add the (leaf or composite) component `ci` to a composite's list of components.
"""
function add_comp!(obj::AbstractCompositeComponentInstance, ci::AbstractComponentInstance)
    obj.comps_dict[nameof(ci)] = ci
    return nothing
end

#
# Setting/getting parameter and variable values
#

# Get the object stored for the given variable, not the value of the variable.
# This is used in the model building process to connect internal parameters.
@inline function get_property_obj(obj::ComponentInstanceParameters{NT}, name::Symbol) where {NT}
    return getproperty(nt(obj), name)
end

@inline function get_property_obj(obj::ComponentInstanceVariables{NT}, name::Symbol) where {NT}
    return getproperty(nt(obj), name)
end

@inline function _get_prop(nt::NT, name::Symbol) where {NT <: NamedTuple}
    obj = getproperty(nt, name)
    return fieldtype(NT, name) <: ScalarModelParameter ? obj.value : obj
end

@inline function Base.getproperty(obj::ComponentInstanceParameters{NT}, name::Symbol) where {NT}
    return _get_prop(nt(obj), name)
end

@inline function Base.getproperty(obj::ComponentInstanceVariables{NT}, name::Symbol) where {NT}
    return _get_prop(nt(obj), name)
end

@inline function Base.setproperty!(obj::ComponentInstanceParameters{NT}, name::Symbol, value::VTYPE) where {NT, VTYPE}
    prop_obj = get_property_obj(obj, name)
    T = fieldtype(NT, name)

    if T <: ScalarModelParameter
        return setproperty!(prop_obj, :value, value)
    else
        error("You cannot override indexed parameter $name::$T.")
    end
end

@inline function Base.setproperty!(obj::ComponentInstanceVariables{NT}, name::Symbol, value::VTYPE) where {NT, VTYPE}
    prop_obj = get_property_obj(obj, name)
    T = fieldtype(NT, name)

    if T <: ScalarModelParameter
        return setproperty!(prop_obj, :value, value)
    else
        error("You cannot override indexed variable $name::$T. Make sure you are using proper indexing syntax in the `run_timestep` function: v.varname[t] = ...")
    end
end

comp_paths(obj::AbstractComponentInstanceData) = getfield(obj, :comp_paths)

"""
    get_param_value(ci::AbstractComponentInstance, name::Symbol)

Return the value of parameter `name` in (leaf or composite) component `ci`.
"""
function get_param_value(ci::AbstractComponentInstance, name::Symbol)
    try
        return getproperty(ci.parameters, name)
    catch err
        if isa(err, KeyError)
            error("Component $(ci.comp_id) has no parameter named $name")
        else
            rethrow(err)
        end
    end
end

"""
    get_var_value(ci::AbstractComponentInstance, name::Symbol)

Return the value of variable `name` in component `ci`.
"""
function get_var_value(ci::AbstractComponentInstance, name::Symbol)
    try
        vars = ci.variables
        # @info ("Getting $name from $vars")
        return getproperty(vars, name)
    catch err
        if isa(err, KeyError)
            error("Component $(ci.comp_id) has no variable named $name")
        else
            rethrow(err)
        end
    end
end

"""
    set_param_value(ci::AbstractComponentInstance, name::Symbol, value)

Set the value of parameter `name` in component `ci` to `value`.
"""
set_param_value(ci::AbstractComponentInstance, name::Symbol, value) = setproperty!(ci.parameters, name, value)

"""
    set_var_value(ci::AbstractComponentInstance, name::Symbol, value)

Set the value of variable `name` in component `ci` to `value`.
"""
set_var_value(ci::AbstractComponentInstance, name::Symbol, value) = setproperty!(ci.variables, name, value)

"""
    variables(obj::AbstractCompositeComponentInstance, comp_name::Symbol)

Return the `ComponentInstanceVariables` for `comp_name` in CompositeComponentInstance `obj`.
"""
variables(obj::AbstractCompositeComponentInstance, comp_name::Symbol) = variables(compinstance(obj, comp_name))

variables(obj::AbstractComponentInstance) = obj.variables

function variables(m::Model)
    if ! is_built(m)
        error("Must build model to access variable instances. Use variables(modeldef(m)) to get variable definitions.")
    end
    return variables(modelinstance(m))
end

"""
    parameters(obj::AbstractComponentInstance, comp_name::Symbol)

Return the `ComponentInstanceParameters` for `comp_name` in CompositeComponentInstance `obj`.
"""
parameters(obj::AbstractCompositeComponentInstance, comp_name::Symbol) = parameters(compinstance(obj, comp_name))

parameters(obj::AbstractComponentInstance) = obj.parameters

function Base.getindex(mi::ModelInstance, names::NTuple{N, Symbol}) where N
    obj = mi

    # skip past first element if same as root node
    if length(names) > 0 && head(obj.comp_path) == names[1]
        names = names[2:end]
    end

    for name in names
        if has_comp(obj, name)
            obj = obj[name]
        else
            error("Component $(obj.comp_path) does not have sub-component :$name")
        end
    end
    return obj
end

Base.getindex(mi::ModelInstance, comp_path::ComponentPath) = getindex(mi, comp_path.names)

Base.getindex(mi::ModelInstance, path_str::AbstractString) = getindex(mi, ComponentPath(mi.md, path_str))

function Base.getindex(obj::AbstractCompositeComponentInstance, comp_name::Symbol)
    if ! has_comp(obj, comp_name)
        error("Component :$comp_name does not exist in the given composite")
    end
    return compinstance(obj, comp_name)
end

# TBD we could probably combine the two _get_datum methods into one that takes a
# ci::AbstractComponentInstance, but for now it seems there are enough differences
# that keeping them separate is cleaner
function _get_datum(ci::CompositeComponentInstance, datum_name::Symbol)
    vars = variables(ci)

    if datum_name in keys(vars) # could merge with method below if made names(NamedTuple) = keys(NamedTuple)
        which = vars
    else
        pars = parameters(ci)
        if datum_name in keys(pars) # could merge with method below if made names(NamedTuple) = keys(NamedTuple)
            which = pars
        else
            error("$datum_name is not a parameter or a variable in component $(ci.comp_path).")
        end
    end

    ref = getproperty(which, datum_name)

    return _get_datum(ci.comps_dict[ref.comp_name], ref.datum_name)
end

function _get_datum(ci::LeafComponentInstance, datum_name::Symbol)
    vars = variables(ci)

    if datum_name in names(vars)
        which = vars
    else
        pars = parameters(ci)
        if datum_name in names(pars)
            which = pars
        else
            error("$datum_name is not a parameter or a variable in component $(ci.comp_path).")
        end
    end

    value = getproperty(which, datum_name)

    if value isa TimestepArray
        return value.data isa SubArray ? parent(value.data) : value.data
    else
        return value
    end
end

function Base.getindex(mi::ModelInstance, key::AbstractString, datum::Symbol)
    _get_datum(mi[key], datum)
end

function Base.getindex(mi::ModelInstance, comp_path::ComponentPath, datum::Symbol)
    _get_datum(mi[comp_path], datum)
end

@delegate Base.getindex(m::Model, comp_path::ComponentPath, datum::Symbol) => mi

function Base.getindex(obj::AbstractCompositeComponentInstance, comp_name::Symbol, datum::Symbol)
    ci = obj[comp_name]
    return _get_datum(ci, datum)
end

"""
    dim_count(mi::ModelInstance, dim_name::Symbol)

Return the size of index `dim_name` in model instance `mi`.
"""
@delegate dim_count(mi::ModelInstance, dim_name::Symbol) => md

function reset_variables(ci::AbstractComponentInstance)
    # @info "reset_variables($(ci.comp_id))"
    vars = ci.variables

    for (name, T) in zip(names(vars), types(vars))
        value = getproperty(vars, name)

        if (T <: AbstractArray || T <: TimestepArray) && eltype(value) <: AbstractFloat
            fill!(value, NaN)

        elseif T <: AbstractFloat || (T <: ScalarModelParameter && T.parameters[1] <: AbstractFloat)
            setproperty!(vars, name, NaN)

        elseif (T <: ScalarModelParameter)    # integer or bool
            setproperty!(vars, name, 0)
        end
    end
end

function reset_variables(obj::AbstractCompositeComponentInstance)
    for ci in components(obj)
        reset_variables(ci)
    end
    return nothing
end

function init(ci::AbstractComponentInstance, dims::DimValueDict)
    # @info "init($(ci.comp_id))"
    reset_variables(ci)

    if ci.init != nothing
        ci.init(parameters(ci), variables(ci), dims)
    end
    return nothing
end

function init(obj::AbstractCompositeComponentInstance, dims::DimValueDict)
    for ci in components(obj)
        init(ci, dims)
    end
    return nothing
end

_runnable(ci::AbstractComponentInstance, clock::Clock) = (ci.first <= gettime(clock) <= ci.last)

function get_shifted_ts(ci, ts::FixedTimestep{FIRST, STEP, LAST}) where {FIRST, STEP, LAST}
    if ci.first == FIRST && ci.last == LAST
        return ts
    else
        # shift the timestep over by (ci.first - FIRST)/STEP
        return FixedTimestep{ci.first,STEP,ci.last}(ts.t - Int((ci.first - FIRST)/STEP))
    end
end

function get_shifted_ts(ci, ts::VariableTimestep{TIMES}) where {TIMES}
    if ci.first == TIMES[1] && ci.last == TIMES[end]
        return ts
    else
        # shift the timestep over by the number of timesteps between the model first and the ts first
        idx_start = findfirst(TIMES .== ci.first)
        idx_finish = findfirst(TIMES .== ci.last)
        return VariableTimestep{TIMES[idx_start:idx_finish]}(ts.t - idx_start + 1)
    end
end

function run_timestep(ci::AbstractComponentInstance, clock::Clock, dims::DimValueDict)
    if ci.run_timestep !== nothing && _runnable(ci, clock)
        ci.run_timestep(parameters(ci), variables(ci), dims, get_shifted_ts(ci, clock.ts))
    end

    return nothing
end

function run_timestep(cci::AbstractCompositeComponentInstance, clock::Clock, dims::DimValueDict)
    if _runnable(cci, clock)
        for ci in components(cci)
            run_timestep(ci, clock, dims)
        end
    end
    return nothing
end

"""
    Base.run(mi::ModelInstance, ntimesteps::Int=typemax(Int),
            dimkeys::Union{Nothing, Dict{Symbol, Vector{T} where T <: DimensionKeyTypes}}=nothing)

Run the `ModelInstance` `mi` once with `ntimesteps` and dimension keys `dimkeys`.
"""
function Base.run(mi::ModelInstance, ntimesteps::Int=typemax(Int),
                  dimkeys::Union{Nothing, Dict{Symbol, Vector{T} where T <: DimensionKeyTypes}}=nothing)

    if length(components(mi)) == 0
        error("Cannot run the model: no components have been created.")
    end

    time_keys::Vector{Int} = dimkeys === nothing ? dim_keys(mi.md, :time) : dimkeys[:time]

    # truncate time_keys if caller so desires
    if ntimesteps < length(time_keys)
        time_keys = time_keys[1:ntimesteps]
    end

    clock = Clock(time_keys)

    # Get the dimensions dictionary
    dim_val_dict = DimValueDict(dim_dict(mi.md), clock)

    # recursively initializes all components
    init(mi, dim_val_dict)

    while ! finished(clock)
        run_timestep(mi, clock, dim_val_dict)
        advance(clock)
    end

    nothing
end


text:
abstract type ConductanceMethod end
struct Thom1972 <: ConductanceMethod end
struct Choudhury1988 <: ConductanceMethod end
struct Su2001 <: ConductanceMethod end
struct ConstantKB1 <: ConductanceMethod end

"""
    compute_Gb!(df::AbstractDataFrame, approach; kwargs...)

Estimate boundary layer conductance.

# Arguments
- `df`       : DataFrame with required variables (depend on approach)
- `approach` : one of
  - `Thom1972()`: see [`Gb_Thom`](@ref)
  - `Choudhury1988()`: see [`Gb_Choudhury`](@ref)
  - `Su2001()`: see [`Gb_Su`](@ref)
  - `ConstantKB1()`: see [`Gb_constant_kB1`](@ref)

The different approaches require different variables to be present in `df` and
different keyword arguments.

# Value
updated DataFrame `df` with the following columns:
- `Gb_h`: Boundary layer conductance for heat transfer (m s-1)

To subsequently derived quantities see
- [`compute_Gb_quantities`](@ref) for Resistance, kB^(-1) constant, and CO2 conductance
- [`add_Gb!`](@ref) for conductances of other species given their Schmidt numbers.

# See also
[`Gb_Thom`](@ref), [`Gb_Choudhury`](@ref), [`Gb_Su`](@ref), [`Gb_constant_kB1`](@ref),
[`aerodynamic_conductance!`](@ref)

```jldoctest; output = false
using DataFrames
df = DataFrame(wind=[3,4,5], ustar=[0.5,0.6,0.65])
compute_Gb!(df, Thom1972())
≈(df.Gb_h[1], 0.102, rtol=1e-2)
# output
true
```
"""
function compute_Gb!(df::AbstractDataFrame, approach::Thom1972; kwargs...)
  fr = (ustar) -> Gb_Thom.(ustar; kwargs...)
  #transform!(df, :ustar => ByRow(fr) => :Gb_h) # does not work with SVector
  transform!(df, :ustar => fr => :Gb_h)
end
function compute_Gb!(df::AbstractDataFrame, approach::ConstantKB1; kB_h, kwargs...)
  # do not use ByRow because kb_H can be a vector
  ft = (ustar) -> Gb_constant_kB1.(ustar, kB_h; kwargs...)
  transform!(df, :ustar => ft => :Gb_h)
end
function compute_Gb!(df::AbstractDataFrame, approach::Choudhury1988; kwargs...)
  Gb_Choudhury!(df; kwargs...)
end
function compute_Gb!(df::AbstractDataFrame, approach::Su2001; kwargs...)
  Gb_Su!(df; kwargs...)
end

"""
    compute_Gb_quantities(Gb_h)
    compute_Gb_quantities!(df:::AbstractDataFrame)

Based on boundary layer conductance for heat, compute derived quantities.

# Arguments
- `Gb_h` : Boundary layer conductance for heat transfer (m s-1)
- `df`   : DataFrame with above columns
- `constants=`[`BigleafConstants`](@ref)`()`: entries `Sc_CO2` and `Pr`

# Value
NamedTuple with entries
- `Gb_h`: Boundary layer conductance for heat transfer (m s-1)
- `Rb_h`: Boundary layer resistance for heat transfer (s m-1)
- `kB_h`: kB^(-1) parameter for heat transfer
- `Gb_CO2`: Boundary layer conductance for CO2 (m s-1).
"""
function compute_Gb_quantities(Gb_h, ustar; constants=BigleafConstants())
  (ismissing(Gb_h) || ismissing(ustar)) && return(
    (Rb_h=missing, Gb_h=missing, kB_h=missing, Gb_CO2=missing))
  Rb_h = 1/Gb_h
  kB_h = Rb_h*oftype(ustar,constants.k)*ustar
  #Gb_CO2 = Gb_h / ((Gb_h,constants.Sc_CO2/constants.Pr)^0.67)
  Gb_CO2 = Gb_h /
    ((oftype(Gb_h,constants.Sc_CO2)/oftype(Gb_h,constants.Pr))^oftype(Gb_h,0.67))
  (;Rb_h, Gb_h, kB_h, Gb_CO2)
end
function compute_Gb_quantities!(df::AbstractDataFrame; constants=BigleafConstants())
  ft = (args...) -> compute_Gb_quantities.(args...; constants)
  transform!(df, SA[:Gb_h, :ustar] => ft => AsTable)
end

"""
    add_Gb(Gb_h::Union{Missing,Number}, Sc::Vararg{Pair,N}; constants)
    add_Gb!(df::AbstractDataFrame, Sc::Vararg{Pair,N}; Gb_h = df.Gb_h, kwargs...)

compute boundary layer conductance for additional quantities for given Schmidt-numbers

# Arguments
- `Gb_h`     : Boundary layer conductance for heat transfer (m s-1)
- `Sc`       : several `Pair{Symbol,Number}` Output name and Schmidt number of
               additional conductances to be calculated
- `df`       : DataFrame to add output columns
optional
- `constants=`[`BigleafConstants`](@ref)`()`: Dictionary with entries
  - `Pr` - Prandtl number

# Details
Boundary layer conductance for other quantities x is calculated
based on boundary layer for heat transfer as (Hicks et al. 1987):

  ``Gb_x = Gb_h / (Sc_x / Pr)^{0.67}``

where `Sc_x` is the Schmidt number of quantity x, and Pr is the Prandtl number (0.71).

# Value
a NameTuple or `df` with keys `Gb_x` where `x` are the keys in `Sc` and
corresponding boundary layer conductances (m s-1).

# Examples
```jldoctest; output=false
using DataFrames
df = DataFrame(Gb_h=[0.02, missing, 0.055])
add_Gb!(df, :O2 => 0.84, :CH4 => 0.99)
propertynames(df)[2:3] == [:Gb_O2, :Gb_CH4]
# output
true
```
"""
function add_Gb!(df::AbstractDataFrame, Sc::Vararg{Pair,N}; Gb_h = df.Gb_h, kwargs...) where N
  N == 0 && return(df)
  Scn, Scv = get_names_and_values("Gb_", Sc...)
  ft = () -> add_Gb_.(Gb_h, Ref(Scn), Ref(Scv); kwargs...)
  transform!(df, [] => ft => AsTable)
end
function add_Gb(Gb_h::Union{Missing,Number}, Sc::Vararg{Pair,N}; kwargs...) where N
  Scn, Scv = get_names_and_values("Gb_", Sc...)
  add_Gb_(Gb_h, Scn, Scv; kwargs...)
end
function add_Gb_(Gb_h::FTM, Scn::NTuple{N,Symbol}, Scv::NTuple{N};
  constants=BigleafConstants()) where {N, FTM<:Union{Missing,Number}}
  #Gbxv = @. Gb_h / ((Scv/constants.Pr)^0.67)
  ismissing(Gb_h) && return(NamedTuple{Scn}(ntuple(_->missing,N)))
  Gbxv = @. Gb_h / ((Scv/FTM(constants.Pr))^FTM(0.67))
  Gbx = NamedTuple{Scn}(Gbxv)
end
function get_names_and_values(prefix::AbstractString, Sc::Vararg{Pair,N}) where N
  Scn = ntuple(i -> Symbol(prefix * string(Sc[i].first)), N)
  Scv = ntuple(i -> Sc[i].second, N)
  Scn, Scv
end

"""
    Gb_Thom(ustar; constants)
    compute_Gb!(df, ::Thom1972)

Boundary Layer Conductance according to Thom 1972, an empirical formulation for the
for heat transfer based on a simple ustar (friction velocity) dependency.

# Arguments
- `ustar`     : Friction velocity (m s-1)
- `df`        : DataFrame with above variables
- `constants=`[`BigleafConstants`](@ref)`()`

# Details
The empirical equation for Rb suggested by Thom 1972 is:

``Rb = 6.2 {u^*}^{-0.67}``

Gb (=1/Rb) for water vapor and heat are assumed to be equal in this package.

# Value
see [`compute_Gb!`](@ref)

# References
- Thom, A., 1972: Momentum, mass and heat exchange of vegetation.
  Quarterly Journal of the Royal Meteorological Society 98, 124-134.
- Hicks, B_B., Baldocchi, D_D., Meyers, T_P., Hosker, J_R., Matt, D_R., 1987:
  A preliminary multiple resistance routine for deriving dry deposition velocities
  from measured quantities. Water, Air, and Soil Pollution 36, 311-330.

```@example; output = false
using DataFrames
df = DataFrame(ustar = SA[0.1,missing,0.3])
compute_Gb!(df, Thom1972())
propertynames(df) == [:ustar, :Gb_h]
```
"""
function Gb_Thom(ustar::Union{Missing,Number}; constants=BigleafConstants())
  Rb_h = 6.2*ustar^-0.667
  Gb_h = 1/Rb_h
end

"""
    Gb_constant_kB1(ustar, kB_h; constants)
    compute_Gb!(df, ::ConstantKB1})

Boundary Layer Conductance using constant kB^(-1) value for heat transfer.

# Arguments
- `ustar`     : Friction velocity (m s-1)
- `df`        : DataFrame with above variables
- `kB_h`      : kB^(-1) value for heat transfer
- `constants=`[`BigleafConstants`](@ref)`()`

# Details
Rb_h computed by ``kB_h/(k * ustar)``, where k is the von Karman constant.
"""
function Gb_constant_kB1(ustar, kB_h; constants=BigleafConstants())
  ismissing(ustar) && return(missing)
  Rb_h = kB_h/(oftype(ustar,constants.k) * ustar)
  Gb_h = 1/Rb_h
  # Gb_CO2 = Gb_h / (constants.Sc_CO2]/constants[:Pr)^0.67
  # (;Rb_h, Gb_h, kB_h, Gb_CO2)
end

"""
    Gb_Choudhury(; leafwidth, LAI, wind_zh, constants)
    Gb_Choudhury!(df; leafwidth, LAI, wind_zh, constants)

Estimate the canopy boundary layer conductance
for heat transfer according to Choudhury & Monteith 1988.

# Arguments
- `df`               : DataFrame where `Gb_h` is to be added/updated
- `leafwidth`        : Leaf width (m)
- `LAI`              : One-sided leaf area index
- `wind_zh`          : Wind speed at canopy heihgt (m s-1), see [`wind_profile`](@ref)
- `constants=`[`BigleafConstants`](@ref)`()`

# Value
see [`compute_Gb!`](@ref)

# Details
Boundary layer conductance according to Choudhury & Monteith 1988 is
given by:

``Gb_h = LAI \\left( 2a/\\alpha \\sqrt{u(z_h)/w} (1-e^{-\\alpha/2})\\right)``

where ``\\alpha`` is modeled as an empirical relation to LAI (McNaughton & van den Hurk 1995):

``\\alpha = 4.39 - 3.97 e^{-0.258 \\, LAI}``

``w`` is leafwidth and ``u(zh)`` is the wind speed at the canopy surface.

It can be approximated from measured wind speed at sensor height zr and a wind extinction
coefficient ``\\alpha``: ``u(z_h) = u(z_r) / e^{\\alpha(z_r/z_h -1)}``.
However, here (if not explicitly given) it is estimated by [`wind_profile`](@ref)

# References
- Choudhury, B. J., Monteith J_L., 1988: A four-layer model for the heat
  budget of homogeneous land surfaces. Q. J. R. Meteorol. Soc. 114, 373-398.
- McNaughton, K. G., Van den Hurk, B_J_J_M., 1995: A 'Lagrangian' revision of
  the resistors in the two-layer model for calculating the energy budget of a
  plant canopy. Boundary-Layer Meteorology 74, 261-288.
- Hicks, B_B., Baldocchi, D_D., Meyers, T_P., Hosker, J_R., Matt, D_R., 1987:
  A preliminary multiple resistance routine for deriving dry deposition velocities
  from measured quantities. Water, Air, and Soil Pollution 36, 311-330.
"""
function Gb_Choudhury(; leafwidth, LAI, wind_zh::FT, constants=BigleafConstants()) where FT
  FT == Missing && return(missing)
  alpha   = FT(4.39) - FT(3.97)*exp(FT(-0.258)*FT(LAI))
  # (ismissing(wind_zh) || isnothing(wind_zh)) && return(
  #   (Rb_h=missing, Gb_h=missing, kB_h=missing, Gb_CO2=missing))
  wind_zh = max(FT(0.01), wind_zh) ## avoid zero windspeed
  Gb_h = FT(LAI)*((FT(0.02)/alpha)*sqrt(wind_zh/FT(leafwidth))*(FT(1)-exp(-FT(alpha)/FT(2))))
end
function Gb_Choudhury!(
  df::AbstractDataFrame; leafwidth, LAI, wind_zh, constants=BigleafConstants()
  )
  # if isnothing(wind_zh)
  #   wind_zh = wind_profile(zh, df, d, z0m; zh, zr, stab_formulation, constants)
  # end
  # Broadcasting does not work over keyword arguments, need to pass as positional
  fwind(wind_zh, leafwidth, LAI; kwargs...) = Gb_Choudhury(;wind_zh, leafwidth, LAI, kwargs...)
  ft = () -> fwind.(wind_zh, leafwidth, LAI; constants)
  transform!(df, [] => ft => :Gb_h)
end

"""
    Gb_Su(Tair,pressure,ustar; wind_zh, Dl, fc, N=2, Cd=0.2, hs=0.01, constants)
    Gb_Su!(df; wind_zh, Dl, fc=nothing, N=2, Cd=0.2, hs=0.01, LAI, constants)

Estimate Boundary Layer Conductance to heat transfer using the physically based
formulation according to Su et al. 2001.

# Arguments
- `Tair`      : Air temperature (degC)
- `pressure`  : Atmospheric pressure (kPa)
- `ustar`     : Friction velocity (m s-1)
- `df`        : DataFrame or matrix containing the above variables
- `Dl`        : Leaf characteristic dimension (m)
- `fc`        : Fractional vegetation cover (0-1), if not provided, calculated from LAI
- `LAI`       : One-sided leaf area index (-) - alternative to `fc`.
- `N`         : Number of leaf sides participating in heat exchange (defaults to 2)
- `Cd`        : Foliage drag coefficient (-)
- `hs`        : Roughness height of the soil (m)
- `constants=`[`BigleafConstants`](@ref)`()`

# Value
see [`compute_Gb!`](@ref)

# Details
The formulation is based on the kB^(-1) model developed by Massman 1999.
Su et al. 2001 derived the following approximation:

``k_{B1} = (k C_d f_c^2) / (4C_t u^*/u(z_h)) + k_{Bs-1}(1 - f_c)^2``

If ``f_c`` (fractional vegetation cover) is missing, it is estimated from LAI:
``f_c = 1 - e^{-LAI/2}``

The wind speed at the top of the canopy is calculated using function
[`wind_profile`](@ref).

Ct is the heat transfer coefficient of the leaf (Massman 1999):

``C_t = P_r^{-2/3} R_{eh}^{-1/2}  N``

where ``P_r`` is the Prandtl number (set to 0.71),
and ``R_{eh}`` is the Reynolds number for leaves:

``R_{eh} = D_l \\, wind(z_h) / v``

k_{Bs-1}, the k_{B-1} value for bare soil surface, is calculated according
to Su et al. 2001:

``k_{Bs-1} = 2.46(Re)^{0.25} - ln(7.4)``

# References
- Su, Z., Schmugge, T., Kustas, W. & Massman, W., 2001: An evaluation of
  two models for estimation of the roughness height for heat transfer between
  the land surface and the atmosphere. Journal of Applied Meteorology 40, 1933-1951.
-  Massman, W., 1999: A model study of kB H- 1 for vegetated surfaces using
  localized near-field' Lagrangian theory. Journal of Hydrology 223, 27-43.
-  Hicks, B_B., Baldocchi, D_D., Meyers, T_P., Hosker, J_R., Matt, D_R., 1987:
   A preliminary multiple resistance routine for deriving dry deposition velocities
  from measured quantities. Water, Air, and Soil Pollution 36, 311-330.

```jldoctest; output = false
using DataFrames
df = DataFrame(
  Tair=25.0,pressure=100.0,wind=[3.0,4,5],ustar=[0.5,0.6,0.65],H=[200.0,230,250])
zh = 25; zr = 40
z0m = roughness_parameters(
  Roughness_wind_profile(), df.ustar, df.wind, df.Tair, df.pressure, df.H; zh, zr).z0m
wind_zh = wind_profile(zh, df, 0.7*zh, z0m)
compute_Gb!(df,Su2001(); wind_zh, Dl=0.01, LAI=5)
# the same meteorological conditions, but larger leaves
compute_Gb!(df,Su2001(); wind_zh, Dl=0.1,LAI=5)
# same conditions, large leaves, and sparse canopy cover (LAI = 1.5)
compute_Gb!(df,Su2001(); wind_zh, Dl=0.1,LAI=1.5)
≈(df.Gb_h[1], 0.0638, rtol=1e-3)
# output
true
```
"""
function Gb_Su(Tair,pressure,ustar::FT; wind_zh, Dl, fc, N=2, Cd=0.2, hs=0.01,
  constants=BigleafConstants()
  ) where FT
  FT == Missing && return(missing)
  wind_zh = max(FT(0.01),FT(wind_zh)) ## avoid zero windspeed
  v   = kinematic_viscosity(Tair,pressure;constants)
  Re  = Reynolds_Number(Tair,pressure,ustar,FT(hs); constants)
  kBs = FT(2.46) * (Re)^FT(0.25) - log(FT(7.4))
  Reh = FT(Dl) * wind_zh / v
  #Ct  = 1*(constants.Pr^-0.6667)*Reh^-0.5*N
  Ct  = 1*(FT(constants.Pr)^-FT(0.6667)) * Reh^-FT(0.5) * FT(N)
  #
  kB_h = (FT(constants.k)*FT(Cd))/(4*Ct*ustar/wind_zh)*FT(fc)^2 + kBs*(1 - FT(fc))^2
  Rb_h = kB_h/(FT(constants.k)*ustar)
  Gb_h = 1/Rb_h
end
function Gb_Su!(df::AbstractDataFrame; wind_zh, Dl, fc=nothing,
  N=2, Cd=0.2, hs=0.01, LAI, constants=BigleafConstants()
  )
  if isnothing(fc)
    isnothing(LAI) && error("one of 'fc' or 'LAI' must be provided")
    fc = length(LAI) == 1 ? (1-exp(-LAI/2)) : @. (1-exp(-LAI/2))
  end
  isnothing(Dl) && error(
    "need to provide keyword argument Dl with :Su2001 method")
  inputcols = SA[:Tair,:pressure,:ustar]
  # Broadcasting does not work over keyword arguments, need to pass as positional
  fwind = (wind_zh, Dl, fc, N, Cd, hs, args...; kwargs...) -> Gb_Su(args...;
    wind_zh, Dl, fc=oftype(wind_zh,fc), N, Cd=oftype(wind_zh,Cd), hs=oftype(wind_zh,hs),
    kwargs...)
  ft = (args...) -> fwind.(wind_zh, Dl, fc, N, Cd, hs, args...; constants)
  transform!(df, inputcols => ft => :Gb_h)
end


text:
module Circuits

using Unitful
import Unitful: ElectricalResistance, Capacitance, Inductance, Voltage, Current

using SparseArrays

using UnicodeFun

include("tikz.jl")
using TikzPictures
using .TikZ
const preamble = read(joinpath(dirname(@__FILE__), "preamble.tex"), String)

if VERSION < v"1.1-DEV"
    isnothing(::Nothing) = true
    isnothing(::Any) = false
end

# * Generic element

abstract type Element end
abstract type NPort{N} <: Element end
const TwoPort = NPort{2}

const Pin = Union{<:Integer,String}

is_subscriptable(::Integer) = true
is_subscriptable(c::Char) = c ∈ "0123456789+-=()aeioruvxβγρφ"
is_subscriptable(s::AbstractString) = all(map(is_subscriptable, collect(s)))
is_subscriptable(::Any) = false

label(e::Element) = is_subscriptable(e.label) ? to_subscript(e.label) : string(e.label)

pins(e::Element) = throw(ArgumentError("Cannot list pins for abstract Element $(e)"))
pins(e::NPort{N}) where N = 1:N

symbolic_pins(e::Element) = throw(ArgumentError("Cannot list symbolic pins for abstract Element $(e)"))

function pin(e::Element, p::String)
    i = findfirst(isequal(p), symbolic_pins(e))
    isnothing(i) && throw(ArgumentError("Pin $i not present for element $(e)"))
    i
end
pin(e::Element, p::Integer) = p

num_pins(e::Element) = length(pins(e))

symbol(::Element) = "Generic element"
value(::Element) = nothing

function Base.show(io::IO, e::Element)
    write(io, "$(symbol(e))$(label(e))")
    v = value(e)
    !isnothing(v) && write(io, "($(v))")
end

# * Physical two-port

mutable struct PhysicalTwoPort{U,S,polar,L} <: TwoPort
    value::U
    label::L
    meta::Dict{Symbol,Any}
end
PhysicalTwoPort{U,S,polar}(value::U; label::L = 1, meta...) where {U,S,polar,L} =
    PhysicalTwoPort{U,S,polar,L}(value, label, meta)

symbol(::PhysicalTwoPort{U,S,polar,L}) where {U,S,polar,L} = S
value(tp::PhysicalTwoPort) = tp.value

Base.copy(t::P) where {P<:PhysicalTwoPort} =
    P(t.value, t.label, copy(t.meta))

# ** Misc physical two-ports

const Resistor{L} = PhysicalTwoPort{ElectricalResistance,:R,false,L}
const Capacitor{L} = PhysicalTwoPort{Capacitance,:C,false,L}
const PolarCapacitor{L} = PhysicalTwoPort{Capacitance,:C⁺,true,L}
const Inductor{L} = PhysicalTwoPort{Inductance,:L,false,L}

const Diode{L} = PhysicalTwoPort{Nothing,:D,true,L}
Diode(;kwargs...) = Diode(nothing; kwargs...)

const Switch{L} = PhysicalTwoPort{Bool,:S,false,L}
Switch(;kwargs...) = Switch(false; kwargs...)
value(s::Switch) = s.value ? "on" : "off"

symbolic_pins(::PhysicalTwoPort{U,S,true,L}) where {U,S,L} = ["+", "-"]

# *** Two-port TikZ symbols
tikz_label(e::PhysicalTwoPort{U,S,polar,L}) where {U,S,polar,L} = "\$$(S)_{$(e.label)}\$"
tikz_info(e::PhysicalTwoPort{U,S,polar,L}) where {U,S,polar,L} = "info sloped={$(tikz_label(e))=$(e.value)}"
tikz_info(e::PhysicalTwoPort{Nothing,S,polar,L}) where {S,polar,L} = "info sloped={$(tikz_label(e))}"
tikz_info(s::Switch{L}) where L = "info sloped={$(tikz_label(s))}"

tikz_symbol(r::Resistor{L}) where L = "resistor"
tikz_symbol(r::Capacitor{L}) where L = "capacitor"
tikz_symbol(c::PolarCapacitor{L}) where L = "capacitor" # No proper symbol in TikZ
tikz_symbol(i::Inductor{L}) where L = "inductor"
tikz_symbol(d::Diode{L}) where L = "diode"
tikz_symbol(s::Switch{L}) where L = s.value ? "break contact" : "make contact"

# * Circuit

mutable struct Circuit
    elements::Vector{<:Element}
    connections::SparseMatrixCSC
    offsets::Vector{Int}
end

Circuit() = Circuit(Vector{Element}(), spzeros(Bool,0,0), Vector{Int}())

findelement(c::Circuit, e::Element) = findfirst(isequal(e), c.elements)

function Base.push!(c::Circuit, e::Element)
    !isnothing(findelement(c, e)) &&
        throw(ArgumentError("Element $e already present in circuit $c"))

    m = size(c.connections, 1)
    n = num_pins(e)
    M = m+n
    o = length(c.offsets) > 0 ? n + c.offsets[end] : n

    push!(c.elements, e)
    push!(c.offsets, o)

    connections = spzeros(Bool,M,M)
    copyto!(view(connections, 1:m, 1:m), c.connections)
    c.connections = connections

    c
end

Base.copy(c::Circuit) = Circuit(copy.(c.elements), copy.(c.connections), copy.(c.offsets))

offset(c::Circuit, i::Integer) = i > 1 ? c.offsets[i-1] : 0
function offset(c::Circuit, e::Element)
    i = findelement(c, e)
    isnothing(i) && throw(ArgumentError("$(e) not present in circuit $(c)"))
    offset(c, i)
end

pin(c::Circuit, e::Element, p::Pin) = offset(c, e) + pin(e, p)
pins(c::Circuit, e::Element) = offset(c, e) .+ pins(e)

function Base.show(io::IO, ::MIME"text/plain", c::Circuit)
    write(io, "Circuit with ")
    show(io, "text/plain", c.elements)
end

function Base.show(io::IO, c::Circuit)
    write(io, "Circuit(")
    show(io, c.elements)
    write(io, ")")
end

function connect!(c::Circuit, a::Int, b::Int)
    c.connections[a,b] = true
    c.connections[b,a] = true
    c
end

function connect!(c::Circuit, a::Element, ap::Pin, b::Element, bp::Pin)
    ai = pin(c, a, ap)
    bi = pin(c, b, bp)
    @debug "Connecting pin $(ap)($(ai)) of $(a) to pin $(bp)($(bi)) of $(b)"
    connect!(c, ai, bi)
end

function unique_labels!(c::Circuit)
    d = Dict{Type{<:Element},Int}()
    for e ∈ c.elements
        e.label isa Integer || throw(ArgumentError("Cannot uniquify symbolic labels"))
        et = typeof(e)
        i = get(d, et, 0) + 1
        e.label = i
        d[et] = i
    end
    c
end

function unique_nodes(c::Circuit)
    connected_nodes = reduce(|, c.connections, dims=1) |> vec |> findall
    @debug "Number of connections: $(length(connected_nodes))"

    nodes = collect(1:size(c.connections,1))

    for m in connected_nodes
        for n in findall(c.connections[:,m])
            n ≤ m && continue
            nodes[n] = nodes[m]
        end
    end

    nodes
end

function node_names(c::Circuit)
    map(unique_nodes(c)) do n
        # TODO If node is connected to ground, return :GND instead
        "Net-$(n)"
    end
end

Base.convert(::Type{Circuit}, c::Circuit) = c

# * Subcircuits

struct SubCircuit
    c::Circuit
    pins::Vector{Int}
end

function SubCircuit(c::Circuit, pins::Pair{<:Element,<:Pin}...)
    pins = map(pins) do (e,p)
        pin(c, e, p)
    end
    SubCircuit(copy(c), [pins...])
end

Base.copy(sc::SubCircuit) = SubCircuit(copy(sc.c), copy(sc.pins))

function Base.append!(c::Circuit, sc::SubCircuit)
    m = size(c.connections, 1)
    els = copy(sc.c.elements)
    cons = copy(sc.c.connections)
    m′ = size(cons, 1)
    for e in els
        push!(c, copy(e))
    end
    for i = 1:m′
        copyto!(view(c.connections, m .+ (1:m′), m+i), view(cons, :, i))
        copyto!(view(c.connections, m+i, m .+ (1:m′)), view(cons, i, :))
    end

    unique_labels!(c)
    m .+ sc.pins
end

function attach!(c::Circuit, sc::SubCircuit,
                 connections::Pair{Int,Int}...)
    m = size(c.connections, 1)
    sc_pins = append!(c, sc)
    for (a,b) in connections
        aa = m+sc.pins[a]
        c.connections[b,aa] = true
        c.connections[aa,b] = true
    end
    sc_pins
end

function attach!(c::Circuit, sc::SubCircuit,
                 connections::Pair{Int,<:Tuple{<:Element,<:Pin}}...)
    connections = map(connections) do (a,(be,bp))
        a => pin(c,be,bp)
    end
    attach!(c, sc, connections...)
end

# * TikZ

function Base.convert(TikzPicture, c::Circuit;
                      print_info::Bool=true,
                      grow_sep=2u"cm",
                      kwargs...)
    all(map(e -> e isa Circuits.TwoPort, c.elements)) ||
        throw(ArgumentError("Don't know how to generate TikZ diagrams consisting of anything but two-ports"))

    nodes = unique_nodes(c)
    connection_count = map(nodes) do n
        length(findall(isequal(n), nodes))
    end
    graph = map(1:2:length(nodes)) do i
        a,b = nodes[i],nodes[i+1]
        a == b && @info("Two-port nodes identical")
        e = c.elements[cld(i,2)]
        style = connection_count[i+1] > 2 ? "[mark]" : ""
        el_sym = tikz_symbol(e)*"={"*(print_info ? tikz_info(e) : "")*"}"
        "n$(a) --[$(el_sym)] n$(b)$(style)"
    end |> n -> join(n, ",\n")

    graph_args = convert(MIME"text/tikz",
                         TikZarg[# "empty nodes",
                                 "grow right sep"=>grow_sep,
                                 "grow left sep"=>grow_sep,
                                 "grow up sep"=>grow_sep,
                                 "grow down sep"=>grow_sep])
    graph = "\\graph$(graph_args){$(graph)};"

    options = join(map(tikz_arg, TikZarg["circuit ee IEC",
                                         "every info/.style"=>"{font=\\footnotesize}",
                                         "small circuit symbols",
                                         "layered layout",
                                         "rotate"=>90,
                                         # "xscale"=>-1,
                                         "mark/.style"=>"{fill,circle}",
                                         "mark/.default"=>"",
                                         "every node/.style"=>"inner sep=0pt,minimum size=0pt"]), ",\n")
    TikzPicture(graph; options=options, kwargs...)
end

TikzPictures.save(f::S, c::Circuit; kwargs...) where {S<:TikzPictures.SaveType} =
    save(f, convert(TikzPicture, c; preamble=preamble, kwargs...))

include("bom.jl")

# * Exports

export Element, Resistor, Capacitor, PolarCapacitor, Inductor, Diode, Switch, pin, pins,
    Circuit, connect!, unique_labels!, unique_nodes, node_names,
    SubCircuit, attach!

end # module


text:
using Luna
import LinearAlgebra: norm
import Test: @test, @testset

@testset "refractive index" begin
@testset "$gi" for gi in PhysData.gas
    @testset "$Pi bar" for Pi in (0.01, 0.1, 1, 10)
        @testset "$(λi*1e9) nm" for λi in (200e-9, 400e-9, 800e-9, 1600e-9)
            for divisor in (2, 3, 4, 5)
                divdens = PhysData.density(gi, Pi)/divisor
                divP = PhysData.pressure(gi, divdens)
                ns = PhysData.ref_index(gi, λi, Pi)
                gt = Tuple([gi for _ in 1:divisor])
                pt = Tuple([divP for _ in 1:divisor])
                nm = PhysData.ref_index(gt, λi, pt)
                nm2 = PhysData.ref_index_fun(gt)(λi, [divdens for _ = 1:divisor])
                @test (nm ≈ ns)
                @test (nm2 ≈ ns)
            end
        end
    end
end
end

@testset "propagation" begin
a = 13e-6
gas = :Ar
pres = 5
τ = 30e-15
λ0 = 800e-9
L = 5e-2

# Pressure at which argon density is exactly half of that at 5 bar
halfpres = PhysData.pressure(:Ar, PhysData.density(:Ar, pres)/2)

# Common setup
grid = Grid.RealGrid(L, λ0, (160e-9, 3000e-9), 0.5e-12)
inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)

# Single gas
m = Capillary.MarcatiliMode(a, gas, pres; loss=false)
aeff(z) = Modes.Aeff(m; z=z)
densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end
responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)

linop, βfun!, β1, αfun = LinearOps.make_const_linop(grid, m, λ0)
Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)
statsfun = Stats.default(grid, Eω, m, linop, transform; gas=gas)
output_single = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output_single)

# Mixture
m = Capillary.MarcatiliMode(a, (gas, gas), (halfpres, halfpres); loss=false)
aeff(z) = Modes.Aeff(m; z=z)
densityfun = let dens0=PhysData.density(gas, halfpres)
    z -> [dens0, dens0]
end
responses = (
    (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),),
    (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
)

linop, βfun!, β1, αfun = LinearOps.make_const_linop(grid, m, λ0)
Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)
statsfun = Stats.default(grid, Eω, m, linop, transform; gas=(gas, gas))
output_mix = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output_mix)

@test all(output_mix.data["Eω"][grid.sidx, :] .== output_single.data["Eω"][grid.sidx, :])
end


text:
function main()
    n = 1
    for i in 1:10
        println(n)
        n *= i
    end
end

main()

max_stars_repo_name: shanedrabing/project-euler

max_issues_repo_name: shanedrabing/project-euler

max_forks_repo_name: shanedrabing/project-euler


text:
export update_y_l_parallel

"""
  update l and compute y . This is a subfunction for PARSDMM.jl
"""
function update_y_l_parallel(
  x                     ::Vector{TF},
  i                     ::Integer,
  Blas_active           ::Bool,
  y                     ::Vector{Vector{TF}},
  y_old                 ::Vector{Vector{TF}},
  l                     ::Vector{Vector{TF}},
  l_old                 ::Vector{Vector{TF}},
  rho                   ::Vector{TF},
  gamma                 ::Vector{TF},
  prox                  ::Vector{Any},
  TD_OP                 ::Vector{<:Union{SparseMatrixCSC{TF,TI},joAbstractLinearOperator{TF,TF}}},
  P_sub,
  x_hat                 ::Vector{Vector{TF}},
  r_pri                 ::Vector{Vector{TF}},
  s                     ::Vector{Vector{TF}},
  set_feas              ::Vector{TF},
  feasibility_only=false::Bool
  ) where {TF<:Real,TI<:Integer}

  rho1 = TF(1.0) ./ rho[1];

  copy!(y_old[1],y[1]);
  copy!(l_old[1],l[1]);

  mul!(s[1],TD_OP[1],x) #s[1] = TD_OP[1]*x;
  if Blas_active
    if gamma[1]==1 #without relaxation
      copyto!(y[1],s[1])
      #y[1]       = prox[1]( BLAS.axpy!(-rho1[1],l[1],y[1]) )
      BLAS.axpy!(-rho1[1],l[1],y[1])
      y[1] = prox[1](y[1])
      @. r_pri[1] = -s[1]+y[1];
      BLAS.axpy!(rho[1],r_pri[1],l[1]);
    else #relaxed iterations
      @. x_hat[1] = gamma[1]*s[1]
      BLAS.axpy!(TF(1.0)-gamma[1],y[1],x_hat[1]);
      #y[1]       = copy(x_hat[1]);
      copyto!(y[1],x_hat[1])
      #y[1]       = prox[1]( BLAS.axpy!(-rho1[1],l[1],y[1]))
      BLAS.axpy!(-rho1[1],l[1],y[1])
      y[1] = prox[1](y[1])
       @. r_pri[1] = -s[1]+y[1];
      BLAS.axpy!(rho[1],y[1]-x_hat[1],l[1]);
    end
  else
    if gamma[1]==1 #without relaxation
      #y[1]       = prox[1]( s[1]-l[1]*rho1[1] );
       @. y[1]       = s[1]-l[1]*rho1[1]
       y[1]          = prox[1](y[1]);
       @. r_pri[1]   = -s[1]+y[1];
       @. l[1]       = l[1]+rho[1]*r_pri[1];
    else #relaxed iterations
       @. x_hat[1]   = gamma[1]*s[1] + ( TF(1.0)-gamma[1] )*y[1]
      #y[1]          = prox[1]( x_hat[1]-l[1]*rho1[1] );
       @. y[1]       = x_hat[1]-l[1]*rho1[1]
       y[1]          = prox[1]( y[1] )
       @. r_pri[1]   = -s[1]+y[1]
       @. l[1]       = l[1]+rho[1]*( -x_hat[1]+y[1] )
    end
  end #end blas

  #log_PARSDMM.r_pri[i,1]  = norm(r_pri[1])#./(norm(y[1])+(100*eps(TF))); #add 1-14, because y may be 0 for certain initializations of x,y,l
  #log r_pri outside this function for parallel use (temporarily)
  # if Blas_active
  #   log_PARSDMM.r_dual[i,1] = norm( BLAS.scal!(lx,rho[1],TD_OP[1]'*(y[1]-y_old[1]),1 ))
  # else
  #  log_PARSDMM.r_dual[i,1] = norm( rho[1]*(TD_OP[1]'*(y[1]-y_old[1])) )
  # end

  #log feasibility
  if (mod(i,10) == 0 && length(P_sub)==1 ) && ((feasibility_only == false && myid()<nprocs()) || feasibility_only == true)

    copy!(x_hat[1],s[1]) #use the already allocated x_hat as a temp placeholder
    P_sub[1](x_hat[1])   #P_sub mutates the input in-place: P(A*x), s=A*x, x_hat = s. The next iteration x_hat will be overwritten anyway before use
    set_feas[1] = norm( x_hat[1] .- s[1]) ./ (norm(s[1])+(100*eps(TF)));

    #equivalent to, but with less allocations
    #log_PARSDMM.set_feasibility[counter,ii] = norm( P_sub[ii](deepcopy(s[ii])) .- s[ii]) ./ (norm(s[ii])+(100*eps(TF))); #use x_hat[ii] as output, because it is already allocated and not used afterwards. The next iteration it will be overwritten anyway before use

  end

#return y,l,r_pri,s,log_PARSDMM,counter,y_old,l_old,x_hat,set_feas_vec

end #end of function

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
module MJCore

using MuJoCo_jll
using StaticArrays, UnsafeArrays, CEnum, MacroTools, Base64
using Base: last

include("exports.jl")
export AbstractMJModel, AbstractMJData, isactivated

abstract type AbstractMJModel end
abstract type AbstractMJData end

const MJArray{T,N} = Union{MArray{<:Tuple,T,N},Array{T,N},UnsafeArray{T,N}} where {T,N}
const MJVec{T} = MJArray{T,1}
const MJMat{T} = MJArray{T,2}

const ERRORSIZE = 1000
const KEY_PATH_ENV_VAR = "MUJOCO_KEY_PATH"
const KEY_ENV_VAR = "MUJOCO_KEY"

const ACTIVATED = Ref(false)
isactivated() = ACTIVATED[]

include("util.jl")

# MuJoCo types, enums, and globals
include("wrapper/mjmodel.jl")
include("wrapper/mjdata.jl")
include("wrapper/mjvisualize.jl")
include("wrapper/mjrender.jl")
include("wrapper/mjui.jl")
include("wrapper/cglobals.jl")

const Model = Union{Ptr{mjModel},mjModel,AbstractMJModel}
const Data = Union{Ptr{mjData},mjData,AbstractMJData}

Base.cconvert(::Type{Ptr{mjModel}}, m::mjModel) = pointer_from_objref(m)
Base.cconvert(::Type{Ptr{mjData}}, d::mjData) = pointer_from_objref(d)

include("wrapper/mjmujoco.jl")

function __init__()
    _init_globals()

    if !isactivated() && haskey(ENV, KEY_PATH_ENV_VAR)
        path = normpath(ENV[KEY_PATH_ENV_VAR])
        isfile(path) || @mjerror "activating with key $path: No such file or directory"
        ACTIVATED[] = Bool(mj_activate(path))
        isactivated() || @mjerror "activating with key $path: invalid key"
    elseif !isactivated() && haskey(ENV, KEY_ENV_VAR)
        mjkey = ENV[KEY_ENV_VAR]
        mktemp() do path, io
            write(io, String(base64decode(mjkey)))
            flush(io)
            ACTIVATED[] = Bool(mj_activate(path))
        end
        isactivated() || @mjerror "Key invalid"
    else
        @warn """

            No key found. Please set the env variable $KEY_PATH_ENV_VAR to point to a valid `mjkey` file.
            Alternatively, use `mj_activate(path_to_mjkey.txt)` to manually activate.
        """
    end
    nothing
end

const MJTypes = Union{
    # mjmodel.jl
    mjModel,
    mjLROpt,
    mjVFS,
    mjOption,

    # mjVisual
    mjVisual,
    mjVisual_global,
    mjVisual_quality,
    mjVisual_headlight,
    mjVisual_map,
    mjVisual_scale,
    mjVisual_rgba,

    mjStatistic,

    # mjdata.jl
    mjData,
    mjContact,
    mjWarningStat,
    mjTimerStat,
    mjSolverStat,

    # mjrender.jl
    mjrRect,
    mjrContext,

    # mjui.jl
    mjuiState,
    mjuiThemeSpacing,
    mjuiThemeColor,
    mjuiItem,
    mjuiSection,
    mjUI,
    mjuiDef,

    # mjvisualize.jl
    mjvPerturb,
    mjvCamera,
    mjvGLCamera,
    mjvGeom,
    mjvLight,
    mjvOption,
    mjvScene,
    mjvFigure,
}

function setup_comparators(T::Type)
    quote
        function Base.:(==)(x::$T, y::$T)
            for name in fieldnames($T)
                !(getfield(x, name) == getfield(y, name)) && return false
            end
            return true
        end
        function Base.isequal(x::$T, y::$T)
            for name in fieldnames($T)
                !isequal(getfield(x, name), getfield(y, name)) && return false
            end
            return true
        end
        function Base.hash(x::$T, h::UInt)
            h = hash($T, h)
            for name in fieldnames($T)
                h = hash(getfield(x, name), h)
            end
            return h
        end
    end
end

for T in Iterators.filter(T->T.mutable, union_types(MJTypes))
    @eval $(setup_comparators(T))
end

function Base.isapprox(x::T, y::T, args...; kwargs...) where {T<:MJTypes}
    for name in fieldnames(T)
        f1 = getfield(x, name)
        f2 = getfield(y, name)
        isapprox(f1, f2, args...; kwargs...) || return false
    end
    true
end

end


text:
#######################################################################
function estimatetemperature1(W, μ, T)
#######################################################################

    g = observedunitgalaxyvector(T)

    function loss(x)

        norm(x - μ - proj(W, x - μ))^2 + norm(x - proj(g, x))^2

    end

    opt = Optim.Options(show_trace = false, iterations = 1_000_000)

    result = optimize(loss, randn(3)*10, LBFGS(), opt, autodiff=:forward)

    result.minimizer, result.minimum

end

#######################################################################
function estimatetemperature2(W, μ, T)
#######################################################################

    g = observedunitgalaxyvector(T)

    function loss(s)

        norm(s*g - μ - proj(W, s*g - μ))^2

    end

    opt = Optim.Options(show_trace = true, iterations = 1_000_000)

    result = optimize(loss, 0.001, 1000.0)

    result.minimizer, result.minimum

end


text:
# This file is part of the HORIZONS.jl package
# The HORIZONS.jl package is licensed under the MIT "Expat" License
# Copyright (c) 2017: Jorge Pérez.

using Test, Dates, DelimitedFiles
using HORIZONS, Expect

import Pkg; Pkg.build("Expect"); Pkg.test("Expect")

@testset "Test connection to HORIZONS machine using Expect.jl" begin
    proc = ExpectProc(`telnet $HORIZONS_MACHINE 6775`, 15)
    @show proc.proc
    # Get main prompt and proceed, turning off paging, specifying I/O model,
    # and sending object look-up from command-line
    idx = expect!(proc, ["unknown host", "Horizons> "])
    if idx == 1
        throw("This system cannot find $HORIZONS_MACHINE")
    end
    @test idx == 2
end

@testset "Vector table generation: vec_tbl" begin
    dt0 = DateTime(2029,4,13)
    dtmax = Date(2029,4,14)
    δt = Hour(1)
    apophisraw = vec_tbl("Apophis", dt0, dtmax, δt; CSV_FORMAT=true);
    @test isa(apophisraw, String)
    @test occursin("\$\$SOE", apophisraw)
    @test occursin("\$\$EOE", apophisraw)
    @test occursin(r"\$\$SOE", apophisraw)
    @test occursin(r"\$\$EOE", apophisraw)
    mSOE = match(r"\$\$SOE", apophisraw)
    mEOE = match(r"\$\$EOE", apophisraw)
    @test mSOE.offsets == Int64[]
    @test mEOE.offsets == Int64[]

    # get everything within SOE and EOE
    apophisste = apophisraw[mSOE.offset+7:mEOE.offset-3]
    # turn into 2-dim array
    apophisarr = readdlm(IOBuffer(apophisste), ',')[:,1:end-1]

    @test typeof(apophisarr) == Array{Any,2}
    @test size(apophisarr) == (25, 11)

    # get table labels
    apophishdr = convert(Array{String,2}, strip.(    readdlm(  IOBuffer( match(r"JDTDB.*,\r\n", apophisraw).match ), ','  )[:,1:end-1]    ));
    @test typeof(apophishdr) == Array{String,2}
    @test size(apophishdr) == (1,11)
    # vcat into common 2-dim array
    apophistable = vcat(apophishdr, apophisarr)

    @test typeof(apophistable) == Array{Any,2}
    @test size(apophistable) == (26, 11)

    dt0 = Date(2000); dtmax = Date(2015); δt = Year(1);

    #test case VEC_LABELS=true
    out_str = vec_tbl("1950 DA", dt0, dtmax, δt; VEC_LABELS=true)
    @test isa(out_str, String)
    @test occursin("\$\$SOE", out_str)
    @test occursin("\$\$EOE", out_str)
    @test occursin(r"\$\$SOE.*\$\$EOE"s, out_str)
    @test occursin(r"JDTDB\r\n.*X.*Y.*Z\r\n.*VX.*VY.*VZ\r\n.*LT.*RG.*RR", out_str)
    @test occursin(r"\*+\r\nJDTDB\r\n.*X.*Y.*Z\r\n.*VX.*VY.*VZ\r\n.*LT.*RG.*RR\r\n\*+", out_str)
    @test occursin(r"\*+\r\n\$\$SOE\r\n.*\r\n\$\$EOE\r\n\*+"s, out_str)
    @test occursin(r"\$\$SOE\r\n[0-9]+.[0-9]+ = A.D. [0-9]{4}-[A-Z][a-z]{2}-[0-9]{2} [0-9]{2}:[0-9]{2}:[0-9]{2}.[0-9]{4} TDB \r\n X = ", out_str)

    #test case VEC_LABELS=false
    out_str = vec_tbl("1950 DA", dt0, dtmax, δt; VEC_LABELS=false, CENTER="Paris", REF_PLANE="E")
    @test isa(out_str, String)
    @test occursin("\$\$SOE", out_str)
    @test occursin("\$\$EOE", out_str)
    @test occursin(r"\$\$SOE.*\$\$EOE"s, out_str)
    @test occursin(r"JDTDB\r\n.*X.*Y.*Z\r\n.*VX.*VY.*VZ\r\n.*LT.*RG.*RR", out_str)
    @test occursin(r"\*+\r\nJDTDB\r\n.*X.*Y.*Z\r\n.*VX.*VY.*VZ\r\n.*LT.*RG.*RR\r\n\*+", out_str)
    @test occursin(r"\*+\r\n\$\$SOE\r\n.*\r\n\$\$EOE\r\n\*+"s, out_str)
    @test occursin(r"\$\$SOE\r\n[0-9]+.[0-9]+ = A.D. [0-9]{4}-[A-Z][a-z]{2}-[0-9]{2} [0-9]{2}:[0-9]{2}:[0-9]{2}.[0-9]{4} TDB \r\n +[0-9]+.[0-9]+E[+,-][0-9]+ +[0-9]+.[0-9]+E[+,-][0-9]+", out_str)
end

@testset "Vector table generation with CSV format: vec_tbl_csv" begin
    dt0 = Date(1950,1,1)
    dtmax = DateTime(1959, 12, 31, 11, 59, 59, 999)
    δt = Year(1)
    # 399 corresponds to Earth
    earth_tbl, earth_csv_str = vec_tbl_csv("399", dt0, dtmax, δt; VEC_TABLE = 2)
    @test typeof(earth_tbl) == Array{Any,2}
    @test size(earth_tbl) == (11, 8)
    dtmax = DateTime(1960, 1, 1, 0, 0, 0, 1)
    earth_tbl, earth_csv_str = vec_tbl_csv("399", dt0, dtmax, δt; VEC_TABLE = 2)
    @test typeof(earth_tbl) == Array{Any,2}
    @test size(earth_tbl) == (12, 8)
    earth_tbl2, earth_csv_str2 = vec_tbl_csv("399", dt0, dtmax, δt; VEC_TABLE = "2")
    @test typeof(earth_tbl2) == Array{Any,2}
    @test size(earth_tbl2) == (12, 8)
    @test earth_tbl == earth_tbl2
    @test earth_csv_str == earth_csv_str2
    # generate table with uncertainties for asteroid 1950 DA
    dt0 = Date(2000)
    dtmax = Date(2015)
    δt = Year(1)
    _1950da_tbl, _1950da_csv_str = vec_tbl_csv("1950 DA", dt0, dtmax, δt;
        VEC_TABLE = "2xa", REF_PLANE="F", CENTER="coord", COORD_TYPE="C", SITE_COORD="1,45,45")
    @test typeof(_1950da_tbl) == Array{Any,2}
    @test size(_1950da_tbl) == (17, 20)
    labels_2xa = ["JDTDB" "Calendar_Date_TDB" "X" "Y" "Z" "VX" "VY" "VZ" "X_s" "Y_s" "Z_s" "VX_s" "VY_s" "VZ_s" "A_s" "C_s" "N_s" "VA_s" "VC_s" "VN_s"]
    for i in eachindex(labels_2xa)
        _1950da_tbl[1, i] = labels_2xa[i]
    end
    #NOTE: the second object returned by vec_tbl_csv may be used to produce a DataFrame!
    # e.g.:
    # using DataFrames
    # mydataframe = readtable(IOBuffer(ea_csv_str))
end

### CI is failing currently for these test in Linux due to ftp issues on travis
### For more details, see https://blog.travis-ci.com/2018-07-23-the-tale-of-ftp-at-travis-ci
@testset "Vector table generation and write output to file: vec_tbl" begin
    dt0 = Date(1836)
    dtmax = Date(1994)
    δt = Year(5)
    # 90000033 corresponds to last Halley's apparition
    file_name = vec_tbl("90000033", "Halley.txt", dt0, dtmax, δt; CSV_FORMAT=true, ftp_verbose=true);
    @test isfile(file_name)
    @test isfile("Halley.txt")
    @test file_name == "Halley.txt"
    file_name = vec_tbl("90000033", "", dt0, dtmax, δt; ftp_verbose=true, CSV_FORMAT=true);
    @test isfile(file_name)
end

@testset "Generation and file download of small-bodies binary SPK files: smb_spk" begin
    ftp_name, local_file = smb_spk("b", "DES= 2099942;", DateTime(2021,Jan,1), DateTime(2029,Apr,13), "joe@your.domain.name")
    @test isfile(local_file)
    smb_spk("b", "DES= 2099942;", DateTime(2021,Jan,1), DateTime(2029,Apr,13), "joe@your.domain.name", "mybinaryspk.apophis")
    smb_spk("b", "DES= 2099942;", DateTime(2021,Jan,1), DateTime(2029,Apr,13), "joe@your.domain.name", "mybinaryspk.apophis", ftp_verbose=true)
    @test isfile("mybinaryspk.apophis")
    smb_spk("b", "DES= 2099942;", "2021-1-1", "2029-4-13T21:46:07.999", "joe@your.domain.name", "2099942_.bsp")
    smb_spk("b", "DES= 2099942;", "2021-1-1", "2029-4-13T21:46:07.999", "joe@your.domain.name", "2099942_.bsp", ftp_verbose=true)
    @test isfile("2099942_.bsp")
end

@testset "Generation and file download of small-bodies binary SPK files: smb_spk_ele" begin
    elements_str = "EPOCH=2449526.5 EC=.6570220840219289 QR=.5559654280797371 TP=2449448.890787227 OM=78.10766874391773 W=77.40198125423228 IN=24.4225258251465"
    smb_spk_ele("b", "1990 MU", DateTime(2021,Jan,1), DateTime(2022,Jan,1), elements_str, "joe@your.domain.name", "1990mu.bsp")
    @test isfile("1990mu.bsp")
    elements_str = "EPOCH= 2454733.5 EC=.1911952942528226 QR=.7460724385331012 TP=2454894.912507658200 OM=204.4460284242489 W=126.401880836064 IN=3.331369228495799"
    ftp_name, local_file = smb_spk_ele("b", "99942", DateTime(2021,Jan,1), DateTime(2022,Jan,1), elements_str, "joe@your.domain.name")
    @test isfile(local_file)
    epoch = 2454733.5
    ec = 0.1911952942528226
    qr = 0.7460724385331012
    tp = 2454894.912507658200
    om = 204.4460284242489
    w = 126.401880836064
    inc = 3.331369228495799
    ftp_name, local_file = smb_spk_ele("b", "99942", DateTime(2021,Jan,1), DateTime(2022,Jan,1), epoch, ec, qr, tp, om, w, inc, "joe@your.domain.name")
    @test isfile(local_file)
    ftp_name, local_file = smb_spk_ele("b", "99942", DateTime(2021,Jan,1), DateTime(2022,Jan,1), elements_str, "joe@your.domain.name", ftp_verbose=true)
    @test isfile(local_file)
end

@testset "Test for erroneous arguments" begin
    @test_throws ArgumentError vec_tbl_csv("erroneous-input", Date(2000), Date(2010), Year(1))
    @test_throws ArgumentError vec_tbl_csv(99942, Date(2000), Date(2010), Year(1); CENTER="nomatch")
    @test_throws ArgumentError vec_tbl_csv(499, Date(2009), Date(2010), Year(1); VEC_TABLE = 1, CENTER="mars")
    dt0 = Date(2000); dtmax = Date(2015); δt = Year(1)
    @test_throws ArgumentError vec_tbl_csv("1950 DA", dt0, dtmax, δt; VEC_TABLE = "2xa", CENTER="coord", COORD_TYPE="w")
    @test_throws ArgumentError vec_tbl_csv("1950 DA", dt0, dtmax, δt; VEC_TABLE = "2xa", CENTER="coord", COORD_TYPE="%")
    @test_throws ArgumentError vec_tbl_csv("1950 DA", dt0, dtmax, δt; VEC_TABLE = "2xa", CENTER="coord", COORD_TYPE="")
    @test_throws ArgumentError vec_tbl_csv("1950 DA", dt0, dtmax, δt; VEC_TABLE = "2xa", CENTER="coord", COORD_TYPE="c", SITE_COORD="a,b,c")
    @test_throws ArgumentError vec_tbl_csv("1950 DA", dt0, dtmax, δt; VEC_TABLE = "2xa", CENTER="coord", COORD_TYPE="c", SITE_COORD="")
    @test_throws ArgumentError vec_tbl_csv("1950 DA", dt0, dtmax, δt; REF_PLANE="T", VEC_TABLE = "2xa", CENTER="coord", COORD_TYPE="c", SITE_COORD="10,1,1")
    @test_throws ArgumentError vec_tbl_csv("1950 DA", "1400-1-1", dtmax, δt)
    @test_throws ArgumentError vec_tbl_csv("1950 DA", "bad-start-time", dtmax, δt)
    @test_throws ArgumentError vec_tbl_csv("1950 DA", dt0, "3500-2-1", δt)
    @test_throws ArgumentError vec_tbl_csv("1950 DA", dt0, "bad-stop-time", δt)
    @test_throws ArgumentError vec_tbl_csv("1950 DA", dt0, dtmax, "1 w")
    @test_throws ArgumentError vec_tbl_csv("1950 DA", dt0, dtmax, "bad-step_size")
    @test_throws TypeError vec_tbl_csv("1950 DA", dt0, dtmax, δt; VEC_LABELS=0)
    @test_throws TypeError vec_tbl("1950 DA", dt0, dtmax, δt; VEC_LABELS=0)
end


text:
using DifferentialEquations
using LinearAlgebra, DiffEqSensitivity, Optim
using Flux: flatten, params
using DiffEqFlux, Flux
using Plots
using Flux: train!
using NNlib
using LaTeXStrings

# Generate Data from subexpontial growth model

function model1(du, u, p, t)
    r, α = p
    du .= r .* u .^ α
end
u_0 = [1.0]
p_data = [0.2, 0.5]
tspan_data = (0.0, 30)
prob_data = ODEProblem(model1, u_0, tspan_data, p_data)
data_solve = solve(prob_data, Tsit5(), abstol = 1e-12, reltol = 1e-12, saveat = 1)
data_withoutnois = Array(data_solve)
data = data_withoutnois #+ Float32(2e-1)*randn(eltype(data_withoutnois), size(data_withoutnois))
tspan_predict = (0.0, 40)
prob_predict = ODEProblem(model1, u_0, tspan_predict, p_data)
test_data = solve(prob_predict, Tsit5(), abstol = 1e-12, reltol = 1e-12, saveat = 1)
plot(data_solve)

ann_node = FastChain(FastDense(1, 10, tanh), FastDense(10, 1))
p = Float64.(initial_params(ann_node))
function model1_nn(du, u, p, t)
    du[1] = 0.1 * ann_node(t, p)[1] * u[1] - 0.1 * u[1]
end
prob_nn = ODEProblem(model1_nn, u_0, tspan_data, p)
function train(θ)
    Array(concrete_solve(prob_nn, Tsit5(), u_0, θ, saveat = 1,
        abstol = 1e-6, reltol = 1e-6))#,sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
end
#println(train(p))
function loss(θ)
    pred = train(θ)
    sum(abs2, (data .- pred)), pred # + 1e-5*sum(sum.(abs, params(ann)))
end
#println(loss(p))
const losses = []
callback(θ, l, pred) = begin
    push!(losses, l)
    if length(losses) % 50 == 0
        println(losses[end])
    end
    false
end
res1_node = DiffEqFlux.sciml_train(loss, p, ADAM(0.02), cb = callback, maxiters = 500)
res2_node = DiffEqFlux.sciml_train(loss, res1_node.minimizer, BFGS(initial_stepnorm = 0.01), cb = callback, maxiters = 2000)
println(res2_node.minimizer)
prob_nn2 = ODEProblem(model1_nn, u_0, tspan_predict, res2_node.minimizer)
s_nn = solve(prob_nn2, Tsit5(), saveat = 1)

# I(t)
scatter(data_solve.t, data[1, :], label = "Training Data")
plot!(test_data, label = "Real Data")
plot!(s_nn, label = "Neural Networks")
xlabel!("t(day)")
ylabel!("I(t)")
title!("Subexpotential Model(I(t))")
savefig("Figures/subexpontialIt.png")
# R(t)

plot((2 ./ sqrt.(test_data) .+ 1)', label = L"R_t = 2/\sqrt{u}+1")
plot!((2 ./ sqrt.(s_nn) .+ 1)', label = L"R_t = NN(t)")
xlabel!("t(day)")
ylabel!("Effective Reproduction Number")
title!("Subexpotential Model(Rt)")
savefig("Figures/subexpontialRt.png")


text:
using DataFrames,Missings

export cursor_dicts,dicts_df

function cursor_dicts(cursor::Mongo.MongoCursor)

arr=[]
for r in cursor
d=Dict()

        d=LibBSON.dict(r)
        try merge!(d,d["_id"]) catch "" end
        try delete!(d,"_id") catch "" end
        push!(arr,d)

end

    return arr
end

function dicts_df(dicts::Array,header=[])

df=DataFrame()

    ###check if header defined
    header==[] ? k=keys(dicts[1]) : k=header

        for rk in k

            ark=[]

                for row in dicts

                    haskey(row,rk)? push!(ark,(row[rk])) : push!(ark,missing)

                end

            #=if typeof(ark[1])==String
            df[Symbol(rk)]=@data(map(String,ark))
          else

          end
          =#
            df[Symbol(rk)]=(map(x->x,ark))

        end

    return df_stab_types(df)
end

############# Function stab types
function df_stab_types(df)
fields=names(df)
rows=size(df,1)
def_type=Any
sample_n_arr=map(x->Int(ceil(x)),(rand(100)*rows))
sample_arr=[]

for f in fields

    for r in sample_n_arr

        push!(sample_arr,typeof(df[f][r]))

    end

sample_arr=unique(sample_arr)

    for s in sample_arr

        if s <: AbstractFloat
            def_type=Float64
            break
        end

        if s == Date
            def_type=Date
            break
        end

        if s == DateTime
            def_type=DateTime
            break
        end

        if s <: Integer
            def_type=Int64
            break
        end

        if s <: String
            def_type=String
            break
        end

    end

try df[f]=convert(Array{def_type,1}, df[f])  catch df[f]=convert(Array{Any,1}, df[f]) end

def_type=Any
sample_arr=[]
end

    return df
end


text:
#!/bin/env julia
# Options apply to all layers, weights and biases.
# (Except for regularization which applies to w, not b)

using HDF5
using ArgParse
using Knet

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "x"
        help = "HDF5 file for input"
        required = true
        "net"
        help = "Comma separated list of HDF5 layer files"
        required = true
        "y"
        help = "HDF5 file for desired output"
        required = true
        "out"
        help = "File prefix for trained output layer files"
        required = true
        "--batch"
        help = "Minibatch size"
        arg_type = Int
        default = 128
        "--adagrad"
        help = "If nonzero apply adagrad using arg as epsilon"
        arg_type = Float32
        default = 0f0
        "--dropout"
        help = "Dropout probability"
        arg_type = Float32
        default = 0f0
        "--iters"
        help = "If nonzero limits number of batches"
        arg_type = Int
        default = 0
        "--l1reg"
        help = "L1 regularization parameter"
        arg_type = Float32
        default = 0f0
        "--l2reg"
        help = "L2 regularization parameter"
        arg_type = Float32
        default = 0f0
        "--learningRate"
        help = "Learning rate"
        arg_type = Float32
        default = 0.01f0
        "--loss"
        help = "Loss function"
        arg_type = String
        default = "softmaxloss"
        "--maxnorm"
        help = "If nonzero upper limit on weight matrix row norms"
        arg_type = Float32
        default = 0f0
        "--momentum"
        help = "Momentum"
        arg_type = Float32
        default = 0f0
        "--nesterov"
        help = "Apply nesterov's accelerated gradient"
        arg_type = Float32
        default = 0f0
        "--nogpu"
        help = "Do not use gpu"
        action = :store_true
    end
    args = parse_args(s)
    for (arg,val) in args
        print("$arg:$val ")
    end
    println("")
    return args
end

function main()
    args = parse_commandline()
    Knet.gpu(!args["nogpu"])
    args["nogpu"] && blas_set_num_threads(20)
    x = h5read(args["x"], "/data");
    y = h5read(args["y"], "/data");
    net = map(l->Op(l), split(args["net"],','))
    for (a,v) in args
        if !in(a, ["x","y","nogpu","net","batch","iters","loss","out"])
            setparam!(net, symbol(a), v)
        end
    end
    @time train(net, x, y; batch=args["batch"], iters=args["iters"], loss=eval(parse(args["loss"])))
    out = args["out"]
    for l=1:length(net)
        h5write("$out$l.h5", net[l]);
    end
end

main()


text:
using Revise
using Test
using NNFEM
using PyCall
using PyPlot
using JLD2
using ADCME
using LinearAlgebra
reset_default_graph()
include("nnutil.jl")

testtype = "NeuralNetwork1D"

prop = Dict("name"=> testtype, "rho"=> 0.1, "E"=> 200.0, "B"=> 10.0,
            "sigmaY"=>0.300, "K"=>1/9*200, "A0"=> 1.0, "eta"=> 10.0, "nn"=>post_nn)

include("NNTrussPull_Domain.jl")

# domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
@load "Data/domain1.jld2" domain
state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)
assembleMassMatrix!(globdat, domain)

nntype = "ae_scaled"
# nntype = "linear"
n_data = [1,2,4,5]
losses = Array{PyObject}(undef, length(n_data))
for (i, ni) in enumerate(n_data)
    state_history, fext_history = read_data("$(@__DIR__)/Data/$i.dat")
    losses[i] = DynamicMatLawLoss(domain, globdat, state_history, fext_history, nn,Δt)
end
loss = sum(losses)

sess = Session(); init(sess)
@show run(sess, loss)
BFGS!(sess, loss, 2000)
ADCME.save(sess, "Data/trained_nn_fem.mat")


text:
using SubmodularMaximization

using PyPlot
using LinearAlgebra

num_sample = 1000
limits = [0 1; 0 1]

function analyze_pdf(p, name)
  println("Evaluating $name")

  println("PDF at center")
  @show evaluate_pdf(p, [0.5, 0.5])

  figure()
  @time visualize_pdf(p)
  title("$name")

  println("Sampling PDF")
  @time samples = hcat(map(x->sample_reject(p, limits), 1:num_sample)...)

  scatter(samples[1,:]', samples[2,:]', color = "k", marker = (5, 2, 0))
end

g = Gaussian([0.5, 0.5], Diagonal([0.01, 0.1]))
analyze_pdf(g, "Gaussian")

bimodal_mixture = GaussianMixture([1.0, 1.5],
                                  [Gaussian([0.2, 0.5], Diagonal([0.001, 0.1])),
                                   Gaussian([0.8, 0.5], Diagonal([0.001, 0.1]))])
analyze_pdf(bimodal_mixture, "Bimodal Mixture")

analyze_pdf(standard_mixture(), "Standard Mixture")

nothing

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing