text:
using Flux

cd(@__DIR__)

include("pointcnn_cls.jl")

model = Classifier(40)|>gpu

include("data.jl")

file1 = h5open(train_files[1], "r")

batch_size = 32
train_data_len = 2048

# batches : Array of minibatches of size (2048, 1, 3, batch_size) = (n_p, 1, dimentions, batch_size)
# y : array of (40, batch_size) multidimentional arrays

batches = [
    cat(
        [reshape(transpose(file1["data"][:,:,i] |> gpu), (2048, 1, 3)) for i=j:j+batch_size-1]...,
        dims=4
    ) for j=1:batch_size:Int(floor(train_data_len / batch_size))*batch_size
];

y = [
    Flux.onehotbatch(file1["label"][1,j:j+batch_size-1], 0:39)|>gpu for j=1:batch_size:Int(floor(train_data_len / batch_size))*batch_size
]

# optimizer
opt = Flux.ADAM()

function loss(p, y)
    logpred = model(p)
    Flux.logitcrossentropy(logpred, y)
end

Flux.train!(loss, Flux.params(model), zip(batches, y), opt)


text:
function test_GaussianRewards()
    means = [x -> x[1], x -> 2 * x[1]]
    rdist = GaussianRewards(means, 2.0)
    seed!(123); r1 = sample([1.0], 2, rdist)
    seed!(123); r2 = rand(Normal(2.0, 2.0))
    @test r1 == r2
end

test_GaussianRewards()

ext: jl


text:
# # Gross-Pitaevskii equation with magnetism

# We solve the 2D Gross-Pitaevskii equation with a magnetic field.
# This is similar to the
# previous example ([Gross-Pitaevskii equation in one dimension](@ref)),
# but with an extra term for the magnetic field.

using DFTK
using StaticArrays
using Plots

# Unit cell. Having one of the lattice vectors as zero means a 2D system
a = 10
lattice = a .* [[1 0 0.]; [0 1 0]; [0 0 0]];

# Confining scalar potential, and magnetic vector potential
pot(x, y, z) = (x - a/2)^2 + (y - a/2)^2
Apot(x, y, z) = .2 * @SVector [y - a/2, -(x - a/2), 0]
Apot(X) = Apot(X...);

# Parameters
Ecut = 20  # Increase this for production
C = 500.0
α = 2
n_electrons = 1;  # Increase this for fun

# Collect all the terms, build and run the model
terms = [Kinetic(),
         ExternalFromReal(X -> pot(X...)),
         PowerNonlinearity(C, α),
         Magnetic(Apot),
]
model = Model(lattice; n_electrons=n_electrons,
              terms=terms, spin_polarization=:spinless)  # "spinless electrons"
basis = PlaneWaveBasis(model, Ecut)
scfres = direct_minimization(basis, tol=1e-5)  # Reduce tol for production
heatmap(scfres.ρ.real[:, :, 1], c=:blues)


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

@doc raw"""IngressRule represents the rules mapping the paths under a specified host to the related backend services. Incoming requests are first evaluated for a host match, then routed to the backend associated with the matching IngressRuleValue.

    IoK8sApiExtensionsV1beta1IngressRule(;
        host=nothing,
        http=nothing,
    )

    - host::String : Host is the fully qualified domain name of a network host, as defined by RFC 3986. Note the following deviations from the \"host\" part of the URI as defined in the RFC: 1. IPs are not allowed. Currently an IngressRuleValue can only apply to the    IP in the Spec of the parent Ingress. 2. The `:` delimiter is not respected because ports are not allowed.    Currently the port of an Ingress is implicitly :80 for http and    :443 for https. Both these may change in the future. Incoming requests are matched against the host before the IngressRuleValue. If the host is unspecified, the Ingress routes all traffic based on the specified IngressRuleValue.
    - http::IoK8sApiExtensionsV1beta1HTTPIngressRuleValue
"""
mutable struct IoK8sApiExtensionsV1beta1IngressRule <: SwaggerModel
    host::Any # spec type: Union{ Nothing, String } # spec name: host
    http::Any # spec type: Union{ Nothing, IoK8sApiExtensionsV1beta1HTTPIngressRuleValue } # spec name: http

    function IoK8sApiExtensionsV1beta1IngressRule(;host=nothing, http=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1IngressRule, Symbol("host"), host)
        setfield!(o, Symbol("host"), host)
        validate_property(IoK8sApiExtensionsV1beta1IngressRule, Symbol("http"), http)
        setfield!(o, Symbol("http"), http)
        o
    end
end # type IoK8sApiExtensionsV1beta1IngressRule

const _property_map_IoK8sApiExtensionsV1beta1IngressRule = Dict{Symbol,Symbol}(Symbol("host")=>Symbol("host"), Symbol("http")=>Symbol("http"))
const _property_types_IoK8sApiExtensionsV1beta1IngressRule = Dict{Symbol,String}(Symbol("host")=>"String", Symbol("http")=>"IoK8sApiExtensionsV1beta1HTTPIngressRuleValue")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1IngressRule }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1IngressRule))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1IngressRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1IngressRule[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1IngressRule }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1IngressRule[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1IngressRule)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1IngressRule }, name::Symbol, val)
end


text:
function add_constraint(obj::AbstractMultipleShooting, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "addConstraint", void, (jint, jint, jdouble), arg0, arg1, arg2)
end

function compute(obj::AbstractMultipleShooting)
    return jcall(obj, "compute", List, ())
end

function set_closed_orbit_constraint(obj::AbstractMultipleShooting, arg0::jboolean)
    return jcall(obj, "setClosedOrbitConstraint", void, (jboolean,), arg0)
end

function set_epoch_freedom(obj::AbstractMultipleShooting, arg0::jint, arg1::jboolean)
    return jcall(obj, "setEpochFreedom", void, (jint, jboolean), arg0, arg1)
end

function set_patch_point_component_freedom(obj::AbstractMultipleShooting, arg0::jint, arg1::jint, arg2::jboolean)
    return jcall(obj, "setPatchPointComponentFreedom", void, (jint, jint, jboolean), arg0, arg1, arg2)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end


text:
cd(@__DIR__)
using Pkg
Pkg.activate("..")
Pkg.instantiate()

##
using LinearAlgebra
using StatsBase
using CSV
using DataFrames
using Pipe
using ProgressMeter
using Random

##

using MCPhylo
Random.seed!(4928370335238343681)

##

wals = CSV.read("../../data/charMtx.csv", DataFrame)
d = CSV.read("../../data/fpairMtx.csv", DataFrame)

famFreqs = combine(groupby(wals, :glot_fam), nrow)
sort!(famFreqs, :nrow, rev=true)

families = filter(x -> x.nrow>1, famFreqs).glot_fam
isolates = filter(x -> x.nrow==1, famFreqs).glot_fam

isoDict = @pipe wals |>
      filter(x -> x.glot_fam ∈ isolates, _) |>
      zip(_.glot_fam, _.longname) |>
      Dict
##

fm2trees = Dict()
@showprogress for fm in families
      fm2trees[fm] = MCPhylo.ParseNewick("../../data/posteriorTrees/$fm.posterior.tree")
end

for fm in isolates
      fm2trees[fm] = repeat([Node(isoDict[fm])], 1000)
end

##

function renumberTree(tree::GeneralNode)
      tips = [nd for nd in post_order(tree) if nd.nchild==0]
      nonTips = [nd for nd in post_order(tree) if nd.nchild > 0]
      sort!(tips, lt= (x,y) -> x.name < y.name)
      for (i,nd) in enumerate(tips)
            nd.num = i
      end
      for (i,nd) in enumerate(nonTips)
            nd.num = i+length(tips)
      end
      tree
end

##

for fm in keys(fm2trees)
      fm2trees[fm] = renumberTree.(fm2trees[fm])
end

##

lineages = vcat(families, isolates)

##

ttrees = [[fm2trees[fm][i] for fm in lineages] for i in 1:1000]

##

states = ["a", "b", "c", "d"]

nsites = 1

nbase = 4

rates = [1.0]

##
taxa = Vector{String}[]
for t in ttrees[1]
      push!(taxa, sort([nd.name for nd in pre_order(t) if nd.nchild==0]))
end

##

max_stars_repo_name: erathorn/phylogeneticTypology

max_issues_repo_name: erathorn/phylogeneticTypology

max_forks_repo_name: erathorn/phylogeneticTypology


text:
# FitzHugh_Nagumo_2D_compute_spiral_wave_PS_v2.jl
# ================================================   u.p. 4.9.21

# Computes and saves spiral wave solution of the Fitzhugh-Nagumo model.
# Used for figures 11.4 and 11.6, but also in the animation spiralwave in animations/11.

using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, OrdinaryDiffEq
using Random
using Statistics

using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

function FitzHugh_Nagumo_ODE(du,u,p,t)
    a, b, d, ε, hsq6 = p
    uu = @view u[1,:,:]
    vv = @view u[2,:,:]

    # Create padded u matrix to incorporate Newman boundary conditions - 9-point stencil
    uuu=[ [uu[2,2] uu[2,:]' uu[2,end-1] ] ; [ uu[:,2] uu uu[:,end-1] ] ; [ uu[end-1,2] uu[end-1,:]' uu[end-1,end-1] ] ]
    diff_term = d .* (
                4 .* (uuu[2:end-1,1:end-2] .+ uuu[2:end-1,3:end] .+
                        uuu[3:end,2:end-1] .+ uuu[1:end-2,2:end-1] ) .+
                        uuu[3:end,3:end] .+ uuu[3:end,1:end-2] .+ uuu[1:end-2,3:end] .+
                        uuu[1:end-2,1:end-2] .- 20 .*uu  )  ./ hsq6

    du[1,:,:] = @. a*uu*(1 - uu)*(uu - b) - vv + diff_term
    du[2,:,:] = @. ε*(uu - vv)
end

L = 300.0  #  size of domain
N = 300   #  no. of grid points
h = L/N # spatial steps size
hsq6 = 6*h*h

# parameters
a = 3.
b = 0.2
d = 1.    #  diffusion
ε = 0.01

#  initial conditions
uu = 0.5 .* ones(N,N)
vv = 0.5 .* ones(N,N)
# This initial condition sets the spiral wave
uu[1:end, 1:(2N÷3)] .= 0
vv[1:(2N÷3 - N÷12), 1:end] .= 0

u0 = zeros(2,N,N)
u0[1,:,:] = uu
u0[2,:,:] = vv

# solve ODE system
p = a, b, d, ε, hsq6

saveat = 0.0:600.0
tspan = (0.0, saveat[end])

prob = ODEProblem(FitzHugh_Nagumo_ODE, u0, tspan, p)
@time sol = solve(prob, Tsit5();
    reltol=1e-6, abstol=1e-6, maxiters = Inf, saveat = saveat,
    progress = true, progress_steps = 10,
)
tvec = sol.t
uout = sol.u

params = @ntuple a b d ε hsq6 N L
data = @strdict tvec uout params

wsave(datadir("FitzHugh", "spiralwave.jld2"), data)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# Default implementations of ONNX operators

import NNlib
import Statistics: mean
import StaticArrays: SVector

using LinearAlgebra

flipweights(w::AbstractArray{T,N}) where {T,N} = w[(size(w, i):-1:1 for i = 1:(N-2))..., :, :]

conv(x, w; kw...) = NNlib.conv(x, flipweights(w); kw...)

function conv(x, w, b; kw...)
    d = ndims(x) - 2
    bias_size = (ntuple(_ -> 1, d)..., :, 1)
    b = reshape(b, bias_size)
    return conv(x, w; kw...) .+ b
end

function onnx_gemm(A, B, C; tA = 0, tB = false, α = 1, β = 1)
    A = Bool(tA) ? A' : A
    B = Bool(tB) ? B' : B
    # note: order of arguments reversed due to row-major layout
    return α * B * A .+ β * C
end

function onnx_gemm(A, B; tA = 0, tB = 0, α = 1)
    A = Bool(tA) ? A' : A
    B = Bool(tB) ? B' : B
    # note: order of arguments reversed due to row-major layout
    return α * B * A
end

# Julia-friendly flatten
function flatten(x; dim = ndims(x) - 1)
    sz = size(x)
    keep = dim < ndims(x) ? sz[dim+1:end] : 1
    return reshape(x, :, keep...)
end

# ONNX-specific flatten
function onnx_flatten(x; axis = 1)
    dim = axis >= 0 ? ndims(x) - axis + 1 : axis + 1
    return flatten(x; dim = dim)
end

add(xs...) = .+(xs...)
mul(xs...) = .*(xs...)
relu(x) = NNlib.relu.(x)
maxpool(x; kernel, pad = 0, stride = 1) = NNlib.maxpool(x, kernel; pad = pad, stride = stride)

# common functional implementation for batch and instance normalization based on
# https://github.com/FluxML/Flux.jl/blob/06970a5fbbb1cb485c5d2cba597a78fb453fc713/src/layers/normalise.jl#L166-L197
function normalization(x::AbstractArray{T,N}, γ, β, μ, σ², reduce_dims, affine_shape;
        ϵ=1f-5, mtm=0.1f0, training=false) where {T, N}
    # init variables in the function scope instead of the if's scope
    μnext = μ
    σ²next = σ²
    if !training  # testmode
        stats_shape = ntuple(i -> i == N - 1 ? size(x, N - 1) : 1, N)
        μ = reshape(μ, stats_shape)
        σ² = reshape(σ², stats_shape)
    else  # trainmode or testmode without tracked stats
        μold = μ
        σ²old = σ²
        μ = mean(x; dims = reduce_dims)
        σ² = mean((x .- μ) .^ 2; dims = reduce_dims)
        m = prod(size(x, i) for i in reduce_dims)  # needed for computing corrected var
        μnew = vec(N ∈ reduce_dims ? μ : mean(μ, dims = N))
        σ²new = vec(N ∈ reduce_dims ? σ² : mean(σ², dims = N))
        μnext = (1 - mtm) .* μold .+ mtm .* μnew
        σ²next = (1 - mtm) .* σ²old .+ mtm .* (m / (m - one(eltype(σ²old)))) .* σ²new
    end
    out = (x .- μ) ./ sqrt.(σ² .+ ϵ)
    if !isnothing(γ) && !isnothing(β)
        γ = reshape(γ, affine_shape)
        β = reshape(β, affine_shape)
        out = γ .* out .+ β
    end
    if training
        return out, μnext, σ²next
    else
        return out
    end
end

function batch_norm(x::AbstractArray{T,N}, γ, β, μ, σ²;
        ϵ=1f-5, mtm=0.1f0, training=false) where {T,N}
    reduce_dims = [1:N-2; N]
    affine_shape = ntuple(i -> i == N-1 ? size(x, N-1) : 1, N)
    return normalization(x, γ, β, μ, σ², reduce_dims, affine_shape;
        ϵ=ϵ, mtm=mtm, training=training)
end

function instance_norm(x::AbstractArray{T,N}, γ, β, μ, σ²;
        ϵ=1f-5, mtm=0.1f0, training=false) where {T,N}
    reduce_dims = 1:N-2
    affine_shape = ntuple(i -> i == N-1 ? size(x, N-1) : 1, N)
    return normalization(x, γ, β, μ, σ², reduce_dims, affine_shape;
        ϵ=ϵ, mtm=mtm, training=training)
end

# implementation from
# https://github.com/FluxML/Flux.jl/blob/f66be896d3d2698ce77ce8b7788b4317285bf0b2/src/layers/conv.jl#L605-L614
function global_average_pool(x)
    # Input size
    x_size = size(x)
    # Kernel size
    k = x_size[1:end-2]
    # Pooling dimensions
    pdims = NNlib.PoolDims(x, k)
    return NNlib.meanpool(x, pdims)
end

size_vector(x) = SVector(size(x))

"""
    take(data, idxs; dim=ndims(data))

Take elements from an array along an axis. For example, for a 4D data
and dim=3, it is roughly equivalent to `data[:, :, idxs, :]`, but allows
multidimensional idxs. See `numpy.take` for a more detailed explanation
of the concept.

In the context of ONNX, `take` is used to implement Gather operation.
We do NOT record this function directly to the tape during loading though,
but instead use a more ONNX-friendly wrapper `onnx_gather()`.

Note: in ONNX, Gather is different from GatherElements, GatherND and
Julia's `NNlib.gather()`.
"""
function take(
        data::AbstractArray{T, N}, idxs::AbstractArray{Int, M};
        dim=ndims(data)) where {T, N, M}
    if length(idxs) == 1 && data isa SVector
        # we use SVector to represent array size, Gather(arr_sz, idx)
        # works as size(arr, idx); but since dimensions are reversed,
        # we need to reverse the index as well
        # see https://github.com/FluxML/ONNX.jl/issues/62 for details
        return data[length(data) .- idxs .+ 1]
    end
    if length(idxs) == 1
        # special case, works as getindex
        return data[idxs]
    end
    # we will take slices of data of this size
    size_before = (size(data)[1:dim-1]...,)
    size_after = (size(data)[dim+1:ndims(data)]...,)
    # and put them into output array at out[:, :, ..., idxs[i, j, ...]]
    out = similar(data, (size_before..., size(idxs)..., size_after...))
    colons_before = [(:) for _=1:dim-1]
    colons_after = [(:) for _=dim+1:ndims(data)]
    # iteration over idxs doesn't depend on data or dimension
    # we iterate over the last index purely due to memory layout
    for i=1:size(idxs, ndims(idxs))
        # R - slice of idxs (not slice of data!)
        R = [[(:) for _=1:ndims(idxs)-1]..., i]
        # ensure I = idxs[R...] is itself an array and not a scalar
        I = [idxs[R...]...,]
        slice = data[colons_before..., I, colons_after...]
        out[colons_before..., R..., colons_after...] = slice
    end
    return out
end

take(data::AbstractArray, idxs::Integer; dim=ndims(data)) =
    take(data, [idxs]; dim=dim)

"""
    onnx_gather(data::AbstractArray, idxs::AbstractArray{Int}; dim=ndims(data))

Implemntation of ONNX's Gather operation with 0-based indices.
For a Julia-friendly version, see `take`.
"""
function onnx_gather(
        data::AbstractArray{T, N}, idxs::AbstractArray{Int, M};
        dim=ndims(data)) where {T, N, M}
    @assert all(idxs .>= 0) "Gather on negative indices is not implemented yet"
    idxs_adjusted = idxs .+ 1
    return take(data, idxs_adjusted; dim=dim)
end

# julia-friendly
function NNlib.unsqueeze(x::AbstractArray, dims)
    new_shape = collect(size(x))
    for d in sort(collect(dims))
        insert!(new_shape, d, 1)
    end
    return reshape(x, new_shape...)
end

# ONNX-friendly, e.g. axes is 0-based, row-major
function onnx_unsqueeze(x::AbstractArray, axes::Vector)
    # ndims(data) + length(axes) => size of the array after unsqueezing
    # .- axes                    => to reverse dimensions
    # .+ 1                       => to convert to 1-based indexing
    # .- 1                       => correction by 1
    dims = ndims(x) + length(axes) .- axes
    return NNlib.unsqueeze(x, dims)
end

function onnx_slice(
        data::AbstractArray,
        starts::VecOrMat{Int}, ends::VecOrMat{Int},
        axes::Vector{Int}=Int[], steps::Vector{Int}=Int[])
    axes = isempty(axes) ? collect(0:ndims(data)-1) : axes
    steps = isempty(steps) ? [1 for i=1:ndims(data)] : steps
    @assert all(starts .>= 0) "Negative indices are not supported yet"
    @assert all(ends .>= 0) "Negative indices are not supported yet"
    # construct ranges, adjusting starts to 1-based indexing
    ranges = [s+1 : st : e for (s, st, e) in zip(starts, steps, ends)]
    # reversed, 1-based dimensions
    dims = ndims(data) .- axes
    # dimension => range mapping
    d2r = Dict(zip(dims, ranges))
    I = [get(d2r, i, (:)) for i=1:ndims(data)]
    return data[I...]
end


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
import Aspasia

import GLPK
oa_solver = MOI.OptimizerWithAttributes(
    GLPK.Optimizer,
    "msg_lev" => 0,
    "tol_int" => 1e-9,
    "tol_bnd" => 1e-9,
    "mip_gap" => 0.0,
)

function _run_moi_tests(use_iterative_method::Bool, oa_solver)
    config = MOI.Test.Config(
        atol = 1e-3,
        rtol = 1e-3,
        exclude = Any[
            MOI.ConstraintDual,
            MOI.ConstraintBasisStatus,
            MOI.DualObjectiveValue,
        ],
    )

    aspasia = Aspasia.Optimizer()
    MOI.set(aspasia, MOI.Silent(), true)
    MOI.set(aspasia, MOI.RawOptimizerAttribute("oa_solver"), oa_solver)
    MOI.set(
        aspasia,
        MOI.RawOptimizerAttribute("use_iterative_method"),
        use_iterative_method,
    )
    opt = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        MOI.Bridges.full_bridge_optimizer(aspasia, Float64),
    )

    excludes = String[
        # not implemented:
        "test_attribute_SolverVersion",
        # invalid model:
        "test_constraint_ZeroOne_bounds_3",
        # CachingOptimizer does not throw if optimizer not attached:
        "test_model_copy_to_UnsupportedAttribute",
        "test_model_copy_to_UnsupportedConstraint",
        # other:
        "test_objective_qp",
        "test_quadratic",
    ]
    if use_iterative_method
        push!(excludes, "test_unbounded")
    end
    # excludes = String[]

    # includes = String["test_conic_SecondOrderCone", "test_conic_PositiveSemidefiniteCone"]
    includes = String[]

    return MOI.Test.runtests(opt, config, include = includes, exclude = excludes)
end

println("starting MOI tests")
@testset "MOI tests" begin
    println("starting iterative method tests")
    @testset "iterative method" _run_moi_tests(true, oa_solver)
    println("starting OA solver driven method tests")
    @testset "OA solver driven method" _run_moi_tests(false, oa_solver)
    println("finished MOI tests")
end


text:
using PyPlot
pygui(true)

M = 3
f_s = 11
d = floor(Int, (f_s-1)/2)
N = f_s*M
center_idx = floor(Int, (N+1)/2)
t = Array(range(-M/2, M/2, length=N))
δt = 3
s = zeros(Int, N)
s[center_idx-d:center_idx+d] .= 1
zp = circshift(s, δt)
ze = circshift(zp, -d)
zl = circshift(zp, d)
S = float.(deepcopy(s))
S[center_idx-2*d:center_idx-1] .= Array(0:2*d-1)/(2*d)
S[center_idx+1:center_idx+2*d] .= reverse(S[center_idx-2*d:center_idx-1])

fig = figure(figsize=(7.5, 3))
ax1 = subplot2grid((4,2), (0,0), rowspan=1, colspan=1)
ax1.step(t, ze, "b-")
ylim([0, 1.2])
ax1.axis("off")
ax2 = subplot2grid((4,2), (1,0), rowspan=1, colspan=1)
ax2.step(t, zp, color="grey", linestyle="-")
ylim([0, 1.2])
ax2.axis("off")
ax3 = subplot2grid((4,2), (2,0), rowspan=1, colspan=1)
ax3.step(t, zl, "r-")
ylim([0, 1.2])
ax3.axis("off")
ax4 = subplot2grid((4,2), (3,0), rowspan=1, colspan=1)
ax4.step(t, s, "k-")
ax4.arrow(t[1], 0, t[end]-t[1], 0., fc="k", ec="k", lw = 1.,
             head_width=8/80*(t[end]-t[1]), head_length=1.75/20*1.2,
             overhang=0.3, length_includes_head=true, clip_on=false)
ax4.arrow(0, 0, 0, 11, fc="k", ec="k", lw = 1.,
             head_width=2/30*1.2, head_length=9/60*(t[end]-t[1]),
             overhang=0.3, length_includes_head=true, clip_on=false)
ylim([0, 1.2])
ax4.axis("off")

ax5 = subplot2grid((4,2), (0,1), rowspan=4, colspan=1)
ax5.plot(t, S, "k-", label=L"s[n]")
ax5.vlines(x=t[center_idx+δt-d], ymin=0, ymax=S[center_idx+δt-d],
           color="grey", linestyle=":")
ax5.vlines(x=t[center_idx+δt], ymin=0, ymax=S[center_idx+δt],
           color="grey", linestyle=":")
ax5.vlines(x=t[center_idx+δt+d], ymin=0, ymax=S[center_idx+δt+d],
           color="grey", linestyle=":")
ax5.scatter(t[center_idx+δt-d], S[center_idx+δt-d], c="b", label=L"Z_E")
ax5.scatter(t[center_idx+δt], S[center_idx+δt], c="grey", label=L"Z_P")
ax5.scatter(t[center_idx+δt+d], S[center_idx+δt+d], c="r", label=L"Z_L")
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax5.spines["top"].set_visible(false)
ax5.spines["right"].set_visible(false)
ax5.spines["bottom"].set_visible(false)
ax5.spines["left"].set_visible(false)
# ax5.spines["left"].set_position("center")
ax5.arrow(t[1], 0, t[end]-t[1], 0., fc="k", ec="k", lw = 1.,
             head_width=1/80*(t[end]-t[1]), head_length=1.5/20*1.2,
             overhang=0.3, length_includes_head=true, clip_on=false)
ax5.arrow(0, 0, 0, 1.285, fc="k", ec="k", lw = 1.,
             head_width=2/30*1.2, head_length=1/60*(t[end]-t[1]),
             overhang=0.3, length_includes_head=true, clip_on=false)
# xlim(t[1], t[end])
ylim([0, 1.2])
legend()

subplots_adjust(wspace=0.075, hspace=1.5, left=0.01, right=0.99)
savefig("figures/ch3_correlators_for_figure.svg", dpi=300)


text:
using LinearAlgebra

# Not sure if this is implemented correctly or is just sensitive to meta learning rate.
mutable struct TIDBDOptimizer{T, TF} <: AbstractAccumulatingTraceOptimizer where {T<:Real,TF<:AbstractFuncApprox}
    f::TF
    η::T         # meta stepsize
    α::Array{T}  # learning rates, i.e., exp.(β)
    β::Array{T}  # ln learning rates
    γ::T         # discount factor
    λ::T         # eligibility trace decay
    e::Array{T}  # eligibility trace vector
    h::Array{T}  # memory vector for stepsizes
    g::Array{T}  # gradient place holder

    function TIDBDOptimizer(fun::AbstractFuncApprox, b::T, η::T, γ::T, λ::T) where {T<:Real}
        tp = eltype(get_params(fun))
        s = size(get_params(fun))

        β = zeros(tp, s)
        h = zeros(tp, s)
        fill!(β, b)
        α = exp.(β)
        new{tp, typeof(fun)}(fun, η, α, β, tp(γ), tp(λ), zeros(tp, s), h, zeros(tp, s))
    end
end

function update!(opt::TIDBDOptimizer{T}, δ::T, ϕ::Array{T}) where {T<:Real}
    @. opt.β += opt.η * δ * ϕ * opt.h
    @. opt.α = exp(opt.β)
    @. opt.e = (opt.γ * opt.λ) * opt.e + ϕ

    @. opt.g = opt.e * (δ * opt.α)
    add_to_params!(opt.f, opt.g)
    @. opt.h = opt.h * clamp(1. - opt.α * ϕ * opt.e, 0., Inf) + opt.α * δ * opt.e
end

function clone(opt::TIDBDOptimizer)
    opt2 = TIDBDOptimizer(opt.f, 1.0, opt.η, opt.γ, opt.λ)
    opt2.e .= deepcopy(opt.e)
    opt2.β .= deepcopy(opt.β)
    opt2.α .= deepcopy(opt.α)
    opt2.h .= deepcopy(opt.h)
    return opt2
end

function clone(opt::TIDBDOptimizer, f::AbstractFuncApprox)
    opt2 = TIDBDOptimizer(f, 1.0, opt.η, opt.γ, opt.λ)
    opt2.e .= deepcopy(opt.e)
    opt2.β .= deepcopy(opt.β)
    opt2.α .= deepcopy(opt.α)
    opt2.h .= deepcopy(opt.h)
    return opt2
end

mutable struct AutoTIDBDOptimizer{T, TF} <: AbstractAccumulatingTraceOptimizer where {T<:Real,TF<:AbstractFuncApprox}
    f::TF
    η::T         # meta stepsize
    α::Array{T}  # learning rates, i.e., exp.(β)
    β::Array{T}  # ln learning rates
    γ::T         # discount factor
    λ::T         # eligibility trace decay
    e::Array{T}  # eligibility trace vector, z in the paper
    h::Array{T}  # memory vector for stepsizes
    g::Array{T}  # gradient place holder
    z::Array{T}  # max observed update trace, η in the paper
    τ::T         # large value to perform normalization of meta learning rate

    function AutoTIDBDOptimizer(fun::AbstractFuncApprox, b::T, η::T, τ::T, γ::T, λ::T) where {T<:Real}
        tp = eltype(get_params(fun))
        s = size(get_params(fun))

        β = zeros(tp, s)
        h = zeros(tp, s)
        z = zeros(tp, s)
        fill!(β, b)
        α = exp.(β)
        new{tp, typeof(fun)}(fun, η, α, β, tp(γ), tp(λ), zeros(tp, s), h, zeros(tp, s), z, τ)
    end
end

function update!(opt::AutoTIDBDOptimizer{T}, δ::T, ϕ::Array{T}, ϕ′::Array{T}) where {T<:Real}
    gpp = (opt.γ .* ϕ′ .- ϕ)
    p1 = δ .* gpp .* opt.h
    p2 = opt.z .- (1. / opt.τ) .* opt.α .* gpp .* opt.e .* (δ .* ϕ .* opt.h .- opt.z)
    opt.z .= max.(abs.(p1), p2)

    @. opt.β -= opt.η / opt.z * δ * gpp * opt.h #* ϕ * opt.h
    M = max.(exp.(opt.β) * dot(vec(gpp), vec(opt.e)), 1.)
    @. opt.β -= log.(M)
    @. opt.α = exp(opt.β)
    @. opt.e = (opt.γ * opt.λ) * opt.e + ϕ

    @. opt.g = opt.e * (δ * opt.α)
    add_to_params!(opt.f, opt.g)
    # @. opt.h = opt.h * clamp(1. - opt.α * ϕ * opt.e, 0., Inf) + opt.α * δ * opt.e
    @. opt.h = opt.h * clamp(1. + opt.α * gpp * opt.e, 0., Inf) + opt.α * δ * opt.e
end

function clone(opt::AutoTIDBDOptimizer)
    opt2 = AutoTIDBDOptimizer(opt.f, 1.0, opt.η, opt.τ, opt.γ, opt.λ)
    opt2.e .= deepcopy(opt.e)
    opt2.β .= deepcopy(opt.β)
    opt2.α .= deepcopy(opt.α)
    opt2.h .= deepcopy(opt.h)
    opt2.z .= deepcopy(opt.z)
    return opt2
end

function clone(opt::AutoTIDBDOptimizer, f::AbstractFuncApprox)
    opt2 = AutoTIDBDOptimizer(f, 1.0, opt.η, opt.τ, opt.γ, opt.λ)
    opt2.e .= deepcopy(opt.e)
    opt2.β .= deepcopy(opt.β)
    opt2.α .= deepcopy(opt.α)
    opt2.h .= deepcopy(opt.h)
    opt2.z .= deepcopy(opt.z)

    return opt2
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
"""
    compute_frames!(actions::Vector{AA}; last_frames=nothing) where AA <: AbstractAction

Set action.frames.frames to the computed frames for each action in actions.
"""
function compute_frames!(
    actions::Vector{AA};
    last_frames = nothing,
) where {AA<:AbstractAction}
    for action in actions
        if last_frames === nothing && get_frames(action) === nothing
            throw(ArgumentError("Frames need to be defined explicitly in the initial
            AbstractAction like Action/BackgroundAction or SubAction."))
        end
        if get_frames(action) === nothing
            set_frames!(action, last_frames)
        end
        last_frames = get_frames(action)
    end
end

"""
    get_current_setting()

Return the current setting of the current action
"""
function get_current_setting()
    action = CURRENT_ACTION[1]
    return action.current_setting
end


text:
module KeyedDistributions

using AutoHashEquals
using AxisKeys
using Distributions
using Distributions: GenericMvTDist
using LinearAlgebra: Symmetric
using PDMatsExtras: submat
using Random: AbstractRNG

export KeyedDistribution, KeyedSampleable
export KeyedMvNormal, KeyedGenericMvTDist, MvNormalLike, MvTLike
export axiskeys, distribution

for T in (:Distribution, :Sampleable)
    KeyedT = Symbol(:Keyed, T)
    @eval begin
        """
            $($KeyedT)(d<:$($T), keys::Tuple{Vararg{AbstractVector}})

        Stores `keys` for each variate alongside the `$($T)` `d`,
        supporting all of the common functions of a `$($T)`.
        Common functions that return an `AbstractArray`, such as `rand`,
        will return a `KeyedArray` with keys derived from the `$($T)`.

        The type of `keys` is restricted to be consistent with
        [AxisKeys.jl](https://github.com/mcabbott/AxisKeys.jl).
        The length of the `keys` tuple must be the number of dimensions, which is 1 for
        univariate and multivariate distributions, and 2 for matrix-variate distributions.
        The length of each key vector in must match the length along each dimension.

        !!! Note
            For distributions that can be marginalized exactly, the $($KeyedT)) can be
            marginalised via the indexing or lookup syntax just like `KeyedArray`s.
            i.e. One can use square or round brackets to retain certain indices or keys and
            marginalise out the others. For example for `D::KeyedMvNormal` over `:a, :b, :c`:
             - `D(:a)` or D(1) will marginalise out `:b, :c` and return a `KeyedMvNormal`
               over `:a`.
             - `D([:a, :b])` or `D[[1, 2]]` will marginalise out `:c` and return a
               `KeyedMvNormal` over `:a, :b`.
        """
        @auto_hash_equals struct $KeyedT{F<:VariateForm, S<:ValueSupport, D<:$T{F, S}} <: $T{F, S}
            d::D
            keys::Tuple{Vararg{AbstractVector}}

            function $KeyedT(d::$T{F, S}, keys::Tuple{Vararg{AbstractVector}}) where {F, S}
                key_lengths = map(length, keys)
                key_lengths == _size(d) || throw(ArgumentError(
                    "lengths of key vectors $key_lengths must match " *
                    "size of distribution $(_size(d))"
                ))

                return new{F, S, typeof(d)}(d, keys)
            end
        end

        """
            $($KeyedT)(d<:$($T), keys::AbstractVector)

        Constructor for [`$($KeyedT)`](@ref) with one dimension of variates.
        The elements of `keys` correspond to the variates of the distribution.
        """
        $KeyedT(d::$T{F, S}, keys::AbstractVector) where {F, S} = $KeyedT(d, (keys, ))

        # Allows marginalisation via lookup syntax, using getindex.
        (d::$KeyedT)(keys...) = d[first(map(AxisKeys.findindex, keys, axiskeys(d)))]
    end
end

_size(d) = (length(d),)
_size(d::Sampleable{<:Matrixvariate}) = size(d)

"""
    KeyedDistribution(d::Distribution)

Constructs a [`KeyedDistribution`](@ref) using the keys of the parameter that matches `size(d)`.
If the parameter has no keys, uses `1:n` for the length `n` of each dimension.
"""
KeyedDistribution(d::Distribution) = KeyedDistribution(d, _keys(d))

_keys(d::Distribution) = _keys(first(filter(x -> size(x) == size(d), params(d))))
_keys(x::KeyedArray) = axiskeys(x)
_keys(x) = map(Base.OneTo, size(x))

const KeyedDistOrSampleable = Union{KeyedDistribution, KeyedSampleable}

# MvNormal and MvT are particularly useful - these make it easier to dispatch
const KeyedMvNormal = KeyedDistribution{Multivariate, Continuous, <:MvNormal}
const MvNormalLike = Union{MvNormal, KeyedMvNormal}

const KeyedGenericMvTDist = KeyedDistribution{Multivariate, Continuous, <:GenericMvTDist}
const MvTLike = Union{GenericMvTDist, KeyedGenericMvTDist}

# Use submat to preserve the covariance matrix PDMat type
function Base.getindex(d::KeyedMvNormal, i::Vector)::KeyedMvNormal
    return KeyedDistribution(MvNormal(d.d.μ[i], submat(d.d.Σ, i)), axiskeys(d)[1][i])
end

function Base.getindex(d::KeyedMvNormal, i::Integer)::KeyedDistribution
    return KeyedDistribution(Normal(d.d.μ[i], d.d.Σ[i, i]), [axiskeys(d)[1][i]])
end

function Base.getindex(d::KeyedGenericMvTDist, i::Vector)::KeyedGenericMvTDist
    return KeyedDistribution(MvTDist(d.d.df, d.d.μ[i], submat(d.d.Σ, i)), axiskeys(d)[1][i])
end

# Access methods

"""
    distribution(::KeyedDistribution) -> Distribution
    distribution(::KeyedSampleable{F, S, D}) -> D

Return the wrapped distribution.
"""
distribution(d::KeyedDistOrSampleable) = d.d

# AxisKeys functionality

"""
    axiskeys(d::Union{KeyedDistribution, KeyedSampleable})

Return the keys for the variates of the `KeyedDistribution` or `KeyedSampleable`.
"""
AxisKeys.axiskeys(d::KeyedDistOrSampleable) = d.keys

# Standard functions to overload for new Distribution and/or Sampleable
# https://juliastats.org/Distributions.jl/latest/extends/#Create-New-Samplers-and-Distributions

function Distributions._rand!(
    rng::AbstractRNG,
    d::KeyedDistOrSampleable,
    x::AbstractVector{<:Real}
)
    sample = Distributions._rand!(rng, distribution(d), x)
    return KeyedArray(sample, axiskeys(d))
end

Base.length(d::KeyedDistOrSampleable) = length(distribution(d))

Distributions.size(d::KeyedDistribution{<:Matrixvariate}) = size(distribution(d))

Distributions.sampler(d::KeyedDistribution) = sampler(distribution(d))

Distributions.params(d::KeyedDistOrSampleable) = params(distribution(d))

Base.eltype(d::KeyedDistribution) = eltype(distribution(d))

function Distributions._logpdf(d::KeyedDistribution, x::AbstractArray)
    # Workaround when KeyedArray is parameter of Distribution
    # https://github.com/mcabbott/AxisKeys.jl/issues/54
    dist = distribution(d)
    T = typeof(dist)
    args = map(_maybe_parent, params(dist))
    unkeyed_dist = T.name.wrapper(args...)

    return Distributions._logpdf(unkeyed_dist, x)
end

_maybe_parent(x) = x
_maybe_parent(x::AbstractArray) = parent(x)

# Also need to overload `rand` methods to return KeyedArrays:

function Base.rand(rng::AbstractRNG, d::KeyedDistOrSampleable)
    sample = rand(rng, distribution(d))
    ndims(sample) == 0 && return sample  # univariate returns a Number
    return KeyedArray(sample, axiskeys(d))
end

function Base.rand(rng::AbstractRNG, d::KeyedDistOrSampleable, n::Int)
    samples = rand(rng, distribution(d), n)
    ndims(samples) == 1 && return KeyedArray(samples, Base.OneTo(n))  # univariate
    return KeyedArray(samples, (first(axiskeys(d)), Base.OneTo(n)))
end

function Base.rand(rng::AbstractRNG, d::KeyedDistribution{<:Matrixvariate}, n::Int)
    # Distributions.rand returns a vector of matrices
    samples = [KeyedArray(x, axiskeys(d)) for x in rand(rng, distribution(d), n)]
    return KeyedArray(samples, Base.OneTo(n))
end

# Statistics functions for Distribution

Distributions.mean(d::KeyedDistribution) = KeyedArray(mean(distribution(d)), axiskeys(d))

Distributions.var(d::KeyedDistribution) = KeyedArray(var(distribution(d)), axiskeys(d))

function Distributions.cov(d::KeyedDistribution)
    keys = vcat(axiskeys(d)...)
    return KeyedArray(cov(distribution(d)), (keys, keys))
end

Distributions.entropy(d::KeyedDistribution) = entropy(distribution(d))
Distributions.entropy(d::KeyedDistribution, b::Real) = entropy(distribution(d), b)

# Univariate Distributions only

for f in (:logpdf, :quantile, :mgf, :cf)
    @eval Distributions.$f(d::KeyedDistribution{<:Univariate}, x) = $f(distribution(d), x)
end

for f in (:minimum, :maximum, :modes, :mode, :skewness, :kurtosis)
    @eval Distributions.$f(d::KeyedDistribution{<:Univariate}) = $f(distribution(d))
end

# Needed to avoid method ambiguity
Distributions.cdf(d::KeyedDistribution{<:Univariate}, x::Real) = cdf(distribution(d), x)

function Distributions.insupport(d::KeyedDistribution{<:Univariate}, x::Real)
    return insupport(distribution(d), x)
end

end

max_stars_repo_name: invenia/KeyedDistributions

max_issues_repo_name: invenia/KeyedDistributions

max_forks_repo_name: invenia/KeyedDistributions


text:
using DataFrames
using DatasToolbox

const NROWS = 5*10^6

df = DataFrame(A=rand(NROWS), B=rand(1:100, NROWS))

const M = 0.5

# testfunc(a, b) = a < M && b % 3 == 0
# testfunc(a) = a < M

# @time ocdf = DatasToolbox.constrain_OLD(df, [:A, :B], testfunc)
# @time cdf = constrain(df, [:A, :B], testfunc)

dict = Dict()
expr = :((:A .< M) && (:B % 3 .== 0))
DatasToolbox._checkConstraintExpr!(expr, dict)

mac = macroexpand(:(@constrain(df, (:A .< M) && (:B % 3 .== 0))))

@time cdf = @constrain(df, (:A .< M) && (:B % 3 .== 0))

@time ccdf = constrain(df, A=(a -> a < M), B=(b -> b % 3 == 0))

max_stars_repo_name: JuliaTagBot/DatasToolbox

max_issues_repo_name: JuliaTagBot/DatasToolbox

max_forks_repo_name: JuliaTagBot/DatasToolbox


text:
"""
    FieldBoundaryConditions

An alias for `NamedTuple{(:x, :y, :z)}` that represents a set of three `CoordinateBoundaryCondition`s
applied to a field along x, y, and z.
"""
const FieldBoundaryConditions = NamedTuple{(:x, :y, :z)}

"""
    FieldBoundaryConditions(x, y, z)

Construct a `FieldBoundaryConditions` using a `CoordinateBoundaryCondition` for each of the
`x`, `y`, and `z` coordinates.
"""
FieldBoundaryConditions(x, y, z) = FieldBoundaryConditions((x, y, z))

"""
    HorizontallyPeriodicBCs(;   top = BoundaryCondition(Flux, nothing),
                             bottom = BoundaryCondition(Flux, nothing))

Construct `FieldBoundaryConditions` with `Periodic` boundary conditions in the x and y
directions and specified `top` (+z) and `bottom` (-z) boundary conditions for u, v,
and tracer fields.

`HorizontallyPeriodicBCs` cannot be applied to the the vertical velocity w.
"""
function HorizontallyPeriodicBCs(;    top = BoundaryCondition(Flux, nothing),
                                   bottom = BoundaryCondition(Flux, nothing))

    x = PeriodicBCs()
    y = PeriodicBCs()
    z = CoordinateBoundaryConditions(bottom, top)

    return FieldBoundaryConditions(x, y, z)
end

"""
    ChannelBCs(; north = BoundaryCondition(Flux, nothing),
                 south = BoundaryCondition(Flux, nothing),
                   top = BoundaryCondition(Flux, nothing),
                bottom = BoundaryCondition(Flux, nothing))

Construct `FieldBoundaryConditions` with `Periodic` boundary conditions in the x
direction and specified `north` (+y), `south` (-y), `top` (+z) and `bottom` (-z)
boundary conditions for u, v, and tracer fields.

`ChannelBCs` cannot be applied to the the vertical velocity w.
"""
function ChannelBCs(;  north = BoundaryCondition(Flux, nothing),
                       south = BoundaryCondition(Flux, nothing),
                         top = BoundaryCondition(Flux, nothing),
                      bottom = BoundaryCondition(Flux, nothing)
                    )

    x = PeriodicBCs()
    y = CoordinateBoundaryConditions(south, north)
    z = CoordinateBoundaryConditions(bottom, top)

    return FieldBoundaryConditions(x, y, z)
end


text:
module LogDensity

export density,
    logdensity

"""
    logdensity(μ, x)

Calculate the density of the measure-like `μ` at `x` in log-space.

`logdensity` should be well-defined over the extension of `μ`'s support to the whole range of the
type of compatible `x` values, returning an equivalent of `-Inf` outside the support (for meaningful
values, ignoring \"undefined\" cases like `NaN`)".
"""
function logdensity end

"""
    density(μ, x)

Calculate the density of the measure-like `μ` at `x`, in normal (i.e., non-log) space.  Should
fulfill the property

```julia
density(μ, x) = exp(logdensity(μ, x))
```

"""
density(args..) = exp(logdensity(args...))

end # module


text:
using MIToS.MSA

read("../data/PF00089.stockholm.gz", Stockholm, MultipleSequenceAlignment, deletefullgaps=false)
read("../data/PF00089.sth", Stockholm, MultipleSequenceAlignment, deletefullgaps=false)
read("../data/PF00089.fasta.gz", FASTA, MultipleSequenceAlignment, deletefullgaps=false)
read("../data/PF00089.fasta", FASTA, MultipleSequenceAlignment, deletefullgaps=false)

println("[BENCH] Read compressed Stockholm PF00089: ", @elapsed read("../data/PF00089.stockholm.gz", Stockholm, MultipleSequenceAlignment, deletefullgaps=false))
println("[BENCH] Read uncompressed Stockholm PF00089: ", @elapsed read("../data/PF00089.sth", Stockholm, MultipleSequenceAlignment, deletefullgaps=false))
println("[BENCH] Read compressed FASTA PF00089: ", @elapsed read("../data/PF00089.fasta.gz", FASTA, MultipleSequenceAlignment, deletefullgaps=false))
println("[BENCH] Read uncompressed FASTA PF00089: ", @elapsed read("../data/PF00089.fasta", FASTA, MultipleSequenceAlignment, deletefullgaps=false))

read("../data/PF16957.stockholm.gz", Stockholm, MultipleSequenceAlignment, deletefullgaps=false)
read("../data/PF16957.sth", Stockholm, MultipleSequenceAlignment, deletefullgaps=false)
read("../data/PF16957.fasta.gz", FASTA, MultipleSequenceAlignment, deletefullgaps=false)
read("../data/PF16957.fasta", FASTA, MultipleSequenceAlignment, deletefullgaps=false)

println("[BENCH] Read compressed Stockholm PF16957: ", @elapsed read("../data/PF16957.stockholm.gz", Stockholm, MultipleSequenceAlignment, deletefullgaps=false))
println("[BENCH] Read uncompressed Stockholm PF16957: ", @elapsed read("../data/PF16957.sth", Stockholm, MultipleSequenceAlignment, deletefullgaps=false))
println("[BENCH] Read compressed FASTA PF16957: ", @elapsed read("../data/PF16957.fasta.gz", FASTA, MultipleSequenceAlignment, deletefullgaps=false))
println("[BENCH] Read uncompressed FASTA PF16957: ", @elapsed read("../data/PF16957.fasta", FASTA, MultipleSequenceAlignment, deletefullgaps=false))

max_stars_repo_name: diegozea/mitos-benchmarks

max_issues_repo_name: diegozea/mitos-benchmarks

max_forks_repo_name: diegozea/mitos-benchmarks


text:
const REGEX = r"([+-]?\d*)\.?(\d*)[Ee]?([+-]?\d+)?"

bigint(s::AbstractString) = isempty(s) ? big"0" : parse(BigInt, s)
bigint(::Nothing) = 0

function parse(::Type{Decimal}, s::AbstractString)
    int, digit, exp = match(REGEX, s).captures
    ndigits = length(digit)
    exp = bigint(exp) - ndigits
    c = bigint(int * digit)
    Decimal(c, exp)
end

Decimal(s::AbstractString) = parse(Decimal, s)
Decimal(x::Real) = Decimal(string(x))

macro d_str(s::String)
    parse(Decimal, s)
end

function print(io::IO, x::Decimal)
    x = strip_trailing_zeros(x)
    isneg = sign(x) < 0
    s = string(abs(x.c))
    len = length(s)
    exp = x.q + len - 1

    if abs(exp) >= PRECISION
        s = len == 1 ? s : s[1] * '.' * s[2:end]
        s *=  'E' * string(exp)
    elseif exp >= 0
        s = x.q >= 0 ? s * '0'^x.q : s[1:exp+1] * '.' * s[exp+2:end]
    else
        n = abs(exp)
        s = '0'^n * s
        s = s[1] * '.' * s[2:end]
    end
    print(io, sign(x) < 0 ? '-' * s : s)
end


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

"""
    groupby(sdata, var)

Partition spatial data `sdata` into groups of constant value
for spatial variable `var`.

### Notes

Missing values are grouped into a separate group.
"""
groupby(sdata::AbstractData, var::Symbol) =
  partition(sdata, VariablePartitioner(var))


text:
import Base.DFT: plan_fft!, plan_ifft!, to1
import Base.DFT.FFTW.cFFTWPlan
import ..Log

typealias GMatrix Matrix{SensitiveFloat{Float64}}
typealias fs0mMatrix Matrix{SensitiveFloat{Float64}}
typealias fs1mMatrix Matrix{SensitiveFloat{Float64}}

const plan_fft_lock = Base.Threads.SpinLock()

function __init__()
    global fft_plans, ifft_plans
    fft_plans = Dict{Tuple{Int64, Int64}, cFFTWPlan{Complex{Float64},-1,true,2}}()
    ifft_type = Base.DFT.ScaledPlan{Complex{Float64},Base.DFT.FFTW.cFFTWPlan{Complex{Float64},1,true,2},Float64}
    ifft_plans = Dict{Tuple{Int64, Int64}, ifft_type}()
end

function safe_fft!(A)
    A1 = to1(A)
    global fft_plans

    if !haskey(fft_plans, size(A))
        lock(plan_fft_lock)
        fft_plans[size(A)] = Base.DFT.plan_fft!(A1)
        unlock(plan_fft_lock)
    end

    fft_plans[size(A)] * A1  # mutates A1
    A1
end

function safe_ifft!(A)
    A1 = to1(A)
    global ifft_plans

    if !haskey(ifft_plans, size(A))
        lock(plan_fft_lock)
        ifft_plans[size(A)] = Base.DFT.plan_ifft!(A1)
        unlock(plan_fft_lock)
    end

    ifft_plans[size(A)] * A1  # mutates A1
    A1
end

type FSMSensitiveFloatMatrices
    # The lower corner of the image (in terms of index values)
    h_lower::Int
    w_lower::Int

    fs1m_image::fs1mMatrix;
    fs1m_image_padded::fs1mMatrix;
    fs1m_conv::fs1mMatrix;
    fs1m_conv_padded::fs1mMatrix;

    # We convolve the star fs0m directly using Lanczos interpolation.
    fs0m_conv::fs0mMatrix;

    E_G::GMatrix
    var_G::GMatrix

    # The PSF and its FFT.
    psf_fft::Matrix{Complex{Float64}}
    psf::Matrix{Float64}

    # The amount of padding introduced by the convolution on one side of the
    # image (the total pixels added in each dimension are twice this)
    pad_pix_h::Int
    pad_pix_w::Int

    # Functions for star convolution.
    # kernel_fun is a kernel function that returns the value, derivative, and second
    # derivative of a univariate kernel, e.g. bspline_kernel_with_derivatives()
    kernel_width::Int64
    kernel_fun

    FSMSensitiveFloatMatrices() = begin
        new(1, 1,
            fs1mMatrix(), fs1mMatrix(), fs1mMatrix(), fs1mMatrix(), fs0mMatrix(),
            GMatrix(), GMatrix(),
            Matrix{Complex{Float64}}(),
            Matrix{Float64}(),
            0, 0,
            2, x -> cubic_kernel_with_derivatives(x, 0.0))
    end
end

function initialize_fsm_sf_matrices_band!(
    fsms::FSMSensitiveFloatMatrices,
    s::Int, b::Int, num_active_sources::Int,
    h_lower::Int, w_lower::Int,
    h_upper::Int, w_upper::Int,
    psf_image::Matrix{Float64})

    psf_size = size(psf_image)

    fsms.h_lower = h_lower
    fsms.w_lower = w_lower

    h_width = h_upper - h_lower + 1
    w_width = w_upper - w_lower + 1

    # An fsm value is only sensitive to one source's parameters.
    fsms.fs1m_image = zero_sensitive_float_array(
        Float64, length(GalaxyPosParams), 1, h_width, w_width);
    fsms.fs1m_conv = zero_sensitive_float_array(
        Float64, length(GalaxyPosParams), 1, h_width, w_width);

    fsms.fs0m_conv = zero_sensitive_float_array(
        Float64, length(StarPosParams), 1, h_width, w_width);

    # The amount of padding introduced by the convolution
    (fft_size1, fft_size2) =
        (h_width + psf_size[1] - 1, w_width + psf_size[2] - 1)
    # Make sure that the PSF has an odd dimension.
    @assert psf_size[1] % 2 == 1
    @assert psf_size[2] % 2 == 1
    fsms.pad_pix_h = Integer((psf_size[1] - 1) / 2)
    fsms.pad_pix_w = Integer((psf_size[2] - 1) / 2)

    fsms.fs1m_image_padded = zero_sensitive_float_array(
        Float64, length(GalaxyPosParams), 1, fft_size1, fft_size2);
    fsms.fs1m_conv_padded = zero_sensitive_float_array(
        Float64, length(GalaxyPosParams), 1, fft_size1, fft_size2);

    # Brightness images
    fsms.E_G = zero_sensitive_float_array(
        Float64, length(CanonicalParams), num_active_sources, h_width, w_width);
    fsms.var_G = zero_sensitive_float_array(
        Float64, length(CanonicalParams), num_active_sources, h_width, w_width);

    # Store the psf image and its FFT.
    fsms.psf = deepcopy(psf_image)
    fsms.psf_fft = zeros(Complex{Float64}, fft_size1, fft_size2)
    fsms.psf_fft[1:psf_size[1], 1:psf_size[2]] = fsms.psf

    safe_fft!(fsms.psf_fft)
end

function initialize_fsm_sf_matrices!(
    fsm_mat::Matrix{FSMSensitiveFloatMatrices},
    ea::ElboArgs{Float64},
    psf_image_mat::Matrix{Matrix{Float64}})

    num_active_sources = length(ea.active_sources)

    for n in 1:ea.N, s in 1:ea.S
        p = ea.patches[s, n]
        apb = p.active_pixel_bitmap

        # some sources don't appear in some images
        sum(apb) > 0 || continue

        active_cols = find([ any(apb[:, col]) for col in 1:size(apb, 2) ])
        active_rows = find([ any(apb[row, :]) for row in 1:size(apb, 1) ])
        h1 = p.bitmap_offset[1] + minimum(active_rows)
        w1 = p.bitmap_offset[2] + minimum(active_cols)
        h2 = h1 + maximum(active_rows) - 1
        w2 = w1 + maximum(active_cols) - 1

        initialize_fsm_sf_matrices_band!(
            fsm_mat[s, n], s, n, num_active_sources,
            h1, w1, h2, w2, psf_image_mat[s, n])
    end
end

# This is useful for debugging and exploration.
# Do not delete until May 2017.
function debug_populate_fsm_mat!(
    ea::ElboArgs,
    fsm_mat::Matrix{FSMSensitiveFloatMatrices})

    sbs = load_source_brightnesses(ea,
        calculate_gradient=ea.elbo_vars.elbo.has_gradient,
        calculate_hessian=ea.elbo_vars.elbo.has_hessian);

    gal_mcs_vec = Array(Array{GalaxyCacheComponent{Float64}, 4}, ea.N);
    for b=1:ea.N
        gal_mcs_vec[b] = load_gal_bvn_mixtures(
                ea.S, ea.patches, ea.vp, ea.active_sources, b,
                calculate_gradient=ea.elbo_vars.elbo.has_gradient,
                calculate_hessian=ea.elbo_vars.elbo.has_hessian);
    end

    for b=1:ea.N
        for s in 1:ea.S
            fsms = fsm_mat[s, b]
            clear_brightness!(fsms)
            populate_star_fsm_image!(
                ea, s, b, fsms.psf, fsms.fs0m_conv,
                fsms.h_lower, fsms.w_lower,
                fsms.kernel_fun, fsms.kernel_width)
            populate_gal_fsm_image!(
                ea, s, b, gal_mcs_vec[b], fsm_mat[s, n])
            accumulate_source_image_brightness!(
                ea, s, b, fsms, sbs[s])
        end
    end
end


text:
module Outputs

using REPL.Terminals: TTYTerminal
using REPL.TerminalMenus: RadioMenu, request

using Crayons: Crayon
using Crayons.Box: GREEN_FG, BLUE_FG, RED_FG
using Parameters: type2dict
using PrettyTables: Highlighter, pretty_table, ft_printf
using QuantumESPRESSOBase.Inputs.PWscf: AtomicPositionsCard, CellParametersCard
using QuantumESPRESSOParsers.Outputs.PWscf

export PWOutput, output_parser

struct PWOutput end

function output_parser(terminal::TTYTerminal, ::Type{T}) where {T<:PWOutput}
    calculations = pairs(("scf", "nscf", "bands", "relax", "md", "vc-relax", "vc-md"))
    print(terminal, GREEN_FG("Please give the absolute path to your output file: "))
    path = abspath(strip(readline(terminal)))
    str = try
        open(path, "r") do io
            read(io, String)
        end
    catch
        @warn("File '$path' not found!")
    end
    if isjobdone(str)
        println(terminal, GREEN_FG("The job is done! Ready to parse!"))
    else
        println(terminal, RED_FG("The job is not finished, be careful!"))
    end
    calculation = calculations[request(
        terminal,
        GREEN_FG("What exact calculation is this output?") |> string,
        RadioMenu(collect(values(calculations))),
    )]
    if calculation ∈ ("relax", "vc-relax")
        choice = request(
            terminal,
            GREEN_FG("Do you want to parse the final or all atomic positions?") |> string,
            RadioMenu(["final", "all"]),
        )
        if choice == 1
            ap = tryparsefinal(AtomicPositionsCard, str)
            println(terminal, BLUE_FG("Print the final atomic positions:"))
            display(ap.data)
        else
            ap = tryparseall(AtomicPositionsCard, str)
            println(terminal, BLUE_FG("Print all atomic positions:"))
            foreach(x -> display(x.data), ap)
        end
        if calculation == "vc-relax"
            choice = request(
                terminal,
                GREEN_FG("Do you want to parse the final or all cell parameters?") |>
                string,
                RadioMenu(["final", "all"]),
            )
            if choice == 1
                cp = tryparsefinal(CellParametersCard, str)
                println(terminal, BLUE_FG("Print the final cell parameters:"))
                println(terminal, cp.data)
            else
                cp = tryparseall(CellParametersCard, str)
                println(terminal, BLUE_FG("Print all cell parameters:"))
                foreach(x -> display(x.data), cp)
            end
            if isrelaxed(str)
                println(terminal, GREEN_FG("The structure is relaxed!"))
            else
                println(terminal, RED_FG("The structure is not well-relaxed!"))
            end
        end
    end
    hl_odd =
        Highlighter(f = (data, i, j) -> (i % 2) == 0, crayon = Crayon(background = :blue))
    choice = request(
        terminal,
        GREEN_FG("Do you want to parse its summary?") |> string,
        RadioMenu(["yes", "no"]),
    )
    if choice == 1
        preamble = parse(Preamble, str)
        pretty_table(preamble |> type2dict; highlighters = hl_odd)
    end
    choice = request(
        terminal,
        GREEN_FG("Do you want to parse the energies?") |> string,
        RadioMenu(["yes", "no"]),
    )
    if choice == 1
        df = parse_electrons_energies(str, :combined)
        pretty_table(df; highlighters = hl_odd, formatter = ft_printf("%10.5f"))
    end
    choice = request(
        terminal,
        GREEN_FG("Do you want to parse the time used?") |> string,
        RadioMenu(["yes", "no"]),
    )
    if choice == 1
        df = parse_clock(str)
        pretty_table(df; highlighters = hl_odd, formatter = ft_printf("%10.5f"))
    end
    choice = request(
        terminal,
        GREEN_FG("Do you want to parse the diagonalization info?") |> string,
        RadioMenu(["yes", "no"]),
    )
    if choice == 1
        df = Outputs.PWscf.parse_diagonalization(str)
        pretty_table(df; highlighters = hl_odd, formatter = ft_printf("%10.5f"))
    end
end # function output_parser

end


text:
"""
    read_nnet(fname::String)

Read in neural net from file and return Network struct
"""
function read_nnet(fname::String)
    f = open(fname)
    line = readline(f)
    while occursin("//", line) #skip comments
        line = readline(f)
    end
    # read in layer sizes and populate array
    record = split(line, ",")
    nLayers = parse(Int64, record[1])
    record = split(readline(f), ",")
    layerSizes = Vector{Int64}(undef, nLayers + 1)
    for i = 1: nLayers + 1
        layerSizes[i] = parse(Int64, record[i])
    end
    # read past additonal information
    for i = 1: 5
        line = readline(f)
    end
    # initialize layers
    layers = Vector{Layer}(undef, nLayers)
    for i = 1:nLayers
        curr_layer = init_layer(i, layerSizes, f)
        layers[i] = curr_layer
    end

    return Network(layers)
end

"""
    init_layer(i::Int64, layerSizes::Array{Int64}, f::IOStream)

Read in layer from nnet file and return a Layer struct containing its weights/biases
"""
function init_layer(i::Int64, layerSizes::Array{Int64}, f::IOStream)
     bias = Vector{Float64}(undef, layerSizes[i+1])
     weights = Matrix{Float64}(undef, layerSizes[i+1], layerSizes[i])
     # first read in weights
     for r = 1: layerSizes[i+1]
        line = readline(f)
        record = split(line, ",")
        token = record[1]
        c = 1
        for c = 1: layerSizes[i]
            weights[r, c] = parse(Float64, token)
            token = record[c]
        end
     end
     # now read in bias
     for r = 1: layerSizes[i+1]
        line = readline(f)
        record = split(line, ",")
        bias[r] = parse(Float64, record[1])
     end
     # activation function is set to ReLU as default
     return Layer(weights, bias, ReLU())
end

"""
    compute_output(nnet::Network, input::Vector{Float64})

Propagate a given vector through a nnet and compute the output.
"""
function compute_output(nnet::Network, input::Vector{Float64})
    curr_value = input
    layers = nnet.layers
    for i = 1:length(layers) # layers does not include input layer (which has no weights/biases)
        curr_value = (layers[i].weights * curr_value) + layers[i].bias
        curr_value = layers[i].activation(curr_value)
    end
    return curr_value # would another name be better?
end

# """
#     compute_output(nnet::Network, input::Vector{Float64}, i, j)

# Returns ouput of neuron j in layer i for a given input.
# NOTE: The was necessary for the sampling methods, but now might not be.
# """
# function compute_output(nnet::Network, input::Vector{Float64}, i, j)
#     layers = nnet.layers
#     @assert 0 <= i <= length(layers)         "number of layers in nnet is $(length(layers)). Got $i for number of layers to compute"
#     @assert 0 <= j <= length(layers[i].bias) "number of neurons in layer is $(length(layers[i].bias)). Got $j for neuron index"
#     curr_value = input
#     for m = 1:i
#         curr_value = (layers[m].weights * curr_value) + layers[m].bias
#         curr_value = layers[m].activation(curr_value)
#     end
#     return curr_value[j]
# end

"""
    get_activation(nnet::Network, x::Vector{Float64})

Given a network, find the activation pattern of all neurons at a given point x.
Assume ReLU.
return Vector{Vector{Bool}}. Each Vector{Bool} refers to the activation pattern of a particular layer.
"""
function get_activation(nnet::Network, x::Vector{Float64})
    act_pattern = Vector{Vector{Bool}}(undef, length(nnet.layers))
    curr_value = x
    for (i, layer) in enumerate(nnet.layers)
        curr_value = layer.weights * curr_value + layer.bias
        act_pattern[i] = curr_value .>= 0.0
        curr_value = layer.activation(curr_value)
    end
    return act_pattern
end

"""
    get_activation(nnet::Network, input::Hyperrectangle)

Given a network, find the activation pattern of all neurons for a given input set.
Assume ReLU.
return Vector{Vector{Int64}}.
- 1: activated
- 0: undetermined
- -1: not activated
"""
function get_activation(nnet::Network, input::Hyperrectangle)
    bounds = get_bounds(nnet, input)
    return get_activation(nnet, bounds)
end

"""
    get_activation(nnet::Network, bounds::Vector{Hyperrectangle})

Given a network, find the activation pattern of all neurons given the node-wise bounds.
Assume ReLU.
return Vector{Vector{Int64}}.
- 1: activated
- 0: undetermined
- -1: not activated
"""
function get_activation(nnet::Network, bounds::Vector{Hyperrectangle})
    act_pattern = Vector{Vector{Int}}(undef, length(nnet.layers))
    for (i, layer) in enumerate(nnet.layers)
        before_act_bound = linear_transformation(layer, bounds[i])
        lower = low(before_act_bound)
        upper = high(before_act_bound)
        act_pattern[i] = fill(0, length(layer.bias))
        for j in 1:length(layer.bias) # For evey node
            if lower[j] > 0.0
                act_pattern[i][j] = 1
            elseif upper[j] < 0.0
                act_pattern[i][j] = -1
            end
        end
    end
    return act_pattern
end

"""
    get_gradient(nnet::Network, x::Vector)

Given a network, find the gradient at the input x
"""
function get_gradient(nnet::Network, x::Vector)
    z = x
    gradient = Matrix(1.0I, length(x), length(x))
    for (i, layer) in enumerate(nnet.layers)
        z_hat = layer.weights * z + layer.bias
        σ_gradient = act_gradient(layer.activation, z_hat)
        gradient = Diagonal(σ_gradient) * layer.weights * gradient
        z = layer.activation(z_hat)
    end
    return gradient
end

"""
    act_gradient(act::ReLU, z_hat::Vector{N}) where N

Computing the gradient of an activation function at point z_hat.
Currently only support ReLU.
"""
function act_gradient(act::ReLU, z_hat::Vector)
    return z_hat .>= 0.0
end

"""
    get_gradient(nnet::Network, input::AbstractPolytope)

Get lower and upper bounds on network gradient for a given input set.
Return:
- `LG::Vector{Matrix}`: lower bounds
- `UG::Vector{Matrix}`: upper bounds
"""
function get_gradient(nnet::Network, input::AbstractPolytope)
    LΛ, UΛ = act_gradient_bounds(nnet, input)
    return get_gradient(nnet, LΛ, UΛ)
end

"""
    act_gradient_bounds(nnet::Network, input::AbstractPolytope)

Computing the bounds on the gradient of all activation functions given an input set.
Currently only support ReLU.
Return:
- `LΛ::Vector{Matrix}`: lower bounds
- `UΛ::Vector{Matrix}`: upper bounds
"""
function act_gradient_bounds(nnet::Network, input::AbstractPolytope)
    bounds = get_bounds(nnet, input)
    LΛ = Vector{Matrix}(undef, 0)
    UΛ = Vector{Matrix}(undef, 0)
    for (i, layer) in enumerate(nnet.layers)
        before_act_bound = linear_transformation(layer, bounds[i])
        lower = low(before_act_bound)
        upper = high(before_act_bound)
        l = act_gradient(layer.activation, lower)
        u = act_gradient(layer.activation, upper)
        push!(LΛ, Diagonal(l))
        push!(UΛ, Diagonal(u))
    end
    return (LΛ, UΛ)
end

"""
    get_gradient(nnet::Network, LΛ::Vector{Matrix}, UΛ::Vector{Matrix})

Get lower and upper bounds on network gradient for given gradient bounds on activations
Inputs:
- `LΛ::Vector{Matrix}`: lower bounds on activation gradients
- `UΛ::Vector{Matrix}`: upper bounds on activation gradients
Return:
- `LG::Vector{Matrix}`: lower bounds
- `UG::Vector{Matrix}`: upper bounds
"""
function get_gradient(nnet::Network, LΛ::Vector{Matrix}, UΛ::Vector{Matrix})
    n_input = size(nnet.layers[1].weights, 2)
    LG = Matrix(1.0I, n_input, n_input)
    UG = Matrix(1.0I, n_input, n_input)
    for (i, layer) in enumerate(nnet.layers)
        LG_hat, UG_hat = interval_map(layer.weights, LG, UG)
        LG = LΛ[i] * max.(LG_hat, 0) + UΛ[i] * min.(LG_hat, 0)
        UG = LΛ[i] * min.(UG_hat, 0) + UΛ[i] * max.(UG_hat, 0)
    end
    return (LG, UG)
end

"""
    get_gradient(nnet::Network, LΛ::Vector{Vector{N}}, UΛ::Vector{Vector{N}}) where N

Get lower and upper bounds on network gradient for given gradient bounds on activations
Inputs:
- `LΛ::Vector{Vector{N}}`: lower bounds on activation gradients
- `UΛ::Vector{Vector{N}}`: upper bounds on activation gradients
Return:
- `LG::Vector{Matrix}`: lower bounds
- `UG::Vector{Matrix}`: upper bounds
"""
function get_gradient(nnet::Network, LΛ::Vector{Vector{N}}, UΛ::Vector{Vector{N}}) where N
    n_input = size(nnet.layers[1].weights, 2)
    LG = Matrix(1.0I, n_input, n_input)
    UG = Matrix(1.0I, n_input, n_input)
    for (i, layer) in enumerate(nnet.layers)
        LG_hat, UG_hat = interval_map(layer.weights, LG, UG)
        LG = Diagonal(LΛ[i]) * max.(LG_hat, 0) + Diagonal(UΛ[i]) * min.(LG_hat, 0)
        UG = Diagonal(LΛ[i]) * min.(UG_hat, 0) + Diagonal(UΛ[i]) * max.(UG_hat, 0)
    end
    return (LG, UG)
end

"""
    interval_map(W::Matrix{N}, l::Vector{N}, u::Vector{N}) where N

Simple linear mapping on intervals
Inputs:
- `W::Matrix{N}`: the linear mapping
- `l::Vector{N}`: the lower bound
- `u::Vector{N}`: the upper bound
Outputs:
- `l_new::Vector{N}`: the lower bound after mapping
- `u_new::Vector{N}`: the upper bound after mapping
"""
function interval_map(W::Matrix{N}, l::Vector{N}, u::Vector{N}) where N
    l_new = max.(W, 0) * l + min.(W, 0) * u
    u_new = max.(W, 0) * u + min.(W, 0) * l
    return (l_new, u_new)
end

"""
    interval_map(W::Matrix{N}, l::Vector{N}, u::Vector{N}) where N

Simple linear mapping on intervals
Inputs:
- `W::Matrix{N}`: the linear mapping
- `l::Matrix{N}`: the lower bound
- `u::Matrix{N}`: the upper bound
Outputs:
- `l_new::Matrix{N}`: the lower bound after mapping
- `u_new::Matrix{N}`: the upper bound after mapping
"""
function interval_map(W::Matrix{N}, l::Matrix{N}, u::Matrix{N}) where N
    l_new = max.(W, 0) * l + min.(W, 0) * u
    u_new = max.(W, 0) * u + min.(W, 0) * l
    return (l_new, u_new)
end

"""
    get_bounds(nnet::Network, input::Hyperrectangle)

This function calls maxSens to compute node-wise bounds given a input set.

Return:
- `bounds::Vector{Hyperrectangle}`: bounds for all nodes AFTER activation. `bounds[1]` is the input set.
"""
function get_bounds(nnet::Network, input::Hyperrectangle) # NOTE there is another function by the same name in convDual. Should reconsider dispatch
    solver = MaxSens(0.0, true)
    bounds = Vector{Hyperrectangle}(undef, length(nnet.layers) + 1)
    bounds[1] = input
    for (i, layer) in enumerate(nnet.layers)
        bounds[i+1] = forward_layer(solver, layer, bounds[i])
    end
    return bounds
end

"""
    get_bounds(problem::Problem)

This function calls maxSens to compute node-wise bounds given a problem.

Return:
- `bounds::Vector{Hyperrectangle}`: bounds for all nodes AFTER activation. `bounds[1]` is the input set.
"""
get_bounds(problem::Problem) = get_bounds(problem.network, problem.input)

"""
    get_bounds(nnet::Network, input::Hyperrectangle, act::Bool)

Compute bounds before or after activation by interval arithmetic. To be implemented.

Inputs:
- `nnet::Network`: network
- `input::Hyperrectangle`: input set
- `act::Bool`: `true` for after activation bound; `false` for before activation bound
Return:
- `bounds::Vector{Hyperrectangle}`: bounds for all nodes AFTER activation. `bounds[1]` is the input set.
"""
function get_bounds(nnet::Network, input::Hyperrectangle, act::Bool)
    if act
        return get_bounds(nnet, input)
    else
        error("before activation bounds not supported yet.")
    end
end

"""
    linear_transformation(layer::Layer, input::Hyperrectangle)

Transformation of a set considering linear mappings in a layer.

Inputs:
- `layer::Layer`: a layer in a network
- `input::Hyperrectangle`: input set
Return:
- `output::Hyperrectangle`: set after transformation.
"""
function linear_transformation(layer::Layer, input::Hyperrectangle)
    (W, b, act) = (layer.weights, layer.bias, layer.activation)
    before_act_center = W * input.center + b
    before_act_radius = abs.(W) * input.radius
    return Hyperrectangle(before_act_center, before_act_radius)
end

"""
    linear_transformation(layer::Layer, input::HPolytope)

Transformation of a set considering linear mappings in a layer.

Inputs:
- `layer::Layer`: a layer in a network
- `input::HPolytope`: input set
Return:
- `output::HPolytope`: set after transformation.
"""
function linear_transformation(layer::Layer, input::HPolytope)
    (W, b) = (layer.weights, layer.bias)
    input_v = tovrep(input)
    output_v = [W * v + b for v in vertices_list(input_v)]
    output = tohrep(VPolytope(output_v))
    return output
end

"""
    linear_transformation(W::Matrix, input::HPolytope)

Transformation of a set considering a linear mapping.

Inputs:
- `W::Matrix`: a linear mapping
- `input::HPolytope`: input set
Return:
- `output::HPolytope`: set after transformation.
"""
function linear_transformation(W::Matrix, input::HPolytope)
    input_v = tovrep(input)
    output_v = [W * v for v in vertices_list(input_v)]
    output = tohrep(VPolytope(output_v))
    return output
end

"""
    split_interval(dom::Hyperrectangle, index::Int64)

Split a set into two at the given index.

Inputs:
- `dom::Hyperrectangle`: the set to be split
- `index`: the index to split at
Return:
- `(left, right)::Tuple{Hyperrectangle, Hyperrectangle}`: two sets after split
"""
function split_interval(dom::Hyperrectangle, index::Int64)
    input_lower, input_upper = low(dom), high(dom)

    input_upper[index] = dom.center[index]
    input_split_left = Hyperrectangle(low = input_lower, high = input_upper)

    input_lower[index] = dom.center[index]
    input_upper[index] = dom.center[index] + dom.radius[index]
    input_split_right = Hyperrectangle(low = input_lower, high = input_upper)
    return (input_split_left, input_split_right)
end


text:
#
# Hooks into AstrodynamicalModels.jl
#

"""
$(SIGNATURES)

Returns the `ModelingToolkit.ODESystem` associated
with `R2BPParameters`, provided by `AstrodynamicalModels`.
"""
model(::R2BPParameters) = AstrodynamicalModels.R2BP

"""
$(SIGNATURES)

Returns the `ModelingToolkit.ODESystem` associated
with `CR3BPParameters`, provided by `AstrodynamicalModels`.
"""
model(::CR3BPParameters) = AstrodynamicalModels.CR3BP

"""
$(SIGNATURES)

Returns the `DifferentialEquations.ODEFunction` associated
with `R2BPParameters`, provided by `AstrodynamicalModels`.
"""
vectorfield(::R2BPParameters) = AstrodynamicalModels.R2BPVectorField

"""
$(SIGNATURES)

Returns the `DifferentialEquations.ODEFunction` associated
with `CR3BPParameters`, provided by `AstrodynamicalModels`.
"""
vectorfield(::CR3BPParameters) = AstrodynamicalModels.CR3BPVectorField


text:
using Test: @test_nowarn, @test
using SHA
using Trixi
using Trixi2Img

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples")

function run_trixi(elixir; parameters...)
  @test_nowarn trixi_include(joinpath(EXAMPLES_DIR, elixir); parameters...)
end

function sha1file(filename)
  open(filename) do f
    bytes2hex(sha1(f))
  end
end

function set_pdf_creation_date(filename, date="20200101080000")
  # Implementation based on https://discourse.julialang.org/t/removing-the-first-line-in-a-text-file/35521/12
  write(filename,
        open(filename) do f
          content = read(f, String)
          replace(content, r"/CreationDate \(D:[0-9]*\)"=>"/CreationDate (D:$date)");
        end)
end

function test_trixi2img(filenames, outdir; hashes=nothing, kwargs...)
  @test_nowarn trixi2img(joinpath(outdir, filenames); output_directory=outdir, kwargs...)

  if !isnothing(hashes)
    for (filename, hash_expected) in hashes
      # If output format is PDF, we need to reset the CreationDate to be able to compare hashes
      if endswith(filename, ".pdf")
        set_pdf_creation_date(joinpath(outdir, filename))
      end

      hash_measured = sha1file(joinpath(outdir, filename))
      @test hash_expected == hash_measured
    end
  end
end


text:
function estimate_value_action(estimator::Union{BasicPOMCP.SolvedPORollout,BasicPOMCP.SolvedFORollout}, pomdp::POMDPs.POMDP, start_state, h::BeliefNode, steps::Int)
    rollout_action(estimator, pomdp, start_state, h, steps)
end

function rollout_action(est::BasicPOMCP.SolvedPORollout, pomdp::POMDPs.POMDP, start_state, h::BeliefNode, steps::Int)
    b = extract_belief(est.updater, h)
    sim = BasicPOMCP.RolloutSimulator(est.rng, steps)
    return simulate_action(sim, pomdp, est.policy, est.updater, b, start_state)
end

function extract_belief end

# some defaults are provided
extract_belief(::BasicPOMCP.NothingUpdater, node::BeliefNode) = nothing

function simulate_action(sim::BasicPOMCP.RolloutSimulator, pomdp::POMDP, policy::Policy, updater::Updater, initial_belief, s)
    if sim.eps == nothing
        eps = 0.0
    else
        eps = sim.eps
    end

    if sim.max_steps == nothing
        max_steps = typemax(Int)
    else
        max_steps = sim.max_steps
    end

    disc = 1.0
    r_total = 0.0
    b = initialize_belief(updater, initial_belief)
    step = 1
    action_list = []

    while disc > eps && !isterminal(pomdp, s) && step <= max_steps
        a = action(policy, b)
        push!(action_list, a)
        sp, o, r = @gen(:sp,:o,:r)(pomdp, s, a, sim.rng)
        r_total += disc*r
        s = sp
        bp = update(updater, b, a, o)
        b = bp
        disc *= discount(pomdp)
        step += 1
    end

    return r_total, action_list
end


text:
module TestNARemover

using Test
using Random
using AutoMLPipeline.NARemovers
using AutoMLPipeline.AbsTypes
using DataFrames

function test_naremover()
  Random.seed!(123)
  df  = DataFrame(a=rand([1:3...,missing],100),b=rand([1:9...,missing],100),c=rand([1:20...,missing],100))
  nara = NARemover(0.25)
  @test fit_transform!(nara,df) |> Matrix |> x->sum(skipmissing(x)) == 1546
  narb = NARemover(0.05)
  @test fit_transform!(narb,df) |> Matrix |> x->sum(skipmissing(x)) == 1086
  narc = NARemover(0.01)
  @test fit_transform!(narc,df)  == DataFrame()
end
@testset "NA Remover" begin
  Random.seed!(123)
  test_naremover()
end

end


text:
#!/usr/bin/env julia

#
# Some wrapper function for external maths functions
#
@inline function mittleff(α,z)
    α>0 ? mittlefforiginal(α,z) : 10000000000.
end

@inline function mittleff(α,β,z)
    α>0 ? mittlefforiginal(α,β,z) : 10000000000.
end

#
#  These are not useful at the moment
#
#invLaplace(f::Function, t::Vector{RheoFloat}) = InverseLaplace.talbotarr(f, t)
#invLaplace(f::Function, t::RheoFloat) = InverseLaplace.talbot(f, t)

models_directory = joinpath(@__DIR__, "models")

# spring-pot, spring and dash-pot
include(joinpath(models_directory, "elements.jl"))

# fractional Maxwell model and specialized forms
include(joinpath(models_directory, "maxwell.jl"))

# fractional Kelvin-Voigt model and specialized forms
include(joinpath(models_directory, "kelvinvoigt.jl"))

# fractional Zener model and specialized forms (equivalent to Standard Linear Solid in Maxwell form)
include(joinpath(models_directory, "zener.jl"))

# poynting-thomson model and specialized forms (equivalent to Standard Linear Solid in Kelvin form)
include(joinpath(models_directory, "poynting-thomson.jl"))

# Four-parameters model
include(joinpath(models_directory, "burgers.jl"))

# Empirical
include(joinpath(models_directory, "empirical.jl"))

# Generalized
include(joinpath(models_directory, "generalized.jl"))


text:
using CautiousEngine

@testset "system composition tests" begin

    minPr_mats = [[0.0 0.1; 0.2 0.3], [0.4 0.5; 0.6 0.7]]
    maxPr_mats = [[0.1 0.2; 0.3 0.4], [0.5 0.6; 0.7 0.8]]

    exp_minPr_mat = [0.0 0.1; 0.4 0.5; 0.2 0.3; 0.6 0.7]
    exp_maxPr_mat = [0.1 0.2; 0.5 0.6; 0.3 0.4; 0.7 0.8]

    minPr_mat, maxPr_mat = CautiousEngine.create_switched_system_matrix(minPr_mats, maxPr_mats)

    @test minPr_mat == exp_minPr_mat
    @test maxPr_mat == exp_maxPr_mat
end


text:
# # Making animated and interactive scientific
# # visualizations in Makie.jl

# Alrighty, in this short tutorial I'll teach you how to make
# animated scientific such as this video, or this application.

# Prerequisitives:
# * Basic familiarity with Julia (obviously)
# * Basic familiarity with Makie.jl

# # The process
# 0. Learn about `Observable`s
# 1. Initialize simulation in a stepping manner
# 2. Initialize the `Observable`s of the animation
# 3. Plot the `Observable`s and any other static elements
# 4. Create the "animation stepping function"
# 5. Test it
# 6. Save animations to videos
# 7. Interactive application

# Let's start!
# Load packages that will be used
using DynamicalSystems
using OrdinaryDiffEq
using GLMakie
using DataStructures: CircularBuffer

# %% 0. Learn about `Observable`s
# An `Observable` is a mutable container of an object of type `T`.
# `T` can be any type. The value of an `Observable` can then be
# changed on the spot, just like updating any mutable container.
# (This is similar to the base structure `Ref`, if you're familiar)

# The important part here is that `Observable`s can be "listened" to.
# What does this mean...?

o = Observable(1) # Observable with values of type `Int`

l1 = on(o) do val # Create a listener `l1` of observable.
    println("Observable now has value $val")
end

# `l1` is triggered each time the value of `o` is updated.
# (demo in REPL, set `o[] = 2`.)

# We care about `Observable`s because Makie.jl is hooked up
# to this "listener" system. If any plotted element is
# initialized as an observable, then Makie.jl understands this.
# Updating the observable would update the plot values.

# For example:
ox = 1:4
oy = Observable(rand(4))
lw = Observable(2)

fig, ax = lines(ox, oy; linewidth = lw)
ylims!(ax, 0, 1)

lw[] = 50
oy[] = rand(4)

# This simple process is the basis of creating
# both animations, as well as interactive apps with Makie.jl.
# So in the following let's begin making a simple visualization
# and interactive application of the double pendulum

# %% 1. Initialize simulation in a stepping manner
# (this can also be done with a pre-run simulation)
# the goal is to create a "step" function which
# once called it progresses the data for one animation frame
const L1 = 1.0
const L2 = 0.9
M = 2
u0 = [π/3, 0, 3π/4, -2]
dp = Systems.double_pendulum(u0; L1, L2)

# Solve diffeq with constant step for smoother curves
diffeq = (alg = Tsit5(), adaptive = false, dt = 0.005)

integ = integrator(dp, u0; diffeq...)

function xycoords(state)
    θ1 = state[1]
    θ2 = state[3]
    x1 = L1 * sin(θ1)
    y1 = -L1 * cos(θ1)
    x2 = x1 + L2 * sin(θ2)
    y2 = y1 - L2 * cos(θ2)
    return x1,x2,y1,y2
end

# `integ` is an integrator. `step!(integ)` progresses the integrator
# for one step. `integ.u` is the system state at current step.
# Then `xycoords` converts the states of the integrator
# to their plottable format. So we can imagine something like
function progress_for_one_step!(integ)
    step!(integ)
    u = integ.u
    return xycoords(u)
end
# to be our stepping function that returns the new data.

# If we had finite data instead of a forever-running animation,
# then the "stepping function" would simply be to progress the index `i`
# of the existing data one step forwards...

# %% 2. Initialize the `Observable`s of the animation
# You need to think of this in advance: what things will to be
# animated, and what other plotting elements will be static?
# Animated elements will need to become `Observable`s.

# Here the animated elements will be: balls and rods making the
# double pendulum, and the tail (trajectory) of the pendulum.
x1,x2,y1,y2 = xycoords(u0)
rod   = Observable([Point2f0(0, 0), Point2f0(x1, y1), Point2f0(x2, y2)])
balls = Observable([Point2f0(x1, y1), Point2f0(x2, y2)])
# (Remember: the most optimal way to plot 2D things in Makie.jl is to
# give it a vector of `Point2f0`, the coordinates for the plot)

# Here we have initialized two _different_ observables, because
# rods and balls will be plotted in a different manner (lines/scatter)

# Next is the observable for the tail
tail = 300 # length of plotted trajectory, in units of `dt`
# The circular buffer datastructure makes making stepping-based
# animations very intuitive
traj = CircularBuffer{Point2f0}(tail)
fill!(traj, Point2f0(x2, y2)) # add correct values to the circular buffer
traj = Observable(traj) # make it an observable

# %% 3. Plot the `Observable`s and any other static elements
# Before plotting we need to initialie a figure
fig = Figure(); display(fig)
# in my experience it leads to cleaner code if we first initialize
# an axis and populate it accordingly.
ax = Axis(fig[1,1])

# Now we plot the observables _directly_! First the pendulum
lines!(ax, rod; linewidth = 4, color = :purple)
scatter!(ax, balls; marker = :circle, strokewidth = 2,
    strokecolor = :purple,
    color = :black, markersize = [8, 12]
)

# then its trajectory, with a nice fadeout color
c = to_color(:purple)
tailcol = [RGBAf0(c.r, c.g, c.b, (i/tail)^2) for i in 1:tail]
lines!(ax, traj; linewidth = 3, color = tailcol)

# We can also plot now any other static elements
ax.title = "double pendulum"
ax.aspect = DataAspect()
l = 1.05(L1+L2)
xlims!(ax, -l, l)
ylims!(ax, -l, 0.5l)

# %% 4. Create the "animation stepping function"
# Using the functions of step 1, we now define a function
# that updates the observables. Makie.jl understands observable
# updates and directly reflects this on the plotted elements.
function animstep!(integ, rod, balls, traj)
    x1,x2,y1,y2 = progress_for_one_step!(integ)
    rod[] = [Point2f0(0, 0), Point2f0(x1, y1), Point2f0(x2, y2)]
    balls[] = [Point2f0(x1, y1), Point2f0(x2, y2)]
    push!(traj[], Point2f0(x2, y2))
    traj[] = traj[] # <- important! Updating in-place the value of an
                    # `Observable` does not trigger an update!
end

# %% 5. Test it
for i in 1:1000
    animstep!(integ, rod, balls, traj)
    sleep(0.001)
end

# cool it works. Let's wrap up the creation of the observables
# and plots in a function (just to re-initialie everything)
function makefig(u0)
    dp = Systems.double_pendulum(u0; L1, L2)
    integ = integrator(dp, u0; diffeq...)
    x1,x2,y1,y2 = xycoords(u0)
    rod   = Observable([Point2f0(0, 0), Point2f0(x1, y1), Point2f0(x2, y2)])
    balls = Observable([Point2f0(x1, y1), Point2f0(x2, y2)])
    traj = CircularBuffer{Point2f0}(tail)
    fill!(traj, Point2f0(x2, y2)) # add correct values to the circular buffer
    traj = Observable(traj) # make it an observable
    fig = Figure(); display(fig)
    ax = Axis(fig[1,1])
    lines!(ax, rod; linewidth = 4, color = :purple)
    scatter!(ax, balls; marker = :circle, strokewidth = 2,
        strokecolor = :purple,
        color = :black, markersize = [8, 12]
    )
    lines!(ax, traj; linewidth = 3, color = tailcol)
    ax.title = "double pendulum"
    ax.aspect = DataAspect()
    l = 1.05(L1+L2)
    xlims!(ax, -l, l)
    ylims!(ax, -l, 0.5l)
    # also return the figure object, we'll ues it!
    return fig, integ, rod, balls, traj
end

# %% 6. Save animations to videos
fig, integ, rod, balls, traj = makefig(u0)
frames = 1:200
record(fig, "video.mp4", frames; framerate = 60) do i # i = frame number
    for j in 1:5 # step 5 times per frame
        animstep!(integ, rod, balls, traj)
    end
    # any other manipulation of the figure here...
end # for each step of this loop, a frame is recorded

# %% 7. Interactive application
# Makie.jl has tremendously strong capabilities for real-time
# interactivity. To learn all of this takes time of course,
# and you'll need to consult the online documentation.
# Here we will do two interactions: 1) a play/stop button
# 2) clicking on the screen and getting a new initial condition!

fig, integ, rod, balls, traj = makefig(u0)
# The run button is actually pretty simple, we'll add it below the plot
run = Button(fig[2,1]; label = "run", tellwidth = false)
# This button will start/stop an animation. It's actually surprisingly
# simple to do this. The magic code is:
isrunning = Observable(false)
on(run.clicks) do clicks; isrunning[] = !isrunning[]; end
on(run.clicks) do clicks
    @async while isrunning[]
        isopen(fig.scene) || break # ensures computations stop if closed window
        animstep!(integ, rod, balls, traj)
        sleep(0.001) # or `yield()` instead
    end
end

# `on` an important Observables function when one starts
# doing advanced stuff. It triggers a piece of code once an observable
# triggers its update.

# We'll add one more interactive feature which will trigger once
# we click on the axis. Notice that by default makie performs a zoom
# once one clicks on the axis, so we'll disable this
ax = content(fig[1,1])
Makie.deactivate_interaction!(ax, :rectanglezoom)
# and we'll add a new trigger using the `select_point` function:
spoint = select_point(ax.scene)

# Now we see that we can click on the screen and the `spoint` updates!

# Okay, let's do something useful when it triggers
function θωcoords(x, y)
    θ = atan(y,x) + π/2
    return SVector(θ,0,0,0)
end

on(spoint) do z
    x, y = z
    u = θωcoords(x, y)
    reinit!(integ, u)
    # Reset tail and balls to new coordinates
    x1,x2,y1,y2 = xycoords(u)
    traj[] .= fill(Point2f0(x2, y2), length(traj[]))
    traj[] = traj[]
    rod[] = [Point2f0(0, 0), Point2f0(x1, y1), Point2f0(x2, y2)]
    balls[] = [Point2f0(x1, y1), Point2f0(x2, y2)]
end

# And that's the end of the tutorial!

max_stars_repo_name: chris-revell/ClockOscillators

max_issues_repo_name: chris-revell/ClockOscillators

max_forks_repo_name: chris-revell/ClockOscillators


text:
##the following code is modified from the paper provided : https://doi.org/10.1038/s41588-020-0687-1

using Distributions, StatsBase, DataFrames, GLM, CSV, DelimitedFiles

function getFitness(n)
    (1 + s*n)
end

mutable struct cancercell
    mutations::Array{Int64,1}
    mut_neop::Array{Bool,1}
    fitness::Float64
    neonumber::Float64
    escaped::Array{Bool,1}
end

function newmutations(cancercell, mutID, p, pesc)
    cancercell.mutations = append!(cancercell.mutations, mutID)
    mutID = mutID + 1

    neoep = rand()<p
    if neoep
        cancercell.neonumber = cancercell.neonumber + 1
        cancercell.fitness = getFitness(cancercell.neonumber) #fitness is affected by the number of mutations
        cancercell.mut_neop = append!(cancercell.mut_neop,true)
    else
        cancercell.mut_neop = append!(cancercell.mut_neop,false)
    end

    mut_escape = rand() < pesc
    if mut_escape
        cancercell.escaped = append!(cancercell.escaped,true)
    else
        cancercell.escaped = append!(cancercell.escaped,false)
    end

    return cancercell, mutID
end

function copycell(cancercellold::cancercell)
  newcancercell::cancercell = cancercell(copy(cancercellold.mutations), copy(cancercellold.mut_neop), copy(cancercellold.fitness), copy(cancercellold.neonumber), copy(cancercellold.escaped))
end

function start_population(p, pesc, initial_mut)
    mutID = 1
    N = 1
    cells = cancercell[]
    #muts = Dict{Int64, Float64}()
    push!(cells,cancercell([],[],1,0,[]))
    for i=1:initial_mut
        cells[1],mutID = newmutations(cells[1],mutID, p, pesc)
        #muts[mutID-1] = cells[1].mut_neop[1]
    end

    return cells, mutID,  N

end

function birthdeath_neoep(b0, d0, Nmax, p, initial_mut, mu, pesc)

    dmax = d0 #dmax is updated throughout, starts from d0

    #initialize arrays and parameters
    cells, mutID, N = start_population(p, pesc, initial_mut )
    Nvec = Int64[]
    push!(Nvec,N)
    t = 0.0
    tvec = Float64[]
    push!(tvec,t)

    while (N < Nmax) & (t < 300) #set so we can exit simulation where there is a lot of death

        #pick a random cell
        randcell = rand(1:N)
        Nt = N

        #a cell's immunogenicity depends on its fitness, i.e. the summed antigenicity of neoepitopes
        d = max(0, (d0 - b0)*cells[randcell].fitness + b0)

        if true in cells[randcell].escaped # discard effect of fitness if the cell escape
            d = d0
        end

        if (d > dmax) #update dmax to keep track of the highest death rate in the whole population
            dmax = d
        end

        Rmax = b0 + dmax

        r = rand(Uniform(0,Rmax)) #Pick which reaction should happen to cell

        # If r < birthrate, a birth event happens: a new cell is created and randcell updated as a new one
        if r < b0

            #population increases by one
            N = N + 1
            #copy cell and mutations for cell that reproduces
            push!(cells, copycell(cells[randcell]))

            #add new mutations to both new cells, the number of mutations is Poisson distributed
            for i=1:(rand(Poisson(mu)))
                cells[randcell],mutID = newmutations(cells[randcell],mutID, p, pesc)
                #muts[mutID-1] = cells[randcell].mut_neop[1]
            end
            for i=1:(rand(Poisson(mu)))
                cells[end],mutID = newmutations(cells[end],mutID, p, pesc)
                #muts[mutID-1] = cells[end].mut_neop[1]
            end

            push!(Nvec, N)
            Δt =  1/(Rmax * Nt) .* - log(rand())
            t = t + Δt
            push!(tvec,t)

        end

        #if r has neither birth or death (only possible if it is a non-immunogenic cell), nothing happens
        if  (b0+d)<= r
          push!(Nvec, N)
          Δt =  1/(Rmax * Nt) .* - log(rand())
          t = t + Δt
          push!(tvec,t)
        end

        #death event if r > b but < d
        if b0 <= r < (b0+d)

            #population decreases by 1, overall fitness score also decreases if it was non-zero
            N = N - 1

            #remove deleted cell
            deleteat!(cells,randcell)
            push!(Nvec,N)
            Δt =  1/(Rmax * Nt) .* - log(rand())
            t = t + Δt
            push!(tvec,t)
        end

        #if every cell dies, restart simulation from a single cell again
        if (N == 0)
            cells, mutID, N = start_population(p, pesc, initial_mut)
            push!(Nvec,N)
            push!(tvec,t)
        end

    end

    return Nvec, tvec, mutID, cells
end

function process_mutations(cells)
    mutVec = []
    mut_is_neo = []
    mut_is_escape = []
    for i=1:length(cells)
        append!(mutVec, cells[i].mutations)
        append!(mut_is_neo, cells[i].mut_neop)
        append!(mut_is_escape, cells[i].escaped)

    end

    #detMutDict = filter((k, v) -> v > detLim, countmap(mutVec))
    detMutDict = countmap(mutVec)

    println("Mutations processed for ", length(cells), " cells.")
    return detMutDict, mutVec, mut_is_neo, mut_is_escape

end

b0 = 1
d0 = 0.1
popSize = 1e5
p = 0.1
initial_mut = 0
mu = 1
pesc = 1e-6

s=<s>
for i=1:200
    Nvec, tvec, mutID, cells = birthdeath_neoep(b0, d0, popSize, p, initial_mut, mu,pesc);
    outNDFsim = DataFrame(t=tvec, N=Nvec)
    CSV.write("/public/slst/home/wutao2/julia_simulation/out/minus_"*string(-s)*"/preIT_"*string(i)*".txt", outNDFsim)

    detMutDict, mutVec, mut_is_neo, mut_is_escape = process_mutations(cells)
    writedlm("/public/slst/home/wutao2/julia_simulation/out/minus_"*string(-s)*"/vaf_preIT_"*string(i)*".txt",detMutDict)

    mut_summ = DataFrame(mut=mutVec,mut_neo=mut_is_neo,mut_escaped=mut_is_escape)
    CSV.write("/public/slst/home/wutao2/julia_simulation/out/minus_"*string(-s)*"/mutsumm_"*string(i)*".txt", mut_summ)

end

max_stars_repo_name: XSLiuLab/Immunoediting

max_issues_repo_name: XSLiuLab/Immunoediting

max_forks_repo_name: XSLiuLab/Immunoediting


text:
function part1()
    lines = readlines("inputD3.txt")
    gammastr = ""
    for idx in 1:length(lines[1])
        zeroes = sum([number[idx] == '0' for number in lines])
        ones = length(lines) - zeroes
        if ones > zeroes
            gammastr *= "1"
        else
            gammastr *= "0"
        end
    end

    gamma = parse(Int, gammastr, base = 2)
    epsilon = Int(2 ^ ceil(log2(gamma)) - 1 - gamma)
    return gamma * epsilon
end

function part2()
    lines = readlines("inputD3.txt")
    oxyfiltered = copy(lines)
    co2filtered = copy(lines)
    for i in 1:length(lines[1])
        if length(oxyfiltered) > 1
            zeroes = sum([number[i] == '0' for number in oxyfiltered])
            ones = length(oxyfiltered) - zeroes
            if ones > zeroes || ones == zeroes
                oxyfiltered = [num for num in oxyfiltered if num[i] == '1']
            elseif ones < zeroes
                oxyfiltered = [num for num in oxyfiltered if num[i] == '0']
            end
        end

        if length(co2filtered) > 1
            co2zeroes = sum([number[i] == '0' for number in co2filtered])
            co2ones = length(co2filtered) - co2zeroes
            if co2ones > co2zeroes || co2ones == co2zeroes
                co2filtered = [num for num in co2filtered if num[i] == '0']
            elseif co2ones < co2zeroes
                co2filtered = [num for num in co2filtered if num[i] == '1']
            end
        end
    end

    if length(co2filtered) != 1 || length(oxyfiltered) != 1
        error("???")
    end

    return parse(Int, oxyfiltered[1], base = 2) * parse(Int, co2filtered[1], base = 2)
end

println(part1())
println(part2())


text:
type DEAData <: AbstractArray{Any,1}
  X::Array
  Y::Array
  gx::Array
  gy::Array
  indexes::Vector{Int}

  function DEAData(X::Array, Y::Array, gx::Array = [], gy::Array = []; indexes::Vector{Int} = Vector{Int}(collect(1:size(X,1))))
    if(size(X,1) != size(Y,1))
			error("Number of DMUs in X and Y are not the same: X has $(size(X,1)) DMUs and Y has $(size(Y,1)) DMUs!")
		end
    if(!isempty(gx) && size(gx) != size(X))
      error("Dimensions of gx and X are not the same: size(gx) = $(size(gx)), size(X) = $(size(X))!")
    end
    if(!isempty(gy) && size(gy) != size(Y))
      error("Dimensions of gy and Y are not the same: size(gy) = $(size(gy)) and size(Y) = $(size(Y))!")
    end
    if(isempty(gx) != isempty(gy))
      error("If gx is empty then gy must also be empty (and vice versa)!")
    end
    if(!isempty(gx) && !isempty(gy) && gx == zeros(size(gx)) && gy == zeros(size(gy)))
      error("Direction vectors cannot both be zero!")
    end
    new(X,Y,gx,gy,indexes)
  end
end

function getnrdmu(Data::DEAData)
  return size(Data,1)
end

function getiodim(Data::DEAData)
  return size(Data.X,2), size(Data.Y,2)
end

function getindexes(Data::DEAData)
  return copy(Data.indexes)
end

function setindexes!(Data::DEAData, newindexes::Vector{Int})
  if(length(newindexes) > size(Data.X,1))
    error("$(length(newindexes)) > $(size(Data.X,1))!")
  end
  Data.indexes = copy(newindexes)
end

Base.size(Data::DEAData) = (length(Data.indexes),)
Base.IndexStyle(::Type{DEAData}) = IndexLinear()

function Base.getindex(Data::DEAData,i::Int)
  1 <= i <= size(Data,1) || throw(BoundsError(Data, i))
  if(!isempty(Data.gx))
    return Data.X[Data.indexes[i],:],Data.Y[Data.indexes[i],:],Data.gx[Data.indexes[i],:],Data.gy[Data.indexes[i],:]
  else
    return Data.X[Data.indexes[i],:],Data.Y[Data.indexes[i],:]
  end
end

function Base.getindex(Data::DEAData, I)
  N,M = getiodim(Data)
  K = length(I)
  X = Array{Float64}(K,N)
  Y = Array{Float64}(K,M)
  if(!isempty(Data.gx))
    gx = Array{Float64}(K,N)
    gy = Array{Float64}(K,M)
    for i in eachindex(I)
      X[i,:],Y[i,:],gx[i,:],gy[i,:] = Data[I[i]]
    end
    return X,Y,gx,gy
  else
    for i in eachindex(I)
      X[i,:],Y[i,:] = Data[I[i]]
    end
    return X,Y
  end
end

max_stars_repo_name: kepiej/NEPA


text:
#Newtonian eqn for impurity
function NewtonImp(Φ::Array, V::Array{<:Array{Float64}},  M::Float64, λ::Float64, Δx::Float64)
    d = length(size(Φ))
    ans = zeros(d)
    for i = 1:d
        ans[i] =1.0 / M *  λ * real(sum(V[i+1] .* conj(Φ) .* Φ) * Δx)
    end
    return ans
end

max_stars_repo_name: phyjonas/ColloidFluctuatingField


text:
#######
# Here we store the parameters corresponding to the non-SQUID-like geometry (`mask = false`)
######

p = Params(scale = 10, Δ=.0, Ln=200, Ls=0, W = 320, U = 0.0005, μSL = .8, μSR=.8,
    μN = 0*-0.0034, V_perp =  0., V_par = 1e-6, B_par = 0, B_perp = 0., Vu = 0.,
    edges = :armchair_zigzag, smooth = false, dir = :x, U_dimer = 0.0, α = 0., smoothμ = false,
    localselfenergy = false, analyticalselfenergy = false, normal=true);

p_λarm  = reconstruct(p, λ = 0.03, edges = :armchair_zigzag) # edge regime
p_λzig  = reconstruct(p, λ = 0.03, edges = :zigzag_armchair) # edge regime

p_λarm_gap  = reconstruct(p, λ = 0.03, α = 0.02, V_par = 0.005,
    edges = :armchair_zigzag) # gapped regime
p_λzig_gap  = reconstruct(p, λ = 0.03, α = 0.02, V_par = 0.005,
    edges = :zigzag_armchair) # gapped regime

p_λarm_bulk = reconstruct(p, λ = 0.03, α = 0.02, V_par = 0.005,
    edges = :armchair_zigzag, μN = .1) # bulk regime
p_λzig_bulk = reconstruct(p, λ = 0.03, α = 0.02, V_par = 0.005,
    edges = :zigzag_armchair, μN = .1) # bulk regime

p_λarm_gap_ls  = reconstruct(p_λarm_gap, μSL = 0.00001, μSR = 0.00001,
    Ls = 100, Δ = 1, normal = false) # gapped regime
p_λarm_bulk_ls  = reconstruct(p_λarm_bulk, μSL = 8, μSR = 8,
    Ls = 100, Ln = 0, Δ = 1, normal = false) # bulk regime
p_λarm_ls  = reconstruct(p_λarm, μSL = 0.00001, μSR = 0.00001,
    Ls = 100, Δ = 1, normal = false) # edge regime
pa_inv_smoothμ = reconstruct(p, smoothμ = true, λ = 0.02,
    Ls = 0, Ln = 200, W = 200, scale = 20, Δxμ = 5, λμ = 5, λμbarrier = 2, Δμbarrier = -2)

p_test2 = reconstruct(p_λarm_gap_ls, scale = 1, Ls = 1, Ln = 2, W = 5)
p_test4 = reconstruct(p_test2, scale = 1, Ls = 2, Ln = 4, W = 20, μN = .5);
p_test5 = reconstruct(p_test2, scale = 1, Ls = 4, Ln = 20, W = 20, μN = .5);

pa_inv_smoothμ_ls = reconstruct(pa_inv_smoothμ, Ls = 50, μSL = .8,
     μSR = .8, Δ = 0.004, normal = false)
pa_inv_smoothμ_ls_trivgap = reconstruct(pa_inv_smoothμ, smoothμ = false, U = 0.04)
pa_inv_smoothμ_ls_bulkconduction =
    reconstruct(pa_inv_smoothμ_ls, smoothμ = false, μN = 0.015)

####### Flux to B conversion
Φ0 = 2067.8338 # T*nm^2
function Bprange(p, Φrange)
    Bperp = collect(Φrange./(p.Ln * p.W)) #conversion from flux to Bfield (size dependent)
    Bcriteria(Bperp, p.scale)
    return Bperp
end

"We need to make sure that we are not in the Hofstadter regime,
    for the scaled and non scaled system"
function Bcriteria(B::Array, scale; threshold = 50.)
    ħc = ustrip(u"eV*nm",ħ*c_0)
    a0 = 0.246    #nm
    Bmax = findmax(abs.(B))[1]
    l_B = ħc/Bmax
    l_scale = scale * a0
    println(l_B/ l_scale)
    println(l_B/a0)
    l_B/ l_scale > threshold ? nothing : @warn("lB ≈ a0*scale, decrease scaling")
    l_B/ a0 > threshold ? nothing : @warn("lB ≈ a0, increase BLG area
        (we are in the Hofstadter regime)")
end


text:
#
# Plot orbits
#

"""
Plots every timestep in `sols` in `3D` space. All keyward
arguments are passed directly to `Plots.jl`.
"""
function orbitplot(sols::Trajectory{<:NBodySystem}; bodies=1:length(sols.step[1].body), kwargs...)

    # Referencing:
    # [1] https://discourse.julialang.org/t/smart-kwargs-dispatch/14571/15

    # Set default kwargs (modified from [1])
    defaults = (;   formatter=:scientific,
                    legend=:topleft,
                    xlabel="X Position (km)",
                    ylabel="Y Position (km)",
                    zlabel="Z Position (km)",
                    title ="NBody Positions vs. Time")
    options = merge(defaults, kwargs)

    fig = Plots.plot()
    for i = bodies

        Plots.plot!(fig, ustrip.(u"km", map(x -> x.body[i].r̅[1], sols.step)),
                         ustrip.(u"km", map(x -> x.body[i].r̅[2], sols.step)),
                         ustrip.(u"km", map(x -> x.body[i].r̅[3], sols.step)),
                         label=string("Body ", i))

    end
    Plots.plot!(fig; options...)

    return fig

end