text:
"""
    TRMMDaily{ST<:AbstractString, DT<:TimeType} <: TRMMDataset

Object containing information on Daily TRMM datasets to be downloaded
"""
struct TRMMDaily{ST<:AbstractString, DT<:TimeType} <: TRMMDataset
	npdID :: ST
	lname :: ST
	doi   :: ST
    dtbeg :: DT
    dtend :: DT
    sroot :: ST
    hroot :: ST
    fpref :: ST
    fsuff :: ST
end

"""
    TRMMDaily(
        ST = String,
        DT = Date;
        dtbeg :: TimeType,
        dtend :: TimeType,
        sroot :: AbstractString = homedir(),
    ) -> npd :: TRMMDaily{ST,DT}

Creates a `TRMMDaily` dataset `npd` to retrieve datasets from the Final post-processing runs for Daily output

Keyword Arguments
=================
- `dtbeg` : Date at which download / analysis of the dataset begins
- `dtend` : Date at which download / analysis of the dataset ends
- `sroot` : The directory in which the folder `trmmdaily` will be created for data downloads, storage and analysis, default is the home directoy called by `homedir()`

The following fields in `npd` will be fixed as below:
- `npdID` : trmmdaily
- `lname` : Final TRMM Daily
- `doi`   : 10.5067/TRMM/TMPA/DAY/7
- `hroot` : https://disc2.gesdisc.eosdis.nasa.gov/opendap/TRMM_L3/TRMM_3B42_Daily.7
- `fpref` : 3B42_Daily
- `fsuff` : 7.nc4
"""
function TRMMDaily(
    ST = String,
    DT = Date;
    dtbeg :: TimeType,
    dtend :: TimeType,
    sroot :: AbstractString = homedir(),
)

	@info "$(modulelog()) - Setting up data structure containing information on Final TRMM Daily data to be downloaded"

    fol = joinpath(sroot,"trmmdaily"); if !isdir(fol); mkpath(fol) end

	dtbeg = Date(year(dtbeg),month(dtbeg),1)
	dtend = Date(year(dtend),month(dtend),daysinmonth(dtend))

    return TRMMDaily{ST,DT}(
		"trmmdaily", "Final TRMM Daily", "10.5067/TRMM/TMPA/DAY/7",
        dtbeg, dtend,
		joinpath(sroot,"trmmdaily"),
        "https://disc2.gesdisc.eosdis.nasa.gov/opendap/TRMM_L3/TRMM_3B42_Daily.7",
        "3B42_Daily", "7.nc4",
    )

end

"""
    TRMMDailyNRT(
        ST = String,
        DT = Date;
        dtbeg :: TimeType,
        dtend :: TimeType,
        sroot :: AbstractString = homedir(),
    ) -> npd :: TRMMDaily{ST,DT}

Creates a `TRMMDaily` dataset `npd` to retrieve datasets from the Near Real-Time processing runs for Daily output

Keyword Arguments
=================
- `dtbeg` : Date at which download / analysis of the dataset begins
- `dtend` : Date at which download / analysis of the dataset ends
- `sroot` : The directory in which the folder `trmmdailynrt` will be created for data downloads, storage and analysis, default is the home directoy called by `homedir()`

The following fields in `npd` will be fixed as below:
- `npdID` : trmmdailynrt
- `lname` : Near Real-Time TRMM Daily
- `doi`   : 10.5067/TRMM/TMPA/DAY-E/7
- `hroot` : https://disc2.gesdisc.eosdis.nasa.gov/opendap/TRMM_RT/TRMM_3B42RT_Daily.7
- `fpref` : 3B42RT_Daily
- `fsuff` : 7.nc4
"""
function TRMMDailyNRT(
    ST = String,
    DT = Date;
    dtbeg :: TimeType,
    dtend :: TimeType,
    sroot :: AbstractString = homedir(),
)

	@info "$(modulelog()) - Setting up data structure containing information on Near Real-Time TRMM Daily data to be downloaded"

    fol = joinpath(sroot,"trmmdailynrt"); if !isdir(fol); mkpath(fol) end

	dtbeg = Date(year(dtbeg),month(dtbeg),1)
	dtend = Date(year(dtend),month(dtend),daysinmonth(dtend))

    return TRMMDaily{ST,DT}(
		"trmmdailynrt", "Near Real-Time TRMM Daily", "10.5067/TRMM/TMPA/DAY-E/7",
        dtbeg, dtend,
		joinpath(sroot,"trmmdailynrt"),
        "https://disc2.gesdisc.eosdis.nasa.gov/opendap/TRMM_RT/TRMM_3B42RT_Daily.7",
        "3B42RT_Daily", "7.nc4",
    )

end

function show(io::IO, npd::TRMMDaily{ST,DT}) where {ST<:AbstractString, DT<:TimeType}
    print(
		io,
		"The NASA Precipitation Dataset {$ST,$DT} has the following properties:\n",
		"    Dataset ID      (npdID) : ", npd.npdID, '\n',
		"    Logging Name    (lname) : ", npd.lname, '\n',
		"    DOI URL          (doi)  : ", npd.doi,   '\n',
		"    Data Directory  (sroot) : ", npd.sroot, '\n',
		"    Date Begin      (dtbeg) : ", npd.dtbeg, '\n',
		"    Date End        (dtend) : ", npd.dtend, '\n',
		"    Timestep                : 1 Day\n",
        "    Data Resolution         : 0.25º\n",
        "    Data Server     (hroot) : ", npd.hroot, '\n',
	)
end

max_stars_repo_name: natgeo-wong/NASAPrecipitation

max_issues_repo_name: natgeo-wong/NASAPrecipitation

max_forks_repo_name: natgeo-wong/NASAPrecipitation


text:
using Pkg
Pkg.add("UrlDownload")
Pkg.add("DataFrames")
Pkg.add("CSV")

using UrlDownload
using DataFrames
using CSV

url = "https://covidtracking.com/api/v1/states/current.csv"
df = urldownload(url) |> DataFrame

ext: jl


text:
# grid
lat=30
lon=0

# atmospheric grid
ΔZA=50; # atmospheric layer thickness
vleva=50; # number of atmospheric vertical levels
ZAF=ΔZA*[0:(vleva);];
ZAC=ZAF[1:end-1].+ΔZA/2

# ocean grid
ΔZO=1; # oceanic layer thickness
vlevo=100; # number of oceanic vertical levels
# edges
ZOF=zeros(vlevo+1);
ZOF[2:vlevo+1]=-[1:vlevo;]*ΔZO;
ZOC=ZOF[1:end-1].+ΔZO/2


text:
# Some **really cool** maths:
#
# $$ \exp(i\pi) + 1 \quad = \quad 0 $$
#
# We can show this with some code:

x = exp(im*π) + 1

# that looks close to zero but

x ≈ 0

# however

abs(x) < eps()

# #### Conclusion
#
# The equation is proven thanks to our very rigorous proof.

max_stars_repo_name: mcognetta/franklindocs

max_issues_repo_name: mcognetta/franklindocs

max_forks_repo_name: mcognetta/franklindocs


text: # Julia Pkg dummy file, created to make the package instantiation process happy

max_forks_repo_name: MillironX/wphl-bioinformatics


text:
#######################################
#
# BSA file
#
#######################################

"""
Read Besa's BSA (.bsa) file

#### Input
* `fname`: Name or path for the BSA file

#### Output
* `bsa`: Dipole object
"""
function read_bsa(fname::AbstractString)
    @info("Reading BSA file = $fname")

    # Open file
    file = open(fname, "r")

    # Read version
    first_line        = readline(file)
    separator         = search(first_line, '|')
    version           = first_line[1:separator-1]
    coordinate_system = first_line[separator+1:end-1]

    # Read title line
    regexp = r"(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)"
    m = match(regexp, readline(file))

    # Read useless line
    readline(file)

    # Read dipoles
    dips = EEG.Dipole[]
    while !eof(file)
        dm = match(regexp, readline(file))

        dip = Dipole(coordinate_system, 1u"m" * float(dm.captures[2])/1000, 1u"m" * float(dm.captures[3])/1000, 1u"m" * float(dm.captures[4])/1000,
                                        float(dm.captures[5]), float(dm.captures[6]), float(dm.captures[7]),
                                        float(dm.captures[8]), float(dm.captures[9]), float(dm.captures[10]))

        push!(dips, dip)
    end

    # Close file
    close(file)

    @debug("Version = $version")
    @debug("Coordinate System  = $coordinate_system")
    @debug("Dipoles = $(length(dips))")

    return dips
end


text:
#=
Low level utilityies for dealing with bid sequences.
=#

import Base: isequal

const bitlen = 8

bid2index(n, r) = UInt8(20n + 2r - 20)
index2bid(i) = i % 2 == 0 ? divrem((i + 20) ÷ 2, 10) : (0,0)
parsebid(b) = divrem(b, 10)

struct UInt256 <: Integer
    h::UInt128
    l::UInt128
    UInt256(x) = new(0, x)
    UInt256(x, y) = new(UInt128(x), UInt128(y))
end

Base.hash(x::UInt256, h::UInt) = hash(x.h, hash(x.l, h))
Base.:(==)(x::UInt256, y::UInt256) = (x.h == y.h) & (x.l == y.l)
isequal(x::UInt256, y::UInt256) = isequal(x.h, y.h) & isequal(x.l, y.l)
Base.isless(x::UInt256, y::UInt256) = isequal(x.h, y.h) ? isless(x.l, y.l) : isless(x.h, y.h)
Base.:<(x::UInt256, y::UInt256) = isless(x, y)
Base.:>(x::UInt256, y::UInt256) = y < x
Base.:<<(x::UInt256, y::Integer) = UInt256(x.h << y | (x.l >> (128 - y)), x.l << y)
Base.:<<(x::UInt256, y::Int) = UInt256(x.h << y | (x.l >> (128 - y)), x.l << y)
Base.:<<(x::UInt256, y::Unsigned) = UInt256(x.h << y | (x.l >> (128 - y)), x.l << y)
Base.:>>(x::UInt256, y::Integer) = UInt256(x.h >> y, x.l >> y | (x.h << (128 - y)))
Base.:>>(x::UInt256, y::Int) = UInt256(x.h >> y, x.l >> y | (x.h << (128 - y)))
Base.:>>(x::UInt256, y::Unsigned) = UInt256(x.h >> y, x.l >> y | (x.h << (128 - y)))
Base.:>>>(x::UInt256, y::Integer) = UInt256(x.h >>> y, x.l >>> y | (x.h << (128 - y)))
Base.:>>>(x::UInt256, y::Int) = UInt256(x.h >>> y, x.l >>> y | (x.h << (128 - y)))
Base.:>>>(x::UInt256, y::Unsigned) = UInt256(x.h >>> y, x.l >>> y | (x.h << (128 - y)))
Base.:|(x::UInt256, y::UInt256) = UInt256(x.h | y.h, x.l | y.l)
Base.:|(x::UInt256, y::Unsigned) = UInt256(x.h, x.l | y)
Base.:⊻(x::UInt256, y::UInt256) = UInt256(x.h ⊻ y.h, x.l ⊻ y.l)
Base.:⊻(x::UInt256, y::Unsigned) = UInt256(x.h, x.l ⊻ y.l)
Base.:&(x::UInt256, y::UInt256) = UInt256(x.h & y.h, x.l & y.l)
Base.bitstring(x::UInt256) = bitstring(x.h) * bitstring(x.l)
Base.zero(x::UInt256) = UInt256(0)
Base.zero(UInt256) = UInt256(0)
Base.one(x::UInt256) = UInt256(1)
Base.one(UInt256) = UInt256(1)

peek(x, b) = isless(zero(x), one(x) << UInt8(b) & x)

togglebit(x::UInt256, b::UInt8) = x ⊻ (UInt256(0,1) << b)

function lastbid(hist::UInt256)
    if hist == zero(UInt256)
        return (0, 0)
    end
    if hist.h > zero(UInt128)
        i = UInt8(255)
        x = hist.h
    else
        i = UInt8(127)
        x = hist.l
    end
    while !peek(x, 127)
        x  = x << 0x1
        i -= 0x1
    end
    index2bid(i)
end

function addbid(hist::UInt256, n, r)
    (q, p) = lastbid(hist)
    if n == 0 && q != 0
        return togglebit(hist, bid2index(q, p) + 0x1)
    elseif n != 0
        return togglebit(hist, bid2index(n, r))
    else
        error("Cannot challenge a challenge")
    end
end

addbid(hist::UInt256, b::Integer) = addbid(hist, parsebid(b)...)

function bidindices(hist::UInt256)
    indices = []
    for idx = 0:255
        if isless(UInt256(0), UInt256(1) & hist)
            push!(indices, idx)
        end
        hist = hist >> 0x1
    end
    indices
end

bidhistory(hist::UInt256) = map(index2bid, bidindices(hist))

max_stars_repo_name: jeffreyrosenbluth/LP

max_issues_repo_name: jeffreyrosenbluth/LP

max_forks_repo_name: jeffreyrosenbluth/LP


text:
__precompile__()
"""
This is the main module and provides interfaces for specifying and solving reachability problems.
"""
module Reachability

using Reexport, RecipesBase, Memento, MathematicalSystems, HybridSystems,
      Compat, Suppressor
@reexport using LazySets
using LazySets.Approximations: overapproximate

import LazySets.use_precise_ρ
import LazySets.LinearMap

include("logging.jl")
include("Options/dictionary.jl")
include("Options/validation.jl")
include("Options/default_options.jl")
include("Utils/Utils.jl")
include("Properties/Properties.jl")
include("ReachSets/ReachSets.jl")

@reexport using Reachability.Utils,
                Reachability.ReachSets,
                Reachability.Properties

export project

include("solve.jl")
include("plot_recipes.jl")

# specify behavior of line search algorithm by dispatching on post operator
global discrete_post_operator = ConcreteDiscretePost()

@suppress_err begin
    function use_precise_ρ(cap::Intersection{N})::Bool where N<:Real
        return use_precise_ρ(discrete_post_operator, cap)
    end
end

end # module


text:
module Integrators

export rk4

function rk4(f::Function, y0::Array{Float64, 1}, t0::Float64,
             t1::Float64, h::Float64; inplace=true)
    y = y0
    n = round(Int, (t1 - t0)/h)
    t = t0
    if ~inplace
        hist = zeros(n, length(y0))
    end
    for i in 1:n
        k1 = h * f(t, y)
        k2 = h * f(t + 0.5*h, y + 0.5*k1)
        k3 = h * f(t + 0.5*h, y + 0.5*k2)
        k4 = h * f(t + h, y + k3)
        y = y + (k1 + 2*k2 + 2*k3 + k4)/6
        if ~inplace
            hist[i, :] = y
        end
        t = t0 + i*h
    end
    if ~inplace
        return hist
    else
        return y
    end
end

end


text:
using FastaIO
using Lumberjack
function ens2ucsc(seq_id::String)
    if seq_id == "MT"
        return "chrM"
    else
        return "chr"*seq_id
    end
end

function ucsc2ens(seq_id::String)
   if seq_id == "chrM"
        return "MT"
    else
        return "chr"*seq_id
    end
end

# yes this could be done via sed and cat
# but here it is.
function make_diploid(fastafile1, id1, fastafile2, id2, output_fasta)
    #run through each assign id to description
    fw = FastaWriter(output_fasta)
    fr1 = FastaReader{Vector{Char}}(fastafile1)
    Lumberjack.info("Start reading fasta file $fastafile1")
    for (desc, seq) in fr1
        desc = desc * "_$id1"
        writeentry(fw, desc, seq)
    end
    close(fr1)

    fr2 = FastaReader{Vector{Char}}(fastafile2)
    Lumberjack.info("Start reading fasta file $fastafile2")
    for (desc, seq) in fr2
        desc = desc * "_$id2"
        writeentry(fw, desc, seq)
    end
    close(fr2)
    close(fw)
    Lumberjack.info("Done writing fasta file $output_fasta")
end


text:
using Mustache, HttpCommon, ThreadTools
using DataStructures: OrderedDict
using StatsBase: RealVector, fweights
using Printf
import Statistics: quantile, mean

# used tableau 10 colormap
HTML_COLORS = ["#4E79A7","#F28E2C","#E15759","#76B7B2","#59A14F", "#EDC949","#AF7AA1","#FF9DA7","#9C755F","#BAB0AB"]

quantile(v::Dict{Int, Int}, p::Number) = quantile(v |> keys |> collect, v |> values |> collect |> fweights, [p])[1]
quantile(v::Dict{Int, Int}, p::RealVector) = quantile(v |> keys |> collect, v |> values |> collect |> fweights, p)
mean(v::Dict{Int, Int}) = mean(v |> keys |> collect, v |> values |> collect |> fweights)
stringify(arg::Pair; max_len=1_000) = stringify(arg.first, arg.second, max_len=max_len)
stringify(arg1::String, arg2; max_len=1_000) = isnothing(max_len) ? "$(escapeHTML(arg1)): $(arg2)" : "$(escapeHTML(first(arg1, max_len))): $(arg2)"
stringify(arg1, arg2; max_len=1_000) = "$(arg1): $(arg2)"
calc_filled_percent(e, parent_updated) = isnothing(parent_updated) ? "" : @sprintf ", filled=%.2f%%" (100 * e.updated / parent_updated)
pad_color(pad) = HTML_COLORS[((length(pad)÷2)%length(HTML_COLORS))+1]

function pair2html(i, key, val, max_len, pad)
	pair_repr = stringify(key, val, max_len=max_len)
	pad*" "^2 * "
$pair_repr
\n"
end

function schema2html(e::Entry; pad = "", max_vals=100, max_len=1_000, parent_updated=nothing, parent_key="")
	c = pad_color(pad)
	filled_percent = calc_filled_percent(e, parent_updated)
	sorted_counts = sort(collect(e.counts), by=x->x[2], rev=true)
	min_repr = stringify(sorted_counts[end], max_len=max_len)
	max_repr = stringify(sorted_counts[1], max_len=max_len)
	ret_str = """
$pad
$pad[Scalar - $(join(types(e)))], $(length(keys(e.counts))) unique values,
(updated=$(e.updated)$filled_percent, min=$min_repr, max=$max_repr)
"""

	counts2process = if isnothing(max_vals)
		enumerate(sorted_counts)
	else
		Iterators.take(enumerate(sorted_counts), max_vals)
	end
	li_part = join(tmap(i->pair2html(i[1], i[2].first, i[2].second, max_len, pad), Threads.nthreads(), counts2process))
	ret_str *= li_part
	if !isnothing(max_vals) && length(sorted_counts) > max_vals
		ret_str *= pad*" "^2 * "
and other $(length(e.counts) - max_vals) values
\n"
	end
	ret_str * "$pad
\n"
end

# because sometimes there is empty list in all jsons, this helps to determine the pruning of such element
function schema2html(e::Nothing; pad = "", max_vals=100, max_len=1_000, parent_updated=nothing, parent_key="")
	c = HTML_COLORS[((length(pad)÷2)%length(HTML_COLORS))+1]
	ret_str = """
$pad[Empty list element], this list is empty in all JSONs, can not infer schema, suggesting to delete key $parent_key
"""
	ret_str
end

function schema2html(e::ArrayEntry; pad = "", max_vals=100, max_len=1_000, parent_updated=nothing, parent_key="")
 	c = pad_color(pad)
	# todo: fix it all so it is different method for array of entries and the rest so only nested things are truly nested
	filled_percent = calc_filled_percent(e, parent_updated)
	quantiles = quantile(e.l, [0.1, 0.5, 0.9])
	min_val = minimum(keys(e.l))
	max_val = maximum(keys(e.l))
	mean_val = round(mean(e.l), digits=2)
	q10_val = round(quantiles[1], digits=2)
	median_val = round(quantiles[2], digits=2)
	q90_val = round(quantiles[3], digits=2)
	ret_str = """
$pad
[List] (updated=$(e.updated)$filled_percent, mean=$mean_val,
min=$min_val, max=$max_val, 10th percentile=$q10_val, median=$median_val, 90th percentile=$q90_val)
$pad
with following frequencies
$pad

"""
 	i = 0
	sorted_el = sort(collect(e.l))
 	for (key, val) in sorted_el
		ret_str *= pad*" "^2 * "
$key: $val
\n"
		i += 1
		if !isnothing(max_vals) && i == max_vals && length(sorted_el) > max_vals
			ret_str *= pad*" "^2 * "
and other $(length(e.l) - i) values
\n"
			break
		end
 	end
	ret_str *= """
$pad

$pad

$pad
and data
"""
	ret_str *= schema2html(e.items, pad=pad*" "^2, max_vals=max_vals, max_len=max_len, parent_key="$parent_key.items")
	ret_str * """
$pad

$pad

"""
end

function schema2html(e::DictEntry; pad = "", max_vals=100, max_len=1_000, parent_updated=nothing, parent_key="")
	c = pad_color(pad)
    if isempty(e.childs)
    	return pad * """
Empty Dict
\n"""
    end
	class = isempty(pad) ? "top_dict" : "nested"
	filled_percent = calc_filled_percent(e, parent_updated)
	ret_str = pad * """
[Dict] (updated=$(e.updated)$filled_percent)\n"""
	i = 0
	ordered_childs = OrderedDict(e.childs)
    for (key, val) in sort!(ordered_childs)
		child_key = """$parent_key[$(repr(Symbol(key)))]"""
		ret_str *= pad*" "^2 * """
$key - $child_key<input type="checkbox" name="$(escapeHTML(child_key))" value="$(escapeHTML(child_key))">\n"""
		ret_str *= schema2html(val, pad=pad*" "^4, max_vals=max_vals, max_len=max_len, parent_updated=e.updated, parent_key=child_key)
		ret_str *= pad*" "^2 * "
\n"
		i += 1
		if !isnothing(max_vals) && i == max_vals && length(ordered_childs) > max_vals
			ret_str *= pad*" "^2 * "
and other $(length(e.childs) - i) values
\n"
			break
		end
    end
	ret_str * pad * "
\n"
end

function schema2html(e::MultiEntry; pad = "", max_vals=100, max_len=1_000, parent_updated=nothing, parent_key="")
	c = pad_color(pad)
	filled_percent = calc_filled_percent(e, parent_updated)
	ret_str = pad * """
[MultiEntry] (updated=$(e.updated)$filled_percent)\n"""
	i = 0
    for (key, val) in enumerate(e.childs)
		child_key = """$parent_key["$key"]"""
		ret_str *= pad*" "^2 * """
$key - $child_key<input type="checkbox" name="$(escapeHTML(child_key))" value="$(escapeHTML(child_key))">\n"""
		ret_str *= schema2html(val, pad=pad*" "^4, max_vals=max_vals, max_len=max_len, parent_updated=e.updated, parent_key=child_key)
		ret_str *= pad*" "^2 * "
\n"
		i += 1
		if !isnothing(max_vals) && i == max_vals && length(e.childs) > max_vals
			ret_str *= pad*" "^2 * "
and other $(length(e.childs) - i) values
\n"
			break
		end
    end
	ret_str * pad * "
\n"
end
# queryselectorall v js místo getelementsbyclassname
# a zkusit minifikovat to html-vyházet odsazení, a kouknout na rozdíl velikostí

"""
	generate_html(sch::DictEntry; max_vals=100, max_len=1_000)
	generate_html(file_name, sch::DictEntry; max_vals=100, max_len=1_000)

exports schema to HTML including CSS style and JS allowing to expand / hide
sub-parts of schema, countmaps, and lengthmaps.

# Arguments
- `max_vals` controls maximum number of exported values in countmap
- `max_len` controls maximum number of exported lengts of arrays
- `file_name` a name of file to save HTML with schema

# Return
If provided filename, it does not return anything.
If not, it returns the generated HTML+CSS+JS as a String.

# Example
You can either open the html file in any browser, or open it directly using ElectronDisplay

```julia
using ElectronDisplay
using ElectronDisplay: newdisplay
generated_html = generate_html(sch, max_vals = 100)
display(newdisplay(), MIME{Symbol("text/html")}(), generated_html)
```
"""
function generate_html(sch::DictEntry; max_vals=100, max_len=1_000)
	tpl = mt"""
	<!DOCTYPE html>

    	<meta charset="UTF-8">
    	<title>Json schema dump</title>

		    <button type="button" id="copy_clipboard">Copy selectors to clipboard</button>

		{{{list_dump}}}

"""

	d = Dict(
		"list_dump" => schema2html(sch, max_vals=max_vals, max_len=max_len),
	)
	return render(tpl, d)
end

function generate_html(file_name, sch::DictEntry; max_vals=100, max_len=1_000)
	s = generate_html(sch; max_vals=max_vals, max_len=max_len)
	open(file_name, "w") do f
 		write(f, s)
	end
end

# Base.show(io, m::MIME{Symbol("text/html")}, sch::JsonGrinder.DictEntry) = print(io, repr(m,sch))

Base.repr(::MIME"text/html", sch::DictEntry; max_vals = 100, max_len=1_000, context = nothing) = generate_html(sch; max_vals=max_vals, max_len=max_len)
Base.showable(::MIME"text/html", ::DictEntry) = true


text:
# Autogenerated wrapper script for Libgcrypt_jll for x86_64-unknown-freebsd
export libgcrypt

using Libgpg_error_jll
JLLWrappers.@generate_wrapper_header("Libgcrypt")
JLLWrappers.@declare_library_product(libgcrypt, "libgcrypt.so.22")
function __init__()
    JLLWrappers.@generate_init_header(Libgpg_error_jll)
    JLLWrappers.@init_library_product(
        libgcrypt,
        "lib/libgcrypt.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using TextParse

import TextParse: tryparsenext, unwrap, failedat, AbstractToken, LocalOpts
using Base.Test

# dumb way to compare two AbstractTokens
Base.:(==){T<:AbstractToken}(a::T, b::T) = string(a) == string(b)

import TextParse: eatnewlines
@testset "eatnewlines" begin
    @test eatnewlines("\n\r\nx") == (4, 2)
    @test eatnewlines("x\n\r\nx") == (1, 0)
end

import TextParse: getlineend
@testset "getlineend" begin
    @test getlineend("\nx") == 0
    @test getlineend("x\nx") == 1
    @test getlineend("x\ny", 2) == 1
    @test getlineend("x\nyz", 3) == 4
end

import TextParse: fromtype, Percentage
@testset "Float parsing" begin

    @test tryparsenext(fromtype(Float64), "1", 1, 1) |> unwrap == (1.0, 2)
    @test tryparsenext(fromtype(Float64), "12", 1, 2) |> unwrap == (12.0, 3)
    @test tryparsenext(fromtype(Float64), ".1", 1, 2) |> unwrap == (0.1, 3)
    @test tryparsenext(fromtype(Float64), "1.1", 1, 3) |> unwrap == (1.1, 4)
    @test tryparsenext(fromtype(Float32), "1.", 1, 2) |> unwrap == (1f0,3)
    @test tryparsenext(fromtype(Float64), "-1.1", 1, 4) |> unwrap == (-1.1,5)
    @test tryparsenext(fromtype(Float64), "-1.0e-12", 1, 8) |> unwrap == (-1.0e-12,9)
    @test tryparsenext(fromtype(Float64), "-1e-12") |> unwrap == (-1.0e-12,7)
    @test tryparsenext(fromtype(Float64), "-1.0E-12", 1, 8) |> unwrap == (-1.0e-12,9)
    @test tryparsenext(Percentage(), "33%") |> unwrap == (.33,4)
    @test tryparsenext(Percentage(), "3.3%") |> unwrap == (.033,5)
end

import TextParse: StringToken
@testset "String parsing" begin

    # default options
    @test tryparsenext(StringToken(String), "") |> unwrap == ("", 1)
    x = "x"
    @test tryparsenext(StringToken(String), "x") |> unwrap == ("x", 2)
    @test tryparsenext(StringToken(String), "x ") |> unwrap == ("x ", 3)
    @test tryparsenext(StringToken(String), " x") |> unwrap == (" x", 3)
    @test tryparsenext(StringToken(String), "x\ny") |> unwrap == ("x", 2)
    @test tryparsenext(StringToken(String), "x,y") |> unwrap == ("x", 2) # test escape

    opts = LocalOpts(',', false, '"', '"', true, true)
    @test tryparsenext(StringToken(String), "", opts) |> unwrap == ("", 1)
    @test tryparsenext(StringToken(String), "\"\"", opts) |> unwrap == ("\"\"", 3)
    @test tryparsenext(StringToken(String), "x", opts) |> unwrap == ("x", 2)
    # test including new lines
    @test tryparsenext(StringToken(String), "x\ny", opts) |> unwrap == ("x\ny", 4)
    @test tryparsenext(StringToken(String), "\"x\ny\"", opts) |> unwrap == ("\"x\ny\"", 6)

    opts = LocalOpts(',', false, '"', '"', false, true)
    # test that includequotes option doesn't affect string
    @test tryparsenext(StringToken(String), "\"\"", opts) |> unwrap == ("\"\"", 3)

    opts = LocalOpts(',', false, '"', '\\', false, false)
    str =  "Owner 2 ”Vicepresident\"\""
    @test tryparsenext(Quoted(String), str) |> unwrap == (str, endof(str)+1)
    str1 =  "\"Owner 2 ”Vicepresident\"\"\""
    @test tryparsenext(Quoted(String,quotechar=Nullable('"'), escapechar=Nullable('"')), str1) |> unwrap == (str, endof(str1)+1)
    @test tryparsenext(Quoted(String), "\"\tx\"") |> unwrap == ("\tx", 5)
    opts = LocalOpts(',', true, '"', '\\', false, false)
    @test tryparsenext(StringToken(String), "x y",1,3, opts) |> unwrap == ("x", 2)
end

import TextParse: Quoted, NAToken, Unknown
@testset "Quoted string parsing" begin
    opts = LocalOpts(',', false, '"', '"', true, true)

    @test tryparsenext(Quoted(String), "\"\"") |> unwrap == ("", 3)
    @test tryparsenext(Quoted(String), "\"\" ", opts) |> unwrap == ("", 3)
    @test tryparsenext(Quoted(String), "\"x\"") |> unwrap == ("x", 4)
    @test tryparsenext(Quoted(String, includequotes=true), "\"x\"") |> unwrap == ("\"x\"", 4)
    str2 =  "\"\"\"\""
    @test tryparsenext(Quoted(String), str2, opts) |> unwrap == ("\"\"", endof(str2)+1)
    str1 =  "\"x”y\"\"\""
    @test tryparsenext(Quoted(StringToken(String), required=true), "x\"y\"") |> failedat == 1

    @test tryparsenext(Quoted(String,quotechar=Nullable('"'), escapechar=Nullable('"')), str1) |> unwrap == ("x”y\"\"", endof(str1)+1)
    @test tryparsenext(Quoted(StringToken(String)), "\"x\\\"yz\"") |> unwrap == ("x\\\"yz", 8)
    @test tryparsenext(Quoted(NAToken(fromtype(Int))), "1") |> unwrap |> unwrap == (1,2)
    @test tryparsenext(Quoted(NAToken(fromtype(Int))), "") |> unwrap |> failedat == 1
    @test tryparsenext(Quoted(NAToken(fromtype(Int))), "\"\"") |> unwrap |> failedat == 3
    @test tryparsenext(Quoted(NAToken(fromtype(Int))), "\"1\"") |> unwrap |> unwrap == (1, 4)

    @test tryparsenext(Quoted(StringToken(String)), "\"abc\"") |> unwrap == ("abc", 6)
    @test tryparsenext(Quoted(StringToken(String)), "x\"abc\"") |> unwrap == ("x\"abc\"", 7)
    @test tryparsenext(Quoted(StringToken(String)), "\"a\nbc\"") |> unwrap == ("a\nbc", 7)
    @test tryparsenext(Quoted(StringToken(String), required=true), "x\"abc\"") |> failedat == 1
    @test tryparsenext(Quoted(fromtype(Int)), "21") |> unwrap == (21,3)
    @test tryparsenext(Quoted(NAToken(fromtype(Int))), "21") |> unwrap |> unwrap == (21,3)
    @test tryparsenext(Quoted(NAToken(fromtype(Int))), "") |> unwrap |> failedat == 1
    @test tryparsenext(Quoted(NAToken(fromtype(Int))), "\"\"") |> unwrap |> failedat == 3
    @test tryparsenext(Quoted(NAToken(fromtype(Int))), "\"21\"") |> unwrap |> unwrap == (21, 5)
    @test isnull(tryparsenext(Quoted(NAToken(Unknown())), " ") |> unwrap |> first)
    opts = LocalOpts(',', false,'"', '"', false, false)
    @test tryparsenext(Quoted(StringToken(String)), "x,", opts) |> unwrap == ("x", 2)

    # stripspaces
    @test tryparsenext(Quoted(Percentage()), "\" 10%\",", opts) |> unwrap == (0.1, 7)
    @test tryparsenext(Quoted(String), "\" 10%\",", opts) |> unwrap == (" 10%", 7)
    opts = LocalOpts(',', true,'"', '"', false, false)
    @test tryparsenext(Quoted(StringToken(String)), "\"x y\" y", opts) |> unwrap == ("x y", 6)
    @test tryparsenext(Quoted(StringToken(String)), "x y", opts) |> unwrap == ("x", 2)
end

@testset "NA parsing" begin
    @test tryparsenext(NAToken(fromtype(Float64)), ",") |> unwrap |> failedat == 1 # is nullable
    @test tryparsenext(NAToken(fromtype(Float64)), "X,") |> failedat == 1
    @test tryparsenext(NAToken(fromtype(Float64)), "NA,") |> unwrap |> failedat == 3
    @test tryparsenext(NAToken(fromtype(Float64)), "1.212,") |> unwrap |> unwrap == (1.212, 6)
end

import TextParse: Field
@testset "Field parsing" begin
    f = fromtype(Int)
    @test tryparsenext(Field(f), "12,3") |> unwrap == (12, 4)
    @test tryparsenext(Field(f), "12 ,3") |> unwrap == (12, 5)
    @test tryparsenext(Field(f), " 12 ,3") |> unwrap == (12, 6)
    opts = LocalOpts('\t', false, 'x','x',true,false)
    @test tryparsenext(Field(f), "12\t3", 1, 4, opts) |> unwrap == (12, 4)
    @test tryparsenext(Field(f), "12 \t3", 1, 5, opts) |> unwrap == (12, 5)
    @test tryparsenext(Field(f), " 12 \t 3", 1, 6, opts) |> unwrap == (12, 6)
    opts = LocalOpts('\t', true, 'x','x',true,false)
    @test tryparsenext(Field(f), " 12 3", 1, 5, opts) |> unwrap == (12, 5)
    @test tryparsenext(Field(f, ignore_end_whitespace=false), " 12 \t 3", 1,6, opts) |> unwrap == (12, 5)
    opts = LocalOpts(' ', false, 'x','x',false, false)
    @test tryparsenext(Field(f,ignore_end_whitespace=false), "12 3", 1,4,opts) |> unwrap == (12, 4)
#    @test tryparsenext(Field(f,ignore_end_whitespace=false), "12 \t3", 1,5,opts) |> failedat == 3
    opts = LocalOpts('\t', false, 'x','x',false, false)
    @test tryparsenext(Field(f,ignore_end_whitespace=false), " 12\t 3", 1, 6, opts) |> unwrap == (12,5)
    @test tryparsenext(Field(f,eoldelim=true), " 12\n", 1, 4, opts) |> unwrap == (12,5)
    @test tryparsenext(Field(f,eoldelim=true), " 12\n\r\n", 1, 5, opts) |> unwrap == (12,6)
    @test tryparsenext(Field(f,eoldelim=true), " 12") |> unwrap == (12,4)
end

import TextParse: Record
@testset "Record parsing" begin
    r=Record((Field(fromtype(Int)), Field(fromtype(UInt)), Field(fromtype(Float64))))
    @test tryparsenext(r, "12,21,21,", 1, 9) |> unwrap == ((12, UInt(21), 21.0), 10)
    @test tryparsenext(r, "12,21.0,21,", 1, 9) |> failedat == 6
    s = "12   ,  21,  21.23,"
    @test tryparsenext(r, s, 1, length(s)) |> unwrap == ((12, UInt(21), 21.23), length(s)+1)
end

import TextParse: UseOne
@testset "UseOne" begin
    f = UseOne((Field(fromtype(Int)), Field(fromtype(Float64)), Field(fromtype(Int), eoldelim=true)), 3)
    @test tryparsenext(f, "1, 33.21, 45", 1, 12) |> unwrap == (45, 13)
end

import TextParse: Repeated
@testset "Repeated" begin
    f = Repeated(Field(fromtype(Int)), 3)
    @test tryparsenext(f, "1, 33, 45,", 1, 12) |> unwrap == ((1,33,45), 11)

    inp = join(map(string, [1:45;]), ", ") * ", "
    out = ntuple(identity, 45)
    f2 = Repeated(Field(fromtype(Int)), 45)
    @test tryparsenext(f2, inp, 1, length(inp)) |> unwrap == (out, length(inp))
    #@benchmark tryparsenext($f2, $inp, 1, length($inp))
end

import TextParse: quotedsplit
@testset "quotedsplit" begin
    opts = LocalOpts(',', false, '"', '\\', false, false)
    @test quotedsplit("x", opts, false, 1, 1) == ["x"]
    @test quotedsplit("x, y", opts, false, 1, 4) == ["x", "y"]
    @test quotedsplit("\"x\", \"y\"", opts,false, 1, 8) == ["x", "y"]
    @test quotedsplit("\"x\", \"y\"", opts,true, 1, 8) == ["\"x\"", "\"y\""]
    str = """x\nx,"s,", "\\",x" """
    @test quotedsplit(str, opts, false, 3, length(str)) == ["x", "s,", "\\\",x"]
    @test quotedsplit(",", opts, true, 1, 1) == ["", ""]
    @test quotedsplit(", ", opts, false, 1, 2) == ["", ""]
    str = "1, \"x \"\"y\"\" z\", 1"
    qopts = LocalOpts(',', false,'"', '"', false, false)
    @test quotedsplit(str, qopts,true, 1, endof(str)) == ["1", "\"x \"\"y\"\" z\"", "1"]
end

import TextParse: LocalOpts, readcolnames
@testset "CSV column names" begin
    str1 = """
     a, b,c d, e
    x,1,1,1
    ,1,1,1
    x,1,1.,1
    x y,1.0,1,
    x,1.0,,1
    """

    str2 = """
     a, " b", "c", "d\\" e "
    """
    opts = LocalOpts(',', false, '"', '\\', false, false)
    @test readcolnames(str1, opts, 1, String[]) == (["a", "b", "c d", "e"], 13)
    @test readcolnames("\n\r$str1", opts, 3, Dict(3=>"x")) == (["a", "b", "x", "e"], 15)
    #@test readcolnames("$str2", opts, 3, Dict(3=>"x")) == (["a", "b", "x", "d\" e"], 24)
end

import TextParse: guesstoken, Unknown, Numeric, DateTimeToken, StrRange
@testset "guesstoken" begin
    # Test null values
    @test guesstoken("", Unknown()) == NAToken(Unknown())
    @test guesstoken("null", Unknown()) == NAToken(Unknown())
    @test guesstoken("", NAToken(Unknown())) == NAToken(Unknown())
    @test guesstoken("null", NAToken(Unknown())) == NAToken(Unknown())

    # Test NA
    @test guesstoken("1", NAToken(Unknown())) == NAToken(Numeric(Int))
    @test guesstoken("1", NAToken(Numeric(Int))) == NAToken(Numeric(Int))
    @test guesstoken("", NAToken(Numeric(Int))) == NAToken(Numeric(Int))
    @test guesstoken("1%", NAToken(Unknown())) == NAToken(Percentage())

    # Test non-null numeric
    @test guesstoken("1", Unknown()) == Numeric(Int)
    @test guesstoken("1", Numeric(Int)) == Numeric(Int)
    @test guesstoken("", Numeric(Int)) == NAToken(Numeric(Int))
    @test guesstoken("1.0", Numeric(Int)) == Numeric(Float64)

    # Test strings
    @test guesstoken("x", Unknown()) == StringToken(StrRange)

    # Test nullable to string
    @test guesstoken("x", NAToken(Unknown())) == StringToken(StrRange)

    # Test string to non-null (short circuit)
    @test guesstoken("1", StringToken(StrRange)) == StringToken(StrRange)

    # Test quoting
    @test guesstoken("\"1\"", Unknown()) == Quoted(Numeric(Int))
    @test guesstoken("\"1\"", Quoted(Numeric(Int))) == Quoted(Numeric(Int))

    # Test quoting with Nullable tokens
    @test guesstoken("\"\"", Quoted(Unknown())) == Quoted(NAToken(Unknown()))
    @test guesstoken("\"\"", Quoted(NAToken(Unknown()))) == Quoted(NAToken(Unknown()))
    @test guesstoken("\"\"", Quoted(Numeric(Int))) == Quoted(NAToken(Numeric(Int)))
    @test guesstoken("\"\"", Unknown()) == Quoted(NAToken(Unknown()))
    @test guesstoken("\"\"", Numeric(Int)) == Quoted(NAToken(Numeric(Int)))
    @test guesstoken("", Quoted(Numeric(Int))) == Quoted(NAToken(Numeric(Int)))
    @test guesstoken("", Quoted(NAToken(Numeric(Int)))) == Quoted(NAToken(Numeric(Int)))
    @test guesstoken("1", Quoted(NAToken(Numeric(Int)))) == Quoted(NAToken(Numeric(Int)))
    @test guesstoken("\"1\"", Quoted(NAToken(Numeric(Int)))) == Quoted(NAToken(Numeric(Int)))

    # Test DateTime detection:
    tok = guesstoken("2016-01-01 10:10:10.10", Unknown())
    @test tok == DateTimeToken(DateTime, dateformat"yyyy-mm-dd HH:MM:SS.s")
    @test guesstoken("2016-01-01 10:10:10.10", tok) == tok
    @test guesstoken("2016-01-01 10:10:10.10", Quoted(NAToken(Unknown()))) == Quoted(NAToken(tok))
end

import TextParse: guesscolparsers
@testset "CSV type detect" begin
    str1 = """
     a, b,c d, e
    x,1,1,1
    x,1,1,1
    x,1,1.,1
    x y,1.0,1,
    ,1.0,,1
    """
    opts = LocalOpts(',', false, '"', '\\', false, false)
    _, pos = readcolnames(str1, opts, 1, String[])
    testtill(i, colparsers=[]) = guesscolparsers(str1, String[], opts, pos, i, colparsers)
    @test testtill(0) |> first == Any[]
    @test testtill(1) |> first == map(fromtype, [StrRange, Int, Int, Int])
    @test testtill(2) |> first == map(fromtype, [StrRange, Int, Int, Int])
    @test testtill(3) |> first == map(fromtype, [StrRange, Int, Float64, Int])
    @test testtill(4) |> first == vcat(map(fromtype, [StrRange, Float64, Float64]),
                                       NAToken(fromtype(Int)))
    @test testtill(5) |> first == vcat(map(fromtype, [StrRange, Float64]),
                                       NAToken(fromtype(Float64)),
                                       NAToken(fromtype(Int)))
end

import TextParse: getlineat
@testset "getlineat" begin
    str = "abc\ndefg"
    @test str[getlineat(str,1)] == "abc\n"
    @test str[getlineat(str,4)] == "abc\n"
    @test str[getlineat(str,5)] == "defg"
    @test str[getlineat(str,endof(str))] == "defg"
    @test getlineat("x", 5) == 1:1
end

import TextParse: guessdateformat
@testset "date detection" begin
    @test guessdateformat("2016") |> typeof == DateTimeToken(Date, dateformat"yyyy-mm-dd") |> typeof
    @test guessdateformat("09/09/2016") |> typeof == DateTimeToken(Date, dateformat"mm/dd/yyyy") |> typeof
    @test guessdateformat("24/09/2016") |> typeof == DateTimeToken(Date, dateformat"dd/mm/yyyy") |> typeof
end

@testset "date parsing" begin
    tok = DateTimeToken(DateTime, dateformat"yyyy-mm-dd HH:MM:SS")
    opts = LocalOpts('y', false, '"', '\\', false, false)
    str = "1970-02-02 02:20:20"
    @test tryparsenext(tok, str, 1, length(str), opts) |> unwrap == (DateTime("1970-02-02T02:20:20"), length(str)+1)
    @test tryparsenext(tok, str*"x", 1, length(str)+1, opts) |> unwrap == (DateTime("1970-02-02T02:20:20"), length(str)+1)
    @test tryparsenext(tok, str[1:end-3]*"x", 1, length(str)-2, opts) |> failedat == length(str)-2
    @test tryparsenext(tok, str[1:end-3]*"y", 1, length(str)-2, opts) |> unwrap == (DateTime("1970-02-02T02:20"), length(str)-2)
end

using DataValues
import TextParse: _csvread
@testset "csvread" begin

    str1 = """
     a, b,c d, e
    x,1,1,1
    ,1,,1
    x,1,1.,1
    x y,1.0,1,
    x,1.0,,1
    """
    data = ((["x", "","x","x y","x"],
              ones(5),
              DataValueArray(ones(5), Bool[0,1,0,0,1]),
              DataValueArray(ones(Int,5), Bool[0,0,0,1,0])),
              ["a", "b", "c d", "e"])
    @test isequal(_csvread(str1, ','), data)
    coltype_test1 = _csvread(str1,
                            colparsers=Dict("b"=>Nullable{Float64},
                                          "e"=>DataValue{Float64}))
    coltype_test2 = _csvread(str1,
                            colparsers=Dict(2=>Nullable{Float64},
                                          4=>DataValue{Float64}))

    str2 = """
    x,1,1,1
    ,1,,1
    x,1,1.,1
    x y,1.0,1,
    x,1.0,,1
    """
    coltype_test3 = _csvread(str2, header_exists=false,
                            colparsers=Dict(2=>Nullable{Float64},
                                          4=>Nullable{Float64}))
    @test eltype(coltype_test1[1][2]) == DataValue{Float64}
    @test eltype(coltype_test1[1][4]) == DataValue{Float64}
    @test eltype(coltype_test2[1][2]) == DataValue{Float64}
    @test eltype(coltype_test2[1][4]) == DataValue{Float64}
    @test eltype(coltype_test3[1][2]) == DataValue{Float64}
    @test eltype(coltype_test3[1][4]) == DataValue{Float64}

    @test isequal(data, _csvread(str1, type_detect_rows=1))
    @test isequal(data, _csvread(str1, type_detect_rows=2))
    @test isequal(data, _csvread(str1, type_detect_rows=3))
    @test isequal(data, _csvread(str1, type_detect_rows=4))
    # But we can't go from a non-string column to a string column
    str3 = """
     a, b,c d, e
    1,1,1,1
    x,1,1,1
    """
    @test_throws TextParse.CSVParseError _csvread(str3, type_detect_rows=1)

    # test growing of columns if prediction is too low
    @test _csvread("x,y\nabcd, defg\n,\n,\n", type_detect_rows=1) ==
        ((String["abcd", "", ""], String["defg", "", ""]), String["x", "y"])

    # #19
    s="""
    x,y,z
    1,1,x
    "2",2,x
    1,2,"x \"\"y\"\""
    """

    res = (([1, 2, 1], [1, 2, 2], String["x", "x", "x \"\"y\"\""]), String["x", "y", "z"])
    @test _csvread(s, type_detect_rows=1, escapechar='"') == res
    @test _csvread(s, type_detect_rows=2, escapechar='"') == res

    @test csvread(IOBuffer("x\n1")) == (([1],),["x"])

    @test _csvread("x\n1\n") == (([1],),["x"])

    # test detection of newlines in fields
    s = """x, y
        abc, def
        g
        hi,jkl
        mno,pqr
        """

    @test _csvread(s, type_detect_rows=1) == ((["abc", "g\nhi", "mno"], ["def", "jkl", "pqr"]), ["x", "y"])
    # test custom na strings
    s = """
    x,y
    1,2
    ?,3
    4,*
    """
    nullness = ([false, true, false], [false, false, true])
    @test map(x->x.isnull, first(_csvread(s, nastrings=["?","*"]))) == nullness
    @test map(x->x.isnull, first(_csvread(s, nastrings=["?","*"], type_detect_rows=1))) == nullness

    @test isequal(csvread(["data/a.csv", "data/b.csv"]),
                  (([1.0, 2.0, 1.0, 2.0, 3.0], DataValue{Int64}[2, 2, nothing, nothing, nothing],
                    DataValue{Int64}[nothing, nothing, nothing, 2, 1]), String["x", "y", "z"], [2, 3]))
    @test isequal(csvread(["data/a.csv", "data/b.csv"], samecols=[("y","z")]),
                  (([1.0, 2.0, 1.0, 2.0, 3.0], DataValue{Int64}[2, 2, nothing, 2, 1]), String["x", "y"], [2,3]))

    # shouldn't fail because y doesn't exist
    @test _csvread("x\n1", colparsers=Dict("y"=>String)) == (([1],), ["x"])

    # Don't try to guess type if it's provided by user. Issue JuliaDB.jl#109
    s="""
    time,value
    "2017-11-09T07:00:07.391101180",0
    """
    @test _csvread(s) == ((String["2017-11-09T07:00:07.391101180"], [0]), String["time", "value"])
    @test _csvread(s, colparsers=Dict(:time=>String)) == ((String["2017-11-09T07:00:07.391101180"], [0]), String["time", "value"])

    @test _csvread("") == ((), String[])
end

@testset "skiplines_begin" begin
    str1 = """
    hello

    world
    x,y,z
    1,1,1
    """
    @test _csvread(str1, skiplines_begin=3) == (([1], [1], [1]), String["x", "y","z"])

    s = """
    x,y z
    a,b 1
    e  	3
    """
    @test _csvread(s, spacedelim=true) == ((["a,b", "e"],[1,3]), ["x,y","z"])
end

using PooledArrays

@testset "pooled array promotion" begin
    # test default behavior
    xs = [randstring(10) for i=1:100]
    col = _csvread(join(xs, "\n"), header_exists=false)[1][1]
    @test isa(col, PooledArray)
    @test eltype(col.refs) == UInt8
    @test xs == col

    # test promotion to a widened type
    xs = [randstring(10) for i=1:300]
    col = _csvread(join(xs, "\n"), header_exists=false)[1][1]
    @test isa(col, PooledArray)
    @test eltype(col.refs) == UInt16
    @test xs == col

    # test promotion to a dense array
    xs = [randstring(10) for i=1:513]
    col = _csvread(join(xs, "\n"), header_exists=false)[1][1]
    @test !isa(col, PooledArray)
    @test isa(col, Array)
    @test xs == col

    # test non-promotion
    xs = [rand(["X", "Y"]) for i=1:500]
    col = _csvread(join(xs, "\n"), header_exists=false)[1][1]
    @test isa(col, PooledArray)
    @test eltype(col.refs) == UInt8
    @test xs == col
end

import TextParse: eatwhitespaces
@testset "custom parser" begin
    const floatparser = Numeric(Float64)
    percentparser = CustomParser(Float64) do str, i, len, opts
        num, ii = tryparsenext(floatparser, str, i, len, opts)
        if isnull(num)
            return num, ii
        else
            # parse away the % char
            ii = eatwhitespaces(str, ii, len)
            c, k = next(str, ii)
            if c != '%'
                return Nullable{Float64}(), ii # failed to parse %
            else
                return num, k # the point after %
            end
        end
    end

    @test tryparsenext(percentparser, "10%")  |> unwrap == (10.0, 4)
    @test tryparsenext(percentparser, "10.32 %") |> unwrap == (10.32, 8)
    @test tryparsenext(percentparser, "2k%") |> failedat ==  2
end


text:
using OffsetArrays, ProgressMeter
import KitBase

begin
    cd(@__DIR__)
    D = KitBase.read_dict("config.txt")
    for key in keys(D)
        s = Symbol(key)
        @eval $s = $(D[key])
    end

    γ = KitBase.heat_capacity_ratio(inK, 1)
    set = KitBase.Setup(
        case,
        space,
        flux,
        collision,
        nSpecies,
        interpOrder,
        limiter,
        cfl,
        maxTime,
    )
    pSpace = KitBase.PSpace1D(x0, x1, nx, pMeshType, nxg)
    vSpace = KitBase.VSpace1D(umin, umax, nu, vMeshType, nug)
    μᵣ = KitBase.ref_vhs_vis(knudsen, alphaRef, omegaRef)
    gas =
        KitBase.Gas(knudsen, mach, prandtl, inK, γ, omega, alphaRef, omegaRef, μᵣ, mass, 0)
    wL, primL, hL, bL, bcL, wR, primR, hR, bR, bcR =
        KitBase.ib_rh(gas.Ma, gas.γ, gas.K, vSpace.u)
    ib = KitBase.IB(wL, primL, bcL, wR, primR, bcR)
    ks = KitBase.SolverSet(set, pSpace, vSpace, gas, ib, pwd())

    ctr = OffsetArray{KitBase.ControlVolumeParticle1D}(undef, eachindex(ks.pSpace.x))
    face = Array{KitBase.Interface1D}(undef, ks.pSpace.nx + 1)

    for i in eachindex(ctr)
        if i <= ks.pSpace.nx ÷ 2
            ctr[i] = KitBase.ControlVolumeParticle1D(
                ks.pSpace.x[i],
                ks.pSpace.dx[i],
                ks.ib.wL,
                ks.ib.primL,
                KitBase.vhs_collision_time(ks.ib.primL, ks.gas.μᵣ, ks.gas.ω),
            )
        else
            ctr[i] = KitBase.ControlVolumeParticle1D(
                ks.pSpace.x[i],
                ks.pSpace.dx[i],
                ks.ib.wR,
                ks.ib.primR,
                #ks.ib.wL,
                #ks.ib.primL,
                KitBase.vhs_collision_time(ks.ib.primR, ks.gas.μᵣ, ks.gas.ω),
            )
        end
    end

    for i = 1:ks.pSpace.nx+1
        face[i] = KitBase.Interface1D(ks.ib.wL)
    end

    ptc = KitBase.init_ptc(ks, ctr)
    ptc_new = deepcopy(ptc)
    ks.gas.np = length(ptc) ÷ 2
end

t = 0.0
dt = KitBase.timestep(ks, ctr, t)
res = zeros(3)

function solve!(ks, ctr, ptc, ptc_new, t, dt, res, nt = 10)
    @showprogress for iter = 1:nt
        #for iter in 1:10
        #KitBase.reconstruct!(ks, ctr)

        @inbounds Threads.@threads for i = 1:ks.pSpace.nx+1
            KitBase.flux_equilibrium!(
                face[i].fw,
                ctr[i-1].w,
                ctr[i].w,
                ks.gas.K,
                ks.gas.γ,
                ks.gas.μᵣ,
                ks.gas.ω,
                ks.gas.Pr,
                dt,
                0.5 * ctr[i-1].dx,
                0.5 * ctr[i].dx,
            )
        end

        #KitBase.update!(ks, ctr, ptc, ptc_new, face, dt, res; coll = :bgk, bc = :fix)

        KitBase.particle_transport!(ks, ctr, ptc, ptc_new, dt)
        KitBase.particle_collision!(ks, ctr, ptc_new, face, res, :bgk)
        KitBase.particle_boundary!(ks, ctr, ptc, ptc_new, face, dt, :bgk, :fix)
        KitBase.duplicate!(ptc, ptc_new, ks.gas.np)

        t += dt
    end
    @show res
    return t
end

t = solve!(ks, ctr, ptc, ptc_new, t, dt, res, 2)

KitBase.plot_line(ks, ctr)

#for i = 1:1
KitBase.particle_transport!(ks, ctr, ptc, ptc_new, dt)
KitBase.particle_collision!(ks, ctr, ptc_new, face, res, :bgk)
KitBase.particle_boundary!(ks, ctr, ptc, ptc_new, face, dt, :bgk, :fix)
KitBase.duplicate!(ptc, ptc_new, ks.gas.np)
#end

begin
    using Plots
    pltx = ks.pSpace.x[1:ks.pSpace.nx]
    plty = zeros(ks.pSpace.nx, 6)
    for i in eachindex(pltx)
        for j = 1:2
            plty[i, j] = ctr[i].wf[j] .+ ctr[i].wp[j]
        end
        plty[i, 3] = ctr[i].wf[end] + ctr[i].wp[end]
    end
    plot(pltx, plty[:, 1:3], lw = 2, legend = :none)
end


text:
# Function that reset the forces and partial vectors used for parallel computation

function resetForces!(f :: Forces)
  @. f.f = 0.
  for i in 1:Threads.nthreads()
    @. f.fpartial[i] = 0.
  end
  @. f.upartial = 0.
  @. f.nenc_partial = 0.
end


text:
mutable struct H5Seis
  seis::PyCall.PyObject
  # types must be compatible with h5geo seis IO data
  ind::Array{UInt64,1}
  data::Array{Float32,2}
end


text:
####### introduction
# the Julia code to run the simulations and produce the plots in Lewis Lehe's paper Downtown Tolling and the Distribution of Trip Lengths, published in the journal Economics of Transportation
# if you have questions email Lewis Lehe at lewis500@berkeley.edu
# in addition to using the julia language, a lot of the syntax comes from the matplotlib module from python. Julia is able to call this function with similar commands to those use in python.

######BEGIN
###### Here we defined the module Simulation, which contains all the code required to run our simulation. The functions are placed in a module so that you don't have to recompile every time you edit something.
module Simulation
using PyCall
using PyPlot
using LaTeXStrings
using Distributions
using Roots:fzero
using Optim:optimize,minimizer
using DataFrames
using Cubature
using QuadGK:quadgk

#the World object that holds the parameters for the simulation. a world object is passed into all the other functions we use. For example, the p(k) function (defined like in the paper) needs to know k0 and pf (k_0 and p_f in the paper)
struct World
  λ::Float64
  pf::Float64
  k0::Float64
  μ::Array{Float64}
  Σ::Array{Float64}
  z::MultivariateDistribution #this is produced from Σ and μ
end
#the generator function for World objects. All it does is take the means and covariance matrix used to produce the population, and creates a bivariate lognormal distribution for it as the z parameter in the World object
function World(;
  λ::Float64 = 20.,
  pf::Float64 = 2.2,
  k0::Float64 = 55.,
  μ::Array{Float64} = [1.0; 1.0],
  Σ::Array{Float64} = [.1 .1; .1 .1]
  )
  return World(λ,pf,k0,μ,Σ,MvLogNormal(μ, Σ))
end

# Equilibrium is defined as type consisting of a density, a pace and a circulation
Equilibrium = Tuple{Float64, Float64, Float64}

# pk and qk are physical functions used in our calculations. pk is p(k) in the paper. qk is q(k) in the paper. these are the mfd functions
pk(w::World, k::Float64)::Float64 = w.pf*exp(.5(k/w.k0)^2)
qk(w::World, k::Float64)::Float64 = k/pk(w,k)

# X() is as defined as in the paper: the disutility of a car trip when pace is p, trip length is l and the given tolls
X(p::Float64,l::Float64;τa::Float64=0.,τd::Float64=0.)::Float64 = (p+τd)*l + τa

##### the following block gives all the functions required to calculate the equilbrium.
##### each takes a "World" (a set of parameters), a pace p and the tolls given by obvious keyword arguments

# qdp gives the cirulation depanded for a given p and set of tolls.
function qdp(w::World,p::Float64;τa::Float64=0.,τd::Float64=0.)
  function f(x)::Float64
    l::Float64 = x[2]/(1-x[2])
    ɛ::Float64 = (p+τd)*l + τa + x[1]/(1-x[1])
    return l*pdf(w.z,[ɛ;l])/((1-x[2])*(1-x[1]))^2
  end
  return w.λ*hcubature(f,[0.,0.],[1.,1.];abstol=1e-4)[1]
end

# adp gives the arrival flow demanded for a given p and set of tolls
function adp(w::World,p::Float64;τa::Float64=0.,τd::Float64=0.)::Float64
  function f(x)::Float64
      l::Float64 = x[2]/(1-x[2])
      cost::Float64 = (p+τd)*l + τa
      ɛ::Float64 = cost+ x[1]/(1-x[1])
      pdf(w.z,[ɛ;l])/((1-x[2])*(1-x[1]))^2
  end
  return w.λ*hcubature(f,[0.,0.],[1.,1.];abstol=1e-4)[1]
end

# TSS gives the total social surplus for a given p and set of tolls
function TSS(w::World,p::Float64;τa::Float64=0.,τd::Float64=0.)::Float64
  function f(x)::Float64
    l::Float64 = x[2]/(1-x[2])
    cost::Float64 = (p+τd)*l + τa
    ɛ::Float64 = cost + x[1]/(1-x[1])
    pdf(w.z,[ɛ;l])*(ɛ-l*p)/((1-x[2])*(1-x[1]))^2
  end
  return w.λ*hcubature(f,[0.,0.],[1.,1.];abstol=1e-3)[1]
end

# CS gives the consumer surplus for a given p and set of tolls
function CS(w::World,p::Float64;τa::Float64=0.,τd::Float64=0.)::Float64
  function f(x)::Float64
    l::Float64 = x[2]/(1-x[2])
    cost::Float64 = (p+τd)*l + τa
    ɛ::Float64 = cost + x[1]/(1-x[1])
    pdf(w.z,[ɛ;l])*(ɛ-cost)/((1-x[2])*(1-x[1]))^2
  end
  return w.λ*hcubature(f,[0.,0.],[1.,1.];abstol=1e-3)[1]
end

#### the functions chooseA and chooseT choose the welfare-maximizing access and trip toll levels
function chooseA(w::World)::Float64
  optimize(τd::Float64->begin
    p= getEq(w;τd=τd)[2]
    -TSS(w,p;τd=τd)
  end,0.,4.0,
  ;abs_tol=1e-4)|>minimizer
end

function chooseT(w::World)::Float64
  optimize(τa::Float64->begin
    -TSS(w,getEq(w;τa=τa)[2];τa=τa)
  end,1.,10.,
  ;abs_tol=1e-4)|>minimizer
end

##### the getEq function finds the equilibrium given tolls for a given "world"
function getEq(w::World;τa::Float64=0.,τd::Float64=0.)::Equilibrium
  ke::Float64 = fzero(
    k::Float64-> qk(w,k) - qdp(w,pk(w,k); τa=τa,τd=τd)
    ,w.k0;order=2)
  return (ke, pk(w,ke), qk(w,ke))
end

###### the plotSupplyAndDemand function takes a "World" and plots the q_d and p_s(q) functions in q/p space (the walters diagram)
function plotSupplyAndDemand(w::World)
  kRange = linspace(w.k0*.2,w.k0*1.6, 30)
  pRange = map(k->pk(w,k),kRange)
  qsRange = map(k->qk(w,k),kRange)
  k,p,q = getEq(w)
  close("all")
  fig,ax = subplots(figsize=(4,3))
  # supply curve
  ax[:plot](qsRange,pRange,label=L"p_s/\tilde{p}_s", ls="-",color="royalblue")
  # untolled
  ax[:plot](map(p->qdp(w,p),pRange),pRange,label=L"p_d(q;0,0)",ls="-.",color="coral")
  # distance toll
  τd = chooseA(w)
  ax[:plot](map(p->qdp(w,p;τd=τd,τa=0.),pRange),pRange,label=L"p_d(q;\tau_d^*,0)",ls=":",color="green")
  # access toll
  τa = chooseT(w)
  ax[:plot](map(p->qdp(w,p;τa=τa,τd=0.),pRange),pRange,label=L"p_d(q;0,\tau_a^*)",ls="--",color="red")
  ax[:set_xlim]((0,qk(w,w.k0)*1.2))
  ax[:set_ylim]((0,pRange[end]))
  ax[:set_xlabel](L"q")
  ax[:set_ylabel](L"p")
  ax[:legend](frameon=false)
  fig[:tight_layout]()
  plt[:show]()
  return fig
end

##### plotTripLengths takes a "World," find the optimal access and distance tolls and then plots the trip length distribution
function plotTripLengths(w)
  k,p,q = getEq(w)
  τd = chooseA(w)
  kd,pd,qd = getEq(w;τd=τd)
  τa = chooseT(w)
  kt,pt,qt = getEq(w;τa = τa)

  lRange = linspace(0,mean(w.z)[2]*2,100)

  function plotEq(w::World,ax;label="o",ls="-",τa=0.,τd=0.,color="blue")
    k,p,q = getEq(w;τa=τa,τd=τd)
    mass = map(l->begin
      quadgk(x-> begin
        ɛ = X(p,l;τa=τa,τd=τd) + x/(1-x)
        pdf(w.z,[ɛ;l])
      end,0.,1.)[1]
    end,lRange)
    ax[:plot](lRange,mass,label=label,ls=ls,color=color)
  end

  fig,ax = subplots(figsize=(4,3))
  ax[:set_xlabel]("km")
  ax[:set_ylabel]("density")
  plotEq(w,ax;label=L"g(l;0,0)",ls="-.",color="coral")
  plotEq(w,ax;label=L"g(l;\tau_d^*,0)",τd=τd,ls=":",color="green")
  plotEq(w,ax;label=L"g(l;0,\tau_a^*)",τa=τa,ls="--",color="red")
  ax[:legend](frameon=false)
  fig[:tight_layout]()
  fig
end

##### createTable creates the table used in the simulation by finding the optimal toll values
function createTable(w::World)
  k,p,q = getEq(w)
  a = adp(w,p)
  l = q/a
  tss = TSS(w,p)
  cs = CS(w,p)
  tr = tss - cs

  τd = chooseA(w)
  kd,pd,qd = getEq(w;τd=τd)
  ad = adp(w,pd;τd=τd)
  tssd = TSS(w,pd;τd=τd)
  csd = CS(w,pd;τd=τd)
  ld = qd/ad
  trd = tssd - csd

  τa = chooseT(w)
  kt,pt,qt = getEq(w;τa = τa)
  at = adp(w,pt;τa=τa)
  tsst = TSS(w,pt;τa=τa)
  cst = CS(w,pt;τa=τa)
  lt = qt/at
  trt = tsst - cst

  # a helper function to round all the results to a few decimal places
  rounder(arr::Array{Float64})::Array{Float64} = map(x-> round(x,1) , arr)
  return DataFrame(
    toll=["none","distance","access"],
    τ=[0,τd,τa]|>rounder,
    k=[k,kd,kt]|>rounder,
    l=[l,ld,lt]|>rounder,
    p=[p,pd,pt]|>rounder,
    q=[q,qd,qt]|>rounder,
    a=[a,ad,at]|>rounder,
    CS=[cs,csd,cst]|>rounder,
    TR=[tr,trd,trt]|>rounder,
    TSS=[tss,tssd,tsst]|>rounder
    )
end
#plot the distribution of trip lengths in the UE, with the TSS-maximizing access toll and with the TSS-maximizing distance toll
function plotDivision(w::World)
  n = 30
  xMax = 1.5*mean(w.z)[2]
  yMax = 1.7*mean(w.z)[1]
  X = linspace(0,xMax,n)
  Y = linspace(0,yMax,n)
  Z = zeros(Float64,(n,n))
  for i=1:n,j=1:n
    Z[j,i] = pdf(w.z,[Y[i];X[j]])
  end
  plt[:close]("all")
  fig,ax = subplots(figsize=(4,3))
  ax[:contourf](X,Y,Z,cmap=:Blues)
  k,p,q = getEq(w)
  ax[:plot](X, map(l->l*p, X),label=L"D_{untolled}",ls="-.",c="coral",lw=2)
  τd = chooseA(w)
  kd,pd,qd = getEq(w;τd=τd)
  ax[:plot](X, map(l->l*(pd + τd), X),label=L"D_{\tau_d^*}",ls=":",c="g",lw=2)
  τa = chooseT(w)
  kt,pt,qt = getEq(w;τa = τa)
  ax[:plot](X, map(l->l*pt + τa, X),label=L"D_{\tau_a^*}",ls="--",c="r",lw=2)
  ax[:set_xlim](0,xMax)
  ax[:set_ylim](0,yMax)
  ax[:legend]()
  ax[:set_xlabel](L"l")
  ax[:set_ylabel](L"\varepsilon")
  ax[:text](0.5, 5.0,L"R_1")
  ax[:text](4.0,19.0,L"R_2")
  fig[:tight_layout]()
  plt[:show]()
  return fig
end

#### plotAdvantage creates the figure showing the advantage of the distance toll over the access toll for various values of covariance between benefit and trip length
function plotAdvantage(w::World)
  n = 20
  cvRange = linspace(-.15,.15,n)
  res = map(cv->begin
    w = World(;
      μ=w.μ,
      Σ=[w.Σ[1,1] cv; cv w.Σ[2,2]],
      λ=w.λ,pf=w.pf,k0=w.k0)
    τd = chooseA(w)
    kd,pd,qd = getEq(w;τd=τd)
    tssd = TSS(w,pd;τd=τd)
    τa = chooseT(w)
    kt,pt,qt = getEq(w;τa=τa)
    tsst = TSS(w,pt;τa=τa)
    tssd/tsst-1
  end,cvRange)
  fig,ax = subplots(figsize=(4,3))
  ax[:plot](cvRange,res*100)
  ax[:set_xlabel](L"\sigma^2_{\varepsilon,l}")
  ax[:set_ylabel]("% welfare gain from distance toll")
  fig[:tight_layout]()
  return fig
end
#end of the Simulation Module
end
###### in this part of the code you define the free parameters and run whichever functions you want for the plot and table, saving them to disk in the images folder

#here, world is a "World" variable holding the parameters used in the paper. you can modify it with your own parameters.
world = Simulation.World(;
  μ=[2.4,1.0],
  Σ=[0.2 0.12; 0.12 0.2],
  λ=20.,
  pf=2.2,
  k0=55.)

# print the table
Simulation.createTable(world)

# make the plots
fig = Simulation.plotTripLengths(world)
fig[:savefig]("img/tripLengths.pdf")

# fig2 = Simulation.plotAdvantage(world)
# fig2[:savefig]("img/advantage.pdf")

fig3 = Simulation.plotSupplyAndDemand(world)
fig3[:savefig]("img/supplyDemand.pdf")

fig4 = Simulation.plotDivision(world)
fig4[:savefig]("img/division.pdf")


text:
## return copy of `a` with midi channel changed to `v`
function channel(v::Int, a::Atom)
	@assert 1 <= v <= 16
	if in(:cha,fieldnames(typeof(a)))
		a = deepcopy(a)
		a.cha = v
	end
	a
end

## return copy of `a` with duration and offset multiplied by `v`
function dilate(v::Rat, a::Atom)
	@assert v > 0
	if in(:dur,fieldnames(typeof(a)))
		a = deepcopy(a)
		a.ofs *= v
		a.dur *= v
	end
	a
end

## return copy of `a` with `v` added to offset
function sshift(v::Rat, a::Atom)
	a = deepcopy(a)
	a.ofs += v
	a
end

## return copy of `a` with velocity mutiplied by `v`
function accel(v::Float64, a::Atom)
	@assert 0 < v
	if in(:vel,fieldnames(typeof(a)))
		a = deepcopy(a)
		a.vel *= v
	end
	a
end

## return copy of `a` with `v` added to note inteval value
function transl(v::Int, a::Atom)
	if in(:itv,fieldnames(typeof(a)))
		a = deepcopy(a)
		a.itv += v
	end
	a
end

## return copy of `a` with `v` added to note octave value
function octave(v::Int, a::Atom)
	if in(:ocv,fieldnames(typeof(a)))
		a = deepcopy(a)
		a.ocv = v
	end
	a
end

## return copy of `a` with note scale changed to `v`
function sscale(v::Scales.Scale, a::Atom)
	if in(:sca,fieldnames(typeof(a)))
		a = deepcopy(a)
		a.sca = v
	end
	a
end

## tests if atom default fields have valid values
function atomTest(a::Atom)
	if in(:cha,fieldnames(typeof(a))) ; @assert 1 <= a.cha <= 16 ; end
	if isa(a,Duratom)
		if in(:dur,fieldnames(typeof(a))) ; @assert 0 < a.dur ; end
	else
		if in(:dur,fieldnames(typeof(a))) ; @assert 0 <= a.dur ; end
	end
	if in(:ocv,fieldnames(typeof(a))) ; @assert 1 <= a.ocv <= 5 ; end
	if in(:vel,fieldnames(typeof(a))) ; @assert 0 < a.vel ; end
	a
end

## play a note (on and off midi events)
type Note <: Duratom
	ofs::Rat			# offset
	dur::Rat			# duration
	cha::Int			# channel
	itv::Int			# interval value
	ocv::Int			# octave
	vel::Float64		# velocity
	sca::Scales.Scale	# scale

	Note(ofs,dur,cha,itv,ocv,vel,sca) = atomTest(new(ofs,dur,cha,itv,ocv,vel,sca))
end

## select instrument bank
type BankSel <: Atom
	ofs::Rat			# offset
	dur::Rat			# duration
	cha::Int			# channel
	bank::Int			# bank number

	BankSel(ofs,cha,bank) = ( @assert 0<=bank<=16383 ; atomTest(new(ofs,0//1,cha,bank)) )
end

## select instrument program
type ProgSel <: Atom
	ofs::Rat			# offset
	dur::Rat			# duration
	cha::Int			# channel
	prog::Int			# program

	ProgSel(ofs,cha,prog) = ( @assert 0<=prog<=127 ; atomTest(new(ofs,0//1,cha,prog)) )
end

## set channel volume
type VolSet <: Atom
	ofs::Rat			# offset
	dur::Rat			# duration
	cha::Int			# channel
	vol::Float64		# volume [0,1]

	VolSet(ofs,cha,vol) = ( @assert 0<=vol<=1 ; atomTest(new(ofs,0//1,cha,vol)) )
end

## pitch wheel bend
type PitchWheel <: Atom
	ofs::Rat			# offset
	dur::Rat			# duration
	cha::Int			# channel
	val::Int			# value {0..2^14-1}

	PitchWheel(ofs,cha,val) = ( @assert 0<=val<=16383 ; atomTest(new(ofs,0//1,cha,val)) )
end

## channel aftertouch
type ChanAfter <: Atom
	ofs::Rat			# offset
	dur::Rat			# duration
	cha::Int			# channel
	val::Int			# value {0..127}

	ChanAfter(ofs,cha,val) = ( @assert 0<=val<=127 ; atomTest(new(ofs,0//1,cha,val)) )
end

## channel contol 7bit
type Control7 <: Atom
	ofs::Rat			# offset
	dur::Rat			# duration
	cha::Int			# channel
	num::Int			# control number {64..127}
	val::Int			# control value {0..127}

	Control7(ofs,cha,num,val) = ( @assert 0<=val<=127 && 64<=num<=127 ; atomTest(new(ofs,0//1,cha,num,val)) )
end

## channel contol 14bit
type Control14 <: Atom
	ofs::Rat			# offset
	dur::Rat			# duration
	cha::Int			# channel
	num::Int			# control number {0..31}
	val::Int			# control value {0..2^14-1}

	Control14(ofs,cha,num,val) = ( @assert 0<=val<=16383 && 0<=num<=31 ; atomTest(new(ofs,0//1,cha,num,val)) )
end

## track marker
type Marker <: Atom
	ofs::Rat			# offset
	dur::Rat			# duration
	name::AbstractString	# marker name

	Marker(ofs,name) = ( @assert length(name) > 0 ; atomTest(new(ofs,0//1,name)) )
end

## set track time signature
type TimeSignature <: Atom
	ofs::Rat			# offset
	dur::Rat			# duration
	num::Int			# numerator
	den::Int			# denominator as logarithm of 2
	clocks::Int			# divisions per quarter note
	sub::Int			# subdivisions

	TimeSignature(ofs,num,den,clocks,sub) = ( @assert 1<=num<=32 && 0<=den<=6 && 1<=clocks<64 && 1<=sub<64 ; atomTest(new(ofs,0//1,num,den,clocks,sub)) )
end

## track tempo
type Tempo <: Atom
	ofs::Rat			# offset
	dur::Rat			# duration
	bpm::Int			# beats per minute

	Tempo(ofs,bpm) = ( @assert 12<=bpm<=1728 ; atomTest(new(ofs,0//1,bpm)) )
end

## functions for writing Atoms to a MidiTrack; uses Jumidi functions; returns a's end time in MidiTrack ticks

function toTrack!(mt::MidiTrack, a::Note)
	t1 = timm(mt,a.ofs) ; t2 = t1 + tim(mt,a.dur) ; @assert t1 < t2-1
	ve = ftoi7(a.vel) ; pch::Int = Scales.getPitch(a.sca,a.ocv,a.itv) ; @assert 0 <= pch < 128
	noteOn!(mt,a.ofs,a.cha,pch,ve)
	noteOff!(mt,a.ofs+a.dur,a.cha,pch,ve)
	t2
end

toTrack!(mt::MidiTrack, a::BankSel) = bankSel!(mt,a.ofs,a.cha,a.bank)

toTrack!(mt::MidiTrack, a::ProgSel) = progSel!(mt,a.ofs,a.cha,a.prog)

toTrack!(mt::MidiTrack, a::VolSet) = control14!(mt,a.ofs,a.cha,7,ftoi14(a.vol))

toTrack!(mt::MidiTrack, a::PitchWheel) = pitchWheel!(mt,a.ofs,a.cha,a.val)

toTrack!(mt::MidiTrack, a::ChanAfter) = chanAfter!(mt,a.ofs,a.cha,a.val)

toTrack!(mt::MidiTrack, a::Control7) = control7!(mt,a.ofs,a.cha,a.num,a.val)

toTrack!(mt::MidiTrack, a::Control14) = control14!(mt,a.ofs,a.cha,a.num,a.val)

toTrack!(mt::MidiTrack, a::Marker) = trackMarker!(mt,a.ofs,a.name)

toTrack!(mt::MidiTrack, a::TimeSignature) = timeSigSet!(mt,a.ofs,a.num,a.den,a.clocks,a.sub)

toTrack!(mt::MidiTrack, a::Tempo) = tempoSet!(mt,a.ofs,a.bpm)

function Base.isequal(x::Atom, y::Atom)
	if typeof(x) != typeof(y) ; return false ; end
	for t in fieldnames(typeof(x))
		v1 = getfield(x,t)
		v2 = getfield(y,t)
		if !isequal(v1,v2) ; return false ; end
	end
	true
end


text:
using MCTSParser
using JLD2, FileIO

if false
    tree = GameTree(TicTacToe(), MCTSParser.NN())
    train(tree, 300, 200, 100)
    save("tree.jld2", "tree", tree)
else
    tree = load("tree.jld2", "tree")
    nodes = selfplay(tree, 100)
    for node in nodes
        println(node.state)
        println(node.p.data)
        println(node.v.data)
    end
    println("finish")
end


text:
using Revise
using HSVI
using POMDPs
using POMDPPolicies
using BeliefUpdaters
using POMDPModelTools
using POMDPModels

pomdp = TigerPOMDP()

spomdp = SparseTabularPOMDP(pomdp)
b0 = POMDPPolicies.beliefvec(pomdp, n_states(pomdp), initialstate(pomdp))
solver = HSVISolver(Vector(b0), 1e-3, 1e-3, SawtoothProjection())

policy = solve(solver, spomdp)

using ProfileView
using Profile

Profile.clear()

@profile solve(solver, spomdp)

ProfileView.view()


text:
using Plots, Colors, ColorSchemes, FlagColors
function displaycolorscheme(k,cs)
    heatmap(repeat(collect(range(0,1,length=20))',2,1),
        showaxis=false, legend=false, ticks=false,
        #size=(100,30),
        color = cgrad(cs),
        title = k,
    )
end

examples = [:fi, :de, :se, :cn, :ca, :us, :gb, :fr, :ag, :eg, Symbol("us-ok"), :dk]
pl = [displaycolorscheme(k,flagcolors[k]) for k in examples]
savefig(plot(pl...), "doc/examples.png")


text:
# This file (`LibWasmer.jl` not `prologue.jl`) is automatically generated using
# Clang.jl and should not be edited manually. Take a look at the `gen/` folder
# if something is to be changed.

using Pkg.Artifacts
using Pkg.BinaryPlatforms

tripletnolibc(platform) = replace(triplet(platform), "-gnu" => "")
get_libwasmer_location(platform) =
    "wasmer-v$release_version-$(tripletnolibc(platform))-c-api"

function get_libwasmer_location()
    artifact_info = artifact_meta("libwasmer", joinpath(@__DIR__, "..", "Artifacts.toml"))
    artifact_info === nothing && return nothing

    parent_path = artifact_path(Base.SHA1(artifact_info["git-tree-sha1"]))
    return joinpath(
        parent_path,
        "lib/libwasmer"
    )
end

const libwasmer_key = "LIBWASMER_LOCATION"
const libwasmer = haskey(ENV, libwasmer_key) ?
    ENV[libwasmer] : get_libwasmer_location()


text:
import SpecialFunctions

function flux_Gxp(Gxp, nx, ny, u1, u2, rho, pr)
    tx = ny
    ty = -nx

    ut = u1*tx + u2*ty
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr
    S1 = ut*sqrt(beta)
    B1 = 0.5*exp(-S1*S1)/sqrt(pi*beta)
    A1pos = 0.5*(1 + SpecialFunctions.erf(S1))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un

    Gxp[1] = (rho*(ut*A1pos + B1))

    temp1 = pr_by_rho + ut*ut
    temp2 = temp1*A1pos + ut*B1
    Gxp[2] = (rho*temp2)

    temp1 = ut*un*A1pos + un*B1
    Gxp[3] = (rho*temp1)

    temp1 = (7*pr_by_rho) + u_sqr
    temp2 = 0.5*ut*temp1*A1pos
    temp1 = (6*pr_by_rho) + u_sqr
    Gxp[4] = (rho*(temp2 + 0.5*temp1*B1))
    return nothing
end

function flux_Gxn(Gxn, nx, ny, u1, u2, rho, pr)

    tx = ny
    ty = -nx

    ut = u1*tx + u2*ty
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr
    S1 = ut*sqrt(beta)
    B1 = 0.5*exp(-S1*S1)/sqrt(pi*beta)
    A1neg = 0.5*(1 - SpecialFunctions.erf(S1))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un

    Gxn[1] = (rho*(ut*A1neg - B1))

    temp1 = pr_by_rho + ut*ut
    temp2 = temp1*A1neg - ut*B1
    Gxn[2] = (rho*temp2)

    temp1 = ut*un*A1neg - un*B1
    Gxn[3] = (rho*temp1)

    temp1 = (7*pr_by_rho) + u_sqr
    temp2 = 0.5*ut*temp1*A1neg
    temp1 = (6*pr_by_rho) + u_sqr
    Gxn[4] = (rho*(temp2 - 0.5*temp1*B1))
    return nothing
end

function flux_Gyp(Gyp, nx, ny, u1, u2, rho, pr)
    tx = ny
    ty = -nx

    ut = u1*tx + u2*ty
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr
    S2 = un*sqrt(beta)
    B2 = 0.5*exp(-S2*S2)/sqrt(pi*beta)
    A2pos = 0.5*(1 + SpecialFunctions.erf(S2))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un

    Gyp[1] = (rho*(un*A2pos + B2))

    temp1 = pr_by_rho + un*un
    temp2 = temp1*A2pos + un*B2

    temp1 = ut*un*A2pos + ut*B2
    Gyp[2] = (rho*temp1)

    Gyp[3] = (rho*temp2)

    temp1 = (7*pr_by_rho) + u_sqr
    temp2 = 0.5*un*temp1*A2pos
    temp1 = (6*pr_by_rho) + u_sqr
    Gyp[4] = (rho*(temp2 + 0.5*temp1*B2))
    return nothing
end

function flux_Gyn(Gyn, nx, ny, u1, u2, rho, pr)
    tx = ny
    ty = -nx

    ut = u1*tx + u2*ty
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr
    S2 = un*sqrt(beta)
    B2 = 0.5*exp(-S2*S2)/sqrt(pi*beta)
    A2neg = 0.5*(1 - SpecialFunctions.erf(S2))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un

    Gyn[1] = (rho*(un*A2neg - B2))

    temp1 = pr_by_rho + un*un
    temp2 = temp1*A2neg - un*B2

    temp1 = ut*un*A2neg - ut*B2
    Gyn[2] = (rho*temp1)

    Gyn[3] = (rho*temp2)

    temp1 = (7*pr_by_rho) + u_sqr
    temp2 = 0.5*un*temp1*A2neg
    temp1 = (6*pr_by_rho) + u_sqr
    Gyn[4] = (rho*(temp2 - 0.5*temp1*B2))
    return nothing
end

function flux_Gx(Gx, nx, ny, u1, u2, rho, pr)
    tx = ny
    ty = -nx

    ut = u1*tx + u2*ty
    un = u1*nx + u2*ny

    Gx[1] = rho*ut

    Gx[2] = pr + rho*ut*ut

    Gx[3] = rho*ut*un

    temp1 = 0.5*(ut*ut + un*un)
    rho_e = 2.5*pr + rho*temp1
    Gx[4] = (pr + rho_e)*ut
    return  nothing
end

function flux_Gy(Gy, nx, ny, u1, u2, rho, pr)
    tx = ny
    ty = -nx

    ut = u1*tx + u2*ty
    un = u1*nx + u2*ny

    Gy[1] = rho*un

    Gy[2] = rho*ut*un

    Gy[3] = pr + rho*un*un

    temp1 = 0.5*(ut*ut + un*un)
    rho_e = 2.5*pr + rho*temp1
    Gy[4] = (pr + rho_e)*un
    return nothing
end

# print(flux_Gy(flux_Gxn(1, 3, 2, 2, 2, 1),5,2,3,4,3,2))

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
######### linear regression 1000 obs x 10 var  ###########

using DataFrames, MCMC

# simulate dataset
srand(1)
n = 1000
nbeta = 10 # number of covariates
X = [ones(n) randn((n, nbeta-1))]
beta0 = randn((nbeta,))
Y = X * beta0 + randn((n,))

# define model
ex = quote
	vars ~ Normal(0, 1.0)  # Normal prior, std 1.0 for predictors
	resid = Y - X * vars
	resid ~ Normal(0, 1.0)
end

m = model(ex, vars=zeros(nbeta))

# run random walk metropolis (10000 steps, 1000 for burnin, thinning 10), no adaptation
mcchain01 = run(m * RWM(0.05) * SerialMC(10000:10:100000))

acceptance(mcchain01) # ~ 3%, too low

# with adaptation (target acceptance = 30%)
mcchain02 = run(m * RAM(1., 0.3) * SerialMC(10000:10:100000))

acceptance(mcchain02) # ~ 29.7%, Ok

[mean(mcchain02) beta0] # mean sample vs original coefs


text:
convert{S<:Sculpt,Q<:Qualia,C<:Clay}(::Type{Flex{S,Q,Float64}}, x::Flex{S,Q,C}) = Flex{S,Q,Float64}(convert(Float64,x.lo),convert(Float64,x.hi))

convert{S<:Sculpt,Q<:Qualia,C<:Clay}(::Type{Flex{CLCL,Q,C}}, x::Flex{S,Q,C}) = Flex{CLCL,Q,C}(closed(x)...)
convert{S<:Sculpt,Q<:Qualia,C<:Clay}(::Type{Flex{CLOP,Q,C}}, x::Flex{S,Q,C}) = Flex{CLCL,Q,C}(clopened(x)...)
convert{S<:Sculpt,Q<:Qualia,C<:Clay}(::Type{Flex{OPCL,Q,C}}, x::Flex{S,Q,C}) = Flex{CLCL,Q,C}(opclosed(x)...)
convert{S<:Sculpt,Q<:Qualia,C<:Clay}(::Type{Flex{OPOP,Q,C}}, x::Flex{S,Q,C}) = Flex{CLCL,Q,C}(opened(x)...)

function convert{S<:Sculpt,Q<:Qualia,C<:Clay}(::Type{Flex{S,Q,C}}, lo::Real, hi::Real)
   low = convert(C, lo)
   hig = convert(C, hi)
   low,hig = minmax(low,hig)
   Flex{S,Q,C}(low,hig)
end

function convert{S<:Sculpt,Q<:Qualia,C<:Clay}(::Type{Tuple{C,C}}, x::Flex{S,Q,C})
   lo,hi = (Q==CLCL) ? value(x) : value(opened(x))
   lo = convert(C,lo)
   hi = convert(C,hi)
   (lo,hi)
end

function convert{S<:Sculpt,Q<:Qualia,C<:Clay}(::Type{Array{C,1}}, x::Flex{S,Q,C})
   lo,hi = (Q==CLCL) ? value(x) : value(opened(x))
   lo = convert(C,lo)
   hi = convert(C,hi)
   [lo,hi]
end

# to eliminate ambiguity
function convert{S<:Sculpt,Q<:Qualia,C<:Clay}(::Type{Bool}, x::Flex{S,Q,C})
   lo,hi = (Q==CLCL) ? value(x) : value(opened(x))
   lo = lo != zero(C)
   hi = hi != zero(C)
   if lo == hi
      lo
   else
      false
   end
end

for T in (:Int64, :Int32, :Int16)
    @eval begin
        function convert{S<:Sculpt,Q<:Qualia,C<:Clay}(::Type{$T}, x::Flex{S,Q,C})
           lo,hi = (Q==CLCL) ? value(x) : value(opened(x))
           lo = trunc($T,trunc(lo))
           hi = trunc($T,trunc(hi))
           if lo == hi
              lo
           else
              trunc($T, trunc(mostRepresentativeValue(x.lo,x.hi)))
           end
        end
    end
end

# to eliminate ambiguity
function convert{S<:Sculpt,Q<:Qualia,C<:Clay}(::Type{Integer}, x::Flex{S,Q,C})
   lo,hi = (Q==CLCL) ? value(x) : value(opened(x))
   lo = trunc(Integer,lo)
   hi = trunc(Integer,hi)
   if lo == hi
      lo
   else
      trunc(Integer, trunc(mostRepresentativeValue(x.lo,x.hi)))
   end
end

function convert{S<:Sculpt,Q<:Qualia,C<:Clay,R<:Real}(::Type{R}, x::Flex{S,Q,C})
   lo,hi = (Q==CLCL) ? value(x) : value(opened(x))
   lo = convert(R,lo)
   hi = convert(R,hi)
   if lo == hi
      lo
   else
       mostRepresentativeValue(lo,hi)
   end
end

for (fn,S) in ( (:ClCl,:CLCL), (:ClOp,:CLOP), (:OpCl,:OPCL), (:OpOp,:OPOP),
                (:clcl,:CLCL), (:clop,:CLOP), (:opcl,:OPCL), (:opop,:OPOP) )
   @eval begin
       ($fn){T<:($S),Q<:Qualia,C<:Clay}(x::Flex{T,Q,C}) = x
   end
end

for (fn,Q,S) in ( (:ClCl,:INEXACT,:CLCL), (:ClOp,:INEXACT,:CLOP), (:OpCl,:INEXACT,:OPCL), (:OpOp,:INEXACT,:OPOP),
                  (:clcl,:EXACT,:CLCL), (:clop,:EXACT,:CLOP), (:opcl,:EXACT,:OPCL), (:opop,:EXACT,:OPOP) )
   @eval begin
       convert{C<:Clay}(::Type{Flex{$S,$Q,C}}, fp::C) = Flex{$S,$Q,C}(fp,fp)
       convert{C<:Clay}(::Type{Flex{$S,$Q,C}}, lo::C, hi::C) = Flex{$S,$Q,C}(minmax(lo,hi)...)
       ($fn){C<:Clay}(fp::C) = Flex{$S,$Q,C}(fp,fp)
       ($fn){C<:Clay}(lo::C, hi::C) = Flex{$S,$Q,C}(minmax(lo,hi)...)
   end
end

# see cvtQualia.jl to access conversion logic when qualia differ

# promoting to a third type

for C in (Float64, Float32)
  for (S,T,U) in ((:Flex{CLCL,EXACT,(C)},:Flex{CLOP,INEXACT,(C)}, :Flex{CLCL,INEXACT,(C)}),
                  (:Flex{CLOP,EXACT,(C)},:Flex{CLCL,INEXACT,(C)}, :Flex{CLCL,INEXACT,(C)}),
                  (:Flex{CLCL,EXACT,(C)},:Flex{OPCL,INEXACT,(C)}, :Flex{CLCL,INEXACT,(C)}),
                  (:Flex{OPCL,EXACT,(C)},:Flex{CLCL,INEXACT,(C)}, :Flex{CLCL,INEXACT,(C)}),
                  (:Flex{CLCL,EXACT,(C)},:Flex{OPOP,INEXACT,(C)}, :Flex{CLCL,INEXACT,(C)}),
                  (:Flex{OPOP,EXACT,(C)},:Flex{CLCL,INEXACT,(C)}, :Flex{CLCL,INEXACT,(C)}),
                  (:Flex{CLOP,EXACT,(C)},:Flex{OPCL,INEXACT,(C)}, :Flex{CLCL,INEXACT,(C)}),
                  (:Flex{OPCL,EXACT,(C)},:Flex{CLOP,INEXACT,(C)}, :Flex{CLCL,INEXACT,(C)}),
                  (:Flex{CLOP,EXACT,(C)},:Flex{OPCL,INEXACT,(C)}, :Flex{CLCL,INEXACT,(C)}),
                  (:Flex{OPCL,EXACT,(C)},:Flex{CLOP,INEXACT,(C)}, :Flex{CLCL,INEXACT,(C)}),
                  (:Flex{CLOP,EXACT,(C)},:Flex{OPOP,INEXACT,(C)}, :Flex{CLOP,INEXACT,(C)}),
                  (:Flex{OPOP,EXACT,(C)},:Flex{CLOP,INEXACT,(C)}, :Flex{CLOP,INEXACT,(C)})
                  )
     @eval begin
       convert(::Type{$U}, x::($S)) = ($U)(x.lo,x.hi)
       promote_rule{::Type{$S},::Type{$T}} = ($U)
     end
   end
end


text:
@test eltype(Vector{Int}) == Int

# we use plain Identity Monads for testing

# Applicative defaults
# --------------------

struct TestApDefault
  x
end

TypeClasses.map(f, x::TestApDefault) = TypeClasses.default_map_having_ap_pure(f, x)
TypeClasses.pure(::Type{TestApDefault}, x) = TestApDefault(x)
TypeClasses.ap(f::TestApDefault, x::TestApDefault) = TestApDefault(f.x(x.x))

@test map(TestApDefault(4)) do x
  x*x
end == TestApDefault(4*4)

# Monad defaults
# --------------

struct TestDefaultFFlattenFunctor
  x
end

TypeClasses.map(f, x::TestDefaultFFlattenFunctor) = TestDefaultFFlattenFunctor(f(x.x))
TypeClasses.flatmap(f, x::TestDefaultFFlattenFunctor) = f(x.x)

@test mapn(TestDefaultFFlattenFunctor(3), TestDefaultFFlattenFunctor(4)) do x, y
  x + y
end == TestDefaultFFlattenFunctor(3 + 4)


text:
using Graphs
using SimpleGraphs
using PyPlot

function plot_deg_hist(g::GenericGraph, out_fldr::ASCIIString,fname::ASCIIString,ac::Float64)
    # degree histogram
    deglst = zeros(num_vertices(g))
    for i in vertices(g)
        deglst[i] = out_degree(i,g)
    end
    plt_deg_hist(deglst,out_fldr,fname,ac)
end

function plot_deg_hist(g::SimpleGraphs.SimpleGraph, out_fldr::ASCIIString,fname::ASCIIString,ac::Float64)
    deglst = zeros(NV(g))
    v_lst = vlist(g)
    for i=1:NV(g)
        deglst[i] = deg(g,v_lst[i])
    end
    # deglst = deg_hist(g)
    plt_deg_hist(deglst,out_fldr,fname,ac)
end

function plot_deg_list(g::SimpleGraphs.SimpleGraph, out_fldr::ASCIIString,fname::ASCIIString,ac::Float64)
    # degree histogram
    deglst = deg(g)
    plt_deg_hist(deglst,out_fldr,fname,ac)
end

function plt_deg_hist(deglst::Array{},out_fldr::ASCIIString,fname::ASCIIString,ac::Float64)
    mu = mean(deglst)
    sig = std(deglst)
    PyPlot.clf()
    plt[:hist](deglst,maximum(deglst),normed=1,alpha=0.75,facecolor="green")
    PyPlot.axvline(mu,color="k",ls="dashed")
    PyPlot.axvline(mu-sig,color="r",ls="dotted")
    PyPlot.axvline(mu+sig,color="r",ls="dotted")
    title("Degree Distribution")
    xlabel("Degree")
    ylabel("Proportion")
    annotate(@sprintf("a(G)=%0.3f",ac),xy=(0.7,0.8),fontsize=20,xycoords="axes fraction")
    annotate(@sprintf("mu=%0.3f",mu),xy=(0.7,0.7),fontsize=20,xycoords="axes fraction")
    annotate(@sprintf("std=%0.3f",sig),xy=(0.7,0.6),fontsize=20,xycoords="axes fraction")
    fname = joinpath(out_fldr,fname)
    fname = string(fname,".pdf")
    savefig(fname)
end

function plot_c_hist(c_hist::Array{},out_fldr::ASCIIString,fname::ASCIIString)
    # connectivity history
    PyPlot.clf()
    PyPlot.plot(collect(1:length(c_hist)),c_hist)
    title("Connectivity History")
    xlabel("Iteration")
    ylabel("Algebraic Connectivity")
    fname = joinpath(out_fldr,fname)
    fname = string(fname,".pdf")
    savefig(fname)
    clf()
end

function plot_c_box(c_vals::Array{},n::Int,c_opt::Float64,out_fldr::ASCIIString,fname::ASCIIString)
    # connectivity distribution vs. optimal
    PyPlot.clf()
    boxplot(c_vals)
    PyPlot.boxplot(c_vals)
    PyPlot.axhline(c_opt,color="r",ls="dotted")
    title("Optimization Performance")
    # xlabel("n")
    xticks([1,2],["Opt","Random"])
    ylabel("Algebraic Connectivity")
    fname = joinpath(out_fldr,fname)
    fname = string(fname,".pdf")
    savefig(fname)
end


text:
module Glu_ED_Selfrep_Model
using JuMP
using KNITRO

function glu_ed_selfrep_model(glc_ext_input)

    glu_ed_selfrep = Model(
        optimizer_with_attributes(
            KNITRO.Optimizer,
            "ms_enable" => 1,
            "opttol" => 1E-16,
            "opttolabs" => 1e-16,
            "honorbnds" => 1,
            "ms_maxsolves" => 10,
        ),
    )
    JuMP.set_silent(glu_ed_selfrep)

    # Bounds
    k = 1e-6

    fmax(x) = (x + sqrt(k + x^2)) / 2                   # max(0, x) approximation
    register(glu_ed_selfrep, :fmax, 1, fmax; autodiff = true)
    fmin(x) = -fmax(-x)                               # min(0, x) approximation
    register(glu_ed_selfrep, :fmin, 1, fmin; autodiff = true)

    LB_enz = 0.0     # [g enz / g DW cell]
    UB_enz = 1.0     # [g enz / g DW cell]
    LB_met = log(1e-6)  # [log(1 uM)]
    UB_met = log(10e-3) # [log(10 mM)]
    LB_v = -50       # [mmol/gDW/h]
    UB_v = 50       # [mmol/gDW/h]
    LB_dg = -100       # [kJ/mol]
    UB_dg = 100       # [kJ/mol]

    @variables glu_ed_selfrep begin
        # biomass [1/h] (g6p + ATP -> ADP)
        0 <= mu <= 10

        # Enzymes concentrations [g enz / g DW]
        LB_enz <= pts <= UB_enz # enzyme for: glc external + ATP -> ADP + g6p
        LB_enz <= ed <= UB_enz # enzyme for: g6p + 2NAD + 2ADP -> 2NADH + 2ATP + 2pep
        LB_enz <= pyk <= UB_enz # enzyme for: pep + ADP -> ATP + pyr
        LB_enz <= ldh <= UB_enz # enzyme for: pyr + NADH -> NAD + lac
        LB_enz <= ppc <= UB_enz # enzyme for: pep + CO2 -> oac + P
        LB_enz <= akgsyn <= UB_enz # enzyme for: pyr + oac -> akg
        LB_enz <= gdhm <= UB_enz # enzyme for: akg + NH4 + NADH -> glu + NAD
        LB_enz <= burn <= UB_enz # enzyme for: ATP -> ADP
        LB_enz <= nadtrdh <= UB_enz # enzyme for: NAD + NADPH -> NADH + NADP
        LB_enz <= lp <= UB_enz # enzyme for: lac internal -> lac external

        LB_enz <= ribo <= UB_enz # ribosomes

        # Metabolite concentrations [log(M)]
        LB_met <= g6p <= UB_met # glucose-6-phosphate
        LB_met <= pyr <= UB_met # pyruvate
        LB_met <= pep <= UB_met # phosphoenolpyrovate
        LB_met <= oac <= UB_met # oxalacetate
        LB_met <= akg <= UB_met # alpha-ketogluterate
        LB_met <= glu <= UB_met # glutamate
        LB_met <= lac <= UB_met # lactate
        LB_met <= atp <= UB_met # ATP
        LB_met <= adp <= UB_met # ADP
        LB_met <= nad <= UB_met # NAD
        LB_met <= nadh <= UB_met # NADH
        LB_met <= nadp <= UB_met # NADP
        LB_met <= nadph <= UB_met # NADPH

        #LB_met <= co2   <= UB_met
        #LB_met <= phos  <= UB_met # phosphate
        #LB_met <= nh4   <= UB_met

        # Fluxes [mmol/gDW/h]
        LB_v <= v_pts <= UB_v
        LB_v <= v_ed <= UB_v
        LB_v <= v_pyk <= UB_v
        LB_v <= v_ldh <= UB_v
        LB_v <= v_ppc <= UB_v
        LB_v <= v_akgsyn <= UB_v
        LB_v <= v_gdhm <= UB_v
        LB_v <= v_burn <= UB_v
        LB_v <= v_nadtrdh <= UB_v
        LB_v <= v_lp <= UB_v

        LB_v <= v_ribo <= UB_v

        # Thermodynamic variables
        LB_dg <= dg_pts <= UB_dg
        LB_dg <= dg_ed <= UB_dg
        LB_dg <= dg_pyk <= UB_dg
        LB_dg <= dg_ldh <= UB_dg
        LB_dg <= dg_ppc <= UB_dg
        LB_dg <= dg_akgsyn <= UB_dg
        LB_dg <= dg_gdhm <= UB_dg
        LB_dg <= dg_burn <= UB_dg
        LB_dg <= dg_nadtrdh <= UB_dg
        LB_dg <= dg_lp <= UB_dg

        LB_dg <= dg_ribo <= UB_dg

        # Ribosome fractions for synthesis
        0 <= used_ribo[1:11] <= 1
    end

    @NLparameters glu_ed_selfrep begin

        #R  == 8.3145e-3 # gas constants [kJ/K/mol]
        #T  == 298.15    # tempreature   [K]
        #RT == R * T
        RT == 8.3145e-3 * 298.15

        # Enzyme rates
        kcat_pts == 213.75
        kcat_ed == 268
        kcat_pyk == 58.4
        kcat_ldh == 31
        kcat_ppc == 540
        kcat_akgsyn == 4
        kcat_gdhm == 20
        kcat_burn == 10
        kcat_nadtrdh == 167.9
        kcat_lp == 100

        kcat_ribo == 80

        # Gibbs energy of reaction
        dG0_pts == -16.7
        dG0_ed == -81.8
        dG0_pyk == -31.7
        dG0_ldh == -23.7
        dG0_ppc == -40.3
        dG0_akgsyn == -60.7
        dG0_gdhm == -33.4
        dG0_burn == -57
        dG0_nadtrdh == 0
        dG0_lp == -10.3

        dG0_ribo == -20

        # media conditions
        glc_ext == log(glc_ext_input) # [log(50 mM)]
        lac_ext == log(10e-6) # [log(10 uM)]

        # capacity constraint
        total_proteome_mass_fraction == 0.26 # [g enz/gDW]

        # maintenance requirement
        min_burn_flux == 1 # [mmol/gDW/h]
    end

    @NLconstraints glu_ed_selfrep begin
        # ΔG at concentrations
        dg_pts == dG0_pts + RT * (adp + g6p - atp - glc_ext)
        dg_ed == dG0_ed + RT * (nadph + nadh + atp + pep + pyr - adp - nadp - nad - g6p)
        dg_pyk == dG0_pyk + RT * (pyr + atp - adp - pep)
        dg_ldh == dG0_ldh + RT * (lac + nad - nadh - pyr)
        dg_ppc == dG0_ppc + RT * (pep - oac)
        dg_akgsyn == dG0_akgsyn + RT * (pyr + oac - akg)
        dg_gdhm == dG0_gdhm + RT * (akg + nadh - glu - nad)
        dg_burn == dG0_burn + RT * (adp - atp)
        dg_nadtrdh == dG0_nadtrdh + RT * (nad + nadph - nadh - nadp)
        dg_lp == dG0_lp + RT * (lac_ext - lac)

        dg_ribo == dG0_ribo + RT * (ribo - glu)

        # flux bounds due to kinetics and thermodynamics
        pts * kcat_pts * fmin(exp(-dg_pts / RT) - 1) <= v_pts
        ed * kcat_ed * fmin(exp(-dg_ed / RT) - 1) <= v_ed
        pyk * kcat_pyk * fmin(exp(-dg_pyk / RT) - 1) <= v_pyk
        ldh * kcat_ldh * fmin(exp(-dg_ldh / RT) - 1) <= v_ldh
        ppc * kcat_ppc * fmin(exp(-dg_ppc / RT) - 1) <= v_ppc
        akgsyn * kcat_akgsyn * fmin(exp(-dg_akgsyn / RT) - 1) <= v_akgsyn
        gdhm * kcat_gdhm * fmin(exp(-dg_gdhm / RT) - 1) <= v_gdhm
        burn * kcat_burn * fmin(exp(-dg_burn / RT) - 1) <= v_burn
        nadtrdh * kcat_nadtrdh * fmin(exp(-dg_nadtrdh / RT) - 1) <= v_nadtrdh
        lp * kcat_lp * fmin(exp(-dg_lp / RT) - 1) <= v_lp

        ribo * kcat_ribo * fmin(exp(-dg_ribo / RT) - 1) <= v_ribo

        v_pts <= pts * kcat_pts * fmax(1 - exp(dg_pts / RT))
        v_ed <= ed * kcat_ed * fmax(1 - exp(dg_ed / RT))
        v_pyk <= pyk * kcat_pyk * fmax(1 - exp(dg_pyk / RT))
        v_ldh <= ldh * kcat_ldh * fmax(1 - exp(dg_ldh / RT))
        v_ppc <= ppc * kcat_ppc * fmax(1 - exp(dg_ppc / RT))
        v_akgsyn <= akgsyn * kcat_akgsyn * fmax(1 - exp(dg_akgsyn / RT))
        v_gdhm <= gdhm * kcat_gdhm * fmax(1 - exp(dg_gdhm / RT))
        v_burn <= burn * kcat_burn * fmax(1 - exp(dg_burn / RT))
        v_nadtrdh <= nadtrdh * kcat_nadtrdh * fmax(1 - exp(dg_nadtrdh / RT))
        v_lp <= lp * kcat_lp * fmax(1 - exp(dg_lp / RT))

        v_ribo <= ribo * kcat_ribo * fmax(1 - exp(dg_ribo / RT))

        # Ribosome fractions
        sum(used_ribo[i] for i = 1:11) == 1

        used_ribo[1] * v_ribo - mu * pts == 0
        used_ribo[2] * v_ribo - mu * ed == 0
        used_ribo[3] * v_ribo - mu * pyk == 0
        used_ribo[4] * v_ribo - mu * ldh == 0
        used_ribo[5] * v_ribo - mu * ppc == 0
        used_ribo[6] * v_ribo - mu * akgsyn == 0
        used_ribo[7] * v_ribo - mu * gdhm == 0
        used_ribo[8] * v_ribo - mu * burn == 0
        used_ribo[9] * v_ribo - mu * nadtrdh == 0
        used_ribo[10] * v_ribo - mu * lp == 0
        used_ribo[11] * v_ribo - mu * ribo == 0

        # mass balance constraints
        v_pts - v_ed == mu # g6p
        v_ed - v_pyk == 0  # pep
        v_ed + v_pyk - v_ldh == 0  # pyr
        v_ldh - v_lp == 0  # lac
        v_ppc - v_akgsyn == 0  # oac
        v_akgsyn - v_gdhm == 0  # akg
        v_gdhm - v_ribo == 0  # glu
        v_ed + v_pyk - v_pts - v_burn == mu # atp
        v_pts + v_burn - v_ed - v_pyk == -mu # adp
        v_ldh + v_gdhm - v_ed - v_nadtrdh - v_akgsyn == 0  # nad
        v_ed + v_akgsyn + v_nadtrdh - v_ldh - v_gdhm == 0  # nadh
        v_gdhm - v_akgsyn - v_nadtrdh - v_ed == 0  # nadp
        v_ed + v_akgsyn + v_nadtrdh - v_gdhm == 0  # nadph

        # density constraint(s) (can add more, e.g. membrane)
        pts + ed + pyk + ldh + ppc + akgsyn + gdhm + burn + nadtrdh + lp + ribo <=
        total_proteome_mass_fraction

        # minimum maintenance
        min_burn_flux <= v_burn

        # measured ratio constraints (not strictly necessary)
        atp == log(10) + adp  # atp/adp = 10, but remember that the concentration variables are logged
        nadh == log(0.2) + nad
        nadph == log(0.2) + nadp
    end

    @objective(glu_ed_selfrep, Max, mu)
    optimize!(glu_ed_selfrep)
    objective_value(glu_ed_selfrep)

    results_ed_selfrep = Dict(
        "mu" => value(mu),
        "pts" => value(pts),
        "ed" => value(ed),
        "pyk" => value(pyk),
        "ldh" => value(ldh),
        "ppc" => value(ppc),
        "akgsyn" => value(akgsyn),
        "gdhm" => value(gdhm),
        "burn" => value(burn),
        "nadtrdh" => value(nadtrdh),
        "lp" => value(lp),
        "ribo" => value(ribo),
        "glc_ext" => value(glc_ext),
        "g6p" => value(g6p),
        "pyr" => value(pyr),
        "pep" => value(pep),
        "oac" => value(oac),
        "akg" => value(akg),
        "lac" => value(lac),
        "glu" => value(glu),
        "atp" => value(atp),
        "adp" => value(adp),
        "nad" => value(nad),
        "nadh" => value(nadh),
        "nadp" => value(nadp),
        "nadph" => value(nadph),
        "v_pts" => value(v_pts),
        "v_ed" => value(v_ed),
        "v_pyk" => value(v_pyk),
        "v_ldh" => value(v_ldh),
        "v_ppc" => value(v_ppc),
        "v_akgsyn" => value(v_akgsyn),
        "v_gdhm" => value(v_gdhm),
        "v_burn" => value(v_burn),
        "v_nadtrdh" => value(v_nadtrdh),
        "v_lp" => value(v_lp),
        "v_ribo" => value(v_ribo),
    )
    return results_ed_selfrep
end

end #module


text:
"""
  generality(N::UnipartiteNetwork)

This function mesure the generality. The generality represents the link density (L/S) divided by the number of predators. The number of predators: (prop. of int + prop. of top)xS
"""
function generality(N::UnipartiteNetwork)
    out_degree = collect(values(degree(N; dims=1)))
    in_degree = collect(values(degree(N; dims=2)))
    degrees = DataFrame(out=out_degree, in=in_degree)
    B = nrow(degrees[degrees.out .== 0, :]) / nrow(degrees)
    T = nrow(degrees[degrees.in .== 0, :]) / nrow(degrees)
    I = nrow(degrees[(degrees.out .!= 0) .& (degrees.in .!= 0), :]) / nrow(degrees)
    L = links(N)
    Generality = L ./ (S.^2 .* (T .+ I))
    return Generality
end


text:
using SimpleTraits
using Test

const trait = SimpleTraits.trait

# @test_throws MethodError trait(4)
@test_throws ErrorException istrait(4)

# definition & adding
@traitdef Tr1{X}
@test trait(Tr1{Int})==Not{Tr1{Int}}
@test !istrait(Tr1{Int})
@traitimpl Tr1{Integer}
@test trait(Tr1{Int})==Tr1{Int}
@test istrait(Tr1{Int})
@test trait(Tr1{Bool})==Tr1{Bool}
@test trait(Tr1{AbstractString})==Not{Tr1{AbstractString}}
@test !istrait(Tr1{AbstractString})

# Logic.  trait(Tr) returns the same trait Tr if it is fulfilled and
# Not{Tr} otherwise.  This is a bit confusing.
@test trait(Tr1{AbstractString})==Not{Tr1{AbstractString}}
@test istrait(Tr1{AbstractString})==false
@test trait(Not{Tr1{AbstractString}})==Not{Tr1{AbstractString}}
@test istrait(Not{Tr1{AbstractString}})==true
@test trait(Not{Not{Tr1{AbstractString}}})==Not{Tr1{AbstractString}}
@test istrait(Not{Not{Tr1{AbstractString}}})==false
@test trait(Not{Not{Not{Tr1{AbstractString}}}})==Not{Tr1{AbstractString}}
@test istrait(Not{Not{Not{Tr1{AbstractString}}}})==true

@test trait(Not{Tr1{Integer}})==Tr1{Integer}
@test istrait(Not{Tr1{Integer}})==false
@test trait(Not{Not{Tr1{Integer}}})==Tr1{Integer}
@test istrait(Not{Not{Tr1{Integer}}})==true
@test trait(Not{Not{Not{Tr1{Integer}}}})==Tr1{Integer}
@test istrait(Not{Not{Not{Tr1{Integer}}}})==false

@traitdef Tr2{X,Y}
@test trait(Tr2{Int,AbstractFloat})==Not{Tr2{Int,AbstractFloat}}
@traitimpl Tr2{Integer, Float64}
@test trait(Tr2{Int, Float64})==Tr2{Int, Float64}
@test trait(Tr2{Int, Float32})==Not{Tr2{Int, Float32}}

# issue 9
abstract type A9 end
struct B9<:A9 end
struct C9<:A9 end
@traitdef Tr9{X}
@traitimpl Tr9{A9}
@traitimpl Not{Tr9{B9}}
@traitimpl !Tr9{C9}
@test istrait(Tr9{A9})==true
@test istrait(Tr9{B9})==false
@test istrait(Tr9{C9})==false

# issue #15 is tested through the base-traits.jl tests.

#################
# Trait functions
#################
# functions with `t` postfix are the same as previous ones but using Traitor syntax
@traitfn f(x::X) where {X; Tr1{X}} = 1  # def 1
@traitfn f(x::X) where {X; !Tr1{X}} = 2
@test f(5)==1
@test f(5.)==2

@traitfn ft(x::::Tr1) = 1  # def 1
@traitfn ft(x::::(!Tr1)) = 2
@test ft(5)==1
@test ft(5.)==2

@traitfn f(x::X,y::Y,z) where {X,Y; Tr2{X,Y}} = 1
@test f(5,5., "a")==1
@test_throws MethodError f(5,5, "a")==2
@traitfn f(x::X,y::Y,z) where {X,Y; !Tr2{X,Y}} = 2
@test f(5,5, "a")==2
# Note, two argument traits have no Traitor style syntax

# This will overwrite the definition def1 above and generate a warning
@traitfn f(x::X) where {X; !Tr2{X,X}} = 10
@traitfn f(x::X) where {X; Tr2{X,X}} = 100
@test f(5)==10
@test f(5.)==10
@traitimpl Tr2{Integer, Integer}
@test f(5.)==10
@test f(5)==100
@test f(5.)==10

# VarArg
@traitfn vara(x::X, y...) where {X; Tr1{X}} = y
@traitfn vara(x::X, y...) where {X; !Tr1{X}} = x
@test vara(5, 7, 8)==(7,8)
@test vara(5.0, 7, 8)==5.0
@traitfn vara2(x::X...) where {X; Tr1{X}} = x
@test vara2(5, 7, 8)==(5, 7, 8)
@test_throws MethodError vara2(5, 7, 8.0)

@traitfn vara3(::X...) where {X; Tr1{X}} = X
@test vara3(5, 7, 8)==Int
@test_throws MethodError vara3(5, 7, 8.0)

@traitfn varat(x::::Tr1, y...) = y
@traitfn varat(x::::(!Tr1), y...) = x
@test varat(5, 7, 8)==(7,8)
@test varat(5.0, 7, 8)==5.0
@traitfn vara2t(x::::Tr1...) = x
@test vara2t(5, 7, 8)==(5, 7, 8)
@test_throws MethodError vara2t(5, 7, 8.0)

@traitfn vara3t(::X::Tr1...) where {X} = X
@test vara3t(5, 7, 8)==Int
@test_throws MethodError vara3t(5, 7, 8.0)

# kwargs (issue 27)
@traitfn kwfn1(x::::Tr1; k=1) = x+k
@traitfn kwfn1(x::::(!Tr1); k=2) = x-k
@test kwfn1(5)==6
@test kwfn1(5.0)==3.0
@test kwfn1(5,k=2)==7
@test kwfn1(5.0,k=3)==2.0

@traitfn kwfn2(x::::Tr1, y...; k=1) = x+y[1]+k
@traitfn kwfn2(x::::(!Tr1), y...; k::Int=2) = x+y[1]-k
@test kwfn2(5,5)==11
@test kwfn2(5.0,5)==8.0
@test kwfn2(5,5,k=2)==12
@test kwfn2(5.0,5,k=3)==7.0
@test_throws TypeError kwfn2(5.0,5,k="sadf")

@traitfn kwfn3(x::::Tr1, y...; kws...) = x+y[1]+length(kws)
@traitfn kwfn3(x::::(!Tr1), y...; kws...) = x+y[1]-length(kws)
@test kwfn3(5,5)==10
@test kwfn3(5.0,5)==10
@test kwfn3(5,5,k=2)==11
@test kwfn3(5.0,5,k=3)==9
@test kwfn3(5.0,5,k=3,kk=9)==8

@traitfn kwfn4(x::::Tr1, y...; kws...) = x+y[1]+length(kws)
@test_throws ErrorException @traitfn kwfn4(x::::(!Tr1), y...) = x+y[1]-length(kws)
@traitfn kwfn5(x::::Tr1, y...) = x+y[1]+length(kws)
@test_throws ErrorException @traitfn kwfn5(x::::(!Tr1), y...; k=1) = x+y[1]-length(kws)

# Default args, issue #32
####
@traitfn defargs1(x::::Tr1, y=2) = x+y
@traitfn defargs1(x::::(!Tr1), y=2) = x-y
@test defargs1(1,3)==4
@test defargs1(1)==3
@test defargs1(1.0,4)==-3
@test defargs1(1.0)==-1

@traitfn defargs2(x::::Tr1, y=2) = x+y
@test_throws ErrorException @traitfn defargs2(x::::(!Tr1), y=3) = x-y
@traitfn defargs3(x::::Tr1, y=2) = x+y
@test_throws ErrorException @traitfn defargs3(x::::(!Tr1), y) = x-y

@traitfn defargs4(x::::Tr1, y=2; k=1) = x+y+k
@traitfn defargs4(x::::(!Tr1), y=2; k=2) = x-y+k
@test defargs4(1,3)==5
@test defargs4(1)==4
@test defargs4(1.0,4)==-1
@test defargs4(1.0)==1
@test defargs4(1.0, k=10)==9

@traitfn defargs5(x::::Tr1, y=2, z...; k=1) = (x+y+k,z)
@traitfn defargs5(x::::(!Tr1), y=2, z...; k=2) = (x-y+k,z)
@test defargs5(1,3)==(5,())
@test defargs5(1,3,4,5)==(5,(4,5))
@test defargs5(1)==(4,())
@test defargs5(1.0,4)==(-1,())
@test defargs5(1.0,k=10)==(9,())

@traitfn defargs6(x::X=1, y=2) where {X;  Tr1{X}} = x+y
@traitfn defargs6(x::X=1, y=2) where {X; !Tr1{X}} = x-y
@test defargs6()==3
@test defargs6(1,3)==4
@test defargs6(1)==3
@test defargs6(1.0,4)==-3
@test defargs6(1.0)==-1
# above does not work with Traitor syntax
@test_broken SimpleTraits.traitfn(:(defargs6a(x::::Tr1=1, y=2) = x+y))

# traitfn with macro
@traitfn @inbounds gg(x::X) where {X; Tr1{X}} = x
@test gg(5)==5
@traitfn @generated ggg(x::X) where {X; Tr1{X}} = X<:AbstractArray ? :(x+x) : :(x)
@test ggg(5)==5
@traitimpl Tr1{AbstractArray}
@test ggg([5])==[10]

@traitfn @inbounds ggt(x::::Tr1) = x
@test ggt(5)==5
@traitfn @generated gggt(x::X::Tr1) where {X} = X<:AbstractArray ? :(x+x) : :(x)
@test gggt(5)==5
@test gggt([5])==[10]

# traitfn with Type
@traitfn ggt(::Type{X}, y) where {X; Tr1{X}} = (X,y)
@test ggt(Array, 5)==(Array, 5)
# no equivalent with Traitor syntax

# traitfn with ::X
@traitfn gg27(::X) where {X; Tr1{X}} = X
@test gg27([1])==Array{Int,1}

@traitfn gg27t(::X::Tr1) where {X} = X
@test gg27t([1])==Array{Int,1}

##
@traitfn f11(x::Dict{T}) where {T<:Number;  Tr1{Dict{T}}} = 1
@traitfn f11(x::Dict{T}) where {T<:Number; !Tr1{Dict{T}}} = 2
@traitimpl Tr1{Dict{Int}}
@test f11(Dict(1=>1))==1
@test f11(Dict(5.5=>1))==2

@traitfn f11t(x::Dict{T}::Tr1) where {T<:Number} = 1
@traitfn f11t(x::Dict{T}::(!Tr1)) where {T<:Number} = 2
@test f11t(Dict(1=>1))==1
@test f11t(Dict(5.5=>1))==2

##
@traitfn f12t(::::Tr1) = 1
@traitfn f12t(::::(!Tr1)) = 2
@test f12t(1)==1
@test f12t(5.5)==2

###
# @traitimpl Tr{X} <- istr(X)
@traitdef TrArrow1{X}
isarrow(X) = eltype(X)<:Integer ? true : false
@traitimpl TrArrow1{X} <- isarrow(X)
@test istrait(TrArrow1{Vector{Int}})
@test !istrait(TrArrow1{Vector{Float64}})

@traitdef TrArrow2{X}
@traitimpl Not{TrArrow2{X}} <- isarrow(X)
@test !istrait(TrArrow2{Vector{Int}})
@test istrait(TrArrow2{Vector{Float64}})

####
# issue with key in dispatch_cache
@traitfn f_dc(::::Tr1) = 1
module Mod
using SimpleTraits
using Compat
using Test
@traitdef Tr1{X}
@traitimpl Tr1{Integer}
@traitfn f_dc(::::Tr1) = 2
@test f_dc(1) == 2
end

######
# Other tests
#####
include("base-traits.jl")
include("base-traits-inference.jl")
include("backtraces.jl")


text:
# This example shows a simple use for a remote channel and a feeder

using Distributed
@everywhere using Logging
@everywhere using Base.CoreLogging

# Set Debug for Logging
@everywhere global_logger(ConsoleLogger(stdout, Logging.Debug))

function read_and_print(channel::RemoteChannel)

    @debug("Running updater.")

    while true

        x, i = try
            take!(channel);
        catch e

            if isa(e, InvalidStateException)
                break
            end

            @warn("Something wrong when reading from channel. Will skip.", e)

            continue

        end

        @debug("Updater has read values from channel. x = $(x) and i = $(i).")

    end

    @debug("Update channel closed. Exiting thread.")

end

c = RemoteChannel(() -> Channel{Tuple{Vector{Int}, Int}}(0))

# Task in main process that waits for the channel to be fed
task1 = @async read_and_print(c)

# Feeding the channel
task2 = @sync @distributed for i = 1:10
    put!(c, (ones(Int, 3), i))
    sleep(rand() / 2)
    @debug("Process $(myid()) has fed the channel.")
end

sleep(1)

# Wrong feeding
put!(c, 1)

# This should close `task1`
close(c)

# Just to be sure
sleep(1)

println("Task1 done = $(istaskdone(task1))")

println("Task2 done = $(istaskdone(task2))")


text:
# NegativeBinomial is the distribution of the number of failures
# before the r-th success in a sequence of Bernoulli trials.
# We do not enforce integer size, as the distribution is well defined
# for non-integers, and this can be useful for e.g. overdispersed
# discrete survival times.

immutable NegativeBinomial <: DiscreteUnivariateDistribution
    r::Float64
    prob::Float64

    function NegativeBinomial(r::Real, p::Real)
        zero(p) < p <= one(p) || error("prob must be in (0, 1].")
        zero(r) < r || error("r must be positive.")
        new(float64(r), float64(p))
    end

    NegativeBinomial() = new(1.0, 0.5)
end

@_jl_dist_2p NegativeBinomial nbinom

isupperbounded(::Union(NegativeBinomial, Type{NegativeBinomial})) = false
islowerbounded(::Union(NegativeBinomial, Type{NegativeBinomial})) = true
isbounded(::Union(NegativeBinomial, Type{NegativeBinomial})) = false

minimum(::Union(NegativeBinomial, Type{NegativeBinomial})) = 0
maximum(::Union(NegativeBinomial, Type{NegativeBinomial})) = Inf

insupport(::NegativeBinomial, x::Real) = isinteger(x) && zero(x) <= x
insupport(::Type{NegativeBinomial}, x::Real) = isinteger(x) && zero(x) <= x

function probs(d::NegativeBinomial, rgn::UnitRange)
    r = d.r
    p0 = 1.0 - d.prob
    f, l = rgn[1], rgn[end]
    0 <= f <= l || throw(BoundsError())
    res = Array(Float64, l - f + 1)
    res[1] = v = pdf(d, f)
    b = f - 1
    for x = f+1:l
        c = (x + r - 1) * p0 / x
        res[x-b] = (v *= c)
    end
    return res
end

function mgf(d::NegativeBinomial, t::Real)
    r, p = d.r, d.prob
    return ((1.0 - p) * exp(t))^r / (1.0 - p * exp(t))^r
end

function cf(d::NegativeBinomial, t::Real)
    r, p = d.r, d.prob
    return ((1.0 - p) * exp(im * t))^r / (1.0 - p * exp(im * t))^r
end

function mean(d::NegativeBinomial)
    p = d.prob
    (1.0 - p) * d.r / p
end

function var(d::NegativeBinomial)
    p = d.prob
    (1.0 - p) * d.r / (p * p)
end

function std(d::NegativeBinomial)
    p = d.prob
    sqrt((1.0 - p) * d.r) / p
end

function skewness(d::NegativeBinomial)
    p = d.prob
    (2.0 - p) / sqrt((1.0 - p) * d.r)
end

function kurtosis(d::NegativeBinomial)
    p = d.prob
    6.0 / d.r + (p * p) / ((1.0 - p) * d.r)
end

function mode(d::NegativeBinomial)
    p = d.prob
    ifloor((1.0 - p) * (d.r - 1.) / p)
end