text:
using MAT
D = matread("TripAdvisor_NAM_Europe_H.mat")
H = D["H"]
classes = D["classes"]

alpha, B1, B2, H1, H2, R1, R2, N = Hypergraph_to_Scores(H,classes,25)
class1 = "N. America"
class2 = "Europe"

## Plot
using Plots
lw = 2.5
ms = 9
gfs = 15
tfs = 15
titlesize = 14
s1 = 350
s2 = 300
for k = 3:3:9
color1 = RGB(27/255,158/255,119/255)
color2 = RGB(217/255,95/255,2/255)
title = "$k-hotel reviews"
title = ""
leg = false
xlab = "Affinity type t"
ylab = "Affinity"
p = Plots.plot(fontfamily = "helvetica",linewidth = 2,legend = false,grid = false, size = (s1,s2),color = :gray)
title!(p,title,titlefontsize = titlesize)
xaxis!(p,xlab, xticks = k/3:k/3:k,tickfontsize=tfs,guidefontsize = gfs)
yaxis!(p,ylab,tickfontsize=tfs,guidefontsize = gfs)

plot!(p,1:k,H2[k,1:k],color = color2, linewidth = lw, markersize = ms,
    markershape = :circle, markerstrokecolor = color2, label = class2)
plot!(p,1:k,H1[k,1:k],color = color1, linewidth = lw, markersize = ms,
    markershape = :circle,markerstrokecolor = color1, label = class1)

savefig("Figures/Affinities_$(k)_TA.pdf")
end

## Next
for k = 3:3:9
s1 = 400
s2 = 300
title = ""
leg = false
xlab = "Affinity type t"
ylab = "Affinity / Baseline"
p = Plots.plot([1; k],[1; 1],yscale = :log10,fontfamily = "helvetica",linewidth = 2,legend = false,grid = false, size = (s1,s2),color = :gray)
title!(p,title,titlefontsize = titlesize)
xaxis!(p,xlab, xticks = k/3:k/3:k,tickfontsize=tfs,guidefontsize = gfs)
yaxis!(p,ylab,tickfontsize=tfs,guidefontsize = gfs)
plot!(p,1:k,R1[k,1:k],color = color1, linewidth = lw, markersize = ms,
    markershape = :circle,markerstrokecolor = color1, label = class1)
plot!(p,1:k,R2[k,1:k],color = color2, linewidth = lw, markersize = ms,
    markershape = :circle, markerstrokecolor = color2, label = class2)

savefig("Figures/Ratios_$(k)_TA.pdf")
end

## Group Homophily Index
r = 15
S1,S2  = GHI(H,classes,r)
print("\n \$k\$ ")
start = 2
for k = start:13
    print("& $(k) ")
end
print("\\\\\n \\midrule \nN. America GHI ")
for k = start:13
    print("& $(round(Int64,S1[k])) ")
end
print("\\\\ \n Europe GHI")
for k = start:13
    print("& $(round(Int64,S2[k])) ")
end
print("\\\\ \n")

## use this table
# {\Large
# \begin{tabular}{lllllllllllllllll}
# \toprule
# group size $k$ & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 \\
#  \midrule
# N. America GHI & 1 & 1 & 1 & 1 & 2 & 2 & 2 & 3 & 3 & 3 & 4 \\
#  Europe GHI& 1 & 1 & 1 & 1 & 1 & 1 & 2 & 3 & 3 & 3 & 3 \\
# \bottomrule
# \end{tabular}
# }

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily


text:
snouter_tree = read_newick("../data/snouters.nwk")[1]
root!(snouter_tree)

dir = SpeciesDirectory(tiplabels(snouter_tree))
snouter_tree.dir = dir
setspecies!(snouter_tree)

@testset "Simulate Brownian Motion" begin
    x₀ = [1.0, 12.0, 3.4, 2.0, 12.0]
    σ² = [1.0, 2.0, 10.0, 100.0, 0.0]

    dm = simulate_bm(snouter_tree)
    @test dm isa SpeciesDataMatrix{Float64}
    @test size(dm) == (n_species(snouter_tree), 1)

    dm = simulate_bm(snouter_tree, 12)
    @test dm isa SpeciesDataMatrix{Float64}
    @test size(dm) == (n_species(snouter_tree), 12)

    dm = simulate_bm(snouter_tree, σ²=3.0)
    @test dm isa SpeciesDataMatrix{Float64}
    @test size(dm) == (n_species(snouter_tree), 1)

    dm = simulate_bm(snouter_tree, σ²=3.0, x₀=x₀)
    @test dm isa SpeciesDataMatrix{Float64}
    @test size(dm) == (n_species(snouter_tree), 5)

    dm = simulate_bm(snouter_tree, σ²=σ²)
    @test dm isa SpeciesDataMatrix{Float64}
    @test size(dm) == (n_species(snouter_tree), 5)
    @test dm[:,5] == zeros(10)

    dm = simulate_bm(snouter_tree, x₀=x₀)
    @test dm isa SpeciesDataMatrix{Float64}
    @test size(dm) == (n_species(snouter_tree), 5)

    dm = simulate_bm(snouter_tree, x₀=x₀, σ²=σ²)
    @test dm isa SpeciesDataMatrix{Float64}
    @test size(dm) == (n_species(snouter_tree), 5)
end

ext: jl


text:
# Path Primitives

# From form.jl:
# Path
# ----

# An implementation of the SVG path mini-language.

abstract type PathOp end

struct MoveAbsPathOp <: PathOp
    to::Vec
end

function assert_pathop_tokens_len(op_type, tokens, i, needed)
    provided = length(tokens) - i + 1
    provided < needed &&
            error("In path $(op_type) requires $(needed) argumens but only $(provided) provided.")
end

function parsepathop(::Type{MoveAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(MoveAbsPathOp, tokens, i, 2)
    op = MoveAbsPathOp((x_measure(tokens[i]), y_measure(tokens[i + 1])))
    return (op, i + 2)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::MoveAbsPathOp) =
        MoveAbsPathOp(resolve(box, units, t, p.to))

struct MoveRelPathOp <: PathOp
    to::Vec
end

function parsepathop(::Type{MoveRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(MoveRelPathOp, tokens, i, 2)
    op = MoveRelPathOp((tokens[i], tokens[i + 1]))
    return (op, i + 2)
end

function resolve_offset(box::AbsoluteBox, units::UnitBox, t::Transform, p::Vec)
    absp = resolve(box, units, t, p)
    zer0 = resolve(box, units, t, (0w, 0h))
    return (absp[1] - zer0[1], absp[2] - zer0[2])
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::MoveRelPathOp) =
        MoveRelPathOp(resolve_offset(box, units, t, p.to))

struct ClosePathOp <: PathOp
end

parsepathop(::Type{ClosePathOp}, tokens::AbstractArray, i) = (ClosePathOp(), i)

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::ClosePathOp) = p

struct LineAbsPathOp <: PathOp
    to::Vec
end

function parsepathop(::Type{LineAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(LineAbsPathOp, tokens, i, 2)
    op = LineAbsPathOp((x_measure(tokens[i]), y_measure(tokens[i + 1])))
    return (op, i + 2)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::LineAbsPathOp) =
        LineAbsPathOp(resolve(box, units, t, p.to))

struct LineRelPathOp <: PathOp
    to::Vec
end

function parsepathop(::Type{LineRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(LineRelPathOp, tokens, i, 2)
    op = LineRelPathOp((x_measure(tokens[i]), y_measure(tokens[i + 1])))
    return (op, i + 2)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::LineRelPathOp) =
        LineRelPathOp(resolve(box, units, t, p.to))

struct HorLineAbsPathOp <: PathOp
    x::Measure
end

function parsepathop(::Type{HorLineAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(HorLineAbsPathOp, tokens, i, 1)
    op = HorLineAbsPathOp(x_measure(tokens[i]))
    return (op, i + 1)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::HorLineAbsPathOp) =
        HorLineAbsPathOp(resolve(box, units, t, (p.x, 0mm))[1])

struct HorLineRelPathOp <: PathOp
    Δx::Measure
end

function parsepathop(::Type{HorLineRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(HorLineRelPathOp, tokens, i, 1)
    op = HorLineRelPathOp(x_measure(tokens[i]))
    return (op, i + 1)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::HorLineRelPathOp) =
        HorLineRelPathOp(resolve(box, units, t, p.Δx))

struct VertLineAbsPathOp <: PathOp
    y::Measure
end

function parsepathop(::Type{VertLineAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(VertLineAbsPathOp, tokens, i, 1)
    op = VertLineAbsPathOp(y_measure(tokens[i]))
    return (op, i + 1)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::VertLineAbsPathOp) =
        VertLineAbsPathOp(resolve(box, units, t, (0mm, p.y))[2])

struct VertLineRelPathOp <: PathOp
    Δy::Measure
end

function parsepathop(::Type{VertLineRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(VertLineRelPathOp, tokens, i, 1)
    op = VertLineRelPathOp(y_measure(tokens[i]))
    return (op, i + 1)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::VertLineRelPathOp) =
        VertLineAbsPathOp(resolve(box, units, t, (0mmm, p.Δy))[2])

struct CubicCurveAbsPathOp <: PathOp
    ctrl1::Vec
    ctrl2::Vec
    to::Vec
end

function parsepathop(::Type{CubicCurveAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(CubicCurveAbsPathOp, tokens, i, 6)
    op = CubicCurveAbsPathOp((tokens[i],     tokens[i + 1]),
                             (tokens[i + 2], tokens[i + 3]),
                             (tokens[i + 4], tokens[i + 5]))
    return (op, i + 6)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::CubicCurveAbsPathOp) =
        CubicCurveAbsPathOp(
            resolve(box, units, t, p.ctrl1),
            resolve(box, units, t, p.ctrl2),
            resolve(box, units, t, p.to))

struct CubicCurveRelPathOp <: PathOp
    ctrl1::Vec
    ctrl2::Vec
    to::Vec
end

function parsepathop(::Type{CubicCurveRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(CubicCurveRelPathOp, tokens, i, 6)
    op = CubicCurveRelPathOp((tokens[i],     tokens[i + 1]),
                             (tokens[i + 2], tokens[i + 3]),
                             (tokens[i + 4], tokens[i + 5]))
    return (op, i + 6)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::CubicCurveRelPathOp) =
        CubicCurveRelPathOp(
            resolve(box, units, t, p.ctrl1),
            resolve(box, units, t, p.ctrl2),
            resolve(box, units, t, p.to))

struct CubicCurveShortAbsPathOp <: PathOp
    ctrl2::Vec
    to::Vec
end

function parsepathop(::Type{CubicCurveShortAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(CubicCurveShortAbsPathOp, tokens, i, 4)
    op = CubicCurveShortAbsPathOp((x_measure(tokens[i]),     y_measure(tokens[i + 1])),
                                  (x_measure(tokens[i + 2]), y_measure(tokens[i + 3])))
    return (op, i + 4)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::CubicCurveShortAbsPathOp) =
        CubicCurveShortAbsPathOp(
            resolve_offset(box, units, t, p.ctrl2),
            resolve_offset(box, units, t, p.to))

struct CubicCurveShortRelPathOp <: PathOp
    ctrl2::Vec
    to::Vec
end

function parsepathop(::Type{CubicCurveShortRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(CubicCurveShortRelPathOp, tokens, i, 4)
    op = CubicCurveShortRelPathOp((x_measure(tokens[i]),     y_measure(tokens[i + 1])),
                                  (x_measure(tokens[i + 2]), y_measure(tokens[i + 3])))
    return (op, i + 4)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::CubicCurveShortRelPathOp) =
        CubicCurveShortRelPathOp(
            resolve(box, units, t, p.ctrl2),
            resolve(box, units, t, p.to))

struct QuadCurveAbsPathOp <: PathOp
    ctrl1::Vec
    to::Vec
end

function parsepathop(::Type{QuadCurveAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(QuadCurveAbsPathOp, tokens, i, 4)
    op = QuadCurveAbsPathOp((tokens[i],     tokens[i + 1]),
                            (tokens[i + 2], tokens[i + 3]))
    return (op, i + 4)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::QuadCurveAbsPathOp) =
        QuadCurveAbsPathOp(
            resolve(box, units, t, p.ctrl1),
            resolve(box, units, t, p.to))

struct QuadCurveRelPathOp <: PathOp
    ctrl1::Vec
    to::Vec
end

function parsepathop(::Type{QuadCurveRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(QuadCurveRelPathOp, tokens, i, 4)
    op = QuadCurveRelPathOp((tokens[i],     tokens[i + 1]),
                            (tokens[i + 2], tokens[i + 3]))
    return (op, i + 4)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::QuadCurveRelPathOp) =
        QuadCurveRelPathOp(
            (resolve(box, units, t, p.ctrl1[1]),
             resolve(box, units, t, p.ctrl1[2])),
            (resolve(box, units, t, p.to[1]),
             resolve(box, units, t, p.to[2])))

struct QuadCurveShortAbsPathOp <: PathOp
    to::Vec
end

function parsepathop(::Type{QuadCurveShortAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(QuadCurveShortAbsPathOp, tokens, i, 2)
    op = QuadCurveShortAbsPathOp((tokens[i], tokens[i + 1]))
    return (op, i + 2)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::QuadCurveShortAbsPathOp) =
        QuadCurveShortAbsPathOp(resolve(box, units, t, p.to))

struct QuadCurveShortRelPathOp <: PathOp
    to::Vec
end

function parsepathop(::Type{QuadCurveShortRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(QuadCurveShortRelPathOp, tokens, i, 2)
    op = QuadCurveShortRelPathOp((tokens[i], tokens[i + 1]))
    return (op, i + 2)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::QuadCurveShortRelPathOp) =
        QuadCurveShortRelPathOp(
            (resolve(box, units, t, p.to[1]),
             resolve(box, units, t, p.to[2])))

struct ArcAbsPathOp <: PathOp
    rx::Measure
    ry::Measure
    rotation::Float64
    largearc::Bool
    sweep::Bool
    to::Vec
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::ArcAbsPathOp) =
        ArcAbsPathOp(
            resolve(box, units, t, p.rx),
            resolve(box, units, t, p.ry),
            p.rotation,
            p.largearc,
            p.sweep,
            resolve(box, units, t, p.to))

struct ArcRelPathOp <: PathOp
    rx::Measure
    ry::Measure
    rotation::Float64
    largearc::Bool
    sweep::Bool
    to::Vec
end

function parsepathop(::Type{T}, tokens::AbstractArray, i) where T <: Union{ArcAbsPathOp, ArcRelPathOp}
    assert_pathop_tokens_len(T, tokens, i, 7)

    if isa(tokens[i + 3], Bool)
        largearc = tokens[i + 3]
    elseif tokens[i + 3] == 0
        largearc = false
    elseif tokens[i + 3] == 1
        largearc = true
    else
        error("largearc argument to the arc path operation must be boolean")
    end

    if isa(tokens[i + 4], Bool)
        sweep = tokens[i + 4]
    elseif tokens[i + 4] == 0
        sweep = false
    elseif tokens[i + 4] == 1
        sweep = true
    else
        error("sweep argument to the arc path operation must be boolean")
    end

    isa(tokens[i + 2], Number) || error("path arc operation requires a numerical rotation")

    op = T(x_measure(tokens[i]),
           y_measure(tokens[i + 1]),
           convert(Float64, tokens[i + 2]),
           largearc, sweep,
           (x_measure(tokens[i + 5]), y_measure(tokens[i + 6])))

    return (op, i + 7)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::ArcRelPathOp) =
        ArcRelPathOp(
            resolve(box, units, t, p.rx),
            resolve(box, units, t, p.ry),
            p.rotation,
            p.largearc,
            p.sweep,
            (resolve(box, units, t, p.to[1]),
             resolve(box, units, t, p.to[2])))

const path_ops = Dict(
     :M => MoveAbsPathOp,
     :m => MoveRelPathOp,
     :Z => ClosePathOp,
     :z => ClosePathOp,
     :L => LineAbsPathOp,
     :l => LineRelPathOp,
     :H => HorLineAbsPathOp,
     :h => HorLineRelPathOp,
     :V => VertLineAbsPathOp,
     :v => VertLineRelPathOp,
     :C => CubicCurveAbsPathOp,
     :c => CubicCurveRelPathOp,
     :S => CubicCurveShortAbsPathOp,
     :s => CubicCurveShortRelPathOp,
     :Q => QuadCurveAbsPathOp,
     :q => QuadCurveRelPathOp,
     :T => QuadCurveShortAbsPathOp,
     :t => QuadCurveShortRelPathOp,
     :A => ArcAbsPathOp,
     :a => ArcRelPathOp
)

# A path is an array of symbols, numbers, and measures following SVGs path
# mini-language.
function parsepath(tokens::AbstractArray)
    ops = PathOp[]
    last_op_type = nothing
    i = 1
    while i <= length(tokens)
        tok = tokens[i]
        strt = i
        if isa(tok, Symbol)
            if !haskey(path_ops, tok)
                error("$(tok) is not a valid path operation")
            else
                op_type = path_ops[tok]
                i += 1
                op, i = parsepathop(op_type, tokens, i)
                push!(ops, op)
                last_op_type = op_type
            end
        else
            op, i = parsepathop(last_op_type, tokens, i)
            push!(ops, op)
        end
    end

    return ops
end

struct PathPrimitive <: FormPrimitive
    ops::Vector{PathOp}
end

const Path = Form{PathPrimitive}

path(tokens::AbstractArray, tag=empty_tag) = Path([PathPrimitive(parsepath(tokens))], tag)

path(tokens::AbstractArray{T}, tag=empty_tag) where T <: AbstractArray =
        Path([PathPrimitive(parsepath(ts)) for ts in tokens], tag)

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::PathPrimitive) =
        PathPrimitive([resolve(box, units, t, op) for op in p.ops])

# TODO: boundingbox

# From cairo_backends.jl:

function draw(img::Image, prim::PathPrimitive)
    for op in prim.ops
        draw_path_op(img, op)
    end
    fillstroke(img)
end

draw_path_op(img::Image, op::MoveAbsPathOp) = move_to(img, op.to)
draw_path_op(img::Image, op::MoveRelPathOp) = rel_move_to(img, op.to)
draw_path_op(img::Image, op::ClosePathOp)   = close_path(img)
draw_path_op(img::Image, op::LineAbsPathOp) = line_to(img, op.to)
draw_path_op(img::Image, op::LineRelPathOp) = rel_line_to(img, op.to)

function draw_path_op(img::Image, op::HorLineAbsPathOp)
    pos = current_point(img)
    line_to(img, (op.x, pos.y))
end

draw_path_op(img::Image, op::HorLineRelPathOp) = rel_line_to(img, (op.Δx, 0.0mm))

function draw_path_op(img::Image, op::VertLineAbsPathOp)
    pos = current_point(img)
    line_to(img, (pos.x, op.y))
end

draw_path_op(img::Image, op::VertLineRelPathOp) = rel_line_to(img, (0.0mm, op.Δy))

function draw_path_op(img::Image, op::CubicCurveAbsPathOp)
    curve_to(img, op.ctrl1, op.ctrl2, op.to)
    img.last_ctrl2_point = op.ctrl2
end

function draw_path_op(img::Image, op::CubicCurveRelPathOp)
    xy = current_point(img)
    rel_curve_to(img, op.ctrl1, op.ctrl2, op.to)
    img.last_ctrl2_point = (op.ctrl2[1] + xy[1], op.ctrl2[2] + xy[2])
end

function draw_path_op(img::Image, op::CubicCurveShortAbsPathOp)
    xy = current_point(img)
    ctrl1 = img.last_ctrl2_point
    if ctrl1 === nothing
        ctrl1 = xy
    else
        ctrl1 = (2*xy[1] - ctrl1[1], 2*xy[2] - ctrl1[2])
    end
    curve_to(img, ctrl1, op.ctrl2, op.to)
    img.last_ctrl2_point = op.ctrl2
end

function draw_path_op(img::Image, op::CubicCurveShortRelPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = op.to[1].value, op.to[2].value

    ctrl1 = img.last_ctrl2_point
    if ctrl1 === nothing
        ctrl1 = xy
    else
        ctrl1 = (Measure(abs=(2*x1 - ctrl1[1].value) - x1),
                 Measure(abs=(2*y1 - ctrl1[2].value) - y1))
    end
    cx, cy = ctrl1[1].value, ctrl1[2].value

    rel_curve_to(img, ctrl1, op.ctrl2, op.to)
    img.last_ctrl2_point =
        (Measure(abs=op.ctrl2[1].value + xy.x.abs),
         Measure(abs=op.ctrl2[2].value + xy.y.abs))
end

function draw_path_op(img::Image, op::QuadCurveAbsPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = op.to[1].value, op.to[2].value
    cx, cy = op.ctrl1[1].value, op.ctrl1[2].value
    curve_to(img,
             (Measure(abs=(x1 + 2*cx)/3),
              Measure(abs=(y1 + 2*cy)/3)),
             (Measure(abs=(x2 + 2*cx)/3),
              Measure(abs=(y2 + 2*cy)/3)),
             op.to)
    img.last_ctrl1_point = op.ctrl1
end

function draw_path_op(img::Image, op::QuadCurveRelPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = op.to[1].value, op.to[2].value
    cx, cy = op.ctrl1[1].value, op.ctrl1[2].value
    rel_curve_to(img,
                 (Measure(abs=(x1 + 2*cx)/3),
                  Measure(abs=(y1 + 2*cy)/3)),
                 (Measure(abs=(x2 + 2*cx)/3),
                  Measure(abs=(y2 + 2*cy)/3)),
             op.to)
    img.last_ctrl1_point =
        (Measure(abs=op.ctrl1[1].value + xy.x.abs),
         Measure(abs=op.ctrl1[2].value + xy.y.abs))
end

function draw_path_op(img::Image, op::QuadCurveShortAbsPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = op.to[1].value, op.to[2].value

    ctrl1 = img.last_ctrl1_point
    if img.last_ctrl1_point === nothing
        ctrl1 = xy
    else
        ctrl1 = (Measure(abs=2*x1 - ctrl1[1].value),
                 Measure(abs=2*y1 - ctrl1[2].value))
    end
    cx, cy = ctrl1[1].value, ctrl1[2].value

    curve_to(img,
             (Measure(abs=(x1 + 2*cx)/3),
              Measure(abs=(y1 + 2*cy)/3)),
             (Measure(abs=(x2 + 2*cx)/3),
              Measure(abs=(y2 + 2*cy)/3)),
             (Measure(abs=x2), Measure(abs=y2)))
    img.last_ctrl1_point = ctrl1
end

function draw_path_op(img::Image, op::QuadCurveShortRelPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = x1 + op.to[1].value, y1 + op.to[2].value

    ctrl1 = img.last_ctrl1_point
    if ctrl1 === nothing
        ctrl1 = xy
    else
        ctrl1 = (Measure(abs=(2*x1 - ctrl1[1].value) - x1),
                 Measure(abs=(2*y1 - ctrl1[2].value) - y1))
    end
    cx, cy = ctrl1[1].value, ctrl1[2].value

    rel_curve_to(img,
                 (Measure(abs=(x1 + 2*cx)/3),
                  Measure(abs=(y1 + 2*cy)/3)),
                 (Measure(abs=(x2 + 2*cx)/3),
                  Measure(abs=(y2 + 2*cy)/3)),
                 (Measure(abs=x2), Measure(abs=y2)))
    img.last_ctrl1_point =
        (Measure(abs=op.ctrl1[1].value + x1),
         Measure(abs=op.ctrl1[2].value + y1))
end

function draw_path_op(img::Image, op::ArcAbsPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = op.to[1].value, op.to[2].value
    rx, ry = op.rx.abs, op.ry.abs
    φ = deg2rad(op.rotation)
    draw_endpoint_arc(img, rx, ry, φ, op.largearc, op.sweep, x1, y1, x2, y2)
end

function draw_path_op(img::Image, op::ArcRelPathOp)
    xy = current_point(img)
    x1, y1 = xy[1].value, xy[2].value
    x2, y2 = x1 + op.to[1].value, y1 + op.to[2].value
    rx, ry = op.rx.abs, op.ry.abs
    φ = deg2rad(op.rotation)
    draw_endpoint_arc(img, rx, ry, φ, op.largearc, op.sweep, x1, y1, x2, y2)
end

# Draw an SVG style elliptical arc
function draw_endpoint_arc(img::Image, rx::Float64, ry::Float64, φ::Float64,
                           largearc::Bool, sweep::Bool,
                           x1::Float64, y1::Float64,
                           x2::Float64, y2::Float64)
    function uvangle(ux, uy, vx, vy)
        t = (ux * vx + uy * vy) / (sqrt(ux^2 + uy^2) * sqrt(vx^2 + vy^2))
        t = max(min(t, 1.0), -1.0)
        return (ux * vy - uy * vx < 0.0 ? -1 : 1.0) * acos(t)
    end

    # From: http://www.w3.org/TR/SVG/implnote.html#ArcConversionEndpointToCenter
    xm, ym = (x1 - x2)/2, (y1 - y2)/2
    x1p =  cos(φ) * xm + sin(φ) * ym
    y1p = -sin(φ) * xm + cos(φ) * ym

    u = (rx^2 * ry^2 - rx^2 * y1p^2 - ry^2 * x1p^2) / (rx^2 * y1p^2 + ry^2 * x1p^2)
    u = u >= 0.0 ? sqrt(u) : 0.0

    cxp =  u * (rx * y1p) / ry
    cyp = -u * (ry * x1p) / rx
    if sweep == largearc
        cxp = -cxp
        cyp = -cyp
    end
    cx = (x1 + x2)/2 + cos(φ) * cxp - sin(φ) * cyp
    cy = (y1 + y2)/2 + sin(φ) * cxp + cos(φ) * cyp

    θ1 = uvangle(1.0, 0.0, (x1p - cxp) / rx, (y1p - cyp) / ry)
    Δθ = uvangle((x1p - cxp) / rx, (y1p - cyp) / ry,
                 (-x1p - cxp) / rx, (-y1p - cyp) / ry) % (2.0*π)
    if Δθ > 0.0 && !sweep
        Δθ -= 2*π
    elseif Δθ < 0.0 && sweep
        Δθ += 2*π
    end

    Cairo.save(img.ctx)
    Cairo.translate(img.ctx,
                    absolute_native_units(img, cx),
                    absolute_native_units(img, cy))
    Cairo.rotate(img.ctx, φ)
    Cairo.scale(img.ctx, rx, ry)
    if sweep
        arc(img, 0.0, 0.0, 1.0, θ1, θ1 + Δθ)
    else
        arc_negative(img, 0.0, 0.0, 1.0, θ1, θ1 + Δθ)
    end
    Cairo.restore(img.ctx)
end

# From svg.jl:

function svg_print_path_op(io::IO, op::MoveAbsPathOp)
    print(io, 'M')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

function svg_print_path_op(io::IO, op::MoveRelPathOp)
    print(io, 'm')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

svg_print_path_op(io::IO, op::ClosePathOp) = print(io, 'z')

function svg_print_path_op(io::IO, op::LineAbsPathOp)
    print(io, 'L')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

function svg_print_path_op(io::IO, op::LineRelPathOp)
    print(io, 'l')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

function svg_print_path_op(io::IO, op::HorLineAbsPathOp)
    print(io, 'H')
    svg_print_float(io, op.x.value)
end

function svg_print_path_op(io::IO, op::HorLineRelPathOp)
    print(io, 'h')
    svg_print_float(io, op.Δx.value)
end

function svg_print_path_op(io::IO, op::VertLineAbsPathOp)
    print(io, 'V')
    svg_print_float(io, op.y.value)
end

function svg_print_path_op(io::IO, op::VertLineRelPathOp)
    print(io, 'v')
    svg_print_float(io, op.Δy.value)
end

function svg_print_path_op(io::IO, op::CubicCurveAbsPathOp)
    print(io, 'C')
    svg_print_float(io, op.ctrl1[1].value)
    print(io, ' ')
    svg_print_float(io, op.ctrl1[2].value)
    print(io, ' ')
    svg_print_float(io, op.ctrl2[1].value)
    print(io, ' ')
    svg_print_float(io, op.ctrl2[2].value)
    print(io, ' ')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

function svg_print_path_op(io::IO, op::CubicCurveRelPathOp)
    print(io, 'c')
    svg_print_float(io, op.ctrl1[1].value)
    print(io, ' ')
    svg_print_float(io, op.ctrl1[2].value)
    print(io, ' ')
    svg_print_float(io, op.ctrl2[1].value)
    print(io, ' ')
    svg_print_float(io, op.ctrl2[2].value)
    print(io, ' ')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

function svg_print_path_op(io::IO, op::CubicCurveShortAbsPathOp)
    print(io, 'S')
    svg_print_float(io, op.ctrl2[1].value)
    print(io, ' ')
    svg_print_float(io, op.ctrl2[2].value)
    print(io, ' ')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

function svg_print_path_op(io::IO, op::CubicCurveShortRelPathOp)
    print(io, 's')
    svg_print_float(io, op.ctrl2[1].value)
    print(io, ' ')
    svg_print_float(io, op.ctrl2[2].value)
    print(io, ' ')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

function svg_print_path_op(io::IO, op::QuadCurveAbsPathOp)
    print(io, 'Q')
    svg_print_float(io, op.ctrl1[1].value)
    print(io, ' ')
    svg_print_float(io, op.ctrl1[2].value)
    print(io, ' ')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

function svg_print_path_op(io::IO, op::QuadCurveRelPathOp)
    print(io, 'q')
    svg_print_float(io, op.ctrl1[1].value)
    print(io, ' ')
    svg_print_float(io, op.ctrl1[2].value)
    print(io, ' ')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

function svg_print_path_op(io::IO, op::QuadCurveShortAbsPathOp)
    print(io, 'T')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

function svg_print_path_op(io::IO, op::QuadCurveShortRelPathOp)
    print(io, 't')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

function svg_print_path_op(io::IO, op::ArcAbsPathOp)
    print(io, 'A')
    svg_print_float(io, op.rx.value)
    print(io, ' ')
    svg_print_float(io, op.ry.value)
    print(io, ' ')
    svg_print_float(io, op.rotation)
    print(io, ' ',
          op.largearc ? 1 : 0, ' ',
          op.sweep ? 1 : 0, ' ')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

function svg_print_path_op(io::IO, op::ArcRelPathOp)
    print(io, 'a')
    svg_print_float(io, op.rx.value)
    print(io, ' ')
    svg_print_float(io, op.ry.value)
    print(io, ' ')
    svg_print_float(io, op.rotation)
    print(io, ' ',
          op.largearc ? 1 : 0, ' ',
          op.sweep ? 1 : 0, ' ')
    svg_print_float(io, op.to[1].value)
    print(io, ' ')
    svg_print_float(io, op.to[2].value)
end

function draw(img::SVG, prim::PathPrimitive, idx::Int)
    indent(img)
    print(img.out, "<path d=\"")
    for op in prim.ops
        svg_print_path_op(img.out, op)
    end
    print(img.out, '"')
    print_vector_properties(img, idx)
    print(img.out, "/>\n")
end


text:
using CUDAdrv, CUDArt
using Base.Test

using Compat

dev = CuDevice(0)
ctx = CuContext(dev)

CUDArt.@compile dev kernel_vadd """
    __global__ void kernel_vadd(const float *a, const float *b, float *c)
    {
        int i = blockIdx.x *blockDim.x + threadIdx.x;
        c[i] = a[i] + b[i];
    }
"""

dims = (3,4)
a = round.(rand(Float32, dims) * 100)
b = round.(rand(Float32, dims) * 100)

d_a = CuArray(a)
d_b = CuArray(b)
d_c = similar(d_a)

len = prod(dims)
cudacall(kernel_vadd, len, 1, Tuple{Ptr{Cfloat},Ptr{Cfloat},Ptr{Cfloat}}, d_a, d_b, d_c)
c = Array(d_c)
@test a+b ≈ c

destroy!(ctx)


text:
replstr(x) = sprint((io,x) -> writemime(io,MIME("text/plain"),x), x)

@test replstr(cell(2)) == "2-element Array{Any,1}:\n #undef\n #undef"
@test replstr(cell(2,2)) == "2x2 Array{Any,2}:\n #undef  #undef\n #undef  #undef"
@test replstr(cell(2,2,2)) == "2x2x2 Array{Any,3}:\n[:, :, 1] =\n #undef  #undef\n #undef  #undef\n\n[:, :, 2] =\n #undef  #undef\n #undef  #undef"

immutable T5589
    names::Vector{UTF8String}
end
@test replstr(T5589(Array(UTF8String,100))) == "T5589(UTF8String[#undef,#undef,#undef,#undef,#undef,#undef,#undef,#undef,#undef,#undef  …  #undef,#undef,#undef,#undef,#undef,#undef,#undef,#undef,#undef,#undef])"

@test replstr(parse("type X end")) == ":(type X\n    end)"
@test replstr(parse("immutable X end")) == ":(immutable X\n    end)"
s = "ccall(:f,Int,(Ptr{Void},),&x)"
@test replstr(parse(s)) == ":($s)"

# expression printing

macro test_repr(x)
    quote
        # Note: We can't just compare x1 and x2 because interpolated
        # strings get converted to string Exprs by the first show().
        # This could produce a few false positives, but until string
        # interpolation works we don't really have a choice.
        local x1 = parse($x)
        local x2 = eval(parse(repr(x1)))
        local x3 = eval(parse(repr(x2)))
        x3 == x2 ? nothing : error(string(
            "repr test failed:",
            "\noriginal: ", $x,
            "\n\nparsed: ", x2, "\n", sprint(dump, x2),
            "\n\nreparsed: ", x3, "\n", sprint(dump, x3)
        ))
    end
end

# primitive types
@test_repr "x"
@test_repr "123"
@test_repr "\"123\""
@test_repr ":()"
@test_repr ":(x, y)"

# basic expressions
@test_repr "x + y"
@test_repr "2e"
@test_repr "!x"
@test_repr "f(1,2,3)"
@test_repr "x = ~y"
@test_repr ":(:x, :y)"
@test_repr ":(:(:(x)))"

# order of operations
@test_repr "x + y * z"
@test_repr "x * y + z"
@test_repr "x * (y + z)"
@test_repr "!x^y"
@test_repr "!x^(y+z)"
@test_repr "!(x^y+z)"
@test_repr "x^-y"
@test_repr "x^-(y+z)"
@test_repr "x^-f(y+z)"
@test_repr "+(w-x)^-f(y+z)"
@test_repr "w = (x = y) = z"
@test_repr "a & b && c"
@test_repr "a & (b && c)"

# precedence tie resolution
@test_repr "(a * b) * (c * d)"
@test_repr "(a / b) / (c / d / e)"
@test_repr "(a == b == c) != (c == d < e)"

# control structures (shamelessly stolen from base/bitarray.jl)
@test_repr """type BitArray{N} <: AbstractArray{Bool, N}
    chunks::Vector{Uint64}
    len::Int
    dims::NTuple{N,Int}
    function BitArray(dims::Int...)
        if length(dims) != N
            error(\"number of dimensions must be \$N (got \$(length(dims)))\")
        end
        n = 1
        for d in dims
            if d < 0
                error(\"dimension size must be nonnegative (got \$d)\")
            end
            n *= d
        end
        nc = num_bit_chunks(n)
        chunks = Array(Uint64, nc)
        if nc > 0
            chunks[end] = uint64(0)
        end
        b = new(chunks, n)
        if N != 1
            b.dims = dims
        end
        return b
    end
end"""

@test_repr """function copy_chunks(dest::Vector{Uint64}, pos_d::Integer, src::Vector{Uint64}, pos_s::Integer, numbits::Integer)
    if numbits == 0
        return
    end
    if dest === src && pos_d > pos_s
        return copy_chunks_rtol(dest, pos_d, pos_s, numbits)
    end

    kd0, ld0 = get_chunks_id(pos_d)
    kd1, ld1 = get_chunks_id(pos_d + numbits - 1)
    ks0, ls0 = get_chunks_id(pos_s)
    ks1, ls1 = get_chunks_id(pos_s + numbits - 1)

    delta_kd = kd1 - kd0
    delta_ks = ks1 - ks0

    u = _msk64
    if delta_kd == 0
        msk_d0 = ~(u << ld0) | (u << ld1 << 1)
    else
        msk_d0 = ~(u << ld0)
        msk_d1 = (u << ld1 << 1)
    end
    if delta_ks == 0
        msk_s0 = (u << ls0) & ~(u << ls1 << 1)
    else
        msk_s0 = (u << ls0)
    end

    chunk_s0 = glue_src_bitchunks(src, ks0, ks1, msk_s0, ls0)

    dest[kd0] = (dest[kd0] & msk_d0) | ((chunk_s0 << ld0) & ~msk_d0)

    if delta_kd == 0
        return
    end

    for i = 1 : kd1 - kd0 - 1
        chunk_s1 = glue_src_bitchunks(src, ks0 + i, ks1, msk_s0, ls0)

        chunk_s = (chunk_s0 >>> (63 - ld0) >>> 1) | (chunk_s1 << ld0)

        dest[kd0 + i] = chunk_s

        chunk_s0 = chunk_s1
    end

    if ks1 >= ks0 + delta_kd
        chunk_s1 = glue_src_bitchunks(src, ks0 + delta_kd, ks1, msk_s0, ls0)
    else
        chunk_s1 = uint64(0)
    end

    chunk_s = (chunk_s0 >>> (63 - ld0) >>> 1) | (chunk_s1 << ld0)

    dest[kd1] = (dest[kd1] & msk_d1) | (chunk_s & ~msk_d1)

    return
end"""

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using JSON, JLD, DSP

# Load and normalize data within [0,1].

function loaddata(fmin=5800, fmax=7300)
    data  = Any[]
    files = readdir(".")
    for f in files
        ismatch(r"^\w+\.json$", f) || continue
        println(f)
        dict = JSON.parsefile(f)
        for (k,v) in dict
            if isa(v,Void)
                dict[k] = Any[]
            else
                for i=1:length(v)
                    v[i] = (convert(Array{Float32},v[i]) - fmin) ./ (fmax-fmin)
                end
            end
        end
        push!(data, dict)
    end
    return data
end

# Save Julia version:
# JLD.save("groupedMelSegData.jld","data",loaddata())

# subsample using DSP.resample (takes long)
function subsample(data, rate)
    sdata = Any[]
    for dict in data
        sdict = Dict()
        for (k,v) in dict
            svecs = Any[]
            for a in v
                push!(svecs, resample(a, rate))
            end
            sdict[k] = svecs
        end
        push!(sdata, sdict)
    end
    return sdata
end

# Generate positive and negative examples: Returns a pair of matrices,
# one for true one for false examples.  The columns of matrices
# contain the pair concatenated.
function pairdata(data)
    tout = Float32[]
    fout = Float32[]
    for d in data
        r = d["RefSegsTrue"]
        t = d["PerSegsTrue"]
        f = d["PerSegsFalse"]
        for i=1:length(r)
            for j=1:length(t)
                append!(tout, r[i])
                append!(tout, t[j])
            end
        end
        for i=1:length(r)
            for j=1:length(f)
                append!(fout, r[i])
                append!(fout, f[j])
            end
        end
    end
    rows = 2*length(data[1]["RefSegsTrue"][1])
    cols = div(length(tout),rows)
    tout = reshape(tout, (rows,cols))
    cols = div(length(fout),rows)
    fout = reshape(fout, (rows,cols))
    return (tout, fout)
end

# Split the pair data into a balanced train, dev and test set with a given minibatch size.
function trntst(pdata; batch=100, splt=((20000,50000),(1700,4300)))
    (xt,xf) = pdata
    (nd,nt) = size(xt)
    (nd,nf) = size(xf)
    rt = randperm(nt)
    rf = randperm(nf)
    nt = nf = 0
    bdata = Any[]
    for (t,f) in splt
        x = hcat(xt[:,rt[nt+1:nt+t]], xf[:,rf[nf+1:nf+f]])
        y = hcat(ones(Float32,1,t), -ones(Float32,1,f))
        nt += t; nf += f
        r = randperm(t+f)
        batches = Any[]
        for i=1:batch:length(r)
            xbatch = x[:,r[i:i+batch-1]]
            ybatch = y[:,r[i:i+batch-1]]
            push!(batches, (xbatch, ybatch))
        end
        push!(bdata, batches)
    end
    return bdata
end

# To load and save KnetArrays:
using Knet
import JLD: writeas, readas
type _KnetArray; a::Array; end
writeas(c::KnetArray) = _KnetArray(Array(c))
readas(d::_KnetArray) = KnetArray(d.a)

# bdata = trntst(pairdata(loaddata()))

# function mat32(x)
#     isa(x,Void) && return x
#     cols = length(x)
#     rows = length(x[1])
#     m = Array(Float32, rows, cols)
#     for i = 1:length(x)
#         m[:,i] = x[i]
#     end
#     return m
# end

# Splitting the melodies:
# @load "groupedMelSegData.jld" # loads data
# trnpairs = pairdata(data[1:30])  # ((3644,17434),(3644,46135))
# tstpairs = pairdata(data[31:40]) # ((3644,4279),(3644,9506))
# trndata = trntst(trnpairs; splt=((17400,46100),))[1]
# tstdata = trntst(tstpairs; splt=((4200,9500),))[1]
# #trndata = map(a->(KnetArray(a[1]),KnetArray(a[2])), trndata)
# #tstdata = map(a->(KnetArray(a[1]),KnetArray(a[2])), tstdata)
# mlprun([trndata,tstdata]; epochs=100)

max_stars_repo_name: JuliaTagBot/melseg

max_issues_repo_name: JuliaTagBot/melseg

max_forks_repo_name: JuliaTagBot/melseg


text:
using Jute

f1 = @global_fixture for x in 1:3
    @produce x
end

f2 = @local_fixture begin
    @produce 1
end

f3 = [1, 2, 3]

for i in 1:100
    @testcase "t$i" for x in f1, y in f2, z in f3
        @test 1 == 1
    end
end


text:
export CERK, cerk!, IERK, ierk!

abstract type AbstractRK{T} end

#### Classic Explicit Runge-Kutta (CERK) ####
struct CERK <: AbstractRK{Float64}
    c2::Float64
    c3::Float64
    c4::Float64
    a2::Float64
    a3::Float64
    a4::Float64
    bv::NTuple{4,Float64}

    CERK(Δt::Real) = new(0.5*Δt, 0.5*Δt, Δt, 0.5*Δt, 0.5*Δt, Δt, (Δt/6.0, Δt/3.0, Δt/3.0, Δt/6.0))
end

function cerk!(ks::MatO{Tk}, up::VecI{Tu}, un::VecI{Tu}, tn::Real, diff!::Function, cerk::CERK; diff_args...) where {Tk<:Real,Tu<:Real}
    #### initialize ####
    @nview ks 4
    one2n = eachindex(un)
    #### generate k1 ####
    copy!(un, up)
    diff!(k1, un, tn; diff_args...)
    #### generate k2 ####
    @simd for i in one2n
        @inbounds un[i] = up[i] + cerk.a2 * k1[i]
    end
    diff!(k2, un, tn + cerk.c2; diff_args...)
    #### generate k3 ####
    @simd for i in one2n
        @inbounds un[i] = up[i] + cerk.a3 * k2[i]
    end
    diff!(k3, un, tn + cerk.c3; diff_args...)
    #### generate k4 ####
    @simd for i in one2n
        @inbounds un[i] = up[i] + cerk.a4 * k3[i]
    end
    diff!(k4, un, tn + cerk.c4; diff_args...)
    #### update un ####
    copy!(un, up)
    gemv!(TXN, 1, ks, cerk.bv, 1, un)
end

function cerk!(ks::MatO{Tk}, up::VecI{Tu}, un::VecI{Tu}, diff!::Function, cerk::CERK; diff_args...) where {Tk<:Real,Tu<:Real}
    #### initialize ####
    @nview ks 4
    one2n = eachindex(un)
    #### generate k1 ####
    copy!(un, up)
    diff!(k1, un; diff_args...)
    #### generate k2 ####
    @simd for i in one2n
        @inbounds un[i] = up[i] + cerk.a2 * k1[i]
    end
    diff!(k2, un; diff_args...)
    #### generate k3 ####
    @simd for i in one2n
        @inbounds un[i] = up[i] + cerk.a3 * k2[i]
    end
    diff!(k3, un; diff_args...)
    #### generate k4 ####
    @simd for i in one2n
        @inbounds un[i] = up[i] + cerk.a4 * k3[i]
    end
    diff!(k4, un; diff_args...)
    #### update un ####
    copy!(un, up)
    gemv!(TXN, 1, ks, cerk.bv, 1, un)
end

#### Improved Explicit Runge-Kutta (IERK) ####
struct IERK <: AbstractRK{Float64}
    c2::Float64
    c3::Float64
    c4::Float64
    a2::Float64
    a3::NTuple{2,Float64}
    a4::NTuple{3,Float64}
    bv::NTuple{4,Float64}

    IERK(Δt::Real) = new(
        0.4*Δt, 0.45573725*Δt, Δt, 0.4*Δt, (0.29697761*Δt,  0.15875964*Δt),
        (0.21810040*Δt, -3.05096516*Δt, 3.83286476*Δt), (0.17476028*Δt, -0.55148066*Δt, 1.20553560*Δt, 0.17118478*Δt)
    )
end

function ierk!(ks::MatO{Tk}, up::VecI{Tu}, un::VecI{Tu}, tn::Real, diff!::Function, ierk::IERK; diff_args...) where {Tk<:Real,Tu<:Real}
    #### initialize ####
    @nview ks 4
    #### generate k1 ####
    copy!(un, up)
    diff!(k1, un, tn; diff_args...)
    #### generate k2 ####
    @simd for i in eachindex(un)
        @inbounds un[i] = up[i] + ierk.a2 * k1[i]
    end
    diff!(k2, un, tn + ierk.c2; diff_args...)
    #### generate k3 ####
    copy!(un, up)
    gemv!(TXN, 1, ks, ierk.a3, 1, un)
    diff!(k3, un, tn + ierk.c3; diff_args...)
    #### generate k4 ####
    copy!(un, up)
    gemv!(TXN, 1, ks, ierk.a4, 1, un)
    diff!(k4, un, tn + ierk.c4; diff_args...)
    #### update un ####
    copy!(un, up)
    gemv!(TXN, 1, ks, ierk.bv, 1, un)
end

function ierk!(ks::MatO{Tk}, up::VecI{Tu}, un::VecI{Tu}, diff!::Function, ierk::IERK; diff_args...) where {Tk<:Real,Tu<:Real}
    #### initialize ####
    @nview ks 4
    #### generate k1 ####
    copy!(un, up)
    diff!(k1, un; diff_args...)
    #### generate k2 ####
    @simd for i in eachindex(un)
        @inbounds un[i] = up[i] + ierk.a2 * k1[i]
    end
    diff!(k2, un; diff_args...)
    #### generate k3 ####
    copy!(un, up)
    gemv!(TXN, 1, ks, ierk.a3, 1, un)
    diff!(k3, un; diff_args...)
    #### generate k4 ####
    copy!(un, up)
    gemv!(TXN, 1, ks, ierk.a4, 1, un)
    diff!(k4, un; diff_args...)
    #### update un ####
    copy!(un, up)
    gemv!(TXN, 1, ks, ierk.bv, 1, un)
end


text:
if joinpath(@__DIR__, "..") ∉ LOAD_PATH
    push!(LOAD_PATH, joinpath(@__DIR__, ".."))
end
using Documenter
using OperatorFlux

makedocs(
    modules = [OperatorFlux],
    sitename = "OperatorFlux.jl",
    authors = "CliMA",
    format = Documenter.HTML(
        prettyurls = !isempty(get(ENV, "CI", "")),
        collapselevel = 1,
        mathengine = MathJax3(),
    ),
    pages = [
        "Home" => "index.md",
        "API" => "api.md",
        # "Installation instructions" => "installation_instructions.md",
        # "Running instructions" => "running_instructions.md",
        # "Contributor Guide" => "contributor_guide.md",
        # "Function Index" => "function_index.md",
    ],
)

deploydocs(repo = "github.com/bischtob/OperatorFlux.jl.git", devbranch = "main")


text:
module AWSCorePrototype

export @service, Services

using AWSCore
using Retry

#include("AWSCorePrototypeServices.jl")
include("AWSMetadata.jl")

macro service(module_name::Symbol)
    service_name = "services/" * lowercase(string(module_name)) * ".jl"

    return Expr(:toplevel,
    :(module ($(esc(module_name)))
        Base.include($(esc(module_name)), $(esc(service_name)))
     end))
end

mutable struct RestXMLService
    name::String
    api_version::String
end

mutable struct QueryService
    name::String
    api_version::String
end

mutable struct JSONService
    name::String
    api_version::String

    json_version::String
    target::String
end

mutable struct RestJSONService
    name::String
    api_version::String
end

function (service::RestXMLService)(aws::AWSConfig, request_method::String, request_uri::String, args=[])
    request = Dict{Symbol,Any}(args)
    args = stringdict(request[:args])

    request[:headers] = Dict{String,String}(get(args, "headers", []))
    delete!(args, "headers")
    request[:content] = get(args, "Body", "")
    delete!(args, "Body")

    request[:resource] = rest_resource(request, args)

    query_str  = HTTP.escapeuri(args)

    if query_str  != ""
        if occursin("?", request[:resource])
            request[:resource] *= "&$query_str"
        else
            request[:resource] *= "?$query_str"
        end
    end

    #FIXME deal with bucket prefix
    request[:url] = service_url(aws, request)

    do_request(merge(request, aws))
end
(service::RestXMLService)(request_method::String, request_uri::String, args=[]) = service(default_aws_config(), request_method, request_uri, args)
(service::RestXMLService)(a...; b...) = service(a..., b)

function (service::QueryService)(aws::AWSConfig, operation, args=[])
    request = Dict{Symbol,Any}(args)

    request[:verb] = "POST"
    request[:resource] = get(aws, :resource, "/")
    request[:url] = service_url(aws, request)
    request[:headers] = Dict("Content-Type" => "application/x-www-form-urlencoded; charset=utf-8")

    request[:query] = aws_args_dict(request[:args])
    request[:query]["Action"] = request[:operation]
    request[:query]["Version"] = request[:version]

    if request[:service] == "iam"
        aws = merge(aws, Dict(:region => "us-east-1"))
    end

    if request[:service] in ["iam", "sts", "sqs", "sns"]
        request[:query]["ContentType"] = "JSON"
    end

    request[:content] = HTTP.escapeuri(flatten_query(request[:service], request[:query]))
    do_request(merge(request, aws))
end
(service::QueryService)(operation, args=[]) = service(default_aws_config(), operation, args)
(service::QueryService)(a...; b...) = service(a..., b)

function (service::JSONService)(aws::AWSConfig, operation, args=[])
    request = Dict{Symbol,Any}(args)

    request[:verb] = "POST"
    request[:resource] = "/"
    request[:url] = service_url(aws, request)

    request[:headers] = Dict(
        "Content-Type" => "application/x-amz-json-$(request[:json_version])",
        "X-Amz-Target" => "$(request[:target]).$(request[:operation])"
    )

    request[:content] = json(aws_args_dict(request[:args]))

    do_request(merge(request, aws))
end
(service::JSONService)(operation, args=[]) = service(default_aws_config(), operation, args)
(service::JSONService)(a...; b...) = service(a..., b)

function (service::RestJSONService)(aws::AWSConfig, request_method::String, request_uri::String, args=[])
    request = Dict{Symbol,Any}(args)
    args = Dict(request[:args])

    request[:resource] = rest_resource(request, args)
    request[:url] = service_url(aws, request)

    request[:headers] = Dict{String,String}(get(args, "headers", []))
    delete!(args, "headers")
    request[:headers]["Content-Type"] = "application/json"
    request[:content] = json(aws_args_dict(args))

    do_request(merge(request, aws))
end
(service::RestJSONService)(request_method::String, request_uri::String, args=[]) = service(default_aws_config(), request_method, request_uri, args)
(service::RestJSONService)(a...; b...) = service(a..., b)

####################
# Helper Functions #
####################
function aws_args_dict(args)
    result = stringdict(args)

    dictlike(t) = (t <: AbstractDict || t <: Vector && t.parameters[1] <: Pair{String})

    for (k, v) in result
        if dictlike(typeof(v))
            result[k] = aws_args_dict(v)
        elseif isa(v, Vector)
            result[k] = [dictlike(typeof(i)) ? aws_args_dict(i) : i for i in v]
        end
    end

    return result
end

function rest_resource(request, args)
    r = request[:resource]

    for (k,v) in args
        if occursin("{$k}", r)
            r = replace(r, "{$k}" => v)
            delete!(args, k)
        elseif occursin("{$k+}", r)
            r = replace(r, "{$k+}" => HTTP.escapepath(v))
            delete!(args, k)
        end
    end

    return r
end

function service_url(aws, request)
    endpoint = get(request, :endpoint, request[:service])
    region = aws[:region]
    regionless_endpoints = ("iam", "route53")

    if endpoint in regionless_endpoints || (endpoint == "sdb" && region == "us-east-1")
        region = ""
    end

    service_host = get(aws, :service_host, "amazonaws.com")

    if !isempty(region)
        return string("https://", endpoint, ".", region, ".", service_host, request[:resource])
    else
        return string("https://", endpoint, ".", service_host, request[:resource])
    end
end

function flatten_query(service, query, prefix="")

    result = Dict{String,String}()

    for (k, v) in query

        if typeof(v) <: AbstractDict

            merge!(result, flatten_query(service, v, "$prefix$k."))

        elseif typeof(v) <: Array

            for (i, x) in enumerate(v)

                suffix = service in ["ec2", "sqs"] ? "" : ".member"
                pk = "$prefix$k$suffix.$i"

                if typeof(x) <: AbstractDict
                    merge!(result, flatten_query(service, x, "$pk."))
                else
                    result[pk] = string(x)
                end
            end
        else
            result["$prefix$k"] = string(v)
        end
    end

    return result
end

function do_request(r::AWSRequest)

    response = nothing

    # Try request 3 times to deal with possible Redirect and ExiredToken...
    @repeat 3 try

        # Default headers...
        if !haskey(r, :headers)
            r[:headers] = Dict{String,String}()
        end
        r[:headers]["User-Agent"] = "AWSCore.jl/0.0.0"
        r[:headers]["Host"]       = HTTP.URI(r[:url]).host

        # If existing, use credentials to sign request...
        r[:creds] === nothing || sign!(r)

        if debug_level > 0
            dump_aws_request(r)
        end

        # Send the request...
        response = http_request(r)

        if response.status in [301, 302, 307] &&
           HTTP.header(response, "Location") != ""
            r[:url] = HTTP.header(response, "Location")
            continue
        end

    catch e
        if e isa HTTP.StatusError
            e = AWSException(e)
        end

        # Handle expired signature...
        @retry if :message in fieldnames(typeof(e)) &&
                  occursin(r"Signature expired", e.message)
            if debug_level > 1
                println("Caught $e during request $(dump_aws_request(r)), retrying due to expired signature...")
            end
        end

        # Handle ExpiredToken...
        # See `credsExpiredCodes` in
        # https://github.com/aws/aws-sdk-go/blob/master/aws/request/retryer.go
        @retry if ecode(e) in ("ExpiredToken",
                               "ExpiredTokenException",
                               "RequestExpired")

            # Reload local system credentials if needed...
            check_credentials(r[:creds], force_refresh=true)
            if debug_level > 1
                println("Caught $e during request $(dump_aws_request(r)), retrying due to expired credentials...")
            end
        end

        # Handle throttling
        # see botocore for list of codes:
        # https://github.com/boto/botocore/blob/master/botocore/data/_retry.json
        # Recommended for SDKs at:
        # https://docs.aws.amazon.com/general/latest/gr/api-retries.html
        @delay_retry if e isa AWSException && (
                  http_status(e.cause) == 429 ||
                  ecode(e) in ("Throttling",
                               "ThrottlingException",
                               "ThrottledException",
                               "RequestThrottledException",
                               "TooManyRequestsException",
                               "ProvisionedThroughputExceededException",
                               "LimitExceededException",
                               "RequestThrottled",
                               "PriorRequestNotComplete"))
            if debug_level > 1
                cause = "throttling"

                if ecode(e) == "PriorRequestNotComplete"
                    cause = ecode(e)
                end
                println("Caught $e during request $(dump_aws_request(r)), retrying due to $cause...")
            end
        end

        # Handle BadDigest error and CRC32 thing
        @retry if e isa AWSException && (
            header(e.cause, "crc32body") == "x-amz-crc32" ||
            ecode(e) in ("BadDigest", "RequestTimeout", "RequestTimeoutException")
        )
            if debug_level > 1
                cause = if header(e.cause, "crc32body") == "x-amz-crc32"
                    "CRC32"
                else
                    ecode(e)
                end
                println("Caught $e during request $(dump_aws_request(r)), retrying due to $cause...")
            end
        end
    end

    if debug_level > 1
        display(response)
        println()
    end

    # For HEAD request, return headers...
    if r[:verb] == "HEAD"
        return Dict(response.headers)
    end

    # Return response stream if requested...
    if get(r, :return_stream, false)
        return r[:response_stream]
    end

    # Return raw data if requested...
    if get(r, :return_raw, false)
        return response.body
    end

    # Parse response data according to mimetype...
    mime = HTTP.header(response, "Content-Type", "")
    if mime == ""
        if length(response.body) > 5 && String(response.body[1:5]) == "<?xml"
            mime = "text/xml"
        end
    end

    body = String(copy(response.body))

    if occursin(r"/xml", mime)
        return parse_xml(body)
    end

    if occursin(r"/x-amz-json-1.[01]$", mime)
        if isempty(response.body)
            return nothing
        end
        if get(r, :ordered_json_dict, true)
            return JSON.parse(body, dicttype=OrderedDict)
        else
            return JSON.parse(body)
        end
    end

    if occursin(r"json$", mime)
        if isempty(response.body)
            return nothing
        end
        if get(r, :ordered_json_dict, true)
            info = JSON.parse(body, dicttype=OrderedDict)
        else
            info = JSON.parse(body)
        end
        @protected try
            action = r[:query]["Action"]
            info = info[action * "Response"]
            info = info[action * "Result"]
        catch e
            @ignore if typeof(e) == KeyError end
        end
        return info
    end

    if occursin(r"^text/", mime)
        return body
    end

    # Return raw data by default...
    return response.body
end

end  # module AWSCorePrototype

max_stars_repo_name: mattBrzezinski/Julia-AWS-Status

max_issues_repo_name: mattBrzezinski/Julia-AWS-Status

max_forks_repo_name: mattBrzezinski/Julia-AWS-Status


text:
# NOTE: Carefully review all the NOTES here.
using Distributed

# This is needed to compile once before loading on all cores.
include(joinpath(@__DIR__, "imports.jl"))

# NOTE: We are using 4 temperatures, so using 4 cores.
addprocs(4)

@everywhere include(joinpath(@__DIR__, "imports.jl"))
@everywhere const rfam = CytofRepFAM.Model

@everywhere function preprocess(y, lower, upper)
  @assert lower < upper
  idx_keep = .!vec(any(y .< lower, dims=2) .|
                   any(y .> upper, dims=2))
  return y[idx_keep, :]
end

@everywhere function run(phi, path_to_data::Vector{String},
                         results_dir, aws_bucket;
                         nburn=10000, nsamps=5000, thin=1,
                         K=20, L=Dict(0=>6, 1=>3),
                         tempers=[1, 1.0003, 1.0006, 1.001],
                         batchsizes=[200, 200, 200], pthin=5,
                         y_lower=-7, y_upper=4)

  @assert length(path_to_data) == length(batchsizes)

  # Print setup
  println("Start time: ", Dates.now())
  println("pid: $(getpid())")
  println("Threads: $(Threads.nthreads())")
  println("pwd: $(pwd())")
  println("Data used:")
  foreach(println, "  " .* path_to_data)
  flush(stdout)

  # TODO: Make y
  function read_and_format(path::String)
    df = CSV.read(path, missingstring="NA")
    mat = Matrix(coalesce.(df, NaN))
    return preprocess(mat, y_lower, y_upper), names(df)
  end
  data = read_and_format.(path_to_data)
  y = [d[1] for d in data]
  println("data sizes: ", size.(y))
  markers = let
    _markers = [d[2] for d in data]
    @assert length(unique(_markers)) == 1
    _markers[1]
  end
  println("markers: $(String.(markers))")

  # Results dir
  mkpath(results_dir)

  # number of temperatures
  ntemps = length(tempers)

  # Repulsive FAM penalty scale
  sim_z = rfam.sim_fn_exp_decay_generator(phi)
  use_repulsive = phi > 0

  function init_state_const_data(y; K, L)
    deltaz_prior = TruncatedNormal(1.0, 0.1, 0.0, Inf)
    d = rfam.Data(y)
    c = rfam.defaultConstants(d, K, L,
                              sig2_prior=InverseGamma(3, 1),
                              delta0_prior=deltaz_prior,
                              delta1_prior=deltaz_prior,
                              alpha_prior=Gamma(1, 1),  # shape, scale
                              # NOTE:  These 3 items are important
                              yQuantiles=[.0, .25, .5],
                              pBounds=[.05, .8, .05],
                              y_grid=collect(range(-10,
                                                   stop=4,
                                                   length=100)),
                              probFlip_Z=1.0,
                              similarity_Z=sim_z)

    # NOTE: To spread weights evenly.
    c.eta_prior = Dict(z => Dirichlet(L[z], 1) for z in 0:1)

    # Initialize state.
    # NOTE: Kmeans init
    # states = [rfam.smartInit(c, d, modelNames="kmeans",
    #                          seed=1 + s, iterMax=100)
    #           for s in 1:ntemps]
    # NOTE: Mclust init
    # VII: spherical, unequal volume
    states = [rfam.smartInit(c, d, modelNames="VII", seed=1 + s)
              for s in 1:ntemps]

    # Create tunenrs
    t = rfam.Tuners(d.y, c.K)

    # Create covriates matrix
    X = rfam.eye(Float64, d.I)

    # Create model constants
    cfs = rfam.ConstantsFS(c)

    # NOTE: Use informative priors for W* and omega, to control feature
    # selection.
    # cfs.W_star_prior = Gamma(1000, 0.1) # shape, scale. Works, but extreme?
    # cfs.W_star_prior = Gamma(100, 0.1) # shape, scale. Not working.
    # cfs.W_star_prior = Gamma(1/c.K, 1) # shape, scale. Not working.
    # cfs.W_star_prior = Gamma(1/c.K, c.K * 10) # shape, scale. Not working.
    # cfs.W_star_prior = Gamma(1, 100) # shape, scale. Not super helpful?
    cfs.W_star_prior = Gamma(10, 1) # shape, scale. Try?
    # similar to p ~ Beta(1, 99)
    cfs.omega_prior = Normal(-5, 1)

    dfs = rfam.DataFS(d, X)
    sfss = [rfam.StateFS{Float64}(s, dfs, omega_prior=cfs.omega_prior,
                                  W_star_prior=cfs.W_star_prior, eps_r=0.05)
            for s in states]
    tfs = rfam.TunersFS(t, states[1], X)

    return Dict(:dfs => dfs, :cfs => cfs, :sfss => sfss, :tfs => tfs, :y => y)
  end

  # Parameters to monitor
  monitor1 = [:theta__Z, :theta__v, :theta__alpha,
              :omega, :r, :theta__lam, :W_star, :theta__eta,
              :theta__W, :theta__delta, :theta__sig2]
  monitor2 = [:theta__y_imputed, :theta__gam]

  # MCMC Specs
  nsamps_to_thin(nsamps::Int, nmcmc::Int) = max(1, div(nmcmc, nsamps))
  thin_samps = thin  # Factor to thin the
                     # primary parameters
  mcmc_iter = nsamps * thin_samps  # Number of MCMC iterations

  # LPML / DIC are computed based on `mcmc_iter` samples

  # Configurations: priors, initial state, data, etc.
  config = init_state_const_data(y, K=K, L=L)

  # Print constants
  println("N: $(config[:dfs].data.N)")
  println("J: $(config[:dfs].data.J)")
  rfam.printConstants(config[:cfs])
  flush(stdout)

  # Fit model
  @time out = rfam.fit_fs_imcmc_pt!(
    config[:cfs], config[:dfs],
    inits=config[:sfss],
    nmcmc=mcmc_iter, nburn=nburn,
    batchsizes=batchsizes,
    prior_thin=pthin,
    tempers=tempers,
    randpair=1.0, swap_freq=1.0,
    Z_marg_lamgam=1.0,
    Z_marg_lamgam_decay_rate=100.0,
    Z_marg_lamgam_min=1.0,
    printFreq=1,
    seed=1,
    computedden=true,
    computeDIC=true,
    computeLPML=true,
    thins=[thin, nsamps_to_thin(10, mcmc_iter)],
    ndden_samps=200,
    time_updates=false,
    save_all_states=true,
    verbose=4)

  # Dump output
  BSON.bson("$(results_dir)/output.bson", out)

  println("Completed! ", Dates.now())

  # Send results to aws
  rfam.s3sync(from=results_dir, to=aws_bucket, tags=`--exclude '*.nfs'`)

  return
end  # simfn

### MAIN ###

# READ COMMAND ARGS
phi = parse(Float64, ARGS[1])
path_to_data = String.(split(ARGS[2], ","))
results_dir = ARGS[3]
aws_bucket = ARGS[4]
istest = Bool(parse(Int, ARGS[5]))

println("phi: $phi")
println("data path: $path_to_data")
println("results dir: $results_dir")
println("AWS bucket: $aws_bucket")
println("is test: $istest")
flush(stdout)

if istest
  nsamps = 20
  nburn = 10
else
  nsamps = 5000
  nburn = 10000
end

@time run(phi, path_to_data, results_dir, aws_bucket,
          nsamps=nsamps, nburn=nburn)


text:
module Journals

include("journalentry.jl")
include("journal.jl")

export
    Journal

end


text:
""""
This function gathers data for making price and other market data predictions.
"""
function gather_prediction_parameters(investor::Investor,
                                       sys_data_dir::String,
                                       iteration_year::Int64)
    investor_name = get_name(investor)
    investor_dir = get_data_dir(investor)

    load_dir = joinpath(investor_dir, "timeseries_data_files", "Load")
    dir_exists(load_dir)
    cp(joinpath(sys_data_dir, "timeseries_data_files", "Load", "rep_load_$(iteration_year - 1).csv"),
    joinpath(load_dir, "load_$(iteration_year - 1).csv"), force = true)

    reserve_dir = joinpath(investor_dir, "timeseries_data_files", "Reserves")
    dir_exists(reserve_dir)
    cp(joinpath(sys_data_dir, "timeseries_data_files", "Reserves", "rep_reserve_up_$(iteration_year - 1).csv"),
    joinpath(reserve_dir, "reserve_up_$(iteration_year - 1).csv"), force = true)

    cp(joinpath(sys_data_dir, "timeseries_data_files", "Reserves", "rep_reserve_down_$(iteration_year - 1).csv"),
    joinpath(reserve_dir, "reserve_down_$(iteration_year - 1).csv"), force = true)

    market_names = get_markets(investor)

    rep_hour_weight = get_rep_hour_weight(investor)

    scenarios = get_scenario_data(get_forecast(investor))

    return investor_name, investor_dir, market_names, rep_hour_weight, scenarios
end

"""
This function runs CEM for price predictions based on user-defined parallelization settings.
"""
function create_investor_predictions(investors::Vector{Investor},
                                          active_projects::Vector{Project},
                                          iteration_year::Int64,
                                          yearly_horizon::Int64,
                                          sys_data_dir::String,
                                          sys_results_dir::String,
                                          average_capital_cost_multiplier::Float64,
                                          zones::Vector{String},
                                          lines::Vector{ZonalLine},
                                          peak_load::Float64,
                                          solver::JuMP.MOI.OptimizerWithAttributes,
                                          parallelize_investors::Bool,
                                          parallelize_scenarios::Bool)

    if parallelize_investors

        if parallelize_scenarios

            scenarios_pmap = Scenario[]
            investor_name_pmap = String[]
            investor_dir_pmap = String[]
            market_names_pmap = Vector{Symbol}[]
            rep_hour_weight_pmap = Vector{Float64}[]
            expected_portfolio_pmap = Vector{Project}[]

            for investor in investors

                investor_name,
                investor_dir,
                market_names,
                rep_hour_weight,
                scenarios = gather_prediction_parameters(investor, sys_data_dir, iteration_year)

                for scenario in scenarios
                    push!(scenarios_pmap, scenario)
                    push!(investor_name_pmap, investor_name)
                    push!(investor_dir_pmap, investor_dir)
                    push!(market_names_pmap, market_names)
                    push!(rep_hour_weight_pmap, rep_hour_weight)
                    push!(expected_portfolio_pmap, active_projects)

                end

            end

            num_tasks = length(scenarios_pmap)
            Distributed.pmap(create_expected_marketdata,
                 investor_dir_pmap,
                 market_names_pmap,
                 expected_portfolio_pmap,
                 repeat([zones], num_tasks),
                 repeat([lines], num_tasks),
                 repeat([peak_load], num_tasks),
                 rep_hour_weight_pmap,
                 repeat([average_capital_cost_multiplier], num_tasks),
                 scenarios_pmap,
                 repeat([iteration_year], num_tasks),
                 repeat([yearly_horizon], num_tasks),
                 repeat([solver], num_tasks),
                 repeat([sys_results_dir], num_tasks),
                 investor_name_pmap)

        else

            num_tasks = length(investors)
            Distributed.pmap(parallelize_only_investors,
                            investors,
                            repeat([sys_data_dir], num_tasks),
                            repeat([active_projects], num_tasks),
                            repeat([zones], num_tasks),
                            repeat([lines], num_tasks),
                            repeat([peak_load], num_tasks),
                            repeat([average_capital_cost_multiplier], num_tasks),
                            repeat([iteration_year], num_tasks),
                            repeat([yearly_horizon], num_tasks),
                            repeat([solver], num_tasks),
                            repeat([sys_results_dir], num_tasks),
                            get_name.(investors))
        end

    else

        for investor in investors

            investor_name,
            investor_dir,
            market_names,
            rep_hour_weight,
            scenarios = gather_prediction_parameters(investor, sys_data_dir, iteration_year)

            if parallelize_scenarios

                num_scenarios = length(scenarios)

                Distributed.pmap(create_expected_marketdata,
                    repeat([investor_dir], num_scenarios),
                    repeat([market_names], num_scenarios),
                    repeat([active_projects], num_scenarios),
                    repeat([zones], num_scenarios),
                    repeat([lines], num_scenarios),
                    repeat([peak_load], num_scenarios),
                    repeat([rep_hour_weight], num_scenarios),
                    repeat([average_capital_cost_multiplier], num_scenarios),
                    scenarios,
                    repeat([iteration_year], num_scenarios),
                    repeat([yearly_horizon], num_scenarios),
                    repeat([solver], num_scenarios),
                    repeat([sys_results_dir], num_scenarios),
                    repeat([investor_name], num_scenarios))

            else

                for scenario in scenarios
                    create_expected_marketdata(investor_dir,
                                            market_names,
                                            active_projects,
                                            zones,
                                            lines,
                                            peak_load,
                                            rep_hour_weight,
                                            average_capital_cost_multiplier,
                                            scenario,
                                            iteration_year,
                                            yearly_horizon,
                                            solver,
                                            sys_results_dir,
                                            investor_name)
                end

            end

        end

    end

    return
end


text:
using ProtoBuf
import ProtoBuf.meta

type BinaryOpReq
    i1::Int64
    i2::Int64
    BinaryOpReq() = (o=new(); fillunset(o); o)
end #type BinaryOpReq
const __req_BinaryOpReq = Symbol[:i1,:i2]
meta(t::Type{BinaryOpReq}) = meta(t, __req_BinaryOpReq, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)

type BinaryOpResp
    result::Int64
    BinaryOpResp() = (o=new(); fillunset(o); o)
end #type BinaryOpResp
const __req_BinaryOpResp = Symbol[:result]
meta(t::Type{BinaryOpResp}) = meta(t, __req_BinaryOpResp, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)

# service methods for TestMath
const _TestMath_methods = MethodDescriptor[
        MethodDescriptor("Mul", 1, BinaryOpReq, BinaryOpResp),
        MethodDescriptor("Add", 2, BinaryOpReq, BinaryOpResp)
    ] # const _TestMath_methods
const _TestMath_desc = ServiceDescriptor("TestMath", 1, _TestMath_methods)

TestMath(impl::Module) = ProtoService(_TestMath_desc, impl)

type TestMathStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    TestMathStub(channel::ProtoRpcChannel) = new(ProtoServiceStub(_TestMath_desc, channel))
end # type TestMathStub

type TestMathBlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    TestMathBlockingStub(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_TestMath_desc, channel))
end # type TestMathBlockingStub

Mul(stub::TestMathStub, controller::ProtoRpcController, inp::BinaryOpReq, done::Function) = call_method(stub.impl, _TestMath_methods[1], controller, inp, done)
Mul(stub::TestMathBlockingStub, controller::ProtoRpcController, inp::BinaryOpReq) = call_method(stub.impl, _TestMath_methods[1], controller, inp)

Add(stub::TestMathStub, controller::ProtoRpcController, inp::BinaryOpReq, done::Function) = call_method(stub.impl, _TestMath_methods[2], controller, inp, done)
Add(stub::TestMathBlockingStub, controller::ProtoRpcController, inp::BinaryOpReq) = call_method(stub.impl, _TestMath_methods[2], controller, inp)

export BinaryOpReq, BinaryOpResp, TestMath, TestMathStub, TestMathBlockingStub, Mul, Add


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## reductions ##

###### Generic (map)reduce functions ######

if Int === Int32
    const SmallSigned = Union{Int8,Int16}
    const SmallUnsigned = Union{UInt8,UInt16}
else
    const SmallSigned = Union{Int8,Int16,Int32}
    const SmallUnsigned = Union{UInt8,UInt16,UInt32}
end

"""
    Base.add_sum(x, y)

The reduction operator used in `sum`. The main difference from [`+`](@ref) is that small
integers are promoted to `Int`/`UInt`.
"""
add_sum(x, y) = x + y
add_sum(x::SmallSigned, y::SmallSigned) = Int(x) + Int(y)
add_sum(x::SmallUnsigned, y::SmallUnsigned) = UInt(x) + UInt(y)
add_sum(x::Real, y::Real)::Real = x + y

"""
    Base.mul_prod(x, y)

The reduction operator used in `prod`. The main difference from [`*`](@ref) is that small
integers are promoted to `Int`/`UInt`.
"""
mul_prod(x, y) = x * y
mul_prod(x::SmallSigned, y::SmallSigned) = Int(x) * Int(y)
mul_prod(x::SmallUnsigned, y::SmallUnsigned) = UInt(x) * UInt(y)
mul_prod(x::Real, y::Real)::Real = x * y

## foldl && mapfoldl

@noinline function mapfoldl_impl(f, op, nt::NamedTuple{(:init,)}, itr, i...)
    init = nt.init
    # Unroll the while loop once; if init is known, the call to op may
    # be evaluated at compile time
    y = iterate(itr, i...)
    y === nothing && return init
    v = op(init, f(y[1]))
    while true
        y = iterate(itr, y[2])
        y === nothing && break
        v = op(v, f(y[1]))
    end
    return v
end

function mapfoldl_impl(f, op, nt::NamedTuple{()}, itr)
    y = iterate(itr)
    if y === nothing
        return Base.mapreduce_empty_iter(f, op, itr, IteratorEltype(itr))
    end
    (x, i) = y
    init = mapreduce_first(f, op, x)
    return mapfoldl_impl(f, op, (init=init,), itr, i)
end

"""
    mapfoldl(f, op, itr; [init])

Like [`mapreduce`](@ref), but with guaranteed left associativity, as in [`foldl`](@ref).
If provided, the keyword argument `init` will be used exactly once. In general, it will be
necessary to provide `init` to work with empty collections.
"""
mapfoldl(f, op, itr; kw...) = mapfoldl_impl(f, op, kw.data, itr)

"""
    foldl(op, itr; [init])

Like [`reduce`](@ref), but with guaranteed left associativity. If provided, the keyword
argument `init` will be used exactly once. In general, it will be necessary to provide
`init` to work with empty collections.

# Examples
```jldoctest
julia> foldl(=>, 1:4)
((1=>2)=>3) => 4

julia> foldl(=>, 1:4; init=0)
(((0=>1)=>2)=>3) => 4
```
"""
foldl(op, itr; kw...) = mapfoldl(identity, op, itr; kw...)

## foldr & mapfoldr

function mapfoldr_impl(f, op, nt::NamedTuple{(:init,)}, itr, i::Integer)
    init = nt.init
    # Unroll the while loop once; if init is known, the call to op may
    # be evaluated at compile time
    if isempty(itr) || i == 0
        return init
    else
        x = itr[i]
        v  = op(f(x), init)
        while i > 1
            x = itr[i -= 1]
            v = op(f(x), v)
        end
        return v
    end
end

function mapfoldr_impl(f, op, ::NamedTuple{()}, itr, i::Integer)
    if isempty(itr)
        return Base.mapreduce_empty_iter(f, op, itr, IteratorEltype(itr))
    end
    return mapfoldr_impl(f, op, (init=mapreduce_first(f, op, itr[i]),), itr, i-1)
end

"""
    mapfoldr(f, op, itr; [init])

Like [`mapreduce`](@ref), but with guaranteed right associativity, as in [`foldr`](@ref). If
provided, the keyword argument `init` will be used exactly once. In general, it will be
necessary to provide `init` to work with empty collections.
"""
mapfoldr(f, op, itr; kw...) = mapfoldr_impl(f, op, kw.data, itr, lastindex(itr))

"""
    foldr(op, itr; [init])

Like [`reduce`](@ref), but with guaranteed right associativity. If provided, the keyword
argument `init` will be used exactly once. In general, it will be necessary to provide
`init` to work with empty collections.

# Examples
```jldoctest
julia> foldr(=>, 1:4)
1 => (2=>(3=>4))

julia> foldr(=>, 1:4; init=0)
1 => (2=>(3=>(4=>0)))
```
"""
foldr(op, itr; kw...) = mapfoldr(identity, op, itr; kw...)

## reduce & mapreduce

# `mapreduce_impl()` is called by `mapreduce()` (via `_mapreduce()`, when `A`
# supports linear indexing) and does actual calculations (for `A[ifirst:ilast]` subset).
# For efficiency, no parameter validity checks are done, it's the caller's responsibility.
# `ifirst:ilast` range is assumed to be a valid non-empty subset of `A` indices.

# This is a generic implementation of `mapreduce_impl()`,
# certain `op` (e.g. `min` and `max`) may have their own specialized versions.
@noinline function mapreduce_impl(f, op, A::AbstractArray, ifirst::Integer, ilast::Integer, blksize::Int)
    if ifirst == ilast
        @inbounds a1 = A[ifirst]
        return mapreduce_first(f, op, a1)
    elseif ifirst + blksize > ilast
        # sequential portion
        @inbounds a1 = A[ifirst]
        @inbounds a2 = A[ifirst+1]
        v = op(f(a1), f(a2))
        @simd for i = ifirst + 2 : ilast
            @inbounds ai = A[i]
            v = op(v, f(ai))
        end
        return v
    else
        # pairwise portion
        imid = (ifirst + ilast) >> 1
        v1 = mapreduce_impl(f, op, A, ifirst, imid, blksize)
        v2 = mapreduce_impl(f, op, A, imid+1, ilast, blksize)
        return op(v1, v2)
    end
end

mapreduce_impl(f, op, A::AbstractArray, ifirst::Integer, ilast::Integer) =
    mapreduce_impl(f, op, A, ifirst, ilast, pairwise_blocksize(f, op))

"""
    mapreduce(f, op, itr; [init])

Apply function `f` to each element in `itr`, and then reduce the result using the binary
function `op`. If provided, `init` must be a neutral element for `op` that will be returned
for empty collections. It is unspecified whether `init` is used for non-empty collections.
In general, it will be necessary to provide `init` to work with empty collections.

[`mapreduce`](@ref) is functionally equivalent to calling
`reduce(op, map(f, itr); init=init)`, but will in general execute faster since no
intermediate collection needs to be created. See documentation for [`reduce`](@ref) and
[`map`](@ref).

# Examples
```jldoctest
julia> mapreduce(x->x^2, +, [1:3;]) # == 1 + 4 + 9
14
```

The associativity of the reduction is implementation-dependent. Additionally, some
implementations may reuse the return value of `f` for elements that appear multiple times in
`itr`. Use [`mapfoldl`](@ref) or [`mapfoldr`](@ref) instead for
guaranteed left or right associativity and invocation of `f` for every value.
"""
mapreduce(f, op, itr; kw...) = mapfoldl(f, op, itr; kw...)

# Note: sum_seq usually uses four or more accumulators after partial
# unrolling, so each accumulator gets at most 256 numbers
pairwise_blocksize(f, op) = 1024

# This combination appears to show a benefit from a larger block size
pairwise_blocksize(::typeof(abs2), ::typeof(+)) = 4096

# handling empty arrays
_empty_reduce_error() = throw(ArgumentError("reducing over an empty collection is not allowed"))

"""
    Base.reduce_empty(op, T)

The value to be returned when calling [`reduce`](@ref), [`foldl`](@ref) or [`foldr`](@ref)
with reduction `op` over an empty array with element type of `T`.

If not defined, this will throw an `ArgumentError`.
"""
reduce_empty(op, T) = _empty_reduce_error()
reduce_empty(::typeof(+), T) = zero(T)
reduce_empty(::typeof(+), ::Type{Bool}) = zero(Int)
reduce_empty(::typeof(*), T) = one(T)
reduce_empty(::typeof(*), ::Type{<:AbstractChar}) = ""
reduce_empty(::typeof(&), ::Type{Bool}) = true
reduce_empty(::typeof(|), ::Type{Bool}) = false

reduce_empty(::typeof(add_sum), T) = reduce_empty(+, T)
reduce_empty(::typeof(add_sum), ::Type{T}) where {T<:SmallSigned}  = zero(Int)
reduce_empty(::typeof(add_sum), ::Type{T}) where {T<:SmallUnsigned} = zero(UInt)
reduce_empty(::typeof(mul_prod), T) = reduce_empty(*, T)
reduce_empty(::typeof(mul_prod), ::Type{T}) where {T<:SmallSigned}  = one(Int)
reduce_empty(::typeof(mul_prod), ::Type{T}) where {T<:SmallUnsigned} = one(UInt)

"""
    Base.mapreduce_empty(f, op, T)

The value to be returned when calling [`mapreduce`](@ref), [`mapfoldl`](@ref`) or
[`mapfoldr`](@ref) with map `f` and reduction `op` over an empty array with element type
of `T`.

If not defined, this will throw an `ArgumentError`.
"""
mapreduce_empty(f, op, T) = _empty_reduce_error()
mapreduce_empty(::typeof(identity), op, T) = reduce_empty(op, T)
mapreduce_empty(::typeof(abs), op, T)      = abs(reduce_empty(op, T))
mapreduce_empty(::typeof(abs2), op, T)     = abs2(reduce_empty(op, T))

mapreduce_empty(f::typeof(abs),  ::typeof(max), T) = abs(zero(T))
mapreduce_empty(f::typeof(abs2), ::typeof(max), T) = abs2(zero(T))

mapreduce_empty_iter(f, op, itr, ::HasEltype) = mapreduce_empty(f, op, eltype(itr))
mapreduce_empty_iter(f, op::typeof(&), itr, ::EltypeUnknown) = true
mapreduce_empty_iter(f, op::typeof(|), itr, ::EltypeUnknown) = false
mapreduce_empty_iter(f, op, itr, ::EltypeUnknown) = _empty_reduce_error()

# handling of single-element iterators
"""
    Base.reduce_first(op, x)

The value to be returned when calling [`reduce`](@ref), [`foldl`](@ref`) or
[`foldr`](@ref) with reduction `op` over an iterator which contains a single element
`x`. This value may also used to initialise the recursion, so that `reduce(op, [x, y])`
may call `op(reduce_first(op, x), y)`.

The default is `x` for most types. The main purpose is to ensure type stability, so
additional methods should only be defined for cases where `op` gives a result with
different types than its inputs.
"""
reduce_first(op, x) = x
reduce_first(::typeof(+), x::Bool) = Int(x)
reduce_first(::typeof(*), x::AbstractChar) = string(x)

reduce_first(::typeof(add_sum), x) = reduce_first(+, x)
reduce_first(::typeof(add_sum), x::SmallSigned)   = Int(x)
reduce_first(::typeof(add_sum), x::SmallUnsigned) = UInt(x)
reduce_first(::typeof(mul_prod), x) = reduce_first(*, x)
reduce_first(::typeof(mul_prod), x::SmallSigned)   = Int(x)
reduce_first(::typeof(mul_prod), x::SmallUnsigned) = UInt(x)

"""
    Base.mapreduce_first(f, op, x)

The value to be returned when calling [`mapreduce`](@ref), [`mapfoldl`](@ref`) or
[`mapfoldr`](@ref) with map `f` and reduction `op` over an iterator which contains a
single element `x`. This value may also used to initialise the recursion, so that
`mapreduce(f, op, [x, y])` may call `op(reduce_first(op, f, x), f(y))`.

The default is `reduce_first(op, f(x))`.
"""
mapreduce_first(f, op, x) = reduce_first(op, f(x))

_mapreduce(f, op, A::AbstractArray) = _mapreduce(f, op, IndexStyle(A), A)

function _mapreduce(f, op, ::IndexLinear, A::AbstractArray{T}) where T
    inds = LinearIndices(A)
    n = length(inds)
    if n == 0
        return mapreduce_empty(f, op, T)
    elseif n == 1
        @inbounds a1 = A[first(inds)]
        return mapreduce_first(f, op, a1)
    elseif n < 16 # process short array here, avoid mapreduce_impl() compilation
        @inbounds i = first(inds)
        @inbounds a1 = A[i]
        @inbounds a2 = A[i+=1]
        s = op(f(a1), f(a2))
        while i < last(inds)
            @inbounds Ai = A[i+=1]
            s = op(s, f(Ai))
        end
        return s
    else
        return mapreduce_impl(f, op, A, first(inds), last(inds))
    end
end

mapreduce(f, op, a::Number) = mapreduce_first(f, op, a)

_mapreduce(f, op, ::IndexCartesian, A::AbstractArray) = mapfoldl(f, op, A)

"""
    reduce(op, itr; [init])

Reduce the given collection `itr` with the given binary operator `op`. If provided, the
initial value `init` must be a neutral element for `op` that will be returned for empty
collections. It is unspecified whether `init` is used for non-empty collections.

For empty collections, providing `init` will be necessary, except for some special cases
(e.g. when `op` is one of `+`, `*`, `max`, `min`, `&`, `|`) when Julia can determine the
neutral element of `op`.

Reductions for certain commonly-used operators may have special implementations, and
should be used instead: `maximum(itr)`, `minimum(itr)`, `sum(itr)`, `prod(itr)`,
 `any(itr)`, `all(itr)`.

The associativity of the reduction is implementation dependent. This means that you can't
use non-associative operations like `-` because it is undefined whether `reduce(-,[1,2,3])`
should be evaluated as `(1-2)-3` or `1-(2-3)`. Use [`foldl`](@ref) or
[`foldr`](@ref) instead for guaranteed left or right associativity.

Some operations accumulate error. Parallelism will be easier if the reduction can be
executed in groups. Future versions of Julia might change the algorithm. Note that the
elements are not reordered if you use an ordered collection.

# Examples
```jldoctest
julia> reduce(*, [2; 3; 4])
24

julia> reduce(*, [2; 3; 4]; init=-1)
-24
```
"""
reduce(op, itr; kw...) = mapreduce(identity, op, itr; kw...)

reduce(op, a::Number) = a  # Do we want this?

###### Specific reduction functions ######

## sum

"""
    sum(f, itr)

Sum the results of calling function `f` on each element of `itr`.

The return type is `Int` for signed integers of less than system word size, and
`UInt` for unsigned integers of less than system word size.  For all other
arguments, a common return type is found to which all arguments are promoted.

# Examples
```jldoctest
julia> sum(abs2, [2; 3; 4])
29
```

Note the important difference between `sum(A)` and `reduce(+, A)` for arrays
with small integer eltype:

```jldoctest
julia> sum(Int8[100, 28])
128

julia> reduce(+, Int8[100, 28])
-128
```

In the former case, the integers are widened to system word size and therefore
the result is 128. In the latter case, no such widening happens and integer
overflow results in -128.
"""
sum(f, a) = mapreduce(f, add_sum, a)

"""
    sum(itr)

Returns the sum of all elements in a collection.

The return type is `Int` for signed integers of less than system word size, and
`UInt` for unsigned integers of less than system word size.  For all other
arguments, a common return type is found to which all arguments are promoted.

# Examples
```jldoctest
julia> sum(1:20)
210
```
"""
sum(a) = sum(identity, a)
sum(a::AbstractArray{Bool}) = count(a)

## prod
"""
    prod(f, itr)

Returns the product of `f` applied to each element of `itr`.

The return type is `Int` for signed integers of less than system word size, and
`UInt` for unsigned integers of less than system word size.  For all other
arguments, a common return type is found to which all arguments are promoted.

# Examples
```jldoctest
julia> prod(abs2, [2; 3; 4])
576
```
"""
prod(f, a) = mapreduce(f, mul_prod, a)

"""
    prod(itr)

Returns the product of all elements of a collection.

The return type is `Int` for signed integers of less than system word size, and
`UInt` for unsigned integers of less than system word size.  For all other
arguments, a common return type is found to which all arguments are promoted.

# Examples
```jldoctest
julia> prod(1:20)
2432902008176640000
```
"""
prod(a) = mapreduce(identity, mul_prod, a)

## maximum & minimum
_fast(::typeof(min),x,y) = min(x,y)
_fast(::typeof(max),x,y) = max(x,y)
function _fast(::typeof(max), x::AbstractFloat, y::AbstractFloat)
    ifelse(isnan(x),
        x,
        ifelse(x > y, x, y))
end

function _fast(::typeof(min),x::AbstractFloat, y::AbstractFloat)
    ifelse(isnan(x),
        x,
        ifelse(x < y, x, y))
end

isbadzero(::typeof(max), x::AbstractFloat) = (x == zero(x)) & signbit(x)
isbadzero(::typeof(min), x::AbstractFloat) = (x == zero(x)) & !signbit(x)
isbadzero(op, x) = false
isgoodzero(::typeof(max), x) = isbadzero(min, x)
isgoodzero(::typeof(min), x) = isbadzero(max, x)

function mapreduce_impl(f, op::Union{typeof(max), typeof(min)},
                        A::AbstractArray, first::Int, last::Int)
    a1 = @inbounds A[first]
    v1 = mapreduce_first(f, op, a1)
    v2 = v3 = v4 = v1
    chunk_len = 256
    start = first + 1
    simdstop  = start + chunk_len - 4
    while simdstop <= last - 3
        # short circuit in case of NaN
        v1 == v1 || return v1
        v2 == v2 || return v2
        v3 == v3 || return v3
        v4 == v4 || return v4
        @inbounds for i in start:4:simdstop
            v1 = _fast(op, v1, f(A[i+0]))
            v2 = _fast(op, v2, f(A[i+1]))
            v3 = _fast(op, v3, f(A[i+2]))
            v4 = _fast(op, v4, f(A[i+3]))
        end
        checkbounds(A, simdstop+3)
        start += chunk_len
        simdstop += chunk_len
    end
    v = op(op(v1,v2),op(v3,v4))
    for i in start:last
        @inbounds ai = A[i]
        v = op(v, f(ai))
    end

    # enforce correct order of 0.0 and -0.0
    # e.g. maximum([0.0, -0.0]) === 0.0
    # should hold
    if isbadzero(op, v)
        for i in first:last
            x = @inbounds A[i]
            isgoodzero(op,x) && return x
        end
    end
    return v
end

maximum(f, a) = mapreduce(f, max, a)
minimum(f, a) = mapreduce(f, min, a)

"""
    maximum(itr)

Returns the largest element in a collection.

# Examples
```jldoctest
julia> maximum(-20.5:10)
9.5

julia> maximum([1,2,3])
3
```
"""
maximum(a) = mapreduce(identity, max, a)

"""
    minimum(itr)

Returns the smallest element in a collection.

# Examples
```jldoctest
julia> minimum(-20.5:10)
-20.5

julia> minimum([1,2,3])
1
```
"""
minimum(a) = mapreduce(identity, min, a)

## all & any

"""
    any(itr) -> Bool

Test whether any elements of a boolean collection are `true`, returning `true` as
soon as the first `true` value in `itr` is encountered (short-circuiting).

If the input contains [`missing`](@ref) values, return `missing` if all non-missing
values are `false` (or equivalently, if the input contains no `true` value), following
[three-valued logic](https://en.wikipedia.org/wiki/Three-valued_logic).

# Examples
```jldoctest
julia> a = [true,false,false,true]
4-element Array{Bool,1}:
  true
 false
 false
  true

julia> any(a)
true

julia> any((println(i); v) for (i, v) in enumerate(a))
1
true

julia> any([missing, true])
true

julia> any([false, missing])
missing
```
"""
any(itr) = any(identity, itr)

"""
    all(itr) -> Bool

Test whether all elements of a boolean collection are `true`, returning `false` as
soon as the first `false` value in `itr` is encountered (short-circuiting).

If the input contains [`missing`](@ref) values, return `missing` if all non-missing
values are `true` (or equivalently, if the input contains no `false` value), following
[three-valued logic](https://en.wikipedia.org/wiki/Three-valued_logic).

# Examples
```jldoctest
julia> a = [true,false,false,true]
4-element Array{Bool,1}:
  true
 false
 false
  true

julia> all(a)
false

julia> all((println(i); v) for (i, v) in enumerate(a))
1
2
false

julia> all([missing, false])
false

julia> all([true, missing])
missing
```
"""
all(itr) = all(identity, itr)

"""
    any(p, itr) -> Bool

Determine whether predicate `p` returns `true` for any elements of `itr`, returning
`true` as soon as the first item in `itr` for which `p` returns `true` is encountered
(short-circuiting).

If the input contains [`missing`](@ref) values, return `missing` if all non-missing
values are `false` (or equivalently, if the input contains no `true` value), following
[three-valued logic](https://en.wikipedia.org/wiki/Three-valued_logic).

# Examples
```jldoctest
julia> any(i->(4<=i<=6), [3,5,7])
true

julia> any(i -> (println(i); i > 3), 1:10)
1
2
3
4
true

julia> any(i -> i > 0, [1, missing])
true

julia> any(i -> i > 0, [-1, missing])
missing

julia> any(i -> i > 0, [-1, 0])
false
```
"""
any(f, itr) = _any(f, itr, :)

function _any(f, itr, ::Colon)
    anymissing = false
    for x in itr
        v = f(x)
        if ismissing(v)
            anymissing = true
        elseif v
            return true
        end
    end
    return anymissing ? missing : false
end

"""
    all(p, itr) -> Bool

Determine whether predicate `p` returns `true` for all elements of `itr`, returning
`false` as soon as the first item in `itr` for which `p` returns `false` is encountered
(short-circuiting).

If the input contains [`missing`](@ref) values, return `missing` if all non-missing
values are `true` (or equivalently, if the input contains no `false` value), following
[three-valued logic](https://en.wikipedia.org/wiki/Three-valued_logic).

# Examples
```jldoctest
julia> all(i->(4<=i<=6), [4,5,6])
true

julia> all(i -> (println(i); i < 3), 1:10)
1
2
3
false

julia> all(i -> i > 0, [1, missing])
missing

julia> all(i -> i > 0, [-1, missing])
false

julia> all(i -> i > 0, [1, 2])
true
```
"""
all(f, itr) = _all(f, itr, :)

function _all(f, itr, ::Colon)
    anymissing = false
    for x in itr
        v = f(x)
        if ismissing(v)
            anymissing = true
        # this syntax allows throwing a TypeError for non-Bool, for consistency with any
        elseif v
            continue
        else
            return false
        end
    end
    return anymissing ? missing : true
end

## count

"""
    count(p, itr) -> Integer
    count(itr) -> Integer

Count the number of elements in `itr` for which predicate `p` returns `true`.
If `p` is omitted, counts the number of `true` elements in `itr` (which
should be a collection of boolean values).

# Examples
```jldoctest
julia> count(i->(4<=i<=6), [2,3,4,5,6])
3

julia> count([true, false, true, true])
3
```
"""
function count(pred, itr)
    n = 0
    for x in itr
        n += pred(x)::Bool
    end
    return n
end
function count(pred, a::AbstractArray)
    n = 0
    for i in eachindex(a)
        @inbounds n += pred(a[i])::Bool
    end
    return n
end
count(itr) = count(identity, itr)

max_stars_repo_name: ashleysommer/julia

max_issues_repo_name: ashleysommer/julia

max_forks_repo_name: ashleysommer/julia


text:
using Adapt
using Statistics

import Oceananigans.BoundaryConditions: fill_halo_regions!

#####
##### AbstractReducedField stuff
#####

abstract type AbstractReducedField{X, Y, Z, A, G, T, N} <: AbstractDataField{X, Y, Z, A, G, T, 3} end

const ARF = AbstractReducedField

fill_halo_regions!(field::AbstractReducedField, arch, args...) =
    fill_halo_regions!(field.data, field.boundary_conditions, arch, field.grid, args...; reduced_dimensions=field.dims)

const DimsType = NTuple{N, Int} where N

function validate_reduced_dims(dims)
    dims = Tuple(dims)

    # Check type
    dims isa DimsType || error("Reduced dims must be an integer or tuple of integers.")

    # Check length
    length(dims) > 3  && error("Models are 3-dimensional. Cannot reduce over 4+ dimensions.")

    # Check values
    all(1 <= d <= 3 for d in dims) || error("Dimensions must be one of 1, 2, 3.")

    return dims
end

function validate_reduced_locations(X, Y, Z, dims)
    loc = (X, Y, Z)

    # Check reduced locations
    for i in dims
        isnothing(loc[i]) || ArgumentError("The location of reduced dimensions must be Nothing")
    end

    return nothing
end

#####
##### Concrete ReducedField
#####

struct ReducedField{X, Y, Z, A, D, G, T, N, B} <: AbstractReducedField{X, Y, Z, A, G, T, N}
                   data :: D
           architecture :: A
                   grid :: G
                   dims :: NTuple{N, Int}
    boundary_conditions :: B

    @doc """
        ReducedField{X, Y, Z}(data, grid, dims)

    Returns a `ReducedField` at location `(X, Y, Z)` with `data` on `grid`
    that is reduced over the dimensions in `dims`.
    """
    function ReducedField{X, Y, Z}(data::D, arch::A, grid::G, dims, bcs::B) where {X, Y, Z, A, D, G, B}

        dims = validate_reduced_dims(dims)
        validate_reduced_locations(X, Y, Z, dims)
        validate_field_data(X, Y, Z, data, grid)

        N = length(dims)
        T = eltype(grid)

        return new{X, Y, Z, A, D, G, T, N, B}(data, arch, grid, dims, bcs)
    end
end

"""
    ReducedField(X, Y, Z, arch, grid; dims, data=nothing, boundary_conditions=nothing)

Returns a `ReducedField` reduced over `dims` on `grid` and `arch`itecture with `boundary_conditions`.

The location `(X, Y, Z)` may be the parent, three-dimension location or the reduced location.

If `data` is specified, it should be an `OffsetArray` with singleton reduced dimensions;
otherwise `data` is allocated.

If `boundary_conditions` are not provided, default boundary conditions are constructed
using the reduced location.
"""
function ReducedField(FT::DataType, Xr, Yr, Zr, arch, grid::AbstractGrid; dims, data=nothing,
                      boundary_conditions=nothing)

    dims = validate_reduced_dims(dims)

    # Reduce non-reduced dimensions
    X, Y, Z = reduced_location((Xr, Yr, Zr); dims=dims)

    if isnothing(data)
        data = new_data(FT, arch, grid, (X, Y, Z))
    end

    if isnothing(boundary_conditions)
        boundary_conditions = FieldBoundaryConditions(grid, (X, Y, Z))
    end

    return ReducedField{X, Y, Z}(data, arch, grid, dims, boundary_conditions)
end

ReducedField(Xr, Yr, Zr, arch, grid::AbstractGrid; dims, kw...) = ReducedField(eltype(grid), Xr, Yr, Zr, arch, grid; dims, kw...)
ReducedField(Lr::Tuple, arch, grid::AbstractGrid; dims, kw...) = ReducedField(Lr..., arch, grid; dims, kw...)
ReducedField(FT::DataType, Lr::Tuple, arch, grid::AbstractGrid; dims, kw...) = ReducedField(FT, Lr..., arch, grid; dims, kw...)

# Canonical `similar` for AbstractReducedField
Base.similar(r::AbstractReducedField{X, Y, Z, Arch}) where {X, Y, Z, Arch} =
    ReducedField(X, Y, Z, Arch(), r.grid; dims=r.dims, boundary_conditions=r.boundary_conditions)

#####
##### ReducedField utils
#####

reduced_location(loc; dims) = Tuple(i ∈ dims ? Nothing : loc[i] for i in 1:3)

Adapt.adapt_structure(to, reduced_field::ReducedField{X, Y, Z}) where {X, Y, Z} =
    ReducedField{X, Y, Z}(adapt(to, reduced_field.data), nothing, adapt(to, reduced_field.grid), reduced_field.dims, nothing)

#####
##### Field reductions
#####

# Risky to use these without tests. Docs would also be nice.
Statistics.norm(a::AbstractDataField) = sqrt(mapreduce(x -> x * x, +, interior(a)))
Statistics.dot(a::AbstractDataField, b::AbstractDataField) = mapreduce((x, y) -> x * y, +, interior(a), interior(b))

# The more general case, for AbstractOperations
function Statistics.norm(a::AbstractField)
    arch = architecture(a)
    grid = a.grid

    r = zeros(arch, grid, 1)

    Base.mapreducedim!(x -> x * x, +, r, a)

    return CUDA.@allowscalar sqrt(r[1])
end

# TODO: In-place allocations with function mappings need to be fixed in Julia Base...
const SumReduction     = typeof(Base.sum!)
const ProdReduction    = typeof(Base.prod!)
const MaximumReduction = typeof(Base.maximum!)
const MinimumReduction = typeof(Base.minimum!)
const AllReduction     = typeof(Base.all!)
const AnyReduction     = typeof(Base.any!)

initialize_reduced_field!(::SumReduction,  f, r::AbstractReducedField, c) = Base.initarray!(interior(r), Base.add_sum, true, interior(c))
initialize_reduced_field!(::ProdReduction, f, r::AbstractReducedField, c) = Base.initarray!(interior(r), Base.mul_prod, true, interior(c))
initialize_reduced_field!(::AllReduction,  f, r::AbstractReducedField, c) = Base.initarray!(interior(r), &, true, interior(c))
initialize_reduced_field!(::AnyReduction,  f, r::AbstractReducedField, c) = Base.initarray!(interior(r), |, true, interior(c))

initialize_reduced_field!(::MaximumReduction, f, r::AbstractReducedField, c) = Base.mapfirst!(f, interior(r), c)
initialize_reduced_field!(::MinimumReduction, f, r::AbstractReducedField, c) = Base.mapfirst!(f, interior(r), c)

filltype(f, grid) = eltype(grid)
filltype(::Union{AllReduction, AnyReduction}, grid) = Bool

# Allocating and in-place reductions
for reduction in (:sum, :maximum, :minimum, :all, :any)

    reduction! = Symbol(reduction, '!')

    @eval begin

        # In-place
        Base.$(reduction!)(f::Function, r::AbstractReducedField, a::AbstractArray; kwargs...) =
            Base.$(reduction!)(f, interior(r), a; kwargs...)

        Base.$(reduction!)(r::AbstractReducedField, a::AbstractArray; kwargs...) =
            Base.$(reduction!)(identity, interior(r), a; kwargs...)

        # Allocating
        function Base.$(reduction)(f::Function, c::AbstractField; dims=:)
            if dims === (:)
                r = zeros(architecture(c), c.grid, 1, 1, 1)
                Base.$(reduction!)(f, r, c)
                return CUDA.@allowscalar r[1, 1, 1]
            else
                FT = filltype(Base.$(reduction!), c.grid)
                r = ReducedField(FT, location(c), architecture(c), c.grid; dims)
                initialize_reduced_field!(Base.$(reduction!), f, r, c)
                Base.$(reduction!)(f, r, c, init=false)
                return r
            end
        end

        Base.$(reduction)(c::AbstractField; dims=:) = Base.$(reduction)(identity, c; dims)
    end
end

Statistics._mean(f, c::AbstractField, ::Colon) = sum(f, c) / length(c)

function Statistics._mean(f, c::AbstractField, dims)
    r = sum(f, c; dims)
    n = mapreduce(i -> size(c, i), *, unique(dims); init=1)
    parent(r) ./= n
    return r
end


text:
using UNet: checksize, minsize, trim, outputsize, adjustsize, padding

@testset "utils" begin
    x1 = rand(32, 32, 4, 6)
    x2 = rand(32, 32, 4, 6)
    @test (chcat(x1, x2) |> size) == (32, 32, 8, 6)

    @test checksize((256, 512), padding = true, nlevels = 4) == (true, true)
    @test checksize((572, 572), padding = false, nlevels = 4) == (true, true)
    @test checksize((256, 500), padding = true, nlevels = 4) == (true, false)
    @test checksize((572, 512), padding = false, nlevels = 4) == (true, false)

    @test minsize(true, 4) == 64
    @test minsize(false, 4) == 188

    @test trim(1, 4) == 88

    @test outputsize((512, 512), true, 4) == (512, 512)
    @test outputsize((572, 572), false, 4) == (388, 388)

    @test adjustsize((256, 256), true, 4) == (256, 256)
    @test adjustsize((388, 388), false, 4) == (572, 572)
    @test adjustsize(4, true, 4) == 64
    @test adjustsize(4, false, 4) == 188

    res = (((92, 92), (92, 92)), ((0, 0), (0, 0)))
    @test padding((388, 388), padding = false, nlevels = 4) == res
    res = (((0, 0), (0, 0)), ((0, 0), (0, 0)))
    @test padding((256, 256), padding = true, nlevels = 4) == res
end