text:
struct Occu <: UnmarkedModel
  data::UmData
  opt::UnmarkedOpt
  submodels::NamedTuple
end

"""
    occu(ψ_form::FormulaTerm, p_form::FormulaTerm, data::Umdata)

Fit single-season occupancy models. Covariates on occupancy and detection
are specified with `ψ_form` and `p_form`, respectively. The `UmData`
object should contain `site_covs` and `obs_covs` `DataFrames` containing
columns with names matching the formulas.
"""
function occu(ψ_form::FormulaTerm, p_form::FormulaTerm, data::UmData)

  occ = UmDesign(:Occupancy, ψ_form, LogitLink(), data.site_covs)
  det = UmDesign(:Detection, p_form, LogitLink(), data.obs_covs)
  add_idx!([occ, det])
  np = get_np([occ, det])

  y = data.y
  N, J = size(y)
  nd = ndetects(y)

  function loglik(β::Array)

    ψ = transform(occ, β)
    p = transform(det, β)

    ll = zeros(eltype(β), N)

    idx = 0
    for i = 1:N
      if all(ismissing.(y[i,:]))
        idx += J
        continue
      end
      cp = 1.0
      for j = 1:J
        idx += 1
        if ismissing(y[i,j])
          continue
        end
        cp *= p[idx]^y[i,j] * (1-p[idx])^(1-y[i,j])
      end
      ll[i] = log(ψ[i] * cp + (1-ψ[i]) * (nd[i]==0))
    end

    return -sum(ll)
  end

  opt = optimize_loglik(loglik, np)
  Occu(data, opt, (occ=UnmarkedSubmodel(occ,opt),
                   det=UnmarkedSubmodel(det,opt)))
end

"""
    occu(ψ_formula::Union{Array,FormulaTerm},
         p_formula::Union{Array,FormulaTerm}, data::UmData)

Provide multiple formulas for ψ and/or p and fit models for all formula
combinations. To fit all model subsets, wrap formula with `allsubs()`.
"""
function occu(ψ_form::Union{Array,FormulaTerm},
              p_form::Union{Array,FormulaTerm}, data::UmData)
  cf = combine_formulas(ψ_form, p_form)
  out = Array{UnmarkedModel}(undef, length(cf))
  p = Progress(length(cf), desc=string("Fitting ", length(cf), " models "))
  Threads.@threads for i in 1:length(cf)
    out[i] = occu(cf[i][1], cf[i][2], data)
    next!(p)
  end
  return UnmarkedModels(out)
end

# Simulations------------------------------------------------------------------

"Simulate new dataset from a fitted model"
function simulate(fit::Occu)

  ψ = predict(occupancy(fit))
  p = predict(detection(fit))

  y = sim_y(Occu, ψ, p)
  rep_missing!(y, fit.data.y)

  UmData(y, fit.data.site_covs, fit.data.obs_covs)
end

"Simulate new dataset from provided formulas"
function simulate(::Type{Occu}, ψ_form::FormulaTerm, p_form::FormulaTerm,
                  ydims::Tuple{Int,Int}, coef::Array)

  sc = gen_covs(ψ_form, ydims[1])
  oc = gen_covs(p_form, prod(ydims))

  occ = UmDesign(:Occ, ψ_form, LogitLink(), sc)
  det = UmDesign(:Det, p_form, LogitLink(), oc)
  add_idx!([occ, det])
  np = get_np([occ, det])

  if np != length(coef) error(string("Coef array must be length ",np)) end

  ψ = transform(occ, coef)
  p = transform(det, coef)

  y = sim_y(Occu, ψ, p)

  UmData(y, sc, oc)
end

#Simulate y matrix from provided psi and p vectors"
function sim_y(::Type{Occu}, ψ::Array, p::Array)

  N = length(ψ)
  J = Integer(length(p) / N)
  z = map(x -> rand(Binomial(1, x))[1], ψ)
  y = Array{Union{Int,Missing}}(undef, N, J)
  idx = 0
  for i = 1:N
    for j = 1:J
      idx += 1
      y[i,j] = rand(Binomial(1, p[idx] * z[i]))[1]
    end
  end

  return y
end

# Update ----------------------------------------------------------------------

function update(fit::Occu, data::UmData)
  occu(occupancy(fit).formula, detection(fit).formula, data)
end

# Goodness-of-fit -------------------------------------------------------------

#Get cohorts with unique NA placement
function get_cohorts(fit::Occu)
  data = fit.data
  na_pattern = mapslices(x -> string(ismissing.(x)*1...), data.y, dims=2)[:,1]
  cohorts = string.(keys(countmap(na_pattern)))

  out = []
  for i in 1:length(cohorts)
    #Discard cohort where all obs are NAs
    if cohorts[i] == repeat("1", size(data.y)[2]) continue end
    idx = findall(na_pattern .== cohorts[i])
    push!(out, data[idx])
  end

  return out
end

#Get probability of a particular encounter history
function get_prob_eh(::Type{Occu}, eh::String, pvec::Array, psi::Float64)
  ehv = map(x -> x == "missing" ? missing : parse(Int, x), split(eh, " "))
  pvals = (ehv .* pvec) + (1 .- ehv) .* (1 .- pvec)
  pvals = collect(skipmissing(pvals))

  if sum(skipmissing(ehv)) == 0
    return psi * prod(pvals) + (1-psi)
  end

  return psi * prod(pvals)
end

#Get counts for each observed encounter history
function get_obs_counts(fit::Occu, data::UmData)
  eh_strings = mapslices(x -> join(x, " "), data.y, dims=2)
  return countmap(eh_strings)
end

#Get expected counts of each observed encounter history
function get_exp_counts(fit::Occu, data::UmData)
  N,J = size(data.y)

  obs_counts = get_obs_counts(fit, data)

  eh_obs = collect(keys(obs_counts))
  nobs = length(eh_obs)

  psi = predict(occupancy(fit), data.site_covs)
  p = predict(detection(fit), data.obs_covs)

  counts_expect = fill(0.0, nobs)
  for i = 1:nobs
    pstart = 1

    for j = 1:N
      pend = pstart + J - 1
      counts_expect[i] += get_prob_eh(Occu, eh_obs[i], p[pstart:pend], psi[j])
      pstart += J
    end
  end

  return Dict(zip(eh_obs, counts_expect))
end

#Compute MacKenzie-Bailey (2004) chi-square statistic
function mb_chisq(fit::Occu)

  #Divide dataset into cohorts based on NA placement
  cohorts = get_cohorts(fit::Occu)

  chisq_all = fill(0.0, length(cohorts))
  for i in 1:length(cohorts)

    N = size(cohorts[i].y)[1]

    obs_dict = get_obs_counts(fit, cohorts[i])
    exp_dict = get_exp_counts(fit, cohorts[i])

    neh = length(obs_dict)

    obs = collect(values(obs_dict))
    exp = collect(values(exp_dict))

    #For observed encounter histories
    chisq_prt1 = sum(map(x -> (obs[x] - exp[x])^2 / exp[x], 1:neh))

    #For unobserved possible encounter histories
    chisq_prt2 = maximum([0, N - sum(exp)])

    chisq_all[i] = chisq_prt1 + chisq_prt2
  end

  return sum(chisq_all)
end

"""
    gof(fit::Occu, nsims::Int=50)

Compute the MacKenzie-Bailey (2004) goodness-of-fit test for a
single-season occupancy model. The distribution of the test statistic
is generated using parametric bootstrapping with `nsims` simulations.
"""
function gof(fit::Occu, nsims::Int=50)
  test_name = "MacKenzie-Bailey Goodness-of-fit"
  t0 = mb_chisq(fit)
  tstar = parboot(fit, nsims, mb_chisq)
  return UmGOF(test_name,t0, tstar)
end

function Base.show(io::IO, gof::UmGOF)
  pval = mean(gof.tstar .> gof.t0)
  chat = gof.t0 / mean(gof.tstar)
  println()
  println(gof.test_name)
  println()
  println(string("  χ2 = ", @sprintf("%.4f",round(gof.t0, digits=4))))
  println(string("  P-value = ", @sprintf("%.4f",round(pval, digits=4))))
  print(string("  Est. c-hat = ",
               @sprintf("%.4f", round(chat, digits=4))))
end

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


text:
platforms_and_sources = Dict(
    Platform("x86_64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.2.0.53/11.3_04222021/cudnn-11.3-linux-x64-v8.2.0.53.tgz",
                      "7a195dc93a7cda2bdd4d9b73958d259c784be422cd941a9a625aab75309f19dc")],
    Platform("powerpc64le", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.2.0.53/11.3_04222021/cudnn-11.3-linux-ppc64le-v8.2.0.53.tgz",
                      "cfe06735671a41a5e25fc7542d740177ac8eab1ab146bd30f19e0fa836895611")],
    Platform("x86_64", "windows") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.2.0.53/11.3_04222021/cudnn-11.3-windows-x64-v8.2.0.53.zip",
                      "c8f108f8b646fd22ab35da7a7d73e2d977f671a1abd3364d372292102e0d7429")],
)

ext: jl

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
# Replace terminal/non-terminal with
# Traits.  Not sure if this is a win.

abstract type LarkSymbol end

abstract type TerminalTrait end
struct IsTerminal <: TerminalTrait end
struct IsNotTerminal <: TerminalTrait end

# == Symbols == #
Base.:(==)(s1::LarkSymbol,s2::LarkSymbol) = begin
    typeof(s1) == typeof(s2) && s1.name == s2.name
end

Base.hash(s1::LarkSymbol,h::UInt64) = begin
    hash(s1.name,h)
end

Base.show(io::IO,ls::LarkSymbol) = print(io,"$(typeof(ls))($(ls.name))")
Base.isless(s1::LarkSymbol,s2::LarkSymbol) = s1.name < s2.name
# == Terminals == #
struct Terminal <: LarkSymbol
    name::String
    filter_out::Bool
end

Terminal(name;filter_out=false) = Terminal(name,filter_out)
TerminalTrait(::Type{Terminal}) = IsTerminal()
Base.show(io::IO,ls::Terminal) = print(io,"$(ls.name)")

# == Non terminals == #
struct NonTerminal <: LarkSymbol
    name::String
    #is_term::Bool
end

TerminalTrait(::Type{NonTerminal}) = IsNotTerminal()

Base.show(io::IO,ls::NonTerminal) = print(io,"$(ls.name)")

is_terminal(l::T) where {T} = is_terminal(TerminalTrait(T),l)
is_terminal(::IsTerminal,l) = true
is_terminal(::IsNotTerminal,l) = false

# == Rule options == #

struct RuleOptions
    keep_all_tokens::Bool
    expand1::Bool
    priority::Union{Int,Nothing}
    template_source
    empty_indices
end

RuleOptions(;keep_all_tokens=false,expand1=false,priority=nothing,template_source=nothing,empty_indices=()) = RuleOptions(keep_all_tokens,expand1,priority,template_source,empty_indices)

# A copy-but-change-keep-all-tokens rule
RuleOptions(ro::RuleOptions,b::Bool) = RuleOptions(b,ro.expand1,ro.priority,ro.template_source,ro.empty_indices)

Base.show(io::IO,ro::RuleOptions) = begin
    print(io,"RuleOptions($(ro.keep_all_tokens),$(ro.expand1), $(ro.priority), $(ro.template_source))")
end

mutable struct Rule
    origin::NonTerminal
    expansion::Array{LarkSymbol}
    order::Int32
    alias::Union{Nothing,String}
    options::RuleOptions
end

Rule(origin,expansion;order=0,alias=nothing,options=nothing) = begin
    Rule(origin,expansion,order,alias,options === nothing ? RuleOptions() : options)
end

Base.show(io::IO,r::Rule) = print(io,"<$(r.origin.name):  $(join(r.expansion," "))>")

# Equality, needed for indexing into rules

Base.:(==)(r1::Rule,r2::Rule) = begin
    return r1.origin == r2.origin && r1.expansion == r2.expansion
end

# The rule name should be unique enough
Base.hash(r1::Rule,h::UInt64) = hash(r1.origin,h)


text:
function jacobian(m::KrusellSmith)

    # Load in endogenous state and eq cond indices
    endo = augment_model_states(m.endogenous_states_unnormalized,
                                n_model_states_unnormalized(m))
    eq   = m.equilibrium_conditions

    # Load in grid settings
    nw::Int = get_setting(m, :nw)
    ns::Int = get_setting(m, :ns)

    zhi::Float64 = get_setting(m, :zhi)
    zlo::Float64 = get_setting(m, :zlo)
    smoother::Float64 = get_setting(m, :smoother)

    wgrid = get_grid(m, :wgrid)
    sgrid = get_grid(m, :sgrid)

    # Make the Jacobian
    JJ = zeros(4*nw+2,8*nw+4)

    # Create auxiliary variables
    g::Array{Float64,1} = trunc_lognpdf(sgrid, m[:μ_s].value, m[:σ_s].value)                   # density g (of skill, s) evaluated on sgrid

    m[:Wstar] = MPL(1.0, m[:Lstar].value, m[:Kstar].value, m[:α].value)            # Steady state wages
    m[:Rstar] = MPK(1.0, m[:Lstar].value, m[:Kstar].value, m[:α].value, m[:δ].value)     # Steady state net return on capital

    # KF equation
    KFmollificand::Array{Float64, 2} = (repeat(wgrid.points, 1, nw*ns) - m[:Rstar]*repeat(wgrid.points'-m[:cstar].value', nw, ns))/m[:Wstar] - kron(sgrid.points',ones(nw, nw))
    μRHS::Array{Float64,1} = mollifier.(KFmollificand, zhi, zlo, smoother)*kron(sgrid.weights.*(g/m[:Wstar]), wgrid.weights.*m[:μstar].value)

    # These matrices (from the Euler equation) correspond to the matrices in the PDF documentation on Dropbox
    # (note that Γ is Λ removing the Σ_i^nx, and the ω^x_i terms from the summation.
    ξ = Array{Float64, 2}(undef, nw, nw) #zeros(nw, nw)
    Γ = Array{Float64, 2}(undef, nw, nw) #zeros(nw, nw)

    # These matrices (from the KF equation) corresp. to the lyx document on dropbox, script A and ζ
    AA = Array{Float64,2}(undef, nw, nw)
    ζ = Array{Float64, 2}(undef, nw, nw)

    # i_w is the ith entry in the nw discretization of x
    # i_wp is i_w' the grid of potential future x's
    for i_w in 1:nw
        for i_wp in 1:nw
            # The following hold the sums over i=1̣̣̣̣, ..., ns
            sum_ns_ξ::Float64 = 0.0
            sum_ns_Γ::Float64 = 0.0
            sum_ns_AA::Float64 = 0.0
            sum_ns_ζ::Float64 = 0.0
            # The first part of ξ and Γ which only relies on i_w' so we calculate outside of 1:ns loop
            front = ((m[:cstar].value[i_wp].^(-m[:γ]))/m[:Wstar])
            first_part_mollificand = (wgrid.points[i_wp] - m[:Rstar] * (wgrid.points[i_w] -
                                                                        m[:cstar].value[i_w]))/m[:Wstar]
            for iss in 1:ns
                mollificand::Float64 = first_part_mollificand - sgrid.points[iss]
                # dm and mm are the mollifiers which correspond to q(mollificand) in the paper
                dm = dmollifier(mollificand, zhi, zlo, smoother)
                mm = mollifier(mollificand, zhi, zlo, smoother)
                skill_distrib_weights = g[iss] * sgrid.weights[iss]

                sum_ns_ξ += front * dm * skill_distrib_weights
                sum_ns_Γ += front * mm * skill_distrib_weights
                sum_ns_AA+= (1.0/m[:Wstar]) * mm * skill_distrib_weights
                sum_ns_ζ += (1.0/m[:Wstar]) * dm * skill_distrib_weights
            end
            ξ[i_w,i_wp] = sum_ns_ξ
            Γ[i_w,i_wp] = sum_ns_Γ
            AA[i_w,i_wp] = sum_ns_AA
            ζ[i_w,i_wp] = sum_ns_ζ
        end
    end

    dRdK::Float64 = -(1-m[:α])*((m[:Rstar]+m[:δ]-1)/m[:Kstar])
    dWdK::Float64 = (m[:α]*m[:Wstar]/m[:Kstar])
    dRdZ::Float64 = m[:Rstar]+m[:δ]-1.0
    dWdZ::Float64 = m[:Wstar].value
    lRHS = m[:β]*m[:Rstar]*Γ*wgrid.weights

    # Fill in Jacobian

    # Euler equation (EE)
    term1_EE::Array{Float64, 1} = m[:β]*Γ*wgrid.weights - m[:β]*(m[:Rstar]/m[:Wstar]) * ((wgrid.points-m[:cstar].value) .* (ξ*wgrid.weights))
    term2_EE::Array{Float64, 1} = -(1.0/m[:Wstar]) * (lRHS+((m[:β]*m[:Rstar])/m[:Wstar]) * ξ.*(repeat(wgrid.points', nw, 1) - m[:Rstar]*repeat(wgrid.points-m[:cstar].value,1,nw))*wgrid.weights)

    JJ[eq[:eq_euler], endo[:K′_t]]  = term1_EE*dRdK + term2_EE*dWdK
    JJ[eq[:eq_euler], endo[:z′_t]] = term1_EE*dRdZ + term2_EE*dWdZ

    JJ[eq[:eq_euler], endo[:l′_t]] = m[:β]*m[:Rstar]*Γ*Matrix(Diagonal(((wgrid.weights.*(m[:lstar].value.^(-(1.0+m[:γ])/m[:γ]))
                                                                         .*(m[:lstar].value.^(-1.0/m[:γ]) .<= wgrid.points))./m[:cstar].value)))

    JJ[eq[:eq_euler], endo[:l_t]]  = -(eye(nw) + Matrix(Diagonal((m[:β]/m[:γ])*(m[:Rstar]*m[:Rstar]/m[:Wstar])*(ξ*wgrid.weights)
                                                                 .*(m[:lstar].value.^(-(1.0+m[:γ])/m[:γ])).*(m[:lstar].value.^(-1.0/m[:γ]) .<= wgrid.points))))

    # KF Equation
    JJ[eq[:eq_kolmogorov_fwd], endo[:μ_t1]]   = AA'*Matrix(Diagonal(wgrid.weights))

    JJ[eq[:eq_kolmogorov_fwd], endo[:l_t1]] = (-(m[:Rstar]/m[:Wstar])*(1.0/m[:γ])*ζ') * Matrix(Diagonal( m[:μstar].value.*wgrid.weights.*(m[:lstar].value.^(-(1.0+m[:γ])/m[:γ])) .*(m[:lstar].value.^(-1.0/m[:γ]) .<= wgrid.points)))

    term1_KF::Array{Float64,1} = -(1.0/m[:Wstar])*(ζ'.*repeat((wgrid.points-m[:cstar].value)',nw,1))*(m[:μstar].value.*wgrid.weights)
    term2_KF::Array{Float64,1} = -(μRHS/m[:Wstar] + (1.0/(m[:Wstar]*m[:Wstar]))*(ζ'.*(repeat(wgrid.points,1,nw) - m[:Rstar]*repeat( (wgrid.points-m[:cstar].value)',nw,1)))*(m[:μstar].value.*wgrid.weights))

    JJ[eq[:eq_kolmogorov_fwd], endo[:K_t]]  = term1_KF*dRdK + term2_KF*dWdK
    JJ[eq[:eq_kolmogorov_fwd], endo[:z_t]]   = term1_KF*dRdZ + term2_KF*dWdZ

    JJ[eq[:eq_kolmogorov_fwd], endo[:μ_t]]    = -eye(nw)

    # DEFN of LM(t+1) = M(t)
    JJ[eq[:eq_μ], endo[:μ′_t1]]  = eye(nw)

    JJ[eq[:eq_μ], endo[:μ_t]]    = -eye(nw)

    # DEFN of LELL(t+1) = ELL(t)
    JJ[eq[:eq_l], endo[:l′_t1]]= eye(nw)

    JJ[eq[:eq_l], endo[:l_t]]  = -eye(nw)

    # LOM K
    JJ[eq[:eq_K_law_of_motion], endo[:l_t]]  = -(1.0/m[:γ])*m[:μstar].value'*(Matrix(Diagonal(wgrid.weights.*(m[:lstar].value.^(-(1.0+m[:γ])/m[:γ]))
                                                                                .*(m[:lstar].value.^(-1.0/m[:γ]) .<= wgrid.points))))

    JJ[eq[:eq_K_law_of_motion], endo[:μ_t]]  .= -(wgrid.points-m[:cstar].value)'*Matrix(Diagonal(wgrid.weights))

    JJ[eq[:eq_K_law_of_motion], endo[:K′_t]]  .= 1.0

    # TFP
    JJ[eq[:eq_TFP], endo[:z′_t]]   .= 1.0

    JJ[eq[:eq_TFP], endo[:z_t]]    .= -m[:ρ_z]

    if !m.testing && get_setting(m, :normalize_distr_variables)
        JJ = normalize(m, JJ)
    end

    return JJ
end

function normalize(m::KrusellSmith, JJ::Matrix{Float64})
    nw = get_setting(m, :nw)

    P = eye(nw)
    P[:,1] = ones(nw)

    Q = Matrix{Float64}(undef, nw, nw)
    Q::Matrix{Float64}, _ = qr(P)

    S = Matrix{Float64}(undef, nw, nw-1)
    S::Matrix{Float64} = Q[:, 2:end]'

    Qf = zeros(4*nw, 4*nw+2)
    Qf[1:nw, 1:nw] = eye(nw)
    Qf[nw+1:2*nw-1, nw+1:2*nw] = S
    Qf[2*nw:3*nw-2, 2*nw+1:3*nw] = S
    Qf[3*nw-1:4*nw, 3*nw+1:4*nw+2] = eye(nw+2)

    # Qx and Qy are for normalizing any variables that represent distributions
    # the S component is for normalizing a distribution, the other identity portions
    # are for non-distributional variables
    # The distribution in the states, x, is the lagged μ
    # The distribution in the jumps, y, is the current μ
    Qx = cat(S, eye(nw), [1], [1], dims = [1 2])   # 161 x 162
    Qy = cat(S, eye(nw), dims = [1 2])             # 159 x 160

    m <= Setting(:n_predetermined_variables, size(Qx, 1))

    # whenever outputs and/or inputs are densities,
    # pre/postmultiply the jacobians from step 2 by an appropriate
    # matrix so things integrate to 1

    Qright = cat(Qx',Qy',Qx',Qy', dims = [1 2])
    Jac1 = Qf*JJ*Qright

    return Jac1
end


text:
#
# Stupid way to recreate a com so that
# can test com extractors
#

lxd(n) = F.LxDef("\\" * n, 1, F.subs(""))

function lxc(mds)
    F.def_GLOBAL_LXDEFS!()
    empty!(F.GLOBAL_LXDEFS)
    name = match(r"\\(.*?)\{", mds).captures[1]
    F.GLOBAL_LXDEFS[name] = lxd(name)
    tokens = F.find_tokens(mds, F.MD_TOKENS, F.MD_1C_TOKENS)
    blocks, tokens = F.find_all_ocblocks(tokens, F.MD_OCB2)
    lxdefs, tokens, braces, _ = F.find_lxdefs(tokens, blocks)
    lxdefs = cat(collect(values(F.GLOBAL_LXDEFS)), lxdefs, dims=1)
    lxcoms, _ = F.find_lxcoms(tokens, lxdefs, braces)
    return lxcoms[1]
end

F.LOCAL_VARS["sections"] = Pair{String,String}[] => (Pair{String,String},)
F.GLOBAL_VARS["section_counter"] = 1 => (Int,)

isapproxstr(s1, s2) =
    isequal(map(s->replace(s, r"\s|\n"=>""), String.((s1, s2)))...)


text:
import Base.Broadcast: _broadcast_getindex, preprocess, preprocess_args, Broadcasted, broadcast_unalias, combine_axes, broadcast_axes, broadcast_shape, check_broadcast_axes, check_broadcast_shape, throwdm, broadcastable, AbstractArrayStyle, DefaultArrayStyle
import Base: copyto!, tail, axes, length, ndims
struct DiffEqBC{T}
    x::T
end
@inline axes(b::DiffEqBC) = axes(b.x)
@inline length(b::DiffEqBC) = length(b.x)
@inline broadcastable(b::DiffEqBC) = b
@inline Base.ndims(b::Type{DiffEqBC{T}}) where T = ndims(T)
Base.@propagate_inbounds _broadcast_getindex(b::DiffEqBC, i) = _broadcast_getindex(b.x, i)
Base.@propagate_inbounds _broadcast_getindex(b::DiffEqBC{<:AbstractArray{<:Any,0}}, i) = b.x[]
Base.@propagate_inbounds _broadcast_getindex(b::DiffEqBC{<:AbstractVector}, i) = b.x[i[1]]
Base.@propagate_inbounds _broadcast_getindex(b::DiffEqBC{<:AbstractArray}, i) = b.x[i]
diffeqbc(x::Array) = DiffEqBC(x)
diffeqbc(x) = x

# Ensure inlining
@inline combine_axes(A, B) = broadcast_shape(broadcast_axes(A), broadcast_axes(B)) # Julia 1.0 compatible
@inline check_broadcast_axes(shp, A::Union{Number, Array, Broadcasted}) = check_broadcast_shape(shp, axes(A))

@inline preprocess(f, dest, bc::Broadcasted{Style}) where {Style} = Broadcasted{Style}(bc.f, preprocess_args(f, dest, bc.args), bc.axes)
preprocess(f, dest, x) = f(broadcast_unalias(dest, x))

@inline preprocess_args(f, dest, args::Tuple) = (preprocess(f, dest, args[1]), preprocess_args(f, dest, tail(args))...)
@inline preprocess_args(f, dest, args::Tuple{Any}) = (preprocess(f, dest, args[1]),)
preprocess_args(f, dest, args::Tuple{}) = ()

# Performance optimization for the common identity scalar case: dest .= val
@inline copyto!(dest::DiffEqBC, bc::Broadcasted{<:AbstractArrayStyle{0}}) = copyto!(dest.x, bc)
@inline function copyto!(dest::DiffEqBC, bc::Broadcasted)
    axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc))
    dest′ = dest.x
    # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match
    if bc.f === identity && bc.args isa Tuple{AbstractArray} # only a single input argument to broadcast!
        A = bc.args[1]
        if axes(dest) == axes(A)
            return copyto!(dest′, A)
        end
    end
    bcs′ = preprocess(diffeqbc, dest, bc)
    @simd ivdep for I in eachindex(bcs′)
        @inbounds dest′[I] = bcs′[I]
    end
    return dest′ # return the original array without the wrapper
end

# Forcing `broadcasted` to inline is not necessary, since `Vern9` plays well
# with the Base implementation, and `Feagin`s do not use broadcasting.
#
#import Base.Broadcast: broadcasted, combine_styles
#map_nostop(f, t::Tuple{})              = ()
#map_nostop(f, t::Tuple{Any,})          = (f(t[1]),)
#map_nostop(f, t::Tuple{Any, Any})      = (f(t[1]), f(t[2]))
#map_nostop(f, t::Tuple{Any, Any, Any}) = (f(t[1]), f(t[2]), f(t[3]))
#map_nostop(f, t::Tuple)                = (Base.@_inline_meta; (f(t[1]), map_nostop(f,tail(t))...))
#@inline function broadcasted(f::Union{typeof(*), typeof(+), typeof(muladd)}, arg1, arg2, args...)
#    arg1′ = broadcastable(arg1)
#    arg2′ = broadcastable(arg2)
#    args′ = map_nostop(broadcastable, args)
#    broadcasted(combine_styles(arg1′, arg2′, args′...), f, arg1′, arg2′, args′...)
#end

macro ..(x)
    expr = Base.Broadcast.__dot__(x)
    if expr.head in (:(.=), :(.+=), :(.-=), :(.*=), :(./=), :(.\=), :(.^=)) # we exclude `÷=` `%=` `&=` `|=` `⊻=` `>>>=` `>>=` `<<=` because they are for integers
      name = gensym()
      dest = :(DiffEqBase.diffeqbc($(expr.args[1])))
      expr.args[1] = name
      return esc(quote
                  $name = $dest
                  $expr
                 end)
    else
      return esc(expr)
    end
end


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

name = "LibpcapBuilder"
version = v"1.9.0"

# Collection of sources required to build LibpcapBuilder
sources = [
    "http://www.tcpdump.org/release/libpcap-1.9.0.tar.gz" =>
    "2edb88808e5913fdaa8e9c1fcaf272e19b2485338742b5074b9fe44d68f37019",

    #"http://www.linuxfromscratch.org/patches/blfs/svn/libpcap-1.9.0-enable_bluetooth-1.patch" =>
    #"5000de2daec00fb31ccc949f7e265bfb621635c8e2daf9bb17770d75997195b5",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libpcap-1.9.0/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain
make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc, :eabihf),
    Linux(:powerpc64le, :glibc),
    FreeBSD(:x86_64)
]

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

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

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jpsamaroo/LibpcapBuilder

max_issues_repo_name: jpsamaroo/LibpcapBuilder

max_forks_repo_name: jpsamaroo/LibpcapBuilder


text:
using StaticArrays
import Base.length

Vec3 = SVector{3};
## treating vector as location in space
x(v::Vec3) = v[1];
y(v::Vec3) = v[2];
z(v::Vec3) = v[3];
## treating vector as color
r(v::Vec3) = v[1];
g(v::Vec3) = v[2];
b(v::Vec3) = v[3];
## special operators
length(v::Vec3) = sqrt(v[1]^2 + v[2]^2 + v[3]^2);
squared_length(v::Vec3) = v[1]^2 + v[2]^2 + v[3]^2;
make_unit_vector(v::Vec3) = v/length(v);


text:
#################### regular (non-updatable) QR structure ######################
struct GivensQR{T, QT, RT<:AbstractMatrix{T}} <: AbstractQR{T}
    Q::QT
    R::RT
    n::Int
    m::Int
end

GivensQR(A::AbstractMatrix) = GivensQR!(copy(A))
# overwrites A
GivensQR!(A::AbstractMatrix) = GivensQR!(A, allocate_rotations(A))
# overwrites A and rotations
function GivensQR!(A::AbstractMatrix, rotations::AbstractVector{<:Givens})
    n, m = size(A)
    n ≥ m || throw(DimensionMismatch("GivensQR only supports systems with full column rank."))
    rot_index = 1
    for j in 1:m
        Aj = @view A[:, j:end] # do not need to treat columns 1:j-1 further
        @inbounds for i in n:-1:j+1
            G, rho = givens(Aj[i-1, 1], Aj[i, 1], i-1, i)
            lmul!(G, Aj)
            rotations[rot_index] = G
            rot_index += 1
        end
    end
    Q = GivensQ(rotations, n, m)
    R = UpperTriangular(@view(A[1:m, 1:m]))
    return GivensQR(Q, R, n, m)
end

Base.copy(F::GivensQR) = GivensQR(copy(F.Q), copy(F.R), F.n, F.m)
function Base.Matrix(F::GivensQR)
    n, m = size(F)
    A = zeros(eltype(F), n, m)
    @. A[1:m, 1:m] = F.R
    return lmul!(F.Q, A) # avoids every explicitly creating F.Q
end

# overwrites the first m elements of x with the solution
function LinearAlgebra.ldiv!(F::GivensQR, x::AbstractVecOrMat)
    ldiv!(F.Q, x)
    xm = typeof(x) <: AbstractVector ? view(x, 1:F.m) : view(x, 1:F.m, :)
    ldiv!(F.R, xm)
    return xm
end


text:
module MarginalTreatmentEffectsWithMultipleInstruments

    using JuMP
    using Clp
    using DataFrames
    using Parameters
    using LinearAlgebra
    using CSV
    using Suppressor
    using NamedArrays
    using DataStructures

    include("basis.jl")
    include("mtr.jl")
    include("dgp.jl")
    include("target_parameters.jl")
    include("ivlike.jl")
    include("mutual_consistency.jl")
    include("bounds.jl")
    include("paper_simulations.jl")
    include("util_save.jl")
    include("toplevel.jl")

end


text:
function ABCDQR_em(y,u,Ai,Bi,Ci,Di,Qi,Ri,m1i,P1i,max_iter::Int,tol::Float64,txo::Bool)
	#
	# estimate A, B, C, D, Q, R, m1, P1 using the EM algorithm for model
	#
	# x_{t+1} = A*x_{t} + B*u_{t} + w_{t}
	# y_{t}   = C*x_{t} + D*u_{t} + v_{t}
	#
	# cov(w_{t},v_{t}) = [Q 0;0 R]
	# x1 -> N(m1,P1)
	#
	# javier.cara@upm.es, 2016-02
	#

	# number of states
	nx = size(Ai,1)

	# data by rows
	y,ny,nt = byrow(y)
	u,nu,nt = byrow(u)

	# initial values
	A = Ai
	B = Bi
	C = Ci
	D = Di
	Q = Qi
	R = Ri
	m1 = m1i
	P1 = P1i

	# log-likelihood values
	loglikv = zeros(max_iter)

	# Syy, Suu and Syu do not depend on the iterations
	Syy = zeros(ny,ny)
	Suu = zeros(nu,nu)
	Syu = zeros(ny,nu)
	for t in 1:nt
		Syy = Syy + y[:,t]*y[:,t]'
		Suu = Suu + u[:,t]*u[:,t]'
		Syu = Syu + y[:,t]*u[:,t]'
	end

	tol1 = 1.0
	iter = 1
	while (iter <= max_iter) && (tol1 > tol)
		tic()

		# E-step
		# ---------------------------------------------------------------------------------
		# Kalmanfilter
		(xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik) = ABCDQR_kfilter(y,u,A,B,C,D,Q,R,m1,P1)
		(xtN,PtN,Pt1tN) = ACQR_ksmoother(A,xtt,Ptt,xtt1,Ptt1)

		loglikv[iter] = loglik
		if iter > 1
			tol1 = abs( (loglikv[iter] - loglikv[iter-1])/loglikv[iter-1] )
		end

		# initial values
		Sxx = zeros(nx,nx)
		Sx1x = zeros(nx,nx)
		Syx = zeros(ny,nx)
		Sx1u = zeros(nx,nu)
		Sxu = zeros(nx,nu)

		# matrices Sxx, Sx1x, Syx, Sx1u, Sxu, Sx1x1
		for t in 1:nt
			Sxx = Sxx + PtN[:,:,t] + xtN[:,t]*xtN[:,t]'
			Sx1x = Sx1x + Pt1tN[:,:,t] + xtN[:,t+1]*xtN[:,t]'
			Syx = Syx + y[:,t]*xtN[:,t]'
		  Sx1u= Sx1u + xtN[:,t+1]*u[:,t]'
		  Sxu = Sxu + xtN[:,t]*u[:,t]'
		end
		Sx1x1 = Sxx - (PtN[:,:,1] + xtN[:,1]*xtN[:,1]') + (PtN[:,:,nt+1] + xtN[:,nt+1]*xtN[:,nt+1]')

		# M-step
		# -------------------------------------------------------------------------------------
		# Matrices x0e y P0e
		m1 = xtN[:,1]
		P1 = PtN[:,:,1]

		# AB
		AB = [Sx1x Sx1u]/[Sxx Sxu;Sxu' Suu]

		# Matrix A
		A = AB[:,1:nx]

		# Matrix B
		B = AB[:,nx+1:nx+nu]

		# Matrix Q
		M1 = Sx1x*A'
		M2 = Sx1u*B'
		M3 = A*Sxu*B'
		Q = Sx1x1 - M1 - M1' - M2 - M2' + M3 + M3' + A*Sxx*A' + B*Suu*B'
		Q = 1/nt*Q
		Q = (Q + Q')/2 # to make sure it's a symmetric matrix

		# CD
		CD = [Syx Syu]/[Sxx Sxu;Sxu' Suu]

		# Matrix Ae
		C = CD[:,1:nx]

		# Matrix Be
		D = CD[:,nx+1:nx+nu]

		# Matrix Re
		M1 = Syx*C'
		M2 = Syu*D'
		M3 = C*Sxu*D'
		R = Syy - M1' - M1 - M2 - M2' + M3 + M3' + C*Sxx*C' + D*Suu*D'
		R = 1/nt*R
		R = (R + R')/2 # to make sure it's a symmetric matrix

		etime = toq()
		if txo
			println( "Iter " * @sprintf("%3d",iter) * ",   @time = " * @sprintf("%5.2f",etime) * ",   logLik = " * @sprintf("%.6E",loglik) *
			         ",   tol = ", @sprintf("%.2E",tol1) )
		end

		iter += 1

	end

	loglikv = loglikv[1:(iter-1)]

	# Akaike Information Criterion
	P = ny*nx + nx*nx + nx*(nx+1)/2 + ny*(ny+1)/2
	aic = -2*loglikv[end] + 2*P

	# output
	return A,B,C,D,Q,R,m1,P1,loglikv,aic

end


text:
module OrdinaryDiffEq

  using Reexport
  @reexport using DiffEqBase

  using Logging

  using MuladdMacro, SparseArrays, FastClosures

  using LinearAlgebra

  # Interfaces
  import DiffEqBase: solve!, step!, initialize!, isadaptive

  # Internal utils
  import DiffEqBase: ODE_DEFAULT_NORM, ODE_DEFAULT_ISOUTOFDOMAIN, ODE_DEFAULT_PROG_MESSAGE, ODE_DEFAULT_UNSTABLE_CHECK

  using DiffEqBase: DiffEqArrayOperator, DEFAULT_UPDATE_FUNC

  using DiffEqBase: TimeGradientWrapper, UJacobianWrapper, TimeDerivativeWrapper, UDerivativeWrapper

  import RecursiveArrayTools: chain, recursivecopy!

  using UnPack, GenericSVD, ForwardDiff, RecursiveArrayTools,
        DataStructures, FiniteDiff, ArrayInterface

  import ForwardDiff.Dual

  using ExponentialUtilities

  using NLsolve
  # Required by temporary fix in not in-place methods with 12+ broadcasts
  # `MVector` is used by Nordsieck forms
  import StaticArrays: SArray, MVector, SVector, @SVector, StaticArray

  # Integrator Interface
  import DiffEqBase: resize!,deleteat!,addat!,full_cache,user_cache,u_cache,du_cache,
                     resize_non_user_cache!,deleteat_non_user_cache!,addat_non_user_cache!,
                     terminate!,get_du, get_dt,get_proposed_dt,set_proposed_dt!,
                     u_modified!,savevalues!,add_tstop!,add_saveat!,set_reltol!,
                     set_abstol!, postamble!, last_step_failed,
                     isautodifferentiable

  using DiffEqBase: check_error!, @def, @.. , _vec, _reshape

  using DiffEqBase: AbstractNLSolverAlgorithm, AbstractNLSolverCache, NLStatus
  using DiffEqBase: nlsolve_f, qrdelete!, qradd!, build_jac_config, resize_jac_config!
  using DiffEqBase: Convergence, Divergence
  const TryAgain = DiffEqBase.SlowConvergence

  import DiffEqBase: calculate_residuals, calculate_residuals!, unwrap_cache, @tight_loop_macros, islinear, timedepentdtmin

  import SparseDiffTools
  import SparseDiffTools: matrix_colors, forwarddiff_color_jacobian!,
                          forwarddiff_color_jacobian, ForwardColorJacCache,
                          default_chunk_size, getsize

  using MacroTools, Adapt

  const CompiledFloats = Union{Float32,Float64,
    ForwardDiff.Dual{ForwardDiff.Tag{T,W},K,3} where {T,W<:Union{Float64,Float32},
                                                        K<:Union{Float64,Float32}}}

  include("misc_utils.jl")
  include("algorithms.jl")
  include("alg_utils.jl")

  include("nlsolve/type.jl")
  include("nlsolve/utils.jl")
  include("nlsolve/nlsolve.jl")
  include("nlsolve/functional.jl")
  include("nlsolve/newton.jl")

  include("generic_rosenbrock.jl")

  include("caches/basic_caches.jl")
  include("caches/low_order_rk_caches.jl")
  include("caches/high_order_rk_caches.jl")
  include("caches/low_storage_rk_caches.jl")
  include("caches/ssprk_caches.jl")
  include("caches/feagin_caches.jl")
  include("caches/verner_caches.jl")
  include("caches/sdirk_caches.jl")
  include("caches/firk_caches.jl")
  include("caches/kencarp_kvaerno_caches.jl")
  include("caches/linear_caches.jl")
  include("caches/linear_nonlinear_caches.jl")
  include("caches/symplectic_caches.jl")
  include("caches/rosenbrock_caches.jl")
  include("caches/rkn_caches.jl")
  include("caches/adams_bashforth_moulton_caches.jl")
  include("caches/nordsieck_caches.jl")
  include("caches/bdf_caches.jl")
  include("caches/rkc_caches.jl")
  include("caches/extrapolation_caches.jl")
  include("caches/prk_caches.jl")
  include("caches/pdirk_caches.jl")
  include("caches/dae_caches.jl")

  include("tableaus/low_order_rk_tableaus.jl")
  include("tableaus/high_order_rk_tableaus.jl")
  include("tableaus/symplectic_tableaus.jl")
  include("tableaus/verner_tableaus.jl")
  include("tableaus/feagin_tableaus.jl")
  include("tableaus/rosenbrock_tableaus.jl")
  include("tableaus/sdirk_tableaus.jl")
  include("tableaus/firk_tableaus.jl")
  include("tableaus/rkn_tableaus.jl")
  include("tableaus/rkc_tableaus.jl")

  include("integrators/type.jl")
  include("integrators/controllers.jl")
  include("integrators/integrator_utils.jl")
  include("cache_utils.jl")
  include("integrators/integrator_interface.jl")
  include("initialize_dae.jl")
  include("wrappers.jl")

  include("perform_step/fixed_timestep_perform_step.jl")
  include("perform_step/symplectic_perform_step.jl")
  include("perform_step/rkn_perform_step.jl")
  include("perform_step/split_perform_step.jl")
  include("perform_step/linear_perform_step.jl")
  include("perform_step/exponential_rk_perform_step.jl")
  include("perform_step/explicit_rk_perform_step.jl")
  include("perform_step/low_order_rk_perform_step.jl")
  include("perform_step/high_order_rk_perform_step.jl")
  include("perform_step/verner_rk_perform_step.jl")
  include("perform_step/feagin_rk_perform_step.jl")
  include("perform_step/low_storage_rk_perform_step.jl")
  include("perform_step/ssprk_perform_step.jl")
  include("perform_step/sdirk_perform_step.jl")
  include("perform_step/kencarp_kvaerno_perform_step.jl")
  include("perform_step/firk_perform_step.jl")
  include("perform_step/rosenbrock_perform_step.jl")
  include("perform_step/composite_perform_step.jl")
  include("perform_step/adams_bashforth_moulton_perform_step.jl")
  include("perform_step/nordsieck_perform_step.jl")
  include("perform_step/bdf_perform_step.jl")
  include("perform_step/rkc_perform_step.jl")
  include("perform_step/extrapolation_perform_step.jl")
  include("perform_step/prk_perform_step.jl")
  include("perform_step/pdirk_perform_step.jl")
  include("perform_step/dae_perform_step.jl")

  include("dense/generic_dense.jl")
  include("dense/interpolants.jl")
  include("dense/rosenbrock_interpolants.jl")
  include("dense/stiff_addsteps.jl")
  include("dense/low_order_rk_addsteps.jl")
  include("dense/verner_addsteps.jl")
  include("dense/high_order_rk_addsteps.jl")

  include("derivative_utils.jl")
  include("nordsieck_utils.jl")
  include("adams_utils.jl")
  include("bdf_utils.jl")
  include("rkc_utils.jl")
  include("derivative_wrappers.jl")
  include("iterator_interface.jl")
  include("constants.jl")
  include("composite_solution.jl")
  include("solve.jl")
  include("initdt.jl")
  include("interp_func.jl")
  include("composite_algs.jl")

  #General Functions
  export solve, solve!, init, step!

  export OrdinaryDiffEqAlgorithm

  #Callback Necessary
  export addsteps!, ode_interpolant, terminate!, savevalues!, copyat_or_push!, isfsal

  export constructDormandPrince

  # Reexport the Alg Types

  export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3, OwrenZen4, OwrenZen5,
         BS3, BS5, RK46NL, DP5, Tsit5, DP8, Vern6, Vern7, Vern8, TanYam7, TsitPap8,
         Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65, PFRK87, RKM

  export SSPRK22, SSPRK33, KYKSSPRK42, SSPRK53, SSPRK53_2N1, SSPRK53_2N2, SSPRK53_H, SSPRK63, SSPRK73, SSPRK83, SSPRK43, SSPRK432,
         SSPRKMSVS32, SSPRKMSVS43, SSPRK932, SSPRK54, SSPRK104

  export ORK256, CarpenterKennedy2N54, SHLDDRK64, HSLDDRK64, DGLDDRK73_C, DGLDDRK84_C, DGLDDRK84_F, NDBLSRK124, NDBLSRK134, NDBLSRK144,
         CFRLDDRK64, TSLDDRK74,SHLDDRK52,SHLDDRK_2N,
         CKLLSRK43_2,CKLLSRK54_3C,CKLLSRK95_4S,CKLLSRK95_4C,CKLLSRK95_4M,
         CKLLSRK54_3C_3R, CKLLSRK54_3M_3R, CKLLSRK54_3N_3R, CKLLSRK85_4C_3R, CKLLSRK85_4M_3R, CKLLSRK85_4P_3R,
         CKLLSRK54_3N_4R, CKLLSRK54_3M_4R, CKLLSRK65_4M_4R, CKLLSRK85_4FM_4R, CKLLSRK75_4M_5R,
         ParsaniKetchesonDeconinck3S32, ParsaniKetchesonDeconinck3S82,
         ParsaniKetchesonDeconinck3S53, ParsaniKetchesonDeconinck3S173,
         ParsaniKetchesonDeconinck3S94, ParsaniKetchesonDeconinck3S184,
         ParsaniKetchesonDeconinck3S105, ParsaniKetchesonDeconinck3S205,
         KYK2014DGSSPRK_3S2

  export RadauIIA3, RadauIIA5

  export ImplicitEuler, ImplicitMidpoint, Trapezoid, TRBDF2, SDIRK2, SDIRK22,
         Kvaerno3, KenCarp3, Cash4, Hairer4, Hairer42, SSPSDIRK2, Kvaerno4,
         Kvaerno5, KenCarp4, KenCarp47, KenCarp5, KenCarp58, ESDIRK54I8L2SA, SFSDIRK4, SFSDIRK5, CFNLIRK3, SFSDIRK6, SFSDIRK7, SFSDIRK8
         Kvaerno5, KenCarp4, KenCarp5, ESDIRK54I8L2SA, SFSDIRK4, SFSDIRK5, CFNLIRK3, SFSDIRK6, SFSDIRK7, SFSDIRK8

  export MagnusMidpoint, LinearExponential, MagnusLeapfrog, LieEuler, CayleyEuler, MagnusGauss4, MagnusNC6, MagnusGL6, MagnusGL8, MagnusNC8, MagnusGL4,
         MagnusAdapt4, RKMK2, RKMK4, LieRK4, CG2, CG3

  export Rosenbrock23, Rosenbrock32, RosShamp4, Veldd4, Velds4, GRK4T, GRK4A,
         Ros4LStab, ROS3P, Rodas3, Rodas4, Rodas42, Rodas4P, Rodas5,
         RosenbrockW6S4OS, ROS34PW1a, ROS34PW1b, ROS34PW2, ROS34PW3

  export LawsonEuler, NorsettEuler, ETD1, ETDRK2, ETDRK3, ETDRK4, HochOst4, Exp4, EPIRK4s3A, EPIRK4s3B,
         EPIRK5s3, EXPRB53s3, EPIRK5P1, EPIRK5P2, ETD2, Exprb32, Exprb43

  export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog,
         McAte2, Ruth3, McAte3, CandyRoz4, McAte4, McAte42, McAte5,
         CalvoSanz4, Yoshida6, KahanLi6, McAte8, KahanLi8, SofSpa10

  export SplitEuler

  export Nystrom4, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent,
         IRKN3, IRKN4, DPRKN6, DPRKN8, DPRKN12, ERKN4, ERKN5

  export ROCK2, ROCK4, RKC, IRKC, ESERK4, ESERK5, SERK2

  export AB3, AB4, AB5, ABM32, ABM43, ABM54

  export VCAB3, VCAB4, VCAB5, VCABM3, VCABM4, VCABM5

  export VCABM

  export IMEXEuler, CNAB2, CNLF2

  export AN5, JVODE, JVODE_Adams, JVODE_BDF

  export ABDF2, QNDF1, QBDF1, QNDF2, QBDF2, QNDF, QBDF

  export SBDF2, SBDF3, SBDF4

  export MEBDF2

  export AutoSwitch, AutoTsit5, AutoDP5,
         AutoVern6, AutoVern7, AutoVern8, AutoVern9

  export AitkenNeville, ExtrapolationMidpointDeuflhard, ExtrapolationMidpointHairerWanner, ImplicitEulerExtrapolation,
         ImplicitDeuflhardExtrapolation, ImplicitHairerWannerExtrapolation, ImplicitEulerBarycentricExtrapolation

  export KuttaPRK2p5, PDIRK44, DImplicitEuler, DABDF2

  export ShampineCollocationInit, BrownFullBasicInit, NoInit
end # module


text:
using Revise
using Soss

m = @model begin
    n = 100 # number of samples
    m = 100 # reviewers for each example

    p_bad ~ Beta(1,3) |> iid(n)

    fpr ~ Beta(2,5) |> iid(m)
    tpr ~ Beta(5,2) |> iid(m)

    y ~ For(n, m) do nj,mj
            Mix([Bernoulli(fpr[mj]), Bernoulli(tpr[mj])]
                , [1 - p_bad[nj], p_bad[nj]])
        end
end;

using Random
truth = rand(m());
symlogpdf(m())

using BenchmarkTools
@btime logpdf(m(),truth)
@btime logpdf(m(),truth, codegen)

f1 = Soss._codegen(m, true);
f2 = Soss._codegen(m,false);

@btime f1((),truth)
@btime f2((),truth)

codegen(m(),truth)

logpdf(m(), merge(truth, (p_bad=shuffle(truth.p_bad),)), codegen)

@time result = dynamicHMC(m(), (y=truth.y,), codegen) ;

# result = @time advancedHMC(m(), (y=truth.y,))

pairs(truth)
result |> particles |> pairs


text:
module psi4

using PyCall
using ..Calculators
using ..Clusters

export Psi4

"""
    gpsi4init=false

Global variable to see if Psi4 was initiated for current process
"""
global gpsi4init=false

"""
    gPsi4 = undef

Hold Psi4 object given by PyCall. Used to do calculation in current process
"""
global gPsi4 = undef

"""
    initpsi4(;memory="500 MiB", quiet=true, nthreads=1)

Used to intialize Psi4 environment
"""
function initpsi4(;memory="500 MiB", quiet=true, nthreads=1)
    global gPsi4 = pyimport("psi4")
    global gpsi4init = true
    quiet && gPsi4.core.be_quiet()
    gPsi4.set_memory(memory)
    nthreads > 1 && gPsi4.set_num_threads(nthreads)
end

"""
    mutable struct Psi4 <: AbstractCalculationProgram

Holds information that calculations are to be done with Psi4

# Fields
- `memory="500MiB"` : memory used by Psi4
- `nthreads=1`  : number of threads used by Psi4
"""
mutable struct Psi4 <: AbstractCalculationProgram
    memory::String
    nthreads::UInt
    function Psi4(;memory="500MiB", nthreads=1)
        initpsi4(memory=memory, nthreads=nthreads)
        @debug gPsi4
        new(memory, nthreads)
    end
end

function Calculators.calculate_energy(cal::Calculator{Psi4}, point::Cluster;
                                      basename="base", ghost=undef, id="", pchannel=undef)
    ! gpsi4init && initpsi4(memory=cal.calculator.memory, nthreads=cal.calculator.nthreads)
    s=sprint( (io, x) -> print_xyz(io,x, printheader=false), point)
    c = gPsi4.geometry(s)
    out = gPsi4.energy(cal.method*"/"*cal.basis, molecule=c)
    pchannel != undef && put!(pchannel,true)
    return out
end

function Calculators.calculate_energy(cal::Calculator{Psi4}, points;
                                      basename="base", ghost=undef, id="", pchannel=undef)
    return map( x -> calculate_energy(cal, x, basename=basename, ghost=ghost, id=id, pchannel=pchannel), points )
end

function Calculators.bsse_corrected_energy(cal::Calculator{Psi4}, c1::Cluster, c2::Cluster;
                                           basename="base", id="", pchannel=undef)
    ! gpsi4init && initpsi4(memory=cal.calculator.memory)
    s1=sprint( (io, x) -> print_xyz(io,x, printheader=false), c1)
    s2=sprint( (io, x) -> print_xyz(io,x, printheader=false), c2)
    c = gPsi4.geometry(s1*"\n--\n"*s2)
    out = gPsi4.energy(cal.method*"/"*cal.basis, molecule=c, bsse_type="cp")
    pchannel != undef && put!(pchannel,true)
    return out
end

function Calculators.bsse_corrected_energy(cal::Calculator{Psi4}, c1, c2;
                                           basename="base", id="", pchannel=undef)
    return map( (x,y) -> bsse_corrected_energy(cal, x, y, basename=basename, id=id, pchannel=pchannel ), c1, c2  )
end

end  # module psi4


text:
# Write function to automate plotting qqplot
# - Plots standard quantile-quantile plot for assessing linear model fits

function plot_qq(model)

    # Extract predicted values from model object
    pred = GLM.predict(model)

    # Extract residuals from model object
    res = GLM.residuals(model)

    # Calculate mean of residuals
    res_mean = StatsBase.mean(res)

    # Calculate length (n) of residuals vector
    n = length(res)

    # Calculate standardised residuals
    stdRes = (res - res_mean * ones(n)) / (Statistics.std(res))

    # Define quantiles
    qx = Distributions.quantile.(Distributions.Normal(),
                                    range(0.5,
                                    stop = (n - 0.5),
                                    length = (n)) / (n + 1))

    # Create plot
    p = Gadfly.plot(
                    # Add points layer
                    Gadfly.layer(x = qx,
                                 y = sort(stdRes),
                                 Gadfly.Geom.point),
                    # Add 1:1 line
                    Gadfly.layer(x = [-3,3],
                                 y = [-3,3],
                                 Gadfly.Geom.line,
                                 Gadfly.style(line_style = [:dot])),
                    # Change plot aesthetics
                    Gadfly.Guide.title("Normal Q-Q plot"),
                    Gadfly.Guide.xlabel("Theoretical Quantiles"),
                    Gadfly.Guide.ylabel("Standardized residuals"))

end

max_stars_repo_name: guysutton/GLMdiagnostics

max_issues_repo_name: guysutton/GLMdiagnostics

max_forks_repo_name: guysutton/GLMdiagnostics


text:
export keeprow_vec, filter_tbl_by_splits

using DataFrames: nrow
using Tables

"""
    Returns the boolean that can used to filter the table

* tbl - A Tables.jl table
* node - A JLBoostTree
"""
function keeprow_vec(tbl, node::Union{Nothing, AbstractJLBoostTree})::BitArray
    @assert Tables.istable(tbl)

    tblc = Tables.columns(tbl)

    # a boolean on whether to keep the row
    if node === nothing
        return keeprow = trues(nrow(tbl))
    end

    if node.parent.parent === nothing
        keeprow = trues(nrow(tbl))
    else
        keeprow = keeprow_vec(tbl, node.parent)
    end

    if is_left_child(node)
        keeprow .&= getproperty(tblc, node.parent.splitfeature) .<= node.parent.split
    else
        keeprow .&= getproperty(tblc, node.parent.splitfeature) .> node.parent.split
    end

    keeprow
end

"""
    Filter the Tables.jl compitable `tbl` by the logic embedded in the node

* tbl - A Tables.jl table
* node - A JLBoostTree
"""
function filter_tbl_by_splits(tbl, node::AbstractJLBoostTree)
    @assert Tables.istable(tbl)

    tblc = Tables.columns(tbl)

    # a boolean on whether to keep the row
    keeprow = keeprow_vec(tbl, node)

    view(tblc, keeprow, :)
end


text:
function PV_savings(consumer_input::Consumer, system::PVSystem; tariff_increase=true)

    consumer = deepcopy(consumer_input)

    netsavings = Array{Float64,2}(undef,12,consumer.decision_timeframe)

    for y in 1:consumer.decision_timeframe

        ebalance = annual_energy_balance(consumer, system, print_output=false);

        for m in 1:12

            bill = monthly_bill(ebalance[m], consumer, print_output=false)

            netsavings[m,y] = bill["total_savings"]

        end

        tariff_increase ? increase_tariff!(consumer) : true

    end

    return netsavings

end

function PV_cashflow(consumer::T, system::PVSystem, financial_terms::Financial; tariff_increase=true) where {T <: Consumer}

    #println(consumer.econsumption, " ", system.capacity)
    if system.capacity > 0.0
        PV = PVCost(system,1.0)
        payment = amortize(financial_terms.apr/12, financial_terms.term*12, PV*(1-financial_terms.downpayment))
    else
        #dirty way to have 0 payments, needs some fixing later
        PV = 0.0
		payment = amortize(financial_terms.apr/12, financial_terms.term*12, 0.0)
    end

    bill_savings = PV_savings(consumer, system; tariff_increase=tariff_increase)

    for y in 1:financial_terms.term

        for m in 1:12

            bill_savings[m,y] = bill_savings[m,y] - payment[4]*financial_terms.XhR

        end

    end

	bill_savings[1,1] -= PV*(financial_terms.downpayment)

    return bill_savings

end

function PV_cashflow(consumer::Residential, system::PVSystem, financial_terms::Financial; tariff_increase=true)

    #println(consumer.econsumption, " ", system.capacity)
    if system.capacity > 0.0
        PV = PVCost(system,financial_terms.XhR)
        payment = amortize(financial_terms.apr/12, financial_terms.term*12, PV*(1-financial_terms.downpayment))
    else
        #dirty way to have 0 payments, needs some fixing later
        payment = amortize(financial_terms.apr/12, financial_terms.term*12, 0.0)
    end

    bill_savings = PV_savings(consumer, system; tariff_increase=tariff_increase)

    for y in 1:financial_terms.term

        for m in 1:12

            bill_savings[m,y] = bill_savings[m,y] - payment[4]

        end

    end

    return bill_savings

end

function optimal_pv(consumer::Residential, system::PVSystem, capacity_range::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, financial_terms::Financial; tariff_increase=true)

    results = Array{Float64,2}(undef, length(capacity_range), 2)

    for (ix,cap) in enumerate(capacity_range)
            system.capacity = cap
            cash_flows = PV_cashflow(consumer, system, financial_terms, tariff_increase=tariff_increase)
            annuities = sum(cash_flows,dims=1)

			if sum(annuities[1:5]) < 100

				results[ix,2] = -99
				results[ix,1] = 0.0;

			else

				results[ix,2] = npv(annuities, consumer.rate_return)
				results[ix,1] = cap;

            end

    end

	fm = findmax(results[:,2])

	if fm[1] == -99
		return (0.0, 0.0), results
	else
		return (fm[1], capacity_range[fm[2]]), results
	end

end

function optimal_pv(consumer::T, system::PVSystem, capacity_range::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, financial_terms::Financial; tariff_increase=true) where {T <: Consumer}

    results = Array{Float64,2}(undef, length(capacity_range), 2)

    for (ix,cap) in enumerate(capacity_range)
            system.capacity = cap
            cash_flows = PV_cashflow(consumer, system, financial_terms, tariff_increase=tariff_increase)
            annuities = sum(cash_flows,dims=1)
            if sum(annuities) < 100

				results[ix,2] = -99
				results[ix,1] = 0.0;

			else

				results[ix,2] = npv(annuities, consumer.rate_return)
				results[ix,1] = cap;

            end

    end

	fm = findmax(results[:,2])

	if fm[1] == -99
		return (0.0, 0.0), results
	else
		return (fm[1], capacity_range[fm[2]]), results
	end

end

function plot_monte_carlo_model_prediction(company, tariff, random_consumption, mc_data; width=50)
    fig = plt.figure()
    ax1 = fig.add_subplot()
    ax1.boxplot([[mc_data[j][i] for j in 1:size(mc_data,1)] for i in 1:size(mc_data[1],1)], positions=random_consumption, manage_ticks=false, widths=width)
#     ax1.set_xticks(random_consumption)
#     ax1.set_xticklabels(random_consumption, rotation="vertical")
    ylabel("Optimal PV System Capacity [kW]")
    xlabel("Consumer Monthly Energy use [kWh]")
    grid("on");
    title(string("Optimal System Choice for ", company, " ", tariff, " Consumer (Various Locations)"))
end

# NOTE: be sure to set the appropriate tariff before calling this function
function retrieve_monte_carlo_model_prediction(company, num_samples, tariff, tariff_obj, random_consumption, cap_range)
    mc_filename = string("data/monte_carlo_data/", company, "_", num_samples, ".txt")
    if isfile(mc_filename)
        mc_pv_output = readdlm(mc_filename, '\t', Float64, '\n')
    else
        println("Can't find NSRDB+SAM output")
        return
    end

    mc_data = []

    for pv_output in eachcol(mc_pv_output)
        mcPV = -1 # Placeholder because Julia's "Nothing" is strange and I don't understand it yet
        if tariff == "R"
            mcPV = newPVRes(pv_output)
        elseif tariff == "CI"
            mcPV = newPVComInd(pv_output)
        else
            mcPV = newPVTMT(pv_output)
        end
        pv_data = []
        for (ix, co) in enumerate(random_consumption)
            tariff_obj.econsumption = co; get_pmax(tariff_obj);
            new_data = optimal_pv(tariff_obj, mcPV, cap_range, BAC1, tariff_increase = true)
            push!(pv_data, new_data[1][2])
        end
        push!(mc_data, pv_data)
    end
    return mc_data
end

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
const Point{T} = AbstractVector{T} where T<:Real
const Points = AbstractVector{<:Point}
const Sigma = AbstractVector{<:Integer}  # Encoded the vertex by the ids of its generators.
const Vertex = Tuple{<:Sigma, <: Point}
const Vertices = Dict{<:Sigma, <:Point}

voronoi(x; kwargs...) = voronoi(vecvec(x); kwargs...)

""" construct the voronoi diagram from `x` through breadth-first search """
function voronoi(xs::Points, searcher = Raycast(xs))
    sig, r = descent(xs, searcher)
    verts, boundary = explore(sig, r, xs, searcher)
    return verts::Vertices, xs, boundary
end

voronoi_random(x, args...; kwargs...) = voronoi_random(vecvec(x), args...; kwargs...)

""" construct a (partial) voronoi diagram from `x` through a random walk """
function voronoi_random(xs::Points, iter=1000; maxstuck=typemax(Int))
    searcher = Raycast(xs)
    sig, r = descent(xs, searcher)
    verts = walk(sig, r, iter, xs, searcher, maxstuck)
    return verts::Vertices, xs
end

vecvec(x::Matrix) = map(SVector{size(x,1)}, eachcol(x))
vecvec(x::Vector{<:SVector}) = x

""" starting at given points, run the ray shooting descent to find vertices """
function descent(xs::Points, searcher, start = 1) :: Vertex
    sig = [start]
    r = xs[start]
    d = length(r)

    for k in d:-1:1  # find an additional generator for each dimension
        u = randray(xs[sig])
        (tau, t) = raycast(sig, r, u, xs, searcher)
        if t == Inf
            u = -u
            (tau, t) = raycast(sig, r, u, xs, searcher)
        end
        if t == Inf
            error("Could not find a vertex in both directions of current point." *
                "Consider increasing search range (tmax)")
        end
        sig = tau
        r = r + t*u
    end

    return (sig, r)
end

""" starting at vertices, walk nsteps along the voronoi graph to find new vertices """
function walk(sig::Sigma, r::Point, nsteps::Int, xs::Points, searcher, maxstuck::Int) :: Vertices
    verts = Dict(sig => r)
    nonew = 0
    prog = Progress(maxstuck, 1, "Voronoi walk")
    progmax = 0

    for s in 1:nsteps
        nonew += 1
        progmax = max(progmax, nonew)
        ProgressMeter.update!(prog, progmax)
        randdir = rand(1:length(sig))
        sig, r = walkray(sig, r, xs, searcher, randdir)
        r == Inf && continue  # we encountered an infinite ray
        get!(verts, sig) do
            nonew = 0
            return r
        end
        nonew > maxstuck && break
    end

    return verts
end

""" find the vertex connected to `v` by moving away from its `i`-th generator """
function walkray(sig::Sigma, r::Point, xs::Points, searcher, i)
    sig_del = deleteat(sig, i)
    u = randray(xs[sig_del])
    if (u' * (xs[sig[i]] - xs[sig_del[1]])) > 0
        u = -u
    end
    sig′, t = raycast(sig_del, r, u, xs, searcher)
    if t < Inf
        r′ = r + t*u
        return sig′, r′
    else
        return sig, r  # if the vertex has an unbounded ray, return the same vertex
    end
end

""" BFS of vertices starting from `S0` """
function explore(sig, r, xs::Points, searcher) # :: Vertices
    verts = Dict(sig=>r)
    queue = copy(verts)
    edgecount = Dict{Vector{Int64}, Int}()
    boundary = []

    cache = true  # Caches if both points along an edge are known. Trades memory for runtime.
    hit = 0
    miss = 0
    new = 0
    unbounded = 0

    while length(queue) > 0
        (sig,r) = pop!(queue)
        for i in 1:length(sig)

            if cache && get(edgecount, deleteat(sig, i), 0) == 2
                hit += 1
                continue
            end

            sig′, r′ = walkray(sig, r, xs, searcher, i)

            if sig′ == sig
                push!(boundary, sig)
                unbounded += 1
                continue
            end

            if !haskey(verts, sig′)
                push!(queue, sig′ => r′)
                push!(verts, sig′ => r′)
                if cache
                    for j in 1:length(sig′)
                        edge = deleteat(sig′, j)
                        edgecount[edge] = get(edgecount, edge, 0) + 1
                    end
                end
                new += 1
            else
                miss += 1
            end
        end

    end

    #@show hit, miss, new, unbounded
    return verts, boundary
end

deleteat(sig, i) = deleteat!(copy(sig), i)

""" generate a random ray orthogonal to the subspace spanned by the given points """
function randray(xs::Points)
    k = length(xs)
    d = length(xs[1])
    v = similar(xs, k-1)

    # Gram Schmidt
    for i in 1:k-1
        v[i] = xs[i] .- xs[k]
        for j in 1:(i-1)
            v[i] = v[i] .- dot(v[i], v[j]) .* v[j]
        end
        v[i] = normalize(v[i])
    end

    u = randn(d)
    for i in 1:k-1
        u = u - dot(u, v[i]) * v[i]
    end
    u = normalize(u)
    return u
end


text:
# cd(@__DIR__)

# using Pkg
# Pkg.activate(@__DIR__)
# Pkg.instantiate()

using Documenter
using GeometricalOptics, Plots

makedocs(
    sitename = "GeometricalOptics.jl",
    format = Documenter.HTML(),
    modules = [GeometricalOptics],
    pages=[ "index.md",
            "Quickstart.md",
            "Options.md",
            "Examples.md",
            "TroubleConv.md",
            "Reference.md"]
)

# 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/airspaced-nk5/GeometricalOptics.jl.git",
       )


text:
using OmegaCore, Distributions
using Test
import Random
Random.seed!(1)

function test_cond!()
  function f(ω)
    x = (1 ~ Normal(0, 1))(ω)
    cond!(ω, x > 0)
    x
  end
  @test all(randsample(f, 100) .>= 0.0)
end

"Test conditioning on predicate of positive measure"
function test_pos_measure()
  x = 1 ~ Normal(0, 1)
  y(ω) = x(ω) > 0
  x_cond = x |ᶜ y
  samples = randsample(x_cond, 100000; alg = RejectionSample)
  samplemean = mean(samples)
  exactmean = mean(truncated(Normal(0, 1), 0, Inf))
  @test samplemean ≈ exactmean atol = 0.01
end

"Test whether the logpdf is correct"
function test_density_cond()
  rng = Random.MersenneTwister(0)
  μ = 1 ~ Normal(0, 1)
  x = Normal(μ, 1.0)

  μ_ = -0.4321
  x_ = 0.1234
  μₓ = μ |ᶜ (x ==ₚ x_)
  ω = LazyΩ(μ => μ_)
  logpdf_ = logpdf(Normal(0, 1), μ_) + logpdf(Normal(μ_, 1), x_)
  @test logpdf(μₓ, ω) == logpdf_
end

function test_out_of_order_condition()
  x = 1 ~ Poisson(1.3)
  function f(ω)
    n = x(ω)
    x_ = 0.0
    for i = 1:n
      x_ += (i + 1 ~ Normal(0, 1))(ω)
    end
    x_
  end

  x_ = x |ᶜ (f ==ₚ 3.0)
  g_ = ω -> (x(ω), f_(ω))
end

function test_parent()
  μ = 1 ~ Normal(0, 1)
  x = 2 ~ Normal(μ, 1)
  x_ = 0.123
  μ_ = 0.1
  μ |ᶜ x == x_
  ω = defΩ()((1,) => μ_)
  @test logpdf(μ_, ω) = logpdf(Normal(0, 1), μ_) + logpdf(Normal(μ_, 1), x_)
end

function test_condition()
  μ = 1 ~ Normal(0, 1)
  x = 2 ~ Normal(μ, 1)
  x_ = 0.123
  μ_ = 0.987

  μc = μ |ᶜ x ==ₚ x_
  ω = defω()
  ω[(1,)] = μ_
  μc(ω)
  @test x(ω) == x_
  logpdf_ = logpdf(Normal(0, 1), μ) + logpdf(Normal(μ, 1), x_)
  # Will ω have a value for (2,)?
end

@testset "Conditions" begin
  test_cond!()
  test_pos_measure()
  # test_density_cond()
  # test_out_of_order_condition()
  # test_parent()
  # test_condition()
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
#*******************************************************/
#* Copyright(c) 2018 by Artelys                        */
#* This source code is subject to the terms of the     */
#* MIT Expat License (see LICENSE.md)                  */
#*******************************************************/

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#  This example demonstrates how to use Knitro to solve the following
#  simple quadratically constrained quadratic programming problem(QCQP).
#
#  min   1000 - x0^2 - 2 x1^2 - x2^2 - x0 x1 - x0 x2 #  s.t.  8 x0 + 14 x1 + 7 x2 = 56
#        x0^2 + x1^2 + x2^2 >= 25
#        x0 >= 0, x1 >= 0, x2 >= 0
#
#  The start point(2, 2, 2) converges to the minimum at(0, 0, 8),
#  with final objective = 936.0.  From a different start point,
#  Knitro may converge to an alternate local solution at(7, 0, 0),
#  with objective = 951.0.
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using KNITRO
using Test

function example_qcqp1(; verbose=true)
    # Create a new Knitro solver instance.
    kc = KNITRO.KN_new()

    # Illustrate how to override default options by reading from
    # the knitro.opt file.
    options = joinpath(dirname(@__FILE__), "..", "examples", "knitro.opt")
    KNITRO.KN_load_param_file(kc, options)
    kn_outlev = verbose ? KNITRO.KN_OUTLEV_ITER : KNITRO.KN_OUTLEV_NONE
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_OUTLEV, kn_outlev)

    # Initialize Knitro with the problem definition.

    # Add the variables and set their bounds and initial values.
    # Note: unset bounds assumed to be infinite.
    KNITRO.KN_add_vars(kc, 3)
    KNITRO.KN_set_var_lobnds_all(kc, [0.0, 0.0, 0.0])
    KNITRO.KN_set_var_primal_init_values_all(kc, [2.0, 2.0, 2.0])

    # Add the constraints and set their bounds.
    KNITRO.KN_add_cons(kc, 2)
    KNITRO.KN_set_con_eqbnd(kc, 0, 56.0)
    KNITRO.KN_set_con_lobnd(kc, 1, 25.0)

    # Add coefficients for linear constraint.
    lconIndexVars = Int32[0, 1, 2]
    lconCoefs = [8.0, 14.0, 7.0]
    KNITRO.KN_add_con_linear_struct_one(kc, 3, 0, lconIndexVars, lconCoefs)

    # Add coefficients for quadratic constraint
    qconIndexVars1 = Int32[0, 1, 2]
    qconIndexVars2 = Int32[0, 1, 2]
    qconCoefs = [1.0, 1.0, 1.0]
    KNITRO.KN_add_con_quadratic_struct_one(
        kc,
        3,
        1,
        qconIndexVars1,
        qconIndexVars2,
        qconCoefs,
    )

    # Set minimize or maximize(if not set, assumed minimize)
    KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MINIMIZE)

    # Add constant value to the objective.
    KNITRO.KN_add_obj_constant(kc, 1000.0)

    # Set quadratic objective structure.
    qobjIndexVars1 = Int32[0, 1, 2, 0, 0]
    qobjIndexVars2 = Int32[0, 1, 2, 1, 2]
    qobjCoefs = [-1.0, -2.0, -1.0, -1.0, -1.0]

    KNITRO.KN_add_obj_quadratic_struct(kc, 5, qobjIndexVars1, qobjIndexVars2, qobjCoefs)

    # Solve the problem.
    #
    # Return status codes are defined in "knitro.h" and described
    # in the Knitro manual.
    nStatus = KNITRO.KN_solve(kc)
    nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)

    # An example of obtaining solution information.
    if verbose
        println("Knitro converged with final status = ", nStatus)
        println("  optimal objective value  = ", objSol)
        println("  optimal primal values x  = ", x)
        println("  feasibility violation    = ", KNITRO.KN_get_abs_feas_error(kc))
        println("  KKT optimality violation = ", KNITRO.KN_get_abs_opt_error(kc))
    end

    # Delete the Knitro solver instance.
    KNITRO.KN_free(kc)

    @testset "Example QCQP1" begin
        @test nStatus == 0
        @test objSol ≈ 936.0
        @test x ≈ [0.0, 0.0, 8.0]
    end
end

example_qcqp1(; verbose=isdefined(Main, :KN_VERBOSE) ? KN_VERBOSE : true)


text:
module TestModelSearch

using Test
using MLJModels
import MLJBase.info

pca = info("PCA", pkg="MultivariateStats")
cnst = info("ConstantRegressor", pkg="MLJModels")

@test_throws ArgumentError info("Julia")

@info(ConstantRegressor) == cnst
@test info(Standardizer()) == info("Standardizer", pkg="MLJModels")

@testset "localmodels" begin
    tree = info("DecisionTreeRegressor")
    @test cnst in localmodels(modl=TestModelSearch)
    @test !(tree in localmodels(modl=TestModelSearch))
    import MLJDecisionTreeInterface.DecisionTreeRegressor
    @test tree in localmodels(modl=TestModelSearch)
end

@testset "models() and localmodels()" begin
    t(model) = model.is_pure_julia
    mods = models(t)
    @test pca in mods
    @test cnst in mods
    @test !(info("SVC") in mods)
    mods = localmodels(t, modl=TestModelSearch)
    @test cnst in mods
    @test !(pca in mods)
    u(model) = !(model.is_supervised)
    @test pca in models(u, t)
    @test !(cnst in models(u, t))
end

@testset "models(needle::Union{AbstractString,Regex}) and localmodels(needle::Union{AbstractString,Regex})" begin
    @test pca in models("PCA")
    @test pca ∉ models("PCA′")

    @test pca in models(r"PCA")
    @test pca in models(r"pca"i)
    @test pca ∉ models(r"PCA′")

    info("DecisionTreeRegressor") in localmodels("Decision"; modl = TestModelSearch)
    info("DecisionTreeRegressor") in localmodels(r"Decision"; modl = TestModelSearch)
end

end
true


text:
#####
##### Reorder
#####

reorder(from::Memory{T}) where {T} = reorder(memorydesc(T, size(from)), from)
function reorder(md::MemoryDesc, from::Memory)
    to = similar(from, eltype(from), size(from), md)
    reorder!(to, from)
    return kernel_exit_hook(to)
end

function reorder!(to::Memory, from::Memory, attributes = noattributes())
    temp_primitive(
        Lib.dnnl_reorder_primitive_desc_create,
        from,
        global_engine(),
        to,
        global_engine(),
        attributes,
    ) do p, _
        execute!(p, @dnnl_args to from)
    end
    return to
end

# Pullback is simply the identity.
# Let downstream kernels decide if they want to reorder the sensitivities.
function ChainRulesCore.rrule(
    ::typeof(reorder), from::Memory
)
    to = reorder(from)
    function reorder_pullback(Δ)
        return (ChainRulesCore.NoTangent(), Δ)
    end
    return to, reorder_pullback
end

function ChainRulesCore.rrule(
    ::typeof(reorder), desc::MemoryDesc, from::Memory
)
    to = reorder(desc, from)
    function reorder_pullback(Δ)
        return (ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent(), Δ)
    end
    return to, reorder_pullback
end

toeltype(::Type{T}, from::OneDNN.Memory{T}) where {T} = from
function toeltype(::Type{T}, from::OneDNN.Memory) where {T}
    to = similar(from, T)
    return reorder!(to, from)
end

function ChainRulesCore.rrule(::typeof(toeltype), ::Type{T}, from::OneDNN.Memory{U}) where {T,U}
    to = toeltype(T, from)
    function toeltype_pullback(Δ)
        return (ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent(), toeltype(U, Δ))
    end
    return to, toeltype_pullback
end

#####
##### Eltwise
#####

eltwise(f::F, src::Memory) where {F} = eltwise(src, forward_expand(f)...)
eltwise(::typeof(identity), src::Memory) = src
function eltwise(
    src::Memory, kind::Lib.dnnl_alg_kind_t, alpha = one(Float32), beta = zero(Float32)
)
    # Keep similar format to source
    dst = similar(src)
    eltwise!(dst, src, kind, alpha, beta)
    return kernel_exit_hook(dst)
end

function eltwise!(
    dst::Memory,
    src::Memory,
    algo::Lib.dnnl_alg_kind_t,
    alpha = one(Float32),
    beta = zero(Float32),
)
    opdesc = Ref{Lib.dnnl_eltwise_desc_t}()
    @apicall dnnl_eltwise_forward_desc_init(opdesc, Inference(), algo, src, alpha, beta)

    temp_primitive(opdesc, noattributes(), global_engine(), noforward()) do primitive, _
        execute!(primitive, @dnnl_args dst src)
    end
    return dst
end

#####
##### Backward Eltwise
#####

eltwise_backward(::typeof(identity), diff_data::Memory, data::Memory) = diff_data
function eltwise_backward(f::F, diff_data::Memory, data::Memory) where {F}
    return eltwise_backward(diff_data, data, backward_expand(f)...)
end

function eltwise_backward(
    diff_data::Memory,
    data::Memory,
    kind::Lib.dnnl_alg_kind_t,
    alpha = one(Float32),
    beta = zero(Float32),
    use_dst_for_bwd::Bool = false,
)
    diff_src = similar(diff_data)
    eltwise_backward!(diff_src, diff_data, data, kind, alpha, beta, use_dst_for_bwd)
    return kernel_exit_hook(diff_src)
end

function eltwise_backward!(
    diff_src::Memory,
    diff_dst::Memory,
    data::Memory,
    kind::Lib.dnnl_alg_kind_t,
    alpha = one(Float32),
    beta = zero(Float32),
    use_dst_for_bwd::Bool = false,
)
    opdesc = Ref{Lib.dnnl_eltwise_desc_t}()
    @apicall dnnl_eltwise_backward_desc_init(opdesc, kind, diff_dst, data, alpha, beta)

    if use_dst_for_bwd
        dst = data
        args = @dnnl_args diff_src diff_dst dst
    else
        src = data
        args = @dnnl_args diff_src diff_dst src
    end

    temp_primitive(opdesc, noattributes(), global_engine(), noforward()) do p, _
        execute!(p, args)
    end
    return diff_src
end

#####
##### Convenience Wrappers.
#####

# Aliases
# N.B.: Must be kept insync with tests!
Base.identity(x::Memory) = x
const ELTWISE_ALIASES = [
    :(Base.abs), :(Flux.sigmoid), :(Base.sqrt), :(Flux.relu), :(Base.log)
]
for op in ELTWISE_ALIASES
    @eval $op(x::Memory) = eltwise($op, x)
end

# Conversion to OneDNN forward eltwise ops
struct Linear
    α::Float32
    β::Float32
end

forward_expand(::typeof(Base.abs)) = (Lib.dnnl_eltwise_abs, zero(Float32), zero(Float32))
forward_expand(x::Linear) = (Lib.dnnl_eltwise_linear, x.α, x.β)
function forward_expand(::typeof(Flux.sigmoid))
    return (Lib.dnnl_eltwise_logistic, zero(Float32), zero(Float32))
end
forward_expand(::typeof(Base.sqrt)) = (Lib.dnnl_eltwise_sqrt, zero(Float32), zero(Float32))
forward_expand(::typeof(Flux.relu)) = (Lib.dnnl_eltwise_relu, zero(Float32), zero(Float32))
forward_expand(::typeof(Base.log)) = (Lib.dnnl_eltwise_log, zero(Float32), zero(Float32))

# Conversion to OneDNN backward eltwise ops
function backward_expand(::typeof(Base.abs))
    return (Lib.dnnl_eltwise_abs, zero(Float32), zero(Float32), false)
end
backward_expand(x::Linear) = (Lib.dnnl_eltwise_linear, x.α, x.β, false)
function backward_expand(::typeof(Flux.sigmoid))
    return (Lib.dnnl_eltwise_logistic_use_dst_for_bwd, zero(Float32), zero(Float32), true)
end
function backward_expand(::typeof(Base.sqrt))
    return (Lib.dnnl_eltwise_sqrt, zero(Float32), zero(Float32), false)
end
function backward_expand(::typeof(Flux.relu))
    return (Lib.dnnl_eltwise_relu_use_dst_for_bwd, zero(Float32), zero(Float32), true)
end
function backward_expand(::typeof(Base.log))
    return (Lib.dnnl_eltwise_log, zero(Float32), zero(Float32), false)
end

# Can this elementwise op be fused into a previous operation?
# Only if the backward op can operate on the destination tensor only.
dst_for_bwd(f::F) where {F} = last(backward_expand(f))
canfuse(f::F) where {F} = dst_for_bwd(f)
canfuse(::typeof(identity)) = true

# Forward backprop definition to the main entry point.
for op in ELTWISE_ALIASES
    @eval function ChainRulesCore.rrule(::typeof($op), x::Memory)
        y, pullback = ChainRulesCore.rrule(eltwise, $op, x)
        return y, dy -> Base.tail(pullback(dy))
    end
end

function ChainRulesCore.rrule(::typeof(eltwise), f::F, from::Memory) where {F}
    to = eltwise(f, from)
    data = dst_for_bwd(f) ? to : from
    function eltwise_pullback(Δ)
        diff_from = eltwise_backward(f, Memory(Δ), data)
        return (ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent(), diff_from)
    end
    return to, eltwise_pullback
end

#####
##### Binary
#####

# Note: No need to define `rrule`s here since higher level definitions should take
# care of this automatically.
function binary(f::F, src_0::Memory, src_1::Memory) where {F}
    return binary(src_0, src_1, binary_forward(f))
end
function binary(src_0::Memory, src_1::Memory, kind)
    dst = similar(src_0)
    binary!(dst, src_0, src_1, kind)
    return kernel_exit_hook(dst)
end

function binary!(f::F, dst::Memory, src_0::Memory, src_1::Memory) where {F}
    return binary!(dst, src_0, src_1, binary_forward(f))
end
function binary!(dst::Memory, src_0::Memory, src_1::Memory, kind)
    opdesc = Ref{Lib.dnnl_binary_desc_t}()
    @apicall dnnl_binary_desc_init(opdesc, kind, src_0, src_1, dst)

    args = @dnnl_args dst src_0 src_1
    temp_primitive(opdesc, noattributes(), global_engine(), noforward()) do p, _
        execute!(p, args)
    end
    return dst
end

Base.:+(a::Memory, b::Memory) = binary(+, a, b)
Base.:-(a::Memory, b::Memory) = binary(-, a, b)
Base.broadcasted(::typeof(+), a::Memory, b::Memory) = binary(+, a, b)
Base.broadcasted(::typeof(-), a::Memory, b::Memory) = binary(-, a, b)
Base.broadcasted(::typeof(*), a::Memory, b::Memory) = binary(*, a, b)
Base.broadcasted(::typeof(/), a::Memory, b::Memory) = binary(/, a, b)
Base.broadcasted(::typeof(min), a::Memory, b::Memory) = binary(min, a, b)
Base.broadcasted(::typeof(max), a::Memory, b::Memory) = binary(max, a, b)

binary_forward(::typeof(+)) = Lib.dnnl_binary_add
binary_forward(::typeof(*)) = Lib.dnnl_binary_mul
binary_forward(::typeof(-)) = Lib.dnnl_binary_sub
binary_forward(::typeof(/)) = Lib.dnnl_binary_div
binary_forward(::typeof(max)) = Lib.dnnl_binary_max
binary_forward(::typeof(min)) = Lib.dnnl_binary_min

#####
##### Optimized Format Conversion.
#####

maybe_reorder(md::MemoryDesc, x::Memory) = (md == memorydesc(x)) ? x : reorder(md, x)
function maybe_reorder(descriptor::PrimitiveDescriptor, x::Memory, key)
    return maybe_reorder(query_md(descriptor, key), x)
end

reorder_if_any(::Ptr{MemoryDesc}, ::PrimitiveDescriptor, x::Memory, key) = x
function reorder_if_any(md::MemoryDesc, desc::PrimitiveDescriptor, x::Memory, key)
    return isany(md) ? maybe_reorder(desc, x, key) : x
end

#####
##### Logsoftmax
#####

function logsoftmax(src::Memory{T,N}, axis::Integer) where {T,N}
    opdesc = Ref{Lib.dnnl_logsoftmax_desc_t}()
    @apicall dnnl_logsoftmax_forward_desc_init(opdesc, Inference(), src, N - axis)

    return temp_primitive(opdesc, noattributes(), global_engine(), noforward()) do p, _
        dst = similar(src)
        execute!(p, @dnnl_args dst src)
        return dst
    end
end

function logsoftmax_backward(diff_dst::Memory{T,N}, dst::Memory{T,N}, axis::Integer) where {T,N}
    opdesc = Ref{Lib.dnnl_logsoftmax_desc_t}()
    @apicall dnnl_logsoftmax_backward_desc_init(opdesc, diff_dst, dst, N - axis)

    return temp_primitive(opdesc, noattributes(), global_engine(), noforward()) do p, _
        diff_src = similar(diff_dst)
        execute!(p, @dnnl_args diff_dst dst diff_src)
        return diff_src
    end
end

function ChainRulesCore.rrule(::typeof(logsoftmax), src::Memory, axis)
    dst = logsoftmax(src, axis)
    function logsoftmax_pullback(diff_dst)
        diff_src = logsoftmax_backward(diff_dst, dst, axis)
        return (
            ChainRulesCore.NoTangent(),
            diff_src,
            ChainRulesCore.NoTangent(),
        )
    end
    return dst, logsoftmax_pullback
end


text:
# ---
# title: 878. Nth Magical Number
# id: problem878
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Binary Search
# link: <https://leetcode.com/problems/nth-magical-number/description/>
# hidden: true
# ---
#
# A positive integer is _magical_  if it is divisible by either A or B.
#
# Return the N-th magical number.  Since the answer may be very large, **return
# it modulo**`10^9 + 7`.
#
#
#
# **Example 1:**
#
#
#
#     Input: N = 1, A = 2, B = 3
#     Output: 2
#
#
# **Example 2:**
#
#
#
#     Input: N = 4, A = 2, B = 3
#     Output: 6
#
#
# **Example 3:**
#
#
#
#     Input: N = 5, A = 2, B = 4
#     Output: 10
#
#
# **Example 4:**
#
#
#
#     Input: N = 3, A = 6, B = 4
#     Output: 8
#
#
#
#
# **Note:**
#
#   1. `1 <= N <= 10^9`
#   2. `2 <= A <= 40000`
#   3. `2 <= B <= 40000`
#
#
## @lc code=start
using LeetCode

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


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

name = "GAP_lib"
version = v"4.11.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/gap-system/gap.git", "069a6497424113dc5e5ffe133be9da5bfe6acb24"),
#    ArchiveSource("https://github.com/gap-system/gap/releases/download/v$(version)/gap-$(version)-core.tar.bz2",
#                  "6637f66409bc91af21eaa38368153270b71b13b55b75cc1550ed867c629901d1"),
    ArchiveSource("https://github.com/gap-system/gap/releases/download/v$(version)/packages-required-v$(version).tar.gz",
                  "29ab6e2752f39d22e3f0a19e5bcbfec1710993b4cfd52337a0ca10fda6a76537";
                  unpack_target="pkg"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/gap*

mv ../pkg .
find pkg -name '._*' -exec rm \{\} \; # unwanted files

# run autogen.sh if compiling from it source and/or if configure was patched
./autogen.sh

# compile a native version of GAP so we can use it to generate the manual
# (the manual is only in FULL gap release tarballs, not in the -core tarball
# nor in git snapshots)
mkdir native-build
cd native-build
#apk add gmp-dev zlib-dev
../configure --build=${MACHTYPE} --host=${MACHTYPE} CC=${CC_BUILD} CXX=${CXX_BUILD} --with-zlib=${prefix} --with-gmp=${prefix}
make -j${nproc}
make html   # build the manual (only HTML and txt; for PDF we'd need LaTeX)
cd ..

# remove the native build, it has done its job
rm -rf native-build

# the license
install_license LICENSE

# "install" most of the files
rm -rf autom4te.cache dev etc extern hpcgap/extern pkg
mkdir -p ${prefix}/share/gap/
mv * ${prefix}/share/gap/
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

# The products that we will ensure are always built
products = Product[
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Zlib_jll"),
    BuildDependency("GMP_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: vchuravy/Yggdrasil

max_issues_repo_name: vchuravy/Yggdrasil

max_forks_repo_name: vchuravy/Yggdrasil