text:
using Base.Meta: parse

#################
# Overload of ~ #
#################

# TODO: Replace this macro, see issue #514
"""
Usage: @VarName x[1,2][1+5][45][3]
  return: (:x,[1,2],6,45,3)
"""
macro VarName(expr::Union{Expr, Symbol})
    ex = deepcopy(expr)
    isa(ex, Symbol) && return var_tuple(ex)
    (ex.head == :ref) || throw("VarName: Mis-formed variable name $(expr)!")
    inds = :(())
    while ex.head == :ref
        if length(ex.args) >= 2
            pushfirst!(inds.args, Expr(:vect, ex.args[2:end]...))
            end
        ex = ex.args[1]
        isa(ex, Symbol) && return var_tuple(ex, inds)
    end
    throw("VarName: Mis-formed variable name $(expr)!")
end
function var_tuple(sym::Symbol, inds::Expr=:(()))
    return esc(:($(QuoteNode(sym)), $inds, $(QuoteNode(gensym()))))
end

"""
    generate_observe(observation, distribution)

Generate an observe expression.
"""
function generate_observe(observation, distribution)
    return esc(
        quote
            vi.logp += Turing.observe(
                sampler,
                $(distribution),
                $(observation),
                vi
            )
        end
    )
end

"""
    generate_assume(variable, distribution, syms)

Generate an assume expression.
"""
function generate_assume(variable, distribution, syms)
    return esc(
        quote
            varname = Turing.VarName(vi, $syms, "")
            ($(variable), _lp) = if isa($(distribution), Vector)
                Turing.assume(
                    sampler,
                    $(distribution),
                    varname,
                    $(variable),
                    vi
                )
            else
                Turing.assume(
                    sampler,
                    $(distribution),
                    varname,
                    vi
                )
            end
            vi.logp += _lp
        end
    )
end

function generate_assume(variable::Expr, distribution)
    return esc(
        quote
            sym, idcs, csym = @VarName $variable
            csym_str = string(Turing._compiler_[:name])*string(csym)
            indexing = mapfoldl(string, *, idcs, init = "")
            varname = Turing.VarName(vi, Symbol(csym_str), sym, indexing)

            $(variable), _lp = Turing.assume(
                sampler,
                $(distribution),
                varname,
                vi
            )
            vi.logp += _lp
        end
    )
end

"""
    macro: @~ var Distribution()

Tilde notation macro. This macro constructs Turing.observe or
Turing.assume calls depending on the left-hand argument.
Note that the macro is interconnected with the @model macro and
assumes that a `compiler` struct is available.

Example:
```julia
@~ x Normal()
```
"""
macro ~(left, right)
    return tilde(left, right)
end

function tilde(left, right)
    return generate_observe(left, right)
end

function tilde(left::Symbol, right)
    # Check if left-hand side is a observation.
    if left in Turing._compiler_[:args]
        if !(left in Turing._compiler_[:dvars])
            @info " Observe - `$(left)` is an observation"
            push!(Turing._compiler_[:dvars], left)
        end

        return generate_observe(left, right)
    else
        # Assume it is a parameter.
        if !(left in Turing._compiler_[:pvars])
            msg = " Assume - `$(left)` is a parameter"
            if isdefined(Main, left)
                msg  *= " (ignoring `$(left)` found in global scope)"
            end

            @info msg
            push!(Turing._compiler_[:pvars], left)
        end

        sym, idcs, csym = @VarName(left)
        csym = Symbol(Turing._compiler_[:name], csym)
        syms = Symbol[csym, left]

        return generate_assume(left, right, syms)
    end
end

function tilde(left::Expr, right)
    vsym = getvsym(left)
    @assert isa(vsym, Symbol)

    if vsym in Turing._compiler_[:args]
        if !(vsym in Turing._compiler_[:dvars])
            @info " Observe - `$(vsym)` is an observation"
            push!(Turing._compiler_[:dvars], vsym)
        end

        return generate_observe(left, right)
    else
        if !(vsym in Turing._compiler_[:pvars])
            msg = " Assume - `$(vsym)` is a parameter"
            if isdefined(Main, vsym)
                msg  *= " (ignoring `$(vsym)` found in global scope)"
            end

            @info msg
            push!(Turing._compiler_[:pvars], vsym)
        end

        return generate_assume(left, right)
    end
end

#################
# Main Compiler #
#################

"""
    @model(name, fbody)

Macro to specify a probabilistic model.

Example:

```julia
@model Gaussian(x) = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0,sqrt.(s))
    for i in 1:length(x)
        x[i] ~ Normal(m, sqrt.(s))
    end
    return (s, m)
end
```

Compiler design: `sample(fname(x,y), sampler)`.
```julia
fname(x=nothing,y=nothing; compiler=compiler) = begin
    ex = quote
        # Pour in kwargs for those args where value != nothing.
        fname_model(vi::VarInfo, sampler::Sampler; x = x, y = y) = begin
            vi.logp = zero(Real)

            # Pour in model definition.
            x ~ Normal(0,1)
            y ~ Normal(x, 1)
            return x, y
        end
    end
    return Main.eval(ex)
end
```
"""
macro model(fexpr)

    # translate all ~ occurences to macro calls
    fexpr = translate(fexpr)

    # extract model name (:name), arguments (:args), (:kwargs) and definition (:body)
    modeldef = MacroTools.splitdef(fexpr)

    # function body of the model is empty
    if all(l -> isa(l, LineNumberNode), modeldef[:body].args)
        @warn("Model definition seems empty, still continue.")
    end

    # construct compiler dictionary
    compiler = Dict(
        :name => modeldef[:name],
        :closure_name => Symbol(modeldef[:name], :_model),
        :args => [],
        :kwargs => modeldef[:kwargs],
        :dvars => Set{Symbol}(),
        :pvars => Set{Symbol}()
    )

    # Manipulate the function arguments.
    fargs = deepcopy(vcat(modeldef[:args], modeldef[:kwargs]))
    for i in 1:length(fargs)
        if isa(fargs[i], Symbol)
            fargs[i] = Expr(:kw, fargs[i], :nothing)
        end
    end

    # Construct closure.
    closure = MacroTools.combinedef(
        Dict(
            :name => compiler[:closure_name],
            :kwargs => [],
            :args => [
                :(vi::Turing.VarInfo),
                :(sampler::Turing.AnySampler)
            ],
            # Initialise logp in VarInfo.
            :body => Expr(:block, :(vi.logp = zero(Real)), modeldef[:body].args...)
        )
    )

    # Construct aliases.
    alias1 = MacroTools.combinedef(
        Dict(
            :name => compiler[:closure_name],
            :args => [:(vi::Turing.VarInfo)],
            :kwargs => [],
            :body => :(return $(compiler[:closure_name])(vi, Turing.SampleFromPrior()))
        )
    )

    alias2 = MacroTools.combinedef(
        Dict(
            :name => compiler[:closure_name],
            :args => [:(sampler::Turing.AnySampler)],
            :kwargs => [],
            :body => :(return $(compiler[:closure_name])(Turing.VarInfo(), Turing.SampleFromPrior()))
        )
    )

    alias3 = MacroTools.combinedef(
        Dict(
            :name => compiler[:closure_name],
            :args => [],
            :kwargs => [],
            :body => :(return $(compiler[:closure_name])(Turing.VarInfo(), Turing.SampleFromPrior()))
        )
    )

    # Add definitions to the compiler.
    compiler[:closure] = closure
    compiler[:alias1] = alias1
    compiler[:alias2] = alias2
    compiler[:alias3] = alias3

    # Construct user function.
    modelfun = MacroTools.combinedef(
        Dict(
            :name => compiler[:name],
            :kwargs => [Expr(:kw, :compiler, compiler)],
            :args => fargs,
            :body => Expr(:block,
                            quote
                                Turing.eval(:(_compiler_ = deepcopy($compiler)))
                                # Copy the expr of function definition and callbacks
                                closure = Turing._compiler_[:closure]
                                alias1 = Turing._compiler_[:alias1]
                                alias2 = Turing._compiler_[:alias2]
                                alias3 = Turing._compiler_[:alias3]
                                modelname = Turing._compiler_[:closure_name]
                            end,
                            # Insert argument values as kwargs to the closure
                            map(data_insertion, fargs)...,
                            # Eval the closure's methods globally and return it
                            quote
                                Main.eval(Expr(:(=), modelname, closure))
                                Main.eval(alias1)
                                Main.eval(alias2)
                                Main.eval(alias3)
                                return $(compiler[:closure_name])
                            end,
                        )
        )
    )

    return esc(modelfun)
end

####################
# Helper functions #
####################

function data_insertion(k)
        if isa(k, Symbol)
            _k = k
        elseif k.head == :kw
            _k = k.args[1]
        else
        return :()
        end

    return  quote
                if $_k == nothing
                    # Notify the user if an argument is missing.
                    @warn("Data `"*$(string(_k))*"` not provided, treating as parameter instead.")
                else
                    if $(QuoteNode(_k)) ∉ Turing._compiler_[:args]
                        push!(Turing._compiler_[:args], $(QuoteNode(_k)))
                    end
                    closure = Turing.setkwargs(closure, $(QuoteNode(_k)), $_k)
                end
            end
end

function setkwargs(fexpr::Expr, kw::Symbol, value)
    # Split up the function definition.
    funcdef = MacroTools.splitdef(fexpr)

    # Add the new keyword argument.
    push!(funcdef[:kwargs], Expr(:kw, kw, value))

    # Recompose the function.
    return MacroTools.combinedef(funcdef)
end

getvsym(s::Symbol) = s
function getvsym(expr::Expr)
    @assert expr.head == :ref "expr needs to be an indexing expression, e.g. :(x[1])"
    return getvsym(expr.args[1])
end

translate!(ex::Any) = ex
function translate!(ex::Expr)
    if ex.head === :call && ex.args[1] === :(~)
        ex.head = :macrocall
        ex.args[1] = Symbol("@~")
        insert!(ex.args, 2, LineNumberNode(@__LINE__))
    else
        map(translate!, ex.args)
    end
    return ex
end
translate(ex::Expr) = translate!(deepcopy(ex))


text:
import OrdinaryDiffEq

"""
    schroedinger(psi0, tspan, H;
    \treltol=1.0e-12, abstol=1.0e-12, fout=nothing, alg=OrdinaryDiffEq.Tsit5())

Integrate Schroedinger equation to evolve states or compute propagators

`` \\frac{d}{d t}\\vert\\psi(t)\\rangle = - \\frac{i}{\\hbar} \\hat H\\vert\\psi(t)\\rangle, \\quad \\hbar = 1``.

# Arguments
* `psi0`: Initial state vector (can be a bra or a ket) or initial propagator.
* `tspan`: Vector specifying the points of time for which output should be displayed.
* `H`: Arbitrary operator specifying the Hamiltonian.
* `reltol`: Relative tolerance for OrdinaryDiffEq solver and its inner states.
* `abstol`: Absolute tolerance for OrdinaryDiffEq solver and its inner states.
* `fout=nothing`: If given, this function `fout(t, psi)` is called every time
        an output should be displayed. ATTENTION: The state `psi` is neither
        normalized nor permanent! It is still in use by the ode solver and
        therefore must not be changed.
* `alg`: Algorithm with which OrdinaryDiffEq will solve Schroedinger equation.
"""
function schroedinger(
    psi0::T,
    tspan::Array,
    H::AbstractOperator{B,B};
    reltol::Float64 = 1.0e-12,
    abstol::Float64 = 1.0e-12,
    alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.Tsit5(),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Union{AbstractOperator{B,B},StateVector{B}}}
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    dschroedinger_(t, psi::T, dpsi::T) = dschroedinger(psi, H, dpsi)
    x0 = psi0.data
    state = copy(psi0)
    dstate = copy(psi0)
    OpenQuantumSystems.integrate(
        tspan_,
        dschroedinger_,
        x0,
        state,
        dstate,
        fout;
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

"""
    schroedinger_dynamic(psi0, tspan, f;
    \treltol=1.0e-12, abstol=1.0e-12, fout=nothing, alg=OrdinaryDiffEq.Tsit5())

Integrate time-dependent Schroedinger equation to evolve states or compute propagators

``\\frac{d}{d t}\\vert\\psi(t)\\rangle = - \\frac{i}{\\hbar} \\hat H(t)\\vert\\rho(t)\\rangle, \\quad \\hbar = 1``.

# Arguments
* `psi0`: Initial state vector (can be a bra or a ket) or initial propagator.
* `tspan`: Vector specifying the points of time for which output should be displayed.
* `f`: Function `f(t, psi) -> H` returning the time and or state dependent Hamiltonian.
* `reltol`: Relative tolerance for OrdinaryDiffEq solver and its inner states.
* `abstol`: Absolute tolerance for OrdinaryDiffEq solver and its inner states.
* `fout=nothing`: If given, this function `fout(t, psi)` is called every time
        an output should be displayed. ATTENTION: The state `psi` is neither
        normalized nor permanent! It is still in use by the ode solver and
        therefore must not be changed.
* `alg`: Algorithm with which OrdinaryDiffEq will solve Schroedinger equation.
"""
function schroedinger_dynamic(
    psi0::T,
    tspan::Array,
    f::Function;
    reltol::Float64 = 1.0e-6,
    abstol::Float64 = 1.0e-8,
    alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.DP5(),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {T<:Union{StateVector,AbstractOperator}}
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    dschroedinger_(t, psi::T, dpsi::T) = dschroedinger_dynamic(t, psi, f, dpsi)
    x0 = psi0.data
    state = copy(psi0)
    dstate = copy(psi0)
    OpenQuantumSystems.integrate(
        tspan_,
        dschroedinger_,
        x0,
        state,
        dstate,
        fout;
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

recast!(x::D, psi::StateVector{B,D}) where {B<:Basis,D} = (psi.data = x);
recast!(psi::StateVector{B,D}, x::D) where {B<:Basis,D} = nothing

function dschroedinger(
    psi::Union{Ket{B},AbstractOperator{B,B}},
    H::AbstractOperator{B,B},
    dpsi::Union{Ket{B},AbstractOperator{B,B}},
) where {B<:Basis}
    OpenQuantumSystems.mul!(dpsi, H, psi, eltype(psi)(-im), zero(eltype(psi)))
    return dpsi
end

function dschroedinger(psi::Bra{B}, H::AbstractOperator{B,B}, dpsi::Bra{B}) where {B<:Basis}
    OpenQuantumSystems.mul!(dpsi, psi, H, eltype(psi)(im), zero(eltype(psi)))
    return dpsi
end

function dschroedinger_dynamic(
    t,
    psi0::T,
    f::Function,
    dpsi::T,
) where {T<:Union{StateVector,AbstractOperator}}
    H = f(t, psi0)
    dschroedinger(psi0, H, dpsi)
end

function check_schroedinger(psi::Ket, H::AbstractOperator)
    check_multiplicable(H, psi)
    check_samebases(H)
end

function check_schroedinger(psi::Bra, H::AbstractOperator)
    check_multiplicable(psi, H)
    check_samebases(H)
end

ext: jl


text:
using MPI
using OrderedCollections
using StaticArrays
using Test

using CLIMAParameters
using CLIMAParameters.Planet: ρ_cloud_liq, ρ_cloud_ice, cp_l, cp_i, T_0, LH_f0
using CLIMAParameters.Atmos.Microphysics: K_therm
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using ClimateMachine
using ClimateMachine.Land.SoilHeatParameterizations
using ClimateMachine.Land

@testset "Land heat parameterizations" begin
    FT = Float64
    # Density of liquid water (kg/m``^3``)
    _ρ_l = FT(ρ_cloud_liq(param_set))
    # Density of ice water (kg/m``^3``)
    _ρ_i = FT(ρ_cloud_ice(param_set))
    # Volum. isobaric heat capacity liquid water (J/m3/K)
    _ρcp_l = FT(cp_l(param_set) * _ρ_l)
    # Volumetric isobaric heat capacity ice (J/m3/K)
    _ρcp_i = FT(cp_i(param_set) * _ρ_i)
    # Reference temperature (K)
    _T_ref = FT(T_0(param_set))
    # Latent heat of fusion at ``T_0`` (J/kg)
    _LH_f0 = FT(LH_f0(param_set))
    # Thermal conductivity of dry air
    κ_air = FT(K_therm(param_set))

    @test temperature_from_ρe_int(5.4e7, 0.05, 2.1415e6, param_set) ==
          FT(_T_ref + (5.4e7 + 0.05 * _ρ_i * _LH_f0) / 2.1415e6)

    @test volumetric_heat_capacity(0.25, 0.05, 1e6, param_set) ==
          FT(1e6 + 0.25 * _ρcp_l + 0.05 * _ρcp_i)

    @test volumetric_internal_energy(0.05, 2.1415e6, 300.0, param_set) ==
          FT(2.1415e6 * (300.0 - _T_ref) - 0.05 * _ρ_i * _LH_f0)

    @test saturated_thermal_conductivity(0.25, 0.05, 0.57, 2.29) ==
          FT(0.57^(0.25 / (0.05 + 0.25)) * 2.29^(0.05 / (0.05 + 0.25)))

    @test saturated_thermal_conductivity(0.0, 0.0, 0.57, 2.29) == FT(0.0)

    @test relative_saturation(0.25, 0.05, 0.4) == FT((0.25 + 0.05) / 0.4)

    # Test branching in kersten_number
    soil_param_functions = SoilParamFunctions(
        FT;
        porosity = 0.2,
        ν_ss_gravel = 0.1,
        ν_ss_om = 0.1,
        ν_ss_quartz = 0.1,
        κ_solid = 0.1,
        ρp = 1.0,
    )
    # ice fraction = 0
    @test kersten_number(0.0, 0.75, soil_param_functions) == FT(
        0.75^((FT(1) + 0.1 - 0.24 * 0.1 - 0.1) / FT(2)) *
        (
            (FT(1) + exp(-18.1 * 0.75))^(-FT(3)) -
            ((FT(1) - 0.75) / FT(2))^FT(3)
        )^(FT(1) - 0.1),
    )

    # ice fraction ~= 0
    @test kersten_number(0.05, 0.75, soil_param_functions) ==
          FT(0.75^(FT(1) + 0.1))

    @test thermal_conductivity(1.5, 0.7287, 0.7187) ==
          FT(0.7287 * 0.7187 + (FT(1) - 0.7287) * 1.5)

    @test volumetric_internal_energy_liq(300.0, param_set) ==
          FT(_ρcp_l * (300.0 - _T_ref))

    @test k_solid(FT(0.5), FT(0.25), FT(2.0), FT(3.0), FT(2.0)) ==
          FT(2)^FT(0.5) * FT(2)^FT(0.25) * FT(3.0)^FT(0.25)

    @test ksat_frozen(FT(0.5), FT(0.1), FT(0.4)) ==
          FT(0.5)^FT(0.9) * FT(0.4)^FT(0.1)

    @test ksat_unfrozen(FT(0.5), FT(0.1), FT(0.4)) ==
          FT(0.5)^FT(0.9) * FT(0.4)^FT(0.1)

    @test k_dry(param_set, soil_param_functions) ==
          ((FT(0.053) * FT(0.1) - κ_air) * FT(0.8) + κ_air * FT(1.0)) /
          (FT(1.0) - (FT(1.0) - FT(0.053)) * FT(0.8))
end


text:
using Printf

mutable struct Point
    x::Float64 # meters
    y::Float64 # meters
    z::Float64 # meters
    vx::Float64 # meters/second
    vy::Float64 # meters/second
    vz::Float64 # meters/second
    mass::Float64 # kg
    fx::Float64 # N - reset every iteration
    fy::Float64 # N
    fz::Float64 # N
end

mutable struct Spring
    k::Float64 # N/m
    p1::Int # Index of first point
    p2::Int # Index of second point
    l0::Float64 # meters
end

function main()
    kNumPerSide::Int64 = 2
    kSpring::Float64 = 500.0
    kGround::Float64 = 100000.0
    kOscillationFrequency::Float64 = 0#10000.0
    kDropHeight::Float64 = 0.2

    t::Float64 = 0.0
    points::Array{Point} = []
    springs::Array{Spring} = []
    cache::Dict{Int64, Dict{Int64, Dict{Int64, Point}}} = Dict()

    # Create the points
    for x = 1:kNumPerSide
        for y = 1:kNumPerSide
            for z = 1:kNumPerSide
                # (0,0,0) or (0.1,0.1,0.1) and all combinations
                p::Point = Point(x / 10.0, kDropHeight + y / 10.0, z / 10.0, 0, 0, 0, 0.1, 0, 0, 0)
                push!(points, p)
                if !haskey(cache, x)
                    cache[x] = Dict()
                end
                if !haskey(cache[x], y)
                    cache[x][y] = Dict()
                end
                cache[x][y][z] = p
            end
        end
    end

    connected::Dict{Int64, Array{Int64}} = Dict()
    #Create the springs
    for x = 1:kNumPerSide
        for y = 1:kNumPerSide
            for z = 1:kNumPerSide
                p1index::Int = z + kNumPerSide * (y - 1) + kNumPerSide * kNumPerSide * (x - 1)
                if !haskey(connected, p1index)
                    connected[p1index] = []
                end
                p1::Point = cache[x][y][z]
                for x1 = (x - 1):(x+1)
                    if x1 == kNumPerSide + 1 || x1 <= 0
                        continue
                    end
                    for y1 = (y - 1):(y+1)
                        if y1 == kNumPerSide + 1 || y1 <= 0
                            continue
                        end
                        for z1 = z:(z+1)
                            if z1 == kNumPerSide + 1 || z1 <= 0 || (x1 == x && y1 == y && z1 == z)
                                continue
                            end
                            p2index::Int = z1 + kNumPerSide * (y1 - 1) + kNumPerSide * kNumPerSide * (x1 - 1)
                            if !haskey(connected, p2index)
                                connected[p2index] = []
                            elseif p2index in connected[p1index]
                                continue
                            end
                            push!(connected[p1index], p2index)
                            push!(connected[p2index], p1index)
                            p2::Point = cache[x1][y1][z1]
                            length::Float64 = sqrt(abs2(p1.x - p2.x) + abs2(p1.y - p2.y) + abs2(p1.z - p2.z))
                            push!(springs, Spring(kSpring, p1index, p2index, length))
                        end
                    end
                end
            end
        end
    end
    staticFriction::Float64 = 0.5
    kineticFriction::Float64 = 0.3
    dt::Float64 = 0.0001
    dampening::Float64 = 1 - (dt * 5)
    gravity::Float64 = -9.81

    limit::Float64 = 5
    println("num springs evaluated: ", size(springs)[1])
    println("time multiplier: ",  limit / dt)

    @time begin
        while t < limit
            adjust::Float64 = 1 + sin(t * kOscillationFrequency) * 0.1
            for l in springs
                p1::Point = points[l.p1]
                p2::Point = points[l.p2]
                xd::Float64 = p1.x - p2.x;
                yd::Float64 = p1.y - p2.y;
                zd::Float64 = p1.z - p2.z;
                dist::Float64 = sqrt(xd * xd + yd * yd + zd * zd)

                # negative if repelling, positive if attracting
                f::Float64 = l.k * (dist - (l.l0 * adjust))
                # distribute force across the axes
                fdist::Float64 = f / dist
                dx::Float64 = fdist * xd
                dy::Float64 = fdist * yd
                dz::Float64 = fdist * zd

                p1.fx -= dx
                p2.fx += dx

                p1.fy -= dy
                p2.fy += dy

                p1.fz -= dz
                p2.fz += dz
            end

            for p in points
                fy = p.fy
                fx = p.fx
                fz = p.fz
                mass = p.mass

                if p.y < 0
                    fy += -kGround * p.y
                    fh::Float64 = sqrt(fx * fx + fz * fz)
                    if fh < abs(fy * staticFriction)
                        fx = 0
                        fz = 0
                    else
                        fyfric = fy * kineticFriction
                        fx = fx - fyfric
                        fz = fz - fyfric
                    end
                end
                ax = fx / mass
                ay = fy / mass + gravity
                az = fz / mass
                # reset the force cache
                p.fx = 0
                p.fy = 0
                p.fz = 0
                p.vx = (ax * dt + p.vx) * dampening
                p.vy = (ay * dt + p.vy) * dampening
                p.vz = (az * dt + p.vz) * dampening
                p.x += p.vx * dt
                p.y += p.vy * dt
                p.z += p.vz * dt
            end
            t += dt
        end
    end
    for p in points
        println(p.x, ",", p.y, ",", p.z)
    end
end

main()
#@profiler main()

max_stars_repo_name: unionfinder/PrimordialOozebot

max_issues_repo_name: unionfinder/PrimordialOozebot

max_forks_repo_name: unionfinder/PrimordialOozebot


text:
export SampleListFormConstraint, LeftProposal, RightProposal

using Random

struct LeftProposal end
struct RightProposal end

"""
    SampleListFormConstraint

One of the form constraint objects. Approximates `DistProduct` with a SampleList object.

See also: [`constrain_form`](@ref), [`DistProduct`](@ref)
"""
struct SampleListFormConstraint{N, R, S, M} <: AbstractFormConstraint
    rng      :: R
    strategy :: S
    method   :: M
end

SampleListFormConstraint(nsamples::Int, strategy::S, method::M = BootstrapImportanceSampling())         where { S, M }                   = SampleListFormConstraint(Random.GLOBAL_RNG, nsamples, strategy, method)
SampleListFormConstraint(rng::R, nsamples::Int, strategy::S, method::M = BootstrapImportanceSampling()) where { R <: AbstractRNG, S, M } = SampleListFormConstraint{nsamples, R, S, M}(rng, strategy, method)

default_form_check_strategy(::SampleListFormConstraint) = FormConstraintCheckLast()

is_point_mass_form_constraint(::SampleListFormConstraint) = false

constrain_form(::SampleListFormConstraint, something) = something

__approximate(constraint::SampleListFormConstraint{N, R, S, M}, left, right) where { N, R, S <: LeftProposal, M }  = approximate_prod_with_sample_list(constraint.rng, constraint.method, left, right, N)
__approximate(constraint::SampleListFormConstraint{N, R, S, M}, left, right) where { N, R, S <: RightProposal, M } = approximate_prod_with_sample_list(constraint.rng, constraint.method, right, left, N)

function constrain_form(constraint::SampleListFormConstraint, something::Message{ <: DistProduct })
    product = getdata(something)
    left    = constrain_form(constraint, getleft(product))
    right   = constrain_form(constraint, getright(product))
    return Message(__approximate(constraint, left, right), is_clamped(something), is_initial(something))
end


text:
export flattened_tup_chain

flattened_tup_chain(::Type{NamedTuple{(), Tuple{}}}; prefix = (Symbol(),)) = ()
flattened_tup_chain(::Type{T}; prefix = (Symbol(),)) where {T <: Real} =
    (prefix,)
flattened_tup_chain(::Type{T}; prefix = (Symbol(),)) where {T <: SArray} =
    (prefix,)
flattened_tup_chain(
    ::Type{T};
    prefix = (Symbol(),),
) where {T <: SHermitianCompact} = (prefix,)

"""
    flattened_tup_chain(::Type{T}) where {T <: Union{NamedTuple,NTuple}}

An array of tuples, containing symbols
and integers for every combination of
each field in the `Vars` array.
"""
function flattened_tup_chain(
    ::Type{T};
    prefix = (Symbol(),),
) where {T <: Union{NamedTuple, NTuple}}
    map(1:fieldcount(T)) do i
        Ti = fieldtype(T, i)
        name = fieldname(T, i)
        sname = name isa Int ? name : Symbol(name)
        flattened_tup_chain(
            Ti,
            prefix = prefix == (Symbol(),) ? (sname,) : (prefix..., sname),
        )
    end |>
    Iterators.flatten |>
    collect
end

max_forks_repo_name: jm-c/CLIMA


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ baff755d-1b7f-42ad-b7c4-8e8098ceadf2
# we write a function that solves Ly = b for y where L is lower triangular and b is a vector
function Lsolve(L, b)

    # we get the sizes of L and b
    LRows, LCols = size(L)
    bRows = size(b)[1]

    # we check that L and b ahve appropriate sizes
    if LRows != LCols || bRows != LRows
        error("The size of L or b is not appropriate")
    end

    # we initialize y to b
    y = copy(b);

    # we loop over the rows of L
    for k = 1:LRows

        # we check that we dont divide by zero
        if L[k,k] == 0
            error("There are zeros on the diagonal of L")
        end

        # we modify y accordingly
        for j = 1:k-1
            y[k] -= L[k,j] * y[j];
        end
        y[k] /= L[k,k];
    end

    # we return y
    return y
end

# ╔═╡ b62739ea-a4df-4516-88a5-59a159053a26
md"Solve $Ly = b$"

# ╔═╡ 15374d6a-ba4a-402b-acd3-3322e6fb189a
L = [1.0 0.0; 2.0 5.0];

# ╔═╡ 48329efb-5c10-4ada-8baf-b5b69a97fc1f
b = [2.0; 3.0];

# ╔═╡ a8948867-db61-4b3d-9e65-f9213783e1c1
y = Lsolve(L, b)

# ╔═╡ 1f47da5b-668b-4471-88c6-0cdfd623e46d
md"Check Lsolve, this should eval to zero:"

# ╔═╡ 7d9856ad-b34f-4277-a7be-a86428c1c56e
L*y-b

# ╔═╡ 3c84b49b-d112-4e8a-ac59-b1961ff2ac1c
md"Test case: b has the wrong size"

# ╔═╡ 4c6ef08d-f462-424f-93b9-d9e4d1425422
bWrong = [2.0; 3.0; 4.0];

# ╔═╡ e1ae3cd5-f6e1-4d63-82ef-e5cf66fc7f64
Lsolve(L, bWrong)

# ╔═╡ eba7f592-fb39-40c6-80b4-931fab81ec96
md"Test case: L is singular"

# ╔═╡ b4a3f8e3-c332-4977-b2c2-e5874d88ee68
Lwrong = [1.0 2.0; 2.0 5.0];

# ╔═╡ 241572c6-3cdc-4147-b2e0-8a6364bc0545
yWrong = Lsolve(Lwrong, b);

# ╔═╡ a7c3d19d-7207-4f24-8c60-1c7d93e2a8b8
Lwrong*yWrong - b

# ╔═╡ 6f6d27f4-c3cd-4ebe-a69a-d522aa5a9503
md"Should be zero but is not... It only looks at the lower triangle of L."

# ╔═╡ 307b2b77-aca1-418b-9084-4816a02fdcf1
# we write a function that solves Ux = y for x where U is upper triangular and y is a vector
function Usolve(U, y)

    # we get the sizes of U and y
    URows, UCols = size(U)
    yRows = size(y)[1]

    # we check that U and y ahve appropriate sizes
    if URows != UCols || yRows != URows
        error("The size of U or y is not appropriate")
    end

    # we initialize x to y
    x = copy(y);

    # we loop backwards over the rows of U
    for k = URows:-1:1

        # we check that we don't divide by zero
        if U[k,k] == 0
            error("There are zeros on the diagonal of U")
        end

        # we modify x accordingly
        for j = k+1:URows
            x[k] -= U[k,j] * x[j];
        end
        x[k] /= U[k,k];
    end

    # we return x
    return x
end

# ╔═╡ 499dd3c6-ab18-47c1-a38f-eaee61e2371f
U = [2.0 1.0; 0.0 3.0];

# ╔═╡ 9aa0e987-a432-431c-9200-9a902de25726
md"solve $Ux=y$, this is now the solution to $LUx = b$."

# ╔═╡ f70a999b-a08f-4622-b018-140ef0d0a3e3
x = Usolve(U, y)

# ╔═╡ a357e6fa-691d-4d49-88b0-f42e6a02c886
md"Check Usolve, this should both eval to zero (up to roundoff):"

# ╔═╡ 0116a06e-1771-4c5f-a73b-174524b25e15
L*U*x - b

# ╔═╡ 5e4a3e9c-659d-46d0-8306-705f9231c0f7
# we write a function that decomposes a square matrix A as A = LU
# where L is lower triangular and U is upper triangular
function LUdecomp(A)

    # we get the size of A
    n, m = size(A)

    # we check if A is indeed a square matrix
    if n != m
        error("Input must be a square matrix")
    end

    # we initialize the output matrices
    L = zeros(n, n)
    U = zeros(n, n)

    # we make a copy of A in order not to modify the original
    B = copy(A)

    for k = 1:n

        # we set the entries of the kth row of U
        for j = k:n
            U[k,j] = B[k,j]
        end

        # we check that we don't divide by zero
        if U[k,k] == 0
            error("** A^(k-1)_{k,k} == 0 in LU decomp")
        end

        # we set the entries of the kth column of L
        for i = k:n
           L[i,k] = B[i,k] / U[k,k]
        end

        # we modify A for the next iteration
        for i = k:n
            for j = k:n
                B[i,j] -= L[i,k] * U[k,j]
            end
        end
    end

    # we return L and U
    return (L, U)
end

# ╔═╡ a5248f09-c39c-4a6a-a327-866ae458b8d0
A = rand(5,5);

# ╔═╡ 3060c850-52ed-499b-b29a-1b00d7276d33
LA, UA = LUdecomp(A);

# ╔═╡ 9ffbfac3-21bc-4a19-a982-d58882a860cc
md"Check LUdecomp, this should eval to zero (up to roundoff):"

# ╔═╡ 040b92ed-bfd9-4a9c-b8e8-82fc955af5a3
LA*UA-A

# ╔═╡ 71e7a466-eea0-4b50-8f16-db774235e704
#we write a function that solves Ax = b for x by LU decomposition
function Asolve(A, b)
	L, U = LUdecomp(A)
	y = Lsolve(L, b)
	x = Usolve(U, y)
	return x
end

# ╔═╡ 1d575359-23c0-4236-bd3f-937400fc57ce
b2 = rand(5);

# ╔═╡ d30c113e-45c6-488d-83ff-c2d557e5f776
x2 = Asolve(A,b2);

# ╔═╡ bf6681d4-03b3-4eaa-99c1-236bf8daaa54
md"Check Asolve, this should eval to zero (up to roundoff):"

# ╔═╡ cc526a1c-8f34-4032-b5ae-548b67cac8bb
 A*x2 - b2

# ╔═╡ 20be458b-7b6c-4aff-822a-71969b013c53
md"This will throw an error..."

# ╔═╡ 5e8286c6-268a-47a2-8b4e-e026b42b5671
C = [0.0 1.0; 1.0 0.0];

# ╔═╡ af0082b7-7ee5-4e7f-9492-7e221a1328ab
LC, UC = LUdecomp(C)

# ╔═╡ 9488d28a-b76c-4863-80da-fcefd5ab2d4d
md"Let's try it with a large matrix."

# ╔═╡ d302f22e-c368-4e17-a31e-abd43db74e30
D = rand(800, 800);

# ╔═╡ 5bb781c8-7af0-4ccc-bb44-fe279127c44d
LD, UD = LUdecomp(D);

# ╔═╡ 1d599c5b-ff5b-42b9-bf79-7b67f5dd0d94
md"Largest divaition from zero in $L_DU_D-D$ is basically zero."

# ╔═╡ db2526a8-1a4a-4c00-a26d-8e55f55f1923
max(abs.(LD*UD-D)...)

# ╔═╡ Cell order:
# ╠═baff755d-1b7f-42ad-b7c4-8e8098ceadf2
# ╟─b62739ea-a4df-4516-88a5-59a159053a26
# ╠═15374d6a-ba4a-402b-acd3-3322e6fb189a
# ╠═48329efb-5c10-4ada-8baf-b5b69a97fc1f
# ╠═a8948867-db61-4b3d-9e65-f9213783e1c1
# ╟─1f47da5b-668b-4471-88c6-0cdfd623e46d
# ╠═7d9856ad-b34f-4277-a7be-a86428c1c56e
# ╟─3c84b49b-d112-4e8a-ac59-b1961ff2ac1c
# ╠═4c6ef08d-f462-424f-93b9-d9e4d1425422
# ╠═e1ae3cd5-f6e1-4d63-82ef-e5cf66fc7f64
# ╟─eba7f592-fb39-40c6-80b4-931fab81ec96
# ╠═b4a3f8e3-c332-4977-b2c2-e5874d88ee68
# ╠═241572c6-3cdc-4147-b2e0-8a6364bc0545
# ╠═a7c3d19d-7207-4f24-8c60-1c7d93e2a8b8
# ╟─6f6d27f4-c3cd-4ebe-a69a-d522aa5a9503
# ╠═307b2b77-aca1-418b-9084-4816a02fdcf1
# ╠═499dd3c6-ab18-47c1-a38f-eaee61e2371f
# ╟─9aa0e987-a432-431c-9200-9a902de25726
# ╠═f70a999b-a08f-4622-b018-140ef0d0a3e3
# ╟─a357e6fa-691d-4d49-88b0-f42e6a02c886
# ╠═0116a06e-1771-4c5f-a73b-174524b25e15
# ╠═5e4a3e9c-659d-46d0-8306-705f9231c0f7
# ╠═a5248f09-c39c-4a6a-a327-866ae458b8d0
# ╠═3060c850-52ed-499b-b29a-1b00d7276d33
# ╟─9ffbfac3-21bc-4a19-a982-d58882a860cc
# ╠═040b92ed-bfd9-4a9c-b8e8-82fc955af5a3
# ╠═71e7a466-eea0-4b50-8f16-db774235e704
# ╠═1d575359-23c0-4236-bd3f-937400fc57ce
# ╠═d30c113e-45c6-488d-83ff-c2d557e5f776
# ╟─bf6681d4-03b3-4eaa-99c1-236bf8daaa54
# ╠═cc526a1c-8f34-4032-b5ae-548b67cac8bb
# ╟─20be458b-7b6c-4aff-822a-71969b013c53
# ╠═5e8286c6-268a-47a2-8b4e-e026b42b5671
# ╠═af0082b7-7ee5-4e7f-9492-7e221a1328ab
# ╟─9488d28a-b76c-4863-80da-fcefd5ab2d4d
# ╠═d302f22e-c368-4e17-a31e-abd43db74e30
# ╠═5bb781c8-7af0-4ccc-bb44-fe279127c44d
# ╟─1d599c5b-ff5b-42b9-bf79-7b67f5dd0d94
# ╠═db2526a8-1a4a-4c00-a26d-8e55f55f1923

max_stars_repo_path: 1a/LU Decomposition/LUtest.jl

max_issues_repo_path: 1a/LU Decomposition/LUtest.jl

max_forks_repo_path: 1a/LU Decomposition/LUtest.jl


text:
planets = ["earth"=>["mu"=>398600.448073446, "j2"=>0.0010826269, "re"=>6378.1363,
    "rp"=>6356.7523, "f"=>0.003352813178, "omega"=>7.29211514670698e-05]]


text:
function jump_max_area_perimeter_Bn(n)
    if mod(log2(n),1)==0 && n>=8
        d = pi/n - asin(sin(2*pi/n)/2);
        Ll = 2*n*sin(pi/(2*n))*cos(d/2);
        Lu = 2*n*sin(pi/(2*n));
        model = Model(Ipopt.Optimizer)
        #model = Model(() -> AmplNLWriter.Optimizer(Couenne_jll.amplexe))
        # Variables
        @variable(model, a[i = 1:Int(n/4+1)] >= 0, start = pi/n+(-1)^(i-1)*d)
        # Perimeter
        @NLobjective(model, Max, 4*sin(a[1]/2)+sum(8*sin(a[i]/2) for i in 2:Int(n/4))+4*sin(a[Int(n/4+1)]/2))
        # Optimal diameter graph
        @constraint(model, a[1] + sum(2*a[k] for k in 2:Int(n/4)) + a[Int(n/4+1)] == pi/2)
        @NLconstraint(model, sin(a[1])+sum((-1)^(i-1)*sin(a[1]+
        sum(2*a[j] for j in 2:i)) for i in 2:Int(n/4)) == -1/2)
        @constraint(model, [i = 1:Int(n/4)], a[i] <= pi/6)
        @constraint(model, a[Int(n/4+1)] <= pi/3)
        # Bounds
        @NLconstraint(model, Ll <= 4*sin(a[1]/2)+sum(8*sin(a[i]/2) for i in 2:Int(n/4))+4*sin(a[Int(n/4+1)]/2) <= Lu)
        # Solving
        optimize!(model)
        L = objective_value(model)
        a = value.(a)
        ct = solve_time(model)
        return(L,a,ct)
    end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
using FromFile
@from "../subsubchain1.jl" import a


text:
"""
$(TYPEDEF)

A system of ordinary differential equations.

# Fields
$(FIELDS)

# Example

```
using ModelingToolkit

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

de = ODESystem(eqs,t,[x,y,z],[σ,ρ,β])
```
"""
struct ODESystem <: AbstractODESystem
    """The ODEs defining the system."""
    eqs::Vector{Equation}
    """Independent variable."""
    iv::Variable
    """Dependent (state) variables."""
    states::Vector{Variable}
    """Parameter variables."""
    ps::Vector{Variable}
    """
    Time-derivative matrix. Note: this field will not be defined until
    [`calculate_tgrad`](@ref) is called on the system.
    """
    tgrad::RefValue{Vector{Expression}}
    """
    Jacobian matrix. Note: this field will not be defined until
    [`calculate_jacobian`](@ref) is called on the system.
    """
    jac::RefValue{Matrix{Expression}}
    """
    `Wfact` matrix. Note: this field will not be defined until
    [`generate_factorized_W`](@ref) is called on the system.
    """
    Wfact::RefValue{Matrix{Expression}}
    """
    `Wfact_t` matrix. Note: this field will not be defined until
    [`generate_factorized_W`](@ref) is called on the system.
    """
    Wfact_t::RefValue{Matrix{Expression}}
    """
    Name: the name of the system
    """
    name::Symbol
    """
    systems: The internal systems
    """
    systems::Vector{ODESystem}
end

function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
                   systems = ODESystem[],
                   name=gensym(:ODESystem))
    iv′ = convert(Variable,iv)
    dvs′ = convert.(Variable,dvs)
    ps′ = convert.(Variable,ps)
    tgrad = RefValue(Vector{Expression}(undef, 0))
    jac = RefValue(Matrix{Expression}(undef, 0, 0))
    Wfact   = RefValue(Matrix{Expression}(undef, 0, 0))
    Wfact_t = RefValue(Matrix{Expression}(undef, 0, 0))
    ODESystem(deqs, iv′, dvs′, ps′, tgrad, jac, Wfact, Wfact_t, name, systems)
end

var_from_nested_derivative(x::Constant) = (missing, missing)
var_from_nested_derivative(x,i=0) = x.op isa Differential ? var_from_nested_derivative(x.args[1],i+1) : (x.op,i)

iv_from_nested_derivative(x) = x.op isa Differential ? iv_from_nested_derivative(x.args[1]) : x.args[1].op
iv_from_nested_derivative(x::Constant) = missing

function ODESystem(eqs; kwargs...)
    # NOTE: this assumes that the order of algebric equations doesn't matter
    diffvars = OrderedSet{Variable}()
    allstates = OrderedSet{Variable}()
    ps = OrderedSet{Variable}()
    # reorder equations such that it is in the form of `diffeq, algeeq`
    diffeq = Equation[]
    algeeq = Equation[]
    # initial loop for finding `iv`
    iv = nothing
    for eq in eqs
        if !(eq.lhs isa Constant) # assume eq.lhs is either Differential or Constant
            iv = iv_from_nested_derivative(eq.lhs)
        end
    end
    iv === nothing && throw(ArgumentError("No differential variable detected."))
    for eq in eqs
        for var in vars(eq.rhs for eq ∈ eqs)
            var isa Variable || continue
            if isparameter(var)
                isequal(var, iv) || push!(ps, var)
            else
                push!(allstates, var)
            end
        end
        if eq.lhs isa Constant
            push!(algeeq, eq)
        else
            diffvar = first(var_from_nested_derivative(eq.lhs))
            iv == iv_from_nested_derivative(eq.lhs) || throw(ArgumentError("An ODESystem can only have one independent variable."))
            diffvar in diffvars && throw(ArgumentError("The differential variable $diffvar is not unique in the system of equations."))
            push!(diffvars, diffvar)
            push!(diffeq, eq)
        end
    end
    algevars = setdiff(allstates, diffvars)
    # the orders here are very important!
    return ODESystem(append!(diffeq, algeeq), iv, vcat(collect(diffvars), collect(algevars)), ps; kwargs...)
end

Base.:(==)(sys1::ODESystem, sys2::ODESystem) =
    _eq_unordered(sys1.eqs, sys2.eqs) && isequal(sys1.iv, sys2.iv) &&
    _eq_unordered(sys1.states, sys2.states) && _eq_unordered(sys1.ps, sys2.ps)
# NOTE: equality does not check cached Jacobian

function rename(sys::ODESystem,name)
    ODESystem(sys.eqs, sys.iv, sys.states, sys.ps, sys.tgrad, sys.jac, sys.Wfact, sys.Wfact_t, name, sys.systems)
end


text:
module YaDeepLearning

using Pkg
Pkg.activate("MLJ_env", shared=true)

using Random, Distributions, Flux.Data.MNIST, Images

export TensorDT
export AbstractLayers
export AbstractOptimizers
export Layers
export Optimizers
export Activations
export Initializations
export Loss
export AA

include("./tensor_dt.jl")
import .TensorDT: Tensor

include("./abstract_layers.jl")
import .AbstractLayers

include("./abstract_optimizers.jl")
import .AbstractOptimizers

include("./optimizers.jl")
import .Optimizers

include("./initializations.jl")
import .Initializations

include("./layers.jl")
import .Layers

include("./activations.jl")
import .Activations

include("./loss.jl")
import .Loss

## ======================================================================
## Const
## ======================================================================
const AA = AbstractArray

end ## Module

max_stars_repo_path: Data Science From Scratch/src/YaDeepLearning.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
#init state
function initstate(domain::Domain, problem::GenericProblem)
    return initstate(domain, problem.objtypes, problem.continuous_inits, problem.init)
end

function initstate(domain::Domain, objtypes::AbstractDict)
    types = Set{Term}([Compound(ty, Term[o]) for (o, ty) in objtypes])
return GenericState(types, Set{Term}(), Dict{Symbol,Any}())
end

function initstate(domain::Domain, objtypes::AbstractDict, cont_inits::Compound,
    fluents::AbstractVector)
    state = initstate(domain, objtypes)
    for t in fluents
        if t.name == :(==) # Non-Boolean fluents
            @assert length(t.args) == 2 "Assignments must have two arguments."
            term, val = t.args[1], t.args[2]
            @assert !isa(term, Var) "Initial terms cannot be unbound variables."
            state[term] = evaluate(domain, state, val)
        else # Boolean fluents
            push!(state.facts, t)
        end
    end
    for var in cont_inits.args
        exp = Meta.parse(var.name)
        state.values[exp.args[2]] = exp.args[3]
    end
    return state
end

function initstate(domain::Domain, objtypes::AbstractDict,
    fluents::AbstractDict)
    state = initstate(domain, objtypes)
    for (name, val) in fluents
        set_fluent!(state, val, name)
    end
    return state
end

#goal state
function goalstate(domain::Domain, problem::GenericProblem)
    types = Set{Term}([Compound(ty, Term[o]) for (o, ty) in problem.objtypes])
    facts = Set{Term}(flatten_conjs(problem.goal))
    cont_inits = problem.continuous_goals
    state = GenericState(types, facts, Dict{Symbol,Any}())
    for var in cont_inits.args
        exp = Meta.parse(var.name)
        state.values[exp.args[2]] = exp.args[3]
    end
    return state
end

function goalstate(domain::Domain, objtypes::AbstractDict, terms)
    types = Set{Term}([Compound(ty, Term[o]) for (o, ty) in objtypes])
    goal = GenericState(types, Set{Term}(), Dict{Symbol,Any}())
    for t in flatten_conjs(terms)
        if t.name == :(==) # Function equality
            @assert length(t.args) == 2 "Assignments must have two arguments."
            term, val = t.args[1], t.args[2]
            @assert(is_func(term, domain) && !is_attached_func(term, domain),
                "Unrecognized function $(term.name).")
            @assert(is_ground(term), "Assigned terms must be ground.")
            @assert(isa(val, Const), "Terms must be equal to constants.")
            goal[term] = val.name
        elseif is_pred(t, domain) # Predicates
            push!(goal.facts, t)
        else
            error("Term $t in $terms cannot be handled.")
        end
    end
    return goal
end


text:
using open62541

running = Ref{Bool}(true)
server = UA_Server_new()
UA_ServerConfig_setMinimalCustomBuffer(UA_Server_getConfig(server), 4840, C_NULL, 0, 0)
retval = UA_Server_run(server, running)