text:
# Basically an augmented LogBinner
struct LightObservable{T, N} <: AbstractObservable
    B::LogBinner{T,N}

    # parameters (external)
    name::String
    outfile::String
    group::String # where to put data in JLD file
end

LightObservable(::Type{T}; kwargs...) where T = LightObservable(zero(T); kwargs...)

function LightObservable(x::T;
                    name::String="unnamed",
                    alloc::Integer=1000,
                    outfile::String="Observables.jld",
                    group::String=name) where T <: Union{Number, AbstractArray}

    B = LogBinner(x, capacity=alloc)
    lo = LightObservable(B,
        name,
        outfile,
        group
    )
    return lo
end

@forward LightObservable.B (Base.length, Base.ndims, Base.push!,Base.append!,
                             Base.empty!, Base.isempty, Base.eltype,
                             Statistics.mean, Statistics.var, BinningAnalysis.std_error)

@inline name(obs::LightObservable) = obs.name
@inline inmemory(obs::LightObservable) = true
@inline Base.size(obs::LightObservable) = @inbounds size(obs.B.x_sum[1])

function _print_header(io::IO, obs::LightObservable{T,N}) where {T,N}
    print(io, "LightObservable{$(T),$(N)}")
    nothing
end

function _println_body(io::IO, obs::LightObservable{T,N}) where {T,N}
    n = length(obs)
    println(io)
    println(io, "| Name: ", name(obs))
    print(io, "| Measurements: ", n)
    if n > 0 && ndims(obs) == 0
        print(io, "\n| Mean: ", round.(mean(obs), digits=5))
        print(io, "\n| StdError: ", round.(std_error(obs), digits=5))
    end
    nothing
end

# short version (shows up in arrays etc.)
Base.show(io::IO, obs::LightObservable{T,N}) where {T,N} = print(io, "LightObservable{$(T),$(N)}()")
# verbose version (shows up in the REPL)
Base.show(io::IO, m::MIME"text/plain", obs::LightObservable) = (_print_header(io, obs); _println_body(io, obs))

import Base: error
@deprecate error(obs::LightObservable) std_error(obs)

BinningAnalysis.tau(obs::LightObservable) = BinningAnalysis.tau(obs.B)

# Assure var(obs) == var(data) and same for std
Statistics.var(obs::LightObservable) = var(obs, 1)
Statistics.std(obs::LightObservable) = map(sqrt, var(obs))

# -------------------------------------------------------------------------
#   Exporting results
# -------------------------------------------------------------------------
"""
    export_results(obs::LightObservable[, filename::AbstractString, group::AbstractString])

Export result for given observable nicely to JLD.

Will export name, number of measurements, estimates for mean and standard error.
"""
function export_result(obs::LightObservable, filename::AbstractString=obs.outfile,
                         group::AbstractString=obs.group*"_export"; error=true)
    grp = endswith(group, "/") ? group : group*"/"

    jldopen(filename, isfile(filename) ? "r+" : "w") do f
        !HDF5.has(f.plain, grp) || delete!(f, grp)
        write(f, joinpath(grp, "name"), name(obs))
        write(f, joinpath(grp, "count"), length(obs))
        m = mean(obs)
        write(f, joinpath(grp, "mean"), m)
        if error
            err = std_error(obs)
            write(f, joinpath(grp, "error"), err)
            write(f, joinpath(grp, "error_rel"), abs.(err ./ m))
        end
    end
    nothing
end

"""
    export_error(obs::LightObservable[, filename::AbstractString, group::AbstractString;])

Export standard error estimate.
"""
function export_error(obs::LightObservable, filename::AbstractString=obs.outfile,
                         group::AbstractString=obs.group)
    grp = endswith(group, "/") ? group : group*"/"

    jldopen(filename, isfile(filename) ? "r+" : "w") do f
        !HDF5.has(f.plain, grp*"error") || delete!(f, grp*"error")
        !HDF5.has(f.plain, grp*"error_rel") || delete!(f, grp*"error_rel")
        err = std_error(obs)
        write(f, joinpath(grp, "error"), err)
        write(f, joinpath(grp, "error_rel"), abs.(err ./ mean(obs)))
    end
    nothing
end


text:
### Abstract variable NStates

abstract VariableNState{F<:VariateForm}

add_dimension(n::Number) = eltype(n)[n]
add_dimension(a::Array, sa::Tuple=size(a)) = reshape(a, sa..., 1)

Base.writemime(io::IO, ::MIME"text/plain", nstate::VariableNState) = show(io, nstate)

### Basic variable NState subtypes

## BasicUnvVariableNState

type BasicUnvVariableNState{N<:Number} <: VariableNState{Univariate}
  value::Vector{N}
  n::Integer
end

BasicUnvVariableNState{N<:Number}(value::Vector{N}) = BasicUnvVariableNState{N}(value, length(value))

BasicUnvVariableNState{N<:Number}(n::Integer, ::Type{N}=Float64) = BasicUnvVariableNState{N}(Array(N, n), n)

Base.eltype{N<:Number}(::Type{BasicUnvVariableNState{N}}) = N
Base.eltype{N<:Number}(::BasicUnvVariableNState{N}) = N

Base.copy!(nstate::BasicUnvVariableNState, state::BasicUnvVariableState, i::Integer) = (nstate.value[i] = state.value)

function Base.show{N<:Number}(io::IO, nstate::BasicUnvVariableNState{N})
  indentation = "  "

  println(io, "BasicUnvVariableNState:")

  println(io, indentation*"eltype: $(eltype(nstate))")
  print(io, indentation*"number of states = $(nstate.n)")
end

Base.writemime{N<:Number}(io::IO, ::MIME"text/plain", nstate::BasicUnvVariableNState{N}) = show(io, nstate)

## BasicMuvVariableNState

type BasicMuvVariableNState{N<:Number} <: VariableNState{Multivariate}
  value::Matrix{N}
  size::Integer
  n::Integer
end

BasicMuvVariableNState{N<:Number}(value::Matrix{N}) = BasicMuvVariableNState{N}(value, size(value)...)

BasicMuvVariableNState{N<:Number}(size::Integer, n::Integer, ::Type{N}=Float64) =
  BasicMuvVariableNState{N}(Array(N, size, n), size, n)

Base.eltype{N<:Number}(::Type{BasicMuvVariableNState{N}}) = N
Base.eltype{N<:Number}(::BasicMuvVariableNState{N}) = N

Base.copy!(nstate::BasicMuvVariableNState, state::BasicMuvVariableState, i::Integer) =
  (nstate.value[1+(i-1)*state.size:i*state.size] = state.value)

function Base.show{N<:Number}(io::IO, nstate::BasicMuvVariableNState{N})
  indentation = "  "

  println(io, "BasicMuvVariableNState:")

  println(io, indentation*"eltype: $(eltype(nstate))")
  println(io, indentation*"state size = $(nstate.size)")
  print(io, indentation*"number of states = $(nstate.n)")
end

Base.writemime{N<:Number}(io::IO, ::MIME"text/plain", nstate::BasicMuvVariableNState{N}) = show(io, nstate)

## BasicMavVariableNState

type BasicMavVariableNState{N<:Number} <: VariableNState{Matrixvariate}
  value::Array{N, 3}
  size::Tuple{Integer, Integer}
  n::Integer
end

BasicMavVariableNState{N<:Number}(value::Array{N, 3}) =
  BasicMavVariableNState{N}(value, (size(value, 1), size(value, 2)), size(value, 3))

BasicMavVariableNState{N<:Number}(size::Tuple, n::Integer, ::Type{N}=Float64) =
  BasicMavVariableNState{N}(Array(N, size..., n), size, n)

Base.eltype{N<:Number}(::Type{BasicMavVariableNState{N}}) = N
Base.eltype{N<:Number}(::BasicMavVariableNState{N}) = N

Base.copy!(nstate::BasicMavVariableNState, state::BasicMavVariableState, i::Integer, statelen::Integer=prod(state.size)) =
  (nstate.value[1+(i-1)*statelen:i*statelen] = state.value)

function Base.show{N<:Number}(io::IO, nstate::BasicMavVariableNState{N})
  indentation = "  "

  println(io, "BasicMavVariableNState:")

  println(io, indentation*"eltype: $(eltype(nstate))")
  println(io, indentation*"state size = $(nstate.size)")
  print(io, indentation*"number of states = $(nstate.n)")
end

ext: jl


text:
using ThreeBodyProblem
using Test
using LinearAlgebra
using SPICE
furnsh("../src/kernels/de440s.bsp")
furnsh("../src/kernels/naif0012.tls")

err_tol = 1e-6

# Altair Meteor data
r_enu = [191.0973115173407; 78.7745585024403; 57.3291938159787] # [km]
v_enu = [-52.199447002937305; -21.517740337533645; -15.659828241680069] #[km/s]

rv_enu = [r_enu; v_enu]

# coordinates of observer (ALTAIR radar)
ϕ = 9.3987 # [deg] latitude
λ = 167.482 # [deg] longitude
h = 0.083532 # [km] altitude

# verified against matlab stuff 4/10/21
rv_ecef = enu2ecef(rv_enu, ϕ, λ, h, geodetic=true, ang_unit=:deg)

# Matlab result
@test_broken norm(rv_ecef - 1.0e+03*[-6.2281, 1.1870, 1.1219, 0.0230, 0.0484,-0.0238]) < err_tol
# Reversible function
@test norm(rv_enu - ecef2enu(rv_ecef, ϕ, λ, h)) < err_tol

# UT1 November 17, 2000 15:00:05.634
date = [2000,11,17,15,0,5.64]
mjd = date2mjd(date)
θ = mjd2gmst(mjd)

rv_eci = ecef2eci(rv_ecef, θ, ang_unit=:rad)
@test norm(rv_ecef - eci2ecef(rv_eci, θ)) < err_tol

sys = sun_earth()
a,e,i,Ω,ω,ν = cart2oe(rv_eci, 3.986004328969392e+05)
@test_broken norm(rv_eci - oe2cart(a,e,i,Ω,ω,ν,G*sec.m)) < err_tol

# Convert the calendar date to ephemeris seconds past J2000
et = utc2et(date2str(date))
# Get the position of Mars at `et` w.r.t. Earth
rv_sun_eci = spkezr("sun", et, "J2000", "none", "earth")[1]
rv_sci = eci2sci(rv_eci, rv_sun_eci)
@test norm(rv_eci - sci2eci(rv_sci, rv_sun_eci)) < err_tol

# Test started on April 22, 2021 I want to see a circle
sys = sun_earth()
Lpts = computeLpts(sys)
tt = 0:π/10:2π
xx = [[Lpts[1];zeros(3)] for t in tt]
# plot(sys)
# plot!(xx)
xx_inert = [rot2inert(xx[i],tt[i],sys) for i = 1:length(tt)]
@test norm(xx - [inert2rot(xx_inert[i],tt[i],sys) for i = 1:length(tt)]) < 1e-12


text:
using Clang.Generators
using SSGraphBLAS_jll

cd(@__DIR__)

const HEADER_BASE = joinpath(SSGraphBLAS_jll.artifact_dir, "include")
const GRAPHBLAS = joinpath(HEADER_BASE, "GraphBLAS.h")

headers = [GRAPHBLAS]

options = load_options(joinpath(@__DIR__, "generator.toml"))
args = get_default_args()
push!(args, "-I$HEADER_BASE")

ctx = create_context(headers, args, options)

build!(ctx)


text:
################################
####### Define parameters ######
################################
import Parameters: @with_kw

"""
Mutable structure containing the variables required to solve the analytical approach. All the functions are solve using the internal values of the structure. You should declare a mutable structure to the perform the analytical estimations.

# Parameters
 - `gamNeg::Int64`: Selection coefficient for deleterious alleles
 - `gL::Int64`: Selection coefficient for weakly benefical alleles
 - `gH::Int64`: Selection coefficient for strongly benefical alleles
 - `alLow::Float64`: Proportion of α due to weak selection
 - `alTot::Float64`: α
 - `thetaF::Float64`: Mutation rate defining BGS strength
 - `thetaMidNeutral::Float64`: Mutation rate on coding region
 - `al::Float64`: DFE shape parameter
 - `be::Float64`: DFE scale parameter
 - `B::Float64`: BGS strength
 - `B_bins::Array{Float64,1}`: BGS values to simulate
 - `pposL::Float64`: Fixation probabily of weakly beneficial alleles
 - `pposH::Float64`: Fixation probabily of strongly beneficial alleles
 - `N::Int64`: Population size
 - `n::Int64`: Sample size
 - `Lf::Int64`: Flanking region length
 - `rho::Float64`: Recombination rate
 - `TE::Float64`

"""
@with_kw mutable struct parameters
	gamNeg::Int64              = -457
	gL::Int64                  = 10
	gH::Int64                  = 500
	alLow::Float64             = 0.2
	alTot::Float64             = 0.4
	thetaF::Float64            = 1e-3
	thetaMidNeutral::Float64   = 1e-3
	al::Float64                = 0.184
	be::Float64                = 0.000402
	B::Float64                 = 0.999
	B_bins::Array{Float64,1}   = push!(collect(0.1:0.025:0.975),0.999)
	pposL::Float64             = 0
	pposH::Float64             = 0.001
	N::Int64                   = 1000
	n::Int64                   = 500
	Lf::Int64                  = 2*10^5
	rho::Float64               = 0.001
	TE::Float64                = 5.0
	diploid::Bool              = false

	NN::Int64 = 2*N
	nn::Int64 = 2*n
	θᵣ::Array{Float64,1} = fill(thetaMidNeutral,nn-1)
	dac::Array{Int64,1} = [2,4,5,10,20,50,200,500,700]
end

"""
Mutable structure containing the downsampled SFS.

# Returns
 - `bn::Dict`: SparseMatrixCSC containing the binomial convolution

"""
@with_kw mutable struct binomialDict
	bn::Dict = Dict{Float64,SparseMatrixCSC{Float64,Int64}}()
end
#=const binomialDict = Dict{Float64,SparseMatrixCSC{Float64,Int64}}()=#

################################
###### Solving parameters ######
################################
"""
	Br(Lmax,theta)

Expected reduction in nucleotide diversity. Explored at [Charlesworth B., 1994](https://doi.org/10.1017/S0016672300032365):

```math
\\frac{\\pi}{\\pi_{0}} = e^{\\frac{4\\muL}{2rL+t}}
```

# Arguments
 - `param::parameters`
 - `theta::Float64`
# Returns
 - `Float64`: expected reduction in diversity given a non-coding length, mutation rate and defined recombination.
"""
function Br(param::parameters,theta::Float64)

	ρ::Float64     = param.rho
	t::Float64     = -1.0*param.gamNeg/(param.NN+0.0)
	μ::Float64     = theta/(2.0*param.NN)
	r::Float64     = ρ/(2.0*param.NN)

	out::Float64  = ℯ^(-4*μ*param.Lf/(2*param.Lf*r+t))
	return out
end

# Set mutation rate given the expected reduction in nucleotide diversity (B value ) in a locus.
"""
	setThetaF!(param)

Find the optimum mutation given the expected reduction in nucleotide diversity (B value) in a locus.

# Returns
 - `adap.thetaF::Float64`: changes adap.thetaF value.
"""
function setThetaF!(param::parameters)

	i(θ,p=param) = Br(p,θ)-p.B
	thetaF       = find_zero(i,0.0)
	param.thetaF = thetaF
end

function alphaExpSimLow(param::parameters,pposL::Float64,pposH::Float64)
	return fixPosSim(param,param.gL,0.5*pposL)/(fixPosSim(param,param.gL,0.5*pposL)+fixPosSim(param,param.gH,0.5*pposH) + fixNegB(param,0.5*pposL+0.5*pposH))
end

function alphaExpSimTot(param::parameters,pposL::Float64,pposH::Float64)
	return (fixPosSim(param,param.gL,0.5*pposL)+fixPosSim(param,param.gH,0.5*pposH))/(fixPosSim(param,param.gL,0.5*pposL)+fixPosSim(param,param.gH,0.5*pposH)+fixNegB(param,0.5*pposL+0.5*pposH))
end

function solvEqns(param,config)

	pposL,pposH = config
	return (alphaExpSimTot(param,pposL,pposH)-param.alTot,alphaExpSimLow(param,pposL,pposH)-param.alLow)
end

"""
	setPpos!(param)

Find the probabilty of positive selected alleles given the model. It solves a equation system taking into account fixations probabilities of weakly and strong beneficial alleles.

# Returns
 - `Tuple{Float64,Float64}`: weakly and strong beneficial alleles probabilites.
"""
function setPpos!(param::parameters)

	function f!(F,x,param=param)
		F[1] = alphaExpSimTot(param,x[1],x[2])-param.alTot
		F[2] = alphaExpSimLow(param,x[1],x[2])-param.alLow
	end

	pposL,pposH = nlsolve(f!,[0.0; 0.0]).zero

	if pposL < 0.0
		 pposL = 0.0
	end
	if pposH < 0.0
		 pposH = 0.0
	end

	param.pposL,param.pposH = pposL, pposH

 end

"""
	binomOp(param)

Binomial convolution to sample the allele frequencies probabilites depending on background selection values and sample size.

# Arguments
 - `param::parameters`
 - `convolutedSamples::binomialDict`

# Returns
 - `Array{Float64,2}`: convoluted SFS for each B value defined in the model (param.B_bins). The estimations are saved at *convolutedBn.bn*.
"""
function binomOp!(param::parameters,convolutedBn::Dict)

	#=bn = Dict(param.B_bins[i] => spzeros(param.nn+1,param.NN) for i in 1:length(param.B_bins))=#

	for bVal in param.B_bins

		NN2          = convert(Int64,ceil(param.NN*bVal))
		samples      = collect(1:(param.nn-1))
		pSize        = collect(0:NN2)
		samplesFreqs = permutedims(pSize/NN2)
		neutralSfs   = @. 1/pSize
		replace!(neutralSfs, Inf => 0.0)

		f(x,y=param.nn) = Binomial(y,x)
		z    = f.(samplesFreqs)

		out  = pdf.(z,samples)
		out  = round.(out,digits=10)
		outS = SparseArrays.dropzeros(SparseArrays.sparse(out))
		convolutedBn[bVal] = outS
		# param.bn[bVal] = outS
		#=param.neut[bVal] = round.(param.B*(param.thetaMidNeutral)*0.25*(outS*neutralSfs),digits=10)=#

	end
end

"""

	phiReduction(param,gamma,ppos)

Reduction in fixation probabilty due to background selection and linkage. The formulas used have been subjected to several theoretical works ([Charlesworth B., 1994](https://doi.org/10.1017/S0016672300032365), [Hudson et al., 1995](https://www.genetics.org/content/141/4/1605), [Nordborg et al. 1995](https://doi.org/10.1017/S0016672300033619), [Barton NH., 1995](https://www.genetics.org/content/140/2/821)).

The fixation probabilty of selected alleles are reduce by a factor ``\\phi``:

```math
\\phi(t,s) = e^{[\\frac{-2\\mu}{t(1+\\frac{rL}{t}+\\frac{2s}{t})}]}
```

Multiplying across all deleterious linkes sites, we find:

```math
\\Phi = \\prod_{1}^{L} = \\phi(t,s) = e^{\\frac{-2t\\mu(\\psi[1,\\frac{r+2s+L}{r}] - \\psi[1,\\frac{r(L+1)+2s+t}{r}])}{r^2}}
```

```math
\\phi(t,s) = e^{\\frac{-2t\\mu(\\psi[1,\\frac{r+2s+L}{r}] - \\psi[1,\\frac{r(L+1)+2s+t}{r}])}{r^2}}
```

# Arguments
 - `gamma::Int64`: selection coefficient.

# Returns
 - `Float64`: expected rate of positive fixations under background selection.

"""
function phiReduction(param::parameters,gammaValue::Int64)
	S::Float64 = abs(param.gamNeg/(1.0*param.NN))
	r::Float64 = param.rho/(2.0*param.NN)
	μ::Float64 = param.thetaF/(2.0*param.NN)
	s::Float64 = gammaValue/(param.NN*1.0)

	Ψ0::Float64 = polygamma(1,(s+S)/r)
	Ψ1::Float64 = polygamma(1,(r+param.Lf*r+s+S)/r)
	CC::Float64 = 1.0

	out::Float64 = (ℯ^(-2.0*S*μ*(Ψ0-Ψ1)/(r^2)))
	return out
end

"""
	analyticalAlpha(param, convolutedSamples)

Analytical α(x) estimation. Solve α(x) generally. We used the expected rates of divergence and polymorphism to approach the asympotic value accouting for background selection, weakly and strong positive selection. α(x) can be estimated taking into account the role of positive selected alleles or not. In this way we explore the role of linkage to deleterious alleles in the coding region.

```math
\\mathbb{E}[\\alpha_{x}] =  1 - \\left(\\frac{\\mathbb{E}[D_{s}]}{\\mathbb{E}[D_{N}]}\\frac{\\mathbb{E}[P_{N}]}{\\mathbb{E}[P_{S}]}\\right)
```

# Arguments
 - `param::parameters`
 - `convolutedSamples::binomialDict`
# Returns
 - `Array{Float64,1}` α(x).
"""
function analyticalAlpha(;param::parameters,convolutedSamples::binomialDict)

	################################################################
		# Solve the model similarly to original python mktest  #
	################################################################
	B = param.B

	setThetaF!(param)
	thetaF = param.thetaF
	# Solve the probabilities of fixations without BGS
	## Set non-bgs
	param.B = 0.999
	## Solve the mutation rate
	setThetaF!(param)
	## Solve the probabilities
	setPpos!(param)
	# Return to the original values
	param.thetaF = thetaF
	param.B = B

	################################################################
	 # Accounting for positive alleles segregating due to linkage #
	################################################################

	# Fixation
	fN     = param.B*fixNeut(param)
	fNeg   = param.B*fixNegB(param,0.5*param.pposH+0.5*param.pposL)
	fPosL  = fixPosSim(param,param.gL,0.5*param.pposL)
	fPosH  = fixPosSim(param,param.gH,0.5*param.pposH)

	ds = fN
	dn = fNeg + fPosL + fPosH

	## Polymorphism
	neut = DiscSFSNeutDown(param,convolutedSamples.bn[param.B])

	selH::Array{Float64,1} = if isinf(exp(param.gH * 2))
		DiscSFSSelPosDownArb(param,param.gH,param.pposH,convolutedSamples.bn[param.B])
	else
		DiscSFSSelPosDown(param,param.gH,param.pposH,convolutedSamples.bn[param.B])
	end

	selL::Array{Float64,1} = DiscSFSSelPosDown(param,param.gL,param.pposL,convolutedSamples.bn[param.B])
	selN::Array{Float64,1} = DiscSFSSelNegDown(param,param.pposH+param.pposL,convolutedSamples.bn[param.B])

	splitColumns(matrix::Array{Float64,2}) = (view(matrix, :, i) for i in 1:size(matrix, 2));
	tmp = cumulativeSfs(hcat(neut,selH,selL,selN),false)

	neut, selH, selL, selN = splitColumns(tmp)
	sel = (selH+selL)+selN

	# ps = @. neut / (sel+neut)
	# pn = @. sel / (sel+neut)

	## Outputs
	α = @. 1 - ((ds/dn) * (sel/neut))

	#############################################################
	#		Accounting for neutral and deleterious alleles		#	#############################################################
	## Fixation
	fN_nopos     = fN*(param.thetaMidNeutral/2.)*param.TE*param.NN
	fNeg_nopos   = fNeg*(param.thetaMidNeutral/2.)*param.TE*param.NN
	fPosL_nopos  = fPosL*(param.thetaMidNeutral/2.)*param.TE*param.NN
	fPosH_nopos  = fPosH*(param.thetaMidNeutral/2.)*param.TE*param.NN

	ds_nopos = fN_nopos
	dn_nopos = fNeg_nopos + fPosL_nopos + fPosH_nopos

	## Polymorphism
	sel_nopos = selN
	ps_nopos = @. neut / (sel_nopos + neut)
	pn_nopos = @. sel_nopos / (sel_nopos + neut)

	α_nopos = 1 .- (ds_nopos/dn_nopos) .* (sel_nopos./neut)

	##########
	# Output #
	##########
	return (α,α_nopos)
end


text:
using Test, MLToolkit.DistributionsX
using MLToolkit: seed!
using Distributions: Normal, MvNormal
using Statistics: mean, var, std
using Flux: gpu
seed!(1)

@testset "Normal" begin
    n_randtests = 3
    n_samples = 10_000
    d, n = 5, 10
    rtol = 0.02

    @testset "BroadcastedNormalStd" begin
        for _ = 1:n_randtests
            # Construction
            m, s = randn(d, n) |> gpu, rand(d, n) |> gpu
            bd1 = BroadcastedNormalStd(m, s)
            # Statistics
            @test mean(bd1) == m
            @test std(bd1)  == s
            # Sample
            xs = [rand(bd1) for _ = 1:n_samples]
            @test mean(xs) ≈ mean(bd1) rtol=rtol
            @test std(xs)  ≈ std(bd1)  rtol=rtol
            # Sample (multi)
            X = rand(bd1, n_samples)
            @test mean(X; dims=3) ≈ m rtol=rtol
            @test  std(X; dims=3) ≈ s rtol=rtol
            # Density
            @test logpdf.(Normal.(m, s), X) ≈ logpdf(bd1, X)
            # Density (multi)
            X1 = X[:,:,1]
            m1, s1 = m[:,1], s[:,1]
            @test logpdf(MvNormal(m1, s1), X1) ≈ sum(logpdf(BroadcastedNormalStd(m1, s1), X1), true; dims=1)
            # KL divergence
            m, s = (m, s) .+ (randn(d, n) |> gpu, rand(d, n) |> gpu)
            bd2 = BroadcastedNormalStd(m, s)
            kl_mc = mean(logpdf.(bd1, xs) - logpdf.(bd2, xs))
            @test kldiv(bd1, bd2) ≈ kl_mc rtol=rtol
        end
    end

    @testset "BroadcastedNormalLogStd" begin
        for _ = 1:n_randtests
            # Construction
            m, s = randn(d, n) |> gpu, rand(d, n) |> gpu
            bd1 = BroadcastedNormalLogStd(m, log.(s))
            # Statistics
            @test mean(bd1) == m
            @test std(bd1)  ≈  s
            # Sample
            xs = [rand(bd1) for _ = 1:n_samples]
            @test mean(xs) ≈ mean(bd1) rtol=rtol
            @test std(xs)  ≈ std(bd1)  rtol=rtol
            # Sample (multi)
            X = rand(bd1, n_samples)
            @test mean(X; dims=3) ≈ m rtol=rtol
            @test  std(X; dims=3) ≈ s rtol=rtol
            # Density
            @test logpdf.(Normal.(m, s), X) ≈ logpdf(bd1, X)
            # Density (multi)
            X1 = X[:,:,1]
            m1, s1 = m[:,1], s[:,1]
            @test logpdf(MvNormal(m1, s1), X1) ≈ dropdims(sum(logpdf(BroadcastedNormalLogStd(m1, log.(s1)), X1); dims=1); dims=1)
            # KL divergence
            m, s = (m, s) .+ (randn(d, n) |> gpu, rand(d, n) |> gpu)
            bd2 = BroadcastedNormalLogStd(m, log.(s))
            kl_mc = mean(logpdf.(bd1, xs) - logpdf.(bd2, xs))
            @test kldiv(bd1, bd2) ≈ kl_mc rtol=rtol
        end
    end

    # using LinearAlgebra: I

    # @testset "MvNormal" begin
    #     for _ = 1:NUM_RANDTESTS
    #         μ1 = zeros(FT, d); Σ1 = Matrix{FT}(I, d, d)
    #         μ2 = μ1 .+ rand(FT); Σ2 = Σ1 .* abs(rand(FT))

    #         mvn1 = MvNormal(μ1, Σ1); mvn2 = MvNormal(μ2, Σ2)
    #         x = rand(mvn1, n)

    #         kl_12 = mean(logpdf(mvn1, x) - logpdf(mvn2, x))
    #         @test kldiv(mvn1, mvn2) ≈ kl_12 atol=(d * ATOL_RAND)
    #     end
    # end
end


text:
"""
    AbstractDH{T}

Abstract base class for data handler objects.
"""
abstract type AbstractDH{T} end
export AbstractDH

"""
    DataHandler{T} <: AbstractDH{T}

Type for handling datasets.  This is basically a wrapper for a dataframe with methods
for splitting it into training and test sets and creating input and output numerical
arrays.  It is intended that most reformatting of the dataframe is done before passing
it to an instance of this type.

The parameter T specifies the datatype of the input, output arrays.
"""
mutable struct DataHandler{T} <: AbstractDH{T}

    # dataframe where all data is kept
    df::DataTable

    colsInput::Vector{Symbol}
    colsOutput::Vector{Symbol}

    colsNormalize::Vector{Symbol}

    mu::Vector{T}
    norm::Vector{T}
    userange::Bool

    # training data, is a subset of df
    dfTrain::DataTable
    # testing data, is a subset of df
    dfTest::DataTable

    # these are the arrays of training, testing data
    X_train::Array{T}
    y_train::Array{T}
    X_test::Array{T}
    y_test::Array{T}

    # this is where predictions are kept
    yhat::Array{T}
    yhat_train::Array{T}

    """
        DataHandler{T}(df::DataTable; testfrac::AbstractFloat=0.0, shuffle::Bool=false,
                    input_cols::Array{Symbol}=Symbol[], output_cols::Array{Symbol}=Symbol[],
                    normalize_cols::Array{Symbol}=Symbol[], assign::Bool=false,
                    userange::Bool=false)

    Constructor for a `DataHandler` object.  The resulting object will produce machine
    learning datasets from the provided dataframe.  The columns for machine learning input
    and output (target) should be provided with keywords `input_cols` and `output_cols`
    respectively.  The constructor will randomly do a train-test split with test fraction
    `testfrac`.

    Note that currently `DataHandler` can't handle dataframes with null values, and
    automatically drops all rows containing null values.
    """
    function DataHandler{T}(df::DataTable; testfrac::AbstractFloat=0.0, shuffle::Bool=false,
                            input_cols::Vector{Symbol}=Symbol[],
                            output_cols::Vector{Symbol}=Symbol[],
                            normalize_cols::Vector{Symbol}=Symbol[],
                            assign::Bool=false,
                            userange::Bool=false) where T
        if sum(!completecases(df)) ≠ 0
            throw(ArgumentError("DataHandler only accepts complete dataframes."))
        end
        ndf = copy(df)
        # TODO convert to all non-nullable arrays!!!
        intersection = input_cols ∩ output_cols
        if length(intersection) ≠ 0
            warn("Estimator input and output both contain $intersection.")
        end
        o = new(ndf, input_cols, output_cols, normalize_cols)
        o.userange = userange
        split!(o, testfrac, shuffle=shuffle, assign=assign)
        computeNormalizeParameters!(o, dataset=:dfTrain)
        if canNormalize(o)
            normalizeTrain!(o)
            size(o.dfTest,1) > 0 && normalizeTest!(o)
        end
        o
    end
end
export DataHandler

# TODO all of this normalization stuff should be on the estimator side.
# get rid of it and make a good framework for doing this with estimators
"""
    computeNormalizeParameters!{T}(dh::AbstractDH{T}; dataset::Symbol=:dfTrain)

Gets the parameters for centering and rescaling from either the training dataset
(`dataset=:dfTrain`) or the test dataset (`dataset=:dfTest`).

Does this using the training dataframe by default, but can be set to use test.
Exits normally if this doesn't need to be done for any columns.

This should always be called before `normalize!`, that way you have control over what
dataset the parameters are computed from.
"""
function computeNormalizeParameters!{T}(dh::AbstractDH{T}; dataset::Symbol=:dfTrain)
    if length(dh.colsNormalize)==0 return end
    df = getfield(dh, dataset)
    mu = Array{T, 1}(length(dh.colsNormalize))
    norm = Array{T, 1}(length(dh.colsNormalize))
    for (i, col) in enumerate(dh.colsNormalize)
        dfcol = dropnull(df[col])
        mu[i] = mean(dfcol)
        if dh.userange
            norm[i] = maximum(dfcol) - minimum(dfcol)
        else
            norm[i] = std(dfcol)
        end
    end
    dh.mu = mu
    dh.norm = norm
    return
end
export computeNormalizeParameters!

"""
    canNormalize(dh::AbstractDH)

Determines whether the data in the datahandler can be normlized, i.e. because
the parameters have or haven't been computed yet.
"""
function canNormalize(dh::AbstractDH)
    can = isdefined(dh, :mu) && isdefined(dh, :norm)
    can = can && length(dh.mu) == length(dh.colsNormalize)
    can = can && length(dh.norm) == length(dh.colsNormalize)
    return can
end

"""
    normalize!{T}(dh::AbstractDH{T}; dataset::Symbol=:dfTrain)
    normalizeTrain!(dh::AbstractDH)
    normalizeTest!(dh::AbstractDH)

Centers and rescales the columns set by `normalize_cols` in the `DataHandler` constructor.
"""
function normalize!{T}(dh::AbstractDH{T}; dataset::Symbol=:dfTrain)
    if !canNormalize(dh)
        error("Trying to normalize before parameters have been computed.")
    end
    df = getfield(dh, dataset)
    # note that this does something because these are not deep copies
    for (i, col) in enumerate(dh.colsNormalize)
        dfcol = dropnull(df[col])
        dfcol = (dfcol - dh.mu[i])./dh.norm[i]
    end
    return
end
export normalize!
normalizeTrain!(dh::AbstractDH) = normalize!(dh, dataset=:dfTrain)
export normalizeTrain!
normalizeTest!(dh::AbstractDH) = normalize!(dh, dataset=:dfTest)
export normalizeTest!

"""
    unnormalize!{T}(dh::AbstractDH{T}, X::Matrix{T}, cols::Vector{Symbol})

Performs the inverse of the centering and rescaling operations on a matrix.
This can also be called on a single column with a `Symbol` as the last argument.
"""
function unnormalize!{T}(dh::AbstractDH{T}, X::Matrix{T}, cols::Vector{Symbol})
    err = "Array must have same number of dimensions as are to be inverted."
    @assert size(X)[2] == length(cols) err
    for (i, col) in enumerate(cols)
        idx = find(dh.colsNormalize .== col)[1]
        X[:, i] = dh.norm[idx]*X[:, i] + dh.mu[idx]
    end
end
export unnormalize!

function unnormalize!{T}(dh::AbstractDH{T}, X::Array{T, 2}, col::Symbol)
    unnormalize!(dh, X, [col])
end
export unnormalize!

"""
    shuffle!(dh::AbstractDH)

Shuffles the main dataframe of the DataHandler.
"""
shuffle!(dh::AbstractDH) = shuffle!(dh.df)
export shuffle!

function _get_assign_data{T}(dataname::Symbol, dh::DataHandler{T})
    data = getfield(dh, dataname)
    if isempty(data)
        throw(ErrorException("Attempting to assign data from empty data set."))
    end
    X = convert(Array{T}, data[:, dh.colsInput])
    y = convert(Array{T}, data[:, dh.colsOutput])
    X, y
end

"""
    assignTrain!(dh::AbstractDH)

Assigns the training data in the data handler so it can be retrieved in proper form.
"""
function assignTrain!{T}(dh::AbstractDH{T})
    X, y = _get_assign_data(:dfTrain, dh)
    dh.X_train = X
    dh.y_train = y
    X, y
end
export assignTrain!

"""
    assignTest!(dh::AbstractDH)

Assigns the test data in the data handler.
"""
function assignTest!{T}(dh::AbstractDH{T})
    X, y = _get_assign_data(:dfTest, dh)
    dh.X_test = X
    dh.y_test = y
    X, y
end
export assignTest!

"""
    assign!(dh::AbstractDH)

Assigns training and test data in the data handler.
"""
function assign!{T}(dh::AbstractDH{T})
    assignTrain!(dh)
    assignTest!(dh)
    nothing
end
export assign!

"""
    getTrainData(dh::AbstractDH; flatten::Bool=false)

Gets the training data input, output tuple `X, y`.

If `flatten`, attempts to flatten `y`.
"""
function getTrainData(dh::AbstractDH; flatten::Bool=false)
    X, y = dh.X_train, dh.y_train
    if flatten
        @assert size(y,2) == 1 "Attempted to flatten rank-2 array."
        y = squeeze(y, 2)
    end
    X, y
end
export getTrainData

"""
    getTestData(dh::AbstractDH; flatten::Bool=false)

Gets the test data input, output tuple `X, y`.

If `flatten`, attempts to flatten `y`.
"""
function getTestData(dh::AbstractDH; flatten::Bool=false)
    X, y = dh.X_test, dh.y_test
    if flatten
        @assert size(y,2) == 1 "Attempted to flatten rank-2 array."
        y = squeeze(y, 2)
    end
    X, y
end
export getTestData

"""
    getData(dh[; flatten=false])

Gets all data (i.e. not split into training and test).  Note that to use this one does
not have to call `assign!` since the data is taken directly from the main dataset.
"""
function getData(dh::AbstractDH; flatten::Bool=false)
    X, y = _get_assign_data(:df, dh)
    if flatten
        @assert size(y,2) == 1 "Attempted to flatten rank-2 array."
        y = squeeze(y,2)
    end
    X, y
end
export getData

"""
    split!(dh::AbstractDH, index::Integer; assign::Bool=true)

Creates a train, test split by index.  The index given is the last index of the training set.
If `assign`, this will assign the training and test data.
"""
function split!(dh::AbstractDH, index::Integer; assign::Bool=true)
    @assert index ≤ size(dh.df)[1] "Index value too large for dataframe."
    dh.dfTrain = dh.df[1:index, :]
    if index == size(dh.df)[1]
        dh.dfTest = DataTable()
    else
        dh.dfTest = dh.df[(index+1):end, :]
    end
    assign && assign!(dh)
    nothing
end
export split!

"""
    split!(dh::AbstractDH, testfrac::AbstractFloat; shuffle::Bool=false,
           assign::Bool=true)

Creates a train, test split by fraction.  The fraction given is the test fraction.
"""
function split!(dh::AbstractDH, testfrac::AbstractFloat; shuffle::Bool=false,
                assign::Bool=true)
    @assert testfrac ≤ 1.0 "Test fraction must be ∈ [0, 1]."
    if shuffle shuffle!(dh) end
    index = convert(Int64, round((1. - testfrac)*size(dh.df)[1]))
    split!(dh, index, assign=assign)
    nothing
end
export split!

"""
    split!(dh::AbstractDH, constraint::BitArray)

Splits the data into training and test sets using a BitArray that must correspond to elements
of dh.df.  The elements of the dataframe for which the BitArray holds 1 will be in the test
set, the remaining elements will be in the training set.
"""
function split!(dh::AbstractDH, constraint::BitArray)
    dh.dfTest = dh.df[constraint, :]
    dh.dfTrain = dh.df[!constraint, :]
end
export split!

"""
    @selectTrain!(dh, expr)

Set the training set to be the subset of the `DataHandler`'s dataframe for which `expr` is
true.  `expr` should be an expression which evaluates to a `Bool` with `Symbol`s corresponding
to column names for values.  See the documentation for `@constrain` and `@split!`
for examples.
"""
macro selectTrain!(dh, expr)
    esc(:($dh.dfTrain = @constrain($dh.df, $expr)))
end
export @selectTrain!

"""
    @selectTrain!(dh, expr)

Set the test set to be the subset of the `DataHandler`'s dataframe for which `expr` is
true.  `expr` should be an expression which evaluates to a `Bool` with `Symbol`s corresponding
to column names for values.  See the documentation for `@constrain` and `@split!`
for examples.
"""
macro selectTest!(dh, expr)
    esc(:($dh.dfTest = @constrain($dh.df, $expr)))
end
export @selectTest!

"""
    @split!(dh, expr)

Splits the `DataHandler`'s DataTable into training in test set, such that the test set is the
set of datapoints for which `expr` is true, and the training set is the set of datapoints for
which `expr` is false.  `expr` should be an expression that evaluates to `Bool` with symbols
in place of column names.  For example, if the columns of the dataframe are `[:x1, :x2, :x3]`
one can do `expr = (:x1 > 0.0) | (tanh(:x3) > e^6)`.

See documentation on the `@constrain` macro which `@split!` calls internally.
"""
macro split!(dh, expr)
    esc(quote
        $dh.dfTest  = @constrain($dh.df, $expr)
        $dh.dfTrain = @constrain($dh.df, !($expr))
    end)
end
export @split!

function getDefaultTestAnalysisColumnNames(dh::AbstractDH)
    [Symbol(string(col, "_hat")) for col ∈ dh.colsOutput]
end

"""
## `DataHandler`

    getTestAnalysisData(dh::AbstractDH, ŷ::Array; names::Vector{Symbol}=Symbol[],
                        squared_error::Bool=true)

Creates a dataframe from the test dataframe and a supplied prediction.

The array names supplies the names for the columns, otherwise will generate default names.

Also generates error columns which are the difference between predictions and test data.
If `squared_error`, will also create a column with squared error.

Note that this currently does nothing to handle transformations of the data.
"""
function getTestAnalysisData(dh::AbstractDH, ŷ::Array; names::Vector{Symbol}=Symbol[],
                             squared_error::Bool=true)
    # convert vectors to matrices
    if length(size(ŷ)) == 1
        ŷ = reshape(ŷ, (length(ŷ), 1))
    end
    @assert size(ŷ,2) == length(dh.colsOutput) ("Supplied array must have same number of
                                                  columns as the handler's output.")
    if length(names) == 0
        names = getDefaultTestAnalysisColumnNames(dh)
    end

    @assert length(dh.colsOutput) == length(names) ("Wrong number of provided names.")
    # if ŷ is short it is assumed to correspond to begining of dataframe, useful for
    # time series where only a partial sequence has been generated
    df = dh.dfTest[1:size(ŷ)[1], :]
    for (idx, name) ∈ enumerate(names)
        df[name] = ŷ[:, idx]
        orig_col = convert(Vector, df[dh.colsOutput[idx]])
        err = ŷ[:, idx] - orig_col
        df[Symbol(string(name, "_Error"))] = err
        if squared_error
            df[Symbol(string(name, "_Error²"))] = err.^2
        end
    end
    df
end
export getTestAnalysisData


text:
"""
ptciBeam(n=63, delta0=1.e-2, PvecKnowsdelta=true, pside = "right";
lsolver="gmres")

Solves the buckling beam problem with ptcsoli. You can play
- left/right preconditioning
- pseudo time step dependent preconditioning
- relationship of delta0 to n (hint, it's not mesh-independent)
"""
function ptciBeam(n = 63, delta0 = 1.e-2, PvecKnowsdelta = true,
pside = "right"; lsolver="gmres")
    lambda = 20.0
    maxit = 1000
    delta0 = 0.01
    PvecKnowsdelta ? Pvec = ptvbeamdelta : Pvec = ptvbeam
    bdata = beaminit(n, 0.0, lambda)
    x = bdata.x
    u0 = x .* (1.0 .- x) .* (2.0 .- x)
    u0 .*= exp.(-10.0 * u0)
    FS = copy(u0)
    FPJV = zeros(n, 20)
    pout = ptcsoli(
        FBeam!,
        u0,
        FS,
        FPJV;
        lsolver=lsolver,
        delta0 = delta0,
        pdata = bdata,
        lmaxit = 19,
        eta = 1.e-2,
        rtol = 1.e-10,
        maxit = maxit,
        Pvec = Pvec,
        PvecKnowsdelta = PvecKnowsdelta,
        pside = pside,
    )
    return pout
end

"""
ptvbeamdelta(v, x, bdata)

Precondition buckling beam problem with delta-aware preconditioner.
"""
function ptvbeamdelta(v, x, bdata)
    delta = bdata.deltaval[1]
    J = bdata.D2 + (1.0 / delta) * I
    ptv = J \ v
end

"""
ptvbeamp(v, x, bdata)

Precondition buckling beam problem with inverse of high-order term.
"""
function ptvbeam(v, x, bdata)
    J = bdata.D2
    ptv = J \ v
end


text:
"""
    Find the best question to ask by iterating over every feature / value
    and calculating the information gain.
"""
function find_best_split(rows)::Tuple{Float64, Union{Question, Nothing}}
    best_gain = 0.0  # keep track of the best information gain
    best_question::Union{Nothing, Question} = nothing  # keep track of the feature / value that produced it
    current_uncertainty = gini(rows)
    n_features = length(rows[1]) - 1  # number of columns

    for col in 1:n_features  # for each feature
        values::Set{Feature} = Set([row[col] for row in rows])  # unique values in the column
        for val in values  # for each value
            question = Question(col, val)

            # try splitting the dataset
            true_rows, false_rows = partition(rows, question)

            # Skip this split if it doesn't divide the dataset.
            if length(true_rows) == 0 || length(false_rows) == 0
                continue
            end
            # Calculate the information gain from this split
            gain = info_gain(true_rows, false_rows, current_uncertainty)

            # You actually can use ">" instead of ">=" here
            # but I wanted the tree to look a certain way for our
            # toy dataset.
            if gain >= best_gain
                best_gain, best_question = gain, question
            end # if
        end # for val in values
    end # for col in 1:n_features
    best_gain, best_question
end


text:
using ProgressMeter # If you didn't install ProgressMeter, do Pkg.add("ProgressMeter")
unirand(a,b) = a + (b-a) * rand(); # define fn returning uniform random number in (a,b)

# make directory to store the positions and the velocities data.
output_dir = "data"
if ! ispath(output_dir); mkdir(output_dir); end

# distane between i th particle and j th particle at periodic boundary condition
function dist(x1::Array{Float64,1} ,x2::Array{Float64,1} ,L::Float64);
    # x[1]: x position of a particle
    # x[2]: y position of a particle
    # L: system size
    dx = abs(x1[1] - x2[1]) % L;
    dx = minimum([dx; L-dx]);
    dy = abs(x1[2] - x2[2]) % L;
    dy = minimum([dy; L-dy]);
    return sqrt(dx^2 + dy^2)
end

# pairwise distance
function ppdist(x::Array{Float64,1} ,y::Array{Float64,1}, num_particle::Int64, L::Float64);
    D = zeros(num_particle, num_particle);
    for i in 1:num_particle-1
            for j in i+1:num_particle
                    D[i,j] = dist([ x[i]; y[i] ], [ x[j]; y[j] ], L);
            end
    end

    return D;
end

# calculate mean angle
function cal_mean_angle(D::Array{Float64,2}, theta::Array{Float64,1}, num_particle::Int64, r::Float64)
    meanangle = zeros(num_particle)
    for iter2 in 1:num_particle
        truth = D[:,iter2] .< r;
        near_particles_angles = theta[truth];
        meanangle[iter2] =
        atan2(
            mean( sin(near_particles_angles) ),
            mean( cos(near_particles_angles) )
        );
    end
    return meanangle;
end

function vicsek()
    # L:             system size
    # num_particle: the number of particles
    # vc:             constant speed
    # dt:             time step
    # steps:         the number of run steps
    # eta:             Order
    # r:             interaction radius
    const L = 5.0; const num_particle = 300; const vc = 1.0;
    const dt = 0.1; const steps = 500; const eta = 0.1; const r = 0.5;

    # initial position of the particles
    # x_i and y_i corespond to the i th positon of particle.
    x = rand(num_particle) * L; y = rand(num_particle) * L;

    # initial angle (-pi, pi) of the particles;
    theta = [unirand(-pi, pi) for i in 1:num_particle];
    #theta = unirand(-pi, pi, num_particle);

    # initial velocity
    vx = vc .* cos(theta); vy = vc .* sin(theta);

    # start fish swimming!
    progress = Progress(steps);
    for iter in 0:steps-1;
        # output configuration
        filename = @sprintf("./data/%05d.csv",iter);
        writecsv(filename, [x y vx vy])

        # pairwise distance
        # D[i,j] is a distance between i th and j th particles.
        D = ppdist(x, y, num_particle, L);
        D += D';

        # calculate mean angle
        meanangle = zeros(num_particle);
        meanangle = cal_mean_angle(D, theta, num_particle, r);

        # update angle with noise
        theta = meanangle .+ 2*pi*eta .* (rand(num_particle) .- 0.5);

        # update positions and velocities;
        vx = vc .* cos(theta); vy = vc .* sin(theta);
        x += vx * dt; y += vy * dt;

        # periodic boundary condition
        #x = (x .+ L) % L; y = (y .+ L) % L;
        x = mod(x,L); y = mod(y,L);

        next!(progress);
    end

end;
@time vicsek();

max_stars_repo_name: goropikari/vicsek-model

max_issues_repo_name: goropikari/vicsek-model

max_forks_repo_name: goropikari/vicsek-model


text:
using DifferentialEquations
using Distributions
using Plots
using DelimitedFiles

#-----------------------------------
# Number of cells
#-----------------------------------

nCell = 3

#-----------------------------------
# Signaling topology
#-----------------------------------

topology = "chain" # "loop" or "chain"

#-------------------------------------------------
# Wait time distribution and pulse time
#-------------------------------------------------

mean = 50.0

pulse = 5.0

#-----------------------------------
# Number of iterations
#-----------------------------------

nits = 10

#-----------------------------------
# lambda
#-----------------------------------

lambda = 1.0

#-----------------------------------
# Range of A0 values
#-----------------------------------

A0_values = 0.7:0.005:0.72

#-----------------------------------
# Noise and the variance
#-----------------------------------

noiseE = "off" # "on" or "off", extracellular noise
noiseP = "off" # "on" or "off", intrinsic noise

etaE = Normal(0.0, 0.01)
etaP = Normal(0.0, 0.1)

#-----------------------------------
# Internal ODE system
#-----------------------------------

function GataPu1(du,u,p,t)
    a = p[1:5]
    b = p[6:8]
    g = p[9:17]
    L = p[18:20]
    du[1] = -b[1]*u[1] + (a[1]*L[1] + a[2]*u[1])/(1 + g[1]*L[1] + g[2]*u[1] + g[3]*u[1]*u[2])
    du[2] = -b[2]*u[2] + (a[3]*L[2] + a[4]*u[2])/(1 + g[4]*L[2] + g[5]*u[2] + g[6]*u[1]*u[2] + g[7]*u[1]*u[3])
    du[3] = -b[3]*u[3] + (a[5]*u[1])/(1 + g[8]*u[1] + g[9]*L[3])
end

#-----------------------------------
# Parameter vector p
#-----------------------------------

p = zeros(0)
a = [1.0,0.25,1.0,0.25,0.01]
b = [0.01,0.01,0.01]
g = [1.0,0.25,1.0,1.0,0.25,1.0,0.13,0.01,10]
L = [1.0,0.5,0] #L = [A,B,C] = [A,0.5,0]
append!(p,a)
append!(p,b)
append!(p,g)
append!(p,L)

#-----------------------------------
# Wait time distribution
#-----------------------------------

function rand_exponential(mean)
    if mean <= 0.0
        error("mean must be positive")
    end
    -mean*log(rand())
end

function wait_time(mean)
    round(rand_exponential(mean))
end

#-----------------------------------
#Simulation Code
#-----------------------------------

GHighCounts = zeros(length(A0_values), nCell+1)

for l = 1:length(A0_values)
    A0 = A0_values[l]
    for j = 1:nits

        CellValues = [[20.0, 1.0, 20.0, A0] for i=1:nCell, n=1:1]

        if noiseE == "on"
            CellValues = [[20.0, 1.0, 20.0, A0 + rand(etaE)] for i=1:nCell, n=1:1]
        end

        CellIts = zeros(nCell)
        CellWaitTimes = [wait_time(mean) for i=1:nCell, n=1:1]
        CellState = ["wait" for i=1:nCell, n=1:1]

        while min(CellIts...) < 1000 + 40*nCell
            t = min(CellWaitTimes...)

            for i = 1:nCell
                u0 = [CellValues[i][1], CellValues[i][2], CellValues[i][3]]
                p[18] = CellValues[i][4]
                tspan = (0.0, t)
                prob = ODEProblem(GataPu1,u0,tspan,p)
                sol = solve(prob)
                CellValues[i][1] = sol[length(sol)][1]
                CellValues[i][2] = sol[length(sol)][2]
                CellValues[i][3] = sol[length(sol)][3]
            end

            CellWaitTimes = CellWaitTimes .- t

            if topology == "loop"
                if CellWaitTimes[1] == 0
                    if CellState[1] == "wait"
                        CellWaitTimes[1] = pulse
                        CellState[1] = "pulse"
                        CellIts[1] = CellIts[1] + 1
                        CellValues[1][4] = A0 * (CellValues[nCell][1]/(lambda + CellValues[nCell][2]))
                        if noiseP == "on"
                            CellValues[1][4] = A0 * (CellValues[nCell][1]/(lambda + CellValues[nCell][2])) + rand(etaP)
                        end
                    else
                        CellWaitTimes[1] = wait_time(mean)
                        CellState[1] = "wait"
                        CellValues[1][4] = A0
                        if noiseE == "on"
                            CellValues[1][4] = A0 + rand(etaE)
                        end
                    end
                end
            else
                CellWaitTimes[1] = 100.0
                CellIts[1] = 2000 + 40*nCell
            end

            for i = 2:nCell
                if CellWaitTimes[i] == 0
                    if CellState[i] == "wait"
                        CellWaitTimes[i] = pulse
                        CellState[i] = "pulse"
                        CellIts[i] = CellIts[i] + 1
                        CellValues[i][4] = A0 * (CellValues[i-1][1]/(lambda + CellValues[i-1][2]))
                        if noiseP == "on"
                            CellValues[i][4] = A0 * (CellValues[i-1][1]/(lambda + CellValues[i-1][2])) + rand(etaP)
                        end
                    else
                        CellWaitTimes[i] = wait_time(mean)
                        CellState[i] = "wait"
                        CellValues[i][4] = A0
                        if noiseE == "on"
                            CellValues[i][4] = A0 + rand(etaE)
                        end
                    end
                end
            end

        end

        for i = 1:nCell
            if CellValues[i][1] > 30
                GHighCounts[l, i+1] = GHighCounts[l, i+1] + 1
                GHighCounts[l, 1] = A0
            end
        end
    end
    for i = 1:nCell
        GHighCounts[l, i+1] = GHighCounts[l, i+1]/nits
    end
end

writedlm("../results/"*topology*"_$nCell"*"_cells.txt", GHighCounts, ", ")

max_stars_repo_name: alavendelm/Cell-Cell-Communication

max_issues_repo_name: alavendelm/Cell-Cell-Communication

max_forks_repo_name: alavendelm/Cell-Cell-Communication


text:
@testset "activations" begin
    x = randn(4, 5);

    @test gradcheck(x -> sum(logistic.(x)), x)
    @test gradcheck(x -> sum(softplus.(x)), x)
    @test gradcheck(x -> sum(softsign.(x)), x)
    @test gradcheck(x -> sum(logsigmoid.(x)), x)
    @test gradcheck(x -> sum(relu.(x)), x)
    @test gradcheck(x -> sum(leakyrelu.(x, 0.01)), x)
    @test gradcheck(x -> sum(elu.(x, 1)), x)

    @test gradcheck(x -> sum(softmax(x)), x)
    @test gradcheck(x -> sum(logsoftmax(x)), x)
end


text:
workspace();
using BoxChopper;

box = [0 1; 0 1; 0 1]
r₀ = [0.5, 0.5, 0.5]
nout = [-2,-1,0]

nr₀ = nout⋅r₀
vcbits = BoxChopper.calc_vcbits(box, nout, nr₀)
nvc = count_ones(vcbits)  # number of vertices contained

assert(BoxChopper.isquadsect(vcbits))
@code_warntype BoxChopper.rvol_quadsect(box, nout, nr₀, vcbits)


text:
function init_constraint!(
    com::CS.CoM,
    constraint::OrConstraint,
    fct,
    set::OrSet;
)
    set_impl_functions!(com,  constraint.lhs)
    set_impl_functions!(com,  constraint.rhs)
    lhs_feasible = true
    if constraint.lhs.impl.init
        lhs_feasible = init_constraint!(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    end
    rhs_feasible = true
    if constraint.rhs.impl.init
        rhs_feasible = init_constraint!(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    end
    return lhs_feasible || rhs_feasible
end

"""
    is_constraint_solved(
        constraint::OrConstraint,
        fct,
        set::OrSet,
        values::Vector{Int},
    )

Checks if at least one of the inner constraints is solved
"""
function is_constraint_solved(
    constraint::OrConstraint,
    fct,
    set::OrSet,
    values::Vector{Int},
)
    lhs_num_vars = get_num_vars(constraint.lhs.fct)
    lhs_solved = is_constraint_solved(constraint.lhs, constraint.lhs.fct, constraint.lhs.set, values[1:lhs_num_vars])
    lhs_solved && return true
    rhs_num_vars = get_num_vars(constraint.rhs.fct)
    rhs_solved = is_constraint_solved(constraint.rhs, constraint.rhs.fct, constraint.rhs.set, values[end-rhs_num_vars+1:end])
    rhs_solved && return true
end

"""
    function is_constraint_violated(
        com::CoM,
        constraint::OrConstraint,
        fct,
        set::OrSet,
    )

Check if both inner constraints are violated already
"""
function is_constraint_violated(
    com::CoM,
    constraint::OrConstraint,
    fct,
    set::OrSet,
)
    lhs_violated = is_constraint_violated(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    rhs_violated = is_constraint_violated(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    return lhs_violated && rhs_violated
end

"""
    still_feasible(com::CoM, constraint::OrConstraint, fct, set::OrSet, vidx::Int, value::Int)

Return whether the constraint can be still fulfilled when setting a variable with index `vidx` to `value`.
**Attention:** This assumes that it isn't violated before.
"""
function still_feasible(
    com::CoM,
    constraint::OrConstraint,
    fct,
    set::OrSet,
    vidx::Int,
    value::Int,
)
    lhs_feasible = !is_constraint_violated(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    if lhs_feasible
        lhs_indices = constraint.lhs.indices
        for i in 1:length(lhs_indices)
            if lhs_indices[i] == vidx
                lhs_feasible = still_feasible(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set, vidx, value)
                lhs_feasible && return true
                break
            end
        end
    end
    rhs_feasible = !is_constraint_violated(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    if rhs_feasible
        rhs_indices = constraint.rhs.indices
        for i in 1:length(rhs_indices)
            if rhs_indices[i] == vidx
                rhs_feasible = still_feasible(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set, vidx, value)
                rhs_feasible && return true
                break
            end
        end
    end
    return rhs_feasible || lhs_feasible
end

"""
    prune_constraint!(com::CS.CoM, constraint::OrConstraint, fct, set::OrSet; logs = true)

Reduce the number of possibilities given the `OrConstraint` by pruning both parts
Return whether still feasible
"""
function prune_constraint!(
    com::CS.CoM,
    constraint::OrConstraint,
    fct,
    set::OrSet;
    logs = true,
)
    lhs_violated = is_constraint_violated(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    rhs_violated = is_constraint_violated(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    if lhs_violated && rhs_violated
        return false
    end
    if lhs_violated
        activate_rhs!(com, constraint)
        return prune_constraint!(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set; logs=logs)
    end
    if rhs_violated
        activate_lhs!(com, constraint)
        return prune_constraint!(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set; logs=logs)
    end
    return true
end


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using Weave

function compile_report(fname, out_type = "pdf")
    if out_type == "pdf"
        Weave.weave(fname, doctype="md2pdf", template="bee4750.tpl", fig_path="figures")
    elseif out_type == "html"
        Weave.weave(fname, doctype="md2html", fig_path="figures")
    else
        throw(ArgumentError("out_type must be pdf or html"))
    end
end

if !isinteractive()
    compile_report(ARGS...)
end


text:
# -------------------------------------------
# set_initQbase
# 基本量の初期値を設定
# リスタートするか初期値を代入する
# -------------------------------------------

function set_initQbase(xmax, ymax, zmax, restart_file, init_rho,
                        init_u, init_v, init_w, init_p,
                        specific_heat_ratio, out_file_front,
                        out_ext, out_dir, restartnum, Rd)
    Qbase=[]
    cellxmax = xmax - 1
    cellymax = ymax - 1
    cellzmax = zmax - 1

    restart_check = 0
    try Qbase = setup_restart_value(cellxmax, cellymax, cellzmax, out_dir, restart_file)
        println("Restart "*restart_file)
        restart_check = 2
    catch
        restart_check = 1
    end

    if restart_check == 1
        Qbase = setup_init_value(cellxmax, cellymax, cellzmax, init_rho, init_u, init_v, init_w, init_p)
        println("Start Initial condition")
        restart_num = 0
        yplus = zeros(cellxmax, cellymax, cellzmax)
        output_result(0, Qbase, cellxmax, cellymax, cellzmax, out_file_front, out_ext, out_dir, Rd, yplus)
    end

    return Qbase, cellxmax, cellymax, cellzmax, restart_num
end

function setup_init_value(cellxmax, cellymax, cellzmax, init_rho, init_u, init_v, init_w, init_p)
    Qbase = zeros(cellxmax, cellymax, cellzmax, 5)
    for i in 1:cellxmax
        for j in 1:cellymax
            for k in 1:cellzmax
                Qbase[i,j,k,1] = init_rho
                Qbase[i,j,k,2] = init_u
                Qbase[i,j,k,3] = init_v
                Qbase[i,j,k,4] = init_w
                Qbase[i,j,k,5] = init_p
            end
        end
    end
    return Qbase
end

function setup_restart_value(cellxmax, cellymax, cellzmax, out_dir, restart_file)
    Qbase = zeros(cellxmax, cellymax, cellzmax, 5)

    skipnum = 1
    fff = []
    open("result/"*restart_file, "r") do f
        fff = read(f,String)
    end
    fff = split(fff,"\n",keepempty=false)

    for i in 1+skipnum:length(fff)
        fff[i] = replace(fff[i]," \r" => "")
    end

    ite = 1
    for i in 2:cellxmax-1
        for j in 2:cellymax-1
            for k in 2:cellzmax-1
                temp = split(fff[ite+skipnum]," ")

                for l in 1:5
                    Qbase[i,j,k,l] = parse(Float64,temp[l])
                end
                ite = ite+1
            end
        end
    end
    return Qbase
end

function quicksort(list, first, last, num_cell_data, standard_num)
    # list = zezros(m,n)
    x = list[Int(floor((first+last)/2)),standard_num]
    i = first
    j = last

    while true
        while list[i,standard_num] < x
            i = i+1
        end
        while x < list[j,standard_num]
            j = j-1
        end

        if (i >= j)
            break
        else
            for k in 1:num_cell_data
                t = list[i,k]
                list[i,k] = list[j,k]
                list[j,k] = t
            end
            i = i+1
            j = j-1
        end
    end

    if (first < i - 1)
        list = quicksort(list, first, i - 1,num_cell_data,standard_num)
    end
    if (j + 1 < last)
        list = quicksort(list, j + 1, last,num_cell_data,standard_num)
    end
    return list
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


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

struct TableSelection{T,C}
  table::T
  cols::C
  names::Vector{Symbol}
  function TableSelection(table::T, names::Vector{Symbol}) where {T}
    cols = Tables.columns(table)
    @assert names ⊆ Tables.columnnames(cols)
    new{T,typeof(cols)}(table, cols, names)
  end
end

function Base.:(==)(a::TableSelection, b::TableSelection)
  a.names != b.names && return false
  all(Tables.getcolumn(a, nm) == Tables.getcolumn(b, nm) for nm in a.names)
end

function Base.show(io::IO, t::TableSelection)
  println(io, "TableSelection")
  pretty_table(io, t, vcrop_mode=:middle)
end

# Tables.jl interface
Tables.istable(::Type{<:TableSelection}) = true
Tables.columnaccess(::Type{<:TableSelection}) = true
Tables.columns(t::TableSelection) = t
Tables.columnnames(t::TableSelection) = t.names
Tables.getcolumn(t::TableSelection, i::Int) =
  Tables.getcolumn(t.cols, t.names[i])
function Tables.getcolumn(t::TableSelection, nm::Symbol)
  nm ∉ t.names && error("Table has no column $nm.")
  Tables.getcolumn(t.cols, nm)
end

function Tables.schema(t::TableSelection)
  schema = Tables.schema(t.table)
  names = schema.names
  types = schema.types
  inds = indexin(t.names, collect(names))
  Tables.Schema(t.names, types[inds])
end

Tables.materializer(t::TableSelection) =
  Tables.materializer(t.table)

"""
    Select(col₁, col₂, ..., colₙ)
    Select([col₁, col₂, ..., colₙ])
    Select((col₁, col₂, ..., colₙ))

The transform that selects columns `col₁`, `col₂`, ..., `colₙ`.

    Select(regex)

Selects the columns that match with `regex`.

# Examples

```julia
Select(1, 3, 5)
Select([:a, :c, :e])
Select(("a", "c", "e"))
Select(r"[ace]")
```
"""
struct Select{S<:ColSpec} <: Stateless
  colspec::S
end

# argument errors
Select(::Tuple{}) = throw(ArgumentError("Cannot create a Select object with empty tuple."))
Select() = throw(ArgumentError("Cannot create a Select object without arguments."))

Select(cols::T...) where {T<:ColSelector} =
  Select(cols)

isrevertible(::Type{<:Select}) = true

function apply(transform::Select, table)
  # original columns
  cols = Tables.columns(table)

  # retrieve relevant column names
  allcols = collect(Tables.columnnames(cols))
  select  = choose(transform.colspec, allcols)
  reject  = setdiff(allcols, select)

  # keep track of indices to revert later
  sinds = indexin(select, allcols)
  rinds = indexin(reject, allcols)

  # sort indices to facilitate reinsertion
  sperm  = sortperm(sinds)
  sorted = sortperm(rinds)
  reject = reject[sorted]
  rinds  = rinds[sorted]

  # rejected columns
  rcols = [Tables.getcolumn(cols, name) for name in reject]

  TableSelection(table, select), (reject, rcols, sperm, rinds)
end

function revert(::Select, newtable, cache)
  # selected columns
  cols   = Tables.columns(newtable)
  select = Tables.columnnames(cols)
  # https://github.com/JuliaML/TableTransforms.jl/issues/76
  scols  = Any[Tables.getcolumn(cols, name) for name in select]

  # rejected columns
  reject, rcols, sperm, rinds = cache

  # restore rejected columns
  anames = collect(select[sperm])
  acols  = collect(scols[sperm])
  for (i, rind) in enumerate(rinds)
    insert!(anames, rind, reject[i])
    insert!(acols, rind, rcols[i])
  end
  𝒯 = (; zip(anames, acols)...)
  𝒯 |> Tables.materializer(newtable)
end

# reverting a single TableSelection is trivial
revert(::Select, newtable::TableSelection, cache) = newtable.table

"""
    Reject(col₁, col₂, ..., colₙ)
    Reject([col₁, col₂, ..., colₙ])
    Reject((col₁, col₂, ..., colₙ))

The transform that discards columns `col₁`, `col₂`, ..., `colₙ`.

    Reject(regex)

Discards the columns that match with `regex`.

# Examples

```julia
Reject(:b, :d, :f)
Reject(["b", "d", "f"])
Reject((2, 4, 6))
Reject(r"[bdf]")
```
"""
struct Reject{S<:ColSpec} <: Stateless
  colspec::S
end

# argumet erros
Reject(::Tuple{}) = throw(ArgumentError("Cannot create a Reject object with empty tuple."))
Reject(::Colon) = throw(ArgumentError("Cannot reject all columns."))
Reject() = throw(ArgumentError("Cannot create a Reject object without arguments."))

Reject(cols::T...) where {T<:ColSelector} =
  Reject(cols)

isrevertible(::Type{<:Reject}) = true

function apply(transform::Reject, table)
  cols = Tables.columns(table)
  allcols = Tables.columnnames(cols)
  reject  = choose(transform.colspec, allcols)
  select  = setdiff(allcols, reject)
  strans  = Select(select)
  newtable, scache = apply(strans, table)
  newtable, (strans, scache)
end

function revert(::Reject, newtable, cache)
  strans, scache = cache
  revert(strans, newtable, scache)
end


text:
@doc doc"""
    SpecialEuclidean(n)

Special Euclidean group $\mathrm{SE}(n)$, the group of rigid motions.

``\mathrm{SE}(n)`` is the semidirect product of the [`TranslationGroup`](@ref) on $ℝ^n$ and
[`SpecialOrthogonal(n)`](@ref)

````math
\mathrm{SE}(n) ≐ \mathrm{T}(n) ⋊_θ \mathrm{SO}(n),
````

where $θ$ is the canonical action of $\mathrm{SO}(n)$ on $\mathrm{T}(n)$ by vector rotation.

This constructor is equivalent to calling

```julia
Tn = TranslationGroup(n)
SOn = SpecialOrthogonal(n)
SemidirectProductGroup(Tn, SOn, RotationAction(Tn, SOn))
```

Points on $\mathrm{SE}(n)$ may be represented as points on the underlying product manifold
$\mathrm{T}(n) \times \mathrm{SO}(n)$ or as affine matrices with size `(n + 1, n + 1)`.
"""
const SpecialEuclidean{N} = SemidirectProductGroup{
    TranslationGroup{Tuple{N},ℝ},
    SpecialOrthogonal{N},
    RotationAction{TranslationGroup{Tuple{N},ℝ},SpecialOrthogonal{N},LeftAction},
}

function SpecialEuclidean(n)
    Tn = TranslationGroup(n)
    SOn = SpecialOrthogonal(n)
    A = RotationAction(Tn, SOn)
    return SemidirectProductGroup(Tn, SOn, A)
end

show(io::IO, ::SpecialEuclidean{n}) where {n} = print(io, "SpecialEuclidean($(n))")

Base.@propagate_inbounds function submanifold_component(
    ::SpecialEuclidean{n},
    x::AbstractMatrix,
    ::Val{1},
) where {n}
    return view(x, 1:n, n + 1)
end
Base.@propagate_inbounds function submanifold_component(
    ::SpecialEuclidean{n},
    x::AbstractMatrix,
    ::Val{2},
) where {n}
    return view(x, 1:n, 1:n)
end

function submanifold_components(G::SpecialEuclidean{n}, x::AbstractMatrix) where {n}
    @assert size(x) == (n + 1, n + 1)
    @inbounds t = submanifold_component(G, x, Val(1))
    @inbounds R = submanifold_component(G, x, Val(2))
    return (t, R)
end

Base.@propagate_inbounds function _padpoint!(
    ::SpecialEuclidean{n},
    y::AbstractMatrix,
) where {n}
    for i ∈ 1:n
        y[n+1, i] = 0
    end
    y[n+1, n+1] = 1
    return y
end

Base.@propagate_inbounds function _padvector!(
    ::SpecialEuclidean{n},
    v::AbstractMatrix,
) where {n}
    for i ∈ 1:n+1
        v[n+1, i] = 0
    end
    return v
end

compose(::SpecialEuclidean, x::AbstractMatrix, y::AbstractMatrix) = x * y

function compose!(
    ::SpecialEuclidean,
    z::AbstractMatrix,
    x::AbstractMatrix,
    y::AbstractMatrix,
)
    return mul!(z, x, y)
end


text:
"""
Reconnection rate postprocessing related to ParaView.

Based on the tests, using concave (1) or convex (2) hull will effect the
potential drop calculation. Method 1 will at most give a value 10% smaller than
method 2.

In the near future, I want to call ParaView directly in Julia.
Hongyang Zhou, hyzhou@umich.edu 11/11/2019
modified on 03/16/2020
"""
module ReconnectionRate

using DelimitedFiles, NumericalIntegration, LinearAlgebra, Random, PyPlot
using ConcaveHull, LazySets
#using MATLAB

export Rg, e
export findBoundary, saveBoundary, createSeeds, integrate_along_boundary,
   integrate_along_boundary_hall

const Rg = 2634000        # [m], radius of Ganymede
const e  = 1.60217662e-19 # [C], electron charge

"""
	findBoundary(filename; DoPlot=false, method=1, nSmooth=3)

Identify boundary points from Paraview polyline tracing output.
`method` 1 is the convex hull of the boundary points, 2 is the concave hull
computed with KNN, and 3 is the Matlab `boundary` function call.  `nSmooth` is
the number of point considered in KNN.
To use 2 or 3,
you need to import the corresponding packages.
"""
function findBoundary(filename; DoPlot=false, method=1, nSmooth=3)

   fieldline = readdlm(filename, ',', Float32, header=true)
   data = fieldline[1]

   # count the total number of field lines
   nlinetotal = 0
   start = Int32[]
   for i = 1:size(data,1)
      if data[i,4] == 0.0 # IntegrationTime starts at 0.0
         nlinetotal += 1
         push!(start, i)
      end
   end

   line = Array{Float32, 2}(undef, 6, nlinetotal)

   for i = 1:nlinetotal-1
      iStart, iEnd = start[i], start[i+1]-1
      line[1,i] = data[iStart,5]
      line[2,i] = data[iStart,6]
      line[3,i] = data[iStart,7]
      line[4,i] = data[iEnd,5]
      line[5,i] = data[iEnd,6]
      line[6,i] = data[iEnd,7]
   end

   line[1,end] = data[start[end],5]
   line[2,end] = data[start[end],6]
   line[3,end] = data[start[end],7]
   line[4,end] = data[end,5]
   line[5,end] = data[end,6]
   line[6,end] = data[end,7]

   # status
   # (0: open, 1: closed along B, 2: closed along -B 3: fully closed)
   # (-1: cells inside body, -2: loop ray within block, -3: strange)
   status = zeros(Int8, nlinetotal)

   # Cluster seed points based on their end points
   for i = 1:nlinetotal
      if line[4,i]^2 + line[5,i]^2 + line[6,i]^2 < 1.1^2
         status[i] = 1
      end
   end

   pts_closed = line[1:2, status .== 1]
   pts_open = line[1:2, status .== 0]

   if method == 1 # LazySets, convex hull
      x = pts_closed[1,:]
      y = pts_closed[2,:]
      v = [[x[i], y[i]] for i in 1:length(x)]
      hull = convex_hull(v) # LazySets
      xy = hcat(hull...)
      xB = [xy[1,:]; xy[1,1]]
      yB = [xy[2,:]; xy[2,1]]
   elseif method == 2 # ConcaveHull
      ϵ = Float32(1e-6)
      Random.seed!(1234)
      x = [i + ϵ * rand(Float32) for i in pts_closed[1,:]]
      y = [i + ϵ * rand(Float32) for i in pts_closed[2,:]]
      v = [[x[i], y[i]] for i in 1:length(x)]
      hull = concave_hull(v, nSmooth) # ConcaveHull
      xB = [p[1] for p in hull.vertices]
      yB = [p[2] for p in hull.vertices]
   elseif method ==3 # Matlab
      x = pts_closed[1,:]
      y = pts_closed[2,:]
      k, area = mxcall(:boundary, 2, x, y)
      k = Int.(k)
      xB, yB = x[k], y[k]
   end
   zB = fill(2.0f0, length(xB))

   if DoPlot
      fig, ax = subplots(1,2)
      ax[1].scatter(pts_closed[1,:], pts_closed[2,:])
      ax[1].scatter(pts_open[1,:], pts_open[2,:])
      ax[2].scatter(xB, yB)
   end

   return xB, yB, zB
end

function saveBoundary(x,y,z, filename="boundary.txt")
   open(filename, "w") do io
      write(io, "\"X\",\"Y\",\"Z\"\n")
      writedlm(io, [x y z], ',')
   end
end

"Create seeds for field line tracing."
function createSeeds(filename="seeds.txt", dir="./")

   xMin, xMax = -2.0, 4.3 # xrange
   yMin, yMax = -3.0, 3.1 # yrange
   z = 2.0 # z cut plane

   dx, dy = 1/32, 1/32 # resolutions

   xRange = xMin:dx:xMax
   yRange = yMin:dy:yMax

   X = [x for x in xRange, y in yRange if 3.5^2 > (x-1.2)^2 + (y-0.05)^2 > 2.7^2]
   Y = [y for x in xRange, y in yRange if 3.5^2 > (x-1.2)^2 + (y-0.05)^2 > 2.7^2]
   Z = fill(z, size(X))

   open(joinpath(dir,filename), "w") do io
      write(io, "\"X\",\"Y\",\"Z\"\n")
      writedlm(io, [X Y Z], ',')
   end
end

"Integrate U cross B along the 2D boundary line."
function integrate_along_boundary(filename, DoPlot=false)
   bc = readdlm(filename, ',', header=true)
   data = bc[1]

   x = @view data[:,7]
   y = @view data[:,8]

   # Get E field: E = -uxB
   Ex = Vector{Float64}(undef,size(data,1))
   Ey = Vector{Float64}(undef,size(data,1))
   Ez = Vector{Float64}(undef,size(data,1))

   for i = 1:size(data,1)
      Ex[i], Ey[i], Ez[i] = -cross(data[i,1:3], data[i,4:6])
   end

   ϕx = cumul_integrate(x, Ex)
   ϕy = cumul_integrate(y, Ey)
   Coef = 1e3 * 1e-9 * Rg * 1e-3 # [km/s]*[nT]*[Rg] --> [kV]
   ϕ  = (ϕx .+ ϕy) .* Coef # [kV]

   if DoPlot
      figure()
      plot(ϕ)
   end

   maximum(ϕ) - minimum(ϕ)
end

"Integrate Ue cross B along the magnetopause boundary."
function integrate_along_boundary_hall(filename, DoPlot=false)
   bc = readdlm(filename, ',', header=true)
   data = bc[1]

   ρ = @view data[:,1]
   ux= @view data[:,2]
   uy= @view data[:,3]
   uz= @view data[:,4]
   bx= @view data[:,5]
   by= @view data[:,6]
   bz= @view data[:,7]
   jx= @view data[:,8]
   jy= @view data[:,9]
   jz= @view data[:,10]
   x = @view data[:,11]
   y = @view data[:,12]

   # Get E field: E = -uxB
   Ex = Vector{Float64}(undef,size(data,1))
   Ey = Vector{Float64}(undef,size(data,1))
   Ez = Vector{Float64}(undef,size(data,1))

   # Calculate electron bulk velocity in Hall MHD
   # Transform the hall velocity into [km/s]
   Coef = 1e-6 /e * 1e-3 * 1e-6
   uex = @. ux - jx/ρ * Coef
   uey = @. uy - jy/ρ * Coef
   uez = @. uz - jz/ρ * Coef

   for i = 1:size(data,1)
      Ex[i], Ey[i], Ez[i] =
         -cross([uex[i], uey[i], uez[i]], [bx[i], by[i], bz[i]])
   end

   ϕx = cumul_integrate(x, Ex)
   ϕy = cumul_integrate(y, Ey)
   Coef = 1e3 * 1e-9 * Rg * 1e-3 # [km/s]*[nT]*[Rg] --> [kV]
   ϕ  = (ϕx .+ ϕy) .* Coef # [kV]

   if DoPlot
      figure()
      plot(ϕ)
   end

   #maximum(ϕ) - minimum(ϕ) # This would return the CPCP on the night side
   maximum(ϕ) - minimum(ϕ) + ϕ[end] # Return the CPCP on the day side
end

end


text:
function portfolioassignments!(df::AbstractDataFrame,groupbyvar::Union{Integer,Symbol,Vector{<:Union{Integer,Symbol}}},var::Union{Integer,Symbol},portname::Symbol;nysebrkpts=false)
   df[portname] = allowmissing(fill(99,size(df,1)))
   for subdf in groupby(df,groupbyvar)
      if nysebrkpts
         dsub = subdf[:EXCHCD .== 1, var]
      else
         dsub = subdf[var]
      end
      filter!(!ismissing,dsub)
      if length(dsub) >= 1
         breaks = quantile(convert(Vector{Float64},dsub),0.0:0.1:1.0)
         subdf[portname] = map(x->searchsortedlast(breaks[1:end-1],x),subdf[var])
         subdf[ismissing.(subdf[var]),portname] = missing
      end
   end
end

# assumes :PERMNO and :DATE exist on both DataFrames
function compustat2crspjoin(crsp::AbstractDataFrame,compustat::AbstractDataFrame,joincols::Union{Integer,Symbol,Vector{<:Union{Integer,Symbol}}})
   df = copy(crsp)
   dfjoin = compustat[vcat(:PERMNO,:DATE,joincols)]

   for j  in joincols
      df[j] = Vector{Union{Float64,Missing}}(fill(missing,(nrow(df),)))
   end
   df[:DDATE] = Vector{Union{Date,Missing}}(fill(missing,(nrow(df),)))
   for subdf in groupby(df,:PERMNO)
      pno = subdf[:PERMNO][1]
      joiners=filter(r->r[:PERMNO] == pno,dfjoin)
      nrj = nrow(joiners)
      if nrj >= 1
         @inbounds for i=1:nrj
            dt = joiners[i,:DATE]
            daterows = dt + Dates.Month(6) .<= subdf[:DATE] .< dt + Dates.Month(18)
            for j in joincols
               subdf[daterows,j] = joiners[i,j]
            end
            subdf[daterows,:DDATE] = dt
         end
      end
   end
   return df
end


text:
using ModelingToolkit, LinearAlgebra, SparseArrays
using Symbolics
using Symbolics: scalarize

# Define the constants for the PDE
const α₂ = 1.0
const α₃ = 1.0
const β₁ = 1.0
const β₂ = 1.0
const β₃ = 1.0
const r₁ = 1.0
const r₂ = 1.0
const _DD = 100.0
const γ₁ = 0.1
const γ₂ = 0.1
const γ₃ = 0.1
const N = 8
const X = reshape([i for i in 1:N for j in 1:N], N, N)
const Y = reshape([j for i in 1:N for j in 1:N], N, N)
const α₁ = 1.0 .* (X .>= 4 * N / 5)

const Mx = Tridiagonal([1.0 for i in 1:(N - 1)], [-2.0 for i in 1:N],
                       [1.0 for i in 1:(N - 1)])
const My = copy(Mx)
Mx[2, 1] = 2.0
Mx[end - 1, end] = 2.0
My[1, 2] = 2.0
My[end, end - 1] = 2.0

# Define the initial condition as normal arrays
@variables du[1:N, 1:N, 1:3] u[1:N, 1:N, 1:3] MyA[1:N, 1:N] AMx[1:N, 1:N] DA[1:N, 1:N]

du, u, MyA, AMx, DA = scalarize.((du, u, MyA, AMx, DA))
@show typeof.((du, u, MyA, AMx, DA))

# Define the discretized PDE as an ODE function
function f(du, u, p, t)
    A = @view u[:, :, 1]
    B = @view u[:, :, 2]
    C = @view u[:, :, 3]
    dA = @view du[:, :, 1]
    dB = @view du[:, :, 2]
    dC = @view du[:, :, 3]
    mul!(MyA, My, A)
    mul!(AMx, A, Mx)
    @. DA = _DD * (MyA + AMx)
    @. dA = DA + α₁ - β₁ * A - r₁ * A * B + r₂ * C
    @. dB = α₂ - β₂ * B - r₁ * A * B + r₂ * C
    @. dC = α₃ - β₃ * C + r₁ * A * B - r₂ * C
end

f(du, u, nothing, 0.0)

multithreadedf = eval(ModelingToolkit.build_function(du, u, fillzeros = true,
                                                     parallel = ModelingToolkit.MultithreadedForm())[2])

MyA = zeros(N, N);
AMx = zeros(N, N);
DA = zeros(N, N);
# Loop to catch syncronization issues
for i in 1:100
    _du = rand(N, N, 3)
    _u = rand(N, N, 3)
    multithreadedf(_du, _u)
    _du2 = copy(_du)
    f(_du2, _u, nothing, 0.0)
    @test _du ≈ _du2
end

#=
jac = sparse(ModelingToolkit.jacobian(vec(du),vec(u)))
fjac = eval(ModelingToolkit.build_function(jac,u,parallel=ModelingToolkit.SerialForm())[2])
multithreadedfjac = eval(ModelingToolkit.build_function(jac,u,parallel=ModelingToolkit.MultithreadedForm())[2])

u = rand(N,N,3)
J = similar(jac,Float64)
fjac(J,u)

J2 = similar(jac,Float64)
multithreadedfjac(J2,u)
@test J ≈ J2

using FiniteDiff
J3 = Array(similar(jac,Float64))
FiniteDiff.finite_difference_jacobian!(J2,(du,u)->f!(du,u,nothing,nothing),u)
maximum(J2 .- Array(J)) < 1e-5
=#

jac = ModelingToolkit.sparsejacobian(vec(du), vec(u))
serialjac = eval(ModelingToolkit.build_function(vec(jac), u)[2])
multithreadedjac = eval(ModelingToolkit.build_function(vec(jac), u,
                                                       parallel = ModelingToolkit.MultithreadedForm())[2])

MyA = zeros(N, N)
AMx = zeros(N, N)
DA = zeros(N, N)
_du = rand(N, N, 3)
_u = rand(N, N, 3)

f(_du, _u, nothing, 0.0)
multithreadedf(_du, _u)

#=
using BenchmarkTools
@btime f(_du,_u,nothing,0.0)
@btime multithreadedf(_du,_u)

_jac = similar(jac,Float64)
@btime serialjac(_jac,_u)
@btime multithreadedjac(_jac,_u)
=#


text:
using HoJBot
using Test

@testset "HoJBot.jl" begin
    include("test_utils.jl")
    include("test_discourse.jl")
    include("test_ig.jl")
    include("test_mod.jl")
end


text:
# -*- coding: utf-8 -*-
using Catlab
using ModelingToolkit
import ModelingToolkit: Constant
using Test
using Petri
using SemanticModels.PetriModels

using Catlab.WiringDiagrams
using Catlab.Doctrines
import Catlab.Doctrines.⊗
import Catlab.Graphics: to_graphviz, to_tikz
import Catlab.Graphics.Graphviz: run_graphviz
⊗(a::WiringDiagram, b::WiringDiagram) = otimes(a,b)
import Base: ∘
∘(a::WiringDiagram, b::WiringDiagram) = compose(b, a)
⊚(a,b) = b ∘ a

# Generators
S, E, I, R, D = Ob(FreeSymmetricMonoidalCategory, :S, :E, :I, :R, :D)

# +
infecting = Hom(:infection, S ⊗ I, I⊗I)

# -
# +
inf  = to_wiring_diagram(infecting)
expo = to_wiring_diagram(Hom(:exposure, S⊗I, E⊗I))
rec  = to_wiring_diagram(Hom(:recovery, I,   R))
wan  = to_wiring_diagram(Hom(:waning,   R,   S))

# -
# +
si    = to_wiring_diagram(Hom(:infection,   S⊗I, I⊗I))
se    = to_wiring_diagram(Hom(:exposure,    S⊗I, E⊗I))
prog  = to_wiring_diagram(Hom(:progression, E,   I))
fatal = to_wiring_diagram(Hom(:die,  I, D))
rip   = to_wiring_diagram(Hom(:rest, D, D))

sir    = si    ⊚ (rec   ⊗ rec)
seir   = se    ⊚ (prog  ⊗ rec)
seirs  = seir  ⊚ (wan   ⊗ wan)
seird  = seir  ⊚ (fatal ⊗ to_wiring_diagram(Hom(:id, R, R)))
seirds = seird ⊚ (rip   ⊗ wan)

model(PetriModel, sir)

models = [sir, seir, seirs, seird, seirds]
nets   = model.(PetriModel, models)

modelnames = ["sir",
         "seir",
         "seirs",
         "seird",
         "seirds",
         ]
map(zip(nets, modelnames)) do (net, name)
    Δ = reverse(net.model.Δ)
    println("Model: $name\n  Connections:")
    println("    $Δ\n")
    Δ
end

function writesvg(f::Union{IO,AbstractString}, d::WiringDiagram)
    write(f, to_graphviz(d, labels=true)
          |>g->run_graphviz(g, format="svg")
          )
end

writesvg(x) = writesvg(x...)

try
    mkpath("img")
    writesvg.(zip(["img/$f.svg"
                   for f in modelnames],
                  models))
catch e
    error("Could not write images, does ./img exist?\n $e")
end
println("Model Drawings are in ./img")


text:
using Documenter, HMCUtilities

makedocs(
    modules = [HMCUtilities],
    format = :html,
    checkdocs = :exports,
    sitename = "HMCUtilities.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/salilab/HMCUtilities.jl.git",
)


text:
import Base: IndexStyle, eachindex, eltype, fill!, first, isempty, iterate, lastindex, length, ndims, ones, similar, size, sizeof, stride, strides, zero
using AutoGrad: AutoGrad, @primitive1, ungetindex
using Base: to_shape, DimOrInd

# AbstractArray interface

eachindex(a::KnetArray) = (1:length(a))
eltype(::KnetArray{T}) where {T}=T
eltype(::Type{KnetArray{T}}) where {T} = T
eltype(::Type{KnetArray{T,n}}) where {T,n} = T
lastindex(a::KnetArray) = length(a)
lastindex(a::KnetArray,d) = size(a,d)
fill!(a::KnetArray{T},x) where {T}=(a[:] .= T(x);a)
first(a::KnetArray) = a[1]
# AutoGrad leaves `first` as a compound proc calling start which doesn't work with KnetArrays
@primitive1  first(x::KnetArray),dy,y  AutoGrad.ungetindex(x,dy,1)
isempty(a::KnetArray) = (0==length(a))
length(a::KnetArray)=prod(size(a))

ndims(a::KnetArray{T,N}) where {T,N}=N
ones(a::KnetArray{T}) where {T}=fill!(similar(a),one(T))

similar(a::KnetArray{T}) where {T}                             = similar(a, T)
similar(a::KnetArray, ::Type{T}) where {T}                     = similar(a, T, size(a))
similar(a::KnetArray{T}, dims::Tuple) where {T}                = similar(a, T, to_shape(dims))
similar(a::KnetArray{T}, dims::DimOrInd...) where {T}          = similar(a, T, to_shape(dims))
similar(a::KnetArray, ::Type{T}, dims::DimOrInd...) where {T}  = similar(a, T, to_shape(dims))
similar(a::KnetArray, ::Type{T}, dims::Dims{N}) where {T,N}    = KnetArray{T,N}(undef,dims)

size(a::KnetArray)=a.dims
size(a::KnetArray{T,N},i::Integer) where {T,N}=(if i>N; 1; else; size(a)[i]; end)
sizeof(a::KnetArray{T}) where {T} = length(a) * sizeof(T)
stride(a::KnetArray{T,N},i::Integer) where {T,N}=(if i>N; length(a); else; s=1; for n=1:(i-1); s*=size(a,n); end; s; end)
strides(a::KnetArray{T,N}) where {T,N}=ntuple(n->stride(a,n), N)
zero(a::KnetArray{T}) where {T}=fill!(similar(a),zero(T))

iterate(A::KnetArray,I...)=error() # throw(MethodError(iterate,A,I...))
IndexStyle(::Type{<:KnetArray})=IndexLinear()


text:
# This file is part of the TaylorIntegration.jl package; MIT licensed

"""
    stabilitymatrix!(eqsdiff!, t, x, δx, dδx, jac, _δv, params[, jacobianfunc!=nothing])

Updates the matrix `jac::Matrix{Taylor1{U}}` (linearized equations of motion)
computed from the equations of motion (`eqsdiff!`), at time `t` at `x`; `x` is
of type `Vector{Taylor1{U}}`, where `U<:Number`. `δx`, `dδx` and `_δv` are
auxiliary arrays of type `Vector{TaylorN{Taylor1{U}}}` to avoid allocations.
Optionally, the user may provide a Jacobian function `jacobianfunc!` to compute
`jac`. Otherwise, `jac` is computed via automatic differentiation using
`TaylorSeries.jl`.

"""
function stabilitymatrix!(eqsdiff!, t::Taylor1{T}, x::Vector{Taylor1{U}},
        δx::Vector{TaylorN{Taylor1{U}}}, dδx::Vector{TaylorN{Taylor1{U}}},
        jac::Matrix{Taylor1{U}}, _δv::Vector{TaylorN{Taylor1{U}}}, params,
        jacobianfunc! =nothing) where {T<:Real, U<:Number}

    if isa(jacobianfunc!, Nothing)
        # Set δx equal to current value of x plus 1st-order variations
        for ind in eachindex(δx)
            @inbounds δx[ind] = x[ind] + _δv[ind]
        end
        # Equations of motion
        eqsdiff!(dδx, δx, params, t)
        TaylorSeries.jacobian!(jac, dδx)
    else
        jacobianfunc!(jac, x, params, t)
    end
    nothing
end

# Modified from `cgs` and `mgs`, obtained from:
# http://nbviewer.jupyter.org/url/math.mit.edu/~stevenj/18.335/Gram-Schmidt.ipynb
# Classical Gram–Schmidt (Trefethen algorithm 7.1), implemented in the simplest way
# (We could make it faster by unrolling loops to avoid temporaries arrays etc.)
function classicalGS!(A, Q, R, aⱼ, qᵢ, vⱼ)
    m,n = size(A)
    fill!(R, zero(eltype(A)))
    for j = 1:n
        # aⱼ = A[:,j]
        for ind = 1:m
            @inbounds aⱼ[ind] = A[ind,j]
            @inbounds vⱼ[ind] = aⱼ[ind]
        end
        # vⱼ = copy(aⱼ) # use copy so that modifying vⱼ doesn't change aⱼ
        for i = 1:j-1
            # qᵢ = Q[:,i]
            for ind = 1:m
                @inbounds qᵢ[ind] = Q[ind,i]
            end
            @inbounds R[i,j] = dot(qᵢ, aⱼ)
            # vⱼ -= R[i,j] * qᵢ
            @inbounds for ind = 1:m
                vⱼ[ind] -= R[i,j] * qᵢ[ind]
            end
        end
        @inbounds R[j,j] = norm(vⱼ)
        # Q[:,j] = vⱼ / R[j,j]
        for ind = 1:m
            @inbounds Q[ind,j] = vⱼ[ind] / R[j,j]
        end
    end
    return nothing
end

# Modified Gram–Schmidt (Trefethen algorithm 8.1); see also
# http://nbviewer.jupyter.org/url/math.mit.edu/~stevenj/18.335/Gram-Schmidt.ipynb
function modifiedGS!(A, Q, R, aⱼ, qᵢ, vⱼ)
    m,n = size(A)
    fill!(R, zero(eltype(A)))
    for j = 1:n
        # aⱼ = A[:,j]
        for ind = 1:m
            @inbounds aⱼ[ind] = A[ind,j]
            @inbounds vⱼ[ind] = aⱼ[ind]
        end
        # vⱼ = copy(aⱼ) # use copy so that modifying vⱼ doesn't change aⱼ
        for i = 1:j-1
            # qᵢ = Q[:,i]
            for ind = 1:m
                @inbounds qᵢ[ind] = Q[ind,i]
            end
            @inbounds R[i,j] = dot(qᵢ, vⱼ) # ⟵ NOTICE: mgs has vⱼ, clgs has aⱼ
            # vⱼ -= R[i,j] * qᵢ
            @inbounds for ind = 1:m
                vⱼ[ind] -= R[i,j] * qᵢ[ind]
            end
        end
        @inbounds R[j,j] = norm(vⱼ)
        # Q[:,j] = vⱼ / R[j,j]
        for ind = 1:m
            @inbounds Q[ind,j] = vⱼ[ind] / R[j,j]
        end
    end
    return nothing
end

@doc doc"""
    lyap_jetcoeffs!(t, x, dx, jac, varsaux)

Similar to [`jetcoeffs!`](@ref) for the calculation of the Lyapunov spectrum.
Updates *only* the elements of `x` which correspond to the solution of the 1st-order
variational equations ``\dot{\xi}=J \cdot \xi``, where ``J`` is the Jacobian
matrix, i.e., the linearization of the equations of motion. `jac` is the Taylor
expansion of ``J`` wrt the independent variable, around the current initial
condition. `varsaux` is an auxiliary array of type `Array{eltype(jac),3}` to
avoid allocations. Calling this method assumes that `jac` has been computed
previously using [`stabilitymatrix!`](@ref).

"""
function lyap_jetcoeffs!(t::Taylor1{T}, x::AbstractVector{Taylor1{S}},
        dx::AbstractVector{Taylor1{S}}, jac::Matrix{Taylor1{S}},
        varsaux::Array{Taylor1{S},3}) where {T <: Real, S <: Number}
    order = t.order
    # `dofrange` behaves like 1:dof, where `dof = size(jac, 1)`. Used to initialize `b` and `li`
    dofrange = axes(jac, 1)
    # The `view` below allows us to obtain CartesianIndices from `varsaux` when calling `eachindex(b)`
    # `varsaux` is an array with size n×n×n
    b = view(varsaux, dofrange, dofrange, dofrange)
    # We use `li` to map between cartesian and linear indices
    li = LinearIndices((dofrange, dofrange))

    # 0-th order evaluation of variational equations `dx = jac * x`
    # Initialize LHS of variational equations at zero
    dx .= zero(x[1])
    # Compute 0-th Taylor coefficients of matrix product `jac * x` and save into `dx`
    for i in eachindex(b)
        varsaux[i] = Taylor1(constant_term(jac[ i[1], i[3] ]) * constant_term(x[ li[i[3], i[2]] ]), order)
        lin_indx = li[ i[1], i[2]] # map from cartesian to linear index
        dx[lin_indx] = Taylor1(constant_term(dx[lin_indx]) + constant_term(varsaux[i]), order)
    end
    # Recursion relations, 0-th order
    for __idx = eachindex(x)
        (x[__idx]).coeffs[2] = (dx[__idx]).coeffs[1]
    end

    # Compute Taylor coefficients of variational equations `dx = jac * x` up to order `order`
    for ord = 1:order - 1
        ordnext = ord + 1
        # Compute `ord`-th Taylor coefficients of matrix product `jac * x` and save into `dx`
        for i in eachindex(b)
            TaylorSeries.mul!(varsaux[i], jac[ i[1], i[3] ], x[ li[i[3], i[2]] ], ord)
            lin_indx = li[ i[1], i[2]] # map from cartesian to linear index
            TaylorSeries.add!(dx[lin_indx], dx[lin_indx], varsaux[i], ord)
        end
        # Recursion relations, `ord`-th order
        for __idx = eachindex(x)
            (x[__idx]).coeffs[ordnext + 1] = (dx[__idx]).coeffs[ordnext] / ordnext
        end
    end

    return nothing
end

"""
    lyap_taylorstep!(f!, t, x, dx, xaux, δx, dδx, jac, t0, t1, order, abstol, _δv, varsaux, params[, jacobianfunc!])

Similar to [`taylorstep!`](@ref) for the calculation of the Lyapunov spectrum.
`jac` is the Taylor expansion (wrt the independent variable) of the
linearization of the equations of motion, i.e, the Jacobian. `xaux`, `δx`, `dδx`,
`varsaux` and `_δv` are auxiliary vectors, and `params` define the parameters
of the ODEs. Optionally, the user may provide a
Jacobian function `jacobianfunc!` to compute `jac`. Otherwise, `jac` is computed
via automatic differentiation using `TaylorSeries.jl`.

"""
function lyap_taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}},
        dx::Vector{Taylor1{U}}, xaux::Vector{Taylor1{U}},
        δx::Array{TaylorN{Taylor1{U}},1}, dδx::Array{TaylorN{Taylor1{U}},1},
        jac::Array{Taylor1{U},2}, t0::T, t1::T, order::Int,
        abstol::T, _δv::Vector{TaylorN{Taylor1{U}}}, varsaux::Array{Taylor1{U},3},
        params, parse_eqs::Bool=true, jacobianfunc! =nothing) where {T<:Real, U<:Number}

    # Dimensions of phase-space: dof
    nx = length(x)
    dof = length(δx)

    # Compute the Taylor coefficients associated to trajectory
    __jetcoeffs!(Val(parse_eqs), f!, t, view(x, 1:dof), view(dx, 1:dof), view(xaux, 1:dof), params)

    # Compute stability matrix
    stabilitymatrix!(f!, t, x, δx, dδx, jac, _δv, params, jacobianfunc!)

    # Compute the Taylor coefficients associated to variational equations
    lyap_jetcoeffs!(t, view(x, dof+1:nx), view(dx, dof+1:nx), jac, varsaux)

    # Compute the step-size of the integration using `abstol`
    δt = stepsize(view(x, 1:dof), abstol)
    δt = min(δt, t1-t0)

    return δt
end

"""
    lyap_taylorinteg(f!, q0, t0, tmax, order, abstol[, f!]; maxsteps::Int=500)

Similar to [`taylorinteg`](@ref) for the calculation of the Lyapunov
spectrum. Note that the number of `TaylorN` variables should be set
previously by the user (e.g., by means of `TaylorSeries.set_variables`) and
should be equal to the length of the vector of initial conditions `q0`.
Otherwise, whenever `length(q0) != TaylorSeries.get_numvars()`, then
`lyap_taylorinteg` throws an `AssertionError`. Optionally, the user may provide
a Jacobian function `jacobianfunc!` to evaluate the current value of the Jacobian.
Otherwise, the current value of the Jacobian is computed via automatic
differentiation using `TaylorSeries.jl`.

"""
function lyap_taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T,
        order::Int, abstol::T, params = nothing, jacobianfunc! =nothing;
        maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number}
    # Allocation
    tv = Array{T}(undef, maxsteps+1)
    dof = length(q0)
    xv = Array{U}(undef, dof, maxsteps+1)
    λ = similar(xv)
    λtsum = similar(q0)
    jt = Matrix{U}(I, dof, dof)
    _δv = Array{TaylorN{Taylor1{U}}}(undef, dof)

    # Initial conditions
    @inbounds tv[1] = t0
    @inbounds for ind in eachindex(q0)
        xv[ind,1] = q0[ind]
        λ[ind,1] = zero(U)
        λtsum[ind] = zero(U)
    end
    x0 = vcat(q0, reshape(jt, dof*dof))
    nx0 = length(x0)
    t00 = t0

    # Initialize the vector of Taylor1 expansions
    t = Taylor1(T, order)
    x = Array{Taylor1{U}}(undef, nx0)
    x .= Taylor1.( x0, order )
    @inbounds t[0] = t0

    # If user does not provide Jacobian, check number of TaylorN variables and initialize _δv
    if isa(jacobianfunc!, Nothing)
        @assert get_numvars() == dof "`length(q0)` must be equal to number of variables set by `TaylorN`"
        for ind in eachindex(q0)
            _δv[ind] = one(x[1])*TaylorN(Taylor1{U}, ind, order=1)
        end
    end

    #Allocate auxiliary arrays
    dx = Array{Taylor1{U}}(undef, nx0)
    xaux = Array{Taylor1{U}}(undef, nx0)
    δx = Array{TaylorN{Taylor1{U}}}(undef, dof)
    dδx = Array{TaylorN{Taylor1{U}}}(undef, dof)
    jac = Array{Taylor1{U}}(undef, dof, dof)
    varsaux = Array{Taylor1{U}}(undef, dof, dof, dof)
    fill!(jac, zero(x[1]))
    QH = Array{U}(undef, dof, dof)
    RH = Array{U}(undef, dof, dof)
    aⱼ = Array{U}(undef, dof )
    qᵢ = similar(aⱼ)
    vⱼ = similar(aⱼ)

    # Determine if specialized jetcoeffs! method exists
    parse_eqs = _determine_parsing!(parse_eqs, f!, t, view(x, 1:dof), view(dx, 1:dof), params)

    # Integration
    nsteps = 1
    while t0 < tmax
        δt = lyap_taylorstep!(f!, t, x, dx, xaux, δx, dδx, jac, t0, tmax,
            order, abstol, _δv, varsaux, params, parse_eqs, jacobianfunc!)
        evaluate!(x, δt, x0) # Update x0
        for ind in eachindex(jt)
            @inbounds jt[ind] = x0[dof+ind]
        end
        modifiedGS!( jt, QH, RH, aⱼ, qᵢ, vⱼ )
        t0 += δt
        @inbounds t[0] = t0
        tspan = t0-t00
        nsteps += 1
        @inbounds tv[nsteps] = t0
        @inbounds for ind in eachindex(q0)
            xv[ind,nsteps] = x0[ind]
            λtsum[ind] += log(RH[ind,ind])
            λ[ind,nsteps] = λtsum[ind]/tspan
        end
        for ind in eachindex(QH)
            @inbounds x0[dof+ind] = QH[ind]
        end
        x .= Taylor1.( x0, order )
        if nsteps > maxsteps
            @warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end
    end

    return view(tv,1:nsteps),  view(transpose(xv),1:nsteps,:),  view(transpose(λ),1:nsteps,:)
end

function lyap_taylorinteg(f!, q0::Array{U,1}, trange::AbstractVector{T},
        order::Int, abstol::T, params = nothing, jacobianfunc! = nothing;
        maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number}
    # Allocation
    nn = length(trange)
    dof = length(q0)
    xv = Array{U}(undef, dof, nn)
    fill!(xv, U(NaN))
    λ = Array{U}(undef, dof, nn)
    fill!(λ, U(NaN))
    λtsum = similar(q0)
    jt = Matrix{U}(I, dof, dof)
    _δv = Array{TaylorN{Taylor1{U}}}(undef, dof)

    # Initial conditions
    @inbounds for ind in eachindex(q0)
        xv[ind,1] = q0[ind]
        λ[ind,1] = zero(U)
        λtsum[ind] = zero(U)
    end

    # Initialize the vector of Taylor1 expansions
    t = Taylor1(T, order)
    x0 = vcat(q0, reshape(jt, dof*dof))
    q1 = similar(q0)
    nx0 = length(x0)
    x = Array{Taylor1{U}}(undef, nx0)
    x .= Taylor1.( x0, order )
    @inbounds t[0] = trange[1]
    @inbounds t0, t1, tmax = trange[1], trange[2], trange[end]
    t00 = trange[1]
    tspan = zero(T)

    # If user does not provide Jacobian, check number of TaylorN variables and initialize _δv
    if isa(jacobianfunc!, Nothing)
        @assert get_numvars() == dof "`length(q0)` must be equal to number of variables set by `TaylorN`"
        for ind in eachindex(q0)
            _δv[ind] = one(x[1])*TaylorN(Taylor1{U}, ind, order=1)
        end
    end

    #Allocate auxiliary arrays
    dx = Array{Taylor1{U}}(undef, nx0)
    xaux = Array{Taylor1{U}}(undef, nx0)
    δx = Array{TaylorN{Taylor1{U}}}(undef, dof)
    dδx = Array{TaylorN{Taylor1{U}}}(undef, dof)
    jac = Array{Taylor1{U}}(undef, dof, dof)
    varsaux = Array{Taylor1{U}}(undef, dof, dof, dof)
    fill!(jac, zero(x[1]))
    QH = Array{U}(undef, dof, dof)
    RH = Array{U}(undef, dof, dof)
    aⱼ = Array{U}(undef, dof )
    qᵢ = similar(aⱼ)
    vⱼ = similar(aⱼ)

    # Determine if specialized jetcoeffs! method exists
    parse_eqs = _determine_parsing!(parse_eqs, f!, t, view(x, 1:dof), view(dx, 1:dof), params)

    # Integration
    iter = 2
    nsteps = 1
    while t0 < tmax
        δt = lyap_taylorstep!(f!, t, x, dx, xaux, δx, dδx, jac, t0, tmax,
            order, abstol, _δv, varsaux, params, parse_eqs, jacobianfunc!)
        evaluate!(x, δt, x0) # Update x0
        tnext = t0+δt
        # # Evaluate solution at times within convergence radius
        while t1 < tnext
            evaluate!(x[1:dof], t1-t0, q1)
            @inbounds xv[:,iter] .= q1
            for ind in eachindex(jt)
                @inbounds jt[ind] = evaluate(x[dof+ind], δt)
            end
            modifiedGS!( jt, QH, RH, aⱼ, qᵢ, vⱼ )
            tspan = t1-t00
            @inbounds for ind in eachindex(q0)
                λ[ind,iter] = (λtsum[ind]+log(RH[ind,ind]))/tspan
            end
            iter += 1
            @inbounds t1 = trange[iter]
        end
        if δt == tmax-t0
            @inbounds xv[:,iter] .= x0[1:dof]
            for ind in eachindex(jt)
                @inbounds jt[ind] = x0[dof+ind]
            end
            modifiedGS!( jt, QH, RH, aⱼ, qᵢ, vⱼ )
            tspan = tmax-t00
            @inbounds for ind in eachindex(q0)
                λ[ind,iter] = (λtsum[ind]+log(RH[ind,ind]))/tspan
            end
            break
        end

        for ind in eachindex(jt)
            @inbounds jt[ind] = x0[dof+ind]
        end
        modifiedGS!( jt, QH, RH, aⱼ, qᵢ, vⱼ )

        t0 = tnext
        @inbounds t[0] = t0
        nsteps += 1
        @inbounds for ind in eachindex(q0)
            λtsum[ind] += log(RH[ind,ind])
        end
        for ind in eachindex(QH)
            @inbounds x0[dof+ind] = QH[ind]
        end
        x .= Taylor1.( x0, order )
        if nsteps > maxsteps
            @warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end
    end

    return transpose(xv), transpose(λ)
end


text:
function streamplot(;x::AbstractVector = Int[], y::AbstractVector = Int[], group::AbstractVector = Int[])

    v = areaplot(x = x, y = y, group = group, stacked = true)
    v.scales[1]._type = "ordinal"
    v.width = 900
    v.height = 400
    v.marks[1].from.transform[1].x["offset"] = "center"

    return v

end

max_stars_repo_name: jpthor/vega


text:
function heat_equation(cur_reactor::AbstractReactor)
  cur_G_T = K_DV(cur_reactor)

  cur_G_T *= cur_reactor.sigma_v ^ (1/3.2)

  cur_equation = Equation(1, 0, -1, cur_G_T)

  cur_equation
end