text:
# Make arrays of field behave like Maps

function return_cache(f::AbstractArray{T},x::Point) where T<:Field
  S = return_type(testitem(f),x)
  cr = CachedArray(zeros(S,size(f)))
  if isconcretetype(T)
    cf = return_cache(testitem(f),x)
  else
    cf = nothing
  end
  cr, cf
end

"""
Implementation of `return_cache` for a array of `Field`.

If the field vector has length `nf` and it is evaluated in one point, it
returns an `nf` vector with the result. If the same array is applied to a
vector of `np` points, it returns a matrix `np` x `nf`.
"""
function evaluate!(c,f::AbstractArray{T},x::Point) where T<:Field
  cr, cf = c
  setsize!(cr,size(f))
  r = cr.array
  if isconcretetype(T)
    for j in eachindex(f)
      @inbounds r[j] = evaluate!(cf,f[j],x)
    end
  else
    for j in eachindex(f)
      @inbounds r[j] = evaluate(f[j],x)
    end
  end
  r
end

function return_cache(f::AbstractArray{T},x::AbstractArray{<:Point}) where T<:Field
  S = return_type(testitem(f),testitem(x))
  cr = CachedArray(zeros(S,(size(x)...,size(f)...)))
  if isconcretetype(T)
    cf = return_cache(f,testitem(x))
  else
    cf = nothing
  end
  cr, cf
end

function evaluate!(c,f::AbstractArray{T},x::AbstractArray{<:Point}) where T<:Field
  cr, cf = c
  setsize!(cr,(size(x)...,size(f)...))
  r = cr.array
  if isconcretetype(T)
    for i in eachindex(x)
      fxi = evaluate!(cf,f,x[i])
      for j in CartesianIndices(f)
        r[i,j] = fxi[j]
      end
    end
  else
    for j in eachindex(f)
      for i in eachindex(x)
        r[i,j] = evaluate(f[j],x[i])
      end
    end
  end
  r
end

function testargs(f::AbstractArray{T},x::Point) where T<:Field
  testargs(testitem(f),x)
end

function testargs(f::AbstractArray{T},x::AbstractArray{<:Point}) where T<:Field
  testargs(testitem(f),x)
end

function test_field_array(f::AbstractArray{<:Field}, x, v, cmp=(==); grad=nothing, gradgrad=nothing)
  test_map(v,f,x;cmp=cmp)
  if grad != nothing
    test_map(grad,Broadcasting(∇)(f),x;cmp=cmp)
  end
  if gradgrad != nothing
    test_map(gradgrad,Broadcasting(∇∇)(f),x;cmp=cmp)
  end
  true
end

# Opening the door to optimize arrays of field gradients

"""
A wrapper that represents the broadcast of `gradient` over an array of fields.
Ng is the number of times the gradient is applyed
"""
struct FieldGradientArray{Ng,A,T,N} <: AbstractArray{T,N}
  fa::A
  function FieldGradientArray{Ng}(f::AbstractArray{<:Field}) where Ng
    T = typeof(gradient(testitem(f),Val(Ng)))
    N = ndims(f)
    A = typeof(f)
    new{Ng,A,T,N}(f)
  end
end

function return_value(k::Broadcasting{typeof(∇)},a::AbstractArray{<:Field})
  evaluate(k,a)
end

function return_value(k::Broadcasting{typeof(∇∇)},a::AbstractArray{<:Field})
  evaluate(k,a)
end

function evaluate!(cache,k::Broadcasting{typeof(∇)},a::AbstractArray{<:Field})
  FieldGradientArray{1}(a)
end

function evaluate!(cache,k::Broadcasting{typeof(∇)},a::FieldGradientArray{N}) where N
  FieldGradientArray{N+1}(a.fa)
end

function evaluate!(cache,k::Broadcasting{typeof(∇∇)},a::AbstractArray{<:Field})
  FieldGradientArray{2}(a)
end

function gradient(a::AbstractArray{<:Field})
  msg =
  """\n
  Function gradient (aka ∇) is not defined for arrays of Field objects.
  Use Broadcasting(∇) instead.
  """
  @unreachable msg
end

function ∇∇(a::AbstractArray{<:Field})
  msg =
  """\n
  Double gradient application (aka ∇∇) is not defined for arrays of Field objects.
  Use Broadcasting(∇∇) instead.
  """
  @unreachable msg
end

Base.size(a::FieldGradientArray) = size(a.fa)
Base.axes(a::FieldGradientArray) = axes(a.fa)
Base.getindex(a::FieldGradientArray{Ng},i::Integer) where Ng = gradient(a.fa[i],Val(Ng))
Base.getindex(
   a::FieldGradientArray{Ng,A,T,N},i::Vararg{Integer,N}) where {Ng,A,T,N} = gradient(a.fa[i...],Val(Ng))
Base.IndexStyle(::Type{<:FieldGradientArray{Ng,A}}) where {Ng,A} = IndexStyle(A)

# Optimizing linear_combination.

function linear_combination(a::AbstractVector{<:Number},b::AbstractVector{<:Field})
  column = 1
  LinearCombinationField(a,b,column)
end

struct LinearCombinationField{V,F} <: Field
  values::V
  fields::F
  column::Int
end

for T in (:(Point),:(AbstractVector{<:Point}))
  @eval begin

  function return_cache(a::LinearCombinationField,x::$T)
    cf = return_cache(a.fields,x)
    fx = return_value(a.fields,x)
    v = a.values
    k = LinearCombinationMap(a.column)
    ck = return_cache(k,v,fx)
    cf, ck
  end

  function evaluate!(cache,a::LinearCombinationField,x::$T)
    cf, ck = cache
    fx = evaluate!(cf,a.fields,x)
    v = a.values
    k = LinearCombinationMap(a.column)
    evaluate!(ck,k,v,fx)
  end

  end
end

for op in (:∇,:∇∇)
  @eval begin
    function $op(a::LinearCombinationField)
      fields = Broadcasting($op)(a.fields)
      LinearCombinationField(a.values,fields,a.column)
    end
  end
end

function linear_combination(a::AbstractMatrix{<:Number},b::AbstractVector{<:Field})
  #[ LinearCombinationField(a,b,i) for i in 1:size(a,2) ]
  LinearCombinationFieldVector(a,b)
end

struct LinearCombinationFieldVector{V,F} <: AbstractVector{LinearCombinationField{V,F}}
  values::V
  fields::F
  function LinearCombinationFieldVector(values::AbstractMatrix{<:Number},fields::AbstractVector{<:Field})
    @check size(values,1) == length(fields) """\n
    Incompatible sizes for performing the linear combination

        linear_combination(values,fields) = transpose(values)*fields

    size(values,1) != length(fields)
    """
    V = typeof(values)
    F = typeof(fields)
    new{V,F}(values,fields)
  end
end

Base.size(a::LinearCombinationFieldVector) = (size(a.values,2),)
Base.getindex(a::LinearCombinationFieldVector,i::Integer) = LinearCombinationField(a.values,a.fields,i)
Base.IndexStyle(::Type{<:LinearCombinationField}) = IndexLinear()

for T in (:(Point),:(AbstractVector{<:Point}))
  @eval begin

    function return_cache(a::LinearCombinationFieldVector,x::$T)
      cf = return_cache(a.fields,x)
      fx = return_value(a.fields,x)
      v = a.values
      k = LinearCombinationMap(:)
      ck = return_cache(k,v,fx)
      cf, ck
    end

    function evaluate!(cache,a::LinearCombinationFieldVector,x::$T)
      cf, ck = cache
      fx = evaluate!(cf,a.fields,x)
      v = a.values
      k = LinearCombinationMap(:)
      evaluate!(ck,k,v,fx)
    end

  end
end

for op in (:∇,:∇∇)
  @eval begin
    function evaluate!(cache,k::Broadcasting{typeof($op)},a::LinearCombinationFieldVector)
      fields = Broadcasting($op)(a.fields)
      LinearCombinationFieldVector(a.values,fields)
    end
  end
end

function get_children(n::TreeNode, a::LinearCombinationFieldVector)
  (similar_tree_node(n,a.values),similar_tree_node(n,a.fields))
end

# This is the map that acts on values
struct LinearCombinationMap{T} <: Map
  column::T
  LinearCombinationMap(column::Integer)  = new{typeof(column)}(column)
  LinearCombinationMap(column::Colon)  = new{typeof(column)}(column)
end

function evaluate!(cache,k::LinearCombinationMap{<:Integer},v::AbstractArray,fx::AbstractVector)
  z = zero(return_type(outer,testitem(fx),testitem(v)))
  @check length(fx) == size(v,1)
  @inbounds for i in eachindex(fx)
    # We need to do the product in this way
    # so that the gradient also works
    z += outer(fx[i],v[i,k.column])
  end
  z
end

function return_value(k::LinearCombinationMap{<:Integer},v::AbstractArray,fx::AbstractMatrix)
  if size(fx,2) == size(v,1)
    evaluate(k,v,fx)
  else
    c = return_cache(k,v,fx)
    c.array
  end
end

function return_cache(k::LinearCombinationMap{<:Integer},v::AbstractArray,fx::AbstractMatrix)
  vf = testitem(fx)
  vv = testitem(v)
  T = typeof( vf⊗vv + vf⊗vv )
  r = zeros(T,size(fx,1))
  CachedArray(r)
end

function evaluate!(cache,k::LinearCombinationMap{<:Integer},v::AbstractArray,fx::AbstractMatrix)
  @check size(fx,2) == size(v,1)
  setsize!(cache,(size(fx,1),))
  r = cache.array
  z = zero(eltype(r))
  @inbounds for p in 1:size(fx,1)
    rp = z
    for i in 1:size(fx,2)
      rp += outer(fx[p,i],v[i,k.column])
    end
    r[p] = rp
  end
  r
end

function evaluate!(cache,k::LinearCombinationMap{Colon},v::AbstractVector,fx::AbstractVector)
  evaluate!(cache,LinearCombinationMap(1),v,fx)
end

function return_value(k::LinearCombinationMap{Colon},v::AbstractVector,fx::AbstractMatrix)
  return_value(LinearCombinationMap(1),v,fx)
end

function return_cache(k::LinearCombinationMap{Colon},v::AbstractVector,fx::AbstractMatrix)
  return_cache(LinearCombinationMap(1),v,fx)
end

function evaluate!(cache,k::LinearCombinationMap{Colon},v::AbstractVector,fx::AbstractMatrix)
  evaluate!(cache,LinearCombinationMap(1),v,fx)
end

function return_cache(k::LinearCombinationMap{Colon},v::AbstractMatrix,fx::AbstractVector)
  vf = testitem(fx)
  vv = testitem(v)
  T = typeof( vf⊗vv + vf⊗vv )
  r = zeros(T,size(v,2))
  CachedArray(r)
end

function evaluate!(cache,k::LinearCombinationMap{Colon},v::AbstractMatrix,fx::AbstractVector)
  @check length(fx) == size(v,1)
  setsize!(cache,(size(v,2),))
  r = cache.array
  @inbounds for j in eachindex(r)
    rj = zero(eltype(r))
    for i in eachindex(fx)
      rj += outer(fx[i],v[i,j])
    end
    r[j] = rj
  end
  r
end

function return_cache(k::LinearCombinationMap{Colon},v::AbstractMatrix,fx::AbstractMatrix)
  vf = testitem(fx)
  vv = testitem(v)
  T = typeof( vf⊗vv + vf⊗vv )
  r = zeros(T,size(fx,1),size(v,2))
  CachedArray(r)
end

function evaluate!(cache,k::LinearCombinationMap{Colon},v::AbstractMatrix,fx::AbstractMatrix)
  @check size(fx,2) == size(v,1)
  setsize!(cache,(size(fx,1),size(v,2)))
  r = cache.array
  @inbounds for p in 1:size(fx,1)
    for j in 1:size(r,2)
      rj = zero(eltype(r))
      for i in 1:size(fx,2)
        rj += outer(fx[p,i],v[i,j])
      end
      r[p,j] = rj
    end
  end
  r
end

# Optimizing transpose
testitem(a::Transpose{<:Field}) = testitem(a.parent)
evaluate!(cache,k::Broadcasting{typeof(∇)},a::Transpose{<:Field}) = transpose(k(a.parent))
evaluate!(cache,k::Broadcasting{typeof(∇∇)},a::Transpose{<:Field}) = transpose(k(a.parent))

return_cache(k::Transpose{<:Field},x::Point) = return_cache(k.parent,x)
evaluate!(cache,k::Transpose{<:Field},x::Point) = transpose(evaluate!(cache,k.parent,x))

return_cache(k::Transpose{<:Field},x::AbstractVector{<:Point}) = return_cache(k.parent,x)
function evaluate!(cache,k::Transpose{<:Field},x::AbstractVector{<:Point})
  TransposeFieldIndices(evaluate!(cache,k.parent,x))
end

struct TransposeMap <: Map end
evaluate!(cache,k::TransposeMap,a::AbstractVector) = transpose(a)
evaluate!(cache,k::TransposeMap,a::AbstractMatrix) = TransposeFieldIndices(a)

"""
Given a matrix `np` x `nf1` x `nf2` result of the evaluation of a field vector
on a vector of points, it returns an array in which the field axes (second and
third axes) are permuted. It is equivalent as `Base.permutedims(A,(1,3,2))`
but more performant, since it does not involve allocations.
"""
struct TransposeFieldIndices{A,T} <: AbstractArray{T,3}
  matrix::A
  function TransposeFieldIndices(matrix::AbstractMatrix{T}) where T
    A = typeof(matrix)
    new{A,T}(matrix)
  end
end

function TransposeFieldIndices{A,T}(::UndefInitializer,shape::NTuple{3,Integer}) where {A,T}
  TransposeFieldIndices(similar(A,(shape[1],shape[3])))
end
Base.size(a::TransposeFieldIndices) = (size(a.matrix,1),1,size(a.matrix,2))
Base.axes(a::TransposeFieldIndices) = (axes(a.matrix,1),Base.OneTo(1),axes(a.matrix,2))
Base.IndexStyle(::Type{<:TransposeFieldIndices{A}}) where A = IndexStyle(A)
Base.getindex(a::TransposeFieldIndices,i::Integer,j::Integer,k::Integer) = a.matrix[i,k]
Base.getindex(a::TransposeFieldIndices,i::Integer) = a.matrix[i]
Base.setindex!(a::TransposeFieldIndices,v,i::Integer,j::Integer,k::Integer) = (a.matrix[i,k] = v)
Base.setindex!(a::TransposeFieldIndices,v,i::Integer) = (a.matrix[i] = v)

# Integration

"""
Integration of a given array of fields in the "physical" space
"""
function integrate(a::AbstractArray{<:Field},x::AbstractVector{<:Point},w::AbstractVector{<:Real})
  cache = return_cache(integrate,a,x,w)
  evaluate!(cache,integrate,a,x,w)
end

"""
Integration of a given array of fields in the "reference" space
"""
function integrate(a::AbstractArray{<:Field},q::AbstractVector{<:Point},w::AbstractVector{<:Real},j::Field)
  cache = return_cache(integrate,a,q,w,j)
  evaluate!(cache,integrate,a,q,w,j)
end

# Broadcast operations

function return_value(k::Broadcasting{<:Operation},args::Union{Field,AbstractArray{<:Field}}...)
  BroadcastOpFieldArray(k.f.op,args...)
end

function evaluate!(cache,k::Broadcasting{<:Operation},args::Union{Field,AbstractArray{<:Field}}...)
  BroadcastOpFieldArray(k.f.op,args...)
end

"""
Type that represents a broadcast operation over a set of `AbstractArray{<:Field}`.
The result is a sub-type of `AbstractArray{<:Field}`
"""
struct BroadcastOpFieldArray{O,T,N,A} <: AbstractArray{T,N}
  op::O
  args::A
  function BroadcastOpFieldArray(op,args::Union{Field,AbstractArray{<:Field}}...)
    fs = map(testitem,args)
    T = return_type(Operation(op),fs...)
    s = map(size,args)
    bs = Base.Broadcast.broadcast_shape(s...)
    N = length(bs)
    A = typeof(args)
    O = typeof(op)
    new{O,T,N,A}(op,args)
  end
end

Base.size(a::BroadcastOpFieldArray) = Base.Broadcast.broadcast_shape(map(size,a.args)...)
Base.axes(a::BroadcastOpFieldArray) = Base.Broadcast.broadcast_shape(map(axes,a.args)...)
Base.IndexStyle(::Type{<:BroadcastOpFieldArray}) = IndexLinear()
Base.getindex(a::BroadcastOpFieldArray,i::Integer) = broadcast(Operation(a.op),a.args...)[i]
function testitem(a::BroadcastOpFieldArray)
  fs = map(testitem,a.args)
  return_value(Operation(a.op),fs...)
end

for T in (:(Point),:(AbstractArray{<:Point}))
  @eval begin

    function return_cache(f::BroadcastOpFieldArray,x::$T)
      cfs = map(fi -> return_cache(fi,x),f.args)
      rs = map(fi -> return_value(fi,x),f.args)
      bm = BroadcastingFieldOpMap(f.op)
      r = return_cache(bm,rs...)
      r, cfs
    end

    function evaluate!(c,f::BroadcastOpFieldArray,x::$T)
      r, cfs = c
      rs = map((ci,fi) -> evaluate!(ci,fi,x),cfs,f.args)
      bm = BroadcastingFieldOpMap(f.op)
      evaluate!(r,bm,rs...)
    end

  end
end

# With this type we mark that we are doing Broadcasting(op) on the result of evaluating Fields/FieldArrays
# This allow us to do some optimizations for block arrays that are only true in this context, not in a
# general Broadcasting operation.
struct BroadcastingFieldOpMap{F} <: Map
  op::F
end

return_value(a::BroadcastingFieldOpMap,args...) = return_value(Broadcasting(a.op),args...)
return_cache(a::BroadcastingFieldOpMap,args...) = return_cache(Broadcasting(a.op),args...)
evaluate!(cache,a::BroadcastingFieldOpMap,args...) = evaluate!(cache,Broadcasting(a.op),args...)

return_value(a::BroadcastingFieldOpMap,args::AbstractArray...) = return_value(Broadcasting(a.op),args...)
return_cache(a::BroadcastingFieldOpMap,args::AbstractArray...) = return_cache(Broadcasting(a.op),args...)
evaluate!(cache,a::BroadcastingFieldOpMap,args::AbstractArray...) = evaluate!(cache,Broadcasting(a.op),args...)

# Follow optimizations are very important to achieve performance

function evaluate!(
  cache,
  f::BroadcastingFieldOpMap,
  a::AbstractArray{T,N},
  b::AbstractArray{S,N}) where {T,S,N}

  @check size(a) == size(b) || (length(a)==0 && length(b)==0)
  setsize!(cache,size(a))
  r = cache.array
  for i in eachindex(a)
    r[i] = f.op(a[i],b[i])
  end
  r
end

function evaluate!(
  cache,
  f::BroadcastingFieldOpMap,
  a::AbstractMatrix,
  b::AbstractArray{S,3} where S)

  @check size(a,1) == size(b,1)
  @check size(b,2) == 1 || size(b,1) == 0
  np, ni = size(a)
  nj = size(b,3)
  setsize!(cache,(np,ni,nj))
  r = cache.array
  for j in 1:nj
    for p in 1:np
      bpj = b[p,1,j]
      for i in 1:ni
        r[p,i,j] = f.op(a[p,i],bpj)
      end
    end
  end
  r
end

function evaluate!(
  cache,
  f::BroadcastingFieldOpMap,
  b::AbstractArray{S,3} where S,
  a::AbstractMatrix)

  @check size(a,1) == size(b,1)
  @check size(b,2) == 1 || size(b,1) == 0
  np, ni = size(a)
  nj = size(b,3)
  setsize!(cache,(np,ni,nj))
  r = cache.array
  for p in 1:np
    for j in 1:nj
      bpj = b[p,1,j]
      for i in 1:ni
        r[p,i,j] = f.op(bpj,a[p,i])
      end
    end
  end
  r
end

function evaluate!(
  cache,
  f::BroadcastingFieldOpMap,
  a::AbstractVector,
  b::AbstractMatrix)

  @check size(a,1) == size(b,1)
  np, ni = size(b)
  setsize!(cache,(np,ni))
  r = cache.array
  for p in 1:np
    ap = a[p]
    for i in 1:ni
      r[p,i] = f.op(ap,b[p,i])
    end
  end
  r
end

function evaluate!(
  cache,
  f::BroadcastingFieldOpMap,
  b::AbstractMatrix,
  a::AbstractVector)

  @check size(a,1) == size(b,1)
  np, ni = size(b)
  setsize!(cache,(np,ni))
  r = cache.array
  for p in 1:np
    ap = a[p]
    for i in 1:ni
      r[p,i] = f.op(b[p,i],ap)
    end
  end
  r
end

function evaluate!(
  cache,
  f::BroadcastingFieldOpMap,
  a::AbstractVector,
  b::AbstractArray{S,3} where S)

  @check size(a,1) == size(b,1)
  np, ni, nj = size(b)
  setsize!(cache,(np,ni,nj))
  r = cache.array
  for p in 1:np
    ap = a[p]
    for j in 1:nj
      for i in 1:ni
        r[p,i,j] = f.op(ap,b[p,i,j])
      end
    end
  end
  r
end

function evaluate!(
  cache,
  f::BroadcastingFieldOpMap,
  b::AbstractArray{S,3} where S,
  a::AbstractVector)

  @check size(a,1) == size(b,1)
  np, ni, nj = size(b)
  setsize!(cache,(np,ni,nj))
  r = cache.array
  for p in 1:np
    ap = a[p]
    for j in 1:nj
      for i in 1:ni
        r[p,i,j] = f.op(b[p,i,j],ap)
      end
    end
  end
  r
end

# Gradient of the sum
for op in (:+,:-)
  @eval begin
    function evaluate!(cache,::Broadcasting{typeof(∇)},a::BroadcastOpFieldArray{typeof($op)})
      f = a.args
      g = map( Broadcasting(∇), f)
      Broadcasting(Operation($op))(g...)
    end
  end
end

# Gradient of the product
for op in (:*,:⋅,:⊙,:⊗)
  @eval begin
    function evaluate!(cache,::Broadcasting{typeof(∇)},a::BroadcastOpFieldArray{typeof($op)})
      f = a.args
      @notimplementedif length(f) != 2
      f1, f2 = f
      g1, g2 = map(Broadcasting(∇), f)
      k(F1,F2,G1,G2) = product_rule($op,F1,F2,G1,G2)
      Broadcasting(Operation(k))(f1,f2,g1,g2)
    end
  end
end


text:
#   This file is part of UnitSystems.jl
#   It is licensed under the MIT license
#   UnitSystems Copyright (C) 2020 Michael Reed
#       _           _                         _
#      | |         | |                       | |
#   ___| |__   __ _| | ___ __ __ ___   ____ _| | __ _
#  / __| '_ \ / _` | |/ / '__/ _` \ \ / / _` | |/ _` |
# | (__| | | | (_| |   <| | | (_| |\ V / (_| | | (_| |
#  \___|_| |_|\__,_|_|\_\_|  \__,_| \_/ \__,_|_|\__,_|
#
#   https://github.com/chakravala
#   https://crucialflow.com

const deka,byte = 𝟐*𝟓,𝟐^3
const hecto, kilo =  deka^2, deka^3
const mega,giga,tera,peta,exa = kilo^2, kilo^3, kilo^4, kilo^5, kilo^6
const deci,centi,milli,micro,nano,pico,femto,atto = inv(deka),inv(hecto),inv(kilo),inv(mega),inv(giga),inv(tera),inv(peta),inv(exa)
const kibi,mebi,gibi,tebi,pebi,exbi,zebi,yobi = 𝟐^10,𝟐^20,𝟐^30,𝟐^40,𝟐^50,𝟐^60,(Constant(1.0)*𝟐)^70,(Constant(1.0)*𝟐)^80

const fur,°R,K,HOUR,k = 𝟐^2*𝟑*𝟓*𝟏𝟏*ft,𝟓/𝟑^2,𝟑^2/𝟓,𝟐^4*𝟑^2*𝟓^2,kG*τ/(𝟐^7*𝟑^4*𝟓^3)
const mₑ,μ₀ = αinv^2*R∞*𝟐*𝘩/𝘤,𝟐*𝘩/𝘤*α/𝘦^2 # ≈ 4π*(1e-7+5.5e-17), exact charge
const ħ,μₚₑ,Rᵤ,αL,αG,Mᵤ = 𝘩/τ,μₚᵤ/μₑᵤ,NA*kB,centi/𝘤,(mₑ/mP)^2,NA*mₑ/μₑᵤ
const pc,G,DAY,nm = au*𝟐^7*𝟑^4*𝟓^3/τ,𝘤*ħ/mP^2,𝟐^7*𝟑^3*𝟓^2,sqrt(GME/g₀)*τ/𝟐^5/𝟑^3/𝟓^2
const GM☉ =au^3*k^2/DAY^2; const th = 𝟏𝟎^3*pc/H0; const ΛC = 𝟑*ΩΛ*(th*𝘤)^-2
const lc,mc,ρΛ = 𝟐*sqrt(τ/ΛC),𝘤^2/(𝟐*sqrt(τ*ΛC*G)),ΛC*𝘤^4/(𝟐^2*τ)/G
const lcq,mcq = sqrt.(sqrt.((𝘤*ħ/ρΛ,ρΛ*ħ^3/𝘤^5)))
const tcq,em,mi = lcq*sqrt(mcq/sqrt(sqrt(ρΛ*(𝘤*ħ)^3))),sqrt(GME/g₀)*τ/𝟐^9/𝟓^7,𝟐^5*𝟑*𝟓*𝟏𝟏

@pure sackurtetrode(U::UnitSystem,P=atm,T=𝟏,m=atomicmass(U)) = normal(log((Constant(exp(5/2))*boltzmann(U)*sqrt(boltzmann(U)/gravity(U)/turn(U)/planckreduced(U)^2)^3)*(T/P*sqrt(m*T)^3)))

const Universe = Coupling(αG,α,μₑᵤ,μₚᵤ,ΩΛ)

Constant(::UnitSystem{a,b,c,d,e,f,g}) where {a,b,c,d,e,f,g} = UnitSystem(Constant(a),Constant(b),Constant(c),Constant(d),Constant(e),Constant(f),Constant(g[1]),Constant(g[2]),Constant(g[3]),Constant(g[4]),Constant(g[5]),Universe,Constant(g[7]),𝟐,𝟑,𝟓,𝟕,𝟏𝟏,𝟏𝟗,𝟒𝟑)

export MetricSystem, ConventionalSystem, RankineSystem
export AstronomicalSystem, ElectricSystem, GaussSystem, EntropySystem

"""
    MetricSystem(Mu=Mᵤ,μ0=μ₀,Ru=Rᵤ,g0=𝟏,h=𝘩)

Constructs new `UnitSystem` from `molarmass` constant, `vacuumpermeability`, `molargas` constant, `gravity` force reference, and `planck` constant.

Examples include `SI2019`, `Metric`, `SI2019Engineering`, `MetricEngineering`, `SI1976`. In addition, the `ConventionalSystem` constructor further builds on `MetricSystem`, resulting in variations.
"""
MetricSystem(Mu=Mᵤ,μ0=μ₀,Ru=Rᵤ,g0=𝟏,h=𝘩,me=αinv^2*R∞*𝟐*h/𝘤) = UnitSystem(Ru*me/Mu/μₑᵤ/g0,h/τ/g0,𝘤,μ0,me,Mu,Kcd*(mₑ/me)^2*(h/𝘩)*g0,𝟏,𝟏,𝟏,g0,Universe,τ,𝟐,𝟑,𝟓,𝟕,𝟏𝟏,𝟏𝟗,𝟒𝟑)

"""
    ConventionalSystem(RK,KJ,Ru=Rᵤ,g0=𝟏) = MetricSystem(milli,𝟐*RK/𝘤*α,Ru,g0,𝟐^2/RK/KJ^2)

Constructs new `UnitSystem` from von `klitzing` constant and `josephson` constant, with an optional specification of `universal` gas constant and `gravity` reference constant.

Examples include `Conventional` (based on 1990) and `CODATA` (based on 2014).
"""
ConventionalSystem(klitz,joseph,Ru=Rᵤ,g0=𝟏) = MetricSystem(milli,𝟐*klitz/𝘤*α,Ru,g0,(𝟐*𝟐)/klitz/(joseph*joseph))

"""
    RankineSystem(U::UnitSystem,l,m,g0=𝟏) = EntropySystem(U,𝟏,l,m,°R,vacuumpermeability(U)/m/l/g0,kilo*molarmass(U),g0)

Constructs new `UnitSystem` from `U` rescaled along `length` and `mass` with optional `gravity` reference constant used to define technical and engineering units.

Examples: `FPS`, `British`, `IPS`, `English`, `Survey`.
"""
RankineSystem(u,l,m,g0=𝟏) = EntropySystem(u,𝟏,l,m,°R,UnitSystems.vacuumpermeability(u)/(m*l)/g0,UnitSystems.unit(kilo*UnitSystems.molarmass(u)),g0)

# historical units

const SI2019 = Quantity(MetricSystem())
const Metric = Quantity(MetricSystem(milli,τ/𝟐^6/𝟓^7))
const SI2019Engineering = Quantity(MetricSystem(Mᵤ,μ₀/g₀,Rᵤ,g₀))
const MetricEngineering = Quantity(MetricSystem(milli,τ/𝟐^6/𝟓^7/g₀,Rᵤ,g₀))
const SI1976 = Quantity(MetricSystem(milli,τ/𝟐^6/𝟓^7,Constant(8.31432)))
const CODATA = Quantity(ConventionalSystem(RK2014,KJ2014,Rᵤ2014))
const Conventional = Quantity(ConventionalSystem(RK1990,KJ1990))
const International = Quantity(ElectricSystem(Metric,Ωᵢₜ,Vᵢₜ))
const InternationalMean = Quantity(ElectricSystem(Metric,Constant(1.00049),Constant(1.00034)))

const EMU = Quantity(GaussSystem(Metric,𝟏,𝟐*τ))
const ESU = Quantity(GaussSystem(Metric,(hecto*𝘤)^-2,𝟐*τ))
const Gauss = Quantity(GaussSystem(Metric,𝟏,𝟐*τ,centi/𝘤))
const LorentzHeaviside = Quantity(GaussSystem(Metric,𝟏,𝟏,centi/𝘤))
#const Thomson = Quantity(GaussSystem(Metric,𝟏,𝟐*τ,inv(𝟐)))
const Kennelly = Quantity(GaussSystem(Metric,inv(𝟏𝟎^7),𝟐*τ,𝟏,𝟏,𝟏))

const British = Quantity(RankineSystem(Metric,ft,lb*g₀/ft))
#const British2019 = Quantity(RankineSystem(SI2019,ft,lb*g₀/ft))
const Survey = Quantity(RankineSystem(Metric,ftUS,lb,g₀/ftUS))
#const Survey2019 = Quantity(RankineSystem(SI2019,ftUS,lb,g₀/ftUS))
const English = Quantity(RankineSystem(Metric,ft,lb,g₀/ft))
#const English2019 = Quantity(RankineSystem(SI2019,ft,lb,g₀/ft))
const FPS = Quantity(RankineSystem(Metric,ft,lb))
#const FPS2019 = Quantity(RankineSystem(SI2019,ft,lb))
const IPS = Quantity(RankineSystem(Metric,ft/𝟐^2/𝟑,lb*g₀*𝟐^2*𝟑/ft))
#const IPS2019 = Quantity(RankineSystem(SI2019,ft/𝟐^2/𝟑,lb*g₀*𝟐^2*𝟑/ft))

#const Astronomical = Quantity(AstronomicalSystem(Metric))
const Hubble = Quantity(EntropySystem(Metric,th,𝘤*th,𝟏))
const Cosmological = Quantity(EntropySystem(Metric,lc/𝘤,lc,mc))
const CosmologicalQuantum = Quantity(EntropySystem(Metric,tcq,lcq,mcq))
#const EMU2019 = Quantity(EntropySystem(SI2019,𝟏,centi,milli))
#const ESU2019 = Quantity(EntropySystem(SI2019,𝟏,centi,milli,𝟏,kilo*μ₀/𝘤^2))
#const Mixed = Quantity(EntropySystem(Metric,𝟏,𝟏,𝟏,𝟏,μ₀))
const Nautical = Quantity(EntropySystem(Metric,HOUR,nm,em^3,𝟏,τ*𝟑^3/𝟐^10/𝟓^12,milli))
const Meridian = Quantity(EntropySystem(Metric,𝟏,em,em^3,𝟏,τ/𝟐^6/𝟓^7,milli))
const MeridianEngineering = Quantity(EntropySystem(Metric,𝟏,em,em^3,𝟏,τ/𝟐^6/𝟓^7*em/g₀,milli,g₀/em))
const GravitationalSI2019 = Quantity(EntropySystem(SI2019,𝟏,𝟏,g₀))
const GravitationalMetric = Quantity(EntropySystem(Metric,𝟏,𝟏,g₀))
const GravitationalMeridian = Quantity(EntropySystem(Metric,𝟏,em,g₀*em^2,𝟏,τ/𝟐^6/𝟓^7*em/g₀,milli))
const IAU☉ = Quantity(EntropySystem(Metric,DAY,au,GM☉/G))
const IAUE = Quantity(EntropySystem(Metric,DAY,LD,GME/G))
const IAUJ = Quantity(EntropySystem(Metric,DAY,JD,GMJ/G))
const MTS = Quantity(EntropySystem(Metric,𝟏,𝟏,kilo))
const KKH = Quantity(EntropySystem(Metric,HOUR,kilo,𝟏))
const MPH = Quantity(EntropySystem(FPS,HOUR,mi,𝟏))
const FFF = Quantity(EntropySystem(Metric,𝟕*𝟐*DAY,fur,(𝟐*𝟑^2*𝟓)*lb,°R,Constant(0.),𝟏))

# natural units

const Planck = Quantity(UnitSystem(𝟏,𝟏,𝟏,𝟏,√(𝟐*τ*αG)))
const PlanckGauss = Quantity(UnitSystem(𝟏,𝟏,𝟏,𝟐*τ,√αG))
const Stoney = Quantity(UnitSystem(𝟏,αinv,𝟏,𝟐*τ,√(αG*αinv)))
const Hartree = Quantity(UnitSystem(𝟏,𝟏,αinv,𝟐*τ*α^2,𝟏))
const Rydberg = Quantity(UnitSystem(𝟏,𝟏,𝟐*αinv,τ/𝟐*α^2,inv(𝟐)))
const Schrodinger = Quantity(UnitSystem(𝟏,𝟏,αinv,𝟐*τ*α^2,√(αG*αinv)))
const Electronic = Quantity(UnitSystem(𝟏,αinv,𝟏,𝟐*τ,𝟏))
const Natural = Quantity(UnitSystem(𝟏,𝟏,𝟏,𝟏,𝟏))
const NaturalGauss = Quantity(UnitSystem(𝟏,𝟏,𝟏,𝟐*τ,𝟏))
const QCD = Quantity(UnitSystem(𝟏,𝟏,𝟏,𝟏,inv(μₚₑ)))
const QCDGauss = Quantity(UnitSystem(𝟏,𝟏,𝟏,𝟐*τ,inv(μₚₑ)))
const QCDoriginal = Quantity(UnitSystem(𝟏,𝟏,𝟏,𝟐*τ*α,inv(μₚₑ)))

export SI, MKS, SIE, ME, GSI2019, GSI, GM, CGS, CGS2019, CGSm, CGSe, HLU, FFF, AE, EE, BG
export EnglishEngineering, BritishGravitational, AbsoluteEnglish, EnglishUS, EE2019, IAU
const SI, MKS, SIE, ME, IAU = SI2019, Metric, SI2019Engineering, MetricEngineering, IAU☉
const GSI2019, GSI, GM = GravitationalSI2019, GravitationalSI2019, GravitationalMetric
const CGS, CGSm, CGSe, HLU = Gauss, EMU, ESU, LorentzHeaviside
const EnglishEngineering, BritishGravitational, BG = English, British, British
const EnglishUS, AbsoluteEnglish, AE, EE = Survey, FPS, FPS, English
#const EnglishEngineering2019,BritishGravitational2019,BG2019 = English2019,British2019,British2019
#const EE2019,AbsoluteEnglish2019,AE2019 = English2019,FPS2019,FPS2019

for u ∈ Systems
    @eval unitname(::typeof(normal($u))) = $(string(u))
end

ext: jl


text:
## FFC for selected failure scenarios
function FFC(GRB_ENV, edges, capacity, flows, demand, scenarios, T, Tf)
    printstyled("\n** solving FFC LP...\n", color=:yellow)

    nedges = size(edges,1)
    nflows = size(flows,1)
    ntunnels = size(T,1)

    tunnel_num = 0
    for x in 1:length(Tf)
        if size(Tf[x],1) > tunnel_num
            tunnel_num = size(Tf[x],1)
        end
    end

    # CREATE RESIDUAL TUNNELS BY SCENARIO BY FLOW (References Tf): fiber cut scenarios
    Tsf = []
    for s in 1:size(scenarios,1)
        sft = []
        for f in 1:size(Tf,1)
            ft = []
            for t in 1:size(Tf[f],1)
                up = true
                if (length(T[Tf[f][t]]) == 0)
                    up = false
                end
                for e in T[Tf[f][t]]
                    if scenarios[s][e] == 0
                        up = false
                    end
                end
                if up
                    push!(ft, t)
                end
            end
            push!(sft,ft)
        end
        push!(Tsf,sft)
    end

    #CREATE TUNNEL EDGE MATRIX
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    model = Model(() -> Gurobi.Optimizer(GRB_ENV))
    set_optimizer_attribute(model, "OutputFlag", 0)
    set_optimizer_attribute(model, "Threads", 32)

    model_def_time = @elapsed begin

    @variable(model, b[1:nflows] >= 0)
    @variable(model, a[1:nflows,1:tunnel_num] >= 0)

    for f in 1:nflows
        @constraint(model, sum(a[f,t] for t in 1:size(Tf[f],1)) >= b[f])   #the sum of all allocated bandwidths on every flow must be >= the total bandwidth for that flow
    end

    for e in 1:nedges
        @constraint(model, sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows for t in 1:size(Tf[f],1)) <= capacity[e])   #overlapping flows cannot add up to the capacity of that link
    end

    for f in 1:nflows
        for s in 1:size(scenarios,1)
            @constraint(model, sum(a[f,t] for t in Tsf[s][f]) >= b[f])   #residual tunnels must be able to carry bandwidth
        end
    end

    for f in 1:nflows
        @constraint(model, b[f] >= 0)
    end

    for f in 1:nflows
        @constraint(model, b[f] <= demand[f])   #all allocated bandwidths must be less than the demand for that flow
        for t in 1:size(Tf[f],1)
            @constraint(model, a[f,t] >= 0)     #each allocated bandwidth on for flow f on tunnel t >= 0
        end
    end

    @objective(model, Max, sum((b[i] for i in 1:size(b,1))))
    end  # model_def_time end

    optimize!(model)
    solve_runtime = solve_time(model)
    opt_runtime = solve_runtime + model_def_time

    return value.(a), value.(b), objective_value(model), solve_runtime, opt_runtime
end

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow


text:
# Run Experiments, display results (and possibly save data) on a Bandit Problem to be specified

using HDF5

# DO YOU WANT TO SAVE RESULTS?
typeExp = "Save"
# typeExp = "NoSave"

# TYPE OF DISTRIBUTION
@everywhere typeDistribution="Bernoulli"
@everywhere include("BAIalgos.jl")

# CHANGE NAME (save mode)
fname="results/Experiment4arms"

# BANDIT PROBLEM
@everywhere mu=vec([0.3 0.25 0.2 0.1])
@everywhere best=find(mu.==maximum(mu))[1]
K=length(mu)

# RISK LEVEL
delta=0.1

# NUMBER OF SIMULATIONS
N=100

# OPTIMAL SOLUTION
@everywhere v,optWeights=OptimalWeights(mu)
@everywhere gammaOpt=optWeights[best]
print("mu=$(mu)\n")
print("Theoretical number of samples: $((1/v)*log(1/delta))\n")
print("Optimal weights: $(optWeights)\n\n")

# POLICIES

@everywhere ChernoffPTSHalf(x,y,z)=ChernoffPTS(x,y,z,0.5)
@everywhere ChernoffPTSOpt(x,y,z)=ChernoffPTS(x,y,z,gammaOpt)

policies=[TrackAndStop,ChernoffBC2,ChernoffPTSHalf,ChernoffPTSOpt,KLLUCB,UGapE]
names=["TrackAndStop","ChernoffBC","ChernoffPTS","ChernoffPTSOpt","KLLUCB","UGapE"]

# EXPLORATION RATES
@everywhere explo(t,n,delta)=log((log(t)+1)/delta)

lP=length(policies)
rates=[explo for i in 1:lP]

# RUN EXPERIMENTS

function MCexp(mu,delta,N)
	for imeth=1:lP
            policy=policies[imeth]
            beta=rates[imeth]
	    startTime=time()
	    a = Array(RemoteRef, N)
	    for j in 1:N
		a[j] =  @spawn policy(mu,delta,beta)
	    end
	    res = collect([fetch(a[j]) for j in 1:N])
            proportion = zeros(N,K)
            for j in 1:N
                n=res[j][2]
                proportion[j,:]=n/sum(n)
            end
            NonTerminated = sum(collect([res[j][1]==0 for j in 1:N]))
            FracNT=(NonTerminated/N)
            FracReco=zeros(K)
            for k in 1:K
                FracReco[k]=sum([(res[j][1]==k)?1:0 for j in 1:N])/(N*(1-FracNT))
            end
            print("Results for $(policy), average on $(N) runs\n")
	    print("proportion of runs that did not terminate: $(FracNT)\n")
	    print("average number of draws: $(sum([sum(res[j][2]) for j in 1:N])/(N-NonTerminated))\n")
            print("average proportion of draws: \n $(mean(proportion,1))\n")
	    print("proportion of errors: $(1-sum([res[j][1]==best for j in 1:N])/(N-NonTerminated))\n")
            print("proportion of recommendation made when termination: $(FracReco)\n")
            print("elapsed time: $(time()-startTime)\n\n")
	    print("")
	end
end

function SaveData(mu,delta,N)
        K=length(mu)
        for imeth=1:lP
            Draws=zeros(N,K)
            policy=policies[imeth]
            namePol=names[imeth]
            startTime=time()
	    a = Array(RemoteRef, N)
            rate=rates[imeth]
	    for j in 1:N
		a[j] =  @spawn policy(mu,delta,rate)
	    end
	    res = [fetch(a[j]) for j in 1:N]
            proportion=zeros(N,K)
            for k in 1:N
               r=res[k][2]
               Draws[k,:]=res[k][2]
               proportion=r/sum(r)
            end
            Reco=[res[j][1] for j in 1:N]
            Error=collect([(r==best)?0:1 for r in Reco])
            FracNT=sum([r==0 for r in Reco])/N
            FracReco=zeros(K)
            for k in 1:K
                FracReco[k]=sum([(r==k)?1:0 for r in Reco])/(N*(1-FracNT))
            end
            print("Results for $(policy), average on $(N) runs\n")
	    print("proportion of runs that did not terminate: $(FracNT)\n")
	    print("average number of draws: $(sum(Draws)/(N*(1-FracNT)))\n")
            print("average proportions of draws: $(mean(proportion,1))\n")
	    print("proportion of errors: $(sum(Error)/(float(N*(1-FracNT))))\n")
            print("proportion of recommendation made when termination: $(FracReco)\n")
            print("elapsed time: $(time()-startTime)\n\n")
            name="$(fname)_$(namePol)_delta_$(delta)_N_$(N).h5"
            h5write(name,"mu",mu)
            h5write(name,"delta",delta)
            h5write(name,"FracNT",collect(FracNT))
            h5write(name,"FracReco",FracReco)
            h5write(name,"Draws",Draws)
            h5write(name,"Error",Error)
	end
end

if (typeExp=="Save")
   SaveData(mu,delta,N)
else
   MCexp(mu,delta,N)
end

max_stars_repo_name: jsfunc/best-arm-identification

max_issues_repo_name: ClaireVern/best-arm-identification

max_forks_repo_name: ClaireVern/best-arm-identification


text:
abstract type AbstractParameterSpace end

mutable struct SetupBatchC{T<:AbstractGenomicData}
    Isolate::T
    n_polymers::Int64
    n_monomers::Int64
    n_microbes::Int64
    n_enzymes::Int64
    n_minerals::Int64
    u0_monomers::Array{Float64,1}
    u0_reserves::Array{Float64,1}
    u0_structures::Array{Float64,1}
    function SetupBatchC{T}(Isolate::T, n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals) where {T<:AbstractGenomicData}
        u0_monomers = zeros(n_monomers)
        u0_reserves = 0.5*1e11.*0.47*cell_volume_to_dry_mass(genome_size_to_cell_volume(Isolate.L_DNA))./12.011
        u0_structures = 0.5*1e11.*0.47*cell_volume_to_dry_mass(genome_size_to_cell_volume(Isolate.L_DNA))./12.011
        new(Isolate, n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals, u0_monomers, u0_reserves, u0_structures)
    end
end

SetupBatchC(Isolate::T, n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals) where {T<:AbstractGenomicData} = SetupBatchC{T}(Isolate::T, n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals)

@units @description mutable struct ParamsBatchC{B<:SetupBatchC, D<:AbstractMonomer, H<:AbstractMetabolism, A<:AbstractAssimilation} <:AbstractParameterSpace
    Setup::B                             | "System dimensions and initial conditions"
    Monomer::D                           | "Monomer defined by BaseMonomer"
    PMetabolismC::H                      | "C metabolism parameters"
    PAssimilation::A                     | "Assimilation parameters"
    #
    n_polymers::Int64
    n_monomers::Int64
    n_microbes::Int64
    n_enzymes::Int64
    n_minerals::Int64
    dim::Int64
    #
    y_DE::Array{Float64,1}   | mol/mol   | "Yield on assimilation"
    N_C::Array{Float64,1}                | "Number of C atoms"
    V_S::Array{Float64,2}    |mol/(mol*d)| "Transporter specific rate"
    N_SB::Array{Float64,2}   |mol/mol    | "Number of uptake sites"
    K_D:: Array{Float64,2}   |mol/m^3    | "Affinity constant"
    #
    γ_X::Array{Float64,1}    | d^-1      | "enzyme turnover rate"
    γ_V::Array{Float64,1}    | d^-1      | "microbial turnover rate"
    β::Array{Float64,1}                  | "density dependent turnover"
    f_E::Array{Float64,1}                | "fractionation of dead reserve into monomers"
    #
    I_D::Array{Float64,1}    |mol/(m^3*d)| "Monomer input"
    #
    function ParamsBatchC{B, D, H, A}(Setup::B, Monomer::D, PMetabolismC::H, PAssimilation::A) where {B<:SetupBatchC, D<:AbstractMonomer, H<:AbstractMetabolism, A<:AbstractAssimilation}
        n_polymers = Setup.n_polymers
        n_monomers = Setup.n_monomers
        n_microbes = Setup.n_microbes
        n_enzymes = Setup.n_enzymes
        n_minerals = Setup.n_minerals
        dim = n_polymers + n_monomers + 3*n_microbes + n_enzymes
        y_DE = Monomer.y_DE
        N_C  = Monomer.N_C
        V_S  = 100.0*24*60*60*ones(n_monomers,n_microbes)
        N_SB = PAssimilation.N_SB
        K_D = PAssimilation.K_D
        γ_X = 1e-4*ones(n_enzymes)
        γ_V = 1e-3*ones(n_microbes)
        β   = ones(n_microbes)
        f_E = ones(n_monomers).*1/n_monomers
        I_D = zeros(n_monomers)
        new(Setup, Monomer, PMetabolismC, PAssimilation, n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals, dim, y_DE, N_C, V_S, N_SB, K_D, γ_X, γ_V, β, f_E, I_D)
    end
end

ParamsBatchC(Setup::B, Monomer::D, PMetabolismC::H, PAssimilation::A) where {B<:SetupBatchC, D<:AbstractMonomer, H<:AbstractMetabolism, A<:AbstractAssimilation} = ParamsBatchC{B, D, H, A}(Setup::B, Monomer::D, PMetabolismC::H, PAssimilation::A)

function rhs_BatchC!(du, u, p::AbstractParameterSpace, t)
    D                            = u[1+p.n_polymers:p.n_polymers+p.n_monomers]
    E                            = u[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes]
    V                            = u[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes]
    X                            = u[1+p.n_polymers+p.n_monomers+2*p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes+p.n_enzymes]
    CO2                          = u[1+p.n_polymers+p.n_monomers+2*p.n_microbes+p.n_enzymes:p.n_polymers+p.n_monomers+2*p.n_microbes+p.n_enzymes+p.n_microbes]
    #
    r                            = growth!(zeros(p.n_microbes), E, V, p.PMetabolismC)
    x, rG_CO2, rM_CO2, rX_CO2    = growth_production(r, E, V, p.PMetabolismC)
    ECA_D                        = DEBTrait.ECA_kinetics!(zeros(p.n_monomers,p.n_microbes), D, V, p.K_D, p.V_S, p.N_SB)
    #
    F_D                          = vcat(sum(ECA_D[:, 1:p.n_microbes], dims=2)...).*p.N_C
    F_DE                         = vcat(sum((1.0./p.y_DE .- 1.0).*p.N_C.*ECA_D[:, 1:p.n_microbes], dims=1)...)
    F_DE_CO2                     = vcat(sum((1.0./p.y_DE).*p.N_C.*ECA_D[:, 1:p.n_microbes], dims=1)...)
    #
    F_X                          = sum(x.*V, dims=1)
    #
    F_ED                         = p.f_E.*sum(p.γ_V.*E.^p.β, dims=1)
    #
    @. du[1+p.n_polymers:p.n_polymers+p.n_monomers] = p.I_D - F_D + F_ED +  p.γ_X*X
    @. du[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes] = F_DE - (p.PMetabolismC.k_E - r)*E - p.γ_V*E^p.β
    @. du[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes] = r*V - p.γ_V*V^p.β
    @. du[1+p.n_polymers+p.n_monomers+2*p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes+p.n_enzymes] = F_X - p.γ_X*X
    @. du[1+p.n_polymers+p.n_monomers+2*p.n_microbes+p.n_enzymes:p.n_polymers+p.n_monomers+2*p.n_microbes+p.n_enzymes+p.n_microbes]  = rG_CO2 + rM_CO2 + rX_CO2 + F_DE_CO2
    return du
end

function calc_BR(du, u, p, t)
    out = rhs_BatchC!(du, u, p, t)
    BR = out[1+p.n_polymers+p.n_monomers+2*p.n_microbes+p.n_enzymes:p.n_polymers+p.n_monomers+2*p.n_microbes+p.n_enzymes+p.n_microbes]
    E = u[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes]
    V = u[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes]
    V_cell = DEBTrait.genome_size_to_cell_volume(p.Setup.Isolate.L_DNA)
    Biomass = @. E + V
    N_cells = Biomass.*12.011*0.47./DEBTrait.cell_volume_to_dry_mass(V_cell)
    BR_norm = BR./N_cells
end

function calc_BP(du, u, p, t)
    out = rhs_BatchC!(du, u, p, t)
    BP_E = out[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes]
    BP_V = out[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes]
    BP = @. BP_E + BP_V
    V_cell = DEBTrait.genome_size_to_cell_volume(p.Setup.Isolate.L_DNA)
    E = u[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes]
    V = u[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes]
    Biomass = @. E + V
    N_cells = Biomass.*12.011*0.47./DEBTrait.cell_volume_to_dry_mass(V_cell)
    BP_norm = BP./N_cells
end

function calc_BGE(du,u,p,t)
    BP = calc_BP(du,u,p,t)
    BR = calc_BR(du,u,p,t)
    denom = max.(1e-15, BP.+BR)
    BGE = @. max.(1e-15, BP/denom)
end

function calc_growth_rate(u, p)
    E = u[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes]
    V = u[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes]
    r = growth!(zeros(p.n_microbes), E, V, p.PMetabolismC)
    return r
end

function calc_growth_production(u, p)
    E = u[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes]
    V = u[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes]
    r = growth!(zeros(p.n_microbes), E, V, p.PMetabolismC)
    x, rG_CO2, rM_CO2, rX_CO2 = growth_production(r, E, V, p.PMetabolismC)
    V_cell = DEBTrait.genome_size_to_cell_volume(p.Setup.Isolate.L_DNA)
    Biomass = @. E + V
    N_cells = Biomass.*12.011*0.47./DEBTrait.cell_volume_to_dry_mass(V_cell)
    return x./N_cells, rG_CO2./N_cells, rM_CO2./N_cells, rX_CO2./N_cells
end

function calc_specific_assimilation_rate(u, p)
    D = u[1+p.n_polymers:p.n_polymers+p.n_monomers]
    V = u[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes]
    ECA_D = DEBTrait.ECA_kinetics!(zeros(p.n_monomers,p.n_microbes), D, V, p.K_D, p.V_S, p.N_SB)
    F_D   = vcat(sum(ECA_D[:, 1:p.n_microbes], dims=2)...).*p.N_C
    F_DE  = vcat(sum((1.0./p.y_DE .- 1.0).*p.N_C.*ECA_D[:, 1:p.n_microbes], dims=1)...)
    E = u[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes]
    V_cell = DEBTrait.genome_size_to_cell_volume(p.Setup.Isolate.L_DNA)
    Biomass = @. E + V
    N_cells = Biomass.*12.011*0.47./DEBTrait.cell_volume_to_dry_mass(V_cell)
    return F_DE./N_cells
end

function calc_assimilation_rate(u, p)
    D = u[1+p.n_polymers:p.n_polymers+p.n_monomers]
    V = u[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes]
    ECA_D = DEBTrait.ECA_kinetics!(zeros(p.n_monomers,p.n_microbes), D, V, p.K_D, p.V_S, p.N_SB)
    F_D   = vcat(sum(ECA_D[:, 1:p.n_microbes], dims=2)...).*p.N_C
    F_DE  = vcat(sum((1.0./p.y_DE .- 1.0).*p.N_C.*ECA_D[:, 1:p.n_microbes], dims=1)...)
end

function calc_specific_assimilation_respiration(u, p)
    D = u[1+p.n_polymers:p.n_polymers+p.n_monomers]
    V = u[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes]
    ECA_D = DEBTrait.ECA_kinetics!(zeros(p.n_monomers,p.n_microbes), D, V, p.K_D, p.V_S, p.N_SB)
    E = u[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes]
    V_cell = DEBTrait.genome_size_to_cell_volume(p.Setup.Isolate.L_DNA)
    Biomass = @. E + V
    N_cells = Biomass.*12.011*0.47./DEBTrait.cell_volume_to_dry_mass(V_cell)
    F_DE_CO2  = vcat(sum((1.0./p.y_DE).*p.N_C.*ECA_D[:, 1:p.n_microbes], dims=1)...)
    return F_DE_CO2./N_cells
end

function calc_assimilation_growth_ratio(u,p)
    F_DE_CO2 = calc_assimilation_rate(u,p)
    r        = calc_growth_rate(u,p)
    return F_DE_CO2./r
end

function calc_reserve_density(u, p)
    E = u[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes]
    V = u[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes]
    m_E = @. E/V
end

function calc_specific_biomass(u, p)
    E = u[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes]
    V = u[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes]
    Bio = @. E+V
    V_cell = DEBTrait.genome_size_to_cell_volume(p.Setup.Isolate.L_DNA)
    n_0 = DEBTrait.cell_volume_to_cell_number_density(V_cell)
    return Bio.*n_0
end


text:
function lintlintpragma(ex::Expr, ctx::LintContext)
    if length(ex.args) >= 2 && isa(ex.args[2], AbstractString)
        m = match(r"^((Print)|(Info)|(Warn)|(Error)) ((type)|(me)|(version)) +(.+)"s, ex.args[2])
        if m != nothing
            action = m.captures[1]
            infotype = m.captures[6]
            rest_str = m.captures[10]
            if infotype == "type"
                v = parse(rest_str)
                if isexpr(v, :incomplete)
                    msg(ctx, :E138, rest_str, "incomplete pragma expression")
                    str = ""
                else
                    str = "typeof($rest_str) == $(guesstype(v, ctx))"
                end
            elseif infotype == "me"
                str = rest_str
            elseif infotype == "version"
                v = convert(VersionNumber, rest_str)
                reachable = ctx.versionreachable(v)
                if reachable
                    str = "Reachable by " * string(v)
                else
                    str = "Unreachable by " * string(v)
                end
            end

            if action == "Print"
                println(str)
            elseif action == "Info"
                msg(ctx, :I271, str)
            elseif action == "Warn"
                msg(ctx, :W241, str)
            else
                msg(ctx, :E221, str)
            end
        else
            if !ctx.versionreachable(VERSION)
                return
            end
            pragma!(ctx.current, ex.args[2], location(ctx))
        end
    else
        msg(ctx, :E137, "lintpragma must be called using only string literals")
    end
end

function pragmaexists(s::String, ctx::_LintContext; deep=true)
    if s in keys(pragmas(ctx))
        true
    elseif istoplevel(ctx)
        false
    elseif deep
        pragmaexists(s, parent(ctx); deep=true)
    else
        false
    end
end

pragmaexists(s::AbstractString, ctx::_LintContext; deep=true) =
    pragmaexists(String(s), ctx; deep = deep)


text:
using BakSneppenSimulation
using BakSneppenSimulation: n, diameter
using Test
using Random

@testset "BakSneppenSimulation.jl" begin

    @testset "3 mutations" begin
        Random.seed!(1234)

        config = SimConfig(10, 0.6, 3)
        state = SimState(rand(0:1, 10))
        sim = SimCircular(config, state)

        @test state.xlist == [0, 1, 1, 1, 0, 0, 0, 1, 0, 1]
        @test state.x0list == [1, 5, 6, 7, 9]

        mutate(sim)

        @test state.xlist == [0, 1, 1, 1, 1, 0, 0, 1, 0, 1]
        @test state.x0list == findall(x->x==0, state.xlist)
        @test diameter(sim) == 6

        mutate(sim)

        @test state.xlist == [0, 1, 1, 1, 1, 1, 0, 1, 0, 1]
        @test state.x0list == findall(x->x==0, state.xlist)
        @test diameter(sim) == 5

        mutate(sim)

        @test state.xlist == [0, 1, 1, 1, 1, 1, 0, 0, 0, 0]
        @test sort!(state.x0list) == findall(x->x==0, state.xlist)
        @test diameter(sim) == 5
    end

    @testset "2 mutations" begin
        Random.seed!(1234)

        config = SimConfig(10, 0.6, 2)
        state = SimState(rand(0:1, 10))
        sim = SimLinear(config, state, true)

        @test state.xlist == [0, 1, 1, 1, 0, 0, 0, 1, 0, 1]
        @test state.x0list == [1, 5, 6, 7, 9]

        mutate(sim)

        @test state.xlist == [0, 1, 1, 1, 1, 1, 0, 1, 0, 1]
        @test state.x0list == findall(x->x==0, state.xlist)
        @test diameter(sim) == 9

        mutate(sim)

        @test state.xlist == [0, 1, 1, 1, 1, 1, 0, 1, 0, 1]
        @test state.x0list == findall(x->x==0, state.xlist)
        @test diameter(sim) == 9

        mutate(sim)

        @test state.xlist == [1, 0, 1, 1, 1, 1, 0, 1, 0, 1]
        @test sort!(state.x0list) == findall(x->x==0, state.xlist)
        @test diameter(sim) == 8
    end
end


text:
using EffectiveWaves, Test

@testset "match purely numerical solution" begin

## high attenuating material
    medium = Acoustic(2; ρ=1.0, c=1.0)
    # cs = [0.1,0.5]
    cs = [0.1,0.5]
    ms = MultipleScattering

    species = [
        Specie(Particle(
            Acoustic(2; ρ=c, c=0.6-c), ms.Circle(2.0*c));
            volume_fraction=c
        )
    for c in cs]

    ω = 1.1
    k = ω/medium.c
    θin = 0.3
    tol = 1e-8
    basis_order = 2

    normal = [-1.0,0.0] # an outward normal to the surface
    materials = [Material(Halfspace(normal),s) for s in species]
    source = PlaneSource(medium, [cos(θin),sin(θin)])

    material = materials[1]

    # import StaticArrays: SVector
    #
    # function unreachable_test(ω::T, source::AbstractSource, material::Material{Dim,S,Sps}; kws...) where {T,Dim,S<:Shape{Dim},Sps<:Species{T,Dim}}
    #
    #     k_effs = rand(2) + rand(2) .* im
    #     wave_effs = map(k_effs) do k_eff
    #         # amp = Complex{Float64}[0.000326919-0.000330135im; -0.000605635-0.0396968im; -0.880755+0.456138im; -0.000605635-0.0396968im; 0.000326919-0.000330135im];
    #         amp = rand(5) + rand(5) .* im;
    #         # wavevector = Complex{Float64}[1.62185+1.33695im, 0.325072+0.0im]
    #         wavevector = rand(2) + rand(2) .* im
    #         basis_order = 2
    #         EffectivePlaneWaveMode(ω, basis_order, SVector(wavevector...), amp)
    #         # WaveMode(ω, k_eff, source, material; kws...)
    #     end
    #
    #     return wave_effs
    # end
    #
    # data = [
    #     unreachable_test(ω, source, materials[i])
    # for i in eachindex(species)]

    wave_effs_arr = Vector{Vector{EffectivePlaneWaveMode{Float64,2}}}(undef,length(species))

    i = 1
    wave_effs_arr[i] = WaveModes(ω, source, materials[i];
        basis_order=basis_order,
        box_k = [[-18.0,18.0],[0.0,20.0]],
        num_wavenumbers = 6,
        tol = tol)

    # k_effs = wavenumbers_bisection(ω, medium, [species[1]];
    #     basis_order = basis_order,
    #     box_k = [[-18.0,18.0],[0.0,20.0]],
    #     num_wavenumbers = 8,
    #     tol = tol)

    i = 2
    wave_effs_arr[i] = WaveModes(ω, source, materials[i];
        basis_order=basis_order,
        num_wavenumbers=6,
        tol = tol);
        # extinction_rescale = false)

    # Causes unreachable error..
    # for i in eachindex(species)
    #     wave_effs_arr[i] = WaveModes(ω, source, materials[i];
    #         basis_order=basis_order,
    #         mesh_points=5,
    #         num_wavenumbers=5,
    #         tol = tol,
    #         extinction_rescale = false)
    # end

   # use only 6 least attenuating
   wave_effs_arr2 = [w[1:6] for w in wave_effs_arr];

    match_ws = [
        MatchPlaneWaveMode(ω, source, materials[i];
            basis_order=basis_order,
            tol = tol,
            wave_effs = wave_effs_arr2[i],
            max_size=300)
            # max_size=480)
    for i in eachindex(species)];

    @test maximum(match_error(match_ws[i],materials[i].shape) for i in eachindex(species)) < 200*tol

    avgs = [
        DiscretePlaneWaveMode(ω, source, materials[i];
                basis_order=basis_order,
                tol = tol, max_size=600,
                wave_effs = wave_effs_arr2[i])
    for i in eachindex(species)]

    R_ms = [reflection_coefficient(ω, match_ws[i], source, materials[i]) for i in eachindex(species)]
    R_ds = [reflection_coefficient(ω, avgs[i], source, materials[i]) for i in eachindex(species)]
    @test maximum(abs.(R_ms - R_ds)) < 5e-4

    avg_eff = DiscretePlaneWaveMode(match_ws[2].x_match[end]:0.002:40, match_ws[2].PlaneWaveModes, material.shape);
    R1 = reflection_coefficient(ω, avg_eff, source, materials[2])
    R2 = reflection_coefficient(ω, match_ws[2].PlaneWaveModes, source, materials[2]; x=avg_eff.x[1])
    @test norm(R1 - R2) < 1e-7

    map(eachindex(species)) do i
        j0 = findmin(abs.(avgs[i].x .- match_ws[i].x_match[1]))[2]
        x0 = avgs[i].x[j0+1:end]
        avg_m = DiscretePlaneWaveMode(x0, match_ws[i].PlaneWaveModes, materials[i].shape)
        maximum(abs.(avgs[i].amplitudes[j0+1:end,:,:][:] - avg_m.amplitudes[:]))
        @test norm(avgs[i].amplitudes[j0+1:end,:,:][:] - avg_m.amplitudes[:])/norm(avg_m.amplitudes[:]) < 1e-2
        @test maximum(abs.(avgs[i].amplitudes[j0+1:end,:,:][:] - avg_m.amplitudes[:])) < 3e-4
        maximum(abs.(avgs[i].amplitudes[j0+1:end,:,:][:] - avg_m.amplitudes[:]))
    end
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
module ArXivTools

include("download.jl")

end


text:
mutable struct PrioritizedExperienceReplay{CB<:CircularBuffer} <: AbstractWeightedReplay
    buffer::CB
    sumtree::SumTree
end

PrioritizedExperienceReplay(size, types, column_names) =
    new(CircularBuffer(size, types, column_names),
        SumTree{Int64}(size))

# size(er::PrioritizedExperienceReplay) = size(er.buffer)
@forward PrioritizedExperienceReplay.buffer getindex

function add!(er::PrioritizedExperienceReplay, experience, weight)
    idx = add!(er.buffer, experience)
    add!(er.sumtree, weight, idx)
    return
end

function update!(er::PrioritizedExperienceReplay, idx, weight)
end

function sample(er::PrioritizedExperienceReplay, batch_size; rng=Random.GLOBAL_RNG)
    batch_idx, batch_priorities, idx = sample(er.sumtree, batch_size; rng=rng)
    return getrow(er.buffer, idx)
end


text:
const PWGRID_VERSION = 3

mutable struct GVectorsW
    Ngwx::Int64
    idx_gw2r::Array{Int}
end

mutable struct GVectors
    Ng::Int64
    G::Array{Float64,2}
    G2::Array{Float64}
end

mutable struct PWGrid
    ecutwfc::Float64
    ecutrho::Float64
    Ns::Array{Int64}
    LatVecs::Array{Float64,2}
    RecVecs::Array{Float64,2}
    Ω::Float64
    r::Array{Float64,2}
    gvec::GVectors
    gvecw::GVectorsW
end

function PWGrid( ecutwfc::Float64, LatVecs::Array{Float64,2} )
    ecutrho = 4.0*ecutwfc
    #
    RecVecs = 2*pi*inv(LatVecs')
    Ω = det(LatVecs)
    #
    LatVecsLen = Array{Float64}(undef,3)
    LatVecsLen[1] = norm(LatVecs[:,1])
    LatVecsLen[2] = norm(LatVecs[:,2])
    LatVecsLen[3] = norm(LatVecs[:,3])

    Ns = Array{Int64}(undef,3)
    Ns[1] = 2*round( Int, sqrt(ecutrho/2)*LatVecsLen[1]/pi ) + 1
    Ns[2] = 2*round( Int, sqrt(ecutrho/2)*LatVecsLen[2]/pi ) + 1
    Ns[3] = 2*round( Int, sqrt(ecutrho/2)*LatVecsLen[3]/pi ) + 1

    # Use even sampling numbers
    Ns[1] = Ns[1] % 2 == 1 ? Ns[1] + 1 : Ns[1]
    Ns[2] = Ns[2] % 2 == 1 ? Ns[2] + 1 : Ns[2]
    Ns[3] = Ns[3] % 2 == 1 ? Ns[3] + 1 : Ns[3]

    Npoints = prod(Ns)
    r = init_grid_R( Ns, LatVecs )

    gvec = init_grid_G( Ns, RecVecs )
    gvecw = init_gvecw( ecutwfc, gvec.G2 )

    return PWGrid( ecutwfc, ecutrho, Ns, LatVecs, RecVecs, Ω, r, gvec, gvecw )
end

function mm_to_nn(mm::Int,S::Int)
    if mm > S/2
        return mm - S
    else
        return mm
    end
end

function init_grid_G( Ns, RecVecs )

    Ng = prod(Ns)

    G  = Array{Float64}(undef,3,Ng)
    G2 = Array{Float64}(undef,Ng)

    ig = 0
    for k in 0:Ns[3]-1
    for j in 0:Ns[2]-1
    for i in 0:Ns[1]-1
        ig = ig + 1
        gi = mm_to_nn( i, Ns[1] )
        gj = mm_to_nn( j, Ns[2] )
        gk = mm_to_nn( k, Ns[3] )
        G[1,ig] = RecVecs[1,1]*gi + RecVecs[1,2]*gj + RecVecs[1,3]*gk
        G[2,ig] = RecVecs[2,1]*gi + RecVecs[2,2]*gj + RecVecs[2,3]*gk
        G[3,ig] = RecVecs[3,1]*gi + RecVecs[3,2]*gj + RecVecs[3,3]*gk
        G2[ig] = G[1,ig]^2 + G[2,ig]^2 + G[3,ig]^2
    end
    end
    end

    return GVectors( Ng, G, G2 )
end

function init_gvecw( ecutwfc, G2 )
    idx_gw2r = findall( 0.5*G2 .< ecutwfc )
    Ngwx = length(idx_gw2r)
    return GVectorsW( Ngwx, idx_gw2r )
end

function find_edges( Ns )
    eS = Ns/2 .+ 0.5
    edges = []
    ip = 1
    for k in 0:Ns[3]-1
    for j in 0:Ns[2]-1
    for i in 0:Ns[1]-1
        ei = abs(i - eS[1])
        ej = abs(j - eS[2])
        ek = abs(k - eS[3])
        # if any of i, j, or k is equal to Ns/2 or Ns/2+1
        if ei < 1.0 || ej < 1.0 || ek < 1.0
            push!(edges,ip)
        end
        ip = ip + 1
    end
    end
    end
    return edges
end

function init_grid_R( Ns, LatVecs )
    #
    Npoints = prod(Ns)
    #
    R = Array{Float64}(undef,3,Npoints)
    ip = 0
    for k in 0:Ns[3]-1
    for j in 0:Ns[2]-1
    for i in 0:Ns[1]-1
        ip = ip + 1
        R[1,ip] = LatVecs[1,1]*i/Ns[1] + LatVecs[1,2]*j/Ns[2] + LatVecs[1,3]*k/Ns[3]
        R[2,ip] = LatVecs[2,1]*i/Ns[1] + LatVecs[2,2]*j/Ns[2] + LatVecs[2,3]*k/Ns[3]
        R[3,ip] = LatVecs[3,1]*i/Ns[1] + LatVecs[3,2]*j/Ns[2] + LatVecs[3,3]*k/Ns[3]
    end
    end
    end
    #
    return R
end


text:
include("imports.jl")

@testset ExtendedTestSet "Training" begin

    @testset ExtendedTestSet "fitonecycle!" begin
        learner = testlearner(Recorder())
        @test_nowarn fitonecycle!(learner, 5)
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/PCRE2_jll.jl
module PCRE2_jll

using Libdl

const PATH_list = String[]
const LIBPATH_list = String[]

export libpcre2_8

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libpcre2_8_handle = C_NULL
libpcre2_8_path = ""

if Sys.iswindows()
    const libpcre2_8 = "libpcre2-8-0.dll"
elseif Sys.isapple()
    const libpcre2_8 = "@rpath/libpcre2-8.0.dylib"
else
    const libpcre2_8 = "libpcre2-8.so.0"
end

function __init__()
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = joinpath(Sys.BINDIR, Base.LIBDIR, "julia")
    global libpcre2_8_handle = dlopen(libpcre2_8)
    global libpcre2_8_path = dlpath(libpcre2_8_handle)
end

is_available() = true

end  # module PCRE2_jll

max_stars_repo_name: saberdy/julia

max_issues_repo_name: saberdy/julia

max_forks_repo_name: saberdy/julia


text:
module TestCompression
using DataFrames, RData, Test

rdata_path = joinpath(dirname(@__FILE__), "data_v3")
df = DataFrame(num = [1.1, 2.2])

function test_missing_compressor(filename::AbstractString, codec::Symbol)
    @test_throws CapturedException load(joinpath(rdata_path, filename))
    caught = nothing
    try
        load(joinpath(rdata_path, filename))
    catch ex
        caught = ex
    end

    @test caught isa CapturedException
    @test caught.ex isa CodecMissingError
    @test caught.ex.formatName == codec
end

@testset "Loading compressed RData files" begin
    @testset "When the optional compression codecs are not loaded" begin
        @test_warn "Error encountered while load" begin
            test_missing_compressor("compressed_bzip2.rda", :BZIP2)
        end
        @test_warn "Error encountered while load" begin
            test_missing_compressor("compressed_xz.rda", :XZ)
        end
    end

    @testset "When the optional compression codecs are loaded" begin
        import CodecBzip2, CodecXz

        @test load(joinpath(rdata_path, "compressed_gzip.rda"))["df"] == df
        @test load(joinpath(rdata_path, "compressed_bzip2.rda"))["df"] == df
        @test load(joinpath(rdata_path, "compressed_xz.rda"))["df"] == df
        @test load(joinpath(rdata_path, "compressed_false.rda"))["df"] == df
    end
end

end


text:
using Test
using LogicCircuits
using DataFrames: DataFrame

@testset "Misc utils" begin

    @test issomething("blah")
    @test !issomething(nothing)

    @test isdisjoint(Set([1,3]),Set([5,2]),Set(["a"]))
    @test !isdisjoint(Set([1,"a"]),Set([5,2]),Set(["a"]))

    # repeat a few times for expected code coverage
    [@test length(pushrand!([1,2],3)) == 3 for i=1:100]

    ab = init_array(Bool,5,4)
    @test size(ab) == (5,4)
    @test all(x -> x == 0, ab)

    ab = init_array(Float32,5,4)
    @test size(ab) == (5,4)

    ab = always(Bool,5,4)
    @test size(ab) == (5,4)
    @test all(x -> x == 1, ab)

    ab = always(Float32,5,4)
    @test size(ab) == (5,4)
    @test all(x -> x == 1, ab)

    ab = never(Bool,5,4)
    @test size(ab) == (5,4)
    @test all(x -> x == 0, ab)

    ab = never(Float32,5,4)
    @test size(ab) == (5,4)
    @test all(x -> x == 0, ab)

    un = uniform(5,4)
    @test size(un) == (5,4)
    @test all(x -> x ≈ 1/(5*4), un)

    logun = log.(un)

    lse = logsumexp(logun, 2)
    @test size(lse) == (5,)
    @test all(x -> x ≈ log(1/5), lse)

    lse = logsumexp(logun, 1)
    @test size(lse) == (4,)
    @test all(x -> x ≈ log(1/4), lse)

    lse = logsumexp(logun, (1,2))
    @test size(lse) == ()
    @test all(x -> x ≈ 0, lse)

end

@testset "Imputations and make_missing" begin

    m = [1.1 2.1;
         3.1 4.1;
         5.1 6.1;
         1.0 2.0;
         3.0 4.0;
         2.0 1.0;
         5.0 10.0;
         ]
    df = DataFrame(m, :auto)
    dfb = DataFrame(BitMatrix([true false true;
                               true true true;
                               false true true;
                               true false true;
                               false false false;
                               true false true;
                               false false true;
                            ]), :auto)

    all_missing = make_missing_mcar(df; keep_prob=0.0)
    @test all(ismissing.(Matrix(all_missing)))

    no_missing = make_missing_mcar(df; keep_prob=1.0)
    @test all(.!(ismissing.(Matrix(no_missing))))

    # Imputations Test
    ms_df = make_missing_mcar(df)
    ms_dfb = make_missing_mcar(dfb)

    function test_no_missing()
        @test all(.!(ismissing.(Matrix(imp1))))
        @test all(.!(ismissing.(Matrix(imp2))))
        @test all(.!(ismissing.(Matrix(imp3))))
        @test all(.!(ismissing.(Matrix(imp4))))
        @test all(.!(ismissing.(Matrix(imp5))))
    end

    ## Median Imputation
    imp1 = impute(ms_df)
    imp2 = impute(ms_df, ms_df)
    imp3 = impute(ms_df, df)
    imp4 = impute(ms_dfb, dfb)
    imp5 = impute(ms_dfb)

    test_no_missing()

    ## Mean Imputation
    imp1 = impute(ms_df; method=:mean)
    imp2 = impute(ms_df, ms_df; method=:mean)
    imp3 = impute(ms_df, df; method=:mean)
    imp4 = impute(ms_dfb, dfb; method=:mean)
    imp5 = impute(ms_dfb; method=:mean)

    test_no_missing()

    ## One imputation
    imp1 = impute(ms_df; method=:one)
    imp2 = impute(ms_df, ms_df; method=:one)
    imp3 = impute(ms_df, df; method=:one)
    imp4 = impute(ms_dfb, dfb; method=:one)
    imp5 = impute(ms_dfb; method=:one)

    test_no_missing()
    mask1 = ismissing.(ms_df[:,:])
    @test all(Matrix(imp3)[Matrix(mask1)] .== 1.0)

    ## Zero imputation
    imp1 = impute(ms_df; method=:zero)
    imp2 = impute(ms_df, ms_df; method=:zero)
    imp3 = impute(ms_df, df; method=:zero)
    imp4 = impute(ms_dfb, dfb; method=:zero)
    imp5 = impute(ms_dfb; method=:zero)

    test_no_missing()

    mask1 = ismissing.(ms_df[:,:])
    @test all(Matrix(imp3)[Matrix(mask1)] .== 0.0)

end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
using SortingAlgorithms

# Fast by for multiple function
function fastby(fns::NTuple{N, Function}, byvec::CategoricalVector, valvec::Tuple) where N
    # TODO generalize for categorical
    # TODO can just copy the code from fastby
    res = fastby(fns, byvec.refs, valvec)
    # (byvec.pool.index[res[1].+1], res[2:length(res)]...)
end

# function fastby2(fns::NTuple{N, Function}, byvec::CategoricalVector, valvec::Tuple) where N
#     res = fastby2(fns, byvec.refs, valvec)
#     (byvec.pool.index, res...)
# end

# assumes x is sorted, count the number of uniques
function _ndistinct_sorted(f, x)
    lastx = f(x[1])
    cnt = 1
    @inbounds for i = 2:length(x)
        newx = f(x[i])
        if newx != lastx
            cnt += 1
            lastx = newx
        end
    end
    return cnt
end

# multiple one parameter function for one factors
# one byvec
# function fastby(fns::NTuple{N, Function}, byvec::Vector{T}, valvec::NTuple{N, Vector}) where {N, T}
#     println("miao")
#     l = length(byvec)
#     @time abc = collect(zip(byvec,valvec...))
#     @time sort!(abc, by=x->x[1], alg=RadixSort)

#     @time ucnt = _ndistinct_sorted(x->x[1], abc)

#     # TODO: return a RLE that is easy to traverse and PARALLELIZE
#     lastby::T = abc[1][1]
#     # res = tuple(Vector{typeof(sum(b[1:1]))}(ucnt), Vector{typeof(mean(c[1:1]))}(ucnt))
#     res = ((Vector{typeof(fn(v[1:1]))}(ucnt) for (fn, v) in zip(fns, valvec))...)
#     u_encountered = 1
#     starti = 1
#     @time @inbounds for i in 2:l
#         newby::T = abc[i][1]
#         if newby != lastby
#             # for k = 1:N
#             #     res[k][u_encountered] = fns[k]([abc[j][1+k] for j=starti:i-1])
#             # end
#             u_encountered += 1
#             starti = i
#             lastby = newby
#         end
#     end
#     for k = 1:N
#         res[k][end] = fns[k]([abc[j][1+k] for j=starti:l])
#     end
#     res
# end

# function fastby2(fns::NTuple{N, Function}, byvec::Vector, valvec::Tuple) where N
#     @time abc = collect(zip(valvec...))
#     l = length(byvec)

#     # A number of alternatives are tested
#     # including sorting sort!(byvec, valvec...)
#     # also sort!(byvec, collect(1:length(byvec))) takes about the same amount of
#     # time
#     @time FastGroupBy.grouptwo!(byvec, abc)

#     # TODO: now that it returns a RLE, just PARALLELIZE the summation loop
#     @time res1 = FastGroupBy._contiguousby_vec(sum, byvec, repeat([1], inner = length(byvec)))
#     ucnt = length(res1[1])

#     lastby = byvec[1]
#     # res = tuple(Vector{typeof(sum(b[1:1]))}(ucnt), Vector{typeof(mean(c[1:1]))}(ucnt))
#     res = ((Vector{typeof(fn(v[1:1]))}(ucnt) for (fn, v) in zip(fns, valvec))...)
#     u_encountered = 1
#     starti = 1
#     @time @inbounds for i in 2:l
#         newby = byvec[i]
#         if newby != lastby
#             for k = 1:N
#                 res[k][u_encountered] = fns[k]([abc[j][k] for j=starti:i-1])
#             end
#             u_encountered += 1
#             starti = i
#             lastby = newby
#         end
#     end
#     for k = 1:N
#         res[k][end] = fns[k]([abc[j][k] for j=starti:l])
#     end
#     (res1[1], res...)
# end

# function fastby4(fns::NTuple{N, Function}, byvec::Vector, valvec::NTuple{N, Vector}) where N
#     println("hello")
#     @time abc = collect(zip(valvec...))
#     l = length(byvec)
#
#     # A number of alternatives are tested
#     # including sorting sort!(byvec, valvec...)
#     # also sort!(byvec, collect(1:length(byvec))) takes about the same amount of
#     # time
#     @time FastGroupBy.grouptwo!(byvec, abc)
#
#     # TODO: now that it returns a RLE, just PARALLELIZE the summation loop
#     @time res1 = FastGroupBy._contiguousby_vec(sum, byvec, repeat([1], inner = length(byvec)))
#     ucnt = length(res1[1])
#
#     lastby = byvec[1]
#     # res = tuple(Vector{typeof(sum(b[1:1]))}(ucnt), Vector{typeof(mean(c[1:1]))}(ucnt))
#     res = ((Vector{typeof(fn(v[1:1]))}(ucnt) for (fn, v) in zip(fns, valvec))...)
#     u_encountered = 1
#     starti = 1
#     @time @inbounds for i in 2:l
#         newby = byvec[i]
#         if newby != lastby
#             for k = 1:N
#                 res[k][u_encountered] = fns[k]([abc[j][k] for j=starti:i-1])
#             end
#             u_encountered += 1
#             starti = i
#             lastby = newby
#         end
#     end
#     for k = 1:N
#         res[k][end] = fns[k]([abc[j][k] for j=starti:l])
#     end
#     (res1[1], res...)
# end

# function fastby(fn::NTuple{N, Function}, byvec::Vector{T}, valvec::NTuple{N, Vector}) where {N, T <: Integer}
#     ab = SortingLab.fsortandperm(byvec)
#     orderx = ab[2]
#     # TODO: make a RLE and parallelize the output
#     byby = ab[1]
#
#     # multi-threaded
#     res = Vector{Vector}(N+1)
#     @threads for j=1:length(valvec)
#         vi = valvec[j]
#         @inbounds viv = @view(vi[orderx])
#         @inbounds res1 = FastGroupBy._contiguousby_vec(fn[j], byby, viv)
#         res[j+1] = res1[2]
#         if j == 1
#             res[1] = res1[1]
#         end
#     end
#     res
# end

function fastby!(
    fn::Function,
    byvec::T,
    valvec::AbstractVector{S}
) where {S, T <: Union{PooledArray, CategoricalVector}}

    # println("woohooo")
    l = length(byvec.pool)

    W = valvec[1:1] |> fn |> typeof

    # count the number of occurences of each ref
    # this is the histogram in counting sort
    counter = zeros(UInt, l)
    for r1 in byvec.refs
        @inbounds counter[r1] += 1
    end

    # number of unique groups in the result
    ngrps = sum(counter .!= zero(UInt))

    # Start: Perform counting sort
    lbyvec = length(byvec)
    r1 = byvec.refs[lbyvec]

    # check for degenerate case where there is only sSet of values
    if counter[r1] == lbyvec
        return ([byvec[1]], [fn(valvec)])
    end

    nonzeropos = fcollect(l)[counter .!= 0]

    counter = cumsum(counter)
    rangelo = vcat(0, counter[1:end-1]) .+ 1
    rangehi = copy(counter)

    simvalvec = similar(valvec)

    ci = counter[r1]
    simvalvec[ci] = valvec[lbyvec]
    counter[r1] -= 1

    @inbounds for i = lbyvec-1:-1:1
        r1 = byvec.refs[i]
        ci = counter[r1]
        simvalvec[ci] = valvec[i]
        counter[r1] -= 1
    end
    # Finish: Perform counting sort

    # the result group
    resgrp = copy(@view(byvec[1:ngrps]))
    resval = Vector{W}(undef, ngrps)

    for (i, nzpos) in enumerate(nonzeropos)
        resgrp.refs[i] = nzpos
        (lo, hi) = (rangelo[nzpos], rangehi[nzpos])
        resval[i] = fn(simvalvec[lo:hi])
    end

    return (resgrp, resval)
end


text:
# Utilities for reading access tokens from files.
#
# When using a script you can specify the team, and the appropriate token will be found.
#
# Token are stored in a ".dandelion" in your home directory.

module Util

import DandelionSlack

export
    write_token,
    find_token

token_dir = "~/.dandelion"

token_file(team::AbstractString) = joinpath(expanduser(token_dir), "$(team).token")
find_token(team::AbstractString) = DandelionSlack.Token(utf8(readall(token_file(team))))

ensure_token_dir() = mkpath(expanduser(token_dir))

function write_token(team::AbstractString, token::AbstractString)
    ensure_token_dir()
    open(s -> write(s, token), token_file(team), "a")
end

end


text:
function get_lane_width(veh::Vehicle, roadway::Roadway)

    lane = roadway[veh.state.posF.roadind.tag]

    if n_lanes_left(lane, roadway) > 0
        footpoint = get_footpoint(veh)
        lane_left = roadway[LaneTag(lane.tag.segment, lane.tag.lane + 1)]
        return -proj(footpoint, lane_left, roadway).curveproj.t
    else
        return lane.width
    end
end

function get_markerdist_left(veh::Vehicle, roadway::Roadway)
    t = veh.state.posF.t
    lane_width = get_lane_width(veh, roadway)
    return lane_width/2 - t
end
function get_markerdist_right(veh::Vehicle, roadway::Roadway)
    t = veh.state.posF.t
    lane_width = get_lane_width(veh, roadway)
    return lane_width/2 + t
end


text:
"""
    function make_CMB_T_map(N::Int, pix_size::Float64,
                    ell::AbstractArray{<:Number},
                    DlTT::AbstractArray{<:Number})

# Brief description
Given a CAMB power spectrum return a random CMB
only-temperature map in the`Array{Float64, 2}` format.

# Development
This functions has not presented weird behaviours.
The function return a random CMB temperature anisotropies
maps, that you can draw using the `heatmap` method in
Plots julia's package.

# Args description
N           # pixel number in linear dimension
pix_size    # size of a pixel in arcminutes
ell         # Multipole value - angular scale
DlTT        # Power spectrum value Dₗ = Cₗ⋅(l⋅(l+1)) (μk)²
              Temperature  Temperature
"""
function make_CMB_T_map(N::Int, pix_size::Float64,
                        ell::AbstractArray{<:Number},
                        DlTT::AbstractArray{<:Number})

    ClTT = DlTT .* 2 .* π ./ (ell .* (ell .+ 1.) )
    ClTT[1] = 0.0
    ClTT[2] = 0.0

    X  = reshape(repeat(range(-0.5, 0.5, length=N), N), N, N)
    Y = X'
    R = (X .^2 .+ Y .^2) .^0.5

    pix_to_rad = pix_size / 60. * π  / 180.  # Scale factor between pix to rad
    ell_scale_factor = 2. * π / pix_to_rad

    ell2d = R .* ell_scale_factor
    max_l_R = Int64(floor(maximum(ell2d))+1)

    ClTT_ex = vcat(ClTT, zeros(max_l_R - size(ClTT)[1]))
    indici = trunc.(Int64, ell2d)

    CLTT2d = ClTT_ex[indici]

    rng = MersenneTwister()
    T_rand_arr = randn(rng, N, N)
    fft_T_rand = fft(T_rand_arr)
    FT_sky = (CLTT2d .^0.5) .* fft_T_rand
    T_cmb_map = ifft(ifftshift(FT_sky))
    T_cmb_map = T_cmb_map ./ (pix_size / 60.0 * π / 180.0)

    return real.(T_cmb_map)
end

"""
function make_CMB_pol_maps(N::Int, pix_size::Float64,
                           ell::AbstractArray{<:Number},
                           DlTT::AbstractArray{<:Number},
                           DlEE::AbstractArray{<:Number},
                           DlTE::AbstractArray{<:Number},
                           DlBB::AbstractArray{<:Number})

# Brief description
Given a CAMB power spectrum, in function of the angular scale
(`ell` scale), return a complete set of maps of a hypothetical
random sky. The maps has return concern: Temperature (T_Map),
Q_Stokes_parameter (Q_Map) and U_Stokes_parametes (U_Map).

Whithin the function, the Q and U maps are turn into E and B maps.
On this last two maps you can use the power spectra analysis tools
to extract the BB and EE power_spectrum after the introduction
of the tlescope systematics errors.

# Development
Actually, the function do not work. I have decided to add it
beacouse I think that the issue that I have noticed on this
function has due to CAMB power spectrum format. More test are needed
in particular testing this function with different power-spectrum set.

# Args description
N           # pixel number in linear dimension
pix_size    # size of a pixel in arcminutes
ell         # Multipole value - angular scale
DlTT        # Power spectrum value Dₗ = Cₗ⋅(l⋅(l+1)) (μk)²
              Temperature  Temperature (μk)²
DlEE        # E-modes auto-power_spectrum  polarization  (μk)²
DlTE        # Cross-power_spectrum between Temperature and E modes (μk)²
DlBB        # B-modes auto-power_spectrum (μk)²

"""
function make_CMB_pol_maps(N::Int, pix_size::Float64,
                           ell::AbstractArray{<:Number},
                           DlTT::AbstractArray{<:Number},
                           DlEE::AbstractArray{<:Number},
                           DlTE::AbstractArray{<:Number},
                           DlBB::AbstractArray{<:Number})

    ClTT = DlTT .* 2 .* π ./ (ell .* (ell .+ 1.) )
    ClEE = DlEE .* 2 .* π ./ (ell .* (ell .+ 1.) )
    ClTE = DlTE .* 2 .* π ./ (ell .* (ell .+ 1.) )
    ClBB = DlBB .* 2 .* π ./ (ell .* (ell .+ 1.) )

    ClTT[1:2] .= 0
    ClEE[1:2] .= 0
    ClTE[1:2] .= 0
    ClBB[1:2] .= 0

    # correlation level between T and E-modes
    corr_level_E = ClTE ./ sqrt.(ClTT)
    uncorr_level_E = ClEE .- ( (ClTE .^2.0) ./ ClTT )

    corr_level_E[1:2] .= 0
    uncorr_level_E[1:2] .= 0

    X  = reshape(repeat(range(-0.5, 0.5, length=N), N), N, N)
    Y = X'
    R = (X .^2 .+ Y .^2) .^0.5
    ang = atan.(Y, X)

    pix_to_rad = pix_size / 60. * π  / 180.  # Scale factor between pix to rad
    ell_scale_factor = 2. * π / pix_to_rad

    ell2d = R .* ell_scale_factor
    max_l_R = Int64(floor(maximum(ell2d))+1)

    ClTT_ex = vcat(ClTT, zeros(max_l_R - size(ClTT)[1]))

    ClEE_uncor_ex = vcat(uncorr_level_E,
                         zeros(max_l_R - size(uncorr_level_E)[1]))

    ClEE_corr_ex = vcat(corr_level_E,
                         zeros(max_l_R - size(corr_level_E)[1]))

    ClBB_ex = vcat(ClBB, zeros(max_l_R - size(ClBB)[1]))

    indici = trunc.(Int64, ell2d)
    CLTT2d = ClTT_ex[indici]
    ClEE_uncor2d = ClEE_uncor_ex[indici]
    ClEE_cor2d = ClEE_corr_ex[indici]
    CLBB2d = ClBB_ex[indici]

    rng = MersenneTwister()
    T_rand_arr = randn(rng, N, N)
    E_rand_arr = randn(rng, N, N)
    B_rand_arr = randn(rng, N, N)

    fft_T_rand = fft(T_rand_arr)
    fft_E_rand = fft(E_rand_arr)
    fft_B_rand = fft(B_rand_arr)

    FT_sky = (CLTT2d .^0.5) .* fft_T_rand
    FE_sky = (real.((ClEE_uncor2d .+ 0im) .^0.5)) .* fft_E_rand .+ # Error!
             (real.((ClEE_cor2d .+ 0im) .^0.5)) .* fft_T_rand    # Error!
    FB_sky = (real.((CLBB2d .+ 0im) .^0.5)) .* fft_B_rand          # Error!

    # Conver the E and B maps into Q and U maps
    FQ_sky = (FE_sky .* cos.(2.0 .* ang)) .- (FB_sky .* sin.(2.0 .* ang))
    FU_sky = (FE_sky .* sin.(2.0 .* ang)) .+ (FB_sky .* cos.(2.0 .* ang))

    T_cmb_map = ifft(ifftshift(FT_sky))
    Q_cmb_map = ifft(ifftshift(FQ_sky))
    U_cmb_map = ifft(ifftshift(FU_sky))

    E_map = ifft(ifftshift(FE_sky))
    B_map = ifft(ifftshift(FB_sky))

    T_cmb_map = T_cmb_map ./ (pix_size / 60.0 * π / 180.0)
    Q_cmb_map = Q_cmb_map ./ (pix_size / 60.0 * π / 180.0)
    U_cmb_map = U_cmb_map ./ (pix_size / 60.0 * π / 180.0)
    E_map = E_map ./ (pix_size / 60.0 * π / 180.0)
    B_map = B_map ./ (pix_size / 60.0 * π / 180.0)

    return real.(T_cmb_map), real.(Q_cmb_map), real.(U_cmb_map), real.(E_map), real.(B_map)

end

"""
    function SZ_source(N::Int, pix_size::Float64, Number_of_SZ_Clusters::Int,
                       Mean_Amplitude_of_SZ_Clusters::Float64,
                       SZ_beta::Float64, SZ_Theta_core::Float64)

# Brief description
SZ_source return a sky map Array{Float64, 2} formatted, of the only
SZ sources. The sources are placed random on the map. The pattern is
the well-know SZ-profile given by the `beta_function`. The insensity
and dimensione are normally distribuited around a given mean-amplitude.

# Development
The function has not presented weird behaviours.

# Args description
N                      # pixel number in linear dimension
pix_size               # size of a pixel in arcminutes
Number_of_SZ_Clusters  # Number of SZ sources on the map
Mean_Amplitude         # Mean intensity of those sources
SZ_beta, SZ_Theta_core # Sz spectral indexes

"""
function SZ_source(N::Int, pix_size::Float64, Number_of_SZ_Clusters::Int,
                   Mean_Amplitude_of_SZ_Clusters::Float64,
                   SZ_beta::Float64, SZ_Theta_core::Float64)

    rng = MersenneTwister()
    SZMap = zeros(N, N)
    SZCat = zeros(3, Number_of_SZ_Clusters)
    i = 0
    for i = 1:Number_of_SZ_Clusters
        pix_x = Int64(trunc((N-1)*rand(rng) + 1))
        pix_y = Int64(trunc((N-1)*rand(rng) + 1))
        pix_amplitude = (-1) * Mean_Amplitude_of_SZ_Clusters * randexp(rng)

        SZCat[1, i] += pix_x
        SZCat[2, i] += pix_y
        SZCat[3, i] += pix_amplitude

        SZMap[pix_x, pix_y] += pix_amplitude
    end

    beta = beta_function(N, pix_size, SZ_beta, SZ_Theta_core)
    FT_beta = fft(fftshift(beta))
    FT_SZMap = fft(fftshift(SZMap))
    conv = FT_beta .* FT_SZMap
    SZMap = fftshift(real.(ifft(conv)))

    return SZMap
end

function beta_function(N::Int, pix_size::Float64, SZ_beta::Float64,
                       SZ_Theta_core::Float64)

    inds = ((0:(N-1)) .+ 0.5 .- (N/2)) * pix_size
    X = reshape(repeat(inds, N), N ,N)
    Y = X'
    R = (X .^2 .+ Y .^2) .^0.5
    beta = (1 .+ (R ./ SZ_Theta_core).^2 ) .^ ((1.0 - 3.0 * SZ_beta) / 2.0)

    return beta

end

"""
    function PS_poisson(N::Int, pix_size::Float64, N_Sources::Int,
                        A_Sources::Float64)
# Brief description
Return a sky-map Array{Float64, 2} formatted of the point sources
distribution on the observe sky. You can chose the number of sources
that are places randomly on the map.

This specific function return N_Sources, with amplitude A_Sources.
The amplitude is assigned using a Poisson extraction.
"""
function PS_poisson(N::Int, pix_size::Float64, N_Sources::Int,
                    A_Sources::Float64)
    PS_Pois = zeros(N, N)
    rng = MersenneTwister()
    for i = 1:N_Sources
        pix_x = Int64(trunc((N-1)*rand(rng) + 1))
        pix_y = Int64(trunc((N-1)*rand(rng) + 1))
        PS_Pois[pix_x, pix_y] += A_Sources * randn(rng)
    end

    return PS_Pois

end

"""
    function PS_exp(N::Int, pix_size::Float64, N_Sources::Int,
                    A_Sources::Float64)
# Brief description
Return a sky-map Array{Float64, 2} formatted of the point sources
distribution on the observe sky. You can chose the number of sources
that are places randomly on the map.

This specific function return N_Sources, with amplitude A_Sources.
The amplitude is assigned using a exponential extraction.
"""
function PS_exp(N::Int, pix_size::Float64, N_Sources::Int,
                    A_Sources::Float64)
    PS_expo = zeros(N, N)
    rng = MersenneTwister()
    for i = 1:N_Sources
        pix_x = Int64(trunc((N-1)*rand(rng) + 1))
        pix_y = Int64(trunc((N-1)*rand(rng) + 1))
        PS_expo[pix_x, pix_y] += A_Sources * randexp(rng)
    end

    return PS_expo

end

"""
    function Plot_Sky(N::Int, pix_size::Float64, sky::Array{Float64, 2})

# Brief description
This function plot the CMB heatmap, using the correct axis (in degree)

# Args description
N         Number of pixel (linear dimension)
pix_size  Pixel angular resulotion in arcminute
sky       the sky map in Array{Float64, 2} format
"""
function Plot_Sky(N::Int, pix_size::Float64, sky::Array{Float64, 2};
                  cmin=-400, cmax=400, lab_title="",
                  gradient=ColorGradient([:blue, :white, :red]))
    map_plot = heatmap( sky, c=gradient, xlabel = "x [deg.]",
                    ylabel = "y [deg.]",
                    clims=(cmin, cmax),
                    size=(590,480),
                    xticks = (round.(range(0, N, length=10)),
                              round.(range(0, N*pix_size, length=10) ./ 60.0 )),
                    yticks = (round.(range(0, N, length=10)),
                              round.(range(0, N*pix_size, length=10) ./ 60.0 )),
                    title = lab_title
                  )
    return map_plot

end


text:
#julia fileio
#function
#   isdir,isfile
path = "C:\\Users\\admin\\Desktop\\Scritps\\julia\\Self.jl\\function.jl"
print(isdir(path))
print(isfile(path))
file = open(path)
typeof(file)

close(file)

max_stars_repo_name: zqsheng/snippet

max_issues_repo_name: zqsheng/snippet

max_forks_repo_name: zqsheng/snippet