text:
module LevelDB

export open_db
export close_db
export create_write_batch
export batch_put
export write_batch
export db_put
export db_get
export db_delete
export db_range
export range_close

function open_db(file_path, create_if_missing)
    options = ccall( (:leveldb_options_create, "libleveldb"), Ptr{Void}, ())
    if create_if_missing
        ccall( (:leveldb_options_set_create_if_missing, "libleveldb"), Void,
              (Ptr{Void}, Uint8), options, 1)
    end
    err = Ptr{Uint8}[0]
    db = ccall( (:leveldb_open, "libleveldb"), Ptr{Void},
               (Ptr{Void}, Ptr{Uint8}, Ptr{Ptr{Uint8}}) , options, file_path, err)

    if db == C_NULL
        error(bytestring(err[1]))
    end
    return db
end

function close_db(db)
    ccall( (:leveldb_close, "libleveldb"), Void, (Ptr{Void},), db)
end

function db_put(db, key, value, val_len)
    options = ccall( (:leveldb_writeoptions_create, "libleveldb"), Ptr{Void}, ())
    err = Ptr{Uint8}[0]
    ccall( (:leveldb_put, "libleveldb"), Void,
          (Ptr{Void}, Ptr{Void}, Ptr{Void}, Uint, Ptr{Uint8}, Uint, Ptr{Ptr{Uint8}} ),
          db, options,key, length(key), value, val_len, err)
    if err[1] != C_NULL
        error(bytestring(err[1]))
    end
end

# return an Uint8 array obj
function db_get(db, key)
    # leveldb_get will allocate the buffer for return value
    options = ccall( (:leveldb_readoptions_create, "libleveldb"), Ptr{Void}, ())
    err = Ptr{Uint8}[0]
    val_len = Csize_t[0]
    value = ccall( (:leveldb_get, "libleveldb"), Ptr{Uint8},
          (Ptr{Void}, Ptr{Void}, Ptr{Uint8}, Uint, Ptr{Csize_t},  Ptr{Ptr{Uint8}} ),
          db, options, key, length(key), val_len, err)
    if err[1] != C_NULL
        error(bytestring(err[1]))
    else
        s = pointer_to_array(value, (val_len[1],), true)
        s
    end
end

function db_delete(db, key)
    options = ccall( (:leveldb_writeoptions_create, "libleveldb"), Ptr{Void}, ())
    err = Ptr{Uint8}[0]
    ccall( (:leveldb_delete, "libleveldb"), Void,
          (Ptr{Void}, Ptr{Void}, Ptr{Void}, Uint, Ptr{Ptr{Uint8}} ),
          db, options, key, length(key), err)
    if err[1] != C_NULL
        error(bytestring(err[1]))
    end
end

function create_write_batch()
    batch = ccall( (:leveldb_writebatch_create, "libleveldb"), Ptr{Void},())
    return batch
end

function batch_put(batch, key, value, val_len)
    ccall( (:leveldb_writebatch_put, "libleveldb"), Void,
          (Ptr{Uint8}, Ptr{Uint8}, Uint, Ptr{Uint8}, Uint),
          batch, key, length(key), value, val_len)
end

function write_batch(db, batch)
    options = ccall( (:leveldb_writeoptions_create, "libleveldb"), Ptr{Void}, ())
    err = Ptr{Uint8}[0]
    ccall( (:leveldb_write, "libleveldb"), Void,
          (Ptr{Void}, Ptr{Void}, Ptr{Void},  Ptr{Ptr{Uint8}} ),
          db, options, batch, err)
    if err[1] != C_NULL
        error(bytestring(err[1])) */
    end
end

function create_iter(db::Ptr{Void}, options::Ptr{Void})
  ccall( (:leveldb_create_iterator, "libleveldb"), Ptr{Void},
              (Ptr{Void}, Ptr{Void}),
              db, options)
end

function iter_valid(it::Ptr{Void})
  ccall( (:leveldb_iter_valid, "libleveldb"), Uint8,
    (Ptr{Void},),
    it) == 1
end

function iter_key(it::Ptr{Void})
  k_len = Csize_t[0]
  key = ccall( (:leveldb_iter_key, "libleveldb"), Ptr{Uint8},
    (Ptr{Void}, Ptr{Csize_t}),
    it, k_len)
  bytestring(key, k_len[1])
end

function iter_value(it::Ptr{Void})
  v_len = Csize_t[0]
  value = ccall( (:leveldb_iter_value, "libleveldb"), Ptr{Uint8},
    (Ptr{Void}, Ptr{Csize_t}),
    it, v_len)
  pointer_to_array(value, (v_len[1],), false)
end

function iter_seek(it::Ptr{Void}, key)
  ccall( (:leveldb_iter_seek, "libleveldb"), Void,
    (Ptr{Void}, Ptr{Uint8}, Uint),
    it, key, length(key))
end

function iter_next(it::Ptr{Void})
  ccall( (:leveldb_iter_next, "libleveldb"), Void,
    (Ptr{Void},),
    it)
end

type Range
  iter::Ptr{Void}
  options::Ptr{Void}
  key_start::String
  key_end::String
  destroyed::Bool
end

function db_range(db, key_start, key_end="\uffff")
  options = ccall( (:leveldb_readoptions_create, "libleveldb"), Ptr{Void}, ())
  iter = create_iter(db, options)
  Range(iter, options, key_start, key_end, false)
end

function range_close(range::Range)
  if !range.destroyed
    range.destroyed = true
    ccall( (:leveldb_iter_destroy, "libleveldb"), Void,
      (Ptr{Void},),
      range.iter)
    ccall( (:leveldb_readoptions_destroy, "libleveldb"), Void,
      (Ptr{Void},),
      range.options)
  end
end

function Base.start(range::Range)
  iter_seek(range.iter, range.key_start)
end

function Base.done(range::Range, state=None)
  if range.destroyed
    return true
  end
  isdone = iter_valid(range.iter) ? iter_key(range.iter) > range.key_end : true
  if isdone
    range_close(range)
  end
  isdone
end

function Base.next(range::Range, state=None)
  k = iter_key(range.iter)
  v = iter_value(range.iter)
  iter_next(range.iter)
  ((k, v), None)
end

end


text:
# returns all variables of a cost function, in terms of appearance
extract_variables(t::TermOrExpr) = variables(t)

function extract_variables(t::NTuple{N,TermOrExpr}) where {N}
  x = variables.(t)
  xAll = x[1]
  for i = 2:length(x)
    for xi in x[i]
      if (xi in xAll) == false
        xAll = (xAll...,xi)
      end
    end
  end
  return xAll
end

# extract functions from terms
function extract_functions(t::Term)
  f = displacement(t) == 0 ? t.f : PrecomposeDiagonal(t.f, 1.0, displacement(t)) #for now I keep this
  f = t.lambda == 1. ? f : Postcompose(f, t.lambda)                                  #for now I keep this
  #TODO change this
  return f
end
extract_functions(t::NTuple{N,Term}) where {N} = SeparableSum(extract_functions.(t))
extract_functions(t::Tuple{Term}) = extract_functions(t[1])

# extract functions from terms without displacement
function extract_functions_nodisp(t::Term)
  f = t.lambda == 1. ? t.f : Postcompose(t.f, t.lambda)
  return f
end
extract_functions_nodisp(t::NTuple{N,Term}) where {N} = SeparableSum(extract_functions_nodisp.(t))
extract_functions_nodisp(t::Tuple{Term}) = extract_functions_nodisp(t[1])

# extract operators from terms

# returns all operators with an order dictated by xAll

#single term, single variable
extract_operators(xAll::Tuple{Variable}, t::TermOrExpr)  = operator(t)
extract_operators(xAll::NTuple{N,Variable}, t::TermOrExpr) where {N} = extract_operators(xAll, (t,))

#multiple terms, multiple variables
function extract_operators(xAll::NTuple{N,Variable}, t::NTuple{M,TermOrExpr}) where {N,M}
  ops = ()
  for ti in t
    tex = expand(xAll,ti)
    ops = (ops...,sort_and_extract_operators(xAll,tex))
  end
  return vcat(ops...)
end

sort_and_extract_operators(xAll::Tuple{Variable}, t::TermOrExpr) = operator(t)

function sort_and_extract_operators(xAll::NTuple{N,Variable}, t::TermOrExpr) where {N}
  p = zeros(Int,N)
  xL = variables(t)
  for i in eachindex(xAll)
    p[i] = findfirst( xi -> xi == xAll[i], xL)
  end
  return operator(t)[p]
end

# extract affines from terms

# returns all affines with an order dictated by xAll

#single term, single variable
extract_affines(xAll::Tuple{Variable}, t::TermOrExpr)  = affine(t)

extract_affines(xAll::NTuple{N,Variable}, t::TermOrExpr) where {N} = extract_affines(xAll, (t,))

#multiple terms, multiple variables
function extract_affines(xAll::NTuple{N,Variable}, t::NTuple{M,TermOrExpr}) where {N,M}
  ops = ()
  for ti in t
    tex = expand(xAll,ti)
    ops = (ops...,sort_and_extract_affines(xAll,tex))
  end
  return vcat(ops...)
end

sort_and_extract_affines(xAll::Tuple{Variable}, t::TermOrExpr) = affine(t)

function sort_and_extract_affines(xAll::NTuple{N,Variable}, t::TermOrExpr) where {N}
  p = zeros(Int,N)
  xL = variables(t)
  for i in eachindex(xAll)
    p[i] = findfirst( xi -> xi == xAll[i], xL)
  end
  return affine(t)[p]
end

# expand term domain dimensions
function expand(xAll::NTuple{N,Variable}, t::Term) where {N}
  xt   = variables(t)
  C    = codomainType(operator(t))
  size_out = size(operator(t),1)
  ex = t.A

  for x in xAll
    if !( x in variables(t) )
      ex += Zeros(eltype(~x),size(x),C,size_out)*x
    end
  end
  return Term(t.lambda, t.f, ex)
end

function expand(xAll::NTuple{N,Variable}, ex::AbstractExpression) where {N}
  ex = convert(Expression,ex)
  xt   = variables(ex)
  C    = codomainType(operator(ex))
  size_out = size(operator(ex),1)

  for x in xAll
    if !( x in variables(ex) )
      ex += Zeros(eltype(~x),size(x),C,size_out)*x
    end
  end
  return ex
end

# extract function and merge operator
function extract_merge_functions(t::Term)
  if is_sliced(t)
    if typeof(operator(t)) <: Compose
      op = operator(t).A[2]
    else
      op = Eye(size(operator(t),1)...)
    end
  else
    op = operator(t)
  end
  if is_eye(op)
    f = displacement(t) == 0 ? t.f : PrecomposeDiagonal(t.f, 1.0, displacement(t))
  elseif is_diagonal(op)
    f = PrecomposeDiagonal(t.f, diag(op), displacement(t))
  elseif is_AAc_diagonal(op)
    f = Precompose(t.f, op, diag_AAc(op), displacement(t))
  end
  f = t.lambda == 1. ? f : Postcompose(f, t.lambda) #for now I keep this
  #TODO change this
  return f
end

function extract_proximable(xAll::NTuple{N,Variable}, t::NTuple{M,Term}) where {N,M}
  fs = ()
  for x in xAll
    tx = () #terms containing x
    for ti in t
      if x in variables(ti)
        tx = (tx...,ti) #collect terms containing x
      end
    end
    if isempty(tx)
      fx = IndFree()
    elseif length(tx) == 1          #only one term per variable
      fx = extract_proximable(x,tx[1])
    else
      #multiple terms per variable
      #currently this happens only with GetIndex
      fxi,idxs = (),()
      for ti in tx
        fxi  = (fxi..., extract_merge_functions(ti))
        idx = typeof(operator(ti)) <: Compose ? operator(ti).A[1].idx : operator(ti).idx
        idxs = (idxs...,  idx   )
      end
      fx = SlicedSeparableSum(fxi,idxs)
    end
    fs = (fs...,fx)
  end
  if length(fs) > 1
    return SeparableSum(fs)  ##probably change constructor in Prox?
  else
    return fs[1]
  end
end

extract_proximable(xAll::Variable, t::Term) =  extract_merge_functions(t)
extract_proximable(xAll::NTuple{N,Variable}, t::Term) where {N} = extract_proximable(xAll,(t,))

ext: jl


text:
__precompile__()

module Projections

export project, project!

abstract type ConvexSet end
abstract type NonconvexSet end

"""
    project!(s, v, y)

Overwrite `v` by projection of point `y` to the box set `b`.
"""
function project!(); end

"""
    project(s, y)

Project point `y` to the box set `b`.
"""
function project(); end

include("box.jl")
include("ball.jl")
include("simplex.jl")
include("affine.jl")

end


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

## Powers
for f in (:unsafe_pow_inner,)
    @eval begin
        ($f)(::EmptyInterval{T}, ::S) where {T, S <: Real} =
            emptyinterval(promote_type(T, S))
        ($f)(::T, ::EmptyInterval{S}) where {T <: Real, S} =
            emptyinterval(promote_type(T, S))
        ($f)(::ValidInterval{T}, ::ValidInterval{S}) where {T, S} =
            emptyinterval(promote_type(T, S))
    end
end

# Integer power:

# overwrite new behaviour for small integer powers from
# https://github.com/JuliaLang/julia/pull/24240:

Base.literal_pow(::typeof(^), x::Interval{T}, ::Val{p}) where {T,p} =
    ^(x, p)

^(x::EmptyInterval, ::Integer) = x

function ^(a::Interval{T}, n::Integer) where T <: AbstractFloat
    n < 0 && return a == zero(a) ? emptyinterval(a) : inv(a^(-n))
    n == 0 && return one(a)
    n == 1 && return a
    a == zero(a) && return a
    n == 2 && return sqr(a)
    # we put an arbitrary limit to recursion at n == 10
    if isodd(n)
        n > 10 && return a^unsafe_atomic(T, n)
        return a * a^(n - 1)
    else
        n > 10 && return nonnegatives(T) ∩ a^unsafe_atomic(T, n)
        return nonnegatives(T) ∩ (a * a^(n - 1))
    end
end

# input a must be positive
sqr(a::AbstractFloat, mode::RoundingMode) = *(a, a, mode)

sqr(a::Interval{T}) where T = zero_in(a) ?
    unsafe_interval(zero(T), @round_up(sqr(mag(a)))) :
    @round(sqr(mig(a)), sqr(mag(a)))

# Interval power of an interval:
function ^(a::ValidInterval{T}, x::ValidInterval{T}) where T
    domain = nonnegatives(T)
    return unsafe_pow_inner(a ∩ domain, x)
end

function unsafe_pow_inner(a::Interval{T}, x::Interval{T}) where T
    a == zero(a) && return !isempty(x ∩ nonnegatives(T)) ? a : emptyinterval(T)
    return exp(x * log(a))
end

function sqrt(a::ValidInterval{T}) where T
    domain = nonnegatives(T)
    return sqrt_inner(a ∩ domain)
end
sqrt_inner(a::Interval) = @round(sqrt(a.lo), sqrt(a.hi))
sqrt_inner(x::EmptyInterval)= x

cbrt(a::Interval{BigFloat}) = @round( cbrt(a.lo), cbrt(a.hi) )
cbrt(x::EmptyInterval)= x

hypot(a::Interval{T}, b::Interval{T}) where T = sqrt(sqr(a) + sqr(b))


text:
using Plots

const serial = readdlm("out-serial.txt")
const darray = readdlm("out-darray.txt")
const mpi = readdlm("out-mpi.txt")
const elemental = readdlm("out-elemental.txt")

# column indices
const procs = 1
const threads = 2
const ddim = 3
const mintime = 4
const maxtime = 5

function filter(col,val,data)
  matching_rows = find(x -> x == val, data[:,col])
  result = similar(data, length(matching_rows), size(data,2))
  for (i,j) in enumerate(matching_rows)
    result[i,:] .= data[j,:]
  end
  return result
end

function plotcols(data, xcol, ycol; kwargs...)
  sorted = sortrows(data[:,[xcol,ycol]])
  x = sorted[:,1]
  y = sorted[:,2]
  plot!(log2.(x), y*1000; xticks=(log2.(x), string.(Int.(x))), kwargs...)
  return (x,y)
end

function plot_comparison(darray_filt, mpi_filt, serial_filt_time, filename, legpos)
  plotcols(filter(ddim, 1, darray_filt), procs, mintime, label="DArray, rows", marker=:square)
  plotcols(filter(ddim, 2, darray_filt), procs, mintime, label="DArray, cols", marker=:square)
  plotcols(filter(ddim, 1, mpi_filt), procs, mintime, label="MPI, rows", marker=:circle)
  (ideal_x, y) = plotcols(filter(ddim, 2, mpi_filt), procs, mintime, label="MPI, cols", marker=:circle)
  hline!([serial_filt_time*1000], label="Base.Array")
  plot!(log2.(ideal_x), serial_filt_time*1000 ./ ideal_x, label="Ideal scaling", linestyle=:dash, xlabel="Number of processes", ylabel="Time (ms)", legend=legpos)
  savefig(filename)
end

plotlyjs()

plotcols(filter(threads,1,elemental), procs, mintime, label="Elemental", marker=:hex)
plot_comparison(filter(threads,1,darray), filter(threads,1,mpi), minimum(filter(threads, 1, serial)[:,mintime]), "singlethread.svg", :top)

plot()
plot_comparison(filter(threads,32,darray), filter(threads,32,mpi), minimum(filter(threads, 32, serial)[:,mintime]), "multithread.svg", :topright)


text:
using SBMLCases
using Test

cases_dict = upload_cases()

@testset "upload test" begin
    @test length(cases_dict) == 1780
end

@testset "filter test" begin
    s1 = filter_cases(
        cases_dict,
        include_test_tags=["LocalParameters"],
        include_component_tags=["EventNoDelay", "CSymbolTime"],
        exclude_test_tags=["FastReaction", "Concentration"],
        exclude_component_tags=["AlgebraicRule"]
    )
    s2 = filter_cases(cases_dict)
    @test length(s1) == 454
    @test length(s2) == 1780
end

max_stars_repo_name: insysbio/sbml-heta-cases

max_issues_repo_name: insysbio/sbml-heta-cases

max_forks_repo_name: insysbio/sbml-heta-cases


text:
###################################
## Simple Variational Inference
## for 2D Gaussian

using PyPlot

function calc_KL(mu1, lambda1, mu2, lambda2)
    D = size(mu1, 1)
    px_lnqx = 0.5 * logdet(lambda2) - 0.5 * ((mu1 - mu2)' * lambda2 * (mu1 - mu2) + trace(lambda2 * inv(lambda1)))
    px_lnpx = 0.5 * logdet(lambda1) - 0.5 * D
    KL = - (px_lnqx - px_lnpx)
    return KL[1]
end

function plot_results(result, truth)
    N = size(result, 1)
    H = Int(ceil(sqrt(N)))
    W = Int(ceil(N / H))
    f, ax = subplots(H, W)
    for i in 1 : H
        for j in 1 : W
            n = (i - 1) * W + j
            if n <= N
                p = ax[i, j]
                p[:set_title](@sprintf("%d of %d", n, N))
                plot_gaussian(p, result[n][1], result[n][2], "b", "p(z)")
                plot_gaussian(p, truth[1], truth[2], "r", "q(z)")
                p[:set_xlim]([-1.5, 1.0])
                p[:set_ylim]([-0.6, 0.8])
            end
        end
    end
end

function plot_gaussian(p, Mu, Sigma, col, label)
    res = 100
    p[:plot](Mu[1], Mu[2], "x", color=col)

    F = eigfact(Sigma)
    vec = F[:vectors]
    val = F[:values]
    dw = 2*pi/res
    w = dw * (0 : res)

    c = 1.0
    a = sqrt(c * val[1])
    b = sqrt(c * val[2])
    P1 = a*cos(w)
    P2 = b*sin(w)
    P = Mu .+ vec'*vcat(P1', P2')
    p[:plot](P[1, :]', P[2, :]', "-", color=col, label=label)
end

function main()
    ## creat truth distribution
    D = 2 # dimension
    theta = 2.0*pi/12 # tilt
    A = reshape([cos(theta), -sin(theta),
                 sin(theta), cos(theta)],
                2, 2)
    mu = [0.0, 0.0]
    lambda = inv(A * inv(reshape([1,0,0,10], 2, 2)) * A')

    ## initialize
    mu_h = randn(D)
    lambda_h = zeros(D,D)

    ## main iteration
    max_iter = 10
    KL = Vector{Float64}()
    result = Vector()
    for i in 1 : max_iter
        ## update
        mu_h[1] = mu[1] - inv(lambda[1,1])*lambda[1,2] * (mu_h[2] - mu[2])
        lambda_h[1,1] = lambda[1,1]
        mu_h[2] = mu[2] - inv(lambda[2,2])*lambda[2,1] * (mu_h[1] - mu[1])
        lambda_h[2,2] = lambda[2,2]

        ## calculate KL divergeince
        push!(KL, calc_KL(mu_h, lambda_h, mu, lambda))
        #println(mu_h)
        #println(lambda_h)

        ## store the results
        push!(result, (deepcopy(mu_h), deepcopy(inv(lambda_h))))
    end

    ## visualize results
    plot_results(result, (mu, inv(lambda)))

    f, ax = subplots(1, 1)
    plot_gaussian(ax, mu, inv(lambda), "b", latexstring("\$p(z)\$"))
    plot_gaussian(ax, result[end][1], result[end][2], "r", latexstring("\$q(z)\$"))
    ax[:set_xlabel](latexstring("\$z_1\$"), fontsize=16)
    ax[:set_ylabel](latexstring("\$z_2\$"), fontsize=16)
    ax[:legend](fontsize=16)

    ## KL divergence
    f, ax = subplots(1)
    ax[:plot](1:max_iter, KL)
    ax[:set_ylabel]("KL divergence")
    ax[:set_xlabel]("iteration")
    show()
end

main()

max_stars_repo_name: sammy-suyama/mlblog

max_issues_repo_name: sammy-suyama/mlblog

max_forks_repo_name: sammy-suyama/mlblog


text:
using LinearAlgebra
using Gadfly
using GeometricalPredicates

# Arguments
# fd - funkcja odległości

struct Scaler
    scale::Any
    transx::Any
    transy::Any

    function Scaler(bbox::Array{Int,2})
        _scale = (abs(1 / (bbox[2] - bbox[1])))
        new(
            round(_scale, digits=4),
            round((1 + (0 - bbox[1] * _scale)), digits=4),
            round((1 + (0 - bbox[1, 2] * _scale)), digits=4),
        )
    end
end

function scaled_point(unscaled_point::Point2D, scaler::Scaler)
    return Point(
        round(((getx(unscaled_point) * scaler.scale) + scaler.transx), digits=4),
        round(((gety(unscaled_point) * scaler.scale) + scaler.transy), digits=4)
    )
end

function unscaled_point(scaled_point::Point2D, scaler::Scaler)
    return Point(
        round((getx(scaled_point) - scaler.transx) / scaler.scale, digits=4),
        round((gety(scaled_point) - scaler.transy) / scaler.scale, digits=4)
    )
end

function fh(x, y)::Real
    #  @info "Value is" s = 0.2
    return 0.2
end

function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where {T}
    m, n = length(vy), length(vx)
    #  @info "Length m = " m
    #  @info "Length n = " n
    gx = reshape(repeat(vx, inner = m, outer = 1), m, n)
    gy = reshape(repeat(vy, inner = 1, outer = n), m, n)
    #  @info "Grid x" gx
    #  @info "Grid y" gy
    return gx, gy
end

#bbox = [0 0; 1 1]

function extractbox(bbox, h0::Real, h1::Real)
    v1 = bbox[1, 1]:h0:bbox[2, 1]
    v2 = bbox[1, 2]:h1:bbox[2, 2]
    return v1, v2
end

function calculateh1(h0::Real)
    return h0 * (sqrt(3) / 2)
end

function shiftevenrows!(x, h0::Real)
    x[2:2:end, :] = x[2:2:end, :] .+ h0 / 2
end

function generate(fd, fh, bbox, h0)

    ttol = 0.1
    geps = 0.001 * h0
    Fscale = 1.2
    dptol = 0.001
    deps = sqrt(eps(Float64)) * h0
    h1 = calculateh1(h0)
    deltat = 0.2
    v1, v2 = extractbox(bbox, h0, h1)
    x, y = meshgrid(v1, v2)
    shiftevenrows!(x, h0)
    pb = [
        Point(row[1], row[2])
        for row in eachrow([x[:] y[:]]) if fd(row) < geps
    ]
    r0 = [1 ./ fh(getx(point), gety(point)) .^ 2 for point in pb]
    r0_max = maximum(r0 ./ maximum(r0))
    pold = Inf
    scaler = Scaler(bbox)
    pb = [
        scaled_point(point, scaler)
        for point in pb if (rand(Float64, size(pb)))[1] < r0_max
    ]
    while true
        tesselation = DelaunayTessellation()
        @info "Push to first tess"
        for p in pb
            push!(tesselation, p)
        end
        interior_points = Array{Point2D,1}()
        for triangle in tesselation
            center = centroid(Primitive(geta(triangle), getb(triangle), getc(triangle)))
            unscaled_p = unscaled_point(Point(getx(center), gety(center)), scaler)
            d = fd([getx(unscaled_p) gety(unscaled_p)])
            if d < -geps
                push!(interior_points, geta(triangle))
                push!(interior_points, getb(triangle))
                push!(interior_points, getc(triangle))
            end
        end
        interior_points = unique(interior_points)
        interior_tesselation = DelaunayTessellation()
        @info "Push to second tess without triangles"
        push!(interior_tesselation, interior_points)
        bars = Array{Point2D,1}()
        barvec = Array{Array{Float64,1},1}()
        points_to_fvces = Dict{Point2D,Array{Float64,1}}()
        @info "First iteration over edges started"

        #Point2D(1.4, 1.3464101615137753)
        #Point2D(1.4002786971987105, 1.3468928792218693)
        for edge in delaunayedges(interior_tesselation)
            b = unscaled_point(getb(edge), scaler)
            a = unscaled_point(geta(edge), scaler)
            # Obliczamy jako wartość fh() w połowie kazdego bara
            push!(
                bars,
                Point(
                    getx(a) + ((getx(b) - getx(a)) / 2),
                    gety(a) + ((gety(b) - gety(a)) / 2),
                ),
            )
            push!(
                barvec,
                [
                    getx(b) - getx(a)
                    gety(b) - gety(a)
                ],
            )
            push!(points_to_fvces, geta(edge) => [0, 0])
            push!(points_to_fvces, getb(edge) => [0, 0])
        end
        @info "First iteration over edges finished"
        L = [sqrt(sum(v_sum .^ 2)) for v_sum in barvec]
        iterator = 1
        hbars = [fh(getx(p), gety(p)) for p in bars]
        L0 = hbars * Fscale * sqrt(sum(L .^ 2) / sum(hbars .^ 2))
        F = maximum.(L0-L)
        println(typeof(F))
        @info "Force is " F
        println("END")
        println(size(pb))
        #Fvec=F./L*[1,1].*barvec
        Fvec = F ./ L .* barvec
        @info "Final iteration started"
        for edge in delaunayedges(interior_tesselation)
            prev_a = points_to_fvces[geta(edge)]
            prev_b = points_to_fvces[getb(edge)]
            push!(points_to_fvces, geta(edge) => prev_a + (-Fvec[iterator]))
            push!(points_to_fvces, getb(edge) => prev_b + (Fvec[iterator]))
            iterator = iterator + 1
        end
        @info "Final iteration finished"
        new_p = Array{Point2D,1}()
        d_points = Array{Array{Float64,1},1}()
        for (point, force) in points_to_fvces
            p = unscaled_point(point, scaler)
            np = [getx(p), gety(p)] + deltat * force
            push!(new_p, Point2D(np[1], np[2]))
            d = fd([getx(p), gety(p)])
            if d < -geps
                push!(d_points, force)
            end
        end
        @info "Mapping points"
        final_p = Array{Array{Float64,1},1}()
        for p in new_p
            x = getx(p)
            y = gety(p)
            d = fd([x, y])
            if d > 0
                dgradx = (fd([x + deps, y])-d)
                dgrady = (fd([x, y + deps])-d)
                rdgradx = dgradx/deps
                rdgrady = dgrady/deps
                res = [x, y] - [d*rdgradx, d*rdgrady]
                push!(final_p, res)
            else
                push!(final_p, [x, y])
            end
        end
        pb = map(p -> scaled_point(Point2D(p[1], p[2]), scaler), final_p)
        d_points = map(row -> sum(deltat * row .^ 2), d_points)
        @info "Checking break condition"
        ee = maximum(sqrt.(d_points)/h0)
        println(ee)
        if ee < dptol
            @info "Breaking out of loop"
            break
        end
   end
   tex = DelaunayTessellation()
   @info "Push to output tess"
   @info "Pb size is" size(pb)
   pb = unique!(pb)
   push!(tex, pb)
   x, y = getplotxy(delaunayedges(tex))
   x = map(p -> (p - scaler.transx) / scaler.scale, x)
   y = map(p -> (p - scaler.transy) / scaler.scale, y)
   @info "Plotting started"
   Gadfly.plot(x = x, y=y, Geom.path, Coord.cartesian(fixed=true))
end

function dc(p)
    #return -0.3 + abs(0.7 - sqrt(sum(p .^ 2)))
    return sqrt(sum(p .^ 2)) - 1
end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


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

import Test

Test.@testset "Package BATBase" begin
    include("my_mv_distribution_density.jl")
    include("basetypes/test_basetypes.jl")
end # testset


text:
using FastTransforms
using FastSphericalHarmonics
using ForwardDiff
using LinearAlgebra
using SpecialFunctions
using StaticArrays
using Test

bitsign(b::Bool) = b ? -1 : 1
bitsign(n::Integer) = bitsign(isodd(n))

function unit(val::T, ind::CartesianIndex{D}, size::NTuple{D,Int}) where {T,D}
    return setindex!(zeros(T, size), val, ind)
end
function unit(::Type{T}, ind::CartesianIndex{D},
              size::NTuple{D,Int}) where {T,D}
    return unit(one(T), ind, size)
end
function unit(ind::CartesianIndex{D}, size::NTuple{D,Int}) where {D}
    return unit(true, ind, size)
end

chop(x) = abs2(x) < 100eps(x) ? zero(x) : x
chop(x::Complex) = Complex(chop(real(x)), chop(imag(x)))

chop2(x) = chop.(x)

# [Generalized binomial coefficient](https://en.wikipedia.org/wiki/Binomial_coefficient#Generalization_and_connection_to_the_binomial_series)
function Base.binomial(α::Number, k::Integer)
    k == 0 && return one(α)
    return prod((α - i) / (k - i) for i in 0:(k - 1))
end

# [Jacobi
# Polynomials](https://en.wikipedia.org/wiki/Jacobi_polynomials)
function JacobiP(α, β, n, z)
    return gamma(α + n + 1) / (factorial(n) * gamma(α + β + n + 1)) *
           sum(binomial(n, m) * gamma(α + β + n + m + 1) / gamma(α + m + 1) *
               ((z - 1) / 2)^m for m in 0:n)
end

# [Associated Legendre
# Polynomials](https://en.wikipedia.org/wiki/Associated_Legendre_polynomials)
function LegendreP(l, m, x)
    return bitsign(m) *
           2^l *
           sqrt(1 - x^2)^m *
           sum(factorial(k) / factorial(k - m) *
               x^(k - m) *
               binomial(l, k) *
               binomial((l + k - 1) / 2, l) for k in m:l)
end

# [Spherical
# harmonics](https://mikaelslevinsky.github.io/FastTransforms/transforms.html)
# (section "sph2fourier")
function Ylm(l, m, θ, ϕ)
    return bitsign(abs(m)) *
           sqrt((l + 1 / 2) * factorial(l - abs(m)) / factorial(l + abs(m))) *
           LegendreP(l, abs(m), cos(θ)) *
           sqrt((2 - (m == 0)) / 2π) *
           (m ≥ 0 ? cos(abs(m) * ϕ) : sin(abs(m) * ϕ))
end

# [Spin-weighted spherical
# harmonics](https://mikaelslevinsky.github.io/FastTransforms/transforms.html)
# (section "spinsph2fourier")
function sYlm(s, l, m, θ, ϕ)
    l0 = max(abs(m), abs(s))
    l1 = min(abs(m), abs(s))
    return cis(m * ϕ) / sqrt(2π) *
           sqrt((l + 1 / 2) * factorial(l + l0) * factorial(l - l0) /
                (factorial(l + l1) * factorial(l - l1))) *
           sin(θ / 2)^abs(m + s) *
           cos(θ / 2)^abs(m - s) *
           JacobiP(abs(m + s), abs(m - s), l - l0, cos(θ))
end
sYlm(::Type{<:Complex}, s, l, m, θ, ϕ) = sYlm(s, l, m, θ, ϕ)
function sYlm(::Type{<:Real}, s, l, m, θ, ϕ)
    @assert s == 0
    if m == 0
        return real(sYlm(s, l, abs(m), θ, ϕ))
    elseif m > 0
        return sqrt(2) * real(sYlm(s, l, abs(m), θ, ϕ))
    else
        return sqrt(2) * imag(sYlm(s, l, abs(m), θ, ϕ))
    end
end

c2a(c) = SVector(real(c), imag(c))
a2c(a) = Complex(a[1], a[2])

∂θYlm(l, m, θ, ϕ) = a2c(ForwardDiff.derivative(θ -> c2a(Ylm(l, m, θ, ϕ)), θ))
∂ϕYlm(l, m, θ, ϕ) = a2c(ForwardDiff.derivative(ϕ -> c2a(Ylm(l, m, θ, ϕ)), ϕ))
function ∂θsYlm(s, l, m, θ, ϕ)
    return a2c(ForwardDiff.derivative(θ -> c2a(sYlm(s, l, m, θ, ϕ)), θ))
end
function ∂ϕsYlm(s, l, m, θ, ϕ)
    return a2c(ForwardDiff.derivative(ϕ -> c2a(sYlm(s, l, m, θ, ϕ)), ϕ))
end
function ∂θsYlm(::Type{T}, s, l, m, θ, ϕ) where {T}
    return a2c(ForwardDiff.derivative(θ -> c2a(sYlm(T, s, l, m, θ, ϕ)), θ))
end
function ∂ϕsYlm(::Type{T}, s, l, m, θ, ϕ) where {T}
    return a2c(ForwardDiff.derivative(ϕ -> c2a(sYlm(T, s, l, m, θ, ϕ)), ϕ))
end

ðYlm(l, m, θ, ϕ) = -(∂θYlm(l, m, θ, ϕ) + im / sin(θ) * ∂ϕYlm(l, m, θ, ϕ))
ð̄Ylm(l, m, θ, ϕ) = -(∂θYlm(l, m, θ, ϕ) - im / sin(θ) * ∂ϕYlm(l, m, θ, ϕ))
function ðsYlm(s, l, m, θ, ϕ)
    return (-(∂θsYlm(s, l, m, θ, ϕ) + im / sin(θ) * ∂ϕsYlm(s, l, m, θ, ϕ)) +
            s * cot(θ) * sYlm(s, l, m, θ, ϕ))
end
function ð̄sYlm(s, l, m, θ, ϕ)
    return (-(∂θsYlm(s, l, m, θ, ϕ) - im / sin(θ) * ∂ϕsYlm(s, l, m, θ, ϕ)) -
            s * cot(θ) * sYlm(s, l, m, θ, ϕ))
end
function ðsYlm(::Type{T}, s, l, m, θ, ϕ) where {T}
    return (-(∂θsYlm(T, s, l, m, θ, ϕ) + im / sin(θ) * ∂ϕsYlm(T, s, l, m, θ, ϕ)) +
            s * cot(θ) * sYlm(T, s, l, m, θ, ϕ))
end
function ð̄sYlm(::Type{T}, s, l, m, θ, ϕ) where {T}
    return (-(∂θsYlm(T, s, l, m, θ, ϕ) - im / sin(θ) * ∂ϕsYlm(T, s, l, m, θ, ϕ)) -
            s * cot(θ) * sYlm(T, s, l, m, θ, ϕ))
end

# Real spherical harmonics
Y_0_0(θ, ϕ) = sqrt(1 / 4π)

Y_1_0(θ, ϕ) = sqrt(3 / 4π) * cos(θ)
Y_1p1(θ, ϕ) = sqrt(3 / 4π) * sin(θ) * cos(ϕ)
Y_1m1(θ, ϕ) = sqrt(3 / 4π) * sin(θ) * sin(ϕ)

# Complex spin-weighted spherical harmonics

# sYlm[s_, l_, m_, q_, f_] =
#  Module[{l0 = Max[Abs[m], Abs[s]], l1 = Min[Abs[m], Abs[s]]},
#      (Exp[I*m*f]/Sqrt[2*Pi])*
#    Sqrt[(l + 1/2)*(l + l0)!*((l - l0)!/((l + l1)!*(l - l1)!))]*
#        Sin[q/2]^Abs[m + s]*Cos[q/2]^Abs[m - s]*
#    JacobiP[l - l0, Abs[m + s], Abs[m - s],
#          Cos[q]]]

Y_s0_0_0(θ, ϕ) = 1 / (2 * sqrt(π))

Y_s0_1_0(θ, ϕ) = 1 / 2 * sqrt(3 / π) * cos(θ)
Y_s0_1p1(θ, ϕ) = cis(ϕ) * sqrt(3 / (2 * π)) * cos(θ / 2) * sin(θ / 2)
Y_s0_1p1(θ, ϕ) = cis(-f) * sqrt(3 / (2 * π)) * cos(θ / 2) * sin(θ / 2)

Y_s1_1_0(θ, ϕ) = sqrt(3 / (2 * π)) * cos(θ / 2) * sin(θ / 2)
Y_s1_1p1(θ, ϕ) = 1 / 2 * cis(f) * sqrt(3 / π) * sin(θ / 2)^2
Y_s1_1m1(θ, ϕ) = 1 / 2 * cis(-f) * sqrt(3 / π) * cos(θ / 2)^2

Y_s1_2_0(θ, ϕ) = sqrt(15 / (2 * π)) * cos(θ / 2) * cos(θ) * sin(θ / 2)
Y_s1_2p1(θ, ϕ) = 1 / 4 * cis(f) * sqrt(5 / π) * (2 + 4 * cos(θ)) * sin(θ / 2)^2
Y_s1_2p2(θ, ϕ) = cis(2 * f) * sqrt(5 / π) * cos(θ / 2) * sin(θ / 2)^3

# # Real gradient spherical harmonics
# # [∂θ, 1/sin(θ) ∂ϕ]
# grad_Y_1_0_θ(θ, ϕ) = sqrt(3 / 4π) * sin(θ)
# grad_Y_1_0_ϕ(θ, ϕ) = 0.0
# grad_Y_1p1_θ(θ, ϕ) = sqrt(3 / 4π) * cos(θ) * cos(ϕ)
# grad_Y_1p1_ϕ(θ, ϕ) = sqrt(3 / 4π) * cos(ϕ)
# grad_Y_1m1_θ(θ, ϕ) = sqrt(3 / 4π) * cos(θ) * sin(ϕ)
# grad_Y_1m1_ϕ(θ, ϕ) = sqrt(3 / 4π) * sin(ϕ)
#
# grad_Y_2_0_θ(θ, ϕ) = sqrt(45 / 4π) * cos(θ) * sin(θ)
# grad_Y_2_0_ϕ(θ, ϕ) = 0.0
# grad_Y_2p1_θ(θ, ϕ) = sqrt(15 / 4π) * cos(2θ) * cos(ϕ)
# grad_Y_2p1_ϕ(θ, ϕ) = sqrt(15 / 4π) * cos(θ) * cos(ϕ)
#
# grad_Y_3_0_θ(θ, ϕ) = sqrt(63 / 256π) * (5 * sin(3θ) + sin(θ))
# grad_Y_3_0_ϕ(θ, ϕ) = 0.0
# grad_Y_3p1_θ(θ, ϕ) = sqrt(42 / 256π) * (15 * cos(2θ) - 7) * cos(θ) * cos(ϕ)
# grad_Y_3p1_ϕ(θ, ϕ) = sqrt(42 / 256π) * (5 * cos(2θ) + 3) * cos(ϕ)
#
# # Real curl spherical harmonics
# # [1/sin(θ) ∂ϕ, -∂θ]
# curl_Y_1_0_θ(θ, ϕ) = grad_Y_1_0_ϕ(θ, ϕ)
# curl_Y_1_0_ϕ(θ, ϕ) = -grad_Y_1_0_θ(θ, ϕ)
# curl_Y_1p1_θ(θ, ϕ) = grad_Y_1p1_ϕ(θ, ϕ)
# curl_Y_1p1_ϕ(θ, ϕ) = -grad_Y_1p1_θ(θ, ϕ)
# curl_Y_1m1_θ(θ, ϕ) = grad_Y_1m1_ϕ(θ, ϕ)
# curl_Y_1m1_ϕ(θ, ϕ) = -grad_Y_1m1_θ(θ, ϕ)

# Integrate over the sphere
function integrate(u::Array{T1,2}, v::Array{T2,2}) where {T1,T2}
    @assert size(u) == size(v)
    N, M = size(u)
    @assert M > 0 && N > 0
    @assert M == 2 * N - 1
    Θ, Φ = sph_points(N)
    s = zero(conj(one(T1)) * one(T2) * one(eltype(Θ)))
    for j in 1:M, i in 1:N
        s += conj(u[i, j]) * v[i, j] * sin(Θ[i])
    end
    return 2π^2 / (N * M) * s
end

function integrate(uθ::Array{T1,2}, uϕ::Array{T1,2}, vθ::Array{T2,2},
                   vϕ::Array{T2,2}) where {T1,T2}
    @assert size(uθ) == size(uϕ) == size(vθ) == size(vϕ)
    N, M = size(uθ)
    @assert M > 0 && N > 0
    @assert M == 2 * N - 1
    Θ, Φ = sph_points(N)
    s = zero(conj(one(T1)) * one(T2) * one(eltype(Θ)))
    for j in 1:M, i in 1:N
        s += (conj(uθ[i, j]) * vθ[i, j] + conj(uϕ[i, j]) * vϕ[i, j]) * sin(Θ[i])
    end
    return 2π^2 / (N * M) * s
end

################################################################################

include("test_fasttransforms.jl")
include("test_points.jl")
include("test_scalar.jl")
# include("test_vector.jl")
include("test_spin.jl")


text:
"""
    EGPbeta(σ, ξ, δ)

*EGPbeta* corresponds to the third extended GP model of Naveau et al. (2016).

It is a three parameters family: δ is a threshold tuning parameter, σ is a scale parameter, and ξ is a shape parameter.

```julia
EGPbeta(σ, ξ, δ)   # EGP of Naveau et al. (2016) (type 3) with parameters σ, ξ, δ.

params(d)           # Get the parameters, i.e. (σ, ξ, δ)
```

Reference :

* Naveau, P., Huser, R., Ribereau, P., and Hannart, A. (2016), Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection, Water Resour. Res., 52, 2753– 2769, doi:10.1002/2015WR018552.
"""
struct EGPbeta{T<:Real} <: ContinuousUnivariateDistribution
    σ::T    # scale parameter
    ξ::T    # shape parameter
    δ::T    # threshold tuning parameter

    function EGPbeta{T}(σ::T, ξ::T, δ::T) where {T <: Real}
        new{T}(σ, ξ, δ)
    end
end

function EGPbeta(σ::T, ξ::T, δ::T; check_args=true) where {T <: Real}
    check_args && @check_args(EGPbeta, σ > zero(σ) && δ > zero(δ))
    return EGPbeta{T}(σ, ξ, δ)
end

EGPbeta(σ::Real, ξ::Real, δ::Real) = EGPbeta(promote(σ, ξ, δ)...)
EGPbeta(σ::Integer, ξ::Integer, δ::Integer) = EGPbeta(float(σ), float(ξ), float(δ))

minimum(d::EGPbeta) = 0.0
maximum(d::EGPbeta) = Inf * (d.ξ >= 0) - ( d.σ / d.ξ ) * (d.ξ < 0)
insupport(d::EGPbeta, x::Real) = minimum(d) <= x <= maximum(d)

#### Parameters

scale(d::EGPbeta) = d.σ
shape(d::EGPbeta) = d.ξ    # noms à vérifier...
thresh(d::EGPbeta) = d.δ

params(d::EGPbeta) = (d.σ, d.ξ, d.δ)
partype(::EGPbeta{T}) where {T} = T

#### Evaluation

function logpdf(d::EGPbeta{T}, x::Real) where T<:Real
    μ = 0
    (σ, ξ, δ) = params(d)

    pd = GeneralizedPareto(μ, 1, ξ)

    #g(v::Real, δ) = δ*(1-v)^(δ-1) * pdf(Beta(1/δ, 2), (1 - v)^δ)
    lg(v::Real, δ) = log(δ) + (δ-1)*log(1-v) + Distributions.logpdf(Beta(1/δ, 2), (1 - v)^δ)

    lf = -log(σ) + Distributions.logpdf.(pd, x/σ) + lg(Distributions.cdf(pd, x/σ), δ)

    return lf
end

pdf(d::EGPbeta, x::Real) = exp(logpdf(d, x))

function logcdf(d::EGPbeta{T}, x::Real) where T<:Real
    μ = 0
    (σ, ξ, δ) = params(d)

    pd = GeneralizedPareto(μ, 1, ξ)

    G(v::Real, δ::Real) = 1 - Distributions.cdf(Beta(1/δ, 2), (1 - v)^δ)

    lf = log(G(Distributions.cdf(pd, x/σ), δ))

    return lf
end

cdf(d::EGPbeta, x::Real) = exp(logcdf(d, x))

function quantile(d::EGPbeta{T}, p::Real) where T<:Real
    @assert zero(p)<p<one(p) "the quantile level should be between 0 and 1."

    (σ, ξ, δ) = params(d)

    invG(v::Real, δ::Real) = 1 - Distributions.quantile(Beta(1/δ, 2), (1 - v))^(1/δ)

    x = (σ/ξ)*(((1 - invG(p,δ))^(-ξ)) - 1)

    return x
end

#### Sampling

function rand(rng::AbstractRNG, d::EGPbeta)
    # Generate a Float64 random number uniformly in (0,1].
    u = 1 - rand(rng)

    invG(v::Real, δ::Real) = 1 - Distributions.quantile(Beta(1/δ, 2), (1 - v))^(1/δ)

    return (d.σ/d.ξ)*(((1 - invG(u,d.δ))^(-d.ξ)) - 1)
end

sampler(d::EGPbeta) = d


text:
using LinearAlgebra
import Base.FastMath.@fastmath

abstract type Side end
struct BothSides <: Side end
struct ToTheLeft <: Side end
struct ToTheRight <: Side end

function determine_input_activations(model::GebModel, organism::GebOrganism)
    all_neighbors = _get_neighbors(BothSides(), model, organism)

    left_neighbors = nothing
    right_neighbors = nothing

    loaded_left_neighbors = false
    loaded_right_neighbors = false

    activations = Float64[]

    distances = Dict()

    for input in organism.network.inputs
        string = input.string

        if length(string) == 0
            push!(activations, 0.0)
            continue
        end

        if string[1] == '0'
            if loaded_left_neighbors == false
                left_neighbors = _get_neighbors(ToTheLeft(), model, organism, all_neighbors)
                loaded_left_neighbors == true
            end
            relevant_neighbors = left_neighbors
        else
            if loaded_right_neighbors == false
                right_neighbors = _get_neighbors(ToTheRight(), model, organism, all_neighbors)
                loaded_right_neighbors == true
            end
            relevant_neighbors = right_neighbors
        end

        activation_sum = 0.0
        for neighbor in relevant_neighbors
            neighbor_sum = 0.0

            for output in neighbor.network.outputs
                if _is_match(string[2:end], output.string)
                    neighbor_sum += output.io_value
                end
            end

            if 0.0 < neighbor_sum
                if haskey(distances, neighbor) == false
                    distances[neighbor] = @fastmath norm(organism.coordinates .- neighbor.coordinates)
                end

                activation_sum += DISTANCE_SCALE_FACTOR * (neighbor_sum / distances[neighbor])
            end
        end

        push!(activations, activation_sum)
    end

    return activations
end

function _get_neighbors(side::Side, model::GebModel, organism::GebOrganism, all_neighbors)
    [neighbor for neighbor in all_neighbors if _is_on_side(side, model, neighbor, organism)]
end

function _get_neighbors(side::Side, model::GebModel, organism::GebOrganism)
    _get_neighbors(side, model, organism, _get_neighbors(BothSides(), model, organism))
end

function _get_neighbors(::BothSides, model::GebModel, organism::GebOrganism)
    neighbors = []

    x_center, y_center = _get_grid_coordinates(organism.coordinates)

    for x in (x_center - INFLUENCE_RADIUS):(x_center + INFLUENCE_RADIUS)
        for y in (y_center - INFLUENCE_RADIUS):(y_center + INFLUENCE_RADIUS)
            x, y = mod_1(x, model.size), mod_1(y, model.size)

            if model.grid[x, y] !== nothing && model.grid[x, y] != organism
                push!(neighbors, model.grid[x, y])
            end
        end
    end

    return neighbors
end

function _is_on_side(side::Side, model::GebModel, target::GebOrganism, source::GebOrganism)
    _is_on_side(side, model, target.coordinates, source.coordinates, source.direction)
end
function _is_on_side(side::Side, model::GebModel, target_coordinates, source_coordinates, source_direction)
    Δx = @fastmath (target_coordinates[1] - source_coordinates[1])
    Δy = @fastmath -(target_coordinates[2] - source_coordinates[2]) # negative as direction of y-axis is different in array!

    # test if wrap-over sides is nearer
    if abs(Δx) >= model.size / 2
        Δx = @fastmath (Δx - sign(Δx)*model.size)
    end
    if abs(Δy) >= model.size / 2
        Δy = @fastmath (Δy - sign(Δy)*model.size)
    end

    if Δx == 0
        target_direction = (0 <= Δy) ? 90 : 270
    elseif Δx < 0
        target_direction = @fastmath (180 + atand(Δy / Δx))
    elseif 0 < Δx
        target_direction = @fastmath (atand(Δy / Δx))
    end

    return _is_on_side(side, target_direction, source_direction)
end
function _is_on_side(::ToTheLeft, target_direction::Number, source_direction::Number)
    mod(target_direction - source_direction, 360) <= 180
end
function _is_on_side(::ToTheRight, target_direction::Number, source_direction::Number)
    mod(target_direction - source_direction, 360) > 180
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
"""
$(TYPEDEF)

A `VertexOrder` represents the order in which the vertices of a [`ConvexHull`](@ref)
are stored.
"""
abstract type VertexOrder end
orientation_comparator(o::VertexOrder) = orientation_comparator(typeof(o))
edge_normal_sign_operator(o::VertexOrder) = orientation_comparator(typeof(o))

"""
$(TYPEDEF)

Counterclockwise vertex order.
"""
struct CCW <: VertexOrder end
orientation_comparator(::Type{CCW}) = >
edge_normal_sign_operator(::Type{CCW}) = +

"""
$(TYPEDEF)

Clockwise vertex order.
"""
struct CW <: VertexOrder end
orientation_comparator(::Type{CW}) = <
edge_normal_sign_operator(::Type{CW}) = -


text:
# ==================================
# CD Player
#
# system type: continuous LTI system
# state dimension: 120
# input dimension: 2
# ==================================
using ReachabilityBenchmarks, MathematicalSystems, LazySets, MAT

function cdplayer_model()
    file = matopen(@relpath "cdplayer.mat")

    # system matrix
    A = read(file, "A")

    # input matrix
    B = read(file, "B")

    # state domain
    X = Universe(120)

    # input domain
    U = BallInf([0.0, 0.0], 1.0)

    # continuous LTI system
    S = @system(x' = Ax + Bu, x ∈ X, u ∈ U)

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
@testset "recv_close.jl" begin
    @testset begin receiver = TCPReceiver(4000)
        @test window_size(receiver) == 4000
        @test ackno(receiver) === nothing
        isn = rand(UInt32(0):typemax(UInt32))
        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn+1)
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0
        segment_arrives(receiver; with_fin=true, seqno=isn+1, result=SegmentArrives_Result_OK)
        @test read!(stream_out(receiver)) == ""
        @test assembled_bytes(receiver) == 0
        @test eof(stream_out(receiver))
    end

    @testset begin receiver = TCPReceiver(4000)
        @test window_size(receiver) == 4000
        @test ackno(receiver) === nothing
        isn = rand(UInt32(0):typemax(UInt32))
        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn+1)
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0
        segment_arrives(receiver; data="a", with_fin=true, seqno=isn+1, result=SegmentArrives_Result_OK)
        @test read!(stream_out(receiver)) == "a"
        @test assembled_bytes(receiver) == 1
        @test eof(stream_out(receiver))
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
s = 0
for i = 1:999
    if i%3 == 0 || i%5 == 0
        global s+=i
    end
end
println(s)


text:
using Documenter, LadeGCM

makedocs()

deploydocs(
    deps   = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    repo   = "github.com/Libbum/LadeGCM.jl.git",
    julia  = "0.6"
)

max_stars_repo_name: Libbum/LadeGCM

max_issues_repo_name: Libbum/LadeGCM

max_forks_repo_name: Libbum/LadeGCM


text:
function hgf_url_template(repo_id, revision, filename)
    global ENDPOINT
    return "$(ENDPOINT[])/$(repo_id)/resolve/$(revision)/$(filename)"
end

abstract type REPO end
struct DATASET_REPO <: REPO end
struct SPACE_REPO <: REPO end

const REPO_TYPE = Union{REPO, Nothing}

is_valid_repo(::REPO_TYPE) = true
is_valid_repo(_) = false

repo_string(::DATASET_REPO) = "datasets"
repo_string(::SPACE_REPO) = "spaces"

build_repo(repo_type::REPO, repo_id) = joinpath(repo_string(repo_type), repo_id)
build_repo(::Nothing, repo_id) = repo_id

"""
    HuggingFaceURL(repo_id, [subfolder], filename,
                   repo_type::REPO_TYPE = nothing,
                   revision::AbstractString = "main")

Construct the real url with the inputs.
"""
struct HuggingFaceURL
    repo_id   :: String
    filename  :: String
    repo_type :: REPO_TYPE
    revision  :: String

    function HuggingFaceURL(repo_id, filename,
                            repo_type::REPO_TYPE = nothing,
                            revision::AbstractString = "main")
        is_valid_repo(repo_type) || error("Invalid repo type: $repo_type")
        return new(repo_id, filename, repo_type, revision)
    end
end

HuggingFaceURL(
    repo_id, subfolder, filename,
    repo_type::REPO_TYPE = nothing,
    revision::AbstractString = "main"
) =
    HuggingFaceURL(repo_id, joinpath(subfoler, filename), repo_type, revision)

repo(hgfurl::HuggingFaceURL) = build_repo(hgfurl.repo_type, hgfurl.repo_id)
id(hgfurl::HuggingFaceURL) = joinpath(repo(hgfurl), hgfurl.filename)

Base.string(hgfurl::HuggingFaceURL) = hgf_url_template(repo(hgfurl), hgfurl.revision, hgfurl.filename)
Base.show(io::IO, hgfurl::HuggingFaceURL) = Base.print(io, "HuggingFaceURL(", string(hgfurl), ')')

function get_etag(
    url;
    timeout :: Real = 10,
    headers :: Union{AbstractVector, AbstractDict} = Pair{String,String}[],
)
    r = request(url; method="HEAD", timeout, headers)
    if r.status < 400
        for (key, val) in r.headers
            if key == "etag" || key == "x-linked-etag"
                return val
            end
        end
    end
    error("Distant resource does not have an ETag, we won't be able to reliably ensure reproducibility.")
end

url_to_filename(hgfurl::HuggingFaceURL, etag=nothing) = url_to_filename(string(hgfurl), etag)
function url_to_filename(url, etag=nothing)
    buf = IOBuffer(maxsize=150)
    bytes2hex(buf, sha256(url))
    isnothing(etag) || (write(buf, '.'); bytes2hex(buf, sha256(etag)))
    # endswith(url, ".h5") && write(buf, ".h5")
    # endswith(url, ".py") && wrtie(buf, ".py")
    return String(take!(buf))
end

hgf_download(
    hgfurl :: HuggingFaceURL, dest :: AbstractString;
    headers :: Union{AbstractVector, AbstractDict} = Pair{String,String}[],
    verbose :: Bool = false, io :: IO = Pkg.DEFAULT_IO[]
) =
    hgf_download(string(hgfurl), dest; name = id(hgfurl), headers, verbose, io)

function hgf_download(
    url :: AbstractString, dest :: AbstractString;
    name :: AbstractString = "",
    headers :: Union{AbstractVector, AbstractDict} = Pair{String,String}[],
    verbose :: Bool = false,
    io :: IO = Pkg.DEFAULT_IO[],
)
    do_fancy = verbose && Pkg.can_fancyprint(io)
    progress =  if do_fancy
        bar = MiniProgressBar(header="Downloading $name", color=Base.info_color())
        start_progress(io, bar)
        (total, now) -> begin
            bar.max = total
            bar.current = now
            show_progress(io, bar)
        end
    else
        (total, now) -> nothing
    end
    try
        Downloads.download(url, dest; headers, progress)
    finally
        do_fancy && end_progress(io, bar)
    end
end

"""
    cached_download(
        hgfurl :: HuggingFaceURL;
        local_files_only :: Bool = false,
        auth_token :: Union{AbstractString, Nothing} = nothing,
    )

Find the local cache of given url or do downloading. If `local_files_only` is set, it will try to
 find the file from cache, and error out when not found. For downloading from private repo,
 `auth_token` need to be set, or do `HuggingFaceApi.login()` beforehand.

See also: [`HuggingFaceURL`](@ref), [`login`](@ref)
"""
function cached_download(
    hgfurl :: HuggingFaceURL;
    local_files_only :: Bool = false,
    auth_token :: Union{AbstractString, Nothing} = get_token(),
)
    # add auth token if provided
    headers = auth_header(auth_token)

    url = string(hgfurl)
    etag = local_files_only ? nothing : get_etag(url; headers)
    name = url_to_filename(url, etag)

    # local file exists
    hash = @my_artifact :hash name
    !isnothing(hash) && my_artifact_exists(hash) && return my_artifact_path(hash)

    # no local file
    local_files_only && error("no cached file found and `local_files_only` is set to true.")

    # download and bind
    hash = @my_artifact :download name hgfurl hgf_download headers verbose=true

    # for future local_files_only, bind name without etag
    @my_artifact :bind split(name, '.')[1] hash force=true

    return my_artifact_path(hash)
end

"""
    remove_cache(hgfurl::HuggingFaceURL; now=false)

Remove files link to the given url. If `now` is set to `true`, cache file will be deleted immediately,
 otherwise waiting `OhMyArtifacts` to do the garbage collection.
"""
function remove_cache(hgfurl::HuggingFaceURL; now=false)
    name = url_to_filename(hgfurl)
    artifacts = OhMyArtifacts.load_my_artifacts_toml(hgf_artifacts())

    for (key, val) in artifacts
        if startswith(key, name)
            @my_artifact :unbind key
        end
    end

    if now
        OhMyArtifacts.find_orphanages(; collect_delay = Dates.Hour(0))
    end
end

"""
    hf_hub_download(
        repo_id :: AbstractString,
        filename :: AbstractString;
        repo_type :: REPO_TYPE = nothing,
        revision :: AbstractString = "main",
        auth_token :: Union{AbstractString, Nothing} = get_token(),
        local_files_only :: Bool  = false,
    )

Construct [`HuggingFaceURL`](@ref) and do [`cached_download`](@ref).
"""
function hf_hub_download(
    repo_id :: AbstractString,
    filename :: AbstractString;
    repo_type :: REPO_TYPE = nothing,
    revision :: AbstractString = "main",
    auth_token :: Union{AbstractString, Nothing} = get_token(),
    local_files_only :: Bool  = false,
)
    hgfurl = HuggingFaceURL(repo_id, filename, repo_type, revision)
    return cached_download(hgfurl; auth_token, local_files_only)
end


text:
function diffeq!(du,u,p,t)
    Complex = p[C.B55tot] - u[V.PP2AB55]
    Wee1p = 1 - u[V.Wee1] - u[V.Wee1pp]
    Cdc25p = 1 - u[V.Cdc25] - u[V.Cdc25pp]
    Vwee = (p[C.kweeS]*(1-u[V.Wee1]) + p[C.kweeF]*u[V.Wee1])
    V25 = p[C.k25S]*(1-u[V.Cdc25pp]) + p[C.k25F]*u[V.Cdc25pp]
    VGwl = p[C.kGwENSA]*u[V.Gwlp]

    du[V.Subp] = p[C.kcBc1Sub]*u[V.CycBCdk1]/((1 + (p[C.InhCDK]/p[C.Kd])))*(p[C.SubT]-u[V.Subp]) - p[C.kB55Sub]*u[V.PP2AB55]*u[V.Subp]
    du[V.CycBCdk1] = V25*(p[C.CycBCdk1T] - u[V.CycBCdk1]) - Vwee*u[V.CycBCdk1]
    du[V.PP1] = (p[C.kapp1] + p[C.kapp1a]*u[V.PP1])*(p[C.PP1T] - u[V.PP1]) - (p[C.kipp1] + p[C.kipp1C]*u[V.CycBCdk1]/((1 + (p[C.InhCDK]/p[C.Kd]))))*u[V.PP1]
    du[V.pENSAt] = VGwl*(p[C.ENSAtot] - u[V.pENSAt]) - p[C.kcatB55]*Complex
    du[V.Gwlp] = (p[C.kcBc1G]*u[V.CycBCdk1]/((1 + (p[C.InhCDK]/p[C.Kd]))) + p[C.kcAc2G]*p[C.CycACdk2T])*(p[C.Gwtot] - u[V.Gwlp]) - (p[C.kB55G]*u[V.PP2AB55] + p[C.kppxGwl] + p[C.kPP1Gw]*u[V.PP1])*u[V.Gwlp]
    du[V.PP2AB55] = p[C.kdis]*Complex + p[C.kcatB55]*Complex - p[C.kass]*u[V.PP2AB55]*(u[V.pENSAt] - Complex)
    du[V.Wee1] = (p[C.kppxY15] + p[C.kB55W1]*u[V.PP2AB55])*Wee1p - (p[C.kcBc1W1]*u[V.CycBCdk1]/((1 + (p[C.InhCDK]/p[C.Kd]))) + p[C.kcAc2W1]*p[C.CycACdk2T])*u[V.Wee1]
    du[V.Wee1pp] = (p[C.kcBc1W1]*u[V.CycBCdk1]/((1 + (p[C.InhCDK]/p[C.Kd]))) + p[C.kcAc2W1]*p[C.CycACdk2T])*Wee1p - (p[C.kppxY15] + p[C.kB55W1]*u[V.PP2AB55])*u[V.Wee1pp]
    du[V.Cdc25] = (p[C.kppxY15] + p[C.kB5525]*u[V.PP2AB55])*Cdc25p - (p[C.kcBc125]*u[V.CycBCdk1]/((1 + (p[C.InhCDK]/p[C.Kd]))) + p[C.kcAc225]*p[C.CycACdk2T])*u[V.Cdc25]
    du[V.Cdc25pp] = (p[C.kcBc125]*u[V.CycBCdk1]/((1 + (p[C.InhCDK]/p[C.Kd]))) + p[C.kcAc225]*p[C.CycACdk2T])*(Cdc25p) - (p[C.kppxY15] + p[C.kB5525]*u[V.PP2AB55])*u[V.Cdc25pp]

end

function param_values()::Vector{Float64}
    p::Vector{Float64} = zeros(C.NUM)

    p[C.InhCDK] = 2.0
    p[C.CycBCdk1T] = 8.1808
    p[C.CycACdk2T] = 1.0000
    p[C.PP1T] = 1.0000
    p[C.kapp1] = 0.0115
    p[C.kapp1a] = 0.7054
    p[C.kipp1] = 0.0018
    p[C.kipp1C] = 0.7549
    p[C.kPP1Gw] = 18.4724
    p[C.ENSAtot] = 1.0000
    p[C.B55tot] = 0.2500
    p[C.SubT] = 1.0000
    p[C.kass] = 617.2807
    p[C.kdis] = 0.0088
    p[C.kcatB55] = 1.0338
    p[C.kGwENSA] = 20.8811
    p[C.kppxGwl] = 0.1560
    p[C.kcBc1Sub] = 0.0080
    p[C.kcBc1G] = 0.2393
    p[C.Gwtot] = 1.0000
    p[C.kB55G] = 496.5636
    p[C.kB55Sub] = 0.0593
    p[C.kcAc2G] = 0.1916
    p[C.k25S] = 0.0050
    p[C.k25F] = 0.9411
    p[C.kweeS] = 0.0050
    p[C.kweeF] = 47.2937
    p[C.kcBc1W1] = 1.3132
    p[C.kcBc125] = 1.3132
    p[C.kppxY15] = 0.0050
    p[C.kcAc2W1] = 0.1096
    p[C.kcAc225] = 0.1096
    p[C.kB55W1] = 0.5511
    p[C.kB5525] = 0.5511
    p[C.Kd] = 0.025

    return p
end

function get_derivatives(u::Vector{Float64},p::Vector{Float64})
    # derivatives: dF/d[bifurcation_param]
    dFdp::Vector{Float64} = zeros(V.NUM)

    Wee1p = 1 - u[V.Wee1] - u[V.Wee1pp]
    Cdc25p = 1 - u[V.Cdc25] - u[V.Cdc25pp]

    dFdp[V.Subp] = - (1/p[C.Kd])*(p[C.kcBc1Sub]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*(p[C.SubT]-u[V.Subp])
    dFdp[V.PP1] = (1/p[C.Kd])*(p[C.kipp1] + p[C.kipp1C]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*u[V.PP1]
    dFdp[V.Gwlp] = - (1/p[C.Kd])*(p[C.kcBc1G]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*(p[C.Gwtot] - u[V.Gwlp])
    dFdp[V.Wee1] = (1/p[C.Kd])*(p[C.kcBc1W1]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*u[V.Wee1]
    dFdp[V.Wee1pp] = - (1/p[C.Kd])*(p[C.kcBc1W1]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*Wee1p
    dFdp[V.Cdc25] = (1/p[C.Kd])*(p[C.kcBc125]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*u[V.Cdc25]
    dFdp[V.Cdc25pp] = - (1/p[C.Kd])*(p[C.kcBc125]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*(Cdc25p)

    return dFdp
end

function get_steady_state(p::Vector{Float64})
    tspan::Tuple{Float64,Float64} = (0.0,Inf)
    u0::Vector{Float64} = zeros(V.NUM)
    u0[V.PP1] = 1.0
    u0[V.PP2AB55] = 0.25
    u0[V.Wee1] = 1.0
    u0[V.Cdc25] = 1.0

    prob = ODEProblem(diffeq!,u0,tspan,p)
    prob = SteadyStateProblem(prob)
    sol = solve(prob,DynamicSS(CVODE_BDF()),dt=1.0)

    return sol.u
end


text:
function collect_commom_DAT!(DAT, Ed, iJR)

    # maps
    # ----------------------------------------------------------------------------------

    # ----------------------------------------------------------------------------------
    max_model = iJR.load_model("max_model"; uncompress = false)
    objider = iJR.BIOMASS_IDER

    # ----------------------------------------------------------------------------------
    for exp in DAT[:EXPS]
        # exp dat
        Ed_biom = Ed.val("D", exp)
        DAT[:exp, :flx, "D", exp] = Ed_biom
        DAT[:exp, :err, "D", exp] = 0.0

        # bounds
        fva_model = iJR.load_model("fva_models", exp; uncompress = false)
        max_lb, max_ub = ChU.bounds(max_model, objider)
        fva_lb, fva_ub = ChU.bounds(fva_model, objider)
        lb = max(max_lb, fva_lb)
        ub = min(max_ub, fva_ub)
        DAT[:bounds, "D", exp] = (lb, ub)

        # exchanges
        for Ed_ider in DAT[:EXCH_FLX_IDERS]

            # exp dat
            Ed_flx = Ed.uval(Ed_ider, exp)
            Ed_err = try; Ed.uerr(Ed_ider, exp) catch; 0.0 end
            DAT[:exp, :flx, Ed_ider, exp] = Ed_flx
            DAT[:exp, :err, Ed_ider, exp] = Ed_err

        end

        # inner
        ider_map = iJR.load_rxns_map()
        for model_ider in DAT[:INNER_FLX_IDERS]

            Ed_ider = ider_map[model_ider]

            # exp dat
            Ed_flx = Ed.val(Ed_ider, exp)
            Ed_err = try; Ed.err(Ed_ider, exp) catch; 0.0 end
            DAT[:exp, :flx, Ed_ider, exp] = Ed_flx
            DAT[:exp, :err, Ed_ider, exp] = Ed_err

        end
    end

    return DAT
end


text:
using SciPy
using Test

@testset "spatial.KDTree" begin
    points1 = rand(15, 2)
    points2 = rand(15, 2)
    kd_tree1 = spatial.KDTree(points1)
    kd_tree2 = spatial.KDTree(points2)
    indexes = kd_tree1.query_ball_tree(kd_tree2, r=0.2)
    @test length(indexes) >= 1
end


text:
using Test
using TreeTools
using TreeKnit

MCCs_ref = [
	["A/NewYork/105/2002"],
	["A/NewYork/177/1999"],
	["A/NewYork/137/1999", "A/NewYork/138/1999"],
	["A/NewYork/52/2004", "A/NewYork/59/2003"],
	["A/NewYork/198/2003", "A/NewYork/199/2003", "A/NewYork/32/2003"],
	["A/NewYork/10/2004", "A/NewYork/11/2003", "A/NewYork/12/2003", "A/NewYork/13/2003", "A/NewYork/14/2003", "A/NewYork/15/2003", "A/NewYork/16/2003", "A/NewYork/17/2003", "A/NewYork/18/2003", "A/NewYork/19/2003", "A/NewYork/2/2003", "A/NewYork/21/2003", "A/NewYork/22/2003", "A/NewYork/23/2003", "A/NewYork/24/2003", "A/NewYork/25/2003", "A/NewYork/26/2003", "A/NewYork/27/2003", "A/NewYork/28/2003", "A/NewYork/29/2003", "A/NewYork/30/2003", "A/NewYork/31/2004", "A/NewYork/33/2004", "A/NewYork/34/2003", "A/NewYork/35/2003", "A/NewYork/36/2003", "A/NewYork/38/2003", "A/NewYork/39/2003", "A/NewYork/4/2003", "A/NewYork/40/2003", "A/NewYork/41/2003", "A/NewYork/42/2003", "A/NewYork/43/2003", "A/NewYork/44/2003", "A/NewYork/45/2003", "A/NewYork/46/2003", "A/NewYork/47/2003", "A/NewYork/48/2003", "A/NewYork/49/2003", "A/NewYork/5/2004", "A/NewYork/50/2003", "A/NewYork/51/2003", "A/NewYork/53/2003", "A/NewYork/54/2003", "A/NewYork/55/2003", "A/NewYork/56/2003", "A/NewYork/6/2004", "A/NewYork/60A/2003", "A/NewYork/61A/2003", "A/NewYork/62A/2003", "A/NewYork/63/2003", "A/NewYork/64/2003", "A/NewYork/65/2003", "A/NewYork/67/2003", "A/NewYork/69/2004", "A/NewYork/7/2003", "A/NewYork/70/2004", "A/NewYork/8/2003"],
]

t1 = read_tree("$(dirname(pathof(TreeKnit)))/../test/NYdata/tree_ha.nwk")
t2 = read_tree("$(dirname(pathof(TreeKnit)))/../test/NYdata/tree_na.nwk")

MCCs = computeMCCs(t1, t2)

@testset "NY" begin
	@test MCCs[1:end-1] == MCCs_ref
end
if MCCs[1:end-1] != MCCs_ref
	@warn "Found different MCCs for the NewYork data. Could indicate a problem..."
end

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
using SimplePCHIP   # WIll bring interpolate into namespace
using Test

function test_interpolation_is_piecewise_monotone(xs, ys, N=10000)
    itp = interpolate(xs, ys)
    for (i,j) in monotone_intervals(ys)
        @test is_monotone([itp(x) for x ∈ range(xs[i], stop=xs[j], length=N)])
    end
end

function monotone_intervals(xs)
    N = length(xs)
    function _monotone_intervals(ch :: Channel)
        up_down = 0
        previous_i = 1
        for i = 1:N-1
            current_x, next_x = xs[i:i+1]
            if up_down != 0  &&  up_down*current_x > up_down*next_x
                put!(ch, (previous_i, i))
                previous_i = i
                up_down = 0
            end
            if up_down == 0  &&  current_x != next_x
                up_down = current_x < next_x ? 1 : -1
            end
        end
        put!(ch, (previous_i, N))
    end

    [interval for interval in Channel(_monotone_intervals)]
end

function is_monotone(zs)
    (all(is_approx_le(z, znext) for (z,znext) = zip(zs[1:end-1], zs[2:end]))
    || all(is_approx_le(znext, z) for (z,znext) = zip(zs[1:end-1], zs[2:end])))
end

function is_approx_le(a, b)
    a < b || a ≈ b
end

xs = [0.0 1.0]
ys = [1.0 1.0]
test_interpolation_is_piecewise_monotone(xs, ys)

xs = [0.0 1.0]
ys = [1.0 0.0]
test_interpolation_is_piecewise_monotone(xs, ys)

xs = [0.0 1.0]
ys = [0.0 1.0]
test_interpolation_is_piecewise_monotone(xs, ys)

xs = [0.0 1.0 2.0 100.0]
ys = [0.0 0.0 1.0   1.0]
test_interpolation_is_piecewise_monotone(xs, ys)

xs = [1.0 2.0 5.0 6.0 10.0 22.0]
ys = [0.0 1.0 0.0 1.0  0.0  1.0]
test_interpolation_is_piecewise_monotone(xs, ys)

xs = [   0.0 10.0 10000.0]
ys = [-100.0  1.2     1.1]
test_interpolation_is_piecewise_monotone(xs, ys)

xs = range(0, stop=2π, length=5000)
ys = sin.(xs)
test_interpolation_is_piecewise_monotone(collect(xs), ys)

# Example data from Fritsch and Carlson, SIAM J. NUMER. ANAL. 17 (1980) 238-246.
xs = [ 0.0  2  3  5  6  8  9   11 12 14 15]
ys = [10.0 10 10 10 10 10 10.5 15 50 60 85]
test_interpolation_is_piecewise_monotone(xs, ys)

# Test the test helper functions
@test is_monotone([0.0 0.0 0.0])
@test is_monotone([0.0 1.0 10.0])
@test is_monotone([10.0 1.0 0.0])
@test ! is_monotone([0.0 1.0 0.0])
@test ! is_monotone([1.0 0.0 1.0])

# Test internal functions

# Make sure the correct interval is identified
p = SimplePCHIP._pchip(3, [1.0 2.0 3.0], [0.0 0.0 0.0], [0.0 0.0 0.0])
for search ∈ (SimplePCHIP._pchip_index_linear_search, SimplePCHIP._pchip_index_bisectional_search)
    @test search(p, 1.0) == 1
    @test search(p, 1.0 + eps(1.0)) == 1
    @test search(p, 2.0 - eps(2.0)) == 1
    @test search(p, 2.0) == 2
    @test search(p, 2.0 + eps(2.0)) == 2
    @test search(p, 3.0 - eps(3.0)) == 2
    @test search(p, 3.0) == 3
end

max_stars_repo_name: genkuroki/SimplePCHIP

max_issues_repo_name: PetrKryslUCSD/SimplePCHIP

max_forks_repo_name: PetrKryslUCSD/SimplePCHIP


text:
module Nanocsv

using DataFrames

export read_csv, write_csv

const QUOTE_CHAR = '"'

include("csvreader.jl")
include("csvwriter.jl")

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""
    ObjsComment(;
        comment=nothing,
        created=nothing,
        id=nothing,
        is_intro=nothing,
        is_starred=nothing,
        num_stars=nothing,
        pinned_info=nothing,
        pinned_to=nothing,
        reactions=nothing,
        timestamp=nothing,
        user=nothing,
    )

    - comment::String
    - created::Int32
    - id::DefsCommentId
    - is_intro::Bool
    - is_starred::Bool
    - num_stars::Int32
    - pinned_info::DefsPinnedInfo
    - pinned_to::Vector{DefsChannel}
    - reactions::Vector{ObjsReaction}
    - timestamp::Int32
    - user::DefsUserId
"""
mutable struct ObjsComment <: SwaggerModel
    comment::Any # spec type: Union{ Nothing, String } # spec name: comment
    created::Any # spec type: Union{ Nothing, Int32 } # spec name: created
    id::Any # spec type: Union{ Nothing, DefsCommentId } # spec name: id
    is_intro::Any # spec type: Union{ Nothing, Bool } # spec name: is_intro
    is_starred::Any # spec type: Union{ Nothing, Bool } # spec name: is_starred
    num_stars::Any # spec type: Union{ Nothing, Int32 } # spec name: num_stars
    pinned_info::Any # spec type: Union{ Nothing, DefsPinnedInfo } # spec name: pinned_info
    pinned_to::Any # spec type: Union{ Nothing, Vector{DefsChannel} } # spec name: pinned_to
    reactions::Any # spec type: Union{ Nothing, Vector{ObjsReaction} } # spec name: reactions
    timestamp::Any # spec type: Union{ Nothing, Int32 } # spec name: timestamp
    user::Any # spec type: Union{ Nothing, DefsUserId } # spec name: user

    function ObjsComment(;comment=nothing, created=nothing, id=nothing, is_intro=nothing, is_starred=nothing, num_stars=nothing, pinned_info=nothing, pinned_to=nothing, reactions=nothing, timestamp=nothing, user=nothing)
        o = new()
        validate_property(ObjsComment, Symbol("comment"), comment)
        setfield!(o, Symbol("comment"), comment)
        validate_property(ObjsComment, Symbol("created"), created)
        setfield!(o, Symbol("created"), created)
        validate_property(ObjsComment, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(ObjsComment, Symbol("is_intro"), is_intro)
        setfield!(o, Symbol("is_intro"), is_intro)
        validate_property(ObjsComment, Symbol("is_starred"), is_starred)
        setfield!(o, Symbol("is_starred"), is_starred)
        validate_property(ObjsComment, Symbol("num_stars"), num_stars)
        setfield!(o, Symbol("num_stars"), num_stars)
        validate_property(ObjsComment, Symbol("pinned_info"), pinned_info)
        setfield!(o, Symbol("pinned_info"), pinned_info)
        validate_property(ObjsComment, Symbol("pinned_to"), pinned_to)
        setfield!(o, Symbol("pinned_to"), pinned_to)
        validate_property(ObjsComment, Symbol("reactions"), reactions)
        setfield!(o, Symbol("reactions"), reactions)
        validate_property(ObjsComment, Symbol("timestamp"), timestamp)
        setfield!(o, Symbol("timestamp"), timestamp)
        validate_property(ObjsComment, Symbol("user"), user)
        setfield!(o, Symbol("user"), user)
        o
    end
end # type ObjsComment

const _property_map_ObjsComment = Dict{Symbol,Symbol}(Symbol("comment")=>Symbol("comment"), Symbol("created")=>Symbol("created"), Symbol("id")=>Symbol("id"), Symbol("is_intro")=>Symbol("is_intro"), Symbol("is_starred")=>Symbol("is_starred"), Symbol("num_stars")=>Symbol("num_stars"), Symbol("pinned_info")=>Symbol("pinned_info"), Symbol("pinned_to")=>Symbol("pinned_to"), Symbol("reactions")=>Symbol("reactions"), Symbol("timestamp")=>Symbol("timestamp"), Symbol("user")=>Symbol("user"))
const _property_types_ObjsComment = Dict{Symbol,String}(Symbol("comment")=>"String", Symbol("created")=>"Int32", Symbol("id")=>"DefsCommentId", Symbol("is_intro")=>"Bool", Symbol("is_starred")=>"Bool", Symbol("num_stars")=>"Int32", Symbol("pinned_info")=>"DefsPinnedInfo", Symbol("pinned_to")=>"Vector{DefsChannel}", Symbol("reactions")=>"Vector{ObjsReaction}", Symbol("timestamp")=>"Int32", Symbol("user")=>"DefsUserId")
Base.propertynames(::Type{ ObjsComment }) = collect(keys(_property_map_ObjsComment))
Swagger.property_type(::Type{ ObjsComment }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsComment[name]))}
Swagger.field_name(::Type{ ObjsComment }, property_name::Symbol) =  _property_map_ObjsComment[property_name]

function check_required(o::ObjsComment)
    (getproperty(o, Symbol("comment")) === nothing) && (return false)
    (getproperty(o, Symbol("created")) === nothing) && (return false)
    (getproperty(o, Symbol("id")) === nothing) && (return false)
    (getproperty(o, Symbol("is_intro")) === nothing) && (return false)
    (getproperty(o, Symbol("timestamp")) === nothing) && (return false)
    (getproperty(o, Symbol("user")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjsComment }, name::Symbol, val)
end


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.

using Base.Meta

_get_name(c::Symbol) = c
_get_name(c::Nothing) = ()
_get_name(c::AbstractString) = c
function _get_name(c::Expr)
    if c.head == :string
        return c
    else
        return c.args[1]
    end
end

"""
    _extract_kw_args(args)

Process the arguments to a macro, separating out the keyword arguments.
Return a tuple of (flat_arguments, keyword arguments, and requested_container),
where `requested_container` is a symbol to be passed to `parse_container`.
"""
function _extract_kw_args(args)
    kw_args =
        filter(x -> isexpr(x, :(=)) && x.args[1] != :container, collect(args))
    flat_args = filter(x -> !isexpr(x, :(=)), collect(args))
    requested_container = :Auto
    for kw in args
        if isexpr(kw, :(=)) && kw.args[1] == :container
            requested_container = kw.args[2]
        end
    end
    return flat_args, kw_args, requested_container
end

function _try_parse_idx_set(arg::Expr)
    # [i=1] and x[i=1] parse as Expr(:vect, Expr(:(=), :i, 1)) and
    # Expr(:ref, :x, Expr(:kw, :i, 1)) respectively.
    if arg.head === :kw || arg.head === :(=)
        @assert length(arg.args) == 2
        return true, arg.args[1], arg.args[2]
    elseif isexpr(arg, :call) && (arg.args[1] === :in || arg.args[1] === :∈)
        return true, arg.args[2], arg.args[3]
    else
        return false, nothing, nothing
    end
end
function _explicit_oneto(index_set)
    s = Meta.isexpr(index_set, :escape) ? index_set.args[1] : index_set
    if Meta.isexpr(s, :call) &&
       length(s.args) == 3 &&
       s.args[1] == :(:) &&
       s.args[2] == 1
        return :(Base.OneTo($index_set))
    else
        return index_set
    end
end

function _expr_is_splat(ex::Expr)
    if ex.head == :(...)
        return true
    elseif ex.head == :escape
        return _expr_is_splat(ex.args[1])
    end
    return false
end
_expr_is_splat(::Any) = false

"""
    _parse_ref_sets(expr::Expr)

Helper function for macros to construct container objects. Takes an `Expr` that
specifies the container, e.g. `:(x[i=1:3,[:red,:blue],k=S; i+k <= 6])`, and
returns:

    1. `idxvars`: Names for the index variables, e.g. `[:i, gensym(), :k]`
    2. `idxsets`: Sets used for indexing, e.g. `[1:3, [:red,:blue], S]`
    3. `condition`: Expr containing any conditional imposed on indexing, or `:()` if none is present
"""
function _parse_ref_sets(_error::Function, expr::Expr)
    c = copy(expr)
    idxvars = Any[]
    idxsets = Any[]
    # `:(t[i,j;k])` is a :ref, while `:(t[i;j])` is a :typed_vcat.
    # In both cases `:t` is the first arg.
    if isexpr(c, :typed_vcat) || isexpr(c, :ref)
        popfirst!(c.args)
    end
    condition = :()
    if isexpr(c, :vcat) || isexpr(c, :typed_vcat)
        # Parameters appear as plain args at the end.
        if length(c.args) > 2
            _error("Unsupported syntax $c.")
        elseif length(c.args) == 2
            condition = pop!(c.args)
        end # else no condition.
    elseif isexpr(c, :ref) || isexpr(c, :vect)
        # Parameters appear at the front.
        if isexpr(c.args[1], :parameters)
            if length(c.args[1].args) != 1
                _error(
                    "Invalid syntax: $c. Multiple semicolons are not " *
                    "supported.",
                )
            end
            condition = popfirst!(c.args).args[1]
        end
    end
    if isexpr(c, :vcat) || isexpr(c, :typed_vcat) || isexpr(c, :ref)
        if isexpr(c.args[1], :parameters)
            @assert length(c.args[1].args) == 1
            condition = popfirst!(c.args).args[1]
        end # else no condition.
    end

    for s in c.args
        parse_done = false
        if isa(s, Expr)
            parse_done, idxvar, _idxset = _try_parse_idx_set(s::Expr)
            if parse_done
                idxset = esc(_idxset)
            end
        end
        if !parse_done # No index variable specified
            idxvar = gensym()
            idxset = esc(s)
        end
        if idxvar in idxvars
            _error(
                "The index $(idxvar) appears more than once. The index " *
                "associated with each set must be unique.",
            )
        end
        push!(idxvars, idxvar)
        push!(idxsets, idxset)
    end
    return idxvars, idxsets, condition
end
_parse_ref_sets(_error::Function, expr) = (Any[], Any[], :())

"""
    _build_ref_sets(_error::Function, expr)

Helper function for macros to construct container objects. Takes an `Expr` that
specifies the container, e.g. `:(x[i=1:3,[:red,:blue],k=S; i+k <= 6])`, and
returns:

    1. `idxvars`: Names for the index variables, e.g. `[:i, gensym(), :k]`
    2. `indices`: Iterators over the indices indexing, e.g.
       `Constainers.NestedIterators((1:3, [:red,:blue], S), (i, ##..., k) -> i + k <= 6)`.
"""
function _build_ref_sets(_error::Function, expr)
    idxvars, idxsets, condition = _parse_ref_sets(_error, expr)
    if any(_expr_is_splat.(idxsets))
        _error(
            "cannot use splatting operator `...` in the definition of an index set.",
        )
    end
    has_dependent = has_dependent_sets(idxvars, idxsets)
    if has_dependent || condition != :()
        esc_idxvars = esc.(idxvars)
        idxfuns = [
            :(($(esc_idxvars[1:(i-1)]...),) -> $(idxsets[i])) for
            i in 1:length(idxvars)
        ]
        if condition == :()
            indices = :(Containers.nested($(idxfuns...)))
        else
            condition_fun = :(($(esc_idxvars...),) -> $(esc(condition)))
            indices =
                :(Containers.nested($(idxfuns...); condition = $condition_fun))
        end
    else
        indices =
            :(Containers.vectorized_product($(_explicit_oneto.(idxsets)...)))
    end
    return idxvars, indices
end

function container_code(idxvars, indices, code, requested_container)
    if isempty(idxvars)
        return code
    end
    if !(
        requested_container in
        [:Auto, :Array, :DenseAxisArray, :SparseAxisArray]
    )
        # We do this two-step interpolation, first into the string, and then
        # into the expression because interpolating into a string inside an
        # expression has scoping issues.
        error_message =
            "Invalid container type $requested_container. Must be " *
            "Auto, Array, DenseAxisArray, or SparseAxisArray."
        return :(error($error_message))
    end
    if requested_container == :DenseAxisArray
        requested_container = :(Containers.DenseAxisArray)
    elseif requested_container == :SparseAxisArray
        requested_container = :(Containers.SparseAxisArray)
    end
    esc_idxvars = esc.(idxvars)
    func = :(($(esc_idxvars...),) -> $code)
    if requested_container == :Auto
        return :(Containers.container($func, $indices))
    else
        return :(Containers.container($func, $indices, $requested_container))
    end
end
function parse_container(_error, var, value, requested_container)
    idxvars, indices = _build_ref_sets(_error, var)
    return container_code(idxvars, indices, value, requested_container)
end

"""
    @container([i=..., j=..., ...], expr)

Create a container with indices `i`, `j`, ... and values given by `expr` that
may depend on the value of the indices.

    @container(ref[i=..., j=..., ...], expr)

Same as above but the container is assigned to the variable of name `ref`.

The type of container can be controlled by the `container` keyword. Note that
when the index set is explicitly given as `1:n` for any expression `n`, it is
transformed to `Base.OneTo(n)` before being given to [`container`](@ref).

## Examples

```jldoctest; setup=:(using JuMP)
julia> Containers.@container([i = 1:3, j = 1:3], i + j)
3×3 Array{Int64,2}:
 2  3  4
 3  4  5
 4  5  6

julia> I = 1:3
1:3

julia> Containers.@container([i = I, j = I], i + j)
2-dimensional DenseAxisArray{Int64,2,...} with index sets:
    Dimension 1, 1:3
    Dimension 2, 1:3
And data, a 3×3 Array{Int64,2}:
 2  3  4
 3  4  5
 4  5  6

julia> Containers.@container([i = 2:3, j = 1:3], i + j)
2-dimensional DenseAxisArray{Int64,2,...} with index sets:
    Dimension 1, 2:3
    Dimension 2, Base.OneTo(3)
And data, a 2×3 Array{Int64,2}:
 3  4  5
 4  5  6

julia> Containers.@container([i = 1:3, j = 1:3; i <= j], i + j)
JuMP.Containers.SparseAxisArray{Int64,2,Tuple{Int64,Int64}} with 6 entries:
  [1, 2]  =  3
  [2, 3]  =  5
  [3, 3]  =  6
  [2, 2]  =  4
  [1, 1]  =  2
  [1, 3]  =  4
```
"""
macro container(args...)
    args, kw_args, requested_container = _extract_kw_args(args)
    @assert length(args) == 2
    @assert isempty(kw_args)
    var, value = args
    code = parse_container(error, var, esc(value), requested_container)
    anonvar = isexpr(var, :vect) || isexpr(var, :vcat)
    if anonvar
        return code
    else
        name = Containers._get_name(var)
        return :($(esc(name)) = $code)
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# misc file input/output functions

Base.adjoint(s::T) where T <: AbstractString = s # not sure where else to put this

function readDlmFileNextLine!(file::IOStream; delim::Char = delimiter)
	try
		return readdlm(IOBuffer(readline(file)), delim)
	catch
		return []
	end
end

function readDlmFile(filename::String; delim::Char = delimiter)
	@assert(isfile(filename), "file not found: $filename")
	file = open(filename, "r")
	data = readdlm(file, delim)
	close(file)
	return data
end

# open file for writing after checking that it does not already exist
function openNewFile(filename::String)
	@assert(!isfile(filename), "file already exists: $filename")
	file = open(filename, "w")
	return file
end

# write a line to file, add delimiter between args
function writeDlmLine!(file::IOStream, args...; delim::Char = delimiter)
	for arg in args
		write(file, string(arg), delim)
	end
	write(file, newline)
	# writedlm(file, hcat(args...), delim) # several times slower
end

# create dict with array entries as keys, indices as values
function arrayDict(array::Vector{Any})
	# return Dict(array[i] => i for i = 1:length(array)) # slower for some reason?
	dict = Dict()
	for i = 1:length(array)
		dict[array[i]] = i
	end
	return dict
end

struct Table
	name::String
	header::Vector{Any} # header names for columns in data, header[i] is for data[:,i]
	headerDict::Dict{Any,Int} # headerDict[h] gives index of h in header
	data::Array{Any,2}
	# columns::Dict{Any,SubArray{Any,1}} # columns[h] gives subarray column data for header h
	columns::Dict{Any,Vector{Any}}
	# columns[h] is same as data[:,headerDict[h]]

	function Table(name, header, data)
		name = convert(String, name)
		header = convert(Vector{Any}, header)
		data = convert(Array{Any,2}, data)

		(numRows, numCols) = size(data)
		@assert(length(header) == numCols)
		@assert(allunique(header))

		headerDict = arrayDict(header)
		columns = Dict{Any,Vector{Any}}()
		for i = 1:numCols
			columns[header[i]] = data[:,i]
			# columns[header[i]] = view(data, :, i)
		end

		return new(name, header, headerDict, data, columns)
	end

	function Table(name, header; rows = [], cols = [])
		@assert(isempty(rows) + isempty(cols) == 1)
		if !isempty(rows)
			return Table(name, header, adjoint(hcat(rows...)))
		elseif !isempty(cols)
			return Table(name, header, hcat(cols...))
		end
	end
end

function writeTablesToFile!(file::IOStream, table::Table;
	writeNumRows::Bool = false, writeNumCols::Bool = false)

	(numRows, numCols) = size(table.data)
	writeDlmLine!(file, table.name, writeNumRows ? numRows : "", writeNumCols ? numCols : "")
	writeDlmLine!(file, table.header...)
	writedlm(file, table.data, delimiter)
	writeDlmLine!(file, "") # newline
end

function writeTablesToFile!(file::IOStream, tables::Vector{Table})
	for table in tables
		writeTablesToFile!(file, table)
	end
end

function writeTablesToFile!(file::IOStream, tables::Dict{String, Table})
	for (name, table) in tables
		writeTablesToFile!(file, table)
	end
end

function writeTablesToFile(filename::String, tables::Vector{Table})
	file = open(filename, "w")
	writeTablesToFile!(file, tables)
	close(file)
end

function writeTablesToFile(filename::String, table::Table)
	writeTablesToFile(filename, [table])
end

function writeTablesToFile(filename::String, tables::Dict{String, Table})
	file = open(filename, "w")
	writeTablesToFile!(file, tables)
	close(file)
end

function readTablesFromFile(filename::String)
	return readTablesFromData(readDlmFile(filename))
end

function readTablesFromData(data::Array{Any,2})
	(numDataRows, numDataCols) = size(data)

	tables = Dict{String, Table}()
	i = 1
	while i <= numDataRows
		# find start of table
		while i <= numDataRows && isempty(data[i,1])
			i += 1
		end
		if i > numDataRows
			break
		end

		tableName = data[i,1]

		# read row and column counts, if given
		numRows = numCols = nullIndex
		(size(data,2) >= 2) && !isempty(data[i,2]) && (numRows = data[i,2])
		(size(data,2) >= 3) && !isempty(data[i,3]) && (numCols = data[i,3])

		i += 1 # move to table header

		# determine column count
		j = 1
		while j <= numDataCols && !isempty(data[i,j])
			j += 1
		end
		numCols == nullIndex ? numCols = j - 1 : @assert(numCols == j - 1)

		# read table header
		tableHeader = data[i,1:numCols]

		i += 1 # move to table data

		# determine row count
		j = i
		while j <= numDataRows && !isempty(data[j,1])
			j += 1
		end
		numRows == nullIndex ? numRows = j - i : @assert(numRows == j - i)

		# create table
		tableData = data[i-1 .+ (1:numRows), 1:numCols]
		table = Table(convert(String,tableName), tableHeader, tableData)

		# add table to tables
		@assert(!haskey(tables, table.name))
		tables[table.name] = table

		# go to next table
		i += numRows # i = j
		@assert(i > numDataRows || isempty(data[i,1]))
	end

	return tables
end

# given a table and names of fields from table header to use,
# return a vector of dicts for each row, with each dict mapping from
# a field name to the value for that row and field
function tableRowsFieldDicts(table::Table, fieldNames::Vector{T}) where T
	fieldNames = convert(Vector{String}, fieldNames)
	fieldsDict = Dict{String,Any}([f => nothing for f in fieldNames])
	numRows = size(table.data, 1)
	rowsFieldDicts = [deepcopy(fieldsDict) for i = 1:numRows] # rowsFieldDicts[i] = dict with fieldName => table.data[i,table.headerDict[fieldName]] for fieldName in fieldNames
	for fieldName in fieldNames
		j = table.headerDict[fieldName]
		for i = 1:numRows
			rowsFieldDicts[i][fieldName] = table.data[i,j]
		end
	end
	return rowsFieldDicts
end

# parse json attributes column from table
function parseAttributesColumn(table::Table)
	if haskey(table.columns, "attributes")
		return [isempty(s) ? Dict{String,Any}() : JSON.parse(s) for s in table.columns["attributes"]]
	else
		return [Dict{String,Any}() for i = 1:size(table.data, 1)]
	end
end

# crc checksum
function fileChecksum(filename::String)
	return CRC32c.crc32c(Mmap.mmap(filename))
end

function serializeToFile(filename::String, data::Any)
	file = open(filename, "w")
	serialize(file, data)
	close(file)
end

function deserializeFile(filename::String)
	file = open(filename, "r")
	data = deserialize(file)
	close(file)
	return data
end

# if path is absolute then return it, otherwise prepend with (assumed) absolute path
joinPathIfNotAbs(absPath::String, path::String) = isabspath(path) ? abspath(path) : joinpath(absPath, path)

function interpolateString(s::String)
	return string("\"", escape_string(s), "\"") |> Meta.parse |> eval # note that eval only works for global vars
end
interpolateString(s::SubString{String}) = interpolateString(String(s))

# some convenient functions for reading xml files
xmlFileRoot(filename::String) = root(parse_file(interpolateString(filename)))
findElt = find_element # shorthand
xName(x) = LightXML.name(x) # to avoid conflict with JuMP.name
containsElt(parentElt::XMLElement, eltString::String) = findElt(parentElt, eltString) != nothing
eltContent(elt::XMLElement) = content(elt)
eltContentVal(elt::XMLElement) = eval(Meta.parse(eltContent(elt)))
eltContentInterpVal(elt::XMLElement) = interpolateString(eltContent(elt))
eltContent(parentElt::XMLElement, eltString::String) = try content(findElt(parentElt, eltString));
	catch; error("Element not found: $eltString"); end
eltContentVal(parentElt::XMLElement, eltString::String) = eval(Meta.parse(eltContent(parentElt, eltString)))
eltContentInterpVal(parentElt::XMLElement, eltString::String) = interpolateString(eltContent(parentElt, eltString))
eltAttr = attribute
eltAttrVal(elt::XMLElement, attrString::String) = eval(Meta.parse(eltAttr(elt, attrString)))

function childrenNodeNames(parentElt::XMLElement)
	childNodes = Vector{String}()
	for childNode in child_nodes(parentElt)
		if is_elementnode(childNode)
			push!(childNodes, xName(childNode))
		end
	end
	return childNodes
end

function selectXmlFile(; message::String = "Enter xml filename: ")
	if Sys.iswindows()
		ps1Filename = "$sourceDir/file/select_xml.ps1"
		str = read(`Powershell.exe -executionpolicy remotesigned -File $ps1Filename`, String)
	else
		print(message)
		str = readline()
	end
	xmlFilename = abspath(String(strip(chomp(str), [' ', '"']))) #"

	# check selected file
	@assert(ispath(xmlFilename), "file does not exist: '$xmlFilename'")
	@assert(split(xmlFilename, '.')[end] == "xml", "file is not xml file: '$xmlFilename'")

	return xmlFilename
end


text:
import Plots.px
using Test
using MechanicalUnits
using MechanicalUnits.Unitfu: DimensionError
using Plots
using MechGluePlots
using Plots: default, px
using RecipesBase
using RecipesBase: debug
using Plots
using MechGluePlots: accumulate_unit_info, sertype
using MechGluePlots: print_prettyln, SeriesUnitInfo

default(width = 5)

if !@isdefined vq
    # 'vector' real, quantity
    vr = range(1, 3; step = 0.2)
    vq = range(1, 3; step = 0.2)s
    # Function_range_domain

    f_r_r(t)           = sin(0.3∙2π∙t)     ∙ 9.81
    f_q_q(t::Quantity) = sin(0.3∙2π∙t / s) ∙ 9.81N
    f_q_r(t::Quantity) = sin(0.3∙2π∙t / s)  ∙ 9.81
    f_r_q(t)           = sin(0.3∙2π∙t)     ∙ 9.81N

    s1x = range(-1, 1, length = 4)
    s1y = [-0.75, -0., 0.2, 1.0]
    s1z = s1x + s1y
    s2x = range(0.9, -0.9, length = 4)
    s2y = [-0.5, -0.1, 0.3, 0.5]
    s2z = s2x + s2y

    mxr = hcat(s1x, s2x)
    myr = hcat(s1y, s2y)
    mxq = hcat(s1x∙m, s2x∙m)
    myq = hcat(s1y∙s, s2y∙s)
    mxqm = hcat(s1x∙m, s2x∙m²)
    myqm = hcat(s1y∙s, s2y∙s²)

    s_1_down  = range(20, 6, length = 8)
    s_1_up = range(6, 20, length = 8)
    s_2_down = range(8, 1, length = 8)
    s_2_up =  range(1, 8, length = 8)

end


text:
using SimilaritySearch
using TextModel
using SnowballStemmer
using JSON

function set_oracle(index, fastlinks::Dict{UInt64,KnnResult})
    function oracle(q::VBOW)::Vector{Int32}
        L = Int32[]
        for term in q.tokens
            if haskey(fastlinks, term.id)
                for p in fastlinks[term.id]
                    push!(L, p.objID)
                end
            end
        end

        if length(L) == 0
            # just link randomly for orthogonal vectors
            n = length(index.db)
            return rand(1:n, floor(Int, log2(n+1)))
        end

        # @show L, fastlinks
        L
    end

    index.options.oracle = oracle
    fastlinks
end

function update_oracle(index, fastlinks, bow::VBOW, num_fast_links)
    for term in bow.tokens
        if !haskey(fastlinks, term.id)
            fastlinks[term.id] = KnnResult(num_fast_links)
        end
        push!(fastlinks[term.id], length(index.db), term.weight)
    end
end

function create_index()
    index = LocalSearchIndex(VBOW, angle_distance, recall=0.90, neighborhood=LogSatNeighborhood(1.5))
    fastlinks = set_oracle(index, Dict{UInt64,KnnResult}())
    # index.options.verbose = false
    return index, fastlinks
end

function create_normalizer(stemming, tabufiles, config)
    tabu = Set{String}()
    if length(tabufiles) > 0
        for filename in split(tabufiles, ':')
            words = split(TextModel.normalize_text(readstring(filename), config) |> join)
            union!(tabu, words)
        end
    end

    emptystring = ""
    if length(stemming) > 0
        S = Stemmer(stemming)
        if length(tabu) > 0
            info("CREATING stemming && tabu $(length(tabu))")
            return (w) -> if w ∈ tabu
                emptystring
            else
                stem(S, w)
            end
        else
            return (w) -> stem(S, w)
        end
    else
        if length(tabu) > 0
            return (w) -> if w ∈ tabu
                emptystring
            else
                w
            end
        else
            return identity
        end
    end
end

function main(filename, dup_threshold, key, num_fast_links; stemming="", verbose=false, tabu="")
    config = TextConfig()
    config.nlist = [1]
    config.qlist = []
    config.skiplist = []
    config.normalize = create_normalizer(stemming, tabu, config)
    index, fastlinks = create_index()
    index.options.verbose = verbose

    lineno = 0

    DB = []

    iterlines(filename) do line
        tweet = TextModel.parsetweet(line)
        lineno += 1
        rawbow = compute_bow(tweet[key], config)
        # info(rawbow)
        bow = VBOW(rawbow)
        knn, N = find_neighborhood(index, bow)
        if length(knn) > 0 && first(knn).dist < dup_threshold
            if verbose
                info("- DISCARDED:", line)
                info("-  ORIGINAL:", DB[first(knn).objID])
                info("-  DISTANCE:", first(knn).dist, " < ", dup_threshold)
            end
        else
            push_neighborhood!(index, bow, N, length(index.db))
            update_oracle(index, fastlinks, bow, num_fast_links)  # adding fast links
            push!(DB, line)
            println(line)
        end

        if lineno % 1000 == 0
            info("selected $(length(DB)) of $(lineno), drop-ratio: $(length(DB) / lineno)")
        end

    end
end

if length(ARGS) == 0
    info("""
    neardup-simsearch.jl removes the number of near duplicated items in text databases. It keeps only the original version of the document.

    Usage: [environment-variables] julia neardup-simsearch.jl files...

    Each input file must contain a json-dictionary per line, the indexes are independent per file

    For simplicity, the arguments are passed as environment-variables

    - dup_threshold: items under this value are considered as duplicated
       valid range: 0 < dup_threshold < pi/2
       default: 1.0

    - key: the keyword containing the text for each json
       default: text

    - num_fast_links: the number of documents stored per token to help localsearch indexes to perform fast lookups
       default: 7

    - stemming: specifies stemming rules to apply.
      e.g. stemming=spanish
       valid values: "", "danish", "dutch", "english", "finnish", "french", "german", "hungarian", "italian", "norwegian", "porter", "portuguese", "romanian", "russian", "spanish", "swedish", "turkish"
    - verbose: if enabled (verbose=true) the output becomes too noisy :S
    """)
else
    for filename in ARGS
        main(
            filename,
            parse(Float64, get(ENV, "dup_threshold", "1.0")),
            get(ENV, "key", "text"),
            parse(Int, get(ENV, "num_fast_links", "7")),
            verbose=parse(Bool, get(ENV, "verbose", "false")),
            tabu=get(ENV, "tabu", ""),
            stemming=get(ENV, "stemming", ""),
        )
    end
end

max_stars_repo_name: INGEOTEC/PatternAnalysisTools

max_issues_repo_name: INGEOTEC/PatternAnalysisTools

max_forks_repo_name: INGEOTEC/PatternAnalysisTools


text:
lo = 2
hi = 10
klim = 16

print("n-step Fibonacci for n = (", lo, ",", hi)
println(") up to k = ", klim, ":")
for i in 2:10
    print(@sprintf("%5d => ", i))
    for j in fib(i, klim)
        print(j, " ")
    end
    println()
end

println()
print("n-step Rosetta Code Lucas for n = (", lo, ",", hi)
println(") up to k = ", klim, ":")
for i in 2:10
    print(@sprintf("%5d => ", i))
    for j in luc_rc(i, klim)
        print(j, " ")
    end
    println()
end

println()
print("n-step MathWorld Lucas for n = (", lo, ",", hi)
println(") up to k = ", klim, ":")
for i in 2:10
    print(@sprintf("%5d => ", i))
    for j in luc(i, klim)
        print(j, " ")
    end
    println()
end