text:
# This file is part of Kpax3. License is MIT.

abstract type State end

mutable struct AminoAcidState <: State
  R::Vector{Int}
  C::Matrix{UInt8}

  emptycluster::BitArray{1}
  cl::Vector{Int}
  k::Int

  v::Vector{Int}
  n1s::Matrix{Float64}
  unit::Vector{Vector{Int}}

  logpR::Float64
  logpC::Vector{Float64}
  loglik::Float64
  logpp::Float64
end

function AminoAcidState(data::Matrix{UInt8},
                        partition::Vector{Int},
                        priorR::PriorRowPartition,
                        priorC::AminoAcidPriorCol,
                        settings::KSettings)
  (m, n) = size(data)

  R = normalizepartition(partition, n)
  k = maximum(R)

  maxclust = min(n, max(k, settings.maxclust))

  C = zeros(UInt8, maxclust, m)

  emptycluster = trues(maxclust)
  cl = zeros(Int, maxclust)

  v = zeros(Int, maxclust)
  n1s = zeros(Float64, maxclust, m)
  unit = Vector{Int}[zeros(Int, settings.maxunit) for g in 1:maxclust]

  g = 0
  for a in 1:n
    g = R[a]

    if emptycluster[g]
      emptycluster[g] = false
    end

    if v[g] == length(unit[g])
      resize!(unit[g], min(n, v[g] + settings.maxunit))
    end

    v[g] += 1
    unit[g][v[g]] = a

    for b in 1:m
      n1s[g, b] += float(data[b, a])
    end
  end

  # we do it here and not in the previous loop because we want them sorted
  # it is not really necessary but they will be sorted during the simulation
  # anyway (we will loop on emptycluster from now on)
  i = 1
  for a in 1:length(emptycluster)
    if !emptycluster[a]
      cl[i] = a
      i += 1
    end
  end

  logpR = logdPriorRow(n, k, v, priorR)

  logpC = zeros(Float64, 2)
  computelocalmode!(v, n1s, C, cl, k, logpC, priorC)

  loglik = loglikelihood(C, cl, k, v, n1s, priorC)

  logpp = logpR + logpC[1] + loglik

  AminoAcidState(R, C, emptycluster, cl, k, v, n1s, unit, logpR, logpC, loglik,
                 logpp)
end

function copystate(x::AminoAcidState)
  R = copy(x.R)
  C = copy(x.C)

  emptycluster = copy(x.emptycluster)
  cl = copy(x.cl)
  k = x.k
  v = copy(x.v)
  n1s = copy(x.n1s)

  unit = Array{Vector{Int}}(undef, length(x.unit))
  for l in 1:length(x.unit)
    unit[l] = copy(x.unit[l])
  end

  logpR = x.logpR
  logpC = copy(x.logpC)
  loglik = x.loglik
  logpp = x.logpp

  AminoAcidState(R, C, emptycluster, cl, k, v, n1s, unit, logpR, logpC, loglik,
                 logpp)
end

function copystate!(dest::AminoAcidState,
                    src::AminoAcidState)
  resizestate!(dest, src.cl[src.k])

  if length(dest.R) == length(src.R)
    copyto!(dest.R, src.R)
  else
    dest.R = copy(src.R)
  end

  if size(dest.C, 2) != size(src.C, 2)
    dest.C = zeros(UInt8, size(dest.C, 1), size(src.C, 2))
    dest.n1s = zeros(Float64, size(dest.n1s, 1), size(src.n1s, 2))
  end

  fill!(dest.emptycluster, true)

  g = 0
  for l in 1:src.k
    g = src.cl[l]

    dest.C[g, 1] = src.C[g, 1]
    dest.emptycluster[g] = false
    dest.cl[l] = src.cl[l]
    dest.v[g] = src.v[g]
    dest.n1s[g, 1] = src.n1s[g, 1]

    if length(dest.unit[g]) < src.v[g]
      resize!(dest.unit[g], src.v[g])
    end

    copyto!(dest.unit[g], 1, src.unit[g], 1, src.v[g])
  end

  for b in 2:size(dest.C, 2)
    for l in 1:src.k
      g = src.cl[l]

      dest.C[g, b] = src.C[g, b]
      dest.n1s[g, b] = src.n1s[g, b]
    end
  end

  dest.k = src.k

  dest.logpR = src.logpR

  dest.logpC[1] = src.logpC[1]
  dest.logpC[2] = src.logpC[2]

  dest.loglik = src.loglik

  dest.logpp = src.logpp

  nothing
end

function initializestate(data::Matrix{UInt8},
                         D::Matrix{Float64},
                         kset::UnitRange{Int},
                         priorR::PriorRowPartition,
                         priorC::PriorColPartition,
                         settings::KSettings)
  n = size(data, 2)

  R = ones(Int, n)

  s = AminoAcidState(data, R, priorR, priorC, settings)

  t1 = copystate(s)
  t2 = copystate(s)

  # TODO: remove the hack once kmedoids is fixed

  if settings.verbose
    Printf.@printf("Log-posterior (plus a constant) for one cluster: %.4f.\n",
                   s.logpp)
    Printf.@printf("Now scanning %d to %d clusters.\n", kset[1], kset[end])
  end

  niter = 0
  for k in kset
    if settings.verbose && (k % 10 == 0)
      Printf.@printf("Total number of clusters = %d.\n", k)
    end

    fill!(R, 0)
    try
      copyto!(R, Clustering.kmedoids(D, k).assignments)
    catch
      StatsBase.sample!(1:k, R, replace=true)
      R[StatsBase.sample(1:n, k, replace=false)] = collect(1:k)
    end

    updatestate!(t1, data, R, priorR, priorC, settings)

    niter = 0
    while niter < 10
      try
        copyto!(R, Clustering.kmedoids(D, k).assignments)
      catch
        StatsBase.sample!(1:k, R, replace=true)
        R[StatsBase.sample(1:n, k, replace=false)] = collect(1:k)
      end

      updatestate!(t2, data, R, priorR, priorC, settings)

      if t2.logpp > t1.logpp
        copystate!(t1, t2)
      end

      niter += 1
    end

    if t1.logpp > s.logpp
      copystate!(s, t1)

      if settings.verbose
        Printf.@printf("Found a better solution! ")
        Printf.@printf("Log-posterior (plus a constant) for %d clusters: %.4f.\n",
                       k, s.logpp)
      end
    end
  end

  s
end

function updatestate!(state::AminoAcidState,
                      data::Matrix{UInt8},
                      R::Vector{Int},
                      priorR::PriorRowPartition,
                      priorC::AminoAcidPriorCol,
                      settings::KSettings)
  (m, n) = size(data)

  copyto!(state.R, normalizepartition(R, n))
  k = maximum(state.R)

  resizestate!(state, k, settings)

  state.k = k

  fill!(state.emptycluster, true)
  fill!(state.v, 0)
  fill!(state.n1s, 0.0)

  g = 0
  for a in 1:n
    g = state.R[a]

    if state.emptycluster[g]
      state.emptycluster[g] = false
    end

    if state.v[g] == length(state.unit[g])
      resize!(state.unit[g], min(n, state.v[g] + settings.maxunit))
    end

    state.v[g] += 1
    state.unit[g][state.v[g]] = a

    for b in 1:m
      state.n1s[g, b] += float(data[b, a])
    end
  end

  i = 1
  for a in 1:length(state.emptycluster)
    if !state.emptycluster[a]
      state.cl[i] = a
      i += 1
    end
  end

  state.logpR = logdPriorRow(n, state.k, state.v, priorR)

  computelocalmode!(state.v, state.n1s, state.C, state.cl, state.k,
                    state.logpC, priorC)

  state.loglik = loglikelihood(state.C, state.cl, state.k, state.v, state.n1s,
                               priorC)

  state.logpp = state.logpR + state.logpC[1] + state.loglik

  nothing
end

function resizestate!(state::AminoAcidState,
                      k::Int,
                      settings::KSettings)
  oldlen = size(state.C, 1)

  if oldlen < k
    m = size(state.C, 2)
    n = length(state.R)

    # we don't want to allocate new resources too often, so allocate double the
    # previous size. this should guarantee a logarithmic number of allocations
    newlen = min(n, max(k, 2 * oldlen))

    resize!(state.emptycluster, newlen)
    resize!(state.cl, newlen)
    resize!(state.v, newlen)
    resize!(state.unit, newlen)

    for l in (oldlen + 1):newlen
      state.emptycluster[l] = true
      state.cl[l] = 0
      state.v[l] = 0
      state.unit[l] = zeros(Int, settings.maxunit)
    end

    C = zeros(UInt8, newlen, m)
    n1s = zeros(Float64, newlen, m)

    for b in 1:m
      for l in 1:state.k
        g = state.cl[l]
        C[g, b] = state.C[g, b]
        n1s[g, b] = state.n1s[g, b]
      end
    end

    state.C = C
    state.n1s = n1s
  end

  nothing
end

function resizestate!(state::AminoAcidState,
                      k::Int)
  oldlen = size(state.C, 1)

  if oldlen < k
    m = size(state.C, 2)
    n = length(state.R)

    # we don't want to allocate new resources too often, so allocate double the
    # previous size. this should guarantee a logarithmic number of allocations
    newlen = min(n, max(k, 2 * oldlen))

    resize!(state.emptycluster, newlen)
    resize!(state.cl, newlen)
    resize!(state.v, newlen)
    resize!(state.unit, newlen)

    for l in (oldlen + 1):newlen
      state.emptycluster[l] = true
      state.cl[l] = 0
      state.v[l] = 0
      state.unit[l] = zeros(Int, 1)
    end

    C = zeros(UInt8, newlen, m)
    n1s = zeros(Float64, newlen, m)

    for b in 1:m
      for l in 1:state.k
        g = state.cl[l]
        C[g, b] = state.C[g, b]
        n1s[g, b] = state.n1s[g, b]
      end
    end

    state.C = C
    state.n1s = n1s
  end

  nothing
end

function optimumstate(x::AminoAcidData,
                      partition::String,
                      settings::KSettings)
  R = normalizepartition(partition, x.id)

  priorR = EwensPitman(settings.α, settings.θ)
  priorC = AminoAcidPriorCol(x.data, settings.γ, settings.r)

  AminoAcidState(x.data, R, priorR, priorC, settings)
end


text:
## superseded in CUDA 11

@checked function cuDevicePrimaryCtxRelease(dev)
    @runtime_ccall((:cuDevicePrimaryCtxRelease, libcuda()), CUresult,
                   (CUdevice,),
                   dev)
end

@checked function cuDevicePrimaryCtxSetFlags(dev, flags)
    @runtime_ccall((:cuDevicePrimaryCtxSetFlags, libcuda()), CUresult,
                   (CUdevice, UInt32),
                   dev, flags)
end

@checked function cuDevicePrimaryCtxReset(dev)
    @runtime_ccall((:cuDevicePrimaryCtxReset, libcuda()), CUresult,
                   (CUdevice,),
                   dev)
end

@checked function cuGraphInstantiate(phGraphExec, hGraph, phErrorNode, logBuffer,
                                        bufferSize)
    initialize_api()
    @runtime_ccall((:cuGraphInstantiate, libcuda()), CUresult,
                   (Ptr{CUgraphExec}, CUgraph, Ptr{CUgraphNode}, Cstring, Csize_t),
                   phGraphExec, hGraph, phErrorNode, logBuffer, bufferSize)
end

ext: jl


text:
#this is a stub function.  It needs to be tested before being officially accepted.
@generated function __round(x::Sigmoid{N, ES, :ubit}, extrabits::(@UInt) = zero(@UInt)) where {N,ES}
  #for convenience, store the type of x as T.
  T = x
  z = zero(@UInt)

  innerbit = one(@UInt) << (__BITS - N)
  checkmask = (-(@UInt)(1)) >> (N)
  blankmask = ~((-(@UInt)(1)) >> N)
  quote
    #take exact zero and exact infinity.
    if (extrabits == zero(@UInt))
      (@u(x) == @signbit) && return reinterpret($T, @signbit)
      (@u(x) == zero(@UInt)) && return reinterpret($T, zero(@UInt))
    end

    truncated_value = @u(x) & $blankmask

    #eliminate values that go towards infinity - as "maxreal"
    (truncated_value == (@signbit) - $innerbit) && return reinterpret($T, truncated_value)
    #eliminate values that go down to infinity.
    (truncated_value == (@signbit)) && return reinterpret($T, truncated_value + $innerbit)
    #eliminate values that go down to zero.
    (truncated_value == zero(@UInt)) && return reinterpret($T, $innerbit)
    #eliminate values that go up to zero.
    (truncated_value == -$innerbit) && return reinterpret($T, -$innerbit)

    #first look to see if the checkbit is set.  If it's zero, round down.
    ($checkmask & @u(x) != $z) && return reinterpret($T, truncated_value | $innerbit)
    return reinterpret($T, truncated_value)
  end
end

#unify the three modes because we'll use the same mode for all of them.
RoundedSigmoid{N,ES} = Union{Sigmoid{N,ES,:upper}, Sigmoid{N,ES,:lower}, Sigmoid{N,ES,:exact}, Sigmoid{N,ES,:outward_exact}, Sigmoid{N,ES,:inward_exact}}
Direct_Sigmoid{N,ES} = Union{Sigmoid{N,ES,:outward_ulp}, Sigmoid{N,ES,:inward_ulp}}

@generated function __round(x::Union{RoundedSigmoid{N, ES}, Direct_Sigmoid{N,ES}}, extrabits::(@UInt) = zero(@UInt)) where {N,ES}
  #for convenience, store the type of x as T.
  T = x
  z = zero(@UInt)

  innerbit = one(@UInt) << (__BITS - N)
  checkmask = (-(@UInt)(1)) >> (N)
  blankmask = ~((-(@UInt)(1)) >> N)
  quote
    #take exact zero and exact infinity.
    if (extrabits == zero(@UInt))
      (@u(x) == @signbit) && return reinterpret($T, @signbit)
      (@u(x) == zero(@UInt)) && return reinterpret($T, zero(@UInt))
    end

    truncated_value = @u(x) & $blankmask

    #eliminate values that go towards infinity - as "maxreal"
    (truncated_value == (@signbit) - $innerbit) && return reinterpret($T, truncated_value)
    #eliminate values that go down to infinity.
    (truncated_value == (@signbit)) && return reinterpret($T, truncated_value + $innerbit)
    #eliminate values that go down to zero.
    (truncated_value == zero(@UInt)) && return reinterpret($T, $innerbit)
    #eliminate values that go up to zero.
    (truncated_value == -$innerbit) && return reinterpret($T, -$innerbit)

    #first look to see if the checkbit is set.  If it's zero, round down.
    ($checkmask & @u(x) != $z) && return reinterpret($T, truncated_value | $innerbit)
    return reinterpret($T, truncated_value)
  end
end

resolve_rounding(x::Sigmoid{N,ES,:lower}) where {N,ES}         = upper_ulp(reinterpret(Sigmoid{N,ES,:ubit}, x))
resolve_rounding(x::Sigmoid{N,ES,:upper}) where {N,ES}         = lower_ulp(reinterpret(Sigmoid{N,ES,:ubit}, x))
resolve_rounding(x::Sigmoid{N,ES,:exact}) where {N,ES}         = reinterpret(Sigmoid{N,ES,:ubit}, x)
resolve_rounding(x::Sigmoid{N,ES,:outward_exact}) where {N,ES} = reinterpret(Sigmoid{N,ES,:ubit}, x)
resolve_rounding(x::Sigmoid{N,ES,:inward_exact}) where {N,ES}  = reinterpret(Sigmoid{N,ES,:ubit}, x)

#special resolve_rounding functions which also must take location into account for the inf/zero values.
@generated function resolve_rounding(x::Sigmoid{N,ES,mode}, ::Type{Val{location}}) where {N,ES,mode,location}
    inf_ulp  = (mode == :outward_ulp) ? :(throw(ArgumentError("an outward ulp cannot be inf"))) : (location == :upper ? :maxpos : :minneg)
    zero_ulp = (mode == :inward_ulp)  ? :(throw(ArgumentError("an inward ulp cannot be zero"))) : (location == :upper ? :maxneg : :minpos)
    negative_ulp = (mode == :outward_ulp) ? :lower_ulp : :upper_ulp
    positive_ulp = (mode == :outward_ulp) ? :upper_ulp : :lower_ulp
    quote
        if !isfinite(x)
            $(inf_ulp)(Vnum{N,ES})
        elseif Base.iszero(x)
            $(zero_ulp)(Vnum{N,ES})
        elseif x < zero(x)
            $(negative_ulp)(reinterpret(Vnum{N,ES},x))
        else
            $(positive_ulp)(reinterpret(Vnum{N,ES},x))
        end
    end
end

→(lower::RoundedSigmoid{N,ES}, upper::RoundedSigmoid{N,ES}) where {N,ES} = resolve_rounding(lower)              → resolve_rounding(upper)
→(lower::RoundedSigmoid{N,ES}, upper::Direct_Sigmoid{N,ES}) where {N,ES} = resolve_rounding(lower)              → resolve_rounding(upper, Val{:upper})
→(lower::Direct_Sigmoid{N,ES}, upper::RoundedSigmoid{N,ES}) where {N,ES} = resolve_rounding(lower, Val{:lower}) → resolve_rounding(upper)
→(lower::Direct_Sigmoid{N,ES}, upper::Direct_Sigmoid{N,ES}) where {N,ES} = resolve_rounding(lower, Val{:lower}) → resolve_rounding(upper, Val{:upper})

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
# Newton with line search
function newton!(mechanism::Mechanism{T,Nl,0}; ε = 1e-10, newtonIter = 100, lineIter = 10, warning::Bool = false) where {T,Nl}
    bodies = mechanism.bodies
    eqcs = mechanism.eqconstraints
    graph = mechanism.graph
    ldu = mechanism.ldu
    Δt = mechanism.Δt

    normf0 = normf(mechanism)
    for n = Base.OneTo(newtonIter)
        setentries!(mechanism)
        factor!(graph, ldu)
        solve!(mechanism) # x̂1 for each body and constraint

        normf1 = lineSearch!(mechanism, normf0;iter = lineIter, warning = warning)
        normΔs1 = normΔs(mechanism)

        foreach(s1tos0!, bodies)
        foreach(s1tos0!, eqcs)

        if normf1 < ε && normΔs1 < ε
            # warning && (@info string("Newton iterations: ",n))
            return
        else
            normf0 = normf1
        end
    end

    warning && (@info string("newton! did not converge. n = ", newtonIter, ", tol = ", normf(mechanism), "."))
    return
end

# Newton on interior point with line search
function newton!(mechanism::Mechanism{T,Nl}; ε = 1e-10, σ = 0.1, μ = 1.0, newtonIter = 100, lineIter = 10, warning::Bool = false) where {T,Nl}
    bodies = mechanism.bodies
    eqcs = mechanism.eqconstraints
    ineqcs = mechanism.ineqconstraints
    graph = mechanism.graph
    ldu = mechanism.ldu
    Δt = mechanism.Δt

    foreach(resetVars!, ineqcs)
    mechanism.μ = μ

    meritf0 = meritf(mechanism)
    for n = Base.OneTo(newtonIter)
        setFrictionForce!(mechanism)
        setentries!(mechanism)
        factor!(graph, ldu)
        solve!(mechanism) # x̂1 for each body and constraint

        meritf1 = lineSearch!(mechanism, meritf0;iter = lineIter, warning = warning)

        foreach(s1tos0!, bodies)
        foreach(s1tos0!, eqcs)
        foreach(s1tos0!, ineqcs)

        if normf(mechanism) < ε
            warning && (@info string("Newton iterations: ", n))
            return
        else
            while meritf1 < mechanism.μ
                mechanism.μ = σ * mechanism.μ
                meritf1 = meritf(mechanism)
            end
            meritf0 = meritf1
        end
    end

    warning && (@info string("newton! did not converge. n = ", newtonIter, ", tol = ", normf(mechanism), "."))
    return
end

function lineSearch!(mechanism::Mechanism{T,N,0}, normf0;iter = 10, warning::Bool = false) where {T,N}
    normf1 = normf0
    scale = 0
    ldu = mechanism.ldu
    bodies = mechanism.bodies
    eqcs = mechanism.eqconstraints

    for n = Base.OneTo(iter + 1)
        for body in bodies
            lineStep!(body, getentry(ldu, body.id), scale)
        end
        for eqc in eqcs
            lineStep!(eqc, getentry(ldu, eqc.id), scale)
        end

        normf1 = normf(mechanism)
        if normf1 >= normf0
            scale += 1
        else
            return normf1
        end
    end

    warning && (@info string("lineSearch! did not converge. n = ", iter, ". Last tol: ", normf1))
    return normf1
end

function lineSearch!(mechanism::Mechanism, meritf0;iter = 10, warning::Bool = false)
    meritf1 = meritf0
    scale = 0
    ldu = mechanism.ldu
    bodies = mechanism.bodies
    eqcs = mechanism.eqconstraints
    ineqcs = mechanism.ineqconstraints

    feasibilityStepLength!(mechanism)

    for n = Base.OneTo(iter)
        for body in bodies
            lineStep!(body, getentry(ldu, body.id), scale, mechanism)
        end
        for eqc in eqcs
            lineStep!(eqc, getentry(ldu, eqc.id), scale, mechanism)
        end
        for ineqc in ineqcs
            lineStep!(ineqc, getineqentry(ldu, ineqc.id), scale, mechanism)
        end

        meritf1 = meritf(mechanism)
        if meritf1 >= meritf0
            scale += 1
        else
            return meritf1
        end
    end

    warning && (@info string("lineSearch! did not converge. n = ", iter, ". Last tol: ", meritf1))
    return meritf1
end

@inline function lineStep!(component::Component, diagonal::DiagonalEntry, scale)
    component.s1 = component.s0 - 1 / (2^scale) * diagonal.Δs
    return
end

@inline function lineStep!(node::Component, diagonal, scale, mechanism)
    node.s1 = node.s0 - 1 / (2^scale) * mechanism.α * diagonal.Δs
    return
end

@inline function lineStep!(ineqc::InequalityConstraint, entry, scale, mechanism)
    ineqc.s1 = ineqc.s0 - 1 / (2^scale) * mechanism.α * entry.Δs
    ineqc.γ1 = ineqc.γ0 - 1 / (2^scale) * mechanism.α * entry.Δγ
    return
end


text:
#Pkg.add("Ipopt");Pkg.build("Ipopt")
Pkg.add("PyPlot");Pkg.build("PyPlot")
#Pkg.add("PGFPlots");Pkg.build("PGFPlots")
#Pkg.add("GR");Pkg.build("GR")

using Documenter,MPCDocs,NLOptControl,PrettyPlots
makedocs(modules=[NLOptControl,PrettyPlots],
        doctest=false, clean=true,
        format =:html,
        authors="Huckleberry Febbo",
        sitename="NLOptControl.jl",
        pages = Any[
        "Home" => "index.md",
        "Tutorials"=>Any[
              "tutorials/Brachistochrone/main.md",
              "tutorials/MoonLander/main.md",
              "tutorials/BrysonDenham/main.md",
              "tutorials/Beam/main.md",
              "tutorials/HyperSensitive/main.md",
              "tutorials/RobotArm/main.md",
              "tutorials/Rocket/main.md",
              "tutorials/Unicycle/main.md",
              "tutorials/KinematicBicycle/main.md"
               ],
       "MPC"=>Any[
           "mpc/index.md"
           ]
               ]
               )

deploydocs(
    deps=Deps.pip("mkdocs","python-markdown-math"),
    repo="github.com/JuliaMPC/MPCDocs.jl.git",
    target="build",
    osname="linux",
    julia="0.6",
    make=nothing)

#=

     "Exported Functions"=>Any[
     "functions/NLOptControl.md",
     "functions/PrettyPlots.md"
     ]
     ]=#


text:
for Op in (:sqrt, :exp, :log,
           :sin, :cos, :tan, :csc, :sec, :cot,
           :asin, :acos, :atan, :acsc, :asec, :acot,
           :sinh, :cosh, :tanh, :csch, :sech, :coth,
           :asinh, :acosh, :atanh, :acsch, :asech, :acoth)
    @eval $Op(x::Single32) = Single32($Op(Float64(x)))
end


text:
using StrTables, Emoji_Entities

@static VERSION < v"0.7.0-DEV" ? (using Base.Test) : (using Test)

# Test the functions lookupname, matches, longestmatches, completions
# Check that characters from all 3 tables (BMP, non-BMP, string) are tested

const def = Emoji_Entities.default

@testset "Emoji_Entities" begin
@testset "lookupname" begin
    @test lookupname(def, SubString("My name is Spock", 12)) == ""
    @test lookupname(def, "foobar")   == ""
    @test lookupname(def, "sailboat") == "\u26f5"
    @test lookupname(def, "ring")     == "\U1f48d"
    @test lookupname(def, "flag-us")  == "\U1f1fa\U1f1f8"
end

@testset "matches" begin
    @test isempty(matches(def, ""))
    @test isempty(matches(def, "\u2020"))
    @test isempty(matches(def, SubString("This is \u2020", 9)))
    for (chrs, exp) in (("\u26f5", ["boat", "sailboat"]),
                        ("\U1f48d", ["ring"]),
                        ("\U1f596", ["spock-hand"]),
                        ("\U1f1fa\U1f1f8", ["flag-us", "us"]))
        res = matches(def, chrs)
        @test length(res) >= length(exp)
        @test intersect(res, exp) == exp
    end
end

@testset "longestmatches" begin
    @test isempty(longestmatches(def, "\u2020 abcd"))
    @test isempty(longestmatches(def, SubString("This is \U2020", 9)))
    for (chrs, exp) in (("\u26f5 abcd", ["boat", "sailboat"]),
                        ("\U1f48d abcd", ["ring"]),
                        ("\U1f1fa\U1f1f8 foo", ["flag-us", "us"]))
        res = longestmatches(def, chrs)
        @test length(res) >= length(exp)
        @test intersect(res, exp) == exp
    end
end

@testset "completions" begin
    @test isempty(completions(def, "ScottPaulJones"))
    @test isempty(completions(def, SubString("My name is Scott", 12)))
    for (chrs, exp) in (("al", ["alarm_clock", "alembic", "alien"]),
                        ("um", ["umbrella", "umbrella_on_ground", "umbrella_with_rain_drops"]))
        res = completions(def, chrs)
        @test length(res) >= length(exp)
        @test intersect(res, exp) == exp
    end
end
end


text:
# mri/z-list.jl

include("kspace.jl")
export ir_mri_kspace_ga_radial

include("sensemap-sim.jl")
export ir_mri_sensemap_sim


text:
include("testGetData2D.jl")
include("testGetData3D.jl")
include("testSensMatVec.jl")
include("testTravelTimeInversion.jl")


text:
module Rotations

using Test

using Gridap
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Arrays
using Gridap.Helpers
using STLCutters

using STLCutters: read_stl
using STLCutters: compute_stl_model
using STLCutters: merge_nodes
using STLCutters: get_bounding_box
using STLCutters: compute_grid
using STLCutters: surface, volume, surfaces, volumes
using STLCutters:  FACE_IN, FACE_OUT, FACE_CUT

function test_stl_cut(model,stl,vol)
  subcells,subfaces,labels = subtriangulate(model,stl,surfacesource=:both)

  grid = get_grid(model)

  bgmesh_in_vol, bgmesh_out_vol, bgmesh_cut_vol = volumes(grid,labels.bgcell_to_ioc)
  submesh_in_vol,submesh_out_vol, = volumes(subcells,labels.cell_to_io)
  in_surf,out_surf,skin_surf = surfaces(subfaces,labels.face_to_ios)

  in_volume = bgmesh_in_vol + submesh_in_vol
  out_volume = bgmesh_out_vol + submesh_out_vol
  cut_volume = bgmesh_cut_vol

  stl_surf = surface(get_grid(stl))

#  celldata = [ "inoutcut" => labels.bgcell_to_ioc ]
#  writevtk( grid, "bgcells"; celldata )
#  celldata = [ "inout" => labels.cell_to_io, "bgcell" => labels.cell_to_bgcell ]
#  writevtk( subcells, "subcells"; celldata )
#  celldata = [ "bgcell" => labels.face_to_bgcell, "inoutskin" => labels.face_to_ios ]
#  writevtk( subfaces, "subfaces"; celldata )
#  writevtk( stl, "stl" )

  println("\t εV = $(in_volume + out_volume - volume(grid))")
  println("\t εVin = $(in_volume-vol)")
  println("\t εΓin = $(in_surf - stl_surf))")
  println("\t εΓout = $(out_surf - stl_surf))")
  println("\t εΓskin = $(skin_surf - stl_surf))")

  tol = 1e-10
  @test abs( submesh_in_vol + submesh_out_vol - cut_volume ) / cut_volume < tol
  @test abs( in_surf - stl_surf ) / stl_surf < tol
  @test abs( out_surf - stl_surf ) / stl_surf < tol
  @test abs( skin_surf - stl_surf ) / stl_surf < tol
  @test abs( in_volume + out_volume - volume(grid) ) / volume(grid) < tol
  @test abs( in_volume - vol ) / vol < tol
end

Rx(ϕ) = TensorValue(
  1,0,0,
  0,cos(ϕ),-sin(ϕ),
  0,sin(ϕ),cos(ϕ))

Ry(θ) = TensorValue(
  cos(θ),0,sin(θ),
  0,1,0,
  -sin(θ),0,cos(θ))

Rz(ψ) = TensorValue(
  cos(ψ),sin(ψ),0,
  -sin(ψ),cos(ψ),0,
  0,0,1)

R(ϕ,θ,ψ) = Rx(ϕ)⋅Ry(θ)⋅Rz(ψ)

R(θ) = R(θ,θ,θ)

#X,T,N = read_stl(joinpath(@__DIR__,"data/Bunny-LowPoly.stl"))
#X,T,N = read_stl(joinpath(@__DIR__,"data/wine_glass.stl"))
X,T,N = read_stl(joinpath(@__DIR__,"data/cube.stl"))

stl0 = compute_stl_model(T,X)
stl0 = merge_nodes(stl0)

X0 = get_node_coordinates(get_grid(stl0))
T0 = get_cell_node_ids(stl0)

## Coindident Grid-Cube
δ = 0.2
n = 7
D = 3

pmin,pmax = get_bounding_box(stl0)
O = (pmin+pmax)/2

Δ = (pmax-pmin)*δ
pmin = pmin - Δ
pmax = pmax + Δ
partition = (n,n,n)

model = CartesianDiscreteModel(pmin,pmax,partition)

θs = [0;exp10.( -17:0 )]
Δx = VectorValue(0.0,0.0,0.1)

for θ in θs
  println("Testing θz = $θ ...")
  Oi = O + Δx
  Xi = map(p-> p + Δx,X0)
  Xi = map(p-> Oi + Rz(θ)⋅(p-Oi),Xi)
  stl = compute_stl_model(Table(T0),Xi)
  #writevtk(stl.grid,"stl")
  test_stl_cut(model,stl,1)
end

## Rational Grid
δ = 0.2
n = 7
D = 3

pmin,pmax = get_bounding_box(stl0)
pmin = Point(rationalize.(Tuple(pmin)))
pmax = Point(rationalize.(Tuple(pmax)))
δ = rationalize(δ)

O = (pmin+pmax)/2

Δ = (pmax-pmin)*δ
pmin = pmin - Δ
pmax = pmax + Δ
partition = (n,n,n)

model = CartesianDiscreteModel(pmin,pmax,partition)

δx = Point(0,0,1)
Δxs = [0;exp10.(-17:-1)]

for Δx in Δxs
  println("Testing Δx = $Δx ...")
  Xi = map(p-> p + Δx,X0)
  stl = compute_stl_model(Table(T0),Xi)
  #writevtk(stl.grid,"stl")
  test_stl_cut(model,stl,1)
end

δ = 0.2
n = 7
D = 3

pmin,pmax = get_bounding_box(stl0)
pmin = Point(rationalize.(Tuple(pmin)))
pmax = Point(rationalize.(Tuple(pmax)))
δ = rationalize(δ)

origins = [pmin, pmax, (pmin+pmax)/2 ]

Δ = (pmax-pmin)*δ
pmin = pmin - Δ
pmax = pmax + Δ
partition = (n,n,n)

model = CartesianDiscreteModel(pmin,pmax,partition)

θs = [0;exp10.( -17:-1 )]

for O in origins
  for θ in θs
    println("Testing θ = $θ over $O ...")
    Xi = map(p-> O + R(θ)⋅(p-O),X0)
    stl = compute_stl_model(Table(T0),Xi)
    #writevtk(stl,"stl")
    test_stl_cut(model,stl,1)
  end
end

end # module


text:
# ----------------------------------------------------------
# ---------------  Auxiliary functionalities  --------------
# ----------------------------------------------------------

"""
Manually set the plan of a FunctionOperator
Arguments:
 - `FO` composite FunctionOperator to be changed
 - `f` manually defined plan (function with two arguments, first is the output buffer, second is the input array)
 - `f_str` (Optional) string representation of the plan
"""
function setPlan(FO::FunctionOperatorComposite, f::Function, f_str::String = "manually defined")
    @assert nargs(f) == 3  "plan must accept two arguments!"
    FO.plan_function = f
    FO.plan_string = f_str
    FO
end

function Base.reshape(FO::FunctionOperator; inDims::Union{Nothing, NTuple{N, Int}} where N = nothing,
        outDims::Union{Nothing, NTuple{N, Int}} where N = nothing)
    inDims = inDims isa Nothing ? FO.inDims : inDims
    outDims = outDims isa Nothing ? FO.outDims : outDims
    (iDims, oDims) = FO.adjoint ? (outDims, inDims) : (inDims, outDims)
    FunctionOperator(FO,
        name = "✠($(FO.name))",
        forw = checkTwoInputs(FO.forw) ?
            (b,x) -> begin FO.forw(reshape(b, FO.outDims), reshape(x, FO.inDims)); b; end :
             x -> reshape(FO.forw(reshape(x, FO.inDims)), oDims),
        backw = checkTwoInputs(FO.backw) ?
            (b,x) -> begin FO.backw(reshape(b, FO.inDims), reshape(x, FO.outDims)); b; end :
             x -> reshape(FO.backw(reshape(x, FO.outDims)), iDims),
        inDims = inDims, outDims = outDims)
end

function Base.reshape(FO::FunctionOperatorComposite; inDims::Union{Nothing, NTuple{N, Int}} where N = nothing,
        outDims::Union{Nothing, NTuple{N, Int}} where N = nothing)
    FOᴴ = FO'
    inDims = inDims isa Nothing ? FO.inDims : inDims
    outDims = outDims isa Nothing ? FO.outDims : outDims
    FunctionOperator{eltype(FO)}(name = "✠($(FO.name))",
        forw = (b,x) -> begin mul!(reshape(b, FO.outDims), FO, reshape(x, FO.inDims)); b; end,
        backw = (b,x) -> begin mul!(reshape(b, FOᴴ.outDims), FOᴴ, reshape(x, FOᴴ.inDims)); b; end,
        inDims = inDims, outDims = outDims)
end

# --------- Equality check -----------

# Expand arithmetic expression trees
function expand(expr)
    sortExprList(ExprList) = sort(ExprList, by = x -> string(x), rev = true)
    adj(e) = (e isa Expr && string(e.head) == "'") ? e.args[1] : Expr(Symbol("'"), e)
    MacroTools.postwalk(x -> begin
        # Adjoint expansion. Eg: (op1 * op2)' -> op2' * op1'
        if @capture(x, (e1_ * e2__)')
            Expr(:call, :*, [adj(e) for e in reverse(e2)]..., adj(e1))
        # Multiplication associativity
        elseif @capture(x, (e1_ * e2__) * e3__)
            Expr(:call, :*, e1, e2..., e3...)
        elseif @capture(x, e1__ * (e2_ * e3__))
            Expr(:call, :*, e1..., e2, e3...)
        # Remove unnecessary parentheses around addition by associativity.
        # Eg: op1 + (op2 + op3) -> op1 + op2 + op3
        elseif @capture(x, e1__ + e2_)
            list = [e1..., e2]
            new_list = []
            for item in list
                if @capture(item, E1__ + E2_)
                    new_list = vcat(new_list, E1, E2)
                else
                    push!(new_list, item)
                end
            end
            Expr(:call, :+, sortExprList(new_list)...)
        # Expandion by distributivity. Eg: (op1 ± op2) * op3 -> op1 * op3 ± op2 * op3
        elseif @capture(x, (e1__ + e2_) * (e3_))
            list = sortExprList([e1..., e2])
            Expr(:call, :+, sortExprList([Expr(:call, :*, e, e3) for e in list])...)
        elseif @capture(x, (e1__ + e2_) * e3__)
            list = sortExprList([e1..., e2])
            Expr(:call, :+, sortExprList([Expr(:call, :*, e, e3...) for e in list])...)
        elseif @capture(x, e0__ * (e1__ + e2_) * (e3_))
            list = sortExprList([e1..., e2])
            Expr(:call, :*, e0...,
                    Expr(:call, :+, sortExprList([Expr(:call, :*, e, e3) for e in list])...))
        elseif @capture(x, e0_ * (e1__ + e2_) * e3__) && length(e3) > 0
            list = sortExprList([e1..., e2])
            Expr(:call, :*, e0,
                Expr(:call, :+, sortExprList([Expr(:call, :*, e, e3...) for e in list])...))
        # Substraction as negation
        elseif @capture(x, -(e1_ * e2__))
            Expr(:call, :*, Expr(:call, :-, e1), e2...)
        elseif @capture(x, -(e1_ + e2__))
            Expr(:call, :+, Expr(:call, :-, e1), [Expr(:call, :-, e) for e in e2]...)
        elseif @capture(x, -(e1_ - e2_))
            Expr(:call, :+, Expr(:call, :-, e1), e2)
        elseif @capture(x, - -e2_)
            e2
        elseif @capture(x, e1_ - e2_)
            Expr(:call, :+, e1, Expr(:call, :-, e2))
        else
            x
        end
    end, expr)
end

# Normalize (expand fully) an arithmetic expression given as a string
function normalizeExpression(str)
    expr = Meta.parse(str) # I use the Julia parser to parse my simple arithmetical expressions
    # Not too elegant solution, I know...
    # But even complicated expressions can be fully expanded within a couple of iterations
    for i = 1:30
        new_expr = expand(expr)
        string(new_expr) == string(expr) && break
        expr = new_expr
    end
    string(expr)
end

"""
Equality check based on name

**HEADS UP!** This works only when operators with the same name has also the same functionality!

It performs basic arithmetic transformations on the expressions, so it recognizes even some less obvious equalities. The rules it uses:
 - Associativity: ``op1 * op2 * op3 = (op1 * op2) * op3 = op1 * (op2 * op3)``, ``op1 + op2 + op3 = (op1 + op2) + op3 = op1 + (op2 + op3)``, ``op1 + (op2 - op3) == (op1 + op2) - op3``, ``op1 - (op2 + op3) == (op1 - op2) - op3``
 - Commutativity: ``op1 - op2 - op3 = op1 - op3 - op2``, ``op1 + op2 = op2 + op1``
 - Distributivity: ``(op1 + op2) * op3 = op1 * op3 + op2 * op3`` (Note that ``op1 * (op2 + op3) ≠ op1 * op2 + op1 * op3``)
"""
Base.:(==)(FM1::FunctionOperatorComposite, FM2::FunctionOperatorComposite) =
    normalizeExpression(FM1.name) == normalizeExpression(FM2.name)

# --------- Recycling macro -----------
# Preallocates arrays for the output of different operations

"""
**Recycling macro**: Reduce the number of allocations inside a for loop by preallocation of arrays for the outputs of marked operations. Markers: `@♻` (`\\:recycle:`), `🔝` (`\\:top:`), `🔃` (`\\:arrows_clockwise:`), and `@🔃`

Macro @♻ should be placed right before a for loop, and then it executes the following substitutions:
 - **Expressions marked by `🔝`:**
They are going to be calculated before the loop, the result is stored in a variable, and the expression will be replaced by that variable. It also can be useful when a constant expression is used in the loop, but the idea behind creating that substitution is to allow caching of composite FunctionMatrices. Eg:
```julia
@♻ for i=1:5
    result = 🔝((FuncOp₁ + 2I) * FuncOp₂) * data
end
```
will be transformed to
```julia
🔝_1 = (FuncOp₁ + 2I) * FuncOp₂
for i = 1:5
    result = 🔝_1 * data
end
```
so that way plan is calculated only once, and also buffers for intermediate results of the composite operator are allocated once.

 - **Expressions marked by `🔃`:**
They are going to be calculated before the loop (to allocate an array to store the result), but the expression is also evaluated in each loop iteration. The difference after the substitution is that the result of the expression is always saved to the preallocated array. Eg:
```julia
@♻ for i=1:5
    result = FuncOp₁ * 🔃(A + B)
end
```
will be transformed to
```julia
🔃_1 = A + B
for i = 1:5
    result = FuncOp₁ * @.(🔃_1 = A + B)
end
```
This transformation first allocates an array named `🔃_1`, and then in every iteration it is recalculated, saved to `🔃_1`, and the this value is used for the rest of the operation (i.e.: `FuncOp₁ * 🔃_1`. Note that `@.` macro is inserted before the inline assignment. This is needed otherwise `A + B` would allocate a new array before it is stored in `🔃_1`. **Warning!** It can break your code, e.g. `@.(🔃_1 = A * B) ≠ (🔃_1 = A .* B)` {matrix multiplication vs. elementwise multiplication}! On the other hand, when the marked expression consists only a multiplication, then it is transformed into a call of `mul!`. Eg:
```julia
@♻ for i=1:5
    result = FuncOp₁ * 🔃(A * B)
end
```
will be transformed to
```julia
🔃_1 = A * B
for i = 1:5
    result = FuncOp₁ * mul!(🔃_1, A, B)
end
```

 - **Lastly, assignments marked by `@🔃`:**
They will be transformed into a call of `mul!`. Of course, it works only if `@🔃` is directly followed by an assignment that has a single multiplication on the right side. Eg:
```julia
@♻ for i=1:5
    @🔃 result = FuncOp₁ * A
end
```
will be transformed to
```julia
result = FuncOp₁ * A
for i = 1:5
    mul!(result, FuncOp₁, A)
end
```

Final note: `🔝` can be arbitrarily nested, and it can be embedded in expressions marked by `🔃`. `🔃` can also be nested, and it can be used in assigments marked by `@🔃` (along with `🔝`, of course).
"""
macro ♻(loop)
    🔝s = Array{Tuple{Symbol, Expr}}(undef, 0)
    🔃s = Array{Tuple{Symbol, Expr}}(undef, 0)
    loop = MacroTools.postwalk(x -> begin
        if @capture(x, 🔝(expr_))
            newSymbol = Symbol(:🔝_, length(🔝s) + 1)
            push!(🔝s, (newSymbol, expr))           # preallocate an array before the loop
            newSymbol                              # replace the expression with name of buffer
        elseif @capture(x, 🔃(expr_))
            newSymbol = Symbol(:🔃_, length(🔃s) + 1)
            push!(🔃s, (newSymbol, expr)) # preallocate an array before the loop
            if @capture(expr, lhs_ * rhs_)
                Expr(:call, :mul!, newSymbol, lhs, rhs)
            else
                transformed = MacroTools.postwalk(x -> begin
                    if @capture(x, lhs_ + rhs__)
                        Expr(:call, :(.+), lhs, rhs...)
                    elseif @capture(x, lhs_ - rhs_)
                        Expr(:call, :(.-), lhs, rhs)
                    elseif @capture(x, lhs_ * rhs_)
                        newSymbol = Symbol(:🔃_, length(🔃s) + 1)
                        push!(🔃s, (newSymbol, x)) # preallocate an array before the loop
                        Expr(:call, :mul!, newSymbol, lhs, rhs)
                    else
                        x
                    end
                end, expr)
                :(($newSymbol .= $transformed))
            end
        elseif @capture(x, @🔃 expr_)
            if @capture(expr, lhs_ = rhs1_ * rhs2_) || @capture(expr, lhs_ = mul!(e_, rhs1_, rhs2_))
                push!(🔃s, (lhs, :($rhs1 * $rhs2))) # preallocate an array before the loop
            elseif @capture(expr, lhs_ .= rhs1_ * rhs2_) || @capture(expr, lhs_ .= mul!(e_, rhs1_, rhs2_))
                # nothing to do here
            else
                @assert false "Macro @🔃 must be followed by an assigment that has a single multiplication on the right side"
            end
            :(mul!($lhs, $(rhs1), $(rhs2)))        # replace with mul!
        else
            x
        end
    end, loop)
    🔝_defs = [:($left = $right) for (left, right) in 🔝s]
    🔃_defs = [:($left = $right) for (left, right) in 🔃s]
    extended_loop = Expr(:block, 🔝_defs..., 🔃_defs..., loop)
    FunctionOperators_global_settings.macro_verbose && println(MacroTools.prewalk(rmlines, extended_loop))
    esc(extended_loop)
end


text:
module EvoDynamics

using Agents
using Random
using Distributions
using StatsBase
using LinearAlgebra
using StaticArrays
using YAML

include("simulation.jl")
include("interactions.jl")
include("migration.jl")
include("reproduction.jl")
include("api.jl")
include("load_params.jl")
include("check_params.jl")

end # module


text:
"""
    Truth Sensor

This sensor can detect all physical effects generated on the vehicle.

"""
function TruthSensor()
    # function init(t, constants, state, draws, inputs, effects, effects_bus)
    #     return state
    # end
    function sense(t, constants, state, draws, inputs, effects, effects_bus)
        body, = find_effect(effects, BodyStateEffect())
        return (state, body, true)
    end
    DynamicalModel("truth_sensor",
                   init = sense,
                   update = sense,
                   outputs = BodyStateEffect(),
                   timing = ModelTiming(0.02))
end

"""
    Star Tracker

This generic star tracker always reports the body's orientation (a quaternion)
with a little noise.

"""
function StarTracker()
    function sense(t, constants, state, draws, inputs, effects, effects_bus)
        errors = constants .* draws
        body, = find_effect(effects, BodyStateEffect())
        return (state, qcomp(mrp2q(errors, 4.), body.q_BI), true)
    end
    DynamicalModel("default_star_tracker",
                   init = sense,
                   update = sense,
                   timing = ModelTiming(0.1), # 10Hz update rate
                   constants = [0.001, 0.001, 0.001], # constants: error magnitudes (rad)
                   outputs = [0.; 0.; 0.; 1.],
                   rand = RandSpec(init=3,update=3)) # rand
end

"""
    Gyro

This is a generic rate-integrating gyro with angular random walk due to noise
and a bias random walk. It doesn't use the state's rotation rate, but actually
differences the quaternion.

"""
mutable struct GyroConstants
    angular_random_walk::Float64
    bias_random_walk::Float64
    initial_bias_magnitude::Float64
    # q_SB::Vector{Float64} # Orientation wrt the body
end
mutable struct GyroState
    bias::Vector{Float64}
    q_SI::Vector{Float64}
    t::Float64
end

function Gyro()
    function init(t, constants, state, draws, inputs, effects, effects_bus)
        initial_bias = constants.initial_bias_magnitude * (2. .* draws[1:3] .- 1.)
        body, = find_effect(effects, BodyStateEffect())
        return (GyroState(initial_bias, body.q_BI, t),
                body.ω_BI_B, # Outputs
                true) # Active?
    end
    function update(t, constants, state, draws, inputs, effects, effects_bus)
        body, found  = find_effect(effects, BodyStateEffect())
        if !found; error("The Gyro model requires a BodyState effect but could not find one."); end
        Δt     = t - state.t # It would be nice if models had access to their own dt.
        σ_rate = Δt > 0. ? constants.angular_random_walk / √(Δt) : 0. # These could be saved to constants for a fixed time step, which would save us from dividing by zero at the first sample.
        σ_bias = constants.bias_random_walk * √(Δt)
        noise  = σ_rate * draws[1:3]
        Δb     = σ_bias * draws[4:6]
        θ, r   = q2aa(qdiff(body.q_BI, state.q_SI))
        Δθ     = θ * r + (state.bias + 0.5*Δb) * Δt + noise
        ω      = Δθ / Δt
        return (GyroState(state.bias + Δb, body.q_BI, t), # Updated state
                ω,                                        # Measurement
                true)                                     # Active?
    end
    DynamicalModel("gyro",
                   init = init,
                   update = update,
                   timing = ModelTiming(0.01),
                   constants = GyroConstants(0.00001, 0.0001, 0.001),
                   state = GyroState([0.; 0.; 0.], [0.; 0.; 0.; 1.], 0.),
                   outputs = [0.; 0.; 0.],
                   rand = RandSpec(init=(rand,3), update=6))
end


text:
# Test conversions ta various normal forms
@test to_nnf(@julog(and(a))) == @julog a
@test to_nnf(@julog(true => not(and(not(!a), b, or(not(c), false))))) ==
    @julog(or(not(a), not(b), c))

@test to_cnf(@julog(and(a))) == @julog and(or(a))
@test to_cnf(@julog(and(and(or(a, and(b, or(c, d))), or(e, f)), and(not(x), or(y, z))))) ==
    @julog(and(or(a, b), or(a, c, d), or(e, f), or(not(x)), or(y, z)))
@test to_cnf(@julog(true => not(and(not(!a), b, or(not(c), false))))) ==
    @julog(and(or(not(a), not(b), c)))

@test to_dnf(@julog(and(a))) == @julog or(and(a))
@test to_dnf(@julog(or(or(and(a, or(b, and(c, d))), and(e, f)), or(not(x), and(y, z))))) ==
    @julog(or(and(a, b), and(a, c, d), and(e, f), and(not(x)), and(y, z)))
@test to_dnf(@julog(true => not(and(not(!a), b, or(not(c), false))))) ==
    @julog(or(and(not(a)), and(not(b)), and(c)))

# Test flattening of conjuctions and disjunctions
@test flatten_conjs(@julog(and(a, and(b, c)))) == @julog Const[a, b, c]
@test flatten_conjs(@julog([a, and(b, c)])) == @julog Const[a, b, c]
@test flatten_disjs(@julog(or(a, or(b, c)))) == @julog Const[a, b, c]
@test flatten_disjs(@julog([a, or(b, c)])) == @julog Const[a, b, c]

# Test instantiation of universal quantifiers
clauses = @julog [
    block(a) <<= true,
    block(b) <<= true,
    block(c) <<= true,
    holding(a) <<= true
]
universal_term = @julog(forall(block(X), not(holding(X))))
@test deuniversalize(universal_term, clauses) ==
    @julog and(not(holding(a)), not(holding(b)), not(holding(c)))
universal_clause = Clause(@julog(handempty), [universal_term])
@test deuniversalize(universal_clause, clauses) ==
    @julog handempty <<= and(not(holding(a)), not(holding(b)), not(holding(c)))

# Test regularization of clause bodies
clauses = @julog [
    binary(X) <<= or(woman(X), man(X)) & not(nonbinary(X)),
    bigender(X) <<= and(woman(X), man(X)),
    nonbinary(X) <<= or(genderqueer(X), agender(X), thirdgender(X))
]

regularized = @julog [
    binary(X) <<= woman(X) & not(nonbinary(X)),
    binary(X) <<= man(X) & not(nonbinary(X)),
    bigender(X) <<= woman(X) & man(X),
    nonbinary(X) <<= genderqueer(X),
    nonbinary(X) <<= agender(X),
    nonbinary(X) <<= thirdgender(X)
]

@test regularize_clauses(clauses) == regularized


text:
using FlightSims
const FS = FlightSims

using Plots
using Transducers
using Test

function get_traj_data(env, Δt, tf)
    x0 = State(env)()
    simulator = Simulator(x0, Dynamics!(env); tf=tf)
    df = solve(simulator; savestep=Δt)
    df
end

function plot_figures!(fig, df, pos_cmd_func)
    ts = df.time
    poss = df.sol |> Map(datum -> datum.multicopter.state.p) |> collect
    poss_cmd = ts |> Map(pos_cmd_func) |> collect
    plot!(fig,
          ts, hcat(poss_cmd...)';
          label="desired",
          color=:red,
         )
    plot!(fig,
          ts, hcat(poss...)';
          label="true",
          color=:blue,
         )
end

function main()
    pos_cmd_func = (t) -> [sin(t), cos(t), t]
    env = FS.BacksteppingPositionController_StaticAllocator_Multicopter(pos_cmd_func)
    Δt = 0.01
    tf = 10.0
    df = get_traj_data(env, Δt, tf)
    # plot
    fig = plot(; legend=:topleft)
    plot_figures!(fig, df, pos_cmd_func)
    savefig("figures/multicopter_position_control.png")
    display(fig)
end

@testset "hexacopter position control example" begin
    main()
end


text:
# Julia wrapper for header: libssh2.h
# Automatically generated using Clang.jl

function libssh2_init(flags)
    ccall((:libssh2_init, libssh2), Cint, (Cint,), flags)
end

function libssh2_exit()
    ccall((:libssh2_exit, libssh2), Cvoid, ())
end

function libssh2_free(session, ptr)
    ccall((:libssh2_free, libssh2), Cvoid, (Ptr{LIBSSH2_SESSION}, Ptr{Cvoid}), session, ptr)
end

function libssh2_session_supported_algs(session, method_type, algs)
    ccall((:libssh2_session_supported_algs, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cint, Ptr{Ptr{Cstring}}), session, method_type, algs)
end

function libssh2_session_init_ex(my_alloc, my_free, my_realloc, abstract)
    ccall((:libssh2_session_init_ex, libssh2), Ptr{LIBSSH2_SESSION}, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), my_alloc, my_free, my_realloc, abstract)
end

function libssh2_session_abstract(session)
    ccall((:libssh2_session_abstract, libssh2), Ptr{Ptr{Cvoid}}, (Ptr{LIBSSH2_SESSION},), session)
end

function libssh2_session_callback_set(session, cbtype, callback)
    ccall((:libssh2_session_callback_set, libssh2), Ptr{Cvoid}, (Ptr{LIBSSH2_SESSION}, Cint, Ptr{Cvoid}), session, cbtype, callback)
end

function libssh2_session_banner_set(session, banner)
    ccall((:libssh2_session_banner_set, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cstring), session, banner)
end

function libssh2_banner_set(session, banner)
    ccall((:libssh2_banner_set, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cstring), session, banner)
end

function libssh2_session_startup(session, sock)
    ccall((:libssh2_session_startup, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cint), session, sock)
end

function libssh2_session_handshake(session, sock)
    ccall((:libssh2_session_handshake, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, libssh2_socket_t), session, sock)
end

function libssh2_session_disconnect_ex(session, reason, description, lang)
    ccall((:libssh2_session_disconnect_ex, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cint, Cstring, Cstring), session, reason, description, lang)
end

function libssh2_session_free(session)
    ccall((:libssh2_session_free, libssh2), Cint, (Ptr{LIBSSH2_SESSION},), session)
end

function libssh2_hostkey_hash(session, hash_type)
    ccall((:libssh2_hostkey_hash, libssh2), Cstring, (Ptr{LIBSSH2_SESSION}, Cint), session, hash_type)
end

function libssh2_session_hostkey(session, len, type)
    ccall((:libssh2_session_hostkey, libssh2), Cstring, (Ptr{LIBSSH2_SESSION}, Ptr{Csize_t}, Ptr{Cint}), session, len, type)
end

function libssh2_session_method_pref(session, method_type, prefs)
    ccall((:libssh2_session_method_pref, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cint, Cstring), session, method_type, prefs)
end

function libssh2_session_methods(session, method_type)
    ccall((:libssh2_session_methods, libssh2), Cstring, (Ptr{LIBSSH2_SESSION}, Cint), session, method_type)
end

function libssh2_session_last_error(session, errmsg, errmsg_len, want_buf)
    ccall((:libssh2_session_last_error, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Ptr{Cstring}, Ptr{Cint}, Cint), session, errmsg, errmsg_len, want_buf)
end

function libssh2_session_last_errno(session)
    ccall((:libssh2_session_last_errno, libssh2), Cint, (Ptr{LIBSSH2_SESSION},), session)
end

function libssh2_session_set_last_error(session, errcode, errmsg)
    ccall((:libssh2_session_set_last_error, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cint, Cstring), session, errcode, errmsg)
end

function libssh2_session_block_directions(session)
    ccall((:libssh2_session_block_directions, libssh2), Cint, (Ptr{LIBSSH2_SESSION},), session)
end

function libssh2_session_flag(session, flag, value)
    ccall((:libssh2_session_flag, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cint, Cint), session, flag, value)
end

function libssh2_session_banner_get(session)
    ccall((:libssh2_session_banner_get, libssh2), Cstring, (Ptr{LIBSSH2_SESSION},), session)
end

function libssh2_userauth_list(session, username, username_len)
    ccall((:libssh2_userauth_list, libssh2), Cstring, (Ptr{LIBSSH2_SESSION}, Cstring, UInt32), session, username, username_len)
end

function libssh2_userauth_authenticated(session)
    ccall((:libssh2_userauth_authenticated, libssh2), Cint, (Ptr{LIBSSH2_SESSION},), session)
end

function libssh2_userauth_password_ex(session, username, username_len, password, password_len, passwd_change_cb)
    ccall((:libssh2_userauth_password_ex, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cstring, UInt32, Cstring, UInt32, Ptr{Cvoid}), session, username, username_len, password, password_len, passwd_change_cb)
end

function libssh2_userauth_publickey_fromfile_ex(session, username, username_len, publickey, privatekey, passphrase)
    ccall((:libssh2_userauth_publickey_fromfile_ex, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cstring, UInt32, Cstring, Cstring, Cstring), session, username, username_len, publickey, privatekey, passphrase)
end

function libssh2_userauth_publickey(session, username, pubkeydata, pubkeydata_len, sign_callback, abstract)
    ccall((:libssh2_userauth_publickey, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cstring, Ptr{Cuchar}, Csize_t, Ptr{Cvoid}, Ptr{Ptr{Cvoid}}), session, username, pubkeydata, pubkeydata_len, sign_callback, abstract)
end

function libssh2_userauth_hostbased_fromfile_ex(session, username, username_len, publickey, privatekey, passphrase, hostname, hostname_len, local_username, local_username_len)
    ccall((:libssh2_userauth_hostbased_fromfile_ex, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cstring, UInt32, Cstring, Cstring, Cstring, Cstring, UInt32, Cstring, UInt32), session, username, username_len, publickey, privatekey, passphrase, hostname, hostname_len, local_username, local_username_len)
end

function libssh2_userauth_publickey_frommemory(session, username, username_len, publickeyfiledata, publickeyfiledata_len, privatekeyfiledata, privatekeyfiledata_len, passphrase)
    ccall((:libssh2_userauth_publickey_frommemory, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cstring, Csize_t, Cstring, Csize_t, Cstring, Csize_t, Cstring), session, username, username_len, publickeyfiledata, publickeyfiledata_len, privatekeyfiledata, privatekeyfiledata_len, passphrase)
end

function libssh2_userauth_keyboard_interactive_ex(session, username, username_len, response_callback)
    ccall((:libssh2_userauth_keyboard_interactive_ex, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cstring, UInt32, Ptr{Cvoid}), session, username, username_len, response_callback)
end

function libssh2_poll(fds, nfds, timeout)
    ccall((:libssh2_poll, libssh2), Cint, (Ptr{LIBSSH2_POLLFD}, UInt32, Clong), fds, nfds, timeout)
end

function libssh2_channel_open_ex(session, channel_type, channel_type_len, window_size, packet_size, message, message_len)
    ccall((:libssh2_channel_open_ex, libssh2), Ptr{LIBSSH2_CHANNEL}, (Ptr{LIBSSH2_SESSION}, Cstring, UInt32, UInt32, UInt32, Cstring, UInt32), session, channel_type, channel_type_len, window_size, packet_size, message, message_len)
end

function libssh2_channel_direct_tcpip_ex(session, host, port, shost, sport)
    ccall((:libssh2_channel_direct_tcpip_ex, libssh2), Ptr{LIBSSH2_CHANNEL}, (Ptr{LIBSSH2_SESSION}, Cstring, Cint, Cstring, Cint), session, host, port, shost, sport)
end

function libssh2_channel_forward_listen_ex(session, host, port, bound_port, queue_maxsize)
    ccall((:libssh2_channel_forward_listen_ex, libssh2), Ptr{LIBSSH2_LISTENER}, (Ptr{LIBSSH2_SESSION}, Cstring, Cint, Ptr{Cint}, Cint), session, host, port, bound_port, queue_maxsize)
end

function libssh2_channel_forward_cancel(listener)
    ccall((:libssh2_channel_forward_cancel, libssh2), Cint, (Ptr{LIBSSH2_LISTENER},), listener)
end

function libssh2_channel_forward_accept(listener)
    ccall((:libssh2_channel_forward_accept, libssh2), Ptr{LIBSSH2_CHANNEL}, (Ptr{LIBSSH2_LISTENER},), listener)
end

function libssh2_channel_setenv_ex(channel, varname, varname_len, value, value_len)
    ccall((:libssh2_channel_setenv_ex, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL}, Cstring, UInt32, Cstring, UInt32), channel, varname, varname_len, value, value_len)
end

function libssh2_channel_request_pty_ex(channel, term, term_len, modes, modes_len, width, height, width_px, height_px)
    ccall((:libssh2_channel_request_pty_ex, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL}, Cstring, UInt32, Cstring, UInt32, Cint, Cint, Cint, Cint), channel, term, term_len, modes, modes_len, width, height, width_px, height_px)
end

function libssh2_channel_request_pty_size_ex(channel, width, height, width_px, height_px)
    ccall((:libssh2_channel_request_pty_size_ex, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL}, Cint, Cint, Cint, Cint), channel, width, height, width_px, height_px)
end

function libssh2_channel_x11_req_ex(channel, single_connection, auth_proto, auth_cookie, screen_number)
    ccall((:libssh2_channel_x11_req_ex, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL}, Cint, Cstring, Cstring, Cint), channel, single_connection, auth_proto, auth_cookie, screen_number)
end

function libssh2_channel_process_startup(channel, request, request_len, message, message_len)
    ccall((:libssh2_channel_process_startup, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL}, Cstring, UInt32, Cstring, UInt32), channel, request, request_len, message, message_len)
end

function libssh2_channel_read_ex(channel, stream_id, buf, buflen)
    ccall((:libssh2_channel_read_ex, libssh2), ssize_t, (Ptr{LIBSSH2_CHANNEL}, Cint, Cstring, Csize_t), channel, stream_id, buf, buflen)
end

function libssh2_poll_channel_read(channel, extended)
    ccall((:libssh2_poll_channel_read, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL}, Cint), channel, extended)
end

function libssh2_channel_window_read_ex(channel, read_avail, window_size_initial)
    ccall((:libssh2_channel_window_read_ex, libssh2), Culong, (Ptr{LIBSSH2_CHANNEL}, Ptr{Culong}, Ptr{Culong}), channel, read_avail, window_size_initial)
end

function libssh2_channel_receive_window_adjust(channel, adjustment, force)
    ccall((:libssh2_channel_receive_window_adjust, libssh2), Culong, (Ptr{LIBSSH2_CHANNEL}, Culong, Cuchar), channel, adjustment, force)
end

function libssh2_channel_receive_window_adjust2(channel, adjustment, force, storewindow)
    ccall((:libssh2_channel_receive_window_adjust2, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL}, Culong, Cuchar, Ptr{UInt32}), channel, adjustment, force, storewindow)
end

function libssh2_channel_write_ex(channel, stream_id, buf, buflen)
    ccall((:libssh2_channel_write_ex, libssh2), ssize_t, (Ptr{LIBSSH2_CHANNEL}, Cint, Cstring, Csize_t), channel, stream_id, buf, buflen)
end

function libssh2_channel_window_write_ex(channel, window_size_initial)
    ccall((:libssh2_channel_window_write_ex, libssh2), Culong, (Ptr{LIBSSH2_CHANNEL}, Ptr{Culong}), channel, window_size_initial)
end

function libssh2_session_set_blocking(session, blocking)
    ccall((:libssh2_session_set_blocking, libssh2), Cvoid, (Ptr{LIBSSH2_SESSION}, Cint), session, blocking)
end

function libssh2_session_get_blocking(session)
    ccall((:libssh2_session_get_blocking, libssh2), Cint, (Ptr{LIBSSH2_SESSION},), session)
end

function libssh2_channel_set_blocking(channel, blocking)
    ccall((:libssh2_channel_set_blocking, libssh2), Cvoid, (Ptr{LIBSSH2_CHANNEL}, Cint), channel, blocking)
end

function libssh2_session_set_timeout(session, timeout)
    ccall((:libssh2_session_set_timeout, libssh2), Cvoid, (Ptr{LIBSSH2_SESSION}, Clong), session, timeout)
end

function libssh2_session_get_timeout(session)
    ccall((:libssh2_session_get_timeout, libssh2), Clong, (Ptr{LIBSSH2_SESSION},), session)
end

function libssh2_channel_handle_extended_data(channel, ignore_mode)
    ccall((:libssh2_channel_handle_extended_data, libssh2), Cvoid, (Ptr{LIBSSH2_CHANNEL}, Cint), channel, ignore_mode)
end

function libssh2_channel_handle_extended_data2(channel, ignore_mode)
    ccall((:libssh2_channel_handle_extended_data2, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL}, Cint), channel, ignore_mode)
end

function libssh2_channel_flush_ex(channel, streamid)
    ccall((:libssh2_channel_flush_ex, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL}, Cint), channel, streamid)
end

function libssh2_channel_get_exit_status(channel)
    ccall((:libssh2_channel_get_exit_status, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL},), channel)
end

function libssh2_channel_get_exit_signal(channel, exitsignal, exitsignal_len, errmsg, errmsg_len, langtag, langtag_len)
    ccall((:libssh2_channel_get_exit_signal, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL}, Ptr{Cstring}, Ptr{Csize_t}, Ptr{Cstring}, Ptr{Csize_t}, Ptr{Cstring}, Ptr{Csize_t}), channel, exitsignal, exitsignal_len, errmsg, errmsg_len, langtag, langtag_len)
end

function libssh2_channel_send_eof(channel)
    ccall((:libssh2_channel_send_eof, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL},), channel)
end

function libssh2_channel_eof(channel)
    ccall((:libssh2_channel_eof, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL},), channel)
end

function libssh2_channel_wait_eof(channel)
    ccall((:libssh2_channel_wait_eof, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL},), channel)
end

function libssh2_channel_close(channel)
    ccall((:libssh2_channel_close, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL},), channel)
end

function libssh2_channel_wait_closed(channel)
    ccall((:libssh2_channel_wait_closed, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL},), channel)
end

function libssh2_channel_free(channel)
    ccall((:libssh2_channel_free, libssh2), Cint, (Ptr{LIBSSH2_CHANNEL},), channel)
end

function libssh2_scp_recv(session, path, sb)
    ccall((:libssh2_scp_recv, libssh2), Ptr{LIBSSH2_CHANNEL}, (Ptr{LIBSSH2_SESSION}, Cstring, Ptr{stat}), session, path, sb)
end

function libssh2_scp_recv2(session, path, sb)
    ccall((:libssh2_scp_recv2, libssh2), Ptr{LIBSSH2_CHANNEL}, (Ptr{LIBSSH2_SESSION}, Cstring, Ptr{libssh2_struct_stat}), session, path, sb)
end

function libssh2_scp_send_ex(session, path, mode, size, mtime, atime)
    ccall((:libssh2_scp_send_ex, libssh2), Ptr{LIBSSH2_CHANNEL}, (Ptr{LIBSSH2_SESSION}, Cstring, Cint, Csize_t, Clong, Clong), session, path, mode, size, mtime, atime)
end

function libssh2_scp_send64(session, path, mode, size, mtime, atime)
    ccall((:libssh2_scp_send64, libssh2), Ptr{LIBSSH2_CHANNEL}, (Ptr{LIBSSH2_SESSION}, Cstring, Cint, libssh2_int64_t, Ctime_t, Ctime_t), session, path, mode, size, mtime, atime)
end

function libssh2_base64_decode(session, dest, dest_len, src, src_len)
    ccall((:libssh2_base64_decode, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Ptr{Cstring}, Ptr{UInt32}, Cstring, UInt32), session, dest, dest_len, src, src_len)
end

function libssh2_version(req_version_num)
    ccall((:libssh2_version, libssh2), Cstring, (Cint,), req_version_num)
end

function libssh2_knownhost_init(session)
    ccall((:libssh2_knownhost_init, libssh2), Ptr{LIBSSH2_KNOWNHOSTS}, (Ptr{LIBSSH2_SESSION},), session)
end

function libssh2_knownhost_add(hosts, host, salt, key, keylen, typemask, store)
    ccall((:libssh2_knownhost_add, libssh2), Cint, (Ptr{LIBSSH2_KNOWNHOSTS}, Cstring, Cstring, Cstring, Csize_t, Cint, Ptr{Ptr{libssh2_knownhost}}), hosts, host, salt, key, keylen, typemask, store)
end

function libssh2_knownhost_addc(hosts, host, salt, key, keylen, comment, commentlen, typemask, store)
    ccall((:libssh2_knownhost_addc, libssh2), Cint, (Ptr{LIBSSH2_KNOWNHOSTS}, Cstring, Cstring, Cstring, Csize_t, Cstring, Csize_t, Cint, Ptr{Ptr{libssh2_knownhost}}), hosts, host, salt, key, keylen, comment, commentlen, typemask, store)
end

function libssh2_knownhost_check(hosts, host, key, keylen, typemask, knownhost)
    ccall((:libssh2_knownhost_check, libssh2), Cint, (Ptr{LIBSSH2_KNOWNHOSTS}, Cstring, Cstring, Csize_t, Cint, Ptr{Ptr{libssh2_knownhost}}), hosts, host, key, keylen, typemask, knownhost)
end

function libssh2_knownhost_checkp(hosts, host, port, key, keylen, typemask, knownhost)
    ccall((:libssh2_knownhost_checkp, libssh2), Cint, (Ptr{LIBSSH2_KNOWNHOSTS}, Cstring, Cint, Cstring, Csize_t, Cint, Ptr{Ptr{libssh2_knownhost}}), hosts, host, port, key, keylen, typemask, knownhost)
end

function libssh2_knownhost_del(hosts, entry)
    ccall((:libssh2_knownhost_del, libssh2), Cint, (Ptr{LIBSSH2_KNOWNHOSTS}, Ptr{libssh2_knownhost}), hosts, entry)
end

function libssh2_knownhost_free(hosts)
    ccall((:libssh2_knownhost_free, libssh2), Cvoid, (Ptr{LIBSSH2_KNOWNHOSTS},), hosts)
end

function libssh2_knownhost_readline(hosts, line, len, type)
    ccall((:libssh2_knownhost_readline, libssh2), Cint, (Ptr{LIBSSH2_KNOWNHOSTS}, Cstring, Csize_t, Cint), hosts, line, len, type)
end

function libssh2_knownhost_readfile(hosts, filename, type)
    ccall((:libssh2_knownhost_readfile, libssh2), Cint, (Ptr{LIBSSH2_KNOWNHOSTS}, Cstring, Cint), hosts, filename, type)
end

function libssh2_knownhost_writeline(hosts, known, buffer, buflen, outlen, type)
    ccall((:libssh2_knownhost_writeline, libssh2), Cint, (Ptr{LIBSSH2_KNOWNHOSTS}, Ptr{libssh2_knownhost}, Cstring, Csize_t, Ptr{Csize_t}, Cint), hosts, known, buffer, buflen, outlen, type)
end

function libssh2_knownhost_writefile(hosts, filename, type)
    ccall((:libssh2_knownhost_writefile, libssh2), Cint, (Ptr{LIBSSH2_KNOWNHOSTS}, Cstring, Cint), hosts, filename, type)
end

function libssh2_knownhost_get(hosts, store, prev)
    ccall((:libssh2_knownhost_get, libssh2), Cint, (Ptr{LIBSSH2_KNOWNHOSTS}, Ptr{Ptr{libssh2_knownhost}}, Ptr{libssh2_knownhost}), hosts, store, prev)
end

function libssh2_agent_init(session)
    ccall((:libssh2_agent_init, libssh2), Ptr{LIBSSH2_AGENT}, (Ptr{LIBSSH2_SESSION},), session)
end

function libssh2_agent_connect(agent)
    ccall((:libssh2_agent_connect, libssh2), Cint, (Ptr{LIBSSH2_AGENT},), agent)
end

function libssh2_agent_list_identities(agent)
    ccall((:libssh2_agent_list_identities, libssh2), Cint, (Ptr{LIBSSH2_AGENT},), agent)
end

function libssh2_agent_get_identity(agent, store, prev)
    ccall((:libssh2_agent_get_identity, libssh2), Cint, (Ptr{LIBSSH2_AGENT}, Ptr{Ptr{libssh2_agent_publickey}}, Ptr{libssh2_agent_publickey}), agent, store, prev)
end

function libssh2_agent_userauth(agent, username, identity)
    ccall((:libssh2_agent_userauth, libssh2), Cint, (Ptr{LIBSSH2_AGENT}, Cstring, Ptr{libssh2_agent_publickey}), agent, username, identity)
end

function libssh2_agent_disconnect(agent)
    ccall((:libssh2_agent_disconnect, libssh2), Cint, (Ptr{LIBSSH2_AGENT},), agent)
end

function libssh2_agent_free(agent)
    ccall((:libssh2_agent_free, libssh2), Cvoid, (Ptr{LIBSSH2_AGENT},), agent)
end

function libssh2_agent_set_identity_path(agent, path)
    ccall((:libssh2_agent_set_identity_path, libssh2), Cvoid, (Ptr{LIBSSH2_AGENT}, Cstring), agent, path)
end

function libssh2_agent_get_identity_path(agent)
    ccall((:libssh2_agent_get_identity_path, libssh2), Cstring, (Ptr{LIBSSH2_AGENT},), agent)
end

function libssh2_keepalive_config(session, want_reply, interval)
    ccall((:libssh2_keepalive_config, libssh2), Cvoid, (Ptr{LIBSSH2_SESSION}, Cint, UInt32), session, want_reply, interval)
end

function libssh2_keepalive_send(session, seconds_to_next)
    ccall((:libssh2_keepalive_send, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Ptr{Cint}), session, seconds_to_next)
end

function libssh2_trace(session, bitmask)
    ccall((:libssh2_trace, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Cint), session, bitmask)
end

function libssh2_trace_sethandler(session, context, callback)
    ccall((:libssh2_trace_sethandler, libssh2), Cint, (Ptr{LIBSSH2_SESSION}, Ptr{Cvoid}, libssh2_trace_handler_func), session, context, callback)
end


text:
# generate a tridiagonal matrix
function tridiagm(a, b, c, n)
    dd, du = ones(n), ones(n - 1)
    b*diagm(0 => dd) + a*diagm(-1 => du) + c*diagm(1 => du)
end

# support function
SUITE["ρ"] = BenchmarkGroup()

# support vector
SUITE["σ"] = BenchmarkGroup()

N = Float64
for set_type in (Ball1, BallInf, Hyperrectangle)
    for n in (2, 10, 100, 1000)
        rng = MersenneTwister(n)
        d = rand(rng, n)
        X = rand(set_type, dim=n, rng=rng)
        SUITE["ρ"][string(set_type), "dense", n] = @benchmarkable ρ($d, $X)
        SUITE["σ"][string(set_type), "dense", n] = @benchmarkable σ($d, $X)
        A = tridiagm(1, -2, 0.5, n)
        Y = A * X
        SUITE["ρ"]["LinearMap of $(string(set_type))", "dense", n] = @benchmarkable ρ($d, $Y)
        SUITE["σ"]["LinearMap of $(string(set_type))", "dense", n] = @benchmarkable σ($d, $Y)
    end
end


text:
export
    apparent_sidereal_time_greenwich,
    mean_sidereal_time_greenwich

"""
    mean_sidereal_time_greenwich(jd)

Return the mean sidereal time at Greenwich.

The Julian Day number must represent Universal Time.

Parameters:
    jd : Julian Day number

Return:
    sidereal time in radians (2pi radians = 24 hrs)

```
mean_sidereal_time_greenwich(Dates.datetime2julian(DateTime(1987, 4, 10, 0, 0, 0))) |>
    radianstime_to_fday |>
    fday_to_hms

# -> (13, 10, 46.366821793410054)
```
"""
function mean_sidereal_time_greenwich(jd)
    T = jd_to_jcent(jd)
    T2 = T * T
    T3 = T2 * T
    theta0 = 280.46061837 + 360.98564736629 * (jd - 2451545.0)  + 0.000387933 * T2 - T3 / 38710000
    result = deg2rad(theta0)
    return mod2pi(result)
end

"""
    apparent_sidereal_time_greenwich(jd)

Return the apparent sidereal time at Greenwich for `jd` Julian day.

The Julian Day number must represent Universal Time.

Parameters:
    jd : Julian Day number

Return:
    sidereal time in radians (2pi radians = 24 hrs)
"""
function apparent_sidereal_time_greenwich(jd)
    # Nutation in right ascension should be computed from the DT julian date
    # we neglect the difference here
    return mod2pi(mean_sidereal_time_greenwich(jd) + nut_in_ra(jd))
end

"""
    local_hour_angle(ra, gst)

Find the local hour angle corresponding to `ra` and `gst`.

H, the local hour angle measured westwards from the South, can be calculated from
`H = θ - α`, or `H = θ0 - L - α`, where:

- θ is the local sidereal time
- θ0 the sidereal time at Greenwich
- L the observer's longitude (positive west from Greenwich)
- α is the right ascension in radians
"""
function local_hour_angle(ra, gst)
    return gst - longitude - ra
end


text:
# Example from Zhou et al., 2016, Stochastic Dual Integer Programming
#

using JuMP, Gurobi, SDDiP

# ===================================================
# P1
m1 = Model(solver=GurobiSolver(OutputFlag=0))
@variables(m1, begin
    x1, Bin
    x2, Bin
end)

@objective(m1, Min, x1 + x2)

# ===================================================
# P2
function Q(z1, z2)
    m2 = Model(solver=GurobiSolver(OutputFlag=0))
    @variables(m2, begin
        y, Int
        0 <= x1 <= 1
        0 <= x2 <= 1
    end)
    @constraints(m2, begin
        y <= 4
        y >= 2.6 - 0.5x2 - 0.25x1
    end)

    c1 = @constraint(m2, x1 == z1)
    c2 = @constraint(m2, x2 == z2)

    @objective(m2, Min, 4y)
    m2, c1, c2

end
# L(π)
function L(π1, π2, z1, z2)
    m3 = Model(solver=GurobiSolver(OutputFlag=0))
    @variables(m3, begin
        y, Int
        0 <= x1 <= 1
        0 <= x2 <= 1
    end)
    @constraints(m3, begin
        y <= 4
        y >= 2.6 - 0.5x2 - 0.25x1
    end)
    @objective(m3, Min, 4y - π1 * x1 - π2 * x2)
    m3
end

# ===================================================
# Solve the first problem
solve(m1)
z1, z2 = getvalue(x1), getvalue(x2)
# ===================================================
# Pass values to the second problem
m1, c1, c2 = Q(z1, z2)
# ===================================================
# Solve the second problem using a Lagrangian solver
method = LevelMethod(initialbound = 100.0, quadsolver=GurobiSolver(OutputFlag=0))
# method = SubgradientMethod(100.0) # gap
lp     = LinearProgramData(m1.obj, [c1; c2], [100.0; 100.0], method=method)
π = [0.0; 0.0]
val, status = lagrangiansolve!(lp, m1, π)

println("With Lagrangian solver")
println("Objective at (0,0): ", getobjectivevalue(m1))
println("Duals are:")
println(m1.linconstrDuals[c1.idx])
println(m1.linconstrDuals[c2.idx])

# Note the primal is degenerate
# ===================================================
# With Benders
println("With Benders")
m2, c1, c2 = Q(z1, z2)
solve(m2, relaxation=true)
println("Objective at (0,0): ", getobjectivevalue(m2))
println("Duals are:")
println(m2.linconstrDuals[c1.idx])
println(m2.linconstrDuals[c2.idx])

# ===================================================
# With strengthened Benders
# z1, z2 = 1, 1
π1, π2 = m2.linconstrDuals[c1.idx], m2.linconstrDuals[c2.idx]
# π1, π2 = 0, -4
m3 = L(π1, π2, z1, z2)
solve(m3)
# println(getobjectivevalue(m3) - (π1 * z1 + π2 * z2))
println("Strengthened benders")
println("Objective at (0,0): ", getobjectivevalue(m3))


text:
# Advent of Code, day 7
#
# Find target position that minimizes fuel use for a population
# of crabs that need to move there, under different assumptions
# about fuel use per each crab's distance moved (optimization
# problem, but used brute force).
#
# AK, 7 Dec 2021

# Read input data, convert to list of numbers
#raw_data = readlines(open("test_input.txt"))
raw_data = readlines(open("input.txt"))
data = map(x -> parse(Int64, x), split(raw_data[1], ","))
maxPos = maximum(data)
#println("Data:", data)
println("Max position =", maximum(data))

# Find the "horizontal" position that minimizes fuel use.  For part 1 (function
# argument true), use one unit of fuel expended per step. For part 2 (false),
# fuel use increases by 1 for each step taken.
function crabs(part1::Bool)

    bestPos = 0
    bestFuel = 0

    # Try each possible position
    for pos in 1:maxPos

        # Add up fuel use
        fuel = 0
        for crab in data
            dist = abs(crab - pos)
            if part1
                fuel += dist
            else
                fuel += sumnums(dist)
            end
        end
        #println("Position $pos: fuel = $fuel")

        # Remember if better than last (i.e., less fuel use)
        if bestPos == 0 || fuel < bestFuel
            bestPos = pos
            bestFuel = fuel
        end
    end

    # Show the best result found
    println("Best position $bestPos: fuel = $bestFuel")

end

# Sum up integers up to n
function sumnums(n)
    nn = 0
    for i in 1:n
        nn += i
    end
    return nn
end

println("--- PART 1 ---")
crabs(true)

println("--- PART 2 ---")
crabs(false)


text:
using Documenter
using ResizingTools
using BenchmarkTools

makedocs(; sitename="ResizingTools.jl", pages=["index.md", "manual.md", "references.md"])

deploydocs(; repo="github.com/wangl-cc/ResizingTools.jl.git")


text:
module VTKTests

using Test
using UnstructuredGrids.VTK
using UnstructuredGrids.Core
using UnstructuredGrids.Factories
using UnstructuredGrids.RefCellGallery: LEX_HEXAHEDRON

d = mktempdir()
f = joinpath(d,"grid")

grid = UGrid(domain=(0,1,-1,0),partition=(2,2))

writevtk(grid,f)

grid = UGrid(domain=(0,1,-1,0,2,4),partition=(2,3,4))

writevtk(grid,f)

fgrid = UGrid(grid,2)

writevtk(fgrid,f)

fgrid = UGrid(grid,1)

writevtk(fgrid,f)

writevtk(LEX_HEXAHEDRON,f)

rm(d,recursive=true)

end # module VTKTests


text:
mutable struct ExactSimulator{M} <: AbstractSimulator
  algorithm::M
  next_jump_index::Int
  next_jump_time::Float64
  t::Float64
end

##### constructors ######
function ExactSimulator(algorithm::M) where M
  return ExactSimulator{M}(algorithm, 0, 0.0, 0.0)
end

##### accessors #####

"""
```
cumulative_intensity(simulator::ExactSimulator)
```

Retrieve the cumulative intensity of the process.
"""
@inline cumulative_intensity(simulator::ExactSimulator) = first(simulator.algorithm.total_rate)

"""
```
jump_rates(simulator::ExactSimulator)
```

Retrieve the current rates for each jump.

If the cache type is `HasRates()`, then `rate[j]` is the rate for reaction channel `j`.
Otherwise if the cache type is `HasSums()`, then `rate[j]` is the cumulative `sum(rate[1:j])`.
"""
@inline jump_rates(simulator::ExactSimulator) = simulator.algorithm.rates
# TODO: Fix this for FirstReactionMethod

"""
```
next_jump_index(simulator::ExactSimulator)
```

Retrieve the index of the next jump.
"""
@inline next_jump_index(simulator::ExactSimulator) = simulator.next_jump_index

"""
```
next_jump_time(simulator::ExactSimulator)
```

Retrieve the time to the next jump.
"""
@inline next_jump_time(simulator::ExactSimulator) = simulator.next_jump_time

# start a new simulation and initialize the stepper with the next event
@inline function initialize!(simulator::ExactSimulator, state, model, tfinal)
  # reset the current system time
  simulator.t = 0.0

  # initialize the simulation algorithm
  initialize!(simulator.algorithm, state, model)

  # seed the stepper with the first event
  generate_next_jump!(simulator)

  return nothing
end

# carry out the next jump event
@inline function step!(simulator::ExactSimulator, state, model)
  # unpack information
  j = simulator.next_jump_index
  s = simulator.next_jump_time

  # update the system time
  simulator.t = get_new_time(simulator)

  # update the system state according to jump j
  @inbounds @fastmath execute_jump!(state, model, j)

  # update statistics
  # TODO

  # update the simulation algorithm given that jump j occured
  update!(simulator.algorithm, state, model, j)

  # generate the next jump
  generate_next_jump!(simulator)

  return nothing
end

# generate the next (j,s) pair by invoking the stepper's event sampler
@inline function generate_next_jump!(simulator::ExactSimulator)
  # ask the algorithm for the next random jump
  j, s = generate_jump(simulator.algorithm)

  # update the stepper's fields
  simulator.next_jump_index = j
  simulator.next_jump_time = s

  return nothing
end

@inline get_new_time(simulator::ExactSimulator) = get_new_time(simulator.algorithm, simulator.t, simulator.next_jump_time)

##### helper functions #####

function no_update_step!(simulator::ExactSimulator, state, model)
  # unpack information
  j = simulator.next_jump_index
  s = simulator.next_jump_time

  # update the system time
  simulator.t = get_new_time(simulator)

  # update the system state according to jump j
  @inbounds @fastmath execute_jump!(state, model, j)

  return nothing
end


text:
using ParserCombinator.Parsers.GML

println("celegansneural")
s = open(readall, "gml/celegansneural.gml")
parse_dict(s; debug=false)
@time x = parse_dict(s; debug=false)
@test length(x) == 2
@test x[:Creator] == "Mark Newman on Thu Aug 31 12:59:09 2006"
@test length(x[:graph][1]) == 3
@test length(x[:graph][1][:node]) == 297
@test length(x[:graph][1][:edge]) == 2359


text:
function growth_bottlenecks(filepath; debug=false, quiet=false)
    if !debug
        parameters = Dict{Symbol, Vector}(
            :num_genes => repeat([5000], inner=[3]),
            :representation => [10, 100, 1000],
            :bottleneck_representation => [10, 100, 1000],
            :seq_depth => [100, 100, 1000],
            :num_bottlenecks => collect(1:20),
            :noise => collect(linspace(0.001, 0.1, 5))
        )
        num_runs = 25
    else
        parameters = Dict{Symbol, Vector}(
            :representation => [100],
            :bottleneck_representation => [100],
            :seq_depth => [100],
            :num_bottlenecks => [5, 4],
            :noise => [0.01, 0.5]
        )
        num_runs = 1
    end

    runs = grouped_param_space(GrowthScreen(), parameters, [:num_bottlenecks, :noise], num_runs)

    methods = [venn, auprc]
    measures = [:inc, :dec, :incdec]
    genetypes = [:sigmoidal, :linear, :all]
    test_method_wrapper = (bc_counts, genes) -> test_methods(genes, methods, measures, genetypes)

    max_phenotype_dists = Dict{Symbol, Tuple{Float64, Sampleable}}(
        :inactive => (0.849, Delta(0.0)),
        :negcontrol => (0.05, Delta(0.0)),
        :increasing => (0.001, TruncatedNormal(0.1, 0.1, 0.025, 1)),
        :decreasing => (0.1, TruncatedNormal(-0.55, 0.2, -1, -0.1))
    )

    before = time()
    lib = Library(max_phenotype_dists, CRISPRi())
    results = pmap(args -> run_exp(args[1], lib, test_method_wrapper; run_idx=args[2]), runs)
    (!quiet) && println("$(time() - before) seconds")
    results = DataFrame(permutedims(hcat(results...), [2, 1]))
    results[:crisprtype] = "CRISPRi"
    before = time()
    lib = Library(max_phenotype_dists, CRISPRn())
    results2 = pmap(args -> run_exp(args[1], lib, test_method_wrapper; run_idx=args[2]), runs)
    (!quiet) && println("$(time() - before) seconds")
    results2 = DataFrame(permutedims(hcat(results2...), [2, 1]))
    results2[:crisprtype] = "CRISPRn"
    results = vcat(results, results2)

    hierarchy = vcat([hcat(item...) for item in IterTools.product(map(Symbol, methods), measures, genetypes)]...)
    new_names = [[:method, :measure, :genetype, :score]...; fieldnames(GrowthScreen)...; :run; :crisprtype]
    results = construct_hierarchical_label(hierarchy, results, new_names)
    CSV.write(filepath, results)
end


text:
# ---
# title: 1498. Number of Subsequences That Satisfy the Given Sum Condition
# id: problem1498
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Sort, Sliding Window
# link: <https://leetcode.com/problems/number-of-subsequences-that-satisfy-the-given-sum-condition/description/>
# hidden: true
# ---
#
# Given an array of integers `nums` and an integer `target`.
#
# Return the number of **non-empty** subsequences of `nums` such that the sum of
# the minimum and maximum element on it is less or equal to `target`. Since the
# answer may be too large, return it modulo `109 + 7`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [3,5,6,7], target = 9
#     Output: 4
#     Explanation: There are 4 subsequences that satisfy the condition.
#     [3] -> Min value + max value <= target (3 + 3 <= 9)
#     [3,5] -> (3 + 5 <= 9)
#     [3,5,6] -> (3 + 6 <= 9)
#     [3,6] -> (3 + 6 <= 9)
#
#
# **Example 2:**
#
#
#
#     Input: nums = [3,3,6,8], target = 10
#     Output: 6
#     Explanation: There are 6 subsequences that satisfy the condition. (nums can have repeated numbers).
#     [3] , [3] , [3,3], [3,6] , [3,6] , [3,3,6]
#
# **Example 3:**
#
#
#
#     Input: nums = [2,3,3,4,6,7], target = 12
#     Output: 61
#     Explanation: There are 63 non-empty subsequences, two of them don't satisfy the condition ([6,7], [7]).
#     Number of valid subsequences (63 - 2 = 61).
#
#
# **Example 4:**
#
#
#
#     Input: nums = [5,2,4,1,7,6,8], target = 16
#     Output: 127
#     Explanation: All non-empty subset satisfy the condition (2^7 - 1) = 127
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 105`
#   * `1 <= nums[i] <= 106`
#   * `1 <= target <= 106`
#
#
## @lc code=start
using LeetCode

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


text:
abstract type AbstractChineseTokenizer end

struct OneGramTokenizer <: AbstractChineseTokenizer
    aca::AhoCorasickAutomaton{UInt32}
    nlf::Vector{Float64}
    INF::Float64
end

function OneGramTokenizer(words::Vector{String}, freqs::Vector{Int})
    aca = AhoCorasickAutomaton{UInt32}(words)
    nlf = -log.(freqs ./ Float64(sum(freqs)))
    INF = maximum(nlf)
    return OneGramTokenizer(aca, nlf, INF)
end

function OneGramTokenizer(;filepath::String = joinpath(Pkg.dir("ChineseTokenizers"), "data", "dict"))
    words = String[]
    freqs = Int[]
    open(filepath, "r") do io
        for line in eachline(io)
            cells = split(line, " ")
            push!(words, String(cells[1]))
            push!(freqs, parse(Int, cells[2]))
        end
    end
    return OneGramTokenizer(words, freqs)
end

function split(text::AbstractString, tokenizer::OneGramTokenizer)::Vector{String}
    codes = codeunits(text)
    res = String[]
    aca = tokenizer.aca; nlf = tokenizer.nlf; INF = tokenizer.INF
    positions = reverse!(eachmatch(aca, text))
    @show sort(positions)
    if (length(positions) == 0) return [String(text)] end
    n = length(codes)
    dp = fill(INF, n + 1); dp = cumsum(dp); reverse!(dp)
    to = collect(1:n)
    for pos in positions
        if dp[pos.t + 1] + nlf[pos.i] < dp[pos.s]
            dp[pos.s] = dp[pos.t + 1] + nlf[pos.i]
            to[pos.s] = pos.t
        end
    end
    println("dp: ", dp)
    println("to: ", to)
    s = 1
    while s <= n
        t = to[s]
        if t == s
            t += 1
            while t <= n && t == to[t] t += 1 end
            t -= 1
        end
        push!(res, String(codes[s:t]))
        s = t + 1
    end
    return res
end


text:
module KalmanAD
using LinearAlgebra
using Transducers
using Transducers: R_, start, next, complete, inner, xform, wrap, unwrap, wrapping

include("utils.jl")

export LinearStateSpaceModel
export Simulator
include("statespacemodel.jl")

export LinearObservation
export observe
include("observer.jl")

export KalmanFilter
include("filter.jl")

export energy
include("energy.jl")
end # module


text:
using ImageMagick, Test

img = load(joinpath(@__DIR__,"images", "bad_exif_orientation.jpg"))
@test size(img) == (512,512)


text:
###############################################################################
#
# Make hexbin plot
#
###############################################################################
"""
    plot_hexbin!(
                ax::PyObject,
                xs::Array,
                ys::Array;
                cmap::String = "Greys",
                logbins::Bool = false,
                gridsize::Number = 25
    )
    plot_hexbin!(
                ax::PyObject,
                xs::Array,
                ys::Array,
                xlims::Array,
                ylims::Array;
                cmap::String = "Greys",
                logbins::Bool = false,
                gridsize::Number = 25
    )

Plot density plot on axis, given
- `ax` Axis to plot on
- `xs` Array of X
- `ys` Array of Y
- `cmap` Optional. Color map scheme
- `logbins` Optional. If true, use log(count) to color the bins
- `gridsize` Number of bins on both directions
- `xlim` Limits of x axis. Used to make plot region equal among subplots
- `ylim` Limits of y axis. Used to make plot region equal among subplots
"""
function plot_hexbin!(
            ax::PyObject,
            xs::Array,
            ys::Array;
            cmap::String = "Greys",
            logbins::Bool = false,
            gridsize::Number = 25
)
    if logbins
        return ax.hexbin(xs, ys; cmap=cmap, bins="log", gridsize=gridsize)
    end

    return ax.hexbin(xs, ys; cmap=cmap, gridsize=gridsize)
end

function plot_hexbin!(
            ax::PyObject,
            xs::Array,
            ys::Array,
            xlims::Array,
            ylims::Array;
            cmap::String = "Greys",
            logbins::Bool = false,
            gridsize::Number = 25
)
    newxs = [collect(xs);
             xlims[1] - (xlims[2]-xlims[1]) / 10;
             xlims[2] + (xlims[2]-xlims[1]) / 10];
    newys = [collect(ys);
             ylims[1] - (ylims[2]-ylims[1]) / 10;
             ylims[2] + (ylims[2]-ylims[1]) / 10];
    return plot_hexbin!(ax, newxs, newys; cmap=cmap, logbins=logbins,
                        gridsize=gridsize)
end


text:
# ------------------------------------------------------------------------------------------
# ## What is learning?
#
# Computers read data, as we saw in notebooks 1 and 2. We can then build functions that
# model that data to make decisions, as we saw in notebooks 3 and 5.
#
# But how do you make sure that the model actually fits the data well? In the last notebook,
# we saw that we can fiddle with the parameters of our function defining the model to reduce
# the loss function. However, we don't want to have to pick the model parameters ourselves.
# Choosing parameters ourselves works *well enough* when we have a simple model and only a
# few data points, but can quickly become extremely complex for more detailed models and
# larger data sets.
#
# Instead, we want our machine to *learn* the parameters that fit the model to our data,
# without needing us to fiddle with the parameters ourselves. In this notebook, we'll talk
# about the "learning" in machine learning.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Motivation: Fitting parameters by hand
#
# Let's go back to our example of fitting parameters from notebook 3. Recall that we looked
# at whether the amount of green in the pictures could distinguish between an apple and a
# banana, and used a sigmoid function to model our choice of "apple or banana" using the
# amount of green in an image.
# ------------------------------------------------------------------------------------------

using Plots; gr()
using Images; using Interact

σ(x,w,b) = 1 / (1 + exp(-w*x+b))

apple =  load("data/10_100.jpg")
banana = load("data/104_100.jpg")
apple_green_amount =  mean(Float64.(green.(apple)))
banana_green_amount = mean(Float64.(green.(banana)));

@manipulate for w in -10:0.01:30, b in 0:0.1:30

    plot(x->σ(x,w,b), 0, 1, label="Model", legend = :topleft, lw=3)
    scatter!([apple_green_amount],  [0.0], label="Apple")
    scatter!([banana_green_amount], [1.0], label="Banana")

end

# ------------------------------------------------------------------------------------------
# Intuitively, how did you tweak the sliders so that way the model sends apples to 0 and
# bananas to 1? Most likely, you did the following:
#
# #### Move the sliders a bit, see whether the curve moves in the right direction, and if it
# did, keep doing it.
#
# For a machine, "learning" is that same process, translated into math!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## "Learning by nudging": The process of descent
#
# Let's start to formalize this idea. In order to push the curve in the "right direction",
# we need some measurement of "how right" and "how wrong" the model is. When we translate
# the idea of a "right direction" into math, we end up with a **loss function**, `L(w, b)`,
# as we saw in notebook 5. We say that the loss function is lowest when the model `σ(x, w,
# b)` performs the best.
#
# Now we want to create a loss function that is the lowest when the apple is at `0` and the
# banana is at `1`. If the data (the amount of green) for our apple is $x_1$, then our model
# will output $σ(x_1,w, b)$ for our apple. So, we want the difference $0 - σ(x_1, w, b)$ to
# be small. Similarly, if our data for our banana (the banana's amount of green) is $x_2$,
# we want the difference $1 - σ(x_2, w, b)$ to be small.
#
# To create our loss function, let's add together the squares of the difference of the
# model's output from the desired output for the apple and the banana. We get
#
# $$ L(w,b) = (0 - σ(x_1, w, b))^2 + (1 - σ(x_2, w, b))^2. $$
#
# $L(w, b)$ is lowest when it outputs `0` for the apple and `1` for the banana, and thus the
# cost is lowest when the model "is correct".
#
# We can visualize this function by plotting it in 3D with the `surface` function.
# ------------------------------------------------------------------------------------------

# plotly()
gr()

L(w, b) = (0 - σ(apple_green_amount,w,b))^2 + (1 - σ(banana_green_amount,w,b))^2

w_range = 10:0.1:13
b_range = 0:1:20

L_values = [L(w,b) for b in b_range, w in w_range]

@manipulate for w in w_range, b in b_range
    p1 = surface(w_range, b_range, L_values, xlabel="w", ylabel="b", cam=(70,40), cbar=false, leg=false)
    scatter!(p1, [w], [b], [L(w,b)+1e-2], markersize=5, color = :blue)
    p2 = plot(x->σ(x,w,b), 0, 1, label="Model", legend = :topleft, lw=3)
    scatter!(p2, [apple_green_amount],  [0.0], label="Apple", markersize=10)
    scatter!(p2, [banana_green_amount], [1.0], label="Banana", markersize=10, xlim=(0,1), ylim=(0,1))
    plot(p1, p2, layout=(2,1))
end

# ------------------------------------------------------------------------------------------
# The blue ball on the 3D plot shows the current parameter choices, plotted as `(w,b)`.
# Shown below the 3D plot is a 2D plot of the corresponding model with those parameters.
# Notice that as the blue ball rolls down the hill, the model becomes a better fit. Our loss
# function gives us a mathematical notion of a "hill", and the process of "learning by
# nudging" is simply rolling the ball down that hill.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# To do this mathematically, we need to know which direction is "downhill". Recall from
# calculus that the derivative of `L` with respect to `b` tells you how `L` changes when `b`
# changes. Thus to roll downhill, we should go in the direction where the derivative is
# negative (the function goes down) for each parameter. This direction is the negative of
# what's called the **gradient**, $\nabla L$. This means that the "learn by nudging method"
# can be rephrased in mathematical terms as:
#
# 1. Calculate the gradient
# 2. Move a little bit in the direction of the negative gradient
# 3. Repeat
#
# This process of rolling the ball in the direction of the negative gradient is called
# **gradient descent**; written mathematically, it is
#
# $$p_{n+1} = p_n - \eta \nabla L(p_n).$$
#
# Here, $p_n$ represents the vector of current parameters $(w, b)$; $\nabla L(p_n)$ is the
# gradient of the loss function, given those parameters. We start from $p_n$ and change it
# by $\eta \nabla L(p_n)$, where $\eta$ is a small step size that determines how far we move
# the parameters in the direction of the negative gradient; notice that if you step too far,
# you'll overshoot the minimum!. The result is $p_{n+1}$, the new vector of parameters.
#
# [Picture of Gradient Descent Vectors]
#
# If we repeat this process, then we will end up at parameters where the model correctly
# labels apples as `0` and bananas as `1`. When this happens, the model has learned from the
# data and can then read pictures and tell you whether they are apples or bananas!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# Use the following terms to fill in the sentences below. Terms may be used more than once
# or not at all:
# > gradient, loss function, derivative, gradient descent, learning.
#
# * We can think of a _(A)_ as a 1D version of a _(B)_.
# * We can visualize a _(C)_ as a hill.
# * In the explanation above, rolling downhill is called _(D)_ and means traveling along the
# _(E)_.
# * To quantify the correctness of a model we use a _(F)_.
# * When our program can minimize a _(G)_ on its own, we say it is _(H)_.
#
#

#
# A)

# B)

# C)

# D)

# E)

# F)

# G)

# H)

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

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/09. ML - What is learning.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials