text:
function defaultfcn(z)
    f = 1e9
    ϵᵣ = 5 - 2im
    μᵣ = 1 - 2im
    d = 1e-2
    c = 3e8
    ω = 2π*f
    k₀ = ω/c
    cc = ϵᵣ^2*(k₀*d)^2*(ϵᵣ*μᵣ - 1)
    w = ϵᵣ^2*z^2 + z^2*tan(z)^2 - cc
end

# Analysis parameters
xb = -2  # real part begin
xe = 2  # real part end
yb = -2  # imag part begin
ye = 2  # imag part end
r = 0.1999  # initial mesh step

origcoords = rectangulardomain(complex(xb, yb), complex(xe, ye), r)

# matlab results from https://github.com/PioKow/GRPF for comparison
matlab_zroots = [-1.624715288135189 + 0.182095877702038im,
                 -1.520192978034417 - 0.173670452237129im,
                 -0.515113098919392 + 0.507111597359180im,
                  0.515113098795215 - 0.507111597284675im,
                  1.520192978034417 + 0.173670452237129im,
                  1.624715288135189 - 0.182095877702037im]

matlab_zpoles = [-1.570796326699632 - 0.000000000206961im,
                  1.570796326699632 - 0.000000000206961im]

zroots, zpoles = grpf(defaultfcn, origcoords)

@test length(zroots) == 6
@test length(zpoles) == 2

@test approxmatch(zroots, matlab_zroots)
@test approxmatch(zpoles, matlab_zpoles)

pzroots, pzpoles, quadrants, phasediffs, tess, g2f = grpf(defaultfcn, origcoords, PlotData())

@test approxmatch(pzroots, matlab_zroots)
@test approxmatch(pzpoles, matlab_zpoles)

tmpr, tmpp = grpf(defaultfcn, origcoords, GRPFParams(8000, 1e-9))

@test approxmatch(tmpr, zroots)
@test approxmatch(tmpp, zpoles)

tmpr, tmpp, quadrants, phasediffs, tess, g2f = grpf(defaultfcn, origcoords, PlotData(), GRPFParams(8000, 1e-9))

@test approxmatch(tmpr, zroots)
@test approxmatch(tmpp, zpoles)

tmpr, tmpp = grpf(defaultfcn, origcoords, GRPFParams(8000, 1e-9, true))

@test approxmatch(tmpr, zroots)
@test approxmatch(tmpp, zpoles)

# Test with very low maxnodes
maxnodes = 10
params = GRPFParams(100, maxnodes, 3, maxnodes-1, 1e-9, false)
@test_logs (:warn,"GRPFParams.maxnodes reached") grpf(defaultfcn, origcoords, PlotData(), params)

# Test with big origcoords
xb = big(xb)  # real part begin
xe = big(xe)  # real part end
yb = big(yb)  # imag part begin
ye = big(ye)  # imag part end
r = big(r)  # initial mesh step

origcoords = rectangulardomain(complex(xb, yb), complex(xe, ye), r)

bzroots, bzpoles = grpf(defaultfcn, origcoords)

@test all(isa.(bzroots, Complex{BigFloat}))
@test all(isa.(bzpoles, Complex{BigFloat}))

@test approxmatch(bzroots, zroots)
@test approxmatch(bzpoles, zpoles)

# Test region with no roots or poles
xb, xe = -1, 1
yb, ye = -2, -1
origcoords = rectangulardomain(complex(xb, yb), complex(xe, ye), r)

nozroots, nozpoles = grpf(defaultfcn, origcoords)

@test length(nozroots) == 0
@test length(nozpoles) == 0

nozroots, nozpoles, quadrants, phasediffs, tess, g2f = grpf(defaultfcn, origcoords, PlotData())

@test length(nozroots) == 0
@test length(nozpoles) == 0

xb, xe = big"-1", big"1"
yb, ye = big"-2", big"-1"
r = big"0.2"

origcoords = rectangulardomain(complex(xb, yb), complex(xe, ye), r)

bnozroots, bnozpoles = grpf(defaultfcn, origcoords)

@test eltype(bnozroots) == Complex{BigFloat}
@test eltype(bnozpoles) == Complex{BigFloat}

@test length(bnozroots) == 0
@test length(bnozpoles) == 0


text:
#=
    Figure 7
=#

using Inference
using Model

using JLD2
using KernelDensity
using Plots
using StatsFuns
using StatsPlots
using Random

include("../figure_defaults.jl")

##############################################################
## Load model and results
##############################################################

include("../../Results/MainResult_Setup.jl")
@load "Results/MainResult.jld2"

include("../../Results/DeterministicModel.jl")

##############################################################
## Fig 6ab: Plot statistics at best fit (log scale)
##############################################################

function Wpdf_fraction(θ,t)
    d = param_dist(θ[5:end-1])
    ξ = sample(d,100_000)
    Z = parameters_to_antibody(t,ξ[2,:],ξ[3,:],θ[end])
    W = Z[2] ./ Z[1]
    return x -> pdf(kde(W),x)
end
function Wpdf_amount(θ,t)
    d = param_dist(θ[5:end-1])
    ξ = sample(d,100_000)
    Z = parameters_to_antibody(t,ξ[2,:],ξ[3,:],θ[end])
    W = Z[2]
    return x -> pdf(kde(W),x)
end

# Plot fraction internalised
lt = range(log10(0.01),log10(200.0),length=500)
x  = range(0.0,1.0,length=200)
f  = hcat([Wpdf_fraction(θ,tᵢ)(x) for tᵢ in 10.0.^(lt)]...)

# Normalise pdf by mode
for i = 1:length(lt)
    f[:,i] /= maximum(f[:,i])
end

# Plot amount internalised
lt = range(log10(0.01),log10(200.0),length=500)
x2  = range(0.0,2.2,length=200)
f2  = hcat([Wpdf2(θ,tᵢ)(x2) for tᵢ in 10.0.^(lt)]...)

# Normalise pdf by mode
for i = 1:length(lt)
    f2[:,i] /= maximum(f2[:,i])
end

# Plot
fig6a = contourf(lt,x,f,c=:RdPu_9,lw=0.0,ylabel="Fraction internal")
fig6b = contourf(lt,x2,f2,c=:RdPu_9,lw=0.0,ylabel="Relative amount internalised",ylim=(0.0,2.2))

# Add deterministic solution
x₁ = [solve_ode_model(t,collect(θdet))[2] / solve_ode_model(t,collect(θdet))[1] for t in 10.0.^(lt)]
x₂ = [solve_ode_model(t,collect(θdet))[2]  for t in 10.0.^(lt)]

plot!(fig6a,lt,x₁,lw=3.0,c=:black,label="")
plot!(fig6a,lt,x₁,lw=2.0,c=:white,label="")
plot!(fig6b,lt,x₂,lw=3.0,c=:black,label="")
plot!(fig6b,lt,x₂,lw=2.0,c=:white,label="")

fig6ab = plot(fig6a,fig6b,xticks=(-2:1:2,10.0.^(-2:1:3)),colorbar_title="Probability density")
plot!(fig6ab,size=(800,200),axis=:x,xlabel="Time [min]",widen=:false,grid=:y)

##############################################################
## Fig 6c: Plot with uncertainty
##############################################################

Θ = sample(C[1:100:end],1000)
x = [range(0.0,0.8,length=200); range(0.8,1.0,length=200)]

# Get m,l,u at various time points
t₁ = [1.0;5.0:5.0:30.0]
t₂ = 60.0:30.0:180.0
t  = [t₁;t₂]
L  = Array{Any,1}(undef,length(t))
U  = similar(L)
M  = similar(L)
for (i,tᵢ) in enumerate(t)
    # Calculate lower and upper ci
    l,u = pdf_ci(θ -> Wpdf(θ,tᵢ),x,Θ)
    m = Wpdf(θ,tᵢ)(x)
    # Normalise by mode
    L[i] = l / maximum(m)
    U[i] = u / maximum(m)
    M[i] = m / maximum(m)
end

# Deterministic model
T₁ = range(0.0,35.0,length=200)
T₂ = range(0.0,210.0,length=200)
X₁ = [solve_ode_model(t,collect(θdet))[2] / solve_ode_model(t,collect(θdet))[1] for t in T₁]
X₂ = [solve_ode_model(t,collect(θdet))[2] / solve_ode_model(t,collect(θdet))[1] for t in T₂]

fig6c1 = plot(X₁,T₁,lw=2.0,c=:red,ls=:dash)
for (i,tᵢ) in Iterators.reverse(enumerate(t₁))
    plot!(fig6c1,x,tᵢ .+ 3.33 * U[i],frange=tᵢ .+ 3.33 * L[i],c=RGBA(0.5,0.5,0.5,0.4),lw=0.0)
    plot!(fig6c1,x,tᵢ .+ 3.33 * M[i],c=:black,lw=2.0)
end
plot!(fig6c1,
    yticks = [1; 5:5:30],
    ylim = (0.0,35.0)
)

fig6c2 = plot(X₂,T₂,lw=2.0,c=:red,ls=:dash)
plot!(fig6c2,x,30.0 .+ 15.0 * U[length(t₁)],frange=30.0 .+ 15.0 * L[length(t₁)],c=RGBA(0.5,0.5,0.5,0.4),lw=0.0)
plot!(fig6c2,x,30.0 .+ 15.0 * M[length(t₁)],c=:black,lw=2.0)
for (i,tᵢ) in Iterators.reverse(enumerate(t₂))
    plot!(fig6c2,x,tᵢ .+ 15.0 * U[i+length(t₁)],frange=tᵢ .+ 15.0 * L[i+length(t₁)],c=RGBA(0.5,0.5,0.5,0.4),lw=0.0)
    plot!(fig6c2,x,tᵢ .+ 15.0 * M[i+length(t₁)],c=:black,lw=2.0)
end
plot!(fig6c2,
    yticks = 30:30:180,
    xlim=(0.9,1.0),
    ylim = (0.0,210.0)
)

fig6c = plot(fig6c1,fig6c2,
    xlabel = "Fraction internal antibody",
    ylabel = "Time [min]",
    xrotation = 90,
    yrotation = 90,
    xguidefontrotation = 180,
    ymirror = true,
    axis = :x,
    grid = :x,
    xflip = true,
    size=(400,600),
    bottom_margin=5Plots.mm,
    legend=:none
)

##############################################################
## Fig 6de: Proportion on surface
##############################################################

# Univariate (with uncertainty)
fun = p -> p[3] / sum(p[2:3])
fig6d = plot_distribution(θ -> param_dist(θ[5:end-1]),fun,C,θ,xlabel="Proportion on surface")
vline!(fig6d,[θdet.β / (θdet.λ + θdet.β)],lw=2.0,c=:red,ls=:dash)
plot!(fig6d,
    axis=:x,
    yticks=[],
    grid=:x,
)

# Bivariate (at best fit)
p = rand(param_dist(θ[5:end-1]),100_000)
S = fun.(eachcol(p))
R = p[1,:]
fig6e = density2d(S,R,
            xlabel = "Proportion on surface",
            ylabel = "R",
            fill   = true,
            lw     = 0.0,
            c      = :RdPu_9,
            levels = 8,
            box    = :on,
            ylim = (0.0,3.0),
            xlim = (0.0,1.0)
        )
scatter!(fig6e,[θdet.β / (θdet.λ + θdet.β)],[1.0],c=:white,widen=false,msw=0.0)

# Fig 6de
fig6de = plot(fig6d,fig6e,layout=grid(2,1),size=(300,500),legend=:none,widen=false)
add_plot_labels!(fig6de,offset=3)

##############################################################
## Save...
##############################################################
savefig(fig6ab,"$(@__DIR__)/fig6ab.svg")
savefig(fig6c,"$(@__DIR__)/fig6c.svg")
savefig(fig6de,"$(@__DIR__)/fig6de.svg")

ext: jl

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
using Documenter #, DocumenterLaTeX
using HMMGradients

makedocs(
    sitename = "HMMGradients",
    format = [Documenter.HTML()],#, DocumenterLaTeX.LaTeX()],
    authors = "Niccolò Antonello",
    modules = [HMMGradients],
    pages = [
        "Home" => "index.md",
        "Theory and Notation" => "1_intro.md",
        "Forward-Backward functions" => "2_fb.md",
        "Computing gradients" => "3_grads.md",
        "Demo" => "4_demo.md",
    ],
    doctest=false,
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/idiap/HMMGradients.jl.git",
    devbranch = "main",
)


text:
function initialPopulation(visitors, rooms, hoursNumber)
  Λ = []
  λ = 0
  k = 10
  while λ <= 1
    push!(Λ,λ)
    λ = λ + 1/k
    k = k + 1
  end # λ <= 1
  push!(Λ,1)
  Population = peek(visitors, rooms, hoursNumber, Λ, 1, [seed])
  for i in 1:length(Λ)
    push!(Population, grasp(visitors, rooms, hoursNumber, 1, 0.001, seed)
  end # i in 1:20
  return Population
end # initialPopulation function

function rank!(Population)
  if length(Population) > 1
    middle = Int(length(Population) / 2)
    left = Population[1:middle]
    right = Population[middle+1:length(Population)]
    rank!(left)
    rank!(right)

    #fusion
    aux = []
    i = 1
    j = 1
    while i <= length(left) && j <= length(right)
      if left[i][1] > right[j][1]
        push!(aux, left[i])
        i += 1
      elseif left[i][1] < right[j][1]
        push!(aux, right[j])
        j += 1
      else
        if left[i][2] <= right[j][2]
          push!(aux, left[i])
          i += 1
        else
          push!(aux, right[j])
          j += 1
        end # if 2
      end # if 1
    end # while fusion
    if i = length(left)
      for l in j:length(right)
        push!(aux,right[l])
      end
    else
      for l in i:length(left)
        push!(left[l])
      end
    end
    Population = aux
  end # if > 1
end # rank function

function genetic(visitors, rooms, hoursNumber)
Population = initialPopulation(visitors, rooms, hoursNumber)

end # genetic function


text:
using DashBootstrapComponents

inline_radioitems = dbc_formgroup([
    dbc_label("Choose one"),
    dbc_radioitems(
        options=[
            Dict("label" => "Option 1", "value" => 1),
            Dict("label" => "Option 2", "value" => 2),
        ],
        value=1,
        id="radioitems-inline-input",
        inline=true,
    ),
]);

inline_checklist = dbc_formgroup([
    dbc_label("Choose a bunch"),
    dbc_checklist(
        options=[
            Dict("label" => "Option 1", "value" => 1),
            Dict("label" => "Option 2", "value" => 2),
        ],
        value=[],
        id="checklist-inline-input",
        inline=true,
    ),
]);

inline_switches = dbc_formgroup([
    dbc_label("Toggle a bunch"),
    dbc_checklist(
        options=[
            Dict("label" => "Option 1", "value" => 1),
            Dict("label" => "Option 2", "value" => 2),
        ],
        value=[],
        id="switches-inline-input",
        inline=true,
        switch=true,
    ),
]);

inline_inputs = dbc_form([inline_radioitems, inline_checklist, inline_switches]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module BoardGamesGames

using BoardGames

module TicTacToeGame
    include("TicTacToe.jl")

    export
        #types
        TicTacToe,
        TicTacToeBoard
end

module OthelloGame
    include("Othello.jl")

    export
        #types
        Othello,
        OthelloBoard
end
end


text:
#
# ITensors.jl extensions
#

# Generalize siteind to n-dimensional lattice
function ITensors.siteind(st::SiteType, N1::Integer, N2::Integer, Ns::Integer...; kwargs...)
  s = siteind(st; kwargs...)
  if !isnothing(s)
    ts = "n1=$N1,n2=$N2"
    for i in eachindex(Ns)
      ts *= ",n$(i + 2)=$(Ns[i])"
    end
    return addtags(s, ts)
  end
  return isnothing(s) && error(space_error_message(st))
end

# Generalize siteinds to n-dimensional lattice
function ITensors.siteinds(
  str::AbstractString, N1::Integer, N2::Integer, Ns::Integer...; kwargs...
)
  st = SiteType(str)
  return [siteind(st, ns...) for ns in Base.product(1:N1, 1:N2, UnitRange.(1, Ns)...)]
end

# Get the promoted type of the Index objects in a collection
# of Index (Tuple, Vector, ITensor, etc.)
indtype(i::Index) = typeof(i)
indtype(T::Type{<:Index}) = T
indtype(is::Tuple{Vararg{<:Index}}) = eltype(is)
indtype(is::Vector{<:Index}) = eltype(is)
indtype(A::ITensor...) = indtype(inds.(A))

indtype(tn1, tn2) = promote_type(indtype(tn1), indtype(tn2))
indtype(tn) = mapreduce(indtype, promote_type, tn)

#
# MPS functionality extensions
#

Base.keytype(m::MPS) = keytype(data(m))

# A version of indexing which returns an empty order-0 ITensor
# when out of bounds
get_itensor(x::MPS, n::Int) = n in 1:length(x) ? x[n] : ITensor()

# Reverse the site ordering of an MPS.
# XXX: also reverse the orthogonality limits.
Base.reverse(x::MPS) = MPS(reverse(x.data))


text:
# Julia wrapper for header: /usr/local/ffmpeg/2.4/include/libavutil/frame.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_frame_get_best_effort_timestamp,
    av_frame_set_best_effort_timestamp,
    av_frame_get_pkt_duration,
    av_frame_set_pkt_duration,
    av_frame_get_pkt_pos,
    av_frame_set_pkt_pos,
    av_frame_get_channel_layout,
    av_frame_set_channel_layout,
    av_frame_get_channels,
    av_frame_set_channels,
    av_frame_get_sample_rate,
    av_frame_set_sample_rate,
    av_frame_get_metadata,
    av_frame_set_metadata,
    av_frame_get_decode_error_flags,
    av_frame_set_decode_error_flags,
    av_frame_get_pkt_size,
    av_frame_set_pkt_size,
    avpriv_frame_get_metadatap,
    av_frame_get_qp_table,
    av_frame_set_qp_table,
    av_frame_get_colorspace,
    av_frame_set_colorspace,
    av_frame_get_color_range,
    av_frame_set_color_range,
    av_get_colorspace_name,
    av_frame_alloc,
    av_frame_free,
    av_frame_ref,
    av_frame_clone,
    av_frame_unref,
    av_frame_move_ref,
    av_frame_get_buffer,
    av_frame_is_writable,
    av_frame_make_writable,
    av_frame_copy,
    av_frame_copy_props,
    av_frame_get_plane_buffer,
    av_frame_new_side_data,
    av_frame_get_side_data,
    av_frame_remove_side_data,
    av_frame_side_data_name

function av_frame_get_best_effort_timestamp(frame)
    ccall((:av_frame_get_best_effort_timestamp,libavutil),Int64,(Ptr{AVFrame},),frame)
end

function av_frame_set_best_effort_timestamp(frame,val::Int64)
    ccall((:av_frame_set_best_effort_timestamp,libavutil),Void,(Ptr{AVFrame},Int64),frame,val)
end

function av_frame_get_pkt_duration(frame)
    ccall((:av_frame_get_pkt_duration,libavutil),Int64,(Ptr{AVFrame},),frame)
end

function av_frame_set_pkt_duration(frame,val::Int64)
    ccall((:av_frame_set_pkt_duration,libavutil),Void,(Ptr{AVFrame},Int64),frame,val)
end

function av_frame_get_pkt_pos(frame)
    ccall((:av_frame_get_pkt_pos,libavutil),Int64,(Ptr{AVFrame},),frame)
end

function av_frame_set_pkt_pos(frame,val::Int64)
    ccall((:av_frame_set_pkt_pos,libavutil),Void,(Ptr{AVFrame},Int64),frame,val)
end

function av_frame_get_channel_layout(frame)
    ccall((:av_frame_get_channel_layout,libavutil),Int64,(Ptr{AVFrame},),frame)
end

function av_frame_set_channel_layout(frame,val::Int64)
    ccall((:av_frame_set_channel_layout,libavutil),Void,(Ptr{AVFrame},Int64),frame,val)
end

function av_frame_get_channels(frame)
    ccall((:av_frame_get_channels,libavutil),Cint,(Ptr{AVFrame},),frame)
end

function av_frame_set_channels(frame,val::Integer)
    ccall((:av_frame_set_channels,libavutil),Void,(Ptr{AVFrame},Cint),frame,val)
end

function av_frame_get_sample_rate(frame)
    ccall((:av_frame_get_sample_rate,libavutil),Cint,(Ptr{AVFrame},),frame)
end

function av_frame_set_sample_rate(frame,val::Integer)
    ccall((:av_frame_set_sample_rate,libavutil),Void,(Ptr{AVFrame},Cint),frame,val)
end

function av_frame_get_metadata(frame)
    ccall((:av_frame_get_metadata,libavutil),Ptr{AVDictionary},(Ptr{AVFrame},),frame)
end

function av_frame_set_metadata(frame,val)
    ccall((:av_frame_set_metadata,libavutil),Void,(Ptr{AVFrame},Ptr{AVDictionary}),frame,val)
end

function av_frame_get_decode_error_flags(frame)
    ccall((:av_frame_get_decode_error_flags,libavutil),Cint,(Ptr{AVFrame},),frame)
end

function av_frame_set_decode_error_flags(frame,val::Integer)
    ccall((:av_frame_set_decode_error_flags,libavutil),Void,(Ptr{AVFrame},Cint),frame,val)
end

function av_frame_get_pkt_size(frame)
    ccall((:av_frame_get_pkt_size,libavutil),Cint,(Ptr{AVFrame},),frame)
end

function av_frame_set_pkt_size(frame,val::Integer)
    ccall((:av_frame_set_pkt_size,libavutil),Void,(Ptr{AVFrame},Cint),frame,val)
end

function avpriv_frame_get_metadatap(frame)
    ccall((:avpriv_frame_get_metadatap,libavutil),Ptr{Ptr{AVDictionary}},(Ptr{AVFrame},),frame)
end

function av_frame_get_qp_table(f,stride,_type)
    ccall((:av_frame_get_qp_table,libavutil),Ptr{Int8},(Ptr{AVFrame},Ptr{Cint},Ptr{Cint}),f,stride,_type)
end

function av_frame_set_qp_table(f,buf,stride::Integer,_type::Integer)
    ccall((:av_frame_set_qp_table,libavutil),Cint,(Ptr{AVFrame},Ptr{AVBufferRef},Cint,Cint),f,buf,stride,_type)
end

function av_frame_get_colorspace(frame)
    ccall((:av_frame_get_colorspace,libavutil),Cint,(Ptr{AVFrame},),frame)
end

function av_frame_set_colorspace(frame,val::AVColorSpace)
    ccall((:av_frame_set_colorspace,libavutil),Void,(Ptr{AVFrame},AVColorSpace),frame,val)
end

function av_frame_get_color_range(frame)
    ccall((:av_frame_get_color_range,libavutil),Cint,(Ptr{AVFrame},),frame)
end

function av_frame_set_color_range(frame,val::AVColorRange)
    ccall((:av_frame_set_color_range,libavutil),Void,(Ptr{AVFrame},AVColorRange),frame,val)
end

function av_get_colorspace_name(val::AVColorSpace)
    ccall((:av_get_colorspace_name,libavutil),Ptr{Uint8},(AVColorSpace,),val)
end

function av_frame_alloc()
    ccall((:av_frame_alloc,libavutil),Ptr{AVFrame},())
end

function av_frame_free(frame)
    ccall((:av_frame_free,libavutil),Void,(Ptr{Ptr{AVFrame}},),frame)
end

function av_frame_ref(dst,src)
    ccall((:av_frame_ref,libavutil),Cint,(Ptr{AVFrame},Ptr{AVFrame}),dst,src)
end

function av_frame_clone(src)
    ccall((:av_frame_clone,libavutil),Ptr{AVFrame},(Ptr{AVFrame},),src)
end

function av_frame_unref(frame)
    ccall((:av_frame_unref,libavutil),Void,(Ptr{AVFrame},),frame)
end

function av_frame_move_ref(dst,src)
    ccall((:av_frame_move_ref,libavutil),Void,(Ptr{AVFrame},Ptr{AVFrame}),dst,src)
end

function av_frame_get_buffer(frame,align::Integer)
    ccall((:av_frame_get_buffer,libavutil),Cint,(Ptr{AVFrame},Cint),frame,align)
end

function av_frame_is_writable(frame)
    ccall((:av_frame_is_writable,libavutil),Cint,(Ptr{AVFrame},),frame)
end

function av_frame_make_writable(frame)
    ccall((:av_frame_make_writable,libavutil),Cint,(Ptr{AVFrame},),frame)
end

function av_frame_copy(dst,src)
    ccall((:av_frame_copy,libavutil),Cint,(Ptr{AVFrame},Ptr{AVFrame}),dst,src)
end

function av_frame_copy_props(dst,src)
    ccall((:av_frame_copy_props,libavutil),Cint,(Ptr{AVFrame},Ptr{AVFrame}),dst,src)
end

function av_frame_get_plane_buffer(frame,plane::Integer)
    ccall((:av_frame_get_plane_buffer,libavutil),Ptr{AVBufferRef},(Ptr{AVFrame},Cint),frame,plane)
end

function av_frame_new_side_data(frame,_type::AVFrameSideDataType,size::Integer)
    ccall((:av_frame_new_side_data,libavutil),Ptr{AVFrameSideData},(Ptr{AVFrame},AVFrameSideDataType,Cint),frame,_type,size)
end

function av_frame_get_side_data(frame,_type::AVFrameSideDataType)
    ccall((:av_frame_get_side_data,libavutil),Ptr{AVFrameSideData},(Ptr{AVFrame},AVFrameSideDataType),frame,_type)
end

function av_frame_remove_side_data(frame,_type::AVFrameSideDataType)
    ccall((:av_frame_remove_side_data,libavutil),Void,(Ptr{AVFrame},AVFrameSideDataType),frame,_type)
end

function av_frame_side_data_name(_type::AVFrameSideDataType)
    ccall((:av_frame_side_data_name,libavutil),Ptr{Uint8},(AVFrameSideDataType,),_type)
end


text:
using Optim
using LinearAlgebra: I, det, tr

# Rosenbrock
f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
x0 = [0.0, 0.0]
# optimize(f, x0, LBFGS(), autodiff = :forward)
# optimize(f, x0, LBFGS(), autodiff = :forward)

# analytical gradient
function g!(G, x)
    G[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
    G[2] = 200.0 * (x[2] - x[1]^2)
end
res = optimize(f, g!, x0)

@show summary(res)
@show Optim.minimizer(res)
@show Optim.minimum(res)
@show Optim.iterations(res)
@show Optim.iteration_limit_reached(res)
@show Optim.f_calls(res)
@show Optim.converged(res)

#
# # analytical Hessian
# function h!(H, x)
#     H[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
#     H[1, 2] = -400.0 * x[1]
#     H[2, 1] = -400.0 * x[1]
#     H[2, 2] = 200.0
# end
# optimize(f, g!, h!, x0, LBFGS())
#
lower = [-0.95, -2.1]
upper = [Inf, Inf]
initial_x = [0.0, 0.0]
# # requires using LineSearches
inner_optimizer = GradientDescent()
# res = optimize(f, g!, lower, upper, initial_x, Fminbox(inner_optimizer))

res = optimize(f, lower, upper, initial_x)

@show summary(res)
@show Optim.minimizer(res)
@show Optim.minimum(res)
@show Optim.iterations(res)
@show Optim.iteration_limit_reached(res)
@show Optim.f_calls(res)
@show Optim.converged(res)

# function fMat(Xm::Matrix{Float64})
#     # Xm = reshape(X,(3,3))
#     return det((Xm-I) * (Xm-I))
# end
# X0 = reshape(collect(1:1.0:9.0), (3,3))
# optimize(fMat, X0, LBFGS())
#
# function sq(x)
#     return x^2
# end
#
# function get_fn(a, b)::Function
#     function calc_fn(x::Vector{Float64})
#         return sum((x .* sq(a) .+ b).^2)
#     end
#     return calc_fn
# end
#
# f = get_fn(2.0,3.0)
# res = optimize(f, [1.0, 2.0], LBFGS())


text:
using Distributed, StanSample, CSV, DataFrames, StatsPlots, Statistics

ProjDir = @__DIR__
#include(joinpath(ProjDir, "redcardstudy.jl"))

function timings(model, nts, ncs; N=6, use_cpp=[false, true])

    df = DataFrame()
    res_t = zeros(N)
    println("\n\n")
    for nt in nts                   # Number of C++ threads
      for n in ncs                # Number of chains
        println("(num_threads=$nt, num_cpp_chains=1, \
          num_julia_chains=$n) runs")
        for i in 1:N
            res_t[i] = @elapsed stan_sample(model; data,
              use_cpp_chains=true, check_num_chains=false,
              num_threads=nt,
              num_cpp_chains=1, num_julia_chains=n)
            println("Iteration=$i, time=$(res_t[i])")
        end
        append!(df, DataFrame(
            num_threads=nt,
            num_cpp_chains=1,
            num_julia_chains=n,
            min=minimum(res_t),
            median=median(res_t),
            max=maximum(res_t))
        )
        if n > 1
          println("(num_threads=$nt, num_cpp_chains=$n, \
            num_julia_chains=1) runs")
          for i in 1:N
              res_t[i] = @elapsed stan_sample(model; data,
                use_cpp_chains=true, check_num_chains=false,
                num_threads=nt,
                num_cpp_chains=n, num_julia_chains=1)
              println("Iteration=$i, time=$(res_t[i])")
          end
          append!(df, DataFrame(
              num_threads=nt,
              num_cpp_chains=n,
              num_julia_chains=1,
              min=minimum(res_t),
              median=median(res_t),
              max=maximum(res_t))
          )
        end
      end
  end

    df
end

nts = [1, 4, 8]
ncs = [1, 2, 4]
N = 2

println("\nlogistic_0 runs\n")
model = logistic_0
arm_log_0_df = timings(model, nts, ncs; N)
arm_log_0_df |> display

CSV.write(joinpath(ProjDir, "results", "arm_log_0_df.csv"), arm_log_0_df)

println("\nlogistic_1 runs\n")
model = logistic_1
arm_log_1_df = timings(model, nts, ncs; N)
arm_log_1_df |> display

CSV.write(joinpath(ProjDir, "results", "arm_log_1_df.csv"), arm_log_1_df)


text:
################################################################################
#
#  Types and memory management for arb
#
################################################################################

# typealias arb_t arb
const arb_ptr = Vector{arb_t}
const arb_srcptr = arb_ptr

###############################################################################
#
#   Basic manipulation
#
###############################################################################

zero(x::arb) = parent(x)(0)
one(x::arb) = parent(x)(1)

arb_rel_accuracy_bits(x::arb_t) = accuracy_bits(x)

arb_bits(x::arb_t) = ccall((:arb_bits, :libarb), Fslong, (Ref{arb_t}, ), x)

################################################################################
#
#  String I/O
#
################################################################################

arb_printn(x::arb_t, digits, flags = 0) = ccall((:arb_printn, :libarb), Void, (Ref{arb_t}, Fslong, Fulong), x, digits, flags)

################################################################################
#
#  Predicates
#
################################################################################

arb_is_zero(x::arb_t) = iszero(x)
arb_is_nonzero(x::arb_t) = isnonzero(x)
arb_is_one(x::arb_t) = isone(x)

arb_is_finite(x::arb_t) = isfinite(x)
arb_is_exact(x::arb_t) = isexact(x)
arb_is_int(x::arb_t) = isint(x)

arb_is_positive(x::arb_t) = ispositive(x)
arb_is_nonnegative(x::arb_t) = isnonnegative(x)
arb_is_negative(x::arb_t) = isnegative(x)
arb_is_nonpositive(x::arb_t) = isnonpositive(x)

arb_contains_zero(x::arb_t) = Bool(ccall((:arb_contains_zero, :libarb), Cint, (Ref{arb_t}, ), x))
arb_contains_negative(x::arb_t) = Bool(ccall((:arb_contains_negative, :libarb), Cint, (Ref{arb_t}, ), x))
arb_contains_nonpositive(x::arb_t) = Bool(ccall((:arb_contains_nonpositive, :libarb), Cint, (Ref{arb_t}, ), x))
arb_contains_positive(x::arb_t) = Bool(ccall((:arb_contains_positive, :libarb), Cint, (Ref{arb_t}, ), x))
arb_contains_nonnegative(x::arb_t) = Bool(ccall((:arb_contains_nonnegative, :libarb), Cint, (Ref{arb_t}, ), x))

################################################################################
#
#  Unsafe setting
#
################################################################################

for (typeofy, passtoc) in ((arb_t, Ref{arb_t}), (Ptr{arb_t}, Ptr{arb_t}))
    @eval begin
        arb_set(y::($typeofy), x::arb_t) = _arb_set(y, x)
        arb_set_round(y::($typeofy), x::arb_t, prec) = _arb_set(y, x, prec)
        arb_set_si(y::($typeofy), x::Fslong) = _arb_set(y, x)
        arb_set_ui(y::($typeofy), x::Fulong) = _arb_set(y, x)
        arb_set_d(y::($typeofy), x::Cdouble) = _arb_set(y, x)
        arb_set_fmpz(y::($typeofy), x::fmpz) = _arb_set(y, x)
        arb_set_round_fmpz(y::($typeofy), x::fmpz, prec) = _arb_set(y, x, prec)
        arb_set_fmpq(y::($typeofy), x::fmpq, prec) = _arb_set(y, x, prec)
        arb_set_str(y::($typeofy), x::AbstractString, prec) = _arb_set(y, x, prec)

        arb_set_arf(y::($typeofy), x::arf_t) =  ccall((:arb_set_arf, :libarb), Void, (($passtoc), Ref{arf_t}), y, x)
        arb_set_interval_arf(x::($typeofy), a::arf_t, b::arf_t, prec) = ccall((:arb_set_interval_arf, :libarb), Void, (($passtoc), Ref{arf_t}, Ref{arf_t}, Fslong), x, a, b, prec)
    end
end


text:
polynorm_l2_n_d_ms = [
    [
    (1, 1, 2), # compile run
    (1, 20, 4),
    (1, 20, 8),
    (1, 20, 16),
    (1, 20, 32),
    (1, 20, 64),
    ],
    [
    (1, 1, 2), # compile run
    (1, 40, 4),
    (1, 40, 8),
    (1, 40, 16),
    (1, 40, 32),
    (1, 40, 64),
    ],
    [
    (3, 1, 2), # compile run
    (4, 2, 4),
    (4, 2, 8),
    (4, 2, 16),
    (4, 2, 32),
    ],
    [
    (3, 1, 2), # compile run
    (4, 4, 4),
    (4, 4, 8),
    (4, 4, 16),
    (4, 4, 32),
    ],
    ]

polynorm_l1_n_d_ms = [
    [
    (1, 1, 2), # compile run
    (1, 40, 8),
    (1, 40, 16),
    (1, 40, 32),
    (1, 40, 64),
    (1, 40, 128),
    ],
    [
    (1, 1, 2), # compile run
    (1, 80, 8),
    (1, 80, 16),
    (1, 80, 32),
    (1, 80, 64),
    (1, 80, 128),
    ],
    [
    (3, 1, 2), # compile run
    (4, 2, 8),
    (4, 2, 16),
    (4, 2, 32),
    (4, 2, 64),
    (4, 2, 128),
    ],
    [
    (3, 1, 2), # compile run
    (4, 4, 4),
    (4, 4, 8),
    (4, 4, 16),
    (4, 4, 32),
    (4, 4, 64),
    ],
    ]

polynorm_insts(
    use_l1::Bool,
    use_norm_cone::Bool,
    use_wsos::Bool,
    d_factors::Vector{Int},
    ) = [
    [(n, d, f * d, m, use_l1, use_norm_cone, use_wsos) for (n, d, m) in ndms]
    for f in d_factors
    for ndms in (use_l1 ? polynorm_l1_n_d_ms : polynorm_l2_n_d_ms)
    ]

insts = OrderedDict()
insts["nat"] = (nothing, vcat(
    polynorm_insts(false, true, false, [1, 2]),
    polynorm_insts(true, true, false, [1,]),
    ))
insts["ext"] = (nothing, vcat(
    polynorm_insts(false, false, true, [1, 2]),
    polynorm_insts(true, false, true, [1,]),
    ))
insts["extmat"] = (nothing, vcat(
    polynorm_insts(false, false, false, [1, 2]),
    ))

return (PolyNormJuMP, insts)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

using Test, ArrowTypes, UUIDs

include("tests.jl")

max_stars_repo_name: jonalm/arrow-julia

max_issues_repo_name: jonalm/arrow-julia

max_forks_repo_name: jonalm/arrow-julia


text:
# ```julia
# using Coad
# pkg"add Coad, Plots"
# ```
push!(LOAD_PATH, pwd())

using Coad
using Plots
using Printf

# Set up grid / model simulation
model = Coad1D(
    n = 6,
    kernel = :long,
)

# Set up initial conditions
r̅=(10.0)*1e-6
dist = ExponentialDist(L̅=(1.0)*1e-3, x̅=mass_from_r(r̅))
set!(model, dist)

println("\n Plotting initial conditions...")
p = plot(
    model.rᵢ*1e6, model.gᵢ*1e3, label="t = 0 min",
    xaxis=:log, xlabel="r (μm)",
    # xlim=(0.5, 5000), xticks=[1, 10, 100, 1000],
    # ylim=(0, 0.9), yticks=0:0.06:0.9,
    ylabel="g (g / m³)",
)
display(p)

t = Timestamp(0)
tmax = 60*60 + 1 # seconds
Δt = 5.0 # s
Δt_diag = 5 * 60
Δt_plot = 10 * 60
nt = ceil(Integer, tmax / Δt)

is_diag_step(t; Δt_diag=Δt_diag) = t % Δt_diag == Timestamp(0)
is_plot_step(t; Δt_plot=Δt_plot) = t % Δt_plot == Timestamp(0)

###
# using Profile
# using ProfileSVG
# Profile.clear()
# @profile step!(model, Δt)
# ProfileSVG.save(joinpath("assets", "prof.svg"))
###

println("\n Beginning simulation...")
println(" -----------------------")

g_diags = []
for i in 1:nt
    global t += Δt
    step!(model, Δt)

    mass_tot = sum(model.gᵢ)
    @printf "t = %s | mass %10.3e\n" t mass_tot

    if is_diag_step(t)
        push!(g_diags, model.gᵢ[:])
    end

    if is_plot_step(t)
        display(plot!(p, model.rᵢ*1e6, model.gᵢ*1e3, label = "t = $(t.minutes) min"))
    end
end

println("End; press any key to close.")
xxx = readline()

# using NPZ
# for (i, g) in enumerate(g_plots)
#     file_fn = @sprintf "output_%03d.npy" i
#     @printf "%03d %s %10.3e \n" i file_fn maximum(g)
#     file_pth = joinpath("output", file_fn)
#     npzwrite(file_pth, g)
# end
# npzwrite("output/r_grid.npy", model.rᵢ)

# using NPZ
# for (i, g) in enumerate(g_diags)
#     file_fn = @sprintf "output_%03d.npy" i
#     @printf "%03d %s %10.3e \n" i file_fn maximum(g)
#     file_pth = joinpath("output", file_fn)
#     npzwrite(file_pth, g)
# end
# npzwrite("output/r_grid.npy", model.rᵢ)

# Set up integrator
# solver = Integrator1D(model, Δt=5.0, tmax=3601)

# Run
# run!(solver)

# Do something with the outputs?


text:
using Random, Distributions

Random.seed!(1234)

description = "Reduced four-state model (C_A, O_R, O_A)"

# indices of states below are relative to the reduced, 3-state model.
# open index
oi = 2
# activated states
ai = [1,3]
# indices to save when solving
states_to_save = [2, 3]
# parameter names
parnames = ["a0", "sa", "b0", "sb", "c0", "sc", "d0", "sd", "g1", "g2", "h1", "h2"]

"""
Notes on how performance can be improved:
* Modify the objective function so that `db_pars` is called only once per evoluation of the objective function.
    i.e. `tm(...)` should not have to call `db_pars`
* Modify `tm(...)` to update a pre-allocated array in-place
* Use a static vector `SVector` for parameters, e.g. the output of `db_pars`
* Use views of the transition matrix to define the ODE, i.e. argument `Q` in `rhs!(...)`

Similar performance tips are available in `main_v4.jl`.
The above apply equally to `six_hum.jl`.
"""

function db_pars(pars::AbstractVector{}, db; fmt="x")::Vector{}
    #pars are ln-transformed
    if db == "x"
        a0, sa, b0, sb, c0, sc, d0, sd, g1, g2, h1, h2 = exp.(pars)
    else
        if db in ["g1", "h1", "g2", "h2"]
            a0, sa, b0, sb, c0, sc, d0, sd = exp.(pars[1:8])

            if db == "g1"
                g2, h1, h2 = exp.(pars[9:end])
                g1 = (h1*b0*g2*c0)/(d0*h2*a0)
            elseif db == "h1"
                g1, g2, h2 = exp.(pars[9:end])
                h1 = (g1*d0*h2*a0)/(b0*g2*c0)
            elseif db == "g2"
                g1, h1, h2 = exp.(pars[9:end])
                g2 = (g1*d0*h2*a0)/(h1*b0*c0)
            elseif db == "h2"
                g1, g2, h1 = exp.(pars[9:end])
                h2 = (g2*c0*h1*b0)/(a0*g1*d0)
            end

        else
            if db == "c"
                a0, sa, b0, sb, d0, sd, g1, g2, h1, h2 = exp.(pars)
                c0 = (a0*d0*g1*h2)/(b0*g2*h1)
                sc = 1/( (1/sb) - (1/sd) + (1/sa) )

            elseif db == "b"
                a0, sa, c0, sc, d0, sd, g1, g2, h1, h2 = exp.(pars)
                b0 = (a0*g1*d0*h2)/(g2*c0*h1)
                sb = 1/( (1/sc) + (1/sd) - (1/sa) )
            elseif db == "d"
                a0, sa, b0, sb, c0, sc, g1, g2, h1, h2 = exp.(pars)
                d0 = (b0*g2*c0*h1)/(a0*g1*h2)
                sd = 1/( (1/sb) + (1/sa) - (1/sc) )
            elseif db == "a"
                b0, sb, c0, sc, d0, sd, g1, g2, h1, h2 = exp.(pars)
                a0 = (b0*g2*c0*h1)/(g1*d0*h2)
                sa = 1/( (1/sc) - (1/sb) + (1/sd) )
            else
                error("   $db not understood as detailed balance condition.")
            end
        end
    end

    if fmt == "log"
        return log.([a0, sa, b0, sb, c0, sc, d0, sd, g1, g2, h1, h2])
    else
        return [a0, sa, b0, sb, c0, sc, d0, sd, g1, g2, h1, h2]
    end
end

function tm(v, p::AbstractVector{}, db)
    a0, sa, b0, sb, c0, sc, d0, sd, g1, g2, h1, h2 = db_pars(p, db)

    # hyperpolarization-activated
    a1 = a0*exp(-v/sa)
    b1 = b0*exp(v/sb)
    a2 = c0*exp(-v/sc)
    b2 = d0*exp(v/sd)

    Q = [a1  -(b1 + g1 + a1)  -a1  h1-a1;
        g2  -g2  -(a2+h2+g2)  b2-g2;
        0  g1  a2  -(b2 + h1)]

    return Q
end

function get_vpars(v, p::AbstractVector{}, db)
    a0, sa, b0, sb, c0, sc, d0, sd, g1, g2, h1, h2 = db_pars(p, db)

    a1 = a0*exp(-v/sa)
    b1 = b0*exp(v/sb)
    a2 = c0*exp(-v/sc)
    b2 = d0*exp(v/sd)

    return [a1, b1, a2, b2, g1, g2, h1, h2, v]
end

function rhs!(du, u, Q, t)
    du .= Q[:,2:end] * u + Q[:,1]
end

function prepare_BoundsAndParams(pars::Array{Float64,1}, bounds, db)
    # remove entries from parameter and bound arrays given detailed balance
    if db == true || db == "x"
        return pars, bounds
    else
        ps = copy(pars)
        bds = copy(bounds)
        if db == "a"
            deleteat!(ps, 1:2)
            deleteat!(bds, 1:2)
        elseif db == "b"
            deleteat!(ps, 3:4)
            deleteat!(bds, 3:4)
        elseif db == "c"
            deleteat!(ps, 5:6)
            deleteat!(bds, 5:6)
        elseif db == "d"
            deleteat!(ps, 7:8)
            deleteat!(bds, 7:8)
        else
            println("$db not understood. Available are: 'a', 'b', 'c', 'd'.")
        end
    end
    return ps, bds
end

bds = [
        (1e-9, 1), #b0
        (7, 200), #sb

        (1e-9, 1e3), #a0
        (10, 4000), #sa

        (1e-12, 1e-1), #d0
        (50, 200), #sd

        (1e-12, 1e-1), #c0
        (10, 200), #sc,

        (1e-9, 10.), #g1
        (1e-8, 10), #g2

        (1e-9, 10), #h1
        (1e-9, 1e3), #h2
]

hnc2_bds = [
        (1e-7, 1), #a0
        (8, 300), #sa

        (20, 1000), #b0
        (25, 500), #sb

        (1e-8, 1), #c0
        (10, 100), #sc,

        (1e-8, 1), #d0
        (30, 100), #sd

        (1e-6, 1.), #g1
        (1e-12, 1e-3), #g2

        (1e-9, 1e-1), #h1
        (1e-9, 1e-1), #h2

]

ln_par_bds = [log.(x) for x in hnc2_bds]

max_stars_repo_name: haganenoneko/FunnyChannelGatingReview

max_issues_repo_name: haganenoneko/FunnyChannelGatingReview

max_forks_repo_name: haganenoneko/FunnyChannelGatingReview


text:
"""
# Module KernelAPI

The module `KernelAPI` is part of the package `Perf` and implements the
interface to the Linux kernel's performance event API.
"""
module KernelAPI

#
#  Types & Constants
#

"""

Directly taken from Linux definition
"""
@enum( PerfTypeId::UInt32,
    PERF_TYPE_HARDWARE   = 0,
    PERF_TYPE_SOFTWARE   = 1,
    PERF_TYPE_TRACEPOINT = 2,
    PERF_TYPE_HW_CACHE   = 3,
    PERF_TYPE_RAW        = 4,
    PERF_TYPE_BREAKPOINT = 5,
    PERF_TYPE_MAX,       # non-ABI
)

"""
Generalized performance event event_id types, used by the
attr.event_id parameter of the sys_perf_event_open()
syscall. Common hardware events, generalized by the kernel.

Direct copy of kernel definition.
"""
@enum( PerfHwId::UInt64,
    PERF_COUNT_HW_CPU_CYCLES                = 0,
    PERF_COUNT_HW_INSTRUCTIONS              = 1,
    PERF_COUNT_HW_CACHE_REFERENCES          = 2,
    PERF_COUNT_HW_CACHE_MISSES              = 3,
    PERF_COUNT_HW_BRANCH_INSTRUCTIONS       = 4,
    PERF_COUNT_HW_BRANCH_MISSES             = 5,
    PERF_COUNT_HW_BUS_CYCLES                = 6,
    PERF_COUNT_HW_STALLED_CYCLES_FRONTEND   = 7,
    PERF_COUNT_HW_STALLED_CYCLES_BACKEND    = 8,
    PERF_COUNT_HW_REF_CPU_CYCLES            = 9,
    PERF_COUNT_HW_MAX,                      # non-ABI
)

"""
Bits that can be set in attr.sample_type to request information
in the overflow packets.

Direct copy of kernel definition.
"""
@enum( PerfEventSampleFormat::UInt32,
    PERF_SAMPLE_IP              = one(UInt32) << 0,
    PERF_SAMPLE_TID             = one(UInt32) << 1,
    PERF_SAMPLE_TIME            = one(UInt32) << 2,
    PERF_SAMPLE_ADDR            = one(UInt32) << 3,
    PERF_SAMPLE_READ            = one(UInt32) << 4,
    PERF_SAMPLE_CALLCHAIN       = one(UInt32) << 5,
    PERF_SAMPLE_ID              = one(UInt32) << 6,
    PERF_SAMPLE_CPU             = one(UInt32) << 7,
    PERF_SAMPLE_PERIOD          = one(UInt32) << 8,
    PERF_SAMPLE_STREAM_ID       = one(UInt32) << 9,
    PERF_SAMPLE_RAW             = one(UInt32) << 10,
    PERF_SAMPLE_BRANCH_STACK    = one(UInt32) << 11,
    PERF_SAMPLE_REGS_USER       = one(UInt32) << 12,
    PERF_SAMPLE_STACK_USER      = one(UInt32) << 13,
    PERF_SAMPLE_WEIGHT          = one(UInt32) << 14,
    PERF_SAMPLE_DATA_SRC        = one(UInt32) << 15,
    PERF_SAMPLE_IDENTIFIER      = one(UInt32) << 16,
    PERF_SAMPLE_TRANSACTION     = one(UInt32) << 17,
    PERF_SAMPLE_REGS_INTR       = one(UInt32) << 18,
    PERF_SAMPLE_MAX             = one(UInt32) << 19, #  non-ABI
)

"""
The format of the data returned by read() on a perf event fd,
as specified by attr.read_format:

struct read_format {
 { u64      value;
   { u64        time_enabled; } && PERF_FORMAT_TOTAL_TIME_ENABLED
   { u64        time_running; } && PERF_FORMAT_TOTAL_TIME_RUNNING
   { u64        id;           } && PERF_FORMAT_ID
 } && !PERF_FORMAT_GROUP

 { u64      nr;
   { u64        time_enabled; } && PERF_FORMAT_TOTAL_TIME_ENABLED
   { u64        time_running; } && PERF_FORMAT_TOTAL_TIME_RUNNING
   { u64        value;
     { u64  id;           } && PERF_FORMAT_ID
   }        cntr[nr];
 } && PERF_FORMAT_GROUP
};

Direct copy of kernel definition.
"""
@enum( PerfEventReadFormat::UInt64,
    PERF_FORMAT_UNKNOWN             = 0,
    PERF_FORMAT_TOTAL_TIME_ENABLED  = one(UInt64) << 0,
    PERF_FORMAT_TOTAL_TIME_RUNNING  = one(UInt64) << 1,
    PERF_FORMAT_ID                  = one(UInt64) << 2,
    PERF_FORMAT_GROUP               = one(UInt64) << 3,
    PERF_FORMAT_MAX                 = one(UInt64) << 4,     # non-ABI
)

"""
Keep in sync with kernel definition of `perf_event_attr` in
`<kernel_headers_dir>/include/uapi/linux/perf_event.h`.

It appears there are several versions of this struct, as defined by
```
#define PERF_ATTR_SIZE_VER0 64  /* sizeof first published struct */
#define PERF_ATTR_SIZE_VER1 72  /* add: config2 */
#define PERF_ATTR_SIZE_VER2 80  /* add: branch_sample_type */
#define PERF_ATTR_SIZE_VER3 96  /* add: sample_regs_user */
                                /* add: sample_stack_user */
#define PERF_ATTR_SIZE_VER4 104 /* add: sample_regs_intr */
#define PERF_ATTR_SIZE_VER5 112 /* add: aux_watermark */
```
"""
struct PerfEventAttr

    # This is taken from Linxu Kernel 4.10.0

    typeid              ::PerfTypeId    # original field name: 'type'
    size                ::UInt32        # set this to sizeof(PerfEventAttr) for the kernel to check compatibility issues
    config              ::UInt64        # Depending on `typeid`, one of the given unions
    sample              ::UInt64
    sample_type         ::UInt64
    read_format         ::PerfEventReadFormat
    flags               ::UInt64
    wakeup              ::UInt32
    bp_type             ::UInt32
    bp_addr             ::UInt64
    bp_len              ::UInt64
    branch_sample_type  ::UInt64
    sample_regs_user    ::UInt64
    sample_stack_user   ::UInt32
    clockid             ::Int32
    sample_regs_intr    ::UInt64
    aux_watermark       ::UInt32
    sample_max_stack    ::UInt16
    __reserved_2        ::UInt16

end

CpuCycleCounter() = PerfEventAttr( PERF_TYPE_HARDWARE
                                 , sizeof(PerfEventAttr)
                                 , PERF_COUNT_HW_CPU_CYCLES
                                 , 0, 0, PERF_FORMAT_UNKNOWN
                                 , 0b0110_0001, 0, 0, 0, 0
                                 , 0, 0, 0, 0, 0, 0, 0, 0)

sizeof(PerfEventAttr) != 112 &&
    warning("Type 'struct Perf.PerfEventAttr' does not meet kernel definition.")

"""
Linux syscall identifier for perf_event_open.

Keep in sync with `<linux-headers>/include/uapi/asm-generic/unistd.h`
"""
const __NR_perf_event_open = 298

#
#  Kernel API Wrapper
#

"""
Open and initialize the Linux kernel's performance monitor device.  This
basically wraps the low-level interface in a way like the kernel's API is
defined.
"""
function perf_event_open( attr::Ref{PerfEventAttr}
                        , pid::Int32, cpu::Int32
                        , group_fd::Int32, flags::UInt32
                        ) ::Int32

    ccall( :syscall, Int32
         , ( Int32, Ref{PerfEventAttr}, Int32, Int32, Int32, UInt32 )
         , __NR_perf_event_open, attr, pid, cpu, group_fd, flags )

end

#function test

end # module KernelAPI


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Model of basic streams
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt and Rafael de P. Soares
#* $Id$
#*---------------------------------------------------------------------
type vapour_stream
	vapour_stream()=begin
		PP=outers.PP
		new(
			stream(),
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			[
				:(h = PP.VapourEnthalpy(T, P, z)),
				:(v = 1),
			],
			[
				"Vapour Enthalpy","Vapour stream",
			],
			[:PP,],
		)
	end
	_base_1::stream
	PP::DanaPlugin
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export vapour_stream
function setEquationFlow(in::vapour_stream)
	addEquation(1)
	addEquation(2)
end
function atributes(in::vapour_stream,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Vapour Material Stream"
	fields[:Info]="Model for vapour material streams.
	This model should be used only when the phase of the stream
	is known ''a priori''."
	drive!(fields,_)
	return fields
end
vapour_stream(_::Dict{Symbol,Any})=begin
	newModel=vapour_stream()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
module Parallization

    using ..DataManager
    using MPI
    using Formatting

    export JobDistributionInfo, getYsplitInfoByRank, syncField!
    export printJobDistributionInfo

    mutable struct YSplitInfo
        pull_fr_rng      :: UnitRange
        pull_to_rng      :: UnitRange
        push_fr_rng      :: UnitRange
        push_to_rng      :: UnitRange

        pull_fr_rng_bnd  :: AbstractArray{Union{UnitRange, Nothing}}
        pull_to_rng_bnd  :: AbstractArray{Union{UnitRange, Nothing}}
        push_fr_rng_bnd  :: AbstractArray{Union{UnitRange, Nothing}}
        push_to_rng_bnd  :: AbstractArray{Union{UnitRange, Nothing}}
    end

    mutable struct JobDistributionInfo

        overlap       :: Int64
        nworkers      :: Int64
        wranks        :: Array{Int64, 1} # worker ranks
        wrank_to_idx  :: Dict{Int64, Int64}
        y_split_infos :: AbstractArray{YSplitInfo, 1}

        function JobDistributionInfo(;
            nworkers :: Int64,
            Ny       :: Int64,
            overlap  :: Int64 = 3,
        )
            if nworkers == 0
                throw(ErrorException("No available workers!"))
            end

            wranks = collect(1:nworkers)

            (
                pull_fr_rngs,
                pull_to_rngs,
                push_fr_rngs,
                push_to_rngs,
                pull_fr_rngs_bnd,
                pull_to_rngs_bnd,
                push_fr_rngs_bnd,
                push_to_rngs_bnd

            ) = calParallizationRange(N=Ny, P=nworkers, L=overlap)

            wrank_to_idx = Dict()
            for i = 1:nworkers
                wrank_to_idx[i] = i
            end

            y_split_infos = Array{YSplitInfo}(undef, nworkers)

            for (i, p) in enumerate(wranks)

                y_split_infos[i] = YSplitInfo(
                    pull_fr_rngs[i],
                    pull_to_rngs[i],
                    push_fr_rngs[i],
                    push_to_rngs[i],
                    pull_fr_rngs_bnd[i, :],
                    pull_to_rngs_bnd[i, :],
                    push_fr_rngs_bnd[i, :],
                    push_to_rngs_bnd[i, :],
                )

            end

            return new(
                overlap,
                nworkers,
                wranks,
                wrank_to_idx,
                y_split_infos,
            )

        end

    end

    function calParallizationRange(;
        N = Integer,     # Total grids
        P = Integer,     # Number of procs
        L = Integer,     # Overlapping grids
    )

        if ! (N >= max(1, L) * P)
            throw(ErrorException("Condition must be satisfied: N >= max(1, L) * P"))
        end

        n̄ = floor(Integer, N / P)
        R = N - n̄ * P

        # "Here" is slave, "there" is master.
        # So "pull" means receiving data from master.
        # and "push" means sending data to master.

        pull_fr_rngs = Array{Union{UnitRange, Nothing}}(undef, P)
        pull_to_rngs = Array{Union{UnitRange, Nothing}}(undef, P)
        push_fr_rngs = Array{Union{UnitRange, Nothing}}(undef, P)
        push_to_rngs = Array{Union{UnitRange, Nothing}}(undef, P)

        # 1: lower latitude side (south), 2: higher latitude side (north)
        pull_fr_rngs_bnd = Array{Union{UnitRange, Nothing}}(undef, P, 2)
        pull_to_rngs_bnd = Array{Union{UnitRange, Nothing}}(undef, P, 2)
        push_fr_rngs_bnd = Array{Union{UnitRange, Nothing}}(undef, P, 2)
        push_to_rngs_bnd = Array{Union{UnitRange, Nothing}}(undef, P, 2)

        cnt = 1
        for p = 1:P
            m = (p <= R) ? n̄ + 1 : n̄  # assigned grids

            pull_fr_rngs[p] = cnt-L:cnt+m-1+L
            pull_to_rngs[p] = 1:length(pull_fr_rngs[p])
            push_fr_rngs[p] = L+1:L+m
            push_to_rngs[p] = cnt:cnt+m-1

            # Boundary
            pull_fr_rngs_bnd[p, 1] = cnt-L:cnt-1
            pull_fr_rngs_bnd[p, 2] = cnt+m:cnt+m+L-1

            pull_to_rngs_bnd[p, 1] = 1:L
            pull_to_rngs_bnd[p, 2] = L+m+1:L+m+L

            push_fr_rngs_bnd[p, 1] = L+1:L+L
            push_fr_rngs_bnd[p, 2] = L+m-L+1:L+m

            push_to_rngs_bnd[p, 1] = cnt:cnt+L-1
            push_to_rngs_bnd[p, 2] = cnt+m-L:cnt+m-1

            cnt += m
        end

        # South pole and north pole do not have boundaries
        pull_fr_rngs_bnd[1, 1] = nothing
        pull_to_rngs_bnd[1, 1] = nothing
        push_fr_rngs_bnd[1, 1] = nothing
        push_to_rngs_bnd[1, 1] = nothing

        pull_fr_rngs_bnd[end, 2] = nothing
        pull_to_rngs_bnd[end, 2] = nothing
        push_fr_rngs_bnd[end, 2] = nothing
        push_to_rngs_bnd[end, 2] = nothing

        # Adjust the first and last range (south pole and north pole)
        pull_fr_rngs[1] = (pull_fr_rngs[1][1]+L):pull_fr_rngs[1][end]
        pull_fr_rngs[end] = pull_fr_rngs[end][1]:(pull_fr_rngs[end][end]-L)

        pull_to_rngs[1] = 1:length(pull_fr_rngs[1])
        pull_to_rngs[end] = 1:length(pull_fr_rngs[end])

        push_fr_rngs[1] = 1:length(push_fr_rngs[1])

        # Change range because to southmost boundary is trimmed
        if P > 1
            pull_to_rngs_bnd[1, 2] = pull_to_rngs_bnd[1, 2] .- L
            push_fr_rngs_bnd[1, 2] = push_fr_rngs_bnd[1, 2] .- L
        end

        return pull_fr_rngs,
               pull_to_rngs,
               push_fr_rngs,
               push_to_rngs,
               pull_fr_rngs_bnd,
               pull_to_rngs_bnd,
               push_fr_rngs_bnd,
               push_to_rngs_bnd

    end

    function getYsplitInfoByRank(
        jdi :: JobDistributionInfo,
        rank :: Integer,
    )
        return jdi.y_split_infos[jdi.wrank_to_idx[rank]]
    end

    mutable struct SyncInfo
        vars         :: AbstractArray{DataUnit, 1}
        y_split_info :: AbstractArray{DataUnit, 1}
    end

    function syncField!(
        vars         :: AbstractArray{DataUnit},
        jdi          :: JobDistributionInfo,
        direction    :: Symbol,
        sync_type    :: Symbol,
    )

        comm = MPI.COMM_WORLD
        rank = MPI.Comm_rank(comm)

        is_master = (rank == 0)

        reqs = Array{MPI.Request}(undef,0)

        if direction == :S2M  # Slave to master
            if sync_type == :BLOCK
                if is_master
                    for (i, _rank) in enumerate(jdi.wranks)
                        for (j, var) in enumerate(vars)
                            v = view(var.data, :, :, getYsplitInfoByRank(jdi, _rank).push_to_rng)
                            push!(reqs, MPI.Irecv!(v, i, j, comm))
                        end
                    end
                else
                    for (j, var) in enumerate(vars)
                        v = view(var.data, :, :, getYsplitInfoByRank(jdi, rank).push_fr_rng)
                        push!(reqs, MPI.Isend(v, 0, j, comm))
                    end
                end
            end

            if sync_type == :BND
                if is_master
                    for (i, _rank) in enumerate(jdi.wranks)
                        for (k, push_to_rng_bnd) in enumerate(getYsplitInfoByRank(jdi, _rank).push_to_rng_bnd)
                            (push_to_rng_bnd == nothing) && continue
                            for (j, var) in enumerate(vars)
                                v = view(var.data, :, :, push_to_rng_bnd)
                                push!(reqs, MPI.Irecv!(v, i, k*length(vars) + j, comm))
                            end
                        end
                    end
                else
                    for (k, push_fr_rng_bnd) in enumerate(getYsplitInfoByRank(jdi, rank).push_fr_rng_bnd)
                        (push_fr_rng_bnd == nothing) && continue
                        for (j, var) in enumerate(vars)
                            v = view(var.data, :, :, push_fr_rng_bnd)
                            push!(reqs, MPI.Isend(v, 0, k*length(vars) + j, comm))
                        end
                    end
                end
            end

        elseif direction == :M2S  # Master to slave

            if sync_type == :BLOCK
                if is_master
                    for (i, _rank) in enumerate(jdi.wranks)
                        for (j, var) in enumerate(vars)
                            v = view(var.data, :, :, getYsplitInfoByRank(jdi, _rank).pull_fr_rng)
                            push!(reqs, MPI.Isend(v, i, j, comm))
                        end
                    end
                else
                    for (j, var) in enumerate(vars)
                        v = view(var.data, :, :, getYsplitInfoByRank(jdi, rank).pull_to_rng)
                        push!(reqs, MPI.Irecv!(v, 0, j, comm))
                    end
                end
            end

            if sync_type == :BND
                if is_master
                    for (i, _rank) in enumerate(jdi.wranks)
                        for (k, pull_fr_rng_bnd) in enumerate(getYsplitInfoByRank(jdi, _rank).pull_fr_rng_bnd)
                            (pull_fr_rng_bnd == nothing) && continue
                            for (j, var) in enumerate(vars)
                                v = view(var.data, :, :, pull_fr_rng_bnd)
                                push!(reqs, MPI.Isend(v, i, k * length(vars) + j, comm))
                            end
                        end
                    end
                else
                    for (k, pull_to_rng_bnd) in enumerate(getYsplitInfoByRank(jdi, rank).pull_to_rng_bnd)
                        (pull_to_rng_bnd == nothing) && continue
                        for (j, var) in enumerate(vars)
                            v = view(var.data, :, :, pull_to_rng_bnd)
                            push!(reqs, MPI.Irecv!(v, 0, k*length(vars) + j, comm))
                        end
                    end
                end
            end
        end

        MPI.Waitall!(reqs)

        MPI.Barrier(comm)
        #=
        if direction == :M2S && rank==1
            for (_, var) in enumerate(vars)
                if var.id == "TEMP"
                    println("detect TEMP.")
                    if var.sdata1[1, 1, end] == 0
                        println("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
                    end
                end
            end
        end  # Master to slave
        =#
    end

    function printJobDistributionInfo(jdi :: JobDistributionInfo)
        println(format("overlap  = {:d}", jdi.overlap))
        println(format("nworkers = {:d}", jdi.nworkers))
        println(format("wranks   = {:s}", string(jdi.wranks)))
        println(format("wrank_to_idx = {:s}", string(jdi.wrank_to_idx)))

        for (i, y_split_info) in enumerate(jdi.y_split_infos)
            println(format("[{:d}] pull_fr_rng = {:s}", i, string(y_split_info.pull_fr_rng)))
            println(format("[{:d}] pull_to_rng = {:s}", i, string(y_split_info.pull_to_rng)))
            println(format("[{:d}] push_fr_rng = {:s}", i, string(y_split_info.push_fr_rng)))
            println(format("[{:d}] push_to_rng = {:s}", i, string(y_split_info.push_to_rng)))

            for j = 1:length(y_split_info.pull_fr_rng_bnd)
                println(format("[{:d}] pull_fr_rng_bnd[{:d}] = {:s}", i, j, string(y_split_info.pull_fr_rng_bnd[j])))
                println(format("[{:d}] pull_to_rng_bnd[{:d}] = {:s}", i, j, string(y_split_info.pull_to_rng_bnd[j])))
                println(format("[{:d}] push_fr_rng_bnd[{:d}] = {:s}", i, j, string(y_split_info.push_fr_rng_bnd[j])))
                println(format("[{:d}] push_to_rng_bnd[{:d}] = {:s}", i, j, string(y_split_info.push_to_rng_bnd[j])))
            end
        end
    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
###
# asserts for constructor of JuDGEModel
###

# Definition of an important set of constraints
function NormalConstraints()
    [
        (GenericAffExpr{Float64,VariableRef}, MathOptInterface.LessThan{Float64}),
        (GenericAffExpr{Float64,VariableRef}, MathOptInterface.GreaterThan{Float64}),
        (GenericAffExpr{Float64,VariableRef}, MathOptInterface.EqualTo{Float64}),
        (VariableRef, MathOptInterface.LessThan{Float64}),
        (VariableRef, MathOptInterface.GreaterThan{Float64}),
        (VariableRef, MathOptInterface.EqualTo{Float64}),
    ]
end

function check_specification_is_legal(sub_problems)
    message = "JuDGE Specification Error: "
    if !same_expansions_at_each_node(sub_problems)
        message *= "Every subproblem must have the same expansion variables"
        return error(message)
    end
    if !same_shutdowns_at_each_node(sub_problems)
        message *= "Every subproblem must have the same shutdown variables"
        return error(message)
    end
    if !check_objectives(sub_problems)
        message *= "The subproblem objective must be minimizing and either be of type VariableRef or GenericAffExpr"
        return error(message)
    end

    check_costs(sub_problems)
    check_constraints(sub_problems) #must be at the end of this subroutine!!
end

#### Expansion variables at every node have to be the "same"
function same_expansions_at_each_node(subproblems)
    expansion_structures = collect(map(get_structure_of_expansions, values(subproblems)))
    all(map(x -> x == expansion_structures[1], expansion_structures))
end

function same_shutdowns_at_each_node(subproblems)
    shutdown_structures = collect(map(get_structure_of_shutdowns, values(subproblems)))
    all(map(x -> x == shutdown_structures[1], shutdown_structures))
end

function check_objectives(subproblems)
    for (n,sp) in subproblems
        if !(objective_function_type(sp) <: GenericAffExpr) && objective_function_type(sp) != VariableRef
            return false
        elseif objective_sense(sp)!=MOI.MIN_SENSE
            return false
        end
    end
    true
end

function check_costs(subproblems)
    for (n,sp) in subproblems
        check_sp_costs(sp)
    end
end

function check_sp_costs(model)
    if haskey(model.ext, :capitalcosts) && model.ext[:capitalcosts][:constant]!=0
        #if typeof(model.ext[:capitalcosts])!=AffExpr
        #    error("@capitalcosts must be provided a linear expression (AffExpr)")
        #elseif model.ext[:capitalcosts][:constant]!=0
        error("@capitalcosts should not contain any constant terms")
        # else
        #     av=all_variables(model)
        #     for v in av
        #         found=false
        #         for e in model.ext[:expansions]
        #             if v==e[2] || v in e[2]
        #                 found=true
        #                 break
        #             end
        #         end
        #         if !found && v in keys(model.ext[:capitalcosts].terms)
        #             error("@capitalcosts should only contain expansion variables"*string(v))
        #         end
        #     end
        # end
    elseif haskey(model.ext, :ongoingcosts) && model.ext[:ongoingcosts][:constant]!=0
        #if typeof(model.ext[:ongoingcosts])!=AffExpr
        #    error("@ongoingcosts must be provided a linear expression (AffExpr)")
        #elseif model.ext[:ongoingcosts].constant!=0
        error("@ongoingcosts should not contain any constant terms")
        # else
        #     av=all_variables(model)
        #     for v in av
        #         found=false
        #         for e in model.ext[:expansions]
        #             if v==e[2] || v in e[2]
        #                 found=true
        #                 break
        #             end
        #         end
        #         if !found && v in keys(model.ext[:ongoingcosts].terms)
        #             error("@ongoingcosts should only contain expansion variables")
        #         end
        #     end
        # end
    end
    nothing
end

function get_structure_of_expansions(subproblem)
    expkeys = get_expansion_keys(subproblem)

    function thinwrapper(var)
        if typeof(var) <: AbstractArray
            collect(eachindex(var))
        else
            ()
        end
    end

    return Dict(i => thinwrapper(subproblem.obj_dict[i]) for i in expkeys)
end

function get_structure_of_shutdowns(subproblem)
    shutkeys = get_shutdown_keys(subproblem)

    function thinwrapper(var)
        if typeof(var) <: AbstractArray
            collect(eachindex(var))
        else
            ()
        end
    end

    return Dict(i => thinwrapper(subproblem.obj_dict[i]) for i in shutkeys)
end

function get_expansion_keys(model)
    filter(keys(model.obj_dict)) do key
        for (exp,var) in model.ext[:expansions]
            if var === model.obj_dict[key] && model.ext[:options][exp][1]==:expansion
                return true
            end
        end
        return false
    end
end

function get_shutdown_keys(model)
    filter(keys(model.obj_dict)) do key
        for (exp,var) in model.ext[:expansions]
            if var === model.obj_dict[key] && model.ext[:options][exp][1]==:shutdown
                return true
            end
        end
        return false
    end
end

function check_constraints(subproblems)
    for (n,sp) in subproblems
        check_sp_constraints(sp)
    end
end

function check_sp_constraints(model)
    loct=list_of_constraint_types(model)
    expansion_keys=get_expansion_keys(model)
    shutdown_keys=get_shutdown_keys(model)

    all_expansion_variables=Array{VariableRef,1}()
    all_shutdown_variables=Array{VariableRef,1}()

    for exp_key in expansion_keys
         var=model[exp_key]
         if model.ext[:options][exp_key][4]!=:Con
             if typeof(var)==VariableRef
                 push!(all_expansion_variables,var)
             elseif typeof(var) <: AbstractArray
                 for i in eachindex(var)
                     push!(all_expansion_variables,var[i])
                 end
             end
         end
    end

    for shut_key in shutdown_keys
        var=model[shut_key]
        if typeof(var)==VariableRef
            push!(all_shutdown_variables,var)
        elseif typeof(var) <: AbstractArray
            for i in eachindex(var)
                push!(all_shutdown_variables,var[i])
            end
        end
    end

    warnings = ["Positive coefficient for shutdown variable on LHS of >= constraint",
                "Positive coefficient for expansion variable on LHS of <= constraint",
                "Expansion or shutdown variable in == constraint",
                "Negative coefficient for expansion variable on LHS of >= constraint",
                "Negative coefficient for shutdown variable on LHS of <= constraint"]

    status = [false,false,false,false,false]

    for ct in loct
        for con in all_constraints(model,ct[1],ct[2])
            con_obj=JuMP.constraint_object(con)
            if typeof(con_obj.func)==VariableRef
                # if typeof(con_obj.set)==MathOptInterface.GreaterThan{Float64}
                #     if con_obj.func in all_shutdown_variables
                #         status[1]=true
                #     end
                # elseif typeof(con_obj.set)==MathOptInterface.LessThan{Float64}
                #     if con_obj.func in all_expansion_variables
                #        status[2]=true
                #     end
                if typeof(con_obj.set)==MathOptInterface.EqualTo{Float64}
                    if con_obj.func in all_expansion_variables || con_obj.func in all_shutdown_variables
                        status[3]=true
                    end
                end
            elseif typeof(con_obj.func) <: GenericAffExpr
                if typeof(con_obj.set)==MathOptInterface.GreaterThan{Float64}
                    for (v,c) in con_obj.func.terms
                        if c>0.0 && v in all_shutdown_variables
                            status[1]=true
                        elseif c<0.0 && v in all_expansion_variables
                            status[4]=true
                        end
                    end
                elseif typeof(con_obj.set)==MathOptInterface.LessThan{Float64}
                    for (v,c) in con_obj.func.terms
                        if c<0.0 && v in all_shutdown_variables
                            status[5]=true
                        elseif c>0.0 && v in all_expansion_variables
                            status[2]=true
                        end
                    end
                elseif typeof(con_obj.set)==MathOptInterface.EqualTo{Float64}
                    for (v,c) in con_obj.func.terms
                        if v in all_shutdown_variables || v in all_expansion_variables
                            status[3]=true
                        end
                    end
                end
            elseif typeof(con_obj.func) <: GenericQuadExpr
                LHS=QuadExpr(AffExpr(0.0),OrderedDict{UnorderedPair{VariableRef},Float64}())
                for (v,c) in con_obj.func.terms
                    if v.a in all_shutdown_variables || v.b in all_shutdown_variables || v.a in all_expansion_variables || v.b in all_expansion_variables
                        error("JuDGE Specification Error: Expansion or shutdown variable in quadratic term")
                    end
                end

                if typeof(con_obj.set)==MathOptInterface.GreaterThan{Float64}
                    for (v,c) in con_obj.func.aff.terms
                        if c>0.0 && v in all_shutdown_variables
                            status[1]=true
                        elseif c<0.0 && v in all_expansion_variables
                            status[4]=true
                        end
                    end
                elseif typeof(con_obj.set)==MathOptInterface.LessThan{Float64}
                    for (v,c) in con_obj.func.aff.terms
                        if c<0.0 && v in all_shutdown_variables
                            status[5]=true
                        elseif c>0.0 && v in all_expansion_variables
                            status[2]=true
                        end
                    end
                elseif typeof(con_obj.set)==MathOptInterface.EqualTo{Float64}
                    for (v,c) in con_obj.func.aff.terms
                        if v in all_shutdown_variables || v in all_expansion_variables
                            status[3]=true
                        end
                    end
                end
            elseif typeof(con_obj.func) == Array{GenericAffExpr{Float64,VariableRef},1}
                for aff in con_obj.func
                    for (v,c) in aff.terms
                        if v in all_shutdown_variables || v in all_expansion_variables
                            error("JuDGE Specification Error: Expansion or shutdown variable in indicator constraint")
                        end
                    end
                end
            else
                @warn("JuDGE Specification Error: Unable to verify constraint of type "*string(typeof(con_obj.func)))
            end
        end
    end

    for i in 1:5
        if status[i]
            @warn(warnings[i])
        end
    end

    nothing
end


text:
@with_kw mutable struct PedestrianAvoidancePOMDPFrenet{B <: Updater} <: DriverModel{LatLonAccel}

    timestep::Float64 = 0.2
    t_current::Float64 = 0
    tick::Int64 = 0

    ego_vehicle::Vehicle = Vehicle(VehicleState(VecSE2(0., 0., 0.), 0.), VehicleDef(), 1)
    a::LatLonAccel = LatLonAccel(0.0, 0)
    desired_velocity::Float64 = 50.0 / 3.6

    obstacles::Vector{ConvexPolygon} = []
    sensor_observations::Vector{Vehicle} = []
    objects_tracked = Dict{Int64, Float64}()

    update_tick_high_level_planner::Int64 = 1

    env::CrosswalkEnv = CrosswalkEnv(CrosswalkParams())
    sensor::AutomotiveSensors.GaussianSensor = AutomotiveSensors.GaussianSensor(AutomotiveSensors.LinearNoise(10, 0., 0.),
                                                               AutomotiveSensors.LinearNoise(10, 0., 0.), 0, 0, MersenneTwister(1))

    pomdp::SingleOCFPOMDP = SingleOCFPOMDP()
    policy::AlphaVectorPolicy{SingleOCFPOMDP,SingleOCFAction} = AlphaVectorPolicy(pomdp, Vector{Vector{Float64}}())
    updater::B = SingleOCFUpdater(pomdp)

    policy_dec::DecPolicy = DecPolicy(policy, pomdp, policy.action_map, (x,y) -> min.(x,y))

    b::SingleOCFBelief = SingleOCFBelief([],[])
    b_dict::Dict{Int64, SingleOCFBelief} = Dict{Int64, SingleOCFBelief}()

    prediction_obstacle::Array{Float64} = []

    brake_request::Bool = false
    collision_rate::Float64 = 0.0
    risk::Float64 = 0.0
    ttc::Float64 = 0.0

end

function AutomotiveDrivingModels.observe!(model::PedestrianAvoidancePOMDPFrenet, scene::Scene, roadway::Roadway, egoid::Int)

    ego = scene[findfirst(egoid, scene)]
    model.pomdp.ego_vehicle = ego
    model.ego_vehicle = ego
    #model.sensor_observations = measure(PerfectSensor(), ego, scene, roadway, model.obstacles)

    # get observations from sensor, only visible objects with sensor noise
    model.sensor_observations = measure(model.sensor, ego, scene, roadway, model.obstacles)

    # add sensor delay 0.2s
    model.sensor_observations = objects_time_delay(model)

    # initialization of the belief for the absent state
    if (model.t_current == 0 )
        model.b_dict[PEDESTRIAN_OFF_KEY] = initBeliefAbsentPedestrian(model.pomdp, ego.state.posF.t, ego.state.v)
    end

    ################ High Level Planner ###################################################
    if (  true ) #model.tick % model.update_tick_high_level_planner == 0 )

        #println("--------------------------POMDP high level planner----------------------- t: ", model.t_current)
        #println("EGO: x/y:", ego.state.posG.x, " / ",  ego.state.posG.y, " v: ", ego.state.v)

        # prepare observation dictionary / remove observations outside state space and add absent state
        observations = get_observations_state_space(model, ego, model.sensor_observations)

        # update belief dictionary
        b_new = update(model.pomdp, model.updater, model.b_dict, SingleOCFAction(model.a.a_lon, model.a.a_lat), observations)
        model.b_dict = deepcopy(b_new)
        #println("dict: ", model.b_dict)

#=
b_dict_tmp = Dict{Int64, SingleOCFBelief}()
for oid in keys(model.b_dict)
    if ( oid == PEDESTRIAN_OFF_KEY )
        b_states = []
        b_prob = []
        for (s, prob) in weighted_iterator(model.b_dict[PEDESTRIAN_OFF_KEY] )
            if ( is_state_absent(model.pomdp, s) )
                (ego_y_state_space, ego_v_state_space) = getEgoDataInStateSpace(model.pomdp, ego.state.posF.t, ego.state.v)
                s = SingleOCFState(ego_y_state_space, ego_v_state_space, -10.0, -10.0, model.pomdp.PED_THETA_RANGE[1] ,0.0)
            end
            push!(b_states, s)
            push!(b_prob, prob)
        end
        normalize!(b_prob, 1)
        b_dict_tmp[PEDESTRIAN_OFF_KEY] =  SingleOCFBelief(b_states, b_prob)
    else
        b_dict_tmp[oid] = model.b_dict[oid]
    end
end
=#
        # use policy and belief dictionary to calculate next action
        act = action(model.policy_dec, model.b_dict)
        model.a = LatLonAccel(act.lateral_movement, act.acc)
        #println("action combined: ", model.a)

        # dummy functionality to test transition function / belief update
        #=
        if (model.tick > 2 )
           model.a = LatLonAccel(1.0, -2.0)
            println("manual intervention")
        end
        =#

    end
    model.tick += 1
    model.t_current = model.t_current + model.timestep
end

function AutomotiveDrivingModels.propagate(veh::Vehicle, action::LatLonAccel, roadway::Roadway, Δt::Float64)

    # new velocity
    v_ = veh.state.v + action.a_lon*Δt
    v_ = clamp(v_, 0., v_)

    # lateral offset
    delta_y = action.a_lat * Δt   # a_lat corresponds to lateral velocity --> a_lat == v_lat
    if v_ <= 0.
        delta_y = 0.
    end
    y_ = veh.state.posG.y + delta_y
    y_ = clamp(y_, -1.0, 1.0)

    s_new = v_ * Δt

    # longitudional distance based on required velocity and lateral offset
    if (abs(delta_y) < abs(s_new) )
        delta_x = sqrt(s_new^2 - delta_y^2 )
    else
        delta_x = 0
    end

    delta_x = clamp(delta_x, 0., delta_x)
    x_ = veh.state.posG.x + delta_x

    return VehicleState(VecSE2(x_, y_, veh.state.posG.θ), roadway, v_)
end

function AutomotiveDrivingModels.get_name(model::PedestrianAvoidancePOMDPFrenet)
    return "Pedestrian Avoidance System (POMDP)"
end

AutomotiveDrivingModels.rand(rng::AbstractRNG, model::PedestrianAvoidancePOMDPFrenet) = model.a

@with_kw mutable struct ObservationCallback

    sensor_observations::Vector{Vector{Vehicle}}
    ego_vehicle::Vector{Vehicle}
    ego_a::Vector{Float64}
    collision::Vector{Bool}

    b_dict::Vector{Dict{Int64, SingleOCFBelief}}
    action::Vector{SingleOCFAction}

    collision_rate::Vector{Float64}
    ttc::Vector{Float64}
    risk::Vector{Float64}
    brake_request::Vector{Bool}
    prediction_obstacle::Vector{Array{Float64}}

end

function AutomotiveDrivingModels.run_callback(
        callback::ObservationCallback,
        scenes::Vector{Scene},
        actions::Union{Nothing, Vector{Frame{A}}},
        roadway::R,
        models::Dict{I,M},
        tick::Int) where {A<:EntityAction, S,D,I,R,M<:DriverModel}

    push!(callback.sensor_observations, models[1].sensor_observations)
    push!(callback.ego_vehicle, models[1].ego_vehicle)
    push!(callback.ego_a, models[1].a.a_lon)

    collision = is_crash(scenes[tick])
    push!(callback.collision, collision)

    push!(callback.b_dict, deepcopy(models[1].b_dict))
    push!(callback.action, SingleOCFAction(models[1].a.a_lon, models[1].a.a_lat))

    push!(callback.collision_rate, models[1].collision_rate)
    push!(callback.ttc, models[1].ttc)
    push!(callback.risk, models[1].risk)
    push!(callback.brake_request, models[1].brake_request)
    push!(callback.prediction_obstacle, models[1].prediction_obstacle)

    return collision
end


text:
"""
	runSimulation(β0,cs0,cd0,dChars,
				outDummy,dV,nCons,nProdTot,nProdCons,nA0;
				seed = 24)

Runs simulation given parameter values. Returns the generated dataset and the results of FI, SO and W estimation.

# Arguments:
- `β0`, `cs0`, `cd0`: True parameter values with which data is generated in SD model
- `dChars`:		Distribution of characteristics in X
- `outDummy`:   Bool whether to include a dummy for the outside option (if true, first element in β0 is preference parameter for it)
- `dV`:		 	Distribution of errors v. If all zeros, set to `Normal(0,0)`.
- `nCons`:		# consumers to simulate
- `nProdTot`:   # of total products (from which a random subsample is shown to each consumer)
- `nProdCons`:  # of products shown to each consumer
- `nA0`:		Size of initial awareness (if nA0 = 2, all consumers initially observe the first two products without having to discover them. But they still have to inspect them.)

# Optional arguments:
- `seed`: Sets seed throughout
- `solver`: Option to set solver through either `Optim` or `NLopt` package. Array with first element string indicating package, and second element algorithm (e.g. `["Optim",Optim.NelderMead()"]` or `["NLopt",:LN_BOBYQA]`)
"""
function runSimulation(
	β0,
	cs0,
	cd0,
	dChars,
	outDummy,
	dV,
	nCons,
	nProdTot,
	nProdCons,
	nA0;
	seed = 24,
	 # uses default estimation options (defined in definitions_struct.jl)
	options = estimOptions()
)

	###########
	# SD model
	SD = SearchDiscovery(
		β = β0,
		cs = cs0,
		cd = cd0,
		dze = 0.0,
		dChars = dChars,
		outDummy = outDummy,
		seed = seed,
		nd = 1,
		nA0 = nA0,
		dV = dV,
		dE = Normal(),
		cfun = "linear",
		zdfun = "linear"
	)
	# Generate data
	d = genData(SD, nCons, nProdTot, nProdCons)

	# # Calculate welfare changes
	wSD = calcWelfareChanges(SD,d; seed = seed)

	###########
	# Full information model
	FI = FullInfo(
		β = β0,
		dChars = dChars,
		outDummy = outDummy,
		seed = seed + 1,
		dV = dV,
		dE = Normal(),
		smo = 3,
		options=options
	)
	# For NLopt algorithms, replace here since need number of parameters
	if typeof(options.algo) <: Opt
		FI.options.algo = Opt(options.algo.algorithm,length(β0))
	end

	resFI = estimateModel(FI, d, startvals = β0)
	# input estimated parameters
	FI.β = resFI[1]
	# Calculate welfare changes
	wFI = calcWelfareChanges(FI,d; seed = seed+4)

	# Calcualte summary stats
	dFI = calcSummaryStats(d,FI)

	###########
	# Estimate search in order model
	RS = RandomSearch(
		β = β0,
		cs = cs0,
		dChars = dChars,
		outDummy = outDummy,
		seed = seed + 2,
		dV = dV,
		dE = Normal(),
		cfun = "linear",
		smo = 10,
		options=options,
	)
	if typeof(options.algo) <: Opt
		RS.options.algo = Opt(options.algo.algorithm,length(β0)+1)
	end

	resRS = estimateModel(RS, d, startvals = vcat(β0, cs0))

	# Input estimated pars
	RS.β = resRS[1]
	RS.cs = resRS[2]
	# Calculate welfare changes
	wRS = calcWelfareChanges(RS,d; seed = seed+4)

	# Calcualte summary stats
	dRS = calcSummaryStats(d,RS)

	###########
	#  Estimate Weitzman model without ranking effects
	Wf = Weitzman(
		β = β0,
		cs = cs0,
		cd = 0,
		nd = 1,
		nA0 = nA0,
		dChars = dChars,
		outDummy = outDummy,
		seed = seed + 3,
		dV = dV,
		dE = Normal(),
		cfun = "linear",
		smo = 10,
		options=options,
	)
	if typeof(options.algo) <: Opt
		Wf.options.algo = Opt(options.algo.algorithm,length(β0)+1)
	end
	# Change structure of data
	# d.pos / d.pid data different in Weitzman
	d1 = (seli=d.seli,pid = vcat(d.pid...), pos = vcat(d.pos...),
			chars = d.chars, path = d.path, purch = d.purch)

	# Run estimation
	resWf = estimateModel(Wf, d1, startvals = vcat(β0, cs0))

	# Input estimated pars
	Wf.β = resWf[1]
	Wf.cs = resWf[2]
	# Calculate welfare changes
	wWf = calcWelfareChanges(Wf,d; seed = seed+4)

	# Calcualte summary stats
	dWf = calcSummaryStats(d,Wf)

		###########
	#  Estimate Weitzman model with ranking effects
	W = Weitzman(
		β = β0,
		cs = cs0,
		cd = cd0,
		nd = 1,
		nA0 = nA0,
		dChars = dChars,
		outDummy = outDummy,
		seed = seed + 3,
		dV = dV,
		dE = Normal(),
		cfun = "linear",
		smo = 10,
		options=options,
	)
	if typeof(options.algo) <: Opt
		W.options.algo = Opt(options.algo.algorithm,length(β0)+2)
	end
	# Change structure of data
	# d.pos / d.pid data different in Weitzman
	d1 = (seli=d.seli,pid = vcat(d.pid...), pos = vcat(d.pos...),
			chars = d.chars, path = d.path, purch = d.purch)

	# Run estimation
	resW = estimateModel(W, d1, startvals = vcat(β0, cs0,cd0))

	# Input estimated pars
	W.β = resW[1]
	W.cs = resW[2]
	W.cd = resW[3]
	# Calculate welfare changes
	wW = calcWelfareChanges(W,d; seed = seed+4)

	# Calcualte summary stats
	dW = calcSummaryStats(d,W)

	return d, wSD , (resFI,wFI,dFI),  (resRS,wRS,dRS), (resWf,wWf,dWf), (resW,wW,dW)
end

"""
	calcWelfareChanges(m::Model, d; seed = 511)
Calculates the welfare changes specified in the paper for a given model specification. Returns matrix with different changes in rows, and different measures in cols.

*Rows*: 1. Baseline welfare | 2. No search and discovery costs | 3. 1% 'price' decrease of product on 5th position | 4. 1% 'price' decrase of product on 15th position

*Columns*: 1. Welfare | 2. Utility | 3. Costs | 4. Demand Outside | 5. Demand Position 1 | 6. Demand Position 2 ....

# Arguments:
- `m`: `Model` definition
- `d`: Data input

# Optional arguments:
- `seed`: Sets seed for random draws in welfare calculation

"""
function calcWelfareChanges(m::Model, d; seed = 511)
	nCons = length(d.seli)
	# Make draws (stays same across changes)
	Random.seed!(seed)
	e = rand(m.dE, size(d.chars, 1), m.options.nDraws*5) # use more draws than in estimation for precision
	v = rand(m.dV, size(d.chars, 1), m.options.nDraws*5)

	d0 = 	if typeof(m) <: Weitzman # Changes structure of data for Weitzman model
				(seli=d.seli,pid = vcat(d.pid...), pos = vcat(d.pos...),
					chars = d.chars, path = d.path, purch = d.purch)
			else
				d
			end

	# Initial welfare
	w0 = calcWelfareDemand(m, d0; e , v)

	# Create array to store values
	out = fill(NaN,4,3+length(w0[4]))

	# Fill in values from welfare
	fillOut!(out,w0,1)

	##############################
	# Remove all search costs
	mn = FullInfo(β=m.β,dChars=m.dChars,outDummy=m.outDummy,
					dV=m.dV,dE=m.dE,smo=0) #smoothing does not matter here

	w1 = calcWelfareDemand(mn, d; e , v)
	fillOut!(out,w1,2)

	##############################
	# 1% 'price' change of product shown on 5th position (for all consumers)
	dn = deepcopy(d0)
	iProd = findfirst( x-> x==5,d.pos[1])

	for ii in eachindex(dn.seli)
		if m.outDummy
			dn.chars[dn.seli[ii][iProd],end-1] = dn.chars[dn.seli[ii][iProd],end-1] * 0.99
		else
			dn.chars[dn.seli[ii][iProd],end] = dn.chars[dn.seli[ii][iProd],end] * 0.99
		end
	end

	w1 = calcWelfareDemand(m, dn; e , v)
	fillOut!(out,w1,3)

	##############################
	# 1% 'price' decrease of product shown on 15th position (for all consumers)
	dn = deepcopy(d0)
	iProd = findfirst( x-> x == 15,d.pos[1])

	for ii in eachindex(dn.seli)
		if m.outDummy
			dn.chars[dn.seli[ii][iProd],end-1] = dn.chars[dn.seli[ii][iProd],end-1] * 0.99
		else
			dn.chars[dn.seli[ii][iProd],end] = dn.chars[dn.seli[ii][iProd],end] * 0.99
		end
	end

	w1 = calcWelfareDemand(m, dn; e , v)
	fillOut!(out,w1,4)

	return out
end

"""
	calcSummaryStats(d,m::Model)
"""
function calcSummaryStats(d,m::Model)
	Random.seed!(m.seed + 2934)
	nCons = length(d.seli)

	nS = zeros(2) # Avg. number of searches
	shareOutside = zeros(2) # Share outside option

	# Calculate summary statistics in data
	nS[1] = mean(length.(d.path))
	shareOutside[1] = count(d.purch .== 1) / nCons

	# Simulate summary statistics (across m.options.nDraws draws)
	nSsim, shareOutsidesim = _calcSummaryStats(d,m)
	nS[2] = nSsim
	shareOutside[2] = shareOutsidesim

	return nS, shareOutside
end

function _calcSummaryStats(d,m::FullInfo)
	Random.seed!(m.seed+120937)
	nCons = length(d.seli)

	nS = NaN
	shareOutside = 0.0
	xβv = zeros(size(d.chars,1))
	u 	= similar(xβv)
	for e in 1:m.options.nDraws
		xβv = d.chars * m.β .+ rand(m.dV,size(d.chars,1))
		u   = xβv .+ rand(m.dE,size(d.chars,1))
		purch = genSearchPaths(m,d,u)
		shareOutside += count(purch .== 1) / nCons / m.options.nDraws
	end

	return nS, shareOutside
end

function _calcSummaryStats(d,m::RandomSearch)
	Random.seed!(m.seed+120937)
	nCons = length(d.seli)
	nProdCons = length(d.seli[1])

	# Get values need to generate search paths
	cfun = getCfun(m.cfun)
	μ  = sum(mean(m.dChars) .* m.β[1:end-m.outDummy]) + mean(m.dV)
	σ2  = sum(var(m.dChars) .* m.β[1:end-m.outDummy] .^2) + var(m.dV)
	zs = calcZs(0.0,Normal(μ,sqrt(σ2)),cfun(m.cs))[1]

	xβv = zeros(size(d.chars,1))
	u 	= similar(xβv)

	nS = 0.0
	shareOutside = 0.0

	for e in 1:m.options.nDraws
		xβv = d.chars * m.β .+ rand(m.dV,size(d.chars,1))
		u   = xβv .+ rand(m.dE,size(d.chars,1))
		path,purch = genSearchPaths(m,d,nProdCons,zs,u)
		nS += mean(length.(path)) / m.options.nDraws
		shareOutside += count(purch .== 1) / nCons / m.options.nDraws
	end

	return nS, shareOutside

end

function _calcSummaryStats(d,m::Weitzman)
	Random.seed!(m.seed+120937)
	nCons = length(d.seli)
	nProdCons = length(d.seli[1])

	# Change data structure for Weitzman model
	d1 = (seli=d.seli,pid = vcat(d.pid...), pos = vcat(d.pos...),
				chars = d.chars, path = d.path, purch = d.purch)

	# Get values need to generate search paths
	cfun = 	getCfun(m.cfun)
	ξ   = 	[calcZs(0.0,m.dE,cfun(m.cs + m.cd * pp))[1] for pp in 0:nProdCons]

	# Pre-allocate
	u = zeros(size(d.chars,1))
	zs = similar(u)
	e = similar(u)
	v = similar(u)

	nS = 0.0
	shareOutside = 0.0

	for k in 1:m.options.nDraws
		rand!(m.dV,v)
		rand!(m.dE,e)
		fillUZS!(u,zs,m,d1,m.β,ξ,e,v)

		path,purch,_ = genSearchPaths(m,d1,nProdCons,zs,u)
		nS +=  countSearches(m,path) / nCons / m.options.nDraws
		shareOutside += count(purch .== 1) / nCons / m.options.nDraws
	end

	return nS, shareOutside
end

function countSearches(m::Weitzman,path)
	nS = 0
	for i in 1:size(path,2)
		_nS = findfirst( x -> x==0, path[:,i])
		nS += isnothing(_nS) ? length(path) -1  : _nS -1
	end
	return nS
end

"""
	gatherResults(β0,cs0,cd0, wSD,outFI,outRS,outW)
Gather results for table and analysis.
"""
function gatherResults(fname::String)
	@load fname wSD outFI outRS outWf outW β0 cs0 cd0
	res = [outWf,outW,outRS,outFI]

	##########################
	# Estimated parameters
	b = [m[1][1] for m in res]
	b = hcat(b...)
	b = hcat(β0,b)

	cs = [m[1][2] for m in res[1:end-1]]
	cs = hcat(cs0,cs',NaN)

	cd = hcat(cd0,NaN,outW[1][3],NaN,NaN)

	L = [NaN]
	for m in res
		push!(L,-m[1][end])
	end
	estimOut = vcat(b,cs,cd,L')

	##########################
	# Welfare changes
	dwU = hcat(wSD[1,2],diffFirst(wSD)[:,2]')
	dwC = hcat(wSD[1,3],diffFirst(wSD)[:,3]')
	dw  = dwU .- dwC
	for m in res
		wdiffd = diffFirst(m[2])
		dwU = vcat(dwU,hcat(m[2][1,2],wdiffd[:,2]'))
		dwC = vcat(dwC,hcat(m[2][1,3],wdiffd[:,3]'))
		dw  = vcat(dw,dwU[end,:]' .- dwC[end,:]')
	end
	# Normalize by second coefficient  (only matters for estimated welfare, not changes)
	for i = 1:size(dw,1), m in [dw,dwU,dwC]
		m[i,:] .= m[i,:] ./ abs(b[2,i])
	end

	# Put into % changes of initial welfare
	for m in (dw,dwU,dwC)
		m[:,2:end] .= m[:,2:end] ./ abs.(m[:,1]) * 100
	end

	welfareOut = (dw,dwU,dwC)

	##########################
	# Demand changes
	dDemand5 = wSD[:,10]'
	dDemand1 = wSD[:,5]'
	for r in res
		dDemand5 = vcat(dDemand5,r[2][:,10]')
		dDemand1 = vcat(dDemand1,r[2][:,5]')
	end

	# Make into percentage changes
	dDemand5 = (dDemand5[:,2:end] .- dDemand5[:,1] ) ./ dDemand5[:,1] .* 100
	dDemand1 = (dDemand1[:,2:end] .- dDemand1[:,1] ) ./ dDemand1[:,1] .* 100

	demandOut = (dDemand1,dDemand5)

	##########################
	# Extract size of search set and share of outside Option

	dStats = [outFI[3][1][1] outFI[3][2][1]]  # Note, first element always from generated data
	for r in res
		dStats = vcat(dStats,[r[3][1][2] r[3][2][2]])
	end

	return estimOut, welfareOut, demandOut, dStats
end

##############################################################

function fillOut!(out,w,k)
	for i in 1:size(out,2)
		out[k,i] = 	if i >= 4
						w[4][i-3] # demand for each product (1. is outside)
					else
						w[i]
					end
	end
end

"""
	diffFirst(a::Array{T,2}) where T<:Real
Takes difference to first row for each row in array `a`.
"""
function diffFirst(a::Array{T,2}) where T<:Real
	return hcat([a[i,:] .- a[1,:] for i = 2:size(a,1)]...)'
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export uipage

"""
    uipage(;kwargs...)
    uipage(children::Any;kwargs...)
    uipage(children_maker::Function;kwargs...)

An UIPage component.
ExampleComponent is an example component.
It takes a property, `label`, and
displays it.
It renders an input with the property `value`
which is editable by the user.
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional): Children
- `id` (String; optional): Element id
- `name` (String; required): Page name
- `value` (String; required)
"""
function uipage(; kwargs...)
        available_props = Symbol[:children, :id, :name, :value]
        wild_props = Symbol[]
        return Component("uipage", "UIPage", "dash_carbon_components", available_props, wild_props; kwargs...)
end

uipage(children::Any; kwargs...) = uipage(;kwargs..., children = children)
uipage(children_maker::Function; kwargs...) = uipage(children_maker(); kwargs...)

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components


text:
using DPC, CairoMakie
#include(joinpath(@__DIR__, "utils.jl"))
include("utils.jl")

config_and = """
refractory_duration: 5.01
inputs:
- id: i11
- id: i12
- id: i13
- id: i14
- id: i15
- id: i21
- id: i22
- id: i23
- id: i24
- id: i25
- id: i31
- id: i32
- id: i33
- id: i34
- id: i35
neurons:
- id: n
  θ_syn: 5
  θ_seg: 1
  branches:
    - id: seg0
      θ_syn: 5
      θ_seg: 2
      branches:
          - id: seg1
            θ_syn: 5
          - id: seg2
            θ_syn: 5
synapses:
- {id: syn11, source: i11, target: seg1}
- {id: syn12, source: i12, target: seg1}
- {id: syn13, source: i13, target: seg1}
- {id: syn14, source: i14, target: seg1}
- {id: syn15, source: i15, target: seg1}
- {id: syn21, source: i21, target: seg2}
- {id: syn22, source: i22, target: seg2}
- {id: syn23, source: i23, target: seg2}
- {id: syn24, source: i24, target: seg2}
- {id: syn25, source: i25, target: seg2}
- {id: syn31, source: i31, target: seg0}
- {id: syn32, source: i32, target: seg0}
- {id: syn33, source: i33, target: seg0}
- {id: syn34, source: i34, target: seg0}
- {id: syn35, source: i35, target: seg0}
"""

config_or = """
refractory_duration: 5.01
inputs:
- id: i11
- id: i12
- id: i13
- id: i14
- id: i15
- id: i21
- id: i22
- id: i23
- id: i24
- id: i25
- id: i31
- id: i32
- id: i33
- id: i34
- id: i35
neurons:
- id: n
  θ_syn: 0
  θ_seg: 1
  branches:
    - id: seg0
      θ_syn: 5
      θ_seg: 1
      branches:
        - id: seg1
          θ_syn: 5
        - id: seg2
          θ_syn: 5
synapses:
- {id: syn11, source: i11, target: seg1}
- {id: syn12, source: i12, target: seg1}
- {id: syn13, source: i13, target: seg1}
- {id: syn14, source: i14, target: seg1}
- {id: syn15, source: i15, target: seg1}
- {id: syn21, source: i21, target: seg2}
- {id: syn22, source: i22, target: seg2}
- {id: syn23, source: i23, target: seg2}
- {id: syn24, source: i24, target: seg2}
- {id: syn25, source: i25, target: seg2}
- {id: syn31, source: i31, target: seg0}
- {id: syn32, source: i32, target: seg0}
- {id: syn33, source: i33, target: seg0}
- {id: syn34, source: i34, target: seg0}
- {id: syn35, source: i35, target: seg0}
"""

config_chain = """
refractory_duration: 5.01
inputs:
- id: i11
- id: i12
- id: i13
- id: i14
- id: i15
- id: i21
- id: i22
- id: i23
- id: i24
- id: i25
- id: i31
- id: i32
- id: i33
- id: i34
- id: i35
neurons:
- id: n
  θ_syn: 0
  θ_seg: 1
  branches:
    - id: seg0
      θ_syn: 5
      θ_seg: 1
      branches:
        - id: seg2
          θ_syn: 5
          branches:
            - id: seg1
              θ_syn: 5
synapses:
- {id: syn11, source: i11, target: seg1}
- {id: syn12, source: i12, target: seg1}
- {id: syn13, source: i13, target: seg1}
- {id: syn14, source: i14, target: seg1}
- {id: syn15, source: i15, target: seg1}
- {id: syn21, source: i21, target: seg2}
- {id: syn22, source: i22, target: seg2}
- {id: syn23, source: i23, target: seg2}
- {id: syn24, source: i24, target: seg2}
- {id: syn25, source: i25, target: seg2}
- {id: syn31, source: i31, target: seg0}
- {id: syn32, source: i32, target: seg0}
- {id: syn33, source: i33, target: seg0}
- {id: syn34, source: i34, target: seg0}
- {id: syn35, source: i35, target: seg0}
"""

volleys = [(15.0, :i1), (75.0, :i2), (140.0, :i3), (215.0, :i2), (275.0, :i1), (300.0, :i3), (415.0, :i2), (475.0, :i3)]
spikes = [ (time=t+τ-0.5, name=Symbol("$(pop)$(i)")) for (t, pop) in volleys for (i,τ) in enumerate(5*rand(5))]

## Set up plot

fig = Figure(resolution = (0.75textwidth, 0.75textwidth))
ax_left  = fig[1,1] = Axis(fig; backgroundcolor=:transparent)
ax_right = fig[1,3] = Axis(fig; backgroundcolor=:transparent)
gl = fig[1,2] = GridLayout(alignmode=Outside())
ax_chain = gl[1,1] = Axis(fig;
    yticks=(2.5:5:15, ["C","B","A"]),
    yminorticks=Makie.IntervalsBetween(5, true),
    yminorticksvisible = true, yminorgridvisible = true,
    titlealign=:left, title = "a.    Sequential segments"
)
ax_or    = gl[2,1] = Axis(fig;
    yticks=(2.5:5:15, ["C","B","A"]),
    yminorticks=Makie.IntervalsBetween(5, true),
    yminorticksvisible = true, yminorgridvisible = true,
    titlealign=:left, title = "b.    Parallel segments (either required)"
)
ax_and   = gl[3,1] = Axis(fig;
    yticks=(2.5:5:15, ["C","B","A"]),
    yminorticks=Makie.IntervalsBetween(5, true),
    yminorticksvisible = true, yminorgridvisible = true,
    titlealign=:left, title = "c.    Parallel segments (both required)", xlabel="time [ms]"
)

linkaxes!(ax_chain,ax_and)
linkaxes!(ax_chain,ax_or)

hidexdecorations!(ax_chain, grid=false)
hidexdecorations!(ax_or, grid=false)
hidedecorations!(ax_left)
hidedecorations!(ax_right)

fig

## Simulate each motive
offset = Dict(:chain => Point2f0(0,-0.5), :and => Point2f0(0,-1.35) , :or => Point2f0(0,-2.7))

for (name, config, ax) in ((:chain,config_chain,ax_chain), (:or,config_or,ax_or), (:and,config_and,ax_and))
    (net,objects) = load_network(YAML_source=config)
    input=[ Event(:input_spikes, 0.0, t, objects[name]) for (t, name) in spikes]
    logger=simulate!(net, input)

    # get traces
    syn11 = get_trace(:syn11, logger.data)
    syn12 = get_trace(:syn12, logger.data)
    syn13 = get_trace(:syn13, logger.data)
    syn14 = get_trace(:syn14, logger.data)
    syn15 = get_trace(:syn15, logger.data)
    syn21 = get_trace(:syn21, logger.data)
    syn22 = get_trace(:syn22, logger.data)
    syn23 = get_trace(:syn23, logger.data)
    syn24 = get_trace(:syn24, logger.data)
    syn25 = get_trace(:syn25, logger.data)
    syn31 = get_trace(:syn31, logger.data)
    syn32 = get_trace(:syn32, logger.data)
    syn33 = get_trace(:syn33, logger.data)
    syn34 = get_trace(:syn34, logger.data)
    syn35 = get_trace(:syn35, logger.data)
    seg0  = get_trace(:seg0, logger.data)
    seg1  = get_trace(:seg1, logger.data)
    seg2  = get_trace(:seg2, logger.data)
    n     = get_trace(:n, logger.data)

    plateau0_starts = filter(x->(x.object == :seg0 && x.event == :plateau_starts), logger.data).t
    plateau1_starts = filter(x->(x.object == :seg1 && x.event == :plateau_starts), logger.data).t
    plateau2_starts = filter(x->(x.object == :seg2 && x.event == :plateau_starts), logger.data).t
    plateau1_ends = filter(x->(x.object == :seg1 && x.event == :plateau_ends), logger.data).t
    plateau2_ends = filter(x->(x.object == :seg2 && x.event == :plateau_ends), logger.data).t
    spike_times = filter(x->(x.object == :n && x.event == :spikes), logger.data).t
    # plot dynamics

    statetrace!(ax, [0;seg0.t;500], [DPC.voltage_low;seg0.state;DPC.voltage_low], Dict(DPC.voltage_low=>:transparent, DPC.voltage_elevated=>color_3_25, DPC.voltage_high=>color_3_50), 0 ,  4.9 )
    statetrace!(ax, [0;seg2.t;500], [DPC.voltage_low;seg2.state;DPC.voltage_low], Dict(DPC.voltage_low=>:transparent, DPC.voltage_elevated=>color_2_25, DPC.voltage_high=>color_2_50), 5 ,   4.9 )
    statetrace!(ax, [0;seg1.t;500], [DPC.voltage_low;seg1.state;DPC.voltage_low], Dict(DPC.voltage_low=>:transparent, DPC.voltage_elevated=>color_1_25, DPC.voltage_high=>color_1_50), 10 ,  4.9 )

    # plot spikes
    for (i, syn) in enumerate([syn11,syn12,syn13,syn14,syn15])
        statetrace!(ax, [0;syn.t;500], [0;syn.state;0], Dict(0=>:transparent, 1=>color_1), 9 .+ i ,  0.9 )
    end

    for (i, syn) in enumerate([syn21,syn22,syn23,syn24,syn25])
        statetrace!(ax, [0;syn.t;500], [0;syn.state;0], Dict(0=>:transparent, 1=>color_2), 4 .+ i ,  0.9 )
    end

    for (i, syn) in enumerate([syn31,syn32,syn33,syn34,syn35])
        statetrace!(ax, [0;syn.t;500], [0;syn.state;0], Dict(0=>:transparent, 1=>color_3), -1 .+ i ,  0.9 )
    end

    # plot output spikes
    # lines!.(ax, Rect.(plateau0_starts.-5.5, -0.25, 11, 5.5), color=:gray10, linewidth=2)

    # plot neurons
    plot!(
        name == :and ? ax_right : ax_left,
        objects[:n],
        angle_between = 20/180*pi,
        branch_length = 0.5,
        branch_width = 0.05,
        root_position = offset[name],
        color=Dict(:n=>color_3_25, :seg0=>color_3, :seg1=>color_1, :seg2=>color_2),
    )
end

fig

## Finish up and save
colsize!(fig.layout, 1, Relative(0.1))
colsize!(fig.layout, 3, Relative(0.1))
ylims!(ax_left, (-3.1,1.1))
ylims!(ax_right, (-3.1,1.1))
xlims!(ax_left, (-0.2,0.2))
xlims!(ax_right, (-0.2,0.2))
# xlims!(ax_input, (0,500))

save(joinpath("figures","motifs.pdf"), fig)
save(joinpath("figures","motifs.svg"), fig)
save(joinpath("figures","motifs.png"), fig)
fig