text:
# helper
colnorms(A::AbstractMatrix) = [norm(a) for a in eachcol(A)]

function samesupport(x::AbstractVector, y::AbstractVector)
    samesupport(sparse(x), sparse(y))
end
function samesupport(x::SparseVector, y::SparseVector)
    sort!(x.nzind) == sort!(y.nzind)
end

########################## synthetic data generators ###########################
# creates random k-sparse vector with ±1 as entries, or gaussian depending on flag
function sparse_vector(m::Int, k::Int, gaussian::Bool = false)
    m ≥ k || throw("m = $m < $k = k")
    x = spzeros(m)
    ind = sort!(sample(1:m, k, replace = false))
    x[ind] .= gaussian ? randn(k) : rand((-1, 1), k)
    return x
end

function sparse_data(;n = 32, m = 64, k = 3, rescaled = true)
    A = randn(n, m)
    if rescaled
        ε = 1e-6
        A .-= ε*mean(A, dims = 1) # this apparently is good practice, better conditioning
        A ./= sqrt.(sum(abs2, A, dims = 1))
    end
	x = sparse_vector(m, k)
    b = A*x
    return A, x, b
end
const gaussian_data = sparse_data
sparse_data(n, m, k) = sparse_data(n = n, m = m, k = k)
function correlated_data(n, m, k; normalized = true)
    U = randn(n, n)
    V = randn(n, m)
    S = Diagonal([1/i^2 for i in 1:n])
    A = U*S*V
    # normalize
    if normalized
        A ./= sqrt.(sum(abs2, A, dims = 1))
    end
    x = sparse_vector(m, k)
    b = A*x
    A, x, b
end
const coherent_data = correlated_data

# random perturbation of b with norm δ
function perturb!(b::AbstractVector, δ::Real)
    e = randn(size(b))
    e *= δ/norm(e)
    b .+= e # perturb
end
perturb(b, δ) = perturb!(copy(b), δ)

########################### dictionary preconditioners #########################
function normalize!(A::AbstractVecOrMat)
    A ./= colnorms(A)'
end

# see "On the Uniqueness of Nonnegative Sparse Solutions to Underdetermined Systems of Equations", Bruckstein 2008
function preconditioner(ε::Real)
    function p!(y, x)
		μ = mean(x, dims = 1)
        @. y = x - (1-ε) * μ
    end
    p!(x) = p!(x, x)
    return p!
end
# precondition!(A::AbstractVecOrMat, ε::Real) = preconditioner(ε)(A)

# related but not identical to the one used in:
# "Preconditioned Multiple Orthogonal Least Squares and Applications in Ghost Imaging via Sparsity Constraint"
# preconditions using U*S⁻¹*U' where U, S, V is svd of A
function preconditioner(A::AbstractMatrix{Float64}, min_σ::Real = 1e-6)
    svdA = svd(A)
    U, S, V = svdA.U, svdA.S, svdA.V
	k = size(U, 2)
	function p!(y, x, z)
		size(y) == size(x) || throw(DimensionMismatch("size(x) ≠ size(y)"))
		mul!(z, U', x)
		z ./= max.(S, min_σ)
		mul!(y, U, z)
	end
	p!(y::AbstractMatrix, x::AbstractMatrix) = p!(y, x, similar(x, k, size(x, 2)))
	p!(y::AbstractVector, x::AbstractVector) = p!(y, x, similar(x, k))
	p!(x::AbstractVecOrMat) = p!(x, x)
	return p!
end
precondition!(A::AbstractMatrix) = preconditioner(A)(A)

############################# dictionary analysis ##############################
# calculates mutual coherence
coherence(A::AbstractMatrix) = babel(A, 1)

# Babel function, see GREED IS GOOD: ALGORITHMIC RESULTS FOR SPARSE APPROXIMATION
babel(A::AbstractMatrix, k::Integer) = cumbabel(A, k)[k]

Base.cumsum!(x::AbstractArray) = cumsum!(x, x)
# calculates all babel function values from 1 to k
function cumbabel(A::AbstractMatrix, k::Integer)
    μ₁ = zeros(eltype(A), k)
    inner = similar(A, size(A, 2))
    for (i, ai) in enumerate(eachcol(A))
        mul!(inner, A', ai)
        @. inner = abs(inner)
        inner[i] = 0 # inner product with self does not count
        partialsort!(inner, 1:k, rev = true) # pick largest k inner products
        innerk = @view inner[1:k]
        μ₁ .= max.(μ₁, cumsum!(innerk))
    end
    return μ₁
end

########################### SparseArrays conveniences ##########################
function addindex!(x::SparseVector, AiQR::UpdatableQR, a::AbstractVector, i::Int)
	if i ∉ x.nzind # if the index is not already in the support set
	    x[i] = NaN   # add non-zero index to active set
	    # efficiently update qr factorization using Givens rotations
	    qr_i = findfirst(==(i), x.nzind) # index in qr where new atom should be added
	    add_column!(AiQR, a, qr_i)
	end
    return x, AiQR
end
const addind! = addindex!

function addindex!(x::SparseVector, AiQR::UpdatableQR, A::AbstractMatrix, indices::AbstractVector)
    for (i, j) in enumerate(indices)
		addindex!(x, AiQR, @view(A[:, i]), j)
	end
	return x, AiQR
end

# drops the ith index-value pair of x, if it is in nzind
function dropindex!(x::SparseVector, i::Int)
	j = findfirst(==(i), x.nzind)
	if !isnothing(j)
		deleteat!(x.nzind, j)
		deleteat!(x.nzval, j)
	end
    return x
end
const dropind! = dropindex!
dropind!(x::SparseVector, i::AbstractVector{<:Int}) = dropind!.((x,), i)
# simultaneously deletes an index of x and removes a column of a UpdatableQR
function dropindex!(x::SparseVector, AiQR::UpdatableQR, i::Int)
    j = findfirst(==(i), x.nzind)
    if !isnothing(j)
		_dropindex!(x, AiQR, j)
    end
    return x, AiQR
end
# WARNING: drops index j into x,nzval, NOT into x, as dropindex!
function _dropindex!(x::SparseVector, AiQR::UpdatableQR, j::Int)
	deleteat!(x.nzind, j)
	deleteat!(x.nzval, j)
	remove_column!(AiQR, j)
	return x, AiQR
end

function SparseArrays.droptol!(x::SparseVector, AiQR::UpdatableQR, tol::Real)
    for i in reverse(eachindex(x.nzind)) # reverse is necessary to not mess up indexing into QR factorization
        if abs(x.nzval[i]) ≤ tol
            remove_column!(P.AiQR, i)
        end
    end
    return droptol!(x, tol), AiQR
end

# TODO: this should be in util somewhere, added in Julia 1.7?
function Base.findmin(f, x::AbstractVector)
    k, m = 0, Inf
    for (i, xi) in enumerate(x)
        fxi = f(xi)
        if fxi < m
            k, m = i, fxi
        end
    end
    return m, k
end
Base.argmin(f, x::AbstractVector) = findmin(f, x)[2]
function Base.findmax(f, x::AbstractVector)
	g(x) = -f(x)
	m, k = findmin(g, x)
	return -m, k
end
Base.argmax(f, x::AbstractVector) = findmax(f, x)[2]

# using StatsBase: mean, std
# # makes x mean zero along dimensions dims
# function center!(x::AbstractVecOrMat, ε = 1e-6; dims = :)
#     μ, σ = mean(x, dims = dims), std(x, dims = dims)
#     @. x = (x - μ) / σ
#     # @. x = (x - (1-ε)*μ) / σ
#     return x, μ, σ
# end
# Base.inv(::typeof(center!)) = uncenter
# function uncenter(x::AbstractVecOrMat, μ, σ)
#     @. x = x * σ + μ
#     x, y
# end


text:
using DataFrames
using CSV
using LearningStrategies: Tracer
using SoftComputingFinal

using Plots
pyplot()

const ValType = Type{Val{T}} where T
const datasets = Dict{String, Tuple{Function, ValType, ValType}}(
    "testdata" => (load_testdata, Val{2}, Val{2}),
    "glass" => (load_glass, Val{9}, Val{7}),
    "ionosphere" => (load_ionosphere, Val{34}, Val{2}),
    "segmentation" => (load_segmentation, Val{19}, Val{7}))

getsampler(::Type{Val{M}}, ::Type{Val{N}}) where {M, N} =
    RampedSplitSampler{M, N}(6, 0.5, 0.5,
                             crange = (-10.0, 10.0),
                             maxvars = min(3, M))

kbest(model, k) = view(model.population, selectperm(model.population_fitnesses, 1:k))
best(model) = model.population[indmax(model.population_fitnesses)]

function evaluatedataset(name, N; folds = 10, repetitions = 3, pareto_sample = 50)
    popsize = 250

    load_data, nvars_t, nclasses_t = datasets[name]
    data = load_data()
    accuracies = DataFrame(generation = Int[],
                           accuracy = Float64[])
    sizes = DataFrame(size = Int[], fold = Int[], run = Int[])

    D = size(data)[1]
    shuffled = randperm(D)
    train_indices, val_indices = kfolds(D, folds)

    rngs = randjump(Base.GLOBAL_RNG, Threads.nthreads()) # we need one rng per thread!
    update_lock = Threads.SpinLock()
    function locked(block)
        lock(update_lock)
        block()
        unlock(update_lock)
    end

    Threads.@threads for f = 1:folds
        training_data = view(data, shuffled[train_indices[f]])
        validation_data = view(data, shuffled[val_indices[f]])

        fitness, _ = create_fitness(training_data, nvars_t, nclasses_t,
                                    depth_factor = 2.0,
                                    size_factor = 0.0)
        _, validation_accuracy = create_fitness(validation_data, nvars_t, nclasses_t,
                                                depth_factor = 2.0,
                                                size_factor = 0.0)

        for r = 1:repetitions
            tracer = Tracer(Float64, (m, i) -> validation_accuracy(best(m)))
            initializer = getsampler(nvars_t, nclasses_t)

            pop, trace = runssgp(fitness, popsize, initializer, N,
                                 max_depth = 30,
                                 tracer = tracer,
                                 verbose = false,
                                 rng = rngs[Threads.threadid()])

            locked() do
                append!(accuracies, DataFrame(generation = 1:N, accuracy = collect(trace)))
            end
            locked() do
                append!(sizes, DataFrame(size = treesize.(pop), fold = f, run = r))
            end
        end
    end

    mean_accuracies = by(accuracies, :generation) do df
        x̄ = mean(df[:accuracy])
        σ̂ = std(df[:accuracy], mean = x̄)
        n = length(df[:accuracy])
        DataFrame(mean = x̄, ci = 1.96σ̂ / √n)
    end

    return mean_accuracies, sizes
end

function main()
    if length(ARGS) == 2
        dataset, N = ARGS[1], parse(Int, ARGS[2])
    elseif length(ARGS) == 1
        dataset, N = ARGS[1], 1000
    else
        error("Usage: $PROGRAM_FILE <dataset> [<generations>]")
    end

    mean_accuracies, sizes = evaluatedataset(dataset, N)
    # display(mean_accuracies)

    plt_acc = plot(mean_accuracies[:generation], mean_accuracies[:mean],
                   ribbon = mean_accuracies[:ci],
                   fillalpha = 0.2,
                   xlabel = "Generation",
                   ylabel = "Validation accuracy", ylims = (0, 1),
                   title = "Average validation accuracy\n and 95% confidence interval\n over time",
                   legend = :none)
    png(plt_acc, "$dataset-accuracies.png")

    plt_sizes = histogram(sizes[:size], normalize = true,
                          xlabel = "Size", ylabel = "Frequency",
                          title = "Size distribution in last generation",
                          legend = :none)
    png(plt_sizes, "$dataset-sizes.png")

    CSV.write("$dataset-accuracies.csv", mean_accuracies)
    CSV.write("$dataset-sizes.csv", sizes)
end

main()

ext: jl

max_stars_repo_name: phipsgabler/soft-computing-final

max_issues_repo_name: phipsgabler/soft-computing-final

max_forks_repo_name: phipsgabler/soft-computing-final


text:
"""
θ-method ODE solver
"""
struct ThetaMethod <: ODESolver
  nls::NonlinearSolver
  dt::Float64
  θ::Float64
  function ThetaMethod(nls,dt,θ)
    if θ > 0.0
      return new(nls,dt,θ)
    else
      return ForwardEuler(nls,dt)
    end
  end
end

BackwardEuler(nls,dt) = ThetaMethod(nls,dt,1.0)
MidPoint(nls,dt) = ThetaMethod(nls,dt,0.5)

function solve_step!(uf::AbstractVector,
                     solver::ThetaMethod,
                     op::ODEOperator,
                     u0::AbstractVector,
                     t0::Real,
                     cache) # -> (uF,tF)

  dt = solver.dt
  solver.θ == 0.0 ? dtθ = dt : dtθ = dt*solver.θ
  tθ = t0+dtθ

  if cache === nothing
    ode_cache = allocate_cache(op)
    vθ = similar(u0)
    nl_cache = nothing
  else
    ode_cache, vθ, nl_cache = cache
  end

  ode_cache = update_cache!(ode_cache,op,tθ)

  nlop = ThetaMethodNonlinearOperator(op,tθ,dtθ,u0,ode_cache,vθ)

  nl_cache = solve!(uf,solver.nls,nlop,nl_cache)

  if 0.0 < solver.θ < 1.0
    uf = uf*(1.0/solver.θ)-u0*((1-solver.θ)/solver.θ)
  end

  cache = (ode_cache, vθ, nl_cache)

  tf = t0+dt
  return (uf,tf,cache)

end

"""
Nonlinear operator that represents the θ-method nonlinear operator at a
given time step, i.e., A(t,u_n+θ,(u_n+θ-u_n)/dt)
"""
struct ThetaMethodNonlinearOperator <: NonlinearOperator
  odeop::ODEOperator
  tθ::Float64
  dtθ::Float64
  u0::AbstractVector
  ode_cache
  vθ::AbstractVector
end

function residual!(b::AbstractVector,op::ThetaMethodNonlinearOperator,x::AbstractVector)
  uθ = x
  vθ = op.vθ
  vθ = (x-op.u0)/op.dtθ
  residual!(b,op.odeop,op.tθ,(uθ,vθ),op.ode_cache)
end

function jacobian!(A::AbstractMatrix,op::ThetaMethodNonlinearOperator,x::AbstractVector)
  uF = x
  vθ = op.vθ
  vθ = (x-op.u0)/op.dtθ
  z = zero(eltype(A))
  fillstored!(A,z)
  jacobians!(A,op.odeop,op.tθ,(uF,vθ),(1.0,1/op.dtθ),op.ode_cache)
end

function allocate_residual(op::ThetaMethodNonlinearOperator,x::AbstractVector)
  allocate_residual(op.odeop,x,op.ode_cache)
end

function allocate_jacobian(op::ThetaMethodNonlinearOperator,x::AbstractVector)
  allocate_jacobian(op.odeop,x,op.ode_cache)
end

function zero_initial_guess(op::ThetaMethodNonlinearOperator)
  x0 = similar(op.u0)
  fill!(x0,zero(eltype(x0)))
  x0
end


text:
using Genie
using Logging, LoggingExtras
using Dates

const date_format = "yyyy-mm-dd HH:MM:SS"

const logger =  if Genie.config.log_to_file
                  isdir(Genie.LOG_PATH) || mkpath(Genie.LOG_PATH)
                  DemuxLogger(
                    FileLogger(joinpath(Genie.LOG_PATH, "$(Genie.config.app_env)-$(Dates.today()).log"), always_flush = true, append = true),
                    ConsoleLogger(stdout, Genie.config.log_level),
                    include_current_global = false
                  )
                else
                  ConsoleLogger(stdout, Genie.config.log_level)
                end

timestamp_logger(logger) = TransformerLogger(logger) do log
  merge(log, (; message = "$(Dates.format(now(), date_format)) $(log.message)"))
end

DemuxLogger(MinLevelLogger(logger, Genie.config.log_level), include_current_global = false) |> timestamp_logger |> global_logger


text:
using DataStructures
using DynamicalBilliards: ismagnetic, find_cyclotron

mutable struct ParticleObservable{P<:AbstractParticle{Float32}}
    # Fields necessary for simulation
    p::P   # particle
    i::Int # index of obstacle to be collided with
    tmin::Float32 # time to next collision
    t::Float32    # current time (resets at each collision)
    n::Int        # number of collisions done so far
    T::Float32    # total time
    # Fields used in plotting
    tail::Observable{CircularBuffer{Point2f0}}
    ξsin::Observable{Point2f0}
end
function ParticleObservable(p, bd, n, ξsin = Point2f0(0, 0)) # initializer
    i, tmin::Float32, cp = next_collision(p, bd)
    cb = CircularBuffer{Point2f0}(n)
    for i in 1:n; push!(cb, Point2f0(p.pos)); end
    ParticleObservable(p, i, tmin, 0f0, 0, 0f0, Observable.((cb, ξsin))...)
end
const ParObs = ParticleObservable

function rebind_partobs!(p::ParticleObservable, p0::AbstractParticle, bd, ξsin = p.ξsin[])
    i, tmin::Float32, cp = next_collision(p0, bd)
    ξ = sφ = 0f0 # TODO: Use boundary map on cp
    p.p.pos = p0.pos
    p.p.vel = p0.vel
    ismagnetic(p.p) && (p.p.center = find_cyclotron(p.p))
    p.i, p.tmin, p.t, p.n, p.T = i, tmin, 0f0, 0, 0f0
    L = length(p.tail[])
    append!(p.tail[], [Point2f0(p0.pos) for i in 1:L])
    p.tail[] = p.tail[]
    if ξsin ≠ nothing
        p.ξsin[] = ξsin # This can only be updated from bmap, which gives selection directly
    end
end

function animstep!(parobs, bd, dt, updateplot = true, intervals = nothing)
    if parobs.t + dt - parobs.tmin > 0
        rt = parobs.tmin - parobs.t # remaining time
        animbounce!(parobs, bd, rt, updateplot, intervals)
    else
        propagate!(parobs.p, dt)
        parobs.t += dt
        push!(parobs.tail[], parobs.p.pos)
        if updateplot
            parobs.tail[] = parobs.tail[] # trigger update
        end
    end
    return
end

function animbounce!(parobs, bd, rt, updateplot = true, intervals = nothing)
    propagate!(parobs.p, rt)
    DynamicalBilliards._correct_pos!(parobs.p, bd[parobs.i])
    DynamicalBilliards.resolvecollision!(parobs.p, bd[parobs.i])
    ismagnetic(parobs.p) && (parobs.p.center = find_cyclotron(parobs.p))
    if intervals ≠ nothing
        ξ, sφ = to_bcoords(parobs.p.pos, parobs.p.vel, bd[parobs.i])
        ξ += intervals[parobs.i]
        parobs.ξsin[] = (ξ, sφ)
    end
    i, tmin::Float32, = next_collision(parobs.p, bd)
    parobs.i = i
    parobs.tmin = tmin
    parobs.t = 0f0
    parobs.T += tmin
    parobs.n += 1
    push!(parobs.tail[], parobs.p.pos)
    if updateplot
        parobs.tail[] = parobs.tail[] # trigger update
    end
    return
end


text:
#
# Contains inforation about jevo project
#
# @author DeadbraiN
#
module About
  function version()
    "jevo v1.0.1"
  end
  function description()
    "jevo - evolutionary progorgamming and genetic algorithm related application,
which manages digital organisms in computer memory. They leave in 2D virtual world,
evolve, compete and reproduce itself according to evolutionary rules. Main goal of
jevo - simple Artificial Intelligence creation."
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
const DATA_DUMMY = GLfloat[0,0,0]

const DATA_CUBE = GLfloat[
# Bottom
-1.0, -1.0, -1.0,
1.0, -1.0, -1.0,
-1.0, -1.0, 1.0,
1.0, -1.0, -1.0,
1.0, -1.0, 1.0,
-1.0, -1.0, 1.0,

# Top
-1.0, 1.0, -1.0,
-1.0, 1.0, 1.0,
1.0, 1.0, -1.0,
1.0, 1.0, -1.0,
-1.0, 1.0, 1.0,
1.0, 1.0, 1.0,

# Front
-1.0, -1.0, 1.0,
1.0, -1.0, 1.0,
-1.0, 1.0, 1.0,
1.0, -1.0, 1.0,
1.0, 1.0, 1.0,
-1.0, 1.0, 1.0,

# Back
-1.0, -1.0, -1.0,
-1.0, 1.0, -1.0,
1.0, -1.0, -1.0,
1.0, -1.0, -1.0,
-1.0, 1.0, -1.0,
1.0, 1.0, -1.0,

# Left
-1.0, -1.0, 1.0,
-1.0, 1.0, -1.0,
-1.0, -1.0, -1.0,
-1.0, -1.0, 1.0,
-1.0, 1.0, 1.0,
-1.0, 1.0, -1.0,

# Right
1.0, -1.0, 1.0,
1.0, -1.0, -1.0,
1.0, 1.0, -1.0,
1.0, -1.0, 1.0,
1.0, 1.0, -1.0,
1.0, 1.0, 1.0,
]

const DATA_CUBE_VERTEX = [
1f0,  1,  1,  # 0
-1,  1,  1,   # 1
-1, -1,  1,   # 2
1, -1,  1,    # 3
1, -1, -1,    # 4
-1, -1, -1,   # 5
-1,  1, -1,   # 6
1,  1, -1,    # 7
]

const DATA_CUBE_INDEX = GLuint[
0, 1, 2, 2, 3, 0, # Front face
7, 4, 5, 5, 6, 7, # Back face
6, 5, 2, 2, 1, 6, # Left face
7, 0, 3, 3, 4, 7, # Right face
7, 6, 1, 1, 0, 7, # Top face
3, 2, 5, 5, 4, 3  # Bottom face
]

const DATA_PLANE_VERTEX  = [
1f0,  0,  1,  # 0
1, 0, -1,     # 1
-1,  0,  1,   # 2
-1, 0, -1     # 3
]

const DATA_PLANE_INDEX = GLuint[
0, 1, 2, 3 #GL_TRIANGLE_STRIP
#0, 1, 2, 2, 1, 3 #GL_TRIANGLES
]

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
#==============================================================================#
# http.jl
#
# HTTP Requests with retry/back-off and HTTPException.
#
# Copyright Sam O'Connor 2014 - All rights reserved
#==============================================================================#

import URIParser: URI, query_params
import Requests: format_query_str, process_response, open_stream,
                 mimetype, text, bytes, OnBody,
                 BodyDone, EarlyEOF
import HttpCommon: Request, Response, STATUS_CODES
import Base: show, UVError

export HTTPException

type HTTPException <: Exception
    request
    response
end

http_status(e::HTTPException) = e.response.status
headers(e::HTTPException) = e.response.headers
http_message(e::HTTPException) = String(e.response.data)
content_type(e::HTTPException) = get(e.response.headers, "Content-Type", "")

function show(io::IO,e::HTTPException)

    println(io, string("HTTP ", http_status(e), " -- ",
                       e.request.method, " ", e.request.uri, " -- ",
                        http_message(e)))
end

function show_data(data)
    if length(data) > 1000
        println(String(data[1:1000]))
        println("...")
    else
        println(String(data))
    end
end

http_ok(r) = r.status in [200, 201, 202, 204, 206]

function http_attempt(request::Request, return_stream=false)

    if debug_level > 1
        println("$(request.method) $(request.uri)")
        @show request.headers
        show_data(request.data)
    end

    # Start HTTP transaction...
    stream = open_stream(request)

    try

        @sync begin

            @async process_response(stream)

            # Send request data...
            if length(request.data) > 0

                @protected try
                    write(stream, request.data)
                catch e
                    @ignore if isa(e, UVError) && stream.state >= OnBody end
                end
            end
        end

        # Read result...
        response = stream.response
        if !return_stream || !http_ok(response)
            response.data = read(stream)

            if debug_level > 1
                println(response.status)
                @show response.headers
                show_data(response.data)
            end
        end

        # Return on success...
        if stream.state == BodyDone
            if http_ok(response)
                return return_stream ? stream : response
            end
        else
            @assert stream.state == EarlyEOF
            if http_ok(response)
                throw(EOFError())
            end
        end

        # Throw error on failure...
        throw(HTTPException(request, response))

    finally
        if !return_stream
            close(stream)
        end
    end

    assert(false) # Unreachable.
end

function http_request(request::Request, return_stream=false)

    request.headers["Content-Length"] = string(length(request.data))

    @repeat 4 try

        return http_attempt(request, return_stream)

    catch e
        @delay_retry if typeof(e) == UVError end
        @delay_retry if http_status(e) < 200 ||
                        http_status(e) >= 500 end
    end

    assert(false) # Unreachable.
end

function http_request(host::String, resource::String)

    http_request(Request("GET", resource,
                         Dict{String,String}(), UInt8[],
                         URI(host,resource)))
end

#==============================================================================#
# End of file.
#==============================================================================#


text:
# This file was generated, do not modify it. # hide
m21 = machine(model2, corestrict(X, f, 1), corestrict(y, f, 1))
m22 = machine(model2, corestrict(X, f, 2), corestrict(y, f, 2))
m23 = machine(model2, corestrict(X, f, 3), corestrict(y, f, 3))
y21 = predict(m21, restrict(X, f, 1));
y22 = predict(m22, restrict(X, f, 2));
y23 = predict(m23, restrict(X, f, 3));


text:
using MultivariateSeries

X = @ring x1 x2
L = monoms(X,5)

Xi0 = rand(2,3)
w0 = rand(3)

sigma = series(w0, Xi0, L)

w, Xi = decompose(sigma, cst_rkf(3))


text:
#=
findpeaks.jl
Julia translation from Matlab code "findPeaks4.m" from
https://sethares.engr.wisc.edu/vocoders/matlabphasevocoder.html

The comments in that file say that the Matlab version was modified
from findPeaks.m by P. Moller-Nielson 28.3.03, pm-n.
See http://www.daimi.au.dk/~pmn/sound
=#

"""
    findpeaks(amp, max_peak::Int = 50, eps_peak::Real = 0.005)
Given a real-valued vector of amplitudes `amp`,
find up to `max_peak` peaks.
Ignores any amplitudes `≤ eps_peak * maxium(amp)`.

Returns a matrix with 3 columns where each row describes a peak.
"""
function findpeaks(
    amp::AbstractVector{<:Real},
    max_peak::Int = 50, # max number of peaks
    eps_peak::Real = 0.005, # minimum relative height of peaks
)

    N = length(amp)
    n = 3:(N-1)
    peakamp = amp[n] .* (amp[n] .> amp[n .- 1]) .* (amp[n] .> amp[n .+ 1])
    peakpos = zeros(Int, max_peak)

    maxamp = maximum(peakamp)
    npeak = 0
    for p in 1:max_peak
        (m, b) = findmax(peakamp)
        if m ≤ (eps_peak * maxamp)
            break
        end
        peakpos[p] = b + 2
        peakamp[b] = 0
        npeak = p
    end

    peakpos = sort(peakpos)
    peaks = zeros(Int, npeak, 3)

    last_b = 1
    for p in 1:npeak
        b = peakpos[max_peak - npeak + p]
        first_b = last_b + 1
        if p == npeak
            last_b = N
        else
            next_b = peakpos[max_peak - npeak + p + 1]
            rel_min = argmin(amp[b:next_b])
            last_b = b + rel_min - 1
        end
        peaks[p,:] = [first_b, b, last_b]
    end

    return peaks
end


text:
# Spectrum calculator

struct SpectraDiagnosticsParams <: DiagnosticsGroupParams
    nor::Float64
end

"""
    setup_dump_spectra_diagnostics(
        ::ClimateMachineConfigType,
        interval::String,
        out_prefix::String;
        writer = NetCDFWriter(),
        interpol = nothing,
        nor = Inf,
    )

Create and return a `DiagnosticsGroup` containing a diagnostic that
dumps the spectrum at the specified
`interval` after the velocity fields have been interpolated, into NetCDF files prefixed by
`out_prefix`.
"""
function setup_dump_spectra_diagnostics(
    ::ClimateMachineConfigType,
    interval::String,
    out_prefix::String,
    nor::Float64;
    writer = NetCDFWriter(),
    interpol = nothing,
)
    @assert !isnothing(interpol)

    return DiagnosticsGroup(
        "Spectra",
        Diagnostics.dump_spectra_init,
        Diagnostics.dump_spectra_fini,
        Diagnostics.dump_spectra_collect,
        interval,
        out_prefix,
        writer,
        interpol,
        SpectraDiagnosticsParams(nor),
    )
end

function get_spectrum(mpicomm, mpirank, Q, bl, interpol, nor)
    FT = eltype(Q)
    istate = similar(Q.data, interpol.Npl, number_states(bl, Prognostic(), FT))
    interpolate_local!(interpol, Q.data, istate)
    all_state_data = accumulate_interpolated_data(mpicomm, interpol, istate)
    if mpirank == 0
        u = all_state_data[:, :, :, 2] ./ all_state_data[:, :, :, 1]
        v = all_state_data[:, :, :, 3] ./ all_state_data[:, :, :, 1]
        w = all_state_data[:, :, :, 4] ./ all_state_data[:, :, :, 1]
        x1 = Array(interpol.x1g)
        d = length(x1)
        s, k = power_spectrum(u, v, w, x1[d] - x1[1], d, nor)
        return s, k
    end
    return nothing, nothing
end

function dump_spectra_init(dgngrp, currtime)
    Q = Settings.Q
    bl = Settings.dg.balance_law
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)
    FT = eltype(Q)
    interpol = dgngrp.interpol
    nor = dgngrp.params.nor

    spectrum, wavenumber = get_spectrum(mpicomm, mpirank, Q, bl, interpol, nor)
    if mpirank == 0
        dims = OrderedDict("k" => (wavenumber, Dict()))
        vars = OrderedDict("spectrum" => (("k",), FT, Dict()))

        dprefix = @sprintf(
            "%s_%s-%s",
            dgngrp.out_prefix,
            dgngrp.name,
            Settings.starttime,
        )
        dfilename = joinpath(Settings.output_dir, dprefix)
        init_data(dgngrp.writer, dfilename, dims, vars)
    end

    return nothing
end

function dump_spectra_collect(dgngrp, currtime)
    Q = Settings.Q
    bl = Settings.dg.balance_law
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)
    FT = eltype(Q)
    interpol = dgngrp.interpol
    nor = dgngrp.params.nor

    spectrum, _ = get_spectrum(mpicomm, mpirank, Q, bl, interpol, nor)

    if mpirank == 0
        varvals = OrderedDict("spectrum" => spectrum)
        append_data(dgngrp.writer, varvals, currtime)
    end

    MPI.Barrier(mpicomm)
    return nothing
end

function dump_spectra_fini(dgngrp, currtime) end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
# set / get Node value

a_node = Tree.Node("some_value")
@test a_node.value == "some_value"

another_node = Tree.Node("different_value")
@test another_node.value == "different_value"

@test a_node.value == "some_value"

# set Node children & parent

parent_node = Tree.Node("parent")
child_node = Tree.Node("child_1")
Tree.add_children!(parent_node, [child_node])

@test length(parent_node.children) == 1

another_child_node = Tree.Node("child_2")
third_child_node = Tree.Node("child_3")
Tree.add_children!(parent_node, [another_child_node, third_child_node])

@test length(parent_node.children) == 3

# get Node children

@test parent_node.children == [child_node, another_child_node, third_child_node]

# get Node parent

@test Tree.get_parent(parent_node) == nothing
@test Tree.get_parent(child_node) == parent_node
@test Tree.get_parent(another_child_node) == parent_node
@test Tree.get_parent(third_child_node) == parent_node

# check if Node is leaf (no children)

@test Tree.is_leaf(child_node) == true
@test Tree.is_leaf(parent_node) == false

# check if Node is root (no parent)

@test Tree.is_root(parent_node) == true
@test Tree.is_root(child_node) == false

# check if Node is ancestor (above/parent) of other Node

grandchild_node = Tree.Node("grandchild")
Tree.add_children!(child_node, [grandchild_node])

@test Tree.is_ancestor(parent_node, child_node) == true
@test Tree.is_ancestor(parent_node, grandchild_node) == true
@test Tree.is_ancestor(child_node, grandchild_node) == true

@test Tree.is_ancestor(child_node, parent_node) == false
@test Tree.is_ancestor(grandchild_node, parent_node) == false
@test Tree.is_ancestor(grandchild_node, child_node) == false

# check if Node is descendant (below/child) of other Node

@test Tree.is_descendant(child_node, parent_node) == true
@test Tree.is_descendant(grandchild_node, parent_node) == true
@test Tree.is_descendant(grandchild_node, child_node) == true

@test Tree.is_descendant(parent_node, child_node) == false
@test Tree.is_descendant(parent_node, grandchild_node) == false
@test Tree.is_descendant(child_node, grandchild_node) == false

# check if Node is sibling (same parent) of other Node

@test Tree.is_sibling(child_node, another_child_node) == true
@test Tree.is_sibling(child_node, third_child_node) == true
@test Tree.is_sibling(another_child_node, third_child_node) == true

@test Tree.is_sibling(child_node, grandchild_node) == false
@test Tree.is_sibling(parent_node, grandchild_node) == false

grand_grandchild_node_1 = Tree.Node("grand_grandchild_node_1")
grand_grandchild_node_2 = Tree.Node("grand_grandchild_node_2")
Tree.add_children!(grandchild_node, [grand_grandchild_node_1, grand_grandchild_node_2])

grand_grand_grandchild_node_1 = Tree.Node("grand_grand_grandchild_node_1")
Tree.add_children!(grand_grandchild_node_1, [grand_grand_grandchild_node_1])

#                       parent
#     child_node, another_child_node, third_child_node
# grandchild_node
# ggc1      ggc2
# gggc1

# get depth (node to top)

@test Tree.depth(parent_node) == 0
@test Tree.depth(child_node) == 1
@test Tree.depth(another_child_node) == 1
@test Tree.depth(grandchild_node) == 2
@test Tree.depth(grand_grandchild_node_1) == 3
@test Tree.depth(grand_grandchild_node_2) == 3
@test Tree.depth(grand_grand_grandchild_node_1) == 4

# get height (node to bottom)

@test Tree.height(parent_node) == 4
@test Tree.height(child_node) == 3
@test Tree.height(another_child_node) == 0
@test Tree.height(third_child_node) == 0
@test Tree.height(grandchild_node) == 2
@test Tree.height(grand_grandchild_node_1) == 1
@test Tree.height(grand_grandchild_node_2) == 0
@test Tree.height(grand_grand_grandchild_node_1) == 0

# get root Node

@test Tree.get_root(parent_node) == parent_node
@test Tree.get_root(child_node) == parent_node
@test Tree.get_root(third_child_node) == parent_node
@test Tree.get_root(grand_grandchild_node_2) == parent_node
@test Tree.get_root(grand_grand_grandchild_node_1) == parent_node


text:
module Bitfinex

    using DataFrames

    import HTTP: get
    import JSON: parse
    import JSON2: read

    export getticker,
           getstats,
           getfundingbook,
           getorderbook,
           gettrades,
           getlends,
           getsymbols,
           getsymboldetails

    include("public.jl")
    include("utils.jl")

    const URL = "https://api.bitfinex.com/v1"

end


text:
module Chair

export chair, chair_system, ChairElem, ChairPTile

using Luxor
import Luxor.transform
using Colors
using StructEquality

using ...CoreDefs

@def_structequal struct ChairElem <: DGroupElem
    angle :: Int # should be 0 ≤ i ≤ 3
    x :: Int
    y :: Int
end

struct ChairPTile end

function chair(angle, x, y)
    return ChairElem(angle, x, y) => ChairPTile()
end

function sin(i :: Integer)
    if i == 0 || i == 2
        return 0
    elseif i == 1
        return 1
    else
        return -1
    end
end

function cos(i :: Integer)
    return sin((i+1) % 4 )
end

function Base.:*(s :: ChairElem, t :: ChairElem)
    return ChairElem(
        (s.angle + t.angle) % 4,
        s.x + cos(s.angle)*t.x + sin(s.angle)*t.y,
        s.y - sin(s.angle)*t.x + cos(s.angle)*t.y
    )
end
function Base.:*(s :: ChairElem, t :: Pair{Integer, Integer})
    return [
        s.x + cos(s.angle)*t[1] + sin(s.angle)*t[2],
        s.y - sin(s.angle)*t[1] + cos(s.angle)*t[2],
    ]
end

function CoreDefs.dilate(λ :: Int, t :: ChairElem)
    return ChairElem(t.angle, λ*t.x, λ*t.y)
end
function CoreDefs.id(::ChairElem)
    return ChairElem(0,0,0)
end
function Base.inv(g::ChairElem)
    return ChairElem(-g.angle,0,0)*ChairElem(0,-g.x, -g.y)
end

chair_subst = Dict([
    (ChairPTile(), [
        ChairElem(0, -1, -1) => ChairPTile(),
        ChairElem(0, 1, 1) => ChairPTile(),
        ChairElem(1, -1, 5) => ChairPTile(),
        ChairElem(3, 5, -1) => ChairPTile(),
    ])
])

chair_system = SubSystem(chair_subst, 2)

function CoreDefs.embed_aff(t :: ChairElem)
    return float([
        cos(t.angle), -sin(t.angle),
        sin(t.angle), cos(t.angle),
        t.x, t.y])
end

function CoreDefs.draw(::ChairPTile, action)
    Luxor.poly([
        Point(-1, -1),
        Point(-1, 3),
        Point(1, 3),
        Point(1, 1),
        Point(3, 1),
        Point(3, -1),
    ], close = true, action)
end

const origin = chair(0,0,0)

function in_bounds(t :: Pair{ChairElem, ChairPTile}, n, window)
    return 2^(n+1)*(abs(t[1].x)) < window.w + 2^(n+1)*3 && 2^(n+1)*(abs(t[1].y)) < window.h + 2^(n+1)*3
end

function squares(g)
    return g*[
        (0,0),
        (0,2),
        (2,0),
    ]
end
function adjacent(xs,ys)
    for x in xs
        for y in ys
            if max(abs.(x-y)...) ≤ 2
                return true
            end
        end
    end
    return false
end
function collar_in(tiling, g::ChairElem)
    result = Pair{ChairElem, ChairPTile}[]
    tiling_dict = Dict(tiling)
    squares_g = tiling_dict[g]
    for (h,l) in tiling_dict
        if adjacent(squares(g), squares(h))
            push!(result, (h,l))
        end
    end
    return result
end

end


text:
module TestKMeans

using ParallelKMeans
using ParallelKMeans: MultiThread
using Test
using Random

@testset "singlethread linear separation" begin
    Random.seed!(2020)

    X = rand(100, 3)
    labels, centroids, sum_squares = kmeans(X, 3; tol = 1e-10, verbose = false)

    @test sum_squares ≈ 14.964882850452971
end

@testset "multithread linear separation" begin
    Random.seed!(2020)

    X = rand(100, 3)
    labels, centroids, sum_squares = kmeans(X, 3, MultiThread(); tol = 1e-10, verbose = false)

    @test sum_squares ≈ 14.964882850452971
end

end # module


text:
##
# Implements contextual dispatch through Cassette.jl
# Goals:
# - Rewrite common CPU functions to appropriate GPU intrinsics
#
# TODO:
# - error (erf, ...)
# - min, max
# - mod, rem
# - gamma
# - bessel
# - distributions
# - unsorted

using Cassette

##
# Important for inference to not be able to see the constant value
##
@inline function unknowably_false()
    Base.llvmcall("ret i8 0", Bool, Tuple{})
end

const INTERACTIVE = haskey(ENV, "GPUIFYLOOPS_INTERACTIVE") && ENV["GPUIFYLOOPS_INTERACTIVE"] == "1"

##
# Forces inlining on everything that is not marked `@noinline`
# Cassette has a #265 issue, let's try to work around that.
##
function transform(ctx, ref)
    CI = ref.code_info
    noinline = any(@nospecialize(x) ->
                       Core.Compiler.isexpr(x, :meta) &&
                       x.args[1] == :noinline,
                   CI.code)
    CI.inlineable = !noinline

    @static if INTERACTIVE
    # 265 fix, insert a call to the original method
    # that we later will remove with LLVM's DCE
    unknowably_false = GlobalRef(@__MODULE__, :unknowably_false)
    Cassette.insert_statements!(CI.code, CI.codelocs,
      (x, i) -> i == 1 ?  4 : nothing,
      (x, i) -> i == 1 ? [
          Expr(:call, Expr(:nooverdub, unknowably_false)),
          Expr(:gotoifnot, Core.SSAValue(i), i+3),
          Expr(:call, Expr(:nooverdub, Core.SlotNumber(1)), (Core.SlotNumber(i) for i in 2:ref.method.nargs)...),
          x] : nothing)
    end
    CI.ssavaluetypes = length(CI.code)
    # Core.Compiler.validate_code(CI)
    return CI
end

const GPUifyPass = Cassette.@pass transform

Cassette.@context Ctx
const ctx = Cassette.disablehooks(Ctx(pass = GPUifyPass))

###
# Cassette fixes
###
@inline Cassette.overdub(::Ctx, ::typeof(Core.kwfunc), f) = return Core.kwfunc(f)
@inline Cassette.overdub(::Ctx, ::typeof(Core.apply_type), args...) = return Core.apply_type(args...)

# the functions below are marked `@pure` and by rewritting them we hide that from
# inference so we leave them alone (see https://github.com/jrevels/Cassette.jl/issues/108).
@inline Cassette.overdub(::Ctx, ::typeof(Base.isimmutable), x) = return Base.isimmutable(x)
@inline Cassette.overdub(::Ctx, ::typeof(Base.isstructtype), t) = return Base.isstructtype(t)
@inline Cassette.overdub(::Ctx, ::typeof(Base.isprimitivetype), t) = return Base.isprimitivetype(t)
@inline Cassette.overdub(::Ctx, ::typeof(Base.isbitstype), t) = return Base.isbitstype(t)
@inline Cassette.overdub(::Ctx, ::typeof(Base.isbits), x) = return Base.isbits(x)
@inline Cassette.overdub(::Ctx, ::typeof(StaticArrays.Size), x::Type{<:AbstractArray{<:Any, N}}) where {N} = return StaticArrays.Size(x)

@init @require CUDAnative="be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
    using .CUDAnative
    @inline Cassette.overdub(::Ctx, ::typeof(CUDAnative.datatype_align), ::Type{T}) where {T} = CUDAnative.datatype_align(T)
end

###
# Rewrite functions
###

# define +, -, * as contract

for (f, T) in Base.Iterators.product((:add, :mul, :sub), (Float32, Float64))
    name = Symbol("$(f)_float_contract")
    if T === Float32
        llvmt = "float"
    elseif T === Float64
        llvmt = "double"
    end

    #XXX Use LLVM.jl
    ir = """
        %x = f$f contract $llvmt %0, %1
        ret $llvmt %x
    """
    @eval begin
        # the @pure is necessary so that we can constant propagate.
        Base.@pure function $name(a::$T, b::$T)
            @Base._inline_meta
            Base.llvmcall($ir, $T, Tuple{$T, $T}, a, b)
        end
    end
end
@inline Cassette.overdub(ctx::Ctx, ::typeof(+), a::T, b::T) where T<:Union{Float32, Float64} = add_float_contract(a, b)
@inline Cassette.overdub(ctx::Ctx, ::typeof(-), a::T, b::T) where T<:Union{Float32, Float64} = sub_float_contract(a, b)
@inline Cassette.overdub(ctx::Ctx, ::typeof(*), a::T, b::T) where T<:Union{Float32, Float64} = mul_float_contract(a, b)
@inline Cassette.overdub(ctx::Ctx, ::typeof(^), x::Float64, y::Float64) = CUDAnative.pow(x, y)
@inline Cassette.overdub(ctx::Ctx, ::typeof(^), x::Float32, y::Float32) = CUDAnative.pow(x, y)
@inline Cassette.overdub(ctx::Ctx, ::typeof(^), x::Float64, y::Int32)   = CUDAnative.pow(x, y)
@inline Cassette.overdub(ctx::Ctx, ::typeof(^), x::Float32, y::Int32)   = CUDAnative.pow(x, y)
@inline Cassette.overdub(ctx::Ctx, ::typeof(^), x::Union{Float32, Float64}, y::Int64) = CUDAnative.pow(x, y)

# libdevice.jl
const cudafuns = (:cos, :cospi, :sin, :sinpi, :tan,
          :acos, :asin, :atan,
          :cosh, :sinh, :tanh,
          :acosh, :asinh, :atanh,
          :log, :log10, :log1p, :log2,
          :exp, :exp2, :exp10, :expm1, :ldexp,
          :isfinite, :isinf, :isnan,
          :signbit, :abs,
          :sqrt, :cbrt,
          :ceil, :floor,)
for f in cudafuns
    @eval function Cassette.overdub(ctx::Ctx, ::typeof(Base.$f), x::Union{Float32, Float64})
        @Base._inline_meta
        return CUDAnative.$f(x)
    end
end

"""
    contextualize(::Dev, f)

This contexualizes the function `f` for a given device type `Dev`.

For the device `CUDA()`, `contextualize` replaces calls to math library
functions.  For example, `cos`, `sin`, are replaced with `CUDAnative.cos`,
`CUDAnative.sin`, respectively.

The full list functions that are replaced is $cudafuns.

# Examples
```julia
function kernel!(::Dev, A, f) where {Dev}
    @setup Dev
    @loop for i in (1:size(A,1); threadIdx().x)
        A[i] = f(A[i])
    end
end

g(x) = sin(x)
kernel!(A::Array) = kernel!(CPU(), A, contextualize(CPU(), g))
kernel!(A::CuArray) =
    @cuda threads=length(A) kernel!(CUDA(), A, contextualize(CUDA(), g))

a = rand(Float32, 1024)
b, c = copy(a), CuArray(a)

kernel!(b)
kernel!(c)

@assert g.(a) ≈ b
@assert g.(a) ≈ c
```
"""
contextualize(f::F) where F = (args...) -> Cassette.overdub(ctx, f, args...)


text:
using WiagDataSetup
using Test

@testset "WiagDataSetup.jl" begin
    # Write your tests here.
end

max_stars_repo_name: ghertko/WiagDataSetup

max_issues_repo_name: ghertko/WiagDataSetup

max_forks_repo_name: ghertko/WiagDataSetup


text:
# abstract type AbstractLazyMul end
struct LazyMulAdd{M,O,T<:NativeTypesV} <: Number
    data::T
    @inline LazyMulAdd{M,O,T}(data::T) where {M,O,T<:NativeTypesV} = new{M,O,T}(data)
end
# O for offset is kind of hard to read next to default of 0?
@inline LazyMulAdd{M,O}(data::T) where {M,O,T<:Union{Base.HWReal,AbstractSIMD}} = LazyMulAdd{M,O,T}(data)
@inline LazyMulAdd{M}(data::T) where {M,T<:NativeTypesV} = LazyMulAdd{M,0,T}(data)
@inline LazyMulAdd{0,O}(data::T) where {O,T<:Union{Base.HWReal,AbstractSIMD}} = @assert false#StaticInt{O}()
@inline LazyMulAdd{0}(data::T) where {T<:NativeTypesV} = @assert false#StaticInt{0}()
@inline LazyMulAdd(data::T, ::StaticInt{M}) where {M,T} = LazyMulAdd{M,0,T}(data)
@inline LazyMulAdd(data::T, ::StaticInt{M}, ::StaticInt{O}) where {M,O,T} = LazyMulAdd{M,O,T}(data)

@inline data(lm::LazyMulAdd) = data(lm.data) # for use with indexing

@inline _materialize(a::LazyMulAdd{M,O,I}) where {M,O,I} = vadd(vmul(StaticInt{M}(), a.data), StaticInt{O}())
@inline Base.convert(::Type{T}, a::LazyMulAdd{M,O,I}) where {M,O,I,T<:Number} = convert(T, _materialize(a))

Base.promote_rule(::Type{LazyMulAdd{M,O,I}}, ::Type{T}) where {M,O,I<:Number,T} = promote_type(I, T)
Base.promote_rule(::Type{LazyMulAdd{M,O,MM{W,X,I}}}, ::Type{T}) where {M,O,W,X,I,T} = promote_type(Vec{W,I}, T)
Base.promote_rule(::Type{LazyMulAdd{M,O,Vec{W,I}}}, ::Type{T}) where {M,O,W,I,T} = promote_type(Vec{W,I}, T)

@inline lazymul(a, b) = vmul(a, b)
@inline lazymul(::StaticInt{M}, b) where {M} = LazyMulAdd{M}(b)
@inline lazymul(::StaticInt{1}, b) = b
@inline lazymul(::StaticInt{0}, b) = StaticInt{0}()
@inline lazymul(a, ::StaticInt{M}) where {M} = LazyMulAdd{M}(a)
@inline lazymul(a, ::StaticInt{1}) = a
@inline lazymul(a, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul(a::LazyMulAdd, ::StaticInt{1}) = a
@inline lazymul(::StaticInt{1}, a::LazyMulAdd) = a
@inline lazymul(a::LazyMulAdd, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul(::StaticInt{0}, a::LazyMulAdd) = StaticInt{0}()
# @inline lazymul(a::LazyMulAdd, ::StaticInt{0}) = StaticInt{0}()
# @inline lazymul(::StaticInt{M}, b::MM{W,X}) where {W,M,X} = LazyMulAdd{M}(MM{W}(b.data, StaticInt{M}()*StaticInt{X}()))
# @inline lazymul(a::MM{W,X}, ::StaticInt{M}) where {W,M,X} = LazyMulAdd{M}(MM{W}(a.data, StaticInt{M}()*StaticInt{X}()))
@inline lazymul(::StaticInt{M}, b::MM{W,X}) where {W,M,X} = MM{W}(vmul(StaticInt{M}(), data(b)), StaticInt{X}() * StaticInt{M}())
@inline lazymul(a::MM{W,X}, ::StaticInt{M}) where {W,M,X} = MM{W}(vmul(StaticInt{M}(), data(a)), StaticInt{X}() * StaticInt{M}())
@inline lazymul(a::MM{W,X}, ::StaticInt{0}) where {W,X} = StaticInt{0}()
@inline lazymul(::StaticInt{0}, a::MM{W,X}) where {W,X} = StaticInt{0}()
@inline lazymul(a::MM{W,X}, ::StaticInt{1}) where {W,X} = a
@inline lazymul(::StaticInt{1}, a::MM{W,X}) where {W,X} = a

@inline lazymul_no_promote(::Type{T}, a, b) where {T} = lazymul_no_promote(a, b)
# @inline lazymul_no_promote(::Type{T}, a::MM{W,X}, b::StaticInt) where {W,X,T} = MM{W}(vmul(a.i, b), mulsizeof(T, StaticInt{X}()))
@inline lazymul_no_promote(::Type{T}, a::MM{W,X}, b::StaticInt{M}) where {M,W,X,T} = LazyMulAdd{M}(a)
@inline lazymul_no_promote(::Type{T}, a::MM{W,X}, b::Integer) where {W,X,T} = MM{W}(vmul(a.i, b), mulsizeof(T, StaticInt{X}()))
@inline lazymul_no_promote(::Type{T}, b::StaticInt, a::MM{W,X}) where {W,X,T} = MM{W}(vmul(b, a.i), mulsizeof(T, StaticInt{X}()))
@inline lazymul_no_promote(::Type{T}, b::Integer, a::MM{W,X}) where {W,X,T} = MM{W}(vmul(b, a.i), mulsizeof(T, StaticInt{X}()))

# @inline lazymul(::StaticInt{M}, b::MM{W,X,StaticInt{N}}) where {W,M,X,N} = MM{W}(StaticInt{N}() * StaticInt{M}(), StaticInt{X}() * StaticInt{M}())
# @inline lazymul(a::MM{W,X,StaticInt{N}}, ::StaticInt{M}) where {W,M,X,N} = MM{W}(StaticInt{N}() * StaticInt{M}(), StaticInt{X}() * StaticInt{M}())

# @inline lazymul(::StaticInt{N}, b::MM) where {N} = LazyMulAdd{N}(Vec(b))
# @inline lazymul(a::MM, ::StaticInt{N}) where {N} = LazyMulAdd{N}(Vec(a))
@inline lazymul(::StaticInt{M}, ::StaticInt{N}) where {M, N} = StaticInt{M}()*StaticInt{N}()
@inline lazymul(::StaticInt{0}, ::StaticInt) = StaticInt{0}()
@inline lazymul(::StaticInt{1}, ::StaticInt{M}) where {M} = StaticInt{M}()
@inline lazymul(::StaticInt, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul(::StaticInt{M}, ::StaticInt{1}) where {M} = StaticInt{M}()
@inline lazymul(::StaticInt{0}, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul(::StaticInt{0}, ::StaticInt{1}) = StaticInt{0}()
@inline lazymul(::StaticInt{1}, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul(::StaticInt{1}, ::StaticInt{1}) = StaticInt{1}()

@inline lazymul(a::LazyMulAdd{M,O}, ::StaticInt{N}) where {M,O,N} = LazyMulAdd(a.data, StaticInt{M}() * StaticInt{N}(), StaticInt{O}() * StaticInt{N}())
@inline lazymul(::StaticInt{M}, b::LazyMulAdd{N,O}) where {M,O,N} = LazyMulAdd(b.data, StaticInt{M}() * StaticInt{N}(), StaticInt{O}() * StaticInt{M}())

@inline lazymul(a::LazyMulAdd{M,<:MM{W,X}}, ::StaticInt{N}) where {M,N,W,X} = LazyMulAdd(MM{W}(a.data, StaticInt{M}()*StaticInt{N}()*StaticInt{X}()), StaticInt{M}()*StaticInt{N}())
@inline lazymul(::StaticInt{M}, b::LazyMulAdd{N,<:MM{W,X}}) where {M,N,W,X} = LazyMulAdd(MM{W}(b.data, StaticInt{M}()*StaticInt{N}()*StaticInt{X}()), StaticInt{M}()*StaticInt{N}())
@inline lazymul(a::LazyMulAdd{M,<:MM{W,X}}, ::StaticInt{0}) where {M,W,X} = StaticInt{0}()
@inline lazymul(::StaticInt{0}, b::LazyMulAdd{N,<:MM{W,X}}) where {N,W,X} = StaticInt{0}()
@inline lazymul(a::LazyMulAdd{M,<:MM{W,X}}, ::StaticInt{1}) where {M,W,X} = a
@inline lazymul(::StaticInt{1}, b::LazyMulAdd{N,<:MM{W,X}}) where {N,W,X} = b
@inline lazymul(a::LazyMulAdd{M}, b::LazyMulAdd{N}) where {M,N} = LazyMulAdd(vmul(a.data, b.data), StaticInt{M}()*StaticInt{N}())

@inline Base.:(>>>)(a::LazyMulAdd{M,O}, ::StaticInt{N}) where {M,O,N} = LazyMulAdd(a.data, StaticInt{M}() >>> StaticInt{N}(), StaticInt{O}() >>> StaticInt{N}())

@inline lazymul_no_promote(a, b) = vmul_no_promote(a, b)
@inline lazymul_no_promote(::StaticInt{N}, b) where {N} = LazyMulAdd{N}(b)
@inline lazymul_no_promote(a, ::StaticInt{N}) where {N} = LazyMulAdd{N}(a)
@inline lazymul_no_promote(::StaticInt{N}, b::MM{W,X,StaticInt{M}}) where {N,W,X,M} = MM{W,X}(StaticInt{M}() * StaticInt{N}())
@inline lazymul_no_promote(a::MM{W,X,StaticInt{M}}, ::StaticInt{N}) where {N,W,X,M} = MM{W,X}(StaticInt{M}() * StaticInt{N}())
@inline lazymul_no_promote(::StaticInt{M}, ::StaticInt{N}) where {M, N} = StaticInt{M}()*StaticInt{N}()
@inline lazymul_no_promote(::StaticInt{1}, b::MM{W,X,StaticInt{M}}) where {W,X,M} = b
@inline lazymul_no_promote(a::MM{W,X,StaticInt{M}}, ::StaticInt{1}) where {W,X,M} = a

@inline lazymul_no_promote(::StaticInt{1}, b) = b
@inline lazymul_no_promote(a, ::StaticInt{1}) = a
@inline lazymul_no_promote(::StaticInt{1}, ::StaticInt{N}) where {N} = StaticInt{N}()
@inline lazymul_no_promote(::StaticInt{M}, ::StaticInt{1}) where {M} = StaticInt{M}()
@inline lazymul_no_promote(::StaticInt{1}, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul_no_promote(::StaticInt{0}, ::StaticInt{1}) = StaticInt{0}()
@inline lazymul_no_promote(::StaticInt{1}, ::StaticInt{1}) = StaticInt{1}()

@inline lazymul_no_promote(a::LazyMulAdd{M}, ::StaticInt{N}) where {M,N} = LazyMulAdd(a.data, StaticInt{M}()*StaticInt{N}())
@inline lazymul_no_promote(::StaticInt{M}, b::LazyMulAdd{N}) where {M,N} = LazyMulAdd(b.data, StaticInt{M}()*StaticInt{N}())
@inline lazymul_no_promote(a::LazyMulAdd{M}, b::LazyMulAdd{N}) where {M,N} = LazyMulAdd(vmul(a.data, b.data), StaticInt{M}()*StaticInt{N}())
@inline lazymul_no_promote(a::LazyMulAdd{M}, ::One) where {M} = a
@inline lazymul_no_promote(::One, b::LazyMulAdd{N}) where {N} = b

@inline lazymul_no_promote(::StaticInt{0}, b) = StaticInt{0}()
@inline lazymul_no_promote(a, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul_no_promote(::StaticInt{0}, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul_no_promote(::StaticInt{0}, ::StaticInt{M}) where {M} = StaticInt{0}()
@inline lazymul_no_promote(::StaticInt{M}, ::StaticInt{0}) where {M} = StaticInt{0}()
@inline lazymul_no_promote(::StaticInt{0}, b::MM{W,X,StaticInt{M}}) where {W,X,M} = StaticInt{0}()
@inline lazymul_no_promote(a::MM{W,X,StaticInt{M}}, ::StaticInt{0}) where {W,X,M} = StaticInt{0}()
@inline lazymul_no_promote(a::LazyMulAdd{M}, ::StaticInt{0}) where {M} = StaticInt{0}()
@inline lazymul_no_promote(::StaticInt{0}, b::LazyMulAdd{N}) where {N} = StaticInt{0}()

# @inline lazyadd(a, b) = vadd(a, b)
@inline vadd(a::LazyMulAdd{M,O,T}, ::StaticInt{A}) where {M,O,T,A} = LazyMulAdd(a.data, StaticInt{M}(), StaticInt{O}()+StaticInt{A}())
@inline vadd(::StaticInt{A}, a::LazyMulAdd{M,O,T}) where {M,O,T,A} = LazyMulAdd(a.data, StaticInt{M}(), StaticInt{O}()+StaticInt{A}())
@inline vadd(a::LazyMulAdd{M,O,T}, ::StaticInt{0}) where {M,O,T} = a
@inline vadd(::StaticInt{0}, a::LazyMulAdd{M,O,T}) where {M,O,T} = a
# @inline vadd(a::LazyMulAdd{M,O,T}, ::StaticInt{A}) where {M,O,T<:MM,A} = LazyMulAdd(a.data, StaticInt{M}(), StaticInt{O}()+StaticInt{A}())
# @inline vadd(::StaticInt{A}, a::LazyMulAdd{M,O,T}) where {M,O,T,A} = LazyMulAdd(a.data, StaticInt{M}(), StaticInt{O}()+StaticInt{A}())

@inline vadd(a::LazyMulAdd{M,O,MM{W,X,I}}, b::Integer) where {M,O,W,X,I} = MM{W}(vadd(vmul(StaticInt{M}(), data(a)), vadd(StaticInt{0}(), b)), StaticInt{X}() * StaticInt{M}())
@inline vadd(::StaticInt{N}, a::LazyMulAdd{M,O,MM{W,X,I}}) where {N,M,O,W,X,I} = LazyMulAdd(a.data, StaticInt{M}(), StaticInt{O}()+StaticInt{N}())
@inline vadd(a::LazyMulAdd{M,O,MM{W,X,I}}, ::StaticInt{N}) where {N,M,O,W,X,I} = LazyMulAdd(a.data, StaticInt{M}(), StaticInt{O}()+StaticInt{N}())
@inline vadd(b::Integer, a::LazyMulAdd{M,O,MM{W,X,I}}) where {M,O,W,X,I} = MM{W}(vadd(vmul(StaticInt{M}(), data(a)), vadd(StaticInt{0}(), b)), StaticInt{X}() * StaticInt{M}())
@inline vadd(a::LazyMulAdd{M,O,MM{W,X,I}}, ::StaticInt{0}) where {M,O,W,X,I} = a
@inline vadd(::StaticInt{0}, a::LazyMulAdd{M,O,MM{W,X,I}}) where {M,O,W,X,I} = a
# @inline vadd(::StaticInt{M}, a::LazyMulAdd{M,O,MM{W,X,I}}) where {M,O,W,X,I} = a

@inline vadd(a::LazyMulAdd{M,O}, b::LazyMulAdd{M,A}) where {M,O,A} = LazyMulAdd(vadd(a.data, b.data), StaticInt{M}(), StaticInt{O}()+StaticInt{A}())

@inline vadd(a::LazyMulAdd{M,O,MM{W,X,I}}, b::LazyMulAdd{N,P,J}) where {M,O,W,X,I,N,P,J<:IntegerTypes} = vadd(a, _materialize(b))
@inline vadd(b::LazyMulAdd{N,P,J}, a::LazyMulAdd{M,O,MM{W,X,I}}) where {M,O,W,X,I,N,P,J<:IntegerTypes} = vadd(a, _materialize(b))
@inline vadd(a::LazyMulAdd{M,O,MM{W,X,I}}, b::LazyMulAdd{M,A,J}) where {M,O,W,X,I,A,J<:IntegerTypes} = LazyMulAdd(vadd(a.data, b.data), StaticInt{M}(), StaticInt{O}()+StaticInt{A}())
@inline vadd(b::LazyMulAdd{M,A,J}, a::LazyMulAdd{M,O,MM{W,X,I}}) where {M,O,W,X,I,A,J<:IntegerTypes} = LazyMulAdd(vadd(a.data, b.data), StaticInt{M}(), StaticInt{O}()+StaticInt{A}())

@inline vadd(::MM{W,X,StaticInt{A}}, a::LazyMulAdd{M,O,T}) where {M,O,T<:Integer,A,W,X} = LazyMulAdd(MM{W,X}(a.data), StaticInt{M}(), StaticInt{O}()+StaticInt{A}())
@inline vadd(a::LazyMulAdd{M,O,T}, ::MM{W,X,StaticInt{A}}) where {M,O,T<:Integer,A,W,X} = LazyMulAdd(MM{W,X}(a.data), StaticInt{M}(), StaticInt{O}()+StaticInt{A}())
@inline vadd(a::MM{W,X,I}, b::LazyMulAdd{N,P,J}) where {W,X,I<:IntegerTypesHW,N,P,J<:Integer} = vadd(a, _materialize(b))
@inline vadd(b::LazyMulAdd{N,P,J}, a::MM{W,X,I}) where {W,X,I<:IntegerTypesHW,N,P,J<:Integer} = vadd(a, _materialize(b))


text:
module day21

using StatsBase: countmap

include("inputs.jl")

ex = inputs.ex_d21()
input = inputs.input_d21()

function parse_input(input)
    foods = split(input, "\n")

    return map(foods) do food
        rgx = r"(.*) \(contains (.*)\)"
        ing_str, alg_str = match(rgx, food).captures
        ingredients = String.(split(ing_str |> strip, " "))
        allergens = String.(split(alg_str |> strip, ", "))
        (; ingredients, allergens)
    end
end

# candidate dict mapping allergens to possible ingredients
function find_candidates(foods)
    candidates = Dict{String, Set{String}}()
    for food in foods
        for alg in food.allergens
            if haskey(candidates, alg)
                # push!(candidates[alg], (food.ingredients)...)
                candidates[alg] = intersect(candidates[alg], Set(food.ingredients))
            else
                candidates[alg] = Set(food.ingredients)
            end
        end
    end
    return candidates
end

# resolve candidates by mapping allergen to exactly one ingredient
function resolve_candidates(candidates)
    result = Dict()
    queue = collect(keys(candidates))

    while length(result) < length(candidates)
        allergen = popfirst!(queue)
        if length(candidates[allergen]) == 1
            ingredient = only(candidates[allergen])
            result[allergen] = ingredient
            for k in keys(candidates)
                if k != allergen
                    delete!(candidates[k], ingredient)
                end
            end
        else
            push!(queue, allergen)
        end
    end
    return result
end

# count non-allergen ingredient occurences
function count_safe_ingredients(foods, result)
    safe = []
    for food in foods, ingredient in food.ingredients
            if ingredient ∉ collect(values(result))
                push!(safe, ingredient)
            end
    end
    sum(values(countmap(safe)))
end

function part1(input)
    foods = parse_input(input)
    candidates = find_candidates(foods)
    result = resolve_candidates(candidates)
    count_safe_ingredients(foods, result)
end

@assert part1(ex) == 5
@info "Day 21, Part 1 answer: $(part1(input))"

function part2(input)
    foods = parse_input(input)
    candidates = find_candidates(foods)
    result = resolve_candidates(candidates)
    join([result[k] for k in sort(collect(keys(result)))], ",")
end

@assert part2(ex) == "mxmxvkd,sqjhc,fvjkl"
@assert part2(input) == "ppdplc,gkcplx,ktlh,msfmt,dqsbql,mvqkdj,ggsz,hbhsx"
@info "Day 21, Part 2 answer: $(part2(input))"

end


text:
using LinearAlgebra, Arpack, IterativeSolvers
using Random, Distributions, Distances
using Zygote
using Zygote: @adjoint
import Base: display, size
import LinearAlgebra: logabsdet, logdet, eltype, mul!

include("./kernel.jl")
include("./transform.jl")
include("./covariate.jl")
include("./opt.jl")

struct Data{T}
    X::Matrix{T} # d x n
    Fu::LU{Float64, Matrix{T}} # p x p
    W::Matrix{T} # p x n - p
    y::Vector{T} # n
end

function display(d::Data)
    println("    X: Observation Locations")
    display(d.X)
    println("\n    Fu: Covariate Basis")
    display(d.Fu)
    println("\n    W: Transformed Covariates")
    display(d.W)
    println("\n    y: Observations")
    display(d.y)
    return nothing
end

function data(f, X, y)
    Fx = covariate(f, X)
    p, n = size(Fx)
    Fu = lu(Fx[:, 1:p])
    W = Fu \ Fx[:, p+1:end]
    return Data(X, Fu, W, y)
end

mutable struct KernelSystem{
    T, K<:AbstractKernel, G<:AbstractTransform, F<:AbstractCovariate
}
    f::F
    k::K
    g::G
    d::Data{T}
    l::LengthScale{T}
    km::KernelMatrix{T}
    ϵ::T
    o::Output{T}
end

function display(ks::KernelSystem)
    println("        f: Covariate Function")
    display(ks.f)
    println("\n        k: Kernel")
    display(ks.k)
    println("\n        g: Transformation")
    display(ks.g)
    println("\n        d: Data")
    display(ks.d)
    println("\n        l: Scaled Distances")
    display(ks.l)
    println("\n        km: Kernel Matrix")
    display(ks.km)
    println("\n        ϵ: Nugget term")
    display(ks.ϵ)
    println("\n        o: Transformed Observations")
    display(ks.o)
    return nothing
end

function kernelsystem(f, k, g, d, ℓ, θ, ϵ, λ)
    p, m = size(d.W)
    l = lengthscale(distance(k), p, ℓ, d.X)
    km = kernelmatrix(k, θ, l)
    o = output(d.W, g, λ, d.y, p)
    return KernelSystem(f, k, g, d, l, km, ϵ, o)
end
schur(ks::KernelSystem) = Schur(ks.d.W, ks.km.Ku, ks.km.Kux, ks.km.Kx, ks.ϵ)
function logabsdet(ks::KernelSystem)
    s = (-1.0) ^ (size(ks.d.Fu, 1) % 2)
    return 2 * logabsdet(ks.d.Fu)[1] * logdet(schur(ks))
end
function logabsdet(ks::KernelSystem, rKx)
    s = (-1.0) ^ (size(ks.d.Fu, 1) % 2)
    return 2 * logabsdet(ks.d.Fu)[1] * logdet(rKx)
end

struct Schur{T}
    W::Matrix{T}
    Ku::Matrix{T}
    Kux::Matrix{T}
    Kx::Matrix{T}
    ϵ::T
end

function Matrix(sc::Schur)
    B = sc.Kux .- sc.Ku * sc.W .- sc.ϵ .* sc.W
    return sc.Kx .- sc.W' * B .- sc.Kux' * sc.W + sc.ϵ * I
end
eltype(::Schur{T}) where T = T
size(sc::Schur, args...) = size(sc.Kx, args...)
function mul!(y, sc::Schur, x, α, β)
    # (Kx - Kux'W - W'Kux + W'(Ku + ϵI)W + ϵI)x * α + y * β
    p, m = size(sc.W)
    v1 = similar(x, p)
    v2 = similar(x, p)

    y .= y .* β .+ x .* (α * sc.ϵ)

    mul!(y, sc.Kx, x, α, 1)

    mul!(v1, sc.W, x)
    mul!(v2, v1, ϵ)
    mul!(v2, sc.Ku, v1, 1, 1)

    mul!(y, sc.W', v2, α, 1)
    mul!(y, sc.Kux', v1, -α, 1)

    mul!(v1, sc.Kux, x)

    mul!(y, sc.W', v1, -α, 1)

    return y
end

function lanczos_arpack(A, k, v; maxiter, tol)
    T = eltype(A)
    n = size(A, 1)
    mulA! = (y, x) -> mul!(y, A, x)
    id = x -> x
    # in: (T, mulA!, mulB, solveSI, n, issym, iscmplx, bmat,
    #            nev, ncv, whichstr, tol, maxiter, mode, v0)
    # out: (resid, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, TOL)
    out = Arpack.aupd_wrapper(T, mulA!, id, id, n, true, false, "I",
                       1, k, "LM", tol, maxiter, 1, v)

    α = out[7][k+1:2*k-1]
    β = out[7][2:k-1]

    return out[2], α, β, out[1]
end

function _lanczos_logdet!(z, acc, sc, k; maxiter, tol, nsamples)
    for i in 1:nsamples
        rand!(Normal(), z)
        z .= sign.(z)
        Q, α, β, resid = lanczos_arpack(sc, k, z; maxiter=maxiter, tol=tol)
        T = SymTridiagonal(α, β)
        Λ = eigen(T)
        wts = Λ.vectors[1, :] .^ 2 .* norm(z) ^ 2
        acc += dot(wts, log.(Λ.values))
    end
    return acc / nsamples
end

function lanczos_logdet(A, k; maxiter, tol, nsamples)
    z = Vector{eltype(A)}(undef, size(A, 1))
    return _lanczos_logdet!(z, 0, A, k;
                            maxiter=maxiter, tol=tol, nsamples=nsamples)
end

logdet(sc::Schur) = lanczos_logdet(sc, max(20, size(sc, 1));
                                   maxiter=1000, tol=0.0, nsamples=50)

function logprob(ks, λ, logprior)
    p, m = size(ks.d.W)
    lnJ = logjacobian(ks.g, λ, ks.o.z)
    rKx = cholesky(Symmetric(Matrix(schur(ks))))
    lnq = dot(ks.o.rz, rKx \ ks.o.rz)
    lnDet = logabsdet(ks, rKx)[1]
    return m / (m + p) * lnJ - m / 2 * lnq - lnDet / 2 + logprior
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
##########
##########
# Simulate to get equilibrium abundances for the webs
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-06-25
##########
##########

# set up =======================================================================

##include("./src/eqn_component_functions.jl")
#include("./src/species_list_creation.jl")
#include("./src/network_metrics.jl")
using DifferentialEquations
using DelimitedFiles

# functions ====================================================================

"""

    biomass(dx, x, p, t)

Calling biomass functions, create the actual equation for each species, to be
passed to the simulation.
"""
function biomass_basic(du, u, p, t)
    species_list, network = p[1:2]
    for i in 1:size(species_list, 1)
        du[i] = growth_rate(i, species_list, u) +
                biomass_gain(species_list, network, u, 2, i) -
                biomass_loss(species_list, network, u, 2, i) -
                mortality_change(i, u)
    end
end

"""
    starting_eq_sims(species_list, network, initial_sim_length)

Function runs through all initial simulations from t=0 to the point of starting
measurement.
"""
function starting_eq_sims!(
    species_list,
    network,
    initial_sim_length,
    network_id
    )

    #initials = species_list, network

    # cut network to size
    network = network[:,2:end]

    #  allocate x0
    x0 = Float64[]
    for i in 1:size(species_list, 1)
        append!(x0, species_list[i, biomass_pos])
    end

    # keep track of how many extinctions there were
    extinct_sp = Float64[]

    # loop through all time values from 0 to 5000
    time_start = 0.0
    while time_start < initial_sim_length

        # put species list and network in same object
        p = [species_list, network]

        # set time and defie ODE problem
        tspan = (time_start, time_start+9.9)
        prob = ODEProblem(biomass!, x0, tspan, p)
        sol = solve(prob)

        # keep only last timestep and check for extinction
        x0 = sol.u[end]
        check_extinctions!(x0)
        species_list[:,biomass_pos] = x0

        # get new network
        network = create_interaction_matrix(species_list)

        # change alive to extinct if species is extinct
        for i in 1:size(species_list, 1)
            if species_list[i, biomass_pos] == 0
                species_list[i, alive_pos] = 0
                append!(extinct_sp, species_list[i, organism_id_pos])
            end
        end

        time_start += 10.0
    end

    # bind back the network_id and write the species list
    species_list = hcat(species_list, fill(network_id, 50, 1))
    writedlm("./data/intermediate-files/species-lists-pre-perturbation/species_list_$network_id.csv", species_list, ",")

    # create and write network
    network = create_interaction_matrix(species_list)
    network = hcat(fill(network_id, 50, 1), network)
    writedlm("./data/intermediate-files/networks-pre-perturbation/network_$network_id.csv", network, ",")

end

"""
    extinction_cond(out, x, t, integrator)

Condition asking if any of the values in the integrator go below the extinction
threshold.
"""
function extinction_threshold_pre_perturb(out,u,t,integrator)
    # loop through all species to make the condition check all of them
    for i in 1:size(species_list, 1)
        out[i] = 1e-06 - u[i]
    end
end

"""
    extinction_affect!(integrator, event_idx)

Affect function to use in callback vector, this will allow us to check
if any of the values of the threshold condition are true, and if so,
set the value to zero in the integrator.
"""
function extinction_affect_pre_perturb!(integrator, event_idx)
    # loop again through all species
    for i in 1:size(species_list, 1)
        if event_idx == i
            integrator.u[i] = 0
        end
    end
end

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
export AbstractAsyncInstrument,
    AbstractSyncInstrument,
    Measurement,
    Counter,
    ObservableCounter,
    Histogram,
    ObservableGauge,
    UpDownCounter,
    ObservableUpDownCounter

struct Measurement{V,T<:StaticAttrs}
    value::V
    attributes::T
end

Measurement(v) = Measurement(v, StaticAttrs())

abstract type AbstractSyncInstrument{T} <: AbstractInstrument{T} end
abstract type AbstractAsyncInstrument{T} <: AbstractInstrument{T} end

function examine_instrument(
    ins::AbstractInstrument;
    max_unit_length = 63,
    max_description_length = 1024,
)
    !isnothing(match(r"[a-zA-Z][_0-9a-zA-Z\\.\\-]{0,62}$", ins.name)) ||
        throw(ArgumentError("invalid name: $(ins.name)"))
    length(ins.unit) <= max_unit_length ||
        throw(ArgumentError("length of unit should be no more than $max_unit_length"))
    length(ins.description) <= max_description_length || throw(
        ArgumentError(
            "length of description should be no more than $max_description_length",
        ),
    )
end

(ins::AbstractSyncInstrument)(amount; kw...) =
    ins(Measurement(amount, StaticAttrs(values(kw))))
(ins::AbstractSyncInstrument)(m::Measurement) = push!(ins.meter.provider, ins => m)

(ins::AbstractAsyncInstrument)() = push!(ins.meter.provider, ins => ins.callback())

#####

struct Counter{T} <: AbstractSyncInstrument{T}
    meter::Meter
    name::String
    unit::String
    description::String
    function Counter{T}(name, meter; unit = "", description = "") where {T}
        c = new{T}(meter, name, unit, description)
        examine_instrument(c)
        push!(meter.provider, c)
        c
    end
end

(c::Counter{T})(; kw...) where {T} = c(one(T); kw...)

function (c::Counter)(m::Measurement)
    if m.value < 0
        throw(ArgumentError("amount must be non-negative, got $(m.value)"))
    end
    push!(c.meter.provider, c => m)
end

struct ObservableCounter{T,F} <: AbstractAsyncInstrument{T}
    callback::F
    meter::Meter
    name::String
    unit::String
    description::String
    function ObservableCounter{T}(
        callback::F,
        name,
        meter;
        unit = "",
        description = "",
    ) where {T,F}
        c = new{T,F}(callback, meter, name, unit, description)
        examine_instrument(c)
        push!(meter.provider, c)
        c
    end
end

struct Histogram{T} <: AbstractSyncInstrument{T}
    meter::Meter
    name::String
    unit::String
    description::String
    function Histogram{T}(name, meter; unit = "", description = "") where {T}
        h = new{T}(meter, name, unit, description)
        examine_instrument(h)
        push!(meter.provider, h)
        h
    end
end

struct ObservableGauge{T,F} <: AbstractAsyncInstrument{T}
    callback::F
    meter::Meter
    name::String
    unit::String
    description::String
    function ObservableGauge{T}(
        callback::F,
        name,
        meter;
        unit = "",
        description = "",
    ) where {T,F}
        g = new{T,F}(callback, meter, name, unit, description)
        examine_instrument(g)
        push!(meter.provider, g)
        g
    end
end

struct UpDownCounter{T} <: AbstractSyncInstrument{T}
    meter::Meter
    name::String
    unit::String
    description::String
    function UpDownCounter{T}(name, meter; unit = "", description = "") where {T}
        c = new{T}(meter, name, unit, description)
        examine_instrument(c)
        push!(meter.provider, c)
        c
    end
end

struct ObservableUpDownCounter{T,F} <: AbstractAsyncInstrument{T}
    callback::F
    meter::Meter
    name::String
    unit::String
    description::String
    function ObservableUpDownCounter{T}(
        callback::F,
        name,
        meter;
        unit = "",
        description = "",
    ) where {T,F}
        c = new{T,F}(callback, meter, name, unit, description)
        examine_instrument(c)
        push!(meter.provider, c)
        c
    end
end


text:
# ---
# title: 1049. Last Stone Weight II
# id: problem1049
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/last-stone-weight-ii/description/>
# hidden: true
# ---
#
# We have a collection of rocks, each rock has a positive integer weight.
#
# Each turn, we choose **any two rocks**  and smash them together.  Suppose the
# stones have weights `x` and `y` with `x <= y`.  The result of this smash is:
#
#   * If `x == y`, both stones are totally destroyed;
#   * If `x != y`, the stone of weight `x` is totally destroyed, and the stone of weight `y` has new weight `y-x`.
#
# At the end, there is at most 1 stone left.  Return the **smallest possible**
# weight of this stone (the weight is 0 if there are no stones left.)
#
#
#
# **Example 1:**
#
#
#
#     Input: [2,7,4,1,8,1]
#     Output: 1
#     Explanation:
#     We can combine 2 and 4 to get 2 so the array converts to [2,7,1,8,1] then,
#     we can combine 7 and 8 to get 1 so the array converts to [2,1,1,1] then,
#     we can combine 2 and 1 to get 1 so the array converts to [1,1,1] then,
#     we can combine 1 and 1 to get 0 so the array converts to [1] then that's the optimal value.
#
#
#
#
# **Note:**
#
#   1. `1 <= stones.length <= 30`
#   2. `1 <= stones[i] <= 100`
#
#
## @lc code=start
using LeetCode

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


text:
# TODO enable scattering points across pr and roc curves based on thresholds, tpr, fpr, or precision

abstract type AbstractCurve end

# by default always compute auc and curve points from all possible points
function auc(::Type{C}, args...; kwargs...) where {C<:AbstractCurve}
    auc_trapezoidal(curve(C, args...; kwargs...))
end
function auc(::Type{C}, enc::TwoClassEncoding, targets::AbstractArray{<:AbstractVector},
               scores::AbstractArray{<:RealVector}, args...; kwargs...) where {C<:AbstractCurve}
    auc_trapezoidal.(curve(C, enc, targets, scores, args...; kwargs...))
end
function auc(::Type{C}, targets::AbstractArray{<:AbstractVector},
               scores::AbstractArray{<:RealVector}, args...; kwargs...) where {C<:AbstractCurve}
    auc_trapezoidal.(curve(C, targets, scores, args...; kwargs...))
end
function curve(::Type{C}, args...; npoints=Inf) where {C<:AbstractCurve}
    apply(C, args...; npoints=npoints)
end
function curve(::Type{C}, cs::CMVector, args...; kwargs...) where {C<:AbstractCurve}
    apply(C, cs, args...)
end

macro curve(name)
    name_lw = Symbol(lowercase(string(name)))
    name_auc = Symbol(lowercase(string("au_", name)))

    quote
        abstract type $(esc(name)) <: AbstractCurve end

        Base.@__doc__  function $(esc(name_lw))(args...; kwargs...)
            curve($(esc(name)), args...; kwargs...)
        end

        function $(esc(name_auc))(args...; kwargs...)
            auc($(esc(name)), args...; kwargs...)
        end
    end
end

@recipe function f(::Type{Val{:mlcurve}}, x, y, z; indexes=Int[], diagonal=false)
    # main curve
    @series begin
        seriestype := :path
        marker     := :none
        x          := x
        y          := y
        ()
    end

    # points on the main curve
    if !isempty(indexes)
        @series begin
            primary           := false
            seriestype        := :scatter
            markerstrokecolor := :auto
            label             := ""
            x                 := x[indexes]
            y                 := y[indexes]
            ()
        end
    end

    # diagonal
    if diagonal && get(plotattributes, :xscale, :identity) === :identity
        @series begin
            primary    := false
            seriestype := :path
            fill       := false
            line       := (:red, :dash, 0.5)
            marker     := :none
            label      := ""
            x          := [0, 1]
            y          := [0, 1]
            ()
        end
    end
end

@shorthands mlcurve

apply(::Type{C}, targets, scores; kwargs...) where {C<:AbstractCurve} =
    apply(C, current_encoding(), targets, scores; kwargs...)
apply(::Type{C}, targets, scores, thres; kwargs...) where {C<:AbstractCurve} =
    apply(C, current_encoding(), targets, scores, thres; kwargs...)

function apply(::Type{C}, enc::TwoClassEncoding,
               targets::AbstractArray{<:AbstractVector},
               scores::AbstractArray{<:RealVector}, args...; kwargs...)  where {C<:AbstractCurve}
    return [apply(C, t, s, args...; kwargs...) for (t,s) in zip(targets, scores)]
end

function apply(::Type{C}, enc::TwoClassEncoding,
               targets::AbstractVector, scores::RealVector;
               npoints::Real=300, xscale::Symbol=:identity,
               xlims=(0, 1), kwargs...) where {C<:AbstractCurve}
    # maximal resolution
    if npoints >= length(targets) + 1
        thres = thresholds(scores)
    else
        if xscale === :identity
            quantils = range(0, stop=1.0, length=npoints)
        else
            # if user didn't provide custom xlims, compute them so that we sample efficiently around
            # the area with lowest nonzero FPR/TPR
            if xlims == (0, 1)
                xlims = (lowest_metric_value(C, enc, targets), 1)
            end
            quantils = exp10.(range(log10(xlims[1]), log10(min(1, xlims[2])), length=npoints))
        end
        thres = sampling_function(C)(enc, targets, scores, quantils)
    end
    return apply(C, enc, targets, scores, thres)
end

function apply(::Type{C},
               enc::TwoClassEncoding,
               targets::AbstractVector,
               scores::RealVector, thres::RealVector; kwargs...)  where {C<:AbstractCurve}
    if !(0 < sum(ispositive.(enc, targets)) < length(targets))
        throw(ArgumentError("Only one class present in `targets` with encoding $enc."))
    end
    return apply(C, ConfusionMatrix(enc, targets, scores, thres))
end

@recipe function f(::Type{C}, args...) where {C<:AbstractCurve}
    points = apply(C, args...; plotattributes...)
    xl, yl = _lims(points, plotattributes)
    xlims --> xl
    ylims --> yl
    delete!(plotattributes, :npoints)
    label := auc_label(plotattributes, auc(C, args...), args...)
    delete!(plotattributes, :aucshow)
    return points
end
@recipe function f(::Type{C}, cs::AbstractArray{<:CMVector}, args...) where {C<:AbstractCurve}
    points = apply.(C, cs)
    xl, yl = _lims(points, plotattributes)
    xlims --> xl
    ylims --> yl
    label := auc_label(plotattributes, auc_trapezoidal(points), args...)
    delete!(plotattributes, :aucshow)
    return points
end

function auc_label(plotattributes, auc_score, args...)
    user_label = get(plotattributes, :label, "AUTO")
    if get(plotattributes, :aucshow, false)
        auc_label = string.("auc: ", round.(100 * auc_score', digits = 2), "%")
        if user_label != "AUTO"
            return string.(user_label, " (", auc_label, ")")
        else
            return auc_label
        end
    else
        user_label
    end
end

function _lims((x,y)::Tuple, plotattributes)
    xscale = get(plotattributes, :xscale, :identity)
    yscale = get(plotattributes, :yscale, :identity)
    xlims = xscale === :identity ? (0, 1.01) : (minimum(x[x .> 0]), 1.01)
    ylims = yscale === :identity ? (0, 1.01) : (minimum(y[y .> 0]), 1.01)
    xlims, ylims
end

function _lims(points::Vector{<:Tuple}, plotattributes)
    limss = [_lims(p, plotattributes) for p in points]
    xlimss, ylimss = zip(limss...) |> collect
    xlimss_low = [x[1] for x in xlimss]
    xlimss_high = [x[2] for x in xlimss]
    ylimss_low = [y[1] for y in ylimss]
    ylimss_high = [y[2] for y in ylimss]
    xlims = (minimum(xlimss_low), maximum(xlimss_high))
    ylims = (minimum(ylimss_low), maximum(ylimss_high))
    xlims, ylims
end

# ROC curve
"""
    $(SIGNATURES)

Returns false positive rates and true positive rates.
"""
@curve ROCCurve
apply(::Type{ROCCurve}, cms::CMVector) = (false_positive_rate(cms), true_positive_rate(cms))
sampling_function(::Type{ROCCurve}) = threshold_at_fpr
# TODO provide a better estimate of these
# smallest possible FPR
function lowest_metric_value(::Type{ROCCurve}, enc::TwoClassEncoding, targets::AbstractVector)
    1.0 / sum(isnegative.(enc, targets))
end

@userplot ROCPlot

@recipe function f(h::ROCPlot)
    seriestype := :mlcurve
    diagonal   --> true
    legend     --> :bottomright
    fillrange  --> 0
    fillalpha  --> 0.15
    title      --> "ROC curve"
    xguide     --> "false positive rate"
    yguide     --> "true positive rate"
    xgrid      --> true
    ygrid      --> true
    aucshow    --> true

    (ROCCurve, h.args...)
end

# Precision-Recall curve
"""
    $(SIGNATURES)

Returns recalls and precisions.
"""
@curve PRCurve
apply(::Type{PRCurve}, cms::CMVector) = (recall(cms), precision(cms))
sampling_function(::Type{PRCurve}) = threshold_at_tpr
# TODO provide a better estimate of these
# smallest possible TPR
function lowest_metric_value(::Type{PRCurve}, enc::TwoClassEncoding, targets::AbstractVector)
    1.0 / sum(ispositive.(enc, targets))
end

@userplot PRPlot

@recipe function f(h::PRPlot)
    seriestype := :mlcurve
    legend     --> :bottomleft
    fillrange  --> 0
    fillalpha  --> 0.15
    title      --> "Precision-Recall curve"
    xguide     --> "recall"
    yguide     --> "precision"
    xgrid      --> true
    ygrid      --> true
    aucshow    --> true

    (PRCurve, h.args...)
end


text:
module MyPkg
# using QuadGK
# Write your package code here.
include("simpson.jl")
export Simpson
export find_num_in_vec

"""
    function find_num_in_vec(A::Vector, position_index::Vector{Int64}, x)

Find the number of values which in each vector elements equal to `x` according to the corresponding index position specified by the element in the `position_index` vector in the given vetcer `A`.

# Examples
```julia-repl
julia> A =  [[0.09,0.09,0.1],[0.3,0.09,0.1],[0.09]]
3-element Vector{Vector{Float64}}:
 [0.09, 0.09, 0.1]
 [0.3, 0.09, 0.1]
 [0.09]

julia> position_index =  [1,2,3]
3-element Vector{Int64}:
 1
 2
 3

julia> find_num_in_vec(A::Vector, position_index::Vector{Int64}, 0.09)
3-element Vector{Int64}:
2
1
1
```
"""
function find_num_in_vec(A::Vector, position_index::Vector{Int64}, x)
    number_in_vec = Vector{Int64}(undef, length(position_index))
    for i in eachindex(position_index)
        count = 0
        for j in eachindex(A[position_index[i]])
            if A[position_index[i]][j] == x
                count += 1
            end
        end
        number_in_vec[i] = count
    end
    return number_in_vec
end

function g(x)
    2*x
end
end


text:
x = [ x^2 for x = 1:5 ]                # Square of the numbers 1 to 5

A = [ x+y for x = 1:3, y = 1:5 ]       # 3x5 matrix with entries A_ij = i + j

B = [ x^y for x = 1:3, y = 1:5 ]       # 3x5 matrix with entries A_ij = i^j

x = [ x^2 for x = 1:10 if x%3 == 0]    # Only for x that are multiples of 3

y = Float64[ rand(-5:5) for i = 1:2, j = 1:10 ]

gen = (i*(i-1) for i = 1:10)        # A generator expression

for x in gen                        # Loop over all the values in the generator expression
    print(x, " ")
end

reshape(1:15, 3, 5)

permutedims(reshape(1:15, 5, 3))

repeat(1:2, 3)

repeat(1:2, 3, 4)

repeat([1 2; 3 4], 1, 2)


text:
const Cmode_t = Cushort

"Generic structure used for passing keys and data in and out of the database."
struct MDBValue
    size::Csize_t   # size of the data item
    data::Ptr{Nothing} # address of the data item
end

MDBValue() = MDBValue(zero(Csize_t), C_NULL)
MDBValue(_::Nothing) = MDBValue()
MDBValue(val::String) = MDBValue(sizeof(val), pointer(val))
function MDBValue(val::T) where {T}
    isbitstype(T) && error("Can not wrap a $T in MDBValue. Use a $T array instead")
    val_size = sizeof(eltype(val))*length(val)
    return MDBValue(val_size, pointer(val))
end

convert(::Type{T}, mdb_val_ref::Ref{MDBValue}) where {T} = _convert(T, mdb_val_ref[])
function _convert(::Type{String}, mdb_val::MDBValue)
    unsafe_string(convert(Ptr{UInt8}, mdb_val.data), mdb_val.size)
end
function _convert(::Type{Vector{T}}, mdb_val::MDBValue) where {T}
    res = unsafe_wrap(Array, convert(Ptr{UInt8}, mdb_val.data), mdb_val.size)
    reinterpret(T,res)
end
function _convert(::Type{T}, mdb_val::MDBValue) where {T}
    unsafe_load(convert(Ptr{T}, mdb_val.data))
end

# Environment Flags
# -----------------
@enum(EnvironmentFlags,
      EMPTY      = 0x00000000, # no flags
      FIXEDMAP   = 0x00000001, # mmap at a fixed address
      NOSUBDIR   = 0x00004000, # no environment directory
      NOSYNC     = 0x00010000, # don't fsync after commit
      RDONLY     = 0x00020000, # read only
      NOMETASYNC = 0x00040000, # don't fsync metapage after commit
      WRITEMAP   = 0x00080000, # use writable mmap
      MAPASYNC   = 0x00100000, # use asynchronous msync when #MDB_WRITEMAP is used
      NOTLS      = 0x00200000, # tie reader locktable slots to #MDB_txn objects instead of to threads
      NOLOCK     = 0x00400000, # don't do any locking, caller must manage their own locks
      NORDAHEAD  = 0x00800000, # don't do readahead (no effect on Windows)
      NOMEMINIT  = 0x01000000  # don't initialize malloc'd memory before writing to datafile
)

# Database Flags
# --------------
@enum(DatabaseFlags,
      REVERSEKEY = 0x00000002, # use reverse string keys
      DUPSORT    = 0x00000004, # use sorted duplicates
      INTEGERKEY = 0x00000008, # numeric keys in native byte order. The keys must all be of the same size.
      DUPFIXED   = 0x00000010, # with MDB_DUPSORT, sorted dup items have fixed size
      INTEGERDUP = 0x00000020, # with MDB_DUPSORT, dups are numeric in native byte order
      REVERSEDUP = 0x00000040, # with #MDB_DUPSORT, use reverse string dups
      CREATE     = 0x00040000  # create DB if not already existing
)

# Write Flags
# -----------
@enum(WriteFlags,
      NOOVERWRITE =0x00000010, # For put: Don't write if the key already exists.
#= Only for #MDB_DUPSORT
 * For put: don't write if the key and data pair already exist.
 * For mdb_cursor_del: remove all duplicate data items.
=#
      NODUPDATA = 0x00000020,
      CURRENT   = 0x00000040, # For mdb_cursor_put: overwrite the current key/data pair
      RESERVE   = 0x00010000, # For put: Just reserve space for data, don't copy it. Return a pointer to the reserved space.
      APPEND    = 0x00020000, # Data is being appended, don't split full pages.
      APPENDDUP = 0x00040000, # Duplicate data is being appended, don't split full pages.
      MULTIPLE  = 0x00080000  # Store multiple data items in one call. Only for DUPFIXED.
)

# Cursor `get` operations
# -----------------------
@enum(CursorOps,
      FIRST,          # Position at first key/data item
      FIRST_DUP,      # Position at first data item of current key. Only for #DUPSORT
      GET_BOTH,       # Position at key/data pair. Only for #MDB_DUPSORT
      GET_BOTH_RANGE, # Position at key, nearest data. Only for #MDB_DUPSORT
      GET_CURRENT,    # Return key/data at current cursor position
      GET_MULTIPLE,   # Return key and up to a page of duplicate data items from current cursor position. Move cursor to prepare for #NEXT_MULTIPLE. Only for #DUPFIXED
      LAST,           # Position at last key/data item
      LAST_DUP,       # Position at last data item of current key. Only for #DUPSORT
      NEXT,           # Position at next data item
      NEXT_DUP,       # Position at next data item of current key. Only for #DUPSORT
      NEXT_MULTIPLE,  # Return key and up to a page of duplicate data items from next cursor position. Move cursor to prepare for #NEXT_MULTIPLE. Only for #DUPFIXED
      NEXT_NODUP,     # Position at first data item of next key
      PREV,           # Position at previous data item
      PREV_DUP,       # Position at previous data item of current key. Only for #MDB_DUPSORT
      PREV_NODUP,     # Position at last data item of previous key
      SET,            # Position at specified key
      SET_KEY,        # Position at specified key, return key + data
      SET_RANGE)      # Position at first key greater than or equal to specified key.

"""Return the LMDB library version and version information

Function returns tuple `(VersionNumber, String)` that contains a library version and a library version string.
"""
function version()
    major = Cint[0]
    minor = Cint[0]
    patch = Cint[0]
    ver_str = ccall( (:mdb_version, liblmdb), Cstring, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), major, minor, patch)
    return VersionNumber(major[1],minor[1],patch[1]), unsafe_string(ver_str)
end

"""Return a string describing a given error code

Function returns description of the error as a string. It accepts following arguments:
* `err::Int32`: An error code.
"""
function errormsg(err::Cint)
    errstr = ccall( (:mdb_strerror, liblmdb), Cstring, (Cint,), err)
    return unsafe_string(errstr)
end

"""LMDB exception type"""
struct LMDBError <: Exception
    code::Cint
    msg::AbstractString
    LMDBError(code::Cint) = new(code, errormsg(code))
end
show(io::IO, err::LMDBError) = print(io, "Code[$(err.code)]: $(err.msg)")

""" Check if binary flag is set in provided value"""
isflagset(value, flag) = (value & flag) == flag


text:
#! /usr/bin/env julia

module problem102
export countTrianglesContainingOrigin

immutable Point
  x :: Int64
  y :: Int64
end

immutable Line
  a :: Point
  b :: Point
end

immutable Triangle
  a :: Point
  b :: Point
  c :: Point
end

function sideOfLine(pt :: Point, ln :: Line)
  return sign((ln.b.x - ln.a.x)*(pt.y - ln.a.y) - (ln.b.y - ln.a.y)*(pt.x - ln.a.x))
end

function pointInTriangle(pt :: Point, tri :: Triangle)
  s1 = sideOfLine(pt, Line(tri.a, tri.b))
  s2 = sideOfLine(pt, Line(tri.b, tri.c))
  s3 = sideOfLine(pt, Line(tri.c, tri.a))
  return (s1 == s2) && (s2 == s3)
end

function triangleContainsOrigin(tri :: Triangle)
  return pointInTriangle(Point(0,0), tri)
end

function readTriangles(filename :: String)
  data = readcsv(filename, Int64)
  return [Triangle(Point(data[i, 1], data[i, 2]), Point(data[i, 3], data[i, 4]), Point(data[i, 5], data[i, 6])) for i = 1:size(data,1)]
end

function countTrianglesContainingOrigin()
  tris = readTriangles("triangles.txt")
  return count(triangleContainsOrigin, tris)
end

using Base.Test

@test triangleContainsOrigin(Triangle(Point(-340,495), Point(-153,-910), Point(835,-947)))
@test !triangleContainsOrigin(Triangle(Point(-175,41), Point(-421,-714), Point(574,-645)))

end

println(problem102.countTrianglesContainingOrigin())

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE

const seg_integration_points = Dict(
    3 => ((5.0/9.0, -sqrt(3.0/5.0)),
          (8.0/9.0, 0.0),
          (5.0/9.0, +sqrt(3.0/5.0))))

"""
    calculate_mortar_matrices(slave_element_id::Int,
                              elements::Dict{Int, Vector{Int}},
                              element_types::Dict{Int, Symbol},
                              coords::Dict{Int, Vector{Float64}},
                              normals::Dict{Int, Vector{Float64}},
                              segmentation:MortarSegmentation)

Calculate mortar matrices De and Me for slave element.

# Example

```jldoctest
elements = Dict(
    1 => [1, 2],
    2 => [3, 4])
element_types = Dict(
    1 => :Seg2,
    2 => :Seg2)
coords = Dict(
    1 => [1.0, 2.0],
    2 => [3.0, 2.0],
    3 => [2.0, 2.0],
    4 => [0.0, 2.0])
normals = Dict(
    1 => [0.0, -1.0],
    2 => [0.0, -1.0])
segmentation = Dict(1 => [(2, [-1.0, 0.0])])
De, Me = calculate_mortar_matrices(1, elements, element_types,
                                   coords, normals, segmentation)

# output

([0.583333 0.166667; 0.166667 0.0833333], Dict(2=>[0.541667 0.208333; 0.208333 0.0416667]))

```

"""
function calculate_mortar_matrices(slave_element_id, elements, element_types,
                                   coords, normals, segmentation)
    @assert element_types[slave_element_id] == :Seg2
    # Initialization + calculate jacobian
    De = zeros(2, 2)
    Me = Dict{Int, Matrix{Float64}}()
    scon = elements[slave_element_id]
    xs1 = coords[scon[1]]
    xs2 = coords[scon[2]]
    ns1 = normals[scon[1]]
    ns2 = normals[scon[2]]
    J = norm(xs2-xs1)/2.0
    # 1. calculate De
    for (mid, (xi1, xi2)) in segmentation[slave_element_id]
        s = abs(xi2-xi1)/2.0
        for (w, ip) in seg_integration_points[3]
            xi_s = (1-ip)/2*xi1 + (1+ip)/2*xi2
            N1 = [(1-xi_s)/2 (1+xi_s)/2]
            De += w*N1'*N1*s*J
        end
    end
    # 2. calculate Me
    for (mid, (xi1, xi2)) in segmentation[slave_element_id]
        @assert element_types[mid] == :Seg2
        Me[mid] = zeros(2, 2)
        mcon = elements[mid]
        xm1 = coords[mcon[1]]
        xm2 = coords[mcon[2]]
        s = abs(xi2-xi1)/2.0
        for (w, xi) in seg_integration_points[3]
            xi_s = (1-xi)/2*xi1 + (1+xi)/2*xi2
            N1 = [(1-xi_s)/2 (1+xi_s)/2]
            n_s = N1[1]*ns1 + N1[2]*ns2
            x_g = N1[1]*xs1 + N1[2]*xs2
            xi_m = project_from_slave_to_master(Val{:Seg2}, x_g, n_s, xm1, xm2)
            N2 = [(1-xi_m)/2 (1+xi_m)/2]
            Me[mid] += w*N1'*N2*s*J
        end
    end
    return De, Me
end


text:
#Written By Johnathan Bizzano
using ..HyperSphere
import ..Error
import ..Layer
import ..Optimizer
import ...AbstractObject
import ...Double
import ..Initializer
import .Processor

export Model

struct Model{StorageArray, InputType, OutputType} <: AbstractObject
    fun::Function
    program::Program
    constants::StorageArray
    constant_range
    error_function::Error.Func
    memsettings::MemorySettings

    Model(designer::ModelDesigner{S}, fun, program, constants, constant_range, InputType, OutputType) where S =
        new{S, InputType, OutputType}(
                  fun,
                  program,
                  constants,
                  constant_range,
                  designer.error_function,
                  designer.memsettings)

    function (f::Model{S, I, O})(args::I)::O where {S, I, O}
        return f.fun(args)
    end

    Functions.trainer(f::Model; optimizer::MemoryWrapper{Optimizer.Func} = Optimizer.blackboxoptimizer()) = ModelTrainer(f, optimizer)

    function Base.string(x::Model)
        layerz = join([string(layer) for layer in x.layers], ", ")
        const_count = length(x.constants)
        err_func = x.error_function
        return "Model:\n\tLayers: $layerz\n\tConstant Count: $const_count\n\tError Function: $err_func\n"
    end
end

struct ModelTrainer{StorageType, InputType, OutputType} <: AbstractTrainer{InputType, -1}
        net::Model{StorageType, InputType, OutputType, Device, DeviceArray}
        optimizer::Optimizer.Func{StorageType, OutputType}

        ModelTrainer(net::Model{S, I, O}, optimizer::MemoryWrapper{Optimizer.Func}) where {S, I, O} = new{S, I, O}(net, optimizer(net.memsettings))

        function Base.string(x::ModelTrainer)
                layerz = join([string(layer) for layer in x.net.layers], ", ")
                const_count = length(x.net.constants)
                err_func = x.net.error_function
                return "ModelTrainer:\n\tLayers: $layerz\n\tConstant Count: $const_count\n\tError Function: $err_func\n"
        end
end

function Functions.train!(f::ModelTrainer, data::AbstractDataSet; epochs::Int = 1, IsVerbose=true, testset::AbstractDataSet=data)
    t = time()
    err_func = Error.Func2(() -> f.net.error_function(data, f.net))
    test_func = Error.Func2(() -> f.net.error_function(testset, f.net))

    if IsVerbose
        println("========MODEL TRAINING STARTED=======")
        println("Initial: Data Error:", err_func(), ". Test Error:", test_func())
        println("Number of Parameters:$(length(f.net.constants))")
    end

    for i in 1:epochs
        res = f.optimizer(f.net.constants, err_func, first(f.net.constant_range), last(f.net.constant_range))
        copy!(f.net.constants, res[1])
        reset(data)
        if IsVerbose
            println("Epoch:$i. Took (", time() - t, ") seconds. Data Error:", res[2], ". Test Error:", test_func())
            t = time()
        end
    end

    if IsVerbose
        println("========MODEL TRAINING ENDED=======")
    end

    (err_func(), test_func())
end


text:
immutable TensorMap{S,P,T}<:AbstractTensorMap{S,P,T}
    tensor::AbstractTensor{S,P,T}
    codom::P
    dom::P
    function TensorMap(t::AbstractTensor{S,P,T},codom::P,dom::P=codom)
        t in codom ⊗ dual(dom).' || throw(SpaceError("tensor cannot represent map from $dom to $codom"))
        new(t,codom,dom)
    end
end
tensormap{S,P,T}(t::AbstractTensor{S,P,T},codom::TensorSpace,dom::TensorSpace=codom)=TensorMap{S,P,T}(t,codom,dom)

# properties
codomain(A::TensorMap)=A.codom
domain(A::TensorMap)=A.dom

Base.isreal(A::TensorMap)=isreal(A.tensor)

Base.transpose(A::TensorMap)=tensormap(A.tensor.',dual(A.dom),dual(A.codom))
Base.ctranspose(A::TensorMap)=tensormap(A.tensor',conj(dual(A.dom)),conj(dual(A.codom)))

# comparison
==(A::TensorMap,B::TensorMap)=A.map==B.map

# multiplication with vector
function Base.A_mul_B!{S,P}(y::AbstractTensor{S,P},A::TensorMap{S,P},x::AbstractTensor{S,P})
    x in domain(A) || throw(SpaceError("x in not in domain of map"))
    y in codomain(A) || throw(SpaceError("y in not in codomain of map"))
    N1=length(codomain(A))
    N2=length(domain(A))
    tensorcontract!(1,A.tensor,[-(1:N2),reverse(1:N1)],'N',x,1:N1,'N',0,y,-(1:N2))
    return y
end
function Base.Ac_mul_B!{S,P}(y::AbstractTensor{S,P},A::TensorMap{S,P},x::AbstractTensor{S,P})
    x in dual(conj(codomain(A))) || throw(SpaceError("x in not in domain of map"))
    y in dual(conj(domain(A))) || throw(SpaceError("y in not in domain of map"))
    N1=length(codomain(A))
    N2=length(domain(A))
    tensorcontract!(1,A.tensor,[(1:N2),-reverse(1:N1)],'C',x,1:N2,'N',0,y,-(1:N1))
    return y
end
function Base.At_mul_B!{S,P}(y::AbstractTensor{S,P},A::TensorMap{S,P},x::AbstractTensor{S,P})
    x in dual(codomain(A)) || throw(SpaceError("x in not in domain of map"))
    y in dual(domain(A)) || throw(SpaceError("y in not in domain of map"))
    N1=length(codomain(A))
    N2=length(domain(A))
    tensorcontract!(1,A.tensor,[(1:N2),-reverse(1:N1)],'N',x,1:N2,'N',0,y,-(1:N1))
    return y
end


text:
InputFolder = './Images/EyeVesselsSpots/';
OutputFolder = './Results/Images/EyeVesselsSpotsFilam2/';
Lbl = 1;
Fill = -1;

@iA = '*.tif';

@fxg_kConnLocThrFilam [iA] > [M1];
params.GRad = 1;
params.LocMeanRad = 9;
params.LocMeanThr1 = 1.7;
params.LocMeanThr2 = 1.35;
params.MinSeedArea = 5;
/endf

@fxkg_kSklClean [M1, iA] > [M2];
params.MinBrchLgth = 5;
params.SearchRad = 9;
params.MinMean = 30;
params.MinArea = 25;
/endf

/show iA > M2;
/keep M2 > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
# Test type relations

using Distributions
using ForwardDiff: Dual

@assert UnivariateDistribution <: Distribution
@assert MultivariateDistribution <: Distribution
@assert MatrixDistribution <: Distribution

@assert DiscreteDistribution <: Distribution
@assert ContinuousDistribution <: Distribution

@assert DiscreteUnivariateDistribution <: DiscreteDistribution
@assert DiscreteUnivariateDistribution <: UnivariateDistribution
@assert ContinuousUnivariateDistribution <: ContinuousDistribution
@assert ContinuousUnivariateDistribution <: UnivariateDistribution
@assert DiscreteMultivariateDistribution <: DiscreteDistribution
@assert DiscreteMultivariateDistribution <: MultivariateDistribution
@assert ContinuousMultivariateDistribution <: ContinuousDistribution
@assert ContinuousMultivariateDistribution <: MultivariateDistribution
@assert DiscreteMatrixDistribution <: DiscreteDistribution
@assert DiscreteMatrixDistribution <: MatrixDistribution
@assert ContinuousMatrixDistribution <: ContinuousDistribution
@assert ContinuousMatrixDistribution <: MatrixDistribution

@testset "Test Sample Type" begin
    for T in (Float64,Float32,Dual{Nothing,Float64,0})
        @testset "Type $T" begin
            for d in (MvNormal,MvLogNormal,MvNormalCanon,Dirichlet)
                dist = d(map(T,ones(2)))
                @test eltype(typeof(dist)) == T
                @test eltype(rand(dist)) == eltype(dist)
            end
            dist = Distributions.mvtdist(map(T,1.0),map(T,[1.0 0.0; 0.0 1.0]))
            @test eltype(typeof(dist)) == T
            @test eltype(rand(dist)) == eltype(dist)
        end
    end
end

@testset "equality" begin
    dist1 = Normal(1, 1)
    dist2 = Normal(1.0, 1.0)

    # Check h is used
    @test hash(dist1, UInt(1)) != hash(dist1, UInt(2))

    @test dist1 == deepcopy(dist1)
    @test hash(dist1) == hash(deepcopy(dist1))
    @test dist1 == dist2
    @test isequal(dist1, dist2)
    @test isapprox(dist1, dist2)
    @test hash(dist1) == hash(dist2)

    dist3 = Normal(1, 0.8)
    @test dist1 != dist3
    @test !isequal(dist1, dist3)
    @test !isapprox(dist1, dist3)
    @test hash(dist1) != hash(dist3)
    @test isapprox(dist1, dist3, atol=0.3)

    dist4 = LogNormal(1, 1)
    @test dist1 != dist4
    @test !isequal(dist1, dist4)
    @test !isapprox(dist1, dist4)
    @test hash(dist1) != hash(dist4)
end