text:
using Makie: ScenePlot
import Makie.MakieLayout: registration_setup!, process_interaction

export NodeHoverHandler, EdgeHoverHandler
export NodeHoverHighlight, EdgeHoverHighlight
export NodeClickHandler, EdgeClickHandler
export NodeDragHandler, EdgeDragHandler
export NodeDrag, EdgeDrag

"""
    convert_selection(element, idx)

`mouse_selection` returns the basic plot type. In case of `Lines` check if it is
part of a `BezierSegments` and convert selection accordingly. In case of `LineSegments`
check if it is part of a `EdgePlot` and convert idx.
"""
convert_selection(element, idx) = (element, idx)
function convert_selection(element::LineSegments, idx)
    if element.parent isa EdgePlot
        return (element.parent, Int(idx / 2))
    end
    return (element, idx)
end
function convert_selection(element::Lines, idx)
    if element.parent isa BezierSegments && element.parent.parent isa EdgePlot
        bezier = element.parent
        idx = findfirst(isequal(element), bezier.plots)
        return (bezier.parent, idx)
    end
    return (element, idx)
end

"""
    abstract type GraphInteraction

This abstract type supertypes other `GraphInteractions`. Some `GraphInteractions` need
references to the scatter and linesegments plot of the `GraphPlot`. By overloding
`set_nodeplot!(T<:GraphInteraction)` and `set_edgeplot!` those references will be set
during `registration_setup!`.
"""
abstract type GraphInteraction end;
set_nodeplot!(::GraphInteraction, plot) = nothing
set_edgeplot!(::GraphInteraction, plot) = nothing

function registration_setup!(parent, inter::GraphInteraction)
    # in case of multiple graph plots in one axis this won't work
    gplots = filter(p -> p isa GraphPlot, parent.scene.plots)
    if length(gplots) !== 1
        @warn "There are multiple GraphPlots, register interaction to first!"
    end
    gplot = gplots[1]
    set_nodeplot!(inter, get_node_plot(gplot))
    set_edgeplot!(inter, get_edge_plot(gplot))
end

####
#### Hover Interaction
####

"""
    mutable struct HoverHandler{P<:ScenePlot, F} <: GraphInteraction

Object to handle hovers on `plot::P`. Calls `fun` on hover.
"""
mutable struct HoverHandler{P<:ScenePlot,F} <: GraphInteraction
    idx::Union{Nothing,Int}
    plot::Union{Nothing,P}
    fun::F
end
set_nodeplot!(h::HoverHandler{Scatter}, plot) = h.plot = plot
set_edgeplot!(h::HoverHandler{EdgePlot}, plot) = h.plot = plot

function process_interaction(handler::HoverHandler, event::MouseEvent, axis)
    if event.type === MouseEventTypes.over
        (element, idx) = convert_selection(mouse_selection(axis.scene)...)
        if element == handler.plot
            if handler.idx === nothing
                handler.idx = idx
                ret = handler.fun(true, handler.idx, event, axis)
                return ret isa Bool ? ret : false
            end
        else
            if handler.idx !== nothing
                ret = handler.fun(false, handler.idx, event, axis)
                handler.idx = nothing
                return ret isa Bool ? ret : false
            end
        end
    end
    return false
end

"""
    NodeHoverHandler(fun)

Initializes `HoverHandler` for nodes. Calls function

    fun(hoverstate, idx, event, axis)

with `hoverstate=true` on hover and `false` at the end of hover. `idx` is the node index.

# Example
```
julia> g = wheel_digraph(10)
julia> f, ax, p = graphplot(g, node_size = [20 for i in 1:nv(g)])
julia> function action(state, idx, event, axis)
           p.node_size[][idx] = state ? 40 : 20
           p.node_size[] = p.node_size[] #trigger observable
       end
julia> register_interaction!(ax, :nodehover, NodeHoverHandler(action))
```
"""
NodeHoverHandler(fun::F) where {F} = HoverHandler{Scatter,F}(nothing, nothing, fun)

"""
    NodeHoverHeighlight(p::GraphPlot, factor=2)

Magnifies the `node_size` of node under cursor by `factor`.

# Example
```
julia> g = wheel_graph(10)
julia> f, ax, p = graphplot(g, node_size = [20 for i in 1:nv(g)])
julia> register_interaction!(ax, :nodehover, NodeHoverHighlight(p))
```
"""
function NodeHoverHighlight(p::GraphPlot, factor=2)
    @assert p.node_size[] isa Vector{<:Real} "`node_size` object needs to be an Vector{<:Real} for this interaction to work!"
    action = (state, idx, _, _) -> begin
        old = p.node_size[][idx]
        p.node_size[][idx] = state ? old * factor : old / factor
        p.node_size[] = p.node_size[] #trigger observable
    end
    return NodeHoverHandler(action)
end

"""
    EdgeHoverHandler(fun)

Initializes `HoverHandler` for edges. Calls function

    fun(hoverstate, idx, event, axis)

with `hoverstate=true` on hover and `false` at the end of hover. `idx` is the edge index.

# Example
```
julia> g = wheel_digraph(10)
julia> f, ax, p = graphplot(g, edge_width = [3.0 for i in 1:ne(g)])
julia> function action(state, idx, event, axis)
           p.edge_width[][idx] = state ? 6.0 : 3.0
           p.edge_width[] = p.edge_width[] #trigger observable
       end
julia> register_interaction!(ax, :edgehover, EdgeHoverHandler(action))
```
"""
EdgeHoverHandler(fun::F) where {F} = HoverHandler{EdgePlot,F}(nothing, nothing, fun)

"""
    EdgeHoverHeighlight(p::GraphPlot, factor=2)

Magnifies the `edge_width` of edge under cursor by `factor`.
If `arrow_size isa Vector{<:Real}` it also magnefies the arrow scatter.

# Example
```
julia> g = wheel_digraph(10)
julia> f, ax, p = graphplot(g, edge_width = [3 for i in 1:ne(g)],
                               arrow_size=[10 for i in 1:ne(g)])
julia> register_interaction!(ax, :nodehover, EdgeHoverHighlight(p))
```
"""
function EdgeHoverHighlight(p::GraphPlot, factor=2)
    @assert p.edge_width[] isa Vector{<:Real} "`edge_width` object needs to be an Vector{<:Real} for this interaction to work!"
    scale_arrows = p.arrow_size[] isa Vector{<:Real}

    if length(p.edge_width[]) == ne(p[:graph][])
        action = (state, idx, _, _) -> begin
            old = p.edge_width[][idx]
            p.edge_width[][idx] = state ? old * factor : old / factor
            p.edge_width[] = p.edge_width[] #trigger observable
            if scale_arrows
                old = p.arrow_size[][idx]
                p.arrow_size[][idx] = state ? old * factor : old / factor
                p.arrow_size[] = p.arrow_size[] #trigger observable
            end
        end
    elseif length(p.edge_width[]) == 2 * ne(p[:graph][])
        action = (state, idx, _, _) -> begin
            oldA = p.edge_width[][2 * idx - 1]
            oldB = p.edge_width[][2 * idx]
            p.edge_width[][2 * idx - 1] = state ? oldA * factor : oldA / factor
            p.edge_width[][2 * idx] = state ? oldB * factor : oldB / factor
            p.edge_width[] = p.edge_width[] #trigger observable
            if scale_arrows
                old = p.arrow_size[][idx]
                p.arrow_size[][idx] = state ? old * factor : old / factor
                p.arrow_size[] = p.arrow_size[] #trigger observable
            end
        end
    else
        error("Can not make sense of `length(p.edge_width)`")
    end

    return EdgeHoverHandler(action)
end

####
#### Drag Interaction
####

"""
    mutable struct DragHandler{P<:ScenePlot, F} <: GraphInteraction

Object to handle left mous drags on `plot::P`.
"""
mutable struct DragHandler{P<:ScenePlot,F} <: GraphInteraction
    dragstate::Bool
    idx::Int
    plot::Union{Nothing,P}
    fun::F
end
set_nodeplot!(h::DragHandler{Scatter}, plot) = h.plot = plot
set_edgeplot!(h::DragHandler{EdgePlot}, plot) = h.plot = plot

function process_interaction(handler::DragHandler, event::MouseEvent, axis)
    if handler.dragstate == false # not in drag state
        if event.type === MouseEventTypes.leftdown
            # on leftdown save idx if happens to be on right element/plot
            (element, idx) = convert_selection(mouse_selection(axis.scene)...)
            handler.idx = element === handler.plot ? idx : 0
        elseif event.type === MouseEventTypes.leftdragstart && handler.idx != 0
            # if idx!=0 the last leftdown was on right element!
            handler.dragstate = true
        end
    elseif handler.dragstate == true # drag state
        if event.type === MouseEventTypes.leftdrag
            ret = handler.fun(true, handler.idx, event, axis)
            return ret isa Bool ? ret : false
        elseif event.type === MouseEventTypes.leftdragstop
            ret = handler.fun(false, handler.idx, event, axis)
            handler.idx = 0
            handler.dragstate = false
            return ret isa Bool ? ret : false
        end
    end
    return false
end

"""
    NodeDragHandler(fun)

Initializes `DragHandler` for Nodes. Calls function

    fun(dragstate, idx, event, axis)

where `dragstate=true` during the drag and `false` at the end of the drag,
the last time `fun` is triggered. `idx` is the node index.

# Example
```
julia> g = wheel_digraph(10)
julia> f, ax, p = graphplot(g, node_size=20)
julia> deregister_interaction!(ax, :rectanglezoom)
julia> function action(state, idx, event, axis)
           p[:node_pos][][idx] = event.data
           p[:node_pos][] = p[:node_pos][]
       end
julia> register_interaction!(ax, :nodedrag, NodeDragHandler(action))
```
"""
NodeDragHandler(fun::F) where {F} = DragHandler{Scatter,F}(false, 0, nothing, fun)

"""
    NodeDrag(p::GraphPlot)

Allows drag and drop of Nodes. Please deregister the `:rectanglezoom` interaction.

# Example
```
julia> g = wheel_graph(10)
julia> f, ax, p = graphplot(g, node_size = [10 for i in 1:nv(g)])
julia> deregister_interaction!(ax, :rectanglezoom)
julia> register_interaction!(ax, :nodehover, NodeHoverHighlight(p))
julia> register_interaction!(ax, :nodedrag, NodeDrag(p))
```
"""
function NodeDrag(p)
    action = (state, idx, event, _) -> begin
        p[:node_pos][][idx] = event.data
        p[:node_pos][] = p[:node_pos][]
    end
    return NodeDragHandler(action)
end

"""
    EdgeDragHandler(fun)

Initializes `DragHandler` for Edges. Calls function

    fun(dragstate, idx, event, axis)

where `dragstate=true` during the drag and `false` at the end of the drag,
the last time `fun` is triggered. `idx` is the edge index.

See [`EdgeDrag`](@ref) for a concrete implementation.
```
"""
EdgeDragHandler(fun::F) where {F} = DragHandler{EdgePlot,F}(false, 0, nothing, fun)

"""
    EdgeDrag(p::GraphPlot)

Allows drag and drop of Edges. Please deregister the `:rectanglezoom` interaction.

# Example
```
julia> g = wheel_graph(10)
julia> f, ax, p = graphplot(g, edge_width = [3 for i in 1:ne(g)])
julia> deregister_interaction!(ax, :rectanglezoom)
julia> register_interaction!(ax, :edgehover, EdgeHoverHighlight(p))
julia> register_interaction!(ax, :edgedrag, EdgeDrag(p))
```
"""
function EdgeDrag(p)
    action = EdgeDragAction(p)
    return EdgeDragHandler(action)
end

mutable struct EdgeDragAction{PT<:GraphPlot}
    p::PT
    init::Union{Nothing,Point2f} # save click position
    src::Union{Nothing,Point2f}  # save src vertex position
    dst::Union{Nothing,Point2f}  # save dst vertex position
    EdgeDragAction(p::T) where {T} = new{T}(p, nothing, nothing, nothing)
end

function (action::EdgeDragAction)(state, idx, event, _)
    edge = collect(edges(action.p[:graph][]))[idx]
    if state == true
        if action.src === action.dst === action.init === nothing
            action.init = event.data
            action.src = action.p[:node_pos][][edge.src]
            action.dst = action.p[:node_pos][][edge.dst]
        end
        offset = event.data - action.init
        action.p[:node_pos][][edge.src] = action.src + offset
        action.p[:node_pos][][edge.dst] = action.dst + offset
        action.p[:node_pos][] = action.p[:node_pos][] # trigger change
    elseif state == false
        action.src = action.dst = action.init = nothing
    end
end

####
#### Click Interaction
####

"""
    mutable struct ClickHandler{P<:ScenePlot, F} <: GraphInteraction

Object to handle left mouse clicks on `plot::P`.
"""
mutable struct ClickHandler{P<:ScenePlot,F} <: GraphInteraction
    plot::Union{Nothing,P}
    fun::F
end
set_nodeplot!(h::ClickHandler{Scatter}, plot) = h.plot = plot
set_edgeplot!(h::ClickHandler{EdgePlot}, plot) = h.plot = plot

function process_interaction(handler::ClickHandler, event::MouseEvent, axis)
    if event.type === MouseEventTypes.leftclick
        (element, idx) = convert_selection(mouse_selection(axis.scene)...)
        if element == handler.plot
            ret = handler.fun(idx, event, axis)
            return ret isa Bool ? ret : false
        end
    end
    return false
end

"""
    NodeClickHandler(fun)

Initializes `ClickHandler` for nodes. Calls function

    fun(idx, event, axis)

on left-click events where `idx` is the node index.

# Example
```
julia> using Makie.Colors
julia> g = wheel_digraph(10)
julia> f, ax, p = graphplot(g, node_size=30, node_color=[colorant"red" for i in 1:nv(g)])
julia> function action(idx, event, axis)
           p.node_color[][idx] = rand(RGB)
           p.node_color[] = p.node_color[]
       end
julia> register_interaction!(ax, :nodeclick, NodeClickHandler(action))
```
"""
NodeClickHandler(fun::F) where {F} = ClickHandler{Scatter,F}(nothing, fun)

"""
    EdgeClickHandler(fun)

Initializes `ClickHandler` for edges. Calls function

    fun(idx, event, axis)

on left-click events where `idx` is the edge index.

# Example
```
julia> using Makie.Colors
julia> g = wheel_digraph(10)
julia> f, ax, p = graphplot(g, edge_width=4, edge_color=[colorant"black" for i in 1:ne(g)])
julia> function action(idx, event, axis)
           p.edge_color[][idx] = rand(RGB)
           p.edge_color[] = p.edge_color[]
       end
julia> register_interaction!(ax, :edgeclick, EdgeClickHandler(action))
```
"""
EdgeClickHandler(fun::F) where {F} = ClickHandler{EdgePlot,F}(nothing, fun)


text:
const c = Center()
const f = Face()

using Oceananigans.AbstractOperations: GridMetricOperation
import Oceananigans.Grids: solid_node, solid_interface

"""
    `solid_node` returns true only if a location is completely immersed
    `solid_interface` returns true if a location is partially immersed

    as an example (in a 1D immersed grid) :

     Immersed      Fluid
    ----------- ...........
   |     ∘     |     ∘
   f     c     f     c
  i-1   i-1    i     i

     `solid_interface(f, c, c, i, 1, 1, grid) = true`
     `solid_node(f, c, c, i, 1, 1, grid)      = false`

     `solid_node(Center(), Center(), Center(), args...) == solid_interface(Center(), Center(), Center(), args...)` as
     `Center(), Center(), Center()` can be only either fully immersed or not at all

     `solid_interface` is used in `GridMetricOperation` to assess the grid metric value.
     `metric = ifelse(solid_interface, 0.0, metric)` eliminating values inside the immersed domain and values lying on the
      immersed boundary from the reduction

      `is_immersed_boundary` returns true only if the interface has a solid and a fluid side (the actual immersed boundary)
      which is true only when `solid_node = false` and `solid_interface = true` (as the case of the face at `i` above)
"""

@inline solid_node(i, j, k, ibg::IBG) = is_immersed(i, j, k, ibg.grid, ibg.immersed_boundary) | solid_node(i, j, k, ibg.grid)

# Defining all the metrics for Immersed Boundaries

for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ)
    for dir in (:x, :y, :z), operator in (:Δ, :A)

        metric = Symbol(operator, dir, LX, LY, LZ)
        @eval begin
            import Oceananigans.Operators: $metric
            @inline $metric(i, j, k, ibg::IBG) = $metric(i, j, k, ibg.grid)
        end
    end

    volume = Symbol(:V, LX, LY, LZ)
    @eval begin
        import Oceananigans.Operators: $volume
        @inline $volume(i, j, k, ibg::IBG) = $volume(i, j, k, ibg.grid)
    end
end

@inline Δzᵃᵃᶜ(i, j, k, ibg::IBG) = Δzᵃᵃᶜ(i, j, k, ibg.grid)
@inline Δzᵃᵃᶠ(i, j, k, ibg::IBG) = Δzᵃᵃᶠ(i, j, k, ibg.grid)

ext: jl


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   General tests.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "stdout" begin
    old_stdout = stdout

    in, out = redirect_stdout()
    pretty_number(1.19861987e6)
    result = String(readavailable(in))
    redirect_stdout(old_stdout)

    expected = "1.19862 × 10⁶"

    @test result == expected
end


text:
function dns_service_browse(service_type)

    result = Channel{DNSService}(0)

    # Browse for service instances and display output in zone file format.
    p = open(`dns-sd -Z $service_type`)

    atexit(()->kill(p, Base.SIGKILL))

    @async try
        while isopen(result) &&
              (process_running(p) || bytesavailable(p) > 0)

            l = readline(p)
            if isempty(l) || startswith(l, ";")
                continue
            end
            l = split(l)
            if l[2] == "SRV"
                # e.g.: FooBar._tinyrpc._tcp  SRV 0 0 2020 foobar.local. ; ..."
                service, SRV, priority, weight, port, target = l
                name = split(service, ".")[1]
                push!(result, name => (target, parse(Int16, port)))
            end
        end
    catch err
        if !isopen(result) &&
           err isa InvalideStateException &&
           err.state == closed
            # Ignore push to closed Channel error.
        else
            exception=(err, catch_backtrace())
            @error "Error reading `dns-sd -Z` output." exception
        end
    finally
        close(result)
        kill(p)
    end

    return result
end

function register_dns_service(name, service_type, port)
    global dns_process
    dns_process[(name, service_type)] =
    open(`dns-sd -R $name $service_type local $port`)
    nothing
end


text:
module BRS
using Parameters
using RecipesBase
using Statistics
using DSP
using Dierckx
using FFTW
using Polynomials

include("sequenceMethod.jl")
include("xBRS.jl")
include("tfBRS.jl")
include("prsaBRS.jl")
include("arBRS.jl")
include("../Utilities/util.jl")

"""
        getBRS(RR::Vector{<:Real}, SBP::Vector{<:Real})

Returns all BRS measures with default settings for individual methods. For more control use each method individually.

# Methods

* SME: Sequence Method
* RMSSDR: RMSSD ratio
* xBRS: Cross-correlation baroreflex sensitivity
* TFBRS: Transfer Function based BRS (fft)
* PRSABRS: Phase Rectified Signal Averaging
* αLF: LF component of the AR based spectral decomposition
* αHF: HF component of the AR based spectral decomposition
"""
function getBRS(RR::Vector{<:Real}, SBP::Vector{<:Real})
        valueSME = sme(RR, SBP).value
        valueRMSSDR = rmssdr(RR, SBP)
        valueXBRS = xbrs(RR, SBP).value
        valueTFBRS = tfbrs(RR, SBP).value
        valuePRSABRS = prsabrs(RR, SBP).value
        resAR = arbrs(RR, SBP)
        αLF = resAR.αLF
        αHF = resAR.αHF

        return (sme = valueSME, rmssdr = valueRMSSDR, xbrs = valueXBRS, tfbrs = valueTFBRS, prsabrs = valuePRSABRS, αLF = αLF, αHF = αHF)
end

export sme, SME,
rmssdr,
xbrs, xBRS,
tfbrs, tfBRS,
prsabrs, prsaBRS,
arbrs, arBRS, arDecomposition, getSpectralComponent,
getBRS

end # end module


text:
struct TURBO{UNROLL,OPS,ARF,AM,LPSYM,LBV,FLBV} <: Function end

# This should call the same `_turbo_!(Val{UNROLL}(), Val{OPS}(), Val{ARF}(), Val{AM}(), Val{LPSYM}(), _vargs)` as normal so that this
# hopefully shouldn't add much to compile time.

function (::TURBO{UNROLL,OPS,ARF,AM,LPSYM,LBV,FLBV})(p::Ptr{UInt}) where {UNROLL,OPS,ARF,AM,LPSYM,K,LBV,FLBV<:Tuple{Vararg{Any,K}}}
  (_, _vargs) = ThreadingUtilities.load(p, FLBV, 2*sizeof(UInt))
  # Main.VARGS[Threads.threadid()] = first(_vargs)
  # Threads.threadid() == 2 && Core.println(typeof(_vargs))
  ret = _turbo_!(Val{UNROLL}(), Val{OPS}(), Val{ARF}(), Val{AM}(), Val{LPSYM}(), Val{LBV}(), _vargs...)
  ThreadingUtilities.store!(p, ret, Int(register_size()))
  ThreadingUtilities._atomic_store!(p, ThreadingUtilities.SPIN)
  nothing
end
@generated function Base.pointer(::TURBO{UNROLL,OPS,ARF,AM,LPSYM,LBV,FLBV}) where {UNROLL,OPS,ARF,AM,LPSYM,K,LBV,FLBV<:Tuple{Vararg{Any,K}}}
    f = TURBO{UNROLL,OPS,ARF,AM,LPSYM,LBV,FLBV}()
    precompile(f, (Ptr{UInt},))
    quote
        $(Expr(:meta,:inline))
        @cfunction($f, Cvoid, (Ptr{UInt},))
    end
end

@inline function setup_turbo_threads!(p::Ptr{UInt}, fptr::Ptr{Cvoid}, args::LBV) where {K,LBV<:Tuple{Vararg{Any,K}}}
    offset = ThreadingUtilities.store!(p, fptr, sizeof(UInt))
    offset = ThreadingUtilities.store!(p, args, offset)
    nothing
end

struct StaticType{T} end
@inline gettype(::StaticType{T}) where {T} = T

@inline function avx_launch(
    ::Val{UNROLL}, ::Val{OPS}, ::Val{ARF}, ::Val{AM}, ::Val{LPSYM}, ::StaticType{LBV}, fargs::FARGS, tid
) where {UNROLL,OPS,ARF,AM,LPSYM,K,LBV<:Tuple{Vararg{Any,K}},FARGS}
    ThreadingUtilities.launch(setup_turbo_threads!, tid, pointer(TURBO{UNROLL,OPS,ARF,AM,LPSYM,LBV,FARGS}()), fargs)
end

# function approx_cbrt(x)
#     s = significand(x)
#     e = exponent(x)

#     # 40 + 0.00020833333333333335*(x-64000)  -2.1701388888888896e-9*(x-64000)^2*0.5 + 5.6514033564814844e-14 * (x-64000)^3/6
# end
lv_max_num_threads() = ifelse(gt(num_threads(), num_cores()), num_cores(), num_threads())

@generated function calc_factors(::StaticInt{nc}) where {nc}
    t = Expr(:tuple)
    for i ∈ nc:-1:1
        d, r = divrem(nc, i)
        iszero(r) && push!(t.args, (i % UInt, d % UInt))
    end
    t
end
@inline cld_fast(x,y) = Base.udiv_int(vsub_nw(vadd_nw(x, y), one(y)), y)
@inline function choose_num_blocks(MoW::UInt, ::StaticInt{U}, ::StaticInt{T}) where {U,T}
    factors = calc_factors(StaticInt{T}())
    for i ∈ 1:length(factors)-1
        # miter decreases in each iteration of factors
        miter, niter = factors[i]

        r = (MoW % miter)
        # if ((miter * W * U * 2) ≤ M - (W+W)) & ((r == 0) | (miter == (r+1)))
        mlarge = (miter * (U * 2)) ≤ MoW - 2
        # we want `mlarge` enough, or there to be no remainder (`r == 0`)
        # or the remainder to be small relative to `M`, e.g. `M ÷ 128`.
        if mlarge && ((r == 0) || ((miter - r) ≤ (cld_fast(MoW, miter) >>> 3)))
            return miter, niter
        end
    end
    last(factors)
end
@inline function choose_num_blocks(::StaticInt{T}) where {T}
    factors = calc_factors(StaticInt{T}())
    @inbounds factors[(length(factors) + 1)>>>1]
end

# struct ChooseNumBlocks{U,C} <: Function end
# function (cnb::ChooseNumBlocks{U,C})(M::UInt) where {U,C}
#     choose_num_blocks(M, StaticInt{U}(), StaticInt{C}())
# end

# @generated function choose_num_block_table(::StaticInt{U}, ::StaticInt{NC}) where {U,NC}
#     t = Expr(:tuple)
#     for n ∈ 1:NC
#         cnb = :(ChooseNumBlocks{$U,$n}())
#         push!(t.args, :(@cfunction($cnb, Tuple{UInt,UInt}, (UInt,))))
#     end
#     t
# end
@generated function choose_num_block_table(::StaticInt{NC}) where {NC}
    t = Expr(:tuple)
    for n ∈ 1:NC
        push!(t.args, :(choose_num_blocks(StaticInt{$n}())))
    end
    t
end

@generated function _choose_num_blocks(M::UInt, ::StaticInt{U}, nt, ::StaticInt{NTMAX}) where {U,NTMAX}
    # valid range for nt: 2 ≤ nt ≤ NTMAX
    # if NTMAX > 8
    #     return quote
    #         $(Expr(:meta,:inline))
    #         choose_num_blocks_table(M, StaticInt{$U}(), nt, StaticInt{$NTMAX}())
    #     end
    # else
    if NTMAX == 2 # `nt` must be `2`
        return quote
            $(Expr(:meta,:inline))
            choose_num_blocks(M, StaticInt{$U}(), StaticInt{$NTMAX}())
        end
    end
    q = Expr(:block)#, Expr(:meta,:inline))
    ifq = Expr(
        :if,
        :(nt == $NTMAX),
        :(choose_num_blocks(M, StaticInt{$U}(), StaticInt{$NTMAX}()))
    )
    add_bisecting_if_branches!(ifq, 2, NTMAX-1, U, false)
    push!(q.args, ifq)
    q
end
function add_bisecting_if_branches!(q, lb, ub, U, isfirst::Bool)
    if lb == ub
        push!(q.args, :(choose_num_blocks(M, StaticInt{$U}(), StaticInt{$lb}())))
        return
    end
    midpoint = (lb + ub) >> 1
    alt = Expr(isfirst ? :if : :elseif, :(nt > $midpoint))
    add_bisecting_if_branches!(alt, midpoint+1, ub, U, true)
    add_bisecting_if_branches!(alt, lb, midpoint, U, false)
    push!(q.args, alt)
    return
end

# @inline function choose_num_blocks_table(M, ::StaticInt{U}, nt, ::StaticInt{NTMAX}) where {U,NTMAX}
#     if nt == NTMAX
#         choose_num_blocks(M % UInt, StaticInt{U}(), StaticInt{NTMAX}())
#     else
#         @inbounds fptr = choose_num_block_table(StaticInt{U}(), StaticInt{NTMAX}())[nt]
#         VectorizationBase.assume(fptr ≠ C_NULL)
#         ccall(fptr, Tuple{UInt,UInt}, (UInt,), M%UInt)
#     end
# end

# if a threaded loop is vectorized, call
@inline function choose_num_blocks(M, ::StaticInt{U}, nt) where {U}
    _choose_num_blocks(M % UInt, StaticInt{U}(), nt, lv_max_num_threads())
end
# otherwise, call
@inline choose_num_blocks(nt, ::StaticInt{NC} = lv_max_num_threads()) where {NC} = @inbounds choose_num_block_table(StaticInt{NC}())[nt]

# The goal is to minimimize the maximum costs...
# But maybe 'relatively even sizes' heuristics are more robust than fancy modeling?
# At least early on, before lots of test cases with different sorts of loops have informed the modeling.
#
# goal is to produce `nblocks` roughly even block sizes (bM, bN), such that `bM % fM == bN % fN == 0`.
# function roughly_even_blocks(M, N, fM, fN, nblocks)
#     M_N_ratio = M / N
#     block_per_m = sqrt(nblocks * M_N_ratio) # obv not even
#     blocks_per_n = block_per_m / M_N_ratio
#     mi = cld(M, fM)
#     ni = cld(N, fN)
#     block_per_m, blocks_per_n
# end
if Sys.ARCH === :x86_64
  @inline function choose_num_threads(C::T, NT::UInt, x::Base.BitInteger) where {T<:Union{Float32,Float64}}
    _choose_num_threads(Base.mul_float_fast(T(C), T(0.0225)), NT, x)
  end
else
  @inline function choose_num_threads(C::T, NT::UInt, x::Base.BitInteger) where {T<:Union{Float32,Float64}}
    _choose_num_threads(Base.mul_float_fast(C, T(0.0225) * T(0.25)), NT, x)
  end
end
@inline function _choose_num_threads(C::T, NT::UInt, x::Base.BitInteger) where {T<:Union{Float32,Float64}}
  min(Base.fptoui(UInt, Base.ceil_llvm(Base.mul_float_fast(C, Base.sqrt_llvm_fast(Base.uitofp(T, x))))), NT)
end
function push_loop_length_expr!(q::Expr, ls::LoopSet)
    l = 1
    ndynamic = 0
    mulexpr = length(ls.loops) == 1 ? q : Expr(:call, lv(:vmul_nw))
    for loop ∈ ls.loops
        if isstaticloop(loop)
            l *= length(loop)
        else
            ndynamic += 1
            if ndynamic < 3
                push!(mulexpr.args, loop.lensym)
            else
                mulexpr = Expr(:call, lv(:vmul_nw), mulexpr, loop.lensym)
            end
        end
    end
    if length(ls.loops) == 1
        ndynamic == 0 && push!(q.args, l)
    elseif l == 1
        push!(q.args, mulexpr)
    elseif ndynamic == 0
        push!(q.args, l)
    elseif ndynamic == 1
        push!(mulexpr.args, l)
        push!(q.args, mulexpr)
    else
        push!(q.args, Expr(:call, :vmul_nw, mulexpr, l))
    end
    nothing
end
@inline function divrem_fast(numerator, denominator)
    d = Base.udiv_int(numerator, denominator)
    r = numerator - denominator*d
    d, r
end

function outer_reduct_combine_expressions(ls::LoopSet, retv)
  gf = GlobalRef(Core, :getfield)
  q = Expr(:block, :(var"#load#thread#ret#" = $gf(ThreadingUtilities.load(var"#thread#ptr#", typeof($retv), $(reg_size(ls))),2,false)))
  # push!(q.args, :(@show var"#load#thread#ret#"))
  for (i,or) ∈ enumerate(ls.outer_reductions)
    op = ls.operations[or]
    var = name(op)
    mvar = mangledvar(op)
    instr = instruction(op)
    out = Symbol(mvar, "##onevec##")
    instrcall::Expr = if instr.instr ≢ :ifelse
      Expr(:call, lv(reduce_to_onevecunroll(instr)))
    else
      reductexpr::Expr = let ls = ls
        ifelse_reduction(:IfElseReduced, op) do opv
          @assert length(ls.outer_reductions) > 1
          j = findfirst(==(identifier(opv)), ls.outer_reductions)
          otherarg = Expr(:call, lv(:vecmemaybe), Expr(:call, GlobalRef(Core,:getfield), Symbol("#load#thread#ret#"), j, false))
          Expr(:call, lv(:vecmemaybe), Symbol(mangledvar(opv), "##onevec##")), (otherarg,)
        end
      end
      Expr(:call, reductexpr)
    end
    push!(instrcall.args, Expr(:call, lv(:vecmemaybe), out))
    if length(ls.outer_reductions) > 1
      push!(instrcall.args, Expr(:call, lv(:vecmemaybe), Expr(:call, gf, Symbol("#load#thread#ret#"), i, false)))
    else
      push!(instrcall.args, Expr(:call, lv(:vecmemaybe), Symbol("#load#thread#ret#")))
    end
    push!(q.args, Expr(:(=), out, Expr(:call, :data, instrcall)))
    # push!(q.args, Expr(:(=), out, :(@show $data($instrcall))))
  end
  q
end

function thread_loop_summary!(ls::LoopSet, ua::UnrollArgs, threadedloop::Loop, issecondthreadloop::Bool)
  W = ls.vector_width
  @unpack u₁loop, u₂loop, vloop, u₁, u₂max = ua
  u₂ = u₂max
  threadloopnumtag = Int(issecondthreadloop)
  lensym = Symbol("#len#thread#$threadloopnumtag#")
  define_len = if isstaticloop(threadedloop)
    :($lensym = $(length(threadedloop)) % UInt)
  else
    :($lensym = $((threadedloop.lensym)) % UInt)
  end
  unroll_factor = Core.ifelse(threadedloop === vloop, W, 1)
  num_unroll_sym = Symbol("#num#unrolls#thread#$threadloopnumtag#")
  define_num_unrolls = if unroll_factor == 1
    :($num_unroll_sym = $lensym)
  else
    # :($num_unroll_sym = Base.udiv_int($lensym, $(UInt(unroll_factor))))
    :($num_unroll_sym = Base.udiv_int(vadd_nw($lensym, $(UInt(unroll_factor-1))), $(UInt(unroll_factor))))
  end
  iterstart_sym = Symbol("#iter#start#$threadloopnumtag#")
  iterstop_sym = Symbol("#iter#stop#$threadloopnumtag#")
  blksz_sym = Symbol("#nblock#size#thread#$threadloopnumtag#")
  loopstart = if isknown(first(threadedloop))
    :($iterstart_sym::Int = $(gethint(first(threadedloop))))
  else
    :($iterstart_sym::Int = $(getsym(first(threadedloop))))
  end
  if isknown(step(threadedloop))
    mf = gethint(step(threadedloop))
    if isone(mf)
      iterstop = :($iterstop_sym::Int = vadd_nsw($iterstart_sym, $blksz_sym))
      looprange = :(CloseOpen($iterstart_sym))
      lastrange = :(CloseOpen($iterstart_sym))
      push_loopbound_ends!(looprange, lastrange, unroll_factor, threadedloop, iterstop_sym, true)
    else
      iterstop = :($iterstop_sym::Int = vadd_nsw($iterstart_sym, vmul_nw($blksz_sym, $mf)))
      looprange = :($iterstart_sym:StaticInt{$mf}())
      lastrange = :($iterstart_sym:StaticInt{$mf}())
      push_loopbound_ends!(looprange, lastrange, unroll_factor, threadedloop, :(vsub_nsw($iterstop_sym,one($iterstop_sym))), false)
    end
  else
    stepthread_sym = Symbol("#step#thread#$threadloopnumtag#")
    pushpreamble!(ls, :($stepthread_sym = $(getsym(step(threadedloop)))))
    iterstop = :($iterstop_sym = vadd_nsw($iterstart_sym, vmul_nw($blksz_sym, $stepthread_sym)))
    looprange = :($iterstart_sym:$stepthread_sym)
    lastrange = :($iterstart_sym:$stepthread_sym)
    push_loopbound_ends!(looprange, lastrange, unroll_factor, threadedloop, :(vsub_nsw($iterstop_sym,one($iterstop_sym))), false)
  end
  define_len, define_num_unrolls, loopstart, iterstop, looprange, lastrange
end
function push_last_bound!(looprange::Expr, lastrange::Expr, lastexpr, iterstop, unroll_factor::Int)
  push!(lastrange.args, lastexpr)
  unroll_factor ≠ 1 && push!(looprange.args, :(min($lastexpr, $iterstop)))
  nothing
end
function push_loopbound_ends!(
  looprange::Expr, lastrange::Expr, unroll_factor::Int,
  threadedloop::Loop, iterstop, offsetlast::Bool
  )
  if unroll_factor == 1
    push!(looprange.args, iterstop)
  end
  if isknown(last(threadedloop))
    push_last_bound!(looprange, lastrange, gethint(last(threadedloop)) + offsetlast, iterstop, unroll_factor)
  else
    lastsym = getsym(last(threadedloop))
    if offsetlast
      push_last_bound!(looprange, lastrange, :(vadd_nsw($lastsym, one($lastsym))), iterstop, unroll_factor)
    else
      push_last_bound!(looprange, lastrange, lastsym, iterstop, unroll_factor)
    end
  end
  nothing
end
function define_block_size(threadedloop, vloop, tn, W)
  baseblocksizeuint = Symbol("#base#block#size#thread#uint#$tn#")
  baseblocksizeint = Symbol("#base#block#size#thread#$tn#")
  nrem = Symbol("#nrem#thread#$tn#")
  remstep = Symbol("#block#rem#step#$tn#")
  num_unroll = Symbol("#num#unrolls#thread#$tn#")
  thread_factor = Symbol("#thread#factor#$tn#")
  if threadedloop === vloop
    quote
      $baseblocksizeuint, $nrem = divrem_fast($num_unroll, $thread_factor)
      $baseblocksizeint = ($baseblocksizeuint << $(VectorizationBase.intlog2(W))) % Int
      $remstep = $(Int(W))
    end
  else
    quote
      $baseblocksizeuint, $nrem = divrem_fast($num_unroll, $thread_factor)
      $baseblocksizeint = $baseblocksizeuint % Int
      $remstep = 1
    end
  end
end
function scale_cost(c, looplen)
  c = 0.05 * c / looplen
  if Sys.ARCH !== :x86_64
    c *= 0.25
  end
  c
end
function thread_one_loops_expr(
  ls::LoopSet, ua::UnrollArgs, valid_thread_loop::Vector{Bool}, ntmax::UInt, c::Float64,
  UNROLL::Tuple{Bool,Int8,Int8,Int8,Bool,Int,Int,Int,Int,Int,Int,Int,UInt}, OPS::Expr, ARF::Expr, AM::Expr, LPSYM::Expr
)
  looplen = looplengthprod(ls)
  c = scale_cost(c, looplen)
  if all(isstaticloop, ls.loops)
    _num_threads = _choose_num_threads(c, ntmax, Int64(looplen))::UInt
    _num_threads > 1 || return avx_body(ls, UNROLL)
    choose_nthread = Expr(:(=), Symbol("#nthreads#"), _num_threads)
  else
    choose_nthread = :(_choose_num_threads($(Float32(c)), $ntmax))
    push_loop_length_expr!(choose_nthread, ls)
    choose_nthread = Expr(:(=), Symbol("#nthreads#"), choose_nthread)
  end
  threadedid = findfirst(valid_thread_loop)::Int
  threadedloop = getloop(ls, threadedid)
  define_len, define_num_unrolls, loopstart, iterstop, looprange, lastrange = thread_loop_summary!(ls, ua, threadedloop, false)
  loopboundexpr = Expr(:tuple) # for launched threads
  lastboundexpr = Expr(:tuple) # remainder, started on main thread
  for (i,loop) ∈ enumerate(ls.loops)
    if loop === threadedloop
      push!(loopboundexpr.args, looprange)
      push!(lastboundexpr.args, lastrange)
    else
      loop_boundary!(loopboundexpr, ls, loop, true)
      loop_boundary!(lastboundexpr, ls, loop, true)
    end
  end
  avxcall_args = Expr(:tuple, lastboundexpr, Symbol("#vargs#"))
  _turbo_call_ = :(_turbo_!(Val{$UNROLL}(), $OPS, $ARF, $AM, $LPSYM, Val(typeof(var"#avx#call#args#")), flatten_to_tuple(var"#avx#call#args#")...))
  update_return_values = if length(ls.outer_reductions) > 0
    retv = loopset_return_value(ls, Val(false))
    _turbo_call_ = Expr(:(=), retv, _turbo_call_)
    outer_reduct_combine_expressions(ls, retv)
  else
    nothing
  end
  retexpr = length(ls.outer_reductions) > 0 ? :(return $retv) : :(return nothing)
  iterdef = define_block_size(threadedloop, ua.vloop, 0, ls.vector_width)
  q = quote
    $choose_nthread # UInt
    $loopstart
    var"##do#thread##" = var"#nthreads#" > one(var"#nthreads#")
    if var"##do#thread##"
      $define_len
      $define_num_unrolls
      var"#nthreads#" = Base.min(var"#nthreads#", var"#num#unrolls#thread#0#")
      var"#nrequest#" = vsub_nw((var"#nthreads#" % UInt32), 0x00000001)
      var"##do#thread##" = var"#nrequest#" ≠ 0x00000000
      if var"##do#thread##"
        var"#threads#tuple#", var"#torelease#tuple#" = PolyesterWeave.request_threads(var"#nrequest#")
        # var"#threads#tuple#", var"#torelease#tuple#" = PolyesterWeave.request_threads(Threads.threadid()%UInt32, var"#nrequest#")

        var"#thread#factor#0#" = var"#nthreads#"
        $iterdef
        var"#thread#id#" = 0x00000000
        var"##do#thread##" = false
        for var"#threads#" in var"#threads#tuple#"

          var"#thread#launch#count#" = 0x00000000
          var"#thread#mask#" = PolyesterWeave.mask(var"#threads#")
          var"#nrequest#" = length(var"#threads#")
          var"#threads#remain#" = var"#thread#launch#count#" ≠ var"#nrequest#"
          var"##do#thread##" |= var"#threads#remain#"
          while var"#threads#remain#"
            VectorizationBase.assume(var"#thread#mask#" ≠ zero(var"#thread#mask#"))
            var"#trailzing#zeros#" = Base.trailing_zeros(var"#thread#mask#") % UInt32
            var"#nblock#size#thread#0#" = Core.ifelse(
              var"#thread#launch#count#" < (var"#nrem#thread#0#" % UInt32),
              vadd_nw(var"#base#block#size#thread#0#", var"#block#rem#step#0#"),
              var"#base#block#size#thread#0#"
            )
            var"#trailzing#zeros#" = vadd_nw(var"#trailzing#zeros#", 0x00000001)
            $iterstop
            var"#thread#id#" = vadd_nw(var"#thread#id#", var"#trailzing#zeros#")

            var"##lbvargs#to_launch##" = ($loopboundexpr, var"#vargs#")
            avx_launch(Val{$UNROLL}(), $OPS, $ARF, $AM, $LPSYM, StaticType{typeof(var"##lbvargs#to_launch##")}(), flatten_to_tuple(var"##lbvargs#to_launch##"), var"#thread#id#")

            var"#thread#mask#" >>>= var"#trailzing#zeros#"

            var"#iter#start#0#" = var"#iter#stop#0#"
            var"#thread#launch#count#" = vadd_nw(var"#thread#launch#count#", 0x00000001)
            var"#threads#remain#" = var"#thread#launch#count#" ≠ var"#nrequest#"
          end
          var"#nrem#thread#0#" -= var"#nrequest#"
        end
      else# eliminate undef var errors that the compiler should be able to figure out are unreachable, but doesn't
        var"#torelease#tuple#" = (zero(PolyesterWeave.worker_type()),)
        var"#threads#tuple#" = (PolyesterWeave.UnsignedIteratorEarlyStop(zero(PolyesterWeave.worker_type()), 0x00000000),)
      end
    end
    var"#avx#call#args#" = $avxcall_args
    $_turbo_call_
    var"##do#thread##" || $retexpr
    var"#thread#id#" = 0x00000000
    for var"#i#" ∈ eachindex(var"#threads#tuple#")
      var"#threads#" = var"#threads#tuple#"[var"#i#"]

      var"#thread#mask#" = PolyesterWeave.mask(var"#threads#")
      var"#threads#remain#" = (length(var"#threads#tuple#") == 1) || (var"#thread#mask#" ≠ zero(var"#thread#mask#"))
      while var"#threads#remain#"
        VectorizationBase.assume(var"#thread#mask#" ≠ zero(var"#thread#mask#"))
        var"#trailzing#zeros#" = vadd_nw(Base.trailing_zeros(var"#thread#mask#") % UInt32, 0x00000001)
        var"#thread#mask#" >>>= var"#trailzing#zeros#"
        var"#thread#id#" = vadd_nw(var"#thread#id#", var"#trailzing#zeros#")
        var"#thread#ptr#" = ThreadingUtilities.taskpointer(var"#thread#id#")
        ThreadingUtilities.wait(var"#thread#ptr#", var"#thread#id#")
        $update_return_values
        var"#threads#remain#" = var"#thread#mask#" ≠ 0x00000000
      end
      PolyesterWeave.free_threads!(var"#torelease#tuple#"[var"#i#"])
    end
    $retexpr
  end
  # Expr(:block, Expr(:meta,:inline), ls.preamble, q)
  Expr(:block, ls.preamble, q)
end
function define_vthread_blocks(vloop, u₁loop, u₂loop, u₁, u₂, ntmax, tn)
  if tn == 0
    loopunrollname = Symbol("#num#unrolls#thread#0#")
    lhs = :((var"#thread#factor#0#", var"#thread#factor#1#"))
  else
    loopunrollname = Symbol("#num#unrolls#thread#1#")
    lhs = :((var"#thread#factor#1#", var"#thread#factor#0#"))
  end
    sntmax = staticexpr(ntmax % Int)
    if vloop === u₁loop
        :($lhs = _choose_num_blocks($loopunrollname, StaticInt{$u₁}(), var"#nthreads#", $sntmax))
    elseif vloop === u₂loop
        :($lhs = _choose_num_blocks($loopunrollname, StaticInt{$u₂}(), var"#nthreads#", $sntmax))
    else
        :($lhs = _choose_num_blocks($loopunrollname, StaticInt{1}(), var"#nthreads#", $sntmax))
    end
end
function define_thread_blocks(threadedloop1, threadedloop2, vloop, u₁loop, u₂loop, u₁, u₂, ntmax)
  if vloop === threadedloop1
    define_vthread_blocks(threadedloop1, u₁loop, u₂loop, u₁, u₂, ntmax, 0)
  elseif vloop === threadedloop2
    define_vthread_blocks(threadedloop2, u₁loop, u₂loop, u₁, u₂, ntmax, 1)
  else
    :((var"#thread#factor#0#", var"#thread#factor#1#") = choose_num_blocks(var"#nthreads#", StaticInt{$(Int(ntmax))}()))
  end
end
function thread_two_loops_expr(
  ls::LoopSet, ua::UnrollArgs, valid_thread_loop::Vector{Bool}, ntmax::UInt, c::Float64,
  UNROLL::Tuple{Bool,Int8,Int8,Int8,Bool,Int,Int,Int,Int,Int,Int,Int,UInt}, OPS::Expr, ARF::Expr, AM::Expr, LPSYM::Expr
)
  looplen = looplengthprod(ls)
  # c = 0.0225 * c / looplen
  c = scale_cost(c, looplen)
  if all(isstaticloop, ls.loops)
    _num_threads = _choose_num_threads(c, ntmax, Int64(looplen))::UInt
    _num_threads > 1 || return avx_body(ls, UNROLL)
    choose_nthread = Expr(:(=), Symbol("#nthreads#"), _num_threads)
  else
    choose_nthread = :(_choose_num_threads($(Float32(c)), $ntmax))
    push_loop_length_expr!(choose_nthread, ls)
    choose_nthread = Expr(:(=), Symbol("#nthreads#"), choose_nthread)
  end
  threadedid1 = threadedid2 = 0
  for (i,v) ∈ enumerate(valid_thread_loop)
    v || continue
    if threadedid1 == 0
      threadedid1 = i
    else#if threadedid2 == 0
      threadedid2 = i
      break
    end
  end
  @unpack u₁loop, u₂loop, vloop, u₁, u₂max = ua
  u₂ = u₂max
  W = ls.vector_width
  threadedloop1 = getloop(ls, threadedid1)
  threadedloop2 = getloop(ls, threadedid2)
  define_len1, define_num_unrolls1, loopstart1, iterstop1, looprange1, lastrange1 = thread_loop_summary!(ls, ua, threadedloop1, false)
  define_len2, define_num_unrolls2, loopstart2, iterstop2, looprange2, lastrange2 = thread_loop_summary!(ls, ua, threadedloop2, true)
  loopboundexpr = Expr(:tuple)
  lastboundexpr = Expr(:tuple)
  for (i,loop) ∈ enumerate(ls.loops)
    if loop === threadedloop1
      push!(loopboundexpr.args, looprange1)
      push!(lastboundexpr.args, lastrange1)
    elseif loop === threadedloop2
      push!(loopboundexpr.args, looprange2)
      push!(lastboundexpr.args, lastrange2)
    else
      loop_boundary!(loopboundexpr, ls, loop, true)
      loop_boundary!(lastboundexpr, ls, loop, true)
    end
  end
  avxcall_args = Expr(:tuple, lastboundexpr, Symbol("#vargs#"))
  _turbo_call_ = :(_turbo_!(Val{$UNROLL}(), $OPS, $ARF, $AM, $LPSYM, Val(typeof(var"#avx#call#args#")), flatten_to_tuple(var"#avx#call#args#")...))
  # _turbo_orig_ = :(_turbo_!(Val{$UNROLL}(), $OPS, $ARF, $AM, $LPSYM, var"#lv#tuple#args#"))
  update_return_values = if length(ls.outer_reductions) > 0
    retv = loopset_return_value(ls, Val(false))
    _turbo_call_ = Expr(:(=), retv, _turbo_call_)
    # _turbo_orig_ = Expr(:(=), retv, _turbo_orig_)
    outer_reduct_combine_expressions(ls, retv)
  else
    nothing
  end
  blockdef = define_thread_blocks(threadedloop1, threadedloop2, vloop, u₁loop, u₂loop, u₁, u₂, ntmax)
  iterdef1 = define_block_size(threadedloop1, vloop, 0, ls.vector_width)
  iterdef2 = define_block_size(threadedloop2, vloop, 1, ls.vector_width)
  retexpr = length(ls.outer_reductions) > 0 ? :(return $retv) : :(return nothing)
  q = quote
    $choose_nthread # UInt
    # @show var"#nthreads#"
    $loopstart1
    $loopstart2
    var"##do#thread##" = var"#nthreads#" > one(var"#nthreads#")
    if var"##do#thread##"
      # if var"#nthreads#" ≤ 1
      #   $_turbo_orig_
      #   return $retexpr
      # end
      $define_len1
      $define_len2
      $define_num_unrolls1
      $define_num_unrolls2
      var"#unroll#prod#" = vmul_nw(var"#num#unrolls#thread#0#", var"#num#unrolls#thread#1#")
      if var"#nthreads#" ≥ var"#unroll#prod#"
        var"#nthreads#" = var"#unroll#prod#"
        var"#thread#factor#0#" = var"#num#unrolls#thread#0#"
        var"#thread#factor#1#" = var"#num#unrolls#thread#1#"
      else
        var"##thread#0##excess##" = var"#num#unrolls#thread#0#" ≥ var"#nthreads#"
        var"##thread#1##excess##" = var"#num#unrolls#thread#1#" ≥ var"#nthreads#"
        if var"##thread#0##excess##" & var"##thread#1##excess##"
          $blockdef
        elseif var"##thread#0##excess##" # var"#num#unrolls#thread#1#" is small but var"#num#unrolls#thread#0#" is not; we want to place a small one in front
          (var"#thread#factor#1#", var"#thread#factor#0#")  = _choose_num_blocks(var"#num#unrolls#thread#1#", StaticInt{1}(), var"#nthreads#", $(staticexpr(ntmax % Int)))
        else # var"#num#unrolls#thread#0#" is small, and var"#num#unrolls#thread#1#" may or may not be
          (var"#thread#factor#0#", var"#thread#factor#1#")  = _choose_num_blocks(var"#num#unrolls#thread#0#", StaticInt{1}(), var"#nthreads#", $(staticexpr(ntmax % Int)))
        end
        var"#thread#factor#0#" = min(var"#thread#factor#0#", var"#num#unrolls#thread#0#")
        var"#thread#factor#1#" = min(var"#thread#factor#1#", var"#num#unrolls#thread#1#")
      end
      # @show (var"#thread#factor#0#", var"#thread#factor#1#")
      var"#nrequest#" = vsub_nsw(vmul_nsw(var"#thread#factor#0#", var"#thread#factor#1#" % UInt32), 0x00000001)
      var"#loop#1#start#init#" = var"#iter#start#0#"
      var"##do#thread##" = var"#nrequest#" ≠ 0x00000000
      if var"##do#thread##"
        # var"#threads#tuple#", var"#torelease#tuple#" = PolyesterWeave.request_threads(Threads.threadid(), var"#nrequest#")
        var"#threads#tuple#", var"#torelease#tuple#" = PolyesterWeave.request_threads(var"#nrequest#")
        var"##do#thread##" = false
        $iterdef1
        $iterdef2
        var"#thread#launch#count#0#" = 0x00000000
        var"#thread#launch#count#1#" = 0x00000000
        var"#thread#id#" = 0x00000000
        for var"#threads#" ∈ var"#threads#tuple#"
          var"#nrequest#" = length(var"#threads#")
          # @show var"#base#block#size#thread#0#", var"#block#rem#step#0#" var"#base#block#size#thread#1#", var"#block#rem#step#1#"
          var"#thread#launch#count#" = 0x00000000
          var"#thread#mask#" = PolyesterWeave.mask(var"#threads#")
          var"#threads#remain#" = var"#thread#launch#count#" ≠ var"#nrequest#"
          var"##do#thread##" |= var"#threads#remain#"
          while var"#threads#remain#"
            VectorizationBase.assume(var"#thread#mask#" ≠ zero(var"#thread#mask#"))
            var"#trailzing#zeros#" = Base.trailing_zeros(var"#thread#mask#") % UInt32
            var"#nblock#size#thread#0#" = Core.ifelse(
              var"#thread#launch#count#0#" < (var"#nrem#thread#0#" % UInt32),
              vadd_nw(var"#base#block#size#thread#0#", var"#block#rem#step#0#"),
              var"#base#block#size#thread#0#"
            )
            var"#nblock#size#thread#1#" = Core.ifelse(
              var"#thread#launch#count#1#" < (var"#nrem#thread#1#" % UInt32),
              vadd_nw(var"#base#block#size#thread#1#", var"#block#rem#step#1#"),
              var"#base#block#size#thread#1#"
            )
            var"#trailzing#zeros#" = vadd_nw(var"#trailzing#zeros#", 0x00000001)
            $iterstop1
            $iterstop2
            var"#thread#id#" = vadd_nw(var"#thread#id#", var"#trailzing#zeros#")
            # @show var"#thread#id#" $loopboundexpr
            var"##lbvargs#to_launch##" = ($loopboundexpr, var"#vargs#")
            avx_launch(Val{$UNROLL}(), $OPS, $ARF, $AM, $LPSYM, StaticType{typeof(var"##lbvargs#to_launch##")}(), flatten_to_tuple(var"##lbvargs#to_launch##"), var"#thread#id#")
            var"#thread#mask#" >>>= var"#trailzing#zeros#"

            var"##end#inner##" = var"#thread#launch#count#0#" == vsub_nw(var"#thread#factor#0#", 0x00000001)
            var"#thread#launch#count#0#" = Core.ifelse(var"##end#inner##", 0x00000000, vadd_nw(var"#thread#launch#count#0#", 0x00000001))
            var"#thread#launch#count#1#" = Core.ifelse(var"##end#inner##", var"#thread#launch#count#1#" + 0x00000001, var"#thread#launch#count#1#")

            var"#iter#start#0#" = Core.ifelse(var"##end#inner##", var"#loop#1#start#init#", var"#iter#stop#0#")
            var"#iter#start#1#" = Core.ifelse(var"##end#inner##", var"#iter#stop#1#", var"#iter#start#1#")

            var"#thread#launch#count#" = vadd_nw(var"#thread#launch#count#", 0x00000001)
            var"#threads#remain#" = var"#thread#launch#count#" ≠ var"#nrequest#"
          end
        end
      else# eliminate undef var errors that the compiler should be able to figure out are unreachable, but doesn't
        var"#torelease#tuple#" = (zero(PolyesterWeave.worker_type()),)
        var"#threads#tuple#" = PolyesterWeave.UnsignedIteratorEarlyStop(zero(PolyesterWeave.worker_type()), 0x00000000)
      end
    end
    # @show $lastboundexpr
    var"#avx#call#args#" = $avxcall_args
    $_turbo_call_
    var"##do#thread##" || $retexpr
    var"#thread#id#" = 0x00000000
    for var"#i#" ∈ eachindex(var"#threads#tuple#")
      var"#threads#" = var"#threads#tuple#"[var"#i#"]
      var"#thread#mask#" = PolyesterWeave.mask(var"#threads#")
      var"#threads#remain#" = (length(var"#threads#tuple#") == 1) || (var"#thread#mask#" ≠ zero(var"#thread#mask#"))
      while var"#threads#remain#"
        VectorizationBase.assume(var"#thread#mask#" ≠ zero(var"#thread#mask#"))
        var"#trailzing#zeros#" = vadd_nw(Base.trailing_zeros(var"#thread#mask#") % UInt32, 0x00000001)
        var"#thread#mask#" >>>= var"#trailzing#zeros#"
        var"#thread#id#" = vadd_nw(var"#thread#id#", var"#trailzing#zeros#")
        var"#thread#ptr#" = ThreadingUtilities.taskpointer(var"#thread#id#")
        ThreadingUtilities.wait(var"#thread#ptr#", var"#thread#id#")
        $update_return_values
        var"#threads#remain#" = var"#thread#mask#" ≠ 0x00000000
      end
      PolyesterWeave.free_threads!(var"#torelease#tuple#"[var"#i#"])
    end
    $retexpr
  end
  # Expr(:block, Expr(:meta,:inline), ls.preamble, q)
  Expr(:block, ls.preamble, q)
end

function valid_thread_loops(ls::LoopSet)
  order, u₁loop, u₂loop, vectorized, u₁, u₂, c, shouldinline = choose_order_cost(ls)
  # NOTE: `names` are being placed in the opposite order here versus normal lowering!
  copyto!(names(ls), order); init_loop_map!(ls)
  u₁loop = getloop(ls, u₁loop)
  _u₂loop = getloopid_or_nothing(ls, u₂loop)
  u₂loop = _u₂loop === nothing ? u₁loop : getloop_from_id(ls, _u₂loop)
  ua = UnrollArgs(u₁loop, u₂loop, getloop(ls, vectorized), u₁, u₂, u₂)
  valid_thread_loop = fill(true, length(order))
  for op ∈ operations(ls)
    if isstore(op) && (length(reduceddependencies(op)) > 0)
      for reduceddep ∈ reduceddependencies(op)
        for (i,o) ∈ enumerate(order)
          if o === reduceddep
            valid_thread_loop[i] = false
          end
        end
      end
    end
  end
  for (i,o) ∈ enumerate(order)
    loop = getloop(ls, o)
    if isstaticloop(loop) & (length(loop) ≤ 1)
      valid_thread_loop[i] = false
    end
  end
  valid_thread_loop, ua, c
end
function avx_threads_expr(
    ls::LoopSet, UNROLL::Tuple{Bool,Int8,Int8,Int8,Bool,Int,Int,Int,Int,Int,Int,Int,UInt},
    nt::UInt, OPS::Expr, ARF::Expr, AM::Expr, LPSYM::Expr
)
    valid_thread_loop, ua, c = valid_thread_loops(ls)
    num_candiates = sum(valid_thread_loop)
    # num_to_thread = min(num_candiates, 2)
    # candidate_ids =
    if (num_candiates == 0) || (nt ≤ 1) # it was called from `avx_body` but now `nt` was set to `1`
        avx_body(ls, UNROLL)
    elseif (num_candiates == 1) || (nt ≤ 3)
        thread_one_loops_expr(ls, ua, valid_thread_loop, nt, c, UNROLL, OPS, ARF, AM, LPSYM)
    else # requires at least 4 threads
        thread_two_loops_expr(ls, ua, valid_thread_loop, nt, c, UNROLL, OPS, ARF, AM, LPSYM)
    end
end


text:
@testset "readrequestline!" begin
    io = IOBuffer(httptest1)
    c = Connection(io)

    Spirit.readrequestline!(c)

    @test c.method == "POST"
    @test c.uri == "/"
    @test c.httpversion == "1.1"
    close(io)

    io = IOBuffer(httptest2)
    c = Connection(io)

    Spirit.readrequestline!(c)

    @test c.method == "GET"
    @test c.uri == "/docs/index.html"
    @test c.httpversion == "1.1"
    close(io)
end

@testset "readheaders!" begin
    io = IOBuffer(httptest1)
    c = Connection(io)

    Spirit.readrequestline!(c)
    Spirit.readheaders!(c)

    rawheaders = c.rawheaders
    headers = c.headers
    @test length(rawheaders) == 2
    @test length(headers) == 2
    @test rawheaders[1] == ("Host" => "www.somewhere.com")
    @test rawheaders[2] == ("Content-Length" => "9")

    close(io)

    io = IOBuffer(httptest2)
    c = Connection(io)

    Spirit.readrequestline!(c)
    Spirit.readheaders!(c)

    rawheaders = c.rawheaders
    headers = c.headers
    @test length(rawheaders) == 5
    @test length(headers) == 5
    @test rawheaders[1] == ("Host" => "www.nowhere123.com")
    @test rawheaders[2] == ("Accept" => "image/gif, image/jpeg, */*")
    @test rawheaders[3] == ("Accept-Language" => "en-us")
    @test rawheaders[4] == ("Accept-Encoding" => "gzip, deflate")
    @test rawheaders[5] == ("User-Agent" => "Mozilla/4.0 (compatible; MSIE 6.0; Windows NT 5.1)")

    close(io)
end

@testset "createbodystream!" begin

    # with content-length
    io = IOBuffer(httptest1)
    c = Connection(io)

    Spirit.readrequestline!(c)
    Spirit.readheaders!(c)
    Spirit.createbodystream!(c)

    @test isa(c.bodystream, SizedStream)
    @test c.bodystream.size == 9
    @test c.bodystartedread == false
    @test c.bodyfullyread == false
    b = read(c.bodystream, UInt8)
    @test b == UInt8('B')
    @test c.bodystartedread == true
    @test c.bodyfullyread == false
    str = read(c.bodystream)
    @test c.bodystartedread == true
    @test c.bodyfullyread == true
    @test String(str) == "ody Body"
    close(io)

    # transfer encoding
    io = IOBuffer(httptest3)
    c = Connection(io)

    Spirit.readrequestline!(c)
    Spirit.readheaders!(c)
    Spirit.createbodystream!(c)

    @test isa(c.bodystream, ChunkedTransferDecodeStream)
    @test c.bodystartedread == false
    @test c.bodyfullyread == false
    b = read(c.bodystream, UInt8)
    @test b == UInt8('\'')
    @test c.bodystartedread == true
    @test c.bodyfullyread == false
    str = read(c.bodystream)
    @test c.bodystartedread == true
    @test c.bodyfullyread == true
    expected = "Twas brillig, and the slithy toves\n" *
        "      Did gyre and gimble in the wabe:\n" *
        "All mimsy were the borogoves,\n" *
        "      And the mome raths outgrabe."
    @test String(str) == expected

    rawtrailers = c.rawtrailers
    @test length(rawtrailers) == 2
    @test rawtrailers[1] == ("Happiness" => "100")
    @test rawtrailers[2] == ("Author" => "Lewis Carroll")

    close(io)
end

@testset "processrequest!" begin
    io = IOBuffer(httptest1)
    c = Connection(io)
    Spirit.processrequest!(c)
    close(io)
end

@testset "writeresponse!" begin
    io = IOBuffer()
    c = Connection(io)

    body = "Safety First"
    bodyio = IOBuffer(body)
    headers = Spirit.Headers("Content-Length" => string(length(body)))
    Spirit.writeresponse!(c, 200, "OK", headers, bodyio)

    response = take!(io)

    expected = "HTTP/1.1 200 OK\r\n" *
        "Content-Length: 12\r\n" *
        "\r\n" *
        "Safety First"

    @test expected == String(response)
end


text:
using BSON
using PyCall
using StatsBase
using Statistics
using OrderedCollections

using .Evaluation: compute_stats, _get_anomaly_class

const BAYES_CACHE = "bayes_cache.bson"
"""
	BayesianHyperOpt

Wrapper around `skotp.optimizer.Optimizer`.
Uses default params from `skopt.gp_minimize`.
Constructed using a list/array of `skotp.space`'s classes.
"""
struct BayesianHyperOpt
	opt
	function BayesianHyperOpt(dimensions)
		py"""
		import numpy as np
		from skopt import Optimizer
		from skopt.utils import cook_estimator, normalize_dimensions
		from sklearn.utils import check_random_state

		def construct_optimizer(dimensions):
			rng = check_random_state(7) 			# fix seed to fit always the same model
			space = normalize_dimensions(dimensions)

			base_estimator = cook_estimator(
							"GP", space=space,
							random_state=rng.randint(0, np.iinfo(np.int32).max),
							noise="gaussian")

			# default from gp_minimize
			acq_optimizer_kwargs = {"n_points": 10000, "n_restarts_optimizer": 5, "n_jobs": 1}
			acq_func_kwargs = {"xi": 0.01, "kappa": 1.96}

			optimizer = Optimizer(
				space,
				base_estimator,
				n_initial_points=0, 				# should be zero in order to skip the random generation
				initial_point_generator="random", 	# it fallbacks into random even if it is None, so here you go
				n_jobs=-1,
				acq_func="gp_hedge", 				# default from gp_minimize
				acq_optimizer="auto", 				# default from gp_minimize
				random_state=rng,
				model_queue_size=None,				# default from gp_minimize
				acq_optimizer_kwargs=acq_optimizer_kwargs,
				acq_func_kwargs=acq_func_kwargs)

			return optimizer
		"""
		new(py"construct_optimizer"(dimensions))
	end
end

"""
	function tell(bho::BayesianHyperOpt, x0, y0)

Fits the optimizer based on multiple observations of x0 (hyperparameter)
and y0 (objective). Be careful about the input types here as conversions
to Python are tricky. By default `x0` is passed as list of lists, however
vector of values `y0` is converted to `numpy.ndarray`, which cannot
be processed by `skopt` (!! What a nice ecosystem !!).
"""
function tell(bho::BayesianHyperOpt, x0, y0)
	bho.opt.tell(x0, Tuple(y0))
end

function ask(bho::BayesianHyperOpt)
	bho.opt.ask()
end

"""
	skopt_parse(dimension, entry::T)
Converts and entry of a hyperparameter from skotp dimension to named tuple.
Dimension's name have to following this convention:

```py
	Integer(1, 3, name="log2_{name}") # will convert to 2^p
	Integer(1, 3, name="log10_{name}") # will convert to 10^p
	Categorical(categories = [:dense, :all, :none], name="sym_{name}") # will convert to and from Symbol
	Categorical(categories = [identity], name="fun_{name}") # will convert to and from a function name
```
Other cases of `T` such as Bool, Real, Int are handled without any specific treatment.
"""
skopt_parse(dimension, entry::Bool) = (;[Symbol(dimension.name) => entry]...)
skopt_parse(dimension, entry::T) where {T <: AbstractFloat} = (;[Symbol(dimension.name) => entry]...)

function skopt_parse(dimension, entry::T) where {T <: Integer}
	if startswith(dimension.name, "log2_")
		return (;[Symbol(dimension.name[6:end]) => 2^entry]...)
	elseif startswith(dimension.name, "log10_")
		return (;[Symbol(dimension.name[7:end]) => 10^entry]...)
	else
		return (;[Symbol(dimension.name) => entry]...)
	end
end

function skopt_parse(dimension, entry::T) where {T <: AbstractString}
	if startswith(dimension.name, "sym_")
		return (;[Symbol(dimension.name[5:end]) => Symbol(entry)]...)
	elseif startswith(dimension.name, "fun_")
		return (;[Symbol(dimension.name[5:end]) => eval(:($(Symbol(entry))))]...)
	else
		return (;[Symbol(dimension.name) => entry]...)
	end
end

"""
	from_skopt(space, p)
Converts an array of hyperparameters `p` from skotp to a named tuple
based on provided metadata in skotp's Space class instance `space`.
The order of entries in `p` coincides with the order of dimensions in `space`,
if that is not the case something broke down.
"""
function from_skopt(space, p)
	merge([skopt_parse(dim, entry) for (dim, entry) in zip(space, p)]...)
end

tuple_parse(dimension, entry::Bool) = entry
tuple_parse(dimension, entry::T) where {T <: AbstractString} = entry
tuple_parse(dimension, entry::T) where {T <: AbstractFloat} = entry
tuple_parse(dimension, entry::Symbol) = string(entry)
tuple_parse(dimension, entry::Function) = string(entry) 	# SPTN has identity function in hyperparameters

function tuple_parse(dimension, entry::T) where {T <: Integer}
	if startswith(dimension.name, "log2")
		return Int(log2(entry))
	elseif startswith(dimension.name, "log10")
		return Int(log10(entry))
	else
		return entry
	end
end

"""
	to_skopt(space, p)
Converts named tuple of hyperparameters `p` from Julia to a list of hyperparameters in skopt
based on provided metadata in skotp's Space class instance `space`.
The order of entries in `p` does not have to coincide with the order of dimensions in `space`,
but it will be imposed by keys of the `space` tuple.
"""
function to_skopt(space, p)
	[tuple_parse(dim, p[name]) for (name, dim) in pairs(space)]
end

"""
	objective_value(runs;; run_agg=mean, score_agg=max, fail_value=0.0)

Validates entry from bayesian cache and returns aggregated score, which is optimized by `BayesianHyperOpt`.
For empty entries (no model trained) we return `fail_value`. In case of images `run_agg` allows to set
the aggregation over anomaly_classes to be also maximum.
"""
function objective_value(runs; run_agg=mean, score_agg=maximum, fail_value=0.0)
	if length(runs) > 0 # should we be more strict?
		return run_agg([score_agg(v) for v in values(runs)])
	else
		return fail_value
	end
end

"""
	function bayes_params(space, folder, fallback; add_model_seed=false, kwargs...)
Tries to load a cache of results to seed the optimizer with. The file should be named by `BAYES_CACHE`
global constant and located in `folder`. If the cache has not been created, the method calls `fallback`.
Additional keyword arguments are passed into `objective_value` function allowing to modify it's behavior.
"""
function bayes_params(space, folder, fallback; add_model_seed=false, kwargs...)
	cache = load_bayes_cache(folder)
	if cache !== nothing
		x0 = [v[:parameters] for v in values(cache)]
		y0 = [objective_value(v[:runs], kwargs...) for v in values(cache)]
		@info("Loaded $(length(cache)) cached results from $(folder).")
		# conversion to tuples this makes skopt happy
		# as it cannot parse list of numpy arrays as valid "2D" input
		x0 = [(to_skopt(space, x)..., ) for x in x0]
		opt = BayesianHyperOpt(collect(space))
		tell(opt, x0, y0)
		@info("Fitted GP with $(length(x0)) samples.")
		p = from_skopt(space, ask(opt))

		# registers only parameters with empty arrays and store cache
		register_run!(cache, (;parameters=p); fit_results=false)
		@info("Updated cache size: $(length(cache)). Saving.")
		save_bayes_cache(folder, cache)

		@info("Fetched new hyperparameters: $(p).")
		return add_model_seed ? merge((;init_seed=rand(1:Int(1e8))), p) : p
	else
		@warn("Could not find bayesian cache file at $(folder), calling fallback.")
		return fallback()
	end
end

"""
	load_bayes_cache(folder)

Tries to load serialized results of Bayesian optimization from `folder`.
"""
function load_bayes_cache(folder)
    file = joinpath(folder, BAYES_CACHE)
	isfile(file) ? BSON.load(file) : OrderedDict{UInt64, Any}()
end

"""
	save_bayes_cache(folder)

Saves serialized results of Bayesian optimization to `folder`.
"""
function save_bayes_cache(folder, cache)
    file = joinpath(folder, BAYES_CACHE)
	wsave(file, cache)
end

"""
	update_bayes_cache(folder, results; kwargs...)

Loads bayesian cache from folder and updates it with results.
Additional keyword arguments are passed into `register_run!`.
"""
function update_bayes_cache(folder, results; kwargs...)
	cache = load_bayes_cache(folder)
	@info("Loaded cache size: $(length(cache)). Updating.")
	for r in results
		register_run!(cache, r; kwargs...)
	end
	@info("Updated cache size: $(length(cache)). Saving.")
	save_bayes_cache(folder, cache)
end

"""
	register_run!(cache, r; fit_results=true, metric=:val_auc, flip_sign=true, ignore=Set([:init_seed]))

Updates `cache` with new result from `r`(named tuple or symbol indexed dictionary).
All basic metrics from `.Evaluation` module are supported and their column name
can be passed in the `metric` argument.
Set of hyperparameters that are not optimized is given by the `ignore` argument,
which by default filters `:init_seed`.
If `fit_results` is false register_run! will write an empty entry under a given parameters. This is used
to indicate that a traing has begun for such parameters and if there are not filled, we assume that it failed.
Flag `flip_sign` is used for aligning sign of the metric with minimization objective.
"""
function register_run!(cache, r; fit_results=true, metric=:val_auc, flip_sign=true, ignore=Set([:init_seed]))
	metric_value = try
		fit_results ? compute_stats(r, top_metrics=false)[metric] : []
	catch e
		@warn("Cache has not been updated. Stats computation failed due to $e")
		return
	end

	# compute hash from all parameters (used for filtering during evaluation)
	phash = hash(r[:parameters])
	# filter only those parameters that are being optimized
	parameters = (;[p for p in pairs(r[:parameters]) if !(p[1] in ignore)]...)
	ophash = hash(parameters)   # corresponds only to optimizable parameters
								# there may be colisions (e.g. parameter being run with two seeds)

	if ophash in keys(cache) && fit_results				# add to existing entry
		@info("Adding run to existing entry $(repr(ophash)) - $(parameters).")
		entry = cache[ophash]
		seed = r[:seed]
		anomaly_class = _get_anomaly_class(r)

		runs = entry[:runs]
		phashes = vcat(entry[:phashes], phash)

		key = (seed, anomaly_class)
		if key in keys(runs)	# this is true if for some `ophash` we have multiple outputs (scores)
			runs[key] = vcat(runs[key], flip_sign ? -metric_value : metric_value)
		else
			runs[key] = flip_sign ? -metric_value : metric_value
		end

		cache[ophash] = merge(entry, (;runs=runs, phashes=phashes))
	elseif !(ophash in keys(cache)) && fit_results		# or create new entry if there are results (mainly used during manual cache creation)
		@info("Creating new non-empty entry $(repr(ophash)) - $(parameters).")
		seed = r[:seed]
		anomaly_class = _get_anomaly_class(r)
		runs = Dict((seed, anomaly_class) => flip_sign ? -metric_value : metric_value)
		phashes = [phash]

		cache[ophash] = (;
			parameters = parameters,
			runs = runs,
			phashes=phashes)
	elseif !(ophash in keys(cache)) && fit_results == false 	# add an empty entry (this is to indicate that such parameters are about to be trained)
		@info("Creating new empty entry $(repr(ophash)) - $(parameters).")
		entry = (;
			parameters = parameters,
			runs = Dict{Tuple{Int,Int}, Any}(),
			phashes = UInt64[])

		cache[ophash] = entry
	else 		# otherwise warn the user that the chosen parameters already exist in the cache
		@warn("Entry $(repr(ophash)) - $(parameters) already exist in the cache.")
	end
end


text:
## Exercise 9-1
## Write a program that reads words.txt and prints only the words with more than 20 characters (not counting whitespace).

println("Ans: ")

for line in eachline("words.txt")
    str = replace(strip(line), " " => "")
    if length(str) > 20
        println(str)
    end
end

println("End.")


text:
module ColumnPartitions

export column_partition, vector_partition, map_elementwise

using ArgCheck: @argcheck

####
#### schemas
####

abstract type AbstractSchema end

function promote_schema_and_eltype(a, b)
    sch_a, T_a = a
    sch_b, T_b = b
    (promote_schema(sch_a, sch_b), promote_type(T_a, T_b))
end

###
### AtomSchema
###

struct AtomSchema <: AbstractSchema end

function schema_and_eltype(x::T) where T
    AtomSchema(), T
end

Base.length(sch::AtomSchema) = 1

flat_get(column, offset, sch::AtomSchema) = column[offset + 1]

function flat_set!(column, offset, sch::AtomSchema, x)
    column[offset + 1] = x
end

promote_schema(::AtomSchema, ::AtomSchema) = AtomSchema()

eltype_schema(::AtomSchema, T) = T

###
### ArraySchema
###

struct ArraySchema{S} <: AbstractSchema end       # FIXME check S

Base.size(sch::ArraySchema{S}) where S = S

Base.length(sch::ArraySchema) = prod(size(sch))

function schema_and_eltype(A::AbstractArray)
    @argcheck !Base.has_offset_axes(A)
    ArraySchema{size(A)}(), eltype(A)
end

function flat_get(column, offset, sch::ArraySchema)
    reshape(view(column, offset .+ (1:length(sch))), size(sch))
end

function flat_set!(column, offset, sch::ArraySchema, A::AbstractArray)
    @argcheck size(A) == size(sch)
    column[offset .+ (1:length(sch))] = vec(A)
end

function promote_schema(a::ArraySchema, b::ArraySchema)
    @argcheck size(a) == size(b)
    a
end

eltype_schema(sch::ArraySchema, T) = Array{T, length(size(sch))}

###
### TupleSchema
###

struct TupleSchema{T <: Union{Tuple,NamedTuple}} <: AbstractSchema
    sub_schemas::T
end

function schema_and_eltype(tup::Union{Tuple,NamedTuple})
    sch_and_Ts = map(schema_and_eltype, tup)
    TupleSchema(map(first, sch_and_Ts)), mapreduce(last, promote_type, sch_and_Ts)
end

Base.length(sch::TupleSchema) = sum(length, sch.sub_schemas)

function flat_get(column, offset, sch::TupleSchema)
    map(sch.sub_schemas) do sub_schema
        ret = flat_get(column, offset, sub_schema)
        offset += length(sub_schema)
        ret
    end
end

function flat_set!(column, offset, sch::TupleSchema{S}, tup::T
                   ) where {S, T <: Union{Tuple, NamedTuple}}
    @argcheck fieldnames(S) == fieldnames(T)
    map(sch.sub_schemas, tup) do sub_schema, elt
        ret = flat_set!(column, offset, sub_schema, elt)
        offset += length(sub_schema)
        ret
    end
end

function promote_schema(a::TupleSchema{S}, b::TupleSchema{T}) where {S,T}
    @argcheck fieldnames(S) == fieldnames(T) "incompatible schemas"
    TupleSchema(map(promote_schema, a.sub_schemas, b.sub_schemas))
end

function eltype_schema(sch::TupleSchema{<:Tuple}, T)
    Tuple{map(s -> eltype_schema(s, T), sch.sub_schemas)...}
end

function eltype_schema(sch::TupleSchema{S}, T) where {N, S <: NamedTuple{N}}
    NamedTuple{N, Tuple{map(s -> eltype_schema(s, T), values(sch.sub_schemas))...}}
end

_unflatten(sch, v) = flat_get(v, 0, sch)

####
#### API
####

###
### Column partitioned matrix
###

struct ColumnPartition{T, S <: AbstractSchema, M <: AbstractMatrix} <: AbstractVector{T}
    sch::S
    matrix::M
end

function ColumnPartition(sch::S, matrix::M) where {M,S}
    ColumnPartition{eltype_schema(sch, eltype(M)),S,M}(sch, matrix)
end

schema(p::ColumnPartition) = p.sch

function column_partition(v::AbstractVector)
    @argcheck !Base.has_offset_axes(v)
    @argcheck !isempty(v)
    sch, T = mapreduce(schema_and_eltype, promote_schema_and_eltype, v)
    matrix = Matrix{T}(undef, length(sch), length(v))
    for (column, elt) in zip(eachcol(matrix), v)
        flat_set!(column, 0, sch, elt)
    end
    ColumnPartition(sch, matrix)
end

Base.size(p::ColumnPartition) = (size(p.matrix, 2), )

Base.getindex(p::ColumnPartition, i) = _unflatten(p.sch, @view p.matrix[:, i])

###
### Partitioned vector
###

struct VectorPartition{S <: AbstractSchema, V <: AbstractVector}
    sch::S
    vector::V
end

Base.getindex(p::VectorPartition) = _unflatten(p.sch, p.vector)

####
#### mapping
####

function schema_and_flatten(x)
    sch, T = schema_and_eltype(x)
    @show sch
    v = Vector{T}(undef, length(sch))
    flat_set!(v, 0, sch, x)
    sch, @show v
end

# FIXME replace by broadcast?
function map_elementwise(f, p_first::ColumnPartition, p_rest::ColumnPartition...)
    sch = p_first.sch
    @assert all(p -> schema(p) == sch, p_rest)
    VectorPartition(sch, map(f, eachrow(p_first.matrix), map(p -> eachrow(p.matrix), p_rest)...))
end

function map_elementwise(f, p_first::VectorPartition, p_rest::VectorPartition...)
    sch = p_first.sch
    @assert all(p -> schema(p) == sch, p_rest)
    VectorPartition(sch, map(f, p_first.vector, map(p -> p.vector, p_rest)...))
end

end # module


text:
using Images
using ImageView

function float2img(a::AbstractArray{T, 3}) where {T <: Real}
    img = permutedims(a, (3, 1, 2))
    img = colorview(RGB, img)
    return img
end

function img2float(img::Array{RGB{T}}) where T
    a = channelview(img)
    a = rawview(a)
    a = permutedims(a, (2,3,1))
    a = a / 255
    return a
end

img2float(img::Array{Gray{T}}) where T = img2float(RGB.(img))

function load_image(path::AbstractString, transform=nothing)
    img = load(path)
    if transform != nothing
        img = transform(img)
    end
    return img2float(img)
end

function view_image(a::AbstractArray{<:Real, 3})
    img = float2img(a)
    imshow(img)
end

struct ImageFolder
    dirpath::AbstractString
    paths::Vector{String}
    labels::Vector{Int}
    classes::Vector{String}
    transform::Function
end

function ImageFolder(dirpath::AbstractString; transform::Function=identity)
    paths = String[]
    classes = String[]
    labels = Int[]
    for subdir in readdir(dirpath)
        subdir_path = joinpath(dirpath, subdir)
        for filename in readdir(subdir_path)
            path = joinpath(subdir_path, filename)
            push!(paths, path)
            if isempty(classes) || classes[end] != subdir
                push!(classes, subdir)
            end
            push!(labels, length(classes))
        end
    end
    return ImageFolder(dirpath, paths, labels, classes, transform)
end

Base.show(io::IO, dataset::ImageFolder) = print(io, "ImageFolder(\"$(dataset.dirpath)\")")
Base.length(dataset::ImageFolder) = length(dataset.labels)

function Base.getindex(dataset::ImageFolder, idxs)
    labels = dataset.labels[idxs]
    imgs = Vector{Array}(undef, length(idxs))
    # Threads.@threads - causes concurrency violation exception from FileIO
    for (i, idx) in collect(enumerate(idxs))
        img = load(dataset.paths[idx])
        img = dataset.transform(img)
        arr = img2float(img)
        imgs[i] = arr
    end
    @assert(all(sz -> sz == size(imgs[1]), [size(img) for img in imgs]),
            "Cannot concat images into a batch since they have different size. " *
            "Use `transform` keyword argument to ImageFolder constructor to resize " *
            "all images to the same size.")
    imgs4d = cat([reshape(img, size(img)..., 1) for img in imgs]...; dims=4)
    return imgs4d, labels
end

function Avalon.getbatch(dataset::ImageFolder, i::Int, sz::Int)
    start = (i-1)*sz + 1
    finish = min(i*sz, length(dataset))
    start > finish && return nothing
    return dataset[start:finish]
end


text:
"""
    vphaseP(θ, VP0, VS0, ϵ, δ)

Compute P wave phase velocity along phase angle θ (with respect to axis of symmetry) based on Vp₀ and Vs₀ (P and S wave speed along axis of symmetry) and Thomsen parameters ϵ and δ.
"""
function vphaseP(θ,VP0,VS0,ϵ,δ)
    f = 1 - (VS0^2)/(VP0^2)
    D = sqrt(1 + 4*sin(θ)^2/f*(2*δ*cos(θ)^2 - ϵ*cos(2θ)) +
             (2*ϵ*sin(θ)^2/f)^2)
    vsq = 1 + ϵ*sin(θ)^2 - f/2 + D*f/2
    return sqrt(vsq)*VP0
end

"""
    vphaseSV(θ, VP0, VS0, ϵ, δ)

Compute SV wave phase velocity along phase angle θ (with respect to axis of symmetry) based on Vp₀ and Vs₀ (P and S wave speed along axis of symmetry) and Thomsen parameters ϵ and δ.
"""
function vphaseSV(θ,VP0,VS0,ϵ,δ)
    f = 1 - (VS0^2)/(VP0^2)
    D = sqrt(1 + 4*sin(θ)^2/f*(2*δ*cos(θ)^2 - ϵ*cos(2θ)) +
             (2*ϵ*sin(θ)^2/f)^2)
    vsq = 1 + ϵ*sin(θ)^2 - f/2 - D*f/2
    return sqrt(vsq)*VP0
end

"""
    vphaseSH(θ, VP0, VS0, γ)

Compute SH wave phase velocity along phase angle θ (with respect to axis of symmetry) based on Vp₀ and Vs₀ (P and S wave speed along axis of symmetry) and Thomsen parameter γ.
"""
function vphaseSH(θ,VS0,γ)
    return VS0*sqrt(1 + 2*γ*sin(θ)^2)
end

"""
    vphaseP(θ, c::AbstractMatrix, ρ)

Compute P wave phase velocity along phase angle θ (with respect to axis of symmetry) based on stiffness matrix c (in Voight notation) and density ρ.
"""
function vphaseP(θ,c::AbstractMatrix,ρ)
    stsq = sin(θ)^2
    ctsq = cos(θ)^2
    D = sqrt(((c[1,1]-c[5,5])*stsq - (c[3,3]-c[5,5])*ctsq)^2 +
             4*(c[1,3]+c[5,5])^2*stsq*ctsq)
    rv = (c[1,1] + c[5,5])*stsq + (c[3,3]+c[5,5])*ctsq + D
    return sqrt(rv/(2ρ))
end

"""
    vphaseSV(θ, c::AbstractMatrix, ρ)

Compute SV wave phase velocity along phase angle θ (with respect to axis of symmetry) based on stiffness matrix c (in Voight notation) and density ρ.
"""
function vphaseSV(θ,c::AbstractMatrix,ρ)
    stsq = sin(θ)^2
    ctsq = cos(θ)^2
    D = sqrt(((c[1,1]-c[5,5])*stsq - (c[3,3]-c[5,5])*ctsq)^2 +
             4*(c[1,3]+c[5,5])^2*stsq*ctsq)
    rv = (c[1,1] + c[5,5])*stsq + (c[3,3]+c[5,5])*ctsq - D
    return sqrt(rv/(2ρ))
end

"""
    vphaseSH(θ, c::AbstractMatrix, ρ)

Compute SH wave phase velocity along phase angle θ (with respect to axis of symmetry) based on stiffness matrix c (in Voight notation) and density ρ.
"""
function vphaseSH(θ,c::AbstractMatrix,ρ)
    return sqrt((c[6,6]*sin(θ)^2 + c[5,5]*cos(θ)^2)/ρ)
end

"""
    moduli2thomsen(c::AbstractMatrix)

Compute Thomsen's parameters (ϵ, δ, γ) based on stiffness matrix c (in Voight notation).
"""
function moduli2thomsen(c::AbstractMatrix)
    epsilon = (c[1,1] - c[3,3])/(2*c[3,3])
    delta = ((c[1,3]+c[5,5])^2 - (c[3,3]-c[5,5])^2)/(2*c[3,3]*(c[3,3]-c[5,5]))
    gamma = (c[6,6]-c[5,5])/(2*c[5,5])
    return epsilon, delta, gamma
end

"""
    velocities0(c::AbstractMatrix, ρ)

Compute phase velocities Vp₀ and Vs₀ (P and S wave speed along symmetry axis) based on stiffness matrix c (in Voight notation) and density ρ.
"""
function velocities0(c::AbstractMatrix,ρ)
    return (sqrt(c[3,3]/ρ), sqrt(c[5,5]/ρ))
end

"""
    thomsen2moduli(VP0, VS0, ρ, ϵ, δ, γ)

Compute 6×6 stiffness matrix from Vp₀, Vs₀ (P and S wave speed along symmetry axis), density ρ, and Thomsen parameters ϵ, δ and γ.
"""
function thomsen2moduli(VP0,VS0,ρ,ϵ,δ,γ)
    c33 = ρ*VP0^2
    c55 = ρ*VS0^2
    c11 = c33*(2*ϵ+1)
    c66 = c55*(2*γ+1)
    c13 = sqrt(2*c33*δ*(c33-c55) + (c33-c55)^2) - c55

    return stiffmatrixTI(c11,c13,c33,c55,c66)
end

"""
    vpvsYoungsPoisson(vp,vs,ρ)

Compute (E, ν), Youngs modulus and Poisson's ratio, from P and S wave speeds, assuming isotropy.
"""
function vpvsYoungsPoisson(vp,vs,ρ)

    nu = (0.5*(vp/vs)^2 - 1)/((vp/vs)^2 - 1)
    E = ρ*vs^2*2*(1+nu)

    return E, nu
end

"""
    stiffmatrixTI(E, nu)

Compute 6×6 stiffness matrix in an isotropic medium characterised by Young's modulus E and Poisson's ratio nu.
"""
function stiffmatrixTI(E,nu)
    c0 = E/(1+nu)/(1-2nu)
    c = zeros(typeof(E),6,6)
    c[1,1] = 1-nu
    c[2,2] = 1-nu
    c[3,3] = 1-nu
    c[1,2] = nu
    c[1,3] = nu
    c[2,3] = nu
    c[4,4] = (1-2nu)/2
    c[5,5] = (1-2nu)/2
    c[6,6] = (1-2nu)/2
    return Symmetric(c).*c0
end

"""
    stiffmatrixTI(c11, c13, c33, c55, c66)

Compute 6×6 stiffness matrix in a VTI medium characterised by components (in Voight notation) c11, c13, c33, c55 and c66.
"""
function stiffmatrixTI(c11,c13,c33,c55,c66)
    c = zeros(typeof(c11),6,6)
    c[1,1] = c11
    c[1,2] = c11-2*c66
    c[1,3] = c13
    c[2,2] = c11
    c[2,3] = c13
    c[3,3] = c33
    c[4,4] = c55
    c[5,5] = c55
    c[6,6] = c66
    return Symmetric(c)
end

"""
    dvPdthetav(θ, VP0VS0, ϵ, δ)

Compute (1/v)∂v/∂θ, derivative of P wave phase velocity  with respect to phase angle θ, normalised by Vp, based on ratio Vp₀/Vs₀ (P and S wave speed along symmetry axis) and Thomsen parameters ϵ and δ.
"""
function dvPdthetav(θ,VP0VS0,ϵ,δ)

    f = 1 - VP0VS0^(-2)

    a = sqrt((δ+f-δ*cos(4θ)-4*ϵ*cos(2θ)*sin(θ)^2+4*ϵ^2*sin(θ)^4)/f)

    vpv = (δ*sin(4θ) + ϵ*sin(2θ)*(1+ϵ-(2+ϵ)*cos(2θ) + a))/(a*(2+2*ϵ*sin(θ)^2 + f*(-1+a)))

    return vpv
end

"""
    dvPdthetav(θ, VP0VS0, ϵ, δ)

Compute (1/v)∂v/∂θ, derivative of SV wave phase velocity  with respect to phase angle θ, normalised by Vsv, based on ratio Vp₀/Vs₀ (P and S wave speed along symmetry axis) and Thomsen parameters ϵ and δ.
"""
function dvSVdthetav(θ,VP0VS0,ϵ,δ)

    f = 1 - VP0VS0^(-2)

    a = sqrt((δ+f-δ*cos(4θ)-4*ϵ*cos(2θ)*sin(θ)^2+4*ϵ^2*sin(θ)^4)/f)

    vpv = (δ*sin(4θ) - ϵ*sin(2θ)*(-1-ϵ+(2+ϵ)*cos(2θ) + a))/(a*(-2-2*ϵ*sin(θ)^2 + f*(1+a)))

    return vpv
end

"""
    vgroup(v,dvv)

Compute group wave velocity based on phase velocity v and normalised derivative of phase velocity with respect to phase angle dvv.
"""
function vgroup(v,dvv)
    return v*sqrt(1 + dvv^2)
end

"""
    groupangle(θ, dvdθv)

Compute group angle from phase angle θ and normalised derivative of phase velocity with respect to phase angle.
"""
function groupangle(θ,dvdθv)
    return atan((sin(θ) + dvdθv*cos(θ))/(cos(θ) - dvdθv*sin(θ)))
end

"""
    vgroupSH(ψ, VS0, γ)

Compute SH group velocity along group angle ψ, based on Vs₀ (S wave speed along symmetry axis) and Thomsen parameter γ.
"""
function vgroupSH(ψ,VS0,γ)
    return VS0*sqrt((1+2*γ)/(1+2*γ*cos(ψ)^2))
end

"""
    groupangleSH(θ,γ)

Compute SH group angle based on phase angle θ and Thomsen parameter γ.
"""
function groupangleSH(θ,γ)
    return atan(tan(θ)*(1+2*γ))
end

"""
    phaseangleP(ψ, VP0VS0, ϵ, δ)

Determine (using root finding method) the P wave phase angle based on group angle ψ, ratio Vp₀/Vs₀ (P and S wave speed along symmetry axis), and Thomsen parameters ϵ and δ.
"""
function phaseangleP(ψ,VP0VS0,ϵ,δ)
    fzero(θ->(ψ - groupangle(θ,dvPdthetav(θ,VP0VS0,ϵ,δ))),ψ)
end

"""
    phaseangleSV(ψ, VP0VS0, ϵ, δ)

Determine (using root finding method) the SV wave phase angle based on group angle ψ, ratio Vp₀/Vs₀ (P and S wave speed along symmetry axis), and Thomsen parameters ϵ and δ.
"""
function phaseangleSV(ψ,VP0VS0,ϵ,δ)
    fzero(θ->(ψ - groupangle(θ,dvSVdthetav(θ,VP0VS0,ϵ,δ))),ψ)
end

"""
    phaseangleSH(ψ, γ)

Compute SH wave phase angle based on group angle ψ and Thomsen parameter γ.
"""
function phaseangleSH(ψ,γ)
    return atan(tan(ψ)/(1+2*γ))
end


text:
module TestDuplicates
	using Base.Test
	using DataArrays
	using DataFrames

	df = DataFrame({"a" => [1, 2, 3, 3, 4]})
	@assert isequal(duplicated(df), [false, false, false, true, false])
	drop_duplicates!(df)
	@assert isequal(df, DataFrame({"a" => [1, 2, 3, 4]}))
end


text:
using TropicalYao
using Test

@testset "reversible" begin
    include("reversible/reversible.jl")
end


text:
using ClimaAtmos
using ClimaCore.Geometry

using ClimaCore: DataLayouts, Operators, Geometry
using ClimaAtmos.Interface: TimeStepper, Simulation
using ClimaAtmos.Interface: PeriodicRectangle, BarotropicFluidModel
using ClimaCore: Fields, Domains, Topologies, Meshes, Spaces

using IntervalSets
using UnPack
using OrdinaryDiffEq: SSPRK33
using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
using Logging: global_logger
using TerminalLoggers: TerminalLogger

using RecursiveArrayTools
using LinearAlgebra
using OrdinaryDiffEq
using Random

global_logger(TerminalLogger())

const CI = !isnothing(get(ENV, "CI", nothing))

# general parameters
const FT = Float64

# define model equations:
include("atmos_rhs.jl")

parameters = (
    # atmos physics (overwriting CliMAParameters defaults?)
    MSLP = 1e5, # mean sea level pressure [Pa]
    grav = 9.8, # gravitational acceleration [m /s^2]
    R_d = 287.058, # R dry (gas constant / mol mass dry air)  [J / K / kg]
    C_p = 287.058 * 1.4 / (1.4 - 1), # heat capacity at constant pressure [J / K / kg]
    C_v = 287.058 / (1.4 - 1), # heat capacity at constant volume [J / K / kg]
    R_m = 87.058, # moist R, assumed to be dry [J / K / kg]
    f = 7.29e-5, # Coriolis parameters [1/s]
    ν = 0.1, #0.01 # viscosity, diffusivity
    Ch = 0.0015, # bulk transfer coefficient for sensible heat
    Cd = 0.01 / (2e2 / 30.0), #drag coeff
    ug = 1.0,
    vg = 0.0,
    d = sqrt(2.0 * 0.01 / 5e-5), #?

    # radiation parameters for DryBulkFormulaWithRadiation() SurfaceFluxType
    τ = 0.9,     # atmospheric transmissivity
    α = 0.5,     # surface albedo
    σ = 5.67e-8, # Steffan Boltzmann constant [kg / s^3 / K^4]
    g_a = 0.06,    # aerodynamic conductance for heat transfer [kg / m^2 / s]
    ϵ = 0.98,    # broadband emissivity / absorptivity of the surface
    F_a = 0.0,     # downward LW flux from the atmosphere [W / m^2]
    F_sol = 1361,   # incoming solar TOA radiation [W / m^2]
    τ_d = 10,     # idealized daily cycle period [s]

    # surface fluxes
    λ = FT(0.01),#FT(1e-5)    # coupling transfer coefficient for LinearRelaxation() SurfaceFluxType
)

function atmos_simulation(
    land_simulation;
    Lz,
    Nz,
    minimum_reference_temperature = 230.0, # K
    start_time = 0.0,
    stop_time = 1.0,
    Δt_min = 0.02,
)

    # Get surface temperature from the land simulation
    land_surface_temperature = land_sim.p[5]

    ########
    # Set up atmos domain
    ########

    domain_atm = Domains.IntervalDomain(0, Lz, x3boundary = (:bottom, :top)) # struct
    mesh_atm = Meshes.IntervalMesh(domain_atm, nelems = Nz) # struct, allocates face boundaries to 5,6
    center_space_atm = Spaces.CenterFiniteDifferenceSpace(mesh_atm) # collection of the above, discretises space into FD and provides coords
    face_space_atm = Spaces.FaceFiniteDifferenceSpace(center_space_atm)

    """ Initialize fields located at cell centers in the vertical. """
    function init_centers(zc, parameters)
        UnPack.@unpack grav, C_p, MSLP, R_d = parameters

        # temperature
        Γ = grav / C_p
        T = max(land_surface_temperature - Γ * zc, minimum_reference_temperature)

        # pressure
        p = MSLP * (T / land_surface_temperature)^(grav / (R_d * Γ))

        if T == minimum_reference_temperature
            z_top = (land_surface_temperature - minimum_reference_temperature) / Γ
            H_min = R_d * minimum_reference_temperature / grav
            p *= exp(-(zc - z_top) / H_min)
        end

        # potential temperature
        θ = land_surface_temperature

        # density
        ρ = p / (R_d * θ * (p / MSLP)^(R_d / C_p))

        # velocties
        u = 1.0
        v = 0.0

        return (ρ = ρ, u = u, v = v, ρθ = ρ * θ)
    end

    """ Initialize fields located at cell interfaces in the vertical. """
    function init_faces(zf, parameters)
        return (; w = 0.0 .* zf)
    end

    # Initialize the atmospheric states Yc and Yf
    z_centers = Fields.coordinate_field(center_space_atm)
    z_faces = Fields.coordinate_field(face_space_atm)
    Yc = init_centers.(z_centers, Ref(parameters))
    Yf = init_faces.(z_faces, Ref(parameters))
    T_atm_0 = (Yc = Yc, Yf = Yf)

    # Put all prognostic variable arrays into a vector and ensure that solve can partition them
    Y_atm = ArrayPartition((T_atm_0.Yc, T_atm_0.Yf, zeros(3)))
    prob_atm = ODEProblem(∑tendencies_atm!, Y_atm, (start_time, stop_time), (parameters, [land_surface_temperature]))
    simulation = init(prob_atm, SSPRK33(), dt = Δt_min, saveat = 1 * Δt_min)

    return simulation
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
function read_image(file, imgpath)
    target = joinpath(imgpath, file)
    img = load(target)
    return img
end

function process_image(
    img, newsize; rgbmean=nothing, crop=true, randomcrop=false)
    scaled = ntuple(i->div(size(img,i)*newsize[i],minimum(size(img))),2)
    a1 = Images.imresize(img, scaled)

    # randomcrop vs. centercrop
    if randomcrop
        offsets = ntuple(i->rand(1:scaled[i]-minimum(scaled)+1),2)
    else
        offsets = ntuple(i->div(size(a1,i)-newsize[i],2)+1,2)
    end

    if crop
        a1 = a1[offsets[1]:offsets[1]+newsize[1]-1,
                offsets[2]:offsets[2]+newsize[2]-1]
    else
        a1 = Images.imresize(a1, newsize)
    end

    b1 = permutedims(collect(channelview(a1)), (3,2,1))
    colordim = size(b1, 3)
    if colordim != 3
        c1 = convert(Array{Float32}, b1)
        c1 = cat(3, cat(3, c1, c1), c1)
    else
        c1 = convert(Array{Float32}, b1) # type conversion
    end
    d1 = reshape(c1[:,:,1:3], (newsize[1],newsize[2],3,1)) # reshape
    e1 = (255 * d1) # 8bit image representation
    if rgbmean != nothing
        e1 = e1 .- rgbmean
    end

    return permutedims(e1, [2,1,3,4]) # transpose
end

function make_average_image(avgimg)
    output = reshape(avgimg, size(avgimg)..., 1)
    output = permutedims(output, [2,1,3,4])
    output = convert(Array{Float32}, output)
end

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
# Autogenerated wrapper script for METIS4_jll for aarch64-linux-gnu
export libmetis_a

JLLWrappers.@generate_wrapper_header("METIS4")
JLLWrappers.@declare_file_product(libmetis_a)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        libmetis_a,
        "lib/libmetis.a",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
########## Vanilla GBP messages: Factor to variable ##########
function messageFactorVariable(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateFactor
        for j in gbp.graph.rowptr[i]
            Mrow = gbp.graph.meanIndirect[j]; Vrow = gbp.graph.varianceIndirect[j]; q = gbp.graph.toVariable.nzval[j]
            for k in gbp.graph.rowptr[i]
                if j != k
                    Mrow -= gbp.graph.coefficient[k] * gbp.inference.meanVariableFactor[k]
                    Vrow += gbp.graph.coefficient[k]^2 * gbp.inference.varianceVariableFactor[k]
                end
            end
            gbp.inference.meanFactorVariable[q] = Mrow / gbp.graph.coefficient[j]
            gbp.inference.varianceFactorVariable[q] = Vrow / (gbp.graph.coefficient[j]^2)
        end
    end
end

########## Vanilla GBP means: Factor to variable ##########
function meanFactorVariable(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateFactor
        for j in gbp.graph.rowptr[i]
            Mrow = gbp.graph.meanIndirect[j]
            for k in gbp.graph.rowptr[i]
                if j != k
                    Mrow -= gbp.graph.coefficient[k] * gbp.inference.meanVariableFactor[k]
                end
            end
            gbp.inference.meanFactorVariable[gbp.graph.toVariable.nzval[j]] = Mrow / gbp.graph.coefficient[j]
        end
    end
end

########## Vanilla GBP variances: Factor to variable ##########
function varianceFactorVariable(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateFactor
        for j in gbp.graph.rowptr[i]
            Vrow = gbp.graph.varianceIndirect[j]
            for k in gbp.graph.rowptr[i]
                if j != k
                    Vrow += gbp.graph.coefficient[k]^2 * gbp.inference.varianceVariableFactor[k]
                end
            end
            gbp.inference.varianceFactorVariable[gbp.graph.toVariable.nzval[j]] = Vrow / (gbp.graph.coefficient[j]^2)
        end
    end
end

########## Vanilla GBP damp messages: Factor to variable ##########
function messageDampFactorVariable(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateFactor
        for j in gbp.graph.rowptr[i]
            Mrow = gbp.graph.meanIndirect[j]; Vrow = gbp.graph.varianceIndirect[j]; q = gbp.graph.toVariable.nzval[j]
            for k in gbp.graph.rowptr[i]
                if j != k
                    Mrow -= gbp.graph.coefficient[k] * gbp.inference.meanVariableFactor[k]
                    Vrow += gbp.graph.coefficient[k]^2 * gbp.inference.varianceVariableFactor[k]
                end
            end
            gbp.inference.meanFactorVariable[q] = gbp.graph.alphaNew[j] * Mrow / gbp.graph.coefficient[j] + gbp.graph.alphaOld[j] * gbp.inference.meanFactorVariable[q]
            gbp.inference.varianceFactorVariable[q] = Vrow / (gbp.graph.coefficient[j]^2)
        end
    end
end

########## Vanilla GBP damp means: Factor to variable ##########
function meanDampFactorVariable(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateFactor
        for j in gbp.graph.rowptr[i]
            Mrow = gbp.graph.meanIndirect[j]; q = gbp.graph.toVariable.nzval[j]
            for k in gbp.graph.rowptr[i]
                if j != k
                    Mrow -= (gbp.graph.coefficient[k] * gbp.inference.meanVariableFactor[k])
                end
            end
            gbp.inference.meanFactorVariable[q] = gbp.graph.alphaNew[j] * Mrow / gbp.graph.coefficient[j] + gbp.graph.alphaOld[j] * gbp.inference.meanFactorVariable[q]
        end
    end
end

############# Vanilla GBP messages: Variable to factor ##########
function messageVariableFactor(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateVariable
        for j in gbp.graph.colptr[i]
            Mcol = gbp.graph.meanDirect[i]; Wcol = gbp.graph.weightDirect[i]; q = gbp.graph.toFactor.nzval[j]
            for k in gbp.graph.colptr[i]
                if j != k
                    Mcol += gbp.inference.meanFactorVariable[k] / gbp.inference.varianceFactorVariable[k]
                    Wcol += 1 / gbp.inference.varianceFactorVariable[k]
                end
            end
            gbp.inference.varianceVariableFactor[q] = 1 / Wcol
            gbp.inference.meanVariableFactor[q] = Mcol / Wcol
        end
    end
end

############# Vanilla GBP means: Variable to factor ##########
function meanVariableFactor(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateVariable
        for j in gbp.graph.colptr[i]
            Mcol = gbp.graph.meanDirect[i]; q = gbp.graph.toFactor.nzval[j]
            for k in gbp.graph.colptr[i]
                if j != k
                    Mcol += gbp.inference.meanFactorVariable[k] / gbp.inference.varianceFactorVariable[k]
                end
            end
            gbp.inference.meanVariableFactor[q] = Mcol * gbp.inference.varianceVariableFactor[q]
        end
    end
end

############# Vanilla GBP variances: Variable to factor ##########
function varianceVariableFactor(gbp::ContinuousModel)
    @inbounds Threads.@threads for i in gbp.graph.iterateVariable
        for j in gbp.graph.colptr[i]
            Wcol = gbp.graph.weightDirect[i]
            for k in gbp.graph.colptr[i]
                if j != k
                    Wcol += 1 / gbp.inference.varianceFactorVariable[k]
                end
            end
            gbp.inference.varianceVariableFactor[gbp.graph.toFactor.nzval[j]] = 1 / Wcol
        end
    end
end