text:
"""
    a, g, b, c = realizeNTF(ntf, form="CRFB", stf=nothing)

Convert a noise transfer function into coefficients for the desired structure.

Supported structures are:
- `CRFB`:    Cascade of resonators, feedback form.
- `CRFF`:    Cascade of resonators, feedforward form.
- `CIFB`:    Cascade of integrators, feedback form.
- `CIFF`:    Cascade of integrators, feedforward form.
- `CRFBD`:   CRFB with delaying quantizer.
- `CRFFD`:   CRFF with delaying quantizer.
- `PFF`:     Parallel feed-forward.
- `Stratos`: A CIFF-like structure with non-delaying resonator feedbacks,
             contributed in 2007 by Jeff Gealow
- `DSFB`:    Cascade of double-sampled integrators, feedback form.

The order of the `ntf` zeros must be (real, complex conj. pairs).
The order of the zeros is used when mapping the `ntf` onto the chosen topology.

`stf` is a TransferFunction object.
"""
function realizeNTF(ntf, form="CRFB", stf=nothing)
    # The basic idea is to equate the loop filter at a set of
    # points in the z-plane to L1 = 1-1/ntf at those points.

    # Code common to all functions
    ntf_z = ntf.matrix[1].z
    ntf_p = ntf.matrix[1].p
    order = length(ntf_p)
    order2 = order ÷ 2
    odd = order - 2*order2

    a = zeros(order)
    g = zeros(order2)
    b = zeros(order+1)
    c = ones(order)

    # Choose a set of points in the z-plane at which to try to make L1 = 1-1/H
    # Author of MATLAB code this is based on did't know how to do this properly,
    # but the code below seems to work for order <= 11
    N = 200
    min_distance = 0.09
    C = zeros(ComplexF64, N)
    j = 1
    for i = 1:N
        z = 1.1 * exp(2im*pi*i/N)
        if all(x -> abs(x - z) > min_distance, ntf_z)
            C[j] = z
            j = j + 1
        end
    end
    zSet = C[1:j-1]
    T = zeros(ComplexF64, order, order*2)

    if form =="CRFB"
        # Find g
        # Assume the roots are ordered, real first, then cx conj. pairs
        for i = 1:order2
            g[i] = 2 * (1 - real(ntf_z[2*i-1+odd]))
        end

        # Form the linear matrix equation a*T = L1
        L1 = zeros(ComplexF64, 1, order*2)
        for i = 1:order*2
            z = zSet[i]

            # L1(z) = 1-1/H(z)
            L1[i] = 1 - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z)

            Dfactor = (z - 1) / z
            product = 1
            for j = order:-2:(odd+1)
                product = z / evalRPoly(ntf_z[j-1:j], z) * product
                T[j,i] = Dfactor * product
                T[j-1,i] = product
            end
            if odd == 1
                T[1,i] = product / (z - 1)
            end
        end

        a = -real(L1 / T)[:]
        if isnothing(stf)
            b[1:order] = a
            b[order+1] = 1
        end

    elseif form == "CRFF"
        # Find g
        # Assume the roots are ordered, real first, then cx conj. pairs
        for i = 1:order2
            g[i] = 2 * (1 - real(ntf_z[2*i-1+odd]))
        end

        # Form the linear matrix equation a*T = L1
        L1 = zeros(ComplexF64, 1, order*2)
        for i = 1:order*2
            z = zSet[i]

            # L1(z) = 1-1/H(z)
            L1[i] = 1 - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z)

            if odd == 1
                Dfactor = z - 1
                product = 1 / Dfactor
                T[1,i] = product
            else
                Dfactor = (z - 1) / z
                product = 1
            end
            for j = (odd+1):2:order
                product = z / evalRPoly(ntf_z[j:j+1], z) * product
                T[j,i] = product * Dfactor
                T[j+1,i] = product
            end
        end

        a = -real(L1 / T)[:]
        if isnothing(stf)
            b = [1; zeros(order-1); 1]
        end

    elseif form == "CIFB"
        # Assume the roots are ordered, real first, then cx conj. pairs
        # Note ones which are moved significantly.
        if any(ntf_z -> abs(real(ntf_z) - 1) > 1e-3, ntf_z)
            @warn "The ntf's zeros have had their real parts set to one."
        end
        ntf_z = 1 .+ 1im*imag(ntf_z)

        for i = 1:order2
            g[i] = imag(ntf_z[2*i-1+odd])^2
        end

        # Form the linear matrix equation a*T = L1
        L1 = zeros(ComplexF64, 1, order*2)
        for i = 1:order*2
            z = zSet[i]

            # L1(z) = 1-1/H(z)
            L1[i] = 1 - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z)

            Dfactor = z - 1
            product = 1
            for j = order:-2:(odd+1)
                product = product / evalRPoly(ntf_z[j-1:j], z)
                T[j,i] = product * Dfactor
                T[j-1,i] = product
            end
            if odd == 1
                T[1,i] = product / (z - 1)
            end
        end

        a = -real(L1 / T)[:]
        if isnothing(stf)
            b[1:order] = a
            b[order+1] = 1
        end

    elseif form == "CIFF"
        # Assume the roots are ordered, real first, then cx conj. pairs
        # Note ones which are moved significantly.
        if any(ntf_z -> abs(real(ntf_z) - 1) > 1e-3, ntf_z)
            @warn "The ntf''s zeros have had their real parts set to one."
        end
        ntf_z = 1 .+ 1im*imag(ntf_z)

        for i = 1:order2
            g[i] = imag(ntf_z[2*i-1+odd])^2
        end

        # Form the linear matrix equation a*T = L1
        L1 = zeros(ComplexF64, 1, order*2)
        for i = 1:order*2
            z = zSet[i]

            # L1(z) = 1-1/H(z)
            L1[i] = 1 - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z)

            Dfactor = z - 1
            if odd == 1
                product = 1 / (z - 1)
                T[1,i] = product
            else
                product = 1
            end
            for  j= (odd+1):2:(order-1)
                product = product / evalRPoly(ntf_z[j:j+1], z)
                T[j,i] = product * Dfactor
                T[j+1,i] = product
            end
        end

        a = -real(L1 / T)[:]
        if isnothing(stf)
            b = [1; zeros(order-1); 1]
        end

    elseif form == "CRFBD"
        # Find g
        # Assume the roots are ordered, real first, then cx conj. pairs
        for i = 1:order2
            g[i] = 2 * (1 - real(ntf_z[2*i-1+odd]))
        end

        # Form the linear matrix equation a*T = L1
        L1 = zeros(ComplexF64, 1, order*2)
        for i = 1:order*2
            z = zSet[i]

            # L1(z) = 1-1/H(z)
            L1[i] = 1 - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z)

            Dfactor = z - 1
            product = 1 / z
            for j = order:-2:(odd+1)
                product = z / evalRPoly(ntf_z[(j-1):j], z) * product
                T[j,i] = product * Dfactor
                T[j-1,i] = product
            end
            if odd == 1
                T[1,i] = product * z / (z - 1)
            end
        end

        a = -real(L1 / T)[:]
        if isnothing(stf)
            b[1:order] = a
            b[order+1] = 1
        end

    elseif form == "CRFFD"
        # Find g
        # Assume the roots are ordered, real first, then cx conj. pairs
        for i = 1:order2
            g[i] = 2 * (1 - real(ntf_z[2*i-1+odd]))
        end

        # Form the linear matrix equation a*T = zL1
        zL1 = zeros(ComplexF64, 1, order*2)
        for i = 1:order*2
            z = zSet[i]

            # zL1 = z*(1-1/H(z))
            zL1[i] = z * (1 - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z))

            if odd == 1
                Dfactor = (z - 1) / z
                product = 1 / Dfactor
                T[1,i] = product
            else
                Dfactor = z - 1
                product = 1
            end
            for j = (odd+1):2:order
                product = z / evalRPoly(ntf_z[j:j+1], z) * product
                T[j,i] = product * Dfactor
                T[j+1,i] = product
            end
        end

        a = -real(zL1 / T)[:]
        if isnothing(stf)
            b = [1; zeros(order-1); 1]
        end

    elseif form == "PFF"
        # Find g
        # Assume the roots are ordered, real first, then cx conj. pairs
        # with the secondary zeros after the primary zeros
        for i = 1:order2
            g[i] = 2 * (1 - real(ntf_z[2*i-1+odd]))
        end

        # Find the dividing line between the zeros
        theta0 = abs(angle(ntf_z[1]))
        # 0.5 radians is an arbitrary separation
        i = findall(x -> abs(abs(angle(x))-theta0) > 0.5, ntf_z)
        order_1 = i[1] - 1
        order_2 = order - order_1
        if length(i) != order_2
            error("For the PFF form, the NTF zeros must be sorted into primary and secondary zeros")
        end
        odd_1 = mod(order_1, 2)
        odd_2 = mod(order_2, 2)

        # Form the linear matrix equation a*T = L1
        L1 = zeros(ComplexF64, 1, order*2)
        for i = 1:order*2
            z = zSet[i]

            # L(z) = 1-1/H(z)
            L1[i] = 1 - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z)

            if odd_1 == 1
                Dfactor = z - 1
                product = 1 / Dfactor
                T[1,i] = product
            else
                Dfactor = (z - 1) / z
                product = 1
            end
            for j = (odd_1+1):2:order_1
                product = z / evalRPoly(ntf_z[j:j+1], z) * product
                T[j,i] = product * Dfactor
                T[j+1,i] = product
            end

            if odd_2 == 1
                Dfactor = z - 1
                product = 1 / Dfactor
                T[order_1+1,i] = product
            else
                Dfactor = (z - 1) / z
                product = 1
            end
            for j = (order_1+1+odd_2):2:order
                product = z / evalRPoly(ntf_z[j:j+1], z) * product
                T[j,i] = product * Dfactor
                T[j+1,i] = product
            end
        end

        a = -real(L1 / T)[:]
        if isnothing(stf)
            b = [1; zeros(order_1-1); 1; zeros(order_2-1); 1]
        end

    elseif form == "Stratos"
        # code copied from case 'CRFF':
        # Find g
        # Assume the roots are ordered, real first, then cx conj. pairs
        for i = 1:order2
            g[i] = 2 * (1 - real(ntf_z[2*i-1+odd]))
        end

        # code copied from case 'CIFF':
        # Form the linear matrix equation a*T = L1
        L1 = zeros(ComplexF64, 1, order*2)
        for i = 1:order*2
            z = zSet[i]

            # L1(z) = 1-1/H(z)
            L1[i] = 1 - evalRPoly(ntf_p,z) / evalRPoly(ntf_z, z)

            Dfactor = z - 1
            if odd == 1
                product = 1 / (z - 1)
                T[1,i] = product
            else
                product = 1
            end
            for  j= (odd+1):2:(order-1)
                product = product / evalRPoly(ntf_z[j:j+1], z)
                T[j,i] = product * Dfactor
                T[j+1,i] = product
            end
        end

        a = -real(L1 / T)[:]
        if isnothing(stf)
            b = [1; zeros(order-1); 1]
        end

    elseif form == "DSFB"
        # Find g
        # Assume the roots are ordered, real first, then cx conj. pairs
        for i = 1:((order+odd)÷2-1)
            g[i] = tan(angle(ntf_z[2*i-1+odd]) / 2) ^ 2
        end
        if odd == 0
            g[order÷2] = 2 * (1 - cos(angle(ntf_z[order])))
        end

        # Form the linear matrix equation a*T = L1
        L1 = zeros(ComplexF64, 1, order*2)
        for i = 1:order*2
            z = zSet[i]

            # L(z) = 1-1/H(z)
            L1[i] = 1 - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z)

            if odd == 1
                product = 1 / (z - 1)
                T[order,i] = product
            else
                Dfactor = (z - 1)^2 + g[order÷2]*z
                T[order,i] = (z - 1) / Dfactor
                product = (z + 1) / Dfactor
                T[order-1,i] = product
            end
            for j = (order+odd-2):-2:1
                Dfactor = (z - 1)^2 + g[j÷2]*(z + 1)^2
                product = product / Dfactor
                T[j,i] = product * (z^2 - 1)
                product = product * (z + 1)^2
                T[j-1,i] = product
            end
        end

        a = -real(L1 / T)[:]
        if isnothing(stf)
            b[1:order] = a
            b[order+1] = 1
        end

    else
        error("Unsupported form $form")
    end

    if !isnothing(stf)
        # Compute the TF from each feed-in to the output
        # and solve for coefficients which yield the best match
        # THIS CODE IS NOT OPTIMAL, in terms of computational efficiency.
        stfList = Array{TransferFunction}(undef, order+1)
        for i = 1:order+1
            bi = zeros(order+1)
            bi[i] = 1
            ABCD = stuffABCD(a, g, bi, c, form)
            if occursin("FF", form)
                ABCD[1,order+2] = -1    # b1 is used to set the feedback
            end
            _, stfList[i] = calculateTF(ABCD)
        end
        # Build the matrix equation b*A = x and solve it.
        A = zeros(ComplexF64, order+1, length(zSet))
        for i = 1:order+1
            A[i,:] = evalTF(stfList[i], zSet)
        end
        x = evalTF(stf, zSet)
        x = transpose(x)
        b = real(x / A)[:]
    end

    return a, g, b, c
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
using DrWatson
@quickactivate "LightningArt"
using LightningArt
using Luxor
using LinearAlgebra
using ForwardDiff
import GeometryBasics
GB = GeometryBasics
using LightGraphs, MetaGraphs
using FileIO
using Revise
using ImageMagick
using ColorTypes, ColorVectorSpace
using Plots, GraphPlot
using Interpolations

guitar_img = load(datadir("exp_raw/guitar_model_1.png"))
plot(guitar_img)
bool_mask_to_blackwhite(bool) = bool .* RGBA(0.0, 0.0, 0.0, 1.0) .+ !bool .* RGBA(1.0, 1.0, 1.0, 1.0)

binary_collision_mask = map(col->norm(col) < 2, guitar_img)

function generate_edge_graph(binary_collision_mask)

    height, width = size(binary_collision_mask)
    pixel_graph = LightGraphs.SimpleGraph(height * width)
    linear_indices_mask = LinearIndices(binary_collision_mask)
    for j in 1:width-1, i in 1:height-1
        if !binary_collision_mask[i, j] && !binary_collision_mask[i+1, j]
            add_edge!(pixel_graph, linear_indices_mask[i, j], linear_indices_mask[i+1, j])
        end

        if !binary_collision_mask[i, j] && !binary_collision_mask[i, j+1]
            add_edge!(pixel_graph, linear_indices_mask[i, j], linear_indices_mask[i, j+1])
        end
    end
    pixel_graph
end

pixel_graph = generate_edge_graph(binary_collision_mask)

connected_components_guitar = LightGraphs.connected_components(pixel_graph)
large_connected_components_guitar = filter(x-> length(x) > 1, connected_components_guitar)
linear_indices_mask = LinearIndices(binary_collision_mask)
index_of_an_interior_pixel = linear_indices_mask[200, 100]
inside_connected_guitar_indices = Set(filter(x->index_of_an_interior_pixel in x, large_connected_components_guitar)[1])
length(inside_connected_guitar_indices)
inside_guitar_mask = map(linearindex->linearindex in inside_connected_guitar_indices, linear_indices_mask)
plot(bool_mask_to_blackwhite.(inside_guitar_mask))

inside_guitar_interpolation = interpolate(inside_guitar_mask, BSpline(Constant()))
function inside_guitar_function(p)
    interp_coords = clamp.(p, 0.0, 1.0) .* (size(inside_guitar_mask) .- 1) .+ 1
    inside_guitar_interpolation(interp_coords...) ≈ 1
end

height, width = size(inside_guitar_mask)
points_box = [GB.Point2(x, y) for x in 0.0:1/height:1.0, y in 0.0:1/width:1.0]
guitar_interp_mask = inside_guitar_function.(points_box)
plot(bool_mask_to_blackwhite.(guitar_interp_mask))

plot(bool_mask_to_blackwhite.(binary_collision_mask))

### TODO: define a function that is 1 inside the guitar and 0 outside

ext: jl

max_stars_repo_name: zoemcc/LightningArt

max_issues_repo_name: zoemcc/LightningArt

max_forks_repo_name: zoemcc/LightningArt


text:
module SumTypes

export @sum_type, @case #@match

using MacroTools: MacroTools

function match end
function parent end
constructors(::Type{T}) where T = throw(error("$T is not a sum type"))

"""
    @sum_type(T, blk)

Create a sum type `T` with constructors in the codeblock `blk`.

Examples:
*) An unparameterized sum type `Foo` with constructors `Bar` and `Baz`

    @sum_type Foo begin
        Bar(::Int)
        Baz(::Float64)
    end

    julia> Bar(1)
    Foo(Bar(1))

*) 'The' `Either` sum type with constructors `Left` and `Right`

    @sum_type Either{A, B} begin
        Left{A, B}(::A)
        Right{A, B}(::B)
    end

    julia> Left{Int, Int}(1)
    Either{Int64,Int64}(Left{Int64,Int64}(1))

    julia> Right{Int, Float64}(1.0)
    Either{Int64,Float64}(Right{Int64,Float64}(1.0))

*) A recursive `List` sum type with constructors `Nil` and `Cons`

    @sum_type List{A} begin
        Nil{A}()
        Cons{A}(::A, ::List{A})
    end

    julia> Nil{Int}()
    List{Int64}(Nil{Int64}())

    julia> Cons{Int}(1, Cons{Int}(1, Nil{Int}()))
    List{Int64}(Cons{Int64}(1, List{Int64}(Cons{Int64}(1, List{Int64}(Nil{Int64}())))))
"""
macro sum_type(T, blk::Expr, recur::Expr=:(recursive=false))
    @assert blk isa Expr && blk.head == :block
    T_name, T_params, T_params_constrained = if T isa Symbol
        T, [], []
    elseif T isa Expr && T.head == :curly
        T.args[1], (x -> x isa Expr && x.head == :(<:) ? x.args[1] : x).(T.args[2:end]), T.args[2:end]
    end
    filter!(x -> !(x isa LineNumberNode), blk.args)
    constructors = map(blk.args) do con_
        @assert isa(con_, Expr) "variants in sum type macro must be typed $(con_)(), not $(con_)"
        con::Expr = con_
        @assert con.head == :call
        con_name = con.args[1] isa Expr && con.args[1].head == :curly ? con.args[1].args[1] : con.args[1]
        con_params = (con.args[1] isa Expr && con.args[1].head == :curly) ? con.args[1].args[2:end] : []
        @assert con_params == T_params "constructors currently must have same parameters as the sum type. Got $T and $(con.args[1])"
        con_nameparam = isempty(con_params) ? con_name : :($con_name{$(con_params...)})
        con_field_names = map(enumerate(con.args[2:end])) do (i, field)
            @assert field isa Symbol || (field isa Expr && field.head == :(::)) "malformed constructor field $field"
            if field isa Symbol
                field
            elseif  length(field.args) == 1
                Symbol(:_, i)
            elseif length(field.args) == 2
                field.args[1]
            end
        end
        con_field_types = map(con.args[2:end]) do field
            @assert field isa Symbol || (field isa Expr && field.head == :(::)) "malformed constructor field $field"
            if field isa Symbol
                Any
            elseif  length(field.args) == 1
                field.args[1]
            elseif length(field.args) == 2
                field.args[2]
            end
        end
        (name=con_name, params=con_params, nameparam = con_nameparam, field_names=con_field_names, types=con_field_types)
    end
    out = Expr(:toplevel)
    foreach(constructors) do (name, params, nameparam, field_names, types)
        nameparam_constrained = isempty(params) ? name : :($name{$(T_params_constrained...)})
        field_names_typed = map(field_names, types) do name, type
            :($name :: $type)
        end
        struct_def = Expr(
            :struct, false, nameparam_constrained,
            Expr(:block,
                 field_names_typed...,
                 :($nameparam($(field_names_typed...)) where {$(T_params_constrained...)} =
                   $(Expr(:new, T, Expr(:new, nameparam, field_names...))))))
        ex = quote
            $struct_def
            @inline $Base.iterate(x::$name, s = 1) = s ≤ fieldcount($name) ? (getfield(x, s), s + 1) : nothing
            $Base.indexed_iterate(x::$name, i::Int, state=1) = (Base.@_inline_meta; (getfield(x, i), i+1))
            $SumTypes.parent(::Type{<:$name}) = $T_name
        end
        push!(out.args, ex)
    end

    con_nameparams = (x -> x.nameparam).(constructors)
    con_names      = (x -> x.name     ).(constructors)

    sum_struct_def = quote
        struct $T
            data::Union{$(con_nameparams...)}
            _1() = nothing
        end
    end
    #
    ex = quote
        $sum_struct_def
        $SumTypes.constructors(::Type{<:$T_name}) = ($(con_names...),)
        function $SumTypes.match(f, x::$T) where {$(T_params...)}
            _x = x.data
            $(_unionsplit(con_nameparams, :(f(_x))))
        end
    end
    push!(out.args, ex)
    MacroTools.@capture(recur, recursive=use_recur_) ||
        throw("Malformed recur option. Expected `recursive=true` or `recursive=false`")
    if use_recur
        pushfirst!(out.args, f(sum_struct_def)) # hack to allow mutually recursive types
    end
    esc(out)
end

f(x) = :(try $x catch; nothing end) #split out to it's own function to to possible parsinf error?

"""
    @case T fdef

Define a pattern matcher `fdef` to deconstruct a `SumType`

Examples:

    @case Either f((x,)::Left)  = x + 1
    @case Either f((x,)::Right) = x - 1

Calling `f` on an `Either` type will use manually unrolled dispatch, rather than julia's automatic dynamic dispatch machinery.That is, it'll emit code that is just a series of if/else calls.
"""
macro case(T, fdef)
    d = MacroTools.splitdef(fdef)
    f = esc(d[:name])
    T = esc(T)
    quote
        $f(x::$T)= $SumTypes.match($f, x)
        $(:($fdef) |> esc)
    end
end

function _unionsplit(thetypes, call)
    MacroTools.@capture(call, f_(arg_))
    first_type, rest_types = Iterators.peel(thetypes)
    code = :(if $arg isa $first_type; $call end)
    the_args = code.args
    for next_type in rest_types
        clause = :(if $arg isa $next_type # use `if` so this parses, then change to `elseif`
                   $call
                   end)
        clause.head = :elseif
        push!(the_args, clause)
        the_args = clause.args
    end
    push!(the_args, :(@assert false))
    return code
end

iscomplete(matcher, ::Type{T}) where {T} = all(constructors(T)) do con
    hasmethod(matcher, (con,))
end

macro foo()
    ex = quote end
    x = :(x = 1)
    pushfirst!(ex.args, :(try $x catch _; end))
    ex
end

end # module


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Grammar-Based Expression Search.
Available algorithms: Simulated Annealing (SA), Monte Carlo (MC), Grammatical Evolution (GE),
Monte Carlo Tree Search (MCTS) (no committing steps).

Usage: using ExprSearch.MC; result = exprsearch(p, problem)
"""
module ExprSearch

export ExprProblem, initialize!, get_grammar, get_fitness
export SearchParams, SearchResult, exprsearch
export get_expr, get_derivtree

const MODULEDIR = joinpath(dirname(@__FILE__), "..", "modules")

using Reexport
using RLESUtils, ModLoader
@reexport using GrammaticalEvolution

load_to_path(MODULEDIR)
const PKGS = readdir(MODULEDIR)

using DerivationTrees

import RLESTypes.SymbolTable
import DerivationTrees: initialize!, get_expr

abstract ExprProblem
abstract SearchParams
abstract SearchResult

exprsearch(p::SearchParams, problem::ExprProblem) = error("Please use a submodule.")

initialize!(problem::ExprProblem) = problem #user should overload this, default to nothing
get_grammar(problem::ExprProblem) = error("ExprSearch::get_grammar() not defined")
get_fitness(problem::ExprProblem, derivtree::DerivationTree,
    userargs::SymbolTable) = error("ExprSearch::get_fitness() not defined")

get_expr(::Void) = ""
get_fitness(::Void) = realmax(Float64)

get_derivtree(result::SearchResult) = error("ExprSearch::get_derivtree() not defined")

"""
Test an individual submodule
"""
function test(pkgs::AbstractString...; coverage::Bool=false)
  cd(() -> Pkg.Entry.test(AbstractString[pkgs...]; coverage=coverage), MODULEDIR)
end

"""
Test all submodules in modules folder.  Don't stop on error.
"""
function testall()
    for pkg in PKGS
        try
            test(pkg)
        catch
            println("Error in $pkg")
        end
    end
end

include("GP/src/GPExprSearch.jl") #Genetic Programming (standard tree-based)
include("GE/src/GEExprSearch.jl") #Grammatical Evolution
include("MC/src/MCExprSearch.jl") #Monte Carlo (random)
include("MCTS/src/MCTSExprSearch.jl") #Monte Carlo Tree Search without committing steps
include("CE/src/CEExprSearch.jl") #Cross-entropy method optimization using PCFG

#include("Ref/src/RefExprSearch.jl") #Ref
#include("PMC/src/PMCExprSearch.jl") #Parallel Monte Carlo
#include("SA/src/SAExprSearch.jl") #SA

end #module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
include("../../../src/histogram_code.jl")

if !@isdefined(CGS)
  include("../../../src/CGS.jl")
  using Main.CGS
end
using PyPlot
using Printf
using JLD2

# Load in the likelihood profile dataset:
@load "../../../data/T1_likelihood_profile_student_all_3.0sig.jld2"

# Load in the Markov chain dataset:
@load "../../../data/T1_hmc_total_02212020.jld2"

using Optim

function normal_data(m0,sm0,mgrid,dchi)
  function fit_normal(x)
    return sum((exp.(-0.5*dchi) .-exp.(-0.5*((mgrid .-x[1])/x[2]).^2)).^2)
  end
  result = Optim.optimize(fit_normal,[m0,sm0])
  return Optim.minimizer(result)
end

;clf()
fig,axes = subplots(3,1,figsize=(8,8))
fac = 0.09*MSUN/MEARTH; dchi = chi_grid_all .-minimum(chi_grid_all)
# Overax.plot Simon's values:
#mgrimm = [1.310,1.307,0.367,0.642,0.950,1.245,0.373]
mgrimm = [1.017,1.156,0.297,0.772,0.934,1.148,0.331]
#sgrimm = [0.082,0.073,0.022,0.027,0.026,0.027,0.034]
sgrimm1 = [0.143,0.131,0.035,0.075,0.078,0.095,0.049]
sgrimm2 = [0.154,0.142,0.039,0.079,0.080,0.098,0.056]
#cplanet = [:blue,:green,cplanet[i],cplanet[i],cplanet[i],:yellow,:black]
cplanet = ["C0","C1","C2","C3","C4","C5","C6"]
planet =["b","c","d","e","f","g","h"]
# Fit Gaussians to each set of parameters:
nplanet = n-1
xfit = zeros(2,nplanet)
for i=1:nplanet
  xfit[:,i] = normal_data(elements_grid_all[5i-4,i+1,1,16]*fac,0.01,elements_grid_all[5i-4,i+1,1,:]*fac,dchi[5i-4,:])
  println("planet: ",i," mass: ",@sprintf("%6.4f",xfit[1,i]),"+-",@sprintf("%6.4f",xfit[2,i]))
end
#i = 5
#xfit[:,i] = normal_data(1.0,sgrimm1[i],elements_grid_all[5i-4,i+1,1,:]*fac,dchi[5i-4,:])

# First, ax.plot planets d, h:
ax = axes[1]
i=3
mgrid = elements_grid_all[5i-4,i+1,1,:].*fac
ax.plot(mgrid,exp.(-.5*dchi[5i-4,:]),label=planet[i],linewidth=1,color=cplanet[i])
#ax.plot(mgrid,exp.(-.5*(mgrid.-xfit[1,i]).^2/xfit[2,i]^2),linestyle=":",color=:black)
ax.plot(mgrid,exp.(-.5*(mgrid.-xfit[1,i]).^2/(cov_save[(i-1)*5+1,(i-1)*5+1]*fac^2)),linestyle=":",color=cplanet[i],linewidth=2,alpha=0.5)
ax.plot(mgrimm[i]*[1,1],[0.4,0.7],color=cplanet[i])
ax.plot((mgrimm[i]-sgrimm1[i])*[1,1],[0.4,0.7],color=cplanet[i],linewidth=2)
ax.plot((mgrimm[i]+sgrimm2[i])*[1,1],[0.4,0.7],color=cplanet[i],linewidth=2)
ax.plot(mgrimm[i] .+[-sgrimm1[i],sgrimm2[i]],[0.55,0.55],color=cplanet[i],linewidth=2)
#ax.errorbar(mgrimm[i],0.5,xerr=[sgrimm1[i],sgrimm2[i]],color=cplanet[i])
# Now, make a histogram:
md_bin,md_hist,md_bin_square,md_hist_square = histogram(state_total[(i-1)*5+1,:].*fac,50)
#ax.plot(md_bin,md_hist ./maximum(md_hist),color=cplanet[i],linewidth=2,linestyle="--")
ax.plot(md_bin_square,md_hist_square ./maximum(md_hist_square),color=cplanet[i],linewidth=2)
i=7
mgrid = elements_grid_all[5i-4,i+1,1,:].*fac
ax.plot(mgrid,exp.(-.5*dchi[5i-4,:]),label=planet[i],linewidth=1,color=cplanet[i])
#ax.plot(mgrid,exp.(-.5*(mgrid .-xfit[1,i]).^2/xfit[2,i]^2),linestyle=":",color=:black)
ax.plot(mgrid,exp.(-.5*(mgrid .-xfit[1,i]).^2/(cov_save[(i-1)*5+1,(i-1)*5+1]*fac^2)),linestyle=":",color=cplanet[i],linewidth=2,alpha=0.5)
ax.plot(mgrimm[i]*[1,1],[0.35,0.65],color=cplanet[i])
ax.plot((mgrimm[i]-sgrimm1[i])*[1,1],[0.35,0.65],color=cplanet[i],linewidth=2)
ax.plot((mgrimm[i]+sgrimm2[i])*[1,1],[0.35,0.65],color=cplanet[i],linewidth=2)
ax.plot(mgrimm[i] .+[-sgrimm1[i],sgrimm2[i]],[0.5,0.5],color=cplanet[i],linewidth=2)
# Now plot HMC posterior:
mh_bin,mh_hist,mh_bin_square,mh_hist_square = histogram(state_total[(i-1)*5+1,:].*fac,50)
#ax.plot(mh_bin,mh_hist./maximum(mh_hist),color=cplanet[i],linewidth=2,linestyle="--")
ax.plot(mh_bin_square,mh_hist_square./maximum(mh_hist_square),color=cplanet[i],linewidth=2)
ax.set_xlabel(L"$(M/M_\oplus)(M_*/0.09 M_\odot)^{-1}$",fontsize=10)
#ax.set_xlabel(L"$\frac{M}{M_\oplus}\frac{0.09 M_\odot}{M_*}$",fontsize=10)
#ax.set_ylabel(L"$\exp(-\frac{1}{2}\Delta \chi^2)$",fontsize=10)
ax.set_ylabel("Probability",fontsize=10)
ax.legend(loc = "upper left")
ax.axis([0.25,0.45,0,1.1])
#ax.plot([0.2,0.5],exp(-0.5)*[1,1],color="k",linestyle=":")
# Next, e, f:
ax = axes[2]
i=4
mgrid = elements_grid_all[5i-4,i+1,1,:].*fac
ax.plot(mgrid,exp.(-.5*dchi[5i-4,:]),label=planet[i],linewidth=1,color=cplanet[i])
#ax.plot(mgrid,exp.(-.5*(mgrid.-xfit[1,i]).^2/xfit[2,i]^2),linestyle=":",color=cplanet[i],linewidth=2)
ax.plot(mgrid,exp.(-.5*(mgrid.-xfit[1,i]).^2/(cov_save[(i-1)*5+1,(i-1)*5+1]*fac^2)),linestyle=":",color=cplanet[i],alpha=0.5,linewidth=2)
ax.plot(mgrimm[i]*[1,1],[0.45,0.75],color=cplanet[i])
ax.plot((mgrimm[i]-sgrimm1[i])*[1,1],[0.45,0.75],color=cplanet[i],linewidth=2)
ax.plot((mgrimm[i]+sgrimm2[i])*[1,1],[0.45,0.75],color=cplanet[i],linewidth=2)
ax.plot(mgrimm[i] .+[-sgrimm1[i],sgrimm2[i]],[0.60,0.60],color=cplanet[i],linewidth=2)
# Now plot HMC posterior:
me_bin,me_hist,me_bin_square,me_hist_square = histogram(state_total[(i-1)*5+1,:].*fac,50)
#ax.plot(me_bin,me_hist./maximum(me_hist),color=cplanet[i],linewidth=2,linestyle="--")
ax.plot(me_bin_square,me_hist_square./maximum(me_hist_square),color=cplanet[i],linewidth=2)
i=5
mgrid = elements_grid_all[5i-4,i+1,1,:].*fac
ax.plot(mgrid,exp.(-.5*dchi[5i-4,:]),label=planet[i],linewidth=1,color=cplanet[i])
#ax.plot(mgrid,exp.(-.5*(mgrid .-xfit[1,i]).^2/xfit[2,i]^2),linestyle=":",color=:black)
ax.plot(mgrid,exp.(-.5*(mgrid .-xfit[1,i]).^2/(cov_save[(i-1)*5+1,(i-1)*5+1]*fac^2)),linestyle=":",color=cplanet[i],alpha=0.5,linewidth=2)
# Now plot HMC posterior:
mf_bin,mf_hist,mf_bin_square,mf_hist_square = histogram(state_total[(i-1)*5+1,:].*fac,50)
ax.plot(mf_bin_square,mf_hist_square./maximum(mf_hist_square),color=cplanet[i],linewidth=2)
#ax.plot(mf_bin,mf_hist./maximum(mf_hist),color=cplanet[i],linewidth=2,linestyle="--")
ax.plot(mgrimm[i]*[1,1],[0.4,0.7],color=cplanet[i])
ax.plot((mgrimm[i]-sgrimm1[i])*[1,1],[0.4,0.7],color=cplanet[i],linewidth=2)
ax.plot((mgrimm[i]+sgrimm2[i])*[1,1],[0.4,0.7],color=cplanet[i],linewidth=2)
ax.plot(mgrimm[i] .+[-sgrimm1[i],sgrimm2[i]],[0.55,0.55],color=cplanet[i],linewidth=2)
ax.set_xlabel(L"$(M/M_\oplus)(M_*/0.09 M_\odot)^{-1}$",fontsize=10)
#ax.set_ylabel(L"$\exp(-\frac{1}{2}\Delta \chi^2)$",fontsize=10)
ax.set_ylabel("Probability",fontsize=10)
ax.legend(loc = "upper left")
ax.axis([0.6,1.1,0,1.1])
#ax.plot([0.5,1.1],exp(-0.5)*[1,1],color="k",linestyle=":")
# Extend errors from b/c/g to this panel:
i=1
ax.plot(mgrimm[i]*[1,1],[0.45,0.75],color=cplanet[i])
ax.plot((mgrimm[i]-sgrimm1[i])*[1,1],[0.45,0.75],color=cplanet[i],linewidth=2)
ax.plot((mgrimm[i]+sgrimm2[i])*[1,1],[0.45,0.75],color=cplanet[i],linewidth=2)
ax.plot(mgrimm[i] .+[-sgrimm1[i],sgrimm2[i]],[0.6,0.6],color=cplanet[i],linewidth=2)
i=2
ax.plot(mgrimm[i]*[1,1],[0.4,0.7],color=cplanet[i])
ax.plot((mgrimm[i]-sgrimm1[i])*[1,1],[0.4,0.7],color=cplanet[i],linewidth=2)
ax.plot((mgrimm[i]+sgrimm2[i])*[1,1],[0.4,0.7],color=cplanet[i],linewidth=2)
ax.plot(mgrimm[i] .+[-sgrimm1[i],sgrimm2[i]],[0.55,0.55],color=cplanet[i],linewidth=2)
i=6
ax.plot(mgrimm[i]*[1,1],[0.35,0.65],color=cplanet[i])
ax.plot((mgrimm[i]-sgrimm1[i])*[1,1],[0.35,0.65],color=cplanet[i],linewidth=2)
ax.plot((mgrimm[i]+sgrimm2[i])*[1,1],[0.35,0.65],color=cplanet[i],linewidth=2)
ax.plot(mgrimm[i] .+[-sgrimm1[i],sgrimm2[i]],[0.5,0.5],color=cplanet[i],linewidth=2)
# Next, b, c, g:
ax = axes[3]
i=1
mgrid = elements_grid_all[5i-4,i+1,1,:].*fac
ax.plot(mgrid,exp.(-.5*dchi[5i-4,:]),label=planet[i],linewidth=1,color=cplanet[i])
#ax.plot(mgrid,exp.(-.5*(mgrid.-xfit[1,i]).^2/xfit[2,i]^2),linestyle=":",color=:black)
ax.plot(mgrid,exp.(-.5*(mgrid.-xfit[1,i]).^2/(cov_save[(i-1)*5+1,(i-1)*5+1]*fac^2)),linestyle=":",color=cplanet[i],alpha=0.5,linewidth=2)
ax.plot(mgrimm[i]*[1,1],[0.45,0.75],color=cplanet[i])
ax.plot((mgrimm[i]-sgrimm1[i])*[1,1],[0.45,0.75],color=cplanet[i],linewidth=2)
ax.plot((mgrimm[i]+sgrimm2[i])*[1,1],[0.45,0.75],color=cplanet[i],linewidth=2)
ax.plot(mgrimm[i] .+[-sgrimm1[i],sgrimm2[i]],[0.6,0.6],color=cplanet[i],linewidth=2)
# Now plot HMC posterior:
mb_bin,mb_hist,mb_bin_square,mb_hist_square = histogram(state_total[(i-1)*5+1,:].*fac,50)
#ax.plot(mb_bin,mb_hist./maximum(mb_hist),color=cplanet[i],linewidth=2,linestyle="--")
ax.plot(mb_bin_square,mb_hist_square./maximum(mb_hist_square),color=cplanet[i],linewidth=2)
i=2
mgrid = elements_grid_all[5i-4,i+1,1,:].*fac
ax.plot(mgrid,exp.(-.5*dchi[5i-4,:]),label=planet[i],linewidth=1,color=cplanet[i])
#ax.plot(mgrid,exp.(-.5*(mgrid.-xfit[1,i]).^2/xfit[2,i]^2),linestyle=":",color=:black)
ax.plot(mgrid,exp.(-.5*(mgrid.-xfit[1,i]).^2/(cov_save[(i-1)*5+1,(i-1)*5+1]*fac^2)),linestyle=":",color=cplanet[i],alpha=0.5,linewidth=2)
# Now plot HMC posterior:
mc_bin,mc_hist,mc_bin_square,mc_hist_square = histogram(state_total[(i-1)*5+1,:].*fac,50)
#ax.plot(mc_bin,mc_hist./maximum(mc_hist),color=cplanet[i],linewidth=2,linestyle="--")
ax.plot(mc_bin_square,mc_hist_square./maximum(mc_hist_square),color=cplanet[i],linewidth=2)
ax.plot(mgrimm[i]*[1,1],[0.4,0.7],color=cplanet[i])
ax.plot((mgrimm[i]-sgrimm1[i])*[1,1],[0.4,0.7],color=cplanet[i],linewidth=2)
ax.plot((mgrimm[i]+sgrimm2[i])*[1,1],[0.4,0.7],color=cplanet[i],linewidth=2)
ax.plot(mgrimm[i] .+[-sgrimm1[i],sgrimm2[i]],[0.55,0.55],color=cplanet[i],linewidth=2)
i=6
mgrid = elements_grid_all[5i-4,i+1,1,:].*fac
ax.plot(mgrid,exp.(-.5*dchi[5i-4,:]),label=planet[i],linewidth=1,color=cplanet[i])
#ax.plot(mgrid,exp.(-.5*(mgrid.-xfit[1,i]).^2/xfit[2,i]^2),linestyle=":",color=:black)
ax.plot(mgrid,exp.(-.5*(mgrid.-xfit[1,i]).^2/(cov_save[(i-1)*5+1,(i-1)*5+1]*fac^2)),linestyle=":",color=cplanet[i],alpha=0.5,linewidth=2)
# Now plot HMC posterior:
mg_bin,mg_hist,mg_bin_square,mg_hist_square = histogram(state_total[(i-1)*5+1,:].*fac,50)
#ax.plot(mg_bin,mg_hist./maximum(mg_hist),color=cplanet[i],linewidth=2,linestyle="--")
ax.plot(mg_bin_square,mg_hist_square./maximum(mg_hist_square),color=cplanet[i],linewidth=2)
ax.plot(mgrimm[i]*[1,1],[0.35,0.65],color=cplanet[i])
ax.plot((mgrimm[i]-sgrimm1[i])*[1,1],[0.35,0.65],color=cplanet[i],linewidth=2)
ax.plot((mgrimm[i]+sgrimm2[i])*[1,1],[0.35,0.65],color=cplanet[i],linewidth=2)
ax.plot(mgrimm[i] .+[-sgrimm1[i],sgrimm2[i]],[0.5,0.5],color=cplanet[i],linewidth=2)
ax.set_xlabel(L"$(M/M_\oplus)(M_*/0.09 M_\odot)^{-1}$",fontsize=10)
#ax.set_ylabel(L"$\exp(-\frac{1}{2}\Delta \chi^2)$",fontsize=10)
ax.set_ylabel("Probability",fontsize=10)
ax.legend()
ax.axis([1.1,1.6,0,1.1])
#ax.plot([1.1,1.6],exp(-0.5)*[1,1],color="k",linestyle=":")

tight_layout()
savefig("../T1_masses_03312020.pdf", bbox_inches="tight")


text:
function show(io::IO, p::Poly)
    print(io,"Poly(")
    print(io,p)
    print(io,")")
end

function printexponent(io,var,i)
    if i == 0
    elseif i == 1
        print(io,var)
    else
        print(io,var,"^",i)
    end
end

function printterm{T}(io::IO,p::Poly{T},j,first)
    pj = p[j]
    if pj == zero(T)
        return false
    end
    neg = pj < zero(T)
    if first
        neg && print(io, "-")    #Prepend - if first and negative
    else
        neg ? print(io, " - ") : print(io," + ")
    end
    pj = abs(pj)
    if pj != one(T) || j == 0
        show(io,pj)
    end
    printexponent(io,p.var,j)
    true
end

function printterm{T<:Complex}(io::IO,p::Poly{T},j,first)
    pj = p[j]
    abs_repj = abs(real(pj))
    abs_impj = abs(imag(pj))
    if abs_repj < 2*eps(T) && abs_impj < 2*eps(T)
        return false
    end

    # We show a negative sign either for any complex number with negative
    # real part (and then negate the immaginary part) of for complex
    # numbers that are pure imaginary with negative imaginary part

    neg = ((abs_repj > 2*eps(T)) && real(pj) < 0) ||
            ((abs_impj > 2*eps(T)) && imag(pj) < 0)

    if first
        neg && print(io, "-")    #Prepend - if first and negative
    else
        neg ? print(io," - ") : print(io," + ")
    end

    if abs_repj > 2*eps(T)    #Real part is not 0
        if abs_impj > 2*eps(T)    #Imag part is not 0
            print(io,'(',neg ? -pj : pj,')')
        else
            print(io, neg ? -real(pj) : real(pj))
        end
    else
        if abs_impj > 2*eps(T)
            print(io,'(', imag(pj),"im)")
        end
    end
    printexponent(io,p.var,j)
    true
end

function print{T}(io::IO, p::Poly{T})
    first = true
    printed_anything = false
    n = length(p)-1
    for i = 0:n
        printed = printterm(io,p,i,first)
        first &= !printed
        printed_anything |= printed
    end
    printed_anything || print(io,zero(T))
end


text:
function _optimisedTauLeap!(model::Model, p::Path, tf::Float32, tauOpt::TauOpt, tauVar::TauVar)

    # ---- Main loop
    a = zeros(Float64, model.nr)
    while p.t < tf
        moleculeLimit!(p, tf) && break

        model.F(p.x, p.θ, a, p.t)
        if sum(a) <= 0.0
            break
        end
        tr = tf - p.t

        τ, switch = computeLeap(model, p, a, tr, tauOpt, tauVar)

        # ---- Switch to direct method?
        if switch == true
            _nDirectMethod!(p, model.F, model.ν, tf, model.nr, tauOpt.nd)

        # ---- Else accept step
        else
            for i in eachindex(p.x)
                p.x[i] += tauVar.dx[i]
            end

            p.t += τ
        end
    end
end

_optimisedTauLeap = _optimisedTauLeap!

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
module iss  #jl

using ReachabilityAnalysis, MAT
using ReachabilityModels: @relpath

file = matopen(@relpath "iss.mat")

# system matrix
A = read(file, "A")

# input matrix
B = read(file, "B")

# state domain
X = Universe(270)

# input domain
U = Hyperrectangle([0.05, 0.9, 0.95], [0.05, 0.1, 0.05])

function model(X0)
    S = @system(x' = Ax + Bu, x ∈ X, u ∈ U)
    return IVP(S, X0)
end

end  #jl

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using TerminalRegressionTests
using Test

const thisdir = dirname(@__FILE__)
TerminalRegressionTests.automated_test(
                joinpath(thisdir,"TRT.multiout"),
                ["Julia\n","Yes!!\n"]) do emuterm
    print(emuterm, "Please enter your name: ")
    name = strip(readline(emuterm))
    print(emuterm, "\nHello $name. Do you like tests? ")
    resp = strip(readline(emuterm))
    @assert resp == "Yes!!"
end

mktemp() do _, io
    redirect_stderr(io) do
        redirect_stdout(io) do
            @test_throws ErrorException TerminalRegressionTests.automated_test(
                            joinpath(thisdir,"TRT2.multiout"),
                            [""]) do emuterm
                println(emuterm, "Hello, world!")   # generate with "wurld" rather than "world"
                readline(emuterm)   # needed to produce output?
            end
        end
    end
end

function compare_replace(em, output; replace=nothing)
    buf = IOBuffer()
    decoratorbuf = IOBuffer()
    TerminalRegressionTests.VT100.dump(buf,decoratorbuf,em)
    outbuf = take!(buf)
    if replace !== nothing
        output = Base.replace(output, replace)
    end
    TerminalRegressionTests._compare(Vector{UInt8}(codeunits(output)), outbuf) || return false
    return true
end
const cmp(a, b, decorator) = compare_replace(a, b; replace="wurld"=>"world")
TerminalRegressionTests.automated_test(cmp,
                joinpath(thisdir,"TRT2.multiout"),
                [""]) do emuterm
    println(emuterm, "Hello, world!")
    readline(emuterm)
end


text:
module Curl

export
    with_handle,
    Easy,
        set_url,
        set_method,
        set_verbose,
        set_body,
        set_upload_size,
        set_seeker,
        set_ca_roots_path,
        set_timeout,
        add_headers,
        enable_upload,
        enable_progress,
        upload_data,
        get_protocol,
        get_effective_url,
        get_response_status,
        get_response_info,
        get_curl_errstr,
        status_ok,
    Multi,
        add_handle,
        remove_handle

using LibCURL
using LibCURL: curl_off_t
# not exported: https://github.com/JuliaWeb/LibCURL.jl/issues/87

using NetworkOptions
using Base: preserve_handle, unpreserve_handle

include("utils.jl")

function __init__()
    @check curl_global_init(CURL_GLOBAL_ALL)
end

const CURL_VERSION_INFO = unsafe_load(curl_version_info(CURLVERSION_NOW))
const SSL_VERSION = unsafe_string(CURL_VERSION_INFO.ssl_version)
const SYSTEM_SSL =
    Sys.isapple() && startswith(SSL_VERSION, "SecureTranspart")
    Sys.iswindows() && startswith(SSL_VERSION, "Schannel")

const CURL_VERSION = unsafe_string(curl_version())
const USER_AGENT = "$CURL_VERSION julia/$VERSION"

include("Easy.jl")
include("Multi.jl")

function with_handle(f, handle::Union{Multi, Easy})
    try f(handle)
    finally
        Curl.done!(handle)
    end
end

end # module


text:
using CompScienceMeshes

m = CompScienceMeshes.tetmeshsphere(1.0, 0.35)
@show numcells(m)

br = barycentric_refinement(m)
@show numcells(br)

@assert numcells(br) == 24*numcells(m)

V = 155
sm = submesh(tet -> V in tet, br.mesh)
@show numcells(sm)

bndry = boundary(sm)
import PlotlyJS
PlotlyJS.plot(patch(bndry))

E = 100
Edges = skeleton(m,1)
sm = submesh(br.mesh) do tet
    cells(Edges)[E][1] in tet && return true
    cells(Edges)[E][2] in tet && return true
    return false
end

PlotlyJS.plot(patch(boundary(sm)))


text:
# Autogenerated wrapper script for ONNXRuntime_jll for x86_64-w64-mingw32
export libonnxruntime

JLLWrappers.@generate_wrapper_header("ONNXRuntime")
JLLWrappers.@declare_library_product(libonnxruntime, "onnxruntime.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libonnxruntime,
        "bin\\onnxruntime.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Feather
include(joinpath("src", "HNScraper.jl"))
data = HNScraper.scrape_topstories(n_timeslices = 60 * 48, interval_minutes = 1)
Feather.write(joinpath("data", "hacker_news_data.feather"), data)

newstories = HNScraper.scrape_newstories(n_timeslices = 60, interval_minutes = 1)
Feather.write(joinpath("data", "hacker_news_newstories.feather"), newstories)

max_stars_repo_name: JohannesNakayama/hacker-news-analysis

max_issues_repo_name: JohannesNakayama/hacker-news-analysis

max_forks_repo_name: JohannesNakayama/hacker-news-analysis


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

using Markdown
using InteractiveUtils

# ╔═╡ 2c0f3d74-8b1e-4648-b0fa-f8050035b0fd
using Pkg, InteractiveUtils

# ╔═╡ ac57d957-0bdd-457a-ac15-9a4f94f0c785
# Remove this cell to use release versions of dependencies
# hideall
let
    docs_dir = dirname(@__DIR__)
    pkg_dir = dirname(docs_dir)

    Pkg.activate(docs_dir)
    Pkg.develop(; path=pkg_dir)
    Pkg.instantiate()
end;

# ╔═╡ 467c2d13-6bfe-4feb-9626-fb14796168aa
begin
    using ArviZ, CmdStan, Distributions, LinearAlgebra, PyPlot, Random, Soss, Turing
    using Soss.MeasureTheory: HalfCauchy
    using SampleChainsDynamicHMC: getchains, dynamichmc
end

# ╔═╡ 56a39a90-0594-48f4-ba04-f7b612019cd1
using PlutoUI

# ╔═╡ a23dfd65-50a8-4872-8c41-661e96585aca
md"""
# [ArviZ.jl Quickstart](#quickstart)

!!! note

    This tutorial is adapted from [ArviZ's quickstart](https://python.arviz.org/en/latest/getting_started/Introduction.html).
"""

# ╔═╡ d2eedd48-48c6-4fcd-b179-6be7fe68d3d6
md"""
## [Setup](#setup)

Here we add the necessary packages for this notebook and load a few we will use throughout.
"""

# ╔═╡ 06b00794-e97f-472b-b526-efe4815103f8
# ArviZ ships with style sheets!
ArviZ.use_style("arviz-darkgrid")

# ╔═╡ 5acbfd9a-cfea-4c3c-a3f0-1744eb7e4e27
md"""
## [Get started with plotting](#Get-started-with-plotting)

ArviZ.jl is designed to be used with libraries like [CmdStan](https://github.com/StanJulia/CmdStan.jl), [Turing.jl](https://turing.ml), and [Soss.jl](https://github.com/cscherrer/Soss.jl) but works fine with raw arrays.
"""

# ╔═╡ efb3f0af-9fac-48d8-bbb2-2dd6ebd5e4f6
rng1 = Random.MersenneTwister(37772);

# ╔═╡ 401e9b91-0bca-4369-8d36-3d9f0b3ad60b
begin
    plot_posterior(randn(rng1, 100_000))
    gcf()
end

# ╔═╡ 2c718ea5-2800-4df6-b62c-e0a9e440a1c3
md"""
Plotting a dictionary of arrays, ArviZ.jl will interpret each key as the name of a different random variable.
Each row of an array is treated as an independent series of draws from the variable, called a _chain_.
Below, we have 10 chains of 50 draws each for four different distributions.
"""

# ╔═╡ 49f19c17-ac1d-46b5-a655-4376b7713244
let
    s = (10, 50)
    plot_forest(
        Dict(
            "normal" => randn(rng1, s),
            "gumbel" => rand(rng1, Gumbel(), s),
            "student t" => rand(rng1, TDist(6), s),
            "exponential" => rand(rng1, Exponential(), s),
        ),
    )
    gcf()
end

# ╔═╡ a9789109-2b90-40f7-926c-e7c87025d15f
md"""
## [Plotting with MCMCChains.jl's `Chains` objects produced by Turing.jl](#Plotting-with-MCMCChains.jl's-Chains-objects-produced-by-Turing.jl)

ArviZ is designed to work well with high dimensional, labelled data.
Consider the [eight schools model](https://statmodeling.stat.columbia.edu/2014/01/21/everything-need-know-bayesian-statistics-learned-eight-schools/), which roughly tries to measure the effectiveness of SAT classes at eight different schools.
To show off ArviZ's labelling, I give the schools the names of [a different eight schools](https://en.wikipedia.org/wiki/Eight_Schools_Association).

This model is small enough to write down, is hierarchical, and uses labelling.
Additionally, a centered parameterization causes [divergences](https://mc-stan.org/users/documentation/case-studies/divergences_and_bias.html) (which are interesting for illustration).

First we create our data and set some sampling parameters.
"""

# ╔═╡ 69124a94-a6aa-43f6-8d4f-fa9a6b1cfef1
begin
    J = 8
    y = [28.0, 8.0, -3.0, 7.0, -1.0, 1.0, 18.0, 12.0]
    σ = [15.0, 10.0, 16.0, 11.0, 9.0, 11.0, 10.0, 18.0]
    schools = [
        "Choate",
        "Deerfield",
        "Phillips Andover",
        "Phillips Exeter",
        "Hotchkiss",
        "Lawrenceville",
        "St. Paul's",
        "Mt. Hermon",
    ]
    ndraws = 1_000
    ndraws_warmup = 1_000
    nchains = 4
end;

# ╔═╡ 10986e0e-55f7-40ea-ab63-274fbab3126d
md"""
Now we write and run the model using Turing:
"""

# ╔═╡ f383d541-e22d-44b4-b8cb-28b3d67944a1
Turing.@model function model_turing(y, σ, J=length(y))
    μ ~ Normal(0, 5)
    τ ~ truncated(Cauchy(0, 5), 0, Inf)
    θ ~ filldist(Normal(μ, τ), J)
    for i in 1:J
        y[i] ~ Normal(θ[i], σ[i])
    end
end

# ╔═╡ 86cb5e19-49e4-4e5e-8b89-e76936932055
rng2 = Random.MersenneTwister(16653);

# ╔═╡ 85bbcba7-c0f9-4c86-9cdf-a27055d3d448
begin
    param_mod_turing = model_turing(y, σ)
    sampler = NUTS(ndraws_warmup, 0.8)

    turing_chns = Turing.sample(
        rng2, model_turing(y, σ), sampler, MCMCThreads(), ndraws, nchains
    )
end;

# ╔═╡ bd4ab044-51ce-4af9-83b2-bd8fc827f810
md"""
Most ArviZ functions work fine with `Chains` objects from Turing:
"""

# ╔═╡ 500f4e0d-0a36-4b5c-8900-667560fbf1d4
begin
    plot_autocorr(turing_chns; var_names=["μ", "τ"])
    gcf()
end

# ╔═╡ 1129ad94-f65a-4332-b354-21bcf7e53541
md"""
### Convert to `InferenceData`

For much more powerful querying, analysis and plotting, we can use built-in ArviZ utilities to convert `Chains` objects to xarray datasets.
Note we are also giving some information about labelling.

ArviZ is built to work with [`InferenceData`](https://julia.arviz.org/stable/reference/#ArviZ.InferenceData) (a netcdf datastore that loads data into `xarray` datasets), and the more *groups* it has access to, the more powerful analyses it can perform.
"""

# ╔═╡ 803efdd8-656e-4e37-ba36-81195d064972
idata_turing_post = from_mcmcchains(
    turing_chns;
    coords=Dict("school" => schools),
    dims=Dict("y" => ["school"], "σ" => ["school"], "θ" => ["school"]),
    library="Turing",
)

# ╔═╡ 79f342c8-0738-432b-bfd7-2da25e50fa91
md"""
Each group is an [`ArviZ.Dataset`](https://julia.arviz.org/stable/reference/#ArviZ.Dataset) (a thinly wrapped `xarray.Dataset`).
We can view a summary of the dataset.
"""

# ╔═╡ 6209f947-5001-4507-b3e8-9747256f3328
idata_turing_post.posterior

# ╔═╡ 26692722-db0d-41de-b58c-339b639f6948
md"""
Here is a plot of the trace. Note the intelligent labels.
"""

# ╔═╡ 14046b83-9c0a-4d33-ae4e-36c7d6f1b2e6
begin
    plot_trace(idata_turing_post)
    gcf()
end

# ╔═╡ 737f319c-1ddd-45f2-8d10-aaecdc1334be
md"We can also generate summary stats..."

# ╔═╡ 5be3131d-db8b-44c8-9fcc-571d68695148
summarystats(idata_turing_post)

# ╔═╡ 0787a842-63ca-407d-9d07-a3d949940f92
md"...and examine the energy distribution of the Hamiltonian sampler."

# ╔═╡ 6e8343c8-bee3-4e1d-82d6-1885bfd1dbec
begin
    plot_energy(idata_turing_post)
    gcf()
end

# ╔═╡ cba6f6c9-82c4-4957-acc3-36e9f1c95d76
md"""
### Additional information in Turing.jl

With a few more steps, we can use Turing to compute additional useful groups to add to the `InferenceData`.

To sample from the prior, one simply calls `sample` but with the `Prior` sampler:
"""

# ╔═╡ 79374a8b-98c1-443c-acdd-0ee00bd42b38
prior = Turing.sample(rng2, param_mod_turing, Prior(), ndraws);

# ╔═╡ c23562ba-3b0b-49d9-bb41-2fedaa9e9500
md"""
To draw from the prior and posterior predictive distributions we can instantiate a "predictive model", i.e. a Turing model but with the observations set to `missing`, and then calling `predict` on the predictive model and the previously drawn samples:
"""

# ╔═╡ d0b447fe-26b2-48e0-bcec-d9f04697973a
begin
    # Instantiate the predictive model
    param_mod_predict = model_turing(similar(y, Missing), σ)
    # and then sample!
    prior_predictive = Turing.predict(rng2, param_mod_predict, prior)
    posterior_predictive = Turing.predict(rng2, param_mod_predict, turing_chns)
end;

# ╔═╡ 4d2fdcbe-c1d4-43b6-b382-f2e956b952a1
md"""
And to extract the pointwise log-likelihoods, which is useful if you want to compute metrics such as [`loo`](https://julia.arviz.org/stable/reference/#ArviZ.loo),
"""

# ╔═╡ 5a075722-232f-40fc-a499-8dc5b0c2424a
log_likelihood = let
    log_likelihood = Turing.pointwise_loglikelihoods(
        param_mod_turing, MCMCChains.get_sections(turing_chns, :parameters)
    )
    # Ensure the ordering of the loglikelihoods matches the ordering of `posterior_predictive`
    ynames = string.(keys(posterior_predictive))
    log_likelihood_y = getindex.(Ref(log_likelihood), ynames)
    # Reshape into `(nchains, ndraws, size(y)...)`
    Dict("y" => permutedims(cat(log_likelihood_y...; dims=3), (2, 1, 3)))
end;

# ╔═╡ 1b5af2c3-f2ce-4e9d-9ad7-ac287a9178e2
md"This can then be included in the [`from_mcmcchains`](https://julia.arviz.org/stable/reference/#ArviZ.from_mcmcchains) call from above:"

# ╔═╡ b38c7a43-f00c-43c0-aa6b-9c581d6d0c73
idata_turing = from_mcmcchains(
    turing_chns;
    posterior_predictive,
    log_likelihood,
    prior,
    prior_predictive,
    observed_data=Dict("y" => y),
    coords=Dict("school" => schools),
    dims=Dict("y" => ["school"], "σ" => ["school"], "θ" => ["school"]),
    library="Turing",
)

# ╔═╡ a3b71e4e-ac1f-4404-b8a2-35d03d326774
md"""
Then we can for example compute the expected *leave-one-out (LOO)* predictive density, which is an estimate of the out-of-distribution predictive fit of the model:
"""

# ╔═╡ f552b5b5-9744-41df-af90-46405367ea0b
loo(idata_turing; pointwise=false) # higher is better

# ╔═╡ 9d3673f5-b57b-432e-944a-70b23643128a
md"""
If the model is well-calibrated, i.e. it replicates the true generative process well, the CDF of the pointwise LOO values should be similarly distributed to a uniform distribution.
This can be inspected visually:
"""

# ╔═╡ 05c9be29-7758-4324-971c-5579f99aaf9d
begin
    plot_loo_pit(idata_turing; y="y", ecdf=true)
    gcf()
end

# ╔═╡ 98acc304-22e3-4e6b-a2f4-d22f6847145b
md"""
## [Plotting with CmdStan.jl outputs](#Plotting-with-CmdStan.jl-outputs)

CmdStan.jl and StanSample.jl also default to producing `Chains` outputs, and we can easily plot these chains.

Here is the same centered eight schools model:
"""

# ╔═╡ b46af168-1ce3-4058-a014-b66c645a6e0d
begin
    schools_code = """
    data {
      int<lower=0> J;
      real y[J];
      real<lower=0> sigma[J];
    }

    parameters {
      real mu;
      real<lower=0> tau;
      real theta[J];
    }

    model {
      mu ~ normal(0, 5);
      tau ~ cauchy(0, 5);
      theta ~ normal(mu, tau);
      y ~ normal(theta, sigma);
    }

    generated quantities {
        vector[J] log_lik;
        vector[J] y_hat;
        for (j in 1:J) {
            log_lik[j] = normal_lpdf(y[j] | theta[j], sigma[j]);
            y_hat[j] = normal_rng(theta[j], sigma[j]);
        }
    }
    """

    schools_data = Dict("J" => J, "y" => y, "sigma" => σ)
    stan_chns = mktempdir() do path
        stan_model = Stanmodel(;
            model=schools_code,
            name="schools",
            nchains,
            num_warmup=ndraws_warmup,
            num_samples=ndraws,
            output_format=:mcmcchains,
            random=CmdStan.Random(28983),
            tmpdir=path,
        )
        _, chns, _ = stan(stan_model, schools_data; summary=false)
        return chns
    end
end;

# ╔═╡ ab145e41-b230-4cad-bef5-f31e0e0770d4
begin
    plot_density(stan_chns; var_names=["mu", "tau"])
    gcf()
end

# ╔═╡ ffc7730c-d861-48e8-b173-b03e0542f32b
md"""
Again, converting to `InferenceData`, we can get much richer labelling and mixing of data.
Note that we're using the same [`from_cmdstan`](https://julia.arviz.org/stable/reference/#ArviZ.from_cmdstan) function used by ArviZ to process cmdstan output files, but through the power of dispatch in Julia, if we pass a `Chains` object, it instead uses ArviZ.jl's overloads, which forward to `from_mcmcchains`.
"""

# ╔═╡ 020cbdc0-a0a2-4d20-838f-c99b541d5832
idata_stan = from_cmdstan(
    stan_chns;
    posterior_predictive="y_hat",
    observed_data=Dict("y" => schools_data["y"]),
    log_likelihood="log_lik",
    coords=Dict("school" => schools),
    dims=Dict(
        "y" => ["school"],
        "sigma" => ["school"],
        "theta" => ["school"],
        "log_lik" => ["school"],
        "y_hat" => ["school"],
    ),
)

# ╔═╡ e44b260c-9d2f-43f8-a64b-04245a0a5658
md"""Here is a plot showing where the Hamiltonian sampler had divergences:"""

# ╔═╡ 5070bbbc-68d2-49b8-bd91-456dc0da4573
begin
    plot_pair(
        idata_stan;
        coords=Dict("school" => ["Choate", "Deerfield", "Phillips Andover"]),
        divergences=true,
    )
    gcf()
end

# ╔═╡ 2674d532-a337-471e-8ba7-02b430f49f12
md"""
## [Plotting with Soss.jl outputs](#Plotting-with-Soss.jl-outputs)

With Soss, we can define our model for the posterior and easily use it to draw samples from the prior, prior predictive, posterior, and posterior predictive distributions.

First we define our model:
"""

# ╔═╡ 14408abe-a16f-4cc0-a6f3-0bb2645653b7
constant_data = (; J, σ);

# ╔═╡ 446341da-902e-474b-b6dc-b085ef74a99b
observed_data = (; y);

# ╔═╡ 9daec35c-3d6e-443c-87f9-213d51964f75
model_soss = Soss.@model (J, σ) begin
    μ ~ Soss.Normal(; μ=0, σ=5)
    τ ~ HalfCauchy(; σ=5)
    θ ~ Soss.Normal(; μ, σ=τ) |> iid(J)
    y ~ For(1:J) do j
        Soss.Normal(; μ=θ[j], σ=σ[j])
    end
end;

# ╔═╡ 6a78e4a8-86c4-4438-b9bb-7c433d2bc8c8
param_mod_soss = model_soss(; constant_data...)

# ╔═╡ a94b3ad6-2305-43f6-881f-d871286b906a
md"Then we draw from the prior and prior predictive distributions."

# ╔═╡ bfdfdff7-6551-4dd7-a7c8-69f8b57272ec
rng3 = MersenneTwister(5298);

# ╔═╡ d385ceea-beb5-4ec7-b50d-be691266440b
prior_priorpred = let Normal = Soss.Normal
    [rand(rng3, param_mod_soss, ndraws) for _ in 1:nchains]
end

# ╔═╡ 152822f9-20f7-4fe9-ad9f-bae31945e981
md"""Next, we draw from the posterior using [DynamicHMC.jl](https://github.com/tpapp/DynamicHMC.jl)."""

# ╔═╡ 2c44f500-c102-4a2e-980e-42afcfbba5bb
post = Soss.sample(rng3, param_mod_soss | observed_data, dynamichmc(), ndraws, nchains)

# ╔═╡ e2260ee1-c792-4b74-97c4-bbf6803a56d7
md"Finally, we update the posterior samples with draws from the posterior predictive distribution."

# ╔═╡ 0219a0d0-13ce-43d0-a048-8e95e55500db
postpred = let
    mod_pred = Soss.predictive(model_soss, :μ, :τ, :θ)
    map(getchains(post)) do chain
        [rand(rng3, mod_pred(; constant_data..., draw...)) for draw in chain]
    end
end

# ╔═╡ 500c4b4d-25eb-4f27-8e8b-7287bfeddf92
md"""
Each Soss draw is a `NamedTuple`.
We can plot the rank order statistics of the posterior to identify poor convergence:
"""

# ╔═╡ eac7b059-129d-472b-a69e-b1611c7cc703
begin
    plot_rank(post; var_names=["μ", "τ"])
    gcf()
end

# ╔═╡ 17d9fff5-d8f6-42c5-8cd1-70ecb34084c7
md"Now we combine all of the samples to an `InferenceData`:"

# ╔═╡ ddb1373a-f180-4efc-9193-905d03be4d8a
idata_soss = from_samplechains(
    post;
    posterior_predictive=postpred,
    prior=prior_priorpred,
    prior_predictive=[:y],
    observed_data,
    constant_data,
    coords=Dict("school" => schools),
    dims=Dict("y" => ["school"], "σ" => ["school"], "θ" => ["school"]),
    library=Soss,
)

# ╔═╡ bfa2c30f-f50c-45de-bd2a-da199ce67cfb
md"We can compare the prior and posterior predictive distributions:"

# ╔═╡ e56ee0c8-bc4a-4b37-bfcc-9ca1c440c1f3
begin
    plot_density(
        [idata_soss.posterior_predictive, idata_soss.prior_predictive];
        data_labels=["Post-pred", "Prior-pred"],
        var_names=["y"],
    )
    gcf()
end

# ╔═╡ ac2b4378-bd1c-4164-af05-d9a35b1bb08f
md"## [Environment](#environment)"

# ╔═╡ fd84237a-7a19-4bbb-a702-faa31075ecbc
with_terminal(Pkg.status; color=false)

# ╔═╡ ad29b7f3-6f5e-4b04-bf70-2308a7d110c0
with_terminal(versioninfo)

# ╔═╡ Cell order:
# ╟─a23dfd65-50a8-4872-8c41-661e96585aca
# ╟─d2eedd48-48c6-4fcd-b179-6be7fe68d3d6
# ╠═ac57d957-0bdd-457a-ac15-9a4f94f0c785
# ╠═467c2d13-6bfe-4feb-9626-fb14796168aa
# ╠═06b00794-e97f-472b-b526-efe4815103f8
# ╟─5acbfd9a-cfea-4c3c-a3f0-1744eb7e4e27
# ╠═efb3f0af-9fac-48d8-bbb2-2dd6ebd5e4f6
# ╠═401e9b91-0bca-4369-8d36-3d9f0b3ad60b
# ╟─2c718ea5-2800-4df6-b62c-e0a9e440a1c3
# ╠═49f19c17-ac1d-46b5-a655-4376b7713244
# ╟─a9789109-2b90-40f7-926c-e7c87025d15f
# ╠═69124a94-a6aa-43f6-8d4f-fa9a6b1cfef1
# ╟─10986e0e-55f7-40ea-ab63-274fbab3126d
# ╠═f383d541-e22d-44b4-b8cb-28b3d67944a1
# ╠═86cb5e19-49e4-4e5e-8b89-e76936932055
# ╠═85bbcba7-c0f9-4c86-9cdf-a27055d3d448
# ╟─bd4ab044-51ce-4af9-83b2-bd8fc827f810
# ╠═500f4e0d-0a36-4b5c-8900-667560fbf1d4
# ╟─1129ad94-f65a-4332-b354-21bcf7e53541
# ╠═803efdd8-656e-4e37-ba36-81195d064972
# ╟─79f342c8-0738-432b-bfd7-2da25e50fa91
# ╠═6209f947-5001-4507-b3e8-9747256f3328
# ╟─26692722-db0d-41de-b58c-339b639f6948
# ╠═14046b83-9c0a-4d33-ae4e-36c7d6f1b2e6
# ╟─737f319c-1ddd-45f2-8d10-aaecdc1334be
# ╠═5be3131d-db8b-44c8-9fcc-571d68695148
# ╟─0787a842-63ca-407d-9d07-a3d949940f92
# ╠═6e8343c8-bee3-4e1d-82d6-1885bfd1dbec
# ╟─cba6f6c9-82c4-4957-acc3-36e9f1c95d76
# ╠═79374a8b-98c1-443c-acdd-0ee00bd42b38
# ╟─c23562ba-3b0b-49d9-bb41-2fedaa9e9500
# ╠═d0b447fe-26b2-48e0-bcec-d9f04697973a
# ╟─4d2fdcbe-c1d4-43b6-b382-f2e956b952a1
# ╠═5a075722-232f-40fc-a499-8dc5b0c2424a
# ╟─1b5af2c3-f2ce-4e9d-9ad7-ac287a9178e2
# ╠═b38c7a43-f00c-43c0-aa6b-9c581d6d0c73
# ╟─a3b71e4e-ac1f-4404-b8a2-35d03d326774
# ╠═f552b5b5-9744-41df-af90-46405367ea0b
# ╟─9d3673f5-b57b-432e-944a-70b23643128a
# ╠═05c9be29-7758-4324-971c-5579f99aaf9d
# ╠═98acc304-22e3-4e6b-a2f4-d22f6847145b
# ╠═b46af168-1ce3-4058-a014-b66c645a6e0d
# ╠═ab145e41-b230-4cad-bef5-f31e0e0770d4
# ╟─ffc7730c-d861-48e8-b173-b03e0542f32b
# ╠═020cbdc0-a0a2-4d20-838f-c99b541d5832
# ╟─e44b260c-9d2f-43f8-a64b-04245a0a5658
# ╠═5070bbbc-68d2-49b8-bd91-456dc0da4573
# ╠═2674d532-a337-471e-8ba7-02b430f49f12
# ╠═14408abe-a16f-4cc0-a6f3-0bb2645653b7
# ╠═446341da-902e-474b-b6dc-b085ef74a99b
# ╠═9daec35c-3d6e-443c-87f9-213d51964f75
# ╠═6a78e4a8-86c4-4438-b9bb-7c433d2bc8c8
# ╟─a94b3ad6-2305-43f6-881f-d871286b906a
# ╠═bfdfdff7-6551-4dd7-a7c8-69f8b57272ec
# ╠═d385ceea-beb5-4ec7-b50d-be691266440b
# ╟─152822f9-20f7-4fe9-ad9f-bae31945e981
# ╠═2c44f500-c102-4a2e-980e-42afcfbba5bb
# ╟─e2260ee1-c792-4b74-97c4-bbf6803a56d7
# ╠═0219a0d0-13ce-43d0-a048-8e95e55500db
# ╟─500c4b4d-25eb-4f27-8e8b-7287bfeddf92
# ╠═eac7b059-129d-472b-a69e-b1611c7cc703
# ╟─17d9fff5-d8f6-42c5-8cd1-70ecb34084c7
# ╠═ddb1373a-f180-4efc-9193-905d03be4d8a
# ╟─bfa2c30f-f50c-45de-bd2a-da199ce67cfb
# ╠═e56ee0c8-bc4a-4b37-bfcc-9ca1c440c1f3
# ╟─ac2b4378-bd1c-4164-af05-d9a35b1bb08f
# ╠═56a39a90-0594-48f4-ba04-f7b612019cd1
# ╠═2c0f3d74-8b1e-4648-b0fa-f8050035b0fd
# ╠═fd84237a-7a19-4bbb-a702-faa31075ecbc
# ╠═ad29b7f3-6f5e-4b04-bf70-2308a7d110c0


text:
using Base: @enum
@enum InputPairs begin
INPUT_PAIR_INVALID = 0
    QT_INPUTS
    PQ_INPUTS
    QSmolar_INPUTS
    QSmass_INPUTS
    HmolarQ_INPUTS
    HmassQ_INPUTS
    DmolarQ_INPUTS
    DmassQ_INPUTS

    PT_INPUTS

    DmassT_INPUTS
    DmolarT_INPUTS
    HmolarT_INPUTS
    HmassT_INPUTS
    SmolarT_INPUTS
    SmassT_INPUTS
    TUmolar_INPUTS
    TUmass_INPUTS

    DmassP_INPUTS
    DmolarP_INPUTS
    HmassP_INPUTS
    HmolarP_INPUTS
    PSmass_INPUTS
    PSmolar_INPUTS
    PUmass_INPUTS
    PUmolar_INPUTS

    HmassSmass_INPUTS
    HmolarSmolar_INPUTS
    SmassUmass_INPUTS
    SmolarUmolar_INPUTS

    DmassHmass_INPUTS
    DmolarHmolar_INPUTS
    DmassSmass_INPUTS
    DmolarSmolar_INPUTS
    DmassUmass_INPUTS
    DmolarUmolar_INPUTS
end


text:
@testset "977.squares-of-a-sorted-array.jl" begin
    @test squares_of_a_sorted_array([-4, -1, 0, 3, 10]) == [0, 1, 9, 16, 100]
    @test squares_of_a_sorted_array([-7, -3, 2, 3, 11]) == [4, 9, 9, 49, 121]
end


text:
function moora(decisionMat::DataFrame, weights::Array{Float64,1})

    w = unitize(weights)

    nalternatives, ncriteria = size(decisionMat)

    normalizedMat = normalize(decisionMat)
    weightednormalizedMat = w * normalizedMat

    cmaxs = colmaxs(weightednormalizedMat)

    refmat = similar(weightednormalizedMat)

    for rowind in 1:nalternatives
        refmat[rowind, :] .= cmaxs - weightednormalizedMat[rowind, :]
    end

    scores = rmaxs = rowmaxs(refmat)

    bestIndex = sortperm(rmaxs) |> first

    result = MooraResult(
       decisionMat,
       w,
       weightednormalizedMat,
       refmat,
       scores,
       bestIndex
   )

    return result
end

max_stars_repo_name: ersagunkuruca/JMcDM

max_issues_repo_name: ersagunkuruca/JMcDM

max_forks_repo_name: ersagunkuruca/JMcDM


text:
# ============================================================
# Quadrature Methods
# ============================================================

"""
$(SIGNATURES)

Gauss-Legendre quadrature

* @args n : quadrature order (MUST be even)
* @return points : quadrature points in 3D coordinate
* @return weights : quadrature weights
"""
function legendre_quadrature(n::T) where {T<:Integer}

    pointsxyz = zeros(n * n, 3)
    weights = zeros(n * n)

    # construct Gauss quadrature
    mu, gaussweights = gausslegendre(n)

    # transform between (mu,phi) and (x,y,z)
    phi = [(k + 0.5) * pi / n for k = 0:2*n-1] # equidistance in z axis
    range = 1:n÷2 # only use upper half of the sphere as quadrature point due to pseudo 3D

    x = sqrt.(1.0 .- mu[range] .^ 2) .* cos.(phi)'
    y = sqrt.(1.0 .- mu[range] .^ 2) .* sin.(phi)'
    z = mu[range] .* ones(size(phi))'
    weights = 2.0 * pi / n * repeat(gaussweights[range], 1, 2 * n)

    # assign
    pointsxyz[:, 1] .= x[:]
    pointsxyz[:, 2] .= y[:]
    pointsxyz[:, 3] .= z[:]
    weights = weights[:]

    return pointsxyz, weights

end

"""
$(SIGNATURES)

Octaeder quadrature

* @args n : quadrature order
* @args slerpflag : flag of spherical linear interpolation
* @return points
* @return triangulation
"""
function octa_quadrature(n::T, slerpflag = true::Bool) where {T<:Integer}

    # integral range
    pt0 = [0.0, 0.0, 1.0]
    pt1 = [0.0, 1.0, 0.0]
    pt2 = [1.0, 0.0, 0.0]

    # slerp / linspace
    if slerpflag
        pts01 = slerp(pt0, pt1, n)
        pts02 = slerp(pt0, pt2, n)
    else
        pts01 = linspace(pt0, pt1, n)
        pts02 = linspace(pt0, pt2, n)
    end

    nptsoctant = Int64(n * (n + 1) / 2)
    pts = zeros(3, nptsoctant)

    # generate points in planar geometry
    counter = 0
    for i = 1:n
        if slerpflag
            if i == 1
                tmp = pts01[:, 1]
            else
                tmp = slerp(pts01[:, i], pts02[:, i], i)
            end
        else
            tmp = linspace(pts01[i], pts02[i], i)
        end
        for j = 1:i
            counter += 1
            if slerpflag
                pts[:, counter] = tmp[:, j]
            else
                pts[:, counter] = tmp[j]
            end
        end

    end

    # project points onto sphere
    for i = 1:nptsoctant
        pts[:, i] = pts[:, i] / norm(pts[:, i])
    end

    # enumerate over points and write their connectivity
    ids = zeros(Int64, n, n) # matrix that assigns an ID to the points
    nTrianglesOctant = Int64(n * n - 2 * n + 1)
    triangles = zeros(Int64, 3, nTrianglesOctant) # matrix that contains all triangles

    counter = 0
    for i = 1:n
        for j = 1:i
            counter += 1
            ids[i, j] = counter
        end
    end

    # create triangles
    counter = 0
    tmp = zeros(Int64, 1)
    for i = 1:n
        for j = 1:i-1
            tmp = [ids[i, j], ids[i, j+1], ids[i-1, j]]
            counter += 1
            triangles[:, counter] = tmp
        end
        if i < n
            for j = 1:i-1
                tmp = [ids[i, j], ids[i, j+1], ids[i+1, j+1]]
                counter += 1
                triangles[:, counter] = tmp
            end
        end
    end

    # now we have the quadrature points and triangles for a single octant
    ptsAll = deepcopy(pts)

    tmp = deepcopy(pts)
    tmp[1, :] *= -1.0
    ptsAll = deepcopy(hcat(ptsAll, tmp))

    tmp = deepcopy(pts)
    tmp[2, :] *= -1.0
    ptsAll = deepcopy(hcat(ptsAll, tmp))

    tmp = deepcopy(pts)
    tmp[1, :] *= -1.0
    tmp[2, :] *= -1.0
    ptsAll = deepcopy(hcat(ptsAll, tmp))

    tmp = deepcopy(ptsAll)
    tmp[3, :] *= -1.0
    ptsAll = deepcopy(hcat(ptsAll, tmp))

    trianglesAll = deepcopy(triangles)
    for i = 1:7
        trianglesAll = (hcat(trianglesAll, triangles .+ i .* nptsoctant))
    end
    ptsAll, triangulation = unique(ptsAll, trianglesAll)

    points = permutedims(ptsAll)
    triangulation = permutedims(triangulation)

    return points, triangulation

end

"""
$(SIGNATURES)

Create quadrature weights from points and triangulation

* @arg xyz : quadrature points
* @arg triangles : triangulation
* @return weights : quadrature weights
"""
function quadrature_weights(xyz::X, triangles::Y) where {X<:AA{<:Real,2},Y<:AA{<:Integer,2}}
    weights = zeros(axes(xyz, 1))
    nTriangles = size(triangles, 1)
    xy = zeros(3)
    yz = zeros(3)
    zx = zeros(3)
    mid = zeros(3)

    for n = 1:nTriangles
        # get three nodes of a triangle
        i, j, k = triangles[n, :]

        # get the corresponding points
        x = xyz[i, :]
        y = xyz[j, :]
        z = xyz[k, :]

        # Now get the midpoint of the triangle and the midpoints along the lines
        mid = (x + y + z) / 3.0
        xy = (x + y) / 2.0
        yz = (y + z) / 2.0
        zx = (z + x) / 2.0

        # These points still have to be projected onto the sphere
        mid = mid / norm(mid, 2)
        xy = xy / norm(xy, 2)
        yz = yz / norm(yz, 2)
        zx = zx / norm(zx, 2)

        # By these four points, plus x,y,z we can span 6 triangles
        # Each of these triangles is assigned to one of the three nodes of the triangle x,y,z.
        # the area = the weight

        # i
        weights[i] += area(x, mid, xy, :sphere)
        weights[i] += area(x, mid, zx, :sphere)

        # j
        weights[j] += area(y, mid, xy, :sphere)
        weights[j] += area(y, mid, yz, :sphere)

        # k
        weights[k] += area(z, mid, yz, :sphere)
        weights[k] += area(z, mid, zx, :sphere)
    end

    return weights
end


text:
struct None <: AbstractSmoothing end

function smooth!(::None,klocal::Integer,kglobal::Integer,x::AbstractVector,gprev::AbstractVector,gcurr::AbstractVector)
    gcurr .= gprev
end


text: #TODO: DO TESTS

max_stars_repo_name: adamwillisXanax/SolvencyPredictor

max_issues_repo_name: adamwillisXanax/SolvencyPredictor

max_forks_repo_name: adamwillisXanax/SolvencyPredictor


text:
using SparseArrays
using LinearAlgebra
using Plots

include("vekit.jl")

let
    # Metrizen definieren
    A=[23 4 3; 5 4 5; -30 -10 -10]
    B=[-1 0 1; 0 1 0; 1 0 -1]
    C=[2 -3 1; 3 1 3; -5 2 -4]
    # y_0 definieren
    y=[8; 3; 5]/norm([8; 3; 5])
    # epsilon definieren
    e=10^-10
    # Eigenwert und eigenvektor berechnen, Fehler auf Eigenwert berechnen
    l_A,y_A=vekit(A,y,e)
    f_A=map(x -> abs((x-maximum(eigen(A).values))/maximum(eigen(A).values)), l_A)
    l_B,y_B=vekit(B,y,e)
    f_B=map(x -> abs((x-maximum(eigen(B).values))/maximum(eigen(B).values)), l_B)
    l_C,y_C=vekit(C,y,e)
    f_C=map(x -> abs((x-maximum(eigen(C).values))/maximum(eigen(C).values)), l_C)

    #TODO: Plot
end

max_stars_repo_name: SoulKathi/vekit

max_issues_repo_name: SoulKathi/vekit

max_forks_repo_name: SoulKathi/vekit


text:
struct InstrumentalResponse{T<:Real}
    amp::Vector{T}
    width::Vector{T}    #in pixels
    location::Vector{T} #in pixels
end

function InstrumentalResponse()
    InstrumentalResponse(Float64[],Float64[],Float64[])
end

function (IR::InstrumentalResponse)(lambda,d=1.0)
    #d converts from pixels to wavelength
    sum([a*exp((-((d*l) - lambda)^2)/(d*w)^2) for (a,l,w) in zip(IR.amp,IR.location,IR.width)])
end

function kernel(lambda,IR::InstrumentalResponse,d=1.0)
    length(IR.location) == 0 && return Float64[]
    dL = lambda[2] - lambda[1]
    lmax = maximum(abs.(d*IR.location)) + 6*maximum(sqrt.(((d*IR.width).^2)/2))
    Lp = 0.0:dL:lmax
    L= vcat(-reverse(Lp)[1:end-1],Lp)
    k = IR.(L,d)
    return reflect(centered(k./sum(k))) #reflect for convolution not correlation
end

function kernel(lambda,IR::T,d=1.0) where T<:Real
    return reflect(Images.Kernel.gaussian((IR/d,), (21,)))
end

struct ExperimentalSpectra{T<:Real}
    lambda::Vector{T}
    data::Vector{T}
    err::Vector{T}
end

function Base.vcat(S::ExperimentalSpectra...)
    lambda = vcat((s.lambda for s in S)...)
    data = vcat((s.data for s in S)...)
    err = vcat((s.err for s in S)...)
    return ExperimentalSpectra(lambda,data,err)
end

struct TheoreticalSpectra{T<:Real}
    lambda::Vector{T}
    spec::Vector{T}
end

function Base.vcat(S::TheoreticalSpectra...)
    lambda = vcat((s.lambda for s in S)...)
    spec = vcat((s.spec for s in S)...)
    return TheoreticalSpectra(lambda,spec)
end


text:
# No idea what kind of joint this actually is...
@inline function getVelocityDelta(joint::Rotational1, body1::AbstractBody, body2::Body{T}, ω) where T
    #TODO define this function
    @error("not defined for rot2")
end

@inline function getPositionDelta(joint::Rotational1, body1::AbstractBody, body2::Body{T}, θ) where T
    #TODO define this function
    @error("not defined for rot2")
end

@inline function setForce!(joint::Rotational1, body1::Body, body2::Body{T}, τ::SVector{2,T}, No) where T
    τ1 = vrotate(joint.V12' * -τ, body1.q[No] * joint.qoff)
    τ2 = -τ1

    body1.τ[No] = τ1
    body2.τ[No] = τ2
    return
end

@inline function setForce!(joint::Rotational1, body1::Origin, body2::Body{T}, τ::SVector{2,T}, No) where T
    body2.τ[No] = vrotate(joint.V12' * τ, joint.qoff)
    return
end

@inline function minimalCoordinates(joint::Rotational1, body1::Body, body2::Body, No)
    joint.V12 * (VLᵀmat(joint.qoff) * Lᵀmat(body1.q[No]) * body2.q[No])
end

@inline function minimalCoordinates(joint::Rotational1, body1::Origin, body2::Body, No)
    joint.V12 * (VLᵀmat(joint.qoff) * body2.q[No])
end

@inline function g(joint::Rotational1, body1::Body, body2::Body, Δt, No)
    joint.V3 * (VLᵀmat(joint.qoff) * Lᵀmat(getq3(body1, Δt)) * getq3(body2, Δt))
end

@inline function g(joint::Rotational1, body1::Origin, body2::Body, Δt, No)
    joint.V3 * VLᵀmat(joint.qoff) * getq3(body2, Δt)
end

@inline function ∂g∂posa(joint::Rotational1{T}, body1::Body, body2::Body, No) where T
    if body2.id == joint.cid
        X = @SMatrix zeros(T, 1, 3)
        R = -joint.V3 * VLᵀmat(joint.qoff) * Rmat(body2.q[No]) * RᵀVᵀmat(body1.q[No])

        return [X R]
    else
        return ∂g∂posa(joint)
    end
end

@inline function ∂g∂posb(joint::Rotational1{T}, body1::Body, body2::Body, No) where T
    if body2.id == joint.cid
        X = @SMatrix zeros(T, 1, 3)
        R = joint.V3 * VLᵀmat(joint.qoff) * Lᵀmat(body1.q[No]) * LVᵀmat(body2.q[No])

        return [X R]
    else
        return ∂g∂posb(joint)
    end
end

@inline function ∂g∂vela(joint::Rotational1{T}, body1::Body, body2::Body, Δt, No) where T
    if body2.id == joint.cid
        V = @SMatrix zeros(T, 1, 3)
        Ω = joint.V3 * VLᵀmat(joint.qoff) * Rmat(ωbar(body2, Δt)) * Rmat(body2.q[No]) * Rᵀmat(body1.q[No]) * Tmat(T) * derivωbar(body1, Δt)

        return [V Ω]
    else
        return ∂g∂vela(joint)
    end
end

@inline function ∂g∂velb(joint::Rotational1{T}, body1::Body, body2::Body, Δt, No) where T
    if body2.id == joint.cid
        V = @SMatrix zeros(T, 1, 3)
        Ω = joint.V3 * VLᵀmat(joint.qoff) * Lᵀmat(ωbar(body1, Δt)) * Lᵀmat(body1.q[No]) * Lmat(body2.q[No]) * derivωbar(body2, Δt)

        return [V Ω]
    else
        return ∂g∂velb(joint)
    end
end

@inline function ∂g∂posb(joint::Rotational1{T}, body1::Origin, body2::Body, No) where T
    if body2.id == joint.cid
        X = @SMatrix zeros(T, 1, 3)
        R = joint.V3 * VLᵀmat(joint.qoff) * LVᵀmat(body2.q[No])

        return [X R]
    else
        return ∂g∂posb(joint)
    end
end

@inline function ∂g∂velb(joint::Rotational1{T}, body1::Origin, body2::Body, Δt, No) where T
    if body2.id == joint.cid
        V = @SMatrix zeros(T, 1, 3)
        Ω = joint.V3 * VLᵀmat(joint.qoff) * Lmat(body2.q[No]) * derivωbar(body2, Δt)

        return [V Ω]
    else
        return ∂g∂velb(joint)
    end
end


text:
path = "/home/user/programming/GraphEvolve.jl/src/"
push!(LOAD_PATH, path)
using Documenter, GraphEvolve

makedocs(
    sitename = "GraphEvolve.jl",
    authors = "Cameron Perot",
    pages = [
        "Home" => "index.md"
        "Manual" => Any[
            "man/getting_started.md",
            # "man/examples.md",
            "man/AbstractGraphs.md",
            "man/evolution_methods.md",
            "man/edge_methods.md",
            "man/cluster_methods.md",
            "man/analysis_methods.md",
        ]
    ],
)


text:
module Utils

using Test

function with_log_on_error(f, name = f)
    try
        return f()
    catch err
        @error "Error from `$name`" exception = (err, catch_backtrace())
        rethrow()
    end
end

concurrently(functions...; kwargs...) = run_concurrently(collect(functions); kwargs...)

"""
    run_concurrently(thunks; [spawn = false])

Run `thunks` concurrently.
"""
function run_concurrently(functions::AbstractVector; spawn = false)
    tasks = Task[]
    Base.Experimental.@sync begin  # OK for testing
        for f in functions[1:end-1]
            t = if spawn
                Threads.@spawn with_log_on_error(f)
            else
                @async with_log_on_error(f)
            end
            push!(tasks, t)
        end
        global TASKS = tasks  # for debugging
        functions[end]()
    end
end

# This use of `schedule` is wrong but it's handy when debugging implementation
# of concurrent program that is wrong.
function cancel_tasks(tasks = TASKS)
    for t in tasks
        istaskdone(t) && continue
        schedule(t, ErrorException("!!! cancel !!!"); error = true)
    end
    return tasks
end

function random_sleep(spawn::Bool = true)
    if rand() < 0.1
        # no sleep
    elseif spawn && rand(Bool)
        nspins = rand(0:10000)
        for _ in 1:nspins
            GC.safepoint()
            ccall(:jl_cpu_pause, Cvoid, ())
        end
    else
        sleep(rand() / 1_000_000)
    end
end

macro spawn_named(name, spawn, ex = undef)
    if ex === undef
        ex = spawn
        spawn = true
    end
    quote
        if $spawn
            $Threads.@spawn $with_log_on_error($name) do
                $ex
            end
        else
            $Base.@async $with_log_on_error($name) do
                $ex
            end
        end
    end |> esc
end

macro test_error(expr)
    @gensym err tmp
    quote
        local $err = nothing
        $Test.@test try
            $expr
            false
        catch $tmp
            $err = $tmp
            true
        end
        $err
    end |> esc
end

const ⊏ = occursin

end  # module


text:
using SchwarzOpt
using Test

@testset "SchwarzOpt.jl" begin
    include("test_optimizer.jl")
end


text:
module REPLCompletions

    export completions, shell_completions

    using Base.Meta

    function completes_global(x, name)
        return beginswith(x, name) && !in('#', x)
    end

    # REPL Symbol Completions
    function complete_symbol(sym,ffunc)
        # Find module
        strs = split(sym,".")
        # Maybe be smarter in the future
        context_module = Main

        mod = context_module
        lookup_module = true
        t = None
        for name in strs[1:(end-1)]
            s = symbol(name)
            if lookup_module
                if isdefined(mod,s)
                    b = mod.(s)
                    if isa(b,Module)
                        mod = b
                    elseif Base.isstructtype(typeof(b))
                        lookup_module = false
                        t = typeof(b)
                    else
                        # A.B.C where B is neither a type nor a
                        # module. Will have to be revisited if
                        # overloading is allowed
                        return ASCIIString[]
                    end
                else
                    # A.B.C where B doesn't exist in A. Give up
                    return ASCIIString[]
                end
            else
                # We're now looking for a type
                fields = t.names
                found = false
                for i in 1:length(fields)
                    if s == fields[i]
                        t = t.types[i]
                        if !Base.isstructtype(t)
                            return ASCIIString[]
                        end
                        found = true
                        break
                    end
                end
                if !found
                    #Same issue as above, but with types instead of modules
                    return ASCIIString[]
                end
            end
        end

        name = strs[end]

        suggestions = String[]
        if lookup_module
            # We will exclude the results that the user does not want, as well
            # as excluding Main.Main.Main, etc., because that's most likely not what
            # the user wants
            p = s->(ffunc(mod,s) && s != module_name(mod))
            # Looking for a binding in a module
            if mod == context_module
                # Also look in modules we got through `using`
                mods = ccall(:jl_module_usings,Any,(Any,),Main)
                for m in mods
                    ssyms = names(m)
                    filter!(p,ssyms)
                    syms = map!(string,Array(UTF8String,length(ssyms)),ssyms)
                    append!(suggestions,syms[map((x)->completes_global(x,name),syms)])
                end
                ssyms = names(mod,true,true)
                filter!(p,ssyms)
                syms = map!(string,Array(UTF8String,length(ssyms)),ssyms)
            else
                ssyms = names(mod,true,false)
                filter!(p,ssyms)
                syms = map!(string,Array(UTF8String,length(ssyms)),ssyms)
            end
            append!(suggestions,syms[map((x)->completes_global(x,name),syms)])
        else
            # Looking for a member of a type
            fields = t.names
            for field in fields
                s = string(field)
                if beginswith(s,name)
                    push!(suggestions,s)
                end
            end
        end
        suggestions
    end

    function complete_path(path::ByteString)
        matches = ByteString[]
        dir,prefix = splitdir(path)
        if length(dir) == 0
            files = readdir()
        elseif isdir(dir)
            files = readdir(dir)
        else
            return matches
        end
        for file in files
            if beginswith(file, prefix)
                p = joinpath(dir, file)
                push!(matches, isdir(p) ? joinpath(p,"") : p)
            end
        end
        matches
    end

    const non_word_chars = " \t\n\"\\'`@\$><=:;|&{}()[].,+-*/?%^~"

    function completions(string,pos)
        startpos = min(pos,1)
        dotpos = 0

        instring = false
        incmd = false
        escaped = false
        nearquote = false
        i = start(string)
        while i <= pos
            c,j = next(string, i)
            if c == '\\'
                if instring
                    escaped $= true
                end
                nearquote = false
            elseif c == '\''
                if !instring
                    nearquote = true
                end
            elseif c == '"'
                if !escaped && !nearquote && !incmd
                    instring $= true
                end
                escaped = nearquote = false
            elseif c == '`'
                if !escaped && !nearquote && !instring
                    incmd $= true
                  end
                escaped = nearquote = false
            else
                escaped = nearquote = false
            end
            if c < 0x80
                if instring || incmd
                    if in(c, " \t\n\"\\'`@\$><=;|&{(")
                        startpos = j
                    end
                elseif in(c, non_word_chars)
                    if c == '.'
                        dotpos = i
                    else
                        startpos = j
                    end
                end
            end
            i = j
        end

        if instring || incmd
            r = startpos:pos
            paths = complete_path(string[r])
            if instring && length(paths) == 1
                paths[1] *= "\""
            end
            return sort(paths), r
        end

        ffunc = (mod,x)->true
        suggestions = UTF8String[]
        r = rsearch(string,"using",startpos)
        if !isempty(r) && all(isspace,string[nextind(string,last(r)):prevind(string,startpos)])
            # We're right after using. Let's look only for packages
            # and modules we can reach from here

            # If there's no dot, we're in toplevel, so we should
            # also search for packages
            s = string[startpos:pos]
            if dotpos == 0
                append!(suggestions,filter(pname->begin
                    pname[1] != '.' &&
                    pname != "METADATA" &&
                    beginswith(pname,s)
                end,readdir(Pkg.dir())))
            end
            ffunc = (mod,x)->(isdefined(mod,x) && isa(mod.(x),Module))
        end
        if startpos == 0
            pos = -1
        end
        if dotpos == 0
            dotpos = startpos-1
        end
        append!(suggestions,complete_symbol(string[startpos:pos],ffunc))
        sort(unique(suggestions)), (dotpos+1):pos
    end

    function shell_completions(string,pos)
        # First parse everything up to the current position
        scs = string[1:pos]
        args, last_parse = Base.shell_parse(scs,true)
        # Now look at the last this we parsed
        arg = args.args[end].args[end]
        if isa(arg,String)
            # Treat this as a path (perhaps give a list of comands in the future as well?)
            dir,name = splitdir(arg)
            if isempty(dir)
                files = readdir()
            else
                if !isdir(dir)
                    return ([],0:-1)
                end
                files = readdir(dir)
            end
            # Filter out files and directories that do not begin with the partial name we were
            # completing and append "/" to directories to simplify further completion
            ret = map(filter(x->beginswith(x,name),files)) do x
                if !isdir(joinpath(dir,x))
                    return x
                else
                    return x*"/"
                end
            end
            r = (nextind(string,pos-sizeof(name))):pos
            return (ret,r,string[r])
        elseif isexpr(arg,:escape) && (isexpr(arg.args[1],VERSION >= v"0.3-" ? :incomplete : :continue) || isexpr(arg.args[1],:error))
            r = first(last_parse):prevind(last_parse,last(last_parse))
            partial = scs[r]
            ret, range = completions(partial,endof(partial))
            range += first(r)-1
            return (ret,range)
        end
    end
end


text:
# This file contains the links to connect together the tools of DsSimulator.
import Base: put!, take!, close, isready, eltype, isopen, isreadable, iswritable, bind, collect, iterate

"""
    $TYPEDEF

A `Link` connects an Outpin and Inpin for data flow. See [`Outpin`](@ref), [`Inpin`](@ref), [`connect!`](@ref)

# Fields

    $TYPEDFIELDS

# Example
```jldoctest
julia> l = Link{Int}(5)
Link(state:open, eltype:Int64, isreadable:false, iswritable:false)

julia> l = Link{Bool}()
Link(state:open, eltype:Bool, isreadable:false, iswritable:false)
```
"""
mutable struct Link{T}
    "Internal buffer to record data flowing throgh link"
    buffer::Buffer{Cyclic, T, 1}
    "Internal channel for data flow"
    channel::Channel{T}
    "Unique identifier of the outpin of the source of the link"
    masterid::UUID
    "Unique identifier fo the inpin of the destination of the link"
    slaveid::UUID
    "Unique identifier"
    id::UUID
    Link{T}(ln::Int=64) where {T} = new{T}(Buffer(T, ln), Channel{T}(0), uuid4(), uuid4(), uuid4())
end
Link(ln::Int=64) = Link{Float64}(ln)

show(io::IO, link::Link) = print(io, "Link(state:$(isopen(link) ? :open : :closed), eltype:$(eltype(link)), ",
    "isreadable:$(isreadable(link)), iswritable:$(iswritable(link)))")

"""
    $SIGNATURES

Returns element type of `link`.
"""
eltype(link::Link{T}) where {T} = T

##### Link reading writing.
"""
    $SIGNATURES

Puts `val` to `link`. `val` is handed over to the `channel` of `link`. `val` is also written in to the `buffer` of `link`.

!!! warning
    `link` must be writable to put `val`. That is, a runnable task that takes items from the link must be bounded to `link`.

# Example
```jldoctest
julia> l = Link();

julia> t  = @async while true
       item = take!(l)
       item === NaN && break
       println("Took " * string(item))
       end;

julia> bind(l, t);

julia> put!(l, 1.)
Took 1.0
1.0

julia> put!(l, 2.)
Took 2.0
2.0

julia> put!(l, NaN)
NaN
```
"""
function put!(link::Link, val)
    write!(link.buffer, val)
    put!(link.channel, val)
end

"""
    $SIGNATURES

Take an element from `link`.

!!! warning
    `link` must be readable to take value. That is, a runnable task that puts items from the link must be bounded to `link`.

# Example
```jldoctest
julia> l = Link(5);

julia> t = @async for item in 1. : 5.
       put!(l, item)
       end;

julia> bind(l, t);

julia> take!(l)
1.0

julia> take!(l)
2.0
```
"""
function take!(link::Link)
    val = take!(link.channel)
    return val
end

"""
    $SIGNATURES

Closes `link`. All the task bound the `link` is also terminated safely. When closed, it is not possible to take and put element from the `link`. See also: [`take!(link::Link)`](@ref), [`put!(link::Link, val)`](@ref)
```
"""
function close(link::Link)
    channel = link.channel
    iswritable(link) && put!(link, NaN)   # Terminate taker task
    isreadable(link) && collect(link.channel)   # Terminater putter task
    isopen(link) && close(link.channel)  # Close link channel if it is open.
    return
end

##### State check of link.
"""
    $SIGNATURES

Returns `true` if `link` is open. A `link` is open if its `channel` is open.
"""
isopen(link::Link) = isopen(link.channel)

"""
    $SIGNATURES

Returns `true` if `link` is readable. When `link` is readable, data can be read from `link` with `take` function.
"""
isreadable(link::Link) = length(link.channel.cond_put.waitq) > 0

"""
    $SIGNATURES

Returns `true` if `link` is writable. When `link` is writable, data can be written into `link` with `put` function.
"""
iswritable(link::Link) = length(link.channel.cond_take.waitq) > 0

"""
    $SIGNATURES

Returns `true` if the `buffer` of `link` is full.
"""
isfull(link::Link) = isfull(link.buffer)

"""
    $SIGNATURES

Returns all the data of the `buffer` of `link`.
"""
snapshot(link::Link) = link.buffer.data

##### Launching links.

##### Auxilary functions to launch links.
### The `taker` and `puter` functions are just used for troubleshooting purpose.
function taker(link::Link)
    while true
        val = take!(link)
        val === NaN && break  # Poison-pill the tasks to terminate safely.
        @info "Took " val
    end
end

function putter(link::Link, vals)
    for val in vals
        put!(link, val)
    end
end

"""
    $SIGNATURES

Binds `task` to `link`. When `task` is done `link` is closed.
"""
bind(link::Link, task::Task) = bind(link.channel, task)

"""
    $SIGNATURES

Collects all the available data on the `link`.

!!! warning
    To collect all available data from `link`, a task must be bounded to it.

# Example
```jldoctest
julia> l = Link();  # Construct a link.

julia> t = @async for item in 1 : 5  # Construct a task
       put!(l, item)
       end;

julia> bind(l, t);  # Bind it to the link.

julia> take!(l)  # Take element from link.
1.0

julia> take!(l)  # Take again ...
2.0

julia> collect(l)  # Collect remaining data.
3-element Array{Float64,1}:
 3.0
 4.0
 5.0
```
"""
collect(link::Link) = collect(link.channel)

"""
    $SIGNATURES

Constructs a `taker` task and binds it to `link`. The `taker` task reads the data and prints an info message until `missing` is read from the `link`.
"""
function launch(link::Link)
    task = @async taker(link)
    bind(link.channel, task)
    task
end

"""
    $SIGNATURES

Constructs a `putter` task and binds it to `link`. `putter` tasks puts the data in `valrange`.
"""
function launch(link::Link, valrange)
    task = @async putter(link, valrange)
    bind(link.channel, task)
    task
end

function launch(link::Link, taskname::Symbol, valrange)
    msg = "`launch(link, taskname, valrange)` has been deprecated."
    msg *= "Use `launch(link)` to launch taker task, `launch(link, valrange)` to launch putter task"
    @warn msg
end

"""
    $SIGNATURES

Reconstructst the channel of `link` is its channel is closed.
"""
refresh!(l::Link{T}) where {T} = (l.channel = Channel{T}(); l)

"""
    $(SIGNATURES)

Cleans the data on `link`.
"""
clean!(link::Link) = (clean!(link.buffer); link)


text:
include("fft.jl")
include("multigrid.jl")
function get_f{T}(u::AbstractArray{T,2}, eps::T)
"""Get the f with left hand side"""
    mx,my=size(u)
    fft=rFFTWrapper(u)
    lap = similar(u)
    laplacian!(u,lap,fft)
    dx=similar(u)
    dy=similar(u)
    set_not_calculated(fft)
    derivx!(u,dx,fft)
    derivy!(u,dy,fft)
    @inbounds for x in 1:mx
        for y in 1:my
            lap[x,y] += eps * (dx[x,y]^2 + dy[x,y]^2 + lap[x,y] * u[x,y])
        end
    end
    return lap
end

function solver{T}(f::AbstractArray{T,2}, eps::T, tol::T, maxiter=50)
"""Solving the equation with laplacian"""
    mx,my=size(f)
    fft=rFFTWrapper(f)
    u=zeros(f)
    newu=similar(u)
    dx=similar(u)
    dy=similar(u)
    rhs=similar(u)
    iter=0
    err=T[]
    while iter < maxiter
        iter += 1
        derivx!(u,dx,fft)
        derivy!(u,dy,fft)
        @inbounds for x in 1:mx
            for y in 1:my
                rhs[x,y]= (f[x,y] - eps * (dx[x,y]^2 + dy[x,y]^2)) / (one(T) + eps * u[x,y])
            end
        end
        solve_laplacian!(rhs,newu,fft)
        e=vecnorm(newu - u) / vecnorm(f)
        push!(err,e)
        if e < tol
            break
        end
        u,newu=newu,u
    end
    return u,err
end

function iteration{T}(grid::gridWrapper, eps::T, err::Array{T,1}, maxiter=4)
    mx,my=size(grid.f)
    iter=0
    while iter < maxiter
        iter += 1
        derivx!(grid.u,grid.dx,grid.fft)
        derivy!(grid.u,grid.dy,grid.fft)
        @inbounds for x in 1:mx
            for y in 1:my
                grid.rhs[x,y]= (grid.f[x,y] - eps * (grid.dx[x,y]^2 + grid.dy[x,y]^2)) / (one(T) + eps * grid.u[x,y])
            end
        end
        solve_laplacian!(grid.rhs,grid.newu,grid.fft)
        grid.u,grid.newu=grid.newu,grid.u
        push!(err,vecnorm(grid.newu - grid.u) / vecnorm(grid.f))
    end
end

function iteration{T}(grid::gridWrapper, eps::T, maxiter=4)
    mx,my=size(grid.f)
    iter=0
    while iter < maxiter
        iter += 1
        derivx!(grid.u,grid.dx,grid.fft)
        derivy!(grid.u,grid.dy,grid.fft)
        @inbounds for x in 1:mx
            for y in 1:my
                grid.rhs[x,y]= (grid.f[x,y] - eps * (grid.dx[x,y]^2 + grid.dy[x,y]^2)) / (one(T) + eps * grid.u[x,y])
            end
        end
        solve_laplacian!(grid.rhs,grid.newu,grid.fft)
        grid.u,grid.newu=grid.newu,grid.u
    end
end

function multigrid_solver{T}(f::AbstractArray{T,2}, eps::T, tol::T, maxiter=32)
"""Solve the equation with multigrid"""
    mx,my=size(f)
    grid1=gridWrapper(f)
    f2=similar(f,div(mx,2),div(my,2))
    restriction!(f,f2)
    grid2=gridWrapper(f2)
    f3=similar(f,div(mx,4),div(my,4))
    restriction!(f2,f3)
    grid3=gridWrapper(f3)
    f4=similar(f,div(mx,8),div(my,8))
    restriction!(f3,f4)
    grid4=gridWrapper(f4)
    err=T[]
    full_cycle(iteration,grid1,grid2,grid3,grid4,eps,tol)
    iter=0
    while iter < maxiter
        iter += 1
        iteration(grid1,eps,err)
        if vecnorm(grid1.u - grid1.newu) < tol * vecnorm(grid1.f)
            break
        end
    end
    return grid1.u,err
end

max_stars_repo_name: JerryLingjieMei/NonlinearPDE

max_issues_repo_name: JerryLingjieMei/NonlinearPDE

max_forks_repo_name: JerryLingjieMei/NonlinearPDE


text:
using RCall, MixedModels, Test
using StatsBase: zscore
const LMM = LinearMixedModel
const GLMM = GeneralizedLinearMixedModel

@testset "merMod" begin
    # this is available in MixedModels.dataset(:sleepstudy) but with different
    # capitalization than in R
    sleepstudy = rcopy(R"sleepstudy")
    jlmm = fit!(LMM(@formula(Reaction ~ 1 + round(Days) + (1|Subject)),sleepstudy), REML=false)
    @testset "bare model" begin
        @test_throws ArgumentError (@rput jlmm)
    end
    @testset "reversed tuple" begin
        jm = (sleepstudy, jlmm);
        @test_throws ArgumentError (@rput jm)
    end
end


text:
@recipe function plot(l::AbstractLattice; bondcolor=:grey)
    markershape --> :circle        # if markershape is unset, make it :auto
    markercolor --> :black
    markersize --> 5
    grid --> false
    axis --> false
    legend --> false
    if ndims(l) == 3
        showaxis --> false
    end

    # plot sites
    allsites = sites(l)
    @series begin
        seriestype := :scatter
        x = [point(s)[1] for s in allsites]
        y = [point(s)[2] for s in allsites]

        if ndims(l) == 3
            z = [point(s)[3] for s in allsites]
            x,y,z
        else
            x,y
        end
    end

    # plot bonds
    for b in bonds(l)
        all(iszero.(b.wrap)) || continue
        @series begin
            seriestype := :line
            linestyle := :solid
            color := bondcolor

            # get the coordinates to where the bond is pointing
            from_point = point(allsites[from(b)])
            x1 = from_point[1]
            y1 = from_point[2]
            # get the coordinates to where the bond is pointing
            to_point = point(allsites[to(b)])
            x2 = to_point[1]
            y2 = to_point[2]

            if ndims(l) == 3
                z1 = from_point[3]
                z2 = to_point[3]
                [x1,x2], [y1,y2], [z1,z2]
            else
                [x1,x2], [y1,y2]
            end

        end
    end
end