text:
# maybe rewrite the 2d generator to the style of the 3d one
function generate_grid2D(cols::Int, rows::Int, v1::Vec{dim,T}, v2::Vec{dim,T}) where {dim,T <: Real}
    nodemat = Matrix{Node{dim,T}}(undef, (cols + 1, rows + 1))
    xstep = 1.0 / cols
    ystep = 1.0 / rows
    for i in 0:cols
        for j in 0:rows
            nodemat[i + 1,j + 1] = Node{dim,T}(i * xstep * v1 + j * ystep * v2)
        end
    end
    nodes = reshape(nodemat, length(nodemat))

    cellmat = Matrix{Element}(undef, (cols, rows))
    for i in 1:cols
        for j in 1:rows
            cellmat[i,j] = Element([nodemat[i, j] , nodemat[i + 1,j] ,nodemat[i + 1, j + 1], nodemat[i, j + 1]], nothing)
        end
    end
    cells = reshape(cellmat, length(cellmat))

    vertfacemat = Matrix{Element}(undef, (cols + 1, rows))
    for i in 2:cols
        for j in 1:rows
            vertfacemat[i,j] = attach_trivially!(cellmat[i - 1,j], cellmat[i,j])
        end
    end
    for j in 1:rows
        vertfacemat[1,j] = cellmat[1,j].boundaries[4]
        vertfacemat[cols + 1,j] = cellmat[cols,j].boundaries[2]
    end

    horfacemat = Matrix{Element}(undef, (cols, rows + 1))
    for i in 1:cols
        for j in 2:rows
            horfacemat[i,j] = attach_trivially!(cellmat[i,j - 1], cellmat[i,j])
        end
    end
    for i in 1:cols
        horfacemat[i,1] = cellmat[i,1].boundaries[1]
        horfacemat[i,rows + 1] = cellmat[i,rows].boundaries[3]
    end
    faces = vcat(reshape(vertfacemat, length(vertfacemat)), reshape(horfacemat, length(horfacemat)))

    cellsets = Dict{String,Set{Element}}()
    cellsets["all"] = Set(cells)
    nodesets = Dict{String,Set{Node}}()
    nodesets["corners"] = Set([nodemat[1,1], nodemat[cols + 1,1], nodemat[cols + 1,rows + 1],nodemat[1,rows + 1]])
    facesets = Dict{String,Set{Element}}()
    facesets["bottom"] = Set(horfacemat[:,1])
    facesets["left"] = Set(vertfacemat[1,:])
    facesets["top"] = Set(horfacemat[:,rows + 1])
    facesets["right"] = Set(vertfacemat[cols + 1,:])
    rootcells = cells
    grid = AdaptiveGrid{dim,T}(rootcells, cells, faces, nodes, cellsets, nodesets, facesets)
    return grid
end

function generate_grid3D(cols::Int, rows::Int, depth::Int, v1::Vec{dim,T}, v2::Vec{dim,T}, v3::Vec{dim,T}) where {dim,T <: Real}
    nodearr = Array{Node{dim,T}}(undef, (cols + 1, rows + 1, depth + 1))
    xstep = 1.0 / cols
    ystep = 1.0 / rows
    zstep = 1.0 / depth
    for i in 0:cols
        for j in 0:rows
            for k in 0:depth
                nodearr[i + 1,j + 1,k + 1] = Node{dim,T}(i * xstep * v1 + j * ystep * v2 + k * zstep * v3)
            end
        end
    end
    nodes = reshape(nodearr, length(nodearr))

    cellarr = Array{Element}(undef, (cols, rows, depth))
    cubenumberings = [[0, 0, 0], [1, 0, 0], [1,1,0], [0,1,0], [0, 0, 1], [1, 0, 1], [1,1,1], [0,1,1]]
    for i in 1:cols
        for j in 1:rows
            for k in 1:depth
                cellarr[i,j,k] = Element((x -> nodearr[CartesianIndex(Tuple([i,j,k] + x))]).(cubenumberings), nothing)
            end
        end
    end
    cells = reshape(cellarr, length(cellarr))
    for i in 1:cols
        for j in 1:rows
            for k in 1:depth
                if i < cols attach_trivially!(cellarr[i,j,k], cellarr[i + 1,j,k]) end
                if j < rows attach_trivially!(cellarr[i,j,k], cellarr[i,j + 1,k]) end
                if k < depth attach_trivially!(cellarr[i,j,k], cellarr[i,j,k + 1]) end
            end
        end
    end

    cellsets = Dict{String,Set{Element}}()
    cellsets["all"] = Set(cells)
    nodesets = Dict{String,Set{Node}}()
    nodesets["corners"] = Set((x -> nodearr[CartesianIndex(Tuple([cols + 1, rows + 1, depth + 1].^x))]).(cubenumberings))
    facesets = Dict{String,Set{Element}}()
    rootcells = cells
    grid = AdaptiveGrid{dim,T}(rootcells, cells, [], nodes, cellsets, nodesets, facesets)
    grid = update_leafmesh!(grid)
    grid.facesets["bottom"] = Set(filter(x -> (length(x.attachedto) == 0 && x.boundaryof.boundaries[1] == x), grid.faces)) # labels possibly misleading
    grid.facesets["left"] = Set(filter(x -> (length(x.attachedto) == 0 && x.boundaryof.boundaries[2] == x), grid.faces))
    grid.facesets["front"] = Set(filter(x -> (length(x.attachedto) == 0 && x.boundaryof.boundaries[3] == x), grid.faces))
    grid.facesets["right"] = Set(filter(x -> (length(x.attachedto) == 0 && x.boundaryof.boundaries[4] == x), grid.faces))
    grid.facesets["back"] = Set(filter(x -> (length(x.attachedto) == 0 && x.boundaryof.boundaries[5] == x), grid.faces))
    grid.facesets["top"] = Set(filter(x -> (length(x.attachedto) == 0 && x.boundaryof.boundaries[6] == x), grid.faces))
    return grid
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestConstraintSOCtoNonConvexQuad

using Test

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_RSOCtoNonConvexQuad()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config(
        exclude = Any[MOI.ConstraintDual, MOI.DualObjectiveValue],
    )
    @test MOI.Bridges.Constraint.RSOCtoNonConvexQuadBridge{Float64} ==
          MOI.Bridges.Constraint.concrete_bridge_type(
        MOI.Bridges.Constraint.RSOCtoNonConvexQuadBridge{Float64},
        MOI.VectorOfVariables,
        MOI.RotatedSecondOrderCone,
    )
    @test MOI.supports_constraint(
        MOI.Bridges.Constraint.RSOCtoNonConvexQuadBridge{Float64},
        MOI.VectorOfVariables,
        MOI.RotatedSecondOrderCone,
    )
    @test !MOI.supports_constraint(
        MOI.Bridges.Constraint.RSOCtoNonConvexQuadBridge{Float64},
        MOI.ScalarAffineFunction{Float64},
        MOI.RotatedSecondOrderCone,
    )
    bridged_mock = MOI.Bridges.Constraint.RSOCtoNonConvexQuad{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = ["test_basic_VectorOfVariables_RotatedSecondOrderCone"],
    )
    MOI.empty!(bridged_mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [0.5, 1.0, 1 / √2, 1 / √2])
    MOI.Test.test_conic_RotatedSecondOrderCone_VectorOfVariables(
        bridged_mock,
        config,
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                MOI.RotatedSecondOrderCone,
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        4,
        (
            (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
        ),
    )
    return
end

function test_SOCtoNonConvexQuad()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config(
        exclude = Any[MOI.ConstraintDual, MOI.DualObjectiveValue],
    )
    @test MOI.Bridges.Constraint.SOCtoNonConvexQuadBridge{Float64} ==
          MOI.Bridges.Constraint.concrete_bridge_type(
        MOI.Bridges.Constraint.SOCtoNonConvexQuadBridge{Float64},
        MOI.VectorOfVariables,
        MOI.SecondOrderCone,
    )
    @test MOI.supports_constraint(
        MOI.Bridges.Constraint.SOCtoNonConvexQuadBridge{Float64},
        MOI.VectorOfVariables,
        MOI.SecondOrderCone,
    )
    @test !MOI.supports_constraint(
        MOI.Bridges.Constraint.SOCtoNonConvexQuadBridge{Float64},
        MOI.ScalarAffineFunction{Float64},
        MOI.SecondOrderCone,
    )
    bridged_mock = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = ["test_basic_VectorOfVariables_SecondOrderCone"],
    )
    MOI.empty!(bridged_mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [1.0, 1 / √2, 1 / √2])
    MOI.Test.test_conic_SecondOrderCone_VectorOfVariables(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                MOI.SecondOrderCone,
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        3,
        (
            (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
        ),
    )
    return
end

end  # module

TestConstraintSOCtoNonConvexQuad.runtests()

ext: jl

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using DimensionalData, Test, Dates

# define dims with both long name and Type name
@dim Lon "Longitude"
@dim Lat "Latitude"

@testset "prettyprinting" begin

    timespan = DateTime(2001):Month(1):DateTime(2001,12)
    t = Ti([timespan...])
    x = Lon(Vector(0.5:10.0:359.5))
    y = Lat(Vector{Union{Float32, Missing}}(-89.5:10.0:89.5))
    n = Dim{:n}(Base.OneTo(10); mode=NoIndex())
    z = Z('a':'d')
    d = (x, y, z, t, n)
    A = DimArray(rand(length.(d)...), d; refdims=(Dim{:refdim}(1),))
    B = DimArray(rand(length(x)), (x,))
    C = DimArray(rand(Bool, 100, 100), (X, Y))

    for (d, str) in ((Ti(), "Ti"), (Lat(), "Lat"), (Lon(), "Lon"), (:n, ":n"), (Z(), "Z")) s1 = sprint(show, MIME("text/plain"), A)
        s2 = sprint(show, MIME("text/plain"), dims(A, d))
        s3 = sprint(show, MIME("text/plain"), dims(A, d))
        @test occursin("DimArray", s1)
        for s in (s1, s2, s3)
            @test occursin(str, s)
            @test occursin(str, s)
        end
    end

    # Test again but now with labelled array A
    A = DimArray(rand(length.(d)...), d; name=:test)
    s1 = sprint(show, MIME("text/plain"), A)
    @test occursin("test", s1)

    # Does it propagate after indexing?
    F = A[Ti(1:4)]
    s2 = sprint(show, MIME("text/plain"), F)
    @test occursin("test", s2)

    # Does it propagate after e.g. reducing operations?
    G = sum(A; dims = Ti)
    s3 = sprint(show, MIME("text/plain"), G)
    @test occursin("test", s3)

    # It should NOT propagate after binary operations
    B = DimArray(rand(length.(d)...), d; name=:test2)
    C = A .+ B
    s4 = sprint(show, MIME("text/plain"), C)
    @test !occursin("test", s4)

    # Test that broadcasted setindex! retains name
    D = DimArray(ones(length.(d)...), d; name=:olo)
    @. D = A + B
    s5 = sprint(show, MIME("text/plain"), D)
    @test occursin("olo", s5)

    # Test zero dim show
    D = DimArray([x for x in 1], (); name=:zero)
    sz = sprint(show, MIME("text/plain"), D)
    @test occursin("zero", sz)

    # Test vector show
    D = DimArray(ones(length(t)), t; name=:vec)
    sv = sprint(show, MIME("text/plain"), D)
    @test occursin("1", sv)
    @test occursin("vec", sv)

    # Test matrix show
    D = ones()
    D = DimArray(ones(X(5), Y(5)); name=:vec)
    sv = sprint(show, MIME("text/plain"), D)
    @test occursin("1", sv)

    sv = sprint(show, MIME("text/plain"), X())
    @test occursin("X", sv)

    @testset "show modes" begin
        sv = sprint(show, MIME("text/plain"), Categorical(Unordered()))
        @test occursin("Categorical", sv)
        @test occursin("Unordered", sv)
        sv = sprint(show, MIME("text/plain"), Sampled(Ordered(), Regular(), Points()))
        @test occursin("Sampled", sv)
        @test occursin("Ordered", sv)
        @test occursin("Regular", sv)
        @test occursin("Points", sv)
        sv = sprint(show, MIME("text/plain"), NoIndex())
        @test occursin("NoIndex", sv)
    end

end


text:
# This file was generated, do not modify it. # hide
@load LDA pkg=MultivariateStats
using Distances

clf = machine(LDA(dist=CosineDist()), X3, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

acc = mean(ŷ .== y[test])

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
using Random, Plots
include_ddp()

# Model
include_model("cartpole")
n, m, d = 4, 1, 4
model = Cartpole{Midpoint, FixedTime}(n, m, d, 1.0, 0.2, 0.5, 9.81)

function fd(model::Cartpole{Midpoint, FixedTime}, x, u, w, h, t)
    x + h * f(model, x + 0.5 * h * f(model, x, u, w), u, w)
end

function fd(model::Cartpole{Midpoint, FixedTime}, x, u, w, h, t)
    rz(z) = fd(model, z, x, u, w, h, t) # implicit Euler integration
	x⁺ = newton(rz, copy(x))
    # x⁺ = levenberg_marquardt(rz, copy(x))
	return x⁺
end

fd(model, zeros(model.n), 0.001 * ones(model.m), zeros(model.d), 0.1, 1)

function fdx(model::Cartpole{Midpoint, FixedTime}, x, u, w, h, t)
	rz(z) = fd(model, z, x, u, w, h, t) # implicit Euler integration
    x⁺ = newton(rz, copy(x))
    # x⁺ = levenberg_marquardt(rz, copy(x))
	∇rz = ForwardDiff.jacobian(rz, x⁺)
	rx(z) = fd(model, x⁺, z, u, w, h, t) # implicit Euler integration
	∇rx = ForwardDiff.jacobian(rx, x)
	return -1.0 * ∇rz \ ∇rx
end

fdx(model, ones(model.n), 0.01 * ones(model.m), zeros(model.d), 0.1, 1)

function fdu(model::Cartpole{Midpoint, FixedTime}, x, u, w, h, t)
	rz(z) = fd(model, z, x, u, w, h, t) # implicit Euler integration
    x⁺ = newton(rz, copy(x))
    # x⁺ = levenberg_marquardt(rz, copy(x))
	∇rz = ForwardDiff.jacobian(rz, x⁺)
	ru(z) = fd(model, x⁺, x, z, w, h, t) # implicit Euler integration
	∇ru = ForwardDiff.jacobian(ru, u)
	return -1.0 * ∇rz \ ∇ru
end

fdu(model, zeros(model.n), zeros(model.m), zeros(model.d), 0.1, 1)

# Time
T = 16
tf = 2.5
h = tf / (T - 1)

# Initial conditions, controls, disturbances
x1 = [0.0, 0.0, 0.0, 0.0]
xT = [0.0, π, 0.0, 0.0] # goal state
ū = [(t == 1 ? 1.0 : 0.0) * ones(m) for t = 1:T-1] # initialize no friciton model this direction (not sure why it goes the oppostive direction...)
w = [zeros(model.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model, x1, ū, w, h, T)
# x̄ = linear_interpolation(x1, xT, T)

# Objective
Q = [(t < T ? Diagonal(h * [1.0; 1.0; 1.0; 1.0])
        : Diagonal(1.0 * ones(model.n))) for t = 1:T]
q = [-2.0 * Q[t] * xT for t = 1:T]

R = [h * Diagonal(1.0 * ones(model.m)) for t = 1:T-1]
r = [zeros(model.m) for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

# Constraints
p_con = [t == T ? model.n : 2 * model.m for t = 1:T]

uL = -10.0 * ones(model.m)
uU = 10.0 * ones(model.m)
info_t = Dict(:uL => uL, :uU => uU, :inequality => collect(1:2 * model.m))
info_T = Dict(:xT => xT)

con_set = [StageConstraint(p_con[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	n = model.n
	m = model.m

	if t < T
		c[1:m] = cons.con[t].info[:uL] - u
		# c[m .+ (1:m)] = u - cons.con[t].info[:uU]
	end

	if t == T
		c[1:n] .= x - cons.con[T].info[:xT]
	end
end

# prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T,
#     analytical_dynamics_derivatives = true)
prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T)

# Solve
@time stats = constrained_ddp_solve!(prob,
    max_iter = 1000, max_al_iter = 10,
	ρ_init = 1.0, ρ_scale = 10.0,
    grad_tol = 1.0e-3,
	con_tol = 1.0e-3,
	cache = false)

@show sum([s[:iters] for s in stats[:stats_ilqr]])

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

# Visualize
plot(π * ones(T),
    width = 2.0, color = :black, linestyle = :dash)
plot!(hcat(x̄...)', width = 2.0, label = "",
	ylims = (-4.0, 8.0),
    color = :cyan)
plot(hcat(ū..., ū[end])',
    width = 2.0, linetype = :steppost)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, x, Δt = h)

# Simulate policy
include(joinpath(@__DIR__, "simulate.jl"))

# Model
model_sim = model
x1_sim = copy(x1)
T_sim = 10 * T

# Time
tf = h * (T - 1)
t = range(0, stop = tf, length = T)
t_sim = range(0, stop = tf, length = T_sim)
dt_sim = tf / (T_sim - 1)

# Policy
K = [K for K in prob.p_data.K]
# plot(vcat(K...))
# K = [K[end] for k in K]
# K = [prob.p_data.K[t] for t = 1:T-1]
# K, _ = tvlqr(model, x̄, ū, h, Q, R)
# # # K = [-k for k in K]
# K = [-K[t] for t = 1:T-1]
# plot(vcat(K...))

# Simulate
N_sim = 1
x_sim = []
u_sim = []
J_sim = []
Random.seed!(1)
for k = 1:N_sim
	wi_sim = 0.0 * rand(range(-1.0, stop = 1.0, length = 1000))

	println("sim: $k")#- w = $(wi_sim[1])")

	x_ddp, u_ddp, J_ddp, Jx_ddp, Ju_ddp = simulate_linear_feedback(
		model_sim,
		K,
	    x̄, ū,
		[xT for t = 1:T], [zeros(model.m) for t = 1:T-1],
		Q, R,
		T_sim, h,
		[wi_sim; 0.0; 0.0; 0.0],
		[zeros(model.d) for t = 1:T-1])

	push!(x_sim, x_ddp)
	push!(u_sim, u_ddp)
	push!(J_sim, J_ddp)
end

# Visualize
idx = (1:2)
plt = plot(t, hcat(x̄...)[idx, :]',
	width = 2.0, color = :black, label = "",
	xlabel = "time (s)", ylabel = "state",
	title = "LQR cartpole (J_avg = $(round(mean(J_sim), digits = 3)), N_sim = $N_sim)")

for xs in x_sim
	plt = plot!(t_sim, hcat(xs...)[idx, :]',
	    width = 2.0, color = :cyan, label = "")
end
display(plt)

plt = plot(t, hcat(ū..., ū[end])',
	width = 2.0, color = :black, label = "",
	xlabel = "time (s)", ylabel = "control",
	title = "cartpole (J_avg = $(round(mean(J_sim), digits = 3)), N_sim = $N_sim)")

for us in u_sim
	plt = plot!(t_sim, hcat(us..., us[end])',
		width = 1.0, color = :magenta, label = "",
		linetype = :steppost)
end
display(plt)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Survival, Documenter, StatsBase

makedocs(
    modules = [Survival],
    sitename = "Survival.jl",
    authors = "Alex Arslan",
    pages = [
        "Home" => "index.md",
        "Event Times" => "events.md",
        "Kaplan-Meier" => "km.md",
        "Nelson-Aalen" => "na.md",
        "Cox" => "cox.md",
    ],
)

deploydocs(
    repo = "github.com/JuliaStats/Survival.jl.git",
    target = "build",
)


text:
using DoubleFloats, LinearAlgebra
using Quadmath
using Test

using  Base: IEEEFloat
import Base: isapprox

using Random

const phi = Base.MathConstants.golden
const gamma = Base.MathConstants.eulergamma

const phi_df64   = Double64(phi)
const phi_df32   = Double32(phi)
const gamma_df64 = Double64(gamma)
const gamma_df32 = Double32(gamma)

const phi_df64hi = HI(phi_df64); const phi_df64lo = LO(phi_df64)
const phi_df32hi = HI(phi_df32); const phi_df32lo = LO(phi_df32)
const gamma_df64hi = HI(gamma_df64); const gamma_df64lo = LO(gamma_df64)
const gamma_df32hi = HI(gamma_df32); const gamma_df32lo = LO(gamma_df32)

function isapprox(a::T, b::T) where {T<:IEEEFloat}
   (iszero(a) && (0 < abs(b) < eps(T)^2) ||
    iszero(b) && (0 < abs(a) < eps(T)^2) ||
    isapprox(a, b, atol=eps(T)^2, rtol = eps(T)))
end

# isapprox(lo,eps(one(Float64)),atol=eps(one(Float64)))

# to cover functions that are not exported

macro df(func)
  :(DoubleFloats.$func)
end
macro df(func, arg)
  :(DoubleFloats.$func($arg))
end
macro df(func, arg1, arg2)
  :(DoubleFloats.$func($arg1, $arg2))
end
macro df(func, arg1, arg2, arg3)
  :(DoubleFloats.$func($arg1, $arg2, $arg3))
end
macro df(func, arg1, arg2, arg3, arg4)
  :(DoubleFloats.$func($arg1, $arg2, $arg3, $arg4))
end

include("specialvalues.jl")
include("construct.jl")
include("compare.jl")
include("promote.jl")

include("prearith.jl")
include("arithmetic.jl")
include("fma.jl")
include("modrem.jl")
include("functions.jl")
include("matmul.jl")
include("complex.jl")

include("linearalgebra.jl")
include("special_functions.jl")

include("concrete_accuracy.jl")
include("function_accuracy.jl")

include("coverage.jl")
include("op_.jl")

include("string.jl")
include("parse.jl")


text:
using Documenter
push!(LOAD_PATH, "..")
using MyTools
using MyTools.Ini
using MyTools.Tool
using MyTools.Atom
using MyTools.Nucleus

makedocs(
    sitename="MyTools",
    pages=[
        "index.md",
        "Base" => "base.md",
        "Ini" => "ini.md",
        "Tool" => "tool.md",
        "Atom" => "atom.md",
        "Nucleus" => "nucleus.md"
    ]
)


text:
"""
$(DocStringExtensions.README)
"""
module Symbolics

using TermInterface

using Metatheory

using DocStringExtensions

using LinearAlgebra

using Reexport

using DomainSets

import DomainSets: Domain
@reexport using SymbolicUtils

import TermInterface: similarterm, istree, operation, arguments, symtype

import SymbolicUtils: Term, Add, Mul, Pow, Sym,
                      FnType, @rule, Rewriters, substitute,
                      promote_symtype

import Metatheory.Rewriters: Chain, Prewalk, Postwalk, Fixpoint

import SymbolicUtils.Code: toexpr

import ArrayInterface

using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)

# re-export

export simplify, substitute

using SciMLBase, IfElse
export Num
using MacroTools
import MacroTools: splitdef, combinedef, postwalk, striplines
include("wrapper-types.jl")

include("num.jl")
include("complex.jl")

"""
    substitute(expr, s)

Performs the substitution on `expr` according to rule(s) `s`.
# Examples
```julia
julia> @variables t x y z(t)
4-element Vector{Num}:
    t
    x
    y
 z(t)
julia> ex = x + y + sin(z)
(x + y) + sin(z(t))
julia> substitute(ex, Dict([x => z, sin(z) => z^2]))
(z(t) + y) + (z(t) ^ 2)
```
"""
substitute

export Equation, ConstrainedEquation
include("equations.jl")

include("utils.jl")
export degree

using ConstructionBase
include("arrays.jl")

export @register, @register_symbolic
include("register.jl")

using TreeViews
export @variables, Variable
include("variable.jl")

include("linearity.jl")

using DiffRules, SpecialFunctions, NaNMath

using SparseArrays

export Differential, expand_derivatives

include("diff.jl")

export Difference, DiscreteUpdate

include("difference.jl")

export infimum, supremum
include("domains.jl")

export Integral
include("integral.jl")

include("array-lib.jl")

include("linear_algebra.jl")

import Libdl
include("build_function.jl")
export build_function

import Distributions
include("extra_functions.jl")

using Latexify
include("latexify_recipes.jl")

using RecipesBase
include("plot_recipes.jl")

using Requires

export symbolics_to_sympy
include("init.jl")

include("semipoly.jl")

# Hacks to make wrappers "nicer"
const NumberTypes = Union{AbstractFloat,Integer,Complex{<:AbstractFloat},Complex{<:Integer}}
(::Type{T})(x::SymbolicUtils.Symbolic) where {T<:NumberTypes} = throw(ArgumentError("Cannot convert Sym to $T since Sym is symbolic and $T is concrete. Use `substitute` to replace the symbolic unwraps."))
for T in [Num, Complex{Num}]
    @eval begin
        #(::Type{S})(x::$T) where {S<:Union{NumberTypes,AbstractArray}} = S(Symbolics.unwrap(x))::S

        SymbolicUtils.simplify(n::$T; kw...) = wrap(SymbolicUtils.simplify(unwrap(n); kw...))
        SymbolicUtils.simplify_fractions(n::$T; kw...) = wrap(SymbolicUtils.simplify_fractions(unwrap(n); kw...))
        SymbolicUtils.expand(n::$T) = wrap(SymbolicUtils.expand(unwrap(n)))
        substitute(expr::$T, s::Pair; kw...) = wrap(substituter(s)(unwrap(expr); kw...)) # backward compat
        substitute(expr::$T, s::Vector; kw...) = wrap(substituter(s)(unwrap(expr); kw...))
        substitute(expr::$T, s::Dict; kw...) = wrap(substituter(s)(unwrap(expr); kw...))

        SymbolicUtils.Code.toexpr(x::$T) = SymbolicUtils.Code.toexpr(unwrap(x))

        SymbolicUtils.setmetadata(x::$T, t, v) = wrap(SymbolicUtils.setmetadata(unwrap(x), t, v))
        SymbolicUtils.getmetadata(x::$T, t) = SymbolicUtils.getmetadata(unwrap(x), t)
        SymbolicUtils.hasmetadata(x::$T, t) = SymbolicUtils.hasmetadata(unwrap(x), t)

        Broadcast.broadcastable(x::$T) = x
    end
    for S in [:(Symbolic{<:FnType}), :CallWithMetadata]
        @eval (f::$S)(x::$T, y...) = wrap(f(unwrap(x), unwrap.(y)...))
    end
end

end # module


text:
@inline function GitHub.authenticate(config::Config)
    return GitHub.authenticate(config.github_token)
end

@inline function GitHub.authenticate(github_token::Credential)
    return GitHub.authenticate(github_token[])
end


text:
__precompile__()

module Convertible

using ItemGraphs
using MacroTools

import Base: ∘

export @convertible, @convert, Converter

const registry = ItemGraph{Any}();

abstract type AbstractConverter end
abstract type Converter{From,To} <: AbstractConverter end

origin(::Converter{From,To}) where {From,To} = From
target(::Converter{From,To}) where {From,To} = To

struct ComposedConverter{C1<:AbstractConverter,C2<:AbstractConverter} <: AbstractConverter
    conv1::C1
    conv2::C2
end

(conv::ComposedConverter)(args...) = conv.conv2(conv.conv1(args...))

origin(conv::ComposedConverter) = origin(conv.conv1)
target(conv::ComposedConverter) = target(conv.conv2)

function ∘(conv1::Converter{From,Via}, conv2::Converter{Via,To}) where {From,Via,To}
    ComposedConverter(conv1, conv2)
end

function ∘(conv1::Converter{From,Via1}, conv2::Converter{Via2,To}) where {From,Via1,Via2,To}
    throw(ArgumentError("The output type of the first converter is '$Via1' while the input type of the
            second converter is '$Via2'. These must match."))
end

function (::Type{T})(from, to, args...) where {T<:Converter}
    path = getpath(registry, from, to)
    _converter(T, args, path...)
end

@generated function _converter(T, args, path...)
    ex = :($T.parameters[1]($(path[1])(), $(path[2])(), args...))
    for i in eachindex(path[2:end-1])
        t1 = path[i+1]
        t2 = path[i+2]
        ex = :(ComposedConverter($ex, $T.parameters[1]($t1(), $t2(), args...)))
    end
    ex
end

function _convert(::Type{T}, obj::S) where {T,S}
    path = getpath(registry, S, T)
    _convert(obj, path[2:end]...)
end

@generated function _convert(obj, path...)
    ex = :(obj)
    for p in path
        ex = :($p.parameters[1]($ex))
    end
    ex
end

macro convert(expr::Expr)
    def = splitdef(expr)
    name = def[:name]
    args = def[:args]
    if length(args) == 1
        length(args) != 1 && throw(ArgumentError("Constructor must have a single argument."))
        _, origin, _ = splitarg(args[1])
        return quote
            add_edge!(registry, $(esc(origin)), $(esc(name)))
            $(esc(expr))
        end
    elseif length(args) >= 2
        _, origin, _ = splitarg(args[1])
        _, target, _ = splitarg(args[2])
        return quote
            if !($(esc(name)) <: Converter)
                name = string($(esc(name)))
                throw(ArgumentError("'$name' is not a subtype of 'Converter'."))
            end
            add_edge!(registry, $(esc(origin))(), $(esc(target))())
            $(esc(expr))
        end
    else
        throw(ArgumentError("@convert must be used on a single-argument or Converter constructor."))
    end
end

macro convertible(expr::Expr)
    @capture(expr, struct T_ fields__ end) || throw(ArgumentError("Expected a struct definition."))
    return quote
        $expr
        $(esc(T))(obj::S) where {S} = _convert($(esc(T)), obj)
        nothing
    end
end

end # module


text:
Base.String(::Type{T}) where T<:TruncatedKernel = "Truncated"
Base.String(::Type{T}) where T<:ParzenKernel = "Parzen"
Base.String(::Type{T}) where T<:TukeyHanningKernel = "Tukey-Hanning"
Base.String(::Type{T}) where T<:BartlettKernel = "Bartlett"
Base.String(::Type{T}) where T<:QuadraticSpectralKernel = "Quadratic Spectral"

Base.size(x::CovarianceMatrix) = size(x.V)
Base.getindex(x::CovarianceMatrix, row, col) = getindex(x.V, row, col)
Base.eltype(x::CovarianceMatrix) = eltype(x.V)

LinearAlgebra.inv(x::CovarianceMatrix) = inv(x.F)

## How should I call this terms? Breads? invfactor?
invfact(x::CovarianceMatrix{T1}, upper::Bool = true) where T1<:Cholesky = upper ? inv(x.F.U) : inv(x.F.L)

function invfact(x::CovarianceMatrix{T1}; regularize::Bool = true, atol::Real = 0.0, rtol::Real = (eps(real(float(one(eltype(x)))))*min(size(x)...))*iszero(atol)) where T1<:SVD
    Sr = similar(x.F.S)
    S  = x.F.S
    U  = x.F.Vt
    if regularize
        maxabs = abs.(first(S))
        tol = max(rtol*maxabs, atol)
        @inbounds for j in eachindex(S)
            if S[j] <= tol
                Sr[j] = zero(eltype(x))
            else
                Sr[j] = 1/sqrt(S[j])
            end
        end
    end
    Diagonal(Sr)*U
end

function quadinv(g::AbstractMatrix, x::CovarianceMatrix)
    LinearAlgebra.checksquare(g)
    @assert size(g, 1) == size(x, 1)
    B = invfact(x)
    gw = B*g
    dot(gw, gw)
end

function quadinv(g::AbstractVector, x::CovarianceMatrix)
    @assert length(g) == size(x, 1)
    B = invfact(x)
    gw = B*g
    dot(gw, gw)
end

maxeig(x::CovarianceMatrix{T1}) where {T1<:SVD} = first(x.F.S)
mineig(x::CovarianceMatrix{T1}) where {T1<:SVD} = last(x.F.S)
LinearAlgebra.Symmetric(x::CovarianceMatrix) = LinearAlgebra.Symmetric(x.V)

function LinearAlgebra.logdet(x::CovarianceMatrix{T1}) where {T1<:Cholesky}
    a = zero(eltype(x))
    @simd for j in 1:size(x, 1)
        @inbounds a += log.(x.F.U[j,j])
    end
    2*a
end

LinearAlgebra.logdet(x::CovarianceMatrix{T1}) where {T1<:SVD} = sum(log.(x.F.S))


text:
module PoissonTest

using Test
using ForneyLab

using ForneyLab: outboundType, isApplicable, unsafeMean, unsafeVar, slug, isProper, FactorNode, SoftFactor, Interface, FactorGraph, step!
using ForneyLab: VBPoissonOut, VBPoissonL, SPPoissonOutNP, SPPoissonLPN

@testset "Poisson ProbabilityDistribution construction" begin
    @test ProbabilityDistribution(Univariate, Poisson, l=2.0) == ProbabilityDistribution{Univariate, Poisson}(Dict(:l=>2.0))
    @test ProbabilityDistribution(Poisson, l=2.0) == ProbabilityDistribution{Univariate, Poisson}(Dict(:l=>2.0))
    @test ProbabilityDistribution(Poisson) == ProbabilityDistribution{Univariate, Poisson}(Dict(:l=>1.0))
    @test_throws MethodError ProbabilityDistribution(Multivariate, Poisson, l=2.0)
end

@testset "Poisson Message construction" begin
    @test Message(Univariate, Poisson, l=2.0) == Message{Poisson, Univariate}(ProbabilityDistribution{Univariate, Poisson}(Dict(:l=>2.0))    )
    @test Message(Poisson, l=2.0) == Message{Poisson, Univariate}(ProbabilityDistribution{Univariate, Poisson}(Dict(:l=>2.0))    )
    @test Message(Poisson) == Message{Poisson, Univariate}(ProbabilityDistribution{Univariate, Poisson}(Dict(:l=>1.0)))
    @test_throws MethodError Message(Multivariate, Poisson, l=2.0)
end

@testset "dims" begin
    @test dims(ProbabilityDistribution(Poisson)) == 1
end

@testset "slug" begin
    @test slug(Poisson) == "Poisson"
end

@testset "vague" begin
    @test vague(Poisson) == ProbabilityDistribution(Poisson, l=huge)
end

@testset "isProper" begin
    @test isProper(ProbabilityDistribution(Poisson, l=1.0))
    @test !isProper(ProbabilityDistribution(Poisson, l=0.0))
    @test !isProper(ProbabilityDistribution(Poisson, l=-1.0))
end

@testset "unsafe mean and variance" begin
    @test unsafeMean(ProbabilityDistribution(Poisson, l=2.0)) == 2.0
    @test unsafeVar(ProbabilityDistribution(Poisson, l=2.0)) == 2.0
end

@testset "log pdf" begin
    @test isapprox(logPdf(ProbabilityDistribution(Poisson, l=2.5),1), -1.583709268125845)
end

#-------------
# Update rules
#-------------

@testset "SPPoissonOutNP" begin
    @test SPPoissonOutNP <: SumProductRule{Poisson}
    @test outboundType(SPPoissonOutNP) == Message{Poisson}
    @test isApplicable(SPPoissonOutNP, [Nothing, Message{PointMass}])

    @test ruleSPPoissonOutNP(nothing, Message(Univariate, PointMass, m=2.0)) == Message(Poisson, l=2.0)
end

@testset "SPPoissonLPN" begin
    @test SPPoissonLPN <: SumProductRule{Poisson}
    @test outboundType(SPPoissonLPN) == Message{Gamma}
    @test isApplicable(SPPoissonLPN, [Message{PointMass}, Nothing])

    @test ruleSPPoissonLPN(Message(Univariate, PointMass, m=2.0), nothing) == Message(Gamma, a=3.0, b=1.0)
end

@testset "VBPoissonOut" begin
    @test VBPoissonOut <: NaiveVariationalRule{Poisson}
    @test outboundType(VBPoissonOut) == Message{Poisson}
    @test isApplicable(VBPoissonOut, [Nothing, ProbabilityDistribution])

    @test ruleVBPoissonOut(nothing, ProbabilityDistribution(Gamma, a=2.0, b=3.0)) == Message(Poisson, l=0.5)
end

@testset "VBPoissonL" begin
    @test VBPoissonL <: NaiveVariationalRule{Poisson}
    @test outboundType(VBPoissonL) == Message{Gamma}
    @test isApplicable(VBPoissonL, [ProbabilityDistribution, Nothing])

    @test ruleVBPoissonL(ProbabilityDistribution(Poisson, l=2.0), nothing) == Message(Gamma, a=3.0, b=1.0)
    @test ruleVBPoissonL(ProbabilityDistribution(Univariate, PointMass, m=2.0), nothing) == Message(Gamma, a=3.0, b=1.0)
end

@testset "averageEnergy and differentialEntropy" begin
    @test isapprox(differentialEntropy(ProbabilityDistribution(Poisson, l=1.0)), averageEnergy(Poisson, ProbabilityDistribution(Poisson, l=1.0), ProbabilityDistribution(Univariate, PointMass, m=1.0)))
    @test isapprox(differentialEntropy(ProbabilityDistribution(Poisson, l=10.0)), averageEnergy(Poisson, ProbabilityDistribution(Poisson, l=10.0), ProbabilityDistribution(Univariate, PointMass, m=10.0)))
    @test isapprox(differentialEntropy(ProbabilityDistribution(Poisson, l=100.0)), averageEnergy(Poisson, ProbabilityDistribution(Poisson, l=100.0), ProbabilityDistribution(Univariate, PointMass, m=100.0)))
end

end # module


text:
# Basic Fourier method without any kind of rest approximation

struct Fourier <: AbstractIteratedIntegralAlgorithm end

convorder(::Fourier) = 1//2
errcoeff(m, h, ::Fourier, ::MaxL2) = h*√3/(√2*π)
norv(m, n, ::Fourier) = 2*m*n

"""
    levyarea(W, n, alg::Fourier)

Simulates an approximation of the iterated Itô-integrals ``\\int_0^1W_s\\otimes dW_s``
of the given ``m``-dimensional increment of a Wiener process with step size 1.
The parameter ``n`` specifies the number of terms in the approximation and thus determines the accuracy.
This algorithm is based on a Fourier expansion of the Wiener process.
The algorithm needs approximately ``m^2+2\\cdot m\\cdot n`` Float's
and ``2\\cdot m\\cdot n`` random numbers.
The time complexity is ``\\mathcal{O}(m^2\\cdot n)``.
"""
function levyarea(W::AbstractVector{T}, n::Integer, alg::Fourier) where {T}
    rng = default_rng()
    m = length(W)
    X = randn(rng, T, n, m) # allocates m*n Floats
    Y = randn(rng, T, m, n) # allocates m*n Floats
    Y .= (Y .- √(T(2)).*W) ./ (1:n)'
    A = Y * X # allocates m*m Floats
    # Antisymmetrize
    # G = inv(2*T(π)).*(A .- A') # allocates m*m Floats
    for i = 1:m
        @inbounds A[i,i] = zero(T)
        for j = (i+1):m
            @inbounds A[i,j] = (A[i,j] - A[j,i])/(2*T(π))
            @inbounds A[j,i] = -A[i,j]
        end
    end
    return A
end

max_stars_repo_name: mrbauff/MasterThesis

max_issues_repo_name: mrbauff/MasterThesis


text:
"""
This struct contains all the data for the simulation to be run.
    case: CaseDefinition for creating and running the simulation.
    results_dir: Path to the directory where results are stored.
    iteration_year: Index of current iteration year.
    system_UC: SIIP PSY System for Unit Commitment.
    system_ED: SIIP PSY System for Economic Dispatch.
    zones: Names of modeled zones.
    lines: Inter-zonal ransmission lines modeled.
    hour_weight: Weight allocated to each hour of actual market clearing.
    peak_load: Annual system-wide peak load for capacity market clearing.
    markets: Boolean for selecting which markets are modeled.
    investors: Vectors of all investors created in the simulation.
    derating_data: Derating data foe each technology type.
    annual_growth: Annual growth rate data of different parameters.
"""

mutable struct AgentSimulation
    case::CaseDefinition
    results_dir::String
    iteration_year::Int64
    system_UC::Union{Nothing, PSY.System}
    system_ED::Union{Nothing, PSY.System}
    zones::Vector{String}
    lines::Vector{ZonalLine}
    hour_weight::Vector{Float64}
    peak_load::Float64
    markets::Dict{Symbol, Bool}
    investors::Vector{Investor}
    derating_data::DataFrames.DataFrame
    annual_growth::AxisArrays.AxisArray{Float64, 2}
end

get_case(sim::AgentSimulation) = sim.case
get_results_dir(sim::AgentSimulation) = sim.results_dir
get_iteration_year(sim::AgentSimulation) = sim.iteration_year
get_system_UC(sim::AgentSimulation) = sim.system_UC
get_system_ED(sim::AgentSimulation) = sim.system_ED
get_zones(sim::AgentSimulation) = sim.zones
get_lines(sim::AgentSimulation) = sim.lines
get_hour_weight(sim::AgentSimulation) = sim.hour_weight
get_peak_load(sim::AgentSimulation) = sim.peak_load
get_markets(sim::AgentSimulation) = sim.markets
get_investors(sim::AgentSimulation) = sim.investors
get_derating_data(sim::AgentSimulation) = sim.derating_data
get_annual_growth(sim::AgentSimulation) = sim.annual_growth

function set_iteration_year!(simulation::AgentSimulation, iteration_year::Int64)
    simulation.iteration_year = iteration_year
end

function set_peak_load!(simulation::AgentSimulation, peak_load::Float64)
    simulation.peak_load = peak_load
end

"""
This struct temporarily holds the data for creating a simulation.
"""
mutable struct AgentSimulationData
    case::CaseDefinition
    iteration_year::Int64
    system_UC::Union{Nothing, PSY.System}
    system_ED::Union{Nothing, PSY.System}
    zones::Vector{String}
    lines::Vector{ZonalLine}
    hour_weight::Vector{Float64}
    rep_hour_weight::Vector{Float64}
    peak_load::Float64
    markets::Dict{Symbol, Bool}
    investors::Union{Nothing, Vector{Investor}}
    queue_cost_data::DataFrames.DataFrame
    derating_data::DataFrames.DataFrame
    annual_growth::AxisArrays.AxisArray{Float64, 2}
end

function AgentSimulationData(case::CaseDefinition,
                        system_UC::Union{Nothing, PSY.System},
                        system_ED::Union{Nothing, PSY.System},
                        zones::Vector{String},
                        lines::Vector{ZonalLine},
                        hour_weight::Vector{Float64},
                        rep_hour_weight::Vector{Float64},
                        peak_load::Float64,
                        markets::Dict{Symbol, Bool},
                        queue_cost_data::DataFrames.DataFrame,
                        derating_data::DataFrames.DataFrame,
                        annual_growth::AxisArrays.AxisArray{Float64, 2})
    return AgentSimulationData(case,
                          1,
                          system_UC,
                          system_ED,
                          zones,
                          lines,
                          hour_weight,
                          rep_hour_weight,
                          peak_load,
                          markets,
                          nothing,
                          queue_cost_data,
                          derating_data,
                          annual_growth)
end

get_case(sim::AgentSimulationData) = sim.case
get_iteration_year(sim::AgentSimulationData) = sim.iteration_year
get_system_UC(sim::AgentSimulationData) = sim.system_UC
get_system_ED(sim::AgentSimulationData) = sim.system_ED
get_zones(sim::AgentSimulationData) = sim.zones
get_lines(sim::AgentSimulationData) = sim.lines
get_hour_weight(sim::AgentSimulationData) = sim.hour_weight
get_rep_hour_weight(sim::AgentSimulationData) = sim.rep_hour_weight
get_peak_load(sim::AgentSimulationData) = sim.peak_load
get_markets(sim::AgentSimulationData) = sim.markets
get_investors(sim::AgentSimulationData) = sim.investors
get_queue_cost_data(sim::AgentSimulationData) = sim.queue_cost_data
get_derating_data(sim::AgentSimulationData) = sim.derating_data
get_annual_growth(sim::AgentSimulationData) = sim.annual_growth

function set_investors!(simulation_data::AgentSimulationData,
                        investors::Vector{Investor})
    simulation_data.investors = investors
    return
end

function get_allprojects(simulation::Union{AgentSimulation, AgentSimulationData})
    projects = Project[]
    append!(projects, vcat(get_projects.(get_investors(simulation))...))
    return projects
end

# Get all investors' projects except those which have been retired.
function get_activeprojects(simulation::Union{AgentSimulation, AgentSimulationData})
    projects = get_allprojects(simulation)
    active_projects = filter(x -> !isa(x, Project{Retired}) , projects)
    return active_projects
end

# Get all investors' projects except options.
function get_decidedprojects(simulation::Union{AgentSimulation, AgentSimulationData})
    projects = get_allprojects(simulation)
    decided_projects = filter(x -> !isa(x, Project{Option}) , projects)
    return decided_projects
end


text:
abstract type AbstractHeap{T} <: AbstractTree{T} end

mutable struct BinaryHeap{T} <: AbstractHeap{T}
    tree::CompleteBinaryTree{T}
    size::Int
    maxsize::Int
    comparator::Function
    function BinaryHeap{T}(data, comparator) where T
        tree = CompleteBinaryTree{T}(data)
        size = tree.size
        maxsize = DEFAULT_MAX_SIZE
        n = new(tree, size, maxsize, comparator)
        buildheap!(n)
    end
end

BinaryHeap(data::Vector, comparator) = BinaryHeap{eltype(data)}(data, comparator)
MinHeap(data::Vector) = BinaryHeap(data, isless)
MaxHeap(data::Vector) = BinaryHeap(data, (x,y)->isless(y,x))

function _siftdown!(h::BinaryHeap, pos)
    let i = pos
        while !isleaf(h.tree, i)
            j = leftchild(h.tree, i)
            rc = rightchild(h.tree, i)
            # check if rightchild is within bounds
            if (rc != 0) && h.comparator(h.tree.data[rc], h.tree.data[j])
                j = rc
            end
            h.comparator(h.tree.data[i], h.tree.data[j]) && break
            swap!(h.tree, i, j)
            i = j
        end
    end
end

function _siftup!(h::BinaryHeap, pos)
    let i = pos
        while (i != 1) && h.comparator(h.tree.data[i], h.tree.data[parent(h.tree, i)])
            swap!(h.tree, i, parent(h.tree, i))
            i = parent(h.tree, i)
        end
    end
end

function insert!(h::BinaryHeap, item)
    (h.size < h.maxsize) || throw(ErrorException("heap capacity exceeded"))
    append!(h.tree, item)
    h.size += 1
    _siftup!(h, h.size)
end

function remove!(h::BinaryHeap)
    (h.size > 0) || throw(ArgumentError("no current element"))
    item = h.tree.data[1]
    swap!(h.tree, 1, h.size)
    h.size -= 1; h.tree.size -= 1
    if h.size != 0
        _siftdown!(h, 1)
    end
    return item
end

function remove!(h::BinaryHeap, pos)
    (1 <= pos <= h.size) || throw(BoundsError(h, pos))
    item = h.tree.data[pos]
    if pos == h.size
        h.size -= 1; h.tree.size -= 1
    else
        swap!(h.tree, pos, h.size)
        h.size -= 1; h.tree.size -= 1
        _siftup!(h, pos)
        if h.size != 0
            _siftdown!(h, pos)
        end
    end
    return item
end

function buildheap!(h::BinaryHeap)
    let i = h.size ÷ 2
        while i>0
            _siftdown!(h, i)
            i -= 1
        end
    end
    return h
end

peek(h::BinaryHeap) =  (h.size > 0) ? h.tree.data[1] : throw(ArgumentError("no current element"))
length(h::BinaryHeap) = h.size
isempty(h::BinaryHeap) = h.size == 0


text:
function inverseSide(side::HorizonSide)
    return HorizonSide(mod(Int(side)+2, 4))
end

function main(r::Robot)
    goToTheCorner(Sud, West)
    goToTheCorner(Nord, West)
    goToTheCorner(Nord, Ost)
    goToTheCorner(Sud, Ost)
end

function goToTheCorner(side1::HorizonSide, side2::HorizonSide)
    pathBack = []
    while isborder(r, side1) != true || isborder(r, side2) != true
        stepsToX = moveSide(r, side1)
        stepsToY = moveSide(r, side2)
        push!(pathBack, [inverseSide(side1), stepsToX])
        push!(pathBack, [inverseSide(side2), stepsToY])
    end
    putmarker!(r)
    for item in pathBack[end:-1:1]
        followByCount(r, item[1], item[2])
    end
end

function moveSide(r::Robot, side::HorizonSide, needMark=false)
    steps = 0
    while isborder(r, side) != true
        moveStep(r, side, needMark)
        steps += 1
    end
    return steps
end

function followByCount(r::Robot, side::HorizonSide, count::Int)
    for i = 1:count
        moveStep(r, side, false)
    end
end

function moveStep(r::Robot, side::HorizonSide, needMark::Bool)
    move!(r, side)
    if needMark
        putmarker!(r)
    end
end

max_stars_repo_name: TheDarkMoon/mirea-progs

max_issues_repo_name: TheDarkMoon/mirea-progs

max_forks_repo_name: TheDarkMoon/mirea-progs


text:
# Load Turing and MCMCChains.
using Turing
using LinearAlgebra
using CosmoCentral

w0waCDMCosmology = CosmoCentral.Flatw0waCDMCosmology()
println("Loaded central cosmology")

@model function gdemo(vecpCℓData, Cov, ConvolvedDensity, EuclidBias, EuclidIA,
    CosmologicalGrid)
    w₀ ~ Uniform(-1.5, -0.5)
    wₐ ~ Uniform(-0.5, 0.5)
    ΩM ~ Uniform(0.2, 0.4)
    ns ~ Uniform(0.8, 1.1)

    w0waCDMCosmology = CosmoCentral.Flatw0waCDMCosmology()
    w0waCDMCosmology.w0 = w₀
    w0waCDMCosmology.wa = wₐ
    w0waCDMCosmology.ΩM = ΩM
    w0waCDMCosmology.ns = ns
    CℓMCMC = CosmoCentral.EvaluateCℓMCMCStep(w0waCDMCosmology, ConvolvedDensity, EuclidBias,
    EuclidIA, CosmologicalGrid)

    ℓnumber = length(CℓMCMC.CℓArray[:,1,1])
    inumber = length(CℓMCMC.CℓArray[1,:,1])
    L = CosmoCentral.EliminationMatrix(inumber)

    vecpCℓMCMC = zeros(1,floor(Int,inumber*0.5*(inumber+1)))
    for ℓ in 1:ℓnumber
        vecCℓMCMC = vec(CℓMCMC.CℓArray[ℓ,:,:])
        vecpCℓMCMC = zeros(floor(Int,inumber*0.5*(inumber+1)))
        LinearAlgebra.mul!(vecpCℓMCMC, L, vecCℓMCMC)

        vecpCℓData[ℓ,:] ~ MvNormal(vecpCℓMCMC, Cov.Covariance[ℓ,:,:])
    end
end
MultipolesArrayTemp = CosmoCentral.LogSpaced(10.,3000., 101)
MultipolesArray = zeros(100)
#MultipolesWidths = vcat(CosmoCentral.Difference(MultipolesArrayTemp), ones(2000))
MultipolesWidths = CosmoCentral.Difference(MultipolesArrayTemp)
for i in 1:100
    MultipolesArray[i] = (MultipolesArrayTemp[i+1]+MultipolesArrayTemp[i])/2
end
CosmologicalGrid = CosmoCentral.CosmologicalGrid(ZArray = LinRange(0.001, 4., 500),
KArray = CosmoCentral.LogSpaced(1e-5, 50., 1000), ℓBinCenters = MultipolesArray,
ℓBinWidths = MultipolesWidths);

AnalitycalDensity = CosmoCentral.AnalitycalDensity()
CosmoCentral.NormalizeAnalitycalDensity!(AnalitycalDensity)

InstrumentResponse = CosmoCentral.InstrumentResponse()

ConvolvedDensity = CosmoCentral.ConvolvedDensity(DensityGridArray = ones(10, length(CosmologicalGrid.ZArray)))
CosmoCentral.NormalizeConvolvedDensity!(ConvolvedDensity, AnalitycalDensity, InstrumentResponse, CosmologicalGrid)

CosmoCentral.ComputeConvolvedDensityGrid!(CosmologicalGrid, ConvolvedDensity, AnalitycalDensity, InstrumentResponse)

EuclidBias = CosmoCentral.EuclidBias()
EuclidIA = CosmoCentral.ExtendedNLIA()

CℓData = CosmoCentral.EvaluateCℓMCMCStep(w0waCDMCosmology, ConvolvedDensity, EuclidBias,
EuclidIA, CosmologicalGrid)
#TODO now only LL, but this need to be more flexible...maybe list with probes?
Covaₗₘ = CosmoCentral.InstantiateEvaluateCovariance(CℓData, ConvolvedDensity, CosmologicalGrid, "Lensing",
"Lensing")
CovCℓ = CosmoCentral.InstantiateEvaluateCovariance(Covaₗₘ)
inumber = 10
L = CosmoCentral.EliminationMatrix(inumber)
vecpCℓData = zeros(100,floor(Int,inumber*0.5*(inumber+1)))
for ℓ in 1:length(CℓData.CℓArray[:,1,1])
    vecCℓData = vec(CℓData.CℓArray[ℓ,:,:])
    vecpCℓDataT = zeros(floor(Int,inumber*0.5*(inumber+1)))
    LinearAlgebra.mul!(vecpCℓDataT, L, vecCℓData)
    vecpCℓData[ℓ,:] = vecpCℓDataT
    Matrix = (CovCℓ.Covariance[ℓ,:,:] .+ transpose(CovCℓ.Covariance[ℓ,:,:])) ./2
    CovCℓ.Covariance[ℓ,:,:] = Matrix
end
println(size(CovCℓ.Covariance[1,:,:]))

model = gdemo(vecpCℓData, CovCℓ, ConvolvedDensity, EuclidBias, EuclidIA,
CosmologicalGrid)

chains = sample(model, MH(), 300; save_state = true)
write("first_chain-file.jls", chains)


text:
using Test
@testset "Testing solution to Exercise 2" begin

@testset "Running ex2.jl" begin
   include("../ex2.jl")
end;

@testset "Testing that the expected functions exist." begin
   @test @isdefined var_one_pass
   @test @isdefined var_two_pass
   @test @isdefined var_online
end;

@testset "Test that the values are close for a test dataset with N=10." begin
   jan1_2019 = 2458484.500000  # Julian Date for Jan 1, 2019
   y_test = generate_sample(10, jan1_2019)
   @test isapprox(var(y_test),var_one_pass(y_test),rtol=1e-3)
   @test !isapprox(var(y_test),var_one_pass(y_test),rtol=1e-4)
   @test isapprox(var(y_test),var_two_pass(y_test),rtol=1e-5)
   @test isapprox(var(y_test),var_online(y_test),rtol=1e-5)
   @test !isapprox(var(y_test),var_online(y_test),rtol=1e-11)
end;

@testset "Test that the values are close for a test dataset with N=10^4" begin
   jan1_2019 = 2458484.500000  # Julian Date for Jan 1, 2019
   y_test = generate_sample(10^4, jan1_2019)
   @test isapprox(var(y_test),var_one_pass(y_test),rtol=1e-1)
   @test !isapprox(var(y_test),var_one_pass(y_test),rtol=1e-2)
   @test isapprox(var(y_test),var_two_pass(y_test),rtol=1e-3)
   @test isapprox(var(y_test),var_online(y_test),rtol=1e-9)
   @test !isapprox(var(y_test),var_online(y_test),rtol=1e-11)
end;

@testset "Test that the values are close for a test dataset with N=10^6" begin
   jan1_2019 = 2458484.500000  # Julian Date for Jan 1, 2019
   y_test = generate_sample(10^6, jan1_2019)
   @test isapprox(var(y_test),var_one_pass(y_test),rtol=0.6)
   @test !isapprox(var(y_test),var_one_pass(y_test),rtol=0.1)
   @test isapprox(var(y_test),var_two_pass(y_test),rtol=1e-3)
   @test isapprox(var(y_test),var_online(y_test),rtol=1e-9)
   @test !isapprox(var(y_test),var_online(y_test),rtol=1e-11)
end;

end; # Exercise 2


text:
abstract type DiscreteFunctionGenerator end
abstract type AbstractResidueGenerator <: DiscreteFunctionGenerator end

const DFGen = DiscreteFunctionGenerator
const ARGen = AbstractResidueGenerator

# other useful functions #
Base.length(gen :: DFGen) = prod(big.(length.(gen.itr.iterators)))
# number of functions A → B
nfuns(dom :: FinSet, codom :: FinSet) = length(codom)^(length(dom))
domain(gen :: DFGen) = domain(gen.f)
codomain(gen :: DFGen) = codomain(gen.f)

# get n'th element of iterator
function Base.getindex(gen :: DFGen, i :: Int)
    return nth(gen, i)
end

# structures -- generators

"""
Generates all function of type
{0, …, from -1} → {0, …, to - 1}
"""
struct ResidueGenerator <: ARGen
    itr :: UnitRange{Int}
    f :: ResidueFunction
    from :: Int
    to :: Int
end

"""
generate all functions of the type:
ℤ_k1 × … × ℤ_km → ℤ_l
"""
struct ExtendedResidueGenerator{T, N} <: ARGen
    itr :: UnitRange{Int}
    f :: ExtendedResidueFunction{T, N}
    from :: NTuple{N, Int}
    to :: Int
end

"""
generate all product-type functions:
D1 × … × Dk → R1 × … × Rk
"""
struct FunProductGenerator{IT, N, T, D, R} <: DFGen
    itr :: IT
    f :: FunProduct{N, T, D, R}
end

"""
generate all tuple-type functions: (f1, ..., fn)
D1 × … × Dk → R1 × … × Rn
fi : D1 × … × Dk → Ri
"""
struct FunTuplingGenerator{IT, N, T, D, R} <: DFGen
    itr :: IT
    f :: FunTupling{N, T, D, R}
end

const ResGen = ResidueGenerator
const ExtGen = ExtendedResidueGenerator
const FTGen = FunTuplingGenerator
const FPGen = FunProductGenerator
const FXGens = Union{FPGen, FTGen}

# construction of generators #

# Residue generators #
# from two integers [0 .. from-1] → [0 .. to-1]
function ResGen(from :: T, to :: T) where T <: Integer
    itr = 0 : (to^from - 1)
    f = ResFun(from, to)
    return ResGen(itr, f, from, to)
end
#from two segments
ResGen(s1 :: T, s2 :: T) where T <: Segment = ResGen(length(s1), length(s2))

# Extended Residue generators #
# from tuple and number
function ExtGen(from :: NTuple{N}, to) where N
    num_from = prod(from)
    itr = 0 : (to^num_from - 1)
    f = ExtendedResidueFunction(from, to)
    return ExtGen(itr, f, from, to)
end
# from domain and codomain
ExtGen(dom, codom) where N = ExtGen(size(dom), length(codom))
# from domain and number
ExtGen(dom, to) where N = ExtGen(dom, segment(to))

# Function Product generators #
# from domain and codomain
function FPGen(dom :: AbstractProduct, codom :: AbstractProduct)
    n = ndims(dom)
    funs = Tuple(ResFun(dom[i], codom[i]) for i in 1 : n)
    tup = Tuple(nfuns(dom[i], codom[i]) for i in 1 : n)
    itr = Iterators.product(Segment.(tup)...)
    f = FunProduct(funs)
    return FPGen(itr, f)
end

# from two tuples
function FPGen(from :: NTuple{N, Int}, to :: NTuple{N, Int}) where N
    dom = DirectProduct(from)
    codom = DirectProduct(to)
    return FPGen(dom, codom)
end

# Function Tupling generators #
# from domain and codomain
function FTGen(dom :: AbstractProduct, codom :: AbstractProduct)
    from = size(dom)
    to = size(codom)
    n = ndims(codom)
    #funs = Tuple(ExtFun(dom, codom[i]) for i in 1 : n)
    tup = Tuple(nfuns(dom, codom[i]) for i in 1 : n)
    itr = product(Segment.(tup)...)
    f = FunTupling(from, to)
    return FTGen(itr, f)
end

# from two tuples
function FTGen(from :: NTuple{N, Int}, to :: NTuple{M, Int}) where {N, M}
    dom = DirectProduct(from)
    codom = DirectProduct(to)
    return FTGen(dom, codom)
end

function itstep(ctr, gen :: ARGen)
    value = ctr[1]
    state = ctr[2]
    fill_numbers!(gen.f.m, value, base = gen.to)
    return (gen.f, state)
end

function itstep(ctr, gen :: FPGen)
    value = ctr[1]
    state = ctr[2]
    mods = size(domain(gen.f))
    for (i, v) in enumerate(value)
        fi = proj(gen.f, i)
        fill_numbers!(fi.m, v, base = mods[i])
    end
    return (gen.f, state)
end

function itstep(ctr, gen :: FTGen)
    value = ctr[1]
    state = ctr[2]
    mods = size(domain(gen))
    for (i, v) in enumerate(value)
        fi = proj(gen.f, i)
        fill_numbers!(fi.m, v, base = mods[i])
    end
    return (gen.f, state)
end

itstep(:: Nothing, gen) = nothing

Base.iterate(gen :: ARGen) = itstep(iterate(gen.itr), gen)
Base.iterate(gen :: ARGen, state) = itstep(iterate(gen.itr, state), gen)
Base.iterate(gen :: FXGens) = itstep(iterate(gen.itr), gen)
Base.iterate(gen :: FXGens, state) = itstep(iterate(gen.itr, state), gen)

#=
itstep(n :: Nothing, gen :: FTGen) = nothing
itstep(n :: Nothing, gen :: FPGen) = nothing
=#

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
function slic(img, K::Int, M::Int, iters=20::Int)

	# Initialize variables
	imgLAB = Lab.(img)
	height, width = size(img)
	S = round(Int, (sqrt((height * width) / K)))
	clusters = Cluster[]
	labels = fill(-1, height, width)
	distance = fill(Inf, height, width)
	nPixels = Integer[]

	#################### Define methods ###########################
	function getGradient(y, x)
		x+1 > width && (x -= 2)
		y+1 > height && (y -= 2)

		return imgLAB[y+1, x+1].l - imgLAB[y, x].l +
               imgLAB[y+1, x+1].a - imgLAB[y, x].a +
               imgLAB[y+1, x+1].b - imgLAB[y, x].b
    end

	function clusterPixels()
        for i = 1:length(clusters)
            for x = (clusters[i].x - 2 * S):(clusters[i].x + 2 * S)
                (x <= 0 || x > width) && continue

				for y = (clusters[i].y - 2 * S):(clusters[i].y + 2 * S)
                    (y <= 0 || y > height) && continue

					# Distance in LAB space
                    Dc = sqrt((imgLAB[y, x].l - clusters[i].l)^2 +
                              (imgLAB[y, x].a - clusters[i].a)^2 +
                              (imgLAB[y, x].b - clusters[i].b)^2)

					# Distance in pixel space
                    Ds = sqrt((y - clusters[i].y)^2 +
                              (x - clusters[i].x)^2)

					# Scaled distance metric
                    D = sqrt((Dc / M)^2 + (Ds / S)^2)

                    if D < distance[y, x]
                        distance[y, x] = D
                        labels[y, x] = i
                    end
                end
            end
        end
    end

	function updateClusters()
        # Reset clusters
        for i = 1:length(clusters)
           clusters[i].y = clusters[i].x = nPixels[i] = 0
        end

		# Recalculate clusters
        for x in 1:width
            for y in 1:height
                labelIdx = labels[y, x]
                labelIdx == -1 && continue

                clusters[labelIdx].y += y
                clusters[labelIdx].x += x
                nPixels[labelIdx] += 1
            end
        end

		# Reassign clusters
        for i = 1:length(clusters)
            newY = div(clusters[i].y, nPixels[i])
            newX = div(clusters[i].x, nPixels[i])
            clusters[i].l = imgLAB[newY, newX].l
            clusters[i].a = imgLAB[newY, newX].a
            clusters[i].b = imgLAB[newY, newX].b
            clusters[i].y = newY
            clusters[i].x = newX
        end
    end

	function iterateClusters(cluster, i)
		if i == 1
			return cluster.l
		elseif i == 2
			return cluster.a
		elseif i == 3
			return cluster.b
		elseif i == 4
			return cluster.x
		elseif i == 5
			return cluster.y
		end
	end
	###################################################################
	# Initialize cluster
	for x = div(S, 2):S:width
        for y = div(S, 2):S:height
            push!(clusters, Cluster(imgLAB[y, x].l, imgLAB[y, x].a, imgLAB[y, x].b,
                                   y, x))
            push!(nPixels, 0)
        end
    end

	# Check gradient against surrounding pixels and reassign if lower
	for i = 1:length(clusters)
        grad = getGradient(clusters[i].y, clusters[i].x)

        for dy = -1:1
            for dx = -1:1
                _y = clusters[i].y + dy
                _x = clusters[i].x + dx
                newGrad = getGradient(_y, _x)
                if newGrad < grad
                    clusters[i].l = imgLAB[_y, _x].l
                    clusters[i].a = imgLAB[_y, _x].a
                    clusters[i].b = imgLAB[_y, _x].b
                    clusters[i].y = _y
                    clusters[i].x = _x
                    grad = newGrad
                end
            end
        end
    end

	# Run SLIC algorithm
	@time for i = 1:iters
        println("Running SLIC iteration $(i) ...")
        clusterPixels()
        updateClusters()
    end

	return clusters, labels
end

max_stars_repo_name: neuro-myoung/CellSegmenter

max_issues_repo_name: neuro-myoung/CellSegmenter

max_forks_repo_name: neuro-myoung/CellSegmenter


text:
using LeastSquaresOptim, Base.Test

# simple factor model
# only problem with "real" optimization
# nice example because J'J is not invertible
# but cholfact in sparse handles this case
function factor()
    name = "factor"
    function f!(x, fvec)
        fvec[1] = 3.0 - x[1] * x[4]
        fvec[2] = 2.0 - x[1] * x[5]
        fvec[3] = 5.0 - x[1] * x[6]

        fvec[4] = 4.5 - x[2] * x[4]
        fvec[5] = 3.2 - x[2] * x[5]
        fvec[6] = 2.0 - x[2] * x[6]

        fvec[7] = 5.0 - x[3] * x[4]
        fvec[8] = 1.3 - x[3] * x[5]
        fvec[9] = 1.5 - x[3] * x[6]
    end

    function g!(x, J)
        fill!(J, 0.0)
        J[1, 1] = -x[4]
        J[1, 4] = -x[1]
        J[2, 1] = -x[5]
        J[2, 5] = -x[1]
        J[3, 1] = -x[6]
        J[3, 6] = -x[1]

        J[4, 2] = -x[4]
        J[4, 4] = -x[2]
        J[5, 2] = -x[5]
        J[5, 5] = -x[2]
        J[6, 2] = -x[6]
        J[6, 6] = -x[2]

        J[7, 3] = -x[4]
        J[7, 4] = -x[3]
        J[8, 3] = -x[5]
        J[8, 5] = -x[3]
        J[9, 3] = -x[6]
        J[9, 6] = -x[3]
    end

    function g!(x, J::SparseMatrixCSC)
        Jvals = nonzeros(J)
        i = 0
        i += 1
        Jvals[i] = -x[4]
        i += 1
        Jvals[i] = -x[1]
        i += 1
        Jvals[i] = -x[5]
        i += 1
        Jvals[i] = -x[1]
        i += 1
        Jvals[i] = -x[6]
        i += 1
        Jvals[i] = -x[1]
        i += 1
        Jvals[i] = -x[4]
        i += 1
        Jvals[i] = -x[2]
        i += 1
        Jvals[i] = -x[5]
        i += 1
        Jvals[i] = -x[2]
        i += 1
        Jvals[i] = -x[6]
        i += 1
        Jvals[i] = -x[2]
        i += 1
        Jvals[i] = -x[4]
        i += 1
        Jvals[i] = -x[3]
        i += 1
        Jvals[i] = -x[5]
        i += 1
        Jvals[i] = -x[3]
        i += 1
        Jvals[i] = -x[6]
        i += 1
        Jvals[i] = -x[3]
    end
    x = ones(6)
    return name, f!, g!, x
end

iter = 0
for matrix in (:dense, :sparse)
    for (optimizer, optimizer_abbr) in ((LeastSquaresOptim.Dogleg(), :dl), (LeastSquaresOptim.LevenbergMarquardt(), :lm))
        factorization = matrix == :dense ? LeastSquaresOptim.QR() : LeastSquaresOptim.Cholesky()
        for (solver, solver_abbr) in ((factorization, :fact), (LeastSquaresOptim.LSMR(), :iter))
            iter += 1
            if matrix == :sparse && optimizer == LeastSquaresOptim.LevenbergMarquardt()
                continue
            else
                name, f!, g!, x = factor()
                fcur = ones(9)
                J = ones(9, 6)
                g!(x, J)
                if matrix == :sparse
                    J = sparse(J)
                end
                nls = LeastSquaresProblem(x = x, y = fcur, f! = f!, J = J, g! = g!)
                r = optimize!(nls, optimizer, solver)
                if iter == 1
                    show(r)
                end
                @printf("%-6s %4s %2s %30s %5d %5d   %5d   %10e\n", matrix, solver_abbr, optimizer_abbr, name, r.iterations, r.f_calls, r.g_calls, r.ssr)
                @test r.ssr <= 12
                @test r.converged
                end
            end
    end
end


text:
mutable struct GridCell
    w::Maybe{GtkWidget}
    r::Int
    c::Int
    GridCell(w::GtkWidget) = new(w, 1, 1)
    GridCell(w::GtkWidget, r::Int, c::Int) = new(w, r, c)
    GridCell(cell::GridCell) = cell
    GridCell(::Tuple) = new(missing, 1, 1)
end

Base.convert(::Type{GridCell}, x::T) where T = GridCell(x)

"""
        span(widget, r, c)

Creates a `GridCell` containing `widget` that spans `r` rows and `c` columns.

# Examples

```julia
julia> span(mybutton, 2, 2)
CellSpan(GtkButton..., 2, 2)
```
"""
span(w::GtkWidget, r::Int, c::Int) = [(x == y == 1 ? GridCell(w, r, c) : ()) for x=1:r, y=1:c]

span(::Tuple, r::Int, c::Int) = [() for x=1:r, y=1:c]

"""
        widget → N::Int

Returns `CellSpan` object for `widget` that spans `N` columns.

# Examples

```julia
julia> label"I span 4 columns" → 4
CellSpan(GtkSlider..., 1, 4)
```
"""
→(w::Union{GtkWidget, Tuple}, N::Int) = span(w, 1, N)
"""
        widget ↓ N::Int

Returns `CellSpan` object for `widget` that spans `N` rows.

# Examples

```julia
julia> Slider(1:100) ↓ 2
CellSpan(GtkSlider..., 2, 1)
```
"""
↓(w::Union{GtkWidget, Tuple}, N::Int) = span(w, N, 1)
→(M::Matrix, N::Int) = (w = M[1, 1]; span(w.w, w.r, w.c + N - 1))
↓(M::Matrix, N::Int) = (w = M[1, 1]; span(w.w, w.r + N - 1, w.c))
"""
        widget → r:c

Returns a `GridCell` object for `widget` that spans `r` rows and `c` columns.

# Examples

```julia
julia> Box((), :h) → 3:4
CellSpan(GtkBox..., 3, 4)
```
"""
→(w::Union{GtkWidget, Tuple}, a::UnitRange) = span(w, a.start, a.stop)

const Cell = Union{GtkWidget, GridCell, Tuple}
const Rows = Vector{Cell}

"""

        Grid(rows::Matrix)

Creates a new `GtkGrid` widget.

The `rows` parameter is a `Matrix` that contains the `GridCells`. A `GridCell` can
be either a `GtkWidget`, a `CellSpan` object or a empty tuple `()`.

The position of each element in the `rows` matrix will determine its position
on the `Grid`.

An `CellSpan` object represents a cell that spans multiple rows or/and columns.
An empty tuple indicates an empty cell.

# Examples

```julia
grid = Grid(@margin(10), [
    button  slider
    ()      label
])

grid = Grid(@margin(10), [
    # the button will span 2 rows
    CellSpan(button, 2, 1)  slider
    ()                      label
])

grid = Grid(@margin(10), [
    # the button will span 2 rows
    button ↓ 2  slider
    # the label will span 2 columns
    label → 2
])

grid = Grid(@margin(10)) do
    [
        # the textarea will span 3 rows and 2 columns
        textarea → 3:2  switch
        ()  ()            ()
        ()  ()            ()
        toolbar → 2     button
    ]
end
```
"""
function Grid(rows::Matrix{GridCell}; props...)
    grid = Gtk.GtkGrid()

    for (index, cell) in pairs(rows)
        if cell.w isa GtkWidget
            i, j = Tuple(index)
            r, c = cell.r, cell.c
            x = j:j + c - 1
            y = i:i + r - 1
            grid[x, y] = cell.w
        end
    end

    length(props) > 0 && set!(grid; props...)
    return grid
end

Grid(rows::Matrix{T}; props...) where {T} = Grid(GridCell.(rows); props...)

Grid(rows::Function; props...) = Grid(rows(); props...)


text:
# Constants corresponding to define in c code
const BLOSC_MIN_HEADER_LENGTH = 16 #!< Minimum header length (Blosc1)
const BLOSC_EXTENDED_HEADER_LENGTH = 32 #!< Extended header length (Blosc2, see README_HEADER)
const BLOSC_MAX_OVERHEAD = BLOSC_EXTENDED_HEADER_LENGTH #!< The maximum overhead during compression in bytes. This equals
                                                        #!< to @ref BLOSC_EXTENDED_HEADER_LENGTH now, but can be higher in future
                                                        #!< implementations.
const BLOSC_MAX_BUFFERSIZE = (typemax(Cint) - BLOSC_MAX_OVERHEAD) #!< Maximum source buffer size to be compressed
const BLOSC_MAX_TYPESIZE = 255 #!< Maximum typesize before considering source buffer as a stream of bytes.
                                #!< Cannot be larger than 255.
const BLOSC_MIN_BUFFERSIZE = 128 #!< Minimum buffer size to be compressed. Cannot be smaller than 66.

const BLOSC2_BDEFINED_FILTERS = 32 #!< Blosc-defined filters must be between 0 - 31.
const BLOSC2_REGISTERED_FILTERS = 160 # !< Blosc-registered filters must be between 32 - 159.
const BLOSC2_UDEFINED_FILTERS = 256 #!< User-defined filters must be between 128 - 255.
const BLOSC2_MAX_FILTERS = 6 #!< Maximum number of filters in the filter pipeline
const BLOSC2_MAX_UDFILTERS = 16 #!< Maximum number of filters that a user can register.

const BLOSC2_IO_FILESYSTEM = 0
const BLOSC_IO_LAST_BLOSC_DEFINED = 1
const BLOSC_IO_LAST_REGISTERED = 32

const BLOSC2_MAX_METALAYERS = 16
const BLOSC2_METALAYER_NAME_MAXLEN = 31

const BLOSC2_MAX_VLMETALAYERS = BLOSC2_MAX_METALAYERS
const BLOSC2_VLMETALAYERS_NAME_MAXLEN = BLOSC2_METALAYER_NAME_MAXLEN

"""
Codes for filters.
"""
const  BLOSC_NOFILTER = 0    #!< No filter.
const  BLOSC_SHUFFLE = 1     #!< Byte-wise shuffle.
const  BLOSC_BITSHUFFLE = 2  #!< Bit-wise shuffle.
const  BLOSC_DELTA = 3       #!< Delta filter.
const  BLOSC_TRUNC_PREC = 4  #!< Truncate precision filter.
const  BLOSC_LAST_FILTER = 5 #!< sentinel

"""
Codes for internal flags (see `blosc_cbuffer_metainfo`)
"""
const  BLOSC_DOSHUFFLE = 0x1     #!< byte-wise shuffle
const  BLOSC_MEMCPYED = 0x2      #!< plain copy
const  BLOSC_DOBITSHUFFLE = 0x4  #!< bit-wise shuffle
const  BLOSC_DODELTA = 0x8       #!< delta coding

"""
Codes for the different compressors shipped with Blosc
"""
const  BLOSC_BLOSCLZ = 0
const  BLOSC_LZ4 = 1
const  BLOSC_LZ4HC = 2
const  BLOSC_ZLIB = 4
const  BLOSC_ZSTD = 5
const  BLOSC_LAST_CODEC = 6 #!< Determine the last codec defined by Blosc.

"""
Codes for compression libraries shipped with Blosc (code must be < 8)
"""
const  BLOSC_BLOSCLZ_LIB = 0
const  BLOSC_LZ4_LIB = 1
const  BLOSC_ZLIB_LIB = 3
const  BLOSC_ZSTD_LIB = 4
const  BLOSC_UDCODEC_LIB = 6
const  BLOSC_SCHUNK_LIB = 7   #!< compressor library in super-chunk header

"""
  SplitMode
Split mode for blocks.
NEVER and ALWAYS are for experimenting with compression ratio.
AUTO for nearly optimal behaviour (based on heuristics).
FORWARD_COMPAT provides best forward compatibility (default).
"""
const  BLOSC_ALWAYS_SPLIT = 1
const  BLOSC_NEVER_SPLIT = 2
const  BLOSC_AUTO_SPLIT = 3
const  BLOSC_FORWARD_COMPAT_SPLIT = 4

"""
The codes for compressor formats shipped with Blosc
"""
const  BLOSC_BLOSCLZ_FORMAT = BLOSC_BLOSCLZ_LIB
const  BLOSC_LZ4_FORMAT = BLOSC_LZ4_LIB # LZ4HC and LZ4 share the same format
const  BLOSC_LZ4HC_FORMAT = BLOSC_LZ4_LIB
const  BLOSC_ZLIB_FORMAT = BLOSC_ZLIB_LIB,
const  BLOSC_ZSTD_FORMAT = BLOSC_ZSTD_LIB,
const  BLOSC_UDCODEC_FORMAT = BLOSC_UDCODEC_LIB


text:
# Use Recursion to find the first callback for type-stability

# Base Case: Only one callback
function find_first_continuous_callback(integrator, callback::AbstractContinuousCallback)
  (find_callback_time(integrator,callback)...,1,1)
end

# Starting Case: Compute on the first callback
function find_first_continuous_callback(integrator, callback::AbstractContinuousCallback, args...)
  find_first_continuous_callback(integrator,find_callback_time(integrator,callback)...,1,1,args...)
end

function find_first_continuous_callback(integrator,tmin::Number,upcrossing::Float64,
                                        event_occured::Bool,idx::Int,counter::Int,
                                        callback2)
  counter += 1 # counter is idx for callback2.
  tmin2,upcrossing2,event_occurred2 = find_callback_time(integrator,callback2)

  if event_occurred2 && (tmin2 < tmin || !event_occured)
    return tmin2,upcrossing2,true,counter,counter
  else
    return tmin,upcrossing,event_occured,idx,counter
  end
end

function find_first_continuous_callback(integrator,tmin::Number,upcrossing::Float64,event_occured::Bool,idx::Int,counter::Int,callback2,args...)
  find_first_continuous_callback(integrator,find_first_continuous_callback(integrator,tmin,upcrossing,event_occured,idx,counter,callback2)...,args...)
end

@inline function determine_event_occurance(integrator,callback)
  event_occurred = false
  Θs = linspace(typeof(integrator.t)(0),typeof(integrator.t)(1),callback.interp_points)
  interp_index = 0
  dt = integrator.t - integrator.tprev
  # Check if the event occured
  if typeof(callback.idxs) <: Void
    previous_condition = callback.condition(integrator.uprev,integrator.tprev,integrator)
  elseif typeof(callback.idxs) <: Number
    previous_condition = callback.condition(integrator.uprev[callback.idxs],integrator.tprev,integrator)
  else
    previous_condition = callback.condition(@view(integrator.uprev[callback.idxs]),integrator.tprev,integrator)
  end
  if integrator.event_last_time
    prev_sign = 0.0
  else
    prev_sign = sign(previous_condition)
  end
  prev_sign_index = 1
  if typeof(callback.idxs) <: Void
    next_sign = sign(callback.condition(integrator.u,integrator.t,integrator))
  elseif typeof(callback.idxs) <: Number
    next_sign = sign(callback.condition(integrator.u[callback.idxs],integrator.t,integrator))
  else
    next_sign = sign(callback.condition(@view(integrator.u[callback.idxs]),integrator.t,integrator))
  end
  if ((prev_sign<0 && !(typeof(callback.affect!)<:Void)) || (prev_sign>0 && !(typeof(callback.affect_neg!)<:Void))) && prev_sign*next_sign<=0
    event_occurred = true
    interp_index = callback.interp_points
  elseif callback.interp_points!=0  # Use the interpolants for safety checking
    tmp = integrator.tmp
    for i in 2:length(Θs)
      if !(typeof(callback.idxs) <: Number)
        integrator(tmp,integrator.tprev+dt*Θs[i])
        callback.idxs == nothing ? _tmp = tmp : _tmp = @view tmp[callback.idxs]
      else
        _tmp = integrator(integrator.tprev+dt*Θs[i])[callback.idxs]
      end
      new_sign = callback.condition(_tmp,integrator.tprev+dt*Θs[i],integrator)
      if prev_sign == 0
        prev_sign = sign(new_sign)
        prev_sign_index = i
      end
      if ((prev_sign<0 && !(typeof(callback.affect!)<:Void)) || (prev_sign>0 && !(typeof(callback.affect_neg!)<:Void))) && prev_sign*new_sign<0
        event_occurred = true
        interp_index = i
        break
      end
    end
  end
  event_occurred,interp_index,Θs,prev_sign,prev_sign_index
end

function find_callback_time(integrator,callback)
  event_occurred,interp_index,Θs,prev_sign,prev_sign_index = determine_event_occurance(integrator,callback)
  dt = integrator.t - integrator.tprev
  if event_occurred
    if typeof(callback.condition) <: Void
      new_t = zero(typeof(integrator.t))
    else
      if callback.interp_points!=0
        top_Θ = Θs[interp_index] # Top at the smallest
        bottom_θ = Θs[prev_sign_index]
      else
        top_Θ = typeof(integrator.t)(1)
        bottom_θ = typeof(integrator.t)(0)
      end
      if callback.rootfind
        tmp = integrator.tmp
        zero_func = (Θ) -> begin
          if !(typeof(callback.idxs) <: Number)
            integrator(tmp,integrator.tprev+Θ*dt)
            callback.idxs == nothing ? _tmp = tmp : _tmp = @view tmp[callback.idxs]
          else
            _tmp = integrator(integrator.tprev+Θ*dt)[callback.idxs]
          end
          out = callback.condition(tmp,integrator.tprev+Θ*dt,integrator)
          out
        end
        Θ = prevfloat(find_zero(zero_func,(bottom_θ,top_Θ),FalsePosition(),abstol = callback.abstol/10))
        #Θ = prevfloat(...)
        # prevfloat guerentees that the new time is either 1 floating point
        # numbers just before the event or directly at zero, but not after.
        # If there's a barrier which is never supposed to be crossed,
        # then this will ensure that
        # The item never leaves the domain. Otherwise Roots.jl can return
        # a float which is slightly after, making it out of the domain, causing
        # havoc.
        new_t = dt*Θ
      elseif interp_index != callback.interp_points
        new_t = dt*Θs[interp_index]
      else
        # If no solve and no interpolants, just use endpoint
        new_t = dt
      end
    end
  else
    new_t = zero(typeof(integrator.t))
  end
  new_t,prev_sign,event_occurred
end

function apply_callback!(integrator,callback::ContinuousCallback,cb_time,prev_sign)
  if cb_time != zero(typeof(integrator.t))
    integrator.t = integrator.tprev+cb_time
    integrator(integrator.u,integrator.t)
  end
  saved_in_cb = false

  @inbounds if callback.save_positions[1]
    savevalues!(integrator,true)
    saved_in_cb = true
  end

  integrator.u_modified = true

  if prev_sign < 0
    if typeof(callback.affect!) <: Void
      integrator.u_modified = false
    else
      callback.affect!(integrator)
    end
  elseif prev_sign > 0
    if typeof(callback.affect_neg!) <: Void
      integrator.u_modified = false
    else
      callback.affect_neg!(integrator)
    end
  end

  if integrator.u_modified
    @inbounds if callback.save_positions[2]
      savevalues!(integrator,true)
      saved_in_cb = true
    end
    return true,saved_in_cb
  end

  false,saved_in_cb
end

#Base Case: Just one
@inline function apply_discrete_callback!(integrator,callback::DiscreteCallback)
  saved_in_cb = false
  if callback.condition(integrator.u,integrator.t,integrator)
    @inbounds if callback.save_positions[1]
      savevalues!(integrator,true)
      saved_in_cb = true
    end
    integrator.u_modified = true
    callback.affect!(integrator)
    @inbounds if callback.save_positions[2]
      savevalues!(integrator,true)
      saved_in_cb = true
    end
  end
  integrator.u_modified,saved_in_cb
end

#Starting: Get bool from first and do next
@inline function apply_discrete_callback!(integrator,callback::DiscreteCallback,args...)
  apply_discrete_callback!(integrator,apply_discrete_callback!(integrator,callback)...,args...)
end

@inline function apply_discrete_callback!(integrator,discrete_modified::Bool,saved_in_cb::Bool,callback::DiscreteCallback,args...)
  bool,saved_in_cb2 = apply_discrete_callback!(integrator,apply_discrete_callback!(integrator,callback)...,args...)
  discrete_modified || bool, saved_in_cb || saved_in_cb2
end

@inline function apply_discrete_callback!(integrator,discrete_modified::Bool,saved_in_cb::Bool,callback::DiscreteCallback)
  bool,saved_in_cb2 = apply_discrete_callback!(integrator,callback)
  discrete_modified || bool, saved_in_cb || saved_in_cb2
end


text:
#=
Created on Tuesday 15 December 2020
Last update: -

@author: Bram De Jaegher
bram.de.jaegher@gmail.com

Launches the Pluto notebooks for the project
=#

cd("notebooks")
import Pkg; Pkg.activate(".");
Pkg.instantiate();
Pkg.add(url="https://github.com/Beramos/DS-Julia2925");

using Pluto
Pluto.run()


text:
using RecipesBase
using ArgCheck

function addline!(xs, ys, from, to)
    x_from, y_from = from
    x_to, y_to = to
    push!(xs, NaN, x_from, x_to, NaN)
    push!(ys, NaN, y_from, y_to, NaN)
    xs, ys
end

struct Grid
    edges
end

@recipe function plot(o::Grid)
    @argcheck length(o.edges) == 2
    xs = Float64[]
    ys = Float64[]
    (x_min, x_max), (y_min, y_max) = map(extrema, o.edges)
    s = 1.1
    xlims --> (s*x_min, s*x_max)
    ylims --> (s*y_min, s*y_max)
    for x in o.edges[1]
        from = [x, y_min]
        to = [x, y_max]
        addline!(xs, ys, from, to)
    end
    for y in o.edges[2]
        from = [x_min, y]
        to = [x_max, y]
        addline!(xs, ys, from, to)
    end
    xs, ys
end

@recipe function plot(o::EachTraversal{2})
    # grid
    legend --> :false
    @series begin
        Grid(o.edges)
    end

    # @series begin
    #     seriestype --> :scatter
    #     marker --> :o
    #     [o.position[1]], [o.position[2]]
    # end

    marker --> :x
    hits = collect(o)
    if isempty(hits)
        xs = [NaN]
        ys = [NaN]
    else
        xs, ys = map(1:2) do axis
            xs = Float64[]
            pos = o.position[axis]
            vel = o.velocity[axis]
            for hit in hits
                t = hit.entry_time
                push!(xs, pos + t*vel)
            end
            if !isempty(hits)
                t = last(hits).exit_time
                push!(xs, pos + t*vel)
            end
        end
    end
    xs, ys
end


text:
export fac, stirlingNumber

caching_fac = [1]
function fac(n)
    if n == 0
        return 1
    elseif size(caching_fac, 1) >= n
        return caching_fac[n]
    end
    fac = caching_fac[size(caching_fac, 1)]
    for i = size(caching_fac, 1)+1:n
        append!(caching_fac, fac * i)
    end
    return caching_fac[n]
end

function stirlingNumber(m, k)
    if k > m || k < 0 || (m > 0 && k == 0)
        return 0
    elseif m == 0 && k == 0
        return 1
    else
        return k * stirlingNumber(m - 1, k) + stirlingNumber(m - 1, k - 1)
    end
end

max_stars_repo_name: EllianCarlos/number-theory-but-in-julia

max_issues_repo_name: EllianCarlos/number-theory-but-in-julia

max_forks_repo_name: EllianCarlos/number-theory-but-in-julia


text:
using Black76
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "normcdf" begin
    z = normcdf(0.6)
    @test z > 0.5
end

@testset "normpdf" begin
    z = normcdf(0.6)
    @test z > 0.0
end

@testset "ntv" begin
    v = ntv(-0.6, 0.2)
    u = ntv( 0.6, 0.2)
    @test v == u  # ntv is symmetric in log-strike argument

    @test ntv(log(0), 0.2) == zero(Float64)  # no tv at infinite log-strike
end

@testset "black formulas" begin

    yf = 0.25
    df = exp(-0.08*yf)
    v0 = df*px(60.0/df, 65.0, yf, 1., 0.3);
    @test v0 ≈ 2.1334 atol = 5e-5  # Haug sample value

    round_trip(v) = impliedvol(100., 100., 1., -1., px(100., 100., 1., -1., v))
    # 20 vol roundtrip
    vol = 0.2
    rt_vol = round_trip(vol)
    @test rt_vol ≈ vol atol = 1e-14

    # 200 vol roundtrip
    vol = 2.
    rt_vol = round_trip(vol)
    @test rt_vol ≈ vol atol = 1e-14

    # 500 vol roundtrip
    vol = 5.
    rt_vol = round_trip(vol)
    @test rt_vol ≈ vol atol = 1e-14

    # zero premium returns zero vol
    @test impliedvol(100., 105., 1., 1., 0.) == zero(Float64)
end


text:
# Julia wrapper for header: MiniFB_enums.h
# Automatically generated using Clang.jl

# Julia wrapper for header: MiniFB.h
# Automatically generated using Clang.jl

function mfb_open(title, width, height)
    ccall((:mfb_open, libminifb), Ptr{mfb_window}, (Cstring, UInt32, UInt32), title, width, height)
end

function mfb_open_ex(title, width, height, flags)
    ccall((:mfb_open_ex, libminifb), Ptr{mfb_window}, (Cstring, UInt32, UInt32, UInt32), title, width, height, flags)
end

function mfb_update(window, buffer)
    ccall((:mfb_update, libminifb), mfb_update_state, (Ptr{mfb_window}, Ptr{Cvoid}), window, buffer)
end

function mfb_update_events(window)
    ccall((:mfb_update_events, libminifb), mfb_update_state, (Ptr{mfb_window},), window)
end

function mfb_close(window)
    ccall((:mfb_close, libminifb), Cvoid, (Ptr{mfb_window},), window)
end

function mfb_set_user_data(window, user_data)
    ccall((:mfb_set_user_data, libminifb), Cvoid, (Ptr{mfb_window}, Ptr{Cvoid}), window, user_data)
end

function mfb_get_user_data(window)
    ccall((:mfb_get_user_data, libminifb), Ptr{Cvoid}, (Ptr{mfb_window},), window)
end

function mfb_set_viewport(window, offset_x, offset_y, width, height)
    ccall((:mfb_set_viewport, libminifb), Bool, (Ptr{mfb_window}, UInt32, UInt32, UInt32, UInt32), window, offset_x, offset_y, width, height)
end

function mfb_set_active_callback(window, callback)
    ccall((:mfb_set_active_callback, libminifb), Cvoid, (Ptr{mfb_window}, mfb_active_func), window, callback)
end

function mfb_set_resize_callback(window, callback)
    ccall((:mfb_set_resize_callback, libminifb), Cvoid, (Ptr{mfb_window}, mfb_resize_func), window, callback)
end

function mfb_set_keyboard_callback(window, callback)
    ccall((:mfb_set_keyboard_callback, libminifb), Cvoid, (Ptr{mfb_window}, mfb_keyboard_func), window, callback)
end

function mfb_set_char_input_callback(window, callback)
    ccall((:mfb_set_char_input_callback, libminifb), Cvoid, (Ptr{mfb_window}, mfb_char_input_func), window, callback)
end

function mfb_set_mouse_button_callback(window, callback)
    ccall((:mfb_set_mouse_button_callback, libminifb), Cvoid, (Ptr{mfb_window}, mfb_mouse_button_func), window, callback)
end

function mfb_set_mouse_move_callback(window, callback)
    ccall((:mfb_set_mouse_move_callback, libminifb), Cvoid, (Ptr{mfb_window}, mfb_mouse_move_func), window, callback)
end

function mfb_set_mouse_scroll_callback(window, callback)
    ccall((:mfb_set_mouse_scroll_callback, libminifb), Cvoid, (Ptr{mfb_window}, mfb_mouse_scroll_func), window, callback)
end

function mfb_get_key_name(key)
    ccall((:mfb_get_key_name, libminifb), Cstring, (mfb_key,), key)
end

function mfb_is_window_active(window)
    ccall((:mfb_is_window_active, libminifb), Bool, (Ptr{mfb_window},), window)
end

function mfb_get_window_width(window)
    ccall((:mfb_get_window_width, libminifb), UInt32, (Ptr{mfb_window},), window)
end

function mfb_get_window_height(window)
    ccall((:mfb_get_window_height, libminifb), UInt32, (Ptr{mfb_window},), window)
end

function mfb_get_mouse_x(window)
    ccall((:mfb_get_mouse_x, libminifb), Cint, (Ptr{mfb_window},), window)
end

function mfb_get_mouse_y(window)
    ccall((:mfb_get_mouse_y, libminifb), Cint, (Ptr{mfb_window},), window)
end

function mfb_get_mouse_scroll_x(window)
    ccall((:mfb_get_mouse_scroll_x, libminifb), Cfloat, (Ptr{mfb_window},), window)
end

function mfb_get_mouse_scroll_y(window)
    ccall((:mfb_get_mouse_scroll_y, libminifb), Cfloat, (Ptr{mfb_window},), window)
end

function mfb_get_mouse_button_buffer(window)
    ccall((:mfb_get_mouse_button_buffer, libminifb), Ptr{UInt8}, (Ptr{mfb_window},), window)
end

function mfb_get_key_buffer(window)
    ccall((:mfb_get_key_buffer, libminifb), Ptr{UInt8}, (Ptr{mfb_window},), window)
end

function mfb_set_target_fps(fps)
    ccall((:mfb_set_target_fps, libminifb), Cvoid, (UInt32,), fps)
end

function mfb_wait_sync(window)
    ccall((:mfb_wait_sync, libminifb), Bool, (Ptr{mfb_window},), window)
end

function mfb_timer_create()
    ccall((:mfb_timer_create, libminifb), Ptr{mfb_timer}, ())
end

function mfb_timer_destroy(tmr)
    ccall((:mfb_timer_destroy, libminifb), Cvoid, (Ptr{mfb_timer},), tmr)
end

function mfb_timer_reset(tmr)
    ccall((:mfb_timer_reset, libminifb), Cvoid, (Ptr{mfb_timer},), tmr)
end

function mfb_timer_now(tmr)
    ccall((:mfb_timer_now, libminifb), Cdouble, (Ptr{mfb_timer},), tmr)
end

function mfb_timer_delta(tmr)
    ccall((:mfb_timer_delta, libminifb), Cdouble, (Ptr{mfb_timer},), tmr)
end

function mfb_timer_get_frequency()
    ccall((:mfb_timer_get_frequency, libminifb), Cdouble, ())
end

function mfb_timer_get_resolution()
    ccall((:mfb_timer_get_resolution, libminifb), Cdouble, ())
end


text:
using StaticRanges: grow_to, grow_to!

@testset "grow" begin
    @testset "grow_end" begin
        m,f,s = MutableRange(UnitRange(1, 10)), 1:10, static(UnitRange(1:10))
        x = @inferred(grow_end(m, 2))
        @test m == 1:10
        @test x == 1:12

        x = @inferred(grow_end(f, 2))
        @test f == 1:10
        @test x == 1:12

        x = @inferred((s -> grow_end(s, 2))(s))

        @test s == 1:10
        @test x == 1:12
    end

    @testset "grow_end!" begin
        m,f,s = MutableRange(UnitRange(1, 10)), 1:10, static(UnitRange(1, 10))
        x = @inferred(grow_end!(m, 2))
        @test m == 1:12
        @test x == 1:12

        #= FIXME These should have proper error messages
        x = @inferred(grow_end!(f, 2))
        @test f == 1:10
        @test x == 1:12

        x = @inferred((s -> grow_end(s, 2))(s))

        @test s == 1:10
        @test x == 1:12
        =#

        v = Vector{Int}(undef, 10)
        grow_end!(v, 2)
        @test length(v) == 12
    end

    @testset "grow_beg" begin
        m,f,s = mutable(UnitRange(1:10)), 1:10, static(UnitRange(1:10))
        x = @inferred(grow_beg(m, 2))
        @test m == 1:10
        @test x == -1:10

        x = @inferred(grow_beg(f, 2))
        @test f == 1:10
        @test x == -1:10

        x = @inferred((s -> grow_beg(s, 2))(s))
        @test s == 1:10
        @test x == -1:10
    end

    @testset "grow_beg!" begin
        m,f,s = MutableRange(UnitRange(1, 10)), 1:10, static(UnitRange(1, 10))
        x = @inferred(grow_beg!(m, 2))
        @test m == -1:10
        @test x == -1:10

        v = Vector{Int}(undef, 10)
        grow_beg!(v, 2)
        @test length(v) == 12
    end
end

@testset "shrink" begin
    @testset "shrink_beg" begin
        m,f,s = mutable(UnitRange(1:10)), 1:10, static(UnitRange(1:10))
        x = @inferred(StaticRanges.shrink_beg(m, 2))
        @test m == 1:10
        @test x == 3:10

        x = @inferred(shrink_beg(f, 2))
        @test f == 1:10
        @test x == 3:10

        x = @inferred((s -> shrink_beg(s, 2))(s))
        @test s == 1:10
        @test x == 3:10
    end

    @testset "shrink_beg!" begin
        m,f,s = mutable(UnitRange(1:10)), 1:10, static(UnitRange(1:10))
        x = @inferred(shrink_beg!(m, 2))
        @test m == 3:10
        @test x == 3:10

        v = Vector{Int}(undef, 10)
        shrink_beg!(v, 2)
        @test length(v) == 8
    end

    @testset "shrink_end" begin
        m,f,s = mutable(UnitRange(1:10)), 1:10, static(UnitRange(1:10))
        x = @inferred(shrink_end(m, 2))
        @test m == 1:10
        @test x == 1:8

        x = @inferred(shrink_end(f, 2))
        @test f == 1:10
        @test x == 1:8

        x = @inferred((s -> shrink_end(s, 2))(s))
        @test s == 1:10
        @test x == 1:8
    end

    @testset "shrink_end!" begin
        m,f,s = mutable(UnitRange(1:10)), 1:10, static(UnitRange(1:10))
        x = @inferred(shrink_end!(m, 2))
        @test m == 1:8
        @test x == 1:8

        v = Vector{Int}(undef, 10)
        shrink_end!(v, 2)
        @test length(v) == 8
    end
end

@testset "grow_to" begin
    m,f,s = mutable(UnitRange(1, 10)), 1:10, static(UnitRange(1, 10))
    x = @inferred(grow_to(m, 11))
    @test m == 1:10
    @test x == 1:11

    x = @inferred(grow_to(f, 11))
    @test f == 1:10
    @test x == 1:11

    x = (s -> grow_to(s, 11))(s)
    @test s == 1:10
    @test x == 1:11

    x = @inferred(grow_to!(m, 11))
    @test m == 1:11
    @test x == 1:11

    x = @inferred(grow_to!(m, 13))
    @test m == 1:13
    @test x == 1:13
end

#=
@testset "resize_first" begin
    m,f,s = UnitMRange(1, 10), 1:10, UnitSRange(1, 10)
    for (n, new_range) in ((11, 0:10), (9, 2:10), (10, 1:10))
        x = @inferred(resize_first(m, n))
        @test m == 1:10
        @test x == new_range

        x = @inferred(resize_first(f, n))
        @test f == 1:10
        @test x == new_range

        # FIXME This doesn't infer the type properly right now
        x = (s -> resize_first(s, n))(s)
        @test s == 1:10
        @test x == new_range
    end

    x = @inferred(resize_first!(m, 10))
    @test m == 1:10
    @test x == 1:10

    x = @inferred(resize_first!(m, 9))
    @test m == 2:10
    @test x == 2:10

    x = @inferred(resize_first!(m, 10))
    @test m == 1:10
    @test x == 1:10
end

@testset "[next/prev]_type" begin
    @test @inferred(StaticRanges.grow_end(["a"], 1)) == ["a", "b"]
    @test @inferred(StaticRanges.grow_beg(["a"], 1)) == ["`", "a"]
end
=#


text:
abstract type AbstractFreestream end

"""
    Freestream(α, β, Ω)
    Freestream(α, β, Ω_x, Ω_y, Ω_z)
    Freestream(U, Ω)
    Freestream(;
                 alpha = 0.,
                 beta  = 0.,
                 omega = [0,0,0]
               )

Define freestream conditions with angle of attack ``α`` (degrees), sideslip angle ``β`` (degrees), and (quasi-steady) rotation vector ``Ω`` for a vortex lattice analysis.

Alternatively, provide the velocity vector ``U``, which is normalized to determine the angles.
"""
struct Freestream{T} <: AbstractFreestream
    alpha :: T
    beta  :: T
    omega :: SVector{3,T}
end

Base.broadcastable(fs :: Freestream) = Ref(fs)

Freestream(α_deg, β_deg, Ω) = let T = promote_type(eltype(α_deg), eltype(β_deg), eltype(Ω)); Freestream{T}(deg2rad(α_deg), deg2rad(β_deg), Ω) end

Freestream(α_deg, β_deg, Ω_x, Ω_y, Ω_z) = Freestream(deg2rad(α_deg), deg2rad(β_deg), SVector(Ω_x, Ω_y, Ω_z))

Freestream(U, Ω) = let (V, α, β) = cartesian_to_freestream(normalize(U)); Freestream{T}(V, α, β, Ω) end

Freestream(; alpha = 0., beta = 0., omega = zeros(3)) = Freestream(alpha, beta, omega)

function Base.show(io :: IO, fs :: Freestream)
    println(io, "Freestream: ")
    for fname in fieldnames(typeof(fs))
        println(io, "    ", fname, " = ", getfield(fs, fname))
    end
end

"""
    velocity(freestream :: Freestream)

Compute the velocity of a `Freestream`.
"""
velocity(fs :: Freestream) = freestream_to_cartesian(1., fs.alpha, fs.beta)

"""
    body_frame_velocity(freestream :: Freestream)

Compute the velocity of Freestream in the body reference frame.
"""
body_frame_velocity(fs :: Freestream) = -velocity(fs)


text:
using Pkg

pkg"activate ."

packages = [
  "Agents",
  "AlgebraicPetri",
  "ApproxBayes",
  "BenchmarkTools",
  "BlackBoxOptim",
  "Bridge",
  "Catalyst",
  "Catlab",
  "DataFrames",
  "DiffEqCallbacks",
  "DiffEqJump",
  "DiffEqParamEstim",
  "DiffEqSensitivity",
  "DifferentialEquations",
  "Distances",
  "Distributions",
  "DrWatson",
  "DynamicalSystems",
  "ForwardDiff",
  "GpABC",
  "IJulia",
  "LabelledArrays",
  "Latexify",
  "LinearAlgebra",
  "MCMCChains",
  "ModelingToolkit",
  "NestedSamplers",
  "NLopt",
  "Optim",
  "OrdinaryDiffEq",
  "Petri",
  "PyCall",
  "Plots",
  "Random",
  "ResumableFunctions",
  "SimJulia",
  "SimpleDiffEq",
  "Soss",
  "SparseArrays",
  "StaticArrays",
  "StatsBase",
  "StatsPlots",
  "StochasticDiffEq",
  "Tables",
  "Turing",
  "Weave"
]

for p in packages
  Pkg.add(p)
end

pkg"instantiate"


text:
using Printf
using Profile
using CUDA
using BenchmarkTools

struct Point
    x::Float32 # meters
    y::Float32 # meters
    z::Float32 # meters
    vx::Float32 # meters/second
    vy::Float32 # meters/second
    vz::Float32 # meters/second
    mass::Float32 # kg
    numSprings::Int32
end

struct Spring
    k::Float32 # N/m
    p1::Int # Index of first point
    p2::Int # Index of second point
    l0::Float32 # meters
    dx::Float32
    dy::Float32
    dz::Float32
end

const kMaxSprings = 28
const kSpring = 500.0
const kGround = 100000.0
const kOscillationFrequency = 0
const kDropHeight = 0.2
const kNumSide = 60
const staticFriction = 0.5
const kineticFriction = 0.3
const dt = 0.0001
const dampening = 0.9995
const gravity = -9.81

function main(limit::Float64)
    t::Float64 = 0.0
    points::Array{Point} = []
    springs::Array{Spring} = []
    springIndices::Array{Int32} = Array{Int32}(undef, kNumSide * kNumSide * kNumSide * kMaxSprings)

    # Create the points
    for x = 1:kNumSide
        for y = 1:kNumSide
            for z = 1:kNumSide
                # (0,0,0) or (0.1,0.1,0.1) and all combinations
                p::Point = Point(x / 10.0, kDropHeight + y / 10.0, z / 10.0, 0, 0, 0, 0.1, 0)
                push!(points, p)
            end
        end
    end

    connected::Dict{Int32, Array{Int32}} = Dict()
    #Create the springs
    for x = 1:kNumSide
        for y = 1:kNumSide
            for z = 1:kNumSide
                p1index::Int = z + kNumSide * (y - 1) + kNumSide * kNumSide * (x - 1)
                if !haskey(connected, p1index)
                    connected[p1index] = []
                end
                p1::Point = points[p1index]
                for x1 = (x - 1):(x+1)
                    if x1 == kNumSide + 1 || x1 <= 0
                        continue
                    end
                    for y1 = (y - 1):(y+1)
                        if y1 == kNumSide + 1 || y1 <= 0
                            continue
                        end
                        for z1 = z:(z+1)
                            if z1 == kNumSide + 1 || z1 <= 0 || (x1 == x && y1 == y && z1 == z)
                                continue
                            end
                            p2index::Int = z1 + kNumSide * (y1 - 1) + kNumSide * kNumSide * (x1 - 1)
                            if !haskey(connected, p2index)
                                connected[p2index] = []
                            elseif p2index in connected[p1index]
                                continue
                            end
                            push!(connected[p1index], p2index)
                            push!(connected[p2index], p1index)
                            p2::Point = points[p2index]
                            length::Float32 = sqrt(abs2(p1.x - p2.x) + abs2(p1.y - p2.y) + abs2(p1.z - p2.z))
                            points[p1index] = Point(p1.x, p1.y, p1.z, p1.vx, p1.vy, p1.vz, p1.mass, p1.numSprings + 1)
                            p1 = points[p1index]
                            points[p2index] = Point(p2.x, p2.y, p2.z, p2.vx, p2.vy, p2.vz, p2.mass, p2.numSprings + 1)
                            p2 = points[p2index]
                            push!(springs, Spring(kSpring, p1index, p2index, length, 0, 0, 0))
                            currentNumSprings::Int32 = size(springs)[1]
                            springIndices[(p1index - 1) * kMaxSprings + p1.numSprings] = currentNumSprings
                            springIndices[(p2index - 1) * kMaxSprings + p2.numSprings] = currentNumSprings
                        end
                    end
                end
            end
        end
    end

    numSprings::Int32 = length(springs)
    numPoints::Int32 = length(points)

    springThreads::Int64 = 25
    springBlocks::Int64 = ceil(numSprings / springThreads)
    pointThreads::Int64 = 12
    pointBlocks::Int64 = ceil(numPoints / pointThreads)

    points_d = CuArray(points)
    springs_d = CuArray(springs)
    springIndices_d = CuArray(springIndices)

    println("num springs evaluated: ", numSprings)
    println("time multiplier: ",  limit / dt)

    println("Warming cache")
    @cuda threads=springThreads blocks=springBlocks updateSprings(points_d, springs_d, Float32(0.0), numSprings)
    @cuda threads=pointThreads blocks=pointBlocks updatePoints(points_d, springs_d, springIndices_d, numPoints)
    CUDA.@time @cuda threads=springThreads blocks=springBlocks updateSprings(points_d, springs_d, Float32(0.0), numSprings)
    CUDA.@time @cuda threads=pointThreads blocks=pointBlocks updatePoints(points_d, springs_d, springIndices_d, numPoints)
    @cuda threads=pointThreads blocks=pointBlocks updatePoints(points_d, springs_d, springIndices_d, numPoints)
    points_d = CuArray(points)
    springs_d = CuArray(springs)
    springIndices_d = CuArray(springIndices)
    synchronize()
    println("Starting")

    @time begin
        while t < limit
            adjust::Float32 = 1 + sin(t * kOscillationFrequency) * 0.1

            @cuda threads=springThreads blocks=springBlocks updateSprings(points_d, springs_d, adjust, numSprings)
            @cuda threads=pointThreads blocks=pointBlocks updatePoints(points_d, springs_d, springIndices_d, numPoints)

            t += dt
        end
        synchronize()
    end
    ps::Array{Point} = Array(points_d)
    println(ps[1].y)
    #for p in ps
    #    println(p.x, ",", p.y, ",", p.z)
    #end
end

function updateSprings(points::CuDeviceArray{Point}, springs::CuDeviceArray{Spring}, adjust::Float32, n::Int32)
    i::Int32 = (blockIdx().x-1) * blockDim().x + threadIdx().x
    if i > n
        return
    end
    l::Spring = springs[i]
    p1index::Int32 = l.p1
    p2index::Int32 = l.p2

    p1::Point = points[p1index]
    p2::Point = points[p2index]

    l0::Float32 = l.l0
    k::Float32 = l.k
    px::Float32 = p1.x - p2.x
    py::Float32 = p1.y - p2.y
    pz::Float32 = p1.z - p2.z

    dist::Float32 = sqrt(px * px + py * py + pz * pz)

    # negative if repelling, positive if attracting
    f::Float32 = k * (dist - (l0 * adjust))
    # distribute force across the axes
    fdist::Float32 = f / dist
    dx::Float32 = fdist * px
    dy::Float32 = fdist * py
    dz::Float32 = fdist * pz

    springs[i] = Spring(k, p1index, p2index, l0, dx, dy, dz)

    return nothing
end

function updatePoints(points::CuDeviceArray{Point}, springs::CuDeviceArray{Spring}, springIndices::CuDeviceArray{Int32}, n::Int32)
    i::Int32 = (blockIdx().x-1) * blockDim().x + threadIdx().x
    if i > n
        return
    end
    p::Point = points[i]

    fy = gravity * p.mass
    fx = 0
    fz = 0

    iTMS::Int32 = (i - 1) * kMaxSprings

    for j in 1:p.numSprings
        s::Spring = springs[springIndices[iTMS + j]]

        if s.p1 == i
            fx -= s.dx
            fy -= s.dy
            fz -= s.dz
        else
            fx += s.dx
            fy += s.dy
            fz += s.dz
        end
    end

    if p.y < 0
        fy += -kGround * p.y
        fh::Float32 = sqrt(fx * fx + fz * fz)
        if fh < abs(fy * staticFriction)
            fx = 0
            fz = 0
        else
            fyfric = fy * kineticFriction
            fx = fx - fyfric
            fz = fz - fyfric
        end
    end
    ax::Float32 = fx / p.mass
    ay::Float32 = fy / p.mass + gravity
    az::Float32 = fz / p.mass
    vx::Float32 = (ax * dt + p.vx) * dampening
    vy::Float32 = (ay * dt + p.vy) * dampening
    vz::Float32 = (az * dt + p.vz) * dampening
    x::Float32 = p.x + p.vx * dt
    y::Float32 = p.y + p.vy * dt
    z::Float32 = p.z + p.vz * dt
    points[i] = Point(x, y, z, vx, vy, vz, p.mass, p.numSprings)

    return nothing
end

main(0.1)
#@profile main()

max_stars_repo_name: unionfinder/PrimordialOozebot

max_issues_repo_name: unionfinder/PrimordialOozebot

max_forks_repo_name: unionfinder/PrimordialOozebot