text:
using RCall, CSV, DataFrames

########################################################################
#Chronic Infected (Figure 6)
########################################################################

Det_Homo = CSV.read("noStochFigure6Data.csv")
Det_Hetero = CSV.read("VirusFigure6_Det_Hetero.csv")
Stoch_Homo = CSV.read("StochFigure6Data.csv")

@rput Det_Homo
@rput Det_Hetero
@rput Stoch_Homo
R"""
library(ggplot2)
library(ggpubr)

 commonFigureOptions <- list(
   theme_pubr(border=TRUE),
   xlab("Time (hours)"),
   labs(fill=expression(paste("IFN", beta, " (nM)"))),
   theme_pubr(border=TRUE,base_size = 14),
   scale_x_continuous(expand=c(0,0)),
   scale_y_continuous(expand=c(0,0)),
   xlab(bquote("Interferon Degradation" ~ tau[7])),
   ylab(bquote("JAK/STAT Activity" ~ k[cat8])),
   theme(plot.title = element_text(hjust = 0.5),aspect.ratio = 1) )

p1 <- ggplot(Det_Homo, aes(kcat8, tau7, fill=IFN)) +
  geom_raster(aes(fill=IFN)) +
  ggtitle("ISD \n Deterministic \n Homogeneous") +
  scale_fill_distiller(palette = "Spectral",
   guide = guide_colorbar(frame.colour = "black",
   ticks.colour = "black")) +
  commonFigureOptions

p2 <- ggplot(Det_Hetero, aes(kcat8, tau7, fill=IFN)) +
  geom_raster(aes(fill=IFN)) +
  ggtitle("ISD \n Deterministic \n Heterogeneous") +
  scale_fill_distiller(palette = "Spectral",
   guide = guide_colorbar(frame.colour = "black",
   ticks.colour = "black")) +
  commonFigureOptions

Stoch_Homo$IFN[Stoch_Homo$IFN<0] <- 0

p3 <- ggplot(Stoch_Homo, aes(kcat8, tau7, fill=IFN)) +
  geom_raster(aes(fill=IFN)) +
  ggtitle("ISD \n Stochastic \n Homogeneous") +
  scale_fill_distiller(palette = "Spectral",
   guide = guide_colorbar(frame.colour = "black",
   ticks.colour = "black"),
   limits = c(0,max(Stoch_Homo$IFN))) +
  commonFigureOptions

figure <- ggarrange(p1, p2, p3,
                    labels = "AUTO",
                    legend = "right",
                    align = "hv",
                    ncol = 3, nrow = 1)

ggsave("../Figures/Figure6_All.pdf",width = 15,height=5)
"""

########################################################################
#Stochastic IFN (Figure 7)
########################################################################

R"""
library(tidyverse)
library(ggpubr)
library(scales)

df <- read.csv("VirusFigure7DataHomo50cells.csv")

#And that's when I discovered how to use tidyverse

df_tidy_median <- df %>%
  pivot_longer(c(Healthy,Infected,Dead),names_to = "CellState",values_to = "CellPercent") %>%
  filter(Percent==0.2) %>%
  group_by(Time,CellState) %>%
  summarise(medianPercentCells = median(CellPercent),
            q25 = quantile(CellPercent,0.25),
            q75 = quantile(CellPercent,0.75))

df_tidy_median$CellState = factor(df_tidy_median$CellState,ordered=TRUE,levels = c("Healthy","Infected","Dead"))

color_list <- c("#377EB8","#4DAF4A","#E41A1C")

p1 <- ggplot(df_tidy_median, aes(x=Time, y=medianPercentCells)) +
  geom_line(aes(x=Time, y=medianPercentCells, color=CellState)) +
  geom_ribbon(aes(ymin=q25, ymax=q75, fill=CellState),alpha=0.3) +
  scale_fill_manual(values=color_list) +
  scale_color_manual(values=color_list) +
  scale_x_continuous(breaks=seq(0, 48, 12)) +
  scale_y_continuous(labels = scales::percent) +
  xlab("Time (hours)") +
  ylab("Cell Percent") +
  theme_pubr(border=TRUE)

df_tidy_box <- df %>%
  filter(Time==48)

p2 <-ggplot(df_tidy_box, aes(x=as.factor(Percent), y=Dead)) +
  stat_boxplot(geom = "errorbar", width = 0.2,size=0.1) +
  geom_boxplot(fatten=1, fill = "#f6b9ba") +
  #geom_jitter(shape=16, position=position_jitter(0.2)) +
  scale_y_continuous(labels = scales::percent,limits=c(0,1)) +
  scale_x_discrete(labels = percent((0:10)/10)) +
  xlab(expression("IFN"*beta*" Producing Cells")) +
  ylab("Dead Cells") +
  theme_pubr(border=TRUE)

df <- read.csv("VirusFigure7DataHetero50cells.csv")

df_tidy_median <- df %>%
  pivot_longer(c(Healthy,Infected,Dead),names_to = "CellState",values_to = "CellPercent") %>%
  filter(Percent==0.2) %>%
  group_by(Time,CellState) %>%
  summarise(medianPercentCells = median(CellPercent),
            q25 = quantile(CellPercent,0.25),
            q75 = quantile(CellPercent,0.75))

df_tidy_median$CellState = factor(df_tidy_median$CellState,ordered=TRUE,levels = c("Healthy","Infected","Dead"))

color_list <- c("#377EB8","#4DAF4A","#E41A1C")

p3 <- ggplot(df_tidy_median, aes(x=Time, y=medianPercentCells)) +
  geom_line(aes(x=Time, y=medianPercentCells, color=CellState)) +
  geom_ribbon(aes(ymin=q25, ymax=q75, fill=CellState),alpha=0.3) +
  scale_fill_manual(values=color_list) +
  scale_color_manual(values=color_list) +
  scale_x_continuous(breaks=seq(0, 48, 12)) +
  scale_y_continuous(labels = scales::percent) +
  xlab("Time (hours)") +
  ylab("Cell Percent") +
  theme_pubr(border=TRUE)

df_tidy_box <- df %>%
  filter(Time==48)

p4 <-ggplot(df_tidy_box, aes(x=as.factor(Percent), y=Dead)) +
  stat_boxplot(geom = "errorbar", width = 0.2,size=0.1) +
  geom_boxplot(fatten=1, fill = "#f6b9ba") +
  scale_y_continuous(labels = scales::percent,limits=c(0,1)) +
  scale_x_discrete(labels = percent((0:10)/10)) +
  xlab(expression("IFN"*beta*" Producing Cells")) +
  ylab("Dead Cells") +
  theme_pubr(border=TRUE)

figure <- ggarrange(p1, p2, p3, p4,
                    labels = "AUTO",
                    common.legend = TRUE, legend = "right",
                    widths = c(3,6),
                    align = "hv",
                    ncol = 2, nrow = 2)

ggsave("../Figures/Figure7.pdf",width = 9,height=6,units="in")
"""


text:
type Classical3dSquareSimulator{T} <: Classical3dSimulator{T}
	dimM::Int
	wholeiteration::Int
	countiteration::Int
	normalizationFactor::Array{T}
	data4energy::Array{T}
	lattice::Classical3dSquareLattice{T}

	function Classical3dSquareSimulator{T}(lattice::Classical3dSquareLattice{T}, dimM::Int, wholeiteration::Int)
		this = new{T}()
		this.dimM = dimM
		this.wholeiteration = wholeiteration
		this.lattice = lattice
		this.normalizationFactor = Array{T}(wholeiteration)
		this.data4energy = Array{T}(wholeiteration)
		initializeCount!(this)
		return this
	end
end

Classical3dSquareSimulator{T}(lattice::Classical3dSquareLattice{T}, dimM::Int, wholeiteration::Int) =
Classical3dSquareSimulator{T}(lattice,dimM,wholeiteration)

#####
function (simulator::Classical3dSquareSimulator){T}(temperature::T = getTemperature(simulator), externalfield::T = getExternalfield(simulator))
	setEnvParameters!(simulator, temperature, externalfield)
	initializeCount!(simulator)
	while true
		whichIsTensorT = 1
		countUp!(simulator)
		renormalizeX!(simulator; whichIsTensorT = whichIsTensorT)
		countUp!(simulator)
		renormalizeY!(simulator; whichIsTensorT = whichIsTensorT)
		countUp!(simulator)
		renormalizeZ!(simulator; whichIsTensorT = whichIsTensorT)
		whichIsTensorT = 2
		countUp!(simulator)
		renormalizeX!(simulator; whichIsTensorT = whichIsTensorT)
		countUp!(simulator)
		renormalizeY!(simulator; whichIsTensorT = whichIsTensorT)
		countUp!(simulator)
		renormalizeZ!(simulator; whichIsTensorT = whichIsTensorT)
		if isDone(simulator)
			break
		end
	end
	freeenergy = getFreeEnergy(simulator)
	magnetization = getExpectationValue(simulator)
	if isPottsModel(simulator)
		NOS = convert(T,getStates(simulator))
		magnetization = (NOS * magnetization - one(T))/(NOS - one(T))
	end

	return freeenergy, magnetization
end

#---
# test Function
function testClassical3dSquareSimulator()
	lattice = testClassical3dSquareLattice()
	dimM = 5
	wholeiteration = 12
	simulator = Classical3dSquareSimulator(lattice, dimM, wholeiteration)
	return simulator
end

ext: jl


text:
module Utils
    using JSON, ..Metrics, ..Mutation, ..Recombination, ..FunctionsCollection
    function get_algo_params()
        file_path = pwd() * "/algo_inputs.json"
        parameters = JSON.parse(read(file_path, String))
        parameters["error_function"] = select_error_function(parameters["error_function"])
        parameters["recombination_type"] = select_recombination_type(parameters["recombination_type"])
        parameters["mutation_type"] = select_mutation_type(parameters["mutation_type"])
        parameters["evaluate_function"] = select_evaluate_function(parameters["evaluate_function"])
        parameters["output_file"] = select_output_file(parameters["output_file"])
        return parameters
    end

    function select_output_file(output_file)
        letters = Set(output_file)
        if length(letters) == 1 && ' ' in letters
            return nothing
        elseif output_file == ""
            return nothing
        else
            return output_file
        end
    end

    function select_evaluate_function(evaluate_function)
        if evaluate_function[end] == '1'
            return FunctionsCollection.function1
        elseif evaluate_function[end] == '2'
            return FunctionsCollection.function2
        elseif evaluate_function[end] == '3'
            return FunctionsCollection.function3
        elseif evaluate_function[end] == '4'
            return FunctionsCollection.function4
        elseif evaluate_function[end] == '5'
            return FunctionsCollection.function5
        else
            throw(ArgumentError("Function to evaluate requested not available. Options are:
                                function1, function2, function3, function4, function5"))
        end
    end

    # Get error function based on possibilities
    function select_error_function(error_function)
        if error_function == "mean_sqr_error"
            return Metrics.mean_sqr_error
        elseif error_function == "mean_absolute_error"
            return Metrics.mean_absolute_error
        elseif error_function == "weighted_mean_abs_error"
            return Metrics.weighted_mean_abs_error
        else
            throw(ArgumentError("Error function requested not available. Options are: mean_sqr_error,
                                mean_absolute_error, weighted_mean_abs_error"))
        end
    end

    # Get recombination based on possibilities
    function select_recombination_type(recombination_type)

        if recombination_type == "intermediary_recombination"
            Recombination.intermediary_recombination
        elseif recombination_type == "discrete_recombination"
            Recombination.discrete_recombination
        else
            throw(ArgumentError("Recombination type requested not available. Options are:
                                intermediary_recombination, discrete_recombination"))
        end
    end

    # Get mutation type based on possibilities
    function select_mutation_type(mutation_type)
        if mutation_type == "uncorr_mutation_onestepsize"
            return Mutation.uncorr_mutation_onestepsize
        elseif mutation_type == "uncorr_mutation_n_stepsize"
            return Mutation.uncorr_mutation_n_stepsize
        else
            throw(ArgumentError("Mutation type requested not available. Options are:
                                uncorr_mutation_onestepsize, uncorr_mutation_n_stepsize"))
        end

    end

    # Linspace like utils
    function generate_linespace(min, max, n_values)
        step = (max-min)/n_values
        return collect(min:step:max)
    end

end


text:
# Invertible network based hyperbolic layers (Lensink et. al, 2019)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: February 2020

export NetworkHyperbolic

"""
    H = NetworkHyperbolic(nx, ny, n_in, batchsize, n_hidden, L, K)

 Create an invertible network based on hyperbolic layers. The network consists of L downsampling
 layers (i.e. scales) using the wavelet transform and `L` upsampling layers to bring the output
 back to the original input dimensions. Each of the `L` scales consists of `K` hyperbolic layers.

 *Input*:

 - `nx`, `ny`, `n_in`, `batchsize`: spatial dimensions, number of channels and batchsize of input tensor

 - `n_hidden`: number of hidden units in residual blocks

 - `L`: number of scales

 - `K`: number of time steps per scale

 *Output*:

 - `H`: invertible hyperbolic network.

 *Usage:*

 - Forward mode: `Y, logdet = H.forward(X)`

 - Inverse mode: `X = H.inverse(Y)`

 - Backward mode: `ΔX, X = H.backward(ΔY, Y)`

 *Trainable parameters:*

 - None in `H` itself

 - Trainable parameters in activation normalization `H.AL` and hyperbolic layers `H.HL[j]`.

 See also: [`AffineLayer`](@ref), [`CouplingLayer!`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct NetworkHyperbolic <: InvertibleNetwork
    AL::AffineLayer
    HL::AbstractArray{HyperbolicLayer, 1}
    logdet::Bool
end

@Flux.functor NetworkHyperbolic

# Constructor
function NetworkHyperbolic(nx::Int64, ny::Int64, n_in::Int64, batchsize::Int64, L::Int64, K::Int64;
        k=3, s=1, p=1, logdet=true, α=1f0, hidden_factor=1, ncenter=1)

    depth = Int(2*(L-1)*K + ncenter)
    AL = AffineLayer(Int(nx/2), Int(ny/2), Int(n_in*4); logdet=logdet)
    HL = Array{HyperbolicLayer}(undef, depth)
    nx = Int(nx/2); ny = Int(ny/2); n_in = Int(n_in*2)  # dimensions after initial wavelet transform

    # Downsampling layers
    count=1
    for i=1:L-1
        if K > 1
            for j=2:K
                HL[count] = HyperbolicLayer(Int(nx/2^(i-1)), Int(ny/2^(i-1)), Int(n_in*4^(i-1)),
                    batchsize, k, s, p; action="same", α=α, hidden_factor=hidden_factor)
                count += 1
            end
        end
        HL[count] = HyperbolicLayer(Int(nx/2^(i-1)), Int(ny/2^(i-1)), Int(n_in*4^(i-1)),
            batchsize, k, s, p; action="down", α=α, hidden_factor=hidden_factor)
        count += 1
    end

    # Middle layers at coarsest scale
    for i=1:ncenter
        HL[count] = HyperbolicLayer(Int(nx/2^(L-1)), Int(ny/2^(L-1)), Int(n_in*4^(L-1)),
            batchsize, k, s, p; action="same", α=α, hidden_factor=hidden_factor)
        count += 1
    end

    # Upsampling layers
    for i=L-1:-1:1
        HL[count] = HyperbolicLayer(Int(nx/2^i), Int(ny/2^i), Int(n_in*4^i),
            batchsize, k, s, p; action="up", α=α, hidden_factor=hidden_factor)
        count += 1
        if K > 1
            for j=2:K
                HL[count] = HyperbolicLayer(Int(nx/2^(i-1)), Int(ny/2^(i-1)), Int(n_in*4^(i-1)), batchsize, k, s, p;
                    action="same", α=α, hidden_factor=hidden_factor)
                count += 1
            end
        end
    end

    return NetworkHyperbolic(AL, HL, logdet)
end

# Forward pass
function forward(X, H::NetworkHyperbolic)
    X = wavelet_squeeze(X)
    X, logdet = H.AL.forward(X)
    X_prev, X_curr = tensor_split(X)
    for j=1:length(H.HL)
        X_prev, X_curr = H.HL[j].forward(X_prev, X_curr)
    end
    X = tensor_cat(X_prev, X_curr)
    X = wavelet_unsqueeze(X)
    return X, logdet
end

# Inverse pass
function inverse(Y, H::NetworkHyperbolic)
    Y = wavelet_squeeze(Y)
    Y_curr, Y_new = tensor_split(Y)
    for j=length(H.HL):-1:1
        Y_curr, Y_new = H.HL[j].inverse(Y_curr, Y_new)
    end
    Y = tensor_cat(Y_curr, Y_new)
    Y = H.AL.inverse(Y)
    Y = wavelet_unsqueeze(Y)
    return Y
end

# Backward pass
function backward(ΔY, Y, H::NetworkHyperbolic)
    ΔY = wavelet_squeeze(ΔY)
    Y = wavelet_squeeze(Y)
    ΔY_curr, ΔY_new = tensor_split(ΔY)
    Y_curr, Y_new = tensor_split(Y)
    for j=length(H.HL):-1:1
        ΔY_curr, ΔY_new, Y_curr, Y_new = H.HL[j].backward(ΔY_curr, ΔY_new, Y_curr, Y_new)
    end
    ΔY = tensor_cat(ΔY_curr, ΔY_new)
    Y = tensor_cat(Y_curr, Y_new)
    ΔY, Y = H.AL.backward(ΔY, Y)
    ΔY = wavelet_unsqueeze(ΔY)
    Y = wavelet_unsqueeze(Y)
    return ΔY, Y
end

# Clear gradients
function clear_grad!(H::NetworkHyperbolic)
    depth = length(H.HL)
    clear_grad!(H.AL)
    for j=1:depth
        clear_grad!(H.HL[j])
    end
end

# Get parameters
function get_params(H::NetworkHyperbolic)
    depth = length(H.HL)
    p = get_params(H.AL)
    for j=1:depth
        p = cat(p, get_params(H.HL[j]); dims=1)
    end
    return p
end

# Put parameters
function put_params!(H::NetworkHyperbolic, Params::Array{Any,1})
    depth = length(H.CL)
    put_params!(H.AL, Params[1:2])
    for j = 1:depth
        put_params!(H.HL[j], Params[2*j+1:2*(j+1)])
    end
end


text:
using StructJuMP

m = StructuredModel()

@variable(m, xh[ASSETS] >= 0)
@variable(m, xb[ASSETS] >= 0)
@variable(m, xs[ASSETS] >= 0)

# Stage 0
sub = StructuredModel(parent=m)
@constraint(sub, (1+Ct)*sum(val[j]*xb[j] for j = ASSETS) <= BUDGET)

for stage = 1:T
    sub = StructuredModel(parent=m) # return a ConstraintRef to a Stage?
    for j in ASSETS
        @constraint(sub, xh[j] == (1+Ret[j])*xh[getParent(sub),j] + xb[j] - xs[j])
        @constraint(sub, (1-Ct)*sum(val[j]*xs[j] for j=ASSETS) == LIABILITY + (1+Ct)*sum(val[j]*xb[j] for j=ASSETS))
    end
end

sub = StructuredModel(parent=m)
@variable(sub, wealth)
@constraint(sub, wealth == sum(val[j]*xh[j] for j=ASSETS))
@constraint(sub, @Exp(wealth) == mu)
@objective(sub, mu - Rho*((wealth*wealth)- mu*mu))

# Exp(x) = sum(ndinNODESET:ndincurrentstage) * Prob[nd] ∗x[nd]


text:
# bayesian inference for gauss-gamma distribution

module GaussGamma
  using DataFrames, CSV, Plots, FreqTables, NamedArrays, Distributions
  using StatsBase, Statistics
  pyplot()

  function main(;test=false)
    # step1: generate prior distribution
    # initial parameters
    mu_0 = 200
    zeta_0 = 1
    alpha_0 = 1
    beta_0 = 2
    # prior gamma distribution of lambda
    lambda_dist = gen_lambda_dist(alpha_0, beta_0)
    draw_dist_pdf(lambda_dist, 0, 5, 0.01,
                  "src/bayesian_inference/prior_lambda_dist.png",
                  test)
    # prior normal distribution of mu
    lambda = rand(lambda_dist)
    mu_dist = gen_mu_dist(mu_0, zeta_0, lambda)
    draw_dist_pdf(mu_dist, 180, 220, 0.1,
                  "src/bayesian_inference/prior_mu_dist.png",
                  test)

    # step2: calculate posterior distribution with not many observation
    data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_200.txt")
    df_200_mm = CSV.read(data_path, DataFrame,
                         header=["date", "time", "ir", "lidar"],
                         delim=' ')
    lidar = df_200_mm.lidar
    samples = sample(lidar, 5)
    println("Samples(N=5): $(samples)")
    println("Mean(N=5): $(sum(samples)/length(samples))")
    println("Std dev(N=5): $(Statistics.std(samples, corrected=false))")
    # calculate parameters
    N = length(samples)
    mu_N = 1.0/(N+beta_0)*sum(samples) + beta_0/(N+beta_0)*mu_0
    zeta_N = N + zeta_0
    alpha_N = N/2 + alpha_0
    beta_N = 0.5*(sum([z^2 for z in samples]) + zeta_0*(mu_0^2) - zeta_N*(mu_N^2)) + beta_0
    println("$(mu_N) $(zeta_N) $(alpha_N) $(beta_N)")
    # posterior gamma distribution of lambda
    lambda_dist = gen_lambda_dist(alpha_N, beta_N)
    draw_dist_pdf(lambda_dist, 0, 5, 0.01,
                  "src/bayesian_inference/posterior_lambda_dist_N5.png",
                  test)
    draw_dist_pdf(lambda_dist, 0, 0.01, 0.00001,
                  "src/bayesian_inference/posterior_lambda_dist_N5_zoom.png",
                  test)
    # posterior normal distribution of mu
    lambda = rand(lambda_dist)
    mu_dist = gen_mu_dist(mu_N, zeta_N, lambda)
    draw_dist_pdf(mu_dist, 180, 220, 0.1,
                  "src/bayesian_inference/posterior_mu_dist_N5.png",
                  test)
    println("Mean(N=5): $(mu_dist.μ)")
    println("Std dev(N=5): $(sqrt(1/lambda))")

    # step3: calculate posterior distribution with a lot of observation
    samples = lidar # all data
    println("Samples(all): $(samples)")
    println("Mean(all): $(sum(samples)/length(samples))")
    println("Std dev(all): $(Statistics.std(samples, corrected=false))")
    # calculate parameters
    N = length(samples)
    mu_N = 1.0/(N+beta_0)*sum(samples) + beta_0/(N+beta_0)*mu_0
    zeta_N = N + zeta_0
    alpha_N = N/2 + alpha_0
    beta_N = 0.5*(sum([z^2 for z in samples]) + zeta_0*(mu_0^2) - zeta_N*(mu_N^2)) + beta_0
    println("$(mu_N) $(zeta_N) $(alpha_N) $(beta_N)")
    # posterior gamma distribution of lambda
    lambda_dist = gen_lambda_dist(alpha_N, beta_N)
    draw_dist_pdf(lambda_dist, 0, 5, 0.01,
                  "src/bayesian_inference/posterior_lambda_dist_all.png",
                  test)
    draw_dist_pdf(lambda_dist, 0.035, 0.05, 0.0001,
                  "src/bayesian_inference/posterior_lambda_dist_all_zoom.png",
                  test)
    # posterior normal distribution of mu
    lambda = rand(lambda_dist)
    mu_dist = gen_mu_dist(mu_N, zeta_N, lambda)
    draw_dist_pdf(mu_dist, 180, 220, 0.1,
                  "src/bayesian_inference/posterior_mu_dist_all.png",
                  test)
    println("Mean(all): $(mu_dist.μ)")
    println("Std dev(all): $(sqrt(1/lambda))")
    # compare with histogram
    # draw histogram
    bin_min_max = maximum(df_200_mm.lidar) - minimum(df_200_mm.lidar)
    histogram(df_200_mm.lidar, bins=bin_min_max, color=:orange,
              label="histogram")
    if test == false
      save_path = joinpath(split(@__FILE__, "src")[1], "src/bayesian_inference/lidar_200mm_histogram.png")
      savefig(save_path)
    end
    # calculated distribution
    z_dist = Normal(mu_dist.μ, sqrt(1/(lambda)))
    xs = range(190, 230, length=Int64(floor(40/0.1)))
    ys = [pdf(z_dist, x) for x in xs]
    plot(xs, ys, label="PDF", color=:blue)
    if test == false
      save_path = joinpath(split(@__FILE__, "src")[1], "src/bayesian_inference/lidar_200mm_calc_dist.png")
      savefig(save_path)
    end
  end

  function gen_lambda_dist(alpha, beta)
    return Gamma(alpha, 1/beta)
  end

  function gen_mu_dist(mu_mean, zeta, lambda)
    return Normal(mu_mean, sqrt(1/(zeta*lambda)))
  end

  function draw_dist_pdf(dist, range_min, range_max, step, save_name, test)
    xs = range(range_min, range_max, length=Int64(floor((range_max-range_min)/step)))
    ys = [pdf(dist, x) for x in xs]
    plot(xs, ys, label="PDF")

    if test == false
      save_path = joinpath(split(@__FILE__, "src")[1], save_name)
      savefig(save_path)
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
abstract type AbstractModel{T} end

function Base.show(io::IO, m::AbstractModel)
    @printf io "Dynamic Stochastic General Equilibrium Model\n"
    @printf io "no. states:             %i\n" n_states(m)
    @printf io "no. anticipated shocks: %i\n" n_anticipated_shocks(m)
    @printf io "data vintage:           %s\n" data_vintage(m)
    @printf io "description:\n %s\n"          description(m)
end

@inline function Base.getindex(m::AbstractModel, i::Integer)
    if i <= (j = length(m.parameters))
        return m.parameters[i]
    else
        return m.steady_state[i-j]
    end
end

# need to define like this so we can disable bounds checking
@inline function Base.getindex(m::AbstractModel, k::Symbol)
    i = m.keys[k]
    @inbounds if i <= (j = length(m.parameters))
        return m.parameters[i]
    else
        return m.steady_state[i-j]
    end
end

@inline function Base.setindex!(m::AbstractModel, value::Number, i::Integer)
    if i <= (j = length(m.parameters))
        param = m.parameters[i]
        param.value = value
        if isa(param, ScaledParameter)
            param.scaledvalue = param.scaling(value)
        end
        return param
    else
        steady_state_param = m.steady_state[i-j]
        steady_state_param.value = value
        return steady_state_param
    end
end

"""
```
setindex!(m::AbstractModel, param::AbstractParameter, i::Integer)
```

If `i`<length(m.parameters), overwrites m.parameters[i] with
param. Otherwise, overwrites m.steady_state[i-length(m.parameters).
"""
@inline function Base.setindex!(m::AbstractModel, param::AbstractParameter, i::Integer)
    if i <= (j = length(m.parameters))
        m.parameters[i] = param
    else
        m.steady_state[i-j] = param
    end
    return param
end

Base.setindex!(m::AbstractModel, value, k::Symbol) = Base.setindex!(m, value, m.keys[k])

"""
```
(<=){T}(m::AbstractModel{T}, p::AbstractParameter{T})
```

Syntax for adding a parameter to a model: m <= parameter.
NOTE: If `p` is added to `m` and length(m.steady_state) > 0, `keys(m)` will not generate the
index of `p` in `m.parameters`.
"""
function (<=){T}(m::AbstractModel{T}, p::AbstractParameter{T})

    if !in(p.key, keys(m.keys))

        new_param_index = length(m.keys) + 1

        # grow parameters and add the parameter
        push!(m.parameters, p)

        # add parameter location to dict
        setindex!(m.keys, new_param_index, p.key)
    else
        # overwrite the previous parameter with the new one
        setindex!(m, p, p.key)
    end
end

"""
```
(<=){T}(m::AbstractModel{T}, ssp::SteadyStateParameter)
```

Add a new steady-state value to the model by appending `ssp` to the `m.steady_state` and
adding `ssp.key` to `m.keys`.
"""
function (<=){T}(m::AbstractModel{T}, ssp::SteadyStateParameter)

    if !in(ssp.key, keys(m.keys))
        new_param_index = length(m.keys) + 1

        # append ssp to steady_state vector
        push!(m.steady_state, ssp)

        # add parameter location to dict
        setindex!(m.keys, new_param_index, ssp.key)
    else
        # overwrite the previous parameter with the new one
        setindex!(m, ssp, ssp.key)
    end
end

Distributions.logpdf(m::AbstractModel) = logpdf(m.parameters)
Distributions.pdf(m::AbstractModel) = exp(logpdf(m))

# Number of anticipated policy shocks
n_anticipated_shocks(m::AbstractModel) = get_setting(m, :n_anticipated_shocks)
n_anticipated_shocks_padding(m::AbstractModel) = get_setting(m, :n_anticipated_shocks_padding)

# Dates, indices, number of periods for each regime
date_presample_start(m::AbstractModel) = get_setting(m, :date_presample_start)
date_mainsample_start(m::AbstractModel) = get_setting(m, :date_mainsample_start)
date_zlb_start(m::AbstractModel) = get_setting(m, :date_zlb_start)

date_presample_end(m::AbstractModel) = Dates.lastdayofquarter(get_setting(m, :date_mainsample_start) - Dates.Month(3))
date_prezlb_end(m::AbstractModel) = Dates.lastdayofquarter(get_setting(m, :date_zlb_start) - Dates.Month(3))
date_mainsample_end(m::AbstractModel) = Dates.lastdayofquarter(get_setting(m, :date_forecast_start) - Dates.Month(3))
date_conditional_end(m::AbstractModel) = get_setting(m, :date_conditional_end)

index_presample_start(m::AbstractModel) = 1
index_mainsample_start(m::AbstractModel) = subtract_quarters(date_mainsample_start(m), date_presample_start(m)) + 1
index_zlb_start(m::AbstractModel) = subtract_quarters(date_zlb_start(m), date_presample_start(m)) + 1
index_forecast_start(m::AbstractModel) = subtract_quarters(date_forecast_start(m), date_presample_start(m)) + 1

"""
```
index_shockdec_start(m::AbstractModel)
```

Returns the index starting from which the shock decomposition is saved, where 1 is the index corresponding to date_mainsample_start(m).
"""
index_shockdec_start(m::AbstractModel) = subtract_quarters(date_shockdec_start(m), date_mainsample_start(m)) + 1

"""
```
index_shockdec_end(m::AbstractModel)
```

Returns the last index for which the shock decomposition is saved, where 1 is the index corresponding to date_mainsample_start(m).
"""
index_shockdec_end(m::AbstractModel) = subtract_quarters(date_shockdec_end(m), date_mainsample_start(m)) + 1

n_presample_periods(m::AbstractModel)   = subtract_quarters(date_mainsample_start(m), date_presample_start(m))
n_prezlb_periods(m::AbstractModel)      = subtract_quarters(date_zlb_start(m), date_mainsample_start(m))
n_zlb_periods(m::AbstractModel)         = subtract_quarters(date_forecast_start(m), date_zlb_start(m))
n_mainsample_periods(m::AbstractModel)  = subtract_quarters(date_forecast_start(m), date_mainsample_start(m))
n_conditional_periods(m::AbstractModel) = subtract_quarters(date_conditional_end(m), date_mainsample_end(m))

inds_presample_periods(m::AbstractModel)  = collect(index_presample_start(m):(index_mainsample_start(m)-1))
inds_prezlb_periods(m::AbstractModel)     = collect(index_mainsample_start(m):(index_zlb_start(m)-1))
inds_zlb_periods(m::AbstractModel)        = collect(index_zlb_start(m):(index_forecast_start(m)-1))
inds_mainsample_periods(m::AbstractModel) = collect(index_mainsample_start(m):(index_forecast_start(m)-1))

# Number of a few things that are useful
n_states(m::AbstractModel)                  = length(m.endogenous_states)
n_states_augmented(m::AbstractModel)        = n_states(m) + length(m.endogenous_states_augmented)
n_shocks_exogenous(m::AbstractModel)        = length(m.exogenous_shocks)
n_shocks_expectational(m::AbstractModel)    = length(m.expected_shocks)
n_observables(m::AbstractModel)             = length(m.observables)
n_pseudo_observables(m::AbstractModel)      = length(m.pseudo_observables)
n_equilibrium_conditions(m::AbstractModel)  = length(m.equilibrium_conditions)
n_parameters(m::AbstractModel)              = length(m.parameters)
n_parameters_steady_state(m::AbstractModel) = length(m.steady_state)
n_parameters_free(m::AbstractModel)         = sum([!α.fixed for α in m.parameters])

"""
```
get_key(m, class, index)
```

Returns the name of the state (`class = :states`), observable (`:obs`),
pseudo-observable (`:pseudo`), or shock (`:shocks` or `:stdshocks`)
corresponding to the given `index`.
"""
function get_key(m::AbstractModel, class::Symbol, index::Int)
    dict = if class == :states
        m.endogenous_states
    elseif class == :obs
        m.observables
    elseif class == :pseudo
        m.pseudo_observables
    elseif class in [:shocks, :stdshocks]
        m.exogenous_shocks
    else
        throw(ArgumentError("Invalid class: $class. Must be :states, :obs, :pseudo, :shocks, or :stdshocks"))
    end

    out = Base.filter(key -> dict[key] == index, collect(keys(dict)))
    if length(out) == 0
        error("Key corresponding to index $index not found for class: $class")
    elseif length(out) > 1
        error("Multiple keys corresponding to index $index found for class: $class")
    else
        return out[1]
    end
end

# Parse population mnemonic into 2 Nullable{Symbol}s from one
function parse_population_mnemonic(m::AbstractModel)
    mnemonic = get_setting(m, :population_mnemonic)
    if isnull(mnemonic)
        return [Nullable{Symbol}(), Nullable{Symbol}()]
    else
        return map(s -> Nullable(Symbol(s)), split(string(get(mnemonic)), DSGE_DATASERIES_DELIM))
    end
end

# From an augmented state space with anticipated policy shocks, get indices
# corresponding to pre-ZLB states, shocks, and observables
function inds_states_no_ant(m::AbstractModel)
    if n_anticipated_shocks(m) > 0
        ind_ant1 = m.endogenous_states[:rm_tl1]
        ind_antn = m.endogenous_states[Symbol("rm_tl$(n_anticipated_shocks(m))")]
        return [1:(ind_ant1-1); (ind_antn+1):n_states_augmented(m)]
    else
        return collect(1:n_states_augmented(m))
    end
end

function inds_shocks_no_ant(m::AbstractModel)
    if n_anticipated_shocks(m) > 0
        ind_ant1 = m.exogenous_shocks[:rm_shl1]
        ind_antn = m.exogenous_shocks[Symbol("rm_shl$(n_anticipated_shocks(m))")]
        return [1:(ind_ant1-1); (ind_antn+1):n_shocks_exogenous(m)]
    else
        return collect(1:n_shocks_exogenous(m))
    end
end

function inds_obs_no_ant(m::AbstractModel)
    if n_anticipated_shocks(m) > 0
        ind_ant1 = m.observables[:obs_nominalrate1]
        ind_antn = m.observables[Symbol("obs_nominalrate$(n_anticipated_shocks(m))")]
        return [1:(ind_ant1-1); (ind_antn+1):n_observables(m)]
    else
        return collect(1:n_observables(m))
    end
end

# Interface for I/O settings
spec(m::AbstractModel)         = m.spec
subspec(m::AbstractModel)      = m.subspec
saveroot(m::AbstractModel)     = get_setting(m, :saveroot)
dataroot(m::AbstractModel)     = get_setting(m, :dataroot)

# Interface for data
data_vintage(m::AbstractModel)    = get_setting(m, :data_vintage)
data_id(m::AbstractModel)         = get_setting(m, :data_id)
cond_vintage(m::AbstractModel)    = get_setting(m, :cond_vintage)
cond_id(m::AbstractModel)         = get_setting(m, :cond_id)
cond_full_names(m::AbstractModel) = get_setting(m, :cond_full_names)
cond_semi_names(m::AbstractModel) = get_setting(m, :cond_semi_names)
use_population_forecast(m::AbstractModel) = get_setting(m, :use_population_forecast)
hpfilter_population(m::AbstractModel)     = get_setting(m, :hpfilter_population)

# Interface for general computation settings
use_parallel_workers(m::AbstractModel)    = get_setting(m, :use_parallel_workers)

# Interface for estimation settings
reoptimize(m::AbstractModel)          = get_setting(m, :reoptimize)
calculate_hessian(m::AbstractModel) = get_setting(m, :calculate_hessian)
hessian_path(m::AbstractModel)      = get_setting(m, :hessian_path)
n_hessian_test_params(m::AbstractModel) = get_setting(m, :n_hessian_test_params)

# Interface for Metropolis-Hastings settings
n_mh_blocks(m::AbstractModel)      =  get_setting(m, :n_mh_blocks)
n_mh_simulations(m::AbstractModel) =  get_setting(m, :n_mh_simulations)
n_mh_burn(m::AbstractModel)        =  get_setting(m, :n_mh_burn)
mh_thin(m::AbstractModel)          =  get_setting(m, :mh_thin)

# Interface for forecast settings
date_forecast_start(m::AbstractModel)   = get_setting(m, :date_forecast_start)
forecast_block_size(m::AbstractModel)   = get_setting(m, :forecast_block_size)
forecast_start_block(m::AbstractModel)  = get_setting(m, :forecast_start_block)
forecast_input_file_overrides(m::AbstractModel) = get_setting(m, :forecast_input_file_overrides)
forecast_uncertainty_override(m::AbstractModel) = get_setting(m, :forecast_uncertainty_override)
forecast_smoother(m::AbstractModel)     = get_setting(m, :forecast_smoother)
forecast_tdist_df_val(m::AbstractModel) = get_setting(m, :forecast_tdist_df_val)
forecast_tdist_shocks(m::AbstractModel) = get_setting(m, :forecast_tdist_shocks)
forecast_zlb_value(m::AbstractModel)    = get_setting(m, :forecast_zlb_value)
impulse_response_horizons(m::AbstractModel) = get_setting(m, :impulse_response_horizons)
n_shockdec_periods(m::AbstractModel)    = index_shockdec_end(m) - index_shockdec_start(m) + 1

# Interface for alternative policy settings
alternative_policy(m::AbstractModel) = get_setting(m, :alternative_policy)

function date_forecast_end(m::AbstractModel)
    date = date_forecast_start(m) + Dates.Month(3 * (forecast_horizons(m)-1))
    return Dates.lastdayofquarter(date)
end

function forecast_horizons(m::AbstractModel; cond_type::Symbol = :none)
    horizons = get_setting(m, :forecast_horizons)
    if cond_type == :none
        return horizons
    else
        return horizons - n_conditional_periods(m)
    end
end

function date_shockdec_start(m::AbstractModel)
    startdate = get_setting(m, :shockdec_startdate)
    if !isnull(startdate)
        return get(startdate)
    else
        return date_mainsample_start(m)
    end
end

function date_shockdec_end(m::AbstractModel)
    enddate =  get_setting(m, :shockdec_enddate)
    if !isnull(enddate)
        return get(enddate)
    else
        return date_forecast_end(m)
    end
end

"""
```
load_parameters_from_file(m::AbstractModel,path::String)
```
Returns a vector of parameters, read from a file, suitable for updating `m`.
"""
function load_parameters_from_file(m::AbstractModel, path::String)

    if isfile(path) && splitext(path)[2] == ".h5"
        x  = h5open(path, "r") do file
            try
                read(file, "params")
            catch
                error("$path does not contain variable params")
            end
        end
    else
        error("$path is not a valid HDF5 file.")
    end

    @assert length(x) == length(m.parameters)
    @assert eltype(x) == typeof(m.parameters[1].value)
    return x
end

"""
```
specify_mode!(m::AbstractModel, mode_file::String=""; verbose=:low)
```

Updates the values of `m.parameters` with the values from
`mode_file`. Sets `reoptimize` setting to `false`.

Usage: should be run before calling `estimate(m)`, e.g.:

    m = Model990()
    specify_mode!(m, modefile)
    estimate(m)
"""
function specify_mode!(m::AbstractModel, mode_file::String = ""; verbose=:low)

    m <= Setting(:reoptimize, false)

    if mode_file == ""
        mode_file = inpath(m, "user", "paramsmode.h5")
    end
    mode_file = normpath(mode_file)

    update!(m,load_parameters_from_file(m,mode_file))

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        println("Loaded previous mode from $mode_file.")
    end

end

"""
```
specify_hessian(m::AbstractModel, path::String=""; verbose=:low)
```

Specify a Hessian matrix calculated at the posterior mode to use in the model estimation. If
no path is provided, will attempt to detect location.
"""
function specify_hessian(m::AbstractModel, path::String=""; verbose=:low)
    if isempty(path)
        path = inpath(m, "user", "hessian.h5")
    end

    if isfile(path) && splitext(path)[2] == ".h5"
        m <= Setting(:hessian_path, normpath(abspath(path)))
    else
        error("Invalid input Hessian file: $path",)
    end

    m <= Setting(:calculate_hessian, false)

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        println("Specified hessian from $path.")
    end
end

#=
Build paths to where input/output/results data are stored.

Description:
Creates the proper directory structure for input and output files, treating the DSGE/save
    directory as the root of a savepath directory subtree. Specifically, the following
    structure is implemented:

    dataroot/

    savepathroot/
                 output_data/<spec>/<subspec>/log/
                 output_data/<spec>/<subspec>/<out_type>/raw/
                 output_data/<spec>/<subspec>/<out_type>/work/
                 output_data/<spec>/<subspec>/<out_type>/tables/
                 output_data/<spec>/<subspec>/<out_type>/figures/

Note: we refer to the savepathroot/output_data/<spec>/<subspec>/ directory as saveroot.
=#
# """
# ```
# logpath(model)
# ```
# Returns path to log file. Path built as
# ```
# <output root>/output_data/<spec>/<subspec>/log/log_<filestring>.log
# ```
# """
# function logpath(m::AbstractModel)
#     return savepath(m, "log", "log.log")
# end

strs = [:work, :raw, :tables, :figures, :log]
fns = [Symbol(x, "path") for x in strs]
for (str, fn) in zip(strs, fns)
    @eval begin
        # First eval function
        function $fn(m::AbstractModel,
                     out_type::String,
                     file_name::String = "",
                     filestring_addl::Vector{String}=Vector{String}())
            return savepath(m, out_type, $(string(str)), file_name, filestring_addl)
        end

        # Then, add docstring to it
        @doc $(
        """
        ```
        $fn(m::AbstractModel, out_type::String, file_name::String="", filestring_addl::Vector{String}=Vector{String}())
        ```

        Returns path to specific $str output file, creating containing directory as needed. If
        `file_name` not specified, creates and returns path to containing directory only. Path built
        as
        ```
        <output root>/output_data/<spec>/<subspec>/<out_type>/$str/<file_name>_<filestring>.<ext>
        ```
        """
        ) $fn
    end
end

# Not exposed to user. Actually create path and insert model string to file name.
function savepath(m::AbstractModel,
                  out_type::String,
                  sub_type::String,
                  file_name::String = "",
                  filestring_addl::Vector{String} = Vector{String}())
    # Containing directory
    dir = String(joinpath(saveroot(m), "output_data", spec(m), subspec(m), out_type, sub_type))

    if !isempty(file_name)
        base = filestring_base(m)
        return savepath(dir, file_name, base, filestring_addl)
    else
        return dir
    end
end

function savepath(dir::String,
                  file_name::String = "",
                  filestring_base::Vector{String} = Vector{String}(),
                  filestring_addl::Vector{String} = Vector{String}())
    if !isdir(dir)
        mkpath(dir)
    end

    if !isempty(file_name)
        (base, ext) = splitext(file_name)
        myfilestring = filestring(filestring_base, filestring_addl)
        file_name_detail = base * myfilestring * ext

        return joinpath(dir, file_name_detail)
    else
        return dir
    end
end

# Input data handled slightly differently, because it is not model-specific.
"""
```
inpath{T<:String}(m::AbstractModel, in_type::T, file_name::T="")
```

Returns path to specific input data file, creating containing directory as needed. If
`file_name` not specified, creates and returns path to containing directory only. Valid
`in_type` includes:

* `\"raw\"`: raw input series
* `\"data\"`: transformed data in model units
* `\"cond\"`: conditional data - nowcasts for the current forecast quarter, or related
* `\"user\"`: user-supplied data for starting parameter vector, hessian, or related
* `\"scenarios\"`: alternative scenarios

Path built as
```
<data root>/<in_type>/<file_name>
```
"""
function inpath(m::AbstractModel, in_type::String, file_name::String="")
    path = dataroot(m)
    # Normal cases.
    if in_type in ["raw", "data", "cond", "scenarios"]
        path = joinpath(path, in_type)
    # User-provided inputs. May treat this differently in the future.
    elseif in_type == "user"
        path = joinpath(path, "user")
    else
        error("Invalid in_type: ", in_type)
    end

    # Containing dir
    if !isdir(path)
        mkpath(path)
    end

    # If file_name provided, return full path
    if !isempty(file_name)
        path = joinpath(path, file_name)
    end

    return path
end

function filestring_base(m::AbstractModel)
    if !m.testing
        base = Vector{String}()
        for (skey, sval) in m.settings
            if sval.print
                push!(base, to_filestring(sval))
            end
        end
        return base
    else
        return ["test"]
    end
end

filestring(m::AbstractModel) = filestring(m, Vector{String}())
filestring(m::AbstractModel, d::String) = filestring(m, [String(d)])
function filestring(m::AbstractModel, d::Vector{String})
    base = filestring_base(m)
    return filestring(base, d)
end

function filestring(base::Vector{String}, d::Vector{String})
    filestrings = vcat(base, d)
    sort!(filestrings)
    return "_" * join(filestrings, "_")
end

function filestring(d::Vector{String})
    sort!(d)
    return "_" * join(d, "_")
end

"""
```
transform_to_model_space!{T<:AbstractFloat}(m::AbstractModel, values::Vector{T})
```

Transforms `values` from the real line to the model space, and assigns `values[i]` to
`m.parameters[i].value` for non-steady-state parameters. Recomputes the steady-state
paramter values.

### Arguments
- `m`: the model object
- `values`: the new values to assign to non-steady-state parameters.
"""
function transform_to_model_space!{T<:AbstractFloat}(m::AbstractModel, values::Vector{T})
    new_values = transform_to_model_space(m.parameters, values)
    update!(m, new_values)
    steadystate!(m)
end

"""
```
update!{T<:AbstractFloat}(m::AbstractModel, values::Vector{T})
```

Update `m.parameters` with `values`, recomputing the steady-state parameter values.

### Arguments:
- `m`: the model object
- `values`: the new values to assign to non-steady-state parameters.
"""
function update!{T<:AbstractFloat}(m::AbstractModel, values::Vector{T})
    update!(m.parameters, values)
    steadystate!(m)
end

"""
```
rand(d::DegenerateMvNormal, m::AbstractModel; cc::AbstractFloat = 1.0)
```

Generate a draw from `d` with variance optionally scaled by `cc^2`.
"""
function rand(d::DegenerateMvNormal, m::AbstractModel; cc::AbstractFloat = 1.0)
    return d.μ + cc*d.σ*randn(m.rng, length(d))
end

"""
`rand_prior(m::AbstractModel; ndraws::Int = 100_000)`

Draw a random sample from the model's prior distribution.
"""
function rand_prior(m::AbstractModel; ndraws::Int = 100_000)
    T = typeof(m.parameters[1].value)
    npara = length(m.parameters)
    priorsim = Array{T}(ndraws, npara)

    for i in 1:ndraws
        priodraw = Array{T}(npara)

        # Parameter draws per particle
        for j in 1:length(m.parameters)

            priodraw[j] = if !m.parameters[j].fixed
                prio = rand(m.parameters[j].prior.value)

                # Resample until all prior draws are within the value bounds
                while !(m.parameters[j].valuebounds[1] < prio < m.parameters[j].valuebounds[2])
                    prio = rand(m.parameters[j].prior.value)
                end

                prio
            else
                m.parameters[j].value
            end
        end
        priorsim[i,:] = priodraw'
    end

    priorsim
end

"""
```
type ShockGroup
```

The `ShockGroup` type is used in `prepare_means_table_shockdec` and
`plot_shock_decompositions`. When plotting shock decompositions, we usually want
to group the shocks into categories (financial, monetary policy, etc.) so that
the resulting grouped bar plot is legible.

### Fields

- `name::String`
- `shocks::Vector{Symbol}`
- `color::Colorant`
"""
type ShockGroup
    name::String
    shocks::Vector{Symbol}
    color::Colorant
end

function ShockGroup(name::String, shocks::Vector{Symbol}, color_name::Symbol)
    color = parse(Colorant, color_name)
    return ShockGroup(name, shocks, color)
end


text:
using JSONL
using Test, DataFrames, RDatasets

full_web = readfile("testfiles/jsonlwebsite.jsonl") |> DataFrame;
nrow_fw = nrow(full_web)

mtcars = dataset("datasets", "mtcars")
full_mtcars = readfile("testfiles/mtcars.jsonl") |> DataFrame;
# Fix R export differences
rename!(full_mtcars, :_row => :model);
rename!(full_mtcars, names(full_mtcars) .=> lowercase.(names(full_mtcars)))
rename!(mtcars, names(mtcars) .=> lowercase.(names(mtcars)))

# Read without promotion
noprom_mtcars = readfile("testfiles/mtcars.jsonl") |> DataFrame;
nrow_mt = nrow(mtcars)

oneline = readfile("testfiles/oneline.jsonl")|> DataFrame;
oneline_plus = readfile("testfiles/oneline_plus.jsonl")|> DataFrame;

escaped = readfile("testfiles/escapedeol.jsonl")|> DataFrame;

@testset "Read" begin
    @test full_web.name == ["Gilbert", "Alexa", "May", "Deloise"]
    @test full_web.wins[1] == [["straight", "7♣"], ["one pair", "10♥"]]
    @test full_web.wins[end] == [["three of a kind", "5♣"]]
    @test full_mtcars.mpg == mtcars.mpg
    @test noprom_mtcars.cyl[32] == 4
    @test noprom_mtcars.wt[30] == 2.77
    @test noprom_mtcars.qsec[16] == 17.82
    @test noprom_mtcars[!, :_row][16] == "Lincoln Continental"
    @test noprom_mtcars[!, :drat] == mtcars[!, :drat]
    @test oneline.name == ["Daniel"]
    @test oneline_plus.name == ["Daniel"]
    @test nrow(escaped) == 4
    @test escaped.name[1] == "Daniel\n"
    @test escaped.age[2] == "}"
end

@testset "Mmap Full File" begin
# full file equal
    @test readfile("testfiles/jsonlwebsite.jsonl", usemmap = true) |> DataFrame == full_web
    @test readfile("testfiles/mtcars.jsonl", usemmap = true) |> DataFrame == noprom_mtcars
    @test readfile("testfiles/oneline.jsonl", usemmap = true) |> DataFrame== oneline
    @test readfile("testfiles/oneline_plus.jsonl", usemmap = true)  |> DataFrame == oneline_plus
    @test readfile("testfiles/escapedeol.jsonl", usemmap = true)  |> DataFrame == escaped
end

@testset "skip & nrows" begin
# skip + nrows = nrow(file)
    @test readfile("testfiles/jsonlwebsite.jsonl", skip = 1, nrows = nrow_fw-1)  |> DataFrame == full_web[2:end, :]
    @test readfile("testfiles/mtcars.jsonl", skip = 2, nrows = nrow_mt-2)  |> DataFrame == noprom_mtcars[3:end, :]

# skip + nrows < nrow(file)
    @test readfile("testfiles/jsonlwebsite.jsonl", skip = 1, nrows = 2)  |> DataFrame == full_web[2:3, :]
    @test readfile("testfiles/mtcars.jsonl", skip = 2, nrows = 2)  |> DataFrame == noprom_mtcars[3:4, :]
    @test readfile("testfiles/escapedeol.jsonl", skip = 2, nrows = 1)  |> DataFrame == escaped[3:3, :]

# skip + nrows > nrow(file) (through nrow)
    @test readfile("testfiles/jsonlwebsite.jsonl", skip = 1, nrows = nrow_fw)  |> DataFrame == full_web[2:end, :]
    @test readfile("testfiles/mtcars.jsonl", skip = 12, nrows = nrow_mt + 10)  |> DataFrame == noprom_mtcars[13:end, :]
    @test readfile("testfiles/oneline.jsonl", skip = 0, nrows = 5)  |> DataFrame == oneline
    @test readfile("testfiles/oneline_plus.jsonl", skip = 0, nrows = 2)  |> DataFrame == oneline_plus
    @test readfile("testfiles/escapedeol.jsonl", skip = 2, nrows = 10)  |> DataFrame == escaped[3:end, :]

# skip + nrows > nrow(file) (through skip)
    @test readfile("testfiles/jsonlwebsite.jsonl", skip = nrow_fw+1, nrows = 1)  |> DataFrame == DataFrame()
    @test readfile("testfiles/mtcars.jsonl", skip = nrow_mt +12, nrows = 120)  |> DataFrame ==  DataFrame()
    @test readfile("testfiles/oneline.jsonl", skip = 2, nrows = 10)  |> DataFrame == DataFrame()
    @test readfile("testfiles/oneline_plus.jsonl", skip = 2, nrows = 123)  |> DataFrame == DataFrame()
    @test readfile("testfiles/escapedeol.jsonl", skip = 5, nrows = 1)  |> DataFrame == DataFrame()

# skip = nrow(file) + nrows > 0
    @test readfile("testfiles/jsonlwebsite.jsonl", skip = nrow_fw, nrows = 10)  |> DataFrame == DataFrame()
    @test readfile("testfiles/mtcars.jsonl", skip = nrow_mt, nrows = 1)  |> DataFrame == DataFrame()
    @test readfile("testfiles/oneline.jsonl", skip = 1, nrows = 12) |> DataFrame == DataFrame()
    @test readfile("testfiles/oneline_plus.jsonl", skip = 1, nrows = 1) |> DataFrame == DataFrame()
    @test readfile("testfiles/escapedeol.jsonl", skip = 4, nrows = 1) |> DataFrame == DataFrame()
end

@testset "skip" begin
# skip = nrow(file)
    @test readfile("testfiles/jsonlwebsite.jsonl", skip = nrow_fw)  |> DataFrame == DataFrame()
    @test readfile("testfiles/mtcars.jsonl", skip = nrow_mt)  |> DataFrame == DataFrame()
    @test readfile("testfiles/oneline.jsonl", skip = 1)  |> DataFrame == DataFrame()
    @test readfile("testfiles/oneline_plus.jsonl", skip = 1)  |> DataFrame == DataFrame()
    @test readfile("testfiles/escapedeol.jsonl", skip = 4)  |> DataFrame == DataFrame()

# skip > nrow(file)
    @test readfile("testfiles/jsonlwebsite.jsonl", skip = nrow_fw + 1)  |> DataFrame == DataFrame()
    @test readfile("testfiles/mtcars.jsonl", skip = nrow_mt + 42) |> DataFrame  == DataFrame()
    @test readfile("testfiles/mtcars.jsonl", skip = typemax(Int))  |> DataFrame == DataFrame()
    @test readfile("testfiles/oneline.jsonl", skip = 2)  |> DataFrame == DataFrame()
    @test readfile("testfiles/oneline_plus.jsonl", skip = 2)  |> DataFrame == DataFrame()
    @test readfile("testfiles/escapedeol.jsonl", skip = 5) |> DataFrame  == DataFrame()

# skip < nrow(file)
    @test readfile("testfiles/jsonlwebsite.jsonl", skip = nrow_fw - 1)  |> DataFrame == full_web[end:end, :]
    @test readfile("testfiles/mtcars.jsonl", skip = nrow_mt - 12)  |> DataFrame == noprom_mtcars[(end-11):end, :]
    @test readfile("testfiles/escapedeol.jsonl", skip = 2)  |> DataFrame == escaped[3:end, :]
end

@testset "nrows" begin
# nrows < nrow(file)
    @test readfile("testfiles/jsonlwebsite.jsonl", nrows = 2)  |> DataFrame == full_web[begin:2, :]
    @test readfile("testfiles/mtcars.jsonl", nrows = 12)  |> DataFrame == noprom_mtcars[begin:12, :]
    @test readfile("testfiles/escapedeol.jsonl", nrows = 3)  |> DataFrame == escaped[begin:3, :]

# nrows = nrow(file)
    @test readfile("testfiles/jsonlwebsite.jsonl", nrows = nrow_fw)  |> DataFrame == full_web
    @test readfile("testfiles/mtcars.jsonl", nrows = nrow_mt)  |> DataFrame == noprom_mtcars
    @test readfile("testfiles/oneline.jsonl", nrows = 1)  |> DataFrame == oneline
    @test readfile("testfiles/oneline_plus.jsonl", nrows = 1)  |> DataFrame == oneline_plus
    @test readfile("testfiles/escapedeol.jsonl", nrows = 4)  |> DataFrame == escaped

# nrows > nrow(file)
    @test readfile("testfiles/jsonlwebsite.jsonl", nrows = nrow_fw+1)  |> DataFrame == full_web
    @test readfile("testfiles/mtcars.jsonl", nrows = nrow_mt+100)  |> DataFrame == noprom_mtcars
    @test readfile("testfiles/oneline.jsonl", nrows = 2)  |> DataFrame == oneline
    @test readfile("testfiles/oneline_plus.jsonl", nrows = 2)  |> DataFrame == oneline_plus
    @test readfile("testfiles/escapedeol.jsonl", nrows = 5)  |> DataFrame == escaped
end

writefile("full_web.jsonl", full_web)
writefile("full_mtcars.jsonl", full_mtcars)
writefile("oneline2.jsonl", oneline)
writefile("oneline_plus2.jsonl", oneline_plus)
writefile("escaped2.jsonl", escaped)
@testset "write" begin
    @test readfile("full_web.jsonl") |> DataFrame == full_web
    @test readfile("full_mtcars.jsonl") |> DataFrame == full_mtcars
    @test readfile("oneline2.jsonl") |> DataFrame == oneline
    @test readfile("oneline_plus2.jsonl") |> DataFrame == oneline_plus
    @test readfile("escaped2.jsonl") |> DataFrame == escaped
end


text:
module contactForceLaw_2_boxes_Prismatic

using Modia3D
using Modia3D.StaticArrays
import Modia3D.ModiaMath

vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)    # material of SolidFileMesh
vmat2 = deepcopy(vmat1)                                           # material of convex decomposition of SolidFileMesh
vmat2.transparency = 0.7

cmat = "Steel"

@assembly TwoBoxes begin
  world = Modia3D.Object3D(visualizeFrame=true)

  boxMoving     = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBox(0.5,0.5,0.5) , "Steel", vmat1; contactMaterial = cmat) )
  helpFrame1 = Modia3D.Object3D(visualizeFrame=false)
  helpFrame2 = Modia3D.Object3D(visualizeFrame=false)
  prisX = Modia3D.Prismatic(world, helpFrame1; axis=1, v_start=-6.0, s_start=0.0, canCollide=true)
  prisY = Modia3D.Prismatic(helpFrame1, helpFrame2, axis=2, v_start=2.0, s_start=0.0, canCollide=true)
  prisZ = Modia3D.Prismatic(helpFrame2, boxMoving, axis=3, v_start=4.0, s_start=0.0, canCollide=true)

  box     = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBox(3.0,2.0,5.0) , "Steel", vmat1; contactMaterial = cmat); r=[-3.0, 0.0, 1.5], fixed=true) # R=ModiaMath.rot2(-pi/3),
end

gravField = Modia3D.UniformGravityField(g=9.81, n=[-1,0,0])
threeD = TwoBoxes(sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=true, defaultFrameLength=0.7,nz_max = 100, enableContactDetection=true, visualizeContactPoints=true, visualizeSupportPoints=true))

# Modia3D.visualizeAssembly!( threeD )

model = Modia3D.SimulationModel( threeD )
ModiaMath.print_ModelVariables(model)
result = ModiaMath.simulate!(model; stopTime=5.0, tolerance=1e-6,interval=0.001, log=false)

ModiaMath.plot(result, ["prisX.r","prisX.v"])

println("... success of contactForceLaw_2_boxes_Prismatic.jl!")

end


text:
# ---
# title: 876. Middle of the Linked List
# id: problem876
# author: zhwang
# date: 2022-01-23
# difficulty: Easy
# categories: Linked List
# link: <https://leetcode.com/problems/middle-of-the-linked-list/description/>
# hidden: true
# ---
#
# Given a non-empty, singly linked list with head node `head`, return a middle
# node of linked list.
#
# If there are two middle nodes, return the second middle node.
#
#
#
# **Example 1:**
#
#
#
#     Input: [1,2,3,4,5]
#     Output: Node 3 from this list (Serialization: [3,4,5])
#     The returned node has value 3.  (The judge's serialization of this node is [3,4,5]).
#     Note that we returned a ListNode object ans, such that:
#     ans.val = 3, ans.next.val = 4, ans.next.next.val = 5, and ans.next.next.next = NULL.
#
#
# **Example 2:**
#
#
#
#     Input: [1,2,3,4,5,6]
#     Output: Node 4 from this list (Serialization: [4,5,6])
#     Since the list has two middle nodes with values 3 and 4, we return the second one.
#
#
#
#
# **Note:**
#
#   * The number of nodes in the given list will be between `1` and `100`.
#
#
## @lc code=start
using LeetCode

## double pointers
function middle_node_by_double_pointers(head::ListNode)::ListNode
    p1 = p2 = head
    while !isnothing(p2)
        p2 = next(p2)
        isnothing(p2) && return p1
        p2, p1 = p2.next, p1.next
    end
    return p1
end

## @lc code=end


text:
function give_rand_p(p_range,p_fixed=nothing,indices=nothing)
    if p_fixed == nothing
        p = [(p_range[j][2] -p_range[j][1])*rand() + p_range[j][1] for j in 1:length(p_range)]
    else
        p =  zeros(length(p_range))
        j = 1
        for i in 1:length(p_range)
            if i in indices
                p[i] = p_fixed[j]
                j += 1
            else
                p[i] = (p_range[i][2] -p_range[i][1])*rand() + p_range[i][1]
            end
        end
    end
    p
end

function calc_mean_var(f,p_range,N)
    y1 = Array(f(give_rand_p(p_range)))
    y0 = zero(y1)
    v = zero(y1)
    for i in 1:N
        y1 = Array(f(give_rand_p(p_range)))
        @. y0 += y1
        @. v += y1^2
    end
    y0 = @. y0/N
    y0_sq = [i.^2 for i in y0]
    v = @. v/N - y0_sq
    y0,v
end

function first_order_var(f,p_range,N,y0)
    ys = Array{typeof(y0)}(undef,length(p_range))
    for i in 1:length(p_range)
        y = zero(y0)
        for j in 1:N
            p2 = give_rand_p(p_range)
            p1 = give_rand_p(p_range,[p2[i]],[i])
            yer =  Array(f(p1)) .* Array(f(p2))
            @. y += yer
        end
        y = @. y/N - y0^2
        ys[i] = copy(y)
    end
    ys
end

function second_order_var(f,p_range,N,y0)
    ys = Array{typeof(y0)}(undef,Int((length(p_range)*(length(p_range)-1))/2))
    curr = 1
    for i in 1:length(p_range)
        for j in i+1:length(p_range)
            y = zero(y0)
            for k in 1:N
                p2 = give_rand_p(p_range)
                p1 = give_rand_p(p_range,[p2[i],p2[j]],[i,j])
                y .+=  Array(f(p1)) .* Array(f(p2))
            end
            y = @. y/N - y0^2
            ys[curr] = copy(y)
            curr += 1
        end
    end
    ys_frst_order = first_order_var(f,p_range,N,y0)
    j = 1
    for i in 1:length(p_range)
        for k in i+1:length(p_range)
            ys[j] = @. ys[j] - ( ys_frst_order[i] + ys_frst_order[k] )
            j += 1
        end
    end
    ys
end

function total_var(f,p_range,N,y0)
    ys = Array{typeof(y0)}(undef,length(p_range))
    for i in 1:length(p_range)
        y = zero(y0)
        for j in 1:N
            p_fixed_all = []
            p_fixed_indices = []
            p2 = give_rand_p(p_range)
            for j in 1:length(p2)
                if j != i
                    push!(p_fixed_all,p2[j])
                    push!(p_fixed_indices,j)
                end
            end
            p1 = give_rand_p(p_range,p_fixed_all,p_fixed_indices)
            yer =  Array(f(p1)) .* Array(f(p2))
            @. y += yer
        end
        y = @. y/N - y0^2
        ys[i] = copy(y)
    end
    ys
end

function sobol_sensitivity(f,p_range,N,order=2)
    y0,v = calc_mean_var(f,p_range,N)
    if order == 1
        first_order = first_order_var(f,p_range,N,y0)
        for i in 1:length(first_order)
            first_order[i] = @. first_order[i] / v
        end
        first_order
    elseif order == 2
        second_order = second_order_var(f,p_range,N,y0)
        for i in 1:length(second_order)
            second_order[i] = @. second_order[i] / v
        end
        second_order
    else
        total_indices = total_var(f,p_range,N,y0)
        for i in 1:length(total_indices)
            total_indices[i] = @. 1 - (total_indices[i] / v)
        end
        total_indices
    end
end

function sobol_sensitivity(prob::DiffEqBase.DEProblem,alg,t,p_range,N,order=2)
    f = function (p)
        prob1 = remake(prob;p=p)
        Array(solve(prob1,alg;saveat=t))
    end
    @assert length(prob.p) == length(p_range)
    sobol_sensitivity(f,p_range,N,order)
end


max_stars_repo_name: sugitanishi/competitive-programming

max_issues_repo_name: sugitanishi/competitive-programming

max_forks_repo_name: sugitanishi/competitive-programming


text:
################################################################################
#
#    Module Potentials.jl
#
#    This module will implement a variety of interatomic potentials and defines abstract types
#    that allow these potentials to be used in other packages (PotentialLearning.jl, PotentialUQ.jl,
#    Atomistic.jl, etc...).
#       'implementation' means:
#           1. Having a defined structure for each potential
#               1.1 The structure holds all of the necessary parameters for evaluating energies, forces, ...
#               1.2 The potential structure should expose the trainable and nontrainable parameters (necessary for learning).
#           2. Having a method to get the potential energy of given configuration, as
#               defined by that potential.
#           3. Having a method to produce the force of a given configuration, as defined
#               by that potential.
#           4. Having a method to produce the stresses of a given configuration, as defined
#               by that potential.
#           5. (For inference) Having a method to produce the gradient of each of the above methods with
#               respect to the potential parameters.
#

################################################################################
module InteratomicPotentials

using StaticArrays
using LAMMPS
using LinearAlgebra
using AtomsBase
using Unitful
using UnitfulAtomic
include("Utilities/utils.jl")
include("PotentialTypes/types.jl")
export NeighborList, neighborlist
export potential_energy, force, virial, virial_stress
export grad_potential_energy, grad_force, grad_virial, grad_virial_stress
export SNAPParams, compute_snap, get_num_coeffs
export EmpiricalPotential, LennardJones
# export BornMayer, Coulomb, GaN, MixedPotential

end


text:
"Create the pattern of the Cholesky factor based on the upper triangle of A"
function cholpattern{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti})
    m,n = size(A)
    m == n || error("A must be square")
    parent = etree(A)
    I = Ti[]
    J = Ti[]
    sizehint!(I, nnz(A))
    sizehint!(J, nnz(A))
    for j in Ti[1:n;]
        cj = Base.SparseMatrix.ereach(A, j, parent)
        append!(J,cj)
        push!(J,j)
        append!(I,fill(j,length(cj) + 1))
    end
    sparse(I,J,one(Tv))
end

"""
Convert sparse to dense if the proportion of nonzeros exceeds a threshold.
A no-op for other matrix types.
"""
function densify(S,threshold=0.3)
    issparse(S) || return S
    nnz(S)/(*(size(S)...)) > threshold || return S
    isbits(eltype(S)) && return full(S)
    nzs = nonzeros(S)
    nz1 = nzs[1]
    T = typeof(nz1)
    isa(nz1,Array) && isbits(eltype(nz1)) || error("Nonzeros must be a bitstype or an array of same")
    sz1 = size(nz1)
    all(x->typeof(x) == T && size(x) == sz1, nzs) || error("Inconsistent dimensions in array nonzeros")
    M,N = size(S)
    m,n = size(nz1,1),size(nz1,2) # this construction allows for nz1 to be a vector
    res = Array(eltype(nz1),M*m,N*n)
    rv = rowvals(S)
    for j in 1:size(S,2)
        jm1 = j - 1
        for k in nzrange(S,j)
            copy!(sub(res,(rv[k]-1)*m+(1:m),jm1*n+(1:n)),nzs[k])
        end
    end
    res
end

function Base.scale!(x::Number,t::UpperTriangular{Float64})
    m,n = size(t)
    td = t.data
    for j in 1:n, i in 1:j
        @inbounds td[i,j] *= x
    end
end

function Base.scale!{T<:Number}(s::T,t::LowerTriangular{T})
    m,n = size(t)
    td = t.data
    for j in 1:n, i in j:m
        @inbounds td[i,j] *= s
    end
    t
end


text:
module TraceUtilities

using ..Utilities, ..JSON

export TraceValue, TraceType, Trace, diff,
        ChartType, Visible, HoverInfo, Stream, GroupNorm, StackGaps, Mode,
        HoverOn, Shape, Line, Fill, MarkerSymbol, SizeMode,
        ThicknessMode, LenMode, ColorBar, GradientType, Marker, Selected,
        Unselected, TextPosition, ErrorType, Error, ConstrainText, BoxPoints,
        BoxMean, CoordinateType, ZSmooth, HistFunc, HistNorm, Direction,
        CurrentBin, Cumulative, Bins, Attributes

abstract type TraceValue <: ComponentValue end
abstract type TraceType <: ComponentType end
abstract type Trace <: Component end

function Utilities.diff(ref::ComponentType,obj::ComponentType,dict=Dict())
    names = propertynames(ref)
    for name in names
        refval = getproperty(ref,name)
        objval = getproperty(obj,name)
        d = Utilities.diff(refval,objval)
        if d != nothing && !isempty(d)
            dict[name] = d
        end
    end
    if typeof(obj) <: Trace && obj.type.value != "scatter"
        dict[:type] = obj.type.value
    end
    if haskey(dict,:attributes)
        return merge(dict[:attributes],delete!(dict,:attributes))
    end
    return dict
end

JSON.lower(t::TraceValue) = t.value

charttypes = (:scatter,:bar,:box,:heatmap,:histogram,:histogram2d,:histogram2dcontour,:pie,:contour,:scatterternary,:violin,:scatter3d,:surface,:mesh3d,:cone,:streamtube,:scattergeo,:choropleth,:scattergl,:splom,:pointcloud,:heatmapgl,:parcoords,:parcats,:scattermapbox,:sankey,:table,:carpet,:scattercarpet,:contourcarpet,:ohlc,:candlestick,:scatterpolar,:scatterpolargl,:barpolar)
const ChartType = Enumerated{:charttype,charttypes}

const Visible = Enumerated{:Visible,(true,false,:legendonly)}

const HoverInfo = FlagList{(:x,:y,:z,:text,:name),(:all,:none,:skip)}

struct Stream <: TraceType
    token::String
    maxpoints::RangeValue{Int,(0,10000),(true,true)}

    Stream(token="",maxpts::Int=500) = new(token,RangeValue((0,10000),maxpts))
end

const GroupNorm = Enumerated{:GroupNorm,(:fraction,:percent,Symbol(""))}

const StackGaps = Enumerated{:StackGaps,(Symbol("infer zero"),:interpolate)}

const Mode = FlagList{(:lines,:markers,:text),(:none,)}

const HoverOn = FlagList{(:points,:fills),()}

const Shape = Enumerated{:Shape,(:linear,:spline,:hv,:vh,:hvh,:vhv)}

struct Line <: TraceType
    color::Color
    width::Size{Int}
    shape::Shape
    smoothing::RangeValue{Float64,(0.0,1.3),(true,true)}
    dash::Dash
    simplify::Bool
    cauto::Bool
    cmin::Float64
    cmax::Float64
    cmid::Float64
    colorscale::Union{ColorScale,ColorScaleSet}
    autocolorscale::Bool
    reversescale::Bool

    function Line(;
        color::String = "#444",
        width::Int = 2,
        shape::String = "linear",
        smoothing::Float64 = 1.0,
        dash::String = "solid",
        simplify::Bool = true,
        cauto::Bool = true,
        cmin::Float64 = 0.0,
        cmax::Float64 = 0.0,
        cmid::Float64 = 0.0,
        colorscale::Union{Vector{Tuple{Float64,String}},String} = "RdBu",
        autocolorscale::Bool = true,
        reversescale::Bool = false)
    return new(
        Color(color),
        Size{Int}(width),
        Shape(shape),
        RangeValue((0.0,1.3),smoothing),
        Dash(dash),
        simplify,
        cauto,
        cmin,
        cmax,
        cmid,
        typeof(colorscale) == String ? ColorScaleSet(colorscale) :
            ColorScale([(Normalized(Float64(n)),Color(s)) for (n,s) in colorscale]),
        autocolorscale,
        reversescale)
    end
end

const Fill = Enumerated{:Fill,(:none,:tozeroy,:tozerox,:tonexty,:tonextx,:toself)}

symbols = Symbol.(("circle","circle-open","circle-dot","circle-open-dot","square","square-open","square-dot","square-open-dot","diamond","diamond-open","diamond-dot","diamond-open-dot","cross","cross-open","cross-dot","cross-open-dot","x","x-open","x-dot","x-open-dot","triangle-up","triangle-up-open","triangle-up-dot","triangle-up-open-dot","triangle-down","triangle-down-open","triangle-down-dot","triangle-down-open-dot","triangle-left","triangle-left-open","triangle-left-dot","triangle-left-open-dot","triangle-right","triangle-right-open","triangle-right-dot","triangle-right-open-dot","triangle-ne","triangle-ne-open","triangle-ne-dot","triangle-ne-open-dot","triangle-se","triangle-se-open","triangle-se-dot","triangle-se-open-dot","triangle-sw","triangle-sw-open","triangle-sw-dot","triangle-sw-open-dot","triangle-nw","triangle-nw-open","triangle-nw-dot","triangle-nw-open-dot","pentagon","pentagon-open","pentagon-dot","pentagon-open-dot","hexagon","hexagon-open","hexagon-dot","hexagon-open-dot","hexagon2","hexagon2-open","hexagon2-dot","hexagon2-open-dot","octagon","octagon-open","octagon-dot","octagon-open-dot","star","star-open","star-dot","star-open-dot","hexagram","hexagram-open","hexagram-dot","hexagram-open-dot","star-triangle-up","star-triangle-up-open","star-triangle-up-dot","star-triangle-up-open-dot","star-triangle-down","star-triangle-down-open","star-triangle-down-dot","star-triangle-down-open-dot","star-square","star-square-open","star-square-dot","star-square-open-dot","star-diamond","star-diamond-open","star-diamond-dot","star-diamond-open-dot","diamond-tall","diamond-tall-open","diamond-tall-dot","diamond-tall-open-dot","diamond-wide","diamond-wide-open","diamond-wide-dot","diamond-wide-open-dot","hourglass","hourglass-open","bowtie","bowtie-open","circle-cross","circle-cross-open","circle-x","circle-x-open","square-cross","square-cross-open","square-x","square-x-open","diamond-cross","diamond-cross-open","diamond-x","diamond-x-open","cross-thin","cross-thin-open","x-thin","x-thin-open","asterisk","asterisk-open","hash","hash-open","hash-dot","hash-open-dot","y-up","y-up-open","y-down","y-down-open","y-left","y-left-open","y-right","y-right-open","line-ew","line-ew-open","line-ns","line-ns-open","line-ne","line-ne-open","line-nw","line-nw-open"))
const MarkerSymbol = Enumerated{:MarkerSymbol,symbols}

const SizeMode = Enumerated{:SizeMode,(:diameter,:area)}

const ThicknessMode = Enumerated{:ThicknessMode,(:fraction,:pixels)}

const LenMode = Enumerated{:LenMode,(:fraction,:pixels)}

struct ColorBar <: TraceType
    thicknessmode::ThicknessMode
    thickness::Size{Int}
    lenmode::LenMode
    len::Size{Int}
    x::RangeValue{Float64,(-2,3),(true,true)}
    xanchor::XAnchor
    xpad::Size{Int}
    y::RangeValue{Float64,(-2,3),(true,true)}
    yanchor::YAnchor
    ypad::Size{Int}
    outlinecolor::Color
    outlinewidth::Size{Int}
    bordercolor::Color
    borderwidth::Size{Int}
    bgcolor::Color
    # =================================================================
    # TODO: Make sure this gets written at the correct level from `diff`
    ticks::Ticks # See also LayoutUtilities.Axis
    # =================================================================
    title::Title

    function ColorBar(;
        thicknessmode::String = "pixels",
        thickness::Int = 30,
        lenmode::String = "fraction",
        len::Int = 1,
        x::Float64 = 1.02,
        xanchor::String = "left",
        xpad::Int = 10,
        y::Float64 = 0.5,
        yanchor::String = "middle",
        ypad::Int = 10,
        outlinecolor::String = "#444",
        outlinewidth::Int = 1,
        bordercolor::String = "#444",
        borderwidth::Int = 0,
        bgcolor::String = "rgba(0,0,0,0)",
        # Ticks child component. Should at the same level.
        tickmode::String = "auto",
        nticks::Int = 0,
        tick0::Any = "",
        dtick::Any = "",
        tickvals::Vector{Float64} = Float64[],
        ticktext::Vector{String} = String[],
        ticks::String = "inside",
        ticklen::Int = 5,
        tickwidth::Int = 1,
        tickcolor::String = "#444",
        showticklabels::Bool = true,
        tickfont::Font = Font(),
        tickangle::Float64 = 0.,
        tickprefix::String = "",
        showtickprefix::String = "all",
        ticksuffix::String = "",
        showticksuffix::String = "all",
        showexponent::String = "all",
        exponentformat::String = "B",
        separatethousands::Bool = false,
        tickformat::String = "",
        tickformatstops::TickFormatStops = TickFormatStops(),
        # ===================================================
        title::Union{String,Title} = "")
        return new(
            ThicknessMode(thicknessmode),
            Size{Int}(thickness),
            LenMode(lenmode),
            Size{Int}(len),
            RangeValue((-2,3),x),
            XAnchor(xanchor),
            Size{Int}(xpad),
            RangeValue((-2,3),y),
            YAnchor(yanchor),
            Size{Int}(ypad),
            Color(outlinecolor),
            Size{Int}(outlinewidth),
            Color(bordercolor),
            Size{Int}(borderwidth),
            Color(bgcolor),
            Ticks(
                tickmode=tickmode,
                nticks=nticks,
                tick0=tick0,
                dtick=dtick,
                tickvals=tickvals,
                ticktext=ticktext,
                ticks=ticks,
                ticklen=ticklen,
                tickwidth=tickwidth,
                tickcolor=tickcolor,
                showticklabels=showticklabels,
                tickfont=tickfont,
                tickangle=tickangle,
                tickprefix=tickprefix,
                showtickprefix=showtickprefix,
                ticksuffix=ticksuffix,
                showticksuffix=showticksuffix,
                exponentformat=exponentformat,
                showexponent=showexponent,
                separatethousands=separatethousands,
                tickformat=tickformat,
                tickformatstops=tickformatstops),
            typeof(title) == String ? Title(text=title) : title
        )
    end
end

const GradientType = Enumerated{:GradientType,(:radial,:horizontal,:vertical,:none)}

struct Gradient <: TraceType
    type::GradientType
    color::Color

    function Gradient(;
        type::String = "none",
        color::String = "#000"
        )
        return new(
            GradientType(type),
            Color(color))
    end
end

struct Marker <: TraceType
    symbol::Union{MarkerSymbol,Vector{MarkerSymbol}}
    opacity::Union{Float64,Vector{Float64}}
    size::Union{Size{Int},Vector{Size{Int}}}
    maxdisplayed::Size{Int}
    sizeref::Float64
    sizemin::Size
    sizemode::SizeMode
    colorbar::ColorBar
    line::Line
    gradient::Gradient
    color::Union{Color,Vector{Color}}
    cauto::Bool
    cmin::Float64
    cmax::Float64
    cmid::Float64
    colorscale::Union{ColorScale,ColorScaleSet}
    autocolorscale::Bool
    reversescale::Bool

    function Marker(;
        symbol::Union{String,Vector{String}} = "circle",
        opacity::Union{Float64,Vector{Float64}} = 1.0,
        size::Union{Int,Vector{Int}} = 6,
        maxdisplayed::Int = 0,
        sizeref = 1.,
        sizemin::Int = 0,
        sizemode::String = "diameter",
        colorbar::ColorBar = ColorBar(),
        line::Line = Line(width=1),
        gradient::Gradient = Gradient(),
        color::Union{String,Vector{String}} = "#444",
        cauto::Bool = true,
        cmin::Float64 = 0.0,
        cmax::Float64 = 0.0,
        cmid::Float64 = 0.0,
        colorscale::Union{Vector{Tuple{Float64,String}},String} = "RdBu",
        autocolorscale::Bool = true,
        reversescale::Bool = false)
        return new(
            MarkerSymbol.(symbol),
            opacity,
            Size{Int}.(size),
            Size{Int}(maxdisplayed),
            sizeref,
            Size{Int}(sizemin),
            SizeMode(sizemode),
            colorbar,
            line,
            gradient,
            Color.(color),
            cauto,
            cmin,
            cmax,
            cmid,
            typeof(colorscale) == String ? ColorScaleSet(colorscale) :
                ColorScale([(Normalized(Float64(n)),Color(s)) for (n,s) in colorscale]),
            autocolorscale,
            reversescale)
    end
end

struct Selected <: TraceType
    marker::Marker
    textfont::Font

    function Selected(;
        marker::Marker = Marker(),
        textfont::Font = Font())
        return new(
            marker,
            textfont)
    end
end

struct Unselected <: TraceType
    marker::Marker
    textfont::Font

    function Unselected(;
        marker::Marker = Marker(),
        textfont::Font = Font())
        return new(
            marker,
            textfont)
    end
end

const TextPosition = Enumerated{:TextPosition,Symbol.((
    "top left","top center","top right","middle left","middle center","middle right",
    "bottom left","bottom center","bottom right"))}

const ErrorType = Enumerated{:ErrorType,(:percent,:constant,:sqrt,:data)}

struct Error <: TraceType
    visible::Bool
    type::ErrorType
    symmetric::Bool
    array::Vector{Float64}
    arrayminus::Vector{Float64}
    value::Size{Float64}
    valueminus::Size{Float64}
    traceref::Size{Int}
    tracerefminus::Size{Int}
    copy_ystyle::Bool
    color::Color
    thickness::Size{Int}
    width::Size{Int}

    function Error(;
        visible::Bool = false,
        type::String = "constant",
        symmetric::Bool = true,
        array::Vector{Float64} = Float64[],
        arrayminus::Vector{Float64} = Float64[],
        value::Float64 = 10.,
        valueminus::Float64 = 10.,
        traceref::Int = 0,
        tracerefminus::Int = 0,
        copy_ystyle::Bool = false,
        color::String = "#000",
        thickness::Int = 2,
        width::Int = 1)
        return new(
            visible,
            ErrorType(type),
            symmetric,
            array,
            arrayminus,
            Size{Float64}(value),
            Size{Float64}(valueminus),
            Size{Int}(traceref),
            Size{Int}(tracerefminus),
            copy_ystyle,
            Color(color),
            Size{Int}(thickness),
            Size{Int}(width)
            )
    end
end

const ConstrainText = Enumerated{:ConstrainText,(:inside,:outside,:both,:none)}

const BoxPoints = Enumerated{:BoxPoints,(:all,:outliers,:suspectedoutliers,false)}

const BoxMean = Enumerated{:BoxMean,(true,:sd,false)}

const CoordinateType = Enumerated{:CoordinateType,(:array,:scaled)}

const ZSmooth = Enumerated{:ZSmooth,(:fast,:best,false)}

const HistFunc = Enumerated{:HistFunc,(:count,:sum,:avg,:min,:max)}

const HistNorm = Enumerated{:HistNorm,(:percent,:probability,:density,:probabilitydensity)}

const Direction = Enumerated{:Direction,(:increasing,:decreasing)}

const CurrentBin = Enumerated{:CurrentBin,(:include,:exclude,:half)}

struct Cumulative <: TraceType
    enabled::Bool
    direction::Direction
    currentbin::CurrentBin

    function Cumulative(;
        enabled::Bool = false,
        direction::String = "increasing",
        currentbin::String = "include")
        return new(
            enabled,
            Direction(direction),
            CurrentBin(currentbin))
    end
end

struct Bins <: TraceType
    start::Any
    stop::Any # plotly.js uses `end`
    size::Any

    function Bins(;
        start = "",
        stop = "",
        size = 0)
        return new(start,stop,size)
    end
end

mutable struct Attributes <: TraceType
    # ====================================
    # Attributes common across trace types
    # ====================================
    x::Vector{X} where X<:Union{String,Number}
    x0::Float64
    dx::Float64
    y::Vector{Y} where Y<:Union{String,Number}
    y0::Float64
    dy::Float64
    orientation::Orientation
    text::Union{String,Vector{String}}
    hovertext::Union{String,Vector{String}}
    hovertemplate::Union{String,Vector{String}}
    marker::Marker
    selected::Selected
    unselected::Unselected
    xcalendar::Calendar
    ycalendar::Calendar
    xaxis::String
    yaxis::String
    # ============================================
    # Attributes contained in the plotly.js object
    # ============================================
    customdata::Vector{Any}
    hoverinfo::HoverInfo
    hoverlabel::HoverLabel
    ids::Vector{String}
    legendgroup::String
    name::String
    opacity::Float64
    selectedpoints::Vector{Int}
    showlegend::Bool
    stream::Stream
    # transforms
    # type
    # uid
    # uirevision::Float64
    visible::Visible

    function Attributes(;
        # ====================================
        # Attributes common across trace types
        # ====================================
        x = Float64[],
        x0 = 0.0,
        dx = 1.0,
        y = Float64[],
        y0 = 0.0,
        dy = 1.0,
        orientation::String = "v",
        text::Union{String,Vector{String}} = "",
        hovertext::Union{String,Vector{String}} = "",
        hovertemplate::Union{String,Vector{String}} = "",
        marker::Marker = Marker(),
        selected::Selected = Selected(),
        unselected::Unselected = Unselected(),
        xcalendar::String = "gregorian",
        ycalendar::String = "gregorian",
        xaxis::String = "x",
        yaxis::String = "y",
        # ============================================
        # Attributes contained in the plotly.js object
        # ============================================
        customdata::Vector{Any} = [],
        hoverinfo::String = "all",
        hoverlabel::HoverLabel = HoverLabel(),
        ids::Vector{String} = String[],
        legendgroup::String = "",
        name::String = "",
        opacity::Float64 = 1.0,
        selectedpoints::Vector{Int} = Int[],
        showlegend::Bool = true,
        stream::Stream = Stream(),
        visible::Union{Bool,String,Symbol} = true)

        return new(
            # ====================================
            # Attributes common across trace types
            # ====================================
            collect(x),
            Float64(x0),
            Float64(dx),
            collect(y),
            Float64(y0),
            Float64(dy),
            Orientation(orientation),
            text,
            hovertext,
            hovertemplate,
            marker,
            selected,
            unselected,
            Calendar(xcalendar),
            Calendar(ycalendar),
            xaxis,
            yaxis,
            # ============================================
            # Attributes contained in the plotly.js object
            # ============================================
            customdata,
            HoverInfo(hoverinfo),
            hoverlabel,
            ids,
            legendgroup,
            name,
            opacity,
            selectedpoints,
            showlegend,
            stream,
            Visible(visible))
    end
end

end


text:
#=
Generate steady states and the stability state

Reference
---------
1. Source code from https://www.nature.com/articles/s41598-017-15889-2
=#

steadies1, jac_ms1, stab_modes1 = fixedpoint_gen(sp.ODE!, sp.u0, sp.p, [-10.:1.0:10.,-10.:1.0:10.])

steadies2, jac_ms2, stab_modes2 = fixedpoint_gen(sp.ODE!, sp.p, [-10.:1.0:10.,-10.:1.0:10.])

@test length(steadies2) == sum(length.([sp.ODE_features["StableStates"]
                                       ,sp.ODE_features["UnstableStates"]]))

@test sum(getfield.(stab_modes2, :stable)) == length(sp.ODE_features["StableStates"])

@test sum(getfield.(stab_modes2, :unstable)) == length(sp.ODE_features["UnstableStates"])


text:
# ---
# title: 1109. Corporate Flight Bookings
# id: problem1109
# author: Indigo
# date: 2021-08-31
# difficulty: Medium
# categories: Array, Math
# link: <https://leetcode.com/problems/corporate-flight-bookings/description/>
# hidden: true
# ---
#
# There are `n` flights, and they are labeled from `1` to `n`.
#
# We have a list of flight bookings.  The `i`-th booking `bookings[i] = [i, j,
# k]` means that we booked `k` seats from flights labeled `i` to `j` inclusive.
#
# Return an array `answer` of length `n`, representing the number of seats
# booked on each flight in order of their label.
#
#
#
# **Example 1:**
#
#
#
#     Input: bookings = [[1,2,10],[2,3,20],[2,5,25]], n = 5
#     Output: [10,55,45,25,25]
#
#
#
#
# **Constraints:**
#
#   * `1 <= bookings.length <= 20000`
#   * `1 <= bookings[i][0] <= bookings[i][1] <= n <= 20000`
#   * `1 <= bookings[i][2] <= 10000`
#
#
## @lc code=start
using LeetCode

function corp_flight_bookings(bookings::Vector{Vector{Int}}, n::Int)
    res = fill(0, n + 1)
    for b in bookings
        res[b[1]] += b[3]
        res[b[2] + 1] -= b[3]
    end
    cumsum!(res, res)
    @view(res[1:n])
end
## @lc code=end


text:
# MethodError (foo(::Int64, ::Int64) is ambiguous)

function foo(x::Int64, y)
    42
end

function foo(x, y::Int64)
    1729
end

function cant_know_which_one_at_this_point(x, y)
    foo(x, y)
end

cant_know_which_one_at_this_point(1, 2)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# XXX this is all kind of hacky

struct MLPOMDPSolver <: Solver
    solver
    updater::Nullable{Any}
end
function set_rng!(s::MLPOMDPSolver, rng::AbstractRNG)
    set_rng!(s.solver, rng)
    set_rng!(get(s.updater), rng)
end

struct MLPOMDPAgent <: Policy
    updater::Updater
    previous_belief::Nullable{Any}
    policy::Policy
    previous_action::Nullable{MLAction}
end

function solve(solver::MLPOMDPSolver, problem::MLMDP, up=get(solver.updater, nothing))
    internal_problem = MLPOMDP{MLState, MLAction, MLPhysicalState, typeof(problem.dmodel), typeof(problem.rmodel)}(problem.dmodel, problem.rmodel, problem.discount)
    policy = solve(solver.solver, internal_problem)
    if up == nothing
        up = POMDPs.updater(policy)
    else
        set_problem!(up, internal_problem)
    end
    return MLPOMDPAgent(up, nothing, policy, nothing)
end

function action(agent::MLPOMDPAgent, state::MLState)
    o = MLPhysicalState(state)
    if isnull(agent.previous_belief)
        belief = initialize_belief(agent.updater, ParticleGenerator(agent.policy.problem, state))
    else
        belief = update(agent.updater, get(agent.previous_belief), get(agent.previous_action), o)
    end
    a = action(agent.policy, belief)
    agent.previous_action = a
    agent.previous_belief = belief
    return a
end

struct ParticleGenerator
    physical::MLPhysicalState
    behaviors::BehaviorGenerator
end
function ParticleGenerator(problem::NoCrashProblem, state::Union{MLState, MLPhysicalState})
    return ParticleGenerator(MLPhysicalState(state),
                problem.dmodel.behaviors)
end

function rand(rng::AbstractRNG, gen::ParticleGenerator,
              full_s::MLState=MLState(gen.physical,
                                      Array{CarState}(length(gen.physical.cars))))
                                      #Array(CarState,
                                        #    length(gen.physical.cars))))
    s = gen.physical
    resize!(full_s.cars, length(s.cars))
    for i in 1:length(s.cars)
        behavior = rand(rng, gen.behaviors)
        full_s.cars[i] = CarState(s.cars[i], behavior)
    end
    return full_s
end

#=
type BasicMLReinvigorator <: POMCP.ParticleReinvigorator
    N::Int # target number of particles
    p_change::Float64 # (independent) probability of the behavior of each car changing
    behaviors::AbstractVector
    weights::WeightVec
end

# ideas for future reinvigorators
# 1) change behaviors of newer cars more often
# 2) change behaviors of cars at the edges of the lanes more often
# 3) change behaviors of the cars that are furthest from their observations most often

function POMCP.reinvigorate!(pc::ParticleCollection, r::BasicMLReinvigorator, old_node::POMCP.BeliefNode, a, o)
    while length(pc.particles) < r.N
        new_s =
    end
    return pc
end

function POMCP.handle_unseen_observation(r::BasicMLReinvigorator, old_node::POMCP.BeliefNode, a, o)
    pc = ParticleCollection{MLState}()

end
=#

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
"""
    timeevolution.schroedinger(tspan, psi0, H; fout)

Integrate Schroedinger equation to evolve states or compute propagators.

# Arguments
* `tspan`: Vector specifying the points of time for which output should be displayed.
* `psi0`: Initial state vector (can be a bra or a ket) or initial propagator.
* `H`: Arbitrary operator specifying the Hamiltonian.
* `fout=nothing`: If given, this function `fout(t, psi)` is called every time
        an output should be displayed. ATTENTION: The state `psi` is neither
        normalized nor permanent! It is still in use by the ode solver and
        therefore must not be changed.
"""
function schroedinger(tspan, psi0::T, H::AbstractOperator{B,B};
                fout::Union{Function,Nothing}=nothing,
                kwargs...) where {B,T<:Union{AbstractOperator{B,B},StateVector{B}}}
    dschroedinger_(t, psi, dpsi) = dschroedinger!(dpsi, H, psi)
    x0 = psi0.data
    state = copy(psi0)
    dstate = copy(psi0)
    integrate(tspan, dschroedinger_, x0, state, dstate, fout; kwargs...)
end

"""
    timeevolution.schroedinger_dynamic(tspan, psi0, f; fout)

Integrate time-dependent Schroedinger equation to evolve states or compute propagators.

# Arguments
* `tspan`: Vector specifying the points of time for which output should be displayed.
* `psi0`: Initial state vector (can be a bra or a ket) or initial propagator.
* `f`: Function `f(t, psi) -> H` returning the time and or state dependent Hamiltonian.
* `fout=nothing`: If given, this function `fout(t, psi)` is called every time
        an output should be displayed. ATTENTION: The state `psi` is neither
        normalized nor permanent! It is still in use by the ode solver and
        therefore must not be changed.
"""
function schroedinger_dynamic(tspan, psi0, f;
                fout::Union{Function,Nothing}=nothing,
                kwargs...)
    dschroedinger_(t, psi, dpsi) = dschroedinger_dynamic!(dpsi, f, psi, t)
    x0 = psi0.data
    state = copy(psi0)
    dstate = copy(psi0)
    integrate(tspan, dschroedinger_, x0, state, dstate, fout; kwargs...)
end

"""
    recast!(x,y)

Write the data stored in `y` into `x`, where either `x` or `y` is a quantum
object such as a [`Ket`](@ref) or an [`Operator`](@ref), and the other one is
a vector or a matrix with a matching size.
"""
recast!(psi::StateVector{B,D},x::D) where {B, D} = (psi.data = x);
recast!(x::D,psi::StateVector{B,D}) where {B, D} = nothing

"""
    dschroedinger!(dpsi, H, psi)

Update the increment `dpsi` in-place according to a Schrödinger equation
as `-im*H*psi`.

See also: [`dschroedinger_dynamic!`](@ref)
"""
function dschroedinger!(dpsi, H, psi)
    QuantumOpticsBase.mul!(dpsi,H,psi,eltype(psi)(-im),zero(eltype(psi)))
    return dpsi
end

function dschroedinger!(dpsi, H, psi::Bra)
    QuantumOpticsBase.mul!(dpsi,psi,H,eltype(psi)(im),zero(eltype(psi)))
    return dpsi
end

"""
    dschroedinger_dynamic!(dpsi, f, psi, t)

Compute the Hamiltonian as `H=f(t, psi)` and update `dpsi` according to a
Schrödinger equation as `-im*H*psi`.

See also: [`dschroedinger!`](@ref)
"""
function dschroedinger_dynamic!(dpsi, f, psi, t)
    H = f(t, psi)
    dschroedinger!(dpsi, H, psi)
end

function check_schroedinger(psi, H)
    check_multiplicable(H, psi)
    check_samebases(H)
end

function check_schroedinger(psi::Bra, H)
    check_multiplicable(psi, H)
    check_samebases(H)
end


text:
"""
Implementation of the FIFO push relabel algorithm with gap heuristic. Takes
approximately O(V^3) time.

Maintains the following auxillary arrays:
- height -> Stores the labels of all vertices
- count  -> Stores the number of vertices at each height
- excess -> Stores the difference between incoming and outgoing flow for all vertices
- active -> Stores the status of all vertices. (e(v)>0 => active[v] = true)
- Q      -> The FIFO queue that stores active vertices waiting to be discharged.

Requires arguments:

- residual_graph::DiGraph                # the input graph
- source::Int                            # the source vertex
- target::Int                            # the target vertex
- capacity_matrix::AbstractArray{T,2}    # edge flow capacities
"""

function push_relabel{T<:Number}(
    residual_graph::DiGraph,               # the input graph
    source::Int,                           # the source vertex
    target::Int,                           # the target vertex
    capacity_matrix::AbstractArray{T,2}    # edge flow capacities
    )

    n = nv(residual_graph)

    flow_matrix = zeros(T, n, n)

    height = zeros(Int, n)
    height[source] = n

    count = zeros(Int, 2*n+1)
    count[0+1] = n-1
    count[n+1] = 1

    excess = zeros(T, n)
    excess[source] = typemax(T)

    active = falses(n)
    active[source] = true
    active[target] = true

    Q = Array{Int,1}()
    sizehint!(Q, n)

    for v in fadj(residual_graph, source)
        push_flow!(residual_graph, source, v, capacity_matrix, flow_matrix, excess, height, active, Q)
    end

    while length(Q) > 0
        v = pop!(Q)
        active[v] = false
        discharge!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q)
    end

    return sum([flow_matrix[v,target] for v in badj(residual_graph, target) ]), flow_matrix
end

"""
Pushes inactive nodes into the queue and activates them.

Requires arguments:

- Q::AbstractArray{Int,1}
- v::Int
- active::AbstractArray{Bool,1}
- excess::AbstractArray{T,1}
"""

function enqueue_vertex!{T<:Number}(
    Q::AbstractArray{Int,1},
    v::Int,                                # input vertex
    active::AbstractArray{Bool,1},
    excess::AbstractArray{T,1}
    )
    if !active[v] && excess[v] > 0
        active[v] = true
        unshift!(Q, v)
    end
    return nothing
end

"""
Pushes as much flow as possible through the given edge.

Requires arguements:

- residual_graph::DiGraph              # the input graph
- u::Int                               # input from-vertex
- v::Int                               # input to-vetex
- capacity_matrix::AbstractArray{T,2}
- flow_matrix::AbstractArray{T,2}
- excess::AbstractArray{T,1}
- height::AbstractArray{Int,1}
- active::AbstractArray{Bool,1}
- Q::AbstractArray{Int,1}
"""

function push_flow!{T<:Number}(
    residual_graph::DiGraph,             # the input graph
    u::Int,                              # input from-vertex
    v::Int,                              # input to-vetex
    capacity_matrix::AbstractArray{T,2},
    flow_matrix::AbstractArray{T,2},
    excess::AbstractArray{T,1},
    height::AbstractArray{Int,1},
    active::AbstractArray{Bool,1},
    Q::AbstractArray{Int,1}
    )
    flow = min(excess[u], capacity_matrix[u,v] - flow_matrix[u,v])

    flow == 0 && return nothing
    height[u] <= height[v] && return nothing

    flow_matrix[u,v] += flow
    flow_matrix[v,u] -= flow

    excess[u] -= flow
    excess[v] += flow

    enqueue_vertex!(Q, v, active, excess)
    nothing
end

"""
Implements the gap heuristic. Relabels all vertices above a cutoff height.
Reduces the number of relabels required.

Requires arguments:

- residual_graph::DiGraph                # the input graph
- h::Int                                 # cutoff height
- excess::AbstractArray{T,1}
- height::AbstractArray{Int,1}
- active::AbstractArray{Bool,1}
- count::AbstractArray{Int,1}
- Q::AbstractArray{Int,1}
"""

function gap!{T<:Number}(
    residual_graph::DiGraph,               # the input graph
    h::Int,                                # cutoff height
    excess::AbstractArray{T,1},
    height::AbstractArray{Int,1},
    active::AbstractArray{Bool,1},
    count::AbstractArray{Int,1},
    Q::AbstractArray{Int,1}                # FIFO queue
    )
    n = nv(residual_graph)
    for v in vertices(residual_graph)
        height[v] < h && continue
        count[height[v]+1] -= 1
        height[v] = max(height[v], n + 1)
        count[height[v]+1] += 1
        enqueue_vertex!(Q, v, active, excess)
    end
    nothing
end

"""
Relabels a vertex with respect to its neighbors, to produce an admissable
edge.

Requires arguments:

- residual_graph::DiGraph                 # the input graph
- v::Int                                  # input vertex to be relabeled
- capacity_matrix::AbstractArray{T,2}
- flow_matrix::AbstractArray{T,2}
- excess::AbstractArray{T,1}
- height::AbstractArray{Int,1}
- active::AbstractArray{Bool,1}
- count::AbstractArray{Int,1}
- Q::AbstractArray{Int,1}
"""

function relabel!{T<:Number}(
    residual_graph::DiGraph,                # the input graph
    v::Int,                                 # input vertex to be relabeled
    capacity_matrix::AbstractArray{T,2},
    flow_matrix::AbstractArray{T,2},
    excess::AbstractArray{T,1},
    height::AbstractArray{Int,1},
    active::AbstractArray{Bool,1},
    count::AbstractArray{Int,1},
    Q::AbstractArray{Int,1}
    )
    n = nv(residual_graph)
    count[height[v]+1] -= 1
    height[v] = 2*n
    for to in fadj(residual_graph, v)
        if capacity_matrix[v,to] > flow_matrix[v,to]
            height[v] = min(height[v], height[to]+1)
        end
    end
    count[height[v]+1] += 1
    enqueue_vertex!(Q, v, active, excess)
    nothing
end

"""
Drains the excess flow out of a vertex. Runs the gap heuristic or relabels the
vertex if the excess remains non-zero.

Requires arguments:

- residual_graph::DiGraph                 # the input graph
- v::Int                                  # vertex to be discharged
- capacity_matrix::AbstractArray{T,2}
- flow_matrix::AbstractArray{T,2}
- excess::AbstractArray{T,1}
- height::AbstractArray{Int,1}
- active::AbstractArray{Bool,1}
- count::AbstractArray{Int,1}
- Q::AbstractArray{Int,1}
"""
function discharge!{T<:Number}(
    residual_graph::DiGraph,                # the input graph
    v::Int,                                 # vertex to be discharged
    capacity_matrix::AbstractArray{T,2},
    flow_matrix::AbstractArray{T,2},
    excess::AbstractArray{T,1},
    height::AbstractArray{Int,1},
    active::AbstractArray{Bool,1},
    count::AbstractArray{Int,1},
    Q::AbstractArray{Int,1}                 # FIFO queue
    )
    for to in fadj(residual_graph, v)
        excess[v] == 0 && break
        push_flow!(residual_graph, v, to, capacity_matrix, flow_matrix, excess, height, active, Q)
    end

    if excess[v] > 0
        if count[height[v]+1] == 1
            gap!(residual_graph, height[v], excess, height, active, count, Q)
        else
            relabel!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q)
        end
    end
    nothing
end


text:
abstract type AbstractConvolutionKernel end

immutable GaussianKernel <: AbstractConvolutionKernel
    Σ::AbstractPDMat
    SquaredExponentialKernel(Σ::AbstractPDMat) = new(Σ)
    SquaredExponentialKernel(Σ::Vector) = new(PDiagMat(Σ))
    SquaredExponentialKernel(Σ::Array) = new(PDMat(Σ))
end

const SquaredExponentialKernel = GaussianKernel

dim(kern::AbstractConvolutionKernel) = kern.Σ.dim

## Calculate weights
function conv_wt(kern::SquaredExponentialKernel, d::Array)
    if size(d, 1) != dim(kern)
        throw(DimensionMismatch)
    end
    exp.(-0.5 * invquad(kern.Σ, d))
end