text:
#################### MCMC Simulation Engine ####################

"""
    mcmc(mc::ModelChains, iters::Integer; verbose::Bool=true, trees::Bool=false)

This function simulates additional draws from a model.

* `mc` is the results of a previous call to the mcmc function.

* `iters` indicates the number of draws to simulate.

* `verbose` controls whether to print porgress statements to the console.

* `trees` indicates if the states of the model nodes describing tree structures should be stored as well.
"""
function mcmc(mc::ModelChains, iters::Integer; verbose::Bool=true, trees::Bool=false)
  thin = step(mc)
  last(mc) == div(mc.model.iter, thin) * thin ||
    throw(ArgumentError("chain is missing its last iteration"))

  mm = deepcopy(mc.model)
  mc2 = mcmc_master!(mm, mm.iter .+ (1:iters), last(mc), thin, mc.chains,
                     verbose, trees)
  if mc2.names != mc.names
    mc2 = mc2[:, mc.names, :]
  end

  ModelChains(vcat(mc, mc2), mc2.model)
end

"""
    mcmc(m::Model, inputs::Dict{Symbol},
              inits::Vector{V} where V<:Dict{Symbol},
              iters::Integer; burnin::Integer=0, thin::Integer=1,
              chains::Integer=1, verbose::Bool=true, trees::Bool=false)

Simulate MCMC draws from the model `m`.

* `inputs` is a dictionary storing the values for input model nodes. Dictionary keys and values should be given for each input node.

* `inits` contains dictionaries with initial values for stochastic model nodes. Dictionary keys and values should be given for each stochastic node. Consecutive runs of the simulator will iterate through the vector’s dictionary elements.

* `iter` Specifies the number of draws to generate for each simulation run

* `burnin` specifies the number of initial draws to discard as a burn-in sequence to allow for convergence.

* `thin` is the step-size between draws to store in the output.

* `chains` specifies the number of simulation runs to perform.

* `verbose` indicates whether the sampler progress should be printed at the console.

* `trees` indicates if the states of the model nodes describing tree structures should be stored as well.
"""
function mcmc(m::Model, inputs::Dict{Symbol},
              inits::Vector{V} where V<:Dict{Symbol},
              iters::Integer; burnin::Integer=0, thin::Integer=1,
              chains::Integer=1, verbose::Bool=true, trees::Bool=false)
  iters > burnin ||
    throw(ArgumentError("burnin is greater than or equal to iters"))
  length(inits) >= chains ||
    throw(ArgumentError("fewer initial values than chains"))

  mm::Model = deepcopy(m)
  setinputs!(mm, inputs)
  setinits!(mm, inits[1:chains])
  mm.burnin = burnin
  mcmc_master!(mm, 1:iters, burnin, thin, 1:chains, verbose, trees)
end

function mcmc_master!(m::Model, window::UnitRange{Int}, burnin::Integer,
                      thin::Integer, chains::AbstractArray{Int}, verbose::Bool, trees::Bool)
  states::Vector{ModelState} = m.states
  m.states = ModelState[]

  N = length(window)
  K = length(chains)

  frame::ChainProgressFrame = ChainProgressFrame(
    "MCMC Simulation of $N Iterations x $K Chain" * "s"^(K > 1), verbose
  )

  lsts = [
    Any[m, states[k], window, burnin, thin, ChainProgress(frame, k, N), trees]
    for k in chains
  ]
  results::Vector{Tuple{Chains, Model, ModelState}} = pmap2(mcmc_worker!, lsts)

  sims::Array{Chains}  = Chains[results[k][1] for k in 1:K]
  model::Model = results[1][2]
  model.states = ModelState[results[k][3] for k in sortperm(chains)]

  ModelChains(cat(sims..., dims=3), model)
end

function mcmc_worker!(args::Vector)::Tuple{Chains, Model, ModelState}
  m::Model, state::ModelState, window::UnitRange{Int}, burnin::Integer, thin::Integer, meter::ChainProgress, store_trees::Bool = args
  llname::AbstractString = "likelihood"
  treeind = 1
  simind = 1
  m.iter = first(window) - 1

  relist!(m, state.value)

  settune!(m, state.tune)
  pnames = vcat(names(m, true), llname)
  treenodes = Symbol[]
  for i in m.nodes
    if isa(i[2], TreeStochastic)
      push!(treenodes ,i[1])
    end
  end

  sim = Chains(last(window), length(pnames), start=burnin + thin, thin=thin,
               names=pnames, ntrees=length(treenodes), tree_names=treenodes)

  reset!(meter)
  for i in window

    sample!(m)
    if i > burnin && (i - burnin) % thin == 0

      sim[i, :, 1] = unlist(m, true)

      if store_trees
       for (ind, tree_node) in enumerate(treenodes)
         sim.trees[treeind, ind, 1] = newick(m[tree_node].value)
       end
       treeind +=1
     end
    end
    next!(meter)
  end

  mv = samparas(m)

  sim.moves[1] = mv
  (sim, m, ModelState(unlist(m), gettune(m)))
end

function track(i::Integer, burnin::Integer, thin::Integer,
               sim::Chains, m::Model, store_trees::Bool, treeind::Integer, treenode)
  if i > burnin && (i - burnin) % thin == 0
    sim[i, :, 1] = unlist(m, true)
    if store_trees
     sim.trees[treeind, 1, 1] = newick(m[treenode].value)
     treeind +=1
   end
  end
end


text:
function print_debug_Z(i::Int, n::Int, j::Int, s::State, c::Constants, d::Data)
  println("y_inj = $(s.y_imputed[i][n, j])")
  println("sig2[i]: $(s.sig2[i])")
  println("mus0: $(mus(false, s, c, d))")
  println("mus1: $(mus(true, s, c, d))")
  println("etaz0_i$(i)_j$(j): $(s.eta[0][i, j, :])")
  println("etaz1_i$(i)_j$(j): $(s.eta[1][i, j, :])")
end

function update_Z!(s::State, c::Constants, d::Data, sb_ibp::Bool;
                   use_repulsive::Bool=false)
  ll0 = zeros(d.J, c.K)
  ll1 = zeros(d.J, c.K)

  for i in 1:d.I
    for n in 1:d.N[i]
      k = s.lam[i][n]
      if k > 0
        for j in 1:d.J
          ll0[j, k] += logdmixture(false, i, n, j, s, c, d)
          ll1[j, k] += logdmixture(true, i, n, j, s, c, d)
        end
      end
    end
  end

  for k in 1:c.K
    v = sb_ibp ? cumprod(s.v) : s.v
    lp0 = log1p(-v[k])
    lp1 = log(v[k])

    for j in 1:d.J
      lfc0 = lp0 + ll0[j, k]
      lfc1 = lp1 + ll1[j, k]

      if use_repulsive
        # Make Z matrix
        Z0 = copy(s.Z)
        Z0[j, k] = false

        Z1 = copy(s.Z)
        Z1[j, k] = true

        # Add penalty terms
        # NOTE: log_penalty_repFAM = log repulsive term
        lfc0 += log_penalty_repFAM(k, Z0, c.log_repulsive_fn)
        lfc1 += log_penalty_repFAM(k, Z1, c.log_repulsive_fn)
      end

      p = 1.0 / (1.0 + exp(lfc0 - lfc1))
      if isnan(p)
        println("WARNING in update_Z: p = NaN.")
      end

      s.Z[j, k] = p > rand() # rand(Bernoulli(p))
    end
  end
end

ext: jl


text:
using StatsBase, Statistics, Dates

# setup
nworkers = 2
windowsize = Second(10)
qlower, qupper = 0.0, 1.0
chin, chout = CodedComputing.setup_profiler_channels()
θs = [0.3, 0.7]
qs = 1 ./ [2, 3]

# put example values into the input channel
## worker 1
timestamps = Time(0) .+ [Second(i) for i in 1:11]
comps = rand(length(timestamps))
comms = rand(length(timestamps))
for i in 1:length(timestamps)
    v = CodedComputing.ProfilerInput(1, θs[1], qs[1], timestamps[i], comps[i], comms[i])
    push!(chin, v)
end
comp_mean = mean(comps) / (θs[1] * qs[1])
comp_var = var(comps ./ (θs[1] * qs[1]), corrected=false)
comm_mean, comm_var = mean(comms), var(comms, corrected=false)
correct1 = CodedComputing.ProfilerOutput(1, θs[1], qs[1], comp_mean, comp_var, comm_mean, comm_var)

## worker 2
timestamps = Time(0) .+ [Second(i) for i in 1:20]
comps = rand(length(timestamps))
comms = rand(length(timestamps))
for i in 1:length(timestamps)
    v = CodedComputing.ProfilerInput(2, θs[2], qs[2], timestamps[i], comps[i], comms[i])
    push!(chin, v)
end
comps = comps[end-10:end]
comms = comms[end-10:end]
comp_mean = mean(comps) / (θs[2] * qs[2])
comp_var = var(comps ./ (θs[2] * qs[2]), corrected=false)
comm_mean, comm_var = mean(comms), var(comms, corrected=false)
correct2 = CodedComputing.ProfilerOutput(2, θs[2], qs[2], comp_mean, comp_var, comm_mean, comm_var)

# start the task and test that it consumes the input values
task = Threads.@spawn CodedComputing.latency_profiler(chin, chout; nworkers, qlower, qupper, windowsize, minsamples=1)

# wait for up to 10 seconds for the input to be consumed
t0 = time_ns()
while (time_ns() - t0)/1e9 < 10 && isready(chin)
    sleep(0.1)
end
@test !isready(chin)

# wait for up to 10 seconds for the subsystem to produce output
t0 = time_ns()
while (time_ns() - t0)/1e9 < 10 && !isready(chout)
    sleep(0.1)
end

# stop the profiler
close(chin)
wait(task)

# there should be 2 values in the output channel
## test the first output value
@test isready(chout)
v = take!(chout)
@test v.worker == 1 || v.worker == 2
correct = v.worker == 1 ? correct1 : correct2
for name in fieldnames(CodedComputing.ProfilerOutput)
    @test getfield(correct, name) ≈ getfield(v, name)
end

# wait for up to 10 seconds for the subsystem to produce output
t0 = time_ns()
while (time_ns() - t0)/1e9 < 10 && !isready(chout)
    sleep(0.1)
end

## test the second output value
@test isready(chout)
v = take!(chout)
correct = v.worker == 1 ? correct1 : correct2
for name in fieldnames(CodedComputing.ProfilerOutput)
    @test getfield(correct, name) ≈ getfield(v, name)
end

close(chout)


text:
"""
This module contains:

  - Data types:
    + `SparseMatrixCSR`:Compressed Sparse Row (CSR) sparse matrix implementation with Bi-based indexing.
    + `SymSparseMatrixCSR`: Symmetric Compressed Sparse Row sparse matrix implementation with Bi-based indexing.

  - Procedures:
    + `push_coo!`: Helper function to build COO arrays for further building a SparseMatrix
    + `finalize_coo!`: Finalization of COO arrays building.
    + `sparsecsr`: Analogous method to [`sparse`](@ref) function to build a SparseMatrixCSR.
    + `symsparsecsr`: Analogous method to [`sparse`](@ref) function to build a SymSparseMatrixCSR.
    + `colvals`: Analogous method to [`rowvals`](@ref) to return `colvals` array.
    + `hasrowmajororder`: Return `true` if matrix values are ordered by row.
    + `hascolmajororder`: Return `true` if matrix values are ordered by column.
    + `getptr`: Return the pointer array of a SparseMatrix (`rowptr` or `colptr` depending on the SparseMatrix type)
    + `getindices`: Return the indices array of a SparseMatrix (`rowval` or `colval` depending on the SparseMatrix type)

  - Overloaded procedures:
    + `*`: SparseMatrix-Vector product.
    + `mul!`: SparseMatrix-Vector product.
    + `nnz`: Return the number of stored (filled) elements in a sparse array.
    + `nonzeros`: Return `nzval` array.
    + `nzrange`: Return the range of indices for a particular row or column (Depending on the SparseMatrix type)
    + `findnz`: Return a tuple (I, J, V) where I and J are the row and column indices.
    + `rowvals`: Return row indices or raises an error (Depending on the SparseMatrix type)
    + `convert`: Type conversion

"""
module SparseMatricesCSR

using SparseArrays
using LinearAlgebra

import Base: convert, size, getindex, show, count, *
import LinearAlgebra: mul!
import SparseArrays: nnz, getnzval, nonzeros, nzrange, findnz, rowvals

export SparseMatrixCSR
export SymSparseMatrixCSR
export push_coo!, finalize_coo!, sparsecsr, symsparsecsr, colvals
export hasrowmajororder, hascolmajororder, getptr, getindices

include("SparseMatrixCSC.jl")
include("SparseMatrixCSR.jl")
include("SymSparseMatrixCSR.jl")

end # module


text:
"""
   update_input(state, agent_idx, config)

Update the incoming directed edges from other agents.

# Arguments
- `state`: a tuple of the current graph and agent_list
- `agent_idx`: agent index
- `config`: Config object as provided by Config

See also: [`Config`](@ref)
"""
function update_input!(
    state::Tuple{AbstractGraph, AbstractArray},
    agent_idx::Integer,
    config::Config
)

    graph, agent_list = state
    this_agent = agent_list[agent_idx]

    pref_attach_list = append!(
        [dst(e) for e in edges(graph)],
        1:nv(graph)
    )

    exclude = union(inneighbors(graph, agent_idx), agent_idx)

    function draw_new_input(pref_attach_list, exclude)
        draw = rand(pref_attach_list)
        if draw in exclude
            draw = draw_new_input(pref_attach_list, exclude)
        end
        return draw
    end

    if (length(exclude) >= (nv(graph) - config.network.new_follows))
        new_input_count = nv(graph) - length(exclude) - 1
    else
        new_input_count = config.network.new_follows
    end

    for _ in 1:new_input_count
        new_neighbor = draw_new_input(pref_attach_list, exclude)
        add_edge!(graph, new_neighbor, agent_idx)
        push!(exclude, new_neighbor)
    end

    unfollow_candidates = Array{Tuple{Int64, Float64}, 1}()

    for neighbor in inneighbors(graph, agent_idx)
        push!(
            unfollow_candidates,
            (neighbor, abs(this_agent.opinion - agent_list[neighbor].opinion))
        )
    end

    sort!(unfollow_candidates, by=last)

    for i in 1:config.network.new_follows
        rem_edge!(graph, popfirst!(unfollow_candidates)[1], agent_idx)
    end

    return state

end


text:
function importSR(SIMurl::String,LOCurl::String)
#1.. Import Data
	#Load Data with colnames:
	#mol_ID|loc_ID|framenumber|pos_x|pos_y
	#[Int64,Int64,Int64,Float64,Float64,Float64,Float64,Float64,Float64,Float64]
	#DataSim = CSV.read("/Users/Patrick/Documents/GitHub/BayesianEstimation/STORM/Simulationen_01-24-19/sim_Tet.csv",
	#DataSim = CSV.read("/Users/Patrick/Documents/GitHub/BayesianEstimation/STORM/ORI_SIM/locs_gamma_conf.csv",
	# print("importing DataSim\n")
	# DataSim = CSV.read(SIMurl,
	# 				types=[Int64,Int64,Int64,Float64,Float64],
	# 				header=["label_ID","oligomer_ID","framenum","x_nm","y_nm"], datarow=2, DataFrame)
	print("importing DataSim\n")
	DataSim = CSV.read(SIMurl,
					types=[Int64,Int64,Float64,Float64,Int64],
					header=["label_ID","framenum","x_nm","y_nm","oligomer_ID"], datarow=2, DataFrame)
	#RealLoc = CSV.read("/Users/Patrick/Documents/GitHub/BayesianEstimation/STORM/Simulationen_01-24-19/sim_Tet_locs.csv",
	#RealLoc = CSV.read("/Users/Patrick/Documents/GitHub/BayesianEstimation/STORM/ORI_SIM/labels_gamma.csv",
	print("importing RealLoc\n")
	RealLoc = CSV.read(LOCurl,
					types=[Int64,Float64,Float64,Int64],
					#header=["label_ID","x_nm","y_nm"], datarow=2, DataFrame)
					header=["label_ID","x_nm","y_nm","oligomer_ID"], datarow=2, DataFrame)
	return DataSim, RealLoc
end

export importSR

max_stars_repo_name: PGSch/SuperResolution

max_issues_repo_name: PGSch/SuperResolution

max_forks_repo_name: PGSch/SuperResolution


text:
using Vectorize

## Use correct version of Base.Test
if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

## Initialize testing globals
println("===== Testing Accelerate =====")
N = 1000
srand(13)

## LOGARITHM
for T in (Float32, Float64)
    @testset "Accelerate: Logarithmic::$T" begin
        X::Array{T} = exp(10*randn(N))
        @testset "Testing $f::$T" for f in [:log,:log2,:log10, :log1p]
            @eval fb = $f
            @eval fa = Vectorize.Accelerate.$f
            @test fa(X) ≈ fb(X)
        end
    end
end

## EXPONENIATION
for T in (Float32, Float64)
    @testset "Accelerate: Exponential::$T" begin
        @testset "Testing $f::$T" for f in [:exp,:exp2,:expm1]
            X::Array{T} = 10*randn(N)
            @eval fb = $f
            @eval fa = Vectorize.Accelerate.$f
            @test fa(X) ≈ fb(X)
        end
    end
end

## TRIGONOMETRIC
for T in (Float32, Float64)
    X::Array{T} = 10*randn(N)
    @testset "Accelerate: Trigonometric::$T" begin
        @testset "Testing $f::$T" for f in [:sin,:sinpi,:cos,:cospi,:tan,:atan] # tanpi not defined in Base
            @eval fb = $f
            @eval fa = Vectorize.Accelerate.$f
            @test fa(X) ≈ fb(X)
        end

        Y::Array{T} = 10*randn(N)
        @testset "Testing $f::$T" for f in [:atan2]
            @eval fb = $f
            @eval fa = Vectorize.Accelerate.$f
            @test fa(X,Y) ≈ fb(X,Y)
        end

        Z::Array{T} = 2*rand(N)-1
        @testset "Testing $f::$T" for f in [:asin,:acos]
            @eval fb = $f
            @eval fa = Vectorize.Accelerate.$f
            @test fa(Z) ≈ fb(Z)
        end
    end
end

## HYPERBOLIC
for T in (Float32, Float64)
    @testset "Accelerate: Hyperbolic::$T" begin
        X = 10*randn(N)
        @testset "Testing $f::$T" for f in [:sinh,:cosh,:tanh,:asinh]
            @eval fb = $f
            @eval fa = Vectorize.Accelerate.$f
            @test fa(X) ≈ fb(X)
        end

        Y = exp(10*randn(N))+1
        @testset "Testing $f::$T" for f in [:acosh]
            @eval fb = $f
            @eval fa = Vectorize.Accelerate.$f
            @test fa(Y) ≈ fb(Y)
        end

        Z = 2*rand(N)-1
        @testset "Testing $f::$T" for f in [:atanh]
            @eval fb = $f
            @eval fa = Vectorize.Accelerate.$f
            @test fa(Z) ≈ fb(Z)
        end
    end
end

## VECTOR OPERATIONS
for T in (Float32, Float64)
    @testset "Accelerate: Vector Operations::$T" begin
        X = 10*randn(N)
        Y = 10*randn(N)
        @testset "Testing $fa::$T" for (f, fa) in [(:+, :add), (:-, :sub), (:.*, :mul), (:./, :div)]
            @eval fb = $f
            @eval facc = Vectorize.Accelerate.$fa
            @test fb(X, Y) ≈ facc(X, Y)
        end
    end
end

println("===== Accelerate Tests Successful =====\n\n")


text:
# Displaying color swatches
# -------------------------

@compat function Base.show(io::IO, ::MIME"image/svg+xml", c::Color)
    write(io,
        """
        <?xml version"1.0" encoding="UTF-8"?>
        <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
         "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
        <svg xmlns="http://www.w3.org/2000/svg" version="1.1"
             width="25mm" height="25mm" viewBox="0 0 1 1">
             <rect width="1" height="1"
                   fill="#$(hex(c))" stroke="none"/>
        </svg>
        """)
end

@compat function Base.show{T <: Color}(io::IO, ::MIME"image/svg+xml",
                                       cs::AbstractVecOrMat{T})
    m,n = ndims(cs) == 2 ? size(cs) : (1,length(cs))

    xsize,xpad = n > 50 ? (250/n,0) : n > 18 ? (5.,1) : n > 12 ? (10.,1) : n > 1 ? (15.,1) : (25.,0)
    ysize,ypad = m > 28 ? (150/m,0) : m > 14 ? (5.,1) : m > 9 ? (10.,1) : m > 1 ? (15.,1) : (25.,0)

    write(io,
        """
        <?xml version="1.0" encoding="UTF-8"?>
        <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
         "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
        <svg xmlns="http://www.w3.org/2000/svg" version="1.1"
             width="$(n*xsize)mm" height="$(m*ysize)mm"
             shape-rendering="crispEdges">
        """)

    for i in 1:m, j in 1:n
        c = ndims(cs) == 2 ? cs[i, j] : cs[j]
        write(io,
            """
            <rect x="$((j-1)*xsize)mm" y="$((i-1)*ysize)mm"
                  width="$(xsize - xpad)mm" height="$(ysize - ypad)mm"
                  fill="#$(hex(c))" stroke="none" />
            """)
    end

    write(io, "</svg>")
end


text:
#=================================
Methods using Taylor Models
=================================#
using .TaylorModels

@inline zeroBox(N) = IntervalBox(0..0, N)
@inline symBox(N) = IntervalBox(-1..1, N)

function _enclose(tme::TaylorModelsEnclosure, f::Function, dom::Interval_or_IntervalBox;
                  kwargs...)

    if tme.normalize
        R = _enclose_TaylorModels_norm(f, dom, tme.order)
    else
        R = _enclose_TaylorModels(f, dom, tme.order)
    end
    return R
end

# univariate
function _enclose_TaylorModels(f::Function, dom::Interval, order::Int)
    x0 = Interval(mid(dom))
    x = TaylorModel1(order, x0, dom)
    return evaluate(f(x-x0), dom)
end

# normalized univariate
function _enclose_TaylorModels_norm(f::Function, dom::Interval, order::Int)
    x0 = Interval(mid(dom))
    x = TaylorModel1(order, x0, dom)
    xnorm = normalize_taylor(x.pol, dom - x0, true)
    xnormTM = TaylorModel1(xnorm, 0..0, 0..0, -1..1)
    return evaluate(f(xnormTM), -1..1)
end

# multivariate
function _enclose_TaylorModels(f::Function, dom::IntervalBox{N}, order::Int) where {N}
    x0 = mid(dom)
    set_variables(Float64, "x", order=2order, numvars=N)
    x = [TaylorModelN(i, order, IntervalBox(x0), dom) for i=1:N]
    return evaluate(f(x...), dom - x0)
end

# normalized multivariate
function _enclose_TaylorModels_norm(f::Function, dom::IntervalBox{N}, order::Int) where {N}
    x0 = mid(dom)
    set_variables(Float64, "x", order=2order, numvars=N)

    zBoxN = zeroBox(N)
    sBoxN = symBox(N)
    x = [TaylorModelN(i, order, IntervalBox(x0), dom) for i=1:N]
    xnorm = [normalize_taylor(xi.pol, dom - x0, true) for xi in x]
    xnormTM = [TaylorModelN(xi_norm, 0..0, zBoxN, sBoxN) for xi_norm in xnorm]
    return evaluate(f(xnormTM...), sBoxN)
end


text:
module AmbulanceDeployment

    import JuMP, Gurobi
    import Base.Collections: PriorityQueue, enqueue!, dequeue!
    import DataFrames: DataFrame, isna, nrow
    import Distributions: Poisson, LogNormal, quantile, sample

    export DeploymentProblem,
           RobustDeployment,
           StochasticDeployment,
           MALPDeployment,
           MEXCLPDeployment,

           DispatchProblem,
           ClosestDispatch,

           NoRedeployModel,
           AssignmentModel,

           solve,
           evaluate,
           deployment,
           convergence_plot,
           compose_neighborhoods,
           compose_locations,
           compose_chloropleth,
           performance,
           test_performance,
           plot_timings,
           simulate_events!,
           initialize!

    include("problem.jl")
    include("model.jl")
    include("simulate.jl")

    include("evaluate.jl")
    include("plot.jl")
end


text:
module Day14

export solvePart1, solvePart2

#= input = 633601 =#
input = 633601

function solvePart1()
    scores = [3, 7]
    elf1, elf2 = 1, 2

    while length(scores) < input + 10
        toAdd = scores[elf1] + scores[elf2]
        if toAdd <= 9
            push!(scores, toAdd)
        else
            push!(scores, div(toAdd, 10))
            push!(scores, toAdd % 10)
        end
        n = length(scores)
        elf1 += 1 + scores[elf1]
        while elf1 > n
            elf1 -= n
        end
        elf2 += 1 + scores[elf2]
        while elf2 > n
            elf2 -= n
        end
    end
    strings = map(score -> string(score), scores[input+1:input+10])
    return join(strings, "")
end

# A local variable / parameter is much cheaper to access than a module-global.
function part2Inner(target)
    scores = [3, 7]
    elf1, elf2 = 1, 2

    digits = Int(floor(log10(target))) + 1

    found = false
    while !found
        toAdd = scores[elf1] + scores[elf2]
        if toAdd <= 9
            push!(scores, toAdd)
        else
            push!(scores, div(toAdd, 10))
            push!(scores, toAdd % 10)
        end
        n = length(scores)
        elf1 += 1 + scores[elf1]
        while elf1 > n
            elf1 -= n
        end
        elf2 += 1 + scores[elf2]
        while elf2 > n
            elf2 -= n
        end

        if n >= digits
            final = 0
            semiFinal = 0
            for i in -digits+1:0
                final = final * 10 + scores[n+i]
            end
            if n >= digits + 1
                for i in -digits:-1
                    semiFinal = semiFinal * 10 + scores[n+i]
                end
            end

            if final == target
                return n - digits
            elseif semiFinal == target
                return n - (digits+1)
            end

        end
    end
end

function solvePart2()
    return part2Inner(input)
end
end


text:
module KnetMetrics

import LinearAlgebra

import Statistics

include("utilities/utilities.jl")
include("classification/classification.jl"); export confusion_matrix, class_confusion, visualize, classification_report, condition_positive, condition_negative, predicted_positive, predicted_negative, correctly_classified, incorrectly_classified, sensitivity_score, true_positive_rate, recall_score, specificity_score, true_negative_rate, precision_score, positive_predictive_value, accuracy_score, balanced_accuracy_score, negative_predictive_value, false_negative_rate, false_positive_rate, false_discovery_rate, false_omission_rate, fbeta_score, f1_score, prevalence_threshold, threat_score, matthews_correlation_coeff, fowlkes_mallows_index, informedness, markedness, cohen_kappa_score, hamming_loss, jaccard_score
include("pairwise/pairwise.jl"); export minkowski_distance, euclidian_distance, manhattan_distance, chebyshev_distance, braycurtis_distance, canberra_distance, cityblock_distance, mahalanobis_distance,
correlation, cosine_distance, cosine_similarity
include("regression/regression.jl"); export max_error, mean_absolute_error, mean_squared_error, mean_squared_log_error, median_absolute_error, mean_absolute_percentage_error

using .classification
using .pairwise
using .regression

end


text:
# This file was generated by JuDoc, do not modify it. # hide
iris2 = coerce(iris, :PetalWidth => OrderedFactor)
first(iris2[[:PetalLength, :PetalWidth]], 1) |> pretty


text:
### A Pluto.jl notebook ###
# v0.16.4

using Markdown
using InteractiveUtils

# ╔═╡ 94075d98-522a-11eb-00bd-8709e96862e9
begin
	import Pkg
	Pkg.activate(Base.current_project())
	using ReinforcementLearning
	using Statistics
	using Flux
	using Plots
end

# ╔═╡ 285f8200-522a-11eb-0260-a3c27fa0762a
md"""
# Chapter 7.2 n-step Sarsa
"""

# ╔═╡ b299f0f4-522a-11eb-040a-e1a9cebe840e
env = RandomWalk1D(N=21)

# ╔═╡ c39deeb6-522a-11eb-1e1e-5193f0ad9fe0
ns, na = length(state_space(env)), length(action_space(env))

# ╔═╡ c7c57a3e-522a-11eb-0fda-91c3aff791d7
true_values = -1:0.1:1

# ╔═╡ dd9e49c6-522a-11eb-2a68-0fc5a404a0cd
md"""
Again, we first define a hook to calculate RMS
"""

# ╔═╡ efe41c64-522a-11eb-20d9-f5b14b5c138c
begin
	struct RecordRMS <: AbstractHook
		rms::Vector{Float64}
		RecordRMS() = new([])
	end
	(f::RecordRMS)(::PostEpisodeStage, agent, env) = push!(f.rms, sqrt(mean((agent.policy.learner.approximator.table[2:end-1] - true_values[2:end-1]).^2)))
end

# ╔═╡ f7443b06-522a-11eb-3801-5525ecada397
function run_once(α, n)
    env = RandomWalk1D(N=21)
    agent = Agent(
        policy=VBasedPolicy(
            learner=TDLearner(
                approximator=TabularVApproximator(;n_state=ns, opt=Descent(α)),
                method=:SRS,
                n=n
                ),
            mapping= (env, V) -> rand(1:na)
            ),
        trajectory=VectorSARTTrajectory()
    )
    hook = RecordRMS()
    run(agent, env, StopAfterEpisode(10; is_show_progress=false), hook)
    mean(hook.rms)
end

# ╔═╡ 0b1b0448-522b-11eb-0264-9131eabe2525
begin
	A = 0.:0.05:1.0
	p = plot(xlabel="α", ylabel="Average RMS error")
	for n in [2^i for i in 0:9]
		avg_rms = Float64[]
		for α in A
			rms = []
			for _ in 1:100
				push!(rms, run_once(α, n))
			end
			push!(avg_rms, mean(rms))
		end
		plot!(p, A, avg_rms, label="n = $n")
	end

	ylims!(p, 0.25, 0.55)
	p
end

# ╔═╡ Cell order:
# ╟─285f8200-522a-11eb-0260-a3c27fa0762a
# ╠═94075d98-522a-11eb-00bd-8709e96862e9
# ╠═b299f0f4-522a-11eb-040a-e1a9cebe840e
# ╠═c39deeb6-522a-11eb-1e1e-5193f0ad9fe0
# ╠═c7c57a3e-522a-11eb-0fda-91c3aff791d7
# ╟─dd9e49c6-522a-11eb-2a68-0fc5a404a0cd
# ╠═efe41c64-522a-11eb-20d9-f5b14b5c138c
# ╠═f7443b06-522a-11eb-3801-5525ecada397
# ╠═0b1b0448-522b-11eb-0264-9131eabe2525


text:
# GeometryTypes doesn't define an Ellipsoid type yet, so we'll make one ourselves!
struct HyperEllipsoid{N, T} <: GeometryPrimitive{N, T}
    center::Point{N, T}
    radii::Vec{N, T}
end

GeometryTypes.origin(geometry::HyperEllipsoid{N, T}) where {N, T} = geometry.center
radii(geometry::HyperEllipsoid{N, T}) where {N, T} = geometry.radii
center(geometry::HyperEllipsoid) = origin(geometry)

@deprecate HyperCylinder(length::T, radius) where {T} Cylinder{3, T}(Point(0., 0., 0.), Point(0, 0, length), radius)

struct PointCloud{T, Point <: StaticVector{3, T}, C <: Colorant} <: AbstractGeometry{3, T}
    position::Vector{Point}
    color::Vector{C}
end

function PointCloud(position::AbstractVector{<:AbstractVector{T}},
           color::AbstractVector{C}=RGB{Float32}[]) where {T, C <: Colorant}
    PointCloud{T, Point{3, T}, C}(position, color)
end

struct Triad <: AbstractGeometry{3, Float64}
    scale::Float64

    Triad(scale=20.0) = new(scale)
end

center(geometry::HyperRectangle) = minimum(geometry) + 0.5 * widths(geometry)
center(geometry::HyperCube) = minimum(geometry) + 0.5 * widths(geometry)
center(geometry::HyperSphere) = origin(geometry)
center(geometry::Cylinder) = origin(geometry) + geometry.extremity / 2

"""
$(SIGNATURES)

Different tools disagree about what various geometric primitives mean. For example,
GeometryTypes.jl considers the "origin" of a cube to be its bottom-left corner, where
DrakeVisualizer and MeshCat consider its origin to be the center. The
intrinsic_transform(g) returns the transform from the GeometryTypes origin to the
MeshCat origin.
"""
intrinsic_transform(g) = IdentityTransformation()
intrinsic_transform(g::HyperRectangle) = Translation(center(g)...)
intrinsic_transform(g::HyperSphere) = Translation(center(g)...)
intrinsic_transform(g::HyperEllipsoid) = Translation(center(g)...) ∘ LinearMap(SDiagonal(radii(g)...))
intrinsic_transform(g::HyperCube) = Translation(center(g)...)

function intrinsic_transform(g::Cylinder{3})
    # Three.js wants a cylinder to lie along the y axis
    R = rotation_between(SVector(0, 1, 0), g.extremity)
    Translation(center(g)) ∘ LinearMap(R)
end


text:
using PGA2D
using Documenter

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

makedocs(;
    modules=[PGA2D],
    authors="Andreas Tell <atell@soundtheory.com> and contributors",
    repo="https://github.com/ATell-SoundTheory/PGA2D.jl/blob/{commit}{path}#{line}",
    sitename="PGA2D.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ATell-SoundTheory.github.io/PGA2D.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/ATell-SoundTheory/PGA2D.jl",
)


text:
sabb(x::String, val::Integer, happ::Integer)::String =
    +(count( c -> c ∈ "sabbatical",x),val,happ) > 22 ? "Sabbatical! Boom!" : "Back to your desk, boy."

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
## -----------  Make deviations for testing
# Names are :d1 etc

function make_matrix_deviation(devNo :: Integer; offset :: Float64 = 0.0)
    sizeV = (5, 4);
    modelV = devNo .+ (1 : sizeV[1]) .+ (1 : sizeV[2])' .+ 0.1 .+ offset;
    idxV = [2:4, 2:3];
    dataV = modelV[idxV...] .+ 0.7;
    wtV = dataV .+ 0.9;
    name, shortStr, longStr, fmtStr = dev_info(devNo);
    d = Deviation{Float64}(name = name,
        modelV = modelV, dataV = dataV, wtV = wtV, idxV = idxV,
        normP = 1.5,
        shortStr = shortStr, longStr = longStr);
    return d;
end

function make_deviation(devNo :: Integer; offset :: Float64 = 0.0)
    dataV = devNo .+ collect(range(2.1, 3.4, length = 5));
    modelV = dataV .+ 0.7 .+ offset;
    wtV = dataV .+ 0.9;
    name, shortStr, longStr, fmtStr = dev_info(devNo);
    d = Deviation{Float64}(name = name,
        modelV = modelV, dataV = dataV, wtV = wtV,
        normP = 1.8,
        shortStr = shortStr, longStr = longStr);
    return d;
end

function make_bounds_deviation(devNo :: Integer, insideBounds :: Bool)
    modelM = collect(1 : 5) .+ collect(2 : 4)' .+ 0.5;
    lbM = modelM .- 0.1;
    ubM = modelM .+ 0.1;
    wtM = modelM .+ 1.3;
    if !insideBounds
        modelM[3] = ubM[3] + 0.01;
    end

    name, shortStr, longStr, fmtStr = dev_info(devNo);
    d = BoundsDeviation{Float64}(name = name,  modelV = modelM,  lbV = lbM,  ubV = ubM,
        wtV = wtM,  shortStr = shortStr, longStr = longStr, fmtStr = fmtStr);
    return d;
end

function make_penalty_deviation(devNo :: Integer)
    modelM = collect(1 : 5) .+ collect(2 : 4)' .+ 0.5;
    name, shortStr, longStr, fmtStr = dev_info(devNo);
    d = PenaltyDeviation{Float64}(name = name,  modelV = modelM,
        scalarDevFct = penalty_dev_fct,
        shortStr = shortStr, longStr = longStr);
    return d;

end

function penalty_dev_fct(modelM)
    return sum(modelM) .- 3.0
end

function make_scalar_deviation(devNo :: Integer)
    name, shortStr, longStr, fmtStr = dev_info(devNo);
    modelV = devNo * 1.1;
    dataV = devNo * 2.2;
    return ScalarDeviation{Float64}(name = name, modelV = modelV,
        dataV = dataV, shortStr = shortStr,
        longStr = longStr)
end

# Make regression deviation for testing.
# No of regressors = devNo + 1
function make_regression_deviation(devNo :: Integer)
    name, shortStr, longStr, fmtStr = dev_info(devNo);
    nc = devNo + 2;
    coeffNameV = Symbol.("beta" .* string.(1 : nc));
    mCoeffV = collect(range(0.1, 0.9, length = nc));
    mSeV = collect(range(0.3, 0.1, length = nc));
    rModel = RegressionTable(coeffNameV, mCoeffV, mSeV);
    rData = RegressionTable(coeffNameV, mCoeffV .+ 0.1, mSeV .+ 0.2);
    return RegressionDeviation{Float64}(name = name,
        shortStr = shortStr, longStr = longStr,
        modelV = rModel, dataV = rData, normP = 1.2)
end

function dev_info(devNo :: Integer)
    name = Symbol("d$devNo");
    shortStr = "dev$devNo";
    longStr = "Deviation $devNo"
    fmtStr = "%.2f";
    return name, shortStr, longStr, fmtStr
end

function make_deviation_vector(n :: Integer; offset :: Float64 = 0.0)
    dv = DevVector();
    for j = 1 : n
        d = make_matrix_deviation(j; offset = offset + 0.2 * j);
        append!(dv, d);
    end
    return dv
end

# ---------------


text:
# Julia 1.3.0

function dist(a, b)
	return sqrt( (a[1] - b[1]) * (a[1] - b[1]) + (a[2] - b[2]) * (a[2] - b[2]) )
end

function sub(a,b)
	return (a[1] - b[1], a[2] - b[2])
end

function angle2(a, b)
	return atan(b[2] - a[2], b[1] - a[1])
end

function part1(coord, pairs)
	rv = []
	for t in pairs
		u1 = sub(t[2], coord)
		x = angle2((0,0), u1)
		md = rad2deg(x)
		if !(md in rv)
			append!(rv, md)
		end
	end
	return rv
end

function part2(coord, pairs)
	rv = []
	for t in pairs
		u1 = sub(t[2], coord)
		x = angle2(coord, t[2])
		md = rad2deg(x)
		push!(rv, (t[2][1], t[2][2], md))
	end
	return rv
end

function los1(m)
	rv = []
	for y in 1:length(m)
		for x in 1:length(m[y])
			if m[y][x] == '.'
				continue
			end
			ck = []
			for y1 in 1:length(m)
				for x1 in 1:length(m[y])
					if m[y1][x1] == '.'
						continue
					end
					if x == x1 && y == y1
						continue
					end
					push!(ck, ((x,y),(x1,y1)))
				end
			end
			angles = part1((x,y), ck)
			push!(rv, (x,y,length(angles)))
		end
	end
	mx = 0
	bs = (0,0)
	for r in rv
		if r[3] > mx
			bs = (r[1], r[2])
			mx = r[3]
		end
	end
	println()
	println("Part 1: ($(bs[1]-1),$(bs[2]-1)) :: $(mx)")
	println()
	return (bs, mx)
end

function m(x, bs, lang)
	rv = ((((x[3] - lang) % 360) + 360) % 360) * 1000 + dist(bs, x)
	#println("    m(",x," ",bs," ",lang," ",dist(bs,x),")=",rv)
	return rv
end

function los2(mm, bs)
	ck = []
	for y1 in 1:length(mm)
		for x1 in 1:length(mm[y1])
			if mm[y1][x1] == '.'
				continue
			end
			if y1 == bs[2] && x1 == bs[1]
				continue
			end
			push!(ck, ((0,0),(x1,y1)))
		end
	end
	aa = part2(bs, ck)
	for i in 1:length(aa)
		aa[i] = (aa[i][1], aa[i][2], (aa[i][3] + 90 + 360) % 360)
	end
	for a in sort(aa)
		#println(a," ",a[3] * 10000 + dist(bs, a))
	end
	shot = []
	lang = 0
	println("left: ",length(aa)," shot: ",length(shot))
	println()

	tgt = (0,0,0)
	lang = 0
	while true
		if length(aa) < 1
			break
		end
		if length(shot) >= 200
			break
		end
		tmpa = sort(aa, by = x -> m(x, bs, lang))
		tgt = tmpa[1]
		idx = 1
		hx = m(tgt, bs, lang)
		println("Target: (",tgt[1],", ",tgt[2],", ",tgt[3],") m: ",hx ," l: ", lang)
		aa = deleteat!(tmpa, idx)
		push!(shot, tgt)
		lang = (tgt[3] + 0.001) % 360
		println("left: ", length(aa), " shot: ", length(shot), " last: ", tgt, " lang: ", lang)
		println()
	end
	println("Part 2: (", tgt[1]-1,"/",tgt[2]-1,") :: ", (tgt[1]-1)*100 + (tgt[2]-1))
	return tgt
end

if length(ARGS) <1
	println("Usage: ", PROGRAM_FILE, " /path/to/file")
	exit()
end
filename = ARGS[1]

println("File: $(filename)")
am = readlines(filename)
(bs, mx) = los1(am)
tgt = los2(am, bs)
println("---")
println("Part 1: ($(bs[1]-1),$(bs[2]-1)) :: $(mx)")
println("Part 2: (", tgt[1]-1,"/",tgt[2]-1,") :: ", (tgt[1]-1)*100 + (tgt[2]-1))

max_stars_repo_name: winks/adventofcode

max_issues_repo_name: winks/adventofcode

max_forks_repo_name: winks/adventofcode


text:
function value_of(::Type{LagrangianPoints}, arg0::JString)
    return jcall(LagrangianPoints, "valueOf", LagrangianPoints, (JString,), arg0)
end

function values(::Type{LagrangianPoints})
    return jcall(LagrangianPoints, "values", Vector{LagrangianPoints}, ())
end


text:
export
    LineNode,
    set!

import ..Nodes.draw

mutable struct LineNode <: Ranger.AbstractNode
    base::Nodes.NodeData
    transform::Nodes.TransformProperties{Float64}

    mesh::Geometry.Mesh

    color::Rendering.Palette

    function LineNode(world::Ranger.World, name::String, parent::Ranger.AbstractNode)
        o = new()

        o.base = Nodes.NodeData(name, parent, world)
        o.transform = Nodes.TransformProperties{Float64}()
        o.mesh = Geometry.Mesh()
        o.color = Rendering.White()

        Geometry.add_vertex!(o.mesh, 0.0, 0.0)
        Geometry.add_vertex!(o.mesh, 0.0, 0.0)

        Geometry.build!(o.mesh)

        o
    end
end

function Nodes.draw(node::LineNode, context::Rendering.RenderContext)
    if Nodes.is_dirty(node)
        Rendering.transform!(context, node.mesh)
        Nodes.set_dirty!(node, false)
    end

    Rendering.set_draw_color(context, node.color)

    v1 = node.mesh.bucket[1]
    v2 = node.mesh.bucket[2]
    Rendering.draw_line(context, v1.x, v1.y, v2.x, v2.y);

end

function set!(node::LineNode, x1::Float64, y1::Float64, x2::Float64, y2::Float64)
    Geometry.set!(node.mesh.vertices[1], x1, y1)
    Geometry.set!(node.mesh.vertices[2], x2, y2)

    Nodes.set_dirty!(node, true);
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
module Bautin

using DiffEqBase: ODEProblem
using Parameters: @with_kw, @unpack
using StaticArrays: SVector
using Setfield: @lens

using ...Bifurcations: BifurcationProblem
using ...Codim2: DiffEqCodim2Problem

@with_kw struct BautinParam{B1, B2, S}
    β₁::B1 = -1.0
    β₂::B2 = -1.0
    s::S = -1
end

function f(u::SVector, p, t)
    @unpack β₁, β₂, s = p
    x, y = u
    #=
    ρ² = x^2 + y^2
    ρ⁴ = ρ²^2
    =#
    ρ² = x * x + y * y
    ρ⁴ = ρ² * ρ²
    a = β₂ * ρ² + s * ρ⁴
    dx = β₁ * x -      y + a * x
    dy =      x + β₁ * y + a * y
    return SVector(dx, dy)
end

# TODO: make complex Bautin normal form work
#=
function f(u::SVector, p, t)
    dz = f(u[1] + u[2] * im, p, t)
    return SVector(real(dz), imag(dz))
end

function f(z::Complex, p, t)
    @unpack β₁, β₂, s = p
    ρ² = abs2(z)
    ρ⁴ = ρ²^2
    return (β₁ + 1im) * z + β₂ * z * ρ² + s * z * ρ⁴
end
=#

f(u::A, p, t) where {A <: AbstractVector} = A(f(SVector{2}(u), p, t))

function f(du, u, p, t)
    du .= f(SVector{2}(u), p, t)
    nothing
end

make_prob(
        p = BautinParam();
        u0 = SVector(0.0, 0.0),
        tspan = (0.0, 30.0),
        ode = ODEProblem{!(u0 isa SVector)}(f, u0, tspan, p),
        param_axis = (@lens _.β₁),
        t_domain = (-2.0, 2.0),
        kwargs...) =
    BifurcationProblem(ode, param_axis, t_domain;
                       kwargs...)

prob = make_prob()
ode = prob.p.de_prob
u0 = ode.u0
param_axis = prob.p.param_axis

function make_codim2_prob(
        p = ode.p;
        u0 = SVector(0.0, 0.0),
        tspan = (0.0, 30.0),
        ode = ODEProblem{!(u0 isa SVector)}(f, u0, tspan, p),
        param_axis1 = (@lens _.β₁),
        param_axis2 = (@lens _.β₂),
        t_domain = ([-2.0, -2.0], [2.0, 2.0]),
        kwargs...)
    return DiffEqCodim2Problem(
        ode,
        param_axis1,
        param_axis2,
        t_domain;
        kwargs...)
end

end  # module


text:
# This file is a part of ROOTFramework.jl, licensed under the MIT License (MIT).

using Cxx

cxxinclude("TCanvas.h")