text:
module UnitfulAngles

import Unitful
import Dates
using Unitful: @unit, Quantity, NoDims, @u_str, uconvert, ustrip
import Base: sin, cos, tan, sec, csc, cot, asin, acos, atan, asec, acsc, acot, convert

######################### Angle units ##########################################
@unit turn          "τ"             Turn          2π*u"rad"     false
@unit doubleTurn    "§"             DoubleTurn    2turn         false
@unit halfTurn      "π"             HalfTurn      turn//2       false
@unit quadrant      "⦜"             Quadrant      turn//4       false
@unit sextant       "sextant"       Sextant       turn//6       false
@unit octant        "octant"        Octant        turn//8       false
@unit clockPosition "clockPosition" ClockPosition turn//12      false
@unit hourAngle     "hourAngle"     HourAngle     turn//24      false
@unit compassPoint  "compassPoint"  CompassPoint  turn//32      false
@unit hexacontade   "hexacontade"   Hexacontade   turn//60      false
@unit brad          "brad"          BinaryRadian  turn//256     false
@unit diameterPart  "diameterPart"  DiameterPart  1u"rad"/60    false # ≈ turn/377
@unit grad          "ᵍ"             Gradian       turn//400     false
@unit arcminute     "′"             Arcminute     u"°"//60      false # = turn/21,600
@unit arcsecond     "″"             Arcsecond     u"°"//3600    false # = turn/1,296,000

######################### Functions ############################################

# cos and sin have *pi versions, and *d versions
for _f in (:cos, :sin)
    @eval $_f(x::Quantity{T, NoDims, typeof(halfTurn)}) where {T} = $(Symbol("$(_f)pi"))(ustrip(x))
    @eval $_f(x::Quantity{T, NoDims, typeof(diameterPart)}) where {T} = $_f(ustrip(uconvert(u"rad", x)))
    for _u in (doubleTurn, turn, quadrant, sextant, octant, clockPosition, hourAngle, compassPoint, hexacontade, brad, grad, arcminute, arcsecond)
        @eval $_f(x::Quantity{T, NoDims, typeof($_u)}) where {T} = $(Symbol("$(_f)pi"))(ustrip(uconvert(halfTurn, x)))
    end
end

# These functions don't have *pi versions, but have *d versions
for _f in (:tan, :sec, :csc, :cot)
    @eval $_f(x::Quantity{T, NoDims, typeof(diameterPart)}) where {T} = $_f(ustrip(uconvert(u"rad", x)))
    for _u in (doubleTurn, turn, halfTurn, quadrant, sextant, octant, clockPosition, hourAngle, compassPoint, hexacontade, brad, grad, arcminute, arcsecond)
        @eval $_f(x::Quantity{T, NoDims, typeof($_u)}) where {T} = $(Symbol("$(_f)d"))(ustrip(uconvert(u"°", x)))
    end
end

# Inverse functions
for _f in (:acos, :asin, :atan, :asec, :acsc, :acot),
    _u in (diameterPart, u"°", u"rad", doubleTurn, turn, halfTurn, quadrant, sextant, octant, clockPosition, hourAngle, compassPoint, hexacontade, brad, grad, arcminute, arcsecond)
    @eval $_f(::typeof($_u), x::Number) = uconvert($_u, $_f(x)*u"rad")
end

for _u in (diameterPart, u"°", u"rad", doubleTurn, turn, halfTurn, quadrant, sextant, octant, clockPosition, hourAngle, compassPoint, hexacontade, brad, grad, arcminute, arcsecond)
    @eval atan(::typeof($_u), y::Number, x::Number) = uconvert($_u, atan(y, x)*u"rad")
end

_tohour(x::T) where {T <: Dates.TimePeriod} = x/convert(T, Dates.Hour(1)) # credit to kristoffer.carlsson see https://discourse.julialang.org/t/convert-time-interval-to-seconds/3806/2

# Fun conversion between time and angles
# NOTE: not sure whether to use `convert` or `uconvert`
for _u in (diameterPart, u"°", u"rad", doubleTurn, turn, halfTurn, quadrant, sextant, octant, clockPosition, hourAngle, compassPoint, hexacontade, brad, grad, arcminute, arcsecond)
    @eval begin
        function convert(::typeof($_u), t::Dates.Time)
            x = t - Dates.Time(0, 0, 0)
            uconvert($_u, _tohour(x)*hourAngle)
        end
    end
    @eval convert(::Type{Dates.Time}, x::Quantity{T, NoDims, typeof($_u)}) where {T} = Dates.Time(0, 0, 0) + Dates.Nanosecond(round(Int, ustrip(uconvert(hourAngle, x))*3600000000000))
end

# Enable precompilation with Unitful extended units
# http://ajkeller34.github.io/Unitful.jl/stable/extending/#precompilation
const localunits = Unitful.basefactors
function __init__()
    merge!(Unitful.basefactors, localunits)
    Unitful.register(UnitfulAngles)
end

end # module


text:
import Clustering
import Distances
import LinearAlgebra: dot

"""
    noselection(algo, data)

Perform no clustering selection.
"""
noselection(algo::Function, data::AbstractMatrix{<:Real}) = Clustering.assignments(algo(data))

"""
Knee point detection
"""
function knee(x)
    n = length(x)
    pts = hcat(1:n, x)'
    # calculate line vector
    lvec = pts[:,end] - pts[:,1]
    lvec ./= sqrt(sum(lvec.^2))
    # find a vector parallel to the line vector
    vecfst = pts .- pts[:,1]
    vecfstpar = lvec * mapslices(v->dot(v, lvec), vecfst, dims=1)
    # calculate a length of a projection to the line vector
    ldist = vec(sqrt.(sum((vecfst - vecfstpar).^2, dims=1)))
    # get maximal projection
    return findmax(ldist)
end

function getmaxk(N::Int; kwargs...)
    N <= 5 && return 2

    # setup parameters
    Kmax = 10 # default
    for (p,v) in kwargs
        if p == :maxk
            Kmax = v
        end
    end
    return min(floor(Int, N/2), Kmax)
end

"""
    elbow(algo, data)

Perform automatic selection of the best clustering of `data` by `algo` using elbow method.
Clustering algorithm `algo` must be a function with one argument which determines number of clusters.
"""
function elbow(algo::Function, data::AbstractMatrix{<:Real}; kwargs...)
    N = size(data,2)
    Kmax = getmaxk(N; kwargs...)

    S = [let C = algo(data, k)
            # println("$k => ", C.centers)
            # println("$k => ", C.costs)
            C.totalcost => Clustering.assignments(C)
         end
    for k in 2:Kmax]
    pushfirst!(S, sum(mapslices(p->sum(p.^2), data .- mean(data, dims=2), dims=1))=>fill(1,N))

    # return a clustering corespondig to an elbow point
    return S[knee(map(first, S)) |> last] |> last
end

"""
    silhouette(algo, data)

Perform automatic selection of the best clustering of `data` by `algo` using silhouette method.
Clustering algorithm `algo` must be a function with one argument which determines number of clusters.
"""
function silhouette(algo::Function, data::AbstractMatrix{<:Real}; kwargs...)
    N = size(data,2)
    Kmax = getmaxk(N; kwargs...)
    S = [let C = algo(data, k),
             D = Distances.pairwise(Distances.Euclidean(), data)
             mean(Clustering.silhouettes(C, D)) => Clustering.assignments(C)
    end
    for k in 2:Kmax]

    # return a clustering corespondig to a maximum average silhouette
    return S[findmax(map(first, S)) |> last] |> last
end

ext: jl


text:
"""
# module Parser

- Julia version:
- Author: bramb
- Date: 2019-02-13

# Examples

```jldoctest
julia>
```
"""
module Parser

export my_parse

function my_parse()
    return "Hello World"
end

end

max_stars_repo_name: rhdekker/HyperGraphTools

max_issues_repo_name: rhdekker/HyperGraphTools

max_forks_repo_name: rhdekker/HyperGraphTools


text:
abstract type Shape end

abstract type ShapeWithParent <: Shape end

macro shape_getter(field, owner_type)
    func_sym = Symbol("shape_get_", field)
    return esc(quote
        $func_sym(s::ShapeWithParent) = $func_sym(s.parent)
        $func_sym(s::$owner_type) = s.$field
    end)
end

struct Rectangle <: Shape
    corner_upper_left  :: Tuple{Float64, Float64}
    corner_lower_right :: Tuple{Float64, Float64}
end
@shape_getter(corner_upper_left, Rectangle)
@shape_getter(corner_lower_right, Rectangle)

struct Circle <: Shape
    center      :: Tuple{Float64, Float64}
    radius      :: Float64
end
@shape_getter(center, Circle)
@shape_getter(radius, Circle)

struct GrayShape{shapeType <: Shape} <: ShapeWithParent
    parent      :: shapeType
    gray_value  :: Float64
end
@shape_getter(gray_value, GrayShape)

struct AlphaShape{shapeType <: Shape} <: ShapeWithParent
    parent      :: shapeType
    alpha_value :: Float64
end
@shape_getter(alhpa_value, AlphaShape)

r1 = AlphaShape(GrayShape(Rectangle((0.0, 1.0), (1.0, 0.0)), 0.5), 0.9)
r2 = GrayShape(AlphaShape(Rectangle((0.0, 1.0), (1.0, 0.0)), 0.9), 0.5)

println("r1: gray_value: ", shape_get_gray_value(r1))
println("r2: gray_value: ", shape_get_gray_value(r2))

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
using Javis
using Unitful
using Plots
import ColorSchemes.deep

function makeRender(iterations::Int, particles::Vector, positions::Array{Quantity,3}; bgcolour="white", lenUnit::Unitful.LengthFreeUnits=u"m", outputFile="circle.gif", size::Int=500)
    print("Beginning animated render.\n")
    @warn "This is extremely buggy. Use of the alternative render methods is suggested."

    # Init video
    function ground(args...)
        background(bgcolour) # canvas background
        sethue("black") # pen color
    end

    function object(p=O, colour="black")
        sethue(colour)
        circle(p, 5, :fill)
        return p
    end

    scaleFactor = size / ustrip(lenUnit, maximum(positions[:,:,:2]))

    myvideo = Video(size, size)
    Background(1:iterations, ground)

    print("Render initialisation complete. Loading points.\n")

    # Get paths followed by particles [iterations, particles]
    point_positions = Matrix{Point}(undef, iterations, length(particles))
    for j in 1:length(particles)
        for i in 1:iterations
            p = ustrip.(lenUnit, positions[i, j, :]) * scaleFactor
            point_positions[i, j] = Point(p[1], p[2])
        end
    end

    # Add particles to vide
    dots = Vector{Object}(undef, length(particles))
    for i in 1:length(particles)
        p = particles[i]
        colour = deep[ustrip(u"C", p.q)]
        pos = ustrip.(lenUnit, p.r) * scaleFactor
        dots[i] = Object((args...) -> object(O, colour), Point(pos[1], pos[2]))
        act!(dots[i], Action(1:iterations, follow_path(point_positions[:, i])))
    end

    print("Points added to video. Beginning render.\n")

    # Render output
    render(
        myvideo;
        pathname=outputFile
    )

    print("Render complete. File saved as $(outputFile).\n")
end

function plotTrajectories(positions)
    positions = ustrip.(u"m", positions)
    plt = plot(positions[:, :, 1], positions[:, :, 2], positions[:, :, 3], legend=false,
    xlabel="x (m)", ylabel="y (m)", zlabel="z(m)",
    title="Particle trajectories")
    display(plt)
end

function plotAtTime(positions, timeIndex, dt)
    positions = ustrip.(u"m", positions)
    plt = scatter(positions[timeIndex, :, 1]', positions[timeIndex, :, 2]', positions[timeIndex, :, 3]',
    legend=false, xlabel="x (m)", ylabel="y (m)", zlabel="z(m)", title="Positions at time $(timeIndex*dt)")
    display(plt)
end

function plotParticle(positions, particleIndex)
    positions = ustrip.(u"m", positions)
    plt = plot(positions[:,particleIndex,1], positions[:,particleIndex,2], positions[:,particleIndex,3],
    legend=false, xlabel="x (m)", ylabel="y (m)", zlabel="z(m)", title="Trajectory of particle $particleIndex")
    display(plt)
end


text:
# ======================
# Setup
# ======================

# The benchmarks in this file are taken from
# [Project Daisy](https://github.com/malyzajko/daisy/blob/master/testcases/).
SUITE["Daisy"] = BenchmarkGroup()
RESULTS = SUITE["Daisy"]

# The following domains are used throughout the tests in Tables 3-5 in
# [1] Althoff, M., Grebenyuk, D., & Kochdumper, N. (2018). Implementation of Taylor models in CORA 2018.
#     In Proc. of the 5th International Workshop on Applied Verification for Continuous and Hybrid Systems.
a = Interval(-4.5, -0.3)
b = Interval(0.4, 0.9)
c = Interval(3.8, 7.8)
d = Interval(8.0, 10.0)
e = Interval(-10.0, 8.0)
f = Interval(1.0, 2.0)

# Dictionary to hold the vector of relative precision intervals for each benchmark
RELPREC = Dict{String, Any}()

# Dictionaries to store data and reference values
DAISY = Dict{String, Any}()

# load functions
include("univariate.jl")
include("multivariate.jl")

# ===============================
# Run benchmarks
# ===============================
for name in keys(DAISY)
    RESULTS[name] = BenchmarkGroup()
    RELPREC[name] = Dict()
    func = DAISY[name]["f"]
    dom = DAISY[name]["dom"]
    ref = DAISY[name]["ref"]
    for ord in [2, 5, 10]
        RESULTS[name]["order $ord"] = BenchmarkGroup()
        RELPREC[name]["order $ord"] = Dict()

        RESULTS[name]["order $ord"]["Taylor Model subs."] = @benchmarkable enclose($func, $dom, :TaylorModels, order=$ord, normalize=false)
        RELPREC[name]["order $ord"]["Taylor Model subs."] = relative_precision(enclose(func, dom, :TaylorModels, order=ord, normalize=false), ref)

        RESULTS[name]["order $ord"]["Normalized Taylor Model subs."] = @benchmarkable enclose($func, $dom, :TaylorModels, order=$ord, normalize=true)
        RELPREC[name]["order $ord"]["Normalized Taylor Model subs."] = relative_precision(enclose(func, dom, :TaylorModels, order=ord, normalize=true), ref)

        RESULTS[name]["order $ord"]["Interval Arithmetic subs."] = @benchmarkable enclose($func, $dom, :IntervalArithmetic)
        RELPREC[name]["order $ord"]["Interval Arithmetic subs."] = relative_precision(enclose(func, dom, :IntervalArithmetic), ref)
    end
end


text:
using CMPFit, GFit, Gnuplot
using QSFit, DataStructures, Statistics, Dierckx

import QSFit: known_spectral_lines, add_qso_continuum!, LineComponent

export T2Recipe

abstract type T2Recipe <: DefaultRecipe end

function default_options(::Type{T}) where T <: T2Recipe
    out = OrderedDict{Symbol, Any}()
    out[:host_template] = "Ell5"
    out[:use_host_template] = true
    out[:use_balmer] = false
    out[:use_ironuv] = false
    out[:use_ironopt] = false
    out[:use_lorentzian_profiles] = false
    out[:n_unk] = 2
    out[:unk_avoid] = [4863 .+ [-1,1] .* 50,
                       6565 .+ [-1,1] .* 150,
                       5008 .+ [-1,1] .* 25]
    return out
end

function known_spectral_lines(source::QSO{T}) where T <: T2Recipe
    list = [
        NarrowLine(                      :Lyb                       ),
        NarrowLine(                      :Lya                       ),
        NarrowLine(                      :NV_1241                   ),
        NarrowLine(                      :CIV_1549                  ),
        NarrowLine(                      :CIII_1909                 ),
        NarrowLine(                      :MgII_2798                 ),
        NarrowLine(                      :NeV_3426                  ),
        NarrowLine(                      :OII_3727                  ),
        NarrowLine(                      :NeIII_3869                ),
        NarrowLine(                      :Hg                        ),
        NarrowLine(                      :Hb                        ),
        NarrowLine(                      :OIII_4959                 ),
        NarrowLine(                      :OIII_5007                 ),
        AsymmTailLine(                   :OIII_5007      ,     :blue),
        NarrowLine(                      :OI_6300                   ),
        NarrowLine(                      :OI_6364                   ),
        NarrowLine(                      :NII_6549                  ),
        NarrowLine(                      :Ha                        ),
        NarrowLine(                      :NII_6583                  ),
        NarrowLine(                      :SII_6716                  ),
        NarrowLine(                      :SII_6731                  )]
    return list
end

function add_qso_continuum(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: T2Recipe
    λ = domain(model)[:]

    comp = QSFit.powerlaw(median(λ))
    comp.alpha.val = -1.8

    model[:qso_cont] = comp
    push!(model[:Continuum].list, :qso_cont)
    evaluate!(model)
end

function LineComponent(source::QSO{T}, line::NarrowLine, multicomp::Bool) where T <: T2Recipe
    lc = LineComponent(parent_recipe(source), line, multicomp) # invoke parent recipe
    lc.comp.fwhm.low  = 10
    lc.comp.fwhm.high = 1000
    lc.comp.voff.high = 500
    @info line.tid lc.comp.fwhm.val
    return lc
end

function LineComponent(source::QSO{T}, line::AsymmTailLine, multicomp::Bool) where T <: T2Recipe
    comp = LineComponent(source, GenericLine(line.tid), multicomp).comp
    comp.fwhm.val  = 500
    comp.fwhm.high = 1e3
    comp.voff.low  = line.side == :blue ?    0  :  720
    comp.voff.high = line.side == :blue ?  720  :    0
    return LineComponent(line, comp, multicomp)
end

function default_unk_line(source::QSO{T}) where T <: T2Recipe
    comp = SpecLineGauss(5e3)
    comp.norm.val = 0.
    comp.center.fixed = false
    comp.center.low = 0
    comp.center.high = Inf
    comp.fwhm.val  = 500
    comp.fwhm.low  = 10
    comp.fwhm.high = 2e3
    comp.voff.fixed = true
    return comp
end

function add_patch_functs!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: T2Recipe
    # Patch parameters
    @try_patch! begin
        # model[:OIII_4959].norm = model[:OIII_5007].norm / 3
        model[:OIII_4959].voff = model[:OIII_5007].voff
    end
    @try_patch! begin
        model[:OIII_5007_bw].voff += model[:OIII_5007].voff
        model[:OIII_5007_bw].fwhm += model[:OIII_5007].fwhm
        model[:OIII_5007_bw].norm += model[:OIII_5007].norm # is this correct way round?
    end
    @try_patch! begin
        # model[:OI_6300].norm = model[:OI_6364].norm / 3
        model[:OI_6300].voff = model[:OI_6364].voff
    end
    @try_patch! begin
        # model[:NII_6549].norm = model[:NII_6583].norm / 3
        model[:NII_6549].voff = model[:NII_6583].voff
    end
    @try_patch! begin
        # model[:SII_6716].norm = model[:SII_6731].norm / 1.5
        model[:SII_6716].voff = model[:SII_6731].voff
    end
end


text:
abstract type SQLQuery <: SQLStatement end

export  Select,
        addqueryfield!,
        setwhereclause!,
        setgroupby!,
        sethavingclause!,
        select

"""
```
Select(
    name::Union{Var, AbstractString, JoinTable}="",
    alias::Union{Var, AbstractString}="";
    distinct::Bool=false )
Select(
    sq::Select,
    alias::Union{Var, AbstractString};
    distinct::Bool=false )
```
This mutable struct represents a `SELECT` SQL statement where `name` is the table to select from, `alias` is the alias for the table, if any is given, and `distinct` determines whether the `DISTINCT` flag is set, that is, identical records are suppressed in the output.

With the second version of the constructor that takes a `Select` as first argument, `alias` is an obligatory argument.
"""
mutable struct Select <: SQLQuery
    isdistinct::Bool
    name::String
    alias::String
    columns::Vector{String}
    aliases::Vector{Var}
    wherec::SQLCondition
    groupby::Vector{Var}
    having::SQLCondition
    orderby::Vector{Var}
    asc::Vector{Bool}
    nulllast::Vector{Bool}
    limit::Int
    offset::Int

    Select( name::SVar="", alias::SVar=""; distinct::Bool=false ) = new( distinct, name, alias |> string, String[], Var[], nocond, Var[], nocond, String[], Bool[], Bool[], -1, 0 )
end

Select( jc::JoinTable, alias::SVar=""; distinct::Bool=false ) = Select( "($jc)", alias, distinct=distinct )
Select( sq::Select, alias::SVar; distinct::Bool=false ) = Select( "($sq)", alias, distinct=distinct )

function Base.show( io::IO, sq::Select )
    isempty(sq.columns) && return
    print( io, "SELECT " )
    sq.isdistinct && print( io, "DISTINCT " )
    colstr = map( eachindex(sq.columns) ) do ii
        isempty(sq.aliases[ii]) && return string(sq.columns[ii])
        "$(sq.columns[ii]) AS $(sq.aliases[ii])"
    end
    print( io, join( colstr, ", " ) )

    if !isempty(sq.name)
        print( io, "\n    FROM $(sq.name)" )
        isempty(sq.alias) || print( io, " AS \"$(sq.alias)\"" )
    end

    isempty(sq.wherec) || print( io, "\n    WHERE $(sq.wherec)" )
    isempty(sq.groupby) || print( io, "\n    GROUP BY $(join( sq.groupby, ", " ))" )
    isempty(sq.having) || print( io, "\n    HAVING $(sq.having)" )

    if !isempty(sq.orderby)
        print( io, "\n    ORDER BY " )
        ostr = map( eachindex(sq.orderby) ) do ii
            string( sq.orderby[ii], ' ', sq.asc[ii] ? "ASC" : "DESC", " NULLS ", sq.nulllast[ii] ? "LAST" : "FIRST" )
        end
        print( io, join( ostr, ", " ) )
    end

    if (sq.limit >= 0) || (sq.offset > 0)
        print( io, "\n   " )
        sq.limit < 0 || print( io, " LIMIT ", sq.limit )
        sq.offset > 0 && print( io, " OFFSET ", sq.offset )
    end
end

"""
```
addqueryfield!(
    column::SVar,
    alias::SVar="" )
```
This function permits the user to add a command that adds the column `column` to the result of a `Struct` struct. This column can have the alias `alias`.

The output of this function is a function that a `Select` as argument and returns the updated argument, such that Julia's function chaining syntax can be used.

Example:
```
sq |> addqueryfield!( "var1" ) |> addqueryfield!( "var2 + var3", "varcomb" )
```
sets the commands to add the columns `var1` and `var2 + var3`, aliased by `varcomb`, to the output of the `Select` statement.
"""
addqueryfield!( column::SVar, alias::SVar="" ) = sq::Select -> begin
    push!( sq.columns, column |> string )
    push!( sq.aliases, alias |> Var )
    sq
end

"""
```
setwhereclause!( cond::SQLCondition )
```
This function permits the user to add a command that sets the `WHERE` clause of a `Update`, `DeleteFrom`, or `Select` struct to `cond`.

The output of this function is a function that takes an argument of the aforementioned types and returns the updated argument, such that Julia's function chaining syntax can be used.

Example:
```
sq |> setwhereclause!( cond )
```
sets the `WHERE` clause of `sq` to `cond` where `cond` ha sbeen defined previously.
"""
function setwhereclause!( cond::SQLCondition )
    dm::Union{Update, DeleteFrom, Select} -> begin
        dm.wherec = cond
        dm
    end
end

"""
```
setgroupby!( columns::Vector{Var} )
setgroupby!( columns::Vector{T} ) where T <: AbstractString
setgroupby!( columns::Union{Var, AbstractString}... )
```
This function permits the user to add a command that sets the `GROUP BY` clause of a `Select` struct to the columns in `columns`. Note that the function doesn't check if these columns exist in the query result.

The output of this function is a function that takes a `Select` as argument and returns the updated argument, such that Julia's function chaining syntax can be used.

Example:
```
sq |> setgroupby!( "var1", "var2" )
```
sets the command to group the result of the `Select` statement by `var1` and `var2`.
"""
setgroupby!( columns::Vector{Var} ) = sq::Select -> begin
    sq.groupby = columns
    sq
end
setgroupby!( columns::Vector{T} ) where T <: AbstractString = Var.(columns) |> setgroupby!
setgroupby!( columns::SVar... ) = Var.( columns |> collect ) |> setgroupby!

"""
```
sethavingclause!( cond::SQLCondition )
```
This function permits the user to add a command that sets the `HAVING` clause of a `Select` struct to `cond`.

The output of this function is a function that takes a `Select` as argument and returns the updated argument, such that Julia's function chaining syntax can be used.

Example:
```
sq |> sethavingclause!( cond )
```
sets the command to consider only the fields that satisfy the condition `cond` when grouping the results of the `Select` statement.
"""
sethavingclause!( cond::SQLCondition ) = sq::Select -> begin
    sq.having = cond
    sq
end

"""
```
select( conn::Conn, df::Bool=true )
```
This function permits the user to execute the SQL equivalent of a `Select` struct on the database that `conn` connects to. The flag `df` determines the format of the output (`DataFrame` or `NamedTuple`) if the statement will return an output.

The output of this function is a function that takes a `Struct` as argument and returns the result of the executed SQL statement.

Example:
```
sq |> select(conn)
```
"""
select( conn::Conn, df::Bool=true ) = sq::Select -> begin
    res = LibPQ.execute( conn, string( sq, ";" ) )
    res |> makeunique! |> tableresult(df)
end


text:
# ---
# title: 169. Majority Element
# id: problem169
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Divide and Conquer, Bit Manipulation
# link: <https://leetcode.com/problems/majority-element/description/>
# hidden: true
# ---
#
# Given an array of size _n_ , find the majority element. The majority element
# is the element that appears **more than** `⌊ n/2 ⌋` times.
#
# You may assume that the array is non-empty and the majority element always
# exist in the array.
#
# **Example 1:**
#
#
#
#     Input: [3,2,3]
#     Output: 3
#
# **Example 2:**
#
#
#
#     Input: [2,2,1,1,1,2,2]
#     Output: 2
#
#
#
## @lc code=start
using LeetCode

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


text:
module Generator

using WAV
using DSP
using Match

export generate_signal, generate_frecs

# Maps digits with its corresponding
# output frecuencies
get_frecs(digit, x, silence_samples) = @match digit begin
	"0" => sen(x, 941) .+ sen(x, 1336)
	"1" => sen(x, 697) .+ sen(x, 1209)
	"2" => sen(x, 697) .+ sen(x, 1336)
	"3" => sen(x, 697) .+ sen(x, 1477)
	"4" => sen(x, 770) .+ sen(x, 1209)
	"5" => sen(x, 770) .+ sen(x, 1336)
	"6" => sen(x, 770) .+ sen(x, 1477)
	"7" => sen(x, 852) .+ sen(x, 1209)
	"8" => sen(x, 852) .+ sen(x, 1336)
	"9" => sen(x, 852) .+ sen(x, 1477)
	"*" => sen(x, 941) .+ sen(x, 1209)
	"#" => sen(x, 941) .+ sen(x, 1477)
	"A" => sen(x, 697) .+ sen(x, 1633)
	"B" => sen(x, 770) .+ sen(x, 1633)
	"C" => sen(x, 852) .+ sen(x, 1633)
	"D" => sen(x, 941) .+ sen(x, 1633)
	"_" => zeros(silence_samples)
end

function sen(x, frec)
	return sin.((2*pi*frec).*x)
end

function parse_sequence(sequence)

	# Starts with a silence
	# for the decoder to
	# work better
	s = ["_"]

	for c in sequence
		push!(s, string(c))
	end

	# Ends with a silence
	# because of the same reason
	# as the start
	push!(s, "_")

	return s
end

# Generates a new dialing signal:
#
# - 'sequence': tones to be generated, it
# 		includes the silences
#
# 		It cames as a 'string' so it has
# 		to be parsed into an array of the
# 		following form
#
# 	example:
# 		1) ["4", "_", "7", "_", "4", "_", "4"] (4_7_4_4)
# 		2) ["4", "_", "_", "7", "_", "4", "_", "_", "4"] (4__7_4__4)
# 		3) ["*", "_", "1", "_", "5", "_", "0", "#"] (*_1_5_0_#)
# - 'time': tones duration
# - 'sample_frec': sample frecuency
#
function generate_signal(sequence; time=70, sample_frec=8000, silence_time=70, noise=false, out_file="tone.wav")

	# Leave them in miliseconds
	time /= 1000
	silence_time /= 1000

	sequence = parse_sequence(sequence)

	res = []

	x = (0:(1/sample_frec):time)

	silence_samples = Integer(sample_frec * silence_time) + 1

	for digit in sequence
		y = get_frecs(digit, x, silence_samples)
		res = [res; y]
	end

	# Add noise
	noise_signal = []
	if noise
		noise_power = 0.00001*sample_frec/2

		noise_signal = gaussian(length(res), sqrt(noise_power))
		exp_drec = exp.((0:(length(noise_signal) - 1)) ./ (-sample_frec*5))

		noise_signal .*= exp_drec

		res .+= noise_signal
	end

	# Transform it in a 64 bit float 'Array'
	res = Array{Float64,1}(res)

	# Store it in a WAV file
	wavwrite(res, out_file, Fs=sample_frec)

	return res
end

# Generates a new signal based on the
# frequencies in the 'tones' array
#
# - freqs: [(f1,f2), (f1,f2), ..]
function generate_frecs(tones; time=70, silence_time=70)

	sample_frec = 8000
	time /= 1000
	silence_time /= 1000

	res = []

	x = (0:(1/sample_frec):time)

	silence_samples = Integer(sample_frec * silence_time) + 1

	for freqs in tones

		# Tone generation
		y = sen(x, freqs[1]) .+ sen(x, freqs[2])
		res = [res; y]

		# Silence generation
		y = zeros(silence_samples)
		res = [res; y]
	end

	res = Array{Float64,1}(res)

	wavwrite(res, "tone-$tones.wav", Fs=sample_frec)
end

end # module

max_stars_repo_name: PatricioIribarneCatella/TP-Telefono

max_issues_repo_name: PatricioIribarneCatella/TP-Telefono

max_forks_repo_name: PatricioIribarneCatella/TP-Telefono


text:
immutable CsvFile
    filename::String
    delim_char::Char
    quote_char::Char
    escape_char::Char
    header::Bool

    function CsvFile(filename, delim_char=',', quote_char='"', escape_char='\\', header=true)
        new(filename, delim_char, quote_char, escape_char, header)
    end
end

function _writevalue(io::IO, value::String, file::CsvFile)
    print(io, file.quote_char)
    for c in value
        if c==file.quote_char
            print(io, file.escape_char)
        end
        print(io, c)
    end
    print(io, file.quote_char)
end

function _writevalue(io::IO, value, file::CsvFile)
    print(io, value)
end

@generated function _writecsv{T}(io::IO, enumerable::Enumerable{T}, csvfile::CsvFile)
    col_names = fieldnames(T)
    n = length(col_names)
    push_exprs = Expr(:block)
    for i in 1:n
        push!(push_exprs.args, :( _writevalue(io, i.$(col_names[i]), csvfile) ))
        if i<n
            push!(push_exprs.args, :( print(io, csvfile.delim_char ) ))
        end
    end
    push!(push_exprs.args, :( println(io) ))

    quote
        for i in enumerable
            $push_exprs
        end
    end
end

function collect{T<:NamedTuple}(enumerable::Enumerable{T}, file::CsvFile)
    open(file.filename, "w") do io
        if file.header
            join(io,["$(file.quote_char)" *replace(string(colname), file.quote_char, "$(file.escape_char)$(file.quote_char)") * "$(file.quote_char)" for colname in fieldnames(T)],file.delim_char)
            println(io)
        end
        _writecsv(io, enumerable, file)
    end
    return nothing
end


text:
# NOTE: This file contains machine-generated code by translate_exported_functions.jl.

# int hts_resize_array_(size_t, size_t, size_t, void *, void **, int, const char *);
# ErrorException("found unnamed argument(s)")

# void hts_free(void *ptr);
@defun hts_free(ptr::Ptr{Cvoid})::Cvoid

# int hts_opt_add(hts_opt **opts, const char *c_arg);
@defun hts_opt_add(opts::Ptr{Ptr{hts_opt}}, c_arg::Cstring)::Cint

# int hts_opt_apply(htsFile *fp, hts_opt *opts);
@defun hts_opt_apply(fp::Ptr{htsFile}, opts::Ptr{hts_opt})::Cint

# void hts_opt_free(hts_opt *opts);
@defun hts_opt_free(opts::Ptr{hts_opt})::Cvoid

# int hts_parse_format(htsFormat *opt, const char *str);
@defun hts_parse_format(opt::Ptr{htsFormat}, str::Cstring)::Cint

# int hts_parse_opt_list(htsFormat *opt, const char *str);
@defun hts_parse_opt_list(opt::Ptr{htsFormat}, str::Cstring)::Cint

# const char *hts_version(void);
@defun hts_version()::Cstring

# int hts_detect_format(struct hFILE *fp, htsFormat *fmt);
@defun hts_detect_format(fp::Ptr{hFILE}, fmt::Ptr{htsFormat})::Cint

# char *hts_format_description(const htsFormat *format);
@defun hts_format_description(format::Ptr{htsFormat})::Ptr{Cchar}

# htsFile *hts_open(const char *fn, const char *mode);
@defun hts_open(fn::Cstring, mode::Cstring)::Ptr{htsFile}

# htsFile *hts_open_format(const char *fn, const char *mode, const htsFormat *fmt);
@defun hts_open_format(fn::Cstring, mode::Cstring, fmt::Ptr{htsFormat})::Ptr{htsFile}

# htsFile *hts_hopen(struct hFILE *fp, const char *fn, const char *mode);
@defun hts_hopen(fp::Ptr{hFILE}, fn::Cstring, mode::Cstring)::Ptr{htsFile}

# int hts_close(htsFile *fp);
@defun hts_close(fp::Ptr{htsFile})::Cint

# const htsFormat *hts_get_format(htsFile *fp);
@defun hts_get_format(fp::Ptr{htsFile})::Ptr{htsFormat}

# const char *hts_format_file_extension(const htsFormat *format);
@defun hts_format_file_extension(format::Ptr{htsFormat})::Cstring

# int hts_set_opt(htsFile *fp, enum hts_fmt_option opt, ...);
# ErrorException("unknown token: '.'")

# int hts_getline(htsFile *fp, int delimiter, kstring_t *str);
@defun hts_getline(fp::Ptr{htsFile}, delimiter::Cint, str::Ptr{kstring_t})::Cint

# char **hts_readlines(const char *fn, int *_n);
@defun hts_readlines(fn::Cstring, _n::Ptr{Cint})::Ptr{Ptr{Cchar}}

# char **hts_readlist(const char *fn, int is_file, int *_n);
@defun hts_readlist(fn::Cstring, is_file::Cint, _n::Ptr{Cint})::Ptr{Ptr{Cchar}}

# int hts_set_threads(htsFile *fp, int n);
@defun hts_set_threads(fp::Ptr{htsFile}, n::Cint)::Cint

# int hts_set_thread_pool(htsFile *fp, htsThreadPool *p);
@defun hts_set_thread_pool(fp::Ptr{htsFile}, p::Ptr{htsThreadPool})::Cint

# void hts_set_cache_size(htsFile *fp, int n);
@defun hts_set_cache_size(fp::Ptr{htsFile}, n::Cint)::Cvoid

# int hts_set_fai_filename(htsFile *fp, const char *fn_aux);
@defun hts_set_fai_filename(fp::Ptr{htsFile}, fn_aux::Cstring)::Cint

# int hts_check_EOF(htsFile *fp);
@defun hts_check_EOF(fp::Ptr{htsFile})::Cint

# hts_idx_t *hts_idx_init(int n, int fmt, uint64_t offset0, int min_shift, int n_lvls);
@defun hts_idx_init(n::Cint, fmt::Cint, offset0::UInt64, min_shift::Cint, n_lvls::Cint)::Ptr{hts_idx_t}

# void hts_idx_destroy(hts_idx_t *idx);
@defun hts_idx_destroy(idx::Ptr{hts_idx_t})::Cvoid

# int hts_idx_push(hts_idx_t *idx, int tid, hts_pos_t beg, hts_pos_t end, uint64_t offset, int is_mapped);
@defun hts_idx_push(idx::Ptr{hts_idx_t}, tid::Cint, beg::hts_pos_t, end_::hts_pos_t, offset::UInt64, is_mapped::Cint)::Cint

# int hts_idx_finish(hts_idx_t *idx, uint64_t final_offset);
@defun hts_idx_finish(idx::Ptr{hts_idx_t}, final_offset::UInt64)::Cint

# int hts_idx_fmt(hts_idx_t *idx);
@defun hts_idx_fmt(idx::Ptr{hts_idx_t})::Cint

# int hts_idx_tbi_name(hts_idx_t *idx, int tid, const char *name);
@defun hts_idx_tbi_name(idx::Ptr{hts_idx_t}, tid::Cint, name::Cstring)::Cint

# int hts_idx_save(const hts_idx_t *idx, const char *fn, int fmt) HTS_RESULT_USED;
@defun hts_idx_save(idx::Ptr{hts_idx_t}, fn::Cstring, fmt::Cint)::Cint

# int hts_idx_save_as(const hts_idx_t *idx, const char *fn, const char *fnidx, int fmt) HTS_RESULT_USED;
@defun hts_idx_save_as(idx::Ptr{hts_idx_t}, fn::Cstring, fnidx::Cstring, fmt::Cint)::Cint

# hts_idx_t *hts_idx_load(const char *fn, int fmt);
@defun hts_idx_load(fn::Cstring, fmt::Cint)::Ptr{hts_idx_t}

# hts_idx_t *hts_idx_load2(const char *fn, const char *fnidx);
@defun hts_idx_load2(fn::Cstring, fnidx::Cstring)::Ptr{hts_idx_t}

# hts_idx_t *hts_idx_load3(const char *fn, const char *fnidx, int fmt, int flags);
@defun hts_idx_load3(fn::Cstring, fnidx::Cstring, fmt::Cint, flags::Cint)::Ptr{hts_idx_t}

# uint8_t *hts_idx_get_meta(hts_idx_t *idx, uint32_t *l_meta);
@defun hts_idx_get_meta(idx::Ptr{hts_idx_t}, l_meta::Ptr{UInt32})::Ptr{UInt8}

# int hts_idx_set_meta(hts_idx_t *idx, uint32_t l_meta, uint8_t *meta, int is_copy);
@defun hts_idx_set_meta(idx::Ptr{hts_idx_t}, l_meta::UInt32, meta::Ptr{UInt8}, is_copy::Cint)::Cint

# int hts_idx_get_stat(const hts_idx_t* idx, int tid, uint64_t* mapped, uint64_t* unmapped);
@defun hts_idx_get_stat(idx::Ptr{hts_idx_t}, tid::Cint, mapped::Ptr{UInt64}, unmapped::Ptr{UInt64})::Cint

# uint64_t hts_idx_get_n_no_coor(const hts_idx_t* idx);
@defun hts_idx_get_n_no_coor(idx::Ptr{hts_idx_t})::UInt64

# long long hts_parse_decimal(const char *str, char **strend, int flags);
@defun hts_parse_decimal(str::Cstring, strend::Ptr{Ptr{Cchar}}, flags::Cint)::Clonglong

# const char *hts_parse_reg64(const char *str, hts_pos_t *beg, hts_pos_t *end);
@defun hts_parse_reg64(str::Cstring, beg::Ptr{hts_pos_t}, end_::Ptr{hts_pos_t})::Cstring

# const char *hts_parse_reg(const char *str, int *beg, int *end);
@defun hts_parse_reg(str::Cstring, beg::Ptr{Cint}, end_::Ptr{Cint})::Cstring

# const char *hts_parse_region(const char *s, int *tid, hts_pos_t *beg, hts_pos_t *end, hts_name2id_f getid, void *hdr, int flags);
@defun hts_parse_region(s::Cstring, tid::Ptr{Cint}, beg::Ptr{hts_pos_t}, end_::Ptr{hts_pos_t}, getid::hts_name2id_f, hdr::Ptr{Cvoid}, flags::Cint)::Cstring

# hts_itr_t *hts_itr_query(const hts_idx_t *idx, int tid, hts_pos_t beg, hts_pos_t end, hts_readrec_func *readrec);
@defun hts_itr_query(idx::Ptr{hts_idx_t}, tid::Cint, beg::hts_pos_t, end_::hts_pos_t, readrec::Ptr{hts_readrec_func})::Ptr{hts_itr_t}

# void hts_itr_destroy(hts_itr_t *iter);
@defun hts_itr_destroy(iter::Ptr{hts_itr_t})::Cvoid

# hts_itr_t *hts_itr_querys(const hts_idx_t *idx, const char *reg, hts_name2id_f getid, void *hdr, hts_itr_query_func *itr_query, hts_readrec_func *readrec);
@defun hts_itr_querys(idx::Ptr{hts_idx_t}, reg::Cstring, getid::hts_name2id_f, hdr::Ptr{Cvoid}, itr_query::Ptr{hts_itr_query_func}, readrec::Ptr{hts_readrec_func})::Ptr{hts_itr_t}

# int hts_itr_next(BGZF *fp, hts_itr_t *iter, void *r, void *data) HTS_RESULT_USED;
@defun hts_itr_next(fp::Ptr{BGZF}, iter::Ptr{hts_itr_t}, r::Ptr{Cvoid}, data::Ptr{Cvoid})::Cint

# const char **hts_idx_seqnames(const hts_idx_t *idx, int *n, hts_id2name_f getid, void *hdr);
@defun hts_idx_seqnames(idx::Ptr{hts_idx_t}, n::Ptr{Cint}, getid::hts_id2name_f, hdr::Ptr{Cvoid})::Ptr{Cstring}

# int hts_itr_multi_bam(const hts_idx_t *idx, hts_itr_t *iter);
@defun hts_itr_multi_bam(idx::Ptr{hts_idx_t}, iter::Ptr{hts_itr_t})::Cint

# int hts_itr_multi_cram(const hts_idx_t *idx, hts_itr_t *iter);
@defun hts_itr_multi_cram(idx::Ptr{hts_idx_t}, iter::Ptr{hts_itr_t})::Cint

# hts_itr_t *hts_itr_regions(const hts_idx_t *idx, hts_reglist_t *reglist, int count, hts_name2id_f getid, void *hdr, hts_itr_multi_query_func *itr_specific, hts_readrec_func *readrec, hts_seek_func *seek, hts_tell_func *tell);
@defun hts_itr_regions(idx::Ptr{hts_idx_t}, reglist::Ptr{hts_reglist_t}, count::Cint, getid::hts_name2id_f, hdr::Ptr{Cvoid}, itr_specific::Ptr{hts_itr_multi_query_func}, readrec::Ptr{hts_readrec_func}, seek::Ptr{hts_seek_func}, tell::Ptr{hts_tell_func})::Ptr{hts_itr_t}

# int hts_itr_multi_next(htsFile *fd, hts_itr_t *iter, void *r);
@defun hts_itr_multi_next(fd::Ptr{htsFile}, iter::Ptr{hts_itr_t}, r::Ptr{Cvoid})::Cint

# hts_reglist_t *hts_reglist_create(char **argv, int argc, int *r_count, void *hdr, hts_name2id_f getid);
@defun hts_reglist_create(argv::Ptr{Ptr{Cchar}}, argc::Cint, r_count::Ptr{Cint}, hdr::Ptr{Cvoid}, getid::hts_name2id_f)::Ptr{hts_reglist_t}

# void hts_reglist_free(hts_reglist_t *reglist, int count);
@defun hts_reglist_free(reglist::Ptr{hts_reglist_t}, count::Cint)::Cvoid

# int hts_file_type(const char *fname);
@defun hts_file_type(fname::Cstring)::Cint

# errmod_t *errmod_init(double depcorr);
@defun errmod_init(depcorr::Cdouble)::Ptr{errmod_t}

# void errmod_destroy(errmod_t *em);
@defun errmod_destroy(em::Ptr{errmod_t})::Cvoid

# int errmod_cal(const errmod_t *em, int n, int m, uint16_t *bases, float *q);
@defun errmod_cal(em::Ptr{errmod_t}, n::Cint, m::Cint, bases::Ptr{UInt16}, q::Ptr{Cfloat})::Cint

# int probaln_glocal(const uint8_t *ref, int l_ref, const uint8_t *query, int l_query, const uint8_t *iqual, const probaln_par_t *c, int *state, uint8_t *q);
@defun probaln_glocal(ref::Ptr{UInt8}, l_ref::Cint, query::Ptr{UInt8}, l_query::Cint, iqual::Ptr{UInt8}, c::Ptr{probaln_par_t}, state::Ptr{Cint}, q::Ptr{UInt8})::Cint

# hts_md5_context *hts_md5_init(void);
@defun hts_md5_init()::Ptr{hts_md5_context}

# void hts_md5_update(hts_md5_context *ctx, const void *data, unsigned long size);
@defun hts_md5_update(ctx::Ptr{hts_md5_context}, data::Ptr{Cvoid}, size::Culong)::Cvoid

# void hts_md5_final(unsigned char *digest, hts_md5_context *ctx);
@defun hts_md5_final(digest::Ptr{Cuchar}, ctx::Ptr{hts_md5_context})::Cvoid

# void hts_md5_reset(hts_md5_context *ctx);
@defun hts_md5_reset(ctx::Ptr{hts_md5_context})::Cvoid

# void hts_md5_hex(char *hex, const unsigned char *digest);
@defun hts_md5_hex(hex::Ptr{Cchar}, digest::Ptr{Cuchar})::Cvoid

# void hts_md5_destroy(hts_md5_context *ctx);
@defun hts_md5_destroy(ctx::Ptr{hts_md5_context})::Cvoid


text:
# Autogenerated wrapper script for msolve_jll for x86_64-apple-darwin
export libmsolve, libneogb

using GMP_jll
using FLINT_jll
using MPFR_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("msolve")
JLLWrappers.@declare_library_product(libmsolve, "@rpath/libmsolve-0.1.7.dylib")
JLLWrappers.@declare_library_product(libneogb, "@rpath/libneogb-0.1.7.dylib")
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll, FLINT_jll, MPFR_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libmsolve,
        "lib/libmsolve.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIB = MathOptInterface.Bridges
const MOIU = MathOptInterface.Utilities

include("sddp_bridge.jl")
include("ddp_bridge.jl")
include("lazy_relax.jl")


text:
using Oceananigans
using Printf
using Statistics
using GLMakie

using Oceananigans.ImmersedBoundaries: mask_immersed_field!
using Oceananigans.Utils: prettysummary

# Monin-Obukhov drag coefficient
z₀ = 1e-2 # Charnock roughness
κ = 0.4 # Von Karman constant
Cᴰ(Δz) = (κ / log(Δz / 2z₀))^2

@inline bottom_drag_u(x, y, t, u, w, Cᴰ) = - Cᴰ * u * sqrt(u^2 + w^2)
@inline bottom_drag_w(x, y, t, u, w, Cᴰ) = - Cᴰ * w * sqrt(u^2 + w^2)
@inline bottom_drag_u(x, y, z, t, u, w, Cᴰ) = - Cᴰ * u * sqrt(u^2 + w^2)
@inline bottom_drag_w(x, y, z, t, u, w, Cᴰ) = - Cᴰ * w * sqrt(u^2 + w^2)

function hilly_simulation(; Nx = 64,
                            Nz = Nx,
                            h = 0.1,
                            Re = 1e3,
                            N² = 1e-2,
                            boundary_condition = :no_slip,
                            stop_time = 1,
                            save_interval = 0.1,
                            architecture = CPU(),
                            filename = "flow_over_hills")

    underlying_grid = RectilinearGrid(architecture, size = (Nx, Nz), halo = (3, 3),
                                      x = (0, 2π), z = (0, 1),
                                      topology = (Periodic, Flat, Bounded))

    min_Δz = 1 / 32
    if h > 0
        hills(x, y) = h * (1 + sin(x)) / 2 + 2min_Δz # ensure that the bottom boundary is immersed
        grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(hills))
    else # no hills
        grid = underlying_grid
    end

    closure = isfinite(Re) ? ScalarDiffusivity(ν=1/Re, κ=1/Re) : nothing

    if boundary_condition == :no_slip
        no_slip = ValueBoundaryCondition(0)
        u_bcs = FieldBoundaryConditions(bottom=no_slip, immersed=no_slip)
        w_bcs = FieldBoundaryConditions(immersed=no_slip)
        boundary_conditions = (; u = u_bcs, w = w_bcs)
    elseif boundary_condition == :bottom_drag
        Δz = 1 / Nz
        Δx = 2π / Nz
        u_drag_bc = FluxBoundaryCondition(bottom_drag_u, field_dependencies=(:u, :w), parameters=Cᴰ(Δz))
        w_drag_bc = FluxBoundaryCondition(bottom_drag_w, field_dependencies=(:u, :w), parameters=Cᴰ(Δx))
        u_bcs = FieldBoundaryConditions(bottom=u_drag_bc, immersed=u_drag_bc)
        w_bcs = FieldBoundaryConditions(immersed=w_drag_bc)
        boundary_conditions = (; u = u_bcs, w = w_bcs)
        @info string("Using a bottom drag with coefficient ", Cᴰ(Δz))
    else
        boundary_conditions = NamedTuple()
    end

    model = NonhydrostaticModel(; grid, closure, boundary_conditions,
                                advection = WENO5(),
                                timestepper = :RungeKutta3,
                                tracers = :b,
                                buoyancy = BuoyancyTracer())

    # Steady flow + perturbations
    δh = 0.1
    ∂z_ψᵋ(x, z) = 4π * sin(4x) * cos(4π * z) * exp(-(z - h)^2 / 2δh^2)
    ∂x_ψᵋ(x, z) = 4  * cos(4x) * sin(4π * z) * exp(-(z - h)^2 / 2δh^2)
    bᵢ(x, y, z) = N² * z + 1e-9 * rand()
    uᵢ(x, y, z) = 1.0 + ∂z_ψᵋ(x, z)
    wᵢ(x, y, z) = - ∂x_ψᵋ(x, z)
    set!(model, b=bᵢ, u=uᵢ, w=wᵢ)

    Δx = 2π / Nx
    Δt = 0.1 * Δx
    simulation = Simulation(model; Δt, stop_time)

    u, v, w = model.velocities
    Uᵢ = mean(u)

    wall_clock = Ref(time_ns())

    function progress(sim)
        δU = mean(u) / Uᵢ
        elapsed = 1e-9 * (time_ns() - wall_clock[])
        @info @sprintf("Iter: %d, time: %.2e, δU: %.2e, max|w|: %.2e, wall time: %s",
                       iteration(sim), time(sim), δU, maximum(abs, w), prettytime(elapsed))
        wall_clock[] = time_ns()
        return nothing
    end

    simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))

    wizard = TimeStepWizard(cfl=0.5)
    simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

    U = Average(u, dims=(1, 2, 3))
    ξ = ∂z(u) - ∂x(w)

    ke = @at (Center, Center, Center) (u^2 + v^2 + w^2) / 2
    KE = Average(ke, dims=(1, 2, 3))

    simulation.output_writers[:fields] =
        JLD2OutputWriter(model, merge(model.velocities, model.tracers, (; ξ, U, KE));
                         schedule = TimeInterval(save_interval),
                         with_halos = true,
                         filename,
                         overwrite_existing = true)

    @info "Made a simulation of"
    @show model

    @info "The grid is"
    @show model.grid

    @info "The x-velocity immersed and bottom boundary conditions are"
    @show boundary_condition
    @show model.velocities.u.boundary_conditions.immersed
    @show model.velocities.u.boundary_conditions.bottom

    return simulation
end

#####
##### Run them!
#####

Nx = 128
stop_time = 20.0

experiments = ["reference", "no_slip", "free_slip", "bottom_drag"]
Nexp = length(experiments)

function filename(exp, Nx, h)
    exp == "reference" && (h = 0.0)
    return @sprintf("hills_%s_%d_h%d", exp, Nx, 10h)
end

for exp in experiments
    h = exp == "reference" ? 0.0 : 0.3
    boundary_condition = exp == "reference" ? :no_slip : Symbol(exp)
    reference_sim = hilly_simulation(; stop_time, Nx, filename=filename(exp, Nx, h), h, boundary_condition)
    run!(reference_sim)
end

#####
##### Plot results
#####

h = 0.3
ξ  = Dict(exp => FieldTimeSeries(filename(exp, Nx, h) * ".jld2", "ξ")  for exp in experiments)
U  = Dict(exp => FieldTimeSeries(filename(exp, Nx, h) * ".jld2", "U")  for exp in experiments)
KE = Dict(exp => FieldTimeSeries(filename(exp, Nx, h) * ".jld2", "KE") for exp in experiments)

Nt = minimum(length(ξ[exp].times) for exp in experiments)
t = ξ["reference"].times[1:Nt]
δU_series(U) = [(U[1, 1, 1, n] - U[1, 1, 1, 1]) / U[1, 1, 1, 1] for n = 1:Nt]
δK_series(K) = [(K[1, 1, 1, n] - K[1, 1, 1, 1]) / K[1, 1, 1, 1] for n = 1:Nt]
δU = Dict(exp => δU_series(u) for (exp, u) in U)
δKE = Dict(exp => δK_series(k) for (exp, k) in KE)

fig = Figure(resolution=(1800, 1200))
ax = Dict(exp => Axis(fig[i+1, 2:4], aspect=2π, xlabel="x", ylabel="z", title=exp)
          for (i, exp) in enumerate(experiments))

axu = Axis(fig[2:3, 5], xlabel="t", ylabel="Total momentum")
axe = Axis(fig[4:5, 5], xlabel="t", ylabel="Total kinetic energy")

slider = Slider(fig[Nexp+2, 2:4], range=1:Nt, startvalue=1)
n = slider.value

# Title
title = @lift string("Flow over hills at t = ", prettysummary(t[$n]))
Label(fig[1, 1:5], title)

# Vorticity heatmaps
ξi(ξ) = @lift begin
    ξn = ξ[$n]
    mask_immersed_field!(ξn, NaN)
    interior(ξn, :, 1, :)
end

ξⁱ = Dict(exp => ξi(ξ[exp]) for exp in experiments)
x, y, z = nodes(ξ["reference"])
ξmax = maximum(abs, ξ["reference"])
ξlim = ξmax / 50

hm = Dict(exp => heatmap!(ax[exp], x, z, ξⁱ[exp], colorrange=(-ξlim, ξlim), colormap=:redblue)
          for exp in experiments)

cb = Colorbar(fig[2:5, 1], hm["reference"], vertical=true, flipaxis=true, label="Vorticity, ∂z(u) - ∂x(w)")

# Momentum and energy plots
for exp in experiments
    lines!(axu, t, δU[exp], label=exp)
    lines!(axe, t, δKE[exp], label=exp)
end

axislegend(axu, position=:lb)
axislegend(axe, position=:rt)

tn = @lift t[$n]
min_δU = minimum(minimum(δ) for δ in values(δU))
min_δK = minimum(minimum(δ) for δ in values(δKE))
vlines!(axu, tn, ymin=min_δU, ymax=1.0)
vlines!(axe, tn, ymin=min_δK, ymax=1.0)

display(fig)

moviename = @sprintf("flow_over_hills_%dd_h%d.mp4", Nx, 10h)
record(fig, moviename, 1:Nt, framerate=24) do nn
    n[] = nn
end


text:
function f(n)
    a = inv(rand(n,n))
    out = 0.0
    for j=1:n
        for i=1:n
            out += a[i,j]
        end
    end
    out
end
include("test2.jl")

f(1)
g(1)

Profile.clear()
@profile f(3*10^3)
@profile g(10^2)
using ProfileFile
ProfileFile.write_to_file()


text:
using DelimitedFiles
using OffsetArrays
using SpecialFunctions
using StatsBase

data = readdlm("coal-mining-disaster.txt")
α = 1
β = 200/365.24
L = data[end] - data[1]
coal = data .- data[1]
λ = 3

# count number of disaster in [s1, s2)
function count_disasters(s1::Float64, s2::Float64)
    return(sum( (coal .>= s1) .& (coal .< s2) ))
end

# accept probability for H move
function alphaH(s::OffsetArray, h::OffsetArray, hjp::Float64, j::Int)
    m = count_disasters(s[j], s[j+1])
    log_ll_r = m * (log(hjp) - log(h[j])) - (hjp - h[j]) * (s[j+1] - s[j])
    log_prior_r = α * (log(hjp) - log(h[j])) - β * (hjp - h[j])
    return log_ll_r + log_prior_r
end

# accept probability for P move
function alphaP(s::OffsetArray, h::OffsetArray, sjp::Float64, j::Int)
    m1 = count_disasters(s[j-1], s[j])
    m1p = count_disasters(s[j-1], sjp)
    m2 = count_disasters(s[j], s[j+1])
    m2p = count_disasters(sjp, s[j+1])
    log_ll_r = (m1p - m1) * log(h[j-1]) + (m2p - m2) * log(h[j]) - (sjp - s[j]) * (h[j-1] - h[j])
    log_prior_r = log(s[j+1] - sjp) + log(sjp - s[j-1]) - log(s[j+1] - s[j]) - log(s[j] - s[j-1])
    return log_ll_r + log_prior_r
end

# accept probability for birth move
function alphaB(s::OffsetArray, h::OffsetArray, sstar, h1p, h2p, j, k)
    m1p = count_disasters(s[j], sstar)
    m2p = count_disasters(sstar, s[j+1])
    mj = m1p + m2p
    log_ll_r = (m1p*log(h1p) + m2p*log(h2p) - (m1p + m2p)*log(h[j])) - (h1p*(sstar-s[j])+h2p*(s[j+1]-sstar)-h[j]*(s[j+1]-s[j]))
    log_prior_r = log(2(2k+3)λ) - 2log(L) + log(sstar-s[j]) + log(s[j+1]-sstar) - log(s[j+1]-s[j]) + α*log(β) - log(gamma(α)) + (α - 1) *(log(h1p) + log(h2p) - log(h[j])) - β*(h1p + h2p - h[j])
    log_proposal_r = log(death(k+1)) + log(L) - log(birth(k)) - log(k+1)
    log_Jacobian = 2 * log(h1p + h2p) - log(h[j])
    return log_ll_r + log_prior_r + log_proposal_r + log_Jacobian
end

# accept probability for death move
function alphaD(s::OffsetArray, h::OffsetArray, hjp, j, k)
    m1 = count_disasters(s[j], s[j+1])
    m2 = count_disasters(s[j+1], s[j+2])
    mjp = m1 + m2
    log_ll_r = ( (m1+m2)*log(hjp) - m1*log(h[j]) - m2*log(h[j+1])) - (hjp*(s[j+2]-s[j]) - (h[j]*(s[j+1]-s[j])) - h[j+1]*(s[j+2]-s[j+1]))
    log_prior_r = - (log(2(2k+3)λ) - 2log(L) + log(s[j+1]-s[j]) + log(s[j+2]-s[j+1]) - log(s[j+2]-s[j]) + α*log(β) - log(gamma(α)) + (α - 1) *(log(h[j]) + log(h[j+1]) - log(hjp)) - β*(h[j] + h[j+1] - hjp))
    log_proposal_r = -( log(death(k+1)) + log(L) - log(birth(k)) - log(k+1) )
    log_Jacobian = log(hjp) - 2 * log(h[j] + h[j+1])
    return log_ll_r + log_prior_r + log_proposal_r + log_Jacobian
end

function birth(k; c = 3.6/7)
    return c * min(1, λ/(k+1))
end

function death(k; c = 3.6/7)
    if k == 0
        return 0
    else
        return c * min(1, (k+1)/λ)
    end
end

function move(n::Int, h::OffsetArray, s::OffsetArray, k::Int)
    H = []
    S = []
    println(s)
    K = ones(n)
    K[1] = k
    push!(H, copy(h))
    push!(S, copy(s))
    for i = 2:n
        println(i)
        probP = 0
        if k > 0
            probP = 0.5 * (1 - birth(k) - death(k))
        end
        probH = 1 - birth(k) - death(k) - probP
        flag = rand()
        println("ok")
        if flag < probH
            println("ok1")
            j = sample(0:k)
            u = rand() - 0.5
            hjp = h[j] * exp(u)
            if log(rand()) < alphaH(s, h, hjp, j)
                h[j] = hjp
            end
        elseif flag < probH + probP
            println("ok2")
            j = sample(1:k)
            sjp = rand() * (s[j+1]-s[j-1]) + s[j-1]
            if log(rand()) < alphaP(s, h, sjp, j)
                s[j] = sjp
            end
        elseif flag < probH + probP + birth(k)
            println("ok3")
            sstar = rand() * L
            j = findfirst(s .>= sstar) - 1
            println(j)
            u = rand()
            h1p = h[j] * exp(-(s[j+1]-sstar)/(s[j+1]-s[j]) * log((1-u)/u))
            h2p = (1-u)/u*h1p
            println("ok31")
            if log(rand()) < alphaB(s, h, sstar, h1p, h2p, j, k)
                push!(s.parent, 0)
                s[(j+2):end] .= s[(j+1):(end-1)]
                s[j+1] = sstar
                push!(h.parent, 0)
                h[(j+2):end] .= h[(j+1):(end-1)]
                h[j+1] = h2p
                h[j] = h1p
                k = k + 1
            end
        else
            println("ok4")
            j = sample(1:k) - 1
            r = ( s[j+1] - s[j] ) / (s[j+2] - s[j] )
            hjp = h[j]^r * h[j+1]^(1-r)
            if log(rand()) < alphaD(s, h, hjp, j, k)
                k = k - 1
                h[j] = hjp
                # the first plus 1 for converting index
                # the seconde plus 1 for the removed one h_{j+1}
                deleteat!(h.parent, j+1+1)
                deleteat!(s.parent, j+1+1)
            end
        end
        K[i] = k
        push!(H, copy(h))
        push!(S, copy(s))
    end
    return K, H, S
end

K, H, S = move(44000, OffsetArray([1.7], -1), OffsetArray([0, L], -1), 0)

# results
using Plots
K = K[4001:end]
H = H[4001:end]
S = S[4001:end]
histogram(K, legend = false, xlab = "Number of change points")

# k = 1
S1 = vcat([parent(x)' for x in S[K.==1]]...)
S2 = vcat([parent(x)' for x in S[K.==2]]...)
S3 = vcat([parent(x)' for x in S[K.==3]]...)

using KernelDensity
y1 = kde(S1[:,2])
y2 = kde(vec(S2[:,2:3]))
y3 = kde(vec(S3[:,2:4]))
mm = extrema(union(S1[:,2], vec(S2[:,2:3]), vec(S3[:,2:4])))
p1 = plot(range(mm..., length=100), z->pdf(y1,z), title="k=1")
p2 = plot(range(mm..., length=100), z->pdf(y2,z), title="k=2")
p3 = plot(range(mm..., length=100), z->pdf(y3,z), title="k=3")
plot(p1, p2, p3, layout=(3,1), legend=false)

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
using SymSpell
using BenchmarkTools
using Statistics
using Revise

cd(@__DIR__)

SUITE = BenchmarkGroup()
SUITE["lookup"]= BenchmarkGroup(["match", "opcodes"])

s1 = "qabxcd"
s2 = "abycdf"
SUITE["lookup"]["opcodes"] = @benchmarkable SymSpell.get_opcodes($s1, $s2)

text_w_casing = "Haaw is the weeather in New York?"
text_wo_casing = "how is the weather in new york?"
SUITE["lookup"]["transfer_casing"] = @benchmarkable SymSpell.transfer_casing_for_similar_text($text_w_casing, $text_wo_casing)

# First time only
tune!(SUITE)
BenchmarkTools.save("params.json", params(SUITE))
# results = run(SUITE, verbose = true)
results = run(SUITE)
BenchmarkTools.save("results1.json", results)

# Second time (different implementation)
res1 = BenchmarkTools.load("results1.json")[1]

loadparams!(SUITE, BenchmarkTools.load("params.json")[1], :evals, :samples)
res2 = run(SUITE)

judge(median(res1), median(res2))

#####
# PkgBenchmark

using PkgBenchmark
import SymSpell
br = benchmarkpkg("SymSpell")

cd(@__DIR__)
PkgBenchmark.export_markdown("SomeTest2.md", br)

# Judge

j1 = PkgBenchmark.judge("SymSpell", "optimize_getopcodes",  "master")
export_markdown("/tmp/xxx.md", j1)


text:
module FrostHelperSpeedRingChallenge

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/SpeedRingChallenge" SpeedRingChallengeFrostHelper(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight, timeLimit::Number=1.0,name::String="fh_test", nodes::Array{Tuple{Integer, Integer}, 1}=Tuple{Integer, Integer}[])

const placements = Ahorn.PlacementDict(
    "Speed Ring Challenge (Frost Helper)" => Ahorn.EntityPlacement(
        SpeedRingChallengeFrostHelper,
        "rectangle"
    ),
)

Ahorn.nodeLimits(entity::SpeedRingChallengeFrostHelper) = 1, 99
Ahorn.minimumSize(entity::SpeedRingChallengeFrostHelper) = 8, 8
Ahorn.resizable(entity::SpeedRingChallengeFrostHelper) = true, true

function Ahorn.selection(entity::SpeedRingChallengeFrostHelper)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    nodes = get(entity.data, "nodes", ())

    if isempty(nodes)
        return Ahorn.Rectangle(x, y, width, height)

    else
        res = [Ahorn.Rectangle(x, y, width, height)]
        for node in nodes
            nx, ny = Int.(node)
            push!(res, Ahorn.Rectangle(nx, ny, width, height))
        end
        return res
    end

end

edgeColor = (38, 0, 0, 128) ./ 255
centerColor = (64, 0, 0, 10) ./ 255

function renderSpaceJam(ctx::Ahorn.Cairo.CairoContext, x::Number, y::Number, width::Number, height::Number)
    Ahorn.Cairo.save(ctx)

    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1)

    Ahorn.drawRectangle(ctx, x, y, width, height, edgeColor, centerColor)

    Ahorn.restore(ctx)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::SpeedRingChallengeFrostHelper)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    if !isempty(nodes)
        for node in nodes
            nx, ny = Int.(node)

            cox, coy = floor(Int, width / 2), floor(Int, height / 2)

            renderSpaceJam(ctx, nx, ny, width, height)
            Ahorn.drawArrow(ctx, x + cox, y + coy, nx + cox, ny + coy, Ahorn.colors.selection_selected_fc, headLength=6)
            x = nx
            y = ny
        end
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::SpeedRingChallengeFrostHelper, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    renderSpaceJam(ctx, 0, 0, width, height)
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
using Test, SparseArrays
using LDPCStorage

@testset "save and load alist" begin
    # Save and load this arbitrary matrix:
    H = Int8[
        0 0 1 1 0 0 0 0 1 0 0 1 1 0
        1 0 0 1 1 0 0 0 0 0 1 0 0 1
        0 1 0 1 0 1 1 0 1 0 0 1 1 0
        1 0 0 1 0 0 0 1 0 1 0 1 0 1
    ]
    file_path = tempname() * "_unit_test.alist"
    save_to_alist(H, file_path)

    H_loaded = load_alist(file_path)

    @test H == H_loaded
end


text:
using Oceananigans.Fields: AbstractField
using Oceananigans.Grids:
    Face, Bounded,
    interior_indices,
    left_halo_indices, right_halo_indices,
    underlying_left_halo_indices, underlying_right_halo_indices

# TODO: Move to Grids/grid_utils.jl

#####
##### Viewing halos
#####

west_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
    include_corners ? view(f.data, left_halo_indices(LX, topology(f, 1), f.grid.Nx, f.grid.Hx), :, :) :
                      view(f.data, left_halo_indices(LX, topology(f, 1), f.grid.Nx, f.grid.Hx),
                                   interior_indices(LY, topology(f, 2), f.grid.Ny),
                                   interior_indices(LZ, topology(f, 3), f.grid.Nz))

east_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
    include_corners ? view(f.data, right_halo_indices(LX, topology(f, 1), f.grid.Nx, f.grid.Hx), :, :) :
                      view(f.data, right_halo_indices(LX, topology(f, 1), f.grid.Nx, f.grid.Hx),
                                   interior_indices(LY, topology(f, 2), f.grid.Ny),
                                   interior_indices(LZ, topology(f, 3), f.grid.Nz))

south_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
    include_corners ? view(f.data, :, left_halo_indices(LY, topology(f, 2), f.grid.Ny, f.grid.Hy), :) :
                      view(f.data, interior_indices(LX, topology(f, 1), f.grid.Nx),
                                   left_halo_indices(LY, topology(f, 2), f.grid.Ny, f.grid.Hy),
                                   interior_indices(LZ, topology(f, 3), f.grid.Nz))

north_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
    include_corners ? view(f.data, :, right_halo_indices(LY, topology(f, 2), f.grid.Ny, f.grid.Hy), :) :
                      view(f.data, interior_indices(LX, topology(f, 1), f.grid.Nx),
                                   right_halo_indices(LY, topology(f, 2), f.grid.Ny, f.grid.Hy),
                                   interior_indices(LZ, topology(f, 3), f.grid.Nz))

bottom_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
    include_corners ? view(f.data, :, :, left_halo_indices(LZ, topology(f, 3), f.grid.Nz, f.grid.Hz)) :
                      view(f.data, interior_indices(LX, topology(f, 1), f.grid.Nx),
                                   interior_indices(LY, topology(f, 2), f.grid.Ny),
                                   left_halo_indices(LZ, topology(f, 3), f.grid.Nz, f.grid.Hz))

top_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
    include_corners ? view(f.data, :, :, right_halo_indices(LZ, topology(f, 3), f.grid.Nz, f.grid.Hz)) :
                      view(f.data, interior_indices(LX, topology(f, 1), f.grid.Nx),
                                   interior_indices(LY, topology(f, 2), f.grid.Ny),
                                   right_halo_indices(LZ, topology(f, 3), f.grid.Nz, f.grid.Hz))

underlying_west_halo(f, grid, location) =
    view(f.parent, underlying_left_halo_indices(location, topology(grid, 1), grid.Nx, grid.Hx), :, :)

underlying_east_halo(f, grid, location) =
    view(f.parent, underlying_right_halo_indices(location, topology(grid, 1), grid.Nx, grid.Hx), :, :)

underlying_south_halo(f, grid, location) =
    view(f.parent, :, underlying_left_halo_indices(location, topology(grid, 2), grid.Ny, grid.Hy), :)

underlying_north_halo(f, grid, location) =
    view(f.parent, :, underlying_right_halo_indices(location, topology(grid, 2), grid.Ny, grid.Hy), :)

underlying_bottom_halo(f, grid, location) =
    view(f.parent, :, :, underlying_left_halo_indices(location, topology(grid, 3), grid.Nz, grid.Hz))

underlying_top_halo(f, grid, location) =
    view(f.parent, :, :, underlying_right_halo_indices(location, topology(grid, 3), grid.Nz, grid.Hz))

#####
##### Viewing boundary grid points (used to fill other halos)
#####

left_boundary_indices(loc, topo, N, H) = 1:H
left_boundary_indices(::Type{Nothing}, topo, N, H) = 1:0 # empty

right_boundary_indices(loc, topo, N, H) = N-H+1:N
right_boundary_indices(::Type{Face}, ::Type{Bounded}, N, H) = N-H:N+1
right_boundary_indices(::Type{Nothing}, topo, N, H) = 1:0 # empty

underlying_left_boundary_indices(loc, topo, N, H) = 1+H:2H
underlying_left_boundary_indices(::Type{Nothing}, topo, N, H) = 1:0 # empty

underlying_right_boundary_indices(loc, topo, N, H) = N+1:N+H
underlying_right_boundary_indices(::Type{Face}, ::Type{Bounded}, N, H) = N+2:N+H+1
underlying_right_boundary_indices(::Type{Nothing}, topo, N, H) = 1:0 # empty

underlying_west_boundary(f, grid, location) =
    view(f.parent, underlying_left_boundary_indices(location, topology(grid, 1), grid.Nx, grid.Hx), :, :)

underlying_east_boundary(f, grid, location) =
    view(f.parent, underlying_right_boundary_indices(location, topology(grid, 1), grid.Nx, grid.Hx), :, :)

underlying_south_boundary(f, grid, location) =
    view(f.parent, :, underlying_left_boundary_indices(location, topology(grid, 2), grid.Ny, grid.Hy), :)

underlying_north_boundary(f, grid, location) =
    view(f.parent, :, underlying_right_boundary_indices(location, topology(grid, 2), grid.Ny, grid.Hy), :)

underlying_bottom_boundary(f, grid, location) =
    view(f.parent, :, :, underlying_left_boundary_indices(location, topology(grid, 3), grid.Nz, grid.Hz))

underlying_top_boundary(f, grid, location) =
    view(f.parent, :, :, underlying_right_boundary_indices(location, topology(grid, 3), grid.Nz, grid.Hz))

#####
##### Convinience functions
#####

function sides_in_the_same_dimension(side1, side2)
    x_sides = (:west, :east)
    y_sides = (:south, :north)
    z_sides = (:bottom, :top)
    side1 in x_sides && side2 in x_sides && return true
    side1 in y_sides && side2 in y_sides && return true
    side1 in z_sides && side2 in z_sides && return true
    return false
end

function cubed_sphere_halo(cubed_sphere_field, location, face_index, side)
    LX, LY, LZ = location
    src_field = get_face(cubed_sphere_field, face_index)
    side == :west  && return  underlying_west_halo(src_field.data, src_field.grid, LX)
    side == :east  && return  underlying_east_halo(src_field.data, src_field.grid, LX)
    side == :south && return underlying_south_halo(src_field.data, src_field.grid, LY)
    side == :north && return underlying_north_halo(src_field.data, src_field.grid, LY)
end

function cubed_sphere_boundary(cubed_sphere_field, location, face_index, side)
    LX, LY, LZ = location
    src_field = get_face(cubed_sphere_field, face_index)
    side == :west  && return  underlying_west_boundary(src_field.data, src_field.grid, LX)
    side == :east  && return  underlying_east_boundary(src_field.data, src_field.grid, LX)
    side == :south && return underlying_south_boundary(src_field.data, src_field.grid, LY)
    side == :north && return underlying_north_boundary(src_field.data, src_field.grid, LY)
end


text: my_f(x,y) = 2x+y


text:
function contraction_graph(As::ITensor...)
  N = length(As)

  # Make edges for the contracted indices
  edge_index_list = Dict{Tuple{Int, Int}, Vector{Index}}()
  for nodeᵢ in 1:N
    Aᵢ = As[nodeᵢ]
    for nodeⱼ in nodeᵢ:N
      if nodeᵢ ≠ nodeⱼ
        Aⱼ = As[nodeⱼ]
        for indAⱼ in inds(Aⱼ)
          pos = findfirst(==(indAⱼ), inds(Aᵢ))
          if !isnothing(pos)
            indAᵢ = ind(Aᵢ, pos)
            @assert dir(indAᵢ) == -dir(indAⱼ)
            if indAⱼ in inds(Aᵢ)
              edge = (nodeᵢ, nodeⱼ)
              if dir(indAⱼ) == ITensors.Out
                edge = reverse(edge)
              else
                indAⱼ = dag(indAⱼ)
              end
              current_indsᵢⱼ = get(edge_index_list, edge, Index[])
              indsᵢⱼ = push!(current_indsᵢⱼ, indAⱼ)
              edge_index_list[edge] = indsᵢⱼ
            end
          end
        end
      end
    end
  end

  # Make nodes out of the uncontracted indices
  uncontracted_inds = noncommoninds(As...)
  for nodeᵢ in N+1:N+length(uncontracted_inds)
    uncontracted_indᵢ = uncontracted_inds[nodeᵢ-N]
    for nodeⱼ in 1:N
      Aⱼ = As[nodeⱼ]
      pos = findfirst(==(uncontracted_indᵢ), inds(Aⱼ))
      if !isnothing(pos)
        edge = (nodeᵢ, nodeⱼ)
        if dir(uncontracted_indᵢ) == ITensors.Out
          edge = reverse(edge)
        else
          uncontracted_indᵢ = dag(uncontracted_indᵢ)
        end
        edge_index_list[edge] = [uncontracted_indᵢ]
      end
    end
  end

  return edge_index_list
end

function get_adjacency_list(edge_index_list::Dict{Tuple{Int, Int}, Vector{Index}})
  # Determine the number of nodes from the edge list
  N = 1
  for edge in keys(edge_index_list)
    N = max(N, edge...)
  end

  # Make an adjacency list from the edge list (`adjlist[i]`: lists the nodes
  # that node `i` connects to)
  adjlist = Vector{Int}[Int[] for n in 1:N]
  for edge in keys(edge_index_list)
    for _ in 1:length(edge_index_list[edge])
      push!(adjlist[first(edge)], last(edge))
    end
  end
  return adjlist
end

function get_source_destiny_weight(source::AbstractVector, destiny::AbstractVector)
  if length(source) != length(destiny)
    throw(ArgumentError("Source and destiny must have the same length."))
  end
  source, destiny, Float64[ 1.0 for i in source ]
end

function get_source_destiny_weight(adjlist::AbstractVector{V}) where V<:AbstractVector{T} where T<:Any
  source = Int[]
  destiny = Int[]
  for (i,l) in enumerate(adjlist)
    for j in l
      push!(source, i)
      push!(destiny, j)
    end
  end
  return get_source_destiny_weight(source, destiny)
end

function get_adjacency_matrix(source::AbstractVector{Int}, destiny::AbstractVector{Int}, weights::AbstractVector)
  n = max(maximum(source), maximum(destiny))
  Matrix(sparse(source, destiny, weights, n, n))
end

# GraphRecipes.get_adjacency_matrix(adjlist)
function get_adjacency_matrix(adjlist::Vector{Vector{Int}})
  s, d, w = get_source_destiny_weight(adjlist)
  return get_adjacency_matrix(s, d, w)
end

function get_adjacency_matrix(edge_index_list::Dict{Tuple{Int, Int}, Vector{Index}})
  return get_adjacency_matrix(get_adjacency_list(edge_index_list))
end

max_stars_repo_name: mtfishman/ITensorsVisualization

max_issues_repo_name: mtfishman/ITensorsVisualization

max_forks_repo_name: mtfishman/ITensorsVisualization


text:
if PROGRAM_FILE == basename(@__FILE__)
    println("please  julia -i sevenstars.jl")
    exit()
end

module Front

export WelcomeController

using Bukdu # ApplicationController Conn
using ..Layout: layout
using ..WASM: WasmController
using Sockets: @ip_str
import Bukdu.Actions: index

struct WelcomeController <: ApplicationController
    conn::Conn
end

function get_server()
    if haskey(ENV, "ON_HEROKU")
        (port=parse(Int, ENV["PORT"]), host=Sockets.IPAddr(0,0,0,0))
    else
        (port=8080, host=ip"127.0.0.1")
    end
end

server = get_server()

end # module Front

import .Layout: layout
function Front.index(::Front.WelcomeController)
    wasm_path = Router.Helpers.url_path(get, WasmController, index)
    uri = "http://$(Front.server.host):$(Front.server.port)"
    title = ""
    script = ""
    style = ""
    body = """

Welcome

1. REST API Demo

GET     /customers           CustomerController  index   :api
GET     /customers/new       CustomerController  new     :api
GET     /customers/:id/edit  CustomerController  edit    :api
GET     /customers/:id       CustomerController  show    :api
POST    /customers           CustomerController  create  :api
DELETE  /customers/:id       CustomerController  delete  :api
PATCH   /customers/:id       CustomerController  update  :api
PUT     /customers/:id       CustomerController  update  :api

curl -X GET    $uri/customers
curl -X GET    $uri/customers/new
curl -X GET    $uri/customers/1/edit
curl -X GET    $uri/customers/1
curl -X POST   $uri/customers
curl -X DELETE $uri/customers/1
curl -X PATCH  $uri/customers/1
curl -X PUT    $uri/customers/1

2. WebAssembly Demo

Go => WASM

"""
    layout(title, script, style, body)
end


text:
using Test
using LinearAlgebra, StaticArrays
using BenchmarkTools
using RobotDynamics
using RobotDynamics: LinearModel
const RD = RobotDynamics
include("random_linear.jl")

## Continuous LTI
n,m = 10,5
dt = 0.1
t = 0.0
A,B = gencontrollable(n,m)
model = LinearModel(A,B)
@test model.use_static == true
@test RD.is_discrete(model) == false
@test RD.is_timevarying(model) == false
@test RD.is_affine(model) == false
@test RD.dims(model) == (n,m,n)

x,u = rand(model)
z = KnotPoint(x,u,dt,t)
@test length(x) == n
@test length(u) == m
@test dynamics(model, x, u) ≈ A*x + B*u
@test_throws AssertionError discrete_dynamics(RD.PassThrough, model, x, u, t, dt)
x2 = x + (A*(x + (A*x + B*u)*dt/2) + B*u)*dt
@test discrete_dynamics(RK2, model, z) ≈ x2

@test @ballocated(discrete_dynamics($model, $z), samples=2, evals=2) == 0
@test @ballocated(dynamics($model, $z), samples=2, evals=2) == 0

F = RD.DynamicsJacobian(n,m)
kp = KnotPoint(x,u,dt,t)
@test_throws AssertionError discrete_jacobian!(PassThrough, F, model, kp)
jacobian!(F, model, kp)
@test RD.get_A(F) == A
@test RD.get_B(F) == B

## Non-static model
model = LinearModel(A,B,use_static=false)
@test model.use_static == false
@test RD.is_discrete(model) == false
@test RD.is_timevarying(model) == false
@test RD.is_affine(model) == false
@test RD.dims(model) == (n,m,n)

@test dynamics(model, z) ≈ A*x + B*u
@test discrete_dynamics(RK2, model, z) ≈ x2

@test @ballocated(dynamics($model, $z), samples=2, evals=2) == 0
@test @ballocated(discrete_dynamics(RK2, $model, $z), samples=2, evals=2) == 0

# affine
d = randn(n)
model = LinearModel(A,B,d)
@test RD.is_affine(model) == true
@test RD.is_timevarying(model) == false
@test RD.is_discrete(model) == false
@test dynamics(model, z) ≈ A*x + B*u + d

# Discrete LTI
model = LinearModel(A,B,d;dt=dt)
@test RD.is_affine(model) == true
@test RD.is_timevarying(model) == false
@test RD.is_discrete(model) == true

@test_throws AssertionError dynamics(model, z)
@test discrete_dynamics(RD.PassThrough, model, x, u, t, dt) ≈ A*x + B*u + d

## Continuous LTV
N = 11
tf = dt*(N-1)
t = 0.5
AB = [gencontrollable(n,m) for k = 1:N]
A = [A for (A,B) in AB]
B = [B for (A,B) in AB]
times = range(0,tf, length=N)
k = findfirst(times .== t)

@test_throws AssertionError LinearModel(A,B)
model = LinearModel(A,B, times=times)
x,u = rand(model)
@test length(x) == n
@test length(u) == m
z = KnotPoint(x,u,dt,t)

@test dynamics(model, z) ≈ A[k]*x + B[k]*u
F = RD.DynamicsJacobian(model)
jacobian!(F, model, z)
@test F ≈ [A[k] B[k]]

@test RD.is_discrete(model) == false
@test RD.is_affine(model) == false
@test RD.is_timevarying(model) == true

x2 = x + (A[k]*(x + (A[k]*x + B[k]*u)*dt/2) + B[k]*u)*dt
@test discrete_dynamics(RK2, model, z) ≈ x2

## Affine Discrete LTV
d = [@SVector rand(n) for k = 1:N]
model = LinearModel(A,B,d, times=times, dt=dt)
@test_throws AssertionError dynamics(model, z)
@test discrete_dynamics(PassThrough, model, z) ≈ A[k]*x + B[k]*u + d[k]

discrete_jacobian!(PassThrough, F, model, z)
@test F ≈ [A[k] B[k]]

@test_throws AssertionError discrete_dynamics(RK2, model, z)

@test RD.is_discrete(model) == true
@test RD.is_affine(model) == true
@test RD.is_timevarying(model) == true

## Varying time steps
dts = normalize(rand(N-1), 1)*tf
times = insert!(cumsum(dts), 1, 0)
linmodel = LinearModel(A, B, times=times, dt=NaN)
@test RD.is_discrete(linmodel)

@test discrete_dynamics(PassThrough, linmodel, x, u, 0.0, dts[1]) ≈ A[1]*x + B[1]*u
@test discrete_dynamics(PassThrough, linmodel, x, u, times[2], dts[2]) ≈ A[2]*x + B[2]*u
@test_throws AssertionError discrete_dynamics(PassThrough, linmodel, x, u, 0.0, dts[2])


text:
let
    # Test first function here
    r1, r2 = your_function(3, 1.5)
    @test r1 == 4.5
    @test r2 == 1.5

    # Test corner case to function here
    r1, r2 = your_function(1, -42)
    @test r1 == 0
    @test r2 == 0
end

let
    # Here's an example of using the function array_isapprox
    array = [1.0 1.001 0.999]
    array_isapprox(array, 1.0, atol=0.002)
end


text:
# Issue #14
stack = @make_stack map(x->2x, 1:10)
state = dummy_state(stack)
execute_command(state, state.stack[1], Val{:finish}(), "finish")
@test isempty(state.stack)
@test state.overall_result == 2 .* [1:10...]

# Issue #12
function complicated_keyword_stuff(args...; kw...)
    args[1] == args[1]
    (args..., kw...)
end
stack = @make_stack complicated_keyword_stuff(1)
state = dummy_state(stack)
execute_command(state, state.stack[1], Val{:n}(), "n")
execute_command(state, state.stack[1], Val{:finish}(), "finish")
@test isempty(state.stack)