text:
using AssociativeArrays
using Base.Test

# write your own tests here
@test 1 == 1
# check that we can instantiate AssociativeArrays
d = Dict{NTuple{2,Int}, Float64}()
twod_itof = AssociativeArray(d, (5,5))
@test typeof(AssociativeArray(Dict{NTuple{2,Int},Int}(), (5,6))) == AssociativeArray{Int, 2, Dict{NTuple{2,Int}, Int}}
Dict{NTuple{2,Int}, Int}()

# check that access and assignment work correctly
a = Dict{NTuple{2,Int}, Int}()
a[(3,4)] = 1
@test a[(3,4)] == 1
sa = AssociativeArray(a, (5,5))
@test typeof(sa) == AssociativeArray{Int, 2, Dict{NTuple{2,Int}, Int}}
@test size(sa) == (5,5)
@show getindex(sa, 3,4)
sa[4,4] = -2
@test sa[4,4] == -2
setindex!(sa, -1, 4,4)
@test sa[4,4] == -1
println(sa)
@show fa = similar(sa, Float64, (3,3))
@show fa[3,3] = 1
@show fa

info("special 1d case")
oned_itoi = AssociativeArray(Dict{Int,Int}(), 5)
@show oned_itoi
oned_itoi[4] = 1
@test oned_itoi[4] == 1
@show oned_itoi
@show similar(oned_itoi, Float64, 3)
@show sim1 = similar(oned_itoi, UInt8, 3)
@test length(sim1.data) == 0
a = AssociativeArray{Float64, 2, Dict{Int, Float64}}(Dict{Int, Float64}(), (3,4));

@test_throws MethodError AssociativeArray(Dict{Int, Int}(), (3,4), Float64)
dat = Dict{NTuple{2,Int}, Int}()
# AssociativeArray(dat, (3,4), Int)


text:
# Examples:
#     julia -i activate.jl             -- from command line
#     include("PATH/TO/activate.jl")   -- in Julia script

import Pkg
Pkg.activate(@__DIR__)
import PyVirtualenv
PyVirtualenv.activate_pipenv(@__DIR__)

ext: jl

max_stars_repo_name: tkf/PyVirtualenvExample

max_issues_repo_name: tkf/PyVirtualenvExample

max_forks_repo_name: tkf/PyVirtualenvExample


text:
using MVDA, Plots, Statistics, Random, LinearAlgebra, LaTeXStrings, StableRNGs

font_main = Plots.font("Computer Modern", 12)
font_small = Plots.font("Computer Modern", 10)
w, h = default(:size)
upscale = 1.5
default(
    markerstrokewidth=0.5,
    titlefont=font_main,
    guidefont=font_main,
    tickfont=font_small,
    legendfont=font_small,
    size=(upscale*w, upscale*h),
    dpi=200
)

function draw_circle!(v, ϵ)
    xs = v[1] .+ range(-ϵ, ϵ, length=201)
    ys1 = sqrt.(ϵ^2 .- (xs .- v[1]) .^2) .+ v[2]
    ys2 = -sqrt.(ϵ^2 .- (xs .- v[1]) .^2) .+ v[2]
    plot!(xs, ys1, color=:black, style=:dash, label="")
    plot!(xs, ys2, color=:black, style=:dash, label="")
end

function run(dir, λs)
    # problem dimensions
    n = 300
    p = 1
    c = 3

    # simulate 3 classes with single predictor as in Wu, Lange 2010
    rng = StableRNG(1234)
    μ = Float64[-4,0,4]
    targets = [ones(Int, 100); 2*ones(Int, 100); 3*ones(Int,100)]
    X = zeros(n, p)
    for j in 1:c
        X[100*(j-1)+1:100*j,:] .= μ[j] .+ randn(rng, 100)
    end

    # visualize the samples
    ccp = palette(:seaborn_colorblind)
    colors = [ccp[1] ccp[2] ccp[3]]
    ytmp = [-0.2*ones(100); 0*ones(100); 0.2*ones(100)]
    fig1 = scatter(X, ytmp; group=targets, ylims=(-0.5,0.5), yticks=nothing, palette=ccp, legend=:outerright, xlabel=L"x_{1}", title="True Classes vs Feature Vectors")

    # create a problem instance
    problem = MVDAProblem(targets, X, intercept=true)

    # set hyperparameters
    ϵ = MVDA.maximal_deadzone(problem)

    for λ in λs
        # compute an initial solution w/ regularization
        @time MVDA.init!(SD(), problem, ϵ, λ, verbose=true, maxiter=10^4, nesterov_threshold=0, gtol=1e-6)

        Y = problem.Y
        B = problem.coeff.all
        V = hcat([problem.label2vertex[j] for j in 1:c]...)
        copyto!(problem.proj.all, B)

        # assign vertices and labels
        Yfit = problem.X*B
        labels = MVDA.classify(problem, X)

        # assess misclassification error
        train_error = round(100 * (1 - sum(labels .== targets)/n), sigdigits=4)

        # get distances and plot with assigned labels
        dist = zeros(n, c)
        fig2 = plot(xlabel=L"x_{1}", ylabel="Distance to Class Vertex", ylims=(0.0, 2.25), title="Misclassification Error: $(train_error) %", legend=nothing)
        for j in 1:c
            dist[:,j] .= [ norm(Yfit[i,:] - V[:,j]) for i in 1:n]
            scatter!(fig2, X, dist[:,j], label=j, palette=ccp, color=colors[j])
        end
        hline!(fig2, [ϵ], label="", lw=3, style=:dash, color=:black)

        # visualize data in vertex space
        vertex_labels = map(v -> problem.vertex2label[v], eachcol(V))
        fig3 = scatter(V[1,:], V[2,:], xlabel=L"y_{1}", ylabel=L"y_{2}", title="Classification in Vertex Space", group=vertex_labels, palette=ccp, legend=:bottomright, marker=:star, markersize=upscale*8, aspect_ratio=1)
        scatter!(Yfit[:,1], Yfit[:,2], label=nothing, palette=ccp, color=colors[labels])
        foreach(v -> draw_circle!(v, ϵ), problem.vertex)

        layout = @layout [[a{0.2h}; b{0.8h}] c{0.6w}]
        fig = plot(fig1, fig2, fig3, layout=layout, size=(2*w,1.1*h), margins=5Plots.mm)

        filename = joinpath(dir, "reg-lambda=$(λ).png")
        @info "Saving figure to $(filename)"
        savefig(fig, filename)
    end
end

# Parse input arguments and run.
length(ARGS) < 2 && error("Script requires a directory and at least one input λ > 0 for regularization.")
dir = ARGS[1]
λs = map(arg -> parse(Float64, arg), ARGS[2:end])
@info "Output directory: $(dir)"
@info "λ values: $(λs)"
run(dir, λs)


text:
sparsity(B) = Float32(nnz(B)/length(B))
sparsity(B::Matrix) = 1.0f0
sparsity(B::Adjoint) = sparsity(B.parent)

function preconditioner(B, s, tol=1e-2)
    upsample(droptol!(sparse(inv(downsample2(B, s)))/s/s, tol), s)
end

#downsample(x, s = 10) = [mean(x[i:i+s-1,j:j+s-1]) for i in 1:s:size(x,1)-s, j in 1:s:size(x, 2)-s]
downsample(x, s = 10) = [mean(x[i:i+s-1,j:j+s-1]) for i in 1:s:size(x,1)-s+1, j in 1:s:size(x, 2)-s+1]
upsample(x, s = 10) = repeat(x, inner=(s, s))


text:
rules = Dict(
  [match(r"(\w+ \w+) bags contain", line).captures[1] => Dict(
    [m.captures[2] => parse(Int, m.captures[1])
      for m in eachmatch(r"(\d+) (\w+ \w+)", line)
    ])
    for line in eachline("a7.txt")
  ])

function addcontainersof!(containers, bag)
  for (colour, contents) in rules
    if !in(containers, colour) && haskey(contents, bag)
      push!(containers, colour)
      addcontainersof!(containers, colour)
    end
  end
  containers
end

length(addcontainersof!(Set(), "shiny gold")) |> println

function contains(bag, multiplier)::Int
  multiplier * (1 + sum([contains(colour, count) for (colour, count) in rules[bag]]))
end

contains("shiny gold", 1) - 1 |> println


text:
using KernelAbstractions
using Oceananigans.Utils
using Oceananigans.Architectures: device, @hascuda, CPU, GPU, array_type
using Oceananigans.Fields: datatuple
using Oceananigans.BoundaryConditions: apply_x_bcs!, apply_y_bcs!, apply_z_bcs!

@kernel function update_total_density!(total_density, grid, gases, tracers)
    i, j, k = @index(Global, NTuple)

    @inbounds total_density[i, j, k] = diagnose_density(i, j, k, grid, gases, tracers)
end

# This is for users. Do not use in time_stepping.jl as it doesn't make use of intermediate fields.
function update_total_density!(model)

    density_update_event =
        launch!(model.architecture, model.grid, :xyz, update_total_density!,
                datatuple(model.total_density), model.grid, model.gases, datatuple(model.tracers),
                dependencies=Event(device(model.architecture)))

    wait(device(model.architecture), density_update_event)

    return nothing
end

#####
##### Computing slow source terms (viscous dissipation, diffusion, and Coriolis terms).
#####

function compute_slow_source_terms!(slow_source_terms, arch, grid, thermodynamic_variable, gases, gravity, coriolis, closure, total_density, momenta, tracers, diffusivities, forcing, clock)

    slow_source_terms, total_density, momenta, tracers, diffusivities =
        datatuples(slow_source_terms, total_density, momenta, tracers, diffusivities)

    workgroup, worksize = work_layout(grid, :xyz)
    barrier = Event(device(arch))

    momentum_kernel! = compute_slow_momentum_source_terms!(device(arch), workgroup, worksize)
    tracer_kernel! = compute_slow_tracer_source_terms!(device(arch), workgroup, worksize)
    thermodynamic_variable_kernel! = compute_slow_thermodynamic_variable_source_terms!(device(arch), workgroup, worksize)

    momentum_event = momentum_kernel!(slow_source_terms, grid, coriolis, closure, total_density, momenta, tracers, diffusivities, forcing, clock, dependencies=barrier)

    events = [momentum_event]

    for (tracer_index, ρc_name) in enumerate(propertynames(tracers))
        ρc   = getproperty(tracers, ρc_name)
        S_ρc = getproperty(slow_source_terms.tracers, ρc_name)
        forcing_ρc = getproperty(forcing, ρc_name)

        tracer_event = tracer_kernel!(S_ρc, grid, closure, tracer_index, total_density, ρc, momenta, tracers, diffusivities, forcing_ρc, clock, dependencies=barrier)
        push!(events, tracer_event)
    end

    thermodynamic_variable_event = thermodynamic_variable_kernel!(slow_source_terms.tracers[1], grid, thermodynamic_variable, gases, gravity, closure, total_density, momenta, tracers, diffusivities, dependencies=barrier)
    push!(events, thermodynamic_variable_event)

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end

@kernel function compute_slow_momentum_source_terms!(slow_source_terms, grid, coriolis, closure, total_density, momenta, tracers, diffusivities, forcing, clock)
    i, j, k = @index(Global, NTuple)

    @inbounds slow_source_terms.ρu[i, j, k] = ρu_slow_source_term(i, j, k, grid, coriolis, closure, total_density, momenta, diffusivities) + forcing.ρu(i, j, k, grid, clock, merge(momenta, tracers))
    @inbounds slow_source_terms.ρv[i, j, k] = ρv_slow_source_term(i, j, k, grid, coriolis, closure, total_density, momenta, diffusivities) + forcing.ρv(i, j, k, grid, clock, merge(momenta, tracers))
    @inbounds slow_source_terms.ρw[i, j, k] = ρw_slow_source_term(i, j, k, grid, coriolis, closure, total_density, momenta, diffusivities) + forcing.ρw(i, j, k, grid, clock, merge(momenta, tracers))
end

@kernel function compute_slow_tracer_source_terms!(S_ρc, grid, closure, tracer_index, total_density, ρc, momenta, tracers, diffusivities, forcing, clock)
    i, j, k = @index(Global, NTuple)

    @inbounds S_ρc[i, j, k] = ρc_slow_source_term(i, j, k, grid, closure, tracer_index, total_density, ρc, diffusivities) + forcing(i, j, k, grid, clock, merge(momenta, tracers))
end

@kernel function compute_slow_thermodynamic_variable_source_terms!(S_ρt, grid, thermodynamic_variable, gases, gravity, closure, total_density, momenta, tracers, diffusivities)
    i, j, k = @index(Global, NTuple)

    @inbounds S_ρt[i, j, k] += ρt_slow_source_term(i, j, k, grid, closure, thermodynamic_variable, gases, gravity, total_density, momenta, tracers, diffusivities)
end

#####
##### Computing fast source terms (advection, pressure gradient, and buoyancy terms).
#####

function compute_fast_source_terms!(fast_source_terms, arch, grid, thermodynamic_variable, gases, gravity, advection_scheme, total_density, momenta, tracers, slow_source_terms)

    fast_source_terms, total_density, momenta, tracers, slow_source_terms =
        datatuples(fast_source_terms, total_density, momenta, tracers, slow_source_terms)

    workgroup, worksize = work_layout(grid, :xyz)
    barrier = Event(device(arch))

    momentum_kernel! = compute_fast_momentum_source_terms!(device(arch), workgroup, worksize)
    tracer_kernel! = compute_fast_tracer_source_terms!(device(arch), workgroup, worksize)
    thermodynamic_variable_kernel! = compute_fast_thermodynamic_variable_source_terms!(device(arch), workgroup, worksize)

    momentum_event = momentum_kernel!(fast_source_terms, grid, thermodynamic_variable, gases, gravity, advection_scheme, total_density, momenta, tracers, slow_source_terms, dependencies=barrier)

    events = [momentum_event]

    for ρc_name in propertynames(tracers)
        ρc   = getproperty(tracers, ρc_name)
        F_ρc = getproperty(fast_source_terms.tracers, ρc_name)
        S_ρc = getproperty(slow_source_terms.tracers, ρc_name)

        tracer_event = tracer_kernel!(F_ρc, grid, advection_scheme, total_density, momenta, ρc, S_ρc, dependencies=barrier)
        push!(events, tracer_event)
    end

    thermodynamic_variable_event = thermodynamic_variable_kernel!(fast_source_terms.tracers[1], grid, thermodynamic_variable, gases, gravity, total_density, momenta, tracers, dependencies=barrier)
    push!(events, thermodynamic_variable_event)

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end

@kernel function compute_fast_momentum_source_terms!(fast_source_terms, grid, thermodynamic_variable, gases, gravity, advection_scheme, total_density, momenta, tracers, slow_source_terms)
    i, j, k = @index(Global, NTuple)

    @inbounds fast_source_terms.ρu[i, j, k] = ρu_fast_source_term(i, j, k, grid, thermodynamic_variable, gases, gravity, advection_scheme, total_density, momenta, tracers, slow_source_terms.ρu)
    @inbounds fast_source_terms.ρv[i, j, k] = ρv_fast_source_term(i, j, k, grid, thermodynamic_variable, gases, gravity, advection_scheme, total_density, momenta, tracers, slow_source_terms.ρv)
    @inbounds fast_source_terms.ρw[i, j, k] = ρw_fast_source_term(i, j, k, grid, thermodynamic_variable, gases, gravity, advection_scheme, total_density, momenta, tracers, slow_source_terms.ρw)
end

@kernel function compute_fast_tracer_source_terms!(F_ρc, grid, advection_scheme, total_density, momenta, ρc, S_ρc)
    i, j, k = @index(Global, NTuple)

    @inbounds F_ρc[i, j, k] = ρc_fast_source_term(i, j, k, grid, advection_scheme, total_density, momenta, ρc, S_ρc)
end

@kernel function compute_fast_thermodynamic_variable_source_terms!(F_ρt, grid, thermodynamic_variable, gases, gravity, total_density, momenta, tracers)
    i, j, k = @index(Global, NTuple)

    @inbounds F_ρt[i, j, k] += ρt_fast_source_term(i, j, k, grid, thermodynamic_variable, gases, gravity, total_density, momenta, tracers)
end

#####
##### Calculating boundary tendency contributions
#####

function calculate_boundary_tendency_contributions!(source_terms, arch, momenta, tracers, clock, model_fields)

    barrier = Event(device(arch))

    events = []

    # Momentum fields
    momentum_source_terms = (source_terms.ρu, source_terms.ρv, source_terms.ρw)

    for (ρϕ_source_term, ρϕ) in zip(momentum_source_terms, momenta)
        x_bcs_event = apply_x_bcs!(ρϕ_source_term, ρϕ, arch, barrier, clock, model_fields)
        y_bcs_event = apply_y_bcs!(ρϕ_source_term, ρϕ, arch, barrier, clock, model_fields)
        z_bcs_event = apply_z_bcs!(ρϕ_source_term, ρϕ, arch, barrier, clock, model_fields)

        push!(events, x_bcs_event, y_bcs_event, z_bcs_event)
    end

    # Tracer fields
    for (ρϕ_source_term, ρϕ) in zip(source_terms.tracers, tracers)
        x_bcs_event = apply_x_bcs!(ρϕ_source_term, ρϕ, arch, barrier, clock, model_fields)
        y_bcs_event = apply_y_bcs!(ρϕ_source_term, ρϕ, arch, barrier, clock, model_fields)
        z_bcs_event = apply_z_bcs!(ρϕ_source_term, ρϕ, arch, barrier, clock, model_fields)

        push!(events, x_bcs_event, y_bcs_event, z_bcs_event)
    end

    events = filter(e -> typeof(e) <: Event, events)

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end

#####
##### Advancing state variables
#####

function advance_state_variables!(state_variables, arch, grid, momenta, tracers, fast_source_terms; Δt)

    state_variables, momenta, tracers, fast_source_terms =
        datatuples(state_variables, momenta, tracers, fast_source_terms )

    workgroup, worksize = work_layout(grid, :xyz)
    barrier = Event(device(arch))

    momentum_kernel! = advance_momentum!(device(arch), workgroup, worksize)
    tracer_kernel! = advance_tracer!(device(arch), workgroup, worksize)

    momentum_event = momentum_kernel!(state_variables, grid, momenta, fast_source_terms, Δt, dependencies=barrier)

    events = [momentum_event]

    for ρc_name in propertynames(tracers)
        ρc   = getproperty(tracers, ρc_name)
        ρc⁺  = getproperty(state_variables.tracers, ρc_name)
        F_ρc = getproperty(fast_source_terms.tracers, ρc_name)

        tracer_event = tracer_kernel!(ρc⁺, grid, ρc, F_ρc, Δt, dependencies=barrier)
        push!(events, tracer_event)
    end

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end

@kernel function advance_momentum!(momenta⁺, grid, momenta, fast_source_terms, Δt)
    i, j, k = @index(Global, NTuple)

    @inbounds momenta⁺.ρu[i, j, k] = momenta.ρu[i, j, k] + Δt * fast_source_terms.ρu[i, j, k]
    @inbounds momenta⁺.ρv[i, j, k] = momenta.ρv[i, j, k] + Δt * fast_source_terms.ρv[i, j, k]
    @inbounds momenta⁺.ρw[i, j, k] = momenta.ρw[i, j, k] + Δt * fast_source_terms.ρw[i, j, k]
end

@kernel function advance_tracer!(ρc⁺, grid, ρc, F_ρc, Δt)
    i, j, k = @index(Global, NTuple)

    @inbounds ρc⁺[i, j, k] = ρc[i, j, k] + Δt * F_ρc[i, j, k]
end


text:
using Pkg
using Plots
using DataFrames
pyplot()

### functions for simulations
function calc_new_state(y::AbstractArray,
    a::Dict,
    b::Dict,
    l::Dict;
    from::Float64=0.01)
    if length(y)==0
        y = [from]
    end
    y[end]+(1-y[end])*(a["walk"]*(1-exp(-b["walk"]*l["walk"]*a["walk"]*y[end]))+a["TTC"]*(1-exp(-b["TTC"]*l["TTC"]*a["TTC"]*y[end])))
end

function simul(s::Dict,
               b::Dict,
               l::Dict;
               to_perc::Bool=true,
               from::Float64=0.01,
               to::Union{Float64, Int}=0.99)

               # initial calculations
                e1 = (exp(-(s["walk"])/min(s["walk"], s["TTC"])))
                e2 = (exp(-(s["TTC"])/min(s["walk"], s["TTC"])))
                a = Dict("walk" => e1/sum([e1,e2]), "TTC" => e2/sum([e1,e2]))
                # simulation loop
                y = []
                append!(y, from)
                if to_perc
                    while y[end] <= to
                        append!(y, calc_new_state(y,a,b,l))
                    end
                else
                    for i in 1:to
                        append!(y, calc_new_state(y,a,b,l))
                    end
                end
                y
end

#####################################
### scenario 1:
s = Dict("walk" => 15, "TTC" => 10) # time of travel
l = Dict("walk" => 0.01, "TTC" => 0.02) # contagion
b = Dict("walk" => 5, "TTC" => 5) # interaction frequency

x1 = simul(s,b,l)
length(x1) - 1

### scenario 2:
s = Dict("walk" => 15, "TTC" => 10) # time of travel
l = Dict("walk" => 0.01, "TTC" => 0.02) # contagion
b = Dict("walk" => 3, "TTC" => 5) # interaction frequency

x2 = simul(s,b,l)
length(x2) - 1

# PLOT1
Plots.plot(1:length(x1), x1, linecolor="purple", lw = 2, legend=false, grid=false)
ylabel!("% infected")
xlabel!("days")
savefig("plot1.eps")

# PLOT2
Plots.plot(1:length(x2), x2, linecolor="purple", lw = 2, legend=false, grid=false)
ylabel!("% infected")
xlabel!("days")
savefig("plot2.eps")

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

#PLOT3
Plots.plot(1:length(x1), x1, linecolor="purple", label="Scenario 1", lw=2, grid = false, legend=:topleft)
Plots.plot!(1:length(x2), x2, linecolor="black", label="Scenario 2", lw=2 #linestyle = :dot, )
)
ylabel!("% infected")
xlabel!("days")
savefig("plot3.eps")

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

# PLOT 4
s = Dict("walk" => 15, "TTC" => 10) # time of travel
l = Dict("walk" => 0.01, "TTC" => 0.02) # contagion
b = Dict("walk" => 3, "TTC" => 5) # interaction frequency

results_full = Dict()
results_short = DataFrame()
for i in 1:50
    s["TTC"] = i
    results_full[i] = simul(s,b,l)
    append!(results_short, DataFrame(s_ttc = [i], days = [length(results_full[i])-1]))
end

results_short # ok

Plots.plot(results_short[!,:s_ttc], results_short[!,:days], linecolor="purple", lw = 2, legend=false, grid=false)
ylabel!("days")
xlabel!("time when taking TTC")
savefig("plot4.eps")

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

# PLOT 5
s = Dict("walk" => 15, "TTC" => 10) # time of travel
l = Dict("walk" => 0.01, "TTC" => 0.02) # contagion
b = Dict("walk" => 3, "TTC" => 5) # interaction frequency

results_full2 = Dict()
results_short2 = DataFrame()
for i in 1:50
    s["TTC"] = i
    results_full2[i] = simul(s,b,l, to_perc=false, to=101)
    append!(results_short2, DataFrame(s_ttc = [i], infected = results_full2[i][101]))
end

results_short2 # ok

Plots.plot(results_short2[!,:s_ttc], results_short2[!,:infected], linecolor="purple", lw = 2, legend=false, grid=false)
ylabel!("% infected after 100 days")
xlabel!("time when taking TTC")
savefig("plot5.eps")

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
function run(s)
    # Your code here
    entries = split(s, "\n")
    entries = map(x -> parse(Int, x), entries)
    for (idx1, el1) in enumerate(entries)
        for (idx2, el2) in enumerate(entries[idx1+1:end])
            if el1 + el2 > 2020
                continue
            end
            for el3 in entries[idx1+idx2+1:end]
                if el1 + el2 + el3 == 2020
                    return el1 * el2 * el3
                end
            end
        end
    end
    return 0
end

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

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time*1000)")
    println(res)
end

main()


text:
using SparseArrays
using LinearAlgebra

@testset "utility" begin
    @testset "random_edge" begin
        @test_throws ArgumentError random_edge(empty_graph(0))
        @test_throws ArgumentError random_edge(empty_graph(5))
        A = sparse([1],[2],1.0,5,5)
        @test random_edge(A)[1:2] == (1,2)
        @test random_edge(A)[1:2] == (1,2)
        @test random_edge(A)[1:2] == (1,2)
        @test random_edge(A)[1:2] == (1,2)
        G = lollipop_graph(5,3)
        A = sparse(G)
        # there are 16 edges, to 100*16*log(16) should randomly generate all of them
        M = zeros(8,8)
        Random.seed!(1) # make it determinstic-ish
        ntrials = 100*16*4
        for i=1:ntrials # log16 = 4 in base 2
            M[random_edge(G)[1:2]...] += 1
        end
        @test sum(M.*A) == ntrials          # this means we got all
        # @test all(findall(!iszero,A.*(M .+ A)) .>= 1) # this means we got all entries
        @test all(x -> iszero(x) || x >= 1, A.*(M .+ A))
        @test std(M[findall(iszero,M)]) <= 100

        n = 16
        A = sparse(1:n-1,2:n,1,n,n)
        # there are 15 edges, to 100*15*log(15) should randomly generate all of them
        M = zeros(n,n)
        Random.seed!(1) # make it determinstic-ish
        ntrials = 100*16*4 # we just use the same one
        for i=1:ntrials # log16 = 4 in base 2
            M[random_edge(A)[1:2]...] += 1
        end
        @test sum(M.*A) == ntrials          # this means we got all
        # @test all(findall(iszero,A.*(M .+ A)) .>= 1) # this means we got all entries
        @test all(x -> iszero(x) || x >= 1, A.*(M .+ A))
        @test std(M[findall(iszero,M)]) <= 100
    end
    @testset "undirected_edges" begin
        @test map(length, undirected_edges(empty_graph(0))) == (0,0)
        @test map(length, undirected_edges(empty_graph(1))) == (0,0)
        @test map(length, undirected_edges(lollipop_graph(5,1))) == (5,5)
        @test undirected_edges(lollipop_graph(5,1)) == ([1,2,3,4,5],[2,3,4,5,6])
        A = sparse([2],[1],1.0,5,5)
        @test map(length, undirected_edges(MatrixNetwork(A))) == (0,0)
    end
    @testset "directed_edges" begin
        @test map(length, directed_edges(empty_graph(0))) == (0,0)
        @test map(length, directed_edges(empty_graph(1))) == (0,0)
        @test map(length, directed_edges(lollipop_graph(5,1))) == (10,10)
        A = sparse([1],[2],1.0,5,5)
        @test map(length, directed_edges(MatrixNetwork(A))) == (1,1)
        @test directed_edges(MatrixNetwork(A)) == ([1],[2])
    end
end


text:
"""
	struct ExtractBranch
		vec::Dict{String,Any}
		other::Dict{String,Any}
	end

	extracts all items in `vec` and in `other` and return them as a TreeNode.
"""
struct ExtractBranch{S,V} <: AbstractExtractor
	vec::S
	other::V
	function ExtractBranch(v::S,o::V) where {S<:Union{Dict,Nothing},V<:Union{Dict,Nothing}}
		v = (v == nothing || isempty(v)) ? nothing : v
		o = (o == nothing || isempty(o)) ? nothing : o
		new{typeof(v),typeof(o)}(v,o)
	end
end

function Base.getindex(m::ExtractBranch, s::String)
	m.vec != nothing && haskey(m.vec, s) && return(m.vec[s])
	m.other != nothing && haskey(m.other, s) && return(m.other[s])
	nothing
end

replacebyspaces(pad) = map(s -> (s[1], " "^length(s[2])), pad)

function printdict(io, d::Dict, ml, c, pad)
  k = sort(collect(keys(d)))
  for i in 1:length(k)-1
	  s = "  ├──"*"─"^(ml-length(k[i]))*" "
		  paddedprint(io, s, color=c, pad=pad)
		  show(io, d[k[i]], pad=[pad; (c, "  │"*" "^(ml-length(k[i])+2))], key = k[i])
  end
  s = "  └──"*"─"^(ml-length(k[end]))*" "
  paddedprint(io, s, color=c, pad=pad)
  show(io, d[k[end]], pad=[pad; (c, " "^(ml-length(k[end])+4))], key = k[end])
end

function Base.show(io::IO, m::ExtractBranch; pad = [], key::String="")
  c = COLORS[(length(pad)%length(COLORS))+1]
  ml = m.vec   != nothing ? maximum(length(k) for k in keys(m.vec)) : 0
  ml = m.other != nothing ? max(ml, maximum(length(k) for k in keys(m.other))) : ml
  key *=": "
  paddedprint(io,"$(key)struct\n", color = c)
	if isnothing(m.vec)
		paddedprint(io, "  Empty vec\n", color = c, pad=pad)
	else
		paddedprint(io, "  Vec:\n", color = c, pad=pad)
		printdict(io, m.vec, ml, c, pad)
	end
	if isnothing(m.other)
		paddedprint(io, "  Empty other\n", color = c, pad=pad)
	else
		paddedprint(io, "  Other:\n", color = c, pad=pad)
		printdict(io, m.other, ml, c, pad)
	end
end

extractsmatrix(s::ExtractBranch) = false
(s::ExtractBranch)(v::V) where {V<:Nothing} = s(Dict{String,Any}())

# function (s::ExtractBranch{S,V})(v::Dict) where {S<:Dict,V<:Dict}
# 	x = vcat([f(get(v,k,nothing)) for (k,f) in s.vec]...)
# 	o = [f(get(v,k,nothing)) for (k,f) in s.other]
# 	data = tuple(x, o...)
# 	TreeNode(data)
# end

function (s::ExtractBranch{S,V})(v::Dict) where {S<:Dict,V<:Dict}
	x = vcat([f(get(v,k,nothing)) for (k,f) in s.vec]...)
	o = [Symbol(k) => f(get(v,k,nothing)) for (k,f) in s.other]
	# o = [f(get(v,k,nothing)) for (k,f) in s.other]
	data = (; :scalars => x,o...)
	# data = tuple(x, o...)
	TreeNode(data)
end

(s::ExtractBranch{S,V})(v::Dict) where {S<:Dict,V<:Nothing} = vcat([f(get(v,k,nothing)) for (k,f) in s.vec]...)

function (s::ExtractBranch{S,V})(v::Dict) where {S<:Nothing,V<:Dict}
	# o = [f(get(v,k,nothing)) for (k,f) in s.other]
	o = [Symbol(k) => f(get(v,k,nothing)) for (k,f) in s.other]
	if length(o) == 1
		# return(o[1])
		return(o[1].second)
	else
		return(TreeNode((;o...)))
	end
end

extractbatch(extractor, samples) = reduce(catobs, map(s-> extractor(s), samples))

Base.hash(e::ExtractBranch, h::UInt) = hash((e.vec, e.other), h)
Base.:(==)(e1::ExtractBranch, e2::ExtractBranch) = e1.vec == e2.vec && e1.other == e2.other


text:
reduce(+, map(x->x^2, 1:100)) #338350
reduce(+, 1:100)^2 #25502500

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
simple_auth = TM1.UsernamePassAuth("admin", "apple");
@test simple_auth.username == "admin"
@test simple_auth.password == "apple"


text:
struct Squares
    count::Int
end

Base.iterate(S::Squares, state=1) =
    state > S.count ? nothing : (state*state, state+1)

for item in Squares(9) # or "for item = Squares(9)"
    # body
    display(item)
end

println(25 in Squares(10))

Base.eltype(::Type{Squares}) = Int # Note that this is defined for the type
Base.length(S::Squares) = S.count
Base.sum(S::Squares) = (n = S.count; return n*(n+1)*(2n+1)÷6)

collect(Squares(4)) |> display


text:
# Parse Primals and Duals

struct primal_dual
    primals
    duals
end

function parsePrimalDualVec(pdVec, primalSize::Int64)
    return primal_dual(pdVec[1:primalSize], pdVec[primalSize + 1:end])
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
#=-------------splitstep.jl----------------------------------------------------#
#
#              split step method for macroscopic wave function generation
#
# Purpose: This file is intended to create a splitstep method for BEC systems.
#          Eventually we may incorporate real-time evolution, but we need to
#          first generate the grounds state and go from there.
#
#   Notes: Units of hbar = c  = 1
#          Initial gaussian guess for wavefunction not normalized because
#              abs2(wave) is described by normalization condition
#-----------------------------------------------------------------------------=#

global xmax = 40
global res = 2^8
global g = 500

# This initializes the potential well and also the gaussian init wavefunction
# All the normalization for the initial gaussian guess is commented out
function initialize()
    gaus = zeros(res)
    pot = zeros(res)
    #gausfile = open("check_gaus.dat", "w")

    for ix = 1:res
        x = ix * (xmax/res) - xmax / 2
        gaus[ix] = exp(-x * x )
        pot[ix] = x * x
        #println(gausfile, gaus[ix])
    end

    return gaus, pot
end

# This will return the kinetic and potential energies at every timestep
# Note unphysical units will be fixed in the future.
function energy(wave, pot, dt)

    KE = complex(zeros(size(wave,1)))
    PE = complex(zeros(size(wave,1)))

    dk = 2pi / xmax

    # note that phase space ordering is piecemeal / different than real.
    for i = 1:size(wave,1)
        if i <= size(wave,1) / 2
            k = dk * (i)
        else
            k = dk * (i - (size(wave,1) / 2))
        end
        KE[i] = exp( -0.5 * (k*k) * dt)
        PE[i] = exp( -0.5 * (pot[i] + g*abs2(wave[i])) *dt)
    end
    return PE, KE
end

# The whole shebang
# Psi * Uv -> fft -> Psi * Uk -> ifft -> renormalization -> cycle
function splitstep(stepnum, dt)

    # initialization
    output = open("out.dat", "w")
    wave, pot = initialize()
    density = zeros(size(wave,1))

    # looping through space (both normal and phase)
    for j = 1:stepnum

        # finding PE and KE terms
        PE, KE = energy(wave, pot, dt)
        norm_const = 0

        for i = 1:res
            wave[i] *= PE[i]
        end

        wave = fft(wave)

        for i = 1:res
            wave[i] *= KE[i]
        end

        wave = ifft(wave)

        # normalizing is necessary to keep wavefunction
        norm_const = 0
        for i = 1:res
            norm_const += abs2(wave[i]) * xmax / res
        end

        for i = 1:res
            wave[i] *= 1/norm_const
        end

        # density is Psi^2
        for i = 1:res
            density[i] = abs2(wave[i])
        end

        # Output
        if j % 10 == 0 || j == 1
            for i = 1:res
                println(output, density[i])
                #println(PE[i], '\t', KE[i])
            end

            print(output, '\n', '\n')
        end

    end
end

# Main

splitstep(10000, 0.0001)

max_issues_repo_name: shiffman/simuleios

max_forks_repo_name: shiffman/simuleios


text:
abstract Geometry

immutable Conductor <: Geometry
    name::String
    outperm::Real
    tran::Tuple{Real,Real,Real}
    shape::Shape
end

function Conductor(name::String,outperm::Real,shape::Shape)
    return Conductor(name,outperm,(0,0,0),shape)
end

immutable Dielectric <: Geometry
    name::String
    outperm::Real
    inperm::Real
    tran::Tuple{Real,Real,Real}
    ref::Tuple{Real,Real,Real}
    shape::Shape
end

function Dielectric(name::String,outperm::Real,inperm::Real,ref::Tuple{Real,Real,Real},shape::Shape)
    return Dielectric(name,outperm,inperm,(0,0,0),ref,shape)
end

immutable ConductorGroup <: Geometry
    name::String
    tran::Tuple{Real,Real,Real}
    geoms::Vector{Conductor}
end

function ConductorGroup(name::String,geoms::Vector{Conductor})
    return ConductorGroup(name,(0,0,0),geoms)
end

immutable DielectricGroup <: Geometry
    name::String
    tran::Tuple{Real,Real,Real}
    geoms::Vector{Dielectric}
end

function DielectricGroup(name::String,geoms::Vector{Conductor})
    return DielectricGroup(name,(0,0,0),geoms)
end

move(tran::Tuple{Real,Real,Real},g::Conductor) = Conductor(g.name,g.outperm,g.tran+tran,g.shape)
move(tran::Tuple{Real,Real,Real},g::Dielectric) = Dielectric(g.name,g.outperm,g.inperm,g.tran+tran,g.ref+tran,g.shape)
move(tran::Tuple{Real,Real,Real},g::ConductorGroup) = ConductorGroup(g.name,g.tran+tran,g.geoms)
move(tran::Tuple{Real,Real,Real},g::DielectricGroup) = DielectricGroup(g.name,g.tran+tran,g.geoms)

function regroup(name::String,tran,groups::Vector{ConductorGroup})
    # Return a new conductor group with the given name
    return ConductorGroup(name,tran,geometries(groups))
end

function regroup(name::String,tran,groups::Vector{DielectricGroup})
    # Return a new dielectric group with the given name
    return ConductorGroup(name,tran,geometries(groups))
end

function geometries{T<:Union{ConductorGroup,DielectricGroup}}(groups::Vector{T})
    # Extract the geometries from each group applying the group
    # translation directly to the extracted geometry
    return mapreduce(group->mapreduce(g->move(group.tran,g),vcat,group.geoms),vcat,groups)
end

function shapes{T<:Geometry}(geoms::Vector{T})
    # Extract all the shapes from a list of geometries
    return mapreduce(shapes,vcat,geoms)
end

function shapes(g::Union{Conductor,Dielectric})
    # Extract shape from geometry
    return g.shape
end

function shapes(g::Union{ConductorGroup,DielectricGroup})
    # Extract shapes from geometry group
    return [e.shape for e in g.geoms]
end

function listfull{T<:Geometry}(geoms::Vector{T})
    listfile   = list(geoms)
    shapefiles = hierfile.(unique(shapes(geoms)))
    return join([listfile;"End";shapefiles],NL)
end

hierfile(shape::Shape) = "File "*shape.name*NL*getfile(shape)*NL*"End"

function list{T<:Geometry}(geoms::Vector{T})
    return join(listline.(geoms),NL)
end

prtu(tu) = @sprintf(" % .5E % .5E % .5E",tu...)

function listline(g::Conductor)
    lg = "G $(g.name)"
    lf = "C $(g.shape.name) $(g.outperm)"*prtu(g.tran)
    return join([lg,lf],NL)
end

function listline(g::Dielectric)
    lg = "G $(g.name)"
    lf = "D $(g.shape.name) $(g.outperm) $(g.inperm)"*prtu(g.tran)*prtu(g.ref)*" - 0x7a49a5ff"
    return join([lg,lf],NL)
end

function listline(g::ConductorGroup)
    lg = "G $(g.name)"
    lf = ["C $(c.shape.name) $(c.outperm)"*prtu(c.tran+g.tran) for c in g.geoms]
    return join([lg,join(lf," +$NL")],NL)
end

function listline(g::DielectricGroup)
    lg = "G $(g.name)"
    lf = ["D $(d.shape.name) $(d.outperm) $(d.inperm)"*prtu(d.tran+g.tran)*prtu(d.ref+g.tran)*" - 0x7a49a5ff" for d in g.geoms]
    return join([lg,join(lf,NL)],NL)
end

import Base.+
+(a::Tuple{Real,Real,Real},b::Tuple{Real,Real,Real}) = (a[1]+b[1],a[2]+b[2],a[3]+b[3])

max_stars_repo_name: jebej/FastCap

max_issues_repo_name: jebej/FastCap

max_forks_repo_name: jebej/FastCap


text:
# TODO: Document this. How is this different from SOCElemConstraint? Why do we need both. How does
# conic form work for SOC constraints.
struct SOCConstraint <: Constraint
    head::Symbol
    id_hash::UInt64
    children::Tuple
    dual::ValueOrNothing

    function SOCConstraint(args::AbstractExpr...)
        children = tuple(args...)
        id_hash = hash((children, :soc))
        return new(:soc, id_hash, children, nothing)
    end
end

function conic_form!(c::SOCConstraint, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, c)
        objectives = Vector{ConicObj}(undef, length(c.children))
        @inbounds for iobj = 1:length(c.children)
            objectives[iobj] = conic_form!(c.children[iobj], unique_conic_forms)
        end
        cache_conic_form!(unique_conic_forms, c, ConicConstr(objectives, :SOC, [length(x) for x in c.children]))
    end
    return get_conic_form(unique_conic_forms, c)
end

# For debugging created this constraint
socp(args::AbstractExpr...) = SOCConstraint(args::AbstractExpr...)

struct SOCElemConstraint <: Constraint
    head::Symbol
    id_hash::UInt64
    children::Tuple
    dual::ValueOrNothing

    function SOCElemConstraint(args::AbstractExpr...)
        children = tuple(args...)
        id_hash = hash((children, :soc_elem))
        return new(:soc_elem, id_hash, children, nothing)
    end
end

function conic_form!(c::SOCElemConstraint, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, c)
        num_constrs = length(c.children[1])
        num_children = length(c.children)
        conic_constrs = Vector{ConicConstr}(undef, num_constrs)
        objectives = Vector{ConicObj}(undef, num_children)
        @inbounds for iobj = 1:num_children
            objectives[iobj] = conic_form!(c.children[iobj], unique_conic_forms)
        end
        @inbounds for row = 1:num_constrs
            objectives_by_row = [get_row(obj, row) for obj in objectives]
            conic_constrs[row] = ConicConstr(objectives_by_row, :SOC, [1, 1, 1])
        end
        cache_conic_form!(unique_conic_forms, c, conic_constrs)
    end
    return get_conic_form(unique_conic_forms, c)
end


text:
type Convex <: AbstractDataEnvelopment
end

type FreeDisposal <: AbstractDataEnvelopment
end

# Directional distance function
immutable DDF{S<:AbstractDataEnvelopment,T<:RS} <: AbstractDEA{S,T}
	Data::DEAData

	function DDF{S,T}(X::Array,Y::Array,gx::Array,gy::Array) where {S<:AbstractDataEnvelopment,T<:RS}
		new(DEAData(X,Y,gx,gy))
	end
end

function getdata(D::DDF)
	return D.Data
end

# Solve convex DDF program with (Xk,Yk) as evaluation point in the direction of (gxk,gyk)
function (D::DDF{Convex,T})(Xk::Array,Yk::Array,gxk::Array,gyk::Array) where T<:RS
	K = getnrdmu(D.Data)
	N,M = getiodim(D.Data)

	# Set the appropriate RTS constraint depending on T
	RSA,RSb,RSsense = RSconstraint(K,T())

  	f = zeros(K+1)
  	f[1] = -1
  	l = [-Inf;zeros(K)]
  	u = [Inf for i=1:K+1]

  	sense = Array{Char}(N+M)
  	sense[1:N+M] = '>'

	X,Y,gx,gy = D[1:end]

	A = [reshape(-gxk',:,1) -X'; reshape(-gyk',:,1) Y'; RSA]
	b = [-Xk; Yk; RSb]
	b = b[:]

	# Solve linear program
	sol = linprog(f,A,[sense;RSsense],b,l,u,ClpSolver())

	if sol.status == :Optimal
		res = DEAResult(sol.sol[1],sol.sol[1].*gxk[1,:],sol.sol[1].*gyk[1,:],sol.sol[2:end])
	else
		res = DEAResult(-Inf)
  	println("Error: solution status $(sol.status)")
	end

	return res
end

# Solve free disposal DDF program under VRS with (Xk,Yk) as evaluation point in the direction of (gxk,gyk)
#FIXME Check that this is correct!
function (D::DDF{FreeDisposal,VRS})(Xk::Array,Yk::Array,gxk::Array,gyk::Array)
	gxk0 = map(isapprox,gxk,0.0.*ones(size(gxk)))
	gyk0 = map(isapprox,gyk,0.0.*ones(size(gyk)))
	if(any(gxk0) || any(gyk0))
		Q = []
		for i in eachindex(D.Data)
			Xi,Yi,gxi,gyi = D[i]
			if(all(Yi[gyk0] .>= Yk[gyk0]) && all(Xi[gxk0] .<= Xk[gxk0]))
				Q = [Q;i]
			end
		end
	else
		Q = getdata(D)
	end

  	ypos = gyk .> 0.0
  	xpos = gxk .> 0.0
	yneg = gyk .< 0.0
	xneg = gxk .< 0.0
  	alpha = -Inf
	alphapeer = 0
	gamma = zeros(getnrdmu(D.Data))
	beta = zeros(getnrdmu(D.Data))
	alpha = zeros(getnrdmu(D.Data))
	for i in eachindex(Q)
    Xi,Yi,gxi,gyi = D.Data[i]

		if(any(xneg) || any(yneg))
			gamma[i] = maximum([((Yi[yneg]-Yk[yneg])./gyk[yneg]); ((Xk[xneg]-Xi[xneg])./gxk[xneg])])
		else
			gamma[i] = 0.0
		end

		beta[i] = minimum([((Yi[ypos]-Yk[ypos])./gyk[ypos]); ((Xk[xpos]-Xi[xpos])./gxk[xpos])])

		if(beta[i] < gamma[i])
			alpha[i] = 0.0
		else
			alpha[i] = beta[i]
		end

  end
	alphapeer = indmax(alpha)
	alpha = maximum(alpha)

	peers = zeros(getnrdmu(D.Data))
	peers[alphapeer] = 1.0
	return DEAResult(alpha,alpha.*gxk,alpha.*gyk,peers)
end

function (D::DDF{FreeDisposal,T})(Xk::Array,Yk::Array,gxk::Array,gyk::Array) where T<:RS
	if(all(gxk .<= 0.0) && all(gyk .<= 0.0))
		return DEAResult(NaN)
	end

	X,Y,gx,gy = D[1:end]

	#Find zero entries in gxk,gyk
	gxk0 = map(isapprox,gxk,0.0.*ones(size(gxk)))
	gyk0 = map(isapprox,gyk,0.0.*ones(size(gyk)))

	#Simar and Vanhems(2012) trick using transformed dataset and computing hyperbolic distance function on the transformed data
	#Note: This assumes gxk,gyk >= 0!
	DGR = Hyperbolic{FreeDisposal,T}(exp.(hcat(X[:,.!gxk0]./gxk[:,.!gxk0],X[:,gxk0])),exp.(hcat(Y[:,.!gyk0]./gyk[:,.!gyk0],Y[:,gyk0])))
	DGRres = DGR(exp.(hcat(Xk[:,.!gxk0]./gxk[:,.!gxk0],Xk[:,gxk0])),exp.(hcat(Yk[:,.!gyk0]./gyk[:,.!gyk0],Yk[:,gyk0])))

	beta = log(1./geteff(DGRres))
	return DEAResult(beta,beta.*gxk,beta.*gyk,getpeers(DGRres))
end

function (D::DDF{S,T})() where {S<:AbstractDataEnvelopment,T<:RS}
  	Data = getdata(D)
  	res = Array{DEAResult}(getnrdmu(Data))

    #@sync @parallel
  	for k in eachindex(Data)
  		res[k] = D(Data[k]...)
  	end

  	return res
end

max_stars_repo_name: kepiej/NEPA


text:
#!/usr/bin/julia

using Pkg
pkg"activate ." # Activate enviroment

using DataFrames, CSV
using Statistics
using Makie, StatsMakie

# Directory of data and files
ECODIR="assets/ECOBICI";
FILE="2019-09.csv";

# Load data from csv to a DataFrame
dt = DataFrame(CSV.File("$ECODIR/$FILE"));
show(dt)

sleep(5)

# Create scene of Makie
scene = Scene()

display(scene)

println("")
println("Scene created!")

# Plot hitograms of "Estacion_Retiro" and "Estacion_Arribo"
plot!(scene, histogram(nbins=200 ), filter(x -> x < 100, sort(dt[!, 4] .+ 0.5 ) ), color = :gray )
plot!(scene, histogram(nbins=200 ), filter(x -> x < 100, sort(dt[!, 7] ) ), color = :black )

display(scene)

# ProgressMeter allows to see a nice graph in terminal
using ProgressMeter

@showprogress for i in 1:20
    sleep(1)
end


text:
module ThreePlayer

import Cribbage: CribbageGame, GameState, Player, nextplayer, UnexpectedPlayerException
import Cribbage.TwoPlayer: incrementscore, numplayers, cardsdealt, cardstocrib, players, reportscore, cardstocrib, isgameover, rotatedealer

export incrementscore, isgameover, players, reportscore, rotatedealer, ThreePlayerGame, ThreePlayerGameState

struct ThreePlayerGame <: CribbageGame
    p₁ :: Player
    p₂ :: Player
    p₃ :: Player
    winningscore :: Int8
    function ThreePlayerGame(
        p₁::Player, p₂::Player, p₃::Player, winningscore::Integer=121)
        @assert(winningscore == 121 || winningscore == 61, "winningscore is either 61 or 121")
        @assert(p₁ ≢ p₂)
        @assert(p₁ ≢ p₃)
        @assert(p₂ ≢ p₃)
        new(p₁, p₂, p₃, winningscore)
    end
end

struct ThreePlayerGameState <: GameState
    game :: ThreePlayerGame
    dealer :: Player
    s₁ :: Int8
    s₂ :: Int8
    s₃ :: Int8
    function ThreePlayerGameState(
            game::ThreePlayerGame,
            dealer::Player,
            s₁::Integer=0,
            s₂::Integer=0,
            s₃::Integer=0
        )
        @assert(dealer ∈ players(game), "dealer is not one of the players of the game")
        @assert(0 ≤ s₁ ≤ game.winningscore, "s₁ is an invalid score")
        @assert(0 ≤ s₂ ≤ game.winningscore, "s₂ is an invalid score")
        @assert(0 ≤ s₃ ≤ game.winningscore, "s₃ is an invalid score")
        new(game, dealer, s₁, s₂, s₃)
    end
end

function numplayers(game::ThreePlayerGame)::Int8
    return 3
end

function cardsdealt(game::ThreePlayerGame)::Int8
    return 5
end

function cardstocrib(game::ThreePlayerGame)::Int8
    return 1
end

function players(game::ThreePlayerGame)::Array{Player,1}
    return [ game.p₁, game.p₂, game.p₃ ]
end

function incrementscore(gamestate::ThreePlayerGameState, player::Player, amount::Integer)
    @assert(amount ≥ 0)
    if player === gamestate.game.p₁
        return ThreePlayerGameState(
            gamestate.game,
            gamestate.dealer,
            min(gamestate.s₁ + amount, gamestate.game.winningscore),
            gamestate.s₂,
            gamestate.s₃
        )
    elseif player === gamestate.game.p₂
        return ThreePlayerGameState(
            gamestate.game,
            gamestate.dealer,
            gamestate.s₁,
            min(gamestate.s₂ + amount, gamestate.game.winningscore),
            gamestate.s₃
        )
    elseif player === gamestate.game.p₃
        return ThreePlayerGameState(
            gamestate.game,
            gamestate.dealer,
            gamestate.s₁,
            gamestate.s₂,
            min(gamestate.s₃ + amount, gamestate.game.winningscore)
        )
    else
        throw(UnexpectedPlayerException("$(player.name) is not one of the players in this game"))
    end
end

function isgameover(gamestate::ThreePlayerGameState)::Bool
    return (
        gamestate.s₁ == gamestate.game.winningscore ||
        gamestate.s₂ == gamestate.game.winningscore ||
        gamestate.s₃ == gamestate.game.winningscore
    )
end

function reportscore(gamestate::ThreePlayerGameState)::String
    return "$(gamestate.s₁) - $(gamestate.s₂) - $(gamestate.s₃)"
end

function rotatedealer(gamestate::ThreePlayerGameState)::ThreePlayerGameState
    return ThreePlayerGameState(
        gamestate.game,
        nextplayer(players(gamestate.game), gamestate.dealer),
        gamestate.s₁,
        gamestate.s₂,
        gamestate.s₃
    )
end

end  # module ThreePlayer


text:
# Day 13: Transparent Origami
function load(::Day{13}, path)
    positions = Set{CartesianIndex}()
    folds = Vector{Tuple{Char,Int}}()
    regex = r"fold along ([xy])=(\d+)"
    for line in readlines(path)
        isempty(line) && continue
        instruction = match(regex, line)
        if isnothing(instruction)
            x, y = parse.(Int, split(line, ","))
            push!(positions, CartesianIndex(x, y))
        else
            push!(folds, (first(instruction[1]), parse(Int, instruction[2])))
        end
    end
    return (positions, folds)
end

# Part 1 – single fold
@inline function foldx(I::CartesianIndex, pos)
    x, y = Tuple(I)
    return CartesianIndex(x - 2 * abs(pos - x), y)
end
@inline function foldy(I::CartesianIndex, pos)
    x, y = Tuple(I)
    return CartesianIndex(x, y - 2 * abs(pos - y))
end

function fold!(positions, fold)
    pfold = last(fold)
    if first(fold) == 'x'
        folded = [
            foldx(pop!(positions, p), pfold) for p in filter(p -> p[1] > pfold, positions)
        ]
    else
        folded = [
            foldy(pop!(positions, p), pfold) for p in filter(p -> p[2] > pfold, positions)
        ]
    end
    union!(positions, folded)
    return nothing
end

function solve1(::Day{13}, data)
    positions, folds = data
    fold!(positions, first(folds))
    return length(positions)
end

# Part 2 – show the final code
function show_code(positions::Set{CartesianIndex})
    pos = Tuple.(positions)
    h = maximum(first(t) + 1 for t in pos)
    w = maximum(last(t) + 1 for t in pos)
    img = zeros(Bool, h, w)
    for p in positions
        img[p + CartesianIndex(1, 1)] = 1
    end
    return transpose(img)
end

function solve2(::Day{13}, data)
    positions, folds = data
    for f in folds
        fold!(positions, f)
    end
    # @show show_code(positions) # worth plugging through ImageCore's `Gray`
    return 0
end

testresult1(::Day{13}) = 17

# No test for part 2, just pass:
testresult2(::Day{13}) = 0
test2(::Day{13}) = 0


text:
module LocalCongruence
	using DataStructures
	using LinearAlgebraicRepresentation
	using NearestNeighbors
	using SparseArrays
	using GraphBLASInterface
	using SuiteSparseGraphBLAS
	using SparseMM
	Lar = LinearAlgebraicRepresentation

	include("./verticesCongruence.jl")
	include("./cea-AA.jl")
	include("./cea-SM.jl")
	include("./cea-GB.jl")

	export chainCongruenceSM, chainCongruenceAA, chainCongruenceGB
end


text:
## initialize module
using Pkg;Pkg.activate(".")
include("../src/RecursiveClusterAlgorithm.jl")
using .RecursiveClusterAlgorithm

using VortexDistributions
using Test, Plots

## set up mock data
N = 10
xi = rand(N)
yi = rand(N)
tree = spanning_tree(xi,yi)

## make cluster using vortices and tree
vort = PointVortex.(xi,yi,one.(xi))
c1 = Cluster(vort,tree)
@test c1.vortices |> length == N

## plot cluster
p1 = plot()
plot_cluster!(p1,c1) # check that the plot method still works

## do some random tests of robustness
N = 30
xi,yi = randn(N),randn(N)
tree = spanning_tree(xi,yi)
c1 = Cluster(PointVortex.(xi,yi,one.(xi)),tree)
p1 = plot(grid=false,axis=false)
plot_cluster!(p1,c1)

xi,yi = randn(N),randn(N)
xi .+= 5
tree = spanning_tree(xi,yi)
c1 = Cluster(PointVortex.(xi,yi,-one.(xi)),tree)
plot_cluster!(p1,c1)


text:
module TestGB

using FactCheck
importall GradientBoost.GB
importall GradientBoost.LossFunctions

type DummyGradientBoost <: GBAlgorithm; end

type StubGradientBoost <: GBAlgorithm
  loss_function::LossFunction
  sampling_rate::FloatingPoint
  learning_rate::FloatingPoint
  num_iterations::Int
end

function build_base_func(
  gb::StubGradientBoost,
  instances,
  labels,
  prev_func_pred,
  psuedo)

  function pred(instances)
    num_instances = size(instances, 1)
    predictions = Array(Float64, num_instances)
    for i = 1:num_instances
      predictions[i] = sum(instances[i,:])
    end
    predictions
  end

  model_const = 0.5
  return (instances) -> model_const .* pred(instances)
end

sgb_instances = [
  2 2;
  2 4
]
sgb_labels = [
  1.0;
  3.0
]

facts("Gradient Boost") do
  context("not implemented functions throw an error") do
    dgb = DummyGradientBoost()
    emp_mat = Array(Any, 1, 1)
    emp_vec = Array[]
    @fact_throws build_base_func(
      dgb,
      emp_mat,
      emp_vec,
      emp_vec,
      emp_vec
    )
  end

  context("stochastic_gradient_boost works") do
    # Sanity check
    sgb = StubGradientBoost(LeastSquares(), 1.0, 0.5, 1)
    model = stochastic_gradient_boost(sgb, sgb_instances, sgb_labels)
    @fact 1 => 1
  end

  context("fit returns model") do
    sgb = StubGradientBoost(LeastSquares(), 1.0, 0.5, 1)
    model = stochastic_gradient_boost(sgb, sgb_instances, sgb_labels)
    @fact typeof(model) <: GBModel => true
  end

  context("predict works") do
    expected = [
      3.0
      3.5
    ]
    sgb = StubGradientBoost(LeastSquares(), 1.0, 0.5, 1)
    model = stochastic_gradient_boost(sgb, sgb_instances, sgb_labels)
    predictions = predict(model, sgb_instances)
    @fact predictions => expected
  end

  context("create_sample_indices works") do
    instances = [1:5 6:10]
    labels = [1:5]

    sgb = StubGradientBoost(LeastSquares(), 1, 1, 1)
    indices = create_sample_indices(sgb, instances, labels)
    @fact length(indices) => 5
    @fact length(unique(indices)) => 5
    @fact minimum(indices) >= 1 => true
    @fact maximum(indices) <= 5 => true

    sgb = StubGradientBoost(LeastSquares(), 0.5, 1, 1)
    indices = create_sample_indices(sgb, instances, labels)
    @fact length(indices) => 3
    @fact length(unique(indices)) => 3
    @fact minimum(indices) >= 1 => true
    @fact maximum(indices) <= 5 => true
  end
end

end # module


text:
using Dates, Learn, BusinessDays

struct HolidaysParsed
    weekday::String
    date::Date
    name::String
    localname::String
    holidaytype::String
    comments::String
end
# HolidaysParsed(weekday, date, name, localname, holidaytype, comments) = begin
#     date = Date(date.children[1].attributes["datetime"])
#     HolidaysParsed(nodeText(weekday), date, nodeText(name), nodeText(localname), nodeText(holidaytype), nodeText(comments))
# end

# How I got the testcases for sweden

# using HTTP, Gumbo, Cascadia
# HolidaysParsed(weekday, date, name, localname, holidaytype, comments) = begin
#     date = Date(date.children[1].attributes["datetime"])
#     HolidaysParsed(nodeText(weekday), date, nodeText(name), nodeText(localname), nodeText(holidaytype), nodeText(comments))
# end

# function fetch_holidays(yr::Year)
#     url = "https://www.officeholidays.com/countries/sweden/$(yr.value)"
#     r = HTTP.get(url; require_ssl_verification=false)
#     page = parsehtml(String(r.body))
#     rows = eachmatch(sel".country-table tr", page.root)[2:end]

#     holidays = [HolidaysParsed(eachmatch(sel"td", row)...) for row in rows]
#     [h for h in holidays if h.holidaytype == "National Holiday"]
# end

# hs = [((fetch_holidays(Year(yr)) for yr in 2015:2022)...)...]

# open(joinpath(@__DIR__,"./testHolidaysSweden.txt"), "w") do io
#     for h in hs
#         println(io, "$(h.weekday),$(h.date),$(h.name),$(h.localname),$(h.holidaytype),$(h.comments)")
#     end
# end

HolidaysParsed(weekday, date, name, localname, holidaytype, comments) = begin
    date = Date(date)
    HolidaysParsed(String(weekday), date, String(name), String(localname), String(holidaytype), String(comments))
end

hs = open(joinpath(@__DIR__,"./testHolidaysSweden.txt")) do io
    readlines(io) .|> x -> split(x,",") |> x->HolidaysParsed(x...)
end

holidays = Dict(h.date=>h for h in hs)

@testset "SwedishHolidays" begin
    for d ∈ Date(2015,1,1):Day(1):Date(2022,12,31)
        @test BusinessDays.isweekend(d) || ((d ∉ holidays.keys) == isbday(SwedenCalendar(),d))
    end
end

max_stars_repo_name: binnisb/JuliaLearn

max_issues_repo_name: binnisb/JuliaLearn

max_forks_repo_name: binnisb/JuliaLearn


text:
rrule(::Type{<:Symmetric}, A::AbstractMatrix) = Symmetric(A), Rule(_symmetric_back)

_symmetric_back(ΔΩ) = UpperTriangular(ΔΩ) + LowerTriangular(ΔΩ)' - Diagonal(ΔΩ)
_symmetric_back(ΔΩ::Union{Diagonal, UpperTriangular}) = ΔΩ


text:
import Polynomials
using Polynomials

#
# given the coefficients of a polynomial p(x) = p(n)*x^n + + p(1)*x + p(0)
#
# determne a scaling factor for x and a linear factor t that the transformed
# polynomial q(z) = t * p(s*z) has coefficient p(n) = 1 and all |p(j)| <= 1
#
# Invariant:
# if q, shift, smax, u = polyscale(p) then q = polytransform(p, shift, smax) * u
# up to numerical precision

FloatOrComplex = Union{AbstractFloat, Complex, Complex{Float32}, Complex{BigFloat}}

function polyscale(P::Poly{T}) where {T<:FloatOrComplex}
	  n = length(P.a) -1
    shift = - P.a[n] / n
    P = polytransform(P, shift)
	  p = P.a
    p[n] = 0
    pni = 1 / p[n+1]
    smaxlog = maximum(log(abs(p[1:n] * pni)) ./ (n+1 - (1:n)))
    smax = exp(smaxlog)
    u = pni / exp(smaxlog * n)
    P = polytransform(Poly(p), 0.0, smax) * u
    P.a[n+1] = T(1)
    P, shift, smax, u
end

"""
given a monic polynomial (x^n + ...) determine a scaling factor s that the transformed
polynomial q(z) = p(s*z) is half of the time < 1, rest >= 1. Ignore zero coefficients.
"""
function polyscale2(P::Poly{T}) where T<:FloatOrComplex

  p = abs.(P.a)
  n = length(p) - 1
  pni = T(1) / p[n+1]
  scale!(p, pni)
  p = view(p,1:n)
  bigg = log(realmax(T))
  tiny = log(realmin(T))
  clip(x::T) = x == -Inf ? -x : x

  llim = exp(maximum(((log.(p) - bigg) ./ (n:-1:1))))
  ulim = exp(minimum(clip, ((log.(p) - tiny) ./ (n:-1:1))))
  med = exp(sort(log.(p) ./ (n:-1:1))[n÷2+1])

  med, llim, ulim

end

"""

Find scaling factor for a polynomial p_1 + p_2*x^1 + ... + p_(n+1)*x^n,
which has the minimal possible quotient of coefficients.
Coefficients which are zero ignored.
"""
function polyscale_mini_quotient(P::Poly{T}) where T<:FloatOrComplex
  p = abs.(P.a)
  n = length(p) - 1
  pni = log(abs(p[n+1]))

  plog = log.(p[1:n]) - pni

  f = Array{real(T)}(2n)
  copy!(f, plog)
  copy!(f, n+1, -f, 1, n)
  g = Array{real(T)}(2n)
  g =
  copy!(g, -(n:-1:1))
  copy!(g, n+1, -g, 1, n)

  xa, xb, val = minimize(f, g)
  exp(xa), exp(xb), exp(val)
end

"""

  `minimize(f::Array, g::Array of same size) -> xa, xb, vmin`

Solve the following special linear optimization problem:
Define `f(x) = max{ f1[k] + g1[k] * x / k = 1..n1}`.
Minimize `f(x)` for real x.
The function f is convex as maximum of linear (convex) functions is convex.
If f is bounded below, it has minmum value `vmin` and a solution interval `[xa, xb]`
with `f(x) = vmin for all x in [xa, xb]`.
If f is not bounded below, return `vmin = -Inf` and `xa = xb = ±Inf`.
"""
function minimize(f::AbstractVector{T}, g::AbstractVector{T}) where T<:Real

  n::Int = length(f)
  n == length(g) || error(ArgumentError("lengths of f1 and g1 are different"))

  Z = zero(T)
  inf = typemax(T)

  # remove unwanted infinities
  f, g = unzip(Iterators.filter(x -> x[1] !=  inf, zip(f, g)))::NTuple{2,AbstractArray{T}}
  n = length(f)

  if n == 0
    return -inf, inf, -inf
  elseif n == 1
    if g[1] == Z
      return -inf, inf, f[1]
    elseif g[1] < Z
      return inf, inf, -inf
    else
      return -inf, -inf, -inf
    end
  end

  eva(x::T)::Tuple{T,T,Int} = maximum(1:n) do k; ((f[k] + g[k] * x), g[k], k) end

  function nextto(x::T)::Tuple{T, T, T}
    v, gk, k = eva(x)
    cross(j::Int) = -(f[j] - f[k]) / (g[j] - gk)

    if isnan(v)
      xa = xb = x
    elseif gk > Z
      xa = xb = maximum(filter(y->y<x, cross.(filter(j->j!=k, 1:n))))
    elseif gk < Z
      xa = xb = minimum(filter(y->y>x, cross.(filter(j->j!=k, 1:n))))
    elseif gk == Z
      xa = maximum(filter(y->y<x, cross.(filter(j->j!=k, 1:n))))
      xb = minimum(filter(y->y>x, cross.(filter(j->j!=k, 1:n))))
    else
      xa, xb = inf, -inf
    end
    xa, xb, gk
  end

  # find extreme slopes
  gmax::T, fmax::T, kmax::Int = maximum(1:n) do k; (g[k], f[k], k) end
  gmin::T, fmin::T, kmin::Int = minimum(1:n) do k; (g[k], f[k], k) end

  xmax::T = maximum(filter(x->!isnan(x), -(f-f[kmax]) ./ (g - gmax)))
  xmin::T = minimum(filter(x->!isnan(x), -(f-f[kmin]) ./ (g - gmin)))
  xa, xb = xmax, xmin

  while xmin < xa || xb < xmax
    xa, xb, gab = nextto((xmin + xmax) / 2)
    gs::Int = cmp(gab, Z)
    if gs >= 0
      xmax = xb
    end
    if gs <= 0
      xmin = xa
    end
  end
  xa, xb, eva((xa + xb)/ 2)...
end

maximum1(f::Function, v0, itr) = mapreduce(f, Base.scalarmax, v0, itr)
maximum1(v0, itr) = mapreduce(identity, Base.scalarmax, v0, itr)
maximum1(f::Function, itr) = mapreduce(f, Base.scalarmax, lowerbound(f, itr), itr)
maximum1(itr) = mapreduce(identity, Base.scalarmax, lowerbound(itr), itr)
lowerbound(itr) = lowerbound(eltype(itr))
lowerbound(f::Function, itr) = lowerbound(promote_type(Base.return_types(f, (eltype(itr),))))
lowerbound(T::Type) = error("no lower bound for type $T")
lowerbound(::Type{T}) where T<:Union{AbstractFloat,Unsigned} = typemin(T)
lowerbound(S::Type{Rational{T}}) where T<:Integer = typemin(S)
lowerbound(::Type{String}) = ""

#Base.mr_empty(::typeof(identity),::typeof(Base.scalarmax),::Type{T}) where T<:Real = typemin(T)
#Base.mr_empty(::typeof(identity),::typeof(Base.scalarmin),::Type{T}) where T<:Real = typemax(T)

"""
  `unzip(itr)`

The iterable object itr must deliver N-tuples.
Convert an array of N-tuples into N-tuple of arrays.
"""
function unzip(aot)
  tuple(map(x->Array(collect(x)), zip(aot...))...)
end


text:
#Residential Consumer Balancing Code
function annual_energy_balance(consumer::C, pvsys::PVSystem; print_output=false) where {C <: Consumer}

    results = Dict{Int64,Any}()

    d_d = Truncated(Normal(0.0, 0.5), -0.10, 0.10)
    d_s = Truncated(Normal(-0.2, 0.5), -0.10, 0.15)
    days_p_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    m_h = 0
    count = 0

    # Initial conditions
    carry_over = 0.0
    global_allowance = 0.0
    global_withdrawal = 0.0
    global_generation = 0.0

  #Loop for every Month of the year
    for (ix,m) in enumerate(days_p_month)

        #Get solar Data with houtly resolution for the month
        solar_month_base = pvsys.time_series[24*(m_h) + 1:24*(m_h+m)]
        solar_month = [(s+rand(d_s, 1)[1])*(s > 0.1) for s in solar_month_base]*pvsys.capacity
        m_h = m

        #Simulate month on an hourly basis
        #initial conditions for the month
        allowance = 0.0
        utility_supplied_energy = 0.0
        grid_energy = 0.0
        withdrawn_energy = 0.0
        consumer_energy = 0.0
        PV_energy = 0.0
        injection_grid = 0.0
        peak_power = 0.0
        peak_demand = 0.0

		peak_power_peak = 0.0
        peak_demand_peak = 0.0

		peak_power_valley = 0.0
        peak_demand_valley = 0.0

		peak_power_night = 0.0
        peak_demand_night = 0.0

        # Initial monthly conditions
        peak = [10, 11, 12, 17, 18, 19]
        valley = [6, 7, 8, 9, 13, 14, 15, 16]
        night = [20, 21, 22, 23, 24, 1, 2, 3, 4, 5]

        consumer_energy_peak = 0.0
        consumer_energy_valley = 0.0
        consumer_energy_night = 0.0

        grid_energy_peak = 0.0
        grid_energy_valley = 0.0
        grid_energy_night = 0.0

        #Loop for every day of the month
        for i in 1:m
            var = rand(d_d, 24).+1.0
            daily_p = (consumer.load_curve.*var)*consumer.peak_power[ix]
            solar_day = solar_month[(((m-1)*24)+1):(((m-1)*24)+24)]

            #Loop for every hour of the day
            for t in 1:24

               consumer_energy += daily_p[t]

               (t in peak) ?   consumer_energy_peak += daily_p[t] : true
               (t in valley) ? consumer_energy_valley += daily_p[t] : true
               (t in night) ?  consumer_energy_night += daily_p[t] : true

               PV_energy +=  solar_day[t]

               balance = daily_p[t] - solar_day[t]

               (t in peak) ? grid_energy_peak += max(0.0,balance) : true
               (t in valley) ? grid_energy_valley += max(0.0,balance) : true
               (t in night) ? grid_energy_night += max(0.0,balance) : true

               grid_energy = grid_energy_peak + grid_energy_valley + grid_energy_night

               injection_grid += max(0.0, -1*balance)

               peak_power = max(peak_power, daily_p[t])

               peak_demand = max(peak_demand, balance)

			   #Peak demand corresponds to the peak as measured by the distribution company

               #(t in peak) ? println("peak d"," ",  peak_demand_peak, " ", balance, " ", max(0.0, peak_demand_peak, max(0.0,balance)) ) : true
               (t in peak) ? peak_demand_peak = max(0.0, peak_demand_peak, max(0.0,balance)) : true

               #(t in valley) ? println("valley d", " ", peak_demand_valley," ",  balance, " ", max(0.0, peak_demand_valley, max(0.0,balance)) ) : true
               (t in valley) ? peak_demand_valley = max(0.0, peak_demand_valley, max(0.0,balance)) : true

               #(t in night) ? println("night d", " ", peak_demand_night, " ", balance, " ", max(0.0, peak_demand_night, max(0.0,balance)) ) : true
               (t in night) ? peak_demand_night = max(0.0, peak_demand_night, max(0.0,balance)) : true

               #Peak power corresponds to the peak of the consumer load

               #(t in peak) ? println("peak p", " ", peak_power_peak, " ", daily_p[t], " ", max(0.0, peak_power_peak, daily_p[t]) ) : true
               (t in peak) ? peak_power_peak = max(0.0, peak_power_peak, daily_p[t]) : true

               #(t in valley) ? println("valley p", " ", peak_power_valley, " ", daily_p[t], " ", max(0.0, peak_power_valley, daily_p[t])) : true
               (t in valley) ? peak_power_valley = max(0.0, peak_power_valley, daily_p[t]) : true

               #(t in night) ? println("night p", " ", peak_power_night, " ", daily_p[t], " ", max(0.0, peak_power_night, daily_p[t]) ) : true
               (t in night) ? peak_power_night = max(0.0, peak_power_night, daily_p[t]) : true

            end #end of hourly loop

        end #end of daily loop

            global_generation += PV_energy #Total Energy from Rooftop PV
            global_allowance = 0.49*global_generation - global_withdrawal # Uses the acumulated withdrawls up to month ix-1
            allowance = min(injection_grid+carry_over, global_allowance)
            withdrawn_energy = min(min(grid_energy,injection_grid+carry_over),allowance) #Withdrawn Energy From the Grid
            utility_supplied_energy = max(grid_energy - injection_grid - carry_over,0.0)

            #Update quantities for next month
            carry_over = max(allowance - withdrawn_energy, 0)
            global_withdrawal += withdrawn_energy #Total Energy withdrawn resulting from injections into the grid

          utility = Dict("consumer_energy_peak" => consumer_energy_peak,
                    "consumer_energy_valley" => consumer_energy_valley,
                    "consumer_energy_night" => consumer_energy_night,
                    "grid_energy_peak" => grid_energy_peak,
                    "grid_energy_valley" => grid_energy_valley,
                    "grid_energy_night" =>  grid_energy_night,
					"consumer_demand_peak" => peak_power_peak,
                    "consumer_demand_valley" => peak_power_valley,
                    "consumer_demand_night" => peak_power_night,
                    "grid_demand_peak" => peak_demand_peak,
                    "grid_demand_valley" => peak_demand_valley,
                    "grid_demand_night" =>  peak_demand_night
               )

           results[ix] = Dict("consumer_energy" => consumer_energy,
                   "PV_energy" => PV_energy,
                    "injection_grid" => injection_grid,
                    "withdrawn_energy" => withdrawn_energy,
                    "utility_supplied_energy" => utility_supplied_energy,
                "grid_energy" => grid_energy,
                "peak_demand" => peak_demand,
                "peak_power" => peak_power,
                "carry_over" => carry_over,
                "allowance" => allowance,
                "global_generation" => global_generation,
                "global_withdrawal" => global_withdrawal,
                "global_allowance" => global_allowance,
                "utility_balance" => utility)

    end #end of monthly loop

    if print_output
        print_base(results)
    end

    return results

end

#TMT Consumer Balancing Code
function annual_energy_balance(consumer::TMT, pvsys::PVSystem; print_output=false)
   results = Dict{Int64,Any}()

    d_d = Truncated(Normal(0.0, 0.5), -0.10, 0.10)
    d_s = Truncated(Normal(-0.2, 0.5), -0.10, 0.15)
    days_p_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    m_h = 0
    count = 0

    # Initial conditions
    carry_over = 0.0
    global_allowance = 0.0
    global_withdrawal = 0.0
    global_generation = 0.0

    #Loop for every Month of the year
    for (ix,m) in enumerate(days_p_month)

        #Simulate month on an hourly basis
        #initial conditions for the month
        allowance = 0.0
        utility_supplied_energy = 0.0
        grid_energy = 0.0
        withdrawn_energy = 0.0
        consumer_energy = 0.0
        PV_energy = 0.0

		peak_power_peak = 0.0
        peak_demand_peak = 0.0

		peak_power_valley = 0.0
        peak_demand_valley = 0.0

		peak_power_night = 0.0
        peak_demand_night = 0.0

		consumer_energy_peak = 0.0
        consumer_energy_valley = 0.0
        consumer_energy_night = 0.0

		grid_energy_peak = 0.0
        grid_energy_valley = 0.0
        grid_energy_night = 0.0

        total_injection = 0.0
		injection_grid_peak = 0.0
		injection_grid_valley = 0.0
		injection_grid_night = 0.0

        # Initial monthly conditions
        peak = [10, 11, 12, 17, 18, 19]
        valley = [6, 7, 8, 9, 13, 14, 15, 16]
        night = [20, 21, 22, 23, 24, 1, 2, 3, 4, 5]

        #Get solar Data with houtly resolution for the month
        solar_month_base = pvsys.time_series[24*(m_h) + 1:24*(m_h+m)]
        solar_month = [(s+rand(d_s, 1)[1])*(s > 0.1) for s in solar_month_base]*pvsys.capacity
        m_h = m

        #Loop for every day of the month
        for i in 1:m
            var = rand(d_d, 24).+1.0
            daily_p = (consumer.load_curve.*var)*consumer.peak_power[ix]
            solar_day = solar_month[(((m-1)*24)+1):(((m-1)*24)+24)]

            #Loop for every hour of the day
            for t in 1:24

			   consumer_energy += daily_p[t]

               (t in peak) ?   consumer_energy_peak += daily_p[t] : true
               (t in valley) ? consumer_energy_valley += daily_p[t] : true
               (t in night) ?  consumer_energy_night += daily_p[t] : true

               PV_energy +=  solar_day[t]

               balance = daily_p[t] - solar_day[t]

               (t in peak) ? grid_energy_peak += max(0.0,balance) : true
               (t in valley) ? grid_energy_valley += max(0.0,balance) : true
               (t in night) ? grid_energy_night += max(0.0,balance) : true

               grid_energy = grid_energy_peak + grid_energy_valley + grid_energy_night

               (t in peak) ?  injection_grid_peak += max(0.0, -1*balance) : true
			   (t in valley) ? injection_grid_valley += max(0.0, -1*balance) : true
			   (t in night) ? injection_grid_night += max(0.0, -1*balance) : true

			   total_injection = injection_grid_peak + injection_grid_valley + injection_grid_night

               #Peak demand corresponds to the peak as measured by the distribution company

               #(t in peak) ? println("peak d"," ",  peak_demand_peak, " ", balance, " ", max(0.0, peak_demand_peak, max(0.0,balance)) ) : true
               (t in peak) ? peak_demand_peak = max(0.0, peak_demand_peak, max(0.0,balance)) : true

               #(t in valley) ? println("valley d", " ", peak_demand_valley," ",  balance, " ", max(0.0, peak_demand_valley, max(0.0,balance)) ) : true
               (t in valley) ? peak_demand_valley = max(0.0, peak_demand_valley, max(0.0,balance)) : true

               #(t in night) ? println("night d", " ", peak_demand_night, " ", balance, " ", max(0.0, peak_demand_night, max(0.0,balance)) ) : true
               (t in night) ? peak_demand_night = max(0.0, peak_demand_night, max(0.0,balance)) : true

               #Peak power corresponds to the peak of the consumer load

               #(t in peak) ? println("peak p", " ", peak_power_peak, " ", daily_p[t], " ", max(0.0, peak_power_peak, daily_p[t]) ) : true
               (t in peak) ? peak_power_peak = max(0.0, peak_power_peak, daily_p[t]) : true

               #(t in valley) ? println("valley p", " ", peak_power_valley, " ", daily_p[t], " ", max(0.0, peak_power_valley, daily_p[t])) : true
               (t in valley) ? peak_power_valley = max(0.0, peak_power_valley, daily_p[t]) : true

               #(t in night) ? println("night p", " ", peak_power_night, " ", daily_p[t], " ", max(0.0, peak_power_night, daily_p[t]) ) : true
               (t in night) ? peak_power_night = max(0.0, peak_power_night, daily_p[t]) : true

            end #end of hourly loop

		end #end of daily loop

            global_generation += PV_energy #Total Energy from Rooftop PV
            global_allowance = 0.49*global_generation - global_withdrawal # Uses the acumulated withdrawls up to month ix-1
            allowance = min(total_injection+carry_over, global_allowance)

			if !(carry_over > 0.0)

				withdrawn_energy = min(min(grid_energy,total_injection),allowance) #Withdrawn Energy From the Grid
				utility_supplied_energy_peak = max(grid_energy_peak - injection_grid_peak,0.0)
				utility_supplied_energy_valley = max(grid_energy_valley- injection_grid_valley,0.0)
				utility_supplied_energy_night = max(grid_energy_night,0.0)
                #Update quantities for next month
				carry_over = max(allowance - withdrawn_energy, 0)
			elseif carry_over > 0.0
                temp_carry_over = 0.0

				withdrawn_energy = min(min(grid_energy,total_injection +carry_over),allowance) #Withdrawn Energy From the Grid

				grid_energy_valley - injection_grid_valley > 0 ? temp_carry_over = min(carry_over, grid_energy_valley - injection_grid_valley) : temp_carry_over = 0.0
				utility_supplied_energy_valley = max((grid_energy_valley - injection_grid_valley - temp_carry_over), 0.0)

				grid_energy_peak - injection_grid_peak > 0 ? temp_carry_over = min(carry_over - temp_carry_over, grid_energy_valley - injection_grid_valley) : temp_carry_over = 0.0
                utility_supplied_energy_peak = max((grid_energy_peak - injection_grid_peak - temp_carry_over), 0.0)

				temp_carry_over = min(carry_over - temp_carry_over, grid_energy_night)
				utility_supplied_energy_night = max(grid_energy_night - temp_carry_over,0.0)

				carry_over = max(allowance - withdrawn_energy, 0)
		   end

            global_withdrawal += withdrawn_energy #Total Energy withdrawn resulting from injections into the grid

          utility = Dict("consumer_energy_peak" => consumer_energy_peak,
                    "consumer_energy_valley" => consumer_energy_valley,
                    "consumer_energy_night" => consumer_energy_night,
                    "grid_energy_peak" => grid_energy_peak,
                    "grid_energy_valley" => grid_energy_valley,
                    "grid_energy_night" =>  grid_energy_night,
					"consumer_demand_peak" => peak_power_peak,
                    "consumer_demand_valley" => peak_power_valley,
                    "consumer_demand_night" => peak_power_night,
                    "grid_demand_peak" => peak_demand_peak,
                    "grid_demand_valley" => peak_demand_valley,
                    "grid_demand_night" =>  peak_demand_night
               )

           results[ix] = Dict("consumer_energy_peak" => consumer_energy_peak,
                              "consumer_energy_valley" => consumer_energy_valley,
                              "consumer_energy_night" => consumer_energy_night,
                              "PV_energy" => PV_energy,
							  "grid_energy" => grid_energy,
                              "total_injection" => total_injection,
                              "withdrawn_energy" => withdrawn_energy,
							  "utility_supplied_energy_peak" => utility_supplied_energy_peak,
                              "utility_supplied_energy_valley" => utility_supplied_energy_valley,
                              "utility_supplied_energy_night" => utility_supplied_energy_night,
                              "grid_energy_peak" => grid_energy_peak,
                              "grid_energy_valley" => grid_energy_valley,
                              "grid_energy_night" => grid_energy_night,
                              "peak_demand_peak" => peak_demand_peak,
                              "peak_demand_valley" => peak_demand_valley,
                              "peak_demand_night" => peak_demand_night,
                              "peak_power_peak" => peak_power_peak,
                              "peak_power_valley" => peak_power_valley,
                              "peak_power_night" => peak_power_night,
                              "carry_over" => carry_over,
                              "global_generation" => global_generation,
                              "global_withdrawal" => global_withdrawal,
                              "global_allowance" => global_allowance,
                              "utility_balance" => utility,
                        )

    end

    if print_output
        print_TMT(results)
    end

    return results

end

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
using PackageCompiler
create_sysimage([:Pluto, :PlutoUI, :DifferentialEquations, :Plots];
                precompile_execution_file = "warmup.jl",
                replace_default = true,
                cpu_target = PackageCompiler.default_app_cpu_target())

max_stars_repo_name: matutenun/projectbinder

max_issues_repo_name: matutenun/projectbinder

max_forks_repo_name: matutenun/projectbinder


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ 538f05be-5053-11eb-19a0-959e34e1c2a1
using Pkg, DrWatson

# ╔═╡ 5d84f90c-5053-11eb-076b-5f30fc9685e3
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample, StanOptimize
	using StatisticalRethinking
end

# ╔═╡ 235c5298-5053-11eb-0608-435d1aa4716c
md" ## Clip-07-15.bs.jl"

# ╔═╡ 861accd4-5053-11eb-0432-81db212f4f38
begin
	sppnames = [:afarensis, :africanus, :hapilis, :boisei, :rudolfensis, :ergaster, :sapiens]
	brainvol = [438, 452, 612, 521, 752, 871, 1350]
	masskg = [37, 35.5, 34.5, 41.5, 55.5, 61, 53.5]
	df = DataFrame(species = sppnames, brain = brainvol, mass = masskg)
	scale!(df, [:mass, :brain])
	df.brain_std = df.brain/maximum(df.brain)
end;

# ╔═╡ 290d19ae-59cf-11eb-1636-772efccc7cb9
df

# ╔═╡ 24328f84-541d-11eb-0ec3-4df5b8a8cb19
begin
	scatter(df.mass, df.brain, xlab="body mass [kg]", ylab="brain vol [cc]",
		lab="Observations")
	for (ind, species) in pairs(df.species)
		annotate!([(df[ind, :mass] + 1, df[ind, :brain] + 30,
			Plots.text(df[ind, :species], 6, :red, :right))])
	end
	plot!()
end

# ╔═╡ ee275e7a-5067-11eb-325b-7760b758e85e
stan7_6 = "
data{
	int N;
    vector[N] brain_std;
    vector[N] mass_std;
}
parameters{
    real a;
    real b1;
    real b2;
    real b3;
    real b4;
    real b5;
    real b6;
}
model{
    vector[N] mu;
    b6 ~ normal( 0 , 10 );
    b5 ~ normal( 0 , 10 );
    b4 ~ normal( 0 , 10 );
    b3 ~ normal( 0 , 10 );
    b2 ~ normal( 0 , 10 );
    b1 ~ normal( 0 , 10 );
    a ~ normal( 0.5 , 1 );
    for ( i in 1:N ) {
        mu[i] = a + b1 * mass_std[i] + b2 * mass_std[i]^2 + b3 * mass_std[i]^3 + b4 * mass_std[i]^4 + b5 * mass_std[i]^5 + b6 * mass_std[i]^6;
    }
    brain_std ~ normal( mu , 0.001 );
}
generated quantities{
    vector[N] log_lik;
    vector[N] mu;
    for ( i in 1:N ) {
        mu[i] = a + b1 * mass_std[i] + b2 * mass_std[i]^2 + b3 * mass_std[i]^3 + b4 * mass_std[i]^4 + b5 * mass_std[i]^5 + b6 * mass_std[i]^6;
    }
    for ( i in 1:N ) log_lik[i] = normal_lpdf( brain_std[i] | mu[i] , 0.001 );
}
";

# ╔═╡ 95cbb456-520e-11eb-2e0a-1d173ddc9dec
begin
	data = (N = 7, brain_std = df.brain_std, mass_std = df.mass_s)
	m7_6s = SampleModel("m7.6s", stan7_6)
	rc7_6s = stan_sample(m7_6s; data=data)

	if success(rc7_6s)
		nt7_6s = read_samples(m7_6s)
	end
end;

# ╔═╡ 19f6114a-6349-11eb-26fa-3dc20e36601c
begin
	log_lik = nt7_6s.log_lik'
	n_sam, n_obs = size(log_lik)
	lppds = reshape(logsumexp(log_lik .- log(n_sam); dims=1), n_obs)
end

# ╔═╡ 244f83fe-66fc-11eb-0cc6-ab326a2f2df7
reshape(lppd(log_lik), n_obs)

# ╔═╡ b5f36364-6348-11eb-288d-a15bb85e31a4
size(log_lik)

# ╔═╡ 675109fc-59e6-11eb-247d-157746d5c626
sum(lppds)

# ╔═╡ 8b9f3f2a-6703-11eb-1bb0-69e52401aedd
waic(log_lik)

# ╔═╡ e8343c86-5a79-11eb-1cd4-7b05ac85f7dd
md"
```
    mean   sd  5.5% 94.5% n_eff Rhat4
a   0.51 0.01  0.50  0.52    47  1.07
b1  0.88 0.01  0.86  0.90    50  1.07
b2  1.71 0.03  1.65  1.76    46  1.07
b3 -0.61 0.03 -0.65 -0.56    49  1.07
b4 -3.48 0.05 -3.56 -3.40    48  1.06
b5 -0.35 0.02 -0.38 -0.32    50  1.06
b6  1.63 0.02  1.59  1.66    49  1.06
```
"

# ╔═╡ bc6292b4-5a75-11eb-2a11-a31ae9ae5c4f
PRECIS(read_samples(m7_6s, :dataframe))

# ╔═╡ 91e9af0a-5065-11eb-212c-f751fd114263
md" ## End of clip-07-15.bs.jl"

# ╔═╡ Cell order:
# ╟─235c5298-5053-11eb-0608-435d1aa4716c
# ╠═538f05be-5053-11eb-19a0-959e34e1c2a1
# ╠═5d84f90c-5053-11eb-076b-5f30fc9685e3
# ╠═861accd4-5053-11eb-0432-81db212f4f38
# ╠═290d19ae-59cf-11eb-1636-772efccc7cb9
# ╠═24328f84-541d-11eb-0ec3-4df5b8a8cb19
# ╠═ee275e7a-5067-11eb-325b-7760b758e85e
# ╠═95cbb456-520e-11eb-2e0a-1d173ddc9dec
# ╠═19f6114a-6349-11eb-26fa-3dc20e36601c
# ╠═244f83fe-66fc-11eb-0cc6-ab326a2f2df7
# ╠═b5f36364-6348-11eb-288d-a15bb85e31a4
# ╠═675109fc-59e6-11eb-247d-157746d5c626
# ╠═8b9f3f2a-6703-11eb-1bb0-69e52401aedd
# ╟─e8343c86-5a79-11eb-1cd4-7b05ac85f7dd
# ╠═bc6292b4-5a75-11eb-2a11-a31ae9ae5c4f
# ╟─91e9af0a-5065-11eb-212c-f751fd114263


text:
using Klara

println("    Testing likelihood_model generator...")

θ = BasicContUnvParameter(:θ, 1, signature=:low)
x = Data(:x, 2)
λ = Hyperparameter(:λ, 3)

likelihood_model([θ, x, λ])


text:
import Bibliography
using Test

for file in ["test.bib"] #, "xampl.bib"] #, "ignace_ref.bib"]
    test_import = Bibliography.import_bibtex("../examples/$file")
    # @info "test_import" test_import
    Bibliography.export_web(test_import)
    println("Test import\n $(Bibliography.export_bibtex("result.bib", test_import))\n")

    if file == "test.bib"
        selection = ["CitekeyArticle", "CitekeyBook"]
        test_select = Bibliography.select(test_import, selection)
        println("Test select\n $(Bibliography.export_bibtex("result.bib", test_select))\n")
    end

    rm("result.bib")
end

testdata = """@inproceedings{demo2020proceedings,
 organization  = {DemoOrg},
 pages         = {1--10},
 doi           = {10.1000/001-1-001-00001-1_001},
 author        = {Demo, D},
 note          = {cited by 0},
 year          = {2020},
 booktitle     = {Demo Booktitle},
 title         = {Demo Title}
}"""

write("demo.bib",testdata)
Bibliography.import_bibtex(testdata)
mybib = Bibliography.import_bibtex("demo.bib")
Bibliography.export_bibtex("demo_export.bib",mybib)
mybib2 = Bibliography.import_bibtex("demo_export.bib")

rm("demo.bib")
rm("demo_export.bib")

include("sort_bibliography.jl")


text:
"""
Upper triangular and lower triangular matrices are the same views of the underlying data.
Upper triangular matrices are column major, and lower triangular row major.

As an implementation detail, indexing with @inbounds into the 0-triangle (ie, lower triangle of an Upper Triangular matrix) is undefined. This is because that check creates a branch.
"""
abstract type AbstractTriangularMatrix{T, N, N2} <: AbstractMatrix{T} end

abstract type AbstractUpperTriangular{T, N, N2} <: AbstractTriangularMatrix{T, N, N2} end
abstract type AbstractLowerTriangular{T, N, N2} <: AbstractTriangularMatrix{T, N, N2} end
abstract type MutableUpperTriangular{T, N, N2} <: AbstractUpperTriangular{T, N, N2} end
abstract type MutableLowerTriangular{T, N, N2} <: AbstractLowerTriangular{T, N, N2} end

const MutableTriangularMatrix{T,N,N2} = Union{MutableUpperTriangular{T, N, N2}, MutableLowerTriangular{T, N, N2}}

struct StaticUpperTriangularMatrix{T, N, N2} <: AbstractUpperTriangular{T, N, N2}
    data::NTuple{N2,T}
end
struct StaticLowerTriangularMatrix{T, N, N2} <: AbstractLowerTriangular{T, N, N2}
    data::NTuple{N2,T}
end
mutable struct UpperTriangularMatrix{T, N, N2} <: MutableUpperTriangular{T, N, N2}
    data::NTuple{N2,T}
    function UpperTriangularMatrix{T, N, N2}(d::NTuple{N2,T}) where {T,N,N2}
        isbits(T) || error("Can only construct mutable isbits matrices.")
        new{T,N,N2}(d)
    end
end
mutable struct LowerTriangularMatrix{T, N, N2} <: MutableLowerTriangular{T, N, N2}
    data::NTuple{N2,T}
    function LowerTriangularMatrix{T, N, N2}(d::NTuple{N2,T}) where {T,N,N2}
        isbits(T) || error("Can only construct mutable isbits matrices.")
        new{T,N,N2}(d)
    end
end
struct PointerUpperTriangularMatrix{T, N, N2} <: MutableUpperTriangular{T, N, N2}
    data::Ptr{T}
end
struct PointerLowerTriangularMatrix{T, N, N2} <: MutableLowerTriangular{T, N, N2}
    data::Ptr{T}
end
Base.size(::AbstractTriangularMatrix{T,N}) where {T,N} = (N,N)

const PointerTriangularMatrix{T,N,N2} = Union{
    PointerUpperTriangularMatrix{T, N, N2},
    PointerLowerTriangularMatrix{T, N, N2}
}

point(A::MutableTriangularMatrix{T}) where T = Base.unsafe_convert(Ptr{T}, pointer_from_objref(A))
point(A::PointerTriangularMatrix) = A.data

function upper_triangle_quote(N, S = :A)
    N2 = big_triangle(N)
    q = quote
        @inbounds begin
        end
    end
    @static if VERSION > v"0.7-"
        qa = q.args[2].args[3].args
    else
        qa = q.args[2].args[2].args
    end
    ind = 0
    for i ∈ 1:N, j ∈ 1:i
        ind += 1
        A_i = Symbol(S, :_, ind)
        push!(qa, :( $A_i = $(S)[$j,$i] ) )
    end
    q, N2
end
function lower_triangle_quote(N, S = :A)
    N2 = big_triangle(N)
    q = quote
        @inbounds begin
        end
    end
    @static if VERSION > v"0.7-"
        qa = q.args[2].args[3].args
    else
        qa = q.args[2].args[2].args
    end
    for i ∈ 1:N, j ∈ i:N
        A_i = Symbol(S, :_, small_triangle(j)+i )
        push!(qa, :( $A_i = $(S)[$j,$i] ) )
    end
    push!(q.args, :($M{$T,$N,$N2}( @ntuple $N2 A )))
    q, N2
end

function recursive_lower_triangle_quote!(qa, N, ::Type{T}) where T

end

function UpperTriangularMatrix(A::AbstractMatrix)
    m, n = size(A)
    @assert m == n
    UpperTriangularMatrix(A, Val(n))
end
@generated function UpperTriangularMatrix(A::AbstractMatrix{T}, ::Val{N}) where {T,N}
    q, N2 = upper_triangle_quote(N, T)
    push!(q.args, :(UpperTriangularMatrix{$T,$N,$N2}( @ntuple $N2 A )))
    q
end
function StaticUpperTriangularMatrix(A::AbstractMatrix)
    m, n = size(A)
    @assert m == n
    StaticUpperTriangularMatrix(A, Val(n))
end
@generated function StaticUpperTriangularMatrix(A::AbstractMatrix{T}, ::Val{N}) where {T,N}
    q, N2 = upper_triangle_quote(N, T)
    push!(q.args, :(StaticUpperTriangularMatrix{$T,$N,$N2}( @ntuple $N2 A )))
    q
end
function LowerTriangularMatrix(A::AbstractMatrix)
    m, n = size(A)
    @assert m == n
    LowerTriangularMatrix(A, Val(n))
end
@generated function LowerTriangularMatrix(A::AbstractMatrix{T}, ::Val{N}) where {T,N}
    q, N2 = lower_triangle_quote(N, T)
    push!(q.args, :(LowerTriangularMatrix{$T,$N,$N2}( @ntuple $N2 A )))
    q
end
function StaticLowerTriangularMatrix(A::AbstractMatrix)
    m, n = size(A)
    @assert m == n
    StaticLowerTriangularMatrix(A, Val(n))
end
@generated function StaticLowerTriangularMatrix(A::AbstractMatrix{T}, ::Val{N}) where {T,N}
    q, N2 = lower_triangle_quote(N, T)
    push!(q.args, :(StaticLowerTriangularMatrix{$T,$N,$N2}( @ntuple $N2 A )))
    q
end
UpperTriangularMatrix(data::NTuple{N2,T}) where {T,N2} = UpperTriangularMatrix(data,ValI(Val{N2}()))
StaticUpperTriangularMatrix(data::NTuple{N2,T}) where {T,N2} = StaticUpperTriangularMatrix(data,ValI(Val{N2}()))
UpperTriangularMatrix(data::SVector{N2,T}) where {T,N2} = UpperTriangularMatrix(data.data,ValI(Val{N2}()))
StaticUpperTriangularMatrix(data::SVector{N2,T}) where {T,N2} = StaticUpperTriangularMatrix(data.data,ValI(Val{N2}()))
UpperTriangularMatrix(data::NTuple{N2,T}, ::Val{N}) where {T,N,N2} = UpperTriangularMatrix{T,N,N2}(data)
StaticUpperTriangularMatrix(data::NTuple{N2,T}, ::Val{N}) where {T,N,N2} = StaticUpperTriangularMatrix{T,N,N2}(data)
UpperTriangularMatrix(data::StaticUpperTriangularMatrix{T,N,N2}) where {T,N,N2} = UpperTriangularMatrix{T,N,N2}(data.data)
StaticUpperTriangularMatrix(data::UpperTriangularMatrix{T,N,N2}) where {T,N,N2} = StaticUpperTriangularMatrix{T,N,N2}(data.data)

LowerTriangularMatrix(data::NTuple{N2,T}) where {T,N2} = LowerTriangularMatrix(data,ValI(Val{N2}()))
StaticLowerTriangularMatrix(data::NTuple{N2,T}) where {T,N2} = StaticLowerTriangularMatrix(data,ValI(Val{N2}()))
LowerTriangularMatrix(data::SVector{N2,T}) where {T,N2} = LowerTriangularMatrix(data.data,ValI(Val{N2}()))
StaticLowerTriangularMatrix(data::SVector{N2,T}) where {T,N2} = StaticLowerTriangularMatrix(data.data,ValI(Val{N2}()))
LowerTriangularMatrix(data::NTuple{N2,T}, ::Val{N}) where {T,N,N2} = LowerTriangularMatrix{T,N,N2}(data)
StaticLowerTriangularMatrix(data::NTuple{N2,T}, ::Val{N}) where {T,N,N2} = StaticLowerTriangularMatrix{T,N,N2}(data)
LowerTriangularMatrix(data::StaticLowerTriangularMatrix{T,N,N2}) where {T,N,N2} = LowerTriangularMatrix{T,N,N2}(data.data)
StaticLowerTriangularMatrix(data::LowerTriangularMatrix{T,N,N2}) where {T,N,N2} = StaticLowerTriangularMatrix{T,N,N2}(data.data)

@inline Base.getindex(A::AbstractTriangularMatrix, i::Integer) = A.data[i]
@inline Base.getindex(A::PointerTriangularMatrix, i::Integer) = unsafe_load(A.data, i)

@inline function Base.getindex(A::StaticUpperTriangularMatrix{T,N,N2}, i::Integer, j::Integer) where {T,N,N2}
    @boundscheck begin
        max(i,j) > N && throw(BoundsError())
        if i > j
            return zero(T)
        end
    end
    @inbounds out = A.data[small_triangle(j)+i]
    out
end
@inline function Base.getindex(A::StaticLowerTriangularMatrix{T,N,N2}, i::Integer, j::Integer) where {T,N,N2}
    @boundscheck begin
        max(i,j) > N && throw(BoundsError())
        if i < j
            return zero(T)
        end
    end
    @inbounds out = A.data[small_triangle(i)+j]
    out
end
@inline function Base.getindex(A::UpperTriangularMatrix{T,N,N2}, i::Integer, j::Integer) where {T,N,N2}
    @boundscheck begin
        max(i,j) > N && throw(BoundsError())
        if i > j
            return zero(T)
        end
    end
    # @inbounds out = A.data[][small_triangle(j)+i]
    # out
    A.data[triangle_sub2ind(Val{N}(), i, j)]
end
@inline function Base.getindex(A::LowerTriangularMatrix{T,N,N2}, i::Integer, j::Integer) where {T,N,N2}
    @boundscheck begin
        max(i,j) > N && throw(BoundsError())
        if i < j
            return zero(T)
        end
    end
    # @inbounds out = A.data[][small_triangle(i)+j]
    # out
    # unsafe_load(Base.unsafe_convert(Ptr{T}, pointer_from_objref(A.data)), small_triangle(i)+j)
    A.data[triangle_sub2ind(Val{N}(), j, i)]
end
@inline function Base.setindex!(A::MutableTriangularMatrix{T,N,N2}, val, i::Integer) where  {T,N,N2}
    @boundscheck i > N2 && throw(BoundsError())
    unsafe_store!(point(A), convert(T, val), i)
    return val
end
@inline function Base.setindex!(A::MutableUpperTriangular{T,N,N2}, val, i::Integer, j::Integer) where  {T,N,N2}
    @boundscheck begin
        max(i,j) > N && throw(BoundsError())
        if i > j
            throw("Cannot set index within lower triangular of UpperTriangularMatrix.")
        end
    end
    unsafe_store!(point(A), convert(T, val), triangle_sub2ind(Val{N}(), i, j))
    return val
end
@inline function Base.setindex!(A::MutableLowerTriangular{T,N,N2}, val, i::Integer, j::Integer) where  {T,N,N2}
    @boundscheck begin
        max(i,j) > N && throw(BoundsError())
        if i < j
            throw("Cannot set index within upper triangular of LowerTriangularMatrix.")
        end
    end
    unsafe_store!(point(A), convert(T, val), triangle_sub2ind(Val{N}(), j, i))
    return val
end

"""
This always returns the static form so that it is non-allocating.
"""
LinearAlgebra.transpose(x::AbstractUpperTriangular{T,N,N2}) where {T,N,N2} = LowerTriangularMatrix{T,N,N2}(x.data)
LinearAlgebra.transpose(x::AbstractLowerTriangular{T,N,N2}) where {T,N,N2} = UpperTriangularMatrix{T,N,N2}(x.data)
@static if VERSION > v"0.6.9"
    """
    Taking the adjoint returns a static matrix so that it is non-allocating.
    """
    LinearAlgebra.adjoint(x::AbstractUpperTriangular{T,N,N2}) where {T<:Real,N,N2} = LowerTriangularMatrix{T,N,N2}(x.data)
    LinearAlgebra.adjoint(x::AbstractLowerTriangular{T,N,N2}) where {T<:Real,N,N2} = UpperTriangularMatrix{T,N,N2}(x.data)
end