text:
struct Mul{StyleA, StyleB, AA, BB}
    A::AA
    B::BB
end

@inline Mul{StyleA,StyleB}(A::AA, B::BB) where {StyleA,StyleB,AA,BB} =
    Mul{StyleA,StyleB,AA,BB}(A,B)

@inline Mul(A, B) = Mul{typeof(MemoryLayout(A)), typeof(MemoryLayout(B))}(A, B)

@inline _mul_eltype(A) = A
@inline _mul_eltype(A, B) = Base.promote_op(*, A, B)
@inline _mul_eltype(A, B, C, D...) = _mul_eltype(Base.promote_op(*, A, B), C, D...)

@inline eltype(M::Mul) = _mul_eltype(eltype(M.A), eltype(M.B))

size(M::Mul, p::Int) = size(M)[p]
axes(M::Mul, p::Int) = axes(M)[p]
length(M::Mul) = prod(size(M))
size(M::Mul) = map(length,axes(M))

_mul_axes(A::Tuple{<:Any,<:Any}, ::Tuple{<:Any}) = (A[1],) # matrix * vector
_mul_axes(A::Tuple{<:Any}, B::Tuple{<:Any,<:Any}) = (A[1],B[2]) # vector * matrix
_mul_axes(A::Tuple{<:Any,<:Any}, B::Tuple{<:Any,<:Any}) = (A[1],B[2]) # matrix * matrix
_mul_axes(::Tuple{}, ::Tuple{}) = () # scalar * scalar
_mul_axes(::Tuple{}, B::Tuple) = B # scalar * B
_mul_axes(A::Tuple, B::Tuple{}) = A # A * scalar

axes(M::Mul) = _mul_axes(axes(M.A), axes(M.B))

# The following design is to support QuasiArrays.jl where indices
# may not be `Int`

zeroeltype(M) = zero(eltype(M)) # allow special casing where we know more about zero
zeroeltype(M::Mul{<:Any,<:Any,<:SubArray}) = zeroeltype(Mul(parent(M.A), M.B))

function _getindex(::Type{Tuple{AA}}, M::Mul, (k,)::Tuple{AA}) where AA
    A,B = M.A, M.B
    ret = zeroeltype(M)
    for j = rowsupport(A, k) ∩ colsupport(B,1)
        ret += A[k,j] * B[j]
    end
    ret
end

function _getindex(::Type{Tuple{AA,BB}}, M::Mul, (k, j)::Tuple{AA,BB}) where {AA,BB}
    A,B = M.A,M.B
    ret = zeroeltype(M)
    @inbounds for ℓ in (rowsupport(A,k) ∩ colsupport(B,j))
        ret += A[k,ℓ] * B[ℓ,j]
    end
    ret
end

# linear indexing
_getindex(::Type{NTuple{2,Int}}, M, k::Tuple{Int}) = M[Base._ind2sub(axes(M), k...)...]

_getindex(::Type{Tuple{Int}}, M, (k,)::Tuple{CartesianIndex{1}}) = M[convert(Int, k)]
_getindex(::Type{NTuple{2,Int}}, M, (kj,)::Tuple{CartesianIndex{2}}) = M[kj[1], kj[2]]

"""
   indextype(A)

gives the expected index type for an array, or array-like object. For example,
if it is vector-like it will return `Tuple{Int}`, or if it is matrix-like it will
return `Tuple{Int,Int}`. Other types may have non-integer based indexing.
"""
indextype(A) = indextype(axes(A))

indextype(axes::Tuple) = Tuple{map(eltype, axes)...}

getindex(M::Mul, k...) = _getindex(indextype(M), M, k)

"""
   mulreduce(M::Mul)

returns a lower level lazy multiplication object such as `MulAdd`, `Lmul` or `Rmul`.
The Default is `MulAdd`. Note that the lowered type must overload `copyto!` and `copy`.
"""
mulreduce(M::Mul) = MulAdd(M)

similar(M::Mul, ::Type{T}, axes) where {T,N} = similar(mulreduce(M), T, axes)
similar(M::Mul, ::Type{T}) where T = similar(mulreduce(M), T)
similar(M::Mul) = similar(mulreduce(M))

check_mul_axes(A) = nothing
_check_mul_axes(::Number, ::Number) = nothing
_check_mul_axes(::Number, _) = nothing
_check_mul_axes(_, ::Number) = nothing
_check_mul_axes(A, B) = axes(A,2) == axes(B,1) || throw(DimensionMismatch("Second axis of A, $(axes(A,2)), and first axis of B, $(axes(B,1)) must match"))
function check_mul_axes(A, B, C...)
    _check_mul_axes(A, B)
    check_mul_axes(B, C...)
end

# we need to special case AbstractQ as it allows non-compatiple multiplication
function check_mul_axes(A::AbstractQ, B, C...)
    axes(A.factors, 1) == axes(B, 1) || axes(A.factors, 2) == axes(B, 1) ||
        throw(DimensionMismatch("First axis of B, $(axes(B,1)) must match either axes of A, $(axes(A))"))
    check_mul_axes(B, C...)
end

function instantiate(M::Mul)
    @boundscheck check_mul_axes(M.A, M.B)
    M
end

materialize(M::Mul) = copy(instantiate(M))
@inline mul(A, B) = materialize(Mul(A,B))

copy(M::Mul) = copy(mulreduce(M))
@inline copyto!(dest, M::Mul) = copyto!(dest, mulreduce(M))
@inline copyto!(dest::AbstractArray, M::Mul) = copyto!(dest, mulreduce(M))
mul!(dest::AbstractArray, A::AbstractArray, B::AbstractArray) = copyto!(dest, Mul(A,B))
mul!(dest::AbstractArray, A::AbstractArray, B::AbstractArray, α::Number, β::Number) = muladd!(α, A, B, β, dest)

broadcastable(M::Mul) = M

macro veclayoutmul(Typ)
    ret = quote
        Base.:*(A::AbstractMatrix, B::$Typ) = ArrayLayouts.mul(A,B)
        Base.:*(A::Adjoint{<:Any,<:AbstractMatrix{T}}, B::$Typ{S}) where {T,S} = ArrayLayouts.mul(A,B)
        Base.:*(A::Transpose{<:Any,<:AbstractMatrix{T}}, B::$Typ{S}) where {T,S} = ArrayLayouts.mul(A,B)
        Base.:*(A::LinearAlgebra.AdjointAbsVec, B::$Typ) = ArrayLayouts.mul(A,B)
        Base.:*(A::LinearAlgebra.TransposeAbsVec, B::$Typ) = ArrayLayouts.mul(A,B)
        Base.:*(A::LinearAlgebra.AdjointAbsVec{<:Number}, B::$Typ{<:Number}) = ArrayLayouts.mul(A,B)
        Base.:*(A::LinearAlgebra.TransposeAbsVec{T}, B::$Typ{T}) where T<:Real = ArrayLayouts.mul(A,B)
        Base.:*(A::LinearAlgebra.TransposeAbsVec{T,<:$Typ{T}}, B::$Typ{T}) where T<:Real = ArrayLayouts.mul(A,B)
        Base.:*(A::LinearAlgebra.TransposeAbsVec{T,<:$Typ{T}}, B::AbstractVector{T}) where T<:Real = ArrayLayouts.mul(A,B)

        Base.:*(A::LinearAlgebra.AbstractQ, B::$Typ) = ArrayLayouts.mul(A,B)
    end
    for Struc in (:AbstractTriangular, :Diagonal)
        ret = quote
            $ret

            Base.:*(A::LinearAlgebra.$Struc, B::$Typ) = ArrayLayouts.mul(A,B)
        end
    end
    for Mod in (:AdjointAbsVec, :TransposeAbsVec)
        ret = quote
            $ret

            LinearAlgebra.mul!(dest::AbstractVector, A::$Mod{<:Any,<:$Typ}, b::AbstractVector) =
                ArrayLayouts.mul!(dest,A,b)

            Base.:*(A::$Mod{<:Any,<:$Typ}, B::ArrayLayouts.LayoutMatrix) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::AbstractMatrix) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::AbstractVector) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Number,<:$Typ}, B::AbstractVector{<:Number}) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::$Typ) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Number,<:$Typ}, B::$Typ{<:Number}) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::Diagonal) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::LinearAlgebra.AbstractTriangular) = ArrayLayouts.mul(A,B)
        end
    end

    esc(ret)
end

macro layoutmul(Typ)
    ret = quote
        LinearAlgebra.mul!(dest::AbstractVector, A::$Typ, b::AbstractVector) =
            ArrayLayouts.mul!(dest,A,b)

        LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, B::AbstractMatrix) =
            ArrayLayouts.mul!(dest,A,B)
        LinearAlgebra.mul!(dest::AbstractMatrix, A::AbstractMatrix, B::$Typ) =
            ArrayLayouts.mul!(dest,A,B)
        LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, B::$Typ) =
            ArrayLayouts.mul!(dest,A,B)
        LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, B::AbstractMatrix, α::Number, β::Number) =
            ArrayLayouts.mul!(dest,A,B,α,β)
        LinearAlgebra.mul!(dest::AbstractMatrix, A::AbstractMatrix, B::$Typ, α::Number, β::Number) =
            ArrayLayouts.mul!(dest,A,B,α,β)
        LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, B::$Typ, α::Number, β::Number) =
            ArrayLayouts.mul!(dest,A,B,α,β)

        Base.:*(A::$Typ, B::$Typ) = ArrayLayouts.mul(A,B)
        Base.:*(A::$Typ, B::AbstractMatrix) = ArrayLayouts.mul(A,B)
        Base.:*(A::$Typ, B::AbstractVector) = ArrayLayouts.mul(A,B)
        Base.:*(A::$Typ, B::ArrayLayouts.LayoutVector) = ArrayLayouts.mul(A,B)
        Base.:*(A::AbstractMatrix, B::$Typ) = ArrayLayouts.mul(A,B)
        Base.:*(A::LinearAlgebra.AdjointAbsVec, B::$Typ) = ArrayLayouts.mul(A,B)
        Base.:*(A::LinearAlgebra.TransposeAbsVec, B::$Typ) = ArrayLayouts.mul(A,B)
        Base.:*(A::LinearAlgebra.AdjointAbsVec{<:Any,<:Zeros{<:Any,1}}, B::$Typ) = ArrayLayouts.mul(A,B)
        Base.:*(A::LinearAlgebra.TransposeAbsVec{<:Any,<:Zeros{<:Any,1}}, B::$Typ) = ArrayLayouts.mul(A,B)
        Base.:*(A::LinearAlgebra.Transpose{T,<:$Typ}, B::Zeros{T,1}) where T<:Real = ArrayLayouts.mul(A,B)

        Base.:*(A::LinearAlgebra.AbstractQ, B::$Typ) = ArrayLayouts.mul(A,B)
        Base.:*(A::$Typ, B::LinearAlgebra.AbstractQ) = ArrayLayouts.mul(A,B)
    end
    for Struc in (:AbstractTriangular, :Diagonal)
        ret = quote
            $ret

            Base.:*(A::LinearAlgebra.$Struc, B::$Typ) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Typ, B::LinearAlgebra.$Struc) = ArrayLayouts.mul(A,B)
        end
    end
    for Mod in (:Adjoint, :Transpose, :Symmetric, :Hermitian)
        ret = quote
            $ret

            LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, b::$Mod{<:Any,<:AbstractMatrix}) =
                ArrayLayouts.mul!(dest,A,b)
            LinearAlgebra.mul!(dest::AbstractVector, A::$Mod{<:Any,<:$Typ}, b::AbstractVector) =
                ArrayLayouts.mul!(dest,A,b)
            LinearAlgebra.mul!(dest::AbstractMatrix, A::$Mod{<:Any,<:$Typ}, B::AbstractMatrix, α::Number, β::Number) =
                ArrayLayouts.mul!(dest,A,B,α,β)
            LinearAlgebra.mul!(dest::AbstractMatrix, A::AbstractMatrix, B::$Mod{<:Any,<:$Typ}, α::Number, β::Number) =
                ArrayLayouts.mul!(dest,A,B,α,β)
            LinearAlgebra.mul!(dest::AbstractMatrix, A::$Mod{<:Any,<:$Typ}, B::$Typ, α::Number, β::Number) =
                ArrayLayouts.mul!(dest,A,B,α,β)
            LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, B::$Mod{<:Any,<:$Typ}, α::Number, β::Number) =
                ArrayLayouts.mul!(dest,A,B,α,β)
            LinearAlgebra.mul!(dest::AbstractMatrix, A::$Mod{<:Any,<:AbstractVecOrMat}, B::$Typ, α::Number, β::Number) =
                ArrayLayouts.mul!(dest,A,B,α,β)
            LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, B::$Mod{<:Any,<:AbstractVecOrMat}, α::Number, β::Number) =
                ArrayLayouts.mul!(dest,A,B,α,β)
            LinearAlgebra.mul!(dest::AbstractMatrix, A::$Mod{<:Any,<:$Typ}, B::$Mod{<:Any,<:$Typ}, α::Number, β::Number) =
                ArrayLayouts.mul!(dest,A,B,α,β)

            Base.:*(A::$Mod{<:Any,<:$Typ}, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::AbstractMatrix) = ArrayLayouts.mul(A,B)
            Base.:*(A::AbstractMatrix, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)
            Base.:*(A::LinearAlgebra.AdjointAbsVec, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)
            Base.:*(A::LinearAlgebra.TransposeAbsVec, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::AbstractVector) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::ArrayLayouts.LayoutVector) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::Zeros{<:Any,1}) = ArrayLayouts.mul(A,B)

            Base.:*(A::$Mod{<:Any,<:$Typ}, B::$Typ) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Typ, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)

            Base.:*(A::$Mod{<:Any,<:$Typ}, B::Diagonal) = ArrayLayouts.mul(A,B)
            Base.:*(A::Diagonal, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)

            Base.:*(A::LinearAlgebra.AbstractTriangular, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::LinearAlgebra.AbstractTriangular) = ArrayLayouts.mul(A,B)
        end
    end

    esc(ret)
end

@veclayoutmul LayoutVector

###
# Dot
###
"""
    Dot(A, B)

is a lazy version of `dot(A, B)`, designed to support
materializing based on `MemoryLayout`.
"""
struct Dot{StyleA,StyleB,ATyp,BTyp}
    A::ATyp
    B::BTyp
end

@inline Dot(A::ATyp,B::BTyp) where {ATyp,BTyp} = Dot{typeof(MemoryLayout(ATyp)), typeof(MemoryLayout(BTyp)), ATyp, BTyp}(A, B)
@inline copy(d::Dot) = Base.invoke(LinearAlgebra.dot, Tuple{AbstractArray,AbstractArray}, d.A, d.B)
@inline materialize(d::Dot) = copy(instantiate(d))
@inline Dot(M::Mul{<:DualLayout,<:Any,<:AbstractMatrix,<:AbstractVector}) = Dot(M.A', M.B)
@inline mulreduce(M::Mul{<:DualLayout,<:Any,<:AbstractMatrix,<:AbstractVector}) = Dot(M)
@inline eltype(D::Dot) = promote_type(eltype(D.A), eltype(D.B))

dot(a, b) = materialize(Dot(a, b))
@inline LinearAlgebra.dot(a::LayoutArray, b::LayoutArray) = dot(a,b)
@inline LinearAlgebra.dot(a::LayoutArray, b::AbstractArray) = dot(a,b)
@inline LinearAlgebra.dot(a::AbstractArray, b::LayoutArray) = dot(a,b)
@inline LinearAlgebra.dot(a::LayoutVector, b::AbstractFill{<:Any,1}) = FillArrays._fill_dot(a,b)
@inline LinearAlgebra.dot(a::AbstractFill{<:Any,1}, b::LayoutVector) = FillArrays._fill_dot(a,b)
@inline LinearAlgebra.dot(a::LayoutArray{<:Number}, b::SparseArrays.SparseVectorUnion{<:Number}) = dot(a,b)
@inline LinearAlgebra.dot(a::SparseArrays.SparseVectorUnion{<:Number}, b::LayoutArray{<:Number}) = dot(a,b)

@inline LinearAlgebra.dot(a::SubArray{<:Any,N,<:LayoutArray}, b::AbstractArray) where N = dot(a,b)
@inline LinearAlgebra.dot(a::SubArray{<:Any,N,<:LayoutArray}, b::LayoutArray) where N = dot(a,b)
@inline LinearAlgebra.dot(a::AbstractArray, b::SubArray{<:Any,N,<:LayoutArray}) where N = dot(a,b)
@inline LinearAlgebra.dot(a::LayoutArray, b::SubArray{<:Any,N,<:LayoutArray}) where N = dot(a,b)
@inline LinearAlgebra.dot(a::SubArray{<:Any,N,<:LayoutArray}, b::SubArray{<:Any,N,<:LayoutArray}) where N = dot(a,b)

# Temporary until layout 3-arg dot is added.
# We go to generic fallback as layout-arrays are structured
dot(x, A, y) = dot(x, mul(A, y))
LinearAlgebra.dot(x::AbstractVector, A::LayoutMatrix, y::AbstractVector) = dot(x, A, y)
LinearAlgebra.dot(x::AbstractVector, A::Symmetric{<:Real,<:LayoutMatrix}, y::AbstractVector) = dot(x, A, y)

# allow overloading for infinite or lazy case
@inline _power_by_squaring(_, _, A, p) = Base.invoke(Base.power_by_squaring, Tuple{AbstractMatrix,Integer}, A, p)
# TODO: Remove unnecessary _apply
_apply(_, _, op, Λ::UniformScaling, A::AbstractMatrix) = op(Diagonal(Fill(Λ.λ,size(A,1))), A)
_apply(_, _, op, A::AbstractMatrix, Λ::UniformScaling) = op(A, Diagonal(Fill(Λ.λ,size(A,1))))

for Typ in (:LayoutMatrix, :(Symmetric{<:Any,<:LayoutMatrix}), :(Hermitian{<:Any,<:LayoutMatrix}),
            :(Adjoint{<:Any,<:LayoutMatrix}), :(Transpose{<:Any,<:LayoutMatrix}))
    @eval begin
        @inline Base.power_by_squaring(A::$Typ, p::Integer) = _power_by_squaring(MemoryLayout(A), size(A), A, p)
        @inline +(A::$Typ, Λ::UniformScaling) = _apply(MemoryLayout(A), size(A), +, A, Λ)
        @inline +(Λ::UniformScaling, A::$Typ) = _apply(MemoryLayout(A), size(A), +, Λ, A)
        @inline -(A::$Typ, Λ::UniformScaling) = _apply(MemoryLayout(A), size(A), -, A, Λ)
        @inline -(Λ::UniformScaling, A::$Typ) = _apply(MemoryLayout(A), size(A), -, Λ, A)
    end
end


text:
using CirruParser
using Base.Test
using JSON

testCases = ["comma" "demo" "folding" "html" "indent" "line" "parentheses" "quote" "spaces" "unfolding"]

for name in testCases
  code = readall(open("../examples/$(name).cirru"))
  excepted = strip(readall(open("../ast/$(name).json")))
  ast = CirruParser.pare(code, "")
  formated = strip(JSON.json(ast, 2))
  if formated == excepted
    println("ok: $(name)")
  else
    println("failed:")
    println(formated)
    println("")
  end
  @test formated == excepted
end

println()

ext: jl


text:
# one-column version
function op_H( LF, Potentials, v::Array{Float64,1} )
  #
  V_loc = Potentials.Ionic + Potentials.Hartree + Potentials.XC
  #
  return -0.5*op_nabla2( LF, v ) + V_loc .* v
end

# multi-column version
function op_H( LF, Potentials, v::Array{Float64,2} )
  #
  Npoints = size(v)[1]
  Ncols   = size(v)[2]
  V_loc = Potentials.Ionic + Potentials.Hartree + Potentials.XC
  #
  Hv = zeros( Float64, Npoints, Ncols )
  #
  for ic = 1:Ncols
    Hv[:,ic] = -0.5*op_nabla2( LF, v[:,ic] ) + V_loc .* v[:,ic]
  end
  return Hv
end

# one-column version, sparse
function op_H( LF, ∇2::SparseMatrixCSC{Float64,Int64},
                            Potentials, v::Array{Float64,1} )
  #
  V_loc = Potentials.Ionic + Potentials.Hartree + Potentials.XC
  return -0.5*∇2*v + V_loc.*v
end

# multi-column version, sparse
function op_H( LF, ∇2, Potentials, v::Array{Float64,2} )
  #
  Npoints = size(v)[1]
  Ncols   = size(v)[2]
  V_loc = Potentials.Ionic + Potentials.Hartree + Potentials.XC
  #
  Hv = zeros( Float64, Npoints, Ncols )
  #
  for ic = 1:Ncols
    Hv[:,ic] = -0.5*∇2*v[:,ic] + V_loc.*v[:,ic]
  end
  return Hv
end


text:
module PowerLaws
  using StatsBase
  using Distributions
  using Optim
  using Compat
  import Distributions:rand,shape,pdf,ccdf,scale,params,cdf,cquantile,quantile,fit_mle,mean,median,var,skewness,mode,kurtosis,logccdf,logpdf,logcdf,entropy,@distr_support
  import Base:minimum,maximum,show

  export con_powerlaw,dis_powerlaw,Kolmogorov_smirnov_test,estimate_xmin,bootstrap,bootstrap_p,compare_distributions

  include("helper_macro.jl")
  include("powerlaw_discrete.jl")
  include("powerlaw_continuous.jl")
  include("estimate_xmin.jl")
  include("Kolmogorov-smirnov_test.jl")
  include("bootstrap.jl")
  include("compare_distributions.jl")

end

using Distributions
import PowerLaws:bootstrap,bootstrap_p,con_powerlaw,dis_powerlaw

#this parallel versions of bootstrap,bootstrap_p cannot be defined in module powerlaw at this time because
#of this issue #https://github.com/JuliaLang/julia/issues/13649

function bootstrap(data::AbstractArray,d::UnivariateDistribution, processes::Int64;no_of_sims::Int64 = 10,xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5),seed::Int64 = 0)
  α,θ = params(d)

  n = length(data)

  seed == 0 ? srand() : srand(seed)
  prcs = addprocs(processes)
  rmref = RemoteRef(1)
  put!(rmref,(d,n,xmins,xmax,data))
  helper_arr = Array(Int8,no_of_sims)
  fill!(helper_arr,0)
  @everywhere using PowerLaws
  @everywhere using Distributions
  helper_dict = Dict{Int64,RemoteRef}()
  for prc in workers() # copy data to all processes
    helper_dict[prc] = @spawnat prc begin
      global d,n,xmins,xmax,data
      d,n,xmins,xmax,data = fetch(rmref)
    end
  end
  for prc in workers() # wait till all data is coppied
    fetch(helper_dict[prc])
  end
  helper_dict = 0
  @everywhere function helper_func(args)
    sim_data = sample(data, n, replace=true)
    estimate_xmin(sim_data,typeof(d),xmins = xmins,xmax = xmax)
  end
  statistic = pmap(helper_func,helper_arr)
  rmprocs(prcs)
  take!(rmref)
  return statistic
end

function bootstrap(data::AbstractArray,distribution::Type{con_powerlaw},processes::Int64;no_of_sims::Int64 = 10,xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5),seed::Int64 = 0)
  d,ks = estimate_xmin(data,distribution,xmins = xmins,xmax =xmax)
  bootstrap(data,d,processes,no_of_sims = no_of_sims,xmins = xmins,xmax = xmax,seed =seed)
end

function bootstrap(data::AbstractArray,distribution::Type{dis_powerlaw},processes::Int64;no_of_sims::Int64 = 10,xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5),seed::Int64 = 0)
  d,ks = estimate_xmin(data,distribution,xmins = xmins,xmax =xmax)
  bootstrap(data,d,processes,no_of_sims = no_of_sims,xmins = xmins,xmax = xmax,seed =seed)
end

function bootstrap_p(data::AbstractArray,d::UnivariateDistribution, processes::Int64;no_of_sims::Int64 = 10,xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5),seed::Int64 = 0)
  sort_data = sort(data)
  α,θ = params(d)
  n = length(sort_data)
  tail_indx = findfirst(sort_data,θ)
  tail_p = length(sort_data[tail_indx:end])/n
  KS_stat = Kolmogorov_smirnov_test(sort_data[tail_indx:end],d)

  seed == 0 ? srand() : srand(seed)
  prcs = addprocs(processes)
  rmref = RemoteRef(1)
  put!(rmref,(d,n,tail_p,tail_indx,xmins,xmax,sort_data))
  helper_arr = Array(Int8,no_of_sims)
  fill!(helper_arr,0)
  @everywhere using PowerLaws
  @everywhere using Distributions

  helper_dict = Dict{Int64,RemoteRef}()
  for prc in workers() # copy data to all processes
    helper_dict[prc] = @spawnat prc begin
      global d,n,tail_p,tail_indx,xmins,xmax,sort_data
      d,n,tail_p,tail_indx,xmins,xmax,sort_data = fetch(rmref)
    end
  end
  for prc in workers() # wait till all data is coppied
    fetch(helper_dict[prc])
  end
  helper_dict = 0

  @everywhere helper_func(args) = begin
    n1 = sum(map(x-> x>tail_p,rand(n)))
    n2 = n - n1
    sim_data = Array(Float64,0)
    append!(sim_data,sample(sort_data[1:tail_indx-1],n1,replace = true))
    append!(sim_data,rand(d,n2))
    estimate_xmin(sim_data,typeof(d),xmins = xmins,xmax = xmax)
  end

  statistic= pmap(helper_func,helper_arr)

  rmprocs(prcs)
  take!(rmref)
  P = 0
  for i=1:no_of_sims
    if (KS_stat <= statistic[i][2])
      P +=1
    end
  end

  return statistic,(P/no_of_sims)
end

function bootstrap_p(data::AbstractArray,distribution::Type{con_powerlaw},processes::Int64;no_of_sims::Int64 = 10,xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5),seed::Int64 = 0)
  d,ks = estimate_xmin(data,distribution,xmins = xmins,xmax =xmax)
  bootstrap_p(data,d,processes,no_of_sims = no_of_sims,xmins = xmins,xmax = xmax,seed =seed)
end

function bootstrap_p(data::AbstractArray,distribution::Type{dis_powerlaw},processes::Int64;no_of_sims::Int64 = 10,xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5),seed::Int64 = 0)
  d,ks = estimate_xmin(data,distribution,xmins = xmins,xmax =xmax)
  bootstrap_p(data,d,processes,no_of_sims = no_of_sims,xmins = xmins,xmax = xmax,seed =seed)
end


text:
abstract type PenaltyTerm end

Base.copy(::PT) where PT <: PenaltyTerm = PT()

"""
    Quadratic

Functor object for using a quadratic 2-norm penalty term. Requires an `AbstractMathProgSolver` capable of solving QP problems. Passed by default through `penalty` where applicable.

"""
mutable struct Quadratic <: PenaltyTerm end

function initialize_penaltyterm!(penalty::Quadratic, solver::LQSolver, ::AbstractVector)
    # Get model handle
    model = solver.lqmodel
    # Check if quadratic objectives can be set
    if !applicable(MPB.setquadobj!, model, Int[], Int[], Float64[])
        error("Using a quadratic penalty term requires a solver that can handle quadratic objectives")
    end
    return nothing
end

function update_penaltyterm!(penalty::Quadratic, solver::LQSolver, c::AbstractVector, α::AbstractFloat, ξ::AbstractVector)
    # Get model handle
    model = solver.lqmodel
    # Cache current linear cost
    c_ = copy(c)
    # Linear part
    c_[1:length(ξ)] -= α*ξ
    MPB.setobj!(model, c_)
    # Quadratic part
    qidx = collect(1:length(c))
    qval = fill(α, length(c))
    qval[length(ξ)+1:end] .= zero(α)
    if applicable(MPB.setquadobj!, model, qidx, qidx, qval)
        MPB.setquadobj!(model, qidx, qidx, qval)
    else
        error("Using a quadratic penalty term requires a solver that can handle quadratic objectives")
    end
    return nothing
end

function solve_penalized!(::Quadratic, solver::LQSolver, X::AbstractVector, ::AbstractVector, ::AbstractVector)
    solver(X)
    return nothing
end

"""
    Linearized

Functor object for using an approximately quadratic penalty term, through linearization. Pass through `penalty` where applicable.

...
# Parameters
- `nbreakpoints::Int`: Number of cutting planes used to approximate quadratic term
...
"""
mutable struct Linearized <: PenaltyTerm
    index::Int
    nbreakpoints::Int

    function Linearized(index::Integer, nbreakpoints::Integer)
        n = nbreakpoints >= 3 ? nbreakpoints : 3
        new(index, n)
    end
end
Linearized(; nbreakpoints = 3) = Linearized(-1, nbreakpoints)
Base.copy(linearized::Linearized) = Linearized(-1, linearized.nbreakpoints)

function initialize_penaltyterm!(penalty::Linearized, solver::LQSolver, ξ::AbstractVector)
    # Set constraint index to -1 to indicate that ∞-norm terms have not yet been added
    penalty.index = -1
    # Get model handle
    model = solver.lqmodel
    # Add auxilliary cost variables
    ncols = length(ξ)
    for i in 1:ncols
        MPB.addvar!(model, 0.0, Inf, 1.0)
    end
    return nothing
end

function update_penaltyterm!(penalty::Linearized, solver::LQSolver, c::AbstractVector, α::AbstractFloat, ξ::AbstractVector)
    # Get model handle
    model = solver.lqmodel
    # Cache current linear cost
    c_ = copy(c)
    # Get number of master variables
    ncols = length(ξ)
    # The aux cost variable is the last variable
    tidx = length(c)+1
    # Linear part
    c_[1:length(ξ)] -= α*ξ
    MPB.setobj!(solver.lqmodel, vcat(c_, fill(1.0, length(ξ))))
    # Get current position of linearization constraints
    j = penalty.index
    if j != -1
        # Delete previous constraints (if they exist)
        MPB.delconstrs!(model, collect(j:j+ncols*penalty.nbreakpoints-1))
    end
    x = MPB.getsolution(model)[1:ncols]
    # Add new (or first) set of linearized cost constraints
    for i in 1:ncols
        breakpoints = LinRange(0.0, 2*ξ[i], penalty.nbreakpoints)
        for x in breakpoints
            f = 0.5*α*x^2
            g = α*x
            MPB.addconstr!(model, [i,tidx+(i-1)], [-g,1.], f-g*x, Inf)
        end
    end
    # Update the position of the constraints
    penalty.index = MPB.numconstr(model)-(ncols*penalty.nbreakpoints-1)
    return nothing
end

function solve_penalized!(penalty::Linearized, solver::LQSolver, X::AbstractVector, x::AbstractVector, ξ::AbstractVector)
    push!(X, norm(x-ξ, 2)^2)
    solver(X)
    pop!(X)
    return nothing
end

"""
    InfNorm

Functor object for using a linear ∞-norm penalty term. Pass through `penalty` where applicable.

"""
mutable struct InfNorm <: PenaltyTerm
    index::Int

    InfNorm() = new(-1)
end

function initialize_penaltyterm!(penalty::InfNorm, solver::LQSolver, ::AbstractVector)
    # Set constraint index to -1 to indicate that ∞-norm terms have not yet been added
    penalty.index = -1
    # Get model handle
    model = solver.lqmodel
    # Add ∞-norm auxilliary variable
    MPB.addvar!(model, -Inf, Inf, 1.0)
    return nothing
end

function update_penaltyterm!(penalty::InfNorm, solver::LQSolver, c::AbstractVector, α::AbstractFloat, ξ::AbstractVector)
    # Get model handle
    model = solver.lqmodel
    # Get number of master variables
    ncols = length(ξ)
    # The ∞-norm aux variable is the last variable
    tidx = length(c)+1
    # Linear part
    MPB.setobj!(solver.lqmodel, vcat(c,1.0))
    # Get current position of ∞-norm constraints
    j = penalty.index
    if j != -1
        # Delete previous constraints (if they exist)
        MPB.delconstrs!(model, collect(j:j+2*ncols-1))
    end
    # Add new (or first) set of ∞-norm constraints
    for i in 1:ncols
        MPB.addconstr!(model, [i,tidx], [-α,1], -α*ξ[i], Inf)
        MPB.addconstr!(model, [i,tidx], [-α,-1], -Inf, -α*ξ[i])
    end
    # Update the position of the constraints
    penalty.index = MPB.numconstr(model)-(2*ncols-1)
    return nothing
end

function solve_penalized!(penalty::InfNorm, solver::LQSolver, X::AbstractVector, x::AbstractVector, ξ::AbstractVector)
    push!(X, norm(x-ξ, Inf))
    solver(X)
    pop!(X)
    return nothing
end

"""
    ManhattanNorm

Functor object for using a linear 1-norm penalty term. Pass through `penalty` where applicable.

"""
mutable struct ManhattanNorm <: PenaltyTerm
    index::Int

    ManhattanNorm() = new(-1)
end

function initialize_penaltyterm!(penalty::ManhattanNorm, solver::LQSolver, ξ::AbstractVector)
    # Set constraint index to -1 to indicate that 1-norm terms have not yet been added
    penalty.index = -1
    # Get model handle
    model = solver.lqmodel
    # Add 1-norm auxilliary variables
    ncols = length(ξ)
    for i in 1:ncols
        MPB.addvar!(model, -Inf, Inf, 1.0)
    end
    return nothing
end

function update_penaltyterm!(penalty::ManhattanNorm, solver::LQSolver, c::AbstractVector, α::AbstractFloat, ξ::AbstractVector)
    # Get model handle
    model = solver.lqmodel
    # Get number of master variables
    ncols = length(ξ)
    # The 1-norm aux variables are the last variables
    tidx = length(c)+1
    # Linear part
    MPB.setobj!(solver.lqmodel, vcat(c, fill(1.0, length(ξ))))
    # Get current position of ∞-norm constraints
    j = penalty.index
    if j != -1
        # Delete previous constraints (if they exist)
        MPB.delconstrs!(model, collect(j:j+2*ncols-1))
    end
    # Add new (or first) set of 1-norm constraints
    for i in 1:ncols
        MPB.addconstr!(model, [i,tidx+(i-1)], [-α,1], -α*ξ[i], Inf)
        MPB.addconstr!(model, [i,tidx+(i-1)], [-α,-1], -Inf, -α*ξ[i])
    end
    # Update the position of the constraints
    penalty.index = MPB.numconstr(model)-(2*ncols-1)
    return nothing
end

function solve_penalized!(penalty::ManhattanNorm, solver::LQSolver, X::AbstractVector, x::AbstractVector, ξ::AbstractVector)
    ts = map(abs, x .- ξ)
    append!(X, ts)
    solver(X)
    [pop!(X) for _ in ts]
    return nothing
end


text:
include("sod_shock.jl")
#include("cr_sod_shock_main.jl")
include("cr_sod_shock_noprepopulation.jl")
include("cr_sod_shock_withprepopulation.jl")

"""
    RiemannParameters(;rhol::Float64=1.0, rhor::Float64=0.125,      # density left and right (L&R)
                       Pl::Float64=0.0,   Pr::Float64=0.0,          # pressure L&R
                       Ul::Float64=0.0,   Ur::Float64=0.0,          # internal energy L&R
                       P_cr_l::Float64=0.0, P_cr_r::Float64=0.0,    # CR pressure L&R
                       E_cr_l::Float64=0.0, E_cr_r::Float64=0.0,    # CR energy L&R
                       Bl::Array{Float64,1} = zeros(3),             # B-field left
                       Br::Array{Float64,1} = zeros(3),             # B-field right
                       Mach::Float64=0.0,                           # target Mach number
                       t::Float64,                                  # time of the solution
                       x_contact::Float64=70.0,                     # position of the contact discontinuity along the tube
                       γ_th::Float64=5.0/3.0,                       # adiabatic index of the gas
                       γ_cr::Float64=4.0/3.0,                       # adiabatic index of CRs
                       Pe_ratio::Float64=0.01,                      # ratio of proton to electron energy in acceleration
                       thetaB::Float64=0.0,                         # angle between magnetic field and shock normal
                       theta_crit::Float64=(π/4.0),                 # critical angle for B/Shock angle efficiency
                       dsa_model::Int64=-1,                         # diffuse shock acceleration model
                       xs_first_guess::Float64=4.7,                 # first guess of the resulting shock compression
                       verbose::Bool=true)

Wrapper function to set up parameter objects for different kinds of Riemann problems.
"""
function RiemannParameters(;rhol::Float64=1.0, rhor::Float64=0.125,      # density left and right (L&R)
                            Pl::Float64=0.0,   Pr::Float64=0.0,          # pressure L&R
                            Ul::Float64=0.0,   Ur::Float64=0.0,          # internal energy L&R
                            P_cr_l::Float64=0.0, P_cr_r::Float64=0.0,    # CR pressure L&R
                            E_cr_l::Float64=0.0, E_cr_r::Float64=0.0,    # CR energy L&R
                            Bl::Array{Float64,1} = zeros(3),             # B-field left
                            Br::Array{Float64,1} = zeros(3),             # B-field right
                            Mach::Float64=0.0,                           # target Mach number
                            t::Float64,                                  # time of the solution
                            x_contact::Float64=70.0,                     # position of the contact discontinuity along the tube
                            γ_th::Float64=5.0/3.0,                       # adiabatic index of the gas
                            γ_cr::Float64=4.0/3.0,                       # adiabatic index of CRs
                            Pe_ratio::Float64=0.01,                      # ratio of proton to electron energy in acceleration
                            thetaB::Float64=0.0,                         # angle between magnetic field and shock normal
                            theta_crit::Float64=(π/4.0),                 # critical angle for B/Shock angle efficiency
                            dsa_model::Int64=-1,                         # diffuse shock acceleration model
                            xs_first_guess::Float64=4.7,                 # first guess of the resulting shock compression
                            verbose::Bool=true)

    # Error handling
    if (Pl == 0.0 && Ul == 0.0)
        error("Both Ul and Pl are zero! Plase supply initital pressure values.")
    end # error handling

    # pure hydro case -> no pre-existing CRs and no DSA model selected
    if ( (P_cr_l == 0.0 && P_cr_r == 0.0) ||
         (E_cr_l == 0.0 && E_cr_r == 0.0) ) &&
         dsa_model == -1

         if verbose
             @info "Setting up parameters for pure hydro Sod-shock."
         end

         return SodParameters(rhol=rhol, rhor=rhor, Pl=Pl, Pr=Pr,
                              Ul=Ul, Ur=Ur, Mach=Mach, t=t,
                              x_contact=x_contact, γ_th=γ_th)

    end

    if dsa_model != -1

        if verbose
            @info "Setting up parameters for Sod-shock with CR acceleration."
        end

        if ( (P_cr_l == 0.0 && P_cr_r == 0.0) &&
             (E_cr_l == 0.0 && E_cr_r == 0.0) )

            if verbose
                @info "No seed CRs."
            end

            return SodCRParameters_noCRs(rhol=rhol, rhor=rhor,
                                         Pl=Pl, Pr=Pr, Ul=Ul, Ur=Ur,
                                         Mach=Mach, t=t, x_contact=x_contact,
                                         Pe_ratio=Pe_ratio, γ_th=γ_th, γ_cr=γ_cr,
                                         thetaB=thetaB, theta_crit=theta_crit,
                                         dsa_model=dsa_model,
                                         xs_first_guess=xs_first_guess)
        else
            if verbose
                @info "With seed CRs."
            end

            #error("Sod shock with seed CRs not implemented yet!")
            return SodCRParameters_withCRs(rhol=rhol, rhor=rhor,
                                           Pl=Pl, Pr=Pr, Ul=Ul, Ur=Ur,
                                           P_cr_l=P_cr_l, P_cr_r=P_cr_r,
                                           E_cr_l=E_cr_l, E_cr_r=E_cr_r,
                                           Mach=Mach, t=t, x_contact=x_contact,
                                           Pe_ratio=Pe_ratio, γ_th=γ_th, γ_cr=γ_cr,
                                           thetaB=thetaB, theta_crit=theta_crit,
                                           dsa_model=dsa_model)
        end #
    end # dsa_model != -1

    if ( (P_cr_l != 0.0 && P_cr_r != 0.0) ||
         (E_cr_l != 0.0 && E_cr_r != 0.0) ) &&
         dsa_model == -1

         if verbose
             @info "Setting up parameters for multicomponent shock without CR acceleration."
         end

         #error("Multicomponent fluid shock not implemented yet!")
         return SodCRParameters_withCRs(rhol=rhol, rhor=rhor, Pl=Pl, Pr=Pr,
                                         Ul=Ul, Ur=Ur, P_cr_l=P_cr_l, P_cr_r=P_cr_r,
                                         E_cr_l=E_cr_l, E_cr_r=E_cr_r,
                                         Mach=Mach, t=t, x_contact=x_contact,
                                         Pe_ratio=Pe_ratio, γ_th=γ_th, γ_cr=γ_cr,
                                         dsa_model=dsa_model)

    end # Seed CRs without acc

end # RiemannParameters

"""
    find_xs_first_guess(Ul::Float64, Mach::Float64, CR_seed::Float64=0.0;
                        xs_start::Float64=3.8, delta_xs::Float64=1.e-4,
                        eff_model::Int64=2, thetaB::Float64=0.0,
                        verbose::Bool=false)

Iterates to a first guess of the compression ratio `xs` that gives a solution.
"""
function find_xs_first_guess(Ul::Float64, Mach::Float64, CR_seed::Float64=0.0;
                             xs_start::Float64=3.8, delta_xs::Float64=1.e-4,
                             eff_model::Int64=2, thetaB::Float64=0.0,
                             verbose::Bool=false)

    xs_first_guess = xs_start
    p = RiemannParameters(Ul=Ul, Mach=Mach,
                          t=1.5)

    while true
        if verbose
            println(xs_first_guess)
        end

        try
            P_cr_l = CR_seed * p.Pl
            P_cr_r = CR_seed * p.Pr

            par = RiemannParameters(Pl=p.Pl, Pr=p.Pr,
                                    P_cr_l=P_cr_l, P_cr_r=P_cr_r,
                                    thetaB = thetaB, Mach = Mach,
                                    xs_first_guess = xs_first_guess,
                                    t=1.5, dsa_model = eff_model,
                                    verbose = verbose)

            solve([86.0], par)

            break
        catch
            xs_first_guess += delta_xs
        end

    end
    return xs_first_guess
end

"""
    Multiple dispatch for solve function
"""
# Pure hydro Sod-shock
"""
    solve(x::Array{Float,1}, par::SodParameters)

Solves a standard Sod shock problem.
"""
solve(x::Array{Float64,1}, par::SodParameters) = solveSodShock(x, par=par)
solve(x::Array{Float32,1}, par::SodParameters) = solveSodShock(Float64.(x), par=par)

# CR Sod shock
"""
    solve(x::Array{Float,1}, par::SodCRParameters_noCRs)

Solves a Sod shock with cosmic ray acceleration, but without a pre-existing CR component.
"""
solve(x::Array{Float64,1}, par::SodCRParameters_noCRs)   = solveSodShockCR_noPrepopulation(x, par=par)
solve(x::Array{Float32,1}, par::SodCRParameters_noCRs)   = solveSodShockCR_noPrepopulation(Float64.(x), par=par)

"""
    solve(x::Array{Float,1}, par::SodCRParameters_withCRs)

Solves a Sod shock with cosmic ray acceleration without a pre-existing CR component.
"""
solve(x::Array{Float64,1}, par::SodCRParameters_withCRs) = solveSodShockCR_withPrepopulation(x, par=par)
solve(x::Array{Float32,1}, par::SodCRParameters_withCRs) = solveSodShockCR_withPrepopulation(Float64.(x), par=par)

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
using SymEngine
export FO,PFO,dag
export Up,Dn,Sl

abstract type AbstractFO end

"""
    Up spin
"""
abstract type Up end

"""
    Down spin
"""
abstract type Dn end

"""
    Spinless
"""
abstract type Sl end
const Spin = Union{Type{Up},Type{Sl},Type{Dn}}

"""
struct FO
    This package assume that each fermionic operator consist of such several parts.
    Example of Fermionic Operator: c_4_(up)^dagger
    1. Symbol. Example: c
    2. Site. Example: 4
    3. Spin. Example: up/dn
    4. Dag. Example: true/false
"""
struct FO <: AbstractFO
    symbol::Symbol
    site::Int
    spin::Spin
    dag::Bool
end

"""
struct PFO
    Powered Ferminoic Operators.
    An abstract FO wich indicate a FO with power.
    Each FO is binded with a abstract symbol to indicate its power.
"""
struct PFO <: AbstractFO
    fo::FO
    pow::Basic
end

export dag
dag(fo::FO) = FO(fo.symbol,fo.site,fo.spin,!(fo.dag))
dag(pfo::PFO) = PFO(dag(pfo.fo),pfo.pow)


text:
abstract type AvgFidelityMetric

end

convert(t::Type{T},x::AvgFidelityMetric) where {T<:Real} = x.val
promote_rule(t::Type{T},r::Type{R}) where {T<:AvgFidelityMetric,R<:Real} = promote_rule(Float64,R)

"""
AvgFidelity(f;dim=2)
"""
struct AvgFidelity <: AvgFidelityMetric
  dim::Int64
  val::Float64
  AvgFidelity(f;dim=2) = new(dim,f)
end

avgfidelity(x::AvgFidelity) = x.val
convert(t::Type{AvgFidelity},f::AvgFidelity) = AvgFidelity(f.val,dim=f.dim)
convert(t::Type{Float64},f::AvgFidelity) = f.val

"""
DepolRate(p;dim=2)

Depolarizing rate of unitarily invariant twirl of some random channel.
If the unitarily invariant twirl of some channel is `D`, the depolarizing rate `p` is
defined by

`D(ρ) = p ⋅ ρ + (1-p) ⋅ I/d`

where `I` is the identity matrix. Note that this implies `p` may be negative (i.e., it
is not a probability).
"""
struct DepolRate <: AvgFidelityMetric
  dim::Int64
  val::Float64
  DepolRate(p;dim=2) = new(dim,p)
end

avgfidelity(x::DepolRate) = ((x.dim-1)*x.val+1)/x.dim
convert(t::Type{DepolRate},f::AvgFidelity) = DepolRate((f.dim*f.val-1)/(f.dim-1),dim=f.dim)

struct RBDecayRate <: AvgFidelityMetric
  dim::Int64
  val::Float64
  RBDecayRate(p;dim=2) = new(dim,p)
end

avgfidelity(x::RBDecayRate) = 1-x.val
convert(t::Type{RBDecayRate},f::AvgFidelity) = RBDecayRate(1-f.val,dim=f.dim)

struct EntanglementFidelity <: AvgFidelityMetric
  dim::Int64
  val::Float64
  EntanglementFidelity(e; dim=2) = new(dim,e)
end

avgfidelity(x::EntanglementFidelity) = (x.dim * x.val + 1)/(x.dim + 1)
convert(t::Type{EntanglementFidelity},f::AvgFidelity) = EntanglementFidelity((f.val*(f.dim+1)-1)/f.dim,dim=f.dim)

struct χ00 <: AvgFidelityMetric
  dim::Int64
  val::Float64
  χ00(e; dim=2) = new(dim,e)
end

avgfidelity(x::χ00) = (x.dim * x.val + 1)/(x.dim + 1)
convert(t::Type{χ00},f::AvgFidelity) = χ00((f.val*(f.dim+1)-1)/f.dim,dim=f.dim)

# convert{T1,T2}(t::Type{T1},f::T2) = convert(T1,avgfidelity(,dim=f.dim))


text:
using GLOQ
using Test
using Random
using Plots
using DifferentialEquations
using LinearAlgebra
pyplot()
#include("../src/GLOQ.jl")

@testset "GLOQ.jl" begin
    GLOQ.hello_world()
    include("forward_solve_constant_lindblad.jl")
    # Write your tests here.
end


text:
# TODO: All of this logging stuff might be better handled with multiple dispatch on parametric functions.

# add! a whole structure to the log.
function add!(log::HDF5Logger.Log, slug::String, t::Real, data, n::Int, keep::Bool = false)
    function add_structured!(log, slug, data::T, n, keep) where {T}
        if typeof(data) <: Enum
            # println("I found an enum for " * slug)
            HDF5Logger.add!(log, slug, Base.Enums.basetype(typeof(data))(data), n, keep) # Convert to the base type for the enum.
        elseif isa(data, Tuple)
            @assert isbitstype(eltype(data)) "HDF5Logger doesn't know what to do with a tuple of non-primitive types."
            # println("I found a tuple type for " * slug)
            HDF5Logger.add!(log, slug, [data...,], n, keep) # Convert to an array.
        elseif isa(data, Array)
            @assert isbitstype(eltype(data)) "HDF5Logger doesn't know what to do with an array of non-primitize types."
            # println("I found an array for " * slug)
            HDF5Logger.add!(log, slug, data, n, keep) # HDF5Loggers only know scalars, vectors, and matrices.
        elseif isstructtype(T) # Tuples have to come before this because tuples are struct types?? So are arrays???
            # println("I found a struct type for " * slug)
            for field in fieldnames(typeof(data))
                add_structured!(log, slug * "/" * string(field), getfield(data, field), n, keep)
            end
        elseif isbits(data)
            # println("I found either a bits type or array for " * slug)
            HDF5Logger.add!(log, slug, data, n, keep) # HDF5Loggers only know scalars, vectors, and matrices.
        else
            error("HDF5Logger doesn't know what to do with a " * slug)
        end
    end
    HDF5Logger.add!(log, slug * "time", t, n, keep)
    add_structured!(log, slug * "data", data, n, keep)
end

# log! a whole structure.
function log!(log::HDF5Logger.Log, slug::String, t::Real, data)
    function log_structured!(log, slug, data::T) where {T}
        if typeof(data) <: Enum
            HDF5Logger.log!(log, slug, Base.Enums.basetype(typeof(data))(data)) # Convert to the base type for the enum.
        elseif isa(data, Tuple)
            @assert isbitstype(eltype(data)) "HDF5Logger doesn't know what to do with a tuple of non-primitive types."
            HDF5Logger.log!(log, slug, [data...,]) # Convert to an array.
        elseif isa(data, Array)
            @assert isbitstype(eltype(data)) "HDF5Logger doesn't know what to do with an array of non-primitize types."
            HDF5Logger.log!(log, slug, data)
        elseif isstructtype(T)
            for field in fieldnames(typeof(data))
                log_structured!(log, slug * "/" * string(field), getfield(data, field))
            end
        elseif isbits(data) || (isa(data, Array) && isbitstype(eltype(data)))
            HDF5Logger.log!(log, slug, data)
        else
            error("HDF5Logger doesn't know what to do with a " * slug)
        end
    end
    HDF5Logger.log!(log, slug * "time", t)
    log_structured!(log, slug * "data", data)
end


text:
# -------------------------------- #
# Custom JSON output for our types #
# -------------------------------- #
JSON.lower(a::HasFields) = a.fields

function _apply_style_axis!(p::Plot, ax, force::Bool=false)
    if haskey(p.style.layout.fields, Symbol(ax, "axis")) || force
        ax_names = Iterators.filter(
            _x-> startswith(string(_x), "$(ax)axis"),
            keys(p.layout.fields)
        )

        for ax_name in ax_names
            cur = p.layout.fields[ax_name]
            cur = merge(p.style.layout[Symbol(ax, "axis")], cur)
            p.layout.fields[ax_name] = cur
        end

        if isempty(ax_names)
            nm = Symbol(ax, "axis")
            p.layout.fields[nm] = deepcopy(p.style.layout[nm])
        end
    end

end

_maybe_set_attr!(hf::HasFields, k::Symbol, v::Any) =
    get(hf, k, nothing) == nothing && setindex!(hf, v, k)

# special case for associative to get nested application
function _maybe_set_attr!(hf::HasFields, k1::Symbol, v::Associative)
    for (k2, v2) in v
        _maybe_set_attr!(hf, Symbol(k1, "_", k2), v2)
    end
end

function _maybe_set_attr!(p::Plot, k1::Symbol, v::Associative)
    for (k2, v2) in v
        _maybe_set_attr!(p, Symbol(k1, "_", k2), v2)
    end
end

function _maybe_set_attr!(p::Plot, k::Symbol, v)
    foreach(t -> _maybe_set_attr!(t, k, v), p.data)
end

function _maybe_set_attr!(p::Plot, k::Symbol, v::Cycler)
    ix = 0
    for t in p.data
        if t[k] == Dict()  # was empty
            t[k] = v[ix+=1]
        end
    end
end

function JSON.lower(p::Plot)
    _is3d = any(_x -> contains(string(_x[:type]), "3d"), p.data)

    # apply layout attrs
    if !isempty(p.style.layout)
        # force xaxis and yaxis if plot is 2d
        _apply_style_axis!(p, "x", !_is3d)
        _apply_style_axis!(p, "y", !_is3d)
        _apply_style_axis!(p, "z")

        # extract this so we can pop! off xaxis and yaxis so they aren't
        # applied again
        la = deepcopy(p.style.layout)
        pop!(la.fields, :xaxis, nothing)
        pop!(la.fields, :yaxis, nothing)
        pop!(la.fields, :zaxis, nothing)
        p.layout = merge(la, p.layout)
    end

    # apply global trace attrs
    if !isempty(p.style.global_trace)
        for (k, v) in p.style.global_trace.fields
            _maybe_set_attr!(p, k, v)
        end
    end

    # apply trace specific attrs
    if !isempty(p.style.trace)
        for t in p.data
            t_type = Symbol(get(t, :type, :scatter))
            for (k, v) in get(p.style.trace, t_type, Dict())
                _maybe_set_attr!(p, k, v)
            end
        end
    end
    Dict(:data => p.data, :layout => p.layout)
end

JSON.lower(sp::SyncPlot) = JSON.lower(sp.plot)
JSON.lower(a::Colors.Colorant) = string("#", hex(a))

# Let string interpolation stringify to JSON format
Base.print(io::IO, a::Union{Shape,GenericTrace,PlotlyAttribute,Layout,Plot}) = print(io, JSON.json(a))
Base.print{T<:GenericTrace}(io::IO, a::Vector{T}) = print(io, JSON.json(a))

GenericTrace(d::Associative{Symbol}) = GenericTrace(pop!(d, :type, "scatter"), d)
GenericTrace{T<:AbstractString}(d::Associative{T}) = GenericTrace(_symbol_dict(d))
Layout{T<:AbstractString}(d::Associative{T}) = Layout(_symbol_dict(d))

function JSON.parse(::Type{Plot}, str::AbstractString)
    d = JSON.parse(str)
    data = GenericTrace[GenericTrace(tr) for tr in d["data"]]
    layout = Layout(d["layout"])
    Plot(data, layout)
end

JSON.parsefile(::Type{Plot}, fn) =
    open(fn, "r") do f; JSON.parse(Plot, readstring(f)) end


text:
include("../Resistances/heat.jl")
include("../../Utilities/summary.jl")

using Measurements, Modia
#ENV["MODIA_PLOT"] = "PyPlot"
Modia.usePlotPackage("PyPlot")
@usingModiaPlot

# Boundary model(s)
T_a = 95.6u"K"
T_b = 273.15u"K"

# Parameters
r_inner = 0.1u"m"#measurement(0.1,0.01)u"m"
th_1 = 0.0025u"m"
th_ins = 0.01u"m"
th_2 = 0.0025u"m"

r_ins_inner = r_inner + th_1
r_ins_outer = r_ins_inner + th_ins
r_outer = r_ins_outer + th_2

model = Model(
    boundaryT_a = Temperature | Map(T=T_a),
    boundaryT_b = Temperature | Map(T=T_b),
    convectionInner = Convection | Map(surfaceArea=4*π*r_inner^2,alpha=150u"W/(m^2*K)"),
    linerInner = Sphere | Map(r_inner=r_inner,r_outer=r_ins_inner,lambda=15u"W/(m*K)"),
    contact_1 = Contact | Map(surfaceArea=4*π*r_ins_inner^2,Rc_pp=0.003u"m^2*K/W"),
    insulation = Sphere | Map(r_inner=r_ins_inner,r_outer=r_ins_outer,lambda=0.033u"W/(m*K)"),
    contact_2 = Contact | Map(surfaceArea=4*π*r_ins_outer^2,Rc_pp=0.003u"m^2*K/W"),
    linerOuter = Sphere | Map(r_inner=r_ins_outer,r_outer=r_outer,lambda=15u"W/(m*K)"),
    radiationOuter = Radiation | Map(surfaceArea=4*π*r_outer^2,epsilon=0.7u"m/m"), #nondim units?
    convectionOuter = Convection | Map(surfaceArea=4*π*r_outer^2,alpha=6.0u"W/(m^2*K)"),
    Q_total=Var(init=0.0u"J"),
    equations =:[
        der(Q_total) = convectionInner.port_a.Q_flow
    ],
    connect =:[
        (convectionInner.port_a,boundaryT_a.port)
        (convectionInner.port_b,linerInner.port_a)
        (linerInner.port_b,contact_1.port_a)
        (contact_1.port_b,insulation.port_a)
        (insulation.port_b,contact_2.port_a)
        (contact_2.port_b,linerOuter.port_a)
        (linerOuter.port_b,radiationOuter.port_a,convectionOuter.port_a)
        (convectionOuter.port_b,radiationOuter.port_b,boundaryT_b.port)
    ]
)

result = @instantiateModel(model, log=false,FloatType=Measurement{Float64})
@time simulate!(result, stopTime=1.0, tolerance=1e-4, log=false)
y_min, t_min, y_max, t_max = reportMinMax("Q_total","Thermal Resistances Model",result)
errorCheck(-69.43455u"J"±0.0u"J", y_min)
plot(result, ["Q_total"], figure=1)


text:
module day11

include("inputs.jl")

ex = inputs.ex_d11()
ex1 = inputs.ex1_d11()
ex2 = inputs.ex2_d11()
input = inputs.input_d11()

build_array(input) = collect.(input) |> x->hcat(x...) |> permutedims

function king(idx, grid, dirs)
    """ Neighbors as the King could reach."""
    return map(dirs) do dir
        checkbounds(Bool, grid, dir + idx) ? grid[idx+dir] : '.'
    end
end

function queen(idx, grid, dirs)
    """ Neighbors as the Queen could reach."""
    seat = idx # needed to reset the index after each inner loop
    neighbors = map(dirs) do dir
        idx = seat
        while checkbounds(Bool, grid, dir+idx)
            grid[dir+idx] ∈ ('L', '#') && return grid[dir+idx]
            idx += dir
        end
        return '.'
    end
end

function step(grid, get_neighbors, seen)
    next_grid = copy(grid)
    dirs = [CartesianIndex(i,j) for i ∈ -1:1, j ∈ -1:1 if !(i==j==0)]
    for idx in CartesianIndices(grid)
        neighbors = get_neighbors(idx, grid, dirs)
        if grid[idx] == 'L' && ('#' ∉ neighbors)
            next_grid[idx] = '#'
        elseif grid[idx] == '#' && count(x->x=='#', neighbors) >= seen
            next_grid[idx] = 'L'
        end
    end
    return next_grid
end

function count_occupied(input, get_neighbors, seen)
    grid = build_array(input)
    grid_ = similar(grid)

    while !(grid_ == grid)
        grid_ = copy(grid)
        grid = step(grid_, get_neighbors, seen)
    end
    return count(x->x=='#', grid)
end

part1(input) = count_occupied(input, king, 4)
part2(input) = count_occupied(input, queen, 5)

@assert part1(ex) == 37
@assert part1(input) == 2273
@info "Day 11, Part 1 answer: $(part1(input))"

@assert part2(ex) == 26
@assert part2(input) == 2064
@info "Day 11, Part 2 answer: $(part2(input))"

end


text:
class Map {
    field table, size
    def new(size) {
        self.table = []
        self.size = size
        self.table:(size-1) <- []
    }

    def map_hash(key) hash_uint32_t_(key) % self.size

    def get(key) {
        hv = self.map_hash(key)
        if ~self.table[hv] return False

        for i=0, i < |self.table[hv]|, i=i+1
            if key == self.table[hv][i][0]
                return self.table[hv][i][1]

        return False
    }

    def put(key, value) {
        hv = self.map_hash(key)
        if self.table[hv] {
            for i=0, i < |self.table[hv]|, i=i+1
                if key == self.table[hv][i][0] {
                    self.table[hv][i][1] = value
                    return
                }
            self.table[hv]:0 <- (key, value)
        } else
            self.table[hv] = [(key, value)]
    }
}

max_stars_repo_name: jeffreymanzione/jeff-vm-lang

max_issues_repo_name: jeffreymanzione/jeff-vm-lang

max_forks_repo_name: jeffreymanzione/jeff-vm-lang


text:
using LinearAlgebra
using QuantumSpin

@doc raw"""
    spin_spin_interaction_term(A, σ)

Compute a spin-spin interaction term:

``\sum_{i} h_{i} \sigma_i^\alpha``,

where ``h_{i}`` is a vector and ``\alpha = x,y,z``.
"""
function on_site_term(h, σ)
    N = length(h)
    Res = spzeros(2^N, 2^N)
    for index in findall(!iszero, h)
        Res += ⊗(set_spins(N, index, σ)...) * h[index]
    end
    Res
end

@doc raw"""
    spin_spin_interaction_term(A, σ)

Compute a spin-spin interaction term:

``\sum_{i>j} A_{ij} \sigma_i^\alpha \sigma_j^\alpha``,

where ``A_{ij}`` is an adjacency matrix and ``\alpha = x,y,z``.
"""
function spin_spin_interaction_term(A, σ)
    @assert issymmetric(A) "The adjacency matrix must be symmetric."
    N = size(A)[1]
    Res = spzeros(2^N, 2^N)
    for indeces in findall(!iszero, tril(A))
        Res += ⊗(set_spins(N, indeces, [σ, σ])...) *A[indeces]
    end
    Res
end

@doc raw"""
    hamiltonian(Jˣ, Jʸ, Jᶻ, Γˣ, Γʸ, Γᶻ)

Compute a Hamiltonian:

``H = \sum_{i>j} \left(
    J^x_{ij} \sigma_i^x \sigma_j^x
    + J^y_{ij} \sigma_i^y \sigma_j^y
    + J^z_{ij} \sigma_i^z \sigma_j^z
    \right)
    +
    \sum_{i} \left(
    Γ^x_i \sigma_i^x
    + Γ^y_i \sigma_i^y
    + Γ^z_i \sigma_i^z
    \right).
    ``
"""
function hamiltonian(Jˣ, Jʸ, Jᶻ, Γˣ, Γʸ, Γᶻ)
    H = spin_spin_interaction_term(Jˣ, σˣ)
    H += spin_spin_interaction_term(Jʸ, σʸ)
    H += spin_spin_interaction_term(Jᶻ, σᶻ)
    H += on_site_term(Γˣ, σˣ)
    H += on_site_term(Γʸ, σʸ)
    H += on_site_term(Γᶻ, σᶻ)
end


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

# must keep the `R` specification to enforce consistent sizes across vector entries
function count_simple_plurality(choices::Vector{RankedChoice{R}})::Vector{Int} where R
  first_counts = countmap(choice[1] for choice in choices)
  count_keys = keys(first_counts)
  candidate_range = 1 : maximum(count_keys) # disregard zero
  [first_counts[index] for index in candidate_range]
end

# todo some actual Condorcet (using all preference info, rather than first) methods. beats(..., ...) would be a cool method name

# knock one out at a time; avoid the complexity of the shortcut to knock out as many from the bottom that add up to less than the next up
function find_instant_runoff_winner(choices::Vector{RankedChoice{R}})::Int where R
  banned = Set{Int}()
  tally = Dict{Int, Int}()
  while length(tally) != 1
    current_votes = Iterators.map(choices) do choice
      for candidate in choice
        if candidate in banned
          continue
        end
        if candidate == 0
          continue
        end
        return candidate
      end
      return 0
    end
    tally = countmap(current_votes)
    delete!(tally, 0)
    min_count, min_candidate = findmin(tally)
    push!(banned, min_candidate)
  end
  tally |> keys |> first
end

function find_condorcet_winner(choices::AbstractVector{RankedChoice{R}},
    allow_ties::Bool = false)::Int where R
  N = maximum(reinterpret(Int, choices))
  differences = matrix_of_differences(N, choices)
  for possibility in 1:N
    winner = true
    for alternative in 1:N
      if alternative == possibility
        continue
      end
      operator = allow_ties ? (<) : (<=)
      if operator(differences[possibility, alternative], 0)
        winner = false
        break
      end
    end
    winner && return possibility # ever care about the chance for more than one
  end
end

# here zero signifies a tie, with indiscernible votes having no influence
matrix_of_differences(N::Int, choices::AbstractVector{RankedChoice{R}}) where R = [ # Val{N} ?
  sum(choices) do choice
    winner = find_pairwise_winner(choice, row, col)
    (winner == row ? +1 : (winner == col ? -1 : 0))
  end
for row in 1:N, col in 1:N ]

function find_pairwise_winner(choice::RankedChoice,
    first_candidate::Int, second_candidate::Int)::Int
  if first_candidate == second_candidate
    return 0
  end
  # for our current purposes, it suffices to find the candidate that appears first in the ranked preferences
  for chosen_candidate in choice
    if chosen_candidate == first_candidate # could use short-circuit shorthand notation
      return first_candidate
    elseif chosen_candidate == second_candidate
      return second_candidate
    end
  end
  return 0 # zero signifies indiscernible with the information available
end


text:
using GraphFlux
using Test

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


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    include("../src/main_model.jl")

    m = getpage(scenario, use_permafrost, use_seaice)
    run(m)

    if scenario == "NDCs"
        while m[:Discontinuity,:occurdis_occurrencedummy] != [0.,0.,0.,0.,0.,0.,1.,1.,1.,1.]
            println(m[:Discontinuity,:occurdis_occurrencedummy])
            run(m)
        end
    else
        while m[:Discontinuity,:occurdis_occurrencedummy] != [0.,0.,0.,0.,0.,0.,0.,0.,0.,0.]
            println(m[:Discontinuity,:occurdis_occurrencedummy])
            run(m)
        end
    end

    #climate component
    temp=m[:ClimateTemperature,:rt_g_globaltemperature]
    temp_compare=readpagedata(m,"test/validationdata/$valdir/rt_g_globaltemperature.csv")
    @test temp ≈ temp_compare rtol=1e-3

    slr=m[:SeaLevelRise,:s_sealevel]
    slr_compare=readpagedata(m,"test/validationdata/$valdir/s_sealevel.csv")
    @test slr ≈ slr_compare rtol=1e-2

    #Socio-Economics
    gdp=m[:GDP,:gdp]
    gdp_compare=readpagedata(m,"test/validationdata/$valdir/gdp.csv")
    @test gdp ≈ gdp_compare rtol=1

    pop=m[:Population,:pop_population]
    pop_compare=readpagedata(m,"test/validationdata/$valdir/pop_population.csv")
    @test pop ≈ pop_compare rtol=0.001

    #Abatement Costs
    abatement=m[:TotalAbatementCosts,:tct_totalcosts]
    abatement_compare=readpagedata(m,"test/validationdata/$valdir/tct_totalcosts.csv")
    @test abatement ≈ abatement_compare rtol=1e-2

    #Adaptation Costs
    adaptation=m[:TotalAdaptationCosts,:act_adaptationcosts_total]
    adaptation_compare=readpagedata(m, "test/validationdata/$valdir/act_adaptationcosts_tot.csv")
    @test adaptation ≈ adaptation_compare rtol=1e-2

    #Damages
    damages=m[:Discontinuity,:rcons_per_cap_DiscRemainConsumption]
    damages_compare=readpagedata(m,"test/validationdata/$valdir/rcons_per_cap_DiscRemainConsumption.csv")
    @test damages ≈ damages_compare rtol=10
    #SLR damages
    slrdamages=m[:SLRDamages,:rcons_per_cap_SLRRemainConsumption]
    slrdamages_compare=readpagedata(m, "test/validationdata/$valdir/rcons_per_cap_SLRRemainConsumption.csv")
    @test slrdamages ≈ slrdamages_compare rtol=0.1
    #Market damages
    mdamages=m[:MarketDamages,:rcons_per_cap_MarketRemainConsumption]
    mdamages_compare=readpagedata(m,"test/validationdata/$valdir/rcons_per_cap_MarketRemainConsumption.csv")
    @test mdamages ≈ mdamages_compare rtol=1
    #NonMarket Damages
    nmdamages=m[:NonMarketDamages,:rcons_per_cap_NonMarketRemainConsumption]
    nmdamages_compare=readpagedata(m,"test/validationdata/$valdir/rcons_per_cap_NonMarketRemainConsumption.csv")
    @test nmdamages ≈ nmdamages_compare rtol=1

    te = m[:EquityWeighting, :te_totaleffect]
    if scenario == "NDCs"
        te_compare = 1.0320923880568126e9
    elseif scenario == "2 degC Target"
        te_compare = 5.0496576104580307e8
    end
    @test te ≈ te_compare rtol=1e4
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
function registerdrivers(f::Function)
    GDAL.allregister()
    try
        f()
    finally
        GDAL.destroydrivermanager()
    end
end

"Fetch driver by index"
getdriver(i::Integer) = GDAL.getdriver(i-1)

"Fetch a driver based on the short name (such as `GTiff`)."
getdriver(name::AbstractString) = GDAL.getdriverbyname(name)

"""
Destroy a `GDALDriver`.

This is roughly equivelent to deleting the driver, but is guaranteed to take
place in the GDAL heap. It is important this that function not be called on a
driver that is registered with the `GDALDriverManager`.
"""
destroy(drv::Driver) = GDAL.destroydriver(drv)

"Register a driver for use."
register(drv::Driver) = GDAL.registerdriver(drv)

"Deregister the passed drv."
deregister(drv::Driver) = GDAL.deregisterdriver(drv)

"Return the list of creation options of the driver [an XML string]"
options(drv::Driver) = GDAL.getdrivercreationoptionlist(drv)

"Return the short name of a driver (e.g. `GTiff`)"
getshortname(drv::Driver) = GDAL.getdrivershortname(drv)

"Return the long name of a driver (e.g. `GeoTIFF`), or empty string."
getlongname(drv::Driver) = GDAL.getdriverlongname(drv)

"Fetch the number of registered drivers."
ndriver() = GDAL.getdrivercount()

"Returns a listing of all registered drivers"
function drivers()
    dlist = Dict{ASCIIString,ASCIIString}()
    for i in 1:ndriver()
        dlist[shortname(getdriver(i))] = longname(getdriver(i))
    end
    dlist
end

"""
Identify the driver that can open a raster file.

This function will try to identify the driver that can open the passed filename
by invoking the Identify method of each registered `GDALDriver` in turn. The
first driver that successful identifies the file name will be returned. If all
drivers fail then `NULL` is returned.
"""
identifydriver(filename::AbstractString) = GDAL.identifydriver(filename, C_NULL)

"""
Validate the list of creation options that are handled by a drv.

This is a helper method primarily used by `Create()` and `CreateCopy()` to
validate that the passed in list of creation options is compatible with the
`GDAL_DMD_CREATIONOPTIONLIST` metadata item defined by some drivers.

See also: `GDALGetDriverCreationOptionList()`

If the `GDAL_DMD_CREATIONOPTIONLIST` metadata item is not defined, this
function will return `TRUE`. Otherwise it will check that the keys and values
in the list of creation options are compatible with the capabilities declared
by the `GDAL_DMD_CREATIONOPTIONLIST` metadata item. In case of incompatibility
a (non fatal) warning will be emited and `FALSE` will be returned.

### Parameters
* `hDriver`     the handle of the driver with whom the lists of creation option
must be validated
* `options`     the list of creation options. An array of strings, whose last
element is a `NULL` pointer

### Returns
`TRUE` if the list of creation options is compatible with the `Create()` and
`CreateCopy()` method of the driver, `FALSE` otherwise.
"""
function validate{T <: AbstractString}(drv::Driver, options::Vector{T})
    Bool(GDAL.validatecreationoptions(drv,
                                      Ptr{Ptr{UInt8}}(pointer(options))))
end


text:
###################################################################################################
#                                        Load Packages
###################################################################################################
cd(@__DIR__)
using Pkg
Pkg.activate("../..")
using Revise, ACTRSimulators
import_gui()
import ACTRSimulators: start!, press_key!, repaint!
include("PVT.jl")
include("PVT_Model.jl")
###################################################################################################
#                                        Run Model
###################################################################################################
scheduler = ACTRScheduler(;model_trace=true)
task = PVT(;scheduler, n_trials=2, visible=true, realtime=true)
model = init_model(scheduler, task)
run!(model, task)


text:
using Base.Iterators
using ExprParsers
using Test

@testset "parse_expr plain values" begin
  @test parse_expr(5, 5) == 5
  @test_throws EP.ParseError parse_expr(6, 7)
  @test_throws EP.ParseError parse_expr("hi", nothing)
end

@testset "parse_expr EP.anything" begin
  @test parse_expr(EP.anything, 5) == 5
  @test parse_expr(EP.anything, "hi") == "hi"
  @test parse_expr(EP.anything, nothing) == nothing
end

@testset "parse_expr vectors" begin
  @test parse_expr([EP.anything, EP.anything], [1,2]) == [1,2]
  @test parse_expr([EP.anything, 2], [1,2]) == [1,2]
  @test_throws EP.ParseError parse_expr([EP.anything, 3], [1,2]) == [1,2]
end

@testset "AnyOf,AllOf" begin
  @test_throws EP.ParseError parse_expr(EP.AnyOf("hi", "ho", "hu"), "aa")
  @test parse_expr(EP.AnyOf("hi", "ho", "hu"), "hu") == "hu"
  @test parse_expr(EP.AllOf("a", EP.anything), "a") == "a"
  @test_throws EP.ParseError parse_expr(EP.AllOf("a", EP.anything, "b"), "a")
end

@testset "Isa" begin
  matcher = EP.Isa(Number)
  @test parse_expr(matcher, 4) == 4
  @test parse_expr(matcher, 4.412) == 4.412
  @test parse_expr(matcher, 3//5) == 3//5
  @test_throws EP.ParseError parse_expr(matcher, "hi")
end

@testset "SatisfiesPredicate" begin
  matcher = EP.SatisfiesPredicate() do x
    x in (1,3,6,7)
  end
  @test_throws EP.ParseError parse_expr(matcher, 2)
  @test parse_expr(matcher, 6) == 6
end

@testset "Iterator" begin
  @test parse_expr(EP.anything, [1,2,3,4]) == parse_expr(EP.Utils.Iterator(repeated(EP.anything)), [1,2,3,4]) == [1,2,3,4]
end


text:
using SupervisedLearning
using RDatasets
using UnicodePlots

data = dataset("datasets", "mtcars")
formula = AM ~ DRat + WT

#-----------------------------------------------------------

msg("Encoding dataframe per hand")

mf = ModelFrame(formula, data)
mm = ModelMatrix(mf)
dfBias = in(0,mm.assign)
X = dfBias ? mm.m[:,2:end]: mm.m
extBias = dfBias * 1.0
t = convert(Vector{Float64}, model_response(mf))

ce = ZeroOneClassEncoding(["No", "Yes"])

ds = dataSource(X', t, ce, bias=extBias)
@test typeof(ds) <: EncodedInMemoryLabeledDataSource

#-----------------------------------------------------------

msg("Encoding dataframe with dataSource")

ds = dataSource(formula, data, SignedClassEncoding)
@test typeof(ds) <: EncodedInMemoryLabeledDataSource

#-----------------------------------------------------------

msg("Encoding dataframe with encodeDataSource")

ds = dataSource(formula, data)
@test typeof(ds) <: DataFrameLabeledDataSource
ds = encodeDataSource(ds, ZeroOneClassEncoding)
@test typeof(ds) <: EncodedInMemoryLabeledDataSource

#barplot(classDistribution(ds)...)


text:
using BenchmarkTools

function solve_day2()
    x, d1, d2, aim = 0, 0, 0, 0
    @inbounds for line in eachline("../data/input2.txt")
        command, vals = split(line)
        val = parse(Int, vals)
        if startswith(command, 'f') # forward
            x += val
            d2 += aim * val
        elseif startswith(command, 'd') # down
            aim += val
            d1 += val
        else # up
            aim -= val
            d1 -= val
        end
    end
    return x * d1, x * d2
end

@show solve_day2()
@btime solve_day2()

max_stars_repo_name: lucaferranti/AoC

max_issues_repo_name: lucaferranti/AoC

max_forks_repo_name: lucaferranti/AoC


text:
using ShowItLikeYouBuildIt
using Base.Test
push!(LOAD_PATH, @__DIR__)
using ZArrays
pop!(LOAD_PATH)

const showit_old = VERSION < v"0.7.0-DEV.1790"
const eltype_string = showit_old ? "element type" : "eltype"

if showit_old
    printdims(io, dims) = print(io, dims)
else
    printdims(io, dims) = _printdims(io, dims...)
    _printdims(io, d1, d...) = (print(io, d1, ", "); _printdims(io, d...))
    _printdims(io, d1) = print(io, d1)
end
function printdims(dims)
    io = IOBuffer()
    printdims(io, dims)
    String(take!(io))
end

if showit_old
    for T in (Float64, Bool, Int8, UInt16, Symbol, String, Union{})
        @test type_complexity(T) == 1
    end
    @test type_complexity(Union{Float32,Int8}) == 2
    @test type_complexity(Array{Int,2}) == 3
    @test type_complexity(Array{Complex{Float32},2}) == 4
    @test type_complexity(Array{Array{Complex{Float32},1},2}) == 6
end

@test ShowItLikeYouBuildIt.dimstring(()) == "0-dimensional"

@test sprint(showarg, 3.0) == "::Float64"
@test sprint(showarg, Float64) == "::Type{Float64}"

# Display of objects

# Note that these next tests alter how SubArrays, ReshapedArrays, and
# PermutedDimsArrays are displayed; it's probably best to run this
# only in a "throwaway" julia session.

function ShowItLikeYouBuildIt.showarg(io::IO, v::SubArray)
    print(io, "view(")
    showarg(io, parent(v))
    print(io, ", ", join(v.indexes, ", "))
    print(io, ')')
end

function ShowItLikeYouBuildIt.showarg(io::IO, A::PermutedDimsArray{T,N,perm}) where {T,N,perm}
    print(io, "PermutedDimsArray(")
    showarg(io, parent(A))
    print(io, ", ", perm, ')')
end

function ShowItLikeYouBuildIt.showarg(io::IO, A::Base.ReshapedArray)
    print(io, "reshape(")
    showarg(io, parent(A))
    print(io, ", ")
    printdims(io, A.dims)
    print(io, ')')
end

Base.summary(A::SubArray) = summary_build(A)
Base.summary(A::Base.PermutedDimsArrays.PermutedDimsArray) = summary_build(A)
Base.summary(A::Base.ReshapedArray) = summary_build(A)

a = rand(3,5,7)
v = view(a, :, 3, 2:5)
view_name = showit_old ? "view(::Array{Float64,3}, $(v.indexes[1]), 3, 2:5)" :
                         "view(::Array{Float64,3}, :, 3, 2:5)"
@test summary(v) == "3×4 $view_name with $eltype_string Float64"

c = reshape(v, 4, 3)
str = summary(c)
dimstr = printdims((4, 3))
@test str == "4×3 reshape($view_name, $dimstr) with $eltype_string Float64"

a = reshape(1:24, 3, 4, 2)
b = PermutedDimsArray(a, (2,3,1))
str = summary(b)
intstr = string("Int", Sys.WORD_SIZE)
dimstr = printdims((3, 4, 2))
@test str == "4×2×3 PermutedDimsArray(reshape(::UnitRange{$intstr}, $dimstr), $((2, 3, 1))) with $eltype_string $intstr"

o = ZArray(rand(3,5))
vo = view(o, 0:2:2, :)
if showit_old
    @test summary(vo) == "Base.OneTo(2)×ZArrays.ZeroRange(5) view(::ZArrays.ZArray{Float64,2}, 0:2:2, $(vo.indexes[end])) with $eltype_string Float64"
else
    @test summary(vo) == "view(::ZArrays.ZArray{Float64,2}, 0:2:2, :) with eltype Float64 with indices Base.OneTo(2)×ZArrays.ZeroRange(5)"
end

if showit_old
    # Adding cthresh is deprecated
    Base.summary(A::SubArray) = summary_build(A,1000)
    @test summary(v) == "3×4 SubArray{Float64,2,Array{Float64,3},Tuple{$(typeof(v.indexes[1])),$Int,UnitRange{$Int}},false}"
end

nothing


text:
# simple program to test the new k-means (not ready yet)

using Clustering

m = 3
n = 1000
k = 10

x = rand(m, n)

println("non-weighted")
r = kmeans(x, k; max_iter=50, display=:iter)
println()

println("weighted")
w = rand(n)
r = kmeans(x, k; max_iter=50, display=:iter, weights=w)
println()


text:
using Patchwork
using Reactive

output_path(file, dir) = begin
    out = abspath(joinpath(dir, file))
    return replace(out, r".jl$", ".html")
end

getvalue(ui::Signal) = value(ui)
getvalue(x) = x

mkoutputdir!(dir) = begin
    if isdir(dir)
        return
    elseif isfile(dir)
        error("$dir is not a directory")
    else
        mkdir(dir)
    end
end

escher_make(file, output_dir; single_file=false, assets_dir="pkg/Escher", copy_assets=false) = begin

    opath = output_path(file, output_dir)
    w = Window()
    Reactive.stop_event_loop()
    assets = foldp(push!, Any[], w.assets) # Accumulate assets
    uifn = include(joinpath(pwd(), file))
    ui = uifn(w)
    Reactive.run_till_now()

    mkoutputdir!(output_dir)

    asset_src  = Pkg.dir("Escher", "assets")
    asset_dest = joinpath(output_dir, assets_dir) |> abspath

    dir = dirname(asset_dest)
    if !isdir(dir)
        mkpath(dir)
    end
    if copy_assets
        cp(asset_src, asset_dest)
    else
        if !isfile(asset_dest) && !isdir(asset_dest)
            symlink(asset_src, asset_dest)
        end
    end

    state = Dict()
    state["embedded_signals"] = Dict()

    open(opath, "w") do io
        write(io, """<!doctype html>
        <html>
        <meta charset="utf-8">
        <head>
           <script> $(Patchwork.js_runtime()) </script>
           <script src="$assets_dir/bower_components/webcomponentsjs/webcomponents.min.js"></script>
        </head>

        $(

        join(map(x -> """<link rel="import" href="$(Escher.resolve_asset(x))">""",
                      vcat("basics", value(assets))), "\n")

        )

        <body>
        <signal-container id="root"></signal-container>
        </div>
        <script>
            new Patchwork.Node("root", $(replace(replace(Patchwork.jsonfmt(Escher.render(getvalue(ui), state)) |> json, "<", "\\<"), ">", "\\>")))
        </script>
        </body>
        </html>
        """)
    end
end


text:
using Test
using SparseArrays
using LinearAlgebra
using RandomProjectionTrees: build_rptree, rp_split, num_leaves, num_nodes, leaves

@testset "RandomProjectionTrees tests" begin

end


text:
using Turing
import Turing: consume, produce, TArray, get, tzeros
using Test

ta1 = TArray{Int}(4);
push!(ta1, 1);
push!(ta1, 2);
@test pop!(ta1) == 2

ta1_2 = TArray{Int, 1}(4); # another constructor
push!(ta1_2, 1);
push!(ta1_2, 2);
@test pop!(ta1_2) == 2

ta2 = TArray{Int}(4, 4);
ta3 = TArray{Int, 4}(4, 3, 2, 1);
ta4 = get(ta3);

@test ta3[3] == ta4[3]

ta5 = TArray{Int}(4);
for i in 1:4 ta5[i] = i end
@test Array(ta5) == [1, 2, 3, 4]

@test Array(tzeros(4)) == zeros(4)

ta6 = TArray{Float64}(4);
for i in 1:4 ta6[i] = i / 10 end
@test ta6[1] == 0.1
@test Array(ta6) == [0.1, 0.2, 0.3, 0.4]

ta7 = TArray{Int, 2}((2, 2));   # TODO: add test for use this multi-dim array


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(Scene, ())
    @assert precompile(update_limits!, (Scene,))
    @assert precompile(update_limits!, (Scene, Automatic))
    @assert precompile(update_limits!, (Scene, FRect3D))

    @assert precompile(boundingbox, (Scene,))
    Ngonf0 = GeometryBasics.Ngon{2, Float32, 3, Point2f0}
    # @assert precompile(boundingbox, (Mesh{Tuple{Vector{GeometryBasics.Mesh{2, Float32, Ngonf0, FaceView{Ngonf0}}}}},))  # doesn't seem to work
    @assert precompile(boundingbox, (Poly{Tuple{Vector{Vector{Point2f0}}}},))
    @assert precompile(boundingbox, (String, Vector{Point3f0}, Vector{Float32}, Vector{FTFont}, Vec2f0, Vector{Quaternionf0}, SMatrix{4, 4, Float32, 16}, Float64, Float64))

    @assert precompile(poly_convert, (Vector{Vector{Point2f0}},))
    @assert precompile(rotatedrect, (FRect2D, Float32))

    @assert precompile(plot!, (Scene, Type{Poly{Tuple{IRect2D}}}, Attributes, Tuple{Observable{IRect2D}}, Observable{Tuple{Vector{Vector{Point2f0}}}}))
    @assert precompile(plot!, (Mesh{Tuple{Vector{GeometryBasics.Mesh{2, Float32, Ngonf0, FaceView{Ngonf0}}}}},))
    @assert precompile(plot!, (Scene, Type{Annotations{Tuple{Vector{Tuple{String, Point2f0}}}}}, Attributes, Tuple{Observable{Vector{Tuple{String, Point2f0}}}}, Observable{Tuple{Vector{Tuple{String, Point2f0}}}}))
end


text:
rules = Dict()
open("day07/day07.in") do io
    global rules
    for line in readlines(io)
        line = replace(line, ("." => ""))
        line = split(line, "contain ")
        key_bag = split(line[1], " ")[1] * " " * split(line[1], " ")[2]
        bags = split(line[2], ", ")
        rules[key_bag] = []
        for bag in bags
            push!(rules[key_bag], split(bag, " ")[2] * " " * split(bag, " ")[3])
        end
    end
end

contains_gold = 0
for (bag_type, bags) in rules
    global contains_gold
    used_bags = Set(); push!(used_bags, bag_type)
    heap_bags = Set(); push!(heap_bags, bag_type)

    while length(heap_bags) > 0
        for bag_type_current in heap_bags
            push!(used_bags, bag_type_current)
            delete!(heap_bags, bag_type_current)
            if get(rules, bag_type_current, 0) != 0
                bags = rules[bag_type_current]
                for bag in bags
                    if !in(bag, used_bags)
                        push!(heap_bags, bag)
                    end
                end
            end
        end
    end
    if bag_type != "shiny gold" && in("shiny gold", used_bags)
        contains_gold += 1
    end
end

println(contains_gold)

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
using Compat: mapreduce

function _treehistory(tree::T, nodename::NL) where
    {RT <: Rooted, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    branches, nodes = _treehistory(tree, getnode(tree, nodename))
    return [getbranchname(tree, b) for b in branches],
           [getnodename(tree, n) for n in nodes]
end

function _treehistory(tree::T, node::N) where
    {RT <: Rooted, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    branches = B[]
    nodestoprocess = N[node]
    nodesprocessed = N[]
    while !isempty(nodestoprocess)
        nextnode = pop!(nodestoprocess)
        push!(nodesprocessed, nextnode)
        if hasinbound(tree, nextnode)
            push!(branches, getinbound(tree, nextnode))
            push!(nodestoprocess, getparent(tree, nextnode))
        end
    end
    return branches, nodesprocessed
end

function _treefuture(tree::T, nodename::NL) where
    {RT <: Rooted, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    branches, nodes = _treefuture(tree, getnode(tree, nodename))
    return [getbranchname(tree, b) for b in branches],
           [getnodename(tree, n) for n in nodes]
end

function _treefuture(tree::T, node::N) where
    {RT <: Rooted, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    branches = B[]
    nodestoprocess = N[node]
    nodesprocessed = N[]
    while !isempty(nodestoprocess)
        nextnode = pop!(nodestoprocess)
        push!(nodesprocessed, nextnode)
        append!(branches, getoutbounds(tree, nextnode))
        append!(nodestoprocess, getchildren(tree, nextnode))
    end
    return branches, nodesprocessed
end

"""
    branchhistory(tree::AbstractTree, node)

Find the branch route between a node on a tree and its root
"""
function branchhistory(tree::T, node::Union{NL, N}) where
    {RT <: Rooted, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    return _treehistory(tree, node)[1]
end

"""
    nodehistory(tree::AbstractTree, node)

Find the node route between a node on a tree and its root
"""
function nodehistory(tree::T, node::Union{NL, N}) where
    {RT <: Rooted, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    return _treehistory(tree, node)[2]
end

"""
    branchfuture(tree::AbstractTree, node)

Find the branches between a node on a tree and its leaves
"""
function branchfuture(tree::T, node::Union{NL, N}) where
    {RT <: Rooted, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    return _treefuture(tree, node)[1]
end

"""
    nodefuture(tree::AbstractTree, node)

Find the nodes between a node on a tree and its leaves
"""
function nodefuture(tree::T, node::Union{NL, N}) where
    {RT <: Rooted, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    return _treefuture(tree, node)[2]
end

"""
    branchfuture(tree::AbstractTree, node)

Find the branches between a node on a tree and its leaves
"""
function branchfuture(tree::T, node::NL) where {NL, BL, T <: AbstractTree{NL, BL}}
    return _treefuture(tree, node)[1]
end

"""
    nodefuture(tree::AbstractTree, node)

Find the nodes between a node on a tree and its leaves
"""
function nodefuture(tree::T, node::NL) where {NL, BL, T <: AbstractTree{NL, BL}}
    return _treefuture(tree, node)[2]
end

"""
    branchroute(tree::AbstractTree, node1, node2)

Find the branch route between two nodes on a tree
"""
function branchroute(tree::T, node1::NL, node2::NL) where
    {RT <: Rooted, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    route = branchroute(tree, getnode(tree, node1), getnode(tree, node2))
    return [getbranchname(tree, b) for b in route]
end
function branchroute(tree::T, node1::N, node2::N) where
    {RT <: Rooted, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    branches1, nodes1 = _treehistory(tree, node1)
    branches2, nodes2 = _treehistory(tree, node2)
    nodes1[end] == nodes2[end] ||
        return error("No route between nodes")
    common = branches1 ∩ branches2
    nodes1[end] == nodes2[end] || error("No route between nodes in tree")
    return append!([b for b in branches1 if b ∉ common],
                   [b for b in reverse(branches2) if b ∉ common])
end

"""
    noderoute(tree::AbstractTree, node1, node2)

Find the node route between two nodes on a tree
"""
function noderoute(tree::T, node1::NL, node2::NL) where
    {RT <: Rooted, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    route = noderoute(tree, getnode(tree, node1), getnode(tree, node2))
    return [getnodename(tree, n) for n in route]
end
function noderoute(tree::T, node1::N, node2::N) where
    {RT <: Rooted, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    branches1, nodes1 = _treehistory(tree, node1)
    branches2, nodes2 = _treehistory(tree, node2)
    nodes1[end] == nodes2[end] || error("No route between nodes in tree")
    common = nodes1[end]
    while min(length(nodes1), length(nodes2)) > 0 && nodes1[end] == nodes2[end]
        common = nodes1[end]
        pop!(nodes1)
        pop!(nodes2)
    end
    push!(nodes1, common)
    return append!(nodes1, reverse(nodes2))
end

"""
    distance(tree::AbstractTree, node1, node2)

Distance between two nodes on a tree
"""
function distance(tree::AbstractTree, node1, node2)
    branches = branchroute(tree, node1, node2)
    return mapreduce(branch -> getlength(tree, branch), +, branches;
    init = 0.0)
end

"""
    distances(tree::AbstractTree)

Pairwise distances between all leaf nodes on a tree
"""
function distances(tree::AbstractTree)
    leaves = getleaves(tree)
    return [distance(tree, li, lj) for li in leaves, lj in leaves]
end

"""
    height(tree::AbstractTree, node)

Height of a node of the tree above the root
"""
function heighttoroot(tree::AbstractTree, node)
    return mapreduce(branch -> getlength(tree, branch), +,
                     branchhistory(tree, node); init = 0.0)
end

"""
    heights(tree::AbstractTree)

Height of all of the leaves of the tree above the root
"""
function heightstoroot(tree::AbstractTree)
    return [heighttoroot(tree, leaf) for leaf in getleaves(tree)]
end


text:
struct Firm
  F::Function #production as a function of capital
  τ::Vector{Float64} #tax
  δ::Float64 #depreciation

  Θ::Vector{Float64} #technology per state
  π::Matrix{CartesianIndex} #policy states
  v::Matrix{Float64} #value states

  R::Array{Float64,3}
end

#gets the cobb-douglas produciton function
function getCD(γ::Float64)::Function
  F(k::Float64)::Float64 = k^γ

  return F
end

function Firm(τ::Vector{Float64} = DEF_τ,
  δ::Float64 = DEF_δ,
  Θ::Vector{Float64} = DEF_Θ,
  γ::Float64 = DEF_γ,
  kstates::Vector{Float64} = DEF_KSTATES,
  ystates::Vector{Float64} = DEF_YSTATES)

  kystates = combine(kstates, ystates)
  N::Int = length(kystates)
  M::Int = length(Θ)
  F = getCD(γ)

  π::Matrix{CartesianIndex} = Matrix{CartesianIndex}(undef, N, M)
  v::Matrix{Float64} = Matrix{Float64}(undef, N, M)
  R::Array{Float64,3} = Array{Float64,3}(undef, N, N, M)

  return Firm(F, τ, δ, Θ, π, v, R)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using ModelingStructs
using DifferentialEquations

# Lorenz system
function lorenz!(du, u, p, t)
    du.x = p.σ*(u.y - u.x)
    du.y = u.x*(p.ρ - u.z) - u.y - p.f
    du.z = u.x*u.y - p.β*u.z
    return nothing
end

lorenz_p = (σ=10.0, ρ=28.0)
lorenz_ic = (x=0.0, y=0.0, z=0.0)

# Lotka-Volterra system
function lotka!(du, u, p, t)
    du.x =  p.α*u.x - p.β*u.x*u.y + p.f
    du.y = -p.γ*u.y + p.δ*u.x*u.y
    return nothing
end

lotka_p = (α=1.0, γ=1.1, δ=0.5)
lotka_ic = (x=1.0, y=1.0)

# Composed Lorenz and Lotka-Volterra system
function composed!(du, u, p, t)
    lorenz!(du.lorenz, u.lorenz, (β=p.β, f=u.lotka.x, p.lorenz...), t)
    lotka!(  du.lotka,  u.lotka, (β=p.β, f=u.lorenz.x, p.lotka...), t)
    return nothing
end

comp_p = (β=8/3, lorenz=lorenz_p, lotka=lotka_p)
comp_ic = MStruct(lorenz=lorenz_ic, lotka=lotka_ic)

# Create and solve problem
prob = ODEProblem(composed!, comp_ic, (0.0, 20.0), comp_p)
sol = solve(prob, Tsit5())

max_stars_repo_name: jonniedie/NamedViewVectors


text:
@testset "mst_prim" begin
    A = load_matrix_network("airports")
    A = -A
    T = mst_prim_matrix(A)

    A = load_matrix_network("clr-24-1")
    A[2,3] = 9
    A[3,2] = 9
    T = mst_prim_matrix(A)
    Ttrueijv = [
     2     8     1     4     6     9     3     5     4     3     7     6     8     1     7     3
     1     1     2     3     3     3     4     4     5     6     6     7     7     8     8     9
     4     8     4     7     4     2     7     9     9     4     2     2     1     8     1     2 ]
    Ttrue = sparse(vec(Ttrueijv[1,:]),  vec(Ttrueijv[2,:]), vec(Ttrueijv[3,:]), 9,9)

    @test nnz(T - Ttrue) == 0

    A = load_matrix_network("clr-24-1")
    T1 = mst_prim_matrix(A)
    T2 = mst_prim_matrix(A,false,5)
    T3 = mst_prim_matrix(MatrixNetwork(A),false)
    T1T2diff = sparse([2,1],[3,8],[8,-8],9,9)
    T3T2diff = sparse([2,1],[3,8],[8,-8],9,9)

    @test nnz(triu(T1-T2) - T1T2diff) == 0
    @test nnz(triu(T3-T2) - T3T2diff) == 0

    @test_throws ErrorException mst_prim(MatrixNetwork(-sparse(1.0I,3,3)))
    @test_throws ErrorException mst_prim(MatrixNetwork(-sparse(1.0I,3,3)), false)
    @test_throws ErrorException mst_prim(MatrixNetwork(-sparse(1.0I,3,3)), false, 1)
    @test_throws ErrorException mst_prim(-sparse(1.0I,3,3))
    @test_throws ErrorException mst_prim(-sparse(1.0I,3,3), false)
    @test_throws ErrorException mst_prim(-sparse(1.0I,3,3), false, 1)
end


text:
module Cones

using LinearAlgebra: I, rank, dot
using MathProgBase, Gurobi
using Printf: @printf

import Base: length

export chernikova, chernikova_general, Cone, FiniteCone, PolyhedralCone, project, Polytope, find_generators, min_max_projections

# package code goes here
include("utils.jl")
include("LCP_julia.jl")
include("chernikova.jl")
include("chernikova_general.jl")
include("types.jl")
include("polytope.jl")

end # module


text:
"""
    Dropout(p; dims=:)

Dropout layer.

## Arguments

* `p`: Probability of Dropout (if `p = 0` then [`NoOpLayer`](@ref) is returned)

## Keyword Arguments

* To apply dropout along certain dimension(s), specify the `dims` keyword. e.g. `Dropout(p; dims = 3)` will randomly zero out entire channels on WHCN input (also called 2D dropout).

## Inputs

* `x`: Must be an AbstractArray

## Returns

* `x` with dropout mask applied if `training=Val(true)` else just `x`
* State with updated `rng`

## States

* `rng`: Pseudo Random Number Generator
* `training`: Used to check if training/inference mode

Call [`Lux.testmode`](@ref) to switch to test mode.

See also [`VariationalHiddenDropout`](@ref)
"""
struct Dropout{T, D} <: AbstractExplicitLayer
    p::T
    dims::D
end

function initialstates(rng::AbstractRNG, ::Dropout)
    # FIXME: Take PRNGs seriously
    randn(rng, 1)
    return (rng=replicate(rng), training=Val(true))
end

function Dropout(p; dims=:)
    @assert 0 ≤ p ≤ 1
    iszero(p) && return NoOpLayer()
    return Dropout(p, dims)
end

function (d::Dropout{T})(x::AbstractArray{T}, ps, st::NamedTuple) where {T}
    y, _, rng = dropout(st.rng, x, d.p, d.dims, st.training)
    return y, merge(st, (rng=rng,))
end

function Base.show(io::IO, d::Dropout)
    print(io, "Dropout(", d.p)
    d.dims != Colon() && print(io, ", dims=", d.dims)
    return print(io, ")")
end

"""
    VariationalHiddenDropout(p; dims=:)

VariationalHiddenDropout layer. The only difference from Dropout is that the `mask` is retained until [`Lux.update_state(l, :update_mask, Val(true))`](@ref) is called.

## Arguments

* `p`: Probability of Dropout (if `p = 0` then [`NoOpLayer`](@ref) is returned)

## Keyword Arguments

* To apply dropout along certain dimension(s), specify the `dims` keyword. e.g. `VariationalHiddenDropout(p; dims = 3)` will randomly zero out entire channels on WHCN input (also called 2D dropout).

## Inputs

* `x`: Must be an AbstractArray

## Returns

* `x` with dropout mask applied if `training=Val(true)` else just `x`
* State with updated `rng`

## States

* `rng`: Pseudo Random Number Generator
* `training`: Used to check if training/inference mode
* `mask`: Dropout mask. Initilly set to nothing. After every run, contains the mask applied in that call
* `update_mask`: Stores whether new mask needs to be generated in the current call

Call [`Lux.testmode`](@ref) to switch to test mode.

See also [`Dropout`](@ref)
"""
struct VariationalHiddenDropout{T, D} <: AbstractExplicitLayer
    p::T
    dims::D
end

function initialstates(rng::AbstractRNG, ::VariationalHiddenDropout)
    # FIXME: Take PRNGs seriously
    randn(rng, 1)
    return (rng=replicate(rng), training=Val(true), update_mask=Val(true),
            mask=nothing)
end

function VariationalHiddenDropout(p; dims=:)
    @assert 0 ≤ p ≤ 1
    iszero(p) && return NoOpLayer()
    return VariationalHiddenDropout(p, dims)
end

function (d::VariationalHiddenDropout{T})(x::AbstractArray{T}, ps, st::NamedTuple) where {T}
    y, mask, rng, update_mask = dropout(st.rng, x, st.mask, d.p, d.dims, st.training,
                                        st.update_mask)
    return y, merge(st, (mask=mask, rng=rng, update_mask=update_mask))
end

function Base.show(io::IO, d::VariationalHiddenDropout)
    print(io, "VariationalHiddenDropout(", d.p)
    d.dims != Colon() && print(io, ", dims=", d.dims)
    return print(io, ")")
end


text:
# This file contans the common interface for LightGraphs.
# TODO 0.7: reevaluate use of errors here.

_NI(m) = error("Not implemented: $m")

"""
    AbstractEdge

An abstract type representing a single edge between two vertices of a graph.
"""
abstract type AbstractEdge{T} end

"""
    AbstractEdgeIter

An abstract type representing an edge iterator.
"""
abstract type AbstractEdgeIter end

"""
    AbstractGraph

An abstract type representing a graph.
"""
abstract type AbstractGraph{T} end

@traitdef IsDirected{G<:AbstractGraph}
@traitimpl IsDirected{G} <- is_directed(G)

#
# Interface for AbstractEdges
#

"""
    src(e)

Return the source vertex of edge `e`.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(2);

julia> add_edge!(g, 1, 2);

julia> src(first(edges(g)))
1
```
"""
src(e::AbstractEdge) = _NI("src")

"""
    dst(e)

Return the destination vertex of edge `e`.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(2);

julia> add_edge!(g, 1, 2);

julia> dst(first(edges(g)))
2
```
"""
dst(e::AbstractEdge) = _NI("dst")

Pair(e::AbstractEdge) = _NI("Pair")
Tuple(e::AbstractEdge) = _NI("Tuple")

"""
    reverse(e)

Create a new edge from `e` with source and destination vertices reversed.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleDiGraph(2);

julia> add_edge!(g, 1, 2);

julia> reverse(first(edges(g)))
Edge 2 => 1
```
"""
reverse(e::AbstractEdge) = _NI("reverse")

==(e1::AbstractEdge, e2::AbstractEdge) = _NI("==")

#
# Interface for AbstractGraphs
#
"""
    edgetype(g)

Return the type of graph `g`'s edge
"""
edgetype(g::AbstractGraph) = _NI("edgetype")

"""
    eltype(g)

Return the type of the graph's vertices (must be <: Integer)
"""
eltype(g::AbstractGraph) = _NI("eltype")

"""
    nv(g)

Return the number of vertices in `g`.

# Examples
```jldoctest
julia> using LightGraphs

julia> nv(SimpleGraph(3))
3
```
"""
nv(g::AbstractGraph) = _NI("nv")

"""
    ne(g)

Return the number of edges in `g`.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = PathGraph(3);

julia> ne(g)
2
```
"""
ne(g::AbstractGraph) = _NI("ne")

"""
    vertices(g)

Return (an iterator to or collection of) the vertices of a graph.

### Implementation Notes
A returned iterator is valid for one pass over the edges, and
is invalidated by changes to `g`.

# Examples
```jldoctest
julia> using LightGraphs

julia> collect(vertices(SimpleGraph(4)))
4-element Array{Int64,1}:
 1
 2
 3
 4
```
"""
vertices(g::AbstractGraph) = _NI("vertices")

"""
    edges(g)

Return (an iterator to or collection of) the edges of a graph.
For `AbstractSimpleGraph`s it returns a `SimpleEdgeIter`.
The expressions `e in edges(g)` and `e ∈ edges(ga)` evaluate as
calls to [`has_edge`](@ref).

### Implementation Notes
A returned iterator is valid for one pass over the edges, and
is invalidated by changes to `g`.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = PathGraph(3);

julia> collect(edges(g))
2-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 2 => 3
```
"""
edges(g) = _NI("edges")

"""
    is_directed(g)

Return `true` if the graph is a directed graph; `false` otherwise.

# Examples
```jldoctest
julia> using LightGraphs

julia> is_directed(SimpleGraph(2))
false

julia> is_directed(SimpleDiGraph(2))
true
```
"""
is_directed(g) = _NI("is_directed")
is_directed(::Type{T}) where T = _NI("is_directed")

"""
    has_vertex(g, v)

Return true if `v` is a vertex of `g`.

# Examples
```jldoctest
julia> using LightGraphs

julia> has_vertex(SimpleGraph(2), 1)
true

julia> has_vertex(SimpleGraph(2), 3)
false
```
"""
has_vertex(x, v) = _NI("has_vertex")

"""
    has_edge(g, s, d)

Return true if the graph `g` has an edge from node `s` to node `d`.

An optional `has_edge(g, e)` can be implemented to check if an edge belongs
to a graph, including any data other than source and destination node.

`e ∈ edges(g)` or `e ∈ edges(g)` evaluate as
calls to `has_edge`, c.f. [`edges`](@ref).

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleDiGraph(2);

julia> add_edge!(g, 1, 2);

julia> has_edge(g, 1, 2)
true

julia> has_edge(g, 2, 1)
false
```
"""
has_edge(g, s, d) = _NI("has_edge")
has_edge(g, e) = has_edge(g, src(e), dst(e))

"""
    inneighbors(g, v)

Return a list of all neighbors connected to vertex `v` by an incoming edge.

### Implementation Notes
Returns a reference, not a copy. Do not modify result.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0 0 0; 0 0 1 0 0; 1 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0]);

julia> inneighbors(g, 4)
2-element Array{Int64,1}:
 3
 5
```
"""
inneighbors(x, v) = _NI("inneighbors")

"""
    outneighbors(g, v)

Return a list of all neighbors connected to vertex `v` by an outgoing edge.

# Implementation Notes
Returns a reference, not a copy. Do not modify result.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0 0 0; 0 0 1 0 0; 1 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0]);

julia> outneighbors(g, 4)
1-element Array{Int64,1}:
 5
```
"""
outneighbors(x, v) = _NI("outneighbors")

"""
    zero(g)

Return a zero-vertex, zero-edge version of the same type of graph as `g`.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0 0 0; 0 0 1 0 0; 1 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0]);

julia> zero(g)
{0, 0} directed simple Int64 graph
```
"""
zero(g::AbstractGraph) = _NI("zero")


text:
module AbstractLayers

export AbstractLayer, AL, forward, backward, parms, ∇parms, idims, odims

include("./tensor_dt.jl")
using .TensorDT: Tensor

## ======================================================================
## The Layer Abstraction
## ======================================================================

abstract type AbstractLayer end
const AL = AbstractLayer

forward(::AbstractLayer, ::Tensor) = throw(ArgumentError("Not Implemented"))
backward(::AbstractLayer, ::Tensor) = throw(ArgumentError("Not Implemented"))
parms(::AbstractLayer) = throw(ArgumentError("Not Implemented"))
∇parms(::AbstractLayer) = throw(ArgumentError("Not Implemented"))

idims(self::AbstractLayer) = throw(ArgumentError("Not Implemented"))
odims(self::AbstractLayer) = throw(ArgumentError("Not Implemented"))

end ## Module

max_stars_repo_path: Data Science From Scratch/src/abstract_layers.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks