text:
###############
# H-OBCA: Hierarchical Optimization-based Collision Avoidance - a path planner for autonomous parking
# Copyright (C) 2018
# Atsushi SAKAI [atsushisakai@global.komatsu; Komatsu Ltd / MPC Lab]
# Alexander LINIGER [liniger@control.ee.ethz.ch; Automatic Control Lab, ETH Zurich]
# Xiaojing ZHANG [xiaojing.zhang@berkeley.edu; MPC Lab, UC Berkeley]
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
###############
# The paper describing the theory can be found here:
# 	X. Zhang, A. Liniger and F. Borrelli; "Optimization-Based Collision Avoidance"; Technical Report, 2017
#   X. Zhang, A. Liniger, A. Sakai and F. Borrelli; "Autonomous  Parking  using  Optimization-Based  Collision  Avoidance"; Technical Report, 2018 [add URL]
###############

###############
# Hybrid A star: Julia implementation of Hybrid A* algorithm
###############

module hybrid_a_star

using PyPlot
using DataFrames
using NearestNeighbors
using DataStructures

include("./reeds_shepp.jl")
include("./a_star.jl")
include("./collision_check.jl")

const VEHICLE_RADIUS = 1.0 #[m]; radius of rear ball;
const BUBBLE_DIST = 1.7 #[m]; distance to "forward bubble";

##### Fast Comp Time values from Alex Liniger ######
const OB_MAP_RESOLUTION = 0.1 #[m]; obstacle resolution
const YAW_GRID_RESOLUTION = deg2rad(5.0) #[m]; 10.0 /// 5.0
const N_STEER = 5.0 # number of steer command; 10.0 seems OK /// 5
## For Backwards Parking
const XY_GRID_RESOLUTION = 0.3 #[m];
const MOTION_RESOLUTION = 0.1 #[m];
###################################################

const USE_HOLONOMIC_WITH_OBSTACLE_HEURISTIC = true
const USE_NONHOLONOMIC_WITHOUT_OBSTACLE_HEURISTIC = false

const SB_COST = 10.0 # switch back penalty cost
const BACK_COST = 0.0 # backward penalty cost
const STEER_CHANGE_COST = 10.0 # steer angle change penalty cost
const STEER_COST = 0.0  # steer angle  penalty cost
const H_COST = 1. # Heuristic cost; higher -> heuristic; 1.0

const WB = 2.7 #[m]; 7.0
const MAX_STEER = 0.6  #[rad]

type Node
    xind::Int64 #x index
    yind::Int64 #y index
    yawind::Int64 #yaw index
    direction::Bool # moving direction forword:true, backword:false
    x::Array{Float64} # x position [m]
    y::Array{Float64} # y position [m]
    yaw::Array{Float64} # yaw angle [rad]
    steer::Float64 # steer input
    cost::Float64 # cost
    pind::Int64 # parent index
end

type Config
    minx::Int64
    miny::Int64
    minyaw::Int64
    maxx::Int64
    maxy::Int64
    maxyaw::Int64
    xw::Int64
    yw::Int64
    yaww::Int64
    xyreso::Float64
    yawreso::Float64
    obminx::Int64
    obminy::Int64
    obmaxx::Int64
    obmaxy::Int64
    obxw::Int64
    obyw::Int64
    obreso::Float64
end

function calc_hybrid_astar_path(sx::Float64, sy::Float64, syaw::Float64,
                                gx::Float64, gy::Float64, gyaw::Float64,
                                ox::Array{Float64}, oy::Array{Float64},
                                xyreso::Float64, yawreso::Float64,
                                obreso::Float64)
    """
    Calc hybrid astar path
    sx: start x position [m]
    sy: start y position [m]
    gx: goal x position [m]
    gx: goal x position [m]
    ox: x position list of Obstacles [m]
    oy: y position list of Obstacles [m]
    xyreso: grid resolution [m]
    yawreso: yaw angle resolution [rad]
    """

    syaw, gyaw = pi_2_pi(syaw), pi_2_pi(gyaw)

    const c = calc_config(ox, oy, xyreso, yawreso, obreso)
    kdtree = KDTree(hcat(ox, oy)')
    obmap, gkdtree = calc_obstacle_map(ox, oy, c)
    nstart = Node(round(Int64,sx/xyreso), round(Int64,sy/xyreso), round(Int64, syaw/yawreso),true,[sx],[sy],[syaw],0.0,0.0, -1)
    ngoal = Node(round(Int64,gx/xyreso), round(Int64,gy/xyreso), round(Int64,gyaw/yawreso),true,[gx],[gy],[gyaw],0.0,0.0, -1)

    if USE_HOLONOMIC_WITH_OBSTACLE_HEURISTIC
        h_dp = calc_holonomic_with_obstacle_heuristic(ngoal, ox, oy, xyreso)
    else
        h_dp = Array{Float64}()
    end
    if USE_NONHOLONOMIC_WITHOUT_OBSTACLE_HEURISTIC
        h_rs = calc_nonholonomic_without_obstacle_heuristic(ngoal, c)
    else
        h_rs = Array{Float64}()
    end

    open, closed = Dict{Int64, Node}(), Dict{Int64, Node}()
    open[calc_index(nstart, c)] = nstart
    pq = PriorityQueue()
    enqueue!(pq, calc_index(nstart, c), calc_cost(nstart, h_rs, h_dp, ngoal, c))

    u, d = calc_motion_inputs()
    nmotion = length(u)

    while true
        if length(open) == 0
            println("Error: Cannot find path, No open set")
            return nothing, nothing, nothing
        end

        c_id = dequeue!(pq)
        current = open[c_id]

        isupdated, current = update_node_with_analystic_expantion(current, ngoal, obmap, c, kdtree, ox, oy)
        if isupdated
            closed[calc_index(ngoal, c)] = current
            break #goal
        end

        #move current node from open to closed
        delete!(open, c_id)
        closed[c_id] = current

        for i in 1:nmotion
            node = calc_next_node(current, c_id, u[i], d[i], c, gkdtree)

            if !verify_index(node, obmap, c, kdtree, ox, oy) continue end

            node_ind = calc_index(node, c)

            # If it is already in the closed set, skip it
            if haskey(closed, node_ind)  continue end

            if !haskey(open, node_ind)
                open[node_ind] = node
                enqueue!(pq, calc_index(node, c), calc_cost(node, h_rs, h_dp, ngoal, c))
            end
        end

    end

    # println("final expand node:", length(open) + length(closed))

    rx, ry, ryaw = get_final_path(closed, ngoal, nstart, c)

    return rx, ry, ryaw
end

function update_node_with_analystic_expantion(current::Node,
                                             ngoal::Node,
                                             obmap::Array{Bool,2},
                                             c::Config,
                                             kdtree::NearestNeighbors.KDTree,
                                             ox::Array{Float64},
                                             oy::Array{Float64}
                                            )

    apath = analystic_expantion(current, ngoal, obmap, c, kdtree, ox, oy)
    if apath != nothing
        # println("Find path! with analystic_expantion")
        current.x = vcat(current.x, apath.x[2:end-1])
        current.y = vcat(current.y, apath.y[2:end-1])
        current.yaw = vcat(current.yaw, apath.yaw[2:end-1])
        current.cost += calc_rs_path_cost(apath)
        return true, current
    end

    return false, current #no update
end

function calc_rs_path_cost(rspath::hybrid_a_star.reeds_shepp.Path)

    cost = 0.0
    for l in rspath.lengths
        if l >= 0 # forward
            cost += l
        else # back
            cost += abs(l) * BACK_COST
        end
    end

    # swich back penalty
    for i in 1:length(rspath.lengths) - 1
        if rspath.lengths[i] * rspath.lengths[i+1] < 0.0 # switch back
            cost += SB_COST
        end
    end

    # steer penalyty
    for ctype in rspath.ctypes
        if ctype != "S" # curve
            cost += STEER_COST*abs(MAX_STEER)
        end
    end

    # ==steer change penalty
    # calc steer profile
    nctypes = length(rspath.ctypes)
    ulist = fill(0.0, nctypes)
    for i in 1:nctypes
        if rspath.ctypes[i] == "R"
            ulist[i] = - MAX_STEER
        elseif rspath.ctypes[i] == "L"
            ulist[i] = MAX_STEER
        end
    end

    for i in 1:length(rspath.ctypes) - 1
        cost += STEER_CHANGE_COST*abs(ulist[i+1] - ulist[i])
    end

    # println("RS cost is ", cost)
    return cost
end

function analystic_expantion(n::Node, ngoal::Node, obmap::Array{Bool,2}, c::Config,
                             kdtree::NearestNeighbors.KDTree,
                             ox::Array{Float64},
                             oy::Array{Float64}
                            )

    sx = n.x[end]
    sy = n.y[end]
    syaw = n.yaw[end]

    max_curvature = tan(MAX_STEER)/WB
    path = reeds_shepp.calc_shortest_path(sx, sy, syaw,
                                           ngoal.x[end], ngoal.y[end], ngoal.yaw[end],
                                           max_curvature, step_size=MOTION_RESOLUTION)

    if path == nothing
        return nothing
    end

    if !collision_check.check_collision(path.x, path.y, path.yaw, kdtree, ox, oy)
        return nothing
    end

    # println(paths)
    return path # find good path
end

function calc_motion_inputs()

    up = [i for i in MAX_STEER/N_STEER:MAX_STEER/N_STEER:MAX_STEER]
    u = vcat([0.0], [i for i in up], [-i for i in up])
    d = vcat([1.0 for i in 1:length(u)], [-1.0 for i in 1:length(u)])
    u = vcat(u,u)

    return u, d
end

function verify_index(node::Node, obmap::Array{Bool,2}, c::Config,
                      kdtree::NearestNeighbors.KDTree,
                      ox::Array{Float64},
                      oy::Array{Float64}
                     )::Bool

    # overflow map
    if (node.xind - c.minx) >= c.xw
        return false
    elseif (node.xind - c.minx) <= 0
        return false
    end
    if (node.yind - c.miny) >= c.yw
        return false
    elseif (node.yind - c.miny) <= 0
        return false
    end

    # check collisiton
    # rectangle check
    if !collision_check.check_collision(node.x, node.y,node.yaw, kdtree, ox, oy)
        return false
    end

    return true #index is ok"
end

function pi_2_pi(iangle::Float64)
    while (iangle > pi)
        iangle -= 2.0 * pi
    end
    while (iangle < -pi)
        iangle += 2.0 * pi
    end

    return iangle
end

function calc_next_node(current::Node, c_id::Int64,
                        u::Float64, d::Float64,
                        c::Config,
                        gkdtree::NearestNeighbors.KDTree)

    arc_l = XY_GRID_RESOLUTION

    nlist = round(Int64, arc_l/MOTION_RESOLUTION)+1
    xlist = fill(0.0, nlist)
    ylist = fill(0.0, nlist)
    yawlist = fill(0.0, nlist)
    xlist[1] = current.x[end] + d * MOTION_RESOLUTION*cos(current.yaw[end])
    ylist[1] = current.y[end] + d * MOTION_RESOLUTION*sin(current.yaw[end])
    yawlist[1] = pi_2_pi(current.yaw[end] + d*MOTION_RESOLUTION/WB * tan(u))

    for i in 1:(nlist - 1)
        xlist[i+1] = xlist[i] + d * MOTION_RESOLUTION*cos(yawlist[i])
        ylist[i+1] = ylist[i] + d * MOTION_RESOLUTION*sin(yawlist[i])
        yawlist[i+1] = pi_2_pi(yawlist[i] + d*MOTION_RESOLUTION/WB * tan(u))
    end

    xind = round(Int64, xlist[end]/c.xyreso)
    yind = round(Int64, ylist[end]/c.xyreso)
    yawind = round(Int64, yawlist[end]/c.yawreso)

    addedcost = 0.0
    if d > 0
        direction = true
        addedcost += abs(arc_l)
    else
        direction = false
        addedcost += abs(arc_l) * BACK_COST
    end

    # swich back penalty
    if direction != current.direction # switch back penalty
        addedcost += SB_COST
    end

    # steer penalyty
    addedcost += STEER_COST*abs(u)

    # steer change penalty
    addedcost += STEER_CHANGE_COST*abs(current.steer - u)

    cost = current.cost + addedcost
    node = Node(xind, yind, yawind, direction, xlist, ylist, yawlist, u, cost, c_id)
    # println(node)

    return node
end

function is_same_grid(node1::Node,node2::Node)

    if node1.xind != node2.xind
        return false
    end
    if node1.yind != node2.yind
        return false
    end
    if node1.yawind != node2.yawind
        return false
    end

    return true

end

function calc_index(node::Node, c::Config)
    ind = (node.yawind - c.minyaw)*c.xw*c.yw+(node.yind - c.miny)*c.xw + (node.xind - c.minx)
    if ind <= 0
        println("Error(calc_index):", ind)
    end
    return ind
end

function calc_holonomic_with_obstacle_heuristic(gnode::Node, ox::Array{Float64}, oy::Array{Float64}, xyreso::Float64)
    # println("Calc distance policy")
    h_dp = a_star.calc_dist_policy(gnode.x[end], gnode.y[end], ox, oy, xyreso, VEHICLE_RADIUS)
    return h_dp
end

function calc_nonholonomic_without_obstacle_heuristic(ngoal::Node,
                                                      c::Config)

    h_rs = fill(0.0, (c.xw,c.yw,c.yaww))
    max_curvature = tan(MAX_STEER)/WB

    for ix in 1:c.xw
        for iy in 1:c.yw
            for iyaw in 1:c.yaww
                sx = (ix + c.minx)*c.xyreso
                sy = (iy + c.miny)*c.xyreso
                syaw = pi_2_pi((iyaw + c.minyaw)*c.yawreso)
                L = reeds_shepp.calc_shortest_path_length(sx, sy, syaw,
                                                          ngoal.x[end], ngoal.y[end], ngoal.yaw[end],
                                                          max_curvature, step_size=MOTION_RESOLUTION)
                h_rs[ix, iy, iyaw] = L
            end
        end
    end

    # println(h_rs[:,:,1])

    return h_rs
end

function calc_config(ox::Array{Float64}, oy::Array{Float64}, xyreso::Float64, yawreso::Float64, obreso::Float64)
    minx = round(Int64, minimum(ox)/xyreso)
    miny = round(Int64, minimum(oy)/xyreso)
    maxx = round(Int64, maximum(ox)/xyreso)
    maxy = round(Int64, maximum(oy)/xyreso)
    obminx = round(Int64, minimum(ox)/obreso)
    obminy = round(Int64, minimum(oy)/obreso)
    obmaxx = round(Int64, maximum(ox)/obreso)
    obmaxy = round(Int64, maximum(oy)/obreso)
    # println("minx:", minx)
    # println("miny:", miny)
    # println("maxx:", maxx)
    # println("maxy:", maxy)

    xw = round(Int64,(maxx - minx))
    yw = round(Int64,(maxy - miny))
    obxw = round(Int64,(obmaxx - obminx))
    obyw = round(Int64,(obmaxy - obminy))

    minyaw = round(Int64, - pi/yawreso) - 1
    maxyaw = round(Int64, pi/yawreso)
    yaww = round(Int64,(maxyaw - minyaw))

    config = Config(minx, miny, minyaw, maxx, maxy, maxyaw, xw, yw, yaww,
                    xyreso, yawreso, obminx, obminy, obmaxx, obmaxy, obxw, obyw, obreso)

    return config
end

function calc_obstacle_map(ox::Array{Float64},
                           oy::Array{Float64},
                           c::Config)

    ox = [iox/c.obreso for iox in ox]
    oy = [ioy/c.obreso for ioy in oy]

    obmap = fill(false, (c.obxw, c.obyw))

    gkdtree = KDTree(hcat(ox, oy)')
    for ix in 1:c.obxw
        x = ix + c.obminx
        for iy in 1:c.obyw
            y = iy + c.obminy
            idxs, onedist = knn(gkdtree, [x, y] , 1)
            if onedist[1] <= VEHICLE_RADIUS/c.obreso
                obmap[ix,iy] = true
            end
        end
    end

    return obmap, gkdtree
end

function get_final_path(closed::Dict{Int64, Node},
                        ngoal::Node,
                        nstart::Node,
                        c::Config)

    rx, ry, ryaw = Array{Float64}(ngoal.x),Array{Float64}(ngoal.y),Array{Float64}(ngoal.yaw)
    nid = calc_index(ngoal, c)
    # println("Final cost is ", closed[nid].cost)
    while true
        n = closed[nid]
        rx = vcat(rx, reverse(n.x))
        ry = vcat(ry, reverse(n.y))
        ryaw = vcat(ryaw, reverse(n.yaw))
        nid = n.pind
        if is_same_grid(n, nstart)
            # println("done")
            break
        end
    end

    rx = reverse(rx)
    ry = reverse(ry)
    ryaw = reverse(ryaw)

    dist = sum([sqrt(idx^2+idy^2) for (idx,idy) in zip(diff(rx), diff(ry))])
    # println("Final path distance is ", dist)

    return rx, ry, ryaw
end

function calc_cost(n::Node, h_rs::Array{Float64}, h_dp::Array{Float64}, ngoal::Node, c::Config)

    if length(h_rs) > 1 && length(h_dp) > 1  # Both heuristic cost are activated
        c_h_dp = h_dp[n.xind - c.minx, n.yind - c.miny]
        c_h_rs = h_rs[n.xind - c.minx, n.yind - c.miny, n.yawind - c.minyaw]
        return (n.cost + H_COST*max(c_h_dp, c_h_rs))
    elseif length(h_dp) > 1 # Distance policy heuristics is activated
        return (n.cost + H_COST*h_dp[n.xind - c.minx, n.yind - c.miny])
    elseif length(h_rs) > 1 # Reed Sheep path heuristics is activated
        return (n.cost + H_COST*h_rs[n.xind - c.minx, n.yind - c.miny, n.yawind - c.minyaw])
    end

    return (n.cost + H_COST*calc_euclid_dist(n.x[end] - ngoal.x[end],n.y[end] - ngoal.y[end], n.yaw[end] - ngoal.yaw[end]))
end

function calc_euclid_dist(x::Float64, y::Float64, yaw::Float64)
    """
    Heuristic cost function
    """
    if yaw >= pi
        yaw -= pi
    else yaw <= -pi
        yaw += pi
    end
    return sqrt(x^2+y^2+yaw^2)
end

function main()
    println(PROGRAM_FILE," start!!")

    sx = 20.0  # [m]
    sy = 20.0  # [m]
    syaw = deg2rad(90.0)
    gx = 180.0  # [m]
    gy = 100.0  # [m]
    gyaw = deg2rad(-90.0)

    ox = Float64[]
    oy = Float64[]

    for i in 0:200
        push!(ox, Float64(i))
        push!(oy, 0.0)
    end
    for i in 0:120
        push!(ox, 200.0)
        push!(oy, Float64(i))
    end
    for i in 0:200
        push!(ox, Float64(i))
        push!(oy, 120.0)
    end
    for i in 0:120
        push!(ox, 0.0)
        push!(oy, Float64(i))
    end
    for i in 0:80
        push!(ox, 40.0)
        push!(oy, Float64(i))
    end
    for i in 0:80
        push!(ox, 80.0)
        push!(oy, 120.0-Float64(i))
    end
    for i in 0:40
        push!(ox, 120.0)
        push!(oy, 120.0-Float64(i))
        push!(ox, 120.0)
        push!(oy, Float64(i))
    end
    for i in 0:80
        push!(ox, 160.0)
        push!(oy, 120.0-Float64(i))
    end

    @time rx, ry, ryaw = calc_hybrid_astar_path(sx, sy, syaw, gx, gy, gyaw, ox, oy, XY_GRID_RESOLUTION, YAW_GRID_RESOLUTION, OB_MAP_RESOLUTION)

    plot(ox, oy, ".k",label="obstacles")
    if rx != nothing
        plot(rx, ry, "-r",label="Hybrid A* path")
    end

    legend()
    grid(true)
    axis("equal")

    show()

    println(PROGRAM_FILE," Done!!")
end

if length(PROGRAM_FILE)!=0 &&
    contains(@__FILE__, PROGRAM_FILE)

    main()
end

end #module

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc


text:
using JSON, CodecZlib, Glob

const VALID_CC = ["AR", "BO", "BR", "CA", "CL", "CO", "CR", "CU", "DO", "EC", "ES", "FR", "GB", "GQ", "GT", "HN", "MX", "NI", "PA", "PE", "PR", "PY", "SV", "US", "UY", "VE"]

function parse_tweet(line)
    if length(line) == 0
        return nothing
    end

    if line[1] != '{'
        line = split(line, '\t') |> last
    end

    try
        JSON.parse(line)
    catch
        nothing
    end
end

function eachtweet(callback, filename)
    open(filename) do f
        stream = endswith(filename, ".gz") ? GzipDecompressorStream(f) : f
        for line in eachline(stream)
            tweet = parse_tweet(line)
            tweet !== nothing && callback(tweet)
        end
    end
end

function text(tweet)
    et = get(tweet, "extended_tweet", nothing)
    if et === nothing
        tweet["text"]
    else
        t = get(et, "full_text", nothing)
        t === nothing ? tweet["text"] : t
    end
end

mentions(text) = [m.match for m in eachmatch(r"(@\w+)", text)]
hashtags(text) = [m.match for m in eachmatch(r"(#\w+)", text)]
urls(text) = [m.match for m in eachmatch(r"(http.?://\S+)", text)]
language(tweet) = tweet["lang"]

function countrycode(tweet)
    place = get(tweet, "place", nothing)
    place === nothing && return nothing
    get(place, "country_code", nothing)
end

function preprocess(text)
    #text = lowercase(text)
    text = replace(text, r"[_\s]+"imx => " ")
    text = replace(text, "&" => "&", ">" => ">", "<" => "<")
    text = replace(text, r"(https?://\S+)" => " _url ")
    text = replace(text, r"#\S+" => " _htag ")
    text = replace(text, r"@\S+" => " _usr ")
    text = replace(text, r"j(a|e|i)[jaei]+"imx => s"j\1j\1")
    replace(text, r"h(a|e|i)[haei]+"imx => s"h\1h\1")
end

function emojis(text_)
    [c for c in text_ if isemoji(c)]
end

function segment_by_countrycode(tweets)
    D = Dict{String,Vector}()

    for tweet in tweets
        cc = countrycode(tweet)
        L = get(D, cc, nothing)

        if L === nothing
            D[cc] = [tweet]
        else
            push!(L, tweet)
        end
    end

    D
end

function load_and_segment_by_countrycode(getdata::Function, filenames, lang, validcc)
    D = Dict{String,Vector}()

    for filename in filenames
       println(stderr, "loading $filename")
        eachtweet(filename) do tweet
            t = text(tweet)
            cc = countrycode(tweet)
            if !startswith(t, "RT") && language(tweet) == lang && cc in validcc
                L = get(D, cc, nothing)
                data = getdata(tweet)
                if data !== nothing
                    if L === nothing
                        D[cc] = [data]
                    else
                        push!(L, data)
                    end
                end
            end
        end
    end

    D
end

load_and_segment_by_countrycode(filenames, lang, validcc) = load_and_segment_by_countrycode(identity, filenames, lang, validcc)

function update!(A, B)
    for (k, v) in B
        append!(A[k], v)
    end

    A
end

ext: jl

max_stars_repo_name: sadit/RegionalSpanish

max_issues_repo_name: sadit/RegionalSpanish

max_forks_repo_name: sadit/RegionalSpanish


text:
"""
    move(agent, model)
Move function for the agents. Alters the agents in the spatial layer.
If the agent has a Cognitive layer, moving makes them not-tired.
If the agent has a Physiological layer, moving makes them not-hungry.
"""
function move(agent, model)
    agent.pos = (rand(1:10), rand(1:10))
    if hasproperty(agent, Cognitive)
        agent.Cognitive.tired = false
    elseif hasproperty(agent, :Physiological)
        agent.Physiological.hungry = false
    end
end

"""
    think(agent, model)
If the agent has a cognitive layer, thinking alters them tired,
If the agent has a Physiological layer, thinking makes them hungry.
"""
function think(agent, model)
    if hasproperty(agent, Cognitive)
        agent.Cognitive.tired = true
    elseif hasproperty(agent, Physiological)
        agent.Physiological.hungry = true
    end
end


text:
#reads in the libaries and Data we need to run the program
using FTPClient
using DataFrames
using CSV
ModelData = CSV.read("Models.csv", silencewarnings=true)

#########Upload Test#########
#                   User Inputs
#-------------------------------------------------------
#Variables for FTP download
#Model Range{1 = Access-10, 2 = Access-13,...}
ModelRange = 1:10
# Loop in the Scenario Folders (1=rcp45 and 2=rcp85)
ScenarioRange = [1,2]
#Loop in the Variables (e.g. tasmax , tasmin , pr)
VariableRange = [1,3]
#Loop for Year
YearRange = 2021:2021
#Sets the directory path for download. Creates folder at current directory by default
Home = pwd()
#Home = "C:\\Users\\Breeze\\Documents\\Downloads\\Example_Path"
#--------------------------------------------------------

#loads the FTP to download Files
global Host = FTP(hostname="gdo-dcp.ucllnl.org",username="",password="")
cd(Host, "/pub/dcp/archive/cmip5/loca/LOCA_2016-04-02")
global BaseDir = pwd(Host)

#Model Loop
for X in ModelRange
    for Y in ScenarioRange
        #This if statement is to account for the change in directory naming depending on Variable rcp45 or rcp85
        if Y == 1
            for Z in VariableRange
                #Creates the Directory Paths for both the FTP and Local Computer
                Filepath = string(Home,"\\NETCDF\\",ModelData[X,4],"\\","16th\\",ModelData[Y,5],"\\",ModelData[X,6],"\\",ModelData[Z,3],"\\")
                ftppath = string(BaseDir,"/",ModelData[X,4],"/","16th/",ModelData[Y,5],"/",ModelData[X,6],"/",ModelData[Z,3],"/")

                #reconnects to FTP
                global Host = FTP(hostname="gdo-dcp.ucllnl.org",username="",password="")
                cd(Host, ftppath)

                #Creates Local Directory Path
                mkpath(Filepath)
                print("Downloading $(ModelData[Y,5]) $(ModelData[Z,3]) from $(ModelData[X,4])\r\n")

                for T in YearRange
                    #Creates File name to pull from FTP
                    Local_File_Name = string(ModelData[Z,3],"_day_",ModelData[X,4],"_",ModelData[Y,5],"_",ModelData[X,6],"_",T,"0101-",T,ModelData[1,7])
                    File_name = string(ModelData[Z,3],"_day_",ModelData[X,4],"_",ModelData[Y,5],"_",ModelData[X,6],"_",T,"0101-",T,ModelData[1,7])
                    End_File = string(Filepath,Local_File_Name)
                    print("Begining Download of Year $T\r\n")

                    #Checks if File is already on Local Drive, otherwise it starts the download
                    if isfile(End_File)
                        print("Download of $(ModelData[Y,5]) $(ModelData[Z,3]) from $(ModelData[X,4]) of Year $T. Already Completed.\r\n\r\n")
                    else
                        FTPClient.download(Host, File_name, End_File)
                        print("Completed Download of $(ModelData[Y,5]) $(ModelData[Z,3]) from $(ModelData[X,4]) of Year $T.\r\n\r\n")
                    end
                end
            end
        end
        if Y == 2
            for Z in VariableRange
                #Creates the Directory Paths for both the FTP and Local Computer
                Filepath = string(Home,"\\NETCDF\\",ModelData[X,4],"\\","16th\\",ModelData[Y,5],"\\",ModelData[X,8],"\\",ModelData[Z,3],"\\")
                ftppath = string(BaseDir,"/",ModelData[X,4],"/","16th/",ModelData[Y,5],"/",ModelData[X,8],"/",ModelData[Z,3],"/")

                #reconnects to FTP
                global Host = FTP(hostname="gdo-dcp.ucllnl.org",username="",password="")
                cd(Host, ftppath)

                #Creates Local Directory Path
                mkpath(Filepath)
                print("Downloading $(ModelData[Y,5]) $(ModelData[Z,3]) from $(ModelData[X,4])\r\n")

                for T in YearRange
                    #Creates File name to pull from FTP
                    Local_File_Name = string(ModelData[Z,3],"_day_",ModelData[X,4],"_",ModelData[Y,5],"_",ModelData[X,8],"_",T,"0101-",T,ModelData[1,7])
                    File_name = string(ModelData[Z,3],"_day_",ModelData[X,4],"_",ModelData[Y,5],"_",ModelData[X,8],"_",T,"0101-",T,ModelData[1,7])
                    End_File = string(Filepath,Local_File_Name)
                    print("Begining Download of Year $T\r\n")

                    #Checks if File is already on Local Drive, otherwise it starts the download
                    if isfile(End_File)
                        print("Download of $(ModelData[Y,5]) $(ModelData[Z,3]) from $(ModelData[X,4]) of Year $T. Already Completed.\r\n\r\n")
                    else
                        FTPClient.download(Host, File_name, End_File)
                        print("Completed Download of $(ModelData[Y,5]) $(ModelData[Z,3]) from $(ModelData[X,4]) of Year $T.\r\n\r\n")
                    end
                end
            end
        end
    end
end
close(Host)
cd(Home)

#print(BaseDir)


text:
using Test
using DeepShiba

function numerical_diff(f::Function, x::Real; e=10e-4)
    return (f(x + e) - f(x - e)) / 2e
end

function numerical_diff(f::Function, xs::AbstractArray; e=1e-4)
    grads = zeros(length(xs))
    for idx in 1:length(xs)
        tmp_val = xs[idx]
        xs[idx] = tmp_val + e
        fxh1 = f(xs...)
        xs[idx] = tmp_val - e
        fxh2 = f(xs...)
        grads[idx] = (fxh1 - fxh2) / 2e
        xs[idx] = tmp_val
    end
    return grads
end

function isAbout(x, y; e=10e-4)
    return y - e <= x <= y + e
end

@testset "NumericalDiffTest" begin
    f(x) = 2x^3 + x^2 + 3x + 2
    f′(x) = 6x^2 + 2x + 3
    for i in 0:0.1:1
        num_diff = numerical_diff(f, i)
        @test isAbout(num_diff, f′(i))
    end
end

@testset "OperatorTest" begin
    operators = [+, -, *, /]
    basic_func(a, b, c) = (a^2 - 4b^2 + 10) / c
    # This simple function contains all the basic operators.
    @testset "OperatorForwardTest" begin
        a, b = rand(2)
        shiba_a, shiba_b = variable.((a, b))
        for op in operators
            y = op(a, b)
            shiba_y = op(shiba_a, shiba_b)
            @test y == shiba_y.data
        end
    end

    @testset "OparatorBackwardTest" begin
        a, b = rand(2)
        shiba_a, shiba_b = variable.((a, b))
        for op in operators
            cleargrad!(shiba_a)
            cleargrad!(shiba_b)
            y = op(a, b)
            shiba_y = op(shiba_a, shiba_b)
            backward!(shiba_y)
            shiba_a_grad = shiba_a.grad.data
            shiba_b_grad = shiba_b.grad.data
            numerical_grads = numerical_diff(op, [a, b])
            @test isAbout(shiba_a_grad, numerical_grads[1])
            @test isAbout(shiba_b_grad, numerical_grads[2])
        end
    end

    @testset "BasicFunBackwardTest" begin
        a, b, c = rand(3)
        shiba_a, shiba_b, shiba_c = variable.((a, b, c))
        grads = numerical_diff(basic_func, [a, b, c])
        shiba_y = basic_func(shiba_a, shiba_b, shiba_c)
        backward!(shiba_y)
        for (shiba_grad, numerical_grad) in zip([shiba_a, shiba_b, shiba_c], [a, b, c])
            @test isAbout(numerical_grad, shiba_grad.data)
        end
    end

end

@testset "MathFuncTest" begin
    math_funcs = [sin, cos, tan, tanh]
    @testset "OneOrgMathFuncForwardTest" begin
        x = rand()
        shiba_x = variable(x)
        for func in math_funcs
            y = func(x)
            shiba_y = func(shiba_x)
            @test isAbout(y, shiba_y.data)
        end
    end
    @testset "OneOrgMathFuncBackwardTest" begin
        x = rand()
        shiba_x = variable(x)
        for func in math_funcs
            cleargrad!(shiba_x)
            shiba_y = func(shiba_x)
            backward!(shiba_y)
            numerical_grad = numerical_diff(func, x)
            @test isAbout(shiba_x.grad.data, numerical_grad)
        end
    end
end


text:
@testset "linear combination" begin

    lc = LinearCombination([1, -1])
    @test lc isa Transform
    @test cardinality(lc) == ManyToOne()

    @testset "default reduces over columns" begin
        M = [1 1; 2 2; 3 5]
        lc = LinearCombination([1, -1, 1])
        @test FeatureTransforms.apply(M, lc) == [2, 4]
        @test lc(M) == [2, 4]
    end

    @testset "dims = :" begin
        M = [1 1; 2 2; 3 5]
        lc = LinearCombination([1, -1, 1])
        @test_deprecated FeatureTransforms.apply(M, lc; dims=:)
    end

    @testset "dimension mismatch" begin
        M = [1 1 1; 2 2 2]
        lc = LinearCombination([1, -1, 1])  # there are only 2 rows
        @test_throws DimensionMismatch FeatureTransforms.apply(M, lc)
    end

    @testset "N-dim Array" begin
        A = reshape(1:27, 3, 3, 3)
        lc = LinearCombination([1, -1, 1])
        @test FeatureTransforms.apply(A, lc) == [2 11 20; 5 14 23; 8 17 26]
        @test lc(A) == [2 11 20; 5 14 23; 8 17 26]
    end
end


text:
# ---
# title: 34. Find First and Last Position of Element in Sorted Array
# id: problem34
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/find-first-and-last-position-of-element-in-sorted-array/description/>
# hidden: true
# ---
#
# Given an array of integers `nums` sorted in ascending order, find the starting
# and ending position of a given `target` value.
#
# If `target` is not found in the array, return `[-1, -1]`.
#
# **Follow up:**  Could you write an algorithm with `O(log n)` runtime
# complexity?
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [5,7,7,8,8,10], target = 8
#     Output: [3,4]
#
#
# **Example 2:**
#
#
#
#     Input: nums = [5,7,7,8,8,10], target = 6
#     Output: [-1,-1]
#
#
# **Example 3:**
#
#
#
#     Input: nums = [], target = 0
#     Output: [-1,-1]
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 105`
#   * `-109 <= nums[i] <= 109`
#   * `nums` is a non-decreasing array.
#   * `-109 <= target <= 109`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
Base.@deprecate_binding NO_FIELDS NoTangent()

const EXTERN_DEPRECATION = "`extern` is deprecated, use `unthunk` or `backing` instead, " *
    "depending on the use case."

"""
    extern(x)

Makes a best effort attempt to convert a differential into a primal value.
This is not always a well-defined operation.
For two reasons:
 - It may not be possible to determine the primal type for a given differential.
 For example, `Zero` is a valid differential for any primal.
 - The primal type might not be a vector space, thus might not be a valid differential type.
 For example, if the primal type is `DateTime`, it's not a valid differential type as two
 `DateTime` can not be added (fun fact: `Milisecond` is a differential for `DateTime`).

Where it is defined the operation of `extern` for a primal type `P` should be
`extern(x) = zero(P) + x`.

!!! note
    Because of its limitations, `extern` should only really be used for testing.
    It can be useful, if you know what you are getting out, as it recursively removes
    thunks, and otherwise makes outputs more consistent with finite differencing.

    The more useful action in general is to call `+`, or in the case of a [`Thunk`](@ref)
    to call [`unthunk`](@ref).

!!! warning
    `extern` may return an alias (not necessarily a copy) to data
    wrapped by `x`, such that mutating `extern(x)` might mutate `x` itself.
"""
@inline function extern(x)
    Base.depwarn(EXTERN_DEPRECATION, :extern)
    return x
end

extern(x::ZeroTangent) = (Base.depwarn(EXTERN_DEPRECATION, :extern); return false)  # false is a strong 0. E.g. `false * NaN = 0.0`

function extern(x::NoTangent)
    Base.depwarn(EXTERN_DEPRECATION, :extern)
    throw(ArgumentError("Derivative does not exit. Cannot be converted to an external type."))
end

extern(comp::Tangent) = (Base.depwarn(EXTERN_DEPRECATION, :extern); return backing(map(extern, comp)))  # gives a NamedTuple or Tuple

extern(x::NotImplemented) = (Base.depwarn(EXTERN_DEPRECATION, :extern); throw(NotImplementedException(x)))

@inline extern(x::AbstractThunk) = (Base.depwarn(EXTERN_DEPRECATION, :extern); return extern(unthunk(x)))

for T in (:Thunk, :InplaceableThunk)
    @eval function (x::$T)()
        Base.depwarn("`(x::" * string($T) * ")()` is deprecated, use `unthunk(x)`", Symbol(:call_, $(T)))
        return unthunk(x)
    end
end


text:
using DataStructures

"""
    abstract type AbstractStateEntry end

Any object that can be stacked into the trailer must be a subtype of this.
"""
abstract type AbstractStateEntry end

Base.show(io::IO, se::AbstractStateEntry) = write(io, "AbstractStateEntry")

"""
    Trailer

The trailer is the structure which makes it possible to memorize the previous State of our model
during the search. It makes it possible to handle the backtrack.
"""
mutable struct Trailer
    current     ::Stack{AbstractStateEntry}
    prior       ::Stack{Stack{AbstractStateEntry}}
    Trailer() = new(Stack{AbstractStateEntry}(), Stack{Stack{AbstractStateEntry}}())
end

Base.show(io::IO, tr::Trailer) = write(io, "Trailer")

"""
    StateObject{T}(value::T, trailer::Trailer)

A reversible object of value `value` that has a type `T`, storing its modification into `trailer`.
"""
mutable struct StateObject{T}
    value       ::T
    trailer     ::Trailer
end

Base.show(io::IO, so::StateObject{T}) where {T} = write(io, "StateObject{", string(T), "}: ", string(so.value))

"""
    StateEntry{T}(value::T, object::StateObject{T})

An entry that can be stacked in the trailer, containing the former `value of the object, and a reference to
the `object` so that it can be restored by the trailer.
"""
struct StateEntry{T} <: AbstractStateEntry
    value       ::T
    object      ::StateObject{T}
end

"""
    trail!(var::StateObject{T})

Store the current value of `var` into its trailer.
"""
function trail!(var::StateObject)
    push!(var.trailer.current, StateEntry(var.value, var))
end

"""
    setValue!(var::StateObject{T}, value::T) where {T}

Change the value of `var`, replacing it with `value`, and if needed, store the
former value into `var`'s trailer.
"""
function setValue!(var::StateObject{T}, value::T) where {T}
    if (value != var.value)
        trail!(var)
        var.value = value
    end
    return var.value
end

"""
    saveState!(trailer::Trailer)

Store the current state into the trailer, replacing the current stack with an empty one.
"""
function saveState!(trailer::Trailer)
    push!(trailer.prior, trailer.current)
    trailer.current = Stack{AbstractStateEntry}()
end

"""
    restoreState!(trailer::Trailer)

Iterate over the last state to restore every former value, used to backtrack every change
made after the last call to [`saveState!`](@ref).
"""
function restoreState!(trailer::Trailer)
    for se in trailer.current
        se.object.value = se.value
    end

    if isempty(trailer.prior)
        trailer.current = Stack{AbstractStateEntry}()
    else
        trailer.current = pop!(trailer.prior)
    end
end

"""
    withNewState!(func, trailer::Trailer)

Call the `func` function with a new state, restoring it after. Aimed to be used with the `do` block syntax.

# Examples
```julia
using SeaPearl
trailer = SeaPearl.Trailer()
reversibleInt = SeaPearl.StateObject{Int}(3, trailer)
SeaPearl.withNewState!(trailer) do
    SeaPearl.setValue!(reversibleInt, 5)
end
reversibleInt.value # 3
```
"""
function withNewState!(func, trailer::Trailer)
    saveState!(trailer)
    func()
    restoreState!(trailer)
end

"""
    restoreInitialState!(trailer::Trailer)

Restore every linked object to its initial state. Basically call [`restoreState!`](@ref) until not possible.
"""
function restoreInitialState!(trailer::Trailer)
    while !isempty(trailer.prior)
        restoreState!(trailer)
    end
    restoreState!(trailer)
end


text:
using Pkg, Pkg.Artifacts, SHA, Scratch, ghr_jll, Base.BinaryPlatforms

# We'll download and install `qemu-static` into an artifact for easy mounting
qemu_static_tag = "v6.0.0-2"
qemu_static_url = "https://github.com/multiarch/qemu-user-static/releases/download/$(qemu_static_tag)"
qemu_host_arch_list = (
    "x86_64",
)
qemu_target_arch_list = (
    "x86_64",
    "i386",
    "aarch64",
    "arm",
    "ppc64le",
)

@info("Creating multiarch support tarball")
artifact_hash = create_artifact() do artifact_dir
    # Download QEMU static builds
    for host_arch in qemu_host_arch_list
        for target_arch in qemu_target_arch_list
            url = "$(qemu_static_url)/$(host_arch)_qemu-$(target_arch)-static.tar.gz"
            mktempdir() do unpack_dir
                rm(unpack_dir)
                Pkg.PlatformEngines.download_verify_unpack(url, nothing, unpack_dir)
                for f in readdir(unpack_dir; join=true)
                    mv(f, joinpath(artifact_dir, basename(f)))
                end
            end
        end
    end
end

@info("Archiving")
tarball_path = joinpath(@get_scratch!("archived"), "multiarch-support.tar.gz")
archive_artifact(artifact_hash, tarball_path)

# Hash the tarball
@info("Hashing tarball")
tarball_hash = open(io -> bytes2hex(sha256(io)), tarball_path)

# Upload it to `staticfloat/Sandbox.jl`
tag_name = "multiarch-support-$(bytes2hex(artifact_hash.bytes[end-3:end]))"
@info("Uploading to staticfloat/Sandbox.jl@$(tag_name)")
run(`$(ghr_jll.ghr()) -replace $(tag_name) $(tarball_path)`)
tarball_url = "https://github.com/staticfloat/Sandbox.jl/releases/download/$(tag_name)/$(basename(tarball_path))"

# Bind it into our Artifacts.toml (advertising support for both glibc and musl)
for libc in ("glibc", "musl")
    bind_artifact!(
        joinpath(@__DIR__, "..", "Artifacts.toml"),
        "multiarch-support",
        artifact_hash;
        download_info=[(tarball_url, tarball_hash)],
        platform = Platform("x86_64", "linux"; libc="glibc"),
        force=true,
        lazy=true,
    )
end

# Eventually, we hope to be able to build this for other host architectures as well!

## Also, download and create our `multiarch-testing` artifact, which contains HelloWorldC_jll for all our architectures.
HWC_version = v"1.1.2+0"
HWC_platforms = (
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="musl"),
    Platform("i686", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="musl"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="musl"),
    Platform("armv7l", "linux"; libc="glibc"),
    Platform("armv7l", "linux"; libc="musl"),
    Platform("powerpc64le", "linux"; libc="glibc"),
    # We don't have this one yet
    #Platform("powerpc64le", "linux"; libc="musl"),
)

artifact_hash = create_artifact() do dir
    for platform in HWC_platforms
        triplet = Base.BinaryPlatforms.triplet(platform)
        url = "https://github.com/JuliaBinaryWrappers/HelloWorldC_jll.jl/releases/download/HelloWorldC-v$(HWC_version)/HelloWorldC.v$(HWC_version.major).$(HWC_version.minor).$(HWC_version.patch).$(triplet).tar.gz"
        mktempdir() do temp_dir
            rm(temp_dir)
            Pkg.PlatformEngines.download_verify_unpack(url, nothing, temp_dir)
            try
                mv(joinpath(temp_dir, "bin", "hello_world"), joinpath(dir, "hello_world.$(triplet)"))
            catch
            end
        end
    end
end

@info("Archiving")
tarball_path = joinpath(@get_scratch!("archived"), "multiarch-testing.tar.gz")
archive_artifact(artifact_hash, tarball_path)

# Hash the tarball
@info("Hashing tarball")
tarball_hash = open(io -> bytes2hex(sha256(io)), tarball_path)

# Upload it to `staticfloat/Sandbox.jl`
tag_name = "multiarch-testing-$(bytes2hex(artifact_hash.bytes[end-3:end]))"
@info("Uploading to staticfloat/Sandbox.jl@$(tag_name)")
run(`$(ghr_jll.ghr()) -replace $(tag_name) $(tarball_path)`)
tarball_url = "https://github.com/staticfloat/Sandbox.jl/releases/download/$(tag_name)/$(basename(tarball_path))"

# Bind it into our Artifacts.toml (advertising support for both glibc and musl)
bind_artifact!(
    joinpath(@__DIR__, "..", "Artifacts.toml"),
    "multiarch-testing",
    artifact_hash;
    download_info=[(tarball_url, tarball_hash)],
    force=true,
    lazy=true,
)


text:
module PGM

import Base: ==, !=
import Base: eltype, ndims, size, length, show, keys, values, parent
import Base: getindex, eachindex, start, next, done

export
    # lists
    GroupedList,
    NamedList,

    grlist_from_counts,
    grlist_by,
    ngroups,
    group,
    group_length,
    namedlist,
    namedlist_from_pairs,
    indexof,

    # variables
    VType,
    Var,
    SVar,

    realv,
    intv,
    dvar,
    rvar,
    vtype,
    cardinality,
    varlist,

    # graphs
    Edge,
    # Tree,
    UGraph,

    nvertices,
    nedges,
    nedgekinds,
    vertices,
    edges,
    edges_of,
    indegree_at,
    outdegree_at,
    degree_at,
    inedges_at,
    outedges_at

    # discrete
    # JointPMF,
    # DiscretePairwiseMRF,
    #
    # nvars,
    # vars,
    # probs,
    # marginal_i,
    # marginal,
    # conditional,
    # pwmrf,
    # has_vertex_potentials,
    # has_edge_potentials,
    # set_vertex_potentials!,
    # set_edge_potentials!,
    # vertex_potentials,
    # edge_potentials,
    # var_cardinalities,
    # edge_pdims,
    # tpotential,
    # tpotentials

# sources

include("utils.jl")
include("lists.jl")
include("variables.jl")
include("graphs.jl")

# include("discrete/jointpmf.jl")
# include("discrete/dmrf.jl")

end # module


text:
using Revise
using Test
using ZonalFlow
using BenchmarkTools

lx = 2.0π
ly = 2.0π
nx = 12
ny = 12

dom = Domain(extent=(lx,ly),res=(nx,ny))
@test dom.lx == lx
@test dom.ly == ly
@test dom.nx == nx
@test dom.ny == ny

Ω = 2.0π
θ = 0.0
μ = 0.01
ν = 0.0
ν₄ = 1.0

coeffs = Coefficients(Ω=Ω,θ=θ,μ=μ,ν=ν,ν₄=ν₄)
@test coeffs.Ω == Ω
@test coeffs.θ == θ
@test coeffs.μ == μ
@test coeffs.ν == ν
@test coeffs.ν₄ == ν₄

Ξ = 1.0
Δθ = 0.05
τ = 10.0
forcing = PointJet(Ξ=Ξ,Δθ=Δθ,τ=τ)
@test forcing.Ξ == Ξ
@test forcing.Δθ == Δθ
@test forcing.τ == τ

A₁ = 1.0
A₄ = 2.0
forcing = Kolmogorov(A₁=A₁,A₄=A₄)
@test forcing.A₁ == A₁
@test forcing.A₄ == A₄

kf = 6
dk = 1
ε = 1.0e-3
forcing = Stochastic(kf=kf,dk=dk,ε=ε)
@test forcing.kf == kf
@test forcing.dk == dk
@test forcing.ε == ε

prob1 = BetaPlane(dom,coeffs,forcing)
prob2 = BetaPlane(coeffs,forcing,extent=(lx,ly),res=(nx,ny))
@test prob1.d == dom
@test prob1.f == forcing
@test prob1.c == coeffs
@test prob2 == prob1

@benchmark zeros(NL(),dom)
zeros(NL(),dom)
@test zeros(GQL(rand(0:dom.nx-1)),dom) == zeros(NL(),dom)

@benchmark zeros(CE2(),dom)
A = zeros(GCE2(2),dom)
@test isequal(A.l,A.x[1])
@test isequal(A.h,A.x[2])

B = zeros(CE2(),dom)
@test isequal(B.l,B.x[1])
@test isequal(B.h,B.x[2])

C = zeros(GCE2(2),dom)
@test isequal(C.l,C.x[1])
@test isequal(C.h,C.x[2])


text:
module RobotSimpleFixed

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export cons_rp_robot_simple_fixed

function cons_rp_robot_simple_fixed(m)

    @variable(m, robot_simple[i=1:length(boxes)-1], Bin)
    @variable(m, robot_simple_x[i=1:length(boxes)-1], Bin)
    @variable(m, robot_simple_y[i=1:length(boxes)-1], Bin)
    @variable(m, robot_simple_z[i=1:length(boxes)-1], Bin)

    @constraint(m, c_robot_simple_x1[i=1:length(boxes)-2],
        L * (1 - robot_simple_x[i+1]) + m[:flbx_front_left_bot_x][i+1] .>=
        [m[:flbx_front_left_bot_x][j] + m[:box_x_axis_cover][j] for j in 1:i]
    )

    @constraint(m, c_robot_simple_y1[i=1:length(boxes)-2],
        W * (1 - robot_simple_y[i+1]) + m[:flby_front_left_bot_y][i+1] .>=
        [m[:flby_front_left_bot_y][j] + m[:box_y_axis_cover][j] for j in 1:i]
    )

    @constraint(m, c_robot_simple_z1[i=1:length(boxes)-2],
        H * (1 - robot_simple_z[i+1]) + m[:flbz_front_left_bot_z][i+1] .>=
        [m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j] for j in 1:i]
    )

    """
    @constraint(m, c_robot_simple_x2[i=1:length(boxes)-2],
        L * robot_simple_x[i+1] .>=
        [m[:flbx_front_left_bot_x][i+1] - m[:flbx_front_left_bot_x][j] + m[:box_x_axis_cover][j] for j in 1:i]
    )

    @constraint(m, c_robot_simple_y2[i=1:length(boxes)-2],
        W * robot_simple_y[i+1] .>=
        [m[:flby_front_left_bot_y][i+1] - m[:flby_front_left_bot_y][j] + m[:box_y_axis_cover][j] for j in 1:i]
    )

    @constraint(m, c_robot_simple_z2[i=1:length(boxes)-2],
        H * robot_simple_z[i+1] .>=
        [m[:flbz_front_left_bot_z][i+1] - m[:flbz_front_left_bot_z][i] + m[:hob_height_of_box][i] for j in 1:i]
    )
    """

    @constraint(m, c_robot_simple[i=1:length(boxes)-2],
        robot_simple_x[i+1] + robot_simple_y[i+1] + robot_simple_z[i+1] >= robot_simple[i+1]
    )

    @constraint(m, c_robot[i=1:length(boxes)-2],
        robot_simple[i+1] >= 1
    )

    return m
end

end


text:
using Documenter, VarianceComponentModels

ENV["DOCUMENTER_DEBUG"] = "true"

makedocs(
    format = Documenter.HTML(),
    sitename = "VarianceComponentModels.jl",
    modules = [VarianceComponentModels],
    pages = Any[
        "Home" => "index.md",
        "Manual" => Any[
            "man/mle_reml.md",
            "man/heritability.md"
        ],
        "API" => "man/api.md"
    ]
)

deploydocs(
    repo   = "github.com/OpenMendel/VarianceComponentModels.jl.git",
    target = "build",
#    osname = "linux",
#    julia = "1.0",
    deps = nothing,
    make = nothing)

max_stars_repo_name: Hua-Zhou/VarianceComponentModels

max_issues_repo_name: Hua-Zhou/VarianceComponentModels

max_forks_repo_name: Hua-Zhou/VarianceComponentModels


text:
__precompile__()
module TWPBVP

export twpbvpc, FInt

FInt = Int32

type TWPBVPCProblem
	fsub :: Function
	dfsub :: Function
	gsub :: Function
	dgsub :: Function
end

function __init__()
	global const libtwpbvpc = joinpath(Pkg.dir("TWPBVP"), "deps", "libtwpbvpc.so")
	global const fsub_ptr = cfunction(unsafe_fsub, Void, (Ref{FInt}, Ref{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{FInt}))

	global const dfsub_ptr = cfunction(unsafe_dfsub, Void, (Ref{FInt}, Ref{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{FInt}))

	global const gsub_ptr = cfunction(unsafe_gsub, Void, (Ref{FInt}, Ref{FInt}, Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{FInt}))

	global const dgsub_ptr = cfunction(unsafe_dgsub, Void, (Ref{FInt}, Ref{FInt}, Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{FInt}))
end

"""
Provides the conversion betweeen Julia's function (including closures)
and Fortran's subroutine. The closure is passed in `rpar` parameter
that Fortran believes to be an array of Float64, it never looks into it
anyway and only passes it around.
"""
function unsafe_fsub(n :: FInt, x :: Float64, py :: Ptr{Float64}, pdy :: Ptr{Float64}, rpar :: Ptr{Void}, ipar :: Ptr{FInt}) :: Void
	y = unsafe_wrap(Array, py, n)
	dy = unsafe_wrap(Array, pdy, n)
	problem = unsafe_pointer_to_objref(rpar) :: TWPBVPCProblem
	problem.fsub(x, y, dy)
	return nothing
end

function unsafe_dfsub(n :: FInt, x :: Float64, py :: Ptr{Float64}, pjac :: Ptr{Float64}, rpar :: Ptr{Void}, ipar :: Ptr{FInt}) :: Void
	y = unsafe_wrap(Array, py, n)
	jac = unsafe_wrap(Array, pjac, (n,n))
	problem = unsafe_pointer_to_objref(rpar) :: TWPBVPCProblem
	problem.dfsub(x, y, jac)
	return nothing
end

function unsafe_gsub(i :: FInt, n :: FInt, py :: Ptr{Float64}, pg :: Ptr{Float64}, rpar :: Ptr{Void}, ipar :: Ptr{FInt}) :: Void
	y = unsafe_wrap(Array, py, n)
	g = unsafe_wrap(Array, pg, 1)
	problem = unsafe_pointer_to_objref(rpar) :: TWPBVPCProblem
	g[1] = problem.gsub(i, y)
	return nothing
end

function unsafe_dgsub(i :: FInt, n :: FInt, py :: Ptr{Float64}, pdg :: Ptr{Float64}, rpar :: Ptr{Void}, ipar :: Ptr{FInt}) :: Void
	y = unsafe_wrap(Array, py, n)
	dg = unsafe_wrap(Array, pdg, n)
	problem = unsafe_pointer_to_objref(rpar) :: TWPBVPCProblem
	problem.dgsub(i, y, dg)
	return nothing
end

"""
Thin layer of compatibility
"""
function twpbvpc(nlbc :: FInt,
		aleft :: Float64, aright :: Float64,
		fixpnt :: Nullable{Vector{Float64}},
		ltol :: Vector{FInt}, tol :: Vector{Float64},
		linear :: Bool, givmsh :: Bool, giveu :: Bool, nmsh :: Ref{FInt},
		xx :: Vector{Float64}, u :: Array{Float64, 2}, nmax :: Ref{FInt},
		wrk :: Vector{Float64}, iwrk :: Vector{FInt},
		fsub :: Function, dfsub :: Function,
		gsub :: Function, dgsub :: Function,
		ckappa1 :: Ref{Float64}, gamma1 :: Ref{Float64},
		ckappa :: Ref{Float64},
		# rpar :: Vector{Float64},
		# ipar :: Vector{FInt},
		iflbvp :: Ref{FInt})

	# Keep problem functions in rpar
	# HACK!
	rpar = TWPBVPCProblem(fsub, dfsub, gsub, dgsub)

	# Fake external parameters
	# Can't have it 0-length as it would be Any[0] and not Float64[0]
	# local rpar :: Vector{Float64} = [0.0]
	local ipar :: Vector{FInt} = [0]

	# No need to pass these parameters
	# u is a matrix for the solution only!
	ncomp, nucol = size(u)
	nudim = ncomp
	# Get the maximum of xx
	nxxdim = length(xx)
	# max for mesh points must be the same as the number of columns of u
	assert(nucol == nxxdim)

	# Sizes of work arrays
	lwrkfl = length(wrk)
	lwrkin = length(iwrk)

	# Number of fixed mesh points
	if isnull(fixpnt)
		nfxpnt = 0
		fixpnt_v = [0.0]
	else
		fixpnt_v = fixpnt[]
		nfxpnt = length(fixpnt_v)
	end

	# Size of tolerance vector ≤ ncomp
	ntol = length(ltol)

	println("tol[1] = $(tol[1])")
	println("xx[2] = $(xx[2])")
	println("aright = $aright")

	varight = [aright]

	ccall((:twpbvpc_, libtwpbvpc), Void,
		(Ref{FInt}, Ref{FInt},					# ncomp, nlbc,
		Ref{Float64}, Ptr{Float64},					# aleft, aright
		Ref{FInt}, Ptr{Float64}, 					# nfxpnt, fixpnt
		Ref{FInt}, Ptr{FInt}, Ptr{Float64},		# ntol, ltol, tol
		Ref{FInt}, Ref{FInt}, Ref{FInt},			# linear, givmsh, giveu
		Ref{FInt}, Ref{FInt},						# nmsh, nxxdim
		Ptr{Float64}, Ref{FInt},					# xx, nudim
		Ptr{Float64}, Ref{FInt},					# u, nmax
		Ref{FInt}, Ptr{Float64},					# lwrkfl, wrk
		Ref{FInt}, Ptr{FInt},						# lwrkin, iwrk
		Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void},	# fsub, dfsub, gsub, dgsub
		Ref{Float64}, Ref{Float64},					# ckappa1, gamma1
		Ref{Float64}, Ptr{Void}, Ptr{FInt},		# ckappa, rpar, ipar
		Ref{FInt}),								# iflbvp
		ncomp, nlbc,
		aleft, & aright,
		nfxpnt, fixpnt_v,
		ntol, ltol, tol,
		linear, givmsh, giveu, nmsh,
		nxxdim, xx,
		nudim, u, nmax,					# nudim = ncomp
		lwrkfl, wrk,
		lwrkin, iwrk,
		fsub_ptr, dfsub_ptr, gsub_ptr, dgsub_ptr,
		ckappa1, gamma1, ckappa,
		pointer_from_objref(rpar), ipar,
		iflbvp)
end

end # module


text:
module PatternFormation

using StaticArrays
using LightLattices
using LoopVectorization
using FLoops
using Symbolics
using Distributions

export init_cond
export GS_Periodic!, GS_Neumann0!, NTF_Periodic!, Box!

include("InitCond.jl")
include("PDEDefs.jl")

end


text:
using Pkg, CORBITS

if !isfile(joinpath(dirname(pathof(CORBITS)),"../libcorbits.so"))
   Pkg.build("CORBITS")
end

if !Sys.iswindows()
   symlink( joinpath(dirname(pathof(CORBITS)),"../libcorbits.so"), joinpath(Pkg.devdir(),"ExoplanetsSysSim","libcorbits.so"))
else
   cp( joinpath(Pkg.devdir(),"CORBITS","libcorbits.so"), joinpath(Pkg.devdir(),"ExoplanetsSysSim","libcorbits.so"))
   cp( joinpath(dirname(pathof(CORBITS)),"../libcorbits.so"), joinpath(Pkg.devdir(),"ExoplanetsSysSim","libcorbits.so"))
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
# Test pix_get_depth()
@testset "pix_get_depth" begin
    local err = "Pix has been freed."
    local pix = pix_with()

    @test pix_get_depth(pix) == 32

    pix_delete!(pix)

    @test (@test_logs (:error, err) pix_get_depth(pix)) == 0
end

# =================================================================================================
# Test pix_get_dimensions()
@testset "pix_get_dimensions" begin
    local err = "Pix has been freed."
    local pix = pix_with()

    @test pix_get_dimensions(pix).w == 300
    @test pix_get_dimensions(pix).h == 600
    @test pix_get_dimensions(pix).d == 32

    pix_delete!(pix)

    @test (@test_logs (:error, err) pix_get_dimensions(pix)) == nothing
end


text:
# Fix EM for Gaussian SPNs... (numerically unstable)

@testset "GSPN parameter learning" begin

    println()
    @info "## Learning Gaussian SPNs"

    # First create model to generate data
    A = Vector([
        SumNode(),     # 1
        ProductNode(), # 2
        ProductNode(), # 3
        ProductNode(), # 4
        GaussianDistribution(1, 2, 18),  # 5
        GaussianDistribution(1, 11, 8), # 6
        GaussianDistribution(2, 3, 10),  # 7
        GaussianDistribution(2, -4, 7)   # 8
    ])
    B = [
        [2, 3, 4], # 1 Sum P(A, B) 0.2,0.5,0.3
        [5, 7],    # 2 Prod P(A, B)=P(A)*P(B)
        [5, 8],    # 3 Prod P(A, B)=P(A)*P(B)
        [6, 8],    # 4 Prod P(A, B)=P(A)*P(B)
    ]
    C = sparse([2,   3,     4],
               [1,   1,     1],
               [4/20, 9/20, 7/20])

    SPN = SumProductNetwork{Int64}(A, B, C)
    #println(GSPN)

    @info " Generating data set..."
    N = 10000
    data = rand(SPN,N)

    # Now learn parameters from data
    SPN2 = SumProductNetwork{Int64}(deepcopy(A), deepcopy(B), deepcopy(C))
    initialize(SPN2) # random parameter initialization

    @info "Running Expectation Maximization until convergence..."

    learner = EMParamLearner()
    learner.minimumvariance = 0.1

    while !converged(learner) && learner.steps < 1000
        SumProductNetworks.step(learner,SPN2,data)
        # mean absolute error
        error = MAE(SPN,SPN2,data)
        println("It: $(learner.steps) \t NLL: $(-learner.score) \t MAE: $error")
    end

    @test learner.score > -60000
    @test MAE(SPN,SPN2,data) < 0.3

    for node in leaves(SPN2)
        println(node)
    end

    println()

 end # End of test set


text:
# NMask
# =====
#
# Compressed bit vector of N nucleotides.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

immutable NMask
    blockmask::IndexableBitVectors.SucVector
    blocks::Vector{UInt64}
    len::Int
end

function NMask()
    return NMask(IndexableBitVectors.SucVector(), UInt64[], 0)
end

Base.copy(nmask::NMask) = NMask(copy(nmask.blockmask), copy(nmask.blocks), nmask.len)

function Base.convert(::Type{NMask}, bv::BitVector)
    n = length(bv.chunks)
    blockmask = BitVector()
    blocks = Vector{UInt64}()
    for i in 1:n
        chunk = bv.chunks[i]
        if chunk == 0
            # no N in this block
            push!(blockmask, false)
        else
            push!(blockmask, true)
            push!(blocks, chunk)
        end
    end
    return NMask(blockmask, blocks, length(bv))
end

function Base.convert(::Type{BitVector}, nmask::NMask)
    bv = BitVector()
    for i in 1:length(nmask)
        push!(bv, nmask[i])
    end
    return bv
end

Base.length(nmask::NMask) = nmask.len

divrem64(i) = i >> 6, i & 0b111111

function block_bit(i)
    d, r = divrem64(Int(i - 1))
    return d + 1, r + 1
end

function block_start(blockid)
    return (blockid - 1) << 6 + 1
end

@inline function Base.getindex(nmask::NMask, i::Integer)
    blockid, bitid = block_bit(i)
    #if !nmask.blockmask[blockid]
    if !IndexableBitVectors.unsafe_getindex(nmask.blockmask, blockid)
        return false
    end
    block = nmask.blocks[IndexableBitVectors.rank1(nmask.blockmask, blockid)]
    return ((block >> (bitid - 1)) & 1) == 1
end

function findnextn(nmask::NMask, i::Integer)
    if i > length(nmask)
        return 0
    end
    blockid, bitid = block_bit(i)
    if nmask.blockmask[blockid]
        # try to find in the current block
        block = nmask.blocks[IndexableBitVectors.rank1(nmask.blockmask, blockid)]
        d = findnext_in_block(block, bitid)
        if d > 0
            # found in the block
            return 64 * (blockid - 1) + d
        end
    end
    # search in the following blocks
    blockid = IndexableBitVectors.search1(nmask.blockmask, blockid + 1)
    if blockid == 0
        return 0
    end
    block = nmask.blocks[IndexableBitVectors.rank1(nmask.blockmask, blockid)]
    d = findnext_in_block(block, 1)
    @assert d > 0
    return 64 * (blockid - 1) + d
end

function findnext_in_block(block, bitid)
    block = block >> (bitid - 1)
    return block == 0 ? 0 : bitid + trailing_zeros(block)
end

function findprevn(nmask::NMask, i::Integer)
    if i ≤ 0
        return 0
    end
    blockid, bitid = block_bit(i)
    if nmask.blockmask[blockid]
        # try to find in the current block
        block = nmask.blocks[IndexableBitVectors.rank1(nmask.blockmask, blockid)]
        d = findprev_in_block(block, bitid)
        if d > 0
            # found in the block
            return 64 * (blockid - 1) + d
        end
    end
    # search in the following blocks
    blockid = IndexableBitVectors.rsearch1(nmask.blockmask, blockid - 1)
    if blockid == 0
        return 0
    end
    block = nmask.blocks[IndexableBitVectors.rank1(nmask.blockmask, blockid)]
    d = findprev_in_block(block, 64)
    @assert d > 0
    return 64 * (blockid - 1) + d
end

function findprev_in_block(block, bitid)
    block = block << (64 - bitid)
    return block == 0 ? 0 : bitid - leading_zeros(block)
end

function hasn_within(nmask::NMask, r::UnitRange{Int})
    # fast heuristic to check that there is no 'N' within the range
    lo, _ = block_bit(first(r))
    hi, _ = block_bit(last(r))
    if (lo     == hi && !nmask.blockmask[lo]) ||
       (lo + 1 == hi && !nmask.blockmask[lo] && !nmask.blockmask[hi])
        return false
    end
    # slower but exact algorithm to check
    return findnextn(nmask, first(r)) in r
end

function make_nbitmap(nmask::NMask, r::UnitRange{Int})
    ns = falses(length(r))
    i = findnextn(nmask, first(r))
    while i in r
        ns[i-first(r)+1] = true
        i = findnextn(nmask, i + 1)
    end
    return ns
end


text:
"""
    Base.size(layer::T) where {T <: SimpleSDMLayer}

Returns the size of the grid.
"""
Base.size(layer::T) where {T <: SimpleSDMLayer} = size(layer.grid)

"""
    Base.size(layer::T, i...) where {T <: SimpleSDMLayer}

Returns the size of the grid alongside a dimension.
"""
Base.size(layer::T, i...) where {T <: SimpleSDMLayer} = size(layer.grid, i...)

Base.length(layer::T) where {T <: SimpleSDMLayer} = count(!isnothing, layer.grid)

"""
    Base.keys(layer::T) where {T <: SimpleSDMLayer}

Returns an array of `Point` where every entry in the array is a non-`nothing`
grid coordinate.
"""
function Base.keys(layer::T) where {T <: SimpleSDMLayer}
    _lon = longitudes(layer)
    _lat = latitudes(layer)
    return [Point(_lon[p[2]], _lat[p[1]]) for p in findall(!isnothing, layer.grid)]
end


text:
struct ColorImages <: Format end

ColorImages(::Type{<:Name}; kwargs...) = ColorImages()
samplesdim(::Type{ColorImages}) = 4
saveformat(::Type{ColorImages}) = "bson"

max_stars_repo_name: VaclavMacha/DatasetProvider

max_issues_repo_name: VaclavMacha/DatasetProvider

max_forks_repo_name: VaclavMacha/DatasetProvider


text:
using BenchmarkPlottingUtilities
using Test

@testset "BenchmarkPlottingUtilities.jl" begin
    # Write your tests here.
end


text:
using Test
using ExtensibleScheduler

using TimeFrames

@testset "sample blocking scheduler" begin

    # Time as Float
    _time = time

    function print_time_noparam()
        println("From print_time_noparam $(_time())")
    end

    function print_time_args(x)
        println("From print_time_args $(_time()) $x")
    end

    function print_time_kwargs(; a="default")
        println("From print_time_kwargs $(_time()) $a")
    end

    function print_some_times()
        sched = BlockingScheduler()

        println(_time())
        #add(sched, Action(print_time_noparam), Trigger(Dates.Second(1), n=3))
        add(sched, Action(print_time_noparam), Trigger(TimeFrame("1s"), n=3))
        #add(sched, Action(print_time_args, ("positional, argument")), TimeFrame("10s"); priority=2)
        #add(sched, Action(print_time_kwargs; Dict(:a=>"keyword")...), TimeFrame("5s"); priority=1)
        #add(sched, Action(print_time_noparam), TimeFrame("10s"); priority=1)
        run(sched)
        println(_time())

        # try
        #     # start(sched)
        # catch  # KeyboardInterrupt, SystemExit
        #     #...
        # end
    end

    print_some_times()

    @test 1 == 1

end


text:
export has_recanting_witness

"""
    has_recanting_witness(g::AbstractGraph, u, v,  blocked_edges::AbstractGraph) -> Bool

In a causal DAG with edges `g`, determine whether path-specific effect from vertex `u` to `v`
with edges in `blocked_edges` can be can be computed uniquely from the data available to the investigator,
which is the case if there is no "recanting witness".
Essentially this means that `blocked_edges` could equivalently be replaced by a blocking only
outgoing edges of `u`.

See Alvin, Shpitser, Pearl (2005): "Identifiability of Path-Specific Effects",
https://ftp.cs.ucla.edu/pub/stat_ser/r321-L.pdf.
"""
function has_recanting_witness(g::AbstractGraph{T}, u::Integer, v::Integer,  blocked_edges::AbstractGraph{T}) where T
    on_blocked = zeros(Bool, nv(g))
    on_open = zeros(Bool, nv(g))
    may_blocked = zeros(Bool, nv(g))
    for e in edges(blocked_edges) # mark sources of blocked edges
        may_blocked[src(e)] = true
    end
    u == v && error("u == v")

    for nbr in outneighbors(g, u) # if there is a recanting witness, a neighbour of u is one too
        on_blocked .= false
        on_open .= false
        may_blocked[v] && has_edge(blocked_edges, v, nbr) && continue
        nbr == v && continue # direct edge

        next_open = Vector{T}()
        next_blocked = Vector{T}()

        push!(next_open, nbr)
        on_open[nbr] = true
        while !isempty(next_open) # find all reachable via direct neighbour
            src = popfirst!(next_open) # get new element from queue
            for vertex in outneighbors(g, src)
                if may_blocked[src] && has_edge(blocked_edges, src, vertex)
                    if !on_blocked[vertex]
                        on_blocked[vertex] = true
                        if vertex != v
                            push!(next_blocked, vertex)  # push onto blocked queue
                        end
                    end
                elseif !on_open[vertex]

                    on_open[vertex] = true
                    if vertex != v
                        push!(next_open, vertex) # push onto queue
                    end
                end
            end
        end
        while !isempty(next_blocked) # find all reachable through the blocked edge
            src = popfirst!(next_blocked) # get new element from queue
            for vertex in outneighbors(g, src)
                if !on_blocked[vertex]
                    on_blocked[vertex] = true
                    if vertex != v
                        push!(next_blocked, vertex)  # push onto blocked queue
                    end
                end
            end
        end
        if on_blocked[v] && on_open[v]
            return true
        end
    end
    return false
end


text:
# basically, writing stuff out in a nice format so that other things can read them

SCALE_FACTOR = 100

function output_pdb(coords::Array{Float64, 2}, out_file::AbstractString)
    # coords: each row is a coord
    n, dim = size(coords)
    out = open(out_file, "w")
    for i in 1:n
        c = coords[i, :]
        write(out, @sprintf("ATOM  "))
        write(out, @sprintf("%5d", i))
        write(out, @sprintf(" "))
        write(out, @sprintf("N   ")) # atom name
        write(out, @sprintf(" "))
        write(out, @sprintf("NOD"))
        write(out, @sprintf(" "))
        write(out, @sprintf("%c", 'A'))
        write(out, @sprintf("    "))
        write(out, @sprintf(" "))
        write(out, @sprintf("   "))
        write(out, @sprintf("%8.3f%8.3f%8.3f", (c[1]+1)*SCALE_FACTOR, (c[2]+1)*SCALE_FACTOR, (c[3]+1)*SCALE_FACTOR))
        write(out, @sprintf("%6.2f", 1.0))
        write(out, @sprintf("%6.2f", 50.0))
        write(out, @sprintf(" "))
        write(out, @sprintf("  "))
        write(out, @sprintf("\n"))
    end
    close(out)
end

function output_txt(coords::Array{Float64, 2}, out_file::AbstractString)
    # coords: each row is a coord
    # outputs a space-separated text file
    n, dim = size(coords)
    out = open(out_file, "w")
    write(out, @sprintf("%d\n", n))
    for i in 1:n
        c = coords[i,:]
        write(out, @sprintf("%8.6f ", c[1]))
        write(out, @sprintf("%8.6f ", c[2]))
        write(out, @sprintf("%8.6f ", c[3]))
        write(out, "\n")
    end
    close(out)
end


text:
# Archimedian spirals
using Luxor

spiraldata = [
  (-2, "Lituus",      130),
  (-1, "Hyperbolic", 180),
  ( 1, "Archimedes",   1),
  ( 2, "Fermat",       7)]

grid = GridRect(O - (200, 0), 130, 50)
@png begin
    background("grey15")
    for aspiral in spiraldata
        @layer begin
            sethue("antiquewhite")
            setcolor(sethue("antiquewhite"))
            setline(.5)
            Luxor.translate(nextgridpoint(grid))
            spiral(last(aspiral), first(aspiral), period=20π, :stroke)
            strokepath()
            fontsize(16)
            label(aspiral[2], :S, offset=80)
        end
    end
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
@testset "AggMode" begin
    @test typeof(LearnBase.AggMode.None()) <: LearnBase.AggregateMode
    @test typeof(LearnBase.AggMode.Sum()) <: LearnBase.AggregateMode
    @test typeof(LearnBase.AggMode.Mean()) <: LearnBase.AggregateMode
    @test typeof(LearnBase.AggMode.WeightedSum([1,2,3])) <: LearnBase.AggregateMode
    @test typeof(LearnBase.AggMode.WeightedMean([1,2,3])) <: LearnBase.AggregateMode
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
    SphericalCap{T} <: ParametricSurface{T}

Spherical cap surface, creates a half-space which is essentially the subtraction of a sphere from an infinite plane.
Only the spherical cap itself is visualized, not the plane.
The positive normal side is outside the surface.

**Can be used as a detector in [`AbstractOpticalSystem`](@ref)s.**

```julia
SphericalCap(radius::T, ϕmax::T, [surfacenormal::SVector{3,T}, centrepoint::SVector{3,T}]; interface::NullOrFresnel{T} = nullinterface(T))
```

The minimal case returns a spherical cap centered at the origin with `surfacenormal = [0, 0, 1]`.
"""
struct SphericalCap{T} <: ParametricSurface{T, 3}
    radius::T
    ϕmax::T
    zmax::T
    centrenormal::SVector{3, T}
    centrepoint::SVector{3, T}
    interface::NullOrFresnel{T}
    uvec::SVector{3, T}
    vvec::SVector{3, T}

    function SphericalCap(
        radius::T,
        ϕmax::T;
        interface::NullOrFresnel{T} = NullInterface(T),
    ) where {T <: Real}
        @assert !isnan(radius)
        @assert radius > zero(T) && zero(T) < ϕmax < T(π)
        zmax = radius * (one(T) + sin(-(π / 2 - ϕmax)))
        return new{T}(
            radius,
            ϕmax,
            zmax,
            SVector{3, T}(0.0, 0.0, 1.0),
            SVector{3, T}(0.0, 0.0, 0.0),
            interface,
            SVector{3, T}(1.0, 0.0, 0.0),
            SVector{3, T}(0.0, 1.0, 0.0),
        )
    end

    function SphericalCap(
        radius::T,
        ϕmax::T,
        centrenormal::SVector{3, T},
        centrepoint::SVector{3, T};
        interface::NullOrFresnel{T} = NullInterface(T),
    ) where {T <: Real}
        @assert radius > zero(T) && zero(T) < ϕmax < T(π)
        n̂ = normalize(centrenormal)
        rotationvec = SVector{3, T}(0.0, 1.0, 0.0)
        if abs(dot(rotationvec, n̂)) == one(T)
            rotationvec = SVector{3, T}(1.0, 0.0, 0.0)
        end
        uvec = normalize(cross(normalize(rotationvec), n̂))
        vvec = normalize(cross(n̂, uvec))
        zmax = radius * (one(T) + sin(-(π / 2 - ϕmax)))
        return new{T}(
            radius,
            ϕmax,
            zmax,
            n̂,
            centrepoint,
            interface,
            uvec,
            vvec,
        )
    end
end

function Base.show(io::IO, a::SphericalCap{T}) where {T <: Real}
    return print(
        io,
        "SphericalCap{$T}($(a.radius), $(a.centrepoint), $(a.centrenormal), $(a.ϕmax), $(interface(a)))",
    )
end

interface(a::SphericalCap{T}) where {T <: Real} = a.interface
radius(a::SphericalCap{T}) where {T <: Real} = a.radius
centroid(r::SphericalCap{T}) where {T <: Real} = r.centrepoint

function uvrange(::Type{SphericalCap{T}}) where {T <: Real}
    return ((-T(π), T(π)), (zero(T), one(T)))
end # θ and ρ

function normal(r::SphericalCap{T}, θ::T, ρ::T) where {T <: Real}
    ϕ = π / 2 - ρ * r.ϕmax
    return -normalize(
        cos(-ϕ) * (cos(-θ) * r.uvec + sin(-θ) * r.vvec) +
        sin(-ϕ) * r.centrenormal,
    )
end

function point(r::SphericalCap{T}, θ::T, ρ::T) where {T <: Real}
    ϕ = π / 2 - ρ * r.ϕmax
    return centroid(r) +
           radius(r) * (
        cos(-ϕ) * (cos(-θ) * r.uvec + sin(-θ) * r.vvec) +
        sin(-ϕ) * r.centrenormal +
        r.centrenormal
    )
end

function uv(r::SphericalCap{T}, x::T, y::T, z::T) where {T <: Real}
    return uv(r, SVector{3, T}(x, y, z))
end
function uv(r::SphericalCap{T}, p::SVector{3, T}) where {T <: Real}
    prel = p - centroid(r) - radius(r) * r.centrenormal
    v = dot(prel, r.vvec) / radius(r)
    u = dot(prel, r.uvec) / radius(r)
    θ = -NaNsafeatan(v, u)
    z = clamp(dot(prel, -r.centrenormal) / radius(r), -one(T), one(T))
    ϕ = NaNsafeasin(z)
    ρ = (π / 2 - ϕ) / r.ϕmax
    return θ, ρ
end

function inside(r::SphericalCap{T}, p::SVector{3, T}) where {T <: Real}
    prel = (p - centroid(r) - radius(r) * r.centrenormal) / radius(r)
    v = dot(prel, r.vvec)
    u = dot(prel, r.uvec)
    if (u == zero(T) && v == zero(T)) || sqrt(u^2 + v^2) < sin(r.ϕmax)
        l = norm(prel)
        if l < one(T)
            return false
        end
    end
    return dot(r.centrenormal, p - centroid(r)) < r.zmax
end

function onsurface(r::SphericalCap{T}, p::SVector{3, T}) where {T <: Real}
    prel = (p - centroid(r) - radius(r) * r.centrenormal) / radius(r)
    v = dot(prel, r.vvec)
    u = dot(prel, r.uvec)
    if (u == zero(T) && v == zero(T)) || sqrt(u^2 + v^2) < sin(r.ϕmax)
        return norm(prel) == one(T) &&
               dot(r.centrenormal, p - centroid(r)) < r.zmax
    else
        false
    end
end

function uvtopix(
    ::SphericalCap{T},
    uv::SVector{2, T},
    imsize::Tuple{Int, Int},
) where {T <: Real}
    θ, ρ = uv
    h, w = imsize
    u = (cos(θ) * ρ + one(T)) / 2
    v = (sin(θ) * ρ + one(T)) / 2
    pixu = Int(floor((w - 1) * u)) + 1
    pixv = Int(floor((h - 1) * v)) + 1
    return pixu, pixv
end

function surfaceintersection(
    sph::SphericalCap{T},
    r::AbstractRay{T, 3},
) where {T <: Real}
    rad = radius(sph)
    orel = origin(r) - centroid(sph) - rad * sph.centrenormal

    ox, oy, oz = orel
    dx, dy, dz = direction(r)

    a = dx^2 + dy^2 + dz^2
    b = 2 * (ox * dx + oy * dy + oz * dz)
    c = (ox^2 + oy^2 + oz^2) - rad^2

    temp = quadraticroots(a, b, c)

    if temp === nothing
        if inside(sph, origin(r))
            return rayorigininterval(Infinity(T))
        else
            return EmptyInterval(T)
        end
    end

    t1, t2 = temp

    if isapprox(t1, t2, rtol = 1e-12, atol = 2 * eps(T))
        if inside(sph, origin(r))
            return rayorigininterval(Infinity(T))
        else
            return EmptyInterval(T)
        end
    end

    if t1 > zero(T)
        pt1 = point(r, t1)
    else
        pt1 = nothing
    end

    if t2 > zero(T)
        pt2 = point(r, t2)
    else
        pt2 = nothing
    end

    let int1 = nothing, int2 = nothing
        if pt1 !== nothing
            θ, ρ = uv(sph, pt1)
            if zero(T) <= ρ <= one(T)
                int1 = Intersection(
                    t1,
                    pt1,
                    normal(sph, θ, ρ),
                    θ,
                    ρ,
                    interface(sph),
                )
            end
        end

        if pt2 !== nothing
            θ, ρ = uv(sph, pt2)
            if zero(T) <= ρ <= one(T)
                int2 = Intersection(
                    t2,
                    pt2,
                    normal(sph, θ, ρ),
                    θ,
                    ρ,
                    interface(sph),
                )
            end
        end

        if int1 !== nothing && int2 !== nothing
            if t1 <= t2
                return Interval(int1, int2)
            else
                return Interval(int2, int1)
            end
        elseif int1 !== nothing
            if inside(sph, origin(r))
                return rayorigininterval(int1)
            else
                return positivehalfspace(int1)
            end
        elseif int2 !== nothing
            if inside(sph, origin(r))
                return rayorigininterval(int2)
            else
                return positivehalfspace(int2)
            end
        else
            if inside(sph, origin(r))
                return rayorigininterval(Infinity(T))
            else
                return EmptyInterval(T)
            end
        end
    end
end

function BoundingBox(a::SphericalCap{T}) where {T <: Real}
    return BoundingBox(
        Plane(a.centrenormal, a.centrepoint + a.radius * a.centrenormal),
    )
end


text:
"""
    IRDumps.tofiles(outputdir, dumpfile; clean = true)

Split `dumpfile` that contains the stderr of `JULIA_LLVM_ARGS='-print-after-all'
julia`) into multiple files in a directory `outputdir`.  It empties `outputdir`
first if `clean = true` (default).
"""
function IRDumps.tofiles(outputdir, dumpfile; clean = true)
    @argcheck isfile(dumpfile)
    if clean
        @argcheck !isfile(outputdir)
        if isdir(outputdir)
            rm(outputdir, recursive = true)
        end
    end
    mkpath(outputdir)
    open(dumpfile) do input
        counter = Ref(0)
        foreachdump(input) do foreachline
            tmppath = joinpath(outputdir, "_tmp.ll")
            name = open(tmppath, write = true) do output
                foreachline() do ln
                    print(output, ln)
                end
            end
            filename = string(counter[]; pad = 5) * "-" * asfilename(name) * ".ll"
            mv(tmppath, joinpath(outputdir, filename))
            counter[] += 1
        end
    end
end


text:
function read_volume(
    io::TranscodingStream{GzipDecompressor},
    ::Type{T},
    dim::Tuple{Vararg{Any,N}},
    map::Bool,
) where {T,N}

    if T === Bool
        A = BitArray{N}
    else
        A = Array{T,N}
    end
    if map
        close(io)
        error("cannot mmap a gzipped NIfTI file")
    else
        return read!(io, A(undef, dim))
    end
end

function read_volume(io, ::Type{T}, dim::Tuple{Vararg{Any,N}}, map::Bool) where {T,N}
    if T === Bool
        A = BitArray{N}
    else
        A = Array{T,N}
    end
    if map
        return Mmap.mmap(io, A, dim)
    else
        return read!(io, A(undef, dim))
    end
end

niopen(file::AbstractString, mode) = niopen(open(file, mode))
function niopen(io)
    if isgz(io)
        return GzipDecompressorStream(io)
    else
        return io
    end
end