text:
function to_xml(doc::Documento)
    xml = XMLDocumentNode("1.0")
    drm = addelement!(xml, "DocDRM")

    addelement!(drm, "IdDocto"      , doc.id_docto)
    addelement!(drm, "IdDoctoVersao", doc.id_docto_versao)
    addelement!(drm, "DataBase"     , doc.data_base)
    addelement!(drm, "IdInstFinanc" , "$(doc.id_inst_financ)")
    addelement!(drm, "TipoArq"      , "$(doc.tipo_arq)")
    addelement!(drm, "NomeContato"  , doc.nome_contato)
    addelement!(drm, "FoneContato"  , doc.fone_contato)

    if length(doc.ativo) > 0
        ativo = addelement!(drm, "Ativo")
        for item in sorted_keys(doc.ativo)
            _add_item_carteira!(ativo, item, doc.ativo[item])
        end
    end

    if length(doc.passivo) > 0
        passivo = addelement!(drm, "Passivo")
        for item in sorted_keys(doc.passivo)
            _add_item_carteira!(passivo, item, doc.passivo[item])
        end
    end

    if length(doc.derivativo) > 0
        derivativo = addelement!(drm, "Derivativo")
        for item in sorted_keys(doc.derivativo)
            _add_item_carteira!(derivativo, item, doc.derivativo[item])
        end
    end

    if length(doc.ativo_fundo) > 0
        ativo_fundo = addelement!(drm, "AtivoFundo")
        for item in sorted_keys(doc.ativo_fundo)
            _add_item_carteira!(ativo_fundo, item, doc.ativo_fundo[item])
        end
    end

    if length(doc.passivo_fundo) > 0
        passivo_fundo = addelement!(drm, "PassivoFundo")
        for item in sorted_keys(doc.passivo_fundo)
            _add_item_carteira!(passivo_fundo, item, doc.passivo_fundo[item])
        end
    end

    if length(doc.derivativo_fundo) > 0
        derivativo_fundo = addelement!(drm, "DerivativoFundo")
        for item in sorted_keys(doc.derivativo_fundo)
            _add_item_carteira!(derivativo_fundo, item, doc.derivativo_fundo[item])
        end
    end

    if length(doc.atividade_financeira) > 0
        atividade_financeira = addelement!(drm, "AtividadeFinanceira")
        for item in sorted_keys(doc.atividade_financeira)
            _add_item_carteira!(atividade_financeira, item, doc.atividade_financeira[item])
        end
    end

    return xml
end

function sorted_keys(d::Dict{ItemCarteira, Fluxos})
    v_keys = sort(collect(keys(d)), by = x -> x.item)
    return v_keys
end

function must_be_written(fluxos::Fluxos)
    # write if any vertice must be written
    must = false
    for codigo in sorted_keys(fluxos)
        vertice = fluxos[codigo]
        must = must || must_be_written(vertice)
    end
    return must
end

function _add_item_carteira!(node::EzXML.Node, item::ItemCarteira, fluxos::Fluxos)
    if !must_be_written(fluxos)
        return nothing
    end
    node_item = addelement!(node, "ItemCarteira")
    # item::Symbol
    # id_posicao::SymbolOrNothing
    # fator_risco::Symbol
    # local_registro::SymbolOrNothing
    # carteira_negoc::Symbol
    # fluxos::Vector{FluxoVertice}
    link!(node_item, AttributeNode("Item", _encode_item(item.item)))
    if item.id_posicao != nothing
        link!(node_item, AttributeNode("IdPosicao", _encode_id_posicao(item.id_posicao)))
    end
    link!(node_item, AttributeNode("FatorRisco", _encode_fator_risco(item.fator_risco)))
    if item.local_registro != nothing
        link!(node_item, AttributeNode("LocalRegistro", _encode_local_registro(item.local_registro)))
    end
    link!(node_item, AttributeNode("CarteiraNegoc", _encode_carteira_negoc(item.carteira_negoc)))
    for codigo in sorted_keys(fluxos)
        _add_fluxo_vertice!(node_item, codigo, fluxos[codigo])
    end
    return node
end

function must_be_written(vertice::Vertice)
    valor = _trunc_to_thousands(vertice.valor_alocado)
    return (valor > 0)
end
function _add_fluxo_vertice!(node::EzXML.Node, codigo_vertice::Symbol, vertice::Vertice)
    # Validate if must be written
    if !must_be_written(vertice)
        return nothing
    end
    valor = _trunc_to_thousands(vertice.valor_alocado)
    valor_mam = _trunc_to_thousands(vertice.valor_mam)
    node_fv = addelement!(node, "FluxoVertice")
    link!(node_fv, AttributeNode("CodVertice", "$codigo_vertice"))
    link!(node_fv, AttributeNode("ValorAlocado", "$valor"))
    valor_mam > 0 && link!(node_fv, AttributeNode("ValorMaM", "$valor_mam"))
    return node
end

_encode_item(x::Symbol)::String = "$x"
_encode_id_posicao(x::Symbol)::String = "$x"
_encode_fator_risco(x::Symbol)::String = "$x"
function _encode_local_registro(x::Symbol)::String
    @assert haskey(DECODE_LOCAIS_REGISTRO, x) "local registro invalido: $x"
    return DECODE_LOCAIS_REGISTRO[x]
end
function _encode_carteira_negoc(x::Symbol)::String
    @assert haskey(DECODE_CARTEIRAS, x) "carteira negoc invalida: $x"
    return DECODE_CARTEIRAS[x]
end

_trunc_to_thousands(x::Float64)::Int64 = Int64(trunc(x/1000))

write_xml(io::IO, doc::Documento) = prettyprint(io, to_xml(doc))
write_xml(path::String, doc::Documento) = open(path, "w+") do io write_xml(io, doc) end


text:
#=
maths:
- Julia version: 1.4.2
- Author: akhil
- Date: 2020-05-28
=#
module Maths
    include("Absolute.jl")
    include("Negative.jl")
end

ext: jl


text:
using Test

include("../src/boardrobots.jl")
include("../src/solver.jl")
include("../src/readchallenges.jl")

@testset "Instances" begin
    for file in readdir((@__DIR__)*"/../instances/")
        @test !isempty(find_solution(read_challenge((@__DIR__)*"/../instances/"*file, autodetectformat=true)))
    end
end

@testset "No solution" begin
    @test isempty(find_solution(read_board("""12\n\n\n\n""")))
    @test isempty(find_solution(read_board("""1 2 3\n\n\n\n""")))
end

max_stars_repo_name: remi-garcia/LunarLockout

max_issues_repo_name: remi-garcia/LunarLockout

max_forks_repo_name: remi-garcia/LunarLockout


text:
"""
    function get_thermal_constprice(args...; kwargs...)

Identify generators which cost function is a constant price, locate their buses and returns
the generators, location in the vector of generators, and name. Skips Slack generators,
generators which Active Power is binding and generators with non zero slope. Note: For a
quadratic cost function of the form: αPg²+ βPg + γ, a constant price generator will have the
parameter α == 0.

# Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)
- `res::PowerSimulations.
    OperationsProblemResults`:              Results of the solved OPF for the system
                                             (from PowerSimulations.jl)

# Keywords
- `dual_gen_tol::Float64 = 1e-1`:           Tolerance to identify any binding generators
"""
function get_thermal_constprice(
    sys::System,
    res::PowerSimulations.OperationsProblemResults,
    dual_gen_tol::Float64
    )

    # Define elements
    gen_loc = 0
    gen_bindPg = 0
    gen_constprice = 0
    gens_thermal_constantprice = Array{Any}(undef, 0)
    gens_thermal_constantprice_id = Array{Int64}(undef, 0)
    gens_thermal_constantprice_name = Array{String}(undef, 0)

    # Get all thermal Generators
    gens_thermal = get_components(ThermalStandard, sys)

    # Get dual variables from the OPF results
    dualPgub = get_duals(res)[:P_ub__ThermalStandard__RangeConstraint]
    dualPglb = get_duals(res)[:P_lb__ThermalStandard__RangeConstraint]

    # Identify Generators and save their information
    for gen_thermal in gens_thermal
        gen_loc = gen_loc +1
        if get_bustype(gen_thermal.bus) ≠ BusTypes.REF  #skips slack
            bindub = !isapprox(dualPgub[1, gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            bindlb = !isapprox(dualPglb[1, gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            if !bindub && !bindlb #skips binding Pg gens Dual = 0
                α, β = get_cost(get_variable(get_operation_cost(gen_thermal)))
                γ = get_fixed(get_operation_cost(gen_thermal))
                if α == 0 && (β ≠ 0 || γ ≠ 0)#skips non zero slope
                    gen_constprice = gen_constprice +1
                    resize!(gens_thermal_constantprice, gen_constprice)
                    resize!(gens_thermal_constantprice_id, gen_constprice)
                    resize!(gens_thermal_constantprice_name, gen_constprice)
                    gens_thermal_constantprice[gen_constprice] =
                        get_component(ThermalStandard, sys, gen_thermal.name)
                    gens_thermal_constantprice_id[gen_constprice] = gen_loc
                    gens_thermal_constantprice_name[gen_constprice] = gen_thermal.name
                end
            end
        end
    end

    return (
        gens_thermal_constantprice,
        gens_thermal_constantprice_id,
        gens_thermal_constantprice_name
    )
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
# This is used on Julia version that have the Base.Grisu module

import Base.Grisu

function grisu(v::AbstractFloat, mode, requested_digits)
    return tuple(Grisu.grisu(v, mode, requested_digits)..., Grisu.DIGITS)
end

function plain_precision_heuristic(xs::AbstractArray{<:AbstractFloat})
    ys = filter(isfinite, xs)
    precision = 0
    for y in ys
        len, point, neg, digits = grisu(convert(Float32, y), Grisu.SHORTEST, 0)
        precision = max(precision, len - point)
    end
    return max(precision, 0)
end

# Print a floating point number at fixed precision. Pretty much equivalent to
# @sprintf("%0.$(precision)f", x), without the macro issues.
function format_fixed(x::AbstractFloat, precision::Integer)
    @assert precision >= 0

    if x == Inf
        return "∞"
    elseif x == -Inf
        return "-∞"
    elseif isnan(x)
        return "NaN"
    end

    len, point, neg, digits = grisu(x, Grisu.FIXED, precision)

    buf = IOBuffer()
    if x < 0
        print(buf, '-')
    end

    for c in digits[1:min(point, len)]
        print(buf, convert(Char, c))
    end

    if point > len
        for _ in len:point-1
            print(buf, '0')
        end
    elseif point < len
        if point <= 0
            print(buf, '0')
        end
        print(buf, '.')
        if point < 0
            for _ in 1:-point
                print(buf, '0')
            end
            for c in digits[1:len]
                print(buf, convert(Char, c))
            end
        else
            for c in digits[point+1:len]
                print(buf, convert(Char, c))
            end
        end
    end

    trailing_zeros = precision - max(0, len - point)
    if trailing_zeros > 0 && point >= len
        print(buf, '.')
    end

    for _ in 1:trailing_zeros
        print(buf, '0')
    end

    String(take!(buf))
end

# Print a floating point number in scientific notation at fixed precision. Sort of equivalent
# to @sprintf("%0.$(precision)e", x), but prettier printing.
function format_fixed_scientific(x::AbstractFloat, precision::Integer,
                                 engineering::Bool)
    if x == 0.0
        return "0"
    elseif x == Inf
        return "∞"
    elseif x == -Inf
        return "-∞"
    elseif isnan(x)
        return "NaN"
    end

    mag = floor(Int, log10(abs(x)))
    grisu_precision = precision

    len, point, neg, digits = grisu((x / 10.0^mag), Grisu.FIXED, grisu_precision)
    point += mag

    @assert len > 0

    buf = IOBuffer()
    if x < 0
        print(buf, '-')
    end

    print(buf, convert(Char, digits[1]))
    nextdigit = 2
    if engineering
        while (point - 1) % 3 != 0
            if nextdigit <= len
                print(buf, convert(Char, digits[nextdigit]))
            else
                print(buf, '0')
            end
            nextdigit += 1
            point -= 1
        end
    end

    if precision > 1
        print(buf, '.')
    end

    for i in nextdigit:len
        print(buf, convert(Char, digits[i]))
    end

    for i in (len+1):precision
        print(buf, '0')
    end

    print(buf, "×10")
    for c in string(point - 1)
        if '0' <= c <= '9'
            print(buf, superscript_numerals[c - '0' + 1])
        elseif c == '-'
            print(buf, '⁻')
        end
    end

    return String(take!(buf))
end


text:
remove_spaces(s::AbstractString) = replace(s, r"\s+" => "")

parse_int(s::AbstractString, default = 0) = let m = match(r"\d+", s)
    isnothing(m) ? default : parse(Int64, m.match)
end


text:
using MLStyle

include("MQuery.ConstantNames.jl")
include("MQuery.DynamicInfer.jl")
include("MQuery.Interfaces.jl")
include("MQuery.MacroProcessor.jl")
include("MQuery.Impl.jl")


text:
# GPU CuArray
function gpu_gen_carrier_replica!(
    carrier_replica::CuArray{ComplexF32},
    carrier_frequency,
    sampling_frequency,
    start_phase,
    carrier_amplitude_power,
    start_sample,
    num_samples
)
    @. @views carrier_replica[start_sample:num_samples] = cis(2pi * (1:num_samples) * carrier_frequency / sampling_frequency + start_phase)
end

# GPU StructArray of CuArrays
function gpu_gen_carrier_replica!(
    carrier_replica::StructOfCuVectors,
    carrier_frequency,
    sampling_frequency,
    start_phase,
    carrier_amplitude_power,
    start_sample,
    num_samples
)where{
    StructOfCuVectors <: StructArray{Complex{Float32},1,NamedTuple{(:re, :im),Tuple{CuArray{Float32,1},CuArray{Float32,1}}},Int64}
}
    @. carrier_replica.re = 2pi * (1:num_samples) * carrier_frequency / sampling_frequency + start_phase
    @. carrier_replica.im = sin(carrier_replica.re)
    @. carrier_replica.re = cos(carrier_replica.re)
    return carrier_replica
end

# CPU StructArray
function cpu_gen_carrier_replica!(
    carrier_replica::StructOfVectors,
    carrier_frequency,
    sampling_frequency,
    start_phase,
    carrier_amplitude_power,
    start_sample,
    num_samples
) where {
    StructOfVectors <: StructArray{Complex{Float32},1,NamedTuple{(:re, :im),Tuple{Array{Float32,1},Array{Float32,1}}},Int64}
}
    @. carrier_replica.re = 2pi * (1:num_samples) * carrier_frequency / sampling_frequency + start_phase
    @avx unroll = 6 for i = start_sample:num_samples + start_sample - 1
        carrier_replica.im[i] = sin(carrier_replica.re[i])
        carrier_replica.re[i] = cos(carrier_replica.re[i])
    end
    return carrier_replica
end


text:
"""
    @tspawnat tid -> task
Mimics `Base.Threads.@spawn`, but assigns the task to thread `tid` (with `sticky = true`).
# Example
```julia
julia> t = @tspawnat 4 Threads.threadid()
Task (runnable) @0x0000000010743c70
julia> fetch(t)
4
```
"""
macro tspawnat(thrdid, expr)
    # Copied from ThreadPools.jl with the change task.sticky = false -> true
    # https://github.com/tro3/ThreadPools.jl/blob/c2c99a260277c918e2a9289819106dd38625f418/src/macros.jl#L244
    letargs = Base._lift_one_interp!(expr)

    thunk = esc(:(() -> ($expr)))
    var = esc(Base.sync_varname)
    tid = esc(thrdid)
    quote
        if $tid < 1 || $tid > Threads.nthreads()
            throw(
                AssertionError(
                    "@tspawnat thread assignment ($($tid)) must be between 1 and Threads.nthreads() (1:$(Threads.nthreads()))",
                ),
            )
        end
        let $(letargs...)
            local task = Task($thunk)
            task.sticky = true
            ccall(:jl_set_task_tid, Cvoid, (Any, Cint), task, $tid - 1)
            if $(Expr(:islocal, var))
                put!($var, task)
            end
            schedule(task)
            task
        end
    end
end

"""
# Examples
```julia
interweave([1,2,3,4], [5,6,7,8]) == [1,5,2,6,3,7,4,8]
```
```julia
interweave(1:4, 5:8, 9:12) == [1, 5, 9, 2, 6, 10, 3, 7, 11, 4, 8, 12]
```
"""
function interweave(arrays::AbstractVector{T}...) where {T}
    # check input args
    narrays = length(arrays)
    narrays > 0 || throw(ArgumentError("No input arguments provided."))
    len = length(first(arrays))
    for a in arrays
        length(a) == len || throw(ArgumentError("Only same length inputs supported."))
    end
    # interweave
    res = zeros(T, len * narrays)
    c = 1
    for i in eachindex(first(arrays))
        for a in arrays
            @inbounds res[c] = a[i]
            c += 1
        end
    end
    return res
end

function maybe_gather_sysinfo(lscpustr=nothing; force=false, verbose=false)
    if !SYSINFO_ATTEMPT[] || force
        SYSINFO_ATTEMPT[] = true
        sysinfo = gather_sysinfo_lscpu(lscpustr; verbose)
        if isnothing(sysinfo)
            SYSINFO_SUCCESS[] = false # redundant
            @warn(
                "Couldn't gather system information via `lscpu` (might not be available?). Some features won't work optimally, others might not work at all."
            )
            SYSINFO[] = SysInfo() # default fallback
        else
            SYSINFO_SUCCESS[] = true
            SYSINFO[] = sysinfo
        end
    end
    return nothing
end

function gather_sysinfo_lscpu(lscpustr=nothing; verbose=false)
    local table
    if isnothing(lscpustr)
        try
            table = readdlm(IOBuffer(read(`lscpu --all --extended`, String)))
        catch
            return nothing
        end
    else
        # for debugging purposes
        table = readdlm(IOBuffer(lscpustr))
    end
    colid_socket = findfirst(isequal("SOCKET"), table[1, :])
    colid_numa = findfirst(isequal("NODE"), table[1, :])
    colid_core = findfirst(isequal("CORE"), table[1, :])
    colid_cpu = findfirst(isequal("CPU"), table[1, :])
    colid_online = findfirst(isequal("ONLINE"), table[1, :])
    online_cpu_tblidcs = findall(isequal("yes"), @view(table[:, colid_online]))
    verbose && @show online_cpu_tblidcs
    if length(online_cpu_tblidcs) != Sys.CPU_THREADS
        @warn(
            "Could read `lscpu --all --extended` but number of online CPUs ($(length(online_cpu_tblidcs))) doesn't match Sys.CPU_THREADS ($(Sys.CPU_THREADS))."
        )
    end
    cpuids = table[online_cpu_tblidcs, colid_cpu]
    verbose && @show cpuids
    # hyperthreading?
    hyperthreading = hasduplicates(@view(table[online_cpu_tblidcs, colid_core]))
    verbose && @show hyperthreading
    # count number of sockets
    nsockets = if isnothing(colid_socket)
        1
    else
        length(unique(@view(table[online_cpu_tblidcs, colid_socket])))
    end
    verbose && @show nsockets
    # count number of numa nodes
    nnuma = if isnothing(colid_numa)
        1
    else
        length(unique(@view(table[online_cpu_tblidcs, colid_numa])))
    end
    verbose && @show nnuma
    # cpuids per socket / numa
    cpuids_sockets = [Int[] for _ in 1:nsockets]
    cpuids_numa = [Int[] for _ in 1:nnuma]
    prev_numa = 0
    prev_socket = 0
    numaidcs = Dict{Int,Int}(0 => 1)
    socketidcs = Dict{Int,Int}(0 => 1)
    for i in online_cpu_tblidcs
        cpuid = table[i, colid_cpu]
        numa = isnothing(colid_numa) ? 0 : table[i, colid_numa]
        socket = isnothing(colid_socket) ? 0 : table[i, colid_socket]
        if numa != prev_numa
            if !haskey(numaidcs, numa)
                numaidcs[numa] = maximum(values(numaidcs)) + 1
            end
            prev_numa = numa
        end
        if socket != prev_socket
            if !haskey(socketidcs, socket)
                socketidcs[socket] = maximum(values(socketidcs)) + 1
            end
            prev_socket = socket
        end
        numaidx = numaidcs[numa]
        socketidx = socketidcs[socket]
        push!(cpuids_sockets[socketidx], cpuid)
        push!(cpuids_numa[numaidx], cpuid)
    end
    # if a coreid is seen for a second time
    # the corresponding cpuid is identified
    # as a hypterthread
    ishyperthread = fill(false, length(online_cpu_tblidcs))
    seen_coreids = Set{Int}()
    for i in online_cpu_tblidcs
        cpuid = table[i, colid_cpu]
        coreid = table[i, colid_core]
        if coreid in seen_coreids
            # mark as hyperthread
            cpuidx = findfirst(isequal(cpuid), cpuids)
            ishyperthread[cpuidx] = true
        end
        push!(seen_coreids, coreid)
    end
    return SysInfo(
        nsockets, nnuma, hyperthreading, cpuids, cpuids_sockets, cpuids_numa, ishyperthread
    )
end

hasduplicates(xs::AbstractVector) = length(xs) != length(Set(xs))

function lscpu()
    run(`lscpu --all --extended`)
    return nothing
end


text:
using FiniteDiff: finite_difference_gradient
using BifurcationInference, Test
using Flux: Optimise

error_tolerance = 0.004
@testset "Gradients" begin

	@testset "Saddle Node" begin
		include("minimal/saddle-node.jl")

		@time for θ ∈ [[5.0,1.0],[5.0,-1.0],[-1.0,-1.0],[2.5,-1.0],[-1.0,2.5],[-4.0,1.0]]
			x, y = finite_difference_gradient(θ->loss(F,θ,X),θ), ∇loss(F,θ,X)[2]
			@test acos(x'y/(norm(x)*norm(y)))/π < error_tolerance
		end
	end

	@testset "Pitchfork" begin
		include("minimal/pitchfork.jl")

		@time for θ ∈ [[1.0,1.0],[3.0,3.0],[1.0,-1.0],[3.0,-3.0],[-1.0,-1.0],[-3.0,-3.0],[-1.0,1.0],[-3.0,3.0]]
			x, y = finite_difference_gradient(θ->loss(F,θ,X),θ), ∇loss(F,θ,X)[2]
			@test acos(x'y/(norm(x)*norm(y)))/π < error_tolerance
		end
	end

	@testset "Two State" begin
		include("applied/two-state.jl")

		@time for θ ∈ [[0.5, 0.5, 0.5470, 2.0, 7.5],[0.5, 0.5, 1.5470, 2.0, 7.5],[0.5, 0.88, 1.5470, 2.0, 7.5]]
			x, y = finite_difference_gradient(θ->loss(F,θ,X),θ), ∇loss(F,θ,X)[2]
			@test acos(x'y/(norm(x)*norm(y)))/π < error_tolerance
		end
	end
end

error_tolerance = 0.1
@testset "Optimisations" begin
include("applied/two-state.jl")

	@testset "Mutual Inhibition" begin
		parameters = ( θ=SizedVector{5}(0.5,0.5,0.5470,1.0,1.5), p=minimum(X.parameter) )
		trajectory = train!( F, parameters, X;  iter=300, optimiser=Optimise.ADAM(0.01) )

		steady_states = deflationContinuation(F,X.roots,(p=minimum(X.parameter),θ=trajectory[end]),getParameters(X))
		bifurcations = [ s.z for branch ∈ steady_states for s ∈ branch if s.bif ]

		@test all( z->any(isapprox.(X.targets,z.p; atol=error_tolerance)), bifurcations)
	end

	@testset "Mutual Activation" begin
		parameters = ( θ=SizedVector{5}(0.246,-1.113,-0.87,-1.36,1.28), p=minimum(X.parameter) )
		trajectory = train!( F, parameters, X;  iter=300, optimiser=Optimise.ADAM(0.01) )

		steady_states = deflationContinuation(F,X.roots,(p=minimum(X.parameter),θ=trajectory[end]),getParameters(X))
		bifurcations = [ s.z for branch ∈ steady_states for s ∈ branch if s.bif ]

		@test all( z->any(isapprox.(X.targets,z.p; atol=error_tolerance)), bifurcations)
	end

	include("applied/double-exclusive.jl")
	@testset "Double Exclusive" begin
		parameters = ( θ=SizedVector{21}(-0.5498107083383627, -0.30983502459648393, 0.631939862866151, -1.3407390261868093, 0.6331343612914606, -0.4049028383285403, -0.8524704190293824, -0.11432918583984555, -1.7299987854365662, -2.079032179004907, 0.4916870522978877, -0.010788745553567142, -0.3878513750140796, -0.08334187371052072, -0.2772816729274815, 0.8285363415116678, 0.054887187754647675, -0.49494957503215126, -0.8693294194842291, 0.9289938806228071, 0.22800165770588346), p=minimum(X.parameter) )
		trajectory = train!( F, parameters, X;  iter=300, optimiser=Optimise.ADAM(0.01) )

		steady_states = deflationContinuation(F,X.roots,(p=minimum(X.parameter),θ=trajectory[end]),getParameters(X))
		bifurcations = [ s.z for branch ∈ steady_states for s ∈ branch if s.bif ]

		@test all( z->any(isapprox.(X.targets,z.p; atol=error_tolerance)), bifurcations)
	end
end


text:
using StaticArrays
using IntervalArithmetic
import IntervalArithmetic: interval

import Base: +, *, ^, -

"""
Affine form with center `c`, affine components `γ` and error `Δ`.
"""
struct AF1{N,T<:AbstractFloat}
    c::T   # mid-point
    γ::SVector{N,T}  # affine terms
    Δ::T   # error term.  Error is Δ.(-1..1)
end

function Base.show(io::IO, C::AF1{N,T}) where {N,T}
    print(io, "⟨", C.c, "; ", C.γ, "; ", C.Δ, "⟩")
end

# ==(C::Affine, D::Affine) = C.c == D.c && C.γ == D.γ

"""
Make an `AF1` based on an interval, which is number `i` of `n` total variables.
"""
function AF1(X::Interval{T}, n, i) where {T}
    c = mid(X)
    r = radius(X)

    γ = SVector(ntuple(j->i==j ? r : zero(r), n))

    return AF1(c, γ, T(0))
end

+(x::AF1{N,T}, y::AF1{N,T}) where {N,T} = AF1(x.c + y.c, x.γ .+ y.γ, x.Δ + y.Δ)

-(x::AF1{N,T}, y::AF1{N,T}) where {N,T} = AF1(x.c - y.c, x.γ .- y.γ, x.Δ - y.Δ)

interval(C::AF1) = C.c + sum(abs.(C.γ))*(-1..1) + (C.Δ) * (-1..1)

function *(x::AF1{N,T}, y::AF1{N,T}) where {N,T}
    c = x.c * y.c

    γ = x.c .* y.γ + y.c .* x.γ

    Δ = abs(x.c) * y.Δ + abs(y.c) * x.Δ
    Δ += ( sum(abs.(x.γ)) + abs(x.Δ) ) * ( sum(abs.(y.γ)) + abs(y.Δ) )

    return AF1(c, γ, Δ)

end

*(x::AF1, α::Real) = AF1(α*x.c, α.*x.γ, α*x.Δ)
*(α::Real, x::AF1) = x * α

+(x::AF1, α::Real) = AF1(α+x.c, x.γ, x.Δ)
+(α::Real, x::AF1) = x + α

-(x::AF1) = AF1(-x.c, .-(x.γ), -x.Δ)
-(x::AF1, α::Real) = AF1(x.c - α, x.γ, x.Δ)
-(α::Real, x::AF1) = α + (-x)

^(x::AF1, n::Integer) = Base.power_by_squaring(x, n)

Base.literal_pow(::typeof(^), x::AF1, ::Val{p}) where {T,p} = x^p

x = AF1(3..5, 2, 1)
y = AF1(2..4, 2, 2)

3-x
interval(3-x)

x + y

interval(x+y)

x * y
interval(x * y)

interval(x * y)
interval(x) * interval(y)

z = AF1(-1..1, 1, 1)
z^2
interval(z^2)

using Polynomials

p = Poly([-3, 1])
p2 = p^8

x = 4 ± 1e-4
y = AF1(x, 1, 1)

interval(y)
interval(p2(x))
interval(p2(y))

@time interval(p2(y))

f( (x, y) ) = [x^3 + y, (x - y)^2]

X = IntervalBox(-1..1, -1..1)

f(X)

xx = AF1(X[1], 2, 1)
yy = AF1(X[2], 2, 2)

interval.(f((xx, yy)))

f(X)

x = AF1(4..6, 1, 1)    # example from Messine
f(x) = x * (10 - x)

f(x)
interval(f(x))

interval(10*x - x^2)

function Base.inv(x::AF1)

    range = interval(x)
    a, b = range.lo, range.hi

    # assumes b > a > 0:
    # min-range approx:

    p = -1 / (b^2)
    q = -p * (a + b)^2 / (2a)
    δ = - p * (a- b)^2 / (2a)

    return AF1(p * x.c + q, p .* x.γ, x.Δ)

end


text:
using Distributions
import Random: rand!, rand
import Distributions: logpdf

struct MALA{T}
    ϵ::T
    MALA(ϵ::T) where T = new{T}(ϵ)
end

function rand!(mala::MALA, θ, ∇ll, ::UpdtIdx) where UpdtIdx
    id = [idx(UpdtIdx())...]
    θ[id] .+= mala.ϵ .* ∇ll .+ sqrt.(2.0.*mala.ϵ).*randn(length(∇ll))
    θ
end

function rand(mala::MALA, θ, ∇ll, ::UpdtIdx) where UpdtIdx
    θc = copy(θ)
    rand!(mala, θc, ∇ll, UpdtIdx())
end

function logpdf(mala::MALA, θ, θᵒ, ∇ll, ::UpdtIdx) where UpdtIdx
    id = [idx(UpdtIdx())...]
    logpdf(MvNormal(θ[id] .+ mala.ϵ.*∇ll, diagm(0=>sqrt.(2.0.*mala.ϵ))), θᵒ[id])
end


text:
using Test
using Forecast
import DataFrames: DataFrame

@testset "co2" begin
    co2_ta = co2()
    co2_df = co2(true)
    @test co2_ta isa TimeArray
    @test co2_df isa DataFrame
    @test size(co2_ta) == (4609,)
    @test size(co2_df) == (17166, 17)
end


text:
module FWI

using jInv.Mesh
using jInv.Utils
using jInv.LinearSolvers
using jInv.InverseSolve
using EikonalInv
using KrylovMethods
using Multigrid
using ForwardHelmholtz
using MAT

import jInv.ForwardShare.getData
import jInv.ForwardShare.getSensTMatVec
import jInv.ForwardShare.getSensMatVec
import jInv.LinearSolvers.copySolver

import jInv.ForwardShare.ForwardProbType

FieldsType = Complex32

useSommerfeldBC = true;

fieldsFilenamePrefix = "tempFWIfields"

function getFieldsFileName(omega::Float64)
	omRound = string(round((omega/(2*pi))*100.0)/100.0);
	tfilename = string(fieldsFilenamePrefix,"_f",omRound,"_worker",myid(),".mat");
end

export clear!
export setSourceSelection,setSourceSelectionRatio,setSourceSelectionNum

export FWIparam, getFWIparam
type FWIparam <: ForwardProbType
    omega					:: Float64     # frequencies
	WaveletCoef				:: Complex128
    gamma					:: Vector{Float64}     # attenuation
    Sources					:: Union{Vector{Float64},SparseMatrixCSC,Array{Float64,2}}   # Sources
    Receivers				:: Union{Vector{Float64},SparseMatrixCSC,Array{Float64,2}}
	Fields					:: Array{FieldsType}
    Mesh      				:: RegularMesh
	ForwardSolver			:: AbstractSolver
	forwardSolveBatchSize	:: Int64
	sourceSelection			:: Array{Int64,1}
	useFilesForFields		:: Bool
end

function getFWIparam(omega::Float64, WaveletCoef::Complex128, gamma::Vector{Float64},
							Sources::Union{Vector{Float64},SparseMatrixCSC,Array{Float64,2}},
							Receivers::Union{Vector{Float64},SparseMatrixCSC,Array{Float64,2}},
							Mesh::RegularMesh, ForwardSolver:: AbstractSolver, workerList::Array{Int64},forwardSolveBatchSize::Int64=size(Sources,2),useFilesForFields::Bool = false)
	return getFWIparam([omega], [WaveletCoef],gamma,Sources,Receivers, Mesh,ForwardSolver, workerList,forwardSolveBatchSize,useFilesForFields);
end

function getFWIparam(omega::Array{Float64}, WaveletCoef::Array{Complex128},gamma::Vector{Float64},
							Sources::Union{Vector{Float64},SparseMatrixCSC,Array{Float64,2}},
							Receivers::Union{Vector{Float64},SparseMatrixCSC,Array{Float64,2}},
							Mesh::RegularMesh, ForwardSolver::AbstractSolver, workerList::Array{Int64},forwardSolveBatchSize::Int64=size(Sources,2),useFilesForFields::Bool = false)

	continuationDivision = zeros(Int64,length(omega)+1);
	continuationDivision[1] = 1;

	if workerList==[]
		ActualWorkers = workers();
	else
		ActualWorkers = intersect(workerList,workers());
		if length(ActualWorkers)<length(workerList)
			warn("FWI: workerList included indices of non-existing workers.")
		end
	end
	numWorkers = length(ActualWorkers);
	pFor   = Array{RemoteChannel}(numWorkers*length(omega));
	SourcesSubInd = Array{Array{Int64,1}}(numWorkers*length(omega));
	for k=1:length(omega)
		getFWIparamInternal(omega[k],WaveletCoef[k], gamma,Sources,Receivers,zeros(FieldsType,0), Mesh,
									ForwardSolver, forwardSolveBatchSize ,ActualWorkers,pFor,(k-1)*numWorkers+1,SourcesSubInd,useFilesForFields);
		continuationDivision[k+1] = k*numWorkers+1;
	end
	return pFor,continuationDivision,SourcesSubInd # Array of Remote Refs
end

function getFWIparamInternal(omega::Float64, WaveletCoef::Complex128,gamma::Vector{Float64},
							Sources::Union{Vector{Float64},SparseMatrixCSC,Array{Float64,2}},
							Receivers::Union{Vector{Float64},SparseMatrixCSC,Array{Float64,2}},
							fields::Array{FieldsType}, Mesh::RegularMesh,
							ForwardSolver:: AbstractSolver, forwardSolveBatchSize::Int64,
							Workers::Array{Int64}, pFor::Array{RemoteChannel},startPF::Int64,SourcesSubInd::Array{Array{Int64,1},1},useFilesForFields::Bool = false)
	i = startPF; nextidx() = (idx=i; i+=1; idx)
	nsrc  = size(Sources,2);
	numWorkers = length(Workers);
	# send out jobs
	@sync begin
		for p=Workers
			@async begin
				while true
					idx = nextidx()
					if idx > startPF + numWorkers - 1
						break
					end
					I_k = getSourcesIndicesOfKthWorker(numWorkers,idx - startPF + 1,nsrc)
					SourcesSubInd[idx] = I_k;
					pFor[idx] = initRemoteChannel(getFWIparamInternal,p, omega,WaveletCoef,  gamma, Sources[:,I_k], Receivers, fields, Mesh,
																			copySolver(ForwardSolver),forwardSolveBatchSize,useFilesForFields);
					wait(pFor[idx]);
				end
			end
		end
	end
	return pFor # Array of Remote Refs
end

function getFWIparamInternal(omega::Float64,WaveletCoef::Complex128,
							gamma::Vector{Float64},
	                        Sources::Union{Vector{Float64},SparseMatrixCSC,Array{Float64,2}},
							Receivers::Union{Vector{Float64},SparseMatrixCSC,Array{Float64,2}},
							Fields::Array{FieldsType},
							Mesh::RegularMesh, ForwardSolver:: AbstractSolver, forwardSolveBatchSize::Int64,useFilesForFields::Bool = false)
	return FWIparam(omega,WaveletCoef,gamma,Sources,Receivers,Fields,Mesh,ForwardSolver,forwardSolveBatchSize,Array{Int64}(0),useFilesForFields)
end

# function setSourceSelection(pForRF::RemoteRef{Channel{Any}}, selection::Array{Int64,1})
	# s = copy(selection);
	# if minimum(s) < 1
		# error("FWI: sources selection out of range");
	# end
	# pFor  = take!(pForRF);
	# s = [];
	# if isa(pFor,FWIparam)
		# Sources = pFor.Sources
		# if maximum(s)>size(Sources,2)
			# s = s[s.<=size(Sources,2)];
			# warn("FWI: reducing selection: s = s[s.<=size(Sources,2)]");
		# end
		# pFor.sourceSelection = s;
	# end
	# put!(pForRF,pFor);
	# return s;
# end

# function setSourceSelectionRatio(pForRF::RemoteRef{Channel{Any}}, selectionRatio::Float64)
	# if selectionRatio <= 0.0 || selectionRatio > 1.0
		# error("selectionRatio has to be between 0 and 1");
	# end
	# if selectionRatio == 1.0
		# return;
	# end
	# pFor  = take!(pForRF);
	# s = [];
	# if isa(pFor,FWIparam)
		# Q = pFor.Sources
		# s = randperm(size(Q,2))[1:ceil(Int64,size(Q,2)*selectionRatio)];
		# pFor.sourceSelection = s;
	# end
	# put!(pForRF,pFor);
	# return s;
# end

# function setSourceSelectionNum(pForRF::RemoteRef{Channel{Any}}, selectionNum::Int64)
	# if selectionNum <= 0
		# error("selectionNum has to be bigger than 1");
	# end
	# pFor  = take!(pForRF);
	# s = [];
	# if isa(pFor,FWIparam)
		# Q = pFor.Sources
		# s = randperm(size(Q,2))[1:selectionNum];
		# pFor.sourceSelection = s;
	# end
	# put!(pForRF,pFor);
	# return s;
# end

import jInv.Utils.clear!
function clear!(pFor::FWIparam)
	clear!(pFor.ForwardSolver);
	pFor.Fields = zeros(FieldsType,0);
	clear!(pFor.Mesh);
	return pFor;
end

include("getData.jl")
include("getSensMatVec.jl")
include("getSensTMatVec.jl")
include("FourthOrderHesPrec.jl")
include("freqCont.jl")
include("timeDomainFWI.jl")
include("firstArivalPicking.jl")
end


text:
module motor  #jl

using ReachabilityAnalysis, SparseArrays
using ReachabilityModels: @relpath

# system matrix
I = [1, 2, 2, 3, 3, 3, 3, 4, 5, 6, 6, 7, 7, 7, 7, 8]
J = [2, 3, 2, 1, 2, 3, 4, 1, 6, 7, 6, 5, 6, 7, 8, 5]
V = [1, 8487.2, -1.0865, -2592.1, -21.119, -698.91, -141399.0, 1.0, 1.0,
     8487.2, -1.0865, -2592.1, -21.119, -698.91, -141399.0, 1.0]
A = sparse(I, J, V)

# input matrix
B = sparse([4, 8], [1, 2], [-1.0, -1.0])

# state domain
X = Universe(8)

# input domain
U = Hyperrectangle([0.23, 0.3], [0.07, 0.1])

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 Test
import ArchGDAL as AG
import GeoFormatTypes as GFT

@testset "test_convert.jl" begin

    # Tests high level convert methods
    @testset "convert point format" begin
        point = AG.createpoint(100, 70)
        json = convert(GFT.GeoJSON, point)
        @test sprint(print, convert(AG.IGeometry, json)) ==
              "Geometry: POINT (100 70)"
        kml = convert(GFT.KML, point)
        gml = convert(GFT.GML, point)
        wkb = convert(GFT.WellKnownBinary, point)
        wkt = convert(GFT.WellKnownText, point)
        @test json.val == AG.toJSON(point)
        @test kml.val == AG.toKML(point)
        @test gml.val == AG.toGML(point)
        @test wkb.val == AG.toWKB(point)
        @test wkt.val == AG.toWKT(point)
        @test convert(GFT.GeoJSON, json) ==
              convert(GFT.GeoJSON, wkb) ==
              convert(GFT.GeoJSON, wkt) ==
              convert(GFT.GeoJSON, gml) ==
              json
        @test convert(GFT.KML, gml) == convert(GFT.KML, wkt)
    end

    @testset "convert crs format" begin
        proj4326 = GFT.ProjString("+proj=longlat +datum=WGS84 +no_defs")
        @test convert(
            GFT.ProjString,
            GFT.CRS(),
            convert(
                GFT.WellKnownText,
                GFT.CRS(),
                convert(GFT.ESRIWellKnownText, GFT.CRS(), GFT.EPSG(4326)),
            ),
        ) == proj4326
        @test convert(GFT.CoordSys, GFT.CRS(), proj4326) isa GFT.CoordSys
        @test convert(GFT.GML, GFT.CRS(), proj4326) isa GFT.GML
    end

    @testset "geometry conversions" begin
        geom1 = AG.createpoint(1, 2)
        @test typeof(geom1) == AG.IGeometry{AG.wkbPoint}
        geom2 = convert(AG.IGeometry{AG.wkbUnknown}, geom1)
        @test typeof(geom2) == AG.IGeometry{AG.wkbUnknown}
        @test AG.toWKT(geom1) == AG.toWKT(geom2)
    end

    @testset "type conversions" begin
        @test convert(AG.OGRFieldType, Int32) == AG.OFTInteger
        @test convert(AG.OGRFieldType, Int16) == AG.OFTInteger
        @test convert(AG.OGRFieldType, Bool) == AG.OFTInteger

        @test convert(AG.OGRFieldType, Float32) == AG.OFTReal
        @test convert(AG.OGRFieldType, Float64) == AG.OFTReal

        # Reverse conversion should result in default type, not subtype
        @test convert(DataType, AG.OFTInteger) == Int32
        @test convert(DataType, AG.OFTReal) == Float64
    end
end


text:
# Event

# constructor
@test length(Event()) == 0
@test length(Event(:foo => 1)) == 1
@test length(Event(:foo => 1, :bar => 2)) == 2
@test length(Event(Dict(:foo => 1))) == 1

# show
@test length("$(Event(:foo => 1, :bar => 1))") > 0

# merge
@test merge(Event(:foo => 2), Event(:foo => 1)) == Event(:foo => 1)
@test merge(Event(:foo => 1), Dict(:bar => 1)) == Event(:foo => 1, :bar => 1)
@test merge(Event(:foo => 1), :bar => 1) == Event(:foo => 1, :bar => 1)

# repeat
@test repeat(Event(:foo => 1), 2) == [Event(:foo => 1), Event(:foo => 1)]

# assoc
@test assoc(Event(), :foo, 1) == Event(:foo => 1)
@test dissoc(Event(:foo => 1), :foo) == Event()

# conversion
@test convert(Event, :foo => 1) == Event(:foo => 1)


text:
using ApproxFun, SingularIntegralEquations, Test
    import SingularIntegralEquations: ⁻, ⁺

@testset "Wiener–Hopf" begin
    # Scalar case

    G=Fun(z->2+cos(z+1/z),Circle()) # the symbol of the Toeplitz operator
    T=ToeplitzOperator(G)

    C  = Cauchy(-1)
    V  = (I+(1-G)*C)\(G-1)

    Φmi = 1+C*V
    Φp = V+Φmi

    L  = ToeplitzOperator(1/Φmi)
    U  = ToeplitzOperator(Φp)

    @test norm((T-U*L)[1:10,1:10]) < 10eps()  # check the accuracy

    # Matrix case

    G=Fun(z->[-1 -3; -3 -1]/z +
             [ 2  2;  1 -3] +
             [ 2 -1;  1  2]*z,Circle())

    H = Hilbert(space(G),1)
    @test domainspace(H) == rangespace(H) == space(G)
    C  = Cauchy(-1)

    A=(I+(I-G)*C)

    F=(G-I)[:,1]
    F̃ = vec(F)

    C2=promotedomainspace(C.ops[2],space(F))

    testblockbandedoperator(Hilbert(space(F)[1]))
    testblockbandedoperator(Hilbert(space(F)))
    testblockbandedoperator(C2)

    @test norm((C*F - [C*F[1];C*F[2]]).coefficients) == 0
    @test norm((C*G - [C*G[1] C*G[3];C*G[2] C*G[4]]).coefficients) == 0
    @test cauchy(F,exp(0.1im)⁻) ≈ (C*F)(exp(0.1im))
    @test cauchy(G,exp(0.1im)⁻) ≈ (C*G)(exp(0.1im))

    V1 = A\F
    Ṽ1 = A\F̃

    A1=choosespaces(A,(G-I)[:,1])
    A2=choosespaces(A,Fun((G-I)[:,1]))

    @test A1\Fun((G-I)[:,1])  == V1
    @test A1\(G-I)[:,1]  == V1

    QR=qrfact(A1)

    @test QR\Fun((G-I)[:,1]) == V1
    @test QR\(G-I)[:,1] == V1

    @test norm((V1-Ṽ1).coefficients) == 0
    @test norm((A*V1-F).coefficients) < 100eps()

    @test norm((F-Fun((G-I)[:,1])).coefficients) == 0
    @test Fun(V1) == V1

    Ṽ = QR\(G-I)
    V  = (I+(I-G)*C)\(G-I)

    @test map(f->f(exp(0.1im)),Array(G-I)) ≈ (G-I)(exp(0.1im))

    @test (G-I)[:,1]==Fun((G-I)[:,1],rangespace(QR))
    @test (G-I)[:,1]==Fun(vec((G-I)[:,1]),rangespace(QR))

    @test norm((V-Ṽ).coefficients) == 0

    @test norm((V1-V[:,1]).coefficients) == 0

    V2  = A\(G-I)[:,2]
    @test norm((V2-V[:,2]).coefficients) == 0

    @test norm((A*V[:,1]-(G[:,1]-[1,0])).coefficients) < 100eps()

    z=exp(0.1im)
    @test V(z)+(I-G(z))*cauchy(V,(z)⁻) ≈ G(z)-I

    @test cauchy(V[1,1],exp(0.1im)⁻) ≈ (C*V[1,1])(exp(0.1im))
    @test cauchy(V[2,1],exp(0.1im)⁻) ≈ (C*V[2,1])(exp(0.1im))
    @test cauchy(V[:,1],exp(0.1im)⁻) ≈ (C*V[:,1])(exp(0.1im))
    @test cauchy(V,exp(0.1im)⁻) ≈ (C*V)(exp(0.1im))

    Φmi = I+C*V
    Φp = V+Φmi

    @test Φmi(z) ≈ (I+cauchy(V,(z)⁻))
    @test Φp(z) ≈ (I+cauchy(V,(z)⁺))
    @test Φp(z)*inv(Φmi(z)) ≈ G(z)

    Φm=inv.(Φmi)
    @test Φm(z) ≈ inv(Φmi(z))

    T=ToeplitzOperator(G)

    L  = ToeplitzOperator(Φm)
    U  = ToeplitzOperator(Φp)

    @test norm((T-U*L)[1:10,1:10]) < 100eps()  # check the accuracy
end


text:
function impulses{S<:Int}(solution::Perturbable_Soln,impulse_length::S,innovation_to_shock::S)

  if isa(solution,Union(Blanchard_Kahn_Soln,Klein_Soln)) == true
    responses = compute_first_order_state_space_re_impulses(solution,impulse_length,innovation_to_shock)
  elseif isa(solution,State_Space_Soln) == true
    responses = compute_state_space_op_impulses(solution,impulse_length,innovation_to_shock)
  elseif isa(solution,Binder_Pesaran_Soln) == true
    responses = compute_structural_re_impulses(solution,impulse_length,innovation_to_shock)
  elseif isa(solution,Structural_Soln) == true
    responses = compute_structural_op_impulses(solution,impulse_length,innovation_to_shock)
  elseif isa(solution,Second_Order_State_Space_Soln) == true
    if isa(solution,Lombardo_Sutherland_Soln) == true
      solution = convert_second_order(solution)
    end
    responses = compute_second_order_state_space_re_impulses(solution,impulse_length,innovation_to_shock)
  end

end

function compute_first_order_state_space_re_impulses(solution,impulse_length,innovation_to_shock)

  p     = copy(solution.p)
  k     = copy(solution.k)
  h     = copy(solution.h)
  sigma = copy(solution.sigma)

  nx = size(p,1)
  ny = size(h,1)

  s = chol(sigma[:,:])'

  responses = zeros(nx+ny,impulse_length)
  responses[1:nx,1]           = (k*s)[:,innovation_to_shock]
  responses[(nx+1):(nx+ny),1] = h*responses[1:nx,1]

  for i = 2:impulse_length

    responses[1:nx,i]           = p*responses[1:nx,i-1]
    responses[(nx+1):(nx+ny),i] = h*responses[1:nx,i]

  end

  return responses

end

function compute_state_space_op_impulses(solution,impulse_length,innovation_to_shock)

  p     = copy(solution.p)
  k     = copy(solution.k)
  h     = copy(solution.h)
  f     = copy(solution.f)
  sigma = copy(solution.sigma)

  nx = size(p,1)
  ny = size(h,1)
  np = size(f,1)

  s = chol(sigma[:,:])'

  responses = zeros(nx+ny+np,impulse_length)
  responses[1:nx,1]                 = (k*s)[:,innovation_to_shock]
  responses[(nx+1):(nx+ny),1]       = h*responses[1:nx,1]
  responses[(nx+ny+1):(nx+ny+np),1] = f*responses[1:nx,1]

  for i = 2:impulse_length

    responses[1:nx,i]                 = p*responses[1:nx,i-1]
    responses[(nx+1):(nx+ny),i]       = h*responses[1:nx,i]
    responses[(nx+ny+1):(nx+ny+np),i] = f*responses[1:nx,i]

  end

  return responses

end

function compute_structural_re_impulses(solution,impulse_length,innovation_to_shock)

  p     = copy(solution.p)
  k     = copy(solution.k)
  sigma = copy(solution.sigma)

  s = chol(sigma[:,:])'

  n = size(p,1)

  responses = zeros(n,impulse_length)
  responses[1:n,1] = (k*s)[:,innovation_to_shock]

  for i = 2:impulse_length

    responses[1:n,i] = p*responses[1:n,i-1]

  end

  return responses

end

function compute_structural_op_impulses(solution,impulse_length,innovation_to_shock)

  p     = copy(solution.p)
  k     = copy(solution.k)
  sigma = copy(solution.sigma)

  s = chol(sigma[:,:])'

  n = size(p,1)

  responses = zeros(n,impulse_length)
  responses[1:n,1] = (k*s)[:,innovation_to_shock]

  for i = 2:impulse_length

    responses[1:n,i] = p*responses[1:n,i-1]

  end

  return responses

end

function compute_second_order_state_space_re_impulses(solution,impulse_length,innovation_to_shock)

  hx    = copy(solution.hx)
  gx    = copy(solution.gx)
  hxx   = copy(solution.hxx)
  gxx   = copy(solution.gxx)
  eta   = copy(solution.eta)
  sigma = copy(solution.sigma)

  s = chol(sigma[:,:])'

  nx = size(hx,1)
  ny = size(gx,1)

  responses = zeros(nx+ny,impulse_length)
  responses[1:nx,1]           = (eta*convert(Array,s))[:,innovation_to_shock]
  responses[(nx+1):(nx+ny),1] = gx*responses[1:nx,1] + 0.5*(kron(eye(ny),responses[1:nx,1]'))*gxx*responses[1:nx,1]

  for i = 2:impulse_length

    responses[1:nx,i]           = hx*responses[1:nx,i-1] + 0.5*(kron(eye(nx),responses[1:nx,i-1]'))*hxx*responses[1:nx,i-1]
    responses[(nx+1):(nx+ny),i] = gx*responses[1:nx,i]   + 0.5*(kron(eye(ny),responses[1:nx,i-1]'))*gxx*responses[1:nx,i-1]

  end

  return responses

end


text:
function plot_country_posterior(chains, grouped_dat, country; tspan = [0.0, 365.0], subsample_size=100)

    country_pop = grouped_dat[country].pop[1]
    country_name = grouped_dat[country].Country_Region[1]

    post_dat = DataFrame(chains);

    init = [(country_pop - 100.0) / country_pop, 100.0 / country_pop, 0.0]
    country_param_cols = [Symbol("β[" * string(country) * "]"), Symbol("γ[" * string(country) * "]")]

    posterior_mean_params = mean.(eachcol(post_dat[:, country_param_cols]))
    posterior_mean_prob = ODEProblem(sir_ode, init, tspan, posterior_mean_params)
    posterior_mean_sol = solve(posterior_mean_prob);

    # plt = plot(posterior_mean_sol;
    #     title = country_name * "\n",
    #     title_location = :left,
    #     top_margin = 5mm,
    #     vars = [1],
    #     xguidefontsize = 10, yguidefontsize = 10,
    #     xlabel = "Days since 100 confirmed cases",
    #     ylabel = "Proportion of population",
    #     label = "S",
    #     color = "blue",
    #     linewidth = 2,
    #     legend = :outerright)
    plt = plot(posterior_mean_sol;
        title = country_name * "\n",
        title_location = :left,
        top_margin = 5mm,
        vars = [2],
        xguidefontsize = 10, yguidefontsize = 10,
        xlabel = "Days since 100 confirmed cases",
        ylabel = "Proportion of population",
        label = "I",
        color = "red",
        linewidth = 2,
        legend = :outerright)
    # plot!(posterior_mean_sol; label="I", vars = [2], color = "red", linewidth = 2)
    plot!(posterior_mean_sol; label="R", vars = [3], color = "green", linewidth = 2)

    # scatter!(grouped_dat[country].susc_prop, color = "blue", label = "Obs. S", markerstrokewidth = 0)
    scatter!(grouped_dat[country].case_prop, color = "red", label = "Obs. I", markerstrokewidth = 0)
    scatter!(grouped_dat[country].recov_prop, color = "green", label = "Obs. R", markerstrokewidth = 0)

    # subsample to avoid a crowded plot
    sub_idx = sample(axes(post_dat, 1), subsample_size; replace = false);
    for i in sub_idx
        prob = ODEProblem(sir_ode, init, tspan, post_dat[i, country_param_cols])
        sol = solve(prob);
        # plot!(sol; vars = [1], label="", color = "blue", alpha = 0.1, xlabel = "Days since 100 confirmed cases")
        plot!(sol; vars = [2], label="", color = "red", alpha = 0.1, xlabel = "Days since 100 confirmed cases")
        plot!(sol; vars = [3], label="", color = "green", alpha = 0.1, xlabel = "Days since 100 confirmed cases")
    end

    return plt
end


text:
import Random # used for example

# Define the new generative support information struct
struct MyGenInfo <: InfiniteOpt.AbstractGenerativeInfo  # REPLACE WITH ACTUAL NAME
    attr::Int # REPLACE WITH DESIRED ATTRIBUTES
    # ADD MORE ATTRIBUTES AS NEEDED
end

# Make a unique support label (recommended)
struct MyGenLabel <: InternalLabel end

# Extend support_label
function InfiniteOpt.support_label(info::MyGenInfo)::Type{MyGenLabel} # REPLACE WITH ACTUAL INFO TYPE
    return MyGenLabel # REPLACE WITH ACTUAL MAPPING
end

# Extend make_generative_supports
function InfiniteOpt.make_generative_supports(
    info::MyGenInfo,
    pref,
    supps
    )::Vector{Float64} # REPLACE WITH ACTUAL INFO TYPE
    # REPLACE BELOW WITH ACTUAL CODE TO CREATE THE GENERATIVE SUPPORTS BASED ON THE EXISTING
    num_existing = length(supps)
    num_existing <= 1 && error("`$pref` doesn't have enough supports.")
    num_internal = info.attr
    gen_supps = Float64[]
    for i = 1:num_existing-1
        lb = supps[i]
        ub = supps[i+1]
        append!(gen_supps, rand(num_internal) * (ub - lb) .+ lb)
    end
    return gen_supps
end


text:
module MyExample
using ForwardDiff

greet() = print("Hello World")

include("extra_file.jl")
export my_f, derivative_of_my_f
end


text:
struct Graph
    n :: Int32 # |V|
    m :: Int32 # |E|
    V :: Array{Int32, 1} # V[i] = Real Index of node i
    E :: Array{Tuple{Int32, Int32, Int32, Float64}, 1} # (ID, u, v, w) in Edge Set
end


text:
using FastRef
using LazySets
import FastRef: forward_network, forward_affine_map, ishull

center = fill(1.0, 400)
radius = fill(1.0, 400)
inputSet = Hyperrectangle(center, radius)

nnet = read_nnet("nnet/mnist-20.nnet")
solver1 = MaxSens(2.0, true)
solver2 = FastGrid(2.0)

(W, b) = (nnet.layers[1].weights, nnet.layers[1].bias)

#outputSet1 = forward_network(solver1, nnet, inputSet)

z = forward_affine_map(solver2, W, b, inputSet)
#outputSet2 = forward_network(solver2, nnet, z)

c = W * center + b
r = z.radius .* 0.03

z1 = Hyperrectangle(c, r)

print("z: ")
println(z)

#println(c)

#print("z1: ")
#println(z1)

println(ishull(z1, inputSet, W, b))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libblastrampoline_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libblastrampoline")
JLLWrappers.@generate_main_file("libblastrampoline", UUID("8e850b90-86db-534c-a0d3-1478176c7d93"))
end  # module libblastrampoline_jll


text:
using Pkg
pkg"activate ."

include("Database.jl")
include("Wikipedia.jl")
include("Gameplay.jl")
include("GameSession.jl")
include("WebApp.jl")

using .Wikipedia, .Gameplay, .GameSession, .WebApp

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
using Test
using SVGPlots

x = collect(-4:0.1:4)
figure(;width=800, height=600)
subplot(2,3,1)
plot(x, x)
hold(true)
plot(x, x.^2)
subplot(2,3,2)
plot(x, x.^2)
subplot(2,3,3)
plot(x, x.^3)
subplot(2,3,4)
plot(x, x.^4)
subplot(2,3,5)
plot(x, x.^5)
subplot(2,3,6)
plot(x, x.^6)
display(gcf())
figure()
plot(x, x)
hold(true)
plot(x, x.^2, linetype = :dotted)
plot(x, x.^3, linetype = :dashed)
plot(x, x.^4, linetype = :dashdotted)
plot(x, x.^5, linetype = :none)
title("This is a title")
xlabel("x")
ylabel("y")
display(gcf())
figure()
subplot(2,2,1)
plot(x, x)
subplot(2,2,2)
plot(x, x.^2)
subplot(2,2,3)
plot(x, x.^3)
subplot(2,2,4)
plot(x, x.^4)
hold(true)
plot(x, x.^5)
title("test")
display(gcf())

max_stars_repo_name: JuliaTagBot/SVGPlots

max_issues_repo_name: JuliaTagBot/SVGPlots

max_forks_repo_name: JuliaTagBot/SVGPlots


text:
# Chebyshev Polynomial of the First Kind
function Tk(τk::AbstractFloat, k::Integer)
    return cos(k*acos(τk))
end

function T(τ::AbstractVector, N::Integer)
    # M + 1 (Number of nodes)
    Mp1 = length(τ)

    # Initialize Chebyshev Interpolation Matrix
    Tτ = SizedMatrix{N, Mp1}(ones(N, Mp1))

    # Fill Matrix
    Tτ[2, :] .= τ
    @inbounds Threads.@threads for col in 1:Mp1
        for row in 3:N
            Tτ[row,col] = Tk(τ[col], row - 1)
        end
    end
    return Tτ
end

function T!(Tτ::AbstractMatrix, τ::AbstractVector)
    dims = size(Tτ)

    # Fill Matrix
    Tτ[1, :] .= 1
    Tτ[2, :] .= τ
    @inbounds Threads.@threads for col in 1:dims[2]
        for row in 3:dims[1]
            Tτ[row,col] = Tk(τ[col], row - 1)
        end
    end
end

function ChebyshevInterpMat!(Cx::AbstractMatrix, τ::AbstractVector)
    dims = size(Cx)

    # Fill Matrix
    Cx[:,1] .= 1
    Cx[:,2] .= τ
    @inbounds Threads.@threads for col in 3:dims[2]
        for row in 1:dims[1]
            Cx[row,col] = Tk(τ[row], col - 1)
        end
    end
end


text:
import Base: iterate, length, keys, values, haskey, delete!, size, getindex,
             IndexStyle

# MinHashing

"""
    minhash(s, k[, seed])

Split `s` into `k`-shingles and hash the shingles using a hash function seeded
by `seed`.

# Examples
```jldoctest
julia> minhash("Hello, world!", 5)
0x12da77c8

julia> minhash("Hello, world!", 5, 1)
0x2011211e
```
"""
@inline function minhash(s::AbstractString, k::Int, seed::Integer)
    k > 1 || throw(DomainError("`k` must be >= 2"))
    seed >= 0 || throw(DomainError("`seed` must be >= 0"))

    h = typemax(UInt32)

    for i = 1:length(s)-k+1
        hᵢ = murmur32(SubString(s, i:i+k-1), seed)

        if hᵢ < h
            h = hᵢ
        end
    end

    h
end
minhash(s::AbstractString, k::Int) = minhash(s, k, 0)

"""
    signature(s, k, n)

Split `s` into `k`-shingles and hash using `n` hash functions. The ith hash
function will be seeded with i.

# Examples
```jldoctest
julia> signature("Hello, world!", 5, 3)
3-element Array{UInt32,1}:
 0x2011211e
 0x1134f986
 0x291d3c45

julia> signature("Hello, world!", 6, 2)
2-element Array{UInt32,1}:
 0x441d7a91
 0x1af9b8f8
```
"""
function signature(s::AbstractString, k::Int, n::Int)
    k > 1 || throw(DomainError("`k` must be >= 2"))
    n > 0 || throw(DomainError("`n` must be >= 1"))

    A = Vector{UInt32}(undef, n)

    for seed = 1:n
        A[seed] = minhash(s, k, seed)
    end

    A
end

"""
    AbstractSignature{T} <: AbstractVector{T}

Abstract MinHash signature type.

Any subtype `S` <: `AbstractSignature` should implement the following methods:

* `identifier`: return the identifier for the signature
* `signature`: return the vector of MinHashes that constitute the signature
"""
abstract type AbstractSignature{T} <: AbstractVector{T} end

identifier(s::AbstractSignature) = s.id
signature(s::AbstractSignature) = s.signature
size(s::AbstractSignature) = size(signature(s))
getindex(s::AbstractSignature, i::Int) = getindex(signature(s), i)
IndexStyle(::AbstractSignature) = IndexLinear()

"""
    identifier(s::AbstractSignature)

Return the identifier for the signature `s`.

# Examples
```jldoctest
julia> s = Signature("ABC", UInt8[1,2,3]);

julia> identifier(s)
"ABC"
```
"""
identifier

"""
    signature(s::AbstractSignature)

Return the vector of MinHashes that constitute the signature `s`.

# Examples
```jldoctest
julia> s = Signature("ABC", UInt8[1,2,3]);

julia> signature(s)
3-element Array{UInt8,1}:
 0x01
 0x02
 0x03
```
"""
signature

"""
    Signature(id, signature)

Generic signature type.

# Examples
```jldoctest
julia> Signature("A", UInt32[1,2,3])
Signature{String,Array{UInt32,1}}("A", UInt8[0x01, 0x02, 0x03])
```
"""
struct Signature{T<:Unsigned} <: AbstractSignature{T}
    id::String
    signature::Vector{T}
end

"""
    hashband(signature, start, stop)

Hashes a band of `signature` from `start` to `stop` into a bucket.

# Examples
```jldoctest
julia> s = Signature("ID123", UInt32[1,2,3,4]);

julia> hashband(s, 1, 2)
0xc3642e86

julia> hashband(s, 3, 4)
0xe68b97c2
```
"""
@inline function hashband(s::AbstractSignature, start::Int, stop::Int)
    murmur32(s[start:stop])
end

# HashTable

"""
    HashTable(band, hashtable)
    HashTable(K, V, band)

Construct a hashtable.

# Examples
```jldoctest
julia> h = HashTable(String, Int, 1)

julia> h["A"] = 1;

julia> h["A"] = 2;

julia> h
HashTable{String,Int64} with 1 entry:
  "A" => Set([2, 1])
```
"""
struct HashTable{K,V} <: AbstractDict{K,Set{V}}
    band::Int
    hashtable::Dict{K,Set{V}}
end
HashTable(K::Type, V::Type, band::Int) = HashTable(band, Dict{K,Set{V}}())

band(h::HashTable) = h.band
hashtable(h::HashTable) = h.hashtable

for f in (:iterate, :length, :keys, :values)
    @eval ($f)(h::HashTable) = ($f)(hashtable(h))
end
for f in (:haskey, :delete!)
    @eval ($f)(h::HashTable, k) = ($f)(hashtable(h), k)
end
function Base.setindex!(h::HashTable, v, k)
    haskey(h, k) ? push!(hashtable(h)[k], v) : hashtable(h)[k] = Set([v])
end
Base.iterate(h::HashTable, i::Int) = iterate(hashtable(h), i)
Base.get(h::HashTable, k, default) = get(hashtable(h), k, default)

# LSH

"""
    lsh(signature, bands)

Splits a `signature` into `bands` bands which are hashed into buckets.

# Examples
```jldoctest
julia> signatures = [Signature("A", UInt32[1,2,3,4]), Signature("B", UInt32[1,2,5,6])];

julia> lsh(signatures, 2)
2-element Array{LSMH.HashTable{UInt64,Set{String}},1}:
 LSMH.HashTable{UInt64,Set{String}}(1, Dict(0x8fd27f36c41781c8=>Set(["B", "A"])))
 LSMH.HashTable{UInt64,Set{String}}(2, Dict{UInt64,Set{String}}())
```
"""
function lsh(signatures::AbstractVector{<:AbstractSignature}, bands::Int)
    hashtables = [HashTable(UInt64, String, b) for b = 1:bands]

    for s = signatures
        buckets = lsh(s, bands)

        for b = 1:bands
            hashtable = hashtables[b]
            bucket = buckets[b]
            hashtable[bucket] = identifier(s)
        end
    end

    map(filter_collisions!, hashtables)
    hashtables
end
function lsh(S::AbstractSignature, bands::Int)::Vector{<:Unsigned}
    step = div(size(signature(S))[1], bands)
    buckets = Vector{UInt64}(undef, bands)

    for i = 1:bands
        start = (i - 1) * step + 1
        buckets[i] = hashband(S, start, start+step-1)
    end

    buckets
end

"""
    filter_collisions!(d)

Removes non-colliding sequences from `d`.

Often `d<:AbstractDict{K,<:AbstractSet{V}}`.

# Examples
```jldoctest
julia> d = Dict("A" => Set([1,2]), "B" => Set([3]));

julia> filter_collisions!(d)

julia> d
Dict{String,Set{Int64}} with 1 entry:
  "A" => Set([2, 1])
```
"""
function filter_collisions!(d::AbstractDict)
    for k = keys(d)
        if length(d[k]) == 1
            delete!(d, k)
        end
    end
end

# Similarity

"""
    jaccard(a, b)

Calculate Jaccard similarity of two sets.

# Examples
```jldoctest
julia> jaccard(1:100, 2:100)
0.99

julia> jaccard([], [])
1.0
```
"""
function jaccard(a::T, b::T) where T<:Union{AbstractVector,AbstractSet}
    length(a) == 0 && length(b) == 0 && return 1.
    a = Set(a)
    b = Set(b)
    length(a ∩ b) / length(a ∪ b)
end

"""
    candidates(d)

Return pairs of candidate matches from `d`.

# Examples
```jldoctest
julia> d = Dict(nothing => Set([1,2,3]));

julia> candidates(d)
3-element Array{Set{Int64},1}:
 Set([2, 3])
 Set([2, 1])
 Set([3, 1])
```
"""
function candidates(d::AbstractDict)
    v = Vector{Set{eltype(eltype(values(hashtable)))}}()
    sizehint!(v, length(d))

    for (bucket, IDs) = d
        for (i, j) = combinations(collect(IDs), 2)
            push!(v, Set([i,j]))
        end
    end

    unique!(v)
end

# """
#     filter_candidates(signatures, candidates, threshold)
#
# Remove `candidates` with Jaccard similarity between their `signatures` less than
# `threshold`.
# """
# function filter_candidates(signatures, candidates, threshold)
#     out = Vector{Tuple{Int,Int}}()
#
#     for (i, j) = candidates
#         println(i,j)
#         similarity = jaccard(signatures[:,i], signatures[:,j])
#         similarity >= threshold ? push!(out, (i,j)) : nothing
#     end
#
#     out
# end


text:
#!/usr/bin/env julia

using Distributed, SlurmClusterManager
addprocs(SlurmManager())

@assert nworkers() == parse(Int, ENV["SLURM_NTASKS"])

hosts = map(workers()) do id
  remotecall_fetch(() -> gethostname(), id)
end
sort!(hosts)

@assert hosts == ["c1", "c1", "c2", "c2"]


text:
using MAT
using LinearAlgebra
using SparseArrays
include("../src/Graph_Louvain.jl")
mat = matread("../data/emailEUcore.mat")
A = mat["A"]            # Largest connected component of email EU
ASim = mat["Asim"]      # Undirected version of same graph: ASim = spones(A+A')
m = sum(A)
label = vec(round.(Int64,mat["labels"]))
include("../src/Motif_Adj_Mat.jl")

Ata = ExtractTriMat(ASim)
Acc, p = largest_component(Ata)
core_inds = findall(x->x==true,p)
ASim = ASim[p,p]

## Performs weighted triangle clique expansion on an already undirected graph
H = ExtractTriangle_List(ASim)
n = size(ASim,1)
m = size(H,1)
open("../data/email.hgr","w") do f
    write(f, "$m $n\n")
    for i = 1:m
        write(f,"$(H[i,1]) $(H[i,2]) $(H[i,3])\n")
    end
end

## Do something similar for the floriday bay food web
mat = matread("../data/Florida_Bay_Dataset.mat")
A125 = mat["A125"]

H = ExtractFanList(A125)
n = size(A125,1)
m = size(H,1)
open("../data/flbay.hgr","w") do f
    write(f, "$m $n\n")
    for i = 1:m
        write(f,"$(H[i,1]) $(H[i,2]) $(H[i,3]) $(H[i,4])\n")
    end
end

max_stars_repo_name: nveldt/ParamCC

max_issues_repo_name: nveldt/ParamCC

max_forks_repo_name: nveldt/ParamCC


text:
using CSV
using Dates
using DataFrames
using DelimitedFiles
using XLSX
using YAML
using InteractiveUtils
const IU = InteractiveUtils

using SLiDE  # see src/SLiDE.jl

READ_DIR = joinpath("data", "readfiles")

files_map = XLSXInput("generate_yaml.xlsx", "map_experiment", "B1:D150", "map_experiment")
files_map = write_yaml(READ_DIR, files_map)
y_read = [read_file(files_map[ii]) for ii in 1:length(files_map)]

ind = Symbol.([m[:ind] for m in match.(r"\_(?<ind>\w*).yml", files_map)])
files_map = run_yaml(files_map)
df = [read_file(joinpath(y_read[ii]["PathOut"]...)) for ii in 1:length(y_read)];

cols = names.(df)

# for ii in 1:length(df)
#     global df[ii] = sort(edit_with(df[ii],
#         Rename.(propertynames(df[ii])[2:end], Symbol.(propertynames(df[ii])[2:end], :_, ii))))
# end

"""
Check 1: Are the two WiNDC BEA maps the same? -- YES
"""
cols_windc = intersect(names.(df[1:2])...)

df_windc = [edit_with(df[ii][:,cols_windc],
    Rename.(cols_windc[2:end], Symbol.(cols_windc[2:end], :_, ind[ii]))) for ii in 1:2]
cols_windc = sort(unique([names.(df_windc)...;]))

df_windc = outerjoin(df_windc[1], df_windc[2], on = :bea_code)
df_windc = df_windc[:,cols_windc]

df_windc_miss = df_windc[ismissing.(df_windc[:,2]), :]
dropmissing!(df_windc)

df_windc[!,:windc_det_equal] .= df_windc[:,2] .== df_windc[:,3]
df_windc[!,:windc_sum_equal] .= df_windc[:,4] .== df_windc[:,5]

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
notenames = Dict(0=>"A",1=>"Bb",2=>"B",3=>"C",4=>"Db",5=>"D",
    6=>"Eb",7=>"E",8=>"F",9=>"Gb",10=>"G",11=>"Ab")

convlist = Dict("C#"=>"Db","D#"=>"Eb","F#"=>"Gb","G#"=>"Ab","A#"=>"Bb")

"""Calculate dBSPL of pressure x"""
dBSPL(x) = 20log10(abs(x/2.0e-5))

"""
Hz to Note name
    hztonote(f,pitch=440.0)
"""
function hztonote(f,pitch=440.0)
    c = 1200log2(f/pitch)
    n = Int(mod(c÷100,12)) # note number
    cc = c%100 # cents remainder
    if cc > 50
        cc = 100-cc
        n+=1
    elseif cc <= -50
        cc = 100+cc
        n-=1
        if n < 0
            n += 12
        end
    end
    octbase = pitch*2^(-95/120) # c'-50cent
    oc =Int(floor(log2(f/octbase)))
    nm = notenames[n]
    if oc > -2
        nm = lowercase(nm)
    end
    if oc > -1
        om = "'"^(oc+1)
    elseif oc < -2
        om = ","^(-2-oc)
    else
        om = ""
    end
    nm *= om
    return nm,cc
end

function notenameconv(name)
    if name in keys(convlist)
        nm = convlist[name]
    else
        nm = name
    end
    nm
end

function notenumber(s)
    ss = notenameconv(uppercasefirst(s))
    for p in notenames
        if p.second == ss
            return p.first
        end
    end
    return -1
end

"""
Convert note and cent to Hz
"""
function notetohz(name, cent, pitch = 440.0)
    if isuppercase(name[1])
        om = replace(name,r"[^,]"=>"")
        nm = replace(name,","=>"")
        oc = -3 - length(om)
    else
        om = replace(name,r"[^']"=>"")
        nm = replace(name,"'"=>"")
        oc = length(om)-2
    end
    nn = notenumber(nm)
    if nn < 3 # lower than c
        oc+=1
    end
    f = pitch * 2^((nn*100+cent)/1200+oc)
#     @show nm,nn,oc
    return f
end


text:
#POTENTIAL ISSUE - line 22 may throw an Inexact Error for large problems

# Parameters:
#	N::Int
#	- N is the size of P
#	P::Array{Int, 1}
#	- P is a Partition
#Return Value:
#	NST::Int
#	- the number of Standard Tableau of the Young Diagram of P
function degree(N::Int, P::Array{Int, 1})
	if N > 20
		NST = 2432902008176640000.0
		for i = 1:length(P)
			for j = 1:P[i]
				NST /= hook_length(P, i, j)
			end
		end
		for n = 21:N
			NST *= n
		end
		NST = round(Int,NST)
		return NST
	else
		NST = factorial(N)
		NST /= hook_product(P)
		NST = round(Int,NST)
		return NST
	end
end

# Parameters:
#	P::Array{Int, 1}
#	- P is a Partition
# Return Values:
#	HP::Int
#	- the Hookproduct of the Young Diagram of P
function hook_product(P::Array{Int, 1})
	HP = 1
	for i = 1:length(P)
		for j = 1:P[i]
			HP *= hook_length(P, i, j)
		end
	end
	return HP
end

# Parameters:
#	P::Array{Int, 1}
#	- P is a Partition
#	i::Int
#	- the row index
#	j::Int
#	- the column index
# Return Value
#	HL::Int
#	- the HookLength of the box in the ith row and jth colomn of the Young Diagram of P
function hook_length(P::Array{Int, 1}, i::Int, j::Int)
	R = P[i] - j
	B = 0
	i += 1
	while i <= length(P) && j <= P[i]
		i += 1
		B += 1
	end
	HL = 1 + R + B
	return HL
end


text:
module LinearClassifiers

export AbstractLinearClassifier, CentroidLinearClassifier, classify

"""
Supertype of all linear classifiers.

Each linear classifier should have the fields:
- `classify::Function`
"""
abstract type AbstractLinearClassifier end

struct CentroidLinearClassifier <: AbstractLinearClassifier
    classify::Function
    # `coeffs` and `offset` are parameters for the linear decision boundary, for
    # a function of the form:
    # f(x) = (coeffs \dot x) + offset
    # where the decision boundary is defined by f(x) = 0
    coeffs::AbstractVector{<:Number}
    offset::Number
end

"""
    classify(point, model)

Classifies `point` using the `classify` function of `model`.
"""
classify(point::AbstractVector{<:Number},
         model::AbstractLinearClassifier) = model.classify(point)

end # module

max_stars_repo_name: brianxie/automaton-edification

max_issues_repo_name: brianxie/automaton-edification

max_forks_repo_name: brianxie/automaton-edification


text:
# Autogenerated wrapper script for GRASS_jll for x86_64-apple-darwin
export libgrass_arraystats, libgrass_bitmap, libgrass_btree, libgrass_btree2, libgrass_cairodriver, libgrass_calc, libgrass_ccmath, libgrass_cdhc, libgrass_cluster, libgrass_datetime, libgrass_dbmibase, libgrass_dbmiclient, libgrass_dbmidriver, libgrass_dbstubs, libgrass_dgl, libgrass_dig2, libgrass_display, libgrass_driver, libgrass_dspf, libgrass_g3d, libgrass_gis, libgrass_gmath, libgrass_gpde, libgrass_gproj, libgrass_htmldriver, libgrass_imagery, libgrass_interpdata, libgrass_interpfl, libgrass_iortho, libgrass_lidar, libgrass_linkm, libgrass_lrs, libgrass_manage, libgrass_neta, libgrass_pngdriver, libgrass_psdriver, libgrass_qtree, libgrass_raster, libgrass_rli, libgrass_rowio, libgrass_rtree, libgrass_segment, libgrass_shape, libgrass_sim, libgrass_sqlp, libgrass_stats, libgrass_symb, libgrass_temporal, libgrass_vector, libgrass_vedit

using GDAL_jll
using PROJ_jll
using Cairo_jll
using FreeType2_jll
using LAPACK_jll
using FFTW_jll
using Libiconv_jll
JLLWrappers.@generate_wrapper_header("GRASS")
JLLWrappers.@declare_library_product(libgrass_arraystats, "@rpath/libgrass_arraystats.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_bitmap, "@rpath/libgrass_bitmap.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_btree, "@rpath/libgrass_btree.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_btree2, "@rpath/libgrass_btree2.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_cairodriver, "@rpath/libgrass_cairodriver.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_calc, "@rpath/libgrass_calc.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_ccmath, "@rpath/libgrass_ccmath.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_cdhc, "@rpath/libgrass_cdhc.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_cluster, "@rpath/libgrass_cluster.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_datetime, "@rpath/libgrass_datetime.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_dbmibase, "@rpath/libgrass_dbmibase.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_dbmiclient, "@rpath/libgrass_dbmiclient.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_dbmidriver, "@rpath/libgrass_dbmidriver.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_dbstubs, "@rpath/libgrass_dbstubs.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_dgl, "@rpath/libgrass_dgl.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_dig2, "@rpath/libgrass_dig2.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_display, "@rpath/libgrass_display.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_driver, "@rpath/libgrass_driver.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_dspf, "@rpath/libgrass_dspf.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_g3d, "@rpath/libgrass_g3d.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_gis, "@rpath/libgrass_gis.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_gmath, "@rpath/libgrass_gmath.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_gpde, "@rpath/libgrass_gpde.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_gproj, "@rpath/libgrass_gproj.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_htmldriver, "@rpath/libgrass_htmldriver.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_imagery, "@rpath/libgrass_imagery.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_interpdata, "@rpath/libgrass_interpdata.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_interpfl, "@rpath/libgrass_interpfl.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_iortho, "@rpath/libgrass_iortho.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_lidar, "@rpath/libgrass_lidar.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_linkm, "@rpath/libgrass_linkm.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_lrs, "@rpath/libgrass_lrs.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_manage, "@rpath/libgrass_manage.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_neta, "@rpath/libgrass_neta.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_pngdriver, "@rpath/libgrass_pngdriver.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_psdriver, "@rpath/libgrass_psdriver.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_qtree, "@rpath/libgrass_qtree.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_raster, "@rpath/libgrass_raster.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_rli, "@rpath/libgrass_rli.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_rowio, "@rpath/libgrass_rowio.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_rtree, "@rpath/libgrass_rtree.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_segment, "@rpath/libgrass_segment.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_shape, "@rpath/libgrass_shape.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_sim, "@rpath/libgrass_sim.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_sqlp, "@rpath/libgrass_sqlp.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_stats, "@rpath/libgrass_stats.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_symb, "@rpath/libgrass_symb.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_temporal, "@rpath/libgrass_temporal.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_vector, "@rpath/libgrass_vector.7.8.dylib")
JLLWrappers.@declare_library_product(libgrass_vedit, "@rpath/libgrass_vedit.7.8.dylib")
function __init__()
    JLLWrappers.@generate_init_header(GDAL_jll, PROJ_jll, Cairo_jll, FreeType2_jll, LAPACK_jll, FFTW_jll, Libiconv_jll)
    JLLWrappers.@init_library_product(
        libgrass_arraystats,
        "lib/libgrass_arraystats.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_bitmap,
        "lib/libgrass_bitmap.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_btree,
        "lib/libgrass_btree.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_btree2,
        "lib/libgrass_btree2.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_cairodriver,
        "lib/libgrass_cairodriver.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_calc,
        "lib/libgrass_calc.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_ccmath,
        "lib/libgrass_ccmath.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_cdhc,
        "lib/libgrass_cdhc.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_cluster,
        "lib/libgrass_cluster.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_datetime,
        "lib/libgrass_datetime.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dbmibase,
        "lib/libgrass_dbmibase.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dbmiclient,
        "lib/libgrass_dbmiclient.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dbmidriver,
        "lib/libgrass_dbmidriver.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dbstubs,
        "lib/libgrass_dbstubs.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dgl,
        "lib/libgrass_dgl.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dig2,
        "lib/libgrass_dig2.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_display,
        "lib/libgrass_display.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_driver,
        "lib/libgrass_driver.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dspf,
        "lib/libgrass_dspf.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_g3d,
        "lib/libgrass_g3d.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_gis,
        "lib/libgrass_gis.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_gmath,
        "lib/libgrass_gmath.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_gpde,
        "lib/libgrass_gpde.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_gproj,
        "lib/libgrass_gproj.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_htmldriver,
        "lib/libgrass_htmldriver.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_imagery,
        "lib/libgrass_imagery.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_interpdata,
        "lib/libgrass_interpdata.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_interpfl,
        "lib/libgrass_interpfl.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_iortho,
        "lib/libgrass_iortho.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_lidar,
        "lib/libgrass_lidar.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_linkm,
        "lib/libgrass_linkm.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_lrs,
        "lib/libgrass_lrs.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_manage,
        "lib/libgrass_manage.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_neta,
        "lib/libgrass_neta.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_pngdriver,
        "lib/libgrass_pngdriver.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_psdriver,
        "lib/libgrass_psdriver.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_qtree,
        "lib/libgrass_qtree.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_raster,
        "lib/libgrass_raster.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_rli,
        "lib/libgrass_rli.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_rowio,
        "lib/libgrass_rowio.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_rtree,
        "lib/libgrass_rtree.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_segment,
        "lib/libgrass_segment.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_shape,
        "lib/libgrass_shape.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_sim,
        "lib/libgrass_sim.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_sqlp,
        "lib/libgrass_sqlp.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_stats,
        "lib/libgrass_stats.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_symb,
        "lib/libgrass_symb.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_temporal,
        "lib/libgrass_temporal.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_vector,
        "lib/libgrass_vector.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_vedit,
        "lib/libgrass_vedit.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct RecurrencePlan{T}
    A::Vector{T}
    B::Vector{T}
    C::Vector{T}
    rf₀::Vector{T}
    rf₁::Vector{T}
    rf₀inv::Vector{T}
    rf₁inv::Vector{T}
    rb₀::Vector{T}
    rb₁::Vector{T}
    rb₀inv::Vector{T}
    rb₁inv::Vector{T}
end

function RecurrencePlan(α::T,β::T,N::Int) where T
    A  = T[recA(α,β,k) for k=0:N-1]
    B  = T[recB(α,β,k) for k=0:N-1]
    C  = T[recC(α,β,k) for k=1:N  ]
    rf₀ = T[rf0(α,β,k) for k=0:N-1]
    rf₁ = T[rf1(α,β,k) for k=0:N-1]
    rf₀inv = T[inv(rf₀[k]) for k=1:N]
    rf₁inv = T[inv(rf₁[k]) for k=1:N]
    rb₀ = T[rb0(α,β,k) for k=1:N-1]
    rb₁ = T[rb1(α,β,k) for k=1:N-1]
    rb₀inv = T[inv(rb₀[k]) for k=1:N-1]
    rb₁inv = T[inv(rb₁[k]) for k=1:N-1]
    RecurrencePlan(A,B,C,rf₀,rf₁,rf₀inv,rf₁inv,rb₀,rb₁,rb₀inv,rb₁inv)
end

recA(α::T,β::T,k::Int) where T = k == 0 && ((α+β==0)||(α+β==-1)) ? (α+β)/2+1 : (2k+α+β+1)*(2k+α+β+2)/(2*(k+1)*(k+α+β+1))
recB(α::T,β::T,k::Int) where T = k==0&&((α+β==0)||(α+β==-1)) ? (α-β)/2 : (α-β)*(α+β)*(2k+α+β+1)/(2*(k+1)*(k+α+β+1)*(2k+α+β))
recC(α::T,β::T,k::Int) where T = (k+α)*(k+β)*(2k+α+β+2)/((k+1)*(k+α+β+1)*(2k+α+β))

rf0(α::T,β::T,k::Int) where T = (k+α+one(T))/(k+one(T))
rf1(α::T,β::T,k::Int) where T = -(k+β+one(T))/(k+one(T))

rb0(α::T,β::T,k::Int) where T = (k+one(T))/k*(α+β+k+1)/(k+β)*(α+β+2k)/(α+β+2k+2)
rb1(α::T,β::T,k::Int) where T = -(k+one(T))/k*(α+β+k+1)/(k+α)*(α+β+2k)/(α+β+2k+2)

# Forward polynomial recurrence. Select modified algorithms for x≈±1.

function forward_recurrence!(p::AbstractVector,N::Int,θ::Number,plan::RecurrencePlan{T}) where T
    if θ ≤ 1/4
        reinsch_f0!(p,N,sinpi(θ/2),plan)
    elseif θ ≥ 3/4
        reinsch_f1!(p,N,cospi(θ/2),plan)
    else#if 1/4 < θ < 3/4
        orthogonal_polynomial_recurrence!(p,N,cospi(θ),plan)
    end
end

forward_recurrence!(p::AbstractVector,θ::Number,plan::RecurrencePlan{T}) where T =
    forward_recurrence!(p,length(p),θ,plan)

function forward_recurrence!(p::AbstractVector,N::Int,θ::Number,cpθ::Number,spθ::Number,plan::RecurrencePlan{T}) where T
    if θ ≤ 1/4
        reinsch_f0!(p,N,spθ,plan)
    elseif θ ≥ 3/4
        reinsch_f1!(p,N,cpθ,plan)
    else#if 1/4 < θ < 3/4
        orthogonal_polynomial_recurrence!(p,N,cospi(θ),plan)
    end
end

# Forward recurrence

function orthogonal_polynomial_recurrence!(p::AbstractVector,N::Int,x::Number,plan::RecurrencePlan{V}) where V
    A=plan.A
    B=plan.B
    C=plan.C

    T = promote_type(eltype(p),typeof(x),V)
    pk1, pk = one(T), muladd(A[1],x,B[1])
    p[1],p[2] = pk1, pk

    @inbounds for k = 2:N-1
        Ak, Bk, Ck = A[k], B[k], C[k-1]
        pk, pk1 = muladd(muladd(Ak,x,Bk),pk,-Ck*pk1), pk
        p[k+1] = pk
    end
end

# Modified forward recurrence algorithm for θ near 0
# spθ = sinpi(θ/2)

function reinsch_f0!(p::AbstractVector,N::Int,spθ::Number,plan::RecurrencePlan{V}) where V
    A=plan.A
    C=plan.C
    r=plan.rf₀
    rinv=plan.rf₀inv

    T = promote_type(eltype(p),typeof(spθ),V)
    xm1 = -2spθ^2
    p[1] = one(T)
    dk = A[1]*xm1*rinv[1]
    p[2] = (p[1]+dk)*r[1]

    @inbounds for k = 2:N-1
        Ak, Ck, rk, rkinv = A[k], C[k-1], r[k], rinv[k]
        temp = muladd(Ak*xm1,p[k],Ck*dk)
        dk, p[k+1] = temp*rkinv, muladd(rk,p[k],temp)
    end
end

# Modified forward recurrence for θ near 1
# cpθ = cospi(θ/2)

function reinsch_f1!(p::AbstractVector,N::Int,cpθ::Number,plan::RecurrencePlan{V}) where V
    A=plan.A
    C=plan.C
    r=plan.rf₁
    rinv=plan.rf₁inv

    T = promote_type(eltype(p),typeof(cpθ),V)
    xp1 = 2cpθ^2
    p[1] = one(T)
    dk = A[1]*xp1*rinv[1]
    p[2] = (p[1]+dk)*r[1]

    @inbounds for k = 2:N-1
        Ak, Ck, rk, rkinv = A[k], C[k-1], r[k], rinv[k]
        temp = muladd(Ak*xp1,p[k],Ck*dk)
        dk, p[k+1] = temp*rkinv, muladd(rk,p[k],temp)
    end
end

# Backward polynomial recurrence. Select modified algorithms for x≈±1.

function backward_recurrence(c::AbstractVector,N::Int,θ::Number,plan::RecurrencePlan{T}) where T
    if θ ≤ 1/4
        reinsch_b0(c,N,sinpi(θ/2),plan)
    elseif θ ≥ 3/4
        reinsch_b1(c,N,cospi(θ/2),plan)
    else#if 1/4 < θ < 3/4
        clenshaw(c,N,cospi(θ),plan)
    end
end

backward_recurrence(c::AbstractVector,θ::Number,plan::RecurrencePlan{T}) where T =
    backward_recurrence(c,length(c),θ,plan)

function backward_recurrence(c::AbstractVector,N::Int,θ::Number,cpθ::Number,spθ::Number,plan::RecurrencePlan{T}) where T
    if θ ≤ 1/4
        reinsch_b0(c,N,spθ,plan)
    elseif θ ≥ 3/4
        reinsch_b1(c,N,cpθ,plan)
    else#if 1/4 < θ < 3/4
        clenshaw(c,N,cospi(θ),plan)
    end
end

# Clenshaw-Smith algorithm

function clenshaw(c::AbstractVector,N::Int,x::Number,plan::RecurrencePlan{V}) where V
    A=plan.A
    B=plan.B
    C=plan.C

    T = promote_type(eltype(c),typeof(x),V)
    bk, bk1, bk2 = zero(T), zero(T), zero(T)

    @inbounds for k = N:-1:2
        ck, Ak, Bk, Ck = c[k], A[k], B[k], C[k]
        bk = muladd(muladd(Ak,x,Bk),bk1,muladd(-Ck,bk2,ck))
        bk2, bk1, bk = bk1, bk, bk2
    end

    @inbounds ck, Ak, Bk, Ck = c[1], A[1], B[1], C[1]
    muladd(muladd(Ak,x,Bk),bk1,muladd(-Ck,bk2,ck))
end

# Modified Clenshaw-Smith algorithm for θ near 0
# spθ = sinpi(θ/2)

function reinsch_b0(c::AbstractVector,N::Int,spθ::Number,plan::RecurrencePlan{V}) where V
    A=plan.A
    C=plan.C
    r=plan.rb₀
    rinv=plan.rb₀inv

    T = promote_type(eltype(c),typeof(spθ),V)
    xm1 = -2spθ^2
    uk, dk = zero(T), zero(T)

    @inbounds for k=N:-1:2
        ck, Ak, Ck, rk, rkinv = c[k], A[k], C[k], r[k-1], rinv[k-1]
        temp = muladd(Ak*xm1,uk,muladd(Ck,dk,ck))
        dk, uk = temp*rk, muladd(rkinv,uk,temp)
    end

    @inbounds ck, Ak, Ck = c[1], A[1], C[1]
    muladd(Ak*xm1,uk,muladd(Ck,dk,ck))
end

# Modified Clenshaw-Smith algorithm for θ near 1
# cpθ = cospi(θ/2)

function reinsch_b1(c::AbstractVector,N::Int,cpθ::Number,plan::RecurrencePlan{V}) where V
    A=plan.A
    C=plan.C
    r=plan.rb₁
    rinv=plan.rb₁inv

    T = promote_type(eltype(c),typeof(cpθ),V)
    xp1 = 2cpθ^2
    uk, dk = zero(T), zero(T)

    @inbounds for k=N:-1:2
        ck, Ak, Ck, rk, rkinv = c[k], A[k], C[k], r[k-1], rinv[k-1]
        temp = muladd(Ak*xp1,uk,muladd(Ck,dk,ck))
        dk, uk = temp*rk, muladd(rkinv,uk,temp)
    end

    @inbounds ck, Ak, Ck = c[1], A[1], C[1]
    muladd(Ak*xp1,uk,muladd(Ck,dk,ck))
end


text:
import Nemo.sub!, Base.gcd
export induce_rational_reconstruction, induce_crt, root, roots,
       is_monic, radical_extension, is_radical_extension,
       is_kummer_extension, cyclotomic_field, wildanger_field,
       compositum

add_verbose_scope(:PolyFactor)
add_assert_scope(:PolyFactor)
add_verbose_scope(:CompactPresentation)
add_assert_scope(:CompactPresentation)

if Int==Int32
  global const p_start = next_prime(2^30)
else
  global const p_start = next_prime(2^60)
end

include("NfAbs/NfAbs.jl")
include("NfAbs/Conjugates.jl")
include("NfAbs/ComplexEmbeddings.jl")
include("NfAbs/CompactRepresentation.jl")
include("NfAbs/Elem.jl")
include("NfAbs/NonSimple.jl")
include("NfAbs/Poly.jl")
include("NfAbs/Simplify.jl")
include("NfAbs/DiscriminantBounds.jl")
include("NfAbs/NormRelation.jl")
include("NfAbs/PolyFact.jl")
include("NfAbs/MPolyGcd.jl")
include("NfAbs/MPolyFactor.jl")
include("NfAbs/MPolyAbsFact.jl")
include("NfAbs/ConjugatesNS.jl")
include("NfAbs/ComplexEmbeddingsNS.jl")


text:
using IVPTestSuite  # imports ODE, DASSL and Sundials
const S = Solvers # for convenience

# choose test-case (all of them are in IVPTestSuite.tc_all)
tc_name = [:vdpol
           :rober
           :threebody
           :bruss1d
           :chemakzo
           :hires][3]
tc = IVPTestSuite.tc_all[tc_name]

# Pick a solver.  S.allsolvers is the list of all solvers.  Generally
# each family of solvers also has a list.
solver = S.allsolvers[ODE.ode45_dp]

###
# make a TestRun which combines a TestCase with a Solver + some extras:
abstols = 10.0.^(-5:-1:-10)
reltols = abstols
dt0 = NaN # size of first step
suite = TestSuite(tc, solver, abstols, reltols, [NaN])

# running it will display some results:
res = run_ode_testsuite(suite)


text:
function my_func1(x::Int64, y::Int64)::Int64
end

bit::BitVec{1}
myu32::BitVec{32}
myu64::BitVec{64}


text:
# Autogenerated wrapper script for SLATEC_jll for x86_64-w64-mingw32-libgfortran3
export libslatec

using LAPACK_jll
using OpenBLAS_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("SLATEC")
JLLWrappers.@declare_library_product(libslatec, "libslatec.dll")
function __init__()
    JLLWrappers.@generate_init_header(LAPACK_jll, OpenBLAS_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libslatec,
        "bin\\libslatec.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using LazyArrays, ArrayLayouts, LinearAlgebra, FillArrays
import LazyArrays: materialize!, MemoryLayout, triangulardata, LazyLayout, UnknownLayout, LazyMatrix

# used to test general matrix backends
struct MyMatrix{T} <: LazyMatrix{T}
    A::Matrix{T}
end

MyMatrix{T}(::UndefInitializer, n::Int, m::Int) where T = MyMatrix{T}(Array{T}(undef, n, m))
MyMatrix(A::AbstractMatrix{T}) where T = MyMatrix{T}(Matrix{T}(A))
Base.convert(::Type{MyMatrix{T}}, A::MyMatrix{T}) where T = A
Base.convert(::Type{MyMatrix{T}}, A::MyMatrix) where T = MyMatrix(convert(AbstractArray{T}, A.A))
Base.convert(::Type{MyMatrix}, A::MyMatrix)= A
Base.convert(::Type{AbstractArray{T}}, A::MyMatrix) where T = MyMatrix(convert(AbstractArray{T}, A.A))
Base.convert(::Type{AbstractMatrix{T}}, A::MyMatrix) where T = MyMatrix(convert(AbstractArray{T}, A.A))
Base.convert(::Type{MyMatrix{T}}, A::AbstractArray{T}) where T = MyMatrix{T}(A)
Base.convert(::Type{MyMatrix{T}}, A::AbstractArray) where T = MyMatrix{T}(convert(AbstractArray{T}, A))
Base.convert(::Type{MyMatrix}, A::AbstractArray{T}) where T = MyMatrix{T}(A)
Base.getindex(A::MyMatrix, kj...) = A.A[kj...]
Base.getindex(A::MyMatrix, ::Colon, j::Integer) = A.A[:,j]
Base.getindex(A::MyMatrix, ::Colon, j::AbstractVector) = MyMatrix(A.A[:,j])
Base.setindex!(A::MyMatrix, v, kj...) = setindex!(A.A, v, kj...)
Base.size(A::MyMatrix) = size(A.A)
Base.similar(::Type{MyMatrix{T}}, m::Int, n::Int) where T = MyMatrix{T}(undef, m, n)
Base.similar(::MyMatrix{T}, m::Int, n::Int) where T = MyMatrix{T}(undef, m, n)
Base.similar(::MyMatrix, ::Type{T}, m::Int, n::Int) where T = MyMatrix{T}(undef, m, n)
LinearAlgebra.factorize(A::MyMatrix) = factorize(A.A)

struct MyLazyArray{T,N} <: LazyArray{T,N}
    data::Array{T,N}
end

Base.size(A::MyLazyArray) = size(A.data)
Base.getindex(A::MyLazyArray, j::Int...) = A.data[j...]
LazyArrays.MemoryLayout(::Type{<:MyLazyArray}) = LazyLayout()
LinearAlgebra.factorize(A::MyLazyArray) = factorize(A.data)

@testset "lazymul/ldiv tests" begin
    @testset "*" begin
        A = randn(5,5)
        B = randn(5,5)
        x = randn(5)
        @test MyMatrix(A)*x ≈ apply(*,MyMatrix(A),x) ≈ A*x
        @test MemoryLayout(MyMatrix(A)) isa LazyLayout
        @test all(MyMatrix(A)*MyMatrix(A) .=== apply(*,MyMatrix(A),MyMatrix(A)))
        @test all(MyMatrix(A)*A .=== apply(*,MyMatrix(A),A))
        @test all(A*MyMatrix(A) .=== apply(*,A,MyMatrix(A)))
        @test MyMatrix(A)*MyMatrix(A) ≈ MyMatrix(A)*A ≈ A*MyMatrix(A) ≈ A^2

        @test MyMatrix(A)*MyMatrix(A)*MyMatrix(A) ≈ apply(*,MyMatrix(A),MyMatrix(A),MyMatrix(A)) ≈ A^3

        @test all(UpperTriangular(A) * MyMatrix(A) .=== apply(*,UpperTriangular(A), MyMatrix(A)))
        @test all(MyMatrix(A) * UpperTriangular(A) .=== apply(*, MyMatrix(A),UpperTriangular(A)))
        @test all(Diagonal(A) * MyMatrix(A) .=== apply(*,Diagonal(A), MyMatrix(A)))
        @test all(MyMatrix(A) * Diagonal(A) .=== apply(*, MyMatrix(A),Diagonal(A)))

        @test all(MyMatrix(A)' * x .=== apply(*,MyMatrix(A)',x))

        @test all(MyMatrix(A)' * MyMatrix(A)' .=== apply(*,MyMatrix(A)', MyMatrix(A)'))
        @test all(MyMatrix(A)' * A' .=== apply(*,MyMatrix(A)', A'))
        @test all(A' * MyMatrix(A)' .=== apply(*,MyMatrix(A)', MyMatrix(A)'))
        @test all(MyMatrix(A)' * MyMatrix(A) .=== apply(*,MyMatrix(A)', MyMatrix(A)))
        @test all(MyMatrix(A)' * A .=== apply(*,MyMatrix(A)', A))
        @test all(MyMatrix(A) * MyMatrix(A)' .=== apply(*,MyMatrix(A), MyMatrix(A)'))
        @test all(A * MyMatrix(A)' .=== apply(*,A, MyMatrix(A)'))

        @test all(UpperTriangular(A) * MyMatrix(A) .=== apply(*,UpperTriangular(A), MyMatrix(A)))
        @test all(MyMatrix(A) * UpperTriangular(A) .=== apply(*, MyMatrix(A),UpperTriangular(A)))

        @test all(Diagonal(A) * MyMatrix(A)' .=== apply(*,Diagonal(A), MyMatrix(A)'))
        @test all(MyMatrix(A)' * Diagonal(A) .=== apply(*,MyMatrix(A)',Diagonal(A)))
        @test all(UpperTriangular(A) * MyMatrix(A)' .=== apply(*,UpperTriangular(A), MyMatrix(A)'))
        @test all(MyMatrix(A)' * UpperTriangular(A) .=== apply(*,MyMatrix(A)',UpperTriangular(A)))

        @test ApplyArray(\, MyMatrix(A), x)[1,1] ≈ (A\x)[1]
        @test MyMatrix(A)\x ≈ apply(\,MyMatrix(A),x) ≈ copyto!(similar(x),Ldiv(A,copy(x))) ≈ A\x
        @test eltype(applied(\,MyMatrix(A),x)) == eltype(apply(\,MyMatrix(A),x)) == eltype(MyMatrix(A)\x) == Float64

        @test MyMatrix(A)\MyMatrix(B) ≈ MyMatrix(A)\B ≈ apply(\,MyMatrix(A),B) ≈ copyto!(similar(B),Ldiv(A,copy(B))) ≈ A\B
        @test eltype(applied(\,MyMatrix(A),B)) == eltype(apply(\,MyMatrix(A),B)) == eltype(MyMatrix(A)\B) == Float64

        @test MyMatrix(A) * ApplyArray(exp,B) ≈ apply(*, MyMatrix(A),ApplyArray(exp,B)) ≈ A*exp(B)
        @test ApplyArray(exp,A) * MyMatrix(B)  ≈ apply(*, ApplyArray(exp,A), MyMatrix(B)) ≈ exp(A)*B
        @test ApplyArray(exp,A) * ApplyArray(exp,B) ≈ apply(*, ApplyArray(exp,A),ApplyArray(exp,B)) ≈ exp(A)*exp(B)

        @test MyMatrix(A) * BroadcastArray(exp,B) ≈ apply(*, MyMatrix(A),BroadcastArray(exp,B)) ≈ A*exp.(B)
        @test BroadcastArray(exp,A) * MyMatrix(B)  ≈ apply(*, BroadcastArray(exp,A), MyMatrix(B)) ≈ exp.(A)*B
        @test BroadcastArray(exp,A) * BroadcastArray(exp,B) ≈ apply(*, BroadcastArray(exp,A),BroadcastArray(exp,B)) ≈ exp.(A)*exp.(B)
    end

    @testset "\\" begin
        A = randn(5,5)
        B = randn(5,5)
        x = randn(5)
        @test MyMatrix(A) \ x == apply(\, MyMatrix(A), x)
        @test ldiv!(MyMatrix(A), copy(x)) == materialize!(Ldiv(MyMatrix(A), copy(x)))
        @test MyMatrix(A) \ x ≈ ldiv!(MyMatrix(A), copy(x)) ≈ A\x
        @test MyMatrix(A) \ B == apply(\, MyMatrix(A), B)
        @test ldiv!(MyMatrix(A), copy(B)) == materialize!(Ldiv(MyMatrix(A), copy(B)))
        @test MyMatrix(A) \ B ≈ MyMatrix(A) \ MyMatrix(B) ≈ ldiv!(MyMatrix(A), copy(B)) ≈  A\B
        @test_broken ldiv!(MyMatrix(A), MyMatrix(copy(B))) ≈ A\B

        C = randn(5,3)
        @test all(MyMatrix(C)\x .=== apply(\,MyMatrix(C),x))
        @test MyMatrix(C)\x ≈ C\x
        @test all(MyMatrix(C)\B .=== apply(\,MyMatrix(C),B))
        @test MyMatrix(C)\B ≈ C\B

        @test_throws DimensionMismatch apply(\,MyMatrix(C),randn(4))
        @test_throws DimensionMismatch apply(\,MyMatrix(C),randn(4,3))
    end

    @testset "Lazy" begin
        A = MyLazyArray(randn(2,2))
        B = MyLazyArray(randn(2,2))
        x = MyLazyArray(randn(2))

        @test apply(*,A,x) isa ApplyVector
        @test apply(*,A,Array(x)) isa ApplyVector
        @test apply(*,Array(A),x) isa ApplyVector
        @test apply(*,A,x) ≈ apply(*,Array(A),x) ≈ apply(*,A,Array(x))  ≈ Array(A)*Array(x)

        @test apply(*,A,B) isa ApplyMatrix
        @test apply(*,A,Array(B)) isa ApplyMatrix
        @test apply(*,Array(A),B) isa ApplyMatrix
        @test apply(*,A,B) ≈ apply(*,Array(A),B) ≈ apply(*,A,Array(B))  ≈ Array(A)*Array(B)

        @test apply(\,A,x) isa ApplyVector
        @test apply(\,A,Array(x)) isa ApplyVector
        @test apply(\,Array(A),x) isa ApplyVector
        @test apply(\,A,x) ≈ apply(\,Array(A),x) ≈ apply(\,A,Array(x)) ≈ Array(A)\Array(x)

        @test apply(\,A,B) isa ApplyMatrix
        @test apply(\,A,Array(B)) isa ApplyMatrix
        @test apply(\,Array(A),B) isa ApplyMatrix
        @test apply(\,A,B) ≈ apply(\,Array(A),B) ≈ apply(\,A,Array(B)) ≈ Array(A)\Array(B)

        Ap = applied(*,A,x)
        @test copyto!(similar(Ap), Ap) == A*x
        @test copyto!(similar(Ap,BigFloat), Ap) ≈ A*x

        @test MemoryLayout(typeof(Diagonal(x))) isa DiagonalLayout{LazyLayout}
        @test MemoryLayout(typeof(Diagonal(ApplyArray(+,x,x)))) isa DiagonalLayout{LazyLayout}
        @test MemoryLayout(typeof(Diagonal(1:6))) isa DiagonalLayout{UnknownLayout}

        @test MemoryLayout(typeof(A')) isa LazyLayout
        @test MemoryLayout(typeof(transpose(A))) isa LazyLayout
        @test MemoryLayout(typeof(view(A,1:2,1:2))) isa LazyLayout
        @test MemoryLayout(typeof(reshape(A,4))) isa LazyLayout
    end

    @testset "QR" begin
        B = MyMatrix(randn(3,3))
        Q = qr(randn(3,3)).Q
        @test Q * B ≈ Q*B.A
    end

    @testset "ambiguities" begin
        A = randn(5,5)
        b = MyLazyArray(randn(5))
        c = randn(5)
        c̃ = complex.(c)
        @test A*b isa ApplyVector{Float64,typeof(*)}
        @test UpperTriangular(A)*b isa ApplyVector{Float64,typeof(*)}
        @test A*b ≈ A*Vector(b)
        @test UpperTriangular(A)*b ≈ UpperTriangular(A)*Vector(b)
        @test c'b ≈ c̃'b ≈ c'Vector(b)
        @test transpose(c)b ≈ transpose(c̃)b ≈ transpose(c)Vector(b)
    end

    @testset "InvMatrix" begin
        A = randn(5,5)
        B = randn(5,5)
        b = MyLazyArray(randn(5))
        M = ApplyArray(*, B, b)
        @test InvMatrix(A) * b ≈ A \ b
        @test InvMatrix(A) * M ≈ A \ B * b

        @test ArrayLayouts.ldiv(MyMatrix(A), M) ≈ A\ B * b
    end

    @testset "Tri/Diagonal" begin
        b = MyLazyArray(randn(5))
        c = MyLazyArray(randn(4))
        d = MyLazyArray(randn(3))
        @test copy(Diagonal(b)) == Diagonal(copy(b))
        @test map(copy, Diagonal(b))  == Diagonal(copy(b))
        @test inv(Diagonal(b)) == inv(Diagonal(b.data))
        @test inv(Diagonal(b)) isa Diagonal{Float64,<:BroadcastVector}
        @test copy(Tridiagonal(c, b, c)) == Tridiagonal(copy(c), copy(b), copy(c))
        @test copy(Tridiagonal(c, b, c, d)) == Tridiagonal(copy(c), copy(b), copy(c), copy(d))
        @test copy(Tridiagonal(c, b, c, d)).du2 == d
        @test map(copy, Tridiagonal(c, b, c)) == Tridiagonal(copy(c), copy(b), copy(c))
        @test map(copy, Tridiagonal(c, b, c, d)) == Tridiagonal(copy(c), copy(b), copy(c), copy(d))
        @test map(copy, Tridiagonal(c, b, c, d)).du2 == d

        @test MemoryLayout(Tridiagonal(c, b, c)) isa TridiagonalLayout{LazyLayout,LazyLayout,LazyLayout}
        @test MemoryLayout(SymTridiagonal(b, c)) isa SymTridiagonalLayout{LazyLayout,LazyLayout}
        @test MemoryLayout(Bidiagonal(b, c, :U)) isa BidiagonalLayout{LazyLayout,LazyLayout}

        @test LazyArrays.tridiagonallayout(UnknownLayout(), UnknownLayout(), LazyLayout()) isa TridiagonalLayout{LazyLayout,LazyLayout,LazyLayout}
        @test LazyArrays.tridiagonallayout(UnknownLayout(), LazyLayout(), UnknownLayout()) isa TridiagonalLayout{LazyLayout,LazyLayout,LazyLayout}
        @test LazyArrays.tridiagonallayout(LazyLayout(), UnknownLayout(), UnknownLayout()) isa TridiagonalLayout{LazyLayout,LazyLayout,LazyLayout}
        @test LazyArrays.tridiagonallayout(UnknownLayout(), LazyLayout(), LazyLayout()) isa TridiagonalLayout{LazyLayout,LazyLayout,LazyLayout}
        @test LazyArrays.tridiagonallayout(LazyLayout(), UnknownLayout(), LazyLayout()) isa TridiagonalLayout{LazyLayout,LazyLayout,LazyLayout}
        @test LazyArrays.tridiagonallayout(LazyLayout(), LazyLayout(), UnknownLayout()) isa TridiagonalLayout{LazyLayout,LazyLayout,LazyLayout}

        @test LazyArrays.symtridiagonallayout(UnknownLayout(), LazyLayout()) isa SymTridiagonalLayout{LazyLayout,LazyLayout}
        @test LazyArrays.symtridiagonallayout(LazyLayout(), UnknownLayout()) isa SymTridiagonalLayout{LazyLayout,LazyLayout}

        @test LazyArrays.bidiagonallayout(UnknownLayout(), LazyLayout()) isa BidiagonalLayout{LazyLayout,LazyLayout}
        @test LazyArrays.bidiagonallayout(LazyLayout(), UnknownLayout()) isa BidiagonalLayout{LazyLayout,LazyLayout}
    end

    @testset "Nested" begin
        a = MyLazyArray(randn(5))
        @test a .\ rand(5) .* Zeros(5) ≡ Zeros(5)
        @test broadcast(*, Zeros(5), Base.broadcasted(\, a, rand(5))) ≡ Zeros(5)
    end

    @testset "inv" begin
        A = randn(5,5)
        B = randn(5,5)
        M = ApplyArray(*, A, B)
        b = randn(5)
        @test M \ MyLazyArray(b) ≈ M \ b
    end

    @testset "Diagonal Fill" begin
        b = randn(5)
        B = randn(5,5)
        @test Eye(5) * MyLazyArray(b) == b
        @test MyLazyArray(B) * Eye(5) == B
    end

    @testset "LazyBroadcast" begin
        a = MyLazyArray(randn(5))
        b = a .^ 2
        @test BroadcastArray(view(b,1:3)) == Vector(a)[1:3] .^2
    end

    @testset "Apply*Broadcast" begin
        A = randn(5,5)
        B = randn(5,5)

        @test ApplyArray(*, A, B) * BroadcastArray(*, A, B) ≈ (A*B) * (A .* B)
        @test BroadcastArray(*, A, B) * ApplyArray(*, A, B) ≈ (A .* B) * (A*B)

        @test ApplyArray(*, A, B) \ BroadcastArray(*, A, B) ≈ (A*B) \ (A .* B)
        @test BroadcastArray(*, A, B) \ ApplyArray(*, A, B) ≈ (A .* B) \ (A * B)
        @test BroadcastArray(*, A, B) \ BroadcastArray(*, A, B) ≈ (A .* B) \ (A .* B)
    end

    @testset "CartesianIndex view" begin
        A = randn(5,5)
        B = randn(5,5)
        M = ApplyArray(*, A, B)
        @test layout_getindex(M,[CartesianIndex(1,2),CartesianIndex(3,3)]) == M[[CartesianIndex(1,2),CartesianIndex(3,3)]] == [M[1,2], M[3,3]]
    end
end