text:
export MvNormalMeanCovariance

import Distributions: logdetcov, distrname, sqmahal, sqmahal!, AbstractMvNormal
import LinearAlgebra: diag, Diagonal, dot
import Base: ndims, precision, length, size, prod

struct MvNormalMeanCovariance{ T <: Real, M <: AbstractVector{T}, P <: AbstractMatrix{T} } <: AbstractMvNormal
    μ :: M
    Σ :: P
end

function MvNormalMeanCovariance(μ::AbstractVector{ <: Real }, Σ::AbstractMatrix{ <: Real })
    T = promote_type(eltype(μ), eltype(Σ))
    return MvNormalMeanCovariance(convert(AbstractArray{T}, μ), convert(AbstractArray{T}, Σ))
end

function MvNormalMeanCovariance(μ::AbstractVector{ <: Integer}, Σ::AbstractMatrix{ <: Integer })
    return MvNormalMeanCovariance(float.(μ), float.(Σ))
end

function MvNormalMeanCovariance(μ::AbstractVector, σ::AbstractVector)
    return MvNormalMeanCovariance(μ, matrix_from_diagonal(promote_type(eltype(μ), eltype(σ)), σ))
end

function MvNormalMeanCovariance(μ::AbstractVector{T}) where T
    return MvNormalMeanCovariance(μ, convert(AbstractArray{T}, ones(length(μ))))
end

Distributions.distrname(::MvNormalMeanCovariance) = "MvNormalMeanCovariance"

weightedmean(dist::MvNormalMeanCovariance) = invcov(dist) * mean(dist)

function weightedmean_invcov(dist::MvNormalMeanCovariance)
    W = invcov(dist)
    xi = W * mean(dist)
    return (xi, W)
end

Distributions.mean(dist::MvNormalMeanCovariance)      = dist.μ
Distributions.var(dist::MvNormalMeanCovariance)       = diag(cov(dist))
Distributions.cov(dist::MvNormalMeanCovariance)       = dist.Σ
Distributions.invcov(dist::MvNormalMeanCovariance)    = cholinv(dist.Σ)
Distributions.std(dist::MvNormalMeanCovariance)       = cholsqrt(cov(dist))
Distributions.logdetcov(dist::MvNormalMeanCovariance) = logdet(cov(dist))

Distributions.sqmahal(dist::MvNormalMeanCovariance, x::AbstractVector) = sqmahal!(similar(x), dist, x)

function Distributions.sqmahal!(r, dist::MvNormalMeanCovariance, x::AbstractVector)
    μ = mean(dist)
    for i in 1:length(r)
        @inbounds r[i] = μ[i] - x[i]
    end
    return xT_A_x(r, invcov(dist))
end

Base.eltype(::MvNormalMeanCovariance{T})     where T = T
Base.precision(dist::MvNormalMeanCovariance)         = invcov(dist)
Base.length(dist::MvNormalMeanCovariance)            = length(mean(dist))
Base.ndims(dist::MvNormalMeanCovariance)             = length(dist)
Base.size(dist::MvNormalMeanCovariance)              = (length(dist), )

Base.convert(::Type{ <: MvNormalMeanCovariance }, μ::AbstractVector, Σ::AbstractMatrix) = MvNormalMeanCovariance(μ, Σ)

function Base.convert(::Type{ <: MvNormalMeanCovariance{T} }, μ::AbstractVector, Σ::AbstractMatrix) where { T <: Real }
    return MvNormalMeanCovariance(convert(AbstractArray{T}, μ), convert(AbstractArray{T}, Σ))
end

vague(::Type{ <: MvNormalMeanCovariance }, dims::Int) = MvNormalMeanCovariance(zeros(dims), fill(huge, dims))

prod_analytical_rule(::Type{ <: MvNormalMeanCovariance }, ::Type{ <: MvNormalMeanCovariance }) = ProdAnalyticalRuleAvailable()

function Base.prod(::ProdPreserveType, left::MvNormalMeanCovariance, right::MvNormalMeanCovariance)
    invcovleft  = invcov(left)
    invcovright = invcov(right)
    Σ = cholinv(invcovleft + invcovright)
    μ = Σ * (invcovleft * mean(left) + invcovright * mean(right))
    return MvNormalMeanCovariance(μ, Σ)
end

function Base.prod(::ProdAnalytical, left::MvNormalMeanCovariance, right::MvNormalMeanCovariance)
    xi_left, W_left = weightedmean_precision(left)
    xi_right, W_right = weightedmean_precision(right)
    return MvNormalWeightedMeanPrecision(xi_left + xi_right, W_left + W_right)
end


text:
using EigenCpp
using Base.Test

const rows = 4
const nb_elems = 1000000

julia_matrices = Array{Matrix{Float64}}((nb_elems,));
julia_vectors = Array{Vector{Float64}}((nb_elems,));
julia_results = Array{Vector{Float64}}((nb_elems,));

function init_julia!(mats::Vector{Matrix{Float64}}, vecs::Vector{Vector{Float64}}, result::Vector{Vector{Float64}})
  for i in 1:nb_elems
    mats[i] = ones(rows,rows)
    vecs[i] = ones(rows)*2.
    result[i] = zeros(rows)
  end
end

function matvec_julia!(mats::Vector{Matrix{Float64}}, vecs::Vector{Vector{Float64}}, result::Vector{Vector{Float64}})
    for i in 1:length(mats)
        #result[i] += mats[i]*vecs[i]
        @inbounds LinAlg.BLAS.gemv!('N', 1., mats[i], vecs[i], 1., result[i])
    end
end

function check_result{T}(result::Vector{T})
  const ref = ones(rows)*24
  all(result) do vec
    vec == ref
  end
end

println("Julia init time:")
@time init_julia!(julia_matrices, julia_vectors, julia_results)

println("Julia matvec timings:")
@time matvec_julia!(julia_matrices, julia_vectors, julia_results)
@time matvec_julia!(julia_matrices, julia_vectors, julia_results)
@time matvec_julia!(julia_matrices, julia_vectors, julia_results)

println("Julia check time:")
@time @test check_result(julia_results)

const EigenMat = EigenCpp.Matrix{Float64,rows,rows}
const EigenVec = EigenCpp.Matrix{Float64,rows,1}

eigen_matrices = Array{EigenMat}((nb_elems,));
eigen_vectors = Array{EigenVec}((nb_elems,));
eigen_results = Array{EigenVec}((nb_elems,));

function init_eigen!(mats::Vector{EigenMat}, vecs::Vector{EigenVec}, result::Vector{EigenVec})
  for i in 1:nb_elems
    mats[i] = EigenMat()
    vecs[i] = EigenVec()
    result[i] = EigenVec()
    EigenCpp.setConstant(mats[i], 1)
    EigenCpp.setConstant(vecs[i], 2)
    EigenCpp.setConstant(result[i], 0)
  end
end

function matvec_eigen!(mats::Vector{EigenMat}, vecs::Vector{EigenVec}, result::Vector{EigenVec})
    for i in 1:length(mats)
        #result[i] += mats[i]*vecs[i]
        EigenCpp.plus_assign_product!(mats[i], vecs[i], result[i])
    end
end

println("Eigen init time:")
@time init_eigen!(eigen_matrices, eigen_vectors, eigen_results)

println("Eigen matvec timings:")
@time matvec_eigen!(eigen_matrices, eigen_vectors, eigen_results)
@time matvec_eigen!(eigen_matrices, eigen_vectors, eigen_results)
@time matvec_eigen!(eigen_matrices, eigen_vectors, eigen_results)

println("Eigen check time:")
@time @test check_result(eigen_results)

ext: jl


text:
"""
    PosixPath()
    PosixPath(str)

Represents any posix path (e.g., `/home/user/docs`)
"""
struct PosixPath <: SystemPath
    segments::Tuple{Vararg{String}}
    root::String
end

PosixPath() = PosixPath(tuple(), "")
PosixPath(segments::Tuple; root="") = PosixPath(segments, root)

function Base.tryparse(::Type{PosixPath}, str::AbstractString)
    str = string(str)
    isempty(str) && return PosixPath(tuple("."))

    tokenized = split(str, POSIX_PATH_SEPARATOR)
    root = isempty(tokenized[1]) ? POSIX_PATH_SEPARATOR : ""
    return PosixPath(tuple(map(String, filter!(!isempty, tokenized))...), root)
end

function Base.expanduser(fp::PosixPath)::PosixPath
    p = fp.segments

    if p[1] == "~"
        return length(p) > 1 ? joinpath(home(), p[2:end]...) : home()
    end

    return fp
end

function Base.Filesystem.contractuser(fp::PosixPath)
    h = home()
    if isdescendant(fp, h)
        if fp == h
            return PosixPath("~")
        else
            n = length(h.segments)
            return PosixPath(("~", fp.segments[n+1:end]...))
        end
    end

    return fp
end


text:
#Careful! needs LinearAlgebra in Reconstruction.jl
"""
 SeisSpitzFX(d;<keyword arguments>)

 First order f-x interpolation using Spitz' method.
 Input size is size(d)= (nt,nh). Output size is (nt,2*nh-1)

# Arguments
-`dt`: sampling interval in secs
-`npf`: prediction filter length
- `pre1` : percentage of pre-whitening for Pef estimation
- `pre2` : percentage of pre-whitening for estimation of missing samples
- `flow` : minimun temporal frequency to reconstruct (Hz)
- `fhigh` : maximun temporal frequency to reconstruct (Hz)

# Example
'''julia
d_interp = SeisSpitzFX(d)
'''
*  Reference: Spitz, S., 1991, Seismic trace interpolation in the F-X domain: Geophysics, 56, 785-794.
"""
function SeisSpitzFX(d;dt=0.004,npf=2,pre1=1,pre2=1,flow=1,fhigh=90)

# 1) Fourier Transform to FX domain

 vect = size(d);
 nt = vect[1]
 nf1 = nextpow(2,vect[1]);
 nf2 = 2*nf1
 nh = vect[2]

 d1,ifl1,ifh1 = pad(d,dt,flow,fhigh)
 D1 = fft(d1,1)

 d2,ifl2,ifh2 = pad(d,dt,flow,fhigh;npad=2)
 D2 = fft(d2,1)

 INTDF = zeros(typeof(D2[1]),nf1, 2*nh-1);

 for ia = ifl1:ifh1;

# 2) Select frequency slices from frequencies f and f/2

    x1 = transpose(D1[ia,:]);
    x2 = transpose(D2[ia,:]);
# 3) Estime the PFs needed for frequency sample ia;

    PF = prediction_filter(x2,npf,pre1);

# 4) Interpolate spatial data at freq. sample ia
    y = interpolate_freq(x1,PF,pre2);

# 5) Replace interpolated frequency slice into final matrix

    INTDF[ia,:] = transpose(y);

end

# 6) Honor symmetries

nf1_2 = convert(Int,round(nf1/2))

INTDF[nf1_2+2:nf1,:]=conj(reverse(INTDF[2:nf1_2,:],dims=1));

# 7) Antitransform from f-x to t-x
 d_interp = real(ifft(INTDF,1));

 d_interp = d_interp[1:nt,:];

 return d_interp
end

#*******************************************************************************
"""
**INTERPOLATE_FREQ**
*fx_interpolation for reconstruction of traces*

**IN**
*   x:   spatial data to be recostruction at a given freq.
*       PF:  prediction error filter
*       pre: prewhitening in percentage

**OUT**
*  y:   reconstructed vector of data at a particular frequency
"""
function interpolate_freq(x,PF,pre)

 np = length(PF);
 nx = length(x);
 ny = 2*nx-1;       # length of interpolated signal

# Forward step using PFs for interpolation
 TMPF1 = hcat(reverse(transpose(PF),dims=2),-1)
 W1 = transpose(Convmtx(TMPF1,ny-np));

# Backward step

 TMPF2 = conj(reverse(TMPF1,dims=2));
 W2 = transpose(Convmtx(TMPF2,ny-np));

# Forward and backward combined in an augmented system

 WT = [W1;W2]

# Separation of  Matrix into  known and unknown parts
 A = WT[:,2:2:ny];

B = -1 .* WT[:,1:2:ny]*transpose(x);
# Least squares solution for missing data, prewhitening
# is added to guarantee stability

 R= A'*A;
 g = A'*B;

 mu = (pre/100.)*tr(R)/(nx-1);

y1 =  (R+mu*Matrix(I,nx-1,nx-1))\g;

 y = zeros(typeof(x[1]),1,ny);

y[1:2:ny]=x;
 y[2:2:ny]=transpose(y1);

 return y
 end

 #******************************************************************************
"""
**PREDICTION_FILTER**
* compute prediction filters

**IN**
*   VEC:   spatial data at a given frequency
*       np:    length of pef
*       pre:   prewhitening in percentage

**OUT**
*  PF:    prediction filter
"""
function prediction_filter(xi,np,pre);

 ns = length(xi);

C = zeros(typeof(xi[1]),ns-np+1,np+1);

 for j=1:ns-np

    C[j,:]=xi[j+np:-1:j];

 end

 # Make kernel for pef estimation (forward and backward prediction included)

 A = [C[:,2:np+1];conj(reverse(C[:,1:np],dims=2))];

# RHS of system of equations

 B = [C[:,1];conj(C[:,np+1])];

 # Solve A.PF = B using least-squares with pre-whitening

 R = A'*A;
 g = A'*B;

 mu = real((pre/100.)*tr(R)/np);

 PF =  (R+mu*Matrix(I,np,np))\g;   # same as inv(R)*g

 return PF
 end

#*****************************************************

function pad(d,dt,flow,fhigh;npad=1)

vect = size(d)
nt = vect[1]

#calculate index frequencies

nf = npad*nextpow(2,nt)

dpad = zeros(nf,vect[2])

dpad[1:nt,:] = d

ilow = convert(Int64,floor(flow*nf*dt)+1)
ihigh = convert(Int64,floor(fhigh*nf*dt)+1)

return dpad,ilow,ihigh
end


text:
# Cartpole
T = Float64
model = TrajectoryOptimization.Dynamics.cartpole
model_d = rk3(model)
n = model.n; m = model.m

# costs
Q = 1.0e-2*Diagonal(I,n)
Qf = 100.0*Diagonal(I,n)
R = 1.0e-1*Diagonal(I,m)
x0 = [0; 0.; 0.; 0.]
xf = [0.; pi; 0.; 0.] # (ie, swing up)

N = 101
tf = 5.
dt = tf/(N-1)
U0 = [0.01*ones(m) for k = 1:N-1]
obj = TrajectoryOptimization.LQRObjective(Q,R,Qf,xf,N)
bnd = BoundConstraint(n,m,u_min=-3.,u_max=3.)
goal = goal_constraint(xf)
constraints = Constraints(N)
for k = 1:N-1
    constraints[k] += bnd
end
constraints[N] += goal

cartpole = TrajectoryOptimization.Problem(model_d, obj, constraints=constraints,x0=x0, xf=xf, N=N, dt=dt)
initial_controls!(cartpole, U0)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
function FieldDSSTPropagator(arg0::Field, arg1::FieldODEIntegrator)
    return FieldDSSTPropagator((Field, FieldODEIntegrator), arg0, arg1)
end

function FieldDSSTPropagator(arg0::Field, arg1::FieldODEIntegrator, arg2::AttitudeProvider)
    return FieldDSSTPropagator((Field, FieldODEIntegrator, AttitudeProvider), arg0, arg1, arg2)
end

function FieldDSSTPropagator(arg0::Field, arg1::FieldODEIntegrator, arg2::PropagationType)
    return FieldDSSTPropagator((Field, FieldODEIntegrator, PropagationType), arg0, arg1, arg2)
end

function FieldDSSTPropagator(arg0::Field, arg1::FieldODEIntegrator, arg2::PropagationType, arg3::AttitudeProvider)
    return FieldDSSTPropagator((Field, FieldODEIntegrator, PropagationType, AttitudeProvider), arg0, arg1, arg2, arg3)
end

function add_additional_equations(obj::FieldAbstractIntegratedPropagator, arg0::FieldAdditionalEquations)
    return jcall(obj, "addAdditionalEquations", void, (FieldAdditionalEquations,), arg0)
end

function add_additional_state_provider(obj::FieldAbstractPropagator, arg0::FieldAdditionalStateProvider)
    return jcall(obj, "addAdditionalStateProvider", void, (FieldAdditionalStateProvider,), arg0)
end

function add_event_detector(obj::FieldAbstractIntegratedPropagator, arg0::FieldEventDetector)
    return jcall(obj, "addEventDetector", void, (FieldEventDetector,), arg0)
end

function add_force_model(obj::FieldDSSTPropagator, arg0::DSSTForceModel)
    return jcall(obj, "addForceModel", void, (DSSTForceModel,), arg0)
end

function clear_events_detectors(obj::FieldAbstractIntegratedPropagator)
    return jcall(obj, "clearEventsDetectors", void, ())
end

function compute_mean_state(obj::FieldDSSTPropagator, arg0::FieldSpacecraftState, arg1::AttitudeProvider, arg2::Collection)
    return jcall(obj, "computeMeanState", FieldSpacecraftState, (FieldSpacecraftState, AttitudeProvider, Collection), arg0, arg1, arg2)
end

function compute_mean_state(obj::FieldDSSTPropagator, arg0::FieldSpacecraftState, arg1::AttitudeProvider, arg2::Collection, arg3::jdouble, arg4::jint)
    return jcall(obj, "computeMeanState", FieldSpacecraftState, (FieldSpacecraftState, AttitudeProvider, Collection, jdouble, jint), arg0, arg1, arg2, arg3, arg4)
end

function compute_osculating_state(obj::FieldDSSTPropagator, arg0::FieldSpacecraftState, arg1::AttitudeProvider, arg2::Collection)
    return jcall(obj, "computeOsculatingState", FieldSpacecraftState, (FieldSpacecraftState, AttitudeProvider, Collection), arg0, arg1, arg2)
end

function get_additional_state_providers(obj::FieldAbstractPropagator)
    return jcall(obj, "getAdditionalStateProviders", List, ())
end

function get_all_force_models(obj::FieldDSSTPropagator)
    return jcall(obj, "getAllForceModels", List, ())
end

function get_attitude_provider(obj::FieldAbstractPropagator)
    return jcall(obj, "getAttitudeProvider", AttitudeProvider, ())
end

function get_basic_dimension(obj::FieldAbstractIntegratedPropagator)
    return jcall(obj, "getBasicDimension", jint, ())
end

function get_calls(obj::FieldAbstractIntegratedPropagator)
    return jcall(obj, "getCalls", jint, ())
end

function get_events_detectors(obj::FieldAbstractIntegratedPropagator)
    return jcall(obj, "getEventsDetectors", Collection, ())
end

function get_field(obj::FieldAbstractPropagator)
    return jcall(obj, "getField", Field, ())
end

function get_frame(obj::FieldAbstractPropagator)
    return jcall(obj, "getFrame", Frame, ())
end

function get_generated_ephemeris(obj::FieldAbstractIntegratedPropagator)
    return jcall(obj, "getGeneratedEphemeris", FieldBoundedPropagator, ())
end

function get_initial_state(obj::FieldAbstractPropagator)
    return jcall(obj, "getInitialState", FieldSpacecraftState, ())
end

function get_managed_additional_states(obj::FieldAbstractIntegratedPropagator)
    return jcall(obj, "getManagedAdditionalStates", Vector{JString}, ())
end

function get_mode(obj::FieldAbstractPropagator)
    return jcall(obj, "getMode", jint, ())
end

function get_mu(obj::FieldAbstractIntegratedPropagator)
    return jcall(obj, "getMu", RealFieldElement, ())
end

function get_orbit_type(obj::FieldDSSTPropagator)
    return jcall(obj, "getOrbitType", OrbitType, ())
end

function get_position_angle_type(obj::FieldDSSTPropagator)
    return jcall(obj, "getPositionAngleType", PositionAngle, ())
end

function get_pv_coordinates(obj::FieldAbstractPropagator, arg0::FieldAbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, (FieldAbsoluteDate, Frame), arg0, arg1)
end

function get_satellite_revolution(obj::FieldDSSTPropagator)
    return jcall(obj, "getSatelliteRevolution", jint, ())
end

function get_selected_coefficients(obj::FieldDSSTPropagator)
    return jcall(obj, "getSelectedCoefficients", Set, ())
end

function initial_is_osculating(obj::FieldDSSTPropagator)
    return jcall(obj, "initialIsOsculating", jboolean, ())
end

function is_additional_state_managed(obj::FieldAbstractIntegratedPropagator, arg0::JString)
    return jcall(obj, "isAdditionalStateManaged", jboolean, (JString,), arg0)
end

function propagate(obj::FieldAbstractIntegratedPropagator, arg0::FieldAbsoluteDate)
    return jcall(obj, "propagate", FieldSpacecraftState, (FieldAbsoluteDate,), arg0)
end

function propagate(obj::FieldAbstractIntegratedPropagator, arg0::FieldAbsoluteDate, arg1::FieldAbsoluteDate)
    return jcall(obj, "propagate", FieldSpacecraftState, (FieldAbsoluteDate, FieldAbsoluteDate), arg0, arg1)
end

function remove_force_models(obj::FieldDSSTPropagator)
    return jcall(obj, "removeForceModels", void, ())
end

function reset_initial_state(obj::FieldDSSTPropagator, arg0::FieldSpacecraftState)
    return jcall(obj, "resetInitialState", void, (FieldSpacecraftState,), arg0)
end

function set_attitude_provider(obj::FieldDSSTPropagator, arg0::AttitudeProvider)
    return jcall(obj, "setAttitudeProvider", void, (AttitudeProvider,), arg0)
end

function set_ephemeris_mode(obj::FieldAbstractIntegratedPropagator)
    return jcall(obj, "setEphemerisMode", void, ())
end

function set_initial_state(obj::FieldDSSTPropagator, arg0::FieldSpacecraftState)
    return jcall(obj, "setInitialState", void, (FieldSpacecraftState,), arg0)
end

function set_initial_state(obj::FieldDSSTPropagator, arg0::FieldSpacecraftState, arg1::PropagationType)
    return jcall(obj, "setInitialState", void, (FieldSpacecraftState, PropagationType), arg0, arg1)
end

function set_interpolation_grid_to_fixed_number_of_points(obj::FieldDSSTPropagator, arg0::jint)
    return jcall(obj, "setInterpolationGridToFixedNumberOfPoints", void, (jint,), arg0)
end

function set_interpolation_grid_to_max_time_gap(obj::FieldDSSTPropagator, arg0::RealFieldElement)
    return jcall(obj, "setInterpolationGridToMaxTimeGap", void, (RealFieldElement,), arg0)
end

function set_master_mode(obj::FieldAbstractIntegratedPropagator, arg0::FieldOrekitStepHandler)
    return jcall(obj, "setMasterMode", void, (FieldOrekitStepHandler,), arg0)
end

function set_master_mode(obj::FieldAbstractPropagator, arg0::RealFieldElement, arg1::FieldOrekitFixedStepHandler)
    return jcall(obj, "setMasterMode", void, (RealFieldElement, FieldOrekitFixedStepHandler), arg0, arg1)
end

function set_mu(obj::FieldDSSTPropagator, arg0::RealFieldElement)
    return jcall(obj, "setMu", void, (RealFieldElement,), arg0)
end

function set_reset_at_end(obj::FieldAbstractIntegratedPropagator, arg0::jboolean)
    return jcall(obj, "setResetAtEnd", void, (jboolean,), arg0)
end

function set_satellite_revolution(obj::FieldDSSTPropagator, arg0::jint)
    return jcall(obj, "setSatelliteRevolution", void, (jint,), arg0)
end

function set_selected_coefficients(obj::FieldDSSTPropagator, arg0::Set)
    return jcall(obj, "setSelectedCoefficients", void, (Set,), arg0)
end

function set_slave_mode(obj::FieldAbstractIntegratedPropagator)
    return jcall(obj, "setSlaveMode", void, ())
end

function tolerances(::Type{FieldDSSTPropagator}, arg0::RealFieldElement, arg1::FieldOrbit)
    return jcall(FieldDSSTPropagator, "tolerances", Vector{Vector{jdouble}}, (RealFieldElement, FieldOrbit), arg0, arg1)
end


text:
# ------------------------------------------------------------------------------
# diffusion operator example
# ------------------------------------------------------------------------------

using LinearAlgebra

# setup
mesh = Mesh2D(1.0, 1.0)
basis = TensorH1LagrangeBasis(4, 4, 1, 2)

# weak form
function diffusionweakform(du::Array{Float64}, w::Array{Float64})
    dv = du * w[1]
    return [dv]
end

# diffusion operator
inputs = [
    OperatorField(basis, [EvaluationMode.gradient]),
    OperatorField(basis, [EvaluationMode.quadratureweights]),
]
outputs = [OperatorField(basis, [EvaluationMode.gradient])]
diffusion = Operator(diffusionweakform, mesh, inputs, outputs)

# compute operator symbols
A = computesymbols(diffusion, [π, π])
eigenvalues = real(eigvals(A))

# ------------------------------------------------------------------------------


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_input(id="input", value = "initial value"),
    html_div(id="output")
end
callback!(app, callid"input.value => output.children") do input
    return html_div() do
        dcc_input(id = "sub-input-1", value = "sub input initial value"),
        html_div(id = "sub-output-1")
    end
end

callback!(app, callid"sub-input-1.value => sub-output-1.children") do value
    return value
end

run_server(app)


text:
# ---
# title: 56. Merge Intervals
# id: problem56
# author: Jimmy Shen
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Sort
# link: <https://leetcode.com/problems/merge-intervals/description/>
# hidden: true
# ---
#
# Given an array of `intervals` where `intervals[i] = [starti, endi]`, merge all
# overlapping intervals, and return _an array of the non-overlapping intervals
# that cover all the intervals in the input_.
#
#
#
# **Example 1:**
#
#
#
#     Input: intervals = [[1,3],[2,6],[8,10],[15,18]]
#     Output: [[1,6],[8,10],[15,18]]
#     Explanation: Since intervals [1,3] and [2,6] overlaps, merge them into [1,6].
#
#
# **Example 2:**
#
#
#
#     Input: intervals = [[1,4],[4,5]]
#     Output: [[1,5]]
#     Explanation: Intervals [1,4] and [4,5] are considered overlapping.
#
#
#
#
# **Constraints:**
#
#   * `1 <= intervals.length <= 104`
#   * `intervals[i].length == 2`
#   * `0 <= starti <= endi <= 104`
#
#
## @lc code=start
using LeetCode

function merge_intervals(intervals::Vector{Vector{Int}})::Vector{Vector{Int}}
    result = Vector{Vector{Int}}([])
    for (a, b) in sort!(intervals, by = i -> i[1])
        if size(result, 1) > 0 && result[end][2] >= a
            result[end][2] = max(result[end][2], b)
        else
            push!(result, [a, b])
        end
    end
    return result
end
## @lc code=end


text:
mutable struct PairwiseTTestResult<:TestResult
    levels::Vector{String}
    correction::Symbol
    n_obs::IntVector
    means::FloatVector
    sds::FloatVector
    statistic::FloatVector
    dof::FloatVector
    pval::FloatVector
    cohen_d::FloatVector
end

function show(io::IO, pttr::PairwiseTTestResult)
end

function coeftable(pttr::PairwiseTTestResult)
end


text:
"""
	$(SIGNATURES)

Constructor for `EndowmentDraws`.
"""
EndowmentDraws() = EndowmentDraws(Vector{Endowment}(), DataFrame())

# EndowmentDraws(meta :: Vector{Endowment{Any}}) =
#     EndowmentDraws(meta, DataFrame());

# EndowmentDraws(meta :: Vector{Endowment{T1}}) where T1 =
#     EndowmentDraws(meta, DataFrame());

"""
	$(SIGNATURES)

Validate endowment draws. Includes bounds checking.
"""
function validate_draws(ed :: EndowmentDraws)
    isValid = true;
    if !isempty(ed)
        for dName in names(ed)
            isValid = isValid  &&
                validate_draws(get_meta(ed, dName), get_draws(ed, dName));
        end
    end
    return isValid
end

"""
	$(SIGNATURES)

Compute the correlation matrix of select endowments. Defaults to all endowments in the order determined by `names(ed)`.
Note that not all endowments are scalar. Correlations are only computed for those that are. The others are set to NaN.
"""
function corr_matrix(ed :: EndowmentDraws, nameV :: Vector{Symbol} = names(ed))
    n = length(nameV);
    @assert n > 1
    corrM = zeros(n, n);
    for i1 = 1 : (n-1)
        d1 = get_draws(ed, nameV[i1]);
        corrM[i1, i1] = 1.0;
        if isa(d1, Vector{<:Real})
            for i2 = (i1+1) : n
                d2 = get_draws(ed, nameV[i2]);
                if isa(d2, Vector{<:Real})
                    corrM[i1, i2] = cor(d1, d2);
                else
                    corrM[i1, i2] = NaN;
                end
            end
        end
    end
    corrM[n, n] = 1.0;
    return corrM
end

"""
	$(SIGNATURES)

Create formatted correlation matrix for selected endowments. Returns a string matrix with headers.
"""
function formatted_corr_matrix(ed :: EndowmentDraws,
    nameV :: Vector{Symbol} = names(ed))

    corrM = corr_matrix(ed, nameV);
    return formatted_corr_matrix(corrM, nameV);
end

"""
	$(SIGNATURES)

Create formatted correlation matrix from a numeric matrix and variable names. Correlation matrix may contain NaNs.
"""
function formatted_corr_matrix(corrM :: AbstractMatrix{F1}, nameV :: Vector{T2}) where {F1 <: Real, T2}

    n = length(nameV);
    @assert size(corrM) == (n, n)

    m = fill("", n+1, n+1);
    for j = 1 : n
        m[1, j+1] = string(nameV[j]);
        m[j+1, 1] = string(nameV[j]);
        m[j+1, j+1] = "1.0";
    end
    for j1 = 1 : (n-1)
        for j2 = (j1 + 1) : n
            c = corrM[j1, j2];
            if isnan(c)
                s = "--";
            else
                s = "$(round(c, digits = 2))";
            end
            m[j1+1, j2+1] = s;
        end
    end
    return m
end

## -------------  Access

Base.show(io :: IO, ed :: EndowmentDraws) =
    print(io,  "EndowmentDraws of length ", length(ed), " with fields ", names(ed));
Base.isempty(ed :: EndowmentDraws) = (length(ed) == 0);
Base.length(ed :: EndowmentDraws) = size(ed.draws, 1);

# Returns names as Vector{Symbol}
function Base.names(ed :: EndowmentDraws)
    if isempty(ed)
        return nothing
    else
        return Symbol.(names(ed.draws));
    end
end

"""
	$(SIGNATURES)

Retrieve endowment draws for one type as a DataFrame row. One can splat this as a tuple into a function:

`foo(type_endowments(e, j)...)` is the same as `foo(x, y)`
"""
type_endowments(e :: EndowmentDraws, j :: Integer) = e.draws[j,:];

"""
	$(SIGNATURES)

Create a copy of an `EndowmentDraws` object with only select cases.
"""
select_rows(e :: EndowmentDraws, idxV) = EndowmentDraws(e.meta, e.draws[idxV, :]);

# """
# 	$(SIGNATURES)

# Select cases in place.
# """
# select_rows!(e :: EndowmentDraws, idxV) = (e.draws = e.draws[idxV, :]);

Base.eltype(e :: EndowmentDraws, eName :: Symbol) =
    eltype(get_meta(e, eName));

"""
	$(SIGNATURES)

Retrieve label for one endowment.
"""
get_label(e :: EndowmentDraws, eName :: Symbol) =
    label(get_meta(e, eName));

"""
	$(SIGNATURES)

Retrieve the meta information about one `Endowment`. Returns an `Endowment` object. Nothing if not found.
"""
function get_meta(e :: EndowmentDraws, eName :: Symbol)
    idx = find_meta(e, eName);
    if !isnothing(idx)
        return e.meta[idx];
    else
        return nothing
    end
end

# Find index into `meta` for an endowment. Nothing if not found.
function find_meta(e :: EndowmentDraws, eName :: Symbol)
    if isempty(e)
        return nothing
    else
        return findfirst(n -> eName == n,  names(e));
    end
end

has_endowment(e :: EndowmentDraws, eName :: Symbol) =
    !isnothing(find_meta(e, eName));

"""
	$(SIGNATURES)

Retrieve draws for one endowment; for all individuals.
"""
get_draws(e :: EndowmentDraws, eName) = e.draws[!, eName];

get_draws(e :: EndowmentDraws, eName, idxV) = e.draws[idxV, eName];

"""
	$(SIGNATURES)

Return DataFrame with all endowments. Useful for running regressions on endowments.
"""
get_dataframe(e :: EndowmentDraws) = e.draws;

"""
	$(SIGNATURES)

Returns quantiles of one endowment from its draws. Returns the same `DataType` as the endowment itself (e.g. nearest Integer).

Quantile = quantile from Statistics, rounded to the `eltype` of the endowment draw.

Only for scalar, numeric endowments.
"""
function endow_quantiles(e :: EndowmentDraws, eName, pctV)
    eType = eltype(e, eName);
    @assert eType <: Real  "Only for scalar numeric endowments. $eName is of $eType"
    qV = quantile(get_draws(e, eName), pctV);
    if eType <: Integer
        eqV = round.(eType, qV);
    else
        eqV = convert.(eType, qV);
    end
    return eqV
end

## ---------  Modify

"""
	$(SIGNATURES)

Add draws for one variable to existing `EndowDraws`.
"""
function add_draws!(e :: EndowmentDraws, endow :: Endowment{T1},
    dV :: AbstractVector{T1}) where T1

    @assert !has_endowment(e, name(endow))
    push!(e.meta, endow);
    setproperty!(e.draws, name(endow), dV);
end

"""
	$(SIGNATURES)

Replace values for one endowment. Type cannot change.
`newVals` is a Vector. One element for each draw. Or it is a scalar.
"""
function replace_draws!(e :: EndowmentDraws, eName :: Symbol, newVals :: T1) where T1
    oldType = eltype(getproperty(e.draws, eName));
    if oldType == T1
        # This is necessary in case each entry is a Vector
        for j = 1 : length(e)
            e.draws[j, eName] = newVals;
        end
    elseif eltype(newVals) == oldType
        # Input is a vector. One draw for each old draw.
        @assert length(newVals) == length(e)
        e.draws[!, eName] .= newVals;
    else
        error("""
        Type mismatch for $eName:
        $(typeof(newVals))  vs  $(oldType)
        """);
    end
end

# """
# 	$(SIGNATURES)

# Replace endowments with their quantiles.
# """
# function replace_with_quantiles!(e :: EndowmentDraws, eNames,
#     pctV :: AbstractVector{F1}) where F1 <: AbstractFloat

#     @assert length(pctV) == length(e)  "Need one percentile for each observation";
#     for eName in eNames
#         replace_draws!(e, eName, endow_quantiles(e, eName, pctV));
#     end
# end

"""
	$(SIGNATURES)

Return `n` endowment draws where all numeric endowments are set to the same quantile in their distributions. Non-scalar endowments are set arbitrarily (because they don't have well-defined quantiles).

# Arguments
- `addNoise`: scale for additive noise. If not `nothing`, create a Vector in (0, addNoise) and add it to the draws. The purpose is to leave some variation in the draws for cases where identical draws cause numerical issues.
"""
function draw_fixed_percentiles(e :: EndowmentDraws, n :: Integer, pct :: F1;
    addNoise :: F2 = 0.0) where {F1 <: AbstractFloat, F2 <: AbstractFloat}

    draws = select_rows(e, 1 : n);
    for eName in names(e)
        eType = eltype(e, eName);
        if eType <: Real
            q = endow_quantiles(e, eName, pct);
            replace_draws!(draws, eName, q .+ gen_noise(eType, addNoise, n));
        end
    end
    return draws
end

function gen_noise(F1 :: Type{<:AbstractFloat}, addNoise, n :: Integer)
    if addNoise > zero(F1)
        nV = F1(addNoise) .* rand(F1, n);
    else
        nV = zero(F1);
    end
end

# For Integers, cannot add noise
gen_noise(I1 :: Type{<:Integer}, addNoise, n :: Integer) = zero(I1);

# -------  For testing

function make_test_endowment_draws(n :: Integer)
    rng = MersenneTwister(43);
    ev = make_test_endowment_vector();
    ed = EndowmentDraws();
    for endow in ev
        dV = draw_test_endowments(endow, n, rng);
        add_draws!(ed, endow, dV);
    end
    @assert validate_draws(ed)
    return ed
end

# -----------------


text:
@testset "Group properties" begin

    σ1 = Permutation([1, 2, 3])
    σ2 = Permutation([2, 3, 1])
    σ3 = Permutation([3, 1, 2])
    τ1 = Permutation([1, 3, 2])
    τ2 = Permutation([3, 2, 1])
    τ3 = Permutation([2, 1, 3])
    cycles = [σ1, σ2, σ3]
    transpositions = [τ1, τ2, τ3]

    @testset "Existence of identity" begin
        for π = [cycles..., transpositions...]
            @test σ1 * π == π * σ1 == π
        end
    end

    @testset "Existence of the inverse" begin
        @test σ2 * σ3 == σ3 * σ2 == σ1
        for τ = transpositions
            @test τ * τ == σ1
        end
    end

    @testset "Closeness under multiplication" begin
        @test σ2 * σ2 == σ3
        @test σ3 * σ3 == σ2
        @test τ1 * τ2 == τ2 * τ3 == τ3 * τ1 == σ2
        @test τ3 * τ2 == τ2 * τ1 == τ1 * τ3 == σ3
    end

    @testset "Associativity" begin
        p_group = [cycles..., transpositions...]
        for π1 = p_group, π2 = p_group, π3 = p_group
            @test (π1 * π2) * π3 == π1 * (π2 * π3)
        end
    end

end


text:
using Sundials, Test, SparseArrays

# Test for Jacobian usage
function Lotka(du, u, p, t)
  du[1] = u[1] - u[1] * u[2] # REPL[7], line 3:
  du[2] = -3 * u[2] + 1 * u[1] * u[2]
  nothing
end

jac_called = false
function Lotka_jac(J,u,p,t)
  global jac_called
  jac_called = true
  J[1,1] = 1.0 - u[2]
  J[1,2] = -u[1]
  J[2,1] = 1 * u[2]
  J[2,2] = -3 + u[1]
  nothing
end

Lotka_f = ODEFunction(Lotka,jac=Lotka_jac)
prob = ODEProblem(Lotka_f,ones(2),(0.0,10.0))
good_sol = solve(prob,CVODE_BDF())
@test jac_called == true

Lotka_f = ODEFunction(Lotka,jac=Lotka_jac,
                      jac_prototype = sparse([1,2,1,2],[1,1,2,2],zeros(4)))

prob = ODEProblem(Lotka_f,ones(2),(0.0,10.0))
jac_called = false
sol9 = solve(prob,CVODE_BDF(linear_solver=:KLU))
@test jac_called == true

function f2!(res, du, u, p, t)
    res[1] = 1.01du[1]
    return
end

jac_called = false
function f2_jac!(out, du, u, p, gamma, t)
    global jac_called
    jac_called = true
    out[1] = 1.01
end

f2_f = DAEFunction(f2!,jac=f2_jac!)

u0 = [0.]
tspan = (0.0, 10.)
du0 = [0.]
dae_prob = DAEProblem(f2_f,du0, u0,tspan, differential_vars=[true])
good_sol = solve(dae_prob,IDA())
@test jac_called == true

function testjac(res, du, u, p, t)
  res[1] = du[1] - 1.5 * u[1] + 1.0 * u[1]*u[2]
  res[2] = du[2] +3 * u[2] - u[1]*u[2]
end

jac_called = false
function testjac_jac(J, du, u, p, gamma, t)
  global jac_called
  jac_called = true
  J[1,1] = gamma - 1.5 + 1.0 * u[2]
  J[1,2] = 1.0 * u[1]
  J[2,1] = - 1 * u[2]
  J[2,2] = gamma + 3 - u[1]
  nothing
end

testjac_f = DAEFunction(testjac,jac=testjac_jac)

prob3 = DAEProblem(testjac_f,[0.5,-2.0],ones(2),(0.0,10.0),differential_vars=[true,true])
sol3 = solve(prob3, IDA())

@test jac_called == true
jac_called = false
prob4 = DAEProblem(testjac,[0.5,-2.0],ones(2),(0.0,10.0))
sol4 = solve(prob4, IDA())

@test jac_called == false

println("Jacobian vs no Jacobian difference:")
println(maximum(sol3-sol4))
@test maximum(sol3-sol4) < 1e-6


text:
using PerformanceRNN
using Test

max_stars_repo_name: VasanthManiVasi/PerformanceRNN

max_issues_repo_name: VasanthManiVasi/PerformanceRNN

max_forks_repo_name: VasanthManiVasi/PerformanceRNN


text:
## Module associated with the MeshArrays package

module MeshArrays

include("Types.jl");
include("Grids.jl");
include("Operations.jl");
include("Exchanges.jl");
include("ReadWrite.jl");
include("Demos.jl");
include("Solvers.jl");

export AbstractMeshArray, MeshArray, gcmgrid
export exchange, gradient, convergence, smooth, mask
export GridSpec, GridLoad, GridOfOnes
export findtiles, LatitudeCircles, ThroughFlow
export ScalarPotential, VectorPotential
#The following exch_UV differs from normal exchange; incl. exch_UV_N.
export exch_UV

#The following codes add dependencies to Plots
#using Plots; include(joinpath(dirname(pathof(MeshArrays)),"Plots.jl"));

#The following codes add dependencies to NetCDF
#include("gcmfaces_nctiles.jl");

#The following codes add being deprecated
#include("gcmfaces_plot.jl");
#include("gcmfaces_convert.jl");

end # module


text:
mutable struct DRAWOutput
    mus
    logsigmas
    sigmas
    cs
end

DRAWOutput() = DRAWOutput([], [], [], [])

function push!(o::DRAWOutput, c)
    push!(o.cs, c)
end

function push!(o::DRAWOutput, mu, logsigma, sigma, c)
    push!(o.mus, mu)
    push!(o.logsigmas, logsigma)
    push!(o.sigmas, sigma)
    push!(o.cs, c)
end

function empty!(o::DRAWOutput)
    empty!(o.mus, mu)
    empty!(o.logsigmas, logsigma)
    empty!(o.sigma, sigma)
    empty!(o.cs, c)
end

function report(epoch, trn, tst=nothing)
    print("[$(now())] epoch=$epoch, trn=$trn")
    if tst != nothing; println(", tst=$tst")
    else; println(); end
    flush(stdout)
end

max_stars_repo_name: ilkerkesen/draw

max_issues_repo_name: ilkerkesen/draw

max_forks_repo_name: ilkerkesen/draw


text:
# Autogenerated wrapper script for Attr_jll for aarch64-linux-musl
export attr

JLLWrappers.@generate_wrapper_header("Attr")
JLLWrappers.@declare_library_product(attr, "libattr.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        attr,
        "lib/libattr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using DataDeps

register(DataDep(
  "dermatology",
  "http://archive.ics.uci.edu/ml/datasets/Dermatology",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/dermatology/dermatology.data",
  "455eba77f72cd087ce54a5a266c514c5f34e85000c8689ad09e796d68ad45742",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=35,
        feature_indices=1:34,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
"""
    assigndist_beta(α, β; trunc_lower = -Inf, trunc_upper = Inf)

Assign parameters to a Beta distribution with parameters `α` and `β`, optionally
truncating the distribution.
"""
function assigndist_beta(α, β; trunc_lower = -Inf, trunc_upper = Inf)
    Truncated(Beta(α, β), trunc_lower, trunc_upper)
end

export assigndist_beta


text:
using Tries
function Either(x::Vector{<:AbstractString})
    P = Trie{Char,Nothing}
    r = P()
    for e in x
        r[e...] = nothing
    end
    Either{String}(r)
end

function Either(x::Dict)
    P = Trie{Char,valtype(x)}
    r=P()
    for (e,v) in pairs(x)
        r[e...] = v
    end
    Either{valtype(x)}(r)
end
either_state_type(T::Type{<:Trie}) = NCodeunitsState{T}

"""
    _iterate(p::Trie{Char}, str, till, posi, next_i, ::Nothing)

Match char path in `p` greedily, recording shorter matches in state.
"""
@inline function _iterate(p::Union{Trie{Char},SubTrie{Char}}, str, till, posi, next_i, state::Nothing)
    ni = ni_ = posi
    st = st_ = p
    while st !== nothing && ni <= till
        @inbounds c, ni = iterate(str,ni)
        @inbounds st = get( Tries.nodes(st),c, nothing)
        st === nothing && break
        if get(st) !== missing
            ni_ = ni
            st_ = st
        end
    end
    if get(st_) !== missing
        return ni_, NCodeunitsState(ni_-posi,st_)
    else
        return nothing
    end
end

@inline function _iterate(p::AbstractTrie{Char}, str, till, posi, next_i, state)
    _iterate(p, str, prevind(str,next_i,2), posi, posi, nothing)
end

function _iterate(p::Either{<:AbstractTrie}, str, till, posi, next_i, state)
    _iterate(p.options, str, till, posi, next_i, state)
end

@inline nextind(str,i::Int,parser::Either{<:AbstractTrie},x::NCodeunitsState)  = i+x.nc
@inline prevind(str,i::Int,parser::Either{<:AbstractTrie},x::NCodeunitsState)  = i-x.nc

children(x::Either{<:AbstractTrie}) =
    children(x.options)

function deepmap_parser(f::typeof(lowercase),mem::AbstractDict,x::Either{<:AbstractTrie},a...;kw...)
    get!(mem,x) do
        Either{result_type(x)}(Trie(lowercase.(k)=>v for (k,v) in pairs(x.options) if !isempty(k)))
    end
end

# struct InternedParser{E, I<:Either{<:AbstractTrie}} <: WrappedParser{E,NCodeunitsState,T} where {P}
#     interned::I
#     parser::E
# end


text:
@testset "Redundancy removal" begin
    A = [9.69462e-5 9.69462e-5 9.69462e-5 9.69462e-5 0.0 0.0 0.0 0.0]
    B = [8.95335e-5 8.95335e-5 8.95335e-5 8.95335e-5 0.0 0.0 0.0 0.0]
    a = [1.]
    @test CutPruners.checkredundancy(A, a, B, a, true, true, 1e-5) == [1]
    @test CutPruners.checkredundancy(A, a, B, a, true, false, 1e-5) == [1]
    @test CutPruners.checkredundancy(A, a, B, a, true, true, 1e-6) == Int[]
    @test CutPruners.checkredundancy(A, a, B, a, true, false, 1e-6) == Int[]
    # A and B will be renormalized so A would be very close to B and even 1e-16 works
    @test CutPruners.checkredundancy(A, a, B, a, false, true, 1e-16) == Int[]
    @test CutPruners.checkredundancy(A, a, B, a, false, false, 1e-16) == [1]
end
@testset "Self redundancy" begin
    A = [9.69462e-5 9.69462e-5;
         8.95335e-5 8.95335e-5;
         0.0 0.0; 0.0 0.0]
    a = ones(Float64, 4)
    @test CutPruners.checkredundancy(A, a, true, true, 1e-5) == [2, 4]
    @test CutPruners.checkredundancy(A, a, true, false, 1e-5) == [2, 4]
    @test CutPruners.checkredundancy(A, a, true, true, 1e-6) == Int[4]
    @test CutPruners.checkredundancy(A, a, true, false, 1e-6) == Int[4]
    # A and B will be renormalized so A would be very close to B and even 1e-16 works
    @test CutPruners.checkredundancy(A, a, false, true, 1e-16) == Int[4]
    @test CutPruners.checkredundancy(A, a, false, false, 1e-6) == [2, 4]
end


text:
using Test
using RingLists

a = RingList(1,2,3)
b = RingList(2,3,1)
@test a==b
@test reverse(a)==reverse(b)
@test hash(a) == hash(b)
@test length(a) == 3
@test sort(collect(keys(a))) == [1,2,3]
@test Set(a) == Set([1,2,3])
@test first(a,true) == 1
@test a[1] == 2
@test Vector(b) == [1,2,3]
@test eltype(b) == Int
a = RingList{Int}()
insert!(a,1)
insert!(a,2)
@test Vector(a) == [1,2]

a = RingList(1,2,3,4)
b = RingList(1,2,4)
delete!(a,3)
@test a==b


text:
# # Chapter 1

# ### Inspect the generated code for the power function.

# For integers:
@code_native 3^2

# And for floats
@code_native 3.5^2


text:
using ArgParse2
using DataFrames
using DataStructures
using BioSequences
using CodecZlib
using BioAlignments
using DataFramesMeta
using Statistics
using st16SseqJuliaTools
using XAM
using BioTools.BLAST
using CSV
using Artifacts

function main()
    parser = ArgumentParser(prog = "analyse_alnself.jl",
                        description = "Detect chimeras",
                        epilog = "Removing potential chimeras and repair snps")
    add_argument!(parser, "--inxml", "-x", help = "input blast xml ",type = String)
    add_argument!(parser, "--ref", "-r", help = "reference fasta",type = String)
    add_argument!(parser, "--out", "-o", help = "Fasta file to write out cleaned up sequences",default="chimerasfree.faSTA",type = String)
    add_argument!(parser, "--min-seqid", "-s", help = "Minimum sequence identity in match",type = Float64, default=0.95)
    add_argument!(parser, "--min-covid", "-c", help = "Minimum sequence length fraction in  match",type = Float64, default=0.96)
    #add_argument!(parser, "--output", "-o", help = "File with names",type = String)

    maxNcount = 5 #max N counts in ref
    max_match_start_pos = 5

    args = parse_args(parser)
    ref = fasta_to_dict(args.ref)
    xmldata = open(args.inxml) do file
        read(file, String)
    end
    records = split_by_query(readblastXML2(xmldata))
    testid = "2_3209;size=1"

    #Limits for real match ... mayve should be options
    minf_len = args.min_covid
    minf_id = args.min_seqid # Identity of first 250 fragments of contigs
    pairs = Array{Tuple{String,String},1}()
    chosenseqs = OrderedDict{String,BioSequences.LongSequence{BioSequences.DNAAlphabet{4}}}()
    ct = 0
    info = DataFrame(contig = Array{String,1}(),
    coverage=Array{Float64,1}(),
    identity=Array{Float64,1}(),
    chimeric=Array{Bool,1}(),
    gaps=Array{Int64,1}(),
    max_indel_length=Array{Int64,1}(),
    good=Array{Bool,1}(),
    refNcounts=Array{Int64,1}(),
    match_from=Array{Int64,1}(),
    match_to=Array{Int64,1}(),
    totall_length=Array{Int64,1}(),
    match_taxid=Array{String,1}())
    for contig_id in keys(records)
        ct+=1
        matching_hits = get_coverage_maximising_hits(contig_id,records)
        r = matching_hits[1]
        contig_seq = ref[contig_id]
        gaps = r.gaps
        aln = r.alignment.aln
        count_amb = count(isambiguous, r.hit)
        if testid == contig_id
            println(r)
            println(r.hit)
            println(count_amb)
            for h in matching_hits
                println(h)
            end
        end
        s = first(aln.anchors).seqpos
        e = last(aln.anchors).seqpos
        l=abs(e-s)
        #count count_deletions
        indel_lengths = Array{Int64,1}()
        push!(indel_lengths,0)
        seqpos, refpos
        if length(aln.anchors) > 1
            for i in 2:length(aln.anchors)
                anc = aln.anchors[i]
                anb = aln.anchors[i-1]
                op = anc.op
                if isdeleteop(op) || isinsertop(op)
                    idl = maximum([abs(anc.refpos - anb.refpos),abs(anc.seqpos - anb.seqpos)])
                    push!(indel_lengths,idl)
                end
            end
        end
        contig_length = length(contig_seq)
        covf = round(l/contig_length ,digits=2)
        idf = round(r.identity/l, digits = 2)
        isgood = false
        ischimeric = (length(matching_hits) >1)
        push!(info,[contig_id,covf,idf,ischimeric,gaps,maximum(indel_lengths), isgood, count_amb, r.queryfrom, r.queryto,contig_length,r.hittaxid])
    end
    println("Matching genus to taxids...")
    species = taxid_to_species(info.match_taxid,artifact"taxdump")
    println("...done")
    info.genus = map(x -> split(x)[1],species)
    info.species_wog = map(x -> join(split(x),"_"),species)
    info.species = species
    info.cluster = map(x -> parse(Int64,split(x,"_")[1]),info.contig)
    info.matchlen = abs.(info.match_to .- info.match_from)
    sort!(info,[:cluster,:genus,:species])
    #Prefilter based on id and coverage for search of a typical species
    infoF  = @where(info, :coverage .>= 0.90 , :identity .>= 0.96, :match_from .<=max_match_start_pos ,:refNcounts.<=maxNcount)
    sub = @select(infoF, :cluster, :genus, :contig,:matchlen)
    cl_s = Dict{Int64,String}() #expected genus per cluster
    gd = groupby(sub, :cluster)
    for cl in unique(sub.cluster)
        ginfo=get(gd, (cluster=cl,), nothing)
        sinfog= groupby(ginfo, :genus)
        s_vs_l = @combine(sinfog, Length = median(:matchlen))
        sort!(s_vs_l,[:Length],rev=true)
        s = s_vs_l.genus[1]
        cl_s[cl] = s
    end

    expected_genus = map(x -> cl_s[x],info.cluster)
    info.Typical_cluster_genus = (info.genus .== expected_genus)

    #test = @where(info, map(x -> occursin("8_246",x),:contig))
    #println(test)
    #exit()

    #good  = @where(info, :coverage .> 0.96 , :identity .> 0.98, .!:chimeric).contig
    good  = @where(info, :gaps .<= 5,:max_indel_length .<=2 , :coverage .>= 0.90 , :identity .>= 0.96, :match_from .<=max_match_start_pos ,:refNcounts.<=maxNcount,:Typical_cluster_genus).contig
    info.good = map(x->(x in good),info.contig)
    CSV.write(args.out*"_blastanalysis.csv",info)
    #collect info about end positions fro write out
    ends4writeout = Dict(zip(info.contig,info.match_to))
    species4writeout = Dict(zip(info.contig,info.species_wog))
    ct_good = 0
    for k in keys(ref)
        if k  in good
            ct_good += 1
            ko = k*"@"*species4writeout[k]
            chosenseqs[ko] = ref[k][1:ends4writeout[k]]
        end
    end
    ct_all = length(keys(ref))

    println("$ct_all in total sequences, $ct_good sequences were  written to ",args.out)
    dict_to_fasta(chosenseqs,args.out)
end

main()


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

name = "HepMC3"
version = v"3.2.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://hepmc.web.cern.ch/hepmc/releases/HepMC3-3.2.2.tar.gz", "0e8cb4f78f804e38f7d29875db66f65e4c77896749d723548cc70fb7965e2d41")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/HepMC3-*/
mkdir build && cd build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DHEPMC3_ENABLE_ROOTIO:BOOL=OFF -DHEPMC3_ENABLE_TEST:BOOL=ON -DHEPMC3_INSTALL_INTERFACES:BOOL=ON -DHEPMC3_ENABLE_PYTHON:BOOL=OFF -DHEPMC3_BUILD_STATIC_LIBS:BOOL=OFF -DHEPMC3_BUILD_DOCS:BOOL=OFF ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("armv7l", "linux"; libc="glibc"),
    Platform("powerpc64le", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="musl"),
    Platform("aarch64", "linux"; libc="musl"),
    Platform("armv7l", "linux"; libc="musl"),
    Platform("x86_64", "linux"; libc="musl"),
    Platform("x86_64", "macos";)
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libHepMC3search", :libHepMC3search),
    LibraryProduct("libHepMC3", :libHepMC3)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; 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


text:
using PyPlot
using DelimitedFiles

if !isdir("figures_summary")
  mkdir("figures_summary")
end

m = 15
n = 30
h = 30.0 # meter
dt = 0.00025
nt = 3001

rc("axes", titlesize=16)
rc("axes", labelsize=12)
rc("xtick", labelsize=12)
rc("ytick", labelsize=12)
rc("legend", fontsize=16)

shot1=read("CO2/Data1/Shot8.bin");
shot1 = reshape(reinterpret(Float32,shot1),(nt,142))
fig,ax = subplots()
imshow(shot1, extent=[0,m*h,(nt-1)*dt,0], cmap="gray", aspect=1.5*(m*h)/((nt-1)*dt));
# imshow(shot1', extent=[0,(nt-1)*dt,m*h,0], cmap="gray", aspect=0.8*(nt-1)*dt/(m*h));
xlabel("Depth (m)")
ylabel("Time (s)")
ax.xaxis.tick_top()
ax.xaxis.set_label_position("top")
savefig("figures_summary/CO2_Data1_Shot8.pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

shot2=read("CO2/Data11/Shot8.bin");
shot2 = reshape(reinterpret(Float32,shot2),(nt,142))
fig,ax = subplots()
imshow(shot2, extent=[0,m*h,(nt-1)*dt,0], cmap="gray", aspect=1.5*(m*h)/((nt-1)*dt));
xlabel("Depth (m)")
ylabel("Time (s)")
ax.xaxis.tick_top()
ax.xaxis.set_label_position("top")
savefig("figures_summary/CO2_Data11_Shot8.pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

fig,ax = subplots()
imshow(shot2-shot1, extent=[0,m*h,(nt-1)*dt,0], cmap="gray", aspect=1.5*(m*h)/((nt-1)*dt));
xlabel("Depth (m)")
ylabel("Time (s)")
ax.xaxis.tick_top()
ax.xaxis.set_label_position("top")
savefig("figures_summary/CO2_Data_diff.pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);


text:
using Test
using LinearAlgebra
using SparseArrays
using CartesianMesh
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using HDGElasticity

function traction_boundary_classifier(x)
    if x[2] ≈ 0.0 || x[1] ≈ 0.0
        return :mixed
    else
        return :traction
    end
end

function displacement_boundary_classifier(x)
    if x[2] ≈ 0.0 || x[1] ≈ 0.0
        return :mixed
    elseif x[1] ≈ 1.0
        return :displacement
    else
        return :traction
    end
end

function mixed_bc_components(x)
    if x[2] ≈ 0.0
        return [0.0,1.0],[1.,0.]
    elseif x[1] ≈ 0.0
        return [1.,0.],[0.,1.]
    end
end

function mixed_bc_values(x)
    return 0.0,0.0
end

function traction_bc_data(x,xR)
    if x[1] ≈ xR
        return [1.0,0.0]
    else
        return [0.0,0.0]
    end
end

function displacement_bc_data(x,u)
    return [u[1],x[2]*u[2]]
end

function allequal(u,v)
    return all(u .== v)
end

function allapprox(u,v)
    return all(u .≈ v)
end

function allapprox(v1, v2, tol)
    flag = length(v1) == length(v2)
    return flag && all([isapprox(v1[i],v2[i],atol=tol) for i = 1:length(v1)])
end

function distance_function(coords,xc)
    return coords[1,:] .- xc
end

polyorder = 1
numqp = 2
levelset = InterpolatingPolynomial(1,2,polyorder)
mesh = UniformMesh([0.,0.],[1.,1.],[1,1])
interface_location = 2.5
levelsetcoeffs = HDGElasticity.levelset_coefficients(
    x->distance_function(x,interface_location),mesh,levelset.basis
)
dgmesh = HDGElasticity.DGMesh(mesh,levelsetcoeffs,levelset)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,polyorder,numqp,
    levelsetcoeffs,levelset)
lambda1,mu1 = 1.,2.
lambda2,mu2 = 1.,2.
stabilization = 0.1
D1 = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(lambda1,mu1))
D2 = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(lambda2,mu2))

e11 = 1.0/((lambda2+2mu2) - lambda2^2/(lambda2+2mu2))
e22 = -lambda2/(lambda2+2mu2)*e11
u1 = e11*1
u2 = e22*1

cellsolvers = HDGElasticity.CellSolvers(dgmesh,ufs,D1,D2,stabilization)
hybrid_element_numbering = HDGElasticity.HybridElementNumbering(dgmesh,ufs)
facetohelid = hybrid_element_numbering.facetohelid

bc_data = Dict(:mixed_components => mixed_bc_components,
               :mixed_values => mixed_bc_values,
               :traction => x->traction_bc_data(x,1.0))

###################################
# Test Traction BC condition
system_matrix = HDGElasticity.SystemMatrix()
system_rhs = HDGElasticity.SystemRHS()

HDGElasticity.assemble_uniform_boundary_faces!(system_matrix,dgmesh,ufs,
    cellsolvers,facetohelid,2,1,traction_boundary_classifier,bc_data)
HDGElasticity.assemble_traction_face!(system_rhs,dgmesh,ufs,[1.,0.],
    facetohelid,2,1,2)

dofsperelement = ufs.dofsperelement
ndofs = dofsperelement*4

K = dropzeros!(HDGElasticity.sparse(system_matrix,ndofs))
R = HDGElasticity.rhs(system_rhs,ndofs)
sol = K\R

U = reshape(sol,2,:)

testH1 = [0.0 u1  u1  u1 0.0 u1 0.  0.
          0.0 0.0 0.0 u2 u2  u2 0.0 u2]
@test allapprox(U,testH1,1e-10)

############################################
# Test direct global assembly of bc rhs
system_matrix = HDGElasticity.SystemMatrix()
system_rhs = HDGElasticity.SystemRHS()

HDGElasticity.assemble_uniform_boundary_faces!(system_matrix,dgmesh,ufs,
    cellsolvers,facetohelid,traction_boundary_classifier,bc_data)
HDGElasticity.assemble_boundary_faces!(system_rhs,dgmesh,ufs,
    facetohelid,traction_boundary_classifier,bc_data)

dofsperelement = ufs.dofsperelement
ndofs = dofsperelement*4

K = dropzeros!(HDGElasticity.sparse(system_matrix,ndofs))
R = HDGElasticity.rhs(system_rhs,ndofs)
sol = K\R

U = reshape(sol,2,:)

testH1 = [0.0 u1  u1  u1 0.0 u1 0.  0.
          0.0 0.0 0.0 u2 u2  u2 0.0 u2]
@test allapprox(U,testH1,1e-10)

###############################################################################
# Two adjacent elements on boundary

polyorder = 1
numqp = 2
levelset = InterpolatingPolynomial(1,2,polyorder)
mesh = UniformMesh([0.,0.],[2.,1.],[2,1])
interface_location = -5.
levelsetcoeffs = HDGElasticity.levelset_coefficients(
    x->distance_function(x,interface_location),mesh,levelset.basis
)
dgmesh = HDGElasticity.DGMesh(mesh,levelsetcoeffs,levelset)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,polyorder,numqp,
    levelsetcoeffs,levelset)
cellsolvers = HDGElasticity.CellSolvers(dgmesh,ufs,D1,D2,stabilization)
hybrid_element_numbering = HDGElasticity.HybridElementNumbering(dgmesh,ufs)
facetohelid = hybrid_element_numbering.facetohelid
system_matrix = HDGElasticity.SystemMatrix()
system_rhs = HDGElasticity.SystemRHS()

bc_data = Dict(:mixed_components => mixed_bc_components,
               :mixed_values => mixed_bc_values,
               :traction => x->traction_bc_data(x,2.0))

HDGElasticity.assemble_uniform_boundary_faces!(system_matrix,dgmesh,
    ufs,cellsolvers,facetohelid,traction_boundary_classifier,bc_data)
HDGElasticity.assemble_boundary_faces!(system_rhs,dgmesh,ufs,
    facetohelid,traction_boundary_classifier,bc_data)

dofsperelement = ufs.dofsperelement
ndofs = dofsperelement*hybrid_element_numbering.number_of_hybrid_elements
K = dropzeros!(HDGElasticity.sparse(system_matrix,ndofs))
R = HDGElasticity.rhs(system_rhs,ndofs)
sol = K\R

U = reshape(sol,2,:)

testH2 = hcat(testH1,testH1)
testH2[1,9:end] .+= u1
testH2 = testH2[1:2,1:end-2]
@test allapprox(U,testH2,1e-10)

# ###############################################################################
# # Two adjacent elements with coherent interface
polyorder = 1
numqp = 2
levelset = InterpolatingPolynomial(1,2,polyorder)
mesh = UniformMesh([0.,0.],[2.,1.],[2,1])
interface_location = 0.5
levelsetcoeffs = HDGElasticity.levelset_coefficients(
    x->distance_function(x,interface_location),mesh,levelset.basis
)
dgmesh = HDGElasticity.DGMesh(mesh,levelsetcoeffs,levelset)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,polyorder,numqp,
    levelsetcoeffs,levelset)
cellsolvers = HDGElasticity.CellSolvers(dgmesh,ufs,D1,D2,stabilization)
hybrid_element_numbering = HDGElasticity.HybridElementNumbering(dgmesh,ufs)
facetohelid = hybrid_element_numbering.facetohelid
interfacehelid = hybrid_element_numbering.interfacehelid
system_matrix = HDGElasticity.SystemMatrix()
system_rhs = HDGElasticity.SystemRHS()

HDGElasticity.assemble_cut_boundary_faces!(system_matrix,dgmesh,ufs,cellsolvers,
    facetohelid,interfacehelid,traction_boundary_classifier,bc_data)
HDGElasticity.assemble_uniform_boundary_faces!(system_matrix,dgmesh,ufs,
    cellsolvers,facetohelid,traction_boundary_classifier,bc_data)
HDGElasticity.assemble_coherent_interface!(system_matrix,dgmesh,ufs,cellsolvers,
    hybrid_element_numbering)
HDGElasticity.assemble_boundary_faces!(system_rhs,dgmesh,ufs,facetohelid,
    traction_boundary_classifier,bc_data)

dofsperelement = ufs.dofsperelement
ndofs = dofsperelement*hybrid_element_numbering.number_of_hybrid_elements
K = dropzeros!(HDGElasticity.sparse(system_matrix,ndofs))
R = HDGElasticity.rhs(system_rhs,ndofs)
sol = K\R

U = reshape(sol,2,:)

testU1 = [0.0 u1  u1  u1 0.0  u1 0.5u1 0.5u1
          0.0 0.  0.  u2  u2  u2  0.    u2]
testU2 = [0.0 u1  0.0 u1  0.0 0.0 0.5u1 0.5u1
          0.0 0.  u2  u2  0.  u2  0.    u2]
testU3 = [u1 2u1 2u1 2u1 u1 2u1
          0. 0.  0.  u2  u2  u2]
testU = hcat(testU1,testU2,testU3)
@test allapprox(testU,U,1e-10)

###############################################################################
# 9 cell grid
polyorder = 1
numqp = 2
levelset = InterpolatingPolynomial(1,2,polyorder)
mesh = UniformMesh([0.,0.],[1.,1.],[3,3])
interface_location = 1/5
levelsetcoeffs = HDGElasticity.levelset_coefficients(
    x->distance_function(x,interface_location),mesh,levelset.basis
)
dgmesh = HDGElasticity.DGMesh(mesh,levelsetcoeffs,levelset)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,polyorder,numqp,
    levelsetcoeffs,levelset)
cellsolvers = HDGElasticity.CellSolvers(dgmesh,ufs,D1,D2,stabilization)
hybrid_element_numbering = HDGElasticity.HybridElementNumbering(dgmesh,ufs)
facetohelid = hybrid_element_numbering.facetohelid
interfacehelid = hybrid_element_numbering.interfacehelid

system_matrix = HDGElasticity.SystemMatrix()
system_rhs = HDGElasticity.SystemRHS()

bc_data = Dict(:mixed_components => mixed_bc_components,
               :mixed_values => mixed_bc_values,
               :traction => x->traction_bc_data(x,1.0))

HDGElasticity.assemble_uniform_interior_faces!(system_matrix,dgmesh,ufs,
    cellsolvers,facetohelid)

HDGElasticity.assemble_cut_boundary_faces!(system_matrix,dgmesh,ufs,cellsolvers,
    facetohelid,interfacehelid,traction_boundary_classifier,bc_data)
HDGElasticity.assemble_uniform_boundary_faces!(system_matrix,dgmesh,ufs,
    cellsolvers,facetohelid,traction_boundary_classifier,bc_data)
HDGElasticity.assemble_coherent_interface!(system_matrix,dgmesh,ufs,cellsolvers,
    hybrid_element_numbering)
HDGElasticity.assemble_boundary_faces!(system_rhs,dgmesh,ufs,facetohelid,
    traction_boundary_classifier,bc_data)

ndofs = dofsperelement*hybrid_element_numbering.number_of_hybrid_elements
K = dropzeros!(HDGElasticity.sparse(system_matrix,ndofs))
R = HDGElasticity.rhs(system_rhs,ndofs)
sol = K\R

U = reshape(sol,2,:)

hcoords = HDGElasticity.hybrid_element_coordinates(dgmesh,ufs,hybrid_element_numbering)

disp = copy(hcoords)
disp[1,:] .*= u1
disp[2,:] .*= u2
@test allapprox(disp,U,1e-10)


text:
filter(x -> iseven(x), 1:10)
# 5-element Array{Int64,1}:
#   2
#   4
#   6
#   8
#  10
mapreduce(x -> sqrt(x), +, 1:10) #> 22.4682781862041
# equivalent to:
sum(map(x -> sqrt(x), [1:10])) #> 22.4682781862041
#
1:10 |> x -> sqrt(x) |> sum #> 22.4682781862041
1:5 |> exp |> sum #> 233.2041839862982

# clipboard:
a = 42 #> 42
clipboard(a)
workspace()
# a #> ERROR: a not defined
a = clipboard() #> "42"

max_stars_repo_path: Module 1/Chapter10/Old/stlib.jl

max_issues_repo_path: Module 1/Chapter10/Old/stlib.jl

max_forks_repo_path: Module 1/Chapter10/Old/stlib.jl


text:
######### StanVariational Bernoulli example  ###########

using StanVariational

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

stanmodel = VariationalModel("bernoulli", bernoulli_model)

rc = stan_variational(stanmodel; data=bernoulli_data)

if success(rc)

  chns, cnames = read_variational(stanmodel)

end


text:
module MaxPlusAlgebra

export ℝᵋ, ε, ⊕, ⊙

struct ℝᵋ <: Real
    value::Float64
end

ℝᵋ(x::ℝᵋ) = ℝᵋ(x.value)

Base.max(x::ℝᵋ, y::ℝᵋ) = ℝᵋ(max(x.value, y.value))
Base.:+(x::ℝᵋ, y::ℝᵋ) = ℝᵋ(x.value + y.value)
Base.typemin(::Type{ℝᵋ}) = ℝᵋ(typemin(Float64))

const global ε = typemin(ℝᵋ)

Base.promote_rule(::Type{ℝᵋ}, ::Type{T}) where {T<:Number} = ℝᵋ

⊙(x::T, y::T) where {T<:Real} = x + y
⊙(x::Number, y::Number) = ⊙(promote(x,y)...)

⊕(x::T, y::T) where {T<:Real} = max(x,y)
⊕(x::Number, y::Number) = ⊕(promote(x,y)...)

ℝᵋ(A::Array) = map(ℝᵋ, A)

function ⊡(A::Matrix{T}, B::Matrix{T}, ⊞ = +, ⊠ = *, ∅ = 0) where T<:Number
    @assert size(A)[end] == size(B)[1]
    S = Matrix{T}(undef, size(A)[1], size(B)[end])
    for i = 1:size(A)[1], j = 1:size(B)[end]
        S[i,j] = ∅
        for k = 1:size(A)[end]
            S[i,j] = ⊞(S[i,j], ⊠(A[i,k], B[k,j]))
        end
    end
    return S
end

⊙(A::Matrix{T}, B::Matrix{T}) where T<:Number = ⊡(A, B, ⊕, ⊙, typemin(T))

end


text:
# mdp step simulator and stepthrough
@testset "gridstepsim" begin
    p = LegacyGridWorld()
    solver = RandomSolver(MersenneTwister(2))
    policy = solve(solver, p)
    sim = StepSimulator("s,sp,r,a,ai", rng=MersenneTwister(3), max_steps=100)
    n_steps = 0
    for (s, sp, r, a, ai) in simulate(sim, p, policy)
        @test isa(s, statetype(p))
        @test isa(sp, statetype(p))
        @test isa(r, Float64)
        @test isa(a, actiontype(p))
        @test isa(ai, Nothing)
        n_steps += 1
    end
    @test n_steps <= 100

    n_steps = 0
    for s in stepthrough(p, policy, "s", rng=MersenneTwister(4), max_steps=100)
        @test isa(s, statetype(p))
        n_steps += 1
    end
    @test n_steps <= 100
end

# pomdp step simulator and stepthrough
@testset "babystepsim" begin
    p = BabyPOMDP()
    policy = FeedWhenCrying()
    up = PreviousObservationUpdater()
    sim = StepSimulator("s,sp,r,a,b,ui,i,ai", rng=MersenneTwister(3), max_steps=100)
    s_init = rand(sim.rng, [true, false])
    b_init = false
    n_steps = 0
    for (s, sp, r, a, b, ui, i, ai) in simulate(sim, p, policy, up, b_init, s_init)
        @test isa(s, statetype(p))
        @test isa(sp, statetype(p))
        @test isa(r, Float64)
        @test isa(a, actiontype(p))
        @test isa(b, Bool)
        @test isa(ui, Nothing)
        @test isa(ai, Nothing)
        n_steps += 1
    end
    @test n_steps == 100

    # test with only two args
    collect(stepthrough(p, Starve(); max_steps=3))
end
@testset "stepthroughfeed" begin
    p = BabyPOMDP()
    policy = FeedWhenCrying()
    up = PreviousObservationUpdater()
    s_init = rand(MersenneTwister(3), [true, false])
    b_init = false
    n_steps = 0
    for r in stepthrough(p, policy, up, b_init, s_init, "r", rng=MersenneTwister(4), max_steps=100)
        @test isa(r, Float64)
        @test r <= 0
        n_steps += 1
    end
    @test n_steps == 100
end

# example from stepthrough documentation
@testset "stepthroughrand" begin
    p = BabyPOMDP()
    policy = RandomPolicy(p)

    for (s, a, o, r) in stepthrough(p, policy, "s,a,o,r", max_steps=10)
        println("in state $s")
        println("took action $o")
        println("received observation $o and reward $r")
    end
end

struct SymPOMDP <: POMDP{Symbol, Symbol, Symbol} end
@testset "stepthrougherr" begin
    m = SymPOMDP()
    @test_throws ErrorException stepthrough(m, RandomPolicy(m), NothingUpdater(), [:init], :init)
end


text:
module UncValue

import Base: +, -, *, /, \, ÷, ^, ==, !=, <, >, <=, >=
import Printf: @printf

struct Value{T <: Real}
    x::T
    σ::T

    Value{T}(x, σ) where {T <: Real} = σ < 0 ? error("Negative uncertainty has no meaning.") : new{T}(x, σ)
end
function Value(x::T, σ::U) where {T <: Real, U <: Real}
    xT, σT = promote(x, σ)
    @warn "Value was initialize with different types, promoting $(typeof(x)) & $(typeof(σ)) to $(typeof(xT))."
    Value{typeof(xT)}(xT, σT)
end
Value(x::T, σ::T) where {T <: Real} = Value{T}(x, σ)
Value(x::T) where {T <: Real} = Value{T}(x, zero(x))

export Value, val, unc, set_unc

## MATHEMATICAL OPERATORS

+(a::Value, b::T) where {T <: Real} = Value(a.x + b, a.σ)
+(a::T, b::Value) where {T <: Real} = Value(a + b.x, b.σ)
+(a::Value, b::Value) = Value(a.x + b.x, hypot(a.σ, b.σ))

-(a::Value) = Value(-a.x, a.σ)
-(a::Value, b::T) where {T <: Real} = Value(a.x - b, a.σ)
-(a::T, b::Value) where {T <: Real} = Value(a - b.x, b.σ)
-(a::Value, b::Value) = Value(a.x - b.x, hypot(a.σ, b.σ))

*(a::Value, b::T) where {T <: Real} = Value(a.x * b, a.σ * abs(b))
*(a::T, b::Value) where {T <: Real} = Value(a * b.x, b.σ * abs(a))
*(a::Value, b::Value) = Value(a.x * b.x, hypot(a.σ * b.x, a.x * b.σ))

/(a::Value, b::T) where {T <: Real} = Value(a.x / b, a.σ / abs(b))
/(a::T, b::Value) where {T <: Real} = Value(a / b.x, abs(a) * b.σ / b.x^2)
/(a::Value, b::Value) = Value(a.x / b.x, hypot(a.σ / b.x, b.σ * a.x / b.x^2))

\(a::Value, b::T) where {T <: Real} = Value(a.x \ b, a.x^2 \ b * a.σ)
\(a::T, b::Value) where {T <: Real} = Value(a \ b.x, a \ b.σ)
\(a::Value, b::Value) = Value(a.x \ b.x, hypot( a.x^2 \ b.x * a.σ, a.x \ b.σ))

÷(a::Value, b::T) where {T <: Real} = Value(a.x ÷ b, a.σ ÷ b)
÷(a::T, b::Value) where {T <: Real} = Value(a ÷ b.x, b.σ * a ÷ b.x^2)
÷(a::Value, b::Value) = Value(a.x ÷ b.x, round(typeof(a.x), hypot(a.σ ÷ b.x, b.σ * a.x ÷ b.x^2)))

^(a::Value, b::T) where {T <: Real} = Value(a.x ^ b, abs(b * a.x^(b-1)) * a.σ)
^(a::T, b::Value) where {T <: Real} = Value(a ^ b.x, abs(a^b.x * log(abs(a))) * b.σ)
^(a::Value, b::Value) where {T <: Real} =
    Value(a.x ^ b.x, hypot(b.x * a.x^(b.x-1) * a.σ, a.x^b.x * log(abs(a.x)) * b.σ))

==(a::Value, b::T) where {T <: Real} = a.x == b
==(a::T, b::Value) where {T <: Real} = a == b.x
==(a::Value, b::Value) = a.x == b.x

!=(a::Value, b::T) where {T <: Real} = a.x != b
!=(a::T, b::Value) where {T <: Real} = a != b.x
!=(a::Value, b::Value) = a.x != b.x

>(a::Value, b::T) where {T <: Real} = a.x > b
>(a::T, b::Value) where {T <: Real} = a > b.x
>(a::Value, b::Value) = a.x > b.x

<(a::Value, b::T) where {T <: Real} = a.x < b
<(a::T, b::Value) where {T <: Real} = a < b.x
<(a::Value, b::Value) = a.x < b.x

>=(a::Value, b::T) where {T <: Real} = a.x >= b
>=(a::T, b::Value) where {T <: Real} = a >= b.x
>=(a::Value, b::Value) = a.x >= b.x

<=(a::Value, b::T) where {T <: Real} = a.x <= b
<=(a::T, b::Value) where {T <: Real} = a <= b.x
<=(a::Value, b::Value) = a.x <= b.x

## MATHEMATICAL FUNCTIONS
Base.abs(a::Value) = Value(abs(a.x), a.σ)
Base.abs2(a::Value) = Value(abs2(a.x), 2 * abs(a.x) * a.σ)

Base.sin(a::Value) = Value(sin(a.x), abs(cos(a.x)) * a.σ)
Base.sind(a::Value) = Value(sind(a.x), abs(cosd(a.x)) * a.σ)
Base.sinh(a::Value) = Value(sinh(a.x), cosh(a.x) * a.σ)
Base.Math.sinpi(a::Value) = Value(sinpi(a.x), abs(π * cospi(a.x)) * a.σ)

Base.cos(a::Value) = Value(cos(a.x), abs(sin(a.x)) * a.σ)
Base.cosd(a::Value) = Value(cosd(a.x), abs(sind(a.x)) * a.σ)
Base.cosh(a::Value) = Value(cosh(a.x), abs(sinh(a.x)) * a.σ)
Base.Math.cospi(a::Value) = Value(cospi(a.x), abs(π * sinpi(a.x)) * a.σ)

Base.tan(a::Value) = Value(tan(a.x),  a.σ * sec(a.x)^2)
Base.tand(a::Value) = Value(tand(a.x),  a.σ * secd(a.x)^2)
Base.tanh(a::Value) = Value(tanh(a.x),  a.σ * sech(a.x)^2)

Base.asin(a::Value) = Value(asin(a.x),  a.σ / sqrt(1-a.x^2))
Base.Math.asind(a::Value) = Value(asind(a.x),  a.σ / sqrt(1-a.x^2))
Base.asinh(a::Value) = Value(asinh(a.x),  a.σ / sqrt(1 + a.x^2))

Base.acos(a::Value) = Value(acos(a.x),  a.σ / sqrt(1-a.x^2))
Base.Math.acosd(a::Value) = Value(acosd(a.x),  a.σ / sqrt(1-a.x^2))
Base.acosh(a::Value) = Value(acosh(a.x), a.σ / sqrt(a.x^2-1))

Base.atan(a::Value) = Value(atan(a.x), a.σ / (1 + a.x^2))
Base.atan(a::Value, b::T) where {T <: Real} = Value(atan(a.x, b), abs(b) * a.σ / (a.x^2 + b^2))
Base.atan(a::T, b::Value) where {T <: Real} = Value(atan(a, b.x), abs(a) * b.σ / (a^2 + b.x^2))
Base.atan(a::Value, b::Value) =
    Value(atan(a.x, b.x), hypot(b.x * a.σ, a.x * b.σ) / (a.x^2 + b.x^2))
Base.Math.atand(a::Value) = Value(atand(a.x),  a.σ / (1 + a.x^2))
Base.Math.atand(a::Value, b::T) where {T <: Real} = Value(atand(a.x, b), abs(b) * a.σ / (a.x^2 + b^2))
Base.Math.atand(a::T, b::Value) where {T <: Real} = Value(atand(a, b.x), abs(a) * b.σ / (a^2 + b.x^2))
Base.Math.atand(a::Value, b::Value) =
    Value(atand(a.x, b.x), hypot(b.x * a.σ, a.x * b.σ) / (a.x^2 + b.x^2))
Base.atanh(a::Value) = Value(atanh(a.x),  a.σ / (1 - a.x^2))

Base.Math.sec(a::Value) = Value(sec(a.x), abs(sec(a.x) * tan(a.x)) * a.σ)
Base.Math.secd(a::Value) = Value(secd(a.x), abs(secd(a.x) * tand(a.x)) * a.σ)
Base.Math.sech(a::Value) = Value(sech(a.x), abs(tanh(a.x) * sech(a.x)) * a.σ)

Base.Math.csc(a::Value) = Value(csc(a.x), abs(cot(a.x) * csc(a.x)) * a.σ)
Base.Math.cscd(a::Value) = Value(cscd(a.x), abs(cotd(a.x) * cscd(a.x)) * a.σ)
Base.Math.csch(a::Value) = Value(csch(a.x), abs(coth(a.x) * csch(a.x)) * a.σ)

Base.Math.cot(a::Value) = Value(cot(a.x), csc(a.x)^2 * a.σ)
Base.Math.cotd(a::Value) = Value(cotd(a.x), cscd(a.x)^2 * a.σ)
Base.Math.coth(a::Value) = Value(coth(a.x), csch(a.x)^2 * a.σ)

Base.Math.asec(a::Value) = Value(asec(a.x), a.σ / sqrt(a.x^2 - 1))
Base.Math.asecd(a::Value) = Value(asecd(a.x), a.σ / sqrt(a.x^2 - 1))
Base.Math.asech(a::Value) = Value(asech(a.x), a.σ / sqrt(1 - a.x^2))

Base.Math.acsc(a::Value) = Value(acsc(a.x), a.σ / sqrt(a.x^2 - 1))
Base.Math.acscd(a::Value) = Value(acscd(a.x), a.σ / sqrt(a.x^2 - 1))
Base.Math.acsch(a::Value) = Value(acsch(a.x), a.σ / sqrt(a.x^2 + 1))

Base.Math.acot(a::Value) = Value(acot(a.x), a.σ / (a.x^2 + 1))
Base.Math.acotd(a::Value) = Value(acotd(a.x), a.σ / (a.x^2 + 1))
Base.Math.acoth(a::Value) = Value(acoth(a.x), a.σ / abs(1 - a.x^2))

Base.Math.sinc(a::Value) = Value(sinc(a.x), a.σ * abs(cosc(a.x)))
Base.Math.cosc(a::Value) = Value(cosc(a.x), a.σ * abs((2 * a.x * cos(a.x) + (a.x^2 - 2) * sin(a.x)) / a.x^3))

Base.Math.deg2rad(a::Value) = Value(deg2rad(a.x), deg2rad(a.σ))
Base.Math.rad2deg(a::Value) = Value(rad2deg(a.x), rad2deg(a.σ))

Base.exp(a::Value) = Value(exp(a.x), exp(a.x) * a.σ)
Base.exp2(a::Value) = Value(exp2(a.x), exp2(a.x) * a.σ * log(2))
Base.exp10(a::Value) = Value(exp10(a.x), exp10(a.x) * a.σ * log(10))

Base.log(a::Value) = Value(log(a.x), a.σ / abs(a.x))
Base.log(b::T, a::Value) where {T <: Real} = Value(log(b, a.x), a.σ / abs(a.x * log(b)))
Base.log2(a::Value) = Value(log2(a.x), a.σ / abs(a.x * log(2)))
Base.log10(a::Value) = Value(log10(a.x), a.σ / abs(a.x * log(10)))
Base.log1p(a::Value) = Value(log1p(a.x), a.σ / abs(1 + a.x))

Base.sqrt(a::Value) = Value(sqrt(a.x), a.σ / (2 * sqrt(a.x)))
Base.hypot(a::Value, b::T) where {T <: Real} = Value(hypot(a.x, b), abs(a.x) * a.σ / hypot(a.x, b))
Base.hypot(a::T, b::Value) where {T <: Real} = Value(hypot(a, b.x), abs(b.x) * b.σ / hypot(a, b.x))
function Base.hypot(a::Value, b::Value)
    h = hypot(a.x, b.x)
    Value(h, hypot(a.x * a.σ, b.x * b.σ) / h)
end

Base.cbrt(a::Value) = Value(cbrt(a.x), a.σ / (3 * cbrt(a.x)^2))

## OTHER FUNCTIONS

Base.sign(a::Value) = sign(a.x)
Base.signbit(a::Value) = signbit(a.x)

Base.inv(a::Value) = Value(inv(a.x), a.σ * inv(a.x)^2)

Base.fld(a::Value, b::T) where {T <: Real} = Value(fld(a.x, b), fld(a.σ, b))
Base.fld(a::T, b::Value) where {T <: Real} = Value(fld(a, b.x), fld(b.σ * a, b.x^2))
Base.fld(a::Value, b::Value) = Value(fld(a.x, b.x), floor(typeof(a.x), hypot(fld(a.σ, b.x), fld(b.σ * a.x, b.x^2))))

Base.cld(a::Value, b::T) where {T <: Real} = Value(cld(a.x, b), cld(a.σ, b))
Base.cld(a::T, b::Value) where {T <: Real} = Value(cld(a, b.x), cld(b.σ * a, b.x^2))
Base.cld(a::Value, b::Value) = Value(cld(a.x, b.x), ceil(typeof(a.x), hypot(cld(a.σ, b.x), cld(b.σ * a.x, b.x^2))))

Base.cmp(a::Value, b::T) where {T <: Real} = cmp(a.x, b)
Base.cmp(a::T, b::Value) where {T <: Real} = cmp(a, b.x)
Base.cmp(a::Value, b::Value) = cmp(a.x, b.x)

Base.isless(a::Value, b::T) where {T <: Real} = cmp(a.x, b) < 0
Base.isless(a::T, b::Value) where {T <: Real} = cmp(a, b.x) < 0
Base.isless(a::Value, b::Value) = cmp(a.x, b.x) < 0

function Base.clamp(x::T, lh::Value) where {T <: Real}
    if x > lh.x + lh.σ
        return lh.x + lh.σ
    elseif x < lh.x - lh.σ
        return lh.x - lh.σ
    end

    x
end

function Base.clamp(rh::Value, lh::Value) where {T <: Real}
    if rh.x > lh.x + lh.σ
        return Value(lh.x + lh.σ, rh.σ)
    elseif rh.x < lh.x - lh.σ
        return Value(lh.x - lh.σ, rh.σ)
    end

    rh
end

#check how to return also the uncertainty!
Base.max(a::Value) = a.x + a.σ
Base.max(a::Value, b::T) where {T <: Real} = cmp(a.x, b) <= 0 ? b : a
Base.max(a::T, b::Value) where {T <: Real} = cmp(a, b.x) <= 0 ? b : a
Base.max(a::Value, b::Value)= cmp(a.x, b.x) <= 0 ? b : a

Base.min(a::Value) = a.x - a.σ
Base.min(a::Value, b::T) where {T <: Real} = cmp(a.x, b) <= 0 ? a : b
Base.min(a::T, b::Value) where {T <: Real} = cmp(a, b.x) <= 0 ? a : b
Base.min(a::Value, b::Value)= cmp(a.x, b.x) <= 0 ? a : b

Base.isapprox(a::Value, b::T; significance::Int = 1) where {T <: Real} =
    a.x - significance * a.σ <= b <= a.x + significance * a.σ
Base.isapprox(a::T, b::Value; significance::Int = 1) where {T <: Real} =
    b.x - significance * b.σ <= a <= b.x + significance * b.σ
Base.isapprox(a::Value, b::Value; significance::Int = 1) =
    (b.x - significance * b.σ <= a.x <= b.x + significance * b.σ) |
    (a.x - significance * a.σ <= b.x <= a.x + significance * a.σ)

Base.length(a::Value) = length(a.x)

Base.iterate(a::Value) = (a, nothing)
Base.iterate(a::Value, b::Nothing) = (a, 2)

## PRINTED RESULTS# https://docs.julialang.org/en/v1/base/math/
# precision(a::T) where {T <: Real} = floor(Int, log10(a)) - 1
Base.precision(a::Value) = floor(Int, log10(a.σ)) - 1

function Base.show(io::IO, a::Value)
    p = precision(a)
    fix_x = round(a.x * 10.0^(-p); digits = 2, base = 10) / 10
    fix_ux = round(a.σ * 10.0^(-p); digits = 2, base = 10) / 10
    p += 1
    if p == 0
        @printf(io, "%.1f ± %.1f\n", fix_x, fix_ux)
    else
        @printf(io, "(%.1f ± %.1f)·10^%d\n", fix_x, fix_ux, p)
    end
end

## NEW FUNCTIONS
val(item::Number) = item
val(item::Value) = item.x
function val(A::Array)
    type = typeof(A[1])
    if type <: Value
        type = typeof(A[1].x)
    end

    B = Array{type, ndims(A)}(undef, size(A))
    for index in eachindex(A)
        B[index] = val(A[index])
    end

    B
end

unc(item::Number) = zero(item)
unc(item::Value) = item.σ
function unc(A::Array)
    type = typeof(A[1])
    if type <: Value
        type = typeof(A[1].x)
    end

    B = Array{type, ndims(A)}(undef, size(A))
    for index in eachindex(A)
        B[index] = unc(A[index])
    end

    B
end

function set_unc(A::Array{T, N}, unc::Real) where {T, N}
    type = typeof(A[1])
    if type <: Value
        type = typeof(A[1].x)
    end

    B = Array{Value{type}, ndims(A)}(undef, size(A))
    for index in eachindex(A)
        B[index] = Value(val(A[index]), unc)
    end

    B
end

function set_unc(A::Array{T, N}, unc::Array{T2, N}) where {T, T2 <: Real, N}
    type = typeof(A[1])
    if type <: Value
        type = typeof(A[1].x)
    end

    B = Array{Value{type}, ndims(A)}(undef, size(A))
    for index in eachindex(A)
        B[index] = Value(val(A[index]), unc[index])
    end

    B
end

end


text:
# Untested work in progres...

packet_count = 1000

data_producers = [

    :large_packet_stream => http -> begin
        for i in 1:packet_count
            send_large_chunk(http)
        end
    end,

    :small_packet_stream => http -> begin
        for i in 1:packet_count
            send_small_chunk(http)
        end
    end,

    :large_packets_with_delays => http -> begin
        for i in 1:packet_count
            sleep(0.1)
            send_large_chunk(http)
        end
    end,

    :small_packets_with_delays => http -> begin
        for i in 1:packet_count
            sleep(0.1)
            send_large_chunk(http)
        end
    end,

    :random_packets_and_delays => http -> begin
        for i in 1:packet_count
            sleep(rand(1:300)/1000)
            send_chunk(http, rand(10:10000)
        end
    end
]

idle_behaviours = [

    :delay_1_minute => http -> sleep(60),

    :no_delay => nothing
]

close_behaviours = [

    :tls_close_notify => http -> begin
        closewrite(http)
        close(http.stream.c::SSLContext)
    end,

    :tcp_close => http -> begin
        closewrite(http)
        close(http.stream.c.bio::TCPSocket)
    end
]

combinations = []
for p1 in data_producers,
    p2 in data_producers,
    i in idle_behaviours,
    c in close_behaviours
    push!(combinations, [p1..., p2..., i..., c...])
end
using Random
combinations = Iterators.Stateful(Iterators.cycle(shuffle(combinations)))

timestamp() = string(round(Int, time() * 1000))

HTTP.listen("127.0.0.1", 8080) do http::HTTP.Stream

    p1n, produce1f,
    p2n, produce2f,
    in, idlef,
    cn, closef = popfirst!(combinations)

    name = "$p1n, $p2n, $in and $cn")
    println(name)

    HTTP.setheader(http, "x-test-cobination" => name)
    HTTP.setheader(http, "x-test-tzero" => timestamp())

    produce1f(http)
    produce2f(http)
    idlef(http)
    closef(http)
end

send_small_chunk(http, 10)
send_large_chunk(http, 1000)

send_chunk(http, size)
    chunk = IOBuffer()
    write(chunk, timestamp(), "\n")
    while size > 0
        write(chunk, "0123456789")
        size -= 10
    end
    write(http, chunk)
end


text:
# Generate low-level api from METADATA

# Types defined by the api
immutable NvimApiObject{N} <: NvimObject
    client::NvimClient
    # TODO: use a fixarray or UInt64
    hnd::Vector{UInt8}
end

import Base.==
=={N}(a::NvimApiObject{N},b::NvimApiObject{N}) = a.hnd == b.hnd

NvimApiObject(c, e::Ext) = NvimApiObject{Int(e.typecode)}(c, e.data)
#Not really module-interface clean, I know...
MsgPack.pack{N}(s, o::NvimApiObject{N}) = MsgPack.pack(s, Ext(N, o.hnd))

symbolize(val::Dict) = Dict{Symbol,Any}([(symbolize(k),symbolize(v)) for (k,v) in val])
symbolize(val::Vector{UInt8}) = symbol(bytestring(val))
symbolize(val::ByteString) = symbol(val)
symbolize(val::Vector) = [symbolize(v) for v in val]
symbolize(val) = val

function _get_metadata()
    data = readall(`nvim --api-info`)
    return symbolize(unpack(data))
end

const _metadata = _get_metadata()
const _types = _metadata[:types]
const _functions = _metadata[:functions]

# will break if the api starts using overloading
const api_methods = [f[:name] => f for f in _functions]

typealias Bytes Union{ByteString, Vector{UInt8}}
const typemap = Dict{Symbol,Type}(
    :Integer => Integer,
    :Boolean => Bool,
    :String => Bytes,
    :Array => Vector{Any},
    :Dictionary => Dict,
    :Object => Any,
    :void => Void
)

for (name, info) in _types
    id = info[:id]
    @eval begin
        typealias $(name) NvimApiObject{$id}
        typemap[$(Meta.quot(name))] = $name
    end
end

function parsetype(t::Symbol)
    components = split(rstrip(string(t), ')'), '(')
    top = symbol(components[1])
    params = length(components) > 1 ? split(components[2], ", ") : []
    (top, params)
end

function gettype(t::Symbol)
    top, params = parsetype(t)
    if haskey(typemap, top)
        typemap[top]
    elseif top == :ArrayOf
        tcontains = gettype(symbol(params[1]))
        if length(params) == 1
            Vector{tcontains}
        else
            Tuple{fill(tcontains, parse(Int, params[2]))...}
        end
    else
        Any
    end
end

# fake covariant Vectors
checkarg{T}(::Type{Vector{T}}, val::Vector{T}) = val
checkarg{T}(::Type{Vector{T}}, val::Vector) = convert(Vector{T}, val)
# I'm pretty sure julia somehow can express subdiagnal dispatch,
# but implement the cases manually for now
checkarg{T<:Bytes}(::Type{Vector{Bytes}}, val::Vector{T}) = val
checkarg{T<:Integer}(::Type{Vector{Integer}}, val::Vector{T}) = val
# this must be specialcased though
checkarg(::Type{Vector{Integer}}, val::Vector{UInt8}) = Int[val...]
# tuples are covariant, so no problem
checkarg{T}(::Type{Vector{T}}, val::(Tuple{Vararg{T}})) = [val...]

retconvert(typ::Union{Type{Any},Type{Bytes}}, c, val::Union{ByteString, Vector{UInt8}}) = bytestring(val)
# needed for disambiguation
retconvert(typ::Type{Any}, c, val::Vector{UInt8}) = bytestring(val)
retconvert(typ::Union{Type{Any}, Type{Bool}}, c, val::Bool) = val
# this assumes the current unpack implementation in MsgPack,
# where all int types get promoted to 64 bit
retconvert(typ::Union{Type{Any},Type{Integer}}, c, val::Int64) = val
retconvert(typ::Union{Type{Any},Type{Void}}, c, val::Void) = nothing
retconvert(typ::Union{Type{Any},Type{Dict}}, c, val::Dict) = Dict{ByteString,Any}([(bytestring(k),retconvert(Any,c,v)) for (k,v) in val])
# we assume Msgpack only generates untyped arrays
retconvert{T}(typ::Type{Vector{T}}, c, val::Vector) = [retconvert(T,c,v) for v in val]
retconvert(typ::Type{Any}, c, val::Vector) = Any[retconvert(Any,c,v) for v in val]

# this is not really strict enough
#retconvert{T<:Tuple}(typ::Type{T}, c, val::Vector) = tuple(val...)::T
# handle the only case, 2-tuples, manually for now
retconvert{T,U}(::Type{Tuple{T,U}}, c, val::Vector) = (retconvert(T,c,val[1]),retconvert(U,c,val[2]))

retconvert{N}(typ::Type{NvimApiObject{N}}, c, val::Ext) = NvimApiObject(c, val)::NvimApiObject{N}
#retconvert{T}(typ::Type{T}, c, val::T) = val

# a stagedfunction will probably be simpler and better
function build_function(f)
    name = f[:name]
    params = f[:parameters]
    tret = f[:return_type]

    parts = split(string(name), "_"; limit=2)
    reciever = parts[1]
    shortname = symbol(parts[2])
    if shortname == :eval; shortname = :vim_eval; end

    body = Expr[]
    args = Symbol[symbol(string("a_",p[2])) for p in params]
    j_args = Expr[]

    for (i,p) in enumerate(params)
        #this is probably too restrictive sometimes,
        # use convert for some types (sequences)?
        t = gettype(p[1])
        # if type is an array
        # we will allow any Vector argument
        # and dynamically check if not a subtype
        erased = t
        if t <: Vector
            erased = Union{Vector, Tuple{Vararg{eltype(t)}}}
            push!(body, :( $(args[i]) = checkarg($t, $(args[i]))))
        end
        push!(j_args, :( $(args[i])::($erased) ))
    end

    if reciever == "vim"
        unshift!(j_args, :( c::NvimClient))
    else
        unshift!(body, :( c = ($(args[1])).client ) )
    end

    #when array constructor non-concatenating, we could drop the Any
    arglist = :( Any[] )
    append!(arglist.args, args)
    push!(body, :( res = send_request(c, $(Meta.quot(name)), $arglist)))

    # "representative" type, not neccessarily the one
    # retconvert will assert
    typ = gettype(tret)
    push!(body, :( retconvert($typ, c, res) ))

    j_call = Expr(:call, shortname, j_args...)
    fun = Expr(:function, j_call, Expr(:block, body...) )
    #println(fun)
    fun
end

#TODO: maybe this should be in a submodule
# then one could do `importall Neovim.API`
for f in _functions
    eval(build_function(f))
end


text:
## Triangular
immutable Triangular{T,S<:AbstractMatrix{T},UpLo,IsUnit} <: AbstractMatrix{T}
    data::S
end
function Triangular{T}(A::AbstractMatrix{T}, uplo::Symbol, isunit::Bool=false)
    uplo != :L && uplo != :U && throw(ArgumentError("uplo argument must be either :U or :L"))
    return Triangular{T,typeof(A),uplo,isunit}(A)
end

######################
# BlasFloat routines #
######################

### Note! the BlasFloat restriction can be removed if generic triangular multiplication methods are written.
for (func1, func2) in ((:*, :A_mul_B!), (:Ac_mul_B, :Ac_mul_B!), (:/, :A_rdiv_B!))
    @eval begin
        ($func1){T<:BlasFloat,S<:StridedMatrix,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, B::Triangular{T,S,UpLo,IsUnit}) = ($func2)(A, full(B))
        ($func1){T<:BlasFloat,S<:StridedMatrix,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, B::StridedVecOrMat) = ($func2)(A, copy(B))
    end
end
for (func1, func2) in ((:A_mul_Bc, :A_mul_Bc!), (:A_rdiv_Bc, :A_rdiv_Bc!))
    @eval begin
        ($func1){T<:BlasFloat,S<:StridedMatrix,UpLo,IsUnit}(A::StridedMatrix, B::Triangular{T,S,UpLo,IsUnit}) = ($func2)(copy(A), B)
    end
end

# Vector multiplication
A_mul_B!{T<:BlasFloat,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, b::StridedVector{T}) = BLAS.trmv!(UpLo == :L ? 'L' : 'U', 'N', IsUnit ? 'U' : 'N', A.data, b)
Ac_mul_B!{T<:BlasReal,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, b::StridedVector{T}) = BLAS.trmv!(UpLo == :L ? 'L' : 'U', 'T', IsUnit ? 'U' : 'N', A.data, b)
Ac_mul_B!{T<:BlasComplex,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, b::StridedVector{T}) = BLAS.trmv!(UpLo == :L ? 'L' : 'U', 'C', IsUnit ? 'U' : 'N', A.data, b)
At_mul_B!{T<:BlasFloat,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, b::StridedVector{T}) = BLAS.trmv!(UpLo == :L ? 'L' : 'U', 'T', IsUnit ? 'U' : 'N', A.data, b)

# Matrix multiplication
A_mul_B!{T<:BlasFloat,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, B::StridedMatrix{T}) = BLAS.trmm!('L', UpLo == :L ? 'L' : 'U', 'N', IsUnit ? 'U' : 'N', one(T), A.data, B)
A_mul_B!{T<:BlasFloat,S,UpLo,IsUnit}(A::StridedMatrix{T}, B::Triangular{T,S,UpLo,IsUnit}) = BLAS.trmm!('R', UpLo == :L ? 'L' : 'U', 'N', IsUnit ? 'U' : 'N', one(T), B.data, A)
Ac_mul_B!{T<:BlasComplex,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, B::StridedMatrix{T}) = BLAS.trmm!('L', UpLo == :L ? 'L' : 'U', 'C', IsUnit ? 'U' : 'N', one(T), A.data, B)
Ac_mul_B!{T<:BlasReal,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, B::StridedMatrix{T}) = BLAS.trmm!('L', UpLo == :L ? 'L' : 'U', 'T', IsUnit ? 'U' : 'N', one(T), A.data, B)
A_mul_Bc!{T<:BlasComplex,S,UpLo,IsUnit}(A::StridedMatrix{T}, B::Triangular{T,S,UpLo,IsUnit}) = BLAS.trmm!('R', UpLo == :L ? 'L' : 'U', 'C', IsUnit ? 'U' : 'N', one(T), B.data, A)
A_mul_Bc!{T<:BlasReal,S,UpLo,IsUnit}(A::StridedMatrix{T}, B::Triangular{T,S,UpLo,IsUnit}) = BLAS.trmm!('R', UpLo == :L ? 'L' : 'U', 'T', IsUnit ? 'U' : 'N', one(T), B.data, A)

A_ldiv_B!{T<:BlasFloat,S<:AbstractMatrix,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, B::StridedVecOrMat{T}) = LAPACK.trtrs!(UpLo == :L ? 'L' : 'U', 'N', IsUnit ? 'U' : 'N', A.data, B)
function \{T<:BlasFloat,S<:AbstractMatrix,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, B::StridedVecOrMat{T})
    x = A_ldiv_B!(A, copy(B))
    errors = LAPACK.trrfs!(UpLo == :L ? 'L' : 'U', 'N', IsUnit ? 'U' : 'N', A.data, B, x)
    all(isfinite, [errors...]) || all([errors...] .< one(T)/eps(T)) || warn("""Unreasonably large error in computed solution:
forward errors:
$(errors[1])
backward errors:
$(errors[2])""")
    x
end
Ac_ldiv_B!{T<:BlasReal,S<:AbstractMatrix,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, B::StridedVecOrMat{T}) = LAPACK.trtrs!(UpLo == :L ? 'L' : 'U', 'T', IsUnit ? 'U' : 'N', A.data, B)
Ac_ldiv_B!{T<:BlasComplex,S<:AbstractMatrix,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, B::StridedVecOrMat{T}) = LAPACK.trtrs!(UpLo == :L ? 'L' : 'U', 'C', IsUnit ? 'U' : 'N', A.data, B)

A_rdiv_B!{T<:BlasFloat,S<:AbstractMatrix,UpLo,IsUnit}(A::StridedVecOrMat{T}, B::Triangular{T,S,UpLo,IsUnit}) = BLAS.trsm!('R', UpLo == :L ? 'L' : 'U', 'N', IsUnit ? 'U' : 'N', one(T), B.data, A)
A_rdiv_Bc!{T<:BlasReal,S<:AbstractMatrix,UpLo,IsUnit}(A::StridedMatrix{T}, B::Triangular{T,S,UpLo,IsUnit}) = BLAS.trsm!('R', UpLo == :L ? 'L' : 'U', 'T', IsUnit ? 'U' : 'N', one(T), B.data, A)
A_rdiv_Bc!{T<:BlasComplex,S<:AbstractMatrix,UpLo,IsUnit}(A::StridedMatrix{T}, B::Triangular{T,S,UpLo,IsUnit}) = BLAS.trsm!('R', UpLo == :L ? 'L' : 'U', 'C', IsUnit ? 'U' : 'N', one(T), B.data, A)

inv{T<:BlasFloat,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T,S,UpLo,IsUnit}(LAPACK.trtri!(UpLo == :L ? 'L' : 'U', IsUnit ? 'U' : 'N', copy(A.data)))

#Eigensystems
function eigvecs{T<:BlasFloat,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit})
    if UpLo == :U
        V = LAPACK.trevc!('R', 'A', Array(Bool,1), A.data)
    else # Uplo == :L
        V = LAPACK.trevc!('L', 'A', Array(Bool,1), A.data')
    end
    for i=1:size(V,2) #Normalize
        V[:,i] /= norm(V[:,i])
    end
    V
end

function cond{T<:BlasFloat,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, p::Real=2)
    chksquare(A)
    if p==1
        return inv(LAPACK.trcon!('O', UpLo == :L ? 'L' : 'U', IsUnit ? 'U' : 'N', A.data))
    elseif p==Inf
        return inv(LAPACK.trcon!('I', UpLo == :L ? 'L' : 'U', IsUnit ? 'U' : 'N', A.data))
    else #use fallback
        return cond(full(A), p)
    end
end

####################
# Generic routines #
####################

size(A::Triangular, args...) = size(A.data, args...)

convert{T,S,UpLo,IsUnit}(::Type{Triangular{T,S,UpLo,IsUnit}}, A::Triangular{T,S,UpLo,IsUnit}) = A
convert{T,S,UpLo,IsUnit}(::Type{Triangular{T,S,UpLo,IsUnit}}, A::Triangular) = Triangular{T,S,UpLo,IsUnit}(convert(AbstractMatrix{T}, A.data))
function convert{T,TA,S,UpLo,IsUnit}(::Type{AbstractMatrix{T}}, A::Triangular{TA,S,UpLo,IsUnit})
    M = convert(AbstractMatrix{T}, A.data)
    Triangular{T,typeof(M),UpLo,IsUnit}(M)
end
function convert{Tret,T,S,UpLo,IsUnit}(::Type{Matrix{Tret}}, A::Triangular{T,S,UpLo,IsUnit})
    B = Array(Tret, size(A, 1), size(A, 1))
    copy!(B, A.data)
    (UpLo == :L ? tril! : triu!)(B)
    if IsUnit
        for i = 1:size(B,1)
            B[i,i] = 1
        end
    end
    B
end
convert{T,S,UpLo,IsUnit}(::Type{Matrix}, A::Triangular{T,S,UpLo,IsUnit}) = convert(Matrix{T}, A)

function full!{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit})
    B = A.data
    (UpLo == :L ? tril! : triu!)(B)
    if IsUnit
        for i = 1:size(A,1)
            B[i,i] = 1
        end
    end
    B
end
full{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = convert(Matrix, A)

fill!(A::Triangular, x) = (fill!(A.data, x); A)

function similar{T,S,UpLo,IsUnit,Tnew}(A::Triangular{T,S,UpLo,IsUnit}, ::Type{Tnew}, dims::Dims)
    dims[1] == dims[2] || throw(ArgumentError("a Triangular matrix must be square"))
    length(dims) == 2 || throw(ArgumentError("a Traigular matrix must have two dimensions"))
    A = similar(A.data, Tnew, dims)
    return Triangular{Tnew, typeof(A), UpLo, IsUnit}(A)
end

getindex{T,S}(A::Triangular{T,S,:L,true}, i::Integer, j::Integer) = i == j ? one(T) : (i > j ? A.data[i,j] : zero(T))
getindex{T,S}(A::Triangular{T,S,:L,false}, i::Integer, j::Integer) = i >= j ? A.data[i,j] : zero(T)
getindex{T,S}(A::Triangular{T,S,:U,true}, i::Integer, j::Integer) = i == j ? one(T) : (i < j ? A.data[i,j] : zero(T))
getindex{T,S}(A::Triangular{T,S,:U,false}, i::Integer, j::Integer) = i <= j ? A.data[i,j] : zero(T)
getindex{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, i::Integer) = ((m, n) = divrem(i - 1, size(A,1)); A[m + 1, n + 1])

istril{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = UpLo == :L
istriu{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = UpLo == :U

transpose{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T, S, UpLo == :U ? :L : :U, IsUnit}(transpose(A.data))
ctranspose{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T, S, UpLo == :U ? :L : :U, IsUnit}(ctranspose(A.data))
transpose!{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T, S, UpLo == :U ? :L : :U, IsUnit}(copytri!(A.data, UpLo == :L ? 'L' : 'U'))
ctranspose!{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T, S, UpLo == :U ? :L : :U, IsUnit}(copytri!(A.data, UpLo == :L ? 'L' : 'U', true))
diag{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = IsUnit ? ones(T, size(A,1)) : diag(A.data)
function big{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit})
    M = big(A.data)
    Triangular{T,typeof(M),UpLo,IsUnit}(M)
end

function (*){T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, x::Number)
    n = size(A,1)
    for j = 1:n
        for i = UpLo == :L ? j:n : 1:j
            A.data[i,j] = i == j & IsUnit ? x : A.data[i,j]*x
        end
    end
    A
end
function (*){T,S,UpLo,IsUnit}(x::Number, A::Triangular{T,S,UpLo,IsUnit})
    n = size(A,1)
    for j = 1:n
        for i = UpLo == :L ? j:n : 1:j
            A.data[i,j] = i == j & IsUnit ? x : x*A.data[i,j]
        end
    end
    A
end

A_mul_B!{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, B::Triangular{T,S,UpLo,IsUnit}) = Triangular{T,S,UpLo,IsUnit}(A*full!(B))
A_mul_B!(A::Tridiagonal, B::Triangular) = A*full!(B)
A_mul_B!(C::AbstractVecOrMat, A::Triangular, B::AbstractVecOrMat) = A_mul_B!(A, copy!(C, B))

A_mul_Bc!(A::Triangular, B::QRCompactWYQ) = A_mul_Bc!(full!(A),B)
A_mul_Bc!(A::Triangular, B::QRPackedQ) = A_mul_Bc!(full!(A),B)
A_mul_Bc!(C::AbstractVecOrMat, A::Triangular, B::AbstractVecOrMat) = A_mul_Bc!(A, copy!(C, B))

#Generic multiplication
*(A::Tridiagonal, B::Triangular) = A_mul_B!(full(A), B)
A_mul_Bc(A::Triangular, B::Union(QRCompactWYQ,QRPackedQ)) = A_mul_Bc(full(A), B)
for func in (:*, :Ac_mul_B, :A_mul_Bc, :/, :A_rdiv_Bc)
    @eval begin
        ($func){TA,TB,SA<:AbstractMatrix,SB<:AbstractMatrix,UpLoA,UpLoB,IsUnitA,IsUnitB}(A::Triangular{TA,SA,UpLoA,IsUnitA}, B::Triangular{TB,SB,UpLoB,IsUnitB}) = ($func)(A, full(B))
        ($func){T,S<:AbstractMatrix,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, B::AbstractVecOrMat) = ($func)(full(A), B)
        ($func){T,S<:AbstractMatrix,UpLo,IsUnit}(A::AbstractMatrix, B::Triangular{T,S,UpLo,IsUnit}) = ($func)(A, full(B))
    end
end

function sqrtm{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit})
    n = size(A, 1)
    R = zeros(T, n, n)
    if UpLo == :U
        for j = 1:n
            (T<:Complex || A[j,j]>=0) ? (R[j,j] = IsUnit ? one(T) : sqrt(A[j,j])) : throw(SingularException(j))
            for i = j-1:-1:1
                r = A[i,j]
                for k = i+1:j-1
                    r -= R[i,k]*R[k,j]
                end
                r==0 || (R[i,j] = r / (R[i,i] + R[j,j]))
            end
        end
        return Triangular{T,S,UpLo,IsUnit}(R)
    else # UpLo == :L #Not the usual case
        return sqrtm(A.').'
    end
end

#Generic solver using naive substitution
function naivesub!{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, b::AbstractVector, x::AbstractVector=b)
    N = size(A, 2)
    N==length(b)==length(x) || throw(DimensionMismatch(""))

    if UpLo == :L #do forward substitution
        for j = 1:N
            x[j] = b[j]
            for k = 1:j-1
                x[j] -= A[j,k] * x[k]
            end
            if !IsUnit
                x[j]/= A[j,j]==0 ? throw(SingularException(j)) : A[j,j]
            end
        end
    elseif UpLo == :U #do backward substitution
        for j = N:-1:1
            x[j] = b[j]
            for k = j+1:1:N
                x[j] -= A[j,k] * x[k]
            end
            if !IsUnit
                x[j]/= A[j,j]==0 ? throw(SingularException(j)) : A[j,j]
            end
        end
    else
        throw(ArgumentError("Unknown UpLo=$(UpLo)"))
    end
    x
end

#Generic eigensystems
eigvals(A::Triangular) = diag(A.data)
det(A::Triangular) = prod(eigvals(A))

function eigvecs{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit})
    N = size(A,1)
    evecs = zeros(T, N, N)
    if IsUnit #Trivial
        return eye(A)
    elseif UpLo == :L #do forward substitution
        for i=1:N
            evecs[i,i] = one(T)
            for j = i+1:N
                for k = i:j-1
                    evecs[j,i] -= A[j,k] * evecs[k,i]
                end
                evecs[j,i] /= A[j,j] - A[i,i]
            end
            evecs[i:N, i] /= norm(evecs[i:N, i])
        end
        evecs
    elseif UpLo == :U #do backward substitution
        for i=1:N
            evecs[i,i] = one(T)
            for j = i-1:-1:1
                for k = j+1:i
                    evecs[j,i] -= A[j,k] * evecs[k,i]
                end
                evecs[j,i] /= A[j,j]-A[i,i]
            end
            evecs[1:i, i] /= norm(evecs[1:i, i])
        end
    else
        throw(ArgumentError("Unknown uplo=$(UpLo)"))
    end
    evecs
end

eigfact(A::Triangular) = Eigen(eigvals(A), eigvecs(A))

#Generic singular systems
for func in (:svd, :svdfact, :svdfact!, :svdvals, :svdvecs)
    @eval begin
        ($func)(A::Triangular) = ($func)(full(A))
    end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# ---
# title: 222. Count Complete Tree Nodes
# id: problem222
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Binary Search, Tree
# link: <https://leetcode.com/problems/count-complete-tree-nodes/description/>
# hidden: true
# ---
#
# Given a **complete** binary tree, count the number of nodes.
#
# **Note:**
#
# **_Definition of a complete binary tree
# from[Wikipedia](http://en.wikipedia.org/wiki/Binary_tree#Types_of_binary_trees):_**
# In a complete binary tree every level, except possibly the last, is completely
# filled, and all nodes in the last level are as far left as possible. It can
# have between 1 and 2h nodes inclusive at the last level h.
#
# **Example:**
#
#
#
#     Input:
#         1
#        / \
#       2   3
#      / \  /
#     4  5 6
#
#     Output: 6
#
#
## @lc code=start
using LeetCode

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


text:
F_0(θ_t) = asin(l_t / l_s * sin(θ_t))
F_1(θ_t) = (l_t / l_s * cos(θ_t)) / (1 - (l_t / l_s * sin(θ_t))^2)^0.5

θ_s(θ_t) = F_0(θ_t)
θ′_s(θ_t, θ′_t) = θ′_t * F_1(θ_t)

x_w(θ_t) = l_t * cos(θ_t) + l_s * cos(θ_s(θ_t)) + x_SB

r_w(θ_t) = SA[x_w(θ_t), 0]

x′_w(θ_t, θ′_t) = -l_t * sin(θ_t) * (θ′_t + θ′_s(θ_t, θ′_t))

r′_w(θ_t, θ′_t) = SA[x′_w(θ_t, θ′_t), 0]

const r_RG = SA[x_OAB, 0]

const r_SB = SA[x_SB, 0]

crossθ(θ, v) = SA[θ*v[2], -θ*v[1]]

r_s_∂l_s(θ_t) = SA[cos(θ_s(θ_t)), l_t / l_s * sin(θ_t)]
const r_s_const = r_SB

r_t_∂l_t(θ_t) = SA[-cos(θ_t), sin(θ_t)]
r_t_const(θ_t) = r_w(θ_t)

r_k_∂l_k(θ_k) = SA[-cos(θ_k), sin(θ_k)]
r_k_const(θ_t) = r_w(θ_t)

x_OA_∂l_OA(γ_OA) = cos(γ_OA)
const y_OA_const = y_OAB
y_OA_∂l_OA(γ_OA) = sin(γ_OA)