text:
using DynamicGrids, DynamicGridsInteract, Test, Colors, ColorSchemes, ImageMagick

# life glider sims

init =  Bool[0 0 0 0 0 0
             0 0 0 0 0 0
             0 0 0 0 0 0
             0 0 0 1 1 1
             0 0 0 0 0 1
             0 0 0 0 1 0]

test3 = [0 0 0 0 0 0
         0 0 0 0 0 0
         0 0 0 0 1 1
         0 0 0 1 0 1
         0 0 0 0 0 1
         0 0 0 0 0 0]

test5 = [0 0 0 0 0 0
         0 0 0 0 0 0
         1 0 0 0 1 1
         1 0 0 0 0 0
         0 0 0 0 0 1
         0 0 0 0 0 0]

processor = ColorProcessor(
    scheme=ColorSchemes.leonardo, zerocolor=nothing, maskcolor=nothing, textconfig=nothing,
)

@testset "InteractOutput" begin
    l0 = ARGB32(get(ColorSchemes.leonardo, 0))
    l1 = ARGB32(get(ColorSchemes.leonardo, 1))

    leonardo2 = [l0 l0 l0 l0 l0 l0
                 l0 l0 l0 l0 l0 l0
                 l1 l0 l0 l0 l1 l1
                 l1 l0 l0 l0 l0 l0
                 l0 l0 l0 l0 l0 l1
                 l0 l0 l0 l0 l0 l0]

    ruleset = Ruleset(Life(); boundary=Wrap())
    output = InteractOutput(init; tspan=1:2, ruleset=ruleset, store=true, processor=processor);
    sim!(output, ruleset)
    sleep(10)
    resume!(output, ruleset; tstop=5)
    sleep(2)

    @test output[1] == init
    @test output[3] == test3
    @test output[5] == test5

    @testset "output image matches colorscheme" begin
        @test output.image_obs[].children.tail[1] == leonardo2
    end

    @testset "output works with store=false" begin
        output = InteractOutput(init; ruleset=ruleset, tspan=1:3, store=false, processor=processor);
        sim!(output, ruleset)
        output.graphicconfig.stoppedframe
        DynamicGrids.stoppedframe(output)
        sleep(10)
        resume!(output, ruleset; tstop=5)
        sleep(2)
        @test output[1] == test5
    end

end

if !Sys.islinux() # No graphic head loaded in CI: TODO add this
    @testset "ElectronOutput" begin
        ruleset = Ruleset(Life(); boundary=Wrap())
        output = ElectronOutput(init; ruleset=ruleset, tspan=1:3, store=true, processor=processor)
        DynamicGrids.setrunning!(output, false)
        sim!(output.interface, ruleset)
        sleep(10)
        DynamicGrids.setrunning!(output, false)
        resume!(output.interface, ruleset; tstop=5)
        sleep(2)
        @test output[3] == test3
        @test output[5] == test5
        close(output.window)
    end
end

@testset "ServerOutput" begin
    ruleset = Ruleset(Life(); boundary=Wrap())
    ServerOutput(init; ruleset=ruleset, port=8080, processor=processor)
    # TODO: test the server somehow
end


text:
using StaticArrays, LinearAlgebra, Plots

# function f(u,t)
#   T = eltype(u)
#   @SVector [u[2], 10*(1-u[1]^2)*u[2] - u[1]]
# end

function f(u,t)
  T = eltype(u)
  @SVector [one(T) + u[1]^2*u[2] - 4*u[1],
            3u[1]-u[1]^2*u[2]]
end

function norm_sc(ys::AbstractVector{T}; sc::AbstractVector{T}) where T
  @assert length(ys) == length(sc)

  maxabs = zero(T)
  for i = 1:length(ys)
    maxabs = max(maxabs,abs(ys[i])/sc[i])
  end
  maxabs
end

function run()
  T = Float64
  A = @SMatrix T[0            0           0          0          0         0     0;
                 1//5         0           0          0          0         0     0;
                 3//40        9//40       0          0          0         0     0;
                 44//45      -56//15      32//9       0          0         0     0;
                 19372//6561 -25360//2187 64448//6561 -212//729     0         0     0;
                 9017//3168   -355//33   46732//5247   49//176 -5103//18656  0     0;
                 35//384       0         500//1113  125//192 -2187//6784  11//84 0]
  b = @SMatrix T[35//384 0 500//1113 125//192 -2187//6784 11//84 0;
                 5179//57600 0 7571//16695 393//640 -92097//339200 187//2100 1//40]
  c = @SVector T[0, 1//5, 3//10, 4//5, 8//9, 1, 1]
  e = b[1,:] - b[2,:]
  p = 5

  y0 = @SVector [T(1.01), T(3)]
  t0 = zero(T)
  e1 = e0 = zero(T)
  timeend = T(20)

  Atol = 1e-6
  Rtol = 1e-5
  sc = Atol.+abs.(y0).*Rtol

  d0 = norm_sc(y0; sc=sc)
  d1 = norm_sc(f(y0, t0); sc=sc)

  h0 = (d0 >= 1e-5 && d1 >= 1e-5) ? 0.01*(d0/d1) : 1e-6
  y1 = y0 + h0 * f(y0, t0)
  t1 = t0 + h0

  d2 = norm_sc(f(y1, t1)-f(y0, t0); sc=sc) / h0
  h1 = (max(d1, d2) <= 1e-15) ?  max(1e-6, h0 * 1e-3) : (0.01 / max(d1, d2))^1/(p+1)

  @show dt = min(100*h0, h1)

  nsteps = ceil(Int64, timeend / dt)
  dt = timeend / nsteps

  ts = Vector{typeof(t0)}(undef,nsteps+1)
  ys = Vector{typeof(y0)}(undef,nsteps+1)
  es = Vector{typeof(e0)}(undef,nsteps+1)

  ts[1] = t0
  ys[1] = copy(y0)
  es[1] = e0

  for n = 1:nsteps

    acceptstep = false

    fac = T(0.9)
    facmin = T(0.1)
    facmax = T(5)

    while !acceptstep
      y1 = y0
      t1 = t0
      K1 = f(y1, t1)

      y1 = y0 + dt*A[2,1]*K1
      t1 = t0 + c[2]*dt
      K2 = f(y1, t1)

      y1 = y0 + dt*(A[3,1]*K1 + A[3,2]*K2)
      t1 = t0 + c[3]*dt
      K3 = f(y1, t1)

      y1 = y0 + dt*(A[4,1]*K1 + A[4,2]*K2 + A[4,3]*K3)
      t1 = t0 + c[4]*dt
      K4 = f(y1, t1)

      y1 = y0 + dt*(A[5,1]*K1 + A[5,2]*K2 + A[5,3]*K3 + A[5,4]*K4)
      t1 = t0 + c[5]*dt
      K5 = f(y1, t1)

      y1 = y0 + dt*(A[6,1]*K1 + A[6,2]*K2 + A[6,3]*K3 + A[6,4]*K4 + A[6,5]*K5)
      t1 = t0 + c[6]*dt
      K6 = f(y1, t1)

      y1 = y0 + dt*(A[7,1]*K1 + A[7,2]*K2 + A[7,3]*K3 + A[7,4]*K4 + A[7,5]*K5 + A[7,6]*K6)
      t1 = t0 + c[7]*dt
      K7 = f(y1, t1)

      ye = dt*(e[1]*K1 + e[2]*K2 + e[3]*K3 + e[4]*K4 + e[5]*K5 + e[6]*K6 + e[7]*K7)

      t1 = t0 + dt
      # y1 is computed in the last stage

      k1 = one(T)
      sc = Atol.+abs.(y1).*Rtol
      e1 = norm_sc(ye, sc=sc)

      beta = e1^(-k1/p)
      dt *= min(facmax, max(facmin, fac*beta))

      if e1 <= 1
        acceptstep = true
      else
        facmax = fac = T(0.9)
      end
    end

    ys[n+1] = y1
    ts[n+1] = t1
    es[n+1] = e1

    t0 = t1
    y0 = y1
    e0 = e1
  end

  (ys, ts, es)
end

(ys, ts, es) = run()

a = [y[1] for y in ys]
b = [y[2] for y in ys]

plot(ts, [a, b])

plot(ts[2:end], es[2:end])
yaxis!("error", :log10)

ext: jl

max_stars_repo_name: lcw/adaptiverk

max_issues_repo_name: lcw/adaptiverk

max_forks_repo_name: lcw/adaptiverk


text:
@testset "Morphism" begin
  L = Zlattice(gram = zero_matrix(ZZ, 0, 0))
  @test_throws ArgumentError shortest_vectors(L)
  @test_throws ArgumentError shortest_vectors(L, Vector{Int})
  @test_throws ArgumentError minimum(L)
  @test_throws ArgumentError kissing_number(L)
  @test (@inferred short_vectors(L, 1)) == []
  @test (@inferred short_vectors(L, 1, 2)) == []
end


text:
# We consider as an illustrative example a simple harmonic oscillator system.
# The set propagation definitions introduced so far are illustrated using a single
# degree of freedom second order problem given by the following equation:

# ```math
# u''(t) + ω^2 u(t) = 0,
# ```
# where $ω$ is a scalar and $u(t)$ is the unknown.
#
# This problem can be associated with a spring-mass system, where $u(t)$ is the
# elongation of the spring at time $t$ and the mass and stiffness set a natural
# frequency $ω$. In this case we consider $ω = 4π$.

# Let us introduce the state variable $v(t) := u'(t)$ and define the vector
# $x(t) = [u(t), v(t)]^T$. Then we can rewrite the system in the following first order form
#
#
# ```math
# x'(t) = \begin{pmatrix} 0 & 1 \\ -\omega^2 & 0 \end{pmatrix} x(t)
#
# Two cases are considered for the initial conditions, first a single initial condition
# given by $u(0)=1$ and $v(0)=0$, and a second case where the initial conditions belong to a set.

# This problem is solved using set propagation and the results obtained for the
# single initial condition are shown in the following figures.
#
# The time step-size used is $\delta = 0.025$.
#
# The set $\Omega_0$ includes the analytic trajectory within $[0, \delta]$, and
# such set is propagated to cover the solution for further time intervals as shown below.
#
# It is worth noting that even when the initial condition is a singleton,
# the set propagation technique returns a sequence of sets, and such sequence guarantees
# an enclosure of the analytic solution at \emph{any} time $t \in [0, T] \subset \mathbb{R}$.
#
# For comparison, we also plot the analytic solution (magenta solid line) and a
# numerical solution obtained using the Newmark's method with time step $\delta/5$
# (red triangle markers).

# The result obtained when a set of initial conditions, together with few dozen
# trajectories with random initial conditions drawn from the initial set
# $\mathcal{X}_0 = \mathcal{U}_0 \times \mathcal{V}_0 = [0.9, 1.1] \times [-0.1, 0.1]$,
# where $\times$ denotes the Cartesian product.

# In both cases, it is shown that a single set-based integration covers infinitely
# many trajectories \emph{in dense time}, i.e. for all intermediate times between
# $0$ and $T$, where $T > 0$ is the time horizon, there is a reach-set that covers
# all the exact solutions.

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
function habitable_exoplanet(L, r)
    habitability = ""

    # Your code goes here!

    return habitability
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
function setupTravelTimeTomography(m,filenamePrefix::String, resultsOutputFolderAndPrefix::String,useFilesForFields::Bool = false)

file = matopen(string(filenamePrefix,"_PARAM.mat"));
n_cells = read(file,"MinvN");
OmegaDomain = read(file,"MinvOmega");
Minv = getRegularMesh(OmegaDomain,n_cells);
HO = read(file,"HO");

boundsLow = read(file,"boundsLow");
boundsHigh = read(file,"boundsHigh");
mref =  read(file,"mref");
close(file);

resultsFilename = "";
if resultsOutputFolderAndPrefix!=""
	resultsFilename = string(resultsOutputFolderAndPrefix,tuple((Minv.n+1)...),".dat");
end
###########################################################

### Read receivers and sources files
RCVfile = string(filenamePrefix,"_rcvMap.dat");
SRCfile = string(filenamePrefix,"_srcMap.dat");

srcNodeMap = readSrcRcvLocationFile(SRCfile,Minv);
rcvNodeMap = readSrcRcvLocationFile(RCVfile,Minv);

Q = generateSrcRcvProjOperators(Minv.n .+ 1,srcNodeMap);
Q = Q.*1/(norm(Minv.h)^2);
P = generateSrcRcvProjOperators(Minv.n .+ 1,rcvNodeMap);

println("Travel time tomography: ",size(Q,2)," sources.");
#############################################################################################

println("Reading data:");

(DobsEik,WdEik) =  readDataFileToDataMat(string(filenamePrefix,"_travelTime.dat"),srcNodeMap,rcvNodeMap);

N = prod(Minv.n .+ 1);

Iact = SparseMatrixCSC{Float16,Int32}(I, N, N);
#Iact = convert(SparseMatrixCSC{Float16,Int32},Iact);
mback   = zeros(Float64,N);

## Setting the sea constant:
mask = zeros(N);
sea = abs.(m[:] .- minimum(m)) .< 5e-2;
mask[sea] .= 1.0;
# # setup active cells
mback = vec(m[:].*mask);
sback = velocityToSlowSquared(mback)[1];
sback[mask .== 0.0] .= 0.0;

Iact = Iact[:,mask .== 0.0];
boundsLow = Iact'*boundsLow;
boundsHigh = Iact'*boundsHigh;
mref = Iact'*mref[:];

########################################################################################################
##### Set up remote workers ############################################################################
########################################################################################################

EikMPIWorkers = nworkers(); # this just set the maximal MPI workers. To activate parallelism, run addprocs()

(pFor,contDiv,SourcesSubInd) = getEikonalInvParam(Minv,Q,P,HO,EikMPIWorkers,useFilesForFields);

misfun = SSDFun

Wd 		= Array{Array{Float64}}(undef, length(pFor));
dobs 	= Array{Array{Float64}}(undef, length(pFor));
for i=1:length(pFor)
	I_i = SourcesSubInd[i];
	Wd[i]   = WdEik[:,I_i];
	dobs[i] = DobsEik[:,I_i];
end
WdEik = 0;

pMisRFs = getMisfitParam(pFor, Wd, dobs, misfun, Iact,sback);

return (Q,P,pMisRFs,SourcesSubInd,contDiv,Iact,sback,mref,boundsHigh,boundsLow,resultsFilename,DobsEik);
end


text:
@testset "operations.jl" begin

end


text:
module Commands

include("../types.jl")
include("../controller/requests.jl")
include("network.jl")

using .Types, Bcrypt, Sockets

include("../commands/create.jl")
include("../commands/help.jl")
include("../commands/join.jl")
include("../commands/leave.jl")
include("../commands/login.jl")
include("../commands/register.jl")
include("../commands/who.jl")
include("../commands/unknown.jl")

function exec_command(command::Vector{SubString{String}}, storage, conn::IO)
    ip_addr = string(first(getpeername(conn)))
    command_name = String(command[1])
    args = command[2:length(command)]
    commandInfos = getCommandInfos(command_name)

    if (Network.is_logged(storage, ip_addr) == false && commandInfos.auth == true)
        return (write(conn, "You must be logged in to use this command\n"))
    end
    if (size(args)[1] < commandInfos.args_amount)
        return (write(conn, "Invalid number of arguments\n"))
    end

    call(commandInfos, args, storage, conn)
end

end # Commands


text:
using O_structs
using B_structs


text:
"""
   Linked List

# Brief:
    A linked list is a data structure where each element is connected to the next one.

# Complexity of some operations
    - insert            - O(N)
    - insert to front   - O(1)
    - delete            - O(N)
    - delete first      - O(1)
    - get element       - O(N)

# Functions
    - create_node(val, next=missing) - Create node with value 'val' and the pointer to the next node ( missing by default )
    - create_list(n::Int=0, val=missing) - Create root node of the list with n elements with value set to 'val'
    - insert(list::Node, new_node::Node, index::Int=1) - Add a new node to the list at the specified index
    - insert(list::Node, val, index::Int=1) - Create a new node with value 'val' in the list at the specified index
    - push_back(list::Node, node::Node) - Add a new node to the end of the list
    - push_back(list::Node, val) - Add a new node with value 'val' to the end of the list
    - return_as_array(list::Node) - Return the array representation of the list
    - clear(list::Node) - Remove all elements from the list
    - remove(list::Node, index::Int) - Remove an element at the specified index
    - remove_all(list::Node, val) Remove all elements with the value 'val'
    - remove_first(list::Node, val) - Remove the first element in the list with value 'val'
    - get_node(list::Node, index::Int) - Get a node at the specified index
    - get(list::Node, index::Int) - Get a value from the node at the specified index
    - indexOf(list::Node, val) - Return the index of the first element with a value 'val'
    - is_empty(list) - Return true if list is empty, false if it has elements

# Contributed by: [Nikola Mircic](https://github.com/Nikola-Mircic)
"""
module LinkedList

mutable struct Node
    val::Any # Value of the node
    next::Any # Pointer to the next node
    Node(val, next = missing) = new(val, next)
end

# Create node with value 'val' and the pointer to the next node ( missing by default )
function create_node(val, next = missing)
    return Node(val, next)
end

# Create root node of the list
function create_list(n::Int = 0, val = missing)
    root_node = create_node("root") # Root node with value "root"

    current_node = root_node
    for i in 1:n
        new_node = create_node(val, missing)
        current_node.next = new_node

        current_node = current_node.next
    end

    return root_node
end

# Add a new node to the list at the specified index
function insert(list::Node, new_node::Node, index::Int = 1)
    current_node = get_node(list, index - 1)

    next = current_node.next
    current_node.next = new_node
    return new_node.next = next
end

# Create a new node with value 'val' in the list at the specified index
function insert(list::Node, val, index::Int = 1)
    new_node = create_node(val)
    return insert(list, new_node, index)
end

# Add a new node to the end of the list
function push_back(list::Node, node::Node)
    current_node = list
    while !ismissing(current_node.next)
        current_node = current_node.next
    end

    return current_node.next = node
end

# Add a new node with value 'val' to the end of the list
function push_back(list::Node, val)
    new_node = create_node(val)
    return push_back(list, new_node)
end

function return_as_array(list::Node)
    arr = Array{Any}(missing, 0)
    node = list.next

    while !ismissing(node)
        push!(arr, node.val)
        node = node.next
    end

    return arr
end

# Remove all elements from the list
function clear(list::Node)
    return list.next = missing
end

# Remove an element at the specified index
function remove(list::Node, index::Int)
    node_before = get_node(list, index - 1)
    to_remove = node_before.next
    return node_before.next = to_remove.next
end

# Remove all elements with the value 'val'
function remove_all(list::Node, val)
    current_node = list
    while !ismissing(current_node.next)
        next_node = current_node.next
        if next_node.val == val
            current_node.next = next_node.next
        else
            current_node = next_node
        end
    end
end

# Remove the first element in the list with value 'val'
function remove_first(list::Node, val)
    index = indexOf(list, val)
    return remove(list, index)
end

# Get a node at the specified index
function get_node(list::Node, index::Int)
    current_node = list
    for i in 1:index
        if !ismissing(current_node)
            current_node = current_node.next
        else
            return missing
        end
    end
    return current_node
end

# Get a value from the node at the specified index
function get(list::Node, index::Int)
    node = get_node(list, index)
    if !ismissing(node)
        return node.val
    end
    return missing
end

# Return the index of the first element with a value 'val'
function indexOf(list::Node, val)
    current_node = list.next
    i = 1
    while !ismissing(current_node)
        if current_node.val == val
            return i
        else
            current_node = current_node.next
            i += 1
        end
    end
    return -1
end

# Return true if list is empty, false if it has elements
function is_empty(list)
    return ismissing(list.next)
end

end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
using Distributed
using ClearStacktrace

#using UnicodePlots
#using Plots
#unicodeplots()

#if nprocs()==1
#	addprocs(8)
#end
#@everywhere include("spike_distance_opt.jl")
include("spike_distance_opt.jl")

P, C = make_net_SNN(26)

E, I = P #, EEA]
EE, EI, IE, II = C
SNN.monitor([E, I], [:fire])
#global E_stim = []#Vector
sim_length = 500
@inbounds for t = 1:sim_length*ms
    E.I = vec([11.5 for i = 1:sim_length*ms])#vec(E_stim[t,:])#[i]#3randn(Ne)
    SNN.sim!(P, C, 1ms)

end
#_,_,_,spkd_ground = raster_synchp(P[1])
spkd_ground = get_trains(P[1])
sgg = [convert(Array{Float32,1}, sg) for sg in spkd_ground]

#MU = 10
#ɛ = MU / 2#0.125
#parallelization = :multi,
#f(x) = 10length(x) + sum( x.^2 - 10cos.(2π*x)  )

#Instantiate the bounds, note that bounds should be a $2\times 10$ Matrix where the first row corresponds to the lower bounds whilst the second row corresponds to the upper bounds.

#D = 10
#nobjectives = length(pf[1].f)
#npartitions = nobjectives == 2 ? 100 : 12

#methods = [
#        SMS_EMOA(N = 50, n_samples=500, options=options),
#        NSGA2(options=options),
#        MOEAD_DE(gen_ref_dirs(nobjectives, npartitions), options=Options( seed = 1, iterations = 500)),
#        NSGA3(options=options),
#      ]

#for method in methods

#f_calls = 0
f(x) = begin
    f_calls += 1
    loss(x)
end
f_calls = 0
#E1, spkd_found = eval_best(params)
#errors = loss(10)
#D = length(loss(10))
#@show(errors)
#bounds = Matrix([3.0; 40.0])
#bounds = [ 3.0; 40.0 ]'
D = 10
bounds = [3ones(D) 40ones(D)]'
a = view(bounds, 1, 1)
b = view(bounds, 1, 2)
#b = problem.bounds[2]
#@show(b)
#@show(a)
#@show(method)
#options = Options( seed = 1, iterations=10, f_calls_limit =10)
#N = 50, n_samples=500,
#method = NSGA2(options=options, f_calls_limit =100)#, information = information)
@show(bounds)
@show(loss)

information = Information(f_optimum = 0.0)
method = ECA(options=options, information = information)
result = optimize(loss, bounds, method)
@show(result)
#Approximate the optimum using the function optimize.

#result = optimize(loss, bounds)
#=
options = GA(
    populationSize = MU,
    ɛ = 4,
    mutationRate = 0.5,
    selection = ranklinear(1.5),#ranklinear(1.5),#ss,
    crossover = intermediate(0.5),#xovr,
    mutation = uniform(0.5),#(.015),#domainrange(fill(1.0,ts)),#ms
)

#Random.seed!(0);
result = Evolutionary.optimize(
    loss,
    lower,
    upper,
    initd,
    options,
    Evolutionary.Options(
        iterations = 50,
        successive_f_tol = 1,
        show_trace = true,
        store_trace = true,
    ),
)
fitness = minimum(result)
=#
#parallelization = :thread,

filename = string("GAsolution.jld")#, py"target_num_spikes")#,py"specimen_id)
params = result.minimizer
E1, spkd_found = eval_best(params)
save(
    filename,
    "spkd_ground",
    spkd_ground,
    "spkd_found",
    spkd_found,
    "Ne",
    Ne,
    "Ni",
    Ni,
    "sim_length",
    sim_length,
)
println("best result")
loss(result.minimizer)
println("σee = 0.5,  pee= 0.8,σei = 0.5,  pei= 0.8")

#println("σee = 0.45,  pee= 0.8,σei = 0.4,  pei= 0.9)")
@show(result.minimizer)

@show(fitness)

@show(result)
@show(result.trace)
trace = result.trace
dir(x) = fieldnames(typeof(x))
dir(trace[1, 1, 1])
trace[1, 1, 1].metadata#["population"]
filename = string("PopulationScatter.jld")#, py"target_num_spikes")#,py"specimen_id)
save(filename, "trace", trace)
#evo_population = [t.metadata[""] for t in trace]
E1, spkd_found = eval_best(params)

evo_loss = [t.value for t in trace]

display(plot(evo_loss))

#first_dim1 = [t.metadata["population"][1][1] for t in trace]
#first_dim2 = [t.metadata["population"][1][2] for t in trace]
#first_dim3 = [t.metadata["population"][1][3] for t in trace]
#first_dim4 = [t.metadata["population"][1][4] for t in trace]

#display(plot(first_dim1))
#display(plot(first_dim1,first_dim2,first_dim3))

#run(`python-jl validate_candidate.py`)
using PyCall
py"""import validate_candidate"""
iter = [t.iteration for t in trace]
data = [ trace[i+1,1,1].metadata["pop"] for i in iter ]

evo_loss

model = Chain(Dense(d, 15, relu), Dense(15, nclasses))

@info "MLP" loss=loss(data, evomodel) accuracy = accuracy(data, evomodel)
#model = Dense(2, 1, σ)
#L(x,y) = Flux.mse(model(x), y)
#opt = SGD(params(model))
#Flux.train!(L, zip(xs, ys), opt)
function loss(model)
    σee = model[1]
    pee = model[2]
    σei = model[3]
    pei = model[4]
    P1, C1 = make_net(Ne, Ni, σee = σee, pee = pee, σei = σei, pei = pei)#,a=a)
    E1, I1 = P1
    SNN.monitor([E1, I1], [:fire])
    sim_length = 500
    @inbounds for t = 1:sim_length*ms
        E1.I = vec([11.5 for i = 1:sim_length])#vec(E_stim[t,:])#[i]#3randn(Ne)
        SNN.sim!(P1, C1, 1ms)
    end

    spkd_found = get_trains(P1[1])
    println("Ground Truth \n")
    SNN.raster([E]) |> display
    println("Best Candidate \n")

    SNN.raster([E1]) |> display

    error = raster_difference(spkd_ground, spkd_found)
    error
end
#loss(model) = (x,y)->logitcrossentropy(model(x), y)
#loss(model,x,y) = loss(model)(x, y)
#loss(xy, model) = loss(model)(hcat(map(first,xy)...), hcat(map(last,xy)...))

opt = ADAM(1e-4)
evalcb = Flux.throttle(() -> @show(loss(data, model), accuracy(data, model)), 5)
for i in 1:500
    Flux.train!(loss(model), params(model), data, opt, cb = evalcb)
end
# ## Visualize the result

#contour(0:.1:1, 0:.1:1, (x, y) -> model([x,y])[].data, fill=true)


text:
#+ =========== parse ============

@inline to_words(l) = split(l,r"[^\S\n\r]",keepempty=false)

function parse_file!(
    shift::Int,
    t::Block{TR},
    lines::Vector{String},
    lines1::Vector{String},
    codeinv
    )::Block{TR} where TR

    @assert is_valid_x(t)
    if is_multi(t)
        Tx1 = concat0(t.C[1].x,t.C[end].x)
        N1  = length(Tx1)
        #@assert N1*t.n == length(t.x)
        #: parse t.n times !!
        for k=1:t.n
            shift_n = shift+(k-1)*N1  #: the accumulated shift
            for i=1:length(t.C)
                t.C[i] = parse_file!(shift_n, t.C[i], lines, lines1, codeinv)
            end
        end
    else
        TX = (first(t.x)+shift):(last(t.x)+shift)
        words1 = to_words.(lines1[TX])
        # extract data
        data = TokenValueLine[]
        for (line,word) ∈ zip(lines[TX],words1)
            ii = 0
            keywords = Pair{String,Any}[w=>((cmp(w,"£")==0 || w ∈ MASK_RULES_DIC_INV_KEYS_STRIPPED_NO_£)
                                            ? MASK_RULES_DIC_INV1[w]
                                            : nothing)
                                        for w ∈ word]
            for (i,(w,kw)) ∈ enumerate(keywords)
                line_ii = line[ii+1:end]
                if kw!==nothing
                    r = findfirst(kw,line_ii) #: note the white space !!!
                    if r!==nothing
                        #push!(dt, w=>string(line_ii[r]))
                        keywords[i] = w=>string(line_ii[r])
                    else
                        throw(error("Incompatible keyword structure in line $(line) at word $(w) and keyword $(kw)."))
                    end
                    ii += last(r)
                end
            end
            push!(data, keywords)
        end
        push!(t.DATA, TX=>data)
    end
    return t
end

# extract DATA according to its structure
function extract_DATA(
    DATA::Vector{Pair{IntRange,Vector{TokenValueLine}}};
    symbol="£",
    parser=identity,
    transformer=identity,
    block_merger=identity
    )::Vector{Pair{IntRange,Any}}
    return [bx=>block_merger([transformer([parser(v) for (k,v) ∈ d if k==symbol]) for d ∈ data])
            for (bx,data) ∈ DATA]
end

max_stars_repo_name: algorithmx/SmartParser

max_issues_repo_name: algorithmx/SmartParser

max_forks_repo_name: algorithmx/SmartParser


text:
#=
test_df:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-14
=#

module TestDf
    # packages
    using Test

    # external modules
    include("create_df.jl")
    include("operate_df.jl")

    # methods
    function test()
        @testset "DataFrames" begin
            @testset "DfCreator" begin
                @test_nowarn DfCreator.main()
            end
            @testset "DfOperator" begin
                @test_nowarn DfOperator.main()
            end
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .TestDf
    TestDf.test()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
function asciiOutput(outfile::String = "")
    asciiRegressandTransform(s::String,colmin::Int64,colmax::Int64) = "$s"
    asciiTableHeader(numberOfResults::Int64, align::String) = ""
    asciiTableFooter(numberOfResults::Int64, align::String) = ""
    asciiHeaderRule(headerCellStartEnd::Vector{Vector{Int64}}) = "-"
    toprule = "-"
    midrule = "-"
    bottomrule = "-"
    headerrule = asciiHeaderRule
    headercolsep = "   "
    colsep = "   "
    linestart = ""
    linebreak = ""

    label_fe_yes = "Yes"
    label_fe_no = ""

    label_statistic_n = "N"
    label_statistic_r2 = "R2"
    label_statistic_adjr2 = "Adjusted R2"
    label_statistic_r2_within = "Within-R2"
    label_statistic_f = "F"
    label_statistic_p = "F-test p value"
    label_statistic_f_kp = "First-stage F statistic"
    label_statistic_p_kp = "First-stage p value"
    label_statistic_dof = "Degrees of Freedom"

    label_estimator = "Estimator"
    label_estimator_ols = "OLS"
    label_estimator_iv = "IV"
    label_estimator_nl = "NL"

    foutfile = outfile
    encapsulateRegressand = asciiRegressandTransform
    header = asciiTableHeader
    footer = asciiTableFooter
    return RenderSettings(toprule, midrule, bottomrule, headerrule, headercolsep, colsep, linestart,
        linebreak, label_fe_yes, label_fe_no,
        label_statistic_n, label_statistic_r2, label_statistic_adjr2, label_statistic_r2_within,
        label_statistic_f, label_statistic_p, label_statistic_f_kp, label_statistic_p_kp, label_statistic_dof,
        label_estimator, label_estimator_ols, label_estimator_iv, label_estimator_nl,
        foutfile, encapsulateRegressand, header, footer)
end


text:
################################################################################
#
#    This file contains methods to calculate the bispectrum of a given Configuration
#       Requires Configuration (see utils.jl). All of the command functions use LAMMPS.jl
#
#       The following methods are useful but require fitted snap.β
#       get_bispectrum calculates the bispectrum components using LAMMPS
#       get_dbispectrum calculuates the derivative of the bispectrum components (used for forces)
#       get_vbispectrum calculuates the virial contribution of the bispectrum components (used for virial stress tensor)
#
#       get_snap produces the A matrix in the SNAP potential fitting paradigm
#           A β = b
#           where β is the coefficients of the bispectrum components to produce energies, forces, stresses.
#           For detailed information, see http://dx.doi.org/10.1016/j.jcp.2014.12.018 (Thompson et al. 2014)
#       get_snap is instrumental in the methods that calculate potential_energy, force, virial for the snap potential.
################################################################################

function get_bispectrum(c::Configuration, snap::SNAP; dim = 3)
    num_coeffs = get_num_coeffs(snap.twojmax)

    radii_string = ""
    weight_string = ""
    for j = 1:c.num_atom_types
        radii_string *= string(c.radii[j])
        radii_string *= " "
        weight_string *= string(c.weights[j])
        weight_string *= " "
    end

    A = LMP(["-screen", "none"]) do lmp
        bispectrum = Array{Float64}(undef, num_coeffs, length(c.Atoms))
        for i = 1
            # Initialize
            command(lmp, "log none")
            command(lmp, "units " * c.units)
            command(lmp, "boundary $(c.domain.bound_type[1]) $(c.domain.bound_type[2]) $(c.domain.bound_type[3])")
            command(lmp, "atom_style atomic")
            command(lmp, "atom_modify map array")

            # Setup box
            command(lmp, "region mybox block $(c.domain.bounds[1][1]) $(c.domain.bounds[1][2]) $(c.domain.bounds[2][1]) $(c.domain.bounds[2][2]) $(c.domain.bounds[3][1]) $(c.domain.bounds[3][2])")
            command(lmp, "create_box $(c.num_atom_types) mybox")

            # Create atoms
            println(c.atom_names)
            for j = 1:length(c.Atoms)
                atom_id = findall(c.atom_names .== c.Atoms[j].Type)[1]
                command(lmp, "create_atoms $atom_id single $(c.Atoms[j].Position[1]) $(c.Atoms[j].Position[2]) $(c.Atoms[j].Position[3])")
            end

            if c.units == "lj"
                command(lmp, "mass 1 1.0")
            else
                for j = 1:c.num_atom_types
                    command(lmp, "mass $j $(c.masses[j])")
                end
            end

            # Setup Forcefield
            cutoff = snap.rcutfac * maximum(c.radii)
            max_bounds = max(2*cutoff, min( (c.domain.bounds[1][2] - c.domain.bounds[1][1]), (c.domain.bounds[2][2] - c.domain.bounds[2][1]), (c.domain.bounds[3][2] - c.domain.bounds[3][1]) ) )
            command(lmp, "pair_style zero $max_bounds")
            command(lmp, "pair_coeff * *")
            command(lmp, "compute PE all pe")
            command(lmp, "compute S all pressure thermo_temp")
            string_command = "compute b all sna/atom $(snap.rcutfac) 0.99363 $(snap.twojmax) " * radii_string * weight_string

            args_command = ""
            fields = fieldnames(typeof(snap.keywords))
            for f in fields
                name = string(f)
                value = getfield(snap.keywords, f)
                if value == 0
                    continue
                else
                    value = string(value)
                    if (name == "chemflag")
                        value = string(c.num_atom_types) * " "
                        for j = 0:(c.num_atom_types-1)
                            value *= string(j) * " "
                        end
                    end
                    args_command *= name * " " * value * " "
                end
            end
            string_command *= args_command
            command(lmp, strip(string_command))
            if dim == 2
                command(lmp, "fix lo all enforce2d")
            end
            command(lmp, "thermo_style custom pe")
            command(lmp, "run 0")

            ## Extract bispectrum
            bs = extract_compute(lmp, "b",  LAMMPS.API.LMP_STYLE_ATOM,
                                            LAMMPS.API.LMP_TYPE_ARRAY)
            bispectrum[:, :] = bs
            command(lmp, "clear")
        end
        return bispectrum
    end
    return A

end

function get_dbispectrum(c::Configuration, snap::SNAP; dim = 3)
    num_coeffs = get_num_coeffs(snap.twojmax)

    radii_string = ""
    weight_string = ""
    for j = 1:c.num_atom_types
        radii_string *= string(c.radii[j])
        radii_string *= " "
        weight_string *= string(c.weights[j])
        weight_string *= " "
    end

    A = LMP(["-screen", "none"]) do lmp
        bispectrum = Array{Float64}(undef, 3*c.num_atom_types*num_coeffs, length(c.Atoms))
        for i = 1
            # Initialize
            command(lmp, "log none")
            command(lmp, "units " * c.units)
            command(lmp, "boundary $(c.domain.bound_type[1]) $(c.domain.bound_type[2]) $(c.domain.bound_type[3])")
            command(lmp, "atom_style atomic")
            command(lmp, "atom_modify map array")

            # Setup box
            command(lmp, "region mybox block $(c.domain.bounds[1][1]) $(c.domain.bounds[1][2]) $(c.domain.bounds[2][1]) $(c.domain.bounds[2][2]) $(c.domain.bounds[3][1]) $(c.domain.bounds[3][2])")
            command(lmp, "create_box $(c.num_atom_types) mybox")

            # Create atoms
            for j = 1:length(c.Atoms)
                atom_id = findall(c.atom_names .== c.Atoms[j].Type)[1]
                command(lmp, "create_atoms $atom_id single $(c.Atoms[j].Position[1]) $(c.Atoms[j].Position[2]) $(c.Atoms[j].Position[3])")
            end

            if c.units == "lj"
                command(lmp, "mass 1 1.0")
            else
                for j = 1:c.num_atom_types
                    command(lmp, "mass $j $(c.masses[j])")
                end
            end

            # Setup Forcefield
            cutoff = snap.rcutfac * maximum(c.radii)
            max_bounds = max(2*cutoff, min( (c.domain.bounds[1][2] - c.domain.bounds[1][1]), (c.domain.bounds[2][2] - c.domain.bounds[2][1]), (c.domain.bounds[3][2] - c.domain.bounds[3][1]) ) )
            command(lmp, "pair_style zero $max_bounds")
            command(lmp, "pair_coeff * *")
            command(lmp, "compute PE all pe")
            command(lmp, "compute S all pressure thermo_temp")
            string_command = "compute db all snad/atom $(snap.rcutfac) 0.99363 $(snap.twojmax) " * radii_string * weight_string
            args_command = ""
            fields = fieldnames(typeof(snap.keywords))
            for f in fields
                name = string(f)
                value = getfield(snap.keywords, f)
                if value == 0
                    continue
                else
                    value = string(value)
                    if (name == "chemflag")
                        value = string(c.num_atom_types) * " "
                        for j = 0:(c.num_atom_types-1)
                            value *= string(j) * " "
                        end
                    end
                    args_command *= name * " " * value * " "
                end
            end
            string_command *= args_command
            command(lmp, strip(string_command))
            if dim == 2
                command(lmp, "fix lo all enforce2d")
            end
            command(lmp, "thermo_style custom pe")
            command(lmp, "run 0")

            ## Extract bispectrum
            bs = extract_compute(lmp, "db",  LAMMPS.API.LMP_STYLE_ATOM,
                                            LAMMPS.API.LMP_TYPE_ARRAY)
            bispectrum[:, :] = bs
            command(lmp, "clear")
        end
        return bispectrum
    end
    return A

end

function get_vbispectrum(c::Configuration, snap::SNAP; dim = 3)
    num_coeffs = get_num_coeffs(snap.twojmax)

    radii_string = ""
    weight_string = ""
    for j = 1:c.num_atom_types
        radii_string *= string(c.radii[j])
        radii_string *= " "
        weight_string *= string(c.weights[j])
        weight_string *= " "
    end

    A = LMP(["-screen", "none"]) do lmp
        bispectrum = Array{Float64}(undef, 6*c.num_atom_types*num_coeffs, length(c.Atoms))
        for i = 1
            # Initialize
            command(lmp, "log none")
            command(lmp, "units " * c.units)
            command(lmp, "boundary $(c.domain.bound_type[1]) $(c.domain.bound_type[2]) $(c.domain.bound_type[3])")
            command(lmp, "atom_style atomic")
            command(lmp, "atom_modify map array")

            # Setup box
            command(lmp, "region mybox block $(c.domain.bounds[1][1]) $(c.domain.bounds[1][2]) $(c.domain.bounds[2][1]) $(c.domain.bounds[2][2]) $(c.domain.bounds[3][1]) $(c.domain.bounds[3][2])")
            command(lmp, "create_box $(c.num_atom_types) mybox")

            # Create atoms
            for j = 1:length(c.Atoms)
                atom_id = findall(c.atom_names .== c.Atoms[j].Type)[1]
                command(lmp, "create_atoms $atom_id single $(c.Atoms[j].Position[1]) $(c.Atoms[j].Position[2]) $(c.Atoms[j].Position[3])")
            end

            if c.units == "lj"
                command(lmp, "mass 1 1.0")
            else
                for j = 1:c.num_atom_types
                    command(lmp, "mass $j $(c.masses[j])")
                end
            end

            # Setup Forcefield
            cutoff = snap.rcutfac * maximum(c.radii)
            max_bounds = max(2*cutoff, min( (c.domain.bounds[1][2] - c.domain.bounds[1][1]), (c.domain.bounds[2][2] - c.domain.bounds[2][1]), (c.domain.bounds[3][2] - c.domain.bounds[3][1]) ) )
            command(lmp, "pair_style zero $max_bounds")
            command(lmp, "pair_coeff * *")
            command(lmp, "compute PE all pe")
            command(lmp, "compute S all pressure thermo_temp")
            string_command = "compute vb all snav/atom $(snap.rcutfac) 0.99363 $(snap.twojmax) " * radii_string * weight_string
            args_command = ""
            fields = fieldnames(typeof(snap.keywords))
            for f in fields
                name = string(f)
                value = getfield(snap.keywords, f)
                if value == 0
                    continue
                else
                    value = string(value)
                    if (name == "chemflag")
                        value = string(c.num_atom_types) * " "
                        for j = 0:(c.num_atom_types-1)
                            value *= string(j) * " "
                        end
                    end
                    args_command *= name * " " * value * " "
                end
            end
            string_command *= args_command
            command(lmp, strip(string_command))
            if dim == 2
                command(lmp, "fix lo all enforce2d")
            end
            command(lmp, "thermo_style custom pe")
            command(lmp, "run 0")

            ## Extract bispectrum
            bs = extract_compute(lmp, "vb",  LAMMPS.API.LMP_STYLE_ATOM,
                                            LAMMPS.API.LMP_TYPE_ARRAY)
            bispectrum[:, :] = bs
            command(lmp, "clear")
        end
        return bispectrum
    end
    return A

end

function get_snap(c::Configuration, snap::SNAP; dim = 3, reference = true)
    J = snap.twojmax
    num_coeffs = length(snap.β)

    radii_string = ""
    weight_string = ""
    for j = 1:c.num_atom_types
        radii_string *= string(c.radii[j])
        radii_string *= " "
        weight_string *= string(c.weights[j])
        weight_string *= " "
    end

    A = LMP(["-screen", "none"]) do lmp
        bispectrum = Array{Float64}(undef, num_coeffs, 1+3*length(c.Atoms)+6)
        for i = 1
            # Initialize
            command(lmp, "log none")
            command(lmp, "units " * c.units)
            command(lmp, "boundary $(c.domain.bound_type[1]) $(c.domain.bound_type[2]) $(c.domain.bound_type[3])")
            command(lmp, "atom_style atomic")
            command(lmp, "atom_modify map array")

            # Setup box
            command(lmp, "region mybox block $(c.domain.bounds[1][1]) $(c.domain.bounds[1][2]) $(c.domain.bounds[2][1]) $(c.domain.bounds[2][2]) $(c.domain.bounds[3][1]) $(c.domain.bounds[3][2])")
            command(lmp, "create_box $(c.num_atom_types) mybox")

            # Create atoms
            for j = 1:length(c.Atoms)
                atom_id = findall(c.atom_names .== c.Atoms[j].Type)[1]
                command(lmp, "create_atoms $atom_id single $(c.Atoms[j].Position[1]) $(c.Atoms[j].Position[2]) $(c.Atoms[j].Position[3])")
            end

            if c.units == "lj"
                command(lmp, "mass 1 1.0")
            else
                for j = 1:c.num_atom_types
                    command(lmp, "mass $j $(c.masses[j])")
                end
            end

            # Setup Forcefield
            cutoff = snap.rcutfac * maximum(c.radii)
            max_bounds = max(2*cutoff, min( (c.domain.bounds[1][2] - c.domain.bounds[1][1]), (c.domain.bounds[2][2] - c.domain.bounds[2][1]), (c.domain.bounds[3][2] - c.domain.bounds[3][1]) ) )
            if reference
                command(lmp, "pair_style hybrid/overlay zero $max_bounds zbl 2.0 2.5")
                command(lmp, "pair_coeff * * zero")
                command(lmp, "pair_coeff * * zbl 1 1")
            else
                command(lmp, "pair_style zero $max_bounds")
                command(lmp, "pair_coeff * *")
            end
            command(lmp, "compute PE all pe")
            command(lmp, "compute S all pressure thermo_temp")
            string_command = "compute snap all snap $(snap.rcutfac) 0.99363 $(snap.twojmax) " * radii_string * weight_string
            args_command = ""
            fields = fieldnames(typeof(snap.keywords))
            for f in fields
                name = string(f)
                value = getfield(snap.keywords, f)
                if value == 0
                    continue
                else
                    value = string(value)
                    if (name == "chemflag")
                        value = string(c.num_atom_types) * " "
                        for j = 0:(c.num_atom_types-1)
                            value *= string(j) * " "
                        end
                    end
                    args_command *= name * " " * value * " "
                end
            end
            string_command *= args_command
            command(lmp, strip(string_command))
            if dim == 2
                command(lmp, "fix lo all enforce2d")
            end
            command(lmp, "thermo_style custom pe")
            command(lmp, "run 0")

            ## Extract bispectrum
            bs = extract_compute(lmp, "snap",  LAMMPS.API.LMP_STYLE_GLOBAL,
                                            LAMMPS.API.LMP_TYPE_ARRAY)
            bispectrum[:, :] = bs
            bispectrum[end, 1] = length(c.Atoms)
            command(lmp, "clear")
        end
        return bispectrum
    end

    return copy(transpose(A))
end

function get_snap(r::Vector{Configuration}, p; dim = 3, reference = true)
    n = length(r)
    l = length(p.β)
    Aenergy = Array{Float64}(undef, 0, l)
    Aforce = Array{Float64}(undef, 0, l)
    Astress = Array{Float64}(undef, 0, l)
    for j = 1:n
       A = get_snap(r[j], p; dim = dim, reference = reference)
       Aenergy = vcat(Aenergy, reshape(A[1, :], 1, l))
       Aforce = vcat(Aforce, A[2:end-6, :])
       Astress = vcat(Astress, A[end-5:end, :])
    end
    A = vcat(Aenergy, Aforce, Astress)
    return A
end


text:
module Utils

using Transducers: DistributedEx, ThreadedEx
using Test

in_ci() = lowercase(get(ENV, "CI", "false")) == "true"

struct NoError end

macro test_error(thrower)
    @gensym err tmp
    quote
        local $err = $NoError()
        $Test.@test try
            $thrower
            false
        catch $tmp
            $Base.@debug(
                "@test_error",
                exception = ($tmp, $catch_backtrace()),
                _module = $__module__,
                _file = $(string(__source__.file)),
                _line = $(__source__.line),
            )
            $err = $tmp
            true
        end
        $err
    end |> esc
end

function foreach_executor(f)
    @testset for ex in [ThreadedEx(), DistributedEx()]
        @debug "`$(parentmodule(f)).$f`" ex
        f(ex)
    end
end

function with_racer(f, label::Symbol; timeout = in_ci() ? 5 * 60 : 5)
    startswith(string(label), "_") && error("label must not start with __")
    endswith(string(label), "_") && error("label must not end with __")
    race = Channel{Symbol}(3) do race
        sleep(timeout)
        put!(race, :__timeout__)
        @error "`with_racer(_, :$label)`: TIMEOUT"
        close(race)
    end
    function win()
        put!(race, label)
    end
    local race_results
    try
        f(win)
    finally
        put!(race, :__done__)
        close(race)
        race_results = collect(race)
        @test race_results == [label, :__done__]
    end
    return race_results == [label, :__done__]
end

function with_racer(f, label::Symbol, rest::Symbol...; kwargs...)
    ok2 = false
    ok1 = with_racer(label; kwargs...) do win
        function wrapper(wins...)
            f(win, wins...)
        end
        ok2 = with_racer(wrapper, rest...; kwargs...)
    end
    return ok1 & ok2
end

print_error(err) = print_error(stdout, err)

function print_error(io, err::CompositeException)
    for (i, e) in enumerate(err)
        println(io)
        println(io, "$i-th error: type = $(typeof(e))")
        showerror(io, e)
    end
    println(io)
end

function print_error(io, err)
    println(io)
    println(io, "Got non-`CompositeException` exception:")
    showerror(io, err)
    println(io)
end

end  # module


text:
inpt = """
1
2
3
7
11
13
17
19
23
31
37
41
43
47
53
59
61
67
71
73
79
83
89
97
101
103
107
109
113"""

packs = Int64[]
for item in split(inpt)
	push!(packs, parse(Int, item))
end
# only need to look for the smallest subset that sums to 520 mass (total/3) and leaves two other subsets that sum to total/3
i = 3
found = false
total = div(sum(packs), 4)
println(total)
while true
	# the input is sorted, so the smallest quantum entanglement is going to be found first
	for best in combinations(packs, i)
		if sum(best) == total
			remainder = setdiff(packs, best)
			for parts in partitions(remainder, 3)
				# now try all possibilities of the remaining combinations
				# if they don't all 3 add up to total, we're not done
				found = true
				for p in parts
					if sum(p) != total
						found = false
						break
					end
				end
				found && break
			end
		end
		if found
			println(best, " ", sum(best), " ", prod(best))
			break
		end
	end
	found && break
	i += 1
end

max_stars_repo_name: jstrube/adventofcode

max_issues_repo_name: jstrube/adventofcode

max_forks_repo_name: jstrube/adventofcode


text:
# This file was generated, do not modify it. # hide
println(read("ejemplos-palabras-knn.md", String))


text:
using Tokens
using Documenter

DocMeta.setdocmeta!(Tokens, :DocTestSetup, :(using Tokens); recursive=true)

makedocs(;
    modules=[Tokens],
    authors="Robert Rudolph",
    repo="https://github.com/rryi/Tokens.jl/blob/{commit}{path}#{line}",
    sitename="Tokens.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://rryi.github.io/Tokens.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/rryi/Tokens.jl",
)


text:
using BinDeps

@BinDeps.setup

ENV["JULIA_ROOT"] = abspath(JULIA_HOME, "../../")

# include alias for WinRPM library
libportaudio = library_dependency("libportaudio", aliases=["libportaudio-2"])

# TODO: add other providers with correct names
provides(AptGet, "libportaudio2", libportaudio)
provides(Pacman, "portaudio", libportaudio)

@static if is_apple()
    using Homebrew
    provides(Homebrew.HB, "portaudio", libportaudio)
end

@static if is_windows()
    using WinRPM
    provides(WinRPM.RPM, "libportaudio2", libportaudio, os = :Windows)
end

@BinDeps.install Dict(:libportaudio => :libportaudio, )


text:
using DashBootstrapComponents

button = html_div(
    [
        dbc_button("Button 1", className = "me-md-2"),
        dbc_button("Button 2", className = "me-md-2"),
        dbc_button("Button 3"),
    ],
    className = "d-grid gap-2 d-md-flex justify-content-md-end",
);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
module WMT
using Fetch
using DataDeps
using InternedStrings
using BytePairEncoding

using ..Datasets: Dataset
import ..Datasets: testfile, trainfile, get_vocab

export GoogleWMT

function __init__()
    googlewmt_init()
end

include("./google_wmt.jl")

end


text:
module Fake

using CairoMakie, MakieLayout, FileIO, AbstractPlotting
CairoMakie.activate!()

set_theme!(font = AbstractPlotting.findfont("Helvetica"))

function bad1()
    scene, layout = layoutscene()
    ax = layout[1,1] = LAxis(scene)
    FileIO.save("a.pdf", scene)
end

function bad2()
    scene = Scene()
    lines!(scene, rand(10))
    linesegments!(scene, rand(10) .+ 1)
    poly!(scene, Circle(Point2f0(0, 4), 1f0))
    FileIO.save("a.pdf", scene)
end

function bad3()
    scene = Scene(; show_axis = false)
    lines!(scene, rand(10))
    linesegments!(scene, rand(10) .+ 1)
    poly!(scene, Circle(Point2f0(0, 4), 1f0))
    FileIO.save("a.pdf", scene)
end

function bad4()
    scene = Scene(; show_axis = false)
    lines!(scene, rand(10))
    linesegments!(scene, rand(10) .+ 1)
    FileIO.save("a.pdf", scene)
end

function good1() # all good
    scene = Scene(; show_axis = false)
    FileIO.save("a.pdf", scene)
end

function good2() # all good
    scene = Scene()
    FileIO.save("a.pdf", scene)
end

function good3()
    scene = Scene(; show_axis = false)
    lines!(scene, rand(10))
    FileIO.save("a.pdf", scene)
end

end


text:
"""Initialize the simulation from the middle of the state space."""
initialize_simple!(vStates0, sol::AbstractDDPSolution) =
	initialize_simple!(vStates0, sol.prob.tStateVectors)
function initialize_simple!(vStates0, tStateVectors::NTuple{N, Vector{Float64}}) where N
	vStates0 .= getindex.(tStateVectors, Int.(floor.(length.(tStateVectors)./2)))
end

initialize_choices(sol::DDPSolutionZero{NS,NC,NE0,NC0}, vExogStates0) where {NS,NC,NE0,NC0} =
	[itp(vExogStates0...) for itp in sol.policy0]
initialize_choices(sol::DDPSolutionZero{NS,NC,NE0,1}, vExogStates0) where {NS,NC,NE0} =
	sol.policy0[1](vExogStates0...)

"""Initialization with shocks and choices"""
function initialize_exact!(vStates0, sol::DDPSolutionZero{NS,NC,NE0,NC0},
		init_opt::InitializationOptions{NC0}, vShocks) where {NS,NC,NE0,NC0}
	vStates0[1+NC0:end] .= init_opt.func(vShocks)
	vStates0[1:NC0] .= initialize_choices(sol, vStates0[1+NC0:end])
end
"""when initial shock is one-dimensional."""
function initialize_exact!(vStates0, sol::DDPSolutionZero{NS,NC,NE0,NC0},
		init_opt::InitializationOptions{NC0,1}, vShocks) where {NS,NC,NE0,NC0}
	vStates0[1+NC0:end] .= init_opt.func(vShocks[1])
	vStates0[1:NC0] .= initialize_choices(sol, vStates0[1+NC0:end])
end
"""Only one choice variable"""
function initialize_exact!(vStates0, sol::DDPSolutionZero{NS,NC,NE0,1},
		init_opt::InitializationOptions{1}, vShocks) where {NS,NC,NE0}
	vStates0[2:end] .= init_opt.func(vShocks)
	vStates0[1] = initialize_choices(sol, vStates0[2:end])
end
"""Only one choice variable and shock is one-dimensional."""
function initialize_exact!(vStates0, sol::DDPSolutionZero{NS,NC,NE0,1},
		init_opt::InitializationOptions{1,1}, vShocks) where {NS,NC,NE0}
	vStates0[2:end] .= init_opt.func(vShocks[1])
	vStates0[1] = initialize_choices(sol, vStates0[2:end])
end

function initialize!(vStates0, initialize_exact::Bool, sol::AbstractDDPSolution,
		init_opt, vShocks)
	if initialize_exact
		initialize_exact!(vStates0, sol, init_opt, vShocks)
	else
		initialize_simple!(vStates0, sol)
	end
end

"""Get optimal choice of a firm at time t."""
function get_choice_t!(vChoice, policy::NTuple{NC,T}, vStates, t) where {NC,T}
	for j = 1 : NC
		vChoice[j,t] = policy[j](vStates[:,t]...)
	end
end

"""Get transition for a firm at time t."""
get_transition_t!(transfunc::Transition,
	vStates, vChoices, vShocks, t) =
	get_transition!(transfunc, @view(vStates[:,t+1]), vStates[:,t],
		vChoices[:,t], vShocks[:,t])
"""Only one shock."""
get_transition_t!(transfunc::Transition{ID,NC,1},
	vStates, vChoices, vShocks, t) where {ID,NC} =
	get_transition!(transfunc, @view(vStates[:,t+1]), vStates[:,t],
		vChoices[:,t], vShocks[:,t][1])
"""Only one choice var."""
get_transition_t!(transfunc::Transition{ID,1},
	vStates, vChoices, vShocks, t) where ID =
	get_transition!(transfunc, @view(vStates[:,t+1]), vStates[:,t],
		vChoices[t], vShocks[:,t])
"""Only one choice var and only one shock."""
get_transition_t!(transfunc::Transition{ID,1,1},
	vStates, vChoices, vShocks, t) where ID =
	get_transition!(transfunc, @view(vStates[:,t+1]), vStates[:,t],
		vChoices[t], vShocks[:,t][1])

# convenience wrappers
simulate(sol::AbstractDDPSolution;
			nPeriods::Int64 = 60,
			nFirms::Int64 = 100, kwargs...) =
	_simulate(sol.prob, sol, drawshocks(sol.prob, nPeriods=nPeriods, nFirms=nFirms); kwargs...)

simulate(sol::AbstractDDPSolution, shocks::AbstractDDPShocks; kwargs...) =
	_simulate(sol.prob, sol, shocks; kwargs...)

simulate(p::DDP; kwargs...) =
	_simulate(p, drawshocks(p); kwargs...)
simulate(p::DDP, shocks::AbstractDDPShocks; kwargs...) =
	_simulate(p, solve(p), shocks; kwargs...)

_simulate(p::DDP, sol::AbstractDDPSolution, shocks::AbstractDDPShocks;
        initialize_exact::Bool = typeof(sol) <: DDPSolutionZero,
		get_value::Bool = false) =
		_simulate(p, sol, shocks, p.transfunc, p.options.initialize,
				initialize_exact, get_value)

function _simulate(p::DDP{NS,NC}, sol::AbstractDDPSolution, shocks::AbstractDDPShocks,
				transfunc::Transition{ID,NC,dimShocks},
				init_opt::Union{InitializationOptions,Nothing},
				initialize_exact::Bool,
				get_value::Bool) where {ID, NS, NC, dimShocks}

	!(initialize_exact && !(typeof(sol) <: DDPSolutionZero)) || error(
		"Solution does not contain intial policy -> specify initializationproblem when defining DDP")

	unused, nPeriods, nFirms = size(shocks.aSim)
	aSim = fill!(zeros(NS, nPeriods+1, nFirms), NaN)
	mVal = fill!(zeros(1, nPeriods+1, nFirms), NaN)
	aChoice = fill!(zeros(NC, nPeriods+1, nFirms), NaN)

	value = sol.value
	policy = sol.policy

	# # for exogenous exit
 	# do_exog_exit = isdefined(params,:π)
	# if do_exog_exit
	# 	do_exog_exit =  do_exog_exit * (params.π>0)
	# end

	for i = 1:nFirms

        vSim_i = @view aSim[:,:,i]
		vVal_i = @view mVal[:,:,i]
		vChoice_i = @view aChoice[:,:,i]
        mShocks_i = @view shocks.aSim[:,:,i]

		if initialize_exact
        	mShocks0_i = @view shocks.aInit[:,1,i]
		else
			mShocks0_i = nothing
		end

		initialize!(@view(vSim_i[:,1]), initialize_exact, sol, init_opt,
				 mShocks0_i)

		for t = 1 : nPeriods

			if get_value
            	vVal_i[t] = value(vSim_i[:,t]...)
			end

			get_choice_t!(vChoice_i, policy, vSim_i, t)

			get_transition_t!(transfunc, vSim_i, vChoice_i, mShocks_i, t)

			# if do_exog_exit
			# 	if i/nFirms > (1-params.π)^max(t-2, 0)
			# 		break
			# 	end
			# end # exit condition

		end # period t

		if get_value
			vVal_i[nPeriods+1] = value(vSim_i[:,nPeriods+1]...)
		end

		get_choice_t!(vChoice_i, policy, vSim_i, nPeriods+1)

	end # firm i

	if get_value
		return DDPSimulation(p, sol, mVal, aSim, aChoice)
	else
		return DDPSimulation(p, sol, nothing, aSim, aChoice)
	end
end # simulatemodel


text:
#   This file is part of Reduce.jl. It is licensed under the MIT license
#   Copyright (C) 2017 Michael Reed

export RExpr, @RExpr, @R_str, rcall, @rcall, convert, ==, string, show, squash, list
import Base: parse, convert, ==, getindex, *, split, string, show, join, String

export ExprSymbol
const ExprSymbol = Union{<:Expr,<:Symbol}

"""
Reduce expression

### Summary:
type RExpr <: Any

### Fields:
str::Array{Compat.String,1}
"""
struct RExpr
    str::Array{String,1}
    RExpr(r::Array{String,1}) = new(r)
end

RExpr(r::Array{SubString{String},1}) = RExpr(convert(Array{String,1},r))
RExpr(str::String) = RExpr([str])
RExpr(t::Bool) = RExpr(t ? 1 : 0)

"""
    RExpr(e::Expr)

Convert Julia expression to Reduce expression

## Examples
```julia-repl
julia> RExpr(:(sin(x*im) + cos(y*ϕ)))

     sqrt(5)*y + y
cos(---------------) + sinh(x)*i
           2
```
"""
RExpr(expr::Expr) = expr |> unparse |> RExpr |> split

function RExpr(r::T) where T <: Union{Array,Adjoint,Tuple,Pair}
    out = IOBuffer()
    show_expr(out,r)
    return out |> take! |> String |> RExpr
end

function RExpr(r::Any)
    typeof(r) <: AbstractFloat && isinf(r) && (return RExpr((r > 0 ? "" : "-")*"infinity"))
    typeof(r) <: Irrational && (return RExpr(unparse_irrational(r)))
    y = "$r"
    for key ∈ keys(repjlr)
        y = replace(y, key => repjlr[key])
    end
    return RExpr(y)
end

macro RExpr(r)
    RExpr(r)
end

macro R_str(str)
    RExpr(str)
end

function rtrim(r::Array{String,1})
    n = deepcopy(r)
    h = 1
    l = length(r)
    while h ≤ l
        isempty(n[h]) ? (deleteat!(n,h); l -= 1) : (h += 1)
    end
    return n
end

function split(r::RExpr)
    n = String[]
    for h ∈ 1:length(r.str)
        p = split(replace(r.str[h],r"(\$)|(;\n)"=>";"),r"(?<!\!#[0-9a-fA-F]{4});")
        for t ∈ 1:length(p)
            push!(n,p[t])
    end; end
    return RExpr(rtrim(n))
end

String(r::RExpr) = convert(String,r)
string(r::RExpr) = convert(String,r)
join(r::RExpr) = RExpr(string(r))
join(r::Array{RExpr,1}) = vcat(convert.(Array{String,1},r)...) |> RExpr
show(io::IO, r::RExpr) = print(io,convert(String,r))

cols = 80

"""
    linelength()

This operator is used with the syntax
```Julia
Reduce.linelength()::Integer
```
and sets the output line length to the integer `tput cols`. It returns the output line length (so that it can be stored for later resetting of the output line if needed).
"""
function linelength()
    c = displaysize(stdout)[2]
    global cols
    if c ≠ cols
        ws = rcall("ws")
        rcall("linelength($c)")
        rcall(ws)
        cols = c
    end
    return c
end

function show(io::IO, ::MIME"text/plain", r::RExpr)
    length(r.str) > 1 && (print(io,string(r),";"); return nothing)
    ColCheck() && linelength()
    print(io,rcall(r;on=[:nat]) |> string |> chomp)
end

function show(io::IO, ::MIME"text/latex", r::RExpr)
    rcall(R"on latex")
    write(rs,r)
    rd = readsp(rs)
    rcall(R"off latex")
    sp = split(join(rd),"\n\n")
    print(io,"\\begin{eqnarray}\n")
    ct = 0 # add enumeration
    for str ∈ sp
        ct += 1
        length(sp) ≠ 1 && print(io,"($ct)"*'&'*"\\,")
        print(io,replace(str,r"(\\begin{displaymath})|(\\end{displaymath})"=>""))
        ct ≠ length(sp) && print(io,"\\\\\\\\")
    end # new line
    print(io,"\n\\end{eqnarray}")
end

function extract end

"""
    Reduce.@subtype

Can be used to create an `RExpr` wrapper type with a subtype relation
```Julia
julia> Reduce.@subtype FakeReal <: Real

julia> FakeReal(:(x+1)) + FakeReal(:y)
y + 1 + x
```
"""
macro subtype(x)
    x.head ≠ :<: && throw(error("$x is not a subtype expression"))
    name = x.args[1]
    Expr(:struct,false,x,Expr(:block,Expr(:(::), :r, :RExpr))) |> eval
    Expr(:block,[:($fun(r::$name) = $fun(r.r)) for fun ∈ [:String,:string,:join,:list,:parse,:mat]]...) |> eval
    @eval begin
        export $name
        $name(x) = $name(RExpr(x))
        RExpr(r::$name) = r.r
        show(io::IO, r::$name) = show(io,r.r)
        ==(a::$name,b::$name) = a.r == b.r
        rcall(r::$name,s...) = rcall(r.r,s...)
        convert(::Type{$name}, r::RExpr) = r.r
        convert(::Type{Array{String,1}}, r::$name) = r.r.str
        convert(::Type{String}, r::$name) = convert(String,r.r)
        Algebra.init_subtype($name)
    end
    nothing
end

const r_to_jl = Dict(
    "i"             =>  "im",
    "infinity"      =>  "Inf"
)

const r_to_jl_utf = Dict(
    "pi"            =>  "π",
    "e"             =>  "ℯ",
    "euler_gamma"   =>  "MathConstants.γ",
    "golden_ratio"  =>  "MathConstants.φ"
)

const r_to_jl_ifx = Dict(
    "**"            =>  "^",
    "/"             =>  "/"
    #"~"             =>  "\""
)

const jl_to_r = Dict(
    "im"            =>  "i",
    "eu"            =>  "euler_gamma",
    "eulergamma"    =>  "euler_gamma",
    "golden"        =>  "golden_ratio",
    "Inf"           =>  "infinity"
)

const jl_to_r_utf = Dict(
    "π"             =>  "pi",
    "ℯ"             =>  "e",
    "MathConstants.γ"=> "euler_gamma",
    "MathConstants.φ"=> "golden_ratio"
)

const jl_to_r_ifx = Dict(
    "//"            =>  "/"
    #"\""            =>  "~"
)

"""
    list(r)

The operator `list` is an alternative to the usage of curly brackets. `list` accepts an arbitrary number of arguments and returns a list of its arguments. This operator is useful in cases where operators have to be passed as arguments. E.g.,
```Julia
julia> list(:a,list(list(:b,:c),:d),:e) == R"{{a},{{b,c},d},e}"
true
```
"""
list(r::T) where T <: Tuple = RExpr(r)
list(r::Array{RExpr,1}) = "{$(replace(join(split(join(r)).str,','),":="=>"="))}" |> RExpr
list(a::T) where T <: Vector = length(a) ≠ 0 ? list(lister.(a)) : R"{}"
list(a::T) where T <: Adjoint = list([a...])
list(a::T) where T <: Matrix = list([a[:,k] for k ∈ 1:size(a)[2]])
list(r...) = list(r)
lister(expr) = typeof(expr) <: Vector ? list(expr) : RExpr(expr)

export sub_list

function sub_list(syme::Dict{String,String})
    str = IOBuffer()
    write(str,"{")
    k = length(keys(syme))
    for key in keys(syme)
        k -= 1
        write(str,"$key => $(syme[key])")
        k > 0 && write(str,",")
    end
    write(str,"}")
    return String(take!(str))
end

_syme(syme::Dict{String,String}) = sub_list(syme)[2:end-1]

function _subst(syme::String,expr::T) where T
    convert(T, "!*hold($expr)\$ ws where $syme" |> rcall)
end

const symrjl = _syme(r_to_jl)
const symjlr = _syme(jl_to_r)
const reprjlu = r_to_jl_utf
const repjlru = jl_to_r_utf
reprjl = r_to_jl_ifx
const repjlr = jl_to_r_ifx
const gexrjl = Regex("($(join(keys(r_to_jl),")|(")))")
const gexjlr = Regex("($(join(keys(jl_to_r),")|(")))")

"""
    Reduce.Rational(::Bool)

Toggle whether to use '/' or '//' for division in julia expressions
"""
Rational = ( () -> begin
        gs = false
        return (tf=gs)->(gs≠tf && (gs=tf; reprjl["/"]=gs ? "//" : "/"); return gs)
    end)()

"""
    Reduce.SubCall(::Bool)

Toggle whether to substitute additional expressions
"""
SubCall = ( () -> begin
        gs = true
        return (tf=gs)->(gs≠tf && (gs=tf); return gs)
    end)()

"""
    Reduce.SubHold(::Real)

Sleep timer in case of clogged Reduce pipe on SubCall
"""
SubHold = ( () -> begin
        gs = 1/17
        return (tf=gs)->(gs≠tf && (gs=tf); return gs)
    end)()

"""
    Reduce.SubFail(::Integer)

Failure limit in case of clogged Reduce pipe on SubCall
"""
SubFail = ( () -> begin
        gs = 17
        return (tf=gs)->(gs≠tf && (gs=tf); return gs)
    end)()

"""
    Reduce.ColCheck(::Bool)

Toggle whether to reset REPL linewidth on each show
"""
ColCheck = ( () -> begin
        gs = true
        return (tf=gs)->(gs≠tf && (gs=tf); return gs)
    end)()

"""
    Reduce.PrintLog(::Bool)

Toggle whether to display the log of REDUCE commands
"""
PrintLog = ( () -> begin
        gs = false
        return (tf=gs)->(gs≠tf && (gs=tf); return gs)
    end)()

"""
    Reduce.ListPrint(::Int)

Toggle whether to translate assignment as `:=` or `=` for list parsing.
"""
ListPrint = ( () -> begin
        gs = 0
        return (tf=gs)->(gs≠tf && (gs=tf); return gs)
    end)()

@inline function SubReplace(sym::Symbol,str::String;utf=false)
    a = collect((m.match for m = eachmatch(r"([^ ()+*\^\/-]+|[()+*\^\/-])",str)))
    for s ∈ 1:length(a)
        if !isinfix(a[s]) && !occursin(r"[()]",a[s])
            if utf
                for key in keys(sym == :r ? reprjlu : repjlru)
                    if a[s] == key
                        a[s] = replace(a[s],key=>(sym == :r ? reprjlu[key] : repjlru[key]))
                    end
                end
            else
                if occursin(sym == :r ? gexrjl : gexjlr,a[s])
                    w = _subst(sym == :r ? symrjl : symjlr, a[s])
                    if w == ""
                        c = 1
                        f = SubFail()
                        h = SubHold()
                        while w == "" && c < f
                            sleep(sqrt(c)*h)
                            w = _subst(sym == :r ? symrjl : symjlr, a[s])
                            c += 1
                        end
                        PipeClogged(w ≠ "", c, "substitution")
                    end
                    w ≠ "" ? (a[s] = w) : (@info "If this is a recurring problem, try with `Reduce.SubCall(false)`.")
                end
            end
        elseif isinfix(a[s]) && utf
            if (s ≠ length(a)) && (a[s+1] == a[s]) && (a[s] ∈ ["*","/"])
                a[s] *= a[s+1]
                a[s+1] = ""
            end
            for key in keys(sym == :r ? reprjl : repjlr)
                if a[s] == key
                    a[s] = replace(a[s],key=>(sym == :r ? reprjl[key] : repjlr[key]))
                end
            end
        end
        if sym == :r
            a[s] == "inf" && (a[s] = "Inf")
            a[s] == " - inf" && (a[s] = "-Inf")
            (a[s] == "nan") | (a[s] == " - nan") && (a[s] = "NaN")
        end
    end
    return join(a)
end

@noinline function JSymReplace(str::String)
    str = SubReplace(:jl,str;utf=true)
    SubCall() && !isinfix(str) && (str = SubReplace(:jl,str;utf=false))
    occursin("!#",str) && (str = replace(rcall(str,:nat),r"\n"=>""))
    return str
end

@noinline function RSymReplace(str::String)
    clean = replace(str,r"[ ;\n]"=>"")
    paren = occursin(r"^\(((?>[^\(\)]+)|(?R))*\)$",clean)
    (isempty(clean)|(clean=="()")) && (return str)
    if SubCall() && !isinfix(str) && !occursin(str,".")
        str = SubReplace(:r,str;utf=false)
    end
    if occursin("!#",str)
        rsp = split(str,';')
        for h in 1:length(rsp)
            psr = split(rsp[h],'.')
            mod = false
            for k in 1:length(psr)
                if occursin("!#",psr[k])
                    sp = split(psr[k],r"!#")
                    psr[k] = join([sp[1],replace(rcall("!#"*sp[end]*";",:nat),r"\n"=>"")])
                    mod = true
                end
            end
            mod && (rsp[h] = join(psr))
        end
        str = join(rsp)
    end
    str = SubReplace(:r,str;utf=true)
    return paren ? "("*str*")" : str
end

RSymReplace(str::SubString{String}) = str |> String |> RSymReplace

convert(::Type{RExpr}, r::RExpr) = r
convert(::Type{Array{String,1}}, r::RExpr) = r.str
convert(::Type{String}, r::RExpr) = join(r.str,";\n")
convert(::Type{T}, r::RExpr) where T = T <: Number ? eval(parse(r)) : parse(r)

"""
    rcall(r::RExpr)

Evaluate a Reduce expression.

## Examples
```julia-repl
julia> R\"int(sin(x), x)\" |> RExpr |> rcall
 - cos(x)
```
"""
@noinline function rcall(r::RExpr;
        on::Union{Array{Symbol,1},Array{String,1}}=Symbol[],
        off::Union{Array{Symbol,1},Array{String,1}}=Symbol[])
    typeof(on) == Array{String,1} ? (ona = Symbol.(on)) : (ona = on)
    typeof(off) == Array{String,1} ? (offa = Symbol.(off)) : (offa = off)
    ons = IOBuffer()
    onr = IOBuffer()
    offs = IOBuffer()
    offr = IOBuffer()
    mode = true
    trim = false
    expo = false
    rlfi = false
    for o in ona
        if o == :expand
            write(ons,"on exp\$ ")
            write(onr,"; off exp ")
        elseif o == :latex
            rcall(R"on latex")
            rlfi = true
        else
            write(ons,"on $o\$ ")
            write(onr,"; off $o ")
        end
        o == :factor && (expo = true)
        o in offa && throw(ReduceError("Invalid: switch on and off at once"))
        o in [:latex,:nat] && (mode = false)
        o == :nat && (trim = true)
    end
    for o in offa
        !(o == :factor) && write(offs,"off $o\$ ")
        !(o in [offlist;[:factor]]) && write(offr,"; on $o")
    end
    wrs = String(UInt8[take!(ons)...,take!(offs)...]) *
          string(r) *
          String(UInt8[take!(onr)...,take!(offr)...])
    PrintLog() && println(wrs)
    write(rs,wrs)
    sp = mode ? readsp(rs) : read(rs)
    expo && rcall(R"off exp")
    mode && for h ∈ 1:length(sp)
        sp[h] = replace(sp[h],'\n' => "")
        sp[h] = replace(sp[h],'\\' => "")
    end
    trim && (return join(split(sp,"\n")[2:end-1],'\n'))
    rlfi && rcall(R"off latex")
    for o in offa
        o == :nat && (return join(sp))
        o == :latex && popfirst!(sp)
    end
    return mode ? sp |> RExpr |> split : sp
end

rcall(r::RExpr,switches...) = rcall(r;on=[switches...])

"""
    rcall{T}(e::T)

Evaluate a Julia expression or string using the Reduce interpretor and convert
output back into the input type

## Examples
```julia-repl
julia> rcall(\"int(sin(y)^2, y)\")
\"( - cos(y)*sin(y) + y)/2\"

julia> rcall(:(int(1/(1+x^2), x)))
:(atan(x))
```
"""
function rcall(expr::T;on::Array{A,1}=Symbol[],off::Array{B,1}=Symbol[]) where T where A <: Union{Symbol,String} where B <: Union{Symbol,String}
    comp = rcall(RExpr(expr);on=on,off=off)
    (:latex in on) | (:nat in on) ? (return comp) : (return convert(T,comp))
end

macro rcall(r,on::Union{Array{Symbol,1},Array{String,1}}=Symbol[],off::Union{Array{Symbol,1},Array{String,1}}=Symbol[])
    return Expr(:quote,rcall(r.head == :quote ? r.args[1] : r; on=on,off=off))
end

macro rcall(r,switches...)
    Expr(:quote,rcall(r.head == :quote ? r.args[1] : r, switches...))
end

rcall(r,switches...) = rcall(r;on=Symbol.([switches...]))

function ==(r::RExpr, s::RExpr)
    n = expand(r).str
    m = expand(s).str
    l=length(n)
    l≠length(m) && (return false)
    b = true
    for j∈1:l
        b &= "if($(n[j]))=($(m[j]))then 1 else 0"|>rcall|>Meta.parse|>eval|>Bool
    end
    return b
end

#getindex(r::RExpr, i) = "$r($i)" |> rcall

#= regroup parens * add feature

@noinline function squash(expr)
    if typeof(expr) == Expr && expr.head ∈ [:block,:function]
        nex = deepcopy(expr)
        k = expr.head ≠ :function ? 1 : 2
        while k ≤ length(nex.args)
            found = false
            if typeof(nex.args[k]) == Expr &&
                    occursin(r"[*\/+-^]=$",string(nex.args[k].head))
                var = nex.args[k].args[1]
                for h ∈ 1:k-1
                    if typeof(nex.args[h]) == Expr &&
                            nex.args[h].head == :(=) &&
                            nex.args[h].args[1] == var
                        nex.args[h].args[2] = eval(Expr(:call,
                            Symbol(match(r"[^(=$)]",string(nex.args[k].head)).match),
                            QuoteNode(nex.args[h].args[2]),
                            QuoteNode(nex.args[k].args[2])))
                        deleteat!(nex.args,k)
                        found = true
                        break
                    end
                end
            else
                typeof(nex.args[k]) == Expr && (nex.args[k] = squash(nex.args[k]))
            end
            !found && (k += 1)
        end
        return length(nex.args) > 1 ? nex : nex.args[end]
    elseif typeof(expr)
    else
        return expr
    end
end=#

"""
    squash(expr)

Reduces an entire program statement block using symbolic rewriting
"""
function squash(expr)
    typeof(expr) == Expr && if expr.head == :block
        return @eval Reduce.Algebra $expr
    elseif expr.head == :function
        out = deepcopy(expr)
        out.args[2] = @eval Reduce.Algebra $(Expr(:block,expr.args[2]))
        return out
    else
        return rcall(expr)
    end
end


text:
"""
    EratosthenesRotations

This module includes several rotations and conversions used by the Eratosthenes
simulation engine. For now, only the basic rotation tools necessary for the job
are included. In time, it will become more comprehensive, documented, and
faster.

"""
module EratosthenesRotations

# Questions:
#
# Is a Quat a type, a vector, or a matrix?
# Use arrays of vectors/quats or 2D arrays? Since nothing's vectorized, it seems perfectly easy to do the former.
# Use StaticArrays for speed (very likely)?

export qrot, aa2q, q2aa, qinv, qpos, qcomp, qdiff, qdot, mrp2q

# Make an in-place, batch-oriented function. We could do this in general for all
# of the functions below, but it's not necessary yet. Is this even the style I
# would want to go forward with?
function qrot!(q::Array{T,N}, v::Array{T,N}) where {T,N}
    @assert(size(q, 1) == 4, "First input must be 4-by-n.")
    @assert(size(v, 1) == 3, "Second input must be 3-by-n.")
    @assert(size(q, 2) == size(v, 2), "Quaternion and vector must have same number of elements.")
    for k = 1:size(v,2)
        c1 =  q[4,k]*v[1,k] + q[3,k]*v[2,k] - q[2,k]*v[3,k]
        c2 = -q[3,k]*v[1,k] + q[4,k]*v[2,k] + q[1,k]*v[3,k]
        c3 =  q[2,k]*v[1,k] - q[1,k]*v[2,k] + q[4,k]*v[3,k]
        c4 = -q[1,k]*v[1,k] - q[2,k]*v[2,k] - q[3,k]*v[3,k]
        v[1,k] = -c4*q[1,k] - c3*q[2,k] + c2*q[3,k] + c1*q[4,k]
        v[2,k] =  c3*q[1,k] - c4*q[2,k] - c1*q[3,k] + c2*q[4,k]
        v[3,k] = -c2*q[1,k] + c1*q[2,k] - c4*q[3,k] + c3*q[4,k]
    end
end

function qrot(q::Array{T,N}, v::Array{T,N}) where {T,N}
    v2 = copy(v)
    qrot!(q, v2)
    return v2
end

function aa2q(theta::T, r_hat::Vector{T}) where {T<:Real}
    ht = 0.5*theta
    q = [sin(ht) * r_hat; cos(ht)]
end

# function aa2q(theta::Vector{T}, r_hat::Array{T,2}) where {T<:Real}
#     q = zeros(T, (4, length(theta)))
#     for k = 1:size(theta,1)
#         q[:,k] = aa2q(theta[k], r_hat[:,k])
#     end
#     return q
# end

function q2aa(q::Vector{T}) where {T<:Real}
    half_theta = acos(max(min(q[4], 1.), -1.))
    if half_theta != 0.
        theta = 2. * half_theta
        r = q[1:3] / sin(half_theta)
        r_mag = sqrt(sum(r.^2))
        if r_mag == 0.
            r = [1.; 0.; 0.]
        else
            r = r ./ r_mag
        end
    else
        theta = 0.
        r = [1., 0., 0.]
    end
    return (theta, r)
end

function qinv(q::Vector{T}) where {T<:Real}
    @assert(length(q) == 4, "Input must be a 4-element quaternion.")
    return [-q[1:3]; q[4]]
end

function qpos(q::Vector{T}) where {T<:Real}
    @assert(length(q) == 4, "Input must be a 4-element quaternion.")
    return q[4] < zero(T) ? -q : q
end

function qcomp(a::Vector{T}, b::Vector{T}) where {T<:Real}
    @assert(length(a) == 4, "First input must be a 4-element quaternion.")
    @assert(length(b) == 4, "Second input must be a 4-element quaternion.")
    return [ a[4]*b[1] + a[3]*b[2] - a[2]*b[3] + a[1]*b[4];
            -a[3]*b[1] + a[4]*b[2] + a[1]*b[3] + a[2]*b[4];
             a[2]*b[1] - a[1]*b[2] + a[4]*b[3] + a[3]*b[4];
            -a[1]*b[1] - a[2]*b[2] - a[3]*b[3] + a[4]*b[4]]

end

function qdiff(a::Vector{T}, b::Vector{T}) where {T<:Real}
    @assert(length(a) == 4, "First input must be a 4-element quaternion.")
    @assert(length(b) == 4, "Second input must be a 4-element quaternion.")
    return [-a[4]*b[1] - a[3]*b[2] + a[2]*b[3] + a[1]*b[4];
             a[3]*b[1] - a[4]*b[2] - a[1]*b[3] + a[2]*b[4];
            -a[2]*b[1] + a[1]*b[2] - a[4]*b[3] + a[3]*b[4];
             a[1]*b[1] + a[2]*b[2] + a[3]*b[3] + a[4]*b[4]]
end

function qdot(q::Vector{T}, w::Vector{T}, k=zero(T)) where {T<:Real}
    @assert(length(q) == 4, "First input must be a 4-element quaternion.")
    @assert(length(w) == 3, "Second input must be a 3-element rotation rate.")
    qd = T(0.5) * [               w[3] * q[2] - w[2] * q[3] + w[1] * q[4];
                   -w[3] * q[1]               + w[1] * q[3] + w[2] * q[4];
                    w[2] * q[1] - w[1] * q[2]               + w[3] * q[4];
                   -w[1] * q[1] - w[2] * q[2] - w[3] * q[3]               ]
    if k != zero(T)
        m2 = q[1]^2 + q[2]^2 + q[3]^2 + q[4]^2
        s = k * (one(T) - m2)
        qd = qd + s * q
    end
    return qd
end

function mrp2q(p::Vector{T}, scale=one(T)) where {T<:Real}
    f2  = scale^2
    pm2 = p[1]^2 + p[2]^2 + p[3]^2
    c0  = one(T) / (f2 + pm2)
    c1  = 2scale * c0
    return [c1 * p[1];
            c1 * p[2];
            c1 * p[3];
            c0 * (f2 - pm2)]
end

end


text:
function inner_product(x::AbstractExpr, y::AbstractExpr)
    if x.size == y.size && x.size[1] == x.size[2]
        return real(tr(x'*y))
    else
        error("Arguments must be square matrix of same dimension")
    end
end

inner_product(x::Value, y::AbstractExpr) = inner_product(Constant(x),y)
inner_product(x::AbstractExpr, y::Value) = inner_product(x,Constant(y))


text:
##########  Configuration   #########

have_winston = false
# # Check whether Winston is installed
# macro usingif(status, sym)
#     if Pkg.installed(string(sym)) != nothing
#         return expr(:toplevel, {expr(:using, {sym}), :($(esc(status)) = true)})
#     end
#     :($(esc(status)) = false)
# end
#
# @usingif have_winston Winston

# Find a system image viewer
imshow_cmd = ""
@unix_only begin
if !have_winston
    imshow_cmd_list = ["feh", "gwenview", "open"]
    for thiscmd in imshow_cmd_list
        _, p = readsfrom(`which $thiscmd`)
        wait(p)
        if success(p)
            imshow_cmd = thiscmd
            break
        end
    end
end
end


text:
using Memoize, Test

@test_throws LoadError eval(:(@memoize))
@test_throws LoadError eval(:(@memoize () = ()))

# you can't use test_throws in macros
arun = 0
@memoize function memadd(x::Int, y::Int)::Int
    global arun += 1
    print("evaluating memadd $x $y\n")
    return x + y
end
@test memadd(1, 2) == 3
@test memadd(1, 4) == 5
@test memadd(1, 2) == 3
@test arun == 2

run = 0
@memoize function simple(a)
    global run += 1
    a
end
@test simple(5) == 5
@test run == 1
@test simple(5) == 5
@test run == 1
@test simple(6) == 6
@test run == 2
@test simple(6) == 6
@test run == 2

empty!(memoize_cache(simple))
@test simple(6) == 6
@test run == 3
@test simple(6) == 6
@test run == 3

run = 0
@memoize function typed(a::Int)
    global run += 1
    a
end
@test typed(5) == 5
@test run == 1
@test typed(5) == 5
@test run == 1
@test typed(6) == 6
@test run == 2
@test typed(6) == 6
@test run == 2

run = 0
@memoize function default(a = 2)
    global run += 1
    a
end
@test default() == 2
@test run == 1
@test default() == 2
@test run == 1
@test default(2) == 2
@test run == 1
@test default(6) == 6
@test run == 2
@test default(6) == 6
@test run == 2

run = 0
@memoize function default_typed(a::Int = 2)
    global run += 1
    a
end
@test default_typed() == 2
@test run == 1
@test default_typed() == 2
@test run == 1
@test default_typed(2) == 2
@test run == 1
@test default_typed(6) == 6
@test run == 2
@test default_typed(6) == 6
@test run == 2

run = 0
@memoize function kw(; a = 2)
    global run += 1
    a
end
@test kw() == 2
@test run == 1
@test kw() == 2
@test run == 1
@test kw(a = 2) == 2
@test run == 1
@test kw(a = 6) == 6
@test run == 2
@test kw(a = 6) == 6
@test run == 2

run = 0
@memoize function kw_typed(; a::Int = 2)
    global run += 1
    a
end
@test kw_typed() == 2
@test run == 1
@test kw_typed() == 2
@test run == 1
@test kw_typed(a = 2) == 2
@test run == 1
@test kw_typed(a = 6) == 6
@test run == 2
@test kw_typed(a = 6) == 6
@test run == 2

run = 0
@memoize function default_kw(a = 1; b = 2)
    global run += 1
    (a, b)
end
@test default_kw() == (1, 2)
@test run == 1
@test default_kw() == (1, 2)
@test run == 1
@test default_kw(1, b = 2) == (1, 2)
@test run == 1
@test default_kw(1, b = 3) == (1, 3)
@test run == 2
@test default_kw(1, b = 3) == (1, 3)
@test run == 2
@test default_kw(2, b = 3) == (2, 3)
@test run == 3
@test default_kw(2, b = 3) == (2, 3)
@test run == 3

run = 0
@memoize function default_kw_typed(a::Int = 1; b::Int = 2)
    global run += 1
    (a, b)
end
@test default_kw_typed() == (1, 2)
@test run == 1
@test default_kw_typed() == (1, 2)
@test run == 1
@test default_kw_typed(1, b = 2) == (1, 2)
@test run == 1
@test default_kw_typed(1, b = 3) == (1, 3)
@test run == 2
@test default_kw_typed(1, b = 3) == (1, 3)
@test run == 2
@test default_kw_typed(2, b = 3) == (2, 3)
@test run == 3
@test default_kw_typed(2, b = 3) == (2, 3)
@test run == 3

run = 0
@memoize function required_kw(; a)
    global run += 1
    a
end
@test required_kw(a = 1) == 1
@test run == 1
@test required_kw(a = 1) == 1
@test run == 1
@test required_kw(a = 2) == 2
@test run == 2
@test required_kw(a = 2) == 2
@test run == 2
@test_throws UndefKeywordError required_kw()

run = 0
@memoize function ellipsis(a, b...)
    global run += 1
    (a, b)
end
@test ellipsis(1) == (1, ())
@test run == 1
@test ellipsis(1) == (1, ())
@test run == 1
@test ellipsis(1, 2, 3) == (1, (2, 3))
@test run == 2
@test ellipsis(1, 2, 3) == (1, (2, 3))
@test run == 2

run = 0
@memoize Dict function kw_ellipsis(;a...)
    global run += 1
    a
end
@test isempty(kw_ellipsis())
@test run == 1
@test isempty(kw_ellipsis())
@test run == 1
@test kw_ellipsis(a = 1) == pairs((a = 1,))
@test run == 2
@test kw_ellipsis(a = 1) == pairs((a = 1,))
@test run == 2
@test kw_ellipsis(a = 1, b = 2) == pairs((a = 1, b = 2))
@test run == 3
@test kw_ellipsis(a = 1, b = 2) == pairs((a = 1, b = 2))
@test run == 3

run = 0
@memoize function multiple_dispatch(a::Int)
    global run += 1
    1
end
@memoize function multiple_dispatch(a::Float64)
    global run += 1
    2
end
@test multiple_dispatch(1) == 1
@test run == 1
@test multiple_dispatch(1) == 1
@test run == 1
@test multiple_dispatch(1.0) == 2
@test run == 2
@test multiple_dispatch(1.0) == 2
@test run == 2

run = 0
@memoize function where_clause(a::T) where T
    global run += 1
    T
end
@test where_clause(1) == Int
@test run == 1
@test where_clause(1) == Int
@test run == 1

function outer()
    run = 0
    @memoize function inner(x)
        run += 1
        x
    end
    @memoize function inner(x, y)
        run += 1
        x + y
    end
    @test inner(5) == 5
    @test run == 1
    @test inner(5) == 5
    @test run == 1
    @test inner(6) == 6
    @test run == 2
    @test inner(6) == 6
    @test run == 2
    @test inner(5, 1) == 6
    @test run == 3
    @test inner(5, 1) == 6
    @test run == 3
end
outer()
@test !@isdefined inner

@memoize function typeinf(x)
    x + 1
end
@inferred typeinf(1)
@inferred typeinf(1.0)

finalized = false
@memoize function method_rewrite()
    x = []
    finalizer(x->(global finalized; finalized = true), x)
    x
end
method_rewrite()
@memoize function method_rewrite() end
GC.gc()
@test finalized

run = 0
""" documented function """
@memoize function documented_function(a)
    global run += 1
    a
end
@test strip(string(@doc documented_function)) == "documented function"
@test documented_function(1) == 1
@test run == 1
@test documented_function(1) == 1
@test run == 1
@test documented_function(2) == 2
@test run == 2
@test documented_function(2) == 2
@test run == 2

run = 0
@memoize function vararg_func(list::Vararg{Tuple{Int64,Int64}})
    global run += 1
    return list[1]
end
@test vararg_func((1,1), (1,1)) == (1,1)
@test run == 1
@test vararg_func((1,1), (1,1)) == (1,1)
@test run == 1
@test vararg_func((1,1), (1,2)) == (1,1)
@test run == 2
@test vararg_func((1,1), (1,2)) == (1,1)
@test run == 2

module MemoizeTest
using Test
using Memoize

const MyDict = Dict

run = 0
@memoize MyDict function custom_dict(a)
    global run += 1
    a
end
@test custom_dict(1) == 1
@test run == 1
@test custom_dict(1) == 1
@test run == 1
@test custom_dict(2) == 2
@test run == 2
@test custom_dict(2) == 2
@test run == 2
@test custom_dict(1) == 1
@test run == 2

end # module

using .MemoizeTest
using .MemoizeTest: custom_dict

empty!(memoize_cache(custom_dict))
@test custom_dict(1) == 1
@test MemoizeTest.run == 3
@test custom_dict(1) == 1
@test MemoizeTest.run == 3

empty!(memoize_cache(MemoizeTest.custom_dict))
@test custom_dict(1) == 1
@test MemoizeTest.run == 4

run = 0
@memoize Dict{Tuple{String},Int}() function dict_call(a::String)::Int
    global run += 1
    length(a)
end
@test dict_call("a") == 1
@test run == 1
@test dict_call("a") == 1
@test run == 1
@test dict_call("bb") == 2
@test run == 2
@test dict_call("bb") == 2
@test run == 2


text:
using HSA
include("helpers.jl")
using FactCheck

@with_codegen facts("The Intrinsics generation logic") do
    context("can parse a mangled name") do
        defi32 = HSA.Builtins.parse_mangled_name("@_Z13get_global_idj")

        @fact defi32.name --> "get_global_id"
        @fact defi32.arg_types --> [Int32]

        defv = HSA.Builtins.parse_mangled_name("@_Z20get_global_linear_idv")

        @fact defv.name --> "get_global_linear_id"
        @fact defv.arg_types --> []
    end
end