text:
"""
    topofile, topotype, ntopo = topodata("simlation/path/_output"::AbstractString)
    topofile, topotype, ntopo = topodata("simlation/path/_output/topo.data"::AbstractString)

read topo.data
"""
function topodata(outdir::AbstractString)
    # filename
    filename = occursin("topo.data", basename(outdir)) ? outdir : joinpath(outdir,"topo.data")
    ## check args
    if !isfile(filename); error("$filename is not found."); end;

    ## read ascdata
    f = open(filename,"r")
    ascdata = readlines(f)
    close(f)

    baseline = findfirst(x->occursin("ntopofiles", x), ascdata)
    ntopo = parse(Int64, split(ascdata[baseline], r"\s+", keepempty=false)[1])

    if ntopo == 1
        topofile = replace(ascdata[baseline+2], r"[\'\s]" => "")
        topotype = parse(Int64, split(ascdata[baseline+3], r"\s+", keepempty=false)[1])
    else
        # preallocate
        topofile = Vector{String}(undef, ntopo)
        topotype = Vector{Int64}(undef, ntopo)
        # filename
        for i = 1:ntopo
            topofile[i] = replace(ascdata[baseline-1+3i], r"[\'\s]" => "")
            topotype[i] = parse(Int64, split(ascdata[baseline+3i], r"\s+", keepempty=false)[1])
        end
    end

    # return
    return topofile, topotype, ntopo
end
#################################

#################################
"""
    dtopofile, dtopotype, ndtopo = dtopodata("simlation/path/_output"::AbstractString)
    dtopofile, dtopotype, ndtopo = dtopodata("simlation/path/_output/dtopo.data"::AbstractString)

read dtopo.data
"""
function dtopodata(outdir::AbstractString)
    # filename
    filename = occursin("dtopo.data", basename(outdir)) ? outdir : joinpath(outdir,"dtopo.data")
    ## check args
    if !isfile(filename); error("$filename is not found."); end;

    # read
    f = open(filename,"r")
    ascdata = readlines(f)
    close(f)

    baseline = findfirst(x->occursin("mdtopofiles", x), ascdata)
    ndtopo = parse(Int64, split(ascdata[baseline], r"\s+", keepempty=false)[1])

    if ndtopo == 0; println("No mdtopofile");  return nothing, nothing, ndtopo

    elseif ndtopo == 1
        dtopofile = replace(ascdata[baseline+3], r"[\'\s]" => "")
        dtopotype = parse(Int64, split(ascdata[baseline+4], r"\s+", keepempty=false)[1])

    else
        # preallocate
        dtopofile = Vector{String}(undef, ndtopo)
        dtopotype = Vector{Int64}(undef, ndtopo)
        # filename
        for i = 1:ndtopo
            dtopofile[i] = replace(ascdata[baseline+3i], r"[\'\s]" => "")
            dtopotype[i] = parse(Int64, split(ascdata[baseline+3i+1], r"\s+", keepempty=false)[1])
        end
    end

    # return
    return dtopofile, dtopotype, ndtopo
end
#################################

#################################
"""
    bathtopo = loadtopo(outdir::AbstractString)
    bathtopo = loadtopo(filename::AbstractString, topotype=3::Integer)

load topography data
"""
function loadtopo(filename::AbstractString, topotype=3::Integer)
    ## from _output directory
    if isdir(filename)
        topofile, topotype, ntopo = VisClaw.topodata(filename)
        return VisClaw.loadtopo.(topofile, topotype)
    end

    ## check args
    isfile(filename) || error("file $filename is not found.")
    any(topotype .== [2,3,4]) || error("unsupported topotype")

    ## NetCDF
    ( filename[end-2:end] == ".nc" && topotype != 4 ) && (topotype=4)
    if topotype == 4
        var_x, var_y, var_z = VisClaw.getvarname_nctopo(filename)
        x = NetCDF.ncread(filename, var_x)
        y = NetCDF.ncread(filename, var_y)
        topo = permutedims(NetCDF.ncread(filename, var_z), [2,1])

        bathtopo = VisClaw.Topo(length(x), length(y), x, y, mean(diff(x)), mean(diff(x)), topo)
        return bathtopo
    end

    ## separator in regular expression
    regex = r"([+-]?(?:\d+\.?\d*|\.\d+)(?:[eE][+-]?\d+)?)"
    ## open topofile
    f = open(filename,"r")
        ## read header
        ncols = parse(Int64, match(regex, readline(f)).captures[1])
        nrows = parse(Int64, match(regex, readline(f)).captures[1])
        xll = parse(Float64, match(regex, readline(f)).captures[1])
        yll = parse(Float64, match(regex, readline(f)).captures[1])
        cellsize = parse(Float64, match(regex, readline(f)).captures[1])
        nodata = match(regex, readline(f)).captures[1]
        ## read topography in asc
        dataorg = readlines(f)
    ## close topofile
    close(f)

    ## meshgrid
    x = collect(Float64, LinRange(xll, round(xll+(ncols-1)*cellsize, digits=3), ncols))
    y = collect(Float64, LinRange(yll, round(yll+(nrows-1)*cellsize, digits=3), nrows))

    # check topotype
    tmp = replace(dataorg[1], r"^\s+|,?\s+$" => "")
    tmp = replace(tmp, "," => " ") # for csv data
    tmp = split(tmp, r"\s+",keepempty=false)
    tmp = parse.(Float64, tmp)
    # topotype 2?
    if length(tmp) == 1
        if topotype == 3; println("topotype 2?"); end
        topotype = 2
    end
    # topotype 3?
    if length(tmp) > 1
        if topotype == 2; println("topotype 3?"); end
        topotype = 3
    end

    ## assign topography
    if topotype == 2
        topo = parse.(Float64, dataorg)
        topo = reshape(topo, (ncols, nrows))
        topo = permutedims(topo,[2 1])
    elseif topotype == 3
        topo = zeros(nrows, ncols)
        for k = 1:nrows
            line = replace(dataorg[k], r"^\s+|,?\s+$" => "")
            line = replace(line, "," => " ") # for csv data
            line = split(line, r"\s+",keepempty=false)
            topo[k,:] = parse.(Float64, line)
        end
    end
    topo[topo.==nodata] .= NaN ## replace nodate to NaN
    topo = reverse(topo, dims=1) ## flip
    bathtopo = VisClaw.Topo(ncols, nrows, x, y, cellsize, cellsize, topo)

    return bathtopo
end
#################################

#########################################
"""
    dtopo = loaddtopo(outdir::AbstractString)
    dtopo = loaddtopo(filename::AbstractString, topotype=3::Integer)
    dtopo = loaddeform(outdir::AbstractString)
    dtopo = loaddeform(filename::AbstractString, topotype=3::Integer)

load spatial distribution of seafloor deformation (dtopo)
"""
function loaddeform(filename::AbstractString, topotype=3::Integer)
    ## from _output directory
    if isdir(filename)
        dtopofile, topotype, ntopo = VisClaw.dtopodata(filename)
        return VisClaw.loaddeform.(dtopofile, topotype)
    end

    ## check args
    if !isfile(filename); error("file $filename is not found."); end
    if (topotype!=2) & (topotype!=3); error("Invalid topotype"); end

    ## separator in regular expression
    regex = r"([+-]?(?:\d+\.?\d*|\.\d+)(?:[eE][+-]?\d+)?)"
    ## open topofile
    f = open(filename,"r")
        ## read header
        mx = parse(Int64, match(regex, readline(f)).captures[1])
        my = parse(Int64, match(regex, readline(f)).captures[1])
        mt = parse(Int64, match(regex, readline(f)).captures[1])
        xlow = parse(Float64, match(regex, readline(f)).captures[1])
        ylow = parse(Float64, match(regex, readline(f)).captures[1])
        t0 = parse(Float64, match(regex, readline(f)).captures[1])
        dx = parse(Float64, match(regex, readline(f)).captures[1])
        dy = parse(Float64, match(regex, readline(f)).captures[1])
        dt = parse(Float64, match(regex, readline(f)).captures[1])
        ## read topography in asc
        dataorg = readlines(f)
    ## close topofile
    close(f)

    ## meshgrid
    x = collect(Float64, LinRange(xlow, round(xlow+(mx-1)*dx, digits=3), mx))
    y = collect(Float64, LinRange(ylow, round(ylow+(my-1)*dy, digits=3), my))

    # check topotype
    tmp = replace(dataorg[1], r"^\s+|,?\s+$" => "") # equivalent to strip?
    tmp = replace(tmp, "," => " ") # for csv data
    tmp = split(tmp, r"\s+", keepempty=false)
    tmp = parse.(Float64, tmp)
    # topotype 2?
    if length(tmp) == 1
        if topotype == 3; println("topotype 2?"); end
        topotype = 2
    end
    # topotype 3?
    if length(tmp) > 1
        if topotype == 2; println("topotype 3?"); end
        topotype = 3
    end

    ## assign topography
    if topotype == 2
        deform = parse.(Float64, dataorg)
        deform = reshape(topo, (mx, my, mt))
        deform = permutedims(topo,[2 1 3])
    elseif topotype == 3
        deform = zeros(my, mx, mt)
        for k = 1:mt
            for i = 1:my
                line = replace(dataorg[i+(k-1)my], r"^\s+|,?\s+$" => "")
                line = replace(line, "," => " ") # for csv data
                line = split(line, r"\s+",keepempty=false)
                deform[i,:,k] = parse.(Float64, line)
            end
        end
    end
    if mt==1; deform = dropdims(deform; dims=3); end
    deform = reverse(deform, dims=1) ## flip
    #deform[abs.(deform).<1e-2] .= NaN
    dtopo = VisClaw.DTopo(mx,my,x,y,dx,dy,mt,t0,dt,deform)

    return dtopo
end
#########################################
const loaddtopo = loaddeform
#########################################


text:
using GLVisualize, GLWindow, GLAbstraction, AxisArrays

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0f0,1f0,360))
end

description = """
Maximum intensity volume Plot.
"""

function volume_data(N)
    volume = Float32[sin(x/15.0)+sin(y/15.0)+sin(z/15.0) for x=1:N, y=1:N, z=1:N]
    max = maximum(volume)
    min = minimum(volume)
    volume = AxisArray((volume .- min) ./ (max .- min), (:y, :x, :z), (1, 1, 3))
end

volumedata = if isfile(joinpath(homedir(), "brain.nii"))
    using NIfTI # If we have, use some nice example data
    niread(joinpath(homedir(), "brain.nii")).raw  # TODO: make it an AxisArray
else
    volume_data(128)
end

volume = visualize(volumedata, :mip)

_view(volume, window)

if !isdefined(:runtests)
    renderloop(window)
end

ext: jl


text:
module Sys

export  CPU_CORES,
        OS_NAME,
        WORD_SIZE,
        ARCH,
        MACHINE,
        cpu_info,
        cpu_summary,
        uptime,
        loadavg,
        free_memory,
        total_memory,
        dlext,
	shlib_ext,
        dllist,
        dlpath

import ..Base: WORD_SIZE, OS_NAME, ARCH, MACHINE
import ..Base: show, uv_error

global CPU_CORES

function init_sysinfo()
    # set CPU core count
    global const CPU_CORES = int(
        haskey(ENV,"JULIA_CPU_CORES") ?
        ENV["JULIA_CPU_CORES"] :
        ccall(:jl_cpu_cores, Int32, ())
    )
    global const SC_CLK_TCK = ccall(:jl_SC_CLK_TCK, Clong, ())
end

type UV_cpu_info_t
    model::Ptr{Uint8}
    speed::Int32
    cpu_times!user::Uint64
    cpu_times!nice::Uint64
    cpu_times!sys::Uint64
    cpu_times!idle::Uint64
    cpu_times!irq::Uint64
end
type CPUinfo
    model::ASCIIString
    speed::Int32
    cpu_times!user::Uint64
    cpu_times!nice::Uint64
    cpu_times!sys::Uint64
    cpu_times!idle::Uint64
    cpu_times!irq::Uint64
    CPUinfo(model,speed,u,n,s,id,ir)=new(model,speed,u,n,s,id,ir)
end
CPUinfo(info::UV_cpu_info_t) = CPUinfo(bytestring(info.model), info.speed,
    info.cpu_times!user, info.cpu_times!nice, info.cpu_times!sys,
    info.cpu_times!idle, info.cpu_times!irq)

function show(io::IO, info::CPUinfo, header::Bool=true, prefix::String="    ")
    tck = SC_CLK_TCK
    if header
        println(io, info.model, ": ")
        print(io, " "^length(prefix))
        if tck > 0
            @printf io "    %5s    %9s    %9s    %9s    %9s    %9s\n" "speed" "user" "nice" "sys" "idle" "irq"
        else
            @printf io "    %5s    %9s  %9s  %9s  %9s  %9s ticks\n" "speed" "user" "nice" "sys" "idle" "irq"
        end
    end
    print(io, prefix)
    if tck > 0
        @printf io "%5d MHz  %9d s  %9d s  %9d s  %9d s  %9d s" info.speed info.cpu_times!user/tck info.cpu_times!nice/tck info.cpu_times!sys/tck info.cpu_times!idle/tck info.cpu_times!irq/tck
    else
        @printf io "%5d MHz  %9d  %9d  %9d  %9d  %9d ticks" info.speed info.cpu_times!user info.cpu_times!nice info.cpu_times!sys info.cpu_times!idle info.cpu_times!irq
    end
end

function _cpu_summary(io::IO, cpu::Array{CPUinfo}, i, j)
    if j-i < 9
        header = true
        for x = i:j
            if header == false println(io) end
            show(io,cpu[x],header,"#$(x-i+1) ")
            header = false
        end
    else
        summary = CPUinfo(cpu[i].model,0,0,0,0,0,0)
        count = j-i+1
        for x = i:j
            summary.speed += cpu[i].speed
            summary.cpu_times!user += cpu[x].cpu_times!user
            summary.cpu_times!nice += cpu[x].cpu_times!nice
            summary.cpu_times!sys += cpu[x].cpu_times!sys
            summary.cpu_times!idle += cpu[x].cpu_times!idle
            summary.cpu_times!irq += cpu[x].cpu_times!irq
        end
        summary.speed = div(summary.speed,count)
        show(io,summary,true,"#1-$(count) ")
    end
    println(io)
end

function cpu_summary(io::IO=STDOUT, cpu::Array{CPUinfo}=cpu_info())
    model = cpu[1].model
    first = 1
    for i = 2:length(cpu)
        if model != cpu[i].model
            _cpu_summary(io,cpu,first,i-1)
            first = i
        end
    end
    _cpu_summary(io,cpu,first,length(cpu))
end

function cpu_info()
    UVcpus = Array(Ptr{UV_cpu_info_t},1)
    count = Array(Int32,1)
    uv_error("uv_cpu_info",ccall(:uv_cpu_info, Int32, (Ptr{Ptr{UV_cpu_info_t}}, Ptr{Int32}), UVcpus, count))
    cpus = Array(CPUinfo,count[1])
    for i = 1:length(cpus)
        cpus[i] = CPUinfo(unsafe_load(UVcpus[1],i))
    end
    ccall(:uv_free_cpu_info, Void, (Ptr{UV_cpu_info_t}, Int32), UVcpus[1], count[1])
    cpus
end

function uptime()
    uptime_ = Array(Float64,1)
    uv_error("uv_uptime",ccall(:uv_uptime, Int32, (Ptr{Float64},), uptime_))
    return uptime_[1]
end

function loadavg()
    loadavg_ = Array(Float64,3)
    ccall(:uv_loadavg, Void, (Ptr{Float64},), loadavg_)
    return loadavg_
end

free_memory() = ccall(:uv_get_free_memory, Uint64, ())
total_memory() = ccall(:uv_get_total_memory, Uint64, ())

if OS_NAME === :Darwin
    const dlext = "dylib"
elseif OS_NAME === :Windows
    const dlext = "dll"
else
    #assume OS_NAME === :Linux, or similar
    const dlext = "so"
end

# This is deprecated!  use dlext instead!
const shlib_ext = dlext

@linux_only begin
    immutable dl_phdr_info
        # Base address of object
        addr::Cuint

        # Null-terminated name of object
        name::Ptr{Uint8}

        # Pointer to array of ELF program headers for this object
        phdr::Ptr{Void}

        # Number of program headers for this object
        phnum::Cshort
    end

    # This callback function called by dl_iterate_phdr() on Linux
    function dl_phdr_info_callback( di_ptr::Ptr{dl_phdr_info}, size::Csize_t, dynamic_libraries_ptr::Ptr{Array{String,1}} )
        di = unsafe_load(di_ptr)

        # Skip over objects without a path (as they represent this own object)
        name = bytestring(di.name)
        if !isempty(name)
            dynamic_libraries = unsafe_pointer_to_objref( dynamic_libraries_ptr )
            push!(dynamic_libraries, name )
        end
        convert(Cint, 0)::Cint
    end
end #@linux_only

function dllist()
    dynamic_libraries = Array(String,0)

    @linux_only begin
        const callback = cfunction(dl_phdr_info_callback, Cint, (Ptr{dl_phdr_info}, Csize_t, Ptr{Array{String,1}} ))
        ccall( cglobal("dl_iterate_phdr"), Cint, (Ptr{Void}, Ptr{Void}), callback, pointer_from_objref(dynamic_libraries) )
    end

    @osx_only begin
        numImages = ccall( cglobal("_dyld_image_count"), Cint, (), )

        # start at 1 instead of 0 to skip self
        for i in 1:numImages-1
            name = bytestring(ccall( cglobal("_dyld_get_image_name"), Ptr{Uint8}, (Uint32,), uint32(i)))
            push!(dynamic_libraries, name)
        end
    end

    dynamic_libraries
end

function dlpath( handle::Ptr{Void} )
    return bytestring(ccall( :jl_pathname_for_handle, Ptr{Uint8}, (Ptr{Void},), handle ))
end

function dlpath( libname::String )
    handle = dlopen(libname)
    path = dlpath( handle )
    dlclose(handle)
    return path
end

end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
module Regression

using Senuk.Solvers.gd

type RegressionModel
  weights :: Vector{Float64}
end

abstract Regularization
type Ridge <: Regularization
  l2 :: Float64
end

type Lasso <: Regularization
  l1 :: Float64
end

function create(X :: Matrix{Float64}, y :: Vector{Float64}; step_size = 0.1,
                max_iter = 100, verbose = false, regularization = Regularization)
  if typeof(regularization) == Ridge
    # minimize Ridge cost function
    error("not implemented yet")
  end

  if  typeof(regularization) == Lasso
    # minimize Lasso cost function
    error("not implemented yet")
  end

  # regression without a regularization term
  X = [ones(size(y)[1])  X]
  w = gd(δf, X, y,  step_size, max_iter, verbose)
  lr = RegressionModel(w)
  return lr
end

function δf(w :: Vector{Float64}, X :: Matrix{Float64}, y :: Vector{Float64})
  cost = (y - X*w)'*(y - X*w)
  Δcost = -2*X'*(y -X*w)
  return (cost[1], Δcost)
end

function δf_ridge(w :: Vector{Float64}, X :: Matrix{Float64}, y :: Vector{Float64})
  cost = (y - X*w)'*(y - X*w)
  Δcost = -2*X'*(y -X*w)
  return (cost[1], Δcost)
end

function δf_lasso(w :: Vector{Float64}, X :: Matrix{Float64}, y :: Vector{Float64})
  cost = (y - X*w)'*(y - X*w)
  Δcost = -2*X'*(y -X*w)
  return (cost[1], Δcost)
end

function predict(model :: RegressionModel, X :: Matrix{Float64})
  X = [ones(size(X)[1]) X]
  return X*model.weights
end

end

max_stars_repo_name: upul/Senuk

max_issues_repo_name: upul/Senuk

max_forks_repo_name: upul/Senuk


text:
using CSyntax
using Quaternions
using GLTF

@static if Sys.isapple()
    const VERSION_MAJOR = 4
    const VERSION_MINOR = 1
end

include(joinpath(@__DIR__, "glutils.jl"))
include(joinpath(@__DIR__, "camera.jl"))

# init window
width, height = fb_width, fb_height = 640, 480
window = startgl(width, height)

glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LESS)

# we use glTF file format because it's OpenGL friendly and we don't need to care about
# how to parse and extract mesh data from other mesh file format like .obj format.
# load glTF file
sphere = GLTF.load(joinpath(@__DIR__, "sphere.gltf"))
sphere_data = [read(joinpath(@__DIR__, b.uri)) for b in sphere.buffers]

# load sphere position metadata
search_name(x, keyword) = x[findfirst(x->occursin(keyword, x.name), x)]
pos = search_name(sphere.accessors, "position")
pos_bv = sphere.bufferViews[pos.bufferView]
# load sphere index metadata
indices = search_name(sphere.accessors, "indices")
idx_bv = sphere.bufferViews[indices.bufferView]

# create VBO and EBO
vbo = GLuint(0)
@c glGenBuffers(1, &vbo)
glBindBuffer(pos_bv.target, vbo)
glBufferData(pos_bv.target, pos_bv.byteLength, C_NULL, GL_STATIC_DRAW)
pos_data = sphere_data[pos_bv.buffer]
@c glBufferSubData(pos_bv.target, 0, pos_bv.byteLength, &pos_data[pos_bv.byteOffset])

ebo = GLuint(0)
@c glGenBuffers(1, &ebo)
glBindBuffer(idx_bv.target, ebo)
glBufferData(idx_bv.target, idx_bv.byteLength, C_NULL, GL_STATIC_DRAW)
idx_data = sphere_data[idx_bv.buffer]
@c glBufferSubData(idx_bv.target, 0, idx_bv.byteLength, &idx_data[idx_bv.byteOffset])

# create VAO
vao = GLuint(0)
@c glGenVertexArrays(1, &vao)
glBindVertexArray(vao)
glBindBuffer(pos_bv.target, vbo)
glVertexAttribPointer(0, 3, pos.componentType, pos.normalized, pos_bv.byteStride, Ptr{Cvoid}(pos.byteOffset))
glEnableVertexAttribArray(0)

# load and compile shaders from file
vert_shader = createshader(joinpath(@__DIR__, "quat.vert"), GL_VERTEX_SHADER)
frag_shader = createshader(joinpath(@__DIR__, "quat.frag"), GL_FRAGMENT_SHADER)

# link program
shader_prog = createprogram(vert_shader, frag_shader)

# enable cull face
glEnable(GL_CULL_FACE)
glCullFace(GL_BACK)
glFrontFace(GL_CCW)
# set background color to gray
glClearColor(0.2, 0.2, 0.2, 1.0)

# camera
camera = PerspectiveCamera()
setposition!(camera, [0.0, 0.0, 5.0])
model_loc = glGetUniformLocation(shader_prog, "model")
view_loc = glGetUniformLocation(shader_prog, "view")
proj_loc = glGetUniformLocation(shader_prog, "proj")
glUseProgram(shader_prog)
glUniformMatrix4fv(view_loc, 1, GL_FALSE, get_view_matrix(camera))
glUniformMatrix4fv(proj_loc, 1, GL_FALSE, get_projective_matrix(window, camera))

# spheres in world
sphere_world = GLfloat[-2.0 0.0  0.0;
                        2.0 0.0  0.0;
                       -2.0 0.0 -2.0;
                        1.5 1.0 -1.0]
model_mats = map(1:4) do i
    GLfloat[ 1.0 0.0 0.0 sphere_world[i,1];
             0.0 1.0 0.0 sphere_world[i,2];
             0.0 0.0 1.0 sphere_world[i,3];
             0.0 0.0 0.0       1.0        ]
end

let
updatefps = FPSCounter()
while !GLFW.WindowShouldClose(window)
    updatefps(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, GLFW.GetFramebufferSize(window)...)
    # drawing
    glUseProgram(shader_prog)
    glBindVertexArray(vao)
    glBindBuffer(idx_bv.target, ebo)
    for i = 1:4
        glUniformMatrix4fv(model_loc, 1, GL_FALSE, model_mats[i])
        glDrawElements(GL_TRIANGLES, indices.count, indices.componentType, Ptr{Cvoid}(0))
    end
    # check and call events
    GLFW.PollEvents()
    # move camera
    updatecamera!(window, camera)
    glUniformMatrix4fv(view_loc, 1, GL_FALSE, get_view_matrix(camera))
    # swap the buffers
    GLFW.SwapBuffers(window)
end
end # let

GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
"""
    STL{T <: Real}
"""
struct STL{T <: Real}
    trend::Vector
    seasonal::Vector
    remainder::Vector
end

function nextodd(x::T)::Integer where T <: Real
    cx = ceil(Integer, x)
    return isodd(cx) ? cx+1 : cx
end

function sma(x, n::Integer; center::Bool = true)
    N = length(x)
    @assert 1 <= n & n <= N
    if n == 1
        return x
    end
    res = Vector{Any}(missing,N)

    # initial and final value positions
    ivp = findfirst(!ismissing, x)
    fvp = findlast(!ismissing, x)

    # using missing values to center ma
    a = center ? div(n,2) : 0

    # initial moving average value
    ma = sum(x[ivp:n+ivp-1])/n
    res[a+ivp] = ma
    for i in 1:(N-n-ivp-(N-fvp)+1)
        resai = ma + (x[n+ivp+i-1] - x[ivp+i-1])/n
        # missing values are imputed
        res[a+ivp+i] = ismissing(resai) ? ma : resai
        ma = res[a+ivp+i]
    end
    return res
end

"""
    stl(Yv, np; robust=false,
                nl=nextodd(np),
                ns=10*length(Yv)+1,
                nt=nextodd(1.5*np/(1-1.5/ns)),
                ni=robust ? 1 : 2,
                no=0,
                spm=false,
                qsmp=max(div(np,7),2))

Decompose a time series into trend, seasonal, and remainder components.
"STL has a simple design that consists of a sequence of applications of the loess smoother; the simplicity
allows analysis of the properties of the procedure and allows fast computation, even for very long time series
and large amounts of trend and seasonal smoothing. Other features of STL  are specification of amounts of seasonal
and trend smoothing that range, in a nearly continous way, from very small amount of smoothing to a very large
amount; robust estimates of the trend and seasonal components that are not distorted by aberrant behavior in
the data; specification of the period of the seasonal component to any intenger multiple of the time sampling
interval greater than one; and the ability to decompose time series with missing values."*

All default values are chosen following the recommendations of the original paper when those were recommended.
`ns` is recommended to be chosen of the basis of knowledge of the time series and on the basis of diagnostic
methods; it must nonethelessbe always odd and at least 7. A default value is not advised on the original
paper, instead the same default value used in the stl implementation in R in usere here.
for `no` the authors advise 5 ("safe value") or 10 ("near certainty of convergence") cycles or a
convergence criterion when robustness is required, in this case when `robust` is true computations
stop when convergence is achieved in trend and seasonality.
for `qsmp` the authors do not adivise a default but they use a value close to div(`np`,7).

# Arguments
- `np`: Seasonality.
- `robust`: Robust stl.
- `nl`: Smoothing parameter of the low-pass filter.
- `ns`: Smoothing parameter for the seasonal component.
- `nt`: Smoothing parameter for the trend decomposition.
- `ni`: Number of inner loop cycles.
- `no`: Number of outer loop cycles.
- `spm`: Seasonal post-smoothing.
- `qsmp`: Loess q window for Seasonal post-smoothing.
- `verbose`: If true shows updates for the Seasonal and Trend convergence.
- `cth`: Corvengence threshold for Seasonal and Trend.

# Returns
An `STL` object with the seasonal, trend and remainder components.
* STL: A Seasonal, Trend Decomposition Procedure Based on Loess" Robert B. Cleveland, William S. Cleveland, Jean E. McRae, and Irma Terpenning. Journal of Official Statistics Vol. 6. No. 1, 1990, pp. 3-73 (c) Statistics Sweden.

# Examples
```julia-repl
julia> stl_co2 = stl(co2(),365; robust=true, spm=true)
[ Info: Dataset used in Cleveland et al. paper
[ Info: Corvengence achieved (< 0.01); Stopping computation...
STL{TimeSeries.TimeArray{Union{Missing, Float64},2,
Dates.Date,Array{Union{Missing, Float64},2}}}
(4609×3 TimeSeries.TimeArray{Union{Missing, Float64},2,
Dates.Date,Array{Union{Missing, Float64},2}} 1974-05-17 to 1986-12-31,
"stl(Yn, np=365; nl=365, ns=46091, nt=549, ni=1, no=0, spm=true, qsmp=52)")
julia> plot(stl_co2)
```
"""
function stl(Yv::Vector{T},
             np::Integer;
             robust::Bool=false,
             nl::Integer=nextodd(np),
             ns::Integer=10*length(Yv)+1,
             nt::Integer=nextodd(1.5*np/(1-1.5/ns)),
             ni::Integer=robust ? 1 : 2,
             no::Integer=robust ? 15 : 0,
             spm::Bool=false,
             qsmp::Integer=max(div(np,7),2),
             verbose::Bool=false,
             cth::Float64 = 0.01) where T

    @assert isodd(ns) & (ns>=7) "`ns` is chosen of the basis of knowledge of the time series and on the basis of diagnostic methods; must always be odd and at least 7"

    function B(u)
        (u < 1) & !ismissing(u) ? (1.0-u^2)^2 : 0.0
    end

    N = length(Yv)
    # initial robustness weights
    rhov = ones(N)
    # intial trend
    Tv = Tv0 = zeros(N)
    Sv = Sv0 = zeros(N)
    Rv = Array{Float64,1}(undef,N)
    Cv = Array{Float64,1}(undef,N+2*np)
    scnv = false # seasonal convergence flag
    tcnv = false # trend convergence flag
    o = 0
    while robust | (o <= no)
        for k in 1:ni
            # Updating sesonal and trend components
            ## 1. Detrending (Yv = Tv + Sv)
            Sv = Yv - Tv

            ### Seasonal convergence criterion
            Md = maximum(abs.(skipmissing(Sv-Sv0)))
            M0 = maximum(skipmissing(Sv0))
            m0 = minimum(skipmissing(Sv0))
            scnv = (Md/(M0-m0) < cth)
            if verbose
                println("Outer loop: " * string(o) * " - " * "Inner loop: " * string(k))
                println("Seasonal Convergence: " * string(Md/(M0-m0)))
            end
            Sv0 = Sv

            # Seasonal Smoothing 2,3,4
            ## 2. Cycle-subseries Smoothing
            for csi in 1:np
                Cv[csi:np:N+2*np] = loess(csi:np:N,
                                          Sv[csi:np:N];
                                          q=ns,d=1,rho=rhov[csi:np:N],
                                          predict=(1.0*csi-np):np:(N+np))
            end
            ## 3. Low-Pass Filtering of Smoothed Cycle-Subseries
            ### centered support instead 1:N to balance out machine error
            Lv = loess(1-ceil(N/2):N-ceil(N/2),
                       collect(skipmissing(sma(sma(sma(Cv,np),np),3))),
                       d=1,q=nl,rho=rhov)
            ## 4. Detreending of Smoothed Cycle-Subseries
            ### Lv is substracted to prevent low-frenquency power
            ### from entering the seasonal component.
            Sv = Cv[np+1:end-np] - Lv

            ## 5. Deseasonalizing
            Dv = Yv - Sv

            # Trend Smoothing
            ## 6. Trend Smoothing
            ### centered support instead 1:N to balance out machine error
            ### (floor isntead ceil like in Lv to balance out even lengths)
            Tv = loess(1-floor(N/2):N-floor(N/2),Dv,q=nt,d=1,rho=rhov)

            ### Trend convergence criterion
            Md = maximum(abs.(skipmissing(Tv-Tv0)))
            M0 = maximum(skipmissing(Tv0))
            m0 = minimum(skipmissing(Tv0))
            tcnv = (Md/(M0-m0) < cth)
            if verbose
                println("Trend    Convergence: " * string(Md/(M0-m0)) * "\n")
            end
            Tv0 = Tv

            scnv & tcnv ? break : nothing
        end
        # Computation of robustness weights
        ## These new weights will reduce the influence of transient behaviour
        ## on the trend and seasonal components.

        # Tv and Sv defined everywhere but Rv is not defined
        # where Yv has missing values

        Rv = Yv - Tv - Sv

        if scnv & tcnv
            verbose &&  println("Corvengence achieved (< " * string(cth) * "); Stopping computation...")
            break
        end

        if 0 < o <= no
            smRv = skipmissing(Rv)
            h = 6*median(abs.(smRv))
            rhov = B.(abs.(Rv)/h)
        end
        o += 1
    end

    if spm
        # Using div(np,7) as default to approximate
        # the q=51 for a np=365 chosen in the original paper
        Sv = loess(1-ceil(N/2):N-ceil(N/2),Sv,q=qsmp)
        Rv = Yv - Tv - Sv
    end

    if !scnv
        @warn "Seasonal convergence not achieved (>= " * string(cth) * ")
         Consider a robust estimation"
    end

    if !tcnv
        @warn "Trend convergence not achieved (>= " * string(cth) * ")
         Consider a robust estimation"
    end

    call = "stl(Yn, np="*string(np)*
        "; nl="*string(nl)*
        ", ns="*string(ns)*
        ", nt="*string(nt)*
        ", ni="*string(ni)*
        ", no="*string(no)*
        ", spm="*string(spm)*
        ", qsmp="*string(qsmp)*")"

    return STL{T}(Tv, Sv, Rv)
end


text:
# see license information at https://github.com/JuliaLang/julia-logo-graphics/LICENSE.md

loadcolorscheme(:julia, [
    colorant"#1f83ff", # blue
    colorant"#CA3C32", # red
    colorant"#399746", # green
    colorant"#9259A3", # purple
], "julia", "sequential")

# Colors from first logo
loadcolorscheme(:julialegacy, [
    colorant"#6a82e0", # blue
    colorant"#d8645f", # red
    colorant"#69ad57", # green
    colorant"#ab7ac1", # purple
], "julia", "sequential")

loadcolorscheme(:jblue, [
    colorant"#d2e6ff",
    colorant"#a5cdff",
    colorant"#79b5ff",
    colorant"#4c9cff",
    colorant"#1f83ff",
    colorant"#1969cc",
    colorant"#134f99",
    colorant"#0c3466",
    colorant"#061a33",
], "julia", "sequential")

loadcolorscheme(:jred, [
    colorant"#f4d8d6",
    colorant"#eab1ad",
    colorant"#df8a84",
    colorant"#d5635b",
    colorant"#ca3c32",
    colorant"#a23028",
    colorant"#79241e",
    colorant"#511814",
    colorant"#280c0a",
], "julia", "sequential")

loadcolorscheme(:jgreen, [
    colorant"#d7eada",
    colorant"#b0d5b5",
    colorant"#88c190",
    colorant"#61ac6b",
    colorant"#399746",
    colorant"#2e7938",
    colorant"#225b2a",
    colorant"#173c1c",
    colorant"#0b1e0e",
], "julia", "sequential")

loadcolorscheme(:jredgreen, [
    colorant"#ca3c32",
    colorant"#d5635b",
    colorant"#df8a84",
    colorant"#eab1ad",
    colorant"#f4d8d6",
    colorant"#f9fafb",
    colorant"#d7eada",
    colorant"#b0d5b5",
    colorant"#88c190",
    colorant"#61ac6b",
    colorant"#399746",
], "julia", "diverging")

loadcolorscheme(:jpurple, [
    colorant"#e9deed",
    colorant"#d3bdda",
    colorant"#be9bc8",
    colorant"#a87ab5",
    colorant"#9259a3",
    colorant"#754782",
    colorant"#583562",
    colorant"#3a2441",
    colorant"#1d1221",
], "julia", "sequential")

max_stars_repo_name: cormullion/ColorSchemes

max_issues_repo_name: cormullion/ColorSchemes

max_forks_repo_name: cormullion/ColorSchemes


text:
# Copyright Christopher Johnstone (2021) under the MIT license.
module ChromatinArt

using Random
using Luxor
using Colors

function idx_to_zero_x(width::Int64, idx::Int64)::Int64
    return (idx - 1) % width
end
function idx_to_zero_y(width::Int64, idx::Int64)::Int64
    return trunc(Int64, (idx - 1) / width)
end
function zero_xy_to_idx(width::Int64, x::Int64, y::Int64)::Int64
    return x + y * width + 1
end

# Inspired by https://aip.scitation.org/doi/10.1063/1.2357935
function hamiltonian_mc_step!(path, width::Int64, height::Int64, dir::Int64, aux1::Int64, rng)
    if rand(rng) > 0.5
        # Reverse the list before we start
        reverse!(path)
    end
    # Try to step away from the head
    dir = rand(rng, 1:4)
    x = idx_to_zero_x(width,path[1])
    y = idx_to_zero_y(width,path[1])
    if dir == 1
        x = x + 1
    elseif dir == 2
        x = x - 1
    elseif dir == 3
        y = y + 1
    elseif dir == 4
        y = y - 1
    end

    if x < 0 || x >= width || y < 0 || y >= height
        return
    end
    # Otherwise, make a step
    aux1 = 0
    for i = 1:length(path)
        if path[i] == zero_xy_to_idx(width, x, y)
            aux1 = i
        end
    end
    reverse!(path, 1, aux1 - 1)
end

function generate_path(width::Int64, height::Int64, inputs::Vector{Int64}, outputs::Vector{Int64}, rng)::Vector{Int64}
    result = Vector{Int64}(1:(width*height))
    generate_path!(result, width, height, inputs, outputs, rng)
    return result
end

function generate_path!(path_array, width::Int64, height::Int64, inputs::Vector{Int64}, outputs::Vector{Int64}, rng)
    n_tiles = width * height
    # Init a Hamiltonian path
    path_array[1:end] = 1:n_tiles
    # Reverse every other row to get a Hamiltonian starting path
    for row in 2:2:height
        reverse!(path_array, 1 + (row - 1) * width, row * width)
    end

    aux1::Int64 = 0
    aux2::Int64 = 0

    # Initialize by doing some random iterations
    for i in 1:(n_tiles * 10)
        hamiltonian_mc_step!(path_array, width, height, aux1, aux2, rng)
    end

    while path_array[1] ∉ inputs  || path_array[end] ∉ outputs
        hamiltonian_mc_step!(path_array, width, height, aux1, aux2, rng)
    end
end

struct FractalStructure
    width::Int64
    height::Int64
    n_levels::Int64
    sides::Vector{Vector{Int64}}
end

function Fractal(width::Int64, height::Int64, levels::Int64)::FractalStructure
    return FractalStructure(
        width,
        height,
        levels,
        [
            Vector(1:width),
            Vector(width:width:(width * height)),
            Vector((width * (height - 1) + 1):(width * height)),
            Vector(1:width:(width * (height - 1) + 1))
        ]
    )
end

struct FractalPath
    structure::FractalStructure
    paths::Vector{Vector{Int64}}
end

function idxes_to_dir(idx1::Int64, idx2::Int64, s::FractalStructure)::Int64
    # Returns the direction that connects two indexes.
    # direction = 1 => up
    # direction = 2 => right
    # direction = 3 => down
    # direction = 4 => left
    d = idx2 - idx1
    if d == 1
        return 2
    elseif d == -1
        return 4
    elseif d == s.width
        return 3
    elseif d == -s.width
        return 1
    end

    # If we are here, we have a cross-index move. Figure out which side each entry is on
    if idx1 ∈ s.sides[1] && idx2 ∈ s.sides[3]
        return 1
    elseif idx1 ∈ s.sides[2] && idx2 ∈ s.sides[4]
        return 2
    elseif idx1 ∈ s.sides[3] && idx2 ∈ s.sides[1]
        return 3
    elseif idx1 ∈ s.sides[4] && idx2 ∈ s.sides[2]
        return 4
    end

     # If we are here, then something unexpected happen. This is likely a
     # cross-index move
    error("Unexpected index-index difference!")
end

function idxes_to_available_outputs(idx1::Int64, idx2::Int64, s::FractalStructure)::Vector{Int64}
    d = idxes_to_dir(idx1, idx2, s)
    return s.sides[d]
end

# Transforms a boundary condition in a certain direction
function reflect_bc(direction::Int64, idx::Int64, s::FractalStructure)::Int64
    if direction == 1
        # Going up! Reflect by adding to width * (height - 1)
        return idx + (s.width * (s.height - 1))
    elseif direction == 2
        # Going right. Subtract (width - 1)
        return idx - (s.width - 1)
    elseif direction == 3
        # Going down. Add width and wrap
        return (idx + s.width)  % (s.width * s.height)
    elseif direction == 4
        # Going left. Add (width - 1)
        return idx + (s.width - 1)
    end
end

function generate_fractal_path(structure::FractalStructure)::FractalPath
    # Generate a fractal path of the desired size and with the specified number of levels.
    result::FractalPath = FractalPath(structure, [])
    n_tiles = structure.width * structure.height
    rng = MersenneTwister();

    # Generate the first level directly
    push!(result.paths, zeros(Int64, n_tiles))
    generate_path!(
        result.paths[1],
        structure.width, structure.height,
        [1], [(structure.width * structure.height)],
        rng
    )
    for level in 2:structure.n_levels
        # We need to iterate over level i-1 and generate paths as we go.
        # Create the large path structure
        push!(result.paths, zeros(Int64,(n_tiles)^level))
        generate_path!(
                  @view(result.paths[level][1:n_tiles]),
                  structure.width, structure.height,
                  [1,structure.width],
                  idxes_to_available_outputs(result.paths[level-1][1], result.paths[level-1][2],structure),
                  rng
        )
        # For the next N-2 entries, reflect the tile
        for i in 2:(length(result.paths[level-1]) - 1)
            generate_path!(
                @view(result.paths[level][(((i-1) * n_tiles) + 1):(i * n_tiles)]),
                structure.width, structure.height,
                [reflect_bc(idxes_to_dir(
                        result.paths[level-1][i-1],
                        result.paths[level-1][i],
                        structure),
                    result.paths[level][(i-1) * n_tiles],structure)],
                idxes_to_available_outputs(
                    result.paths[level-1][i],
                    result.paths[level-1][i+1],
                    structure),
                rng
            )
        end
        # Push the last subpath
        generate_path!(
            @view(result.paths[level][(end-n_tiles+1):end]),
            structure.width, structure.height,
            [reflect_bc(idxes_to_dir(
                    result.paths[level-1][end-1],
                    result.paths[level-1][end],
                    structure),
                result.paths[level][end-n_tiles],structure)],
            [(structure.width * structure.height)],
            rng
        )
    end
    return result
end

struct FractalViz
    level_offsets::Vector{Int64}
    jitter_magnitude::Vector{Float64}
    polymerase_variance::Vector{Float64}
    polymerase_base_rate::Float64
    swirl_strength::Float64
    swirl_dropoff::Float64
    n_chromosomes::Int64
    colors
    skip_loops::Bool
end

function swirl_transform(strength::Float64, dropoff::Float64, total_width::Int64, total_height::Int64, x::Float64, y::Float64)::Tuple{Float64,Float64}
    # Recenter x/y
    recentered_x = x - (total_width / 2.0)
    recentered_y = y - (total_height / 2.0)
    max_r = sqrt((total_width/2)^2+(total_height/2)^2)
    # Convert to polar coordinates
    r, theta = sqrt(recentered_x^2 + recentered_y^2), atan(recentered_y, recentered_x)
    # Transform with a swirl pattern
    t_theta = theta + strength * sin(4.0*π*(r/(max_r))) * exp(-dropoff * r/max_r)
    return (r * sin(t_theta) + (total_width / 2.0), r * cos(t_theta) + (total_height / 2.0))
end

function get_transformed_loc(idx::Int64, path::FractalPath, viz::FractalViz, plot_width, plot_height, jitter)
    indicies = [path.paths[level][
    ceil(Int64,idx / (
        path.structure.width *
        path.structure.height)^(
            path.structure.n_levels - level
        ))] for level in 1:path.structure.n_levels]
    x_offsets = map(idx->idx_to_zero_x(path.structure.width,idx),indicies) .* viz.level_offsets + map(level->jitter[level][indicies[level],1], 1:path.structure.n_levels)
    y_offsets = map(idx->idx_to_zero_y(path.structure.width,idx),indicies) .* viz.level_offsets + map(level->jitter[level][indicies[level],2], 1:path.structure.n_levels)
    x_loc = sum(x_offsets)
    y_loc = sum(y_offsets)
    return swirl_transform(viz.swirl_strength, viz.swirl_dropoff, plot_width, plot_height, x_loc, y_loc)
end

function draw_fractal_path(path::FractalPath, viz::FractalViz, filename)
    # Generate jitter offsets
    jitter = [clamp.(randn(
        (path.structure.width * path.structure.height)^level,2) * viz.jitter_magnitude[level],
        -viz.jitter_magnitude[level], viz.jitter_magnitude[level])
        for level in 1:path.structure.n_levels]
    # Generate polymerase propensities
    rnap_propensity = [
        clamp.(1.0 .+ (randn((path.structure.width * path.structure.height)^level) * viz.polymerase_variance[level]),0.0, 5.0)
        for level in 1:path.structure.n_levels
    ]
    plot_width = sum([viz.level_offsets[level] * path.structure.width for level in 1:path.structure.n_levels])
    plot_height = sum([viz.level_offsets[level] * path.structure.height for level in 1:path.structure.n_levels])

    n_per_chromosome = ceil(Int64, ((path.structure.width * path.structure.height)^path.structure.n_levels) / viz.n_chromosomes)
    Drawing(plot_width + 5, plot_height + 5, filename)
    setline(0.5)
    newpath()
    move(0,0)
    setcolor(viz.colors[1])
    polymerase_accum = []
    for i in 1:length(path.paths[end])
        if i % n_per_chromosome == 1
            # New chromosome!
            do_action(:stroke)
            setcolor(viz.colors[ceil(Int64, i / n_per_chromosome)])
            newpath()
        end
        transformed = get_transformed_loc(i, path, viz, plot_width, plot_height, jitter)
        line(Point(transformed...))
        indicies = [path.paths[level][
        ceil(Int64,i / (
            path.structure.width *
            path.structure.height)^(
                path.structure.n_levels - level
        ))] for level in 1:path.structure.n_levels]

        polymerase_rate = viz.polymerase_base_rate * prod(map(level->rnap_propensity[level][indicies[level]], 1:path.structure.n_levels))
        if rand() < polymerase_rate
            push!(polymerase_accum, transformed)
        end
    end
    do_action(:stroke)

    # Add loops at level n-1
    if ~viz.skip_loops
        gsave()
        setcolor("#9e52e5")
        setline(0.8)
        potential_loop_idx = -1
        for i in 2:(length(path.paths[end-1])-1)
            # Reset if we are in a new subpath
            if i % (path.structure.width * path.structure.height) == 0
                potential_loop_idx = -1
                continue
            end
            # Locate things like:
            # ->
            # <-
            # within each segment
            segment_diff = path.paths[end-1][i+1] - path.paths[end-1][i]
            if potential_loop_idx == -1 && segment_diff == 1
                potential_loop_idx = i
            end

            if potential_loop_idx != -1 && segment_diff == -1 && abs(path.paths[end-1][i+1] - path.paths[end-1][potential_loop_idx]) == path.structure.width
                p1 = Point(get_transformed_loc(
                    ((potential_loop_idx) * path.structure.width * path.structure.height),
                    path, viz, plot_width, plot_height, jitter))
                p2 = Point(get_transformed_loc(
                    ((potential_loop_idx) * path.structure.width * path.structure.height) + 1,
                    path, viz, plot_width, plot_height, jitter))
                p3 = Point(get_transformed_loc(
                    (i * path.structure.width * path.structure.height),
                    path, viz, plot_width, plot_height, jitter))
                p4 = Point(get_transformed_loc(
                    (i * path.structure.width * path.structure.height) + 1,
                    path, viz, plot_width, plot_height, jitter))
                # Skip if the distance between points is out of the desired range
                if sqrt((p1[1] - p4[1])^2 + (p1[2] - p4[2])^2) > 3
                    potential_loop_idx = -1
                    continue
                end
                # We got a loop! Generate it!
                ellipseinquad([p1, p2, p3, p4], :stroke)
                potential_loop_idx = -1
            end
        end
        grestore()
    end

    # Place polymerases
    gsave()
    print("# polymerases:", length(polymerase_accum), "\n")
    for polymerase_loc in polymerase_accum
        setcolor(0.4,0.4,0.4,0.6)
        circle(polymerase_loc...,0.7,:fill)
        setcolor(0.4,0.4,0.4,0.9)
        setline(0.2)
        circle(polymerase_loc...,0.7,:stroke)
        setcolor(0.6, 0.6, 0.6, 1.0)
        setline(0.5)
        dx = randn() * 3
        dy = randn() * 3
        drawbezierpath(bezierfrompoints(
            Point(polymerase_loc...),
            Point(polymerase_loc[1] + (0.25 * dx) + randn() * 0.3, polymerase_loc[2] + (0.25 * dy)),
            Point(polymerase_loc[1] + (0.75 * dx) + randn() * 0.3, polymerase_loc[2] + (0.75 * dy)),
            Point(polymerase_loc[1] + dx, polymerase_loc[2] + dy)), :stroke, close=false)
    end
    grestore()

    finish()
end

# Precompile hints
generate_path(4,4,[1],[2], MersenneTwister());
end # module


text:
module TestStatic

using Test
import ThreadPools: StaticPool
using ThreadPools

include("util.jl")

@testset "StaticPool" begin

    @testset "tforeach" begin
        @testset "foreground" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            primary = Threads.nthreads() == 1
            fn! = (x) -> begin
                Threads.threadid() == 1 && (primary = true)
                x.data += 1
            end
            pool = StaticPool()
            tforeach(fn!, pool, objs)
            close(pool)
            @test [x.data for x in objs] == collect(2:N+1)
            @test primary
            @inferred tforeach(fn!, pool, objs)
        end

        @testset "background" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            fn! = (x) -> begin
                Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
                x.data += 1
            end
            pool = StaticPool(2)
            tforeach(fn!, pool, objs)
            close(pool)
            @test [x.data for x in objs] == collect(2:N+1)
            @inferred tforeach(fn!, pool, objs)
        end
    end

    @testset "tmap" begin
        @testset "foreground" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            primary = Threads.nthreads() == 1
            fn! = (x) -> begin
                Threads.threadid() == 1 && (primary = true)
                x.data
            end
            pool = StaticPool()
            @test tmap(fn!, pool, objs) == collect(1:N)
            close(pool)
            @test primary
            @inferred tmap(fn!, pool, objs)
        end

        @testset "background" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            fn! = (x) -> begin
                Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
                x.data
            end
            pool = StaticPool(2)
            @test tmap(fn!, pool, objs) == collect(1:N)
            close(pool)
            @inferred tmap(fn!, pool, objs)
        end

        if Threads.nthreads() > 2
            t = @tspawnat 2 begin
                @testset "schedule from background" begin
                    N = 2 * Threads.nthreads()
                    objs = [TestObj(x) for x in 1:N]
                    fn! = (x) -> begin
                        Threads.threadid() == 1 && error("Task on primary")
                        Threads.threadid() == 2 && error("Task on scheduler")
                        (x.data, Threads.threadid())
                    end
                    pool = StaticPool(3)
                    result = tmap(fn!, pool, objs)
                    @test [x[1] for x in result] == collect(1:N)
                    @test sort(collect(Set([x[2] for x in result]))) == collect(3:Threads.nthreads())
                    close(pool)
                    @inferred tmap(fn!, pool, objs)
                end
            end
            wait(t)
        end
    end

    @testset "twith" begin
        @testset "foreground" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            primary = Threads.nthreads() == 1
            twith(StaticPool()) do pool
                tforeach(pool, objs) do x
                    Threads.threadid() == 1 && (primary = true)
                    x.data += 1
                end
            end
            @test [x.data for x in objs] == collect(2:N+1)
            @test primary
        end

        @testset "background" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            twith(StaticPool(2)) do pool
                tforeach(pool, objs) do x
                    Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
                    x.data += 1
                end
            end
            @test [x.data for x in objs] == collect(2:N+1)
        end
    end

    @testset "@tthreads" begin
        @testset "foreground" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            primary = Threads.nthreads() == 1
            pool = StaticPool()
            @tthreads pool for obj in objs
                Threads.threadid() == 1 && (primary = true)
                obj.data += 1
            end
            close(pool)
            @test primary
            @test [x.data for x in objs] == collect(2:N+1)
        end

        @testset "background" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            pool = StaticPool(2)
            @tthreads pool for obj in objs
                Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
                obj.data += 1
            end
            close(pool)
            @test [x.data for x in objs] == collect(2:N+1)
        end
    end

end

end # module


text:
using HyperHessians: HyperHessians, hessian, hessian!, Chunk, HessianConfig, HessianConfigThreaded, HyperDual
using DiffTests
using ForwardDiff
using Test
using SIMD

if Threads.nthreads() == 1
    @warn "Not properly testing threaded behavior"
end

@testset "scalar" begin
f(x) = exp(x) / sqrt(sin(x)^3 + cos(x)^3);
x = rand()
@test hessian(f, x) ≈  ForwardDiff.derivative(x->ForwardDiff.derivative(f, x), x)
end # testset

@testset "correctness" begin
for f in (DiffTests.rosenbrock_1,
          DiffTests.rosenbrock_2,
          DiffTests.rosenbrock_3,
          #DiffTests.rosenbrock_4,
          DiffTests.ackley,
          DiffTests.self_weighted_logit,
          DiffTests.nested_array_mul,)
    for n in (1, 4, HyperHessians.DEFAULT_CHUNK_THRESHOLD, HyperHessians.DEFAULT_CHUNK_THRESHOLD + 7)
        if n == 1 && (f == DiffTests.rosenbrock_4 || f == DiffTests.nested_array_mul)
            continue
        end
        x = rand(n)
        H  = ForwardDiff.hessian(f, x)

        for chunk in (1, n÷2, n)
            if chunk <= 0
                continue
            end
            @info "f=$f, n=$n, chunk=$chunk"

            cfg = HessianConfig(x, Chunk{chunk}())
            @test H ≈ hessian(f, x, cfg)
            @test H ≈ hessian(f, x, cfg)

            cfg_threaded = HessianConfigThreaded(x, Chunk{chunk}())
            @test H ≈ hessian(f, x, cfg_threaded)
            @test H ≈ hessian(f, x, cfg_threaded)
        end
    end
end
end # testset

function ackley_stable(x::AbstractVector{T}) where {T}
    a, b, c = T(20.0), T(-0.2), T(2.0*π)
    len_recip = T(inv(length(x)))
    sum_sqrs = sum(i->i^2, x; init=zero(eltype(x)))
    sum_cos = sum(i->cos(c*i), x; init=zero(eltype(x)))
    return (-a * exp(b * sqrt(len_recip*sum_sqrs)) -
            exp(len_recip*sum_cos) + a + exp(T(1)))
end

@testset "Float32" begin
x = rand(Float32, 8)
@test hessian(ackley_stable, x) isa Matrix{Float32}
@test hessian(ackley_stable, x)  ≈ ForwardDiff.hessian(ackley_stable, x)
@test hessian(DiffTests.ackley, x) ≈ ForwardDiff.hessian(DiffTests.ackley, x)
end

using StaticArrays
@testset "StaticArrays" begin
x = rand(SVector{4})
f = DiffTests.rosenbrock_3
@test ForwardDiff.hessian(f, x) ≈ hessian(f, x)
@test hessian(f, x) isa MMatrix
@test_broken hessian(f, x) isa SMatrix
end # testset

@testset "No spurious promotions primitives" begin
h = HyperDual(0.8f0, Vec(0.7f0, 0.7f0), Vec(0.7f0, 0.7f0))
for (fsym, _, _) in HyperHessians.DIFF_RULES
    if fsym in (:asec, :acsc, :asecd)
        hv = inv(h)
    else
        hv = h
    end
    f = @eval $fsym
    try
        v = f(hv)
        @test v isa HyperDual{2, Float32}
    catch e
        e isa DomainError || rethrow()
    end
end
end #testset


text:
using DimensionalData, Test, Dates
using DimensionalData.LookupArrays, DimensionalData.Dimensions
using .LookupArrays: shiftlocus, maybeshiftlocus

@testset "reverse" begin
    @testset "dimension" begin
        revdima = reverse(X(Sampled(10:10:20; order=ForwardOrdered(), span=Regular(10))))
        @test index(revdima) == 20:-10:10
        @test order(revdima) === ReverseOrdered()
        @test span(revdima) === Regular(-10)
    end

    A = [1 2 3; 4 5 6]
    da = DimArray(A, (X(10:10:20), Y(300:-100:100)); name=:test)
    s = DimStack(da)

    rev = reverse(da; dims=Y)
    @test rev == [3 2 1; 6 5 4]
    @test index(rev, X) == 10:10:20
    @test index(rev, Y) == 100:100:300
    @test span(rev, Y) == Regular(100)
    @test order(rev, Y) == ForwardOrdered()
    @test order(rev, X) == ForwardOrdered()

    @testset "NoLookup dim index is not reversed" begin
        da = DimArray(A, (X(), Y()))
        revd = reverse(da)
        @test index(revd) == axes(da)
    end

    @testset "stack" begin
        dimz = (X([:a, :b]), Y(10.0:10.0:30.0))
        da1 = DimArray(A, dimz; name=:one)
        da2 = DimArray(Float32.(2A), dimz; name=:two)
        da3 = DimArray(Int.(3A), dimz; name=:three)

        s = DimStack((da1, da2, da3))
        rev_s = reverse(s; dims=X)
        @test rev_s[:one] == [4 5 6; 1 2 3]
        @test rev_s[:three] == [12 15 18; 3 6 9]
    end
end

@testset "reorder" begin
    A = [1 2 3; 4 5 6]
    da = DimArray(A, (X(10:10:20), Y(300:-100:100)); name=:test)
    s = DimStack(da)

    reo = reorder(da, ReverseOrdered())
    @test reo == [4 5 6; 1 2 3]
    @test index(reo, X) == 20:-10:10
    @test index(reo, Y) == 300:-100:100
    @test order(reo, X) == ReverseOrdered()
    @test order(reo, Y) == ReverseOrdered()

    reo = reorder(da, X=>ForwardOrdered(), Y=>ReverseOrdered())
    @test reo == A
    @test index(reo, X) == 10:10:20
    @test index(reo, Y) == 300:-100:100
    @test order(reo, X) == ForwardOrdered()
    @test order(reo, Y) == ReverseOrdered()

    reo = reorder(da, X=>ReverseOrdered(), Y=>ForwardOrdered())
    @test reo == [6 5 4; 3 2 1]
    @test index(reo, X) == 20:-10:10
    @test index(reo, Y) == 100:100:300
    @test order(reo, X) == ReverseOrdered()
    @test order(reo, Y) == ForwardOrdered()

    revallis = reverse(da; dims=(X, Y))
    @test index(revallis) == (20:-10:10, 100:100:300)
    @test order(revallis) == (ReverseOrdered(), ForwardOrdered())

    d = reorder(dims(da, Y), ForwardOrdered())
    @test order(d) isa ForwardOrdered
    @test index(d) == 100:100:300
end

@testset "modify" begin
    A = [1 2 3; 4 5 6]
    dimz = (X(10:10:20), Y(300:-100:100))
    @testset "array" begin
        da = DimArray(A, dimz)
        mda = modify(A -> A .> 3, da)
        @test dims(mda) === dims(da)
        @test mda == [false false false; true true true]
        @test typeof(parent(mda)) == BitArray{2}
        @test_throws ErrorException modify(A -> A[1, :], da)
    end
    @testset "dataset" begin
        da1 = DimArray(A, dimz; name=:da1)
        da2 = DimArray(2A, dimz; name=:da2)
        s = DimStack(da1, da2)
        ms = modify(A -> A .> 3, s)
        @test parent(ms) == (da1=[false false false; true true true],
                              da2=[false true  true ; true true true])
        @test typeof(parent(ms[:da2])) == BitArray{2}
    end
    @testset "dimension" begin
        dim = X(Sampled(10:10:20))
        mdim = modify(x -> 3 .* x, dim)
        @test index(mdim) == 30:30:60 # in Julia 1.6: typeof(30:30:60)==StepRange ; in Julia 1.7 typeof(30:30:60)==StepRangeLen
        da = DimArray(A, dimz)
        mda = modify(y -> vec(4 .* y), da, Y)
        @test index(mda, Y) == [1200.0, 800.0, 400.0]
    end
end

@testset "dimwise" begin
    A2 = [1 2 3; 4 5 6]
    B1 = [1, 2, 3]
    da2 = DimArray(A2, (X([20, 30]), Y([:a, :b, :c])))
    db1 = DimArray(B1, (Y([:a, :b, :c]),))
    dc1 = dimwise(+, db1, db1)
    @test dc1 == [2, 4, 6]
    dc2 = dimwise(+, da2, db1)
    @test dc2 == [2 4 6; 5 7 9]
    dc4 = dimwise(+, da2, db1)

    A3 = cat([1 2 3; 4 5 6], [11 12 13; 14 15 16]; dims=3)
    da3 = DimArray(A3, (X, Y, Z))
    db2 = DimArray(A2, (X, Y))
    dc3 = dimwise(+, da3, db2)
    @test dc3 == cat([2 4 6; 8 10 12], [12 14 16; 18 20 22]; dims=3)
    dc3 = dimwise!(+, da3, da3, db2)

    A3 = cat([1 2 3; 4 5 6], [11 12 13; 14 15 16]; dims=3)
    da3 = DimArray(A3, (X([20, 30]), Y([:a, :b, :c]), Z(10:10:20)))
    db1 = DimArray(B1, (Y([:a, :b, :c]),))
    dc3 = dimwise(+, da3, db1)
    @test dc3 == cat([2 4 6; 5 7 9], [12 14 16; 15 17 19]; dims=3)

    @testset "works with permuted dims" begin
        db2p = permutedims(da2)
        dc3p = dimwise(+, da3, db2p)
        @test dc3p == cat([2 4 6; 8 10 12], [12 14 16; 18 20 22]; dims=3)
    end

end

@testset "shiftlocus" begin
    dim = X(Sampled(1.0:3.0, ForwardOrdered(), Regular(1.0), Intervals(Center()), NoMetadata()))
    @test index(shiftlocus(Start(), dim)) === 0.5:1.0:2.5
    @test index(shiftlocus(End(), dim)) === 1.5:1.0:3.5
    @test index(shiftlocus(Center(), dim)) === 1.0:1.0:3.0
    @test locus(shiftlocus(Start(), dim)) === Start()
    @test locus(shiftlocus(End(), dim)) === End()
    @test locus(shiftlocus(Center(), dim)) === Center()
    dim = X(Sampled([3, 4, 5], ForwardOrdered(), Regular(1), Intervals(Start()), NoMetadata()))
    @test val(shiftlocus(End(), dim)) == [4, 5, 6]
    @test val(shiftlocus(Center(), dim)) == [3.5, 4.5, 5.5]
    @test val(shiftlocus(Start(), dim)) == [3, 4, 5]
    dim = X(Sampled([3, 4, 5], ForwardOrdered(), Regular(1), Intervals(End()), NoMetadata()))
    @test val(shiftlocus(End(), dim)) == [3, 4, 5]
    @test val(shiftlocus(Center(), dim)) == [2.5, 3.5, 4.5]
    @test val(shiftlocus(Start(), dim)) == [2, 3, 4]

    dates = DateTime(2000):Month(1):DateTime(2000, 12)
    ti = Ti(Sampled(dates, ForwardOrdered(), Regular(Month(1)), Intervals(Start()), NoMetadata()))
    @test index(shiftlocus(Center(), ti)) == dates .+ (dates .+ Month(1) .- dates) ./ 2

    bnds = vcat((0.5:2.5)', (1.5:3.5)')
    dim = X(Sampled(1.0:3.0, ForwardOrdered(), Explicit(bnds), Intervals(Center()), NoMetadata()))
    start_dim = shiftlocus(Start(), dim)
    @test index(start_dim) == [0.5, 1.5, 2.5]
    @test locus(start_dim) == Start()
    end_dim = shiftlocus(End(), start_dim)
    @test index(end_dim) == [1.5, 2.5, 3.5]
    @test locus(end_dim) == End()
    center_dim = shiftlocus(Center(), end_dim)
    @test index(center_dim) == index(dim)
    @test locus(center_dim) == Center()
end

@testset "maybeshiftlocus" begin
    dim = X(Sampled(1.0:3.0, ForwardOrdered(), Regular(1.0), Intervals(Center()), NoMetadata()))
    @test val(maybeshiftlocus(Start(), dim)) == 0.5:1.0:2.5
    dim = X(Sampled(1.0:3.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))
    @test val(maybeshiftlocus(Start(), dim)) == 1.0:3.0
end

@testset "dim2boundsmatrix" begin
    @testset "Regular span" begin
        dim = X(Sampled(1.0:3.0, ForwardOrdered(), Regular(1.0), Intervals(Center()), NoMetadata()))
        @test Dimensions.dim2boundsmatrix(dim) == [0.5 1.5 2.5
                                                        1.5 2.5 3.5]
        dim = X(Sampled(1.0:3.0, ForwardOrdered(), Regular(1.0), Intervals(Start()), NoMetadata()))
        @test Dimensions.dim2boundsmatrix(dim) == [1.0 2.0 3.0
                                                        2.0 3.0 4.0]
        dim = X(Sampled(1.0:3.0, ReverseOrdered(), Regular(1.0), Intervals(End()), NoMetadata()))
        @test Dimensions.dim2boundsmatrix(dim) == [0.0 1.0 2.0
                                                        1.0 2.0 3.0]
        dim = X(Sampled(3.0:-1:1.0, ReverseOrdered(), Regular(1.0), Intervals(Center()), NoMetadata()))
        @test Dimensions.dim2boundsmatrix(dim) == [2.5 1.5 0.5
                                                        3.5 2.5 1.5]
    end
    @testset "Explicit span" begin
        dim = X(Sampled(1.0:3.0, ForwardOrdered(),
                Explicit([0.0 1.0 2.0; 1.0 2.0 3.0]), Intervals(End()), NoMetadata()))
        @test Dimensions.dim2boundsmatrix(dim) == [0.0 1.0 2.0
                                                        1.0 2.0 3.0]
    end
end


text:
using Knet: KnetArray, RNN, Param, relu, sumabs2, mat, conv4, pool, dropout,
    nll, gpu, xavier
using Statistics: mean

struct LinearLayer; w; b; f; pdrop; end

LinearLayer(i::Int, o::Int, f=relu; pdrop=0.5,
            atype=gpu()>=0 ? KnetArray{Float64} : Array{Float64}) = LinearLayer(
                Param(atype(xavier(o, i))), Param(atype(zeros(o))), f, pdrop)

(l::LinearLayer)(x) = l.f.(l.w * mat(dropout(x, l.pdrop)) .+ l.b)
(l::LinearLayer)(x, y) = sumabs2(y - l(x)) / size(y,2)

struct LinearChain
    layers
    LinearChain(layers...) = new(layers)
end

(c::LinearChain)(x) = (for l in c.layers; x = l(x); end; x)
(c::LinearChain)(x, y) = sumabs2(y - c(x)) / size(y,2)

struct CategoricalChain
    layers
    CategoricalChain(layers...) = new(layers)
end

(c::CategoricalChain)(x) = (for l in c.layers; x = l(x); end; x)
(c::CategoricalChain)(x, y) = nll(c(x), y)

struct Conv; w; b; f; pdrop; end

Conv(w1::Int, w2::Int, cx::Int, cy::Int, f=relu; pdrop=0,
     atype=gpu()>=0 ? KnetArray{Float64} : Array{Float64}) =
         Conv(Param(atype(xavier(w1, w2, cx, cy))),
              Param(atype(zeros(1, 1, cy, 1))),
              f,
              pdrop)

(c::Conv)(x) = c.f.(pool(conv4(c.w, dropout(x, c.pdrop)) .+ c.b))

struct OneLayerBiRNN; input; rnn; output; b; pdrop; end

OneLayerBiRNN(input::Int, embed::Int, hidden::Int, output::Int; pdrop=0,
              scale=0.01, atype=gpu()>=0 ? KnetArray{Float64} : Array{Float64},
              rnnType=:gru, bidirectional=false) =
    OneLayerBiRNN(Param(atype(randn(embed, input) * scale)),
                  RNN(embed, hidden, rnnType=rnnType, dataType=Float64,
                      bidirectional=bidirectional),
                  Param(atype(randn(output, hidden) * scale)),
                  Param(atype(zeros(output))),
                  pdrop)

function (c::OneLayerBiRNN)(input)
    embed = c.input[:, permutedims(hcat(input...))]
    embed = dropout(embed, c.pdrop)
    hiddenoutput = c.rnn(embed)
    hiddenoutput = dropout(hiddenoutput, c.pdrop)

    return c.output * hiddenoutput[:,:,end] .+ c.b
end

(c::OneLayerBiRNN)(input,output) = nll(c(input), output)

struct TwoTextsClassifier; input; rnn; output; b; pdrop; end

TwoTextsClassifier(input::Int, embed::Int, hidden::Int, output::Int; pdrop=0,
                   scale=0.01, atype=gpu()>=0 ? KnetArray{Float64} :
                   Array{Float64}, rnnType=:gru, bidirectional=true) =
    TwoTextsClassifier(Param(atype(randn(embed, input) * scale)),
                       RNN(embed, hidden, rnnType=rnnType, dataType=Float64,
                           bidirectional=bidirectional),
                       Param(atype(randn(output, 4hidden) * scale)),
                       Param(atype(zeros(output))),
                       pdrop)

function (c::TwoTextsClassifier)(input)
    input1 = hcat(input...)[1, :]
    input2 = hcat(input...)[2, :]

    embed1 = c.input[:, permutedims(hcat(input1...))]
    embed1 = dropout(embed1, c.pdrop)

    embed2 = c.input[:, permutedims(hcat(input1...))]
    embed2 = dropout(embed2, c.pdrop)

    hiddenoutput1 = c.rnn(embed1)
    hiddenoutput2 = c.rnn(embed2)

    H, B, W = size(hiddenoutput1)

    hiddenoutput1 = reshape(hiddenoutput1, H, :)
    hiddenoutput2 = reshape(hiddenoutput2, H, :)

    hiddenoutput = vcat(hiddenoutput1, hiddenoutput2)
    hiddenoutput = reshape(hiddenoutput, 2H, B, W)
    hiddenoutput = sum(hiddenoutput, dims=3)[:, :, end]
    hiddenoutput = dropout(hiddenoutput, c.pdrop)

    return c.output * hiddenoutput .+ c.b
end

(c::TwoTextsClassifier)(input,output) = nll(c(input), output)

struct TwoLayerBiRNN; input; rnn; rnn2; output; b; pdrop; end

TwoLayerBiRNN(input::Int, embed::Int, hidden1::Int, hidden2::Int,
              output::Int; pdrop=0, scale=0.01,
              atype=gpu()>=0 ? KnetArray{Float64} : Array{Float64},
              rnnType=:gru, bidirectional=false) = TwoLayerBiRNN(
                  Param(atype(randn(embed, input) * scale)),
                  RNN(embed, hidden1, rnnType=rnnType, dataType=Float64,
                      bidirectional=bidirectional),
                  RNN(2hidden1, hidden2, rnnType=rnnType, dataType=Float64,
                      bidirectional=bidirectional),
                  Param(atype(randn(output, 2hidden2) * scale)),
                  Param(atype(zeros(output))),
                  pdrop)

function (c::TwoLayerBiRNN)(input)
    embed = c.input[:, permutedims(hcat(input...))]
    embed = dropout(embed, c.pdrop)
    hiddenoutput = c.rnn(embed)
    hiddenoutput = dropout(hiddenoutput, c.pdrop)
    hiddenoutput2 = c.rnn2(hiddenoutput)
    hiddenoutput2 = dropout(hiddenoutput2, c.pdrop)

    return c.output * hiddenoutput2[:,:,end] .+ c.b
end

(c::TwoLayerBiRNN)(input,output) = nll(c(input), output)

predict(model, x) = map(i->i[1], findmax(Array(model(x)),dims=1)[2])
accuracy(model, x, y) = mean(y .== predict(model, x))
accuracy(model, data) = mean(accuracy(model,x,y) for (x,y) in data)

max_stars_repo_name: alperkesen/AutoML

max_issues_repo_name: alperkesen/AutoML


text:
"""
    op = JopSft(dom,freqs,dt)

Polychromatic slow Fourier transforms for a specified list of frequencies. Tranform
along the fast dimension of `dom::JetSpace{<:Real}`.

Notes:
    - it is expected that the domain is 2D real of size [nt,ntrace]
    - the range will be 2D complex of size [length(freqs),ntrace]
    - regarding the factor of (2/n): the factor "1/n" is from the Fourier transform, and the "2" is from Hermittian symmetry.
"""
function JopSft(dom::JetSpace{T}, freqs::Vector, dt) where {T<:Real}
    @assert length(size(dom)) == 2
    nf = length(freqs)
    n2 = size(dom,2)
    JopLn(dom = dom, rng = JetSpace(Complex{T}, nf, n2), df! = JopSft_df!, df′! = JopSft_df′!, s = (freqs=freqs, dt=dt))
end
export JopSft

function JopSft_df!(rngvec::AbstractArray{Complex{T},N}, domvec::AbstractArray{T,N}; freqs, dt, kwargs...) where {T,N}
    nf = length(freqs)
    n1 = size(domvec,1)
    _domvec = reshape(domvec, n1, :)
    _rngvec = reshape(rngvec, nf, :)
    @assert size(_domvec,2) == size(rngvec,2)
    n2 = size(_domvec,2)

    _rngvec .= 0
    for kfreq = 1:nf
        for k2 = 1:n2
            for k1 = 1:n1
                t = dt * (k1 - 1)
                phs = exp(-im * 2 * pi * freqs[kfreq] * t)
                _rngvec[kfreq,k2] += _domvec[k1,k2] * phs
            end
        end
    end

    _rngvec .*= (2/n1)
end

function JopSft_df′!(domvec::AbstractArray{T,N}, rngvec::AbstractArray{Complex{T},N}; freqs, dt, kwargs...) where {T,N}
    nf = length(freqs)
    n1 = size(domvec,1)
    _domvec = reshape(domvec, n1, :)
    _rngvec = reshape(rngvec, nf, :)
    @assert size(_domvec,2) == size(rngvec,2)
    n2 = size(_domvec,2)

    _domvec .= 0
    for kfreq = 1:nf
        for k2 = 1:n2
            for k1 = 1:n1
                t = dt * (k1 - 1)
                phs = exp(-im * 2 * pi * freqs[kfreq] * t)
                _domvec[k1,k2] += real(_rngvec[kfreq,k2] * conj(phs))
            end
        end
    end

    _domvec .*= (2/n1)
end


text:
@testset "points.jl" begin
    p1 = rad2deg(Point(0.5, 1.0))
    p2 = rad2deg(Point(1.0, 0.5))
    p3 = rad2deg(Point(0.6, 1.0))  # above p1

    rs = RouteSection(p1,p2)

    @test midpoint(p1, p2).λ ≈ rad2deg(0.8106650896663075) atol = 0.0001
    @test midpoint(p1, p2).ϕ ≈ rad2deg(0.76484691145) atol = 0.0001
    @test midpoint(rs).λ ≈ rad2deg(0.8106650896663075) atol = 0.0001
    @test midpoint(rs).ϕ ≈ rad2deg(0.76484691145) atol = 0.0001

    @test intermediate_point(p1, p2, 0.5).λ ≈ rad2deg(0.8106650896663075) atol = 0.0001
    @test intermediate_point(p1, p2, 0.5).ϕ ≈ rad2deg(0.76484691145) atol = 0.0001
    @test intermediate_point(rs, 0.5).λ ≈ rad2deg(0.8106650896663075) atol = 0.0001
    @test intermediate_point(rs, 0.5).ϕ ≈ rad2deg(0.76484691145) atol = 0.0001

    dst = 3888000.0
    brg = 5.81412807071
    @test destination_point(p1, dst, rad2deg(brg)).ϕ ≈ rad2deg(1.0) atol = 0.01
    @test destination_point(p1, dst, rad2deg(brg)).λ ≈ rad2deg(0.5) atol = 0.01

    brg13 = 0.5
    brg23 = 1.5
    brg13_ambiguous = 2.61799387799
    @test isnan(intersection_point(p2, p1, rad2deg(brg13_ambiguous), rad2deg(brg23)).ϕ) == true

    @test intersection_point(p1, p2, rad2deg(brg13), rad2deg(brg23)).ϕ ≈ rad2deg(0.86121333) atol = 0.001
    @test intersection_point(p1, p2, rad2deg(brg13), rad2deg(brg23)).λ ≈ rad2deg(1.31431535) atol = 0.001

    p1ni = Point(50.0, 10.0)
    p2ni = Point(50.0, 30.0)
    p3ni = Point(40.0, 20.0)
    p4ni = Point(60.0, 20.0)

    @test intersection_point(p1ni, p2ni, p3ni, p4ni).ϕ ≈ 50.4313888888 atol = 0.0001
    @test intersection_point(p1ni, p2ni, p3ni, p4ni).λ ≈ 20.0 atol = 0.0001
    @test intersection_point(p1ni, p2ni, p3ni, p2ni).λ ≈ 30.0 atol = 0.0001

    p4ni2 = Point(49.0, 20.0)
    @test isnan(intersection_point(p1ni, p2ni, p3ni, p4ni2).λ)

    p1a = Point(10.0374, 0.0)
    p2a = Point(-90.0, 0.0)
    p3a = Point(0.0, 0.0)
    @test isinf(intersection_point(p3a, p2a, p1a, p3a).ϕ)

    p7 = Point(5.0, 0.0)
    p8 = Point(-5.0, 0.0)
    brg7 = 0.0
    brg8 = 0.0
    @test isinf(intersection_point(p7, p8, brg7, brg8).λ)

    #TODO Check if intersection_point also works when crossing the datum line

    #TODO Uncomment after network.jl has been added
    # ap = Airspace("test", [Point(50.0, 10.0), Point(50.0, 30.0),
    # Point(60.0, 30.0), Point(60.0, 10.0), Point(50.0, 10.0)])
    #
    # bx = box_spherical_polygon(ap.polygon)
    # @test bx[3].ϕ ≈ 60.378348124804496 atol = 0.00001
    # # @test isnan(intersection_point(Point(40.0, 20.0), Point(50.0, 20.0), ap).ϕ)

    @test opposite_point(p1).ϕ ≈ -28.6478897565 atol = 0.0001
    @test opposite_point(p1).λ ≈ -122.704220487 atol = 0.0001

    @test max_latitude(0.0, 45.0) ≈ 45.0 atol = 0.0001
    @test max_latitude(rad2deg(1.0), 90.0) ≈ rad2deg(1.0) atol = 0.0001
    @test max_latitude(rad2deg(1.0), 0.0) ≈ 90.0 atol = 0.0001
    @test max_latitude(rad2deg(-1.0), 0.0) ≈ 90.0 atol = 0.0001
    @test max_latitude(rad2deg(-1.0), 180.0) ≈ 90.0 atol = 0.0001

    p5 = Point(45.0, 45.0)
    @test closest_point_to_pole(p5, 45.0).ϕ ≈ 60.0 atol = 0.001
    @test closest_point_to_pole(p5, 45.0).λ ≈ 99.73561 atol = 0.001
    @test closest_point_to_pole(p5, 89.0).ϕ ≈ 45.00873 atol = 0.001
    @test closest_point_to_pole(p5, 89.0).λ ≈ 46.41407 atol = 0.001
    @test closest_point_to_pole(p5, 135.0).ϕ ≈ -60.0 atol = 0.001
    @test closest_point_to_pole(p5, 135.0).λ ≈ 170.2644 atol = 0.001
    @test closest_point_to_pole(p5, -135.0).ϕ ≈ -60.0 atol = 0.001
    @test closest_point_to_pole(p5, -135.0).λ ≈ -80.2644 atol = 0.001
    @test closest_point_to_pole(p5, -89.0).ϕ ≈ 45.00873 atol = 0.001
    @test closest_point_to_pole(p5, -89.0).λ ≈ 43.58593 atol = 0.001
    @test closest_point_to_pole(p5, -45.0).ϕ ≈ 60.0 atol = 0.001
    @test closest_point_to_pole(p5, -45.0).λ ≈ -9.73561 atol = 0.001

    p6 = Point(-45.0, 45.0)
    @test closest_point_to_pole(p6, 45.0).ϕ ≈ 60.0 atol = 0.001
    @test closest_point_to_pole(p6, 45.0).λ ≈ 170.2644 atol = 0.001
    @test closest_point_to_pole(p6, 91.0).ϕ ≈ -45.00873 atol = 0.001
    @test closest_point_to_pole(p6, 91.0).λ ≈ 46.41407 atol = 0.001
    @test closest_point_to_pole(p6, 135.0).ϕ ≈ -60.0 atol = 0.001
    @test closest_point_to_pole(p6, 135.0).λ ≈ 99.73561 atol = 0.001
    @test closest_point_to_pole(p6, -135.0).ϕ ≈ -60.0 atol = 0.001
    @test closest_point_to_pole(p6, -135.0).λ ≈ -9.73561 atol = 0.001
    @test closest_point_to_pole(p6, -91.0).ϕ ≈ -45.00873 atol = 0.001
    @test closest_point_to_pole(p6, -91.0).λ ≈ 43.58593 atol = 0.001
    @test closest_point_to_pole(p6, -45.0).ϕ ≈ 60.0 atol = 0.001
    @test closest_point_to_pole(p6, -45.0).λ ≈ -80.2644 atol = 0.001
end


text:
# A basic project can be created and accessed like so:

using REDCap

# Create the config object for project creation
super_config = REDCap.Config("<URL>", "<S-API>")

# `create_project()` returns the new projects config object, ready to use.
config = create_project(super_config, "Test Project", 1; purpose_other="Testing REDCap.jl Functionality", project_notes="This is not an actual REDCap Database.")

# ### Importing- NOTE: Records may be incomplete. Only provided fields will be updated
record=[Dict("sex"=>"0",
      "age"=>"56",
      "address"=>"168 Anderson Blvd. Quincy MA 01227",
      "height"=>"80",
      "dob"=>"1962-04-08",
      "record_id"=>"1",
      "bmi"=>"125",
      "comments"=>"Randomly Generated - Demographics",
      "email"=>"ALin@aol.com",
      "first_name"=>"Alexia",
      "demographics_complete"=>"0",
      "telephone"=>"(617) 882-6049",
      "weight"=>"80",
      "last_name"=>"Lin",
      "ethnicity"=>"1",
      "race"=>"1")]

import_records(config, record)

# Create new user with basic import/export permissions
user=[Dict("username" => "john_smith@email.com",
         "email" => "john_smith@email.com",
         "lastname" => "Smith",
         "api_export"=>"1",
         "api_import"=>"1")]

import_users(config, user)

# ### Exporting
records = export_records(config)

# `.pdf` summary of the project
export_pdf(config, "/<path>/export.pdf", allRecords=true)


text:
"Return first true `f` of `v` or `nothing`."
function getfirst(f::Function, v)
    i = findfirst(f,v)
    isnothing(i) ? nothing : v[i]
end


text:
using QuantumEvolution
using Test, Random

@testset "QuantumEvolution.jl" begin
    sol = optimize(x -> (x - 1)^2, i -> randn(), 50)
    @test sol.x ≈ 1 rtol = 1e-4
    sol =
        optimize(x -> (x[1] - 1)^2 + (x[2] - x[1]^2)^2, i -> randn(2), 50, max_iters = 500)
    map(sol.x) do x
        @test x ≈ 1 rtol = 1e-4
    end
    sol = optimize(x -> ifelse(1.2 < x, Inf, (x - 1)^2), i -> rand() * 0.2 + 0.9, 50)
    @test sol.x ≈ 1 rtol = 1e-4
    @test_throws ErrorException optimize(x -> Inf, i -> nothing, 50)
end


text:
using Dash, DashHtmlComponents, DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_input(id = "input-1-state", type = "text", value = "Montreal"),
    dcc_input(id = "input-2-state", type = "text", value = "Canada"),
    html_button(id = "submit-button-state", children = "submit", n_clicks = 0),
    html_div(id = "output-state")
end

callback!(
    app,
    Output("output-state", "children"),
    Input("submit-button-state", "n_clicks"),
    State("input-1-state", "value"),
    State("input-2-state", "value"),
) do clicks, input_1, input_2
    return "The Button has been pressed \"$clicks\" times, Input 1 is \"$input_1\" and Input 2 is \"$input_2\""
end

run_server(app, "0.0.0.0", debug=true)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
# This file contains the definition for the Renewables struct.

function get_res_dict(renewables_df::DataFrame)

   renewables_stacked = stack(renewables_df, Not(:node), variable_name=:fuel,
			           value_name=:capacity)

	renewables_stacked[!, :unit] = Symbol.(string.(renewables_stacked[!, :fuel]) .* "_"
		                      .* string.(renewables_stacked[!, :node]))

	renewables_stacked = renewables_stacked[:, [:unit, :node, :fuel, :capacity]]
	renewables_stacked = renewables_stacked[renewables_stacked[!, :capacity] .> 0, :]
	renewables_dict = df_to_dict_with_id(renewables_stacked)

	return renewables_stacked[!, :unit], renewables_dict
end

function get_avail_dict(renewables_avail_dict::Dict, renewables_dict)

	fuel = renewables_dict[:fuel]
	node = renewables_dict[:node]

	avail_dict = Dict{Symbol,Array{Float64, 1}}()
	for unit in keys(fuel)

		f = fuel[unit]
		n = Symbol(node[unit])

		if haskey(renewables_avail_dict, f)
			if Symbol(node[unit]) in names(renewables_avail_dict[f])
				avail_dict[unit] = renewables_avail_dict[f][!, n]
			else
				avail_dict[unit] = fill(0.0, 8760)
			end
		else
		 	avail_dict[unit] = renewables_avail_dict[:other][!, f]
		end
	end

	return avail_dict
end

function calc_infeed(cap_dict, avail_dict)

	infeed = map(collect(keys(cap_dict))) do unit
		cap = cap_dict[unit]
		avail = avail_dict[unit]

		unit => avail * cap
	end |> Dict

	return infeed
end

mutable struct Renewables
	unit
	fuel
	node
	capacity
	infeed

	function Renewables(renewables_df::DataFrame, renewables_avail_dict::Dict)

		units, renewables_dict = get_res_dict(renewables_df)
		avail_dict = get_avail_dict(renewables_avail_dict, renewables_dict)
		infeed = calc_infeed(renewables_dict[:capacity], avail_dict)

		return new(units, renewables_dict[:fuel], renewables_dict[:node], renewables_dict[:capacity], infeed)
	end
end

max_issues_repo_name: bobbyxiong/redispatch-ptg

max_forks_repo_name: bobbyxiong/redispatch-ptg


text:
#jl #! format: off
# # ACOPF with [PowerSimulations.jl](https://github.com/NREL-SIIP/PowerSimulations.jl) using [PowerModels.jl](https://github.com/lanl-ansi/PowerModels.jl)

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# PowerSimulations.jl supports non-linear AC optimal power flow through a deep integration
# with [PowerModels.jl](https://github.com/lanl-ansi/PowerModels.jl). This example shows a
# single multi-period optimization of economic dispatch with a full representation of
# AC optimal power flow. However, since we use a case where generators are subject to
# minimum operating points, we need to also execute a unit commitment problem to provide the
# ACOPF with a valid commitment pattern. This example uses a `Simulation` with two
# `DecisionModels` to execute the UC-ACOPF workflow for a single period.

# ## Dependencies
using SIIPExamples
using PowerSystems
using PowerSimulations
using PowerSystemCaseBuilder
using Dates
sim_folder = mktempdir(".", cleanup = true)

# We'll just use a suitable `System` that contains valid AC power flow parameters
sys = build_system(PSITestSystems, "modified_RTS_GMLC_DA_sys")
transform_single_time_series!(sys, 1, Hour(1))

# Since we'll be doing non-linear optimization, we need a solver that supports non-linear
# problems. Ipopt is quite good. And, we'll need a separate solver that can handle integer variables.
# So, we'll use HiGHS for the UC problem.
using Ipopt
solver = optimizer_with_attributes(Ipopt.Optimizer)
using HiGHS # mip solver
uc_solver = optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.5)

# Here, we want do define an economic dispatch (linear generation decisions) with an ACOPF
# network representation.
# So, starting with the network, we can select from _almost_ any of the endpoints on this
# tree:
print_tree(PowerSimulations.PM.AbstractPowerModel)

# First, we can setup a template with a suitable ACOPF network formulation, and formulations
# that represent each of the relevant device categories
ed_template = ProblemTemplate()
set_device_model!(ed_template, ThermalStandard, ThermalStandardDispatch)
set_device_model!(ed_template, PowerLoad, StaticPowerLoad)
set_device_model!(ed_template, Line, StaticBranch)
set_device_model!(ed_template, TapTransformer, StaticBranch)
set_device_model!(ed_template, Transformer2W, StaticBranch)
set_device_model!(ed_template, HVDCLine, HVDCDispatch)
set_network_model!(ed_template, NetworkModel(ACPPowerModel, use_slacks = true))

# We also need to setup a UC template with a simplified network representation
uc_template = ProblemTemplate(DCPPowerModel)
set_device_model!(uc_template, ThermalStandard, ThermalBasicUnitCommitment)
set_device_model!(uc_template, PowerLoad, StaticPowerLoad)
set_device_model!(uc_template, Line, StaticBranch)
set_device_model!(uc_template, TapTransformer, StaticBranch)
set_device_model!(uc_template, Transformer2W, StaticBranch)
set_device_model!(uc_template, HVDCLine, HVDCDispatch)
set_service_model!(uc_template, VariableReserve{ReserveUp}, RangeReserve)

# Now we can build a simulation to solve the UC, pass the commitment pattern to the ACOPF
# and then solve the ACOPF.
models = SimulationModels(
    decision_models = [
        DecisionModel(uc_template, sys, name = "UC", optimizer = uc_solver),
        DecisionModel(
            ed_template,
            sys,
            name = "ACOPF",
            optimizer = solver,
            initialize_model = false,
        ),
    ],
)
sequence = SimulationSequence(
    models = models,
    feedforwards = Dict(
        "ACOPF" => [
            SemiContinuousFeedforward(
                component_type = ThermalStandard,
                source = OnVariable,
                affected_values = [ActivePowerVariable, ReactivePowerVariable],
            ),
        ],
    ),
    ini_cond_chronology = InterProblemChronology(),
)

# Note that in the above feedforward definition, the `OnVariable` for the `ThermalStandard`
# components is affecting both the `ActivePowerVariable` and the `ReactivePowerVariable`.
# This is the connection that restricts the ACOPF to only represent active and reactive
# power injections from the units that are committed in the UC problem. It's not guaranteed
# that the UC result generates an AC feasible initial condition, so this problem selects
# a particular period (`initial_time`) where the conditions are suitable.
sim = Simulation(
    name = "UC-ACOPF",
    steps = 12,
    models = models,
    sequence = sequence,
    simulation_folder = sim_folder,
)

build!(sim)

# And solve it ...
execute!(sim, enable_progress_bar = false)

# And extract some results

results = SimulationResults(sim)
ac_results = get_decision_problem_results(results, "ACOPF")

slack_keys = [
    k for k in list_variable_keys(ac_results) if
    PSI.get_entry_type(k) ∈ [SystemBalanceSlackDown, SystemBalanceSlackUp]
]
slack_vars = Dict([k => read_realized_variables(ac_results, k) for k in slack_keys])

#nb # Plot the slack values
#nb plot_results(slack_vars);

max_stars_repo_name: NREL-SIIP/Examples

max_issues_repo_name: NREL-SIIP/Examples

max_forks_repo_name: NREL-SIIP/Examples


text:
using Oceananigans.Operators

"""
    struct MixingLength{FT}

Contains mixing length parameters for CATKE vertical diffusivity.

The mixing length is the maximum of three mixing lengths:

```math
ℓᵩ = max(ℓᵟᵩ, ℓᴬᵩ, ℓ⋆ᵩ)
```

where ``ℓᵟ`` is a grid-spacing-proportional length,
``ℓʰᵩ`` is a convective mixing length, and ``ℓ⋆ᵩ`` is
a mixing length in stably-stratified mixing.

Grid-spacing-proportional mixing length
=======================================

```math
ℓᵟᵩ = Cᵟϕ * Δz
```

Convective mixing length
========================

```math
ℓᴬᵩ = convecting ? Cᴬϕ * e^3/2 / Qᵇ : 0
```

Stably-stratified mixing length
===============================

```math
ℓ⋆ᵩ = σᵩ * min(d, Cᴸᵇ * √e / N)
```

where ``σᵩ(Ri)`` is a stability function that depends on
the local Richardson number, ...

The Richardson number is

    ``Ri = ∂z B / ( (∂z U)² + (∂z V)² )`` ,

where ``B`` is buoyancy and ``∂z`` denotes a vertical derviative.
The Richardson-number dependent diffusivities are multiplied by the stability
function

    1. ``σ(Ri) = σ⁻ * (1 + rσ * step(Ri, Riᶜ, Riʷ))``
    3. ``σ(Ri) = σ⁻ + (σ⁺ - σ⁻) * step(Ri, Riᶜ, Riʷ)``

σ⁻ = σ₀
rσ = (σ⁺ - σ⁻) / σ₀

where ``σ₀``, ``Δσ``, ``Riᶜ``, and ``Riʷ`` are free parameters,
and ``step`` is a smooth step function defined by

    ``step(x, c, w) = (1 + \tanh((x - c) / w)) / 2``.

The 8 free parameters in `RiDependentDiffusivityScaling` have been _experimentally_ calibrated
against large eddy simulations of ocean surface boundary layer turbulence in idealized
scenarios involving monotonic boundary layer deepening into variable stratification
due to constant surface momentum fluxes and/or destabilizing surface buoyancy flux.
See https://github.com/CliMA/LESbrary.jl for more information about the large eddy simulations.
The calibration was performed using a combination of Markov Chain Monte Carlo (MCMC)-based simulated
annealing and noisy Ensemble Kalman Inversion methods.
"""
Base.@kwdef struct MixingLength{FT}
    Cᴸᵇ   :: FT = 1.16
    Cᵟu   :: FT = 0.5
    Cᵟc   :: FT = 0.5
    Cᵟe   :: FT = 0.5
    Cᴬu   :: FT = 0.0
    Cᴬc   :: FT = 0.0
    Cᴬe   :: FT = 0.0
    Cᴷu⁻  :: FT = 0.15
    Cᴷuʳ  :: FT = 3.87
    Cᴷc⁻  :: FT = 0.40
    Cᴷcʳ  :: FT = 0.77
    Cᴷe⁻  :: FT = 0.13
    Cᴷeʳ  :: FT = 1.11
    CᴷRiʷ :: FT = 0.72
    CᴷRiᶜ :: FT = 0.76
end

#####
##### Mixing length
#####

@inline surface(i, j, k, grid)                = znode(Center(), Center(), Face(), i, j, grid.Nz+1, grid)
@inline bottom(i, j, k, grid)                 = znode(Center(), Center(), Face(), i, j, 1, grid)
@inline depthᶜᶜᶜ(i, j, k, grid)               = surface(i, j, k, grid) - znode(Center(), Center(), Center(), i, j, k, grid)
@inline height_above_bottomᶜᶜᶜ(i, j, k, grid) = znode(Center(), Center(), Center(), i, j, k, grid) - bottom(i, j, k, grid)

@inline wall_vertical_distanceᶜᶜᶜ(i, j, k, grid) = min(depthᶜᶜᶜ(i, j, k, grid), height_above_bottomᶜᶜᶜ(i, j, k, grid))

@inline function sqrt_∂z_b(i, j, k, grid, buoyancy, tracers)
    FT = eltype(grid)
    N² = ∂z_b(i, j, k, grid, buoyancy, tracers)
    N²⁺ = max(zero(FT), N²)
    return sqrt(N²⁺)
end

@inline function buoyancy_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy)
    FT = eltype(grid)
    Cᴸᵇ = closure.mixing_length.Cᴸᵇ
    N⁺ = ℑzᵃᵃᶜ(i, j, k, grid, sqrt_∂z_b, buoyancy, tracers)

    @inbounds e⁺ = max(zero(FT), e[i, j, k])

    return @inbounds ifelse(N⁺ == 0, FT(Inf), Cᴸᵇ * sqrt(e⁺) / N⁺)
end

@inline function unscaled_stable_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy)
    d = wall_vertical_distanceᶜᶜᶜ(i, j, k, grid)
    ℓᵇ = buoyancy_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy)
    return min(d, ℓᵇ)
end

@inline function unscaled_convective_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs)
    Qᵇ = top_buoyancy_flux(i, j, grid, buoyancy, tracer_bcs, clock, merge(velocities, tracers))
    N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers)
    convecting = (N² < 0) & (Qᵇ > 0)
    @inbounds eijk = max(zero(eltype(grid)), tracers.e[i, j, k])
    return ifelse(convecting, sqrt(eijk^3) / Qᵇ, zero(eltype(grid)))
end

@inline ϕ²(i, j, k, grid, ϕ, args...) = ϕ(i, j, k, grid, args...)^2

@inline function Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy)
    FT = eltype(grid)
    ∂z_u² = ℑxzᶜᵃᶜ(i, j, k, grid, ϕ², ∂zᵃᵃᶠ, velocities.u)
    ∂z_v² = ℑyzᵃᶜᶜ(i, j, k, grid, ϕ², ∂zᵃᵃᶠ, velocities.v)
    N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers)
    return ifelse(N² == 0, zero(FT), N² / (∂z_u² + ∂z_v²))
end

@inline step(x, c, w) = (1 + tanh((x - c) / w)) / 2

@inline scale(Ri, σ⁻, rσ, c, w) = σ⁻ * (1 + rσ * step(Ri, c, w))

@inline function momentum_stable_mixing_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    Ri = Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy)
    return scale(Ri,
                 closure.mixing_length.Cᴷu⁻,
                 closure.mixing_length.Cᴷuʳ,
                 closure.mixing_length.CᴷRiᶜ,
                 closure.mixing_length.CᴷRiʷ)
end

@inline function tracer_stable_mixing_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    Ri = Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy)
    return scale(Ri,
                 closure.mixing_length.Cᴷc⁻,
                 closure.mixing_length.Cᴷcʳ,
                 closure.mixing_length.CᴷRiᶜ,
                 closure.mixing_length.CᴷRiʷ)
end

@inline function TKE_stable_mixing_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    Ri = Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy)
    return scale(Ri,
                 closure.mixing_length.Cᴷe⁻,
                 closure.mixing_length.Cᴷeʳ,
                 closure.mixing_length.CᴷRiᶜ,
                 closure.mixing_length.CᴷRiʷ)
end

@inline function momentum_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs)
    Cᴬu = closure.mixing_length.Cᴬu
    ℓᴬ = unscaled_convective_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs)

    Cᵟu = closure.mixing_length.Cᵟu
    ℓᵟ = Δzᵃᵃᶜ(i, j, k, grid)

    σu = momentum_stable_mixing_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    ℓ★ = unscaled_stable_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, tracers.e, tracers, buoyancy)

    return max(Cᴬu * ℓᴬ, Cᵟu * ℓᵟ, σu * ℓ★)
end

@inline function tracer_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs)
    Cᴬc = closure.mixing_length.Cᴬc
    ℓᴬ = unscaled_convective_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs)

    Cᵟc = closure.mixing_length.Cᵟc
    ℓᵟ = Δzᵃᵃᶜ(i, j, k, grid)

    σc = tracer_stable_mixing_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    ℓ★ = unscaled_stable_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, tracers.e, tracers, buoyancy)

    return max(Cᴬc * ℓᴬ, Cᵟc * ℓᵟ, σc * ℓ★)
end

@inline function TKE_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs)
    Cᴬe = closure.mixing_length.Cᴬe
    ℓᴬ = unscaled_convective_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs)

    Cᵟe = closure.mixing_length.Cᵟe
    ℓᵟ = Δzᵃᵃᶜ(i, j, k, grid)

    σe = TKE_stable_mixing_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    ℓ★ = unscaled_stable_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, tracers.e, tracers, buoyancy)

    return max(Cᴬe * ℓᴬ, Cᵟe * ℓᵟ, σe * ℓ★)
end


text:
abstract type AbstractFiniteDifferenceSpace <: AbstractSpace end

abstract type Staggering end

""" Cell center location """
struct CellCenter <: Staggering end

""" Cell face location """
struct CellFace <: Staggering end

struct FiniteDifferenceSpace{
    S <: Staggering,
    T <: Topologies.IntervalTopology,
    GG,
    LG,
} <: AbstractFiniteDifferenceSpace
    staggering::S
    topology::T
    global_geometry::GG
    center_local_geometry::LG
    face_local_geometry::LG
end

function Base.show(io::IO, space::FiniteDifferenceSpace)
    indent = get(io, :indent, 0)
    iio = IOContext(io, :indent => indent + 2)
    println(
        io,
        space isa CenterFiniteDifferenceSpace ? "CenterFiniteDifferenceSpace" :
        "FaceFiniteDifferenceSpace",
        ":",
    )
    print(iio, " "^(indent + 2), space.topology)
end

function FiniteDifferenceSpace{S}(
    topology::Topologies.IntervalTopology,
) where {S <: Staggering}
    global_geometry = Geometry.CartesianGlobalGeometry()
    mesh = topology.mesh
    CT = Meshes.coordinate_type(mesh)
    AIdx = Geometry.coordinate_axis(CT)
    # TODO: FD operators  hardcoded to work over the 3-axis, need to generalize
    # similar to spectral operators
    @assert AIdx == (3,) "FiniteDifference operations only work over the 3-axis (ZPoint) domain"
    FT = eltype(CT)
    face_coordinates = collect(mesh.faces)
    LG = Geometry.LocalGeometry{AIdx, CT, FT, SMatrix{1, 1, FT, 1}}
    nface = length(face_coordinates)
    ncent = nface - 1
    center_local_geometry = DataLayouts.VF{LG}(Array{FT}, ncent)
    face_local_geometry = DataLayouts.VF{LG}(Array{FT}, nface)
    for i in 1:ncent
        # centers
        coord⁻ = Geometry.component(face_coordinates[i], 1)
        coord⁺ = Geometry.component(face_coordinates[i + 1], 1)
        # at the moment we use a "discrete Jacobian"
        # ideally we should use the continuous quantity via the derivative of the warp function
        # could we just define this then as deriv on the mesh element coordinates?
        coord = (coord⁺ + coord⁻) / 2
        Δcoord = coord⁺ - coord⁻
        J = Δcoord
        WJ = Δcoord
        ∂x∂ξ = SMatrix{1, 1}(J)
        center_local_geometry[i] = Geometry.LocalGeometry(
            CT(coord),
            J,
            WJ,
            Geometry.AxisTensor(
                (Geometry.LocalAxis{AIdx}(), Geometry.CovariantAxis{AIdx}()),
                ∂x∂ξ,
            ),
        )
    end
    for i in 1:nface
        coord = Geometry.component(face_coordinates[i], 1)
        if i == 1
            # bottom face
            coord⁺ = Geometry.component(face_coordinates[2], 1)
            J = coord⁺ - coord
            WJ = J / 2
        elseif i == nface
            # top face
            coord⁻ = Geometry.component(face_coordinates[i - 1], 1)
            J = coord - coord⁻
            WJ = J / 2
        else
            coord⁺ = Geometry.component(face_coordinates[i + 1], 1)
            coord⁻ = Geometry.component(face_coordinates[i - 1], 1)
            J = (coord⁺ - coord⁻) / 2
            WJ = J
        end
        ∂x∂ξ = SMatrix{1, 1}(J)
        ∂ξ∂x = SMatrix{1, 1}(inv(J))
        face_local_geometry[i] = Geometry.LocalGeometry(
            CT(coord),
            J,
            WJ,
            Geometry.AxisTensor(
                (Geometry.LocalAxis{AIdx}(), Geometry.CovariantAxis{AIdx}()),
                ∂x∂ξ,
            ),
        )
    end
    return FiniteDifferenceSpace(
        S(),
        topology,
        global_geometry,
        center_local_geometry,
        face_local_geometry,
    )
end

FiniteDifferenceSpace{S}(mesh::Meshes.IntervalMesh) where {S <: Staggering} =
    FiniteDifferenceSpace{S}(Topologies.IntervalTopology(mesh))

const CenterFiniteDifferenceSpace = FiniteDifferenceSpace{CellCenter}
const FaceFiniteDifferenceSpace = FiniteDifferenceSpace{CellFace}

function FiniteDifferenceSpace{S}(
    space::FiniteDifferenceSpace,
) where {S <: Staggering}
    FiniteDifferenceSpace(
        S(),
        space.topology,
        space.global_geometry,
        space.center_local_geometry,
        space.face_local_geometry,
    )
end

Base.length(space::FiniteDifferenceSpace) = length(coordinates_data(space))

nlevels(space::FiniteDifferenceSpace) = length(space)

local_geometry_data(space::CenterFiniteDifferenceSpace) =
    space.center_local_geometry

local_geometry_data(space::FaceFiniteDifferenceSpace) =
    space.face_local_geometry

left_boundary_name(space::FiniteDifferenceSpace) =
    propertynames(space.topology.boundaries)[1]

right_boundary_name(space::FiniteDifferenceSpace) =
    propertynames(space.topology.boundaries)[2]


text:
using Documenter, DocPkg

makedocs(;
    modules=[DocPkg],
    format=Documenter.HTML(
        # prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = Documenter.MathJax()
    ),
    pages=[
        "Home" => "index.md",
        "Updated" => "updated.md",
    ],
    strict = true,  # See https://github.com/JuliaOpt/JuMP.jl/issues/1576
    repo="https://github.com/AKS1996/DocPkg.jl/blob/{commit}{path}#L{line}",
    sitename="DocPkg.jl",
    authors="Akshay Sharma"
)

deploydocs(
    repo   = "github.com/AKS1996/DocPkg.jl.git",
)


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Axis Aligned Bounding Box

"""
    AABB

A type representing an Axis-Aligned Bounding Box.
"""
Base.@kwdef struct AABB
    p_M::Point = Point(1f0, 1f0, 1f0)
    p_m::Point = Point(0f0, 0f0, 0f0)
end

"""
    get_t(ray::Ray, aabb::AABB)

Return the parameter `t` at which [`Ray`](@ref) first hits the [`AABB`](@ref). If no hit exists, return `Inf32`.
"""
function get_t(ray::Ray, aabb::AABB)
    dir = ray.dir
    o = ray.origin
    min_point = replace((aabb.p_m - o) ./ dir, NaN32 => 0)
    max_point = replace((aabb.p_M - o) ./ dir, NaN32 => 0)
    overlap = reduce(intersect, map(t -> Interval(extrema(t)...), zip(min_point, max_point)))
    isempty(overlap) && return Inf32
    t1, t2 = overlap.first, overlap.last
    ray.tmin < t1 < ray.tmax && return t1
    ray.tmin < t2 < ray.tmax && return t2
    return Inf32
end

function get_all_ts(ray::Ray, aabb::AABB)
    dir = ray.dir
    o = ray.origin
    min_point = replace((aabb.p_m - o) ./ dir, NaN32 => 0)
    max_point = replace((aabb.p_M - o) ./ dir, NaN32 => 0)
    overlap = reduce(intersect, map(t -> Interval(extrema(t)...), zip(min_point, max_point)))
    isempty(overlap) && return Vector{Float32}()
    t1, t2 = overlap.first, overlap.last
    isfinite(t1) && (@assert isfinite(t2); return [t1, t2])
    @assert !isfinite(t2)
    return Vector{Float32}()
end


text:
export rydberg_tebd!

function cterm(C::Real)
    return C*kron(P1, P1)
end

function zterm(Δ::Real)
    return Δ*Z
end

function xterm(Ω::Real)
    return Ω/2*X
end

function rydberg_tebd_sweep!(peps, g, C::Real, Δ::Real, Ω::Real, δt)
    for edge in edges(g)
        di = degree(g, edge.src)
        dj = degree(g, edge.dst)
        hi = cterm(C) + put(2, 1=>xterm(Ω/di)) + put(2, 2=>xterm(Ω/dj)) + put(2, 1=>zterm(Δ/di)) + put(2, 2=>zterm(Δ/dj))
        apply!(peps, put(nqubits(peps), (edge.src, edge.dst)=>time_evolve(hi, δt)))
    end
    return peps
end

function rydberg_tebd!(reg, g::SimpleGraph; t::Real, C::Real, Δ::Real, Ω::Real, nstep::Int)
    δt = t/nstep
    for _ = 1:nstep-1
        rydberg_tebd_sweep!(reg, g, C, Δ, Ω, δt)
    end
    return reg
end

function rydberg_tebd!(peps::PEPS; t::Real, C::Real, Δ::Real, Ω::Real, nstep::Int)
    g = SimpleGraph(nqubits(peps))
    for (i, j) in virtualbonds(peps)
        add_edge!(g, i, j)
    end
    rydberg_tebd!(peps, g; t=t, C=C, Δ=Δ, Ω=Ω, nstep=nstep)
end