text:
###############################################################
# Comparison
###############################################################

(==)(a::MultiValue,b::MultiValue) = a.data == b.data
(≈)(a::MultiValue,b::MultiValue) = isapprox(get_array(a), get_array(b))
(≈)(a::MultiValue{S,T1,N,0} where T1,b::MultiValue{S,T2,N,0} where T2) where {S,N} = true

function (≈)(
  a::AbstractArray{<:MultiValue}, b::AbstractArray{<:MultiValue})
  if size(a) != size(b); return false; end
  for (ai,bi) in zip(a,b)
    if !(ai≈bi); return false; end
  end
  true
end

function isless(a::MultiValue{Tuple{L}},b::MultiValue{Tuple{L}}) where L
  for d in L:-1:1
    if a[d] < b[d]
      return true
    elseif a[d] > b[d]
      return false
    else
      continue
    end
  end
  false
end

isless(a::Number,b::MultiValue) where {D,T} = all(a .< b.data)

###############################################################
# Addition / subtraction
###############################################################

for op in (:+,:-)
  @eval begin

    function ($op)(a::T) where {T<:MultiValue}
      r = map($op, a.data)
      T(r)
    end

    function ($op)(a::MultiValue{S},b::MultiValue{S})  where S
      r = broadcast(($op), a.data, b.data)
      T = _eltype($op,r,a,b)
      M = change_eltype(a,T)
      M(r)
    end

    function ($op)(a::TensorValue,b::SymTensorValue)
      @notimplemented
    end

    function ($op)(a::SymTensorValue,b::TensorValue)
      @notimplemented
    end

  end
end

###############################################################
# Matrix Division
###############################################################

function (\)(a::MultiValue{Tuple{D,D}} where D, b::MultiValue)
  r = get_array(a) \ get_array(b)
  T = change_eltype(b,eltype(r))
  T(r)
end

###############################################################
# Operations with other numbers
###############################################################

for op in (:+,:-,:*)
  @eval begin
    function ($op)(a::MultiValue,b::Number)
        r = broadcast($op,a.data,b)
        T = _eltype($op,r,a,b)
        M  = change_eltype(a,T)
        M(r)
    end

    function ($op)(a::Number,b::MultiValue)
        r = broadcast($op,a,b.data)
        T = _eltype($op,r,a,b)
        M  = change_eltype(b,T)
        M(r)
    end
  end
end

function (/)(a::MultiValue,b::Number)
    r = broadcast(/,a.data,b)
    T = _eltype(/,r,a,b)
    P  = change_eltype(a,T)
    P(r)
end

@inline function _eltype(op,r,a,b)
  eltype(r)
end

@inline function _eltype(op,r::Tuple{},a,b)
  typeof(op(zero(eltype(a)),zero(eltype(b))))
end

###############################################################
# Dot product (simple contraction)
###############################################################

function (*)(a::MultiValue, b::MultiValue)
  msg = """
  Method (*)(::$(typeof(a)),::$(typeof(b))) has been removed.
  Depending the case, use simple contraction dot aka ⋅ (\\cdot) or full contraction inner aka ⊙ (\\odot) instead.
  """
  error(msg)
  #dot(a,b)
end

dot(a::MultiValue{Tuple{D}}, b::MultiValue{Tuple{D}}) where D = inner(a,b)

dot(a::MultiValue,b::MultiValue) = @notimplemented

@generated function dot(a::A,b::B) where {A<:MultiValue{Tuple{D1}},B<:MultiValue{Tuple{D1,D2}}} where {D1,D2}
    ss = String[]
    for j in 1:D2
      s = ""
      for i in 1:D1
        ak = data_index(A,i)
        bk = data_index(B,i,j)
        s *= "a.data[$ak]*b.data[$bk]+"
      end
        push!(ss,s[1:(end-1)]*", ")
    end
    str = join(ss)
    Meta.parse("VectorValue{$D2}($str)")
end

@generated function dot(a::A,b::B) where {A<:MultiValue{Tuple{D1,D2}},B<:MultiValue{Tuple{D2}}} where {D1,D2}
    ss = String[]
    for i in 1:D1
      s = ""
      for j in 1:D2
        ak = data_index(A,i,j)
        bk = data_index(B,j)
        s *= "a.data[$ak]*b.data[$bk]+"
      end
        push!(ss,s[1:(end-1)]*", ")
    end
    str = join(ss)
    Meta.parse("VectorValue{$D1}($str)")
end

@generated function dot(a::MultiValue{Tuple{D1,D3}}, b::MultiValue{Tuple{D3,D2}}) where {D1,D2,D3}
    ss = String[]
    for j in 1:D2
        for i in 1:D1
            s = join([ "a[$i,$k]*b[$k,$j]+" for k in 1:D3])
            push!(ss,s[1:(end-1)]*", ")
        end
    end
    str = join(ss)
    Meta.parse("TensorValue{$D1,$D2}(($str))")
end

# Double contraction

#(::Colon)(a::MultiValue{Tuple{D1,D2}},b::MultiValue{Tuple{D1,D2}}) where {D1,D2} = inner(a,b)
#(::Colon)(a::MultiValue{Tuple{D1,D2}},b::MultiValue{Tuple{D1,D2,D3,D4}}) where {D1,D2,D3,D4} = inner(a,b)
#(::Colon)(a::MultiValue{Tuple{D1,D2,D3,D4}},b::MultiValue{Tuple{D1,D2}}) where {D1,D2,D3,D4} = inner(a,b)

###############################################################
# Inner product (full contraction)
###############################################################

inner(a::Real,b::Real) = a*b

function inner(a::MultiValue, b::MultiValue)
  @notimplemented
end

@generated function inner(a::MultiValue{S}, b::MultiValue{S}) where S
    str = join([" a[$i]*b[$i] +" for i in 1:length(a) ])
    Meta.parse(str[1:(end-1)])
end

@generated function inner(a::SymTensorValue{D}, b::SymTensorValue{D}) where D
  str = ""
  for i in 1:D
    for j in 1:D
      k = data_index(a,i,j)
      str *= " a.data[$k]*b.data[$k] +"
    end
  end
  Meta.parse(str[1:(end-1)])
end

@generated function inner(a::SymFourthOrderTensorValue{D}, b::SymTensorValue{D}) where D
  str = ""
  for i in 1:D
    for j in i:D
      s = ""
      for k in 1:D
        for l in 1:D
          ak = data_index(a,i,j,k,l)
          bk = data_index(b,k,l)
          s *= " a.data[$ak]*b.data[$bk] +"
        end
      end
      str *= s[1:(end-1)]*", "
    end
  end
  Meta.parse("SymTensorValue{D}($str)")
end

function inner(a::SymFourthOrderTensorValue{D},b::MultiValue{Tuple{D,D}}) where D
  inner(a,symmetric_part(b))
end

const ⊙ = inner

###############################################################
# Reductions
###############################################################

for op in (:sum,:maximum,:minimum)
    @eval begin
        $op(a::MultiValue) = $op(a.data)
    end
end

# Outer product (aka dyadic product)

"""
"""
outer(a::Real,b::Real) = a*b
outer(a::MultiValue,b::Real) = a*b
outer(a::Real,b::MultiValue) = a*b

function outer(a::MultiValue,b::MultiValue)
   @notimplemented
end

@generated function outer(a::MultiValue{Tuple{D}},b::MultiValue{Tuple{Z}}) where {D,Z}
    str = join(["a[$i]*b[$j], " for j in 1:Z for i in 1:D])
    Meta.parse("TensorValue{$D,$Z}($str)")
end

function outer(a::VectorValue{0,Ta},b::VectorValue{1,Tb}) where {Ta,Tb}
  T = promote_type(Ta,Tb)
  TensorValue{0,1,T}()
end

function outer(a::VectorValue{0,Ta},b::Tb) where {Ta,Tb<:Real}
  T = promote_type(Ta,Tb)
  VectorValue{0,T}()
end

@generated function outer(a::MultiValue{Tuple{D}},b::MultiValue{Tuple{D1,D2}}) where {D,D1,D2}
  str = join(["a[$i]*b[$j,$k], "  for k in 1:D2 for j in 1:D1 for i in 1:D])
  Meta.parse("ThirdOrderTensorValue{D,D1,D2}($str)")
end

@generated function outer(a::SymTensorValue{D},b::SymTensorValue{D}) where D
  str = ""
  for i in 1:D
    for j in i:D
      ak = data_index(a,i,j)
      for k in 1:D
        for l in k:D
          bk = data_index(b,k,l)
          str *= "a.data[$ak]*b.data[$bk], "
        end
      end
    end
  end
  Meta.parse("SymFourthOrderTensorValue{D}($str)")
end

const ⊗ = outer

###############################################################
# Cross Product
###############################################################

function cross(a::MultiValue{Tuple{3}}, b::MultiValue{Tuple{3}})
  VectorValue{3}(a[2]b[3]-a[3]b[2], a[3]b[1]-a[1]b[3], a[1]b[2]-a[2]b[1])
end

function cross(a::MultiValue{Tuple{2}}, b::MultiValue{Tuple{2}})
  a[1]b[2]-a[2]b[1]
end

cross(a::MultiValue,b::MultiValue) = error("Cross product only defined for R2 and R3 vectors")

###############################################################
# Linear Algebra
###############################################################

det(a::MultiValue{Tuple{D1,D2}}) where {D1,D2} = det(get_array(a))

det(a::MultiValue{Tuple{1,1}}) = a[1]

function det(a::MultiValue{Tuple{2,2}})
  a_11 = a[1,1]; a_12 = a[1,2]
  a_21 = a[2,1]; a_22 = a[2,2]
  a_11*a_22 - a_12*a_21
end

function det(a::MultiValue{Tuple{3,3}})
  a_11 = a[1,1]; a_12 = a[1,2]; a_13 = a[1,3]
  a_21 = a[2,1]; a_22 = a[2,2]; a_23 = a[2,3]
  a_31 = a[3,1]; a_32 = a[3,2]; a_33 = a[3,3]
  a_11*a_22*a_33 + a_12*a_23*a_31 + a_13*a_21*a_32 -
    (a_11*a_23*a_32 + a_12*a_21*a_33 + a_13*a_22*a_31)
end

inv(a::MultiValue{Tuple{D1,D2}}) where {D1,D2} = TensorValue(inv(get_array(a)))

function inv(a::MultiValue{Tuple{1,1}})
  r = 1/a[1]
  T = change_eltype(a,typeof(r))
  T(r)
end

function inv(a::MultiValue{Tuple{2,2}})
 c = 1/det(a)
 data = (a[2,2]*c, -a[2,1]*c, -a[1,2]*c,  a[1,1]*c)
 TensorValue{2}(data)
end

function inv(a::MultiValue{Tuple{3,3}})
 a_11 = a[1,1]; a_12 = a[1,2]; a_13 = a[1,3]
 a_21 = a[2,1]; a_22 = a[2,2]; a_23 = a[2,3]
 a_31 = a[3,1]; a_32 = a[3,2]; a_33 = a[3,3]
 c = 1/det(a)
 data = (
    ( a_22*a_33 - a_23*a_32 )*c,
   -( a_21*a_33 - a_23*a_31 )*c,
    ( a_21*a_32 - a_22*a_31 )*c,
   -( a_12*a_33 - a_13*a_32 )*c,
    ( a_11*a_33 - a_13*a_31 )*c,
   -( a_11*a_32 - a_12*a_31 )*c,
    ( a_12*a_23 - a_13*a_22 )*c,
   -( a_11*a_23 - a_13*a_21 )*c,
    ( a_11*a_22 - a_12*a_21 )*c)
 TensorValue{3}(data)
end

###############################################################
# Measure
###############################################################

"""
"""
meas(a::MultiValue{Tuple{D}}) where D = sqrt(inner(a,a))
meas(a::MultiValue{Tuple{D,D}}) where D = abs(det(a))
meas(a::TensorValue{0,1,T}) where T = one(T)

function meas(v::MultiValue{Tuple{1,2}})
  n1 = v[1,2]
  n2 = -1*v[1,1]
  n = VectorValue(n1,n2)
  sqrt(n ⋅ n)
end

function meas(v::MultiValue{Tuple{2,3}})
  n1 = v[1,2]*v[2,3] - v[1,3]*v[2,2]
  n2 = v[1,3]*v[2,1] - v[1,1]*v[2,3]
  n3 = v[1,1]*v[2,2] - v[1,2]*v[2,1]
  n = VectorValue(n1,n2,n3)
  sqrt(n ⋅ n)
end

@inline norm(u::MultiValue{Tuple{D}}) where D = sqrt(inner(u,u))
@inline norm(u::MultiValue{Tuple{0},T}) where T = sqrt(zero(T))

###############################################################
# conj
###############################################################

function conj(a::T) where {T<:MultiValue}
  r = map(conj, a.data)
  T(r)
end

###############################################################
# Trace
###############################################################

@generated function tr(v::MultiValue{Tuple{D,D}}) where D
    str = join([" v[$i,$i] +" for i in 1:D ])
    Meta.parse(str[1:(end-1)])
end

@generated function tr(v::MultiValue{Tuple{A,A,B}}) where {A,B}
  lis = LinearIndices((A,A,B))
  str = ""
  for k in 1:B
    for i in 1:A
      if i !=1
        str *= " + "
      end
      p = lis[i,i,k]
      str *= " v.data[$p]"
    end
    str *= ", "
  end
  Meta.parse("VectorValue($str)")
end

###############################################################
# Adjoint and transpose
###############################################################

adjoint(a::MultiValue{Tuple{D,D}}) where D = @notimplemented
transpose(a::MultiValue{Tuple{D,D}}) where D = @notimplemented

@generated function adjoint(a::TensorValue{D1,D2}) where {D1,D2}
  str = ""
  for i in 1:D1
    for j in 1:D2
      k = (j-1)*D1 + i
      str *= "conj(a.data[$k]), "
    end
  end
  Meta.parse("TensorValue{D2,D1}($str)")
end

@generated function transpose(a::TensorValue{D1,D2}) where {D1,D2}
  str = ""
  for i in 1:D1
    for j in 1:D2
      k = (j-1)*D1 + i
      str *= "a.data[$k], "
    end
  end
  Meta.parse("TensorValue{D2,D1}($str)")
end

@inline function adjoint(a::TensorValue{D1,D2,T}) where {D1,D2,T<:Real}
  transpose(a)
end

adjoint(a::SymTensorValue) = conj(a)

@inline adjoint(a::SymTensorValue{D,T} where {D,T<:Real}) = transpose(a)

transpose(a::SymTensorValue) = a

###############################################################
# Symmetric part
###############################################################

"""
"""
@generated function symmetric_part(v::MultiValue{Tuple{D,D}}) where D
    str = "("
    for j in 1:D
        for i in j:D
            str *= "0.5*v[$i,$j] + 0.5*v[$j,$i], "
        end
    end
    str *= ")"
    Meta.parse("SymTensorValue{D}($str)")
end

###############################################################
# Define new operations for Gridap types
###############################################################

for op in (:symmetric_part,)
    @eval begin
        ($op)(a::GridapType) = operate($op,a)
    end
end

for op in (:inner,:outer)#,:(:))
    @eval begin
        ($op)(a::GridapType,b::GridapType) = operate($op,a,b)
        ($op)(a::GridapType,b::Number)     = operate($op,a,b)
        ($op)(a::Number,    b::GridapType) = operate($op,a,b)
        ($op)(a::GridapType,b::Function)   = operate($op,a,b)
        ($op)(a::Function,  b::GridapType) = operate($op,a,b)
    end
end


text:
module demo

export mandelbrot, julia

function mandelbrot(x, y, n=100)
  c = complex(x, y)
  z = zero(c)
  for j ∈ 0:n-1
    z = z^2 + c
    abs(z) > 2.0 && (return j/n)
  end
  return 1.0
end

function julia(x, y, n=100; c = -0.54 + 0.54im)
  z = complex(x, y)
  for j ∈ 0:n-1
    z = z^2 + c
    abs(z) > 2.0 && (return j/n)
  end
  return 1.0
end

end # module

ext: jl


text:
# using DrWatson
# quickactivate(@__DIR__)
# @show projectdir()
# using Plots
##
using CovidSurvey
using UnPack
using Test
using Turing
using Optim

θ = (
	num_rt_steps      = 100,
	lockdown_index    = 30,
	num_time_steps    = 130,
	num_observations  = 130,
	link              = KLogistic(3.),
    invlink           = KLogit(3.),
	population        = 5_837_213,
	num_impute        = 6,
	serial_interval   = National.serialinterval(15),
	num_si            = 15,
)

rt  = let
	n = θ.num_rt_steps
	N = θ.lockdown_index+n
	t = 1:N

	frt(t) = t <= θ.lockdown_index ? 2. : frt1(t)
	frt1(t) = 0.5 + (t-θ.lockdown_index-1)*1/100
	rt  = frt.(t)
	# plot(rt)|>display
	rt
end

## =============================================================================
# test RandomWalk
@testset "latent Rt" begin

	@model function test_rt(Y, θ, predict, ::Type{TV} = Vector{Float64}) where TV
		@unpack invlink, num_rt_steps, num_time_steps = θ
		R0         ~ truncated(Normal(2., .1), 1.5, 2.5)
		R1         ~ truncated(Normal(.5, .1), .3, .7)
		σ_rt       ~ truncated(Normal(0.1, .08), 0, .25)
	    latent_Rt  ~ National.RandomWalk(num_rt_steps, σ_rt, invlink(R1))

		Rt = TV(undef, num_time_steps)
		National.random_walks!(Rt, θ, predict, latent_Rt, R0, σ_rt)
		Y ~ MvNormal(Rt, 1e-5)
		return Rt
	end

	m = test_rt(rt, θ, false)
	@test length(m()) == θ.num_time_steps
	mle = optimize(m, MLE(), LBFGS(), Optim.Options(iterations = 1000, allow_f_increases = true))
	@test mle.optim_result.ls_success

	estimated_Rt = let #generated quantity
		v = mle.values
		latent_Rt = [v[Symbol("latent_Rt[$i]")] for i in 1:θ.num_rt_steps]
		Rt = Vector{Float64}(undef, θ.num_time_steps)
		National.random_walks!(Rt, θ, false, latent_Rt, v[:R0], v[:σ_rt])
		Rt
	end
	# plot!(estimated_Rt)
	@test isapprox(rt, estimated_Rt, atol=0.1)
end
## ============================================================================
# test latent infections
y = 100.
newly_infected       = Vector{Float64}(undef, θ.num_time_steps)
cumulative_infected  = Vector{Float64}(undef, θ.num_time_steps)
effective_Rt         = Vector{Float64}(undef, θ.num_time_steps)
National.infections!(newly_infected, cumulative_infected, effective_Rt, θ, y, rt)

@testset "latent infections with given Rt" begin
	@model function test_latent_infections(Y, θ, Rt, predict, ::Type{TV} = Vector{Float64}) where TV
		@unpack num_time_steps = θ

		############ infection dynamics
		τ  ~ Exponential(100/3)
		T  = typeof(τ)
		y  ~ truncated(Exponential(τ),T(0),T(1000))

		newly_infected       = TV(undef, num_time_steps)
		cumulative_infected  = TV(undef, num_time_steps)
		effective_Rt         = TV(undef, num_time_steps)

		National.infections!(newly_infected, cumulative_infected, effective_Rt, θ, y, Rt)
		Y ~ MvNormal(newly_infected, 1e-1)
		return newly_infected
	end

	m = test_latent_infections(newly_infected, θ, rt, false)
	@test length(m()) == θ.num_time_steps
	mle = optimize(m, MLE(), LBFGS(), Optim.Options(iterations = 1000, allow_f_increases = true))
	@test mle.optim_result.ls_success

	estimated_newly_infected  = let
		y = mle.values[:y]
		newly_infected       = Vector{Float64}(undef, θ.num_time_steps)
		cumulative_infected  = Vector{Float64}(undef, θ.num_time_steps)
		effective_Rt         = Vector{Float64}(undef, θ.num_time_steps)
		National.infections!(newly_infected, cumulative_infected, effective_Rt, θ, y, rt)
		newly_infected
	end
	# p = plot(newly_infected)
	# plot!(estimated_newly_infected)
	# @test all( estimated_newly_infected .≈ newly_infected )
	@test isapprox(estimated_newly_infected, newly_infected, atol=0.1)
end

## ============================================================================
# test latent infections #2
@testset "latent infections without Rt" begin
	@model function test_latent_infections_v2(Y, θ, predict, ::Type{TV} = Vector{Float64}) where TV
		@unpack invlink, num_rt_steps, num_time_steps = θ

		############# time varying reproduction number
	    R0         ~ truncated(Normal(2., .1), 1.5, 2.5)
		R1         ~ truncated(Normal(.5, .1), .3, .7)
		σ_rt       ~ truncated(Normal(0.1, .08), 0, .25)
	    latent_Rt  ~ National.RandomWalk(num_rt_steps, σ_rt, invlink(R1))

		Rt = TV(undef, num_time_steps)
		National.random_walks!(Rt, θ, predict, latent_Rt, R0, σ_rt)

		############ infection dynamics
		# τ  ~ Exponential(100/3)
		# T  = typeof(τ)
		y  ~ InverseGamma2(100, .1)#truncated(Normal(500, 200),T(10),T(1000))

		newly_infected       = TV(undef, num_time_steps)
		cumulative_infected  = TV(undef, num_time_steps)
		effective_Rt         = TV(undef, num_time_steps)

		National.infections!(newly_infected, cumulative_infected, effective_Rt, θ, y, Rt)
		Y ~ MvNormal(newly_infected, 1e-5)
		return newly_infected
	end

	m = test_latent_infections_v2(newly_infected, θ, false)
	@test length(m()) == θ.num_time_steps
	mle = optimize(m, MLE(), LBFGS(), Optim.Options(iterations = 1000, allow_f_increases = true))
	@test mle.optim_result.ls_success

	estimated_newly_infected, estimated_Rt  = let
		v = mle.values
		latent_Rt = [v[Symbol("latent_Rt[$i]")] for i in 1:θ.num_rt_steps]
		Rt = Vector{Float64}(undef, θ.num_time_steps)
		National.random_walks!(Rt, θ, false, latent_Rt, v[:R0], v[:σ_rt])

		y = v[:y]
		newly_infected       = Vector{Float64}(undef, θ.num_time_steps)
		cumulative_infected  = Vector{Float64}(undef, θ.num_time_steps)
		effective_Rt         = Vector{Float64}(undef, θ.num_time_steps)
		National.infections!(newly_infected, cumulative_infected, effective_Rt, θ, y, Rt)
		newly_infected, Rt
	end

	# plot(rt)
	# plot!(estimated_Rt)
	@test isapprox(rt, estimated_Rt, atol=0.1)

	# p = plot(newly_infected)
	# plot!(estimated_newly_infected)
	# @test all( estimated_newly_infected .≈ newly_infected )
	@test isapprox(estimated_newly_infected, newly_infected, atol=0.1)
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
export MI

#=
Bias correction techniques

=#

#=
Quadratic extrapolation
Strong et al 1998

I=I_raw + a/n_s + b/n_s^2

Compute MI with number of repetitions (n_s) equal to total, total/2, and total/4
Plot result against 1/n_s, and fit quadratic polynomial by least squares
Intercept is true information
=#

function MI{B<:QE, T<:transformation}(info::Information{B,T},response::Array{Float64,2},stimulus::Array{Float64,1})

    N_s=length(stimulus)

    myinds=1:N_s

    samplepoints=round(Int64,[N_s; N_s/2; N_s/4])

    MI=zeros(Float64,3)

    (conmat, p_s, p_r)=transform(info,response,stimulus)

    MI[1]=mutual_information(conmat,p_s,p_r)

    for i=2:3

        theseinds=sample(myinds, samplepoints[i], replace=false)

        stimulus_2=stimulus[theseinds]
        response_2=response[theseinds,:]

        (conmat2, p_s2, p_r2)=transform(info,response_2,stimulus_2)

        MI[i]=mutual_information(conmat2,p_s2,p_r2)

    end

    I_corrected=polyfit(1./samplepoints, MI, 2)[1]

    (I_corrected, conmat)

end

#=
Panzeri-Treves (PT)
Panzeri and Treves 1996
=#

#=
Nemenman-Shafee-Bialek (NSB)
Nemenman et al 2002
=#

#=
James-Stein shrinkage
Hausser and Strimmer 2008
=#

function mutual_information(p_rs::Array{Float64,2},p_s::Array{Float64,1},p_r::Array{Float64,1})

    MI=0.0

    for i=1:size(p_rs,1)
        for j=1:size(p_rs,2)
            if (p_rs[i,j]==0.0) || (isnan(p_rs[i,j]))
            elseif p_r[j]==0.0
            else
                MI+=p_rs[i,j]*p_s[i]*log2(p_rs[i,j]/(p_r[j]))
            end

        end
    end

    MI

end

function transform{B<:bias,T<:decoder}(info::Information{B,T},response::Array{Float64,2},stimulus::Array{Float64,1})
    predict=project(info.t,response,stimulus)
    (conmat, p_s, p_r)=conmatrix(info.t,predict,stimulus)
end

function polyfit(x::Array{Float64,1}, y::Array{Float64,1}, n::Int64)
  A = [ x[i]^p for i = 1:length(x), p = 0:n ]
  A \ y
end


text:
# indexnotation/tensormacro.jl
#
# Defines the @tensor macro which switches to an index-notation environment.
const prime = Symbol("'")

macro tensor(arg)
    tensorify(arg)
end

function tensorify(ex::Expr)
    if ex.head == :(=) || ex.head == :(:=) || ex.head == :(+=) || ex.head == :(-=)
        lhs = ex.args[1]
        rhs = ex.args[2]
        if isa(lhs, Expr) && lhs.head == :ref
            dst = tensorify(lhs.args[1])
            src = ex.head == :(-=) ? tensorify(Expr(:call,:-,rhs)) : tensorify(rhs)
            indices = makeindex_expr(lhs)
            if ex.head == :(:=)
                return :($dst = deindexify($src, $indices))
            else
                value = ex.head == :(=) ? 0 : +1
                return :(deindexify!($dst, $src, $indices, $value))
            end
        end
    end
    if ex.head == :ref
        indices = makeindex_expr(ex)
        t = tensorify(ex.args[1])
        return :(indexify($t,$indices))
    end
    if ex.head == :call && ex.args[1] == :scalar
        if length(ex.args) != 2
            error("scalar accepts only a single argument")
        end
        src = tensorify(ex.args[2])
        indices = :(Indices{()}())
        return :(scalar(deindexify($src, $indices)))
    end
    return Expr(ex.head,map(tensorify,ex.args)...)
end
tensorify(ex::Symbol) = esc(ex)
tensorify(ex) = ex

function makeindex_expr(ex::Expr)
    if ex.head == :ref
        for i = 2:length(ex.args)
            if isa(ex.args[i],Expr) && ex.args[i].head == prime
                ex.args[i] = makesymbolprime(ex.args[i])
            end
            isa(ex.args[i],Int) || isa(ex.args[i],Symbol) || isa(ex.args[i],Char) || error("cannot make indices from $ex")
        end
    else
        error("cannot make indices from $ex")
    end
    return :(Indices{$(tuple(ex.args[2:end]...))}())
end

function makesymbolprime(ex::Expr)
    if isa(ex,Expr) && ex.head == prime && length(ex.args) == 1
        if isa(ex.args[1],Symbol) || isa(ex.args[1],Int)
            return Symbol(ex.args[1],prime)
        elseif isa(ex.args[1],Expr) && ex.args[1].head == prime
            return Symbol(makesymbolprime(ex.args[1]),prime)
        end
    end
    error("cannot make indices from $ex")
end


text:
# harbor-model.jl
# SimLynx.jl Harbor Model

using SimLynx
SimLynx.greet()

using Distributions: Exponential, Uniform
using Random

const MEAN_INTERARRIVAL_TIME = 4.0 / 3.0 # days
const MIN_UNLOADING_TIME = 1.0 # days
const MAX_UNLOADING_TIME = 2.5 # days

cycle_time = nothing

dock = nothing
queue = nothing

@process scheduler() begin
    dist = Exponential(MEAN_INTERARRIVAL_TIME)
    i = 1
    while true
        @schedule in 0.0 ship(i)
        work(rand(dist))
        i += 1
    end
end

@process ship(i::Integer) begin
    dist = Uniform(MIN_UNLOADING_TIME, MAX_UNLOADING_TIME)
    arrival_time = current_time()
    current_process_store!(:unloading_time, rand(dist))
    if !harbor_master(current_process(), :arriving)
        enqueue!(queue, current_process())
        suspend()
    end
    work(current_process_store(:unloading_time))
    remove!(dock, current_process())
    cycle_time.value = current_time() - arrival_time
    harbor_master(current_process(), :leaving)
    return nothing
end

function harbor_master(ship::Process, action::Symbol)
    if action == :arriving
        if length(dock) < 2
            # The dock is not full
            if isempty(dock)
                process_store!(ship,
                               :unloading_time,
                               process_store(ship, :unloading_time) / 2)
            else
                other_ship = first(dock)
                notice = interrupt(other_ship)
               notice.time = current_time() + 2*notice.time
                resume(other_ship, notice)
            end
            add!(dock, ship)
            return true
        else
            # The dock is full
            return false
        end
    elseif action == :leaving
        if isempty(queue)
            if !isempty(dock)
                other_ship = first(dock)
                notice = interrupt(other_ship)
                notice.time = current_time() + 2/notice.time
                resume(other_ship, notice)
            end
        else
            next_ship = dequeue!(queue)
            enqueue!(dock, next_ship)
            resume(next_ship, Notice(current_time(), next_ship))
        end
        return true
    else
        error("harbor_master: illegal action value $action")
    end
end

@event stop_sim() begin
    println("Harbor Model - report after $(current_time()) days - $(cycle_time.stats.n) ships unloaded")
    println("Minimum unload time was $(cycle_time.stats.min)")
    println("Maximum unload time was $(cycle_time.stats.max)")
    println("Average unload time was $(cycle_time.stats.mean)")
    println("Average queue of ships waiting to be unloaded was $(queue.n.stats.mean)")
    println("Maximum queue of ships waiting to be unloaded was $(queue.n.stats.max)")
    plot_history(queue.n, title="Queue of Ships History")
    stop_simulation()
end

function run_simulation()
    @simulation begin
        current_trace!(true)
        global cycle_time = Variable{Float64}(data=:tally)
        global dock = Queue{Process}()
        global queue = FifoQueue{Process}()
        @schedule at 0.0 scheduler()
        @schedule at 80.0 stop_sim()
        start_simulation()
    end
end

run_simulation()


text:
using SpWord2Vec

init = SpWord2Vec.SpSkipGramA.init
sampleFactory = SpWord2Vec.SpSkipGramA.sampleFactory

init("/home/cyc/Data/project_data/SpSkipGram/pairs-wiki-1based-20k.idx")
factory = get(sampleFactory)

const producer = factory.producer
batch_set = Task(() -> producer(factory, 1024, 15))

st = time()
for i=1:10000
    consume(batch_set)
end
en = time()

speed = factory.consumed[]/(en-st)
println("speed: $speed bytes/s")

init("/home/cyc/Data/project_data/SpSkipGram/wiki-sg-norm-lc-drop.bz2")

factory = get(sampleFactory)

const producer_ = factory.producer
batch_set = Task(() -> producer_(factory, 1024, 15))

st = time()
for i=1:1000
    consume(batch_set)
end
en = time()

speed = factory.consumed[]/(en-st)
println("speed: $speed lines/s")

init("/home/cyc/Data/project_data/SpSkipGram/wiki-sg-norm-lc-drop.bz2")
factory = get(sampleFactory)
const producer__ = factory.producer
batch_set = Task(() -> producer__(factory, 1024, 7))

st = time()
for i=1:1000
    consume(batch_set)
end
en = time()

speed = factory.consumed[]/(en-st)
println("speed: $speed lines/s")


text:
using Distributed

const SUPPORTED_HDF5 = (Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Float32, Float64)

mutable struct HDF5Writer{C,T} <: Writer{C,T}
    filename::String
    groupname::String
    writemode::String
    writeoffset::Int
    collector::C # The collector used by the Writer.
    localbuffer::TypedBuffer{T} # Stores data entries until the buffer is full.
    remotequeue::RemoteChannel{Channel{TypedBuffer{T}}} # References the channel on the master
                                                        # to which the buffer should be flushed.
end

"""
    HDF5Writer(collector; chunksz = 15000, channels_per_worker = 2)

Return a `Writer` that writes the `collector`'s data to an HDF5 file.

# Arguments
* `filename::String`: the HDF5 file in which to store the data.
* `groupname::String`: the group in the HDF5 file in which to store the data.
* `collector::Collector`: the Collector to use to collect data.
* `chunksz::Integer`: specifies the number of sets of data that the collector's buffer can store
    before sending it off to the master process to write it to disk.
* `channels_per_worker::Integer`: specifies the number of channels that the master process should
    allocate for each worker.
* `writemode::String`: the writemode to use for `filename`. Can be either "r+" (read-write,
    preserving any existing contents) or "w" (read-write, destroying any existing contents, if any).
"""
function HDF5Writer(filename::String, groupname::String, collector::Collector{C};
                    chunksz::Integer = 15000, channels_per_worker::Integer = 2,
                    writemode::String = "w") where C
    for datatype in C.types
        @assert(in(datatype, SUPPORTED_HDF5),
                "collector contains a datatype that isn't supported by HDF5Writer, got: $datatype")
    end
    localbuffer = TypedBuffer(collector, chunksz)
    buffertype = typeof(localbuffer)
    remotequeue = RemoteChannel(() -> Channel{buffertype}(channels_per_worker * nworkers()), 1)
    return HDF5Writer(filename, groupname, writemode, 0, collector, localbuffer, remotequeue)
end

function open_sink(writer::HDF5Writer)
    hdf5file = h5open(writer.filename, writer.writemode)
    if exists(hdf5file, writer.groupname)
        close(hdf5file)
        error("the specified groupname already exists")
    else
        hdf5group = g_create(hdf5file, writer.groupname)
    end
    return hdf5group, hdf5file
end

function close_sink(sink::Tuple{HDF5Group,HDF5File}, writer::HDF5Writer)
    foreach(close, sink)
    return nothing
end

function get_return(writer::HDF5Writer)
    return nothing
end

# Uses a generated function to unroll the loop over the datavector tuple. This avoids the type
# instability that results from iterating over a heterogenous tuple.
@generated function insert!(sink::Tuple{HDF5Group,HDF5File}, writer::HDF5Writer,
                            buffer::TypedBuffer{T}) where T
    n_datavectors = length(T.parameters)

    ex = :()
    for datavector_idx in 1:n_datavectors
        datatype = T.parameters[datavector_idx].parameters[1]

        ex = :($(ex.args...);
            name = names[$datavector_idx];
            if exists(sink[1], name);
                dset = d_open(sink[1], name);
                set_dims!(dset, (newoffset, ));
            else;
                dset = d_create(sink[1], name, $datatype, ((size, ), (-1, )), "chunk",
                                (maxsize, ));
            end;
            dset[(writer.writeoffset + 1):newoffset] = view(data[$datavector_idx], 1:size);
            close(dset);)
    end

    return quote
        if buffer.size > 0
            newoffset = writer.writeoffset + buffer.size

            # Hoist from the loop.
            names = writer.collector.names
            data = buffer.data
            size = buffer.size
            maxsize = buffer.maxsize

            $ex

            writer.writeoffset = newoffset
        end
        return nothing
    end
end


text:
module LPC

using ..DSP: xcorr

using LinearAlgebra: dot

export lpc, LPCBurg, LPCLevinson

# Dispatch types for lpc()
mutable struct LPCBurg; end
mutable struct LPCLevinson; end

"""
    lpc(x::AbstractVector, p::Int, [LPCBurg()])

Given input signal `x` and prediction order `p`, returns IIR coefficients `a`
and average reconstruction error `prediction_err`. Note that this method does
NOT return the leading `1` present in the true autocorrelative estimate; it
omits it as it is implicit in every LPC estimate, and must be manually
reintroduced if the returned vector should be treated as a polynomial.

The algorithm used is determined by the last optional parameter, and can be
either `LPCBurg` or `LPCLevinson`.
"""
function lpc end

"""
    lpc(x::AbstractVector, p::Int, LPCBurg())

LPC (Linear-Predictive-Code) estimation, using the Burg method. This function
implements the mathematics published in [1].

[1] - Enhanced Partial Tracking Using Linear Prediction
(DAFX 2003 article, Lagrange et al)
http://www.sylvain-marchand.info/Publications/dafx03.pdf
"""
function lpc(x::AbstractVector{T}, p::Int, ::LPCBurg) where T <: Number
    ef = x                      # forward error
    eb = x                      # backwards error
    a = [1; zeros(T, p)]        # prediction coefficients

    # Initialize prediction error wtih the variance of the signal
    prediction_err = (sum(abs2, x) ./ length(x))[1]

    for m in 1:p
        efp = ef[1:end-1]
        ebp = eb[2:end]
        k = -2 * dot(ebp,efp) / (sum(abs2, ebp) + sum(abs2,efp))
        ef = efp + k .* ebp
        eb = ebp + k .* efp
        a[1:m+1] = [a[1:m]; 0] + k .* [0; a[m:-1:1]]
        prediction_err *= (1 - k*k)
    end

    return a[2:end], prediction_err
end

"""
    lpc(x::AbstractVector, p::Int, LPCLevinson())

LPC (Linear-Predictive-Code) estimation, using the Levinson method. This
function implements the mathematics described in [1].

[1] - The Wiener (RMS) Error Criterion in Filter Design and Prediction
(N. Levinson, Studies in Applied Mathematics 25(1946), 261-278,
https://doi.org/10.1002/sapm1946251261)
"""
function lpc(x::AbstractVector{<:Number}, p::Int, ::LPCLevinson)
    R_xx = xcorr(x,x)[length(x):end]
    a = zeros(p,p)
    prediction_err = zeros(1,p)

    # for m = 1
    a[1,1] = -R_xx[2]/R_xx[1]
    prediction_err[1] = R_xx[1]*(1-a[1,1]^2)

    # for m = 2,3,4,..p
    for m = 2:p
        a[m,m] = (-(R_xx[m+1] + dot(a[m-1,1:m-1],R_xx[m:-1:2]))/prediction_err[m-1])
        a[m,1:m-1] = a[m-1,1:m-1] + a[m,m] * a[m-1,m-1:-1:1]
        prediction_err[m] = prediction_err[m-1]*(1-a[m,m]^2)
    end

    # Return autocorrelation coefficients and error estimate
    a[p,:], prediction_err[p]
end

# Default users to using Burg estimation as it is in general more stable
lpc(x, p) = lpc(x, p, LPCBurg())

end # module LPC


text:
# # Two dimensional turbulence example
#
# In this example, we initialize a random velocity field and observe its turbulent decay
# in a two-dimensional domain. This example demonstrates:
#
#   * How to run a model with no tracers and no buoyancy model.
#   * How to use `AbstractOperations`.
#   * How to use `ComputedField`s to generate output.

# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, JLD2, Plots"
# ```

# ## Model setup

# We instantiate the model with an isotropic diffusivity. We use a grid with 128² points,
# a fifth-order advection scheme, third-order Runge-Kutta time-stepping,
# and a small isotropic viscosity.  Note that we assign `Flat` to the `z` direction.

using Oceananigans

grid = RegularRectilinearGrid(size=(128, 128), extent=(2π, 2π),
                              topology=(Periodic, Periodic, Flat))

model = NonhydrostaticModel(timestepper = :RungeKutta3,
                              advection = UpwindBiasedFifthOrder(),
                                   grid = grid,
                               buoyancy = nothing,
                                tracers = nothing,
                                closure = IsotropicDiffusivity(ν=1e-5)
                           )

# ## Random initial conditions
#
# Our initial condition randomizes `model.velocities.u` and `model.velocities.v`.
# We ensure that both have zero mean for aesthetic reasons.

using Statistics

u, v, w = model.velocities

uᵢ = rand(size(u)...)
vᵢ = rand(size(v)...)

uᵢ .-= mean(uᵢ)
vᵢ .-= mean(vᵢ)

set!(model, u=uᵢ, v=vᵢ)

# ## Computing vorticity and speed

# To make our equations prettier, we unpack `u`, `v`, and `w` from
# the `NamedTuple` model.velocities:
u, v, w = model.velocities

# Next we create two objects called `ComputedField`s that calculate
# _(i)_ vorticity that measures the rate at which the fluid rotates
# and is defined as
#
# ```math
# ω = ∂_x v - ∂_y u \, ,
# ```

ω = ∂x(v) - ∂y(u)

ω_field = ComputedField(ω)

# We also calculate _(ii)_ the _speed_ of the flow,
#
# ```math
# s = \sqrt{u^2 + v^2} \, .
# ```

s = sqrt(u^2 + v^2)

s_field = ComputedField(s)

# We'll pass these `ComputedField`s to an output writer below to calculate them during the simulation.
# Now we construct a simulation that prints out the iteration and model time as it runs.

progress(sim) = @info "Iteration: $(sim.model.clock.iteration), time: $(round(Int, sim.model.clock.time))"

simulation = Simulation(model, Δt=0.2, stop_time=50, iteration_interval=100, progress=progress)

# ## Output
#
# We set up an output writer for the simulation that saves the vorticity every 20 iterations.

simulation.output_writers[:fields] = JLD2OutputWriter(model, (ω=ω_field, s=s_field),
                                                      schedule = TimeInterval(2),
                                                      prefix = "two_dimensional_turbulence",
                                                      force = true)

# ## Running the simulation
#
# Pretty much just

run!(simulation)

# ## Visualizing the results
#
# We load the output and make a movie.

using JLD2

file = jldopen(simulation.output_writers[:fields].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

# Construct the ``x, y`` grid for plotting purposes,

xω, yω, zω = nodes(ω_field)
xs, ys, zs = nodes(s_field)
nothing # hide

# and animate the vorticity and fluid speed.

using Plots

@info "Making a neat movie of vorticity and speed..."

anim = @animate for (i, iteration) in enumerate(iterations)

    @info "Plotting frame $i from iteration $iteration..."

    t = file["timeseries/t/$iteration"]
    ω_snapshot = file["timeseries/ω/$iteration"][:, :, 1]
    s_snapshot = file["timeseries/s/$iteration"][:, :, 1]

    ω_lim = 2.0
    ω_levels = range(-ω_lim, stop=ω_lim, length=20)

    s_lim = 0.2
    s_levels = range(0, stop=s_lim, length=20)

    kwargs = (xlabel="x", ylabel="y", aspectratio=1, linewidth=0, colorbar=true,
              xlims=(0, model.grid.Lx), ylims=(0, model.grid.Ly))

    ω_plot = contourf(xω, yω, clamp.(ω_snapshot', -ω_lim, ω_lim);
                       color = :balance,
                      levels = ω_levels,
                       clims = (-ω_lim, ω_lim),
                      kwargs...)

    s_plot = contourf(xs, ys, clamp.(s_snapshot', 0, s_lim);
                       color = :thermal,
                      levels = s_levels,
                       clims = (0., s_lim),
                      kwargs...)

    plot(ω_plot, s_plot, title=["Vorticity" "Speed"], layout=(1, 2), size=(1200, 500))
end

mp4(anim, "two_dimensional_turbulence.mp4", fps = 8) # hide


text:
## Deps
using Plots
using StatsBase
using Statistics
using StatsPlots
using LinearAlgebra
using DelimitedFiles

# Default plotting
default(titlefont=(14,"arial"),guidefont=(16,"arial"),tickfont=(12,"arial"))

# Read in data
data_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Real-Data/GM12878/chr22/"
cpel_params = readdlm("$(data_dir)/GM12878_wgbs_cpelnano_theta.txt")[:,4]
cpel_params = map(x->parse.(Float64,split(x,",")),cpel_params)

# Get parameters
a = map(x->x[1],cpel_params)
b = map(x->x[2],cpel_params)
c = map(x->x[3],cpel_params)

a_range = -10.0:0.5:10.0
b_range = -99.5:10.0:99.5
c_range = -20.0:2.0:20.0

a_histo = fit(Histogram,a,a_range)
b_histo = fit(Histogram,b,b_range)
c_histo = fit(Histogram,c,c_range)

a_histo = normalize(a_histo, mode=:probability)
b_histo = normalize(b_histo, mode=:probability)
c_histo = normalize(c_histo, mode=:probability)

sts_a = collect(a_range)
sts_b = collect(b_range)
sts_c = collect(c_range)

cent_a = sts_a.+step(a_range)/2.0
cent_b = sts_b.+step(b_range)/2.0
cent_c = sts_c.+step(c_range)/2.0

p_a = plot(cent_a,a_histo.weights*100,seriestype=:bar,alpha=0.25,xlabel="a",ylabel="Percentage (%)",label="");
p_b = plot(cent_b,b_histo.weights*100,seriestype=:bar,alpha=0.25,xlabel="b",ylabel="Percentage (%)",label="");
p_c = plot(cent_c,c_histo.weights*100,seriestype=:bar,alpha=0.25,xlabel="c",ylabel="Percentage (%)",label="");

p = plot(p_a,p_b,p_c,layout=(3,1),size=(500,800));
savefig(p,"$(data_dir)/Histograms.png")
savefig(p,"$(data_dir)/Histograms.pdf")


text:
import JSON
using Sockets

function main()
  text = open("1.json") do file
    read(file, String)
  end
  jobj = JSON.parse(text)
  coordinates = jobj["coordinates"]
  len = length(coordinates)
  x = y = z = 0

  for coord in coordinates
    x += coord["x"]
    y += coord["y"]
    z += coord["z"]
  end

  println(x / len)
  println(y / len)
  println(z / len)
end

function test()
  x = @timed main()
  println("Elapsed: $(x[2]), Allocated: $(x[3]), GC Time: $(x[4])")
end

for i in 1:4
  test()
end

try
  socket = connect("localhost", 9001)
  write(socket, "Julia")
  close(socket)
catch
  # standalone usage
end

test()

max_stars_repo_name: wmeddie/benchmarks

max_issues_repo_name: wmeddie/benchmarks

max_forks_repo_name: wmeddie/benchmarks


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

# messages
struct Checkpoint{T,L}
    level::Int
    key::Symbol
    var::T
    from::L
end
struct Interval{T}
    x::T
end
struct Load
    from::AbstractString
end
struct Register{L}
    lk::L
end
struct Parent{L}
    lk::L
end
struct Reset
    key::Symbol
end
struct Restore{L}
    level::Int
    key::Symbol
    from::L
end
struct Save end
struct Snapshot end
struct Start end
struct Stop end

# Checkpointing actor
mutable struct Checkpointing{S,T}
    store::S
    level::Int
    file::AbstractString
    parent::Union{Link,Nothing}
    children::Vector{Link}
    interval::T
    timer::Vector{Timer}
end

function (cp::Checkpointing)(msg::Checkpoint)
    cp.store[msg.key] = msg.var
    if cp.level < msg.level && !isnothing(cp.parent)
        send(cp.parent, Checkpoint(msg.level, msg.key, msg.var, self()))
    elseif cp.level ≥ 2
        foreach(cp.children) do child
            if child != msg.from
                foreach(call(child)) do check
                    cp.store[check[1]] = check[2]
                end
            end
        end
    end
end
function (cp::Checkpointing)(msg::Interval)
    cp.interval = msg.x
    foreach(close, cp.timer)
    empty!(cp.timer)
    cp.interval.on && cp(Start())
end
function (cp::Checkpointing)(msg::Load)
    !isempty(cp.file) && open(msg.from, "r") do io
        cp.store = deserialize(io)
    end
end
(cp::Checkpointing)(msg::Parent) = cp.parent = msg.lk
(cp::Checkpointing)(msg::Register) = unique!(push!(cp.children, msg.lk))
function (cp::Checkpointing)(msg::Reset)
    msg.key == :all ?
        empty!(cp.store) :
        delete!(cp.store, msg.key)
end
function (cp::Checkpointing)(msg::Restore)
    if cp.level ≤ msg.level || isempty(cp.parent)
        send(msg.from, get(cp.store, msg.key, nothing))
    else
        send(cp.parent, msg) # propagate upwards
    end
end
function (cp::Checkpointing)(::Save)
    !isempty(cp.file) && open(cp.file, "w") do io
        serialize(io, cp.store)
    end
end
function (cp::Checkpointing)(::Start)
    !cp.interval.on && (cp.interval = merge(cp.interval, (; on=true)))
    if cp.interval.update ≥ 1
        push!(cp.timer, Timer(0, interval=cp.interval.update) do timer
            send(self(), Update())
        end)
    end
    if cp.interval.save ≥ 1 && !isempty(cp.file)
        push!(cp.timer, Timer(0, interval=cp.interval.save) do timer
            send(self(), Save())
        end)
    end
end
function (cp::Checkpointing)(::Stop)
    cp.interval = merge(cp.interval, (; on=false))
    foreach(close, cp.timer)
    empty!(cp.timer)
end
function (cp::Checkpointing)(::Snapshot)
    if cp.level ≥ 2
        foreach(cp.children) do child
            foreach(call(child)) do check
                cp.store[check[1]] = check[2]
            end
        end
    end
end
(cp::Checkpointing)() = cp.store

# Checkpointing API
"""
    checkpointing(level=1, filename::AbstractString=""; kwargs...)

Start a checkpointing actor and return a [`Link`](@ref) to it.

# Arguments
- `filename::AbstractString`: filename where the actor should save its
    checkpointing buffer,
- `on::Bool=false`: should checkpoints automatically be updated and saved,
- `update::Int=10`: update interval in seconds ≥ 1,
- `save::Int=60`: saving interval in seconds ≥ 1,
- `kwargs...`: keyword arguments to [`spawn`](@ref).

!!! note
    Updating or saving is only done automatically if their
    intervals are ≥ 1 seconds.
"""
function checkpointing(level::Int=1 ,filename::AbstractString="";
                        on::Bool=false, update::Int=10, save::Int=60, kwargs...)
    cp = spawn(Checkpointing(Dict{Symbol, Any}(),
                             level, filename, nothing, Link[],
                             (; on, update, save), Timer[]);
            kwargs...)
    on && send(cp, Start())
    return cp
end

"""
    @chkey a b c 123

Build a checkpointing key from the surrounding module and function
names and the given arguments. This is intended for easy construction
of checkpointing keys.

# Example
```julia
module MyModule

using Actors

function mykey()
    a = 1
    b = 2
    c = 3
    # do something
    return @chkey a b c 123
end

export mykey
end # MyModule

julia> using .MyModule

julia> mykey()
:MyModule_mykey_a_b_c_123
```
"""
macro chkey(exs...)
    return quote
        st = stacktrace(backtrace())
        func = +
        for frm in st
            if frm.func != :backtrace && frm.func != Symbol("macro expansion")
                func = frm.func
                break
            end
        end
        ret = Symbol($(__module__),"_",func)
        for ex in $exs
            ret = Symbol(ret, "_", sprint(Base.show_unquoted, ex))
        end
        ret
    end
end

"""
    checkpoint(cp::Link, key::Symbol, args...)

Tell a checkpointing actor to take a checkpoint.

# Arguments
- `cp::Link`: [`Link`](@ref) to the checkpointing actor,
- `key::Symbol`: key for the checkpoint,
- `args...`: variables to save in the checkpoint,
- `level::Int=1`: checkpoint level.
"""
checkpoint(cp::Link, key::Symbol, args...; level::Int=1) =
    send(cp, Checkpoint(level, key, args, nothing))

"""
    restore(cp::Link; key::Symbol)

Tell a checkpointing actor to restore the last taken checkpoint.

# Arguments
- `cp::Link`:  [`Link`](@ref) to the checkpointing actor,
- `key::Symbol`: checkpoint key,
- `level::Int=1`: checkpoint level to restore from.
"""
function restore(cp::Link, key::Symbol; level::Int=1)
    me = newLink()
    send(cp, Restore(level, key, me))
    return receive(me)
end

"Tell the checkpointing actor `cp` to return its stored checkpoint data."
get_checkpoints(cp::Link) = call(cp)

"Tell a checkpointing actor `cp` to save its stored checkpoint data."
save_checkpoints(cp::Link) = send(cp, Save())

"Tell a checkpointing actor `cp` to load checkpointing data from a file"
load_checkpoints(cp::Link, fromfile::AbstractString) = send(cp, Load(fromfile))

"""
    set_interval(cp::Link; kwargs...)

Set the checkpointing intervals of a checkpointing actor `cp`.

# Arguments
- `cp::Link`: [`Link`](@ref) to a checkpointing actor,
- `kwargs...`: allowed keyword arguments:
    - `on:Bool`: should checkpoints automatically be updated and saved,
    - `update::Int`: update interval in seconds ≥ 1,
    - `save::Int`: saving interval in seconds ≥ 1.
"""
function set_interval(cp::Link; kwargs...)
    intv = get_interval(cp)
    kw = NamedTuple(kwargs)
    res = merge(interval, kw)
    typeof(res) == typeof(interval) ?
        send(cp, Interval(res)) :
        throw(AssertionError("forbidden key/value in $kw"))
end

"""
Get the checkpointing intervals of a checkpointing actor.
"""
get_interval(cp::Link) = diag(cp, :bhv).interval

"""
Register a checkpointing actor to a higher level actor.
"""
register_checkpoint(cp::Link) = send(cp, Register(self()))

"""
Start periodic checkpointing.
"""
start_checkpointing(cp::Link) = send(cp, Start())

"""
Stop periodic checkpointing.
"""
stop_checkpointing(cp::Link) = send(cp, Stop())


text:
# boot_biaslp.m
# Luca Brugnolini 2014
# Based on Kilian code

# This function gives bootstrap bias estimates for each horizon of LP

function boot_biaslp(ys, yt, x, rfirf, block, sub)

# Input : dependent variables ys, lagged RHS variables (yt and x)
#         vec(reduced form IRF by LP) for a certain horizon
#         sub=1: sub-bootstrap (need to draw from pre-selected blocks)
#         sub=0: 1st round bootstrap (standard blocks-of-blocks bootstrap)
# Output: bias estimates for each horizon (K^2 by 1)

(T,K) = size(yt);

#  Bootstrap simulation for bias estimation
if sub == 0 nrep = 100                                # number of bootstrap replications sub==0 1000

elseif sub ==1  nrep = 100;                          # number of bootstrap replications sub==1 500

end

IRFrmat = zeros(nrep,K^2);               # IRF for each bootstrap replication

for j = 1:nrep
    pos = ceil(rand()*(size(ys,1)-block))+1;
    ysr = ys[pos:pos+block-1,:];
    ytr = yt[pos:pos+block-1,:];
    xr  = x[pos:pos+block-1,:];

# I  DO NOT UNDERSTAND THE TWO FOR CYCLE for ii BECAUSE THERE ARE NO ii IN THE CODE

if sub == 1

   for ii = 1:ceil(T/block)
   pos = block*ceil(rand()*(T/block-1))+1;
   ysr = [ysr; ys[pos:pos+block-1,:]];
   ytr = [ytr; yt[pos:pos+block-1,:]];
   xr  = [xr ;  x[pos:pos+block-1,:]];

   end

elseif sub == 0

       for ii = 1:ceil(size(ys,1)/block)
       pos = ceil(rand()*(size(ys,1)-block))+1;
       ysr = [ysr; ys[pos:pos+block-1,:]];
       ytr = [ytr; yt[pos:pos+block-1,:]];
       xr  = [xr ;  x[pos:pos+block-1,:]];

       end

end

    IRFrmat[j,:] = vec(irflpest(ysr, ytr, xr));         # IRF by LP for 1~H

end

bias = (mean(IRFrmat,1))'-rfirf;                        # bias estimates

return bias

end


text:
"""
    lebedev_points(n)

Return the lebedev points on the unit sphere.
"""
function lebedev_points(n::Int)
    pts = Vector{SVector{3,Float64}}()
    if n<=6
        push!(pts,_sph_pt(0,90))
        push!(pts,_sph_pt(180,90))
        push!(pts,_sph_pt(90,90))
        push!(pts,_sph_pt(-90,90))
        push!(pts,_sph_pt(90,0))
        push!(pts,_sph_pt(90,180))
        return pts
    elseif n<=14
        push!(pts,lebedev_points(6)...)
        push!(pts,_sph_pt(45,   54.735610317245346))
        push!(pts,_sph_pt(45,   125.264389682754654))
        push!(pts,_sph_pt(-45,  54.735610317245346))
        push!(pts,_sph_pt(-45,  125.264389682754654))
        push!(pts,_sph_pt(135,  54.735610317245346))
        push!(pts,_sph_pt(135,  125.264389682754654))
        push!(pts,_sph_pt(-135, 54.735610317245346))
        push!(pts,_sph_pt(-135, 125.264389682754654))
        return pts
    elseif n<=26
        push!(pts,_sph_pt(0,90))
        push!(pts,_sph_pt(180,90))
        push!(pts,_sph_pt(90,90))
        push!(pts,_sph_pt(-90,90))
        push!(pts,_sph_pt(90,0))
        push!(pts,_sph_pt(90,180))
        push!(pts,_sph_pt(90,45))
        push!(pts,_sph_pt(90,135))
        push!(pts,_sph_pt(-90,45))
        push!(pts,_sph_pt(-90,135))
        push!(pts,_sph_pt(0,45))
        push!(pts,_sph_pt(0,135))
        push!(pts,_sph_pt(180,45))
        push!(pts,_sph_pt(180,135))
        push!(pts,_sph_pt(45,90))
        push!(pts,_sph_pt(-45,90))
        push!(pts,_sph_pt(135,90))
        push!(pts,_sph_pt(-135,90))
        push!(pts,_sph_pt(45,54.735610317245346))
        push!(pts,_sph_pt(45,125.264389682754654))
        push!(pts,_sph_pt(-45,54.735610317245346))
        push!(pts,_sph_pt(-45,125.264389682754654))
        push!(pts,_sph_pt(135,54.735610317245346))
        push!(pts,_sph_pt(135,125.264389682754654))
        push!(pts,_sph_pt(-135,54.735610317245346))
        push!(pts,_sph_pt(-135,125.264389682754654))
        return pts
    elseif n<=38
        push!(pts,  _sph_pt(0.000000000000000,    90.000000000000000))
        push!(pts,_sph_pt( 180.000000000000000,    90.000000000000000))
        push!(pts,_sph_pt(  90.000000000000000,    90.000000000000000))
        push!(pts,_sph_pt( -90.000000000000000,    90.000000000000000))
        push!(pts,_sph_pt(  90.000000000000000,     0.000000000000000))
        push!(pts,_sph_pt(  90.000000000000000,   180.000000000000000))
        push!(pts,_sph_pt(  45.000000000000000,    54.735610317245346))
        push!(pts,_sph_pt(  45.000000000000000,   125.264389682754654))
        push!(pts,_sph_pt( -45.000000000000000,    54.735610317245346))
        push!(pts,_sph_pt( -45.000000000000000,   125.264389682754654))
        push!(pts,_sph_pt( 135.000000000000000,    54.735610317245346))
        push!(pts,_sph_pt( 135.000000000000000,   125.264389682754654))
        push!(pts,_sph_pt(-135.000000000000000,    54.735610317245346))
        push!(pts,_sph_pt(-135.000000000000000,   125.264389682754654))
        push!(pts,_sph_pt(  62.632194841377327,    90.000000000000000))
        push!(pts,_sph_pt( -62.632194841377327,    90.000000000000000))
        push!(pts,_sph_pt( 117.367805158622687,    90.000000000000000))
        push!(pts,_sph_pt(-117.367805158622687,    90.000000000000000))
        push!(pts,_sph_pt(  27.367805158622673,    90.000000000000000))
        push!(pts,_sph_pt( -27.367805158622673,    90.000000000000000))
        push!(pts,_sph_pt( 152.632194841377355,    90.000000000000000))
        push!(pts,_sph_pt(-152.632194841377355,    90.000000000000000))
        push!(pts,_sph_pt(   0.000000000000000,    27.367805158622673))
        push!(pts,_sph_pt(   0.000000000000000,   152.632194841377355))
        push!(pts,_sph_pt( 180.000000000000000,    27.367805158622673))
        push!(pts,_sph_pt( 180.000000000000000,   152.632194841377355))
        push!(pts,_sph_pt(   0.000000000000000,    62.632194841377327))
        push!(pts,_sph_pt(   0.000000000000000,   117.367805158622687))
        push!(pts,_sph_pt( 180.000000000000000,    62.632194841377327))
        push!(pts,_sph_pt( 180.000000000000000,   117.367805158622687))
        push!(pts,_sph_pt(  90.000000000000000,    27.367805158622673))
        push!(pts,_sph_pt(  90.000000000000000,   152.632194841377355))
        push!(pts,_sph_pt( -90.000000000000000,    27.367805158622673))
        push!(pts,_sph_pt( -90.000000000000000,   152.632194841377355))
        push!(pts,_sph_pt(  90.000000000000000,    62.632194841377327))
        push!(pts,_sph_pt(  90.000000000000000,   117.367805158622687))
        push!(pts,_sph_pt( -90.000000000000000,    62.632194841377327))
        push!(pts,_sph_pt( -90.000000000000000,   117.367805158622687))
        return pts
    elseif n<=50
        push!(pts,_sph_pt(  0.000000000000000,    90.000000000000000))
        push!(pts,_sph_pt( 180.000000000000000,    90.000000000000000))
        push!(pts,_sph_pt(  90.000000000000000,    90.000000000000000))
        push!(pts,_sph_pt( -90.000000000000000,    90.000000000000000))
        push!(pts,_sph_pt(  90.000000000000000,     0.000000000000000))
        push!(pts,_sph_pt(  90.000000000000000,   180.000000000000000))
        push!(pts,_sph_pt(  90.000000000000000,    45.000000000000000))
        push!(pts,_sph_pt(  90.000000000000000,   135.000000000000000))
        push!(pts,_sph_pt( -90.000000000000000,    45.000000000000000))
        push!(pts,_sph_pt( -90.000000000000000,   135.000000000000000))
        push!(pts,_sph_pt(   0.000000000000000,    45.000000000000000))
        push!(pts,_sph_pt(   0.000000000000000,   135.000000000000000))
        push!(pts,_sph_pt( 180.000000000000000,    45.000000000000000))
        push!(pts,_sph_pt( 180.000000000000000,   135.000000000000000))
        push!(pts,_sph_pt(  45.000000000000000,    90.000000000000000))
        push!(pts,_sph_pt( -45.000000000000000,    90.000000000000000))
        push!(pts,_sph_pt( 135.000000000000000,    90.000000000000000))
        push!(pts,_sph_pt(-135.000000000000000,    90.000000000000000))
        push!(pts,_sph_pt(  45.000000000000000,    54.735610317245346))
        push!(pts,_sph_pt(  45.000000000000000,   125.264389682754654))
        push!(pts,_sph_pt( -45.000000000000000,    54.735610317245346))
        push!(pts,_sph_pt( -45.000000000000000,   125.264389682754654))
        push!(pts,_sph_pt( 135.000000000000000,    54.735610317245346))
        push!(pts,_sph_pt( 135.000000000000000,   125.264389682754654))
        push!(pts,_sph_pt(-135.000000000000000,    54.735610317245346))
        push!(pts,_sph_pt(-135.000000000000000,   125.264389682754654))
        push!(pts,_sph_pt(  45.000000000000000,    25.239401820678911))
        push!(pts,_sph_pt(  45.000000000000000,   154.760598179321079))
        push!(pts,_sph_pt( -45.000000000000000,    25.239401820678911))
        push!(pts,_sph_pt( -45.000000000000000,   154.760598179321079))
        push!(pts,_sph_pt( 135.000000000000000,    25.239401820678911))
        push!(pts,_sph_pt( 135.000000000000000,   154.760598179321079))
        push!(pts,_sph_pt(-135.000000000000000,    25.239401820678911))
        push!(pts,_sph_pt(-135.000000000000000,   154.760598179321079))
        push!(pts,_sph_pt(  71.565051177077990,    72.451599386207704))
        push!(pts,_sph_pt( -71.565051177077990,    72.451599386207704))
        push!(pts,_sph_pt(  71.565051177077990,   107.548400613792296))
        push!(pts,_sph_pt( -71.565051177077990,   107.548400613792296))
        push!(pts,_sph_pt( 108.434948822922010,    72.451599386207704))
        push!(pts,_sph_pt(-108.434948822922010,    72.451599386207704))
        push!(pts,_sph_pt( 108.434948822922010,   107.548400613792296))
        push!(pts,_sph_pt(-108.434948822922010,   107.548400613792296))
        push!(pts,_sph_pt(  18.434948822922017,    72.451599386207704))
        push!(pts,_sph_pt( 161.565051177078004,    72.451599386207704))
        push!(pts,_sph_pt(  18.434948822922017,   107.548400613792296))
        push!(pts,_sph_pt( 161.565051177078004,   107.548400613792296))
        push!(pts,_sph_pt( -18.434948822922017,    72.451599386207704))
        push!(pts,_sph_pt(-161.565051177078004,    72.451599386207704))
        push!(pts,_sph_pt( -18.434948822922017,   107.548400613792296))
        push!(pts,_sph_pt(-161.565051177078004,   107.548400613792296))
        return pts
    else
        @error "Unable to return $n Lebedev points"
    end
end

function _sph_pt(theta,phi,r=1, center=[0 0 0])
    # convert to radians
    theta = theta * pi/180
    phi   = phi   * pi/180
    # perform translation and dilation
    x = center[1] + r*sin(phi)*cos(theta)
    y = center[2] + r*sin(phi)*sin(theta)
    z = center[3] + r*cos(phi)
    return SVector{3,Float64}(x,y,z)
end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
using Test

include("_dummy_model.jl")
using .DummyModel

@testset "utilities" begin
    include("utilities.jl")
end

@testset "ic_model" begin
    include("ic_model.jl")
end

@testset "controls" begin
    include("controls.jl")
end

@testset "constructors" begin
    include("constructors.jl")
end

@testset "core" begin
    include("core.jl")
end

@testset "traits" begin
    include("traits.jl")
end


text:
#GracePlot code generation tools
#-------------------------------------------------------------------------------

#==Main data structures
===============================================================================#

#A type whose elements simply list valid attributes of a given object:
abstract type AttributeList end

#==Helper functions
===============================================================================#

#==Object builders
===============================================================================#

#Generate Expr for an "AttributeList" object constructor function.
#NOTE:
#   -Constructor function supports optional parameters.
#   -Assumes optional parameters can be "nothing".
#Inputs:
#   fnname: Object constructor name
#   t: data type to build
#   reqfieldcnt: # of required field/arguments (must be first object fields)
#TODO: Find way to restrict t::DataTypes <: AttributeList *without* wasting
#      space generating duplicate code
function genexpr_attriblistbuilder(fnname::Symbol, t::DataType; reqfieldcnt::Int=0)
	fnames = fieldnames(t)
	ftypes = [fieldtype(t,n) for n in fnames]

	#Build list of required parameters
	reqlist = Expr[]
	for i in 1:reqfieldcnt
		(name, _type) = (fnames[i], ftypes[i])
		push!(reqlist, :($name::$_type))
	end

	#Build list of optional parameters
	optlist = Expr[]
	for i in reqfieldcnt+1:length(fnames)
		(name, _type) = (fnames[i], ftypes[i])
		#Ignore argument types for now (assumes ::Any)...
		push!(optlist, Expr(:kw, name, nothing)) #Keyword assignment
	end

	#Build entire list of constructed object fields
	constlist = Symbol[]
	for i in 1:length(fnames)
		(name, _type) = (fnames[i], ftypes[i])
		push!(constlist, :($name))
	end

	constructorcall = Expr(:call, t, constlist...)
	return :($fnname($(reqlist...); $(optlist...)) = $constructorcall)
end

#Last line


text:
###julia --machine-file <host-file> distributed.jl
###julia -p <number-of-processors> distributed.jl

# using Distributed # not necessary when doing `julia -p ...` or `julia --machine-file...`
@everywhere begin
    import Pkg
    Pkg.activate("../")
    # Pkg.instantiate()
end
@everywhere begin
    using Apophis
    using Dates
    using TaylorSeries

    #script parameters (TODO: use ArgParse.jl instead)
    const objname = "Apophis"
    const maxsteps = 10000
    const nyears = 5.0
    const dense = true#false
    const apophisjlpath = pkgdir(Apophis)
    const radarobsfile = joinpath(apophisjlpath, "Apophis_JPL_data_2012_2013.dat")
    const dynamics = RNp1BP_pN_A_J23E_J2S_ng_eph!
    const t0 = datetime2julian(DateTime(2008,9,24,0,0,0)) #starting time of integration
    const tmax = t0+365.25nyears #final time of integration
    @show t0 == 2454733.5
    @show tmax

    function earth_et(et)
        return ss16asteph( Apophis.etsecs2julian(et) )[union(3*4-2:3*4,3*(27+4)-2:3*(27+4))]
    end
    function sun_et(et)
        return ss16asteph( Apophis.etsecs2julian(et) )[union(3*1-2:3*1,3*(27+1)-2:3*(27+1))]
    end
end

# path to local Solar System ephemeris file
# ss_eph_file = joinpath(pkgdir(Apophis), "jldeph", "ss16ast343_eph_24yr_tx.jld")
ss_eph_file = joinpath(pkgdir(Apophis), "jldeph", "ss16ast343_eph_5yr_tx.jld")

ss16asteph, acc_eph, newtonianNb_Potential = Apophis.loadeph(ss_eph_file)

aux = (ss16asteph, acc_eph, newtonianNb_Potential, earth_et, sun_et)
for i in 1:nworkers()
    @spawnat i+1 aux = (ss16asteph, acc_eph, newtonianNb_Potential, earth_et, sun_et)
end

#warmup (compilation) short run on all processes
parallel_run(objname, dynamics, 1, t0, tmax, aux, output=false)
println("*** Finished warmup")

#Full jet transport integration until ~2038: about 8,000 steps
# parallel_run(objname, dynamics, maxsteps, t0, tmax, aux, radarobsfile=radarobsfile)
# println("*** Finished full jet transport integration")


text:
"""
    propagate(veh, action, roadway, Δt)

Take an entity of type {S,D,I} and move it over Δt seconds to produce a new
entity based on the action on the given roadway.
"""
propagate(veh::Entity{S,D,I}, action::A, roadway::R, Δt::Float64) where {S,D,I,A,R} = error("propagate not implemented for Entity{$S, $D, $I}, actions $A, and roadway $R")

propagate(veh::Entity{S,D,I}, state::S, roadway::R, Δt::Float64) where {S,D,I,R}= state

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


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

"""
    FIST(shuffle=true)

Fast Industrial-Strength Triangulation (FIST) of polygons.

This triangulation method is the method behind the famous Mapbox's
Earcut library. It is based on a ear clipping algorithm adapted
for complex n-gons with holes. It has O(n²) time complexity where
n is the number of vertices. In practice it is very efficient due
to heuristics implemented in the algorithm.

The option `shuffle` is used to shuffle the order in which ears
are clipped. It improves the quality of the triangles, which can
be very sliver otherwise.

## References

* Held, M. 1998. [FIST: Fast Industrial-Strength Triangulation of Polygons]
  (https://link.springer.com/article/10.1007/s00453-001-0028-4)
* Eder et al. 2018. [Parallelized ear clipping for the triangulation and
  constrained Delaunay triangulation of polygons]
  (https://www.sciencedirect.com/science/article/pii/S092577211830004X)
"""
struct FIST <: DiscretizationMethod
  shuffle::Bool
end

FIST() = FIST(true)

function discretize(𝒫::Chain, method::FIST)
  # points of resulting mesh
  points = vertices(𝒫)

  # keep track of global indices
  inds = CircularVector(1:nvertices(𝒫))

  # perform ear clipping
  𝒬 = ears(𝒫); method.shuffle && shuffle!(𝒬)
  n = nvertices(𝒫)
  𝒯 = Connectivity{Triangle,3}[]
  clipped = false
  while n > 3
    if !isempty(𝒬) # clip an ear
      # 0. select candidate ear
      i = pop!(𝒬); 𝒬[𝒬.>i] .-= 1
      # 1. push a new triangle to 𝒯
      push!(𝒯, connect((inds[i-1], inds[i], inds[i+1]), Triangle))
      # 2. remove the vertex from 𝒫
      inds = inds[setdiff(1:n, mod1(i,n))]
      𝒫 = Chain(points[inds])
      n = nvertices(𝒫)
      # 3. update 𝒬 near clipped ear
      for j in (i-1, i)
        if isear(𝒫, j)
          𝒬 = 𝒬 ∪ [mod1(j,n)]
        else
          setdiff!(𝒬, [mod1(j,n)])
        end
      end
      clipped = true
    elseif clipped # recompute all ears
      𝒬 = ears(𝒫); method.shuffle && shuffle!(𝒬)
      clipped = false
    else # recovery process
      # check if consecutive edges vᵢ-1 -- vᵢ and vᵢ+1 -- vᵢ+2
      # intersect and fix the issue by clipping ear (vᵢ, vᵢ+1, vᵢ+2)
      v = vertices(𝒫)
      for i in 1:n
        s1 = Segment(v[i-1], v[i])
        s2 = Segment(v[i+1], v[i+2])
        if intersecttype(s1, s2) isa CrossingSegments
          # 1. push a new triangle to 𝒯
          push!(𝒯, connect((inds[i], inds[i+1], inds[i+2]), Triangle))
          # 2. remove the vertex from 𝒫
          inds = inds[setdiff(1:n, mod1(i+1,n))]
          𝒫 = Chain(points[inds])
          n = nvertices(𝒫)
          clipped = true
          break
        end
      end
    end
  end
  # remaining polygonal area is the last triangle
  push!(𝒯, connect((inds[1], inds[2], inds[3]), Triangle))

  SimpleMesh(collect(points), 𝒯)
end

# return index of all ears of 𝒫
ears(𝒫) = filter(i -> isear(𝒫, i), 1:nvertices(𝒫))

# tell whether or not vertex i is an ear of 𝒫
function isear(𝒫, i)
  O = orientation(𝒫, TriangleOrientation())
  if O == :CCW
    isearccw(𝒫, i)
  else
    # reverse chain and index
    n = nvertices(𝒫)
    ℛ = reverse(𝒫)
    j = n - i - 1
    isearccw(ℛ, j)
  end
end

# tells whether or not vertex i is an ear of 𝒫
# assuming that 𝒫 has counter-clockwise orientation
function isearccw(𝒫::Chain{Dim,T}, i) where {Dim,T}
  v = vertices(𝒫)

  # helper function to compute the vexity of vertex i
  function vexity(i)
    α = ∠(v[i-1], v[i], v[i+1]) # oriented angle
    θ = α > 0 ? 2*T(π) - α : -α # inner angle
    θ < π ? :CONVEX : :REFLEX
  end

  # helper function to check if vertex j is inside cone i
  function incone(j, i)
    s1 = sideof(v[j], Segment(v[i], v[i-1]))
    s2 = sideof(v[j], Segment(v[i], v[i+1]))
    if vexity(i) == :CONVEX
      s1 != :LEFT && s2 != :RIGHT
    else
      s1 != :LEFT || s2 != :RIGHT
    end
  end

  # CE1.1: classify angle as convex vs. reflex
  isconvex = vexity(i) == :CONVEX

  # CE1.2: check if segment vᵢ-₁ -- vᵢ+₁ intersects 𝒫
  sᵢ = Segment(v[i-1], v[i+1])
  intersects = false
  for j in 1:nvertices(𝒫)
    sⱼ = Segment(v[j], v[j+1])
    I = intersecttype(sᵢ, sⱼ)
    if !(I isa CornerTouchingSegments || I isa NoIntersection)
      intersects = true
      break
    end
  end

  # CE1.3: check if vᵢ-1 ∈ C(vᵢ, vᵢ+1, vᵢ+2) and vᵢ+1 ∈ C(vᵢ-2, vᵢ-1, vᵢ)
  incones = incone(i-1, i+1) && incone(i+1, i-1)

  isconvex && !intersects && incones
end


text:
include("../src/cel_bulirsch.jl")
include("../src/transit_structure.jl")
include("../src/transit_poly_struct.jl")
include("../src/Mn_compute.jl")
using QuadGK

function Mn_num(k2::T,m::T) where {T <: Real}
# Numerically integrates M_n(k^2)/(4br)^m (note: m can be a half-integer).
# See 11/10/2018 notes.
f(x) = (k2-sin(x)^2)^m
g(x) = (k2-sin(x)^2)^m*sin(x)^2
if k2 < 1.0
  kap2 = convert(T,asin(sqrt(big(k2))))
  Mn,error = quadgk(f,-kap2,kap2,rtol=1e-15)
  Nn,error = quadgk(g,-kap2,kap2,rtol=1e-15)
else
  pi2 = 0.5*pi
  Mn,error = quadgk(f,-pi2,pi2,rtol=1e-15)
  Nn,error = quadgk(g,-pi2,pi2,rtol=1e-15)
end
return Mn,Nn
end

n = 30
t= transit_init(0.1,0.5,ones(n)/n,false)
nb = 7
rgrid = [0.01,0.1,1.0,10.0,100.0]
for r in rgrid
Mnr = zeros(t.n_max+1)
Nnr = zeros(t.n_max+1)
Mnl = zeros(t.n_max+1)
Nnl = zeros(t.n_max+1)
for i=1:nb
#r = 10.^(rand()*3-2); ; b=rand()*((1+r)-b0)+b0; fourbr = 4*b*r
  b0 = maximum([1e-8,r-1])
  b = b0+i/nb*(1+r-b0)
  t.r = r; t.b = b; t.fourbr = 4*b*r; t.sqbr = sqrt(b*r)
  t.sqarea = sqarea_triangle(1.0,r,b)
  t.onembmr2 = (1.0-(r-b)^2); t.onembpr2 = 1.0-(b+r)^2
  t.sqonembmr2 = sqrt(t.onembmr2)
  t.k2 = t.onembmr2/t.fourbr;
  t.onemr2mb2 = 1-r^2-b^2
  if t.k2 < 1.0
    println("k2: ",t.k2)
    t.k2inv=inv(t.k2); t.kc = sqrt(1.0-t.k2); t.k = sqrt(t.k2); t.kap0 = 2*asin(t.k)
    t.Eofk = cel_bulirsch(t.k2,t.kc,1.0,1.0,1.0-t.k2)
    t.Em1mKdm = cel_bulirsch(t.k2,t.kc,1.0,1.0,0.0)
    t.kite_area2 = sqrt(abs(t.sqarea))
  else
    println("k2: ",t.k2)
    t.k2inv=inv(t.k2); t.kc = sqrt(1.0-t.k2inv); t.k = sqrt(t.k2); t.kap0 = pi
    t.Eofk = cel_bulirsch(t.k2inv,t.kc,1.0,1.0,1.0-t.k2inv)
    t.Em1mKdm = cel_bulirsch(t.k2inv,t.kc,1.0,1.0,0.0)
  end
  Mn_raise!(t)
  Mnr .= t.Mn
  Nn_raise!(t)
  Nnr .= t.Nn
  Mn_lower!(t)
  Nn_lower!(t)
  Mnl .= t.Mn
  Nnl .= t.Nn
  Mnn = zeros(t.n_max+1)
  # Now, compute the N_n functions: N_n = (4br)^{n/2} \int_{-\kappa/2}^{\kappa/2} (k^2-\sin^2\phi)^{n/2} \sin^2\phi d\phi
  for m=0:t.n_max
    Mn_n,Nn_num = Mn_num(t.k2,float(m/2))
    Mnn[m+1] = Mn_n*sqrt(t.fourbr)^m
    Nn_num *= sqrt(t.fourbr)^m
    println(Mnr[m+1]," ",Mnl[m+1]," ",Mnn[m+1]," ",Mnr[m+1]-Mnn[m+1]," ",Mnl[m+1]-Mnn[m+1]," ",Mnl[m+1]-Mnr[m+1])
    # Now use function which doesn't involve division by b:
    if t.k2 < 0.5
      println("n: ",m," N_n: ",Nnl[m+1]," N_n num: ",Nn_num)
    else
      println("n: ",m," N_n: ",Nnr[m+1]," N_n num: ",Nn_num)
    end
  end
  println("k2: ",t.k2," r: ",r," b: ",b," max(M_n-M_{n,num}): ",maximum(abs,Mnr-Mnn)," ",maximum(abs,Mnl-Mnn))
  read(STDIN,Char)
end
end


text:
using JSON

function pre()
  # data = JSON.parsefile("training_data1.json")
  data = JSON.parsefile("training_data2.json")

  n = length(data["Labels"])

  labelled_data = Array{Tuple{Array{Float64,1},Int64}}(undef, n)

  for i in 1:n
    labelled_data[i] = ( data["Features"][i], data["Labels"][i] )
  end

  return labelled_data
end

# @time pre()

max_stars_repo_name: dandanua/QML

max_issues_repo_name: dandanua/QML

max_forks_repo_name: dandanua/QML


text:
module PolynomialAsymptotics

using FastGaussQuadrature
using SpecialFunctions
using Roots
using LinearAlgebra

import FastGaussQuadrature: besselroots

include("asy_laguerre.jl")

include("gauss/gausslaguerre.jl")
include("gauss/gaussfreud.jl")
include("gauss/gaussjacobi.jl")

export asy_gausslaguerre
export asy_gaussfreud
export asy_gaussjacobi

end # module


text:
module MySample

export naiveform!, ipform!, opform!, mpform!
export unroll_naive!, unroll_i!, unroll_j!, unroll_k!, unroll_ij!
include("productform.jl")
include("unrolling.jl")

end # module


text:
"""
    merge(left::Array{Int64,1}, right::Array{Int64,1}))

The merge algorithms (MA) merge the subarrays `left` and `right` to produce new sorted
subarrays until there is only one subarray remaining, which will be the sorted array. For
more information see: [https://en.wikipedia.org/wiki/Merge_sort](https://en.wikipedia.org/wiki/Merge_sort)

# Arguments
- `left::Array{Int64,1}`: Unsorted left part of the array
- `right::Array{Int64,1}`: Unsorted right part of the array

"""
function merge(left::Array{Int64,1}, right::Array{Int64,1})
    result = zeros(Int64, false)
    # Left and right may have elements left; consume them.
    while left != [] && right != []
        # Check if left is smaller than right, if yes append
        if left[1] <= right[1]
            append!(result, left[1])
            left = left[2:length(left)]
        else
            # If not mix the order of left and right
            append!(result, right[1])
            right = right[2:length(right)]
        end
    end

    # Either left or right may have elements left; consume them.
    # (Only one of the following loops will actually be entered.)
    while left != []
        append!(result, left[1])
        left = left[2:length(left)]
    end
    while right != []
        append!(result, right[1])
        right = right[2:length(right)]
    end
    return result
end

"""
    merge_sorting(array::Array{Int64,1})

The merge sort algorithms (MSA) are a comparison-based sorting algorithm, which is referred
to as the divide and conquer algorithms.  The stable sort implementation is a widely used
method for the MSA, which means that the order of equal elements is the same in the input
and output. In the current implementation, a top-down implementation is used; however, a
Bottom-up implementation can be used, too. In the top-down implementation, the MSA
recursively splits the array into subarrays until the subarray size is < 2, merging those
subarrays to produce a sorted array by using a new function `merge`. The back copying is
blocked by alternating the direction of the merge with each recursion. For more information
see: [https://en.wikipedia.org/wiki/Merge_sort](https://en.wikipedia.org/wiki/Merge_sort)

# Arguments
- `array::Array{Int64,1}`: Unsorted array of integers

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr = [64, 34, 25, 12, 22, 11, 90]
julia> ClassicAlgorithmsCollections.merge_sorting(arr)
[11, 12, 22, 25, 34, 64, 90]
```
"""
function merge_sorting(array::Array{Int64,1})
    # Base case. A array of zero or one elements is sorted, by definition.
    n = length(array)
    if n < 2
        return array
    end

    # First, divide the array into equal-sized subarrays
    mid = fld(n, 2)
    left = array[1:mid]
    right = array[mid+1:n]
    # Recursively sorting the left and right half
    left = merge_sorting(left)
    right = merge_sorting(right)
    # Then merge the now-sorted left and right subarrays.
    array = merge(left, right)
    # Finally return the sorted array
    return array
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
using Echidna
# using Plots
using Test

# gr()

println("Setting up ZDT1 test problem")
nobjs = 2
nvars = 30
zdt1_problem = Problem(
    nobjs, [false for i in 1:nobjs], nvars,
    [MOGA_Real(0.0, 1.0) for i in 1:30],
    ZDT1
)

println("Setting up NSGAIII config")
niters = 200
ndivs = 25
# refpts = generate_regular_reference_points(2, ndivs)
algo = NSGAIII(zdt1_problem, ZDT1, ndivs, niters)
result = garun(algo)

@show result

# rx = [s.objectives[1] for s in result]
# ry = [s.objectives[2] for s in result]
# scatter(rx, ry)

# seedgenomes = [(n/10.0)*ones(nvars) for n in 1:9]
# seedpop = [candidate_from_genome(zdt1_problem, g) for g in seedgenomes]

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna