text:
module DSP

importall Base.FFTW
import Base.FFTW.normalization

export FFTW, filt, deconv, conv, conv2, xcorr, fftshift, ifftshift,
       dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!,
       # the rest are defined imported from FFTW:
       fft, bfft, ifft, rfft, brfft, irfft,
       plan_fft, plan_bfft, plan_ifft, plan_rfft, plan_brfft, plan_irfft,
       fft!, bfft!, ifft!, plan_fft!, plan_bfft!, plan_ifft!

function filt{T<:Number}(b::Union(AbstractVector{T}, T),a::Union(AbstractVector{T}, T),x::AbstractVector{T})
    if isempty(b); error("b must be non-empty"); end
    if isempty(a); error("a must be non-empty"); end
    if a[1]==0; error("a[1] must be nonzero"); end

    as = length(a)
    bs = length(b)
    sz = max(as, bs)

    if sz == 1
        return (b[1]/a[1]).*x
    end

    if bs<sz
        newb = zeros(T,sz)
        newb[1:bs] = b
        b = newb
    end

    xs = size(x,1)
    y = Array(T, xs)
    silen = sz-1
    si = zeros(T, silen)

    if a[1] != 1
        norml = a[1]
        a ./= norml
        b ./= norml
    end

    if as > 1
        if as<sz
            newa = zeros(T,sz)
            newa[1:as] = a
            a = newa
        end

        @inbounds begin
            for i=1:xs
                y[i] = si[1] + b[1]*x[i]
                for j=1:(silen-1)
                    si[j] = si[j+1] + b[j+1]*x[i] - a[j+1]*y[i]
                end
                si[silen] = b[silen+1]*x[i] - a[silen+1]*y[i]
            end
        end
    else
        @inbounds begin
            for i=1:xs
                y[i] = si[1] + b[1]*x[i]
                for j=1:(silen-1)
                    si[j] = si[j+1] + b[j+1]*x[i]
                end
                si[silen] = b[silen+1]*x[i]
            end
        end
    end
    return y
end

function deconv{T}(b::StridedVector{T}, a::StridedVector{T})
    lb = size(b,1)
    la = size(a,1)
    if lb < la
        return [zero(T)]
    end
    lx = lb-la+1
    x = zeros(T, lx)
    x[1] = 1
    filt(b, a, x)
end

function conv{T<:Base.LinAlg.BlasFloat}(u::StridedVector{T}, v::StridedVector{T})
    nu = length(u)
    nv = length(v)
    n = nu + nv - 1
    np2 = n > 1024 ? nextprod([2,3,5], n) : nextpow2(n)
    upad = [u, zeros(T, np2 - nu)]
    vpad = [v, zeros(T, np2 - nv)]
    if T <: Real
        p = plan_rfft(upad)
        y = irfft(p(upad).*p(vpad), np2)
    else
        p = plan_fft!(upad)
        y = ifft!(p(upad).*p(vpad))
    end
    return y[1:n]
end
conv{T<:Integer}(u::StridedVector{T}, v::StridedVector{T}) = int(conv(float(u), float(v)))
conv{T<:Integer, S<:Base.LinAlg.BlasFloat}(u::StridedVector{T}, v::StridedVector{S}) = conv(float(u), v)
conv{T<:Integer, S<:Base.LinAlg.BlasFloat}(u::StridedVector{S}, v::StridedVector{T}) = conv(u, float(v))

function conv2{T}(u::StridedVector{T}, v::StridedVector{T}, A::StridedMatrix{T})
    m = length(u)+size(A,1)-1
    n = length(v)+size(A,2)-1
    B = zeros(T, m, n)
    B[1:size(A,1),1:size(A,2)] = A
    u = fft([u;zeros(T,m-length(u))])
    v = fft([v;zeros(T,n-length(v))])
    C = ifft(fft(B) .* (u * v.'))
    if T <: Real
        return real(C)
    end
    return C
end

function conv2{T}(A::StridedMatrix{T}, B::StridedMatrix{T})
    sa, sb = size(A), size(B)
    At = zeros(T, sa[1]+sb[1]-1, sa[2]+sb[2]-1)
    Bt = zeros(T, sa[1]+sb[1]-1, sa[2]+sb[2]-1)
    At[1:sa[1], 1:sa[2]] = A
    Bt[1:sb[1], 1:sb[2]] = B
    p = plan_fft(At)
    C = ifft(p(At).*p(Bt))
    if T <: Real
        return real(C)
    end
    return C
end
conv2{T<:Integer}(A::StridedMatrix{T}, B::StridedMatrix{T}) = int(conv2(float(A), float(B)))
conv2{T<:Integer}(u::StridedVector{T}, v::StridedVector{T}, A::StridedMatrix{T}) = int(conv2(float(u), float(v), float(A)))

function xcorr(u, v)
    su = size(u,1); sv = size(v,1)
    if su < sv
        u = [u;zeros(eltype(u),sv-su)]
    elseif sv < su
        v = [v;zeros(eltype(v),su-sv)]
    end
    flipud(conv(flipud(u), v))
end

fftshift(x) = circshift(x, div([size(x)...],2))

function fftshift(x,dim)
    s = zeros(Int,ndims(x))
    s[dim] = div(size(x,dim),2)
    circshift(x, s)
end

ifftshift(x) = circshift(x, div([size(x)...],-2))

function ifftshift(x,dim)
    s = zeros(Int,ndims(x))
    s[dim] = -div(size(x,dim),2)
    circshift(x, s)
end

# Discrete cosine and sine transforms via FFTW's r2r transforms;
# we follow the Matlab convention and adopt a unitary normalization here.
# Unlike Matlab we compute the multidimensional transform by default,
# similar to the Julia fft functions.

fftwcopy{T<:fftwNumber}(X::StridedArray{T}) = copy(X)
fftwcopy{T<:Real}(X::StridedArray{T}) = float(X)
fftwcopy{T<:Complex}(X::StridedArray{T}) = complex128(X)

for (f, fr2r, Y, Tx) in ((:dct, :r2r, :Y, :Number),
                         (:dct!, :r2r!, :X, :fftwNumber))
    plan_f = symbol(string("plan_",f))
    plan_fr2r = symbol(string("plan_",fr2r))
    fi = symbol(string("i",f))
    plan_fi = symbol(string("plan_",fi))
    Ycopy = Y == :X ? 0 : :(Y = fftwcopy(X))
    @eval begin
        function $f{T<:$Tx}(X::StridedArray{T}, region)
            $Y = $fr2r(X, REDFT10, region)
            scale!($Y, sqrt(0.5^length(region) * normalization(X,region)))
            sqrthalf = sqrt(0.5)
            r = map(n -> 1:n, [size(X)...])
            for d in region
                r[d] = 1:1
                $Y[r...] *= sqrthalf
                r[d] = 1:size(X,d)
            end
            return $Y
        end

        function $plan_f{T<:$Tx}(X::StridedArray{T}, region,
                                 flags::Unsigned, timelimit::Real)
            p = $plan_fr2r(X, REDFT10, region, flags, timelimit)
            sqrthalf = sqrt(0.5)
            r = map(n -> 1:n, [size(X)...])
            nrm = sqrt(0.5^length(region) * normalization(X,region))
            return X::StridedArray{T} -> begin
                $Y = p(X)
                scale!($Y, nrm)
                for d in region
                    r[d] = 1:1
                    $Y[r...] *= sqrthalf
                    r[d] = 1:size(X,d)
                end
                return $Y
            end
        end

        function $fi{T<:$Tx}(X::StridedArray{T}, region)
            $Ycopy
            scale!($Y, sqrt(0.5^length(region) * normalization(X, region)))
            sqrt2 = sqrt(2)
            r = map(n -> 1:n, [size(X)...])
            for d in region
                r[d] = 1:1
                $Y[r...] *= sqrt2
                r[d] = 1:size(X,d)
            end
            return r2r!($Y, REDFT01, region)
        end

        function $plan_fi{T<:$Tx}(X::StridedArray{T}, region,
                                 flags::Unsigned, timelimit::Real)
            p = $plan_fr2r(X, REDFT01, region, flags, timelimit)
            sqrt2 = sqrt(2)
            r = map(n -> 1:n, [size(X)...])
            nrm = sqrt(0.5^length(region) * normalization(X,region))
            return X::StridedArray{T} -> begin
                $Ycopy
                scale!($Y, nrm)
                for d in region
                    r[d] = 1:1
                    $Y[r...] *= sqrt2
                    r[d] = 1:size(X,d)
                end
                return p($Y)
            end
        end

    end
    for (g,plan_g) in ((f,plan_f), (fi, plan_fi))
        @eval begin
            $g{T<:$Tx}(X::StridedArray{T}) = $g(X, 1:ndims(X))

            $plan_g(X, region, flags::Unsigned) =
              $plan_g(X, region, flags, NO_TIMELIMIT)
            $plan_g(X, region) =
              $plan_g(X, region, ESTIMATE, NO_TIMELIMIT)
            $plan_g{T<:$Tx}(X::StridedArray{T}) =
              $plan_g(X, 1:ndims(X), ESTIMATE, NO_TIMELIMIT)
        end
    end
end

# DCT of scalar is just the identity:
dct(x::Number, dims) = length(dims) == 0 || dims[1] == 1 ? x : throw(BoundsError())
idct(x::Number, dims) = dct(x, dims)
dct(x::Number) = x
idct(x::Number) = x
plan_dct(x::Number, dims, flags, tlim) = length(dims) == 0 || dims[1] == 1 ? y::Number -> y : throw(BoundsError())
plan_idct(x::Number, dims, flags, tlim) = plan_dct(x, dims)
plan_dct(x::Number) = y::Number -> y
plan_idct(x::Number) = y::Number -> y

end # module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# ---
# title: Amplitude swapping
# id: demo_swapphase
# cover: assets/swapphase.gif
# author: Peter Kovesi
# date: 2018-10-26
# ---

# A demonstration of the importance of phase information in images. Given two
# images`swapphase()` takes their Fourier transforms and constructs two new, synthetic,
# images formed from the swapped phase and amplitude imformation.  In general it is the
# phase information that dominates.  However, for textures where the amplitude spectra can
# be concentrated in a limited set of locations, the reverse can apply.

# See [Oppenheim and Lim's paper "The importance of phase in signals". Proceedings of the
# IEEE. Volume: 69 , Issue: 5 , May 1981](https://ieeexplore.ieee.org/document/1456290)

using TestImages
using Images
using ImagePhaseCongruency

img1 = centered(Float64.(Gray.(restrict(testimage("lighthouse")))))[-127:128, -127:128]
img2 = restrict(Float64.(testimage("mandril_gray")))[1:256, 1:256]

(newimg1, newimg2) = swapphase(img1, img2)

mosaic(Gray.(img1), newimg1, img2, newimg2; nrow=2)
# Bottom 1) left: phase of lighthouse, amplitude of Mandrill 2) right: amplitude of lighthouse, phase of Mandrill

# save cover image #src
isdir("assets") || mkdir("assets") #src
save(joinpath("assets", "swapphase.gif"), Images.gif([img1, newimg1]); fps=1) #src

ext: jl

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
# Types to do a HF calculation
abstract type BasisFunction end

#= 1s Gaussian-type function
Smallest possible basis set is called the minimal basis set, and it contais one orbital.
For example, hydrogen has just one orbital, but carbon has 5 (1s, 2s, 2px, 2py, and 2pz) even though one of the p orbitals for carbon atom will be unoccupied.
=#
struct Gaussian1s <: BasisFunction
  #Gaussian orbital exponent
  alpha :: Float64
  #nuclear coordinates
  center::Float64
end

# Slater Type Orbital fit with N gaussians (STO-NG) type basis
struct STONG <: BasisFunction
  n::Int64
  #contraction coeffiecents
  d::Array{Float64}
  #primative gaussians
  g::Array{Gaussian1s}
end

#=
The STO-3G basis is a very well-known minimal basis set which contracts 3 Gaussian functions to approximate the more accurate (but more difficult to compute) Slater type orbitals. Although a contracted GTO might give a good approximation to an atomic orbital, it lacks any flexibility to expand or shrink in the presence of other atoms in a molecule. Hence, a minimal basis set such as STO-3G is not capable of giving highly accurate results.
=#


text:
struct IAmaSingleton end

x = IAmaSingleton()
y = IAmaSingleton()

println("x == y is $(x == y) and x === y is $(x === y).")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Pkg

Pkg.activate(".")

using BeerTest


text:
if VERSION >= v"0.7.0-"
    using Pkg
    Pkg.status()
else
    for line in split(strip(readstring("REQUIRE")), '\n')[2:end]
        name = split(line)[1]
        println(name, "\t", Pkg.installed(name))
    end
end


text:
module ARCTR

# stdlib
using LinearAlgebra, Logging, Printf, SparseArrays
# JSO
using HSL, Krylov, LinearOperators, NLPModels, SolverCore, SolverTools
# Stopping
using Stopping, StoppingInterface

# Selective includes.
include("hessian_rep.jl")

include("./utils/show.jl")
include("./utils/pdata_struct.jl")
include("./utils/utils.jl")
include("./utils/increase_decrease.jl")
include("./utils/ldlt_symm.jl")

path = joinpath(dirname(@__FILE__), "SolveModel")
files = filter(x -> x[(end-2):end] == ".jl", readdir(path))
for file in files
    if file in ["krylov_aux.jl"]
        continue
    end
    include("SolveModel/" * file)
end

pathHO = joinpath(dirname(@__FILE__), "SolveModel", "high-order-correction")
files = filter(x -> x[(end-2):end] == ".jl", readdir(pathHO))
for file in files
    if file in []
        continue
    end
    include("SolveModel/high-order-correction/" * file)
end

path = joinpath(dirname(@__FILE__), "PreProcess")
files = filter(x -> x[(end-2):end] == ".jl", readdir(path))
for file in files
    if occursin(r"MA", file)
        continue # Remove MA57 and MA97 solvers for now
    end
    include("PreProcess/" * file)
end

####################################################################
## Model, temporary, shoudln't be here
# include("autodiff_high_order_model.jl")
####################################################################

include("TRARC.jl")

export ALL_solvers

ALL_solvers = Symbol[]

include("solvers.jl")

for fun in keys(solvers_const)
    push!(ALL_solvers, fun)

    ht, pt, sm, ka = ARCTR.solvers_const[fun]
    @eval begin
        function $fun(nlpstop::NLPStopping; kwargs...)
            kw_list = Dict{Symbol,Any}()
            if $ka != ()
                for t in $ka
                    push!(kw_list, t)
                end
            end
            merge!(kw_list, Dict(kwargs))
            TRARC(nlpstop; hess_type = $ht, pdata_type = $pt, solve_model = $sm, kw_list...)
        end
    end
    @eval begin
        function $fun(nlp::AbstractNLPModel; kwargs...)
            # https://github.com/vepiteski/Stopping.jl/blob/4036e996be4e9521bb40725a659a03d23fd70e03/src/Stopping/nlp_admissible_functions.jl#L12
            nlpstop = NLPStopping(
                nlp;
                optimality_check = (pb, state) ->
                    unconstrained_check(pb, state, pnorm = 2.0),
                kwargs...,
            )
            nlpstop = $fun(nlpstop; kwargs...)
            return stopping_to_stats(nlpstop)
        end
    end
    @eval export $fun
end

end # module


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

using Test, Libdl, nghttp2_jll

struct nghttp2_info
    age::Cint
    version_num::Cint
    version_str::Cstring
    proto_str::Cstring
end

@testset "nghttp2_jll" begin
    info = unsafe_load(ccall((:nghttp2_version,libnghttp2), Ptr{nghttp2_info}, (Cint,), 0))
    @test VersionNumber(unsafe_string(info.version_str)) == v"1.41.0"
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
################################################################################

using HiddenMarkovModelReaders
using Test
using DelimitedFiles

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

tests = [

  # private functions
  "reset.jl",
  "feed.jl",
  "backTrace.jl",
  "sort.jl",

  # export functions
  "distance.jl",

  "setup.jl",

  "process.jl",

]

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

@testset verbose = true "Hidden Markov Model Readers" begin

  for τ ∈ tests
    include(τ)
  end

end

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

max_stars_repo_name: DanielRivasMD/HiddenMarkovModelReaders

max_issues_repo_name: DanielRivasMD/HiddenMarkovModelReaders


text:
# Simulation of Linear Model Data

# Starting Simrel Code
function simrel{T}(
    n::T         = 100,
    p::T         = 10,
    q::T         = 5,
    relpos::T    = [1, 2, 3],
    gamma::T     = 0.2,
    R2::T        = 0.9;
    ntest::T     = nothing,
    muY::T       = nothing,
    muX::T       = nothing,
    lambdaMin::T = 10e-4,
    sim::T       = nothing
)
    m = length(relpos) ## Number of relevant components
    if (q < m)
        error("Number of relevant predictors must at least equal to $(m).")
    end
    irrelPosition = setdiff(collect(1:p), relpos)
    extraDim = q - m
    predPosition = vcat(relpos, sample(irrelPosition, extraDim,
                                    replace = false, ordered = true))
    nu = lambdaMin * exp(-gamma)/(1 - lambdaMin)
    if (lambdaMin < 0 || lambdaMin >= 1)
        error("Parameter lambdaMin must be in between 0 and 1")
    end
    lambdas = (exp(-gamma * vec(1:p)) + nu) / (exp(-gamma) + nu)
    sigmaZ = diagm(lambdas)
    sigmaZinv = inv(sigmaZ)
    sigmaZy = zeros(p, 1)
    uniformSample = rand(Uniform(-1, 1), m)
    covarianceSample = sign(uniformSample) .* sqrt(R2 * abs(uniformSample) /
             sum(abs(uniformSample)) .* lambdas[relpos])
    sigmaZy[relpos] = covarianceSample
    sigmaY = 1
    sigma = vcat(hcat(sigmaY, sigmaZy'), hcat(sigmaZy, sigmaZ))

    # Rotation Matrix

    Q = randn(q, q) |> x -> x .- mean(x, 1) # Centering columns
    Rq = qr(Q)[:1]
    R = eye(p)
    R[predPosition, predPosition] = Rq
    if (q < (p - 1))
      Q = randn(p - q, p - q) |> x -> x .- mean(x, 1)
      Rnq = qr(Q)[:1]
      notPredPos = setdiff(1:p, predPosition)
      R[notPredPos, notPredPos] = Rnq
    end

    # Regression coefficients
    betaZ = sigmaZinv * sigmaZy
    betaX = R * betaZ
    beta0 = 0

    if !(muY == nothing)
      beta0 = beta0 + muY
    end
    if !(muX == nothing)
      beta0 = beta0 - transpose(betaX) * muX
    end

    R2 = transpose(sigmaZy) * betaZ
    minerror = sigmaY - R2

    if all(eig(sigma)[:1] .> 0)
      sigmaRot = chol(sigma)
      ucal = randn(n, p + 1)
      ucal_1 = ucal * sigmaRot
      Y = ucal_1[:, 1]
      if !(muY == nothing)
        Y = Y .+ muY
      end
      Z = ucal_1[:, 2:end]
      X = Z * transpose(R)
      if !(muX == nothing)
        X = X .+ muX
      end

      # Test Observations
      if !(ntest == nothing)
        utest = randn(ntest, p + 1)
        utest_1 = utest * sigmaRot
        testY = utest_1[:, 1]
        if !(muY == nothing)
          testY = testY .+ muY
          testZ = utest_1[:, 2:end]
          testX = testZ * transpose(R)
        end
        if !(muX == nothing)
          testX = testX .+ muX
        end
      else
        testX = nothing
        testY = nothing
      end
    else
      return error("Correlation matrix is not positive definite")
    end

    return @NT(
      beta0    = beta0,
      beta     = betaX,
      sigma    = sigma,
      R2       = R2,
      minerror = minerror,
      X        = X,
      Y        = Y,
      testX    = testX,
      testY    = testY
    )
end

max_stars_repo_name: ExtraSpace/JuliaSimulatr

max_issues_repo_name: ExtraSpace/JuliaSimulatr

max_forks_repo_name: ExtraSpace/JuliaSimulatr


text:
"""
    tuple_to_size(::Type{S}) where S<:Tuple
Converts a tuple given by `(N, M, ...)` into `Tuple{N, M, ...}`.
"""
Base.@pure tuple_to_size(S::Tuple{Vararg{Int}}) = Tuple{S...}

# Copied (with modifications) from StaticArrays.jl
@generated function check_shape_params(
    ::Type{S},
    ::Type{T},
    ::Type{Val{N}},
    ::Type{Val{L}},
) where {S,T,N,L}
    if !all(x -> isa(x, Int), S.parameters)
        return :(throw(ArgumentError("Shapes parameter S must be a tuple of Ints (e.g. `Tuple{3,3}`")))
    end

    if L != tuple_prod(S) || L < 0 || tuple_minimum(S) < 0 || tuple_length(S) != N
        return :(throw(ArgumentError("Size mismatch in Shapes parameters. Got size $S, dimension $N and length $L.")))
    end

    return nothing
end

@generated function promote_namedtuple_eltype(::Union{NT,Type{NT}}) where {NT<:NamedTuple}
    #eltypes = map(eltype, shapes.parameters)
    eltypes = Tuple(eltype(shape) for shape in NT.parameters[2].parameters)
    T = promote_tuple_eltype(Tuple{eltypes...})
    return quote
        @_inline_meta
        $T
    end
end

# Copied from ValueShapes.jl
"""
    Shapes.default_datatype(T::Type)
Return a default specific type U that is more specific than T, with U <: T.
e.g.
    Shapes.default_datatype(Real) === Float64
    Shapes.default_datatype(Integer) === Int
"""
function default_datatype end

@inline default_datatype(T::Type) = throw(ArgumentError("Type must be <: Real. Got $T."))
@inline default_datatype(T::Type{<:Real}) = _default_datatype(T)
@inline _default_datatype(::Type{>:Int}) = Int
@inline _default_datatype(::Type{>:Float64}) = Float64
@inline _default_datatype(::Type{>:Real}) = Float64
@inline _default_datatype(T::Type) = T

has_unit_axes(A) = all(ax->ax isa AbstractUnitRange{Int}, axes(A))
function check_has_unit_axes(A)
    has_unit_axes(A) || throw(ArgumentError("The axes of data must be <: AbstractUnitRange{Int}"))
end


text:
#######################################################################################
# GPU-crutch solve method
# Makes the simplest possible method for GPU-compatibility
# Out of place only
#######################################################################################
struct GPUSimpleRK4 <: DiffEqBase.AbstractODEAlgorithm end
export GPUSimpleRK4

@muladd function DiffEqBase.solve(prob::ODEProblem,
                          alg::GPUSimpleRK4;
                          dt = error("dt is required for this algorithm"))
  @assert !isinplace(prob)
  u0 = prob.u0
  tspan = prob.tspan
  f = prob.f
  p = prob.p
  t = tspan[1]
  tf = prob.tspan[2]
  ts = tspan[1]:dt:tspan[2]
  us = MVector{length(ts),typeof(u0)}(undef)
  us[1] = u0
  u = u0
  half = convert(eltype(u0),1//2)
  sixth = convert(eltype(u0),1//6)

  for i in 2:length(ts)
      uprev = u; t = ts[i]
      k1 = f(u,p,t)
      tmp = uprev+dt*half*k1
      k2 = f(tmp, p, t+half*dt)
      tmp = uprev+dt*half*k2
      k3 = f(tmp, p, t+half*dt)
      tmp = uprev+dt*k3
      k4 = f(tmp, p, t+dt)
      u = uprev+dt*sixth*(k1+2k2+2k3+k4)
      us[i] = u
  end

  sol = DiffEqBase.build_solution(prob,alg,ts,SArray(us),
                                  k = nothing, destats = nothing,
                                  calculate_error = false)
  DiffEqBase.has_analytic(prob.f) && DiffEqBase.calculate_solution_errors!(sol;timeseries_errors=true,dense_errors=false)
  sol
end


text:
# This file was generated, do not modify it. # hide
ℓλ = symlogdensity(post; noinline=(:λ,))

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
########### Defines the parameters & plot for Fig 3 a,c,e,g,i ###########
c = 0.5; # set the host contact rate as a constant 0.5
# import required variables/modules common to all simulations
include("../Model_scripts/fig3params.jl");
include("../Model_scripts/Functions.jl");
# now that we have all the requirements, let's run the simulation
run_simulation();
Np = N'
########### Figure 3a ###########
# plot the number of infected/uninfected hosts
labels = vcat("Infected", ["" for i in 1:1:n_parasites]);
plot_population_numbers(
    hcat(labels...),
    "c=0.5\n \nNumber of infected and uninfected hosts",
    100,
    size(Np)[1],
    "$figure_directory/graph_3a.png"
)
########### Figure 3c ###########
# calculate the weighted average mortality
ui_w_avg = calculate_mortality()
# calculate the weighted average β
βy_w_avg = calculate_horizontal_transmission()
# calculate weighted H0
H0_w = calculate_H0()
# calculate weighted V0
V0_w = calculate_V0()
# plot the average weighted R0
plot_reproductive_rate(5, size(Np)[1], "$figure_directory/graph_3c.png")
########### Figure 3e ###########
# plot the average V0
plot_vertical_reproductive_ratio(1, size(Np)[1], "$figure_directory/graph_3e.png")
########### Figure 3g ###########
# calculate the weighted virulence (for noise reduction)
vir_w_avg = calculate_average_virulence()
# plot the average horizontal transmission & weighted virulence
plot_horizontal_and_virulence(1, size(Np)[1], "$figure_directory/graph_3g.png")
########### Figure 3i ###########
# calculate & plot the population evenness through time
plot_evenness(1, size(Np)[1], "$figure_directory/graph_3i.png")


text:
let
 tot2 = 0
 tot3 = 0
 open("d2data.txt") do file
    for line in eachline(file)
           check = Dict()
                        has2 = false
                        has3 = false
                        for c in line
                                val = get(check, c, 0)
                                check[c] = val +1
                        end
                println(line)
                    for k in keys(check)
                            if check[k] == 2 && !has2
                                    has2 = true
                                    tot2 += 1
                                    print(k)
                                    print(check[k])
                            end
                            if check[k] == 3 && !has3
                                    has3 = true
                                    tot3 += 1
                                    print(k)
                                    print(check[k])

                            end
                    end
println("")
            end
    end

        println("-----")
        println(tot2)
        println(tot3)
        println(tot2 * tot3)

end


text:
using Term
import Term.renderables: Renderable
import Term: int
import Term.color: hsl2rgb

function make_julia_circles()
    circle = """
    oooo
 oooooooooo
oooooooooooo
oooooooooooo
 oooooooooo
    oooo    """

    # create circles
    green = Panel(
        RenderableText(circle; style="#389826 bold");
        fit=true,
        style = "dim #389826",
        justify = :center,
        padding=(2, 2, 0, 0),
        title = "Made",
        title_style = "default bold red",
    )
    red = Panel(
        RenderableText(circle; style="#CB3C33 bold");
        fit=true,
        style = "dim #CB3C33",
        justify = :center,
        padding=(2, 2, 0, 0),
        subtitle = "with",
        subtitle_style = "default bold #b656e3",
        subtitle_justify = :right,
    )
    purple = Panel(
        RenderableText(circle; style="#9558B2 bold");
        fit=true,
        style = "dim #9558B2",
        justify = :center,
        padding=(2, 2, 0, 0),
        subtitle = "Term",
        subtitle_style = "default bold #389826",
    )

    hspacer = Spacer(green.measure.w / 2 + 1, green.measure.h)
    line = Spacer(green.measure.w * 2 + 6, 1)
    circles =  line / (hspacer * green * hspacer) / (red * Spacer(2, purple.measure.h) * purple)
    return circles

end

function make_rgb_colors(; max_width=88)
    colors = ""

    for y in 0:5
        for x in 0:max_width
            h = x / max_width
            l = 0.1 + (((5-y) / 5) * 0.7)
            color = hsl2rgb(h*360, .9, l)
            bg = hsl2rgb(h*360, .9, l + 0.7/10)

            colors *= "[$color on_$bg]▄[/$color on_$bg]"
            # colors *= "[$color on_$bg]▬▄x[/$color on_$bg]"
        end
        colors *= "\n"
    end

    return colors
end

function rainbow_maker(N)
    R = range(30, 200, length=N)
    G = range(200, 60, length=N)
    B = range(50, 200, length=N)
    colors = []
    for n in 1:N
        r, g, b = int(R[n]), int(G[n]), int(B[n])
        push!(colors, "($r,$g,$b)")
    end
    return colors
end

circles = make_julia_circles()  # 42 x 17

_code_style = "yellow italic bold"

bfc = rainbow_maker(9)
basic_features = Panel("""
    [bright_red bold underline]Features[/bright_red bold underline]

[bold $(bfc[1])]✔[/bold $(bfc[1])][white] [blue]Colored text[/blue][/white]
[bold $(bfc[2])]✔[/bold $(bfc[2])][white] [italic]italic[/italic], [bold]bold[/bold], [underline]underline[/underline], [striked]striked[/striked], [inverse]inverse[/inverse][/white]
[bold $(bfc[3])]✔[/bold $(bfc[3])][white] styling [$_code_style]@macros[/$_code_style][/white]
[bold $(bfc[4])]✔[/bold $(bfc[4])][white] [italic white]markup[/italic white] style syntax[/white]
[bold $(bfc[5])]✔[/bold $(bfc[5])][white] progress bars[/white]
[bold $(bfc[6])]✔[/bold $(bfc[6])][white] [$_code_style]`Expr`[/$_code_style] and [$_code_style]`Type`[/$_code_style] introspection[/white]
[bold $(bfc[7])]✔[/bold $(bfc[7])][white] logging[/white]
[bold $(bfc[8])]✔[/bold $(bfc[8])][white] stacktraces[/white]
[bold $(bfc[9])]✔[/bold $(bfc[9])][white] syntax highlighting[/white]
"""; width=70, padding=(2, 2, 1, 2), justify=:center, style="default blue dim",
title="Term.jl", title_style="bright_blue bold",
subtitle="https://github.com/FedeClaudi/Term.jl", subtitle_style="default dim", subtitle_justify=:right
)

colors_info = TextBox(
"""[bold bright_green]Colors[/bold bright_green]
[bright_green bold]✔[/bright_green bold] 8-bit colors
[bright_green bold]✔[/bright_green bold] 16-bit colors
[bright_green bold]✔[/bright_green bold] hex colors
[bright_green bold]✔[/bright_green bold] rgb colors
[bright_green bold]✔[/bright_green bold] colors conversion
"""; width=25, padding=(0, 0, 0, 0))
colors = make_rgb_colors(;max_width=103)

_lorem = """朗眠裕安無際集正聞進士健音社野件草売規作独特認権価官家複入豚末告設悟自職遠氷育教載最週場仕踪持白炎組特曲強真雅立覧自価宰身訴側善論住理案者券真犯著避銀楽験館稿告
"""
lorem_description = TextBox(
"""[bold bright_yellow]Text reshaping[/bold bright_yellow]
[bright_yellow bold]✔[/bright_yellow bold] reshaping.
[bright_yellow bold]✔[/bright_yellow bold] justification
[bright_yellow bold]✔[/bright_yellow bold] Asian languages support
"""; width=28, padding=(0, 0, 0, 0))
lorem1 = TextBox(_lorem; width=62, padding=(0, 0, 0, 0))
lorem2 = TextBox(_lorem; width=42, padding=(0, 0, 0, 0))

expr = :(2x + 2π/θ)
tree = Renderable(sprint(typestree, Float64))
dendo = Renderable(sprint(inspect, expr))
renderables_info = TextBox(
"""[bold bright_blue]Renderables types[/bold bright_blue]
[bold bright_blue]✔[/bold bright_blue] Panel
[bold bright_blue]✔[/bold bright_blue] TextBox
[bold bright_blue]✔[/bold bright_blue] hLine
[bold bright_blue]✔[/bold bright_blue] vLine
[bold bright_blue]✔[/bold bright_blue] Tree
[bold bright_blue]✔[/bold bright_blue] Dendogram

[bold bright_cyan]Renderables layout[/bold bright_cyan]
[bold bright_cyan]✔[/bold bright_cyan] Horizontal stacking
[bold bright_cyan]✔[/bold bright_cyan] Vertical stacking
[bold bright_cyan]✔[/bold bright_cyan] Shorthand syntax: [bold bright_red]* /[/bold bright_red]

"""; width=28, padding=(0, 0, 0, 0)
)

line = hLine(140; style="bold dim grey35", box=:HEAVY)

# ----------------------------------- print ---------------------------------- #
# width is 140
print("\n"^3)

readme = (Spacer(10, circles.measure.h ) * circles * Spacer(8, circles.measure.h ) * basic_features) /
(Spacer(140, 2) / line) /
(Spacer(3, colors_info.measure.h ) * colors_info * Spacer(3, colors_info.measure.h ) * (Spacer(100, 1 ) / colors)) /
line /
(Spacer(3, lorem_description.measure.h ) * lorem_description * lorem1 * lorem2 ) /
line /
(Spacer(3, tree.measure.h) * renderables_info * tree * Spacer(5, tree.measure.h) * dendo)

print(Spacer(10, readme.measure.h) * readme)


text:
using MCPhylo
using Test

@testset "find_common_clusters" begin
    ref_tree = MCPhylo.parsing_newick_string("((A,(B,(C,(D,E)))),(F,(G,H)))")
    tree2 = MCPhylo.parsing_newick_string("((G,(C,(A,(F,E)))),(B,(D,H)))")
    MCPhylo.number_nodes!.([ref_tree, tree2])
    MCPhylo.set_binary!.([ref_tree, tree2])
    A, B, C, D, E, F, G, H = MCPhylo.find_by_name(tree2, "A"), MCPhylo.find_by_name(tree2, "B"),
                             MCPhylo.find_by_name(tree2, "C"), MCPhylo.find_by_name(tree2, "D"),
                             MCPhylo.find_by_name(tree2, "E"), MCPhylo.find_by_name(tree2, "F"),
                             MCPhylo.find_by_name(tree2, "G"), MCPhylo.find_by_name(tree2, "H")
     expected_dict = Dict{Int64,Tuple{Bool,Union{Missing, Float64}}}([(D.mother.num, (false, missing)), (C.mother.num, (false, missing)), (B.mother.num, (false, missing)),
                           (A.mother.num, (false, missing)), (G.mother.num, (false, missing)), (F.mother.num, (false, missing)),
                           (G.mother.mother.num, (true, 1.0)), (A.num, (true, 1.0)), (B.num, (true, 1.0)), (C.num, (true, 1.0)),
                           (D.num, (true, 1.0)), (E.num, (true, 1.0)), (F.num, (true, 1.0)), (G.num, (true, 1.0)), (H.num, (true, 1.0))])
    @test isequal(MCPhylo.find_common_clusters(ref_tree, tree2) , expected_dict)

    tree3 = MCPhylo.parsing_newick_string("((A,(C,(D,(B,E)))),(G,(F,H)))")
    MCPhylo.number_nodes!(tree3)
    MCPhylo.set_binary!(tree3)
    A, B, C, D, E, F, G, H = MCPhylo.find_by_name(tree3, "A"), MCPhylo.find_by_name(tree3, "B"),
                             MCPhylo.find_by_name(tree3, "C"), MCPhylo.find_by_name(tree3, "D"),
                             MCPhylo.find_by_name(tree3, "E"), MCPhylo.find_by_name(tree3, "F"),
                             MCPhylo.find_by_name(tree3, "G"), MCPhylo.find_by_name(tree3, "H")
     expected_dict = Dict([(D.mother.num, (false, missing)), (C.mother.num, (true, 1.0)), (B.mother.num, (false, missing)),
                           (A.mother.num, (true, 1.0)), (G.mother.num, (true, 1.0)), (F.mother.num, (false, missing)),
                           (A.mother.mother.num, (true, 1.0)), (A.num, (true, 1.0)), (B.num, (true, 1.0)), (C.num, (true, 1.0)),
                           (D.num, (true, 1.0)), (E.num, (true, 1.0)), (F.num, (true, 1.0)), (G.num, (true, 1.0)), (H.num, (true, 1.0))])
    @test isequal(MCPhylo.find_common_clusters(ref_tree, tree3), expected_dict)

    tree4 = MCPhylo.parsing_newick_string("((A,(B,(C,(D,E)))),(F,(G,H)))")
    MCPhylo.number_nodes!(tree4)
    MCPhylo.set_binary!(tree4)
    A, B, C, D, E, F, G, H = MCPhylo.find_by_name(tree4, "A"), MCPhylo.find_by_name(tree4, "B"),
                             MCPhylo.find_by_name(tree4, "C"), MCPhylo.find_by_name(tree4, "D"),
                             MCPhylo.find_by_name(tree4, "E"), MCPhylo.find_by_name(tree4, "F"),
                             MCPhylo.find_by_name(tree4, "G"), MCPhylo.find_by_name(tree4, "H")
    expected_dict = Dict([(D.mother.num, (true, 1.0)), (C.mother.num, (true, 1.0)), (B.mother.num, (true, 1.0)),
                          (A.mother.num, (true, 1.0)), (G.mother.num, (true, 1.0)), (F.mother.num, (true, 1.0)),
                          (A.mother.mother.num, (true, 1.0)), (A.num, (true, 1.0)), (B.num, (true, 1.0)), (C.num, (true, 1.0)),
                          (D.num, (true, 1.0)), (E.num, (true, 1.0)), (F.num, (true, 1.0)), (G.num, (true, 1.0)), (H.num, (true, 1.0))])
    @test isequal(MCPhylo.find_common_clusters(ref_tree, tree4), expected_dict)

    tree5 = MCPhylo.parsing_newick_string("((G,(X,(A,(F,E)))),(B,(D,H)))")
    MCPhylo.number_nodes!(tree5)
    MCPhylo.set_binary!(tree5)
    @test_throws ArgumentError MCPhylo.find_common_clusters(ref_tree, tree5)

    tree6 = MCPhylo.parsing_newick_string("(X,(G,(C,(A,(F,E)))),(B,(D,H)))))")
    MCPhylo.number_nodes!(tree5)
    MCPhylo.set_binary!(tree5)
    @test_throws ArgumentError MCPhylo.find_common_clusters(ref_tree, tree6)
end

@testset "one_way_compatible" begin
    tree = MCPhylo.parsing_newick_string("((A,C,E),(B,D))")
    tree2 = MCPhylo.parsing_newick_string("((A,C),(B,D,E))")
    expected_tree = MCPhylo.parsing_newick_string("(A,C,E,(B,D))")
    MCPhylo.number_nodes!.([tree, tree2, expected_tree])
    MCPhylo.set_binary!.([tree, tree2, expected_tree])
    @test MCPhylo.newick(MCPhylo.one_way_compatible(tree, tree2)) == MCPhylo.newick(expected_tree)
end

@testset "get_leaf_ranks" begin
    ref_tree = MCPhylo.parsing_newick_string("((A,(B,(C,(D,E)))),(F,(G,H)))")
    nodes = post_order(ref_tree)
    @test MCPhylo.get_leaf_ranks(nodes) == Dict([("A", 1), ("B", 2), ("C", 3),
                                                 ("D", 4), ("E", 5), ("F", 6),
                                                 ("G", 7), ("H", 8)])
end

@testset "order_tree!" begin
    tree = MCPhylo.parsing_newick_string("(A,B,(C,(D,E)F)G)H;")
    MCPhylo.number_nodes!(tree)
    MCPhylo.set_binary!(tree)
    A, B, C, D, E, F, G, H = MCPhylo.find_by_name(tree, "A"), MCPhylo.find_by_name(tree, "B"),
                             MCPhylo.find_by_name(tree, "C"), MCPhylo.find_by_name(tree, "D"),
                             MCPhylo.find_by_name(tree, "E"), MCPhylo.find_by_name(tree, "F"),
                             MCPhylo.find_by_name(tree, "G"), MCPhylo.find_by_name(tree, "H")
    cluster_start_indeces = Dict([(A, 3), (B, 7), (C, 2), (D, 8),
                                  (E, 5), (F, 1), (G, 4), (H, 6)])
    ordered_tree = MCPhylo.parsing_newick_string("(A,((E,D)F,C)G,B)H;")
    MCPhylo.number_nodes!(ordered_tree)
    MCPhylo.set_binary!(ordered_tree)
    @test MCPhylo.order_tree!(tree, cluster_start_indeces) == [A, E, D, C, B]
    MCPhylo.number_nodes!(tree)
    MCPhylo.set_binary!(tree)
    @test MCPhylo.newick(tree) == MCPhylo.newick(ordered_tree)
end

@testset "max/min_leaf_rank" begin
    tree = MCPhylo.parsing_newick_string("(A,B,(C,(D,E)F)G)H;")
    F, G, H, A = MCPhylo.find_by_name(tree, "F"), MCPhylo.find_by_name(tree, "G"),
                 MCPhylo.find_by_name(tree, "H"), MCPhylo.find_by_name(tree, "A")
    leaf_ranks = Dict([("A", 1), ("B", 5), ("C", 2), ("D", 4), ("E", 3)])

    @testset "min_leaf_rank" begin
        @test MCPhylo.min_leaf_rank(leaf_ranks, F) == 3
        @test MCPhylo.min_leaf_rank(leaf_ranks, G) == 2
        @test MCPhylo.min_leaf_rank(leaf_ranks, H) == 1
        @test MCPhylo.min_leaf_rank(leaf_ranks, A) == 1
    end

    @testset "max_leaf_rank" begin
        @test MCPhylo.max_leaf_rank(leaf_ranks, F) == 4
        @test MCPhylo.max_leaf_rank(leaf_ranks, G) == 4
        @test MCPhylo.max_leaf_rank(leaf_ranks, H) == 5
        @test MCPhylo.max_leaf_rank(leaf_ranks, A) == 1
    end
end

@testset "x_left_right" begin
    tree = MCPhylo.parsing_newick_string("(A,B,(C,(D,E)F)G)H;")
    MCPhylo.set_binary!(tree)
    MCPhylo.number_nodes!(tree)
    A, B, C, D, E, F, G, H = MCPhylo.find_by_name(tree, "A"), MCPhylo.find_by_name(tree, "B"),
                             MCPhylo.find_by_name(tree, "C"), MCPhylo.find_by_name(tree, "D"),
                             MCPhylo.find_by_name(tree, "E"), MCPhylo.find_by_name(tree, "F"),
                             MCPhylo.find_by_name(tree, "G"), MCPhylo.find_by_name(tree, "H")

    @testset "x_left" begin
        @test MCPhylo.x_left(A) == (H, [A,H])
        @test MCPhylo.x_left(B) == (B, [B,H])
        @test MCPhylo.x_left(C) == (G, [C,G,H])
        @test MCPhylo.x_left(D) == (F, [D,F,G])
        @test MCPhylo.x_left(E) == (E, [E,F])
    end

    @testset "x_right" begin
        @test MCPhylo.x_right(A) == (A, [A,H])
        @test MCPhylo.x_right(B) == (B, [B,H])
        @test MCPhylo.x_right(C) == (C, [C,G])
        @test MCPhylo.x_right(D) == (D, [D,F])
        @test MCPhylo.x_right(E) == (H, [E,F,G,H])
    end
end

@testset "majority_consensus_tree" begin
    tree1 = MCPhylo.parsing_newick_string("(((A,B),C),(D,E))")
    tree2 = MCPhylo.parsing_newick_string("((A,C),(B,D,E))")
    tree3 = MCPhylo.parsing_newick_string("(((B,C),A),D,E)")
    trees = [tree1, tree2, tree3]
    MCPhylo.number_nodes!.(trees)
    MCPhylo.set_binary!.(trees)
    result = MCPhylo.newick(MCPhylo.parsing_newick_string("((A,B,C),D,E)"))
    @test MCPhylo.newick(MCPhylo.majority_consensus_tree(trees)) == result

    tree4 = MCPhylo.parsing_newick_string("(((B,C),A),D,F)")
    push!(trees, tree4)
    @test_throws ArgumentError MCPhylo.majority_consensus_tree(trees)

    tree1 = MCPhylo.parsing_newick_string("(((A,B)AB,C)ABC,(D,E)DE)no_name;")
    tree2 = MCPhylo.parsing_newick_string("((A,C)AC,(B,D,E)BDE)no_name;")
    tree3 = MCPhylo.parsing_newick_string("(((B,C)BC,A)BCA,D,E)no_name;")
    trees = [tree1, tree2, tree3]
    MCPhylo.number_nodes!.(trees)
    MCPhylo.set_binary!.(trees)
    result = newick(MCPhylo.parsing_newick_string("((A,B,C),D,E)"))
    @test newick(MCPhylo.majority_consensus_tree(trees)) == result
end

@testset "loose_consensus_tree" begin
    tree1 = MCPhylo.parsing_newick_string("(((A,B),C),(D,E))")
    tree2 = MCPhylo.parsing_newick_string("((A,C),(B,D,E))")
    tree3 = MCPhylo.parsing_newick_string("(((B,C),A),D,E)")
    trees = [tree1, tree2, tree3]
    MCPhylo.number_nodes!.(trees)
    MCPhylo.set_binary!.(trees)
    result = newick(MCPhylo.parsing_newick_string("(A,B,C,(D,E))"))
    @test newick(MCPhylo.loose_consensus_tree(trees)) == result

    tree4 = MCPhylo.parsing_newick_string("(((B,C),A),D,F)")
    push!(trees, tree4)
    @test_throws ArgumentError MCPhylo.loose_consensus_tree(trees)

    tree1 = MCPhylo.parsing_newick_string("(((A,B)AB,C)ABC,(D,E)DE)no_name;")
    tree2 = MCPhylo.parsing_newick_string("((A,C)AC,(B,D,E)BDE)no_name;")
    tree3 = MCPhylo.parsing_newick_string("(((B,C)BC,A)BCA,D,E)no_name;")
    trees = [tree1, tree2, tree3]
    MCPhylo.number_nodes!.(trees)
    MCPhylo.set_binary!.(trees)
    result = newick(MCPhylo.parsing_newick_string("(A,B,C,(D,E)no_name)"))
    @test newick(MCPhylo.loose_consensus_tree(trees)) == result
end

@testset "greedy_consensus_tree" begin
    tree1 = MCPhylo.parsing_newick_string("(((A,B),C),(D,E))")
    tree2 = MCPhylo.parsing_newick_string("((A,C),(B,D,E))")
    tree3 = MCPhylo.parsing_newick_string("(((B,C),A),D,E)")
    trees = [tree1, tree2, tree3]
    MCPhylo.number_nodes!.(trees)
    MCPhylo.set_binary!.(trees)
    result = newick(MCPhylo.parsing_newick_string("(((A,B),C),(D,E))"))
    # result = newick(MCPhylo.parsing_newick_string("(((A,C),B),(D,E))"))
    # result = newick(MCPhylo.parsing_newick_string("(((B,C),A),(D,E))"))
    @test newick(MCPhylo.greedy_consensus_tree(trees)) == result
end

@testset "count_cluster_occurences" begin
    bit_vec = BitArray.([[1,0,0,0,0], [0,1,0,0,0], [0,0,1,0,0], [0,0,0,1,0],
              [0,0,0,0,1], [1,1,0,0,0], [1,1,1,0,0], [0,0,0,1,1], [1,0,0,0,0],
              [0,1,0,0,0], [0,0,1,0,0], [0,0,0,1,0], [0,0,0,0,1], [1,0,1,0,0],
              [0,1,0,1,1], [0,1,1,0,0], [1,1,1,0,0], [0,0,0,1,1], [1,0,0,0,0],
              [0,1,0,0,0], [0,0,1,0,0], [0,0,0,1,0], [0,0,0,0,1]])
    sort!(bit_vec, alg=QuickSort)
    queue = MCPhylo.PriorityQueue{BitVector, Int64}()
    queue[BitArray([0,0,0,0,1])] = -3
    queue[BitArray([0,0,0,1,0])] = -3
    queue[BitArray([0,1,0,0,0])] = -3
    queue[BitArray([0,0,1,0,0])] = -3
    queue[BitArray([1,0,0,0,0])] = -3
    queue[BitArray([0,0,0,1,1])] = -2
    queue[BitArray([1,1,0,0,0])] = -1
    queue[BitArray([1,0,1,0,0])] = -1
    queue[BitArray([0,1,0,1,1])] = -1
    queue[BitArray([0,1,1,0,0])] = -1
    queue[BitArray([1,1,1,0,0])] = -2

    @test MCPhylo.count_cluster_occurences(bit_vec) == queue
end

"""
# Additional test with a big file of trees
trees = MCPhylo.ParseNewick("./doc/Tree/Drav_mytrees_1.nwk")
MCPhylo.set_binary!.(trees)
MCPhylo.number_nodes!.(trees)
majority_tree = MCPhylo.majority_consensus_tree(trees)
println(newick(majority_tree))
"""


text:
# Qualified names → objects

function getthing(mod::Module, name::Vector, default = nothing)
  thing = mod
  for sym in name
    sym = Symbol(sym)
    if isdefined(thing, sym)
      thing = getfield(thing, sym)
    else
      return default
    end
  end
  return thing
end

getthing(mod::Module, name::AbstractString, default = nothing) =
  name == "" ?
    default :
    @as _ name split(_, ".", keep=false) map(Symbol, _) getthing(mod, _, default)

getthing(mod::Module, ::Void, default) = default

getthing(args...) = getthing(Main, args...)

# include_string with line numbers

function Base.include_string(s::AbstractString, fname::AbstractString, line::Integer)
  include_string("\n"^(line-1)*s, fname)
end

function Base.include_string(mod::Module, args...)
  eval(mod, :(include_string($(args...))))
end

# Get the current module for a file/pos

function getmodule(code, pos; filemod = nothing)
  codem = codemodule(code, pos)
  modstr = (codem != "" && filemod != nothing) ? "$filemod.$codem" :
           codem == "" ? filemod : codem
  getthing(modstr, Main)
end


text:
using Printf
using Plots; pyplot(fmt = :png)
using SolidStateDetectors
using Unitful
using StatsPlots

include("WaveformAnalysis.jl")

figSize = (1400, 1000)
extension = ".png"

function GetDetectorFacePosition(sim::Simulation)
    #thickness = sim.detector.semiconductors[1].geometry.h
    thickness = 2e-3
    #z = sim.detector.semiconductors[1].geometry.translate.z
    z = thickness/2
    zBottom = z - thickness/2
    zTop = z + thickness/2

    return (zBottom, zTop)
end

function PlotEvents(events::Vector{Event}, dir::String="", prefix::String="", suffix::String="")::Nothing
    if !isempty(events)
        try
            for i in 1:size(events[1].waveforms, 1) # Loop over contacts
                # Plot integrated charge
                pIC = plot(events[1].waveforms[1].time / 1u"ns", [x.waveforms[i].value for x in events],
                xlabel="Time [ns]", ylabel="Integrated charge [arb.]", linecolor="red",
                linealpha=1/sqrt(length(events)), legend=false)

                #savefig(dir * prefix * "integrated_charge_contact_$i" * suffix * extension)

                # Histogram rise times

                riseTimes = CollectRiseTimes(events)
                pRT = histogram(riseTimes / 1u"ns", bins=:scott, xlabel="Rise time [ns]", ylabel="Counts [arb.]")

                #savefig(dir * prefix * "risetime_hist_$i" * suffix * extension)

                # Histogram maximal integrated current

                maxAmplitude = [maximum(abs.(x.waveforms[i].value)) for x in events]
                # println(maxAmplitude)
                # println(length(riseTimes), length(maxAmplitude))
                pMA = histogram(maxAmplitude, bins=:scott, xlabel="Maximal amplitude", ylabel="Counts [arb.]")

                #savefig(dir * prefix * "amplitude_hist_$i" * suffix * extension)

                # 2D Amplitude - Rise time interpolated scatter plot
                pKDE = marginalkde(riseTimes / 1u"ns", maxAmplitude)

                #savefig(dir * prefix * "rise_time_amplitude_kde_$i" * suffix * extension)

                plot(pRT, pIC, pKDE, pMA, layout=(2, 2), size=figSize)

                savefig(dir * prefix * "event_contact_$i" * suffix * extension)
            end
        catch

        end
    end
end

function PlotGeometry(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    @info "Plotting detector geometry"
    plot(sim.detector, size=figSize)

    savefig(dir * prefix * "geometry" * suffix * extension)
end

function PlotElectricPotential(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    zBottom, thickness = GetDetectorFacePosition(sim)
    zTop = zBottom + thickness
    if !ismissing(sim.electric_potential)
        plot(plot(sim.electric_potential, y=0, aspect_ratio = :none), plot(sim.electric_potential, x=0, aspect_ratio = :none),
        plot(sim.electric_potential, z=zBottom), plot(sim.electric_potential, z=zTop),
        layout=(2, 2), size=figSize)

        savefig(dir * prefix * "electric_potential" * suffix * extension)
    end
end

function Plot2DWeightingPotential(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    zBottom, zTop = GetDetectorFacePosition(sim)
    for i in 1:size(sim.weighting_potentials, 1)
        if !ismissing(sim.weighting_potentials[i])
            plot(plot(sim.weighting_potentials[i], y=0, aspect_ratio = :none), plot(sim.weighting_potentials[i], x=0, aspect_ratio = :none),
            plot(sim.weighting_potentials[i], z=zBottom), plot(sim.weighting_potentials[i], z=zTop),
            layout=(2, 2), size=figSize)

            p = @sprintf "weighting_potential_%d" i

            savefig(dir * prefix * p * suffix * extension)
        end
    end
end

function PlotElectricField(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    zBottom, zTop = GetDetectorFacePosition(sim)
    if !ismissing(sim.electric_field)
        p_xz = plot(sim.electric_field, y=0, aspect_ratio = :none)
        plot_electric_fieldlines!(sim)
        p_yz = plot(sim.electric_field, x=0, aspect_ratio = :none)
        plot_electric_fieldlines!(sim)
        p_xy_bottom = plot(sim.electric_field, z=zBottom)
        p_xy_top = plot(sim.electric_field, z=zTop)
        plot(p_xz, p_yz, p_xy_bottom, p_xy_top, layout=(2, 2), size=figSize)

        savefig(dir * prefix * "electric_field" * suffix * extension)
    end
end

function PlotFields(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    @info "Plotting fields"
    PlotElectricPotential(sim, dir, prefix, suffix)
    PlotWeightingPotential(sim, dir, prefix, suffix)
    PlotElectricField(sim, dir, prefix, suffix)
end

function PlotChargeDensity(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    zBottom, zTop = GetDetectorFacePosition(sim)
    if !ismissing(sim.ρ)
        plot(plot(sim.ρ, x=0, aspect_ratio = :none), plot(sim.ρ, y=0, aspect_ratio = :none), plot(sim.ρ, z=zBottom), plot(sim.ρ, z=zTop), layout=(2,2), size=figSize)
        savefig(dir * prefix * "rho" * suffix * extension)
    end
end

function PlotDielectric(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    zBottom, zTop = GetDetectorFacePosition(sim)
    if !ismissing(sim.ϵ)
        plot(plot(sim.ϵ, x=0, aspect_ratio = :none), plot(sim.ϵ, y=0, aspect_ratio = :none), plot(sim.ϵ, z=zBottom), plot(sim.ϵ, z=zTop), layout=(2,2), size=figSize)
        savefig(dir * prefix * "epsilon" * suffix * extension)
    end
end

function PlotPointType(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    zBottom, zTop = GetDetectorFacePosition(sim)
    if !ismissing(sim.point_types)
        plot(plot(sim.point_types, x=0, aspect_ratio = :none), plot(sim.point_types, y=0, aspect_ratio = :none), plot(sim.point_types, z=zBottom), plot(sim.point_types, z=zTop), layout=(2,2), size=figSize)
        savefig(dir * prefix * "point_types" * suffix * extension)
    end
end

function PlotMaterialProperties(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    @info "Plotting material properties"

    PlotChargeDensity(sim, dir, prefix, suffix)
    PlotDielectric(sim, dir, prefix, suffix)
    PlotPointType(sim, dir, prefix, suffix)

    zBottom, zTop = GetDetectorFacePosition(sim)

    if !(ismissing(sim.ρ) || ismissing(sim.ϵ) || ismissing(sim.point_types))
        plot(plot(sim.ρ, x=0), plot(sim.ϵ, x=0), plot(sim.point_types, x=0), layout=(1, 3), size=figSize)
        savefig(dir * prefix * "material_properties_yz" * suffix * extension)
        plot(plot(sim.ρ, y=0), plot(sim.ϵ, y=0), plot(sim.point_types, y=0), layout=(1, 3), size=figSize)
        savefig(dir * prefix * "material_properties_xz" * suffix * extension)
        plot(plot(sim.ρ, z=zBottom), plot(sim.ϵ, z=zBottom), plot(sim.point_types, z=zBottom), layout=(1, 3), size=figSize)
        savefig(dir * prefix * "material_properties_xy_bottom" * suffix * extension)
        plot(plot(sim.ρ, z=zTop), plot(sim.ϵ, z=zTop), plot(sim.point_types, z=zTop), layout=(1, 3), size=figSize)
        savefig(dir * prefix * "material_properties_xy_top" * suffix * extension)
    end
end

function PlotCV(sims::AbstractVector{Simulation}, dir::String="", prefix::String="", suffix::String="")::Nothing
    cv = zeros(size(sims, 1), 2)
    for i in 1:size(sims, 1)
        cv[i, 1] = sims[i].detector.contacts[1].potential
        cv[i, 2] = calculate_capacitance(sims[i]) / 1u"pF"
    end

    sort!(cv, dims=1)

    pCV = plot(cv[:, 1], cv[:, 2], xlabel="Voltage [V]", ylabel="Capacitance [pF]", size=figSize)

    savefig(dir * prefix * "cv" * suffix * extension)

    dCm2dV = zeros(size(sims, 1)-1, 2)
    for i in 1:size(dCm2dV, 1)
        dCm2dV[i, 1] = cv[i, 1]
        dCm2dV[i, 2] = ((cv[i+1, 2])^(-2)-cv[i, 2]^(-2))/(cv[i+1, 1] - cv[i, 1])
    end

    pDerivCV = plot(dCm2dV[:, 1], dCm2dV[:, 2], xlabel="Voltage [V]", ylabel="d(1/C^2)/dV", size=figSize)

    savefig(dir * prefix * "deriv_cv" * suffix * extension)

    plot(pCV, pDerivCV, layout=(1, 2), size=figSize)

    savefig(dir * prefix * "cv_deriv_cv" * suffix * extension)
end

function Plot1DWeightingPotential(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    # Edge and corner plots only work for hexagons
    @info "Plotting weighting potentials vs z"

    for i in 1:length(sim.weighting_potentials)
        grid = sim.weighting_potentials[i].grid
        gridX = collect(grid.x)
        gridY = collect(grid.y)
        # R = sim.detector.contacts[i].geometry.rOuter
        # r = sqrt(3)*R/2
        # edge = 0.9.*(r, 0.0) # only works for center pixel now, expand
        # corner = 0.9.*(r, R/2)

        edge = (0, 6e-3)
        corner = (6e-3, 0)

        thresholdX = maximum(gridX[2:end]-gridX[1:end-1])/2
        thresholdY = maximum(gridY[2:end]-gridY[1:end-1])/2

        println(thresholdX, " ", thresholdY)
        x0 = findlast(x->x==0, collect(grid.x))
        y0 = findlast(x->x==0, collect(grid.y))
        xc = findlast(x->abs(x-corner[1]) < thresholdX, collect(gridX))
        yc = findlast(x->abs(x-corner[2]) < thresholdY, collect(gridY))
        xe = findlast(x->abs(x-edge[1]) < thresholdX, collect(gridX))
        ye = findlast(x->abs(x-edge[2]) < thresholdY, collect(gridY))

        println(grid[x0, y0, 1])
        println(grid[xe, ye, 1])
        println(grid[xc, yc, 1])

        wp0 = sim.weighting_potentials[i][x0, y0, :]
        wpe = sim.weighting_potentials[i][x0, ye, :]
        wpc = sim.weighting_potentials[i][xc, y0, :]

        plot(collect(grid.z), [wp0, wpe, wpc], labels = ["Center" "Edge" "Corner"],
        xlabel = "Z [m]", ylabel="Weighting potential", size=figSize)

        savefig(dir * prefix * "weighting_potential_vs_z_$i" * suffix * extension)
    end
end

function PlotIDP(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    @info "Plotting IDP along z and top of detector"

    bulk = sim.detector.semiconductors[1]
    bulk_geom = bulk.geometry
    center = bulk_geom.translate === missing ? (0,0) : bulk_geom.translate
    R = bulk_geom.r_interval.right
    H = bulk_geom.z_interval.right

    zs = 0:1e-8:H
    charge_vecz = similar(zs)
    for iz in eachindex(zs)
        pt = CartesianPoint{T}(0,0,zs[iz])
        charge_vecz[iz] = pt in bulk?
                            SolidStateDetectors.get_charge_density(bulk, pt) : 0
    end

    xs = -R:1e-7:R
    charge_vecx = similar(xs)
    for ix in eachindex(xs)
        pt = CartesianPoint{T}(xs[ix]+center[1], 0, H-1e-6)
        charge_vecx[ix] = pt in bulk?
                            SolidStateDetectors.get_charge_density(bulk, pt) : 0
    end

    plot(plot(zs, charge_vecz), plot(xs, charge_vecx), layout=(1, 2), size=figSize)
    savefig(dir * prefix * "IDP_z_and_top" * suffix * extension)
end

max_stars_repo_name: rjtayl/NDP

max_issues_repo_name: rjtayl/NDP

max_forks_repo_name: rjtayl/NDP


text:
using Documenter

# Build the docs
makedocs(
    sitename = "Documenter Testing",
    pages = [
        "Home" => "index.md",
    ],
    format = Documenter.HTML(
    ),
)

deploydocs(
    push_preview = true,
    versions = nothing,
    repo = "github.com/fredrikekre/documenter-testing.git",
)

max_stars_repo_name: fredrikekre/documenter-testing

max_issues_repo_name: fredrikekre/documenter-testing

max_forks_repo_name: fredrikekre/documenter-testing


text:
###############################################################################
#
#    Copyright (C) 2015  VoxForge
#
#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
###############################################################################

function mktrihed(monophones, triphones, mktri)
  monophones_arr=open(readlines, monophones)  # automatically closes file handle

  hed=open(mktri, "w")
  write(hed, "CL $triphones\n")
  for phoneln=monophones_arr
    phone=chomp(phoneln)
    if length(phone)>0
      write(hed,"TI T_$phone {(*-$phone+*,$phone+*,*-$phone).transP}\n")
    end
  end
  close(hed)
end

# if called from command line
if length(ARGS) > 0
  if ! isfile(ARGS[1])
    error("can't find monophones file: $ARGS[1]")
  end
  if ! isfile(ARGS[2])
    error("can't find triphones file: $ARGS[2]")
  end
  if length(ARGS) > 3
    error("prompts2list: too many arguments for call from command line")
  end

  mktrihed(ARGS[1], ARGS[2], ARGS[3])
end

max_stars_repo_name: Hyperparticle/lct-master

max_issues_repo_name: Hyperparticle/lct-master

max_forks_repo_name: Hyperparticle/lct-master


text:
# An IntervalCollection is an efficiently stored and indexed set of annotated
# genomic intervals. It looks something like this.
#
#                                      ┌─────┐
#                                      │trees│
#                                      └─────┘
#                                         │
#                              ┌──────────┴──────────┬────────────┐
#                              ▼                     ▼            │
#  Each sequence has       ┌──────┐              ┌──────┐         ▼
#    an associated         │ chr1 │              │ chr2 │
#    IntervalTree     ┌────┴──────┴────┐    ┌────┴──────┴────┐    ...
#                     │ IntervalTree 1 │    │ IntervalTree 2 │
#                     └────────────────┘    └────────────────┘
#                              │
#                    ┌─────────┴─────────────────────────┬────────────────────────┐
#                    │                                   │                        │
#                    ▼                                   ▼                        ▼
#   ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓  ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
#   ┃ Interval{T}(10000, 20000, ...) ┃  ┃ Interval{T}(35000, 40000, ...) ┃      ...
#   ┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛  ┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛
#
#
#
#
#                               ┌────────────────┐
#       ordered_trees holds     │ IntervalTree 1 │
#         an array of the       ├────────────────┤
#       some IntervalTrees,     │ IntervalTree 2 │
#           ordered by          └────────────────┘
#       chromosome for fast
#       ordered iteration.             ...
#

typealias IntervalCollectionTree{T} IntervalTrees.IntervalTree{Int64,Interval{T}}

type IntervalCollection{T} <: IntervalStream{T}
    # Sequence name mapped to IntervalTree, which in turn maps intervals to
    # a list of metadata.
    trees::Dict{StringField, IntervalCollectionTree{T}}

    # Keep track of the number of stored intervals
    length::Int

    # A vector of values(trees) sorted on sequence name.
    # This is used to iterate intervals as efficiently as possible, but is only
    # updated as needed, indicated by the ordered_trees_outdated flag.
    ordered_trees::Vector{IntervalCollectionTree{T}}
    ordered_trees_outdated::Bool

    function IntervalCollection()
        return new(Dict{StringField, IntervalCollectionTree{T}}(), 0,
                   IntervalCollectionTree{T}[], false)
    end

    # bulk insertion
    function IntervalCollection(intervals::AbstractVector{Interval{T}}, sort=false)
        if sort
            sort!(intervals)
        else
            if !issorted(intervals)
                error("Intervals must be sorted, or `sort=true` set, to construct an IntervalCollection")
            end
        end

        n = length(intervals)
        trees = Dict{StringField, IntervalCollectionTree{T}}()
        i = 1
        while i <= n
            j = i
            while j <= n && intervals[i].seqname == intervals[j].seqname
                j += 1
            end
            trees[intervals[i].seqname] = IntervalCollectionTree{T}(view(intervals, i:j-1))
            i = j
        end
        return new(trees, n, IntervalCollectionTree{T}[], true)
    end
end

function IntervalCollection{T}(intervals::AbstractVector{Interval{T}}, sort=false)
    return IntervalCollection{T}(intervals, sort)
end

function IntervalCollection{T}(interval_stream::IntervalStream{T})
    intervals = collect(Interval{T}, interval_stream)
    return IntervalCollection{T}(intervals, true)
end

function IntervalCollection(intervals)
    return IntervalCollection(collect(intervals), true)
end

function update_ordered_trees!{T}(ic::IntervalCollection{T})
    if ic.ordered_trees_outdated
        ic.ordered_trees = collect(IntervalCollectionTree{T}, values(ic.trees))
        p = sortperm(collect(AbstractString, keys(ic.trees)), lt=isless)
        ic.ordered_trees = ic.ordered_trees[p]
        ic.ordered_trees_outdated = false
    end
end

function Base.push!{T}(ic::IntervalCollection{T}, i::Interval{T})
    if !haskey(ic.trees, i.seqname)
        tree = IntervalCollectionTree{T}()
        ic.trees[i.seqname] = tree
        ic.ordered_trees_outdated = true
    else
        tree = ic.trees[i.seqname]
    end
    push!(tree, i)
    ic.length += 1
    return ic
end

function Base.show(io::IO, ic::IntervalCollection)
    const max_entries = 8
    n_entries = length(ic)
    println(io, "IntervalCollection with $(n_entries) intervals:")
    if n_entries > 0
        for (k, i) in enumerate(ic)
            if k > 8
                break
            end
            println(io, "  ", i)
        end
        if n_entries > max_entries
            print(io, "  ⋮")
        end
    end
end

Base.length(ic::IntervalCollection) = ic.length

# Iterators
# --------

typealias IntervalCollectionTreeIteratorState{T} IntervalTrees.IntervalBTreeIteratorState{Int64,Interval{T},64}

immutable IntervalCollectionIteratorState{T}
    i::Int # index into ordered_trees
    tree_state::IntervalCollectionTreeIteratorState{T}

    function IntervalCollectionIteratorState(i::Int)
        return new(i)
    end

    function IntervalCollectionIteratorState(i::Int, tree_state)
        return new(i, tree_state)
    end
end

function Base.start{T}(ic::IntervalCollection{T})
    update_ordered_trees!(ic)
    i = 1
    while i <= length(ic.ordered_trees)
        tree_state = start(ic.ordered_trees[i])
        if !done(ic.ordered_trees[i], tree_state)
            return IntervalCollectionIteratorState{T}(i, tree_state)
        end
        i += 1
    end

    return IntervalCollectionIteratorState{T}(i)
end

function Base.next{T}(ic::IntervalCollection{T},
                      state::IntervalCollectionIteratorState{T})
    i = state.i
    value, tree_state = next(ic.ordered_trees[i], state.tree_state)

    if done(ic.ordered_trees[i], tree_state)
        i += 1
        while i <= length(ic.ordered_trees)
            tree_state = start(ic.ordered_trees[i])
            if !done(ic.ordered_trees[i], tree_state)
                break
            end
            i += 1
        end
    end

    return value, IntervalCollectionIteratorState{T}(i, tree_state)
end

function Base.done{T}(ic::IntervalCollection{T},
                      state::IntervalCollectionIteratorState{T})
    return state.i > length(ic.ordered_trees)
end

"""
Return an iterator over all intervals in the collection that overlap the query interval.
"""
function Base.intersect{T}(a::IntervalCollection{T}, b::Interval)
    if haskey(a.trees, b.seqname)
        return intersect(a.trees[b.seqname], b)
    else
        return IntervalTrees.IntervalIntersectionIterator{Int64, Interval{T}, 64}()
    end
end

type IntersectIterator{S, T}
    a_trees::Vector{IntervalTrees.IntervalBTree{Int64, Interval{S}, 64}}
    b_trees::Vector{IntervalTrees.IntervalBTree{Int64, Interval{T}, 64}}

    i::Int # index into a_trees/b_trees.
    intersect_iterator::IntervalTrees.IntersectionIterator{Int64, Interval{S}, 64, Interval{T}, 64}

    function IntersectIterator(a_trees, b_trees)
        return new(a_trees, b_trees)
    end
end

Base.iteratorsize(::IntersectIterator) = Base.SizeUnknown()

"Iterate over pairs of intersecting intervals in two IntervalCollections"
function Base.intersect{S, T}(a::IntervalCollection{S}, b::IntervalCollection{T})
    seqnames = collect(AbstractString, keys(a.trees) ∩ keys(b.trees))
    sort!(seqnames, lt=isless)

    a_trees = IntervalCollectionTree{S}[a.trees[seqname] for seqname in seqnames]
    b_trees = IntervalCollectionTree{T}[b.trees[seqname] for seqname in seqnames]

    return IntersectIterator{S, T}(a_trees, b_trees)
end

function Base.start{S,T}(it::IntersectIterator{S,T})
    i = 1
    while i <= length(it.a_trees)
        intersect_iterator = intersect(it.a_trees[i], it.b_trees[i])
        intersect_iterator_state = start(intersect_iterator)
        if !done(intersect_iterator, intersect_iterator_state)
            it.i = i
            it.intersect_iterator = intersect_iterator
            return nothing
        end
        i += 1
    end
    it.i = i

    return nothing
end

function Base.next{S,T}(it::IntersectIterator{S, T}, ::Void)
    intersect_iterator = it.intersect_iterator
    value, intersect_iterator_state = next(intersect_iterator, nothing)
    i = it.i
    if done(intersect_iterator, intersect_iterator_state)
        i += 1
        while i <= length(it.a_trees)
            intersect_iterator = intersect(it.a_trees[i], it.b_trees[i])
            intersect_iterator_state = start(intersect_iterator)
            if !done(intersect_iterator, intersect_iterator_state)
                break
            end
            i += 1
        end
    end
    it.i = i
    it.intersect_iterator = intersect_iterator

    return value, nothing
end

function Base.done{S, T}(it::IntersectIterator{S, T}, state)
    return it.i > length(it.a_trees)
end

"""
Iterate over pairs of intersections in an IntervalCollection versus an
IntervalStream.
"""
function Base.intersect(a::IntervalCollection, b::IntervalStreamOrArray)
    return IntervalCollectionStreamIterator(metadatatype(b), a, b)
end

immutable IntervalCollectionStreamIterator{S, T, TV}
    a::IntervalCollection{S}
    b::T
end

function IntervalCollectionStreamIterator{S, T}(TV::Type, a::IntervalCollection{S},
                                                b::T)
    return IntervalCollectionStreamIterator{S, T, TV}(a, b)
end

immutable IntervalCollectionStreamIteratorState{S, TS, TV}
    intersection::IntervalTrees.Intersection{Int64, Interval{S}, 64}
    b_state::TS
    b_value::Interval{TV}

    function IntervalCollectionStreamIteratorState(intersection, b_state, b_value)
        return new(intersection, b_state, b_value)
    end

    function IntervalCollectionStreamIteratorState()
        return new(IntervalTrees.Intersection{Int64, Interval{S}, 64}())
    end
end

function Base.:(==){T}(a::IntervalCollection{T}, b::IntervalCollection{T})
    if length(a) != length(b)
        return false
    end

    for (i, j) in zip(a, b)
        if i != j
            return false
        end
    end

    return true
end

# This mostly follows from SuccessiveTreeIntersectionIterator in IntervalTrees
function Base.start{S,T,TV}(it::IntervalCollectionStreamIterator{S,T,TV})
    b_state = start(it.b)

    # TODO: We need to figure out a way to know these types at compile time.
    TS = typeof(b_state)

    intersection = IntervalTrees.Intersection{Int64, Interval{S}, 64}()
    while !done(it.b, b_state)
        b_value, b_state = next(it.b, b_state)
        if haskey(it.a.trees, b_value.seqname)
            tree = it.a.trees[b_value.seqname]
            IntervalTrees.firstintersection!(tree, b_value, intersection)
            if intersection.index != 0
                return IntervalCollectionStreamIteratorState{S, TS, TV}(
                    intersection, b_state, b_value)
            end
        end
    end

    return IntervalCollectionStreamIteratorState{S, TS, TV}()
end

function Base.next{S,T,TS,TV}(it::IntervalCollectionStreamIterator{S,T,TV},
                              state::IntervalCollectionStreamIteratorState{S,TS,TV})
    intersection = state.intersection
    entry = intersection.node.entries[intersection.index]
    return_value = (entry, state.b_value)
    b_state = state.b_state
    b_value = state.b_value
    IntervalTrees.nextintersection!(
        intersection.node, intersection.index, b_value, intersection)
    while intersection.index == 0 && !done(it.b, b_state)
        b_value, b_state = next(it.b, b_state)
        if haskey(it.a.trees, b_value.seqname)
            tree = it.a.trees[b_value.seqname]
            IntervalTrees.firstintersection!(tree, b_value, intersection)
        end
    end

    return return_value, IntervalCollectionStreamIteratorState{S, TS, TV}(
                            intersection, b_state, b_value)
end

function Base.done{S,T,TS,TV}(it::IntervalCollectionStreamIterator{S,T,TV},
                              state::IntervalCollectionStreamIteratorState{S,TS,TV})
    return state.intersection.index == 0
end


text:
using LabelledArrays

export lorenz_lorenz_lorenz_bidir_forced

@inline @inbounds function eom_lorenz_lorenz_lorenz_bidir_forced(u, p, t)
    c_xy, c_yx, c_zx, c_zy, a₁, a₂, a₃, b₁, b₂, b₃, a₃, c₁, c₂, c₃ = (p...,)
    x₁, x₂, x₃, y₁, y₂, y₃, z₁, z₂, z₃ = (u...,)

    dx₁ = -a₁*(x₁ - x₂) + c_yx*(y₁ - x₁) + c_zx*(z₁ - x₁)
    dx₂ = -x₁*x₃ + a₂*x₁ - x₂
    dx₃ = x₁*x₂ - a₃*x₃

    dy₁ = -b₁*(y₁ - y₂) + c_xy*(x₁ - y₁) + c_zy*(z₁ - y₁)
    dy₂ = -y₁*y₃ + b₂*y₁ - y₂
    dy₃ = y₁*y₂ - b₃*y₃

    dz₁ = -c₁*(z₁ - z₂)
    dz₂ = -z₁*z₃ + c₂*z₁ - z₂
    dz₃ = z₁*z₂ - c₃*z₃

    return SVector{9}(dx₁, dx₂, dx₃, dy₁, dy₁, dy₃, dz₁, dz₂, dz₃)
end

"""
    lorenz_lorenz_lorenz_bidir_forced(; u0 = rand(9),
        c_xy = 0.1, c_yx = 0.1,
        c_zx = 0.05, c_zy = 0.05,
        a₁ = 10, a₂ = 28, a₃ = 8/3,
        b₁ = 10, b₂ = 28, b₃ = 8/3,
        c₁ = 10, c₂ = 28, c₃ = 8/3)

Initialise a system consisting of two bidirectionally coupled 3D Lorenz
systems forced by an external 3D Lorenz system, giving a 9D system.

## Equations of motion

The dynamics is generated by the following vector field

```math
\\begin{aligned}
\\dot{x_1} &= - a_1 (x_1 - x_2) + c_{yx}(y_1 - x_1) + c_{zx}(z_1 - x_1) \\\\
\\dot{x_2} &= - x_1 x_3 + a_2 x_1 - x_2 \\\\
\\dot{x_3} &= x_1 x_2 - a_3 x_3 \\\\
\\dot{y_1} &= -b_1 (y_1 - y_2) + c_{xy} (x_1 - y_1) + c_{zy}(z_1 - y_1) \\\\
\\dot{y_2} &= - y_1 y_3 + b_2 y_1 - y_2 \\\\
\\dot{y_3} &= y_1 y_2 - b_3 y_3 \\\\
\\dot{z_1} &= - c_1 (z_1 - z_2) \\\\
\\dot{z_2} &= - z_1 z_3 + c_2 z_1 - z_2 \\\\
\\dot{z_3} &= z_1 z_2 - c_3 z_3
\\end{aligned}
```

## References

1. Amigó, José M., and Yoshito Hirata. "Detecting directional couplings from
    multivariate flows by the joint distance distribution." Chaos: An
    Interdisciplinary Journal of Nonlinear Science 28.7 (2018): 075302.
"""
function lorenz_lorenz_lorenz_bidir_forced(; u0 = rand(9),
        c_xy = 0.1, c_yx = 0.1,
        c_zx = 0.05, c_zy = 0.05,
        a₁ = 10, a₂ = 28, a₃ = 8/3,
        b₁ = 10, b₂ = 28, b₃ = 8/3,
        c₁ = 10, c₂ = 28, c₃ = 8/3)

    p = @LArray [c_xy, c_yx, c_zx, c_zy, a₁, a₂, a₃, b₁, b₂, b₃, a₃, c₁, c₂, c₃] (:c_xy, :c_yx, :c_zx, :c_zy, :a₁, :a₂, :a₃, :b₁, :b₂, :b₃, :a₃, :c₁, :c₂, :c₃)
    ContinuousDynamicalSystem(eom_lorenz_lorenz_lorenz_bidir_forced, u0, p)
end

function lorenz_lorenz_lorenz_bidir_forced_trajectory(npts;
        n_transient = 2000, dt = 0.1, sample_dt = 1,
        u0 = rand(9),
        c_xy = 1.0, c_yx = 1.0, c_zx = 1.0, c_zy = 1.0, # beyond c = 2, systems syncronize
        a₁ = 10, a₂ = 28, a₃ = 8/3,
        b₁ = 10, b₂ = 28, b₃ = 8/3,
        c₁ = 10, c₂ = 28, c₃ = 8/3)

    s = lorenz_lorenz_lorenz_bidir_forced(u0 = u0,
        c_xy = c_xy, c_yx = c_yx,
        c_zx = c_zx, c_zy = c_zy,
        a₁ = a₁, a₂ = a₂, a₃ = a₃,
        b₁ = b₁, b₂ = b₂, b₃ = b₃,
        c₁ = c₁, c₂ = c₂, c₃ = c₃)

    # the system is recorded at times t0:dt:T
    T = npts*dt*sample_dt
    o = trajectory(s, T, dt = dt, Ttr = n_transient*dt, alg = SimpleDiffEq.SimpleATsit5())[1:sample_dt:end-1, :] #alg = SimpleDiffEq.SimpleATsit5()
end

function good_lorenz_lorenz_lorenz_bidir_forced_trajectory(npts;
        sample_dt = 1,  Ttr = 5000, dt = 0.1,
        Da₁ = Uniform(9.5, 10.5),
        Da₂ = Uniform(27.5, 28.5),
        Da₃ = Uniform(7.5/3, 8.5/3),
        Db₁ = Uniform(9.5, 10.5),
        Db₂ = Uniform(27.5, 28.5),
        Db₃ = Uniform(7.5/3, 8.5/3),
        Dc₁ = Uniform(9.5, 10.5),
        Dc₂ = Uniform(27.5, 28.5),
        Dc₃ = Uniform(7.5/3, 8.5/3),

        a₁ = nothing,
        a₂ = nothing,
        a₃ = nothing,
        b₁ = nothing,
        b₂ = nothing,
        b₃ = nothing,
        c₁ = nothing,
        c₂ = nothing,
        c₃ = nothing,
        c_xy = 0.2,  c_yx = 0.2,
        c_zx = 0.05, c_zy = 0.05,
        u0 = [rand(Uniform(0, 10)) for i = 1:9],
        n_maxtries = 300)

    n_tries = 0

    while n_tries <= n_maxtries
        a₁ == nothing ? a₁ = rand(Da₁) : a₁ = a₁
        a₂ == nothing ? a₂ = rand(Da₂) : a₂ = a₂
        a₃ == nothing ? a₃ = rand(Da₃) : a₃ = a₃
        b₁ == nothing ? b₁ = rand(Db₁) : b₁ = b₁
        b₂ == nothing ? b₂ = rand(Db₂) : b₂ = b₂
        b₃ == nothing ? b₃ = rand(Db₃) : b₃ = b₃
        c₁ == nothing ? c₁ = rand(Dc₁) : c₁ = c₁
        c₂ == nothing ? c₂ = rand(Dc₂) : c₂ = c₂
        c₃ == nothing ? c₃ = rand(Dc₃) : c₃ = c₃
        pts = lorenz_lorenz_lorenz_bidir_forced(npts,
            sample_dt = sample_dt, dt = dt, n_transient = Ttr,
            c_xy = c_xy,  c_yx = c_yx,
            c_zx = c_zx, c_zy = c_zy,
            a₁ = a₁, a₂ = a₂, a₃ = a₃,
            b₁ = b₁, b₂ = b₂, b₃ = b₃,
            c₁ = c₁, c₂ = c₂, c₃ = c₃)

        if all(Matrix(pts) .< 1e9) #&& length(unique(pts)) < length(pts)*0.8
            return pts
        end
        println("no attractor found. trying with new initial condition and parameters")
        n_tries += 1
    end
end


text:
using PackageCompiler
create_sysimage([:ExaPF, :Ipopt, :MadNLP, :ProxAL]; sysimage_path="proxal.so", precompile_execution_file="test/distributed64.jl")


text:
@testset "Flatten/unflatten" begin
    mlp_parameters = params(Chain(Dense(4, 16, σ), Dense(16, 4), softmax))
    mlp_parameters_count = mapreduce(length, +, mlp_parameters)

    flattened = mlp_parameters |> flatten
    @test length(flattened) == mlp_parameters_count

    unflattened = unflatten_set(flattened, mlp_parameters)
    for (p, u) in zip(mlp_parameters, unflattened)
        @test p == u
    end

    mlp_parameters_smol = params(Dense(4, 16, σ))
    @test_throws AssertionError unflatten_set(flattened, mlp_parameters_smol)

    conv_parameters = params(Chain(
        Conv((3, 3), 3=>16, pad=(1,1), relu), BatchNorm(16),
        Conv((3, 3), 16=>32, pad=(1,1), relu), BatchNorm(32),
    ))
    conv_parameters_count = mapreduce(length, +, conv_parameters)
    flattened = conv_parameters |> flatten
    @test length(flattened) == conv_parameters_count

    unflattened = unflatten_set(flattened, conv_parameters)
    for (p, u) in zip(conv_parameters, unflattened)
        @test p == u
    end
end


text:
####
#### PhysicsState
####

resettime!(phys::PhysicsState) = (reset!(phys.timer); phys.elapsedsim = 0; phys)

####
#### UIState
####

function alignscale!(ui::UIState, sim::MJSim)
    ui.cam[].lookat = sim.m.stat.center
    ui.cam[].distance = 1.5 * sim.m.stat.extent
    ui.cam[].type = MJCore.mjCAMERA_FREE
    return ui
end

####
#### Engine
####

@inline mode(e::Engine, idx::Integer = e.curmodeidx) = e.modes[idx]

function switchmode!(e::Engine, idx::Integer)
    teardown!(e.ui, e.phys, mode(e))
    deregister!(e.mngr, e.modehandlers...)

    e.curmodeidx = idx
    e.modehandlers = handlers(e.ui, e.phys, mode(e))
    setup!(e.ui, e.phys, mode(e))
    register!(e.mngr, e.modehandlers...)

    return e
end

function printhelp(e::Engine)
    io = e.ui.io1

    writedescription(io, e.handlers)
    handlerdescription = String(take!(io))

    writedescription(io, e.modehandlers)
    modehandlerdescription = String(take!(io))

    println("Standard Commands:")
    print(handlerdescription)
    if !isempty(modehandlerdescription)
        println("$(nameof(mode(e))) Mode Commands:")
        print(modehandlerdescription)
    end
    println()
    println()

    return
end

function writedescription(io, hs::Vector{EventHandler})
    if !isempty(hs)
        whens = String[]
        whats = String[]
        for h in hs
            if h.when !== nothing && h.what !== nothing
                push!(whens, h.when)
                push!(whats, h.what)
            elseif h.what !== nothing
                push!(whens, "----")
                push!(whats, h.what)
            end
        end

        header = ["Command", "Description"]
        _, ncols = get_terminalsize()
        w1max = max(maximum(length, whens), length(first(header)))
        w2max = max(maximum(length, whats), length(first(header)))

        w1 = min(w1max, div(ncols, 2))
        w2 = min(w2max, ncols - w1 - 4 * length(header)) # each column is padded by 4 spaces
        pretty_table(io, hcat(whens, whats), ["Command", "Description"], alignment = :L, linebreaks = true, autowrap = true, columns_width = [w1, w2])
    end

    return
end

function overlay_info(rect::MJCore.mjrRect, e::Engine)
    ui = e.ui
    io1 = ui.io1
    io2 = ui.io2
    phys = e.phys
    sim = getsim(phys.model)

    seekstart(io1)
    seekstart(io2)

    println(io1, "Mode")
    println(io2, nameof(mode(e)))

    println(io1, "Status")
    if ui.paused
        println(io2, "Paused")
    elseif ui.reversed
        println(io2, "Reverse Simulation")
    else
        println(io2, "Forward Simulation")
    end

    println(io1, "Time")
    @printf io2 "%.3f s\n" time(sim)

    println(io1, "Refresh Rate")
    @printf io2 "%d Hz\n" ui.refreshrate

    println(io1, "Resolution")
    @printf io2 "(%d, %d)\n" e.mngr.state.width e.mngr.state.height

    println(io1, "Sim Speed")
    if ui.speedmode
        if ui.speedfactor < 1
            @printf io2 "%.5gx (slower)\n" 1 / ui.speedfactor
        else
            @printf io2 "%.5gx (faster)\n" ui.speedfactor
        end
    else
        println(io2, "1")
    end

    println(io1, "Frame")
    println(io2, MJCore.mjFRAMESTRING[e.ui.vopt[].frame+1])

    println(io1, "Label")
    println(io2, MJCore.mjLABELSTRING[e.ui.vopt[].label+1])

    # env specific info
    if phys.model isa AbstractMuJoCoEnvironment
        name = string(Base.nameof(typeof(phys.model)))
        println(io1, "Env")
        println(io2, name)

        println(io1, "Reward")
        @printf io2 "%.5g\n" ui.reward
    end

    # mode specific info
    println(io1, "Mode Info")
    println(io2)
    modeinfo(io1, io2, ui, phys, mode(e))

    info1 = string(chomp(String(take!(io1))))
    info2 = string(chomp(String(take!(io2))))

    mjr_overlay(
        MJCore.FONT_NORMAL,
        MJCore.GRID_BOTTOMLEFT,
        rect,
        info1,
        info2,
        ui.con,
    )

    return
end

function startrecord!(e::Engine)
    window = e.mngr.state.window
    SetWindowAttrib(window, GLFW.RESIZABLE, 0)
    w, h = GLFW.GetFramebufferSize(window)
    resize!(e.framebuf, 3 * w * h)
    e.ffmpeghandle, e.videodst = startffmpeg(w, h, e.min_refreshrate)
    @info "Recording video. Window resizing temporarily disabled"
    return e
end

function recordframe(e::Engine)
    w, h = GLFW.GetFramebufferSize(e.mngr.state.window)
    rect = mjrRect(Cint(0), Cint(0), Cint(w), Cint(h))
    mjr_readPixels(e.framebuf, C_NULL, rect, e.ui.con)
    write(e.ffmpeghandle, e.framebuf)
    return nothing
end

function stoprecord!(e::Engine)
    SetWindowAttrib(e.mngr.state.window, GLFW.RESIZABLE, 1)
    @info "Recording finished, window resizing re-enabled. Waiting for transcoding to finish."
    ispaused = e.ui.paused
    setpause!(e.ui, e.phys, true)
    close(e.ffmpeghandle)
    @info "Finished recording! Video saved to $(e.videodst)"
    e.ffmpeghandle = e.videodst = nothing
    setpause!(e.ui, e.phys, ispaused)
    return e
end

function setpause!(ui::UIState, p::PhysicsState, status::Bool)
    status ? stop!(p.timer) : start!(p.timer)
    ui.paused = status
    ui
end


text:
using Interpolations # To interpolate value function
using Expectations # To easily find expected values
using Distributions # To create normal stock returns =)

# Load files

include("model/Structs.jl")
include("model/Fundamentals.jl")
include("model/MainFunctions.jl")


text:
@inline ldof(S::Real) = Int(2S+1)
function Sp(S::Real)
    S2 = Int(2S)
    ret = [0.5sqrt((S2-m2)*(S2+m2+2) ) for m2 in (S2-2):-2:(-S2)]
    return diagm(1=>ret)
end
function Sp(S::Real, L::Integer, i::Integer)
    leftN = ldof(S)^(i-1)
    rightN = ldof(S)^(L-i)
    kron(eye(leftN), kron(Sp(S), eye(rightN)))
end

function Sm(S::Real)
    S2 = Int(2S)
    ret = [0.5sqrt((S2+m2)*(S2-m2+2) ) for m2 in (S2):-2:(-S2+2)]
    return diagm(-1=>ret)
end
function Sm(S::Real, L::Integer, i::Integer)
    leftN = ldof(S)^(i-1)
    rightN = ldof(S)^(L-i)
    kron(eye(leftN), kron(Sm(S), eye(rightN)))
end

function Sz(S::Real)
    S2 = Int(2S)
    ret = [0.5m2 for m2 in S2:-2:-S2]
    return diagm(0=>ret)
end
function Sz(S::Real, L::Integer, i::Integer)
    leftN = ldof(S)^(i-1)
    rightN = ldof(S)^(L-i)
    kron(eye(leftN), kron(Sz(S), eye(rightN)))
end

function magnetization(S::Real,L::Integer,k::Integer=0)
    sz = Sz(S)
    S2 = Int(2S)
    ld = ldof(S)
    N = ld^L
    res = zeros(N)
    invV = 1.0/L
    for i in 1:N
        for (j,m) in enumerate(digits(i-1,base=ld,pad=L))
            res[i] += sz[m+1,m+1] * cospi(k*invV*(j-1))
        end
    end
    res .*= invV
    return diagm(0=>res)
end

function spinchain(S::Real, L::Integer, Jz::Real, Jxy::Real, h::Real, Guni::Real, Gstag::Real)
    S2 = Int(2S)
    ld = ldof(S)
    sz = Sz(S)
    sp = Sp(S)
    sm = Sm(S)
    sx = 0.5.*(sp.+sm)
    jxy = 0.5Jxy

    G = [-(Guni+Gstag), -(Guni-Gstag)]
    onsite = [ (-h) .* sz .+ G[i] .* sx for i in 1:2]

    N = ld^L
    H = zeros(N,N)
    for i=1:L
        H .+ kron(eye(ld^(i-1)), kron(onsite[mod1(i,2)], eye(ld^(L-i))))
    end

    bondH = Jz.*kron(sz,sz) .+ jxy.*(kron(sp,sm) .+ kron(sm,sp))
    H .+= Jz.*kron(sz, kron(eye(div(N,ld*ld)), sz))
    H .+= jxy .* kron(sm, kron(eye(div(N,ld*ld)), sp))
    H .+= jxy .* kron(sp, kron(eye(div(N,ld*ld)), sm))

    leftN = 1
    rightN = ld^(L-2)
    for i in 1:L-1
        H .+= kron(eye(leftN), kron(bondH, eye(rightN)))
        leftN *= ld
        rightN = div(rightN,ld)
    end
    return H
end

doc"""
\mathcal{H} = Jz \sum_i(S^z_i S^z_{i+1}) + 0.5Jxy \sum_i(S^+_i S^-_{i+1} + h.c.) - h \sum_i S^z_i - G \sum_i S^x_i
"""
struct SpinChainSolver <: Solver
    ef :: Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
    S  :: Float64
    L  :: Int
    function SpinChainSolver(S::Real, L::Integer
                             ;
                             Jz::Real = 1.0, Jxy::Real = 1.0,
                             h::Real = 0.0,
                             Guni::Real = 0.0, Gstag::Real = 0.0)
        new(eigen(spinchain(S, L, Jz, Jxy, h, Guni, Gstag)), S, L)
    end
end
creator(solver::SpinChainSolver) = Sp(solver.S)
creator(solver::SpinChainSolver, i::Integer) = Sp(solver.S, solver.L, i)
annihilator(solver::SpinChainSolver) = Sm(solver.S)
annihilator(solver::SpinChainSolver, i::Integer) = Sm(solver.S, solver.L, i)
basis(solver::SpinChainSolver) = Sz(solver.S)
basis(solver::SpinChainSolver, i::Integer) = Sz(solver.S, solver.L, i)
orderparameter(solver::SpinChainSolver, k::Integer = 0) = magnetization(solver.S, solver.L, k)


text:
const PyJuliaAnyValue_Type__ref = Ref(PyNULL)
PyJuliaAnyValue_Type() = begin
    ptr = PyJuliaAnyValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaBaseValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.AnyValue",
                base = base,
                repr = pyjlany_repr,
                str = pyjlany_str,
                getattro = pyjlany_getattro,
                setattro = pyjlany_setattro,
                call = pyjlany_call,
                iter = pyjlany_iter,
                richcompare = pyjlany_richcompare,
                as_mapping = (
                    length = pyjlany_length,
                    subscript = pyjlany_getitem,
                    ass_subscript = pyjlany_setitem,
                ),
                as_sequence = (contains = pyjlany_contains,),
                as_number = (
                    positive = pyjlany_positive,
                    negative = pyjlany_negative,
                    absolute = pyjlany_absolute,
                    power = pyjlany_power,
                    add = pyjlany_binop(+),
                    subtract = pyjlany_binop(-),
                    multiply = pyjlany_binop(*),
                    truedivide = pyjlany_binop(/),
                    divmod = pyjlany_binop((x, y) -> (fld(x, y), mod(x, y))),
                    floordivide = pyjlany_binop(fld),
                    remainder = pyjlany_binop(mod),
                    lshift = pyjlany_binop(<<),
                    rshift = pyjlany_binop(>>),
                    and = pyjlany_binop(&),
                    xor = pyjlany_binop(xor),
                    or = pyjlany_binop(|),
                ),
                methods = [
                    (name = "__dir__", flags = Py_METH_NOARGS, meth = pyjlany_dir),
                    (name = "_repr_mimebundle_", flags = Py_METH_VARARGS | Py_METH_KEYWORDS, meth = pyjlany_repr_mimebundle),
                    (name = "_jl_raw", flags = Py_METH_NOARGS, meth = pyjlany_toraw),
                    (name = "_jl_display", flags = Py_METH_NOARGS, meth = pyjlany_display),
                    (name = "_jl_help", flags = Py_METH_NOARGS, meth = pyjlany_help),
                ],
                getset = [(name = "__name__", get = pyjlany_name)],
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaAnyValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaAnyValue_New(x) = PyJuliaValue_New(PyJuliaAnyValue_Type(), x)
PyJuliaValue_From(x) = PyJuliaAnyValue_New(x)

pyjlany_repr(xo::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)
        s = "<jl $(repr(x))>"
        PyUnicode_From(s)
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlany_str(xo::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)
        s = string(x)
        PyUnicode_From(s)
    catch err
        PyErr_SetJuliaError(err)
        return PyNULL
    end

pyjlany_getattro(xo::PyPtr, ko::PyPtr) = begin
    # Try generic lookup first
    ro = PyObject_GenericGetAttr(xo, ko)
    if isnull(ro) && PyErr_IsSet(PyExc_AttributeError())
        PyErr_Clear()
    else
        return ro
    end
    # Now try to get the corresponding property
    x = PyJuliaValue_GetValue(xo)
    k = PyUnicode_AsString(ko)
    isempty(k) && PyErr_IsSet() && return PyNULL
    k = pyjl_attr_py2jl(k)
    try
        v = getproperty(x, Symbol(k))
        PyObject_From(v)
    catch err
        if !hasproperty(x, Symbol(k)) ||
           (err isa UndefVarError && err.var === Symbol(k)) ||
           (err isa ErrorException && occursin("has no field", err.msg))
            PyErr_SetStringFromJuliaError(PyExc_AttributeError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

propertytype(x, k) =
    propertiesarefields(typeof(x)) && hasfield(typeof(x), k) ? fieldtype(typeof(x), k) : Any
@generated propertiesarefields(::Type{T}) where {T} =
    which(getproperty, Tuple{T,Symbol}) == which(getproperty, Tuple{Any,Symbol})

pyjlany_setattro(xo::PyPtr, ko::PyPtr, vo::PyPtr) = begin
    # Try generic lookup first
    ro = PyObject_GenericSetAttr(xo, ko, vo)
    if ism1(ro) && PyErr_IsSet(PyExc_AttributeError())
        PyErr_Clear()
    else
        return ro
    end
    if isnull(vo)
        PyErr_SetString(PyExc_TypeError(), "attribute deletion not supported")
        return Cint(-1)
    end
    # Now try to set the corresponding property
    x = PyJuliaValue_GetValue(xo)
    k = PyUnicode_AsString(ko)
    isempty(k) && PyErr_IsSet() && return Cint(-1)
    k = pyjl_attr_py2jl(k)
    try
        V = propertytype(x, Symbol(k))
        ism1(PyObject_Convert(vo, V)) && return Cint(-1)
        v = takeresult(V)
        setproperty!(x, Symbol(k), v)
        Cint(0)
    catch err
        if !hasproperty(x, Symbol(k)) ||
           (err isa UndefVarError && err.var === Symbol(k)) ||
           (err isa ErrorException && occursin("has no field", err.msg))
            PyErr_SetStringFromJuliaError(PyExc_AttributeError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        Cint(-1)
    end
end

pyjl_dir(x) = propertynames(x)
pyjl_dir(x::Module) = begin
    r = Symbol[]
    append!(r, names(x, all = true, imported = true))
    for m in ccall(:jl_module_usings, Any, (Any,), x)::Vector
        append!(r, names(m))
    end
    r
end

pyjlany_dir(xo::PyPtr, _::PyPtr) = begin
    fo = PyObject_GetAttrString(PyJuliaBaseValue_Type(), "__dir__")
    isnull(fo) && return PyNULL
    ro = PyObject_CallNice(fo, PyObjectRef(xo))
    Py_DecRef(fo)
    isnull(ro) && return PyNULL
    x = PyJuliaValue_GetValue(xo)
    ks = try
        collect(map(string, pyjl_dir(x)))
    catch err
        Py_DecRef(ro)
        PyErr_SetJuliaError(err)
        return PyNULL
    end
    for k in ks
        ko = PyUnicode_From(pyjl_attr_jl2py(k))
        isnull(ko) && (Py_DecRef(ro); return PyNULL)
        err = PyList_Append(ro, ko)
        Py_DecRef(ko)
        ism1(err) && (Py_DecRef(ro); return PyNULL)
    end
    return ro
end

pyjlany_call(fo::PyPtr, argso::PyPtr, kwargso::PyPtr) = begin
    f = PyJuliaValue_GetValue(fo)
    if isnull(argso)
        args = Vector{Any}()
    else
        ism1(PyObject_Convert(argso, Vector{Any})) && return PyNULL
        args = takeresult(Vector{Any})
    end
    if isnull(kwargso)
        kwargs = Dict{Symbol,Any}()
    else
        ism1(PyObject_Convert(kwargso, Dict{Symbol,Any})) && return PyNULL
        kwargs = takeresult(Dict{Symbol,Any})
    end
    try
        x = f(args...; kwargs...)
        PyObject_From(x)
    catch err
        if err isa MethodError && err.f === f
            PyErr_SetStringFromJuliaError(PyExc_TypeError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

pyjlany_length(xo::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)
        Py_ssize_t(length(x))
    catch err
        if err isa MethodError && err.f === length
            PyErr_SetStringFromJuliaError(PyExc_TypeError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        Py_ssize_t(-1)
    end

@generated pyjl_keytype(::Type{T}) where {T} =
    try
        keytype(T)
    catch
        nothing
    end;
pyjl_hasvarindices(::Type) = true
pyjl_hasvarindices(::Type{<:AbstractDict}) = false

pyjl_getindices(x, ko) =
    if (K = pyjl_keytype(typeof(x))) !== nothing
        ism1(PyObject_Convert(ko, K)) ? PYERR() : (takeresult(K),)
    elseif pyjl_hasvarindices(typeof(x)) && PyTuple_Check(ko)
        ism1(PyObject_TryConvert(ko, Tuple)) ? PYERR() : takeresult(Tuple)
    else
        ism1(PyObject_TryConvert(ko, Any)) ? PYERR() : takeresult(Any)
    end

pyjlany_getitem(xo::PyPtr, ko::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)
    k = pyjl_getindices(x, ko)
    k === PYERR() && return PyNULL
    try
        PyObject_From(x[k...])
    catch err
        if err isa BoundsError && err.a === x
            PyErr_SetStringFromJuliaError(PyExc_IndexError(), err)
        elseif err isa KeyError && (err.key === k || (err.key,) === k)
            PyErr_SetStringFromJuliaError(PyExc_KeyError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end
end

@generated pyjl_valtype(::Type{T}) where {T} =
    try
        valtype(T)
    catch
        try
            eltype(T)
        catch
            nothing
        end
    end;

pyjl_getvalue(x, vo) =
    if (V = pyjl_valtype(typeof(x))) !== nothing
        ism1(PyObject_Convert(vo, V)) ? PYERR() : takeresult(V)
    else
        ism1(PyObject_Convert(vo, Any)) ? PYERR() : takeresult(Any)
    end

pyjlany_setitem(xo::PyPtr, ko::PyPtr, vo::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)
    k = pyjl_getindices(x, ko)
    k === PYERR() && return Cint(-1)
    try
        if isnull(vo)
            delete!(x, k...)
            Cint(0)
        else
            v = pyjl_getvalue(x, vo)
            v === PYERR() && return Cint(-1)
            x[k...] = v
            Cint(0)
        end
    catch err
        if err isa BoundsError && err.a === x
            PyErr_SetStringFromJuliaError(PyExc_IndexError(), err)
        elseif err isa KeyError && (err.key === k || (err.key,) === k)
            PyErr_SetStringFromJuliaError(PyExc_KeyError(), err)
        elseif err isa MethodError && err.f === delete!
            PyErr_SetStringFromJuliaError(PyExc_TypeError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        Cint(-1)
    end
end

pyjlany_iter(xo::PyPtr) = PyJuliaIteratorValue_New(Iterator(PyJuliaValue_GetValue(xo)))

pyjlany_contains(xo::PyPtr, vo::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)
    r = PyObject_TryConvert(vo, eltype(x))
    r == -1 && return Cint(-1)
    r == 0 && return Cint(0)
    v = takeresult(eltype(x))
    try
        Cint(v in x)
    catch err
        if err isa MethodError && err.f === :in
            PyErr_SetStringFromJuliaError(PyExc_TypeError(), err)
        else
            PyErr_SetJuliaError(err)
        end
        Cint(-1)
    end
end

pyjlany_richcompare(xo::PyPtr, yo::PyPtr, op::Cint) = begin
    x = PyJuliaValue_GetValue(xo)
    r = PyObject_TryConvert(yo, Any)
    r == -1 && return PyNULL
    r == 0 && return PyNotImplemented_New()
    y = takeresult()
    try
        if op == Py_EQ
            PyObject_From(x == y)
        elseif op == Py_NE
            PyObject_From(x != y)
        elseif op == Py_LE
            PyObject_From(x <= y)
        elseif op == Py_LT
            PyObject_From(x < y)
        elseif op == Py_GE
            PyObject_From(x >= y)
        elseif op == Py_GT
            PyObject_From(x > y)
        else
            PyErr_SetString(PyExc_ValueError(), "bad op given to richcompare: $op")
            PyNULL
        end
    catch err
        if err isa MethodError && err.f in (==, !=, <=, <, >=, >)
            PyNotImplemented_New()
        else
            PyErr_SetJuliaError(err)
            PyNULL
        end
    end
end

const ALL_MIMES = [
    "text/plain",
    "text/html",
    "text/markdown",
    "text/json",
    "text/latex",
    "text/xml",
    "text/csv",
    "application/javascript",
    "application/pdf",
    "application/ogg",
    "image/jpeg",
    "image/png",
    "image/svg+xml",
    "image/gif",
    "image/webp",
    "image/tiff",
    "image/bmp",
    "audio/aac",
    "audio/mpeg",
    "audio/ogg",
    "audio/opus",
    "audio/webm",
    "audio/wav",
    "audio/midi",
    "audio/x-midi",
    "video/mpeg",
    "video/ogg",
    "video/webm",
]

pyjlany_repr_mimebundle(xo::PyPtr, args::PyPtr, kwargs::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)
    ism1(PyArg_CheckNumArgsEq("_repr_mimebundle_", args, 0)) && return PyNULL
    ism1(PyArg_GetArg(Union{Set{String},Nothing}, "_repr_mimebundle_", kwargs, "include", nothing)) && return PyNULL
    inc = takeresult(Union{Set{String},Nothing})
    ism1(PyArg_GetArg(Union{Set{String},Nothing}, "_repr_mimebundle_", kwargs, "exclude", nothing)) && return PyNULL
    exc = takeresult(Union{Set{String},Nothing})
    # decide which mimes to include
    if inc === nothing
        # default set of mimes to try
        mimes = copy(ALL_MIMES)
        # looks for mimes on show methods for the type
        for meth in methods(show, Tuple{IO, MIME, typeof(x)}).ms
            mimetype = meth.sig.parameters[3]
            mimetype isa DataType || continue
            mime = string(mimetype.parameters[1])
            push!(mimes, mime)
        end
    else
        mimes = push!(inc, "text/plain")
    end
    exc === nothing || setdiff!(mimes, exc)
    # make the bundle
    bundle = PyDict_New()
    isnull(bundle) && return PyNULL
    for m in mimes
        try
            io = IOBuffer()
            show(io, MIME(m), x)
            v = take!(io)
            mo = PyUnicode_From(m)
            isnull(mo) && (Py_DecRef(bundle); return PyNULL)
            vo = istextmime(m) ? PyUnicode_From(v) : PyBytes_From(v)
            isnull(vo) && (Py_DecRef(mo); Py_DecRef(bundle); return PyNULL)
            err = PyDict_SetItem(bundle, mo, vo)
            Py_DecRef(mo)
            Py_DecRef(vo)
            ism1(err) && (Py_DecRef(bundle); return PyNULL)
        catch err
            # silently skip anything that didn't work
        end
    end
    bundle
end

pyjlany_name(xo::PyPtr, ::Ptr{Cvoid}) =
    try
        PyObject_From(string(nameof(PyJuliaValue_GetValue(xo))))
    catch err
        if err isa MethodError && err.f === nameof
            PyErr_SetString(PyExc_AttributeError(), "__name__")
        else
            PyErr_SetJuliaError(err)
        end
        PyNULL
    end

pyjlany_toraw(xo::PyPtr, ::PyPtr) = PyJuliaRawValue_New(PyJuliaValue_GetValue(xo))

struct ExtraNewline{T}
    value :: T
end
Base.show(io::IO, m::MIME, x::ExtraNewline) = show(io, m, x.value)
Base.show(io::IO, m::MIME"text/plain", x::ExtraNewline) = (show(io, m, x.value); println(io))
Base.showable(m, x::ExtraNewline) = showable(m, x.value)

pyjlany_display(xo::PyPtr, ::PyPtr) = try
    x = PyJuliaValue_GetValue(xo)
    display(ExtraNewline(x))
    PyNone_New()
catch err
    PyErr_SetJuliaError(err)
    PyNULL
end

pyjlany_help(xo::PyPtr, ::PyPtr) =
    try
        x = Docs.doc(PyJuliaValue_GetValue(xo))
        display(ExtraNewline(x))
        PyNone_New()
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlany_positive(xo::PyPtr) =
    try
        PyObject_From(+(PyJuliaValue_GetValue(xo)))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlany_negative(xo::PyPtr) =
    try
        PyObject_From(-(PyJuliaValue_GetValue(xo)))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlany_absolute(xo::PyPtr) =
    try
        PyObject_From(abs(PyJuliaValue_GetValue(xo)))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

struct pyjlany_binop{F}
    f::F
end
(f::pyjlany_binop)(xo::PyPtr, yo::PyPtr) = begin
    PyJuliaValue_Check(xo) || return PyNotImplemented_New()
    PyJuliaValue_Check(yo) || return PyNotImplemented_New()
    x = PyJuliaValue_GetValue(xo)
    y = PyJuliaValue_GetValue(yo)
    try
        PyObject_From(f.f(x, y))
    catch err
        if err isa MethodError && err.f === f.f
            PyNotImplemented_New()
        else
            PyErr_SetJuliaError(err)
            PyNULL
        end
    end
end

pyjlany_power(xo::PyPtr, yo::PyPtr, zo::PyPtr) = begin
    PyJuliaValue_Check(xo) || return PyNotImplemented_New()
    PyJuliaValue_Check(yo) || return PyNotImplemented_New()
    x = PyJuliaValue_GetValue(xo)
    y = PyJuliaValue_GetValue(yo)
    if PyNone_Check(zo)
        try
            PyObject_From(x^y)
        catch err
            if err isa MethodError && err.f === ^
                PyNotImplemented_New()
            else
                PyErr_SetJuliaError(err)
                PyNULL
            end
        end
    else
        PyJuliaValue_Check(zo) || return PyNotImplemented_New()
        z = PyJuliaValue_GetValue(zo)
        try
            PyObject_From(powermod(x, y, z))
        catch err
            if err isa MethodError && err.f === powermod
                PyNotImplemented_New()
            else
                PyErr_SetJuliaError(err)
                PyNULL
            end
        end
    end
end