text:
using .Gadfly

export plotICX, plotICXRatio, plotICXRatio2

function plotICX(z::Integer)
    colors = [ "red", "green", "blue", "yellow", "lightgreen", "purple", "brown", "orange", "lightblue"]
    name = [ "K", "L1", "L2", "L3", "M1", "M2", "M3", "M4", "M5" ]
    lyr(bs, i, color) =
        Gadfly.layer(loge->log10(max(0.1, ionizationcrosssection(bs.z, i, 10.0^loge) / 1.0e-24)), log10(bs.edge[i]), log10(1.0e9), Gadfly.Theme(default_color = color))
    bs = BoteSalvatElectron[z]
    idx = eachindex(bs.edge)
    layers = [ lyr(bs, i, colors[i]) for i in idx]
    Gadfly.plot(layers...,
      Gadfly.Guide.manual_color_key("Shell", name[idx], colors[idx]),
      Gadfly.Guide.xlabel("log(E) (eV)"), Guide.ylabel("log(σ) (barns)"), Guide.title("Z = $z"),
      Gadfly.Coord.cartesian(xmin = round(log10(minimum(bs.edge)) - 0.5), xmax = 9.0, ymin = 0.0))
end

function plotICXRatio(z::Integer)
    colors = [ "red", "green", "blue", "yellow", "lightgreen", "purple", "brown", "orange", "lightblue"]
    name = [ "K", "L1", "L2", "L3", "M1", "M2", "M3", "M4", "M5" ]
    lyr(i, ip, color) =
        Gadfly.layer(u->ionizationcrosssection(z, i, u*edgeenergy(z,i)) / ionizationcrosssection(z, ip, u*edgeenergy(z,ip)), 1.2, 10.0, Gadfly.Theme(default_color = color))
    layers = ( lyr(i, i>4 ? 5 : (i>1 ? 2 : 1), colors[i]) for i in eachindex(name) )
    Gadfly.plot(layers...,
      Gadfly.Guide.manual_color_key("Shell", name[eachindex(name)], colors[eachindex(name)]),
      Gadfly.Guide.xlabel("U (eV/eV)"), Guide.ylabel("σ/σ' (barns/barns)"), Guide.title("Z = $z"),
      Gadfly.Coord.cartesian(xmin = 1.0, xmax = 10.0, ymin = 0.0))
end

function plotICXRatio2(u::AbstractFloat)
    colors = [ "red", "green", "blue", "yellow", "lightgreen", "purple", "brown", "orange", "lightblue"]
    name = [ "K", "L1", "L2", "L3", "M1", "M2", "M3", "M4", "M5" ]
    function lyr(i, ip)
        zs = filter(z->( hasedge(z,i) && hasedge(z, ip)),1:92)
        xs = map(z->ionizationcrosssection(z, i, u*edgeenergy(z,i)) / ionizationcrosssection(z, ip, u*edgeenergy(z,ip)),zs)
        Gadfly.layer(x=zs, y=xs, Gadfly.Theme(default_color = colors[i]))
    end
    layers = ( lyr(i, i>4 ? 5 : (i>1 ? 2 : 1)) for i in eachindex(name) )
    Gadfly.plot(layers..., Guide.title("U = $u"),
      Gadfly.Guide.manual_color_key("Shell", name[eachindex(name)], colors[eachindex(name)]),
      Gadfly.Guide.xlabel("Atomic Number"), Guide.ylabel("σ/σ' (barns/barns)"),
      Gadfly.Coord.cartesian(xmin = 1.0, xmax = 92.0, ymin = 0.0))
end


text:
using Conda

# Install Jupyter kernel-spec file.
include("kspec.jl")
kernelpath = installkernel("Julia")

# make it easier to get more debugging output by setting JULIA_DEBUG=1
# when building.
IJULIA_DEBUG = lowercase(get(ENV, "IJULIA_DEBUG", "0"))
IJULIA_DEBUG = IJULIA_DEBUG in ("1", "true", "yes")

# remember the user's Jupyter preference, if any; empty == Conda
prefsfile = joinpath(first(DEPOT_PATH), "prefs", "IJulia")
mkpath(dirname(prefsfile))
jupyter = get(ENV, "JUPYTER", isfile(prefsfile) ? readchomp(prefsfile) : Sys.isunix() && !Sys.isapple() ? "jupyter" : "")
condajupyter = normpath(Conda.SCRIPTDIR, exe("jupyter"))
if isempty(jupyter) || dirname(jupyter) == abspath(Conda.SCRIPTDIR)
    jupyter = condajupyter # will be installed if needed
elseif isabspath(jupyter)
    if !Sys.isexecutable(jupyter)
        @warn("ignoring non-executable JUPYTER=$jupyter")
        jupyter = condajupyter
    end
elseif jupyter != basename(jupyter) # relative path
    @warn("ignoring relative path JUPYTER=$jupyter")
    jupyter = condajupyter
elseif Sys.which(jupyter) === nothing
    @warn("JUPYTER=$jupyter not found in PATH")
end

function write_if_changed(filename, contents)
    if !isfile(filename) || read(filename, String) != contents
        write(filename, contents)
    end
end

# Install the deps.jl file:
deps = """
    const IJULIA_DEBUG = $(IJULIA_DEBUG)
    const JUPYTER = $(repr(jupyter))
"""
write_if_changed("deps.jl", deps)
write_if_changed(prefsfile, jupyter)

ext: jl


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

function rhs!(du, u, t,
              mesh::StructuredMesh{2}, equations,
              initial_condition, boundary_conditions, source_terms,
              dg::DG, cache)
  # Reset du
  @trixi_timeit timer() "reset ∂u/∂t" du .= zero(eltype(du))

  # Calculate volume integral
  @trixi_timeit timer() "volume integral" calc_volume_integral!(
    du, u, mesh,
    have_nonconservative_terms(equations), equations,
    dg.volume_integral, dg, cache)

  # Calculate interface fluxes
  @trixi_timeit timer() "interface flux" calc_interface_flux!(
    cache, u, mesh,
    have_nonconservative_terms(equations), equations,
    dg.surface_integral, dg)

  # Calculate boundary fluxes
  @trixi_timeit timer() "boundary flux" calc_boundary_flux!(
    cache, u, t, boundary_conditions, mesh, equations, dg.surface_integral, dg)

  # Calculate surface integrals
  @trixi_timeit timer() "surface integral" calc_surface_integral!(
    du, u, mesh, equations, dg.surface_integral, dg, cache)

  # Apply Jacobian from mapping to reference element
  @trixi_timeit timer() "Jacobian" apply_jacobian!(
    du, mesh, equations, dg, cache)

  # Calculate source terms
  @trixi_timeit timer() "source terms" calc_sources!(
    du, u, t, source_terms, equations, dg, cache)

  return nothing
end

function calc_volume_integral!(du, u,
                               mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}},
                               nonconservative_terms::Val{false}, equations,
                               volume_integral::VolumeIntegralWeakForm,
                               dg::DGSEM, cache)
  @unpack derivative_dhat = dg.basis
  @unpack contravariant_vectors = cache.elements

  @threaded for element in eachelement(dg, cache)
    for j in eachnode(dg), i in eachnode(dg)
      u_node = get_node_vars(u, equations, dg, i, j, element)

      flux1 = flux(u_node, 1, equations)
      flux2 = flux(u_node, 2, equations)

      # Compute the contravariant flux by taking the scalar product of the
      # first contravariant vector Ja^1 and the flux vector
      Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
      contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2
      for ii in eachnode(dg)
        multiply_add_to_node_vars!(du, derivative_dhat[ii, i], contravariant_flux1, equations, dg, ii, j, element)
      end

      # Compute the contravariant flux by taking the scalar product of the
      # second contravariant vector Ja^2 and the flux vector
      Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
      contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2
      for jj in eachnode(dg)
        multiply_add_to_node_vars!(du, derivative_dhat[jj, j], contravariant_flux2, equations, dg, i, jj, element)
      end
    end
  end

  return nothing
end

function calc_interface_flux!(cache, u,
                              mesh::StructuredMesh{2},
                              nonconservative_terms, # can be Val{true}/Val{false}
                              equations, surface_integral, dg::DG)
  @unpack elements = cache

  @threaded for element in eachelement(dg, cache)
    # Interfaces in negative directions
    # Faster version of "for orientation in (1, 2)"

    # Interfaces in x-direction (`orientation` = 1)
    calc_interface_flux!(elements.surface_flux_values,
                         elements.left_neighbors[1, element],
                         element, 1, u, mesh,
                         nonconservative_terms, equations,
                         surface_integral, dg, cache)

    # Interfaces in y-direction (`orientation` = 2)
    calc_interface_flux!(elements.surface_flux_values,
                         elements.left_neighbors[2, element],
                         element, 2, u, mesh,
                         nonconservative_terms, equations,
                         surface_integral, dg, cache)
  end

  return nothing
end

@inline function calc_interface_flux!(surface_flux_values, left_element, right_element,
                                      orientation, u,
                                      mesh::StructuredMesh{2},
                                      nonconservative_terms::Val{false}, equations,
                                      surface_integral, dg::DG, cache)
  # This is slow for LSA, but for some reason faster for Euler (see #519)
  if left_element <= 0 # left_element = 0 at boundaries
    return nothing
  end

  @unpack surface_flux = surface_integral
  @unpack contravariant_vectors, inverse_jacobian = cache.elements

  right_direction = 2 * orientation
  left_direction = right_direction - 1

  for i in eachnode(dg)
    if orientation == 1
      u_ll = get_node_vars(u, equations, dg, nnodes(dg), i, left_element)
      u_rr = get_node_vars(u, equations, dg, 1,          i, right_element)

      # If the mapping is orientation-reversing, the contravariant vectors' orientation
      # is reversed as well. The normal vector must be oriented in the direction
      # from `left_element` to `right_element`, or the numerical flux will be computed
      # incorrectly (downwind direction).
      sign_jacobian = sign(inverse_jacobian[1, i, right_element])

      # First contravariant vector Ja^1 as SVector
      normal_direction = sign_jacobian * get_contravariant_vector(1, contravariant_vectors,
                                                                  1, i, right_element)
    else # orientation == 2
      u_ll = get_node_vars(u, equations, dg, i, nnodes(dg), left_element)
      u_rr = get_node_vars(u, equations, dg, i, 1,          right_element)

      # See above
      sign_jacobian = sign(inverse_jacobian[i, 1, right_element])

      # Second contravariant vector Ja^2 as SVector
      normal_direction = sign_jacobian * get_contravariant_vector(2, contravariant_vectors,
                                                                  i, 1, right_element)
    end

    # If the mapping is orientation-reversing, the normal vector will be reversed (see above).
    # However, the flux now has the wrong sign, since we need the physical flux in normal direction.
    flux = sign_jacobian * surface_flux(u_ll, u_rr, normal_direction, equations)

    for v in eachvariable(equations)
      surface_flux_values[v, i, right_direction, left_element] = flux[v]
      surface_flux_values[v, i, left_direction, right_element] = flux[v]
    end
  end

  return nothing
end

@inline function calc_interface_flux!(surface_flux_values, left_element, right_element,
                                      orientation, u,
                                      mesh::StructuredMesh{2},
                                      nonconservative_terms::Val{true}, equations,
                                      surface_integral, dg::DG, cache)
  # See comment on `calc_interface_flux!` with `nonconservative_terms::Val{false}`
  if left_element <= 0 # left_element = 0 at boundaries
    return nothing
  end

  @unpack surface_flux = surface_integral
  @unpack contravariant_vectors, inverse_jacobian = cache.elements

  right_direction = 2 * orientation
  left_direction = right_direction - 1

  for i in eachnode(dg)
    if orientation == 1
      u_ll = get_node_vars(u, equations, dg, nnodes(dg), i, left_element)
      u_rr = get_node_vars(u, equations, dg, 1,          i, right_element)

      # If the mapping is orientation-reversing, the contravariant vectors' orientation
      # is reversed as well. The normal vector must be oriented in the direction
      # from `left_element` to `right_element`, or the numerical flux will be computed
      # incorrectly (downwind direction).
      sign_jacobian = sign(inverse_jacobian[1, i, right_element])

      # First contravariant vector Ja^1 as SVector
      normal_direction = sign_jacobian * get_contravariant_vector(1, contravariant_vectors,
                                                                  1, i, right_element)
    else # orientation == 2
      u_ll = get_node_vars(u, equations, dg, i, nnodes(dg), left_element)
      u_rr = get_node_vars(u, equations, dg, i, 1,          right_element)

      # See above
      sign_jacobian = sign(inverse_jacobian[i, 1, right_element])

      # Second contravariant vector Ja^2 as SVector
      normal_direction = sign_jacobian * get_contravariant_vector(2, contravariant_vectors,
                                                                  i, 1, right_element)
    end

    # If the mapping is orientation-reversing, the normal vector will be reversed (see above).
    # However, the flux now has the wrong sign, since we need the physical flux in normal direction.
    flux = sign_jacobian * surface_flux(u_ll, u_rr, normal_direction, equations)

    # Call pointwise nonconservative term; Done twice because left/right orientation matters
    # See Bohm et al. 2018 for details on the nonconservative diamond "flux"
    # Scale with sign_jacobian to ensure that the normal_direction matches that from the flux above
    noncons_primary   = sign_jacobian * noncons_interface_flux(u_ll, u_rr, normal_direction, :weak, equations)
    noncons_secondary = sign_jacobian * noncons_interface_flux(u_rr, u_ll, normal_direction, :weak, equations)

    for v in eachvariable(equations)
      surface_flux_values[v, i, right_direction, left_element] = flux[v] + noncons_primary[v]
      surface_flux_values[v, i, left_direction, right_element] = flux[v] + noncons_secondary[v]
    end
  end

  return nothing
end

# TODO: Taal dimension agnostic
function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic,
                             mesh::StructuredMesh{2}, equations, surface_integral, dg::DG)
  @assert isperiodic(mesh)
end

function calc_boundary_flux!(cache, u, t, boundary_condition,
                             mesh::StructuredMesh{2}, equations, surface_integral, dg::DG)
  calc_boundary_flux!(cache, u, t,
                      (boundary_condition, boundary_condition,
                       boundary_condition, boundary_condition),
                      mesh, equations, surface_integral, dg)
end

function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple,Tuple},
                             mesh::StructuredMesh{2}, equations, surface_integral, dg::DG)
  @unpack surface_flux_values = cache.elements
  linear_indices = LinearIndices(size(mesh))

  for cell_y in axes(mesh, 2)
    # Negative x-direction
    direction = 1
    element = linear_indices[begin, cell_y]

    for j in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1,
                                       boundary_conditions[direction],
                                       mesh, equations, surface_integral, dg, cache,
                                       direction, (1, j), (j,), element)
    end

    # Positive x-direction
    direction = 2
    element = linear_indices[end, cell_y]

    for j in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1,
                                       boundary_conditions[direction],
                                       mesh, equations, surface_integral, dg, cache,
                                       direction, (nnodes(dg), j), (j,), element)
    end
  end

  for cell_x in axes(mesh, 1)
    # Negative y-direction
    direction = 3
    element = linear_indices[cell_x, begin]

    for i in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2,
                                       boundary_conditions[direction],
                                       mesh, equations, surface_integral, dg, cache,
                                       direction, (i, 1), (i,), element)
    end

    # Positive y-direction
    direction = 4
    element = linear_indices[cell_x, end]

    for i in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2,
                                       boundary_conditions[direction],
                                       mesh, equations, surface_integral, dg, cache,
                                       direction, (i, nnodes(dg)), (i,), element)
    end
  end
end

function apply_jacobian!(du,
                         mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}},
                         equations, dg::DG, cache)
  @unpack inverse_jacobian = cache.elements

  @threaded for element in eachelement(dg, cache)
    for j in eachnode(dg), i in eachnode(dg)
      factor = -inverse_jacobian[i, j, element]

      for v in eachvariable(equations)
        du[v, i, j, element] *= factor
      end
    end
  end

  return nothing
end

end # @muladd


text:
# Methods to more easily handle financial data

# Check for various key financial field names
has_open(x::TS)::Bool = any(occursin.(r"(op)"i, String.(x.fields)))
has_high(x::TS)::Bool = any(occursin.(r"(hi)"i, String.(x.fields)))
has_low(x::TS)::Bool = any(occursin.(r"(lo)"i, String.(x.fields)))
has_volume(x::TS)::Bool = any(occursin.(r"(vo)"i, String.(x.fields)))
function has_close(x::TS; allow_settle::Bool=true, allow_last::Bool=true)::Bool
    columns = String.(x.fields)
    if allow_settle && allow_last
        return any([occursin(r"(cl)|(last)|(settle)"i, column) for column in columns])
    end
    if allow_last && !allow_settle
		return any([occursin(r"(cl)|(last)"i, column) for column in columns])
    end
    if allow_settle && !allow_last
		return any([occursin(r"(cl)|(settle)"i, column) for column in columns])
    end
    if !allow_last && !allow_settle
		return any([occursin(r"(cl)"i, column) for column in columns])
    end
    return false
end

# Identify OHLC(V) formats
is_ohlc(x::TS)::Bool = has_open(x) && has_high(x) && has_low(x) && has_close(x)
is_ohlcv(x::TS)::Bool = is_ohlc(x) && has_volume(x)

# Extractor functions
op(x::TS)::TS = x[:,findfirst([occursin(r"(op)"i, String(field)) for field in x.fields])]
hi(x::TS)::TS = x[:,findfirst([occursin(r"(hi)"i, String(field)) for field in x.fields])]
lo(x::TS)::TS = x[:,findfirst([occursin(r"(lo)"i, String(field)) for field in x.fields])]
vo(x::TS)::TS = x[:,findfirst([occursin(r"(vo)"i, String(field)) for field in x.fields])]
function cl(x::TS; use_adj::Bool=true, allow_settle::Bool=true, allow_last::Bool=true)::TS
    columns = String.(x.fields)
    if use_adj
        j = findfirst([occursin(r"(adj((usted)|\s|)+)?(cl)"i, column) for column in columns])
        if !isa(j, Nothing)
            return x[:,j]
        end
    else
        j = findfirst([occursin(r"(?!adj)*(cl(ose|))"i, column) for column in columns])
        if !isa(j, Nothing)
            return x[:,j]
        end
    end
    if allow_settle
        j = findfirst([occursin(r"(settle)"i, column) for column in columns])
        if !isa(j, Nothing)
            return x[:,j]
        end
    end
    if allow_last
        j = findfirst([occursin(r"(last)"i, column) for column in columns])
        if !isa(j, Nothing)
            return x[:,j]
        end
    end
    error("No closing prices found.")
end
ohlc(x::TS)::TS = [op(x) hi(x) lo(x) cl(x)]
ohlcv(x::TS)::TS = [op(x) hi(x) lo(x) cl(x) vo(x)]
hlc(x::TS)::TS = [hi(x) lo(x) cl(x)]
hl(x::TS)::TS = [hi(x) lo(x)]
hl2(x::TS)::TS = (hi(x) + lo(x)) * 0.5
hlc3(x::TS; args...)::TS = (hi(x) + lo(x) + cl(x; args...)) * 0.3333333333333333
ohlc4(x::TS; args...)::TS = (op(x) + hi(x) + lo(x) + cl(x; args...)) * 0.25


text:
using Test, Convex, Cliffords, SCS, SchattenNorms, Distributions, QuantumInfo
import Random

import Base.kron

set_default_solver(SCSSolver(verbose=0, eps=1e-6, max_iters=5_000))

function paulichannel(p::Vector{T}) where T <: Real
    n_ = log(4,length(p))
    if !isinteger(n_)
        error("Probability vector must have length 4^n for some integer n")
    end
    n = round(Int,n_)
    allpaulisops = map(m->liou(complex(m)),allpaulis(n))
    return reduce(+,map(*,p,allpaulisops))
end

function randp(n)
    return rand(Dirichlet(ones(4^n)))
end

function dnormp(p1,p2)
    return norm(p1-p2,1)
end

Random.seed!(123456)

p1 = randp(2)
p2 = randp(2)

pc1 = paulichannel(p1)
pc2 = paulichannel(p2)

dnormcptp(pc1,pc2)
SchattenNorms.dnormcptp2(pc1,pc2)
dnorm(pc1-pc2)

times = Vector[]
results = Vector[]
for i in 1:100
    println("Iteration ",i)

    p1 = randp(2)
    p2 = randp(2)

    pc1 = paulichannel(p1)
    pc2 = paulichannel(p2)

    tic()
    r1 = dnormcptp(pc1,pc2)
    te1 = toc()

    tic()
    r2 = SchattenNorms.dnormcptp2(pc1,pc2)
    te2 = toc()

    tic()
    r3 = dnorm(pc1-pc2)
    te3 = toc()

    push!(times, [te1, te2, te3])
    push!(results, [r1, r2, r3, dnormp(p1,p2)])
end
times = reduce(hcat,times)'
results = reduce(hcat,results)'
results = results .- results[:,4]


text:
using KernelRidgeRegression
using StatsBase
using MLKernels

N = 5000
x = rand(1, N) * 4π - 2π
yy = sinc.(x) # vec(sinc.(4 .* x) .+ 0.2 .* sin.(30 .* x))
y = squeeze(yy + 0.1randn(1, N), 1)
xnew = collect(-2.5π:0.01:2.5π)'

@time mykrr = fit(KRR, x, y, 1e-3/5000, GaussianKernel(1.0))
showcompact(mykrr)
show(mykrr)

@time mynystkrr = fit(KernelRidgeRegression.NystromKRR, x, y, 1e-3/5000, 280, MLKernels.GaussianKernel(100.0))
showcompact(mynystkrr)
show(mynystkrr)

@time myfastkrr = fit(KernelRidgeRegression.FastKRR, x, y, 4/5000, 11, MLKernels.GaussianKernel(100.0))
showcompact(myfastkrr)
show(myfastkrr)

@time mytnkrr = fit(KernelRidgeRegression.TruncatedNewtonKRR, x, y, 4/5000, MLKernels.GaussianKernel(100.0), 0.5, 200)
showcompact(mytnkrr)
show(mytnkrr)

@time myrandkrr = fit(KernelRidgeRegression.RandomFourierFeatures, x, y, 1e-3/5000, 200 , 1.0)
showcompact(myrandkrr)
show(myrandkrr)


text:
using CSV

const defdir = joinpath(dirname(@__FILE__), "..", "datasets")

function get1cdtdata(dir)
       mkpath(joinpath(defdir, "synthetic"))
       path = download("https://raw.githubusercontent.com/Conradox/datastreams/master/sinthetic/1CDT.csv")
    mv(path, joinpath(defdir, "synthetic/1CDT.csv"))
end

function getug2c5ddata(dir)
    mkpath(joinpath(defdir, "synthetic"))
    path = download("https://raw.githubusercontent.com/Conradox/datastreams/master/sinthetic/UG_2C_5D.csv")
    mv(path, joinpath(defdir, "synthetic/UG_2C_5D.csv"))
end

function Dataset1CDT(batch::Int)::BatchStream
    filename = "$(defdir)/synthetic/1CDT.csv"

    isfile(filename) || get1cdtdata(defdir)

    data = CSV.read(filename; header = false)

    conn = EasyStream.TablesConnector(data)

    stream = BatchStream(conn; batch = batch)

    return stream
end

Dataset1CDT() = Dataset1CDT(1)

function DatasetUG_2C_5D(batch::Int)::BatchStream
    filename = "$(defdir)/synthetic/UG_2C_5D.csv"

    isfile(filename) || getug2c5ddata(defdir)

    data = CSV.read(filename; header = false)

    conn = EasyStream.TablesConnector(data)

    stream = BatchStream(conn, batch)

    return stream
end

DatasetUG_2C_5D() = DatasetUG_2C_5D(1)


text:
# https://www.codewars.com/kata/50654ddff44f800200000004

module Solution
  export multiply
  function multiply(a, b)
    a * b
  end
end

max_stars_repo_name: hiljusti/codewars-solutions

max_issues_repo_name: hiljusti/codewars-solutions

max_forks_repo_name: hiljusti/codewars-solutions


text:
# Configuration for nearly neutral simuation
# Run with small values of N, ngens for testing
# Example run from command line:  julia run.jl examples/ia_example0
function dfemod(x)
  dfe_mod(x,fit_inc=1.2)
end
const nn_simtype = 1    # nearly neutral infinite alleles
const type_str = "MD"
@everywhere const N_list     = [8]     # sample size n list, popsize N may be larger
const N_mu_list= [2.0]   # Mutation rate as a multiple of 1.0/N
const L = 1           # number of loci
const ngens = 6      # Generations after burn-in
const burn_in= 0.0    # generations of burn_in as a multiple of 1/mu
const use_poplist=false   # save a list of populations as a check of correctness
dfe=dfemod
dfe_str = "dfe_mod fit_inc: 1.2"

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
using BAT
using Test
using Statistics
using StatsBase
using Distributions
using HypothesisTests

@testset "Funnel Distribution" begin
    #Tests different forms of instantiate Funnel Distribution
    @test @inferred(BAT.FunnelDistribution(1., 2., 3)) isa BAT.FunnelDistribution
    @test @inferred(BAT.FunnelDistribution()) isa BAT.FunnelDistribution
    @test @inferred(BAT.FunnelDistribution(a = 1.0, b = 0.5, n = 3)) isa BAT.FunnelDistribution

    funnel = BAT.FunnelDistribution(a = Float32(1.0), b = Float32(0.5), n = Int32(3))

    #Check Parameters
    @test @inferred(StatsBase.params(funnel)) == (Float32(1.0), Float32(0.5), Int32(3))

    @test @inferred(Base.length(funnel)) == 3 #Number of dimensions
    @test @inferred(Base.eltype(funnel)) == Float32 #Element Type

    #Check mean and covariance
    @test @inferred(Statistics.mean(funnel)) == [0.0, 0.0, 0.0] #Check mean
    @test isapprox(@inferred(Statistics.cov(funnel)), [1.00256 -0.0124585 -0.00373376;
    -0.0124585 7.04822 -0.165097;  -0.00373376  -0.165097  7.1126], rtol = 1e-1, atol = 0)

    #logpdf
    @test isapprox(@inferred(Distributions._logpdf(funnel, [0., 0., 0.])), -2.75681, atol = 1e-5)

    #KS Test
    #Test the constant-variance Gaussian
    funnel = BAT.FunnelDistribution(a = 1., b = 0., n = 1)
    ks_test = HypothesisTests.ExactOneSampleKSTest(rand(funnel, 10^6)[:], Normal(0., 1.))
    @test pvalue(ks_test) > 0.05

    #Test the variable-variance Gaussian
    funnel = BAT.FunnelDistribution(a = 0., b = 0., n = 2)
    ks_test = HypothesisTests.ExactOneSampleKSTest(rand(funnel, 10^6)[2,:], Normal(0., 1.))
    @test pvalue(ks_test) > 0.05

end


text:
"""
LDA Slater exchange (DOI: 10.1017/S0305004100016108 and 10.1007/BF01340281)
"""
energy_per_particle(::Val{:lda_x}, ρ) = -3/4 * cbrt(3/π * ρ)


text:
"""
    ConjugateModel <: ProbabilisticModel

`ConjugateModel` is a `ProbabilisticModel` with a conjugate prior of the
corresponding likelihood function.
"""
abstract type ConjugateModel <: ProbabilisticModel end

function Base.show(io::IO, model::ConjugateModel)
    modelinfo = replace(string(model.dist), "Distributions."=>"", count=1)
    message = "$(typeof(model)) : $(modelinfo)"
    print(io, message)
end

"""
    ConjugateBernoulli <: ConjugateModel

Bernoulli likelihood with Beta distribution as the conjugate prior.

```julia
ConjugateBernoulli(α, β)              # construct a ConjugateBernoulli

update!(model, stats)             # update model with statistics from data
samplepost(model, numsamples)    # sampling from the posterior distribution
samplestats(model, numsamples)   # sampling statistics from the data generating distribution
```

"""
mutable struct ConjugateBernoulli <: ConjugateModel
    dist::Beta
    function ConjugateBernoulli(α, β)
        return new(Beta(α, β))
    end
end

defaultparams(::ConjugateBernoulli) = [:θ]

function update!(model::ConjugateBernoulli, stats::BetaStatistics)
    numsuccesses = stats.s
    numtrials = stats.n
    α = model.dist.α + numsuccesses
    β = model.dist.β + numtrials - numsuccesses
    model.dist = Beta(α, β)
    return nothing
end

"""
    ConjugateExponential <: ConjugateModel

Exponential likelihood with Gamma distribution as the conjugate prior.

```julia
ConjugateExponential(α, β)            # construct a ConjugateExponential

update!(model, stats)             # update model with statistics from data
samplepost(model, numsamples)    # sampling from the posterior distribution
samplestats(model, numsamples)   # sampling statistics from the data generating distribution
```
"""
mutable struct ConjugateExponential <: ConjugateModel
    dist::Gamma
    function ConjugateExponential(α, θ)
        return new(Gamma(α, θ))
    end
end

defaultparams(::ConjugateExponential) = [:θ]

mutable struct ConjugatePoisson <: ConjugateModel
    dist::Gamma
    function ConjugatePoisson(α, θ)
        return new(Gamma(α, θ))
    end
end

defaultparams(::ConjugatePoisson) = [:λ]

function update!(model::T, stats::S) where
        {T<:Union{ConjugateExponential, ConjugatePoisson},
         S<:Union{GammaStatistics}}
    n = stats.n
    x̄ = stats.x̄
    α = model.dist.α + n
    θ = model.dist.θ / (1 + model.dist.θ * n * x̄)
    model.dist = Gamma(α, θ)
    return nothing
end

"""
    ConjugateNormal <: ConjugateModel

Normal likelihood and Normal Inverse Gamma distribution as the
conjugate prior.

## Parameters

- μ: mean of normal distribution
- v:  scale variance of Normal
- α:  shape of Gamma distribution
- θ:  scale of Gamma distribution

```julia
ConjugateNormal(μ, v, α, θ)           # construct a ConjugateNormal

update!(model, stats)             # update model with statistics from data
samplepost(model, numsamples)    # sampling from the posterior distribution
samplestats(model, numsamples)   # sampling statistics from the data generating distribution
```

## References

- The update rule for Normal distribution is based on this
  [lecture notes](https://people.eecs.berkeley.edu/~jordan/courses/260-spring10/other-readings/chapter9.pdf).
"""
mutable struct ConjugateNormal{S} <: ConjugateModel
    dist::S
    function ConjugateNormal{Normal}(μ::Real, σ::Real)
        T = promote_type(typeof(μ), typeof(σ))
        return new(Normal{T}(μ, σ))
    end
    function ConjugateNormal{NormalInverseGamma}(μ::Real, v::Real, α::Real, θ::Real)
        return new(NormalInverseGamma(μ, v, α, θ))
    end
end

ConjugateNormal(μ::T, v::T, α::T, θ::T) where T <: Real = ConjugateNormal{NormalInverseGamma}(μ, v, α, θ)
ConjugateNormal(μ::T, σ::T) where T <: Real = ConjugateNormal{Normal}(μ, σ)

defaultparams(::ConjugateNormal{NormalInverseGamma}) = [:μ]

"""
    ConjugateLogNormal(μ, v, α, θ)

A model with Normal likelihood and Normal Inverse distribution with log transformed data.
Notice `LogNormal` in `Distributions.jl` takes mean and standard deviation of \$\\log(x)\$
instead of \$x\$ as the input parameters.

```julia
ConjugateLogNormal(μ, v, α, θ) # construct a ConjugateLogNormal

lognormalparams(μ_logx, σ²_logx) # convert normal parameters to log-normal parameters
update!(model, stats)              # update model with statistics from data
samplepost(model, numsamples)     # sampling from the posterior distribution
samplestats(model, numsamples)    # sampling statistics from the data generating distributio
```

"""
mutable struct ConjugateLogNormal{S} <: ConjugateModel
    dist::S
    function ConjugateLogNormal{NormalInverseGamma}(μ::Real, v::Real, α::Real, θ::Real)
        return new(NormalInverseGamma(μ, v, α, θ))
    end
end

ConjugateLogNormal(μ::T, v::T, α::T, θ::T) where T <: Real = ConjugateLogNormal{NormalInverseGamma}(μ, v, α, θ)

function lognormalparams(μ_logx, σ²_logx)
    μ_x  = @. exp(μ_logx + σ²_logx / 2)
    σ²_x = @. (exp(σ²_logx) - 1) * exp(2 * μ_logx + σ²_logx)
    return (μ_x, σ²_x)
end

defaultparams(::ConjugateLogNormal{NormalInverseGamma}) = [:μ_x]

function update!(model::Union{ConjugateNormal{NormalInverseGamma},ConjugateLogNormal{NormalInverseGamma}},
        stats::Union{NormalStatistics,LogNormalStatistics})
    n, x̄, sdx = getall(stats)
    ΣΔx² = sdx^2 * (n - 1)

    μ0 = model.dist.μ
    v0 = model.dist.v
    α0 = model.dist.α
    θ0 = model.dist.θ
    inv_v0 = 1.0 / v0
    inv_v = inv_v0 + n

    μ = (inv_v0 * μ0 + n * x̄) / inv_v
    v = 1 / inv_v
    α = α0 + n / 2
    θ = θ0 + 0.5 * (ΣΔx² + (n * inv_v0) * (x̄ - μ0)^2 * v)
    model.dist = NormalInverseGamma(μ, v, α, θ)
    return nothing
end

abstract type ChainOperator end

struct MultiplyOperator <: ChainOperator end
(m::MultiplyOperator)(a, b) = .*(a, b)

"""
    ChainedModel <: ProbabilisticModel

`ChainedModel` is a combination of `ConjugateModel`s chained by the specified operator.
It can be used to model a multiple step process.
"""
struct ChainedModel <: ProbabilisticModel
    models::Vector{ConjugateModel}
    operators::Vector{ChainOperator}
    function ChainedModel(models::Vector{ConjugateModel}, operators)
        length(models) - 1 == length(operators) || error("need to specify (number of model - 1) chaining operators")
        return new(models, operators)
    end
end

function ChainedModel(models::Vector{ConjugateModel})
    operators = [MultiplyOperator() for _ in 1:(length(models) -1)]
    return ChainedModel(models, operators)
end

function update!(chainedmodel::ChainedModel, listofstats::Vector{T}) where T <: ModelStatistics
    @assert(length(chainedmodel.models) == length(listofstats),
            "Number of model should be equal to number of statistics.")
    for (model, stats) = zip(chainedmodel.models, listofstats)
        update!(model, stats)
    end
    return nothing
end

function defaultparams(chainedmodel::ChainedModel)
    params = Symbol[]
    for model in chainedmodel.models
        push!(params, defaultparams(model)[1])
    end
    return params
end

"""
    samplepost(model, numsamples)

Sample from the posterior distribution of the model.
"""
function samplepost(model::ConjugateBernoulli, numsamples::Int)
    return BernoulliPosteriorSample(rand(model.dist, numsamples))
end

function samplepost(model::ConjugateExponential, numsamples::Int)
    # Gamma distribution generates the rate (λ) of Exponential distribution
    # convert it to scale θ to make it consistent with Distributions.jl
    λs = rand(model.dist, numsamples)
    θs = 1 ./ λs
    return ExponentialPosteriorSample(θs)
end

function samplepost(model::ConjugatePoisson, numsamples::Int)
    inv_rates = rand(model.dist, numsamples)
    λs = 1 ./ inv_rates
    return PoissonPosteriorSample(λs)
end

function samplepost(model::ConjugateNormal{NormalInverseGamma}, numsamples::Int)
    μ, σ² = rand(model.dist, numsamples)
    return NormalPosteriorSample(μ, σ²)
end

function samplepost(model::ConjugateLogNormal{NormalInverseGamma}, numsamples::Int)
    # sample for normal means and variances
    μ_logx, σ²_logx = rand(model.dist, numsamples)
    μ_x, σ²_x = lognormalparams(μ_logx, σ²_logx)
    return LogNormalPosteriorSample(μ_logx, σ²_logx, μ_x, σ²_x)
end

function samplepost(model::T, parameter::Symbol, numsamples::Int) where T <: ConjugateModel
    samples = samplepost(model, numsamples)
    return getfield(samples, parameter)
end

function samplepost(model::T, parameters::Vector{Symbol}, numsamples::Int)  where T <: ConjugateModel
    parameter = toparameter(parameters)
    return samplepost(model, parameter, numsamples)
end

function samplepost(models::Vector{T}, parameters::Vector{Symbol}, numsamples::Int) where T <: ProbabilisticModel
    samples = Array{Float64,2}(undef, numsamples, length(models))
    for (modelindex, model) in enumerate(models)
        sample = samplepost(model, parameters, numsamples)
        samples[:, modelindex] = sample
    end
    return samples
end

function samplepost(chainedmodel::ChainedModel, parameters::Vector{Symbol}, numsamples::Int)
    length(parameters) == length(chainedmodel.models) ||
        throw(ArgumentError("Number of parameters must be equal to number of chained models"))
    models = chainedmodel.models
    operators = chainedmodel.operators
    post_samples = [samplepost(model, parameter, numsamples) for (model, parameter) in zip(models, parameters)]
    chainedsample = post_samples[1]
    for i = 2:length(post_samples)
        chainedsample = operators[i - 1](chainedsample, post_samples[i])
    end
    return chainedsample
end

"""
    samplestats(model, dist, numsamples)

Sample from the distribution of the data generating process, and
calculate the corresponding statistics for the model.
"""
function samplestats(::ConjugateBernoulli, dist::Bernoulli, numsamples::Integer)
    return BetaStatistics(rand(dist, numsamples))
end

function samplestats(::ConjugateExponential, dist::Exponential, numsamples::Integer)
    return GammaStatistics(rand(dist, numsamples))
end

function samplestats(::ConjugatePoisson, dist::Poisson, numsamples::Integer)
    return GammaStatistics(rand(dist, numsamples))
end

function samplestats(::ConjugateNormal{NormalInverseGamma}, dist::Normal, numsamples::Integer)
    return NormalStatistics(rand(dist, numsamples))
end

function samplestats(::ConjugateLogNormal{NormalInverseGamma}, dist::LogNormal, numsamples::Integer)
    return LogNormalStatistics(rand(dist, numsamples))
end

function samplestats(chainedmodel::ChainedModel, dists::Vector{T}, listofnumsamples::Vector{Int}) where T
    listofstats = Vector{ModelStatistics}()
    for (model, dist, numsamples) in zip(chainedmodel.models, dists, listofnumsamples)
        push!(listofstats, samplestats(model, dist, numsamples))
    end
    return listofstats
end


text:
module MatOp

using Blocks
using Compat

importall Blocks
importall Base

export MatOpBlock, Block, RandomMatrix, op, convert, *

##
# RamdomMatrix: A way to represent distributed memory random matrix that is easier to work with blocks
# Actually just a placeholder for eltype, dims and rng seeds, with a few AbstractMatrix methods defined on it
type RandomMatrix{T} <: AbstractMatrix{T}
    t::Type{T}
    dims::NTuple{2,Int}
    seed::Int
end
eltype{T}(m::RandomMatrix{T}) = T
size(m::RandomMatrix) = m.dims

##
# MatrixSplits hold remote refs of parts of the matrix, along with the matrix definition
type MatrixSplits
    m::AbstractMatrix
    splitrefs::Dict
end

##
# Block definition on DenseMatrix.
# Each part contains the part dimensions.
# Parts are fetched from master node on first use, and their resusable references are stored at master node.
# This is not a terribly useful type, just used as an example implementation.
type DenseMatrixPart
    r::RemoteRef
    split::Tuple
end

function Block(A::DenseMatrix, splits::Tuple)
    msplits = MatrixSplits(A, Dict())
    r = RemoteRef()
    put!(r, msplits)
    blks = [DenseMatrixPart(r, split) for split in splits]
    Block(A, blks, Blocks.no_affinity, as_it_is, as_it_is)
end

##
# Block definition on RandomMatrix.
# Each part contains the type, part dimensions, and rng seed.
# Parts are created locally on first use, and their resusable references are stored at master node
# This is not a terribly useful type, just used as an example implementation.
type RandomMatrixPart
    r::RemoteRef
    t::Type
    split::Tuple
    splitseed::Int
end

function Block(A::RandomMatrix, splits::Tuple)
    msplits = MatrixSplits(A, Dict())
    r = RemoteRef()
    put!(r, msplits)
    splitseeds = A.seed + (1:length(splits))
    blks = [RandomMatrixPart(r, eltype(A), splits[idx], splitseeds[idx]) for idx in 1:length(splits)]
    Block(A, blks, Blocks.no_affinity, as_it_is, as_it_is)
end

function matrixpart_get(blk)
    refpid = blk.r.where # get the pid of the master process
    # TODO: need a way to avoid unnecessary remotecalls after the initial fetch
    remotecall_fetch(()->matrixpart_get(blk, myid()), refpid)
end

function matrixpart_get(blk, pid)
    msplits = fetch(blk.r)
    splitrefs = msplits.splitrefs
    key = (pid, blk.split)
    get(splitrefs, key, nothing)
end

function matrixpart_set(blk, part_ref)
    refpid = blk.r.where # get the pid of the master process
    remotecall_wait(()->matrixpart_set(blk, myid(), part_ref), refpid)
end

function matrixpart_set(blk, pid, part_ref)
    msplits = fetch(blk.r)
    splitrefs = msplits.splitrefs
    key = (pid, blk.split)
    splitrefs[key] = part_ref
    nothing
end

function matrixpart(blk)
    # check at the master process if we already have this chunk, and get its ref
    chunk_ref = matrixpart_get(blk)
    (chunk_ref === nothing) || (return fetch(chunk_ref))

    # create the chunk
    part = matrixpart_create(blk)
    part_ref = RemoteRef()
    put!(part_ref, part)
    # update its reference
    matrixpart_set(blk, part_ref)
    # return the chunk
    part
end

function matrixpart_create(blk::RandomMatrixPart)
    part_size = map(length, blk.split)
    srand(blk.splitseed)
    rand(blk.t, part_size...)
end

function matrixpart_create(blk::DenseMatrixPart)
    splits = fetch(blk.r)
    A = splits.m
    part_range = blk.split
    A[part_range...]
end

convert{T}(::Type{DenseMatrix}, r::Block{RandomMatrix{T}}) = convert(DenseMatrix{T}, r)
function convert{T}(::Type{DenseMatrix{T}}, r::Block{RandomMatrix{T}})
    A = r.source
    ret = zeros(eltype(A), size(A))
    for blk in blocks(r)
        part_range = blk.split
        ret[part_range...] = matrixpart_create(blk)
    end
    ret
end

# Blocked operations on matrices
type MatOpBlock
    mb1::Block
    mb2::Block
    oper::Symbol

    function MatOpBlock(m1::AbstractMatrix, m2::AbstractMatrix, oper::Symbol, np::Int=0)
        s1 = size(m1)
        s2 = size(m2)

        (0 == np) && (np = nworkers())
        np = min(max(s1...), max(s2...), np)

        (splits1, splits2) = (oper == :*) ? matop_block_mul(s1, s2, np) : error("operation $oper not supported")

        mb1 = Block(m1, splits1)
        mb2 = Block(m2, splits2)
        new(mb1, mb2, oper)
    end
end

function Block(mb::MatOpBlock)
    (mb.oper == :*) && return block_mul(mb)
    error("operation $(mb.oper) not supported")
end

op{T<:MatOpBlock}(blk::Block{T}) = (eval(blk.source.oper))(blk)

##
# internal methods
function common_factor_around(around::Int, nums::Int...)
    gf = gcd(nums...)
    ((gf == 1) || (gf < around)) && return gf

    factors = Int[]

    n = @compat(Int(floor(gf/around)))+1
    while n < gf
        (0 == (gf%n)) && (push!(factors, @compat(round(Int,gf/n))); break)
        n += 1
    end

    n = @compat(Int(floor(gf/around)))
    while n > 0
        (0 == (gf%n)) && (push!(factors, @compat(round(Int,gf/n))); break)
        n -= 1
    end
    (length(factors) == 1) && (return factors[1])
    ((factors[2]-around) > (around-factors[1])) ? factors[1] : factors[2]
end

function mat_split_ranges(dims::Tuple, nrsplits::Int, ncsplits::Int)
    row_splits = Base.splitrange(dims[1], nrsplits)
    col_splits = Base.splitrange(dims[2], ncsplits)
    splits = Array(Tuple,0)
    for cidx in 1:ncsplits
        for ridx in 1:nrsplits
            push!(splits, (row_splits[ridx],col_splits[cidx]))
        end
    end
    splits
end

# Input:
#   - two matrices (or their dimensions)
#   - number of processors to split over
# Output:
#   - tuples of ranges to split each matrix into, suitable for block multiplication
function matop_block_mul(s1::Tuple, s2::Tuple, np::Int)
    fc = common_factor_around(round(Int,min(s1[2],s2[1])/np), s1[2], s2[1])
    f1 = common_factor_around(round(Int,s1[1]/np), s1[1])
    f2 = common_factor_around(round(Int,s2[2]/np), s2[2])

    splits1 = mat_split_ranges(s1, round(Int,s1[1]/f1), round(Int,s1[2]/fc)) # splits1 is f1 x fc blocks
    splits2 = mat_split_ranges(s2, round(Int,s2[1]/fc), round(Int,s2[2]/f2)) # splits2 is fc x f2 blocks
    (tuple(splits1...), tuple(splits2...))
end

function block_mul(mb::MatOpBlock)
    blklist = Any[]
    afflist = Any[]
    proclist = workers()
    # distribute by splits on m1
    for blk1 in blocks(mb.mb1)
        split1 = blk1.split
        proc = shift!(proclist)
        push!(proclist, proc)
        for blk2 in blocks(mb.mb2)
            split2 = blk2.split
            (split1[2] != split2[1]) && continue
            result_range = (split1[1], split2[2])
            push!(blklist, (result_range, blk1, blk2, proc))
            # set affinities of all blocks each split of m1 needs to same processor
            push!(afflist, proc)
        end
    end

    Block(mb, blklist, afflist, as_it_is, as_it_is)
end

##
# operations
function *{T<:MatOpBlock}(blk::Block{T})
    mb = blk.source
    m1 = mb.mb1.source
    m2 = mb.mb2.source
    m1size = size(m1)
    m2size = size(m2)
    restype = promote_type(eltype(m1), eltype(m2))
    res = zeros(restype, m1size[1], m2size[2])

    pmapreduce((t)->(t[1], matrixpart(t[2])*matrixpart(t[3])), (v,t)->begin v[t[1]...] += t[2]; v; end, res, blk)
    res
end

end # module MatOP


text:
@testset "Test OpCode" begin
    @testset "Test Make" begin
        for (op, operands, expected) in [
            (m.OpConstant, [65534], [UInt8(m.OpConstant), 255, 254]),
            (m.OpAdd, [], [UInt8(m.OpAdd)]),
            (m.OpSub, [], [UInt8(m.OpSub)]),
            (m.OpMul, [], [UInt8(m.OpMul)]),
            (m.OpDiv, [], [UInt8(m.OpDiv)]),
            (m.OpEqual, [], [UInt8(m.OpEqual)]),
            (m.OpNotEqual, [], [UInt8(m.OpNotEqual)]),
            (m.OpLessThan, [], [UInt8(m.OpLessThan)]),
            (m.OpGreaterThan, [], [UInt8(m.OpGreaterThan)]),
            (m.OpGetLocal, [255], [UInt8(m.OpGetLocal), 255]),
            (m.OpGetFree, [1], [UInt8(m.OpGetFree), 1]),
            (m.OpGetOuter, [1, 2, 3], [UInt8(m.OpGetOuter), 1, 2, 3]),
            (m.OpClosure, [65534, 255], [UInt8(m.OpClosure), 255, 254, 255]),
            (m.OpIllegal, [], []),
        ]
            @test begin
                instruction = m.make(op, operands...)

                @assert length(instruction) == length(expected) "Instruction has wrong length. Expected $(length(expected)), got $(length(instruction)) instead."

                for (i, e) in enumerate(expected)
                    @assert instruction[i] == e "Instruction has wrong value at index $(i). Expected $(e), got $(instruction[i]) instead."
                end

                true
            end
        end
    end

    @testset "Test Stringifying Instructions" begin
        for (instructions, expected) in [
            (
                [
                    m.make(m.OpConstant, 1),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpConstant, 65535),
                ],
                "0000 OpConstant 1\n0003 OpConstant 2\n0006 OpConstant 65535\n",
            ),
            (
                [m.make(m.OpAdd), m.make(m.OpConstant, 2), m.make(m.OpConstant, 65535)],
                "0000 OpAdd\n0001 OpConstant 2\n0004 OpConstant 65535\n",
            ),
            (
                [m.make(m.OpSub), m.make(m.OpConstant, 2), m.make(m.OpConstant, 65535)],
                "0000 OpSub\n0001 OpConstant 2\n0004 OpConstant 65535\n",
            ),
            (
                [m.make(m.OpMul), m.make(m.OpConstant, 2), m.make(m.OpConstant, 65535)],
                "0000 OpMul\n0001 OpConstant 2\n0004 OpConstant 65535\n",
            ),
            (
                [m.make(m.OpDiv), m.make(m.OpConstant, 2), m.make(m.OpConstant, 65535)],
                "0000 OpDiv\n0001 OpConstant 2\n0004 OpConstant 65535\n",
            ),
            (
                [m.make(m.OpEqual), m.make(m.OpConstant, 2), m.make(m.OpConstant, 65535)],
                "0000 OpEqual\n0001 OpConstant 2\n0004 OpConstant 65535\n",
            ),
            (
                [
                    m.make(m.OpNotEqual),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpConstant, 65535),
                ],
                "0000 OpNotEqual\n0001 OpConstant 2\n0004 OpConstant 65535\n",
            ),
            (
                [
                    m.make(m.OpLessThan),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpConstant, 65535),
                ],
                "0000 OpLessThan\n0001 OpConstant 2\n0004 OpConstant 65535\n",
            ),
            (
                [
                    m.make(m.OpGreaterThan),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpConstant, 65535),
                ],
                "0000 OpGreaterThan\n0001 OpConstant 2\n0004 OpConstant 65535\n",
            ),
            (
                [
                    m.make(m.OpAdd),
                    m.make(m.OpGetLocal, 1),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpConstant, 65535),
                ],
                "0000 OpAdd\n0001 OpGetLocal 1\n0003 OpConstant 2\n0006 OpConstant 65535\n",
            ),
            (
                [
                    m.make(m.OpAdd),
                    m.make(m.OpGetLocal, 1),
                    m.make(m.OpConstant, 2),
                    m.make(m.OpConstant, 65535),
                    m.make(m.OpClosure, 65535, 255),
                ],
                "0000 OpAdd\n0001 OpGetLocal 1\n0003 OpConstant 2\n0006 OpConstant 65535\n0009 OpClosure 65535 255\n",
            ),
            ([m.Instructions([UInt8(m.OpIllegal)])], "ERROR: unknown opcode: 36\n"),
        ]
            @test string(vcat(instructions...)) == expected
        end
    end

    @testset "Test Reading Operands" begin
        for (op, operands, bytes_read) in [
            (m.OpConstant, [65535], 2),
            (m.OpGetLocal, [255], 1),
            (m.OpClosure, [65535, 255], 3),
        ]
            @test begin
                instruction = m.make(op, operands...)
                def = m.lookup(Integer(op))

                @assert !isnothing(def) "Definition for $op not found."

                operands_read, n = m.read_operands(def, instruction[2:end])
                @assert n == bytes_read "Wrong number of bytes read. Expected $bytes_read, got $n instead."
                for (i, want) in enumerate(operands)
                    @assert want == operands_read[i] "Wrong operand at index $i. Expected $(want), got $(operands_read[i]) instead."
                end

                true
            end
        end
    end
end


text:
# ---
# title: 29. Divide Two Integers
# id: problem29
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Binary Search
# link: <https://leetcode.com/problems/divide-two-integers/description/>
# hidden: true
# ---
#
# Given two integers `dividend` and `divisor`, divide two integers without using
# multiplication, division, and mod operator.
#
# Return the quotient after dividing `dividend` by `divisor`.
#
# The integer division should truncate toward zero, which means losing its
# fractional part. For example, `truncate(8.345) = 8` and `truncate(-2.7335) =
# -2`.
#
# **Note:**
#
#   * Assume we are dealing with an environment that could only store integers within the 32-bit signed integer range: [−231,  231 − 1]. For this problem, assume that your function **returns 2 31 − 1 when the division result overflows**.
#
#
#
# **Example 1:**
#
#
#
#     Input: dividend = 10, divisor = 3
#     Output: 3
#     Explanation: 10/3 = truncate(3.33333..) = 3.
#
#
# **Example 2:**
#
#
#
#     Input: dividend = 7, divisor = -3
#     Output: -2
#     Explanation: 7/-3 = truncate(-2.33333..) = -2.
#
#
# **Example 3:**
#
#
#
#     Input: dividend = 0, divisor = 1
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: dividend = 1, divisor = 1
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `-231 <= dividend, divisor <= 231 - 1`
#   * `divisor != 0`
#
#
## @lc code=start
using LeetCode

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


text:
function LogPredict(θ;μ=0,σ=1,ξ=0.5)
       P=function (X)

      m=size(X,1)
      s=zeros(m)
      X=(X.-μ)./σ
      X=[ones(m) X]
      p=Sigmoid(X*θ)
      s[p.>=ξ] .=1
        return s
    end

        return P
end


text:
function [d_time]=GH_greens_funct(rd,dt,nt,c,src,opt)

% Inputs:
% f: frequency [Hz]
% rd: zeros(ns,nr);src and rev distance [m].
% src: zeros(1,nt);usually Ricker.
%
% Outputs:
% d_time: data (Greens function) in the time domain.
%
%  Copyright (C) 2019, Exploration and Environmental Geophysics (EEG) at
%  ETH Zurich

src=reshape(src,[1,nt]);

ff = 1/dt;    % [Hz] Sampling rate
fs = linspace(0,ff,nt); % [Hz] frequency vector

k = 2*pi*fs/c;
if strcmp(opt,'greens1d')
    greens = @(r,f) -1i./(2*k) .* exp(-1i*r*k);
elseif strcmp(opt,'greens2d_1st')
    greens = @(r,f) 1i/4 * besselh(0,2, r*k) .* (-2*pi*fs*1i); % first-order system
elseif strcmp(opt,'greens2d_2nd')
    greens = @(r,f) 1i/4 * besselh(0,2, r*k); % second-order system
elseif strcmp(opt,'greens3d')
    greens = @(r,f) exp(-1i*r*k)./(4*pi*r);
else
    error('opt not supported.\n');
end

%% ---------
SRC=fft(src);

d_f = zeros(length(fs),length(rd));
for l=1:length(rd)
    d_f(:,l) = SRC .* sum( greens(rd(:,l),fs) , 1 );
end
d_f( isnan(d_f) ) = 0;

%% 6. Display in time domain
d_time = ifft( d_f, [], 1, 'symmetric' );

end


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function writeregion(reg,filename)

    tmp = [reg.triclinic];

    nsize = zeros(4,1);
    nsize[1] = length(tmp);
    nsize[2] = length(reg.boxlo[:]);
    nsize[3] = length(reg.boxhi[:]);
    nsize[4] = length(reg.boxtilt[:]);

    fileID = open(filename,"w");
    write(fileID,Float64(length(nsize[:])));
    write(fileID,Float64.(nsize[:]));
    write(fileID,Float64.(tmp[:]));
    write(fileID,Float64.(reg.boxlo[:]));
    write(fileID,Float64.(reg.boxhi[:]));
    write(fileID,Float64.(reg.boxtilt[:]));
    close(fileID);

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


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

const GPU = nothing

function get_imfilter_impl(GPU)
  if GPU ≠ nothing
    imfilter_gpu
  else
    imfilter_cpu
  end
end

function convdist(img::AbstractArray, kern::AbstractArray;
                  weights::AbstractArray=fill(1.0, size(kern)))
  imfilter_impl = get_imfilter_impl(GPU)

  wkern = weights.*kern

  A² = imfilter_impl(img.^2, weights)
  AB = imfilter_impl(img, wkern)
  B² = sum(wkern .* kern)

  parent(@. abs(A² - 2AB + B²))
end

cart2lin(dims, ind) = LinearIndices(dims)[ind]
lin2cart(dims, ind) = CartesianIndices(dims)[ind]

function event!(buff, hard::Dict, tile::CartesianIndices)
  @inbounds for (i, coord) in enumerate(tile)
    h = get(hard, coord, NaN)
    buff[i] = ifelse(isnan(h), 0.0, h)
  end
end

function event(hard::Dict, tile::CartesianIndices)
  buff = Array{Float64}(undef, size(tile))
  event!(buff, hard, tile, def)
  buff
end

function indicator!(buff, hard::Dict, tile::CartesianIndices)
  @inbounds for (i, coord) in enumerate(tile)
    nan = isnan(get(hard, coord, NaN))
    buff[i] = ifelse(nan, false, true)
  end
end

function indicator(hard::Dict, tile::CartesianIndices)
  buff = Array{Bool}(undef, size(tile))
  indicator!(buff, hard, tile)
  buff
end

function activation!(buff, hard::Dict, tile::CartesianIndices)
  @inbounds for (i, coord) in enumerate(tile)
    cond = coord ∈ keys(hard) && isnan(hard[coord])
    buff[i] = ifelse(cond, false, true)
  end
end

function activation(hard::Dict, tile::CartesianIndices)
  buff = Array{Bool}(undef, size(tile))
  activation!(buff, hard, tile)
  buff
end

function preprocess_images(trainimg::AbstractArray{T,N}, soft::AbstractVector,
                           geoconfig::NamedTuple) where {T,N}
  padsize = geoconfig.padsize

  TI = Float64.(trainimg)
  replace!(TI, NaN => 0.)

  SOFT = map(soft) do (aux, auxTI)
    prepend = ntuple(i->0, N)
    append  = padsize .- min.(padsize, size(aux))
    padding = Pad(:symmetric, prepend, append)

    AUX   = Float64.(padarray(aux, padding))
    AUXTI = Float64.(auxTI)
    replace!(AUX, NaN => 0.)
    replace!(AUXTI, NaN => 0.)

    AUX, AUXTI
  end

  TI, SOFT
end

function find_disabled(trainimg::AbstractArray{T,N}, geoconfig::NamedTuple) where {T,N}
  TIsize   = geoconfig.TIsize
  tilesize = geoconfig.tilesize
  distsize = geoconfig.distsize

  disabled = falses(distsize)
  for ind in findall(isnan, trainimg)
    start  = @. max(ind.I - tilesize + 1, 1)
    finish = @. min(ind.I, distsize)
    tile   = CartesianIndex(start):CartesianIndex(finish)
    disabled[tile] .= true
  end

  disabled
end

function find_skipped(hard::Dict, geoconfig::NamedTuple)
  ntiles   = geoconfig.ntiles
  tilesize = geoconfig.tilesize
  spacing  = geoconfig.spacing
  simsize  = geoconfig.simsize

  skipped = Set{Int}()
  datainds = Vector{Int}()
  for tileind in CartesianIndices(ntiles)
    # tile corners are given by start and finish
    start  = @. (tileind.I - 1)*spacing + 1
    finish = @. start + tilesize - 1
    tile   = CartesianIndex(start):CartesianIndex(finish)

    # skip tile if either
    #   1) tile is beyond true simulation size
    #   2) all values in the tile are NaN
    if any(start .> simsize) || !any(activation(hard, tile))
      push!(skipped, cart2lin(ntiles, tileind))
    elseif any(indicator(hard, tile))
      push!(datainds, cart2lin(ntiles, tileind))
    end
  end

  skipped, datainds
end

function genpath(rng::AbstractRNG, extent::Dims{N}, kind::Symbol, datainds::AbstractVector{Int}) where {N}
  path = Vector{Int}()

  if isempty(datainds)
    if kind == :raster
      for ind in LinearIndices(extent)
        push!(path, ind)
      end
    end

    if kind == :random
      path = randperm(rng, prod(extent))
    end

    if kind == :dilation
      nelm = prod(extent)
      pivot = rand(rng, 1:nelm)

      grid = falses(extent)
      grid[pivot] = true
      push!(path, pivot)

      while !all(grid)
        dilated = dilate(grid)
        append!(path, findall(vec(dilated .& .!grid)))
        grid = dilated
      end
    end
  else
    # data-first path
    shuffle!(rng, datainds)

    grid = falses(extent)
    for pivot in datainds
      grid[pivot] = true
      push!(path, pivot)
    end

    while !all(grid)
      dilated = dilate(grid)
      append!(path, findall(vec(dilated .& .!grid)))
      grid = dilated
    end
  end

  path
end


text:
makepath(p) = joinpath("subdir", p)

@testset "included" begin
    @test true
end

# test that ReTest can handle non-literal arguments in include
SUB="sub"
include(makepath("$SUB.jl"))

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
module HistogramThresholding

using LinearAlgebra
using ImageBase

# TODO: port ThresholdAPI to ImagesAPI
include("ThresholdAPI/ThresholdAPI.jl")
import .ThresholdAPI: AbstractThresholdAlgorithm,
                      find_threshold, build_histogram

include("common.jl")
include("algorithms/otsu.jl")
include("algorithms/yen.jl")
include("algorithms/minimum_error.jl")
include("algorithms/unimodal.jl")
include("algorithms/moments.jl")
include("algorithms/minimum.jl")
include("algorithms/intermodes.jl")
include("algorithms/balancedthreshold.jl")
include("algorithms/entropy_thresholding.jl")
include("deprecations.jl")

export
  # main functions
  find_threshold,
  build_histogram,
  Otsu,
  UnimodalRosin,
  Moments,
  MinimumIntermodes,
  Intermodes,
  MinimumError,
  Balanced,
  Yen,
  Entropy
end # module

max_stars_repo_name: zygmuntszpak/HistogramThresholding

max_forks_repo_name: zygmuntszpak/HistogramThresholding


text:
using TSML
using TSMLextra
using TSML.ArgumentParsers

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
  tsmlmain()
end

max_stars_repo_name: ppalmes/TSMLextra

max_issues_repo_name: ppalmes/TSMLextra

max_forks_repo_name: ppalmes/TSMLextra


text:
# TODO: I cound that "utils.jl" has been `include`ed in the main `VectorSphericalWaves` module. Is there a more neat way? when I `include("utils.jl")` here, I get lots of warnings.

# exporting functions that calculate for all m,n combinations
export M_N_wave_all_m_n
export B_C_P_mn_of_θ_ϕ_for_all_m_n
export πₘₙ_τₘₙ_all_all_m_n

#############################################################################################
# This version is different from "VectorSphericalHarmonics.jl", as it calculate VSWF for all m,n and all kr, θ, ϕ in one call
# This should be more stable and faster, as it employs recurrence relations in calculating πₘₙ and τₘₙ
# using recurrence relations in calculating πₘₙ and τₘₙ

#############################################################################################
# calculate π(θ) and τ(θ) using recurrence relations
function πₘₙ_τₘₙ_all_all_m_n(n_max::I, θ::R; verbose=false) where {R <: Real, I <: Integer}
    # calculate A with recurrence relation
    A = SortedDict(0 => 1.0)
    for m = 0:n_max - 1
        A[m + 1] = A[m] * sqrt((2m + 1) / (2 * (m + 1)))
    end

    # calculate πₘₙ with recurrence relation
    πₘₙ_all = zeros(get_max_single_index_from_n_max(n_max))
    for m = 1:n_max   # TODO: I think I need to start m from 0, not 1.
        n = m
        πₘₙ_all[single_index_from_m_n(m, n)] = m * A[m] * sin(θ)^(m - 1)
        if verbose; println("m=$m, n=$n, πₘₙ_all=$(πₘₙ_all[single_index_from_m_n(m, n)])"); end

        # calculate π₋ₘₙ from πₘₙ
        if m != 0
            πₘₙ_all[single_index_from_m_n(-m, n)] = (-1)^(m + 1) * πₘₙ_all[single_index_from_m_n(m, n)]
            if verbose; println("m=$(-m), n=$n, πₘₙ_all=$(πₘₙ_all[single_index_from_m_n(-m, n)])"); end
        end

        for n = (m + 1):n_max
            if n == m + 1
                πₘₙ_all[single_index_from_m_n(m, n)] = 1 / sqrt(n^2 - m^2) * ((2n - 1) * cos(θ) * πₘₙ_all[single_index_from_m_n(m, n - 1)])
                if verbose; println("m=$m, n=$n, πₘₙ_all=$(πₘₙ_all[single_index_from_m_n(m, n)])"); end
            else
                πₘₙ_all[single_index_from_m_n(m, n)] = 1 / sqrt(n^2 - m^2) * ((2n - 1) * cos(θ) * πₘₙ_all[single_index_from_m_n(m, n - 1)]) - sqrt((n - 1)^2 - m^2) * πₘₙ_all[single_index_from_m_n(m, n - 2)]
                if verbose; println("m=$m, n=$n, πₘₙ_all=$(πₘₙ_all[single_index_from_m_n(m, n)])"); end
            end

            # calculate π₋ₘₙ from πₘₙ
            if m != 0
                πₘₙ_all[single_index_from_m_n(-m, n)] = (-1)^(m + 1) * πₘₙ_all[single_index_from_m_n(m, n)]
                if verbose; println("m=$(-m), n=$n, πₘₙ_all=$(πₘₙ_all[single_index_from_m_n(-m, n)])"); end
            end
        end
    end
    τₘₙ_all = 0 # TODO: add the code for it
    return πₘₙ_all, τₘₙ_all
end

#############################################################################################
# Legendre and Associated Legendre
function Legendre_polynomials_Pn_array(n_max::I, x::NN) where {I <: Integer, NN <: Number}
    """
    Calculate all Legendre polynomials Pₙ(x) from n = 1 up to n=n_max using recurrence relation
    https://en.wikipedia.org/wiki/Legendre_polynomials

    returns
    dictionary. TODO: find a better way if this cause adjoint calculation trouble
    """
    P = SortedDict(
        0 => 1,
        1 => x,
    )

    for n = 1:(n_max - 1)
        P[n + 1] = ((2n + 1) * x * P[n] - n * P[n - 1]) / (n + 1)
    end

    return P
end

function Legendre_polynomials_Pn(n::I, x::NN) where {I <: Integer, NN <: Number}
    """
    Calculate Legendre polynomials Pₙ(x) at a given n
    """
    return Legendre_polynomials_Pn_array(n, x)[n]
end

function Associated_Legendre_polynomials_Pmn_array(m::I, n_max::I, x) where {I <: Integer}
    """
    Associated Legendre polynomials Pᵐₙ(x) from n = m up to n=n_max using recurrence relation
    https://en.wikipedia.org/wiki/Associated_Legendre_polynomials#Recurrence_formula

    returns
    dictionary. TODO: find a better way if this cause adjoint calculation trouble
    """
    n = m
    P = SortedDict(
        n => (-1)^n * factorial(factorial(2n - 1)) * (1 - x^2)^(n / 2),
    )
    P[n + 1] = x * (2n + 1) * P[n]

    for n = (m + 1):(n_max - 1)
        P[n + 1] = ( (2n + 1) * x * P[n] - (n + m) * P[n - 1] ) / (n - m + 1)
    end

    return P
end

function Associated_Legendre_polynomials_Pmn(m::I, n::I, x) where {I <: Integer}
    """
    Associated Legendre polynomials Pᵐₙ(x) at a given n
    """
    return Associated_Legendre_polynomials_Pmn_array(m, n, x)[n]
end

#############################################################################################
# Wigner-d, using recurrence
# TODO: I think this will not work for large s,m,n values, try to fix it as in `wignerdjmn_ELZOUKA`
function wignerd_and_∂wignerd_for_all_s(s_max::I, m::I, n::I, θ::R; get_derivatives=true, verbose=false) where {R <: Real, I <: Integer}
    """
    Calculate dˢₘₙ(θ) and ∂(dˢₘₙ(θ))/∂θ for all values of s, where s starts from s_min up to s_max. s_min is the maximum of |m| and |n|

    """
    # TODO: special case of m=0, n=0 (eq. B.27)
    s_min = max(abs(m), abs(n))

    # calculate ξ from eq. B.16
    if n >= m
        ξ_mn = 1
    else
        ξ_mn = (-1)^(m - n)
    end

    x = cos(θ)

    # calculate d^s_min__m_n(θ) from eq. B.24
    d_smin_m_n =
        ξ_mn * 2.0^(-s_min) * sqrt(
            factorial(BigInt(2s_min)) / (factorial(BigInt(abs(m - n))) * factorial(BigInt(abs(m + n))))
        ) *
        (1 - x)^(abs(m - n) / 2) *
        (1 + x)^(abs(m + n) / 2)

    d = SortedDict(
        s_min - 1 => 0.0,
        s_min   => d_smin_m_n
    )

    # applying the recurrence relation
    if n == 0
        if m == 0
            if verbose; println("m=$m, n=$n, I will use Legendre_polynomials_Pn_array"); end
            P_s = Legendre_polynomials_Pn_array(s_max + 1, x)
            for s = s_min:s_max + 1
                d[s] = P_s[s]
            end
        else
            if verbose; println("m=$m, n=$n, I will use Associated_Legendre_polynomials_Pmn_array"); end
            P_m_s = Associated_Legendre_polynomials_Pmn_array(m, s_max + 1, x)
            for s = s_min:s_max + 1
                d[s] = sqrt(factorial(s - m) / factorial(s + m)) * P_m_s[s]
            end
        end
    else
        for s = s_min:s_max + 1
            if verbose; println("m=$m, n=$n, I will use the general recurrence"); end
            d[s + 1] = 1 / (s * sqrt((s + 1)^2 - m^2) * sqrt((s + 1)^2 - n^2)) * (
                (2s + 1) * (s * (s + 1) * x - m * n) * d[s]
                - 1 * (s + 1) * sqrt(s^2 - m^2) * sqrt(s^2 - n^2) * d[s - 1]
            ) # eq. B.22
        end
    end

    # calculate the derivative ∂(dˢₘₙ(θ))/∂θ
    if get_derivatives
        ∂d_∂θ = SortedDict(
            s_min - 1 => 0.0,
            s_min   => 0.0
        )

        sin_theta = sin(θ)

        for s = s_min:s_max
            if verbose; println("s=$s"); end
            ∂d_∂θ[s] = 1 / sin_theta * (
                -1 * ((s + 1) * sqrt((s^2 - m^2) * (s^2 - n^2))) / (s * (2s + 1)) * d[s - 1]
                - 1 * (m * n) / (s * (s + 1)) * d[s]
                + 1 * (
                        s *
                        sqrt((s + 1)^2 - m^2) *
                        sqrt((s + 1)^2 - n^2)
                    ) /
                    ((s + 1) * (2s + 1)) * d[s + 1]
            )
        end
        return d, ∂d_∂θ
    else
        return d
    end

end

#############################################################################################
# calculate π(θ) and τ(θ) using Wigner-d that was calculated using recurrence relations
function πₘₙ_τₘₙ_all_all_m_n_using_wigner(n_max::I, θ::R; verbose=false) where {R <: Real, I <: Integer}
    πₘₙ_all = zeros(get_max_single_index_from_n_max(n_max))
    τₘₙ_all = zeros(get_max_single_index_from_n_max(n_max))

    for m = 0:n_max  # TODO: special case of m=0
        d_rec_all_n, ∂d_∂θ_rec_all_n  = wignerd_and_∂wignerd_for_all_s(n_max, 0, m, θ)
        for n = m:n_max
            if n != 0
                if verbose; println("m=$m, n=$n, calculate πₘₙ_all, τₘₙ_all"); end
                πₘₙ_all[single_index_from_m_n(m, n)] = m / sin(θ) * d_rec_all_n[n]
                τₘₙ_all[single_index_from_m_n(m, n)] = ∂d_∂θ_rec_all_n[n]

                if m != 0
                    πₘₙ_all[single_index_from_m_n(-m, n)] = (-1)^(m + 1) * πₘₙ_all[single_index_from_m_n(m, n)]
                    τₘₙ_all[single_index_from_m_n(-m, n)] = (-1)^(m)   * τₘₙ_all[single_index_from_m_n(m, n)]
                end
            end
        end
    end
    return πₘₙ_all, τₘₙ_all
end

function B_C_mn_of_θ_for_all_m_n(n_max::I, θ::R) where {R <: Real, I <: Integer}
    """
    Calculate Bₙₘ(θ), Cₙₘ(θ), Pₙₘ(θ) equations C.19, C.20, C.21
    The order of m,n is according to the function "single_index_from_m_n"
    """
    πₘₙ_all, τₘₙ_all = πₘₙ_τₘₙ_all_all_m_n_using_wigner(n_max, θ)
    B = (_ -> zero(SVector{3,Complex})).(πₘₙ_all)
    C = (_ -> zero(SVector{3,Complex})).(πₘₙ_all)
    for idx in eachindex(πₘₙ_all)
        B[idx] = [0,      τₘₙ_all[idx], im * πₘₙ_all[idx] ]
        C[idx] = [0, im * πₘₙ_all[idx], -1 * τₘₙ_all[idx] ]
    end

    return B, C
end

function P_mn_of_θ_for_all_m_n(n_max::I, θ::R) where {R <: Real, I <: Integer}
    """
    Calculate Pₙₘ(θ) using equations C.19, C.20
    The order of m,n is according to the function "single_index_from_m_n"
    """
    P = fill(zero(SVector{3,Complex}), get_max_single_index_from_n_max(n_max))
    for m = 0:n_max
        d_rec_all_n  = wignerd_and_∂wignerd_for_all_s(n_max, 0, m, θ; get_derivatives=false)
        for n = m:n_max
            if n != 0
                P[single_index_from_m_n(m, n)] = [d_rec_all_n[n], 0, 0]
                P[single_index_from_m_n(-m, n)] = [(-1)^m * d_rec_all_n[n], 0, 0]  # using symmetry relation B.5, we can get d_(0,-m) from d_(0,m)
            end
        end
    end
    return P
end

function B_C_P_mn_of_θ_ϕ_for_all_m_n(n_max::I, θ::R, ϕ::R) where {R <: Real, I <: Integer}
    """
    Calculate Bₙₘ(θ,ϕ), Cₙₘ(θ,ϕ), Pₙₘ(θ,ϕ) for all m and n
    """
    B_of_θ, C_of_θ = B_C_mn_of_θ_for_all_m_n(n_max, θ)
    P_of_θ = P_mn_of_θ_for_all_m_n(n_max, θ)

    for n = 1:n_max
        for m = -n:n
            factor = (-1)^m * convert(R, sqrt_factorial_n_plus_m_over_factorial_n_minus_m(m,n)) * exp(im * m * ϕ)
            B_of_θ[single_index_from_m_n(m, n)] *= factor
            C_of_θ[single_index_from_m_n(m, n)] *= factor
            P_of_θ[single_index_from_m_n(m, n)] *= factor
        end
    end

    return B_of_θ, C_of_θ, P_of_θ
end

function M_N_wave_all_m_n(n_max::I, kr::NN, θ::R, ϕ::R; kind="regular") where {R <: Real, I <: Integer, NN <: Number}
    """
    Parameters
    ==========
    kind: string, either ["regular" or "incoming"] or ["irregular" or "outgoing"]
    """
    radial_function, radial_function_special_derivative  = get_radial_function_and_special_derivative_given_kind(kind)

    B_of_θ_ϕ, C_of_θ_ϕ, P_of_θ_ϕ  = B_C_P_mn_of_θ_ϕ_for_all_m_n(n_max, θ, ϕ)

    M = 0 .* B_of_θ_ϕ
    N = 0 .* B_of_θ_ϕ

    for n = 1:n_max
        for m = -n:n
            M[single_index_from_m_n(m, n)] = convert(R, γ_mn(m, n)) * radial_function(n, kr) * C_of_θ_ϕ[single_index_from_m_n(m, n)]
            N[single_index_from_m_n(m, n)] = convert(R, γ_mn(m, n)) * (
                n * (n + 1) / kr * radial_function(n, kr)    * P_of_θ_ϕ[single_index_from_m_n(m, n)]
                + (radial_function_special_derivative(n, kr) * B_of_θ_ϕ[single_index_from_m_n(m, n)])
            )
        end
    end
    return M, N
end


text:
module FrostHelperFallingBlockIgnoreSolids

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/FallingBlockIgnoreSolids" FallingBlockIgnoreSolids(x::Integer, y::Integer, width::Integer=16, height::Integer=16, tiletype::String="3", climbFall::Bool=true, behind::Bool=false)

const placements = Ahorn.PlacementDict(
    "Falling Block (Ignore Solids, Frost Helper)" => Ahorn.EntityPlacement(
        FallingBlockIgnoreSolids,
        "rectangle",
        Dict{String, Any}(),
        Ahorn.tileEntityFinalizer
    ),
)

Ahorn.editingOptions(entity::FallingBlockIgnoreSolids) = Dict{String, Any}(
    "tiletype" => Ahorn.tiletypeEditingOptions()
)

Ahorn.minimumSize(entity::FallingBlockIgnoreSolids) = 8, 8
Ahorn.resizable(entity::FallingBlockIgnoreSolids) = true, true

Ahorn.selection(entity::FallingBlockIgnoreSolids) = Ahorn.getEntityRectangle(entity)

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::FallingBlockIgnoreSolids, room::Maple.Room) = Ahorn.drawTileEntity(ctx, room, entity)

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
using Printf
using Dates: AbstractTime

using Oceananigans.Units

maybe_int(t) = isinteger(t) ? Int(t) : t

"""
    prettytime(t, longform=true)

Convert a floating point value `t` representing an amount of time in
SI units of seconds to a human-friendly string with three decimal places.
Depending on the value of `t` the string will be formatted to show `t` in
nanoseconds (ns), microseconds (μs), milliseconds (ms),
seconds, minutes, hours, days, or years.

With `longform=false`, we use s, m, hrs, d, and yrs in place of seconds,
minutes, hours, and years.
"""
function prettytime(t, longform=true)
    # Modified from: https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/src/trials.jl

    # Some shortcuts
    s = longform ? "seconds" : "s"
    iszero(t) && return "0 $s"
    t < 1e-9 && return @sprintf("%.3e %s", t, s) # yah that's small

    t = maybe_int(t)
    value, units = prettytimeunits(t, longform)

    if isinteger(value)
        return @sprintf("%d %s", value, units)
    else
        return @sprintf("%.3f %s", value, units)
    end
end

function prettytimeunits(t, longform=true)
    t < 1e-9 && return t, "" # just _forget_ picoseconds!
    t < 1e-6 && return t * 1e9, "ns"
    t < 1e-3 && return t * 1e6, "μs"
    t < 1    && return t * 1e3, "ms"
    if t < minute
        value = t
        !longform && return value, "s"
        units = value == 1 ? "second" : "seconds"
        return value, units
    elseif t < hour
        value = maybe_int(t / minute)
        !longform && return value, "m"
        units = value == 1 ? "minute" : "minutes"
        return value, units
    elseif t < day
        value = maybe_int(t / hour)
        units = value == 1 ? (longform ? "hour" : "hr") :
                             (longform ? "hours" : "hrs")
        return value, units
    elseif t < year
        value = maybe_int(t / day)
        !longform && return value, "d"
        units = value == 1 ? "day" : "days"
        return value, units
    else
        value = maybe_int(t / year)
        units = value == 1 ? (longform ? "year" : "yr") :
                             (longform ? "years" : "yrs")
        return value, units
    end
end

prettytime(dt::AbstractTime) = "$dt"


text:
"""
    dirichlet_beta()

 Calculates Dirichlet beta function,
     [https://en.wikipedia.org/wiki/Dirichlet_beta_function](https://en.wikipedia.org/wiki/Dirichlet_beta_function)

## Arguments
* ``s`` `::Number`: it should work for any type of number, but mainly tested for `Complex{Float64}`

## Examples
```jldoctest; setup = :(using Polylogarithms)
julia> dirichlet_beta(1.5)
0.8645026534612017
```
"""
function dirichlet_beta(s::Number)
    β = 4.0^(-s) * ( SpecialFunctions.zeta(s,0.25) - SpecialFunctions.zeta(s,0.75) )
end


text:
using ReachabilityBenchmarks, MathematicalSystems, Reachability, Plots

include(@relpath "linearSwitching_model.jl")
include(@relpath "linearSwitching_specifications.jl")

S = linearSwitching_model()
X0, options = linearSwitching_specification()
options = Options(options)

# initial value problem
problem = InitialValueProblem(S, X0)

# reachability algorithm
algorithm_continuous = BFFPS19(:δ => 1e-4, :partition=>[1:2, 3:3, 4:4, 5:5])
algorithm_hybrid = DecomposedDiscretePost(:out_vars=>[1, 2], :clustering=>:none)

# compute flowpipe
options[:mode] = "reach"
solution = solve(problem, options, algorithm_continuous, algorithm_hybrid)

# project flowpipe to x₁ and x₂
solution.options[:plot_vars] = [1, 2]
solution_proj = project(solution)

# plot projection
plot(solution_proj)
savefig("linearSwitching")

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
function str_quote(str::AbstractString, delim::Char, na::AbstractString)::String
    do_quote = false
    for c in str
        if c in [QUOTE_CHAR, delim, '\n', '\r']
            do_quote = true
            break
        end
    end

    do_quote || str == na || return str

    io = IOBuffer(sizehint=sizeof(str)+8)
    write(io, QUOTE_CHAR)
    for c in str
        write(io, c)
        if c == QUOTE_CHAR
            write(io, QUOTE_CHAR)
            had_q = true
        end
    end
    write(io, QUOTE_CHAR)
    String(take!(io))
end

quote_with_missing(val, delim, na) =
    ismissing(val) ? na : str_quote(string(val), delim, na)

"""
    write_csv(filename::AbstractString, df::DataFrame;
              delim::Char=',', na::AbstractString="")

    Writes `df` to disk as `filename` CSV file using `delim` separator
    and `na` as string for representing missing value.

    `"` character is used for quoting fields. In quoted field use `""` to
    represent `"`.

    `na` is written to CSV verbatim when `DataFrame` contains `missing`.
    If a string is equal to `na` then it is written to CSV quoted as `"na"`.
"""
function write_csv(filename::AbstractString, df::DataFrame;
                   delim::Char=',', na::AbstractString="")
    if delim in [QUOTE_CHAR, '\n', '\r']
        throw(ArgumentError("delim is a quote char or a newline"))
    end
    if any(occursin.([QUOTE_CHAR, delim, '\n', '\r'], na))
        throw(ArgumentError("na contains quote, separator or a newline"))
    end
    open(filename, "w") do io
        if length(names(df)) > 0
            println(io, join(str_quote.(string.(names(df)), delim, na), delim))
            for i in 1:nrow(df)
                line = [quote_with_missing(df[i,j], delim, na) for j in 1:ncol(df)]
                println(io, join(line, delim))
            end
        end
    end
end


text:
__precompile__()
module BufferedStream

import Base: getindex, length, push!, eltype, linearindexing, size, similar, show
export LinkedStream, StreamView
export shift_origin!

# The foundation of a linked-list of chunks
type StreamChunk{T}
    data::Array{T,1}
    next::Nullable{StreamChunk{T}}

    StreamChunk() = new(T[], Nullable{StreamChunk{T}}())
    StreamChunk(data::Array{T,1}) = new(data, Nullable{StreamChunk{T}}())
end

# Just return the length of this chunk.
length(chunk::StreamChunk) = length(chunk.data)

# Show something well-behaved
function show(io::IO, chunk::StreamChunk)
    write(io, "StreamChunk")
    Base.showlimited(io, chunk.data)
    if isnull(chunk.next)
        write(io, " -> <null>")
    else
        write(io, " -> StreamChunk")
    end
end

# The object that keeps track of the last chunk in the linked list
type LinkedStream{T}
    # This is where our lovely loving data gets stored
    latest::StreamChunk{T}

    # How many samples have we discarded?
    samples_past::Int64

    # We can only create one that initializes with an empty list
    LinkedStream() = new(StreamChunk{T}(), 0)
end

# Pushes a new chunk onto the end of the stream.  Julia should automatically garbage collect
# the chunks any views have moved past, since they will no longer be accessible.
function push!{T}(stream::LinkedStream{T}, data::Array{T,1})
    # Create a new chunk based off of this data
    chunk = StreamChunk{T}(data)

    # Account for the samples we're putting behind us
    stream.samples_past += length(stream.latest)

    # Move the new chunk into the lead position
    stream.latest.next = Nullable(chunk)
    stream.latest = chunk
    return
end

function show{T}(io::IO, stream::LinkedStream{T})
    write(io, "LinkedStream{$(string(T))}, $(stream.samples_past) samples post-origin")
end

# A view that acts like an array, but that you can move forward in the stream
type StreamView{T} <: DenseArray{T, 1}
    # Where is our origin?
    chunk::StreamChunk{T}
    sample_offset::Int64

    # What is our origin's absolute index?
    abs_idx::Int64
end

# Define stuff for AbstractArray
linearindexing(::StreamView) = Base.LinearFast()
# Calculates how many samples after our origin there are buffered up
function length(view::StreamView)
    # Count total length of all chunks after the origin of this guy
    len = length(view.chunk)
    chunk = view.chunk
    while !isnull(chunk.next)
        chunk = chunk.next.value
        len += length(chunk)
    end

    # Subtract sample_offset, and ensure that zero is the smallest we will ever return
    return max(len - view.sample_offset, 0)
end

size(view::StreamView) = (length(view),)
eltype{T}(view::StreamView{T}) = T

# Construct a new view off of a stream, with an optional offset
function StreamView{T}(stream::LinkedStream{T})
    offset = length(stream.latest)
    return StreamView{T}(stream.latest, offset, stream.samples_past + offset)
end

# Construct a new view off of another view (easy peasey; just a straight copy)
StreamView{T}(v::StreamView{T}) = StreamView{T}(v.chunk, v.sample_offset, v.abs_idx)

# Move this view forward `offset` samples; we only support positive offsets!
function shift_origin!(view::StreamView, offset::Integer)
    if offset < 0
        throw(DomainError())
    end

    # Bump up the absolute index immediately
    view.abs_idx += offset

    # Move through chunks until offset lands within this current chunk, or we run out of chunks
    while (offset >= length(view.chunk) - view.sample_offset) && !isnull(view.chunk.next)
        offset -= length(view.chunk) - view.sample_offset
        # sample_offset can be greater than length(view.chunk), if we had shifted past the end of
        # buffered data previously; this will move us forward as far as we can
        view.sample_offset = max(view.sample_offset - length(view.chunk), 0)
        view.chunk = view.chunk.next.value
    end

    # Once our offset lands within this chunk, update the sample_offset and call it good!  If we are
    # past the end of the current offset (e.g. we are shifting into the great unknown) do the same,
    # simply leaving sample_offset greater than length(view.chunk)
    view.sample_offset += offset
    return
end

# If someone asks for an index, serve it if we can!
function getindex(view::StreamView, idx::Integer)
    # Scoot us forward as much as the sample_offset requires
    idx += view.sample_offset

    # Zoom through chunks as much as we need
    chunk = view.chunk
    while idx > length(chunk)
        # If we can't move forward anymore, throw a BoundsError!
        if isnull(chunk.next)
            throw(BoundsError())
        end
        idx -= length(chunk)
        chunk = chunk.next.value
    end

    # If we got to the correct chunk, index into it and yield the data!
    return chunk.data[idx]
end

# Don't define setindex!() for now; that may or may not be a great idea.

end # module