text:
module Deprecations
    using CSTParser
    using CSTParser: EXPR, MacroCall, IDENTIFIER, LITERAL
    using AbstractTrees
    using Tokenize: Tokens

    export edit_text, edit_file, edit_markdown

    struct Deprecation
        description::String
        package::String
        minsupported::Base.VersionNumber
        minrequired::Base.VersionNumber
        # Don't apply this deprecation if all versions are above
        # this version - generally because this syntax is now available
        # for re-use, and we want to avoid false positives for such syntax
        maxver::Base.VersionNumber
    end

    include("CSTUtils/CSTUtils.jl")
    using .CSTUtils
    import .CSTUtils: span_text, fullspan_text, leading_trivia, trailing_trivia
    include("CSTAnalyzer/CSTAnalyzer.jl")
    using .CSTAnalyzer
    using .CSTAnalyzer: State, FileSystem, Scope, Location, File
    include("formatting.jl")
    include("astmatching.jl")
    include("resolutions.jl")

    all_templates = Any[]
    all_replacements = Any[]
    deprecation = Any[]

    function applicable_dep(dep, vers)
        i = info(dep)
        haskey(vers, i.package) || return false
        pvers = vers[i.package]
        all(interval->i.minsupported <= interval.lower, pvers.intervals) || return false
        all(interval->i.maxver <= interval.lower, pvers.intervals) && return false
        true
    end

    all_deprecations = Dict()
    templates = Dict{Any,Vector{Tuple{Any, Any, Any, Any}}}()
    custom_resolutions = Dict{Any, Any}()
    register(dep, info) = all_deprecations[dep] = info
    info(dep) = all_deprecations[dep]

    dep_for_vers(dep, vers) = dep()

    function applicable_deprecations(vers)
        deps = Any[]
        for dep in keys(all_deprecations)
            applicable_dep(dep, vers) || continue
            push!(deps, dep_for_vers(dep, vers))
        end
        deps
    end

    identity_formatter(orig_tree, tree, matches) = tree

    function match(dep, template::String, replacement::String, formatter = identity_formatter; filter = (args...)->return true)
        haskey(templates, dep) || (templates[dep] = Vector{Tuple{Any, Any, Any}}())
        push!(templates[dep], (template, replacement, formatter, filter))
    end

    function match_macrocall(mod, name, dep, template::String, replacement::String)
        # TODO: Make sure the macro in the current module is actually the one from `mod`
        haskey(templates, dep) || (templates[dep] = Vector{Tuple{Any, Any, Any}}())
        push!(templates[dep], (template, replacement, identity_formatter))
    end

    function match(f::Function, dep, expr_kind)
        haskey(custom_resolutions, dep) || (custom_resolutions[dep] = Vector{Any}())
        push!(custom_resolutions[dep], (expr_kind, f))
    end

    include("database.jl")

    struct Context
        in_macrocall::Bool
        top_macrocall
    end

    # By default, don't apply rewrites in macro context
    applies_in_macrocall(dep, context) = false

    struct TextReplacement
        dep::Union{Deprecation, Nothing}
        range::UnitRange{Int}
        text::String
    end
    function TextReplacement(dep, range::UnitRange{Int}, text::String)
        dep = all_deprecations[typeof(dep)]
        TextReplacement(dep, range, text)
    end

    function overlay_parse(text, cont = true)
        p = CSTParser.parse(text, cont)
        OverlayNode(p, text)
    end

    function changed_text(text, resolutions)
        sort!(resolutions, by=x->first(x.range))
        buf = IOBuffer()
        lastoffset = 1
        for r in resolutions
            if lastoffset > first(r.range)
                # Overlapping replacements. Only apply the first one for now.
                # We'll re-run this to convergence
                continue
            end
            write(buf, text[(lastoffset:first(r.range)-1)])
            write(buf, r.text)
            lastoffset = last(r.range)+1
        end
        write(buf, text[lastoffset:end])
        (length(resolutions) != 0, String(take!(buf)))
    end

    struct ParsedReplacement
        dep
        t::CSTUtils.OverlayNode
        r::CSTUtils.OverlayNode
        formatter
        filter
    end

    struct CustomReplacement
        dep
        matchT::Type
        f
    end

    function text_replacements(text, deps; analysis = nothing)
        match = overlay_parse(text)
        # Re-construct analysis by assuming top-level
        if analysis === nothing
            S = State{FileSystem}(Scope(), Location("top", 0), "", [], [], 0:0, false, Dict(), FileSystem())
            Scop = S.current_scope
            CSTAnalyzer.trav(CSTAnalyzer.DefaultWalker(), match, Scop, S)
        else
            (S, Scop) = analysis
        end
        # CST Matchers
        replacements = Any[]
        customs = Any[]
        for dep in deps
            haskey(templates, typeof(dep)) && append!(replacements, collect(Iterators.product((dep,), templates[typeof(dep)])))
            haskey(custom_resolutions, typeof(dep)) && append!(customs,
                map(x->CustomReplacement(x[1], x[2]...),
                    collect(Iterators.product((dep,), custom_resolutions[typeof(dep)])))
            )
        end
        parsed_replacementes::Vector{ParsedReplacement} = map(x->ParsedReplacement(x[1], overlay_parse(x[2][1],false),overlay_parse(x[2][2],false),x[2][3:end]...), replacements)
        function find_replacements(x, results, context=Context(false, nothing))
            for (i, pr) in enumerate(parsed_replacementes)
                if matches_template2(x, pr.t)
                    (!context.in_macrocall || applies_in_macrocall(pr.dep, context)) || continue
                    result = Dict{Any,Any}()
                    match_parameters(pr.t, x, result)[1] || continue
                    pr.filter((S, Scop), pr.dep, x, result) || continue
                    rtree = reassemble_tree(pr.r, result)
                    buf = IOBuffer()
                    print_replacement(buf, pr.formatter(x, rtree, result))
                    push!(results, TextReplacement(pr.dep, x.span, String(take!(buf))))
                end
            end
            for (i, custom) in enumerate(customs)
                if isexpr(x, custom.matchT)
                    custom.f((custom.dep, x, results, context, (S, Scop)))
                end
            end
            if isexpr(x, MacroCall) && !context.in_macrocall
                # Don't consider doc macros here. We know they don't
                # affect the meaning of the parsed code.
                if !isexpr(children(x)[1], CSTParser.GlobalRefDoc)
                    macroname = children(children(x)[1])[2]
                    if !(isexpr(macroname, IDENTIFIER) && Expr(macroname) in (
                           Symbol("eval"), Symbol("inline"), Symbol("views"), Symbol("test"),
                           Symbol("testset"), Symbol("pure"), Symbol("test_throws"),
                           Symbol("inferred")
                            ))
                        context = Context(true, x)
                    end
                end
                # Recurse into code examples in documentation
                if isexpr(children(x)[1], CSTParser.GlobalRefDoc)
                    docnode = children(x)[2]
                    doc = span_text(docnode)
                    nquotes = doc[1:3] == "\"\"\"" ? 3 : 1
                    doc = doc[(1+nquotes):(end-nquotes)]
                    if !isa(doc, Expr)
                        changed, new_doc = edit_markdown(doc, deps)
                        if changed
                            # TODO: More fine grained deprecation messages inside
                            # doc strings
                            push!(results, TextReplacement(nothing, (first(docnode.span)+nquotes):(last(docnode.span)-nquotes), new_doc))
                        end
                    end
                end
            end
            for arg in children(x)
                find_replacements(arg, results, context)
            end
        end
        results = TextReplacement[]
        find_replacements(match, results)
        results
    end

    function edit_text(text, deps = map(x->x(), keys(all_deprecations)); kwargs...)
        changed_text(text, text_replacements(text, deps; kwargs...))
    end

    function edit_file(fname, deps = map(x->x(), keys(all_deprecations)), edit=edit_text; kwargs...)
        text = readstring(fname)
        any_changed, new_text = edit(text, deps; kwargs...)
        if any_changed
            open(fname, "w") do io
                write(io, new_text)
            end
        end
        any_changed
    end

    function edit_markdown(text, deps = map(x->x(), keys(all_deprecations)); kwargs...)
        content = Base.Markdown.parse(text)
        changed_any = false
        new_text = text
        map(content.content) do x
            isa(x, Base.Markdown.Code) || return
            if !(x.language =="julia" || ismatch(r"jldoctest[ ]?(.*)$", x.language))
                return
            end
            any_matches = false
            new_code = x.code
            for m in eachmatch(r"^julia> "m, x.code)
                any_matches = true
                curidx = startidx = m.offset
                while (curidx = search(x.code, '\n', curidx)) != 0
                    # In order to be a line continuation, needs to be at least
                    # 7 spaces after a newline.
                    is_continuation = true
                    idx = curidx
                    i = 1
                    while i <= 7
                        idx = nextind(x.code, idx)
                        idx > endof(x.code) && break
                        if !isspace(x.code[idx])
                            is_continuation = false
                            break
                        end
                        i += 1
                    end
                    is_continuation || break
                    curidx = nextind(x.code, curidx)
                end
                curidx == 0 && (curidx = endof(x.code))
                text = x.code[(startidx+7):curidx]
                try
                    new_code = replace(new_code, text, edit_text(text, deps)[2])
                end
            end
            if !any_matches
                # Parse the whole thing
                try
                    new_code = edit_text(new_code, deps)[2]
                end
            end
            if x.code != new_code
                changed_any = true
                new_text = replace(new_text, x.code, new_code)
            end
        end
        return (changed_any, new_text)
    end

    struct IncludeWalker <: CSTAnalyzer.Walker
        include_map::Vector{Pair{Scope, String}}
    end
    IncludeWalker() = IncludeWalker(Vector{Pair{Scope, String}}())

    function CSTAnalyzer.lint_call(w::IncludeWalker, x::CSTParser.EXPR{CSTParser.Call}, s, S)
        fname = CSTParser.get_name(x)
        if CSTParser.str_value(fname) == "include"
            path = CSTAnalyzer.get_path(x)
            isempty(path) && return
            path = isabspath(path) ? path : joinpath(dirname(s.loc.path), path)
            push!(w.include_map, s=>path)
        end
    end

    function scope_include!(parent, scope)
        push!(parent.children, scope)
        scope.parent = parent
        merge!(parent.names, scope.names)
        empty!(scope.names)
    end

    # Given a set of files in a package, process all of them
    # using CSTAnalyzer
    function process_all(files)
        parsed = [ file=>overlay_parse(readstring(file)) for file in files ]
        S = State{FileSystem}(Scope(), Location("top", 0), "", [], [], 0:0, false, Dict(), FileSystem())
        walker = IncludeWalker()
        file_scopes = Dict(file => Scope() for file in files)
        for (file, p) in parsed
            s = file_scopes[file]
            s.loc = Location(file, 0)
            S.current_scope = s
            S.loc = s.loc
            CSTAnalyzer.trav(walker, p, s, S)
        end
        for (scope, incl) in walker.include_map
            haskey(file_scopes, incl) || continue
            scope_include!(scope, file_scopes[incl])
        end
        S, file_scopes
    end

end # module


text:
function preconditioner(d, q)
    P_preallocated = Diagonal(zeros(d*(q+1), d*(q+1)))
    Pinv_preallocated = Diagonal(zeros(d*(q+1), d*(q+1)))

    @fastmath @inbounds function P(h)
        @simd for i in 1:d
            @simd for j in 0:q
                P_preallocated[j*d + i,j*d + i] = h^(j-q)
            end
        end
        return P_preallocated
    end

    @fastmath @inbounds function P_inv(h)
        @simd for i in 1:d
            @simd for j in 0:q
                Pinv_preallocated[j*d + i,j*d + i] = h^(q-j)
            end
        end
        return Pinv_preallocated
    end

    return (P=P, P_inv=P_inv)
end

ext: jl

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
# Alias of Model{K,Kₑ,Kₘ, K², K₊₁,K₊₂, AK₊₁,AK₊₂} for 1D TEM Maxwell's equations.
# The type parameters AK₊₁ and AK₊₂ specify device-specific arrays types (e.g., CuArray) and
# are user-defined in the constructor.
const ModelTEM{AK₊₁,AK₊₂} = Model{1,1,1, 1, 2,3, AK₊₁,AK₊₂}

# Convenience constructor
function ModelTEM(grid::Grid; Atype::Type=Array)
    cmpₛ = SInt(3)  # shapes in z-axis
    cmpₑ = SInt(1)  # E-field with x-component
    cmpₘ = SInt(2)  # H-field with y-component

    return ModelTEM{Atype{ComplexF,2},Atype{ComplexF,3}}(;grid, cmpₛ, cmpₑ, cmpₘ)
end

# Assign the material parameters on the grid and smooth them.
function calc_matparams!(mdl::ModelTEM)
    # Take necessary fields of mdl.
    grid = mdl.grid

    isbloch = grid.isbloch
    N = grid.N
    l = grid.l
    σ = grid.σ

    lg = grid.ghosted.l
    τl = grid.ghosted.τl
    ∆τ = grid.ghosted.∆τ

    boundft = mdl.boundft
    gₑ = ft2gt.(EE,boundft)
    gₘ = ft2gt.(HH,boundft)

    εarr = mdl.εarr
    µarr = mdl.μarr

    oind2shp = mdl.oind2shp
    oind2εind = mdl.oind2εind
    oind2μind = mdl.oind2μind

    cmpₑ = mdl.cmpₑ
    cmpₘ = mdl.cmpₘ

    εind2ε = sub_pind2matprm(mdl.εind2ε, cmpₑ)
    μind2μ = sub_pind2matprm(mdl.μind2μ, cmpₘ)

    # Create temporary storages.
    εxx_oind1d = create_oind_array(N)
    μyy_oind1d = create_oind_array(N)

    # Assign the material parameters.
    assign_param!(εarr, tuple(μyy_oind1d), oind2shp, oind2εind, εind2ε, gₑ, τl, isbloch)  # ε tensors (rand-0, so scalars)
    assign_param!(μarr, tuple(εxx_oind1d), oind2shp, oind2μind, μind2μ, gₘ, τl, isbloch)  # μ tensors (rand-0, so scalars)

    # Smooth the material parameters.
    smooth_param!(εarr, tuple(εxx_oind1d), oind2shp, oind2εind, εind2ε, gₑ, l, lg, σ, ∆τ, mdl.ise˔shp)  # ε tensors (rand-0, so scalars)
    smooth_param!(μarr, tuple(μyy_oind1d), oind2shp, oind2μind, μind2μ, gₘ, l, lg, σ, ∆τ, mdl.ish˔shp)  # μ tensors (rand-0, so scalars)

    return nothing
end


text: abstract AbstractScatteredLayerState <: Mocha.LayerState


text:
include("helper_functions.jl")

begin
    λ1 = 1
    λ2 =
    α = 1
    μ(α) = 1 / (1.5 * (α))
    ν(α) = 1 / (0.7 * (α))
    parameters = (λ1, λ2, α, μ, ν)
    opponent_parameters = (λ2, λ1, α, μ, ν)
    p, q = equilibrium_probability(parameters...)
    qs = 0:0.001:1
end
αs = collect(0.5:0.01:3)
λ1 = 0.7 ./ αs
λ2 = 1.2 ./ αs

plot(αs, first.(equilibrium_probability.(λ1, λ2, αs, μ, ν)))


text:
module FlightSims

using Reexport
@reexport using FSimBase
using OrdinaryDiffEq: Tsit5  # ODEProblem default solver
using OrdinaryDiffEq: FunctionMap  # DiscreteProblem default solver
@reexport using FSimZoo

include("APIs/APIs.jl")

end


text:
function find_energy(img)

    energy_x = imfilter(brightness.(img), Kernel.sobel()[2])
    energy_y = imfilter(brightness.(img), Kernel.sobel()[1])

    return sqrt.(energy_x.^2 + energy_y.^2)
end

function find_energy_map(energy)

    sz = size(energy)
    energy_map = zeros(size(sz))
    energy_map[end,:] = energy[end,:]

    next_elements = zeros(Int, sz)

    for i in sz[1]-1:-1:1, j in sz[2]
        left = max(1, j-1)
        right = min(j+1, sz[2])
        local_energy, next_element = findmin(@view energy_map[i+1, left:right])

        # Minimal energy to bottom
        energy_map[i,j] += local_energy + energy[i,j]

        # Element that minimizes the energy
        next_elements[i,j] = next_element - 2

        if left == 1
            next_elements[i,j] += 1
        end
    end

    return energy_map, next_elements
end

"""
Just looks up the minimal path from the `next_elements` array.
"""
function find_seam_at(next_elements, element)

    seam = zeros(Int, size(next_elements)[1])
    seam[1] = element
    for i = 2:length(seam)
        seam[i] = seam[i-1] + next_element[i, seam[i-1]]
    end

    return seam
end

function find_seam(energy)

    energy_map, next_elements = find_energy_map(energy)

    _, min_element = findmin(@view energy_map[1,:])

    return find_seam_at(next_elements, min_element)
end

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
"""
    Manifold
A manifold type. The `Manifold` is used to dispatch to different exponential
and logarithmic maps as well as other function on manifold.
"""
abstract type Manifold end

abstract type AbstractRetractionMethod end

"""
    ExponentialRetraction
Retraction using the exponential map.
"""
struct ExponentialRetraction <: AbstractRetractionMethod end

"""
    retract!(M::Manifold, y, x, v, [t=1], [method::AbstractRetractionMethod=ExponentialRetraction()])
Retraction (cheaper, approximate version of exponential map) of tangent
vector `t*v` at point `x` from manifold `M`.
Result is saved to `y`.
Retraction method can be specified by the last argument. Please look at the
documentation of respective manifolds for available methods.
"""
retract!(M::Manifold, y, x, v, method::ExponentialRetraction) = exp!(M, y, x, v)

retract!(M::Manifold, y, x, v) = retract!(M, y, x, v, ExponentialRetraction())

retract!(M::Manifold, y, x, v, t::Real) = retract!(M, y, x, t*v)

retract!(M::Manifold, y, x, v, t::Real, method::AbstractRetractionMethod) = retract!(M, y, x, t*v, method)
struct Euclidean{T<:Tuple} <: Manifold where {T} end

Euclidean(n::Int) = Euclidean{Tuple{n}}()
Euclidean(m::Int, n::Int) = Euclidean{Tuple{m,n}}()

function representation_size(::Euclidean{Tuple{n}}) where {n}
    return (n,)
end

function representation_size(::Euclidean{Tuple{m,n}}) where {m,n}
    return (m,n)
end

@generated manifold_dimension(::Euclidean{T}) where {T} = *(T.parameters...)

exp!(M::Euclidean, y, x, v) = (y .= x .+ v)

log!(M::Euclidean, v, x, y) = (v .= y .- x)

function zero_tangent_vector!(M::Euclidean, v, x)
    fill!(v, 0)
    return v
end

project_point!(M::Euclidean, x) = x

function project_tangent!(M::Euclidean, w, x, v)
    w .= v
    return w
end
retract(M::Euclidean, x, v, t::Real) = x+t*v

function retract(M::Manifold, x, v, method::AbstractRetractionMethod)
    xr = similar_result(M, retract, x, v)
    retract!(M, xr, x, v, method)
    return xr
end

function retract(M::Manifold, x, v)
    xr = similar_result(M, retract, x, v)
    retract!(M, xr, x, v)
    return xr
end

retract(M::Manifold, x, v, t::Real) = retract(M, x, t*v)

retract(M::Manifold, x, v, t::Real, method::AbstractRetractionMethod) = retract(M, x, t*v, method)


text:
```Write a function that takes a number and returns a list of its digits. So for
 2342 it should return [2,3,4,2].```

 num = 123456

 function split_num(x::Int64)
     digits(x)
 end

 split_num(num)
 digits(num)

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
using Documenter, PtFEM

makedocs(
    format = :html,
    sitename = "PtFEM",
    pages = Any[
        "Introduction" => "INTRO.md",
        "Getting started" => "GETTINGSTARTED.md",
        "Changes w.r.t. PtFEM" => "CHANGES.md",
        "PtFEM.jl documentation" => "index.md",
        "Versions" => "VERSIONS.md",
        "Todo" => "TODO.md",
        "References" => "REFERENCES.md",
    ]
)

deploydocs(
    repo = "github.com/PtFEM/PtFEM.jl.git",
    target = "build",
    julia = "0.6",
    osname = "linux",
    latest = "v0.0.4",
    deps = nothing,
    make = nothing
)


text:
#!/usr/bin/env julia
using Base
using Logging
using Dates
using SHA

# includeList = [ "C:\\Users\\mandi\\OneDrive\\Temp", "C:\\Users\\mandi\\OneDrive\\Scans" ]
includeList = [ "/home/root", "/home/vlado", ]
excludeList = [ Regex("/Camera Roll"), Regex(".lrcat"), Regex(".lrdata"), Regex("/Prints"), Regex("System Volume Information") ]

filesExcluded = 0
filesIncluded = 0
filesSkipped = 0
filesDuplicate = 0
filesError = 0
foldersWalked = 0
files = Dict{}() # structure will be hash of the content plus a tuple of filename/size/ctime/mtime
buf = Vector{UInt8}(undef, 1024*1024)

function dt(t)
    return Dates.format(unix2datetime(t), "mm/dd/YYYY-HH:MM:SS")
end

function checkFile(fileName)
    try
        if filesIncluded % 1000 == 0
#            print(".") # print . for progress every 1000 files
        end
        st = lstat(fileName)
        global filesIncluded += 1
        if (st.size > 64*1024) # only check files larger than 64k
            @debug "checkFile: $fileName"
            fd = open(fileName, "r")
            bytes = readbytes!(fd, buf, 1024*1024, all=false) # read first 1MB of the file to calculate sha2 hash
            close(fd)
            hash = bytes2hex(sha256(buf))
            if haskey(files, hash)
                global filesDuplicate += 1
                @info "duplicate: $fileName size=$(st.size) ctime=$(dt(st.ctime)) mtime=$(dt(st.mtime))"
                @info "original:  $(files[hash].name) size=$(files[hash].size) ctime=$(dt(files[hash].ctime)) mtime=$(dt(files[hash].mtime))"
            end
            global files[hash] = (name = fileName, size = st.size, ctime = st.ctime, mtime = st.mtime)
        else
            global filesSkipped += 1
        end
        return true
    catch err
        global filesError += 1
        @error "$fileName: $err"
        return false
    end
end

function excludeFile(file)
    for regex in excludeList
        if match(regex, file) !== nothing
            global filesExcluded += 1
            return true
        end
    end
    return false
end

function walkDir(dir)
    try
        global foldersWalked += 1
        dirFiles = readdir(dir)
        for fileName in dirFiles
            f = joinpath(dir, fileName)
            if excludeFile(f)
                continue
            end
            if isdir(f)
                walkDir(f)
            elseif isfile(f)
                checkFile(f)
            end
        end
    catch err
        global filesError += 1
        @error "$dir: $err"
    end
end

function startDirScan()
    for dir in includeList
        @info "processing dir: $dir"
        try
            walkDir(dir)
        catch err
            @error "$dir $err"
        end
        @info "finished dir: $dir"
    end
    @info "folders walked: $foldersWalked"
    @info "objects included: $filesIncluded"
    @info "objects excluded: $filesExcluded"
    @info "objects skipped: $filesSkipped"
    @info "objects with duplicates: $filesDuplicate"
    @info "objects with error: $filesError"
end

function init()
#    global logFile = open("vDupes.log", "w+")
#    logger = SimpleLogger(logFile, Logging.Info, Dict{Any,Int64}())
    logger = ConsoleLogger(stdout, Logging.Info; meta_formatter=Logging.default_metafmt, show_limited=true, right_justify=0)
    global_logger(logger)
end

function finish()
#    flush(logFile)
#    close(logFile)
end

function main()
    @info "vdupes starting"
    init()
    time = @elapsed startDirScan()
    speed = round(filesIncluded / time, RoundUp)
    @info "vdupes finished in: $time sec ($speed files/sec)"
    finish()
end

main()

max_stars_repo_name: vladmandic/vdupes

max_issues_repo_name: vladmandic/vdupes

max_forks_repo_name: vladmandic/vdupes


text:
# Programs testing the RPT Instruction.

# Just sandwich an add between RPT instructions.
# If we repeat 10 times, final result should be 10
@asap4asm function repeat_1()
    RPT(10, nop3)
    ADD(dmem[0], 1, bypass[1])
    END_RPT()
end

# Try inserting 2 instructions between start of repeat and end of repeat.
# Use Bypass to to basically interleave these two instructions.
#
# Results in dmem[0] and dmem[1] should both be 10
@asap4asm function repeat_2()
    RPT(10, nop3)
    ADD(dmem[0], 1, bypass[2])
    ADD(dmem[1], 1, bypass[2])
    END_RPT()
end

# Try inserting 3 instructions. Use bypass register 3 to do the same trick
# that we did in repeat_2, but with 3 ADD instrucitons.
@asap4asm function repeat_3()
    RPT(10, nop3)
    ADD(dmem[0], 1, bypass[3])
    ADD(dmem[1], 1, bypass[3])
    ADD(dmem[2], 1, bypass[3])
    END_RPT()
end

# Use the RPT instruction to repeatedly call a function.
@asap4asm function repeat_4()
    # Move 1 into DMEM 0.
    #
    # Our silly little function is going to add this value to itself and then
    # return.
    MOVI(dmem[0], 1)
    RPT(10, nop3)
        BRL(add_dmem_0, j)
    END_RPT()
    BRL(exit)

    @label add_dmem_0
        ADD(dmem[0], 1, dmem[0], nop3)
        BRL(back)

    @label exit
        MOVI(dmem[1], 1)
end

# Jump out of one repeat loop into another repeat loop. Make sure that the
# other repeat loop gets initialized properly.
@asap4asm function repeat_5()
    RPT(5, nop3)
        BRL(second_loop)
    END_RPT()
    BRL(exit)

    @label second_loop
        MOVI(dmem[0], 0, nop3)
        RPT(10, nop3)
            ADD(dmem[0], 1, bypass[1])
        END_RPT()
        BRL(exit)

    @label exit
        MOVI(dmem[1], 10)
end


text:
module test_suite

using DomainSets, BasisFunctions, FrameFun, StaticArrays, FastTransforms
using Test, Printf, LinearAlgebra, Random

FE = FrameFun
BA = BasisFunctions

## Settings

# Test fourier extensions for all parameters
Extensive = false

# Show matrix vector product timings
const show_mv_times = false
const verbose = true
total_mv_allocs = 0
total_mv_time = 0.0

# const include_1d_tests = false
# const include_2d_tests = false
# const include_3d_tests = false
# const include_bigfloat_tests = false
const include_1d_tests = true
const include_2d_tests = true
const include_3d_tests = true
const include_bigfloat_tests = true

########
# Auxiliary functions
########

function delimit(s::AbstractString)
    println()
    println("############")
    println("# ",s)
    println("############")
end

show_mv_timings(basis, domain, fun; options...) = show_mv_times && show_mv_timings(ExtensionFrame(domain, basis), fun; options...)

function show_mv_timings(dict::Dictionary, fun; options...)
    A, B = full_discretization(FrameFun.promote_dictionary(dict, fun), fun; options...)
    op = A
    c1 = zeros(eltype(op), size(src(op)))
    c2 = zeros(eltype(op), size(dest(op)))
    apply!(op, c2, c1)
    t = @timed apply!(op, c2, c1)
    if t[3] > 700
        printstyled("\"\tMV product: ", color=:red)
    elseif t[3] == 0
        printstyled("\"\tMV product: ", color=:green)
    else
        printstyled("\"\tMV product: ", color=:blue)
    end
    println(t[3], " bytes allocated, in ", t[2], " seconds")

    global total_mv_allocs += t[3]
    global total_mv_time += t[2]
    return nothing
end

#######
# Testing -- Accuracy
#######

function test_1d_cases()
    delimit("1D")

    ## The tridiagonal solver for Fourier extensions
    solverstyle = TridiagonalProlateStyle()

    n = 61
    T = 1.9
    domain = -1.0..1.0
    basis = Fourier(n) → -T..T

    println()
    println("Testing \t solver = $solverstyle, \n\t\t Domain = $domain, \n\t\t Basis = Fourier,\n\t\t ELT = Float64 ")
    verbose && println("N\t T\t Complex?\t abserror\t time\t\t memory   ")

    # Complex and real functions / Float64
    f(x) = cos(x^2-1//2)-1
    g(x) = 1im*cos(x^2-1//2)-1
    for func in (f,g)
        F = @timed Fun(func, basis, domain; solverstyle=solverstyle, normalizedsampling=false)
        error = abserror(func, F[1])
        if verbose
            print("$n\t $T\t\t")
            func==g ? print("Y\t") : print("N\t")
            @printf("%3.2e\t %3.2e s\t %3.2e bytes \n",error, F[2],F[3])
        end
        @test  (error < sqrt(eps(Float64))*10)
        show_mv_timings(basis, domain, func; solverstyle=solverstyle, normalizedsampling=false)
    end

    @testset "result" for ELT in (Float32,Float64),
            Basis in (Fourier, ChebyshevT),
            domain in [Interval(ELT(-1.5),ELT(0.7)), DomainSets.UnionDomain(Interval(ELT(-1.5),ELT(-0.5)),Interval(ELT(0.5),ELT(1.5)))],
            solverstyle in (AZStyle(), DirectStyle())

        println()
        println("Testing \t solver = $solverstyle, \n\t\t Domain = $domain, \n\t\t Basis = $(repr(Basis(10))),\n\t\t ELT = $ELT ")
        verbose && println("N\t T\t Complex?\t abserror\t time\t\t memory   ")

        for n in (61,)

            # There is some symmetry around T=2, test smaller and larger values
            for T in (ELT(1.9),)
                for func in (f,g)

                    basis = Basis{ELT}(n) → -T..T
                    F = @timed Fun(func, basis, domain; solverstyle=solverstyle)
                    error = abserror(func, F[1])
                    if verbose
                        print("$n\t $T\t\t")
                        func==g ? print("Y\t") : print("N\t")
                        @printf("%3.2e\t %3.2e s\t %3.2e bytes \n",error, F[2],F[3])
                    end
                    @test  (error < sqrt(eps(ELT))*10)
                    show_mv_timings(basis, domain, func; solverstyle=solverstyle)
                end
            end
        end
    end
end

function test_bigfloat()
    f(x) = cos(x.^2) - big(1.0)
    g(x) = big(1.0)im * cos(x.^2) - big(1.0)
    @testset "result" for Basis in (Fourier, ChebyshevT),
            domain in [Interval(big(-3.0)/2, big(7.0)/10)]
        println()
        println("Testing \t solver = QR_solver\n\t\t Domain = $domain, \n\t\t Basis = $(repr(Basis(10))),\n\t\t ELT = BigFloat ")
        verbose && println("N\t T\t Complex?\t abserror\t time\t\t memory   ")
        for T in (big(17.0)/10,)
            for func in (f,g)
                basis = Basis(91) → -T..T
                F = @timed Fun(func, basis, domain; solverstyle = DirectStyle(), directsolver=:qr)
                error = abserror(func, F[1])
                if verbose
                    @printf("91\t %3.2e\t\t",T)
                    func==g ? print("Y\t") : print("N\t")
                    @printf("%3.2e\t %3.2e s\t %3.2e bytes \n",error, F[2],F[3])
                end
                @test error < 1e-20
                show_mv_timings(basis, domain, func; solverstyle = DirectStyle(), directsolver=:qr)
            end
        end
    end
end

function test_2d_cases()
    delimit("2D")

    f(x,y) = cos(0.5*x)+2*sin(0.2*y)-1.0*x*y
    g(x,y) = 1im*cos(0.5*x)+2*sin(0.2*y)-1.0im*x*y
    @testset "result" for Basis in (Fourier, ChebyshevT),
            domain in [Disk(1.2,SA[-0.1,-0.2]), Rectangle(SA[-1.0,-1.5],SA[0.5,0.7])],
            solverstyle in (AZStyle(), DirectStyle())

        println()
        println("Testing \t solver = $solverstyle \n\t\t Domain = $domain, \n\t\t Basis = $(repr(Basis(10)⊗Basis(10))),\n\t\t ELT = Float64 ")
        verbose && println("N\t\t T\t\t Complex?\t abserror\t time\t\t memory   ")

        for n in ((11,11),)
            for T in ((1.7,2.3),)
                basis = (Basis(n[1]) → -T[1]..T[1]) ⊗ (Basis(n[2]) → -T[2]..T[2])
                for func in (f,g)
                    F = @timed Fun(func, basis, domain; solverstyle=solverstyle)
                    error = abserror(func, F[1])
                    if verbose
                        print("$n \t $T \t\t")
                        func==g ? print("Y\t") : print("N\t")
                        @printf("%3.2e\t %3.2e s\t %3.2e bytes \n",error, F[2],F[3])
                    end
                    @test error < 1e-3
                    show_mv_timings(basis, domain, func; solverstyle=solverstyle)
                end

            end
        end
    end
end

function test_3d_cases()
    delimit("3D")

    f(x,y,z) = cos(x)+sin(y)-x*z
    # @testset "result" for Basis in (Fourier, ChebyshevT), D in (Cube((-1.2,-1.0,-0.9),(1.0,0.9,1.2)),FrameFun.tensorproduct(Interval(-1.0,1.0),Disk(1.05)), FrameFun.Ball(1.2,[-0.3,0.25,0.1])), solver in (AZSolver, )
    #             show(solver); println()
    @testset "result" for Basis in (Fourier, ChebyshevT),
                domain in (Rectangle(SA[-1.2,-1.0,-0.9],SA[1.0,0.9,1.2]), Ball(1.2,SA[-0.3,0.25,0.1])),
                solverstyle in (AZStyle(), )

        println()
        println("Testing \t solver = $solverstyle \n\t\t Domain = $domain, \n\t\t Basis = $(repr(Basis(10)⊗Basis(10)⊗Basis(10))),\n\t\t ELT = Float64 ")
        verbose && println("N\t\t T\t\t Complex?\t abserror\t time\t\t memory   ")

        n = (12,12,12)

        for T in ((1.7,1.2,1.3),)
            basis = (Basis(n[1]) → -T[1]..T[1]) ⊗ (Basis(n[2]) → -T[2]..T[2]) ⊗ (Basis(n[3]) → -T[3]..T[3])
            F = @timed Fun(f, basis, domain; solverstyle=solverstyle, threshold=10.0^(3/4*log10(eps(real(codomaintype(basis))))), oversamplingfactor=1.5)
            error = abs(f(0.4,0.4,0.4)-F[1](0.4,0.4,0.4))
            if verbose
                print("$n \t $T \t\t")
                print("N\t")
                @printf("%3.2e\t %3.2e s\t %3.2e bytes \n",error, F[2],F[3])
            end
            @test  error < 1e-1
            show_mv_timings(basis, domain, f; solverstyle=solverstyle, threshold=10.0^(3/4*log10(eps(real(codomaintype(basis))))), oversamplingfactor=1.5)
        end
    end
end

delimit("Algorithm Implementation and Accuracy")
# FFTW.set_num_threads(Sys.CPU_CORES)

if include_1d_tests
    test_1d_cases()
end

if include_bigfloat_tests
    test_bigfloat()
end

if include_2d_tests
    test_2d_cases()
end

if include_3d_tests
    test_3d_cases()
end

# delimit("Random circles")
# dom = FrameFun.randomcircles(10)
# #    b = Fourier(21) ⊗ Fourier(21)
# b = Fourier(31) ⊗ ChebyshevT(31)
# f(x,y) = cos(20*x+22*y)
# @time F = Fun(f,b,dom,verbose=true)

if show_mv_times
    println("Total bytes in MV products:\t$total_mv_allocs")
    println("Total time in MV products:\t$total_mv_time")
end

end


text:
using LazySets, Polyhedra, Makie, LinearAlgebra

inputSet = Hyperrectangle(low=zeros(3), high=ones(3))

W = [1.8 -4 -2; 1.7 6 -8]

eigenvector = [1.0 0.25 0.4]
#big box
Cb = vcat(Array(Diagonal(ones(3))), Array(-Diagonal(ones(3))))
db = vcat(2*ones(3), 1*ones(3))

#small box / input LazySets
di = vcat(ones(3), zeros(3))

#project plane
C = vcat(Cb, eigenvector, -eigenvector)
d = vcat(db, 0.825, -0.825)
proj_plane = HPolytope(C, d)
#inputSet in front of project plane
infr = HPolytope(vcat(Cb, -eigenvector), vcat(di, -0.825)) #lightpink

#inputSet behind project plane
inbe = HPolytope(vcat(Cb, eigenvector), vcat(di, 0.825)) #lightpink3

#inter proj input
d1 = vcat(di, 0.825, -0.825)
inter = HPolytope(C, d1)

#q1
Cq = vcat(Array(Diagonal(ones(2))), Array(-Diagonal(ones(2))))
dq1 = [1, 6, 0, -5]

#q2
dq2 = [1, 4, 0, -3]

#q3
dq3 = [1, 1, 0, 0]

#p1
Cp = Cq * W

p1if = HPolytope(vcat(Cb, Cp, -eigenvector), vcat(db, dq1, -0.825)) #cyan
p1be = HPolytope(vcat(Cb, Cp, eigenvector), vcat(db, dq1, 0.825)) #cyan3
p1inter = HPolytope(vcat(Cb, Cp, eigenvector, -eigenvector), vcat(db, dq1, 0.825, -0.825))

p2if = HPolytope(vcat(Cb, Cp, -eigenvector), vcat(db, dq2, -0.825))
p2be = HPolytope(vcat(Cb, Cp, eigenvector), vcat(db, dq2, 0.825))
p2inter = HPolytope(vcat(Cb, Cp, eigenvector, -eigenvector), vcat(db, dq2, 0.825, -0.825))
p2in = HPolytope(vcat(Cb, Cp), vcat(di, dq2))

p3if = HPolytope(vcat(Cb, Cp, -eigenvector), vcat(db, dq3, -0.825))
p3be = HPolytope(vcat(Cb, Cp, eigenvector), vcat(db, dq3, 0.825))
p3inter = HPolytope(vcat(Cb, Cp, eigenvector, -eigenvector), vcat(db, dq3, 0.825, -0.825))
p3inf = HPolytope(vcat(Cb, Cp, -eigenvector), vcat(di, dq3, -0.825))
p3inb = HPolytope(vcat(Cb, Cp, eigenvector), vcat(di, dq3, 0.825))

scene1 = plot3d(proj_plane, color=:dodgerblue)

scene1 = plot3d!(inbe, color=:lightpink3)
scene1 = plot3d!(infr, color=:lightpink)
scene1 = plot3d!(inter, color=:dodgerblue4)
scene1 = plot3d!(p1be, color=:cyan3)
scene1 = plot3d!(p1if, color=:cyan)
scene1 = plot3d!(p1inter, color=:red)

scene1 = plot3d(proj_plane, color=:dodgerblue)

scene1 = plot3d!(inbe, color=:lightpink3)
scene1 = plot3d!(infr, color=:lightpink)
scene1 = plot3d!(inter, color=:dodgerblue4)
scene1 = plot3d!(p2be, color=:cyan3)
scene1 = plot3d!(p2if, color=:cyan)
scene1 = plot3d!(p2inter, color=:red)
scene1 = plot3d!(p2in, color=:yellow)

#

scene1 = plot3d(proj_plane, color=:dodgerblue)

scene1 = plot3d!(inbe, color=:lightpink3)
scene1 = plot3d!(p3be, color=:cyan3)
scene1 = plot3d!(p3inb, color=:yellow3)

scene1 = plot3d!(infr, color=:lightpink)
scene1 = plot3d!(inter, color=:dodgerblue4)
scene1 = plot3d!(p3if, color=:cyan)
scene1 = plot3d!(p3inf, color=:yellow)
scene1 = plot3d!(p3inter, color=:red)

#save("plot.png", scene)


text:
module ShallowWater

export SWModel, SWProblem

using StaticArrays
using ..VariableTemplates
using LinearAlgebra: I, dot
using CLIMAParameters.Planet: grav

import ClimateMachine.DGMethods:
    BalanceLaw,
    vars_state_auxiliary,
    vars_state_conservative,
    vars_state_gradient,
    vars_state_gradient_flux,
    vars_integrals,
    flux_first_order!,
    flux_second_order!,
    source!,
    wavespeed,
    boundary_state!,
    compute_gradient_argument!,
    init_state_auxiliary!,
    init_state_conservative!,
    LocalGeometry,
    compute_gradient_flux!

using ..DGMethods.NumericalFluxes

×(a::SVector, b::SVector) = StaticArrays.cross(a, b)
⋅(a::SVector, b::SVector) = StaticArrays.dot(a, b)
⊗(a::SVector, b::SVector) = a * b'

abstract type SWProblem end

abstract type TurbulenceClosure end
struct LinearDrag{L} <: TurbulenceClosure
    λ::L
end
struct ConstantViscosity{L} <: TurbulenceClosure
    ν::L
end

abstract type AdvectionTerm end
struct NonLinearAdvection <: AdvectionTerm end

struct SWModel{PS, P, T, A, S} <: BalanceLaw
    param_set::PS
    problem::P
    turbulence::T
    advection::A
    c::S
end

function vars_state_conservative(m::SWModel, T)
    @vars begin
        U::SVector{3, T}
        η::T
    end
end

function vars_state_auxiliary(m::SWModel, T)
    @vars begin
        f::SVector{3, T}
        τ::SVector{3, T}  # value includes τₒ, g, and ρ
    end
end

function vars_state_gradient(m::SWModel, T)
    @vars begin
        U::SVector{3, T}
    end
end

function vars_state_gradient_flux(m::SWModel, T)
    @vars begin
        ν∇U::SMatrix{3, 3, T, 9}
    end
end

@inline function flux_first_order!(
    m::SWModel,
    F::Grad,
    q::Vars,
    α::Vars,
    t::Real,
)
    FT = eltype(q)
    _grav::FT = grav(m.param_set)
    U = q.U
    η = q.η
    H = m.problem.H

    F.η += U
    F.U += _grav * H * η * I

    advective_flux!(m, m.advection, F, q, α, t)

    return nothing
end

advective_flux!(::SWModel, ::Nothing, _...) = nothing

@inline function advective_flux!(
    m::SWModel,
    A::NonLinearAdvection,
    F::Grad,
    q::Vars,
    α::Vars,
    t::Real,
)
    U = q.U
    H = m.problem.H

    F.U += 1 / H * U ⊗ U

    return nothing
end

function compute_gradient_argument!(
    m::SWModel,
    f::Vars,
    q::Vars,
    α::Vars,
    t::Real,
)
    compute_gradient_argument!(m.turbulence, f, q, α, t)
end

compute_gradient_argument!(::LinearDrag, _...) = nothing

@inline function compute_gradient_argument!(
    T::ConstantViscosity,
    f::Vars,
    q::Vars,
    α::Vars,
    t::Real,
)
    f.∇U = q.U

    return nothing
end

function compute_gradient_flux!(
    m::SWModel,
    σ::Vars,
    δ::Grad,
    q::Vars,
    α::Vars,
    t::Real,
)
    compute_gradient_flux!(m.turbulence, σ, δ, q, α, t)
end

compute_gradient_flux!(::LinearDrag, _...) = nothing

@inline function compute_gradient_flux!(
    T::ConstantViscosity,
    σ::Vars,
    δ::Grad,
    q::Vars,
    α::Vars,
    t::Real,
)
    ν = T.ν
    ∇U = δ.∇U

    σ.ν∇U = ν * ∇U

    return nothing
end

function flux_second_order!(
    m::SWModel,
    G::Grad,
    q::Vars,
    σ::Vars,
    ::Vars,
    α::Vars,
    t::Real,
)
    flux_second_order!(m.turbulence, G, q, σ, α, t)
end

flux_second_order!(::LinearDrag, _...) = nothing

@inline function flux_second_order!(
    ::ConstantViscosity,
    G::Grad,
    q::Vars,
    σ::Vars,
    α::Vars,
    t::Real,
)
    G.U -= σ.ν∇U

    return nothing
end

@inline wavespeed(m::SWModel, n⁻, q::Vars, α::Vars, t::Real) = m.c

@inline function source!(
    m::SWModel{P},
    S::Vars,
    q::Vars,
    diffusive::Vars,
    α::Vars,
    t::Real,
    direction,
) where {P}
    τ = α.τ
    f = α.f
    U = q.U
    S.U += τ - f × U

    linear_drag!(m.turbulence, S, q, α, t)

    return nothing
end

linear_drag!(::ConstantViscosity, _...) = nothing

@inline function linear_drag!(T::LinearDrag, S::Vars, q::Vars, α::Vars, t::Real)
    λ = T.λ
    U = q.U

    S.U -= λ * U

    return nothing
end

function shallow_init_aux! end
function init_state_auxiliary!(m::SWModel, aux::Vars, geom::LocalGeometry)
    shallow_init_aux!(m.problem, aux, geom)
end

function shallow_init_state! end
function init_state_conservative!(m::SWModel, state::Vars, aux::Vars, coords, t)
    shallow_init_state!(m.problem, m.turbulence, state, aux, coords, t)
end

function boundary_state!(
    nf,
    m::SWModel,
    state⁺::Vars,
    aux⁺::Vars,
    n⁻,
    state⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
    _...,
)
    shallow_boundary_state!(
        nf,
        m,
        m.turbulence,
        state⁺,
        aux⁺,
        n⁻,
        state⁻,
        aux⁻,
        t,
    )
end

@inline function shallow_boundary_state!(
    ::RusanovNumericalFlux,
    m::SWModel,
    ::LinearDrag,
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
)
    U⁻ = state⁻.U
    n⁻ = SVector(n⁻)

    state⁺.η = state⁻.η
    state⁺.U = U⁻ - 2 * (n⁻ ⋅ U⁻) * n⁻

    return nothing
end

shallow_boundary_state!(
    ::CentralNumericalFluxGradient,
    m::SWModel,
    ::LinearDrag,
    _...,
) = nothing

shallow_boundary_state!(
    ::CentralNumericalFluxSecondOrder,
    m::SWModel,
    ::LinearDrag,
    _...,
) = nothing

function boundary_state!(
    nf,
    m::SWModel,
    q⁺::Vars,
    σ⁺::Vars,
    α⁺::Vars,
    n⁻,
    q⁻::Vars,
    σ⁻::Vars,
    α⁻::Vars,
    bctype,
    t,
    _...,
)
    shallow_boundary_state!(nf, m, m.turbulence, q⁺, σ⁺, α⁺, n⁻, q⁻, σ⁻, α⁻, t)
end

@inline function shallow_boundary_state!(
    ::RusanovNumericalFlux,
    m::SWModel,
    ::ConstantViscosity,
    q⁺,
    α⁺,
    n⁻,
    q⁻,
    α⁻,
    t,
)
    q⁺.η = q⁻.η
    q⁺.U = -q⁻.U

    return nothing
end

@inline function shallow_boundary_state!(
    ::CentralNumericalFluxGradient,
    m::SWModel,
    ::ConstantViscosity,
    q⁺,
    α⁺,
    n⁻,
    q⁻,
    α⁻,
    t,
)
    q⁺.U = 0

    return nothing
end

@inline function shallow_boundary_state!(
    ::CentralNumericalFluxSecondOrder,
    m::SWModel,
    ::ConstantViscosity,
    q⁺,
    σ⁺,
    α⁺,
    n⁻,
    q⁻,
    σ⁻,
    α⁻,
    t,
)
    q⁺.U = -q⁻.U
    σ⁺.ν∇U = σ⁻.ν∇U

    return nothing
end

end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
"""Validates the contents of a Branch."""
function validate(branch::Branch)::Bool
    is_valid = true

    @debug "Branch validation" branch.name is_valid
    return is_valid
end


text:
using ConstrainedDynamics
using ConstrainedDynamicsVis

# Parameters
joint_axis = [1.0;0.0;0.0]

length1 = 0.5
width, depth = 0.5, 0.5
box1 = Box(width, depth, length1, 1., color = RGBA(1., 1., 0.))
b1 = Box(0.1, 0.1, .1, .1, color = RGBA(1., 0., 0.))

#Corner Vector v
v = [length1 / 2; length1 / 2; -length1 / 2]

# Initial orientation
ϕ1 = 0;
q1 = UnitQuaternion(RotX(ϕ1))

# Links
origin = Origin{Float64}()
link1 = Body(box1)

# Constraints

joint0to1 = EqualityConstraint(Floating(origin, link1))
impact = InequalityConstraint(Impact(link1,[0;0;1.0]; p = v))

eqcs = [joint0to1]
ineqcs = [impact]
shapes = [box1]

mech = Mechanism(origin, [link1], eqcs, ineqcs, shapes = shapes)

setPosition!(link1,x = [0.;-2;1.5])

ωtemp = (rand(3) .- 0.5) * 100
setVelocity!(link1,v = [0;3;7.],ω = ωtemp)

storage = simulate!(mech, 10., record = true)
visualize(mech, storage, shapes)


text:
using Distributions
using NHST

srand(1)

x = randn(20)
y = 1.0 + randn(20)

@assert isapprox(test_statistic(TTest, x, y), -2.8204)
@assert isapprox(p_value(TTest, x, y), 0.00758227)
@assert isapprox(ci(TTest, x, y)[1], -1.5176904)
@assert isapprox(ci(TTest, x, y)[2], -0.2493534)

@assert isapprox(test_statistic(run(TTest, x, y)), -2.8204)
@assert isapprox(p_value(run(TTest, x, y)), 0.00758227)
@assert isapprox(ci(run(TTest, x, y))[1], -1.5176904)
@assert isapprox(ci(run(TTest, x, y))[2], -0.2493534)

run(TTest, x, y)


text:
export TwoDMaxwell

"""
    TwoDMaxwell( mesh, degree)

- mesh       : cartesian mesh
- s_deg_0    : spline degree 0-forms
- s_deg_1    : spline degree 1-forms

"""
struct TwoDMaxwell
    s_deg_0::Int
    s_deg_1::Int

    mesh::TwoDGrid

    mass_line_0::Vector{Array{Float64,1}}
    mass_line_1::Vector{Array{Float64,1}}
    mass_line_mixed::Vector{Array{Float64,1}}
    inv_mass_1::Vector{TwoDLinearSolverSplineMass}
    inv_mass_2::Vector{TwoDLinearSolverSplineMass}
    poisson::TwoDPoisson
    wk1d::Array{Float64,1}
    wk2d::Array{Float64,2}

    function TwoDMaxwell(mesh, degree)
        nx, ny = mesh.nx, mesh.ny
        dx, dy = mesh.dx, mesh.dy
        s_deg_0 = degree
        s_deg_1 = degree - 1

        # Sparse matrices
        # Assemble the mass matrices
        # First assemble a mass line for both degrees

        mass_line_0 = [
            spline_fem_mass_line(s_deg_0) .* dx, spline_fem_mass_line(s_deg_0) .* dy
        ]

        mass_line_1 = [
            spline_fem_mass_line(s_deg_1) .* dx, spline_fem_mass_line(s_deg_1) .* dy
        ]

        mass_line_mixed = [
            spline_fem_mixedmass_line(s_deg_0) .* dx,
            spline_fem_mixedmass_line(s_deg_0) .* dy,
        ]

        # Next put together the 1d parts of the 2d Kronecker product

        eig_values_mass_0_1 = spline_fem_compute_mass_eig(nx, s_deg_0, mass_line_0[1])
        eig_values_mass_0_2 = spline_fem_compute_mass_eig(ny, s_deg_0, mass_line_0[2])
        eig_values_mass_1_1 = spline_fem_compute_mass_eig(nx, s_deg_1, mass_line_1[1])
        eig_values_mass_1_2 = spline_fem_compute_mass_eig(ny, s_deg_1, mass_line_1[2])

        inv_mass_1 = [
            TwoDLinearSolverSplineMass(nx, ny, eig_values_mass_1_1, eig_values_mass_0_2),
            TwoDLinearSolverSplineMass(nx, ny, eig_values_mass_0_1, eig_values_mass_1_2),
            TwoDLinearSolverSplineMass(nx, ny, eig_values_mass_0_1, eig_values_mass_0_2),
        ]

        inv_mass_2 = [
            TwoDLinearSolverSplineMass(nx, ny, eig_values_mass_0_1, eig_values_mass_1_2),
            TwoDLinearSolverSplineMass(nx, ny, eig_values_mass_1_1, eig_values_mass_0_2),
            TwoDLinearSolverSplineMass(nx, ny, eig_values_mass_1_1, eig_values_mass_1_2),
        ]

        poisson = TwoDPoisson(mesh, s_deg_0)

        wk1d = zeros(nx * ny)
        wk2d = zeros(nx, ny)

        return new(
            s_deg_0,
            s_deg_1,
            mesh,
            mass_line_0,
            mass_line_1,
            mass_line_mixed,
            inv_mass_1,
            inv_mass_2,
            poisson,
            wk1d,
            wk2d,
        )
    end
end

function spline_fem_mixedmass_line(deg)
    n = min(3 * deg + 1, 10)

    spline_val_0 = zeros(deg + 1, n)
    spline_val_1 = zeros(deg, n)

    x, w = gausslegendre(n)
    x .+= 1
    x ./= 2
    w ./= 2

    for j in 1:n
        spline_val_0[:, j] .= uniform_bsplines_eval_basis(deg, x[j])
        spline_val_1[:, j] .= uniform_bsplines_eval_basis(deg - 1, x[j])
    end

    mass_line = zeros(2deg)
    for j in 2:(deg + 1), i in j:(deg + 1), k in 1:n
        mass_line[j + deg - 1] += spline_val_0[i, k] * spline_val_1[i - j + 1, k] * w[k]
    end

    for j in (-deg + 1):0, i in 1:(deg + j), k in 1:n
        mass_line[j + deg] += spline_val_0[i, k] * spline_val_1[i - j, k] * w[k]
    end

    return mass_line
end

"""
    compute_rhs_from_function(solver, func, component, form)

Compute the FEM right-hand-side for a given function f and periodic splines of given degree
Its components are ``\\int f N_i dx`` where ``N_i`` is the B-spline starting at ``x_i``
"""
function compute_rhs_from_function(
    solver::TwoDMaxwell, f::Function, component::Int, form::Int
)
    nx, ny = solver.mesh.nx, solver.mesh.ny
    dx, dy = solver.mesh.dx, solver.mesh.dy
    deg_0 = solver.s_deg_0
    deg_1 = solver.s_deg_1

    coefs_dofs = zeros(nx * ny)  # Finite Element right-hand-side

    # Define the spline degree in the 3 dimensions, depending on form and component of the form
    if form == 0
        degree = [deg_0, deg_0]
    elseif (form == 1)
        degree = [deg_0, deg_0]
        component < 3 && (degree[component] = deg_1)
    elseif form == 2
        degree = [deg_1, deg_1]
        component < 3 && (degree[component] = deg_0)
    elseif form == 3
        degree = [deg_1, deg_1]
    else
        @error " Wrong form "
    end

    d1, d2 = degree

    # take enough Gauss points so that projection is exact for splines of degree deg
    # rescale on [0,1] for compatibility with B-splines
    bspl_d1 = zeros(d1 + 1, d1 + 1)
    x1, w1 = gausslegendre(d1 + 1)
    x1 .+= 1
    x1 ./= 2
    w1 ./= 2
    # Compute bsplines at gauss_points
    for k in 1:(d1 + 1)
        bspl_d1[k, :] .= uniform_bsplines_eval_basis(d1, x1[k])
    end

    bspl_d2 = zeros(d2 + 1, d2 + 1)
    x2, w2 = gausslegendre(d2 + 1)
    x2 .+= 1
    x2 ./= 2
    w2 ./= 2
    # Compute bsplines at gauss_points
    for k in 1:(d2 + 1)
        bspl_d2[k, :] .= uniform_bsplines_eval_basis(d2, x2[k])
    end

    counter = 1
    # Compute coefs_dofs = int f(x)N_i(x)
    for i2 in 1:ny, i1 in 1:nx
        coef = 0.0
        # loop over support of B spline
        for j1 in 1:(d1 + 1), j2 in 1:(d2 + 1)
            # loop over Gauss points
            for k1 in 1:(d1 + 1), k2 in 1:(d2 + 1)
                x = dx * (x1[k1] + i1 + j1 - 2)
                y = dy * (x2[k2] + i2 + j2 - 2)
                coef +=
                    w1[k1] * w2[k2] * f(x, y) *
                    bspl_d1[k1, d1 + 2 - j1] *
                    bspl_d2[k2, d2 + 2 - j2]
            end
        end

        # rescale by cell size
        coefs_dofs[counter] = coef * dx * dy
        counter = counter + 1
    end

    return coefs_dofs
end

function compute_e_from_rho!(efield, solver::TwoDMaxwell, rho)
    return compute_e_from_rho!(efield, solver.poisson, rho)
end

"""
    compute_fem_rhs(solver, func , component, form, coefs_dofs)

Compute the FEM right-hand-side for a given function f and periodic splines of given degree
Its components are ``\\int f N_i dx`` where ``N_i`` is the B-spline starting at ``x_i``
- component : Specify the component
- form : Specify 0,1,2 or 3 form
coefs_dofs - Finite Element right-hand-side
"""
function compute_fem_rhs!(
    coefs_dofs::Vector{Float64}, solver::TwoDMaxwell, f, component::Int, form::Int
)
    n1, n2 = solver.mesh.nx, solver.mesh.ny
    δ1, δ2 = solver.mesh.dx, solver.mesh.dy
    s_deg_0, s_deg_1 = solver.s_deg_0, solver.s_deg_1

    if form == 0
        degree = [s_deg_0, s_deg_0]
    elseif form == 1
        degree = [s_deg_0, s_deg_0]
        if component < 3
            degree[component] = s_deg_1
        end
    elseif form == 2
        degree = [s_deg_1, s_deg_1]
        if component < 3
            degree[component] = s_deg_0
        end
    elseif form == 3
        degree = [s_deg_1, s_deg_1]
    end

    d1, d2 = degree

    # take enough Gauss points so that projection is exact for splines of degree deg
    # rescale on [0,1] for compatibility with B-splines

    bspl_d1 = zeros(d1 + 1, d1 + 1)
    x1, w1 = gausslegendre(d1 + 1)
    x1 .+= 1
    x1 ./= 2
    w1 ./= 2
    # Compute bsplines at gauss_points
    for k in 1:(d1 + 1)
        bspl_d1[k, :] .= uniform_bsplines_eval_basis(d1, x1[k])
    end

    bspl_d2 = zeros(d2 + 1, d2 + 1)
    x2, w2 = gausslegendre(d2 + 1)
    x2 .+= 1
    x2 ./= 2
    w2 ./= 2
    # Compute bsplines at gauss_points
    for k in 1:(d2 + 1)
        bspl_d2[k, :] .= uniform_bsplines_eval_basis(d2, x2[k])
    end

    # Compute coefs_dofs = int f(x)N_i(x)
    ind = 0
    for i2 in 1:n2, i1 in 1:n1
        ind += 1
        coef = 0.0
        # loop over support of B spline
        for j1 in 1:(d1 + 1), j2 in 1:(d2 + 1)
            # loop over Gauss points
            for k1 in 1:(d1 + 1), k2 in 1:(d2 + 1)
                xg1 = δ1 * (x1[k1] + i1 + j1 - 2)
                xg2 = δ2 * (x2[k2] + i2 + j2 - 2)
                coef +=
                    w1[k1] * w2[k2] * f(xg1, xg2) *
                    bspl_d1[k1, d1 + 2 - j1] *
                    bspl_d2[k2, d2 + 2 - j2]
            end
        end
        # rescale by cell size
        coefs_dofs[ind] = coef * δ1 * δ2
    end
end

export l2projection
"""
Compute the L2 projection of a given function f on periodic splines of given degree
"""
function l2projection(solver::TwoDMaxwell, f, component, form)
    n1, n2 = solver.mesh.nx, solver.mesh.ny

    compute_fem_rhs!(solver.wk1d, solver, f, component, form)

    if form == 1
        solve(solver.inv_mass_1[component], solver.wk1d)
    elseif form == 2
        solve(solver.inv_mass_2[component], solver.wk1d)
    end
end

function spline_fem_multiply_mass(n_cells, degree, mass, invec)
    outvec = copy(invec)

    ind = 1
    # For the first degree rows we need to put the first part to the back due to periodic boundaries
    for row in 1:degree
        outvec[row] = mass[1] * invec[row]
        for column in 1:(row - 1)
            outvec[row] += mass[column + 1] * (invec[row + column] + invec[row - column])
        end
        for column in row:degree
            outvec[row] +=
                mass[column + 1] * (invec[row + column] + invec[row - column + n_cells])
        end
    end

    for row in (degree + 1):(n_cells - degree)
        outvec[row] = mass[1] * invec[row]
        for column in 1:degree
            outvec[row] += mass[column + 1] * (invec[row + column] + invec[row - column])
        end
    end

    # For the last degree rows, we need to put the second part to the front due to periodic boundaries
    for row in (n_cells - degree + 1):n_cells
        outvec[row] = mass[1] * invec[row]
        for column in 1:(n_cells - row)
            outvec[row] += mass[column + 1] * (invec[row + column] + invec[row - column])
        end
        for column in (n_cells - row + 1):degree
            outvec[row] +=
                mass[column + 1] * (invec[row + column - n_cells] + invec[row - column])
        end
    end

    return outvec
end

"""
    multiply_mass_2dkron!( c_out, solver, mass_line_1, mass_line_2,  c_in )

Multiply by the mass matrix
"""
function multiply_mass_2dkron!(c_out, solver, mass_line_1, mass_line_2, c_in)
    nx1 = solver.mesh.nx
    nx2 = solver.mesh.ny
    deg1 = size(mass_line_1)[1] - 1
    deg2 = size(mass_line_2)[1] - 1

    solver.wk2d .= reshape(c_in, nx1, nx2)

    for j in 1:nx2
        solver.wk2d[:, j] .= spline_fem_multiply_mass(
            nx1, deg1, mass_line_1, solver.wk2d[:, j]
        )
    end
    for i in 1:nx1
        solver.wk2d[i, :] .= spline_fem_multiply_mass(
            nx2, deg2, mass_line_2, solver.wk2d[i, :]
        )
    end

    return c_out .= vec(solver.wk2d)
end

"""
    compute_e_from_b!(e, solver, delta_t, b)

compute Ey from Bz using weak Ampere formulation
"""
function compute_e_from_b!(e, solver, dt, b)
    nx1, nx2 = solver.mesh.nx, solver.mesh.ny
    dx1, dx2 = solver.mesh.dx, solver.mesh.dy

    work = deepcopy(b)

    multiply_mass_2dkron!(
        work[1], solver, solver.mass_line_0[1], solver.mass_line_1[2], b[1]
    )
    multiply_mass_2dkron!(
        work[2], solver, solver.mass_line_1[1], solver.mass_line_0[2], b[2]
    )
    multiply_mass_2dkron!(
        work[3], solver, solver.mass_line_1[1], solver.mass_line_1[2], b[3]
    )

    curl_b = deepcopy(work)

    # Compute curl with periodic boundary conditions
    for j in 1:nx2
        indp2 = j == nx2 ? nx1 * (1 - nx2) : nx1
        for i in 1:nx1
            indp1 = i == nx1 ? 1 - nx1 : 1
            ind2d = (j - 1) * nx1 + i
            ind2d_1 = ind2d + indp1
            ind2d_2 = ind2d + indp2
            curl_b[1][ind2d] = -(work[3][ind2d] - work[3][ind2d_2]) / dx2
            curl_b[2][ind2d] = (work[3][ind2d] - work[3][ind2d_1]) / dx1
            curl_b[3][ind2d] =
                (work[1][ind2d] - work[1][ind2d_2]) / dx2 -
                (work[2][ind2d] - work[2][ind2d_1]) / dx1
        end
    end

    curl_b[1] .= solve(solver.inv_mass_1[1], curl_b[1])
    curl_b[2] .= solve(solver.inv_mass_1[2], curl_b[2])
    curl_b[3] .= solve(solver.inv_mass_1[3], curl_b[3])

    # Update b from solver value
    e[1] .= e[1] .+ dt .* curl_b[1]
    e[2] .= e[2] .+ dt .* curl_b[2]
    return e[3] .= e[3] .+ dt .* curl_b[3]
end

"""
    compute_b_from_e!(b, solver, delta_t, e)

Compute Bz from Ey using strong 1D Faraday equation for spline coefficients
``B_z^{new}(x_j) = B_z^{old}(x_j) - \\frac{\\Delta t}{\\Delta x} (E_y(x_j) - E_y(x_{j-1}) ``

multiplying by discrete curl matrix

"""
function compute_b_from_e!(b, solver, dt, e)
    nx1 = solver.mesh.nx
    nx2 = solver.mesh.ny
    dx1 = solver.mesh.dx
    dx2 = solver.mesh.dy

    for j in 1:nx2
        indp2 = j == 1 ? nx1 * (nx2 - 1) : -nx1
        for i in 1:nx1
            indp1 = i == 1 ? nx1 - 1 : -1
            ind2d = (j - 1) * nx1 + i
            ind2d_1 = ind2d + indp1
            ind2d_2 = ind2d + indp2

            b[1][ind2d] += -dt * (e[3][ind2d] - e[3][ind2d_2]) / dx2
            b[2][ind2d] += dt * (e[3][ind2d] - e[3][ind2d_1]) / dx1
            b[3][ind2d] +=
                -dt *
                ((e[2][ind2d] - e[2][ind2d_1]) / dx1 - (e[1][ind2d] - e[1][ind2d_2]) / dx2)
        end
    end
end

"""
    compute_e_from_j(e, solver, current, component)

Compute E_i from j_i integrated over the time interval using weak Ampere formulation
- solver : Maxwell solver class
- current : Component of the current integrated over time interval
- component : Component of the Efield to be computed
- e : Updated electric field
"""
function compute_e_from_j!(e, solver::TwoDMaxwell, current, component)
    work = solve(solver.inv_mass_1[component], current)

    return e .= e .- work
end

export compute_rho_from_e!

"""
    compute_rho_from_e(rho, solver, efield)

compute rho from e using weak Gauss law ( rho = G^T M_1 e )
"""
function compute_rho_from_e!(rho, solver, efield)
    work = deepcopy(efield)
    multiply_mass_2dkron!(
        work[1], solver, solver.mass_line_1[1], solver.mass_line_0[2], efield[1]
    )
    multiply_mass_2dkron!(
        work[2], solver, solver.mass_line_0[1], solver.mass_line_1[2], efield[2]
    )
    multiply_mass_2dkron!(
        work[3], solver, solver.mass_line_0[1], solver.mass_line_0[2], efield[3]
    )

    dx1, dx2 = solver.mesh.dx, solver.mesh.dy
    nx1, nx2 = solver.mesh.nx, solver.mesh.ny

    # Compute div with periodic boundary conditions
    for j in 1:nx2
        indp2 = j == nx2 ? -nx1 * (nx2 - 1) : nx1
        for i in 1:nx1
            indp1 = i == nx1 ? 1 - nx1 : 1
            ind2d = (j - 1) * nx1 + i
            ind2d_1 = ind2d + indp1
            ind2d_2 = ind2d + indp2
            rho[ind2d] = (
                (work[1][ind2d] - work[1][ind2d_1]) / dx1 +
                (work[2][ind2d] - work[2][ind2d_2]) / dx2
            )
        end
    end

    return rho .*= -1
end

"""
    inner_product( solver, coefs1_dofs, coefs2_dofs, component, form)

- solver : Maxwell solver object
- coefs1_dofs : Coefficient for each DoF
- coefs2_dofs : Coefficient for each DoF
- component : Specify the component
- form : Specify 0,1,2 or 3 form
- r : Result: squared L2 norm
"""
function inner_product(solver, coefs1_dofs, coefs2_dofs, component, form)
    if form == 0
        multiply_mass_2dkron!(
            solver.wk1d, solver, solver.mass_line_0[1], solver.mass_line_0[2], coefs2_dofs
        )
    elseif form == 1
        if component == 1
            multiply_mass_2dkron!(
                solver.wk1d,
                solver,
                solver.mass_line_1[1],
                solver.mass_line_0[2],
                coefs2_dofs,
            )
        elseif component == 2
            multiply_mass_2dkron!(
                solver.wk1d,
                solver,
                solver.mass_line_0[1],
                solver.mass_line_1[2],
                coefs2_dofs,
            )
        elseif component == 3
            multiply_mass_2dkron!(
                solver.wk1d,
                solver,
                solver.mass_line_0[1],
                solver.mass_line_0[2],
                coefs2_dofs,
            )
        end
    elseif form == 2
        if component == 1
            multiply_mass_2dkron!(
                solver.wk1d,
                solver,
                solver.mass_line_0[1],
                solver.mass_line_1[2],
                coefs2_dofs,
            )
        elseif component == 2
            multiply_mass_2dkron!(
                solver.wk1d,
                solver,
                solver.mass_line_1[1],
                solver.mass_line_0[2],
                coefs2_dofs,
            )
        elseif component == 3
            multiply_mass_2dkron!(
                solver.wk1d,
                solver,
                solver.mass_line_1[1],
                solver.mass_line_1[2],
                coefs2_dofs,
            )
        end
    elseif form == 3
        multiply_mass_2dkron!(
            solver.wk1d, solver, solver.mass_line_1[1], solver.mass_line_1[2], coefs2_dofs
        )
    end

    return sum(coefs1_dofs .* solver.wk1d)
end


text:
include("definitions.jl")

function readData(dataFilename::AbstractString, dictionaryFilename::AbstractString)
#read a dataset in non-sparse format. One line per document, word indices (with one-based indexing), space separated.
    f = open(dictionaryFilename);
        dictionary = readlines(f); #an array of strings
    close(f);

    f = open(dataFilename);
    documents = Array{Document}(0);
    for ln in eachline(f)
        splitLn = split(ln);
        wordVector = Array{Int64}(length(splitLn));
        for i = 1:length(splitLn)
           wordVector[i] = parse(Int64, splitLn[i]);
        end
        push!(documents, Document(wordVector));
    end
    close(f);

    documents, dictionary
end


text:
@testset "unimodal" begin

    img = testimage("moonsurface")
    img2 = copy(img)
    img_binary = binarize(UnimodalRosin(), img2)

    # check that the original image has not been changed
    @test img == img2

    # Check that the ! version of the function modifies the parameters
    img3 = copy(img)
    binarize!(UnimodalRosin(),img3)
    @test img3 != img

    # Check that the only entries in the image are 1 or 0.
    non_zeros = findall(x-> x != 0.0 && x != 1.0, img_binary)
    @test length(non_zeros) == 0

    # Check that output and input image are same type
    @test typeof(img_binary) == typeof(img)

    # Check that ones and zeroes have benn assigned to correct sides.
    max, maxpos = findmax(img)
    @test img_binary[maxpos] == 1
    min, minpos = findmin(img)
    @test img_binary[minpos] == 0

end


text:
using HDF5
using Random
using Unitful

include("atmosphere.jl")

"""
    function read_quadrature(fname::String)

Read quarature weights and angles from file. Returns weights, horizontal angle,
azimuthal angle, and number of quadrature points. Quadratures found in
https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/645/A101#/browse
from Bestard & Bueno (2021)
"""
function read_quadrature(fname::String)
    n_points = ""
    switch = false
    for (i, char) in enumerate(fname)
        if char == 'n'
            switch = true
        elseif switch == true
            try parse(Int, char)
                n_points = string(n_points, char)
            catch
                break
            end
        end
    end

    n_points = parse(Int, n_points)

    weights = zeros(n_points)
    θ_array = zeros(n_points)
    ϕ_array = zeros(n_points)

    open(fname, "r") do io
        for (i, line) in enumerate(eachline(fname))
            weights[i] = parse(Float64, split(line)[1])
            θ_array[i] = parse(Float64, split(line)[2])
            ϕ_array[i] = parse(Float64, split(line)[3])
        end
    end

    return weights::AbstractArray, θ_array::AbstractArray, ϕ_array::AbstractArray, n_points::Int
end

"""
    rejection_sampling(n_sites::Int)

Sample 3D vectors according to hydrogen distribution. Samples are obtained
through the acceptance-rejection method. Read
https://towardsdatascience.com/what-is-rejection-sampling-1f6aff92330d
for a nice explanation. Takes number of samples to generate. These are the
voronoi sites.

Make this method better by comparing with another distribution??? This will
increase the rate for acceptance.
Right now the reference distribution is the uniform pdf, scaled between N_H_max
and N_H_min.
"""
function rejection_sampling(n_sites::Int, atmos::Atmosphere, quantity::AbstractArray)
    # Find max and min to convert random number between 0 and 1 to coordinate
    println("---Sampling new sites---")
    z_min = minimum(atmos.z); z_max = maximum(atmos.z)
    x_min = minimum(atmos.x); x_max = maximum(atmos.x)
    y_min = minimum(atmos.y); y_max = maximum(atmos.y)

    Δz = z_max - z_min
    Δx = x_max - x_min
    Δy = y_max - y_min

    # Find max and min populations to scale uniform distribution
    q_min = minimum(quantity)
    q_max = maximum(quantity)

    Δq = q_max - q_min

    # allocate arrays for new sites
    p_vec = Matrix{Float64}(undef, (3, n_sites))u"m"

    for i in 1:n_sites
        print("site $i/$n_sites \r")
        while true
            ref_vec = rand(Float64, 3)
            z_ref = ref_vec[1]*Δz + z_min
            x_ref = ref_vec[2]*Δx + x_min
            y_ref = ref_vec[3]*Δy + y_min

            # acceptance criterion, "reference"
            density_ref = trilinear(z_ref, x_ref, y_ref, atmos, quantity)
            # random sample, compare to reference
            density_ran = rand(Float64)*Δq + q_min
            if density_ref > density_ran
                # a point is accepted, store position and move on
                p_vec[:, i] .= (z_ref, x_ref, y_ref)
                # break to find next site
                break
            end
        end
    end
    print("                                                                 \r")
    return p_vec
end

function rejection_sampling(n_sites::Int, boundaries::Matrix, quantity::AbstractArray)
    # Find max and min to convert random number between 0 and 1 to coordinate
    println("---Sampling new sites---")
    z_min = boundaries[1,1]; z_max = boundaries[1,2]
    x_min = boundaries[2,1]; x_max = boundaries[2,2]
    y_min = boundaries[3,1]; y_max = boundaries[3,2]

    Δz = z_max - z_min
    Δx = x_max - x_min
    Δy = y_max - y_min

    nz = length(quantity[:,1,1])
    nx = length(quantity[1,:,1])
    ny = length(quantity[1,1,:])

    z = Vector{Float64}(undef, nz)
    x = Vector{Float64}(undef, nx)
    y = Vector{Float64}(undef, ny)

    δz = (z_max - z_min)/(nz + 1)
    δx = (x_max - x_min)/(nx + 1)
    δy = (y_max - y_min)/(ny + 1)

    z[1] = z_min
    x[1] = x_min
    y[1] = y_min

    for i in 2:nz
        z[i] = z[i-1] + δz
    end

    for i in 2:nx
        x[i] = x[i-1] + δx
    end

    for i in 2:ny
        y[i] = y[i-1] + δy
    end

    z[end] = z_max
    x[end] = x_max
    y[end] = y_max

    # Find max and min populations to scale uniform distribution
    q_min = minimum(quantity)
    q_max = maximum(quantity)

    Δq = q_max - q_min

    # allocate arrays for new sites
    p_vec = Matrix{Float64}(undef, (3, n_sites))

    for i in 1:n_sites
        print("site $i/$n_sites \r")
        while true
            ref_vec = rand(Float64, 3)
            z_ref = ref_vec[1]*Δz + z_min
            x_ref = ref_vec[2]*Δx + x_min
            y_ref = ref_vec[3]*Δy + y_min

            # acceptance criterion, "reference"
            density_ref = trilinear(z_ref, x_ref, y_ref, z, x, y, quantity)
            # random sample, compare to reference
            density_ran = rand(Float64)*Δq + q_min
            if density_ref > density_ran
                # a point is accepted, store position and move on
                p_vec[:, i] .= (z_ref, x_ref, y_ref)
                # break to find next site
                break
            end
        end
    end
    print("                                                                 \r")
    return p_vec*1u"m"
end

@doc raw"""
    trilinear(x_mrk, y_mrk, z_mrk, hydrogen_populations)

Three-dimensional linear interpolation. Takes a function
$f: \mathbb{R}^3 -> \mathbb{R}$ and returns the trilinear interpolation in the
coordinates (x, y, z). Assumes an original cartesian grid (x, y, z), and
values for each grid point (hydrogen_populations) are defined
"""
function trilinear(z_mrk::Unitful.Length, x_mrk::Unitful.Length, y_mrk::Unitful.Length,
                   atmos::Atmosphere, vals::AbstractArray)
    # Returns the index of the first value in a greater than or equal to x
    # Subtract by 1 to get coordinate of lower corner
    idz = searchsortedfirst(atmos.z, z_mrk) - 1
    idx = searchsortedfirst(atmos.x, x_mrk) - 1
    idy = searchsortedfirst(atmos.y, y_mrk) - 1

    # bounding corner coordinates
    z0 = atmos.z[idz]; z1 = atmos.z[idz+1]
    x0 = atmos.x[idx]; x1 = atmos.x[idx+1]
    y0 = atmos.y[idy]; y1 = atmos.y[idy+1]

    # difference between coordinates and interpolation point
    x_d = (x_mrk - x0)/(x1 - x0)
    y_d = (y_mrk - y0)/(y1 - y0)
    z_d = (z_mrk - z0)/(z1 - z0)

    # values at each corner (z is first index in data array)
    c000 = vals[idz, idx, idy]
    c010 = vals[idz, idx, idy+1]
    c100 = vals[idz, idx+1, idy]
    c110 = vals[idz, idx+1, idy+1]
    c001 = vals[idz+1, idx, idy]
    c011 = vals[idz+1, idx, idy+1]
    c101 = vals[idz+1, idx+1, idy]
    c111 = vals[idz+1, idx+1, idy+1]

    # interpolate in x direction
    c00 = c000*(1 - x_d) + c100*x_d
    c01 = c001*(1 - x_d) + c101*x_d
    c10 = c010*(1 - x_d) + c110*x_d
    c11 = c011*(1 - x_d) + c111*x_d

    # interpolate in y direction
    c0 = c00*(1 - y_d) + c10*y_d
    c1 = c01*(1 - y_d) + c11*y_d

    # intepolate in z direction
    c = c0*(1 - z_d) + c1*z_d
    return c
end

function trilinear(z_mrk::Unitful.Length, x_mrk::Unitful.Length, y_mrk::Unitful.Length,
                   z::Vector{<:Unitful.Length}, x::Vector{<:Unitful.Length}, y::Vector{<:Unitful.Length},
                   vals::AbstractArray)
    # Returns the index of the first value in a greater than or equal to x
    # Subtract by 1 to get coordinate of lower corner
    idz = searchsortedfirst(z, z_mrk) - 1
    idx = searchsortedfirst(x, x_mrk) - 1
    idy = searchsortedfirst(y, y_mrk) - 1

    # bounding corner coordinates
    z0 = atmos.z[idz]; z1 = atmos.z[idz+1]
    x0 = atmos.x[idx]; x1 = atmos.x[idx+1]
    y0 = atmos.y[idy]; y1 = atmos.y[idy+1]

    # difference between coordinates and interpolation point
    x_d = (x_mrk - x0)/(x1 - x0)
    y_d = (y_mrk - y0)/(y1 - y0)
    z_d = (z_mrk - z0)/(z1 - z0)

    # values at each corner (z is first index in data array)
    c000 = vals[idz, idx, idy]
    c010 = vals[idz, idx, idy+1]
    c100 = vals[idz, idx+1, idy]
    c110 = vals[idz, idx+1, idy+1]
    c001 = vals[idz+1, idx, idy]
    c011 = vals[idz+1, idx, idy+1]
    c101 = vals[idz+1, idx+1, idy]
    c111 = vals[idz+1, idx+1, idy+1]

    # interpolate in x direction
    c00 = c000*(1 - x_d) + c100*x_d
    c01 = c001*(1 - x_d) + c101*x_d
    c10 = c010*(1 - x_d) + c110*x_d
    c11 = c011*(1 - x_d) + c111*x_d

    # interpolate in y direction
    c0 = c00*(1 - y_d) + c10*y_d
    c1 = c01*(1 - y_d) + c11*y_d

    # intepolate in z direction
    c = c0*(1 - z_d) + c1*z_d
    return c
end

@doc raw"""
    bilinear(x_mrk, y_mrk, x_bounds, y_bounds, vals)

Two-dimensional linear interpolation. Takes a function
$f: \mathbb{R}^2 -> \mathbb{R}$ and returns the biilinear interpolation in the
coordinates (x_mrk, y_mrk). Assumes an underlying cartesian grid (x_i, y_i)
Values are defined on the corners of the rectangle spanned by x_bounds and
y_bounds. x_mrk and y_mrk have to lie inside this rectangle.
"""
function bilinear(x_mrk::Unitful.Length, y_mrk::Unitful.Length,
                  x_bounds::Tuple, y_bounds::Tuple,
                  vals::AbstractMatrix)

    # corner coordinates
    x1 = x_bounds[1]; x2 = x_bounds[2]
    y1 = y_bounds[1]; y2 = y_bounds[2]

    # function value for each corner. 11 is lower left corner, 22 is upper right
    Q11 = vals[1,1]     # (x1, y1)
    Q12 = vals[1,2]     # (x1, y2)
    Q21 = vals[2,1]     # (x2, y1)
    Q22 = vals[2,2]     # (x2, y2)

    # Rectangle side length
    dx = x2 - x1
    dy = y2 - y1

    # Interpolate in x-direction
    f1 = ((x2 - x_mrk)*Q11 + (x_mrk - x1)*Q21)/dx
    f2 = ((x2 - x_mrk)*Q12 + (x_mrk - x1)*Q22)/dx

    # Interpolate in y-direction
    f = ((y2 - y_mrk)*f1 + (y_mrk - y1)*f2)/dy
    return f
end

function bilinear(x_mrk::Float64, y_mrk::Float64,
                  x_bounds::Tuple, y_bounds::Tuple,
                  vals::AbstractMatrix)

    # corner coordinates
    x1 = x_bounds[1]; x2 = x_bounds[2]
    y1 = y_bounds[1]; y2 = y_bounds[2]

    # function value for each corner. 11 is lower left corner, 22 is upper right
    Q11 = vals[1,1]     # (x1, y1)
    Q12 = vals[1,2]     # (x1, y2)
    Q21 = vals[2,1]     # (x2, y1)
    Q22 = vals[2,2]     # (x2, y2)

    # Rectangle side length
    dx = x2 - x1
    dy = y2 - y1

    # Interpolate in x-direction
    f1 = ((x2 - x_mrk)*Q11 + (x_mrk - x1)*Q21)/dx
    f2 = ((x2 - x_mrk)*Q12 + (x_mrk - x1)*Q22)/dx

    # Interpolate in y-direction
    f = ((y2 - y_mrk)*f1 + (y_mrk - y1)*f2)/dy
    return f
end

"""
    trapezoidal(Δx::AbstractFloat, a::AbstractFloat, b::AbstractFloat)

Trapezoidal rule.
"""
function trapezoidal(Δx::Unitful.Length, a::PerLength, b::PerLength)
    area = Δx*(a + b)/2
    return area
end

"""
    xy_intersect(ϕ::AbstractFloat)

Finds x and y direction ray is moving in defined by the azimuthal angle ϕ
"""
function xy_intersect(ϕ::Float64)
    local sign_x, sign_y
    @assert 0 <= ϕ < 2π "Bad angle. Expects ϕ ∈ [0, 2π), but ϕ = $ϕ"
    if 0 <= ϕ < π/2
        # 1st quadrant. Negative x, negative y
        sign_x = -1
        sign_y = -1
    elseif π/2 < ϕ < π
        # 2nd quadrant. Positive x, negative y
        sign_x = 1
        sign_y = -1
    elseif π < ϕ < 3π/2
        # 3rd quadrant. Positive x, positive y
        sign_x = 1
        sign_y = 1
    elseif 3π/2 < ϕ < 2π
        # 4th quadrant. Negative x, positive y
        sign_x = -1
        sign_y = 1
    end
    return sign_x::Int, sign_y::Int
end

"""
    xy_intersect(k::Vector{Float64})

Finds x and y direction ray is moving in defined by direction vector k
"""
function xy_intersect(k::Vector{Float64})
    local sign_x, sign_y
    @assert norm(k) ≈ 1 "Bad direction vector. Expects length to be 1"
    if k[2] > 0 && k[3] > 0
        # 1st quadrant. Negative x, negative y
        sign_x = -1
        sign_y = -1
    elseif k[2] < 0 && k[3] > 0
        # 2nd quadrant. Positive x, negative y
        sign_x = 1
        sign_y = -1
    elseif k[2] < 0 && k[3] < 0
        # 3rd quadrant. Positive x, positive y
        sign_x = 1
        sign_y = 1
    elseif k[2] > 0 && k[3] < 0
        # 4th quadrant. Negative x, positive y
        sign_x = -1
        sign_y = 1
    else
        # theta is 0 or 180, ray either straight up or down, doesn't really matter what
        # the sign is here, interpolation is exact anyways as the ray hits the
        # grid points...
        sign_x = 1
        sign_y = 1
    end
    return sign_x::Int, sign_y::Int
end

"""
    function range_bounds(sign::Int, bound::Int)

Given a quadrant from sign_x and sign_x from xy_intersect(), this function
determines the loop start and stop point for the short characteristics rays.
"""
function range_bounds(sign::Int, bound::Int)
    if sign == 1
        start = 2
        stop = bound-1
    elseif sign == -1
        start = bound-1
        stop = 2
    end
    return start::Int, stop::Int
end

"""
    linear_weights(Δτ::AbstractFloat)

Compute weights for linear integration of Source function in formal solution,
also return weight for I_0. Approximates expressions with Taylor expansion for
large or small optical depth
"""
function linear_weights(Δτ::AbstractFloat)
    local α, β, expΔτ
    if Δτ < 5e-4
        expΔτ = 1 - Δτ + 0.5*Δτ^2
        α = Δτ*(1/2 - Δτ/3)
        β = Δτ*(1/2 - Δτ/6)
    elseif Δτ > 50
        expΔτ = 0.0
        α = 1/Δτ
        β = 1.0 - α
    else
        expΔτ = exp(-Δτ)
        α = (1 - expΔτ)/Δτ - expΔτ
        β = 1 - α - expΔτ
    end
    return α, β, expΔτ
end

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
abstract type DataOneV2 <: DataRepo end

export KNB, TERN
include("KNB.jl")
include("TERN.jl")

description(repo::DataOneV2, mainpage) = try desc_(repo, mainpage) catch missing end

author(repo::DataOneV2, mainpage) = try authors_(repo, mainpage) catch missing end

published_date(repo::DataOneV2, mainpage) = try pub_date(repo, mainpage) catch missing end

function license(repo::DataOneV2, mainpage)
    try
        license_ele = eachmatch(sel"intellectualRights", mainpage.root)
        if length(license_ele) > 0
            text_only(first(license_ele))
        else
            ""
        end
    catch
        missing
    end
end

function paper_title_(repo::DataOneV2, mainpage)
    try
        paper_title_ele = eachmatch(sel"title", mainpage.root)
        if length(paper_title_ele) > 0
            text_only(first(paper_title_ele))
        else
            ""
        end
    catch
        missing
    end
end


text:
# Initialize data structures in element container
function init_elements!(elements, mesh::P4estMesh{3}, basis::LobattoLegendreBasis)
  @unpack node_coordinates, jacobian_matrix,
          contravariant_vectors, inverse_jacobian = elements

  calc_node_coordinates!(node_coordinates, mesh, basis.nodes)

  for element in 1:ncells(mesh)
    calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis)

    calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix,
                                node_coordinates, basis)

    calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix, basis)
  end

  return nothing
end

# Interpolate tree_node_coordinates to each quadrant at the nodes of the specified basis
function calc_node_coordinates!(node_coordinates,
                                mesh::P4estMesh{3},
                                basis::LobattoLegendreBasis)
  # Hanging nodes will cause holes in the mesh if its polydeg is higher
  # than the polydeg of the solver.
  @assert length(basis.nodes) >= length(mesh.nodes) "The solver can't have a lower polydeg than the mesh"

  calc_node_coordinates!(node_coordinates, mesh, basis.nodes)
end

# Interpolate tree_node_coordinates to each quadrant at the specified nodes
function calc_node_coordinates!(node_coordinates,
                                mesh::P4estMesh{3},
                                nodes::AbstractVector)
  # Macros from p4est
  p4est_root_len = 1 << P4EST_MAXLEVEL
  p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l)

  trees = unsafe_wrap_sc(p8est_tree_t, mesh.p4est.trees)

  for tree in eachindex(trees)
    offset = trees[tree].quadrants_offset
    quadrants = unsafe_wrap_sc(p8est_quadrant_t, trees[tree].quadrants)

    for i in eachindex(quadrants)
      element = offset + i
      quad = quadrants[i]

      quad_length = p4est_quadrant_len(quad.level) / p4est_root_len

      nodes_out_x = 2 * (quad_length * 1/2 * (nodes .+ 1) .+ quad.x / p4est_root_len) .- 1
      nodes_out_y = 2 * (quad_length * 1/2 * (nodes .+ 1) .+ quad.y / p4est_root_len) .- 1
      nodes_out_z = 2 * (quad_length * 1/2 * (nodes .+ 1) .+ quad.z / p4est_root_len) .- 1

      matrix1 = polynomial_interpolation_matrix(mesh.nodes, nodes_out_x)
      matrix2 = polynomial_interpolation_matrix(mesh.nodes, nodes_out_y)
      matrix3 = polynomial_interpolation_matrix(mesh.nodes, nodes_out_z)

      multiply_dimensionwise!(
        view(node_coordinates, :, :, :, :, element),
        matrix1, matrix2, matrix3,
        view(mesh.tree_node_coordinates, :, :, :, :, tree)
      )
    end
  end

  return node_coordinates
end

# Initialize node_indices of interface container
@inline function init_interface_node_indices!(interfaces::P4estInterfaceContainer{3},
                                              faces, orientation, interface_id)
  # Iterate over primary and secondary element
  for side in 1:2
    # Align interface at the primary element (primary element has surface indices (:i, :j)).
    # The secondary element needs to be indexed differently.
    if side == 1
      surface_index1 = :i
      surface_index2 = :j
    else
      surface_index1, surface_index2 = orientation_to_indices_p4est(faces[2], faces[1], orientation)
    end

    if faces[side] == 0
      # Index face in negative x-direction
      interfaces.node_indices[side, interface_id] = (:one, surface_index1, surface_index2)
    elseif faces[side] == 1
      # Index face in positive x-direction
      interfaces.node_indices[side, interface_id] = (:end, surface_index1, surface_index2)
    elseif faces[side] == 2
      # Index face in negative y-direction
      interfaces.node_indices[side, interface_id] = (surface_index1, :one, surface_index2)
    elseif faces[side] == 3
      # Index face in positive y-direction
      interfaces.node_indices[side, interface_id] = (surface_index1, :end, surface_index2)
    elseif faces[side] == 4
      # Index face in negative z-direction
      interfaces.node_indices[side, interface_id] = (surface_index1, surface_index2, :one)
    else # faces[side] == 5
      # Index face in positive z-direction
      interfaces.node_indices[side, interface_id] = (surface_index1, surface_index2, :end)
    end
  end

  return interfaces
end

# Initialize node_indices of boundary container
@inline function init_boundary_node_indices!(boundaries::P4estBoundaryContainer{3},
                                             face, boundary_id)
  if face == 0
    # Index face in negative x-direction
    boundaries.node_indices[boundary_id] = (:one, :i, :j)
  elseif face == 1
    # Index face in positive x-direction
    boundaries.node_indices[boundary_id] = (:end, :i, :j)
  elseif face == 2
    # Index face in negative y-direction
    boundaries.node_indices[boundary_id] = (:i, :one, :j)
  elseif face == 3
    # Index face in positive y-direction
    boundaries.node_indices[boundary_id] = (:i, :end, :j)
  elseif face == 4
    # Index face in negative z-direction
    boundaries.node_indices[boundary_id] = (:i, :j, :one)
  else # face == 5
    # Index face in positive z-direction
    boundaries.node_indices[boundary_id] = (:i, :j, :end)
  end

  return boundaries
end

# Initialize node_indices of mortar container
# faces[1] is expected to be the face of the small side.
@inline function init_mortar_node_indices!(mortars::P4estMortarContainer{3},
                                           faces, orientation, mortar_id)
  for side in 1:2
    # Align mortar at small side.
    # The large side needs to be indexed differently.
    if side == 1
      surface_index1 = :i
      surface_index2 = :j
    else
      surface_index1, surface_index2 = orientation_to_indices_p4est(faces[2], faces[1], orientation)
    end

    if faces[side] == 0
      # Index face in negative x-direction
      mortars.node_indices[side, mortar_id] = (:one, surface_index1, surface_index2)
    elseif faces[side] == 1
      # Index face in positive x-direction
      mortars.node_indices[side, mortar_id] = (:end, surface_index1, surface_index2)
    elseif faces[side] == 2
      # Index face in negative y-direction
      mortars.node_indices[side, mortar_id] = (surface_index1, :one, surface_index2)
    elseif faces[side] == 3
      # Index face in positive y-direction
      mortars.node_indices[side, mortar_id] = (surface_index1, :end, surface_index2)
    elseif faces[side] == 4
      # Index face in negative z-direction
      mortars.node_indices[side, mortar_id] = (surface_index1, surface_index2, :one)
    else # faces[side] == 5
      # Index face in positive z-direction
      mortars.node_indices[side, mortar_id] = (surface_index1, surface_index2, :end)
    end
  end

  return mortars
end

# Convert p4est orientation code to node indices.
# Return node indices that index "my side" wrt "other side",
# i.e., i and j are indices of other side.
function orientation_to_indices_p4est(my_face, other_face, orientation_code)
  # my_face and other_face are the face directions (zero-based)
  # of "my side" and "other side" respectively.
  # Face corner 0 of the face with the lower face direction connects to a corner of the other face.
  # The number of this corner is the orientation code in p4est.
  lower = my_face <= other_face

  # x_pos, y_neg, and z_pos are the directions in which the face has right-handed coordinates
  # when looked at from the outside.
  my_right_handed = my_face in (1, 2, 5)
  other_right_handed = other_face in (1, 2, 5)

  # If both or none are right-handed when looked at from the outside, they will have different
  # orientations when looked at from the same side of the interface.
  flipped = my_right_handed == other_right_handed

  # In the folowing illustrations, p4est's face corner numbering is shown.
  # ξ and η are the local coordinates of the respective face.
  # We're looking at both faces from the same side of the interface, so that "other side"
  # (in the illustrations on the left) has right-handed coordinates.
  if !flipped
    if orientation_code == 0
      # Corner 0 of other side matches corner 0 of my side
      #   2┌──────┐3   2┌──────┐3
      #    │      │     │      │
      #    │      │     │      │
      #   0└──────┘1   0└──────┘1
      #     η            η
      #     ↑            ↑
      #     │            │
      #     └───> ξ      └───> ξ
      surface_index1 = :i
      surface_index2 = :j
    elseif ((lower && orientation_code == 2) # Corner 0 of my side matches corner 2 of other side
        || (!lower && orientation_code == 1)) # Corner 0 of other side matches corner 1 of my side
      #   2┌──────┐3   0┌──────┐2
      #    │      │     │      │
      #    │      │     │      │
      #   0└──────┘1   1└──────┘3
      #     η            ┌───> η
      #     ↑            │
      #     │            ↓
      #     └───> ξ      ξ
      surface_index1 = :j_backwards
      surface_index2 = :i
    elseif ((lower && orientation_code == 1) # Corner 0 of my side matches corner 1 of other side
        || (!lower && orientation_code == 2)) # Corner 0 of other side matches corner 2 of my side
      #   2┌──────┐3   3┌──────┐1
      #    │      │     │      │
      #    │      │     │      │
      #   0└──────┘1   2└──────┘0
      #     η                 ξ
      #     ↑                 ↑
      #     │                 │
      #     └───> ξ     η <───┘
      surface_index1 = :j
      surface_index2 = :i_backwards
    else # orientation_code == 3
      # Corner 0 of my side matches corner 3 of other side and
      # corner 0 of other side matches corner 3 of my side.
      #   2┌──────┐3   1┌──────┐0
      #    │      │     │      │
      #    │      │     │      │
      #   0└──────┘1   3└──────┘2
      #     η           ξ <───┐
      #     ↑                 │
      #     │                 ↓
      #     └───> ξ           η
      surface_index1 = :i_backwards
      surface_index2 = :j_backwards
    end
  else # flipped
    if orientation_code == 0
      # Corner 0 of other side matches corner 0 of my side
      #   2┌──────┐3   1┌──────┐3
      #    │      │     │      │
      #    │      │     │      │
      #   0└──────┘1   0└──────┘2
      #     η            ξ
      #     ↑            ↑
      #     │            │
      #     └───> ξ      └───> η
      surface_index1 = :j
      surface_index2 = :i
    elseif orientation_code == 2
      # Corner 0 of my side matches corner 2 of other side and
      # corner 0 of other side matches corner 2 of my side.
      #   2┌──────┐3   0┌──────┐1
      #    │      │     │      │
      #    │      │     │      │
      #   0└──────┘1   2└──────┘3
      #     η            ┌───> ξ
      #     ↑            │
      #     │            ↓
      #     └───> ξ      η
      surface_index1 = :i
      surface_index2 = :j_backwards
    elseif orientation_code == 1
      # Corner 0 of my side matches corner 1 of other side and
      # corner 0 of other side matches corner 1 of my side.
      #   2┌──────┐3   3┌──────┐2
      #    │      │     │      │
      #    │      │     │      │
      #   0└──────┘1   1└──────┘0
      #     η                 η
      #     ↑                 ↑
      #     │                 │
      #     └───> ξ     ξ <───┘
      surface_index1 = :i_backwards
      surface_index2 = :j
    else # orientation_code == 3
      # Corner 0 of my side matches corner 3 of other side and
      # corner 0 of other side matches corner 3 of my side.
      #   2┌──────┐3   2┌──────┐0
      #    │      │     │      │
      #    │      │     │      │
      #   0└──────┘1   3└──────┘1
      #     η           η <───┐
      #     ↑                 │
      #     │                 ↓
      #     └───> ξ           ξ
      surface_index1 = :j_backwards
      surface_index2 = :i_backwards
    end
  end

  return surface_index1, surface_index2
end


text:
type PlaceOccupied <: Exception end

type Position
     x::Int
     y::Int
end

abstract Zwierze

type NULL <: Zwierze end

type Drapieznik <: Zwierze
     name::String
     attack::Float16
     position::Position
end

type Ofiara <: Zwierze
     name::String
     defense::Float16
     position::Position
end

type Swiat
     width::Int
     height::Int
     mapa
end

function initWorld(world::Swiat)
	 for i = 1:n
	     for j = 1:n
	         oneWorld.mapa[i,j] = NULL()
	     end
	 end
end

function printWorld(world::Swiat)
	 println("Swiat: D - drapieznik, O - ofiara")
	 for i = 1:world.width
	     for j = 1:world.height
	         if (world.mapa[i,j] == NULL())
		     print("N")
		 elseif (typeof(world.mapa[i,j]) == Ofiara)
		     print("O")
		 elseif (typeof(world.mapa[i,j]) == Drapieznik)
		     print("D")
		 end
		 print(" ")
	     end
	     println()
	 end
end

function addAnimal(animal::Zwierze, world::Swiat)
	 if world.mapa[animal.position.x, animal.position.y] != NULL()
	     throw(PlaceOccupied)
	 else
	     world.mapa[animal.position.x, animal.position.y] = animal
	 end
end

function seek(animal::Zwierze, world::Swiat)
	 for i = 1:world.width
	     for j = 1:world.height
	     	 if world.mapa[i,j] == NULL()
		    return Position(i,j)
		 end
	     end
	 end

	 return Position(0,0)
end

function odleglosc{T1 <: Zwierze, T2 <: Zwierze}(a1::T1, a2::T2)
	 return abs(a1.position.x - a2.position.x) + abs(a1.position.y - a2.position.y)
end

function interakcja(o::Ofiara, d::Drapieznik, world::Swiat)
  local newPositionVictim = seek(o,world)
  local newPositionPredator = o.position

  world.mapa[d.position.x, d.position.y] = NULL()
  world.mapa[o.position.x, o.position.y] = d
  world.mapa[newPositionVictim.x, newPositionVictim.y] = o

  d.position = newPositionPredator
  o.position = newPositionVictim

end

function interakcja(d::Drapieznik, o::Ofiara, world::Swiat)
  local newPositionVictim = o.position

  world.mapa[d.position.x,d.position.y] = NULL()
  world.mapa[o.position.x,o.position.y] = d

  d.position = newPositionVictim
  o = NULL()
end

function interakcja(d1::Drapieznik, d2::Drapieznik, world::Swiat)
	 "Wrrrr"
end

function interakcja(o1::Ofiara, o2::Ofiara, world::Swiat)
	 "Beeee"
end

#run section

n = 3

oneWorld = Swiat(n,n,Array(Zwierze,n,n))

initWorld(oneWorld)

predator = Drapieznik("lew",1.32,Position(1,1))
predator2 = Drapieznik("nosorożec",1.22,Position(1,2))

victim = Ofiara("antylopa",0.59,Position(2,3))
victim2 = Ofiara("zebra",1.09,Position(3,3))

addAnimal(predator,oneWorld)
addAnimal(predator2,oneWorld)
#addAnimal(predator2,oneWorld)

addAnimal(victim,oneWorld)
addAnimal(victim2,oneWorld)

printWorld(oneWorld)

println("Odleglosc w metryce taksowkowej miedzy (1,1) a (1,2): ",odleglosc(predator,victim))

println(interakcja(predator,predator2,oneWorld))
println(interakcja(victim, victim2, oneWorld))

println("Interakcja ofiara-drapieznik")
interakcja(victim,predator,oneWorld)
printWorld(oneWorld)

println("Interakcja drapieznik-ofiara")
interakcja(predator,victim,oneWorld)
printWorld(oneWorld)

max_stars_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_issues_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_forks_repo_name: wchmielarz/NobodyDoesItLikeJuliet


text:
"""
controls mask visibility responds to keyboard input
"""
module KeyboardVisibility
using ModernGL, ..DisplayWords, ..StructsManag, Setfield, ..PrepareWindow,   ..DataStructs , Rocket, GLFW,Dictionaries,  ..ForDisplayStructs, ..TextureManag,  ..OpenGLDisplayUtils,  ..Uniforms, Match, Parameters,DataTypesBasic
export processKeysInfo, setVisAndRender

"""
processing information from keys - the instance of this function will be chosen on
the basis mainly of multiple dispatch
"""
function processKeysInfo(textSpecObs::Identity{TextureSpec{T}},actor::SyncActor{Any, ActorWithOpenGlObjects},keyInfo::KeyboardStruct )where T
    textSpec =  textSpecObs.value
    if(keyInfo.isCtrlPressed)
        setVisAndRender(false,actor.actor,textSpec.uniforms )
        @info " set visibility of $(textSpec.name) to false"
        #to enabling undoing it
        addToforUndoVector(actor, ()-> setVisAndRender(true,actor.actor,textSpec.uniforms )    )
    elseif(keyInfo.isShiftPressed)
        setVisAndRender(true,actor.actor,textSpec.uniforms )
        @info " set visibility of $(textSpec.name) to true"
       #to enabling undoing it
       addToforUndoVector(actor, ()->setVisAndRender(false,actor.actor,textSpec.uniforms )   )

    elseif(keyInfo.isAltPressed)
        oldTex = actor.actor.textureToModifyVec
        actor.actor.textureToModifyVec= [textSpec]
        @info " set texture for manual modifications to  $(textSpec.name)"
       if(!isempty(oldTex))
       addToforUndoVector(actor, ()->begin  @info actor.actor.textureToModifyVec=[oldTex[1]] end)
       end
    end #if
end #processKeysInfo

"""
sets  visibility and render the result to the screen
"""
function setVisAndRender(isVis::Bool,actor::ActorWithOpenGlObjects,unifs::TextureUniforms )
    setTextureVisibility(isVis,unifs )
    basicRender(actor.mainForDisplayObjects.window)

end#setVisAndRender

end#module


text:
reload("gencond.jl")

function genba( kappa, m, n, kind )
  kappa = log10(kappa)
  sqrt_kappa = sqrt(kappa)
  V = qr(randn(m,m))[1]
  D = logspace(0, -kappa, m)
  B = scale(V, D)*V'
  W = qr(randn(n,n))[1]
  if kind == 1
    A = scale(V[:,m-n+1:m], logspace(0,-sqrt_kappa,n))*W'
  elseif kind == 2
    A = scale(V[:,1:n], logspace(0,-sqrt_kappa,n))*W'
  elseif kind == 3
    A = scale([V[:,m-cld(n,2)+1:m] V[:,1:fld(n,2)]],logspace(0,-sqrt_kappa,n))*W'
  elseif kind == 4
    A = gencond(10^sqrt_kappa, m, n)
  else
    U = [V[:,m-cld(n,2)+1:m] V[:,1:fld(n,2)]]
    S = [D[m-cld(n,2)+1:m]; D[1:fld(n,2)]]
    A = scale(U,1./sqrt(S))*W'
  end
  return sparse(triu(B)),A
end


text:
using ParquetFiles
using Parquet
using IteratorInterfaceExtensions
using TableTraits
using Test

@testset "ParquetFiles" begin

parquet_pkg_dir = joinpath(dirname(pathof(Parquet)), "..")

include(joinpath(parquet_pkg_dir, "test", "get_parcompat.jl"))

test_filename = joinpath(parquet_pkg_dir, "test", "parquet-compatibility", "parquet-testdata", "impala", "1.1.1-NONE", "nation.impala.parquet")

pqf = load(test_filename)

@test IteratorInterfaceExtensions.isiterable(pqf) == true
@test TableTraits.isiterabletable(pqf) == true

it = IteratorInterfaceExtensions.getiterator(pqf)

ar = collect(it)

@test length(ar) == 25
@test ar[1] == (n_nationkey = 0, n_name = "ALGERIA", n_regionkey = 0, n_comment = " haggle. carefully final deposits detect slyly agai")

@test sprint((stream,data)->show(stream, "text/html", data), pqf)[1:100] == "

n_nationkey	n_name	n_regionkey	n_comment
show(stream, "application/vnd.dataresource+json", data), pqf)[1:100] == "{\"schema\":{\"fields\":[{\"name\":\"n_nationkey\",\"type\":\"integer\"},{\"name\":\"n_name\",\"type\":\"string\"},{\"nam"
@test sprint(show, pqf)[1:100] == "25x4 Parquet file\nn_nationkey │ n_name    │ n_regionkey\n────────────┼─"
@test showable("text/html", pqf) == true
@test showable("application/vnd.dataresource+json", pqf) == true

end


text:
using Temporal, Dates

function find_index_col(fields)::Int
    regex = r"date|time|index"i
    @inbounds for (j, field) in enumerate(fields)
        if occursin(regex, String(field))
            return j
        end
    end
    return 0
end

function TS(named_tuple::NamedTuple, index_element::Int=find_index_col(keys(named_tuple)))::TS
    index = named_tuple[index_element]
    array = zeros(Float64, (length(index), length(named_tuple)-1))
    value_elements = setdiff(1:length(named_tuple), index_element)
    fields = [f for f in keys(named_tuple)[value_elements]]
    @inbounds for (array_col, tuple_element) in enumerate(value_elements)
        array[:,array_col] = Real.(named_tuple[tuple_element])
    end
    return TS(array, index, fields)
end