text:
"""
    mce(lxβ, y, mask_n [, n=length(y), mn=length(mask_n)])

For statistical design matrix `x`, model `b`, and binary response `y`, compute the mean classification error between a vector `lxβ = logistic(x*b)` and `y` with potentially masked entries indicated by the `0`s of the `Int` vector `mask_n`.
"""
function mce{T <: Float}(
    lxβ    :: DenseVector{T},
    y      :: DenseVector{T},
    mask_n :: DenseVector{Int},
    mn     :: Int = sum(mask_n)
)
    s = zero(T)
    @inbounds for i in eachindex(y)
        if mask_n[i] == 1
            s += abs(round(lxβ[i]) - round(y[i]))
        end
    end
    return s / mn
end

"""
    update_x!(x, x0, y, idxs, a)

If fed an `Int` vector `idxs` that sorts `x` by `abs()`, then `x` is updated as `x[idxs[1:k]] = x0[idxs[1:k]] - a*y[idxs[1:k]]`.
The output is comparable to

    copy!(x,x0)
    BLAS.axpy!(p,-a,y,1,x,1)

but `update_x!()` exploits the implied sparsity of `x` by only updating the top `k` components in magnitude.
"""
function update_x!{T <: Float}(
    x    :: DenseVector{T},
    x0   :: DenseVector{T},
    y    :: DenseVector{T},
    idxs :: DenseVector{Int},
    a    :: T
)
    fill!(x, zero(T))
    @inbounds for i in eachindex(idxs)
        idx    = idxs[i]
        x[idx] = x0[idx] - a*y[idx]
    end
    return nothing
end

"""
    df_norm(df, idxs, k0, k)

This subroutine efficiently computes `norm(df[idxs[k0:k]])`.
"""
function df_norm{T <: Float}(
    df   :: DenseVector{T},
    idxs :: DenseVector{Int},
    k0   :: Int,
    k    :: Int
)
    s = zero(T)
    @inbounds for i = k0:k
        idx = idxs[i]
        s += df[idx]*df[idx]
    end
    return sqrt(s)
end

"""
    logistic!(y, x)

This subroutine computes `y = logistic(x)` in-place, modifying `y`."
"""
function logistic!{T <: Float}(
	y :: DenseVecOrMat{T},
	x :: DenseVecOrMat{T}
)
    @assert size(y) == size(x)
    @inbounds for i in eachindex(x)
		y[i] = one(T) / (one(T) + exp(-x[i]))
	end
	return nothing
end

"""
    logistic_loglik(xβ, y, b, idxs, λ, k)

If fed a vector `idxs` of indices that sort `b`, then `logistic_loglik` computes the loglikelihood over the `k` active components of `b`.
"""
function logistic_loglik{T <: Float}(
	xβ   :: DenseVector{T},
    y    :: DenseVector{T},
	β    :: DenseVector{T},
	idxs :: DenseVector{Int},
	λ    :: T,
    k    :: Int
)
    n = length(xβ)
	s = zero(T)
	@inbounds for i = 1:n
		s += log(one(T) + exp(xβ[i])) - xβ[i]*y[i]
	end
    s /= n
	@inbounds for i = 1:k
		idx = idxs[i]
        s  += λ*β[idx]*β[idx] / 2
	end
	return s
end

function logistic_loglik{T <: Float}(
	xβ   :: DenseVector{T},
    y    :: DenseVector{T},
	β    :: DenseVector{T},
	idxs :: BitArray{1},
	λ    :: T,
    k    :: Int
)
    n = length(xβ)
	s = zero(T)
	@inbounds for i = 1:n
		s += log(one(T) + exp(xβ[i])) - xβ[i]*y[i]
	end
    s /= n
    ntrue = 0
    @inbounds for i in eachindex(idxs)
		idx = idxs[i]
        if idx
            ntrue += 1
            s  += λ*β[idx]*β[idx] / 2
        end
        ntrue >= k && break
	end
	return s
end

"""
    logistic_loglik(xβ, y, b, λ)

Compute the regularized negative logistic loglikelihood

   ( sum(log(1 + exp(x*β))) - y'*x*β) / n + 0.5*λ*sumabs2(β)

The unregularized loglikelihood is computed by setting `λ = 0.0`.

Arguments:

- `xβ = x*β`. Note that this must be precomputed; the function does not compute `x*b`.
- `y` is the response vector of `0`s and `1`s.
- `β` is the statistical model.
- `λ` is the regularization parameter.
"""
function logistic_loglik{T <: Float}(
	xβ :: DenseVector{T},
    y  :: DenseVector{T},
	β  :: DenseVector{T},
	λ  :: T
)
    n = length(xβ)
	s = zero(T)
	@inbounds for i = 1:n
		s += log(one(T) + exp(xβ[i])) - xβ[i]*y[i]
	end
    s /= n
    @inbounds for i in eachindex(β)
        s += λ*β[i]*β[i] / 2
	end
	return s
end

function logistic_loglik{T <: Float}(
	xβ     :: DenseVector{T},
    y      :: DenseVector{T},
	β      :: DenseVector{T},
    mask_n :: DenseVector{Int},
    m_val  :: Int,
	λ      :: T,
	k      :: Int;
    mn     :: Int = sum(mask_n),
)
	n = length(xβ),
	s = zero(T)
	@inbounds for i = 1:n
        if mask_n[i] != m_val
            s += log(one(T) + exp(xβ[i])) - xβ[i]*y[i]
        end
	end
#    s /= n
    s /= mn
	@inbounds for i = 1:k
        s += λ*β[i]*β[i] / 2
	end
	return s
end

function logistic_loglik{T <: Float}(
	xβ      :: DenseVector{T},
    y       :: DenseVector{T},
	b       :: DenseVector{T},
    idxs    :: DenseVector{Int},
    mask_n  :: DenseVector{Int},
    m_val   :: Int,
	λ  :: T,
	k       :: Int;
	n       :: Int = length(xb),
    mn      :: Int = sum(mask_n),
)
	s = zero(T)
	@inbounds for i = 1:n
        if mask_n[i] != m_val
            s += log(one(T) + exp(xβ[i])) - xβ[i]*y[i]
        end
	end
#    s /= n
    s /= mn
	@inbounds for i = 1:k
        idx = idxs[i]
        s += λ*b[idx]*b[idx] / 2
	end
	return s
end

function logistic_loglik{T <: Float}(
	xβ      :: DenseVector{T},
    y       :: DenseVector{T},
	b       :: DenseVector{T},
    idxs    :: BitArray{1},
    mask_n  :: DenseVector{Int},
    m_val   :: Int,
	λ  :: T,
	k       :: Int;
	n       :: Int = length(xb),
    mn      :: Int = sum(mask_n),
)
	s = zero(T)
	@inbounds for i = 1:n
        if mask_n[i] != m_val
            s += log(one(T) + exp(xβ[i])) - xβ[i]*y[i]
        end
	end
#    s /= n
    s /= mn
    nk = 0
	@inbounds for i = 1:length(idxs)
        idx = idxs[i]
        if idx
            s += λ*b[idx]*b[idx] / 2
            nk += 1
        end
        nk >= k && break
	end
	return s
end

"""
    logistic_grad!(df, lxβ, x, y, b, xβ, idxs, k, λ)

Compute the *partial* gradient of the regularized negative logistic loglikelihood as indexed by the `Int` vector `idxs`.
"""
function logistic_grad!{T <: Float}(
    df   :: DenseVector{T},
    lxβ  :: DenseVector{T},
    x    :: DenseMatrix{T},
    y    :: DenseVector{T},
    b    :: DenseVector{T},
    xβ   :: DenseVector{T},
    idxs :: DenseVector{Int},
    k    :: Int,
    λ    :: T;
)
    n = length(xβ)
    logistic!(lxβ,xβ)
    fill!(df, zero(T))
    @inbounds for i = 1:k
        idx = idxs[i]
#        df[idx] = zero(T)
        @inbounds for j = 1:n
            df[idx] += x[j,idx] * (lxβ[j] - y[j])
        end
        df[idx] /= n
        df[idx] += λ*b[idx]
    end
    return nothing
end

function logistic_grad!{T <: Float}(
    df   :: DenseVector{T},
    lxβ  :: DenseVector{T},
    x    :: DenseMatrix{T},
    y    :: DenseVector{T},
    b    :: DenseVector{T},
    xβ   :: DenseVector{T},
    idxs :: BitArray{1},
    k    :: Int,
    λ    :: T;
)
    n = length(xβ)
    logistic!(lxβ,xβ)
    fill!(df, zero(T))
    ntrue = 0
    @inbounds for i in eachindex(idxs)
        idx = idxs[i]
#        df[idx] = zero(T)
        if idx
            ntrue += 1
            @inbounds for j = 1:n
                df[idx] += x[j,idx] * (lxβ[j] - y[j])
            end
            df[idx] /= n
            df[idx] += λ*b[idx]
        end
        ntrue >= k && break
    end
    return nothing
end

### 18 OCT 2016:
### must redo logistic_grad! for revamped BEDFiles

"""
    logistic_grad!(df, lxβ, x::BEDFile, y, b, xβ, means, invstds, idxs, mask_n, k, λ [, n=length(xβ)])

Compute the *partial* gradient of the regularized negative logistic loglikelihood, as indexed by the `Int` vector `idxs`, with a `BEDFile` object `x`.
"""
function logistic_grad!{T <: Float}(
    df      :: DenseVector{T},
    lxβ     :: DenseVector{T},
    x       :: BEDFile,
    y       :: DenseVector{T},
    b       :: DenseVector{T},
    xβ      :: DenseVector{T},
    idxs    :: DenseVector{Int},
    mask_n  :: BitArray{1},
    k       :: Int,
    λ  :: T;
    mn      :: Int = sum(mask_n),
)
    n = length(y)
    logistic!(lxβ, xβ,n=n)
    fill!(df, zero(T))
    @inbounds for i = 1:k
        idx = idxs[i]
        m   = x.means[idx]
        d   = x.invstds[idx]
#        df[idx] = zero(T)
        @inbounds for j = 1:n
            if mask_n[i]
                df[idx] += (x[j,idx] - m) * d * (lxβ[j] - y[j])
            end
        end
#        df[idx] /= n
        df[idx] /= mn
        df[idx] += λ*b[idx]
    end
    return nothing
end

"""
    logistic_grad!(df, lxβ, x, y, β, xβ, λ)

Computes the full gradient of the negative logistic loglikelihood.
"""
function logistic_grad!{T <: Float}(
    df  :: DenseVector{T},
    lxβ :: DenseVector{T},
    x   :: DenseMatrix{T},
    y   :: DenseVector{T},
    β   :: DenseVector{T},
    xβ  :: DenseVector{T},
    λ   :: T;
)
    n = length(y)
    logistic!(lxβ, xβ)
    BLAS.axpy!(-one(T), y, lxβ)
    At_mul_B!(df, x, lxβ)
    scale!(df, 1/n)
    BLAS.axpy!(λ, β, df)
    return nothing
end

"""
    logistic_grad!(df, lxβ, x::BEDFile, y, b, xβ, means, invstds, λ, [, n=length(xβ), p=length(df)])

Computes the full gradient of the negative logistic loglikelihood with a `BEDFile` object `x`.
"""
function logistic_grad!{T <: Float}(
    df      :: SharedVector{T},
    lxβ     :: SharedVector{T},
    x       :: BEDFile{T},
    y       :: SharedVector{T},
    b       :: SharedVector{T},
    xβ      :: SharedVector{T},
    mask_n  :: BitArray{1},
    λ       :: T;
    pids    :: Vector{Int} = procs(x)
)
    n,p = size(x)
    mn  = sum(mask_n)
    logistic!(lxβ, xβ)
    BLAS.axpy!(-one(T), sdata(y), sdata(lxβ))
    At_mul_B!(df, x, lxβ, mask_n, pids=pids)
    scale!(df, 1/mn)
    BLAS.axpy!(λ, sdata(b), sdata(df))
    return nothing
end

"""
    log2xb!(lxβ, l2xb, xβ)

For `n`-vectors lxβ, l2xβ, and xβ, Compute `lxβ = logistic(xβ)` and `l2xb = lxβ*(1 - lxβ) / n` in-place and in one pass, overwriting `lxβ` and `l2xb`.
"""
function log2xb!{T <: Float}(
    lxβ  :: DenseVector{T},
    l2xβ :: DenseVector{T},
    xβ   :: DenseVector{T}
)
    n = length(lxβ)
    @inbounds for i = 1:n
        lxβ[i]  = one(T) / (one(T) + exp(-xβ[i]))
        l2xβ[i] = lxβ[i]*(one(T) - lxβ[i]) / n
    end
    return nothing
end


text:
module ImageComponentAnalysis

using Images

abstract type ComponentAnalysisAlgorithm end
struct ContourTracing <: ComponentAnalysisAlgorithm end

include("contour_tracing.jl")

export
    label_components,
	ContourTracing
end # module

ext: jl


text:
module mymodule
include("union_find.jl")

using Test

function test1()
    uf = UnionFind(5)
    @show uf
    # uf = UnionFind([-1, -1, -1, -1, -1], [1, 1, 1, 1, 1])
    @show unite!(uf, 2, 3)
    # unite!(uf, 2, 3) = true
    @show unite!(uf, 2, 3)
    # unite!(uf, 2, 3) = false
    @show uf
    # uf = UnionFind([-1, -1, 2, -1, -1], [1, 2, 2, 1, 1])
    @show root!(uf, 2)
    # root!(uf, 2) = 2
    @show issame(uf, 1, 2)
    # issame(uf, 1, 2) = false
    @show issame(uf, 2, 3)
    # issame(uf, 2, 3) = true
    @show size(uf, 2)
    # size(uf, 2) = 2
    @show size(uf, 3)
    # size(uf, 3) = 2
    @show members(uf, 1)
    # members(uf, 1) = [1]
    @show members(uf, 2)
    # members(uf, 2) = [2, 3]
    @show members(uf, 3)
    # members(uf, 3) = [2, 3]
    @show members(uf, 4)
    # members(uf, 4) = [4]
    @show roots(uf)
    # roots(uf) = [1, 2, 4, 5]
    @show group_size(uf)
    # group_size(uf) = 4
    @show all_group_members(uf)
    # all_group_members(uf) = Dict(5 => [5], 4 => [4], 2 => [2, 3], 1 => [1])

end

test1()

end


text:
using TriMeshes
using Test

@testset "TriMeshes.jl" begin
    @test areaXY(Vertex(0,0,0), Vertex(0,1,0), Vertex(1,0,0)) == 0.5
    @test magnitude(Vertex(0,1,0)) == magnitude(Vertex(1,0,0))
end


text:
using PyPlot
partial = true
include("pendulum.jl")

clf()
plot(X.tt, first.(X.yy), label=L"X_1", color=:darkviolet)
plot(X.tt, last.(X.yy), label=L"X_2", linewidth=0.2, color=:darkviolet)

plot(X.tt, first.(Xsm.yy), label=L"x_1", color=:darkgreen)
plot(X.tt, last.(Xsm.yy), ":", label=L"x_2", color=:darkgreen)

if partial
    plot(V.tt, V.yy, "*", label=L"V",  color=:darkviolet)
else
    plot(V.tt, first.(V.yy), ",", label=L"V_1",  color=:darkviolet)
    plot(V.tt, last.(V.yy), "*", label=L"V_2",  color=:darkviolet)
end
legend()


text:
println("Hi Web World from Julia!")

# Output:
# Hi Web World from Julia!


text:
# An axon takes an input from the soma and routes it to a destination without delay.
# a) Just poisson input
# b) Poisson and Stimulus

# ----------------------------------------------------------------
# Direct
# ----------------------------------------------------------------
mutable struct DirectAxon <: AbstractAxon
    base::BaseData{UInt8}
    len::Integer

    function DirectAxon()
        o = new()
        o.base = BaseData{UInt8}()
        o.len = 0.0
        o
    end
end

function set!(axon::DirectAxon, v::UInt8)
    axon.base.input = v
    axon.base.output = v
end

function input(axon::DirectAxon)
    axon.base.input
end

function output(axon::DirectAxon)
    axon.base.output
end

function reset!(axon::DirectAxon)
    axon.base.input = 0
    axon.base.output = 0
end

function step!(axon::DirectAxon)
    # Instantanious traversal down axon.
    axon.base.output = axon.base.input
end


text:
# Process any ASCII escape sequences in a raw string captured by string macro
function _unsafe_unescape!(c)
    n = ncodeunits(c)
    a = Base.unsafe_wrap(Array, pointer(c)::Ptr{UInt8}, n)
    for i = 1:n
        if a[i] == 0x5c # \
            if a[i+1] == 0x30 # \0
                a[i] = 0x00
                n = _advance!(a, i+1, n)
            elseif a[i+1] == 0x61 # \a
                a[i] = 0x07
                n = _advance!(a, i+1, n)
            elseif a[i+1] == 0x62 # \b
                a[i] = 0x08
                n = _advance!(a, i+1, n)
            elseif a[i+1] == 0x66 # \f
                a[i] = 0x0c
                n = _advance!(a, i+1, n)
            elseif a[i+1] == 0x6e # \n
                a[i] = 0x0a
                n = _advance!(a, i+1, n)
            elseif a[i+1] == 0x72 # \r
                a[i] = 0x0d
                n = _advance!(a, i+1, n)
            elseif a[i+1] == 0x74 # \t
                a[i] = 0x09
                n = _advance!(a, i+1, n)
            elseif a[i+1] == 0x76 # \v
                a[i] = 0x0b
                n = _advance!(a, i+1, n)
            elseif a[i+1] == 0x5c # \\
                a[i] = 0x5c
                n = _advance!(a, i+1, n)
            elseif a[i+1] == 0x27 # \'
                a[i] = 0x27
                n = _advance!(a, i+1, n)
            elseif a[i+1] == 0x22 # \"
                a[i] = 0x22
                n = _advance!(a, i+1, n)
            end
        end
    end
    return n
end

@inline function _advance!(a::AbstractArray{UInt8}, i::Int, n::Int)
    copyto!(a, i, a, i+1, n-i)
    a[n] = 0x00
    n -= 1
end


text:
using GeometricBase
using Test

@testset "$(rpad("Interface Definition",80))" begin
    struct TestSolution{dType, tType, N} <: Solution{dType, tType, N} end
    test_sol = TestSolution{Float64, Float64, 2}()

    @test_throws ErrorException nsave(test_sol)
    @test_throws ErrorException ntime(test_sol)
    @test_throws ErrorException nsamples(test_sol)

    @test_throws ErrorException counter(test_sol)
    @test_throws ErrorException offset(test_sol)
    @test_throws ErrorException lastentry(test_sol)
    @test_throws ErrorException timesteps(test_sol)

    @test_throws ErrorException eachtimestep(test_sol)
    @test_throws ErrorException eachsample(test_sol)
end


text:
module ArrayFire

using Compat

export AFArray

include("config.jl")

type AFArray{T,N} <: AbstractArray{T,N}
    ptr::Ptr{Void}
    function AFArray(ptr::Ptr{Void})
        a = new(ptr)
        finalizer(a, x -> af_release_array(x))
        a
    end
end

typealias AFVector{T} AFArray{T,1}
typealias AFMatrix{T} AFArray{T,2}

include("error.jl")
include("defns.jl")
include("indexing.jl")
include("wrap.jl")
include("utils.jl")
include("create.jl")
include("math.jl")
include("vector.jl")
include("backend.jl")
include("stats.jl")
include("linalg.jl")
include("image.jl")
include("cv.jl")
include("signal.jl")
include("device.jl")

end


text:
function coupleproximal!(system::LegSystem,n::Int64,h::Float64,uprox::Float64)
    # update right-running invariant
    system.branches.W1root = 2*uprox-system.branches.W2root;
    # update proximal A, Q w/ invariants
    system.branches.A[1][1,n+2] = ((2*system.solverparams.rho/
        system.branches.beta[1][end])^2*(0.125*(system.branches.W1root-
        system.branches.W2root) + system.branches.c0[1][end])^4);
    system.branches.Q[1][1,n+2] = (system.branches.A[1][1,n+2]*0.5*(
        system.branches.W1root + system.branches.W2root));
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
module Matdescra

using MKLSparse
using Test, SparseArrays, LinearAlgebra

sA = sprand(5, 5, 0.01)
sS = sA'sA
sTl = tril(sS)
sTu = triu(sS)

@test MKLSparse.BLAS.matdescra(Symmetric(sTl,:L)) == "SLNF"
@test MKLSparse.BLAS.matdescra(Symmetric(sTu,:U)) == "SUNF"
@test MKLSparse.BLAS.matdescra(LowerTriangular(sTl)) == "TLNF"
@test MKLSparse.BLAS.matdescra(UpperTriangular(sTu)) == "TUNF"
@test MKLSparse.BLAS.matdescra(UnitLowerTriangular(sTl)) == "TLUF"
@test MKLSparse.BLAS.matdescra(UnitUpperTriangular(sTu)) == "TUUF"
@test MKLSparse.BLAS.matdescra(sA) == "GUUF"

end  # module

module BLAS

using MKLSparse
using Test, LinearAlgebra, SparseArrays

macro test_blas(ex)
    return quote
        MKLSparse.BLAS.__counter[] = 0
        @test $(esc(ex))
        @test MKLSparse.BLAS.__counter[] == 1
    end
end

@testset "matrix-vector multiplication (non-square)" begin
    for i = 1:5
        a = sprand(10, 5, 0.5)
        b = rand(5)
        @test_blas maximum(abs.(a*b - Array(a)*b)) < 100*eps()
        b = rand(5, 5)
        @test_blas maximum(abs.(a*b - Array(a)*b)) < 100*eps()
        b = rand(10)
        @test_blas maximum(abs.(a'*b - Array(a)'*b)) < 100*eps()
        @test_blas maximum(abs.(transpose(a)*b - Array(a)'*b)) < 100*eps()
        b = rand(10,10)
        @test_blas maximum(abs.(a'*b - Array(a)'*b)) < 100*eps()
        @test_blas maximum(abs.(transpose(a)*b - Array(a)'*b)) < 100*eps()
    end
end

#?
@testset "complex matrix-vector multiplication" begin
    for i = 1:5
        a = I + im * 0.1*sprandn(5, 5, 0.2)
        b = randn(5,3) + im*randn(5,3)
        c = randn(5) + im*randn(5)
        d = randn(5) + im*randn(5)
        α = rand(ComplexF64)
        β = rand(ComplexF64)
        @test_blas (maximum(abs.(a*b - Array(a)*b)) < 100*eps())
        @test_blas (maximum(abs.(a'*b - Array(a)'*b)) < 100*eps())
        @test_blas (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*eps())
        @test_blas (maximum(abs.(mul!(similar(b), a, b) - Array(a)*b)) < 100*eps())
        @test_blas (maximum(abs.(mul!(similar(c), a, c) - Array(a)*c)) < 100*eps())
        @test_blas (maximum(abs.(mul!(similar(b), transpose(a), b) - transpose(Array(a))*b)) < 100*eps())
        @test_blas (maximum(abs.(mul!(similar(c), transpose(a), c) - transpose(Array(a))*c)) < 100*eps())

        c = randn(6) + im*randn(6)
        @test_throws DimensionMismatch transpose(a)*c
        @test_throws DimensionMismatch a.*c
        @test_throws DimensionMismatch a.*c
    end
end

@testset "triangular" begin
    n = 100
    A = sprandn(n, n, 0.5) + sqrt(n)*I
    b = rand(n)
    symA = A + transpose(A)
    trilA = tril(A)
    triuA = triu(A)
    trilUA = tril(A, -1) + I
    triuUA = triu(A, 1)  + I

    @test_blas LowerTriangular(trilA) \ b ≈ Array(LowerTriangular(trilA)) \ b
    @test_blas LowerTriangular(trilA) * b ≈ Array(LowerTriangular(trilA)) * b

    @test_blas UpperTriangular(triuA) \ b ≈ Array(UpperTriangular(triuA)) \ b
    @test_blas UpperTriangular(triuA) * b ≈ Array(UpperTriangular(triuA)) * b

    @test_blas UnitLowerTriangular(trilUA) \ b ≈ Array(UnitLowerTriangular(trilUA)) \ b
    @test_blas UnitLowerTriangular(trilUA) * b ≈ Array(UnitLowerTriangular(trilUA)) * b

    @test_blas UnitUpperTriangular(triuUA) \ b ≈ Array(UnitUpperTriangular(triuUA)) \ b
    @test_blas UnitUpperTriangular(triuUA) * b ≈ Array(UnitUpperTriangular(triuUA)) * b

    @test_blas Symmetric(symA) * b ≈ Array(Symmetric(symA)) * b
end

end # module


text:
ENV["GKSwstype"] = "nul"
using Plots
using SpikingNeuralNetworks
SNN.@load_units

S = SNN.Rate(;N = 200)
SS = SNN.PINningSynapse(S, S; σ = 1.5, p = 1.0)
P, C = [S], [SS]

SNN.monitor(SS, [(:g, [1])])

A = 1.3 / 1.5; fr = 1 / 60ms
f(t) = (A /1.0) * sin(1π * fr * t) +  (A / 2.0) * sin(2π * fr * t) +
(A / 6.0) * sin(3π * fr * t) +   (A / 3.0) * sin(4π * fr * t)

ts = 0:0.1ms:1440ms
for (i, t) in enumerate(ts)
    SS.f .= [f(t); SS.g[2:end]]
    SpikingNeuralNetworks.train!(P, C, 0.1ms, t)
end

SNN.vecplot(SS, :g); plot!(f.(ts))
png("results/pinning_net")


text:
"""
    linesfill(xs, ys; lower, upper, kwargs...)

Line plot with a shaded area between `lower` and `upper`. If `lower` and `upper`
are not given, shaded area is between `0` and `ys`.
## Attributes
$(ATTRIBUTES)
"""
@recipe(LinesFill) do scene
    l_theme = default_theme(scene, Lines)
    Attributes(
        color=l_theme.color,
        colormap=l_theme.colormap,
        colorrange=get(l_theme.attributes, :colorrange, automatic),
        linestyle=l_theme.linestyle,
        linewidth=l_theme.linewidth,
        fillalpha=0.15,
        lower=automatic,
        upper=automatic,
    )
end

function Makie.plot!(p::LinesFill)
    lines!(p, p[1:2]...;
        color = p.color,
        linestyle = p.linestyle,
        linewidth = p.linewidth,
        colormap = p.colormap,
        colorrange = p.colorrange,
    )
    lower = lift(p[:lower], p[2]) do lower, y
        return lower === automatic ? zero(y) : lower
    end
    upper = lift(p[:upper], p[2]) do upper, y
        return upper === automatic ? y : upper
    end
    meshcolor = lift(to_color∘tuple, p.color, p.fillalpha)
    band!(p, p[1], lower, upper; color = meshcolor)
end


text:
"""
    SinkFreeGraph{T<:LG.SimpleDiGraph{Int64}} <: AbstractGraphPointProcess{T}

Concrete type representing a point process on the edges of a `graph` characterizing the uniform distribution on the orientations of the edges conditioned on the absence of sinks.

**See also**

- Section 4.1 of [GuJeLi19](@cite).

# Example

A realization from a ``5\\times 5`` grid graph

![assets/sink_free_graph.png](assets/sink_free_graph.png)
"""
struct SinkFreeGraph{T<:LG.SimpleDiGraph{Int64}} <: AbstractGraphPointProcess{T}
    "Graph"
    graph::LG.SimpleGraph{Int64}
end

function Base.show(io::IO, pp::SinkFreeGraph{T}) where {T}
    print(io, "SinkFreeGraph{$T}\n- graph = $(pp.graph)")
end

"""
    SinkFreeGraph(graph::LG.SimpleGraph{T}) where {T<:Int}

Construct a [`PRS.SinkFreeGraph`](@ref).

```jldoctest; output = true
using PartialRejectionSampling
using LightGraphs; const LG = LightGraphs

sfg = PRS.SinkFreeGraph(LG.grid([5, 5]))

# output

SinkFreeGraph{LightGraphs.SimpleGraphs.SimpleDiGraph{Int64}}
- graph = {25, 40} undirected simple Int64 graph
```
"""
SinkFreeGraph(graph::LG.SimpleGraph{Int64}) = SinkFreeGraph{LG.SimpleDiGraph{Int64}}(graph)

"""
    generate_sample([rng::Random.AbstractRNG,] pp::SinkFreeGraph)

Generate an exact sample from the [`PRS.SinkFreeGraph`](@ref).

Default sampler is [`PRS.generate_sample_prs`](@ref).
"""
generate_sample(rng::Random.AbstractRNG, pp::SinkFreeGraph) = generate_sample_prs(rng, pp)
# generate_sample(pp::SinkFreeGraph) = generate_sample_prs(Random.default_rng(), pp)

"""
    generate_sample_prs(
        [rng::Random.AbstractRNG,]
        pp::SinkFreeGraph{T}
    )::T where {T}

Generate an orientated version of `pp.graph` uniformly at random among all possible orientations conditioned on the absence of sinks, using Partial Rejection Sampling (PRS).

**See also**

- Section 4.1 of [GuJeLi19](@cite).

# Example

A illustration of the procedure on a ``5 \\times 5`` grid grah.

- red: variables involved in constraints violation (orientation of the edges pointing to a sink node)
- orange: variables resampled (orientation of the previous edges)

![assets/sink_free_graph_prs.gif](assets/sink_free_graph_prs.gif)
"""
function generate_sample_prs(
    rng::Random.AbstractRNG,
    pp::SinkFreeGraph{T}
)::T where {T}
    g = random_edge_orientation(rng, pp.graph)
    while true
        sinks = sink_nodes(g)
        isempty(sinks) && break
        # Resample i.e. flip edges forming a sink uniformly at random (Bernoulli(0.5))
        for v in sinks
            in_neighbors_v = copy(LG.inneighbors(g, v))
            for w in in_neighbors_v
                if rand(rng) < 0.5
                    LG.rem_edge!(g, w, v)
                    LG.add_edge!(g, v, w)
                end
            end
        end
    end
    return g
end

# generate_sample_prs(pp::SinkFreeGraph) = generate_sample_prs(Random.default_rng(), pp)


text:
include("utils.jl")

m = 1
n = 21

a = constant(0.01*Array(LinRange(1,2,n+1)))
b = constant(0.02*Array(LinRange(1,2,n+1)))
σn = 50*Array(LinRange(1,2,n+1))
τ = 30*Array(LinRange(1,2,n+1))
η = 4.7
v0 = 1e-6
f0 = 0.6
Ψ0 = Array(LinRange(1,2,n+1)) * 1.01 * f0
Dc = 0.01

u_left_bd_S2 = Array(LinRange(1,2,n+1))
v_bd_S2 = rate_state_friction(a, u_left_bd_S2, v0, Ψ0, σn, τ, η, 1.0) # t+ Δt + δΔt

# du is defined on the left side
function StateRate(Ψ, v)
    v = abs(v)
    s = v/(2v0)
    o = exp(-Ψ/a)
    f = a * (Ψ/a + log(s+sqrt(s^2+o^2)))
    # f = a* tf.asinh(v/(2v0)*exp(Ψ/a))
    Δψ = -v/Dc*(f-f0+(b-a)*log(v/v0))
end

dpsi_S3  = StateRate(Ψ0, v_bd_S2) # t + Δt

sess = Session();
run(sess, v_bd_S2)

# run(sess, dpsi_S3)

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
import SimpleDrawing: newdraw, finish
export newdraw, finish, draw

"""
`draw(X)` draws the hyperbolic object `X` in a graphics window.

`draw` may be applied to a list of hyperbolic objects.

The typical sequence of drawing starts by clearing the screen with the `plot()` function
(from the `Plots` module), then various to calls to `draw` and then
concludes with a call to `finish()` (see the help message for that function).
"""
function SimpleDrawing.draw(P::HPoint)
    x,y = reim(getz(P))
    draw_point(x,y;P.attr...)
end

function SimpleDrawing.draw(S::HSegment)
    P,Q = endpoints(S)

    if P==Q
        return
    end

    M = midpoint(S)

    z = getz(P)
    zz = getz(Q)
    w = getz(M)

    if abs(imag( (z-w)/(zz-w) )) <= THRESHOLD*eps(1.0)   # they're linear
        draw_segment(z,zz;S.attr...)
    else
        draw_arc(z,w,zz; S.attr...)
    end
end

function SimpleDrawing.draw(L::HLine)
    t1 = L.s
    t2 = L.t

    if t1==t2
        return
    end

    x = cos(t1)
    y = sin(t1)
    xx = cos(t2)
    yy = sin(t2)

    if abs(abs(t1-t2)-pi) <= THRESHOLD*eps(1.0)
        draw_segment(x,y,xx,yy;L.attr...)
    else
        P = point_on_line(L)
        m = getz(P)
        draw_arc(x+im*y, m, xx+im*yy; L.attr...)
    end
end

function SimpleDrawing.draw(R::HRay)
    t = R.t
    A = get_vertex(R)
    B = point_on_ray(R)

    a = getz(A)
    b = getz(B)
    c = exp(im*t)

    z = (b-a)/(c-a)  # if this is real, we draw a segment
    if abs(imag(z)) <= THRESHOLD*eps(1.0)
        draw_segment(a,c;R.attr...)
    else
        draw_arc(a,b,c;R.attr...)
    end
end

function SimpleDrawing.draw(HP::HPlane)
    draw_circle(0,0,1; HP.attr...)
end

function SimpleDrawing.draw(list::Array{T,1}) where T <: HObject
    g = plot!()
    for X in list
        g = draw(X)
    end
    return g
end

# SimpleDrawing.draw(args...) = draw(collect(args))

function SimpleDrawing.draw(HO::HObject, args...)
    draw(HO)
    draw(collect(args))
end

function SimpleDrawing.draw(C::HCircle)
    X,Y,Z= points_on_circle(C)
    draw_circle(getz(X),getz(Y),getz(Z);C.attr...)
end

function SimpleDrawing.draw(HC::Horocycle)
    c = euclidean_center(HC)
    r = abs(getz(HC.pt)-c)
    draw_circle(c,r;HC.attr...)
end

function SimpleDrawing.draw(X::Union{HPolygon,HTriangle})
    for S in sides(X)
        copy_attr(S,X)
        draw(S)
    end
    draw()
end

SimpleDrawing.draw(C::HContainer) = draw(collect(C))


text:
using Test
using SafeTestsets

@safetestset "Example Tests" begin
    include("example_test.jl")
end


text:
module MyMeshIO

f(x) = 2x

end

max_stars_repo_name: powerzf/MyMeshIO

max_issues_repo_name: powerzf/MyMeshIO

max_forks_repo_name: powerzf/MyMeshIO


text:
# baiaiyu
# 2021,Jul,15
# Many-Body Localization Transition https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.119.075702
# programe for application
# 请君洗耳 听我之言

# bitstring() 获取二进制结构
# 稀疏矩阵
# ↓ 实对称矩阵在julia
# https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/12
# ↓ LAPACK在julia
# LinearAlgebra.LAPACK.syevr!

module MBL_ED

export cal_one_task, primitive_run, give_tasks, auto_manger, give_empty_dataframe_for_format

using SparseArrays  # 稀疏矩阵支持
using LinearAlgebra # QR法得特征值和特征向量
# using KrylovKit     # lanczos 直接得到特征值和特征向量
# using Arpack
using Statistics    # 平均值与标准差
# using Gadfly        # 绘图
using Distributed   #   并行
using JLD   # 保存数据
using DataFrames    # 数据整合
using Dates # 时间

######### GPU ##########
# using CUDA.CUSOLVER
# using CUDA
########################

include("./my_hamiltonion.jl")
include("./phy_quanlities.jl")
include("./mymath.jl")
include("./parallel_manager.jl")    # 并行计算管理器
include("./mydata.jl")    # 数据整合

# 封装计算任务
struct QPTask
    H::Matrix{Float64}
    states::Vector{Int}
    left_chain::Vector{Int}
    right_chain::Vector{Int}
    nhil::Int
    n_tol::Int
    W::Float64
    omega::Float64
    phi::Float64
end

struct RTask
    H::Matrix{Float64}
    states::Vector{Int}
    left_chain::Vector{Int}
    right_chain::Vector{Int}
    nhil::Int
    n_tol::Int
    W::Float64
    omega::Float64
end

function cal_one_task(one_R_task::RTask)
    St = 0.5*(one_R_task.n_tol*log(2) - 1)
    # HH = Array{Float64}(undef, one_R_task.nhil, one_R_task.nhil)
    # HH .= one_R_task.H
    HH = copy(one_R_task.H)
    wave = fill_Sz(one_R_task.states, one_R_task.n_tol, one_R_task.nhil, one_R_task.W, one_R_task.omega)
    for i in 1:one_R_task.nhil
        @inbounds HH[i, i] = HH[i, i] + wave[i]
    end

    empty_df = give_empty_dataframe_for_format()

    #################################
    # CPU
    eig_sol = eigen(Symmetric(HH))
    push!( empty_df, ["R", string(one_R_task.n_tol), one_R_task.W, lsr(eig_sol.values), EE(eig_sol.vectors, one_R_task.left_chain, one_R_task.right_chain, one_R_task.nhil)/St] )
    #################################

    #################################
    # GPU
    # eig_sol = CUSOLVER.syevd!( 'V', 'U', CuArray(Symmetric(HH)) )
    # #"R", one_R_task.n_tol, one_R_task.W, lsr(Array(eig_sol[1])), EE(Array(eig_sol[2]), one_R_task.left_chain, one_R_task.right_chain, one_R_task.nhil)/St
    # push!( empty_df, ["R", string(one_R_task.n_tol), one_R_task.W, lsr(eig_sol.values), EE(eig_sol.vectors, one_R_task.left_chain, one_R_task.right_chain, one_R_task.nhil)/St] )
    #################################
end

function cal_one_task(one_QP_task::QPTask)
    St = 0.5*(one_QP_task.n_tol*log(2) - 1)
    # HH = Array{Float64}(undef, one_QP_task.nhil, one_QP_task.nhil)
    # HH .= one_QP_task.H
    HH = copy(one_QP_task.H)
    wave = fill_Sz(one_QP_task.states, one_QP_task.n_tol, one_QP_task.nhil, one_QP_task.W, one_QP_task.omega, one_QP_task.phi)
    for i in 1:one_QP_task.nhil
        @inbounds HH[i, i] = HH[i, i] + wave[i]
    end

    empty_df = give_empty_dataframe_for_format()

    #################################
    # CPU
    eig_sol = eigen(Symmetric(HH))
    push!(empty_df, ["QP", string(one_QP_task.n_tol), one_QP_task.W, lsr(eig_sol.values), EE(eig_sol.vectors, one_QP_task.left_chain, one_QP_task.right_chain, one_QP_task.nhil)/St] )
    #################################

    #################################
    # GPU
    # eig_sol = CUSOLVER.syevd!( 'V', 'U', CuArray(Symmetric(HH)) )
    # #"QP", one_QP_task.n_tol, one_QP_task.W, lsr(Array(eig_sol[1])), EE(Array(eig_sol[2]), one_QP_task.left_chain, one_QP_task.right_chain, one_QP_task.nhil)/St
    # push!(empty_df, ["QP", one_QP_task.n_tol, one_QP_task.W, lsr(eig_sol.values), EE(eig_sol.vectors, one_QP_task.left_chain, one_QP_task.right_chain, one_QP_task.nhil)/St] )
    #################################
end

function primitive_run(
            ; tol_n::Int = 10
            , up::Int = 5
            , J = 1.0
            , Jp = 1.0
            , Jz = 1.0
            , W = 0.5:0.5:8.0
            , k = (sqrt(5) - 1)/2
            , samples::Int = 10
            )

    omega = 2*pi*k
    # Hilbert space 维数
    nhil::Int = div(factorial(big(tol_n)), (factorial(big(up))*factorial(big(tol_n - up))))
    states = Array{Int}(undef, nhil)
    # 得到所有要求的状态
    all_required_states!(states, nhil, tol_n, up)
    # 得到左半链和右半链的状态
    left_states = Array{Int}(undef, nhil)
    right_states = Array{Int}(undef, nhil)
    for i in 1:nhil
        left_states[i], right_states[i] = left_right_half_chain(states[i], tol_n)
    end
    # 得到不含随机项的哈密顿量
    H = fill_H(nhil, states, tol_n, J, Jp, Jz)

    cal_tasks_array = []
    for w in W
        for i in 1:samples
            push!(cal_tasks_array, RTask(H, states, left_states, right_states, nhil, tol_n, w, omega))
            push!(cal_tasks_array, QPTask(H, states, left_states, right_states, nhil, tol_n, w, omega, (2.0*pi*(rand() - 0.5))))
        end
    end

    # println("getalltasks")
    auto_manger(cal_tasks_array)

    # one_R_task = RTask(H, states, left_states, right_states, nhil, tol_n, W, omega)
    # cal_one_task(cal_tasks_array[2])
    # cal_tasks_array
end

"""

    give_tasks(
    ; tol_n::Int = 10
    , up::Int = 5
    , J = 1.0
    , Jp = 1.0
    , Jz = 1.0
    , W = 0.5:0.5:8.0
    , k = (sqrt(5) - 1)/2
    , samples::Int = 10
    )

DOI: 10.1103/PhysRevLett.119.075702
给出如下哈密顿量对应的计算任务, 一半为R模型一半为QP模型, 非周期性边界条件

``
\\begin{aligned}
H^{\\mathrm{QP} / R}=& J \\sum_{i=1}^{L-1}\\left(S_{i}^{x} S_{i+1}^{x}+S_{i}^{y} S_{i+1}^{y}\\right)+J_{z} \\sum_{i=1}^{L-1} S_{i}^{z} S_{i+1}^{z} \\\\
&+\\sum_{i=1}^{L} W \\cos \\left(2 \\pi k i+\\phi_{i}^{\\mathrm{QP} / R}\\right) S_{i}^{z} \\\\
&+J^{\\prime} \\sum_{i=1}^{L-2}\\left(S_{i}^{x} S_{i+2}^{x}+S_{i}^{y} S_{i+2}^{y}\\right)
\\end{aligned}
``
# Arguments
- `tol_n::Int`: 一维链长度
- `up::Int`: Sz向上的数量
- `samples::Int`: 对于每组参数的采样数

# Examples

```julia-repl
julia> tasks = give_tasks();
julia> tasks[1]
MBL_ED.RTask([2.75 0.5 … 0.0 0.0; 0.5 1.75 … 0.0 0.0; … ; 0.0 0.0 … 1.75 0.5; 0.0 0.0 … 0.5 2.75], [31, 47, 55, 59, 61, 62, 79, 87, 91, 93  …  15361, 15362, 15364, 15368, 15376, 15392, 15424, 15488, 15616, 15872], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  121, 121, 121, 121, 121, 121, 121, 122, 123, 125], [32, 48, 56, 60, 62, 63, 80, 88, 92, 94  …  2, 3, 5, 9, 17, 33, 65, 1, 1, 1], 2002, 14, 0.5, 3.883222077450933)

julia> tasks[2]
MBL_ED.QPTask([2.75 0.5 … 0.0 0.0; 0.5 1.75 … 0.0 0.0; … ; 0.0 0.0 … 1.75 0.5; 0.0 0.0 … 0.5 2.75], [31, 47, 55, 59, 61, 62, 79, 87, 91, 93  …  15361, 15362, 15364, 15368, 15376, 15392, 15424, 15488, 15616, 15872], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  121, 121, 121, 121, 121, 121, 121, 122, 123, 125], [32, 48, 56, 60, 62, 63, 80, 88, 92, 94  …  2, 3, 5, 9, 17, 33, 65, 1, 1, 1], 2002, 14, 0.5, 3.883222077450933, 2.0325482732260935)
```
"""
function give_tasks(
        ; tol_n::Int = 10
        , up::Int = 5
        , J = 1.0
        , Jp = 1.0
        , Jz = 1.0
        , W = 0.5:0.5:8.0
        , k = (sqrt(5) - 1)/2
        , samples::Int = 10
        )

    omega = 2*pi*k
    # Hilbert space 维数
    nhil::Int = div(factorial(big(tol_n)), (factorial(big(up))*factorial(big(tol_n - up))))
    states = Array{Int}(undef, nhil)
    # 得到所有要求的状态
    all_required_states!(states, nhil, tol_n, up)
    # 得到左半链和右半链的状态
    left_states = Array{Int}(undef, nhil)
    right_states = Array{Int}(undef, nhil)
    for i in 1:nhil
    left_states[i], right_states[i] = left_right_half_chain(states[i], tol_n)
    end
    # 得到不含随机项的哈密顿量
    H = fill_H(nhil, states, tol_n, J, Jp, Jz)

    cal_tasks_array = []
    for w in W
        for i in 1:samples
            push!(cal_tasks_array, RTask(H, states, left_states, right_states, nhil, tol_n, w, omega))
            push!(cal_tasks_array, QPTask(H, states, left_states, right_states, nhil, tol_n, w, omega, (2.0*pi*(rand() - 0.5))))
        end
    end

    # auto_manger(cal_tasks_array)

    # one_R_task = RTask(H, states, left_states, right_states, nhil, tol_n, W, omega)
    # cal_one_task(cal_tasks_array[2])
    cal_tasks_array
end

# 准备不含随机项的哈密顿矩阵等计算所需数据与模型参数
# 包括tol_n, up, nhil
function prepare_H(
    ; tol_n::Int = 10
    , up::Int = 5
    , J = 1.0
    , Jp = 1.0
    , Jz = 1.0
    )

    # Hilbert space 维数
    nhil::Int = div(factorial(big(tol_n)), (factorial(big(up))*factorial(big(tol_n - up))))
    states = Array{Int}(undef, nhil)
    # 得到所有要求的状态
    all_required_states!(states, nhil, tol_n, up)
    # 得到左半链和右半链的状态
    left_states = Array{Int}(undef, nhil)
    right_states = Array{Int}(undef, nhil)
    for i in 1:nhil
    left_states[i], right_states[i] = left_right_half_chain(states[i], tol_n)
    end
    # 得到不含随机项的哈密顿量
    H = fill_H(nhil, states, tol_n, J, Jp, Jz)

    file_name = "$tol_n" * "_$up"* "_.jld"
    save(file_name, "H", H, "states", states, "left_chain", left_states,"right_chain", right_states,"tol_n", tol_n, "up", up, "nhil", nhil)

    nothing
    # H
end

# 使用以保存的哈密顿量计算
struct spTask
    Hjld::String
    omega::Float64
    W::Float64
    phi::Float64
end

# 使用已经储存的哈密顿量计算
function cal_one_sptask(one_sptask::spTask)

    omega = 2*pi*k

    # 载入数据
    data = load(one_sptask.Hjld)

    cal_one_task(
        QPTask(data["H"], data["states"], data["left_chain"], data["right_chain"], data["nhil"], data["n_tol"], W, one_sptask.omega, (2.0*pi*(rand() - 0.5)))
    )

end

function run_sptasks(::Vector{String},)

end

end # module end


text:
module CellwiseEditingDifferentiation

using Distributions
using Lora

export
  logit,
  inv_logit,
  beta_a_from_mv,
  beta_b_from_mv,
  beta_pars_from_mv,
  vpcprior,
  logdvdw,
  wpcprior

include("logit.jl")
include("beta_moments.jl")
include("vpcprior.jl")

end


text:
using TestsetScope
using Test

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


text:
module Agents

using Random
using StatsBase: mean, sample, pweights
using Ramnet.Encoders

using ..Environments: AbstractEnvironment, observation_type, observation_length, observation_extrema, action_type, action_set, action_length
import ..reset!, ..update!, ..select_action

export AbstractAgent, select_action, select_action!, observe!, update!, agent, encoding, make_agent

# TODO: Make this type parametric. All agents should know the types of its observations,
#       actions and encoder
abstract type AbstractAgent{O<:AbstractVector,A} end

# I can make use of multiple dispatch here so that both observe_first! and observe!
# are called observe!. The "observe_first!" equivalent only takes a observation while
# the "observe!" equivalent takes an action, reward and next observation
# function observe_first! end
# function select_action end
function observe! end
# function update! end

# NOTE: This version assumes discriminators take in as input binary patterns, already encoded
# function q_values!(agent::G, observation::O, dest::AbstractVector{Float64}) where {O <: AbstractVector,A <: Real,G <: AbstractAgent{O,A}}
#     q_max = typemin(Float64)
#     q = 0.0
#     encoded_obs = encode(agent.encoder, observation) # TODO: Don't have to create a bunch of temp arrays here
#     for (i, action) in Iterators.zip(eachindex(dest), agent.actions)
#         q = predict(agent.Q̂[action], encoded_obs)
#         dest[i] = q
#         q_max = q > q_max ? q : q_max
#     end

#     return q_max
# end

# TODO: Make a version that takes a destination vector
# TODO: This implementation is wrong. It must take into account the epsilon prob.
#       Refer to the Alberta MOOC exercise.
function expected_q_value(agent::G, observation::O) where {O<:AbstractVector,A<:Real,G<:AbstractAgent{O,A}}
    values = similar(agent.actions, Float64)
    q_max = q_values!(agent, observation, values)

    return mean(values, pweights((values .== q_max) .+ agent.ϵ))
end

include("RandomAgent.jl")
export RandomAgent

# include("WNN/MonteCarloDiscountedDiscriminatorAgent.jl")
# export MonteCarloDiscountedDiscriminatorAgent

include("WNN/SARSADiscountedDiscriminatorAgent.jl")
export SARSADiscountedDiscriminatorAgent

# include("WNN/ExpectedSARSADiscriminatorAgent.jl")
# export ExpectedSARSADiscriminatorAgent

# include("WNN/QLearningDiscountedDiscriminatorAgent.jl")
# export QLearningDiscountedDiscriminatorAgent

include("WNN/MCDiscriminatorAgent.jl")
export MCDiscriminatorAgent

# include("WNN/MCDifferentialDiscriminatorAgent.jl")
# export MCDifferentialDiscriminatorAgent

include("WNN/PG/FunctionalPG.jl")
export FunctionalPG

include("WNN/PG/ContinousFunctionalPG.jl")
export ContinousFunctionalPG

include("WNN/PG/DiscreteFunctionalPG.jl")
export DiscreteFunctionalPG

include("WNN/PG/FunctionalAC.jl")
export FunctionalAC

include("WNN/PG/BinaryRegularizedAC.jl")
export BinaryRegularizedAC

include("WNN/PG/BinaryAdaptiveAC.jl")
export BinaryAdaptiveAC

include("WNN/PG/DiscreteFunctionalAC.jl")
export DiscreteFunctionalAC

include("WNN/PG/ContinousFunctionalAC.jl")
export ContinousFunctionalAC

include("WNN/PG/ContinousRegularizedAC.jl")
export ContinousRegularizedAC

# include("WNN/PG/AltContinousFunctionalAC.jl")
# export AltContinousFunctionalAC

include("WNN/PG/ContinousAdaptiveAC.jl")
export ContinousAdaptiveAC

const agent_table = Dict{Symbol,Type{<:AbstractAgent}}(
    :MCDiscriminatorAgent => MCDiscriminatorAgent,
    # :MonteCarloDiscountedDiscriminatorAgent => MonteCarloDiscountedDiscriminatorAgent,
    :SARSADiscountedDiscriminatorAgent => SARSADiscountedDiscriminatorAgent,
    # :ExpectedSARSADiscriminatorAgent => ExpectedSARSADiscriminatorAgent,
    # :QLearningDiscountedDiscriminatorAgent => QLearningDiscountedDiscriminatorAgent
    :FunctionalPG => FunctionalPG,
    :FunctionalAC => FunctionalAC,
    :BinaryRegularizedAC => BinaryRegularizedAC,
    :BinaryAdaptiveAC => BinaryAdaptiveAC,
    :DiscreteFunctionalAC => DiscreteFunctionalAC,
    :ContinousFunctionalPG => ContinousFunctionalPG,
    :DiscreteFunctionalPG => DiscreteFunctionalPG,
    :ContinousFunctionalAC => ContinousFunctionalAC,
    # :ContinousFunctionalAC => ContinousFunctionalAC,
    :ContinousRegularizedAC => ContinousRegularizedAC,
    # :AltContinousFunctionalAC => AltContinousFunctionalAC,
    :ContinousAdaptiveAC => ContinousAdaptiveAC
)

const encoding_table = Dict{Symbol,Type{<:AbstractEncoder}}(
    :Thermometer => Thermometer,
    :CircularThermometer => CircularThermometer
)

function encoding(name::String, minimum, maximum, resolution)
    encoding_table[Symbol(name)](eltype(minimum), minimum, maximum, resolution)
end

function agent(name::String, ::Type{O}, ::Type{A}, actions, obs_size, encoder::E; kargs...) where {O,A,E<:AbstractEncoder}
    agent_table[Symbol(name)]{O,A,E}(actions, obs_size, encoder; kargs...)
end

function agent(agent_spec::Dict{Symbol,Any}, env::AbstractEnvironment)
    enc_spec = agent_spec[:encoding]
    enc = encoding(
        enc_spec[:name],
        observation_extrema(env)...,
        enc_spec[:resolution]
    )

    agent(
        agent_spec[:name],
        observation_type(env),
        action_type(env),
        action_set(env),
        observation_length(env),
        enc;
        agent_spec[:args]...
    )
end

function make_agent(env_type::Type{<:AbstractEnvironment}, agentspec::Dict{Symbol,Any}; seed::Union{Nothing,UInt} = nothing)
    encoder_name = agentspec[:encoding][:name] |> Symbol
    encoder_resolution = agentspec[:encoding][:resolution]

    min, max = observation_extrema(env_type)

    if haskey(agentspec[:encoding], :min)
        spec_min = agentspec[:encoding][:min]

        if spec_min isa Real # Scalar?
            min = similar(min)
            fill!(min, spec_min)
        else
            # Should check if spec_min is the right shape and type. Raise an
            # exception if not
            min = spec_min
        end
    end

    if haskey(agentspec[:encoding], :max)
        spec_max = agentspec[:encoding][:max]

        if spec_max isa Real # Scalar?
            max = similar(max)
            fill!(max, spec_max)
        else
            # Should check if spec_max is the right shape and type. Raise an
            # exception if not
            max = spec_max
        end
    end

    enc = encoding_table[encoder_name](min, max, encoder_resolution)

    agent_name = agentspec[:name] |> Symbol

    agent_table[agent_name](env_type, enc; agentspec[:args]..., seed)
end

# TODO: Make this a function of the Q-value approximation, not the agent.
#       (Define somewhere else a DiscreteQDiscriminator. Maybe in a Approximators
#        module)
# function select_action!(agent::MonteCarloDiscountedDiscriminatorAgent{O,A,E}, observation::O; ϵ::Float64=agent.ϵ) where {O <: AbstractVector,A <: Real,E <: AbstractEncoder}
#     # TODO: I could avoid creating these temporary arrays here
#     values = similar(agent.actions, Float64)
#     q_max = q_values!(agent, observation, values)

#     if rand(agent.rng) < ϵ # Random action taken with ϵ probability
#         agent.action = rand(agent.rng, agent.actions)
#     else
#         # Selecting action that maximizes Q̂, with ties broken randomly
#         agent.action = sample(agent.rng, agent.actions, pweights(values .== q_max))
#     end

#     return agent.action::A
# end

end


text:
module Wavefronts

export intensity, field, Wavefront, AddWavefront, amplitude
export Piston, Tilt, Tip, Defocus, Astigmatism, Coma, AddAberration
export Trefoil, Spherical, VerticalAstigmatism, ObliqueAstigmatism
export Aberration, HorizontalComa, VerticalComa, ObliqueTrefoil, VerticalTrefoil
export phase, project, correct
export circ
export rawphase, unwrapphase, retrievephase

cartesian_to_polar(x,y) = √(x^2+y^2), atan(x,y)
polar_to_cartesian(r,θ) = r*cos(θ),r*sin(θ)

include("coordinates.jl")
include("aberrations.jl")
include("wavefronts.jl")
include("phasestep.jl")

end # module


text:
"""
CSV provides fast, flexible reader & writer for delimited text files in various formats.

Reading:
- `CSV.File` reads delimited data and returns a `CSV.File` object, which allows dot-access to columns and iterating rows.
- `CSV.read` is similar to `CSV.File` but used when the input will be passed directly to a sink function such as a `DataFrame`.
- `CSV.Rows` reads delimited data and returns a `CSV.Rows` object, which allows "streaming" the data by iterating and thereby has a lower memory footprint than `CSV.File`.
- `CSV.Chunks` allows processing extremely large files in "batches" or "chunks".
Writing:
- `CSV.write` writes a [Tables.jl interface input](https://github.com/JuliaData/Tables.jl) such as a `DataFrame` to a csv file or an in-memory IOBuffer.
- `CSV.RowWriter` creates an iterator that produces csv-formatted strings for each row in the input table.
Here is an example of reading a csv file and passing the input to a `DataFrame`:
```julia
using CSV, DataFrames
ExampleInputDF = CSV.read("ExampleInputFile.csv", DataFrame)
```
Here is an example of writing out a `DataFrame` to a csv file:
```julia
using CSV, DataFrames
ExampleOutputDF = DataFrame(rand(10,10), :auto)
CSV.write("ExampleOutputFile.csv", ExampleOutputDF)
```
"""
module CSV

if !isdefined(Base, :contains)
    contains(haystack, needle) = occursin(needle, haystack)
end

# stdlib
using Mmap, Dates, Unicode
# Parsers.jl is used for core type parsing from byte buffers
# and all other parsing options (quoted fields, delimiters, dateformats etc.)
using Parsers
# Tables.jl allows integration with all other table/data file formats
using Tables
# PooledArrays.jl is used for materializing pooled columns
using PooledArrays
# SentinelArrays.jl allow efficient conversion from Vector{Union{T, Missing}} to Vector{T}
# it also provides the MissingVector and ChainedVector array types
using SentinelArrays
# InlineStrings and WeakRefStrings provide the InlineString and PosLenString types for more gc-friendly string types
using InlineStrings, WeakRefStrings
export PosLenString, InlineString
export String1, String3, String7, String15, String31, String63, String127, String255
# CodecZlib is used for unzipping gzip files
using CodecZlib
# FilePathsBase allows more structured file path types
using FilePathsBase

struct Error <: Exception
    msg::String
end

Base.showerror(io::IO, e::Error) = println(io, e.msg)

# constants
const DEFAULT_STRINGTYPE = InlineString
const DEFAULT_POOL = 0.25
const DEFAULT_ROWS_TO_CHECK = 30
const DEFAULT_MAX_WARNINGS = 100
const DEFAULT_MAX_INLINE_STRING_LENGTH = 32
const TRUE_STRINGS = ["true", "True", "TRUE", "T", "1"]
const FALSE_STRINGS = ["false", "False", "FALSE", "F", "0"]
const StringCodeUnits = Base.CodeUnits{UInt8, String}
const ValidSources = Union{Vector{UInt8}, SubArray{UInt8, 1, Vector{UInt8}}, StringCodeUnits, IO, Cmd, AbstractString, AbstractPath}
const MAX_INPUT_SIZE = Int === Int64 ? 2^42 : typemax(Int32)
const EMPTY_INT_ARRAY = Int[]

include("keyworddocs.jl")
include("utils.jl")
include("detection.jl")
include("context.jl")
include("file.jl")
include("chunks.jl")
include("rows.jl")
include("write.jl")

"""
`CSV.read(source, sink::T; kwargs...)` => T

Read and parses a delimited file or files, materializing directly using the `sink` function. Allows avoiding excessive copies
of columns for certain sinks like `DataFrame`.

$KEYWORD_DOCS
"""
function read(source, sink=nothing; copycols::Bool=false, kwargs...)
    if sink === nothing
        throw(ArgumentError("provide a valid sink argument, like `using DataFrames; CSV.read(source, DataFrame)`"))
    end
    Tables.CopiedColumns(CSV.File(source; kwargs...)) |> sink
end

include("precompile.jl")
_precompile_()

function __init__()
    CSV.Context(IOBuffer(CSV.PRECOMPILE_DATA))
    # CSV.File(IOBuffer(CSV.PRECOMPILE_DATA))
    # foreach(row -> row, CSV.Rows(IOBuffer(PRECOMPILE_DATA)))
    # CSV.File(joinpath(dirname(pathof(CSV)), "..", "test", "testfiles", "promotions.csv"))
end

end # module


text:
# Zygote fill has issues with non-numbers
ZygoteRules.@adjoint function fill(x::T, dims...) where {T}
    return ZygoteRules.pullback(x, dims...) do x, dims...
        return reshape([x for i in 1:prod(dims)], dims)
    end
end

## Uniform ##

ZygoteRules.@adjoint function Distributions.Uniform(args...)
    return ZygoteRules.pullback(TuringUniform, args...)
end

## Product

# Tests with `Kolmogorov` seem to fail otherwise?!
ZygoteRules.@adjoint function Distributions._logpdf(d::Product, x::AbstractVector{<:Real})
    return ZygoteRules.pullback(d, x) do d, x
        sum(map(logpdf, d.v, x))
    end
end
ZygoteRules.@adjoint function Distributions._logpdf(
    d::FillVectorOfUnivariate,
    x::AbstractVector{<:Real},
)
    return ZygoteRules.pullback(d, x) do d, x
        _flat_logpdf(d.v.value, x)
    end
end

# Loglikelihood of multi- and matrixvariate distributions: multiple samples
# workaround for Zygote issues discussed in
# https://github.com/TuringLang/DistributionsAD.jl/pull/198
ZygoteRules.@adjoint function Distributions.loglikelihood(
    d::MultivariateDistribution, x::AbstractMatrix{<:Real}
)
    return ZygoteRules.pullback(d, x) do d, x
        return sum(xi -> Distributions._logpdf(d, xi), eachcol(x))
    end
end
ZygoteRules.@adjoint function Distributions.loglikelihood(
    d::MatrixDistribution, x::AbstractArray{<:Real,3}
)
    return ZygoteRules.pullback(d, x) do d, x
        return sum(xi -> Distributions._logpdf(d, xi), eachslice(x; dims=3))
    end
end

## Wishart ##

# Custom adjoint since Zygote can't differentiate through `@warn`
# TODO: Remove when fixed upstream in Distributions
ZygoteRules.@adjoint function Wishart(df::T, S::AbstractPDMat{T}, warn::Bool = true) where T<:Real
    function _Wishart(df::T, S::AbstractPDMat{T}, warn::Bool = true) where T
        df > 0 || throw(ArgumentError("df must be positive. got $(df)."))
        p = dim(S)
        rnk = p
        singular = df <= p - 1
        if singular
            isinteger(df) || throw(ArgumentError("singular df must be an integer. got $(df)."))
            rnk = convert(Integer, df)
            warn && _warn("got df <= dim - 1; returning a singular Wishart")
        end
        logc0 = Distributions.wishart_logc0(df, S, rnk)
        R = Base.promote_eltype(T, logc0)
        prom_S = convert(AbstractArray{T}, S)
        Wishart{R, typeof(prom_S), typeof(rnk)}(R(df), prom_S, R(logc0), rnk, singular)
    end
    return ZygoteRules.pullback(_Wishart, df, S, warn)
end

_warn(msg) = @warn(msg)
@non_differentiable _warn(msg)

ZygoteRules.@adjoint function Distributions.Wishart(df::Real, S::AbstractMatrix{<:Real})
    return ZygoteRules.pullback(TuringWishart, df, S)
end
ZygoteRules.@adjoint function Distributions.InverseWishart(
    df::Real,
    S::AbstractMatrix{<:Real}
)
    return ZygoteRules.pullback(TuringInverseWishart, df, S)
end

## General definitions of `logpdf` for arrays

ZygoteRules.@adjoint function Distributions.logpdf(
    dist::MultivariateDistribution,
    X::AbstractMatrix{<:Real},
)
    size(X, 1) == length(dist) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    return ZygoteRules.pullback(dist, X) do dist, X
        return map(xi -> Distributions._logpdf(dist, xi), eachcol(X))
    end
end

ZygoteRules.@adjoint function Distributions.logpdf(
    dist::MatrixDistribution,
    X::AbstractArray{<:Real,3},
)
    (size(X, 1), size(X, 2)) == size(dist) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    return ZygoteRules.pullback(dist, X) do dist, X
        return map(xi -> Distributions._logpdf(dist, xi), eachslice(X; dims=3))
    end
end

ZygoteRules.@adjoint function Distributions.logpdf(
    dist::MatrixDistribution,
    X::AbstractArray{<:AbstractMatrix{<:Real}},
)
    return ZygoteRules.pullback(dist, X) do dist, X
        return map(x -> logpdf(dist, x), X)
    end
end


text:
__precompile__(true)
module FWF

@static if VERSION < v"0.7.0-DEV.2005"
    using DataStreams, Missings, DataFrames
else
    using DataStreams, Missings, DataFrames, Dates, Mmap
end

struct ParsingException <: Exception
    msg::String
end

"""
Configuration Settings for fixed width file parsing.

  * `skip`        : integer of the number of lines to skip; default `0`
  * `trimstrings` : true if strings should be trimmed; default `true`
  * `usemissings` : true if fields should be checked for null values; default `true`
  * `errorlevel`  : if `:parse` then as much as possible is parsed and missing data is replaced by `missing`;
                    if `:skip` then malformed line is skipped on error;
                    if `:error` then an exception is thrown on error; default `:parse`
  * `countnybytes`: true if field parsing should happen by bytes, false for character based parsing; default `true`
  * `missingvals` : Set{String} for values that equal missing
  * `dateformats` : Dictionary in the form of Int=>DateFormat to specify date formats for a column
  * `columnrange` : Vector of UnitRanges that specifcy the widths of each column.
  * ``
"""

struct Options
    usemissings::Bool
    trimstrings::Bool
    errorlevel::Symbol
    unitbytes::Bool
    skip::Int
    missingvals::Set{String}
    dateformats::Dict{Int, DateFormat}
    columnrange::Vector{UnitRange{Int}}
end

function Options(;usemissings=true, trimstrings=true, errorlevel=:parse, unitbytes=true,
                 skip=0, missingvals=Set{String}(),
                 dateformats=Dict{Int, DateFormat}(), columnrange=Vector{UnitRange{Int}}())
    if !usemissings && (errorlevel == :parse)
        println(STDERR, "Warning: Combination of usemissings==false and errorlevel==:parse\n"*
               "will lead to an error when malformed lines are present in the data.\n"*
               "In order to avoid this set usemissings to true or errorlevel to :skip.")
    end
    if !(errorlevel in [:parse, :skip, :error])
        throw(ArgumentError("Allowed values for errorlevel are :parse, :skip or :error"))
    end
    Options(usemissings, trimstrings, errorlevel, unitbytes,
            skip, missingvals, dateformats, columnrange)
end

function Base.show(io::IO, op::Options)
    println(io, "   FWF.Options:")
    println(io, "     nullcheck: ", op.usemissings)
    println(io, "   trimstrings: ", op.trimstrings)
    println(io, "   errorlevel: ", op.errorlevel)
    println(io, "  countbybytes: ", op.unitbytes)
    println(io, "          skip: ", op.skip)
    println(io, "   missingvals:", )
    show(io, op.missingvals)
    println(io)
    println(io, "   dateformats:")
    show(io, op.dateformats)
    println(io)
    println(io, "   columnranges:")
    show(io, op.columnrange)
    println(io)
end

include("Source.jl")
include("parsefields.jl")
include("io.jl")
include("scan.jl")
end


text:
using Revise
using SoapySDR
using GLMakie
using LinearAlgebra
using FFTW
using DSP

# Don't forget to add/import a device-specific plugin package!
# using xtrx_jll
# using SoapyLMS7_jll
using SoapyRTLSDR_jll
# using SoapyPlutoSDR_jll
# using SoapyUHD_jll

# Here we want to test the behavior in a loop, to ensure
# that we can block on buffer overflow conditions, and
# handle partial reads, measure latnecy, etc

freq = 104.3e6u"Hz"

function makie_viz(freq=freq)
    dev = open(Devices()[1])
    rx_chan = dev.rx[1]
    @show rx_chan
    rx_chan.gain_mode = true
    rx_chan.frequency = freq
    rx_stream = SoapySDR.Stream(ComplexF32, [rx_chan])
    @show rx_stream.mtu
    buf = Vector{ComplexF32}(undef, rx_stream.mtu)
    @show typeof(fft)
    sum = zeros(Float32, length(buf))
    norms = Observable(Vector{Float32}(undef,length(buf)))

    println("Makie Setup...")
    fig = Figure(); display(fig)
    ax = Axis(fig[1,1])
    xlims!(ax, 0, 100)
    ylims!(ax, -500, 500)
    lines!(ax, norms)

    h = Float32.(DSP.hamming(length(buf)))

    println("Starting fft loop..")
    Base.atexit(()->SoapySDR.deactivate!(rx_stream))

    SoapySDR.activate!(rx_stream)
    k = 1
    while true
        read!(rx_stream, (buf,))
        buf = fftshift(FFTW.fft!(buf.*h))
        norms[] .= (norm.(buf))
        #norms[] .= sum ./ k
        norms[] = norms[]
        println("fft") # makie needs a slight delay here
        k += 1
    end

end


text:
## ExoplanetsSysSim/examples/tess/abc_run.jl
## (c) 2019 Danley C. Hsu & Eric B. Ford
# Script for producing TESS planet candidate occurrence rate estimates

include("abc_setup.jl")

#using SysSimABC
using ExoplanetsSysSim
using JLD
using StatsBase

out2txt = false # Write occurrence rates & densities to text files
expandpart = false # Expand final generation for robust posteriors
prior_choice = "uniform"
bin_size_factor = 2.0

println("Setting up simulation...")
@time abc_plan = SysSimABC.setup_abc(prior_choice = prior_choice, bin_size_factor = bin_size_factor)
println("")
println("Running simulation...")
@time output = SysSimABC.run_abc(abc_plan)
println("")
println("Running simulation (part 2)...")
@time abc_plan = SysSimABC.setup_abc_p2(abc_plan)
@time output = SysSimABC.run_abc(abc_plan, output)
#@time abc_plan = change_distance()
#@time output = run_abc(abc_plan, output)
println("")

save(string("test-pop-out.jld"), "output", output, "ss_true", EvalSysSimModel.get_ss_obs())

if expandpart
    println("Expanding to large generation...")
    @time theta_largegen, weights_largegen = SysSimABC.run_abc_largegen(abc_plan, output, EvalSysSimModel.get_ss_obs(), output.accept_log.epsilon[end-1], npart=200)
    println("")

    save(string("test-pop-out.jld"), "output", output, "ss_true", EvalSysSimModel.get_ss_obs(), "theta_largegen", theta_largegen, "weights_largegen", weights_largegen)
end

if out2txt
    file_rate = open("rate_output.txt", "w")
    file_dens = open("dens_output.txt", "w")
end

limitP = get_any(EvalSysSimModel.sim_param_closure, "p_lim_arr", Array{Float64,1})
limitR = get_any(EvalSysSimModel.sim_param_closure, "r_lim_arr", Array{Float64,1})
const r_dim = length(limitR)-1

if expandpart
    weight_vec = aweights(weights_largegen)
    #weight_vec = aweights(fill(1.0, length(weights_largegen)))
else
    weight_vec = aweights(output.weights)
    #weight_vec = aweights(fill(1.0, length(output.weights)))
end

for p_ind = 1:(length(limitP)-1)
    col_ind = (p_ind-1)*(r_dim)+1
    for r_ind = 1:r_dim
        bin_ind = (p_ind-1)*(r_dim)+r_ind
        dens_denom = 1.0/log(limitP[p_ind+1]/limitP[p_ind])/log(limitR[r_ind+1]/limitR[r_ind])

        if prior_choice == "dirichlet" && r_dim > 1
            col_ind = (p_ind-1)*(r_dim+1)+1
            bin_ind = (p_ind-1)*(r_dim+1)+r_ind+1
            if expandpart
                quant_arr = quantile(theta_largegen[bin_ind,:].*theta_largegen[col_ind,:], weight_vec, [0.1587, 0.5, 0.8413])
            else
                quant_arr = quantile(output.theta[bin_ind,:].*output.theta[col_ind,:], weight_vec, [0.1587, 0.5, 0.8413])
            end
        elseif prior_choice == "beta"
            col_lambda = bin_size_factor * 3 * log(limitP[p_ind+1]/limitP[p_ind])/log(2)
            if expandpart
                quant_arr = quantile(theta_largegen[bin_ind,:]*col_lambda, weight_vec, [0.1587, 0.5, 0.8413])
            else
                quant_arr = quantile(output.theta[bin_ind,:]*col_lambda, weight_vec, [0.1587, 0.5, 0.8413])
            end
        else
            if expandpart
                quant_arr = quantile(theta_largegen[bin_ind,:], weight_vec, [0.1587, 0.5, 0.8413])
            else
                quant_arr = quantile(output.theta[bin_ind,:], weight_vec, [0.1587, 0.5, 0.8413])
            end
        end

        println("-----------------------------")
        println("Orbital Period (day) = ", string(limitP[p_ind:p_ind+1]), " / Planet Radius (R_earth) = ", string(limitR[r_ind:r_ind+1]/ExoplanetsSysSim.earth_radius))
        println("")
        println("Rate = ", string(quant_arr[2], " + ", quant_arr[3]-quant_arr[2], " - ", quant_arr[2]-quant_arr[1]))
        println("Density = ", string(quant_arr[2]*dens_denom, " + ", (quant_arr[3]-quant_arr[2])*dens_denom, " - ", (quant_arr[2]-quant_arr[1])*dens_denom))

        if out2txt
            write(file_rate, string(quant_arr[2], " + ", quant_arr[3]-quant_arr[2], " - ", quant_arr[2]-quant_arr[1], "\n"))
            write(file_dens, string(quant_arr[2]*dens_denom, " + ", (quant_arr[3]-quant_arr[2])*dens_denom, " - ", (quant_arr[2]-quant_arr[1])*dens_denom, "\n"))
        end
    end
end

if out2txt
    close(file_rate)
    close(file_dens)
end

println("-----------------------------")
println("")
println(EvalSysSimModel.get_ss_obs())


text:
pkgs = ["OhMyREPL", "Revise", "Pluto", "JET"]

using Pkg
Pkg.add(pkgs)

# pkg"add LanguageServer#master"
# pkg"add StaticLint#master"
# pkg"add SymbolServer#master"

max_stars_repo_name: mgarbiso/markus-config

max_issues_repo_name: mgarbiso/markus-config

max_forks_repo_name: mgarbiso/markus-config


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

mutable struct DataLakeStoreAccountPropertiesBasic <: SwaggerModel
    accountId::Any # spec type: Union{ Nothing, String } # spec name: accountId
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    state::Any # spec type: Union{ Nothing, String } # spec name: state
    creationTime::Any # spec type: Union{ Nothing, DateTime } # spec name: creationTime
    lastModifiedTime::Any # spec type: Union{ Nothing, DateTime } # spec name: lastModifiedTime
    endpoint::Any # spec type: Union{ Nothing, String } # spec name: endpoint

    function DataLakeStoreAccountPropertiesBasic(;accountId=nothing, provisioningState=nothing, state=nothing, creationTime=nothing, lastModifiedTime=nothing, endpoint=nothing)
        o = new()
        validate_property(DataLakeStoreAccountPropertiesBasic, Symbol("accountId"), accountId)
        setfield!(o, Symbol("accountId"), accountId)
        validate_property(DataLakeStoreAccountPropertiesBasic, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(DataLakeStoreAccountPropertiesBasic, Symbol("state"), state)
        setfield!(o, Symbol("state"), state)
        validate_property(DataLakeStoreAccountPropertiesBasic, Symbol("creationTime"), creationTime)
        setfield!(o, Symbol("creationTime"), creationTime)
        validate_property(DataLakeStoreAccountPropertiesBasic, Symbol("lastModifiedTime"), lastModifiedTime)
        setfield!(o, Symbol("lastModifiedTime"), lastModifiedTime)
        validate_property(DataLakeStoreAccountPropertiesBasic, Symbol("endpoint"), endpoint)
        setfield!(o, Symbol("endpoint"), endpoint)
        o
    end
end # type DataLakeStoreAccountPropertiesBasic

const _property_map_DataLakeStoreAccountPropertiesBasic = Dict{Symbol,Symbol}(Symbol("accountId")=>Symbol("accountId"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("state")=>Symbol("state"), Symbol("creationTime")=>Symbol("creationTime"), Symbol("lastModifiedTime")=>Symbol("lastModifiedTime"), Symbol("endpoint")=>Symbol("endpoint"))
const _property_types_DataLakeStoreAccountPropertiesBasic = Dict{Symbol,String}(Symbol("accountId")=>"String", Symbol("provisioningState")=>"String", Symbol("state")=>"String", Symbol("creationTime")=>"DateTime", Symbol("lastModifiedTime")=>"DateTime", Symbol("endpoint")=>"String")
Base.propertynames(::Type{ DataLakeStoreAccountPropertiesBasic }) = collect(keys(_property_map_DataLakeStoreAccountPropertiesBasic))
Swagger.property_type(::Type{ DataLakeStoreAccountPropertiesBasic }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DataLakeStoreAccountPropertiesBasic[name]))}
Swagger.field_name(::Type{ DataLakeStoreAccountPropertiesBasic }, property_name::Symbol) =  _property_map_DataLakeStoreAccountPropertiesBasic[property_name]

const _allowed_DataLakeStoreAccountPropertiesBasic_provisioningState = ["Failed", "Creating", "Running", "Succeeded", "Patching", "Suspending", "Resuming", "Deleting", "Deleted", "Undeleting", "Canceled"]

const _allowed_DataLakeStoreAccountPropertiesBasic_state = ["Active", "Suspended"]

function check_required(o::DataLakeStoreAccountPropertiesBasic)
    true
end

function validate_property(::Type{ DataLakeStoreAccountPropertiesBasic }, name::Symbol, val)
    if name === Symbol("provisioningState")
        Swagger.validate_param(name, "DataLakeStoreAccountPropertiesBasic", :enum, val, _allowed_DataLakeStoreAccountPropertiesBasic_provisioningState)
    end
    if name === Symbol("state")
        Swagger.validate_param(name, "DataLakeStoreAccountPropertiesBasic", :enum, val, _allowed_DataLakeStoreAccountPropertiesBasic_state)
    end
end


text:
module StatsMakie

using Observables
using AbstractPlotting
import AbstractPlotting: convert_arguments, used_attributes, plot!, combine, to_plotspec
using AbstractPlotting: plottype, Plot, PlotFunc, to_tuple
using AbstractPlotting: node_pairs, extrema_nan, automatic, default_theme
using AbstractPlotting: GeometryTypes

# Moved in https://github.com/JuliaGizmos/Observables.jl/pull/40
if isdefined(Observables, :to_value)
    using Observables: to_value
else
    using AbstractPlotting: to_value
end

using Statistics, KernelDensity
import StatsBase
using Distributions
using IntervalSets
using Tables, StructArrays
using StructArrays: uniquesorted, finduniquesorted, GroupPerm
using IntervalSets: Interval, endpoints
using Loess

export Data, Group, Style
export Position
export bycolumn
export frequency

include(joinpath("group", "analysis.jl"))
include(joinpath("group", "scales.jl"))
include(joinpath("group", "group.jl"))
include(joinpath("group", "tables.jl"))
include(joinpath("group", "dodge.jl"))

include(joinpath("typerecipes", "density.jl"))
include(joinpath("typerecipes", "histogram.jl"))
include(joinpath("typerecipes", "distribution.jl"))
include(joinpath("typerecipes", "smooth.jl"))
include(joinpath("typerecipes", "namedarray.jl"))

# include(joinpath("recipes", "corrplot.jl"))
include(joinpath("recipes", "boxplot.jl"))
include(joinpath("recipes", "violin.jl"))
include(joinpath("recipes", "ribbon.jl"))
include(joinpath("recipes", "errorbar.jl"))

end