text:
"""
    E = kepler(M, e; atol=1e-12, maxit=1000)

Solves [Kepler's equation](https://en.wikipedia.org/wiki/Kepler%27s_equation) for eccentric anomaly or
hyperbolic anomaly using Newton's method.
"""
function kepler(M::Real, e::Real; atol = 1e-12, maxit = 1000)
    if e == 1           # Parabolic case
        throw(ArgumentError("Parabolic case (e=1) has not been implemented"))
    elseif e < 1        # Elliptical case
        # # from Seppo Mikkola (slower, but better initial guess)
        # α = (1-e)/(4*e+0.5)
        # β = 0.5*M/(4*e+0.5)
        # z = copysign((abs(β) + √(β^2 + α^3))^(1/3), β)
        # s = z - α/z
        # s = s - 0.078*s^5/(1+e)
        # E = M + e*(3*s - 4*s^3)
        E = M
        # iterate until convergence
        err = 1+atol
        it = 1
        while (err > atol) && (it < maxit)
            it = it + 1
            Ep = E
            errp = err
            E = E - (E - e*sin(E) - M)/(1-e*cos(E))       # Newton's method
            # f = E - e*sin(E) - M
            # df = 1 - e*cos(E)
            # ddf = e*sin(E)
            # E = E - 2*f*df/(2*df^2 - f*ddf)               # Halley's method (slower, no accuracy improvement)
            err = abs(E-Ep)
            # if error gets worse, use bisection (makes up for poor initial guess)
            if err > errp
                E = (E+Ep)/2
            end
        end
        if it == maxit
            @warn "Kepler solver failed to converge: err = $err"
        end
        return E
    else                # Hyperbolic case
        if abs(M) > 4π
            H = sign(M) * 4π
        else
            H = M
        end
        err = 1+atol
        it = 1
        while (err > atol) && (it < maxit)
            it = it + 1
            Hp = H
            errp = err
            H = H + (M - e*sinh(H) + H)/(e*cosh(H) - 1)     # Newton's method
            err = abs(H-Hp)
            # if error gets worse, use bisection (makes up for poor initial guess)
            if err > errp
                H = (H+Hp)/2
            end
        end
        if it == maxit
            @warn "Kepler solver failed to converge: err = $err"
        end
        return H
    end
end


text:
module kuramoto_model

using DifferentialEquations, Flux, Optim, DiffEqFlux, Plots, Random

rng=MersenneTwister(123) #make the results consistent across runs
N=10000
ω=randn(rng,N)
norder=2 #OA=1, versus Hannay Ansatz=2

#ω0,γ,K,β
model_params=[2π, 0.1,0.8,0.1]
ic=rand(rng,N)*2π;

tend=20.0
tstart=0.0
sampling=0.05

function full_model(dϕ,ϕ,p,t)
    ω0, γ, K, β=p
    ψ,R=OrderParameter(ϕ)
    for i in 1:N
        dϕ[i]=ω0+γ*ω[i]+R*K*sin(ψ-ϕ[i]+β)
    end
end

function reduced_model(du,u,p,t)
    ω0, γ, K, β=p
    ψ,R=u
    n=norder
    du[1]=ω0+K/2.0*sin(β)*(1+R^(2*norder))
    du[2]=-γ*R+K/2.0*cos(β)*R*(1-R^(2*norder))
end

function OrderParameter(ϕ;m=1.0)
    Z=1.0/N*sum(exp.(im.*ϕ*m))
    R=abs(Z)
    ψ=angle(Z)
    return ψ,R
end

function OrderParameterComplex(ϕ;m=1.0)
    Z=1.0/N*sum(exp.(im.*ϕ*m))
    return Z
end

icr=[angle(OrderParameterComplex(ic)), abs(OrderParameterComplex(ic))]

function InitPlot()

    lay = @layout [ a{0.6w} [b; c]]

    t=tstart:0.01:tend
    sol=solve(ODEProblem(full_model,ic,(tstart,tend),model_params), Tsit5(), saveat=t)

    num_tps=length(t)
    R=zeros(num_tps)
    ψ=zeros(num_tps)
    for tp in 1:num_tps
        psi, rr=OrderParameter(sol[:,tp])
        R[tp]=rr
        ψ[tp]=psi
    end

    pl1=plot(ψ,R,lw=2, legend=false, color=:blue, proj = :polar)

    #radius plot
    pl2=plot(t,R,lw=2, legend=false, color=:blue)

    #angle plot
    pl3=plot(t,sin.(ψ), lw=2, legend=false, color=:blue)

    pl=plot(pl1, pl2, pl3, layout=lay)
    display(pl)

end

function predict_adjoint(param) # Our 1-layer neural network
    prob=ODEProblem(reduced_model,icr,(tstart,tend), model_params)
    Array(concrete_solve(prob,Tsit5(),icr,param,saveat=tstart:sampling:tend,abstol=1e-8,reltol=1e-6))
end

function generate_data(ic; param=model_params)
    prob=ODEProblem(full_model,ic,(tstart,tend), param)
    tps=tstart:sampling:tend
    num_tps=length(tps)
    d1=Array(concrete_solve(prob,Tsit5(),ic,param,saveat=tps))
    # Get a time series of the order parameter
    d2=mapslices(OrderParameterComplex, d1, dims=1)
    data=Float64[angle.(d2); abs.(d2)]
    data=reshape(data,2,num_tps)
    return data
end

# Generate some data to fit, and add some noise to it
data=generate_data(ic)
#σN=0.05
#data+=σN*randn(size(data))

function loss_adjoint(param)
    prediction = predict_adjoint(param)
    loss=0.0
    for k in 1:length(prediction[1,:])
        Rp=prediction[2,k]
        Rd=data[2,k]
        θp=prediction[1,k]
        θd=data[1,k]
        loss+=sqrt((Rp*cos(θp)-Rd*cos(θd))^2+(Rp*sin(θp)-Rd*sin(θd))^2)
    end
    loss,prediction
end

function train_model(;pguess=[6.1, 0.15,1.1,0.15])
    println("The initial loss function is $(loss_adjoint(pguess)[1])")
    res = DiffEqFlux.sciml_train(loss_adjoint,pguess,BFGS(initial_stepnorm=0.0001), maxiters=4000)
    println("The parameters are $(res.minimizer) with final loss value $(res.minimum)")
    return(res)
end

function plotFit(param; tend=20.0)

    validationPlot(param,ic; tend=tend)

end

function validationPlot(param,ic; tend=10.0)

    lay = @layout [ a{0.6w} [b; c]]

    t=tstart:0.01:tend
    icr=[angle(OrderParameterComplex(ic)), abs(OrderParameterComplex(ic))]

    sol=solve(ODEProblem(full_model,ic,(tstart,tend),model_params), Tsit5(), saveat=t)
    sol_fit=solve(ODEProblem(reduced_model,icr,(tstart,tend),param), Tsit5(), saveat=t)

    num_tps=length(t)
    R=zeros(num_tps)
    ψ=zeros(num_tps)
    for tp in 1:num_tps
        psi, rr=OrderParameter(sol[:,tp])
        R[tp]=rr
        ψ[tp]=psi
    end

    pl1=plot(ψ,R,lw=2, legend=false, color=:blue, proj = :polar)
    plot!(pl1, sol_fit, color=:green, proj=:polar, vars=(1,2))

    #radius plot
    pl2=plot(t,R,lw=2, legend=false, color=:blue)
    plot!(pl2, sol_fit, color=:green, lw=2.0, vars=(0,2))

    #angle plot
    pl3=plot(t,sin.(ψ), lw=2, legend=false, color=:blue)
    plot!(pl3, sol_fit.t, sin.(sol_fit[1,:]), color=:green)

    pl=plot(pl1, pl2, pl3, layout=lay)
    display(pl)

end

end #module

#d=kuramoto_model.generate_data(randn(kuramoto_model.N))
#kuramoto_model.predict_adjoint(kuramoto_model.model_params)

kuramoto_model.InitPlot()
results_kuramoto=kuramoto_model.train_model()
kuramoto_model.plotFit(results_kuramoto.minimizer)
kuramoto_model.validationPlot(results_kuramoto.minimizer, π*rand(kuramoto_model.N))
kuramoto_model.validationPlot(results_kuramoto.minimizer, 1.5*π*rand(kuramoto_model.N))

ext: jl

max_stars_repo_name: khannay/FittingParamsDiffEqFlux

max_issues_repo_name: khannay/FittingParamsDiffEqFlux

max_forks_repo_name: khannay/FittingParamsDiffEqFlux


text:
@testset "JDF.jl jdf_str" begin
df = DataFrame([collect(1:100) for i =1:3000])
df[!, :int_missing] =
    rand(
        [rand(rand([UInt, Int, Float64, Float32, Bool])), missing],
        nrow(df)
    )

df[!, :missing] .= missing
df[!, :strs] = [randstring(8) for i in 1:nrow(df)]
df[!, :stringarray] = StringVector([randstring(8) for i in 1:nrow(df)])

df[!, :strs_missing] = [rand([missing, randstring(8)]) for i in 1:nrow(df)]
df[!, :stringarray_missing] = StringVector([rand([missing, randstring(8)]) for i in 1:nrow(df)])
df[!, :symbol_missing] = [rand([missing, Symbol(randstring(8))]) for i in 1:nrow(df)]
df[!, :char] = getindex.(df[!, :strs], 1)
df[!, :char_missing] = allowmissing(df[!, :char])
df[rand(1:nrow(df), 10), :char_missing] .= missing

savejdf("a.jdf", df)

df2 = loadjdf(jdf"a.jdf")

@test ncol(df2) == 3009
@test nrow(df2) == 100

@test all(all(isequal(df[!,n], df2[!, n])) for n in names(df))

df3 = loadjdf(jdf"a.jdf", cols = [:missing, :strs])

@test ncol(df3) == 2
@test nrow(df3) == 100
# clean up
rm("a.jdf", force=true, recursive=true)
end


text:
module Containers

using PyCall: PyObject

using DockerPy: Collection, docker, @pyinterface
using DockerPy.Images: Image

export Container, DockerPath,
    attach, attach_socket, exec_run, pause, reload, rename, restart, start, stats, stop

mutable struct Container
    o::PyObject
end
Container() = Container(docker.models.Container())

const ContainerCollection = Collection{Container}

Container(x::ContainerCollection, image::Image; kwargs...) =
    Container(PyObject(x).create(PyObject(image); kwargs...))

ExecResult(exit_code, output) = (exit_code = exit_code, output = output)
ExecResult(x) = ExecResult(x...)

struct DockerPath
    container::Container
    path::AbstractString
end

"""
    attach(container::Container; stdout = true, stderr = true, stream = false, logs = false, demux = false)

Attach to this container.

# Arguments
- `stdout::Bool`: include stdout.
- `stderr::Bool`: include stderr.
- `stream::Bool`: return container output progressively as an iterator of strings, rather than a single string.
- `logs::Bool`: include the container’s previous output.
"""
attach(x::Container; kwargs...) = PyObject(x).attach(kwargs...)
"""
    attach_socket(container::Container, params = nothing, ws = false)

Like `attach`, but returns the underlying socket-like object for the HTTP request.
"""
attach_socket(x::Container; kwargs...) = PyObject(x).attach_socket(kwargs...)
function exec_run(
    x::Container,
    cmd;
    stdout = true,
    stderr = true,
    stdin = false,
    tty = false,
    privileged = false,
    user = "",
    detach = false,
    stream = false,
    socket = false,
    environment = nothing,
    workdir = nothing,
    demux = false,
)
    exec_result = ExecResult(PyObject(x).exec_run(
        cmd,
        stdout,
        stderr,
        stdin,
        tty,
        privileged,
        user,
        detach,
        stream,
        socket,
        environment,
        workdir,
        demux,
    ))
    if stream || socket || detach
        @warn "`detach`, `stream` or `socket` is `true`! Unable to determine if there is an error!"
    else
        if !iszero(exec_result.exit_code)
            @error "the exit code is not zero! Something wrong happened!"
        end
    end
    return exec_result
end # function exec_run
pause(x::Container) = PyObject(x).pause()
reload(x::Container) = PyObject(x).reload()
rename(x::Container, name::AbstractString) = PyObject(x).rename(name)
restart(x::Container; timeout = nothing) = PyObject(x).restart(timeout = timeout)
start(x::Container) = PyObject(x).start()
function start(f::Function, x::Container; remove = false, kwargs...)
    start(x)
    try
        f(x; kwargs...)
    finally
        stop(x)
        if remove
            rm(x)  # TODO: so slow
        end
    end
end # function start
stats(x::Container; decode::Bool = false, stream::Bool = true) =
    PyObject(x).stats(decode = decode, stream = stream)
stop(x::Container, timeout = nothing) = PyObject(x).stop(timeout = timeout)
Base.kill(x::Container, signal = nothing) = PyObject(x).kill(signal)
Base.rm(x::Container; kwargs...) = PyObject(x).remove(kwargs...)
Base.show(io::IO, x::Container) =
    print(io, "Container(name = \"$(x.name)\", short_id = \"$(x.short_id)\")")

Base.collect(x::ContainerCollection) = map(Container, PyObject(x).list())
Base.get(x::ContainerCollection, container_id) = Container(PyObject(x).get(container_id))
Base.empty!(x::ContainerCollection, filters = nothing) =
    PyObject(x).prune(filters = filters)

@pyinterface Container

end # module Containers


text:
###########################################################################################
#                           Hierarchical Navigable Small World                            #
###########################################################################################
mutable struct HierarchicalNSW{T, F, V, M}
    lgraph::LayeredGraph{T}
    added::Vector{Bool}
    data::V
    ep::T
    entry_level::Int
    vlp::VisitedListPool
    metric::M
    efConstruction::Int #size of dynamic candidate list
    ef::Int
end

###########################################################################################
#                                 Creating Graphs / Struct                                #
###########################################################################################
function HierarchicalNSW(data;
        metric=Euclidean(),
        M = 10, #5 to 48
        M0 = 2M,
        m_L = 1/log(M),
        efConstruction = 100,
        ef = 10,
        max_elements=length(data))
    T = max_elements <= typemax(UInt32) ? UInt32 : UInt64
    lg = LayeredGraph{T}(max_elements, M0, M, m_L)
    ep = T(0)
    F = eltype(data[1])
    vlp = VisitedListPool(1,max_elements)
    return HierarchicalNSW{T,F,typeof(data),typeof(metric)}(
        lg, fill(false, max_elements), data, ep, 0, vlp, metric, efConstruction, ef)
end

"""
    add_to_graph!(notify_func, hnsw, indices, multithreading=false)

Add `i ∈ indices` referring to `data[i]` into the graph.

`notify_func(i)` provides an interface for a progress notification by current index.

Indices already added previously will be ignored.
"""
function add_to_graph!(notify_func::Function,
                       hnsw::HierarchicalNSW{T}, indices) where {T}
    any(hnsw.added[indices]) && @warn "Some of the points have already been added!"
    for i ∈ indices
        hnsw.added[i] || insert_point!(hnsw, T(i))
        hnsw.added[i] = true
        notify_func(i)
    end
    hnsw
end

"""
    add_to_graph!(hnsw, indices)

short form of `add_to_graph!(notify_func, hnsw, indices)`
"""
add_to_graph!(hnsw::HierarchicalNSW{T}, indices) where {T} =
    add_to_graph!(identity, hnsw, indices)

add_to_graph!(notify_func::Function, hnsw::HierarchicalNSW; kwargs...) =
    add_to_graph!(notify_func, hnsw, eachindex(hnsw.data); kwargs...)

add_to_graph!(hnsw::HierarchicalNSW; kwargs...) =
    add_to_graph!(identity, hnsw; kwargs...)

set_ef!(hnsw::HierarchicalNSW, ef) = hnsw.ef = ef

###########################################################################################
#                                    Utility Functions                                    #
###########################################################################################
get_enter_point(hnsw::HierarchicalNSW) = hnsw.ep
function set_enter_point!(hnsw::HierarchicalNSW, ep)
    hnsw.ep = ep
    hnsw.entry_level = levelof(hnsw.lgraph, ep)
end
get_entry_level(hnsw::HierarchicalNSW) = hnsw.entry_level

@inline distance(hnsw, i, j) = @inbounds evaluate(hnsw.metric, hnsw.data[i], hnsw.data[j])
@inline distance(hnsw, i, q::AbstractVector) = @inbounds evaluate(hnsw.metric, hnsw.data[i], q)
@inline distance(hnsw, q::AbstractVector, j) = @inbounds evaluate(hnsw.metric, hnsw.data[j], q)

function Base.show(io::IO, hnsw::HierarchicalNSW)
    lg = hnsw.lgraph
    maxpoints = length(lg.linklist)
    addedpoints = count(isassigned.(Ref(lg.linklist), 1:maxpoints))
    println(io, "Hierarchical Navigable Small World with $(get_entry_level(hnsw)) layers")
    println(io, "$addedpoints of $maxpoints have been added to the index")
    println(io, "Index parameters are:")
    println(io, "\t M0 = $(lg.M0), M = $(lg.M), m_L ≈ $(round(lg.m_L,digits=1))")
    println(io, "\t metric = $(hnsw.metric)")
    println(io, "\t efConstruction = $(hnsw.efConstruction), ef = $(hnsw.ef)")

end


text:
using Revise
using Dice
using Dice: num_flips, num_nodes, ifelse

# Number of nodes SBK needs to model a distribution on n bits
sbk_num_nodes(n) = 2^n * (n - 1) + 3

function generate_code_sbk(p::Vector{Float64})
    @dice begin
        function helper(i)
            if i == length(p)
                DistInt(i - 1)
            else
                ifelse(flip(p[i] / sum(p[i:length(p)])),
                    DistInt(dicecontext(), i - 1),
                    helper(i + 1))
            end
        end
        helper(1)
    end
end

function test_sbk_num_nodes(p::Vector{Float64})
    code = generate_code_sbk(p)
    bdd = compile(code)
    @assert infer(code, :bdd) ≈ p  # Verify correctness

    bdd_nodes = num_nodes(bdd)
    num_bits = round(Int, log2(length(p)))
    @assert bdd_nodes == sbk_num_nodes(num_bits)  # Verify num bits

    println("SBK for a $(num_bits)-bit distribution: $(bdd_nodes) bdd nodes, $(num_flips(bdd)) flips")
end

dist_2 = [0.5, 0.5]
dist_4 = [0.1, 0.2, 0.3, 0.4]
dist_8 = [0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14, 0.44]
dist_16 = [0.005, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.001, 0.11, 0.12, 0.13, 0.14, 0.022, 0.032]

for dist in [dist_2, dist_4, dist_8, dist_16]
    test_sbk_num_nodes(dist)
end

# To investigate: the below discrete() fuction, which only meaningfully changes
# the order of flip introduction, takes fewer nodes (but still O(n2^n)).

#==                                       vvvvvvvvvvvvvvv diff
sbk_num_nodes(n) = 2^n * (n - 1) + 3      - 2^(n - 1) + 1

function generate_code_sbk(p::Vector{Float64})
    @dice begin
        function discrete(p::Vector{Float64})
            mb = length(p)
            v = Vector(undef, mb)
            sum = 1
            for i=1:mb
                v[i] = p[i]/sum
                sum = sum - p[i]
            end

            ans = DistInt(dicecontext(), mb-1)
            for i=mb-1:-1:1
                ans = if flip(v[i]) DistInt(dicecontext(), i-1) else ans end
            end
            return ans
        end
        discrete(p)
    end
end
==#


text:
## Deprecations for Languages

function tokenize(::Type{S}, s::T) where {S <: Language, T <: AbstractString}
    depwarn("Use of Languages as types is deprecated. Use instances.",  Symbol(S))
    tokenize(S(), s)
end

function ngramize(::Type{S}, words::Vector{T}, n::Int) where {S <: Language, T <: AbstractString}
    depwarn("Use of Languages as types is deprecated. Use instances.",  Symbol(S))
    ngramize(S(), words, n)
end

function onegramize(::Type{S}, words::Vector{T}) where {S <: Language, T <: AbstractString}
    depwarn("Use of Languages as types is deprecated. Use instances.",  Symbol(S))
    onegramize(S(), words)
end

function stem_all(stemmer::Stemmer, lang::Type{S}, sentence::AbstractString) where S <: Language
    depwarn("Use of Languages as types is deprecated. Use instances.",  Symbol(S))
    stem_all(stemmer, S(), sentence)
end

#pre-processing functions

function remove_whitespace!(entity::(Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_whitespace! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_whitespace)
end
function remove_punctuation!(entity::(Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_punctuation! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_punctuation)
end
function remove_nonletters!(entity::(Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_nonletters! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_non_letters)
end
function remove_numbers!(entity::(Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_numbers! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_numbers)
end
function remove_articles!(entity::(Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_articles! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_articles)
end
function remove_indefinite_articles!(entity::(Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_indefinite_articles! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_indefinite_articles)
end
function remove_definite_articles!(entity::(Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_definite_articles! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_definite_articles)
end
function remove_prepositions!(entity::(Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_prepositions! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_prepositions)
end
function remove_pronouns!(entity::(Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_pronouns! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_pronouns)
end
function remove_stop_words!(entity::(Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_stop_words! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_stopwords)
end


text:
@testset "High-level operations on distributed data" begin

    W = addprocs(2)
    @everywhere using GigaSOM

    Random.seed!(1)
    dd = rand(11111, 5)
    buckets = rand([1, 2, 3], 11111)

    di1 = scatter_array(:test1, dd, W[1:1])
    di2 = scatter_array(:test2, dd, W)
    buckets1 = scatter_array(:buckets1, buckets, W[1:1])
    buckets2 = scatter_array(:buckets2, buckets, W)

    dc = gather_array(di1)
    dc[:, 1:2] = asinh.(dc[:, 1:2] ./ 1.23)
    dtransform_asinh(di1, [1, 2], 1.23)
    dtransform_asinh(di2, [1, 2], 1.23)

    @testset "dtransform_asinh" begin
        @test isapprox(gather_array(di1), dc)
        @test isapprox(gather_array(di2), dc)
    end

    unscatter(di1)
    unscatter(di2)
    rmprocs(W)
end


text:
function retrieve!(quad::Quadtree, shape::Shape)
    global points_shape = Point[]
    level::Int64 = 0
    list_of_objects = Set{Shape}()

    for i = 1:length(shape.bounds[:,1])
        push!(points_shape, Point(shape.bounds[i,:]))
    end
    push!(points_shape, Point(shape.bounds[1,:]))

    function retrieve_inside!(quad::Quadtree, shape::Shape)
        Area = Float64
        s = Float64
        t = Float64

        ## Get all the objects in that level
        for i = 1:length(quad.objects)
            push!(list_of_objects, quad.objects[i])
        end

        ## Center of quadtree
        xx = (quad.x0 + quad.x1) / 2
        yy = (quad.y0 + quad.y1) / 2
        ## Test whether the half-quad is inside the shape
        for i = 0:length(shape.bounds[:,1])-3
            Area = 1/2*(-points_shape[2+i].y*points_shape[3+i].x + points_shape[1].y*(-points_shape[2+i].x + points_shape[3+i].x) + points_shape[1].x*(points_shape[2+i].y - points_shape[3+i].y) + points_shape[2+i].x*points_shape[3+i].y)

            s = 1/(2*Area)*(points_shape[1].y*points_shape[3+i].x - points_shape[1].x*points_shape[3+i].y + (points_shape[3+i].y - points_shape[1].y)*xx + (points_shape[1].x - points_shape[3+i].x)*yy)

            t = 1/(2*Area)*(points_shape[1].x*points_shape[2+i].y - points_shape[1].y*points_shape[2+i].x + (points_shape[1].y - points_shape[2+i].y)*xx + (points_shape[2+i].x - points_shape[1].x)*yy);

            # Test whether half-quad is inside the shape_test
            # if [ s>=0 && t>=0 && 1-s-t>=0 ] then it is inside
            if s>=0 && t>=0 && 1-s-t>=0
                if level <= max_level_quadtree[1]
                    if isdefined(quad, :ne)
                        level += 1
                        retrieve_inside!(quad.ne, shape)
                        level -= 1
                    end
                    if isdefined(quad, :nw)
                        level += 1
                        retrieve_inside!(quad.nw, shape)
                        level -= 1
                    end
                    if isdefined(quad, :sw)
                        level += 1
                        retrieve_inside!(quad.sw, shape)
                        level -= 1
                    end
                    if isdefined(quad, :se)
                        level += 1
                        retrieve_inside!(quad.se, shape)
                        level -= 1
                    end
                end
            end
        end

        ## Test whether the shape has points inside the quad
        if isdefined(quad, :ne)
            if isinside(quad.ne, points_shape)
                level += 1
                retrieve_inside!(quad.ne, shape)
                level -= 1
            end
        end
        if isdefined(quad, :nw)
            if isinside(quad.nw, points_shape)
                level += 1
                retrieve_inside!(quad.nw, shape)
                level -= 1
            end
        end
        if isdefined(quad, :sw)
            if isinside(quad.sw, points_shape)
                level += 1
                retrieve_inside!(quad.sw, shape)
                level -= 1
            end
        end
        if isdefined(quad, :se)
            if isinside(quad.se, points_shape)
                level += 1
                retrieve_inside!(quad.se, shape)
                level -= 1
            end
        end

        return

    end

    retrieve_inside!(quad, shape)
    return list_of_objects
end

function isinside(quad::Quadtree, points::Vector{Point})
    ## Test whether the shape has points inside the quad
    for i=1:length(points)
        xxx = points[i].x > quad.x0 && points[i].x < quad.x1
        yyy = points[i].y > quad.y0 && points[i].y < quad.y1

        if xxx && yyy
            return true
        end
    end

    ## Test line intersections
    l = quad.x0
    r = quad.x1
    t = quad.y0
    b = quad.y1
    for i=1:length(points)-1
        x1 = points_shape[i].x
        y1 = points_shape[i].y
        x2 = points_shape[i+1].x
        y2 = points_shape[i+1].y

        m = (x2-x1) != 0 ? (y2-y1)/(x2-x1) : pi
        c = y1 - (m*x1)

        if m>0 && m!=pi
            top_intersection = (m*l  + c)
            bottom_intersection = (m*r  + c)
        else
            top_intersection = (m*r  + c)
            bottom_intersection = (m*l  + c)
        end

        if y1 < y2
            toptrianglepoint = y1
            bottomtrianglepoint = y2
        else
            toptrianglepoint = y2
            bottomtrianglepoint = y1
        end

        topoverlap = top_intersection > toptrianglepoint ? top_intersection : toptrianglepoint
        botoverlap = bottom_intersection < bottomtrianglepoint ? bottom_intersection : bottomtrianglepoint

        isinside2 = (topoverlap < botoverlap) && (!((botoverlap<t)||(topoverlap>b)))
        if isinside2
            return true
        end
    end

    return false
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

#=
Notes for CSG bounding boxes
    The intersection nodes could compute bounding boxes as the CSGTree is being constructed because the bbox can only be equal or smaller than either of the children. Subtrees of union nodes should be processed as a whole so they can be properly subdivided or the bounding boxes are likely to be very inefficient.
    Need to cluster nodes geometrically to subdivide well but the order in which the unions are called by the user's code is unlikely to have the best geometric distribution.
    This makes the recursion a bit more complicated than conventional BVH because of the interleaved union and intersection operations.
=#

"""if fewer than this many primitives in a bounding box then do no subdivide further"""
const minprims = 4

abstract type BVHNode{T<:Real,S<:Primitive{T}} end

struct SplitBVH{T<:Real,S<:Primitive{T}} <: BVHNode{T,S}
    left::BVHNode{T,S}
    right::BVHNode{T,S}
    axis::Int64
end

struct LeafBVH{T<:Real,S<:Primitive{T}} <: BVHNode{T,S}
    primitives::Vector{S}
end

struct BVHData{T<:Real,S<:Primitive{T}}
    primitive::S
    centroid::SVector{3,T}
    boundingbox::BoundingBox{T}
end

centroid(a::BVHData) = a.centroid
boundingbox(a::BVHData) = a.boundingbox
primitive(a::BVHData) = a.primitive

struct PrimitiveData{T<:Real,S<:Primitive{T}}
    data::Vector{BVHData{T,S}}
    boundingbox::BoundingBox{T}

    function PrimitiveData(primitives::Vector{S}) where {T<:Real,S<:Primitive{T}}
        combined = Vector{BVHData{T,S}}(undef, length(primitives))
        totalbox = nothing

        for i in eachindex(primitives)
            box = BoundingBox(primitives[i])
            combined[i] = BVHData{T,S}(primitives[i], centroid(primitives[i]), box)
            totalbox = union(totalbox, box)
        end

        return new{T,S}(combined, totalbox)
    end
end

function computebvh(a::PrimitiveData)
    subdivide(a.boundingbox, a)

end

function subdivide(parentbounds::BoundingBox{T}, primitives::AbstractVector{S}) where {T<:Real,S<:Primitive{T}}
    if length(primitives <= minprims)
        return LeafBVH(primitives)
    else
        # test subdivision in 3 axes.
        #partition primitives
        # compute costs using parentbounds
        # return SplitBVH(subdivide(reduce(union,leftprims),leftprims),subdivide(reduce(union,rightprims,rightprims))
    end
end

cost(enclosingbox::BoundingBox{T}, childboxes::AbstractVector{BoundingBox{T}}) where {T<:Real} = sum(area.(childboxes)) / area(enclosingbox)

cost(boxes::Vector{BVHData{T,S}}) where {T,S} = sum(area.(boundingbox.(boxes)))

function partition!(a::Vector{S}, valfunc::Function, value::T) where {S,T<:Real}
    lower = 0
    upper = lastindex(a) + 1

    while true
        while upper - 1 >= 1 && valfunc(a[upper - 1]) > value
            upper = upper - 1
        end
        while lower + 1 <= lastindex(a) && valfunc(a[lower + 1]) <= value
            lower = lower + 1
        end

        if upper - lower == 1
            if lower == 0
                return (nothing, a)
            else
                if upper == lastindex(a) + 1
                    return (a, nothing)
                else
                    return (view(a, 1:lower), view(a, upper:lastindex(a)))
                end
            end
        end

        temp = a[upper - 1]
        a[upper - 1] = a[lower + 1]
        a[lower + 1] = temp
        upper = upper - 1
        lower = lower + 1
    end
end

"""compute compacted version of bounding volume hierarchy that can be traversed more quickly than the original tree structure, without recursion"""
function linearizebvh(bvhtree::BVHNode{T,S}) where {T<:Real,S<:Primitive{T}} end

"""traverses BVH and computes surface intersection"""
function surfaceintersection(a::BVHNode{T,S}, r::Ray{T,N}) where {T<:Real,S<:Primitive{T},N} end

# function testpartition()
#     split = .5

#     for i in 1:10000000
#         a = rand(5)
#         lower,upper = partition!(a,(x)->x,split)

#         if lower !== nothing
#             for i in lower
#                 @assert i < split
#             end
#         end

#         if upper !== nothing
#             for i in upper
#                 @assert i > split
#             end
#         end
#     end
# end

export PrimitiveData


text:
## ---------------------------------------------------------
# Ping
const _PING_SIGNAL_FILE_NAME = "gl-ping-signal"
_ping_signal_file(gl::GitLink) =
    joinpath(global_state_dir(gl), _PING_SIGNAL_FILE_NAME)

function _write_ping_signal(gl::GitLink, vtime = 60.0)
    fn = _ping_signal_file(gl)
    _mkdir(fn)
    extime = time() + vtime
    write(fn, string(extime))
    return extime
end

function _read_ping_signal(gl::GitLink)
    fn = _ping_signal_file(gl)
    !isfile(fn) && return -1.0
    txt = read(fn, String)
    extime = tryparse(Float64, txt)
    extime = isnothing(extime) ? -1.0 : extime
    return extime
end

_do_ping(gl::GitLink) = _read_ping_signal(gl) > time()


text:
@testset "Embeddings" begin
  Qx, x = QQ["x"]
  K, a = NumberField(x^3 - 2, "a")
  K2, (a2, ) = NumberField([x^3 - 2], "a")
  QQx, x = Hecke.rationals_as_number_field()[1]["x"]
  K3, a3 = NumberField(x^3 - 2, "a")
  K4, (a4, ) = NumberField([x^3 - 2], "a")
  # Let's go a bit crazy with the type of the base field
  for (K, a) in [(K, a), (K2, a2), (K3, a3), (K4, a4)]
    Kt, t = K["t"]
    g = t^3 + a*t + 5
    L, b = NumberField(g, "b")
    r, s = @inferred signature(L)
    @test (r, s) == (1, 4)

    emb = @inferred complex_embeddings(L)
    @test emb === complex_embeddings(L)
    @test (@inferred number_field(emb[1])) === L

    # isreal
    @test length(real_embeddings(L)) == 1
    @test @inferred isreal(emb[1])
    @test count(isreal, emb) == 1
    @test count(!isreal, emb) == 8
    @test count(isimaginary, emb) == 8

    for e in emb
      if isreal(e)
        @test (@inferred conj(e)) == e
      else
        @test conj(e) in emb && conj(e) != e
      end
    end

    # print
    sprint(show, emb[1])
    sprint(show, MIME"text/plain"(), emb[1])

    # restrict
    @test number_field(@inferred restrict(emb[1], K)) === K
    Set([restrict(e, K) for e in emb]) == Set(complex_embeddings(K))
    Set([restrict(e, L) for e in emb]) == Set(emb)
    @test_throws ArgumentError restrict(emb[1], number_field(x^2 - 2, "a", cached = false)[1])

    # compare with complex_embeddings of absolute field
    Labs, LabstoL = @inferred absolute_simple_field(L)
    embabs = complex_embeddings(Labs)
    ainLabs = LabstoL\(L(a))
    binLabs = LabstoL\b

    for e in emb
      imageofb = @inferred e(b)
      @test overlaps(conj(e)(b), conj(e(b)))
      imageofa = e(L(a))
      @test count(ee -> overlaps(ee(binLabs), imageofb) && overlaps(ee(ainLabs), imageofa), embabs) == 1
      @test Hecke.radiuslttwopower(e(b, 128), -128)
    end

    # restrict along morphism
    LtoLabs = inv(LabstoL)
    for e in embabs
      ee = @inferred restrict(e, LtoLabs)
      @test overlaps(ee(LabstoL(gen(Labs))), e(gen(Labs)))
      @test ee == LtoLabs * e
      @test_throws ArgumentError restrict(e, LabstoL)
    end

    for e in embabs
      ee = @inferred extend(e, LabstoL)
      @test length(ee) == 1
      @test LabstoL * ee[1] == e
    end

    for e in emb
      ee = @inferred extend(e, LtoLabs)
      @test length(ee) == 1
      @test LtoLabs * ee[1] == e
    end

    z = (a + b)^100
    for p in [32, 128, 256, 1024]
      for e in emb
        fn = evaluation_function(e, p)
        @test Hecke.radiuslttwopower(fn(z), -p)
        c = e(z, p)
        @test Hecke.radiuslttwopower(c, -p)
      end
    end
  end

  # More complicated restriction

  Qx, x = QQ["x"]
  K, a = number_field(x^2 - 2, "a")
  L, b = number_field(polynomial(K, [-3, 0, 0, 1], cached = false))
  M, c = number_field(polynomial(L, [-5, 0, 1], cached = false))
  emb = complex_embeddings(M)
  @test Set([restrict(e, L) for e in emb]) == Set(complex_embeddings(L))
  @test Set([restrict(e, K) for e in emb]) == Set(complex_embeddings(K))

end


text:
import ProximalOperators: ProximableFunction, prox!, IndPSD

# TODO zero cone
const conemap = Dict{Symbol, ProximableFunction}(
    :Free => IndFree(),
    :Zero => IndZero(),
    :NonNeg => IndNonnegative(),
    :NonPos => IndNonpositive(),
    :SOC => IndSOC(),
    :SOCRotated => IndRotatedSOC(),
    :SDP => IndPSD(scaling=true),
    :ExpPrimal => IndExpPrimal(),
    :ExpDual => IndExpDual()
)
const badcones = ProximableFunction[]

# type DualCone{T<:ProximableFunction} <: ProximableFunction
#     C::T
# end

# dual{T<:ProximableFunction}(C::T) = DualCone{T}(C)
# dual(C::DualCone) = C.C

# function prox!(cone::DualCone, x, y)
#     prox!(cone.C, -x, y)
#     for i = 1:length(x)
#         y[i] = x[i] + y[i]
#     end
# end

mutable struct ConeProduct{N,T} <: ProximableFunction
    ranges::NTuple{N, UnitRange{Int64}}
    cones::T
end
#
ConeProduct() = ConeProduct{0,Any}((),())

# function ConeProduct(ranges::NTuple{N, UnitRange{Int64}}, cones::T) where {N,T}
#     return ConeProduct{N,T}(ranges, cones)
# end

toRanges(rangesIn::NTuple{N, UnitRange{Int64}}) where N = rangesIn

function toRanges(rangesIn::NTuple{N,Array{T,1}}) where {N, T<:Integer}
    ranges = Array{UnitRange{Int64},1}(undef, N)
    for j in 1:N
        range = rangesIn[j]
        for (i,el) in enumerate(range[1]:range[end])
            if el != range[i]
                @error "Invalid range in input"
            end
        end
        ranges[j] = range[1]:range[end]
    end
    return tuple(ranges...)
end

#Should accept tro tuples with UnitRange{Int64} and ProximableFunction
#Cant enforce types or ambigous with default constructor?
function ConeProduct(rangesIn, cones)
    ranges = toRanges(rangesIn)
    N = length(cones)
    @assert typeof(ranges) == NTuple{N, UnitRange{Int64}}
    @assert length(ranges) == N
    @assert typeof(cones) <: Tuple
    #Verify that ranges covers the whole range
    prevRangeEnd = 0
    for i = 1:N
      @assert ranges[i][1] == prevRangeEnd + 1
      prevRangeEnd = ranges[i][end]
      @assert typeof(cones[i]) <: ProximableFunction
    end
    T = typeof(cones)
    #println("Dumping cones input")
    #dump(cones)
    ConeProduct(ranges, cones)
end

#Wrapper for dual prox to avoid double duals
function proxDual!(y::AbstractArray, C::ProximableFunction, x::AbstractArray)
    prox!(y, C, -x)
    @simd for i = 1:length(x)
        y[i] = x[i] + y[i]
    end
end
# proxDual!(y, C::DualCone, x) = prox!(y, C.C, x)

function prox!(y::AbstractArray, C::ConeProduct{N,T}, x::AbstractArray) where {N,T}
    #TODO Paralell implementation
    for i = 1:N
        prox!(view(y, C.ranges[i]), C.cones[i], view(x, C.ranges[i]))
    end
end

## Some better dual projections
const myIndFree = IndFree()
proxDual!(y::AbstractArray, C::IndZero, x::AbstractArray) = prox!(y, myIndFree, x)
const myIndZero = IndPoint()
proxDual!(y::AbstractArray, C::IndFree, x::AbstractArray)          = prox!(y, myIndZero, x)
proxDual!(y::AbstractArray, C::IndNonnegative, x::AbstractArray)   = prox!(y, C, x)
proxDual!(y::AbstractArray, C::IndNonpositive, x::AbstractArray)   = prox!(y, C, x)
# TODO figure out if self dual PSD
#proxDual!(y::AbstractArray, C::IndPSD, x::AbstractArray)           = prox!(y, C, x)

function proxDual!(y::AbstractArray, C::ConeProduct{N,T}, x::AbstractArray) where {N,T}
    #TODO Paralell implementation
    for i = 1:N
        proxDual!(view(y, C.ranges[i]), C.cones[i], view(x, C.ranges[i]))
    end
end

""" Indicator of K2×K1*×R+ × K2*×K1×R+ ∈ (R^n,R^m,R)^2"""
mutable struct DualConeProduct{T1<:ConeProduct,T2<:ConeProduct} <: ProximableFunction
    K1::T1
    K2::T2
    m::Int64
    n::Int64
end

DualConeProduct(K1::T1,K2::T2) where {T1,T2} = DualConeProduct{T1,T2}(K1, K2, K1.ranges[end][end], K2.ranges[end][end])
function prox!(y::AbstractArray, K::DualConeProduct, x::AbstractArray)
    m, n = K.m, K.n
    K1, K2 = K.K1, K.K2
    nu = n+m+1
    xx = view(x, 1:n)
    xy = view(x, (n+1):(n+m))
    xr = view(x, (nu+1):(nu+n))
    xs = view(x, (nu+n+1):(nu+n+m))

    yx = view(y, 1:n)
    yy = view(y, (n+1):(n+m))
    yr = view(y, (nu+1):(nu+n))
    ys = view(y, (nu+n+1):(nu+n+m))

    prox!(    yx, K2, xx)
    proxDual!(yy, K1, xy)
    y[nu] = max(x[nu], 0)
    proxDual!(yr, K2, xr)
    prox!(    ys, K1, xs)
    y[2nu] = max(x[2nu], 0)
end


text:
function check(i, j)
    id, im = div(i, 9), mod(i, 9)
    jd, jm = div(j, 9), mod(j, 9)

    jd == id && return true
    jm == im && return true

    div(id, 3) == div(jd, 3) &&
    div(jm, 3) == div(im, 3)
end

const lookup = zeros(Bool, 81, 81)

for i in 1:81
    for j in 1:81
        lookup[i,j] = check(i-1, j-1)
    end
end

function solve_sudoku(callback::Function, grid::Array{Int64})
    (function solve()
        for i in 1:81
            if grid[i] == 0
                t = Dict{Int64, Void}()

                for j in 1:81
                    if lookup[i,j]
                        t[grid[j]] = nothing
                    end
                end

                for k in 1:9
                    if !haskey(t, k)
                        grid[i] = k
                        solve()
                    end
                end

                grid[i] = 0
                return
            end
        end

        callback(grid)
    end)()
end

function display(grid)
    for i in 1:length(grid)
        print(grid[i], " ")
        i %  3 == 0 && print(" ")
        i %  9 == 0 && print("\n")
        i % 27 == 0 && print("\n")
    end
end

grid = Int64[5, 3, 0, 0, 2, 4, 7, 0, 0,
             0, 0, 2, 0, 0, 0, 8, 0, 0,
             1, 0, 0, 7, 0, 3, 9, 0, 2,
             0, 0, 8, 0, 7, 2, 0, 4, 9,
             0, 2, 0, 9, 8, 0, 0, 7, 0,
             7, 9, 0, 0, 0, 0, 0, 8, 0,
             0, 0, 0, 0, 3, 0, 5, 0, 6,
             9, 6, 0, 0, 1, 0, 3, 0, 0,
             0, 5, 0, 6, 9, 0, 0, 1, 0]

solve_sudoku(display, grid)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ---
# title: 96. Unique Binary Search Trees
# id: problem96
# author: zhwang
# date: 2022-02-23
# difficulty: Medium
# categories: Dynamic Programming, Tree
# link: <https://leetcode.com/problems/unique-binary-search-trees/description/>
# hidden: true
# ---
#
# Given _n_ , how many structurally unique **BST 's** (binary search trees) that
# store values 1 ...  _n_?
#
# **Example:**
#
#
#
#     Input: 3
#     Output: 5
#     Explanation: Given _n_ = 3, there are a total of 5 unique BST's:
#
#        1         3     3      2      1
#         \       /     /      / \      \
#          3     2     1      1   3      2
#         /     /       \                 \
#        2     1         2                 3
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 19`
#
#
## @lc code=start
using LeetCode

num_of_trees(n::Int) = binomial(2 * n, n) - binomial(2 * n, n + 1)
## @lc code=end


text:
using Plots, GeometryTypes

@recipe function f(r::HyperRectangle)
    points = decompose(Point{2,Float64}, r)
    rectpoints = points[[1,2,4,3],:]
    xs = [pt[1] for pt in rectpoints];
    ys = [pt[2] for pt in rectpoints];
    seriestype := :shape
    color = :orange
	#m = (:black, stroke(0))
	s = Shape(xs[:], ys[:])
end

@recipe function f(pt::Point)
    xs = [pt[1]]
    ys = [pt[2]]
    seriestype --> :scatter
    color := :orange
    markersize := 6
    xs, ys
end

@recipe function f(l::LineSegment)
    xs = [ l[1][1], l[2][1] ]
    ys = [ l[1][2], l[2][2] ]
    # seriestype = :line
    color = :red
    lw := 3
    xs, ys
end

@recipe function f(rectList::Vector{<:HyperRectangle})
    for r in rectList
	  @series begin
         r
	  end
   end
end

@recipe function f(ptList::Vector{<:Point})
    for p in ptList
	  @series begin
         p
	  end
   end
end

@recipe function f(lineList::Vector{<:LineSegment})
    for l in lineList
	  @series begin
         l
	  end
   end
end
sizePlot = (400,400)
plot(opacity=0.5, legend=false, size=sizePlot, yaxis=( (0,10), 0:1:10), xaxis=( (0,10), 0:1:10), foreground_color_grid= :cyan)

# defined in terms of SW x,y and then w,h
rect1 = HyperRectangle(Vec(1,1), Vec(1.,1))
rect2 = HyperRectangle(Vec(3,3), Vec(1., 2))
rect3 = HyperRectangle(Vec(1,1), Vec(4.,4))
rect4 = HyperRectangle(Vec(3,3), Vec(4.,4))
rect5 = HyperRectangle(Vec(2,2), Vec(1.,1))

vecR = Vector{HyperRectangle}()
push!(vecR,rect1, rect2, rect3, rect4, rect5)

vecPt = Vector{Point{2, Float64}}()

##########TEST POINTS################
nl = LineSegment(Point(6,8), Point(-1,0))
np = Point(6,8) # This doesnt work with plot (error,ERROR: LoadError: No user recipe defined for Float32)
mp = Point(8,8)
zp = Point(1,1)
xp = Point(3,3)
#mp = Point(9,9)

push!(vecPt, mp, zp, np, xp)

n = rand(1:100)
footitle = "update check: $n"
plot(opacity=0.5, legend=false, size=sizePlot, yaxis=( (0,10), 0:1:10), xaxis=( (0,10), 0:1:10), foreground_color_grid= :lightgray, title=footitle)

# opacity

a = minmax_euclidean(rect1, rect2)
b = minmax_euclidean(rect2, rect3)
c = overlaps(rect1, rect2)
d = overlaps(rect2, rect3)

print("test $a \n $b \n $c \n $d \n ")

e = contains(rect3, rect1) #rect 1 is tentirely inside rect 3
i = overlaps(rect3, rect1) #rect 1 is tentirely inside rect 3
h = in(rect3, rect1) #rect 1 is tentirely inside rect 3
f = in(rect3, rect5) #rect 3 and 5 do not overlap at all
g = overlaps(rect3, rect4) #rec 3 and 4 definitely overlap
print("test $e, $i, $h, $f, $g")

#ashape = Shape([0,2,3,0],[0,0,3,4])
gui()

# Check collision of node and rectangle
# works for points (but overlaps not defined)
# p = Point(1,1)
gp = in(rect3, zp) #rec 3 and 4 definitely overlap

###############TEST LINES#################
# Check collision of edge and rectangle
aline = LineSegment(np, mp) #horz at top, line in middle of nowhere
bline = LineSegment(zp,xp) #short diagonal
cline = LineSegment(zp,np) #long diagonal

vecLine = Vector{LineSegment}()
push!(vecLine, aline, bline, cline)

stmt = (rect1, bline)
a = minmax_euclidean(rect1, rect2)
b = minmax_euclidean(rect2, rect3)
c = overlaps(rect1, rect2)
d = overlaps(rect2, rect3)

## Bline: collides with 2,4,5, but not 1 and 3?

#print("\nAre the rectangle and line colliding? $stmt")

# aline shouldn't collide with anything
# bline should collide with rect4, rect5,
# cline should collide with rect4, but not rect1

function decompRect(r::HyperRectangle) #GeometryTypes.HyperRectangle{2,Float64}
    corners = decompose(Point{2, Float64}, r)
    corners = [Point(pt) for pt in corners]
    lineBottom  = LineSegment(corners[1], corners[2])
    lineTop     = LineSegment(corners[3], corners[4])
    lineLeft    = LineSegment(corners[1], corners[3])
    lineRight   = LineSegment(corners[2], corners[4])
    lines = Vector{LineSegment}()
    push!(lines, lineTop, lineRight, lineBottom, lineLeft)
    return lines
end

oppa = 0.5
plot!(vecR, opacity=oppa)
plot!(mp, seriestype=:scatter, markersize=3, color=:blue)
plot!(vecPt)
plot!(vecLine)

lines = decompRect(rect1);

# test an edge now.

# let's assume nodes not in rectangle for now -- or double-check again, just in case?
function isCollidingEdge(edge::LineSegment, obsList::Vector{HyperRectangle})
    for obs in obsList
        rectLines = decompRect(obs)
        for line in rectLines
            if intersects(edge, line)[1]
                return true
            end
        end
    end
    return false
end

#testE = LineSegment(Point(1.,0), Point(1.,10)) #collinear seems to collide correctly!
testE = LineSegment(Point(0.,0), Point(1.,2))
iscol = isCollidingEdge(testE, vecR)
print("\nAre the rectangles and line colliding? $iscol\n")

plot!(lines)

testl = LineSegment(Point(0.,0), Point(5.,5)) #remember to use period for float64 instead of int64
testl2 = LineSegment(Point(1.,1), Point(3.,3))
testl3 = LineSegment(Point(1.,3), Point(3.,1))
plot!(testE, color=:red)
plot!(testl)
plot!(testl2)
plot!(testl3)

#TODO: test for collinear lines....

# treeNode
# id
# prevID #either that, or store edges?
# state
#
# graphNode
# id
# state
#
#
# Edge
# linsegement
# startid
# endied
#
# function isCollidingNode(n::GraphNode, obsList::Vector{HyperRectangle})
    # for rect in obsList
        # if in(rect, n.state)
            # return true
        # end
    # end
    # return false
# end
#
# function isCollidingEdge(e::Edge, obsList::Vector{HyperRectangle})
    # for rect in obsList
        # if overlaps
#
#plot!(s)

#dist = min_euclidean(rect1, rect2)
#print("Dist $dist")
#x = 1:10; y = rand(10,2) # 2 columns means two lines
#plot(x,y)


text:
let
    biomass_f = 0.01

    # Feed jobs
    Ch = Channel(1) do ch
        cGLCs = Hd.val("cGLC")
        for (exp, cGLC)  in enumerate(cGLCs)
            put!(ch, (exp, cGLC))
        end
    end

    @threads for thid in 1:nthreads()
        for (exp, cGLC) in Ch
            # handle cache
            datfile = dat_file(string(DAT_FILE_PREFFIX, ME_Z_FIXXED_G_BOUNDED); exp)
            if isfile(datfile)
                lock(WLOCK) do
                    INDEX[ME_Z_FIXXED_G_BOUNDED, :DFILE, exp] = datfile
                    @info("Cached loaded (skipping)",
                        exp, cGLC, datfile, thid
                    ); println()
                end
                continue
            end

            # setup
            model = load_model(exp)
            objidx = ChU.rxnindex(model, iJR.BIOMASS_IDER)
            M, N = size(model)
            exp_growth = Hd.val("D", exp)
            ChU.ub!(model, iJR.BIOMASS_IDER, exp_growth * (1.0 + biomass_f))
            ChU.lb!(model, iJR.BIOMASS_IDER, exp_growth * (1.0 - biomass_f))
            model = ChLP.fva_preprocess(model,
                batchlen = 50,
                check_obj = iJR.BIOMASS_IDER,
                verbose = true
            )

            lock(WLOCK) do
                @info("Doing $(ME_Z_FIXXED_G_BOUNDED)",
                    exp, cGLC, thid
                ); println()
            end

            # maxent
            epout = ChEP.maxent_ep(model; alpha = Inf, damp = 0.985, epsconv = 1e-4,
                        verbose = false, maxiter = 5000
                    )

            # storing
            lock(WLOCK) do
                # Storing
                dat = Dict()
                dat[:exp_beta] = 0.0
                dat[:epouts] = Dict(0.0 => epout)
                dat[:model] = model |> ChU.compressed_model

                # caching
                serialize(datfile, dat)
                INDEX[ME_Z_FIXXED_G_BOUNDED, :DFILE, exp] = datfile

                @info("Finished ", exp, thid)
                println()
            end
        end # for (exp, cGLC) in Ch
    end # for thid in 1:nthreads()
end


text:
using DataStructures

function solve(collisions=false)
    d = hcat(map(x -> parse.(Int, x), matchall.(r"(-?\d+)", readlines("input.txt")))...)'
    distance = zeros(Int, (size(d, 1), 1))
    for n in 1:1000
        d[:, 4:6] += d[:, 7:9]
        d[:, 1:3] += d[:, 4:6]
        distance += sum(abs.(d[:, 1:3]), 2)
        if collisions
            positions = DefaultDict(() -> Int[])
            foreach(x -> push!(positions[d[x, 1:3]], x), 1:size(d, 1))
            for (k, v) in positions
                if k != [0, 0, 0] && length(v) > 1
                    foreach(x -> d[x,:] = zeros(Int, 9), v)
                end
            end
        end
    end
    return distance, d
end

println(findmin(solve()[1])[2] - 1)
println(sum(solve(true)[2][:, 1] .!= 0))


text:
# This file contains the plotting function

# load packages
using PyPlot
using KernelDensity
using Distributions

include(pwd()*"/utilities/posteriorinference.jl")

# plotting function

doc"""
    analyse_results(problem, res::Result)

Prints the results of the mcmc alg. and also plots the parameter posterior
and diagnostics plots.
"""
function analyse_results(Theta, loglik, accept_vec, prior_vec,
  theta_true, burn_in, Theta_parameters, lower_q_int_limit::Real=2.5,upper_q_int_limit::Real=97.5)

  # calc acceptance rate
  acceptance_rate = sum(accept_vec[burn_in:end])/length(accept_vec[burn_in:end])
  nbr_out_side_prior = sum(prior_vec)

  # calc acf
  acf_c1 = autocor(Theta[1,burn_in:end],1:50)
  acf_c2 = autocor(Theta[2,burn_in:end],1:50)
  acf_c3 = autocor(Theta[3,burn_in:end],1:50)

  # print results
  @printf "True value for r : %.4f, estimated value %.4f \n" exp(theta_true[1]) exp(mean(Theta[1,burn_in+1:end]))
  @printf "True value for phi: %.4f, estimated value %.4f \n" exp(theta_true[2]) exp(mean(Theta[2,burn_in+1:end]))
  @printf "True value for sigma: %.4f, estimated value %.4f \n" exp(theta_true[3]) exp(mean(Theta[3,burn_in+1:end]))

  @printf "Accept rate: %.4f %% \n" acceptance_rate*100
  @printf "Nbr outside of prior: %d  \n" nbr_out_side_prior

  @printf "Posterior mean:\n"
  Base.showarray(STDOUT,round(mean(Theta[:,burn_in+1:end],2),2),false)
  @printf "\n"

  @printf "Posterior standard deviation:\n"
  Base.showarray(STDOUT,std(Theta[:,burn_in+1:end],2),false)
  @printf "\n"

  @printf "Posterior quantile intervals (2.5th and 97.5th quantiles as default):\n"
  Base.showarray(STDOUT,round(calcquantileint(Theta[:,burn_in+1:end],lower_q_int_limit,upper_q_int_limit),2),false)
  @printf "\n"

  @printf "Loss for parameter estimations:\n"
  Base.showarray(STDOUT,loss(theta_true, Theta[:,burn_in+1:end]),false)
  @printf "\n"

  # text and lable size
  text_size = 15
  label_size = 10

  # plot chains
  PyPlot.figure(figsize=(10,20))
  ax = axes()
  ax1 = PyPlot.subplot(311)
  PyPlot.plot(Theta[1,:])
  PyPlot.plot(ones(size(Theta,2),1)*theta_true[1], "k")
  PyPlot.ylabel(L"$\log r$",fontsize=text_size)
  ax1[:axes][:get_xaxis]()[:set_ticks]([])
  ax2 = PyPlot.subplot(312)
  PyPlot.plot(Theta[2,:])
  PyPlot.plot(ones(size(Theta,2),1)*theta_true[2], "k")
  PyPlot.ylabel(L"$\log \phi$",fontsize=text_size)
  ax2[:axes][:get_xaxis]()[:set_ticks]([])
  ax3 = PyPlot.subplot(313)
  PyPlot.plot(Theta[3,:])
  PyPlot.plot(ones(size(Theta,2),1)*theta_true[3], "k")
  PyPlot.ylabel(L"$\log \sigma$",fontsize=text_size)
  PyPlot.xlabel("Iteration",fontsize=text_size)
  ax[:tick_params]("both",labelsize = label_size)

  # plot chains after burn in
  PyPlot.figure(figsize=(10,20))
  ax = axes()
  ax1 = PyPlot.subplot(311)
  PyPlot.plot(Theta[1,burn_in:end])
  PyPlot.plot(ones(size(Theta[:,burn_in:end],2),1)*theta_true[1], "k")
  PyPlot.ylabel(L"$\log r$",fontsize=text_size)
  ax1[:axes][:get_xaxis]()[:set_ticks]([])
  ax2 = PyPlot.subplot(312)
  PyPlot.plot(Theta[2,burn_in:end])
  PyPlot.plot(ones(size(Theta[:,burn_in:end],2),1)*theta_true[2], "k")
  PyPlot.ylabel(L"$\log \phi$",fontsize=text_size)
  ax2[:axes][:get_xaxis]()[:set_ticks]([])
  ax1 = PyPlot.subplot(313)
  PyPlot.plot(Theta[3,burn_in:end])
  PyPlot.plot(ones(size(Theta[:,burn_in:end],2),1)*theta_true[3], "k")
  PyPlot.ylabel(L"$\log \sigma$",fontsize=text_size)
  PyPlot.xlabel("Iteration",fontsize=text_size)
  ax[:tick_params]("both",labelsize = label_size)

    # plot acf
  PyPlot.figure()
  PyPlot.subplot(311)
  PyPlot.plot(acf_c1,"--r")
  PyPlot.ylabel(L"log $r$")
  PyPlot.subplot(312)
  PyPlot.plot(acf_c2,"--r")
  PyPlot.ylabel(L"log $\phi$")
  PyPlot.subplot(313)
  PyPlot.plot(acf_c3,"--r")
  PyPlot.ylabel(L"log $\sigma$")

  # plot posterior dist

  # calc grid for prior dist

x_c1 = Theta_parameters[1,1]-0.5:0.01:Theta_parameters[1,2]+0.5
x_c2 = Theta_parameters[2,1]-0.5:0.01:Theta_parameters[2,2]+0.5
x_c3 = Theta_parameters[3,1]-0.5:0.01:Theta_parameters[3,2]+0.5

# calc grid for kernel dens est of marginal posterior
#x_r_kerneldens = linspace(minimum(Theta[1,burn_in:end]),maximum(Theta[1,burn_in:end]),500)
#x_phi_kerneldens = linspace(minimum(Theta[2,burn_in:end]),maximum(Theta[2,burn_in:end]),500)
x_c1_kerneldens = x_c1
x_c2_kerneldens = x_c2
x_c3_kerneldens = x_c3

# calc prior dist
priordens_c1 = pdf.(Uniform(Theta_parameters[1,1], Theta_parameters[1,2]), x_c1)
priordens_c2 = pdf.(Uniform(Theta_parameters[2,1], Theta_parameters[2,2]), x_c2)
priordens_c3 = pdf.(Uniform(Theta_parameters[3,1], Theta_parameters[3,2]), x_c3)

h1 = kde(Theta[1,burn_in:end])
h2 = kde(Theta[2,burn_in:end])
h3 = kde(Theta[3,burn_in:end])

PyPlot.figure()
ax = axes()

subplot(311)
PyPlot.plot(h1.x,h1.density, "b")
PyPlot.plot(x_c1,priordens_c1, "g")
PyPlot.plot((theta_true[1], theta_true[1]), (0, maximum(h1.density)), "k")
PyPlot.ylabel(L"log $r$",fontsize=text_size)

subplot(312)
PyPlot.plot(h2.x,h2.density, "b")
PyPlot.plot(x_c2,priordens_c2, "g")
PyPlot.plot((theta_true[2], theta_true[2]), (0, maximum(h2.density)), "k")
PyPlot.ylabel(L"log $\phi$",fontsize=text_size)

subplot(313)
PyPlot.plot(h3.x,h3.density, "b")
PyPlot.plot(x_c3,priordens_c3, "g")
PyPlot.plot((theta_true[3], theta_true[3]), (0, maximum(h3.density)), "k")
PyPlot.ylabel(L"log $\sigma$",fontsize=text_size)
ax[:tick_params]("both",labelsize=label_size)

# print loglik
PyPlot.figure()
PyPlot.plot(loglik)
PyPlot.ylabel("Loglik")

posterior_samples = Theta[:,burn_in:end]

PyPlot.figure()
counter = 1
for i = 1:3
  for j = 1:3
    PyPlot.subplot(3,3,counter) # set subplot
    if i == j
      # plot hist on diag. elements
      h = PyPlot.plt[:hist](posterior_samples[i,:],50)
      PyPlot.plot((theta_true[i], theta_true[i]), (0, maximum(h[1])), "k")
    else
      # plot scatters on non-diagonal elements
      PyPlot.scatter(posterior_samples[j,:], posterior_samples[i,:])
    end
    counter += 1

    if i == 1 && j == 1
      PyPlot.ylabel(L"log $r$")
    elseif i == 3 && j == 1
      PyPlot.xlabel(L"log $r$")
    end

    if j == 1 && i == 2
      PyPlot.ylabel(L"log $\phi$")
    elseif i == 3 && j == 2
      PyPlot.xlabel(L"log $\phi$")
    end

    if j == 1 && i == 3
      PyPlot.ylabel(L"log $\sigma$")
    elseif i == 3 && j == 3
      PyPlot.xlabel(L"log $\sigma$")
    end

  end
end

end

max_stars_repo_path: Ricker model/Results/plotting.jl

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_path: Ricker model/Results/plotting.jl

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_path: Ricker model/Results/plotting.jl

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
"""
    tear_graph(sys) -> sys

Tear the bipartite graph in a system. End users are encouraged to call [`structural_simplify`](@ref)
instead, which calls this function internally.
"""
function tear_graph(sys)
    find_solvables!(sys)
    s = structure(sys)
    @unpack graph, solvable_graph, var_eq_matching, scc = s

    @set! sys.structure.partitions = map(scc) do c
        ieqs = filter(eq->isalgeq(s, eq), c)
        vars = Int[var for var in invview(var_eq_matching)[ieqs] if var !== unassigned]

        ict = IncrementalCycleTracker(DiCMOBiGraph{true}(graph); dir=:in)
        SystemPartition(tearEquations!(ict, solvable_graph.fadjlist, ieqs, vars)...)
    end
    return sys
end

"""
    algebraic_equations_scc(sys)

Find strongly connected components of algebraic equations in a system.
"""
function algebraic_equations_scc(sys)
    s = get_structure(sys)
    if !(s isa SystemStructure)
        sys = initialize_system_structure(sys)
        s = structure(sys)
    end

    # skip over differential equations
    algvars = isalgvar.(Ref(s), 1:ndsts(s.graph))
    var_eq_matching = complete(matching(s, algvars, s.algeqs))
    components = find_scc(s.graph, var_eq_matching)

    @set! sys.structure.var_eq_matching = var_eq_matching
    @set! sys.structure.scc = components
    return sys
end


text:
module NiLang

using Reexport
@reexport using NiLangCore
import NiLangCore: ⊕, ⊖

include("utils.jl")
include("vars.jl")
include("instructs.jl")
include("functions.jl")
include("nnlib.jl")

include("autodiff/autodiff.jl")

export AD
end # module


text:
# An acceleration index that helps do a `join` over the listed columns
abstract type Index{names}; end

@inline labels(i::Index{names}) = _labels(names...)

# Similar to above but the values are unique
abstract type UniqueIndex{names} <: Index{names}; end

# There is no acceleration but the assertion of uniqueness may make joins, etc, faster.
struct Unique{names} <: UniqueIndex{names}
end