text:
function search_template(model::TierraModel, organism::TierrianOrganism, start_address, direction::SearchDirection)
    template = _get_template(model, start_address)
    template_length = length(template)

    if 0 < template_length
        for distance in 1:MAX_SEARCH_DISTANCE
            template_address = _get_address_if_match(direction, model, start_address, template, distance)

            if (template_address != -1)
                return template_address, template_length
            end
        end
    end

    # no match found -> set error flag and return address after template
    organism.error_flag = true

    template_address = normalize_address(start_address + template_length)
    return template_address, template_length
end

function _get_address_if_match(direction::Union{SearchBackward, SearchForward}, model, template_address, template, distance)
    template_length = length(template)
    start_address = _get_address_within_distance(direction, template_address, template_length, distance)

    for i in 0:(template_length - 1)
        instruction = read_memory(model, start_address + i)

        if instruction != template[i+1]
            return -1
        end
    end

    next_address = normalize_address(start_address + template_length)
    return next_address
end

function _get_address_if_match(::SearchBoth, model, template_address, template, distance)
    forward_address = _get_address_if_match(SearchForward(), model, template_address, template, distance)

    if forward_address != -1
        return forward_address
    else
        return _get_address_if_match(SearchBackward(), model, template_address, template, distance)
    end
end

function _get_address_within_distance(::SearchBackward, template_address, template_length, distance)
    return normalize_address(template_address - distance - template_length + one(UInt16))
end
function _get_address_within_distance(::SearchForward, template_address, template_length, distance)
    return normalize_address(template_address + template_length + distance - one(UInt16))
end

function _get_template(model, address)
    template = []

    current_address = address

    while length(template) < MAX_TEMPLATE_LENGTH
        instruction_int = read_memory(model, current_address)
        instruction = convert_instruction(instruction_int)

        if (instruction != NOP_0() && instruction != NOP_1()) break end

        push!(template, 1-instruction_int)
        current_address += one(UInt16)
    end

    return template
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
struct LiquidProperties{P, U}
    η::P # volume fraction
    T::P # reduced temperature
    potential::U
end

LiquidProperties(η::P, potential::U) where {P, U<:Potential} =
    LiquidProperties{P, U}(η, one(P), potential)

potential(p::LiquidProperties) = p.potential
temperature(p::LiquidProperties) = p.T
volume_fraction(p::LiquidProperties) = p.η

### Liquids types
abstract type AbstractLiquid{U} end

const LIQUIDS_NAMES = (
    :AttractiveHardSpheres,
    :DipolarHardSpheres,
    :HardDisks,
    :HardSpheres
)

for liquid in LIQUIDS_NAMES
    @eval begin
        struct $liquid{P, U} <: AbstractLiquid{U}
            properties::LiquidProperties{P, U}
        end
    end
end

### Liquids constructors
function HardDisks(η::Real)
    properties = LiquidProperties(η, HardCore{2}())
    return HardDisks(properties)
end

function HardSpheres(η::Real)
    properties = LiquidProperties(η, HardCore{3}())
    return HardSpheres(properties)
end

function AttractiveHardSpheres(η::Real, T::Real, potential::CompositePotential)
    η′, T′ = promote(η, T)
    properties = LiquidProperties(η′, T′, potential)
    return AttractiveHardSpheres(properties)
end
#
function AttractiveHardSpheres(η::Real, T::Real, potential::AttractivePotential)
    T′ = T / energy_scale(potential)
    potential′ = CompositePotential(HardCore{3}(), potential)
    return AttractiveHardSpheres(η, T′, potential′)
end

function DipolarHardSpheres(η::Real, T::Real)
    η′, T′ = promote(η, T)
    potential = CompositePotential(HardCore{3}(), DipoleDipole())
    properties = LiquidProperties(η′, T′, potential)
    return DipolarHardSpheres(properties)
end

potential(l::AbstractLiquid) = potential(l.properties)
temperature(l::AbstractLiquid) = temperature(l.properties)
volume_fraction(l::AbstractLiquid) = volume_fraction(l.properties)

ext: jl


text:
using SafeTestsets

@time begin
    using eFEMpart
#    @time @safetestset "Cell List Tests" begin include("cl_test.jl") end
end

max_stars_repo_name: pseastham/eFEMpart


text:
@propagate_inbounds find_all(f, x) = find_all(f, x, order(x))
find_all(f, x, xo) = findall(f, x)

find_all(f::Fix2{typeof(in)}, y, yo) = _findin(f.x, order(f.x), y, yo)

_fallback_find_all(f, a) = collect(first(p) for p in pairs(a) if f(last(p)))

# <, <=
@propagate_inbounds find_all(f::F2Lt, r, ro::ForwardOrdering) = _find_all(keytype(r), firstindex(r), find_last(f, r, ro))
@propagate_inbounds find_all(f::F2Lt, r, ro::ReverseOrdering) = _find_all(keytype(r), find_first(f, r, ro), lastindex(r))
find_all(f::F2Lt, r::AbstractRange, ::UnorderedOrdering) = _empty_ur(keytype(r))
find_all(f::F2Lt, r::AbstractArray, ::UnorderedOrdering) = _fallback_find_all(f, a)

# >, >=
@propagate_inbounds find_all(f::F2Gt, r, ro::ForwardOrdering) = _find_all(keytype(r), find_first(f, r, ro), lastindex(r))
@propagate_inbounds find_all(f::F2Gt, r, ro::ReverseOrdering) = _find_all(keytype(r), firstindex(r), find_last(f, r, ro))
find_all(f::F2Gt, r::AbstractRange, ::UnorderedOrdering) = _empty_ur(keytype(r))
find_all(f::F2Gt, r::AbstractArray, ::UnorderedOrdering) = _fallback_find_all(f, a)

# find_all(==(x), r)
find_all(f::F2Eq, r::AbstractRange, ::UnorderedOrdering) = _empty_ur(keytype(r))
find_all(f::F2Eq, r, ro::Ordering) = _find_all(keytype(r), find_first(f, r, ro), find_last(f, r, ro))

# only really applies to ordered vectors
_find_all(::Type{T},        fi,        li) where {T} = fi:li
_find_all(::Type{T}, ::Nothing,        li) where {T} = _empty_ur(T)
_find_all(::Type{T},        fi, ::Nothing) where {T} = _empty_ur(T)
_find_all(::Type{T}, ::Nothing, ::Nothing) where {T} = _empty_ur(T)

_empty_ur(::Type{T}) where {T} = one(T):zero(T)

@propagate_inbounds find_all(f::BitAnd, r, ro) = _bit_and(r, ro, find_all(f.f1, r, ro), find_all(f.f2, r, ro))
@propagate_inbounds find_all(f::BitOr, r, ro) = _bit_or(r, ro, find_all(f.f1, r, ro), find_all(f.f2, r, ro))

function _bit_and(r, ro, inds1, inds2)
    if isempty(inds1)
        return isempty(inds2) ? inds2 : inds1
    else
        return isempty(inds2) ? inds2 : intersect(inds1, inds2)
    end
end

function _bit_or(r, ro, inds1, inds2)
    if isempty(inds1)
        return isempty(inds2) ? inds1 : inds2
    else
        return isempty(inds2) ? inds1 : _merge_bit_find(inds1, inds2)
    end
end

# we assume that neither of these are empty at this point
function _merge_bit_find(inds1::AbstractRange{T}, inds2::AbstractRange{T}) where {T}
    if is_after(inds1, inds2)
        return vcat(inds2, inds1)
    elseif is_before(inds1, inds2)
        return vcat(inds1, inds2)
    else
        return _group_min(inds1, Forward, inds2, Forward):_group_max(inds1, Forward, inds2, Forward)
    end
end

@propagate_inbounds find_all(f::Fix2{typeof(!=)}, r, ro) = find_all(<(f.x) | >(f.x), r, ro)


text:
# This is the script used to generate the results comparing mix-SQP
# vs. the EM algorithm. The results for the projected gradient method
# are generated separately with run_pg.m.
csvfile         = "simdata-n=20000-m=20.csv"
outfile_em      = "em-n=20000-m=20.csv";
outfile_mixsqp1 = "mixsqp-exact-n=20000-m=20.csv";
outfile_mixsqp2 = "mixsqp-approx-n=20000-m=20.csv";

using Printf
using DelimitedFiles
using LinearAlgebra
using SparseArrays
using LowRankApprox
using Statistics
include("../code/mixEM.jl");
include("../code/mixSQP.jl");

# Load data..
@printf("Loading data.\n")
L = readdlm(csvfile,',');
n = size(L,1);
m = size(L,2);
@printf "Loaded %d x %d matrix.\n" n m

# Run the methods for a small number of iterations to precompile the
# code.
@printf "Precompiling EM and mixSQP code.\n"
xem, fem, tem = mixEM(L,maxiter = 10);
outsqp1 = mixSQP(L,lowrank = "none",maxiter = 4,maxqpiter = 4,verbose = false);
outsqp2 = mixSQP(L,lowrank = "qr",maxiter = 4,maxqpiter = 4,verbose = false);

# Run the EM algorithm.
@printf "Fitting model using EM.\n"
@time xem, fem, dem, tem = mixEM(L,maxiter = 1000,tol = 1e-6);
fem = fem/n;
writedlm(outfile_em,[fem tem],',');

# Run mix-SQP with no approximation to the input matrix.
@printf "Fitting model using mix-SQP with exact L.\n"
@time outsqp1 = mixSQP(L,lowrank = "none",eps = 1e-8,sptol = 1e-6,
                       maxqpiter = m,maxiter = 200,verbose = false);
outsqp1["obj"] = outsqp1["obj"]/n;
writedlm(outfile_mixsqp1,[outsqp1["obj"] outsqp1["timing"]],',');

# Run the mix-SQP with a low rank (truncated QR) approximation to the
# input matrix.
@printf "Fitting model using mix-SQP with approximate L.\n"
@time outsqp2 = mixSQP(L,lowrank = "qr",eps = 1e-8,sptol = 1e-6,
                       maxqpiter = m,maxiter = 200,verbose = false);
outsqp2["obj"] = outsqp2["obj"]/n;
writedlm(outfile_mixsqp2,[outsqp2["obj"] outsqp2["timing"]],',');

# Compare the quality of the solutions. Note that we need to divide
# the objective by n to match the objective used in the manuscript.
@printf "Objective at EM sol.:   %0.12f\n" mixobjective(L,xem,1e-8)/n
@printf "Objective at SQP1 sol.: %0.12f\n" mixobjective(L,outsqp1["x"],1e-8)/n
@printf "Objective at SQP2 sol.: %0.12f\n" mixobjective(L,outsqp2["x"],1e-8)/n

max_stars_repo_name: stephenslab/mixsqp-paper

max_issues_repo_name: stephenslab/mixsqp-paper

max_forks_repo_name: stephenslab/mixsqp-paper


text:
test_f(x, y) = x^2 + 3y + 5

function test_f_2()

    @syms a b
    expr = simplify(a+(b+1)+4+12^2)
    println(expr)
    println(a+1)

end


text:
# save to disk/read from disk
savfile1  = "test.seis"
savfile2  = "test.hdr"
savfile3  = "test.evt"
legf_050  = path * "/SampleFiles/SEIS/2011.314.7D.J34A..HH1.seis"
legf_052  = path * "/SampleFiles/SEIS/seisio_testfile_v0.52.seis"
uw_file   = path * "/SampleFiles/UW/02062915175o"

# Changing this test to guarantee at least one campaign-style measurement ... and test splat notation ... and something with no notes
printstyled("  SeisData\n", color=:light_green)
S = breaking_seis()
wseis(savfile1, S)
R = rseis(savfile1)[1]
@test(R==S)
S.notes[2][1] = string(String(Char.(0x00:0x7e)), String(Char.(0x80:0xff)))
wseis(savfile1, S)
R = rseis(savfile1)[1]
@test(R!==S)

printstyled("  SeisHdr\n", color=:light_green)
H = randSeisHdr()
wseis(savfile2, H)
H2 = rseis(savfile2)[1]
@test(H==H2)

printstyled("  SeisEvent\n", color=:light_green)
EV = SeisEvent(hdr=H, data=convert(EventTraceData, S))
EV.data.misc[1] = breaking_dict
wseis(savfile3, EV)
R = rseis(savfile3)[1]
@test R == EV

printstyled("  read/write of each type to same file\n", color=:light_green)
Ch = randSeisChannel(s=true)
Ch.x = rand(24)
Ch.t = vcat(Ch.t[1], [24 0])
wseis(savfile3, EV, S, H, Ch)

R = rseis(savfile3)
@test(R[1]==EV)
@test(R[2]==S)
@test(R[3]==H)
@test(R[4]==Ch)
@test(S.misc[1] == R[1].data.misc[1] ==  R[2].misc[1])

# read one file with one record number
printstyled("  read file with integer record number\n", color=:light_green)
R = rseis("test.seis", c=1, v=1)
@test R[1] == S

# read everything
printstyled("  read a multi-record file\n", color=:light_green)
R = rseis("test*", v=1)
@test R[3] == R[5] # Header is read twice, test.evt (file 1) record 3, test.hdr (file 2) record 1
@test R[2] == R[6] # Seis is read twice, test.evt (file 1) record 2, test.seis (file 3) record 1

# read when some files have record 3 but not all
printstyled("  read file list with list of record numbers\n", color=:light_green)
R = rseis("test.*", c = [1,3], v=1)
@test(R[3]==R[2])
@test(R[1].data.misc[1]==R[4].misc[1])

# read nothing as each target file has one record
printstyled("  read nothing due to an intentionally poor choice of record numbers\n", color=:light_green)
R = rseis(["test.seis", "test.h*"], c=[2, 3], v=1)
@test isempty(R)

# read the first record of each file
printstyled("  read first record from each SeisIO file using a wildcard list\n", color=:light_green)
R = rseis("test*", c=1, v=1)
@test R[1] == EV
@test R[2] == H
@test R[3] == S

printstyled("  test that every custom Type can be written and read faithfully\n", color=:light_green)
redirect_stdout(out) do
  A = Array{Any,1}(undef, 0)
  for T in SeisIO.TNames
    println("testing ", T)
    if T == PhaseCat
      push!(A, randPhaseCat())
    elseif T == MultiStageResp
      push!(A, MultiStageResp(6))
    else
      push!(A, getfield(SeisIO, Symbol(T))())
    end
  end
  wseis(savfile1, A...)
  R = rseis(savfile1, v=2)

  for i = 1:length(R)
    @test R[i] == A[i]
  end

  # add an incompatible type; should throw a warning, and another when read back in
  push!(A, rand(Float64,3))
  wseis(savfile1, A...)
  R = rseis(savfile1, v=2)
end

printstyled("  test read/write with data compression\n", color=:light_green)
SeisIO.KW.comp = 0x02
S = randSeisData()
nx = 4
S.t[1] = [1 0; nx 0]
S.x[1] = randn(eltype(S.x[1]), nx)
wseis(savfile1, S)
R = rseis(savfile1, v=2)[1]
@test R == S

SeisIO.KW.comp = 0x01
S = randSeisEvent()
C = SeisChannel()
nx = SeisIO.KW.n_zip*2
C.t = [1 0; nx 0]
C.x = randn(nx)
n = S.data.n
push!(S.data, C)
@test S.data.n == n+1
C = SeisChannel()
nx = 4
C.t = [1 0; nx 0]
C.x = randn(nx)
push!(S.data, C)
wseis(savfile1, S)
R = rseis(savfile1, v=2)[1]
@test R == S

# read_data("seisio", ...)
S1 = read_data(savfile1)
@test convert(SeisData, S.data) == S1
S1 = read_data("seisio", savfile1)
@test convert(SeisData, S.data) == S1
S2 = rseis(savfile3)
S1 = read_data(savfile3)
Tr = S2[1].data
@test Tr.n == S1.n
for i = 1:Tr.n
  if any(isnan, Tr.x[i]) == false
    @test Tr.x[i] == S1.x[i]
  end
end

# Type unit tests with read_data("seisio", ...)
C = randSeisChannel()
wseis(savfile1, C)
S1 = read_data("seisio", savfile1)
@test S1[1] == C
C = convert(EventChannel, randSeisChannel())
wseis(savfile1, C)
S1 = read_data("seisio", savfile1)
@test S1[1] == convert(SeisChannel, C)
S = randSeisData()
wseis(savfile1, S)
S1 = read_data("seisio", savfile1)
@test S1 == S
S = convert(EventTraceData, S)
wseis(savfile1, S)
S1 = read_data("seisio", savfile1)
@test S1 == convert(SeisData, S)
Ev = randSeisEvent()
L = GeoLoc(lat=45.560504, lon=-122.617020, el=51.816, az=180.0, inc=0.0)
wseis(savfile1, L, Ev)
S1 = read_data("seisio", savfile1)
@test S1 == convert(SeisData, Ev.data)

rm(savfile1)
rm(savfile2)
rm(savfile3)

# Legacy file reading ========================================================
# 0.50  all custom types can use write(); rseis, wseis no longer required
#       String arrays and :misc are written in a completely different way
#       Type codes for :misc changed
#       deprecated BigFloat/BigInt support in :misc
#       :n is no longer stored as a UInt32
#       :x compression no longer automatic and changed from Blosc to lz4

set_file_ver(legf_050, 0.4)
@test get_file_ver(legf_050) == 0.4f0
@test_throws ErrorException rseis(legf_050)
@test_throws ErrorException get_file_ver(uw_file)
set_file_ver(legf_050, 0.50)
S = rseis(legf_050)[1]
@test S.n == 3
@test S.id == ["7D.J34A..HH1", "7D.J34A..HH2", "7D.J34A..HHZ"]

# 0.52  2019-09-03    added Types: CoeffResp, MultiStageResp
# 0.51  2019-08-01    added :f0 to PZResp, PZResp64
# (legacy reader is the same for both)

S = rseis(legf_052)[1]
@test S.n == 3
@test S.id == ["IU.COR.00.BH1", "IU.COR.00.BH2", "IU.COR.00.BHZ"]
@test typeof(S.resp[1]) == MultiStageResp
@test typeof(S.resp[1].stage[1]) == PZResp64
@test isapprox(S.resp[1].stage[1].p, [ -0.0173949 + 0.01234im,
                                       -0.0173949 - 0.01234im,
                                       -0.0175489 + 0.0im    ,
                                       -0.0591783 + 0.0im    ,
                                           -39.18 + 49.12im  ,
                                           -39.18 - 49.12im ])
@test isapprox(S.resp[1].gain, [2660.0, 1.67772e6, 1.0])
@test S.resp[1].fac == [0,1,1]
@test S.resp[1].o[1] == "m/s"
@test S.resp[1].i[1] == "v" # note: this was wrong at the time in absolute sense. Should be "V" in UCUM.


text:
export Fun, VFun, RFun, NFun, NVFun, NRFun, EFun, NEFun,
             @Fun, @VFun, @RFun, @NFun, @NVFun, @NRFun, @EFun, @NEFun, FunDef, @CFun

struct Fun{R, A, H}
    _function::Function

    "Pass Through"
    Fun(_function::Fun) = _function
    Fun{R, A, H}(_function::Fun) where {R, A, H} = Fun{R, A, H} == typeof(_function) ? _function : new{R, A, H}(_function._function)
    Fun(R, A, H, _function::Function) = new{R, A, H}(_function)
    Fun{R, A, H}(_function::Function) where {R, A, H} = Fun(R, A, H, _function)

    function Fun(_function::Function)
        r = Base.return_types(_function)[1]
        a = methods(_function)[1].sig.parameters[2:end]
        if length(a) > 1
            a = Tuple{a...,}
        end
        Fun(r, a, _function, Symbol(_function))
    end

    Base.string(f::Fun) = string(typeof(f))
    Base.isequal(f::Fun, f2::Fun) = (typeof(f) == typeof(f2)) && f._function == f2._function
    Base.show(io::IO, f::Fun) = print(io, string(f))
    Base.print(io::IO, f::Fun) = print(io, string(f))
    Docs.getdoc(f::Fun) = string(f)

    checkeq(::Type{Fun{R, A, H}}, ::Type{Fun{R2, A2, H2}}) where {R, A, H, R2, A2, H2} = (R == R2 && A == A2 && H == H2)
    Base.isequal(t::Type{F1}, t2::Type{F2}) where {F1 <: Fun, F2 <: Fun} = checkeq(t, t2)
    Base.print(io::IO, f::Type{F}) where F <: Fun = print(io, string(f))
    Base.show(io::IO, f::Type{F}) where F <: Fun = show(io, string(f))
    Base.string(f::Type{Fun{R, A, H}}) where {R, A, H} = "Fun{$R, $A}"
    Docs.getdoc(f::Type{F}) where F <: Fun = string(f)

end

function fun_str(alias, msg, generic_values)
    idx = 1
    while alias isa UnionAll
        msg = replace(msg, (string(alias.var.name) * "~") => generic_values[idx])
        alias = alias.body
        idx += 1
    end
    replace(msg, "~" => "")
end

fun_str(msg) = replace(msg, "~" => "")

"Default Function TypeDef Generator"
function FunDef(mod, line, fun, ret, args, msg)
    partialTypeDef = nothing
    partialStringDef = nothing
    partialCallDef = nothing
    partialStringDef2 = nothing

    try
        #Strip Return Label
        ret = (ret isa Expr && ret.head == :(::)) ? ret.args[2] : ret
        argtypes = nothing
        argnames = nothing
        args = Any[arg for arg in args]
        aliasName = fun isa Expr ? fun.args[1] : fun
        generics = fun isa Expr ? fun.args[2:end] : []

        if args isa Expr || args isa Symbol
            #Strip Arg Labels
            if args isa Expr && args.head == :(::)
                argtypes = [args.args[2]]
                argnames = [args.args[1]]
            else
                #Add Argument Label
                argnames = [:arg1]
                argtypes = [args]
                args = Expr(:(::), :arg1, args)
            end
        else
            argtypes = Any[arg for arg in args]
            argnames = [Symbol("arg$i") for i in 1:length(args)]
            for i in 1:length(args)
                if argtypes[i] isa Expr && args[i].head == :(::)
                    argnames[i] = argtypes[i].args[1]
                    argtypes[i] = argtypes[i].args[2]
                else
                    args[i] = Expr(:(::), argnames[i], args[i])
                end
            end
            argtypes = Expr(:curly, :Tuple, argtypes...)
        end

        #Create Alias (Use Hash/Rand for Equality Checking)
        curlyBlock = Expr(:curly, GlobalRef(@__MODULE__, :Fun), ret, argtypes, msg === nothing ? rand() : hash(msg))
        partialTypeDef = Expr(:const, Expr(:(=), fun, curlyBlock))

        aliasResult = Core.eval(mod, partialTypeDef)

        #Define most of the specialized function call expr
        partialCallDef = :(@inline function (f::$fun)()::$ret where {}
                                        f._function()
                                   end)

        partialCallDef.args[3].args[2].args[1] = LineNumberNode(line.line, line.file)
        partialCallDef.args[3].args[2].args[2] = LineNumberNode(line.line, line.file)

        #Append Function Generics
        whereExpr = partialCallDef.args[3].args[1].args
        append!(whereExpr, generics)

        #Append the arguments to the function Call
        funArgs = partialCallDef.args[3].args[1].args[1].args[1].args
        append!(funArgs, args)

        #Append the arguments to the inner function
        funBlock = partialCallDef.args[3].args[2].args[3].args
        append!(funBlock, argnames)

        Core.eval(mod, partialCallDef)
        msg === nothing && return aliasResult

        #Eval No Generic Message
        if length(generics) > 0
            #Define most of the function string expr
            partialStringDef2 = :(Base.string(x::Type{$aliasName}) = fun_str($msg))

            #Replace with safer reference
            fun_str_call = partialStringDef2.args[2].args[2].args
            fun_str_call[1] = GlobalRef(@__MODULE__, :fun_str)

            Core.eval(mod, partialStringDef2)
        end

        #Define most of the function string expr
        partialStringDef = :(Base.string(x::Type{$fun}) where {} = fun_str($aliasName, $msg))

        #Append Function Generics
        whereExpr = partialStringDef.args[1]
        append!(whereExpr.args, generics)

        fun_call = partialStringDef.args[2].args[2].args
        #Replace with safer reference
        fun_call[1] = GlobalRef(@__MODULE__, :fun_str)

        #Create Tuple instead of curly to pass to function
        tupleExpr = Expr(:tuple)
        append!(tupleExpr.args, generics)
        push!(fun_call, tupleExpr)

        #Set Line Information
        partialStringDef.args[2].args[1] = LineNumberNode(line.line, line.file)
        Core.eval(mod, partialStringDef)
        return aliasResult
    catch e
        println("Error in Function Generating $fun in $mod")
        println("PartialTypeDef:$partialTypeDef")
        println("PartialCallDef:$partialCallDef")
        println("PartialStringDef2:$partialStringDef")
        println("PartialStringDef:$partialStringDef")
        throw(e)
    end
end

function deffhash(mod, expr, ret, args)
    str = "$mod.$expr("
    startswith(str, "Main.") && (str = SubString(str, 6, length(str)))
    if args != :Nothing
        str *= join(args, "~, ") * "~"
    end

    (ret == :Nothing) && return str * ")"
    return string(str * ")::($ret~)")
end

"Default Message Function (Args & Return) Type Definition Generator"
macro Fun(expr, ret, args...)
    FunDef(__module__, __source__,  expr, ret, args, deffhash(__module__,expr, ret, args))
end

"Custom FuncType And Additional Function (Args & Return) Type Definition Generator"
macro CFun(expr, ret, msg, args...)
    FunDef(__module__, __source__,  expr, ret, args, deffhash(msg, expr, ret, args))
end

"Void Message Function (Args Only)"
const VFun{A} = Fun{Nothing, A}

"Default Void Message Function (Args Only) Type Definition Generator"
macro VFun(expr, args...)
    return FunDef(__module__, __source__, expr, :Nothing, args, deffhash(__module__, expr, :Nothing, args))
end

"Void Message Function"
const EFun{MSG} = Fun{Nothing, Nothing}

"Void Default Message Function Type Definition Generator"
macro EFun(expr)
    return FunDef(__module__, __source__, expr, :Nothing, :Nothing, deffhash(__module__,expr, :Nothing, args))
end

"Default Return Message Function (Return Only)"
const RFun{R} = Fun{R, Nothing}

"Default Return Message Function (Return Only) Type Definition Generator"
macro RFun(expr, ret)
    return FunDef(__module__, __source__, expr, ret, :Nothing, deffhash(__module__,expr, ret, :Nothing))
end

"No Message Function (Args & Return)"
const NFun{R, A} = Fun{R, A, Nothing}

"No Message Function (Args & Return) Type Definition Generator"
macro NFun(expr, ret, args...)
    return FunDef(__module__, __source__, expr, ret, args, nothing)
end

"No Message Void Function (Args only)"
const NVFun{A} = Fun{Nothing, A, Nothing}

"No Message Void Function (Args only) Type Definition Generator"
macro NFun(expr, args...)
    FunDef(__module__, __source__, expr, :Nothing, args, nothing)
end

"No Message Return Function (Return only)"
const NRFun{R} = Fun{R, Nothing, Nothing}

"No Message Return Function (Return only) Type Definition Generator"

macro NFun(expr, ret,  args...)
    FunDef(__module__, __source__, expr, ret, :Nothing, nothing)
end

"Void No Message Function"
const NEFun = Fun{Nothing, Nothing, Nothing}

"Void No Message Function Type Definition Generator"
macro NEFun(expr)
    FunDef(__module__, __source__, expr, :Nothing, :Nothing, nothing)
end


text:
@testset "150.evaluate-reverse-polish-notation.jl" begin
    @test eval_rpn(["4", "13", "5", "/", "+"]) == 6
    @test eval_rpn(["2", "1", "+", "3", "*"]) == 9
    @test eval_rpn(["10", "6", "9", "3", "+", "-11", "*", "/", "*", "17", "+", "5", "+"]) == 22
    @test eval_rpn(["2", "1", "-", "3", "*"]) == 3
end


text:
using Revise
using Bolt
using ForwardDiff
#using PyPlot
using Plots
pyplot()
using BenchmarkTools
using Printf
#using ThreadPools

#input ingredients
par = CosmoParams()
bg = Background(par)
ih = IonizationHistory(Peebles(), par, bg)
k_grid = quadratic_k(0.1bg.H₀, 1000bg.H₀, 1000) #quadratically spaced k points
println(bg.H₀) #the units on this are in H0?
 #must need to account for evoln somehow?
#numerical parameters
ℓᵧ=8 #cutoff
ℓ_ν=10 #not used except for size here, should pass
reltol=1e-11

#solve hierarchy at single x to check
x=-8 #just picking a number
pertlen = 2(ℓᵧ+1)+(ℓ_ν+1)+5
println("pert vector length=",pertlen)
results=zeros(pertlen,length(k_grid))

#@btime @qthreads
#seems to take like a minute (without threads)
for (i_k, k) in enumerate(k_grid)
    hierarchy = Hierarchy(BasicNewtonian(), par, bg, ih, k, ℓᵧ)
    perturb = boltsolve(hierarchy; reltol=reltol)
    u = perturb(x)  #z this can be optimized away, save timesteps at the grid!
    results[:,i_k] = u #z should use unpack somehow
end
#'Need to figure out the units of k here'
labels = [raw"$\Phi$",raw"$\delta$",raw"$v$",raw"$|\delta_{b}|$",raw"$|v_{b}|$"]
plot()
for i in 1:5
    plot!(log10.(k_grid/ bg.H₀), log10.(abs.(results[2(ℓᵧ+1)+(ℓ_ν+1)+i,:])),label=labels[i])
end
plot!(log10.(k_grid/ bg.H₀), log10.(abs.(results[1,:])),label=raw"$|\Theta_{0}|$")

plot!(log10.(k_grid/ bg.H₀), log10.(abs.(results[2(ℓᵧ+1)+1,:])),
      label=raw"$|\mathcal{N}_{0}|$")
# plot ℓ=1 perts:
# plot!(log10.(k_grid/ bg.H₀), log10.(abs.(results[2,:])),label=raw"$|\Theta_{1}|$")
# plot!(log10.(k_grid/ bg.H₀), log10.(abs.(results[2(ℓᵧ+1)+2,:])),
#     label=raw"$|\mathcal{N}_{1}|$")
#xlims!(1,3)

ylabel!(raw"$\delta_{i}(k)$")
xlabel!(raw"$k \ H₀$")
title!("z=$(@sprintf("%.0f", exp(-x)-1))")
savefig("../compare/bolt_perts_k_z$(@sprintf("%.0f", exp(-x)-1)).png")


text:
using Printf

function readData(filename::AbstractString)

    !isfile(filename) && error("Input file not found")

	phylo_profile = open(filename)
	list_species = String[]
	list_domains = String[]
	for line in eachline(phylo_profile)
		data = (split(chomp(line)))
		N = length(data)
		for k = 1:N
			if(k == 1)
				if(in(data[k],list_species) == false)
					push!(list_species,data[k])
				end

			else
				if(in(data[k],list_domains) == false)
					push!(list_domains,data[k])
				end
			end
		end
	end
	close(phylo_profile)
	sort!(list_species)
	sort!(list_domains)

	return list_species,list_domains
end

function makePhyloMatrix(filename::AbstractString,
                         list_species_ord,
                         list_domains_ord)

	#N=number of domains, M=number of species
	M=length(list_species_ord)
	N=length(list_domains_ord)

	#make phylogenetic profile
	#P=phylogenetic matrix, M times N matrix
	P=zeros(Int8,M,N)

    #map species/domains to their index in the list
    dict_species=Dict(list_species_ord[i] => i for i=1:length(list_species_ord))
    dict_domains=Dict(list_domains_ord[i] => i for i=1:length(list_domains_ord))

	phylo_profile=open(filename)
	for line in eachline(phylo_profile)
		data=(split(chomp(line)))
		len=length(data)
		index_spec=0
		index_dom=0
        index_spec=dict_species[data[1]]
		for k=2:len
			index_dom=dict_domains[data[k]]
			P[index_spec,index_dom]=1 #binary matrix
		end
	end
	return P
end

#############################################################
#TO be implemented well, that's just for now
function read_benchmark(filename::AbstractString)

	!isfile(filename) && error("Benchamark file not found")

	bench_matr=readdlm(filename)
	return bench_matr
end
##################################################

function print_result(final_matrix::Matrix)
	result_file=open("results.dat", "w")

	lungh,num=size(final_matrix)
    for i=1:lungh
        @printf(result_file,"%s %s %f ", final_matrix[i,1], final_matrix[i,2], final_matrix[i,3])
        for k=4:num
            @printf(result_file, " %d ", final_matrix[i,k])
        end
        @printf(result_file, "\n")
    end
    close(result_file)
end

max_stars_repo_name: clayfos/PhyDCA

max_issues_repo_name: clayfos/PhyDCA

max_forks_repo_name: clayfos/PhyDCA


text:
immutable BandedQR{T} <: Factorization{T}
    H::Matrix{T}  # Represents the orthogonal matrix Q
    R::BandedMatrix{T}
end

Base.size(QR::BandedQR,k) = size(QR.R,k)

immutable BandedQ{T} <: AbstractMatrix{T}
    H::Matrix{T}
    m::Int   # number of rows/cols
end

size(A::BandedQ) = (A.m, A.m)
size(A::BandedQ,i::Integer) = i <= 0 ? error("dimension out of range") :
                                i == 1 ? A.m :
                                i == 2 ? A.m : 1

Base.At_mul_B{T<:Real}(A::BandedQ{T},B::Union{Vector{T},Matrix{T}}) = Ac_mul_B(A,B)
Base.At_mul_B!{T<:Real}(Y,A::BandedQ{T},B::Union{Vector{T},Matrix{T}}) = Ac_mul_B!(Y,A,B)

Base.Ac_mul_B{T<:Real,V<:Complex}(A::BandedQ{T},B::Vector{V}) =
    Ac_mul_B(A,real(B))+im*Ac_mul_B(A,imag(B))
Base.Ac_mul_B{T<:Real,V<:Real}(A::BandedQ{T},B::Vector{V}) =
    Ac_mul_B(A,Vector{T}(B))
Base.Ac_mul_B{T<:Complex}(A::BandedQ{T},B::Vector) =
    Ac_mul_B(A,Vector{T}(B))
Base.Ac_mul_B{T<:Real}(A::BandedQ{T},B::Vector{T}) = Ac_mul_B!(similar(B),A,B)
Base.Ac_mul_B{T<:Complex}(A::BandedQ{T},B::Vector{T}) = Ac_mul_B!(similar(B),A,B)

function Base.Ac_mul_B!{T<:BlasFloat}(Y::Vector{T},A::BandedQ{T},B::Vector{T})
    if length(Y) != size(A,1) || length(B) != size(A,2)
        throw(DimensionMismatch("Matrices have wrong dimensions"))
    end

    H=A.H
    m=A.m

    M=size(H,1)
    b=pointer(B)
    y=pointer(Y)
    h=pointer(H)
    st=stride(H,2)

    sz=sizeof(T)

    BLAS.blascopy!(m,B,1,Y,1)

    for k=1:min(size(H,2),m-M+1)
        wp=h+sz*st*(k-1)
        yp=y+sz*(k-1)

        dt=dot(M,wp,1,yp,1)
        BLAS.axpy!(M,-2*dt,wp,1,yp,1)
    end

    for k=m-M+2:size(H,2)
        p=k-m+M-1

        wp=h+sz*st*(k-1)
        yp=y+sz*(k-1)

        dt=dot(M-p,wp,1,yp,1)
        BLAS.axpy!(M-p,-2*dt,wp,1,yp,1)
    end
    Y
end

# Each householder is symmetyric, this just reverses the order of application
function Base.A_mul_B!{T<:BlasFloat}(Y::Vector{T},A::BandedQ{T},B::Vector{T})
    if length(Y) != size(A,1) || length(B) != size(A,2)
        throw(DimensionMismatch("Matrices have wrong dimensions"))
    end

    H=A.H
    m=A.m

    M=size(H,1)
    b=pointer(B)
    y=pointer(Y)
    h=pointer(H)
    st=stride(H,2)

    sz=sizeof(T)

    BLAS.blascopy!(m,B,1,Y,1)

    for k=size(H,2):-1:m-M+2
        p=k-m+M-1

        wp=h+sz*st*(k-1)
        yp=y+sz*(k-1)

        dt=dot(M-p,wp,1,yp,1)
        BLAS.axpy!(M-p,-2*dt,wp,1,yp,1)
    end

    for k=min(size(H,2),m-M+1):-1:1
        wp=h+sz*st*(k-1)
        yp=y+sz*(k-1)

        dt=dot(M,wp,1,yp,1)
        BLAS.axpy!(M,-2*dt,wp,1,yp,1)
    end

    Y
end

function Base.A_mul_B!(Y::Matrix,A::BandedQ,B::Matrix)
    for j=1:size(A,2)
        Y[:,j]=A*B[:,j]
    end
    Y
end

function Base.Ac_mul_B!(Y::Matrix,A::BandedQ,B::Matrix)
    for j=1:size(A,2)
        Y[:,j]=A'*B[:,j]
    end
    Y
end

Base.full(A::BandedQ) = A*eye(eltype(A),size(A,1))

Base.linearindexing{T}(::Type{BandedQ{T}}) = Base.LinearSlow()
Base.getindex(A::BandedQ,k::Int,j::Int) = (A*eltype(A)[zeros(j-1);1.0;zeros(size(A,2)-j)])[k]

function Base.getindex(QR::BandedQR,d::Symbol)
    d == :Q && return BandedQ(QR.H,size(QR,1))
    d == :R && return QR.R

    throw(KeyError(d))
end

function Base.qr(A::BandedMatrix)
    QR = qrfact(A)
    QR[:Q],QR[:R]
end

function Base.qrfact(A::BandedMatrix)
    R=bzeros(eltype(A),size(A,1),size(A,2),A.l,A.l+A.u)
    R.data[A.l+1:end,:]=A.data
    banded_qrfact!(R)
end

flipsign(x,y) = Base.flipsign(x,y)
flipsign(x::BigFloat,y::BigFloat) = sign(y)==1?x:(-x)
flipsign(x,y::Complex) = y==0?x:x*sign(y)

function banded_qrfact!(R::BandedMatrix)
    T=eltype(R)
    M=R.l+1   # number of diag+subdiagonal bands
    m,n=size(R)
    W=Array(T,M,(n<m?n:m-1))
    w=pointer(W)
    r=pointer(R.data)
    sz=sizeof(T)
    st=stride(R.data,2)
    stw=stride(W,2)

    for k=1:min(size(R,1)-R.l,n)
        v=r+sz*(R.u + (k-1)*st)    # diagonal entry
        wp=w+stw*sz*(k-1)          # k-th column of W
        BLAS.blascopy!(M,v,1,wp,1)
        W[1,k]+= flipsign(BLAS.nrm2(M,wp,1),W[1,k])
        normalize!(M,wp)

        for j=k:min(k+R.u,n)
            v=r+sz*(R.u + (k-1)*st + (j-k)*(st-1))
            dt=dot(M,wp,1,v,1)
            BLAS.axpy!(M,-2*dt,wp,1,v,1)
        end
    end

    for k=m-R.l+1:(n<m?n:m-1)
        p=k-m+R.l
        v=r+sz*(R.u + (k-1)*st)    # diagonal entry
        wp=w+stw*sz*(k-1)          # k-th column of W
        BLAS.blascopy!(M-p,v,1,wp,1)
        W[1,k]+= flipsign(BLAS.nrm2(M-p,wp,1),W[1,k])
        normalize!(M-p,wp)

        for j=k:min(k+R.u,n)
            v=r+sz*(R.u + (k-1)*st + (j-k)*(st-1))
            dt=dot(M-p,wp,1,v,1)
            BLAS.axpy!(M-p,-2*dt,wp,1,v,1)
        end
    end

    BandedQR(W,BandedMatrix(R.data[1:R.u+1,:],m,0,R.u))
end


text:
_tofile(f::AbstractString) = h5open(f, "r+")
_tofile(f::DataFile) = f

get_chunk(arr::HDF5Array{<:Any, <:Any, <:Any, Nothing}) = error("given HDF5Array is not chunked")
get_chunk(arr::HDF5Array{<:Any, N, HDF5Dataset, NTuple{N, Int}}) where {N} = arr.chunk

function HDF5Array(dset::HDF5Dataset)
    try
        c = get_chunk(dset)
        return HDF5Array{eltype(dset), ndims(dset), HDF5Dataset, typeof(c)}(dset, c)
    catch ErrorException
        return HDF5Array{eltype(dset), ndims(dset), HDF5Dataset, Nothing}(dset)
    end
end
HDF5Array(file::Union{AbstractString, DataFile}, path::AbstractString) =
    HDF5Array(_tofile(file)[path])

HDF5Array{T, N, D, Nothing}(data) where {T, N, D} = HDF5Array{T, N, D, Nothing}(data, nothing)

_make_dim_args(
    T::Type, dims::NTuple{N, Integer},
    max_dims::NTuple{N, Integer}, chunk::NTuple{N, Integer}
) where {N} = (T, (dims, max_dims), "chunk", chunk)

_make_dim_args(
    T::Type, dims::NTuple{N, Integer},
    max_dims::Nothing, chunk::Nothing
) where {N} = (datatype(T), dataspace(dims...))

_make_dim_args(
    T::Type, dims::NTuple{N, Integer},
    max_dims::Nothing, chunk::NTuple{N, Integer}
) where {N} = (datatype(T), dataspace(dims...), "chunk", chunk)

_convert_chunk_type(c::NTuple{N, Integer}) where {N} = NTuple{N, Int}(c)
_convert_chunk_type(::Nothing) = nothing

_make_filter_args(compress::Integer) = ("compress", compress)
_make_filter_args(::Nothing) = ()

function HDF5Array{T, N}(
    file::Union{AbstractString, DataFile}, path::AbstractString,
    ::UndefInitializer, dims::Vararg{Integer, N};
    chunk::Union{NTuple{N, Integer}, Nothing} = nothing,
    max_dims::Union{NTuple{N, Integer}, Nothing} = nothing,
    compress::Union{Integer, Nothing} = nothing
) where {T, N}
    c = _convert_chunk_type(chunk)
    return HDF5Array{T, N, HDF5Dataset, typeof(c)}(
        d_create(
            _tofile(file), path,
            _make_dim_args(T, dims, max_dims, chunk)...,
            _make_filter_args(compress)...
        ),
        c
    )
end


text:
using LinearAlgebra
using Test

using QD

function approxeq(::Type{T}, tolfac, xs, r0, r1) where {T}
    tol = max(abs.(big.(xs))..., abs(r0), abs(r1))
    iseq = abs(r0 - r1) <= tolfac * tol * big(eps(T))
    if !iseq
        @show T tolfac xs r0 r1 abs(r0 - r1)
    end
    iseq
end

include("test128.jl")
include("test256.jl")


text:
struct Basis
    Ψ::Matrix
    r::Int
end

function VandermondeBasis(x,r::Int;verbose::Bool=false)
    """Makes a vandermonde basis given a range and a number of basis modes"""
    Ψ = zeros((length(x),r))
    for (j,col) in enumerate(eachcol(Ψ))
        ψ = make_vander_column(j,x)
        if verbose
            println("ψ","_"*string(j)," size: ",size(ψ))
        end
        Ψ[:,j] = make_vander_column(j,x) #Each col is x^(k-1)
    end
    if verbose
        print("Final Ψ size: ",size(Ψ))
    end
    return Basis(Ψ,r) #Truncated basis matrix
end

function make_vander_column(k,col_k)
    ψ = [x^(k-1) for x in col_k]
    return ψ #Column of basis matrix
end


text:
using Revise, WaltonForaging

## Leaving times learning rate
gr(size = (600,600))
patch = ["poor", "medium","rich"]
travel = ["long", "short"]
stim = [false, true]
LR_SimulatedStim  =  DataFrame(patch = String[], travel = String[], stim = Bool[], time = Float64[],
    color = Symbol[], line = Symbol[], group = String[])
for p in patch
    pv = p == "rich" ? 3 : p == "medium" ? 2 : 1
    for t in travel
        tv = t == "short" ? -1 : 1
        l = t == "short" ? :solid : :dash
        tgroup = t =="short" ? "Short Travel" : "Long Travel"
        for s in stim
            if s
                pre = p == "rich" ? 3 : p == "medium" ? 2 : 1
                sv = 1.3 - pre*(0.4)
                c = :blue
                sgroup = " Stim"
            else
                sv = 0
                c = :black
                sgroup = " No stim"
            end
            value = 3 * pv + tv + sv
            label = tgroup * sgroup
            push!(LR_SimulatedStim, (patch = p, travel = t, stim = s, time = value, color = c, line = l, group = label))
        end
    end
end
@df LR_SimulatedStim plot(:patch,:time, group = :group,
    legend = :topleft, color = :color, linestyle = :line,
    ylabel = "Patch leaving time (simulated)",
    xlabel = "Patch Richness", ylims = (0,11))
##
folder = "/Users/dariosarra/Dropbox (Personal)/2022_BBSRC_Foraging/Figures/Illustrative_Serotonin_effects"
savefig(joinpath(folder,"Bhv_learningrate1.pdf"))
##
lr_df1 = combine(groupby(LR_SimulatedStim,[:patch, :stim]), :time => mean)
change_lr_patch=unstack(lr_df1, :stim, :time_mean)
change_lr_patch[!,:change] = change_lr_patch[:,3] - change_lr_patch[:,2]
change_lr_patch[!,:color] .= :blue
change_lr_patch[!,:Hypothesis] .= "Local"
@df change_lr_patch bar(:patch, :change,
    ylabel = "Stimulated change in patch leving time", xlabel = "Patch Richness",
    color = :blue, legend = false, ylims = (0,1))
    # @df change_Stim_patch plot!(:patch, :change,linestyle = :dot, color = :grey)
##
savefig(joinpath(folder,"Bhv_learningrate2.pdf"))
##
lr_df2 = combine(groupby(LR_SimulatedStim,[:travel, :stim]), :time => mean)
change_lr_trav=unstack(lr_df2, :stim, :time_mean)
change_lr_trav[!,:change] = change_lr_trav[:,3] - change_lr_trav[:,2]
change_lr_trav[!,:color] .= :blue
change_lr_trav[!,:Hypothesis] .= "Local"
@df change_lr_trav bar(:travel, :change,
    ylabel = "Stimulated change in patch leving time",
    xlabel = "Travel distance", color = :blue, legend = false, ylims = (0,1))
    # @df change_Stim_patch plot!(:travel, :change,linestyle = :dot, color = :grey)
##
savefig(joinpath(folder,"Bhv_learningrate3.pdf"))
##
patch = ["poor","rich"]
time = [0,1,2,3]
stim = [false, true]
SimulatedPhoto  =  DataFrame(time = Float64[], patch = String[], stim = Bool[], RPE = Float64[],
    color = Symbol[], line = Symbol[], group = String[])
for t in time
    for p in patch
        pv = p == "poor" ? 1.2 : 0.8
        c = p == "poor" ? :brown : :green
        pgroup = p == "poor" ? "Poor patch" : "Rich patch"
        for s in stim
            sv = s ? 0.85 : 1
            l = s ? :dash : :solid
            sgroup = s ? " Stim" : " No stim"
            value = (pv*sv)*t + 1
            label = pgroup * sgroup
            push!(SimulatedPhoto,(time = t, patch = p, stim = s, RPE = value,
                color = c, line = l, group = label))
        end
    end
end
@df SimulatedPhoto plot(:time, :RPE, group = :group,
    legend = :topleft, color = :color, linestyle = :line,
    ylabel = "RPE response at reward (simulated)",
    xlabel = "Elapsed time")
#"u expect reward to have not reduce as much"
##
savefig(joinpath(folder,"DA_learningrate1.pdf"))
##
df1 = combine(groupby(SimulatedPhoto,[:patch, :stim]), :RPE => mean)
change_Stim_patch = unstack(df1, :stim, :RPE_mean)
change_Stim_patch[!,:change] = change_Stim_patch[:,3] - change_Stim_patch[:,2]
@df change_Stim_patch bar(:patch, abs.(:change),
    ylabel = "Stimulated RPE absolute magnitude change", xlabel = "Patch Richness",
    color = :blue, legend = false)
    # @df change_Stim_patch plot!(:patch, :change,linestyle = :dot, color = :grey)
##
savefig(joinpath(folder,"DA_learningrate2.pdf"))
##


text:
module DirichletProcess

using Distributions

export DP, dp_sample!
export HDP, hdp_sample!

# Write your package code here.
mutable struct DP
    base_measure::Function
    alpha::Float64
    cache
    weights
    total_stick_used::Float64
    DP(base_measure, alpha) = new(base_measure, alpha, [], [], 0.0)
end

mutable struct HDP
    base_measure::Function
    alpha1::Float64
    alpha2::Float64
    dp1::DP
    dp2::DP
    function HDP(base_measure, alpha1, alpha2)
        dp1 = DP(base_measure, alpha1)
        dp2 = DP(() -> dp_sample!(dp1), alpha2)
        return new(base_measure, alpha1, alpha2, dp1, dp2)
    end
end

function dp_sample!(dp::DP)
    remaining = 1.0 - dp.total_stick_used
    i = rand(Categorical([dp.weights..., remaining]))
    if i <= length(dp.weights)
        return dp.cache[i]
    else
        stick_piece = rand(Beta(1, dp.alpha)) * remaining
        dp.total_stick_used += stick_piece
        push!(dp.weights, stick_piece)
        new_value = dp.base_measure()
        push!(dp.cache, new_value)
        return new_value
    end
end

function hdp_sample!(hdp::HDP)
    return dp_sample!(hdp.dp2)
end

end


text:
@testset "ilog2 $U(2^$e)" for e in (0, 7, 63, 64, 128, 200), U in (UInt8, UInt16, UInt32, UInt64, UInt128)

    a = U(2)^e
    @test iszero(a) || ilog2(a) == e
    @test iszero(a) || ilog2(a - 1) == e - 1
end

@testset "ilog2 $U(2^$e)" for e in (0, 7, 63, 64, 128, 200), U in (Int8, Int16, Int32, Int64, Int128, BigInt)

    a = U(2)^e
    @test a <= 0 || ilog2(a) == e
    @test a <= 0 || ilog2(a - 1) == e - 1
    @test a <= 0 || ilog2(-a) == e
    @test a <= 0 || ilog2(-a + 1) == e - 1
end

@testset "ilog2 $U(0)" for U in (UInt8, UInt16, UInt32, UInt64, UInt128, BigInt)
    a = zero(U)
    @test ilog2(a) == -1
    @test ilog2(signed(a)) == - 1
end


text:
using DistStat, Random, LinearAlgebra

mutable struct MultUpdate
    maxiter::Int
    step::Int
    verbose::Bool
    tol::Real
    lambda::Real
    function MultUpdate(; maxiter::Int=100, step::Int=10, verbose::Bool=false, tol::Real=1e-10, lambda::Real=2e-38)
        maxiter > 0 || throw(ArgumentError("maxiter must be greater than 0."))
        tol > 0 || throw(ArgumentError("tol must be positive."))
        lambda >= 0 || throw(ArgumentError("lambda must be nonnegative."))
        new(maxiter, step, verbose, tol, lambda)
    end
end

mutable struct NMFVariables{T, A}
    # NMF variables corresponding to X::MPIMatrix{T,A} # m x [n], X is external.
    m::Int # rows
    n::Int # cols
    r::Int # intermediate
    Vt::MPIMatrix{T,A} # r x [m]
    Vt_prev::MPIMatrix{T,A} # r x [m]
    W::MPIMatrix{T,A} # r x [n]
    W_prev::MPIMatrix{T,A} # r x [n]
    tmp_rr::A # broadcasted
    tmp_rm1::MPIMatrix{T,A}
    tmp_rm2::MPIMatrix{T,A}
    tmp_rn1::MPIMatrix{T,A}
    tmp_rn2::MPIMatrix{T,A}
    tmp_mn::Union{MPIMatrix, Nothing}
    tmp_rm_local::A
    function NMFVariables(X::MPIMatrix{T,A}, r; eval_obj=false, seed=nothing) where {T,A}
        m, n    = size(X)
        Vt      = MPIMatrix{T,A}(undef, r, m)
        Vt_prev = MPIMatrix{T,A}(undef, r, m)
        W       = MPIMatrix{T,A}(undef, r, n)
        W_prev  = MPIMatrix{T,A}(undef, r, n)
        rand!(Vt; seed=seed, common_init=true)
        if seed ≠ nothing
            seed += 1
        end
        rand!(W; seed=seed, common_init=true)
        tmp_rr  = A{T}(undef, r, r)
        tmp_rm1 = MPIMatrix{T,A}(undef, r, m)
        tmp_rm2 = MPIMatrix{T,A}(undef, r, m)
        tmp_rn1 = MPIMatrix{T,A}(undef, r, n)
        tmp_rn2 = MPIMatrix{T,A}(undef, r, n)
        if eval_obj
            tmp_mn = MPIMatrix{T,A}(undef, m, n)
        else
            tmp_mn = nothing
        end
        tmp_rm_local = A{T}(undef, r, m)
        new{T,A}(m, n, r, Vt, Vt_prev, W, W_prev, tmp_rr, tmp_rm1, tmp_rm2, tmp_rn1, tmp_rn2, tmp_mn, tmp_rm_local)
    end
end

function reset!(v::NMFVariables{T,A}; seed=nothing) where {T,A}
    rand!(v.Vt; seed=seed, common_init=true)
    if seed ≠ nothing
        seed += 1
    end
    rand!(v.W; seed=seed, common_init=true)
end

function update_V!(X::MPIArray, u::MultUpdate, v::NMFVariables)
    WXt = mul!(v.tmp_rm1, v.W, transpose(X); tmp=v.tmp_rm_local)
    WWt = mul!(v.tmp_rr, v.W, transpose(v.W))
    WWtVt = mul!(v.tmp_rm2, WWt, v.Vt)
    v.Vt .= v.Vt .* WXt ./ (WWtVt .+ u.lambda)
end

function update_W!(X::MPIArray, u::MultUpdate, v::NMFVariables)
    VtX = mul!(v.tmp_rn1, v.Vt, X; tmp=v.tmp_rm_local)
    VtV = mul!(v.tmp_rr, v.Vt, transpose(v.Vt))
    VtVW = mul!(v.tmp_rn2, VtV, v.W)
    v.W .= v.W .* VtX ./ (VtVW .+ u.lambda)
end

function nmf_get_objective!(X::MPIArray, u::MultUpdate, v::NMFVariables{T,A}) where {T,A}
    if v.tmp_mn != nothing
        mul!(v.tmp_mn, transpose(v.Vt), v.W; tmp=v.tmp_rm_local) # TODO: improve: print amount of update, etc.
        v.tmp_mn .= (v.tmp_mn .- X).^ 2
        return false, (sum(v.tmp_mn))
    else
        v.Vt_prev .= abs.(v.Vt_prev .- v.Vt)
        v.W_prev  .= abs.(v.W_prev  .- v.W)
        return false, (maximum(v.Vt_prev), maximum(v.W_prev))
    end
end

function nmf_mult_one_iter!(X::MPIArray, u::MultUpdate, v::NMFVariables)
    copyto!(v.W_prev, v.W)
    copyto!(v.Vt_prev, v.Vt)
    update_V!(X, u, v)
    update_W!(X, u, v)
end

function loop!(X::MPIArray, u, iterfun, evalfun, args...)
    converged = false
    t = 0
    while !converged && t < u.maxiter
        t += 1
        iterfun(X, u, args...)
        if t % u.step == 0
            converged, monitor = evalfun(X, u, args...)
            if DistStat.Rank() == 0
                println(t, ' ', monitor)
            end
        end
    end
end

function nmf!(X::MPIArray, u::MultUpdate, v::NMFVariables)
    loop!(X, u, nmf_mult_one_iter!, nmf_get_objective!, v)
end

include("cmdline.jl")
opts = parse_commandline_nmf()
if DistStat.Rank() == 0
    println("world size: ", DistStat.Size())
    println(opts)
end

m = opts["rows"]
n = opts["cols"]
r = opts["r"]
iter = opts["iter"]
interval = opts["step"]
T = Float64
A = Array
if opts["gpu"]
    using CUDA
    A = CuArray
end
if opts["Float32"]
    T = Float32
end
init_opt = opts["init_from_master"]
seed = opts["seed"]
eval_obj = opts["eval_obj"]

X = MPIMatrix{T,A}(undef, m, n)
rand!(X; common_init=init_opt, seed=0)
uquick = MultUpdate(;maxiter=2, step=1, verbose=true)
u = MultUpdate(;maxiter=iter, step=interval, verbose=true)
v = NMFVariables(X, r; eval_obj=eval_obj, seed=seed)
nmf!(X, uquick, v)
reset!(v; seed=seed)
if DistStat.Rank() == 0
    @time nmf!(X, u, v)
else
    nmf!(X, u, v)
end


text:
module CUFFT

using CUDAapi

using ..CuArrays
using ..CuArrays: libcufft
import ..CuArrays: unsafe_free!

using CUDAdrv
using CUDAdrv: CUstream

using CEnum

include("libcufft_common.jl")
include("error.jl")

version() = VersionNumber(cufftGetProperty(CUDAapi.MAJOR_VERSION),
                          cufftGetProperty(CUDAapi.MINOR_VERSION),
                          cufftGetProperty(CUDAapi.PATCH_LEVEL))

include("libcufft.jl")

include("util.jl")
include("wrappers.jl")
include("fft.jl")

end


text:
function optkeywords()
  #this will hold the arguments
  optkwargs::Dict{Symbol,Real} = Dict()
  kws = [:f_tol, :g_tol, :x_tol, :iterations, :Δ]
  sizehint!(optkwargs, length(kws))

  paramkw(s) = Symbol(PARAM[:optprefix], :_ , s)

  for kw ∈ kws
    arg = PARAM[paramkw(kw)]
    if !(arg===nothing)
      push!(optkwargs, kw=>arg)
    end
  end

  return optkwargs
end

function getoptbounds(zs::ZSpec,
  Θ::VolumePartsDLSq, dims::NamedTuple,
  PredictionType::Val = Val(PARAM[:predictiontype]);
  minvol::Float64=10^-6)

  Abounds = PARAM[:abounds]
  V₀bounds = PARAM[:v0bounds]

  #this willl hold the bounds
  ZAbounds::Matrix{Tuple{Float64,Float64}} = Matrix{Tuple{Float64,Float64}}(undef, dims.K, dims.T)
  ZV₀bounds::Vector{Tuple{Float64,Float64}}= Vector{Tuple{Float64,Float64}}(undef, dims.T-1)

  function no0log(x)
    (x<0) && error("no0log cannot have a negative input (did you use a negative V₀ bound?)")
    return log(max(minvol,x))
  end

  #need to scale the bounds at the same level as the coefficients
  σᵥ::Float64 = zs.scale[zs.Zms.Fvol]
  scaledV₀bounds = (t-> no0log(t / σᵥ)).(V₀bounds)

  #scale the intercept term
  for t ∈ 1:(dims.T-1)
    ZV₀bounds[t] = scaledV₀bounds
  end

  #scale the A matrix coefficeint bounts
  for (k, FZcoef) ∈ enumerate(zs.Zms.Fξs)
    FZw::Symbol = zs.Zms[FZcoef].X[:Fw]
    σ_kw = zs.scale[FZw]
    scaledAbounds = (t->t*σ_kw/σᵥ).(Abounds)
    for t ∈ 1:dims.T
      ZAbounds[k,t] =  scaledAbounds
    end
  end

  #combine the bounds
  bounds=[vec(ZAbounds); ZV₀bounds]

  return bounds
end

#generates an in-place prediction function
function leastsquaresobj(Θ::AbstractVolumeParts, X::NamedTuple, v::AbstractVector,
  PredictionType::Val)

  function obj!(out::AbstractVector, Θ::AbstractVolumeParts)
    out .= Θ(X,v,PredictionType) .- v
    return nothing
  end

  return obj!
end

function optmodel(panel::AbstractDataFrame,  zs::ZSpec, ::Val{:levellsq},
    ::Type{T}=PARAM[:opttype],
    ::Type{TV}=Vector{T},
    ::Type{TM}=Matrix{T},
    PredictionType::Val = Val(PARAM[:predictiontype]),
    optkwargs = optkeywords(),
    ) where {
      T<:Real, TV<:AbstractVector{T}, TM<:AbstractMatrix{T}}
  PARAM[:predictiontype] ∈ [:level, :levellog] || error(
    "invalid predictiontype $PredictionType for model type $(PARAM[:fluxmodeltype])")
  (T == Float64) || @warn("T=$T, which may lead to numerical instability.
    Is this really what you want?")

  Zms::MeasureSpec = zs.Zms
  Fw::Vector{Symbol} = Fweights(Zms, :Fw)
  FRLw::Vector{Symbol} = Fweights(Zms, :FRLw)
  essentialfields = [Fw; FRLw; Zms.Fvol; ]

  #form the data into a more usable format and expand the arrays
  p = modelparts(panel, Zms,T,TV,TM,
    weightdata = (Fw=Fw, FRLw=FRLw),
    Fvol=Zms.Fvol)
  dims = p.dims
  tidx::Dict = p.tidx
  cleanpanel::SubDataFrame = p.cleanpanel
  Ncleanpanel::Int = nrow(cleanpanel)
  ws::TM = p.dat[:Fw]
  RLws::TM = p.dat[:FRLw]
  ts::Vector{Int} = p.ts
  X = (ws=ws, RLws=RLws)
  v::TV = p.dat[:Fvol]

  #construct the parameter
  Θ = VolumePartsDLSq(dims.T, dims.K, ts, T)

  J, J! = JVolumePartsDLSq(Θ, X, v, T, PredictionType)
  testjacobian(Θ, X, v, T, J = J, J! = J!, dims=dims)
  obj! = leastsquaresobj(Θ, X, v, PredictionType)

  out = similar(v)
  obj!(out,Θ)
  vhatver = Θ(X,v,PredictionType)
  @assert out .+ v ≈ vhatver

  #this is not used in the optimization, but keeps track of our objective
  function loss(Θ::AbstractVolumeParts, X::NamedTuple,y::AbstractVector) where TX<:AbstractMatrix

    #keep count of each iteration
    absμ = Θ(X, y, PredictionType)

    sumsq = sum((absμ .- y) .^ 2)

    return sumsq
  end

  @info "Beginning optimization with $(optkwargs)"
  opt = optimize!(LeastSquaresOptim.LeastSquaresProblem(x=Θ, f! = obj!, g! = J!, J=J, output_length=Ncleanpanel),
    LevenbergMarquardt(LeastSquaresOptim.LSMR()); optkwargs...)

  @info "opt: \n$opt"
  λ⁺ = loss(Θ, X,v)

  #compute the full jacobian with AD
  Θver = VolumePartsLSq(dims.T, dims.K, ts, T)
  Θver .= Θ
  ∇Θver = gradient(()->loss(Θver,X,v),Flux.params(Θver))
  rmsver = sqrt(mean(vec∇(Θver, ∇Θver).^2))
  @info "lst squares gradient rms: $rmsver"

  #=J!(J,Θ)

  Jcolsums = vec(sum(J, dims=1))
  @info "opt:\n$opt"
  @info "Jcolsums: $Jcolsums"
  @info "L₁(Jcolsums)=$(sum(abs.(Jcolsums)))"=#

  results::DataFrame = formresults(panel, Zms, Θ, PredictionType)

  return Int(round(λ⁺)), results
end

function optmodel(panel::AbstractDataFrame,  zs::ZSpec,::Val{:levellsqbb},
    ::Type{T}=PARAM[:opttype],
    ::Type{TV}=PARAM[:optgpu] ? CUDA.CuVector{T} : Vector{T},
    ::Type{TM}=PARAM[:optgpu] ? CUDA.CuMatrix{T} : Matrix{T},
    PredictionType::Val = Val(PARAM[:predictiontype]);
    evalmodels::Bool = PARAM[:optbbcomparemodels],
    ) where {
      T<:Real, TV<:AbstractVector{T}, TM<:AbstractMatrix{T}}
  PARAM[:predictiontype] ∈ [:level, :levellog] || error(
    "invalid predictiontype $PredictionType for model type $(PARAM[:fluxmodeltype])")
  (T == Float64) || @warn("T=$T, which may lead to numerical instability.
    Is this really what you want?")
  CUDA.allowscalar(false)

  Zms::MeasureSpec = zs.Zms
  Fw::Vector{Symbol} = Fweights(Zms, :Fw)
  FRLw::Vector{Symbol} = Fweights(Zms, :FRLw)
  essentialfields = [Fw; FRLw; Zms.Fvol; ]

  #form the data into a more usable format and expand the arrays
  p = modelparts(panel, Zms,T,TV,TM,
    weightdata = (Fw=Fw, FRLw=FRLw),
    Fvol=Zms.Fvol)
  dims = p.dims
  tidx::Dict = p.tidx
  cleanpanel::SubDataFrame = p.cleanpanel
  Ncleanpanel::Int = nrow(cleanpanel)
  ws::TM = p.dat[:Fw]
  RLws::TM = p.dat[:FRLw]
  ts::Vector{Int} = p.ts
  X = (ws=ws, RLws=RLws)
  v::TV = p.dat[:Fvol]

  #construct the parameter
  Θ = VolumePartsDLSq(dims.T, dims.K, ts, T)

  testprediction(Θ, X, v, PredictionType)

  #least squares objective funciton
  function loss(Θ::AbstractVolumeParts, X::NamedTuple,y::AbstractVector)

    #keep count of each iteration
    absμ = Θ(X, y, PredictionType)

    sumsq = sum((absμ .- y) .^ 2)

    return sumsq
  end

  function loss(Θ::AbstractVolumeParts, X::NamedTuple,y::CuVector{T})

    #keep count of each iteration
    absμ = Θ(X, y, PredictionType)

    krn(absμᵢ,yᵢ) = (absμᵢ - yᵢ)*(absμᵢ - yᵢ)
    sumsq = CUDA.sum(krn.(absμ,y))

    return sumsq
  end

  #this just puts the problem in the form used by bboptimize
  function loss(Π::AbstractVector)
    Θ.Π .= TV(Π)
    return loss(Θ, X, v)
  end

  bounds = getoptbounds(zs, Θ, dims, PredictionType)

  if evalmodels
    @warn "Testing models! This could take a while...."
    optbbcomparemodels(
      bounds=bounds,
      loss=loss,
      Θ=Θ,
      X=X,
      y=v,)
  end

  opt, Ψ = optbbestimate(
    bounds=bounds,
    loss=loss,
    Θ=Θ,
    X=X,
    y=v,
  )

  λ⁺ = loss(Θ, X,v)

  #const a differentiable version with nromal data types of parameters
  Θver = VolumePartsLSq(dims.T, dims.K, Θ.ts, T)
  Θver .= Vector(Θ.Π)
  @assert loss(Θver,X,v) ≈ λ⁺

  optbbresults(loss=loss, opt=opt, Ψ=Ψ, Θ=Θ, Θver=Θver, X=X, y=v)
  results::DataFrame = formresults(panel, Zms, Θver, PredictionType)

  return Int(round(min(λ⁺, typemax(Int)÷10))), results
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Test
using ACL

tests = ["floorsum"]

@testset "ACL" begin

for testname in tests
    @testset "$testname" begin
        include("$testname.jl")
    end
end

end


text:
struct Callback{F, S}
    func :: F
    schedule :: S
end

(callback::Callback)(sim) = callback.func(sim)

"""
    Callback(func; schedule=IterationInterval(1))

Return `Callback` that executes `func(sim::Simulation)` on `schedule`.
"""
Callback(func; schedule=IterationInterval(1)) = Callback(func, schedule)

#####
##### Utilities for run!
#####

default_callback_name(n) = Symbol(:callback, n)

function add_callbacks!(sim, callbacks)
    n_existing_callbacks = length(sim.callbacks)

    for (i, callback) in enumerate(callbacks)
        name = default_callback_name(n_existing_callbacks + i)
        sim.callbacks[name] = callback
    end

    return nothing
end


text:
"""
    get_d_hsml_3D( dx::Real, dy::Real, dz::Real,
                   hsml_inv::Real )

Computes the distance in 3D to the pixel center in units of the kernel support.
"""
function get_d_hsml( dx::T, dy::T, dz::T,
                     hsml_inv::T) where T
    √( dx*dx + dy*dy + dz*dz ) * hsml_inv
end

"""
    function calculate_weights_3D(  wk::Array{<:Real,1},
                                    iMin::Integer, iMax::Integer,
                                    jMin::Integer, jMax::Integer,
                                    kMin::Integer, kMax::Integer,
                                    x::Real, y::Real, z::Real,
                                    hsml::Real, hsml_inv::Real,
                                    kernel::SPHKernel,
                                    x_pixels::Integer, y_pixels::Integer )

Calculates the kernel- and geometric weights of the pixels a particle contributes to.
"""
@fastmath function calculate_weights( wk::Vector{Float64},
                                    iMin::Integer, iMax::Integer,
                                    jMin::Integer, jMax::Integer,
                                    kMin::Integer, kMax::Integer,
                                    x::T, y::T, z::T,
                                    hsml::T, hsml_inv::T,
                                    kernel::SPHKernel,
                                    x_pixels::Integer, y_pixels::Integer ) where T

    is_undersampled = false

    if hsml <= 1.0
        is_undersampled = true
    end

    distr_weight = 0.0
    n_distr_pix  = 0

    @inbounds for i = iMin:iMax
        x_dist, dx = get_x_dx(x, hsml, i)

        for j = jMin:jMax
            y_dist, dy = get_x_dx(y, hsml, j)

            for k = kMin:kMax
                z_dist, dz = get_x_dx(z, hsml, k)

                idx = calculate_index(i, j, k, x_pixels, y_pixels)

                dxdydz        = dx * dy * dz

                if is_undersampled

                    wk[idx]       = dxdydz
                    distr_weight += dxdydz
                    n_distr_pix  += 1

                    continue

                else # is_undersampled

                    u = get_d_hsml(x_dist, y_dist, z_dist, hsml_inv)

                    if u <= 1.0

                        wk[idx]       =  𝒲₃(kernel, u, hsml_inv)
                        wk[idx]      *= dxdydz
                        distr_weight += wk[idx]
                        n_distr_pix  += 1

                    else
                        wk[idx] = 0.0
                    end # u < 1.0

                end # is_undersampled
            end # k
        end # j
    end # i

    return wk, n_distr_pix, distr_weight
end

"""
    function calculate_index(i::Integer, j::Integer, x_pixels::Integer)

Calculates the index of a flattened 3D image array.
"""
function calculate_index(  i::T, j::T, k::T, x_pixels::T, y_pixels::T) where T
    return floor(T, i * x_pixels + j * y_pixels + k) + 1
end

"""
    get_quantities_3D( pos, weight, hsml,
                            rho, m, len2pix::T) where T

Helper function to get
"""
function get_quantities_3D( pos, weight, hsml,
                            rho, m, len2pix::T) where T

        hsml     = T(hsml * len2pix)
        hsml_inv = T(1.0/hsml)

        volume  = T(4π/3) * hsml^3
        rho    /= T(len2pix)^3
        dz      = T(m / rho / volume)

    return T.(pos), T(weight), hsml, hsml_inv, volume, dz
end

"""
   sphMapping_3D( Pos::Array{<:Real}, HSML::Array{<:Real},
                  M::Array{<:Real}, Rho::Array{<:Real},
                  Bin_Q::Array{<:Real}, Weights::Array{<:Real}=ones(length(Rho));
                  param::mappingParameters, kernel::SPHKernel,
                  show_progress::Bool=false )

Underlying function to map SPH data to a 3D grid.
"""
function sphMapping_3D( Pos::Array{T}, HSML::Array{T},
                        M::Array{T}, Rho::Array{T},
                        Bin_Q::Array{T}, Weights::Array{T};
                        param::mappingParameters, kernel::SPHKernel,
                        show_progress::Bool=false ) where T

    N = size(M,1)  # number of particles

    # max number of pixels over which the particle can be distributed
    N_distr = param.Npixels[1] * param.Npixels[2] * param.Npixels[3]

    image = zeros(Float64, N_distr,2)

    halfXsize = param.halfsize[1]
    halfYsize = param.halfsize[2]
    halfZsize = param.halfsize[3]

    if param.periodic
        k_start = 0
    else
        k_start = 7
    end

    # allocate arrays for weights
    wk = zeros(Float64, N_distr)

    # allocate array for positions
    #pos = Vector{eltype(Pos[1,1])}(undef, 3)

    if show_progress
        P = Progress(N)
    end

    # loop over all particles
    @inbounds for p = 1:N

        bin_q = Float64(Bin_Q[p])

        if bin_q == 0.0
            continue
        end

        _pos, weight, hsml, hsml_inv, volume, dz = get_quantities_3D(Pos[:,p], Weights[p], HSML[p], Rho[p], M[p], param.len2pix)

        for k_periodic = k_start:7

            x, y, z, skip_k = get_xyz( _pos, HSML[p], k_periodic, param)

            if skip_k
                continue
            end

            # calculate relevant pixel range
            iMin, iMax = pix_index_min_max( x, hsml, param.Npixels[1] )
            jMin, jMax = pix_index_min_max( y, hsml, param.Npixels[2] )
            kMin, kMax = pix_index_min_max( z, hsml, param.Npixels[3] )

            wk, n_distr_pix, distr_weight = calculate_weights( wk,
                                                               iMin, iMax,
                                                               jMin, jMax,
                                                               kMin, kMax,
                                                               x, y, z,
                                                               hsml, hsml_inv,
                                                               kernel,
                                                               param.Npixels[1],
                                                               param.Npixels[2])

            # skip if the particle is not contained in the image
            # ( should only happen with periodic boundary conditions )
            if n_distr_pix == 0
                continue
            end

            weight_per_pix = n_distr_pix / distr_weight
            kernel_norm    = volume / n_distr_pix
            volume_norm    = kernel_norm * weight_per_pix * weight * dz

            @inbounds for i = iMin:iMax, j = jMin:jMax, k = kMin:kMax

                idx = calculate_index( i, j, k,
                                       param.Npixels[1],
                                       param.Npixels[2] )

                image[idx,1], image[idx,2] = update_image( image[idx,1], image[idx,2],
                                                           wk[idx], bin_q, volume_norm)

            end # i, j, k

        end # k_periodic

         # update for ProgressMeter
        if show_progress
            next!(P)
        end
    end # p

    return image

end # function


text:
using AWS.EC2
using AWS

include("config.jl")

instances = ec2_instances_by_tag("ClusterName", ec2_clustername; running_only=false)
println("Starting previously stopped EC2 nodes : ", instances)

ec2_start(instances)
ec2_show_status(instances)


text:
function main()
    scanners = []
    for line in readlines()
        if occursin('s', line)
            push!(scanners, [])
        elseif occursin(',', line)
            push!(scanners[end], [parse.(Int, split(line, ','))...])
        end
    end
    queue = [scanners[1]]
    coords_lookup = Dict(scanners[1] => [0, 0, 0])
    while !isempty(queue)
        curr_scanner = pop!(queue)
        for other_scanner in scanners
            other_scanner ∈ keys(coords_lookup) && continue
            rotations = [(a, b, c) for a=1:3, b=1:3, c=1:3 if a != b != c != a]
            flips = [(nx, ny, nz) for nx=[-1, 1], ny=[-1, 1], nz=[-1, 1]]
            function try_all_perms()
                for (r1, r2, r3) ∈ rotations, flip ∈ flips
                    other_permuted = other_scanner .|> (v -> v .* flip) .|> (v -> [v[r1], v[r2], v[r3]])
                    diff_counter = Dict()
                    for other_vec in other_permuted
                        differences = curr_scanner .|> curr_vec -> curr_vec .- other_vec
                        for diff in differences
                            if (diff_counter[diff] = get!(diff_counter, diff, 0) + 1) >= 12
                                coords_lookup[other_scanner] = diff
                                for i in eachindex(other_scanner)
                                    other_scanner[i] = other_permuted[i] .+ diff
                                end
                                other_scanner ∉ queue && push!(queue, other_scanner)
                                return
                            end
                        end
                    end
                end
            end
            try_all_perms()
        end
    end
    println(length(Set(vcat(scanners...))))
    println([sum(abs.(a .- b)) for a=values(coords_lookup), b=values(coords_lookup)] |> maximum)
end
main()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
file = open(f->read(f, String), "day_12.txt")

regex = r"([A-Z])(\d+)"

input =[x.captures for x in eachmatch(regex,file)]

#converting to cartezian format
@NamedTuple{x::Int64,y::Int64,angle::Int64}

cartezinaize_comand = Dict()
cartezinaize_comand["N"] = (north(i) = (x=0,y=i,angle=0))
cartezinaize_comand["S"] = (south(i) = (x=0,y=-i,angle=0))
cartezinaize_comand["E"] = (east(i) =  (x=i,y=0,angle=0))
cartezinaize_comand["W"] = (west(i) =  (x=-i,y=0,angle=0))
cartezinaize_comand["L"] = (left(i) =  (x=0,y=0,angle=i))
cartezinaize_comand["R"] = (right(i) = (x=0,y=0,angle=-i))

angle_to_card = Dict([(0,"E"),(90,"N"),(180,"W"),(270,"S")])
function change_foward(letter,angle)
    if (letter != "F")
        return letter
    else

        angle = (angle % 360)#overflow angles
        angle += 360*(angle < 0)#fixing negative angle
        return angle_to_card[angle]
    end
end

#start moving
pos = (x=0,y=0,angle=0)

for (letter,value) in input
    #change F to equivalent command
    letter = change_foward(letter,copy(pos.angle))

    movement = cartezinaize_comand[letter](parse(Int64,value))

    newX = pos.x + movement.x
    newY = pos.y + movement.y
    newAngle = pos.angle + movement.angle
    global pos = (x=newX,y=newY,angle=newAngle)
end

print(abs(pos.x)+abs(pos.y))

##Question 2

matrix_rotation = Dict()
matrix_rotation[0] = [1 0;0 1]
matrix_rotation[90] = [0 -1;1 0]
matrix_rotation[180] = [-1 0;0 -1]
matrix_rotation[270] = [0 1;-1 0]

waypoint = [10;1]
pos = (x=0,y=0,angle=0)
for (letter,value) in input
    value = parse(Int64,value)
    if letter == "F"
        newX = pos.x + (waypoint[1] * value)
        newY = pos.y + (waypoint[2] * value)
        global pos = (x=newX,y=newY,angle=0)
    elseif letter == "N"
        global waypoint += [0;value]
    elseif letter == "S"
        global waypoint += [0;-value]
    elseif letter == "E"
        global waypoint += [value;0]
    elseif letter == "W"
        global waypoint += [-value;0]
    else #angle commands
        angle = if (letter == "L") value else -value end# left or right
        angle = (angle % 360)#overflow angles
        angle += 360*(angle < 0)#fixing negative angle

        #applying angle with matrix of rotation
        global waypoint = matrix_rotation[angle] * waypoint
    end

    print("\n $(pos)")
    print("\n $(waypoint)")
end

print("\n$(abs(pos.x)+abs(pos.y))")


text:
# Structure Types
# ---------------
"""Supertype for character structures."""
abstract type AbstractCharStructure end

"""Leaf node of a structure tree; this contains a single component."""
struct Component <: AbstractCharStructure
    component::Char
end
struct LeftRightStructure{L<:AbstractCharStructure, R<:AbstractCharStructure} <: AbstractCharStructure
    left::L
    right::R
end
struct TopBottomStructure{T<:AbstractCharStructure, B<:AbstractCharStructure} <: AbstractCharStructure
    top::T
    bottom::B
end
struct VerticalThirdsStructure{L<:AbstractCharStructure, M<:AbstractCharStructure, R<:AbstractCharStructure} <:
AbstractCharStructure
    left::L
    middle::M
    right::R
end
struct HorizontalThirdsStructure{T<:AbstractCharStructure, M<:AbstractCharStructure, B<:AbstractCharStructure} <:
AbstractCharStructure
    top::T
    middle::M
    bottom::B
end
abstract type NestedStructure{I<:AbstractCharStructure, O<:AbstractCharStructure} <: AbstractCharStructure end
struct CenterNestedStructure{I, O} <: NestedStructure{I, O}
    inside::I
    outside::O
end
struct BottomNestedStructure{I, O} <: NestedStructure{I, O}
    inside::I
    outside::O
end
struct TopNestedStructure{I, O} <: NestedStructure{I, O}
    inside::I
    outside::O
end
struct RightNestedStructure{I, O} <: NestedStructure{I, O}
    inside::I
    outside::O
end
struct BottomRightNestedStructure{I, O} <: NestedStructure{I, O}
    inside::I
    outside::O
end
struct BottomLeftNestedStructure{I, O} <: NestedStructure{I, O}
    inside::I
    outside::O
end
struct TopRightNestedStructure{I, O} <: NestedStructure{I, O}
    inside::I
    outside::O
end
struct OverlapStructure{F<:AbstractCharStructure, B<:AbstractCharStructure} <: AbstractCharStructure
    front::F
    back::B
end

# Structure Parsing
# -----------------
function parse(tokens)::AbstractCharStructure
    next = take!(tokens)
    if next == '⿰'
        LeftRightStructure(parse(tokens), parse(tokens))
    elseif next == '⿱'
        TopBottomStructure(parse(tokens), parse(tokens))
    elseif next == '⿲'
        VerticalThirdsStructure(parse(tokens), parse(tokens), parse(tokens))
    elseif next == '⿳'
        HorizontalThirdsStructure(parse(tokens), parse(tokens), parse(tokens))
    elseif next == '⿴'
        CenterNestedStructure(parse(tokens), parse(tokens))
    elseif next == '⿵'
        BottomNestedStructure(parse(tokens), parse(tokens))
    elseif next == '⿶'
        TopNestedStructure(parse(tokens), parse(tokens))
    elseif next == '⿷'
        RightNestedStructure(parse(tokens), parse(tokens))
    elseif next == '⿸'
        BottomRightNestedStructure(parse(tokens), parse(tokens))
    elseif next == '⿹'
        BottomLeftNestedStructure(parse(tokens), parse(tokens))
    elseif next == '⿺'
        TopRightNestedStructure(parse(tokens), parse(tokens))
    elseif next == '⿻'
        OverlapStructure(parse(tokens), parse(tokens))
    else
        Component(next)
    end
end

"""
    IDSFiles.parse(::AbstractString)

Parse a single IDS string that describes a character/component's decomposition. (This function
does not do any input validation to verify that the string can actually represent a
decomposition; in this case, the method will fail.)

This is a relatively low-level method; you typically do not need to call this directly.

# Examples
We can take the IDS string itself and parse it as a nested structure.
```jldoctest
julia> parse("⿱此二")
CharStructure{TopBottomStructure}(TopBottomStructure(Component('此'), Component('二')))
```

"""
function parse(str::AbstractString)
    tokens = Channel{Char}(length(str))
    for char in Iterators.reverse(str)
        put!(tokens, char)
    end

    return parse(tokens)
end


text:
using NURBS.Bases
import LinearAlgebra: dot, cross, norm

struct Surface{BU<:Basis1D, BV<:Basis1D, T<:Real} <: Object2D
    basisU::BU
    basisV::BV
    d1bU::BU
    d2bU::BU
    d3bU::BU
    d1bV::BV
    d2bV::BV
    d3bV::BV
    points::Array{T, 3}

    function Surface(basisU::BU, basisV::BV, points::Array{T, 3}) where {BU<:Basis1D, BV<:Basis1D, T<:Real}
        d1bU = deriv(basisU)
        d2bU = deriv(d1bU)
        d3bU = deriv(d2bU)
        d1bV = deriv(basisV)
        d2bV = deriv(d1bV)
        d3bV = deriv(d2bV)
        new{BU, BV, T}(basisU, basisV, d1bU, d2bU, d3bU, d1bV, d2bV, d3bV, points)
    end

    function Surface(srf::Surface{BU, BV, T}, moves::Array{T, 3}) where {BU<:Basis1D, BV<:Basis1D, T<:Real}
        points = srf.points + moves
        new{BU, BV, T}(
            srf.basisU, srf.basisV, srf.d1bU, srf.d2bU, srf.d3bU, srf.d1bV, srf.d2bV, srf.d3bV, points
        )
    end
end

const BSurface = Surface{BSplineBasis, BSplineBasis}
const NSurface = Surface{NURBSBasis, NURBSBasis}

function evalpoint(srf::Surface, u, v)
    return srf.basisV(v, srf.basisU(u, srf.points))'
end

function evalnormal(srf::Surface, u, v)
    (du, dv, duu, duv, dvv) = evalderiv2(srf, u, v)
    E = dot(du, du)
    F = dot(du, dv)
    G = dot(dv, dv)
    n = cross(du, dv) ./ norm(cross(du, dv))
    # n ./= norm(n)
    L = dot(duu, n)
    M = dot(duv, n)
    N = dot(dvv, n)
    dS2 = E*G - F^2
    K = (L*N - M^2) / dS2
    H = (E*N - 2*F*M + G*L) / (2*dS2)
    Δ = sqrt(H^2 - K)
    if H > 0  # Rhino orders principals by absolute value
        κ1 = H + Δ
        κ2 = H - Δ
    else
        κ1 = H - Δ
        κ2 = H + Δ
    end
    λ1 = (κ1*F - M) / (N - κ1*G)
    λ2 = (κ2*F - M) / (N - κ2*G)
    d1 = (du + λ1*dv) ./ norm(du + λ1*dv)
    # d1 ./= norm(d1)
    d2 = (du + λ2*dv) ./ norm(du + λ2*dv)
    # d2 ./= norm(d2)
    (n, κ1, κ2, du, dv, d1, d2, dS2)
end

function evalderiv1(srf::Surface, u, v)
    du = srf.basisV(v, srf.d1bU(u, srf.points))
    dv = srf.d1bV(v, srf.basisU(u, srf.points))
    (du', dv')
end

function evalderiv2(srf::Surface, u, v)
    duu = srf.basisV(v, srf.d2bU(u, srf.points))
    du = srf.basisV(v, srf.d1bU(u, srf.points))
    duv = srf.d1bV(v, srf.d1bU(u, srf.points))
    dv = srf.d1bV(v, srf.basisU(u, srf.points))
    dvv = srf.d2bV(v, srf.basisU(u, srf.points))
    (du', dv', duu', duv', dvv')
end

function evalderiv3(srf::Surface, u, v)
    duuu = srf.basisV(v, srf.d3bU(u, srf.points))
    duu = srf.basisV(v, srf.d2bU(u, srf.points))
    duuv = srf.d1bV(v, srf.d2bU(u, srf.points))
    duv = srf.d1bV(v, srf.d1bU(u, srf.points))
    du = srf.basisV(v, srf.d1bU(u, srf.points))
    duvv = srf.d2bV(v, srf.d1bU(u, srf.points))
    dv = srf.d1bV(v, srf.basisU(u, srf.points))
    dvv = srf.d2bV(v, srf.basisU(u, srf.points))
    dvvv = srf.d3bV(v, srf.basisU(u, srf.points))
    (du', dv', duu', duv', dvv', duuu', duuv', duvv', dvvv')
end

function BSurface(orderU::Int, orderV::Int, points::Array{T, 3}) where {T<:Real}
    bu = BSplineBasis(range(-1, stop=1, length=size(points, 1)-orderU+2), orderU)
    bv = BSplineBasis(range(-1, stop=1, length=size(points, 2)-orderV+2), orderV)
    Surface(bu, bv, points)
end

function BSurface(KnotsU::Vector{T}, KnotsV::Vector{T}, orderU::Int, orderV::Int, points) where {T<:Real}
    bu = BSplineBasis(2 .* (KnotsU .- minimum(KnotsU)) ./ (maximum(KnotsU) - minimum(KnotsU)) .- 1, orderU)
    bv = BSplineBasis(2 .* (KnotsV .- minimum(KnotsV)) ./ (maximum(KnotsV) - minimum(KnotsV)) .- 1, orderV)
    Surface(bu, bv, points)
end

function isocurve(srf::Surface, dir::Bool, u)
    b = dir ? srf.basisV : srf.basisU
    p = dir ? srf.basisU(u, srf.points) : srf.basisV(u, permutedims(srf.points, [2, 1, 3]))
    return Curve(b, p)
end

NSurface(srf::BSurface) = Surface(NURBSBasis(srf.basisU), NURBSBasis(srf.basisV))
NSurface(points::Array{T, 3}, degU::Int, degV::Int) where {T<:Real} = NSurface(BSurface(points, degU, degV))


text:
using BFloat16s

const XLAScalar = Union{Bool, Int8, Int16, Int32, Int64,
                        UInt8, UInt16, UInt32, UInt64,
                        BFloat16, Float16, Float32, Float64,
                        Complex{Float32}}

struct XlaType
    which::Int32
end

"""
    The runtime representation of an HLO `token` type. This only matters
    for creating happens-before chains in the compiled code. In dynamic
    execution, happens-before is guaranteed by the statement order of
    side effecting instructions.
"""
struct HloToken
    # Tokens don't actually have any data, but we need to prevent
    # the compiler from trying to constant prop it, so we give is
    # some dummy data. If control was returned to the host, control
    # flow order is already clear.
    unused::Int32
end

const xla_type_mapping = Dict{Int32, Type}(
    xla.PrimitiveType.PRED => Bool,
    xla.PrimitiveType.S8   => Int8,
    xla.PrimitiveType.S16  => Int16,
    xla.PrimitiveType.S32  => Int32,
    xla.PrimitiveType.S64  => Int64,
    xla.PrimitiveType.U8   => UInt8,
    xla.PrimitiveType.U16  => UInt16,
    xla.PrimitiveType.U32  => UInt32,
    xla.PrimitiveType.U64  => UInt64,
    xla.PrimitiveType.BF16  => BFloat16,
    xla.PrimitiveType.F16  => Float16,
    xla.PrimitiveType.F32  => Float32,
    xla.PrimitiveType.F64  => Float64,
    xla.PrimitiveType.C64  => Complex{Float32},
    xla.PrimitiveType.TOKEN => HloToken,
)
const reverse_xla_type_mapping = Dict(v=>k for (k,v) in xla_type_mapping)

function Base.convert(::Type{Type}, t::XlaType)
    xla_type_mapping[t.which]
end

function Base.convert(::Type{XlaType}, t::Type)
    XlaType(reverse_xla_type_mapping[t])
end

function assign_data!(lp::LiteralProto, a::Vector{<:XLAScalar})
    if eltype(a) == Bool
        lp.preds = vec(a)
    elseif eltype(a) == Int8
        lp.s8s = vec(a)
    elseif eltype(a) == Int16
        lp.s16s = collect(reinterpret(UInt8, vec(a)))
    elseif eltype(a) == Int32
        lp.s32s = vec(a)
    elseif eltype(a) == Int64
        lp.s64s = vec(a)
    elseif eltype(a) == UInt8
        lp.u8s = vec(a)
    elseif eltype(a) == UInt16
        lp.u16s = collect(reinterpret(UInt8, vec(a)))
    elseif eltype(a) == UInt32
        lp.u32s = vec(a)
    elseif eltype(a) == UInt64
        lp.u64s = vec(a)
    elseif eltype(a) == Float16
        lp.f16s = vec(a)
    elseif eltype(a) == BFloat16
        lp.bf16s = collect(reinterpret(UInt8, vec(a)))
    elseif eltype(a) == Float32
        lp.f32s = vec(a)
    elseif eltype(a) == Float64
        lp.f64s = vec(a)
    elseif eltype(a) == Complex{Float32}
        lp.c64s = collect(reinterpret(Float32, a))
    end
    nothing
end

function Base.convert(::Type{LiteralProto}, a::XLAScalar)
    xlat = convert(XlaType, typeof(a))
    shape = Shape(
        element_type = xlat.which,
        dimensions = Int64[],
        layout = Layout(
            format = Format.DENSE,
        )
    )
    lp = LiteralProto()
    lp.shape = shape
    assign_data!(lp, typeof(a)[a])
    return lp
end

function Base.convert(::Type{LiteralProto}, a::Array{<:XLAScalar})
    # For performance on TPUs, we need to permute the dimensions to
    # have smaller dimensions be more major
    xlat = convert(XlaType, eltype(a))
    perm = sortperm(collect(size(a)); rev=true)
    shape = Shape(xlat, size(a))
    shape.layout.minor_to_major = perm .- 1
    lp = LiteralProto()
    lp.shape = shape
    assign_data!(lp, perm == [] ? vec(a) : vec(permutedims(a, perm)))
    return lp
end

function to_data(lp::LiteralProto)
    elt = lp.shape.element_type
    if elt == xla.PrimitiveType.PRED
        return lp.preds
    elseif elt == xla.PrimitiveType.U8
        return lp.u8s
    elseif elt == xla.PrimitiveType.U16
        return collect(reinterpret(UInt16, lp.u16s))
    elseif elt == xla.PrimitiveType.S16
        return collect(reinterpret(UInt16, lp.s16s))
    elseif elt == xla.PrimitiveType.S32
        return lp.s32s
    elseif elt == xla.PrimitiveType.S64
        return lp.s64s
    elseif elt == xla.PrimitiveType.U32
        return lp.u32s
    elseif elt == xla.PrimitiveType.U64
        return lp.u64s
    elseif elt == xla.PrimitiveType.BF16
        return collect(reinterpret(BFloat16, lp.bf16s))
    elseif elt == xla.PrimitiveType.F16
        return lp.f16s
    elseif elt == xla.PrimitiveType.F32
        return lp.f32s
    elseif elt == xla.PrimitiveType.F64
        return lp.f64s
    elseif elt == xla.PrimitiveType.C64
        return collect(reinterpret(Compplex{Float32}, lp.c64s))
    end
end

function Base.convert(::Type{Array}, proto::LiteralProto)
    s = proto.shape
    if s.layout.minor_to_major != 0:(length(s.dimensions)-1)
        b = to_data(proto)
        perm = tuple((x+1 for x in s.layout.minor_to_major)...)
        r = reshape(b, tuple(collect(s.dimensions)[collect(perm)]...))
        return convert(Array, permutedims(r, Base.invperm(perm)))
    else
        b = to_data(proto)
        return reshape(b, tuple(s.dimensions...))
    end
end


text:
module REPLCompletions

export completions, shell_completions, latex_completions

using Base.Meta

function completes_global(x, name)
    return beginswith(x, name) && !('#' in x)
end

function filtered_mod_names(ffunc::Function, mod::Module, name::String, all::Bool=false, imported::Bool=false)
    ssyms = names(mod, all, imported)
    filter!(ffunc, ssyms)
    syms = UTF8String[string(s) for s in ssyms]
    filter!(x->completes_global(x, name), syms)
end

# REPL Symbol Completions
function complete_symbol(sym, ffunc)
    # Find module
    strs = split(sym, '.')
    # Maybe be smarter in the future
    context_module = Main

    mod = context_module
    lookup_module = true
    t = None
    for name in strs[1:(end-1)]
        s = symbol(name)
        if lookup_module
            # If we're considering A.B.C where B doesn't exist in A, give up
            isdefined(mod, s) || return UTF8String[]
            b = mod.(s)
            if isa(b, Module)
                mod = b
            elseif Base.isstructtype(typeof(b))
                lookup_module = false
                t = typeof(b)
            else
                # A.B.C where B is neither a type nor a
                # module. Will have to be revisited if
                # overloading is allowed
                return UTF8String[]
            end
        else
            # We're now looking for a type
            fields = t.names
            found = false
            for i in 1:length(fields)
                s == fields[i] || continue
                t = t.types[i]
                Base.isstructtype(t) || return UTF8String[]
                found = true
                break
            end
            #Same issue as above, but with types instead of modules
            found || return UTF8String[]
        end
    end

    name = strs[end]

    suggestions = UTF8String[]
    if lookup_module
        # We will exclude the results that the user does not want, as well
        # as excluding Main.Main.Main, etc., because that's most likely not what
        # the user wants
        p = s->(ffunc(mod, s) && s != module_name(mod))
        # Looking for a binding in a module
        if mod == context_module
            # Also look in modules we got through `using`
            mods = ccall(:jl_module_usings, Any, (Any,), Main)
            for m in mods
                append!(suggestions, filtered_mod_names(p, m, name))
            end
            append!(suggestions, filtered_mod_names(p, mod, name, true, true))
        else
            append!(suggestions, filtered_mod_names(p, mod, name, true, false))
        end
    else
        # Looking for a member of a type
        fields = t.names
        for field in fields
            s = string(field)
            if beginswith(s, name)
                push!(suggestions, s)
            end
        end
    end
    suggestions
end

function complete_keyword(s::ByteString)
    const sorted_keywords = [
        "abstract", "baremodule", "begin", "bitstype", "break", "catch", "ccall",
        "const", "continue", "do", "else", "elseif", "end", "export", "finally",
        "for", "function", "global", "if", "immutable", "import", "importall",
        "let", "local", "macro", "module", "quote", "return", "try", "type",
        "typealias", "using", "while"]
    r = searchsorted(sorted_keywords, s)
    i = first(r)
    n = length(sorted_keywords)
    while i <= n && beginswith(sorted_keywords[i],s)
        r = first(r):i
        i += 1
    end
    sorted_keywords[r]
end

function complete_path(path::ByteString)
    matches = ByteString[]
    dir, prefix = splitdir(path)
    if length(dir) == 0
        files = readdir()
    elseif isdir(dir)
        files = readdir(dir)
    else
        return matches
    end
    for file in files
        if beginswith(file, prefix)
            p = joinpath(dir, file)
            push!(matches, isdir(p) ? joinpath(p,"") : p)
        end
    end
    matches
end

function complete_methods(input::String)
    tokens = split(input, '.')
    fn = Main
    for token in tokens
        sym = symbol(token)
        isdefined(fn, sym) || return UTF8String[]
        fn = fn.(sym)
    end
    isgeneric(fn) || return UTF8String[]
    UTF8String[string(m) for m in methods(fn)]
end

include("latex_symbols.jl")

const non_identifier_chars = [" \t\n\r\"\\'`\$><=:;|&{}()[],+-*/?%^~"...]
const non_filename_chars = [" \t\n\r\"\\'`@\$><=;|&{("...]
const whitespace_chars = [" \t\n\r"...]

# Aux function to detect whether we're right after a
# using or import keyword
function afterusing(string::ByteString, startpos::Int)
    (isempty(string) || startpos == 0) && return false
    str = string[1:prevind(string,startpos)]
    isempty(str) && return false
    rstr = reverse(str)
    r = search(rstr, r"\s(gnisu|tropmi)\b")
    isempty(r) && return false
    fr = chr2ind(str, length(str)-ind2chr(rstr, last(r))+1)
    return ismatch(r"^\b(using|import)\s*(\w+\s*,\s*)*\w*$", str[fr:end])
end

function latex_completions(string, pos)
    slashpos = rsearch(string, '\\', pos)
    if rsearch(string, whitespace_chars, pos) < slashpos
        # latex symbol substitution
        s = string[slashpos:pos]
        latex = get(latex_symbols, s, "")
        if !isempty(latex) # complete an exact match
            return (true, ([latex], slashpos:pos, true))
        else
            # return possible matches; these cannot be mixed with regular
            # Julian completions as only latex symbols contain the leading \
            latex_names = filter(k -> beginswith(k, s), keys(latex_symbols))
            return (true, (sort!(collect(latex_names)), slashpos:pos, true))
        end
    end
    return (false, (UTF8String[], 0:-1, false))
end

function completions(string, pos)
    inc_tag = Base.incomplete_tag(parse(string[1:pos], raise=false))
    if inc_tag in [:cmd, :string]
        startpos = nextind(string, rsearch(string, non_filename_chars, pos))
        r = startpos:pos
        paths = complete_path(string[r])
        if inc_tag == :string && length(paths) == 1 && !isdir(paths[1])
            paths[1] *= "\""
        end
        return sort(paths), r, true
    end

    ok, ret = latex_completions(string,pos)
    ok && return ret

    if inc_tag == :other && string[pos] == '('
        endpos = prevind(string, pos)
        startpos = nextind(string, rsearch(string, non_identifier_chars, endpos))
        return complete_methods(string[startpos:endpos]), startpos:endpos, false
    elseif inc_tag == :comment
        return UTF8String[], 0:-1, false
    end

    dotpos = rsearch(string, '.', pos)
    startpos = nextind(string, rsearch(string, non_identifier_chars, pos))

    ffunc = (mod,x)->true
    suggestions = UTF8String[]
    comp_keywords = true
    if afterusing(string, startpos)
        # We're right after using or import. Let's look only for packages
        # and modules we can reach from here

        # If there's no dot, we're in toplevel, so we should
        # also search for packages
        s = string[startpos:pos]
        if dotpos <= startpos
            append!(suggestions, filter(readdir(Pkg.dir())) do pname
                pname[1] != '.' &&
                pname != "METADATA" &&
                pname != "REQUIRE" &&
                beginswith(pname, s)
            end)
        end
        ffunc = (mod,x)->(isdefined(mod, x) && isa(mod.(x), Module))
        comp_keywords = false
    end
    startpos == 0 && (pos = -1)
    dotpos <= startpos && (dotpos = startpos - 1)
    s = string[startpos:pos]
    comp_keywords && append!(suggestions, complete_keyword(s))
    append!(suggestions, complete_symbol(s, ffunc))
    return sort(unique(suggestions)), (dotpos+1):pos, true
end

function shell_completions(string,pos)
    # First parse everything up to the current position
    scs = string[1:pos]
    local args, last_parse
    try
        args, last_parse = Base.shell_parse(scs, true)
    catch
        return UTF8String[], 0:-1, false
    end
    # Now look at the last this we parsed
    isempty(args.args[end].args) && return UTF8String[], 0:-1, false
    arg = args.args[end].args[end]
    if isa(arg,String)
        # Treat this as a path (perhaps give a list of comands in the future as well?)
        dir,name = splitdir(arg)
        if isempty(dir)
            files = readdir()
        else
            isdir(dir) || return UTF8String[], 0:-1, false
            files = readdir(dir)
        end
        # Filter out files and directories that do not begin with the partial name we were
        # completing and append "/" to directories to simplify further completion
        ret = map(filter(x->beginswith(x, name), files)) do x
            if !isdir(joinpath(dir, x))
                return x
            else
                return x*"/"
            end
        end
        r = (nextind(string, pos-sizeof(name))):pos
        return ret, r, true
    elseif isexpr(arg, :escape) && (isexpr(arg.args[1], :incomplete) || isexpr(arg.args[1], :error))
        r = first(last_parse):prevind(last_parse, last(last_parse))
        partial = scs[r]
        ret, range = completions(partial, endof(partial))
        range += first(r) - 1
        return ret, range, true
    end
end

end # module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
export factor_graph, rank_reveal, projectors, split_into_clusters, decode_factor_graph_state

function split_into_clusters(ig::LabelledGraph{S, T}, assignment_rule) where {S, T}
    cluster_id_to_verts = Dict(
        i => T[] for i in values(assignment_rule)
    )

    for v in vertices(ig)
        push!(cluster_id_to_verts[assignment_rule[v]], v)
    end

    Dict(
        i => first(cluster(ig, verts)) for (i, verts) ∈ cluster_id_to_verts
    )
end

function factor_graph(
    ig::IsingGraph,
    num_states_cl::Int;
    spectrum::Function=full_spectrum,
    cluster_assignment_rule::Dict{Int, T} # e.g. square lattice
) where {T}
    ns = Dict(i => num_states_cl for i ∈ Set(values(cluster_assignment_rule)))
    factor_graph(
        ig,
        ns,
        spectrum=spectrum,
        cluster_assignment_rule=cluster_assignment_rule
    )
end

function factor_graph(
    ig::IsingGraph,
    num_states_cl::Dict{T, Int};
    spectrum::Function=full_spectrum,
    cluster_assignment_rule::Dict{Int, T} # e.g. square lattice
) where {T}
    L = maximum(values(cluster_assignment_rule))
    fg = LabelledGraph{MetaDiGraph}(sort(unique(values(cluster_assignment_rule))))

    for (v, cl) ∈ split_into_clusters(ig, cluster_assignment_rule)
        sp = spectrum(cl, num_states=get(num_states_cl, v, basis_size(cl)))
        set_props!(fg, v, Dict(:cluster => cl, :spectrum => sp))
    end

    for (i, v) ∈ enumerate(vertices(fg)), w ∈ vertices(fg)[i+1:end]
        cl1, cl2 = get_prop(fg, v, :cluster), get_prop(fg, w, :cluster)

        outer_edges, J = inter_cluster_edges(ig, cl1, cl2)

        if !isempty(outer_edges)
            en = inter_cluster_energy(
                get_prop(fg, v, :spectrum).states, J, get_prop(fg, w, :spectrum).states
            )
            pl, en = rank_reveal(en, :PE)
            en, pr = rank_reveal(en, :EP)
            add_edge!(fg, v, w)
            set_props!(
                fg, v, w, Dict(:outer_edges => outer_edges, :pl => pl, :en => en, :pr => pr)
            )
        end
    end
    fg
end

function factor_graph(
    ig::IsingGraph;
    spectrum::Function=full_spectrum,
    cluster_assignment_rule::Dict{Int, T}
) where {T}
    factor_graph(ig, Dict{T, Int}(), spectrum=spectrum, cluster_assignment_rule=cluster_assignment_rule)
end

function rank_reveal(energy, order=:PE)
    @assert order ∈ (:PE, :EP)
    dim = order == :PE ? 1 : 2

    E, idx = unique_dims(energy, dim)

    if order == :PE
        P = zeros(size(energy, 1), size(E, 1))
    else
        P = zeros(size(E, 2), size(energy, 2))
    end

    for (i, elements) ∈ enumerate(eachslice(P, dims=dim))
        elements[idx[i]] = 1
    end

    order == :PE ? (P, E) : (E, P)
end

function decode_factor_graph_state(fg, state::Vector{Int})
    ret = Dict{Int, Int}()
    for (i, vert) ∈ zip(state, vertices(fg))
        spins = get_prop(fg, vert, :cluster).labels
        states = get_prop(fg, vert, :spectrum).states
        if length(states) > 0
            curr_state = states[i]
            merge!(ret, Dict(k => v for (k, v) ∈ zip(spins, curr_state)))
        end
    end
    ret
end