text:
function load_constrained(op::Operation, u₁loop::Symbol, u₂loop::Symbol, innermost_loop_or_vloop::Symbol, forprefetch::Bool = false)
    dependsonu₁ = isu₁unrolled(op)
    dependsonu₂ = isu₂unrolled(op)
    if forprefetch
        (dependsonu₁ & dependsonu₂) || return false
    end
    unrolleddeps = Symbol[]
    if forprefetch # innermost_loop_or_vloop is innermost_loop
        push!(unrolleddeps, innermost_loop_or_vloop)
    else# if this is for `sub_fmas` instead of prefetch, then `innermost_loop_or_vloop` is `vloop`, and we only care if it isn't vectorized
        # if it is vectorized, it can't be broadcasted on AVX512 anyway, so no point trying to enable that optimization
        dependsonu₁ &= u₁loop !== innermost_loop_or_vloop
        dependsonu₂ &= u₂loop !== innermost_loop_or_vloop
    end
    dependsonu₁ && push!(unrolleddeps, u₁loop)
    dependsonu₂ && push!(unrolleddeps, u₂loop)
    length(unrolleddeps) > 0 || return false
    any(parents(op)) do opp
        isload(opp) && all(in(loopdependencies(opp)), unrolleddeps)
    end
end
function check_if_remfirst(ls::LoopSet, ua::UnrollArgs)
    usorig = ls.unrollspecification
    @unpack u₁, u₁loopsym, u₂loopsym, u₂max = ua
    u₁loop = getloop(ls, u₁loopsym)
    u₂loop = getloop(ls, u₂loopsym)
    if isstaticloop(u₁loop) && (usorig.u₁ != u₁)
        return true
    end
    if isstaticloop(u₂loop) && (usorig.u₂ != u₂max)
        return true
    end
    false
end
function sub_fmas(ls::LoopSet, op::Operation, ua::UnrollArgs)
    @unpack u₁, u₁loopsym, u₂loopsym, vloopsym, u₂max = ua
    !(load_constrained(op, u₁loopsym, u₂loopsym, vloopsym) || check_if_remfirst(ls, ua))
end

function parent_unroll_status(op::Operation, u₁loop::Symbol, us::UnrollSpecification)
    parentsop = parents(op)
    u2 = fill(false, length(parentsop))
    u1 = similar(u2)
    for i ∈ eachindex(parentsop)
        u1[i] = isunrolled_sym(parentsop[i], u₁loop, us)
    end
    u1, u2
end
function parent_unroll_status(op::Operation, u₁loop::Symbol, u₂loop::Symbol, vloop::Symbol, u₂max::Int, us::UnrollSpecification)
    u₂max == -1 && return parent_unroll_status(op, u₁loop, us)
    vparents = parents(op);
    # parent_names = Vector{Symbol}(undef, length(vparents))
    parents_u₁syms = Vector{Bool}(undef, length(vparents))
    parents_u₂syms = Vector{Bool}(undef, length(vparents))
    for i ∈ eachindex(vparents)
        parents_u₁syms[i], parents_u₂syms[i] = isunrolled_sym(vparents[i], u₁loop, u₂loop, vloop, us)#, u₂max)
    end
    # parent_names, parents_u₁syms, parents_u₂syms
    parents_u₁syms, parents_u₂syms
end

function _add_loopvalue!(ex::Expr, loopval::Symbol, vloop::Loop, u::Int)
  vloopsym = vloop.itersymbol
  if loopval === vloopsym
    if iszero(u)
      push!(ex.args, _MMind(loopval, step(vloop)))
    else
      vstep = step(vloop)
      mm = _MMind(loopval, vstep)
      if isone(u) & isone(vstep)
        push!(ex.args, Expr(:call, lv(:vadd_nsw), VECTORWIDTHSYMBOL, mm))
      else
        push!(ex.args, Expr(:call, lv(:vadd_nsw), mulexpr(VECTORWIDTHSYMBOL, u, vstep), mm))
      end
    end
  elseif u == 0
    push!(ex.args, loopval)
  else
    push!(ex.args, Expr(:call, lv(:vadd_nsw), loopval, staticexpr(u)))
  end
end
function add_loopvalue!(instrcall::Expr, loopval, ua::UnrollArgs, u₁::Int)
  @unpack u₁loopsym, u₂loopsym, vloopsym, vloop, suffix = ua
  if loopval === u₁loopsym #parentsunrolled[n]
    if isone(u₁)
      _add_loopvalue!(instrcall, loopval, vloop, 0)
    else
      t = Expr(:tuple)
      for u ∈ 0:u₁-1
        _add_loopvalue!(t, loopval, vloop, u)
      end
      push!(instrcall.args, Expr(:call, lv(:VecUnroll), t))
    end
  elseif suffix > 0 && loopval === u₂loopsym
    _add_loopvalue!(instrcall, loopval, vloop, suffix)
  elseif loopval === vloopsym
    push!(instrcall.args, _MMind(loopval, step(vloop)))
  else
    push!(instrcall.args, loopval)
  end
end

vecunrolllen(::Type{VecUnroll{N,W,T,V}}) where {N,W,T,V} = (N::Int + 1)
vecunrolllen(_) = -1
function ifelselastexpr(hasf::Bool, M::Int, vargtypes, K::Int, S::Int, maskearly::Bool)
    q = Expr(:block, Expr(:meta,:inline))
    vargs = Vector{Symbol}(undef, K)
    for k ∈ 1:K
        vargs[k] = Symbol(:varg_,k)
    end
    lengths = Vector{Int}(undef, K);
    for k ∈ 1:K
        lengths[k] = l = vecunrolllen(vargtypes[k])
        if hasf
            gfvarg = Expr(:call, GlobalRef(Core, :getfield), :vargs, k, false)
            if l ≠ -1 # VecUnroll
                gfvarg = Expr(:call, GlobalRef(Core, :getfield), gfvarg, 1, false)
            end
            push!(q.args, Expr(:(=), vargs[k], gfvarg))
        elseif l ≠ -1
            varg = vargs[k]
            vargs[k] = dvarg = Symbol(:d, varg)
            push!(q.args, :($dvarg = data($varg)))
        end
    end
    N = last(lengths)
    start = (hasf | maskearly) ? 1 : M
    Sreduced = (S > 0) && (lengths[S] == -1)
    if Sreduced
        maxlen = maximum(lengths)
        if maxlen == -1
            Sreduced = false
            t = Expr(:tuple)
        else
            hasf || throw(ArgumentError("Argument reduction only supported for `ifelse(last/partial)(f::Function, args...)`"))
            M = maxlen
            t = q
        end
    else
        t = Expr(:tuple)
    end
    for m ∈ 1:start-1
        push!(t.args, :(getfield($(vargs[1]), $m, false)))
    end
    for m ∈ start:M
        call = if hasf
            (maskearly | (m == M)) ? Expr(:call, :ifelse, :f, :m) : Expr(:call, :f)
        else# m == M because !hasf
            Expr(:call, :ifelse, :m)
        end
        for k ∈ 1:K
            if lengths[k] == -1
                push!(call.args, vargs[k])
            else
                # @assert (k == K) || (lengths[k] == M)
                push!(call.args, :(getfield($(vargs[k]), $m, false)))
            end
        end
        if Sreduced
            push!(t.args, Expr(:(=), vargs[S], call))
        elseif N == -1
            push!(q.args, call)
            return q
        else
            push!(t.args, call)
        end
    end
    Sreduced && return q
    for m ∈ M+1:N
        push!(t.args, :(getfield($(vargs[K]), $m, false)))
    end
    # push!(q.args, :(VecUnroll($t)::VecUnroll{$N,$W,$T,$V}))
    # push!(q.args, Expr(:call, lv(:VecUnroll), t))
    push!(q.args, :(VecUnroll($t)))
    q
end
@generated function ifelselast(f::F, m::AbstractMask{W}, ::StaticInt{M}, ::StaticInt{S}, vargs::Vararg{Any,K}) where {F,W,K,M,S}
    ifelselastexpr(true, M, vargs, K, S, false)
end
@generated function ifelselast(m::AbstractMask{W}, ::StaticInt{M}, ::StaticInt{S}, varg_1::V1, varg_2::V2) where {W,V1,V2,M,S}
    ifelselastexpr(false, M, (V1,V2), 2, S, false)
end
@generated function ifelsepartial(f::F, m::AbstractMask{W}, ::StaticInt{M}, ::StaticInt{S}, vargs::Vararg{Any,K}) where {F,W,K,M,S}
    ifelselastexpr(true, M, vargs, K, S, true)
end
@generated function ifelsepartial(m::AbstractMask{W}, ::StaticInt{M}, ::StaticInt{S}, varg_1::V1, varg_2::V2) where {W,V1,V2,M,S}
    ifelselastexpr(false, M, (V1,V2), 2, S, true)
end
# @inline ifelselast(f::F, m::AbstractMask{W}, ::StaticInt{M}, ::StaticInt{S}, vargs::Vararg{NativeTypes,K}) where {F,W,K,M,S} = f(vargs...)
# @inline ifelsepartial(f::F, m::AbstractMask{W}, ::StaticInt{M}, ::StaticInt{S}, vargs::Vararg{NativeTypes,K}) where {F,W,K,M,S} = f(vargs...)
@generated function subset_vec_unroll(vu::VecUnroll{N}, ::StaticInt{S}) where {N,S}
    (1 ≤ S ≤ N + 1) || throw(ArgumentError("`vu` isa `VecUnroll` of `$(N+1)` elements, but trying to subset $S of them."))
    t = Expr(:tuple)
    gf = GlobalRef(Core,:getfield)
    S == 1 && return Expr(:block, Expr(:meta,:inline), :($gf($gf(vu,1),1,false)))
    for s ∈ 1:S
        push!(t.args, Expr(:call, gf, :vud, s, false))
    end
    quote
        $(Expr(:meta,:inline))
        vud = $gf(vu, 1)
        VecUnroll($t)
    end
end
# `S` is the ind to replace with the return value of previous invocation ("S" for "self") if reducing
@generated function partialmap(f::F, default::D, ::StaticInt{M}, ::StaticInt{S}, vargs::Vararg{Any,K}) where {F,M,K,D,S}
    lengths = Vector{Int}(undef, K);
    q = Expr(:block, Expr(:meta,:inline))
    syms = Vector{Symbol}(undef, K)
    isnotpartial = true
    gf = GlobalRef(Core, :getfield)
    for k ∈ 1:K
        l = vecunrolllen(vargs[k])
        # if l
        kisnotpartial = ((l ≡ -1) & (k ≢ S)) | (l ≡ M)
        isnotpartial &= kisnotpartial
        lengths[k] = l
        @assert (l == -1) || (l ≥ M)
        syms[k] = symk = Symbol(:vargs_,k)
        extractq = :($gf(vargs, $k, false))
        if l != -1
            extractq = :(data($extractq))
        end
        push!(q.args, :($symk = $extractq))
    end
    Dlen = vecunrolllen(D)
    N = maximum(lengths)
    Sreduced = (S > 0) && (lengths[S] == -1) && N != -1
    if isnotpartial & (Sreduced | (Dlen == N))
        q =  Expr(:call, :f)
        for k ∈ 1:K
            push!(q.args, :($gf(vargs, $k, false)))
        end
        return Expr(:block, Expr(:meta, :inline), q)
    end
    if Sreduced
        M = N
        t = q
    else
        @assert (N ≤ Dlen)
        if Dlen == -1
            @assert (M == 1)
        else
            push!(q.args, :(dd = data(default)))
        end
        t = Expr(:tuple)
    end
    for m ∈ 1:M
        call = Expr(:call, :f)
        for k ∈ 1:K
            if lengths[k] == -1
                push!(call.args, syms[k])
            else
                push!(call.args, Expr(:call, gf, syms[k], m, false))
            end
        end
        # minimal change in behavior to fix case when !Sreduced by N -> Dlen; TODO: what should Dlen be here?
        if Sreduced ? (N == -1) : (Dlen == -1)
            push!(q.args, call)
            return q
        end
        if Sreduced
            push!(t.args, Expr(:(=), syms[S], call))
        else
            push!(t.args, call)
        end
    end
    Sreduced && return q
    for m ∈ M+1:max(N,Dlen)
        push!(t.args, :($gf(dd, $m, false)))
    end
    push!(q.args, :(VecUnroll($t)))
    q
end

function parent_op_name!(
  q, ls::LoopSet, parents_op::Vector{Operation}, n::Int, modsuffix, suffix, parents_u₁syms, parents_u₂syms, u₁, u₂max, u₂unrolledsym, op, tiledouterreduction
)
  opp = parents_op[n]
  opisvectorized = isvectorized(op)
  parent = mangledvar(opp)
  u = 0
  if n == tiledouterreduction# && isvectorized(opp)
    parent = if u₂unrolledsym #!parents_u₁syms[n]
      Symbol(parent, modsuffix)
    else
      Symbol(parent, '_', modsuffix)
    end
  else
    isouterreduct = isouterreduction(ls, opp) ≠ -1
    u = if !parents_u₁syms[n]
      1
    elseif isouterreduct
      getu₁full(ls, u₁)
    else
      getu₁forreduct(ls, opp, u₁)
    end
    if parents_u₂syms[n]
      if isu₂unrolled(op) # u₂unrolledsym ||
        parent = isouterreduct ? Symbol(parent, suffix) : Symbol(parent, suffix, '_', '_', u)
      elseif u₂max > 1
        t = Expr(:tuple)
        reduction = Expr(:call, GlobalRef(ArrayInterface, :reduce_tup), reduce_to_onevecunroll(opp), t)
        for u₂ ∈ 0:u₂max-1
          push!(t.args, Symbol(parent, u₂, '_', '_', u))
        end
        parent = gensym!(ls, parent)
        push!(q.args, Expr(:(=), parent, reduction))
        parent
      else
        # parent = Symbol(parent, '_', u)
        parent = Symbol(parent, 0, "__", u)
      end
    else
      parent = Symbol(parent, '_', u)
    end
  end
  if opisvectorized && isload(opp) && (!isvectorized(opp))
    parent = Symbol(parent, "##broadcasted##")
  end
  parent, u
end
function getuouterreduct(ls::LoopSet, op::Operation, suffix)
    us = ls.unrollspecification
    if us.vloopnum === us.u₁loopnum # unroll u₂
        suffix
    else # unroll u₁
        us.u₁
    end
end

function getu₁full(ls::LoopSet, u₁::Int)
  Ureduct = ureduct(ls)
  ufull = ls.unrollspecification.u₁
  ufull = Core.ifelse(Ureduct == -1, ufull, min(ufull, Ureduct))
    # ufull = if Ureduct == -1 # no reducing
    #     ls.unrollspecification.u₁
    # else
    #     Ureduct
    # end
    # max is because we may be in the extended (non-reduct) region
  return max(u₁, ufull)
end
function getu₁forreduct(ls::LoopSet, op::Operation, u₁::Int)
    !isu₁unrolled(op) && return 1
    # if `op` is u₁unrolled, we must then find out if the initialization is `u₁unrolled`
    # if it is, then op's `u₁` will be the current `u₁`
    # if it is not, then the opp is initialized once per full u₁
    while true
        selfparentid = findfirst(Base.Fix2(===,name(op)) ∘ name, parents(op))
        selfparentid === nothing && return u₁
        op = parents(op)[selfparentid]
        isreduction(op) || break
    end
    if isu₁unrolled(op)
        return u₁
    elseif (ls.unrollspecification.u₂ != -1) && length(ls.outer_reductions) > 0
        # then `ureduct` doesn't tell us what we need, so
        return ls.unrollspecification.u₁
    else # we need to find u₁-full
        return getu₁full(ls, u₁)
    end
end
isidentityop(op::Operation) = iscompute(op) && (instruction(op).instr === :identity) && (length(parents(op)) == 1)
function reduce_parent!(q::Expr, ls::LoopSet, op::Operation, opp::Operation, parent::Symbol)
  isvectorized(op) && return parent
  if isvectorized(opp)
    oppt = opp
  elseif isidentityop(opp)
    oppt = parents(opp)[1]
    isvectorized(oppt) || return parent
  else
    return parent
  end
  reduct_class = reduction_instruction_class(oppt.instruction)
  if (instruction(op).instr === :mul_fast) & (reduct_class == ADDITIVE_IN_REDUCTIONS)
    op.vectorized = true
    return parent
  end
  newp = gensym(parent)
  if instruction(op).instr ≢ :ifelse
    push!(q.args, Expr(:(=), newp, Expr(:call, lv(reduction_to_scalar(reduct_class)), parent)))#IfElseReducer
  else
    reductexpr = ifelse_reduction(:IfElseReducer,op) do opv
      throw(LoopError("Does not support storing mirrored ifelse-reductions yet"))
    end
    push!(q.args, Expr(:(=), newp, Expr(:call, reductexpr, parent)))
  end
  newp
end
function lower_compute!(
  q::Expr, op::Operation, ls::LoopSet, ua::UnrollArgs, mask::Bool
  )
  @unpack u₁, u₁loopsym, u₂loopsym, vloopsym, u₂max, suffix = ua
  var = name(op)
  instr = instruction(op)
  parents_op = parents(op)
  nparents = length(parents_op)
  # __u₂max = ls.unrollspecification.u₂
  # TODO: perhaps allos for swithcing unrolled axis again
  mvar, u₁unrolledsym, u₂unrolledsym = variable_name_and_unrolled(op, u₁loopsym, u₂loopsym, vloopsym, suffix, ls)
  opunrolled = u₁unrolledsym || isu₁unrolled(op)
  us = ls.unrollspecification
  parents_u₁syms, parents_u₂syms = parent_unroll_status(op, u₁loopsym, u₂loopsym, vloopsym, u₂max, us)
  # tiledouterreduction = if num_loops(ls) == 1# (suffix == -1)# || (vloopsym === u₂loopsym)
  tiledouterreduction = if (suffix == -1)# || (vloopsym === u₂loopsym)
    suffix_ = Symbol("")
    -1
  else
    suffix_ = Symbol(suffix, :_)
    isouterreduction(ls, op)
  end
  if !opunrolled && any(parents_u₁syms) # TODO: Clean up this mess, refactor the naming code, putting it in one place and have everywhere else use it for easy equivalence.
    parents_op = copy(parents_op) # don't mutate the original!
    for i ∈ eachindex(parents_u₁syms)
      parents_u₁syms[i] || continue
      parents_u₁syms[i] = false
      parentop = parents_op[i]
      i == tiledouterreduction && isconstant(parentop) && continue
      newparentop = Operation(
        parentop.identifier, gensym(parentop.variable), parentop.elementbytes, parentop.instruction, parentop.node_type,
        parentop.dependencies, parentop.reduced_deps, parentop.parents, parentop.ref, parentop.reduced_children
      )
      newparentop.vectorized = false
      newparentop.u₁unrolled = false
      newparentop.u₂unrolled = parents_u₂syms[i]
      parentname = mangledvar(parentop)
      newparentname = mangledvar(newparentop)
      parents_op[i] = newparentop
      if parents_u₂syms[i]
        parentname = Symbol(parentname, suffix_)
        newparentname = Symbol(newparentname, suffix_)
      end
      if isconstant(newparentop)
        push!(q.args, Expr(:(=), Symbol(newparentname, '_', 1), Symbol(parentname, '_', 1)))
      else
        newpname = Symbol(newparentname, '_', u₁)
        push!(q.args, Expr(:(=), newpname, Symbol(parentname, '_', u₁)))
        reduce_expr!(q, newparentname, newparentop, u₁, -1, true, false)
        push!(q.args, Expr(:(=), Symbol(newparentname, '_', 1), Symbol(newparentname, "##onevec##")))
      end
    end
  end
  # if suffix === nothing# &&
  # end
  # cache unroll and tiling check of parents
  # not broadcasted, because we use frequent checks of individual bools
  # making BitArrays inefficient.
  # parentsyms = [opp.variable for opp ∈ parents(op)]
  Uiter = opunrolled ? u₁ - 1 : 0
  isreduct = isreduction(op)
  if Base.libllvm_version < v"11.0.0" && (suffix ≠ -1) && isreduct# && (iszero(suffix) || (ls.unrollspecification.u₂ - 1 == suffix))
    # if (length(reduceddependencies(op)) > 0) | (length(reducedchildren(op)) > 0)# && (iszero(suffix) || (ls.unrollspecification.u₂ - 1 == suffix))
    # instrfid = findfirst(isequal(instr.instr), (:vfmadd, :vfnmadd, :vfmsub, :vfnmsub))
    instrfid = findfirst(Base.Fix2(===,instr.instr), (:vfmadd_fast, :vfnmadd_fast, :vfmsub_fast, :vfnmsub_fast))
    # instrfid = findfirst(isequal(instr.instr), (:vfnmadd_fast, :vfmsub_fast, :vfnmsub_fast))
    # want to instcombine when parent load's deps are superset
    # also make sure opp is unrolled
    if !(instrfid === nothing) && (opunrolled && u₁ > 1) && sub_fmas(ls, op, ua)
      specific_fmas = Base.libllvm_version >= v"11.0.0" ? (:vfmadd, :vfnmadd, :vfmsub, :vfnmsub) : (:vfmadd231, :vfnmadd231, :vfmsub231, :vfnmsub231)
      # specific_fmas = Base.libllvm_version >= v"11.0.0" ? (:vfnmadd, :vfmsub, :vfnmsub) : (:vfnmadd231, :vfmsub231, :vfnmsub231)
      # specific_fmas = (:vfmadd231, :vfnmadd231, :vfmsub231, :vfnmsub231)
      instr = Instruction(specific_fmas[instrfid])
    end
  end
  reduceddeps = reduceddependencies(op)
  vecinreduceddeps = isreduct && vloopsym ∈ reduceddeps
  maskreduct = mask & vecinreduceddeps #any(opp -> opp.variable === var, parents_op)
  # if vecinreduceddeps && vectorized ∉ loopdependencies(op) # screen parent opps for those needing a reduction to scalar
  #     # parents_op = reduce_vectorized_parents!(q, op, parents_op, U, u₁loopsym, u₂loopsym, vectorized, suffix)
  #     isreducingidentity!(q, op, parents_op, U, u₁loopsym, u₂loopsym, vectorized, suffix) && return
  # end
  # if a parent is not unrolled, the compiler should handle broadcasting CSE.
  # because unrolled/tiled parents result in an unrolled/tiled dependendency,
  # we handle both the tiled and untiled case here.
  # bajillion branches that go the same way on each iteration
  # but smaller function is probably worthwhile. Compiler could theoreically split anyway
  # but I suspect that the branches are so cheap compared to the cost of everything else going on
  # that smaller size is more advantageous.
  opisvectorized = isvectorized(op)
  modsuffix = 0
  # for u ∈ 0:Uiter
  isouterreduct = false
  instrcall = callexpr(instr)
  dopartialmap = false
  varsym = if tiledouterreduction > 0 # then suffix ≠ -1
    # modsuffix = ((u + suffix*(Uiter + 1)) & 7)
    isouterreduct = true
    # if u₁unrolledsym
    #     modsuffix = ls.unrollspecification.u₁#getu₁full(ls, u₁)#u₁
    #     Symbol(mangledvar(op), '_', modsuffix)
    # else
    if u₁unrolledsym
      # modsuffix = 0
      modsuffix = ls.unrollspecification.u₁
      Symbol(mangledvar(op), '_', modsuffix)
    else
      modsuffix = suffix % ls.ureduct
      Symbol(mangledvar(op), modsuffix)
    end
  elseif u₁unrolledsym
    if isreduct #(isanouterreduction(ls, op))
      # isouterreduct = true
      isouterreduct = isanouterreduction(ls, op)
      u₁reduct = isouterreduct ? getu₁full(ls, u₁) : getu₁forreduct(ls, op, u₁)
      dopartialmap = u₁reduct ≠ u₁
      Symbol(mvar, '_', u₁reduct)
    else
      Symbol(mvar, '_', u₁)
    end
  else
    Symbol(mvar, '_', 1)
  end
  selfopname = varsym
  selfdep = 0
  for n ∈ 1:nparents
    opp = parents_op[n]
    if isloopvalue(opp)
      loopval = first(loopdependencies(opp))
      add_loopvalue!(instrcall, loopval, ua, u₁)
    elseif name(opp) === name(op)
      selfdep = n
      if ((isvectorized(opp) && !isvectorized(op))) ||
        (parents_u₁syms[n] != u₁unrolledsym) || (parents_u₂syms[n] != u₂unrolledsym)

        selfopname, uₚ = parent_op_name!(q, ls, parents_op, n, modsuffix, suffix, parents_u₁syms, parents_u₂syms, u₁, u₂max, u₂unrolledsym, op, tiledouterreduction)
        push!(instrcall.args, selfopname)
      else
        push!(instrcall.args, varsym)
      end
    elseif ((!isu₂unrolled(op)) & isu₂unrolled(opp)) && (parents_u₂syms[n] & (!u₂unrolledsym))
      # elseif parents_u₂syms[n] & (!u₂unrolledsym)
      #&& (isouterreduction(ls, opp) != -1)
      # this checks if the parent is u₂ unrolled but this operation is not, in which case we need to reduce it.
      reduced_u₂ = reduce_expr_u₂(mangledvar(opp), opp, u₂max, Symbol("__", u₁))#ureduct(ls))
      reducedparentname = gensym!(ls, "reducedop")
      push!(q.args, Expr(:(=), reducedparentname, reduced_u₂))
      reduced_u₂ = reduce_parent!(q, ls, op, opp, reducedparentname)
      push!(instrcall.args, reduced_u₂)
    elseif isconstant(opp) && instruction(opp).mod === GLOBALCONSTANT
      push!(instrcall.args, GlobalRef(Base, instruction(opp).instr))
    else
      # 0 for tiledouterreduction, because `name(opp) ≢ name(op)`
      parent, uₚ = parent_op_name!(q, ls, parents_op, n, modsuffix, suffix, parents_u₁syms, parents_u₂syms, u₁, u₂max, u₂unrolledsym, op, 0)
      parent = reduce_parent!(q, ls, op, opp, parent)
      if (selfdep == 0) && search_tree(parents(opp), name(op))
        selfdep = n
        push!(instrcall.args, parent)
      elseif (uₚ ≠ 0) & (uₚ > u₁)
        push!(instrcall.args, :(subset_vec_unroll($parent, StaticInt{$u₁}())))
      else
        push!(instrcall.args, parent)
      end
    end
  end
  selfdepreduce = ifelse(((!u₁unrolledsym) & isu₁unrolled(op)) & (u₁ > 1), selfdep, 0)
  if maskreduct
    ifelsefunc = if us.u₁ == 1
      :ifelse # don't need to be fancy
    elseif (u₁loopsym !== vloopsym)
      :ifelsepartial # ifelse all the early ones
    else# mask last u₁
      :ifelselast # ifelse only the last one
    end
    if last(instrcall.args) === varsym
      pushfirst!(instrcall.args, lv(ifelsefunc))
      # showexpr = true
      insert!(instrcall.args, 3, MASKSYMBOL)
      if !(ifelsefunc === :ifelse)
        insert!(instrcall.args, 4, staticexpr(u₁))
        insert!(instrcall.args, 5, staticexpr(selfdepreduce))
      end
    elseif all(in(loopdependencies(op)), reduceddeps) || selfdep ≠ 0#any(opp -> mangledvar(opp) === mangledvar(op), parents_op)
      # Here, we are evaluating the function, and then `ifelse`-ing it with `hasf == false`.
      # That means we still need to adjust the `instrcall` in case we're reducing/accumulating across the unroll
      if ifelsefunc ≡ :ifelse # ifelse means it's unrolled by 1, no need
        # push!(q.args, Expr(:(=), varsym, Expr(:call, lv(ifelsefunc), MASKSYMBOL, instrcall, selfopname)))
        push!(q.args, Expr(:(=), varsym, Expr(:call, lv(ifelsefunc), MASKSYMBOL, instrcall, selfopname)))
      elseif ((u₁ ≡ 1) | (selfdepreduce ≡ 0))
        # if the current unroll is 1, no need to accumulate. Same if there is no selfdepreduce, but there has to be if we're here?
        # push!(q.args, Expr(:(=), varsym, Expr(:call, lv(ifelsefunc), MASKSYMBOL, staticexpr(u₁), staticexpr(selfdepreduce), instrcall, selfopname)))
        push!(q.args, Expr(:(=), varsym, Expr(:call, lv(ifelsefunc), MASKSYMBOL, staticexpr(u₁), staticexpr(selfdepreduce), instrcall, selfopname)))
      else
        make_partial_map!(instrcall, selfopname, u₁, selfdepreduce)
        # partialmap accumulates
        push!(q.args, Expr(:(=), varsym, Expr(:call, lv(:ifelse), MASKSYMBOL, instrcall, selfopname)))
      end
      return
      # elseif selfdep != 0
      #   make_partial_map!(instrcall, selfopname, u₁, selfdepreduce)
    end
  elseif selfdep != 0 && (dopartialmap ||
    (isouterreduct && (opunrolled) && (u₁ < us.u₁)) ||
    (isreduct & (u₁ > 1) & (!u₁unrolledsym) & isu₁unrolled(op))) # TODO: DRY `selfdepreduce` definition
    # first possibility (`isouterreduct && opunrolled && (u₁ < ls.unrollspecification.u₁)`):
    # checks if we're in the "reduct" part of an outer reduction
    #
    # second possibility (`(isreduct & (u₁ > 1) & (!u₁unrolledsym) & isu₁unrolled(op))`):
    # if the operation is repeated across u₁ (indicated by `isu₁unrolled(op)`) but
    # the variables are not correspondingly replicated across u₁ (indicated by `!u₁unrolledsym`)
    # then we need to accumulate it.
    make_partial_map!(instrcall, selfopname, u₁, selfdepreduce)
    # elseif
  end
  if instr.instr === :identity && isone(length(parents_op))
    if instrcall.args[2] !== varsym
      push!(q.args, Expr(:(=), varsym, instrcall.args[2]))
    end
  elseif identifier(op) ∉ ls.outer_reductions && should_broadcast_op(op)
    push!(q.args, Expr(:(=), varsym, Expr(:call, lv(:vbroadcast), VECTORWIDTHSYMBOL, instrcall)))
  else
    push!(q.args, Expr(:(=), varsym, instrcall))
  end
  # end
end
function make_partial_map!(instrcall, selfopname, u₁, selfdep)
  pushfirst!(instrcall.args, lv(:partialmap))
  insert!(instrcall.args, 3, selfopname)
  insert!(instrcall.args, 4, staticexpr(u₁))
  insert!(instrcall.args, 5, staticexpr(selfdep))
  nothing
end


text:
v = (1, 2, 3)
@show v

ext: jl

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/iup/include/iupgl.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function IupGLCanvasOpen()
  ccall( (:IupGLCanvasOpen, iup), Void, (), )
end
function IupGLCanvas()
  ccall( (:IupGLCanvas, iup), Ptr{Cint}, (), )
end
function IupGLMakeCurrent(ih::Ptr{Cint})
  ccall( (:IupGLMakeCurrent, iup), Void, (Ptr{Cint},), ih)
end
function IupGLIsCurrent(ih::Ptr{Cint})
  ccall( (:IupGLIsCurrent, iup), Cint, (Ptr{Cint},), ih)
end
function IupGLSwapBuffers(ih::Ptr{Cint})
  ccall( (:IupGLSwapBuffers, iup), Void, (Ptr{Cint},), ih)
end
function IupGLPalette(ih::Ptr{Cint}, index::Cint, r::Cfloat, g::Cfloat, b::Cfloat)
  ccall( (:IupGLPalette, iup), Void, (Ptr{Cint}, Cint, Cfloat, Cfloat, Cfloat), ih, index, r, g, b)
end
function IupGLUseFont(ih::Ptr{Cint}, first::Cint, count::Cint, list_base::Cint)
  ccall( (:IupGLUseFont, iup), Void, (Ptr{Cint}, Cint, Cint, Cint), ih, first, count, list_base)
end
function IupGLWait(gl::Cint)
  ccall( (:IupGLWait, iup), Void, (Cint,), gl)
end


text:
"""
kalman_ll(bk,y,B)

Derive the log-likelihood (up to constant) of time series y given
DSGE solution bk. B is a vector containing the indexes of the
observed variables.

Based on the MXNet function LDS_forward, which can be found here:
https://gluon.mxnet.io/chapter12_time-series/lds-scratch.html
"""

function kalman_ll(bk,y,B,σ_obs)

    Zygote.ignore() do
        global nx,nε = size(bk.D)
        global n     = nx+nε
        global Σₕ    = Diagonal([zeros(nx);ones(nε)]) # shocks
        global Bm    = I(nx+nε)[B,:] # observation equation
        global T,nobs = size(y)      # num of observations
        global Σᵥ     = σ_obs        # observation errors
    end

    A  = [[bk.C  bk.D]; zeros(nε,nx+nε)]

    # Priors (updated throughout the loop)
    μₕ  = zeros(n)
    Σₕₕ = Matrix(I(n))

    # initialise variables
    ll  = 0.0

    for t in 1:T

        Σᵥᵥ = update_Σᵥᵥ(Σₕₕ,B,Σᵥ)
        Σch = cholesky(Σᵥᵥ;check=false)
        if Σch.info != 0
            Zygote.ignore() do
                @warn "Cholesky decomposition failed. Covariance matrix not positive definite."
            end
            ll = -Inf
            break
        end
        ΣL = Σch.L

        K   = update_K(Σₕₕ,B,ΣL)
        δ   = update_δ(μₕ,y,t,B)
        μₕ  = update_μ(A,K,μₕ,δ)
        Σₕₕ = update_Σ(A,K,Σₕₕ,Bm,Σᵥ,Σₕ)
        ll  = update_ll(ll,ΣL,δ)

    end

  return ll

end

"""
# Matrix rules with pullbacks to avoid repeat calculations.
"""

"""
sumlogdiag(A::AbstractMatrix{T})

Take logs of the diagonal values of A, then sum the log values.
"""

function sumlogdiag(A::AbstractMatrix{T}) where T
    return sum(log.(diag(A)))
end

function ChainRulesCore.rrule(::typeof(sumlogdiag), A)
    Σ = sumlogdiag(A)
    function sumlogdiag_pullback(ΔΣ)
        ∂A =  Diagonal(ΔΣ ./diag(A))
        return NoTangent(), ∂A
    end
    return Σ, sumlogdiag_pullback
end

diag(F::CHOLMOD.FactorComponent) = diag(CHOLMOD.sparse(F))

function update_K(Σₕₕ,B,ΣL)
    K   = (Σₕₕ[:,B]/(ΣL'))/ΣL
    return K
end

# function ChainRulesCore.rrule(::typeof(update_K),Σₕₕ,B,ΣL)
#     A₁  = Σₕₕ[:,B]
#     A₂  = A₁/(ΣL')
#     K   = A₂/ΣL
#     function update_K_pullback(ΔK)
#         ∂A₂   = ΔK/(ΣL')
#         ∂A₁   = ∂A₂/ΣL
#         ∂ΣL   = - tril(K'*∂A₂ + ∂A₁'*K)
#         ∂Σₕₕ  = zeros(size(Σₕₕ))
#         ∂Σₕₕ[:,B]  += ∂A₁
#         return NoTangent(),∂Σₕₕ,NoTangent(),∂ΣL
#     end
#     return K, update_K_pullback
# end

function update_Σ(A,K,Σₕₕ,Bm,Σᵥ,Σₕ)
    ImKB = I - K*Bm
    X = A*(ImKB*Σₕₕ*ImKB' + sym_square(Σᵥ,K))*A' + Σₕ
    return X
end

function ChainRulesCore.rrule(::typeof(update_Σ),A,K,Σₕₕ,Bm,Σᵥ,Σₕ)
    ImKB = I - K*Bm
    B₁ = ImKB*Σₕₕ*ImKB'
    B₂ = sym_square(Σᵥ,K)
    B₃ = B₁ + B₂
    B₄ = A*B₃*A'
    X  = B₄ + Σₕ
    function update_Σ_pullback(ΔX)
        ∂B₄   = @thunk(ΔX)
        ∂B₃   = @thunk(A'*∂B₄*A)
        ∂A    = @thunk((∂B₄+∂B₄')*A*B₃)
        ∂B₂   = @thunk(∂B₃)
        ∂B₁   = @thunk(∂B₃)
        ∂Σₕₕ  = @thunk(ImKB'*∂B₁*ImKB)
        ∂ImKB = @thunk((∂B₁+∂B₁')*ImKB*Σₕₕ)
        ∂K    = @thunk(Σᵥ*(∂B₂+∂B₂')*K - ∂ImKB * Bm')
        return NoTangent(),∂A,∂K,∂Σₕₕ,NoTangent(),NoTangent(),NoTangent()
    end
    return X, update_Σ_pullback
end

function update_μ(A,K,μₕ,δ)
    μₕ_ = A*(μₕ + K*δ)
    return μₕ_
end

function ChainRulesCore.rrule(::typeof(update_μ),A,K,μₕ,δ)
    M   = (μₕ + K*δ)
    μₕ_ = A*M
    function update_μ_pullback(Δμₕ_)
        ∂A    = @thunk(Δμₕ_*M')
        ∂M    = A'*Δμₕ_
        ∂μₕ   = @thunk(∂M)
        ∂K    = @thunk(∂M*δ')
        ∂δ    = @thunk(K'*∂M)
        return (NoTangent(),∂A,∂K,∂μₕ,∂δ)
    end
    return μₕ_, update_μ_pullback
end

function update_ll(ll,Σᵥᵥ_chol,δ)
    Z   = Σᵥᵥ_chol\ δ
    ll_ = ll -0.5*Z'*Z - sumlogdiag(Σᵥᵥ_chol)
    return ll_
end

function ChainRulesCore.rrule(::typeof(update_ll),ll,Σᵥᵥ_chol,δ)
    Z   = Σᵥᵥ_chol\δ
    ll_ = ll - 0.5*dot(Z,Z) - sumlogdiag(Σᵥᵥ_chol)
    function update_ll_pullback(Δll_)
        ∂ll       =  @thunk(Δll_)
        ∂Z        = -Δll_*Z
        ∂δ        =  Σᵥᵥ_chol' \ ∂Z
        ∂Σᵥᵥ_chol =  @thunk(Diagonal(-Δll_ ./diag(Σᵥᵥ_chol)) - ∂δ *Z')
        return (NoTangent(),∂ll,∂Σᵥᵥ_chol,∂δ)
    end
    return ll_, update_ll_pullback
end

function update_δ(μₕ,ys,t,B)
    δ = ys[t,:] - μₕ[B]
    return δ
end

function ChainRulesCore.rrule(::typeof(update_δ),μₕ,ys,t,B)
    δ = ys[t,:] - μₕ[B]
    function update_δ_pullback(Δδ)
        ∂μₕ     =  zeros(size(μₕ))
        ∂μₕ[B] += -Δδ
        return (NoTangent(),∂μₕ,NoTangent(),NoTangent(),NoTangent())
    end
    return δ, update_δ_pullback
end

function update_Σᵥᵥ(Σₕₕ,B,Σᵥ)
    Σᵥᵥ  = Symmetric(Σₕₕ[B,B],:U) + Σᵥ*I
    return Σᵥᵥ
end

function ChainRulesCore.rrule(::typeof(update_Σᵥᵥ),Σₕₕ,B,Σᵥ)
    Σᵥᵥ = update_Σᵥᵥ(Σₕₕ,B,Σᵥ)
    function update_Σᵥᵥ_pullback(ΔΣᵥᵥ)
        ∂Σₕₕ       =  zeros(size(Σₕₕ))
        ∂Σₕₕ[B,B] +=  ΔΣᵥᵥ
        return (NoTangent(),∂Σₕₕ,NoTangent(),NoTangent())
    end
    return Σᵥᵥ, update_Σᵥᵥ_pullback
end

"""
sym_square(A,B)

Implements the following product:
    A*B*A'
where B is symmetric.
"""

function sym_square(A::AbstractMatrix{Float64},B::Symmetric{Float64})
    return Symmetric(A*B*A', LinearAlgebra.sym_uplo(B.uplo))
end

function sym_square(A::AbstractMatrix{Float64},B::Diagonal{Float64})
    return Symmetric(A*B*A', :U)
end

function ChainRulesCore.rrule(::typeof(sym_square),A,B)
    X = A*B*A'
    function sym_square_pullback(ΔX)
        ∂A = @thunk((ΔX+ΔX')*A*B)
        ∂B = @thunk(A'*ΔX*A)
        return NoTangent(), ∂A, ∂B
    end
    return X, sym_square_pullback
end

function sym_square(α::Float64,K::AbstractMatrix)
    return Symmetric(BLAS.syrk('U','N',α,K))
end

function ChainRulesCore.rrule(::typeof(sym_square),α::Float64,K::AbstractMatrix)
   X = sym_square(α,K)
   function sym_square_pullback(ΔX)
       ∂K = @thunk(α*(ΔX+ΔX')*K)
       return NoTangent(), NoTangent(),∂K
   end
   return X, sym_square_pullback
end


text:
"""
    @declarezero(T, args...)

To create scalar `Float64` variables `x` and `y` that equal 0, run `@declarezero Float64 x y`

taken from <https://stackoverflow.com/questions/37229630/define-multiple-variables-of-same-type-in-one-line-in-julia>

see also <https://stackoverflow.com/questions/31313040/julia-automatically-generate-functions-and-export-them>
"""
macro declarezero(T, args...)
    r = quote end
    for var in args
        # push!(r.args, :( $(esc(var))::$T = zero($T)) )
        push!(r.args, :( $(esc(var)) = zero($T)) )
    end
    r
end

"""
    @addStructFields(T, flds...)

Adds `flds` of type `T` to struct. Note that `T` can have parameters.

# Example

```julia
struct myStruct{T<:Number}
    @addStructFields Vector{T} x y
end
```
"""
macro addStructFields(T, flds::Symbol...)
    r = quote end
    for f in flds
        e = quote $(esc(f))::$(esc(T)) end
        push!(r.args, e)
    end
    return r
end

"declare `vars...` equal to `T`"
macro declareVariables(T, vars::Symbol...)
    r = quote end
    for v in vars
        e = quote $(esc(v)) = $(esc(T)) end
        push!(r.args, e)
    end
    return r
end

"Fill column `j` of `obj.f` with zeros for `f` in `flds"
macro zeroColumn_j(obj, j, flds::Symbol...)
    r = quote end
    for f in flds
        e = quote ShaleDrillingModel.zero!(
            view(
                $(esc(obj)).$(f),
                :, $(esc(j))
            )
        )
        end
        push!(r.args, e)
    end
    return r
end

"creates function f(x::T) = x.f for f in flds"
macro getFieldFunction(T, flds::Symbol...)
    r = quote end
    for f in flds
        e = esc( quote $f( x::$T ) = x.$f end )
        push!(r.args, e)
    end
    return r
end


text:
using AWS.EC2
using AWS

include("config.jl")

instances = ec2_instances_by_tag("ClusterName", ec2_clustername)
hosts = ec2_hostnames(instances)

println("Creating addprocs_external.jl")
f = open("addprocs_external.jl", "w")
machines = ec2_addprocs(instances, ec2_sshkey_file;
                        use_public_dnsname=true,
                        num_workers=ec2_num_workers,
                        workers_per_instance=ec2_workers_per_instance,
                        machines_only=true)

if ec2_install_julia
    ec2_julia_dir = "/usr/bin"
end

for mset in machines
    print(f, "addprocs([")
    lenmset = length(mset)
    for (i,m) in enumerate(mset)
        print(f, "\"$m", i == lenmset ? "\"" : "\"," )
    end
    print(f, "]; tunnel=true, dir=\"$(ec2_julia_dir)\",  sshflags=`-i $ec2_sshkey_file -o UserKnownHostsFile=/dev/null -o StrictHostKeyChecking=no -o LogLevel=ERROR`)\n\n")
end
println(f, "println(\"nworkers = \$(nworkers())\")")
flush(f)
close(f)

println()
println("Use the following commands to execute from a external node on EC2")
println()
println("  julia compute_external.jl")
println("OR")
println("  julia")
println("  julia> include(\"compute_external.jl\")")
println()


text:
push!(LOAD_PATH,joinpath(@__DIR__, ".."))
using Documenter, ICI_simulateur

makedocs(
    modules = [ICI_simulateur],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Nicolas GILET and Maxime COLOMB and ASCII team",
    sitename = "ICI_simulateur.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/icimodel/ici_simulateur.jl.git",
)


text:
module SpkIO

include("./Events.jl")
include("./IO.jl")

end # end module


text:
# benchmark broadcasting

using NumericExtensions

macro bench_bsx(name, rp, a, b1, b2, f1, f2, dims)
    quote
        println("Benchmark on $($name) ...")

        ac = copy($a)
        b1c = copy($b1)
        b2c = copy($b2)

        broadcast!($f1, ac, ac, b1c)
        ($f2)(ac, b2c, $dims)

        t1 = @elapsed for i in 1 : ($rp)
            broadcast!($f1, ac, ac, b1c)
        end

        t2 = @elapsed for i in 1 : ($rp)
            ($f2)(ac, b2c, $dims)
        end

        add_row!(perftable, $name, [t1, t2, t1 / t2])
    end
end

# benchmarks

const perftable = BenchmarkTable("Comparison of broadcasting", ["broadcast", "vbroadcast", "gain"])

# matrices

a = rand(200, 200)
b = rand(200)
b1m = reshape(b, 200, 1)
b2m = reshape(b, 1, 200)

@bench_bsx "matrix-dim(1)" 5000 a b1m b (.+) badd! 1
@bench_bsx "matrix-dim(2)" 5000 a b2m b (.+) badd! 2

# cubes

a = rand(100, 100, 4)
b1 = rand(100)
b2 = rand(100)
b3 = rand(4)
b1m = reshape(b1, 100, 1, 1)
b2m = reshape(b2, 1, 100, 1)
b3m = reshape(b3, 1, 1, 4)

@bench_bsx "cube-dim(1)" 5000 a b1m b1 (.+) badd! 1
@bench_bsx "cube-dim(2)" 5000 a b2m b2 (.+) badd! 2
@bench_bsx "cube-dim(3)" 5000 a b3m b3 (.+) badd! 3

# b12 = rand(100, 100)
# b13 = rand(100, 4)
# b23 = rand(100, 4)
# b12m = reshape(b12, 100, 100, 1)
# b13m = reshape(b13, 100, 1, 4)
# b23m = reshape(b23, 1, 100, 4)

# @bench_bsx "cube-dim(1,2)" 5000 a b12m b12 (.+) badd! (1,2)
# @bench_bsx "cube-dim(1,3)" 5000 a b13m b13 (.+) badd! (1,3)
# @bench_bsx "cube-dim(2,3)" 5000 a b23m b23 (.+) badd! (2,3)

# print results

println()
println(perftable)


text:
function primal_objective(x::Array{Float64,2}, y::Array{Float64,1}; C=1.0)
    s,d = size(x)
    f = 0
    for i in 1:s
        xi = vec(x[i,:])
        f += max(0, 1-y[i]*vecdot(w,xi))^2
    end
    f *= C
    f += 0.5*vecdot(w,w)
    return f
end

"""
`dualcd(x, y; <keyword arguments>)`

Solves the linear SVM problem on the given training set of examples `x`
and labels `y` with the Dual Coordinate Descent
method by Cho-Jui Hsieh, Kai-Wei Chang, Chih-Jen Lin, S. Sathiya Keerthi and S. Sundararajan.

See https://www.csie.ntu.edu.tw/~cjlin/papers/cddual.pdf

### Keyword arguments
* `C=1.0`:
the parameter C of the SVM problem
* `shrinking=false`: the method uses the shrinking heuristic by Thorsten Joachims if true. See https://www.cs.cornell.edu/people/tj/publications/joachims_99a.pdf
* `iterations=1000`: the number of iterations performed until `w` is returned
* `verbose=false`: the differences in the maximum and minimum projected gradient is printed if true
* `eps=0.001`: the algorithm terminates once the difference in the maximum and minimum projected gradient is smaller than `eps`. Typically `eps` has to be large for this stopping criterion to take effect.

### Returns
* weight vector `w`
"""
function dualcd(x::Array{Float64,2}, y::Array{Float64,1}; C=1.0, shrinking=false, iterations=1000, verbose=false, eps=0.001)

    s,d = size(x) # s: number of data points
    alpha = zeros(s)

    # compute initial w:
    w = At_mul_B(x, y.*alpha)

    # Precompute (Q+D)_ii:
    QD = Array{Float64,1}(s)
    for i in 1:s
        xi = vec(x[i,:])
        QD[i] = (0.5/C)
        QD[i] += vecdot(xi,xi) # doesnt have to be multiplied by y[i]*y[i] since this product is always 1
    end

    # the active set as a set of booleans.
    # if an entry is true, then the corresponding point
    # is active
    active = trues(s)
    num_active = s

    PGmax_old = Inf
    PGmin_old = -Inf

    t=0 # number of outer iterations

    while true
        t+=1

        PGmax_new = -Inf
        PGmin_new = Inf

        #TODO: permutate data in active set

        for i in 1:s
            if shrinking == false || active[i] == true
                xi = vec(x[i,:])
                G = y[i]*vecdot(w,xi) - 1.0 + (0.5/C)*alpha[i]

                # project gradient into feasible set:
                PG = 0
                if alpha[i] == 0 # exact zero check is fine here
                    if G > PGmax_old
                        active[i] = false
                        num_active-=1
                        continue
                    elseif G < 0
                        PG = G
                        #PG = min(0,G)
                    end
                else
                    PG = G
                end

                PGmax_new = max(PGmax_new, PG)
                PGmin_new = min(PGmin_new, PG)

                # perform update of alpha and w:
                if abs(PG) > 1e-6
                    alpha_old = alpha[i]
                    alpha[i] = max(0, alpha[i] - G/QD[i])
                    # w update:
                    w = w + (alpha[i] - alpha_old)*y[i]*xi
                end

            end
        end # end inner for loop

        if t >= iterations
            @printf "dualcd stopped after %d iterations.\n" t
            return w
        end

        if verbose
            println("PGmax_new - PGmin_new = ", PGmax_new - PGmin_new)
        end

        if PGmax_new - PGmin_new < eps
            if num_active == s
                @printf "dualcd stopped after %d iterations.\n" t
                return w
            else # reset active set to see if we have a solution for the full set of data points
                num_active = s
                PGmax_old = Inf
                PGmin_old = -Inf
                continue # this missing might have caused problems before
            end
        end

        PGmax_old = PGmax_new
        PGmin_old = PGmin_new

        if PGmax_old <= 0
            PGmax_old = Inf
        end

        if PGmin_old >= 0
            PGmin_old = -Inf
        end

    end # end outer while loop
end


text:
using VoronoiCells
using GeometryBasics
import Deldir
using Test

@testset "Area" begin
    Random.seed!(1)
    points = [Point2(rand(), rand()) for _ in 1:5]

    @testset "Area of Voronoi cells sum to area of window" begin
        @testset "Unit rectangle" begin
            rect = Rectangle(Point2(0, 0), Point2(1, 1))
            tess = voronoicells(points, rect)

            cell_area = voronoiarea(tess)

            @test sum(cell_area) ≈ 1
        end

        @testset "Random rectangle" begin
            rect = Rectangle(Point2(-rand(), -rand()), Point2(1 + rand(), 1 + rand()))
            tess = voronoicells(points, rect)

            cell_area = voronoiarea(tess)

            @test sum(cell_area) ≈ VoronoiCells.area(rect)
        end
    end

    @testset "Area of Voronoi cells are the same as in Deldir" begin
        rect = Rectangle(Point2(0, 0), Point2(1, 1))
        tess = voronoicells(points, rect)

        cell_area1 = VoronoiCells.voronoiarea(tess)

        x = map(getx, points)
        y = map(gety, points)
        rw = [left(rect), right(rect), lower(rect), upper(rect)]
        cell_area2 = Deldir.voronoiarea(x, y, rw)

        @test cell_area1 ≈ cell_area2
    end
end


text:
"""
    function bp_summary(prd, snp, qtl)
---
Summarize simulation results of the `b`reeding `p`rogram.
The figures and tables will be saved with an initial 'bps-'.
"""
function bp_summary(prd, snp, qtl)
    # construct pedigree
    n₀ = length(
        unique(
            Matrix(
                select(
                    filter(row ->row.g8n ==1, prd), :sir, :dam))))
    ped = begin                 # this was the way I recorded the pedigree
        n = n₀ + nrow(prd)      # total ID number, including base
        ped = zeros(Int32, n, 2) # the first `n₀` are automatically 0s
        # parents of generation 1 refer to base
        pm₁ = Matrix(select(filter(row -> row.g8n == 1, prd), :sir, :dam))
        n₁ = size(pm₁)[1]
        ped[n₀+1:n₀+n₁, :] = pm₁
        # parents of generation 2+ are about descendants
        pm₂ = Matrix(select(filter(row -> row.g8n > 1, prd), :sir, :dam))
        pm₂ .+= n₀
        ped[n₀+n₁+1:end, :] = pm₂
        ped
    end

    @info "Inbreeding using A matrix"
    begin                 # generation average inbreeding coefficients
        # only need to record generation, and one ID number of its many full sibs.
        n = size(ped)[1]
        g8n = zeros(Int, n)
        g8n[n₀+1:end] = prd.g8n
        df = DataFrame(g8n=Int[], id=Int[], cnt=Int[], ic=Float64[])
        a = b = 0
        for (g, id, pa, ma) in eachrow([g8n 1:n ped])
            g == 0 && continue
            if pa == a && ma == b
                df.cnt[end] += 1
                continue
            end
            push!(df, [g, id, 1, 0])
            a, b = pa, ma
        end
        Threads.@threads for i in 1:nrow(df)
            df.ic[i] = kinship(ped, df.id[i], df.id[i]) - 1
        end
        boxplot(df.g8n, df.ic, label="Inbreeding measured by A")
        savefig("doc/fig/inbreeding-changes-with-A.pdf")
    end

    @info "Inbreeding using homozygous loci"
    begin
        m = size(snp)[1]
        df = DataFrame(g8n = prd.g8n, ic = zeros(nrow(prd)))
        Threads.@threads for k in 1:nrow(prd)
            df.ic[k] = m - sum(snp[:, 2k-1] .⊻ snp[:, 2k])
        end
        df.ic ./= m
        t = mean(filter(row -> row.g8n == 1, df).ic)
        df.ic ./= t
        df.ic .-= 1
        boxplot(df.g8n, df.ic, label="Inbreeding measured with homozygotes ratio")
        savefig("doc/fig/inbreeding-changes-with-homozygous-loci.pdf")
    end

    # Something not right below.  Will visit this later.
    #@info "Inbreeding using G matrix"
    #begin
    #    p = mean(snp, dims = 2)
    #    tp = 2p
    #    s2pq = sum((1 .- p) .* tp)
    #    r2pq = 1 / s2pq
    #    df = DataFrame(g8n = prd.g8n, ic = zeros(nrow(prd)))
    #    for i in 1:nrow(prd)
    #        z = (snp[:, 2i-1] + snp[:, 2i]) - tp
    #        df.ic[i] = (z'z)[1]*r2pq - 1
    #    end
    #    boxplot(df.g8n, df.ic, label="Inbreeding measured with G")
    #    savefig("doc/fig/inbreeding-changes-with-G.pdf")
    #end

    gp = groupby(prd, :g8n)
    @info "Genetic variance and inbreeding over generations"
    begin                                 # plot vₐ² trends of 2 traits
        p = combine(gp, :tbv => var => :vprd)
        b = combine(gp, :t2c => var => :vbin)
        plot(p.vprd, label="Production", dpi=300)
        plot!(b.vbin, label="Binary")
        savefig("doc/fig/Va-changes.pdf")
    end
    @info join(["",
                "QTL frequency changes",
                "  - if negtive a, switch p to 1-p",
                "  - bubble plot, sizes are with respect to initial frq"], '\n')
    begin
        qg = GEAS.alleles2gt(snp[qtl[2].pos, :])
        n1 = nrow(filter(row -> row.g8n == 1, prd))
        p1 = mean(qg[:, 1:n1], dims=2)./2
        n2 = nrow(filter(row -> row.g8n == prd.g8n[end], prd))
        p2 = mean(qg[:, end-n2:end], dims=2)./2
        scatter(p1, p2 .- p1, ms=p1.*10,
                label="Frequency changes",
                xlabel="MAF",
                ylabel="frequency change",
                dpi=300)
        savefig("doc/fig/frequency-changes.pdf")
    end
end


text:
"""
    mulaipraktikum()
untuk membuka folder praktikum via IJulia.
"""
function mulaipraktikum()
    repo_directory = joinpath(@__DIR__,"..")
    path = joinpath(repo_directory, "notebookpraktikum")
    IJulia.notebook(;dir=path)
end


text:
using LinearAlgebra

# Parameters in regression model
beta = [148.73, -18.85]; v = 54.40;
y_hat(x) = x'*beta + v;
# Evaluate regression model prediction
x = [0.846, 1]; y = 115;
y_hat(x), y                     # (161.37557999999999, 115)
x = [1.324,2];  y = 234.50;
y_hat(x), y

using VMLS                      # see below to install package
D = house_sales_data();
price = D["price"]
area = D["area"]
beds = D["beds"]

v = 54.4017;
beta = [ 147.7251, -18.8534 ];
predicted = v .+ beta[1] * area + beta[2] * beds;

using Plots
scatter(price, predicted, lims = (0,800));
plot!([0, 800], [0, 800], linestyle = :dash);
# make axes equal and add labels
plot!(xlims = (0,800), ylims = (0,800), size = (500,500));
plot!(xlabel = "Actual price", ylabel = "Predicted price");

# Install VMLS to JuliaBox
# Go to git, enter 'https:////github.com/VMLS-book/VMLS.jl.git'
# Go to package manager, add 'VMLS'
# Wait up to 15 minutes (for package build)

# --OR--

# Install VMLS to Julia runtime (local machine)
# Go to Julia prompt, type ']' to open package manager
# Get library from 'pkg>' prompt
#    type 'add https:////github.com/VMLS-book/VMLS.jl'
#    'ctrl + C' to exit pakage manager prompt

max_stars_repo_name: lynnlangit/julia-linear-algebra

max_issues_repo_name: lynnlangit/julia-linear-algebra

max_forks_repo_name: lynnlangit/julia-linear-algebra


text:
abstract type ModifiedManifold{D, C} <: AbstractManifold{D, C} end

# reparametrize ========================================================================== #
struct ReparametrizedManifold{D, C, F, M<:AbstractManifold{D, C}} <: ModifiedManifold{D, C}
    base ::M
    par  ::F
end

Base.show(io::IO, rm::ReparametrizedManifold) =
    print(io, "reparametrized ($(rm.base))")

(sm::ReparametrizedManifold{D})(args::Vararg{T, D}) where {D, T} =
    sm.base(sm.par(args...)...)

reparametrized(man::AbstractManifold{D, C}, f) where {D, C} =
    ReparametrizedManifold{D, C, typeof(f), typeof(man)}(man, f)
reparametrized(man::AbstractManifold, t::Number) =
    reparametrized(man, x->x.*t)

# scaling ================================================================================ #
# TODO? split out const scaling?
struct ScaledManifold{D, C, F, M<:AbstractManifold{D, C}} <: ModifiedManifold{D, C}
    base    ::M
    scaling ::F
end

Base.show(io::IO, sm::ScaledManifold) =
    print(io, "scaled ($(sm.base))")

scaling(::AbstractManifold) = one
scaling(mm::ModifiedManifold) = scaling(mm.base)
scaling(sm::ScaledManifold) = sm.scaling
scaling(sm::ScaledManifold{<:Any, <:Any, <:Number}) = one

scaled(man::AbstractManifold{D, C}, f) where {D, C} =
    ScaledManifold{D, C, typeof(f), typeof(man)}(man, f)
Base.:*(r::Number, man::AbstractManifold) =
    scaled(man, r)
Base.:*(man::AbstractManifold, r::Number) =
    scaled(man, r)

(sm::ScaledManifold{D, <:Any, <:Number})(args::Vararg{T, D}) where {T, D} =
    sm.scaling .* sm.base(args...)
(sm::ScaledManifold{D, <:Any})(args::Vararg{T, D}) where {T, D} =
    sm.base(args...)
Base.rand(rng::AbstractRNG, sm::ScaledManifold{D, <:Any, <:Number}) where D =
    sm.scaling .* rand(rng, sm.base)
Base.rand(rng::AbstractRNG, sm::ScaledManifold) =
    rand(rng, sm.base)

# translate ============================================================================== #
struct TranslatedManifold{D, C, T<:AbstractArray,
                          M<:AbstractManifold{D, C}} <: ModifiedManifold{D, C}
    base        ::M
    translation ::T
end

Base.show(io::IO, tm::TranslatedManifold) =
    print(io, "translated ($(tm.base))")

translation(::AbstractManifold) = false
translation(mm::ModifiedManifold) = translation(mm.base)
translation(tm::TranslatedManifold) = tm.translation

Base.:+(man::M, r::T) where {M<:AbstractManifold{D, C}, T<:AbstractArray} where {D, C} =
    TranslatedManifold{D, C, T, M}(man, r)
Base.:+(r, man::AbstractManifold{D}) where D = r + man

(tm::TranslatedManifold{D})(args::Vararg{T, D}) where {T, D} =
    promotedadd(tm.base(args...), tm.translation)
Base.rand(rng::AbstractRNG, tm::TranslatedManifold) =
    promotedadd(rand(rng, tm.base), tm.translation)

# flattened ============================================================================== #
struct FlattenedManifold{D, C, M<:AbstractManifold{D, C}} <: AbstractManifold{D, C}
    base ::M
end

"""
    flattened(man)

Use to make a product flat, e.g. while `Sphere(1) × Sphere(1)` forms a torus embedded in ℝ³,
`flattened(Sphere(1)) × Sphere(1)` forms a flat torus, embedded in ℝ⁴.
"""
flattened(man::AbstractManifold{D, C}) where {D, C} =
    FlattenedManifold{D, C, typeof(man)}(man)

offsetframe(fm::FlattenedManifold{D, C}, args::Vararg{T, D}) where {D, C, T} =
    fm.base(args...), vcat(zeros(T, D+C, 1), one(T))
dimincrease(::FlattenedManifold{D, C}) where {D, C} =
    D + C

(fm::FlattenedManifold{D})(args::Vararg{T, D}) where {D, T} =
    fm.base(args...)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# ToDo: Add `density_logvalgrad!` to support HMC, etc.

@doc doc"""
    AbstractDensity

Subtypes of `AbstractDensity` must implement the function

* `BAT.density_logval`

For likelihood densities this is typically sufficient, since shape, and
bounds of parameters will be inferred from the prior.

Densities with known parameters bounds may also implement

* `BAT.param_bounds`

If the parameter bounds are unkown, but the number of parameters is known,
the function

* `BAT.nparams`

may be implemented directly (usually it is inferred from the bounds).

!!! note

    The functions `BAT.param_bounds` and `BAT.params_shape` are currently not
    part of the stable public BAT-API, their names and arguments may change
    without notice.

Densities that support named parameters should also implement

* `BAT.params_shape`
"""
abstract type AbstractDensity end
export AbstractDensity

@doc doc"""
    BAT.density_logval(density::AbstractDensity, params::Any)

Compute log of value of a multivariate density function at the given
parameter values.

Input:

* `density`: density function
* `params`: parameter values

Note: If `density_logval` is called with out-of-bounds parameters (see
`param_bounds`), the behaviour is undefined. The result for parameters that
are not within bounds is *implicitly* `-Inf`, but it is the caller's
responsibility to handle these cases.
"""
function density_logval end

@doc doc"""
    param_bounds(
        density::AbstractDensity
    )::Union{AbstractParamBounds,Missing}

Get the parameter bounds of `density`. See `density_logval` for the
implications and handling of the bounds. If the bounds are missing,
`density_logval` must be prepared to handle any parameter values.
"""
param_bounds(density::AbstractDensity) = missing

@doc doc"""
    nparams(density::AbstractDensity)::Union{Int,Missing}

Get the number of parameters of `density`. May return `missing`, if the
density supports a variable number of parameters.
"""
function nparams(density::AbstractDensity)
    bounds = param_bounds(density)
    ismissing(bounds) ? missing : nparams(bounds)
end
export nparams

@doc doc"""
    params_shape(
        density::AbstractDensity
    )::Union{ValueShapes.AbstractValueShape,Missing}

    params_shape(
        density::DistLikeDensity
    )::ValueShapes.AbstractValueShape

    params_shape(
        distribution::Distributions.Distribution
    )::ValueShapes.AbstractValueShape

*BAT-internal, not part of stable public API.*

Get the shapes of parameters of `density`.

For prior densities, the result must not be `missing`, but may be `nothing` if
the prior only supports flat parameter vectors.
"""
function params_shape end

params_shape(density::AbstractDensity) = missing

params_shape(dist::Distribution) = valshape(dist)

@doc doc"""
    eval_density_logval(
        density::AbstractDensity,
        params::AbstractVector{<:Real},
        parshapes::ValueShapes.AbstractValueShape
    )

*BAT-internal, not part of stable public API.*

Evaluates density log-value via `density_logval`.

`parshapes` *must* be compatible with `params_shape(density)`.

Checks that:

* The number of parameters of `density` (if known) matches the length of
  `params`.
* The return value of `density_logval` is not `NaN`.
* The return value of `density_logval` is less than `+Inf`.
"""
function eval_density_logval(
    density::AbstractDensity,
    params::AbstractVector{<:Real},
    parshapes::ValueShapes.AbstractValueShape
)
    npars = nparams(density)
    ismissing(npars) || (length(eachindex(params)) == npars) || throw(ArgumentError("Invalid length of parameter vector"))

    r = float(density_logval(density, _apply_parshapes(params, parshapes)))
    isnan(r) && throw(ErrorException("Return value of density_logval must not be NaN, density has type $(typeof(density))"))
    r < convert(typeof(r), +Inf) || throw(ErrorException("Return value of density_logval must not be posivite infinite, density has type $(typeof(density))"))

    r
end

_apply_parshapes(params::AbstractVector{<:Real}, parshapes::AbstractValueShape) = parshapes(params)

@doc doc"""
    DistLikeDensity <: AbstractDensity

A density that implements part of the `Distributions.Distribution` interface.
Such densities are suitable to be used as a priors.

Typically, custom priors should be implemented as subtypes of
`Distributions.Distribution`. BAT will automatically wrap them in a subtype of
`DistLikeDensity`.

Subtypes of `DistLikeDensity` are required to support more functionality
than a `AbstractDensity`, but less than a
`Distribution{Multivariate,Continuous}`.

The following functions must be implemented for subtypes:

* `BAT.density_logval`

* `BAT.param_bounds`

* `BAT.params_shape`

* `Distributions.sampler`

* `Statistics.cov`

!!! note

    The functions `BAT.param_bounds` and `BAT.params_shape` are currently not
    part of the stable public BAT-API, their names and arguments may change
    without notice.

Prior densities that support non-flat parameters should also implement

* `BAT.params_shape`

A `d::Distribution{Multivariate,Continuous}` can be converted into (wrapped
in) an `DistLikeDensity` via `conv(DistLikeDensity, d)`.
"""
abstract type DistLikeDensity <: AbstractDensity end
export DistLikeDensity

@doc doc"""
    param_bounds(density::DistLikeDensity)::AbstractParamBounds

*BAT-internal, not part of stable public API.*

Get the parameter bounds of `density`. Must not be `missing`.
"""
function param_bounds end

@doc doc"""
    nparams(density::DistLikeDensity)::Int

Get the number of parameters of prior density `density`. Must not be
`missing`, prior densities must have a fixed number of parameters. By default,
the number of parameters is inferred from the parameter bounds.
"""
nparams(density::DistLikeDensity) = nparams(param_bounds(density))


text:
# Parameterize by T so that way it can be Vector{Expression} which is defined after
struct Operation <: AbstractOperation
    op::Function
    args::Vector{Expression}
end

# Recursive ==
function Base.:(==)(x::Operation,y::Operation)
    x.op == y.op && length(x.args) == length(y.args) && all(isequal.(x.args,y.args))
end
Base.:(==)(x::Operation,y::Number) = false
Base.:(==)(x::Number,y::Operation) = false
Base.:(==)(x::Operation,y::Void) = false
Base.:(==)(x::Void,y::Operation) = false
Base.:(==)(x::Variable,y::Operation) = y == Operation(identity,Expression[x])
Base.:(==)(x::Operation,y::Variable) = x == Operation(identity,Expression[y])

# Don't recurse inversion for Jacobians
Base.inv(x::Operation) = x

function Base.Expr(O::Operation)
    Expr(:call, Symbol(O.op), Expr.(O.args)...)
end

Base.show(io::IO,O::Operation) = print(io,string(Expr(O)))

# Bigger printing
# Is there a way to just not have this as the default?
function Base.parse(::Type{Operation},ex::Expr)
    f = ex.args[1]
    operands = ex.args[2:end]
    if ex.head == :call && any(x -> x isa Expr, ex.args)
        args = Expression[parse(Operation,o) for o in operands]
        parse(Operation, f, args)
    else
        parse(Operation, f, Expression[parse(Operation,o) for o in operands])
    end
end
Base.parse(::Type{Operation},x::Expression) = x
Base.parse(::Type{Operation},sym::Symbol,args) = Operation(eval(sym), args)
Base.parse(::Type{Operation},x::Union{Symbol, Number}) = x

"""
find_replace(O::Operation,x::Variable,y::Expression)

Finds the variable `x` in Operation `O` and replaces it with the Expression `y`
"""
function find_replace!(O::Operation,x::Variable,y::Expression)
    for i in eachindex(O.args)
        if isequal(O.args[i],x)
            O.args[i] = y
        elseif typeof(O.args[i]) <: Operation
            find_replace!(O.args[i],x,y)
        end
    end
end

# For inv
Base.convert(::Type{Operation},x::Int) = Operation(identity,Expression[Constant(x)])
Base.convert(::Type{Operation},x::Variable) = Operation(identity,Expression[x])


text:
# domain.jl

create_key(name::String) = string(hash(name * string(time())))

struct Organization
    id::String
    name::String
    address::String
    zip::String
    city::String
    country::String
    # constructors
    Organization(name, address, zip, city, country) = new(create_key(name), name, address, zip, city, country)
end

struct Training
    id::String
    date::DateTime
    days::Int64
    name::String
    price::Float64
end

struct Order
    id::String
    org::Organization
    training::Training
    order_ref::String
    contact_name::String
    contact_email::String
    students::Array{String,1}
    # constructors
    Order(org, training, order_ref, contact_name, contact_email, students) =
        new(create_key(contact_name), org, training, order_ref, contact_name, contact_email, students)
end


text:
@testset "4.6.1.4 (d cot)^n (a+b csc)^m" begin
    (a, b, x, ) = @variables a b x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(e*tan(c+d*x))^m*(a+a*csc(c+d*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*tan(c+d*x))^m*(a+a*csc(c+d*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tan(c+d*x)^m*(a+a*csc(c+d*x))^n=#

    @test_int [tan(x)^4/(a + a*csc(x)), x, 5, x/a - ((15 - 8*csc(x))*tan(x))/(15*a) + ((5 - 4*csc(x))*tan(x)^3)/(15*a) - ((1 - csc(x))*tan(x)^5)/(5*a)]
    @test_int [tan(x)^3/(a + a*csc(x)), x, 3, (5*log(1 - sin(x)))/(16*a) + (11*log(1 + sin(x)))/(16*a) + 1/(8*a*(1 - sin(x))) - 1/(8*a*(1 + sin(x))^2) + 3/(4*a*(1 + sin(x)))]
    @test_int [tan(x)^2/(a + a*csc(x)), x, 4, -(x/a) + ((3 - 2*csc(x))*tan(x))/(3*a) - ((1 - csc(x))*tan(x)^3)/(3*a)]
    @test_int [tan(x)^1/(a + a*csc(x)), x, 3, -(log(1 - sin(x))/(4*a)) - (3*log(1 + sin(x)))/(4*a) - 1/(2*a*(1 + sin(x)))]
    @test_int [cot(x)^1/(a + a*csc(x)), x, 2, log(1 + sin(x))/a]
    @test_int [cot(x)^2/(a + a*csc(x)), x, 3, -(x/a) - atanh(cos(x))/a]
    @test_int [cot(x)^3/(a + a*csc(x)), x, 3, -(csc(x)/a) - log(sin(x))/a]
    @test_int [cot(x)^4/(a + a*csc(x)), x, 4, x/a + atanh(cos(x))/(2*a) + (cot(x)*(2 - csc(x)))/(2*a)]
    @test_int [cot(x)^5/(a + a*csc(x)), x, 3, csc(x)/a + csc(x)^2/(2*a) - csc(x)^3/(3*a) + log(sin(x))/a]
    @test_int [cot(x)^6/(a + a*csc(x)), x, 5, -(x/a) - (3*atanh(cos(x)))/(8*a) + (cot(x)^3*(4 - 3*csc(x)))/(12*a) - (cot(x)*(8 - 3*csc(x)))/(8*a)]
    @test_int [cot(x)^7/(a + a*csc(x)), x, 3, -(csc(x)/a) - csc(x)^2/a + (2*csc(x)^3)/(3*a) + csc(x)^4/(4*a) - csc(x)^5/(5*a) - log(sin(x))/a]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*tan(c+d*x)^m*(a+a*csc(c+d*x))^(n/2)=#

    #= ::Title:: =#
    #=Integrands*of*the*form*(e*tan(c+d*x))^m*(a+b*csc(c+d*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*tan(c+d*x))^m*(a+b*csc(c+d*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tan(c+d*x)^m*(a+b*csc(c+d*x))^n=#

    @test_int [tan(x)^5/(a + b*csc(x)), x, 3, 1/(16*(a + b)*(1 - csc(x))^2) + (5*a + 7*b)/(16*(a + b)^2*(1 - csc(x))) + 1/(16*(a - b)*(1 + csc(x))^2) + (5*a - 7*b)/(16*(a - b)^2*(1 + csc(x))) - ((8*a^2 + 21*a*b + 15*b^2)*log(1 - csc(x)))/(16*(a + b)^3) - ((8*a^2 - 21*a*b + 15*b^2)*log(1 + csc(x)))/(16*(a - b)^3) + (b^6*log(a + b*csc(x)))/(a*(a^2 - b^2)^3) - log(sin(x))/a]
    @test_int [tan(x)^3/(a + b*csc(x)), x, 3, -(1/(4*(a + b)*(1 - csc(x)))) - 1/(4*(a - b)*(1 + csc(x))) + ((2*a + 3*b)*log(1 - csc(x)))/(4*(a + b)^2) + ((2*a - 3*b)*log(1 + csc(x)))/(4*(a - b)^2) + (b^4*log(a + b*csc(x)))/(a*(a^2 - b^2)^2) + log(sin(x))/a]
    @test_int [tan(x)^1/(a + b*csc(x)), x, 3, -(log(1 - csc(x))/(2*(a + b))) - log(1 + csc(x))/(2*(a - b)) + (b^2*log(a + b*csc(x)))/(a*(a^2 - b^2)) - log(sin(x))/a]
    @test_int [cot(x)^1/(a + b*csc(x)), x, 4, log(a + b*csc(x))/a + log(sin(x))/a]
    @test_int [cot(x)^3/(a + b*csc(x)), x, 3, -(csc(x)/b) - ((1 - a^2/b^2)*log(a + b*csc(x)))/a - log(sin(x))/a]
    @test_int [cot(x)^5/(a + b*csc(x)), x, 3, -(((a^2 - 2*b^2)*csc(x))/b^3) + (a*csc(x)^2)/(2*b^2) - csc(x)^3/(3*b) + ((a^2 - b^2)^2*log(a + b*csc(x)))/(a*b^4) + log(sin(x))/a]
    @test_int [cot(x)^7/(a + b*csc(x)), x, 3, -(((a^4 - 3*a^2*b^2 + 3*b^4)*csc(x))/b^5) + (a*(a^2 - 3*b^2)*csc(x)^2)/(2*b^4) - ((a^2 - 3*b^2)*csc(x)^3)/(3*b^3) + (a*csc(x)^4)/(4*b^2) - csc(x)^5/(5*b) + ((a^2 - b^2)^3*log(a + b*csc(x)))/(a*b^6) - log(sin(x))/a]

    @test_int [tan(x)^4/(a + b*csc(x)), x, 16, x/a + (2*b^5*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(a*(a^2 - b^2)^(5/2)) - (b^3*sec(x))/(a^2 - b^2)^2 + (b*sec(x))/(a^2 - b^2) - (b*sec(x)^3)/(3*(a^2 - b^2)) + (a*b^2*tan(x))/(a^2 - b^2)^2 - (a*tan(x))/(a^2 - b^2) + (a*tan(x)^3)/(3*(a^2 - b^2)), -((a*b^2*x)/(a^2 - b^2)^2) + (b^4*x)/(a*(a^2 - b^2)^2) + (a*x)/(a^2 - b^2) + (2*b^5*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(a*(a^2 - b^2)^(5/2)) - (b^3*sec(x))/(a^2 - b^2)^2 + (b*sec(x))/(a^2 - b^2) - (b*sec(x)^3)/(3*(a^2 - b^2)) + (a*b^2*tan(x))/(a^2 - b^2)^2 - (a*tan(x))/(a^2 - b^2) + (a*tan(x)^3)/(3*(a^2 - b^2))]
    @test_int [tan(x)^2/(a + b*csc(x)), x, 10, -(x/a) + (2*b^3*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(a*(a^2 - b^2)^(3/2)) - (b*sec(x))/(a^2 - b^2) + (a*tan(x))/(a^2 - b^2), -((a*x)/(a^2 - b^2)) + (b^2*x)/(a*(a^2 - b^2)) + (2*b^3*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(a*(a^2 - b^2)^(3/2)) - (b*sec(x))/(a^2 - b^2) + (a*tan(x))/(a^2 - b^2)]
    @test_int [cot(x)^2/(a + b*csc(x)), x, 8, -(x/a) - atanh(cos(x))/b + (2*sqrt(a^2 - b^2)*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(a*b)]
    @test_int [cot(x)^4/(a + b*csc(x)), x, 7, x/a - ((2*a^2 - 3*b^2)*atanh(cos(x)))/(2*b^3) + (2*(a^2 - b^2)^(3/2)*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(a*b^3) + (a*cot(x))/b^2 - (cot(x)*csc(x))/(2*b)]
    @test_int [cot(x)^6/(a + b*csc(x)), x, 16, -(x/a) - (3*atanh(cos(x)))/(8*b) - ((a^2 - 3*b^2)*atanh(cos(x)))/(2*b^3) - ((a^4 - 3*a^2*b^2 + 3*b^4)*atanh(cos(x)))/b^5 + (2*(a^2 - b^2)^(5/2)*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(a*b^5) + (a*cot(x))/b^2 + (a*(a^2 - 3*b^2)*cot(x))/b^4 + (a*cot(x)^3)/(3*b^2) - (3*cot(x)*csc(x))/(8*b) - ((a^2 - 3*b^2)*cot(x)*csc(x))/(2*b^3) - (cot(x)*csc(x)^3)/(4*b)]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*tan(c+d*x)^m*(a+b*csc(c+d*x))^(n/2)=#
    end

max_stars_repo_path: test/4 Trig functions/4.6 Cosecant/4.6.1.4 (d cot)^n (a+b csc)^m.jl


text:
module TestPerformanceSpecializations

using Test
using Trixi

include("test_trixi.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "Performance specializations" begin
  @timed_testset "TreeMesh2D, flux_shima_etal_turbo" begin
    trixi_include(@__MODULE__,
      joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_ec.jl"),
      initial_refinement_level=0, tspan=(0.0, 0.0), polydeg=3,
      volume_flux=flux_shima_etal_turbo, surface_flux=flux_shima_etal_turbo)
    u_ode = copy(sol.u[end])
    du_ode = zero(u_ode)

    # Preserve original memory since it will be `unsafe_wrap`ped and might
    # thus otherwise be garbage collected
    GC.@preserve u_ode du_ode begin
      u = Trixi.wrap_array(u_ode, semi)
      du = Trixi.wrap_array(du_ode, semi)
      nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations)

      # Call the optimized default version
      du .= 0
      Trixi.split_form_kernel!(
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_specialized = du[:, :, :, :, 1]

      # Call the plain version - note the argument type `Function` of
      # `semi.solver.volume_integral.volume_flux`
      du .= 0
      invoke(Trixi.split_form_kernel!,
        Tuple{typeof(du), typeof(u), Integer, typeof(semi.mesh),
        typeof(nonconservative_terms), typeof(semi.equations),
        Function, typeof(semi.solver), typeof(semi.cache), Bool},
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_baseline = du[:, :, :, :, 1]

      @test du_specialized ≈ du_baseline
    end
  end

  @timed_testset "TreeMesh2D, flux_ranocha_turbo" begin
    trixi_include(@__MODULE__,
      joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_ec.jl"),
      initial_refinement_level=0, tspan=(0.0, 0.0), polydeg=3,
      volume_flux=flux_ranocha_turbo, surface_flux=flux_ranocha_turbo)
    u_ode = copy(sol.u[end])
    du_ode = zero(u_ode)

    # Preserve original memory since it will be `unsafe_wrap`ped and might
    # thus otherwise be garbage collected
    GC.@preserve u_ode du_ode begin
      u = Trixi.wrap_array(u_ode, semi)
      du = Trixi.wrap_array(du_ode, semi)
      nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations)

      # Call the optimized default version
      du .= 0
      Trixi.split_form_kernel!(
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_specialized = du[:, :, :, :, 1]

      # Call the plain version - note the argument type `Function` of
      # `semi.solver.volume_integral.volume_flux`
      du .= 0
      invoke(Trixi.split_form_kernel!,
        Tuple{typeof(du), typeof(u), Integer, typeof(semi.mesh),
        typeof(nonconservative_terms), typeof(semi.equations),
        Function, typeof(semi.solver), typeof(semi.cache), Bool},
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_baseline = du[:, :, :, :, 1]

      @test du_specialized ≈ du_baseline
    end
  end

  @timed_testset "StructuredMesh2D, flux_shima_etal_turbo" begin
    trixi_include(@__MODULE__,
      joinpath(examples_dir(), "structured_2d_dgsem", "elixir_euler_ec.jl"),
      cells_per_dimension=(1, 1), tspan=(0.0, 0.0), polydeg=3,
      volume_flux=flux_shima_etal_turbo, surface_flux=flux_shima_etal_turbo)
    u_ode = copy(sol.u[end])
    du_ode = zero(u_ode)

    # Preserve original memory since it will be `unsafe_wrap`ped and might
    # thus otherwise be garbage collected
    GC.@preserve u_ode du_ode begin
      u = Trixi.wrap_array(u_ode, semi)
      du = Trixi.wrap_array(du_ode, semi)
      nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations)

      # Call the optimized default version
      du .= 0
      Trixi.split_form_kernel!(
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_specialized = du[:, :, :, :, 1]

      # Call the plain version - note the argument type `Function` of
      # `semi.solver.volume_integral.volume_flux`
      du .= 0
      invoke(Trixi.split_form_kernel!,
        Tuple{typeof(du), typeof(u), Integer, typeof(semi.mesh),
        typeof(nonconservative_terms), typeof(semi.equations),
        Function, typeof(semi.solver), typeof(semi.cache), Bool},
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_baseline = du[:, :, :, :, 1]

      @test du_specialized ≈ du_baseline
    end
  end

  @timed_testset "StructuredMesh2D, flux_ranocha_turbo" begin
    trixi_include(@__MODULE__,
      joinpath(examples_dir(), "structured_2d_dgsem", "elixir_euler_ec.jl"),
      cells_per_dimension=(1, 1), tspan=(0.0, 0.0), polydeg=3,
      volume_flux=flux_ranocha_turbo, surface_flux=flux_ranocha_turbo)
    u_ode = copy(sol.u[end])
    du_ode = zero(u_ode)

    # Preserve original memory since it will be `unsafe_wrap`ped and might
    # thus otherwise be garbage collected
    GC.@preserve u_ode du_ode begin
      u = Trixi.wrap_array(u_ode, semi)
      du = Trixi.wrap_array(du_ode, semi)
      nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations)

      # Call the optimized default version
      du .= 0
      Trixi.split_form_kernel!(
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_specialized = du[:, :, :, :, 1]

      # Call the plain version - note the argument type `Function` of
      # `semi.solver.volume_integral.volume_flux`
      du .= 0
      invoke(Trixi.split_form_kernel!,
        Tuple{typeof(du), typeof(u), Integer, typeof(semi.mesh),
        typeof(nonconservative_terms), typeof(semi.equations),
        Function, typeof(semi.solver), typeof(semi.cache), Bool},
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_baseline = du[:, :, :, :, 1]

      @test du_specialized ≈ du_baseline
    end
  end

  @timed_testset "TreeMesh3D, flux_shima_etal_turbo" begin
    trixi_include(@__MODULE__,
      joinpath(examples_dir(), "tree_3d_dgsem", "elixir_euler_ec.jl"),
      initial_refinement_level=0, tspan=(0.0, 0.0), polydeg=3,
      volume_flux=flux_shima_etal_turbo, surface_flux=flux_shima_etal_turbo)
    u_ode = copy(sol.u[end])
    du_ode = zero(u_ode)

    # Preserve original memory since it will be `unsafe_wrap`ped and might
    # thus otherwise be garbage collected
    GC.@preserve u_ode du_ode begin
      u = Trixi.wrap_array(u_ode, semi)
      du = Trixi.wrap_array(du_ode, semi)
      nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations)

      # Call the optimized default version
      du .= 0
      Trixi.split_form_kernel!(
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_specialized = du[:, :, :, :, 1]

      # Call the plain version - note the argument type `Function` of
      # `semi.solver.volume_integral.volume_flux`
      du .= 0
      invoke(Trixi.split_form_kernel!,
        Tuple{typeof(du), typeof(u), Integer, typeof(semi.mesh),
        typeof(nonconservative_terms), typeof(semi.equations),
        Function, typeof(semi.solver), typeof(semi.cache), Bool},
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_baseline = du[:, :, :, :, 1]

      @test du_specialized ≈ du_baseline
    end
  end

  @timed_testset "TreeMesh3D, flux_ranocha_turbo" begin
    trixi_include(@__MODULE__,
      joinpath(examples_dir(), "tree_3d_dgsem", "elixir_euler_ec.jl"),
      initial_refinement_level=0, tspan=(0.0, 0.0), polydeg=3,
      volume_flux=flux_ranocha_turbo, surface_flux=flux_ranocha_turbo)
    u_ode = copy(sol.u[end])
    du_ode = zero(u_ode)

    # Preserve original memory since it will be `unsafe_wrap`ped and might
    # thus otherwise be garbage collected
    GC.@preserve u_ode du_ode begin
      u = Trixi.wrap_array(u_ode, semi)
      du = Trixi.wrap_array(du_ode, semi)
      nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations)

      # Call the optimized default version
      du .= 0
      Trixi.split_form_kernel!(
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_specialized = du[:, :, :, :, 1]

      # Call the plain version - note the argument type `Function` of
      # `semi.solver.volume_integral.volume_flux`
      du .= 0
      invoke(Trixi.split_form_kernel!,
        Tuple{typeof(du), typeof(u), Integer, typeof(semi.mesh),
        typeof(nonconservative_terms), typeof(semi.equations),
        Function, typeof(semi.solver), typeof(semi.cache), Bool},
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_baseline = du[:, :, :, :, 1]

      @test du_specialized ≈ du_baseline
    end
  end

  @timed_testset "StructuredMesh3D, flux_shima_etal_turbo" begin
    trixi_include(@__MODULE__,
      joinpath(examples_dir(), "structured_3d_dgsem", "elixir_euler_ec.jl"),
      cells_per_dimension=(1, 1, 1), tspan=(0.0, 0.0), polydeg=3,
      volume_flux=flux_shima_etal_turbo, surface_flux=flux_shima_etal_turbo)
    u_ode = copy(sol.u[end])
    du_ode = zero(u_ode)

    # Preserve original memory since it will be `unsafe_wrap`ped and might
    # thus otherwise be garbage collected
    GC.@preserve u_ode du_ode begin
      u = Trixi.wrap_array(u_ode, semi)
      du = Trixi.wrap_array(du_ode, semi)
      nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations)

      # Call the optimized default version
      du .= 0
      Trixi.split_form_kernel!(
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_specialized = du[:, :, :, :, 1]

      # Call the plain version - note the argument type `Function` of
      # `semi.solver.volume_integral.volume_flux`
      du .= 0
      invoke(Trixi.split_form_kernel!,
        Tuple{typeof(du), typeof(u), Integer, typeof(semi.mesh),
        typeof(nonconservative_terms), typeof(semi.equations),
        Function, typeof(semi.solver), typeof(semi.cache), Bool},
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_baseline = du[:, :, :, :, 1]

      @test du_specialized ≈ du_baseline
    end
  end

  @timed_testset "StructuredMesh3D, flux_ranocha_turbo" begin
    trixi_include(@__MODULE__,
      joinpath(examples_dir(), "structured_3d_dgsem", "elixir_euler_ec.jl"),
      cells_per_dimension=(1, 1, 1), tspan=(0.0, 0.0), polydeg=3,
      volume_flux=flux_ranocha_turbo, surface_flux=flux_ranocha_turbo)
    u_ode = copy(sol.u[end])
    du_ode = zero(u_ode)

    # Preserve original memory since it will be `unsafe_wrap`ped and might
    # thus otherwise be garbage collected
    GC.@preserve u_ode du_ode begin
      u = Trixi.wrap_array(u_ode, semi)
      du = Trixi.wrap_array(du_ode, semi)
      nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations)

      # Call the optimized default version
      du .= 0
      Trixi.split_form_kernel!(
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_specialized = du[:, :, :, :, 1]

      # Call the plain version - note the argument type `Function` of
      # `semi.solver.volume_integral.volume_flux`
      du .= 0
      invoke(Trixi.split_form_kernel!,
        Tuple{typeof(du), typeof(u), Integer, typeof(semi.mesh),
        typeof(nonconservative_terms), typeof(semi.equations),
        Function, typeof(semi.solver), typeof(semi.cache), Bool},
        du, u, 1, semi.mesh,
        nonconservative_terms, semi.equations,
        semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
      du_baseline = du[:, :, :, :, 1]

      @test du_specialized ≈ du_baseline
    end
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end #module


text:
# This file is the successor to "images/GenerateImages.jl" (v0.5.1).
module Colormaps

using Colors

struct ColormapSVG <: Main.SVG
    buf::IOBuffer
end

function ColormapSVG(cm::AbstractVector{T}, w="96mm", h="4mm") where T<:Color
    io = IOBuffer()
    n = length(cm)
    write(io,
        """
        <svg xmlns="http://www.w3.org/2000/svg" version="1.1"
             width="$w" height="$h"
             viewBox="0 0 $n 1" preserveAspectRatio="none"
             shape-rendering="crispEdges" stroke="none">
        """)
    for i in 1:n
        c = hex(cm[i])
        write(io,
            """
            <rect width="$(n+1-i)" height="1" x="$(i-1)" y="0" fill="#$c" />
            """)
    end
    write(io, "</svg>")
    ColormapSVG(io)
end

end


text:
include("ProductDomain.jl")

## boundary

boundary(d::ProductDomain{Tuple{A,B}}) where {A<:IntervalOrSegment,B<:IntervalOrSegment} =
    PiecewiseSegment([Vec(leftendpoint(factor(d,1)),leftendpoint(factor(d,2))),
                      Vec(rightendpoint(factor(d,1)),leftendpoint(factor(d,2))),
                      Vec(rightendpoint(factor(d,1)),rightendpoint(factor(d,2))),
                      Vec(leftendpoint(factor(d,1)),rightendpoint(factor(d,2))),
                      Vec(leftendpoint(factor(d,1)),leftendpoint(factor(d,2)))])
boundary(d::ProductDomain{Tuple{A,B}}) where {A<:IntervalOrSegment,B<:PeriodicSegment} =
    UnionDomain((PeriodicSegment(Vec(rightendpoint(factor(d,1)),leftendpoint(factor(d,2))),Vec(rightendpoint(factor(d,1)),rightendpoint(factor(d,2)))),
        PeriodicSegment(Vec(leftendpoint(factor(d,1)),rightendpoint(factor(d,2))),Vec(leftendpoint(factor(d,1)),leftendpoint(factor(d,2))))))
boundary(d::ProductDomain{Tuple{A,B}}) where {A<:PeriodicSegment,B<:IntervalOrSegment} =
    UnionDomain((PeriodicSegment(Vec(leftendpoint(factor(d,1)),leftendpoint(factor(d,2))),Vec(rightendpoint(factor(d,1)),leftendpoint(factor(d,2)))),
        PeriodicSegment(Vec(rightendpoint(factor(d,1)),rightendpoint(factor(d,2))),Vec(leftendpoint(factor(d,1)),rightendpoint(factor(d,2))))))
boundary(d::ProductDomain{Tuple{A,B}}) where {A<:PeriodicSegment,B<:PeriodicSegment} = EmptyDomain()

## Union
function Base.join(p1::AbstractVector{IT},p2::AbstractVector{IT}) where IT<:IntervalOrSegment
    for k=length(p1):-1:1,j=length(p2):-1:1
        if p1[k]==reverse(p2[j])
            deleteat!(p1,k)
            deleteat!(p2,j)
            break
        end
    end
    [p1;p2]
end

function boundary(d::UnionDomain{Tuple{PD1,PD2}}) where {PD1<:ProductDomain,PD2<:ProductDomain}
    bnd=map(d->components(∂(d)),d.domains)
    PiecewiseSegment(reduce(join,bnd))
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 = "nanoflann"
version = v"1.3.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/jlblancoc/nanoflann/archive/refs/tags/v$(version).tar.gz", "e100b5fc8d72e9426a80312d852a62c05ddefd23f17cbb22ccd8b458b11d0bea")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/nanoflann-*/
mkdir build && cd build

cmake .. \
-DCMAKE_INSTALL_PREFIX=${prefix} \
-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
-DCMAKE_BUILD_TYPE=Release \
-DBUILD_EXAMPLES=OFF \
-DBUILD_BENCHMARKS=OFF \
-DBUILD_TESTS=OFF

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 = [AnyPlatform()]

# The products that we will ensure are always built
#nanoflann is a header only library, so no binary products
products = Product[]

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

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#TODO: Apply actions according to load type
function _get_load_data(sys::System, b::Bus)
    activepower = 0.0
    reactivepower = 0.0
    for l in get_components(ElectricLoad, sys)
        if !isa(l,FixedAdmittance) && (l.bus == b)
            activepower += get_activepower(l)
            reactivepower += get_reactivepower(l)
        end
    end
    return activepower, reactivepower
end

"""
    make_pf(sys)

Create the objects needed to solve an powerflow case using NLsolve.jl solvers. Returns
an anonymous function with the powerflow equations, initial conditions and a dict to link the
solutions to the original system. Only supports systems with a single generator per bus and
currently doesn't support distributed slack buses and doesn't enforce reactive power limits.

## Example
```julia
pf!, x0 = make_pf(sys)
res = NLsolve.nlsolve(pf!, x0)
```

# Arguments
    * `sys`::System : a PowerSystems.jl system

"""
function make_pf(system)
    buses = sort(collect(get_components(Bus, system)), by = x -> get_number(x))
    N_BUS = length(buses)

    # assumes the ordering in YBus is the same as in the buses.
    Yb = Ybus(system)
    x0 = zeros(N_BUS * 2)

    # Use vectors to cache data for closure
    # These should be read only
    P_GEN_BUS  = fill(0.0, N_BUS)
    Q_GEN_BUS  = fill(0.0, N_BUS)
    P_LOAD_BUS = fill(0.0, N_BUS)
    Q_LOAD_BUS = fill(0.0, N_BUS)

    BUSES = (b for b in enumerate(buses))

    state_variable_count = 1

    for (ix, b) in BUSES
        bus_number = get_number(b)::Int
        bus_angle = get_angle(b)::Float64
        bus_voltage = get_voltage(b)::Float64
        generator = nothing
        for gen in get_components(Generator, system)
            if gen.bus == b
                !isnothing(generator) && throw(DataFormatError("There is more than one generator connected to Bus $b.name"))
                generator = gen
            end
        end
        P_GEN_BUS[ix] = isnothing(generator) ? 0.0 : get_activepower(generator)
        Q_GEN_BUS[ix] = isnothing(generator) ? 0.0 : get_reactivepower(generator)
        P_LOAD_BUS[ix], Q_LOAD_BUS[ix] = _get_load_data(system, b)

        if b.bustype == REF::BusType
            x0[state_variable_count]     = P_GEN_BUS[ix]
            x0[state_variable_count + 1] = Q_GEN_BUS[ix]
            state_variable_count += 2
        elseif b.bustype == PV::BusType
            x0[state_variable_count]     = Q_GEN_BUS[ix]
            x0[state_variable_count + 1] = bus_angle
            state_variable_count += 2
        elseif b.bustype == PQ::BusType
            x0[state_variable_count]     = bus_voltage
            x0[state_variable_count + 1] = bus_angle
            state_variable_count += 2
        end
    end

    @assert state_variable_count - 1 == N_BUS * 2

    function pf!(F, X)

        P_net = Vector{Float64}(undef, N_BUS)
        Q_net = Vector{Float64}(undef, N_BUS)
        Vm = Vector{Float64}(undef, N_BUS)
        θ = Vector{Float64}(undef, N_BUS)

        for (ix, b) in BUSES
            bus_voltage = get_voltage(b)
            bus_angle = get_angle(b)
            if b.bustype == REF::BusType
                # When bustype == REFERENCE Bus, state variables are Active and Reactive Power Generated
                P_net[ix] = X[2 * ix - 1] - P_LOAD_BUS[ix]
                Q_net[ix] = X[2 * ix] - Q_LOAD_BUS[ix]
                Vm[ix]    = bus_voltage
                θ[ix]     = bus_angle
            elseif b.bustype == PV::BusType
                # When bustype == PV Bus, state variables are Reactive Power Generated and Voltage Angle
                P_net[ix] = P_GEN_BUS[ix] - P_LOAD_BUS[ix]
                Q_net[ix] = X[2 * ix - 1] - Q_LOAD_BUS[ix]
                Vm[ix]    = bus_voltage
                θ[ix]     = X[2 * ix]
            elseif b.bustype == PQ::BusType
                # When bustype == PQ Bus, state variables are Voltage Magnitude and Voltage Angle
                P_net[ix] = P_GEN_BUS[ix] - P_LOAD_BUS[ix]
                Q_net[ix] = Q_GEN_BUS[ix] - Q_LOAD_BUS[ix]
                Vm[ix]    = X[2 * ix - 1]
                θ[ix]     = X[2 * ix]
            end
        end

        # F is active and reactive power balance equations at all buses
        state_count = 1
        for (ix_f, bf) in BUSES
            S = -P_net[ix_f] + -Q_net[ix_f]im
            V_f = Vm[ix_f] * ( cos(θ[ix_f]) + sin(θ[ix_f])im )
            for (ix_t, bt) in BUSES
                iszero(Yb[ix_f, ix_t]::ComplexF64) && continue
                V_t = Vm[ix_t] * ( cos(θ[ix_t]) + sin(θ[ix_t])im )
                S += V_f * conj(V_t) * conj(Yb[ix_f, ix_t]::ComplexF64)
            end
            F[state_count] = real(S)
            F[state_count + 1] = imag(S)
            state_count += 2
        end

        return F

    end

    return ( pf!, x0 )
end


text:
# Solving LP with AFBA
#
#   minimize    c'x             -> f = <c,.>
#   subject to  Ax = b          -> h = ind{b}, L = A
#               x >= 0          -> g = ind{>=0}
#
# Dual LP
#
#   maximize    b'y
#   subject to  A'y <= c
#
# Optimality conditions
#
#   x >= 0              [primal feasibility 1]
#   Ax = b              [primal feasibility 2]
#   A'y <= c            [dual feasibility]
#   X'(c - A'y) = 0     [complementarity slackness]
#

Random.seed!(0)

n = 100 # primal dimension
m = 80 # dual dimension (i.e. number of linear equalities)
k = 50 # number of active dual constraints (must be 0 <= k <= n)
x_star = vcat(rand(k), zeros(n-k)) # primal optimal point
s_star = vcat(zeros(k), rand(n-k)) # dual optimal slack variable
y_star = randn(m) # dual optimal point

A = randn(m, n)
b = A*x_star
c = A'*y_star + s_star

using ProximalOperators

f = Linear(c)
g = IndNonnegative()
h = IndPoint(b)

using ProximalAlgorithms

x0 = zeros(n)
y0 = zeros(m)

it_afba, x_afba, y_afba, solver_afba = ProximalAlgorithms.AFBA(x0, y0; f=f, g=g, h=h, L=A, tol=1e-8, maxit=10000)

# Check and print solution quality measures
# (for some reason the returned dual iterate is the negative of the dual LP variable y)

TOL_ASSERT = 1e-6

meas_nonneg = -minimum(min.(0.0, x_afba))
println("Nonnegativity          : ", meas_nonneg)
@test meas_nonneg <= TOL_ASSERT

meas_feas_primal = norm(A*x_afba - b)
println("Primal feasibility     : ", meas_feas_primal)
@test meas_feas_primal <= TOL_ASSERT

meas_feas_dual = maximum(max.(0.0, -A'*y_afba - c))
println("Dual feasibility       : ", meas_feas_dual)
@test meas_feas_dual <= TOL_ASSERT

meas_compl = abs(dot(c + A'*y_afba, x_afba))
println("Complementarity        : ", meas_compl)
@test meas_compl <= TOL_ASSERT

println("Primal objective       : ", dot(c, x_afba))
println("Dual objective         : ", dot(b, -y_afba))


text:
#-----------------------------------------------------------------------# Kahan Sum
"""
    KahanSum(T::Type = Float64)

Track the overall sum. Includes a compensation term that effectively doubles
precision, see
[Wikipedia](https://en.wikipedia.org/wiki/Kahan_summation_algorithm) for
details.

# Example

    fit!(KahanSum(Float64), fill(1, 100))
"""
mutable struct KahanSum{T<:Number} <: OnlineStat{Number}
    sum::T
    c::T
    n::Int
end
KahanSum(T::Type = Float64) = KahanSum(T(0), T(0), 0)
Base.sum(o::KahanSum) = o.sum
function _fit!(o::KahanSum{T}, x::Number) where {T<:Number}
    y = convert(T, x) - o.c
    t = o.sum + y
    o.c = (t - o.sum) - y
    o.sum = t
    o.n += 1
end

function _merge!(o::T, o2::T) where {T <: KahanSum}
    # correct both sums symmetrically if o and o2 are of different magnitude,
    # the correction factor of the smaller one will be lost.
    y = o2.sum - o.c
    y2 = o.sum - o2.c

    t = y2 + y

    # Here we find the case improved by Neumaier upon Kahan's algorithm, i.e. the
    # number added is orders of magnitude larger than the sum.
    if abs(y) < abs(y2)
        o.c = (t - y2) - y
    else
        o.c = (t - y) - y2
    end

    o.sum = t
    o.n += o2.n
    o
end

#-----------------------------------------------------------------------# KahanMean
"""
    KahanMean(; T=Float64, weight=EqualWeight())

Track a univariate mean. Uses a compensation term for the update.

#Note

This should be more accurate as [`Mean`](@ref) in most cases but the guarantees
of [`KahanSum`](@ref) do not apply. [`merge!`](@ref) can have some accuracy
issues.

# Update

``μ = (1 - γ) * μ + γ * x``

# Example

    @time fit!(KahanMean(), randn(10^6))
"""
mutable struct KahanMean{W, T<:Number} <: OnlineStat{Number}
    μ::T
    c::T
    weight::W
    n::Int
end
KahanMean(T::Type) = KahanMean(T(0), T(0), EqualWeight(), 0)
KahanMean(;T::Type = Float64, weight = EqualWeight()) = KahanMean(T(0.0), T(0.0), weight, 0)
function _fit!(o::KahanMean{W, T}, x) where {W, T}
    o.n += 1

    # This acts under the assumption that the mean and all values are
    # approximately of the same order
    # o.μ = o.μ + T(o.weight(o.n)) * (convert(T, x) - o.μ) - o.c
    y = T(o.weight(o.n)) * (convert(T, x) - o.μ) - o.c
    t = o.μ + y
    o.c = (t - o.μ) - y
    o.μ = t
end
function _merge!(o::KahanMean{W1, T}, o2::KahanMean{W2, T}) where {W1, W2, T}

    o.n += o2.n
    # y = (T(o2.n / o.n) * (o2.μ - o.μ) - o.c) - o2.c
    x1 = o.μ - o2.c
    y = T(o2.n / o.n) * (o2.μ - x1) - o.c
    t = x1 + y

    if abs(x1) < abs(y)
        o.c = (t - y) - x1
    else
        o.c = (t - x1) - y
    end

    o.μ = t

    o
end
Statistics.mean(o::KahanMean) = o.μ
Base.copy(o::KahanMean) = KahanMean(o.μ, o.c, o.weight, o.n)

#-----------------------------------------------------------------------# KahanVariance
"""
    KahanVariance(; T=Float64, weight=EqualWeight())

Track the univariate variance. Uses compensation terms for a higher accuracy.

#Note

This should be more accurate as [`Variance`](@ref) in most cases but the
guarantees of [`KahanSum`](@ref) do not apply. [`merge!`](@ref) can have
accuracy issues.

# Example

    o = fit!(KahanVariance(), randn(10^6))
    mean(o)
    var(o)
    std(o)
"""
mutable struct KahanVariance{W, T<:Number} <: OnlineStat{Number}
    σ2::T
    μ::T
    cμ::T
    cσ2::T
    weight::W
    n::Int
end
KahanVariance(T::Type) =
    KahanVariance(T(0.0), T(0.0), T(0.0), T(0.0), EqualWeight(), 0)
KahanVariance(;T::Type = Float64, weight = EqualWeight()) =
    KahanVariance(T(0.0), T(0.0), T(0.0), T(0.0), weight, 0)
Base.copy(o::KahanVariance) =
    KahanVariance(o.σ2, o.μ, o.cμ, o.cσ2, o.weight, o.n)
function _fit!(o::KahanVariance{W, T}, x) where {W, T}

    o.n += 1
    γ = T(o.weight(o.n))
    μ = o.μ

    xx = convert(T, x)
    # o.μ = μ + γ * (xx - μ)
    y = γ * (xx - μ) - o.cμ
    t = μ + y
    o.cμ = (t - μ) - y
    o.μ = t

    # o.σ2 = o.σ2 + γ * ((xx - μ) * (xx - o.μ) - o.σ2)
    y = γ * ((xx - μ) * (xx - o.μ) - o.σ2) - o.cσ2
    t = o.σ2 + y
    o.cσ2 = (t - o.σ2) - y
    o.σ2 = t

    return nothing
end
function _merge!(o::KahanVariance{W1, T}, o2::KahanVariance{W2, T}) where {W1, W2, T}
    o.n += o2.n

    γ = T(o2.n / o.n)

    μ1 = o.μ - o2.cμ
    μ2 = o2.μ - o.cμ
    # Δμ = μ2 - μ1
    Δμ = μ2 - o.μ

    σ1 = o.σ2 - o2.cσ2
    σ2 = o2.σ2 - o.cσ2
    # Δσ = σ2 - σ1
    Δσ = σ2 - o.σ2

    # o.σ2 = γ * (o2.σ2 - o.σ2) + δ ^ 2 * γ * (1.0 - γ)
    # xx = (γ * Δσ) + ((Δμ ^ 2) * γ * (T(1.0) - γ))
    xx = γ * (Δσ + (Δμ ^ 2) * (T(1.0) - γ))
    y = xx - o.cσ2
    t = o.σ2 + y

    if abs(σ1) < abs(y)
        o.cσ2 = (t - y) - σ1
    else
        o.cσ2 = (t - σ1)  - y
    end

    o.σ2 = t

    # o.μ = o.μ + γ * (o2.μ - o.μ)
    # xx = γ * δ
    xx = γ * Δμ
    y = xx - o.cμ
    t = o.μ + y

    if abs(μ1) < abs(y)
        o.cμ = (t - y) - μ1
    else
        o.cμ = (t - μ1) - y
    end

    o.μ = t

    o
end
value(o::KahanVariance{W, T}) where {W, T} =
    o.n > 1 ? o.σ2 * bessel(o) : T(1.0)
Statistics.var(o::KahanVariance) = value(o)
Statistics.mean(o::KahanVariance) = o.μ


text:
using Plots, Iterators
gr()

function circle(x,y,r,c_color, myplot)
    th = 0:pi/50:2*pi;
    xunit = r * cos.(th) + x;
    yunit = r * sin.(th) + y;
    plot!(myplot, xunit, yunit,color=c_color,linewidth=3.0);
end

plt = plot()
circle(5,5, 5,:green);

function circleObs(obstacle)
    x1, y1 = obstacle.SW
    x2,y2 = obstacle.NE
    r = 0.2
    obsColor =  :black

    circle(x1,y1,r,obsColor)
    circle(x1,y2,r,obsColor)
    circle(x2,y1,r,obsColor)
    circle(x2,y2,r,obsColor)
end

obs1 = Obstacle((1,0),(2,1))
obs1 = Obstacle((1,0),(2,1))

circleObs(obs1)

############
print("\n\nend of function. \n")

module plotHelpers
    function line(startpt, finishpt, l_color)
        x1,y1 = startpt
        x2,y2 = finishpt
        plot([x1, x2], [y1,y2], color=l_color, linewidth=2)
    end

    function goalline(startpt, finishpt)
        x1,y1 = startpt
        x2,y2 = finishpt
        plot([x1, x2], [y1,y2], color="g", linewidth=2)
    end

    ### Some helper functions
    function circle(x,y,r,c_color)
        th = 0:pi/50:2*pi;
        xunit = r * cos.(th) + x;
        yunit = r * sin.(th) + y;
        h = plot(xunit, yunit,color=c_color,linewidth=3.0);
    end

    function circleObs(obstacle)
        x1, y1 = obstacle.SW
        x2,y2 = obstacle.NE
        r = 0.2
        obsColor =  :black

        circle(x1,y1,r,obsColor)
        circle(x1,y2,r,obsColor)
        circle(x2,y1,r,obsColor)
        circle(x2,y2,r,obsColor)
    end

end

###########

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

#foo = Node(0,10,(1,2))
#print(foo, "\n Type ", typeof(foo))

function inCollisionObstacle(rob,obstacle)
end

function planPath(start, goal, ballradius, maxiters)
end

function AddNode(rrt, id, iPrev)
end

function inCollision(newPoint)
end

function inCollision(Node)
end

function distance()
end

function nn() #nearest neighbors
end
# robotRad = 0.5
# rob = Node(0,0,(startx, starty)) #!
#
# rrt = rrtPath(Node[],Edge[])
#
# ## okay -- iterations
# #    iterations = 3
# #    for k = 1:iterations
# #        xy_rand = rand(0:dim)
# #
# #        newNode = Node(i, iprev, xy_rand)
# #        if inCollision(newNode)
# #            continue
# #        end
# # #        newEdge = Edge(nodePrev, newNode)

# robotRad = 0.5
# rob = Node(0,0,(startx, starty)) #!
#
# rrt = rrtPath(Node[],Edge[])
#
# ## okay -- iterations
# #    iterations = 3
# #    for k = 1:iterations
# #        xy_rand = rand(0:dim)
# #
# #        newNode = Node(i, iprev, xy_rand)
# #        if inCollision(newNode)
# #            continue
# #        end
#
# #        newEdge = Edge(nodePrev, newNode)
# #        if inCollision(newEdge):
# #            continue
# #        end
#
#
# #        push!(rrt.nodeslist,newNode)
# #        push!(rrt.edgeslist,newEdge)
# #
#
# #        addNode(rrt,
# #end
# ##
# push!(rrt.nodeslist,rob)
# push!(rrt.nodeslist,rob)
# testgoal = Node(1,0,(goalx, goaly))
# push!(rrt.nodeslist,testgoal)
#
    # line((1,1),(2,2), "red")
#
# ##
#
# #        if inCollision(newEdge):
# #            continue
# #        end
#
#
# #        push!(rrt.nodeslist,newNode)
# #        push!(rrt.edgeslist,newEdge)
# #
#
# #        addNode(rrt,
# #end
# ##
# push!(rrt.nodeslist,rob)
# push!(rrt.nodeslist,rob)
# testgoal = Node(1,0,(goalx, goaly))
# push!(rrt.nodeslist,testgoal)
#
    # line((1,1),(2,2), "red")
#
# ##
#

function plotRRT(room, paths)
    foo = rand(1)
    title("A rrt visualization $(foo)")

    dim = 6
    roomx = [0,0,dim,dim];
    roomy = [0,dim,dim,0];
    plot!(roomx, roomy, color=:black, linewidth=5)
    fig = gcf();

    # Plot the path
    for i=2:length(P)
        plot([P(1,i);P(1,i-1)],[P(2,i);P(2,i-1)],"g","LineWidth",3);
    end
end

module rrtFunctions
    struct rrtPath
        nodeslist
        edgeslist
        #nodeslist::Array{Node}
        #edgeslist::Array{Edge}
    end
end


text:
using Test

include("str_tests.jl")
include("matrix_tests.jl")
include("linkedlist_tests.jl")

max_stars_repo_name: bmoretz/Daily-Coding-Problem

max_issues_repo_name: bmoretz/Daily-Coding-Problem

max_forks_repo_name: bmoretz/Daily-Coding-Problem


text:
#=
runtests:
- Julia version:
- Author: bramb
- Date: 2019-02-20
=#
using Test

@testset "all tests" begin
    include("test_hypercollate.jl")
    include("test_collater.jl")
    include("test_xmlparser.jl")
    include("test_parser2.jl")
end