text:
using StatisticalRethinking: sr_datadir, PI
using StatisticalRethinkingCommon
using Statistics: mean, median
using StatsBase: sample
using Plots: plot, display, scatter!, plot!
using DataFrames: DataFrame
using Gen: @gen, (~), normal, gamma, get_args, get_retval
import Gen
import CSV

@info "Reading Howell dataset"
@time howell1 = CSV.read(sr_datadir("Howell1.csv"), DataFrame)

howell1_adults = howell1[howell1.age .>= 18, :]
howell1_children = howell1[howell1.age .< 13, :]
howell1_boys = howell1_children[howell1_children.male .== 1, :]
howell1_girls = howell1_children[howell1_children.male .== 0, :]

function plot_howell1()

    height_weight_plot = plot(howell1.height, howell1.weight, seriestype=:scatter, color=:blue, title="Height/Weight", xlabel="height (cm)", ylabel="weight (kg)", labels="data")
    return height_weight_plot

end

function plot_howell1_adults()

    height_weight_plot = plot(howell1_adults.height, howell1_adults.weight, seriestype=:scatter, color=:blue, title="Height/Weight Adults", xlabel="height (cm)", ylabel="weight (kg)", labels="data")
    return height_weight_plot

end

function plot_howell1_children()

    height_weight_plot = plot(howell1_children.height, howell1_children.weight, seriestype=:scatter, color=:blue, title="Height/Weight Children", xlabel="height (cm)", ylabel="weight (kg)", labels="data")
    age_weight_plot = plot(howell1_children.age, howell1_children.weight, seriestype=:scatter, color=:blue, title="Age/Weight Children", xlabel="age (years)", ylabel="weight (kg)", labels="data")
    return plot(height_weight_plot, age_weight_plot, layout=2)

end

function plot_howell1_age_weight_boys_girls()
    age_weight_plot = plot([howell1_boys.age, howell1_girls.age],
                           [howell1_boys.weight, howell1_girls.weight],
                           seriestype=:scatter,
                           color=[:blue :red],
                           title="Age/Weight Children",
                           xlabel="age (years)",
                           ylabel="weight (kg)",
                           labels=["Boys" "Girls"])
    return age_weight_plot
end

function plot_howell1_height_weight_boys_girls()
    height_weight_plot = plot([howell1_boys.height, howell1_girls.height],
                              [howell1_boys.weight, howell1_girls.weight],
                              seriestype=:scatter,
                              color=[:blue :red],
                              title="Height/Weight Children",
                              xlabel="height (cm)",
                              ylabel="weight (kg)",
                              labels=["Boys" "Girls"])
    return height_weight_plot
end

function plot_howell1_boys_girls()
    return plot(plot_howell1_height_weight_boys_girls(), plot_howell1_boys_girls(), layout=2)

end

@gen function pos_line_model(xs,
                             alpha_mu=0.,
                             alpha_sigma=1.,
                             beta_shape=1.,
                             beta_scale=1.,
                             sigma_shape=1.,
                             sigma_scale=1.)

    alpha ~ normal(alpha_mu, alpha_sigma)
    beta ~ gamma(beta_shape, beta_scale)
    sigma ~ gamma(sigma_shape, sigma_scale)

    function f(x)
        return alpha .+ beta .* x
    end

    for (i, x) in enumerate(xs)
        ({(:y, i)} ~ normal(f(x), sigma))
    end

    return f

end

@gen function problem3_model(xs, male)
    alpha0  ~ normal(5., 1.0)
    beta0   ~  gamma(1., 1.5)
    sigma0  ~  gamma(1., 1.0)
    alpha1  ~ normal(5., 1.0)
    beta1   ~  gamma(1., 1.5)
    sigma1  ~  gamma(1., 1.0)

    function f(x, m)
        return (alpha0 * (1 - m) + alpha1 * m) + (beta0  * (1 - m) + beta1 * m) * x
    end

    for (i, x) in enumerate(xs)
        m = male[i]
        {(:y,  i)} ~ normal(f(x,m), sigma0 * (1-m) + sigma1 * m)
    end

    return f
end

function render_linear_trace(trace; show_data=true)

    # Pull out xs from the trace
    xs, = get_args(trace)

    xmin = minimum(xs)
    xmax = maximum(xs)

    # Pull out the return value, useful for plotting
    y = get_retval(trace)

    # Draw the line
    test_xs = collect(range(xmin - 5, stop=xmax + 5, length=1000))
    test_ys = map(y, test_xs)

    ymin = minimum(test_ys)
    ymax = maximum(test_ys)

    fig = plot(test_xs, map(y, test_xs), color="black", alpha=0.5, label=nothing,
                xlim=(xmin-0.5, xmax+0.5), ylim=(ymin-0.5, ymax+0.5))

    if show_data
        ys = [trace[(:y, i)] for i=1:length(xs)]

        # Plot the data set
        scatter!(xs, ys, c="black", label=nothing)
    end

    return fig
end

function do_importance_inference(model, xs, ys, steps, modelargs...)

    # Create a choice map that maps model addresses (:y, i)
    # to observed values ys[i]. We leave :slope and :intercept
    # unconstrained, because we want them to be inferred.
    observations = Gen.choicemap()
    for (i, y) in enumerate(ys)
        observations[(:y, i)] = y
    end

    # Call importance_resampling to obtain a likely trace consistent
    # with our observations.
    (trace, _) = Gen.importance_resampling(model, (xs, modelargs...), observations, steps)
    return trace
end

function overlay(renderer, traces; same_data=true, args...)
    fig = renderer(traces[1], show_data=true, args...)

    xs, = get_args(traces[1])
    xmin = minimum(xs)
    xmax = maximum(xs)

    # Draw the line

    for i=2:length(traces)
        y = get_retval(traces[i])
        test_xs = collect(range(xmin - 5, stop=xmax + 5, length=1000))
        test_ys = map(y, test_xs)
        ymin = minimum(test_ys)
        ymax = maximum(test_ys)
        fig = plot!(test_xs, map(y, test_xs), color="black", alpha=0.5, label=nothing,
                    xlim=(xmin - 0.5 , xmax + 0.5), ylim=(ymin - 0.5, ymax + 0.5))
    end
    return fig
end

function predict_new_data(model, trace, new_xs::Vector{Float64}, param_addrs)

    # Copy parameter values from the inferred trace (`trace`)
    # into a fresh set of constraints.
    constraints = Gen.choicemap()
    for addr in param_addrs
        constraints[addr] = trace[addr]
    end

    # Run the model with new x coordinates, and with parameters
    # fixed to be the inferred values.
    (new_trace, _) = Gen.generate(model, (new_xs,), constraints)

    # Pull out the y-values and return them.
    ys = [new_trace[(:y, i)] for i=1:length(new_xs)]
    return ys
end

function infer_and_predict(model, xs, ys, new_xs, param_addrs, num_traces, steps, modelargs...)
    pred_ys = []
    for _ in 1:num_traces
        trace = do_importance_inference(model, xs, ys, steps, modelargs...)
        push!(pred_ys, predict_new_data(model, trace, new_xs, param_addrs))
    end
    return pred_ys
end

"""
Construct a linear regression of weight as predicted by height, using the
adults (age 18 or greater) from the Howell1 dataset. The heights listed below
were recorded in the !Kung census, but weights were not recorded for these
individuals. Provide predicted weights and 89% compatibility intervals for
each of these individuals. That is, fill in the table below, using model-based
predictions.

"""
function problem1(samples=100)

    alpha_mu = 50.
    alpha_sigma= 10.

    beta_shape=1.
    beta_scale=5.

    sigma_shape=1.
    sigma_scale=5.

    xbar = mean(howell1_adults.height)
    xs = howell1_adults.height .- xbar
    ys = howell1_adults.weight

    @info "Inferring:"
    @time all_predictions = infer_and_predict(pos_line_model, xs, ys, [140, 160, 175] .- xbar, (:alpha, :beta, :sigma), samples, 1000, alpha_mu, alpha_sigma, beta_shape, beta_scale, sigma_shape, sigma_scale)
    all_predictions = vectorlist_to_matrix(all_predictions)
    @info "Calculating intervals:"
    @time intervals = [PI(all_predictions[:,i]) for i in (1,2,3)]
    return mean(all_predictions, dims=1), intervals
end

function problem1_support()

    alpha_mu = 50.
    alpha_sigma= 10.

    beta_shape=1.
    beta_scale=5.

    sigma_shape=1.
    sigma_scale=5.

    fs = [pos_line_model((), alpha_mu, alpha_sigma, beta_shape, beta_scale, sigma_shape, sigma_scale) for _ in 1:15]

    prior_plot = plot(xs, [f(xs) for f in fs], legend=false, title="Priors")

    traces = [do_importance_inference(pos_line_model, xs, ys, 100, alpha_mu, alpha_sigma, beta_shape, beta_scale, sigma_shape, sigma_scale) for _ in 1:10]

    posterior = overlay(render_linear_trace, traces)

    problem1_plot = plot(prior_plot, posterior, layout=2)

    display(problem1_plot)

end

function problem2(steps=100, samples=100)

    alpha_mu = 5.
    alpha_sigma = 1.

    beta_shape = 1.
    beta_scale = 3.

    sigma_shape = 1.
    sigma_scale = 1.

    priors = [pos_line_model((), alpha_mu, alpha_sigma, beta_shape, beta_scale, sigma_shape, sigma_scale) for _ in 1:15]

    xs_test = Vector(range(minimum(howell1_children.age) - 0.5, maximum(howell1_children.age) + 0.5, length=1000))

    prior_plot = plot(xs_test, [f(xs_test) for f in priors], legend=false, title="Priors")

    xs = howell1_children.age
    ys = howell1_children.weight

    traces = [do_importance_inference(pos_line_model, xs, ys, steps, alpha_mu, alpha_sigma, beta_shape, beta_scale, sigma_shape, sigma_scale) for _ in 1:samples]

    posteriors = overlay(render_linear_trace, sample(traces, 15, replace=false))
    plot!(posteriors, title="Posteriors")

    problem2_plot = plot(prior_plot, posteriors, layout=2)
    display(problem2_plot)

end

function problem3(steps=100, samples=100)

    priors_boys = [problem3_model((), ()) for _ in 1:15]
    priors_girls = [problem3_model((), ()) for _ in 1:15]
    xs_test = Vector(range(minimum(howell1_children.age) - 0.5, maximum(howell1_children.age) + 0.5, length=1000))

    prior_plot = plot(xs_test, [f.(xs_test, ones(1000)) for f in priors_boys], color=:blue, labels=nothing)
    plot!(prior_plot, xs_test, [f.(xs_test, zeros(1000)) for f in priors_girls], color=:red, labels=nothing)
    plot!(prior_plot, title="Priors")

    xs = howell1_children.age

    ys = howell1_children.weight

    traces = [do_importance_inference(problem3_model, xs, ys, steps, howell1_children.male) for _ in 1:samples]
    posterior_funs = Gen.get_retval.(traces)

    posterior_plot = plot_howell1_age_weight_boys_girls()
    plot!(posterior_plot, xs_test, [f.(xs_test, ones(1000)) for f in sample(posterior_funs, 7, replace=false)], alpha=1/2, color=:red, legend=false)
    plot!(posterior_plot, xs_test, [f.(xs_test, zeros(1000)) for f in sample(posterior_funs, 8, replace=false)], alpha=1/2, color=:blue, legend=false)
    plot!(posterior_plot, title="Posteriors")

    problem3_plot = plot(prior_plot, posterior_plot, layout=2)
    display(problem3_plot)

end


text:
abstract type AbstractGameModel <: AbstractModel
end

function Base.size(model::AbstractGameModel)
	return model.n, model.m, model.pu, model.p
end

ext: jl


text:
import JSON, Glob

function unliterate(path, output)
    if isdir(path)
        !isdir(output) && mkdir(output)
        for file in Glob.glob(joinpath(path, "*.ipynb"))
            unliterate(file, joinpath(output, splitext(splitdir(file)[2])[1] * ".jl"))
        end
        return
    end

    endswith(path, ".ipynb") || throw(ArgumentError("only Jupyter notebooks are supported"))
    nb = open(JSON.parse, path, "r")
    out = open(output, "w")
    prev_type = ""
    for cell in nb["cells"]
        type = cell["cell_type"]
        type == prev_type && print(out, "#-\n\n")
        prev_type = type

        if type == "markdown"
            for line in cell["source"]
                print(out, "#")
                if !isempty(strip(line))
                    print(out, " ")
                    print(out, line)
                else
                    print(out, "\n")
                end
            end
            print(out, "\n\n")
        elseif type == "code"
            for line in cell["source"]
                startswith(line, "# ") && print(out, "#")
                print(out, line)
            end
            print(out, "\n\n")
        else
            error("unknown cell type $type")
        end
    end
    close(out)
    return
end

function main()
    if length(ARGS) != 2
        println("USAGE: unliterate source.ipynb output.jl")
        exit(1)
    end
    unliterate(ARGS[1], ARGS[2])
    return
end

main()

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
mutable struct PatternCache
    table::Dict{BitArray{1}, Int}
    count::Int
    non_uniform_pattern::BitArray{1}
end

function PatternCache(points::Integer)
    temp_pattern = zeros(Bool, points)
    temp_pattern[1:2:points] .= true
    table = Dict{BitArray{1}, Int}()
    table[temp_pattern] = points * (points - 1) + 2
    PatternCache(table, 0, temp_pattern)
end

function lbp_original(bit_pattern::BitArray{1}, uniform_params::PatternCache)
    sum((b * 1) << (length(bit_pattern) - i) for (i, b) in enumerate(bit_pattern)), uniform_params
end

function lbp_uniform(bit_pattern::BitArray{1}, uniform_params::PatternCache)
    variations = sum(bit_pattern[i] != bit_pattern[i + 1] for i in 1:length(bit_pattern) - 1)
    if variations <= 2
        haskey(uniform_params.table, bit_pattern) && return uniform_params.table[bit_pattern], uniform_params
        uniform_params.count += 1
        uniform_params.table[copy(bit_pattern)] = uniform_params.count
        return uniform_params.count, uniform_params
    else
        return uniform_params.table[uniform_params.non_uniform_pattern], uniform_params
    end
end

function lbp_rotation_invariant(bit_pattern::BitArray{1}, uniform_params::PatternCache)
    mini, _ = lbp_original(bit_pattern, uniform_params)
    for i in 2:length(bit_pattern)
          mini = min(mini, lbp_original(vcat(bit_pattern[i:end], bit_pattern[1:i-1]), uniform_params)[1])
    end
    mini, uniform_params
end

function _lbp(img::AbstractArray{T, 2}, points::Integer, offsets::Array, method::Function = lbp_original) where T<:Gray
    uniform_params = PatternCache(points)
    lbp_image = zeros(UInt, size(img))
    R = CartesianIndices(size(img))
    bit_pattern = falses(length(offsets))
    for I in R
        for (i, o) in enumerate(offsets) bit_pattern[i] = img[I] >= bilinear_interpolation(img, I[1] + o[1], I[2] + o[2]) end
        lbp_image[I], uniform_params = method(bit_pattern, uniform_params)
    end
    lbp_image
end

const original_offsets = [[- 1, - 1], [- 1, 0], [- 1, 1], [0, 1], [1, 1], [1, 0], [1, - 1], [0, - 1]]

function circular_offsets(points::Integer, radius::Number)

    return [(round(- radius * sin(2 * pi * i / points), digits=5), round(radius * cos(2 * pi * i / points), digits=5)) for i = 0:points - 1]
end

lbp(img::AbstractArray{T, 2}, method::Function = lbp_original) where {T<:Gray} = _lbp(img, 8, original_offsets, method)

lbp(img::AbstractArray{T, 2}, points::Integer, radius::Number, method::Function = lbp_original) where {T<:Gray} = _lbp(img, points, circular_offsets(points, radius), method)

function _modified_lbp(img::AbstractArray{T, 2}, points::Integer, offsets::Array, method::Function = lbp_original) where T<:Gray
    uniform_params = PatternCache(points)
    lbp_image = zeros(UInt, size(img))
    R = CartesianIndices(size(img))
    bit_pattern = falses(length(offsets))
    for I in R
        avg = (sum(bilinear_interpolation(img, I[1] + o[1], I[2] + o[2]) for o in offsets) + img[I]) / (points + 1)
        for (i, o) in enumerate(offsets) bit_pattern[i] = avg >= bilinear_interpolation(img, I[1] + o[1], I[2] + o[2]) end
        lbp_image[I], uniform_params = method(bit_pattern, uniform_params)
    end
    lbp_image
end

modified_lbp(img::AbstractArray{T, 2}, method::Function = lbp_original) where {T<:Gray} = _modified_lbp(img, 8, original_offsets, method)

modified_lbp(img::AbstractArray{T, 2}, points::Integer, radius::Number, method::Function = lbp_original) where {T<:Gray} = _modified_lbp(img, points, circular_offsets(points, radius), method)

function _direction_coded_lbp(img::AbstractArray{T, 2}, offsets::Array) where T
    lbp_image = zeros(UInt, size(img))
    R = CartesianIndices(size(img))
    p = Int(length(offsets) / 2)
    raw_img = convert(Array{Int}, rawview(channelview(img)))
    neighbours = zeros(Int, length(offsets))
    for I in R
        for (i, o) in enumerate(offsets)
            neighbours[i] = Int(bilinear_interpolation(img, I[1] + o[1], I[2] + o[2]).val.i)
        end
        lbp_image[I] = sum(((neighbours[j] - raw_img[I]) * (neighbours[j + p] - raw_img[I]) >= 0) * (2 ^ (2 * p - 2 * j + 1)) +
                            (abs(neighbours[j] - raw_img[I]) >= abs(neighbours[j + p] - raw_img[I])) * (2 ^ (2 * p - 2 * j)) for j in 1:p)
    end
    lbp_image
end

direction_coded_lbp(img::AbstractArray{T, 2}, method::Function = lbp_original) where {T<:Gray} = _direction_coded_lbp(img, original_offsets)

function direction_coded_lbp(img::AbstractArray{T, 2}, points::Integer, radius::Number, method::Function = lbp_original) where T<:Gray
    @assert points % 2 == 0 "For Direction Coded LBP, the number of points must be an even number."
    _direction_coded_lbp(img, circular_offsets(points, radius))
end

function multi_block_lbp(img::AbstractArray{T, 2}, tl_y::Integer, tl_x::Integer, height::Integer, width::Integer) where T<:Gray
    int_img = integral_image(img)
    h, w = size(img)

    @assert (tl_y + 3 * height - 1 <= h) && (tl_x + 3 * width -1 <= w) "Rectangle Grid exceeds image dimensions."

    center = [tl_y + height, tl_x + width]
    central_sum = boxdiff(int_img, tl_y + height : tl_y + 2 * height - 1, tl_x + width : tl_x + 2 * width - 1)
    lbp_code = 0

    for (i, o) in enumerate(original_offsets)
        cur_tl_y = center[1] + o[1] * height
        cur_tl_x = center[2] + o[2] * width
        cur_window_sum = boxdiff(int_img, cur_tl_y : cur_tl_y + height - 1, cur_tl_x : cur_tl_x + height - 1)
        lbp_code += (cur_window_sum > central_sum ? 1 : 0) * 2 ^ (8 - i)
    end
    lbp_code
end

function _create_descriptor(img::AbstractArray{Gray{T}, 2}, yblocks::Integer = 4, xblocks = 4, lbp_type::Function = lbp, args...) where T<:Normed
    h, w = size(img)
    blockh = ceil(Int, h / (yblocks))
    blockw = ceil(Int, w / (xblocks))
    el_max = typemax(FixedPointNumbers.rawtype(eltype(img[1])))
    edges = 0:Int((el_max+1)^0.5):el_max+1
    descriptor = Int[]
    for i in 1:xblocks
        for j in 1:yblocks
            lbp_image = lbp_type(img[(j-1)*blockh+1 : j*blockh, (i-1)*blockw+1 : i*blockw], args...)
            lbp_norm = lbp_image
            _, hist = imhist(lbp_image, edges)
            append!(descriptor, hist[2 : end - 1])
        end
    end
    descriptor
end

function create_descriptor(img::AbstractArray{Gray{T}, 2}, yblocks::Integer = 4, xblocks = 4; lbp_type::Function = lbp, args...) where T<:Normed
    h, w = size(img)
    y_padded = ceil(Int, h / (yblocks)) * yblocks
    x_padded = ceil(Int, w / (xblocks)) * xblocks

    img_padded = Images.imresize(img, (y_padded, x_padded))
    _create_descriptor(img_padded, yblocks, xblocks, lbp_type, args...)
end


text:
export connect_remote, run_server

struct RemoteServer
    host::String
    port::Int
    conn::TCPSocket
end

const CONNS_CTR = Ref{Int}(0)
const CONNS = Dict{Int,RemoteServer}()

function connect_remote(host, port)
    conn = connect(host, port)
    return RemoteServer(string(host), port, conn)
end

run_server(port=rand(1024:65535)) = run_server(Sockets.localhost, port)
function run_server(host, port)
    srv = listen(host, port)
    while true
        conn = accept(srv)
        @async _run_server(conn)
    end
end
function _run_server(conn)
    while true
        try
            blob = deserialize(conn)
            cmd, data = blob.cmd, blob.data
            handle_cmd(conn, Val(Symbol(cmd)), data)
        catch err
            showerror(stderr, err)
            Base.show_backtrace(stderr, catch_backtrace())
            # TODO: Write stacktrace
            serialize(conn, (err=err))
            try
                close(conn)
            finally
                break
            end
        end
    end
end

function handle_cmd(conn, ::Val{:eval}, ex)
    robj = remote(ex)
    serialize(conn, (result=(uuid=robj.uuid,rtype=robj.rtype),))
end
function handle_cmd(conn, ::Val{:get}, uuid)
    robj = _unwrap(uuid)
    serialize(conn, (result=robj,))
end
function handle_cmd(conn, ::Val{:show}, uuid)
    iob = IOBuffer()
    robj = _unwrap(uuid)
    Base.show(iob, robj)
    serialize(conn, (result=String(take!(iob)),))
end

macro remote(conn, ex)
    quote
        serialize($(esc(conn)).conn, (cmd=:eval, data=$ex))
        blob = deserialize($(esc(conn)).conn)
        if haskey(blob, :result)
            T = eval(Symbol(blob.result.rtype))
            uuid = blob.result.uuid
            id = CONNS_CTR[]+1
            CONNS[id] = $(esc(conn))
            RemoteObject(T, uuid, -id)
        else
            err = blob.err
            # FIXME: Wrap with own exception type
            throw(err)
        end
    end
end


text:
"""
    Arrow.Table(io::IO)
    Arrow.Table(file::String)
    Arrow.Table(bytes::Vector{UInt8}, pos=1, len=nothing)

Read an arrow formatted table, from:
 * `io`, bytes come from `read(io)`
 * `file`, bytes come from `Mmap.mmap(file)`
 # `bytes`, a byte vector directly, optionally allowing specifying the starting byte position `pos` and `len`

Returns a `Arrow.Table` object that allows column access via `table.col1`, `table[:col1]`, or `table[1]`.

`Arrow.Table` also satisfies the Tables.jl interface, and so can easily be materialied via any supporting
sink function: e.g. `DataFrame(Arrow.Table(file))`, `SQLite.load!(db, "table", Arrow.Table(file))`, etc.
"""
struct Table <: Tables.AbstractColumns
    names::Vector{Symbol}
    types::Vector{Type}
    columns::Vector{AbstractVector}
    lookup::Dict{Symbol, AbstractVector}
    schema::Ref{Any}
end

Table() = Table(Symbol[], Type[], AbstractVector[], Dict{Symbol, AbstractVector}(), Ref{Meta.Schema}())

names(t::Table) = getfield(t, :names)
types(t::Table) = getfield(t, :types)
columns(t::Table) = getfield(t, :columns)
lookup(t::Table) = getfield(t, :lookup)
schema(t::Table) = getfield(t, :schema)

Tables.istable(::Table) = true
Tables.columnaccess(::Table) = true
Tables.columns(t::Table) = t
Tables.schema(t::Table) = Tables.Schema(names(t), types(t))
Tables.columnnames(t::Table) = names(t)
Tables.getcolumn(t::Table, i::Int) = columns(t)[i]
Tables.getcolumn(t::Table, nm::Symbol) = lookup(t)[nm]

# high-level user API functions
Table(io::IO, pos::Integer=1, len=nothing; debug::Bool=false, convert::Bool=true) = Table(Base.read(io), pos, len; debug=debug, convert=convert)
Table(str::String, pos::Integer=1, len=nothing; debug::Bool=false, convert::Bool=true) = Table(Mmap.mmap(str), pos, len; debug=debug, convert=convert)

# will detect whether we're reading a Table from a file or stream
function Table(bytes::Vector{UInt8}, off::Integer=1, tlen::Union{Integer, Nothing}=nothing; debug::Bool=false, convert::Bool=true)
    len = something(tlen, length(bytes))
    if len > 24 &&
        _startswith(bytes, off, FILE_FORMAT_MAGIC_BYTES) &&
        _endswith(bytes, off + len - 1, FILE_FORMAT_MAGIC_BYTES)
        off += 8 # skip past magic bytes + padding
    end
    t = Table()
    sch = nothing
    dictencodings = Dict{Int64, DictEncoding}()
    dictencoded = Dict{Int64, Tuple{Bool, Type, Meta.Field}}()
    fieldmetadata = Dict{Int, Dict{String, String}}()
    for batch in BatchIterator{debug}(bytes, off)
        # store custom_metadata of batch.msg?
        header = batch.msg.header
        if header isa Meta.Schema
            debug && println("parsing schema message")
            # assert endianness?
            # store custom_metadata?
            for (i, field) in enumerate(header.fields)
                push!(names(t), Symbol(field.name))
                T, metadata = juliaeltype(field)
                if metadata !== nothing
                    fieldmetadata[i] = metadata
                end
                push!(types(t), T)
                d = field.dictionary
                isencoded = false
                if d !== nothing
                    dictencoded[d.id] = (d.isOrdered, types(t)[end], field)
                    isencoded = true
                end
                debug && println("parsed column from schema: name=$(names(t)[end]), type=$(types(t)[end])$(isencoded ? " dictencoded" : "")")
            end
            sch = header
            schema(t)[] = sch
        elseif header isa Meta.DictionaryBatch
            debug && println("parsing dictionary batch message")
            id = header.id
            recordbatch = header.data
            if haskey(dictencodings, id) && header.isDelta
                # delta
                isOrdered, T, field = dictencoded[id]
                values, _, _ = build(T, field, batch, recordbatch, 1, 1, debug)
                dictencoding = dictencodings[id]
                append!(dictencoding.data, values)
                continue
            end
            # new dictencoding or replace
            isOrdered, T, field = dictencoded[id]
            values, _, _ = build(T, field, batch, recordbatch, 1, 1, debug)
            dictencodings[id] = DictEncoding(id, ChainedVector([values]), isOrdered)
            debug && println("parsed parsed dictionary batch message: id=$id, data=$values\n")
        elseif header isa Meta.RecordBatch
            debug && println("parsing record batch message")
            if isempty(columns(t))
                # first RecordBatch
                for vec in VectorIterator{debug}(types(t), sch, batch, dictencodings)
                    push!(columns(t), vec)
                end
                debug && println("parsed 1st record batch")
            elseif !(columns(t)[1] isa ChainedVector)
                # second RecordBatch
                for (i, vec) in enumerate(VectorIterator{debug}(types(t), sch, batch, dictencodings))
                    columns(t)[i] = ChainedVector([columns(t)[i], vec])
                end
                debug && println("parsed 2nd record batch")
            else
                # 2+ RecordBatch
                for (i, vec) in enumerate(VectorIterator{debug}(types(t), sch, batch, dictencodings))
                    append!(columns(t)[i], vec)
                end
                debug && println("parsed additional record batch")
            end
        else
            throw(ArgumentError("unsupported arrow message type: $(typeof(header))"))
        end
    end
    lu = lookup(t)
    for (i, (k, T, col)) in enumerate(zip(names(t), types(t), columns(t)))
        if haskey(fieldmetadata, i) && haskey(fieldmetadata[i], "ARROW:extension:name")
            if fieldmetadata[i]["ARROW:extension:name"] == "JuliaLang.Symbol"
                TT = finaljuliatype(Symbol)
                types(t)[i] = TT
                col = converter(TT, col)
                columns(t)[i] = col
            elseif fieldmetadata[i]["ARROW:extension:name"] == "JuliaLang.Char"
                TT = finaljuliatype(Char)
                types(t)[i] = TT
                col = converter(TT, col)
                columns(t)[i] = col
            end
        end
        if convert
            TT = finaljuliatype(T)
            if TT !== T
                types(t)[i] = TT
                col = converter(TT, col)
                columns(t)[i] = col
            end
        end
        lu[k] = col
    end
    return t
end

struct BatchIterator{debug}
    bytes::Vector{UInt8}
    startpos::Int
end

struct Batch
    msg::Meta.Message
    bytes::Vector{UInt8}
    pos::Int
end

function Base.iterate(x::BatchIterator{debug}, pos=x.startpos) where {debug}
    if pos + 3 > length(x.bytes)
        debug && println("not enough bytes left for another batch message")
        return nothing
    end
    if readbuffer(x.bytes, pos, UInt32) != CONTINUATION_INDICATOR_BYTES
        debug && println("didn't find continuation byte to keep parsing messages: $(readbuffer(x.bytes, pos, UInt32))")
        return nothing
    end
    pos += 4
    if pos + 3 > length(x.bytes)
        debug && println("not enough bytes left to read length of another batch message")
        return nothing
    end
    msglen = readbuffer(x.bytes, pos, Int32)
    if msglen == 0
        debug && println("message has 0 length; terminating message parsing")
        return nothing
    end
    pos += 4
    msg = FlatBuffers.getrootas(Meta.Message, x.bytes, pos-1)
    pos += msglen
    # pos now points to message body
    return Batch(msg, x.bytes, pos), pos + msg.bodyLength
end

struct VectorIterator{debug}
    types::Vector{Type}
    schema::Meta.Schema
    batch::Batch # batch.msg.header MUST BE RecordBatch
    dictencodings::Dict{Int64, DictEncoding}
end

function Base.iterate(x::VectorIterator{debug}, (columnidx, nodeidx, bufferidx)=(1, 1, 1)) where {debug}
    columnidx > length(x.schema.fields) && return nothing
    field = x.schema.fields[columnidx]
    d = field.dictionary
    if d !== nothing
        validity = buildbitmap(x.batch, x.batch.msg.header, bufferidx, debug)
        bufferidx += 1
        buffer = x.batch.msg.header.buffers[bufferidx]
        debug && @show columnidx, buffer.offset, buffer.length
        S = d.indexType === nothing ? Int32 : juliaeltype(field, d.indexType)
        ptr = convert(Ptr{S}, pointer(x.batch.bytes, x.batch.pos + buffer.offset))
        indices = unsafe_wrap(Array, ptr, div(buffer.length, sizeof(S)))
        encoding = x.dictencodings[d.id]
        A = DictEncoded{x.types[columnidx], eltype(indices)}(validity, indices, encoding)
        nodeidx += 1
        bufferidx += 1
    else
        debug && println("parsing column=$columnidx, T=$(x.types[columnidx]), len=$(x.batch.msg.header.nodes[nodeidx].length)")
        A, nodeidx, bufferidx = build(x.types[columnidx], field, x.batch, x.batch.msg.header, nodeidx, bufferidx, debug)
    end
    return A, (columnidx + 1, nodeidx, bufferidx)
end

const ListTypes = Union{Meta.Utf8, Meta.LargeUtf8, Meta.Binary, Meta.LargeBinary, Meta.List, Meta.LargeList}
const LargeLists = Union{Meta.LargeUtf8, Meta.LargeBinary, Meta.LargeList}

build(T, f::Meta.Field, batch, rb, nodeidx, bufferidx, debug) =
    build(T, f, f.type, batch, rb, nodeidx, bufferidx, debug)

function buildbitmap(batch, rb, bufferidx, debug)
    buffer = rb.buffers[bufferidx]
    debug && @show :validity_bitmap, buffer.offset, buffer.length
    voff = batch.pos + buffer.offset
    return ValidityBitmap(batch.bytes, voff, buffer.length)
end

function build(T, f::Meta.Field, L::ListTypes, batch, rb, nodeidx, bufferidx, debug)
    validity = buildbitmap(batch, rb, bufferidx, debug)
    bufferidx += 1
    buffer = rb.buffers[bufferidx]
    debug && @show T, nodeidx, bufferidx, buffer.offset, buffer.length
    ooff = batch.pos + buffer.offset
    OT = L isa LargeLists ? Int64 : Int32
    ptr = convert(Ptr{OT}, pointer(batch.bytes, ooff))
    offsets = Offsets(unsafe_wrap(Array, ptr, div(buffer.length, sizeof(OT))))
    bufferidx += 1
    len = rb.nodes[nodeidx].length
    nodeidx += 1
    if L isa Meta.Utf8 || L isa Meta.LargeUtf8 || L isa Meta.Binary || L isa Meta.LargeBinary
        buffer = rb.buffers[bufferidx]
        debug && @show T, nodeidx, bufferidx, buffer.offset, buffer.length
        bytesptr = pointer(batch.bytes, batch.pos + buffer.offset)
        A = unsafe_wrap(Array, bytesptr, buffer.length)
        bufferidx += 1
    else
        A, nodeidx, bufferidx = build(eltype(Base.nonmissingtype(T)), f.children[1], batch, rb, nodeidx, bufferidx, debug)
    end
    return List{T, OT, typeof(A)}(validity, offsets, A, len), nodeidx, bufferidx
end

function build(T, f::Meta.Field, L::Union{Meta.FixedSizeBinary, Meta.FixedSizeList}, batch, rb, nodeidx, bufferidx, debug)
    validity = buildbitmap(batch, rb, bufferidx, debug)
    bufferidx += 1
    len = rb.nodes[nodeidx].length
    nodeidx += 1
    if L isa Meta.FixedSizeBinary
        buffer = rb.buffers[bufferidx]
        debug && @show T, nodeidx, bufferidx, buffer.offset, buffer.length
        bytesptr = pointer(batch.bytes, batch.pos + buffer.offset)
        A = unsafe_wrap(Array, bytesptr, buffer.length)
        bufferidx += 1
    else
        A, nodeidx, bufferidx = build(eltype(Base.nonmissingtype(T)), f.children[1], batch, rb, nodeidx, bufferidx, debug)
    end
    return FixedSizeList{T, typeof(A)}(validity, A, len), nodeidx, bufferidx
end

function build(S, f::Meta.Field, L::Meta.Map, batch, rb, nodeidx, bufferidx, debug)
    T = Base.nonmissingtype(S)
    # Map is an alias for List<Struct<K, V>>, but the validity/offset buffers for
    # the List/Struct are reduntant/irrelevant, so we skip them
    # skip the entries list & struct nodes
    nodeidx += 2
    # skip the entries list validity + offsets & struct validity
    bufferidx += 3
    A, nodeidx, bufferidx = build(_keytype(T), f.children[1].children[1], batch, rb, nodeidx, bufferidx, debug)
    B, nodeidx, bufferidx = build(_valtype(T), f.children[1].children[2], batch, rb, nodeidx, bufferidx, debug)
    return Map{_keytype(T), _valtype(T), typeof(A), typeof(B)}(A, B), nodeidx, bufferidx
end

function build(T, f::Meta.Field, L::Meta.Struct, batch, rb, nodeidx, bufferidx, debug)
    validity = buildbitmap(batch, rb, bufferidx, debug)
    bufferidx += 1
    len = rb.nodes[nodeidx].length
    NT = Base.nonmissingtype(T)
    N = getn(NT)
    vecs = []
    nodeidx += 1
    for i = 1:N
        A, nodeidx, bufferidx = build(fieldtype(NT, i), f.children[i], batch, rb, nodeidx, bufferidx, debug)
        push!(vecs, A)
    end
    data = Tuple(vecs)
    return Struct{T, typeof(data)}(validity, data, len), nodeidx, bufferidx
end

function build(T, f::Meta.Field, L::Meta.Union, batch, rb, nodeidx, bufferidx, debug)
    buffer = rb.buffers[bufferidx]
    debug && @show T, nodeidx, bufferidx, buffer.offset, buffer.length
    typeidsptr = pointer(batch.bytes, batch.pos + buffer.offset)
    typeIds = unsafe_wrap(Array, typeidsptr, buffer.length)
    debug && @show typeIds
    bufferidx += 1
    if L.mode == Meta.UnionMode.Dense
        buffer = rb.buffers[bufferidx]
        debug && @show T, nodeidx, bufferidx, buffer.offset, buffer.length
        ooff = batch.pos + buffer.offset
        ptr = convert(Ptr{Int32}, pointer(batch.bytes, ooff))
        offsets = unsafe_wrap(Array, ptr, div(buffer.length, 4))
        bufferidx += 1
    end
    vecs = []
    nodeidx += 1
    types = eltype(T)
    for i = 1:fieldcount(types)
        A, nodeidx, bufferidx = build(fieldtype(types, i), f.children[i], batch, rb, nodeidx, bufferidx, debug)
        push!(vecs, A)
    end
    data = Tuple(vecs)
    if L.mode == Meta.UnionMode.Dense
        return DenseUnion{T, typeof(data)}(typeIds, offsets, data), nodeidx, bufferidx
    else
        return SparseUnion{T, typeof(data)}(typeIds, data), nodeidx, bufferidx
    end
end

function build(T, f::Meta.Field, L::Meta.Null, batch, rb, nodeidx, bufferidx, debug)
    return MissingVector(rb.nodes[nodeidx].length), nodeidx + 1, bufferidx
end

# primitives
function build(T, f::Meta.Field, ::L, batch, rb, nodeidx, bufferidx, debug) where {L}
    validity = buildbitmap(batch, rb, bufferidx, debug)
    bufferidx += 1
    buffer = rb.buffers[bufferidx]
    debug && @show T, nodeidx, bufferidx, buffer.offset, buffer.length
    S = Base.nonmissingtype(T)
    ptr = convert(Ptr{S}, pointer(batch.bytes, batch.pos + buffer.offset))
    A = unsafe_wrap(Array, ptr, div(buffer.length, sizeof(S)))
    len = rb.nodes[nodeidx].length
    return Primitive{T, eltype(A)}(validity, A, len), nodeidx + 1, bufferidx + 1
end


text:
deal_new(deck) = reverse(deck)
cut(deck, i) = (i < 0) ? [deck[end+i+1:end]; deck[1:end+i]] : [deck[i+1:end]; deck[1:i]]
function deal_inc(deck, increment)
    result = copy(deck)
    for i in 1:length(deck)
        result[((i-1)*increment % length(deck)) + 1] = deck[i]
    end
    return result
end

"""Parse an operation, returning a function that implements the required technique"""
function parse_operation(str)
    if str == "deal into new stack"
        return deal_new
    elseif startswith(str, "cut")
        return (x -> cut(x, parse(Int32, split(str)[end])))
    elseif startswith(str, "deal with increment")
        return (x -> deal_inc(x, parse(Int32, split(str)[end])))
    end
end

"""Parse an operation, returning (α, β) for the equation of the form αk + β mod N"""
function parse_operation_αβ(str)
    if str == "deal into new stack"
        return (-1, -1)
    elseif startswith(str, "cut")
        return (1, -parse(Int32, split(str)[end]))
    elseif startswith(str, "deal with increment")
        return (parse(Int32, split(str)[end]), 0)
    end
end

"""Compute repeated product Π(αk + β) m times, via power of two repeated multiplication"""
function fast_multiply_mod(α, β, m, N)
    result = (1, 0)
    power_two = (α, β)
    while m > 0
        if m & 1 == 1
            result = (mod(power_two[1] * result[1], N), mod(power_two[1] * result[2] + power_two[2], N))
        end
        m >>= 1
        power_two = (mod(power_two[1] * power_two[1], N), mod(power_two[1]*power_two[2] + power_two[2], N))
    end
    return result
end

function main()
    # Part 1: Apply all of the operations as described in the input file
    operations = parse_operation.(open(f -> readlines(f), "day22_input.txt"))
    deck = foldl((x, op) -> op(x), operations, init=collect(0:10006))
    println("Part 1: Card 2019 is at position $(indexin(2019, deck)[1] - 1)")

    # Part 2: massive number of cards, shuffled a massive number of times
    N = Int128(119315717514047)
    shuffle_repeats = Int128(101741582076661)
    # We need to figure out which card ends up in position 2020
    # What do the operations do to card k?
    # Deal into new stack:
    #   Forward:  D(C_k) = C_(N-1-k) for k ∈ 0:N-1
    # or alternatively
    #   Forward:  D(C_k) = C_(-k-1 mod N)
    # Cut:
    #   Forward:  U(C_k, c) = C_((k-c) mod N)
    # Deal with increment:
    #   Forward:  I(C_k, i) = C_(k*i mod N)
    #
    operations = parse_operation_αβ.(open(f -> readlines(f), "day22_input.txt"))
    α, β = foldl((x, op) -> (mod(x[1]*op[1], N), mod(x[2]*op[1] + op[2], N)), operations, init=(1, 0))
    # The shuffle repeats complicate things! After one round, k -> (αk + β) mod N.
    # After two rounds we have k -> (α^2k + αβ + β) mod N, etc
    # We need to compute the bit pattern of shuffle_repeats and apply it via power of two exponentiation
    α, β = fast_multiply_mod(α, β, shuffle_repeats, N)
    # Now we need to find the number on the card in position 2020. We have the sum:
    # (αk + β) mod N = 2020
    # To solve for k, we can say (αk) mod N = (2020 - β) mod N
    # So we need to solve linear Diophantine equation αk + Ny = (2020 - β) mod N.
    # Use extended Euclidean algorithm which solves related equation αu + Nv = gcd(α, N).
    d, u, v = gcdx(α, N)
    final = 2020
    k = mod(u * mod(final - β, N) ÷ d, N)
    println("Part 2: Card in pos $final has number $k on.")
end
main()

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
@recipe function f(df::DataFrame,
                   titlefont  = font(8, "Courier"),
                   xtickfont  = font(3, "Courier"),
                   ytickfont  = font(3, "Courier"),
                   yguidefont = font(5, "Courier"),
                   xguidefont = font(5, "Courier"),
                   legendfont = font(5, "Courier"),
                   markershape = :circle,
                   markersize = 2,
                   markerstrokewidth = 0,
                   legend = nothing)

    # Default values
    titlefont  := titlefont
    xtickfont  := xtickfont
    ytickfont  := ytickfont
    yguidefont := yguidefont
    xguidefont := xguidefont
    legendfont := legendfont
    markershape := markershape
    markersize := markersize
    markerstrokewidth := markerstrokewidth

    dfts = tots(df)
    names_x = names(dfts)[2:end]
    x = Array(dfts[:,2:end])
    ts = dfts[:,1]
    name_ts = names(dfts)[1]

    legend := isnothing(legend) ? (size(x,2) <= 4 ? :outertop : :outerright) : legend

    @series begin
        label := permutedims(names_x)
        color_palette := :tab10
        seriescolor := reshape(collect(1:size(x,2)),1,size(x,2))
        seriestype := :line

        ts,x
    end
end


text:
module ReadOnlyArrays

export ReadOnlyArray

"""
    ReadOnlyArray(a)

Return a read-only view into the parent array `a`. Most iteration, indexing,
abstract arrays and strided arrays interface functions defined for array `a`
are transparently defined for `ReadOnlyArray(a)`.
The exceptions are `setindex!` which is not allowed for `ReadOnlyArray`
and `similar` which uses the default definitions for `AbstractArray` from Base.
Also when used in broadcasting `ReadOnlyArray` uses default broadcast machinery.

Use `parent` function to access the parent array of `ReadOnlyArray`.

# Examples
```jldoctest
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> r = ReadOnlyArray(a)
2×2 ReadOnlyArray{Int64,2,Array{Int64,2}}:
 1  2
 3  4

julia> r[1]
1

julia> r[1] = 10
ERROR: setindex! not defined for ReadOnlyArray{Int64,2,Array{Int64,2}}
[...]
```
"""
struct ReadOnlyArray{T,N,P} <: AbstractArray{T,N}
    parent::P
    ReadOnlyArray(parent::AbstractArray{T,N}) where{T,N} =
        new{T, N, typeof(parent)}(parent)
end

Base.IteratorSize(::Type{<:ReadOnlyArray{T,N,P}}) where {T,N,P} =
    Base.IteratorSize(P)
Base.IteratorEltype(::Type{<:ReadOnlyArray{T,N,P}}) where {T,N,P} =
    Base.IteratorEltype(P)
Base.eltype(::Type{<:ReadOnlyArray{T,N,P}}) where {T,N,P} =
    eltype(P)
Base.size(roa::ReadOnlyArray, args...) = size(roa.parent, args...)
Base.@propagate_inbounds Base.getindex(roa::ReadOnlyArray, I...) =
    getindex(roa.parent, I...)
Base.firstindex(roa::ReadOnlyArray) = firstindex(roa.parent)
Base.lastindex(roa::ReadOnlyArray) = lastindex(roa.parent)
Base.IndexStyle(::Type{<:ReadOnlyArray{T,N,P}}) where {T,N,P} = IndexStyle(P)
Base.iterate(roa::ReadOnlyArray, args...) = iterate(roa.parent, args...)
Base.length(roa::ReadOnlyArray) = length(roa.parent)

Base.axes(roa::ReadOnlyArray) = axes(roa.parent)
Base.strides(roa::ReadOnlyArray) = strides(roa.parent)
Base.unsafe_convert(p::Type{Ptr{T}}, roa::ReadOnlyArray) where {T} =
    Base.unsafe_convert(p, roa.parent)
Base.stride(roa::ReadOnlyArray, i::Int) = stride(roa.parent, i)
Base.parent(roa::ReadOnlyArray) = roa.parent

end # module


text:
using Giac
using Test

using SpecialFunctions

#For elementary functions we check whether the Giac/Function of a
#numerical value computes (approximately) the same result as the
#corresponding Julia function

for (ex, res) in [
    (sqrt(giac(0.9)), sqrt(0.9)),
    (exp(giac(0.9)), exp(0.9)),
    (log(giac(0.9)), log(0.9)),
    (log(3,giac(0.9)), log(3,0.9)),
    (log10(giac(0.9)), log10(0.9)), # log10 not available in Giac
    (sin(giac(0.9)), sin(0.9)),
    (cos(giac(0.9)), cos(0.9)),
    (tan(giac(0.9)), tan(0.9)),
    (cot(giac(0.9)), cot(0.9)),
    (sec(giac(0.9)), sec(0.9)),
    (csc(giac(0.9)), csc(0.9)),
    (sinh(giac(0.9)), sinh(0.9)),
    (cosh(giac(0.9)), cosh(0.9)),
    (tanh(giac(0.9)), tanh(0.9)),
    (asin(giac(0.9)), asin(0.9)),
    (acos(giac(0.9)), acos(0.9)),
    (atan(giac(0.9)), atan(0.9)),
    (acot(giac(0.9)), acot(0.9)),
    (asec(giac(1.9)), asec(1.9)), # |arg|>=1
    (acsc(giac(1.9)), acsc(1.9)), # |arg|>=1
    (asinh(giac(0.9)), asinh(0.9)),
    (acosh(giac(1.9)), acosh(1.9)), #arg>=1
    (atanh(giac(0.9)), atanh(0.9)),
    (expm1(giac(0.9)), expm1(0.9)),
    (log1p(giac(0.9)), log1p(0.9)),

    (erf(giac(0.9)), erf(0.9)),
    (erfc(giac(0.9)), erfc(0.9)),
    (gamma(giac(0.9)), gamma(0.9)),
    (beta(giac(0.9),giac(1.9)), beta(0.9,1.9)),
    (zeta(giac(0.9)), zeta(0.9)),
    (factorial(giac(13)), factorial(13)),
    (binomial(giac(13),giac(7)), binomial(13,7)),
    (airyai(giac(0.9)), airyai(0.9)), # Giac/Airy_Ai does not work as expected
    (airybi(giac(0.9)), airybi(0.9)), # Giac/Airy_Bi does not work as expected
    (besselj(2,giac(0.9)), besselj(2,0.9)),
    (bessely(2,giac(0.9)), bessely(2,0.9)),
    (digamma(giac(0.9)), digamma(0.9)),
    (polygamma(2,giac(0.9)), polygamma(2,0.9)),
    (sinint(giac(0.9)), sinint(0.9)),
    (cosint(giac(0.9)), cosint(0.9)),
    (Ei(giac(0.9)), 1.6228117136968674413),
]
    @test convert(Float64,ex) ≈ res
end

@test convert(Float64, sign(giac(-12.3))) == -1.0
@test convert(Float64, sign(giac(0.0))) == 0.0
@test convert(Float64, sign(giac(1.3))) == 1.0
@test convert(Float64, round(giac(-1.7))) == -2.0
@test convert(Float64, round(giac(1.4999))) == 1.0
@test convert(Float64, ceil(giac(1.1))) == 2.0
@test convert(Float64, ceil(giac(-1.1))) == -1.0
@test convert(Float64, floor(giac(1.1))) == 1.0
@test convert(Float64, floor(giac(-1.1))) == -2.0
@test convert(Float64, trunc(giac(1.1))) == 1.0
@test convert(Float64, trunc(giac(-1.1))) == -1.0

#Note: We don't check whether the underlying Giac C++ library computes
#correct results. We only test if the Julia functions in Giac.jl/src/library.jl
#are correctly implemented. We assume that this is the case, if the Julia
#function computes the same (or at least an equivalent) result as the
#corresponding Giac function called via the Giac parser.

@giac x y z n

ex1 = simplify(4*atan(giac(1//5))-atan(giac(1//239)))
ex2 = giac("simplify(4*atan(1/5)-atan(1/239))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = collect(x+2*x+1-4)
ex2 = giac("collect(x+2*x+1-4)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = expand((x+y)*(z+1))
ex2 = giac("expand((x+y)*(z+1))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = factor(x^4-4)
ex2 = giac("factor(x^4-4)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = partfrac((x^2-2*x+3)/(x^2-3*x+2))
ex2 = giac("partfrac((x^2-2*x+3)/(x^2-3*x+2))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = numerator((x^3-1)/(x^2-1))
ex2 = giac("numer((x^3-1)/(x^2-1))")
#Note: Giac/numer -> Julia/numerator
@test evaluate(ex1) == evaluate(ex2)

ex1 = denominator((x^3-1)/(x^2-1))
ex2 = giac("denom((x^3-1)/(x^2-1))")
#Note: Giac/denom -> Julia/denominator
@test evaluate(ex1) == evaluate(ex2)

ex1 = left(x^2-1 ⩦ 2*x+3)
ex2 = giac("left(x^2-1 = 2*x+3)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = right(x^2-1 ⩦ 2*x+3)
ex2 = giac("right(x^2-1 = 2*x+3)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = diff(x^3-x,x)
ex2 = giac("diff(x^3-x,x)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = diff(exp(x*y),x,x,x,y,y)
ex2 = giac("diff(exp(x*y),x,x,x,y,y)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = integrate(1/(4+z^2),z)
ex2 = giac("integrate(1/(4+z^2),z)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = integrate(1/(1-x^4),x,2,3)
ex2 = giac("integrate(1/(1-x^4),x,2,3)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = limit(sin(x)/(x^2-3*x),x=>0)
ex2 = giac("limit(sin(x)/(x^2-3*x),x=0)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = limit((2*x-1)/exp(1/(x-1)),x=>plus_inf)
ex2 = giac("limit((2*x-1)/exp(1/(x-1)),x=+infinity)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = limit(sign(x),x=>0,1)
ex2 = giac("limit(sign(x),x=0,1)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = series((x^4+x+2)/(x^2+1),x=>0, 5)
ex2 = giac("series((x^4+x+2)/(x^2+1),x=0, 5)")
@test evaluate(ex1) == evaluate(ex2)

#ex1 = series(1/(1+x+y),[x,y],[0,0],5)
#ex2 = giac("series(1/(1+x+y),[x,y],[0,0],5)")
#@test evaluate(ex1) == evaluate(ex2)

ex1 = sum(1/n^2,n,1,17)
ex2 = giac("sum(1/n^2,n,1,17)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = sum(1/n^2,n,1,17,2)
ex2 = giac("sum(1/n^2,n,1,17,2)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = sum(1/n^2,n,1,infinity)
ex2 = giac("sum(1/n^2,n,1,infinity)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = sum(cos(n*x),n)
ex2 = giac("sum(cos(n*x),n)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = curl([2*x*y,x*z,y*z],[x,y,z])
ex2 = giac("curl([2*x*y,x*z,y*z],[x,y,z])")
@test evaluate(ex1) == evaluate(ex2)

ex1 = divergence([x^2+y,x+z+y,z^3+x^2],[x,y,z])
ex2 = giac("divergence([x^2+y,x+z+y,z^3+x^2],[x,y,z])")
@test evaluate(ex1) == evaluate(ex2)

ex1 = grad(2*x^2*y-x*z^3,[x,y,z])
ex2 = giac("grad(2*x^2*y-x*z^3,[x,y,z])")
@test evaluate(ex1) == evaluate(ex2)

ex1 = hessian(2*x^2*y-x*z,[x,y,z])
ex2 = giac("hessian(2*x^2*y-x*z,[x,y,z])")
@test evaluate(ex1) == evaluate(ex2)

ex1 = taylor(sin(x)/x,x)
ex2 = giac("taylor(sin(x)/x,x)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = taylor(sin(x)/x,x=>0, 7)
ex2 = giac("taylor(sin(x)/x,x=0, 7)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = solve(x^2-3⩦1)
ex2 = giac("solve(x^2-3=1)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = solve([y-z⩦0,z-x⩦0,x-y⩦0,x-1+y+z⩦0],[x,y,z])
ex2 = giac("solve([y-z=0,z-x=0,x-y=0,x-1+y+z=0],[x,y,z])")
@test evaluate(ex1) == evaluate(ex2)

ex1 = zeros(x^2-4,x)
ex2 = giac("zeros(x^2-4,x)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = zeros([x^2-1,x^2-y^2],[x,y])
ex2 = giac("zeros([x^2-1,x^2-y^2],[x,y])")
@test evaluate(ex1) == evaluate(ex2)

ex1 = cSolve(x^4-1⩦0, x)
ex2 = giac("cSolve(x^4-1=0, x)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = cSolve([x^4-y^4,x+y⩦0,x^2⩦2*x],[x,y])
ex2 = giac("cSolve([x^4-y^4,x+y=0,x^2=2*x],[x,y])")
@test evaluate(ex1) == evaluate(ex2)

ex1 = cZeros(x^4-1, x)
ex2 = giac("cZeros(x^4-1, x)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = cZeros([x^2-1,x^2-y^2],[x,y])
ex2 = giac("cZeros([x^2-1,x^2-y^2],[x,y])")
@test evaluate(ex1) == evaluate(ex2)

ex1 = fSolve(cos(x)⩦x,x=>0)
ex2 = giac("fsolve(cos(x)=x,x=0)")
#Note: renamed fsolve->fSolve, because that is the name the HP Prime uses
@test evaluate(ex1) == evaluate(ex2)

ex1 = deSolve(diff(y(x),x,x)+y(x)⩦0, y)
ex2 = giac("deSolve(diff(y(x),x,x)+y(x)=0, y)")
@test evaluate(ex1) == evaluate(ex2)

#ex1 = deSolve(y''(x)+y(x)⩦0, y)
#ex2 = giac("deSolve(y''+y=0, y)")
##Note: standalone y,y',y'' in this context not yet defined
#@test evaluate(ex1) == evaluate(ex2)

ex1 = linsolve([x+y+z⩦1,x-y⩦2,2*x-z⩦3],[x,y,z])
ex2 = giac("linsolve([x+y+z=1,x-y=2,2*x-z=3],[x,y,z])")
@test evaluate(ex1) == evaluate(ex2)

ex1 = logcollect(log(x)+2*log(y))
ex2 = giac("lncollect(log(x)+2*log(y))")
#Note: renamed lncollect->logcollect because Giac/ln is Julia/log
@test evaluate(ex1) == evaluate(ex2)

ex1 = powexpand(2^(x+y))
ex2 = giac("powexpand(2^(x+y))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = texpand(sin(2*x)+exp(x+y))
ex2 = giac("texpand(sin(2*x)+exp(x+y))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = exp2pow(exp(3*log(x)))
ex2 = giac("exp2pow(exp(3*log(x)))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = pow2exp(x^y)
ex2 = giac("pow2exp(x^y)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = exp2trig(exp(im*x))
ex2 = giac("exp2trig(exp(i*x))")
#Note: Julia/im corresponds to Giac/i
@test evaluate(ex1) == evaluate(ex2)

ex1 = expexpand(exp(3*x+2*y))
ex2 = giac("expexpand(exp(3*x+2*y))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = asin2acos(acos(x)+asin(x))
ex2 = giac("asin2acos(acos(x)+asin(x))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = asin2atan(2*asin(x))
ex2 = giac("asin2atan(2*asin(x))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = sin2costan(sin(x))
ex2 = giac("sin2costan(sin(x))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = acos2asin(acos(x)+asin(x))
ex2 = giac("acos2asin(acos(x)+asin(x))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = acos2atan(2*acos(x))
ex2 = giac("acos2atan(2*acos(x))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = cos2sintan(cos(x))
ex2 = giac("cos2sintan(cos(x))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = atan2asin(atan(2*x))
ex2 = giac("atan2asin(atan(2*x))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = atan2acos(atan(2*x))
ex2 = giac("atan2acos(atan(2*x))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = tan2sincos(tan(x))
ex2 = giac("tan2sincos(tan(x))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = halftan(sin(x))
ex2 = giac("halftan(sin(x))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = trigsin(cos(x)^4+sin(x)^2)
ex2 = giac("trigsin(cos(x)^4+sin(x)^2)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = trigcos(cos(x)^4+sin(x)^2)
ex2 = giac("trigcos(cos(x)^4+sin(x)^2)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = trigtan(cos(x)^4+sin(x)^2)
ex2 = giac("trigtan(cos(x)^4+sin(x)^2)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = atrig2log(atan(x))
ex2 = giac("atrig2ln(atan(x))")
#Note: renamed atrig2ln->atrig2log because Giac/ln is Julia/log
@test evaluate(ex1) == evaluate(ex2)

ex1 = tlin(sin(x)^3)
ex2 = giac("tlin(sin(x)^3)")
@test evaluate(ex1) == evaluate(ex2)

ex1 = tcollect(sin(x)+cos(x))
ex2 = giac("tcollect(sin(x)+cos(x))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = trigexpand(sin(3*x))
ex2 = giac("trigexpand(sin(3*x))")
@test evaluate(ex1) == evaluate(ex2)

ex1 = trig2exp(sin(x))
ex2 = giac("trig2exp(sin(x))")
@test evaluate(ex1) == evaluate(ex2)

# to be continued ...


text:
# Load X and y variable
using JLD
data = load("nonLinear.jld")
(X,y,Xtest,ytest) = (data["X"],data["y"],data["Xtest"],data["ytest"])

# Compute number of training examples and number of features
(n,d) = size(X)

display(size(X))

ntrain = Int(floor(n/2))

Xtrain = X[1:ntrain, :]
ytrain = y[1:ntrain, :]
Xvalidate = X[ntrain+1:n, :]
yvalidate = y[ntrain+1:n, :]

display(size(Xtrain))
display(size(Xvalidate))

# Fit least squares model
#include("leastSquares.jl")
#model = leastSquares(X,y)

bestError = 10000000000000000
bestl = 0
bestsigma = 0

include("leastSquaresRBFL2.jl")

# Find best l, sigma values
for l = 0:0.1:3
  for sigma = 0.1:0.1:3
		model = leastSquaresRBFL2(Xtrain,ytrain,l, sigma)

		# Report the error on the validation set
		t = size(Xvalidate,1)
		yhat = model.predict(Xvalidate)
		validationError = sum((yhat - yvalidate).^2)/t
		@printf("l = %f, sigma = %f, ValidationError = %.2f\n", l, sigma, validationError)
    if validationError < bestError
      bestError = validationError
      bestl = l
      bestsigma = sigma
    end
  end
end

@printf("Best validation l = %f, sigma = %f, error = %f", bestl, bestsigma, bestError)

model = leastSquaresRBFL2(X,y,bestl,bestsigma)

# Report the error on the test set
t = size(Xtest,1)
yhat = model.predict(Xtest)
testError = sum((yhat - ytest).^2)/t
@printf("l = %f, sigma = %f, TestError = %.2f\n", bestl, bestsigma, testError)

# Plot model
using PyPlot
figure()
plot(X,y,"b.")
plot(Xtest,ytest,"g.")
Xhat = minimum(X):.1:maximum(X)
Xhat = reshape(Xhat,length(Xhat),1) # Make into an n by 1 matrix
yhat = model.predict(Xhat)
plot(Xhat,yhat,"r")
ylim((-300,400))
show()


text:
Base.@kwdef struct FieldInfo
    name::String
    ltype::Union{DataType,String,Symbol}
    rtype::Union{DataType,String,Symbol} = rtype
    get_returns_reference::Bool = true
end

struct TypeInfo
    name::String
    fields::Vector{FieldInfo}
end


text:
#EDAData: base types & core functions
#-------------------------------------------------------------------------------

#==Useful constants
===============================================================================#

#==Main data structures
===============================================================================#

#Last Line


text:
export
    istriangular,
    nthtriangular, ntriangular,
    alltriangular, sometriangular, exacttriangular,
    primitive_pythagorean_triples

istriangular(n::Integer) = isperfectsquare(8*n + 1)

nthtriangular(n::Integer) = (n * (n+1)) ÷ 2
ntriangular(n::Int, T::Type = Int) = collect(exacttriangular(n, T))

alltriangular(T::Type = Int) = TriangularIterator{T}()
sometriangular(xmax::T) where {T<:Integer} = @pipe alltriangular(T) |> takewhile(x -> x <= xmax, _)
exacttriangular(n::Int, T::Type = Int) = @pipe alltriangular(T) |> take(_, n)

struct TriangularIterator{T<:Integer}
end

Base.start(::TriangularIterator{T}) where {T<:Integer} = zero(T), one(T)
Base.next(::TriangularIterator{T}, state) where {T<:Integer} = begin
    s, n = state
    s += n
    s, (s, n + 1)
end
Base.done(::TriangularIterator, _) = false

Base.eltype(it::TriangularIterator) = Base.eltype(typeof(it))
Base.eltype(::Type{TriangularIterator{T}}) where {T<:Integer} = T

Base.iteratorsize(::TriangularIterator) = Base.IsInfinite()

primitive_pythagorean_triples(T::Type = Int) = PrimitivePythagoreanTripleIterator{T}()
primitive_pythagorean_triples(maxperimeter::T) where {T} = PrimitivePythagoreanTriplePerimeterIterator{T}(maxperimeter)

abstract type PrimitivePythagoreanTripleIteratorBase{T<:Integer} end

struct PrimitivePythagoreanTripleIterator{T<:Integer} <: PrimitivePythagoreanTripleIteratorBase{T}
end

struct PrimitivePythagoreanTriplePerimeterIterator{T<:Integer} <: PrimitivePythagoreanTripleIteratorBase{T}
    maximum::T
end

Base.start(::PrimitivePythagoreanTripleIterator{T}) where {T} = begin
    q = List{Array{T,1}}()
    n = push!(q, [3,4,5])
    q, n
end
Base.start(it::PrimitivePythagoreanTriplePerimeterIterator{T}) where {T} = begin
    q = List{Array{T,1}}()
    n = it.maximum >= 12 ? push!(q, [3,4,5]) : q.node
    q, n
end

Base.next(::PrimitivePythagoreanTripleIterator, state) = begin
    q, n = state
    n.data, (q, next!(q, n))
end
Base.next(it::PrimitivePythagoreanTriplePerimeterIterator, state) = begin
    q, n = state
    n.data, (q, next!(q, n, it.maximum))
end

Base.done(::PrimitivePythagoreanTripleIteratorBase, state) = begin
    q, n = state
    iseol(q, n)
end

Base.eltype(it::PrimitivePythagoreanTripleIteratorBase) = Base.eltype(typeof(it))
Base.eltype(::Type{PrimitivePythagoreanTripleIterator{T}}) where {T} = Array{T,1}
Base.eltype(::Type{PrimitivePythagoreanTriplePerimeterIterator{T}}) where {T} = Array{T,1}

Base.iteratorsize(::PrimitivePythagoreanTripleIterator) = Base.IsInfinite()
Base.iteratorsize(::PrimitivePythagoreanTriplePerimeterIterator) = Base.SizeUnknown()

iseol(q::List{T}, n::ListNode{T}) where {T} = n == q.node

next!(q::List{T}, n::ListNode{T}) where {T} = begin
    n = n.next
    iseol(q, n) || return n

    x, y, z = calculate(shift!(q)...)
    n = push!(q, x)
    push!(q, y)
    push!(q, z)
    n
end

next!(q::List{T}, n::ListNode{T}, maxsum) where {T} = begin
    n = n.next
    iseol(q, n) || return n

    while iseol(q, n) && !isempty(q)
        x, y, z = calculate(shift!(q)...)
        sum(x) <= maxsum && (nx = push!(q, x); n = nx)
        sum(y) <= maxsum && (ny = push!(q, y); n = iseol(q, n) ? ny : n)
        sum(z) <= maxsum && (nz = push!(q, z); n = iseol(q, n) ? nz : n)
    end
    n
end

calculate(a, b, c) = begin
    x = [ a - 2b + 2c,  2a - b + 2c,  2a - 2b + 3c] # const U = [ 1  2  2, -2 -1 -2, 2 2 3], [a b c] * U
    y = [ a + 2b + 2c,  2a + b + 2c,  2a + 2b + 3c] # const A = [ 1  2  2,  2  1  2, 2 2 3], [a b c] * A
    z = [-a + 2b + 2c, -2a + b + 2c, -2a + 2b + 3c] # const D = [-1 -2 -2,  2  1  2, 2 2 3], [a b c] * D
    x, y, z
end


text:
using PyCall
@pyimport math
@show math.cos(1) # 0.5403023058681398

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/UMAT.jl/blob/master/LICENSE

using UMAT, Test, Materials, Libdl, Tensors
const abaqus = [1 4 6; 0 2 5; 0 0 3]

pkg_dir = dirname(Base.find_package("UMAT"))
if Sys.iswindows()
    lib_dir = joinpath(pkg_dir,"..","deps","usr","bin")
else
    lib_dir = joinpath(pkg_dir,"..","deps","usr","lib")
end

material = UmatMaterial(NTENS=4, NSTATV=0, NPROPS=0,
                        lib_path=joinpath(lib_dir,"libelastic." * dlext))

strains = fromvoigt(SymmetricTensor{2,2},[1.0e-3, 0.0, 0.0, 0.0]; offdiagscale = 2.0, order=abaqus)
ref_stress = [282.6923171355886, 121.1538570784259, 121.1538570784259, 0.0]
material.ddrivers = UmatDriverState(NTENS=4, strain=strains)
UMAT.integrate_material!(material)
update_material!(material)

for i in 1:3
    @test_broken isapprox(tovoigt(material.variables.stress, order=abaqus)[i],ref_stress[i],atol=sqrt(eps()))
end


text:
export  excelFluxReport

function excelFluxReport( mpSim::MPsim, timeGrid::Vector{T},
    transitions::MP.TransitionType...; filename::AbstractString = "fluxReport",
    overwrite::Bool = true, timeFactor::Real = 12.0 ) where T <: Real

    filename, overwrite = setupFile( filename, overwrite )
    tStart = now()
    fluxReport = transitionFluxReport( mpSim, timeGrid, transitions... )
    reportGenerationTime = ( now() - tStart ).value / 1000.0

    if isempty( fluxReport )
        return
    end  # if isempty( fluxReport )

    writeTransFluxData( fluxReport, mpSim, timeFactor, reportGenerationTime,
        filename, overwrite )

end  # excelFluxReport( mpSim, timeGrid, transitions, filename, overwrite, timeFactor )

excelFluxReport( mpSim::MPsim, timeRes::Real,
    transitions::MP.TransitionType...; filename::AbstractString = "fluxReport",
    overwrite::Bool = true, timeFactor::Real = 12.0 ) =
    excelFluxReport( mpSim, MP.generateTimeGrid( mpSim, timeRes ),
    transitions...; filename = filename, overwrite = overwrite,
    timeFactor = timeFactor )

include( "private/transitionflux.jl" )


text:
"""
    BlackKarasinski()

[Black-Karasinski model](https://en.wikipedia.org/wiki/Black%E2%80%93Karasinski_model).
"""
struct BlackKarasinski <: Model end

"""
    evaluate(IRD, BlackKarasinski(), n = 12)

Evaluate interest rate derivative using `BlackKarasinski` model.

# Arguments
- `IRD::InterestRateDerivative`: interest rate derivative
- `n`: number of paths to simulate
"""
function evaluate(IRD::InterestRateDerivative, m::BlackKarasinski, n::Int64 = 12)
    Δt = IRD.t / n
    rates = [IRD.r]
    for i = 0:n
        dr = IRD.k * (IRD.θ - log(rates[end])) * Δt + IRD.σ * rates[end] * randn()
        append!(rates, rates[end] + dr)
    end
    return rates
end


text:
module Stiffness
export StiffnessData

mutable struct StiffnessData
grads::Array{Float64,3}
gq::Array{Float64,3}#grads of 6 basis functions at 3 quadrature points.
StiffnessData()=
    new(zeros(2,3,6),
        reshape([-1.0,-1.0,1.0,1.0,-1.0,-1.0,1.0,0.0,1.0,0.0,-1.0,0.0,0.0,
                 -1.0,0.0,1.0,0.0,1.0,0.0,-2.0,-2.0,-2.0,2.0,0.0,0.0,2.0,2.0,
                 2.0,2.0,0.0,0.0,2.0,-2.0,-2.0,-2.0,0],2,3,6)
        )
end

function op!(S::StiffnessData,x::Array{Float64,1},y::Array{Float64,1},
            m::Array{Float64,1})
    a11=-y[1] + y[3]
    a12= y[1] - y[2]
    a21= x[1] - x[3]
    a22=-x[1] + x[2]

    for f=1:6
        @simd for p=1:3
            S.grads[1,p,f]=a11*S.gq[1,p,f]+a12*S.gq[2,p,f]
            S.grads[2,p,f]=a21*S.gq[1,p,f]+a22*S.gq[2,p,f]
        end

        # this seems slower:
        # S.grads[1,1:3,f]=a11*S.gq[1,1:3,f]+a12*S.gq[2,1:3,f]
        # S.grads[2,1:3,f]=a21*S.gq[1,1:3,f]+a22*S.gq[2,1:3,f]
    end

    det= -(x[2] - x[3])*y[1] + (x[1] - x[3])*y[2] - (x[1] - x[2])*y[3]
    dv=1.0/(6.0*det)
    ii=1
    for i=1:6
        for j=1:i
            s=0.0
            @simd for k=1:3
                s+=S.grads[1,k,i]*S.grads[1,k,j]+S.grads[2,k,i]*S.grads[2,k,j]
            end
            m[ii]=dv*s
            # this seems slower:
            # m[ii]=dv*(dot(S.grads[1,:,i],S.grads[1,:,j])+
            #           dot(S.grads[2,:,i],S.grads[2,:,j]))
            ii+=1
        end
    end

    m
end
end

max_stars_repo_name: fluiddyn/BenchmarksPythonJuliaAndCo

max_issues_repo_name: fluiddyn/BenchmarksPythonJuliaAndCo

max_forks_repo_name: fluiddyn/BenchmarksPythonJuliaAndCo


text:
# ---------------------------------------------------
# Save the data for the figures in Chapter 2
# The figures themselves are created
# from GMT scripts in the ./GMT/ directory
# ---------------------------------------------------
module SaveFigureData

using NetCDF
using DelimitedFiles
using Statistics
using NCDatasets
function RunSaveFigureData(settings)
    println("Saving figure data...")
    fig_1_gmsl_usa(settings)
    fig_2_regional(settings)
    fig_4_gmsl(settings)
    fig_5_divergence(settings)
    return
end

function fig_1_gmsl_usa(settings)
    # Comparison of observed global-mean (Panel a) and Contiguous US (Panel b), Observed sea level, trajectory and scenarios
    # Panel a: GMSL
    GMSL_20c = ncread(settings["fn_proj_glb"],"GMSL_20c")
    GMSL_Altimetry = ncread(settings["fn_proj_glb"],"GMSL_Altimetry")
    GMSL_Trajectory = ncread(settings["fn_proj_glb"],"GMSL_Trajectory")
    save_ts_unc_gmt(settings["years"],GMSL_20c./1000,settings["dir_fig_1_gmsl_usa"]*"GMSL_20c.txt")
    save_ts_gmt(settings["years"],GMSL_Altimetry./1000,settings["dir_fig_1_gmsl_usa"]*"GMSL_Altimetry.txt")
    save_ts_unc_gmt(settings["years"],GMSL_Trajectory./1000,settings["dir_fig_1_gmsl_usa"]*"GMSL_Trajectory.txt")
    for scenario ∈ settings["NCA5_scenarios"]
        GMSL_scn = ncread(settings["fn_proj_glb"],"GMSL_"*scenario)
        save_ts_unc_gmt(settings["years"],GMSL_scn./1000,settings["dir_fig_1_gmsl_usa"]*"GMSL_"*scenario*".txt")
    end

    # Panel b: Contiguous US
    USA_rsl = ncread(settings["fn_proj_reg"],"rsl_obs")[1,:]
    USA_alt = ncread(settings["fn_proj_reg"],"alt_obs")[1,:]
    USA_trajectory = ncread(settings["fn_proj_reg"],"rsl_trajectory")[1,:,:]
    save_ts_gmt(settings["years"],USA_rsl./1000,settings["dir_fig_1_gmsl_usa"]*"USA_20c.txt")
    save_ts_gmt(settings["years"],USA_alt./1000,settings["dir_fig_1_gmsl_usa"]*"USA_Altimetry.txt")
    save_ts_unc_gmt(settings["years"],USA_trajectory./1000,settings["dir_fig_1_gmsl_usa"]*"USA_Trajectory.txt")

    for scenario ∈ settings["NCA5_scenarios"]
        USA_scn = ncread(settings["fn_proj_reg"],"rsl_total_"*scenario)[1,:,:]
        save_ts_unc_gmt(settings["years"],USA_scn./1000,settings["dir_fig_1_gmsl_usa"]*"USA_"*scenario*".txt")
    end
    return nothing
end

function fig_2_regional(settings)
    # Plot observations, trajectory and scenarios for each individual region
    years = ncread(settings["fn_proj_glb"],"years")
    reg_rsl = ncread(settings["fn_proj_reg"],"rsl_obs")
    reg_alt = ncread(settings["fn_proj_reg"],"alt_obs")
    reg_trajectory = ncread(settings["fn_proj_reg"],"rsl_trajectory")
    for (region_idx,region) ∈ enumerate(settings["regions"])
        save_ts_gmt(years,reg_rsl[region_idx,:]./1000,settings["dir_fig_2_regional"]*region*"_20c.txt")
        save_ts_gmt(years,reg_alt[region_idx,:]./1000,settings["dir_fig_2_regional"]*region*"_Altimetry.txt")
        save_ts_unc_gmt(years,reg_trajectory[region_idx,:,:]./1000,settings["dir_fig_2_regional"]*region*"_Trajectory.txt")
        for scenario ∈ settings["NCA5_scenarios"]
            reg_scn = ncread(settings["fn_proj_reg"],"rsl_total_"*scenario)[region_idx,:,:]
            save_ts_unc_gmt(years,reg_scn./1000,settings["dir_fig_2_regional"]*region*"_"*scenario*".txt")
        end
    end
    return nothing
end

function fig_4_gmsl(settings)
    # Observed global sea level change and its components (Panel a) and GMSL versus contiguous US sea level (Panel b)
    t = [1900:2018...]
    glac = zeros(length(t),3)
    GrIS = zeros(length(t),3)
    AIS = zeros(length(t),3)
    tws = zeros(length(t),3)
    steric = zeros(length(t),3)
    budget = zeros(length(t),3)
    GMSL = zeros(length(t),3)
    read_gmsl_var!("glac",glac,settings)
    read_gmsl_var!("GrIS",GrIS,settings)
    read_gmsl_var!("AIS",AIS,settings)
    read_gmsl_var!("tws",tws,settings)
    read_gmsl_var!("sum_of_contrib_processes",budget,settings)
    read_gmsl_var!("sum_of_contrib_processes",budget,settings)
    read_gmsl_var2!("global_average_sea_level_change",GMSL,settings)
    read_gmsl_var2!("global_average_thermosteric_sea_level_change",steric,settings)
    USA_20c = ncread(settings["fn_proj_reg"],"rsl_obs")[1,:]

    save_ts_unc_gmt(t,0.001.*(glac.-glac[101,2]),settings["dir_fig_4_gmsl"]*"glac.txt")
    save_ts_unc_gmt(t,0.001.*(GrIS.-GrIS[101,2]),settings["dir_fig_4_gmsl"]*"GrIS.txt")
    save_ts_unc_gmt(t,0.001.*(AIS.-AIS[101,2]),settings["dir_fig_4_gmsl"]*"AIS.txt")
    save_ts_unc_gmt(t,0.001.*(tws.-tws[101,2]),settings["dir_fig_4_gmsl"]*"tws.txt")
    save_ts_unc_gmt(t,0.001.*(steric.-steric[101,2]),settings["dir_fig_4_gmsl"]*"steric.txt")
    save_ts_unc_gmt(t,0.001.*(budget.-budget[101,2]),settings["dir_fig_4_gmsl"]*"budget.txt")
    save_ts_unc_gmt(t,0.001.*(GMSL.-GMSL[101,2]) ,settings["dir_fig_4_gmsl"]*"GMSL.txt")
    save_ts_gmt(settings["years"],0.001.*(USA_20c.-USA_20c[101]),settings["dir_fig_4_gmsl"]*"USA_MSL.txt")
    return nothing
end

function fig_5_divergence(settings)
    # Plot where each scenario diverges from the trajectory
    # This script uses 2σ instead of 1σ in the rest of the figures

    # Observed and trajectory time series
    GMSL_20c = ncread(settings["fn_proj_glb"],"GMSL_20c")
    GMSL_Trajectory = ncread(settings["fn_proj_glb"],"GMSL_Trajectory")
    save_ts_unc_gmt(settings["years"],0.001.*(GMSL_20c),settings["dir_fig_5_divergence"]*"GMSL_20c.txt")
    save_ts_unc_gmt(settings["years"],0.001.*(GMSL_Trajectory),settings["dir_fig_5_divergence"]*"GMSL_Trajectory.txt")

    # Panel a: relative to trajectory
    for scenario ∈ settings["NCA5_scenarios"]
        # For each scenario:
        # Read and save time series
        # Determine where this scenario starts diverging
        # from trajectory (Both above and below the trajectory)
        GMSL_scn = ncread(settings["fn_proj_glb"],"GMSL_"*scenario)
        save_ts_unc_gmt(settings["years"],GMSL_scn./1000,settings["dir_fig_5_divergence"]*"GMSL_"*scenario*".txt")
        diff_mean = @. (GMSL_scn[:,2] - GMSL_Trajectory[:,2])
        diff_err = @. 2 * sqrt((0.5*(GMSL_scn[:,3] - GMSL_scn[:,1]))^2 + (0.5*(GMSL_Trajectory[:,3] - GMSL_Trajectory[:,1]))^2)
        scn_above_traj = findfirst(@. (diff_mean - diff_err) >= 0)
        scn_below_traj = findfirst(@. (diff_mean + diff_err) <= 0)
        if !isnothing(scn_above_traj)
            p_plot = [settings["years"][scn_above_traj] 0.001*GMSL_scn[scn_above_traj,2] 1]
            l_plot = [settings["years"][scn_above_traj] 0.001*GMSL_scn[scn_above_traj,2] ; settings["years"][scn_above_traj] -0.05]
            writedlm(settings["dir_fig_5_divergence"]*"div_p_"*scenario*".txt",p_plot,";")
            writedlm(settings["dir_fig_5_divergence"]*"div_l_"*scenario*".txt",l_plot,";")
        elseif !isnothing(scn_below_traj)
            p_plot = [settings["years"][scn_below_traj] 0.001*GMSL_scn[scn_below_traj,2] 1]
            l_plot = [settings["years"][scn_below_traj] 0.001*GMSL_scn[scn_below_traj,2] ; settings["years"][scn_below_traj] -0.05 ]
            writedlm(settings["dir_fig_5_divergence"]*"div_p_"*scenario*".txt",p_plot,";")
            writedlm(settings["dir_fig_5_divergence"]*"div_l_"*scenario*".txt",l_plot,";")
        else
            println("ERROR: no divergence for "*scenario*" from trajectory!")
        end
    end

    # Panel b: relative to Intermediate
    GMSL_int = ncread(settings["fn_proj_glb"],"GMSL_Int")
    for scenario ∈ settings["NCA5_scenarios"]
        # For each scenario:
        # Read and save time series
        # Determine where this scenario starts diverging
        # from trajectory (Both above and below the trajectory)
        GMSL_scn = ncread(settings["fn_proj_glb"],"GMSL_"*scenario)
        diff_mean = @. (GMSL_scn[:,2] - GMSL_int[:,2])
        diff_err = @. 2 * sqrt((0.5*(GMSL_scn[:,3] - GMSL_scn[:,1]))^2 + (0.5*(GMSL_int[:,3] - GMSL_int[:,1]))^2)
        scn_above_traj = findfirst(@. (diff_mean - diff_err) >= 0)
        scn_below_traj = findfirst(@. (diff_mean + diff_err) <= 0)
        if !isnothing(scn_above_traj)
            p_plot = [settings["years"][scn_above_traj] 0.001*GMSL_scn[scn_above_traj,2] 1]
            l_plot = [settings["years"][scn_above_traj] 0.001*GMSL_scn[scn_above_traj,2] ; settings["years"][scn_above_traj] -0.05]
            writedlm(settings["dir_fig_5_divergence"]*"div_p_Int_"*scenario*".txt",p_plot,";")
            writedlm(settings["dir_fig_5_divergence"]*"div_l_Int_"*scenario*".txt",l_plot,";")
        elseif !isnothing(scn_below_traj)
            p_plot = [settings["years"][scn_below_traj] 0.001*GMSL_scn[scn_below_traj,2] 1]
            l_plot = [settings["years"][scn_below_traj] 0.001*GMSL_scn[scn_below_traj,2] ; settings["years"][scn_below_traj] -0.05 ]
            writedlm(settings["dir_fig_5_divergence"]*"div_p_Int_"*scenario*".txt",p_plot,";")
            writedlm(settings["dir_fig_5_divergence"]*"div_l_Int_"*scenario*".txt",l_plot,";")
        else
            println("ERROR: no divergence for "*scenario*" from intermediate!")
        end
    end
    return nothing
end

function figure_8_US_map(settings)
    fh = Dataset(settings["fn_proj_gri"],"r")
    year_idx = findfirst(fh["years"][:].==2100)
    rsl = convert.(Float32,fh["rsl_total_Int"][:,:,year_idx,2])
    @. rsl[rsl <-31900] = NaN32
    @. rsl *= 0.00328084 ## millimeter to feet
    lon = fh["lon"][:]
    lat = fh["lat"][:]
    close(fh)

    fh = Dataset(settings["dir_gmt"] * "fig_8_US_map/fig_8_US_map.nc","c")
    defDim(fh,"lon", length(lon))
    defDim(fh,"lat", length(lat))
    defVar(fh,"lon",lon,("lon",),deflatelevel=5)
    defVar(fh,"lat",lat,("lat",),deflatelevel=5)
    defVar(fh,"rsl",rsl,("lon","lat",),deflatelevel=5)
    close(fh)

end

# -----------------
# General functions
# -----------------
function read_gmsl_var!(varname,store,settings)
    store[:,1] = ncread(settings["fn_gmsl_20c_mean"],varname*"_lower")
    store[:,2] = ncread(settings["fn_gmsl_20c_mean"],varname*"_mean")
    store[:,3] = ncread(settings["fn_gmsl_20c_mean"],varname*"_upper")
    return nothing
end

function read_gmsl_var2!(varname,store,settings)
    store[:,1] = ncread(settings["fn_gmsl_20c_mean"],varname*"_lower")
    store[:,2] = ncread(settings["fn_gmsl_20c_mean"],varname)
    store[:,3] = ncread(settings["fn_gmsl_20c_mean"],varname*"_upper")
    return nothing
end

function save_ts_unc_gmt(years,ts,fn)
    acc_idx = isfinite.(ts[:,2])
    sdata = [years[acc_idx] ts[acc_idx,2] ts[acc_idx,1] ts[acc_idx,3]]
    writedlm(fn,sdata,";")
    return nothing
end

function save_ts_gmt(years,ts,fn)
    acc_idx = isfinite.(ts)
    sdata = [years[acc_idx] ts[acc_idx]]
    writedlm(fn,sdata,";")
    return nothing
end

end


text:
@testset "Edmonds Karp" begin
    # Construct DiGraph
    flow_graph = lg.DiGraph(8)

    # Load custom dataset
    flow_edges = [
        (1, 2, 10), (1, 3, 5), (1, 4, 15), (2, 3, 4), (2, 5, 9),
        (2, 6, 15), (3, 4, 4), (3, 6, 8), (4, 7, 16), (5, 6, 15),
        (5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6), (7, 8, 10)
    ]

    for fg in (flow_graph, lg.DiGraph{UInt8}(flow_graph), lg.DiGraph{Int16}(flow_graph))
        capacity_matrix = zeros(Int, 8, 8)
        for e in flow_edges
            u, v, f = e
            lg.add_edge!(fg, u, v)
            capacity_matrix[u, v] = f
        end
        residual_graph = @inferred(LightGraphsFlows.residual(fg))

        # Test with default distances
        @test @inferred(LightGraphsFlows.edmonds_karp_impl(residual_graph, 1, 8, LightGraphsFlows.DefaultCapacity(residual_graph)))[1] == 3

        # Test with capacity matrix
        @test @inferred(LightGraphsFlows.edmonds_karp_impl(residual_graph, 1, 8, capacity_matrix))[1] == 28

        # Test the types of the values returned by fetch_path
        function test_find_path_types(residual_graph, s, t, flow_matrix, capacity_matrix)
            v, P, S, flag = LightGraphsFlows.fetch_path(residual_graph, s, t, flow_matrix, capacity_matrix)
            @test typeof(P) == Vector{Int}
            @test typeof(S) == Vector{Int}
            @test typeof(flag) == Int
            @test typeof(v) == eltype(residual_graph)
        end

        # Test the value of the flags returned.
        function test_find_path_disconnected(residual_graph, s, t, flow_matrix, capacity_matrix)
            h = copy(residual_graph)
            for dst in collect(lg.neighbors(residual_graph, s))
                lg.rem_edge!(residual_graph, s, dst)
            end
            v, P, S, flag = LightGraphsFlows.fetch_path(residual_graph, s, t, flow_matrix, capacity_matrix)
            @test flag == 1
            for dst in collect(lg.neighbors(h, t))
                lg.rem_edge!(h, t, dst)
            end
            v, P, S, flag = LightGraphsFlows.fetch_path(h, s, t, flow_matrix, capacity_matrix)
            @test flag == 0
            for i in collect(lg.inneighbors(h, t))
                lg.rem_edge!(h, i, t)
            end
            v, P, S, flag = LightGraphsFlows.fetch_path(h, s, t, flow_matrix, capacity_matrix)
            @test flag == 2
        end
        flow_matrix = zeros(Int, lg.nv(residual_graph), lg.nv(residual_graph))
        test_find_path_types(residual_graph, 1, 8, flow_matrix, capacity_matrix)
        test_find_path_disconnected(residual_graph, 1, 8, flow_matrix, capacity_matrix)
    end
end


text:
import StableRNGs: StableRNG

# create scattered data
npts = 50
xmin = 0.0
xspan = 3/2*π
x = collect(range(xmin, xmin+xspan, length=npts))
rng = StableRNG(655)
x = x + xspan/npts*(rand(rng,npts) .- 0.5)
# select a subset randomly
idx = unique(rand(rng, collect(eachindex(x)), 20))
t = x[unique(idx)]
npts = length(t)
ut = sin.(t)
stdev = 1e-1*maximum(ut)
u = ut + stdev*randn(rng, npts)
# data must be ordered if t̂ is not provided
idx = sortperm(t)
tₒ = t[idx]
uₒ = u[idx]

tolerance = 1e-6

@testset "Direct smoothing" begin
    # fixed with default λ = 1.0
    A = RegularizationSmooth(uₒ,tₒ; alg=:fixed)
    ans = [0.6456173647252937 0.663974701324226 0.7631218523665086 0.778654700697601 0.7489958320589535 0.7319087707475104 0.6807082599508811 0.6372557895089508 0.5832859790765743 0.5021274805916013 0.3065928203396211 0.1353332321156384 -0.3260000640060584 -0.6557906092739154 -0.9204882447932498]'
    @test isapprox(A.û, ans, rtol=tolerance)
    # non-default d and λ
    A = RegularizationSmooth(uₒ,tₒ, 4; λ=1e-2, alg=:fixed)
    ans = [0.19865190868740357 0.2885349151737291 0.6756699442978945 0.9165887141895426 0.9936113717653254 1.0042825002191034 0.9768118192829827 0.9184595331808411 0.8214983284892922 0.6538356458824783 0.28295521578898 0.018060767871253963 -0.5301723647977373 -0.8349855890541111 -1.1085048455468356]'
    @test isapprox(A.û, ans, rtol=tolerance)
    # GCV (default) to determine λ
    A = RegularizationSmooth(uₒ,tₒ)
    @test isapprox(A.λ, 0.12788440382063268, rtol=tolerance)
    ans = [0.21974931848164914 0.2973284508009968 0.6908546278415386 0.9300465474303226 0.9741453042418977 0.9767572556868123 0.9432951659303452 0.8889834700087442 0.804842790047182 0.6603217445567791 0.30341652659101737 0.05924456463634589 -0.5239939779242144 -0.8421768233191822 -1.107517099580091]'
    @test isapprox(A.û, ans, rtol=tolerance)
    # L-curve to determine λ
    A = RegularizationSmooth(uₒ,tₒ; alg=:L_curve)
    @test isapprox(A.λ, 0.3339750913481082, rtol=tolerance)
    ans = [0.3038273346447746 0.3695574797289306 0.7261740581164446 0.8945169297047363 0.9100773683859941 0.9048101771171335 0.8663722295718835 0.8188893880450858 0.7520030529566242 0.641893306140736 0.35912131844587486 0.12553463033078033 -0.45446129364457777 -0.8246966349034558 -1.118321479210823]'
    @test isapprox(A.û, ans, rtol=tolerance)
end

@testset "Smoothing with weights" begin
    # midpoint rule integration
    A = RegularizationSmooth(uₒ,tₒ, nothing, :midpoint)
    @test isapprox(A.λ, 0.10768759624567165, rtol=tolerance)
    ans = [0.30685734533678666 0.3637106676306693 0.6653926790689593 0.9056721139798669 0.9738743511381565 0.9822029085679189 0.9503276381232628 0.8954395464903986 0.8024964523151614 0.6415834595513483 0.2833951659794013 0.052823783087081515 -0.5334066052968335 -0.8406485796257146 -1.0983153546407178]'
    @test isapprox(A.û, ans, rtol=tolerance)
    # arbitrary weights for wls
    A = RegularizationSmooth(uₒ,tₒ, nothing, collect(1:npts))
    @test isapprox(A.λ, 0.1447077651169621, rtol=tolerance)
    ans = [0.25185799867326253 0.32130333169345704 0.6679088195845279 0.927260505998722 0.9914594632495953 1.0028090695917116 0.9749355187531624 0.9254679051053563 0.8291064785367862 0.6493316082309926 0.2324329561930743 0.07635651700236243 -0.5564454824473795 -0.817911933888329 -1.1130193817730065]'
    @test isapprox(A.û, ans, rtol=tolerance)
    # arbitrary weights for wls and wr
    nhalf = Int(floor(npts/2))
    wls = vcat(ones(nhalf),10*ones(npts-nhalf))
    wr = collect(1:npts-2)
    A = RegularizationSmooth(uₒ,tₒ, nothing, wls, wr)
    @test isapprox(A.λ, 0.01920769986569919, rtol=tolerance)
    ans = [0.16161677980889166 0.34630920002600707 0.6469896059918103 0.9219255811754361 1.0001382485019457 1.0260666135098766 1.0278670301547874 0.9823839714037528 0.8512856706983483 0.6416080323126064 0.1892473613483678 0.10264529505906245 -0.5680977857976282 -0.8075007645180948 -1.1168524120952026]'
    @test isapprox(A.û, ans, rtol=tolerance)
end

@testset "Smoothing with t̂ provided" begin  ##### TBD
    N̂ = 20
    t̂ = collect(range(xmin, xmin+xspan, length=N̂))
    # with t̂, no weights
    A = RegularizationSmooth(u,t, t̂)
    @test isapprox(A.λ, 0.138273889585313, rtol=tolerance)
    ans = [0.21626377852882872 0.39235926952322575 0.5573848799950002 0.7072474496656729 0.8361906119247042 0.9313473799797176 0.9809844353757837 0.9750833208625507 0.9096038940899813 0.7816929736202427 0.6052694276527628 0.4015903497629387 0.1913719025253403 -0.01979786871512895 -0.23400354001942947 -0.44481229967011127 -0.6457913359497256 -0.8405146928672158 -1.0367229293434395 -1.2334090099343238]'
    @test isapprox(A.û, ans, rtol=tolerance)
    # t̂ and wls
    A = RegularizationSmooth(u,t, t̂, collect(1:npts))
    @test isapprox(A.λ, 0.26746430253489195, rtol=tolerance)
    ans = [0.3118247878815087 0.44275860852897864 0.5705834985506882 0.6979119448253899 0.8234189540704866 0.9289458273102476 0.9970803470992273 1.0071205506077525 0.9443157518324818 0.7954860908242515 0.5847385548859145 0.34813493129868633 0.1237494751337505 -0.0823517516424196 -0.28265170846635246 -0.4760833187699964 -0.6615795059853024 -0.844779821396189 -1.0341162283806349 -1.225270266213379]'
    @test isapprox(A.û, ans, rtol=tolerance)
    # t̂, wls, and wr
    nhalf = Int(floor(npts/2))
    wls = vcat(ones(nhalf),10*ones(npts-nhalf))
    wr = collect(1:N̂-2)
    A = RegularizationSmooth(u,t, t̂, wls, wr)
    @test isapprox(A.λ, 0.04542202002453003, rtol=tolerance)
    ans = [0.27982685119190176 0.4629485181120683 0.5612727224498559 0.6647076034312198 0.790945652904372 0.9094680579558058 1.0067349357353885 1.0543162288524837 0.9890994274368918 0.8088681588493614 0.5676793701053269 0.3129632437767546 0.08578111296361303 -0.11476994629546472 -0.30745805411723526 -0.49422228232466525 -0.6767431598309912 -0.858727390770621 -1.0443219211010253 -1.2309211555227055]'
    @test isapprox(A.û, ans, rtol=tolerance)
end


text:
@doc doc"""
    usinc(θ::Real)

Unnormalized version of `sinc` function, i.e.
$\operatorname{usinc}(\theta) = \frac{\sin(\theta)}{\theta}$. This is
equivalent to `sinc(θ/π)`.
"""
@inline usinc(θ::Real) = θ == 0 ? one(θ) : isinf(θ) ? zero(θ) : sin(θ) / θ

@doc doc"""
    usinc_from_cos(x::Real)

Unnormalized version of `sinc` function, i.e.
$\operatorname{usinc}(\theta) = \frac{\sin(\theta)}{\theta}$, computed from
$x = cos(\theta)$.
"""
@inline function usinc_from_cos(x::Real)
    if x >= 1
        return one(x)
    elseif x <= -1
        return zero(x)
    else
        return sqrt(1 - x^2) / acos(x)
    end
end

"""
    eigen_safe(x)

Compute the eigendecomposition of `x`. If `x` is a `StaticMatrix`, it is
converted to a `Matrix` before the decomposition.
"""
@inline eigen_safe(x; kwargs...) = eigen(x; kwargs...)

@inline function eigen_safe(x::StaticMatrix; kwargs...)
    s = size(x)
    E = eigen!(Matrix(parent(x)); kwargs...)
    return Eigen(SizedVector{s[1]}(E.values), SizedMatrix{s...}(E.vectors))
end

"""
    log_safe(x)

Compute the matrix logarithm of `x`. If `x` is a `StaticMatrix`, it is
converted to a `Matrix` before computing the log.
"""
@inline log_safe(x) = log(x)

@inline function log_safe(x::StaticMatrix)
    s = Size(x)
    return SizedMatrix{s[1],s[2]}(log(Matrix(parent(x))))
end

"""
    select_from_tuple(t::NTuple{N, Any}, positions::Val{P})

Selects elements of tuple `t` at positions specified by the second argument.
For example `select_from_tuple(("a", "b", "c"), Val((3, 1, 1)))` returns
`("c", "a", "a")`.
"""
@generated function select_from_tuple(t::NTuple{N, Any}, positions::Val{P}) where {N, P}
    for k in P
        (k < 0 || k > N) && error("positions must be between 1 and $N")
    end
    return Expr(:tuple, [Expr(:ref, :t, k) for k in P]...)
end

"""
    ziptuples(a, b)

Zips tuples `a` and `b` in a fast, type-stable way. If they have different
lengths, the result is trimmed to the length of the shorter tuple.
"""
@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}) where {N,M}
    ex = Expr(:tuple)
    for i = 1:min(N, M)
        push!(ex.args, :((a[$i], b[$i])))
    end
    ex
end

"""
    ziptuples(a, b, c)

Zips tuples `a`, `b` and `c` in a fast, type-stable way. If they have different
lengths, the result is trimmed to the length of the shorter tuple.
"""
@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}, c::NTuple{L,Any}) where {N,M,L}
    ex = Expr(:tuple)
    for i = 1:min(N, M, L)
        push!(ex.args, :((a[$i], b[$i], c[$i])))
    end
    ex
end

"""
    ziptuples(a, b, c, d)

Zips tuples `a`, `b`, `c` and `d` in a fast, type-stable way. If they have
different lengths, the result is trimmed to the length of the shorter tuple.
"""
@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}, c::NTuple{L,Any}, d::NTuple{K,Any}) where {N,M,L,K}
    ex = Expr(:tuple)
    for i = 1:min(N, M, L, K)
        push!(ex.args, :((a[$i], b[$i], c[$i], d[$i])))
    end
    ex
end

"""
    ziptuples(a, b, c, d, e)

Zips tuples `a`, `b`, `c`, `d` and `e` in a fast, type-stable way. If they have
different lengths, the result is trimmed to the length of the shorter tuple.
"""
@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}, c::NTuple{L,Any}, d::NTuple{K,Any}, e::NTuple{J,Any}) where {N,M,L,K,J}
    ex = Expr(:tuple)
    for i = 1:min(N, M, L, K, J)
        push!(ex.args, :((a[$i], b[$i], c[$i], d[$i], e[$i])))
    end
    ex
end

# conversion of SizedArray from StaticArrays.jl
"""
    SizedAbstractArray{Tuple{dims...}}(array)

Wraps an `Array` with a static size, so to take advantage of the (faster)
methods defined by the static array package. The size is checked once upon
construction to determine if the number of elements (`length`) match, but the
array may be reshaped.
(Also, `Size(dims...)(array)` acheives the same thing)
"""
struct SizedAbstractArray{S<:Tuple, T, N, M, TData<:AbstractArray{T,M}} <: StaticArray{S, T, N}
    data::TData

    function SizedAbstractArray{S, T, N, M, TData}(a::TData) where {S, T, N, M, TData<:AbstractArray}
        if length(a) != StaticArrays.tuple_prod(S)
            error("Dimensions $(size(a)) don't match static size $S")
        end
        new{S,T,N,M,TData}(a)
    end

    function SizedAbstractArray{S, T, N, 1}(::UndefInitializer) where {S, T, N}
        new{S, T, N, 1, Array{T, 1}}(Array{T, 1}(undef, StaticArrays.tuple_prod(S)))
    end
    function SizedAbstractArray{S, T, N, N}(::UndefInitializer) where {S, T, N}
        new{S, T, N, N, Array{T, N}}(Array{T, N}(undef, size_to_tuple(S)...))
    end
end

@inline SizedAbstractArray{S,T,N}(a::TData) where {S,T,N,M,TData<:AbstractArray{T,M}} = SizedAbstractArray{S,T,N,M,TData}(a)
@inline SizedAbstractArray{S,T}(a::TData) where {S,T,M,TData<:AbstractArray{T,M}} = SizedAbstractArray{S,T,StaticArrays.tuple_length(S),M,TData}(a)
@inline SizedAbstractArray{S}(a::TData) where {S,T,M,TData<:AbstractArray{T,M}} = SizedAbstractArray{S,T,StaticArrays.tuple_length(S),M,TData}(a)

@inline SizedAbstractArray{S,T,N}(::UndefInitializer) where {S,T,N} = SizedAbstractArray{S,T,N,N}(undef)
@inline SizedAbstractArray{S,T}(::UndefInitializer) where {S,T} = SizedAbstractArray{S,T,StaticArrays.tuple_length(S),StaticArrays.tuple_length(S)}(undef)
@generated function (::Type{SizedAbstractArray{S,T,N,M,TData}})(x::NTuple{L,Any}) where {S,T,N,M,TData,L}
    if L != StaticArrays.tuple_prod(S)
        error("Dimension mismatch")
    end
    exprs = [:(a[$i] = x[$i]) for i = 1:L]
    return quote
        $(Expr(:meta, :inline))
        a = SizedAbstractArray{S,T,N,M}(undef)
        @inbounds $(Expr(:block, exprs...))
        return a
    end
end

@inline SizedAbstractArray{S,T,N}(x::Tuple) where {S,T,N} = SizedAbstractArray{S,T,N,N}(x)
@inline SizedAbstractArray{S,T}(x::Tuple) where {S,T} = SizedAbstractArray{S,T,StaticArrays.tuple_length(S),StaticArrays.tuple_length(S)}(x)
@inline SizedAbstractArray{S}(x::NTuple{L,T}) where {S,T,L} = SizedAbstractArray{S,T,StaticArrays.tuple_length(S),StaticArrays.tuple_length(S)}(x)

# Overide some problematic default behaviour
@inline convert(::Type{SA}, sa::SizedAbstractArray) where {SA<:SizedAbstractArray} = SA(sa.data)
@inline convert(::Type{SA}, sa::SA) where {SA<:SizedAbstractArray} = sa

# Back to Array (unfortunately need both convert and construct to overide other methods)
import Base.Array
@inline Array(sa::SizedAbstractArray) = Array(sa.data)
@inline Array{T}(sa::SizedAbstractArray{S,T}) where {T,S} = Array(sa.data)
@inline Array{T,N}(sa::SizedAbstractArray{S,T,N}) where {T,S,N} = Array(sa.data)

@inline convert(::Type{AbstractArray}, sa::SizedAbstractArray) = sa.data
@inline convert(::Type{AbstractArray{T}}, sa::SizedAbstractArray{S,T}) where {T,S} = sa.data
@inline convert(::Type{AbstractArray{T,N}}, sa::SizedAbstractArray{S,T,N}) where {T,S,N} = sa.data

Base.@propagate_inbounds getindex(a::SizedAbstractArray, i::Int) = getindex(a.data, i)
Base.@propagate_inbounds setindex!(a::SizedAbstractArray, v, i::Int) = setindex!(a.data, v, i)

SizedAbstractVector{S,T,M} = SizedAbstractArray{Tuple{S},T,1,M}
@inline SizedAbstractVector{S}(a::TData) where {S,T,M,TData<:AbstractArray{T,M}} = SizedAbstractArray{Tuple{S},T,1,M,TData}(a)
@inline SizedAbstractVector{S}(x::NTuple{L,T}) where {S,T,L} = SizedAbstractArray{Tuple{S},T,1,1,Vector{T}}(x)

SizedAbstractMatrix{S1,S2,T,M} = SizedAbstractArray{Tuple{S1,S2},T,2,M}
@inline SizedAbstractMatrix{S1,S2}(a::TData) where {S1,S2,T,M,TData<:AbstractArray{T,M}} = SizedAbstractArray{Tuple{S1,S2},T,2,M,TData}(a)
@inline SizedAbstractMatrix{S1,S2}(x::NTuple{L,T}) where {S1,S2,T,L} = SizedAbstractArray{Tuple{S1,S2},T,2,2,Matrix{2}}(x)

"""
    Size(dims)(array)

Creates a `SizedAbstractArray` wrapping `array` with the specified statically-known
`dims`, so to take advantage of the (faster) methods defined by the static array
package.
"""
(::Size{S})(a::AbstractArray) where {S} = SizedAbstractArray{Tuple{S...}}(a)

function promote_rule(::Type{<:SizedAbstractArray{S,T,N,M,TDataA}}, ::Type{<:SizedAbstractArray{S,U,N,M,TDataB}}) where {S,T,U,N,M,TDataA,TDataB}
    SizedAbstractArray{S,promote_type(T,U),N,M,T,promote_type(TDataA, TDataB)}
end

@inline copy(a::SizedAbstractArray) = typeof(a)(copy(a.data))

similar(::Type{<:SizedAbstractArray{S,T,N,M}},::Type{T2}) where {S,T,N,M,T2} = SizedAbstractArray{S,T2,N,M}(undef)
similar(::Type{SA},::Type{T},s::Size{S}) where {SA<:SizedAbstractArray,T,S} = sizedabstractarray_similar_type(T,s,StaticArrays.length_val(s))(undef)
sizedabstractarray_similar_type(::Type{T},s::Size{S},::Type{Val{D}}) where {T,S,D} = SizedAbstractArray{Tuple{S...},T,D,length(s)}

"""
    size_to_tuple(::Type{S}) where S<:Tuple

Converts a size given by `Tuple{N, M, ...}` into a tuple `(N, M, ...)`.
"""
Base.@pure function size_to_tuple(::Type{S}) where S<:Tuple
    return tuple(S.parameters...)
end


text:
const NODES_VISITED = Ref{Int}(0)
const POINTS_VISITED = Ref{Int}(0)
const POINTS_VISITED_UNCHECKED = Ref{Int}(0)

add_node_visited(n::Int) = NODES_VISITED[] += n
add_point_visited(n::Int) = POINTS_VISITED[] += n
function add_point_visited_unchecked(n::Int)
    POINTS_VISITED[] += n
    POINTS_VISITED_UNCHECKED[] += n
end

function reset_stats()
    NODES_VISITED[] = 0
    POINTS_VISITED[] = 0
    return
end

function print_stats()
    println("Nodes visited: $(NODES_VISITED[])")
    println("Points visited: $(POINTS_VISITED[]), out of these: $(POINTS_VISITED_UNCHECKED[]) unchecked.")
end

macro POINT(n) end
macro POINT_UNCHECKED(n) end
macro NODE(n) end

if DEBUG
    println("Debugging for NearestNeighbors is active")
    @eval begin
        macro POINT(n)
            :(add_point_visited($n))
        end

        macro NODE(n)
            :(add_node_visited($n))
        end

        macro POINT_UNCHECK(n)
            :(add_point_visited_unchecked($n))
        end
    end
end


text:
module SimpleTest
export simple_operation
export start_of_opt
using JuMP
using Ipopt
using Clp
using Dates
using MathOptInterface
using CSV
using Documenter
using Literate
include("testLP.jl")
include("testLPMain.jl")
include("simple_operation.jl")
end


text:
module aoc_2019_01

ingest(path) = [parse(Int, l) for l in readlines(path)]

function calculate_fuel_cost(mass)
    return (mass ÷ 3) - 2
end

function calculate_fuel_cost_recursive(mass)
    fuel = (mass ÷ 3) - 2
    fuel > 0 || return 0
    return fuel + calculate_fuel_cost_recursive(fuel)
end

function solve(path, calculator)
    input_values = ingest(path)
    return sum(calculator(m) for m in input_values)
end

solve_part_one(path) = solve(path, calculate_fuel_cost)
solve_part_two(path) = solve(path, calculate_fuel_cost_recursive)

end # module

max_stars_repo_name: BsideAlgorithms/CodingDojo

max_issues_repo_name: BsideAlgorithms/CodingDojo

max_forks_repo_name: BsideAlgorithms/CodingDojo


text:
#File: fStd.jl
#Author:
#Date: 29-June-2020

#STANDARD DEVIATION
function fStd(x,flag,f)
#--------------------------------------------------------------------------
#x:     returns vector
#flag:  0 = sample, 1 = population
#f:     reporting frequency e.g. 12 (monthly)
#sd:    standard deviation (sample or population)
#--------------------------------------------------------------------------

    fStd = std(x,flag).*sqrt(f);
    sd = fStd;
    return sd
end


text:
using Measures
include("helper_functions.jl")

begin
    λ1 = 1
    λ2 = 1
    α = 1
    μ(α) = 1 / (1.5 * (α))
    ν(α) = 1 / (0.7 * (α))
    parameters = (λ1, λ2, α, μ, ν)
    opponent_parameters = (λ2, λ1, α, μ, ν)
    p, q = equilibrium_probability(parameters...)
    iter = 0:0.001:1
end

cost(p, q, parameters...) ≈ cost(q, p, opponent_parameters...)
costs = cost.(iter, q, parameters...)

function find_asymptotes(vec, iter)
    sv = sign.(vec)
    iter[sv .!= circshift(sv, 1)]
end

asymptotes = find_asymptotes(costs, iter)

nash_plot = plot(iter[costs .> 0], costs[costs .> 0], label = "Response Time")
plot!(ylims = (0, 25), xlims = (0, 1), size = (600, 250), palette = :tab10)
plot!([p, p], [0, cost(p, q, parameters...)], ls = :dash, lc = :black, label = "Minimum")
plot!([0, p], [cost(p, q, parameters...), cost(p, q, parameters...)], ls = :dash, lc = :black, label = "")
vline!([asymptotes], ls = :dash, label = "Asymptotes")
plot!(xlabel = "Proportion Sent to Public Server", ylabel = "Expected Response Time", legend = :topleft, margin = 2.5mm)
savefig("outputs/nash_plot_poster.pdf")


text:
using ConditionalDists
using DistributionsAD
#using Distributions: logpdf
using Flux

batch = 10000
xlength = 300
zlength = 150

f = Dense(zlength,xlength*2)
d = TuringDiagMvNormal(zeros(Float32,xlength), ones(Float32,xlength))
p = ConditionalMeanVarMvNormal(d,f)

z = randn(Float32, zlength, batch)
x = randn(Float32, xlength, batch)

@time logpdf(p, x, z)
@time logpdf(p, x, z)
@time logpdf(p, x, z)

@info "done"


text:
const initialcondition_types = Dict(
    :class=>String, :region=>String, :built=>Int, :optioned=>Int
)

const period_types = Dict(
    :name=>String, :weight=>Float64
)

const scenarios_tree_types = Dict(
    :name=>String,
    :parent=>String,
    :probability=>Float64
)

# TODO: Thermal / variable / storage technical types?

const scenarios_resource_param_types = Dict(
    :class=>String, :region=>String,
    :variablecost=>Float64, :fixedcost=>Float64,
    :startupcost=>Float64, :shutdowncost=>Float64,
    :optioncost=>Float64, :buildcost=>Float64, :retirementcost=>Float64,
    :buildleadtime=>Int, :optionleadtime=>Int,
    :newoptionslimit=>Int, :newbuildslimit=>Int, :rec=>Float64
)

const scenarios_transmission_param_types = Dict(
    :interface=>String, :limit=>Float64)

const scenarios_market_param_types = Dict(
    :region=>String, :pricecap=>Float64)

const scenarios_capacitymarket_param_types = Dict(
    :maxprice=>Float64, :midprice=>Float64,
    :maxpricepeakscalar=>Float64, :midpricepeakscalar=>Float64,
    :zeropricepeakscalar=>Float64)


text:
# TODO Instead of saving to a file, send content via this message
function view_profile(data::Union{Nothing,Vector{UInt}} = nothing, period::Union{Nothing,UInt64} = nothing; kwargs...)
    filename = string(tempname(), ".cpuprofile")

    ChromeProfileFormat.save_cpuprofile(filename, data, period, kwargs...)

    JSONRPC.send(conn_endpoint[], repl_showprofileresult_file_notification_type, filename)
    # JSONRPC.send(conn_endpoint[], repl_showprofileresult_notification_type, filename)
end

"""
    @profview f(args...)

Clear the Profile buffer, profile `f(args...)`, and view the result graphically.
"""
macro profview(ex, args...)
    return quote
        Profile.clear()
        Profile.@profile $(esc(ex))
        view_profile(;$(esc.(args)...))
    end
end


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

name = "CMakeBuilder"
version = v"0.1.0"

# Collection of sources required to build CMakeBuilder
sources = [
    "https://cmake.org/files/v3.15/cmake-3.15.0.tar.gz" =>
    "0678d74a45832cacaea053d85a5685f3ed8352475e6ddf9fcb742ffca00199b5",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
ls
cd cmake-3.15.0/
ls
mkdir build
ls
cd build
ls
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DCMAKE_USE_SYSTEM_LIBRARIES=Off -DBUILD_CursesDialog=Off ..
make
make install
exit

"""

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

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

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

]

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

max_stars_repo_name: jalving/CMakeBuilder

max_issues_repo_name: jalving/CMakeBuilder

max_forks_repo_name: jalving/CMakeBuilder