text:
# exgraph.jl - expression graph as a list of primitive expression nodes

abstract type AbstractExGraph end

mutable struct ExGraph <: AbstractExGraph
    tape::Vector{ExNode}           # list of ExNode-s
    idx::Dict{Symbol, ExNode}      # map from var name to its node in the graph
    ctx::Dict{Any,Any}             # settings and caches
end

function ExGraph(; ctx=Dict(), inputs...)
    ctx = to_context(ctx)
    # @get_or_create(ctx, :mod, @__MODULE__)
    @get_or_create(ctx, :mod, get_caller_module())
    g = ExGraph(ExNode[], Dict(), ctx)
    for (var, val) in inputs
        push!(g, :input, var, var; val=val)
    end
    return g
end

function ExGraph(ex::Expr; fuse=true, ctx=Dict(), inputs...)
    ex = preprocess(ex)
    ctx = to_context(ctx)
    g = ExGraph(;ctx=ctx, inputs...)
    g.ctx[:expr] = ex
    method = @get(g.ctx, :method, :parse)
    if method == :parse
        parse!(g, ex)
    elseif method == :track
        eval_tracked!(g, ex, inputs...)
    else
        error("Method $method is not supported")
    end
    if fuse
        g = fuse_assigned(g)
    end
    return g
end

function Base.deepcopy(g::AbstractExGraph)
    ctx_copy = to_context(Dict())
    for (k, v) in g.ctx
        if isa(v, Module)
            ctx_copy[k] = v
        else
            ctx_copy[k] = deepcopy(v)
        end
    end
    return typeof(g)(deepcopy(g.tape), deepcopy(g.idx), ctx_copy)
end

function Base.show(io::IO, g::ExGraph)
    print(io, "ExGraph\n")
    for node in g.tape
        print(io, "  $node\n")
    end
end

Base.haskey(g::AbstractExGraph, var::Symbol) = haskey(g.idx, var)
Base.in(var::Symbol, g::AbstractExGraph) = haskey(g, var)
Base.lastindex(g::AbstractExGraph) = lastindex(g.tape)
Base.length(g::AbstractExGraph) = length(g.tape)
Base.get(g::AbstractExGraph, var::Symbol) = g.idx[var]
Base.getindex(g::AbstractExGraph, var::Symbol) = g.idx[var]
Base.getindex(g::AbstractExGraph, var::String) = g.idx[Symbol(var)]
Base.getindex(g::AbstractExGraph, i::Integer) = g.tape[i]
Base.setindex!(g::AbstractExGraph, nd::ExNode, i::Integer) =
    (g.tape[i] = nd; g.idx[varname(nd)] = nd)

# indexof(g::AbstractExGraph, vname::Symbol) = findfirst(map(varname, g.tape), vname)
indexof(g::AbstractExGraph, vname::Symbol) = findfirst(isequal(vname), map(varname, g.tape))
getsize(g::AbstractExGraph, vname::Symbol) = g.ctx[:rsizes][vname]
getsize(g::AbstractExGraph, nd::ExNode) = getsize(g, varname(nd))

Base.iterate(g::ExGraph) = length(g) >= 1 ? (g[1], 2) : nothing
Base.iterate(g::ExGraph, state) = length(g) >= state ? (g[state], state + 1) : nothing

function Base.cat(g1::T, g2::T) where T<:AbstractExGraph
    return T(vcat(g1.tape, g2.tape), merge(g1.idx, g2.idx), merge(g1.ctx, g2.ctx))
end

function to_expr(g::AbstractExGraph)
    res = quote end
    for nd in g.tape
        if !isa(nd, ExNode{:input})
            push!(res.args, to_expr(nd))
        end
    end
    return res
end

function to_expr_kw(g::AbstractExGraph)
    res = quote end
    for nd in g.tape
        if !isa(nd, ExNode{:input})
            push!(res.args, to_expr_kw(nd))
        end
    end
    return res
end

function eval_tracked!(g::ExGraph, ex::Expr, inputs...)
    og = swap_default_graph!(g)
    evex = ex.head == :block ? deepcopy(ex) : Expr(:block, deepcopy(ex))
    for (var, val) in inputs
        tv = isa(val, AbstractArray) ? TrackedArray(g, var, val) : TrackedReal(g, var, val)
        pushfirst!(evex.args, :($var = $tv))
    end
    Core.eval(@get(g.ctx, :mod, Main), evex)
    new_g = reparse(swap_default_graph!(og); all=true)
    # copy nodes, but keep old context; maybe we will introduce special copyto! for this later
    g.tape = new_g.tape
    g.idx = new_g.idx
    return g
end

"""Generate a new unique name for intermediate variable in graph"""
function genname()
    s = String(gensym())
    return Symbol(replace(s, "##" => "tmp"))
end

function genname(prefix::String)
    s = String(gensym())
    return Symbol(replace(s, "##" => prefix))
end

function genname(prefix::Symbol)
    s = String(gensym())
    return Symbol(replace(s, "##" => prefix))
end

function gennames(count::Int)
    return [genname() for _=1:count]
end

function rename_repeated(g::AbstractExGraph, ex)
    st = @get_or_create(g.ctx, :renamings, Dict())
    st = prop_subs(st)
    return subs(ex, st)
end
rename_repeated(g::AbstractExGraph, ex::ExH) = rename_repeated(g, Expr(ex))

function add_renaming!(g::AbstractExGraph, var::Symbol)
    old = var
    while var in g
        var = inc_var_name(var)
    end
    st = @get_or_create(g.ctx, :renamings, Dict())
    st[old] = var
    return var
end

## push!, insert!, delete!

"""
Add a new node to a graph. Expression should be simple, e.g.
nested calls or blocks are not allowed (use parse!() for it).
"""
function Base.push!(g::AbstractExGraph, nd::ExNode)
    @assert(!haskey(g, varname(nd)),
            "Graph already contains a node with name $(varname(nd))!")
    push!(g.tape, nd)
    g.idx[varname(nd)] = nd
    return varname(nd)
end

function Base.push!(g::AbstractExGraph, C::Symbol, var::Union{Symbol,Expr}, ex::Any;
                    val=nothing, meta=Dict())
    @assert(!haskey(g, split_indexed(var)[1]),
            "Graph already contains a node with name $(var)!")
    nd = ExNode{C}(var, ex; val=val, meta=meta)
    push!(g, nd)
    return var
end

function Base.insert!(g::AbstractExGraph, i::Integer, nd::ExNode)
    g.idx[varname(nd)] = nd
    insert!(g.tape, i, nd)
end

function Base.insert!(g::AbstractExGraph, i::Integer, nds::Vector{ExNode})
    for (j, nd) in enumerate(nds)
        insert!(g, i + j - 1, nd)
    end
end

function Base.delete!(g::AbstractExGraph, i::Integer)
    delete!(g.idx, varname(g[i]))
    deleteat!(g.tape, i)
end

function Base.delete!(g::AbstractExGraph, vname::Symbol)
    delete!(g.idx, vname)
    i = findall(nd -> varname(nd) == vname, g.tape)
    deleteat!(g.tape, i)
end

## parse!

"""
Parse Julia expression and build ExGraph in-place.
Return the the output variable.
"""
parse!(g::AbstractExGraph, ex::Expr) = parse!(g, ExH(ex))
parse!(g::AbstractExGraph, ::LineNumberNode) = :nil
parse!(g::AbstractExGraph, ::ExH{:line}) = :nil
parse!(g::AbstractExGraph, s::Symbol) = rename_repeated(g, s)

function parse!(g::AbstractExGraph, x::Number)
    if haskey(g.ctx, :bitness)
        x = force_bitness(x, Val(g.ctx[:bitness]))
    end
    var = push!(g, :constant, genname(), x; val=x)
    return var
end

function parse!(g::AbstractExGraph, x::AbstractArray)
    if haskey(g.ctx, :bitness)
        x = force_bitness(x, Val(g.ctx[:bitness]))
    end
    var = push!(g, :constant, genname(), x; val=x)
    return var
end

function parse!(g::AbstractExGraph, x::ExH{:vect})
    @assert all(e -> e isa Number, x.args) "Can only create vector literal node with numbers"
    nums = convert(Vector{Float64}, x.args)
    if haskey(g.ctx, :bitness)
        nums = force_bitness(nums, Val(g.ctx[:bitness]))
    end
    var = push!(g, :constant, genname(), nums; val=nums)
    return var
end

function parse!(g::ExGraph, ex::ExH{:(=)})
    lhs, rhs = ex.args
    rhs = rename_repeated(g, rhs)
    dep = parse!(g, rhs)
    if isa(lhs, Symbol)
        if haskey(g, lhs)
            lhs = add_renaming!(g, lhs)
        end
        push!(g, :(=), lhs, dep)
        return lhs
    elseif isa(lhs, Expr) && lhs.head == :tuple
        last_vname = nothing
        for (i, vname) in enumerate(lhs.args)
            parse!(g, :($vname = $dep[$i]))
        end
        return last_vname
    else
        error("LHS of $(Expr(ex)) is neither variable, nor tuple")
    end
end

function parse!(g::ExGraph, ex::ExH{:ref})
    ex = rename_repeated(g, ex)
    # @assert isa(ex.args[2], Number) "Currently only constant indices are supported"
    idx_vars = [parse!(g, v) for v in ex.args[2:end]]
    base = isa(ex.args[1], Symbol) ? ex.args[1] : parse!(g, ex.args[1])
    vname = push!(g, :ref, genname(), Expr(:ref, base, idx_vars...))
    return vname
end

# operations that are guaranted to never change their value
const CONST_OPS = Set([:length])

function parse!(g::ExGraph, ex::ExH{:call})
    ex = rename_repeated(g, ex)
    op = canonical(g.ctx[:mod], ex.args[1])
    args, kw_args = parse_call_args(ex)
    meta = isempty(kw_args) ? Dict() : Dict(:kw => kw_args)
    deps = [parse!(g, arg) for arg in args]
    pex = Expr(:call, op, deps...)
    if op in CONST_OPS && isempty(meta)
        var = push!(g, :constant, genname(), pex)
    else
        var = push!(g, :call, genname(), pex; meta=meta)
    end
    return var
end

function parse!(g::ExGraph, ex::ExH{:.})
    ex = rename_repeated(g, ex)
    if isa(ex.args[2], Expr) && ex.args[2].head == :tuple
        # broadcasting
        op = canonical(g.ctx[:mod], ex.args[1])
        deps = [parse!(g, arg) for arg in ex.args[2].args]
        # pex = Expr(:call, op, deps...)
        pex = Expr(:., op, Expr(:tuple, deps...))
        var = push!(g, :bcast, genname(), pex)
        return var
    elseif isa(ex.args[2], QuoteNode)
        # field
        var = push!(g, :field, genname(), ex)
        return var
    end
end

function parse!(g::ExGraph, ex::ExH{Symbol("'")})
    ex = rename_repeated(g, ex)
    dep = parse!(g, ex.args[1])
    pex = :(transpose($dep))
    var = push!(g, :call, genname(), pex)
    return var
end

function parse!(g::AbstractExGraph, ex::Union{ExH{:block}, ExH{:body}})
    deps = [parse!(g, arg) for arg in ex.args]
    return deps[end]
end

function parse!(g::ExGraph, ex::ExH{:tuple})
    ex = rename_repeated(g, ex)
    deps = [parse!(g, arg) for arg in ex.args]
    pex = Expr(:tuple, deps...)
    vname = push!(g, :tuple, genname(), pex)
end

## reparsing of opaque nodes

function reparse(g::AbstractExGraph; all=false)
    new_g = reset_tape(g)
    for nd in g.tape
        if isa(nd, ExNode{:input}) || isa(nd, ExNode{:constant})
            push!(new_g, nd)
        elseif all || isa(nd, ExNode{:opaque})
            parse!(new_g, to_expr(nd))
        else
            push!(new_g, nd)
        end
    end
    return fuse_assigned(new_g)
end

## graph simlification

istemp(var::Symbol) = startswith(string(var), "tmp")

function external_vars(g::AbstractExGraph)
    ext_vnames = Set{Symbol}()
    for nd in g.tape
        for dep in dependencies(nd)
            if !haskey(g, dep)
                push!(ext_vnames, dep)
            end
        end
    end
    return ext_vnames
end

function assign_chain!(g::AbstractExGraph, nd::ExNode{:(=)},
                       guards::Vector{Expr}, chain::Vector{Symbol})
    if getguards(nd) == guards #  && !is_special_expr(getexpr(nd))
        push!(chain, varname(nd))
        dep = dependencies(nd)[1]
        if haskey(g, dep) && !isa(g[dep], ExNode{:input}) && !isa(g[dep], ExNode{:constant})
            dep_nd = g[dep]
            assign_chain!(g, dep_nd, guards, chain)
        end
    end
    return chain
end

function assign_chain!(g::AbstractExGraph, nd::ExNode{C},
                       guards::Vector{Expr}, chain::Vector{Symbol}) where C
    if getguards(nd) == guards
        push!(chain, varname(nd))
    end
end

"""
Collect all replacable variables from a chain of assignments in a graph.
Variables `y` and `x` are considered replacable if there's a node `y = x`
and both variables have the same set of guards.
Note that this allows nodes to have different sets of indices.
"""
assign_chain(g::AbstractExGraph, nd::ExNode{C}) where {C} =
    assign_chain!(g, nd, getguards(nd), Vector{Symbol}())

function assign_chain_index_replacements(g::AbstractExGraph, chain::Vector{Symbol})
    nd = g[chain[1]]
    root_nd = g[chain[end]]
    st = Dict(zip(varidxs(nd), varidxs(nd)))
    for i=2:length(chain)
        prev_idxs = get_indices(getexpr(g[chain[i-1]]))[1]
        cur_idxs = varidxs(g[chain[i]])
        pair_st = Dict(zip(prev_idxs, cur_idxs))
        new_st = Dict()
        for (k, v) in st
            if haskey(pair_st, v)
                # propagate replacements
                new_st[k] = pair_st[v]
            end
        end
        st = new_st
    end
    # indices that occur in RHS of the root assignment node, but not on its LHS
    root_free_idxs_ = find_indices(to_expr(root_nd)) |> flatten |> Set
    root_free_idxs = Any[idx for idx in root_free_idxs_ if !in(idx, values(st))]
    free_st = index_replacements(Set(keys(st)), root_free_idxs)
    rev_st = Dict(zip(values(st), keys(st)))
    return merge(rev_st, free_st)
end

"""
Collapse unnecessary assignment nodes, rewriting all affected nodes. Example:

    tmp1 = x * y
    z = tmp1

will be rewritten to

    z = x * y
"""
function fuse_assigned(g::AbstractExGraph; outvars=nothing)
    new_g = reset_tape(g)
    for nd in g.tape
        if isa(nd, ExNode{:(=)})
            chain = assign_chain(g, nd)
            root_assign_nd = g[chain[end]]
            new_ex = getexpr(root_assign_nd)
            # st = assign_chain_index_replacements(g, chain)
            # new_ex = subs(new_ex_, st)
            new_nd = copy(root_assign_nd; var=getvar(nd), ex=new_ex)
            push!(new_g, new_nd)
        else
            push!(new_g, copy(nd))
        end
    end
    new_g = remove_unused(new_g, outvars == nothing ? [varname(new_g[end])] : outvars)
    return new_g
end

function fuse_assigned(g::ExGraph; outvars=nothing)
    new_g = reset_tape(g)
    for nd in g.tape
        if isa(nd, ExNode{:(=)})
            chain = assign_chain(g, nd)
            root_assign_nd = g[chain[end]]
            new_ex = getexpr(root_assign_nd)
            new_nd = copy(root_assign_nd; var=getvar(nd), ex=new_ex)
            push!(new_g, new_nd)
        else
            push!(new_g, copy(nd))
        end
    end
    new_g = remove_unused(new_g, outvars == nothing ? [varname(new_g[end])] : outvars)
    return new_g
end


text:
@testset "Lines" begin
    @testset "linesegment intersect" begin
        a = LineSegment(Point2f0(0,0), Point2f0(1,0))
        b = LineSegment(Point2f0(0.5,0.5), Point2f0(0.5,-0.5))
        result = intersects(a,b)
        @test result[1] == true
        @test result[2] ≈ Point2f0(0.5,0.0)

        a = LineSegment(Point2f0(0,0), Point2f0(0.499,0))
        b = LineSegment(Point2f0(0.5,0.5), Point2f0(0.5,-0.5))
        result = intersects(a,b)
        @test intersects(a,b) == (false, zero(Point2f0))

        l = LineSegment(Point(6.0,11), Point(12.0,11))
        r1 = LineSegment(Point(9.0,10.0), Point(9.0, 13.0))
        r2 = LineSegment(Point(9.34365, 10.890), Point(9.34365, 13.402))
        intersect, p = intersects(l, r2)
        @test intersect
        @test p ≈ Point(9.34365, 11.0)
        intersect, p = intersects(l, r1)
        @test intersect
        @test p ≈ Point(9.0, 11.0)
    end

    @testset "linesegment intersect: #133" begin
        s1 = LineSegment(Point{2,Float64}[[0.88033, 1.28211], [3.525, 1.28211]]...)
        s2 = LineSegment(Point{2,Float64}[[2.0375, 1.1625], [2.0375, 1.2875]]...)
        s1_old = deepcopy(s1)
        s2_old = deepcopy(s2)
        result = GeometryTypes.intersects(s1, s2)
        @test s1 == s1_old # make sure it doesn't mutate inputs
        @test s2 == s2_old # make sure it doesn't mutate inputs
        @test result[1] == true
        @test result[2] ≈ Point(2.0375, 1.28211)
    end

    @testset "poly self intersections" begin
        p1 = Point2f0[(x, sin(x)) for x in 0:.5:6]
        p2 = Point2f0[(p1[i][1], 0) for i=length(p1):-1:1]
        points = vcat(p1, p2)
        inds, intersects = self_intersections(points)
        @test inds == [7,19]
        @test intersects ≈ Point2f0[(3.1434426,0.0)]
        polys = split_intersections(points)
        @test length(polys) == 2
        @test length(polys[1]) == 12
        @test length(polys[2]) == 15

        u = range(0, stop=2pi, length=200)
        points = Point2f0[(sin(x), sin(2x)) for x in u]

        inds, intersects = self_intersections(points)
        @test inds == [1,100]
        @test intersects == Point2f0[(0.0,0.0)]
        polys = split_intersections(points)
        @test length(polys) == 2
        @test length(polys[1]) == 99
        @test length(polys[2]) == 102
    end
end

ext: jl


text:
export TeardownLogic, UnsubscribableTeardownLogic, CallableTeardownLogic, VoidTeardownLogic, InvalidTeardownLogic
export Teardown, as_teardown
export unsubscribe!, teardown!, on_unsubscribe!

export InvalidTeardownLogicTraitUsageError, MissingOnUnsubscribeImplementationError

import Base: show

"""
Abstract type for all possible teardown logic traits.

See also: [`UnsubscribableTeardownLogic`](@ref), [`CallableTeardownLogic`](@ref), [`VoidTeardownLogic`](@ref), [`InvalidTeardownLogic`](@ref)
"""
abstract type TeardownLogic end

"""
Unsubscribable teardown logic trait behavior. Unsubscribable teardown object must define its own method
for `on_unsubscribe!()` function which will be invoked when actor decides to `unsubscribe!` from Observable.

See also: [`TeardownLogic`](@ref), [`on_unsubscribe!`](@ref), [`unsubscribe!`](@ref)
"""
struct UnsubscribableTeardownLogic <: TeardownLogic end

"""
Callable teardown logic trait behavior. Callable teardown object must be callable (insert meme with a surprised Pikachu here).

See also: [`TeardownLogic`](@ref), [`on_unsubscribe!`](@ref), [`unsubscribe!`](@ref)
"""
struct CallableTeardownLogic       <: TeardownLogic end

"""
Void teardown logic trait behavior. Void teardown object does nothing in `unsubscribe!` and may not define any additional methods.

See also: [`TeardownLogic`](@ref), [`on_unsubscribe!`](@ref), [`unsubscribe!`](@ref)
"""
struct VoidTeardownLogic           <: TeardownLogic end

"""
Default teardown logic trait behavour. Invalid teardwon object cannot be used in `unsubscribe!` function. Doing so will raise an error.

See also: [`TeardownLogic`](@ref), [`on_unsubscribe!`](@ref), [`unsubscribe!`](@ref)
"""
struct InvalidTeardownLogic        <: TeardownLogic end

"""
Abstract type for any teardown object. Each teardown object must be a subtype of `Teardown`.

See also: [`TeardownLogic`](@ref)
"""
abstract type Teardown end

"""
    as_teardown(::Type)

This function checks teardown trait behavior specification. Should be used explicitly to specify teardown logic trait behavior for any object.

# Examples

```jldoctest
using Rocket

struct MySubscription <: Teardown end

Rocket.as_teardown(::Type{<:MySubscription}) = UnsubscribableTeardownLogic()
Rocket.on_unsubscribe!(s::MySubscription)    = println("Unsubscribed!")

subscription = MySubscription()
unsubscribe!(subscription)
;

# output

Unsubscribed!
```

See also: [`Teardown`](@ref), [`TeardownLogic`](@ref)
"""
as_teardown(::Type)             = InvalidTeardownLogic()
as_teardown(::Type{<:Function}) = CallableTeardownLogic()

"""
    unsubscribe!(o::T) where T

`unsubscribe!` function is used to cancel Observable execution and to dispose any kind of resources used during an Observable execution.

See also: [`Teardown`](@ref), [`TeardownLogic`](@ref), [`on_unsubscribe!`](@ref)
"""
unsubscribe!(teardown::T) where T = teardown!(as_teardown(T), teardown)

teardown!(::UnsubscribableTeardownLogic, teardown) = on_unsubscribe!(teardown)
teardown!(::CallableTeardownLogic,       teardown) = teardown()
teardown!(::VoidTeardownLogic,           teardown) = begin end
teardown!(::InvalidTeardownLogic,        teardown) = throw(InvalidTeardownLogicTraitUsageError(teardown))

"""
    on_unsubscribe!(teardown)

Each valid teardown object with UnsubscribableTeardownLogic trait behavior must implement its own method
for `on_unsubscribe!()` function which will be invoked when actor decides to `unsubscribe!` from Observable.

See also: [`Teardown`](@ref), [`TeardownLogic`](@ref), [`UnsubscribableTeardownLogic`](@ref)
"""
on_unsubscribe!(teardown) = throw(MissingOnUnsubscribeImplementationError(teardown))

# -------------------------------- #
# Errors                           #
# -------------------------------- #

"""
This error will be thrown if `unsubscribe!` function is called with invalid teardown object.

See also: [`unsubscribe!`](@ref)
"""
struct InvalidTeardownLogicTraitUsageError
    teardown
end

function Base.show(io::IO, err::InvalidTeardownLogicTraitUsageError)
    print(io, "Type $(typeof(err.teardown)) has undefined teardown behavior. \nConsider implement as_teardown(::Type{<:$(typeof(err.teardown))}).")
end

"""
This error will be thrown if Julia cannot find specific method of `on_unsubscribe!()` function for given teardown object.

See also: [`on_unsubscribe!`](@ref)
"""
struct MissingOnUnsubscribeImplementationError
    teardown
end

function Base.show(io::IO, err::MissingOnUnsubscribeImplementationError)
    print(io, "You probably forgot to implement on_unsubscribe!(unsubscribable::$(typeof(err.teardown))).")
end


text:
# Cut oracles for MathOptInterface cones
module Cuts

using LinearAlgebra
import LinearAlgebra.copytri!

import MathOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

include("arrayutilities.jl")
include("secondordercone.jl")
include("positivesemidefiniteconetriangle.jl")

end


text:
module TestZulipSnippetBot
using Test

for file in sort([file for file in readdir(@__DIR__) if
                                   occursin(r"^test[0-9]+_.*\.jl$", file)])
    m = match(r"test([0-9]+)_(.*).jl", file)
    filename = String(m[2])
    testnum = string(parse(Int, m[1]))

    # with this test one can run only specific tests, for example
    # Pkg.test("MyPackage", test_args = ["xxx"])
    # or
    # Pkg.test("MyPackage", test_args = ["1"])
    if isempty(ARGS) || (filename in ARGS) || (testnum in ARGS) || (m[1] in ARGS)
        @testset "$filename" begin
            # Here you can optionally exclude some test files
            # VERSION < v"1.1" && file == "test_xxx.jl" && continue

            include(file)
        end
    end
end

end # module

max_stars_repo_name: Arkoniak/ZulipSnippetBot

max_issues_repo_name: Arkoniak/ZulipSnippetBot

max_forks_repo_name: Arkoniak/ZulipSnippetBot


text:
import AxisArrays: axisdim, axisnames, axisname, axisvalues

# All ArrayTypes must define `axisnames` function
const ArrayTypes{T,M} = Union{AxisArray{T,M}} # JuMP.JuMPArray{T,M}

struct WrappedArray{T,N,M,D,Ax} <: AbstractArray{T,N}
    data::D
    axes::Ax
    # Duck-typing, we assume axisnames is defined for data
    function WrappedArray(data::T, axes::NTuple{N,Axis}) where {T<:Union{Number,Symbol}, N}
        new{T,N,0,T,typeof(axes)}(data, axes)
    end
    function WrappedArray(data::ArrayTypes{T,M}, axes::NTuple{N,Axis}) where {T,N,M}
        @assert(issubset(axisnames(data), axisnames(axes...)),
                "WrappedArray is missing axes: $(setdiff(axisnames(data), axisnames(axes...)))")
        new{T,N,M,typeof(data),typeof(axes)}(data, axes)
    end
end
WrappedArray(data, axes::Axis...) = WrappedArray(data, axes)
# Not all axes are Axis, since we insist, you have to provide a function that
# convert one for your iterator type by adding methods to this recursive
# _to_axis
WrappedArray(data, axes...) = WrappedArray(data, _to_axis(axes...)...)
_to_axis(ax::Axis, axes...) = (ax, _to_axis(axes...)...)
_to_axis() = ()

@inline Base.size(A::WrappedArray) = length.(A.axes)
@inline Base.size(A::WrappedArray, Ax::Axis) = length(A.axes[axisdim(A, Ax)])

# Slightly adapted axisdim, axisnames and axisvalues from AxisArrays

axisdim(A::WrappedArray, ax::Axis) = axisdim(A, typeof(ax))
@generated function axisdim(A::WrappedArray, ax::Type{Ax}) where Ax<:Axis
    dim = axisdim(A, Ax)
    :($dim)
end

# The actual computation is done in the type domain, which is a little tricky
# due to type invariance.
function axisdim(::Type{WrappedArray{T,N,M,D,Ax}}, ::Type{<:Axis{name,S} where S}) where {T,N,M,D,Ax,name}
    isa(name, Int) && return name <= N ? name : error("axis $name greater than array dimensionality $N")
    names = axisnames(Ax)
    idx = findfirst(isequal(name), names)
    isnothing(idx) && error("axis $name not found in array axes $names")
    idx
end

axisnames(::WrappedArray{T,N,M,D,Ax}) where {T,N,M,D,Ax}       = AxisArrays._axisnames(Ax)
axisnames(::Type{WrappedArray{T,N,M,D,Ax}}) where {T,N,M,D,Ax} = AxisArrays._axisnames(Ax)
axisvalues(A::WrappedArray) = axisvalues(A.axes...)

# Probably if there are colons or ranges on any of the virtual axes,
# we would instead want to give back a new WrappedArray, but this works well enough.
@generated function Base.getindex(A::WrappedArray{T,N,M,D,Ax}, idxs...) where {T,N,M,D,Ax}
    meta = Expr(:meta, :inline, :propagate_inbounds)
    if M == 0
        :(A.data)
    else
        names = axisnames(A)
        :($meta; A.data[$((Expr(:ref, :idxs, something(findfirst(isequal(n), names))) for n=axisnames(D))...)])
    end
end

# # RelIterators allow to use + Int, and - Int to get to adjacent indices, BUT
# # they are about a 100x slower, so should be used sparingly

# struct RelativeIndex{Ax}
#     ax::Ax
#     i::Int64
# end

# struct RelativeIterator{Ax}
#     ax::Ax
# end

# Base.start(it::RelativeIterator) = 1
# Base.done(it::RelativeIterator, state) = state == length(it.ax) + 1
# Base.next(it::RelativeIterator, state) = RelativeIndex(it.ax, state), state + 1
# Base.length(it::RelativeIterator) = length(it.ax)
# Base.:+(r::RelativeIndex, d::Int) = RelativeIndex(r.ax, r.i+d)
# Base.:-(r::RelativeIndex, d::Int) = RelativeIndex(r.ax, r.i-d)

# # We intercept the getindex calls of JuMPArray and AxisArray to resolve to the
# # actual index value
# for T = (:(JuMP.JuMPArray), :AxisArray)
#     @eval Base.getindex(a::$T, r::RelativeIndex) = a[r.ax[r.i]]
# end

# @inline unrelative(it) = it
# @inline unrelative(it::RelativeIterator) = it.ax

# _to_axis(it::RelativeIterator, axes...) = (unrelative(it), _to_axis(axes...)...)

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
#!/usr/bin/env julia6
using ReadGlobal, InplaceRealFFT

function resize(input::PaddedArray{T,3,L},newdim::NTuple{3,Int64}) where {T,L}
    u1_big = similar(input,newdim);
    fill!(complex(u1_big),0.0+0.0im)
    nx,ny,nz = size(input)
    nxb,nyb,nzb = size(u1_big)
    ratio = prod(newdim)/prod(size(real(input)))

    #get smallest one
    nxs = min(nx,nxb)
    nys = min(ny,nyb)
    nzs = min(nz,nzb)

    rfft!(input)
    @views begin
        u1_big[1:nxs,1:(div(nys,2)),1:(div(nzs,2))] .= ratio .* input[1:nxs,1:(div(nys,2)),1:(div(nzs,2))]
        u1_big[1:nxs,(nyb-div(nys,2)):end,1:(div(nzs,2))] .= ratio .* input[1:nxs,(ny-div(nys,2)):end,1:(div(nzs,2))]
        u1_big[1:nxs,(nyb-div(nys,2)):end,(nzb-div(nzs,2)):end] .= ratio .* input[1:nxs,(ny-div(nys,2)):end,(nz-div(nzs,2)):end]
        u1_big[1:nxs,1:(div(nys,2)),(nzb-div(nzs,2)):end] .= ratio .* input[1:nxs,1:div(nys,2),(nz-div(nzs,2)):end]
    end
    irfft!(u1_big)
    return u1_big
end

function resize(input::String,nnx::Int,nny::Int,nnz::Int)
    nx,ny,nz,_,_,_ = getdimsize()
    dtype,padded = checkinput(input,nx,ny,nz)
    us = PaddedArray{dtype}(input,(nx,ny,nz),padded)
    out = resize(us,(nnx,nny,nnz))
    write(input*"_resized",out)
end

resize(ARGS[1],parse(Int,ARGS[2]),parse(Int,ARGS[3]),parse(Int,ARGS[4]))


text:
module OmegaSoftPredicates

using SoftPredicates
using SoftPredicates: AbstractSoftBool
import OmegaCore
using OmegaCore: pw

export softconstraints,
       pwsofteq,
       pwsoftgt,
       pwsoftlt,
       ==̃ₚ,
       >̃ₚ,
       >=̃ₚ,
       <̃ₚ,
       <=̃ₚ

OmegaCore.initerror(::Type{T}) where T <: AbstractSoftBool = one(T)

softconstraints(x) = OmegaCore.condvar(x, DualSoftBool{Float64}) # zt: should this be fized to float64?

pwsofteq(x, y) = pw(softeq, x, y)
pwsoftgt(x, y) = pw(softgt, x, y)
pwsoftlt(x, y) = pw(softlt, x, y)

const ==̃ₚ = pwsofteq
const >̃ₚ = pwsoftgt
const >=̃ₚ = pwsoftgt
const <=̃ₚ = pwsoftlt
const <̃ₚ = pwsoftlt

export ≊
const ≊ = ==̃ₚ

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
export
    detailindex, detailrange, degree,
    scales, scalerange, levelmax, overlaplevelmax,
    dyadiclevelmax, todyadiclevel, hasdyadicfactor, isdyadic

#####
##### Detail and scale coefficient utilities
#####

"""
    detailindex(size, level, offset)

Return the detail coefficient for a `size`-length signal at `level` with `offset`.
"""
detailindex(size, level, offset) = round(Int, (size / degree(level)) + offset)

"""
    detailindex(signal::AbstractArray, level, offset)

Return the detail coefficient for `signal` at `level` with `offset`.
"""
detailindex(signal::AbstractArray, level, offset) = detailindex(first(size(signal)), level, offset)

"""
    detailindex(level, offset)

Return the detail coefficient for a dyadic signal at `level` with `offset`.
"""
detailindex(level, offset) = degree(level) + offset

"""
    detailrange(size::T, level::U) where {T<:Integer,U<:Integer}

Return the range of possible detail coefficients for
a `size`-length signal at `level`.
"""
function detailrange(size, level)
    lo = round(Int, size / degree(level) + one(level))
    hi = round(Int, size / degree(level - one(level)))
    return lo:hi
end

"""
    detailrange(signal::AbstractArray, level)

Return the range of possible detail coefficients for
the first dimension of a `signal` at `level`.
"""
detailrange(signal::AbstractArray, level) = detailrange(first(size(signal)), level)

"""
    detailrange(level)

Return the range of detail coefficients for a dyadic signal at a given `level`.
"""
detailrange(level) = (degree(level) + one(level)):(degree(level + one(level)))

"""
    degree(size, level)

Return the number of detail coefficients for
a `size`-length signal at `level`.
"""
degree(size, level) = round(Int, size / degree(level))

"""
    degree(signal::AbstractArray, level)

Return the number of detail coefficients for
the first dimension of a `signal` at `level`.
"""
degree(signal::AbstractArray, level) = degree(first(size(signal)), level)

"""
    degree(level)

Return the number of detail coefficients for a dyadic signal at a given `level`.
"""
degree(level) = one(level) << level

"""
    scales(signal::AbstractArray)

Return the number of scales in the first dimension of a dyadic `signal`.
"""
scales(signal::AbstractArray) = scales(first(size(signal)))

"""
    scales(level::Integer)

Return the number of scales for a dyadic signal at `level`.
"""
scales(level::Integer) = round(Int, log2(level))

"""
    scalerange(level)

Return the range of scaling coefficients for a dyadic signal at a given `level`.
"""
scalerange(level) = one(level):degree(level)

#####
##### Transformation level utilities
#####

"""
    levelmax(size)

Return the largest possible number of transformations
that can be performed on a `size`-length signal.
"""
function levelmax(size)
    size > 1 || return 0
    for power in Iterators.countfrom()
        hasdyadicfactor(size, power) || return power - 1
    end
end

"""
    levelmax(signal::AbstractArray)

Return the largest possible transformation level along the smallest axis of a `signal`.
"""
levelmax(signal::AbstractArray) = levelmax(minimum(size(signal)))

"""
    overlap_level_maximum(signal::AbstractArray)

Return the largest possible transformation level of a `signal`
when performing a maximal overlap transform.
"""
overlaplevelmax(signal::AbstractArray) = overlaplevelmax(length(signal))

"""
    overlap_level_maximum(size)

Return the largest possible transformation level of a
`size`-length signal when performing a maximal overlap transform.
"""
overlaplevelmax(size) = floor(Int, log2(size))

"""
    dyadiclevelmax(signal::AbstractArray)

Return the largest possible transformation level for a dyadic `signal`.
"""
dyadiclevelmax(signal::AbstractArray) = dyadiclevelmax(first(size(signal)))

"""
    dyadiclevelmax(size::Integer)

Return the largest possible transformation level for a dyadic `size`-length signal.
"""
dyadiclevelmax(size::Integer) = scales(size) - one(size)

"""
    todyadiclevel(size, level)

Given a `size`-length signal, return the corresponding
number of scales for a dyadic signal at `level`.
"""
todyadiclevel(size, level) = scales(size) - level

"""
    todyadiclevel(signal::AbstractArray, level)

Given a `signal`, return the corresponding
number of scales for a dyadic signal at `level`.
"""
todyadiclevel(signal::AbstractArray, level) = todyadiclevel(length(signal), level)

#####
##### Dyadic utilities
#####

"""
    isdyadic(signal::AbstractArray)

Return `true` if all dimensions of a `signal` are dyadic.
"""
isdyadic(signal::AbstractArray) = all(isdyadic, size(signal))

"""
    isdyadic(size::T) where {T<:Integer}

Return `true` if a `size`-length signal is a power of two.
"""
isdyadic(size) = size == one(size) << scales(size)

"""
    hasdyadicfactor(signal::AbstractArray, power)

Return `true` if all dimensions of a `signal` are multiples of the `level`-th power of two.
"""
hasdyadicfactor(signal::AbstractArray, power) = all(hasdyadicfactor(axis, power) for axis in size(signal))

"""
    hasdyadicfactor(size, power)

Return `true` if `size` is a multiple of the `level`-th power of two.
"""
hasdyadicfactor(size, power) = iszero(size % degree(power))


text:
"""
*load_graph_from_csv*

This function loads a graph from csv file containing the adjacency matrix of the graph. The file location and name is passed as argument (location/name.csv) together with the separator used in the csv e.g. ';' or ','.
"""
function load_graph_from_csv(csv_adj_m::String;sep::Char=';')

    #adjacency matrix
    am = readdlm(csv_adj_m, sep, Int, '\n')

    #create weighted graph
    wg = SimpleWeightedGraph(am)
    #create metagraph
    g = MetaGraph(am)

    [ set_prop!(g, r, c, :weight, SimpleWeightedGraphs.weights(wg)[r,c]) for r=1:size(SimpleWeightedGraphs.weights(wg),1),c=1:size(SimpleWeightedGraphs.weights(wg),2) if SimpleWeightedGraphs.weights(wg)[r,c] >0]
    return g
end

"""
It annotates the graph with the coordinates obtained from the passed file name and location (location/name.csv). Annotating the graph is just adding attributes to the MetaGraph object passed.
"""
function add_geo_coordinates_to_graph!(g::MetaGraph,csv_coordinates::String)
    #coordinates
    df_net = CSV.File(csv_coordinates) |> DataFrame
    df_net = sort(df_net,[:id])
    #populate metagraph
    vi = 0
    for r in eachrow(df_net)
        vi += 1
        set_props!(g,vi,Dict(:lon=>r.lon,:lat=>r.lat))
    end
    return g
end

"""
Create a packet using arguments
"""
function create_pkt(src::Int64,dst::Int64,model)
    model.pkt_id += 1
    return DPacket(model.pkt_id,src,dst,model.:pkt_size,model.:ticks,100)
end

"""
Traffic generation per tick
#TODO Implement this function as per required
"""
function generate_traffic!(model)
    #random pkts
    traffic_μ = first(model.traffic_dist_params)
    traffic_sd = last(model.traffic_dist_params)
    q_pkts = abs(round(model.traffic_packets*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd))))
    # q_pkts: A percentage of the model.pkt_per_tick so NEs are able to process traffic coming from different nodes (NEs)
    #src,dst = samplepair(1:nv(model.ntw_graph)) # can be replaced for random pair
    pairs = model.ntw_services
    ##**For testing only**
    #fixed pkts
    #q_pkts = 5
    # pairs =[(5,14)]
    tpkts = 400
    println("$(model.seed)-->[$(model.ticks)] - generating traffic btwn $pairs")
    for p in pairs
        src,dst = p
        #for bipartite sim
        if src == 13
            q_pkts = @match model.traffic_packets begin
                #model.traffic_packets
                100 => abs(round(tpkts÷4*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) # 75% drop
                200 => abs(round(tpkts÷2*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) # 50% drop
                300 => abs(round((tpkts÷4)*3*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) #25% drop

                400 => abs(round(tpkts*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) #baseline
                500 => abs(round((tpkts÷4)*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) #25% drop
                600 => abs(round((tpkts÷2)*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) #50% drop
                700 => abs(round((tpkts÷4)*3*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) # 75% drop
                800 => abs(round(tpkts*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) #100% drop
            end

        elseif src == 7

            q_pkts = @match model.traffic_packets begin
                100 => abs(round(tpkts*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) #irrelevant
                200 => abs(round(tpkts*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) #irrelevant
                300 => abs(round(tpkts*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) #irrelevant

                400 => abs(round(tpkts*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) #baseline
                500 => abs(round(tpkts*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) #25% drop
                600 => abs(round(tpkts*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) #50% drop
                700 => abs(round(tpkts*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) #75% drop
                800 => abs(round(tpkts*get_random(model.seed,model.ticks,Normal(traffic_μ,traffic_sd)))) #100% drop
            end

        end

        sne_src = getindex(model,src)
        sne_dst = getindex(model,dst)
        if is_up(sne_src) && is_up(sne_dst)
            # log_info(model.ticks,sne_src.id,"pkts generate: $(q_pkts) --> $(length(sne_src.queue.data))")
            if q_pkts > sne_src.queue.sz_max - length(sne_src.queue.data)  - Int(round(sne_src.queue.sz_max * .2))
                log_info(model.ticks,sne_src.id,"Traffic generated in origin ($(q_pkts)) greater than current node queue capacity ($(sne_src.queue.sz_max - sne_src.queue.sz_max - length(sne_src.queue.data))) reducing pkts in origin to max. ")
                q_pkts = sne_src.queue.sz_max - length(sne_src.queue.data) - Int(round(sne_src.queue.sz_max * .2))
            end
            for i =1:q_pkts
                pkt = create_pkt(src,dst,model)
                # if model.ticks >= 80
                    # log_info(model.ticks,src,3,"Sending src: $src - dst: $dst -> q_pkts: $q_pkts ==> $pkt packets ")
                # end
                push_msg!(sne_src,OFMessage(next_ofmid!(model), model.ticks,src,0,pkt)) # always from port 0
            end
        end
    end
    # log_info(model.ticks,"end generated ==> $(length(getindex(model,17).queue.data))")
   # log_info("[$(model.ticks)] $(q_pkts) pkts generated")
end

"""
 Get a random number of packets to be processed by a sne on a given tick
 defined by the sequence passed (tick + sne_id)
"""
function get_random_packets_to_process(seed,sequence,max_ppt)
    #max pkts processed per tick
    #1/10 of ppt
    ppt_u = Int(round(max_ppt/10))
    return get_random(seed,sequence,((max_ppt-(4*ppt_u)):ppt_u:max_ppt))
end

"""
Rejoin node to network according to
- base_network
- current_network
- rejoining node id
"""
function rejoin_node!(model,rjn_id::Int64)

   base_g = model.base_ntw_graph
   g = model.ntw_graph

   nbs = neighbors(base_g,rjn_id)
   lsnes = get_live_snes(model)
   #those neighbour nodes that are up
   sne_ids = intersect(nbs,lsnes)

   for sne_id in sne_ids
       add_edge!(g,sne_id,rjn_id)
       add_edge!(g,rjn_id,sne_id)
       sne = getindex(model,sne_id)
       link_up!(sne,rjn_id,model)
   end

   rjn_ag = getindex(model,rjn_id)
   set_up!(rjn_ag)

   #re init ports
   nbs = all_neighbors(model.base_ntw_graph,get_address(rjn_ag.id,model.base_ntw_graph))

   if !in((0,"h$(rjn_ag.id)"),get_state(rjn_ag).port_edge_list)
        push_ep_entry!(rjn_ag,(0,"h$(rjn_ag.id)")) # link to a host of the same id
   end

#    log_info(model.ticks,rjn_id," Rejoining adding $(get_state(rjn_ag).port_edge_list)")
   # creates entries for all nb from current ntw?
   for i in 1:size(nbs,1)
       if is_up(getindex(model,nbs[i]))
          push_ep_entry!(rjn_ag,(i,"s$(nbs[i])"))
       end
   end

   #it simulates control detects sne up:
   # when down, controller id is multiplied by -1, so we do the opposite now
   aid = abs(get_control_agent(rjn_id,model))
   a = getindex(model,aid)
   controlled_sne_up!(a,rjn_id,sne_ids,model)
end

export load_graph_from_csv


text:
# AUTO GENERATED FILE - DO NOT EDIT

export breadcrumbs

"""
    breadcrumbs(;kwargs...)

A Breadcrumbs component.
Separate list of react nodes with given separator. For more information, see: https://mantine.dev/core/breadcrumbs/
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `items` (optional): Link items. items has the following type: Array of lists containing elements 'title', 'href'.
Those elements have the following types:
  - `title` (String; required)
  - `href` (String; required)s
- `separator` (String; optional): Separator between breadcrumbs
- `style` (Dict; optional): Inline style override
"""
function breadcrumbs(; kwargs...)
        available_props = Symbol[:id, :items, :separator, :style]
        wild_props = Symbol[]
        return Component("breadcrumbs", "Breadcrumbs", "dash_mantine_components", available_props, wild_props; kwargs...)
end


text:
using Glutils.Math

title = "Basic Example"
background = "#112"
vsync = false

struct Vertex
   position::Vec{2}
   color::Vec{3}
end

const vertices = @SVector [
  Vertex([ 0.5,  0.5], [0, 0, 1]),
  Vertex([ 0.5, -0.5], [1, 1, 0]),
  Vertex([-0.5, -0.5], [1, 0, 0]),
  Vertex([-0.5,  0.5], [0, 1, 0])
]

const indices = @SVector Cuint[
   0, 1, 3,
   1, 2, 3
]

vao = createvao()
ibo, vbo = createbuffers(2)

glVertexArrayElementBuffer(vao, ibo)
glVertexArrayVertexBuffer(vao, 0, vbo, 0, sizeof(Vertex))

bufferdata!(ibo, indices,  :static_draw)
bufferdata!(vbo, vertices, :static_draw)
attrformat!(vao, 0, Vertex)

program = create_program(
   create_shader("shaders/simple.vert"),
   create_shader("shaders/simple.frag"),
)

@show active_uniforms(program)

function onupdate(dt)
   bindvao!(vao)
   bindprogram!(program)

   left, right  = -10, 10
   bottom, top  = -10, 10
   near, far = 0.0001, 10000

   eye = Vec(0, 0, 20)
   center = Vec(0, 0, 0)
   up = Vec(0, 1, 0)

   @uniforms program {
      uTime  = t
      uView  = lookAt(eye, center, up)
      uModel = scale(10) * rotateZ(2t / π)
      uProj  = ortho(left, right, bottom, top, near, far)
   }

   draw_elements(6, Cuint)
end

max_stars_repo_name: jorge-brito/Glutils

max_issues_repo_name: jorge-brito/Glutils

max_forks_repo_name: jorge-brito/Glutils


text:
using VisualCircuits, Compose

begin
    notQ1 = Symbol("~Q1")
    notQ2 = Symbol("~Q2")

    # shift register
    IC7496(x, y) = IC(x, y, "7496", [
        :CLR    :QA     :QB     :QC     :GND    :QD     :QE     :SER;
        :CLK    :A      :B      :C      :VCC    :D      :E      :PRE
    ])
    # full adder
    IC7483(x, y) = IC(x, y, "7483", [
        :B4     :Q4     :C4     :C0     :GND    :B1     :A1     :Q1;
        :A4     :Q3     :A3     :B3     :VCC    :Q2     :B2     :A2
    ])
    # D-flipflop
    IC7474(x, y) = IC(x, y, "7474", [
        :VCC    :CLR2   :D2     :CLK2   :PRE2   :Q2     notQ2;
        :CLR1   :D1     :CLK1   :PRE1   :Q1     notQ1   :GND
    ])

    width, height = 18, 9.6
    le, re = 0, width
    reg1 = IC7496(4, 2.7)
    reg2 = IC7496(5.65, 7.1)
    add  = IC7483(9.2, 3.1)
    dff  = IC7474(10.45, 7.1)
    reg3 = IC7496(14.5, 3.1)

    # guide (not displayed)
    box = compose(
        context(),
        rectangle(0, 0, width, height),
        stroke("cyan"), fill(nothing)
    )

    # power
    y1 = 0.2
    y2 = 4.6
    y3 = 9.2
    power = compose(
        context(),
        connect((le, y1), (re, y1)), connect((le, y2), (re, y2)),
        connect((le, y3), (re, y3)),
        connect(reg1[:VCC], y2, :y), connect(reg2[:VCC], y3, :y),
        connect(reg3[:VCC], y2, :y), connect(add[:VCC], y2, :y),
        connect(dff[:VCC], y2, :y),
        connect(dff[:PRE2], y2, :y),
        dots(
            (reg1[:VCC][1], y2), (reg2[:VCC][1], y3),
            (reg3[:VCC][1], y2), (add[:VCC][1], y2),
            (dff[:VCC][1], y2),  (dff[:PRE2][1], y2),
            color = "darkred"
        ),
        stroke("darkred")
    )

    y1 = 0.4
    y2 = 4.8
    y3 = 9.4
    y4 = 4.2
    ground = compose(
        context(),
        connect((le, y1), (re, y1)), connect((le, y2), (re, y2)),
        connect((le, y3), (re, y3)),
        connect(reg1[:GND], y1, :y), connect(reg2[:GND], y2, :y),
        connect(reg3[:GND], y1, :y), connect(add[:GND], y1, :y),
        connect(dff[:GND], y3, :y),
        connect(add[:B2], y2, :y),
        connect(add[:A2], (add[:A2][1], y4), (add[:B2][1], y4)),
        connect(reg1[:SER], y1, :y),
        connect(reg2[:SER], y2, :y),
        connect(reg3[:PRE], y2, :y),
        dots(
            (reg1[:GND][1], y1), (reg2[:GND][1], y2),
            (reg3[:GND][1], y1), (add[:GND][1], y1),
            (dff[:GND][1], y3),  (add[:B2][1], y4),
            (add[:B2][1], y2),   (reg1[:SER][1], y1),
            (reg2[:SER][1], y2), (reg3[:PRE][1], y2),
            color = "darkblue"
        ),
        stroke("darkblue")
    )

    # input
    y1 = 3.8
    y2 = 8.2
    inputs = compose(
        context(),
        connect(reg1[:A], y1, :y), connect(reg2[:A], y2, :y),
        connect(reg1[:B], y1, :y), connect(reg2[:B], y2, :y),
        connect(reg1[:C], y1, :y), connect(reg2[:C], y2, :y),
        connect(reg1[:D], y1, :y), connect(reg2[:D], y2, :y),
        connect(reg1[:PRE], y1, :y), connect(reg2[:PRE], y2, :y),
        dots(
            (reg1[:A][1], y1),  (reg2[:A][1], y2),
            (reg1[:B][1], y1),  (reg2[:B][1], y2),
            (reg1[:C][1], y1),  (reg2[:C][1], y2),
            (reg1[:D][1], y1),  (reg2[:D][1], y2),
            (reg1[:PRE][1], y1),  (reg2[:PRE][1], y2),
            color = "darkgreen"
        ),
        stroke("darkgreen")
    )

    # input leds
    y1 = 1.6
    y2 = 6
    inputLEDs = compose(
        context(),
        connect(reg1[:QA], y1, :y), connect(reg2[:QA], y2, :y),
        connect(reg1[:QB], y1, :y), connect(reg2[:QB], y2, :y),
        connect(reg1[:QC], y1, :y), connect(reg2[:QC], y2, :y),
        connect(reg1[:QD], y1, :y), connect(reg2[:QD], y2, :y),
        dots(
            (reg1[:QA][1], y1),  (reg2[:QA][1], y2),
            (reg1[:QB][1], y1),  (reg2[:QB][1], y2),
            (reg1[:QC][1], y1),  (reg2[:QC][1], y2),
            (reg1[:QD][1], y1),  (reg2[:QD][1], y2),
            color = "orange"
        ),
        stroke("orange")
    )

    # input text
    y1 = 1.6
    y2 = 3.8
    y3 = 6
    y4 = 8.2
    input_labels = compose(
        context(),
        text(reg1[:QA][1], y1-0.2, "b1", hcenter, vbottom),
        text(reg1[:QB][1], y1-0.2, "b2", hcenter, vbottom),
        text(reg1[:QC][1], y1-0.2, "b3", hcenter, vbottom),
        text(reg1[:QD][1]+0.2, y1, "b4", hleft, vcenter),

        text(reg1[:A][1], y2+0.2, "b1", hcenter, vtop),
        text(reg1[:B][1], y2+0.2, "b2", hcenter, vtop),
        text(reg1[:C][1], y2+0.2, "b3", hcenter, vtop),
        text(reg1[:D][1], y2+0.2, "b4", hcenter, vtop),

        text(reg2[:QA][1], y3-0.2, "b'1", hcenter, vbottom),
        text(reg2[:QB][1], y3-0.2, "b'2", hcenter, vbottom),
        text(reg2[:QC][1], y3-0.2, "b'3", hcenter, vbottom),
        text(reg2[:QD][1]+0.2, y3, "b'4", hleft, vcenter),

        text(reg2[:A][1], y4+0.2, "b'1", hcenter, vtop),
        text(reg2[:B][1], y4+0.2, "b'2", hcenter, vtop),
        text(reg2[:C][1], y4+0.2, "b'3", hcenter, vtop),
        text(reg2[:D][1], y4+0.2, "b'4", hcenter, vtop),

        text(reg1[:PRE][1], y2+0.2, "SET", hcenter, vtop),
        text(reg2[:PRE][1], y4+0.2, "SET", hcenter, vtop),
        fontsize(20pt), font("Helvetica-Bold"), linewidth(0mm)
    )

    # output
    y1 = 2.0
    outputs = compose(
        context(),
        connect(reg3[:QA], y1, :y), connect(reg3[:QB], y1, :y),
        connect(reg3[:QC], y1, :y), connect(reg3[:QD], y1, :y),
        connect(reg3[:QE], y1, :y),
        dots(
            (reg3[:QA][1], y1), (reg3[:QB][1], y1),
            (reg3[:QC][1], y1), (reg3[:QD][1], y1),
            (reg3[:QE][1], y1), color = "orange"
        ),
        stroke("orange")
    )

    output_labels = compose(
        context(),
        text(reg3[:QA][1], y1-0.2, "B1", hcenter, vbottom),
        text(reg3[:QB][1], y1-0.2, "B2", hcenter, vbottom),
        text(reg3[:QC][1], y1-0.2, "B3", hcenter, vbottom),
        text(reg3[:QD][1], y1-0.2, "B4", hcenter, vbottom),
        text(reg3[:QE][1], y1-0.2, "B5", hcenter, vbottom),
        fontsize(20pt), font("Helvetica-Bold"), linewidth(0mm)
    )

    # clock
    x1 = 1
    x2 = 8.2
    y1 = 8.8
    y2 = 6
    clock = compose(
        context(),
        connect((x1, y1), (reg3[:CLK][1], y1), reg3[:CLK]),
        connect(dff[:CLK2], (dff[:CLK2][1], y2), (x2, y2), (x2, y1)),
        connect(reg1[:CLK], y1, :y),
        connect(reg2[:CLK], y1, :y),
        stroke("darkcyan")
    )

    # clear
    x1 = 1
    x2 = 1.75
    y1 = 0.8
    y2 = 5.2
    y3 = 5.6
    clear = compose(
        context(),
        connect((x1, y1), (reg3[:CLR][1], y1), reg3[:CLR]),
        connect(reg1[:CLR], y1, :y),
        connect(reg2[:CLR], (reg2[:CLR][1], y2), (x2, y2), (x2, y1)),
        connect(dff[:CLR2], (dff[:CLR2][1], y2), (reg2[:CLR][1], y2)),
        dots((reg2[:CLR][1], y2), color = "purple"),
        stroke("purple")
    )

    # d flip flop
    x1 = 6.8
    x2 = 6.4
    y1 = 5.6
    y2 = 2
    y3 = 1.2
    y4 = 1.6
    y5 = 6
    wires = compose(
        context(),
        connect(dff[:D2], add[:Q2]),
        connect(
            dff[:Q2], (dff[:Q2][1], y1), (x1, y1),
            (x1, y2), (add[:C0][1], y2), add[:C0]
        ),
        connect(add[:Q1], (add[:Q1][1], y3), (reg3[:SER][1], y3), reg3[:SER]),
        connect(add[:A1], (add[:A1][1], y3), (reg1[:QD][1], y3), reg1[:QD]),
        connect(add[:B1], (add[:B1][1], y4), (x2, y4), (x2, y5), (reg2[:QD][1], y5))
    )

    x = 1
    y1 = 0.8
    y2 = 8.8
    labels = compose(
        context(),
        text(x-0.2, y1, "CLR", hright, vcenter),
        text(x-0.2, y2, "CLK", hright, vcenter),
        fontsize(20pt), font("Helvetica-Bold"), linewidth(0mm)
    )

    circuit = compose(
        context(units = UnitBox(0, 0, width, height)),
        render(reg1), render(reg2), render(reg3), render(add), render(dff),

        # box,
        power, ground,
        inputs, inputLEDs, input_labels,
        outputs, output_labels,
        clock, clear,
        wires, labels,

        stroke("black"), linewidth(2mm)
    )

    circuit |> SVG("4bit_adder.svg", autosize(width, height)...)
    # circuit |> PDF("4bit_adder.pdf", autosize(width, height)...)
end


text:
"""
    Normal(μ,σ)

The *Normal distribution* with mean `μ` and standard deviation `σ≥0` has probability density function

```math
f(x; \\mu, \\sigma) = \\frac{1}{\\sqrt{2 \\pi \\sigma^2}}
\\exp \\left( - \\frac{(x - \\mu)^2}{2 \\sigma^2} \\right)
```

Note that if `σ == 0`, then the distribution is a point mass concentrated at `μ`.
Though not technically a continuous distribution, it is allowed so as to account for cases where `σ` may have underflowed,
and the functions are defined by taking the pointwise limit as ``σ → 0``.

```julia
Normal()          # standard Normal distribution with zero mean and unit variance
Normal(μ)         # Normal distribution with mean μ and unit variance
Normal(μ, σ)      # Normal distribution with mean μ and variance σ^2

params(d)         # Get the parameters, i.e. (μ, σ)
mean(d)           # Get the mean, i.e. μ
std(d)            # Get the standard deviation, i.e. σ
```

External links

* [Normal distribution on Wikipedia](http://en.wikipedia.org/wiki/Normal_distribution)

"""
struct Normal{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    σ::T
    Normal{T}(µ::T, σ::T) where {T<:Real} = new{T}(µ, σ)
end

function Normal(μ::T, σ::T; check_args=true) where {T <: Real}
    check_args && @check_args(Normal, σ >= zero(σ))
    return Normal{T}(μ, σ)
end

#### Outer constructors
Normal(μ::Real, σ::Real) = Normal(promote(μ, σ)...)
Normal(μ::Integer, σ::Integer) = Normal(float(μ), float(σ))
Normal(μ::T) where {T <: Real} = Normal(μ, one(T))
Normal() = Normal(0.0, 1.0, check_args=false)

const Gaussian = Normal

# #### Conversions
convert(::Type{Normal{T}}, μ::S, σ::S) where {T <: Real, S <: Real} = Normal(T(μ), T(σ))
convert(::Type{Normal{T}}, d::Normal{S}) where {T <: Real, S <: Real} = Normal(T(d.μ), T(d.σ), check_args=false)

@distr_support Normal -Inf Inf

#### Parameters

params(d::Normal) = (d.μ, d.σ)
@inline partype(d::Normal{T}) where {T<:Real} = T

location(d::Normal) = d.μ
scale(d::Normal) = d.σ

Base.eltype(::Type{Normal{T}}) where {T} = T

#### Statistics

mean(d::Normal) = d.μ
median(d::Normal) = d.μ
mode(d::Normal) = d.μ

var(d::Normal) = abs2(d.σ)
std(d::Normal) = d.σ
skewness(d::Normal{T}) where {T<:Real} = zero(T)
kurtosis(d::Normal{T}) where {T<:Real} = zero(T)

entropy(d::Normal) = (log2π + 1)/2 + log(d.σ)

function kldivergence(p::Normal, q::Normal)
    μp = mean(p)
    σ²p = var(p)
    μq = mean(q)
    σ²q = var(q)
    σ²p_over_σ²q = σ²p / σ²q
    return (abs2(μp - μq) / σ²q - logmxp1(σ²p_over_σ²q)) / 2
end

#### Evaluation

# Helpers
"""
    xval(d::Normal, z::Real)

Computes the x-value based on a Normal distribution and a z-value.
"""
function xval(d::Normal, z::Real)
    if isinf(z) && iszero(d.σ)
        d.μ + one(d.σ) * z
    else
        d.μ + d.σ * z
    end
end
"""
    zval(d::Normal, x::Real)

Computes the z-value based on a Normal distribution and a x-value.
"""
zval(d::Normal, x::Real) = (x - d.μ) / d.σ

gradlogpdf(d::Normal, x::Real) = -zval(d, x) / d.σ

# logpdf
_normlogpdf(z::Real) = -(abs2(z) + log2π)/2

function logpdf(d::Normal, x::Real)
    μ, σ = d.μ, d.σ
    if iszero(d.σ)
        if x == μ
            z = zval(Normal(μ, one(σ)), x)
        else
            z = zval(d, x)
            σ = one(σ)
        end
    else
        z = zval(Normal(μ, σ), x)
    end
    return _normlogpdf(z) - log(σ)
end

# pdf
_normpdf(z::Real) = exp(-abs2(z)/2) * invsqrt2π

function pdf(d::Normal, x::Real)
    μ, σ = d.μ, d.σ
    if iszero(σ)
        if x == μ
            z = zval(Normal(μ, one(σ)), x)
        else
            z = zval(d, x)
            σ = one(σ)
        end
    else
        z = zval(Normal(μ, σ), x)
    end
    return _normpdf(z) / σ
end

# logcdf
function _normlogcdf(z::Real)
    if z < -one(z)
        return log(erfcx(-z * invsqrt2)/2) - abs2(z)/2
    else
        return log1p(-erfc(z * invsqrt2)/2)
    end
end

function logcdf(d::Normal, x::Real)
    if iszero(d.σ) && x == d.μ
        z = zval(Normal(zero(d.μ), d.σ), one(x))
    else
        z = zval(d, x)
    end
    return _normlogcdf(z)
end

# logccdf
function _normlogccdf(z::Real)
    if z > one(z)
        return log(erfcx(z * invsqrt2)/2) - abs2(z)/2
    else
        return log1p(-erfc(-z * invsqrt2)/2)
    end
end

function logccdf(d::Normal, x::Real)
    if iszero(d.σ) && x == d.μ
        z = zval(Normal(zero(d.μ), d.σ), one(x))
    else
        z = zval(d, x)
    end
    return _normlogccdf(z)
end

# cdf
_normcdf(z::Real) = erfc(-z * invsqrt2)/2

function cdf(d::Normal, x::Real)
    if iszero(d.σ) && x == d.μ
        z = zval(Normal(zero(d.μ), d.σ), one(x))
    else
        z = zval(d, x)
    end
    return _normcdf(z)
end

# ccdf
_normccdf(z::Real) = erfc(z * invsqrt2)/2

function ccdf(d::Normal, x::Real)
    if iszero(d.σ) && x == d.μ
        z = zval(Normal(zero(d.μ), d.σ), one(x))
    else
        z = zval(d, x)
    end
    return _normccdf(z)
end

# quantile
function quantile(d::Normal, p::Real)
    # Promote to ensure that we don't compute erfcinv in low precision and then promote
    _p, _μ, _σ = map(float, promote(p, d.μ, d.σ))
    q = xval(d, -erfcinv(2*_p) * sqrt2)
    if isnan(_p)
        return oftype(q, _p)
    elseif iszero(_σ)
        # Quantile not uniquely defined at p=0 and p=1 when σ=0
        if iszero(_p)
            return oftype(q, -Inf)
        elseif isone(_p)
            return oftype(q, Inf)
        else
            return oftype(q, _μ)
        end
    end
    return q
end

# cquantile
function cquantile(d::Normal, p::Real)
    # Promote to ensure that we don't compute erfcinv in low precision and then promote
    _p, _μ, _σ = map(float, promote(p, d.μ, d.σ))
    q = xval(d, erfcinv(2*_p) * sqrt2)
    if isnan(_p)
        return oftype(q, _p)
    elseif iszero(_σ)
        # Quantile not uniquely defined at p=0 and p=1 when σ=0
        if iszero(_p)
            return oftype(q, Inf)
        elseif isone(_p)
            return oftype(q, -Inf)
        else
            return oftype(q, _μ)
        end
    end
    return q
end

mgf(d::Normal, t::Real) = exp(t * d.μ + d.σ^2 / 2 * t^2)
cf(d::Normal, t::Real) = exp(im * t * d.μ - d.σ^2 / 2 * t^2)

#### Affine transformations

Base.:+(d::Normal, c::Real) = Normal(d.μ + c, d.σ)
Base.:*(c::Real, d::Normal) = Normal(c * d.μ, abs(c) * d.σ)

#### Sampling

rand(rng::AbstractRNG, d::Normal{T}) where {T} = d.μ + d.σ * randn(rng, float(T))

#### Fitting

struct NormalStats <: SufficientStats
    s::Float64    # (weighted) sum of x
    m::Float64    # (weighted) mean of x
    s2::Float64   # (weighted) sum of (x - μ)^2
    tw::Float64    # total sample weight
end

function suffstats(::Type{<:Normal}, x::AbstractArray{T}) where T<:Real
    n = length(x)

    # compute s
    s = x[1]
    for i = 2:n
        @inbounds s += x[i]
    end
    m = s / n

    # compute s2
    s2 = abs2(x[1] - m)
    for i = 2:n
        @inbounds s2 += abs2(x[i] - m)
    end

    NormalStats(s, m, s2, n)
end

function suffstats(::Type{<:Normal}, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Real
    n = length(x)

    # compute s
    tw = w[1]
    s = w[1] * x[1]
    for i = 2:n
        @inbounds wi = w[i]
        @inbounds s += wi * x[i]
        tw += wi
    end
    m = s / tw

    # compute s2
    s2 = w[1] * abs2(x[1] - m)
    for i = 2:n
        @inbounds s2 += w[i] * abs2(x[i] - m)
    end

    NormalStats(s, m, s2, tw)
end

# Cases where μ or σ is known

struct NormalKnownMu <: IncompleteDistribution
    μ::Float64
end

struct NormalKnownMuStats <: SufficientStats
    μ::Float64      # known mean
    s2::Float64     # (weighted) sum of (x - μ)^2
    tw::Float64     # total sample weight
end

function suffstats(g::NormalKnownMu, x::AbstractArray{T}) where T<:Real
    μ = g.μ
    s2 = abs2(x[1] - μ)
    for i = 2:length(x)
        @inbounds s2 += abs2(x[i] - μ)
    end
    NormalKnownMuStats(g.μ, s2, length(x))
end

function suffstats(g::NormalKnownMu, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Real
    μ = g.μ
    s2 = abs2(x[1] - μ) * w[1]
    tw = w[1]
    for i = 2:length(x)
        @inbounds wi = w[i]
        @inbounds s2 += abs2(x[i] - μ) * wi
        tw += wi
    end
    NormalKnownMuStats(g.μ, s2, tw)
end

struct NormalKnownSigma <: IncompleteDistribution
    σ::Float64

    function NormalKnownSigma(σ::Float64)
        σ > 0 || throw(ArgumentError("σ must be a positive value."))
        new(σ)
    end
end

struct NormalKnownSigmaStats <: SufficientStats
    σ::Float64      # known std.dev
    sx::Float64      # (weighted) sum of x
    tw::Float64     # total sample weight
end

function suffstats(g::NormalKnownSigma, x::AbstractArray{T}) where T<:Real
    NormalKnownSigmaStats(g.σ, sum(x), Float64(length(x)))
end

function suffstats(g::NormalKnownSigma, x::AbstractArray{T}, w::AbstractArray{T}) where T<:Real
    NormalKnownSigmaStats(g.σ, dot(x, w), sum(w))
end

# fit_mle based on sufficient statistics

fit_mle(::Type{<:Normal}, ss::NormalStats) = Normal(ss.m, sqrt(ss.s2 / ss.tw))
fit_mle(g::NormalKnownMu, ss::NormalKnownMuStats) = Normal(g.μ, sqrt(ss.s2 / ss.tw))
fit_mle(g::NormalKnownSigma, ss::NormalKnownSigmaStats) = Normal(ss.sx / ss.tw, g.σ)

# generic fit_mle methods

function fit_mle(::Type{<:Normal}, x::AbstractArray{T}; mu::Float64=NaN, sigma::Float64=NaN) where T<:Real
    if isnan(mu)
        if isnan(sigma)
            fit_mle(Normal, suffstats(Normal, x))
        else
            g = NormalKnownSigma(sigma)
            fit_mle(g, suffstats(g, x))
        end
    else
        if isnan(sigma)
            g = NormalKnownMu(mu)
            fit_mle(g, suffstats(g, x))
        else
            Normal(mu, sigma)
        end
    end
end

function fit_mle(::Type{<:Normal}, x::AbstractArray{T}, w::AbstractArray{Float64}; mu::Float64=NaN, sigma::Float64=NaN) where T<:Real
    if isnan(mu)
        if isnan(sigma)
            fit_mle(Normal, suffstats(Normal, x, w))
        else
            g = NormalKnownSigma(sigma)
            fit_mle(g, suffstats(g, x, w))
        end
    else
        if isnan(sigma)
            g = NormalKnownMu(mu)
            fit_mle(g, suffstats(g, x, w))
        else
            Normal(mu, sigma)
        end
    end
end


text:
using Test
using RayTraceEllipsoids, LinearAlgebra, CoordinateTransformations

@test all(1:10^5) do _
    oz, l, dz = rand(4) .- 0.5
    positive = rand(Bool)
    a = RayTraceEllipsoids.isgoodz(oz, l, dz, positive)
    z = oz + l*dz
    b = z > 0 ? positive : !positive
    a == b
end

@test all(1:10^5) do _
    z = rand()
    orig = V3(0.,0,z + 1)
    dir = V3(0.,0,-1)
    pos = true
    RayTraceEllipsoids.distance(orig, dir, pos) == z
end

@test RayTraceEllipsoids.distance(V3(2,2,2.), V3(0,0,-1.), true) ≡ nothing

@test all(1:10^5) do _
    z = rand()
    c = V3(1.,1,1)
    r = V3(1.,2,3)
    h = true
    s = Ellipsoid(c, r, h)
    orig = V3(1,1,c[3] + r[3] + z)
    dir = V3(0.,0,-1)
    r = Ray(orig, dir)
    RayTraceEllipsoids.distance(r, s) ≈ z
end

@test all(1:10^5) do _
    z = rand()
    c = V3(1.0,1,1)
    r = V3(1.0,2,3)
    h = true
    s = Ellipsoid(c, r, h)
    orig = V3(1.0,1,c[3] + r[3] + z)
    dir = V3(0.0,0,-1)
    r = Ray(orig, dir)
    RayTraceEllipsoids.distance(r, s) ≈ z
end

@test all(1:10^5) do _
    orig = V3(0.,0,0)
    l = rand()
    dir = V3(0.,0,-1)
    o = V3(rand(3)...)
    ucs = Translation(o)
    rorig = o
    RayTraceEllipsoids._getl(orig, l, dir, ucs, rorig) == l
end

#1,3
r = Ray(V3(2.,2,2), V3(0.,0,-1))
s = Cylinder(1.5,-2.,2.)
L = 4.
@test RayTraceEllipsoids.incylinder(r, s, L) == 0
#2
r = Ray(V3(1.,1,2), V3(0.,0,-1))
s = Cylinder(1.5,-2.,2.)
L = 4.
@test RayTraceEllipsoids.incylinder(r, s, L) == 4
#4
r = Ray(V3(0.,0,2), V3(sqrt(1/2),0,-sqrt(1/2)))
s = Cylinder(1.,-2.,2.)
L = sqrt(4^2+4^2)
@test RayTraceEllipsoids.incylinder(r, s, L) == sqrt(2)
#5
r = Ray(V3(1.,-1,2), V3(0.,1,-1))
s = Cylinder(1.,-2.,2.)
L = sqrt(4^2+4^2)
@test RayTraceEllipsoids.incylinder(r, s, L) == 0
#6
r = Ray(V3(-4tand(45),0,2), V3(sqrt(1/2),0,-sqrt(1/2)))
s = Cylinder(1.,-2.,2.)
L = sqrt(4^2+4^2)
@test RayTraceEllipsoids.incylinder(r, s, L) ≈ sqrt(2)
#7
r = Ray(V3(-2tand(45),0,2), V3(sqrt(1/2),0,-sqrt(1/2)))
s = Cylinder(1.,-2.,2.)
L = sqrt(4^2+4^2)
@test RayTraceEllipsoids.incylinder(r, s, L) ≈ sqrt(2^2+2^2)
#8
r = Ray(V3(-6tand(45),0,2), V3(sqrt(1/2),0,-sqrt(1/2)))
s = Cylinder(1.,-2.,2.)
L = sqrt(4^2+4^2)
@test RayTraceEllipsoids.incylinder(r, s, L) == 0

@test RayTraceEllipsoids.absorption(0, 0) == 0
@test RayTraceEllipsoids.absorption(1, 0) == 0
@test RayTraceEllipsoids.absorption(0, 1) == 0
@test RayTraceEllipsoids.absorption(1, 1) == 1/(1 + RayTraceEllipsoids.LOG10EXP)

#=s = Cylinder(1.,-2.,2.)
m = Retina(s, 0.1)
r = Ray(V3(-2tand(45),0,2), V3(sqrt(1/2),0,-sqrt(1/2)))
l = sqrt(4^2+4^2)
L = sqrt(2^2+2^2)
w = RayTraceEllipsoids.absorption(L, 0.1)
W = rand()
RayTraceEllipsoids.register!(m, r, l, W)
@test m.signal.photoreceptor == w
@test m.signal.retina == W

s = Cylinder(1.,-2.,2.)
m = Retina(s, 0.1)
r = Ray(V3(-2tand(45),0,2), V3(sqrt(1/2),0,-sqrt(1/2)))
l = sqrt(4^2+4^2)
o2 = r.orig + r.dir*l
RayTraceEllipsoids.absorbmove!(r, l, m)
i = RayTraceEllipsoids.absorption(l, 0.1)
@test r.int == i
@test r.orig == o2

s = Cylinder(1.,-2.,2.)
m = Retina(s, 1 - RayTraceEllipsoids.ENERGYTHRESHOLD)
r = Ray(V3(-2tand(45),0,2), V3(sqrt(1/2),0,-sqrt(1/2)))
l = sqrt(4^2+4^2)
o2 = r.orig + r.dir*l
@test_throws RayTraceEllipsoids.DeadRay() RayTraceEllipsoids.absorbmove!(r, l, m)=#

#=@test all(1:10^5) do _
    z = rand()
    c = V3(1.0,1,1)
    r = V3(2.0,2,2)
    orig = V3(c[1] + 1, c[2] + 1, c[3] + r[3] + z)
    dir = V3(0.0,0,-1)
    p1 = copy(dir)
    ray = Ray(orig, dir)
    h = true
    ss = Ellipsoid(c, r, h)
    l = RayTraceEllipsoids.distance(ray, ss)
    s = Cylinder(r[1], c[2] - r[2], c[3] + r[3])
    m = Retina(s, 0.5)
    RayTraceEllipsoids.absorbmove!(ray, l, m)
    ou = OpticUnit(ss, false, 0.9, m)
    RayTraceEllipsoids.bend!(ray, ou.interface)
    p2 = ray.dir
    acosd(dot(p1, p2) / sqrt(norm(p1)*norm(p2))) ≈ 45 - asind(sind(45)*0.9)
end=#
#=
@test all(1:10^5) do _
    d = 10*rand()
    r = rand()
    cz = rand()
    d += cz
    l = Light(d, 2r, r, r, cz, 0)
    ray = l(1)
    θ = atan(r/d)
    ray.dir ≈ V3(sin(θ), 0, -cos(θ)) && ray.orig == V3(0, 0, d + cz)
end

@test all(1:10^5) do _
    d = 10*rand()
    r = rand()
    cz = rand()
    d += cz
    l = Light(d, 2r, r, r, cz, 0)
    ray = l(1, 0.5)
    ray.dir ≈ V3(0, 0, -1) && ray.orig == V3(0, 0, d + cz)
end

@test all(1:10^5) do _
    z = rand()
    c = V3(1.0,1,1)
    r = V3(2.0,2,2)
    orig = V3(c[1] + 1, c[2] + 1, c[3] + r[3] + z)
    dir = V3(0.0,0,-1)
    p1 = copy(dir)
    ray = Ray(orig, dir)
    h = true
    ss = Ellipsoid(c, r, h)
    l = RayTraceEllipsoids.distance(ray, ss)
    s = Cylinder(r[1], c[2] - r[2], c[3] + r[3])
    m = Retina(s, 0.5)
    ou = OpticUnit(ss, false, 0.9, m)
    raytrace!(ray, ou)
    p2 = ray.dir
    acosd(dot(p1, p2) / sqrt(norm(p1)*norm(p2))) ≈ 45 - asind(sind(45)*0.9) && norm(ray.orig - orig) == l
end=#

orig = V3(0,0,1)
dir = V3(0,0,-1)
r = Ray(orig, dir)
s = Cylinder(1, -1, 1)
L = 2.0
@test RayTraceEllipsoids.incylinder(r, s, L) == L

orig = V3(1,2,1)
dir = V3(0,0,-1)
r = Ray(orig, dir)
s = Cylinder(1, -1, 1)
L = 2.0
@test RayTraceEllipsoids.incylinder(r, s, L) == 0

orig = V3(2,2,1)
dir = normalize(V3(0,1,-1))
r = Ray(orig, dir)
s = Cylinder(1, -1, 1)
L = 2.0
@test RayTraceEllipsoids.incylinder(r, s, L) == 0

orig = V3(1,1,1)
dir = normalize(V3(1,1,-1))
r = Ray(orig, dir)
s = Cylinder(1, -1, 1)
L = 2.0
@test RayTraceEllipsoids.incylinder(r, s, L) == 0

orig = V3(0,2,0)
dir = normalize(V3(0,1,0))
r = Ray(orig, dir)
s = Cylinder(1, -1, 1)
L = 2.0
@test RayTraceEllipsoids.incylinder(r, s, L) == 0

orig = V3(0,0,1)
dir = normalize(V3(0,-1,-1))
r = Ray(orig, dir)
s = Cylinder(1, -1, 1)
L = 2.0
@test RayTraceEllipsoids.incylinder(r, s, L) == sqrt(2)

orig = V3(0,0,1)
dir = normalize(V3(0,-0.1,-1))
r = Ray(orig, dir)
s = Cylinder(1, -1, 1)
L = 2.0
@test RayTraceEllipsoids.incylinder(r, s, L) == L

orig = V3(0,2,1)
dir = normalize(V3(0,-1,0))
r = Ray(orig, dir)
s = Cylinder(1, -1, 1)
L = 2.0
@test RayTraceEllipsoids.incylinder(r, s, L) == 1
L = 9.0
@test RayTraceEllipsoids.incylinder(r, s, L) == 2
L = 0.3
@test RayTraceEllipsoids.incylinder(r, s, L) == 0


text:
using Simplify: Symbolic, Match
using SpecialSets

@testset "Patterns" begin

    @syms f g h
    @syms a b c d e
    @vars x y z

    @testset "Symbolic" begin
        @test f == f
        @test f ≠ g
        @test f ≠ @term f
        @test f ≠ @term :f

        @test get(@term 1) == 1
        @test get(@term π) == pi
        @test get(@term im) == im
        @test get(@term sin) == sin

        @test_throws UndefVarError @term qwerty
    end

    @testset "Variable" begin
        @test x == x
        @test x == Variable(:x)
        @test x ≠ y

        @test match(@term(x), @term(y)) == Match(x => y)
        @test @term(y) ⊆ @term(x)

        @test replace(@term(x), Dict(x => y)) == @term(y)
        @test replace(@term(x), Dict(y => x)) == @term(x)
        @test replace(@term(x), Dict(y => z)) == @term(x)
    end

    @testset "Constant" begin
        _1, _2 = Term.([1, 2])
        @vars x y

        @test _1 == _1
        @test _1 ≠ _2

        @test _1 == @term 1

        @test match(_1, _1) == one(Match)
        @test _1 ⊆ _1
        @test match(_2, _1) == zero(Match)
        @test _1 ⊈ _2

        @test match(@term(:x), @term(:x)) == one(Match)
        @test match(@term(:x), @term(:y)) == zero(Match)

        @test match(@term(f(x, 0)), @term(f(y, 0))) == Match(x => y)

        @test replace(_1, Dict(x => y)) == _1
    end

    @testset "Function" begin

        @testset "standard" begin

            @test @term(f(x)) == @term(f(x))
            @test @term(f(x)) ≠ @term(f(y))

            @test match(@term(f()), @term(f())) == one(Match)
            @test match(@term(f(x)), @term(f(y))) == Match(x => y)
            @test @term(f(y)) ⊆ @term(f(x))
            @test match(@term(f(x)), @term(g(x))) == zero(Match)
            @test match(@term(f(f(), x)), @term(f(g(), y))) == zero(Match)
            @test match(@term(f(x, x)), @term(f(y, z))) == zero(Match)
            @test match(@term(f(x)), @term(g(x, y))) == zero(Match)
            @test match(@term(x - x), @term(-y)) == zero(Match)
            @test @term(g(x, y)) ⊈ @term(f(x))
            @test @term(f(:a, 2, :a)) ⊆ @term(f(x, 2, x))
            @test @term(f(:a, 2, :b)) ⊈ @term(f(x, 2, x))

            @test match(@term(x(y)), @term(sin(3))) == Match(Dict(x => sin, y => 3))
            @test match(@term(map(x, [])), @term(map(iseven, []))) == Match(x => iseven)
            @test match(@term(map(x, [])), @term(map(iseven, [1]))) == zero(Match)

            @test replace(@term(f(x)), Dict(x => y)) == @term(f(y))
            @test replace(@term(f(x)), Dict(y => x)) == @term(f(x))
        end

        @testset "associative" begin

            with_context(Context(Associative(f))) do
                @test match(@term(f(1, 2, f(3, 4))), @term(f(1, 2, 3, 4))) == zero(Match)

                @test match(@term(f(g(x), g(y), z)), @term(f(g(a), g(b), g(c), g(d), g(e)))) ==
                    Match(Dict(x => a, y => b, z => f(g(c), g(d), g(e))))

                @test match(@term(f(x, y)), @term(f(1, b))) ==
                    Match(Dict(x => 1, y => b))
            end

            with_context(Context(Associative(*))) do
                @test match(@term(x), @term(a * b)) ==
                    Match(x => get(@term a * b))

                @test match(@term(x * y), @term(a * b))::Match ==
                    Match(Dict(x => a, y => b))

                @test match(@term(x * y), @term(a * b * c)) == Match(
                    Dict(x => get(@term a * b), y => c),
                    Dict(x => a, y => get(@term b * c)),
                )

                @test match(@term(x * "_" * y), @term(a * b * "_" * c)) ==
                    Match(Dict(x => get(@term a * b), y => c))

                with_context(Context(Associative(*))) do
                    @test match(@term(a * b), @term(a * b)) == one(Match)
                    @test match(@term(a * b), @term(b * a)) == zero(Match)
                end

                @test replace(@term(a * b * (c * b)), Dict(b => d)) == @term(a * d * (c * d))
                @test_broken replace(@term(a * b * c), Dict(get(@term(b * c)) => 2)) == @term(a * 2)
            end

        end

        @testset "commutative" begin

            @testset "standard" begin

                with_context(Context(Commutative(f))) do
                    @test match(@term(f(x, 1)), @term(f(1, y))) == Match(x => y)

                    @test match(@term(f(g(x), g(y), z)), @term(f(h(a), g(b), g(c))))::Match == Match(
                        Dict(x => b, y => c, z => h(a)),
                        Dict(x => c, y => b, z => h(a)),
                    )

                    @test replace(@term(f(x, y)), Dict(x => 1)) == @term(f(1, y))
                    @test replace(@term(f(f(x, y), z)), Dict(f(x, y) => 2)) == @term(f(2, z))
                end

            end

            @testset "partial" begin
                @test_broken match(@term(f() * g()), @term(g() * f())) ==
                    one(Match)

                @test_broken match(@term(x * y), @term(1 * b)) == Match(
                    Dict(x => 1, y => b),
                    Dict(x => b, y => 1),
                )

                @test_broken normalize(@term(3 * (x * 2) * y)) == @term(2 * 3 * x * y)
            end

            @testset "associative" begin

                with_context(Context(Associative(f), Commutative(f))) do
                    @test replace(@term(f(x, y, z)), Dict(y => get(@term(x^3)))) == @term(f(x, x^3, z))
                    @test_skip replace(@term(f(x, y, z)), Dict(f(x, z) => 1)) == @term(f(1, y))
                end

                with_context(Context(Associative(+), Commutative(+), Associative(*), Commutative(*))) do
                    @test normalize(@term((x+y+b*a))) == normalize(@term((a*b+x+y)))
                end

                @test match(@term(x), @term(a + b)) ==
                    Match(x => get(@term(a + b)))

                @test match(@term(x + y), @term(a() + b()))::Match == Match(
                    Dict(x => a(), y => b()),
                    Dict(x => b(), y => a()),
                )

                @test match(@term(x + y), @term(1 + b)) == Match(
                    Dict(x => 1, y => b),
                    Dict(x => b, y => 1),
                )

                @test match(@term(x + y), @term(a + b)) == Match(
                    Dict(x => a, y => b),
                    Dict(x => b, y => a),
                )

                @test length(match(@term(x + y), @term(a + b + c))) == 12
                @test match(@term(a + b + c), @term(x + y)) == zero(Match)

                @test match(@term(x + 0), @term(f() + 0 + g())) == Match(
                    Dict(x => get(@term(f() + g()))),
                    Dict(x => get(@term(g() + f()))),
                )

                @test match(@term(x + y + 1), @term(a + 1 + b)) == Match(
                    Dict(x => a, y => b),
                    Dict(x => b, y => a),
                )

                @test match(@term(f() + g()), @term(f() + g())) ==
                    one(Match)

                @test match(@term(g() + f()), @term(f() + g())) ==
                    one(Match)

                @test match(@term(f() + f()), @term(f() + g())) ==
                    zero(Match)

            end

        end

    end

    @testset "show" begin
        let
            x = 1
            @test sprint(show, @term([x, :x])) == "@term([1, :x])"
        end
        let
            @syms x
            @test sprint(show, @term(2x + 3)) == "@term(2x + 3)"
        end
        @test sprint(show, @term(Base.Broadcast.materialize(1, 2))) ==
            "@term(Broadcast.materialize(1, 2))"
        let
            @syms a
            @vars x
            @test sprint(show, @term(2^3a+5x)) == "@term($(:(2^3a+5x)))"
        end
    end

    @testset "miscellaneous" begin
        @test match(@term([x, y, z]), @term([a, b, c])) == Match(Dict(
            x => a, y => b, z => c,
        ))
        @test match(@term([x, y, z]), @term([2a, a-1, b])) == Match(Dict(
            x => get(@term 2a), y => get(@term a-1), z => b
        ))
        @test match(@term([x, y, z]), @term([a, b])) == zero(Match)

        @test match(@term(x ? y : z), @term(a ? b : c)) == Match(Dict(
            x => a, y => b, z => c,
        ))
        @test match(@term(x ? y : y), @term(a ? b : b)) == Match(Dict(
            x => a, y => b,
        ))
        @test match(@term(x ? y : y), @term(a ? b : c)) == zero(Match)
    end

end


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 3, Zadanie 4.1.
#

include("iterative-methods.jl")
using .IterativeMethods

fl = IterativeMethods.FL

r, fr, k, err =  IterativeMethods.mbisekcji(
    x -> sin(x) - (x/2)^2,
    fl(1.5),
    fl(2),
    fl((10^(-5))/2),
    fl((10^(-5))/2)
)

println(
    "r = ", r, "\n",
    "f(r) = ", fr, "\n",
    "k = ", k, "\n",
    err == 0 ? "" : string("error code: ", err)
)

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook