text:
module PortableMinskyMachineNotation
export @pmmn_str, compile, load
using SparseArrays

struct PMMNError <: Exception
    msg::String
end

function inc(memory::SparseVector{I}, n, i::Integer = oneunit(I)) where {I<:Integer}
    memory[n+oneunit(n)] += i
    nothing
end

function dec(memory::SparseVector{I}, n, i::Integer = oneunit(I)) where {I<:Integer}
    n += oneunit(n)
    before = memory[n]
    memory[n] = max(zero(I), memory[n] - i)
    memory[n] != before
end

function input(io::IO, memory, n::Integer)
    if !eof(io)
        memory[n+oneunit(n)] += read(io, UInt8) + 1
    end
    nothing
end

function output(io::IO, memory, n::Integer)
    n += oneunit(n)
    if iswritable(io) && !iszero(memory[n])
        write(io, Char((memory[n] - 1) & 255))
    end
    memory[n] = zero(memory[n])
    nothing
end

function scan(source)
    (m.match for m in eachmatch(
        r"\d+|}|end|if|else|while|inc|dec|input|output",
        replace(source, r"(
            /\*(?:[^/]*|(?<!\*)/)*\*/|      # C style multiline comments
            \#=(?:[^#]*|(?<!\=)\#)=\#|      # Julia style multiline comments
            (\#|//).*                       # single line comments
        )"x => "")
    ))
end

parse(tokens) = parse(Iterators.Stateful(tokens))
parse(tokens::Iterators.Stateful) = commandlist(tokens)
parse(tokens::AbstractString) = parse(scan(tokens))

function expect(value::AbstractString, token::AbstractString)
    if token != value
        throw(PMMNError("Expected '$(value)', but got '$(token)'."))
    end
    token
end

function expect(value::AbstractVector, token::AbstractString)
    if !(token in value)
        throw(PMMNError("Expected '$(first(value))', but got '$(token)'."))
    end
    token
end

expect(value, tokens) = expect(value, popfirst!(tokens))

function commandlist(tokens)
    result = Expr[]
    while !isempty(tokens) && peek(tokens) != "}" && peek(tokens) != "end"
        push!(result, command(tokens))
    end
    if isempty(result)
        throw(PMMNError("Expected a command list, but found nothing."))
    end
    Expr(:block, result...)
end

function command(tokens, token = popfirst!(tokens))
    if token == "inc" || token == "dec"
        expr = Expr(:call, Symbol(token), :memory, Base.parse(Int, popfirst!(tokens)))
        if all(isdigit, peek(tokens))
            push!(expr.args, Base.parse(Int, popfirst!(tokens)))
        end
        expr
    elseif token == "if"
        condition = command(tokens, expect("dec", token))
        ifblock = block(tokens)
        if peek(tokens) == "else"
            Expr(:if, condition, ifblock, block(tokens))
        else
            Expr(:if, condition, ifblock)
        end
    elseif token == "while"
        condition = command(tokens, expect("dec", tokens))
        whileblock = block(tokens)
        Expr(:while, condition, whileblock)
    elseif token == "input" || token == "output"
        Expr(
            :call,
            Symbol(token),
            Symbol(replace(token, "put" => "")),
            :memory,
            Base.parse(Int, popfirst!(tokens))
        )
    else
        PMMNError("Unexpected token: $(token)")
    end
end

function block(tokens)
    expr = commandlist(tokens)
    expect(["}", "end"], tokens)
    expr
end

function compile(commands::Expr)
    fn = :(
        function (in = stdin; out = IOBuffer())
            if in isa AbstractString
                in = IOBuffer(in)
            end
            memory = spzeros(BigInt, 2000000000)
        end
    )
    push!(
        last(fn.args).args,
        commands,
        :(
            if out != stdout
                readchomp(seekstart(out))
            end
        )
    )
    eval(fn)
end

"""

    compile(str::AbstractString)

Compiles a string containing a Portable Minsky Machine Notation program.
"""
compile(str::AbstractString) = compile(parse(str))

"""

    pmmn""

A convenience macro for writing Portable Minsky Machine Notation.
"""
macro pmmn_str(str)
    :(compile($(str)))
end

"""

    load(filepath::AbstractString)

Load a Minsky machine from the specified file.
"""
load(filepath::AbstractString) = compile(readchomp(filepath))

end # module

max_stars_repo_name: stellartux/Esolangs

max_issues_repo_name: stellartux/Esolangs

max_forks_repo_name: stellartux/Esolangs


text:
import SolverTraces: scalar_format

@testset "Scalar column" begin
    @test format(scalar_format(Float64, false), 1.0) == " 1.000e+00"
    @test format(scalar_format(Float64, false), 100.0) == " 1.000e+02"
    @test format(scalar_format(Float64, true), 1.0) == "+1.000e+00"
    @test format(scalar_format(Float64, true), -1.0) == "-1.000e+00"

    @test format(scalar_format(Int, false), 1) == "   1 "
    @test format(scalar_format(Int, false), 100) == " 100 "
    @test format(scalar_format(Int, true), 1) == "  +1 "
    @test format(scalar_format(Int, true), -1) == "  -1 "
    @test format(scalar_format(Int, true), 100) == "+100 "
    @test format(scalar_format(Int, true), -100) == "-100 "

    scf = ScalarColumn(identity, "", Float64)
    @test format(scf, 1) == " 1.000e+00"
    @test format(scf, 100) == " 1.000e+02"

    sci = ScalarColumn(i -> i^2, "", Int)
    @test format(sci, 1) == "   1 "
    @test format(sci, 10) == " 100 "

    sci = ScalarColumn(5, "")
    @test format(sci, 1) == "   5 "
    @test format(sci, 10) == "   5 "
    sci.n = 6
    @test format(sci, 1) == "   6 "
    @test format(sci, 10) == "   6 "

    sci = ScalarColumn(Inf, "")
    @test format(sci, 1) == "       Inf"
    @test format(sci, 10) == "       Inf"
    sci.n = 6
    @test format(sci, 1) == " 6.000e+00"
    @test format(sci, 10) == " 6.000e+00"
end

ext: jl


text:
#g=g_u(uvec)
f_target(x,Amp,h1,Lt,Ht) = Amp^2*((abs(x[1])<Lt/2)&&(abs(x[2]-h1/2)<Ht/2))

function g_u(u_vec;Amp,h1,Lt,Ht,U,V,dΩ_t)
  It(x) = f_target(x,Amp,h1,Lt,Ht)
  uh_t = FEFunction(U,u_vec)
  sum(∫(abs2(abs2(uh_t)-It))dΩ_t)
end

#uvec = u_pf(pf)
function u_pf(pf;x0,δ,Amp,P,Pf,β,η,flag_t,flag_f,ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,hd,U,V,dΩ,dΓ)
  if (flag_f)
      pfh = FEFunction(Pf,pf)
  else
      pfh = FEFunction(P,pf)
  end
  ph = (pf->Threshold(β,η,flag_t,pf))∘pfh
  A_mat = MatrixA(ph,ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,hd,U,V,dΩ)
  B_vec = MatrixB(x0,δ,2*π*Amp,V,dΩ,dΓ)
  u_vec = A_mat\B_vec
  u_vec
end

#pf = pf_p(p)
function pf_p(p;r,flag_f,P,Pf,Qf,dΩ,dΓ_d,tags,design_tag)
  pvec = p_vec(p,P,tags,design_tag)
  pf = Filter(pvec,r,flag_f,P,Pf,Qf,dΩ,dΓ_d)
  pf
end
# Chain Rule : dg/dp = dg/dg*dg/du*du/dpf*dpf/dp
# dg/du=dg/dg*dg/du
function rrule(::typeof(g_u),u_vec;Amp,h1,Lt,Ht,U,V,dΩ_t)
function g_pullback(dgdg)
  NO_FIELDS, dgdg*Dgdu(u_vec,Amp,h1,Lt,Ht,U,V,dΩ_t)
end
g_u(u_vec;Amp,h1,Lt,Ht,U,V,dΩ_t), g_pullback
end

function Dgdu(u_vec,Amp,h1,Lt,Ht,U,V,dΩ_t)
It(x) = f_target(x,Amp,h1,Lt,Ht)
uh_t = FEFunction(U,u_vec)
l_temp(du)=∫(4*uh_t*(abs2(uh_t)-It)*du)dΩ_t
assemble_vector(l_temp,V)
end

# dg/dpf=dg/du*du/dpf
function rrule(::typeof(u_pf),pf;x0,δ,Amp,P,Pf,β,η,flag_t,flag_f,ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,hd,U,V,dΩ,dΓ)
u_vec = u_pf(pf;x0,δ,Amp,P,Pf,β,η,flag_t,flag_f,ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,hd,U,V,dΩ,dΓ)
function u_pullback(dgdu)
  NO_FIELDS, Dgdpf(dgdu,u_vec,pf,P,Pf,β,η,flag_t,flag_f,ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,hd,U,V,dΩ)
end
u_vec, u_pullback
end

Dξdp(pf,ϵmin,ϵmax,β,η,flag_t)=(ϵmax-ϵmin)*(!flag_t+flag_t*β*(1.0-tanh(β*(pf-η))^2)/(tanh(β*η)+tanh(β*(1.0-η))))
dG(pfh,u,v,dp,ϵmin,ϵmax,k,β,η,flag_t) = real(k^2*((pf->Dξdp(pf,ϵmin,ϵmax,β,η,flag_t))∘pfh)*(v*u)*dp)

function Dgdpf(dgdu,u_vec,pf,P,Pf,β,η,flag_t,flag_f,ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,hd,U,V,dΩ)
  if (flag_f)
      pfh = FEFunction(Pf,pf)
      ph = (pf->Threshold(β,η,flag_t,pf))∘pfh
      A_mat = MatrixA(ph,ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,hd,U,V,dΩ)
      λ_vec = A_mat'\dgdu

      uh = FEFunction(U,u_vec)
      λh = FEFunction(V,conj(λ_vec))
      l_temp(dp) = ∫(dG(pfh,uh,λh,dp,ϵ1,ϵ2,k,β,η,flag_t))*dΩ
      dgdpf = assemble_vector(l_temp,Pf)
      return dgdpf
  else
      pfh = FEFunction(P,pf)
      ph = (pf->Threshold(β,η,flag_t,pf))∘pfh
      A_mat = MatrixA(ph,ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,hd,U,V,dΩ)
      λ_vec = A_mat'\dgdu

      uh = FEFunction(U,u_vec)
      λh = FEFunction(V,conj(λ_vec))
      l_temp(dp) = ∫(dG(pfh,uh,λh,dp,ϵ1,ϵ2,k,β,η,flag_t))*dΩ
      dgdpf = assemble_vector(l_temp,P)
      return dgdpf
  end
end

# dg/dp=dg/dpf*dpf/dp
function rrule(::typeof(pf_p),p;r,flag_f,P,Pf,Qf,dΩ,dΓ_d,tags,design_tag)
function pf_pullback(dgdpf)
  NO_FIELDS, Dgdp(dgdpf,p,r,flag_f,P,Pf,Qf,dΩ,dΓ_d,tags,design_tag)
end
pf_p(p;r,flag_f,P,Pf,Qf,dΩ,dΓ_d,tags,design_tag), pf_pullback
end

function Dgdp(dgdpf,p,r,flag_f,P,Pf,Qf,dΩ,dΓ_d,tags,design_tag)
np = length(p)
if (flag_f)
  A = assemble_matrix(Pf,Qf) do u, v
          ∫( a_f(r,u,v))dΩ
      end
  λvec = A'\dgdpf
  λh = FEFunction(Pf,λvec)
  l_temp(dp) = ∫(λh*dp)*dΩ
  return extract_design(assemble_vector(l_temp,P),np,tags,design_tag)
else
  return extract_design(dgdpf,np,tags,design_tag)
end
end

# Final objective function
function g_p(p::Vector;x0,δ,Amp,r,flag_f,P,Pf,Qf,β,η,flag_t,
      ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,hd,h1,Lt,Ht,U,V,dΩ,dΓ,dΩ_t,dΓ_d,tags,design_tag)
  pf = pf_p(p;r,flag_f,P,Pf,Qf,dΩ,dΓ_d,tags,design_tag)
  u_vec=u_pf(pf;x0,δ,Amp,P,Pf,β,η,flag_t,flag_f,ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,hd,U,V,dΩ,dΓ)
  g_u(u_vec;Amp,h1,Lt,Ht,U,V,dΩ_t)
end

function g_p(p::Vector,grad::Vector;x0,δ,Amp,r,flag_f,P,Pf,Qf,β,η,flag_t,
      ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,hd,h1,Lt,Ht,U,V,dΩ,dΓ,dΩ_t,dΓ_d,tags,design_tag)
  if length(grad) > 0
      dgdp, = Zygote.gradient(p->g_p(p;x0,δ,Amp,r,flag_f,P,Pf,Qf,β,η,flag_t,
      ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,hd,h1,Lt,Ht,U,V,dΩ,dΓ,dΩ_t,dΓ_d,tags,design_tag),p)
      grad[:] = dgdp
  end
  g_value = g_p(p;x0,δ,Amp,r,flag_f,P,Pf,Qf,β,η,flag_t,
      ϵ1,ϵ2,μ,σs,k,LHp,LHn,dpml,hd,h1,Lt,Ht,U,V,dΩ,dΓ,dΩ_t,dΓ_d,tags,design_tag)
  @show g_value
  return g_value
end

max_stars_repo_name: Tiamat-Tech/hpinn

max_issues_repo_name: louhz/hpinn

max_forks_repo_name: louhz/hpinn


text:
include("ED.jl")

@testset "ED checks" begin
    void = state_from_integer(0, 1, 2)
    up = state_from_integer(1, 1, 2)
    down = state_from_integer(2, 1, 2)
    updown = state_from_integer(3, 1, 2)

    # create up at site 1
    _sign, s = create(void, 1, 1)
    @test _sign == 1.0 &&  s == up
    _sign, s = create(up, 1, 1)
    @test _sign == 0 &&  s == 0
    _sign, s = create(down, 1, 1)
    @test _sign == 1.0 &&  s == updown
    _sign, s = create(updown, 1, 1)
    @test _sign == 0 &&  s == 0

    # create down at site 1
    _sign, s = create(void, 1, 2)
    @test _sign == 1.0 &&  s == down
    _sign, s = create(up, 1, 2)
    @test _sign == -1.0 &&  s == updown
    _sign, s = create(down, 1, 2)
    @test _sign == 0 &&  s == 0
    _sign, s = create(updown, 1, 2)
    @test _sign == 0 &&  s == 0

    # annihilate up at site 1
    _sign, s = annihilate(void, 1, 1)
    @test _sign == 0 &&  s == 0
    _sign, s = annihilate(up, 1, 1)
    @test _sign == 1.0 &&  s == void
    _sign, s = annihilate(down, 1, 1)
    @test _sign == 0 &&  s == 0
    _sign, s = annihilate(updown, 1, 1)
    @test _sign == 1.0 &&  s == down

    # annihilate down at site 1
    _sign, s = annihilate(void, 1, 2)
    @test _sign == 0 &&  s == 0
    _sign, s = annihilate(up, 1, 2)
    @test _sign == 0 &&  s == 0
    _sign, s = annihilate(down, 1, 2)
    @test _sign == 1.0 &&  s == void
    _sign, s = annihilate(updown, 1, 2)
    @test _sign == -1.0 &&  s == up

    # check Greens consistency
    model = HubbardModelAttractive(
        L = 2, dims = 2,
        U = rand(), mu = rand(), t = rand()
    )
    H = HamiltonMatrix(model)
    for substate1 in 1:2, substate2 in 1:2
        for site1 in 1:model.l.sites, site2 in 1:model.l.sites
            G = expectation_value(
                Greens(site1, site2, substate1, substate2),
                H,
                N_sites = model.l.sites,
            )
            G_perm = expectation_value(
                Greens_permuted(site1, site2, substate1, substate2),
                H,
                N_sites = model.l.sites,
            )
            @test G ≈ G_perm
        end
    end
end

@testset "Attractive Hubbard Model (ED)" begin
    model = HubbardModelAttractive(
        L = 2,
        dims = 2,
        U = 1.0,
        mu = 1.0,
        t = 1.0
    )

    @info "Running DQMC β=1.0, 10k + 20k sweeps, ≈1min"
    dqmc = DQMC(model, beta=1.0)
    run!(dqmc, thermalization = 10_000, sweeps = 20_000, verbose=false)
    G_DQMC = mean(dqmc.measurements[:Greens].obs)

    @info "Running ED"
    H = HamiltonMatrix(model)
    G_ED = calculate_Greens_matrix(H, model.l, beta=1.0)

    # G_DQMC is smaller because it doesn't differentiate between spin up/down
    for i in 1:size(G_DQMC, 1), j in 1:size(G_DQMC, 2)
        @test isapprox(G_DQMC[i, j], G_ED[i, j], atol=0.025, rtol=0.1)
        # @test isapprox(G_DQMC[i, j], G_ED[i, j])
    end
end


text:
function log_network(state::Tuple{AbstractGraph, AbstractArray}, tick_nr::Int64)
    graph, agent_list = state
    agent_opinion = [a.opinion for a in agent_list]
    agent_perceiv_publ_opinion = [a.perceiv_publ_opinion for a in agent_list]
    agent_inclin_interact = [a.inclin_interact for a in agent_list]
    agent_inactive_ticks = [a.inactive_ticks for a in agent_list]
    agent_active_state = [a.active for a in agent_list]
    agent_indegree = indegree(graph)
    return DataFrame(
        TickNr = tick_nr,
        AgentID = 1:length(agent_list),
        Opinion = agent_opinion,
        PerceivPublOpinion = agent_perceiv_publ_opinion,
        InclinInteract = agent_inclin_interact,
        InactiveTicks = agent_inactive_ticks,
        Indegree = agent_indegree,
        ActiveState = agent_active_state
    )
end

max_stars_repo_name: digitalemuendigkeit/twitter-abm

max_issues_repo_name: digitalemuendigkeit/twitter-abm

max_forks_repo_name: digitalemuendigkeit/twitter-abm


text:
using StatsBase: sample, Weights

const PACKET_SIZE = 8
const WATER_WASTING_FACTOR = 2.0

struct WaterDrop
    "The index of the first block, doubled as RNG seed"
    index::UInt32

    "The number of blocks"
    degree::UInt32

    "The block data XORed together"
    data::Vector{UInt8}
end

function get_block(data::AbstractVector{UInt8}, block_size::Integer, i::Integer)::SubArray{UInt8, 1}
    view(data, (block_size * (i - 1) + 1 : block_size * i))
end

function generate_degrees(soliton::Soliton, n_samples::Int)::Vector{Int}
    prob = get_probabilities(soliton)
    sample(1:length(prob), Weights(prob), n_samples)
end

"`n`: maximum index number."
function generate_indices(seed::UInt32, n::Int, n_samples::Int)::Vector{Int}
    lehmer = Lehmer(seed)
    Int.(random_sample(lehmer, UInt32(n), UInt32(n_samples)))
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
using GeneralizedSampling
using Base.Test

println("Testing NFFT...")

#=
Test if the (N)FFT part of a Freq2Wave type behaves like a DFT matrix by
removing the parts related to wavelets.

For integer samples [0:N-1] the DFT matrix is Hermitian.

The special version of NFFT with uniform sampling locations should agree
with the general version.
=#

# ------------------------------------------------------------
# 1D

begin
	N = 16.0
	samples = collect( 0:N-1 )

	J = Int(log2(N))
	T = Freq2Wave(samples, "haar", J)
	for idx in eachindex(T.internal)
		T.internal[idx] = 1.0 + 0.0*im
	end

	F = collect(T)
	B = F'*F / N

	@test_approx_eq B eye(B)
end

# ------------------------------------------------------------
# 2D

begin
	M = 4
	samples = grid(M)

	J = Int(log2(M))
	T = Freq2Wave(samples, "haar", J)
	for idx in eachindex(T.internal)
		T.internal[idx] = 1.0 + 0.0*im
	end

	F = collect(T)
	N = size(T,2)
	B = F'*F / N

	@test_approx_eq B eye(B)
end

# ------------------------------------------------------------
# Uniform NFFT

J = rand(2:5)
M = 2*rand( 2^(J-1):2^J )
epsilon = 1/rand(2:5)
xi = GeneralizedSampling.grid(M, epsilon)

using NFFT

N1 = 2*rand( 2:2^(J-1) )
chi = xi * 2.0^-J
P1 = NFFTPlan( frac(chi), N1 )

f = rand(N1) + rand(N1)*im
y1 = nfft(P1, f)

P2 = GeneralizedSampling.FFTPlan( xi, J, N1 )
y2 = similar(y1)
GeneralizedSampling.myfft!(y2, P2, f)

@test_approx_eq_eps y1 y2 1e-5


text:
pages = [
    "MethodOfLines.jl: Automated Finite Difference for Phyiscs-Informed Learning" => "index.md",
    "Tutorials" => ["tutorials/brusselator.md", "tutorials/heat.md", "tutorials/params.md", "tutorials/heatss.md", "tutorials/sispde.md", "tutorials/icbc_sampled.md"],
    "MOLFiniteDifference" => "MOLFiniteDifference.md",
    "Grid and Solution Retrieval" => "get_grid.md",
    "Boundary Conditions" => "boundary_conditions.md",
    "Non-Uniform Rectilinear Grids" => "nonuniform.md",
    "Curvilinear Grids" => "curvilinear_grids.md",
    "How it Works" => "howitworks.md",
    "Notes for Developers: Implement a Scheme" => "devnotes.md",
    "Generated Examples" => ["generated/bruss_code.md", "generated/bruss_ode_eqs.md"]
    # "Tutorial: Burgers" => "tutorials/burgers.md",
    # "Tutorial: 1D Linear Diffusion" => "tutorials/1d_linear_diffusion.md",
    # "Tutorial: 1D Non-Linear Diffusion" => "tutorials/1d_nonlinear_diffusion.md",
    # "Tutorial: 2D Diffusion" => "tutorials/2d_diffusion.md",
    # "Tutorial: 1D Linear Convection" => "tutorials/1d_linear_convection.md",
    # "Tutorial: 1D Higher Order" => "tutorials/1d_higher_order.md",
    # "Tutorial: Stationary Nonlinear Problems" => "tutorials/stationary_nonlinear_problems.md",
    # "Tutorial: 1D Partial DAE" => "tutorials/1d_partial_DAE.md",

]


text:
const Point = SVector{2,Float64}

struct HalfPlane
    p::Point
    pq::Point
    angle::Float64

    HalfPlane(p::Point, pq::Point) = begin
        angle = atan(pq[2], pq[1])
        new(p, pq, angle)
    end
end

function out(h::HalfPlane, r::Point)
    cross(h.pq, r - h.p) < -1e-9
end

function Base.isless(x::HalfPlane, y::HalfPlane)
    if abs(x.angle - y.angle) < 1e-6
        return cross(x.pq, y.p - x.p) < 0
    end
    x.angle < y.angle
end

function isect(s::HalfPlane, t::HalfPlane)
    alpha = cross((t.p - s.p), t.pq) / cross(s.pq, t.pq)
    s.p + s.pq * alpha
end

function remove_duplicates(s::Vector{HalfPlane})
    res = HalfPlane[]
    for hp in s
        if length(res) > 0
            if abs(res[end].angle - hp.angle) > 1e-12
                push!(res, hp)
            end
        else
            push!(res, hp)
        end
    end
    res
end

"""
Computes the polygonal intersection of two-dimensional halfplanes. If there is
no intersection, returns `nothing`.

Requires that the halfplanes being intersected have finite area.

Warning: This version of the code does not handle parallel lines correctly.
"""
function intersect_halfplanes(halfplanes::Vector{HalfPlane})
    halfplanes = copy(halfplanes)
    # append!(halfplanes, [
    #     HalfPlane(Point(-1e9, -1e9), Point(1.0, 0.0)),
    #     HalfPlane(Point(1e9, -1e9), Point(0.0, 1.0)),
    #     HalfPlane(Point(1e9, 1e9), Point(-1.0, 0.0)),
    #     HalfPlane(Point(-1e9, 1e9), Point(0.0, -1.0)),
    # ])
    sort!(halfplanes)
    halfplanes = remove_duplicates(halfplanes)
    dq = HalfPlane[]
    for halfplane in halfplanes
        while length(dq) >= 2 && out(halfplane, isect(dq[end], dq[end-1]))
            pop!(dq)
        end

        while length(dq) >= 2 && out(halfplane, isect(dq[1], dq[2]))
            popfirst!(dq)
        end
        push!(dq, halfplane)
    end

    while length(dq) >= 3 && out(dq[1], isect(dq[end], dq[end-1]))
        pop!(dq)
    end
    while length(dq) >= 3 && out(dq[end], isect(dq[1], dq[2]))
        popfirst!(dq)
    end

    len = length(dq)
    if len < 3
        return nothing
    end
    hcat([vec(isect(dq[i], dq[mod1(i + 1, len)])) for i = 1:len]...)
end

function intersect_halfplanes_slow(s::Vector{HalfPlane})
    # s is a vector of half planes in two dimensions.
    A = zeros(Float64, (length(s), 2))
    b = Vector{Float64}(undef, length(s))
    ones = BitSet()
    sign = 1.0
    for i = 1:length(s)
        # maybe we need to flip all the signs here?
        A[i, 1] = sign * s[i].pq[2]
        A[i, 2] = sign * -s[i].pq[1]
        b[i] = sign * (s[i].p[1] * s[i].pq[2] - s[i].p[2] * s[i].pq[1])
    end
    poly = polyhedron(hrep(A, b, ones))
    points = Polyhedra.points(vrep(poly))
    if length(points) < 3
        return nothing
    end
    hcat(points...)
end

"""
Returns the equi-pressure intersection between two tetrahedra as a polygon, as
well as the normal vector to the plane of intersection.

Takes in objects `o1`, `o2`, followed by the indices of the tets in each object.
"""
function intersect_tets(o1::Object, o2::Object, a_face_idx::Int64, b_face_idx::Int64)
    m1 = o1.mesh
    m2 = o2.mesh
    coords_A = transform(m1.verts[:, m1.tets[:, a_face_idx]], o1.pose) # 3 x 4 matrix
    coords_B = transform(m2.verts[:, m2.tets[:, b_face_idx]], o2.pose) # 3 x 4 matrix

    # First, do an initial check for bounding boxes of tetrahedra
    bbox_A = bounding_box(coords_A)
    bbox_B = bounding_box(coords_B)
    if any(bbox_A[2] .< bbox_B[1]) || any(bbox_B[2] .< bbox_A[1])
        return nothing
    end

    function get_equation(coords::Matrix{Float64}, potentials::Vector{Float64})
        ones_arr = ones(size(coords, 2))
        mat = [coords' ones_arr]
        mat \ potentials
    end

    a_pot = get_equation(coords_A, m1.potentials[m1.tets[:, a_face_idx]])
    b_pot = get_equation(coords_B, m2.potentials[m2.tets[:, b_face_idx]])
    plane = a_pot - b_pot  # dot(plane, [x; 1.0]) == 0.0
    if norm(plane[1:3]) < 1e-8
        return nothing
    end

    # Check if the tets actually intersect with the plane.
    values_A = plane[1:3]' * coords_A .+ plane[4]
    values_B = plane[1:3]' * coords_B .+ plane[4]
    if !(
        minimum(values_A) < -1e-6 &&
        maximum(values_A) > 1e-6 &&
        minimum(values_B) < -1e-6 &&
        maximum(values_B) > 1e-6
    )
        return nothing
    end

    # Find a 2x3 projection from the plane onto two dimensions, along with an
    # inverse 3x2 projection that has the following properties:
    #   1. proj_mat * inv_proj_mat = I
    #   2. plane[1:3]' * (inv_proj_mat * x + inv_proj_offset) + plane[4] = 0
    if abs(plane[1]) / norm(plane[1:3]) > 1e-3
        proj_mat = [0.0 1.0 0.0; 0.0 0.0 1.0]
        inv_proj_mat = [(-plane[2]/plane[1]) (-plane[3]/plane[1]); 1.0 0.0; 0.0 1.0]
        inv_proj_offset = [-plane[4] / plane[1], 0.0, 0.0]
    elseif abs(plane[2]) / norm(plane[1:3]) > 1e-3
        proj_mat = [1.0 0.0 0.0; 0.0 0.0 1.0]
        inv_proj_mat = [1.0 0.0; (-plane[1]/plane[2]) (-plane[3]/plane[2]); 0.0 1.0]
        inv_proj_offset = [0.0, -plane[4] / plane[2], 0.0]
    elseif abs(plane[3]) / norm(plane[1:3]) > 1e-3
        proj_mat = [1.0 0.0 0.0; 0.0 1.0 0.0]
        inv_proj_mat = [1.0 0.0; 0.0 1.0; (-plane[1]/plane[3]) (-plane[2]/plane[3])]
        inv_proj_offset = [0.0, 0.0, -plane[4] / plane[3]]
    else
        error("no valid projection found, aborting...")
    end

    # Convert each tetrahedron into four 2D halfplanes, on the projection of
    # the plane to two dimensions along some orthogonal axis.
    halfplanes = HalfPlane[]
    for coords in (coords_A, coords_B)
        for i = 1:4
            potentials = zeros(4)
            potentials[i] = 1.0
            halfspace = get_equation(coords, potentials)
            normal_2d = (halfspace[1:3]' * inv_proj_mat)'
            if norm(normal_2d) > 1e-9
                p =
                    normal_2d * (-halfspace[4] - halfspace[1:3]' * inv_proj_offset) /
                    dot(normal_2d, normal_2d)
                pq = normalize(Point(normal_2d[2], -normal_2d[1]))
                push!(halfplanes, HalfPlane(Point(p), pq))
            end
        end
    end

    all_points = intersect_halfplanes_slow(halfplanes)
    if isnothing(all_points)
        return nothing
    end

    # for p in eachcol(all_points)
    #     print(p)
    #     for hp in halfplanes
    #         if out(hp, p)
    #             println("Halfplanes: ", halfplanes)
    #             println(hp)
    #             println(p)
    #             #error("bad!!")
    #         end
    #     end
    # end

    @assert size(all_points, 2) >= 3 "sanity check length"
    hcat([inv_proj_mat * p + inv_proj_offset for p in eachcol(all_points)]...), plane[1:3]
end

function intersect_polygons(polygonA::Matrix{Float64}, polygonB::Matrix{Float64})
    n, m = size(polygonA, 2), size(polygonB, 2)
    halfplanes::Array{HalfPlane} = []
    for i = 1:n
        push!(
            halfplanes,
            HalfPlane(
                Point(polygonA[:, i]),
                Point(polygonA[:, mod1(i + 1, n)]) - Point(polygonA[:, i]),
            ),
        )
    end
    for i = 1:m
        push!(
            halfplanes,
            HalfPlane(
                Point(polygonB[:, i]),
                Point(polygonB[:, mod1(i + 1, m)]) - Point(polygonB[:, i]),
            ),
        )
    end
    intersect_halfplanes_slow(halfplanes)
end

"""
Result of a force calculation, expressed as vectors in the world frame.
"""
struct ForceResult
    F_AB::SVector{3,Float64} # The force applied to object A from B.
    F_BA::SVector{3,Float64} # The force applied to object B from A.
    τ_AB::SVector{3,Float64} # The torque applied to object A from B.
    τ_BA::SVector{3,Float64} # The torque applied to object B from A.
end

"""
Compute overall force between two tets A[i], B[j] of objects A, B.

Returns the force applied to A (in the direction of A.com - B.com), along with
the location of force applied to A and B.
"""
function tet_force(
    A::Object,
    B::Object,
    i::Int64,
    j::Int64,
)::Union{Tuple{Vector{Float64},Vector{Float64}},Nothing}
    result = intersect_tets(A, B, i, j)
    if isnothing(result)
        return nothing
    end

    intersection_polygon, normal = result
    normalize!(normal)
    if dot(normal, transform(A.mesh.com, A.pose) - transform(B.mesh.com, B.pose)) < 0
        normal *= -1
    end

    total_force = 0.0
    intersection_com = zeros(3)
    total_area = 0.0

    vtx_inds = A.mesh.tets[:, i]
    vtx_coords = transform(A.mesh.verts[:, vtx_inds], A.pose)
    vtx_coords = vcat(vtx_coords, ones(1, 4))  # 4x4 matrix of vtxs padded w 1s
    for i = 3:size(intersection_polygon, 2)
        # calculate barycentric coordinates of each point in the triangle.
        vtxs = intersection_polygon[:, [1, i - 1, i]]
        com = push!(mean(eachcol(vtxs)), 1)
        res = vtx_coords \ com
        pressure = sum(res .* A.mesh.potentials[vtx_inds])
        area = 0.5 * norm(cross(vtxs[:, 1] - vtxs[:, 2], vtxs[:, 1] - vtxs[:, 3]))
        total_force += pressure * area
        total_area += area
        intersection_com += area * com[1:3]
    end

    if total_area < 1e-9
        return nothing
    end

    intersection_com /= total_area
    total_force * normal, intersection_com
end

"""
Computes force on object A due to contact with object B, along with the net
torques on both objects.
"""
function compute_force(A::Object, B::Object)::ForceResult
    if A.mesh.m > B.mesh.m
        # Swap order for efficiency
        result = compute_force(B, A)
        return ForceResult(result.F_BA, result.F_AB, result.τ_BA, result.τ_AB)
    end

    X_BA = inv(B.pose) * A.pose
    force = zeros(3)
    τ_AB = zeros(3)
    τ_BA = zeros(3)

    for i = 1:A.mesh.m
        # Coordinates of tet i in object A, from the frame of object B
        coords_BA = A.mesh.verts[:, A.mesh.tets[:, i]]
        coords_BA = transform(coords_BA, X_BA)
        bb_min, bb_max = bounding_box(coords_BA)
        rect = SI.Rect((bb_min[1], bb_min[2], bb_min[3]), (bb_max[1], bb_max[2], bb_max[3]))

        for elem::SI.SpatialElem in intersects_with(B.mesh.rtree.t, rect)
            j = elem.val
            tets_result = tet_force(A, B, i, j)
            if !isnothing(tets_result)
                force += tets_result[1]
                τ_AB += cross(tets_result[2] - A.mesh.com, tets_result[1])
                τ_BA += cross(tets_result[2] - B.mesh.com, -tets_result[1])
            end
        end
    end
    ForceResult(force, -1 * force, τ_AB, τ_BA)
end

export compute_force

max_stars_repo_name: ekzhang/hydroelastics

max_issues_repo_name: ekzhang/hydroelastics

max_forks_repo_name: ekzhang/hydroelastics


text:
module CITIfile

export rdciti

include("rdciti.jl")

end # module


text:
using MatchCore

comm_monoid = @commutative_monoid (*) 1

fold_mul = @theory begin
	a::Number * b::Number |> a * b
end

t = comm_monoid ∪ fold_mul

@testset "Extraction 1 - Commutative Monoid" begin
    G = EGraph(:(3 * 4))
    saturate!(G, t)
    @test (12 == extract!(G, astsize))

    ex = :(a * 3 * b * 4)
    G = EGraph(cleanast(ex))
	params=SaturationParams(timeout=15)
    saturate!(G, t, params)
    extr = extract!(G, astsize)
	println(extr)
	@test extr == :((12 * a) * b) || extr == :(12 * (a * b)) || extr == :(a * (b * 12)) ||
		extr == :((a * b) * 12) || extr == :((12a) * b) || extr == :(a * (12b)) ||
		extr == :((b * (12a))) || extr == :((b * 12) * a) || extr == :((b * a) * 12) ||
		extr == :(b * (a * 12)) || extr == :((12b) * a)
end

fold_add = @theory begin
	a::Number + b::Number |> a + b
end

@testset "Extraction 2" begin
	comm_group = @abelian_group (+) 0 inv

    t = comm_monoid ∪ comm_group ∪ distrib(:(*), :(+)) ∪ fold_mul ∪ fold_add

    # for i ∈ 1:20
    # sleep(0.3)
    ex = cleanast(:((x*(a+b)) + (y*(a+b))))
    G = EGraph(ex)
    saturate!(G, t)
    # end

    extract!(G, astsize) == :((x+y) * (b+a))
end

@testset "Lazy Extraction 2" begin
	comm_group = @abelian_group (+) 0 inv

    t = comm_monoid ∪ comm_group ∪ distrib(:(*), :(+)) ∪ fold_mul ∪ fold_add

    # for i ∈ 1:20
    # sleep(0.3)
    ex = cleanast(:((x*(a+b)) + (y*(a+b))))
    G = EGraph(ex)
    saturate!(G, t)
    # end

    extract!(G, astsize) == :((x+y) * (b+a))
end

@testset "Extraction - Adding analysis after saturation" begin
    G = EGraph(:(3 * 4))
    addexpr!(G, 12)
    saturate!(G, t)
    addexpr!(G, :(a * 2))
    saturate!(G, t)

    saturate!(G, t)

    @test (12 == extract!(G, astsize))

    # for i ∈ 1:100
    ex = :(a * 3 * b * 4)
    G = EGraph(cleanast(ex))
    analyze!(G, NumberFold)
	params=SaturationParams(timeout=15)
    saturate!(G, comm_monoid, params)

    extr = extract!(G, astsize)

	@test extr == :((12 * a) * b) || extr == :(12 * (a * b)) || extr == :(a * (b * 12)) ||
		extr == :((a * b) * 12) || extr == :((12a) * b) || extr == :(a * (12b)) ||
		extr == :((b * (12a))) || extr == :((b * 12) * a) || extr == :((b * a) * 12)
end

comm_monoid = @commutative_monoid (*) 1

comm_group = @abelian_group (+) 0 inv

powers = @theory begin
	a * a => a^2
	a => a^1
	a^n * a^m => a^(n+m)
end
logids = @theory begin
	log(a^n) => n * log(a)
	log(x * y) => log(x) + log(y)
	log(1) => 0
	log(:e) => 1
	:e^(log(x)) => x
end

t = comm_monoid ∪ comm_group ∪ distrib(:(*), :(+)) ∪ powers ∪ logids  ∪ fold_mul ∪ fold_add

@testset "Complex Extraction" begin
	G = EGraph(:(log(e) * log(e)))
	params=SaturationParams(timeout=8)
	saturate!(G, t, params)
	# display(G.classes);println()
	@test extract!(G, astsize) == 1

	G = EGraph(:(log(e) * (log(e) * e^(log(3)))  ))
	params=SaturationParams(timeout=7)
	saturate!(G, t, params)
	@test extract!(G, astsize) == 3

	@time begin
		G = EGraph(:(a^3 * a^2))
		saturate!(G, t)
		ex = extract!(G, astsize)
	end
	@test ex == :(a^5)

	@time begin
		G = EGraph(:(a^3 * a^2))
		saturate!(G, t)
		ex = extract!(G, astsize)
	end
	@test ex == :(a^5)

	# TODO could serve as example for more advanced
	# symbolic mathematics simplification based on the computation cost
	# of the expressions
	function cust_astsize(n::ENode, g::EGraph, an::Type{<:AbstractAnalysis})
		cost = 1 + arity(n)

		if n.head == :^
			cost += 2
		end

		for id ∈ n.args
	        eclass = geteclass(g, id)
	        !hasdata(eclass, an) && (cost += Inf; break)
	        cost += last(getdata(eclass, an))
	    end
	    return cost
	end

	@time begin
		G = EGraph(:((log(e) * log(e)) * (log(a^3 * a^2))))
		saturate!(G, t)
		ex = extract!(G, cust_astsize)
	end
	@test ex == :(5*log(a)) || ex == :(log(a)*5)
end

function costfun(n::ENode, g::EGraph, an)
	arity(n) != 3 && (return 1)
	left = n.args[2]
	left_class = geteclass(g, left)
	ENode(:a) ∈ left_class.nodes ? 1 : 100
end

moveright = @theory begin
    (:b * (:a * c)) => (:a * (:b * c))
end

expr = :(a * (a * (b * (a * b))))
res = rewrite( expr , moveright)

g = EGraph(expr)
saturate!(g, moveright)
resg = extract!(g, costfun)

@testset "Symbols in Right hand" begin
    @test resg == res == :(a * (a * (a * (b * b))))
end

co = @theory begin
	foo(x ⋅ :bazoo ⋅ :woo) => Σ(:n * x)
end
@testset "Consistency with Matchcore backend" begin
	ex = :(foo(wa(rio) ⋅ bazoo ⋅ woo))
	g = EGraph(ex)
	saturate!(g, co)

	res = extract!(g, astsize)

	resclassic = rewrite(ex, co)

	@test res == resclassic
end

@testset "No arguments" begin
	ex = :(f());
	g = EGraph(ex);
	@test :(f()) == extract!(g, astsize)

	ex = :(f() + g())

	t = @theory begin
		f() + g() => h()
	end;

	_, ex2 = @extract ex t astsize;

    @test ex2 == :(h())
end


text:
using FiniteVolumes
using Printf
using Plots

mesh = PeriodicCartesianMesh(40, 40)
model = LinearAdvectionFlux([1.0, 1.0])

is_in_square(x, side=0.5) = all(0.5-side/2 .<= x .<= 0.5+side/2)

w₀ = [is_in_square(x) ? 1.0 : 0.0 for x in cell_centers(mesh)]

anim = Animation()
function plot_callback(i, t, w)
    if i % 2 == 0
        title = @sprintf "time = %.2f" t
        plot(mesh, w, 1, title=title, clims=(0, 1))
        frame(anim)
    end
end

t, w = FiniteVolumes.run(directional_splitting(model), mesh, w₀,
                         time_step=FixedCourant(0.2),
                         nb_time_steps=200,
                         callback=plot_callback)
gif(anim, "advection.gif")


text:
using Distances

# Linear scan using PQ codebooks no rotation
function linscan_pq(
  B::Matrix{UInt8},           # m-by-n. The database, encoded
  X::Matrix{Cfloat},         # d-by-nq. The queries.
  C::Vector{Matrix{Cfloat}}, # The cluster centers
  b:: Int,                    # Number of bits per code -- log2(h) * m
  k:: Int = 10000)            # Number of knn results to return

  m, n  = size( B )
  d, nq = size( X )

  @show k, nq
  dists = zeros( Cfloat, k, nq )
  res   = zeros(  Cuint, k, nq )

  ccall(("linscan_aqd_query", "src/linscan/cpp/linscan_aqd.so"), Void,
    (Ptr{Cfloat}, Ptr{Cuint}, Ptr{Cuchar}, Ptr{Cfloat},
    Ptr{Cfloat}, Cint, Cuint, Cint, Cint, Cint, Cint, Cint),
    dists, res, B, cat(3,C...), X, Cint(n), Cuint(nq),
    Cint(b), Cint(k), Cint(m), Cint(d), Cint(d/m) )

  return dists, (res.+=1)

end

# Linear scan using OPQ
function linscan_opq(
  B::Matrix{UInt8},           # m-by-n. The database, encoded
  X::Matrix{Cfloat},         # d-by-nq. The queries.
  C::Vector{Matrix{Cfloat}}, # The cluster centers
  b::Int,                     # Number of bits per code -- log2(h) * m
  R::Matrix{Cfloat},         # Rotation matrix
  k::Int = 10000)             # Number of knn results to return

  # Simply rotate
  RX = R'*X

  # And call the function as usual
  return linscan_pq( B, RX, C, b, k )
end

# Linear scan using LSQ, with dbnorms not encoded.
function linscan_lsq(
  B::Matrix{UInt8},           # m-by-n. The database, encoded
  X::Matrix{Cfloat},         # d-by-nq. The queries.
  C::Vector{Matrix{Cfloat}}, # (m-1)-long. The cluster centers
  dbnorms::Vector{Cfloat},   # h-long. Table with database norms
  R::Matrix{Cfloat},         # Rotation matrix
  k::Int = 10000)             # Number of knn results to return

  RX = R' * X;

  m, n  = size( B );
  d, nq = size( RX );
  _, h  = size( C[1] );

  dists = zeros( Cfloat, k, nq );
  res   = zeros(  Cint,  k, nq  );

  ccall(("linscan_aqd_query_extra_byte", "src/linscan/cpp/linscan_aqd_pairwise_byte.so"), Void,
    (Ptr{Cfloat}, Ptr{Cint},
    Ptr{Cuchar}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat},
    Cuint, Cint, Cint, Cint, Cint, Cint),
    dists, res,
    B, RX, hcat(C...), dbnorms,
    Cint(nq), Cint(n), Cint(m), Cint(h), Cint(d), Cint(k) );

  return dists, res

end

# Evaluate ANN search vs ground truth. Produces a recall@N curve.
function eval_recall{T <: Integer}(
  ids_gnd::Vector{T},
  ids_predicted::Matrix{T},
  k::Integer)

  # Modified from Hervé Jégou's test_compute_stats.m code

  nquery = size( ids_predicted, 2 );
  assert( nquery == length( ids_gnd) );

  nn_ranks = zeros( nquery );
  hist_pqc = zeros( k+1 );

  for i = 1:nquery
    gnd_ids = ids_gnd[i];
    nn_pos = find( ids_predicted[:,i] .== gnd_ids );

    if length(nn_pos) == 1
      nn_ranks[i] = nn_pos[1];
    else
      nn_ranks[i] = k+1;
    end
  end

  nn_ranks = sort( nn_ranks );

  recall_at_i = zeros( k );

  for i = [1 2 5 10 20 50 100 200 500 1000 2000 5000 10000]
    if i <= k
      recall_at_i[i] = length( find( (nn_ranks .<= i) .& (nn_ranks .<= k) )) ./ nquery * 100;
      println("r@$(i) = $(recall_at_i[i])");
    end
  end

  for i = 1:k
    recall_at_i[i] = length(find( (nn_ranks .<= i) .& (nn_ranks .<= k) )) ./ nquery;
  end

  return recall_at_i

end

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
## distribution.jl : collection of probability density functions
#
# Representation of some commonly used distributions.
#
# Note: ideally, we want to link with Distributions.jl. However, this is troublesome when
# coupling with QMC points. Minimal requirements are the implementation of rng_native_52 and
# rand, see [RandomNumbers.jl](https://github.com/sunoru/RandomNumbers.jl) and the Julia
# documentation [Random Numbers](https://docs.julialang.org/en/v1/stdlib/Random/index.html).
# As a temporary fix, we will define some commonly used probability distributions, and
# define the required transformations ourselves (with less efficiency).
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

## AbstractDistribution ##
abstract type AbstractDistribution end

## Uniform ##
"""
    Uniform([a = 0, [b = 1]])

Return a uniform distribution on [`a`, `b`].

# Examples
```jldoctest; setup = :(using MultilevelEstimators)
julia> Uniform()
Uniform{Int64}(a = 0, b = 1)
```
See also: [`Normal`](@ref), [`TruncatedNormal`](@ref), [`Weibull`](@ref)
"""
struct Uniform{T} <: AbstractDistribution
    a::T
    b::T
end

Uniform(a::T1=0, b::T2=1) where {T1<:Real, T2<:Real} =
check_finite(Uniform, a, "a") &&
check_finite(Uniform, b, "b") &&
check_ordered(Uniform, a, b, "a", "b") &&
Uniform{promote_type(T1, T2)}(a, b)

show(io::IO, distr::Uniform{T}) where T = print(io, string("Uniform{", T, "}(a = ", distr.a, ", b = ", distr.b, ")"))

"""
    transform(D::AbstractDistribution, x::Real)

Apply a transformation to the uniformly distributed number `x` such that it is distributed according to `D`.

The transform is implemented by applying the inverse CDF of the distribution `D`.

!!! note

    This is currently limited only to distributions with an analytic expression for the inverse CDF.

# Examples
```jldoctest; setup = :(using MultilevelEstimators; import Random; Random.seed!(1))
julia> x = rand(10)
10-element Array{Float64,1}:
 0.23603334566204692
 0.34651701419196046
 0.3127069683360675
 0.00790928339056074
 0.4886128300795012
 0.21096820215853596
 0.951916339835734
 0.9999046588986136
 0.25166218303197185
 0.9866663668987996

julia> broadcast(i -> transform(Normal(), i), x)
10-element Array{Float64,1}:
 -0.7191204773101185
 -0.39474101736576256
 -0.4881918784059056
 -2.413074939355834
 -0.02854727903738184
 -0.8030663180302484
  1.6637253880181684
  3.7310515426407296
 -0.6692682660790877
  2.2163540181904553
```
"""
transform(T::Uniform, x::Number) = T.a + (T.b - T.a)*x

## Normal ##
"""
    Normal([μ = 0, [σ = 1]])

Return a normal distribution with mean `μ` and standard deviation `σ`.

# Examples
```jldoctest; setup = :(using MultilevelEstimators)
julia> Normal()
Normal{Int64}(μ = 0, σ = 1)
```
See also: [`Uniform`](@ref), [`TruncatedNormal`](@ref), [`Weibull`](@ref)
"""
struct Normal{T} <: AbstractDistribution
    μ::T
    σ::T
end

Normal(μ::T1=0, σ::T2=1) where {T1<:Real, T2<:Real} =
check_finite(Normal, μ, "μ") &&
check_finite(Normal, σ, "σ") &&
check_larger_than(Normal, σ, "σ", 0) &&
Normal{promote_type(T1, T2)}(μ, σ)

show(io::IO, distr::Normal{T}) where T = print(io, string("Normal{", T, "}(μ = ", distr.μ, ", σ = ", distr.σ, ")"))

transform(T::Normal, x::Number) = T.μ + T.σ*Φ⁻¹(x)

## TruncatedNormal ##
"""
    TruncatedNormal([μ = 0, [σ = 1, [a = -2, [b = 2]]]])

Return a truncated normal distribution with mean `μ` and standard deviation `σ`.

# Examples
```jldoctest; setup = :(using MultilevelEstimators)
julia> TruncatedNormal()
TruncatedNormal{Int64}(μ = 0, σ = 1, a = -2, b = 2)
```
See also: [`Uniform`](@ref), [`Normal`](@ref), [`Weibull`](@ref)
"""
struct TruncatedNormal{T} <: AbstractDistribution
    μ::T
    σ::T
    a::T
    b::T
end

TruncatedNormal(μ::T1=0, σ::T2=1, a::T3=-2, b::T4=2) where {T1<:Real, T2<:Real, T3<:Real, T4<:Real} =
check_finite(TruncatedNormal, μ, "μ") &&
check_finite(TruncatedNormal, σ, "σ") &&
check_finite(TruncatedNormal, a, "a") &&
check_finite(TruncatedNormal, b, "b") &&
check_larger_than(TruncatedNormal, σ, "σ", 0) &&
check_ordered(TruncatedNormal, a, b, "a", "b") &&
TruncatedNormal{promote_type(T1, T2, T3, T4)}(μ, σ, a, b)

show(io::IO, distr::TruncatedNormal{T}) where T = print(io, string("TruncatedNormal{", T, "}(μ = ", distr.μ, ", σ = ", distr.σ, ", a = ", distr.a, ", b = ", distr.b, ")"))

function transform(T::TruncatedNormal, x::Number)
    α = (T.a - T.μ)/T.σ
    β = (T.b - T.μ)/T.σ
    T.μ + T.σ*Φ⁻¹(Φ(α) + (Φ(β) - Φ(α))*x)
end

Φ⁻¹(x::T where {T<:Real}) = √2*erfinv(2*x-1)

Φ(x::T where {T<:Real}) = 1/2*(1+erf(x/√2))

## Weibull ##
"""
    Weibull([k = 2, [λ = √2]])

Return a Weibull distribution with shape parameter `k` and scale parameter `λ`.

# Examples
```jldoctest; setup = :(using MultilevelEstimators)
julia> Weibull()
Weibull{Float64}(k = 2.0, λ = 1.4142135623730951)
```
See also: [`Uniform`](@ref), [`Normal`](@ref), [`TruncatedNormal`](@ref)
"""
struct Weibull{T} <: AbstractDistribution
    k::T
    λ::T
end

Weibull(k::T1=2, λ::T2=√2) where {T1<:Real, T2<:Real} =
check_finite(Weibull, k, "k") &&
check_finite(Weibull, λ, "λ") &&
check_larger_than(Weibull, k, "k", 0) &&
check_larger_than(Weibull, λ, "λ", 0) &&
Weibull{promote_type(T1, T2)}(k, λ)

show(io::IO, distr::Weibull{T}) where T = print(io, string("Weibull{", T, "}(k = ", distr.k, ", λ = ", distr.λ, ")"))

transform(T::Weibull, x::Number) = T.λ * ( -log(1 - x))^(1/T.k)


text:
module BitSAD

using DataStructures
using UUIDs, PearsonHash
using Random: MersenneTwister
using MacroTools
using Reexport

export AbstractBit, AbstractBitstream
export SBit, SBitstream
export DBit, DBitstream
export float, pos, neg
export zero, one
export +, -, *, /, ÷, sqrt, decorrelate, norm
export generate, generate!, estimate!
export push!, pop!, observe, length
export SDM

include("types/bitstream.jl")
include("types/sbitstream.jl")
include("types/dbitstream.jl")
include("modules/sdm.jl")

include("hardware/HW.jl")

@reexport using .HW

end # module


text:
struct DataClientError <: Exception
  msg::String
end

struct ContentParserError <: Exception
  msg::String
end

struct APIResponseError <: Exception
  msg::String
end

struct RequestBuilderError <: Exception
  msg::String
end

struct APIRateExceeded <: Exception
  msg::String
end

struct SchemaValidationError <: Exception
  msg::String
end


text:
export Reaction

"""
A [`Message`](@ref) reaction.
More details [here](https://discordapp.com/developers/docs/resources/channel#reaction-object).
"""
struct Reaction
    count::Int
    me::Bool
    emoji::Emoji
end
@boilerplate Reaction :dict :docs :lower :merge


text:
using Parameters
using Base.Test
using Compat
import Compat.String

# parameters.jl
###############

"Test documentation"
@with_kw immutable MT1
    r::Int = 4
    c = "sdaf"
end
MT1()
@test "Test documentation\n" == Markdown.plain(@doc MT1)

# parameter-less
@with_kw type MT2
    r::Int
    c
    a::Float64
end
MT2(r=4, a=5., c=6)
MT2(r=4, a=5, c="asdf")
MT2(4, "dsaf", 5)
@test_throws ErrorException MT2(r=4)
@test_throws ErrorException MT2()

@with_kw type MT3
    r::Int=5
    a::Float64
end
MT3(r=4, a=5.)
MT3(r=4, a=5)
MT3(a=5)
MT3(4,5)
@test_throws ErrorException MT3(r=4)
@test_throws ErrorException MT3()

# parameter-less with supertype
abstract MM1
abstract MM2{T}
@with_kw type MT3_1 <: MM1
    r::Int=5
    a::Float64
end
MT3_1(r=4, a=5.)
MT3_1(r=4, a=5)
MT3_1(a=5)
MT3_1(4,5)
@test_throws ErrorException MT3_1(r=4)
@test_throws ErrorException MT3_1()

@with_kw type MT3_2 <: MM2{Int}
    r::Int=5
    a::Float64
end
MT3_2(r=4, a=5.)
MT3_2(r=4, a=5)
MT3_2(a=5)
MT3_2(4,5)
@test_throws ErrorException MT3_2(r=4)
@test_throws ErrorException MT3_2()

# with type-parameters
@with_kw immutable MT4{R,I}
    r::R=5
    a::I
end
@test MT4(r=4, a=5.0)==MT4{Int,Float64}(4,5.0)
MT4{Float32, Int}(r=4, a=5.)
MT4{Float32, Int}(a=5.)
MT4{Float32, Int}(5.4, 4)  # inner positional
mt4=MT4(5.4, 4) # outer positional
@test MT4(mt4)==mt4 # outer reconstruct
@test MT4(mt4; a=77)==MT4(5.4, 77)
@test_throws ErrorException MT4{Float32, Int}()
@test_throws InexactError MT4{Float32,Int}(a=5.5)
@test_throws InexactError MT4{Float32,Int}(5.5, 5.5)

# with type-parameters 2
abstract AMT{R<:Real}
@with_kw type MT5{R,I<:Integer} <: AMT{R}
    r::R=5
    a::I
end
@test_throws MethodError MT5(r=4, a=5.) # need to specify type parameters
MT5{Float32, Int}(r=4, a=5.)
MT5{Float32, Int}(a=5.)
MT5{Float32, Int}(5.4, 4)  # inner positional
mt5=MT5(5.4, 4) # outer positional
@test MT5(mt5).r==mt5.r # outer reconstruct
@test MT5(mt5).a==mt5.a # outer reconstruct
@test MT5(mt5; a=77).a==MT5(5.4, 77).a
@test MT5(mt5; a=77).r==MT5(5.4, 77).r
@test_throws ErrorException MT5{Float32, Int}()
@test_throws InexactError MT5{Float32,Int}(a=5.5)
@test_throws InexactError MT5{Float32,Int}(5.5, 5.5)
@test_throws  MethodError MT5(5., "asdf")
@test_throws  TypeError MT5( "asdf", 5)
@test_throws  TypeError MT5{Float64, @compat String}(5., "asdf")
@test_throws  TypeError MT5{@compat String, Int}("asdf", 6)

# with type parameters and supertype
@with_kw type MT4_1{T} <: MM1
    r::Int=5
    a::T
end
MT4_1{Float64}(r=4, a=5.)
MT4_1{Float64}(r=4, a=5)
MT4_1{Float64}(a=5)
MT4_1{Float64}(4,5)
@test_throws ErrorException MT4_1{Float64}(r=4)
@test_throws ErrorException MT4_1{Float64}()

@with_kw type MT4_2{T} <: MM2{T}
    r::Int=5
    a::T
end
MT4_2{Float64}(r=4, a=5.)
MT4_2{Float64}(r=4, a=5)
MT4_2{Float64}(a=5)
MT4_2{Float64}(4,5)
@test_throws ErrorException MT4_2{Float64}(r=4)
@test_throws ErrorException MT4_2{Float64}()

# user defined inner positional constructor
@with_kw immutable MT6{R,I<:Integer} <: AMT{R}
    r::R=5
    a::I
    MT6(r) = new(r,r)
    MT6(r,a) = (@assert a>r; new(r,a))
end
@test_throws MethodError MT6(r=4, a=5.) # need to specify type parameters
MT6{Float32, Int}(r=4, a=5.)
MT6{Float32, Int}(a=6.)
MT6{Float32, Int}(5.4, 6)  # inner positional
mt6=MT6(5.4, 6) # outer positional
@test MT6(mt6)==mt6 # outer reconstruct
@test MT6(mt6; a=77)==MT6(5.4, 77)
@test_throws ErrorException MT6{Float32, Int}()
@test_throws InexactError MT6{Float32,Int}(a=5.5)
@test_throws InexactError MT6{Float32,Int}(5.5, 6.5)
@test_throws  MethodError MT6(5., "asdf")
@test_throws  TypeError MT6( "asdf", 5)
@test_throws  TypeError MT6{Float64, @compat String}(5., "asdf")
@test_throws  TypeError MT6{@compat String, Int}("asdf", 6)

# user defined BAD inner positional constructor
@with_kw type MT7{R,I<:Integer} <: AMT{R}
    r::R=5
    a::I
    MT7(r) = new(r,r+8)
    # no MT7(r,a)
end
@test_throws MethodError MT7{Float32, Int}(r=4, a=5.)

# user defined BAD inner positional constructor
@test_throws ErrorException Parameters.with_kw(:(immutable MT8{R,I<:Integer} <: AMT{R}
    r::R=5
    a::I
    MT8() = new(5,6) # this would shadow the keyword constructor!
    MT8(r,a) = new(r,a)
end))
@test_throws ErrorException Parameters.with_kw(:(type MT8{R,I<:Integer} <: AMT{R}
    r::R=5
    a::I
    MT8(;a=7) = new(5,a) # this would shadow the keyword constructor!
    MT8(r,a) = new(r,a)
end))

# default type annotation (adapted from MT6 test above)
@with_kw immutable MT8{R,I<:Integer} <: AMT{R} @deftype R
    r=5
    a::I
    MT8(r) = new(r,r)
    MT8(r,a) = (@assert a>r; new(r,a))
end
@test_throws MethodError MT8(r=4, a=5.) # need to specify type parameters
MT8{Float32, Int}(r=4, a=5.)
MT8{Float32, Int}(a=6.)
MT8{Float32, Int}(5.4, 6)  # inner positional
mt6=MT8(5.4, 6) # outer positional
@test MT8(mt6)==mt6 # outer reconstruct
@test MT8(mt6; a=77)==MT8(5.4, 77)
@test_throws ErrorException MT8{Float32, Int}()
@test_throws InexactError MT8{Float32,Int}(a=5.5)
@test_throws InexactError MT8{Float32,Int}(5.5, 6.5)
@test_throws  MethodError MT8(5., "asdf")
@test_throws  TypeError MT8( "asdf", 5)
@test_throws  TypeError MT8{Float64, @compat String}(5., "asdf")
@test_throws  TypeError MT8{@compat String, Int}("asdf", 6)
@test MT8.types[1].name==:R
@test MT8.types[2].name==:I
@test MT8.types[2].ub==Integer
@test MT8{Float32,Int32}.types[1]==Float32
@test MT8{Float32,Int32}.types[2]==Int32

# parameter interdependence
@with_kw immutable MT9{R<:Real}
    a::R = 5
    b::R
    c::R = a+b
end
@test MT9{Float64}(b=1).c==6
@test MT9{Float64}(b=1, c=1).c==1

@with_kw immutable MT10{R<:Real}
    b::R = 6
    c::R = a+b
    a::R = 5
end
@test_throws UndefVarError MT10{Float64}() # defaults are evaluated in order
@test MT10{Float64}(b=1).c==6  # this shouldn't work but does: https://github.com/JuliaLang/julia/issues/9535#issuecomment-73717708
@test MT10{Float64}(b=1, c=1).c==1

# binding outside variables
a = 7
b = [1,2]
@with_kw immutable MT11
    aa::Int=a  # a::Int=a is not possible as the outside a gets shadowed
    bb::Vector{Int}=b
end
m = MT11()
@test m.aa===a
@test m.bb===b
a = 5
b[1] = 2
@test m.aa===7
@test m.bb[1]==2

## (Un)pack
@with_kw immutable P1
    r::Int
    c
    a::Float64
end

let
    mt = P1(r=4, a=5, c=6)
    @unpack_P1 mt
    @test r===4
    @test c===6
    @test a===5.
    r = 1
    a = 2
    c = 3
    @pack_P1 mt
    @test mt===P1(r=1, a=2, c=3)
    @test string(mt) == "P1\n  r: 1\n  c: 3\n  a: 2.0"
end

### Assertions
@with_kw immutable MT12
    a=5; @assert a>=5
    b
    @assert b>a
end

@test_throws AssertionError MT12(b=2)
@test_throws AssertionError MT12(a=1,b=2)
@test MT12(b=6)==MT12(5,6)

# only asserts allowed if no inner constructors
@test_throws ErrorException Parameters.with_kw(:(immutable MT13
                                                 a=5; @assert a>=5
                                                 MT13(a) = new(8)
                                                 end))

####
# issue 10: infer type parameters from kw-args

@with_kw immutable I10{T} @deftype Int
    a::T
    b = 10
    c::T="aaa"
end
@test_throws ErrorException I10()
@test I10(1,2,3)==I10{Int}(1,2,3)
@test_throws MethodError I10(a=10) # typeof(a)!=typeof(c)
a =  I10(a="asd")
b = I10{String}("asd",10,"aaa")
for fn in fieldnames(a)
    # complicated testing because of mutable T
    @test getfield(a, fn)==getfield(b, fn)
end

# issue #12: only one-liner inner constructors were parsed correctly.
@with_kw immutable T9867
    r::Int
    function T9867(r)
        new(r)
    end
end
@test_throws ErrorException T9867()
@test T9867(r=2).r == 2

###########################
# Packing and unpacking @unpack, @pack
##########################
# Example with dict:
d = Dict{Symbol,Any}(:a=>5.0,:b=>2,:c=>"Hi!")
@unpack a, c = d
@test a == 5.0 #true
@test c == "Hi!" #true

# Example with type:
type A; a; b; c; end
d = A(4,7.0,"Hi!")
@unpack a, c = d
@test a == 4 #true
@test c == "Hi!" #true

## Packing

# Example with dict:
a = 5.0
c = "Hi!"
d = Dict{Symbol,Any}()
@pack d = a, c
@test d==Dict{Symbol,Any}(:a=>5.0,:c=>"Hi!")

# Example with type:
a = 99
c = "HaHa"
d = A(4,7.0,"Hi")
@pack d = a, c
@test d.a == 99
@test d.c == "HaHa"

# older tests ported
type UP1
    a
    b
end
uu = UP1(1,2)
@test_throws ErrorException @unpack c = uu
@test_throws ErrorException @unpack a, c = uu

a, b = 0, 0
@unpack a = uu
@test a==1
@test b==0
a, b = 0, 0
@unpack a, b = uu
@test a==1
@test b==2

vv = uu
a = 99
@pack uu = a
@test uu==vv
@test uu.a==99

immutable UP2
    a
    b
end
uu = UP2(1,2)
@test_throws ErrorException @unpack c = uu
@test_throws ErrorException @unpack a, c = uu

a, b = 0, 0
@unpack a = uu
@test a==1
@test b==0
a, b = 0, 0
@unpack a,b = uu
@test a==1
@test b==2

a = 99
@test_throws ErrorException @pack uu = a

# check that inference works
immutable UP3
    a::Float64
    b::Int
end
function f(u::UP3)
    @unpack a,b = u
    a,b
end
if VERSION>v"0.5-"
    @inferred f(UP3(1,2))
end


text:
@testset "tosecond" begin
    @test_throws MethodError Wflow.tosecond(Month(2))
    @test Wflow.tosecond(Week(2)) == 86400 * 14
    @test Wflow.tosecond(Day(2)) == 86400 * 2
    @test Wflow.tosecond(Hour(2)) == 3600 * 2
    @test Wflow.tosecond(Minute(2)) == 60 * 2
    @test Wflow.tosecond(Second(2)) == 2
    @test Wflow.tosecond(Millisecond(2)) == 2e-3
    @test Wflow.tosecond(Microsecond(2)) == 2e-6
    @test Wflow.tosecond(Nanosecond(2)) == 2e-9
end


text:
module ReportMetrics

import Pkg
import PrettyTables
import Coverage

"""
    report_allocs(;
        job_name::String,
        run_cmd::String = "",
        deps_to_monitor::Vector{Module} = Module[],
        dirs_to_monitor::Vector{String} = String[],
        pkg_name::Union{Nothing, String} = nothing,
        n_unique_allocs::Int = 10,
        write_csv::Bool = false,
        csv_prefix_path::String = "",
    )

Reports allocations
 - `job_name` name of job
 - `run_cmd` a `Base.Cmd` to run script
 - `deps_to_monitor` a `Vector` of modules to monitor
 - `dirs_to_monitor` a `Vector` of directories to monitor
 - `pkg_name` name of package being tested (for helping with string processing)
 - `n_unique_allocs` limits number of unique allocation sites to report (to avoid large tables)
 - `suppress_url` (` = true`) suppress trying to use URLs in the output table
 - `write_csv` Bool to write a csv output
 - `csv_prefix_path` prefix path to write the csv file

## Notest
 - `deps_to_monitor` and `dirs_to_monitor` are merged together.
"""
function report_allocs(;
        job_name::String,
        run_cmd::Union{Nothing, Base.Cmd} = nothing,
        deps_to_monitor::Vector{Module} = Module[],
        dirs_to_monitor::Vector{String} = String[],
        process_filename::Function = process_filename_default,
        is_loading_pkg::Function = (fn, ln) -> false,
        pkg_name::Union{Nothing, String} = nothing,
        n_unique_allocs::Int = 10,
        suppress_url::Bool = true,
        write_csv::Bool = false,
        csv_prefix_path::String = "",
        format_locid::Function = x -> x,
    )

    ##### Collect deps
    dep_dirs = map(deps_to_monitor) do dep
        pkgdir(dep)
    end
    all_dirs_to_monitor = [dirs_to_monitor..., dep_dirs...]

    ##### Run representative work & track allocations
    local allocs
    try
        run(run_cmd)
        allocs = Coverage.analyze_malloc(all_dirs_to_monitor)
    finally

        ##### Clean up files
        for d in all_dirs_to_monitor
            all_files = [
                joinpath(root, f) for
                (root, dirs, files) in Base.Filesystem.walkdir(d) for f in files
            ]
            all_mem_files = filter(x -> endswith(x, ".mem"), all_files)
            for f in all_mem_files
                rm(f)
            end
        end
    end

    ##### Process and plot results
    filter!(x -> x.bytes ≠ 0, allocs)
    filter!(x -> !is_loading_pkg(x.filename, x.linenumber), allocs) # Try to skip loading module if pkg_name is included

    n_alloc_sites = length(allocs)
    if n_alloc_sites == 0
        @info "$(job_name): Number of sampled allocating sites: $n_alloc_sites 🎉"
        return
    end

    all_bytes = reverse(getproperty.(allocs, :bytes))
    all_filenames = reverse(getproperty.(allocs, :filename))
    all_linenumbers = reverse(getproperty.(allocs, :linenumber))
    process_fn(fn) = post_process_fn(process_filename(fn))

    bytes_subset = Int[]
    filenames_subset = String[]
    linenumbers_subset = Int[]
    loc_ids_subset = String[]
    truncated_allocs = false
    for (bytes, filename, linenumber) in zip(all_bytes, all_filenames, all_linenumbers)
        is_loading_pkg(filename, linenumber) && continue # Try to skip loading module if pkg_name is included
        loc_id = "$(process_fn(filename)):$(linenumber)"
        if !(bytes in bytes_subset) && !(loc_id in loc_ids_subset)
            push!(bytes_subset, bytes)
            push!(filenames_subset, filename)
            push!(linenumbers_subset, linenumber)
            push!(loc_ids_subset, loc_id)
            if length(bytes_subset) ≥ n_unique_allocs
                truncated_allocs = true
                break
            end
        end
    end
    sum_bytes = sum(bytes_subset)
    trunc_msg = truncated_allocs ? " (truncated) " : ""
    @info "$(job_name): $(length(bytes_subset)) unique allocating sites, $sum_bytes total bytes$trunc_msg"
    xtick_name(filename, linenumber) = "$filename:$linenumber"
    labels = xtick_name.(process_fn.(filenames_subset), linenumbers_subset)

    # TODO: get urls for hypertext
    pkg_urls = Dict(map(all_dirs_to_monitor) do dep_dir
        proj = Pkg.Types.read_project(joinpath(dep_dir, "Project.toml"))
        if proj.uuid ≠ nothing
            url = Pkg.Operations.find_urls(Pkg.Types.Context().registries, proj.uuid)
            Pair(proj.name, url)
        else
            Pair(proj.name, "https://www.google.com")
        end
    end...)

    fileinfo = map(zip(filenames_subset, linenumbers_subset)) do (filename, linenumber)
        label = xtick_name(process_fn(filename), linenumber)
        if suppress_url
            label
        else
            url = ""
            # name = basename(pkg_dir_from_file(dirname(filename)))
            # TODO: incorporate URLS into table
            # if haskey(pkg_urls, name)
            #     url = pkg_urls[name]
            # else
            #     url = "https://www.google.com"
            # end
            PrettyTables.URLTextCell(label, url)
        end
    end

    alloc_percent = map(bytes_subset) do bytes
        alloc_perc = bytes / sum_bytes
        Int(round(alloc_perc*100, digits = 0))
    end
    header = (
        ["<file>:<line number>", "Allocations", "Allocations %"],
        ["", "(bytes)", "(xᵢ/∑x)"],
    )

    table_data = hcat(
        fileinfo,
        bytes_subset,
        alloc_percent,
    )

    PrettyTables.pretty_table(
        table_data;
        header,
        formatters = PrettyTables.ft_printf("%s", 2:2),
        header_crayon = PrettyTables.crayon"yellow bold",
        subheader_crayon = PrettyTables.crayon"green bold",
        crop = :none,
        alignment = [:l, :c, :c],
    )
    if write_csv
        mkpath(csv_prefix_path)
        open(joinpath(csv_prefix_path, "$job_name" * ".csv"), "w") do fh
            for (bytes, loc) in zip(bytes_subset, loc_ids_subset)
                println(fh, "$bytes\t|\t$(format_locid(loc))")
                # println(fh, "$bytes,$loc")
            end
        end
    end
end

# Try to find a package directory from a file
function pkg_dir_from_file(filename, candidates = String[])
    _pkg_dir_from_file!(dirname(filename), candidates)
    if length(candidates) > 1
        @debug "Multiple Project.toml files found recursively through $filename.\n$(first(candidates)) used."
    end
    return first(candidates)
end

function _pkg_dir_from_file!(dir, candidates)
    @assert ispath(dir)
    if isfile(joinpath(dir, "Project.toml"))
        push!(candidates, dir)
        _pkg_dir_from_file!(dirname(dir), candidates)
    end
end

function post_process_fn(fn)
    # Remove ###.mem.
    fn = join(split(fn, ".jl")[1:(end - 1)], ".jl") * ".jl"
    if startswith(fn, Base.Filesystem.path_separator)
        fn = fn[2:end]
    end
    return fn
end

function process_filename_default(fn)
    # TODO: make this work for Windows
    if occursin(".julia/packages/", fn)
        fn = last(split(fn, ".julia/packages/"))
        pkg_name = first(split(fn, "/"))
        if occursin("$pkg_name/src", fn)
            return fn
        else
            fn = join(split(fn, pkg_name)[2:end], pkg_name)
            sha = split(fn, "/")[2]
            fn = replace(fn, sha*"/" => "")
            fn = pkg_name*fn
            return fn
        end
    end
    return fn
end

#=
# name helper functions
function filename_only(fn)
    if occursin(".jl", fn)
        fn = join(split(fn, ".jl")[1:(end - 1)], ".jl") * ".jl"
    end
    fn = replace(fn, "\\" => "/") # for windows...
    isempty(deps_to_monitor) && return fn

    for dep_name in string.(deps_to_monitor)
        if occursin(dep_name, fn)
            fn = dep_name * last(split(fn, dep_name))
        end
    end
    isempty(splitbys) && return fn
    for splitby in splitbys
        if occursin(splitby, fn)
            fn = splitby * last(split(fn, splitby))
        end
    end
    return fn
end
function compile_pkg(fn, linenumber)
    c1 = if pkg_name ≠ nothing
        endswith(filename_only(fn), pkg_name)
    else
        false
    end
    c2 = linenumber == 1
    return c1 && c2
end
=#

#=
For integrating with buildkite

function format_locid(locid)
    relpth, lineno = split(locid, ':')
    if haskey(ENV, "BUILDKITE") && ENV["BUILDKITE"] == "true"
        return join(
            [PKG_ORG, PKG_NAME, "blob", ENV["BUILDKITE_COMMIT"], relpth],
            '/',
        ) * "#$(lineno)"
    end
    return join([PKG_NAME, locid], '/')
end
=#

import SnoopCompile

"""
    report_invalidations(;
        job_name::String,
        invalidations,
        n_rows::Int = 10,
        process_filename::Function = x -> x,
        write_csv::Bool = false,
        csv_prefix_path::String = "",
        format_locid::Function = x -> x,
    )

Report a table of invalidations.

For documentation, see https://timholy.github.io/SnoopCompile.jl/stable/snoopr/

## Usage example
```julia
import SnoopCompileCore
invalidations = SnoopCompileCore.@snoopr begin

    # load packages & do representative work

end;
import ReportMetrics
ReportMetrics.report_invalidations(;
    job_name = "MyWork",
    invalidations,
    process_filename = x -> last(split(x, "packages/")),
)
```
"""
function report_invalidations(;
        job_name::String,
        invalidations,
        n_rows::Int = 10,
        process_filename::Function = x -> x,
        write_csv::Bool = false,
        csv_prefix_path::String = "",
        format_locid::Function = x -> x,
    )

    trees = reverse(SnoopCompile.invalidation_trees(invalidations))

    n_total_invalidations = length(SnoopCompile.uinvalidated(invalidations))
    @info "Number of invalidations for $job_name: $n_total_invalidations"

    invs_per_method = map(trees) do methinvs
        SnoopCompile.countchildren(methinvs)
    end
    n_invs_total = length(invs_per_method)
    truncated_invs = n_rows < n_invs_total
    sum_invs = sum(invs_per_method)
    invs_per_method = invs_per_method[1:min(n_rows, n_invs_total)]
    trees = trees[1:min(n_rows, n_invs_total)]
    trunc_msg = truncated_invs ? " (showing $n_rows) " : ""
    @info "$(job_name): $n_invs_total total method invalidations$trunc_msg"

    n_invalidations_percent = map(invs_per_method) do inv
        inv_perc = inv / sum_invs
        Int(round(inv_perc*100, digits = 0))
    end
    meth_name = map(trees) do inv
        "$(inv.method.name)"
    end
    fileinfo = map(trees) do inv
        "$(process_filename(string(inv.method.file))):$(inv.method.line)"
    end

    header = (
        ["<file name>:<line number>", "Method Name", "Invalidations", "Invalidations %"],
        ["", "", "Number", "(xᵢ/∑x)"],
    )

    table_data = hcat(
        fileinfo,
        meth_name,
        invs_per_method,
        n_invalidations_percent,
    )

    PrettyTables.pretty_table(
        table_data;
        header,
        formatters = PrettyTables.ft_printf("%s", 2:2),
        header_crayon = PrettyTables.crayon"yellow bold",
        subheader_crayon = PrettyTables.crayon"green bold",
        crop = :none,
        alignment = [:l, :c, :c, :c],
    )
    if write_csv
        mkpath(csv_prefix_path)
        open(joinpath(csv_prefix_path, "$job_name" * ".csv"), "w") do fh
            for (name, finfo, inv) in zip(meth_name, fileinfo, invs_per_method)
                println(fh, "$(format_locid(finfo))\t|\t$name\t|\t$inv")
            end
        end
    end

end

end # module


text:
using Base.Test

N = 10

C = rand(N, N)
A = C - C'

T, Z, λ = schur(A)
@test Z*T*Z' ≈ A

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

A = Antisymmetric(A)
T, Z, λ = schur(A)

@test Z*T*Z' ≈ A

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

A = AntisymTridiagonal(rand(N))
@show A

T, Z, λ = schur(A)
@test Z*T*Z' ≈ A


text:
export DefaultStateStyleEnv

struct DefaultStateStyleEnv{S,E} <: AbstractEnvWrapper
    env::E
end

"""
    DefaultStateStyleEnv{S}(env::E)

Reset the result of `DefaultStateStyle` without changing the original behavior.
"""
DefaultStateStyleEnv{S}(env::E) where {S,E} = DefaultStateStyleEnv{S,E}(env)

RLBase.DefaultStateStyle(::DefaultStateStyleEnv{S}) where {S} = S

RLBase.state(env::DefaultStateStyleEnv{S}) where S = state(env.env, S)
RLBase.state(env::DefaultStateStyleEnv, ss::RLBase.AbstractStateStyle) = state(env.env, ss)
RLBase.state(env::DefaultStateStyleEnv{S}, player) where S = state(env.env, S, player)

RLBase.state_space(env::DefaultStateStyleEnv{S}) where S = state_space(env.env, S)
RLBase.state_space(env::DefaultStateStyleEnv, ss::RLBase.AbstractStateStyle) = state_space(env.env, ss)


text:
using Thebes, Luxor

@draw begin
    background("black")
    setopacity(0.5)
    sethue("gold")
    c = pin.([Point3D(randn() * 50, randn() * 50, randn() * 50) for x in 1:75, y in 1:75, z in 1:75], gfunction = (pt3, pt2) -> begin
        circle(pt2, 1, :fill)
    end)
    axes3D()
    @show length(c)
end


text:
/* Testing that main is allowed not to return. */

int main() {
        return 0 ;

}

max_stars_repo_name: lukaszcz/javalette

max_issues_repo_name: lukaszcz/javalette

max_forks_repo_name: lukaszcz/javalette


text:
module RegisterMismatchCommon

using RegisterCore, CenterIndexedArrays, ImageCore

export correctbias!, nanpad, mismatch0, aperture_grid, allocate_mmarrays, default_aperture_width, truncatenoise!
export DimsLike, WidthLike, each_point, aperture_range, assertsamesize, tovec, mismatch, padsize, set_FFTPROD
export padranges, shiftrange, checksize_maxshift, register_translate

const DimsLike = Union{AbstractVector{Int}, Dims}   # try to avoid these and just use Dims tuples for sizes
const WidthLike = Union{AbstractVector,Tuple}
FFTPROD = [2,3]
set_FFTPROD(v) = global FFTPROD = v

mismatch(fixed::AbstractArray{T}, moving::AbstractArray{T}, maxshift::DimsLike; normalization = :intensity) where {T<:AbstractFloat} = mismatch(T, fixed, moving, maxshift; normalization=normalization)
mismatch(fixed::AbstractArray, moving::AbstractArray, maxshift::DimsLike; normalization = :intensity) = mismatch(Float32, fixed, moving, maxshift; normalization=normalization)

mismatch_apertures(fixed::AbstractArray{T}, moving::AbstractArray{T}, args...; kwargs...) where {T<:AbstractFloat} = mismatch_apertures(T, fixed, moving, args...; kwargs...)
mismatch_apertures(fixed::AbstractArray, moving::AbstractArray, args...; kwargs...) = mismatch_apertures(Float32, fixed, moving, args...; kwargs...)

function mismatch_apertures(::Type{T}, fixed::AbstractArray, moving::AbstractArray, gridsize::DimsLike, maxshift::DimsLike; kwargs...) where T
    cs = coords_spatial(fixed)
    aperture_centers = aperture_grid(map(d->size(fixed, d), cs), gridsize)
    aperture_width = default_aperture_width(fixed, gridsize)
    mismatch_apertures(T, fixed, moving, aperture_centers, aperture_width, maxshift; kwargs...)
end

"""
`correctbias!(mm::MismatchArray)` replaces "suspect" mismatch
data with imputed data.  If each pixel in your camera has a different
bias, then matching that bias becomes an incentive to avoid
shifts.  Likewise, CMOS cameras tend to have correlated row/column
noise. These two factors combine to imply that `mm[i,j,...]` is unreliable
whenever `i` or `j` is zero.

Data are imputed by averaging the adjacent non-suspect values.  This
function works in-place, overwriting the original `mm`.
"""
function correctbias!(mm::MismatchArray{ND,N}, w = correctbias_weight(mm)) where {ND,N}
    T = eltype(ND)
    mxshift = maxshift(mm)
    Imax = CartesianIndex(mxshift)
    Imin = CartesianIndex(map(x->-x,mxshift)::NTuple{N,Int})
    I1 = CartesianIndex(ntuple(d->d>2 ? 0 : 1, N)::NTuple{N,Int})  # only first 2 dims
    for I in eachindex(mm)
        if w[I] == 0
            mms = NumDenom{T}(0,0)
            ws = zero(T)
            strt = max(Imin, I-I1)
            stop = min(Imax, I+I1)
            for J in CartesianIndices(ntuple(d->strt[d]:stop[d],N))
                wJ = w[J]
                if wJ != 0
                    mms += wJ*mm[J]
                    ws += wJ
                end
            end
            mm[I] = mms/ws
        end
    end
    mm
end

"`correctbias!(mms)` runs `correctbias!` on each element of an array-of-MismatchArrays."
function correctbias!(mms::AbstractArray{M}) where M<:MismatchArray
    for mm in mms
        correctbias!(mm)
    end
    mms
end

function correctbias_weight(mm::MismatchArray{ND,N}) where {ND,N}
    T = eltype(ND)
    w = CenterIndexedArray(ones(T, size(mm)))
    for I in eachindex(mm)
        anyzero = false
        for d = 1:min(N,2)   # only first 2 dims
            anyzero |= I[d] == 0
        end
        if anyzero
            w[I] = 0
        end
    end
    w
end

"""
`fixedpad, movingpad = nanpad(fixed, moving)` will pad `fixed` and/or
`moving` with NaN as needed to ensure that `fixedpad` and `movingpad`
have the same size.
"""
function nanpad(fixed, moving)
    ndims(fixed) == ndims(moving) || error("fixed and moving must have the same dimensionality")
    if size(fixed) == size(moving)
        return fixed, moving
    end
    rng = map(d->1:max(size(fixed,d), size(moving,d)), 1:ndims(fixed))
    T = promote_type(eltype(fixed), eltype(moving))
    get(fixed, rng, nanval(T)), get(moving, rng, nanval(T))
end

nanval(::Type{T}) where {T<:AbstractFloat} = convert(T, NaN)
nanval(::Type{T}) where {T} = convert(Float32, NaN)

"""
`mm0 = mismatch0(fixed, moving, [normalization])` computes the
"as-is" mismatch between `fixed` and `moving`, without any shift.
`normalization` may be either `:intensity` (the default) or `:pixels`.
"""
function mismatch0(fixed::AbstractArray{Tf,N}, moving::AbstractArray{Tm,N}; normalization = :intensity) where {Tf,Tm,N}
    size(fixed) == size(moving) || throw(DimensionMismatch("Size $(size(fixed)) of fixed is not equal to size $(size(moving)) of moving"))
    _mismatch0(zero(Float64), zero(Float64), fixed, moving; normalization=normalization)
end

function _mismatch0(num::T, denom::T, fixed::AbstractArray{Tf,N}, moving::AbstractArray{Tm,N}; normalization = :intensity) where {T,Tf,Tm,N}
    if normalization == :intensity
        for i in eachindex(fixed, moving)
            vf = T(fixed[i])
            vm = T(moving[i])
            if isfinite(vf) && isfinite(vm)
                num += (vf-vm)^2
                denom += vf^2 + vm^2
            end
        end
    elseif normalization == :pixels
        for i in eachindex(fixed, moving)
            vf = T(fixed[i])
            vm = T(moving[i])
            if isfinite(vf) && isfinite(vm)
                num += (vf-vm)^2
                denom += 1
            end
        end
    else
        error("Normalization $normalization not recognized")
    end
    NumDenom(num, denom)
end

"""
`mm0 = mismatch0(mms)` computes the "as-is"
mismatch between `fixed` and `moving`, without any shift.  The
mismatch is represented in `mms` as an aperture-wise
Arrays-of-MismatchArrays.
"""
function mismatch0(mms::AbstractArray{M}) where M<:MismatchArray
    mm0 = eltype(M)(0, 0)
    cr = eachindex(first(mms))
    z = first(cr)+last(cr)  # all-zeros CartesianIndex
    for mm in mms
        mm0 += mm[z]
    end
    mm0
end

"""
`ag = aperture_grid(ssize, gridsize)` constructs a uniformly-spaced
grid of aperture centers.  The grid has size `gridsize`, and is
constructed for an image of spatial size `ssize`.  Along each
dimension the first and last elements are at the image corners.
"""
function aperture_grid(ssize::Dims{N}, gridsize) where N
    if length(gridsize) != N
        if length(gridsize) == N-1
            @info("ssize and gridsize disagree; possible fix is to use a :time axis (AxisArrays) for the image")
        end
        error("ssize and gridsize must have the same length, got $ssize and $gridsize")
    end
    grid = Array{NTuple{N,Float64},N}(undef, (gridsize...,))
    centers = map(i-> gridsize[i] > 1 ? collect(range(1, stop=ssize[i], length=gridsize[i])) : [(ssize[i]+1)/2], 1:N)
    for I in CartesianIndices(size(grid))
        grid[I] = ntuple(i->centers[i][I[i]], N)
    end
    grid
end

"""
`mms = allocate_mmarrays(T, gridsize, maxshift)` allocates storage for
aperture-wise mismatch computation. `mms` will be an
Array-of-MismatchArrays with element type `NumDenom{T}` and half-size
`maxshift`. `mms` will be an array of size `gridsize`. This syntax is
recommended when your apertures are centered at points of a grid.

`mms = allocate_mmarrays(T, aperture_centers, maxshift)` returns `mms`
with a shape that matches that of `aperture_centers`. The centers can
in general be provided as an vector-of-tuples, vector-of-vectors, or a
matrix with each point in a column.  If your centers are arranged in a
rectangular grid, you can use an `N`-dimensional array-of-tuples (or
array-of-vectors) or an `N+1`-dimensional array with the center
positions specified along the first dimension.  (But you may find the
`gridsize` syntax to be simpler.)
"""
function allocate_mmarrays(::Type{T}, aperture_centers::AbstractArray{C}, maxshift) where {T,C<:Union{AbstractVector,Tuple}}
    isempty(aperture_centers) && error("aperture_centers is empty")
    N = length(first(aperture_centers))
    sz = map(x->2*x+1, maxshift)
    mm = MismatchArray(T, sz...)
    mms = Array{typeof(mm)}(undef, size(aperture_centers))
    f = true
    for i in eachindex(mms)
        if f
            mms[i] = mm
            f = false
        else
            mms[i] = MismatchArray(T, sz...)
        end
    end
    mms
end

function allocate_mmarrays(::Type{T}, aperture_centers::AbstractArray{R}, maxshift) where {T,R<:Real}
    N = ndims(aperture_centers)-1
    mms = Array{MismatchArray{T,N}}(undef, size(aperture_centers)[2:end])
    sz = map(x->2*x+1, maxshift)
    for i in eachindex(mms)
        mms[i] = MismatchArray(T, sz...)
    end
    mms
end

function allocate_mmarrays(::Type{T}, gridsize::NTuple{N,Int}, maxshift) where {T<:Real,N}
    mms = Array{MismatchArray{NumDenom{T},N}}(undef, gridsize)
    sz = map(x->2*x+1, maxshift)
    for i in eachindex(mms)
        mms[i] = MismatchArray(T, sz...)
    end
    mms
end

struct ContainerIterator{C}
    data::C
end

Base.iterate(iter::ContainerIterator) = iterate(iter.data)
Base.iterate(iter::ContainerIterator, state) = iterate(iter.data, state)

struct FirstDimIterator{A<:AbstractArray,R<:CartesianIndices}
    data::A
    rng::R

    FirstDimIterator{A,R}(data::A) where {A,R} = new{A,R}(data, CartesianIndices(Base.tail(size(data))))
end
FirstDimIterator(A::AbstractArray) = FirstDimIterator{typeof(A),typeof(CartesianIndices(Base.tail(size(A))))}(A)
function Base.iterate(iter::FirstDimIterator)
    isempty(iter.rng) && return nothing
    index, state = iterate(iter.rng)
    iter.data[:, index], state
end
function Base.iterate(iter::FirstDimIterator, state)
    state == last(iter.rng) && return nothing
    index, state = iterate(iter.rng, state)
    iter.data[:, index], state
end

"""
`iter = each_point(points)` yields an iterator `iter` over all the
points in `points`. `points` may be represented as an
AbstractArray-of-tuples or -AbstractVectors, or may be an
`AbstractArray` where each point is represented along the first
dimension (e.g., columns of a matrix).
"""
each_point(aperture_centers::AbstractArray{C}) where {C<:Union{AbstractVector,Tuple}} = ContainerIterator(aperture_centers)

each_point(aperture_centers::AbstractArray{R}) where {R<:Real} = FirstDimIterator(aperture_centers)

"""
`rng = aperture_range(center, width)` returns a tuple of
`UnitRange{Int}`s that, for dimension `d`, is centered on `center[d]`
and has width `width[d]`.
"""
function aperture_range(center, width)
    length(center) == length(width) || error("center and width must have the same length")
    ntuple(d->leftedge(center[d], width[d]):rightedge(center[d], width[d]), length(center))
end

"""
`aperturesize = default_aperture_width(img, gridsize, [overlap])`
calculates the aperture width for a regularly-spaced grid of aperture
centers with size `gridsize`.  Apertures that are adjacent along
dimension `d` may overlap by a number pixels specified by
`overlap[d]`; the default value is 0.  For non-negative `overlap`, the
collection of apertures will yield full coverage of the image.
"""
function default_aperture_width(img, gridsize::DimsLike, overlap::DimsLike = zeros(Int, sdims(img)))
    sc = coords_spatial(img)
    length(sc) == length(gridsize) == length(overlap) || error("gridsize and overlap must have length equal to the number of spatial dimensions in img")
    for i = 1:length(sc)
        if gridsize[i] > size(img, sc[i])
            error("gridsize $gridsize is too large, given the size $(size(img)[sc]) of the image")
        end
    end
    gsz1 = max.(1, [gridsize...].-1)
    gflag = [gridsize...].>1
    tuple((([map(d->size(img,d),sc)...]-gflag)./gsz1+2*[overlap...].*gflag)...)
end

"""
`truncatenoise!(mm, thresh)` zeros out any entries of the
MismatchArray `mm` whose `denom` values are less than `thresh`.
"""
function truncatenoise!(mm::AbstractArray{NumDenom{T}}, thresh::Real) where T<:Real
    for I in eachindex(mm)
        if mm[I].denom <= thresh
            mm[I] = NumDenom{T}(0,0)
        end
    end
    mm
end

function truncatenoise!(mms::AbstractArray{A}, thresh::Real) where A<:MismatchArray
    for i = 1:length(denoms)
        truncatenoise!(mms[i], thresh)
    end
    nothing
end

"""
`shift = register_translate(fixed, moving, maxshift, [thresh])`
computes the integer-valued translation which best aligns images
`fixed` and `moving`. All shifts up to size `maxshift` are considered.
Optionally specify `thresh`, the fraction (0<=thresh<=1) of overlap
required between `fixed` and `moving` (default 0.25).
"""
function register_translate(fixed, moving, maxshift, thresh=nothing)
    mm = mismatch(fixed, moving, maxshift)
    _, denom = separate(mm)
    if thresh==nothing
        thresh = 0.25maximum(denom)
    end
    indmin_mismatch(mm, thresh)
end

function checksize_maxshift(A::AbstractArray, maxshift)
    ndims(A) == length(maxshift) || error("Array is $(ndims(A))-dimensional, but maxshift has length $(length(maxshift))")
    for i = 1:ndims(A)
        size(A,i) == 2*maxshift[i]+1 || error("Along dimension $i, the output size $(size(A,i)) does not agree with maxshift[$i] = $(maxshift[i])")
    end
    nothing
end

function padranges(blocksize, maxshift)
    padright = [maxshift...]
    transformdims = findall(padright.>0)
    paddedsz = [blocksize...] + 2*padright
    for i in transformdims
        # Pick a size for which one can efficiently calculate ffts
        padright[i] += padsize(blocksize, maxshift, i) - paddedsz[i]
    end
    rng = UnitRange{Int}[ 1-maxshift[i]:blocksize[i]+padright[i] for i = 1:length(blocksize) ]
end

padsize(blocksize::Dims{N}, maxshift::Dims{N}) where N = map(padsize, blocksize, maxshift, ntuple(identity, Val(N)))

function padsize(blocksize::Int, maxshift::Int, dim::Int)
    p = blocksize + 2maxshift
    return maxshift > 0 ? (dim == 1 ? nextpow(2, p) : nextprod(FFTPROD, p)) : p   # we won't FFT along dimensions with maxshift 0
end

function assertsamesize(A, B)
    if !issamesize(A,B)
        error("Arrays are not the same size")
    end
end

function issamesize(A::AbstractArray, B::AbstractArray)
    n = ndims(A)
    ndims(B) == n || return false
    for i = 1:n
        axes(A, i) == axes(B, i) || return false
    end
    true
end

function issamesize(A::AbstractArray, indices)
    n = ndims(A)
    length(indices) == n || return false
    for i = 1:n
        size(A, i) == length(indices[i]) || return false
    end
    true
end

# This yields the _effective_ overlap, i.e., sets to zero if gridsize==1 along a coordinate
# imgssz = image spatial size
function computeoverlap(imgssz, blocksize, gridsize)
    gsz1 = max(1,[gridsize...].-1)
    tmp = [imgssz...]./gsz1
    blocksize - [ceil(Int, x) for x in tmp]
end

leftedge(center, width) = ceil(Int, center-width/2)
rightedge(center, width) = leftedge(center+width, width) - 1

# These avoid making a copy if it's not necessary
tovec(v::AbstractVector) = v
tovec(v::Tuple) = [v...]

shiftrange(r, s) = r .+ s

### Utilities for unsafe indexing of views
# TODO: redesign this whole thing to be safer?
using Base: ViewIndex, to_indices, unsafe_length, index_shape, tail

@inline function extraunsafe_view(V::SubArray{T,N}, I::Vararg{ViewIndex,N}) where {T,N}
    idxs = unsafe_reindex(V, V.indices, to_indices(V, I))
    SubArray(V.parent, idxs)
end

function get_index_wo_boundcheck(r::AbstractUnitRange, s::AbstractUnitRange{<:Integer}) # BlockRegistration issue #36
    f = first(r)
    st = oftype(f, f + first(s)-1)
    range(st, length=length(s))
end

unsafe_reindex(V, idxs::Tuple{UnitRange, Vararg{Any}}, subidxs::Tuple{UnitRange, Vararg{Any}}) =
    (Base.@_propagate_inbounds_meta; @inbounds new1 = get_index_wo_boundcheck(idxs[1],subidxs[1]);
     (new1, unsafe_reindex(V, tail(idxs), tail(subidxs))...))

unsafe_reindex(V, idxs, subidxs) = Base.reindex(V, idxs, subidxs)

### Deprecations

function padsize(blocksize, maxshift)
    Base.depwarn("padsize(::$(typeof(blocksize)), ::$(typeof(maxshift)) is deprecated, use Dims-tuples instead", :padsize)
    sz = Vector{Int}(undef, length(blocksize))
    padsize!(sz, blocksize, maxshift)
end

function padsize!(sz::Vector, blocksize, maxshift)
    n = length(blocksize)
    for i = 1:n
        sz[i] = padsize(blocksize, maxshift, i)
    end
    sz
end

function padsize(blocksize, maxshift, dim)
    m = maxshift[dim]
    p = blocksize[dim] + 2m
    return m > 0 ? (dim == 1 ? nextpow(2, p) : nextprod(FFTPROD, p)) : p   # we won't FFT along dimensions with maxshift[i]==0
end

end #module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

include("geoops/cat.jl")
include("geoops/split.jl")
include("geoops/slice.jl")
include("geoops/sample.jl")
include("geoops/filter.jl")
include("geoops/groupby.jl")
include("geoops/boundbox.jl")
include("geoops/inside.jl")
include("geoops/unique.jl")


text:
"""
    printExperimentInfo(genj::GeneticJuliaStructure = GenJ)

documentation
"""
function printExperimentInfo(outputFile::String = "", genj::GeneticJuliaStructure = GenJ)

    io = outputFile != "" ? open(outputFile, "a") : Base.stdout

    # Experiment Information
    flush(io)
    println(io, "\n\tExperiment information:")
    println(io, "\t\t· Individual Type → $(genj._experimentInfo._individualType)")
    if isdefined(genj._experimentInfo, :_GPExperimentInfo)
        println(io, "\t\t· Genetic Programming Experiment Information:")

        if typeof(genj._experimentInfo._GPExperimentInfo) <: GEInfo
            #verbose para ge
            print(io, "\t\t\t· Variables → ")
            for i=1:(length(genj._experimentInfo._GPExperimentInfo._variables)-1)
                print(io, "$(genj._experimentInfo._GPExperimentInfo._variables[i]), ")
            end
            (length(genj._experimentInfo._GPExperimentInfo._variables) != 0) ? println(io, "$(genj._experimentInfo._GPExperimentInfo._variables[end])") : println(io, "None")
        elseif typeof(genj._experimentInfo._GPExperimentInfo) <: CGPInfo ||
               typeof(genj._experimentInfo._GPExperimentInfo) <: STGPInfo ||
               typeof(genj._experimentInfo._GPExperimentInfo) <: GEPInfo

            println(io, "\t\t\t· Function Set → $(typeof(genj._experimentInfo._GPExperimentInfo._functionSet)) $(size(genj._experimentInfo._GPExperimentInfo._functionSet))")
            println(io, "\t\t\t· Terminal Set → $(typeof(genj._experimentInfo._GPExperimentInfo._terminalSet)) $(size(genj._experimentInfo._GPExperimentInfo._terminalSet))")

            if typeof(genj._experimentInfo._GPExperimentInfo) <: GEPInfo
                println(io, "\t\t\t· Head Size → $(genj._experimentInfo._GPExperimentInfo._head)")
                println(io, "\t\t\t· Tail Size → $(genj._experimentInfo._GPExperimentInfo._tail)")
            else
                println(io, "\t\t\t· Maximum Tree Depth → $(genj._experimentInfo._GPExperimentInfo._maxTreeDepth)")
            end

            print(io, "\t\t\t· Variables → ")
            for i=1:(length(genj._experimentInfo._GPExperimentInfo._variables)-1)
                print(io, "$(genj._experimentInfo._GPExperimentInfo._variables[i]._name), ")
            end
            (length(genj._experimentInfo._GPExperimentInfo._variables) != 0) ? println(io, "$(genj._experimentInfo._GPExperimentInfo._variables[end]._name)") : println(io, "None")
        end
    end
    println(io, "\t\t· Random Seed → $(genj._experimentInfo._randomSeed)")
    print(io, "\t\t· Evolutionary Algorithm → $(genj._experimentInfo._algorithm)(")
    for i=1:length(genj._experimentInfo._algorithmArgs)-1
        print(io, "$(genj._experimentInfo._algorithmArgs[i]), ")
    end
    length(genj._experimentInfo._algorithmArgs) != 0 ? println(io, "$(genj._experimentInfo._algorithmArgs[end]))") : println(io, ")")

    # Stop Conditions
    println(io, "\tStop Conditions:")
    genj._stopCondition._maxEvaluations != -1 && println(io, "\t\t· Maximum Evaluations → $(genj._stopCondition._maxEvaluations)")
    genj._stopCondition._maxIterations != -1 && println(io, "\t\t· Maximum Iterations → $(genj._stopCondition._maxIterations)")
    genj._stopCondition._maxIterNotImproving != -1 && println(io, "\t\t· Maximum Iterations without improving → $(genj._stopCondition._maxIterNotImproving)")
    genj._stopCondition._maxTime != Inf && println(io, "\t\t· Maximum Time → $(genj._stopCondition._maxTime)")

    # Evaluation
    println(io, "\tEvaluation:")
    for ff in genj._evaluator._fitnessFunctions
        println(io, "\t\t· Fitness function → $(ff._fitnessFunction)")
        println(io, "\t\t· Weight → $(ff._weight)")
    end
    if genj._evaluator._globalFitnessFunction._fitnessFunction != noFunc
        print(io, "\t\t· Global Fitness Function → $(genj._evaluator._globalFitnessFunction._fitnessFunction)(")
        for i=1:length(genj._evaluator._globalFitnessFunction._varArgs)-1
            print(io, "$(genj._evaluator._globalFitnessFunction._varArgs[i]), ")
        end
        length(genj._evaluator._globalFitnessFunction._varArgs) != 0 ? println(io, "$(genj._evaluator._globalFitnessFunction._varArgs[end]))") : println(io, ")")
    end
    println(io, "\t\t· Comparison Function → $(genj._evaluator._compareFunction)")
    println(io, "\t\t· Fitness Comparison Mode → $(genj._evaluator._compareFunctionArgs)")

    # Generation
    println(io, "\tGeneration:")
    println(io, "\t\t· Population Size → $(genj._generator._popSize)")
    println(io, "\t\t· Individuals Generated One By One → $(genj._generator._generateOneByOne)")
    print(io, "\t\t· Generation Method → $(genj._generator._method)(")
    for i=1:length(genj._generator._varArgs)-1
        print(io, "$(genj._generator._varArgs[i]), ")
    end
    length(genj._generator._varArgs) != 0 ? println(io, "$(genj._generator._varArgs[end]))") : println(io, ")")

    # Selection
    println(io, "\tSelection:")
    print(io, "\t\t· Selection Method → $(genj._selector._method)(")
    for i=1:length(genj._selector._varArgs)-1
        print(io, "$(genj._selector._varArgs[i]), ")
    end
    length(genj._selector._varArgs) != 0 ? println(io, "$(genj._selector._varArgs[end]))") : println(io, ")")
    println(io, typeof(genj._selector._nSelectedParents) == Integer ? "\t\t· # of Selected Parents → $(genj._selector._nSelectedParents)" : "\t\t· # of Selected Parents → $(round(genj._selector._nSelectedParents*100, digits=2))%")
    println(io, "\t\t· Method needs to compare between individuals → $(genj._selector._needsComparison)")
    #genj._selector._needsComparison && println(io, "\t\t· Comparison Mode → $(genj._selector._individualMode)")
    println(io, "\t\t· Sampling with Replacement → $(genj._selector._samplingWithRep)")

    # Crossover
    println(io, "\tCrossover")
    print(io, "\t\t· Crossover Method → $(genj._crossoverOp._method)(")
    for i=1:length(genj._crossoverOp._varArgs)-1
        print(io, "$(genj._crossoverOp._varArgs[i]), ")
    end
    length(genj._crossoverOp._varArgs) != 0 ? println(io, "$(genj._crossoverOp._varArgs[end]))") : println(io, ")")
    println(io, "\t\t· Crossover Probability → $(round(genj._crossoverOp._probability*100, digits=2))%")
    println(io, "\t\t· Number of Parents → $(genj._crossoverOp._nParents)")
    println(io, "\t\t· Number of Children → $(genj._crossoverOp._nChildren)")

    # Mutation
    println(io, "\tMutation")
    print(io, "\t\t· Mutation Method → $(genj._mutationOp._method)(")
    for i=1:length(genj._mutationOp._varArgs)-1
        print(io, "$(genj._mutationOp._varArgs[i]), ")
    end
    length(genj._mutationOp._varArgs) != 0 ? println(io, "$(genj._mutationOp._varArgs[end]))") : println(io, ")")
    println(io, "\t\t· Mutation Probability → $(round(genj._mutationOp._probability*100, digits=2))%")

    # Replacement
    println(io, "\tReplacement:")
    print(io, "\t\t· Replacement Method → $(genj._replacementOp._method)(")
    for i=1:length(genj._replacementOp._varArgs)-1
        print(io, "$(genj._replacementOp._varArgs[i]), ")
    end
    length(genj._replacementOp._varArgs) != 0 ? println(io, "$(genj._replacementOp._varArgs[end]))") : println(io, ")")
    println(io, "\t\t· Method needs to compare between individuals → $(genj._replacementOp._needsComparison)")
    #genj._replacementOp._needsComparison && println(io, "\t\t· Comparison Mode → $(genj._replacementOp._individualMode)")

    io == Base.stdout || close(io)
    return nothing
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using DiffWrappers
using Base.Test

using ContinuousTransformations
import DiffResults
import ForwardDiff
using StaticArrays

struct QuadForm{TΣ}
    Σ::TΣ
end

(q::QuadForm)(x) = x' * q.Σ * x

@testset "quadratic form gradient" begin
    for _ in 1:100
        N = rand(3:10)
        static = rand() < 0.5
        A = randn(N, N)
        q = QuadForm(Symmetric(static ? SMatrix{N, N}(A + A') : (A + A')))
        x = static ? SVector{N}(randn(N)) : randn(N)
        qq = ForwardGradientWrapper(q, x)
        qqx = qq(x)
        qq(randn(N))           # another call to test that structure is unshared
        @inferred qq(x)
        @test q(x) == DiffResults.value(qqx)
        @test ForwardDiff.gradient(q, x) == DiffResults.gradient(qqx)
        @test length(qq) == N
    end
end

@testset "TransformLogLikelihood constructor" begin
    ℓ(x) = x[1] + log(x[2])
    t = TransformationTuple((IDENTITY, bridge(ℝ, ℝ⁺)))
    ℓt = TransformLogLikelihood(ℓ, t)
    q = ForwardGradientWrapper(ℓt)
    x = randn(2)
    qx = q(x)
    @test DiffResults.value(qx) == ℓt(x)
    @test DiffResults.gradient(qx) == ForwardDiff.gradient(ℓt, x)
end


text:
@assert ARGS[1] in ["chimera", "square"]
@assert ARGS[2] in ["cpu", "gpu"]

const DEVICE = ARGS[2]

using CUDAnative, CuArrays
device!(parse(Int, get(ARGS, 3, "0")))
using CuYao
using DelimitedFiles, TropicalTensors
using BenchmarkTools

if ARGS[1] == "chimera"
    const LATTICE = ChimeraLattice
else
    const LATTICE = SquareLattice
end

suite = BenchmarkGroup()

function bfunc(::Type{T}, L::Int) where T
    lt = LATTICE(L, L)
    if T <: Integer
        Js = rand(T[-1, 1], length(sgbonds(lt)))
    else
        Js = randn(T, length(sgbonds(lt)))
    end
    hs = zeros(T, length(lt))
    sg = Spinglass(lt, Js, hs)
    if DEVICE == "gpu"
        CuArrays.@sync solve(sg; usecuda=true)
    else
        solve(sg; usecuda=false)
    end
end

if LATTICE === SquareLattice
    const RANGE = 4:2:32
    const TYPE = Float32
else
    const RANGE = 2:8
    const TYPE = Int16
end

for L = RANGE
    suite[L] = @benchmarkable bfunc($TYPE, $L)
    #suite["CPU"][L] = @benchmarkable (CuArrays.@sync spinglass_yao(Float32, $L, $(load_J(L, Val(:randn))); usecuda=false))
end

println("loaded")

tune!(suite)
res = run(suite)

function analyze_res(res)
    times = zeros(length(res))
    for (k,L) = enumerate(RANGE)
        times[k] = minimum(res[L].times)
    end
    return times
end

times = analyze_res(res)
fname = joinpath(@__DIR__, "data", "bench_$(LATTICE)_$DEVICE.dat")
println("Writing benchmark results to file: $fname.")
writedlm(fname, times)


text:
#=Copyright (c) 2017 Gabriel Goh

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
the Software, and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.=#

using Plots

# A = ones(18,67)

# A[2:end-1,2:end-1] = 0
# A[1:7,32-5:34+5] = 1
# A[end-7:end,32-5:34+5] = 1

Box = zeros(16,16)
Cavity = ones(16,16)
Cavity[8:9,:] = 0
R = [Box Cavity Box Cavity Box Cavity Box]
A = ones(size(R,1)+2, size(R,2) + 2)
A[2:end-1,2:end-1] = R

n = size(A,1) - 2
m = size(A,2) - 2

R = reshape(1:n*m, n,m)
Id(i,j) = R[i-1,j-1]

Q = zeros(n*m,n*m)
for i = 2:(size(A,1)-1), j = 2:(size(A,2)-1)
  k = Id(i,j)
  if A[i,j] == 0
    A[i,j+1] == 0 ? Q[Id(i,j+1),k] = 1 : nothing
    A[i+1,j] == 0 ? Q[Id(i+1,j),k] = 1 : nothing
    A[i,j-1] == 0 ? Q[Id(i,j-1),k] = 1 : nothing
    A[i-1,j] == 0 ? Q[Id(i-1,j),k] = 1 : nothing
  end
end

Q = (Q + Q')/2

v = sum(Q,1)[:]
v[v .== 0] = 1
L = diagm(v) - Q
L[1,1] = 2.1

Σ,U = eig(L)

pyplot(leg=false, ticks=nothing)
@gif for i in 2:500
  plot(size = (880,140)); heatmap!(reshape(U[:,224*3+i],n,m))
end

e1 = zeros(n*m)
e1[8:9] = 1

α = 1/max(eig(L)[1]...)
iters = 5000
X = zeros(n*m,iters)
x = 10*zeros(n*m)
for i = 1:iters
  X[:,i] = x
  for i = 1:100
    x = x - (1*α)*(L*x - e1)
  end
  println(norm(L*x - e1))
end

plot(size = (880,140)); heatmap!(reshape(X[:,100],n,m))
pyplot(leg=false, ticks=nothing)

@gif for i in 10:10:2000
  println(i)
  xi = reshape(X[:,i],n,m)
  plot(size = (880,140)); heatmap!(xi)
end

lambdak(k,λ) = α*(1 - (1 - α*λ).^k)./(α*λ)
c = U'e1

@gif for i in 10:500:300000
  println(i)
  z = U*(lambdak(i,Σ).*(c)
  plot(size = (880,140)); heatmap!(reshape(z, n, m))
end

function getiterk(Q,b)
  Σ,U = eig(Q)
  c = U'b
  lambdak(k,λ) = α*(1 - (1 - α*λ).^k)./(α*λ)
  return U*(lambdak(k,Σ).*c)
end

max_stars_repo_name: sbeleidy/post--momentum

max_issues_repo_name: sbeleidy/post--momentum

max_forks_repo_name: sbeleidy/post--momentum