text:
using Test
using BBSPI

@testset "BBSPI" begin

    # Input pin that produces a test pattern
    mutable struct TestOutPin
        bit
        bits
    end
    Base.getindex(b::TestOutPin) = b.bit
    Base.length(b::TestOutPin) = 1

    # Output pins that reocrd the output value
    mutable struct TestInPin
        v
    end
    Base.getindex(c::TestInPin) = c.v

    function Base.setindex!(pin::TestInPin, v)

        if pin == clk
            # On rising clock edge, load a bit into MISO...
            if pin.v == 0 && v == 1
                if miso isa Vector
                    for p in miso
                        p.bit = popfirst!(p.bits)
                    end
                else
                    miso.bit = popfirst!(miso.bits)
                end
            end
#            println(  "clock = ", v,
#                    ", cs = ", cs.v,
#                    ", mosi = ", mosi.v,
#                    ", miso = ",  miso.bit)
        end

        # Record transmitted bits
        if pin == mosi
            push!(txbuf, v)
        end
    end

    rxbuf = zeros(UInt8, 4)

    BBSPI.delay(s::BBSPI.SPISlave) = nothing

    cs = TestInPin(0)
    clk = TestInPin(0)
    mosi = TestInPin(0)
    txbuf = Int[]
    miso = TestOutPin(-1, [0,1,0,1,0,1,0,1,
                           1,0,1,0,1,0,1,0,
                           0,1,0,1,0,1,0,1,
                           1,0,1,0,1,0,1,0])

    spi = BBSPI.SPISlave(cs=cs, clk=clk, mosi=mosi, miso=miso)

    BBSPI.transfer(spi, [0xAA, 0x55], rxbuf)
    @test rxbuf == [0x55, 0xaa, 0x55, 0xaa]
    @test txbuf == [1,0,1,0,1,0,1,0,
                    0,1,0,1,0,1,0,1,
                    0,0,0,0,0,0,0,0,
                    0,0,0,0,0,0,0,0]

    cs = TestInPin(0)
    clk = TestInPin(0)
    mosi = TestInPin(0)
    txbuf = Int[]
    miso1 = TestOutPin(-1, [0,1,0,1,0,1,0,1,
                            1,0,1,0,1,0,1,0,
                            0,1,0,1,0,1,0,1,
                            1,0,1,0,1,0,1,0])

    miso2 = TestOutPin(-1, [1,0,1,0,1,0,1,0,
                            0,1,0,1,0,1,0,1,
                            1,0,1,0,1,0,1,0,
                            0,1,0,1,0,1,0,1])

    miso = [miso1, miso2]

    spi = BBSPI.SPISlave(cs=cs, clk=clk, mosi=mosi, miso=miso)

    rxbuf = zeros(UInt8, 4, 2)

    BBSPI.transfer(spi, [0xAA, 0x55], rxbuf)

    @test rxbuf[:,1] == [0x55, 0xaa, 0x55, 0xaa]
    @test rxbuf[:,2] == [0xaa, 0x55, 0xaa, 0x55]
end


text:
using CSV
using Plots
using DataFrames

headers = ["timestamp", "cambrian", "info", "gen", "max", "avg", "std"]
df = CSV.read("path/to/csv", DataFrame, header=headers)

plot(df.gen, df.avg, xlabel="Generation", ylabel="Negative fitness", label="label")

ext: jl

max_stars_repo_name: mahmoud-al-najar/cgp-morphcast

max_issues_repo_name: mahmoud-al-najar/cgp-morphcast

max_forks_repo_name: mahmoud-al-najar/cgp-morphcast


text:
struct PhaseRecord
    phase_name::String
    obj::Num
    grad::Vector{Num}
    hess::Matrix{Num}
    mass::Vector{Num}
    mass_jac::Matrix{Num}
    state_variables::Vector{Num}
    site_fractions::Vector{Num}
    phase_matrix::Matrix{Num}
    inv_phase_matrix::Matrix{Num}
    ℵ::Num
end

function phase_matrix(hess, num_potentials, num_site_fractions, sublattice_dof)
    phase_mat = Matrix{Num}(undef, num_site_fractions + length(sublattice_dof), num_site_fractions + length(sublattice_dof))
    # Fill the upper-left corner with the Hessian terms
    for i in 1:num_site_fractions
        for j in 1:num_site_fractions
            phase_mat[i,j] = hess[num_potentials+i,num_potentials+j]
        end
    end
    # Fill the remaining rows (and the symmetric terms in the columns)
    # using the condition that the sum of site fractions is one
    subl_offset = 0
    for s in 1:length(sublattice_dof)
        # The first columns until the offset are zero
        phase_mat[num_site_fractions+s,1:subl_offset] .= 0.0
        # rows in symmetric column
        phase_mat[1:subl_offset,num_site_fractions+s] .= 0.0
        # The columns for the current sublattice are all one
        phase_mat[num_site_fractions+s,subl_offset+1:subl_offset+sublattice_dof[s]] .= 1.0
        # rows in symmetric column
        phase_mat[subl_offset+1:subl_offset+sublattice_dof[s],num_site_fractions+s] .= 1.0
        # Remaining columns are zero
        phase_mat[num_site_fractions+s,subl_offset+sublattice_dof[s]+1:end] .= 0.0
        # rows in symmetric columns,
        # some repeated work here, but it's consistent and logically equal
        phase_mat[subl_offset+sublattice_dof[s]+1:end,num_site_fractions+s] .= 0.0
        subl_offset += sublattice_dof[s]
    end
    return phase_mat
end

function PhaseRecord(phase_name, obj, mass, potentials, site_fractions, sublattice_dof, unique_phase_name=nothing)
    ℵ = Num(Variable(Symbol("ℵ_$(something(unique_phase_name, phase_name))")))
    grad = Symbolics.gradient(obj, [potentials..., site_fractions...])
    hess = Symbolics.hessian(obj, [potentials..., site_fractions...])
    mass_jac = Symbolics.jacobian(mass, [potentials..., site_fractions...])
    pm = phase_matrix(hess, length(potentials), length(site_fractions), sublattice_dof)

    return PhaseRecord(phase_name,
                       obj,
                       grad,
                       hess,
                       mass,
                       mass_jac,
                       potentials,
                       site_fractions,
                       pm,
                       LinearAlgebra.inv(pm),
                       ℵ,
                       )
end

mutable struct CompSet
    phase_rec::PhaseRecord
    Y::Vector{Float64}
    ℵ::Float64
end

function Base.show(io::IO, composition_set::CompSet)
    print(string(typeof(composition_set)))
    print("($(composition_set.phase_rec.phase_name), $(composition_set.Y), $(composition_set.ℵ))")
end


text:
# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
  Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using Documenter, ImageQuilting

makedocs(
  format = :html,
  sitename = "ImageQuilting.jl",
  authors = "Júlio Hoffimann Mendes",
  pages = [
    "Home" => "index.md",
    "Concepts" => "concepts.md",
    "Examples" => "examples.md",
    "Voxel reuse" => "voxel-reuse.md",
    "GPU support" => "gpu-support.md",
    "About" => [
      "Author" => "about/author.md",
      "License" => "about/license.md",
      "Citation" => "about/citation.md"
    ]
  ]
)

deploydocs(
  repo  = "github.com/juliohm/ImageQuilting.jl.git",
  target = "build",
  deps = nothing,
  make = nothing,
  julia = "1.0"
)


text:
abstract StochasticApproximationOptimizer <: Optimizer

SPSADefaultParameters = {
  :Alpha => 0.602,  # The optimal value is 1.0 but values down to 0.602 often can give faster convergence
  :Gamma => 0.101,  # The optimal value is 1/6 but values down to 0.101 often can give faster convergence
  :a     => 0.0017,
  :c     => 1.9, # Recommendation is value 1.9 but that assumes noisy function, otherwise should be low
  :A     => 10
}

type SimultaneousPerturbationSA2 <: StochasticApproximationOptimizer
  search_space::SearchSpace
  parameters::Parameters
  k::Int64
  n::Int64
  theta::Array{Float64, 2}
  delta_ck::Array{Float64, 2}

  SimultaneousPerturbationSA2(parameters) = begin
    ss = parameters[:SearchSpace]
    n = numdims(ss)
    new(ss, Parameters(parameters, SPSADefaultParameters),
      0, n, rand_individual(ss), zeros(Float64, n, 2))
  end
end

name(spsa::SimultaneousPerturbationSA2) = "SPSA2 (Simultaneous Perturbation Stochastic Approximation, 1st order, 2 samples)"

sample_bernoulli_vector(n::Int64) = 2.0 * round(rand(n,1)) - 1.0

function ask(spsa::SimultaneousPerturbationSA2)

  delta = sample_bernoulli_vector(spsa.n)
  ck = spsa.parameters[:c]/(spsa.k + 1)^spsa.parameters[:Gamma]
  spsa.delta_ck = ck * delta

  theta_plus = spsa.theta + spsa.delta_ck
  theta_minus = spsa.theta - spsa.delta_ck

  [(theta_plus, 1), (theta_minus, 2)]

end

function tell!(spsa::SimultaneousPerturbationSA2, rankedCandidates)

  # Use index of rank to get right values for yplus and yminus, respectively.
  if rankedCandidates[1][2] == 1
    yplus = rankedCandidates[1][3]
    yminus = rankedCandidates[2][3]
  else
    yplus = rankedCandidates[2][3]
    yminus = rankedCandidates[1][3]
  end

  # Estimate gradient.
  ghat = (yplus - yminus) ./ (2.0 * spsa.delta_ck)

  # Calc new estimate of theta based on estimate of gradient, ghat.
  ak = spsa.parameters[:a]/(spsa.k + 1 + spsa.parameters[:A])^spsa.parameters[:Alpha]
  theta_new = spsa.theta - ak * ghat
  rand_bound_from_target!(theta_new, spsa.theta, spsa.search_space)
  spsa.theta = theta_new
  spsa.k += 1

  return 0
end


text:
module Revealables

using Interact
using Reactive

import Base.writemime, Base.display

export Revealable
export revealable
export encode!
export encode
export decode!
export decode

type Revealable
    markdown::ASCIIString
    label::ASCIIString
    show::Bool
    encoded::Bool
end

Revealable(markdown::ASCIIString, label::ASCIIString = "Show/Hide", show = false, encoded = false) = Revealable(markdown, label, show, encoded)

function revealable(markdown::ASCIIString, label::ASCIIString, show::Bool)
    x = Revealable(markdown, label, show)
    revealable(x)
end

function revealable(markdown::ASCIIString, label::ASCIIString = "Show/Hide")
    x = Revealable(markdown, label, false)
    revealable(x)
end

function revealable(x::Revealable)
    @manipulate for n in togglebutton(; label = x.label, value=x.show, signal=Input(x.show))
        x.show = n
        x
    end
end

function rot(text::ASCIIString, key::Int)
    map(text) do c
        if 'a' <= c <= 'z'
            char(mod(c - 'a' + key, 26) + 'a')
        elseif 'A' <= c <= 'Z'
            char(mod(c - 'A' + key, 26) + 'A')
        elseif '&' <= c <= '?'
            char(mod(c - '&' + key, 26) + '&')
        else
            c
        end
    end
end

function encode!(x::Revealable, password::ASCIIString)
    if x.encoded
        x
    else
        key = sum(password)
        x.markdown = rot(x.markdown, key)
        x.encoded = true
        x
    end
end

function encode(x::Revealable, password::ASCIIString)
    r = Revealable(x.markdown, x.label, x.show)
    encode!(r, password)
end

function decode!(x::Revealable, password::ASCIIString)
    if x.encoded
        key = 26 - sum(password)
        x.markdown = x.markdown = rot(x.markdown, key)
        x.encoded = false
        x
    else
        x
    end
end

function decode(x::Revealable, password::ASCIIString)
    r = Revealable(x.markdown, x.label, x.show)
    decode!(r, password)
end

function writemime(io::IO, ::MIME"text/markdown", x::Revealable)
    if x.show
        println(io, x.markdown)
    else
        println(io, "")
    end
end

function display(io::IO, x::Revealable)
    if x.show
        x
    else
        println(io, "")
    end
end

end # module


text:
import Base: show, length, collect
import RigidBodyTools: RigidTransform

export PotentialFlowBody, subtractcirculation!, getU, getΩ, getΓ, setU, setΩ, setΓ

"""
$TYPEDEF

Defines a rigid body for `VortexModel` in `GridPotentialFlow.jl`.
"""
mutable struct PotentialFlowBody
    """body: """
    points::Body
    """U: Linear velocity of the body."""
    U::Tuple{Float64,Float64}
    """Ω: Rotational velocity of the body."""
    Ω::Float64
    """Γ: Current circulation about the body."""
    Γ::Float64
    """regidx: Array of the indices of any regularized edges."""
    edges::Vector{Int}
    """σ: Array of suction parameters or suction parameter ranges associated with the entries of `edges`."""
    σ::Vector{Union{SuctionParameter,SuctionParameterRange}}
end

"""
$(TYPEDSIGNATURES)

Constructs a potential flow body with shape `b`, linear velocity `U`, rotational velocity `Ω`, and initial circulation `Γ`. Edges that should be regularized can be specified in `edges` and their suction parameter ranges in `σ`, which defaults to zero for each regularized edge.
"""
function PotentialFlowBody(b::Body{N,C}; U = (0.0,0.0), Ω = 0.0, Γ = 0.0, edges = Int64[], σ = SuctionParameterRange[]) where {N,C}

    Ne = length(edges)

    if Ne > 0 && isempty(σ)
        σ = [SuctionParameterRange(0.0,0.0) for e in 1:Ne]
    elseif Ne > length(σ)
        println("Not enough suction parameters provided. Setting all suction parameters to zero.")
    end

    return PotentialFlowBody(b,U,Ω,Γ,edges,σ)
end

"""
$(TYPEDSIGNATURES)

Returns the number of surface points in the body `b`.
"""
function Base.length(b::PotentialFlowBody)
    return length(b.points)
end

"""
$(TYPEDSIGNATURES)

Returns the subrange of indices in the global set of surface point data corresponding to body `i` in `b`.
"""
function getrange(b::AbstractVector{PotentialFlowBody},i::Int)
    first = 1
    j = 1
    while j < i
        first += length(b[j])
        j += 1
    end
    last = first+length(b[i])-1
    return first:last
end

"""
$(TYPEDSIGNATURES)

Returns the linear velocity of the body `b`.
"""
function getU(b::PotentialFlowBody)
    return b.U
end

"""
$(TYPEDSIGNATURES)

Returns the rotational velocity of the body `b`.
"""
function getΩ(b::PotentialFlowBody)
    return b.Ω
end

"""
$(TYPEDSIGNATURES)

Returns the bound circulation of the body `b`.
"""
function getΓ(b::PotentialFlowBody)
    return b.Γ
end

"""
$(TYPEDSIGNATURES)

Sets the linear velocity of the body `b` to `U`.
"""
function setU(b::PotentialFlowBody, U::Tuple{Float64,Float64})
    b.U = U
end

"""
$(TYPEDSIGNATURES)

Sets the rotational velocity of the body `b` to `Ω`.
"""
function setΩ(b::PotentialFlowBody, Ω::Float64)
    b.Ω = Ω
end

"""
$(TYPEDSIGNATURES)

Sets the bound circulation of the body `b` to `Γ`.
"""
function setΓ(b::PotentialFlowBody, Γ::Float64)
    b.Γ = Γ
end

"""
$(TYPEDSIGNATURES)

Returns the indices in the global set of surface point data of the regularized points of body `i` in `b`.
"""
function getregularizededges(b::AbstractVector{PotentialFlowBody},i::Int)
    r = getrange(b,i)
    idx = deepcopy(b[i].edges) .+ r[1] .- 1
    return idx
end

"""
$(TYPEDSIGNATURES)

Returns the indices in the global set of surface point data of all regularized points in `b`.
"""
function getregularizededges(b::AbstractVector{PotentialFlowBody})
    idx = Int64[]
    for i in 1:length(b)
        append!(idx,getregularizededges(b,i))
    end
    return idx
end

"""
$(TYPEDSIGNATURES)

Returns the indices in the global set of surface point data of all regularized points in `b`.
"""
function subtractcirculation!(b::AbstractVector{PotentialFlowBody}, δΓ_vec::AbstractVector)
    i = 0
    for j in 1:length(b)
        for k in 1:length(b[j].edges)
            i += 1
            b[j].Γ -= δΓ_vec[i]
        end
    end
end

function bodypointsvelocity!(v::AbstractVector, b::AbstractVector{PotentialFlowBody}, dir::Int)
    for i in 1:length(b)
        vi = view(v,getrange(b,i))
        vi .= b[i].U[dir]
    end
end

function Base.collect(b::AbstractVector{PotentialFlowBody})
    x = Vector{Float64}()
    y = Vector{Float64}()
    for i in 1:length(b)
        xi, yi = RigidBodyTools.collect(b[i].points)
        append!(x,deepcopy(xi))
        append!(y,deepcopy(yi))
    end

    return x, y
end

function (T::RigidBodyTools.RigidTransform)(b::PotentialFlowBody)
    T(b.points)
    return b
end

function RigidBodyTools.dlength(b::PotentialFlowBody)
    dlength(b.points)
end

function RigidBodyTools.dlengthmid(b::PotentialFlowBody)
    dlengthmid(b.points)
end

function Base.show(io::IO, b::PotentialFlowBody)
    Ne = length(b.edges)
    iobody = IOBuffer()
    Base.show(iobody, b.points)
    sbody = String(take!(iobody))
    println(io,"Potential flow body with $Ne regularized $(Ne == 1 ? "edge" : "edges")")
    println(io,"Shape and position: $(sbody[1:end-1])")
    println(io,"Linear velocity: $(b.U)")
    println(io,"Angular velocity: $(b.Ω)")
    println(io,"Current circulation: $(b.Γ)")
    if Ne > 0
        println(io,"Suction parameter limits:")
        for e in 1:Ne
            println(io,"   Edge $e: $(b.σ[e])")
        end
    end
end

@recipe f(::Type{PotentialFlowBody}, pfb::PotentialFlowBody) = pfb.points


text:
using Test
using Zygote, CUDA
import SpecialFunctions
import StatsFuns
import CUDAExtensions

CUDA.allowscalar(false)

function test_ad(f, args...)
    y, ∂_cpu = Zygote.pullback(f, args...)
    _, ∂_gpu = Zygote.pullback(f, map(cu, args)...)

    ∇_cpu = ∂_cpu(one(y))
    ∇_gpu = ∂_gpu(one(y))

    # Convert back to CPU
    ∇_gpu_cpu = map(Array, ∇_gpu)

    for i = 1:length(∇_cpu)
        # Ensure that the values align
        @test ∇_cpu[i] ≈ ∇_gpu_cpu[i] atol=1e-3

        # Ensure that eltype is preserved
        @test eltype(∇_gpu_cpu[i]) == eltype(args[i])
    end
end

test_ad_broadcast(f, args...) = test_ad(args...) do (inner_args...)
    sum(f.(inner_args...))
end

@testset "SpecialFunctions.jl" begin
    x = Float32.(randn(10))

    test_ad_broadcast(SpecialFunctions.erf, x)
    test_ad_broadcast(SpecialFunctions.erfinv, SpecialFunctions.erf.(x))
    test_ad_broadcast(SpecialFunctions.erfc, x)
    test_ad_broadcast(SpecialFunctions.erfcinv, SpecialFunctions.erfc.(x))
    # test_ad_broadcast(SpecialFunctions.erfi, x) # Apparently doesn't have a CUDA equivalent
    test_ad_broadcast(SpecialFunctions.erfcx, x)
end


text:
using Crayons.Box
using DataStructures
using NumericalDataManipulation.CommonMath
using NumericalDataManipulation.Common

function piecewise_function(
  definitions::OrderedDict{Tuple{Float64, Float64}, Tuple{Function, Bool}}; ϵᵇᵖ::Float64=1e-2)
  @assert all(d -> d.first[1] < d.first[2], definitions)
  keypoints = flatten(keys(definitions))
  @assert issorted(keypoints)
  inf_point_pos = sort(findall(keypoint -> abs(keypoint) == Inf, keypoints))
  @assert length(inf_point_pos) ≤ 2
  @assert all(pos -> pos ∈ [1, length(keypoints)], inf_point_pos)
  @assert ifelse(length(inf_point_pos) == 2, keypoints[1] == -Inf, true)
  @assert ifelse(length(inf_point_pos) == 2, keypoints[end] == Inf, true)
  @assert ifelse(length(inf_point_pos) == 1 && inf_point_pos[1] == 1, keypoints[1] == -Inf, true)
  @assert ifelse(length(inf_point_pos) == 1 && inf_point_pos[1] == length(keypoints), keypoints[end] == Inf, true)
  @assert ifelse(length(inf_point_pos) == 0, abs(keypoints[1]) ≠ Inf, true)
  @assert ifelse(length(inf_point_pos) == 0, abs(keypoints[end]) ≠ Inf, true)

  lfuncs = map(d -> d.first[2] => d.second, collect(definitions)[1:end-1])
  rfuncs = map(d -> d.first[1] => d.second, collect(definitions)[2:end])
  @assert length(lfuncs) == length(rfuncs)
  Nₚ = length(lfuncs)
  point_pair_funcs = Vector{Pair{Tuple{Float64,Float64}, Tuple{Tuple{Function, Bool},Tuple{Function, Bool}}}}()
  for k = 1:Nₚ
    ppf = (lfuncs[k].first, rfuncs[k].first) => (lfuncs[k].second, rfuncs[k].second)
    push!(point_pair_funcs, ppf)
  end
  @assert all(ppf -> abs(ppf.first[1]) ≠ Inf, point_pair_funcs)
  @assert all(ppf -> abs(ppf.first[2]) ≠ Inf, point_pair_funcs)
  @assert all(ppf -> abs(ppf.first[1] - ppf.first[2]) ≤ ϵᵇᵖ, point_pair_funcs)
  breakpoint_funcs = map(ppf -> (ppf.first[1] + ppf.first[2]) / 2 => ppf.second, point_pair_funcs)
  breakpoint_func_gaps = map(bf -> bf.first => (bf.second[2][1](bf.first + ϵᵇᵖ) - bf.second[1][1](bf.first - ϵᵇᵖ), bf.second[1][2], bf.second[2][2]), breakpoint_funcs)
  @assert length(breakpoint_funcs) == length(breakpoint_func_gaps)
  breakpoint_parameters = map(
  bp_entry -> begin
      k = bp_entry[1]
      func_pars = bp_entry[2]
      func_gap = func_pars.second[1]
      keep_left = func_pars.second[2]
      keep_right = func_pars.second[3]
      Δy = func_gap; Δx = abs(𝐆 * Δy)
      Δx = Δx ≥ ϵᵇᵖ ? Δx : ϵᵇᵖ
      x₀ = func_pars.first
      f1 = breakpoint_funcs[k].second[1][1]
      f2 = breakpoint_funcs[k].second[2][1]
      # @assert typeof(f1) == Function
      # @assert typeof(f2) == Function
      return (x₀, Δx, Δy, (f1, f2), keep_left, keep_right)
    end,
    enumerate(breakpoint_func_gaps))
  return make_piecewise_sigmoid_function(definitions, breakpoint_funcs, breakpoint_parameters)
end

function make_piecewise_sigmoid_function(
  definitions::OrderedDict{Tuple{Float64, Float64},Tuple{Function, Bool}},
  breakpoint_funcs, breakpoint_parameters;
  ϵᵇᵖ::Float64=1e-2)
  @assert length(definitions) ≥ 2

  breakpoints = Vector{Float64}(undef, 0)
  functions = Vector{Function}(undef, 0)
  definition_list = collect(definitions)
  for k = 1:length(definition_list)-1
    defₖ_left = definition_list[k]
    xₖ_a_left = defₖ_left.first[1]; xₖ_b_left = defₖ_left.first[2]
    @assert xₖ_a_left < xₖ_b_left
    @assert ifelse(xₖ_a_left == -Inf, k == 1, true)
    keep_leftₖ = defₖ_left.second[2]
    Δxₖ_left = xₖ_b_left - xₖ_a_left
    if keep_leftₖ
      @info LIGHT_RED_FG("Interval [$xₖ_a_left, $xₖ_b_left] will be kept (untouched function).")
    end

    defₖ_right = definition_list[k + 1]
    xₖ_a_right = defₖ_right.first[1]; xₖ_b_right = defₖ_right.first[2]
    @assert xₖ_a_right < xₖ_b_right
    @assert ifelse(xₖ_b_right == Inf, k + 1 == length(definition_list), true)
    keep_rightₖ = defₖ_right.second[2]
    Δxₖ_right = xₖ_b_right - xₖ_a_right
    if keep_rightₖ
      @info LIGHT_RED_FG("Interval [$xₖ_a_right, $xₖ_b_right] will be kept (untouched function).")
    end

    parametersₖ = breakpoint_parameters[k]
    x₀ₖ = parametersₖ[1]
    Δxₛ = parametersₖ[2]
    Δyₛ = parametersₖ[3]
    f1, f2 = parametersₖ[4]
    kf1 = parametersₖ[5]
    kf2 = parametersₖ[6]
    @assert kf1 == keep_leftₖ
    @assert kf2 == keep_rightₖ

    δxₖ_left = ifelse(xₖ_a_left == -Inf,
      10 * 𝐆 * Δxₛ,
      ifelse(𝐆⁻¹ * Δxₛ < Δxₖ_left, 𝐆⁻¹ * Δxₛ/π, 𝐆1⁻¹ * Δxₖ_left / 2))

    δxₖ_right = ifelse(xₖ_b_right == Inf,
      10 * 𝐆 * Δxₛ,
      ifelse(𝐆⁻¹ * Δxₛ < Δxₖ_right, 𝐆⁻¹ * Δxₛ/π, 𝐆1⁻¹ * Δxₖ_right / 2))
    #@assert(δxₖ_left ≥ 𝐆⁻¹ * Δxₛ / 2, "δxₖ_left ≥ 𝐆⁻¹Δxₛ/2: δxₖ_left = $(δxₖ_left), 𝐆⁻¹Δxₛ/2 = $(𝐆⁻¹*Δxₛ/2)")
    #@assert(δxₖ_right ≥ 𝐆⁻¹ * Δxₛ / 2, "δxₖ_right ≥ 𝐆⁻¹Δxₛ/2: δxₖ_right = $(δxₖ_right),  𝐆⁻¹Δxₛ/2 = $(𝐆⁻¹*Δxₛ/2)")

    α = 0.2 * (δxₖ_left + δxₖ_right)
    δxₖ_left_new = keep_leftₖ ? δxₖ_left / 100.0 : δxₖ_left
    δxₖ_right_new = keep_rightₖ ? δxₖ_right / 100.0 : δxₖ_right
    if δxₖ_left_new < δxₖ_left
      x₀ₖ = x₀ₖ + δxₖ_left - δxₖ_left_new
    end
    if δxₖ_right_new < δxₖ_right
      x₀ₖ = x₀ₖ - δxₖ_right + δxₖ_right_new
    end
    @info "Blending function $f1; $f2; $x₀ₖ; $α"
    blending_functionₖ = sigmoid_of_name(f1, f2, x₀ₖ, α, "")

    breakpoint_left = x₀ₖ - δxₖ_left_new
    breakpoint_right = x₀ₖ + δxₖ_right_new
    function_left = defₖ_left.second[1]
    function_middle = blending_functionₖ
    function_right = defₖ_right.second[1]

    @info "Made blending function in interval:
    [$(breakpoint_left), $(breakpoint_right)],
    α = $α, x₀ₖ = $x₀ₖ, δxₖ_left = $δxₖ_left_new, δxₖ_right = $δxₖ_right_new,
    keep_leftₖ = $keep_leftₖ, keep_rightₖ = $keep_rightₖ,
    δxₖ_left_old = $δxₖ_left, δxₖ_right_old = $δxₖ_right"

    push!(functions, function_left)
    push!(breakpoints, breakpoint_left)
    push!(functions, function_middle)
    push!(breakpoints, breakpoint_right)
    if k + 1 == length(definition_list)
      push!(functions, function_right)
    end
  end
  pwf = make_piecewise_function(breakpoints, functions)
  return pwf, breakpoints
end

function make_piecewise_function(
  breakpoints::Vector{Float64}, functions::Vector{Function})
  pwf(R) = begin
    nᶠ = searchsortedfirst(breakpoints, R)
    f_piecewise = functions[nᶠ]
    return f_piecewise(R)
  end
end

const sigmoid = (x, x₀, α) -> 1/(1 + exp(-(x - x₀)/α))
function sigmoid_of_name(f1::Function, f2::Function, x₀, α, name::AbstractString)
  return x -> begin
    sf = (1 - sigmoid(x, x₀, α))*f1(x) + sigmoid(x, x₀, α)*f2(x)
    if sf === NaN || sf === NaN64 || sf === NaN32 || sf === NaN16
      @error "SIGMOID NaN: x=$x, x₀=$x₀, α=$α, f1(x)=$(f1(x)), f2(x)=$(f2(x)), σ(x)=$(sigmoid(x, x₀, α)), name=$name"
    end
    return sf
  end
end


text:
module Collision_3Elements

using  Modia3D
import Modia3D.ModiaMath

vmat1 = Modia3D.Material(color="Red", transparency=0.5)
vmat2 = Modia3D.Material(color="Blue", transparency=0.5)
vmat3 = Modia3D.Material(color="Green", transparency=0.5)
cmat  = "Steel"

@signal Sine(;A=1.0, freqHz = 1.0) begin
   y = ModiaMath.RealScalar(causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::Sine, sim::ModiaMath.SimulationState)
    signal.y.value  = signal.A*sin(2*pi*signal.freqHz*sim.time)
end

@assembly CollisionWithoutJoint(;length=1.0, width=0.2) begin
   world = Modia3D.Object3D(Modia3D.CoordinateSystem(0.6))
   #world_f1 = Modia3D.Object3D(world , r=[0.5, 0.0, 0.5])

   red = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), "Aluminium", vmat1;contactMaterial = cmat),
                        [ [-length/2,0,0],
                          [ length/2,0,0] ])
   blue = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), "Aluminium", vmat2;contactMaterial = cmat),
                        [ [-length/2,0,0],
                          [ length/2,0,0] ])
   green = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), "Aluminium", vmat3;contactMaterial = cmat),
                        [ [-length/2,0,0] ])

   Modia3D.connect(world, red.frames[1])
   rev2 = Modia3D.Revolute(red.frames[2], blue.frames[1]) #; canCollide = true)
   rev3 = Modia3D.Revolute(blue.frames[2], green.frames[1]) #; canCollide = true)

   sig1 = Sine(A=1.5,freqHz=0.5)
   sig2 = Sine(A=2.5,freqHz=1.0)
   sig3 = Sine(A=3.5,freqHz=1.5)

   sine2     = Modia3D.SignalToFlangeAngle(sig1.y)
   sine3     = Modia3D.SignalToFlangeAngle(sig1.y)

   Modia3D.connect(sine2,   rev2)
   Modia3D.connect(sine3,   rev3)
end

@assembly CollisionWithJoint(;length=1.0, width=0.2) begin
   world = Modia3D.Object3D(Modia3D.CoordinateSystem(0.6))
   world_f1 = Modia3D.Object3D(world , r=[0.5, 0.0, 0.5])

   red = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), "Aluminium", vmat1;contactMaterial = cmat),
                        [ [-length/2,0,0],
                          [ length/2,0,0] ])
   #blue = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), "Aluminium", vmat2;contactMaterial = cmat),
   #                     [ [-length/2,0,0],
   #                       [ length/2,0,0] ])
   blue = Modia3D.Part(Modia3D.Box(length,width,width, material=vmat2),
                                               [ [-length/2,0,0],
                                                 [ length/2,0,0] ])

   green = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), "Aluminium", vmat3;contactMaterial = cmat),
                        [ [-length/2,0,0] ])

   rev1 = Modia3D.Revolute(world_f1     , red.frames[1]) #; canCollide = true)
   rev2 = Modia3D.Revolute(red.frames[2], blue.frames[1]) #; canCollide = true)
   rev3 = Modia3D.Revolute(blue.frames[2], green.frames[1]) #; canCollide = true)

   sig1 = Sine(A=1.5,freqHz=0.5)
   sig2 = Sine(A=2.5,freqHz=1.0)
   sig3 = Sine(A=3.5,freqHz=1.5)

   sine1     = Modia3D.SignalToFlangeAngle(sig1.y)
   sine2     = Modia3D.SignalToFlangeAngle(sig1.y)
   sine3     = Modia3D.SignalToFlangeAngle(sig1.y)

   Modia3D.connect(sine1,   rev1)
   Modia3D.connect(sine2,   rev2)
   Modia3D.connect(sine3,   rev3)
end

collWithoutJoint=CollisionWithoutJoint(sceneOptions=Modia3D.SceneOptions(visualizeFrames=true,defaultFrameLength=0.3,enableContactDetection=true))

collWithJoint=CollisionWithJoint(sceneOptions=Modia3D.SceneOptions(visualizeFrames=true,defaultFrameLength=0.3,enableContactDetection=true))

#model = Modia3D.SimulationModel(collWithoutJoint, analysis=ModiaMath.KinematicAnalysis )
model = Modia3D.SimulationModel(collWithJoint, analysis=ModiaMath.KinematicAnalysis )
result = ModiaMath.simulate!(model, stopTime=0.01)

ModiaMath.plot(result, ("rev1.phi", "rev2.phi") )

println("... success of Collision_3Elements.jl!")
end


text:
# this is the 2nd attempt; in this attempt I will try use MCTS and gradients
using StatsBase, DataFrames, Flux, JLD2, Plots
include("2048_fn.jl")

model = Flux.Chain(
  Dense(14*16+4, 128, relu),
  Dense(128, 64, relu),
  Dense(64, 32, relu),
  Dense(32, 1, relu));

loss(x,y) = sum((model(x) .- y).^2)
opt = ADAM(params(model))

# play one game to see how it goes
function assess_model(model)
  init_grid = init_game()
  @time (init_grid, fnl, seq, cart, tf, reward) = simulate_game!(init_grid, model, 0.0)

  # obtain the discounted rewards
  y = discounted_rewards(reward, 0.99)[1:end-1];

  # score every move using  model
  updated_model = [Tracker.data(model(grid_move_to_feature(sim_seq(init_grid, seq[1:i+1], cart[1:i+1], tf[1:i+1]), seq[i])))[1] for i=1:length(seq)-1]
  display(sim_seq(init_grid, seq, cart, tf))
  g = plot(updated_model, label="model")
  plot!(y, label="discounted rewards")
  (init_grid, fnl, seq, cart, tf, reward, g)
end

ϵ = 0.5

model = nothing
@load "no_neg_model.jld2"
model

@time for i=1:10000
  global ϵ
  ϵ -= 0.0001
  simulate_one_game_and_train!(model, opt, ϵ)
  if i % 100 == 0
    @save "no_neg_model.jld2" model;
    println("epsilon: $ϵ")
    init_grid = init_game()
    res = simulate_game!(init_grid, model, 0)
    display(res[2])
    @save "no_neg_model.jld2" model;
  end
end

@time (init, fnl, seq, cart, tf, reward, g) = assess_model(model);
g

last_n_move = length(seq)-20
grid_new = sim_seq(init, seq[1:end-last_n_move], cart[1:end-last_n_move], tf[1:end - last_n_move])
grid_try = copy(grid_new)
display(grid_new)
vm = valid_moves(grid_new)
x = vec(Tracker.data(model(hcat(grid_move_to_feature.([grid_new], vm)...))))
y = x .- minimum(x)

using DataFrames
DataFrame(vm=vm, chosen=vm.== seq[end-last_n_move], x=vec(x),  prob_exp = round.(exp.(y)/sum(exp.(y)),digits=2), prob_square =((y.+1).^2)./sum((y.+1).^2))


text:
using BinaryBuilder, Pkg
using Base.BinaryPlatforms: arch, os

include("../../fancy_toys.jl")

name = "NCCL"
version = v"2.11.4"

script = raw"""
mkdir -p ${libdir} ${prefix}/include

cd ${WORKSPACE}/srcdir/nccl*
find .

install_license LICENSE.txt

mv lib/libnccl*.so* ${libdir}
mv include/* ${prefix}/include
"""

products = [
    LibraryProduct("libnccl", :libnccl, dont_dlopen = true),
]

# XXX: CUDA_loader_jll's CUDA tag should match the library's CUDA version compatibility.
#      lacking that, we can't currently dlopen the library

dependencies = [Dependency(PackageSpec(name="CUDA_loader_jll"))]

cuda_versions = [v"10.2", v"11.0", v"11.4"]
for cuda_version in cuda_versions
    cuda_tag = "$(cuda_version.major).$(cuda_version.minor)"
    include("build_$(cuda_tag).jl")

    for (platform, sources) in platforms_and_sources
        augmented_platform = Platform(arch(platform), os(platform); cuda=cuda_tag)
        should_build_platform(triplet(augmented_platform)) || continue
        build_tarballs(ARGS, name, version, sources, script, [augmented_platform],
                       products, dependencies; lazy_artifacts=true)
    end
end

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module ArrayLisp

using ReplMaker: initrepl

export expand, parse, SExpr, enable_repl, include_lisp, format, format_file
export @sexpr_str

const KEYWORDS = Set([
    :function,
    :parameters,
    :kw,
    :(=),
    :call,
    :macrocall,
    :block,
    :begin,
    :vect,
    :vcat,
    :ref,
    :curly,
    :using,
    :import,
    :(.),
    :for,
    :if,
    :elseif,
    :else,
    :(&&),
    :(||),
])

struct Atom
    val::Any
end

function Base.show(io::IO, a::Atom)
    write(io, string(a.val))
end

function Base.parse(::Type{Atom}, x::AbstractString)
    val = something(
        tryparse(Int, x),
        tryparse(Float64, x),
        x == "true" ? true : nothing,
        x == "false" ? false : nothing,
        first(x) == '"' && last(x) == '"' ? strip1(x, '"') : nothing,
        first(x) == '\'' && last(x) == '\'' && length(x) == 3 ? strip1(x, '\'') :
            nothing,
        Symbol(x),
    )
    return Atom(val)
end

function strip1(x, char)
    left_ind = 1
    right_ind = lastindex(x)
    if x[left_ind] == char
        left_ind = nextind(x, left_ind)
    end
    if x[right_ind] == char
        right_ind = prevind(x, right_ind)
    end
    return x[left_ind:right_ind]
end

struct SExpr
    terms::Any
end
Base.length(s::SExpr) = length(s.terms)
Base.push!(s::SExpr, x) = push!(s.terms, x)
Base.getindex(s::SExpr, i) = s.terms[i]
Base.lastindex(s::SExpr) = lastindex(s.terms)
Base.iterate(s::SExpr) = iterate(s.terms)
Base.iterate(s::SExpr, i::Int) = iterate(s.terms, i)

function Base.show(io::IO, s::SExpr)
    write(io, "(")
    for (i, e) in enumerate(s.terms)
        show(io, e)
        i < length(s) && write(io, " ")
    end
    write(io, ")")
end

const DEFAULT_FORMAT_LINELENGTH = 92
const DEFAULT_FORMAT_INDENT = 0

function format(
    io::IO,
    s::SExpr;
    linelength = DEFAULT_FORMAT_LINELENGTH,
    indent = DEFAULT_FORMAT_INDENT,
)
    head = expand(s[1])
    if head === :function
        format_function(io, s; linelength = linelength, indent = indent)
    elseif head === :for
        format_for(io, s; linelength = linelength, indent = indent)
    elseif head === :begin
        format_begin(io, s; linelength = linelength, indent = indent)
    else
        format_call(io, s; linelength = linelength, indent = indent)
    end
end

function format(
    s::SExpr;
    linelength = DEFAULT_FORMAT_LINELENGTH,
    indent = DEFAULT_FORMAT_INDENT,
)
    format(stdout, s; linelength = linelength, indent = indent)
end

function format(
    io::IO,
    s::Atom;
    linelength = DEFAULT_FORMAT_LINELENGTH,
    indent = DEFAULT_FORMAT_INDENT,
)
    show(io, s)
end

function format(
    s::Atom;
    linelength = DEFAULT_FORMAT_LINELENGTH,
    indent = DEFAULT_FORMAT_INDENT,
)
    format(stdout, s; linelength = linelength, indent = indent)
end

function format_call(
    io::IO,
    s::SExpr;
    linelength = DEFAULT_FORMAT_LINELENGTH,
    indent = DEFAULT_FORMAT_INDENT,
)
    _format(io, s; indent = indent) =
        format(io, s; linelength = linelength, indent = indent)

    operator_str = let io = IOBuffer()
        print(io, "(")
        _format(io, s[1])
        length(s) > 1 && print(io, " ")
        String(take!(io))
    end
    print(io, operator_str)
    length(s) == 1 && (print(io, ")"); return)

    operands_str = let io = IOBuffer()
        for (i, expr) in enumerate(s[2:end])
            i > 1 && print(io, " ")
            _format(io, expr)
        end
        print(io, ")")
        String(take!(io))
    end
    max_operand_linelength = maximum(length.(split(operands_str, "\n")))
    if max_operand_linelength + indent <= linelength
        print(io, operands_str)
    else
        let indent = indent + length(operator_str)
            for (i, expr) in enumerate(s[2:end])
                i > 1 && print(io, "\n", repeat(" ", indent))
                _format(io, expr)
            end
        end
        print(io, ")")
    end
end

function format_function(
    io::IO,
    s::SExpr;
    linelength = DEFAULT_FORMAT_LINELENGTH,
    indent = DEFAULT_FORMAT_INDENT,
)
    _format(io, s; indent = indent) =
        format(io, s; linelength = linelength, indent = indent)

    func_start_str = let io = IOBuffer()
        print(io, "(")
        _format(io, s[1])
        print(io, " ")
        _format(io, s[2])
        String(take!(io))
    end
    if length(func_start_str) <= linelength - indent
        print(io, func_start_str)
    else
        func_start_str = let io = IOBuffer()
            print(io, "(")
            _format(io, s[1])
            print(io, " (")
            func_formals = s[2]
            _format(io, func_formals[1])
            String(take!(io))
        end
        if length(func_formals) > 1
            let indent = indent + length(func_start_str)
                for v in func_formals[2:end]
                    print(io, "\n")
                    print(io, repeat(" ", indent))
                    _format(io, v; indent = indent)
                end
            end
        end
        print(io, ")")
    end
    let indent = indent + 2
        for expr in s[3:end]
            print(io, "\n")
            print(io, repeat(" ", indent))
            _format(io, expr; indent = indent)
        end
    end
    print(io, ")")
end

function format_for(
    io::IO,
    s::SExpr;
    linelength = DEFAULT_FORMAT_LINELENGTH,
    indent = DEFAULT_FORMAT_INDENT,
)
    _format(io, s; indent = indent) =
        format(io, s; linelength = linelength, indent = indent)

    for_head = let io = IOBuffer()
        print(io, "(")
        _format(io, s[1])
        print(io, " (")
        String(take!(io))
    end
    print(io, for_head)
    let indent = indent + length(for_head)
        for (i, iter_expr) in enumerate(s[2])
            i > 1 && print(io, "\n", repeat(" ", indent))
            _format(io, iter_expr)
        end
    end
    print(io, ")")
    let indent = indent + 2
        for body_expr in s[3:end]
            print(io, "\n", repeat(" ", indent))
            _format(io, body_expr)
        end
    end
    print(io, ")")
end

function format_begin(
    io::IO,
    s::SExpr;
    linelength = DEFAULT_FORMAT_LINELENGTH,
    indent = DEFAULT_FORMAT_INDENT,
)
    _format(io, s; indent = indent) =
        format(io, s; linelength = linelength, indent = indent)

    print(io, "(")
    _format(io, s[1])
    let indent = indent + 2
        for body_expr in s[2:end]
            print(io, "\n", repeat(" ", indent))
            _format(io, body_expr)
        end
    end
    print(io, ")")
end

function format_file(
    io::IO,
    filename;
    linelength = DEFAULT_FORMAT_LINELENGTH,
    indent = DEFAULT_FORMAT_INDENT,
)
    for expr in parse(SExpr, read(filename, String))
        format(expr; linelength = linelength, indent = indent)
        println()
        println()
    end
end

function format_file(
    filename;
    linelength = DEFAULT_FORMAT_LINELENGTH,
    indent = DEFAULT_FORMAT_INDENT,
)
    format_file(stdout, filename; linelength = linelength, indent = indent)
end

const DELIM_REGEX = r"\s|\(|\)|\[|\]|\{|\}|;"

function Base.parse(::Type{SExpr}, x::AbstractString)
    s = SExpr([])
    s_stack = [s]
    i = firstindex(x)
    j = i
    while true
        # @show i x[i:j]
        delim = findfirst(DELIM_REGEX, x[i:end])
        delim === nothing && break
        j = last(delim) - 1 + i
        c = x[i:j]
        if c == "("
            push!(s, SExpr([]))
            push!(s_stack, s)
            s = s[end]
            i = nextind(x, j)
            j = i
        elseif c == "["
            push!(s, SExpr(Any[:ref]))
            push!(s_stack, s)
            s = s[end]
            i = nextind(x, j)
            j = i
        elseif c == "{"
            push!(s, SExpr(Any[:curly]))
            push!(s_stack, s)
            s = s[end]
            i = nextind(x, j)
            j = i
        elseif c == ")" || c == "]" || c == "}"
            s = pop!(s_stack)
            i = nextind(x, j)
            j = i
        elseif all(isspace, c)
            i = nextind(x, j)
            j = i
        elseif c == ";"
            next_newline = findfirst("\n", x[i:end])
            i = last(next_newline) - 1 + i
            j = i
        else
            token_start = i
            if x[token_start] == '"'
                quote_end = findfirst(r"[^\\]\"", x[nextind(x, i):end])
                if quote_end === nothing
                    throw(ArgumentError("Unterminated string found at index $(i)"))
                end
                token_end = last(quote_end)
            end
            token_end = prevind(x, j)
            token = x[token_start:token_end]
            push!(s, parse(Atom, token))
            i = nextind(x, token_end)
            j = i
        end
        i > lastindex(x) && break
    end
    return first(s_stack)
end

expand(x) = x
expand(x::Atom) = x.val
expand(x::SExpr) = expand(x.terms)
expand(x::LineNumberNode) = x

struct ModuleExpr
    module_spec::Any
end

function expand(x::ModuleExpr)
    x_str = expand(x.module_spec) |> string
    leading_dots_end = findfirst(r"[^.]", x_str)
    x_name = Symbol(x_str[last(leading_dots_end):end])
    n_dots = 1 + length(1:(prevind(x_str, first(leading_dots_end))))
    return Expr(vcat(repeat([:(.)], n_dots), [x_name])...)
end

struct FunctionExpr
    function_spec::Any
end

function expand(x::FunctionExpr)
    # (function (f <args>...) <body>...)
    # <args> <- <pos_arg> | <kw_arg>
    # <pos_arg> <- <symbol>
    # <kw_arg> <- (<symbol> [<default>])
    # ->
    # (function (call f (parameters [<kw_args>...]) <pos_args>...) (block <body>...))
    function_name = expand(x.function_spec[2][1])
    pos_args = []
    kw_args = []
    for a in x.function_spec[2][2:end]
        if a isa Atom
            push!(pos_args, expand(a))
        elseif a isa SExpr
            ae = a.terms
            if length(ae) == 1
                push!(kw_args, (ae[1], nothing))
            elseif length(ae) == 2
                push!(kw_args, (ae[1], ae[2]))
            else
                throw(ArgumentError("Keyword argument must have form `(var)` or `(var default)`, got $(x.function_spec)"))
            end
        end
    end

    if length(kw_args) > 0
        @show function_name pos_args kw_args
        function_def = Expr(
            :call,
            function_name,
            Expr(
                :parameters,
                [
                    kw[2] !== nothing ? Expr(:kw, expand.(kw)...) : expand(kw[1])
                    for kw in kw_args
                ]...,
            ),
            pos_args...,
        )
    else
        function_def = Expr(:call, function_name, pos_args...)
    end

    body = expand.(x.function_spec[3:end])

    return Expr(:function, function_def, Expr(:block, body...))
end

function expand(x::Array)
    @show x
    head = expand(x[1])
    @show head
    if head ∉ KEYWORDS || (head isa Array && x[1] == :curly)
        # check for macros
        if first(string(head)) == '@'
            return expand(vcat([:macrocall, x[1], LineNumberNode(1, :none)], x[2:end]))
        else
            return expand(vcat([:call], x))
        end
    elseif head === :function
        return expand(FunctionExpr(x))
    elseif head === :begin
        return Expr(:block, expand.(x[2:end])...)
    elseif head === :(=)
        return Expr(head, expand.(x[2:end])...)
    elseif head === :using || head === :import
        module_exprs = expand.(ModuleExpr.(x[2:end]))
        return Expr(head, module_exprs...)
    elseif head === :(.)
        return Expr(head, expand(x[2]), QuoteNode(expand(x[3])))
    elseif head === :for
        # (for (<iter>...) <body>)
        # <iter> <- (<symbol> <expr>)
        # ->
        # (for (block (= <iter>...)) (block <body>)))
        iteration_expr = [Expr(:(=), expand.(it)...) for it in x[2]]
        @show iteration_expr
        body = Expr(:block, expand.(x[3:end])...)
        if length(iteration_expr) == 1
            return Expr(head, iteration_expr[1], body)
        else
            return Expr(head, Expr(:block, iteration_expr...), body)
        end
    elseif head === :if
        if length(x) == 3
            # (if <pred> <value>)
            # ->
            # (if <pred> (block <value>))
            return Expr(head, expand(x[2]), Expr(:block, expand(x[3])))
        elseif length(x) == 4
            # (if <pred> <true-value> <false-value>)
            # ->
            # (if <pred> (block <true-value>) (block <false-value>))
            return Expr(
                head,
                expand(x[2]),
                Expr(:block, expand(x[3])),
                Expr(:block, expand(x[4])),
            )
        elseif length(x) > 4
            # (if <pred> <body>... [elseif <pred> <body>...]... [else <body>...])
            # ->
            # (if <pred> (block <body>...)
            #     ([elseif <pred>]/else (block <body>) [elseif/else...]))
            i_else = findall(a -> a === :else, x)
            i_elseif = findall(a -> a === :elseif, x)
            if length(i_else) == 1
                tail_i = i_else[1]
                expr = Expr(:block, expand.(x[tail_i+1:end])...)
            elseif length(i_elseif) >= 1
                tail_i = pop!(i_elseif)
                expr = Expr(
                    :elseif,
                    Expr(:block, expand(x[tail_i+1])),
                    Expr(:block, expand.(x[tail_i+2:end])...),
                )
            else
                # No elseif or else blocks.
                return Expr(
                    :if,
                    Expr(:block, expand(x[2])),
                    Expr(:block, expand.(x[3:end])...),
                )
            end
            # Now walk backward for any elseif statements.
            while length(i_elseif) >= 1
                i = pop!(i_elseif)
                expr = Expr(
                    :elseif,
                    Expr(:block, expand(x[i+1])),
                    Expr(:block, expand.(x[i+2:tail_i-1])...),
                    expr,
                )
                tail_i = i
            end
            #Finally add the initial if block
            return Expr(
                :if,
                Expr(:block, expand(x[2])),
                Expr(:block, expand.(x[3:tail_i-1])...),
                expr,
            )
        end
    else
        return Expr(head, expand.(x[2:end])...)
    end
end

### REPL functionality
function enable_repl()
    initrepl(
        (x) -> map(eval ∘ expand, parse(SExpr, x));
        prompt_text = "arraylisp> ",
        prompt_color = :blue,
        start_key = ')',
        mode_name = "ArrayLisp_mode",
    )
end

function _sexpr_str(sexpr_str)
    s = parse(SExpr, sexpr_str)
    terms = []
    for si in s
        @show si
        e = si |> expand
        @show e
        push!(terms, :(@eval($(e))))
    end
    return Expr(:block, terms...)
end

macro sexpr_str(sexpr_str)
    return _sexpr_str(sexpr_str)
end

function include_lisp(mod, names...)
    for name in names
        s = read(name, String)
        println(s)
        mod.eval(_sexpr_str(s))
    end
end

### Testing utilities
macro dump_eval(expr)
    Meta.dump(expr)
    return :(@eval($(expr)))
end

end # module


text:
include("../LoadGraph/LoadGraph.jl")
using .LoadGraph: load_adjacency_list, AdjacencyList

function main(dir_path="../../cleaned_data", dtype::Type=Int32)
    for file in filter(x->endswith(x, ".txt"), readdir(dir_path))
        if file != "friendster.txt"
        println(file)
        graph = load_adjacency_list("$dir_path/$file", dtype)
        # @time triangles = list_triangles(graph, dtype)
        # println(length(triangles))
        @time n_triangles = count_triangles(graph, dtype)
        println(n_triangles)
        println()
        end
    end
end

function main2()
    G = AdjacencyList{Int}(()->Vector{Int}())
    G[1] = [2, 3, 4]
    G[2] = [1, 3, 4]
    G[3] = [1, 2]
    G[4] = [1, 2]
    triangles = list_triangles(G,Int)
    println(triangles)
    println()
end

function list_triangles(graph::AdjacencyList, dtype::Type=UInt32)
    truncated = reindex_truncate(graph, dtype)
    triangles = find_triangles(truncated)
    return triangles
end

function count_triangles(graph::AdjacencyList, dtype::Type=UInt32)
    truncated = reindex_truncate(graph, dtype)
    n_triangles = count_triangles(truncated)
    return n_triangles
end

function reindex_truncate(graph::AdjacencyList, dtype::Type=UInt32)
    sorted_keys = sort(collect(keys(graph)),
                       by=k->length(graph[k]),
                       rev=true)
    new_keys = collect(dtype, 1:length(graph))
    key_map = Dict(zip(sorted_keys, new_keys))

    reindexed = [sort([key_map[key]
                       for key in graph[old_key]
                       if key_map[key] > new_key])
                 for (old_key, new_key) in zip(sorted_keys, new_keys)]
    return reindexed
end

function find_triangles(graph::Vector{Vector{T}}) where T<:Real
    triangles = Vector{Vector{T}}()
    for (node, neighbours) in enumerate(graph)
        for (i, neighbour) in enumerate(neighbours[1:end-1])
            for second_neighbour in intersect(neighbours[i+1:end], graph[neighbour])
                push!(triangles, [node, neighbour, second_neighbour])
            end
        end
    end
    return triangles
end

function count_triangles(graph::Vector{Vector{T}}) where T<:Real
    n_triangles = 0
    for (node, neighbours) in enumerate(graph)
        for (i, neighbour) in enumerate(neighbours[1:end-1])
            n_triangles += length(intersect(neighbours[i+1:end], graph[neighbour]))
        end
    end
    return n_triangles
end

function degree(graph::AdjacencyList, node)
    return length(graph[node])
end

main2()
main()

max_stars_repo_name: quentin-burthier/GrandsGraphes

max_issues_repo_name: quentin-burthier/GrandsGraphes

max_forks_repo_name: quentin-burthier/GrandsGraphes


text:
# This file constains testset for Printer

@testset "PrinterTestSet" begin
    @info "Running PrinterTestSet ..."

    # Printer construction
    printer = Printer(input=Inport(2), buflen=100)
    @test typeof(printer.trigger) == Inpin{Float64}
    @test typeof(printer.handshake) == Outpin{Bool}
    @test size(printer.timebuf) == (100,)
    @test size(printer.databuf) == (2, 100)
    @test isa(printer.input, Inport)
    @test printer.plugin === nothing
    @test typeof(printer.callbacks) <: Nothing
    @test typeof(printer.sinkcallback) <: Callback

    # Driving Printer
    oport, iport, trg, hnd, tsk, tsk2 = prepare(printer)
    for t in 1 : 200
        put!(trg, t)
        put!(oport, ones(2) * t)
        take!(hnd)
        @test read(printer.timebuf) == t
        @test [read(pin.links[1].buffer) for pin in oport] == ones(2) * t
    end
    put!(trg, NaN)
    sleep(0.1)
    @test istaskdone(tsk)

    @info "Done PrinterTestSet ..."
end  # testset


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module QuadGK
export gauss, kronrod, quadgk
using Base.Collections
import Base: isless, Order.Reverse, AnyDict

# Adaptive Gauss-Kronrod quadrature routines (arbitrary precision),
# written and contributed to Julia by Steven G. Johnson, 2013.

###########################################################################

# cache of (T,n) -> (x,w,gw) Kronrod rules, to avoid recomputing them
# unnecessarily for repeated integration.   We initialize it with the
# default n=7 rule for double-precision calculations.
const rulecache = AnyDict( (Float64,7) => # precomputed in 100-bit arith.
  ([-9.9145537112081263920685469752598e-01,
    -9.4910791234275852452618968404809e-01,
    -8.6486442335976907278971278864098e-01,
    -7.415311855993944398638647732811e-01,
    -5.8608723546769113029414483825842e-01,
    -4.0584515137739716690660641207707e-01,
    -2.0778495500789846760068940377309e-01,
    0.0],
   [2.2935322010529224963732008059913e-02,
    6.3092092629978553290700663189093e-02,
    1.0479001032225018383987632254189e-01,
    1.4065325971552591874518959051021e-01,
    1.6900472663926790282658342659795e-01,
    1.9035057806478540991325640242055e-01,
    2.0443294007529889241416199923466e-01,
    2.0948214108472782801299917489173e-01],
    [1.2948496616886969327061143267787e-01,
     2.797053914892766679014677714229e-01,
     3.8183005050511894495036977548818e-01,
     4.1795918367346938775510204081658e-01]) )

# integration segment (a,b), estimated integral I, and estimated error E
immutable Segment
    a::Number
    b::Number
    I
    E::Real
end
isless(i::Segment, j::Segment) = isless(i.E, j.E)

# Internal routine: approximately integrate f(x) over the interval (a,b)
# by evaluating the integration rule (x,w,gw). Return a Segment.
function evalrule(f, a,b, x,w,gw, nrm)
    # Ik and Ig are integrals via Kronrod and Gauss rules, respectively
    s = convert(eltype(x), 0.5) * (b-a)
    n1 = 1 - (length(x) & 1) # 0 if even order, 1 if odd order
    # unroll first iterationof loop to get correct type of Ik and Ig
    fg = f(a + (1+x[2])*s) + f(a + (1-x[2])*s)
    fk = f(a + (1+x[1])*s) + f(a + (1-x[1])*s)
    Ig = fg * gw[1]
    Ik = fg * w[2] + fk * w[1]
    for i = 2:length(gw)-n1
        fg = f(a + (1+x[2i])*s) + f(a + (1-x[2i])*s)
        fk = f(a + (1+x[2i-1])*s) + f(a + (1-x[2i-1])*s)
        Ig += fg * gw[i]
        Ik += fg * w[2i] + fk * w[2i-1]
    end
    if n1 == 0 # even: Gauss rule does not include x == 0
        Ik += f(a + s) * w[end]
    else # odd: don't count x==0 twice in Gauss rule
        f0 = f(a + s)
        Ig += f0 * gw[end]
        Ik += f0 * w[end] +
              (f(a + (1+x[end-1])*s) + f(a + (1-x[end-1])*s)) * w[end-1]
    end
    Ik *= s
    Ig *= s
    E = nrm(Ik - Ig)
    if isnan(E) || isinf(E)
        throw(DomainError())
    end
    return Segment(a, b, Ik, E)
end

rulekey(::Type{BigFloat}, n) = (BigFloat, precision(BigFloat), n)
rulekey(T,n) = (T,n)

# Internal routine: integrate f over the union of the open intervals
# (s[1],s[2]), (s[2],s[3]), ..., (s[end-1],s[end]), using h-adaptive
# integration with the order-n Kronrod rule and weights of type Tw,
# with absolute tolerance abstol and relative tolerance reltol,
# with maxevals an approximate maximum number of f evaluations.
function do_quadgk{Tw}(f, s, n, ::Type{Tw}, abstol, reltol, maxevals, nrm)
    if eltype(s) <: Real # check for infinite or semi-infinite intervals
        s1 = s[1]; s2 = s[end]; inf1 = isinf(s1); inf2 = isinf(s2)
        if inf1 || inf2
            if inf1 && inf2 # x = t/(1-t^2) coordinate transformation
                return do_quadgk(t -> begin t2 = t*t; den = 1 / (1 - t2);
                                            f(t*den) * (1+t2)*den*den; end,
                                 map(x -> isinf(x) ? copysign(one(x), x) : 2x / (1+hypot(1,2x)), s),
                                 n, Tw, abstol, reltol, maxevals, nrm)
            end
            s0,si = inf1 ? (s2,s1) : (s1,s2)
            if si < 0 # x = s0 - t/(1-t)
                return do_quadgk(t -> begin den = 1 / (1 - t);
                                            f(s0 - t*den) * den*den; end,
                                 reverse!(map(x -> 1 / (1 + 1 / (s0 - x)), s)),
                                 n, Tw, abstol, reltol, maxevals, nrm)
            else # x = s0 + t/(1-t)
                return do_quadgk(t -> begin den = 1 / (1 - t);
                                            f(s0 + t*den) * den*den; end,
                                 map(x -> 1 / (1 + 1 / (x - s0)), s),
                                 n, Tw, abstol, reltol, maxevals, nrm)
            end
        end
    end

    key = rulekey(Tw,n)
    x,w,gw = haskey(rulecache, key) ? rulecache[key] :
       (rulecache[key] = kronrod(Tw, n))
    segs = Segment[]
    for i in 1:length(s) - 1
        heappush!(segs, evalrule(f, s[i],s[i+1], x,w,gw, nrm), Reverse)
    end
    numevals = (2n+1) * length(segs)
    I = segs[1].I
    E = segs[1].E
    for i in 2:length(segs)
        I += segs[i].I
        E += segs[i].E
    end
    # Pop the biggest-error segment and subdivide (h-adaptation)
    # until convergence is achieved or maxevals is exceeded.
    while E > abstol && E > reltol * nrm(I) && numevals < maxevals
        s = heappop!(segs, Reverse)
        mid = (s.a + s.b) * 0.5
        s1 = evalrule(f, s.a, mid, x,w,gw, nrm)
        s2 = evalrule(f, mid, s.b, x,w,gw, nrm)
        heappush!(segs, s1, Reverse)
        heappush!(segs, s2, Reverse)
        I = (I - s.I) + s1.I + s2.I
        E = (E - s.E) + s1.E + s2.E
        numevals += 4n+2
    end
    # re-sum (paranoia about accumulated roundoff)
    I = segs[1].I
    E = segs[1].E
    for i in 2:length(segs)
        I += segs[i].I
        E += segs[i].E
    end
    return (I, E)
end

# Gauss-Kronrod quadrature of f from a to b to c...

function quadgk{T<:AbstractFloat}(f, a::T,b::T,c::T...;
                                  abstol=zero(T), reltol=sqrt(eps(T)),
                                  maxevals=10^7, order=7, norm=vecnorm)
    do_quadgk(f, [a, b, c...], order, T, abstol, reltol, maxevals, norm)
end

function quadgk{T<:AbstractFloat}(f, a::Complex{T},
                                  b::Complex{T},c::Complex{T}...;
                                  abstol=zero(T), reltol=sqrt(eps(T)),
                                  maxevals=10^7, order=7, norm=vecnorm)
    do_quadgk(f, [a, b, c...], order, T, abstol, reltol, maxevals, norm)
end

"""
    quadgk(f, a,b,c...; reltol=sqrt(eps), abstol=0, maxevals=10^7, order=7, norm=vecnorm)

Numerically integrate the function `f(x)` from `a` to `b`, and optionally over additional
intervals `b` to `c` and so on. Keyword options include a relative error tolerance `reltol`
(defaults to `sqrt(eps)` in the precision of the endpoints), an absolute error tolerance
`abstol` (defaults to 0), a maximum number of function evaluations `maxevals` (defaults to
`10^7`), and the `order` of the integration rule (defaults to 7).

Returns a pair `(I,E)` of the estimated integral `I` and an estimated upper bound on the
absolute error `E`. If `maxevals` is not exceeded then `E <= max(abstol, reltol*norm(I))`
will hold. (Note that it is useful to specify a positive `abstol` in cases where `norm(I)`
may be zero.)

The endpoints `a` et cetera can also be complex (in which case the integral is performed over
straight-line segments in the complex plane). If the endpoints are `BigFloat`, then the
integration will be performed in `BigFloat` precision as well.

!!! note
    It is advisable to increase the integration `order` in rough proportion to the
    precision, for smooth integrands.

More generally, the precision is set by the precision of the integration
endpoints (promoted to floating-point types).

The integrand `f(x)` can return any numeric scalar, vector, or matrix type, or in fact any
type supporting `+`, `-`, multiplication by real values, and a `norm` (i.e., any normed
vector space). Alternatively, a different norm can be specified by passing a `norm`-like
function as the `norm` keyword argument (which defaults to `vecnorm`).

!!! note
    Only one-dimensional integrals are provided by this function. For multi-dimensional
    integration (cubature), there are many different algorithms (often much better than simple
    nested 1d integrals) and the optimal choice tends to be very problem-dependent. See the
    Julia external-package listing for available algorithms for multidimensional integration or
    other specialized tasks (such as integrals of highly oscillatory or singular functions).

The algorithm is an adaptive Gauss-Kronrod integration technique: the integral in each
interval is estimated using a Kronrod rule (`2*order+1` points) and the error is estimated
using an embedded Gauss rule (`order` points). The interval with the largest error is then
subdivided into two intervals and the process is repeated until the desired error tolerance
is achieved.

These quadrature rules work best for smooth functions within each interval, so if your
function has a known discontinuity or other singularity, it is best to subdivide your
interval to put the singularity at an endpoint. For example, if `f` has a discontinuity at
`x=0.7` and you want to integrate from 0 to 1, you should use `quadgk(f, 0,0.7,1)` to
subdivide the interval at the point of discontinuity. The integrand is never evaluated
exactly at the endpoints of the intervals, so it is possible to integrate functions that
diverge at the endpoints as long as the singularity is integrable (for example, a `log(x)`
or `1/sqrt(x)` singularity).

For real-valued endpoints, the starting and/or ending points may be infinite. (A coordinate
transformation is performed internally to map the infinite interval to a finite one.)
"""
# generic version: determine precision from a combination of
# all the integration-segment endpoints
function quadgk(f, a, b, c...; kws...)
    T = promote_type(typeof(float(a)), typeof(b))
    for x in c
        T = promote_type(T, typeof(x))
    end
    cT = map(T, c)
    quadgk(f, convert(T, a), convert(T, b), cT...; kws...)
end

###########################################################################
# Gauss-Kronrod integration-weight computation for arbitrary floating-point
# types and precision, implemented based on the description in:
#
#    Dirk P. Laurie, "Calculation of Gauss-Kronrod quadrature rules,"
#    Mathematics of Computation, vol. 66, no. 219, pp. 1133-1145 (1997).
#
# for the Kronrod rule, and for the Gauss rule from the description in
#
#    Lloyd N. Trefethen and David Bau, Numerical Linear Algebra (SIAM, 1997).
#
# Arbitrary-precision eigenvalue (eignewt & eigpoly) and eigenvector
# (eigvec1) routines are written by SGJ, independent of the above sources.
#
# Since we only implement Gauss-Kronrod rules for the unit weight function,
# the diagonals of the Jacobi matrices are zero and certain things simplify
# compared to the general case of an arbitrary weight function.

# Given a symmetric tridiagonal matrix H with H[i,i] = 0 and
# H[i-1,i] = H[i,i-1] = b[i-1], compute p(z) = det(z I - H) and its
# derivative p'(z), returning (p,p').
function eigpoly(b,z,m=length(b)+1)
    d1 = z
    d1deriv = d2 = one(z)
    d2deriv = zero(z)
    for i = 2:m
        b2 = b[i-1]^2
        d = z * d1 - b2 * d2
        dderiv = d1 + z * d1deriv - b2 * d2deriv
        d2 = d1
        d1 = d
        d2deriv = d1deriv
        d1deriv = dderiv
    end
    return (d1, d1deriv)
end

# compute the n smallest eigenvalues of the symmetric tridiagonal matrix H
# (defined from b as in eigpoly) using a Newton iteration
# on det(H - lambda I).  Unlike eig, handles BigFloat.
function eignewt(b,m,n)
    # get initial guess from eig on Float64 matrix
    H = SymTridiagonal(zeros(m), Float64[ b[i] for i in 1:m-1 ])
    lambda0 = sort(eigvals(H))

    lambda = Array{eltype(b)}(n)
    for i = 1:n
        lambda[i] = lambda0[i]
        for k = 1:1000
            (p,pderiv) = eigpoly(b,lambda[i],m)
            lambda[i] = (lamold = lambda[i]) - p / pderiv
            if abs(lambda[i] - lamold) < 10 * eps(lambda[i]) * abs(lambda[i])
                break
            end
        end
        # do one final Newton iteration for luck and profit:
        (p,pderiv) = eigpoly(b,lambda[i],m)
        lambda[i] = lambda[i] - p / pderiv
    end
    return lambda
end

# given an eigenvalue z and the matrix H(b) from above, return
# the corresponding eigenvector, normalized to 1.
function eigvec1(b,z::Number,m=length(b)+1)
    # "cheat" and use the fact that our eigenvector v must have a
    # nonzero first entries (since it is a quadrature weight), so we
    # can set v[1] = 1 to solve for the rest of the components:.
    v = Array{eltype(b)}(m)
    v[1] = 1
    if m > 1
        s = v[1]
        v[2] = z * v[1] / b[1]
        s += v[2]^2
        for i = 3:m
            v[i] = - (b[i-2]*v[i-2] - z*v[i-1]) / b[i-1]
            s += v[i]^2
        end
        scale!(v, 1 / sqrt(s))
    end
    return v
end

# return (x, w) pair of N quadrature points x[i] and weights w[i] to
# integrate functions on the interval (-1, 1).  i.e. dot(f(x), w)
# approximates the integral.  Uses the method described in Trefethen &
# Bau, Numerical Linear Algebra, to find the N-point Gaussian quadrature
function gauss{T<:AbstractFloat}(::Type{T}, N::Integer)
    if N < 1
        throw(ArgumentError("Gauss rules require positive order"))
    end
    o = one(T)
    b = T[ n / sqrt(4n^2 - o) for n = 1:N-1 ]
    x = eignewt(b,N,N)
    w = T[ 2*eigvec1(b,x[i])[1]^2 for i = 1:N ]
    return (x, w)
end

# Compute 2n+1 Kronrod points x and weights w based on the description in
# Laurie (1997), appendix A, simplified for a=0, for integrating on [-1,1].
# Since the rule is symmetric only returns the n+1 points with x <= 0.
# Also computes the embedded n-point Gauss quadrature weights gw (again
# for x <= 0), corresponding to the points x[2:2:end].  Returns (x,w,wg).
function kronrod{T<:AbstractFloat}(::Type{T}, n::Integer)
    if n < 1
        throw(ArgumentError("Kronrod rules require positive order"))
    end
    o = one(T)
    b = zeros(T, 2n+1)
    b[1] = 2*o
    for j = 1:div(3n+1,2)
        b[j+1] = j^2 / (4j^2 - o)
    end
    s = zeros(T, div(n,2) + 2)
    t = zeros(T, div(n,2) + 2)
    t[2] = b[n+2]
    for m = 0:n-2
        u = zero(T)
        for k = div(m+1,2):-1:0
            l = m - k + 1
            k1 = k + n + 2
            u += b[k1]*s[k+1] - b[l]*s[k+2]
            s[k+2] = u
        end
        s,t = t,s
    end
    for j = div(n,2):-1:0
        s[j+2] = s[j+1]
    end
    for m = n-1:2n-3
        u = zero(T)
        for k = m+1-n:div(m-1,2)
            l = m - k + 1
            j = n - l
            k1 = k + n + 2
            u -= b[k1]*s[j+2] - b[l]*s[j+3]
            s[j+2] = u
        end
        k = div(m+1,2)
        if 2k != m
            j = n - (m - k + 2)
            b[k+n+2] = s[j+2] / s[j+3]
        end
        s,t = t,s
    end
    for j = 1:2n
        b[j] = sqrt(b[j+1])
    end

    # get negative quadrature points x
    x = eignewt(b,2n+1,n+1) # x <= 0

    # get quadrature weights
    w = T[ 2*eigvec1(b,x[i],2n+1)[1]^2 for i in 1:n+1 ]

    # Get embedded Gauss rule from even-indexed points, using
    # the method described in Trefethen and Bau.
    for j = 1:n-1
        b[j] = j / sqrt(4j^2 - o)
    end
    gw = T[ 2*eigvec1(b,x[i],n)[1]^2 for i = 2:2:n+1 ]

    return (x, w, gw)
end

###########################################################################

end # module QuadGK

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using Test, Random, FillArrays
import LuxurySparse: IMatrix, PermMatrix

Random.seed!(2)

p1 = IMatrix{4}()
sp = sprand(ComplexF64, 4,4, 0.5)
ds = rand(ComplexF64, 4,4)
pm = PermMatrix([2,3,4,1], randn(4))
v = [0.5, 0.3im, 0.2, 1.0]
dv = Diagonal(v)

@testset "basic" begin
    @test p1==copy(p1)
    @test eltype(p1) == Bool
    @test size(p1) == (4, 4)
    @test size(p1, 1) == size(p1, 2) == 4
    @test Matrix(p1) == [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
end

@testset "conversion" begin
    for mat in [p1, pm, dv]
        @test mat == SparseMatrixCSC(mat)
        @test mat == Matrix(mat)
    end
    for mat in [p1, pm, dv]
        @test mat == PermMatrix(mat)
    end
    @test Diagonal(p1) == p1

    @test SparseMatrixCSC(Diagonal(Fill(2, 4))) ≈ Diagonal(Fill(2, 4))
end

@testset "sparse" begin
    @test nnz(p1) == 4
    @test nonzeros(p1) == ones(4)
end

@testset "linalg" begin
    for op in [conj, real, transpose, copy, inv]
        @test op(p1) == Matrix(I, 4, 4)
        @test typeof(op(p1)) == typeof(p1)
    end
    @test imag(p1) == zeros(4, 4)
    @test p1' == Matrix(I, 4, 4)

    # This will be lazy evaluated in 0.7+
    @static if VERSION < v"0.7-"
        @test typeof(p1') == typeof(p1)
    end

    @test ishermitian(p1)
end

@testset "elementary" begin
    @test all(isapprox.(conj(p1), conj(Matrix(p1))))
    @test all(isapprox.(real(p1), real(Matrix(p1))))
    @test all(isapprox.(imag(p1), imag(Matrix(p1))))
end

@testset "basicmath" begin
    @test p1*2im == Matrix(p1)*2im
    @test p1/2.0 == Matrix(p1)/2.0
end


text:
# replstr(x) = sprint((io,x) -> show(IOContext(io, :limit => true, :displaysize => (24, 80)), MIME("text/plain"), x), x)
# showstr(x) = sprint((io,x) -> show(IOContext(io, :limit => true, :displaysize => (24, 80)), x), x)

@testset "Bias correction" begin

d = DateTime(1961,1,1):Day(1):DateTime(1990,12,31)
Random.seed!(42)
data = randn(2, 2, 10957)
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
dimension_dict = Dict(["lon" => "lon", "lat" => "lat"])
obs = ClimateTools.ClimGrid(axisdata, variable = "tasmax", dimension_dict=dimension_dict)
ref = obs * 1.05
fut = obs * 1.05

D = qqmap(obs, ref, fut, method = "Additive", detrend=false)
@test round(D[1][1, 1, 1], digits=5) == -0.55603#68761463861
D = qqmap(obs, ref, fut, method = "Additive", detrend=true)
@test round(D[1][1, 1, 1], digits=2) == -0.56#5603#74620422795

# ===================================
# Ensure we have similar statistics
# ===================================
d = DateTime(1961,1,1):Day(1):DateTime(1990,12,31)
Random.seed!(42)
data = randn(2, 2, 10957)
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
dimension_dict = Dict(["lon" => "lon", "lat" => "lat"])
obs = ClimateTools.ClimGrid(axisdata, variable = "tasmax", dimension_dict=dimension_dict)
ref = obs + 2.0
fut = obs + 2.0

D = qqmap(obs, ref, fut, method="Additive", detrend=false)
@test mean(obs[1][1,1,:]) - mean(fut[1][1,1,:]) ≈ -1.99999999999
@test mean(obs[1][1,1,:]) - mean(D[1][1,1,:]) ≈ -0.00029395108

ref = obs * 2.0
fut = obs * 2.0
D = qqmap(obs, ref, fut, method="Multiplicative", detrend=false)
@test std(obs[1][1,1,:]) - std(fut[1][1,1,:]) ≈ -1.012096874
@test round(std(obs[1][1,1,:]) - std(D[1][1,1,:]), digits=3) ≈ -0.0#5.233448404e-5

# ============================
Random.seed!(42)
data = randn(2, 2, 10957)
axisdata = AxisArray(data .- minimum(data), Axis{:lon}(1:2), Axis{:lat}(1:2),Axis{:time}(d))
obs = ClimateTools.ClimGrid(axisdata, variable = "tasmax", dimension_dict=dimension_dict)
ref = obs * 1.05
fut = obs * 1.05

D = qqmap(obs, ref, fut, method = "Multiplicative", detrend=false)
@test round(D[1][1, 1, 1], digits=2) == 3.78

# Bias extremes
d = DateTime(1961,1,1):Day(1):DateTime(1990,12,31)
lat = 1.0:2.0
lon = 1.0:2.0
Random.seed!(42)
data = randn(2, 2, 10957)
longrid, latgrid = ndgrid(lon, lat)
axisdata = AxisArray(data, Axis{:lon}(lon), Axis{:lat}(lat), Axis{:time}(d))
dimension_dict = Dict(["lon" => "lon", "lat" => "lat"])
varattribs = Dict(["standard_name" => "random noise"])
obs = ClimateTools.ClimGrid(axisdata, variable = "psl", typeofvar="psl", longrid=longrid, latgrid=latgrid, dimension_dict=dimension_dict, varattribs=varattribs)
ref = obs * 1.1
fut = obs * 1.05

mu = [36.6208, 37.423]
sigma = [11.555, 10.381]
xi = [0.08393, 0.08393]
gevparams = DataFrame([lat, lon, mu, sigma, xi], [:lat, :lon, :mu, :sigma, :xi])
Dext = biascorrect_extremes(obs, ref, fut, detrend=true, gevparams=gevparams, method="multiplicative")
@test round(maximum(Dext), digits=5) == 103.03357

Dext = biascorrect_extremes(obs, ref, fut, detrend=false, gevparams=gevparams, method="multiplicative")
@test round(maximum(Dext), digits=5) == 103.03378

# Create a ClimGrid with a clear trend
x = 1:10957
x = x
d = DateTime(1961,1,1):Day(1):DateTime(1990,12,31)
data = zeros(2, 2, 10957)
data[1,1,:] .= 0.0 .+ x
data[1,2,:] .= 0.0 .+ 0.2.*x .+ x.^2
data[2,1,:] .= 0.0 .+ 0.04.*x - 0.01 .* x.^3
data[2,2,:] .= 0.0 .+ 0.14.*x.^2 .- 0.032 .* x .^4
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimGrid(axisdata)
poly = polyfit(C)
val = polyval(C, poly)
D = C - val
@test D[1] == (C - val)[1]

end


text:
###########################
### Basic Pull Requests ###
###########################

listRepoPullRequests(c::Client,owner::String,repo::String,page::Int,state::String) =  getParsedResponse(Vector{PullRequest},c,Requests.get,"/repos/$(owner)/$(repo)/pulls"; json = Dict("page" => page, "state" => state))

getPullRequest(c::Client,owner::String,repo::String,index::Int) =  getParsedResponse(PullRequest,c,Requests.get,"/repos/$(owner)/$(repo)/pulls/$(index)")

FieldTags.@tag immutable CreatePullRequestOption
	title::String
	head::String
	base::String
	body::Nullable{String}
	assignee::Nullable{String}
	milestone::Nullable{Int64}
	labels::Nullable{Vector{Int64}}
end

function CreatePullRequestOption(title, head, base; kwargs...)
	kwargsDict = Dict{Symbol,Any}(k=>v for (k,v) in kwargs)

	args = Any[title,head,base]

	for field in fieldnames(CreatePullRequestOption)[4:end]
		push!(args,get(kwargsDict,field, fieldtype(CreatePullRequestOption,field)()))
	end

	CreatePullRequestOption(args...)
end

createPullRequest(c::Client,owner::String,repo::String,opt::CreatePullRequestOption) =  getParsedResponse(PullRequest,c,Requests.post,"/repos/$(owner)/$(repo)/pulls"; json = marshalJSON(opt))

FieldTags.@tag immutable EditPullRequestOption
	title::Nullable{String}
	body::Nullable{String}
	assignee::Nullable{String}
	milestone::Nullable{Int64}
	labels::Nullable{Vector{Int64}}
	state::Nullable{String}
end

function EditPullRequestOption(;kwargs...)
	kwargsDict = Dict{Symbol,Any}(k=>v for (k,v) in kwargs)

	args = Any[]

	for field in fieldnames(EditPullRequestOption)[1:end]
		push!(args,get(kwargsDict,field, fieldtype(EditPullRequestOption,field)()))
	end

	EditPullRequestOption(args...)
end

editPullRequest(c::Client,owner::String,repo::String,index::Int, opt::EditPullRequestOption) =  getParsedResponse(PullRequest,c,Requests.patch,"/repos/$(owner)/$(repo)/pulls/$(index)"; json = marshalJSON(opt))

mergePullRequest(c::Client,owner::String,repo::String,index::Int) =  getParsedResponse(PullRequest,c,Requests.post,"/repos/$(owner)/$(repo)/pulls/$(index)/merge")

function isPullRequestMerged(c::Client, user::String, repo::String, index::String)
	code = getStatusCode(c,Requests.get,"/repos/$(user)/$(repo)/pulls/$(index)/merge")
	return code == 204
end


text:
# ---
# title: 1226. The Dining Philosophers
# id: problem1226
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/the-dining-philosophers/description/>
# hidden: true
# ---
#
# Five silent philosophers sit at a round table with bowls of spaghetti. Forks
# are placed between each pair of adjacent philosophers.
#
# Each philosopher must alternately think and eat. However, a philosopher can
# only eat spaghetti when they have both left and right forks. Each fork can be
# held by only one philosopher and so a philosopher can use the fork only if it
# is not being used by another philosopher. After an individual philosopher
# finishes eating, they need to put down both forks so that the forks become
# available to others. A philosopher can take the fork on their right or the one
# on their left as they become available, but cannot start eating before getting
# both forks.
#
# Eating is not limited by the remaining amounts of spaghetti or stomach space;
# an infinite supply and an infinite demand are assumed.
#
# Design a discipline of behaviour (a concurrent algorithm) such that no
# philosopher will starve;  _i.e._ , each can forever continue to alternate
# between eating and thinking, assuming that no philosopher can know when others
# may want to eat or think.
#
# ![](https://assets.leetcode.com/uploads/2019/09/24/an_illustration_of_the_dining_philosophers_problem.png)
#
# _The problem statement and the image above are taken
# from[wikipedia.org](https://en.wikipedia.org/wiki/Dining_philosophers_problem)_
#
#
#
# The philosophers' ids are numbered from **0** to **4** in a **clockwise**
# order. Implement the function `void wantsToEat(philosopher, pickLeftFork,
# pickRightFork, eat, putLeftFork, putRightFork)` where:
#
#   * `philosopher` is the id of the philosopher who wants to eat.
#   * `pickLeftFork` and `pickRightFork` are functions you can call to pick the corresponding forks of that philosopher.
#   * `eat` is a function you can call to let the philosopher eat once he has picked both forks.
#   * `putLeftFork` and `putRightFork` are functions you can call to put down the corresponding forks of that philosopher.
#   * The philosophers are assumed to be thinking as long as they are not asking to eat (the function is not being called with their number).
#
# Five threads, each representing a philosopher, will simultaneously use one
# object of your class to simulate the process. The function may be called for
# the same philosopher more than once, even before the last call ends.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 1
#     Output: [[4,2,1],[4,1,1],[0,1,1],[2,2,1],[2,1,1],[2,0,3],[2,1,2],[2,2,2],[4,0,3],[4,1,2],[0,2,1],[4,2,2],[3,2,1],[3,1,1],[0,0,3],[0,1,2],[0,2,2],[1,2,1],[1,1,1],[3,0,3],[3,1,2],[3,2,2],[1,0,3],[1,1,2],[1,2,2]]
#     Explanation:
#     n is the number of times each philosopher will call the function.
#     The output array describes the calls you made to the functions controlling the forks and the eat function, its format is:
#     output[i] = [a, b, c] (three integers)
#     - a is the id of a philosopher.
#     - b specifies the fork: {1 : left, 2 : right}.
#     - c specifies the operation: {1 : pick, 2 : put, 3 : eat}.
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 60`
#
#
## @lc code=start
using LeetCode

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


text:
doc=
"""
usage:
stencil.jl <outputdir> <builddir> <refine> <scheme> <vtk>

     builddir path      Directory to write discretization.
     outputdir path     Directory to write output data to.
     refine int         Level of grid refinement to use. refine = 0 (no grid
                                                                  refinement)
     scheme string      Choose between 'staggered' or 'collocated'
     vtk     int        How often to write vtk files

Generate input data for stencil-based computation.

"""
if length(ARGS) != 5
        println(doc)
        exit(-1)
end

builddir=string(ARGS[1], "/")
outputdir=string(ARGS[2], "/")
refine=parse(Int32, ARGS[3])
scheme=ARGS[4]
vtk_stride=parse(Int32, ARGS[5])

if scheme == "collocated"
        scheme_num = 0
elseif scheme == "staggered"
        scheme_num = 1
else
        error("Unknown scheme selected")
end

import sbp

using Printf
using LinearAlgebra
using SparseArrays

num_moment=4
num_smoothness=4
ref=2^refine
mx=64*ref
my=32*ref
nx=mx + 1
ny=my + 1
r1_s=0.5
r2_s=0.5
r1_r=0.4
r2_r=0.4
dt=2.0/mx
num_steps=250*ref
scheme=0
println(dt)

if scheme_num == 1
        println("Using staggered discretization")
        ac = sbp.StaggeredAcoustic
else
        println("Using collocated discretization")
        ac = sbp.CollocatedAcoustic
end

# Number of grid points and grids
if scheme_num == 1
        nxp, nyp = ac.grid_points("p", nx, ny)
        nx1, ny1 = ac.grid_points("v1", nx, ny)
        nx2, ny2 = ac.grid_points("v2", nx, ny)

        r1_p, r2_p = ac.grids("p", nx, ny, pad=false)
        r1_1, r2_1 = ac.grids("v1", nx, ny, pad=false)
        r1_2, r2_2 = ac.grids("v2", nx, ny, pad=false)
        r1, r2 = ac.grids("node", nx, ny, pad=false)

        pr1_p, pr2_p = ac.grids("p", nx, ny, pad=true)
        pr1_1, pr2_1 = ac.grids("v1", nx, ny, pad=true)
        pr1_2, pr2_2 = ac.grids("v2", nx, ny, pad=true)
        pr1, pr2 = ac.grids("node", nx, ny, pad=true)
else
        nxp = nx
        nyp = ny
        nx1 = nx
        ny1 = ny
        nx2 = nx
        ny2 = ny
        r1_p, r2_p = ac.grid(nx, ny)
        r1_1, r2_1 = ac.grid(nx, ny)
        r1_2, r2_2 = ac.grid(nx, ny)
        r1, r2 = ac.grid(nx, ny)
end

function mapping(r1, r2)
        x = 10.0 .* r1
        sigma = 10/1.05
        y = (5 .+ 1 .* exp.(-sigma^2*(r1 .- 0.5).^2)) .* r2
        return x, y
end

function initial_conditions(xp, yp, x1, y1, x2, y2)
        p = zeros(size(xp))
        v1 = zeros(size(x1))
        v2 = zeros(size(x2))
        return p, v1, v2
end

function zero_pad_x(a, qx, qy)
        b = zeros(qx* qy)

        for i=1:(qx-1)
        for j=1:qy
                b[j + (i - 1) * qy] = a[j + (i - 1) * qy]
        end
        end
        return b
end

function zero_pad_y(a, qx, qy)
        b = zeros(qx* qy)

        for i=1:qx
        for j=1:(qy-1)
                b[j + (i - 1) * qy] = a[j + (i - 1) * (qy - 1)]
        end
        end
        return b
end

x, y = mapping(r1, r2)
xp, yp = mapping(r1_p, r2_p)
x1, y1 = mapping(r1_1, r2_1)
x2, y2 = mapping(r1_2, r2_2)

if scheme_num == 1
        pxp, pyp = mapping(pr1_p, pr2_p)
        px1, py1 = mapping(pr1_1, pr2_1)
        px2, py2 = mapping(pr1_2, pr2_2)
        p, v1, v2 = initial_conditions(pxp, pyp, px1, py1, px2, py2)
        x1 = zero_pad_x(x1, nx + 1, ny + 1)
        y1 = zero_pad_x(y1, nx + 1, ny + 1)
        x2 = zero_pad_y(x2, nx + 1, ny + 1)
        y2 = zero_pad_y(y2, nx + 1, ny + 1)

else
        p, v1, v2 = initial_conditions(xp, yp, x1, y1, x2, y2)
end

xs, ys = mapping(r1_s, r2_s)
xr, yr = mapping(r1_r, r2_r)
@printf("Source location: (%f, %f) \n", xs, ys)
@printf("Receiver location: (%f, %f) \n", xr, yr)

sizes = [0, size(p,1), size(v1,1), size(v2,1)]
ranges = accumulate(+, sizes)

@time begin
println("Discretizing in space...")
@printf("refine = %d mx = %d my = %d dt = %g \n", refine, mx, my, dt)
if scheme_num == 1
        ops = ac.init_operators(nx, ny, sbp.OP2019.build_operators, x, y)
else
        st = sbp.Strand1994
        desc_D, desc_H = sbp.Strand1994.load_description(order=4)
        ops = ac.init_operators(nx, ny, st, desc_D, desc_H, x, y)
end
end

"""
Export to C++ code

"""
h1 = 1 / mx
h2 = 1 / my
cfg = Dict("nx" => nxp,
           "ny" => nyp,
           "dt" => dt,
           "hi1" => 1/h1,
           "hi2" => 1/h2,
           "nvtk" => vtk_stride,
           "nt" => num_steps,
           "outputdir" => outputdir,
           "scheme" => scheme_num
                   )

if scheme_num == 0
        r1, h1 = sbp.Grid.grid_xp(nx)
        r2, h2 = sbp.Grid.grid_xp(ny)
        else
        r1, h1 = sbp.Grid.grid_xm(nx + 1)
        r2, h2 = sbp.Grid.grid_xm(ny + 1)
end
d = Array(sbp.Source.source_discretize_2d(r1_s, r2_s, num_moment,
                                    num_smoothness, r1, r2, h1,
                                    h2))
file = x -> string(builddir, x)

@time begin
println("Writing discretization to ", builddir)
sbp.lmv.write_config(file("config.txt"), cfg, skip_types=true)

arr = x -> convert(Array{Float32}, Array(diag(x)))
conv = x -> convert(Array{Float32}, Array(x))

if scheme_num == 0
        sbp.lmv.write_vector(file("J.bin"), arr(ops.J))
        sbp.lmv.write_vector(file("g11.bin"), arr(ops.G.b11))
        sbp.lmv.write_vector(file("g12.bin"), arr(ops.G.b12))
        sbp.lmv.write_vector(file("g22.bin"), arr(ops.G.b22))
else
        Jp = arr(ops.Jp)
        J1 = conv(zero_pad_x(arr(ops.J1), nx+1, ny+1))
        J2 = conv(zero_pad_y(arr(ops.J2), nx+1, ny+1))
        g11 = conv(zero_pad_x(arr(ops.G1.b11), nx+1, ny+1))
        g12 = arr(ops.Gp.b12)
        g22 = conv(zero_pad_y(arr(ops.G2.b22), nx+1, ny+1))

        sbp.lmv.write_vector(file("Jp.bin"), Jp)
        sbp.lmv.write_vector(file("J1.bin"), J1)
        sbp.lmv.write_vector(file("J2.bin"), J2)
        sbp.lmv.write_vector(file("g1_11.bin"), g11)
        sbp.lmv.write_vector(file("gp_12.bin"), g12)
        sbp.lmv.write_vector(file("g2_22.bin"), g22)
end
        sbp.lmv.write_vector(file("xp.bin"), conv(xp))
        sbp.lmv.write_vector(file("yp.bin"), conv(yp))
        sbp.lmv.write_vector(file("d.bin"), conv(d))
end

sbp.lmv.write_vector(file("v1.bin"), conv(v1))
sbp.lmv.write_vector(file("v2.bin"), conv(v2))
sbp.lmv.write_vector(file("p.bin"), conv(p))

max_stars_repo_name: ooreilly/sbp

max_issues_repo_name: ooreilly/sbp

max_forks_repo_name: ooreilly/sbp


text:
using SmithNormalForm
using BenchmarkTools

srand(2903872398473)
rows, cols = 64, 93
TMP = rand(1:10, rows*cols)
D = reshape(TMP, rows, cols)
i, j = 2, 3
a, b, c, d = 1, 2, 3, 4

suite = BenchmarkGroup()

suite["internal"] = BenchmarkGroup(["swap", "elimination"])
# suite["snf"] = BenchmarkGroup(["decomposition"])

suite["internal"]["col_swap"] = @benchmarkable SmithNormalForm.cswap!(D, i, j) gctrial=true evals=1000 samples=100000
suite["internal"]["row_swap"] = @benchmarkable SmithNormalForm.rswap!(D, i, j) gctrial=true evals=1000 samples=100000

suite["internal"]["col_elim"] = @benchmarkable SmithNormalForm.colelimination(D, a, b, c, d, i, j) evals=500 gctrial=true
suite["internal"]["row_elim"] = @benchmarkable SmithNormalForm.rowelimination(D, a, b, c, d, i, j) evals=500 gctrial=true

Z = rand(1:100, 500, 520)
Z[Z .> 66] = 1
Z[Z .> 33] = -1
Z[Z .> 0] = 0
suite["snf"]["dense"] = @benchmarkable snf(Z)

Z = sprand(Int, 500, 520, 0.3)
Z[Z .< 0] = -1
Z[Z .> 0] = 1
suite["snf"]["sparse"] = @benchmarkable snf(Z)