text:
__precompile__()
# See readme for module description.
module FilmDataStructures
export lyrDsc, dptDsc
# Data structure definition for calculating heat transfer between two layered infinite half spaces.
"""
Data structure for permittivity profile of layered structure.
Each provided permittivity in the permittivity list is presumed to be constant within an interval beginning on the previous interface depth, and terminating on the next interface depth (infinity when no other boundary exists).

# Arguments
.bdrLoc: relative positions of interfaces, units of microns.
.tmpLst: assumed (constant) temperature of each layer in Kelvin.
.rspPrf: permittivity values for each depth.
.tfrFac: layer transfer factor, imaginary part of electric susceptibility for heat transfer.
"""
struct lyrDsc

	bdrLoc::Array{Float64,1}
	tmpLst::Array{Float64,1}
	rspPrf::Any
	tfrFac::Any
end
"""
Data structure describing properties of supposed dopants.

# Arguments
.enr: offset of impurity energy from band edge in electron volts.
.ccn: impurity concentration in electron volts.
"""
# Antimonide (n-type)
# enr = 0.039
# Phosphorous, most common n-type
# enr = 0.044
# Boron, most common p-type
# enr = 0.045
struct dptDsc

	enr::Float64
	ccn::Float64
end
end


text:
module SandwichLava

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Sandwich Lava" => Ahorn.EntityPlacement(
        Maple.SandwichLava
    )
)

function Ahorn.selection(entity::Maple.SandwichLava)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

edgeColor = (246, 98, 18, 255) ./ 255
centerColor = (209, 9, 1, 102) ./ 255

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.SandwichLava, room::Maple.Room) = Ahorn.drawImage(ctx, Ahorn.Assets.lavaSanwitch, -12, -12)

end

ext: jl

max_stars_repo_name: zyowo/Ahorn

max_issues_repo_name: zyowo/Ahorn

max_forks_repo_name: zyowo/Ahorn


text:
using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic

# demonstrate the various window flags. typically you would just use the default!
control_no_titlebar = false
control_no_scrollbar = true
control_no_menu = true
control_no_move = false
control_no_resize = false
control_no_collapse = false
control_no_close = true
control_no_nav = true
control_no_background = false
control_no_bring_to_front = false

"""
    ShowControlWindow(p_open::Ref{Bool})

Show control window with camera info and controls.
"""
function ShowControlWindow(p_open::Ref{Bool})
    global cam, camImage
    global camSettings, camGPIO
    global sessionStat

    camIsRunning = isrunning(cam)

    # Control window
    # demonstrate the various window flags. typically you would just use the default!
    window_flags = CImGui.ImGuiWindowFlags(0)
    control_no_titlebar       && (window_flags |= CImGui.ImGuiWindowFlags_NoTitleBar;)
    control_no_scrollbar      && (window_flags |= CImGui.ImGuiWindowFlags_NoScrollbar;)
    !control_no_menu          && (window_flags |= CImGui.ImGuiWindowFlags_MenuBar;)
    control_no_move           && (window_flags |= CImGui.ImGuiWindowFlags_NoMove;)
    control_no_resize         && (window_flags |= CImGui.ImGuiWindowFlags_NoResize;)
    control_no_collapse       && (window_flags |= CImGui.ImGuiWindowFlags_NoCollapse;)
    control_no_nav            && (window_flags |= CImGui.ImGuiWindowFlags_NoNav;)
    control_no_background     && (window_flags |= CImGui.ImGuiWindowFlags_NoBackground;)
    control_no_bring_to_front && (window_flags |= CImGui.ImGuiWindowFlags_NoBringToFrontOnFocus;)
    control_no_close && (p_open = C_NULL;) # don't pass our bool* to Begin

    # specify a default position/size in case there's no data in the .ini file.
    # typically this isn't required! we only do it to make the Demo applications a little more welcoming.
    CImGui.SetNextWindowPos((1280-550-20, 20), CImGui.ImGuiCond_FirstUseEver)
    CImGui.SetNextWindowSize((550, 310), CImGui.ImGuiCond_FirstUseEver)

    CImGui.Begin("Control", p_open, window_flags) || (CImGui.End(); return)
    CImGui.Text("ImGui $(CImGui.IMGUI_VERSION)")

    # most "big" widgets share a common width settings by default.
    # CImGui.PushItemWidth(CImGui.GetWindowWidth() * 0.65)    # use 2/3 of the space for widgets and 1/3 for labels (default)
    CImGui.PushItemWidth(CImGui.GetFontSize() * -12)        # use fixed width for labels (by passing a negative value), the rest goes to widgets. We choose a width proportional to our font size.

    CImGui.Text(@sprintf("Display: %.1f FPS", CImGui.GetIO().Framerate))
    CImGui.Text(@sprintf("Framegrab: %.1f FPS", perfGrabFramerate))
    CImGui.Spacing()
    CImGui.Text("Settings")
    ## FRAMERATE
    @cstatic ctrl_framerate=Cfloat(12.00) begin
        ctrl_framerate = Cfloat(camSettings.acquisitionFramerate)
        CImGui.PushItemWidth(200)
        @c CImGui.SliderFloat("Framerate (fps)", &ctrl_framerate,
            camSettingsLimits.acquisitionFramerate[1],
            camSettingsLimits.acquisitionFramerate[2], "%.3f")
        camSettings.acquisitionFramerate = ctrl_framerate
    end

    ## EXPOSURE
    @cstatic ctrl_mode=0 ctrl_exposure=Cfloat(12.00) ctrl_auto=false begin
        ctrl_exposure = Cfloat(camSettings.exposureTime)
        ctrl_auto = (camSettings.exposureAuto == :continuous)
        CImGui.PushItemWidth(200)
        @c CImGui.SliderFloat("Exposure (ms)", &ctrl_exposure,
            camSettingsLimits.exposureTime[1],
            camSettingsLimits.exposureTime[2], "%.3f")
        CImGui.SameLine()
        @c CImGui.Checkbox("Auto exp.", &ctrl_auto)
        #=
        items = map(x->string(x),camSettingsLimits.exposureMode)
        CImGui.PushItemWidth(50); CImGui.SameLine()
        @cstatic item_current="timed" begin #NOT LINKED!!
            # here our selection is a single pointer stored outside the object.
            if CImGui.BeginCombo("Mode", item_current) # the second parameter is the label previewed before opening the combo.
                for n = 0:length(items)-1
                    is_selected = item_current == items[n+1]
                    CImGui.Selectable(items[n+1], is_selected) && (item_current = items[n+1];)
                    is_selected && CImGui.SetItemDefaultFocus() # set the initial focus when opening the combo (scrolling + for keyboard navigation support in the upcoming navigation branch)
                end
                CImGui.EndCombo()
            end
        end
        =#
        if ctrl_auto
            camSettings.exposureAuto = :continuous
        else
            camSettings.exposureAuto = :off
            camSettings.exposureTime = ctrl_exposure
        end
    end

    ## GAIN
    @cstatic ctrl_gain=Cfloat(0.00) ctrl_autogain=false begin
        ctrl_gain = Cfloat(camSettings.gain)
        ctrl_autogain = (camSettings.gainAuto == :continuous)
        CImGui.PushItemWidth(200)
        @c CImGui.SliderFloat("Gain (dB)", &ctrl_gain,
            camSettingsLimits.gain[1],
            camSettingsLimits.gain[2], "%.3f")
        CImGui.SameLine()
        @c CImGui.Checkbox("Auto Gain", &ctrl_autogain)

        if ctrl_autogain
            camSettings.gainAuto = :continuous
        else
            camSettings.gainAuto = :off
            camSettings.gain = ctrl_gain
        end
    end

    #==
    ## GAMMA
    @cstatic ctrl_gamma=Cfloat(12.00) begin
        ctrl_gamma = Cfloat(camSettings.gamma)
        CImGui.PushItemWidth(200)
        @c CImGui.SliderFloat("Gamma", &ctrl_gamma,
            camSettingsLimits.gamma[1],
            camSettingsLimits.gamma[2], "%.3f")
        camSettings.gamma = ctrl_gamma
    end

    ## BLACK LEVEL
    @cstatic ctrl_blacklevel=Cfloat(12.00) begin
        ctrl_blacklevel = Cfloat(camSettings.blackLevel)
        CImGui.PushItemWidth(200)
        @c CImGui.SliderFloat("Black Level", &ctrl_blacklevel,
            camSettingsLimits.blackLevel[1],
            camSettingsLimits.blackLevel[2], "%.3f")
        camSettings.blackLevel = ctrl_blacklevel
    end
    =#

    CImGui.Spacing()
    CImGui.Text("Image Format")

    ## IMAGE OFFSET X
    @cstatic ctrl_val=Cint(0) begin
        ctrl_val = Cint(camSettings.offsetX)
        CImGui.PushItemWidth(200)
        @c CImGui.SliderInt("Offset X", &ctrl_val,
            camSettingsLimits.offsetX[1],
            camSettingsLimits.offsetX[2], "%i")
        camSettings.offsetX = 2*(round(ctrl_val/2,RoundToZero))
    end
    if !camIsRunning
        CImGui.SameLine()
        ## IMAGE WIDTH
        @cstatic ctrl_val=Cint(100) begin
            ctrl_val = Cint(camSettings.width)
            CImGui.PushItemWidth(200)
            @c CImGui.SliderInt("Width", &ctrl_val,
                camSettingsLimits.width[1],
                camSettingsLimits.width[2], "%i")
            camSettings.width = 2*(round(ctrl_val/2,RoundToZero))
        end
    end
    ## IMAGE OFFSET Y
    @cstatic ctrl_val=Cint(0) begin
        ctrl_val = Cint(camSettings.offsetY)
        CImGui.PushItemWidth(200)
        @c CImGui.SliderInt("Offset Y", &ctrl_val,
            camSettingsLimits.offsetY[1],
            camSettingsLimits.offsetY[2], "%i")
        camSettings.offsetY = 2*(round(ctrl_val/2,RoundToZero))
    end
    if !camIsRunning
        CImGui.SameLine()
        ## IMAGE HEIGHT
        @cstatic ctrl_val=Cint(100) begin
            ctrl_val = Cint(camSettings.height)
            CImGui.PushItemWidth(200)
            @c CImGui.SliderInt("Height", &ctrl_val,
                camSettingsLimits.height[1],
                camSettingsLimits.height[2], "%i")
            camSettings.height = 2*(round(ctrl_val/2,RoundToZero))
        end
    end

    ## PLAY/PAUSE
    @cstatic begin
        if camIsRunning
            if CImGui.Button("Pause",(200,25))
                stop!(cam)
                @info "Camera stopped"
            end
        else
            if CImGui.Button("Run",(200,25))
                startcheckrunningfix!(cam)
            end
        end
    end

    ## RECORDING
    @cstatic begin
        if sessionStat.recording
            if CImGui.Button("Stop Recording",(200,25))
                sessionStat.recording = false
                buffermode!(cam,"NewestOnly")
            end
        else
            if CImGui.Button("Record",(200,25))
                sessionStat.recording = true
                buffermode!(cam,"OldestFirst")
            end
        end
        CImGui.SameLine()
        CImGui.Text(@sprintf("Frames: [%i buffered. %i dropped. %i saved]", sessionStat.bufferedframes,sessionStat.droppedframes,sessionStat.savedframes))
    end

    CImGui.End() # demo
end


text:
using DifferentialEquations, LinearAlgebra
using Plots
plotly()

#=
m*x`` + b*x` + k*x + a*x^3 = -mg
y = x`
=#

function fun(du,u,p,t)
    m,k,g,a,b = p
    du[1] = u[2]
    du[2] = -1/m*(b*u[2]+ k*u[1]+ a*u[1]^3+ m*g)
end

m, k, g, a = 1220.0, 35600.0, 17.5, 450000.0
b= 1000.0

tspan = (0.0, 10.0)
u0 = [0.0, 5.0]

prob1 = ODEProblem(fun, u0, tspan, (m, k, g, a, b))
sol = solve(prob1, Tsit5(),reltol = 1e-4, abstol=1e-10)

plot(sol)

max_stars_repo_path: src/solvind second order diffEq.jl

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_path: src/solvind second order diffEq.jl

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_path: src/solvind second order diffEq.jl

max_forks_repo_name: SabbahMohammed/Dev


text:
import ACEatoms
import ACEatoms: AtomState, ACESitePotential

const AtomStateU1{T} = ACE.State{(:mu, :rr, :u),
                                 Tuple{AtomicNumber, SVector{3, T}, T}}

const AtEnvDef = AtomicEnvironment{AtomState{Float64}}
const AtEnvU1 = AtomicEnvironment{AtomStateU1{Float64}}

struct I2LModel{TL1, TL2} <: AbstractACEModel
   L1::TL1
   L2::TL2
end

# ------------ Model construction codes

I2LModel(L1::Dict, L2::Dict) =
      I2LModel( ACESitePotential(L1), ACESitePotential(L2) )

# ------------ Parameter wrangling

nparams(m::I2LModel) = nparams(m.L1) + nparams(m.L2)

params(m::I2LModel) = [ params(m.L1); params(m.L2) ]

function set_params!(m::I2LModel, c)
   n1 = nparams(m.L1)
   set_params!(m.L1, c[1:n1])
   set_params!(m.L2, c[(n1+1):end])
end

# ----------- evaluation codes

function AtEnvU1(at, nlist, i, U)
   Js, Rs, Zs = JuLIP.Potentials.neigsz(nlist, at, i)
   z0 = at.Z[i]
   X0 = AtomStateU1{Float64}( (mu = z0, u = U[i]) )
   Xs = [ AtomStateU1{Float64}( (mu = Zs[t], rr = Rs[t], u = U[Js[t]]) )
          for t = 1:length(Js) ]
   return AtomicEnvironment(X0, Xs)
end

alloc_temp(m::I2LModel, args...) =
      ( tmpL1 = alloc_temp(m.L1, args...),
        tmpL2 = alloc_temp(m.L2, args...)
      )

function JuLIP.energy(m::I2LModel, at::Atoms)
   # evaluate the first layer
   nlist = JuLIP.neighbourlist(at, cutoff(m.L1))
   U = zeros(Float64, length(at))
   for i = 1:length(at)
      env, Js = environment(m.L1, at, nlist, i)
      U[i] = evaluate(m.L1.models[at.Z[i]], env)
   end

   # filter the U values to go between 0, 1
   map!(u -> 1 / (1 + exp(u)), U, U)

   # evaluate the second layer
   E = zero(fltype(m.L2))
   nlist = JuLIP.neighbourlist(at, cutoff(m.L2 ))  # could technically be different rcut
   for i = 1:length(at)
      envU = AtEnvU1(at, nlist, i, U)
      E += evaluate(m.L2.models[at.Z[i]], envU)
   end

   return E.val
end


text:
type Variable{T,N} <: Node{T}
    name::Symbol
    tensor::Array{T, N}
end
Variable{T}(name::Symbol, ::Type{T}, shape::Int...) = Variable(name, Array{T}(shape...))

name(a::Variable) = a.name
output(a::Variable) = a.tensor

function Base.push!{T}(net::ForwardNet{T}, ::Type{Variable},
    name::Symbol,
    tensor::Array{T},
    )

    push!(net, Variable(name, tensor))
end
function Base.push!{T}(net::ForwardNet{T}, ::Type{Variable},
    name::Symbol,
    shape::Int...
    )

    push!(net, Variable(name, Array{T}(shape...)))
end
function Base.push!{T}(net::ForwardNet{T}, ::Type{Variable},
    name::Symbol,
    parent::NameOrIndex,
    ::Symbol
    )#

    node = Variable(name, output(net[parent]))
    push!(net, node, parent)
end


text:
export ModelKitHomotopy

"""
    ModelKitHomotopy(system::ModelKit.System) <: AbstractSystem
"""
struct ModelKitHomotopy{T,P<:Union{Nothing,AbstractVector}} <: AbstractHomotopy
    homotopy::ModelKit.CompiledHomotopy{T}
    parameters::P
end

function ModelKitHomotopy(H::ModelKit.Homotopy; parameters = nothing)
    ModelKitHomotopy(ModelKit.compile(H), parameters)
end

Base.size(H::ModelKitHomotopy) = size(H.homotopy)

function set_parameters!(H::ModelKitHomotopy, p::Tuple, γ::Union{Nothing,Tuple} = nothing)
    p1, p2 = p
    n = length(p1)
    for i = 1:n
        H.parameters[i] = p1[i]
    end
    for i = 1:n
        H.parameters[i+n] = p2[i]
    end

    if γ !== nothing
        γ₁, γ₂ = γ
        H.parameters[2n+1] = γ₁
        H.parameters[2n+2] = γ₂
    end

    H
end

"""
    set_start_parameters!(H::ModelKitHomotopy, p)

Update the start parameters of `H`. This assumes that the first half of the parameters
are the start parameters.
"""
function set_start_parameters!(H::ModelKitHomotopy, p)
    for i = 1:div(length(p), 2)
        H.parameters[i] = p[i]
    end
    H
end

"""
    set_target_parameters!(H::ModelKitHomotopy, p)

Update the target parameters of `H`. This assumes that the second half of the parameters
are the start parameters.
"""
function set_target_parameters!(H::ModelKitHomotopy, p)
    n = div(length(p), 2)
    for i = 1:n
        H.parameters[i+n] = p[i]
    end
    H
end

cache(H::ModelKitHomotopy, x, t, p = nothing) = HomotopyNullCache()

function evaluate!(u, H::ModelKitHomotopy, x, t, ::HomotopyNullCache)
    ModelKit.evaluate!(u, H.homotopy, x, t, H.parameters)
end
function evaluate(H::ModelKitHomotopy, x, t, ::HomotopyNullCache)
    convert(Vector, ModelKit.evaluate(H.homotopy, x, t, H.parameters))
end
(H::ModelKitHomotopy)(x, t) = evaluate(F, x, t, H.parameters)

function jacobian!(U, H::ModelKitHomotopy, x, t, ::HomotopyNullCache)
    ModelKit.jacobian!(U, H.homotopy, x, t, H.parameters)
end
function jacobian(H::ModelKitHomotopy, x, t, ::HomotopyNullCache)
    convert(Matrix, ModelKit.jacobian(H.homotopy, x, t, H.parameters))
end

function evaluate_and_jacobian!(u, U, H::ModelKitHomotopy, x, t, ::HomotopyNullCache)
    ModelKit.evaluate_and_jacobian!(u, U, H.homotopy, x, t, H.parameters)
end

function dt!(u, H::ModelKitHomotopy, x, t, ::HomotopyNullCache)
    ModelKit.dt!(u, H.homotopy, x, t, H.parameters)
end
function dt(H::ModelKitHomotopy, x, t, ::HomotopyNullCache)
    convert(Vector, ModelKit.dt(H.homotopy, x, t, H.parameters))
end

function jacobian_and_dt!(U, u, H::ModelKitHomotopy, x, t, ::HomotopyNullCache)
    ModelKit.jacobian_and_dt!(U, u, H.homotopy, x, t, H.parameters)
end
function jacobian_and_dt(H::ModelKitHomotopy, x, t, ::HomotopyNullCache)
    ModelKit.jacobian_and_dt(H.homotopy, x, t, H.parameters)
end


text:
# module model_fit

using DifferentialEquations, Flux, Optim, DiffEqFlux, DiffEqSensitivity, Plots
using Zygote, ForwardDiff, ReverseDiff
using BlackBoxOptim, GalacticOptim
using BenchmarkTools
using RecursiveArrayTools
# -------------------------------------------------------------------------

mutable struct M
    model   # model function
    u₀      # initial conditions
    tspan   # time span
    dens    # site densities
    cons    # constant parameters
end

mutable struct DissData
    y
    t
    f
end

model = bimolecular_diss!
u₀ = [1.0]
tspan =  (0.0,10.0)
dens = []
cons = [1]
γδM = M(model,u₀,tspan,dens,cons)
# γδM.model = bimolecular!

function optimization(M,data,options)
# fit function - given model, data, densities, options get best parameters
    unique_t = [unique(times) for times in data.t]
    loss_in = (p) -> loss(p,data,M,unique_t)
    callback_in = (p,l) -> callback(p,l,data,M,unique_t)

    lb,ub = options

    @show loss_in(pinit)

    sr = [(0.0,1.0),(0.0,100.0)]
    res_bbo = bboptimize(loss_in; SearchRange = sr,
    # res_bbo = bboptimize(loss_in; SearchRange = options,
                                NumDimensions = length(pinit),
                                Method = :adaptive_de_rand_1_bin_radiuslimited,
                                # Method = :separable_nes,
                                NThreads=Threads.nthreads()-1,
                                MaxSteps = 20000)#,

    p₀_bbo = best_candidate(res_bbo)
    @show p₀_bbo

    res = optimize(loss_in,lb,ub,p₀_bbo, Fminbox(BFGS()), autodiff = :forward,
                                Optim.Options(show_trace=true,
                                              f_tol = 1e-5,
                                              outer_iterations = 10))

    rates = res.minimizer
    l = res.minimum

    # bootstrap_rates = bootstrapping(rates,data,model,dens,cons,t)

    @show M.model, l, rates

    hes_zyg = []
    bootstrap_rates = []

    return rates, l
end

function run_model(M,p,dens,cons,t) # version of run function with multiple models
    # f = (du,u,p,t) -> model(du,u,10 .^ p,t,densities,constants)
    f = (du,u,p,t) -> M.model(du,u,p,t,dens,cons)
    tspan = (0.0,maximum(t))
    # # tmp_prob = ODEProblem(f,M.u₀,M.tspan,p) #allows to specify tspan
    # tmp_prob = ODEProblem(f,M.u₀,tspan,p)  #automatic tspan

    # # uncomment to fit the initial conditions as well...
    # # #u₀ = [(1-p[1])/2, (1-p[1])/2, p[1]]
    u₀ = [(1-p[1])*0.610, (1-p[1])*0.390, p[1]] #ratio depends on γ and δ affinities
    p = p[2:end]
    tmp_prob = ODEProblem(f,u₀,tspan,p)
    # # end fit initial conditions
    tmp_sol = solve(tmp_prob,AutoVern7(KenCarp4()),saveat=t, abstol=1e-8,reltol=1e-8,verbose=false)
end

function loss(p,data,M,unique_t)
    Δy_all = [] # keep track of errors for all curves (change this to array of arrays!)
    # Δy_all = Zygote.Buffer([[0.,0.],[0.,0.]],length(data))
    for i in 1:length(data.y) # for each bin
      # run model with densities and parameter/dens/cons combinations
      tmp_sol = run_model(M, p, M.dens[i], M.cons[i], unique_t[i])

      if tmp_sol.retcode != :Success #if there is an error return infitite loss
        Δy_sol = fill(Inf, length(data.t[i]))
      else
        Δy_sol = find_diffs(tmp_sol,data.t[i],unique_t[i],data.y[i])
      end
      Δy_all = vcat(Δy_all,Δy_sol)
    end

    t_all = reduce(vcat,data.t)

    Losses = Δy_all .^2 .* t_all .^2
    L = sum(Losses) # calculate error across all differences from all curves

    return L
end

function find_diffs(tmp_sol,t,unique_t,y) # solutions use unique time points,
# need to evaluate the solution at each experimental time point as well as
# calculating the difference between fit and experiment

    Σ_sol = sum(Array(tmp_sol),dims=1)
    Δy = []
    for (i, tᵢ) in enumerate(t)
        ind = isequal.(unique_t,tᵢ)
        yᵢ_exp = Σ_sol[ind][1]
        yᵢ_obs = y[i]
        Δyᵢ = yᵢ_obs .- yᵢ_exp
        Δy=vcat(Δy,Δyᵢ)
    end
    return copy(Δy) #Δy[1:end] # copy(Δy)
end

# function callback(p,l,tmp_sol,data,model,dens,cons,t)
function callback(p,l,data,M,unique_t)
  # @show l, p
  @show l

  false
end


text:
"""
    Get the rhs for a strictly less than set such that it cna be used as a <= constraint.

"""
function get_rhs_from_strictly(com::CS.CoM, constraint::LinearConstraint,
        fct, set::Strictly{T, MOI.LessThan{T}}) where {T<:Real}

    constraint.is_rhs_strong && return constraint.rhs

    constraint.is_rhs_strong = true
    constraint.rhs = set.set.upper - fct.constant

    # change the rhs in such a way that <= can be used
    # if all coefficients are 1 or -1
    if all(abs(term.coefficient) == 1 for term in fct.terms)
        # => just subtract one from the rhs if rhs is discrete
        if isapprox_discrete(com, constraint.rhs)
            return constraint.rhs - 1
        else # otherwise round down so 9.5 => 9
            return floor(constraint.rhs)
        end
    else
        # use lp solver if it exists and if it supports MIP
        fallback_rhs = constraint.rhs - com.options.atol
        com.options.lp_optimizer === nothing && return fallback_rhs

        if !MOI.supports_constraint(
            com.options.lp_optimizer.optimizer_constructor(),
            SVF, MOI.Integer
        ) || !MOI.supports_constraint(
            com.options.lp_optimizer.optimizer_constructor(),
            typeof(constraint.fct),
            typeof(constraint.set.set),
        ) || !MOI.supports(com.options.lp_optimizer.optimizer_constructor(),
              MOI.ObjectiveFunction{typeof(constraint.fct)}()
        ) || !MOI.supports(com.options.lp_optimizer.optimizer_constructor(), MOI.ObjectiveSense())
            return fallback_rhs
        end

        # supports MIP and <=
        mip_model = Model()
        mip_backend = backend(mip_model)

        set_optimizer(mip_model, com.options.lp_optimizer)
        lp_x = Vector{VariableRef}(undef, length(com.search_space))
        for variable in com.search_space
            lp_x[variable.idx] = @variable(
                mip_model,
                lower_bound = variable.lower_bound,
                upper_bound = variable.upper_bound,
                integer = true
            )
        end
        MOI.add_constraint(mip_backend, constraint.fct, typeof(constraint.set.set)(set.set.upper - com.options.atol))
        MOI.set(mip_backend, MOI.ObjectiveFunction{typeof(constraint.fct)}(), constraint.fct)
        MOI.set(mip_backend, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        optimize!(mip_model)

        if termination_status(mip_model) == MOI.OPTIMAL
            return objective_value(mip_model)
        else
            # TODO: return that it's not feasible
            return fallback_rhs
        end
    end
end

"""
    init_constraint!(com::CS.CoM, constraint::CS.LinearConstraint,fct::SAF{T}, set::Strictly{MOI.LessThan{T}})

Initialize the LinearConstraint by checking whether it might be an unfillable constraint
without variable i.e x == x-1 => x -x == -1 => 0 == -1 => return false
"""
function init_constraint!(
    com::CS.CoM,
    constraint::CS.LinearConstraint,
    fct::SAF{T},
    set::Strictly{T, MOI.LessThan{T}};
) where {T<:Real}
    # rhs will be changed to use as <=
    constraint.rhs = get_rhs_from_strictly(com, constraint, fct, set)
    constraint.strict_rhs = set.set.upper - fct.constant
    constraint.is_strict = true
    is_no_variable_constraint(constraint) && return fct.constant < set.set.upper

    return true
end

"""
    init_constraint!(com::CS.CoM, constraint::CS.LinearConstraint,fct::SAF{T}, set::MOI.LessThan{T})

Initialize the LinearConstraint by checking whether it might be an unfillable constraint
without variable i.e x == x-1 => x -x == -1 => 0 == -1 => return false
"""
function init_constraint!(
    com::CS.CoM,
    constraint::CS.LinearConstraint,
    fct::SAF{T},
    set::MOI.LessThan{T};
) where {T<:Real}
    constraint.rhs = set.upper - fct.constant
    is_no_variable_constraint(constraint) && return fct.constant <= set.upper
    return true
end

function init_constraint!(
    com::CS.CoM,
    constraint::CS.LinearConstraint,
    fct::SAF{T},
    set::MOI.EqualTo{T};
) where {T<:Real}
    constraint.rhs = set.value - fct.constant
    is_no_variable_constraint(constraint) && return fct.constant == set.value
    return true
end

"""
    get_new_extrema_and_sum(search_space, vidx, i, terms, full_min, full_max, pre_mins, pre_maxs)

Get the updated full_min, full_max as well as updated pre_mins[i] and pre_maxs[i]
after values got removed from search_space[vidx]
Return full_min, full_max, pre_mins[i], pre_maxs[i]
"""
function get_new_extrema_and_sum(
    search_space,
    vidx,
    i,
    terms,
    full_min,
    full_max,
    pre_mins,
    pre_maxs,
)
    new_min = pre_mins[i]
    new_max = pre_maxs[i]
    if terms[i].coefficient > 0
        coeff_min = search_space[vidx].min * terms[i].coefficient
        coeff_max = search_space[vidx].max * terms[i].coefficient
        full_max -= (coeff_max - pre_maxs[i])
        full_min += (coeff_min - pre_mins[i])
        new_min = coeff_min
        new_max = coeff_max
    else
        coeff_min = search_space[vidx].max * terms[i].coefficient
        coeff_max = search_space[vidx].min * terms[i].coefficient
        full_max -= (coeff_max - pre_maxs[i])
        full_min += (coeff_min - pre_mins[i])
        new_min = coeff_min
        new_max = coeff_max
    end
    return full_min, full_max, new_min, new_max
end

function recompute_lc_extrema!(
    com::CS.CoM,
    constraint::LinearConstraint,
    fct::SAF{T},
) where {T<:Real}
    indices = constraint.indices
    search_space = com.search_space
    maxs = constraint.maxs
    mins = constraint.mins
    pre_maxs = constraint.pre_maxs
    pre_mins = constraint.pre_mins

    for (i, vidx) in enumerate(indices)
        if fct.terms[i].coefficient >= 0
            max_val = search_space[vidx].max * fct.terms[i].coefficient
            min_val = search_space[vidx].min * fct.terms[i].coefficient
        else
            min_val = search_space[vidx].max * fct.terms[i].coefficient
            max_val = search_space[vidx].min * fct.terms[i].coefficient
        end
        maxs[i] = max_val
        mins[i] = min_val
        pre_maxs[i] = max_val
        pre_mins[i] = min_val
    end
end

"""
    get_fixed_rhs(com::CS.CoM, constraint::Constraint)

Compute the fixed rhs based on all already fixed variables
"""
function get_fixed_rhs(com::CS.CoM, constraint::Constraint)
    rhs = constraint.rhs
    fct = constraint.fct
    search_space = com.search_space
    @inbounds for li in 1:length(constraint.indices)
        vidx = constraint.indices[li]
        var = search_space[vidx]
        !isfixed(var) && continue
        rhs -= CS.value(var) * fct.terms[li].coefficient
    end
    return rhs
end

"""
    set_new_extrema(i, pre_mins, pre_maxs, new_min, new_max)

Update pre_mins and pre_maxs if new_min or new_max is smaller or bigger.
Return wheter a value was updated
"""
function set_new_extrema(i, pre_mins, pre_maxs, new_min, new_max)
    changed = false
    if new_min != pre_mins[i]
        changed = true
        pre_mins[i] = new_min
    end
    if new_max != pre_maxs[i]
        changed = true
        pre_maxs[i] = new_max
    end
    return changed
end

"""
    prune_constraint!(com::CS.CoM, constraint::LinearConstraint, fct::SAF{T}, set::MOI.EqualTo{T}; logs = true) where T <: Real

Reduce the number of possibilities given the equality `LinearConstraint` .
Return if still feasible and throw a warning if infeasible and `logs` is set to `true`
"""
function prune_constraint!(
    com::CS.CoM,
    constraint::LinearConstraint,
    fct::SAF{T},
    set::Union{MOI.LessThan, MOI.EqualTo, Strictly{T, MOI.LessThan{T}}};
    logs = true,
) where {T<:Real}
    indices = constraint.indices
    search_space = com.search_space
    rhs = constraint.rhs

    # compute max and min values for each index
    recompute_lc_extrema!(com, constraint, fct)
    maxs = constraint.maxs
    mins = constraint.mins
    pre_maxs = constraint.pre_maxs
    pre_mins = constraint.pre_mins

    # for each index compute the maximum and minimum value possible
    # to fulfill the constraint
    full_max = sum(maxs) - rhs
    full_min = sum(mins) - rhs

    # if the maximum is smaller than 0 (and not even near zero)
    # or if the minimum is bigger than 0 (and not even near zero)
    # the equation can't sum to 0 => infeasible
    if full_min > com.options.atol
        com.bt_infeasible[indices] .+= 1
        return false
    end
    if full_max < -com.options.atol && constraint.is_equal
        com.bt_infeasible[indices] .+= 1
        return false
    end

    changed = true
    while changed
        changed = false
        for (i, vidx) in enumerate(indices)
            if isfixed(search_space[vidx])
                continue
            end
            # minimum without current index
            c_min = full_min - mins[i]

            # maximum without current index
            c_max = full_max - maxs[i]

            p_max = -c_min
            if p_max < maxs[i]
                maxs[i] = p_max
            end

            if constraint.is_equal
                p_min = -c_max
                if p_min > mins[i]
                    mins[i] = p_min
                end
            end
        end

        # update all
        for (i, vidx) in enumerate(indices)
            # if the maximum of coefficient * variable got reduced
            # get a safe threshold because of floating point errors
            if maxs[i] < pre_maxs[i]
                if fct.terms[i].coefficient > 0
                    threshold =
                        get_safe_upper_threshold(com, maxs[i], fct.terms[i].coefficient)
                    still_feasible = remove_above!(com, search_space[vidx], threshold)
                else
                    threshold =
                        get_safe_lower_threshold(com, maxs[i], fct.terms[i].coefficient)
                    still_feasible = remove_below!(com, search_space[vidx], threshold)
                end
                full_min, full_max, new_min, new_max = get_new_extrema_and_sum(
                    search_space,
                    vidx,
                    i,
                    fct.terms,
                    full_min,
                    full_max,
                    pre_mins,
                    pre_maxs,
                )
                changed = set_new_extrema(i, pre_mins, pre_maxs, new_min, new_max)
                mins[i] = pre_mins[i]
                maxs[i] = pre_maxs[i]
                !still_feasible && return false
            end
            # same if a better minimum value could be achieved
            if mins[i] > pre_mins[i]
                new_min = pre_mins[i]
                new_max = pre_maxs[i]
                if fct.terms[i].coefficient > 0
                    threshold =
                        get_safe_lower_threshold(com, mins[i], fct.terms[i].coefficient)
                    still_feasible = remove_below!(com, search_space[vidx], threshold)
                else
                    threshold =
                        get_safe_upper_threshold(com, mins[i], fct.terms[i].coefficient)
                    still_feasible = remove_above!(com, search_space[vidx], threshold)
                end
                full_min, full_max, new_min, new_max = get_new_extrema_and_sum(
                    search_space,
                    vidx,
                    i,
                    fct.terms,
                    full_min,
                    full_max,
                    pre_mins,
                    pre_maxs,
                )
                changed = set_new_extrema(i, pre_mins, pre_maxs, new_min, new_max)
                mins[i] = pre_mins[i]
                maxs[i] = pre_maxs[i]

                !still_feasible && return false
            end
        end
    end

    # the following is only for an equal constraint
    !constraint.is_equal && return true
    n_unfixed = count_unfixed(com, constraint)
    n_unfixed != 2 && return true

    return prune_is_equal_two_var!(com, constraint, fct)
end

"""
    prune_is_equal_two_var!(com::CS.CoM,
        constraint::CS.LinearConstraint,
        fct::SAF{T}) where T

Prune something like ax+by+cz == d with a,b,c,d constants and x,y,z variables
where only two variables aren't fixed yet. One should make sure with [`count_unfixed`](@ref)
that this is the case.
Return whether feasible
"""
function prune_is_equal_two_var!(com::CS.CoM,
    constraint::CS.LinearConstraint,
    fct::SAF{T}) where T

    search_space = com.search_space

    fixed_rhs = get_fixed_rhs(com, constraint)

    local_vidx_1, vidx_1, local_vidx_2, vidx_2 = get_two_unfixed(com, constraint)
    var1 = com.search_space[vidx_1]
    var2 = com.search_space[vidx_2]
    both_in_same_all_different = any(var1.in_all_different .& var2.in_all_different)
    for ((this_local_vidx, this_vidx), (other_local_vidx, other_vidx)) in zip(
        ((local_vidx_1, vidx_1), (local_vidx_2, vidx_2)),
        ((local_vidx_2, vidx_2), (local_vidx_1, vidx_1))
    )
        this_var = search_space[this_vidx]
        other_var = search_space[other_vidx]
        for val in values(this_var)
            # if we choose this value but the other wouldn't be an integer
            # => remove this value
            if !isapprox_divisible(
                com,
                (fixed_rhs - val * fct.terms[this_local_vidx].coefficient),
                fct.terms[other_local_vidx].coefficient,
            )
                if !rm!(com, this_var, val)
                    return false
                end
                continue
            end

            remainder = fixed_rhs - val * fct.terms[this_local_vidx].coefficient
            remainder /= fct.terms[other_local_vidx].coefficient
            remainder_int = get_approx_discrete(remainder)
            if !has(other_var, remainder_int)
                !rm!(com, this_var, val) && return false
                continue
            end
            # if in same all different and the remainder equals the value
            # like x+y = 2x
            if both_in_same_all_different && remainder_int == val
                !rm!(com, this_var, val) && return false
                !rm!(com, other_var, val) && return false
            end
        end
    end

    # if only one is unfixed => fix it
    var_1 = search_space[vidx_1]
    var_2 = search_space[vidx_2]
    if isfixed(var_1) || isfixed(var_2)
        if isfixed(var_1)
            unfixed_var = var_2
            unfixed_local_idx = local_vidx_2
            fixed_local_idx = local_vidx_1
        else
            unfixed_var = var_1
            unfixed_local_idx = local_vidx_1
            fixed_local_idx = local_vidx_2
        end
        fixed_coeff = fct.terms[fixed_local_idx].coefficient
        unfixed_coeff = fct.terms[unfixed_local_idx].coefficient
        if isfixed(var_1)
            fixed_rhs -= fixed_coeff*value(var_1)
        else
            fixed_rhs -= fixed_coeff*value(var_2)
        end
        !isapprox_divisible(com, fixed_rhs, unfixed_coeff) && return false
        remainder = fixed_rhs
        remainder /= unfixed_coeff
        remainder_int = get_approx_discrete(remainder)
        !has(unfixed_var, remainder_int) && return false
        !fix!(com, unfixed_var, remainder_int) && return false
    end
    return true
end

"""
    still_feasible(com::CoM, constraint::LinearConstraint, fct::SAF{T},
                   set::MOI.EqualTo{T}, vidx::Int, val::Int) where T <: Real

Return whether setting `search_space[vidx]` to `val` is still feasible given `constraint`.
"""
function still_feasible(
    com::CoM,
    constraint::LinearConstraint,
    fct::SAF{T},
    set::Union{MOI.LessThan, MOI.EqualTo, Strictly{T, MOI.LessThan{T}}},
    vidx::Int,
    val::Int,
) where {T<:Real}
    search_space = com.search_space
    rhs = constraint.rhs
    csum = 0
    num_not_fixed = 0
    not_fixed_idx = 0
    not_fixed_i = 0
    max_extra = 0
    min_extra = 0
    for (i, cvidx) in enumerate(constraint.indices)
        if cvidx == vidx
            csum += val * fct.terms[i].coefficient
            continue
        end
        if isfixed(search_space[cvidx])
            csum += CS.value(search_space[cvidx]) * fct.terms[i].coefficient
        else
            num_not_fixed += 1
            not_fixed_idx = cvidx
            not_fixed_i = i
            if fct.terms[i].coefficient >= 0
                max_extra += search_space[cvidx].max * fct.terms[i].coefficient
                min_extra += search_space[cvidx].min * fct.terms[i].coefficient
            else
                min_extra += search_space[cvidx].max * fct.terms[i].coefficient
                max_extra += search_space[cvidx].min * fct.terms[i].coefficient
            end
        end
    end
    if num_not_fixed == 0
        if constraint.is_equal
            return isapprox(csum, rhs; atol = com.options.atol, rtol = com.options.rtol)
        else
            if constraint.is_strict
                return csum < constraint.strict_rhs
            else
                return csum <= rhs
            end
        end
    end
    if num_not_fixed == 1 && constraint.is_equal
        if isapprox_divisible(com, rhs - csum, fct.terms[not_fixed_i].coefficient)
            return has(
                search_space[not_fixed_idx],
                get_approx_discrete((rhs - csum) / fct.terms[not_fixed_i].coefficient),
            )
        else
            return false
        end
    end

    if csum + min_extra > rhs + com.options.atol
        return false
    end

    if constraint.is_equal && csum + max_extra < rhs - com.options.atol
        return false
    end

    return true
end

function is_constraint_solved(
    constraint::LinearConstraint,
    fct::SAF{T},
    set::MOI.EqualTo{T},
    values::Vector{Int}
) where {T<:Real}
    return sum(values .* constraint.coeffs) + fct.constant ≈ set.value
end

function is_constraint_solved(
    constraint::LinearConstraint,
    fct::SAF{T},
    set::MOI.LessThan{T},
    values::Vector{Int}
) where {T<:Real}
    return sum(values .* constraint.coeffs) + fct.constant <= set.upper
end

function is_constraint_solved(
    constraint::LinearConstraint,
    fct::SAF{T},
    set::Strictly{T, MOI.LessThan{T}},
    values::Vector{Int}
) where {T<:Real}
    return sum(values .* constraint.coeffs) + fct.constant < set.set.upper
end

"""
    is_constraint_violated(
        com::CoM,
        constraint::LinearConstraint,
        fct::SAF{T},
        set::MOI.EqualTo{T}
    ) where {T<:Real}

Checks if the constraint is violated as it is currently set. This can happen inside an
inactive reified or indicator constraint.
"""
function is_constraint_violated(
    com::CoM,
    constraint::LinearConstraint,
    fct::SAF{T},
    set::Union{MOI.LessThan, MOI.EqualTo, Strictly{T, MOI.LessThan{T}}},
) where {T<:Real}
    if all(isfixed(var) for var in com.search_space[constraint.indices])
        return !is_constraint_solved(
            constraint,
            fct,
            set,
            CS.value.(com.search_space[constraint.indices]),
        )
    end
    return false
end


text:
@testset "Promotions" begin
    @test 1v1 + 1.0v1 == 2.0v1
    @test 1v1 + 1.0v2 == 1.0v1 + 1.0v2
    @test 1v1 - 1.0v2 == 1.0v1 - 1.0v2
    @test 1v1 - 1.0v1 == zero(KVector{1}, Float64)
    @test 1v1 ∧ 1.0v2 == 1v1 * 1.0v2 == 1.0v12
end


text:
"""
    epoch!(learner, phase[, dataiter])

Train `learner` for one epoch on `dataiter`. Iterates through
`dataiter` and [`step!`](#)s for each batch/item.

If no data iterator is passed in, use `learner.data[phasedataiter(phase)]`.

## Extending

The default implementation iterates over every batch in `dataiter`
and calls [`step!`](#) for each. This behavior can be overloaded
by implementing `epoch!(learner, ::MyPhase, dataiter)`.

If you're implementing a custom `epoch!` method, it is recommended
you make use of [`runepoch`](#) to get begin and end events as well
as proper handling of [`CancelEpochException`](#)s.

See the default implementation for reference.
"""
function epoch!(learner, phase::Phase, dataiter=learner.data[phasedataiter(phase)])
    runepoch(learner, phase) do _
        for batch in dataiter
            step!(learner, phase, batch)
        end
    end
end

"""
    step!(learner, phase::Phase, batch)

Run one step of training for `learner` on batch.
Behavior is customized through `phase`.

## Extending

This is a required method for custom [`Phase`](#)s to implement.
To implement `step!`, it is recommended you make use of [`runstep`](#)
to get begin and end events as well as proper handling of
[`CancelStepException`](#)s.

See the implementations of [`TrainingPhase`](#) and [`ValidationPhase`](#)
for reference.
"""
function step! end

function step!(learner, phase::TrainingPhase, batch)
    xs, ys = batch
    runstep(learner, phase, (; xs=xs, ys=ys)) do handle, state
        state.grads = gradient(learner.params) do
            state.ŷs = learner.model(state.xs)
            handle(LossBegin())
            state.loss = learner.lossfn(state.ŷs, state.ys)
            handle(BackwardBegin())
            return state.loss
        end
        handle(BackwardEnd())
        update!(learner.optimizer, learner.params, state.grads)
    end
end

function step!(learner, phase::ValidationPhase, batch)
    xs, ys = batch
    runstep(learner, phase, (;xs=xs, ys=ys)) do _, state
        state.ŷs = learner.model(state.xs)
        state.loss = learner.lossfn(state.ŷs, state.ys)
    end
end

"""
    runepoch(epochfn, learner, phase)

Run `epochfn` inside the context of an epoch. Calls `epochfn(handle)`
where `handle(e)` can be called to dispatch events.

Takes care of dispatching [`EpochBegin`](#) and [`EpochEnd`](#)
events as well as handling [`CancelEpochException`](#)s.

"""
function runepoch(epochfn, learner, phase::Phase)
    handlefn(e) = handle(learner.callbacks.runner, e, phase, learner)
    try
        handlefn(EpochBegin())
        epochfn(handlefn)
        handlefn(EpochEnd())
    catch e
        if e isa CancelEpochException
            @debug "Epoch skipped" error = e
            handlefn(EpochEnd())
        else
            rethrow()
        end
    end
end

"""
    runstep(stepfn, learner, phase) -> state

Run `stepfn` inside the context of a step. Calls `stepfn(handle, state)`
where `handle(e)` can be called to dispatch events and `state` is a [`PropDict`](#)
which step data, gradients and losses can be written to. Return `state`.

Takes care of dispatching [`StepBegin`](#) and [`StepEnd`](#)
events as well as handling [`CancelStepException`](#)s.
"""
function runstep(stepfn, learner, phase::Phase, initialstate = (;))
    state = PropDict(pairs(initialstate))
    handlefn(e) = handle(learner.callbacks.runner, e, phase, learner)
    try
        learner.step = state
        handlefn(StepBegin())
        stepfn(handlefn, state)
        handlefn(StepEnd())
        return state
    catch e
        if e isa CancelStepException
            @debug "Step skipped" error = e
        else
            rethrow()
        end
    end
    return state
end

# Utilities

"""
    fit!(learner, nepochs)
    fit!(learner, nepochs, (trainiter, validiter))

Train `learner` for `nepochs` of training and validation each. Use data
iterators that are passed in. If none are given, use `learner.data.training`
and `learner.data.validation`.

## Examples

```julia
fit!(learner, 10)
fit!(learner, 10, (traindl, valdl))
```

"""
function fit!(learner, nepochs::Int, (trainiter, validiter))
    for i in 1:nepochs
        epoch!(learner, TrainingPhase(), trainiter)
        epoch!(learner, ValidationPhase(), validiter)
    end
end

function fit!(learner, nepochs::Int)
    fit!(learner, nepochs, (learner.data.training, learner.data.validation))
end


text:
using ProbabilityModels, LinearAlgebra, DistributionParameters
sg = quote
    σ ~ Gamma(1.0, 0.05)
    L ~ LKJ(2.0)
    β ~ Normal(10.0) # μ = 0
    μ ~ Normal(100.0)

    σL = Diagonal(σ) * L
    Y₁ ~ Normal( X₁, β,       μ',  σL )
    Y₂ ~ Normal( X₂, β[1:7,:], μ', σL )
    # Y₁ ~ Normal( μ' .+ X₁ .*ˡ β,        σL )
    # Y₂ ~ Normal( μ' .+ X₂ .*ˡ β[1:7,:], σL )
end;
m = ProbabilityModels.read_model(sg, Main);

N, K₁, K₂, P = 100, 10, 7, 12;
μ = rand(P);
σU = rand(P, (3P)>>1) |> x -> cholesky(Hermitian(x * x')).U;
X₁ = rand(N, K₁);
X₂ = rand(N, K₂);
β = randn(K₁, P);
Y₁ = mul!(rand(N, P) * σU, X₁, β, 1.0, 1.0);
Y₂ = mul!(rand(N, P) * σU, X₂, view(β, 1:K₂, :), 1.0, 1.0);

datant = (
    Y₁ = Y₁, Y₂ = Y₂, X₁ = X₁, X₂ = X₂,
    μ = RealVector{12}(), β = RealMatrix{10,12}(),
    σ = RealVector{12,0.0,Inf}(), L = CorrelationMatrixCholesyFactor{12}()
);

ProbabilityModels.preprocess!(m, typeof(datant));

ProbabilityModels.ReverseDiffExpressions.lower(m)

dm = ProbabilityModels.ReverseDiffExpressions.differentiate(m);
ProbabilityModels.ReverseDiffExpressions.lower(dm)

@model SubGroups begin
    σ ~ Gamma(1.0, 0.05)
    L ~ LKJ(2.0)
    β ~ Normal(10.0) # μ = 0
    μ ~ Normal(100.0)

    σL = Diagonal(σ) * L
    Y₁ ~ Normal( X₁, β,       μ',  σL )
    Y₂ ~ Normal( X₂, β[1:7,:], μ', σL )
end

@model LogisticRegression begin
    β ~ Normal(10.0) # μ = 0
    y ~ BernoulliLogit( X * β )
end

@model SubGroups begin
    σ ~ Gamma(1.0, 0.05)
    L ~ LKJ(2.0)
    β ~ Normal(10.0) # μ = 0
    μ ~ Normal(100.0)

    σL = Diagonal(σ) * L
    Y₁ ~ Normal( X₁, β,       μ',  σL )
    Y₂ ~ Normal( X₂, β[1:7,:], μ', σL )
end

@model SubGroup begin
    σ ~ Gamma(1.0, 0.05)
    L ~ LKJ(2.0)
    β ~ Normal(10.0) # μ = 0
    μ ~ Normal(100.0)

    σL = Diagonal(σ) * L
    Y₁ ~ Normal( μ' .+ X₁ .*ˡ β,        σL )
    Y₂ ~ Normal( μ' .+ X₂ .*ˡ β[1:7,:], σL )
end

using PaddedMatrices, StructuredMatrices, DistributionParameters, LoopVectorization
using InplaceDHMC, VectorizedRNG
using Random, SpecialFunctions, MCMCChainSummaries, LinearAlgebra
using ProbabilityModels, LoopVectorization, SLEEFPirates, SIMDPirates, ProbabilityDistributions, PaddedMatrices
using ProbabilityModels: HierarchicalCentering, ∂HierarchicalCentering, ITPExpectedValue, ∂ITPExpectedValue
using DistributionParameters: CovarianceMatrix, MissingDataVector#, add
using PaddedMatrices: vexp
BLAS.set_num_threads(1)

@model ITPModel begin
    # Non-hierarchical Priors
    ρ ~ Beta(3, 1)
    lκ ~ lsgg(8.5, 1.5, 3.0, 1.5271796258079011) # μ = 1, σ² = 10
    σ ~ Gamma(1.5, 0.25) # μ = 6, σ² = 2.4
    θ ~ Normal(10)
    L ~ LKJ(2.0)
    # Hierarchical Priors.
    # h subscript, for highest in the hierarhcy.
    μₕ₁ ~ Normal(10) # μ = 0
    μₕ₂ ~ Normal(10) # μ = 0
    σₕ ~ Normal(10) # μ = 0
    # Raw μs; non-cenetered parameterization
    μᵣ₁ ~ Normal() # μ = 0, σ = 1
    μᵣ₂ ~ Normal() # μ = 0, σ = 1
    # Center the μs
    μᵦ₁ = HierarchicalCentering(μᵣ₁, μₕ₁, σₕ)
    μᵦ₂ = HierarchicalCentering(μᵣ₂, μₕ₂, σₕ)
    σᵦ ~ Normal(10) # μ = 0
    # Raw βs; non-cenetered parameterization
    βᵣ₁ ~ Normal()
    βᵣ₂ ~ Normal()
    # Center the βs.
    β₁ = HierarchicalCentering(βᵣ₁, μᵦ₁, σᵦ, domains)
    β₂ = HierarchicalCentering(βᵣ₂, μᵦ₂, σᵦ, domains)
    # Likelihood
    κ = vexp(lκ)
    μ₁ = vec(ITPExpectedValue(time, β₁, κ, θ))
    μ₂ = vec(ITPExpectedValue(time, β₂, κ, θ))
    Σ = CovarianceMatrix(ρ, Diagonal(σ) * L, time)
    # Tuple (Y₁, Y₂)
    (Y₁, Y₂) ~ Normal((μ₁, μ₂)[AvailableData], Σ[AvailableData])
end

function rinvscaledgamma(::Val{N},a::T,b::T,c::T) where {N,T}
    rg = MutableFixedSizeVector{N,T}(undef)
    log100 = log(100)
    @inbounds for n ∈ 1:N
        rg[n] = log100 / exp(log(VectorizedRNG.randgamma(a/c)) / c + b)
    end
    rg
end

const domains = ProbabilityModels.Domains(2,2,3)

const n_endpoints = sum(domains)

const times = MutableFixedSizeVector{36,Float64,36}(undef); times .= 0:35;

structured_missing_pattern = push!(vcat(([1,0,0,0,0] for i ∈ 1:7)...), 1);
missing_pattern = vcat(
    structured_missing_pattern, fill(1, 4length(times)), structured_missing_pattern, structured_missing_pattern
);

const availabledata = MissingDataVector{Float64}(missing_pattern);

const κ₀ = (8.5, 1.5, 3.0)

AR1(ρ, t) = @. ρ ^ abs(t - t')

function generate_true_parameters(domains, times, κ₀)
    K = sum(domains)
    D = length(domains)
    T = length(times)
    μ = MutableFixedSizeVector{K,Float64}(undef)
    offset = 0
    for i ∈ domains
        domain_mean = 5randn()
        for j ∈ 1+offset:i+offset
            μ[j] = domain_mean + 5randn()
        end
        offset += i
    end
    σ = VectorizedRNG.randgamma( 6.0, 1/6.0)
    ρ = VectorizedRNG.randbeta(4.0,4.0)
    κ = rinvscaledgamma(Val(K), κ₀...)
    lt = last(times)
    β = MutableFixedSizeVector{7,Float64,7}(( 0.0625,  0.0575,  0.0525,  0.0475,  0.0425,  0.04,  0.0375))
    θ₁ = @. μ' - β' * ( 1.0 - exp( - κ' * times) ) / (1.0 - exp( - κ' * lt) )
    θ₂ = @. μ' + β' * ( 1.0 - exp( - κ' * times) ) / (1.0 - exp( - κ' * lt) )

    L_T, info = LAPACK.potrf!('L', AR1(ρ, times))
    @inbounds for tc ∈ 2:T, tr ∈ 1:tc-1
        L_T[tr,tc] = 0.0
    end
    X = PaddedMatrices.MutableFixedSizeMatrix{K,K+3,Float64,K}(undef); randn!(X)
    U_K, info = LAPACK.potrf!('U', BLAS.syrk!('U', 'N', σ, X, 0.0, zero(MutableFixedSizeMatrix{K,K,Float64,K})))
    (
        U_K = U_K, L_T = L_T, μ = μ, θ₁ = θ₁, θ₂ = θ₂, domains = domains, time = times
    )
end

@generated function randomize!(
    sp::PaddedMatrices.StackPointer,
    A::AbstractArray{T,P},
    B::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,M,T},
    C::PaddedMatrices.AbstractMutableFixedSizeMatrix{N,N,T},
    D::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,T}
) where {M,N,T,P}
    quote
        @boundscheck begin
            d = size(A,1)
            for p ∈ 2:$(P-1)
                d *= size(A,p)
            end
            d == M*N || PaddedMatrices.ThrowBoundsError("Earlier dims size(A) == $(size(A)) does not match size(D) == ($M,$N)")
        end
        ptr = pointer(sp, $T)
        E = PtrMatrix{$M,$N,$T,$M}( ptr )
        F = PtrMatrix{$M,$N,$T,$M}( ptr + $(sizeof(T) * M * N) )
        ptr_A = pointer(A)
        GC.@preserve A begin
            for n ∈ 0:size(A,$P)-1
                Aₙ = PtrMatrix{$M,$N,$T,$M}( ptr_A + n*$(sizeof(T)*M*N) )
                randn!(ProbabilityModels.GLOBAL_PCGs[1], E)
                mul!(F, B, E)
                Aₙ .= D
                PaddedMatrices.gemm!(Aₙ, F, C)
            end
        end
        sp
    end
end

sample_data( N, truth, missingness, missingvals = (Val{0}(),Val{0}()) ) = sample_data( N, truth, missingness, missingvals, truth.domains )
@generated function sample_data(
    N::Tuple{Int,Int},
	truth, missingness,
	::Tuple{Val{M1},Val{M2}},
	::ProbabilityModels.Domains{S}
) where {S,M1,M2}
    K = sum(S)
    D = length(S)
    quote
        N₁, N₂ = N
        L_T = truth.L_T
        U_K = truth.U_K
        T = size(L_T,1)

        sp = ProbabilityModels.STACK_POINTER_REF[]
        (sp,Y₁) = PaddedMatrices.DynamicPtrArray{Float64,3}(sp, (T, $K, N₁), T)
        (sp,Y₂) = PaddedMatrices.DynamicPtrArray{Float64,3}(sp, (T, $K, N₂), T)
        randomize!(sp, Y₁, L_T, U_K, truth.θ₁)
        randomize!(sp, Y₂, L_T, U_K, truth.θ₂)

        c = length(missingness.indices)
        inds = missingness.indices

        Y₁sub = reshape(Y₁, (T * $K, N₁))[inds, :]
        $(M1 > 0 ? quote
            Y₁union = Array{Union{Missing,Float64}}(Y₁sub)
            perm = randperm(length(Y₁sub))
            @inbounds for m in 1:$M1
                Y₁union[perm[m]] = Base.missing
	    	end
    		Y₁ = convert(MissingDataArray{$M1,Bounds(-Inf,Inf)}, Y₁union)
		end : quote
			Y₁ = Y₁sub
		end)

        Y₂sub = reshape(Y₂, (T * $K, N₂))[inds, :]
		$(M2 > 0 ? quote
            Y₂union = Array{Union{Missing,Float64}}(Y₂sub)
            perm = randperm(length(Y₂sub))
			@inbounds for m in 1:$M2
                Y₂union[perm[m]] = Base.missing
            end
            Y₂ = convert(MissingDataArray{$M2,Bounds(-Inf,Inf)},Y₂union)
		end : quote
		    Y₂ = Y₂sub
		end)

        ITPModel(
            domains = truth.domains,
	        AvailableData = missingness,
            Y₁ = Y₁,
            Y₂ = Y₂,
            time = truth.time,
            L = LKJCorrCholesky{$K},
            ρ = RealVector{$K,Bounds(0,1)},
            lκ = RealVector{$K},
            θ = RealVector{$K},
            μₕ₁ = RealFloat,
            μₕ₂ = RealFloat,
            μᵣ₁ = RealVector{$D},
            μᵣ₂ = RealVector{$D},
            βᵣ₁ = RealVector{$K},
            βᵣ₂ = RealVector{$K},
            σᵦ = RealFloat{Bounds(0,Inf)},
            σₕ = RealFloat{Bounds(0,Inf)},
            σ = RealVector{$K,Bounds(0,Inf)}
        )
    end
end

truth = generate_true_parameters(domains, times, κ₀);

data = sample_data((100,100), truth, availabledata);
# mdata = sample_data((100,100), truth, availabledata, ((Val(10),Val(9))));

a = randn(length(data)); g = similar(a);

@time logdensity(data, a)
@time logdensity_and_gradient(data, a)

@time ProbabilityModels.check_gradient(data, a)

using BenchmarkTools

@benchmark logdensity($data, $a)
@benchmark logdensity_and_gradient!($g, $data, $a)

pg = PtrVector{73,Float64,73}(pointer(g));
pa = PtrVector{73,Float64,73}(pointer(a));
@code_warntype logdensity_and_gradient!(pg, data, pa, ProbabilityModels.STACK_POINTER_REF[])

using MacroTools, PaddedMatrices, DiffRules, VectorizationBase, SLEEFPirates
using MacroTools: striplines, @capture, prewalk, postwalk, @q

using LogDensityProblems, DynamicHMC, VectorizedRNG

includet("/home/chriselrod/.julia/dev/DistributionParameters/src/uniform_mapped_parameters.jl")
includet("/home/chriselrod/Documents/progwork/julia/autodiff_work/distribution_functions.jl")
includet("/home/chriselrod/Documents/progwork/julia/autodiff_work/special_diff_rules.jl")
includet("/home/chriselrod/Documents/progwork/julia/autodiff_work/model_macro_passes.jl")

# q = @q begin
#     β₀ ~ Normal(μ₀, σ₀)
#     β₁ ~ Normal(μ₁, σ₁)
#     y ~ Bernoulli_logit(β₀ + x * β₁)
# end
#
# q2 = translate_sampling_statements(q)
# q3 = flatten_expression(q2)

model_parameters = Set([:β₀,  :β₁])
ProbabilityModels.first_updates_to_assignemnts(q3, varnames)

q = @q begin
    β₀ ~ Normal(μ₀, σ₀)
    β₁ ~ Normal(μ₁, σ₁)
    y ~ Bernoulli_logit(β₀ + x * β₁)
end
model_name = :BernoulliLogitModel
variable_set = ProbabilityModels.determine_variables(q)
variables = [v for v ∈ variable_set] # ensure order is constant
variable_type_names = [Symbol("##Type##", v) for v ∈ variables]
struct_quote = quote
    struct $model_name{$(variable_type_names...)} <: LogDensityProblems.AbstractLogDensityProblem
        $([:( $(variables[i])::$(variable_type_names[i]) ) for i ∈ eachindex(variables)]...)
    end
end
struct_kwarg_quote = quote
    function $model_name(; $(variables...))
        $model_name($([:(ProbabilityModels.types_to_vals($v)) for v ∈ variables]...))
    end
end

# Translate the sampling statements, and then flatten the expression to remove nesting.
expr = ProbabilityModels.translate_sampling_statements(q) |>
            ProbabilityModels.flatten_expression
first_pass = quote end
second_pass = quote end
tracked_vars = Set(model_parameters)
first_pass, name_dict = ProbabilityModels.rename_assignments(first_pass)
expr, name_dict = ProbabilityModels.rename_assignments(expr, name_dict)

θ_sym = $(QuoteNode(θ)) # This creates our symbol θ

second_pass, name_dict = ProbabilityModels.rename_assignments(second_pass, name_dict)
ProbabilityModels.reverse_diff_pass!(first_pass, second_pass, expr, tracked_vars)
# variable renaming rather than incrementing makes initiazing
# target to an integer okay.
expr_out = quote
    # target = 0
    $first_pass
    $(Symbol("###seed###", name_dict[:target])) = ProbabilityModels.One()
    $second_pass
    (
        $(name_dict[:target]),
        $(Expr(:tuple, [Symbol("###seed###", mp) for mp ∈ model_parameters]...))
    )
end

using ProbabilityModels, MacroTools, PaddedMatrices, DiffRules, VectorizationBase, SLEEFPirates
using DistributionParameters
using MacroTools: striplines, @capture, prewalk, postwalk, @q

using LogDensityProblems, DynamicHMC, VectorizedRNG

# model_parameters = Set([:β₀,  :β₁])
# ProbabilityModels.first_updates_to_assignemnts(q3, varnames)

q = @q begin
    β₀ ~ Normal(μ₀, σ₀)
    β₁ ~ Normal(μ₁, σ₁)
    y ~ Bernoulli_logit(β₀ + x * β₁)
end

# variable_set = ProbabilityModels.determine_variables(expr)
# variables = [v for v ∈ variable_set] # ensure order is constant
# variable_type_names = [Symbol("##Type##", v) for v ∈ variables]

q2 = ProbabilityModels.translate_sampling_statements(q)
expr = ProbabilityModels.flatten_expression(q2)
return_partials = true
model_parameters = Symbol[]
first_pass = quote end
second_pass = quote end

push!(model_parameters, :β₀)
DistributionParameters.load_parameter(first_pass.args, second_pass.args, :β₀, RealFloat, return_partials)
push!(model_parameters, :β₁)
DistributionParameters.load_parameter(first_pass.args, second_pass.args, :β₁, RealVector{4,Float64}, return_partials)

tracked_vars = Set(model_parameters)
first_pass, name_dict = ProbabilityModels.rename_assignments(first_pass)
expr, name_dict = ProbabilityModels.rename_assignments(expr, name_dict)

second_pass, name_dict = ProbabilityModels.rename_assignments(second_pass, name_dict)
TLθ = 5 #type_length($θ) # This refers to the type of the input
ProbabilityModels.reverse_diff_pass!(first_pass, second_pass, expr, tracked_vars)
first_pass
second_pass

combined_q = quote
    $first_pass
    $(Symbol("###seed###", name_dict[:target])) = ProbabilityModels.One()
    $second_pass
end;

ProbabilityModels.first_updates_to_assignemnts(combined_q, model_parameters)

T_sym = gensym(:T); θ_sym = gensym(:θ);
expr_out = quote
    # target = zero($T_sym)
    $(Symbol("##θparameter##")) = VectorizationBase.vectorizable($θ_sym)
    $first_pass
    $(Symbol("##∂θparameter##m")) = PaddedMatrices.MutableFixedSizeVector{$TLθ,$T_sym}(undef)
    $(Symbol("##∂θparameter##")) = VectorizationBase.vectorizable($(Symbol("##∂θparameter##m")))
    $(Symbol("###adjoint###", name_dict[:target])) = ProbabilityModels.One()
    $second_pass
    LogDensityProblems.ValueGradient(
        $(name_dict[:target]),
        PaddedMatrices.ConstantFixedSizeVector($(Symbol("##∂θparameter##m")))
    )
end;
final_expr = first_updates_to_assignemnts(expr_out, model_parameters)

using MacroTools, PaddedMatrices, DiffRules, VectorizationBase, SLEEF, Random, SIMDPirates
using MacroTools: striplines, @capture, prewalk, postwalk, @q
using LogDensityProblems, DynamicHMC, VectorizedRNG, RandomNumbers

includet("/home/chriselrod/.julia/dev/DistributionParameters/src/uniform_mapped_parameters.jl")
includet("/home/chriselrod/Documents/progwork/julia/autodiff_work/distribution_functions.jl")
includet("/home/chriselrod/Documents/progwork/julia/autodiff_work/special_diff_rules.jl")
includet("/home/chriselrod/Documents/progwork/julia/autodiff_work/model_macro_passes.jl")
includet("/home/chriselrod/.julia/dev/ProbabilityModels/src/dynamic_hmc_interface.jl")

using ProbabilityModels, DistributionParameters, LogDensityProblems, DynamicHMC
using VectorizationBase, LoopVectorization
@model BernoulliLogitModel begin
    β₀ ~ Normal(μ₀, σ₀)
    β₁ ~ Normal(μ₁, σ₁)
    y ~ Bernoulli_logit(β₀ + x * β₁)
end

# using StaticArrays
N = 800; N_β = 4;
X = randn(N, N_β);
β₁ = [-1.595526740808615, -1.737875659746032, -0.26107993378119343, 0.6500851571519769];
β₀ = -0.05
Xβ₁ = X * β₁;
p = rand(N);
y = @. p < 1 / (1 + exp( - Xβ₁ - β₀));
sum(y)

ℓ = BernoulliLogitModel(
    σ₀ = 10.0, σ₁ = 5.0, μ₀ = 0.0, μ₁ = 0.0,
    β₀ = RealFloat, β₁ = RealVector{4},
    y = y, x = X
);

# eval(θquote)
dimension(ℓ)
a = fill(1.0, dimension(ℓ));
logdensity(LogDensityProblems.ValueGradient, ℓ, a)

using LogDensityProblems, DynamicHMC
@time mcmc_chain, tuned_sampler = NUTS_init_tune_mcmc_default(ℓ, 4000);
sample_mean(mcmc_chain)
sample_cov(mcmc_chain)

using BenchmarkTools
@benchmark logdensity(LogDensityProblems.ValueGradient, $ℓ, $a)
@benchmark NUTS_init_tune_mcmc_default($ℓ, 4000, report = DynamicHMC.ReportSilent())

struct_quote, struct_kwarg_quote, directquote, thetaquote, dimquote = generate_generated_funcs_expressions(:BernoulliLogitModel, q);
using LogDensityProblems, DynamicHMC
eval(struct_quote)
eval(struct_kwarg_quote)
eval(thetaquote)
eval(dimquote)
@time mcmc_chain, tuned_sampler = NUTS_init_tune_mcmc_default(rng, ℓ1, 4000);
sample_mean(mcmc_chain)
sample_cov(mcmc_chain)

@benchmark NUTS_init_tune_mcmc($rng, $ℓ1, 4000, report = DynamicHMC.ReportSilent())

mcmc_chainv = reshape(reinterpret(Float64, get_position.(mcmc_chain)), (8,length(mcmc_chain)))[1:5,:];

@code_warntype NUTS_init_tune_mcmc(rng, ℓ1, 4000, report = DynamicHMC.ReportSilent())
# Not type stable?

@code_warntype DynamicHMC.NUTS_init(rng, ℓ1, report = DynamicHMC.ReportSilent())
sampler_init = DynamicHMC.NUTS_init(rng, ℓ1, report = DynamicHMC.ReportSilent())
tuners = DynamicHMC.bracketed_doubling_tuner()
@code_warntype DynamicHMC.tune(sampler_init, tuners)
sampler_tuned = DynamicHMC.tune(sampler_init, DynamicHMC.bracketed_doubling_tuner());
@code_warntype DynamicHMC.mcmc(sampler_tuned, 4000)
mcmcres = DynamicHMC.mcmc(sampler_tuned, 4000)

function iterator_test(N)
    out = 0.0
    for i ∈ 1:4:N
        out += 2.3i
    end
    out
end

GK1 = GaussianKE(dimension(ℓ1));
H = DynamicHMC.Hamiltonian(ℓ1, GK1);
q = randn(rng, dimension(ℓ1));
p = rand(rng, GK1);
z = DynamicHMC.phasepoint_in(H, q, p)
DynamicHMC.find_initial_stepsize(DynamicHMC.InitialStepsizeSearch(), H, z)
@code_warntype DynamicHMC.NUTS_init(rng, ℓ1, report = false)

sampler_init = DynamicHMC.NUTS_init(rng, ℓ1);
tuners = DynamicHMC.bracketed_doubling_tuner()
@which DynamicHMC.tune(sampler_init,tuners.tuners[1])

sample, A = DynamicHMC.mcmc_adapting_ϵ(sampler_init, tuners.tuners[1].N)

function sample_sum(sample)
    N = length(sample)
    x̄ = DynamicHMC.get_position(sample[1])
    @inbounds for n ∈ 2:N
        x̄ += DynamicHMC.get_position(sample[n])
    end
    x̄
end
sample_v = get_position.(sample);
sample_mat = reshape(reinterpret(Float64, sample_v), (8,75))[1:5,:];

using ForwardDiff, Optim
inv_logit(x) = 1 / (1 + exp(-x))
function Bernoulli_logit_fmadd_simple(y, X, β, α)
    p = inv_logit.( X * β .+ α )
    T = eltype(p)
    target = zero(T)
    @inbounds @simd ivdep for i ∈ eachindex(p, y)
        OmP = one(T) - p[i]
        target += y[i] ? log(p[i]) : log(OmP)
    end
    target
end
opt = optimize(b -> -Bernoulli_logit_fmadd_simple(y, X, b, 1.0), ones(4), BFGS())

βones =  fill(1.0, PaddedMatrices.Static{4}())
∂Bernoulli_logit_fmadd_logeval_dropconst(y, X, βones, 1.0, Val{(false,false,true,true)}())

Bernoulli_logit_fmadd_simple(y, X, ones(4), 1.0)
ForwardDiff.derivative(a -> Bernoulli_logit_fmadd_simple(y, X, ones(4), a), 1.0)
ForwardDiff.gradient(b -> Bernoulli_logit_fmadd_simple(y, X, b, 1.0), ones(4))
Normal_logeval_dropconst(βones, 0.0, 5.0, Val{(true,false,false)}())

using ProbabilityModels, DistributionParameters, ProbabilityDistributions,
    LoopVectorization, LinearAlgebra,
    LogDensityProblems, SLEEFPirates, SIMDPirates, StructuredMatrices, ScatteredArrays, PaddedMatrices
using DistributionParameters: LKJ_Correlation_Cholesky, RealFloat, PositiveFloat, UnitFloat, RealVector, PositiveVector
using ProbabilityModels: HierarchicalCentering, ∂HierarchicalCentering, ITPExpectedValue, ∂ITPExpectedValue

@model ITPModel begin

    # Non-hierarchical Priors
    (0.5ρ + 0.5) ~ Beta(2, 2)
    κ ~ Gamma(0.1, 0.1) # μ = 1, σ² = 10
    σ ~ Gamma(1.5, 0.25) # μ = 6, σ² = 2.4
    θ ~ Normal(10)
    L ~ LKJ(2.0)

    # Hierarchical Priors.
    # h subscript, for highest in the hierarhcy.
    μₕ₁ ~ Normal(10) # μ = 0
    μₕ₂ ~ Normal(10) # μ = 0
    σₕ ~ Normal(10) # μ = 0
    # Raw μs; non-cenetered parameterization
    μᵣ₁ ~ Normal() # μ = 0, σ = 1
    μᵣ₂ ~ Normal() # μ = 0, σ = 1
    # Center the μs
    μᵦ₁ = HierarchicalCentering(μᵣ₁, μₕ₁, σₕ)
    μᵦ₂ = HierarchicalCentering(μᵣ₂, μₕ₂, σₕ)
    σᵦ ~ Normal(10) # μ = 0
    # Raw βs; non-cenetered parameterization
    βᵣ₁ ~ Normal()
    βᵣ₂ ~ Normal()
    # Center the βs.
    β₁ = HierarchicalCentering(βᵣ₁, μᵦ₁, σᵦ, domains)
    β₂ = HierarchicalCentering(βᵣ₂, μᵦ₂, σᵦ, domains)

    U = inv′(Diagonal(σ) * L)

    # Likelihood
    μ₁ = ITPExpectedValue(t, β₁, κ, θ)
    μ₂ = ITPExpectedValue(t, β₂, κ, θ)
    #AR = AutoregressiveMatrix(2ρ-1, δₜ)
    AR = AutoregressiveMatrix(ρ, δₜ)
    Y₁ ~ Normal(μ₁, AR, U)
    Y₂ ~ Normal(μ₂, AR, U)

end

# mub1 = HierarchicalCentering(xc.μᵣ₁, xc.μₕ₁, xc.σₕ);
# mub2 = HierarchicalCentering(xc.μᵣ₂, xc.μₕ₂, xc.σₕ);
# b1 = HierarchicalCentering(xc.βᵣ₁, mub1, xc.σᵦ, domains);
# b2 = HierarchicalCentering(xc.βᵣ₂, mub2, xc.σᵦ, domains);
# mu1 = ITPExpectedValue(times, b1, xc.κ, xc.θ);
# mu2 = ITPExpectedValue(times, b2, xc.κ, xc.θ);
#
# ProbabilityDistributions.Normal(Y₁, mu1, AutoregressiveMatrix(xc.ρ, δₜ), sLinv, Val((false,true,true,true)))
# ProbabilityDistributions.Normal(Y₂, mu2, AutoregressiveMatrix(xc.ρ, δₜ), sLinv, Val((false,true,true,true)))

#     Defined model: ITPModel.
#     Unknowns: Y₂, domains, μₕ₂, δt, μᵣ₁, βₖ, η₀, αₖ, σᵦ, θ, μᵣ₂, ρ, σₕ, μₕ₁, κ, L, Y₁, βᵣ₂, t, βᵣ₁.

using ProbabilityModels, DistributionParameters, ProbabilityDistributions,
    LoopVectorization, LinearAlgebra, Random, LogDensityProblems, SLEEFPirates, SIMDPirates, StructuredMatrices, ScatteredArrays, PaddedMatrices
using DistributionParameters: LKJ_Correlation_Cholesky, RealFloat, PositiveFloat, UnitFloat, RealVector, PositiveVector
using ProbabilityModels: HierarchicalCentering, ∂HierarchicalCentering, ITPExpectedValue, ∂ITPExpectedValue
BLAS.set_num_threads(1)

domains = ProbabilityModels.Domains(2,2,2,3);
# domains = ProbabilityModels.Domains(2,2,3);
T = 24; K = sum(domains); D = length(domains);

κ = (1/32) * reduce(+, (@Constant randexp(K)) for i ∈ 1:8) # κ ~ Gamma(8, 32)
σd = sum(@Constant randexp(4)) / 16 # σd ~ Gamma(4,16)
θ = 2.0 * (@Constant randn(K)) # θ ~ Normal(0,2)
S = (@Constant randn(K,4K)) |> x -> x * x'
S *= (1/16)
pS = StructuredMatrices.SymmetricMatrixL(S)
L = PaddedMatrices.chol(S); U = PaddedMatrices.invchol(pS)
muh1, muh2 = -3.0, 9.0
m01 = muh1 + (@Constant randn(D)); # placebo
m02 = muh2 + (@Constant randn(D)); #treatment
b1 = HierarchicalCentering((@Constant randn(K)), m01, σd, domains); # placebo
b2 = HierarchicalCentering((@Constant randn(K)), m02, σd, domains); # treatment

δₜ = 0.06125 *  reduce(+, (@Constant randexp(T-1)) for i ∈ 1:8);
t = vcat(zero(ConstantFixedSizeVector{1,Float64}), cumsum(δₜ));
mu1 = ProbabilityModels.ITPExpectedValue(t, b1, κ, θ);
mu2 = ProbabilityModels.ITPExpectedValue(t, b2, κ, θ);

ρ = 0.7; ARmat = StructuredMatrices.AutoregressiveMatrix(ρ, δₜ);
# ARcholinv = ConstantFixedSizeMatrix(ARmat);
ARchol = PaddedMatrices.chol(ConstantFixedSizeMatrix(ARmat));

pσU = MutableFixedSizeMatrix{K,K,Float64}(undef); pσU .= L'; cσU = ConstantFixedSizeMatrix(pσU);
Y1 = [ARchol * (@Constant randn(T, K)) * cσU + mu1 for n in 1:56];
Y2 = [ARchol * (@Constant randn(T, K)) * cσU + mu2 for n in 1:56];
Y1c = ChunkedArray(Y1);
Y2c = ChunkedArray(Y2);

mu1a = Array(mu1);
ARcholinva = inv(Array(ARchol)); Ua = Array(U);
-0.5*sum(A -> sum(abs2,ARcholinva * (Array(A) .- mu1a) * Ua), Y1)
ProbabilityDistributions.Normal(Y1c, mu1, ARmat, U, Val((false,true,false,false)))
ProbabilityDistributions.∂Normal(Y1c, mu1, ARmat, U, Val((false,true,false,false)))
ProbabilityDistributions.Normal(Y1c, mu1, ARmat, U, Val((false,true,true,true)))
ProbabilityDistributions.∂Normal(Y1c, mu1, ARmat, U, Val((false,true,true,true)))

using BenchmarkTools
@benchmark ProbabilityDistributions.Normal($Y1c, $mu1, $ARmat, $U, Val((false,true,true,true)))
@benchmark ProbabilityDistributions.∂Normal($Y1c, $mu1, $ARmat, $U, Val((false,true,true,true)))

Y1r = reshape(reinterpret(Float64, Y1), (T, K, length(Y1)));
Y1p = permutedims(Y1r, (2,3,1));
Y1p2 = Array{Float64,3}(undef, K,length(Y1),T);
Y1p3 = Array{Float64,3}(undef, K,length(Y1),T);
mu1a2 = reshape(Array(mu1)', (K,1,T));
ARmata = Bidiagonal(vcat(1.0,Array(ARmat.spacing.ρᵗ)), -1.0 .* ARmat.spacing.rinvOmρ²ᵗ .* ARmat.spacing.ρᵗ, :U)
function lpdfn(Y1p3::AbstractArray{Tf,3}, Y1p2::AbstractArray{Tf,2}, Y1p, mu1, ARmat, U) where {Tf}
    @inbounds Y1p3 .= Y1p .- mu1
    K,N,T = size(Y1p3)
    ρᵗ = ARmat.spacing.ρᵗ
    rinvOmρ²ᵗ = ARmat.spacing.rinvOmρ²ᵗ
    Y1p3_2 = reshape(Y1p3, (K,N*T))
    Y1p2_3 = reshape(Y1p2, (K,N,T))
    @inbounds Y1p2_3[:,:,1] .= Y1p3[:,:,1]
    Y1p3_t2 = reshape(Y1p3, (K*N,T))
    Y1p2_t2 = reshape(Y1p2, (K*N,T))
    @inbounds for t ∈ 1:T-1
        rt = ρᵗ[t]
        rit = rinvOmρ²ᵗ[t]
        @fastmath @simd ivdep for i ∈ 1:K*N
            Y1p2_t2[i,t+1] =  rit * (Y1p3_t2[i,t+1] - rt* Y1p3_t2[i,t])
        end
    end
    mul!(Y1p3_2, U', Y1p2)
    -0.5sum(abs2, Y1p3_2)
end
lpdfn(Y1p3, Y1p2, Y1p, mu1a2, ARmat, Ua)

Y1r = reshape(reinterpret(Float64, Y1), (T, K, length(Y1)));
Y1p = permutedims(Y1r, (1,3,2));
Y1p2 = Array{Float64,3}(undef, T,length(Y1),K);
Y1p3 = Array{Float64,3}(undef, T,length(Y1),K);
mu1a2 = reshape(Array(mu1), (T,1,K));
Ua = Array(U);
ARmata = Bidiagonal(vcat(1.0,Array(ARmat.spacing.rinvOmρ²ᵗ)), -1.0 .* ARmat.spacing.rinvOmρ²ᵗ .* ARmat.spacing.ρᵗ, :L)
function lpdfn2(Y1p3::AbstractArray{Tf,3}, Y1p2::AbstractArray{Tf,3}, Y1p, mu1, ARmat, U) where {Tf}
    T,N,K = size(Y1p3)
    @inbounds Y1p3 .= Y1p .- mu1
    mul!(reshape(Y1p2,(T,N*K)), ARmat, reshape(Y1p3,(T,N*K)))
    mul!(reshape(Y1p3,(T*N,K)), reshape(Y1p2,(T*N,K)), U)
    N*(T*logdet(UpperTriangular(U)) + K*sum(log, ARmat.dv)) - 0.5dot(Y1p3, Y1p3)
end
lpdfn2(Y1p3, Y1p2, Y1p, mu1a2, ARmata, Ua)

Y1c2 = ChunkedArray(Y1);
using VectorizationBase
function mul2!(Y1c2, Y1c, U)
    vY1c = VectorizationBase.vectorizable(Y1c)
    vY1c2 = VectorizationBase.vectorizable(Y1c2)
    V = Vec{8,Float64}
    for i ∈ 0:(length(Y1c)>>3)-1
        mul!(vload(V, vY1c2+8i), vload(V,vY1c + 8i), U)
    end
end
Y1pm = reshape(Y1p, (size(Y1p,1)*size(Y1p,2),size(Y1p,3)));
Y1pmUa = similar(Y1pm);
@benchmark mul!($Y1pmUa, $Y1pm, $Ua)
@benchmark mul2!($Y1c2, $Y1c, $U)

ℓ = ITPModel(
    domains = domains, Y₁ = Y1c, Y₂ = Y2c, t = t, δₜ = δₜ,
    L = LKJ_Correlation_Cholesky{K}, ρ = UnitFloat,
    κ = PositiveVector{K}, θ = RealVector{K},
    μₕ₁ = RealFloat, μₕ₂ = RealFloat,
    μᵣ₁ = RealVector{D}, μᵣ₂ = RealVector{D},
    βᵣ₁ = RealVector{K}, βᵣ₂ = RealVector{K},
    σᵦ = PositiveFloat, σₕ = PositiveFloat,
    σ = PositiveVector{K}
    #η₀ = 2.0, αₖ = 0.1, βₖ = 0.1
);
## domains = (2,2,3)
# Unknowns: Y₂, domains, μₕ₂, μᵣ₁, σ, σᵦ, θ, μᵣ₂, ρ, σₕ, μₕ₁, κ, L, δₜ, Y₁, βᵣ₂, t, βᵣ₁.
# D = 3, K = 7
# μₕ₂ = 1, μᵣ = 2:4, σ = 5:11, σᵦ = 12, θ = 13:19, μᵣ₂ = 20:22, ρ = 23, σₕ = 24, μₕ₁ = 25, κ = 26:32
# L = 33:53,  βᵣ₂ = 54:60, βᵣ₁ = 61:67

## domains = (2,2,2,3)
# Unknowns: Y₂, domains, μₕ₂, μᵣ₁, σ, σᵦ, θ, μᵣ₂, ρ, σₕ, μₕ₁, κ, L, δₜ, Y₁, βᵣ₂, t, βᵣ₁.
# D = 4, K = 9
# μₕ₂ = 1, μᵣ = 2:5, σ = 6:14, σᵦ = 15, θ = 16:24, μᵣ₂ = 25:28, ρ = 29, σₕ = 30, μₕ₁ = 31, κ = 32:40
# L = 41:76,  βᵣ₂ = 77:85, βᵣ₁ = 86:94
## 33 through 39: wrong

##

using BenchmarkTools
dimension(ℓ)
# a = fill(1.0, dimension(ℓ));
a = randn(dimension(ℓ));
logdensity(LogDensityProblems.Value, ℓ, a)
vg = logdensity(LogDensityProblems.ValueGradient, ℓ, a)
vgg = vg.gradient;

@benchmark logdensity(LogDensityProblems.ValueGradient, $ℓ, $a)

function gradi(ℓ, a, a2, i)
    step = sqrt(eps(a[i]))
    # step = cbrt(eps(a[i]))
    v1 = logdensity(LogDensityProblems.Value, ℓ, a)
    a2[i] += step
    v2 = logdensity(LogDensityProblems.Value, ℓ, a2)
    a2[i] = a[i]
    (v2.value - v1.value) / step
end
function gradicompare(ℓ, vgg, a, a2, i)
    d1, d2 = (gradi(ℓ, a, a2, i), vgg[i])
    d1, d2, round((d1-d2)/d1, digits = 4)
end
function compare_grads(ℓ, a)
    a2 = MutableFixedSizeVector(a);
    vgg = logdensity(LogDensityProblems.ValueGradient, ℓ, a).gradient
    for i ∈ 1:dimension(ℓ)
       # ℓ, vgg, a, a2
       @show i, gradicompare(ℓ, vgg, a, a2, i)
    end
end
compare_grads(ℓ, a)
compare_grads(ℓ, randn(dimension(ℓ)))

using LogDensityProblems, DynamicHMC, MCMCDiagnostics
# @time mcmc_chain, tuned_sampler = NUTS_init_tune_mcmc_default(ℓ, 500);#, max_depth = 15);#, ϵ = 0.007);#, δ = 0.99);
@time mcmc_chain, tuned_sampler = NUTS_init_tune_mcmc_default(ℓ, 2000, δ = 0.75);#, max_depth = 15);#, ϵ = 0.007);#, δ = 0.99);
NUTS_statistics(mcmc_chain)
tuned_sampler
chain_matrix = get_position_matrix(mcmc_chain);
[effective_sample_size(chain_matrix[:,i]) for i in 1:20]'

@time mcmc_chain2, tuned_sampler2 = NUTS_init_tune_mcmc_default(ℓ, 2000, δ = 0.75);#, max_depth = 15);#, ϵ = 0.007);#, δ = 0.99);
NUTS_statistics(mcmc_chain2)
tuned_sampler2
chain_matrix2 = get_position_matrix(mcmc_chain2);
[effective_sample_size(chain_matrix2[:,i]) for i in 1:20]'

@time mcmc_chain3, tuned_sampler3 = NUTS_init_tune_mcmc_default(ℓ, 2000, δ = 0.75);#, max_depth = 15);#, ϵ = 0.007);#, δ = 0.99);
NUTS_statistics(mcmc_chain3)
tuned_sampler3
chain_matrix3 = get_position_matrix(mcmc_chain3);
[effective_sample_size(chain_matrix3[:,i]) for i in 1:20]'

using VectorizationBase, PaddedMatrices, StructuredMatrices, BenchmarkTools
δ = MutableFixedSizeMatrix{T,K,NTuple{8,Core.VecElement{Float64}},T,(T*K)}(undef);
δU = MutableFixedSizeMatrix{T,K,NTuple{8,Core.VecElement{Float64}},T,(T*K)}(undef);
Yᵥ = VectorizationBase.vectorizable(Y1c);
remmask = Y1c.mask;
@benchmark PaddedMatrices.vload!($δ, $Yᵥ + 0)
@benchmark PaddedMatrices.diff!($δ, $mu1, $δ)
@benchmark PaddedMatrices.mask!($δ, $remmask)
@benchmark mul!($δU, $δ, $U)

function chain_to_array(chain::Array{NUTS_Transition{ConstantFixedSizeArray{Tuple{N},Float64,1,Tuple{1},L},Float64},1}) where {N,L}
    pos = get_position.(chain)
    s = sizeof(eltype(pos)) >> 3
    reshape(reinterpret(Float64, pos), (s,length(pos)))[1:N,:]
end
function process_chain(chain)
    copy(chain_to_array(chain)')
end
c = process_chain(mcmc_chain);
ntuple(i -> effective_sample_size(c[:,i]), Val(16))
NUTS_statistics(mcmc_chain)
tuned_sampler

@time mcmc_chains = ProbabilityModels.NUTS_init_tune_threaded(ℓ, 1000, max_depth = 12);

@time mcmc_chains = ProbabilityModels.NUTS_init_tune_threaded(ℓ, 2000, max_depth = 12);

@time mcmc_chain, tuned_sampler = NUTS_init_tune_mcmc_default(ℓ, 1000, max_depth = 12);

report = ReportIO(countΔ = -1, time_nsΔ = -1);
@time mcmc_chain, tuned_sampler = NUTS_init_tune_mcmc_default(ℓ, 4000, report = report);

c1 = process_chain(mcmc_chains[2]);
# for i ∈ eachindex(mcmc_chains)
effective_sample_size(c1[:,1])
effective_sample_size(c1[:,2])
effective_sample_size(c1[:,3])
effective_sample_size(c1[:,4])

μₕ₁ = @Constant randn(4); # placebo
μₕ₂ = 2.0 + @Constant randn(4); #treatment
β₁ = vcat(ntuple(i -> (@Constant randn(domains[i])) + μₕ₁[i], Val(4))...); # placebo
β₂ = vcat(ntuple(i -> (@Constant randn(domains[i])) + μₕ₂[i], Val(4))...); # treatment

δt = @Constant rand(23); lastt = Ref(0.0);
t = ConstantFixedSizeVector{24,Float64}(ntuple(i -> i == 1 ? 0.0 : lastt[] += δt[i-1], Val(24)));

μ₁ = ITPExpectedValue(t, β₁, κ)
μ₂ = ITPExpectedValue(t, β₂, κ)

β = @CFixedSize randn(16);
κ = @CFixedSize randexp(16);

μ = ProbabilityModels.ITPExpectedValue(t, β, κ)

ℓitp = BernoulliLogitModel(
    α_κ = 1.0, δt = δtv, μ_β, η₀, σ_β, Y, β, ρ, α_ρ, κ, L, β_κ, t, β_ρ
);

# y =


text:
# Logistic Equation (test incomplete)

# Packages and inclusions
using ModelingToolkit,DiffEqOperators,DiffEqBase,LinearAlgebra,Test

# Tests
@testset "Test 00: Dt(u(t)) ~ u(t)*(1.0-u(t))" begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ 0.0*Dx(u(t,x))+u(t,x)*(1.0-u(t,x))
    bcs = [ u(0,x) ~ 0.5+x*0.0,
            u(t,0) ~ 0.5,
            u(t,1) ~ 0.5]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,5.0),
               x ∈ IntervalDomain(0.0,1.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])

    # Method of lines discretization
    dx = 0.1
    order = 1
    discretization = MOLFiniteDifference(dx,order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    using OrdinaryDiffEq

    sol = solve(prob,Euler(),dt=0.01,saveat=0.1)

    # Plot and save results
    #using Plots
    #time = domains[1].domain.lower:0.1:domains[1].domain.upper

    #plot(time,sol[5,1,:])
    #savefig("MOL_0D_Logistic.png")

    # Test
    # x_interval = domains[2].domain.lower+dx:dx:domains[2].domain.upper-dx
    # u = @. (0.5/(0.2*sqrt(2.0*3.1415)))*exp(-(x_interval-(0.75+0.6))^2/(2.0*0.2^2))
    # t_f = size(sol,3)
    # @test sol[t_f] ≈ u atol = 0.1;

end


text:
module test_julia_typealiases

using Test

struct G
end

const GVec = Vector{G}

if VERSION >= v"1.4"
    @test endswith(sprint(show, GVec([G()])), ".test_julia_typealiases.G()]")
elseif VERSION >= v"1.3"
    @test endswith(sprint(show, GVec([G()])), ".test_julia_typealiases.G()]")
end

function Base.show(io::IO, v::GVec)
    print(io, "GVec", '(')
    Base.show_delim_array(io, v, '[', ",", ']', false)
    print(io, ')')
end

if VERSION >= v"1.2"
    @test endswith(sprint(show, GVec([G()])), ".test_julia_typealiases.G()])")
else
    @test sprint(show, GVec([G()])) == "GVec([G()])"
end

const A = Vector{Int64}
if VERSION >= v"1.6"
    @test sprint(io -> show(io, MIME"text/plain"(), A)) == "Vector{Int64} (alias for Array{Int64, 1})"
end

end # module test_julia_typealiases


text:
function validate(obj::ZulipRequest, opts)
    obj.data = strip(obj.data)
    if obj.message.sender_id < 0 || isempty(obj.data) || isempty(obj.token)
        return false, "Wrong message, contact bot maintainer"
    end
    if obj.token != opts.token
        return false, "Incorrect token, verify ReminderBot server configuration"
    end

    if obj.message.type == "stream"
        if !startswith(obj.data, "@**")
            return false, ""
        end
    end
    if obj.data[1] == '@'
        m = match(r"^@\*\*[^\s]+\s+(.*)$"ms, obj.data)
        m === nothing && return false, "Wrong message. Refer to `help` on the usage of the Reminder Bot."
        obj.data = m[1]
    end

    return true, ""
end

########################################
# Parse time setup
########################################

function zparse_mehere(msg)
    rexp = r"^(me|here)\s+(.*)"msi
    m = match(rexp, msg)
    m === nothing && return :me, msg
    if lowercase(m[1]) == "me"
        return :me, m[2]
    else
        return :here, m[2]
    end
end

function zparse_relative(msg, exects)
    rexps = [
             r"^\s*([+-]?[0-9]+)\s+(months?)\s*(.*)"msi,
             r"^\s*([+-]?[0-9]+)\s+(weeks?)\s*(.*)"msi,
             r"^\s*([+-]?[0-9]+)\s+(days?)\s*(.*)"msi,
             r"^\s*([+-]?[0-9]+)\s+(hours?)\s*(.*)"msi,
             r"^\s*([+-]?[0-9]+)\s+(minutes?|min)\s*(.*)"msi,
             r"^\s*([+-]?[0-9]+)\s+(seconds?|sec)\s*(.*)"msi,
            ]
    matched = false
    while true
        ended = true
        for (i, r) in pairs(rexps)
            m = match(r, msg)
            m === nothing && continue
            matched = true
            ended = false
            msg = m[3]
            delta = parse(Int, m[1])
            if i == 1
                exects += Month(delta)
            elseif i == 2
                exects += Week(delta)
            elseif i == 3
                exects += Day(delta)
            elseif i == 4
                exects += Hour(delta)
            elseif i == 5
                exects += Minute(delta)
            else
                exects += Second(delta)
            end
        end
        ended && break
    end

    return matched, msg, exects
end

function zparse_absolute(msg, exects)
    rexpdt = r"^\s*([0-9]{4}-[0-9]{2}-[0-9]{2})(.*)"ms
    m = match(rexpdt, msg)
    m === nothing && return false, msg, exects
    msg = m[2]
    try
        exects = DateTime(m[1])
    catch
        return false, msg, exects
    end
    if msg[1] == '\n'
        msg = strip(msg)
        return true, msg, exects
    end
    if msg[1] == ' ' || msg[1] == 'T'
        msg = msg[2:end]
    else
        return false, msg, exects
    end

    rexptm = r"^([0-2][0-9]):([0-5][0-9]):([0-5][0-9])\s+(.*)"ms
    m = match(rexptm, msg)
    if m !== nothing
        hh = parse(Int, m[1])
        0 <= hh <= 23 || return false, msg, exects
        mm = parse(Int, m[2])
        0 <= mm <= 59 || return false, msg, exects
        ss = parse(Int, m[3])
        0 <= ss <= 59 || return false, msg, exects
        exects += Hour(hh) + Minute(mm) + Second(ss)
        return true, m[4], exects
    end

    rexptm = r"^([0-2][0-9]):([0-5][0-9])\s+(.*)"ms
    m = match(rexptm, msg)
    if m !== nothing
        hh = parse(Int, m[1])
        0 <= hh <= 23 || return false, msg, exects
        mm = parse(Int, m[2])
        0 <= mm <= 59 || return false, msg, exects
        exects += Hour(hh) + Minute(mm)
        return true, m[3], exects
    end

    rexptm = r"^([0-2][0-9])\s+(.*)"ms
    m = match(rexptm, msg)
    if m !== nothing
        hh = parse(Int, m[1])
        0 <= hh <= 23 || return false, msg, exects
        exects += Hour(hh)
        return true, m[2], exects
    end

    return true, strip(msg), exects
end

function zparse(msg, exects = Dates.now())
    gde, msg = zparse_mehere(msg)
    matched, msg, exects = zparse_relative(msg, exects)
    matched && return gde, :relative, msg, exects
    matched, msg, exects = zparse_absolute(msg, exects)
    matched && return gde, :absolute, msg, exects

    return gde, :unknown, msg, exects
end

########################################
# Narrow
########################################

function narrow(msg)
    return ZulipOpts[].baseep * "/#narrow/stream/$(msg.stream_id)-$(HTTP.escapeuri(msg.display_recipient))/topic/$(HTTP.escapeuri(msg.subject))/near/$(msg.id)"
end

function startswithnarrow(txt)
    contains(txt, r"^https?://[^\s]+/#narrow/stream/[^\s]+/near/[0-9]+")
end


text:
# Point in 1D space

export Point1D

mutable struct Point1D <: Point
    x :: Float64
end

point(x) = Point1D(x)
zero(::Type{Point1D}) = Point1D(0.0)
rand(r::AbstractRNG, ::Random.SamplerType{Point1D}) = Point1D(rand(r, [1.0, -1.0]))

+(lhs :: Point1D, rhs :: Point1D) = Point1D(lhs.x + rhs.x)
-(lhs :: Point1D, rhs :: Point1D) = Point1D(lhs.x - rhs.x)
-(p :: Point1D) = Point1D(-p.x)
*(a :: Real, p :: Point1D) = Point1D(a*p.x)
/(p :: Point1D, a :: Real) = Point1D(p.x/a)

abs(p :: Point1D) = abs(p.x)
abs2(p :: Point1D) = p.x*p.x

plot_str(p::Point1D) = string(p.x)

max_stars_repo_name: JuliaTagBot/JuliaBook-Samples

max_issues_repo_name: JuliaTagBot/JuliaBook-Samples

max_forks_repo_name: JuliaTagBot/JuliaBook-Samples


text:
#####
##### getindex
#####

function frule((_, ẋ), ::typeof(getindex), x::AbstractArray, inds...)
    return x[inds...], ẋ[inds...]
end

function rrule(::typeof(getindex), x::Array{<:Number}, inds...)
    # removes any logical indexing, CartesianIndex etc
    # leaving us just with a tuple of Int, Arrays of Int and Ranges of Int
    plain_inds = Base.to_indices(x, inds)
    y = getindex(x, plain_inds...)
    function getindex_pullback(ȳ)
        function getindex_add!(Δ)
            # this a optimizes away for simple cases
            for (ȳ_ii, ii) in zip(ȳ, Iterators.product(plain_inds...))
                Δ[ii...] += ȳ_ii
            end
            return Δ
        end

        x̄ = InplaceableThunk(
            getindex_add!,
            @thunk(getindex_add!(zero(x))),
        )
        īnds = broadcast(Returns(NoTangent()), inds)
        return (NoTangent(), x̄, īnds...)
    end

    return y, getindex_pullback
end

#####
##### view
#####

function frule((_, ẋ), ::typeof(view), x::AbstractArray, inds...)
    return view(x, inds...), view(ẋ, inds...)
end

#####
##### setindex!
#####

function frule((_, ẋ, v̇), ::typeof(setindex!), x::AbstractArray, v, inds...)
    return setindex!(x, v, inds...), setindex!(ẋ, v̇, inds...)
end

#####
##### `eachslice` and friends
#####

function rrule(::typeof(eachrow), x::AbstractVecOrMat)
    allrows(dy) = (NoTangent(), ∇eachslice(unthunk(dy), x, Val(1)))
    return collect(eachrow(x)), allrows
end

function rrule(::typeof(eachcol), x::AbstractVecOrMat)
    allcols(dy) = (NoTangent(), ∇eachslice(unthunk(dy), x, Val(2)))
    return collect(eachcol(x)), allcols
end

function rrule(::typeof(eachslice), x::AbstractArray; dims)
    y = collect(eachslice(x; dims=dims))
    @assert length(dims) == 1 """That's amazing, after many years JuliaLang/julia#32310
        actually landed. Sadly, the gradient rule for `eachslice` is unable to handle this
        case right now, please make an issue at https://github.com/JuliaDiff/ChainRules.jl"""
    dim = only(dims)
    allslices(dy) = (NoTangent(), ∇eachslice(unthunk(dy), x, Val(dim)))
    return y, allslices
end

# Using Val(dim) here is worth a factor of 2 in this, on Julia 1.8-
# @btime rrule(eachcol, $([1 2; 3 4]))[2]($([[10, 20], [30, 40]]))
function ∇eachslice(dys_raw, x::AbstractArray, vd::Val{dim}) where {dim}
    dys = unthunk(dys_raw)
    i1 = findfirst(dy -> dy isa AbstractArray, dys)
    if i1 === nothing  # all slices are Zero!
        return _zero_fill!(similar(x, float(eltype(x)), axes(x)))
    end
    T = promote_type(eltype(dys[i1]), eltype(x))
    # The whole point of this gradient is that we can allocate one `dx` array:
    dx = similar(x, T, axes(x))
    for i in axes(x, dim)
        slice = selectdim(dx, dim, i)
        if dys[i] isa AbstractZero
            _zero_fill!(slice)  # Avoids this: copyto!([1,2,3], ZeroTangent()) == [0,2,3]
        else
            copyto!(slice, dys[i])
        end
    end
    return ProjectTo(x)(dx)
end

_zero_fill!(dx::AbstractArray{<:Number}) = fill!(dx, zero(eltype(dx)))
_zero_fill!(dx::AbstractArray) = map!(zero, dx, dx)

function rrule(::typeof(∇eachslice), dys, x, vd::Val)
    function ∇∇eachslice(dz_raw)
        dz = unthunk(dz_raw)
        # eachslice(dz; dims=_val(vd)) does not make @code_warntype happy
        iter = vd == Val(1) ? eachrow(dz) : vd == Val(2) ? eachcol(dz) : eachslice(dz; dims=_val(vd))
        return (NoTangent(), collect(iter), NoTangent(), NoTangent())
    end
    return ∇eachslice(dys, x, vd), ∇∇eachslice
end

# These rules help with testing, and won't hurt:
# They are correct as we always `collect` the primal result as we need that
# information for the reverse pass
ChainRules.rrule(::typeof(collect∘eachrow), x) = rrule(eachrow, x)
ChainRules.rrule(::typeof(collect∘eachcol), x) = rrule(eachcol, x)
ChainRules.rrule(::typeof(collect∘eachslice), x; dims) = rrule(eachslice, x; dims=dims)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/solving-difference-equation-part-2/67057

# %%
using DynamicalSystems, DelimitedFiles, .Threads,  BenchmarkTools
## Components of a test DiscreteDynamicalSystem

function dds_constructor(u0 = [0.5, 0.7];  r=1.0, k=2.0)
    return DiscreteDynamicalSystem(dds_rule, u0, [r, k], dds_jac)
end
## equations of motion:
function dds_rule(x, par, n)
   r, k = par
   a, mu, d = 5.0, 0.5, 0.2
   dx = x[1]*exp(r*(1-x[1]/k)-x[2]/(a+x[1]^2))
   dy = x[2]*exp(mu*x[1]/(a+x[1]^2)-d)
   return @SVector [dx, dy]
end
## Jacobian:
function dds_jac( x, par, n)
    r, k = par;
    a = 5.0; mu = 0.5; d = 0.2;
    J11 = exp(- x[2]/(x[1]^2 + a) - r*(x[1]/k - 1)) - x[1]*exp(- x[2]/(x[1]^2 + a) - r*(x[1]/k - 1))*(r/k - (2*x[1]*x[2])/(x[1]^2 + a)^2)
    J12 = -(x[1]*exp(- x[2]/(x[1]^2 + a) - r*(x[1]/k - 1)))/(x[1]^2 + a)
    J21 = x[2]*exp((mu*x[1])/(x[1]^2 + a) - d)*(mu/(x[1]^2 + a) - (2*mu*x[1]^2)/(x[1]^2 + a)^2)
    J22 = exp((mu*x[1])/(x[1]^2 + a) - d)
    return @SMatrix [J11  J12; J21 J22]
end

# %%
function seqper_new(x; tol=1e-3)       # function to calculate periodicity
    @inbounds for k in 2:length(x)
        if abs(x[k] - x[1]) ≤ tol
            all(j -> abs(x[j] - x[j-k+1]) ≤ tol, k:lastindex(x)) && return k - 1
        end
    end
    return length(x)
end

# %%
function meshh(x,y)    # create a meshgrid of x and y
    len_x = length(x);
    len_y = length(y);
    xmesh = Float64[]
    ymesh = Float64[]
    for i in 1:len_x
        for j in 1:len_y
            push!(xmesh, x[i])
            push!(ymesh, y[j])
        end
    end
    return  xmesh, ymesh
end

# %%
#function isoperiodic_test(ds, par_area::Matrix{Float64}, NIter::Int64, filename::String,
function isoperiodic_test(ds, par_area::Matrix{Float64}, NIter::Int64, iso_filename::String,
                         nxblock::Int64, nyblock::Int64, NTr::Int64, xpts::Int64, ypts::Int64)

    io = open(iso_filename, "a")   # opeing file to write output
    p1_st =  par_area[1];          # beginning of first parameter
    p1_nd =  par_area[2];          # end of first parameter
    p2_st =  par_area[3];          # beginning of second parameter
    p2_nd =  par_area[4];          # end of second parameter

    p1_block = range(p1_st, p1_nd, length=nxblock + 1);    # divide par_area in blocks
    p2_block = range(p2_st, p2_nd, length=nyblock + 1);    # divide par_area in blocks
    l_bipar = xpts * ypts;                                 # total number of pts in each block
    total_blocks = nxblock * nyblock                       # total number of blocks
    sol_last = Array{Float64,2}(undef, NIter + 1, l_bipar);   # stores x values of solution for l_bipar (r,k) paris

    for ii in 1:nxblock
         step1size = (p1_block[ii + 1] - p1_block[ii]) / xpts;
         par1range = range(p1_block[ii], p1_block[ii + 1] - step1size, length=xpts);
        for jj in 1:nyblock
            #display("Progress: " * string(nxblock * (ii - 1) + jj) * " out of " * string(total_blocks) * " steps.")
            display("Progress: " * string(nyblock * (ii - 1) + jj) * " out of " * string(total_blocks) * " steps.")
            step2size = (p2_block[jj + 1] - p2_block[jj]) / ypts;
            par2range = range(p2_block[jj], p2_block[jj + 1] - step2size, length=ypts);
            Threads.@threads for i = 1:xpts
                for j = 1:ypts
                     set_parameter!(ds, [par1range[i], par2range[j]]) # change the parameter values <-- thread unsafe! bug!
                     tr = trajectory(ds, NIter; Ttr=NTr);             # find the solution
                     #sol_last[:,(i - 1) * xpts + j] =  tr[:,1];       # store x-values
                     sol_last[:,(i - 1) * ypts + j] =  tr[:,1];       # store x-values
                end
            end
            periods = seqper_new.(eachcol(sol_last), tol=0.001)       # seqper_new calculates periodicity
            par1mesh, par2mesh = meshh(par1range, par2range)          # meshgrid of parameter values
            writedlm(io, [par1mesh par2mesh periods])                 # write data in txt file
        end
    end
    close(io)
    return iso_filename
end

# %%
##  Test

parameter_area = [1.0 5.0 2.0 5.0]
nxblock = 1;
nyblock = 1;
xpts = 100;
ypts = 100;
NIter = 2000;
NTr = 50000;
init = [0.4, 0.5];
ds = dds_constructor(init);
##
iso_filename = "isoperiodic_test_data.txt"
isfile(iso_filename) && rm(iso_filename)
@time isoperiodic_test(ds, parameter_area, NIter, iso_filename, nxblock, nyblock, NTr, xpts, ypts)

# %%
##  Test

parameter_area = [1.0 5.0 2.0 5.0]
nxblock = 2;
nyblock = 2;
xpts = 100;
ypts = 100;
NIter = 2000;
NTr = 50000;
init = [0.4, 0.5];
ds = dds_constructor(init);
##
iso_filename = "isoperiodic_test_data2x2.txt"
isfile(iso_filename) && rm(iso_filename)
@time isoperiodic_test(ds, parameter_area, NIter, iso_filename, nxblock, nyblock, NTr, xpts, ypts)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public