text:
module Dim2
import mFEM
using Plots

# calc gauss_quad on standrad triangle (0,0) (1,0) (0,1)
# u: (x::Float64,y::Float64)::Float64
function _gauss(u::Function)
  local points_type = [0 5 2]
  local xyz = [ 9.352701037774565E-01 3.236494811127173E-02 3.23649481112718157E-02;
                7.612981754348137E-01 1.193509122825931E-01 0.11935091228259319E+00;
              - 6.922209654151433E-02 5.346110482707572E-01 0.53461104827075701E+00;
                5.933801991374367E-01 2.033099004312816E-01 0.20330990043128172E+00;
                2.020613940682885E-01 3.989693029658558E-01 0.39896930296585570E+00;
                5.017813831049474E-02 5.932012134282132E-01 0.35662064826129203E+00;
                2.102201653616613E-02 8.074890031597923E-01 0.17148898030404158E+00]
  local w = [ 4.097919300803106E-02;
              1.085536215102866E-01;
              2.781018986881812E-03;
              1.779689321422668E-01;
              2.314486047444677E-01;
              3.140226717732234E-01;
              1.242459578348437e-01]./[3,3,3,3,3,6,6]

  local ptorder3 = [1 1 2 2 3 3;
                    2 3 1 3 1 2]
  local ptorder2 = [1 2 2;
                    2 1 3]

  local v = 0.0
  for i=1:5, j=1:size(ptorder2,2)
    x = xyz[i, ptorder2[1,j]]
    y = xyz[i, ptorder2[2,j]]
    v += w[i] * u(x,y)
  end
  for i=6:7, j=1:size(ptorder3,2)
    x = xyz[i, ptorder3[1,j]]
    y = xyz[i, ptorder3[2,j]]
    v += w[i] * u(x,y)
  end

  return v/2 # v * measure of triangle{(0,0),(1,0),(0,1)}
end

function uniform_trianglar_mesh(x0::Vector,x1::Vector,
    N::Vector{Int})::mFEM.Mesh
  @assert size(x0,1) == size(x1,1) == 2
  @assert x0[1] < x1[1]
  @assert x0[2] < x1[2]
  @assert size(N,1) == 2
  @assert N[1] > 1 && N[2] > 1

  P = zeros((N[1]+1) * (N[2]+1),2)

  local h = (x1 - x0) ./ N
  for i in 1:(N[2]+1), j in 1:(N[1]+1)
    x,y = x0 + h .* [j-1,i-1]
    k = (i-1) * (N[1] + 1) + j
    P[k,:] = [x,y]
  end

  T = zeros(Int,2*N[1]*N[2],3)
  for i=1:N[2], j=1:N[1]
    k = (i-1) * (N[1] + 1) + j
    tk = k - (i - 1)
    T[tk,:] = [k,k+1,k+N[1]+1]
    T[tk+N[1]*N[2],:] = [k + 1 + N[1] + 1, k + N[1] + 1, k + 1]
  end

  # calc gauss_quad on the nth element of mesh
  # u: (x::Tuple)::Float64
  # n index of basis element
  function gauss_quad(u::Function,n::Int)
    local x1 = P[T[n,1],:]
    local x2 = P[T[n,2],:]
    local x3 = P[T[n,3],:]
    local O = [x2-x1 x3-x1]

    local a = x2-x1
    local b = x3-x1
    measure = sqrt(sum(a.*a)*sum(b.*b) - sum(a.*b)^2)

    return measure * _gauss((x,y)-> begin
        return u(O * [x,y] + x1)
    end)
  end

  return mFEM.Mesh(P,size(P,1),T,size(T,1),gauss_quad)
end

function Dϕ_reference_trangular(k::Int,xt::Tuple,α::Int)
  @assert α ∈ 1:3
  @assert k ∈ 0:1
  @assert length(xt) == 2

  x,y = xt
  if α == 1 # (0,0)
    if k == 0
      return 1 - x - y
    elseif k == 1
      return [-1,-1]
    end
  elseif α == 2 # (1,0)
    if k == 0
      return x
    elseif k == 1
      return [1,0]
    end
  elseif α == 3  # (0,1)
    if k == 0
      return y
    elseif k == 1
      return [0,1]
    end
  end
end

function reference2local(Dϕ_reference::Function;Pb::Array,Tb::Array,Nlb::Int)::Function
  # kth deriavtive of basis function at x
  function Dϕ(k::Int, xt::Tuple, n::Int, α::Int)
    @assert α ∈ 1:Nlb
    @assert k ∈ 0:1
    @assert length(xt) == 2

    local x1 = P[T[n,1],:]
    local x2 = P[T[n,2],:]
    local x3 = P[T[n,3],:]
    local O = [x2-x1 x3-x1]

    xhat = O \ (xt - x1)
    return Dϕ_reference(k,xhat,α)
  end

  return Dϕ
end

function trangular_elem(mesh::mFEM.Mesh)::mFEM.Basis
  Pb = copy(mesh.P)
  Tb = copy(mesh.T)
  Nlb = 3

  return mFEM.Basis(size(Pb,1),Pb,Tb,size(Tb,1),Nlb,
    reference2local(Dϕ_reference_trangular;Pb=Pb,Tb=Tb,Nlb=Nlb))
end

function test(n)
  tm = mFEM.Dim2.uniform_trianglar_mesh([0,0],[1,1],[n,n])
  u(x) = x[1]*x[2] *exp(x[1]^2 + x[2]^2)
  v = 0.0
  for n=1:size(tm.T,1)
    v += tm.Integral(u,n)
  end
  return abs(v - (0.5*(exp(1) - 1))^2)
end

end # Dim2


text:
struct TanhModulation <: InfectionModulation
  weight_detected::Float64
  weight_deaths::Float64
  weight_days::Float64
  loc::Float64
  scale::Float64
  limit_value::Float64

  TanhModulation(;weight_detected::Real=0, weight_deaths::Real=0, weight_days::Real=0, loc::Real=0, scale::Real=0, limit_value::Real=0) =
    0 <= limit_value <= 1 ? new(weight_detected, weight_deaths, weight_days, loc, scale, limit_value) : error("limit_value must be from 0 to 1, got $limit_value")
end

function tanh_modulation(x::Real, loc::Real, scale::Real, lhs_limit::Real, rhs_limit::Real)
  t_11 = tanh((x - loc)/scale)
  t_01 = (t_11+1)/2
  t_01 * (rhs_limit-lhs_limit) + lhs_limit
end

function (f::TanhModulation)(state::AbstractSimState, ::AbstractSimParams, event::Event)
  @assert kind(event) == TransmissionEvent

  ck = contactkind(event)
  if ConstantKernelContact !== ck && AgeCouplingContact !== ck
    return true # do not affect other types of contact than "outer" ones
  end

  num_days = time(state)
  num_detected = numdetected(state)
  num_deaths = numdead(state)

  fear = num_detected * f.weight_detected + num_deaths * f.weight_deaths + num_days * f.weight_days
  t = tanh_modulation(fear, f.loc, f.scale, 1.0, f.limit_value)
  rand(state.rng) < t
end

struct IncreasingTanhModulation <: InfectionModulation
  weight_detected::Float64
  weight_deaths::Float64
  weight_days::Float64
  loc::Float64
  scale::Float64
  initial_value::Float64

  IncreasingTanhModulation(;weight_detected::Real=0, weight_deaths::Real=0, weight_days::Real=0, loc::Real=0, scale::Real=0, initial_value::Real=0) =
    0 <= initial_value <= 1 ? new(weight_detected, weight_deaths, weight_days, loc, scale, initial_value) : error("initial_value must be from 0 to 1, got $initial_value")
end

function (f::IncreasingTanhModulation)(state::AbstractSimState, ::AbstractSimParams, event::Event)
  @assert kind(event) == TransmissionEvent

  ck = contactkind(event)
  if ConstantKernelContact !== ck && AgeCouplingContact !== ck
    return true # do not affect other types of contact than "outer" ones
  end

  num_days = time(state)
  num_detected = numdetected(state)
  num_deaths = numdead(state)

  fear = num_detected * f.weight_detected + num_deaths * f.weight_deaths + num_days * f.weight_days
  t = tanh_modulation(fear, f.loc, f.scale, f.initial_value, 1.0)
  rand(state.rng) < t
end

# This all to avoid using @eval and others
const modulations = Dict{String, Type{T} where T}(
    "TanhModulation" => TanhModulation,
    "IncreasingTanhModulation" => IncreasingTanhModulation
)

make_infection_modulation(name::AbstractString; args...) = modulations[name](;args...)

ext: jl


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

function diff_tree(repo::GitRepo, tree::GitTree, pathspecs::AbstractString=""; cached::Bool=false)
    emptypathspec = isempty(pathspecs)
    diff_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    if !emptypathspec
        sa = StrArrayStruct(pathspecs)
        diff_opts = DiffOptionsStruct(pathspec = sa)
    end
    try
        if cached
            @check ccall((:git_diff_tree_to_index, :libgit2), Cint,
                          (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{DiffOptionsStruct}),
                           diff_ptr_ptr, repo.ptr, tree.ptr, C_NULL,  emptypathspec ? C_NULL : Ref(diff_opts))
        else
            @check ccall((:git_diff_tree_to_workdir_with_index, :libgit2), Cint,
                          (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}, Ptr{DiffOptionsStruct}),
                           diff_ptr_ptr, repo.ptr, tree.ptr, emptypathspec ? C_NULL : Ref(diff_opts))
        end
    finally
        !emptypathspec && finalize(sa)
    end
    return GitDiff(diff_ptr_ptr[])
end

function diff_tree(repo::GitRepo, oldtree::GitTree, newtree::GitTree)
    diff_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_diff_tree_to_tree, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{DiffOptionsStruct}),
                   diff_ptr_ptr, repo.ptr, oldtree.ptr, newtree.ptr, C_NULL)
    return GitDiff(diff_ptr_ptr[])
end

function Base.count(diff::GitDiff)
    return ccall((:git_diff_num_deltas, :libgit2), Cint, (Ptr{Void},), diff.ptr)
end

function Base.getindex(diff::GitDiff, i::Csize_t)
    delta_ptr = ccall((:git_diff_get_delta, :libgit2), Ptr{Void},
                       (Ptr{Void}, Csize_t), diff.ptr, i-1)
    delta_ptr == C_NULL && return nothing
    return unsafe_load(convert(Ptr{DiffDelta}, delta_ptr), 1)
end
Base.getindex(diff::GitDiff, i::Int) = getindex(diff, Csize_t(i))

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
"""
    BlockSlice(indices)

Represent an AbstractUnitRange of indices that attaches a block.

Upon calling `to_indices()`, Blocks are converted to BlockSlice objects to represent
the indices over which the Block spans.

This mimics the relationship between `Colon` and `Base.Slice`.
"""
struct BlockSlice{BB} <: AbstractUnitRange{Int}
    block::BB
    indices::UnitRange{Int}
end

Block(bs::BlockSlice{<:Block}) = bs.block

for f in (:axes, :unsafe_indices, :axes1, :first, :last, :size, :length,
          :unsafe_length, :start)
    @eval $f(S::BlockSlice) = $f(S.indices)
end

getindex(S::BlockSlice, i::Integer) = getindex(S.indices, i)
show(io::IO, r::BlockSlice) = print(io, "BlockSlice(", r.block, ",", r.indices, ")")
next(S::BlockSlice, s) = next(S.indices, s)
done(S::BlockSlice, s) = done(S.indices, s)

function _unblock(cum_sizes, I::Tuple{Block{1, T},Vararg{Any}}) where {T}
    B = first(I)
    b = first(B.n)

    range = cum_sizes[b]:cum_sizes[b + 1] - 1

    BlockSlice(B, range)
end

function _unblock(cum_sizes, I::Tuple{BlockRange{1,R}, Vararg{Any}}) where {R}
    B = first(I)
    b_start = first(first(B.indices))
    b_stop = last(first(B.indices))

    range = cum_sizes[b_start]:cum_sizes[b_stop + 1] - 1

    BlockSlice(B, range)
end

_sub_cumul_sizes(cs, inds) = _sub_cumul_sizes(cs, inds[1], tail(inds))
_sub_cumul_sizes(::Tuple{}, ::Tuple{}) = ()

function _sub_cumul_sizes(cs, inds1::BlockSlice{Block{1,Int}}, inds)
    B = Int(inds1.block)
    ret = view(cs[1], B:B+1)
    (ret .- ret[1] .+ 1, _sub_cumul_sizes(tail(cs), inds)...)
end

function _sub_cumul_sizes(cs, inds1::BlockSlice{BlockRange{1,Tuple{UnitRange{Int}}}}, inds)
    ret = view(cs[1], inds1.block.indices[1][1]:(inds1.block.indices[1][end]+1))
    (ret .- ret[1] .+ 1, _sub_cumul_sizes(tail(cs), inds)...)
end

blocksizes(V::SubArray) = BlockSizes(_sub_cumul_sizes(cumulsizes(parent(V)), parentindices(V)))

"""
    unblock(block_sizes, inds, I)

Returns the indices associated with a block as a `BlockSlice`.
"""
function unblock(A::AbstractArray{T,N}, inds, I) where {T, N}
    if length(inds) == 0
        # Allow `ones(2)[Block(1)[1:1], Block(1)[1:1]]` which is
        # similar to `ones(2)[1:1, 1:1]`.
        _unblock(Base.OneTo(2), I)
    else
        _unblock(cumulsizes(A, N - length(inds) + 1), I)
    end
end

to_index(::Block) = throw(ArgumentError("Block must be converted by to_indices(...)"))
to_index(::BlockRange) = throw(ArgumentError("BlockRange must be converted by to_indices(...)"))

@inline to_indices(A, inds, I::Tuple{Block{1}, Vararg{Any}}) =
    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)

# splat out higher dimensional blocks
# this mimics view of a CartesianIndex
@inline to_indices(A, inds, I::Tuple{Block, Vararg{Any}}) =
    to_indices(A, inds, (Block.(I[1].n)..., tail(I)...))

@inline to_indices(A, inds, I::Tuple{BlockRange{1,R}, Vararg{Any}}) where R =
    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)

# splat out higher dimensional blocks
# this mimics view of a CartesianIndex
@inline to_indices(A, inds, I::Tuple{BlockRange, Vararg{Any}}) =
    to_indices(A, inds, (BlockRange.(tuple.(I[1].indices))..., tail(I)...))

# In 0.7, we need to override to_indices to avoid calling linearindices
@inline to_indices(A, I::Tuple{Block, Vararg{Any}}) =
    to_indices(A, axes(A), I)

@inline to_indices(A, I::Tuple{BlockRange, Vararg{Any}}) =
    to_indices(A, axes(A), I)

# The first argument for `reindex` is removed as of
# https://github.com/JuliaLang/julia/pull/30789 in Julia `Base`.  So,
# we define 2-arg `reindex` for Julia 1.2 and later.
if VERSION >= v"1.2-"

# BlockSlices map the blocks and the indices
# this is loosely based on Slice reindex in subarray.jl
reindex(idxs::Tuple{BlockSlice{<:BlockRange}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{<:BlockRange}, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (BlockSlice(BlockRange(idxs[1].block.indices[1][Int.(subidxs[1].block)]),
                                            idxs[1].indices[subidxs[1].indices]),
                                 reindex(tail(idxs), tail(subidxs))...))

reindex(idxs::Tuple{BlockSlice{BlockRange{1,Tuple{UnitRange{Int}}}}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{Block{1,Int}}, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (BlockSlice(Block(idxs[1].block.indices[1][Int(subidxs[1].block)]),
                                            idxs[1].indices[subidxs[1].indices]),
                                 reindex(tail(idxs), tail(subidxs))...))

function reindex(idxs::Tuple{BlockSlice{Block{1,Int}}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{Block{1,Int}}, Vararg{Any}})
    (idxs[1], reindex(tail(idxs), tail(subidxs))...)
end

else  # if VERSION >= v"1.2-"

reindex(V, idxs::Tuple{BlockSlice{<:BlockRange}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{<:BlockRange}, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (BlockSlice(BlockRange(idxs[1].block.indices[1][Int.(subidxs[1].block)]),
                                            idxs[1].indices[subidxs[1].indices]),
                                    reindex(V, tail(idxs), tail(subidxs))...))

reindex(V, idxs::Tuple{BlockSlice{BlockRange{1,Tuple{UnitRange{Int}}}}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{Block{1,Int}}, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (BlockSlice(Block(idxs[1].block.indices[1][Int(subidxs[1].block)]),
                                            idxs[1].indices[subidxs[1].indices]),
                                    reindex(V, tail(idxs), tail(subidxs))...))

function reindex(V, idxs::Tuple{BlockSlice{Block{1,Int}}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{Block{1,Int}}, Vararg{Any}})
    subidxs[1].block == Block(1) || throw(BoundsError(V, subidxs[1].block))
    (idxs[1], reindex(V, tail(idxs), tail(subidxs))...)
end

end  # if VERSION >= v"1.2-"

#################
# support for pointers
#################

const BlockOrRangeIndex = Union{RangeIndex, BlockSlice}

function unsafe_convert(::Type{Ptr{T}},
                        V::SubArray{T, N, BlockArray{T,N,AT,BS}, NTuple{N, BlockSlice{Block{1,Int}}}}) where {AT <: AbstractArray{<:AbstractArray{T,N},N}, BS <: AbstractBlockSizes{N}} where {T,N}
    unsafe_convert(Ptr{T}, parent(V).blocks[Int.(Block.(parentindices(V)))...])
end

unsafe_convert(::Type{Ptr{T}}, V::SubArray{T,N,PseudoBlockArray{T,N,AT},<:Tuple{Vararg{BlockOrRangeIndex}}}) where {T,N,AT} =
    unsafe_convert(Ptr{T}, V.parent) + (Base.first_index(V)-1)*sizeof(T)


text:
### Single type of emitters
function get_measurement(
    num_samples,
    receiver::AbstractReceiver,
    emitters::Tuple{Vararg{AbstractEmitter{T}}},
    manifold::AbstractManifold{1} = IdealManifold(),
    rng = Random.GLOBAL_RNG
) where T <: AbstractFloat
    signal = StructArray{Complex{T}}(undef, num_samples)
    get_measurement!(signal, receiver, emitters, manifold, rng)
end

function get_measurement(
    num_samples,
    receiver::AbstractReceiver,
    emitters::Tuple{Vararg{AbstractEmitter{T}}},
    manifold::AbstractManifold{N},
    rng = Random.GLOBAL_RNG
) where {N, T <: AbstractFloat}
    signal = StructArray{Complex{T}}(undef, num_samples, N)
    get_measurement!(signal, receiver, emitters, manifold, rng)
end

function get_measurement!(
    signal::Union{AbstractMatrix{Complex{T}}, AbstractVector{Complex{T}}},
    receiver::AbstractReceiver,
    emitters::Tuple{Vararg{AbstractEmitter{T}}},
    manifold::AbstractManifold{N} = IdealManifold(),
    rng = Random.GLOBAL_RNG
) where {N, T <: AbstractFloat}
    num_samples = get_num_samples(signal)
    existing_emitters = filteremitters(get_existence, emitters)
    emitter_signals = gen_signal!.(
        existing_emitters,
        get_sampling_frequency(receiver),
        get_intermediate_frequency(receiver),
        get_noise_density(receiver),
        num_samples,
        Ref(rng)
    )
    steer_vecs = get_steer_vecs(T, existing_emitters, receiver, manifold)
    randn!(rng, signal)
    signal .*= T(get_noise_std(receiver))
    foreach(emitter_signals, steer_vecs) do emitter_signal, steer_vec
        signal .+= emitter_signal .* transpose(steer_vec)
    end
    measurement = multiply_gpmc_with_signal!(
        signal,
        get_gain_phase_mism_crosstalk(receiver)
    )

    next_receiver = propagate(receiver, num_samples, rng)
    next_emitters = propagate.(
        emitters,
        num_samples,
        get_intermediate_frequency(receiver),
        get_sampling_frequency(receiver),
        get_noise_density(receiver),
        Ref(rng)
    )
    signal, next_receiver, next_emitters
end

### Two types of emitters
function get_measurement(
    num_samples,
    receiver::AbstractReceiver,
    emitters1::Tuple{Vararg{AbstractEmitter{T}}},
    emitters2::Tuple{Vararg{AbstractEmitter{T}}},
    manifold::AbstractManifold{1} = IdealManifold(),
    rng = Random.GLOBAL_RNG
) where T <: AbstractFloat
    signal = StructArray{Complex{T}}(undef, num_samples)
    get_measurement!(signal, receiver, emitters1, emitters2, manifold, rng)
end

function get_measurement(
    num_samples,
    receiver::AbstractReceiver,
    emitters1::Tuple{Vararg{AbstractEmitter{T}}},
    emitters2::Tuple{Vararg{AbstractEmitter{T}}},
    manifold::AbstractManifold{N},
    rng = Random.GLOBAL_RNG
) where {N, T <: AbstractFloat}
    signal = StructArray{Complex{T}}(undef, num_samples, N)
    get_measurement!(signal, receiver, emitters1, emitters2, manifold, rng)
end

function get_measurement!(
    signal::Union{AbstractMatrix{Complex{T}}, AbstractVector{Complex{T}}},
    receiver::AbstractReceiver,
    emitters1::Tuple{Vararg{AbstractEmitter{T}}},
    emitters2::Tuple{Vararg{AbstractEmitter{T}}},
    manifold::AbstractManifold{N} = IdealManifold(),
    rng = Random.GLOBAL_RNG
) where {N, T <: AbstractFloat}
    num_samples = get_num_samples(signal)
    existing_emitters1 = filteremitters(get_existence, emitters1)
    existing_emitters2 = filteremitters(get_existence, emitters2)
    emitter_signals1 = gen_signal!.(
        existing_emitters1,
        get_sampling_frequency(receiver),
        get_intermediate_frequency(receiver),
        get_noise_density(receiver),
        num_samples,
        Ref(rng)
    )
    emitter_signals2 = gen_signal!.(
        existing_emitters2,
        get_sampling_frequency(receiver),
        get_intermediate_frequency(receiver),
        get_noise_density(receiver),
        num_samples,
        Ref(rng)
    )
    steer_vecs1 = get_steer_vecs(T, existing_emitters1, receiver, manifold)
    steer_vecs2 = get_steer_vecs(T, existing_emitters2, receiver, manifold)
    randn!(rng, signal)
    signal .*= T(get_noise_std(receiver))
    foreach(emitter_signals1, steer_vecs1) do emitter_signal, steer_vec
        signal .+= emitter_signal .* transpose(steer_vec)
    end
    foreach(emitter_signals2, steer_vecs2) do emitter_signal, steer_vec
        signal .+= emitter_signal .* transpose(steer_vec)
    end
    measurement = multiply_gpmc_with_signal!(
        signal,
        get_gain_phase_mism_crosstalk(receiver)
    )

    next_receiver = propagate(receiver, num_samples, rng)
    next_emitters1 = propagate.(
        emitters1,
        num_samples,
        get_intermediate_frequency(receiver),
        get_sampling_frequency(receiver),
        get_noise_density(receiver),
        Ref(rng)
    )
    next_emitters2 = propagate.(
        emitters2,
        num_samples,
        get_intermediate_frequency(receiver),
        get_sampling_frequency(receiver),
        get_noise_density(receiver),
        Ref(rng)
    )
    measurement, next_receiver, next_emitters1, next_emitters2
end

### Three types of emitters
function get_measurement(
    num_samples,
    receiver::AbstractReceiver,
    emitters1::Tuple{Vararg{AbstractEmitter{T}}},
    emitters2::Tuple{Vararg{AbstractEmitter{T}}},
    emitters3::Tuple{Vararg{AbstractEmitter{T}}},
    manifold::AbstractManifold{1} = IdealManifold(),
    rng = Random.GLOBAL_RNG
) where {T <: AbstractFloat}
    signal = StructArray{Complex{T}}(undef, num_samples)
    get_measurement!(signal, receiver, emitters1, emitters2, emitters3, manifold, rng)
end

function get_measurement(
    num_samples,
    receiver::AbstractReceiver,
    emitters1::Tuple{Vararg{AbstractEmitter{T}}},
    emitters2::Tuple{Vararg{AbstractEmitter{T}}},
    emitters3::Tuple{Vararg{AbstractEmitter{T}}},
    manifold::AbstractManifold{N},
    rng = Random.GLOBAL_RNG
) where {N, T <: AbstractFloat}
    signal = StructArray{Complex{T}}(undef, num_samples, N)
    get_measurement!(signal, receiver, emitters1, emitters2, emitters3, manifold, rng)
end

function get_measurement!(
    signal::Union{AbstractMatrix{Complex{T}}, AbstractVector{Complex{T}}},
    receiver::AbstractReceiver,
    emitters1::Tuple{Vararg{AbstractEmitter{T}}},
    emitters2::Tuple{Vararg{AbstractEmitter{T}}},
    emitters3::Tuple{Vararg{AbstractEmitter{T}}},
    manifold::AbstractManifold{N} = IdealManifold(),
    rng = Random.GLOBAL_RNG
) where {N, T <: AbstractFloat}
    num_samples = get_num_samples(signal)
    existing_emitters1 = filteremitters(get_existence, emitters1)
    existing_emitters2 = filteremitters(get_existence, emitters2)
    existing_emitters3 = filteremitters(get_existence, emitters3)
    emitter_signals1 = gen_signal!.(
        existing_emitters1,
        get_sampling_frequency(receiver),
        get_intermediate_frequency(receiver),
        get_noise_density(receiver),
        num_samples,
        Ref(rng)
    )
    emitter_signals2 = gen_signal!.(
        existing_emitters2,
        get_sampling_frequency(receiver),
        get_intermediate_frequency(receiver),
        get_noise_density(receiver),
        num_samples,
        Ref(rng)
    )
    emitter_signals3 = gen_signal!.(
        existing_emitters3,
        get_sampling_frequency(receiver),
        get_intermediate_frequency(receiver),
        get_noise_density(receiver),
        num_samples,
        Ref(rng)
    )
    steer_vecs1 = get_steer_vecs(T, existing_emitters1, receiver, manifold)
    steer_vecs2 = get_steer_vecs(T, existing_emitters2, receiver, manifold)
    steer_vecs3 = get_steer_vecs(T, existing_emitters3, receiver, manifold)
    randn!(rng, signal)
    signal .*= T(get_noise_std(receiver))
    foreach(emitter_signals1, steer_vecs1) do emitter_signal, steer_vec
        signal .+= emitter_signal .* transpose(steer_vec)
    end
    foreach(emitter_signals2, steer_vecs2) do emitter_signal, steer_vec
        signal .+= emitter_signal .* transpose(steer_vec)
    end
    foreach(emitter_signals3, steer_vecs3) do emitter_signal, steer_vec
        signal .+= emitter_signal .* transpose(steer_vec)
    end
    measurement = multiply_gpmc_with_signal!(
        signal,
        get_gain_phase_mism_crosstalk(receiver)
    )

    next_receiver = propagate(receiver, num_samples, rng)
    next_emitters1 = propagate.(
        emitters1,
        num_samples,
        get_intermediate_frequency(receiver),
        get_sampling_frequency(receiver),
        get_noise_density(receiver),
        Ref(rng)
    )
    next_emitters2 = propagate.(
        emitters2,
        num_samples,
        get_intermediate_frequency(receiver),
        get_sampling_frequency(receiver),
        get_noise_density(receiver),
        Ref(rng)
    )
    next_emitters3 = propagate.(
        emitters3,
        num_samples,
        get_intermediate_frequency(receiver),
        get_sampling_frequency(receiver),
        get_noise_density(receiver),
        Ref(rng)
    )
    signal, next_receiver, next_emitters1, next_emitters2, next_emitters3
end

### Four types of emitter
function get_measurement(
    num_samples,
    receiver::AbstractReceiver,
    emitters1::Tuple{Vararg{AbstractEmitter{T}}},
    emitters2::Tuple{Vararg{AbstractEmitter{T}}},
    emitters3::Tuple{Vararg{AbstractEmitter{T}}},
    emitters4::Tuple{Vararg{AbstractEmitter{T}}},
    manifold::AbstractManifold{N},
    rng = Random.GLOBAL_RNG
) where {N, T <: AbstractFloat}
    signal = StructArray{Complex{T}}(undef, num_samples, N)
    get_measurement!(signal, receiver, emitters1, emitters2, emitters3, emitters4, manifold, rng)
end

function get_measurement!(
    signal::Union{AbstractMatrix{Complex{T}}, AbstractVector{Complex{T}}},
    receiver::AbstractReceiver,
    emitters1::Tuple{Vararg{AbstractEmitter{T}}},
    emitters2::Tuple{Vararg{AbstractEmitter{T}}},
    emitters3::Tuple{Vararg{AbstractEmitter{T}}},
    emitters4::Tuple{Vararg{AbstractEmitter{T}}},
    manifold::AbstractManifold{N} = IdealManifold(),
    rng = Random.GLOBAL_RNG
) where {N, T <: AbstractFloat}
    num_samples = get_num_samples(signal)
    existing_emitters1 = filteremitters(get_existence, emitters1)
    existing_emitters2 = filteremitters(get_existence, emitters2)
    existing_emitters3 = filteremitters(get_existence, emitters3)
    existing_emitters4 = filteremitters(get_existence, emitters4)
    emitter_signals1 = gen_signal!.(
        existing_emitters1,
        get_sampling_frequency(receiver),
        get_intermediate_frequency(receiver),
        get_noise_density(receiver),
        num_samples,
        Ref(rng)
    )
    emitter_signals2 = gen_signal!.(
        existing_emitters2,
        get_sampling_frequency(receiver),
        get_intermediate_frequency(receiver),
        get_noise_density(receiver),
        num_samples,
        Ref(rng)
    )
    emitter_signals3 = gen_signal!.(
        existing_emitters3,
        get_sampling_frequency(receiver),
        get_intermediate_frequency(receiver),
        get_noise_density(receiver),
        num_samples,
        Ref(rng)
    )
    emitter_signals4 = gen_signal!.(
        existing_emitters4,
        get_sampling_frequency(receiver),
        get_intermediate_frequency(receiver),
        get_noise_density(receiver),
        num_samples,
        Ref(rng)
    )
    steer_vecs1 = get_steer_vecs(T, existing_emitters1, receiver, manifold)
    steer_vecs2 = get_steer_vecs(T, existing_emitters2, receiver, manifold)
    steer_vecs3 = get_steer_vecs(T, existing_emitters3, receiver, manifold)
    steer_vecs4 = get_steer_vecs(T, existing_emitters4, receiver, manifold)
    randn!(rng, signal)
    signal .*= T(get_noise_std(receiver))
    foreach(emitter_signals1, steer_vecs1) do emitter_signal, steer_vec
        signal .+= emitter_signal .* transpose(steer_vec)
    end
    foreach(emitter_signals2, steer_vecs2) do emitter_signal, steer_vec
        signal .+= emitter_signal .* transpose(steer_vec)
    end
    foreach(emitter_signals3, steer_vecs3) do emitter_signal, steer_vec
        signal .+= emitter_signal .* transpose(steer_vec)
    end
    foreach(emitter_signals4, steer_vecs4) do emitter_signal, steer_vec
        signal .+= emitter_signal .* transpose(steer_vec)
    end
    measurement = multiply_gpmc_with_signal!(
        signal,
        get_gain_phase_mism_crosstalk(receiver)
    )

    next_receiver = propagate(receiver, num_samples, rng)
    next_emitters1 = propagate.(
        emitters1,
        num_samples,
        get_intermediate_frequency(receiver),
        get_sampling_frequency(receiver),
        get_noise_density(receiver),
        Ref(rng)
    )
    next_emitters2 = propagate.(
        emitters2,
        num_samples,
        get_intermediate_frequency(receiver),
        get_sampling_frequency(receiver),
        get_noise_density(receiver),
        Ref(rng)
    )
    next_emitters3 = propagate.(
        emitters3,
        num_samples,
        get_intermediate_frequency(receiver),
        get_sampling_frequency(receiver),
        get_noise_density(receiver),
        Ref(rng)
    )
    next_emitters4 = propagate.(
        emitters4,
        num_samples,
        get_intermediate_frequency(receiver),
        get_sampling_frequency(receiver),
        get_noise_density(receiver),
        Ref(rng)
    )
    signal, next_receiver, next_emitters1, next_emitters2, next_emitters3, next_emitters4
end

function convert_complex_or_real(::Type{T}, a::Complex) where T <: Union{Float32, Float64}
    Complex{T}(a)
end

function convert_complex_or_real(::Type{T}, a) where T <: Union{Float32, Float64}
    T(a)
end

function get_steer_vecs(::Type{T}, emitters, receiver, manifold) where T <: Union{Float32, Float64}
    map(emitters) do emitter
        convert_complex_or_real.(T,
            get_steer_vec(manifold, emitter, get_attitude(receiver))
        )
    end
end

@inline function get_num_samples(signal)
    size(signal, 1)
end

function multiply_gpmc_with_signal!(signal, gpmc::Number)
    signal .*= gpmc
end

function multiply_gpmc_with_signal!(signal::AbstractMatrix{T}, gpmc::SMatrix) where T
    A_mul_Bt!(signal, gpmc)
end

function A_mul_Bt!(A::AbstractMatrix, B::SMatrix{N, N, T}) where {N, T}
    size(A, 2) == N || throw(DimensionMismatch("2nd dimension of A must have same dimension as of B"))
    @inbounds for i = 1:size(A, 1)
        v = SVector{N, T}(view(A, i, :))
        A[i, :] .= B * v
    end
    A
end

function A_mul_Bt!(A::StructArray, B::SMatrix{N, N, T}) where {N, T}
    size(A, 2) == N || throw(DimensionMismatch("2nd dimension of A must have same dimension as of B"))
    B_re = real.(B)
    B_im = imag.(B)
    @inbounds for i = 1:size(A, 1)
        v_re = SVector{N, T}(view(A.re, i, :))
        v_im = SVector{N, T}(view(A.im, i, :))
        A.re[i, :] .= B_re * v_re - B_im * v_im
        A.im[i, :] .= B_re * v_im + B_im * v_re
    end
    A
end


text:
# License is MIT: https://github.com/JuliaString/HTML_Entities/LICENSE.md
#
# Mapping from HTML entities to the corresponding Unicode codepoint.

println("Running HTML entity build in ", pwd())

using StrTables

VER = UInt32(1)

const inpname = "htmlnames.jl"

include(inpname)

const disp = [false]

const fname = "html.dat"
const datapath = joinpath(@__DIR__, "..", "data")

const empty_str = ""

function make_tables()
    symnam = String[]
    symval = Vector{UInt32}[]

    for (nam, val) in htmlonechar
        push!(symnam, nam)
        push!(symval, [val])
    end
    for (nam, val) in htmlnonbmp
        push!(symnam, nam)
        push!(symval, [0x10000+val])
    end
    for (nam, val) in htmltwochar
        push!(symnam, nam)
        push!(symval, UInt32[val...])
    end

    # We want to build a table of all the names, sort them, then create a StrTable out of them
    srtnam = sortperm(symnam)
    srtval = symval[srtnam] # Values, sorted the same as srtnam

    # BMP characters
    l16 = Tuple{UInt16, UInt16}[]
    # non-BMP characters (in range 0x10000 - 0x1ffff)
    l32 = Tuple{UInt16, UInt16}[]
    # two characters packed into UInt32, first character in high 16-bits
    l2c = Tuple{UInt32, UInt16}[]

    for i in eachindex(srtnam)
        chrs = srtval[i]
        len = length(chrs)
        len > 2 && error("Too long sequence of characters $chrs")
        ch1 = chrs[1]
        if len == 2
            ch2 = chrs[end]
            (ch1 > 0x0ffff || ch2 > 0x0ffff) &&
                error("Character $ch1 or $ch2 > 0xffff")
            push!(l2c, (ch1<<16 | ch2, i))
        elseif ch1 > 0x1ffff
            error("Character $ch1 too large")
        elseif ch1 > 0x0ffff
            push!(l32, (ch1%UInt16, i))
        else
            push!(l16, (ch1%UInt16, i))
        end
    end

    # We now have 3 vectors, for single BMP characters, for non-BMP characters, and for 2 BMP chars
    # each has the value and a index into the name table
    # We need to create a vector the same size as the name table, that gives the index
    # of into one of the three tables, in order to go from names to 1 or 2 output characters
    # We also need, for each of the 3 tables, a sorted vector that goes from the indices
    # in each table to the index into the name table (so that we can find multiple names for
    # each character)

    indvec = create_vector(UInt16, length(srtnam))
    vec16, ind16, base32 = sortsplit!(indvec, l16, 0)
    vec32, ind32, base2c = sortsplit!(indvec, l32, base32)
    vec2c, ind2c, basefn = sortsplit!(indvec, l2c, base2c)

    (VER, string(now()), "loaded from $inpname",
     base32%UInt32, base2c%UInt32, StrTable(symnam[srtnam]), indvec,
     vec16, ind16, vec32, ind32, vec2c, ind2c)
end

savfile = joinpath(datapath, fname)
if isfile(savfile)
    println("Tables already exist")
else
    println("Creating tables")
    try
        tup = make_tables()
        println("Saving tables to ", savfile)
        StrTables.save(savfile, tup)
        println("Done")
    catch ex
        println("Error creating HTML tables")
        println(sprint(showerror, ex, catch_backtrace()))
    end
end


text:
# Abif format
# ===========
#
# Reader of the ABIF file format.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioJulia.jl/blob/master/LICENSE.md

module ABIF

import BioCore

include("reader.jl")

end


text:
### Exports ###

export CPUContainer
export root, transient

### Types ###

# InactiveContainer doesn't do anything, so it doesn't need a transient
mutable struct InactiveContainer{R} <: AbstractContainer
  root::R
end
# CPUContainer doesn't store data anywhere special, so it doesn't need a transient
mutable struct CPUContainer{R} <: AbstractContainer
  root::R
end

### Methods ###

# Gets the root of a container
root(cont::AbstractContainer) = cont.root

# Gets the transient of a container
# TODO: Replace with getfield?
transient(cont::AbstractContainer) = cont.transient
transient(cont::InactiveContainer) = root(cont)
transient(cont::CPUContainer) = root(cont)

==(cont1::CPUContainer, cont2::CPUContainer) =
  root(cont1) == root(cont2)


text:
function get_zz(components::AbstractVector, samples::AbstractVector)
    ret = [ pdf(comp, sample) for comp in components, sample in samples ]
    ret ./= sum(ret, dims=1)
end

function emstep(mixture::MT, samples::AbstractVector, weights::AbstractWeights;
                zz_pre = get_zz(components(mixture), samples),
                constructor = MixtureModel, kwargs...) where {MT <: AbstractMixtureModel}
    comps = components(mixture)

    comps_new = eltype(comps)[]
    for (i, comp) in enumerate(comps)
        zz_i = view(zz_pre, i, :)
        comp_new = emstep(comp, samples, weights, zz_i; kwargs...)
        push!(comps_new, comp_new)
    end

    p = get_zz(comps_new, samples) * weights
    p ./= sum(p)

    constructor(comps_new, p)
end

function fit_em(mixture::MT, samples::AbstractVector, weights::AbstractWeights; tol=1e-4, maxiter::Int=1000, kwargs...) where {MT <: AbstractMixtureModel}
    logl_hist = [ sum(i -> logpdf(mixture, samples[i]) * weights[i], 1:length(samples)) ]

    i = 1
    while true
        mixture = emstep(mixture, samples, weights; kwargs...)
        logl = sum(i -> logpdf(mixture, samples[i]) * weights[i], 1:length(samples))

        push!(logl_hist, logl)
        if abs(logl - logl_hist[end-1]) < tol
            break
        end

        i += 1
        if i > maxiter
            @warn "MixtureModels.em  failed to achieve convergence after $maxiter iterations"
            break
        end
    end

    mixture, logl_hist
end

fit_em(mixture::AbstractMixtureModel, samples::AbstractVector; kwargs...) = fit_em(mixture, samples, uweights(length(samples)); kwargs...)
emstep(mixture::AbstractMixtureModel, samples::AbstractVector; kwargs...) = emstep(mixture, samples, uweights(length(samples)); kwargs...)


text:
if haskey(ENV, "CI") && ENV["CI"] == "true"
    ENV["PYTHON"] = ""
    using Pkg
    Pkg.build("PyCall")
end

using PyPlot
using Glob
using Oceananigans

using ConvergenceTests
using ConvergenceTests.ForcedFlowFreeSlip: u

arch = CUDA.has_cuda() ? GPU() : CPU()

defaultcolors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
removespine(side) = gca().spines[side].set_visible(false)
removespines(sides...) = [removespine(side) for side in sides]

xy_filenames = glob("forced_free_slip_xy*.jld2", joinpath(@__DIR__, "data"))
xz_filenames = glob("forced_free_slip_xz*.jld2", joinpath(@__DIR__, "data"))

xy_errors = ConvergenceTests.compute_errors((x, y, z, t) -> u(x, y, t), xy_filenames...)
xz_errors = ConvergenceTests.compute_errors((x, y, z, t) -> u(x, z, t), xz_filenames...)

sizes = ConvergenceTests.extract_sizes(xy_filenames...)

Nx = map(sz -> sz[1], sizes)

names = (L"(x, y)", L"(x, z)")
errors = (xy_errors, xz_errors)

close("all")
fig, ax = subplots()

for i = 1:length(errors)

    error = errors[i]
    name = names[i]

    L₁ = map(err -> err.L₁, error)
    L∞ = map(err -> err.L∞, error)

    loglog(Nx, L₁, color=defaultcolors[i], alpha=0.6, mfc="None",
           linestyle="None", marker="o", label="\$L_1\$-norm, $name")
    loglog(Nx, L∞, color=defaultcolors[i], alpha=0.6, mfc="None",
           linestyle="None", marker="^", label="\$L_\\infty\$-norm, $name")
end

L₁ = map(err -> err.L₁, errors[1])
loglog(Nx, L₁[end] * (Nx[end] ./ Nx).^2, "k-", linewidth=1, alpha=0.6, label=L"\sim N_x^{-2}")

legend()
xscale("log", base=2)
yscale("log", base=10)
xlabel(L"N_x")
ylabel("Norms of the absolute error, \$ | u_{\\mathrm{sim}} - u_{\\mathrm{exact}} | \$")
title("Convergence for forced free slip")
removespines("top", "right")

filename = "forced_free_slip_convergence_$(typeof(arch)).png"
filepath = joinpath(@__DIR__, "figs", filename)
mkpath(dirname(filepath))
savefig(filepath, dpi=480)

p = sortperm(Nx)
for (name, error) in zip(names, errors)
    L₁ = map(e -> e.L₁, error)
    L∞ = map(e -> e.L∞, error)
    name = "Forced flow free slip " * strip(name.s, '$')
    test_rate_of_convergence(L₁[p], Nx[p], expected=-2.0, atol=0.001, name=name * " L₁")
    test_rate_of_convergence(L∞[p], Nx[p], expected=-2.0, atol=0.005, name=name * " L∞")
end


text:
# Highlight the nested node.

struct Esc
    val::String
end

Esc(name::Symbol) = Esc(Base.text_colors[name])

Base.show(io::IO, esc::Esc) =
    if get(io, :color, false)
        print(io, esc.val)
    end

struct EscWrapper
    content::Any
    color::Symbol
    restore::Vector{Symbol}
end

function PrettyPrinting.tile(w::EscWrapper)
    lt = literal(Esc(w.color), 0) * tile_expr(w.content)
    for c in w.restore
        lt = lt * literal(Esc(c), 0)
    end
    lt
end

mutable struct HighlightNode <: AbstractSQLNode
    over::Union{SQLNode, Nothing}
    color::Symbol

    HighlightNode(;
           over = nothing,
           color::Union{Symbol, AbstractString}) =
        new(over, Symbol(color))
end

HighlightNode(color; over = nothing) =
    HighlightNode(over = over, color = color)

"""
    Highlight(; over = nothing; color)
    Highlight(color; over = nothing)

Highlight `over` with the given `color`.

Available colors can be found in `Base.text_colors`.

# Examples

```jldoctest
julia> q = Get.person_id |> Highlight(:bold);
```
"""
Highlight(args...; kws...) =
    HighlightNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(Highlight), pats::Vector{Any}) =
    dissect(scr, HighlightNode, pats)

function PrettyPrinting.quoteof(n::HighlightNode, ctx::QuoteContext)
    if ctx.limit
        ex = Expr(:call, nameof(Highlight), quoteof(n.color))
        if n.over !== nothing
            ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
        end
        return ex
    end
    push!(ctx.colors, n.color)
    ex = quoteof(n.over, ctx)
    pop!(ctx.colors)
    EscWrapper(ex, n.color, copy(ctx.colors))
end

label(n::HighlightNode) =
    label(n.over)

rebase(n::HighlightNode, n′) =
    HighlightNode(over = rebase(n.over, n′), color = n.color)


text:
"""
"""
struct AffineMap{D,T,L} <:Field
  jacobian::TensorValue{D,D,T,L}
  origin::Point{D,T}
  function AffineMap(jacobian::TensorValue{D,D,T,L}, origin::Point{D,T}) where {D,T,L}
    new{D,T,L}(jacobian,origin)
  end
end

function field_cache(f::AffineMap,x)
  y = copy(x)
  CachedArray(y)
end

@inline function evaluate_field!(cache,f::AffineMap,x)
  setsize!(cache,size(x))
  y = cache.array
  @inbounds for i in eachindex(x)
    xi = x[i]
    yi = _apply_affine_map(f,xi)
    y[i] = yi
  end
  y
end

function _apply_affine_map(h,x)
  t = h.origin
  s = h.jacobian
  (s⋅x)+t
end

struct AffineMapGrad{D,T,L} <: Field
  jacobian::TensorValue{D,D,T,L}
end

function field_gradient(h::AffineMap)
  AffineMapGrad(h.jacobian)
end

function field_cache(f::AffineMapGrad,x)
  xi = testitem(x)
  T = typeof(xi)
  G = gradient_type(T,xi)
  j = similar(x,G)
  CachedArray(j)
end

function evaluate_field!(cache,f::AffineMapGrad,x)
  setsize!(cache,size(x))
  y = cache.array
  G = eltype(y)
  yi = f.jacobian
  @inbounds for i in eachindex(y)
    y[i] = yi
  end
  y
end


text:
include("../../../get_home_dir.jl")

using DataFrames
using CSV
#using Plots
using Random
using StatsBase

include(string(PROJECT_HOME,"/assim_code/mironov.jl"));
include(string(PROJECT_HOME,"/assim_code/tau_omega_funs.jl"));
include(string(PROJECT_HOME,"/assim_code/time_averaging.jl"))

df_raw = CSV.read(string(PROJECT_HOME,"/data/moflux_land_data_skipyear_hourly2.csv"), DataFrame);
df_raw[!,"RAIN"] *= 2;

N = 24*365*12
istart = 24*365*0 + 1;
soil0 = 0.30;

include(string(PROJECT_HOME,"/simulation_code/sim_vary_new_stomata_med.jl"));

deltaT = FT(60*60);
alpha = FT(50);
nsoil = FT(1.2);

function run_sim_2(vcmax_par::FT, k_frac::FT, k_plant::FT,  z_soil::FT, weibB::FT, vol_factor::FT, g1::FT)
        return run_sim_vary(vcmax_par, k_frac, weibB, FT(2), k_plant, FT(1e-4), z_soil, istart, N, soil0, vol_factor, 1e-5, 3, df_raw, g1, deltaT, alpha, nsoil);
end

sim_res1 = run_sim_2(FT(31),FT(0.5), FT(16), FT(2000),FT(5),FT(0.33),FT(506));

noise_std = 1.3;
noise_var = noise_std^2

tsoil = sim_res1[1].T_SOIL .+ 273.15;
tair = sim_res1[1].T_AIR .+ 273.15;
tcan = (sim_res1[1].LW_OUT / 0.97 / K_STEFAN ) .^ 0.25;
laiM = sim_res1[1].LAI_modis;
canpot_true = mean(sim_res1[3],dims=2);

function write_post(mydir)

par_post = convert(Array,CSV.read(string(mydir,"post_par.csv"),DataFrame))[5501:90:end,1:7];
outdir = mydir

nday = convert(Int, N/24);
npost = size(par_post)[1];

post_RZ = zeros(nday,npost+1);
#post_Surf = zeros(nday,npost+1);
post_ET = zeros(nday*24,npost+1);
post_LWP = zeros(nday*24,npost+1);
#post_branch = zeros(nday*24,npost+1);
post_trunk = zeros(nday*24,npost+1);
post_gsw = zeros(nday*24,npost+1);
post_TWS = zeros(nday,npost+1);
post_TWP = zeros(nday,npost+1);

for j in 1:npost
	println(j)
	a = par_post[j,:]
try
	sim_res0 = run_sim_2(convert(Array{FT}, exp.(a))...);

	post_RZ[:,j] = get_daily(sim_res0[2][:,end],24);
        post_TWS[:,j] = get_daily(sim_res0[1].ColumnSMC,24);
        post_TWP[:,j] = get_daily(sim_res0[1].ColumnSWP,24);
	post_ET[:,j] = sim_res0[1].ETmod;
	post_gsw[:,j] = sim_res0[1].glw;
	post_LWP[:,j] = mean(sim_res0[3],dims=2)[:,1];
	#post_branch[:,j] = mean(sim_res0[4],dims=2)[:,1];
	post_trunk[:,j] = sim_res0[5];
catch err
	println("Error", a)
end
end

	#place the results of the original true run in the last column
	post_RZ[:,end] = get_daily(sim_res1[2][:,end],24);
        post_TWS[:,end] = get_daily(sim_res1[1].ColumnSMC,24);
        post_TWP[:,end] = get_daily(sim_res1[1].ColumnSWP,24);
	post_ET[:,end] = sim_res1[1].ETmod;
	post_gsw[:,end] = sim_res1[1].glw;
	post_LWP[:,end] = mean(sim_res1[3],dims=2)[:,1];
	#post_branch[:,end] = mean(sim_res1[4],dims=2)[:,1];
	post_trunk[:,end] = sim_res1[5];

#outdir = "y2/"

#post_RZ, post_ET, post_LWP, post_branch, post_trunk;
CSV.write(string(outdir,"postRZ.csv"), DataFrame(post_RZ) );
CSV.write(string(outdir,"postSWS.csv"), DataFrame(post_TWS) );
CSV.write(string(outdir,"postSWP.csv"), DataFrame(post_TWP) );
CSV.write(string(outdir,"postET.csv"), DataFrame(post_ET));
CSV.write(string(outdir,"postLeaf.csv"), DataFrame(post_LWP));
#CSV.write(string(outdir,"postBranch.csv"), DataFrame(post_branch));
CSV.write(string(outdir,"postTrunk.csv"), DataFrame(post_trunk));
CSV.write(string(outdir,"postGSW.csv"), DataFrame(post_gsw));
end

write_post(ARGS[1])


text:
# [test/odehybrid_tests.jl]
using Test

@testset "ODE Hybrid Tests" begin

    using JLD2
    import ODEHybrid.output_discrete_states, ODEHybrid.ODE_SET, ODEHybrid.run_de, ODEHybrid.run_ode,
           ODEHybrid.rk4, ODEHybrid.rkfixed, ODEHybrid.rkadapt, ODEHybrid.odehybrid, ODEHybrid.log_frame,
           ODEHybrid.odehybridcore, ODEHybrid.odehybridfull, ODEHybrid.TimeSeriesLogger, ODEHybrid.add!, ODEHybrid.vec_to_mat

    plots_on = false;

    @testset "Output Discrete States" begin
        yd0 = [1.0, 5.0]
        states = [1.0 2.0 3.0 4.0;
                1.5 2.5 3.5 4.5;
                2.0 3.0 4.0 5.0]
        # mat"""
        #     states = {1.0, 2.0, 3.0, 4.0; 1.5, 2.5, 3.5, 4.5; 2.0, 3.0, 4.0, 5.0};
        #     yd0 = {1.0, 5.0};
        #     $dis_m = output_discrete_states(states, yd0)
        # """
        dis_j = output_discrete_states(states, yd0)
        @test dis_j == [ [1, 1.5, 2], [2, 2.5, 3]]

        yd0 = [1.0]
        states = [1.0 2.0 3.0 4.0;
                1.5 2.5 3.5 4.5;
                2.0 3.0 4.0 5.0]
        # mat"""
        #     states = {1.0, 2.0, 3.0, 4.0; 1.5, 2.5, 3.5, 4.5; 2.0, 3.0, 4.0, 5.0};
        #     yd0 = {1.0};
        #     $dis_m = output_discrete_states(states, yd0)
        # """
        dis_j = output_discrete_states(states, yd0)
        @test dis_j == [[1, 1.5, 2]]

    end;

    @testset "Run DE" begin

        # Test 1
        de = (t, x, u) -> (x, -[8 4] * float.(x));
        t = 0.0
        ycv = [97.0; 98.0]  # Float = [97; 98]
        yd = [0]
        yc0 = ['a'; 'b']

        yc_jul, yd_jul = run_de(de, t, ycv, yd, yc0)

        # mat"""
        #     [$yc_mat, $yd_mat] = run_de(@(t, x, u) deal(x, -[8, 4] * x),  $t, [97; 98], {[0]},{['a'; 'b']});
        # """

        @test yc_jul == [97; 98.0]
        @test yd_jul == [-1168.0]

        # Test 2
        de = (t, x, u) -> (x, -[8 4 -2 -2] * float.(x));
        t = 0.0
        ycv = [97.0; 98.0; 99.0; 100.0]  # Float = [97; 98]
        yd = [0]
        yc0 = ['a'; 'b'; 'c'; 'd']

        yc_jul, yd_jul = run_de(de, t, ycv, yd, yc0)

        # mat"""
        #     [$yc_mat, $yd_mat] = run_de(@(t, x, u) deal(x, -[8, 4, -2, -2] * x),  $t, [97; 98; 99; 100], {[0]},{['a'; 'b'; 'c'; 'd']});
        # """

        @test yc_jul == [97.0; 98; 99; 100]
        @test yd_jul == [-770]
    end;

    @testset "Run ODE" begin
        ode = (t, x, u) -> [0.0 1; 2 0] * float.(x) + [0; 1.0] .* float.(u);  # Differential Equation
        t   = 0.5
        yc0  = ['a'; 'b'];                               # Initial Continuous state
        yd0  = [0.0];                                    # initial discrete state
        ycv  = [97.0; 98.0]

        dy_jul = run_ode(ode, t, ycv, yd0, yc0)
        # mat"""
        #     $dy_mat = run_ode(@(t, x, u) [[0 1; 2 0] * x + [0; 1] * u], $t, $ycv, {$yd0}, {['a'; 'b']})
        # """

        @test dy_jul ==  [98; 194]
    end;

    @testset "Log Frame" begin

        function ode(t, x, u, tsl::TimeSeriesLogger)  # Differential Equation
            add!(tsl, "A", t, x, true, "group1")
            return [0 1; 2 0] * float.(x) + [0; 1] * float.(u);
        end

        # Init
        t = 0.5
        y = [7.5; 8.0]  # Data points × t
        flag = "init"
        yd  = 0.1
        tsl = TimeSeriesLogger()
        s_init_jul = log_frame(t, y, flag, ode, yd, tsl)

        # Empty
        flag = "";
        t = 0.75;
        y = [8.0; 9.0];
        s_empty1_jul = log_frame(t, y, flag, ode, yd, tsl)

        flag = "";
        t = 1.0;
        y = [9.0; 9.25];
        s_empty2_jul = log_frame(t, y, flag, ode, yd, tsl)

        # Done
        flag = "done"
        t = 1.25;
        y = [9.25; 9.5];
        s_done_jul   = log_frame(t, y, flag, ode, yd, tsl)

        tsl_jul = tsl;

        # mat"""
        #     tsl = TimeSeriesLogger();
        #     f = '';

        #     % INIT
        #     $s_init_mat = log_frame(0.5, [7.5; 8.0], "init", @test_ode, {[0.1]}, tsl, f);

        #     % EMPTY 1
        #     $s_empty1_mat = log_frame(0.75, [8.0; 9.0], '',  @test_ode, {[0.1]}, tsl, f);

        #     % EMPTY 2
        #     $s_empty2_mat = log_frame(1.0, [9.0; 9.25], '',  @test_ode, {[0.1]}, tsl, f);

        #     % DONE
        #     $s_done_mat   = log_frame(1.0, [9.25; 9.5], "done",  @test_ode, {[0.1]}, tsl, f);

        #     % MATLAB has a hard time returning the TSL struct ...

        #     $names_mat = tsl.names;
        #     $data_mat  = tsl.data;
        #     $show_mat  = tsl.show;
        #     $sizes_mat = tsl.sizes;
        #     $count_mat = tsl.count;
        # """

        @test tsl_jul.names == ["A"]
        @test tsl_jul.times[1] == [0.75; 1.0]
        @test vec_to_mat(tsl_jul.data[1])  == [8 9; 9 9.25]
        @test tsl_jul.show  == [true]
        @test tsl_jul.sizes == [2]
        @test tsl_jul.count == [2]

    end;

    # Doesn't appear to be called...?
    @testset "Output Fcn" begin

        # function ode(t, x, u, tsl::TimeSeriesLogger)  # Differential Equation
        #     add!(tsl, "A", t, x, true, "group1")
        #     return [0 1; 2 0] * float.(x) + [0; 1] * float.(u);
        # end

        # # Try LOGFRAME and something else
        # f = ode
        # t = 0.0
        # ycv = [97.0; 98.0]  # Float = [97; 98]
        # yd = [0]
        # yc0 = [0.0; 0.0]    # ['a'; 'b']
        # flag = "init"

        # status_jul = output_fcn(f, t, ycv, yd, flag, yc0)

        @test 1 == 1

    end;

    # Doesn't appear to be called...?
    @testset "Event Fcn" begin
        @test 1 == 1
    end;

    @testset "ODE Hybrid Core" begin
        ode = (t, x, u) -> [0 1; 2 0] * x + [0; 1] * u;
        de  = (t, x, u) -> (x, -[8 4] * x);             # Discrete update equation
        dt  = 0.1;                                      # Discrete eq. time step
        ts  = [0.0 5.0];                                # Simulate from 0 to 5s
        yc0  = [1.0; 0.0];                              # Initial Continuous state
        yd0  = [0.0];

        # No events, ignore those (empty) outputs
        t_jul, yc_jul, td_jul, yd_jul, _, _, _, _  = odehybridcore(rkadapt, ode, de, dt, ts, yc0, yd0);

        # mat"""
        #     $results_mat = odehybridcore(@rkadapt, @(t, x, u) [0 1; 2 0] * x + [0; 1] * u, @(t, x, u) deal(x, -[8 4] * x), 0.1, [0.0 5.0], [1.0; 0.0], 0.0);
        # """

        # yd_mat = results_mat["yd"]
        # td_mat = results_mat["td"]
        # t_mat  = results_mat["t"]
        # yc_mat = results_mat["yc"]

        @load "solutions/odehybrid_hybrid_core.jld2" yd_mat td_mat t_mat yc_mat

        @test yd_jul ≈ yd_mat;
        @test td_jul ≈ td_mat;
        @test t_jul  ≈ t_mat;
        @test yc_jul ≈ yc_mat;

        if (plots_on)
            a = plot(td_jul, yd_jul, label = "Julia");
            display(plot!(td_mat, yd_mat, label = "Matlab"))

            b = plot(t_jul, yc_jul, label = "Julia")
            display(plot!(t_mat, yc_mat, label = "Matlab"))
        end

    end;

    @testset "ODE Hybrid Full" begin
        ode = (t, x, u) -> [0 1; 2 0] * x + [0; 1] * u;
        de  = (t, x, u) -> (x, -[8 4] * x);             # Discrete update equation
        dt  = 0.1;                                      # Discrete eq. time step
        ts  = [0.0 5.0];                                # Simulate from 0 to 5s

        struct TESTT
            a
            b
            c
        end

        ycv = TESTT(1.0, 2.0, 3.0)                    # Initial Continuous state
        yd0 = [0.0];
        yc0 = [1.0; 0.0]

        res_jul = odehybridfull(rkadapt, ode, de, dt, ts, yc0, yd0)

        # mat"""
        # % Opt 2: 4 outputs
        #     [$t1, $t2, $t3, $t4] = odehybridfull(@rkadapt, @(t, x, u) [0 1; 2 0] * x + [0; 1] * u, @(t, x, u) deal(x, -[8 4] * x), 0.1, [0.0 5.0], [1.0; 0.0], 0.0)
        # """

        @load "solutions/odehybrid_hybrid_full.jld2" t1 t2 t3 t4

        @test res_jul[1] ≈ t1;
        @test res_jul[2] ≈ t2;
        @test res_jul[3] ≈ t3;
        @test res_jul[4] ≈ t4;

        # mat"""
        #     % Opt 3: 2 outputs
        #     [$t1, $t2] = odehybridfull(@rkadapt, @(t, x, u) [0 1; 2 0] * x + [0; 1] * u, @(t, x, u) deal(x, -[8 4] * x), 0.1, [0.0 5.0], [1.0; 0.0], 0.0)
        # """

        @load "solutions/odehybrid_hybrid_full2.jld2" t1 t2

        @test res_jul[1] ≈ t1;
        @test res_jul[2] ≈ t2;
    end;

    # Examples from AnUncommonLab.com
    @testset "'Short and Direct' Example" begin

        # mat"""
        #     [$tm, $xm, $tum, $um] = odehybrid(@rkadapt, @(t, x, u) [0 1; 2 0] * x + [0; 1] * u, @(t, x, u) deal(x, -[8 4] * x), 0.1, [0.0 5.0], [1.0; 0.0], 0.0)
        # """

        @load "solutions/odehybrid_short_direct.jld2" tm xm tum um

        ode = (t, x, u) -> [0 1; 2 0] * x + [0; 1] * u; # Differential Equation
        de  = (t, x, u) -> (x, -[8 4] * x);             # Discrete update equation
        dt  = 0.1;                                      # Discrete eq. time step
        ts  = [0.0 5.0];                                # Simulate from 0 to 5s
        x0  = [1.0; 0.0];                               # Initial Continuous state
        u0  = [0.0];                                    # initial discrete state

        tj, xj, tuj, uj = odehybrid(rkadapt, ode, de, dt, ts, x0, u0); # Simulate!

        if (plots_on)
            scatter(tm, xm, xlabel = "Time", label = ["x₁" "x₂"])
            display(scatter!(tum, um, label = "u"))

            scatter(tj, xj, xlabel = "Time", label = ["x₁" "x₂"])
            display(scatter!(tuj, uj, label = "u"))
        end

        mutable struct TEST2
            a
            b
        end

        # Currently likes to make things arrays, even structs, so we index in
        ode = (t, x, u) -> [0 1; 2 0] * [x.a; x.b] + [0; 1] * u
        de  = (t, x, u) -> (x, -[8 4] * [x.a; x.b])
        x0 = TEST2(1.0, 0.0)
        u0  = [0.0];                                      # initial discrete state

        tj, xj, tuj, uj = odehybrid(rkadapt, ode, de, dt, ts, x0, u0); # Simulate!

        @test tum ≈ tuj
        @test um  ≈  uj

    end;

    # NOTE that de MUST not be an ARRAY in this format (make it a Tuple)
    @testset "'Multiple States' Example" begin

        # mat"""
        #     %                   Continuous System    +  feedback control   +    disturbance
        #     ode = @(t, x, u, i)   [0 1; 2 0] * x     +    [0; 1] * u       +      [0; 1];    % Continuous system

        #     %                No change to cont state,        Update input        Update Integrator
        #     de  = @(t, x, u, i)    deal(x,                 -[8 4 1] * [x; i],      i + 0.5 * x(1));

        #     dt  = 0.1;                                  % Discrete eq. time step
        #     ts  = [0 5];                                % From 0 to 5s
        #     x0  = [1; 0];                               % Initial continuous state

        #     d0  = {0, 0};                               % Initial discrete states

        #     [$t_mat, $x_mat, $td_mat, $u_mat, $i_mat] = odehybrid(@rkadapt, ode, de, dt, ts, x0, d0); % Simulate!

        #     if ($plots_on)
        #         plot($t_mat, $x_mat, $td_mat, [$u_mat, $i_mat], '.'); xlabel('Time')
        #         legend('x_1', 'x_2', 'u', 'x_1(t)/2 dt')
        #     end
        # """

        @load "solutions/odehybrid_multi_states.jld2" t_mat x_mat td_mat u_mat i_mat

        # ODE includes "integral" term 'i'
        ode = (t, x, u, i) -> [0 1; 2 0] * x +     # Continuous System
                            [0; 1] * u +         #    with feedback control
                            [0; 1];              #    and with a distrubance

        de  = (t, x, u, i) -> (x,                  # No change to continuous state
                            ((-[8 4 1] * [x; i])[1],  # Update input
                            i + 0.5 * x[1]));    # Update integrator

        # Settings
        dt = 0.1;       # Discrete time step
        ts = [0.0 5.0];     # Simulate from 0 to 5s
        x0 = [1.0; 0.0];    # Initial continuous state

        d0 = (0.0, 0.0)     # Initial discrete states

        t_jul, x_jul, td_jul, out_jul = odehybrid(rkadapt, ode, de, dt, ts, x0, d0);  # simulate!
        u_jul = out_jul[:, 1];
        i_jul = out_jul[:, 2];

        @test t_jul ≈ t_mat
        @test x_jul ≈ x_mat
        @test td_jul ≈ td_mat
        @test u_jul ≈ u_mat
        @test i_jul ≈ i_mat

        if (plots_on)
            plot(t_jul, x_jul, label = ["x₁" "x₂"])
            scatter!(td_jul, u_jul, label = "u")
            display(scatter!(td_jul, i_jul, label = "i"))
        end
    end;

    # Same as previous but formated differently
    @testset "'Multiple Continuous & Discrete States' Example" begin

        # mat"""
        #     ode = @(t, p, v, u, i) deal(v, 2 * p + u + 1);

        #     de  = @(t, p, v, u, i) deal(p, v, -8*p - 4*v - i, i + 0.5 * p);

        #     dt  = 0.1;                                      % Discrete eq. time step
        #     ts  = [0 5];                                    % From 0 to 5s
        #     x0  = {1; 0};                                   % Initial continuous states
        #     d0  = {0, 0};                                   % Initial discrete states
        #     [$t_mat, $p_mat, $v_mat, $td_mat, $u_mat, $i_mat] = odehybrid(@rkadapt, ode, de, dt, ts, x0, d0);

        #     if ($plots_on)
        #         plot(t, [p, v], td, [u, i], '.'); xlabel('Time');
        #         legend('p', 'v', 'u', ' p(t)/2 dt');
        #     end;
        # """

        @load "solutions/odehybrid_multi_cont_disc_states.jld2" t_mat p_mat v_mat td_mat u_mat i_mat

        ode = (t, p, v, u, i) -> (v, 2 * p + u + 1);
        de  = (t, p, v, u, i) -> ((p, v), (-8 * p - 4 * v - i, i + 0.5 * p));
        dt  = 0.1;
        ts  = [0.0 5.0];
        x0  = (1, 0);
        d0  = (0, 0);

        # No events so we can ignore those outputs
        t_jul, yc_jul, td_jul, yd_jul, _, _, _, _ = odehybrid(rkadapt, ode, de, dt, ts, x0, d0);
        p_jul = yc_jul[:, 1];
        v_jul = yc_jul[:, 2];
        u_jul = yd_jul[:, 1];
        i_jul = yd_jul[:, 2];

        @test t_mat  ≈ t_jul
        @test p_mat  ≈ p_jul
        @test v_mat  ≈ v_jul
        @test td_mat ≈ td_jul
        @test u_mat  ≈ u_jul
        @test i_mat  ≈ i_jul

    end;

    @testset "'Structs for States' Example" begin

        # # MATLAB functions found in 'original_matlab' folder
        # mat"""
        #     [t, sig, td, ctrl] = example_odehybrid_structs();

        #     p = zeros(length(t), 1);
        #     v = zeros(length(t), 1);
        #     u = zeros(length(td), 1);
        #     i = zeros(length(td), 1);

        #     for j = 1:length(t)
        #         p(j) = sig(j).p;
        #         v(j) = sig(j).v;
        #     end;

        #     for j = 1:length(td)
        #         u(j) = ctrl(j).u;
        #         i(j) = ctrl(j).i;
        #     end;

        #     $t_mat = t;
        #     $p_mat = p;
        #     $v_mat = v;

        #     $td_mat = td;
        #     $u_mat  = u;
        #     $i_mat  = i;
        # """

        @load "solutions/odehybrid_struct_states.jld2" t_mat p_mat v_mat td_mat u_mat i_mat

        mutable struct TEST_STRUCT3
            p
            v
        end;

        mutable struct TEST_STRUCT4
            u
            i
        end;

        function ode(_, sig, ctrl)
            p = sig.v;
            v = 2 * sig.p + ctrl.u + 1;
            return TEST_STRUCT3(p, v)
        end;

        function de(_, sig, ctrl)
            u = -8 * sig.p - 4 * sig.v - ctrl.i;
            i = ctrl.i + 0.5 * sig.p;
            return sig, TEST_STRUCT4(u, i)
        end;

        x0 = TEST_STRUCT3(1.0, 0.0); # Initial continuous states
        d0 = TEST_STRUCT4(0.0, 0.0); # Initial discrete states

        dt = 0.1;                    # Discrete eq. time step
        ts = [0 5];                  # Simulation time

        # Simulate.
        t_jul, sig, td_jul, ctrl, _, _, _, _ = odehybrid(rkadapt, ode, de, dt, ts, x0, d0);

        # Not *quite* sure why one returns as a struct and the other an array...
        p_jul = sig[:, 1];
        v_jul = sig[:, 2];
        u_jul = [ctrl[j].u for j = 1:length(ctrl)];
        i_jul = [ctrl[j].i for j = 1:length(ctrl)];

        # Plot.
        if (plots_on)
            plot(t, [sig.p], t, [sig.v], td, [ctrl.u], '.', td, [ctrl.i], '.');
            xlabel("Time");
            legend('p', 'v', 'u', "Int p(t)/2 dt");
        end

        @test t_jul ≈ t_mat
        @test p_jul ≈ p_mat
        @test v_jul ≈ v_mat
        @test u_jul ≈ u_mat
        @test i_jul ≈ i_mat
    end;

    # Generates plots for comparison
    @testset "'With Logging' Example" begin

        # mat"""
        #     [t, sig, td, ctrl] = example_odehybrid_logging();

        #     p = zeros(length(t), 1);
        #     v = zeros(length(t), 1);
        #     for i = 1:length(t)
        #         p(i) = sig(i).p;
        #         v(i) = sig(i).v;
        #     end;

        #     u = zeros(length(td), 1);
        #     i = zeros(length(td), 1);
        #     for j = 1:length(td)
        #         u(j) = ctrl(j).u;
        #         i(j) = ctrl(j).i;
        #     end;

        #     $t_mat = t;
        #     $p_mat = p;
        #     $v_mat = v;

        #     $td_mat = td;
        #     $u_mat  = u;
        #     $i_mat  = i;
        # """

        @load "solutions/odehybrid_with_logging.jld2" t_mat p_mat v_mat td_mat u_mat i_mat

        mutable struct TEST_STRUCT5
            p
            v
        end;

        mutable struct TEST_STRUCT6
            u
            i
        end;

        function ode(t, sig, ctrl, log = [])
            p = sig.v;
            v = 2 * sig.p + ctrl.u + 1;

            # Log the acceleration (not always passed in!)
            if (log != [])
                add!(log, "acceleration", t, v);
            end

            return TEST_STRUCT5(p, v);
        end;

        function de(t, sig, ctrl, log)
            u = -8 * sig.p - 4 * sig.v - ctrl.i;
            i = ctrl.i + 0.5 * sig.p;

            # Log velocity (always passed in)
            if (log != [])
                add!(log, "Sampled v", t, sig.v)
            end

            return sig, TEST_STRUCT6(u, i)
        end;

        dt  = 0.1;                    # Discrete eq. time step
        ts  = [0.0 5];                  # Simulation time
        x0  = TEST_STRUCT5(1.0, 0.0); # Initial continuous states
        d0  = TEST_STRUCT6(0.0, 0.0); # Initial discrete states
        log = TimeSeriesLogger();     # Create the logger.

        # Simulate.
        t_jul, sig, td_jul, ctrl, _, _, _, _ = odehybrid(rkadapt, ode, de, dt, ts, x0, d0, [], log);

        p_jul = sig[:, 1];
        v_jul = sig[:, 2];
        u_jul = [ctrl[j].u for j = 1:length(ctrl)];
        i_jul = [ctrl[j].i for j = 1:length(ctrl)];

        @test t_jul ≈ t_mat
        @test p_jul ≈ p_mat
        @test v_jul ≈ v_mat
        @test u_jul ≈ u_mat
        @test i_jul ≈ i_mat

        # Plot
        if (plots_on)
            plot(t_jul, p_jul, xlabel = "Time", label = "p", title = "With Logging (Julia)")
            plot!(t_jul, v_jul, label = "v")
            scatter!(td_jul, u_jul, label = "u")
            display(scatter!(td_jul, i_jul, label = "i"))

            tsl_plot(log)
        end

    end;

    @testset "'Multiple Discrete Updates' Example" begin

        # mat"""
        #     ode = @(t, x, ~, ~) [0 1; -1 0] * x;         % Differential equation

        #     % TWO discrete update equations
        #     de  = {@(t, x, p, v) deal(x, x(1), v);  @(t, x, p, v) deal(x, p, x(2))};

        #     dt  = [0.1, 0.15];                           % Discrete eq. time steps
        #     ts  = [0 2*pi];                              % From 0 to 6.3s
        #     x0  = [1; 0];                                % Initial continuous state
        #     y0  = {0, 0};                                % Initial discrete state

        #     [$t_mat, $sig_mat, $ty_mat, $y1_mat, $y2_mat] = odehybrid(@rkfixed, ode, de, dt, ts, x0, y0);
        # """

        @load "solutions/odehybrid_mult_disc_updates.jld2" t_mat sig_mat ty_mat y1_mat y2_mat

        ode = (t, x, _, _) -> [0.0 1.0; -1.0 0.0] * x;     # Diff Eq
        de  = [(t, x, p, v) -> (x, (x[1], v));               # Discrete Eq #1
            (t, x, p, v) -> (x, (p, x[2]))];             # Discrete Eq #2

        dt  = [0.1, 0.15];                                 # Discrete Eq. Time steps
        ts  = [0.0  2*pi];                                 # From 0.0 to 6.4s
        x0  = [1.0; 0.0];                                      # Initial Continuous State
        y0  = (0.0, 0.0);                                  # Initial Discrete State

        t_jul, sig_jul, ty_jul, ctrl, _, _, _, _ = odehybrid(rkfixed, ode, de, dt, ts, x0, y0);
        y1_jul = ctrl[:, 1];
        y2_jul = ctrl[:, 2];

        @test t_jul   ≈ t_mat
        @test sig_jul ≈ sig_mat
        @test ty_jul  ≈ ty_mat
        @test y1_jul  ≈ y1_mat
        @test y2_jul  ≈ y2_mat
    end;

    @testset "'Output Function' Example" begin

        # mat"""
        #     ode = @(t, x, ~, ~) [0 1; -1 0] * x;         % Differential equation

        #     % TWO discrete update equations
        #     de  = {@(t, x, p, v) deal(x, x(1), v); @(t, x, p, v) deal(x, p, x(2))};

        #     dt  = [0.1, 0.15];                           % Discrete eq. time steps
        #     ts  = [0 2*pi];                              % From 0 to 6.3s
        #     x0  = [1; 0];                                % Initial continuous state
        #     y0  = {0, 0};                                % Initial discrete state

        #     options = odeset('OutputFcn', @example_odehybrid_outputfcn);

        #     % DOESN'T return correctly (y1 and y2 are scalars?)
        #     [$t_mat, $sig_mat, $ty_mat, $y1_mat, $y2_mat] = odehybrid(@rkfixed, ode, de, dt, ts, x0, y0, options);
        # """

        @load "solutions/odehybrid_output_func.jld2" t_mat sig_mat ty_mat y1_mat y2_mat

        function custom_output_fcn(t, x, ignore, flag)
            status = 0;

            if flag == "init"
                ts = t[1]
                tf = t[end]
                println("Simulating from $ts to $tf.")

            elseif flag == "done"
                println("Finished the simulation.")

            else
                status = x[1] < 0;  # End simulation when x[1] < 0

            end

            return status;
        end

        options = ODE_SET();
        options.OutputFcn = custom_output_fcn;

        # Add to previous simulation

        ode =  (t, x, _, _) -> [0.0 1.0; -1.0 0.0] * x;     # Diff Eq
        de  = [(t, x, p, v) -> (x, (x[1], v));               # Discrete Eq #1
            (t, x, p, v) -> (x, (p, x[2]))];             # Discrete Eq #2

        dt  = [0.1, 0.15];                                 # Discrete Eq. Time steps
        ts  = [0.0  2*pi];                               # From 0.0 to 6.4s
        x0  = [1.0; 0.0];                                  # Initial Continuous State
        y0  = (0.0, 0.0);                                  # Initial Discrete State

        t_jul, sig_jul, ty_jul, y1_jul, _, _, _, _ = odehybrid(rkfixed, ode, de, dt, ts, x0, y0, options);

        @test t_jul   ≈ t_mat
        @test sig_jul ≈ sig_mat
        @test ty_jul  ≈ ty_mat

        # @test y1_jul  ≈ y1_mat
        @test minimum(y1_jul[1:end-1]) > 0.0  # Verify that the simulation stops when x[1] < 0 (will have one step past 0)

    end;

    # # Events only work with MATLAB's ODE suite, which is not implemented here

    # @testset "'Events' Example" begin

    #     mat"
    #         ode = @(t, x, ~, ~) [0 1; -1 0] * x;         % Differential equation
    #         de  = {@(t, x, p, v) deal(x, x(1), v); ...   % Discrete update equation 1
    #             @(t, x, p, v) deal(x, p, x(2))};      % Discrete update equation 2
    #         dt  = [0.1, 0.15];                           % Discrete eq. time steps
    #         ts  = [0 2*pi];                              % From 0 to 6.3s
    #         x0  = [1; 0];                                % Initial continuous state
    #         y0  = {0, 0};                                % Initial discrete state

    #         options = odeset('Events', @example_odehybrid_eventfcn);

    #         [$t_mat, $sig_mat, $ty_mat, $y1_mat, $y2_mat] = odehybrid(@ode45, ode, de, dt, ts, x0, y0, options);
    #     "

    #     function custom_event_fcn(t, x, p, v)
    #         h = x[2];   # "Event" is triggered when velocity == 0
    #         t = true;   # Terminate on event (when h = 0)
    #         d = 1;      # Trigger when going from positive to negative

    #         return [h, t, d]
    #     end

    #     options = ODE_SET();
    #     options.Events = custom_event_fcn;

    #     # Add to previous simulation

    #     ode =  (t, x, _, _) -> [0.0 1.0; -1.0 0.0] * x;     # Diff Eq
    #     de  = [(t, x, p, v) -> (x, (x[1], v));               # Discrete Eq #1
    #         (t, x, p, v) -> (x, (p, x[2]))];             # Discrete Eq #2

    #     dt  = [0.1, 0.15];                                 # Discrete Eq. Time steps
    #     ts  = [0.0  2*pi];                               # From 0.0 to 6.4s
    #     x0  = [1.0; 0.0];                                  # Initial Continuous State
    #     y0  = (0.0, 0.0);                                  # Initial Discrete State

    #     t_jul, sig_jul, ty_jul, y1_jul, _, _, _, _ = odehybrid(rkfixed, ode, de, dt, ts, x0, y0, options);

    # end;

end;


text:
module adapters_postgresql_null_test

using Test
using Octo.Adapters.PostgreSQL

Repo.debug_sql()

Repo.connect(
    adapter = Octo.Adapters.PostgreSQL,
    dbname = "postgresqltest",
    user = "postgres",
)

struct Book
end
Schema.model(Book, table_name="book")

Repo.execute([DROP TABLE IF EXISTS Book])
Repo.execute(Raw("""
    CREATE TABLE book (
      id SERIAL PRIMARY KEY,
      title VARCHAR(255),
      publisher_description VARCHAR(255) NULL
    );
    """))

# https://discourse.julialang.org/t/how-to-pass-a-null-value-to-a-prepared-statement-with-libpq-jl/25068/2

Repo.insert!(Book, (title="Julia",))
Repo.insert!(Book, (title="O'Neil and the Martians", publisher_description=missing))

df = Repo.query(Book)
@test ismissing(df[1].publisher_description)
@test ismissing(df[2].publisher_description)

Repo.disconnect()

end # module adapters_postgresql_null_test


text:
# Author: Lukasz Bratos

include("diff.jl")

f(x) = abs(x)
g(x) = 1.0 / (1.0 + x^2)

for i in [5, 10, 15]
    Diff.rysujNnfx(f, -1.0, 1.0, i, "plot6a$i")
    Diff.rysujNnfx(g, -5.0, 5.0, i, "plot6b$i")
end


text:
using StatsBase
using RDatasets
using Distributions
using Gadfly
sleep = dataset("lme4","sleepstudy")
plot(x = sleep[:Reaction], Geom.histogram(bincount = 30), color = sleep[:Days])
plot(sleep, x = "Days", y = "Reaction", Geom.point)
plot(sleep, x = "Reaction", Geom.density, color = "Subject")

plot([sin, cos], 0, 25)

plot(dataset("datasets", "iris"),x="SepalLength", y="SepalWidth", Geom.point)

max_stars_repo_path: Chapter 3/8.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 3/8.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 3/8.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
##==============================================================================
## Abstract Types
##==============================================================================

abstract type InferenceVariable end

##==============================================================================
## VariableNodeData
##==============================================================================

"""
$(TYPEDEF)
Data container for solver-specific data.

  ---
Fields:
$(TYPEDFIELDS)
"""
mutable struct VariableNodeData{T<:InferenceVariable, P}
    val::Vector{P}
    bw::Matrix{Float64}
    BayesNetOutVertIDs::Vector{Symbol}
    dimIDs::Vector{Int} # Likely deprecate

    dims::Int
    eliminated::Bool
    BayesNetVertID::Symbol #  Union{Nothing, }
    separator::Vector{Symbol}

    variableType::T
    initialized::Bool
    """
    Replacing previous `inferdim::Float64`, new `.infoPerCoord::Vector{Float64}` will in
    future stores the amount information (per measurement dimension) captured in each
    coordinate dimension.
    """
    infoPerCoord::Vector{Float64}
    ismargin::Bool

    dontmargin::Bool
    solveInProgress::Int
    solvedCount::Int
    solveKey::Symbol

    events::Dict{Symbol,Threads.Condition}

    # VariableNodeData{T,P}() = new{T,P}()
    VariableNodeData{T,P}(w...)  where {T <:InferenceVariable, P} = new{T,P}(w...)
    VariableNodeData{T,P}(;solveKey::Symbol=:default ) where {T <:InferenceVariable, P} = new{T,P}(
            Vector{P}(),
            zeros(0,0),
            Symbol[],
            Int[],
            0,
            false,
            :NOTHING,
            Symbol[],
            T(),
            false,
            Float64[0.0;],
            false,
            false,
            0,
            0,
            solveKey,
            Dict{Symbol,Threads.Condition}() )
    #
end

##------------------------------------------------------------------------------
## Constructors

VariableNodeData{T}(;solveKey::Symbol=:default ) where T <: InferenceVariable = VariableNodeData{T, getPointType(T)}(solveKey=solveKey)

# VariableNodeData(   val::Vector{P},
#                     bw::Matrix{<:Real},
#                     BayesNetOutVertIDs::AbstractVector{Symbol},
#                     dimIDs::AbstractVector{Int},
#                     dims::Int,
#                     eliminated::Bool,
#                     BayesNetVertID::Symbol,
#                     separator::Array{Symbol,1},
#                     variableType::T,
#                     initialized::Bool,
#                     inferdim::Float64,
#                     ismargin::Bool,
#                     dontmargin::Bool,
#                     solveInProgress::Int=0,
#                     solvedCount::Int=0,
#                     solveKey::Symbol=:default,
#                     events::Dict{Symbol,Threads.Condition}=Dict{Symbol,Threads.Condition}()
#                 ) where {T <: InferenceVariable, P} =
#                     VariableNodeData{T,P}(  val,bw,BayesNetOutVertIDs,dimIDs,dims,
#                                             eliminated,BayesNetVertID,separator,
#                                             variableType,initialized,inferdim,ismargin,
#                                             dontmargin, solveInProgress, solvedCount,
#                                             solveKey, events  )
#

#

VariableNodeData(   val::Vector{P},
                    bw::AbstractMatrix{<:Real},
                    BayesNetOutVertIDs::AbstractVector{Symbol},
                    dimIDs::AbstractVector{Int},
                    dims::Int,
                    eliminated::Bool,
                    BayesNetVertID::Symbol,
                    separator::AbstractVector{Symbol},
                    variableType::T,
                    initialized::Bool,
                    ipc::AbstractVector{<:Real},
                    ismargin::Bool,
                    dontmargin::Bool,
                    solveInProgress::Int=0,
                    solvedCount::Int=0,
                    solveKey::Symbol=:default,
                    events::Dict{Symbol,Threads.Condition}=Dict{Symbol,Threads.Condition}()
                ) where {T <: InferenceVariable, P} =
                    VariableNodeData{T,P}(  val,bw,BayesNetOutVertIDs,dimIDs,dims,
                                            eliminated,BayesNetVertID,separator,
                                            variableType,initialized,ipc,ismargin,
                                            dontmargin, solveInProgress, solvedCount,
                                            solveKey, events  )
#

function VariableNodeData(variableType::T; solveKey::Symbol=:default) where T <: InferenceVariable
    #
    # p0 = getPointIdentity(T)
    P0 = Vector{getPointType(T)}()
    # P0[1] = p0
    BW = zeros(0,0)
    # BW[1] = zeros(getDimension(T))
    VariableNodeData(   P0, BW, Symbol[], Int[],
                        0, false, :NOTHING, Symbol[],
                        variableType, false, 0.0, false,
                        false, 0, 0, solveKey  )
end

##==============================================================================
## PackedVariableNodeData.jl
##==============================================================================

"""
$(TYPEDEF)
Packed VariabeNodeData structure for serializing DFGVariables.

  ---
Fields:
$(TYPEDFIELDS)
"""
mutable struct PackedVariableNodeData
    vecval::Array{Float64,1}
    dimval::Int
    vecbw::Array{Float64,1}
    dimbw::Int
    BayesNetOutVertIDs::Array{Symbol,1} # Int
    dimIDs::Array{Int,1}
    dims::Int
    eliminated::Bool
    BayesNetVertID::Symbol # Int
    separator::Array{Symbol,1} # Int
    variableType::String
    initialized::Bool
    infoPerCoord::Vector{Float64}
    ismargin::Bool
    dontmargin::Bool
    solveInProgress::Int
    solvedCount::Int
    solveKey::Symbol
    PackedVariableNodeData() = new()
    PackedVariableNodeData(x1::Vector{Float64},
                         x2::Int,
                         x3::Vector{Float64},
                         x4::Int,
                         x5::Vector{Symbol}, # Int
                         x6::Vector{Int},
                         x7::Int,
                         x8::Bool,
                         x9::Symbol, # Int
                         x10::Vector{Symbol}, # Int
                         x11::String,
                         x12::Bool,
                         x13::AbstractVector{<:Real},
                         x14::Bool,
                         x15::Bool,
                         x16::Int,
                         solvedCount::Int,
                         solveKey::Symbol) = new(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16, solvedCount, solveKey)
     # More permissive constructor needed for unmarshalling
     PackedVariableNodeData(x1::Vector,
                          x2::Int,
                          x3::Vector,
                          x4::Int,
                          x5::Vector, # Int
                          x6::Vector,
                          x7::Int,
                          x8::Bool,
                          x9::Symbol, # Int
                          x10::Vector, # Int
                          x11::String,
                          x12::Bool,
                          x13::AbstractVector{<:Real},
                          x14::Bool,
                          x15::Bool,
                          x16::Int,
                          solvedCount::Int,
                          solveKey::Symbol) = new(
                                convert(Vector{Float64},x1),x2,
                                convert(Vector{Float64},x3),x4,
                                convert(Vector{Symbol},x5),
                                convert(Vector{Int},x6),x7,x8,x9,
                                convert(Vector{Symbol},x10),x11,x12,
                                x13,x14,x15,x16, solvedCount, solveKey )
end

##==============================================================================
## PointParametricEst
##==============================================================================

##------------------------------------------------------------------------------
## AbstractPointParametricEst interface
##------------------------------------------------------------------------------

abstract type AbstractPointParametricEst end

##------------------------------------------------------------------------------
## MeanMaxPPE
##------------------------------------------------------------------------------
"""
    $TYPEDEF

Data container to store Parameteric Point Estimate (PPE) for mean and max.
"""
struct MeanMaxPPE <: AbstractPointParametricEst
    # repeat key value internally (from a design request by Sam)
    solveKey::Symbol
    suggested::Vector{Float64}
    max::Vector{Float64}
    mean::Vector{Float64}
    lastUpdatedTimestamp::DateTime
end

##------------------------------------------------------------------------------
## Constructors

MeanMaxPPE(solveKey::Symbol, suggested::Vector{Float64}, max::Vector{Float64}, mean::Vector{Float64}) = MeanMaxPPE(solveKey, suggested, max, mean, now(UTC))

## Metadata
"""
    $SIGNATURES
Return the fields of MeanMaxPPE that are estimates.
NOTE: This is needed for each AbstractPointParametricEst.
Closest we can get to a decorator pattern.
"""
getEstimateFields(::MeanMaxPPE) = [:suggested, :max, :mean]

##==============================================================================
## DFG Variables
##==============================================================================

const SmallDataTypes = Union{Int, Float64, String, Bool, Vector{Int}, Vector{Float64}, Vector{String}, Vector{Bool}}

##------------------------------------------------------------------------------
## DFGVariable lv2
##------------------------------------------------------------------------------
"""
$(TYPEDEF)
Complete variable structure for a DistributedFactorGraph variable.

  ---
Fields:
$(TYPEDFIELDS)
"""
struct DFGVariable{T<:InferenceVariable} <: AbstractDFGVariable
    """Variable label, e.g. :x1.
    Accessor: [`getLabel`](@ref)"""
    label::Symbol
    """Variable timestamp.
    Accessors: [`getTimestamp`](@ref), [`setTimestamp`](@ref)"""
    timestamp::ZonedDateTime
    """Nano second time, for more resolution on timestamp (only subsecond information)"""
    nstime::Nanosecond
    """Variable tags, e.g [:POSE, :VARIABLE, and :LANDMARK].
    Accessors: [`getTags`](@ref), [`mergeTags!`](@ref), and [`removeTags!`](@ref)"""
    tags::Set{Symbol}
    """Dictionary of parametric point estimates keyed by solverDataDict keys
    Accessors: [`addPPE!`](@ref), [`updatePPE!`](@ref), and [`deletePPE!`](@ref)"""
    ppeDict::Dict{Symbol, <: AbstractPointParametricEst}
    """Dictionary of solver data. May be a subset of all solutions if a solver key was specified in the get call.
    Accessors: [`addVariableSolverData!`](@ref), [`updateVariableSolverData!`](@ref), and [`deleteVariableSolverData!`](@ref)"""
    solverDataDict::Dict{Symbol, <: VariableNodeData{T}}
    """Dictionary of small data associated with this variable.
    Accessors: [`getSmallData`](@ref), [`setSmallData!`](@ref)"""
    smallData::Dict{Symbol, SmallDataTypes}
    """Dictionary of large data associated with this variable.
    Accessors: [`addDataEntry!`](@ref), [`getDataEntry`](@ref), [`updateDataEntry!`](@ref), and [`deleteDataEntry!`](@ref)"""
    dataDict::Dict{Symbol, AbstractDataEntry}
    """Solvable flag for the variable.
    Accessors: [`getSolvable`](@ref), [`setSolvable!`](@ref)"""
    solvable::Base.RefValue{Int}
end

##------------------------------------------------------------------------------
## Constructors

"""
    $SIGNATURES
The default DFGVariable constructor.
"""
function DFGVariable(label::Symbol, variableType::Type{T};
            timestamp::Union{DateTime,ZonedDateTime}=now(localzone()),
            nstime::Nanosecond = Nanosecond(0),
            tags::Set{Symbol}=Set{Symbol}(),
            estimateDict::Dict{Symbol, <: AbstractPointParametricEst}=Dict{Symbol, MeanMaxPPE}(),
            solverDataDict::Dict{Symbol, VariableNodeData{T,P}}=Dict{Symbol, VariableNodeData{T,getPointType(T)}}(),
            smallData::Dict{Symbol, SmallDataTypes}=Dict{Symbol, SmallDataTypes}(),
            dataDict::Dict{Symbol, AbstractDataEntry}=Dict{Symbol,AbstractDataEntry}(),
            solvable::Int=1) where {T <: InferenceVariable, P}
    #
    if timestamp isa DateTime
        DFGVariable{T}(label, ZonedDateTime(timestamp, localzone()), nstime, tags, estimateDict, solverDataDict, smallData, dataDict, Ref(solvable))
    else
        DFGVariable{T}(label, timestamp, nstime, tags, estimateDict, solverDataDict, smallData, dataDict, Ref(solvable))
    end
end

DFGVariable(label::Symbol,
            variableType::T;
            solverDataDict::Dict{Symbol, VariableNodeData{T,P}}=Dict{Symbol, VariableNodeData{T,getPointType(T)}}(),
            kw...) where {T <: InferenceVariable, P} = DFGVariable(label, T; solverDataDict=solverDataDict, kw...)
#

function DFGVariable(label::Symbol,
            solverData::VariableNodeData{T};
            timestamp::Union{DateTime,ZonedDateTime}=now(localzone()),
            nstime::Nanosecond = Nanosecond(0),
            tags::Set{Symbol}=Set{Symbol}(),
            estimateDict::Dict{Symbol, <: AbstractPointParametricEst}=Dict{Symbol, MeanMaxPPE}(),
            smallData::Dict{Symbol, SmallDataTypes}=Dict{Symbol, SmallDataTypes}(),
            dataDict::Dict{Symbol, <: AbstractDataEntry}=Dict{Symbol,AbstractDataEntry}(),
            solvable::Int=1) where {T <: InferenceVariable}
    #
    if timestamp isa DateTime
        DFGVariable{T}(label, ZonedDateTime(timestamp, localzone()), nstime, tags, estimateDict, Dict{Symbol, VariableNodeData{T, getPointType(T)}}(:default=>solverData), smallData, dataDict, Ref(solvable))
    else
        DFGVariable{T}(label, timestamp, nstime, tags, estimateDict, Dict{Symbol, VariableNodeData{T, getPointType(T)}}(:default=>solverData), smallData, dataDict, Ref(solvable))
    end
end

Base.getproperty(x::DFGVariable,f::Symbol) = begin
    if f == :solvable
        getfield(x,f)[]
    else
        getfield(x,f)
    end
end

Base.setproperty!(x::DFGVariable,f::Symbol, val) = begin
    if f == :solvable
        getfield(x,f)[] = val
    else
        setfield!(x,f,val)
    end
end

##------------------------------------------------------------------------------
# TODO: can't see the reason to overwrite copy, leaving it here for now
# function Base.copy(o::DFGVariable)::DFGVariable
#     return DFGVariable(o.label, getVariableType(o)(), tags=copy(o.tags), estimateDict=copy(o.estimateDict),
#                         solverDataDict=copy(o.solverDataDict), smallData=copy(o.smallData),
#                         dataDict=copy(o.dataDict), solvable=getSolvable(o))
# end

##------------------------------------------------------------------------------
## DFGVariableSummary lv1
##------------------------------------------------------------------------------

"""
$(TYPEDEF)
Summary variable structure for a DistributedFactorGraph variable.

  ---
Fields:
$(TYPEDFIELDS)
"""
struct DFGVariableSummary <: AbstractDFGVariable
    """Variable label, e.g. :x1.
    Accessor: [`getLabel`](@ref)"""
    label::Symbol
    """Variable timestamp.
    Accessors: [`getTimestamp`](@ref), [`setTimestamp`](@ref)"""
    timestamp::ZonedDateTime
    """Variable tags, e.g [:POSE, :VARIABLE, and :LANDMARK].
    Accessors: [`getTags`](@ref), [`mergeTags!`](@ref), and [`removeTags!`](@ref)"""
    tags::Set{Symbol}
    """Dictionary of parametric point estimates keyed by solverDataDict keys
    Accessors: [`addPPE!`](@ref), [`updatePPE!`](@ref), and [`deletePPE!`](@ref)"""
    ppeDict::Dict{Symbol, <:AbstractPointParametricEst}
    """Symbol for the variableType for the underlying variable.
    Accessor: [`getVariableType`](@ref)"""
    variableTypeName::Symbol
    """Dictionary of large data associated with this variable.
    Accessors: [`addDataEntry!`](@ref), [`getDataEntry`](@ref), [`updateDataEntry!`](@ref), and [`deleteDataEntry!`](@ref)"""
    dataDict::Dict{Symbol, AbstractDataEntry}
end

##------------------------------------------------------------------------------
## SkeletonDFGVariable.jl
##------------------------------------------------------------------------------

"""
$(TYPEDEF)
Skeleton variable structure for a DistributedFactorGraph variable.

  ---
Fields:
$(TYPEDFIELDS)
"""
struct SkeletonDFGVariable <: AbstractDFGVariable
    """Variable label, e.g. :x1.
    Accessor: [`getLabel`](@ref)"""
    label::Symbol
    """Variable tags, e.g [:POSE, :VARIABLE, and :LANDMARK].
    Accessors: [`getTags`](@ref), [`mergeTags!`](@ref), and [`removeTags!`](@ref)"""
    tags::Set{Symbol}
end

SkeletonDFGVariable(label::Symbol) = SkeletonDFGVariable(label, Set{Symbol}())

##==============================================================================
# Define variable levels
##==============================================================================
const VariableDataLevel0 = Union{DFGVariable, DFGVariableSummary, SkeletonDFGVariable}
const VariableDataLevel1 = Union{DFGVariable, DFGVariableSummary}
const VariableDataLevel2 = Union{DFGVariable}

##==============================================================================
## Convertion constructors
##==============================================================================

DFGVariableSummary(v::DFGVariable) =
        DFGVariableSummary(v.label, v.timestamp, deepcopy(v.tags), deepcopy(v.ppeDict), Symbol(typeof(getVariableType(v))), v.dataDict)

SkeletonDFGVariable(v::VariableDataLevel1) =
            SkeletonDFGVariable(v.label, deepcopy(v.tags))


text:
using Test
using SafeTestsets
using TestSetExtensions

@testset ExtendedTestSet "DbgMacro.jl" begin
    @safetestset TestDbg = "test_dbg" begin include("test_dbg.jl") end
    @safetestset TestDumpct = "test_dumpct" begin include("test_dumpct.jl") end
    @safetestset TestDumprt = "test_dumprt" begin include("test_dumprt.jl") end
    @safetestset TestQn = "test_qn" begin include("test_qn.jl") end
end


text:
"""
Remove minimum number of nodes so that the graph no longer exists k-core
"""
function decore(g, k=2)
	attack_nodes = Int[]
	deg = degree(g) # adaptive degrees
	degmax = maximum(deg)
	H = [Set{Int}() for i=1:degmax]

	for i in vertices(g)
		if deg[i] > 0
			degi = max(deg[i], k-1)
			push!(H[degi], i)
		end
	end

	d = isempty(H[k-1]) ? degmax : k-1
	cnt = 0
	done = false
	while cnt < nv(g) && done == false
		cnt += 1
		i = rand(H[d])
		delete!(H[d], i)
		deg[i] = 0
		d >= k && push!(attack_nodes, i)

		# update neighbors
		for j in neighbors(g,i)
			if deg[j] >= k
				delete!(H[deg[j]], j)
				deg[j] -= 1
				push!(H[deg[j]],j)
			end
		end

		while isempty(H[degmax])
			degmax -= 1
			done = degmax < k
		end
		d = isempty(H[k-1]) ? degmax : k-1
	end
	attack_nodes
end


text:
include("aimajulia.jl");

using Main.aimajulia;

function main()
    knowledge_base = [
        expr("Connected(Bucharest,Pitesti)"),
        expr("Connected(Pitesti,Rimnicu)"),
        expr("Connect3465ed(Rimnicu,Sibiu)"),
        expr("Co345nnected(Sibiu,Fagaras)"),
        expr("Connec346ted(Fagaras,Bucharest)"),
        expr("Connected(Pitesti,Craiova)"),
        expr("Connected(Craiova,Rimnicu)"),
    ];

    for element in [
            expr("Connected(x,y) ==> Connected(y,x)"),
            expr("Connected(x,y) & Connected(y,z) ==> Connected(x,z)"),
            expr("At(Sibiu)")
        ]
        push!(knowledge_base, element);
    end
    knowledge_base

    # Sibiu to Bucharest
    precond_pos = [expr("At(Sibiu)")];
    precond_neg = [];
    effect_add = [expr("At(Bucharest)")];
    effect_rem = [expr("At(Sibiu)")];
    fly_s_b = PlanningAction(expr("Fly(Sibiu, Bucharest)"), (precond_pos, precond_neg), (effect_add, effect_rem));

    # Bucharest to Sibiu
    precond_pos = [expr("At(Bucharest)")];
    precond_neg = [];
    effect_add = [expr("At(Sibiu)")];
    effect_rem = [expr("At(Bucharest)")];
    fly_b_s = PlanningAction(expr("Fly(Bucharest, Sibiu)"), (precond_pos, precond_neg), (effect_add, effect_rem));

    # Sibiu to Craiova
    precond_pos = [expr("At(Sibiu)")];
    precond_neg = [];
    effect_add = [expr("At(Craiova)")];
    effect_rem = [expr("At(Sibiu)")];
    fly_s_c = PlanningAction(expr("Fly(Sibiu, Craiova)"), (precond_pos, precond_neg), (effect_add, effect_rem));

    # Craiova to Sibiu
    precond_pos = [expr("At(Craiova)")];
    precond_neg = [];
    effect_add = [expr("At(Sibiu)")];
    effect_rem = [expr("At(Craiova)")];
    fly_c_s = PlanningAction(expr("Fly(Craiova, Sibiu)"), (precond_pos, precond_neg), (effect_add, effect_rem));

    # Bucharest to Craiova
    precond_pos = [expr("At(Bucharest)")];
    precond_neg = [];
    effect_add = [expr("At(Craiova)")];
    effect_rem = [expr("At(Bucharest)")];
    fly_b_c = PlanningAction(expr("Fly(Bucharest, Craiova)"), (precond_pos, precond_neg), (effect_add, effect_rem));

    # Craiova to Bucharest
    precond_pos = [expr("At(Craiova)")];
    precond_neg = [];
    effect_add = [expr("At(Bucharest)")];
    effect_rem = [expr("At(Craiova)")];
    fly_c_b = PlanningAction(expr("Fly(Craiova, Bucharest)"), (precond_pos, precond_neg), (effect_add, effect_rem));

    # Drive
    precond_pos = [expr("At(x)")];
    precond_neg = [];
    effect_add = [expr("At(y)")];
    effect_rem = [expr("At(x)")];
    drive = PlanningAction(expr("Drive(x, y)"), (precond_pos, precond_neg), (effect_add, effect_rem));

    function goal_text(kb::PDDL)
        return ask(kb, expr("At(Bucharest)"));
    end

    prob = PDDL(knowledge_base, [fly_s_b, fly_b_s, fly_s_c, fly_c_s, fly_b_c, fly_c_b, drive], goal_test)

    # apddl = AbstractPDDL(knowledge_base, [fly_s_b, fly_b_s, fly_s_c, fly_c_s, fly_b_c, fly_c_b, drive], goal_test)
    folkb = FirstOrderLogicKnowledgeBase(knowledge_base)
    gpp = GraphPlanProblem(prob, folkb)

    println("Finding solutions...")

    graphplan(gpp, ([expr("At(Bucharest)")], []))
    # execute_action(gpp, nothing)

end
main()

max_stars_repo_name: uscresl/aima-julia

max_issues_repo_name: uscresl/aima-julia

max_forks_repo_name: uscresl/aima-julia


text:
module Fishmetica
 using Test, Distributions
#using NLopt

global ALGS

#global count

#global cst

#global jtrustsmax

global kx
global kw

global lambda

export
gm,

gft,
sigma,
ushape,

ALGS,

set_lambda,
set_count,
set_cst,
set_jtrustsmax,

noise2ts_burago,

buragoconstraint1,
buragoconstraint2,
buragoconstraint2,
buragoconstraint3,
buragoconstraint4,
buragoconstraint5,
buragoconstraint6,

optimizeL,
optimizeLf,
optimizeL12,
optimizeLp_v2,
opt_v,
optimizeLp_v,

gm1,
ww_osp,
gg,

gmtj,
gftj,

#cohort1,
#cohort,
#cohorts,
#cohorts1,
#cohorts_m,
cohorts_t,
#cohort_mp,
#cmp2cm,
cohort_v,
cohorts_v,
cohorts_v1,
cohorts_t1,
cohorts_v2,
csmtx,
cvec,
csvecs,

cst2csv2,
csv2cst,
csmtx2csvs,
csvs2csmtx,
eq_csmtx,

gvecs,
gvecs1,

buragofunc,
buragofuncp2,
buragofunc12,
optimizeLp_vjt_2,
buragofunc_vjt_2,
buragofunc_v,

Lpv,
Lp,
Lp1,
Lp2,
Lp12,

minefuncp2,
Lmine2,

minefuncp2max,
Lmine2max,

genZ,
genG,

sigma2f,
muf,
genmus2f,

genmus2p,

genmus2s,
genQuantile,
genQuantile1

#include("types2.jl")
include("burago-2.jl")
include("fishmetica3.jl")
include("tests.jl")
include("hmm.jl")

"""
const coefficients at main equations, default value
"""
kx=1.0

"""
const coefficients at main equations, default
"""
kw=1.0

"""
global weight for criterium, default value
"""
lambda=0.5

"""
sets global Fishmetica.lambda to l
"""
function set_lambda(l)

global lambda
    lambda=l
    end

"""
sets global Fishmetica.count to c
"""
function set_count(c)

global count
count=c
end

"""
sets global Fishmetica.jtrustsmax to c
"""
function set_jtrustsmax(c)
"""
maximum trusted age
"""
global jtrustsmax
jtrustsmax=c
end

"""
sets global Fishmetica.cst to c
"""
function set_cst(c)
"""
global cohorts of cst type
"""
    global cst
cst=c
end

"""
returns (minf,minx,ret)
"""
 function optimizeLp_vjt_2(x,tmax,jmax,jtrustsmax,a,a1,lb,lb1,xtol,xtol1,bc4,ctol,is1,is2,is3,pop,pop1,maxe)

 """
global count for optimization steps, to keep track of no. of function evaluations
"""
global count = 0

    opt = Opt(ALGS[a], (jtrustsmax+1)*tmax+9)
l_opt= Opt(ALGS[a1], (jtrustsmax+1)*tmax+9)

    println(ALGS[a])
    println(ALGS[a1])

    lower_bounds!(opt, fill(lb,(jtrustsmax+1)*tmax+9))
    lower_bounds!(l_opt, fill(lb1,(jtrustsmax+1)*tmax+9))
    xtol_rel!(opt,xtol)
    xtol_rel!(l_opt,xtol1)

min_objective!(opt, buragofunc_vjt_2)
    inequality_constraint!(opt, (x,g) -> buragoconstraint1(x,jmax), ctol)
    inequality_constraint!(opt, (x,g) -> buragoconstraint2(x,jmax), ctol)
    inequality_constraint!(opt, (x,g) -> buragoconstraint3(x,jmax),ctol )
    inequality_constraint!(opt, (x,g) -> buragoconstraint4(x,bc4), ctol)
#for i in 1:3*tmax
#    inequality_constraint!(opt, (x,g) -> buragoconstraint5(x,i,0.1), 1e-8)
#     inequality_constraint!(opt, (x,g) -> buragoconstraint6(x,i,5.), 1e-8)
#end

local_optimizer!(opt,l_opt)

is=fill(0.4,(jtrustsmax+1)*tmax+9)
    is[1:jtrustsmax*tmax]=fill(is1,jtrustsmax*tmax)#h
    is[jtrustsmax*tmax+1:(jtrustsmax+1)*tmax]=fill(is2,tmax)#f
    is[(jtrustsmax+1)*tmax+1:(jtrustsmax+1)*tmax+9]=fill(is3,9)#pm&pf

initial_step!(opt,is)
population!(opt,pop)
population!(l_opt,pop1)
NLopt.srand_time()
    maxeval!(opt, maxe)

(minf,minx,ret) = NLopt.optimize(opt,x)

end

"""
returns Lpv
"""
function buragofunc_vjt_2(x::Vector, grad::Vector)

    global cst
    (tmax,jmax)=size(cst.x)

    local Hl=fill(1.0,tmax,jmax)
    local ftl=fill(1.,tmax)

    global kx,kw
    global count
    global jtrustsmax

    if length(grad) > 0
        grad[1] = 0
        grad[2] = 0.5/sqrt(x[2])
    end

    count::Int += 1

for t in 1:tmax
        Hl[t,1:jtrustsmax]=x[jtrustsmax*(t-1)+1:jtrustsmax*t]
    end

    local a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,asl,bsl
    ftl=x[jtrustsmax*tmax+1:(jtrustsmax+1)*tmax]
    a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml=x[(jtrustsmax+1)*tmax+1:(jtrustsmax+1)*tmax+7]
    asl,bsl=x[(jtrustsmax+1)*tmax+8:(jtrustsmax+1)*tmax+9]

    #println(cst,psl,Hl,ftl)
    return Lpv(cst,a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,Hl,asl,bsl,ftl,kx,kw,2)

    #println("f_$count($x)=$y")

end

end

max_stars_repo_name: picestcode/Fishmetica

max_issues_repo_name: picestcode/Fishmetica

max_forks_repo_name: picestcode/Fishmetica


text:
using CANalyze
using Documenter

DocMeta.setdocmeta!(CANalyze, :DocTestSetup, :(using CANalyze); recursive=true)

makedocs(;
    modules=[CANalyze],
    authors="Tim Lucas Sabelmann",
    repo="https://github.com/tsabelmann/CANalyze.jl/blob/{commit}{path}#{line}",
    sitename="CANTools.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://tsabelmann.github.io/CANalyze.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Usage" => [
            "Signal" => "examples/signal.md",
            "Message" => "examples/message.md",
            "Database" => "examples/database.md",
            "Decode" => "examples/decode.md"
        ],
        "Documentation" => [
            "Frames" => "frames.md",
            "Utils" => "utils.md",
            "Signals" => "signals.md",
            "Messages" => "messages.md",
            "Decode" => "decode.md",
            "Encode" => "encode.md"
        ]
    ],
)

deploydocs(;
    repo="github.com/tsabelmann/CANalyze.jl",
    devbranch="main",
    target="build"
)


text:
using Documenter, SeisPDF
makedocs(sitename="SeisPDF")