text:
import Base: copy, reset

abstract type DynamicsType end
abstract type Continuous <: DynamicsType end
abstract type Discrete <: DynamicsType end

abstract type ModelType end
abstract type Nominal <: ModelType end
abstract type Uncertain <: ModelType end

abstract type AbstractModel end
abstract type Model{M<:ModelType,D<:DynamicsType} <: AbstractModel end

export
    Nominal,
    Uncertain,
    AnalyticalModel,
    discretize_model,
    dynamics

"""
$(TYPEDEF)
Dynamics model
Holds all information required to uniquely describe a dynamic system, including
a general nonlinear dynamics function of the form `ẋ = f(x,u)`, where x ∈ ℜⁿ are
the states and u ∈ ℜᵐ are the controls.
Dynamics function, f, should be of the form
    f(ẋ,x,u,p) for Continuous models, where ẋ is the state derivative
    f(ẋ,x,u,p,dt) for Discrete models, where ẋ is the state at the next time step
    and x is the state vector, u is the control input vector, and p is an optional `NamedTuple` of static parameters (mass, gravity, etc.)
Dynamics jacobians, ∇f, should be of the form
    ∇f(Z,x,u,p) for Continuous models, and
    ∇f(Z,x,u,,p,dt) for discrete models
    where p is the same `NamedTuple` of parameters used in the dynamics
"""
struct AnalyticalModel{M,D} <: Model{M,D}
    f::Function   # dynamics f(ẋ,x,u)
    ∇f::Function  # dynamics jacobian
    n::Int        # number of states
    m::Int        # number of controls
    r::Int       # number of uncertain parameters
    params::NamedTuple
    evals::Vector{Int}
    info::Dict{Symbol,Any}

    """ $(SIGNATURES)
    Create a dynamics model given a dynamics function and Jacobian, with n states and m controls.
    Dynamics function should be of the form
        f(ẋ,x,u,p) for Continuous models, where ẋ is the state derivative
        f(ẋ,x,u,p,dt) for Discrete models, where ẋ is the state at the next time step
        and x is the state vector, u is the control input vector, and p is a `NamedTuple` of static parameters (mass, gravity, etc.)
    Optionally pass in a dictionary `d` with model information.
    `check_functions` option runs verification checks on the dynamics function and Jacobian to make sure they have the correct forms.
    """
    function AnalyticalModel{M,D}(f::Function, ∇f::Function, n::Int, m::Int,r::Int,
            p::NamedTuple=NamedTuple(), d::Dict{Symbol,Any}=Dict{Symbol,Any}();
            check_functions::Bool=false) where {M<:ModelType,D<:DynamicsType}
        d[:evals] = 0
        evals = [0;0]
        if check_functions
            # Make dynamics inplace
            if is_inplace_dynamics(f,n,m,r)
                f! = f
            else
                f! = wrap_inplace(f)
            end
            # ∇f! = _check_jacobian(D,f,∇f,n,m)
            new{M,D}(f!,∇f,n,m,r,p,evals,d)
        else
            new{M,D}(f,∇f,n,m,r,p,evals,d)
        end
    end
end

"Docstring for this guy"
function AnalyticalModel{M,D}(f::Function, n::Int, m::Int, r::Int, d::Dict{Symbol,Any}=Dict{Symbol,Any}()) where {M<:ModelType,D<:DynamicsType}
    p = NamedTuple()
    ∇f, = generate_jacobian(M,D,f,n,m)
    AnalyticalModel{M,D}(f,∇f,n,m,r,p,d)
end

function AnalyticalModel{Uncertain,D}(f::Function, n::Int, m::Int, r::Int, d::Dict{Symbol,Any}=Dict{Symbol,Any}()) where {D<:DynamicsType}
    p = NamedTuple()
    ∇f, = generate_jacobian(Uncertain,D,f,n,m,r)
    AnalyticalModel{Uncertain,D}(f,∇f,n,m,r,p,d)
end

function AnalyticalModel{M,D}(f::Function, n::Int, m::Int, r::Int, p::NamedTuple, d::Dict{Symbol,Any}=Dict{Symbol,Any}()) where {M<:ModelType,D<:DynamicsType}
    f_p(ẋ,x,u) = f(ẋ,x,u,p)
    f_p(ẋ,x,u,p) = f(ẋ,x,u,p)
    ∇f, = generate_jacobian(M,D,f_p,n,m)
    AnalyticalModel{M,D}(f_p,∇f,n,m,r,p,d)
end

""" $(TYPEDSIGNATURES)
Create a dynamics model, using ForwardDiff to generate the dynamics jacobian, with parameters
Dynamics function passes in parameters:
    f(ẋ,x,u,p)
    where p in NamedTuple of parameters
"""
Model(f::Function, n::Int, m::Int, d::Dict{Symbol,Any}=Dict{Symbol,Any}()) =
    AnalyticalModel{Nominal,Continuous}(f,n,m,0,d)

UncertainModel(f::Function, n::Int, m::Int, r::Int, d::Dict{Symbol,Any}=Dict{Symbol,Any}()) =
    AnalyticalModel{Uncertain,Continuous}(f,n,m,r,d)

""" $(TYPEDSIGNATURES)
Create a dynamics model, using ForwardDiff to generate the dynamics jacobian, without parameters
Dynamics function of the form:
    f(ẋ,x,u)
"""
Model(f::Function, n::Int, m::Int, p::NamedTuple, d::Dict{Symbol,Any}=Dict{Symbol,Any}()) =
    AnalyticalModel{Nominal,Continuous}(f,n,m,0,p,d)

UncertainModel(f::Function, n::Int, m::Int, r::Int, p::NamedTuple, d::Dict{Symbol,Any}=Dict{Symbol,Any}()) =
    AnalyticalModel{Uncertain,Continuous}(f,n,m,p,r,d)

""" $(TYPEDSIGNATURES)
Create a dynamics model with an analytical Jacobian, with parameters
Dynamics functions pass in parameters:
    f(ẋ,x,u,p)
    ∇f(Z,x,u,p)
    where p in NamedTuple of parameters
"""
Model(f::Function, ∇f::Function, n::Int, m::Int, p::NamedTuple, d::Dict{Symbol,Any}=Dict{Symbol,Any}()) = begin
    f_p(ẋ,x,u) = f(ẋ,x,u,p)
    f_p(ẋ,x,u,p) = f(ẋ,x,u,p)
    ∇f_p(Z,x,u) = ∇f(Z,x,u,p)
    AnalyticalModel{Nominal,Continuous}(f_p,∇f_p,n,m,0,p,d, check_functions=true); end

UncertainModel(f::Function, ∇f::Function, n::Int, m::Int, r::Int, p::NamedTuple, d::Dict{Symbol,Any}=Dict{Symbol,Any}()) = begin
    f_p(ẋ,x,u,w) = f(ẋ,x,u,w,p)
    f_p(ẋ,x,u,w,p) = f(ẋ,x,u,w,p)
    ∇f_p(Z,x,u,w) = ∇f(Z,x,u,w,p)
    AnalyticalModel{Uncertain,Continuous}(f_p,∇f_p,n,m,r,p,d, check_functions=true); end

""" $(TYPEDSIGNATURES)
Create a dynamics model with an analytical Jacobian, without parameters
Dynamics functions pass of the form:
    f(ẋ,x,u)
    ∇f(Z,x,u)
"""
Model(f::Function, ∇f::Function, n::Int, m::Int, d::Dict{Symbol,Any}=Dict{Symbol,Any}()) = begin
    p = NamedTuple()
    AnalyticalModel{Nominal,Continuous}(f,∇f,n,m,0,p,d, check_functions=true); end

UncertainModel(f::Function, ∇f::Function, n::Int, m::Int, r::Int, d::Dict{Symbol,Any}=Dict{Symbol,Any}()) = begin
    p = NamedTuple()
    AnalyticalModel{Uncertain,Continuous}(f,∇f,n,m,r,p,d, check_functions=true); end

""" ```julia
evaluate!(ẋ, model::Model{M,Continuous}, x, u)
```
Evaluate the continuous dynamics at state `x` and control `u`
Keeps track of the number of evaluations
"""
function evaluate!(ẋ::AbstractVector,model::Model{M,Continuous},x::AbstractVector,u::AbstractVector) where M <: ModelType
    model.f(ẋ,x,u)
    model.evals[1] += 1
end

""" ```julia
evaluate!(ẋ, model::Model{M,Discrete}, x, u, dt)
```
 Evaluate the discrete dynamics at state `x` and control `u` and time step `dt`
Keeps track of the number of evaluations
"""
function evaluate!(ẋ::AbstractVector,model::Model{M,Discrete},x::AbstractVector,u::AbstractVector,dt::T) where {M <: ModelType,T}
    model.f(ẋ,x,u,dt)
    model.evals[1] += 1
end

""" ```julia
jacobian!(Z, model::Model{M,Continuous}, x, u)
```
Evaluate the dynamics Jacobian simultaneously at state `x` and control `x`
Keeps track of the number of evaluations
"""
function jacobian!(Z::AbstractMatrix,model::Model{M,Continuous},x::AbstractVector,u::AbstractVector) where M <: ModelType
    model.∇f(Z,x,u)
    model.evals[2] += 1
end

""" ```julia
jacobian!(Z, model::Model{M,Discrete}, x, u, dt)
```
Evaluate the dynamics Jacobian simultaneously at state `x` and control `x`
Keeps track of the number of evaluations
"""
function jacobian!(Z::AbstractArray{T},model::Model{M,Discrete},x::AbstractVector,u::AbstractVector,dt::T) where {M <: ModelType,T <: AbstractFloat}
    model.∇f(Z,x,u,dt)
    model.evals[2] += 1
end

""" ```julia
jacobian!(Z::PartedVecTrajectory, model::Model{M,Discrete}, X, U, dt)
```
Evaluate discrete dynamics Jacobian along entire trajectory
"""
function jacobian!(Z::PartedMatTrajectory{T},model::Model{M,Discrete},X::VectorTrajectory{T},U::VectorTrajectory{T},dt::Vector{T}) where {M<:ModelType,T}
    N = length(X)
    for k = 1:N-1
        jacobian!(Z[k],model,X[k],U[k],dt[k])
    end
end

# Uncertain Dynamics
function evaluate!(ẋ::AbstractVector,model::Model{Uncertain,Continuous},x::AbstractVector,u::AbstractVector)
    model.f(view(ẋ,1:model.n),x[1:model.n],u[1:model.m],zeros(model.r))
    model.evals[1] += 1
end

function evaluate_uncertain!(ẋ::AbstractVector,model::Model{Uncertain,Continuous},x::AbstractVector,u::AbstractVector,w::AbstractVector)
    model.f(view(ẋ,1:model.n),x[1:model.n],u[1:model.m],w)
    model.evals[1] += 1
end

function evaluate!(ẋ::AbstractVector,model::Model{Uncertain,Discrete},x::AbstractVector,u::AbstractVector)
    # model.f(view(ẋ,1:model.n),x[1:model.n],u[1:model.m],zeros(model.r))
    model.f(ẋ,x[1:model.n],u[1:model.m],zeros(model.r))

    model.evals[1] += 1
end

function evaluate_uncertain!(ẋ::AbstractVector,model::Model{Uncertain,Discrete},x::AbstractVector,u::AbstractVector,w::AbstractVector)
    model.f(view(ẋ,1:model.n),x[1:model.n],u[1:model.m],w)
    model.evals[1] += 1
end

function evaluate!(ẋ::AbstractVector,model::Model{Uncertain,Discrete},x::AbstractVector,u::AbstractVector,dt::T) where T
    # model.f(view(ẋ,1:model.n),x[1:model.n],u[1:model.m],zeros(model.r),dt)
    model.f(ẋ,x[1:model.n],u[1:model.m],zeros(model.r),dt)
    model.evals[1] += 1
end

function evaluate_uncertain!(ẋ::AbstractVector,model::Model{Uncertain,Discrete},x::AbstractVector,u::AbstractVector,w::AbstractVector,dt::T) where T
    model.f(view(ẋ,1:model.n),x[1:model.n],u[1:model.m],w,dt)
    model.evals[1] += 1
end

function evaluate!(Z::AbstractMatrix,ẋ::AbstractVector,model::Model{Uncertain,Continuous},x::AbstractVector,u::AbstractVector)
    idx = NamedTuple{(:x,:u)}((1:model.n,1:model.m))
    model.∇f(Z,ẋ[idx.x],x[idx.x],u[idx.u],zeros(model.r))
    model.evals[1] += 1
    model.evals[2] += 1
end

function evaluate_uncertain!(Z::AbstractMatrix,ẋ::AbstractVector,model::Model{Uncertain,Continuous},x::AbstractVector,u::AbstractVector,w::AbstractVector)
    idx = NamedTuple{(:x,:u)}((1:model.n,1:model.m))
    model.∇f(Z,ẋ[idx.x],x[idx.x],u[idx.u],w)
    model.evals[1] += 1
    model.evals[2] += 1
end

function evaluate!(Z::AbstractMatrix,ẋ::AbstractVector,model::Model{Uncertain,Discrete},x::AbstractVector,u::AbstractVector,dt::T) where {M <: ModelType,T}
    idx = NamedTuple{(:x,:u)}((1:model.n,1:model.m))
    model.∇f(Z,ẋ[idx.x],x[idx.x],u[idx.u],zeros(model.r),dt)
    model.evals[1] += 1
    model.evals[2] += 1
end

function evaluate_uncertain!(Z::AbstractMatrix,ẋ::AbstractVector,model::Model{Uncertain,Discrete},x::AbstractVector,u::AbstractVector,dt::T) where {M <: ModelType,T}
    idx = NamedTuple{(:x,:u)}((1:model.n,1:model.m))
    model.∇f(Z,ẋ[idx.x],x[idx.x],u[idx.u],w,dt)
    model.evals[1] += 1
    model.evals[2] += 1
end
jacobian!(Z::AbstractArray,ẋ::AbstractVector,model::Model{Uncertain,Continuous},x::AbstractVector,u::AbstractVector) = evaluate!(Z,ẋ,model,x,u,zeros(model.r))
jacobian_uncertain!(Z::AbstractArray,ẋ::AbstractVector,model::Model{Uncertain,Continuous},x::AbstractVector,u::AbstractVector,w::AbstractVector) = evaluate!(Z,ẋ,model,x,u,w)
jacobian!(Z::AbstractArray,ẋ::AbstractVector,model::Model{Uncertain,Discrete},x::AbstractVector,u::AbstractVector,dt::T) where T = evaluate!(Z,ẋ,model,x,u,zeros(model.r),dt)
jacobian_uncertain!(Z::AbstractArray,ẋ::AbstractVector,model::Model{Uncertain,Discrete},x::AbstractVector,u::AbstractVector,w::AbstractVector,dt::T) where T = evaluate!(Z,ẋ,model,x,u,w,dt)

function jacobian!(Z::AbstractMatrix,model::Model{Uncertain,Continuous},x::AbstractVector,u::AbstractVector)
    model.∇f(Z,x,u,zeros(model.r))
    model.evals[2] += 1
end

function jacobian_uncertain!(Z::AbstractMatrix,model::Model{Uncertain,Continuous},x::AbstractVector,u::AbstractVector,w::AbstractVector)
    model.∇f(Z,x,u,w)
    model.evals[2] += 1
end

function jacobian!(Z::AbstractArray{T},model::Model{Uncertain,Discrete},x::AbstractVector{T},u::AbstractVector{T},dt::T) where T <: AbstractFloat
    model.∇f(Z,x,u,zeros(model.r),dt)
    model.evals[2] += 1
end

function jacobian_uncertain!(Z::AbstractArray,model::Model{Uncertain,Discrete},x::AbstractVector,u::AbstractVector,w::AbstractVector,dt::T) where T
    model.∇f(Z,x,u,w,dt)
    model.evals[2] += 1
end

function jacobian!(Z::PartedMatTrajectory{T},model::Model{M,Discrete},X::VectorTrajectory{T},U::VectorTrajectory{T},dt) where {M<:ModelType,T}
    N = length(X)
    for k = 1:N-1
        jacobian!(Z[k],model,X[k],U[k],dt[k])
    end
end

function jacobian!(Z::PartedMatTrajectory{T},model::Model{Uncertain,Discrete},X::VectorTrajectory{T},U::VectorTrajectory{T},dt) where T
    N = length(X)
    n = model.n; m = model.m; r = model.r
    m_e = length(U[1]) - length(U[2])
    idx = [(1:n)...,(n .+ (1:m))...,((n+m+m_e) .+ (1:r))...,(n+m+m_e+r+1)]
    jacobian!(view(Z[1],1:n,idx),model,X[1],U[1][1:m],dt[1])
    for k = 2:N-1
        jacobian!(Z[k],model,X[k],U[k],dt[k])
    end
end

function jacobian_uncertain!(Z::PartedMatTrajectory{T},model::Model{Uncertain,Discrete},X::VectorTrajectory{T},U::VectorTrajectory{T},W::VectorTrajectory{T},dt) where T
    N = length(X)
    for k = 1:N-1
        jacobian!(Z[k],model,X[k],U[k],W[k],dt[k])
    end
end

""" $(SIGNATURES) Return the number of dynamics evaluations """
evals(model::Model)::Int = model.evals[1]

""" $(SIGNATURES) Reset the evaluation counts for the model """
reset(model::Model) = begin model.evals[1] = 0; return nothing end

Base.length(model::Model{M,Discrete}) where M = model.n + model.m + 1
Base.length(model::Model{M,Continuous}) where M = model.n + model.m

Base.length(model::Model{Uncertain,Discrete}) = model.n + model.m + model.r + 1
Base.length(model::Model{Uncertain,Continuous}) = model.n + model.m + model.r

PartedArrays.create_partition(model::Model{M,Discrete}) where M = create_partition((model.n,model.m,1),(:x,:u,:dt))
PartedArrays.create_partition(model::Model{Uncertain,Discrete}) = create_partition((model.n,model.m,model.r,1),(:x,:u,:w,:dt))

PartedArrays.create_partition2(model::Model{M,Discrete}) where M = create_partition2((model.n,),(model.n,model.m,1),Val((:xx,:xu,:xdt)))
PartedArrays.create_partition2(model::Model{Uncertain,Discrete}) = create_partition2((model.n,),(model.n,model.m,model.r,1),Val((:xx,:xu,:xw,:xdt)))

PartedArrays.create_partition(model::Model{M,Continuous}) where M = create_partition((model.n,model.m),(:x,:u))
PartedArrays.create_partition(model::Model{Uncertain,Continuous}) = create_partition((model.n,model.m,model.r),(:x,:u,:w))

PartedArrays.create_partition2(model::Model{M,Continuous}) where M= create_partition2((model.n,),(model.n,model.m),Val((:xx,:xu)))
PartedArrays.create_partition2(model::Model{Uncertain,Continuous}) = create_partition2((model.n,),(model.n,model.m,model.r),Val((:xx,:xu,:xw)))

PartedArrays.PartedVector(model::Model) = PartedArray(zeros(length(model)),create_partition(model))
PartedArrays.PartedVector(T::Type,model::Model) = PartedArray(zeros(T,length(model)),create_partition(model))
PartedArrays.PartedMatrix(model::Model) = PartedArray(zeros(model.n,length(model)),create_partition2(model))
PartedArrays.PartedMatrix(T::Type,model::Model) = PartedArray(zeros(T,model.n,length(model)),create_partition2(model))

function dynamics(model::Model{M,D},xdot,x,u) where {M,D}
    model.f(xdot,x,u)
    model.evals[1] += 1
end

#TODO
function dynamics(model::Model{Uncertain,D},xdot::AbstractVector,x::AbstractVector,u::AbstractVector) where D <: DynamicsType
    model.f(xdot,x,u,zeros(model.r))
    model.evals[1] += 1
end

function dynamics(model::Model{Uncertain,D},xdot::AbstractVector,x::AbstractVector,u::AbstractVector,w::AbstractVector) where D <: DynamicsType
    model.f(xdot,x,u,w)
    model.evals[1] += 1
end

"""
$(TYPEDEF)
RigidBodyDynamics model. Wrapper for a RigidBodyDynamics Mechanism
"""
struct RBDModel{M,D} <: Model{M,D}
    f::Function # continuous dynamics (ie, differential equation)
    ∇f::Function  # dynamics jacobian
    n::Int # number of states
    m::Int # number of controls
    r::Int # number of uncertain parameters
    params::NamedTuple
    mech::Mechanism  # RigidBodyDynamics Mechanism
    evals::Vector{Int}
    info::Dict{Symbol,Any}
end

"""
$(SIGNATURES)
Model(mech::Mechanism, torques::Array{Bool, 1}) Constructor for an underactuated mechanism, where torques is a binary array
that specifies whether a joint is actuated.
"""
function Model(mech::Mechanism, torques::Array)

    # construct a model using robot dynamics equation assembed from URDF file
    n = num_positions(mech) + num_velocities(mech) + num_additional_states(mech)
    num_joints = length(joints(mech))  # subtract off joint to world

    if length(torques) != num_joints
        error("Torque underactuation specified does not match mechanism dimensions")
    end

    m = convert(Int,sum(torques)) # number of actuated (ie, controllable) joints
    torque_matrix = 1.0*Matrix(I,num_joints,num_joints)[:,torques.== 1] # matrix to convert from control inputs to mechanism joints

    statecache = StateCache(mech)
    dynamicsresultscache = DynamicsResultCache(mech)

    function f(ẋ::AbstractVector{T},x::AbstractVector{T},u::AbstractVector{T}) where T
        state = statecache[T]
        dyn = dynamicsresultscache[T]
        dynamics!(view(ẋ,1:n), dyn, state, x, torque_matrix*u)
        return nothing
    end

    f_wrap(ẋ,z) = f(ẋ,z[1:n],z[n .+ (1:m)])
    _∇f(Z,z) = ForwardDiff.jacobian!(Z,f_wrap,zeros(n),z)
    z0 = zeros(n+m)

    function ∇f(Z::AbstractArray{T},x::AbstractVector{T},u::AbstractVector{T}) where T
        z0[1:n] = x
        z0[n .+ (1:m)] = u
        _∇f(Z,z0)
    end

    d = Dict{Symbol,Any}()

    evals = [0;0]
    RBDModel{Nominal,Continuous}(f, ∇f, n, m, 0, NamedTuple(), mech, evals, d)
end

"""
$(SIGNATURES)
 Construct model from a `Mechanism` type from `RigidBodyDynamics`
 """
function Model(mech::Mechanism)
    m = length(joints(mech))  # subtract off joint to world
    Model(mech,ones(m))
end

"""$(SIGNATURES) Construct a fully actuated model from a string to a urdf file"""
function Model(urdf::String)
    # construct model using string to urdf file
    mech = parse_urdf(urdf)
    Model(mech)
end

"""$(SIGNATURES) Construct a partially actuated model from a string to a urdf file, where torques is a binary array that specifies whether a joint is actuated."""
function Model(urdf::String,torques::Array{Float64,1})
    # underactuated system (potentially)
    mech = parse_urdf(urdf)
    Model(mech,torques)
end

generate_jacobian(f!::Function,n::Int,m::Int,p::Int=n) = generate_jacobian(Nominal,Continuous,f!,n,m,p)

function generate_jacobian(::Type{Nominal},::Type{Continuous},f!::Function,n::Int,m::Int,p::Int=n)
    inds = (x=1:n,u=n .+ (1:m), px=(1:p,1:n),pu=(1:p,n .+ (1:m)))
    Z = PartedArray(zeros(p,n+m),inds)
    z = zeros(n+m)
    v0 = zeros(p)
    f_aug(dZ::AbstractVector,z::AbstractVector) = f!(dZ,view(z,inds.x), view(z,inds.u))
    ∇fz(Z::AbstractMatrix,v::AbstractVector,z::AbstractVector) = ForwardDiff.jacobian!(Z,f_aug,v,z)

    ∇f!(Z::AbstractMatrix,v::AbstractVector,x::AbstractVector,u::AbstractVector) = begin
        z[inds.x] = x
        z[inds.u] = u
        ∇fz(Z,v,z)
        return nothing
    end
    ∇f!(Z::AbstractMatrix,x::AbstractVector,u::AbstractVector) = begin
        z[inds.x] = x
        z[inds.u] = u
        ∇fz(Z,zeros(eltype(x),p),z)
        return nothing
    end
    ∇f!(x::AbstractVector,u::AbstractVector) = begin
        z[inds.x] = x
        z[inds.u] = u
        Z = PartedArray(zeros(eltype(x),p,n+m),inds)
        ∇fz(Z,zeros(eltype(x),p),z)
        return Z
    end
    ∇f!(Z::AbstractMatrix,x::AbstractVector) = ForwardDiff.jacobian!(Z,f!,v0,x)
    return ∇f!, f_aug
end

function generate_jacobian(::Type{Nominal},::Type{Discrete},fd!::Function,n::Int,m::Int)
    inds = (x=1:n,u=n .+ (1:m), dt=n+m+1, xx=(1:n,1:n),xu=(1:n,n .+ (1:m)), xdt=(1:n,n+m.+(1:1)))
    S0 = zeros(n,n+m+1)
    s = zeros(n+m+1)
    ẋ0 = zeros(n)

    fd_aug!(xdot,s) = fd!(xdot,s[inds.x],s[inds.u],s[inds.dt])
    Fd!(S,xdot,s) = ForwardDiff.jacobian!(S,fd_aug!,xdot,s)
    ∇fd!(S::AbstractMatrix,ẋ::AbstractVector,x::AbstractVector,u::AbstractVector,dt::Float64) = begin
        s[inds.x] = x
        s[inds.u] = u
        s[inds.dt] = dt
        Fd!(S,ẋ,s)
        return nothing
    end
    ∇fd!(S::AbstractMatrix,x::AbstractVector,u::AbstractVector,dt::Float64) = begin
        s[inds.x] = x
        s[inds.u] = u
        s[inds.dt] = dt
        Fd!(S,zero(x),s)
        return nothing
    end
    ∇fd!(x::AbstractVector,u::AbstractVector,dt::Float64) = begin
        s[inds.x] = x
        s[inds.u] = u
        s[inds.dt] = dt
        S0 = zeros(eltype(x),n,n+m+1)
        Fd!(S0,zero(x),s)
        return S0
    end
    return ∇fd!, fd_aug!
end

function generate_jacobian(::Type{Uncertain},::Type{Continuous},f!::Function,n::Int,m::Int,r::Int)
    inds = (x=1:n,u=n .+ (1:m), w=(n+m) .+ (1:r), xx=(1:n,1:n),xu=(1:n,n .+ (1:m)), xw=(1:n,(n+m) .+ (1:r)))
    Z = PartedArray(zeros(n,n+m+r),inds)
    z = zeros(n+m+r)
    v0 = zeros(n)
    f_aug(dZ::AbstractVector,z::AbstractVector) = f!(dZ,view(z,inds.x), view(z,inds.u), view(z,inds.w))
    ∇fz(Z::AbstractMatrix,v::AbstractVector,z::AbstractVector) = ForwardDiff.jacobian!(Z,f_aug,v,z)

    ∇f!(Z::AbstractMatrix,v::AbstractVector,x::AbstractVector,u::AbstractVector,w::AbstractVector) = begin
        z[inds.x] = x
        z[inds.u] = u
        z[inds.w] = w
        ∇fz(Z,v,z)
        return nothing
    end
    ∇f!(Z::AbstractMatrix,x::AbstractVector,u::AbstractVector,w::AbstractVector) = begin
        z[inds.x] = x
        z[inds.u] = u
        z[inds.w] = w
        ∇fz(Z,zeros(eltype(x),n),z)
        return nothing
    end
    ∇f!(x::AbstractVector,u::AbstractVector,w::AbstractVector) = begin
        z[inds.x] = x
        z[inds.u] = u
        z[inds.w] = w
        Z = PartedArray(zeros(eltype(x),n,n+m+r),inds)
        ∇fz(Z,zeros(eltype(x),n),z)
        return Z
    end
    return ∇f!, f_aug
end

function generate_jacobian(::Type{Uncertain},::Type{Discrete},fd!::Function,n::Int,m::Int,r::Int)
    inds = (x=1:n,u=n .+ (1:m),w=(n+m) .+ (1:r), dt=n+m+1, xx=(1:n,1:n),xu=(1:n,n .+ (1:m)), xw=(1:n, (n+m) .+ (1:r)), xdt=(1:n,(n+m+r) .+ (1:1)))
    S0 = zeros(n,n+m+r+1)
    s = zeros(n+m+r+1)
    ẋ0 = zeros(n)

    fd_aug!(xdot,s) = fd!(xdot,s[inds.x],s[inds.u],s[inds.w],s[inds.dt])
    Fd!(S,xdot,s) = ForwardDiff.jacobian!(S,fd_aug!,xdot,s)
    ∇fd!(S::AbstractMatrix,ẋ::AbstractVector,x::AbstractVector,u::AbstractVector,w::AbstractVector,dt::Float64) = begin
        s[inds.x] = x
        s[inds.u] = u
        s[inds.w] = w
        s[inds.dt] = dt
        Fd!(S,ẋ,s)
        return nothing
    end
    ∇fd!(S::AbstractMatrix,x::AbstractVector,u::AbstractVector,w::AbstractVector,dt::Float64) = begin
        s[inds.x] = x
        s[inds.u] = u
        s[inds.w] = w
        s[inds.dt] = dt
        Fd!(S,zero(x),s)
        return nothing
    end
    ∇fd!(x::AbstractVector,u::AbstractVector,w::AbstractVector,dt::Float64) = begin
        s[inds.x] = x
        s[inds.u] = u
        s[inds.w] = w
        s[inds.dt] = dt
        S0 = zeros(eltype(x),n,n+m+r+1)
        Fd!(S0,zero(x),s)
        return S0
    end
    return ∇fd!, fd_aug!
end

"""$(SIGNATURES)
Convert a continuous dynamics model into a discrete one using the given discretization function.
    The discretization function can either be one of the currently supported functions (midpoint, rk3, rk4) or a custom method that has the following form
    ```
    function discretizer(f::Function,dt::Float64)
        function fd!(xdot,x,u,dt)
            # Your code
            return nothing
        end
        return fd!
    end
    ```
"""
function discretize_model(model::Model{M,Continuous},discretizer::Symbol, dt::T=1.0) where {M<:ModelType,T}
    fd!,∇fd! = discretize(model.f,discretizer,dt,model.n,model.m)
    info_d = deepcopy(model.info)
    integration = string(discretizer)
    info_d[:integration] = Symbol(replace(integration,"TrajectoryOptimization." => ""))
    info_d[:fc] = model.f
    info_d[:∇fc] = model.∇f
    AnalyticalModel{M,Discrete}(fd!, ∇fd!, model.n, model.m, model.r, model.params, info_d)
end

function discretize_model(model::Model{Uncertain,Continuous},discretizer::Symbol,dt::T=1.0) where T
    fd!,∇fd! = discretize_uncertain(model.f,discretizer,dt,model.n,model.m,model.r)
    info_d = deepcopy(model.info)
    integration = string(discretizer)
    info_d[:integration] = Symbol(replace(integration,"TrajectoryOptimization." => ""))
    info_d[:fc] = model.f
    info_d[:∇fc] = model.∇f
    AnalyticalModel{Uncertain,Discrete}(fd!, ∇fd!, model.n, model.m, model.r, model.params, info_d)
end

function continuous(model::Model{M,Discrete}) where M<:ModelType
    n,m = model.n, model.m
    if haskey(model.info, :fc)
        f = model.info[:fc]
        if haskey(model.info, :∇fc)
            ∇f = model.info[:∇fc]
            Model(f,∇f,n,m)
        else
            Model(f,n,m)
        end
    else
        error("Continuous time model not stashed in model.info")
    end
end

midpoint(model::Model{M,Continuous},dt::T=1.0) where {M <: ModelType,T} = discretize_model(model, :midpoint, dt)
rk3(model::Model{M,Continuous},dt::T=1.0) where {M <: ModelType,T} = discretize_model(model, :rk3, dt)
rk4(model::Model{M,Continuous},dt::T=1.0) where {M <:ModelType,T} = discretize_model(model, :rk4, dt)

rk3_implicit(model::Model{M,Continuous},dt::T=1.0) where {M,T} = discretize_model(model,:rk3_implicit,dt)
midpoint_implicit(model::Model{M,Continuous},dt::T=1.0) where {M,T} = discretize_model(model,:midpoint_implicit,dt)

function discretize(f::Function,discretization::Symbol,dt::T,n::Int,m::Int) where T

    if discretization in [:rk3,:rk4,:midpoint] # ie, explicit
            fd! = eval(discretization)(f,dt)
    elseif discretization in [:rk3_implicit,:midpoint_implicit] # ie, implicit
        fd! = eval(discretization)(f,n,m,dt)
    elseif String(discretization)[1:6] == "DiffEq"
            fd! = DiffEqIntegrator(f,dt,Symbol(split(String(discretization),"_")[2]),n,m)
    else
        error("Integration not defined")
    end
    ∇fd!, = generate_jacobian(Nominal,Discrete,fd!,n,m)
    return fd!,∇fd!
end

function discretize_uncertain(f::Function,discretization::Symbol,dt::T,n::Int,m::Int,r::Int) where T

    if discretization in [:rk3,:rk4,:midpoint] # ie, explicit
            fd! = eval(Symbol(String(discretization) * "_uncertain"))(f,dt)
    elseif discretization in [:rk3_implicit,:midpoint_implicit] # ie, implicit
        fd! = eval(Symbol(String(discretization) * "_uncertain"))(f,n,m,r,dt)
    elseif String(discretization)[1:6] == "DiffEq"
        fd! = DiffEqIntegratorUncertain(f,dt,Symbol(split(String(discretization),"_")[2]),n,m,r)
    else
        error("Integration not defined")
    end

    ∇fd!, = generate_jacobian(Uncertain,Discrete,fd!,n,m,r)
    return fd!,∇fd!
end

"""
$(SIGNATURES)
Determine if the dynamics in model are in place. i.e. the function call is of
the form `f!(xdot,x,u)`, where `xdot` is modified in place. Returns a boolean.
"""
function is_inplace_dynamics(model::Model)::Bool
    is_inplace_dynamics(model.f,model.n,model.m,model.r)
end

function is_inplace_dynamics(f::Function,n::Int,m::Int,r::Int)::Bool
    x = rand(n)
    u = rand(m)
    xdot = rand(n)
    if r == 0
        try
            f(xdot,x,u)
        catch x
            if x isa MethodError
                return false
            end
        end
    else
        w = rand(r)
        try
            f(xdot,x,u,w)
        catch x
            if x isa MethodError
                return false
            end
        end
    end
    return true
end
"""
$(SIGNATURES)
Makes the dynamics function `f(x,u)` appear to operate as an inplace operation of the
form `f!(xdot,x,u)`.
"""
function wrap_inplace(f::Function)
    f!(xdot,x,u) = copyto!(xdot, f(x,u))
    f!(xdot,x) = copyto!(xdot, f(x))
end

"""$(SIGNATURES)
Checks jacobians of functions of the form `f(v,x,u)` to make sure they have the correct forms.
Jacobians should have the following three forms
```∇f(x,u)
∇f(Z,x,u)
∇f(Z,v,x,u)```
"""
function _test_jacobian(::Type{Continuous},∇f::Function)
    form = [true,true,true]
    form[1] = hasmethod(∇f,(AbstractVector,AbstractVector))
    form[2] = hasmethod(∇f,(AbstractMatrix,AbstractVector,AbstractVector))
    form[3] = hasmethod(∇f,(AbstractMatrix,AbstractVector,AbstractVector,AbstractVector))
    return form
end

"""$(SIGNATURES)
Checks jacobians of functions of the form `f(v,x,u,dt)` to make sure they have the correct forms.
Jacobians should have the following three forms
```∇f(x,u,dt)
∇f(S,x,u,dt)
∇f(S,v,x,u,dt)```
"""
function _test_jacobian(::Type{Discrete},∇f::Function)
    form = [true,true,true]
    form[1] = hasmethod(∇f,(AbstractVector,AbstractVector,Float64))
    form[2] = hasmethod(∇f,(AbstractMatrix,AbstractVector,AbstractVector,Float64))
    form[3] = hasmethod(∇f,(AbstractMatrix,AbstractVector,AbstractVector,AbstractVector,Float64))
    return form
end

function _check_dynamics(f::Function,n::Int,m::Int)
    no_params = hasmethod(f,(AbstractVector,AbstractVector,AbstractVector))
    with_params = hasmethod(f,(AbstractVector,AbstractVector,AbstractVector,Any))
    return [no_params,with_params]
end

"Add slack controls to dynamics to make artificially fully actuated"
function add_slack_controls(model::Model{Nominal,D}) where D<:Discrete
    n = model.n; m = model.m
    nm = n+m

    idx = merge(create_partition((m,n),(:u,:inf)),(x=1:n,))
    idx2 = [(1:nm)...,2n+m+1]

    function f!(x₊::AbstractVector{T},x::AbstractVector{T},u::AbstractVector{T},dt::T) where T
        model.f(x₊,x,u[idx.u],dt)
        x₊ .+= u[idx.inf]
    end

    function ∇f!(Z::AbstractMatrix{T},x::AbstractVector{T},u::AbstractVector{T},dt::T) where T
        model.∇f(view(Z,idx.x,idx2),x[idx.x],u[idx.u],dt)
        view(Z,idx.x,(idx.x) .+ nm) .= Diagonal(1.0I,n)
    end

    AnalyticalModel{Nominal,Discrete}(f!,∇f!,n,nm,model.r,model.params,model.info)
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
module ComradeDynesty

using Comrade

using AbstractMCMC
using TupleVectors
using Reexport

@reexport using Dynesty

Comrade.samplertype(::Type{<:NestedSampler}) = Comrade.IsCube()
Comrade.samplertype(::Type{<:DynamicNestedSampler}) = Comrade.IsCube()

function AbstractMCMC.sample(post::Comrade.TransformedPosterior,
                             sampler::Union{NestedSampler, DynamicNestedSampler},
                             args...;
                             kwargs...)
    ℓ = logdensityof(post)
    kw = delete!(Dict(kwargs), :init_params)
    res = sample(ℓ, identity, sampler, args...; kw...)
    samples, weights = res["samples"].T, exp.(res["logwt"].T .- res["logz"][end])
    chain = transform.(Ref(post), eachcol(samples)) |> TupleVector
    stats = (logl = res["logl"].T,
             logz = res["logz"][end],
             logzerr = res["logz"][end],
             weights = weights,
            )
    return TupleVector(chain), stats
end

end

ext: jl


text:
function solve_modelARCDiag(nlp_stop, PData :: PDataFact, α:: T) where T
    # Solve the ARC subproblem once diagonalized into Δ
    # printstyled("On est dans solve_modelARCDiag ⇊ \n", color = :red)
    # Setup the problem
    M = ones(T, size(PData.Δ))

    ϵ1 = sqrt(eps(T)) ./ T(100.0)
    ϵ = T.(ϵ1 .* (1.0 .+ PData.λ))
    λ = max(ϵ, PData.λ .+ ϵ)

    d̃ = -(PData.Δ .+ λ .* M) .\ PData.g̃
    normd̃ = sqrt(d̃⋅d̃)
    seuil = normd̃ ./ α

    # Solve the subproblem (Δ + λ I) d̃ = -g̃ such that λ = ||d̃||/α
    d̃,λ = solve_diag(λ, PData.Δ, PData.g̃, seuil, α, ϵ)

    PData.λ = λ
    # Transform back d̃ into d

    d = AInv(PData,d̃)

    return d, NaN * rand(length(d)), λ
end


text:
import ..TopOpt: whichdevice

whichdevice(s::SIMPResult) = whichdevice(s.topology)
whichdevice(s::SIMP) = whichdevice(s.optimizer)
whichdevice(c::ContinuationSIMP) = whichdevice(c.simp)
whichdevice(o::Optimizer) = whichdevice(o.model)
whichdevice(o::BoxOptimizer) = whichdevice(o.obj)


text:
module Tirocinio
#=

import Pluto
Pluto.run()
=#

#= TEST ORIGINALI

DAF
 ---
    contaminante: Tetracloroetilene (PCE)
    tessitura: sand
    metodo: Domenico/Schwartz
    estensione: 10
    densità: non messo
    spessore: non messo
    concentrazione: 100
    profondità: 1000
    darcy: 0.000025
    profondità mixed: 1580
    dispersione x: non messo
    dispersione y: non messo
    indice 1° decadimento: non messo
 ---
    python daf.py -s 5 -t sand -c0 10000 -x 1 -y 1 -sw 10 -l 1

    python daf.py -s 5 -t sand -c0 10000 -x 1 -y 1 -sw 10 -l 1 -g 0.5 -T 5

    nuovo daf

    python daf.py -s 5 -t sand -c0 5 -x 2 -y 0 -T 8 -v 0.09 -opz 0

    python daf.py -s 5 -t sand -c0 5 -x 30 -y 0 -T 8 -v 0.09 -opz 1

    python daf.py -s 5 -t sand -c0 5 -x 30 -y 0 -T 8 -v 0.09 -opz 1 -lev 0

    python daf.py -s 56 -t sand -c0 5 -x 1000 -y 1580 -T 8 -v 0.09 -opz 0 -sw 100

Leaching
    python leaching.py -s 57 -t sand -lf 1 -v 0.000025 -dgw 1587

Plume
 ---
    wind_dir: E (Est == 0?)
    stabilità: a
    outdoor: country (c?)
    concentrazione: 10000
    temperatura: 18
    h stack: 80
    diametro: 1
    temperatura gas: 150
    velocità gas:0.1
    velocità vento: 1
 ---
    python3 plume.py      -q 100 -u   5 -x    4 -y 1 -z 1 -c A -o c -h_s 10 -x_w 45

    python3 plume.py      -q  10 -u 4.5 -x 8000 -y 0 -z 0 -c A -o c -h_s 25 -x_w 90

    python3 plume_loop.py -q  10 -u 4.5 -x 8000 -y 0 -z 0 -c A -o c -h_s 25 -x_w 90

    python3 plume_loop.py -q 720 -u 4.5 -x  500 -y 1 -z 0 -c A -o c -h_s 62 -x_w 90 -d_s 3 -v_s 11.8 -t_s 159

    concentration 100 - 10 - 10 - 720
    wind_speed 5 - 4.5 - 4.5 - 4.5
    coordinates (4, 1, 1) - (8000, 0, 0) - (8000, 0, 0) - (500, 1, 0)
    stability a - a - a - a
    outdoor c - c - c - c
    stack_height 10 - 25 - 25 - 62
    wind_direction 45 - 90 - 90 - 90
    stack_diameter ? - ? - ? - 3
    gas_speed ? - ? - ? - 11.8
    temperature ? - ? - ? - 159

Lake
 ---
    direzione corrente: E (0?)
    concentrazione: 2000
    velocità corrente: 0.03
    coeff. fickian x: 4
    coeff. fickian y: 3
    coeff. lambda: non messo
    tempo: 10
 ---
    python lake.py -Cs 5 -t 2 -x 5 -y 5 -Vx 1 -Vy 1 -Dx 2 -Dy 40

    python lake.py -Cs 2000 -t 20 -x 120 -y 120 -Vx 1 -Vy 1 -Dx 1000 -Dy 1000

River
 ---
    concentrazione 2000
    sezione bagnata: 1
    raggio idraulico: 0.1
    coeff. dispersione: 15
    coeff. scabrezza: 14
    coeff. decadimento: non messo
    tempo: 2
 ---
    python river.py -Cs 20000 -t 60 -x 120 -v 2 -Dl 100 -w 10

Sediment
    python3 sediment.py -Q 4 -t 9999 -hh 13 -dt 1 -dx 1 -dy 10 -x0 0 -y0 0 -x 10 -y 10 -V 1 -w 0.0359
    dredged_mass 4
    time 9999
    mean_depth 13
    time_interval 1
    dispersion x 1
    dispersion y 10
    x0 0
    y0 0
    x 10
    y 10
    mean_flow_speed 1
    mean_sedimentation_velocity 0.0359

Noise
 ---
    wind dir: NW (135?)
    frequenza: 500
    velocità vento: 6
    umidità: 70
    temperatura: 25
    intensità: 110
 ---

=#

using Revise
using ArchGDAL

include(".\\Analysis\\DiluitionAttenuationFactor.jl")
include(".\\Analysis\\Lakes.jl")
include(".\\Analysis\\Plumes.jl")
include(".\\Analysis\\Sediments.jl")

const agd = ArchGDAL
const daf = DiluitionAttenuationFactor
const lks = Lakes
const plm = Plumes
const sed = Sediments

src = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\source_shapefile\\source_32.shp"
dtm = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\DTM_32.tiff"
#   source = agd.read(src)

#=
    contaminante: Tetracloroetilene (PCE)
    concentrazione: 100
    profondità: 1000
    direzione flusso: non si vede
    pioggia: non si vede
    tessitura: sand
    tempo: non si vede
    estensione: 10
    densità: non messo
    spessore: non messo
    darcy: 0.000025
    profondità mixed: 1580
    indice 1° decadimento: non messo
    metodo: Domenico/Schwartz
=#
daf.run_leaching(
    dem_file = src,
    source_file = src,
    contaminants = ["Tetracloroetilene (PCE)"],
    concentrations = [100.0],
    aquifier_depth = 1000.0,
    aquifier_flow_direction = angle,
    mean_rainfall = rain,
    texture = "sand",
    resolution = 25,
    time = time,
    orthogonal_extension = 10.0,
    darcy_velocity = 0.000025,
    mixed_zone_depth = 1580.0,
    option = :domenico,
    output_path = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis results\\test_daf.tiff"
)

#=
    direzione corrente: E (0?)
    concentrazione: 2000
    velocità corrente: 0.03
    tempo: 10
    coeff. fickian x: 4
    coeff. fickian y: 3
    coeff. lambda: non messo
=#
lks.run_lake(
    source_file = src,
    wind_direction = 0,
    pollutant_mass = 2000.0,
    flow_mean_speed = 0.03,
    resolution = 25,
    hours = 10,
    fickian_x = 4.0,
    fickian_y = 3.0,
    output_path = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis results\\test_lake.tiff"
)

#=
    stabilità: a
    outdoor: country (c?)
    concentrazione: 10000
    wind_dir: E (niente vento? Est? Est == 0?)
    velocità vento: 1
    h stack: 80
    velocità gas:0.1
    diametro: 1
    temperatura gas: 150
    temperatura: 18
=#
plm.run_plume(
    dem_file = dtm,
    source_file = src,
    stability = "a",
    outdoor = "c",
    concentration = 10000.0,
    resolution = 25,
    wind_direction = 0,
    wind_speed = 1.0,
    stack_height = 80.0,
    gas_speed = 0.1,
    stack_diameter = 1.0,
    gas_temperature = 180.0,
    temperature = 18.0,
    output_path = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis results\\test_plume.tiff"
)

using ArchGDAL
include(".\\Analysis\\Utils\\Functions.jl")
const agd = ArchGDAL
src = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\source_shapefile\\source_32.shp"
dtm = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\DTM_32.tiff"
d = agd.read(dtm)
point = agd.getgeom( collect(agd.getlayer(agd.read(src), 0))[1], 0 )
src_i = Functions.toIndexes(d, agd.getx(point, 0), agd.gety(point, 0))
src_c = Functions.toCoords(d, src_i...)

using Revise
using ArchGDAL
include(".\\Analysis\\Plumes.jl")
const plm = Plumes
src = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\source_shapefile\\source_32.shp"
dtm = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\DTM_32.tiff"
out = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis results\\test_plume.tiff"
pr = plm.run_plume(dem_file=dtm,source_file=src,stability="a",outdoor="c",concentration=10000.0,resolution=25,wind_direction=0,wind_speed=1.0,stack_height=80.0,gas_speed=0.1,stack_diameter=1.0,gas_temperature=180.0,temperature=18.0,output_path=out)

#=  PROBLEMI
- Il risultato è spostato verso il basso rispetto alla sorgente [V]
- source è sull'angolo della cella source del risultato [V]
- Il risultato sembra essere girato (Il vento va verso Est ma nell'immagine il fumo va reso Ovest)
- L'angolo sembra non fare nulla (sembra dipendere dall'intesità del vento) [V]
- Quando il vento è particolarmente forte  si ottengono risultati strani (le polveri vngono portate a grande distanza quindi
    si ha il non verificarsi della condizione nelle immediate vicinanze della sorgente)
=#
#= POSSIBILI SOLUZIONI
Posizione
- Sottrarre alla Y del geotransform invece che sommare [X]
- Usare minrow mincolumn [V]
- Usare maxcolumn maxrow
- Usare mincolumn maxrow
- Invertire minRow e maxColumn sempre nel geotransform
Angolo
- cambiare l'angolo come se lo 0 fosse al posto di 270 [X]
- cambiare l'angolo come se lo 0 fosse al posto di 180 [X]
=#

using Revise
using ArchGDAL
include(".\\Analysis\\Lakes.jl")
const lks = Lakes
src = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\source_shapefile\\source_32.shp"
dtm = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\DTM_32.tiff"
out = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis results\\test_lake.tiff"
lr = lks.run_lake(dem_file=dtm,source_file=src,wind_direction=0,pollutant_mass=2000.0,flow_mean_speed=0.03,resolution=25,hours=10,fickian_x=4.0,fickian_y=3.0,output_path=out)

d = agd.read(dtm)
s = agd.read(src)

geom = agd.getgeom( collect(agd.getlayer(s, 0))[1] )

sr = agd.toWKT(agd.getspatialref(geom))
dr = agd.getproj(d)

end # module

include(".\\Analysis\\Utils\\Functions.jl")
using ArchGDAL
const agd = ArchGDAL
out = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis results\\test_plume.tiff"
dtm = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\DTM_32.tiff"
src = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\source_shapefile\\source_32.shp"
d = agd.read(dtm)
shp = agd.read(src)
p = agd.read(out)
coords = Functions.toCoords(p, 52, 12)
point = agd.getgeom(collect(agd.getlayer(shp, 0))[1])

#=

using GRASS_jll

ccal(
    # (<FunctionName>, <Lybrary>)
    (:viewshed, :spgrass7),
)

=#

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
# ---
# title: 788. Rotated Digits
# id: problem788
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/rotated-digits/description/>
# hidden: true
# ---
#
# X is a good number if after rotating each digit individually by 180 degrees,
# we get a valid number that is different from X.  Each digit must be rotated -
# we cannot choose to leave it alone.
#
# A number is valid if each digit remains a digit after rotation. 0, 1, and 8
# rotate to themselves; 2 and 5 rotate to each other (on this case they are
# rotated in a different direction, in other words 2 or 5 gets mirrored); 6 and
# 9 rotate to each other, and the rest of the numbers do not rotate to any other
# number and become invalid.
#
# Now given a positive number `N`, how many numbers X from `1` to `N` are good?
#
#
#
#     **Example:**
#     Input: 10
#     Output: 4
#     Explanation:
#     There are four good numbers in the range [1, 10] : 2, 5, 6, 9.
#     Note that 1 and 10 are not good numbers, since they remain unchanged after rotating.
#
#
# **Note:**
#
#   * N  will be in range `[1, 10000]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# ----------------------------------------------------------------------------------------------- #
# Default lean-project things

using LeanProject

se = joinpath(dirname(@__DIR__), "input", "setting.json")

PAR, PAI, PAC, PAO = LeanProject.get_project_path(se)

SE = LeanProject.read_setting(se)

# ----------------------------------------------------------------------------------------------- #
# Project-specific things

using OnePiece


text:
"""
    struct SampledDensity <: AbstractDensity

Stores a density and samples drawn from it.

A report on the density, samples, properties of variates/parameters, etc. can
be generated via `Base.show`.

Constructors:

* ```SampledDensity(density::AbstractPosteriorDensity, samples::DensitySampleVector)```

Fields:

* ```density::AbstractDensity```

* ```samples::DensitySamplesVector```

* ```_stats::BAT.MCMCBasicStats```

* ```_generator::AbstractSampleGenerator```

!!! note

    Fields `_stats` and `_generator` do not form part of the stable public
    API and are subject to change without deprecation.

This type is likely to evolve into a subtype of `AbstractDensity` in future
versions.
"""
struct SampledDensity{D<:AbstractDensity,S<:DensitySampleVector, G<:AbstractSampleGenerator} <: AbstractDensity
    density::D
    samples::S
    _stats::MCMCBasicStats
    _generator::G
end
export SampledDensity

function SampledDensity(density::AbstractDensity, samples::DensitySampleVector)
    stats = MCMCBasicStats(samples)
    return SampledDensity(density, samples, stats, UnknownSampleGenerator())
end

function SampledDensity(density::SampledDensity, samples::DensitySampleVector)
    return SampledDensity(density.density, samples)
end

function SampledDensity(target::AnyDensityLike, samples::DensitySampleVector)
    return SampledDensity(convert(AbstractDensity, target), samples)
end

eval_logval(density::SampledDensity, v::Any, T::Type{<:Real}) = eval_logval(density.density, v, T)

eval_logval_unchecked(density::SampledDensity, v::Any) = eval_logval_unchecked(density.density, v)

ValueShapes.varshape(density::SampledDensity) = varshape(density.density)

var_bounds(density::SampledDensity) = var_bounds(density.density)

# ToDo: Distributions.sampler(density::SampledDensity)
# ToDo: bat_sampler(density::SampledDensity)

get_initsrc_from_target(target::SampledDensity) = target.samples

_get_deep_prior_for_trafo(density::SampledDensity) = _get_deep_prior_for_trafo(density.density)

function bat_transform_impl(target::Union{PriorToUniform,PriorToGaussian}, density::SampledDensity, algorithm::PriorSubstitution)
    new_parent_density, trafo = bat_transform_impl(target, density.density, algorithm)
    new_samples = trafo.(density.samples)
    (result = SampledDensity(new_parent_density, new_samples), trafo = trafo)
end

# ToDo: truncate_density(density::SampledDensity, bounds::AbstractArray{<:Interval})

_approx_cov(target::SampledDensity) = cov(target.samples)

function getdensity(sd::SampledDensity)
    return sd.density
end

function nfreeparams(sd::SampledDensity)
    return totalndof(varshape(sd))
end

function freeparams(sd::SampledDensity)
    return active_keys(sd)
end

function nfixedparams(sd::SampledDensity)
    return length(get_fixed_names(varshape(sd)))
end

function fixedparams(sd::SampledDensity)
    param_shape = varshape(sd)
    fixed_param_keys = Symbol.(get_fixed_names(param_shape))
    fixed_values = [getproperty(param_shape, f).shape.value for f in fixed_param_keys]
    return (; zip(fixed_param_keys, fixed_values)...,)
end

function numberofsamples(sd::SampledDensity)
    return length(sd.samples)
end

function eff_sample_size(sd::SampledDensity)
    effsize = bat_eff_sample_size(unshaped.(sd.samples)).result
    return (; zip(active_keys(sd), effsize)...,)

end

function active_keys(sd::SampledDensity)
    return Symbol.(all_active_names(elshape(sd.samples.v)))
end

function Statistics.mean(sd::SampledDensity)
    means = sd._stats.param_stats.mean
    return (; zip(active_keys(sd), means)...,)
end

function Statistics.std(sd::SampledDensity)
    covm = collect(sd._stats.param_stats.cov)
    stds = sqrt.(LinearAlgebra.diag(covm))

    return (; zip(active_keys(sd), stds)...,)
end

function Distributions.mode(sd::SampledDensity)
    modes = sd._stats.mode
    return (; zip(active_keys(sd), modes)...,)
end

function marginalmode(sd::SampledDensity)
    modes = bat_marginalmode(unshaped.(sd.samples)).result
    return (; zip(active_keys(sd), modes)...,)
end

function Statistics.cov(sd::SampledDensity)
    covm = collect(sd._stats.param_stats.cov)
    names = string.(active_keys(sd))

    return NamedArrays.NamedArray(covm, (names, names), ("cov",""))
end

function Statistics.cor(sd::SampledDensity)
    covm = collect(sd._stats.param_stats.cov)
    corm = cov2cor(covm, sqrt.(diag(covm)))
    names = string.(active_keys(sd))

    return NamedArrays.NamedArray(corm, (names, names), ("cor",""))
end

function parameter_table(sd::SampledDensity)
    tab = TypedTables.Table(
        parameter = active_keys(sd),
        mean = collect(mean(sd)),
        std = collect(std(sd)),
        global_mode = collect(mode(sd)),
        marginal_mode = collect(marginalmode(sd))
        )
    return tab
end

function fixed_parameter_table(sd::SampledDensity)
    fixed = fixedparams(sd)
    freekeys = collect(keys(fixed))

    tab = TypedTables.Table(
        parameter = freekeys,
        value = collect(fixed)
        )
    return tab
end

function Base.show(io::IO, mime::MIME"text/plain", sd::SampledDensity)
    println(io, "BAT.jl - SampledDensity")
    _line(io, length=30)

    println(io, "\nSampling:")
    _line(io, length=25)
    _print_generator(io, sd._generator)
    _print_sampling(io, sd)

    fpt = fixed_parameter_table(sd)
    if !isempty(fpt)
        println(io, "\n\nFixed parameters:")
        _line(io, length=25)
        println(io, "number of fixed parameters: ", nfixedparams(sd))
        println(io, "fixed parameters: ",fixedparams(sd))
    end

    println(io, "\n\nParameter estimates:")
    _line(io, length=25)
    println(io, "number of free parameters: ", nfreeparams(sd), "\n")
    println(io, parameter_table(sd))

    println(io, "\n\nCovariance matrix:")
    _line(io, length=25)
    println(io, cov(sd))
end

function _line(io::IO; length=25, indent=0)
    println(io, repeat(' ', indent), repeat('─', length))
end

function _print_sampling(io::IO, sd::SampledDensity)
    println(io, "total number of samples:", repeat(' ', 6), numberofsamples(sd))
    println(io, "effective number of samples: ", eff_sample_size(sd))
end

function _print_generator(io::IO, generator::GenericSampleGenerator)
    print(io, "algorithm: ")
    show(io, "text/plain", getalgorithm(generator))
    print(io, "\n")
end

function _print_generator(io::IO, generator::UnknownSampleGenerator)
end

function _print_generator(io::IO, generator::MCMCSampleGenerator)
    chains = generator._chains
    nchains = length(chains)
    n_tuned_chains = count(c -> c.info.tuned, chains)
    n_converged_chains = count(c -> c.info.converged, chains)
    print(io, "algorithm: ")
    show(io, "text/plain", getalgorithm(generator))
    print(io, "\n")

    println(io, "number of chains:", repeat(' ', 13), nchains)
    println(io, "number of chains tuned:", repeat(' ', 7), n_tuned_chains)
    println(io, "number of chains converged:", repeat(' ', 3), n_converged_chains)
    println(io, "number of samples per chain:", repeat(' ', 2), nsamples(chains[1]), "\n")
end


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

using FEMBase, Test, LinearAlgebra, SparseArrays, Statistics

include(joinpath("..", "docs", "make.jl"))

@testset "FEMBase.jl" begin
    @testset "test_assembly" begin include("test_assembly.jl") end
    @testset "test_elements" begin include("test_elements.jl") end
    @testset "test_fields" begin include("test_fields.jl") end
    @testset "test_integrate" begin include("test_integrate.jl") end
    @testset "test_problems" begin include("test_problems.jl") end
    @testset "test_sparse" begin include("test_sparse.jl") end
    @testset "test_solvers" begin include("test_solvers.jl") end
    @testset "test_analysis" begin include("test_analysis.jl") end
    @testset "test_test" begin include("test_test.jl") end
    @testset "test_types" begin include("test_types.jl") end
end

include(joinpath("..", "docs", "deploy.jl"))


text:
export exampledense, examplesparse
import HDF5

"""
    exampledense(x, name, dims; density = 1, compress = 6, chunkdim = nothing)

Create an example dense dataset with `name` inside a HDF5 file/group.
`x` may be either a string containing the name of a new HDF5 file, or a handle into an existing HDf5 file or group.
`dims` should be a tuple specifying the dimensions of the dataset.

The HDF5 dataset is created with random double-precision values at the specified `density`.
Lowering the `density` will introduce zero values randomly into the dataset.
Note that this setting does not affect the structural density of the dataset.

`compress` specifies the DEFLATE compression level - if set to zero, no compression is used.
`chunkdim` should be a tuple of length equal to `dims`, specifying the chunk dimensions to be used for compression.
This is heuristically determined if not supplied, and is ignored if `compress = 0`.

Nothing is returned by this function.

# Examples
```jldoctest
julia> using HDF5Arrays

julia> tmp = tempname();

julia> exampledense(tmp, "stuff", (10, 5, 2))

julia> y = DenseHDF5Array(tmp, "stuff");

julia> size(y)
(10, 5, 2)
```
"""
function exampledense(x, name::String, dims::Tuple{Vararg{Integer}}; density = 1, compress = 6, chunkdim = nothing)
    val = Array{Float64,length(dims)}(undef, dims)
    for i in 1:length(val)
        if rand() <= density
            val[i] = rand()
        else
            val[i] = 0
        end
    end

    if compress > 0
        if chunkdim == nothing
            x[name,compress=compress] = val
        else
            x[name,compress=compress,chunk=chunkdim] = val
        end
    else
        x[name] = val
    end

    return
end

function exampledense(x::String, name::String, dims::Tuple{Vararg{Integer}}; density = 1, compress = 6, chunkdim = nothing)
    handle = HDF5.h5open(x, "w");
    exampledense(handle, name, dims; density = density, compress = compress, chunkdim = chunkdim)
    return
end

"""
    examplesparse(x, name, dims, density; compress = 6, chunkdim = 10000)

Create an example sparse dataset with `name` inside a HDF5 file/group.
`x` may be either a string containing the name of a new HDF5 file, or a handle into an existing HDf5 file or group.
`dims` should be a tuple specifying the dimensions of the dataset.

The sparse matrix is created with random double-precision values at the specified `density`.
Contents are written to file inside group `name` in the 10X Genomics format, which uses the compressed sparse column layout.
Briefly: values are stored in `data`, row indices in `indices`, column pointers in `indptr` and the dimensions in `shape`.

`compress` specifies the DEFLATE compression level - if set to zero, no compression is used.
`chunkdim` should be a tuple of length equal to `dims`, specifying the chunk dimensions to be used for compression.
This is heuristically determined if not supplied, and is ignored if `compress = 0`.

Nothing is returned by this function.

# Examples
```jldoctest
julia> using HDF5Arrays

julia> tmp = tempname();

julia> examplesparse(tmp, "stuff", (10, 20), 0.1)

julia> y = SparseHDF5Matrix(tmp, "stuff");

julia> size(y)
(10, 20)
```
"""
function examplesparse(x, name::String, dims::Tuple{Integer, Integer}, density::AbstractFloat; compress = 6, chunkdim = 10000)
    ref = SparseArrays.sprand(dims[1], dims[2], density)
    handle = HDF5.create_group(x, name)

    vals = SparseArrays.getnzval(ref)
    for i in 1:length(vals)
        vals[i] -= 1
    end

    ind = SparseArrays.getrowval(ref)
    for i in 1:length(ind)
        ind[i] -= 1
    end

    ptr = SparseArrays.getcolptr(ref)
    for i in 1:length(ptr)
        ptr[i] -= 1
    end

    if compress > 0
        if chunkdim == nothing
            handle["data",compress=compress] = vals
            handle["indices",compress=compress] = ind
            handle["indptr",compress=compress] = ptr
        else
            handle["data",compress=compress,chunkdim=chunkdim] = vals
            handle["indices",compress=compress,chunkdim=chunkdim] = ind
            handle["indptr",compress=compress,chunkdim=chunkdim] = ptr
        end
    else
        handle["data"] = vals
        handle["indices"] = ind
        handle["indptr"] = ptr
    end
    handle["shape"] = [dims...]

    return
end

function examplesparse(x::String, name::String, dims::Tuple{Integer, Integer}, density::AbstractFloat; compress = 6, chunkdim = nothing)
    handle = HDF5.h5open(x, "w");
    examplesparse(handle, name, dims, density; compress = compress, chunkdim = chunkdim)
    return
end


text:
include("../envs/RLEnvironments.jl")
using .RLEnvironments
using Plots

act(W, s) = argmax(W * s)

function run!(W, env; epochs=200, maxt=100, noise=1f-2, γ=99f-2)
    rewards = zeros(Float32, epochs)
    steps = zeros(Float32, epochs)

    bestR = typemin(Float32)
    bestW = copy(W)
    for i in 1:epochs
        reset!(env)
        done = false
        Ri = 0f0
        j = 0f0
        while !done && j < maxt
            _, r, done = step!(env, act(W, env.current))
            Ri += γ^j * r
            j += 1f0
        end

        if Ri >= bestR
            bestR = Ri
            bestW = copy(W)
            noise = max(1f-3, noise/2f0)
            W += noise .* rand(Float32, size(W))
        else
            noise = min(2f0, noise*2f0)
            W = bestW + noise .* rand(Float32, size(W))
        end
        rewards[i] = Ri
        steps[i] = j
    end
    W .= bestW
    return rewards, steps
end

env = CartPole{Float32}()
W = 1f-4 * randn(Float32, length(env.actionspace), length(env.observationspace))

r, s = run!(W, env)
plt = plot([r s], labels=["Reward" "Steps"], color=[:red :blue], layout=(2, 1))
display(plt)
#savefig(plt, "hillclimbing+adaptivenoisescaling")


text:
using SymPy.Physics
using SymPy: Sym, sqrt, symbols, doit, PI
using SymPy.SpecialFuncs: Ynm
using PyCall: PyError
using Test

@testset "Physics" begin
    @test clebsch_gordan(Sym(3)/2, Sym(1)/2, 2, Sym(3)/2, Sym(1)/2, 2) == 1
    @test clebsch_gordan(Sym(3)/2, Sym(1)/2, 1, Sym(3)/2, -Sym(1)/2, 1) == sqrt(Sym(3))/2
    @test clebsch_gordan(Sym(3)/2, Sym(1)/2, 1, -Sym(1)/2, Sym(1)/2, 0) == -sqrt(Sym(2))/2

    θ, ϕ = symbols("theta, phi")
    @test doit(dot_rot_grad_Ynm(3, 2, 2, 0, θ, ϕ)) == 3*sqrt(Sym(55))*Ynm(5, 2, θ, ϕ)/(11*sqrt(PI))

    @test gaunt(1,0,1,1,0,-1) == -1/(2*sqrt(PI))
    #@test N(gaunt(1000,1000,1200,9,3,-12)) ≈ 0.00689500421922113448 # takes forever to compute

    @test_throws PyError gaunt(1.2,0,1.2,0,0,0)
    @test_throws PyError gaunt(1,0,1,1.1,0,-1.1)

    @test racah(3,3,3,3,3,3) == -1//14

    @test wigner_3j(2, 6, 4, 0, 0, 0) == sqrt(Sym(715))/143
    @test wigner_3j(2, 6, 4, 0, 0, 1) == 0

    @test wigner_6j(3,3,3,3,3,3) == -1//14
    @test wigner_6j(5,5,5,5,5,5) == 1//52

    # Optics
    @test RayTransferMatrix(1,2,3,4) == Sym[1 2; 3 4]
    @test FlatMirror() == Sym[1 0; 0 1]

    # Spin
    @test doit(WignerD(Sym(3), Sym(2), Sym(1), PI, PI/2, -PI)) == -sqrt(Sym(10))/8
    @test doit(WignerD(Sym(1), Sym(1), Sym(0), 0, θ, 0)) == -sin(θ)/sqrt(Sym(2))
end


text:
using Distributions: Uniform
using LinearAlgebra
using StaticArrays
using Random: rand
using Test

using CLIMA
using CLIMA.Atmos
using CLIMA.ConfigTypes
using CLIMA.Diagnostics
using CLIMA.GenericCallbacks
using CLIMA.ODESolvers
using CLIMA.ColumnwiseLUSolver: ManyColumnLU
using CLIMA.Mesh.Filters
using CLIMA.Mesh.Grids
using CLIMA.MoistThermodynamics
using CLIMA.VariableTemplates

using CLIMAParameters
using CLIMAParameters.Planet: R_d, day, grav, cp_d, cv_d, planet_radius
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

struct HeldSuarezDataConfig{FT}
    T_ref::FT
end

function init_heldsuarez!(bl, state, aux, coords, t)
    FT = eltype(state)

    # Parameters need to set initial state
    T_ref = bl.data_config.T_ref
    temp_profile = IsothermalProfile(T_ref)

    # Calculate the initial state variables
    T, p = temp_profile(bl.orientation, bl.param_set, aux)
    thermo_state = PhaseDry_given_pT(p, T, bl.param_set)
    ρ = air_density(thermo_state)
    e_int = internal_energy(thermo_state)
    e_pot = gravitational_potential(bl.orientation, aux)

    # Set initial state with random perturbation
    rnd = FT(1.0 + rand(Uniform(-1e-3, 1e-3)))
    state.ρ = ρ
    state.ρu = SVector{3, FT}(0, 0, 0)
    state.ρe = rnd * state.ρ * (e_int + e_pot)

    nothing
end

function config_heldsuarez(FT, poly_order, resolution)
    exp_name = "HeldSuarez"

    # Parameters
    T_ref::FT = 255
    Rh_ref::FT = 0
    domain_height::FT = 30e3
    turb_visc::FT = 0 # no visc. here

    # Set up a reference state for linearization
    temp_profile_ref = IsothermalProfile(T_ref)
    ref_state = HydrostaticState(temp_profile_ref, Rh_ref)

    # Rayleigh sponge to dampen flow at the top of the domain
    z_sponge = FT(15e3) # height at which sponge begins
    α_relax = FT(1 / 60 / 15) # sponge relaxation rate in (1/seconds)
    u_relax = SVector(FT(0), FT(0), FT(0)) # relaxation velocity
    exp_sponge = 2 # sponge exponent for squared-sinusoid profile
    sponge = RayleighSponge{FT}(
        domain_height,
        z_sponge,
        α_relax,
        u_relax,
        exp_sponge,
    )

    # Set up the atmosphere model
    model = AtmosModel{FT}(
        AtmosGCMConfigType,
        param_set;
        ref_state = ref_state,
        turbulence = ConstantViscosityWithDivergence(turb_visc),
        moisture = DryModel(),
        source = (Gravity(), Coriolis(), held_suarez_forcing!, sponge),
        init_state = init_heldsuarez!,
        data_config = HeldSuarezDataConfig(T_ref),
    )

    config = CLIMA.AtmosGCMConfiguration(
        exp_name,
        poly_order,
        resolution,
        domain_height,
        param_set,
        init_heldsuarez!;
        model = model,
    )

    return config
end

function held_suarez_forcing!(
    bl,
    source,
    state,
    diffusive,
    aux,
    t::Real,
    direction,
)
    FT = eltype(state)

    # Parameters
    T_ref = bl.data_config.T_ref

    # Extract the state
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe

    coord = aux.coord
    e_int = internal_energy(bl.moisture, bl.orientation, state, aux)
    T = air_temperature(e_int, bl.param_set)
    _R_d = FT(R_d(bl.param_set))
    _day = FT(day(bl.param_set))
    _grav = FT(grav(bl.param_set))
    _cp_d = FT(cp_d(bl.param_set))
    _cv_d = FT(cv_d(bl.param_set))

    # Held-Suarez parameters
    k_a = FT(1 / (40 * _day))
    k_f = FT(1 / _day)
    k_s = FT(1 / (4 * _day))
    ΔT_y = FT(60)
    Δθ_z = FT(10)
    T_equator = FT(315)
    T_min = FT(200)
    σ_b = FT(7 / 10)
    λ = longitude(bl, aux)
    φ = latitude(bl, aux)
    z = altitude(bl, aux)
    scale_height = _R_d * T_ref / _grav
    σ = exp(-z / scale_height)

    # TODO: use
    #  p = air_pressure(T, ρ)
    #  σ = p/p0
    exner_p = σ^(_R_d / _cp_d)
    Δσ = (σ - σ_b) / (1 - σ_b)
    height_factor = max(0, Δσ)
    T_equil = (T_equator - ΔT_y * sin(φ)^2 - Δθ_z * log(σ) * cos(φ)^2) * exner_p
    T_equil = max(T_min, T_equil)
    k_T = k_a + (k_s - k_a) * height_factor * cos(φ)^4
    k_v = k_f * height_factor

    # Apply Held-Suarez forcing
    source.ρu -= k_v * projection_tangential(bl, aux, ρu)
    source.ρe -= k_T * ρ * _cv_d * (T - T_equil)
    return nothing
end

function config_diagnostics(FT, driver_config)
    interval = 1000 # in time steps

    _planet_radius = FT(planet_radius(param_set))

    info = driver_config.config_info
    boundaries = [
        FT(-90.0) FT(-180.0) _planet_radius
        FT(90.0) FT(180.0) FT(_planet_radius + info.domain_height)
    ]
    resolution = (FT(10), FT(10), FT(1000))
    interpol =
        CLIMA.InterpolationConfiguration(driver_config, boundaries, resolution)

    dgngrp = setup_dump_state_and_aux_diagnostics(
        interval,
        driver_config.name,
        interpol = interpol,
        project = true,
    )
    return CLIMA.DiagnosticsConfiguration([dgngrp])
end

function main()
    CLIMA.init()

    # Driver configuration parameters
    FT = Float32                        # floating type precision
    poly_order = 5                      # discontinuous Galerkin polynomial order
    n_horz = 5                          # horizontal element number
    n_vert = 5                          # vertical element number
    days = 100                          # experiment day number
    timestart = FT(0)                   # start time (seconds)
    timeend = FT(days * 24 * 60 * 60)   # end time (seconds)

    # Set up driver configuration
    driver_config = config_heldsuarez(FT, poly_order, (n_horz, n_vert))

    # Set up ODE solver configuration
    ode_solver_type = CLIMA.IMEXSolverType(
        linear_solver = ManyColumnLU,
        solver_method = ARK2GiraldoKellyConstantinescu,
    )

    # Set up experiment
    solver_config = CLIMA.SolverConfiguration(
        timestart,
        timeend,
        driver_config,
        ode_solver_type = ode_solver_type,
        Courant_number = 0.14,
        init_on_cpu = true,
        CFL_direction = HorizontalDirection(),
    )

    # Set up diagnostics
    dgn_config = config_diagnostics(FT, driver_config)

    # Set up user-defined callbacks
    # TODO: This callback needs to live somewhere else
    filterorder = 10
    filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            1:size(solver_config.Q, 2),
            solver_config.dg.grid,
            filter,
        )
        nothing
    end

    # Run the model
    result = CLIMA.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbfilter,),
        check_euclidean_distance = true,
    )
end

main()

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
const J0_ROOTS(::Type{Float64}) = (
    (2.4048255576957730, -1.176691651530894e-16),
    (3.8317059702075125, -1.5269184090088067e-16),
    (5.5200781102863110, 8.088597146146722e-17),
    (7.0155866698156190, -9.414165653410389e-17),
    (8.6537279129110130, -2.92812607320779e-16),
    (10.173468135062722, 4.482162274768888e-16),
    (11.791534439014281, 2.812956912778735e-16),
    (13.323691936314223, 2.600408064718813e-16),
    (14.930917708487787, -7.070514505983074e-16),
    (16.470630050877634, -1.619019544798128e-15),
    (18.071063967910924, -9.658048089426209e-16),
    (19.615858510468243, -1.004445634526616e-15),
    (21.211636629879260, 4.947077428784068e-16),
    (22.760084380592772, -4.925749373614922e-16),
    (24.352471530749302, 9.169067133951066e-16),
    (25.903672087618382, 4.894530726419825e-16)
)
const J0_POLYS(::Type{Float64}) = (
    (0.0000000000000000000, -0.5191474972894669, 0.10793870175491979, 0.05660177443794914, -0.008657669593292222, -0.0021942003590739974, 0.0002643770365942964, 4.37291931443113e-5, -4.338825419759815e-6, -5.304927679598406e-7, 4.469819175606667e-8, 4.3284827345621115e-9, -3.135111000732148e-10, -2.628876489517534e-11),
    (-0.402759395702553, 2.476919088072758e-16, 0.20137969785127532, -0.017518715285670765, -0.013352611033152278, 0.0010359438492839443, 0.00037218755624680334, -2.4952042421113972e-5, -5.776086353844158e-6, 3.374317129436161e-7, 5.727482259215452e-8, -2.9561880489355444e-9, -3.905845672635605e-10, 1.971332566705736e-11),
    (0.000000000000000000, 0.34026480655836816, -0.030820651425593214, -0.05298855286760721, 0.004631042145890388, 0.002257440229081131, -0.00017518572879518415, -4.6521091062878115e-5, 3.199785909661533e-6, 5.716500268232186e-7, -3.5112898510841466e-8, -4.684643389757727e-9, 2.562685034682206e-10, 2.7958958795750104e-11),
    (0.30011575252613254, -1.6640272822046001e-16, -0.15005787626306408, 0.007129737603121546, 0.011742619737383848, -0.0006260583453094324, -0.00035093119008693753, 1.7929701912295164e-5, 5.6239324892485754e-6, -2.668437501970219e-7, -5.6648488273749086e-8, 2.48117399780498e-9, 3.8876537586241154e-10, -1.6657136713437192e-11),
    (0.000000000000000000, -0.27145229992838193, 0.015684124960953488, 0.044033774963413, -0.0025093022271948434, -0.0020603351551475315, 0.00011243486771159352, 4.482303558813413e-5, -2.288390108003442e-6, -5.679383588459768e-7, 2.693939234375692e-8, 4.737285529934781e-9, -2.0612709555352797e-10, -2.8163166483726606e-11),
    (-0.2497048770578432, 1.1807897766765572e-16, 0.12485243852891914, -0.0040907858517059345, -0.010102792347641438, 0.00038536375952213334, 0.00031859711440332953, -1.2373899600646271e-5, -5.3013932979548665e-6, 2.001098153528186e-7, 5.4711629662471434e-8, -1.9724572531751518e-9, -3.8121398193699247e-10, 1.3667679743782715e-11),
    (0.000000000000000000, 0.23245983136472478, -0.009857064513825458, -0.03818600911162367, 0.0016073972920762946, 0.0018420433388794816, -7.581358465623415e-5, -4.159284549011371e-5, 1.650645590334915e-6, 5.425453494592871e-7, -2.0556207467977526e-8, -4.620018928884712e-9, 1.642028058414746e-10, 2.7701605444102412e-11),
    (0.21835940724787295, -8.89726402965429e-17, -0.10917970362393398, 0.0027314677279632535, 0.008944552393700088, -0.00026391472261453583, -0.00028847875053074687, 8.858193371737123e-6, 4.9233776180403375e-6, -1.5077786827161215e-7, -5.190218733666561e-8, 1.5539413886301204e-9, 3.674809363354973e-10, -1.1113645791216594e-11),
    (0.00000000000000000, -0.20654643307799603, 0.006916736034268416, 0.034115572697347704, -0.001137276252948717, -0.0016680057255530109, 5.4841792064182565e-5, 3.837965853474541e-5, -1.2335804050962046e-6, -5.106259295634553e-7, 1.592333632709497e-8, 4.423517565793139e-9, -1.3138837384184105e-10, -2.6809397212536384e-11),
    (-0.19646537146865717, 6.979167865106427e-17, 0.09823268573432613, -0.001988037402152532, -0.008095530671166083, 0.00019440675128712672, 0.0002640383898036336, -6.666777683303928e-6, -4.5715696772304925e-6, 1.1666296153560847e-7, 4.8913639696764225e-8, -1.2379867207945651e-9, -3.508930968415813e-10, 9.07632091591013e-12),
    (0.00000000000000000, 0.18772880304043943, -0.005194182350684612, -0.031096513233785917, 0.0008577442641273341, 0.0015312251534677639, -4.184307585284775e-5, -3.5603170534217916e-5, 9.58002109234601e-7, 4.795250964600283e-7, -1.263434500625308e-8, -4.20550167685701e-9, 1.065791122326672e-10, 2.5727417675684577e-11),
    (0.18006337534431555, -5.638484737644332e-17, -0.09003168767215539, 0.0015299132863046024, 0.0074441453680234426, -0.00015060569680378768, -0.0002443398416353605, 5.227001519013193e-6, 4.267152607972633e-6, -9.295966007495808e-8, -4.610438011417262e-8, 1.0049092632275165e-9, 3.339442682325105e-10, -7.4991079301099e-12),
    (0.00000000000000000, -0.17326589422922986, 0.004084217951979124, 0.028749284970146657, -0.0006761643016121907, -0.0014215899173758441, 3.320978125391802e-5, 3.3264379323815026e-5, -7.684444100376941e-7, -4.515479818833484e-7, 1.0271599456634145e-8, 3.993903280527723e-9, -8.793975528824604e-11, -2.4610374225652004e-11),
    (-0.16718460047381803, 4.662597138876655e-17, 0.08359230023690671, -0.0012242529339116864, -0.006925682915280748, 0.00012100729852042988, 0.00022821854128498174, -4.230799709959437e-6, -4.007618360337058e-6, 7.601217108010105e-8, 4.358483157250522e-8, -8.317621452517008e-10, -3.178805429572483e-10, 6.284538399593043e-12),
    (0.000000000000000000, 0.16170155068925002, -0.0033200234006036146, -0.02685937038656159, 0.0005505380905909195, 0.0013316994659124243, -2.715683205388875e-5, -3.1289544794362e-5, 6.326900360777323e-7, 4.2697141781883964e-7, -8.532447325590315e-9, -3.799009445641295e-9, 7.379552811590934e-11, 2.353654960834717e-11),
    (0.15672498625285222, 1.1464880342445208e-19, -0.07836249312642612, 0.0010083833270351796, 0.006501011610557426, -9.993664895655577e-5, -0.00021478298462967253, 3.511086890959515e-6, 3.7857022791122945e-6, -6.350944888510364e-8, -4.135588012575766e-8, 7.135988717747828e-10, 3.2075281131621564e-10, 2.208506585533542e-12)
)
J0_POLY_PIO2(::Type{Float64}) = (
1.000000000000000, -0.25, 0.01562499999999994, -0.00043402777777725544, 6.781684026082576e-6,
-6.781683757550061e-8, 4.709479394601058e-10, -2.4016837144506874e-12, 9.104258208703104e-15
)

J1_ROOTS(::Type{Float64}) = (
        (1.8411837813406593, 4.7898393919093694e-18),
        (3.8317059702075125, -1.5269184090088067e-16),
        (5.3314427735250325, 1.5109105349471405e-16),
        (7.015586669815619, -9.414165653410389e-17),
        (8.536316366346286, -1.5433871213307537e-16),
        (10.173468135062722, 4.482162274768888e-16),
        (11.706004902592063, 7.121366942298246e-16),
        (13.323691936314223, 2.600408064718813e-16),
        (14.863588633909034, -6.265788988781879e-16),
        (16.470630050877634, -1.619019544798128e-15),
        (18.015527862681804, -1.1196999448424267e-16),
        (19.615858510468243, -1.004445634526616e-15),
        (21.16436985918879, 1.7024131380423588e-15),
        (22.760084380592772, -4.925749373614922e-16),
        (24.311326857210776, -2.614798558537172e-16),
        (25.903672087618382, 4.894530726419825e-16)
        )

@inline J1_POLYS(::Type{Float64}) = (
         (0.5818652242815964, 8.973834293036876e-17, -0.20511071214777496, 0.006058948324597746, 0.013801769808047723, -0.0003723170971916689, -0.0003949590750416506, 9.202953658715798e-6, 6.267292697927347e-6, -1.2681571675674473e-7, -6.319326736018798e-8, 1.057852721352679e-9, 4.839653511789343e-10, -2.1534683937074676e-11),
         (0.0, -0.402759395702553, 0.05255614585697725, 0.053410444132722765, -0.005179719245639115, -0.00223312533910116, 0.00017466429071755533, 4.6208700909637e-5, -3.0368632747366735e-6, -5.727790851165856e-7, 3.248220586398651e-8, 4.732693354840469e-9, -2.346792988667587e-10, -2.607612190048509e-11),
         (-0.34612620185379156, 2.6002760076991505e-17, 0.16697453550109428, -0.009678268542879972, -0.012099225779175138, 0.0006654009006841076, 0.0003541389011182841, -1.742720356958137e-5, -5.655294960447268e-6, 2.484316221112215e-7, 5.7101369594782676e-8, -2.2587104440771904e-9, -3.9854619701676605e-10, 1.854508641982197e-11),
         (0.0, 0.30011575252613254, -0.02138921280934146, -0.04697047894974007, 0.0031302917260424344, 0.0021055871432169443, -0.00012550790947431806, -4.499147505171648e-5, 2.401580334889493e-6, 5.665260872870135e-7, -2.7272058756924113e-8, -4.718921752439826e-9, 2.045523137027535e-10, 2.6923755198436e-11),
         (0.27329994163319987, -1.3511933704624036e-16, -0.13477468037992396, 0.005116340346495344, 0.010631861751990702, -0.000448743683849691, -0.0003268000185531749, 1.3382556752353113e-5, 5.363177082527157e-6, -2.0647461181851255e-7, -5.499892476485602e-8, 1.9780335410646812e-9, 3.8456897851545063e-10, -1.5543946040865833e-11),
         (0.0, -0.2497048770578432, 0.012272357555101388, 0.04041116939079181, -0.0019268187972545589, -0.0019115826893618788, 8.661729445925766e-5, 4.241116264477775e-5, -1.800978985430059e-6, -5.47159694982648e-7, 2.1682608288127194e-8, 4.629875997238198e-9, -1.6946465290961774e-10, -2.7001263147077087e-11),
         (-0.23330441717143408, 1.3347222236161232e-16, 0.11580092244607755, -0.0032489977328297747, -0.009372527206041913, 0.0003036138212749204, 0.00029804555521007126, -9.813819271763006e-6, -5.0242293237237345e-6, 1.6136482548346424e-7, 5.251755806094281e-8, -1.6213723779041528e-9, -3.7130918416720704e-10, 1.2840797195120732e-11),
         (0.0, 0.21835940724787295, -0.008194403183877394, -0.03577820957502993, 0.001319573612804675, 0.0017308725061602815, -6.200735153951766e-5, -3.938703730297954e-5, 1.3569938030086903e-6, 5.190652330962483e-7, -1.7082573736581586e-8, -4.464580500726025e-9, 1.3835552663370753e-10, 2.6508049073606443e-11),
         (0.20701265272531905, -1.1572365385200866e-16, -0.10303781563402706, 0.002289729546293356, 0.008432435058846974, -0.00021965941971569214, -0.000272721543802475, 7.36909520166486e-6, 4.682127415369e-6, -1.2611406291992532e-7, -4.9733737497256667e-8, 1.313333637862963e-9, 3.5584811981981077e-10, -1.0581432980865016e-11),
         (0.0, -0.19646537146865717, 0.005964112206447895, 0.032382122684890324, -0.0009720337562254532, -0.0015842303417456071, 4.6667442166602986e-5, 3.6572573450276185e-5, -1.0499611641790313e-6, -4.891787571238003e-7, 1.3609422867092438e-8, 4.2640592578446675e-9, -1.1321195999108497e-10, -2.5685037344174917e-11),
         (-0.18801748852581776, 9.806595982257605e-17, 0.09371909377243136, -0.0017233243363111925, -0.007714266760106843, 0.00016755862384546863, 0.000251827174474777, -5.7335177376628445e-6, -4.372425030172481e-6, 1.0045839852031556e-7, 4.698184959516747e-8, -1.0717704296017814e-9, -3.3952344362582674e-10, 8.773851724781324e-12),
         (0.0, 0.18006337534431555, -0.004589739858905968, -0.029776581472313136, 0.0007530284838491313, 0.0014660390526475475, -3.658900932751018e-5, -3.4137236390211275e-5, 8.366325091910564e-7, 4.6108475352995054e-7, -1.1047200388967438e-8, -4.058793425595929e-9, 9.363513647344476e-11, 2.4720051360744732e-11),
         (0.1734590492857464, -8.319859584484218e-17, -0.08653590193876039, 0.0013568189856187354, 0.007147216520714834, -0.00013295776598760658, -0.00023462959289705953, 4.6031437800947915e-6, 4.103196271673978e-6, -8.18373301294018e-8, -4.44403211204062e-8, 8.872617322304839e-10, 3.2359730989751983e-10, -7.354373969368619e-12),
         (0.0, -0.16718460047381806, 0.0036727588017285762, 0.02770273166133465, -0.0006050364924617241, -0.0013693112504440895, 2.9615596890234806e-5, 3.206096184277923e-5, -6.841058767433935e-7, -4.358877331737043e-7, 9.143763763164119e-9, 3.864037532317281e-9, -7.851667555762691e-11, -2.3728757802646805e-11),
         (-0.1618382095526585, 7.115722920799028e-17, 0.08078219522574913, -0.0011038527651152986, -0.006686444726790514, 0.00010870535083026479, 0.00022030197347228358, -3.791805814546456e-6, -3.871147931227952e-6, 6.80607764766372e-8, 4.2160129307780366e-8, -7.460427055677985e-10, -3.0872139357605097e-10, 6.24158293622982e-12),
         (0.0, 0.15672498625285222, -0.0030251499811055965, -0.026004046442225856, 0.0004996832448037306, 0.001288697907658241, -2.457760943647847e-5, -3.02856199504406e-5, 5.715983781675481e-7, 4.136275513873783e-7, -7.703846567878601e-9, -3.6853569596797437e-9, 6.676439988396374e-11, 2.277101500922319e-11))#for i in 2:2:18

J1_POLY_PIO2(::Type{Float64}) = (0.5, -0.0624999999999989, 0.002604166666657291, -5.42534721917933e-5, 6.781683542660179e-7, -5.651361336587487e-9, 3.36191211106159e-11, -1.4511302591871352e-13)

#for i in 1:2:18
#    display(((i+1)*pi/2-pts[i+1][1], (i+2)*pi/2-pts[i+1][1]))
#    poly = Float64.(Tuple(ratfn_minimax(x->r2(x,pts[i]), (i*pi/2-pts[i][1], (i+1)*pi/2-pts[i][1]), 13, 0)[1]))
#    push!(polys, poly)
#    poly = Float64.(Tuple(ratfn_minimax(x->r(x,pts[i+1]), ((i+1)*pi/2-pts[i+1][1], (i+2)*pi/2-pts[i+1][1]), 12, 0)[1]))
#    push!(polys, (0.0, poly...))
#end


text:
__precompile__()
"""
Simple String Table and load/save functions

Author:     Scott Paul Jones
Copyright:  2017-2018 Gandalf Software, Inc (and any future contributors)
License:    MIT (see https://github.com/JuliaString/StrTables.jl/blob/master/LICENSE.md)

## Public API:

# Exported:
* StrTable takes a Vector of some `AbstractString` type
  and returns a compact string table, that acts as a `AbstractVector{String}`

# Not exported:
* matchfirst takes a sorted string table and a string,
  returns a vector of elements that start with the string
* matchfirstrng takes a sorted string table and a string,
  returns the range of indexes of elements start with the string
* save takes a filename and a collection (of simple types) and
  saves them to the file in a simple format
* load takes a filename and returns a Vector of the values stored by save in the file
"""
module StrTables
export StrTable, PackedTable, AbstractPackedTable, AbstractEntityTable

const V6_COMPAT = VERSION < v"0.7-"

# Utility functions for building tables
export create_vector, sortsplit!, _contains, _replace, _codeunits, cvt_char, parse_hex

parse_hex(T, s) = @static V6_COMPAT ? parse(T, s, 16) : parse(T, s, base=16)
cvt_char(s)     = @static V6_COMPAT ? convert(Vector{Char}, s) : Vector{Char}(s)
_codeunits(s)   = Vector{UInt8}(@static V6_COMPAT ? s : codeunits(s))
_contains(s, r) = @static V6_COMPAT ? ismatch(r, s) : occursin(r, s)
_replace(s, p)  = @static V6_COMPAT ? replace(s, p.first, p.second) : replace(s, p)
create_vector(T, len)  = @static V6_COMPAT ? Vector{T}(len) : Vector{T}(undef, len)
read_vector(s, T, len) = @static V6_COMPAT ? read(s, T, len) : read!(s, create_vector(T, len))

@static if V6_COMPAT
    const copyto! = copy!
    export copyto!
else
    using Dates
    export now
end

function sortsplit!(index::Vector{UInt16}, vec::Vector{Tuple{T, UInt16}}, base) where {T}
    sort!(vec)
    len = length(vec)
    valvec = create_vector(T, len)
    indvec = create_vector(UInt16, len)
    for (i, val) in enumerate(vec)
        valvec[i], ind = val
        indvec[i] = ind
        index[ind] = UInt16(base + i)
    end
    base += len
    valvec, indvec, base
end

abstract type AbstractPackedTable{T} <: AbstractVector{T} end

"""
Compact table
Designed to save memory compared to a `Vector{Vector{S}}`
Allows for fast lookup of ranges when input was sorted
Can be saved/loaded to/from a file quickly
"""
struct PackedTable{T,S,O} <: AbstractPackedTable{T}
    offsetvec::Vector{O}
    namtab::Vector{S}
end

PackedTable(::Type{T}, offvec::Vector{O}, namtab::Vector{S}) where {T,S,O} =
    PackedTable{T,S,O}(offvec, namtab)

_getsize(el::String) = sizeof(el)
_getsize(el::Vector{<:Any}) = length(el)

function pack_table(::Type{T}, ::Type{S}, strvec) where {T,S}
    namvec = Vector{S}()
    offvec = create_vector(UInt32, length(strvec)+1)
    offvec[1] = 0%UInt32
    offset = 0%UInt32
    for (i, str) in enumerate(strvec)
        offset += _getsize(str)
        offvec[i+1] = offset
        append!(namvec, T == String ? _codeunits(str) : str)
    end
    (offset > 0x0ffff
     ? PackedTable{T,S,UInt32}(offvec, namvec)
     : PackedTable{T,S,UInt16}(copyto!(create_vector(UInt16, length(strvec)+1), offvec), namvec))
end

"""Make a single table of a vector of elements of type T"""
PackedTable(strvec::Vector{T}) where {T} = pack_table(T, eltype(T), strvec)
PackedTable(strvec::Vector{String}) = pack_table(String, UInt8, strvec)

const StrTable = PackedTable{T, UInt8} where {T}

"""
Compact string table
Designed to save memory compared to a `Vector{String}`
Allows for fast lookup of ranges when input was sorted
Can be saved/loaded to/from a file quickly
"""
StrTable

StrTable(strvec::Vector{<:AbstractString}) = pack_table(String, UInt8, strvec)

Base.getindex(str::AbstractPackedTable{T}, ind::Integer) where {T} =
    T(str.namtab[str.offsetvec[ind]+1:str.offsetvec[ind+1]])
Base.size(str::AbstractPackedTable) = (length(str.offsetvec)-1,)
Base.IndexStyle(::Type{<:AbstractPackedTable}) = Base.IndexLinear()
@static if VERSION >= v"0.7.0-DEV.5127"
    (Base.iterate(str::AbstractPackedTable, state::Integer=1) =
     state < length(str.offsetvec) ? (getindex(str, state), state+1) : nothing)
else
    Base.start(str::AbstractPackedTable) = 1
    Base.next(str::AbstractPackedTable, state) = (getindex(str, state), state+1)
    Base.done(str::AbstractPackedTable, state) = state >= length(str.offsetvec)
end

# Get all indices that start with a string

_memcmp(v1, v2, l) = ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt), v1, v2, l)
_veccmp(l1, l2, v1, v2) = _memcmp(v1, v2, min(l1, l2))

@inline function _lexcmp(l1, l2, v1, v2)
    c = _veccmp(l1, l2, v1, v2)
    c < 0 ? -1 : c > 0 ? +1 : cmp(l1, l2)
end

_ltvec(v1, v2) = _lexcmp(sizeof(v1), sizeof(v2), v1, v2) < 0

"""Return the range of indices of values that whose beginning matches the string"""
matchfirstrng(tab::AbstractPackedTable, str::AbstractString) = matchfirstrng(tab, String(str))
matchfirstrng(tab::AbstractPackedTable, str::String) = matchfirstrng(tab, _codeunits(str))
function matchfirstrng(tab::AbstractPackedTable, str::Vector{T}) where {T}
    pos = searchsortedfirst(tab, str, lt=_ltvec)
    len = length(tab)
    pos > len && return pos:pos-1
    beg = pos
    l1 = sizeof(str)
    prevoff = tab.offsetvec[pos]
    while pos <= len
        curoff = tab.offsetvec[pos+1]
        l1 > curoff - prevoff && break
        _memcmp(str, pointer(tab.namtab, prevoff+1), l1) != 0 && break
        prevoff = curoff
        pos += 1
    end
    beg:pos-1
end

"""Return a vector of values that whose beginning matches the string"""
matchfirst(tab::AbstractPackedTable, str) = tab[matchfirstrng(tab, str)]

## Support for AbstractEntityTables

export lookupname, matchchar, matches, longestmatches, completions

abstract type AbstractEntityTable <: AbstractVector{String} end

"""
Abstract type for Entity tables:
Supports lookupname, matchchar, matches, longestmatches, completions
"""
AbstractEntityTable

Base.getindex(str::AbstractEntityTable, ind::Integer) = String(str.nam[ind])
Base.size(str::AbstractEntityTable) = (length(str.ind),)
@static if VERSION < v"0.7.0-DEV.5127"
Base.done(str::AbstractEntityTable, state) = state == length(str.ind)
end

const _empty_str = ""
const _empty_str_vec = Vector{String}()

"""Given an entity name, return the string it represents, or an empty string if not found"""
function lookupname end

"""Given a character, return all exact matches to the character as a vector"""
function matchchar end

"""Given a string, return all exact matches to the string as a vector"""
function matches end

"""Given a string, return all of the longest matches to the beginning of the string as a vector"""
function longestmatches end

"""Given a string, return all of the entity names that start with that string, if any"""
function completions end

matchchar(tab::AbstractEntityTable, ch) = matchchar(tab, UInt32(ch))
matches(tab::AbstractEntityTable, str::AbstractString) = matches(tab, cvt_char(str))
longestmatches(tab::AbstractEntityTable, str::AbstractString) = longestmatches(tab, cvt_char(str))

completions(tab::AbstractEntityTable, str::AbstractString) =
    matchfirst(tab.nam, convert(String, str))

## default methods for most common packed string tables

_get_val2c(tab::AbstractEntityTable, val) = string(Char(val>>>16), Char(val&0xffff))

function _get_str(tab::AbstractEntityTable, ind)
    ind <= tab.base32 && return string(Char(tab.val16[ind]))
    ind <= tab.base2c && return string(Char(tab.val32[ind - tab.base32] + 0x10000))
    _get_val2c(tab, tab.val2c[ind - tab.base2c])
end

function _get_strings(tab::AbstractEntityTable, val::T,
                      vec::AbstractVector{T}, ind::Vector{UInt16}) where {T}
    rng = searchsorted(vec, val)
    isempty(rng) ? _empty_str_vec : tab.nam[ind[rng]]
end

function lookupname(tab::AbstractEntityTable, str::AbstractString)
    rng = searchsorted(tab.nam, str)
    isempty(rng) ? _empty_str : _get_str(tab, tab.ind[rng.start])
end

function matchchar(tab::AbstractEntityTable, ch::UInt32)
    (ch <= 0x0ffff
     ? _get_strings(tab, ch%UInt16, tab.val16, tab.ind16)
     : (ch <= 0x1ffff
        ? _get_strings(tab, ch%UInt16, tab.val32, tab.ind32)
        : _empty_str_vec))
end

function matches(tab::AbstractEntityTable, vec::Vector{T}) where {T}
    if length(vec) == 1
        matchchar(tab, vec[1])
    elseif length(vec) == 2 && ((v1 = vec[1]%UInt32) <= 0x0ffff && (v2 = vec[2]%UInt32) <= 0x0ffff)
        _get_strings(tab, v1<<16 | v2, tab.val2c, tab.ind2c)
    else
        _empty_str_vec
    end
end

function longestmatches(tab::AbstractEntityTable, vec::Vector{T}) where {T}
    isempty(vec) && return _empty_str_vec
    if length(vec) >= 2 && (v1 = vec[1]%UInt32) <= 0x0ffff && (v2 = vec[2]%UInt32) <= 0x0ffff
        res = _get_strings(tab, v1<<16 | v2, tab.val2c, tab.ind2c)
        isempty(res) || return res
        # Fall through and check only the first character
    end
    matchchar(tab, vec[1])
end

# Support for saving and loading

# This use a simple format to store single Unsigned values (UInt8-UInt64),
# Vectors of Unsigned values, and StrTable values

const VECTOR_CODE = 0x0
const STRTAB_CODE = 0x1
const PACKED_CODE = 0x2
const STRING1_CODE = 0x3
const STRING2_CODE = 0x4
const STRING4_CODE = 0x5
const BASE_CODE = 0x5

const type_tab = (UInt8, UInt16, UInt32, UInt64, UInt128,
                  Int8, Int16, Int32, Int64, Int128,
                  Float16, Float32, Float64)
const SupTypes = Union{type_tab...}

for (i, typ) in enumerate(type_tab)
    @eval _get_code(::Type{$typ})  = $((i+BASE_CODE)%UInt8)
end

const MAX_CODE = (length(type_tab)+BASE_CODE)%UInt8

const VER = 0x00000001

write_value(io::IO, val::T) where {T<:SupTypes} = write(io, _get_code(T), val)

write_value(io::IO, val::AbstractString) = write_value(io, String(val))
function write_value(io::IO, str::String)
    siz = sizeof(str)
    if siz < 256
        write(io, STRING1_CODE, sizeof(str)%UInt8, str)
    elseif siz < 65536
        write(io, STRING2_CODE, sizeof(str)%UInt16, str)
    elseif siz <= 0xffffffff
        write(io, STRING4_CODE, sizeof(str)%UInt32, str)
    else
        error("String too large: $siz")
    end
end

write_value(io::IO, val::Vector{T}) where {T<:SupTypes} =
    (write(io, _get_code(T) | 0x80, length(val)%UInt32, val))

function write_value(io::IO, tab::PackedTable{T,S,O}) where {T,S,O}
    write(io, (T == String && S == UInt8) ? STRTAB_CODE : PACKED_CODE)
    write_value(io, tab.offsetvec)
    write_value(io, tab.namtab)
end

"""Save a collection of values (StrTable, String, Float*, UInt*, Int*) into an IO object"""
function save(io::IO, values)
    write(io, VER)
    for val in values
        write_value(io, val)
    end
end

"""Save a collection of values (StrTable, String, Float*, UInt*, Int*) into a file"""
function save(filename::AbstractString, values)
    open(filename, "w") do io
        save(io, values)
    end
end

"""Read a single StrTable, String, UInt, or Vector of UInt value from the io device"""
function read_value(io::IO)
    typ = read(io, UInt8)
    # Check for vector
    if (BASE_CODE|0x80) < typ <= (MAX_CODE|0x80)
        len = read(io, UInt32)
        res = read_vector(io, type_tab[typ-0x80-BASE_CODE], len)
    elseif typ == STRTAB_CODE
        off = read_value(io)
        nam = read_value(io)
        res = PackedTable(String, off, nam)
    elseif typ == PACKED_CODE
        off = read_value(io)
        nam = read_value(io)
        res = PackedTable(typeof(nam), off, nam)
    elseif typ == STRING1_CODE
        len = read(io, UInt8)
        res = String(read(io, len))
    elseif typ == STRING2_CODE
        len = read(io, UInt16)
        res = String(read(io, len))
    elseif typ == STRING4_CODE
        len = read(io, UInt32)
        res = String(read(io, len))
    elseif typ <= MAX_CODE
        res = read(io, type_tab[typ-BASE_CODE])
    else
        error("Unsupported type code $typ")
    end
end

"""Return a Vector{Any} with data tables loaded from IO"""
function load(io::IO)
    res = Vector{Any}()
    ver = read(io, UInt32)
    ver != VER && error("File version $ver doesn't match current version $VER")
    while !eof(io)
        push!(res, read_value(io))
    end
    res
end

"""Return a Vector{Any} with data tables loaded from file"""
function load(filename::AbstractString)
    open(filename) do io
        load(io)
    end
end

end # module StrTables


text:
# testing discretization procedures
# example is based on Walter Gautschi "Orthogonal Polynomials: Computation and Approximation"
# examples 2.36

using PolyChaos, Test, StaticArrays
import LinearAlgebra: norm

nodes = [40, 80, 160, 320]
tol = 1e-14

@testset "Stieltjes procedure" begin
    for n in nodes
        myfile = open("fej$n.txt")
        βref = parse.(Float64,readlines(myfile))
        βcom = stieltjes(n,fejer(n)...)[2]
        nodes, weights = fejer(n)
        @test stieltjes(n,SVector(nodes...),SVector(weights...)) == stieltjes(n,nodes,weights) == stieltjes(n,nodes,SVector(weights...)) == stieltjes(n,SVector(nodes...),weights)
        @test isapprox(norm(βref-βcom,Inf),0.;atol=tol)
        close(myfile)
    end
end

@testset "Lanczos procedure" begin
    for n in nodes
        myfile = open("fej$n.txt")
        βref = parse.(Float64,readlines(myfile))
        nodes, weights = fejer(n)
        βcom = lanczos(n,nodes,weights)[2]
        @test lanczos(n,SVector(nodes...),SVector(weights...)) == lanczos(n,nodes,weights) == lanczos(n,nodes,SVector(weights...)) == lanczos(n,SVector(nodes...),weights)
        @test isapprox(norm(βref-βcom,Inf),0.;atol=tol)
        close(myfile)
    end
end


text:
include("colors.jl")
include("points.jl")
include("quaternion.jl")
include("sarray.jl")
include("statistic.jl")
include("triangle.jl")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

julia_version = v"1.6.0"

name = "scopehal"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/azonenberg/scopehal.git", "db29f07e396667135ccf02b37e27f68b8d4998a1"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/scopehal/
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/0001-replace-color-with-string.patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/0001-avx2-fma-on-llvm.patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/0001-rigol-function-generator.patch
git submodule update --init
cd xptools
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/0001-musl.patch
cd ../..
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DJulia_PREFIX=${prefix} \
    .
make -j${nproc}
make install
install_license scopehal/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")
platforms = libjulia_platforms(julia_version)
# ARM targets are not supported
platforms = filter!(p -> BinaryBuilder.proc_family(p) == "intel", platforms)
# LLVM does not yet support OpenMP on BinaryBuilder
 filter!(p -> !Sys.isapple(p), platforms)
platforms = expand_cxxstring_abis(platforms; skip=p->false)

# The products that we will ensure are always built
products = Product[
    LibraryProduct("libscopehal", :scopehal)
    LibraryProduct("libscopehalwrapper", :scopehalwrapper)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="yaml_cpp_jll", uuid="01fea8cc-7d33-533a-824e-56a766f4ffe8"))
    Dependency(PackageSpec(name="ffts_jll", uuid="84ad2fbc-dacc-5b40-99f5-8db9d02a0a8a"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
    Dependency("libcxxwrap_julia_jll")
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"7.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil