text:
function device!(
    x,
    output_ode,
    voltage_r,
    voltage_i,
    current_r,
    current_i,
    ix_range::UnitRange{Int},
    ode_range::UnitRange{Int},
    dynamic_device::DynG,
    inputs::SimulationInputs,
) where {DynG <: PSY.DynamicGenerator}
    #Obtain local device states
    device_states = @view x[ix_range]

    #Obtain references
    sys_Sbase = get_base_power(inputs)
    sys_f0 = get_base_frequency(inputs)
    sys_ω = get_ω_sys(inputs)

    #Update Voltage data
    get_inner_vars(dynamic_device)[VR_gen_var] = voltage_r[1]
    get_inner_vars(dynamic_device)[VI_gen_var] = voltage_i[1]

    #Obtain ODEs and Mechanical Power for Turbine Governor
    mdl_tg_ode!(device_states, view(output_ode, ode_range), sys_ω, dynamic_device)

    #Obtain ODEs for AVR
    mdl_pss_ode!(device_states, view(output_ode, ode_range), sys_ω, dynamic_device)

    #Obtain ODEs for AVR
    mdl_avr_ode!(device_states, view(output_ode, ode_range), dynamic_device)

    #Obtain ODEs for Machine
    mdl_machine_ode!(
        device_states,
        view(output_ode, ode_range),
        current_r,
        current_i,
        sys_Sbase,
        sys_f0,
        dynamic_device,
    )

    #Obtain ODEs for PSY.Shaft
    mdl_shaft_ode!(
        device_states,
        view(output_ode, ode_range),
        sys_f0,
        sys_ω,
        dynamic_device,
    )

    return
end

function device!(
    voltage_r,
    voltage_i,
    current_r,
    current_i,
    device::PSY.Source,
    ::SimulationInputs,
)
    mdl_source!(voltage_r, voltage_i, current_r, current_i, device)
    return
end

function device!(
    voltage_r,
    voltage_i,
    current_r,
    current_i,
    device::PSY.PowerLoad,
    ::SimulationInputs,
)
    mdl_Zload!(voltage_r, voltage_i, current_r, current_i, device)
    return
end

function device!(
    x,
    output_ode::Vector{T},
    voltage_r,
    voltage_i,
    current_r,
    current_i,
    ix_range::UnitRange{Int},
    ode_range::UnitRange{Int},
    dynamic_device::DynI,
    inputs::SimulationInputs,
) where {DynI <: PSY.DynamicInverter, T <: Real}
    #Obtain local device states
    device_states = @view x[ix_range]

    #Obtain references
    Sbase = get_base_power(inputs)
    sys_f0 = get_base_frequency(inputs)
    sys_ω = get_ω_sys(inputs)

    #Update Voltage data
    get_inner_vars(dynamic_device)[VR_inv_var] = voltage_r[1]
    get_inner_vars(dynamic_device)[VI_inv_var] = voltage_i[1]

    #Update V_ref
    V_ref = PSY.get_ext(dynamic_device)[CONTROL_REFS][V_ref_index]
    get_inner_vars(dynamic_device)[V_oc_var] = V_ref

    #Obtain ODES for DC side
    mdl_DCside_ode!(
        device_states,
        view(output_ode, ode_range),
        sys_f0,
        sys_ω,
        dynamic_device,
    )

    #Obtain ODEs for PLL
    mdl_freq_estimator_ode!(
        device_states,
        view(output_ode, ode_range),
        sys_f0,
        sys_ω,
        dynamic_device,
    )

    #Obtain ODEs for OuterLoop
    mdl_outer_ode!(
        device_states,
        view(output_ode, ode_range),
        sys_f0,
        sys_ω,
        dynamic_device,
    )

    #Obtain inner controller ODEs and modulation commands
    mdl_inner_ode!(device_states, view(output_ode, ode_range), dynamic_device)

    #Obtain converter relations
    mdl_converter_ode!(dynamic_device)

    #Obtain ODEs for output filter
    mdl_filter_ode!(
        device_states,
        view(output_ode, ode_range),
        current_r,
        current_i,
        Sbase,
        sys_f0,
        sys_ω,
        dynamic_device,
    )

    return
end


text:
using StochasticIntegrators
using Test

@testset "$(rpad("Stochastic Runge-Kutta Tableaus",80))" begin

    @test typeof(TableauStochasticEuler()) <: TableauSERK
    @test typeof(TableauStochasticHeun()) <: TableauSERK
    @test typeof(TableauPlaten()) <: TableauSERK
    @test typeof(TableauBurrageR2()) <: TableauSERK
    @test typeof(TableauBurrageCL()) <: TableauSERK
    @test typeof(TableauBurrageE1()) <: TableauSERK
    @test typeof(TableauBurrageG5()) <: TableauSERK

    @test typeof(TableauStochasticGLRK(1)) <: TableauSIRK
    @test typeof(TableauStochasticGLRK(2)) <: TableauSIRK
    @test typeof(TableauStochasticDIRK()) <: TableauSIRK

    @test typeof(TableauStochasticSymplecticEuler()) <: TableauSIPRK
    @test typeof(TableauStochasticStoermerVerlet()) <: TableauSIPRK

    @test typeof(TableauStochasticLobattoIIIABD2()) <: TableauSISPRK
    @test typeof(TableauModifiedStochasticStoermerVerlet()) <: TableauSISPRK

    @test typeof(TableauRoesslerRS1()) <: TableauWERK
    @test typeof(TableauRoesslerRS2()) <: TableauWERK

    @test typeof(TableauSRKw1()) <: TableauWIRK
    @test typeof(TableauSRKw2()) <: TableauWIRK

end

ext: jl


text:
import PyPlot

function plot(r::basic_rod)

    # Setup quiver points
    xd = similar(r.d)
    xd[1,:] = r.x[1,:]
    xd[end,:] = r.x[end,:]
    xd[2:end-1,:] = 0.5*(r.x[2:end-2,:] + r.x[3:end-1,:])

    # Create a figure and axes
    fig = PyPlot.figure()
    ax = PyPlot.axes(projection="3d")
    ax.plot(r.x[:,1], r.x[:,2], r.x[:,3], "-ok")
    ax.quiver(xd[:,1], xd[:,2], xd[:,3],
                  r.d[:,1], r.d[:,2], r.d[:,3],
                  length=0.1, normalize=true)
    PyPlot.xlabel("x")
    PyPlot.ylabel("y")
    PyPlot.zlabel("z")
    PyPlot.pygui(true)
    display(fig)
end
#
# import Plots
# function plot(r::basic_rod)
#
#      # Setup quiver points
#      xd = similar(r.d)
#      xd[:,1] = r.x[:,1]
#      xd[:,end] = r.x[:,end]
#      xd[:,2:end-1] = 0.5*(r.x[:,2:end-2] + r.x[:,3:end-1])
#
#      # Create a figure and axes
#      q = Plots.quiver(xd[1,:], xd[2,:], xd[3,:],
#                    quiver = (r.d[1,:], r.d[2,:], r.d[3,:]), arrow=:arrow)
#
#      Plots.plot!(q, xlabel="x", ylabel = "x", zlabel = "z")
#      q
#  end


text:
#=
downsample.jl
Copyright 2019-03-05, Jeff Fessler, University of Michigan

see also MIRT/time/downsample.jl
=#

export downsample_dim1
export downsample1
export downsample2
export downsample3
export downsample

"""
    y = downsample_dim1(x, down ; warn::Bool)
downsample by factor `down` along first dimension by averaging

in
- `x [n1 (Nd)]`
- `down::Int` downsampling factor

option
- `warn::Bool` warn if noninteger multiple; default `isinteractive()`

out
- `y [n1/down (Nd)]`

Copyright 2019-03-05, Jeff Fessler, University of Michigan
"""
function downsample_dim1(
	x::AbstractArray{<:Number},
	down::Int ;
	warn::Bool = isinteractive(),
)

	dim = size(x)
	dim1 = dim[1]
	x = reshape(x, dim1, :) # [n1 *Nd]
	m1 = dim1 ÷ down
	if m1 * down < dim1
		warn && @warn("truncating input size $dim1 to $(m1 * down)")
		x = x[1:(m1*down),:]
	end
	y = reshape(x, down, :)
	y = sum(y, dims=1) / down # mean
	y = reshape(y, m1, dim[2:end]...)
#	out = similar(x, (m1, dim[2:end]...))
#	vec(out) .= vec(y) # failed attempt to help @inferred for 2D input arrays
	return y
end

"""
    y = downsample1(x, down ; warn=true)
downsample 1D vector by factor `down`

in
- `x [n1]`
- `down::Int` downsampling factor

option
- `warn::Bool` warn if noninteger multiple; default `isinteractive()`

out
- `y [n1/down]`

Copyright 2019-03-05, Jeff Fessler, University of Michigan
"""
function downsample1(
	x::AbstractVector{<:Number},
	down::Int ;
	warn::Bool = isinteractive(),
)

	dim = size(x)
	dim1 = dim[1]
	m1 = dim1 ÷ down
	if m1 * down < dim1
		warn && @warn("truncating input size $dim1 to $(m1 * down)")
		y = reshape((@view x[1:(m1*down)]), down, :)
	else
		y = reshape(x, down, :)
	end
	y = sum(y, dims=1) / down # mean
	y = reshape(y, m1, dim[2:end]...)
	return vec(y)
end

"""
    y = downsample2(x, down ; warn=true, T)

downsample by averaging by integer factors
in
- `x [nx ny]`
- `down` can be a scalar (same factor for both dimensions) or a `NTuple{2,Int}`

option
- `warn::Bool` warn if noninteger multiple; default `isinteractive()`
- `T::DataType` specify output eltype; default `eltype(x[1] / down[1])`

out
- `y [nx/down ny/down]`
"""
function downsample2(
	x::AbstractMatrix{<:Number},
	down::NTuple{2,Int} ;
	warn::Bool = isinteractive(),
	T::DataType = eltype(x[1] / down[1])
)

	idim = size(x)
	odim = idim .÷ down

	warn && any(odim .* down .!= idim) && @warn("truncating to $odim")

#	y = similar(x, odim) # fails!?
	y = Array{T}(undef, odim[1], odim[2])
	d1 = down[1]
	d2 = down[2]
	for i2=1:odim[2]
		for i1=1:odim[1]
			y[i1,i2] =
				sum(@view x[(i1-1)*d1 .+ (1:d1), (i2-1)*d2 .+ (1:d2)]) / d1 / d2
		end
	end

#=
	fun = (x, d) -> downsample_dim1(x, d, warn=warn)

#=
#	this old way returned an Adjoint type:
	y = fun(x, down[1])
	y = fun(y', down[end])'
=#

#	this way avoids the Adjoint:
	y = fun(x', down[end])' # doing adjoint first
	y = fun(y, down[1])

#=	failed attempt at helping @inferred:
	f1 = (x, d) -> downsample1(x, d, warn=warn)
	tmp = x'
	tmp = hcat([f1(tmp[:,n], down[end]) for n=1:size(tmp,2)]...)'
	y = hcat([f1(tmp[:,n], down[1]) for n=1:size(tmp,2)]...)
@shows y
=#
=#

	return y
end

downsample2(x::AbstractArray{<:Number,2}, down::Int ; args...) =
	downsample2(x, (down, down) ; args...)

"""
    y = downsample3(x, down ; warn=true, T)

downsample by averaging by integer factors
in
- `x [nx ny nz]`
- `down` can be a scalar (same factor for all dimensions) or a `NTuple{3,Int}`

option
- `warn::Bool` warn if noninteger multiple; default true
- `T::DataType` specify output eltype; default `eltype(x[1] / down[1])`

out
- `y [nx/down ny/down nz/down]`
"""
function downsample3(
	x::AbstractArray{<:Number,3},
	down::NTuple{3,Int} ;
	warn::Bool = isinteractive(),
	T::DataType = eltype(x[1] / down[1]),
)

	idim = size(x)
	odim = idim .÷ down

	warn && any(odim .* down .!= idim) && @warn("truncating to $odim")

	return downsample3_perm(x, Tuple(down)) # because it is faster
end

downsample3(x::AbstractArray{<:Number,3}, down::Int ; args...) =
	downsample3(x, (down, down, down) ; args...)

# this method is good for @inferred but is slower!
function downsample3_loop(
	x::AbstractArray{<:Number,3},
	down::NTuple{3,Int} ;
	T::DataType = eltype(x[1] / down[1]),
)

	odim = size(x) .÷ down

	y = Array{T}(undef, odim[1], odim[2], odim[3])
	d1 = down[1]
	d2 = down[2]
	d3 = down[3]
	d123 = d1 * d2 * d3

	for i3=1:odim[3]
		for i2=1:odim[2]
			for i1=1:odim[1]
				y[i1,i2,i3] = sum(@view x[(i1-1)*d1 .+ (1:d1),
					(i2-1)*d2 .+ (1:d2), (i3-1)*d3 .+ (1:d3)]) / d123
			end
		end
	end
	return y
end

# this method fails @inferred but is faster!
function downsample3_perm(x::AbstractArray{<:Number,3}, down::NTuple{3,Int})

	# down sample along each dimension
	y = downsample_dim1(x, down[1])
	y = downsample_dim1(permutedims(y, [2, 1, 3]), down[2])
	y = downsample_dim1(permutedims(y, [3, 2, 1]), down[3]) # [3 1 2] order
	y = permutedims(y, [2, 3, 1])

	return y
end


text:
using LinearAlgebra
using Test
using DistributedArrays
using Distributed
# # add at least 3 worker processes
if nworkers() < 3
    n = max(3, min(8, Sys.CPU_THREADS))
    addprocs(n; exeflags=`--check-bounds=yes`)
end
@assert nprocs() > 3
@assert nworkers() >= 3

@everywhere using DistributedArrays, SparseArrays, LinearAlgebra, Random, TimerOutputs
@everywhere using JOLI
@everywhere using SetIntersectionProjection

@everywhere mutable struct compgrid
  d :: Tuple
  n :: Tuple
end

@testset "SetIntersectionProjection" begin

  include("test_projectors.jl")
  include("test_setup_constraints.jl")

  include("test_TD_OPs.jl")
  include("test_prox_l2s!.jl")
  include("test_argmin_x.jl")
  include("test_update_y_l.jl")

  #parallel scripts
  include("test_update_y_l_parallel.jl")
  include("test_adapt_rho_gamma_parallel.jl")

  #linear algebra subroutines
  include("test_cg.jl")
  include("test_CDS_Mvp.jl")
  include("test_CDS_scaled_add.jl")
  include("test_Q_update.jl")

  #test full algorithms
  include("test_PARSDMM.jl")
  include("test_PARSDMM_parallel.jl")
  #include("test_PARSDMM_multilevel.jl")

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
module CalcDist

import StatsFuns: normpdf, normcdf, norminvcdf,
                  tdistpdf, tdistcdf, tdistinvcdf,
                  chisqpdf, chisqcdf, chisqinvcdf,
                  fdistpdf, fdistcdf, fdistinvcdf

export normalpdf, normalcdf, invNorm,
       tpdf, tcdf, invT,
       X2pdf, X2cdf, invX2,
       Fpdf, Fcdf, invF,
       binompdf, binomcdf,
       poissonpdf, poissoncdf,
       geometpdf, geometcdf,
       nCr, nPr

##
# Continuous
##

"""
   normalpdf(x, μ = 0, σ = 1) \n
      x: value
      μ: mean
      σ: standard deviation
   probability density function (Normal)
"""
function normalpdf(x, μ = 0, σ = 1)
   return normpdf(μ, σ, x)
end

"""
   normalcdf(lower, upper, μ = 0, σ = 1) \n
      lower: lower boundary
      upper: upper boundary
          μ: mean
          σ: standard deviation
   cumulative distribution function (Normal)
"""
function normalcdf(lower, upper, μ = 0, σ = 1)
   return normcdf(μ, σ, upper) - normcdf(μ, σ, lower)
end

"""
   invNorm(area, μ = 0, σ = 1) \n
      area: probability area
         μ: mean
         σ: standard deviation
   inverse cumulative distribution function (Normal)
"""
function invNorm(area, μ = 0, σ = 1)
   return norminvcdf(μ, σ, area)
end

"""
   tpdf(x, ν) \n
      x: value
      ν: degrees of freedom
   probability density function (Student's t)
"""
function tpdf(x, ν)
   return tdistpdf(ν, x)
end

"""
   tcdf(lower, upper, ν) \n
      lower: lower boundary
      upper: upper boundary
          ν: degrees of freedom
   cumulative distribution function (Student's t)
"""
function tcdf(lower, upper, ν)
   return tdistcdf(ν, upper) - tdistcdf(ν, lower)
end

"""
   invT(area, ν) \n
      area: probability area
         ν: degrees of freedom
   inverse cumulative distribution function (Student's t)
"""
function invT(area, ν)
   return tdistinvcdf(ν, area)
end

"""
   X2pdf(x, ν) \n
      x: value
      ν: degrees of freedom
   probability density function (Chi-squared)
"""
function X2pdf(x, ν)
   return chisqpdf(ν, x)
end

"""
   X2cdf(lower, upper, ν) \n
      lower: lower boundary
      upper: upper boundary
          ν: degrees of freedom
   cumulative distribution function (Chi-squared)
"""
function X2cdf(lower, upper, ν)
   return chisqcdf(ν, upper) - chisqcdf(ν, lower)
end

"""
   invX2(area, ν) \n
      area: probability area
         ν: degrees of freedom
   inverse cumulative distribution function (Chi-squared)
"""
function invX2(area, ν)
   return chisqinvcdf(ν, area)
end

"""
   Fpdf(x, nν, dν) \n
       x: value
      nν: numerator degrees of freedom
      dν: denominator degrees of freedom
   probability density function (F)
"""
function Fpdf(x, nν, dν)
   return fdistpdf(nν, dν, x)
end

"""
   Fcdf(lower, upper, nν, dν) \n
      lower: lower boundary
      upper: upper boundary
         nν: numerator degrees of freedom
         dν: denominator degrees of freedom
   cumulative distribution function (F)
"""
function Fcdf(lower, upper, nν, dν)
   return fdistcdf(nν, dν, upper) - fdistcdf(nν, dν, lower)
end

"""
   invF(area, nν, dν) \n
      area: probability area
        nν: numerator degrees of freedom
        dν: denominator degrees of freedom
   inverse cumulative distribution function (F)
"""
function invF(area, nν, dν)
   return fdistinvcdf(nν, dν, area)
end

##
# Discrete
##

"""
   binompdf(n, p, x = nothing) \n
      n: trials
      p: probability of success
      x: value
   probability mass function (Binomial)
"""
function binompdf(n, p, x = nothing)
   if n ≤ 0
      throw(DomainError(n, "\t[n > 0] must be true."))
   end
   if p < 0 || p > 1
      throw(DomainError(p, "\t[0 ≤ p ≤ 1] must be true."))
   end

   if isnothing(x)
      x = range(0, stop=n)
   end

   probs = Vector{Real}()

   for i in x
      push!(probs, binomial(n, i) * p^i * (1-p)^(n-i))
   end

   return probs
end

"""
   binomcdf(n, p, x = nothing) \n
      n: trials
      p: probability of success
      x: value
   cumulative distribution function (Binomial)
"""
function binomcdf(n, p, x = nothing)
   if n ≤ 0
      throw(DomainError(n, "\t[n > 0] must be true."))
   end
   if p < 0 || p > 1
      throw(DomainError(p, "\t[0 ≤ p ≤ 1] must be true."))
   end

   if isnothing(x)
      x = range(0, stop=n)
   end

   probs = Vector{Real}()

   for i in x
      push!(probs, sum([ binomial(n, i) * p^i * (1-p)^(n-i) for i in 0:i ]))
   end

   return probs
end

"""
   poissonpdf(λ, x) \n
      λ: mean, variance
      x: value
   probability mass function (Poisson)
"""
function poissonpdf(λ, x)
   if λ ≤ 0
      throw(DomainError(λ, "\t[λ > 0] must be true."))
   end

   probs = Vector{Real}()

   for i in x
      push!(probs, ℯ^(-λ) * (λ^i / factorial(i)))
   end

   return probs
end

"""
   poissoncdf(λ, x) \n
      λ: mean, variance
      x: value
   cumulative distribution function (Poisson)
"""
function poissoncdf(λ, x)
   if λ ≤ 0
      throw(DomainError(λ, "\t[λ > 0] must be true."))
   end

   probs = Vector{Real}()

   for i in x
      push!(probs, ℯ^(-λ) * sum([ λ^i / factorial(i) for i in 0:i ]))
   end

   return probs
end

"""
   geometpdf(p, x) \n
      p: probability of success
      x: trial on which 1st success occurs
   probability density function (Geometric)
"""
function geometpdf(p, x)
   if p < 0 || p > 1
      throw(DomainError(p, "\t[0 ≤ p ≤ 1] must be true."))
   end

   probs = Vector{Real}()

   for i in x
      push!(probs, p * (1-p)^(i-1))
   end

   return probs
end

"""
   geometcdf(p, x) \n
      p: probability of success
      x: trial on which 1st success occurs
   cumulative distribution function (Geometric)
"""
function geometcdf(p, x)
   if p < 0 || p > 1
      throw(DomainError(p, "\t[0 ≤ p ≤ 1] must be true."))
   end

   probs = Vector{Real}()

   for i in x
      push!(probs, 1 - (1-p)^i)
   end

   return probs
end

##
# Counting
##

"""
   nCr(n, r) \n
      return number of combinations of n items taken r at a time. n and r must be
      non-negative. both n and r can be vectors.
"""
function nCr(n, r)
   if n < 0
      throw(DomainError(n, "\t[n ≥ 0] must be true."))
   end

   combs = Vector{Real}()

   for i in r
      if i < 0
         throw(DomainError(i, "\t[r ≥ 0] must be true."))
      end
      push!(combs, factorial(n) / ( factorial(i) * factorial(n-i) ))
   end

   return combs
end

function nCr(n::Vector{Int}, r)
   if r < 0
      throw(DomainError(r, "\t[r ≥ 0] must be true."))
   end

   combs = Vector{Real}()

   for i in n
      if i < 0
         throw(DomainError(i, "\t[n ≥ 0] must be true."))
      end
      push!(combs, factorial(i) / ( factorial(r) * factorial(i-r) ))
   end

   return combs
end

"""
   nPr(n, r)
      return number of permutations of n items taken r at a time. n and r must be
      non-negative. both n and r can be vectors.
"""
function nPr(n, r)
   if n < 0
      throw(DomainError(n, "\t[n ≥ 0] must be true."))
   end

   perms = Vector{Real}()

   for i in r
      if i < 0
         throw(DomainError(i, "\t[r ≥ 0] must be true."))
      end
      push!(perms, factorial(n) / factorial(n-i))
   end

   return perms
end

function nPr(n::Vector{Int}, r)
   if r < 0
      throw(DomainError(r, "\t[r ≥ 0] must be true."))
   end

   perms = Vector{Real}()

   for i in n
      if i < 0
         throw(DomainError(i, "\t[n ≥ 0] must be true."))
      end
      push!(perms, factorial(i) / factorial(i-r))
   end

   return perms
end

end # module


text:
## ----------------------------------------------------------------------------
# _GLOBAL_CONF
const _GLOBAL_CONF = Dict()

function _init_globals()
    empty!(_GLOBAL_CONF)
    _GLOBAL_CONF[:CACHE_DIR] = pwd()
    _GLOBAL_CONF[:VERBOSE] = false
    _GLOBAL_CONF[:PRINT_FUN] = Base.println
    _GLOBAL_CONF[:LOAD_FUN] = _load
    _GLOBAL_CONF[:SAVE_FUN] = _save
    _GLOBAL_CONF[:ADD_TAG] = false
    _GLOBAL_CONF[:MK_DIR] = true
end

## ----------------------------------------------------------------------------
_toupper(s::Symbol) = Symbol(uppercase(string(s)))
global_conf() = _GLOBAL_CONF
global_conf(id::Symbol) = getindex(_GLOBAL_CONF, _toupper(id))
global_conf(id::Symbol, val) = (setindex!(_GLOBAL_CONF, val, _toupper(id)); val)


text:
using DataFrames
using Random

struct Individual
    id # ::Any
    amt
    dosing
    observations # ::Vector
    covariates
    function Individual(;id = rand(UInt16), amt, dosing, observations, covariates...)
        return new(id, amt, dosing, observations, covariates)
    end
end

struct Cohort
    population
    length
end

function Cohort(population)
    return Cohort(population, length(population))
end

function cohort(individuals, f)
    population = map(f, individuals)
    return Cohort(population)
end

function cohort!(individuals, f)
    population = map!(f, individuals)
    return Cohort(population)
end

# struct Run
#     id::Symbol
#     amt::Symbol
#     obs::Symbol
#     covariates::Tuple{Union{Symbol, String}}
#     data::DataFrame
# end

# @inline function getter(run::Run, field::Union{Symbol, String})
#     return run.data[getfield(run, field)]
# end

# @inline function id(run::Run)
#     return getter(run, :id)
# end

# @inline function amt(run::Run)
#     return getter(run, :amt)
# end

# @inline function obs(run::Run)
#     return getter(run, :obs)
# end

function DataFrames.DataFrame(cohort::Cohort)
    individuals = map(ind -> DataFrame(ind), cohort.population)
    return vcat(individuals...)
end

function DataFrames.DataFrame(individual::Individual)
    time = vcat(individual.dosing, individual.observations)
    entries = length(time)
    amt = zeros(entries)
    amt[1:length(individual.dosing)] .= coalesce(individual.amt, 0)
    indices = sortperm(time)
    covs = map((key, value) -> key => fill(value, entries), keys(individual.covariates), values(individual.covariates))
    id = fill(individual.id, entries)
    time = time[indices]
    amt = amt[indices]
    dv =  zeros(length(time))
    return DataFrame(
        :ID => id,
        :TIME => time,
        :EVID => Vector{Int64}(amt .> 0), # What ioamtis zero for EVID == 1
        :AMT => amt,
        :DV => dv,
        covs...
    )
end


text:
module TestIO
    using Base.Test
    using DataFrames, Compat
    using LaTeXStrings

    #test_group("We can read various file types.")

    data = joinpath(dirname(@__FILE__), "data")

    filenames = ["$data/blanklines/blanklines.csv",
                 "$data/compressed/movies.csv.gz",
                 "$data/headeronly/headeronly.csv",
                 "$data/newlines/embedded_os9.csv",
                 "$data/newlines/embedded_osx.csv",
                 "$data/newlines/embedded_windows.csv",
                 "$data/padding/space_after_delimiter.csv",
                 "$data/padding/space_around_delimiter.csv",
                 "$data/padding/space_before_delimiter.csv",
                 "$data/quoting/empty.csv",
                 "$data/quoting/escaping.csv",
                 "$data/quoting/quotedcommas.csv",
                 "$data/scaling/10000rows.csv",
                 "$data/utf8/corrupt_utf8.csv",
                 "$data/utf8/short_corrupt_utf8.csv",
                 "$data/utf8/utf8.csv"]

    for filename in filenames
        try
            df = readtable(filename)
        catch
            error(@sprintf "Failed to read %s\n" filename)
        end
    end

    #test_group("We get the right size, types, values for a basic csv.")

    filename = "$data/scaling/movies.csv"
    df = readtable(filename)

    @test size(df) == (58788, 25)

    @test df[1, 1] === Nullable(1)
    @test isequal(df[1, 2], Nullable("\$"))
    @test df[1, 3] === Nullable(1971)
    @test df[1, 4] === Nullable(121)
    @test isnull(df[1, 5])
    @test df[1, 6] === Nullable(6.4)
    @test df[1, 7] === Nullable(348)
    @test df[1, 8] === Nullable(4.5)
    @test df[1, 9] === Nullable(4.5)
    @test df[1, 10] === Nullable(4.5)
    @test df[1, 11] === Nullable(4.5)
    @test df[1, 12] === Nullable(14.5)
    @test df[1, 13] === Nullable(24.5)
    @test df[1, 14] === Nullable(24.5)
    @test df[1, 15] === Nullable(14.5)
    @test df[1, 16] === Nullable(4.5)
    @test df[1, 17] === Nullable(4.5)
    @test isequal(df[1, 18], Nullable(""))
    @test df[1, 19] === Nullable(0)
    @test df[1, 20] === Nullable(0)
    @test df[1, 21] === Nullable(1)
    @test df[1, 22] === Nullable(1)
    @test df[1, 23] === Nullable(0)
    @test df[1, 24] === Nullable(0)
    @test df[1, 25] === Nullable(0)

    @test df[end, 1] === Nullable(58788)
    @test isequal(df[end, 2], Nullable("xXx: State of the Union"))
    @test df[end, 3] === Nullable(2005)
    @test df[end, 4] === Nullable(101)
    @test df[end, 5] === Nullable(87000000)
    @test df[end, 6] === Nullable(3.9)
    @test df[end, 7] === Nullable(1584)
    @test df[end, 8] === Nullable(24.5)
    @test df[end, 9] === Nullable(4.5)
    @test df[end, 10] === Nullable(4.5)
    @test df[end, 11] === Nullable(4.5)
    @test df[end, 12] === Nullable(4.5)
    @test df[end, 13] === Nullable(14.5)
    @test df[end, 14] === Nullable(4.5)
    @test df[end, 15] === Nullable(4.5)
    @test df[end, 16] === Nullable(4.5)
    @test df[end, 17] === Nullable(14.5)
    @test isequal(df[end, 18], Nullable("PG-13"))
    @test df[end, 19] === Nullable(1)
    @test df[end, 20] === Nullable(0)
    @test df[end, 21] === Nullable(0)
    @test df[end, 22] === Nullable(0)
    @test df[end, 23] === Nullable(0)
    @test df[end, 24] === Nullable(0)
    @test df[end, 25] === Nullable(0)

    #test_group("readtable handles common separators and infers them from extensions.")

    df1 = readtable("$data/separators/sample_data.csv")
    df2 = readtable("$data/separators/sample_data.tsv")
    df3 = readtable("$data/separators/sample_data.wsv")
    df4 = readtable("$data/separators/sample_data_white.txt", separator = ' ')

    @test isequal(df1, df2)
    @test isequal(df2, df3)
    @test isequal(df3, df4)

    readtable("$data/quoting/quotedwhitespace.txt", separator = ' ')

    #test_group("readtable handles common newlines.")

    df = readtable("$data/newlines/os9.csv")
    @test isequal(readtable("$data/newlines/osx.csv"), df)
    @test isequal(readtable("$data/newlines/windows.csv"), df)

    @test isequal(df, readtable("$data/newlines/os9.csv", skipblanks = false))
    @test isequal(df, readtable("$data/newlines/osx.csv", skipblanks = false))
    @test isequal(df, readtable("$data/newlines/windows.csv", skipblanks = false))

    #test_group("readtable treats rows as specified.")

    df1 = readtable("$data/comments/before_after_data.csv", allowcomments = true)
    df2 = readtable("$data/comments/middata.csv", allowcomments = true)
    df3 = readtable("$data/skiplines/skipfront.csv", skipstart = 3)
    df4 = readtable("$data/skiplines/skipfront.csv", skipstart = 4, header = false)
    names!(df4, names(df1))
    df5 = readtable("$data/comments/before_after_data_windows.csv", allowcomments = true)
    df6 = readtable("$data/comments/middata_windows.csv", allowcomments = true)
    df7 = readtable("$data/skiplines/skipfront_windows.csv", skipstart = 3)
    df8 = readtable("$data/skiplines/skipfront_windows.csv", skipstart = 4, header = false)
    names!(df8, names(df1))
    # df9 = readtable("$data/skiplines/skipfront.csv", skipstart = 3, skiprows = 5:6)
    # df10 = readtable("$data/skiplines/skipfront.csv", skipstart = 3, header = false, skiprows = [4, 6])
    # names!(df10, names(df1))

    @test isequal(df2, df1)
    @test isequal(df3, df1)
    @test isequal(df4, df1)

    # Windows EOLS
    @test isequal(df5, df1)
    @test isequal(df6, df1)
    @test isequal(df7, df1)
    @test isequal(df8, df1)

    # @test isequal(df9, df1[3:end])
    # @test isequal(df10, df1[[1, 3:end]])

    function normalize_eol!(df)
        for (name, col) in eachcol(df)
            if eltype(col) <: AbstractString
                df[name] = map(s -> replace(s, "\r\n", "\n"), col)
            elseif eltype(col) <: Nullable && eltype(eltype(col)) <: AbstractString
                df[name] = map(s -> replace(get(s), "\r\n", "\n"), col)
            end
        end
        df
    end

    osxpath = "$data/skiplines/complex_osx.csv"
    winpath = "$data/skiplines/complex_windows.csv"

    opts1 = @compat Dict{Any,Any}(:allowcomments => true)
    opts2 = @compat Dict{Any,Any}(:skipstart => 4, :skiprows => [6, 7, 12, 14, 17], :skipblanks => false)

    df1 = readtable(osxpath; opts1...)
    # df2 = readtable(osxpath; opts2...)
    df1w = readtable(winpath; opts1...)
    # df2w = readtable(winpath; opts2...)

    # Normalize line endings in both and test equality
    @test isequal(normalize_eol!(df1w), normalize_eol!(df1))
    # @test normalize_eol!(df2w) == df1

    opts1[:nrows] = 3
    opts2[:nrows] = 3

    @test isequal(normalize_eol!(readtable(osxpath; opts1...)), df1[1:3, :])
    # @test isequalreadtable(osxpath; opts2...), df1[1:3, :]
    @test isequal(normalize_eol!(readtable(winpath; opts1...)), df1[1:3, :])
    # @test isequalreadtable(winpath; opts2...), df1[1:3, :])

    #test_group("readtable handles custom delimiters.")

    readtable("$data/skiplines/skipfront.csv", allowcomments = true, commentmark = '%')

    readtable("$data/separators/sample_data.csv", quotemark = Char[])
    @test_throws ErrorException readtable("$data/newlines/embedded_osx.csv", quotemark = Char[])
    df = readtable("$data/quoting/single.csv", quotemark = ['\''])
    @test isequal(df, readtable("$data/quoting/mixed.csv", quotemark = ['\'', '"']))

    # df = readtable("$data/decimal/period.csv")
    # @test isequaldf[2, :A], 0.3)
    # @test isequaldf[2, :B], 4.0)

    # @test isequal(df, readtable("$data/decimal/comma.tsv", decimal = ','))

    #test_group("readtable column names.")

    ns = [:Var1, :Var2, :Var3, :Var4, :Var5]
    df = readtable("$data/typeinference/mixedtypes.csv")
    names!(df, ns)
    @test isequal(df, readtable("$data/typeinference/mixedtypes.csv", names = ns))

    df = readtable("$data/separators/sample_data.csv", header = false, names = ns[1:3])
    @test isequal(df[1, :Var1], Nullable(0))
    df = readtable("$data/separators/sample_data.csv", names = ns[1:3])
    @test isequal(df[1, :Var1], Nullable(1))

    #test_group("Properties of data frames returned by readtable method.")

    # Readtable ignorepadding
    io = IOBuffer("A , \tB  , C\n1 , \t2, 3\n")
    @test isequal(readtable(io, ignorepadding = true), DataFrame(A = 1, B = 2, C = 3))

    # Readtable c-style escape options

    df = readtable("$data/escapes/escapes.csv", allowescapes = true)
    @test isequal(df[1, :V], Nullable("\t\r\n"))
    @test isequal(df[2, :V], Nullable("\\\\t"))
    @test isequal(df[3, :V], Nullable("\\\""))

    df = readtable("$data/escapes/escapes.csv")
    @test isequal(df[1, :V], Nullable("\\t\\r\\n"))
    @test isequal(df[2, :V], Nullable("\\\\t"))
    @test isequal(df[3, :V], Nullable("\\\""))

    # df = readtable("$data/escapes/escapes.csv", escapechars = ['"'], nrows = 2)
    # @test isequal(df[1, :V], "\\t\\r\\n")
    # @test isequal(df[2, :V], "\\\\\\\\t")

    # Readtable with makefactors active should only make factors from columns
    # of strings.
    filename = "$data/factors/mixedvartypes.csv"
    df = readtable(filename, makefactors = true)

    @test isa(df[:factorvar], NullableCategoricalArray{Compat.UTF8String,1})
    @test isa(df[:floatvar], NullableArray{Float64,1})

    # Readtable shouldn't silently drop data when reading highly compressed gz.
    df = readtable("$data/compressed/1000x2.csv.gz")
    @test size(df) == (1000, 2)

    # Readtable type inference
    filename = "$data/typeinference/bool.csv"
    df = readtable(filename)
    @test isa(df[:Name], NullableArray{Compat.UTF8String,1})
    @test isa(df[:IsMale], NullableArray{Bool,1})
    @test get(df[:IsMale][1])
    @test !get(df[:IsMale][4])

    filename = "$data/typeinference/standardtypes.csv"
    df = readtable(filename)
    @test isa(df[:IntColumn], NullableArray{Int,1})
    @test isa(df[:IntlikeColumn], NullableArray{Float64,1})
    @test isa(df[:FloatColumn], NullableArray{Float64,1})
    @test isa(df[:BoolColumn], NullableArray{Bool,1})
    @test isa(df[:StringColumn], NullableArray{Compat.UTF8String,1})

    filename = "$data/typeinference/mixedtypes.csv"
    df = readtable(filename)
    @test isa(df[:c1], NullableArray{Compat.UTF8String,1})
    @test isequal(df[:c1][1], Nullable("1"))
    @test isequal(df[:c1][2], Nullable("2.0"))
    @test isequal(df[:c1][3], Nullable("true"))
    @test isa(df[:c2], NullableArray{Float64,1})
    @test isequal(df[:c2][1], Nullable(1.0))
    @test isequal(df[:c2][2], Nullable(3.0))
    @test isequal(df[:c2][3], Nullable(4.5))
    @test isa(df[:c3], NullableArray{Compat.UTF8String,1})
    @test isequal(df[:c3][1], Nullable("0"))
    @test isequal(df[:c3][2], Nullable("1"))
    @test isequal(df[:c3][3], Nullable("f"))
    @test isa(df[:c4], NullableArray{Bool,1})
    @test isequal(df[:c4][1], Nullable(true))
    @test isequal(df[:c4][2], Nullable(false))
    @test isequal(df[:c4][3], Nullable(true))
    @test isa(df[:c5], NullableArray{Compat.UTF8String,1})
    @test isequal(df[:c5][1], Nullable("False"))
    @test isequal(df[:c5][2], Nullable("true"))
    @test isequal(df[:c5][3], Nullable("true"))

    # Readtable defining column types
    filename = "$data/definedtypes/mixedvartypes.csv"

    df = readtable(filename)
    @test isa(df[:n], NullableArray{Int,1})
    @test isequal(df[:n][1], Nullable(1))
    @test isa(df[:s], NullableArray{Compat.UTF8String,1})
    @test isequal(df[:s][1], Nullable("text"))
    @test isa(df[:f], NullableArray{Float64,1})
    @test isequal(df[:f][1], Nullable(2.3))
    @test isa(df[:b], NullableArray{Bool,1})
    @test isequal(df[:b][1], Nullable(true))

    df = readtable(filename, eltypes = [Int64, Compat.UTF8String, Float64, Bool])
    @test isa(df[:n], NullableArray{Int64,1})
    @test isequal(df[:n][1], Nullable(1))
    @test isa(df[:s], NullableArray{Compat.UTF8String,1})
    @test isequal(df[:s][1], Nullable("text"))
    @test isequal(df[:s][4], Nullable("text ole"))
    @test isa(df[:f], NullableArray{Float64,1})
    @test isequal(df[:f][1], Nullable(2.3))
    @test isa(df[:b], NullableArray{Bool,1})
    @test isequal(df[:b][1], Nullable(true))
    @test isequal(df[:b][2], Nullable(false))

    df = readtable(filename, eltypes = [Int64, Compat.UTF8String, Float64, Compat.UTF8String])
    @test isa(df[:n], NullableArray{Int64,1})
    @test isequal(df[:n][1], Nullable(1.0))
    @test isnull(df[:s][3])
    @test isa(df[:f], NullableArray{Float64,1})
    # Float are not converted to int
    @test isequal(df[:f][1], Nullable(2.3))
    @test isequal(df[:f][2], Nullable(0.2))
    @test isequal(df[:f][3], Nullable(5.7))
    @test isa(df[:b], NullableArray{Compat.UTF8String,1})
    @test isequal(df[:b][1], Nullable("T"))
    @test isequal(df[:b][2], Nullable("FALSE"))

    # Readtable name normalization
    abnormal = "\u212b"
    ns = [:Å, :_B_C_, :_end]
    @test !in(Symbol(abnormal), ns)

    io = IOBuffer(abnormal*",%_B*\tC*,end\n1,2,3\n")
    @test names(readtable(io)) == ns

    # With normalization disabled
    io = IOBuffer(abnormal*",%_B*\tC*,end\n1,2,3\n")
    @test names(readtable(io, normalizenames=false)) == [Symbol(abnormal),Symbol("%_B*\tC*"),:end]

    # Test writetable with Nullable() and compare to the results
    tf = tempname()
    isfile(tf) && rm(tf)
    df = DataFrame(A = NullableArray(Nullable{Int}[1,Nullable()]),
                   B = NullableArray(Nullable{String}["b", Nullable()]))
    writetable(tf, df)
    @test readcsv(tf) == ["A" "B"; 1 "b"; "NULL" "NULL"]

    # Test writetable with nastring set and compare to the results
    isfile(tf) && rm(tf)
    writetable(tf, df, nastring="none")
    @test readcsv(tf) == ["A" "B"; 1 "b"; "none" "none"]
    rm(tf)

    # Test writetable with append
    df1 = DataFrame(a = NullableArray([1, 2, 3]), b = NullableArray([4, 5, 6]))
    df2 = DataFrame(a = NullableArray([1, 2, 3]), b = NullableArray([4, 5, 6]))
    df3 = DataFrame(a = NullableArray([1, 2, 3]), c = NullableArray([4, 5, 6])) # 2nd column mismatch
    df3b = DataFrame(a = NullableArray([1, 2, 3]), b = NullableArray([4, 5, 6]), c = NullableArray([4, 5, 6])) # number of columns mismatch

    # Would use joinpath(tempdir(), randstring()) to get around tempname
    # creating a file on Windows, but Julia 0.3 has no srand() to unset the
    # seed set in test/data.jl -- annoying for local testing.
    tf = tempname()
    isfile(tf) && rm(tf)

    # Written as normal if file doesn't exist
    writetable(tf, df1, append = true)
    @test isequal(readtable(tf), df1)

    # Written as normal if file is empty
    open(io -> print(io, ""), tf, "w")
    writetable(tf, df1, append = true)
    @test isequal(readtable(tf), df1)

    # Appends to existing file if append == true
    writetable(tf, df1)
    writetable(tf, df2, header = false, append = true)
    @test isequal(readtable(tf), vcat(df1, df2))

    # Overwrites file if append == false
    writetable(tf, df1)
    writetable(tf, df2)
    @test isequal(readtable(tf), df2)

    # Enforces matching column names iff append == true && header == true
    writetable(tf, df2)
    @test_throws KeyError writetable(tf, df3, append = true)
    writetable(tf, df3, header = false, append = true)

    # Enforces matching column count if append == true
    writetable(tf, df3)
    @test_throws DimensionMismatch writetable(tf, df3b, header = false, append = true)
    rm(tf)

    # Quotemarks are escaped
    tf = tempname()
    isfile(tf) && rm(tf)

    df = DataFrame(a = ["who's"]) # We have a ' in our string

    # Make sure the ' doesn't get escaped for no reason
    writetable(tf, df)
    @test isequal(readtable(tf), df)

    # Make sure the ' does get escaped when needed
    writetable(tf, df, quotemark='\'')
    @test readstring(tf) == "'a'\n'who\\'s'\n"
    rm(tf)

    ### Tests for nonstandard string literals
    # Test basic @csv_str usage
    df1 = csv"""
        name,  age, squidPerWeek
        Alice,  36,         3.14
        Bob,    24,         0
        Carol,  58,         2.71
        Eve,    49,         7.77
        """
    @test size(df1) == (4, 3)
    @test names(df1) == [:name, :age, :squidPerWeek]
    @test isequal(df1[1], NullableArray(["Alice","Bob","Carol","Eve"]))
    @test isequal(df1[2], NullableArray([36,24,58,49]))
    @test isequal(df1[3], NullableArray([3.14,0,2.71,7.77]))
    @test isa(df1[1], NullableArray{Compat.UTF8String,1})

    # Test @wsv_str
    df2 = wsv"""
        name  age squidPerWeek
        Alice  36         3.14
        Bob    24         0
        Carol  58         2.71
        Eve    49         7.77
        """
    @test isequal(df2, df1)

    # Test @tsv_str
    df3 = tsv"""
        name	age	squidPerWeek
        Alice	36	3.14
        Bob	24	0
        Carol	58	2.71
        Eve	49	7.77
        """
    @test isequal(df3, df1)

    # csv2 can't be tested until non-'.' decimals are implemented
    #df4 = csv2"""
    #    name;  age; squidPerWeek
    #    Alice;  36;         3,14
    #    Bob;    24;         0
    #    Carol;  58;         2,71
    #    Eve;    49;         7,77
    #    """
    #@test isequal(df4, df1)

    # Test 'f' flag
    df5 = csv"""
        name,  age, squidPerWeek
        Alice,  36,         3.14
        Bob,    24,         0
        Carol,  58,         2.71
        Eve,    49,         7.77
        """f
    @test isa(df5[1], NullableCategoricalArray{Compat.UTF8String,1})

    # Test 'c' flag
    df6 = csv"""
        name,  age, squidPerWeek
        Alice,  36,         3.14
        Bob,    24,         0
        #Carol,  58,         2.71
        Eve,    49,         7.77
        """c
    @test isequal(df6, df1[[1,2,4],:])

    # Test 'H' flag
    df7 = csv"""
        Alice,  36,         3.14
        Bob,    24,         0
        Carol,  58,         2.71
        Eve,    49,         7.77
        """H
    @test names(df7) == [:x1,:x2,:x3]
    names!(df7, names(df1))
    @test isequal(df7, df1)

    # Test multiple flags at once
    df8 = csv"""
        Alice,  36,         3.14
        Bob,    24,         0
        #Carol,  58,         2.71
        Eve,    49,         7.77
        """fcH
    @test isa(df8[1], NullableCategoricalArray{Compat.UTF8String,1})
    @test names(df8) == [:x1,:x2,:x3]
    names!(df8, names(df1))
    @test isequal(df8, df1[[1,2,4],:])

    # Test invalid flag
    # Need to wrap macro call inside eval to prevent the error from being
    # thrown prematurely
    @test_throws ArgumentError eval(:(csv"foo,bar"a))

    # Test LaTeX export
    df = DataFrame(A = 1:4,
                   B = ["\$10.0", "M&F", "A~B", "\\alpha"],
                   C = [L"\alpha", L"\beta", L"\gamma", L"\sum_{i=1}^n \delta_i"],
                   D = [1.0, 2.0, Nullable(), 3.0]
                   )
    str = """
        \\begin{tabular}{r|cccc}
        \t& A & B & C & D\\\\
        \t\\hline
        \t1 & 1 & \\\$10.0 & \$\\alpha\$ & 1.0 \\\\
        \t2 & 2 & M\\&F & \$\\beta\$ & 2.0 \\\\
        \t3 & 3 & A\\textasciitilde{}B & \$\\gamma\$ &  \\\\
        \t4 & 4 & \\textbackslash{}alpha & \$\\sum_{i=1}^n \\delta_i\$ & 3.0 \\\\
        \\end{tabular}
        """
    @test reprmime(MIME("text/latex"), df) == str

    #Test HTML output for IJulia and similar
    df = DataFrame(Fish = ["Suzy", "Amir"], Mass = [1.5, Nullable()])
    io = IOBuffer()
    show(io, "text/html", df)
    str = takebuf_string(io)
    @test str == "

	Fish	Mass

1	Suzy	1.5

2	Amir	#NULL

"

    # test limit attribute of IOContext is used
    df = DataFrame(a=collect(1:1000))
    ioc = IOContext(IOBuffer(), displaysize=(10, 10), limit=false)
    show(ioc, "text/html", df)
    @test length(takebuf_string(ioc.io)) > 10000

    io = IOBuffer()
    show(io, "text/html", df)
    @test length(takebuf_string(io)) < 10000

end


text:
export H2oModel, H2oRegressor, H2oClassifier

@with_kw mutable struct H2oModel <: BaseEstimator
    pyo::PyObject = PyNULL()
    isclf::Bool = true
    name::String = "glm"
    family::String = "gaussian"
    lambda_search::Bool = true
    alpha::Float32 = 0
    gamma::Float32 = -1
    hyper_param::Float32 = 1
    ntrees::Int = 50
    sample_rate::Float32 = 1
    col_sample_rate::Float32 = 1
    col_sample_rate_per_tree::Float32 = 1
    max_depth::Int = 5
    min_rows::Int = 1
    score_tree_interval::Int = 5
    stopping_rounds::Int = 5
    learn_rate::Float32 = 0.1
    epochs::Int = 10
    adaptive_rate::Bool = true
    activation::String = "rectifier_with_dropout"
    rho::Float32 = 0.99
    epsilon::Float32 = 1e-8
    input_dropout_ratio::Float32 = 0
    hidden::Vector = [100, 100]
    hidden_dropout_ratios::Vector = [0.5, 0.5]
    max_models::Int = 10
end

is_classifier(m::H2oModel) = m.isclf

support_multiclass(m::H2oModel) = m.name != "psvm"

function fit!(m::H2oModel, x, y, w = nothing; columns = nothing)
    columns = something(columns, string.(1:size(x, 1)))
    @unpack name, isclf = m
    ENV["COLS"] = join(columns, '|')
    # init
    max_mem_size = @sprintf("%dG", 0.8 * Sys.total_memory() / 2^30)
    @imports h2o; h2o.init(max_mem_size = max_mem_size)
    @from h2o.estimators.glm imports H2OGeneralizedLinearEstimator
    @from h2o.estimators.psvm imports H2OSupportVectorMachineEstimator
    @from h2o.estimators.random_forest imports H2ORandomForestEstimator
    @from h2o.estimators.gbm imports H2OGradientBoostingEstimator
    @from h2o.estimators.xgboost imports H2OXGBoostEstimator
    @from h2o.estimators.deeplearning imports H2ODeepLearningEstimator
    @from h2o.automl imports H2OAutoML
    h2o.remove_all()
    h2o.api("POST /3/GarbageCollect")
    # create data
    dfx = DataFrame(pymat(x), columns = columns)
    dfy = DataFrame(vec(y), columns = ["label"])
    if isnothing(w)
        hdf = to_h2o(pdhcat(dfx, dfy))
    else
        dfw = DataFrame(vec(w), columns = ["weight"])
        hdf = to_h2o(pdhcat(dfx, dfy, dfw))
    end
    isclf && (hdf["label"] = hdf["label"].asfactor())
    # create model
    pyo = if name == "glm"
        @unpack family, lambda_search, alpha = m
        num_class = isclf ? length(unique(y)) : 1
        family = num_class == 1 ? family : num_class == 2 ? "binomial" : "multinomial"
        H2OGeneralizedLinearEstimator(;@NT(family, lambda_search, alpha)...)
    elseif name == "psvm"
        @unpack gamma, hyper_param = m
        H2OSupportVectorMachineEstimator(gamma = gamma, hyper_param = hyper_param)
    elseif name == "random_forest"
        H2ORandomForestEstimator()
    elseif name == "gbm"
        @unpack ntrees, sample_rate, col_sample_rate, col_sample_rate_per_tree, max_depth, min_rows = m
        H2OGradientBoostingEstimator(;@NT(ntrees, sample_rate, col_sample_rate, col_sample_rate_per_tree, max_depth, min_rows)...)
    elseif name == "xgboost"
        @unpack score_tree_interval, stopping_rounds, ntrees, learn_rate, col_sample_rate = m
        @unpack col_sample_rate_per_tree, max_depth, min_rows, sample_rate = m
        H2OXGBoostEstimator(;@NT(score_tree_interval, stopping_rounds, ntrees, learn_rate,
            col_sample_rate, col_sample_rate_per_tree, max_depth, min_rows, sample_rate)...)
    elseif name == "deeplearning"
        @unpack epochs, adaptive_rate, activation, rho, epsilon, input_dropout_ratio, hidden = m
        H2ODeepLearningEstimator(;@NT(epochs, adaptive_rate, activation, rho, epsilon, input_dropout_ratio)...)
    elseif name == "autoglm"
        @unpack max_models = m
        H2OAutoML(max_models = max_models, max_runtime_secs = 0, include_algos = ["GLM"])
    elseif name == "automl"
        @unpack max_models = m
        H2OAutoML(max_models = max_models, max_runtime_secs = 0, exclude_algos = ["XGBoost"])
    end
    # train
    weights_column = isnothing(w) ? nothing : "weight"
    pyo.train(x = columns, y = "label", weights_column = weights_column, training_frame = hdf)
    if occursin("auto", name)
        pyo, lb = pyo.leader, pyo.leaderboard
        @assert !isnothing(pyo)
        @redirect "summary.txt" begin
            println('*'^100, '\n')
            println("leader summary\n", '='^20)
            @trys println(pyo.summary().as_data_frame().to_string(), '\n')
            println("leaderboard table\n", '='^20)
            @trys println(lb.as_data_frame().to_string(), '\n')
            println("variable importance\n", '='^20)
            @trys println(pyo.varimp(true).to_string(), '\n')
            println('*'^100, '\n')
        end
    end
    h2o.save_model(pyo, force = true)
    # cleanup
    h2o.remove(hdf)
    h2o.api("POST /3/GarbageCollect")
    @pack! m = pyo
    BSON.bson("model.bson", model = m)
    return m
end

function predict_h2o(m::H2oModel, x)
    @unpack pyo = m
    columns = split(ENV["COLS"], '|')
    df = DataFrame(pymat(x), columns = columns)
    hdf = to_h2o(df)
    index = vec(hdf["index"].as_data_frame().values)
    pred = from_h2o(pyo.predict(hdf)).set_index(index).sort_index()
    @imports h2o
    h2o.remove(hdf)
    h2o.api("POST /3/GarbageCollect")
    return pred
end

predict_proba(m::H2oModel, x) = Array(predict_h2o(m, x).drop(columns = "predict"))

predict(m::H2oModel, x) = Array(predict_h2o(m, x)["predict"])

function gridparams(m::H2oModel)
    glm_comm = Dict("name" => "glm", "lambda_search" => true)
    glm_grid = Dict("alpha" => [0.0, 0.2, 0.4, 0.6, 0.8, 1.0])
    dl_comm = Dict("name" => "deeplearning", "epochs" => 100, "adaptive_rate" => true,
                    "activation" => "rectifier_with_dropout")
    dl_grid1 = Dict("rho" => [0.9, 0.95, 0.99], "epsilon" => [1e-6, 1e-7, 1e-8, 1e-9],
                    "input_dropout_ratio" => [0.0, 0.05, 0.1, 0.15, 0.2])
    dl_grid2 = Dict("hidden" => [[50], [200], [500]])
    dl_grid3 = Dict("hidden" => [[50, 50], [200, 200], [500, 500]], "hidden_dropout_ratios" =>
                    [[0.0, 0.0], [0.1, 0.1], [0.2, 0.2], [0.3, 0.3], [0.4, 0.4], [0.5, 0.5]])
    gbm_comm = Dict("name" => "gbm", "ntrees" => 100, "sample_rate" => 0.8,
                "col_sample_rate" => 0.8, "col_sample_rate_per_tree" => 0.8)
    gbm_grid = Dict("max_depth" => [6, 7, 8, 10, 15], "min_rows" => [1, 10, 10, 10, 100])
    xgb_comm = Dict("name" => "xgboost", "score_tree_interval" => 5, "stopping_rounds" => 5,
                        "ntrees" => 100, "learn_rate" => 0.05, "col_sample_rate" => 0.8,
                        "col_sample_rate_per_tree" => 0.8)
    xgb_grid = Dict("max_depth" => [10, 20, 5], "min_rows" => [5, 10, 3], "sample_rate" => [0.6, 0.6, 0.8])
    params = []
    append!(params, [merge(glm_comm, p) for p in gridparams(glm_grid, zip)])
    append!(params, [merge(gbm_comm, p) for p in gridparams(gbm_grid, zip)])
    append!(params, [merge(xgb_comm, p) for p in gridparams(xgb_grid, zip)])
    append!(params, [merge(dl_comm, p1, p2) for p1 in gridparams(xgb_grid, zip)
            for p2 in [gridparams(dl_grid2)..., gridparams(dl_grid3)...]])
    return params
end

modelhash(m::H2oModel) = hash(m.pyo)

function from_h2o(hdf)
    @from h2o imports export_file
    dst = tempname() * ".csv"
    for (c, t) in hdf.types
        if occursin(r"string|enum", t)
            hdf[c] = hdf[c].gsub("\"", "|")
        end
    end
    export_file(hdf, dst, force = true)
    df = pd.read_csv(dst)
    rm(dst, force = true)
    for (c, t) in hdf.types
        if t == "time"
            df[c] = pd.to_datetime(df[c], unit = "ms")
        end
    end
    return df
end

function to_h2o(df)
    @imports h2o
    dst = mkpath(joinpath("/dev/shm", tempname()))
    ncpu = h2o.cluster().nodes[1]["nthreads"]
    npart = ceil(Int, length(df) / ncpu)
    df["part"] = repeat(1:ncpu, inner = npart)[1:length(df)]
    df.reset_index(inplace = true)
    if Sys.ARCH == :x86_64
        @imports pyarrow as pa
        @imports pyarrow.parquet as pq
        table = pa.Table.from_pandas(df, preserve_index = false)
        pq.write_to_dataset(table, root_path = dst, partition_cols = ["part"])
    else
        @imports fastparquet as pq
        pq.write(dst, df, partition_on = ["part"], file_scheme = "hive")
        rm(joinpath(dst, "_metadata"), force = true)
        rm(joinpath(dst, "_common_metadata"), force = true)
    end
    hdf = h2o.import_file(dst)
    rm(dst, force = true, recursive = true)
    return hdf
end

H2oRegressor(;ka...) = H2oModel(;isclf = false, ka...)

H2oClassifier(;ka...) = H2oModel(;isclf = true, ka...)


text:
function harmonic(n::Integer)
    H = 0.0
    for i in n:-1:1
        H += 1 / i
    end
    return H
end

function harmonic(n::Integer, m::Integer)
    H = 0.0
    for i in n:-1:1
        H += 1 / (i ^ m)
    end
    return H
end

td_a1(n::Integer) = harmonic(n - 1)

td_a2(n::Integer) = harmonic(n - 1, 2)

td_b1(n::Integer) = (n + 1) / (3 * (n - 1))

function td_b2(n::Integer)
    return (2 * ((n ^ 2) + n + 3)) / ((9 * n) * (n - 1))
end

td_c1(a1::AbstractFloat, b1::AbstractFloat) = b1 - 1 / a1

function td_c2(n::Integer, a1::AbstractFloat, a2::AbstractFloat, b2::AbstractFloat)
    return b2 - ((n + 2) / (a1 * n)) + a2 / (a1 ^ 2)
end

td_e1(a1::AbstractFloat, c1::AbstractFloat) = c1 / a1

td_e2(a1::AbstractFloat, a2::AbstractFloat, c2::AbstractFloat) = c2 / (a1 ^ 2 + a2)

function tajimad(π::AbstractFloat, S::Integer, a1::AbstractFloat, e1::AbstractFloat, e2::AbstractFloat)
    return (π - S / a1) / sqrt(e1 * S + e2 * S * (S - 1))
end

"""
    tajimad(π::AbstractFloat, S::Integer, n::Integer)

Compute Tajima's D from:

* `π`: The average number of SNPs found in (n choose 2) pairwise comparisons of
       a sample of sequences.

* `S`: The number of segregating sites in a sample of sequences.

* `n`: The number of sequences in your sample.

*Example*

```julia
tajimad(3.88888, 16, 10)
```

"""
function tajimad(π::AbstractFloat, S::Integer, n::Integer)
    a1 = td_a1(n)
    a2 = td_a2(n)
    b1 = td_b1(n)
    b2 = td_b2(n)
    c1 = td_c1(a1, b1)
    c2 = td_c2(n, a1, a2, b2)
    e1 = td_e1(a1, c1)
    e2 = td_e2(a1, a2, c2)
    return tajimad(π, S, a1, e1, e2)
end

"""
    tajimad(seqs)

Compute Tajima's D from a collection of BioSequences{DNAAlphabet{n}} (n = 2 or 4).

This will estimate the `π`, `S`, and `n` parameters from the sequences and use
those parameters to estimate Tajima's D.

*Example*

```julia

sample = [dna"ATAATAAAAAAATAATAAAAAAATAAAAAAAATAAAAAAAA",
          dna"AAAAAAAATAAATAATAAAAAAATAAAAAAAAAAAAAAAAA",
          dna"AAAATAAAAATATAATAAAAAAATATAAAAAAAAAAAAAAA",
          dna"AAAAAAAAAAAATAATAAAAAAATAAATAAATAAAAAAAAA",
          dna"AAAATAAAAAAAATATAAAAAAATAAAAAAAAAAAAAAAAA",
          dna"AAAATAAAAAAAAAATAAAAAAAAAAAAAAAAAAATAAAAA",
          dna"AAAAAATAAAAATAATAAAAAAATAAAAAAAAAAAAAAAAA",
          dna"AAAAAAAAAAAAAAATAAAAAAATAAAAAAAAAAAAAAATA",
          dna"AAAAAAAAAAAAAAAAAAAAAAATAAAAAAAAAAAAAAAAA",
          dna"AAAAAAAAAAAAAAATAAAAAAATAATAAAAAAAAAAAAAA"]

tajimad(sample)
```

"""
function tajimad(seqs)
    π = avg_mut(seqs)
    S = count(Segregating, seqs)
    n = length(seqs)
    return tajimad(π, S[1], n)
end


text:
using Zygote, StaticArrays

# Naive Matrix Multiply

function mul(a::AbstractMatrix, b::AbstractMatrix)
  [sum(a[i, k]*b[k,j] for k = 1:size(a, 2))
   for i = 1:size(a, 1), j = 1:size(b, 2)]
end

a = [1 2; 3 4]

mul(a, a)

# Static Matrix Multiple

_sum(xs) = reduce((a, b) -> :($a+$b), xs)

@generated function mul(a::SMatrix{M,N}, b::SMatrix{N,O}) where {M,N,O}
  :(SMatrix{M,O}(
      $([_sum(:(a[$i, $k]*b[$k,$j]) for k = 1:size(a, 2))
         for i = 1:size(a, 1), j = 1:size(b, 2)]...),))
end

@generated function mul(a::SMatrix{M,N}, b::SMatrix{N,O}) where {M,N,O}
  quote
    Base.@_inline_meta
    SMatrix{M,O}(
      $([_sum(:(a[$i, $k]*b[$k,$j]) for k = 1:size(a, 2))
        for i = 1:size(a, 1), j = 1:size(b, 2)]...),)
  end
end

a = @SArray [1 2; 3 4]

mul(a, a)

sqmul(a) = sum(mul(a, a))

@code_llvm sqmul(a)

# See the code output

function _mul(a::SMatrix{M,N}, b::SMatrix{N,O}) where {M,N,O}
  :(SMatrix{M,O}(
      $([_sum(:(a[$i, $k]*b[$k,$j]) for k = 1:size(a, 2))
         for i = 1:size(a, 1), j = 1:size(b, 2)]...),))
end

_mul(a, a)

max_stars_repo_name: MikeInnes/zygote-paper

max_issues_repo_name: MikeInnes/zygote-paper

max_forks_repo_name: MikeInnes/zygote-paper


text:
# This only works if v is a n*1 matrix and the op is *.
function Base.broadcast{Tv,Ti}(op, v::Array{Tv,2}, A::SparseMatrixCSC{Tv,Ti})
  I, J = findn(A)
  V = zeros(nnz(A))
  vn, vm = size(v)
  if vn >= 1 && vm == 1
    for(l in 1:nnz(A))
      row = I[l]
      V[l] = op(v[row], A.nzval[l])
    end
  elseif vn == 1 && vm >= 1
    for(l in 1:nnz(A))
      col = J[l]
      V[l] = op(v[col], A.nzval[l])
    end
  else
    throw(ArgumentError("invalid dimensions"))
  end
  sparse(I, J, V)
end

n = 100
m = 42

v = randn(n,1)
s = sprandn(n,m,0.75)
@assert full(broadcast(*, v, s)) == broadcast(*, v, full(s))

v = randn(1,n)
s = sprandn(m,n,0.42)
res = broadcast(*, v, s)
@assert full(broadcast(*, v, s)) == broadcast(*, v, full(s))


text:
## low-level pcre interface ##

module PCRE

include("pcre_h.jl")

const VERSION = bytestring(ccall((:pcre_version, :libpcre), Ptr{Uint8}, ()))

# supported options for different use cases

const COMPILE_MASK      =
      ANCHORED          |
      CASELESS          |
      DOLLAR_ENDONLY    |
      DOTALL            |
      EXTENDED          |
      FIRSTLINE         |
      JAVASCRIPT_COMPAT |
      MULTILINE         |
      NEWLINE_ANY       |
      NEWLINE_ANYCRLF   |
      NEWLINE_CR        |
      NEWLINE_CRLF      |
      NEWLINE_LF        |
      NO_AUTO_CAPTURE   |
      NO_START_OPTIMIZE |
      NO_UTF8_CHECK     |
      UNGREEDY          |
      UTF8

const EXECUTE_MASK      =
      NEWLINE_ANY       |
      NEWLINE_ANYCRLF   |
      NEWLINE_CR        |
      NEWLINE_CRLF      |
      NEWLINE_LF        |
      NOTBOL            |
      NOTEMPTY          |
      NOTEMPTY_ATSTART  |
      NOTEOL            |
      NO_START_OPTIMIZE |
      NO_UTF8_CHECK     |
      PARTIAL_HARD      |
      PARTIAL_SOFT

const OPTIONS_MASK = COMPILE_MASK | EXECUTE_MASK

function info{T}(
    regex::Ptr{Void},
    extra::Ptr{Void}, what::Integer, ::Type{T}
)
    buf = Array(Uint8,sizeof(T))
    ret = ccall((:pcre_fullinfo, :libpcre), Int32,
                (Ptr{Void}, Ptr{Void}, Int32, Ptr{Uint8}),
                regex, extra, what, buf)
    if ret != 0
        error(ret == ERROR_NULL      ? "NULL regex object" :
              ret == ERROR_BADMAGIC  ? "invalid regex object" :
              ret == ERROR_BADOPTION ? "invalid option flags" :
                                       "unknown error $ret")
    end
    reinterpret(T,buf)[1]
end

function config{T}(what::Integer, ::Type{T})
    buf = Array(Uint8, sizeof(T))
    ret = ccall((:pcre_config, :libpcre), Int32,
                (Int32, Ptr{Uint8}),
                what, buf)

    if ret != 0
        error("error $n")
    end
    reinterpret(T,buf)[1]
end

function compile(pattern::String, options::Integer)
    errstr = Array(Ptr{Uint8},1)
    erroff = Array(Int32,1)
    re_ptr = ccall((:pcre_compile, :libpcre), Ptr{Void},
                    (Ptr{Uint8}, Int32, Ptr{Ptr{Uint8}}, Ptr{Int32}, Ptr{Uint8}),
                    pattern, options, errstr, erroff, C_NULL)
    if re_ptr == C_NULL
        error("$(bytestring(errstr[1]))",
              " at position $(erroff[1]+1)",
              " in $(repr(pattern))")
    end
    re_ptr
end

function study(regex::Ptr{Void}, options::Integer)
    # NOTE: options should always be zero in current PCRE
    errstr = Array(Ptr{Uint8},1)
    extra = ccall((:pcre_study, :libpcre), Ptr{Void},
                  (Ptr{Void}, Int32, Ptr{Ptr{Uint8}}),
                  regex, options, errstr)
    if errstr[1] != C_NULL
        error("$(bytestring(errstr[1]))")
    end
    extra
end
study(re::Ptr{Void}) = study(re, int32(0))

free_study(extra::Ptr{Void}) =
    ccall((:pcre_free_study, :libpcre), Void, (Ptr{Void},), extra)
free(regex::Ptr{Void}) =
    ccall(unsafe_load(cglobal((:pcre_free, :libpcre),Ptr{Void})), Void, (Ptr{Void},), regex)

exec(regex::Ptr{Void}, extra::Ptr{Void}, str::SubString, offset::Integer, options::Integer, cap::Bool) =
    exec(regex, extra, str.string, str.offset, offset, sizeof(str), options, cap)
exec(regex::Ptr{Void}, extra::Ptr{Void}, str::ByteString, offset::Integer, options::Integer, cap::Bool) =
    exec(regex, extra, str, 0, offset, sizeof(str), options, cap)
function exec(regex::Ptr{Void}, extra::Ptr{Void},
              str::ByteString, shift::Integer, offset::Integer, len::Integer, options::Integer, cap::Bool)
    if offset < 0 || len < offset || len+shift > sizeof(str)
        error(BoundsError)
    end
    ncap = info(regex, extra, INFO_CAPTURECOUNT, Int32)
    ovec = Array(Int32, 3(ncap+1))
    n = ccall((:pcre_exec, :libpcre), Int32,
              (Ptr{Void}, Ptr{Void}, Ptr{Uint8}, Int32,
               Int32, Int32, Ptr{Int32}, Int32),
              regex, extra, pointer(str.data,shift+1), len,
              offset, options, ovec, length(ovec))
    if n < -1
        error("error $n")
    end
    cap ? ((n > -1 ? ovec[1:2(ncap+1)] : Array(Int32,0)), ncap) : n > -1
end

end # module

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
############
# Indexing #
############

@inline function compute_index{dim}(::Type{SymmetricTensor{2, dim}}, i::Int, j::Int)
    if i < j
        i, j = j, i
    end
    # We are skipping triangle over the diagonal = (j-1) * j / 2 indices
    skipped_indicies = div((j-1) * j, 2)
    return dim*(j-1) + i - skipped_indicies
end

@inline function compute_index{dim}(::Type{Tensor{2, dim}}, i::Int, j::Int)
    return dim*(j-1) + i
end

@inline function compute_index{dim}(T::Type{Tensor{4, dim}}, i::Int, j::Int, k::Int, l::Int)
    lower_order = Tensor{2, dim}
    I = compute_index(lower_order, i, j)
    J = compute_index(lower_order, k, l)
    n = n_components(lower_order)
    return (J-1) * n + I
end

@inline function compute_index{dim}(T::Type{SymmetricTensor{4, dim}}, i::Int, j::Int, k::Int, l::Int)
    lower_order = SymmetricTensor{2, dim}
    I = compute_index(lower_order, i, j)
    J = compute_index(lower_order, k, l)
    n = n_components(lower_order)
    return (J-1) * n + I
end

###########################
# getindex general tensor #
###########################
@inline function Base.getindex(S::Tensor, i::Int)
    @boundscheck checkbounds(S, i)
    @inbounds v = get_data(S)[i]
    return v
end

@inline function Base.getindex{dim}(S::SymmetricTensor{2, dim}, i::Int, j::Int)
    @boundscheck checkbounds(S, i, j)
    @inbounds v = get_data(S)[compute_index(SymmetricTensor{2, dim}, i, j)]
    return v
end

@inline function Base.getindex{dim}(S::SymmetricTensor{4, dim}, i::Int, j::Int, k::Int, l::Int)
    @boundscheck checkbounds(S, i, j, k, l)
    @inbounds v = get_data(S)[compute_index(SymmetricTensor{4, dim}, i, j, k, l)]
    return v
end

# Slice
@inline Base.getindex(v::Vec, ::Colon) = v

function Base.getindex(S::Union{SecondOrderTensor, FourthOrderTensor}, ::Colon)
    throw(ArgumentError("S[:] not defined for S of order 2 or 4, use Array(S) to convert to an Array"))
end

@inline @generated function Base.getindex{dim, T}(S::SecondOrderTensor{dim, T}, ::Colon, j::Int)
    idx2(i,j) = compute_index(get_base(S), i, j)
    ex1 = Expr(:tuple, [:(get_data(S)[$(idx2(i,1))]) for i in 1:dim]...)
    ex2 = Expr(:tuple, [:(get_data(S)[$(idx2(i,2))]) for i in 1:dim]...)
    ex3 = Expr(:tuple, [:(get_data(S)[$(idx2(i,3))]) for i in 1:dim]...)
    return quote
        @boundscheck checkbounds(S,Colon(),j)
        if     j == 1 return Vec{dim}($ex1)
        elseif j == 2 return Vec{dim}($ex2)
        else          return Vec{dim}($ex3)
        end
    end
end
@inline @generated function Base.getindex{dim, T}(S::SecondOrderTensor{dim, T}, i::Int, ::Colon)
    idx2(i,j) = compute_index(get_base(S), i, j)
    ex1 = Expr(:tuple, [:(get_data(S)[$(idx2(1,j))]) for j in 1:dim]...)
    ex2 = Expr(:tuple, [:(get_data(S)[$(idx2(2,j))]) for j in 1:dim]...)
    ex3 = Expr(:tuple, [:(get_data(S)[$(idx2(3,j))]) for j in 1:dim]...)
    return quote
        @boundscheck checkbounds(S,i,Colon())
        if     i == 1 return Vec{dim}($ex1)
        elseif i == 2 return Vec{dim}($ex2)
        else          return Vec{dim}($ex3)
        end
    end
end


text:
########################################################################################################
## Details
# This file will create the structure for the UQ inference problems.
# This package is designed to be used with a family of julia packages,
#           - Potential.jl (*)
#               - Creates potential types, has abstract 'Potential' type (*), let p <: 'Potential'
#               - Assume mutable struct Potential{T<:AbstractFloat}
#                            x       # Potential trainable parameters
#                            p       # Potential nontrainable parameters
#                        end
#               - Assume mutable struct SNAP{T<:AbstractFloat}
#                            A
#                            b
#                            β
#                        end
#           - PotentialLearning.jl
#               - Fits data-driven potential
#               - We assume that potential learning updates parameters p.x
#           - PotentialUQ.jl (this package)
#               - Defines distribution over potential parameters.
#               - Initially must provide distribution over p.x (assume user provides distribution, in form of log likelihood)
#               - In case of SNAP, potential is predetermined.
#               - Allow users to modify distribution with priors by adding to log-likelihood.
#
#########################################################################################################
module PotentialUQ

using Base: NamedTuple
using LinearAlgebra: tr, logdet, transpose, Symmetric, diag, I, pinv, diagm, svd, dot
using Statistics: cov, std, mean
using TransformVariables, Bijectors
using Distributions
using Random
import InteratomicPotentials as Potentials

include("types/types.jl")
include("samplers/samplers.jl")
include("utilities/utils.jl")
# include("optimizers/map.jl")

export state, chain, chain_diagnostics, mcmc_step, mcmc, hmc_step, hmc, Sample

end


text:
# ## Function plot
# ### simple function

#md # ![FILE_NAME.png](images/FILE_NAME.png)

f(x) = 5exp(-x^2)
g(x) = x^2

plot([f, g], -3, 3)


text:
function  removeedge(C, C_age)
global amax
row, col = findnz(C_age .> amax);
a = size(row,2);
if !isempty(row)
    for i = 1:a
        C_age[row[i],col[i]] = 0;
        C_age[col[i],row[i]] = 0;
        C[row[i],col[i]] = 0;
        C[col[i],row[i]] = 0;
    end
end
return C, C_age
end


text:
"""
Julia port of MPF's SPGL1
"""
module GenSPGL

# Modules
using Printf
using Nullables
using LinearAlgebra
using MAT
using JLD
using Arpack
using JOLI

# Types
include("spgOptions.jl")
include("spgExitCondition.jl")
include("spgInit.jl")
include("spgInfo.jl")

# Methods
include("spgl1.jl")
include("spglcore.jl")
include("spg_lasso.jl")
include("NormL1_project.jl")
include("NormL1_dual.jl")
include("NormL1_primal.jl")
include("oneprojector.jl")
include("spgline.jl")
include("spglinecurvy.jl")
include("NLfunForward.jl")
include("TraceNorm_project.jl")
include("TraceNorm_dual.jl")
include("TraceNorm_primal.jl")
include("afun.jl")

# Penalties
include("Penalty/funLS.jl")

# Composite
include("Composite/funCompR1.jl")
include("Composite/funCompR2.jl")

# Examples
#include("Examples/spgdemo.jl")
include("Examples/compare/jl_cs.jl")
include("Examples/pres/ex_cs.jl")
include("Examples/scaletest/scale_cs.jl")
include("Examples/compare/jl_complex.jl")

end # module


text:
################################################################################################
#
# demo 5: set video properties
#
################################################################################################

# Julia variables
winname = "Welcome"
fname = joinpath(Pkg.dir("OpenCV"), "./test/images/movie.avi")

# C++ OpenCV code
cxx"""
void setVideoProperties(const char *winname, const char *fname)
{
    // convert const char* => string
    std::string wname = winname;
    std::string fvideo = fname;

    // Open video file
    cv::VideoCapture capture(fvideo);

    if (!capture.isOpened())
    {
        std::cout << "Cannot open the video file!" << std::endl;
        exit(0);
    }

//    Properties to control:
//    int propID
//    cv::CAP_PROP_POS_MSEC - current position of the video in milliseconds
//    cv::CAP_PROP_POS_FRAMES - current position of the video in frames
//    cv::CAP_PROP_FRAME_WIDTH - width of the frame of the video stream
//    cv::CAP_PROP_FRAME_HEIGHT - height of the frame of the video stream
//    cv::CAP_PROP_FPS - frame rate (frames per second)
//    cv::CAP_PROP_FOURCC - four character code  of codec

    capture.set(cv::CAP_PROP_POS_MSEC, 300); //300 ms into the video
    capture.set(cv::CAP_PROP_FRAME_WIDTH, 600);
    capture.set(cv::CAP_PROP_FRAME_HEIGHT, 400);
    capture.set(cv::CAP_PROP_FPS, 5);

    // Get the fps
    double fps = capture.get(cv::CAP_PROP_FPS);

    std::cout << "Frame per second: "<< fps << std::endl;
    cv::namedWindow("Welcome", cv::WINDOW_AUTOSIZE);

    cv::Mat frame;

    while(true)
    {
        bool bSuccess = capture.read(frame);
        if (!bSuccess)
        {
             std::cout << "Can not read the video frame!"<< std::endl;
             break;
        }

        cv::imshow("Welcome", frame);

        if (cv::waitKey(30)== 27)  // wait for esc key(27) press for (30) msec
        {
            std::cout << "Pressed esc to stop."<< std::endl;
            cv::destroyWindow("Welcome");
            break;
        }
    }
}
"""

@cxx setVideoProperties(pointer(winname), pointer(fname))


text:
push!(LOAD_PATH,"../")

using CubicEoSDatabase
using Documenter

makedocs(;
    sitename="CubicEoSDatabase.jl",
    modules=[CubicEoSDatabase],
    pages=[
        "index.md",
        "outofboxdb.md",
        "apireference.md",
    ],
)

deploydocs(;
    repo = "github.com/stepanzh/CubicEoSDatabase.jl.git",
    devbranch = "main",
)


text:
using ProfileView
using PolynomialBasis
using ImplicitDomainQuadrature
using BenchmarkTools
# using Revise
using CutCellDG
include("../../../test/useful_routines.jl")
include("test-bounds.jl")

interfacecenter = [1.5, 0.5]
interfaceradius = 1.0

levelsetbasis = HermiteTensorProductBasis(2)
interpgrad = InterpolatingPolynomial(2, 2, 3)
numqp = required_quadrature_order(3) + 2
dim, numpts = size(interpolation_points(levelsetbasis))
quad = tensor_product_quadrature(2, 4)
cgmesh = CutCellDG.CGMesh([0.0, 0.0], [1.0, 1.0], [1, 1], numpts)

levelset = CutCellDG.LevelSet(
    x -> -circle_distance_function(x, interfacecenter, interfaceradius)[1],
    cgmesh,
    levelsetbasis,
    quad,
)

CutCellDG.load_coefficients!(levelset, 1)
update_interpolating_gradient!(interpgrad, CutCellDG.interpolater(levelset))

function run_area_quadrature(levelset, interpgrad, numqp; numiter = 100)
    for i = 1:numiter
        pquad = area_quadrature(
            CutCellDG.interpolater(levelset),
            interpgrad,
            +1,
            [-1.0, -1.0],
            [1.0, 1.0],
            numqp,
        )
    end
end

# Run once to compile
# run_area_quadrature(levelset,interpgrad,numqp)

# @profview run_area_quadrature(levelset, interpgrad, numqp)

using IntervalArithmetic
xL,xR = [-1.,-1.],[1.,1.]
box = IntervalBox(xL,xR)
# @btime ImplicitDomainQuadrature.is_suitable(1,interpgrad,box)

curve(x) = ImplicitDomainQuadrature.curvature_measure(interpgrad(x),1,4)
# @btime sign(curve,box,tol=1e-3)
interval_arithmetic_sign_search(curve,box,1e-3,0.0,2)

@btime area_quadrature(CutCellDG.interpolater(levelset),interpgrad)


text:
using  WebSockets
import WebSockets:Response, Request

using JSON
using Sockets
using Dates
using Random

include("game/abstract-game.jl")
include("game/yes-please-logic.jl")
include("game/yes-please-server.jl")

include("server/output.jl")
include("server/spectation.jl")
include("server/room.jl")
include("server/utils.jl")
include("server/disconnection.jl")

const LOCALIP = string(Sockets.getipaddr())
const HTTPPORT = 8080

getFile(str) = read(joinpath(@__DIR__,str), String)

const validSuffix = Dict("/"              => getFile("web/index.html"),
                         "/index.html"    => getFile("web/index.html"),
                         "/style.css"     => getFile("web/style.css"),
                         "/comms.js"      => getFile("web/comms.js"),
                         "/myprompts.js"  => getFile("web/myprompts.js"),
                         "/login.js"      => getFile("web/login.js"),
                         "/room.js"       => getFile("web/room.js"),
                         "/winstack.js"   => getFile("web/winstack.js"))

function HTMLHandler(req::Request)
    if req.method == "GET" && haskey(validSuffix, req.target)
        return validSuffix[req.target] |> Response
    end
    "" |> Response ## 404 - Page Not Found Message
end

function WSGatekeeper(req, ws)
    orig = WebSockets.origin(req)
    # To be replaced with External IP on deployment
    if occursin(LOCALIP, orig)
        # Run Main Coroutine
        coroutine(ws)
    else
        @warn("Unauthorized websocket connection, $orig not approved by gatekeeper, expected $LOCALIP")
    end
    nothing
end

function coroutine(currws)
    while isopen(currws)
        data, success = readguarded(currws)
        !success && break
        textmsg = String(data)
        println("Recieved: $textmsg")
        msg = JSON.parse(textmsg)
        !validJSON(msg, "querytype") && break
        if msg["querytype"] == "join-room"
            !validJSON(msg, "particulars") && break
            joinroom(currws, msg["particulars"])
        elseif msg["querytype"] == "make-room"
            !validJSON(msg, "particulars") && break
            makeroom(currws, msg["particulars"])
        elseif msg["querytype"] == "leave-room"
            leaveroom(currws)
        elseif msg["querytype"] == "change-state-spectate"
            changestatespec(currws)
        elseif msg["querytype"] == "change-state-play"
            changestateplay(currws)

        ######
        # Change this line to modify the server game
        ######
        elseif msg["querytype"] in YPCommandsList
            YPServer(msg, currws)
        end
    end

    # Clean up External References
    removereferences(currws)
    nothing
end

global SERVER = WebSockets.ServerWS(HTMLHandler, WSGatekeeper)
@async WebSockets.serve(SERVER, LOCALIP, HTTPPORT)
println("HTTP Server listening on $LOCALIP:$HTTPPORT.")


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct CleanEnergyReserve{T <: ReserveDirection} <: Reserve{T}
        name::String
        available::Bool
        time_frame::Float64
        requirement::Float64
        ext::Dict{String, Any}
        time_series_container::InfrastructureSystems.TimeSeriesContainer
        internal::InfrastructureSystemsInternal
    end

Data Structure for the procurement products for system simulations.

# Arguments
- `name::String`
- `available::Bool`
- `time_frame::Float64`: the relative saturation time_frame, validation range: `(0, nothing)`, action if invalid: `error`
- `requirement::Float64`: the required quantity of the product should be scaled by a TimeSeriesData
- `ext::Dict{String, Any}`
- `time_series_container::InfrastructureSystems.TimeSeriesContainer`: internal time_series storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct CleanEnergyReserve{T <: ReserveDirection} <: Reserve{T}
    name::String
    available::Bool
    "the relative saturation time_frame"
    time_frame::Float64
    "the required quantity of the product should be scaled by a TimeSeriesData"
    requirement::Float64
    ext::Dict{String, Any}
    "internal time_series storage"
    time_series_container::InfrastructureSystems.TimeSeriesContainer
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function CleanEnergyReserve{T}(name, available, time_frame, requirement, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), ) where T <: ReserveDirection
    CleanEnergyReserve{T}(name, available, time_frame, requirement, ext, time_series_container, InfrastructureSystemsInternal(), )
end

function CleanEnergyReserve{T}(; name, available, time_frame, requirement, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), internal=InfrastructureSystemsInternal(), ) where T <: ReserveDirection
    CleanEnergyReserve{T}(name, available, time_frame, requirement, ext, time_series_container, internal, )
end

# Constructor for demo purposes; non-functional.
function CleanEnergyReserve{T}(::Nothing) where T <: ReserveDirection
    CleanEnergyReserve{T}(;
        name="init",
        available=false,
        time_frame=0.0,
        requirement=0.0,
        ext=Dict{String, Any}(),
        time_series_container=InfrastructureSystems.TimeSeriesContainer(),
    )
end

"""Get [`CleanEnergyReserve`](@ref) `name`."""
PSY.get_name(value::CleanEnergyReserve) = value.name
"""Get [`CleanEnergyReserve`](@ref) `available`."""
PSY.get_available(value::CleanEnergyReserve) = value.available
"""Get [`CleanEnergyReserve`](@ref) `time_frame`."""
PSY.get_time_frame(value::CleanEnergyReserve) = value.time_frame
"""Get [`CleanEnergyReserve`](@ref) `requirement`."""
PSY.get_requirement(value::CleanEnergyReserve) = value.requirement
"""Get [`CleanEnergyReserve`](@ref) `ext`."""
PSY.get_ext(value::CleanEnergyReserve) = value.ext
"""Get [`CleanEnergyReserve`](@ref) `time_series_container`."""
PSY.get_time_series_container(value::CleanEnergyReserve) = value.time_series_container
"""Get [`CleanEnergyReserve`](@ref) `internal`."""
PSY.get_internal(value::CleanEnergyReserve) = value.internal

"""Set [`CleanEnergyReserve`](@ref) `name`."""
PSY.set_name!(value::CleanEnergyReserve, val) = value.name = val
"""Set [`CleanEnergyReserve`](@ref) `available`."""
PSY.set_available!(value::CleanEnergyReserve, val) = value.available = val
"""Set [`CleanEnergyReserve`](@ref) `time_frame`."""
PSY.set_time_frame!(value::CleanEnergyReserve, val) = value.time_frame = val
"""Set [`CleanEnergyReserve`](@ref) `requirement`."""
PSY.set_requirement!(value::CleanEnergyReserve, val) = value.requirement = val
"""Set [`CleanEnergyReserve`](@ref) `ext`."""
PSY.set_ext!(value::CleanEnergyReserve, val) = value.ext = val
"""Set [`CleanEnergyReserve`](@ref) `time_series_container`."""
PSY.set_time_series_container!(value::CleanEnergyReserve, val) = value.time_series_container = val

max_stars_repo_name: sourabhdalvi/PowerSystemExtensions

max_issues_repo_name: sourabhdalvi/PowerSystemExtensions

max_forks_repo_name: sourabhdalvi/PowerSystemExtensions


text:
struct SimpleFunctionMap end
export SimpleFunctionMap

# ConstantCache version
function DiffEqBase.__solve(prob::DiffEqBase.DiscreteProblem{uType,tupType,false},
                          alg::SimpleFunctionMap;
                          calculate_values=true) where {uType,tupType}
    tType = eltype(tupType)
    tspan = prob.tspan
    f = prob.f
    u0 = prob.u0
    p = prob.p
    dt = 1

    n = Int((tspan[2] - tspan[1])/ dt + 1)
    u = [u0 for i in 1:n]
    t = tspan[1]:dt:tspan[2]
    if calculate_values
        for i in 2:n
            u[i] = f(u[i-1],p,t[i])
        end
    end
    sol = DiffEqBase.build_solution(prob,alg,t,u,dense=false,
                         interp = DiffEqBase.ConstantInterpolation(t,u),
                         calculate_error = false)
end

# Cache version
function DiffEqBase.__solve(prob::DiscreteProblem{uType,tupType,true},
                          alg::SimpleFunctionMap;
                          calculate_values=true) where {uType,tupType}
    tType = eltype(tupType)
    tspan = prob.tspan
    f = prob.f
    u0 = prob.u0
    p = prob.p
    dt = 1

    n = Int((tspan[2] - tspan[1])/ dt + 1)
    u = [similar(u0) for i in 1:n]
    u[1] .= u0
    t = tspan[1]:dt:tspan[2]
    if calculate_values
        for i in 2:n
            f(u[i],u[i-1],p,t[i])
        end
    end
    sol = DiffEqBase.build_solution(prob,alg,t,u,dense=false,
                         interp = DiffEqBase.ConstantInterpolation(t,u),
                         calculate_error = false)
end

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

# Integrator version

mutable struct DiscreteIntegrator{F,uType,tType,P,S} <: DiffEqBase.DEIntegrator
    f::F
    u::uType
    t::tType
    uprev::uType
    p::P
    sol::S
    i::Int
end

function DiffEqBase.__init(prob::DiscreteProblem,
                         alg::SimpleFunctionMap)
    sol = solve(prob,alg;calculate_values=false)
    DiscreteIntegrator(prob.f,prob.u0,prob.tspan[1],copy(prob.u0),prob.p,sol,1)
end

function step!(integrator::DiscreteIntegrator)
    integrator.t = integrator.i
    integrator.i += 1
    u = integrator.u
    uprev = integrator.uprev
    p = integrator.p
    f = integrator.f
    i = integrator.i

    if isinplace(integrator.sol.prob)
        f(integrator.sol.u[i],uprev,p,i)
        integrator.uprev = integrator.u
        integrator.u = integrator.sol.u[i]
    else
        u = f(uprev,p,i)
        integrator.sol.u[i] = u
        integrator.uprev = integrator.u
        integrator.u = u
    end
end


text:
"""
    error_occurred(log::String)

Return `true` if an error occurred.
It would be more stable if Weave would have a fail on error option or something similar.
"""
function error_occurred(log::String)
    weave_error = contains(log, "ERROR")
end

const WEAVE_LOG_FILE = "weave.log"

log_path(folder) = joinpath(REPO_DIR, "tutorials", folder, WEAVE_LOG_FILE)

folder2filename(folder) = replace(folder, '-' => '_'; count=1)

"""
    markdown_output(folder)

Returns the Markdown output for a folder.
The output seems to be the only place where Weave prints the full stacktrace.
"""
function markdown_output(folder)
    file = folder2filename(folder)
    file = "$file.md"
    path = joinpath(REPO_DIR, "markdown", folder, file)
    text = read(path, String)
    """
    Markdown output (contains stacktrace):
    $text
    """
end

"""
    build_folder(folder)

It seems that Weave has no option to fail on error, so we handle errors ourselves.
Also, this method only shows the necessary information in the CI logs.
If something crashes, then show the logs immediately.
If all goes well, then store the logs in a file, but don't show them.
"""
function build_folder(folder)
    println("$folder - Starting build")
    cache = :all
    c = IOCapture.capture() do
        @timed weave_folder(folder; cache)
    end
    stats = c.value
    gib = round(stats.bytes / 1024^3, digits=2)
    min = round(stats.time / 60, digits=2)
    println("$folder - Build took $min minutes and allocated $gib GiB:")
    log = c.output
    md_out = markdown_output(folder)
    if error_occurred(log)
        @error """
        $folder - Error occured:
        $log

        $md_out
        """
    end
    path = log_path(folder)
    println("$folder - Writing log to $path")
    write(path, log)
end

"""
    safe_instantiate(folders)

Install all packages sequentially since Pkg.jl is thread-unsafe.
See https://github.com/JuliaLang/Pkg.jl/issues/2219 for details.
"""
function safe_instantiate(folders)
    script = "import Pkg; Pkg.activate(ARGS[1]); Pkg.instantiate()"
    for folder in folders
        folder == "99-test" && return nothing
        dir = tutorial_path(folder)
        @info "Instantiating project environment in $folder"
        cmd = `$(Base.julia_cmd()) -e $script $dir`
        p = run(cmd)
        if !success(p)
            error("Couldn't instantiate project environment of $folder")
        end
    end
    return nothing
end

"""
    parallel_build(folders)

Build `folders` in parallel inside new Julia instances.
This has two benefits, namely that it ensures that globals are reset and reduces the
running time.
"""
function parallel_build(folders)
    safe_instantiate(folders)
    # The static schedule creates one task per thread.
    Threads.@threads :static for folder in folders
        ex = """using TuringTutorials; build_folder("$folder")"""
        cmd = `$(Base.julia_cmd()) --project -e $ex`
        run(cmd)
    end
end

function log_has_error(folder)::Bool
    path = log_path(folder)
    if isfile(path)
        println("$folder - Verifying the log")
        log = read(path, String)
        has_error = error_occurred(log)
        println("""$folder - Log contains $(has_error ? "an" : "no") error""")
        return has_error
    else
        println("$folder - No file found to verify")
        return false
    end
end

"""
    verify_logs(T::Vector)::Bool

Return `true` if logs for the tutorials `T` contain an error.
This method is used at the end of the CI in order to allow the CI to fail only after
running all the tutorials (similar to `Pkg.test()`).
"""
verify_logs(T::Vector)::Bool = !any(log_has_error.(T))

"""
    build(T::Vector=changed_tutorials())::Bool

Build all changed outputs.
For example, pass `tutorials()` to build all tutorials or `["00-introduction"]` to build
only the first.
"""
function build(T::Vector=changed_tutorials())::Bool
    "CI" in keys(ENV) && download_artifacts()
    clean_weave_cache()
    parallel_build(T)
    out = verify_logs(T)
    # Avoid committing cache files to the artifacts branch.
    clean_weave_cache()
    return out
end
build(tutorial::AbstractString)::Bool = build([tutorial])

"""
    build_and_exit(T::Union{Vector,AbstractString})

Build tutorial(s) `T` and exit with 1 if an error occurred during build.
This method is used in the CI job.
"""
function build_and_exit(T)
    success = build(T)
    if !success
        println("One of the logs contains an error. Exiting with `exit(1)`")
    end
    code = success ? 0 : 1
    exit(code)
end

max_stars_repo_name: staticfloat/TuringTutorials

max_issues_repo_name: staticfloat/TuringTutorials

max_forks_repo_name: staticfloat/TuringTutorials


text:
using Liblet
using Test

test_derivation_repr() = begin
    G = Grammar("""
        S -> A B
        A -> a
        B -> b
    """, false)
    d = Derivation(G)
    for (prod, pos) ∈ [(1, 1), (2, 1), (3, 2)]
        d = next(d, prod, pos)
    end
    d.repr == "S -> A B -> a B -> a b"
end

test_derivation_eq() = begin
    G = Grammar("""
            S -> A B
            A -> a
            B -> b
        """, false)
    d0 = Derivation(G)
    for (prod, pos) ∈ [(1, 1), (2, 1), (3, 2)]
        d0 = next(d0, prod, pos)
    end
    d1 = Derivation(G)
    for (prod, pos) ∈ [(1, 1), (2, 1), (3, 2)]
        d1 = next(d1, prod, pos)
    end
    d0 == d1
end

test_derivation_hash() = begin
    G = Grammar("""
            S -> A B
            A -> a
            B -> b
        """, false)
    d0 = Derivation(G)
    for (prod, pos) ∈ [(1, 1), (2, 1), (3, 2)]
        d0 = next(d0, prod, pos)
    end
    d1 = Derivation(G)
    for (prod, pos) ∈ [(1, 1), (2, 1), (3, 2)]
        d1 = next(d1, prod, pos)
    end
    S = Dict( d0 => 0, d1 => 1)
    length(S) == 1
end

test_derivation_steps() = begin
    G = Grammar("""
            S -> A B
            A -> a
            B -> b
        """, false)
    d = Derivation(G)
    steps = [(1, 1), (2, 1), (3, 2)]
    for (prod, pos) ∈ steps
        d = next(d, prod, pos)
    end
    steps == d.steps
end

test_derivation_byprods() = begin
    G = Grammar("""
            S -> A B
            A -> a
            B -> b
        """)
    d = Derivation(G)
    p = [Production("S", ["A", "B"]), Production("A", ["a"])]
    leftmost(d, p) |> sententialform == ["a", "B"]
end

test_derivation_byprod() = begin
    G = Grammar("""
        S -> A B
        A -> a
        B -> b
    """)
    d = Derivation(G)
    p = Production("S", ["A", "B"])
    leftmost(d, p) |> sententialform == ["A", "B"]
end

test_derivation_steps_list() = begin
    G = Grammar("""
        S -> A B
        A -> a
        B -> b
    """, false)
    k = [(1, 1), (2, 1), (3, 2)]
    d = Derivation(G)
    d = next(d,k)
    k == steps(d)
end

test_derivation_wrong_step() = begin
    G = Grammar("""
        S -> A B
        A -> a
        B -> b
    """, false)
    d = Derivation(G)
    steps = [(1, 1), (2, 2)]
    for (prod, pos) ∈ steps
        d = next(d, prod, pos)
    end
end

test_derivation_leftmost() = begin
    G = Grammar("""
        S -> A B
        A -> a
        B -> b
    """)
    d = Derivation(G) |>
        (x -> leftmost(x, 1)) |>
        (x -> leftmost(x, 2)) |>
        (x -> leftmost(x, 3))
    s = [(1, 1), (2, 1), (3, 2)]
    s == steps(d)
end

test_derivation_leftmost_list() = begin
    G = Grammar("""
        S -> A B
        A -> a
        B -> b
    """)
    d = leftmost(Derivation(G),[1, 2, 3])
    s = [(1, 1), (2, 1), (3, 2)]
    s == steps(d)
end

test_derivation_leftmost_ncf() = Grammar("S -> s\nT U ->s", false) |> Derivation |> (x -> leftmost(x, 1))

test_derivation_leftmost_allterminals() = begin
    G = Grammar("""
    E -> M | A | n
    M -> E * E
    A -> E + E
    """)
    d = Derivation(G) |>
        (x -> leftmost(x, 2)) |>
        (x -> leftmost(x, 5)) |>
        (x -> leftmost(x, 1)) |>
        (x -> leftmost(x, 4)) |>
        (x -> leftmost(x, 3)) |>
        (x -> leftmost(x, 3)) |>
        (x -> leftmost(x, 3))
    leftmost(d, 3)
end

test_derivation_leftmost_wrongsymbol() = begin
    G = Grammar("""
    E -> M | A | n
    M -> E * E
    A -> E + E
    """)
    d = Derivation(G) |>
        (x -> leftmost(x, 2)) |>
        (x -> leftmost(x, 5)) |>
        (x -> leftmost(x, 1)) |>
        (x -> leftmost(x, 4)) |>
        (x -> leftmost(x, 3))
    leftmost(d, 4)
end

test_derivation_rightmost() = begin
    G = Grammar("""
        S -> A B
        A -> a
        B -> b
    """)
    d = Derivation(G) |>
        (x -> rightmost(x, 1)) |>
        (x -> rightmost(x, 3)) |>
        (x -> rightmost(x, 2))
    s = [(1, 1), (3, 2), (2, 1)]
    s == steps(d)
end

test_derivation_rightmost_list() = begin
    G = Grammar("""
        S -> A B
        A -> a
        B -> b
    """)
    d = Derivation(G) |> x -> rightmost(x,[1, 3, 2])
    s = [(1, 1), (3, 2), (2, 1)]
    s == steps(d)
end

test_derivation_rightmost_ncf() = Grammar("S -> s\nT U ->s", false) |> Derivation |> x -> rightmost(x, 1)

test_derivation_rightmost_allterminals() = begin
    G = Grammar("""
    E -> M | A | n
    M -> E * E
    A -> E + E
    """)
    d = Derivation(G) |>
        (x -> rightmost(x, 1)) |>
        (x -> rightmost(x, 4)) |>
        (x -> rightmost(x, 3)) |>
        (x -> rightmost(x, 3))
        rightmost(d, 3)
end

test_derivation_rightmost_wrongsymbol() = begin
    G = Grammar("""
    E -> M | A | n
    M -> E * E
    A -> E + E
    """)
    d = Derivation(G) |>
        (x -> rightmost(x, 1)) |>
        (x -> rightmost(x, 4))
    rightmost(d, 4)
end

test_derivation_possible_steps() = begin
    G = Grammar("""
        S -> A B
        A -> a
        B -> b
    """, false)
    expected = [(2, 1), (3, 2)]
    actual = Derivation(G) |> (x -> next(x, 1, 1)) |> possiblesteps |> collect
    actual == expected
end

test_derivation_possible_steps_prod() = begin
    G = Grammar("""
        S -> A B
        A -> a
        B -> b
    """, false)
    expected = [(2, 1)]
    actual = Derivation(G) |> (x -> next(x, 1, 1)) |> (x-> possiblesteps(x, prod = 2)) |> collect
    expected == actual
end

test_derivation_possible_steps_pos() = begin
    G = Grammar("""
        S -> A B
        A -> a
        B -> b
    """, false)
    expected = [(3, 2)]
    actual = Derivation(G) |> (x -> next(x, 1, 1)) |> (x-> possiblesteps(x, pos = 2)) |> collect
    expected == actual
end

# Tests

function runderivationtests()
    @testset "derivation_tests" begin
        @test test_derivation_repr()
        @test test_derivation_eq()
        @test test_derivation_hash()
        @test test_derivation_steps()
        @test test_derivation_byprods()
        @test test_derivation_byprod()
        @test test_derivation_steps_list()
        @test_throws ArgumentError test_derivation_wrong_step()
        @test test_derivation_leftmost()
        @test test_derivation_leftmost_list()
        @test_throws ArgumentError test_derivation_leftmost_ncf()
        @test_throws ArgumentError test_derivation_leftmost_allterminals()
        @test_throws ArgumentError test_derivation_leftmost_wrongsymbol()
        @test test_derivation_rightmost()
        @test test_derivation_rightmost_list()
        @test_throws ArgumentError test_derivation_rightmost_ncf()
        @test_throws ArgumentError test_derivation_rightmost_allterminals()
        @test_throws ArgumentError test_derivation_rightmost_wrongsymbol()
        @test test_derivation_possible_steps()
        @test test_derivation_possible_steps_prod()
        @test test_derivation_possible_steps_pos()
    end
end