text:
#######################
### Document Struct ###
#######################

mutable struct Document
	terms::Vector{Int}
	counts::Vector{Int}
	readers::Vector{Int}
	ratings::Vector{Int}
	stamp::Float64
	title::String

	function Document(terms; counts=ones(length(terms)), readers=Int[], ratings=ones(length(readers)), stamp=-Inf, title="")
		doc = new(terms, counts, readers, ratings, stamp[1], title)
		checkdoc(doc)
		return doc
	end
end

#####################
### Corpus Struct ###
#####################

mutable struct Corpus
	docs::Vector{Document}
	lex::Dict{Int, String}
	users::Dict{Int, String}

	function Corpus(;docs=Document[], lex=[], users=[])
		isa(lex, Dict) || (lex = Dict(lkey => term for (lkey, term) in enumerate(lex)))
		isa(users, Dict) || (users = Dict(ukey => user for (ukey, user) in enumerate(users)))

		corp = new(docs, lex, users)
		checkcorp(corp)
		return corp
	end
end

#################################################
### Functions for Reading and Writing Corpora ###
#################################################

function readcorp(;docfile::AbstractString="", lexfile::AbstractString="", userfile::AbstractString="", titlefile::AbstractString="", delim::Char=',', counts::Bool=false, readers::Bool=false, ratings::Bool=false, stamps::Bool=false)
	(ratings <= readers) || (ratings = false; warn("Ratings require readers, ratings switch set to false."))
	(!isempty(docfile) | isempty(titlefile)) || warn("No docfile, titles will not be assigned.")
	stamp = stamps

	corp = Corpus()
	if !isempty(docfile)
		docs = open(docfile)
		dockwargs = [:counts, :readers, :ratings, :stamp]
		for (d, docblock) in enumerate(partition(readlines(docs), counts + readers + ratings + stamp + 1))
			try
			doclines = Vector{Float64}[[parse(Float64, p) for p in split(line, delim)] for line in docblock]
			docinput = zip(dockwargs[[counts, readers, ratings, stamp]], doclines[2:end])
			push!(corp, Document(doclines[1]; docinput...))
			catch error("Document $d beginning on line $((d - 1) * (counts + readers + ratings + stamp) + d) failed to load.")
			end
		end
	else warn("No docfile, topic models cannot be trained without documents.")
	end

	if !isempty(lexfile)
		lex = readdlm(lexfile, '\t', comments=false)
		lkeys = lex[:,1]
		terms = [string(j) for j in lex[:,2]]
		corp.lex = Dict{Int, String}(zip(lkeys, terms))
		@assert all(ispositive.(collect(keys(corp.lex))))
	end

	if !isempty(userfile)
		users = readdlm(userfile, '\t', comments=false)
		ukeys = users[:,1]
		users = [string(u) for u in users[:,2]]
		corp.users = Dict{Int, String}(zip(ukeys, users))
		@assert all(ispositive.(collect(keys(corp.users))))
	end

	if !isempty(titlefile)
		titles = readdlm(titlefile, '\n', String)
		for (d, doc) in enumerate(corp)
			doc.title = titles[d]
		end
	end

	return corp
end

function writecorp(corp::Corpus; docfile::AbstractString="", lexfile::AbstractString="", userfile::AbstractString="", titlefile::AbstractString="", delim::Char=',', counts::Bool=false, readers::Bool=false, ratings::Bool=false, stamps::Bool=false)
	(ratings <= readers) || (ratings = false; warn("Ratings require readers, ratings switch set to false."))
	stamp = stamps

	if !isempty(docfile)
		dockwargs = [:counts, :readers, :ratings, :stamp]
		dockwargs = dockwargs[[counts, readers, ratings, stamp]]
		docfile = open(docfile, "w")
		for doc in corp
			write(docfile, join(doc.terms, delim), '\n')
			for arg in dockwargs
				write(docfile, join(doc.(arg), delim), '\n')
			end
		end
		close(docfile)
	end

	if !isempty(lexfile)
		lkeys = sort(collect(keys(corp.lex)))
		lex = zip(lkeys, [corp.lex[lkey] for lkey in lkeys])
		writedlm(lexfile, lex)
	end

	if !isempty(userfile)
		ukeys = sort(collect(keys(corp.users)))
		users = zip(ukeys, [corp.users[ukey] for ukey in ukeys])
		writedlm(userfile, users)
	end

	if !isempty(titlefile)
		writedlm(titlefile, [doc.title for doc in corp])
	end
	nothing
end

#######################
### corp! Functions ###
#######################

function abridgecorp!(corp::Corpus; stop::Bool=false, order::Bool=true, abr::Integer=1)
	v = "v$(VERSION.major).$(VERSION.minor)"

	if stop
		stopwords = vec(readdlm(pwd() * "/.julia/$v/topicmodelsvb/datasets/stopwords.txt", String))
		stopkeys = filter(j -> lowercase(corp.lex[j]) in stopwords, collect(keys(corp.lex)))
		for doc in corp
			keep = Bool[!(j in stopkeys) for j in doc.terms]
			doc.terms = doc.terms[keep]
			doc.counts = doc.counts[keep]
		end
	end

	if !order
		for doc in corp
			docdict = Dict(Int(j) => 0 for j in doc.terms)
			for (j, c) in zip(doc.terms, doc.counts)
				docdict[j] += c
			end
			doc.terms = collect(keys(docdict))
			doc.counts = collect(values(docdict))
		end
	end

	if abr > 1
		doclkeys = Set(vcat([doc.terms for doc in corp]...))
		lexcount = Dict(Int(j) => 0 for j in doclkeys)
		for doc in corp, (j, c) in zip(doc.terms, doc.counts)
			lexcount[j] += c
		end

		for doc in corp
			keep = Bool[lexcount[j] >= abr for j in doc.terms]
			doc.terms = doc.terms[keep]
			doc.counts = doc.counts[keep]
		end

		keep = Bool[length(doc.terms) > 0 for doc in corp]
		corp.docs = corp[keep]
	end
	nothing
end

function trimcorp!(corp::Corpus; lex::Bool=true, terms::Bool=true, users::Bool=true, readers::Bool=true)
	if lex
		doclkeys = Set(vcat([doc.terms for doc in corp]...))
		corp.lex = Dict(lkey => corp.lex[lkey] for lkey in intersect(keys(corp.lex), doclkeys))
	end

	if terms
		doclkeys = Set(vcat([doc.terms for doc in corp]...))
		boguslkeys = setdiff(doclkeys, keys(corp.lex))
		for doc in corp
			keep = Bool[!(j in boguslkeys) for j in doc.terms]
			doc.terms = doc.terms[keep]
			doc.counts = doc.counts[keep]
		end
	end

	if users
		docukeys = Set(vcat([doc.readers for doc in corp]...))
		corp.users = Dict(ukey => corp.users[ukey] for ukey in intersect(keys(corp.users), docukeys))
	end

	if readers
		docukeys = Set(vcat([doc.readers for doc in corp]...))
		bogusukeys = setdiff(docukeys, keys(corp.users))
		for doc in corp
			keep = Bool[!(u in bogusukeys) for u in doc.readers]
			doc.readers = doc.readers[keep]
			doc.ratings = doc.ratings[keep]
		end
	end
	nothing
end

function compactcorp!(corp::Corpus; lex::Bool=true, users::Bool=true, alphabetize::Bool=true)
	if lex
		lkeys = sort(collect(keys(corp.lex)))
		lkeymap = zip(lkeys, 1:length(corp.lex))
		if alphabetize
			alphabetdict = Dict(j => lkey for (j, lkey) in zip(sortperm([corp.lex[lkey] for lkey in lkeys]), 1:length(corp.lex)))
			lkeydict = Dict(lkey => alphabetdict[j] for (lkey, j) in lkeymap)
		else
			lkeydict = Dict(lkey => j for (lkey, j) in lkeymap)
		end

		corp.lex = Dict(lkeydict[lkey] => corp.lex[lkey] for lkey in keys(corp.lex))
		for lkey in keys(corp.lex)
			try if corp.lex[lkey][1:5] == "#term"; corp.lex[lkey] = "#term$(lkey)"; end
			end
		end
	end

	if users
		ukeys = sort(collect(keys(corp.users)))
		ukeymap = zip(ukeys, 1:length(corp.users))
		if alphabetize
			alphabetdict = Dict(r => ukey for (r, ukey) in zip(sortperm([corp.users[ukey] for ukey in ukeys]), 1:length(corp.users)))
			ukeydict = Dict(ukey => alphabetdict[r] for (ukey, r) in ukeymap)
		else
			ukeydict = Dict(ukey => r for (ukey, r) in ukeymap)
		end

		corp.users = Dict(ukeydict[ukey] => corp.users[ukey] for ukey in keys(corp.users))
		for ukey in keys(corp.users)
			try if corp.users[ukey][1:5] == "#user"; corp.users[ukey] = "#user$(ukey)"; end
			end
		end
	end

	for doc in unique(corp)
		if lex; doc.terms = [lkeydict[lkey] for lkey in doc.terms]; end
		if users; doc.readers = [ukeydict[ukey] for ukey in doc.readers]; end
	end
	nothing
end

function padcorp!(corp::Corpus; lex::Bool=true, users::Bool=true)
	if lex
		doclkeys = Set(vcat([doc.terms for doc in corp]...))
		for lkey in setdiff(doclkeys, keys(corp.lex))
			corp.lex[lkey] = string(join(["#term",lkey]))
		end
	end
	if users
		docukeys = Set(vcat([doc.readers for doc in corp]...))
		for ukey in setdiff(docukeys, keys(corp.users))
			corp.users[ukey] = string(join(["#user",ukey]))
		end
	end
	nothing
end

function cullcorp!(corp::Corpus; lex::Bool=false, users::Bool=false, len::Integer=1)
	lexkeys = keys(corp.lex)
	userkeys = keys(corp.users)
	bogusdocs = Int[]
	for (d, doc) in enumerate(corp)
		if lex
			for j in doc.terms
				if !(j in lexkeys)
					push!(bogusdocs, d)
					break
				end
			end
		end

		if users & !(d in bogusdocs)
			for u in doc.readers
				if !(u in userkeys)
					push!(bogusdocs, d)
					break
				end
			end
		end

		if (len > 1) & !(d in bogusdocs)
			if length(doc) < len
				push!(bogusdocs, d)
			end
		end

		if isempty(doc.terms) & !(d in bogusdocs)
			push!(bogusdocs, d)
		end
	end

	deleteat!(corp, bogusdocs)
	nothing
end

function fixcorp!(corp::Corpus; lex::Bool=true, terms::Bool=true, users::Bool=true, readers::Bool=true, stop::Bool=false, order::Bool=true, abr::Integer=1, len::Integer=1, alphabetize::Bool=true)
	println("Abridging corpus..."); abridgecorp!(corp, stop=stop, order=order, abr=abr)
	println("Trimming corpus..."); trimcorp!(corp, lex=lex, terms=terms, users=users, readers=readers)
	println("Culling corpus..."); cullcorp!(corp, len=len)
	println("Compacting corpus..."); compactcorp!(corp, lex=lex, users=users, alphabetize=alphabetize)
	nothing
end

########################
### Corpus Utilities ###
########################

function Base.show(io::IO, doc::Document)
	print(io, "Document with:\n * $(length(doc.terms)) terms\n * $(length(doc.readers)) readers")
	if isfinite(doc.stamp)
		print(io, "\n * $(doc.stamp) stamp")
	end
end

Base.length(doc::Document) = length(doc.terms)
Base.size(doc::Document) = sum(doc.counts)
Base.in(doc::Document, corp::Corpus) = in(doc, corp.docs)

Base.show(io::IO, corp::Corpus) = print(io, "Corpus with:\n * $(length(corp)) docs\n * $(length(corp.lex)) lex\n * $(length(corp.users)) users")
Base.start(corp::Corpus) = 1
Base.next(corp::Corpus, d::Int) = corp.docs[d], d + 1
Base.done(corp::Corpus, d::Int) = length(corp.docs) == d - 1
Base.push!(corp::Corpus, doc::Document) = push!(corp.docs, doc)
Base.pop!(corp::Corpus) = pop!(corp.docs)
Base.unshift!(corp::Corpus, doc::Document) = unshift!(corp.docs, doc)
Base.unshift!(corp::Corpus, docs::Vector{Document}) = unshift!(corp.docs, docs)
Base.shift!(corp::Corpus) = shift!(corp.docs)
Base.insert!(corp::Corpus, d::Int, doc::Document) = insert!(corp.docs, d, doc)
Base.deleteat!(corp::Corpus, d::Int) = deleteat!(corp.docs, d)
Base.deleteat!(corp::Corpus, ds::Vector{Int}) = deleteat!(corp.docs, ds)
Base.deleteat!(corp::Corpus, ds::UnitRange{Int}) = deleteat!(corp.docs, ds)
Base.getindex(corp::Corpus, d::Int) = getindex(corp.docs, d)
Base.getindex(corp::Corpus, ds::Vector{Int}) = getindex(corp.docs, ds)
Base.getindex(corp::Corpus, ds::UnitRange{Int}) = getindex(corp.docs, ds)
Base.getindex(corp::Corpus, ds::Vector{Bool}) = getindex(corp.docs, find(ds))
Base.setindex!(corp::Corpus, doc::Document, d::Int) = setindex!(corp.docs, doc, d)
Base.setindex!(corp::Corpus, docs::Vector{Document}, ds::Vector{Int}) = setindex!(corp.docs, docs, ds)
Base.setindex!(corp::Corpus, docs::Vector{Document}, ds::UnitRange{Int}) = setindex!(corp.docs, docs, ds)
Base.findfirst(corp::Corpus, doc::Document) = findfirst(corp.docs, doc)
Base.findin(corp::Corpus, docs::Vector{Document}) = findin(corp.docs, docs)
Base.findin(corp::Corpus, doc::Document) = findin(corp, [doc])
Base.length(corp::Corpus) = length(corp.docs)
Base.size(corp::Corpus) = (length(corp), length(corp.lex), length(corp.users))
Base.copy(corp::Corpus) = Corpus(docs=copy(corp.docs), lex=copy(corp.lex), users=copy(corp.users))
Base.endof(corp::Corpus) = length(corp)

function checkdoc(doc::Document)
	pass =
	(!isempty(doc.terms)
	& all(ispositive.(doc.terms))
	& all(ispositive.(doc.counts))
	& isequal(length(doc.terms), length(doc.counts))
	& all(ispositive.(doc.readers))
	& all(ispositive.(doc.ratings))
	& isequal(length(doc.readers), length(doc.ratings)))
	return pass
end

function checkcorp(corp::Corpus)
	pass = true
	for (d, doc) in enumerate(corp)
		checkdoc(doc) || (println("Document $d failed check."); pass=false)
	end
	@assert all(ispositive.(collect(keys(corp.lex))))
	@assert all(ispositive.(collect(keys(corp.users))))
	return pass
end

function showdocs{T<:Integer}(corp::Corpus, ds::Vector{T})
	@assert checkbounds(Bool, 1:length(corp), ds) "Some document indices outside docs range."

	for d in ds
		doc = corp[d]
		@juliadots "doc $d\n"
		if !isempty(doc.title)
			@juliadots "$(doc.title)\n"
		end
		println(join([corp.lex[lkey] for lkey in corp[d].terms], " "), '\n')
	end
end

function showdocs(corp::Corpus, docs::Vector{Document})
	ds = findin(corp, docs)
	showdocs(corp, ds)
end

showdocs{T<:Integer}(corp::Corpus, ds::UnitRange{T}) = showdocs(corp, collect(ds))
showdocs(corp::Corpus, d::Integer) = showdocs(corp, [d])
showdocs(corp::Corpus, doc::Document) = showdocs(corp, [doc])

getlex(corp::Corpus) = sort(collect(values(corp.lex)))
getusers(corp::Corpus) = sort(collect(values(corp.users)))

######################################
### Pre-packaged Dataset Shortcuts ###
######################################

function readcorp(corpsym::Symbol)
	v = "v$(VERSION.major).$(VERSION.minor)"

	if corpsym == :nsf
		docfile = homedir() * "/.julia/$v/topicmodelsvb/datasets/nsf/nsfdocs.txt"
		lexfile = homedir() * "/.julia/$v/topicmodelsvb/datasets/nsf/nsflex.txt"
		titlefile = homedir() * "/.julia/$v/topicmodelsvb/datasets/nsf/nsftitles.txt"
		corp = readcorp(docfile=docfile, lexfile=lexfile, titlefile=titlefile, counts=true, stamps=true)

	elseif corpsym == :citeu
		docfile = homedir() * "/.julia/$v/topicmodelsvb/datasets/citeu/citeudocs.txt"
		lexfile = homedir() * "/.julia/$v/topicmodelsvb/datasets/citeu/citeulex.txt"
		userfile = homedir() * "/.julia/$v/topicmodelsvb/datasets/citeu/citeuusers.txt"
		titlefile = homedir() * "/.julia/$v/topicmodelsvb/datasets/citeu/citeutitles.txt"
		corp = readcorp(docfile=docfile, lexfile=lexfile, userfile=userfile, titlefile=titlefile, counts=true, readers=true)
		padcorp!(corp)

	elseif corpsym == :mac
		docfile = homedir() * "/.julia/$v/topicmodelsvb/datasets/mac/macdocs.txt"
		lexfile = homedir() * "/.julia/$v/topicmodelsvb/datasets/mac/maclex.txt"
		titlefile = homedir() * "/.julia/$v/topicmodelsvb/datasets/mac/mactitles.txt"
		corp = readcorp(docfile=docfile, lexfile=lexfile, titlefile=titlefile, counts=true, stamps=true)

	else
		println("Included corpora:\n:nsf\n:citeu\n:mac")
		corp = nothing
	end

	return corp
end


text:
module PDSampler

using ApproxFun
using StatsBase: autocov
using Polynomials:
        Poly,
        roots,
        polyint,
        polyval
using DataStructures:
        PriorityQueue,
        enqueue!,
        dequeue!,
        peek
using ProgressMeter
using Distributions: Beta
using LinearAlgebra
using Random
using Statistics

export
    # MODELS
    loglik,
    gradloglik,
    gradloglik_cv

### source files (keep the order)

include("models/mvgaussian.jl")
include("models/tmvgaussian.jl")
include("models/logreg.jl")
include("models/pmf.jl")

include("geometry.jl")
include("ippsampler.jl")
include("path.jl")
include("kernels.jl")
include("simulate.jl")

include("local/event.jl")
include("local/factorgraph.jl")
include("local/simulate.jl")

end # module

ext: jl


text:
module Batchifier

import ..MCTS, ..Util

export BatchedOracle

"""
    launch_server(f, num_workers)

Launch an inference requests server.

# Arguments

  - `f` is the function to be evaluated (e.g. a neural network).
    It takes a batch of inputs and returns a batch of results
    of similar size.
  - `num_workers` is the number of workers that are expected to
    query the server. It also corresponds to the batch size that is used to
    evaluate `f`.

# How to use

This function returns a channel along which workers can send queries.
A query can be either:

    - `:none` when a worker is done and about to terminate
    - a named tuple with fields `query` (the input to be given to `f`) and
      `answer_channel` (the channel the sever must use to return its answer).

The server stops automatically after all workers send `:none`.
"""
function launch_server(f, num_workers)
  channel = Channel(num_workers)
  Threads.@spawn Util.@printing_errors begin
    num_active = num_workers
    pending = []
    while num_active > 0
      req = take!(channel)
      if req == :done
        num_active -= 1
      else
        push!(pending, req)
      end
      @assert length(pending) <= num_active
      if length(pending) == num_active && num_active > 0
        batch = [p.query for p in pending]
        results = f(batch)
        for i in eachindex(pending)
          put!(pending[i].answer_channel, results[i])
        end
        empty!(pending)
      end
    end
  end
  return channel
end

client_done!(reqc) = put!(reqc, :done)

struct BatchedOracle{Game, F} <: MCTS.Oracle{Game}
  make_query :: F # turn state into a query (this is usually the identity)
  reqchan :: Channel
  anschan :: Channel
  function BatchedOracle{G}(f, reqchan) where G
    return new{G, typeof(f)}(f, reqchan, Channel(1))
  end

end

function BatchedOracle{G}(reqchan) where G
  return BatchedOracle{G}(x -> x, reqchan)
end

function (oracle::BatchedOracle)(state)
  query = oracle.make_query(state)
  put!(oracle.reqchan, (query=query, answer_channel=oracle.anschan))
  answer = take!(oracle.anschan)
  return answer
end

end


text:
const CR_BYTE = 0x0d # \r
const LF_BYTE = 0x0a # \n
const DASH_BYTE = 0x2d # -
const HTAB_BYTE = 0x09 # \t
const SPACE_BYTE = 0x20
const RETURN_BYTES = [CR_BYTE, LF_BYTE]

"""
    find_multipart_boundary(bytes, boundaryDelimiter; start::Int=1)

Find the first and last index of the next boundary delimiting a part, and if
the discovered boundary is the terminating boundary.
"""
function find_multipart_boundary(bytes::AbstractVector{UInt8}, boundaryDelimiter::AbstractVector{UInt8}; start::Int=1)
    # The boundary delimiter line is prepended with two '-' characters
    # The boundary delimiter line starts on a new line, so must be preceded by a \r\n.
    # The boundary delimiter line ends with \r\n, and can have "optional linear whitespace" between
    # the end of the boundary delimiter, and the \r\n.
    # The last boundary delimiter line has an additional '--' at the end of the boundary delimiter
    # [RFC2046 5.1.1](https://tools.ietf.org/html/rfc2046#section-5.1.1)

    i = start
    end_index = i+length(boundaryDelimiter)-1
    while end_index <= length(bytes)
        if bytes[i:end_index] == boundaryDelimiter
            # boundary delimiter line start on a new line ...
            if i > 1
                (i == 2 || bytes[i-2] != CR_BYTE || bytes[i-1] != LF_BYTE) && error("boundary delimiter found, but it was not the start of a line")
                # the CRLF preceding the boundary delimiter is "conceptually attached
                # to the boundary", so account for this with the index
                i-=2
            end

            # need to check if there are enough characters for the CRLF or for two dashes
            end_index < length(bytes)-1 || error("boundary delimiter found, but did not end with new line")

            is_terminating_delimiter = bytes[end_index+1] == DASH_BYTE && bytes[end_index+2] == DASH_BYTE
            is_terminating_delimiter && (end_index+=2)

            # ... there can be arbitrary SP and HTAB space between the boundary delimiter ...
            while (end_index < length(bytes) && (bytes[end_index+1] in [HTAB_BYTE, SPACE_BYTE]))
                end_index+=1
            end
            # ... and ends with a new line
            newlineEnd = end_index < length(bytes)-1 &&
                         bytes[end_index+1] == CR_BYTE &&
                         bytes[end_index+2] == LF_BYTE
            if !newlineEnd
                error("boundary delimiter found, but did not end with new line")
            end

            end_index += 2

            return (is_terminating_delimiter, i, end_index)
        end

        i += 1
        end_index += 1
    end

    error("boundary delimiter not found")
end

"""
    find_multipart_boundaries(bytes, boundary; start=1)

Find the start and end indexes of all the parts of the multipart object.  Ultimately this method is
looking for the data between the boundary delimiters in the byte array.  A vector containing all
the start/end pairs is returned.
"""
function find_multipart_boundaries(bytes::AbstractVector{UInt8}, boundary::AbstractVector{UInt8}; start=1)
    idxs = Tuple{Int, Int}[]
    while true
        (is_terminating_delimiter, i, end_index) = find_multipart_boundary(bytes, boundary; start = start)
        push!(idxs, (i, end_index))
        is_terminating_delimiter && break
        start = end_index + 1
    end
    return idxs
end

"""
    find_header_boundary(bytes)

Find the end of the multipart header in the byte array.  Returns a Tuple with the
start index(1) and the end index. Headers are separated from the body by CRLFCRLF.

[RFC2046 5.1](https://tools.ietf.org/html/rfc2046#section-5.1)
[RFC822 3.1](https://tools.ietf.org/html/rfc822#section-3.1)
"""
function find_header_boundary(bytes::AbstractVector{UInt8})
    delimiter = UInt8[CR_BYTE, LF_BYTE, CR_BYTE, LF_BYTE]
    length(delimiter) > length(bytes) && (return nothing)

    l = length(bytes) - length(delimiter) + 1
    i = 1
    end_index = length(delimiter)
    while (i <= l)
        bytes[i:end_index] == delimiter && (return (1, end_index))
        i += 1
        end_index += 1
    end
    error("no delimiter found separating header from multipart body")
end

"""
    content_disposition_tokenize(str)

Tokenize the "arguments" for the Content-Disposition declaration.  A vector of
strings is returned that contains each token and separator found in the source
string. Tokens are separated by either an equal sign(=) or a semi-colon(;) and
may be quoted or escaped with a backslash(\\). All tokens returned are stripped
of whitespace at the beginning and end of the string, quotes are retained.
"""
function content_disposition_tokenize(str)
    retval = Vector{SubString}()
    start = 1
    quotes = false
    escaped = false

    for offset in eachindex(str)
        if escaped == false
            if quotes == true && str[offset] == '"'
                quotes = false
            elseif str[offset] == '\\'
                escaped = true
            elseif str[offset] == '"'
                quotes = true
            elseif quotes == false && (str[offset] == ';' || str[offset] == '=')
                prev = prevind(str, offset)
                if prev > start
                    push!(retval, strip(SubString(str, start, prev)))
                end
                push!(retval, SubString(str, offset, offset))
                start = nextind(str, offset)
            end
        else
            escaped = false
        end
    end

    if start != lastindex(str)
        push!(retval, strip(SubString(str, start)))
    end

    return retval
end

"""
    content_disposition_extract(str)

Extract all the flags and key/value arguments from the Content-Disposition
line.  The result is returned as an array of tuples.

In the case of a flag the first value of the tuple is false the second value
is the flag and the third value is nothing.

In the case of a key/value argument the first value is true, the second is the
key, and the third is the value (or nothing if no value was specified).
"""
function content_disposition_extract(str)
    retval = Vector{Tuple{Bool, SubString, Union{SubString,Nothing}}}()
    tokens = content_disposition_tokenize(str)
    total  = length(tokens)

    function strip_quotes(val)
        if val[1] == '"' && val[end] == '"'
            SubString(val, 2, lastindex(val) - 1)
        else
            val
        end
    end

    i = 1
    while i < total
        if tokens[i] != ';'
            pair  = (i + 1 <= total && tokens[i + 1] == "=")
            key   = strip_quotes(tokens[i])
            value = (pair && i + 2 <= total && tokens[i + 2] != ";" ? strip_quotes(tokens[i + 2]) : nothing)

            push!(retval, (pair, key, value))

            if pair
                i += 3
            else
                i += 1
            end
        else
            i += 1
        end
    end
    return retval
end

"""
    parse_multipart_chunk(chunk)

Parse a single multi-part chunk into a Multipart object.  This will decode
the header and extract the contents from the byte array.
"""
function parse_multipart_chunk(chunk)
    (startIndex, end_index) = find_header_boundary(chunk)

    headers = String(view(chunk, startIndex:end_index))
    content = view(chunk, end_index+1:lastindex(chunk))

    disposition = match(r"(?i)Content-Disposition: form-data(.*)\r\n", headers)

    if disposition === nothing
        @warn "Content disposition is not specified dropping the chunk." chunk
        return # Specifying content disposition is mandatory
    end

    name = nothing
    filename = nothing

    for (pair, key, value) in content_disposition_extract(disposition[1])
        if pair && key == "name"
            name = value
        elseif pair && key == "filename"
            filename = value
        end
    end

    name === nothing && return

    match_contenttype = match(r"(?i)Content-Type: (\S*[^;\s])", headers)
    contenttype = match_contenttype !== nothing ? match_contenttype[1] : "text/plain" # if content_type is not specified, the default text/plain is assumed

    return Multipart(filename, IOBuffer(content), contenttype, "", name)
end

"""
    parse_multipart_body(body, boundary)::Vector{Multipart}

Parse the multipart body received from the client breaking it into the various
chunks which are returned as an array of Multipart objects.
"""
function parse_multipart_body(body::AbstractVector{UInt8}, boundary::AbstractString)::Vector{Multipart}
    multiparts = Multipart[]
    idxs = find_multipart_boundaries(body, Vector{UInt8}("--$(boundary)"))
    length(idxs) > 1 || (return multiparts)

    for i in 1:length(idxs)-1
        chunk = view(body, idxs[i][2]+1:idxs[i+1][1]-1)
        push!(multiparts, parse_multipart_chunk(chunk))
    end
    return multiparts
end

"""
    parse_multipart_form(req::Request)::Vector{Multipart}

Parse the full mutipart form submission from the client returning and
array of Multipart objects containing all the data.

The order of the multipart form data in the request should be preserved.
[RFC7578 5.2](https://tools.ietf.org/html/rfc7578#section-5.2).

The boundary delimiter MUST NOT appear inside any of the encapsulated parts. Note
that the boundary delimiter does not need to have '-' characters, but a line using
the boundary delimiter will start with '--' and end in \r\n.
[RFC2046 5.1](https://tools.ietf.org/html/rfc2046#section-5.1.1)
"""
function parse_multipart_form(req::Request)::Union{Vector{Multipart}, Nothing}
    # parse boundary from Content-Type
    m = match(r"multipart/form-data; boundary=(.*)$", req["Content-Type"])
    m === nothing && return nothing

    boundary_delimiter = m[1]

    # [RFC2046 5.1.1](https://tools.ietf.org/html/rfc2046#section-5.1.1)
    length(boundary_delimiter) > 70 && error("boundary delimiter must not be greater than 70 characters")

    return parse_multipart_body(payload(req), boundary_delimiter)
end


text:
# DBL_EPSILON - jest to dokładność dla liczb zmiennoprzecinkowych
# dla Float64 jest to 2.2204460492503131e-16
DBL_EPSILON = 2.2204460492503131e-16

function DIFFSIGN(x, y)
    if (x <=0 && y >= 0) || (x >= 0 && y <= 0)
        return true
    else
        return false
    end
end

function fun(x)
    return 1 / (x - 3) - 6
end

function between(x, a, b)
    if b > a
        return (x >= a && x <= b)
    else
        return (x >= b && x <= a)
    end
end

function lfun(b, a, fb, fa)
    if fb != fa
        return b - fb*(b - a) / (fb - fa) # metoda siecznych
    elseif fa != 0
        return Inf
    else
        return b
    end
end

function hfun(b, c)
    if c > b
        return b + abs(b*DBL_EPSILON)
    else
        return b - abs(b*DBL_EPSILON)
    end
end

function mfun(b, c)
    return 0.5*(b + c) # metoda bisekcji
end

function wfun(l, b, c)
    h = hfun(b, c)
    m = mfun(b, c)

    if between(l, h, m) == true
         return l
    elseif (fun(abs(l - b)) <= abs(b*DBL_EPSILON)) && (between(l, b, m) == false)
        return h
    else
        return m
    end
end

function ffun(a, b, fa, fb)
     return (fa - fb) / (a - b)
end

function rfun(b, a, d, fb, fa, fd)
    alpha = ffun(b, d, fb, fd)*fa
    beta = ffun(a, d, fa, fd)*fb

    if beta != alpha
        return b - beta*(b - a) / (beta - alpha)
    elseif alpha != 0
        return Inf
    else
        return 0 # beta == alpha == 0
    end
end

function DekkerM(x0, x1, Eps)
    """
    b - ostatnie przyblizenie wyniku
    c - kontrapunkt b, punkt w którym funkcja f ma przeciwny znak niż w punkcie b
    a – poprzednia wartość punktu a, używana do wyliczania następnego punktu metodą siecznych

    Metoda bisekcji z punktów b i c tworzy punkt m pomiędzy nimi na środku przedziału.

    Wyliczany jest ciąg xi, którego ostatni element oznaczany jest przez x, a poprzedni przez xp.

    xk - ostatni punkt w ciągu, który ma różny znak niż x.

    Punkt x wyliczany jest dwoma metodami:
        - siecznych
        - bisekcji
    i wybierany jest ten obliczony z metody siecznych jeśli leży pomiędzy punktem b
    (ze względów dokładnościowych z pewną poprawką) a punktem m wyliczonym z bisekcji.

    Jeżeli f(x) czy f(xk) leży bliżej zera i jeśli f(x) leży bliżej zera
    wtedy b ma wartość x, c ma wartość xk, w przeciwnym razie zamiana.
    """
    d = NaN
    fd = NaN

    fxp = fun(x0)
    fx = fun(x1)

    if x0 == x1
        return fx
    end

    if DIFFSIGN(fx, fxp) == false
        return 0
    end

    if abs(fx) <= abs(fxp)
        b = x1
        a = c = x0
        fa = fxp
        fb = fx
    else
        b = x0
        a = c = x1
        fa = fx
        fb = fxp
    end

    xk = xp = x0
    fxk = fxp
    x = x1
    iter = 0
    age = 0
    bp = b
    cp = c
    ap = a

    while abs(b - c) > 2*Eps
        iter = iter + 1
        age = age + 1

        if abs(b - c) <= (0.5 + 2 * DBL_EPSILON)*(abs(bp - cp) + abs(b*DBL_EPSILON))
            age = 1
        end

        xp = x

        if iter == 2
            lambda = lfun(b, a, fb, fa)

            if abs(lambda - b) < abs(b*DBL_EPSILON)
                break
            end
            x = wfun(lambda, b, c)
        elseif iter >= 3 && age <= 3
            rho = rfun(b, a, d, fb, fa, fd)

            if abs(rho - b) < abs(b*DBL_EPSILON)
                break
            end
            x = wfun(rho, b, c)
        elseif iter >=3 && age == 4
            rho = rfun(b, a, d, fb, fa, fd)

            if abs(rho - b) < abs(b*DBL_EPSILON)
                break
            end
            x = wfun(2 * rho - b, b, c)
        else
            x = mfun(b, c)
        end

        fxp = fx
        fx = fun(x)

        if DIFFSIGN(fxp, fx) == true
            xk = xp
            fxk = fxp
        end

        bp = b
        fbp = fb
        ap = a
        fap = fa
        cp = c

        if abs(fx) <= abs(fxk)
            a = b
            fa = fb
            b = x
            fb = fx
            c = xk
        else
            b = xk
            fb = fxk
            a = c = x
            fa = fx
        end

        if b == x || b == bp
            d = ap
            fd = fap
        else
            d = ap
            fd = fdp
        end
    end
    println("Number of iterations: ", iter)
    return b
end

println("f(x) = 1/(x-3)-6\n")
result = DekkerM(3.01, 4, 1e-12)
println("x0 = ", result)

max_stars_repo_path: Dekker Method R.jl

max_stars_repo_name: radoslawczapp/DekkerMethod

max_issues_repo_path: Dekker Method R.jl

max_issues_repo_name: radoslawczapp/DekkerMethod

max_forks_repo_path: Dekker Method R.jl

max_forks_repo_name: radoslawczapp/DekkerMethod


text:
@require SQLite begin
using SQLite
export table

immutable SQLiteTable
    db::SQLite.DB
    tablename::String
end

function table{T<:AbstractString}(db::SQLite.DB,tablename::T)
    return SQLiteTable(db,tablename)
end

immutable QueryableSQLite{T,Provider} <: Queryable{T,Provider}
    db::SQLite.DB
    tablename::AbstractString
end

type QueryProviderSQLite <: QueryProvider
end

function query(table::SQLiteTable)
	columns = SQLite.columns(table.db, table.tablename)

    col_expressions = Array{Expr,1}()
    for i in 1:size(columns,1)
    	type_string = get(columns[i,:type])
    	eltype = type_string=="INTEGER" ? Int64 :
    	         type_string[1:8]=="NVARCHAR" ? String :
    	         type_string=="DATETIME" ? DateTime : error("Unsupported column type")
    	if get(columns[i,:notnull])==0
    		eltype=Nullable{eltype}
    	end
        push!(col_expressions, Expr(:(::), get(columns[i,:name]), Type(eltype)))
    end
    t_expr = NamedTuples.make_tuple(col_expressions)
    #t_expr.args[1] = t_expr.args[1].args[1]

    t2 = :(Query.QueryableSQLite{Float64,Query.QueryProviderSQLite})
    t2.args[2] = t_expr

    eval(NamedTuples, :(import Query))
    t = eval(NamedTuples, t2)

    q = t(table.db, table.tablename)

    return q
end

# This must be a contender for most hacky, unstable function ever...
function collect(::Type{QueryProviderSQLite}, source)
    query_elements = []
    current_source = source
    while !isa(current_source, QueryableSQLite)
        push!(query_elements, current_source)
        current_source = current_source.source
    end

    where_count = count(i->isa(i,QueryableWhere),query_elements)
    where_clause = ""
    if where_count==1
    	where_o = filter(i->isa(i,QueryableWhere),query_elements)[1]
    	filter_el = where_o.filter
    	if filter_el.head!=:(->)
    		error()
    	end
    	name_of_index_var = filter_el.args[1]
    	operator = where_o.filter.args[2].args[1]
    	if operator!=:(==)
    		error()
    	end
    	if where_o.filter.args[2].args[2].head!=:(.) || where_o.filter.args[2].args[2].args[1]!=name_of_index_var
    		error()
    	end
    	left_hand_name = string(where_o.filter.args[2].args[2].args[2].value)
    	right_hand_value = where_o.filter.args[2].args[3]
    	where_clause = "WHERE $left_hand_name=\"$right_hand_value\""
    elseif where_count>1
    	error("At most one where clause is supported.")
    end

    select_count = count(i->isa(i,QueryableSelect),query_elements)
    select_clause = "*"
    if select_count==1
		select_o = filter(i->isa(i,QueryableSelect),query_elements)[1]
		f_el = select_o.f
		# Get columns we need
		col_names = Array(String,0)
		for i in 2:length(f_el.args[2].args)
			col_name = string(f_el.args[2].args[i].args[2].args[2].value)
			desired_name = string(f_el.args[2].args[i].args[1])
			push!(col_names, "$col_name AS $desired_name")
		end
		select_clause = join(col_names, ", ")
    elseif select_count>1
    	error("At most one select clause is supported.")
    end

    root = current_source
    sql_query = "SELECT $select_clause FROM $(root.tablename) $where_clause"
    println(sql_query)
    df = SQLite.query(root.db, sql_query)

    collect(query(df))
end

end


text:
function compute_nde_solution_history(datasets, NDEType, algorithm, nn_filepath, nn_history_filepath; gc_interval=50)
    final_nn = jldopen(nn_filepath, "r")

    Nz = final_nn["grid_points"]
    NN = final_nn["neural_network"]
    T_scaling = final_nn["T_scaling"]
    wT_scaling = final_nn["wT_scaling"]

    history = jldopen(nn_history_filepath, "r")
    epochs = keys(history["neural_network"]) |> length

    ids = keys(datasets)
    solution_history = Dict(id => [] for id in ids)

    for e in 1:epochs
        @info "Computing NDE solutions for epoch $e/$epochs..."

        NN = history["neural_network/$e"]

        for (id, ds) in datasets
            nde_sol = solve_nde(ds, NN, NDEType, algorithm, T_scaling, wT_scaling)
            push!(solution_history[id], nde_sol)
        end

        e % gc_interval == 0 && GC.gc() # Mercy if you're running this on a little laptop.
    end

    close(final_nn)
    close(history)

    return solution_history
end

function plot_epoch_loss(ids_train, ids_test, nde_solutions, true_solutions, T_scaling; title, filepath)
    ids = (ids_train..., ids_test...)
    epochs = length(nde_solutions[first(ids)])

    p = plot(dpi=200)

    ylims = (1e-1, 1e-1)
    kwargs = (linewidth=2, linealpha=0.8, yaxis=:log, ylabel="Mean squared error",
              xlabel="Epochs", title=title, grid=false, legend=:outertopright, framestyle=:box,
              foreground_color_legend=nothing, background_color_legend=nothing)

    for id in ids
        true_sol_scaled = T_scaling.(true_solutions[id].T)
        loss_history = [Flux.mse(true_sol_scaled, T_scaling.(nde_solutions[id][e].T)) for e in 1:epochs]

        min_loss, max_loss = extrema(loss_history)
        ylims_min = min(ylims[1], 10^floor(log10(min_loss)))
        ylims_max = max(ylims[2], 10^ceil(log10(max_loss)))
        ylims = (ylims_min, ylims_max)
        label = @sprintf("id=%d (%s)", id, id in ids_train ? "train" : "test")
        linestyle = id in ids_train ? :solid : :dot

        plot!(p, 1:epochs, loss_history, label=label, linestyle=linestyle, ylims=ylims, xlims=(1, epochs); kwargs...)
    end

    savefig(filepath)

    return p
end

function animate_nde_loss(datasets, ids_train, ids_test, nde_solutions, true_solutions, T_scaling; title, filepath, fps=15)
    ids = (ids_train..., ids_test...)
    epochs = length(nde_solutions[first(ids)])
    times = datasets[first(ids)]["T"].times

    ylims=(1e-6, 1)
    kwargs = (linewidth=2, linealpha=0.8, yaxis=:log, xlabel="Simulation time (days)",
              ylabel="Mean squared error", grid=false, legend=:outertopright, framestyle=:box,
              foreground_color_legend=nothing, background_color_legend=nothing)

    anim = @animate for e in 1:epochs
        @info "Plotting NDE loss evolution... epoch $e/$epochs"
        title_epoch = title * ": epoch $e"

        p = plot(dpi=200)

        for id in ids
            true_sol_scaled = T_scaling.(true_solutions[id].T)
            nde_sol_scaled = T_scaling.(nde_solutions[id][e].T)
            nde_loss = Flux.mse(true_sol_scaled, nde_sol_scaled, agg=x->mean(x, dims=1))[:]

            min_loss, max_loss = extrema(filter(!iszero, nde_loss))
            ylims = (min(ylims[1], 10^floor(log10(min_loss))), max(ylims[2], 10^ceil(log10(max_loss))))

            label = @sprintf("id=%d (%s)", id, id in ids_train ? "train" : "test")
            linestyle = id in ids_train ? :solid : :dot

            plot!(p, times, nde_loss, label=label, title=title_epoch, linestyle=linestyle,
                  xlims=extrema(times), ylims=ylims; kwargs...)
        end
    end

    @info "Saving $filepath..."
    mp4(anim, filepath * ".mp4", fps=fps)
    gif(anim, filepath * ".gif", fps=fps)

    return nothing
end

minimum_nonzero(xs...) = min([minimum(filter(!iszero, x)) for x in xs]...)
maximum_nonzero(xs...) = max([maximum(filter(!iszero, x)) for x in xs]...)

# TODO: Allow to selectively turn off different parameterizations.
function plot_comparisons(ds, id, ids_train, nde_sol, kpp_sol, tke_sol, convective_adjustment_sol, oceananigans_sol, T_scaling; filepath, frameskip=1, fps=15)
    T = ds["T"]
    wT = ds["wT"]
    FT = eltype(T)
    Nz = size(T, 3)
    Nt = size(T, 4)
    zc = znodes(T)
    zf = znodes(wT)
    times = T.times ./ days

    kwargs = (linewidth=2, linealpha=0.8, grid=false, framestyle=:box,
              foreground_color_legend=nothing, background_color_legend=nothing)

    T_lims = interior(T) |> extrema
    wT_lims = interior(wT) |> extrema

    loss(T, T̂) = Flux.mse(T_scaling.(T), T_scaling.(T̂))

    loss_nde = [loss(interior(T)[1, 1, :, n], nde_sol.T[:, n]) for n in 1:Nt]
    loss_kpp = [loss(interior(T)[1, 1, :, n], kpp_sol.T[:, n]) for n in 1:Nt]
    loss_tke = [loss(interior(T)[1, 1, :, n], tke_sol.T[:, n]) for n in 1:Nt]
    loss_ca = [loss(interior(T)[1, 1, :, n], convective_adjustment_sol.T[:, n]) for n in 1:Nt]
    loss_emb = [loss(interior(T)[1, 1, :, n], oceananigans_sol.T[:, n]) for n in 1:Nt]

    loss_min = minimum_nonzero(loss_nde, loss_kpp, loss_tke, loss_ca, loss_emb)
    loss_max = maximum_nonzero(loss_nde, loss_kpp, loss_tke, loss_ca, loss_emb)
    loss_extrema = (loss_min, loss_max)

    anim = @animate for n=1:frameskip:Nt
        @info "Plotting comparisons [frame $n/$Nt]: $filepath ..."

        time_str = @sprintf("%.2f days", times[n])
        title = @sprintf("Free convection (Q = %d W/m², %s): %s", -4e6 * ds.metadata["temperature_flux"], id in ids_train ? "train" : "test", time_str)

        wT_plot = plot(margin=5Plots.mm)
        T_plot = plot(margin=5Plots.mm)

        plot!(wT_plot, interior(wT)[1, 1, :, n], zf, xlims=wT_lims, label="LES", color="dodgerblue", ylims=extrema(zf),
              xlabel="Heat flux w′T′ (m/s ⋅ K)", ylabel="Depth z (meters)", legend=nothing; kwargs...)

        plot!(T_plot, interior(T)[1, 1, :, n], zc, label="LES", color="dodgerblue", xlabel="Temperature T (°C)",
              xlims=T_lims, ylims=extrema(zf), title=title, legend=:bottomright; kwargs...)

        plot!(T_plot, convective_adjustment_sol.T[:, n], zc, label="Convective adjustment", color="gray"; kwargs...)

        plot!(wT_plot, nde_sol.wT[:, n], zf, label="NDE", color="forestgreen"; kwargs...)
        plot!(T_plot, nde_sol.T[:, n], zc, label="NDE", color="forestgreen"; kwargs...)

        plot!(wT_plot, oceananigans_sol.wT[:, n], zf, label="Embedded", color="darkorange", linestyle=:dot; kwargs...)
        plot!(T_plot, oceananigans_sol.T[:, n], zc, label="Embedded", color="darkorange", linestyle=:dot; kwargs...)

        plot!(wT_plot, kpp_sol.wT[:, n], zf, label="KPP", color="crimson"; kwargs...)
        plot!(T_plot, kpp_sol.T[:, n], zc, label="KPP", color="crimson"; kwargs...)

        plot!(wT_plot, tke_sol.wT[:, n], zf, label="TKE mass flux", color="darkmagenta"; kwargs...)
        plot!(T_plot, tke_sol.T[:, n], zc, label="TKE mass flux", color="darkmagenta"; kwargs...)

        loss_plot = plot(margin=5Plots.mm)

        time_in_days = times[1:n]

        plot!(loss_plot, time_in_days, loss_nde[1:n], label="NDE", color="forestgreen",
              yaxis=(:log10, (1e-6, 1e-1)), xlims=extrema(times), ylims=(1e-6, 1e-1),
              xlabel="Time (days)", ylabel="Mean squared error", legend=nothing; kwargs...)

        plot!(loss_plot, time_in_days, loss_nde[1:n], label="NDE", color="forestgreen"; kwargs...)
        plot!(loss_plot, time_in_days, loss_emb[1:n], label="Embedded", color="darkorange", linestyle=:dot; kwargs...)
        plot!(loss_plot, time_in_days, loss_kpp[1:n], label="KPP", color="crimson"; kwargs...)
        plot!(loss_plot, time_in_days, loss_tke[1:n], label="TKE mass flux", color="darkmagenta"; kwargs...)

        plot(wT_plot, T_plot, loss_plot, layout=(1, 3), size=(1000, 400), dpi=200)
    end

    @info "Saving $filepath..."
    mp4(anim, filepath * ".mp4", fps=fps)
    gif(anim, filepath * ".gif", fps=fps)
end

# TODO: Allow to selectively turn off different parameterizations.
function plot_loss_matrix(datasets, ids_train, nde_sols, kpp_sols, tke_sols, convective_adjustment_sols, oceananigans_sols, T_scaling; filepath)

    loss(T, T̂) = Flux.mse(T_scaling.(T), T_scaling.(T̂))

    plots = []

    for (id, ds) in datasets

        T = ds["T"]
        wT = ds["wT"]
        Nz = size(T, 3)
        zc = znodes(T)
        zf = znodes(wT)
        H = abs(zf[1])
        Nt = size(T, 4)
        times = T.times ./ days

        kwargs = (linewidth=2, linealpha=0.8, grid=false, framestyle=:box, yaxis=(:log10, (1e-6, 1e-1)),
                  xlabel="Time (days)", ylabel="Mean squared error", xlims=extrema(times), ylims=(1e-6, 1e-1),
                  title = @sprintf("Q = %d W/m² (%s)", -4e6 * ds.metadata["temperature_flux"], id in ids_train ? "train" : "test"),
                  legend=isempty(plots) ? :bottomright : nothing,
                  foreground_color_legend=nothing, background_color_legend=nothing)

        loss_nde = [loss(interior(ds["T"])[1, 1, :, n], nde_sols[id].T[:, n]) for n in 1:Nt]
        loss_kpp = [loss(interior(ds["T"])[1, 1, :, n], kpp_sols[id].T[:, n]) for n in 1:Nt]
        loss_tke = [loss(interior(ds["T"])[1, 1, :, n], tke_sols[id].T[:, n]) for n in 1:Nt]
        loss_ca = [loss(interior(ds["T"])[1, 1, :, n], convective_adjustment_sols[id].T[:, n]) for n in 1:Nt]
        loss_emb = [loss(interior(ds["T"])[1, 1, :, n], oceananigans_sols[id].T[:, n]) for n in 1:Nt]

        p = plot()

        plot!(p, times, loss_ca, label="Convective adjustment", color="gray"; kwargs...)
        plot!(p, times, loss_nde, label="NDE", color="forestgreen"; kwargs...)
        plot!(p, times, loss_emb, label="Embedded", color="darkorange", linestyle=:dot; kwargs...)
        plot!(p, times, loss_kpp, label="KPP", color="crimson"; kwargs...)
        plot!(p, times, loss_tke, label="TKE mass flux", color="darkmagenta"; kwargs...)

        push!(plots, p)
    end

    P = plot(plots..., size=(800, 600), dpi=200)

    savefig(filepath)

    return P
end


text:
using VLKeggSDK

# reaction -
reaction_phrase = "C00533 + C00014 + 3 C00126 <=> C05361 + C00001 + 3 C00125"

# get symbols -
symbol_array = extract_metabolite_symbols(reaction_phrase)


text:
using Base.Test
import Ogre

@testset "Value sets" begin
    vs_notemp = Ogre.ValueSet(1,2,3)
    vs_full = Ogre.ValueSet(1,2,3,4)
    @test_throws MethodError Ogre.ValueSet(1,2)

    @testset "Type hierarchy" begin
        tempdep = Ogre.ValueSet{Ogre.WithTemp}
        notemp = Ogre.ValueSet{Ogre.NoTemp}
        @test !isa(vs_notemp, tempdep)
        @test isa(vs_notemp, notemp)
        @test isa(vs_notemp, Ogre.ValueSet)
        @test isa(vs_full, tempdep)
        @test !isa(vs_full, notemp)
        @test isa(vs_full, Ogre.ValueSet)
    end

    @testset "Get values from them" begin
        @test Ogre.mass(vs_notemp) == 1
        @test Ogre.nonmass(vs_notemp) == [2,3]
        @test Ogre.mass(vs_full) == 1
        @test Ogre.nonmass(vs_full) == [2,3,4]
        @test Ogre.radius(vs_notemp) == 2
        @test Ogre.pressure(vs_full) == 3
        @test_throws MethodError Ogre.temperature(vs_notemp)
        @test Ogre.temperature(vs_full) == 4
    end

    @testset "Zero values or less are considered 'unphysical'" begin
        @test !Ogre.isphysical(Ogre.zero(Ogre.ValueSet{Ogre.NoTemp}))
        @test !Ogre.isphysical(Ogre.zero(Ogre.ValueSet{Ogre.WithTemp}))
        @test !Ogre.isphysical(Ogre.ValueSet(0, 0, -1))
        @test !Ogre.isphysical(Ogre.ValueSet(0, 0, 0, -1))
    end
end


text:
# Experiment that should test J4PR functionality
module experiment3

reload("j4pr")
using DataFrames
import j4pr
# Generate data and write it to file
# TODO

# Read data and preproces
# -----------------------

# Read from csv
data = permutedims(readcsv("/home/zgornel/projects/j4pr/trunk/data/sample.csv"),(2,1));

label_col = 7

# Apply data transforms
# ---------------------

# 1. Process independently missing values on columns
data_filt = data |> (data |> j4pr.filterg(Dict(1=>"majority",
					2=>(x,y)->"missing",
					3=>(x,y)->zero(typeof(x)),
					4:5=>"majority",
					6=>(x,y)->"missing")
				   )
		 )

# 2. define a domain for a variable
dom_x1 = Dict(2=>x->x in ["0","1","a"] ? x : 0);
data_filt |> j4pr.filterdomain!(dom_x1);

# 3. Run one-hot encoder
data_ohenc = data_filt[1:6,:] |> j4pr.ohenc(j4pr._variable_(data_filt,1:6), Dict(1:6=>"binary"))
labels = Vector{String}(data_filt[7,:])
dc = j4pr.datacell(data_ohenc, labels)

end # module end


text:
"""
Created on Fri May 04 2020

@author: Yoann Pradat

    CentraleSupelec
    MICS laboratory
    9 rue Juliot Curie, Gif-Sur-Yvette, 91190 France

Algorithm for evaluating factors stability in rank selection method in NMF based on average stability of latent factor
across perturbed versions of the input matrix V.
"""

abstract type AbstractRSCluNMFResults <: AbstractNMFResults end

"""
    RSCluNMFResults(cluster_W_avg, cluster_W_std, cluster_H_avg, cluster_H_std, cluster_indices, cluster_stab_avg, cluster_intra_dist)

Struct for storing results returned by the clustering algorithm `evaluate_stability`.
"""
mutable struct RSCluNMFResults{T} <: AbstractRSCluNMFResults
    W_avg::Matrix{T}
    W_std::Matrix{T}
    H_avg::Matrix{T}
    H_std::Matrix{T}
    labels::Vector{Int64}
    stab_avg::Vector{Float64}
    stab_std::Vector{Float64}

    function RSCluNMFResults{T}(
        W_avg::Matrix{T},
        W_std::Matrix{T},
        H_avg::Matrix{T},
        H_std::Matrix{T},
        labels::Vector{Int64},
        stab_avg::Vector{Float64},
        stab_std::Vector{Float64}) where T <: Real

        new{T}(W_avg, W_std, H_avg, H_std, labels, stab_avg, stab_std)
    end
end

"""
    evaluate_stability(list_nmf_results[, dist_func, n_init, max_iter, stopping_tol, stopping_iter, rng])

Evaluates the stability of multiple runs of a NMF factorisation by clustering in a as many clusters
as the rank of the factorisation and compute the average silhouette index within each cluster.

#Arguments
- list_nmf_results::Vector{NMFResults{T}} : vector containing results of each NMF run
- dist_func::SemiMetric : SemiMetric from the `Distances` library used for computing distances between vectors.
- rng::AbstractRNG : random number generator, for the random initialization of the centroids positions.
- n_init::Integer : number of random initialization
- n_repl::Integer : maxium number of centroids update for each initialization.
"""

function evaluate_stability(list_nmf_results::Array{NMFResults{T},1}; dist_func::SemiMetric=CosineDist(), n_init::Integer=5, max_iter::Integer=100, stopping_crit::Symbol=:dist, stopping_tol::Real=0.005, stopping_iter::Integer=10, rng::AbstractRNG=MersenneTwister(0)) where T <: Real

    W_all = hcat([nmf_results.W for nmf_results in list_nmf_results]...)
    H_all = vcat([nmf_results.H for nmf_results in list_nmf_results]...)

    #### recover parameters of the algorithm
    n_nmf = length(list_nmf_results)
    n_tot = size(W_all, 2)
    n_fac = n_tot ÷ n_nmf

    list_centroids   = Array{Matrix{T}, 1}()
    list_labels      = Array{Array{Int64}, 1}()
    list_intra_dists = Array{Array{Float64,1}, 1}()
    list_cost        = Array{Float64,1}()

    perm_nmf = randperm(rng, n_nmf)

    for i_init = 1:min(n_init, n_nmf)
        i_iter = perm_nmf[i_init]
        init   = convert(Array{Int64,1}, ((i_iter-1) * n_fac + 1):(i_iter * n_fac))

        centroids, labels, cost = kmeans_init(
            X             = W_all,
            init          = init,
            max_iter      = max_iter,
            dist_func     = dist_func,
            stopping_crit = stopping_crit,
            stopping_tol  = stopping_tol,
            stopping_iter = stopping_iter
        )

        push!(list_centroids, centroids)
        push!(list_labels, labels)
        push!(list_cost, cost)
    end

    #### index of best initialization
    i_init_best = argmin(list_cost)

    centroids   = list_centroids[i_init_best]
    labels      = list_labels[i_init_best]

    #### compute silhouette indices per cluster
    if n_fac == 1
        stab_avg = [1.0]
        stab_std = [0.0]
    else
        stab_all = silhouettes(labels, pairwise(dist_func, W_all, dims=2))
        stab_avg = [mean(stab_all[findall(x -> x == label, labels)]) for label in sort(unique(labels))]
        stab_std = [std(stab_all[findall(x -> x == label, labels)]) for label in sort(unique(labels))]
    end

    #### rename labels in decreasing stab avg
    sorted_order = sortperm(stab_avg, rev=true)
    stab_avg = stab_avg[sorted_order]
    stab_std = stab_std[sorted_order]
    centroids = centroids[:,sorted_order]

    sorting_dict  = Dict( j => i for (i,j) in enumerate(sortperm(stab_avg, rev=true)))
    labels = map(x -> sorting_dict[x], labels)

    #### compute cluster standard deviations
    W_avg = zeros(Float64, size(centroids))
    W_std = zeros(Float64, size(W_avg))

    for label = 1:n_fac
        W_avg[:, label] = centroids[:, label]
        W_std[:, label] = std(W_all[:, labels .== label], dims=2)
    end

    #### compute cluster H avg and std
    H_avg = zeros(Float64, n_fac, size(H_all, 2))
    H_std = zeros(Float64, n_fac, size(H_all, 2))

    for label = 1:n_fac
        H_avg[label, :] = mean(H_all[labels .== label, :], dims = 1)
        H_std[label, :] = std(H_all[labels .== label, :], dims = 1)
    end

    RSCluNMFResults{T}(
        W_avg,
        W_std,
        H_avg,
        H_std,
        labels,
        stab_avg,
        stab_std
    )
end

#function evaluate_stability(list_nmf_results::Array{NMFResults{T},1}; dist_func::SemiMetric=CosineDist(),
#    rng::AbstractRNG=MersenneTwister(0), n_init::Integer=5, n_repl::Integer=100) where T <: Real
#    #### constants
#    convergence_cutoff = 0.005
#    convergence_iter = 10
#
#    W_all = hcat([nmf_results.W for nmf_results in list_nmf_results]...)
#    H_all = vcat([nmf_results.H for nmf_results in list_nmf_results]...)
#
#    #### recover parameters of the algorithm
#    n_clus = list_nmf_results[1].global_params.rank
#    n_iter = size(list_nmf_results, 1)
#    n_tot  = n_clus * n_iter
#
#    centroids_final = zeros(T, size(W_all, 1), n_clus)
#    cluster_indices = zeros(Int64, n_tot)
#    cluster_indices_final = zeros(Int64, n_tot)
#    cluster_intra_dist = zeros(Float64, n_iter, n_clus)
#    cluster_intra_dist_final = zeros(Float64, n_iter, n_clus)
#    min_cluster_intra_dist = Inf
#
#    #### random permutation for selection of initial centroids
#    #### i.e random selection of one matrix W among the n_iter
#    #### W matrices.
#    perm_iter = randperm(rng, n_iter)
#
#    for i_init = 1:min(n_init, n_iter)
#        centroids_cur = list_nmf_results[perm_iter[i_init]].W
#        centroids_pre = rand(rng, T, size(centroids_cur))
#
#        count_i_repl = 0
#
#        for i_repl = 1:n_repl
#            dist_all = pairwise(dist_func, hcat(centroids_cur, W_all), dims=2)
#            centroids_dist = transpose(dist_all[1:n_clus, (n_clus+1):(n_clus+n_tot)])
#
#            #### update cluster labels
#            perm_n_clus = randperm(rng, n_clus)
#            for j_n_clus in perm_n_clus
#                for i = 1:n_clus:n_tot
#                    i_range =  i:(i+n_clus-1)
#                    _, index_min = findmin(centroids_dist[i_range, j_n_clus], dims=1)
#                    centroids_dist[i_range[index_min],:] .= Inf
#                    cluster_indices[i_range[index_min]] .= j_n_clus
#                end
#            end
#
#            #### update centroids
#            dist_max_to_new_centroids = 0
#            for j_n_clus = 1:n_clus
#                centroids_cur[:, j_n_clus] = mean(W_all[:, cluster_indices .== j_n_clus], dims=2)
#                dist_to_new_centroids = dist_func(centroids_cur[:, j_n_clus], centroids_pre[:, j_n_clus])
#                dist_max_to_new_centroids = max(dist_max_to_new_centroids, dist_to_new_centroids)
#            end
#
#            #### check convergence
#            if dist_max_to_new_centroids < convergence_cutoff
#                count_i_repl += 1
#            else
#                count_i_repl = 0
#                centroids_pre .= centroids_cur
#            end
#
#            if count_i_repl == convergence_iter
#                break
#            end
#        end
#
#        #### update intra cluster distances
#        for j = 1:n_clus
#            cluster_intra_dist_one = pairwise(dist_func, hcat(centroids_cur[:,j], W_all[:, cluster_indices .== j]), dims=2)
#            cluster_intra_dist[:, j] = cluster_intra_dist_one[1, 2:size(cluster_intra_dist_one, 2)]
#        end
#
#        #### update final cluster centroids, indices and dist if lower
#        #### mean intra cluster distance is achieved
#        if min_cluster_intra_dist > mean(cluster_intra_dist)
#            centroids_final .= centroids_cur
#            cluster_indices_final .= cluster_indices
#            cluster_intra_dist_final .= cluster_intra_dist
#
#            min_cluster_intra_dist = mean(cluster_intra_dist)
#        end
#    end
#
#    centroids = centroids_final
#    cluster_indices = cluster_indices_final
#    cluster_intra_dist = cluster_intra_dist_final
#
#    #### reorder according to mean intra cluster distances
#    cluster_intra_dist_avg = vec(mean(cluster_intra_dist, dims=1))
#    cluster_intra_dist_avg_indices = sortperm(cluster_intra_dist_avg)
#
#    centroids = centroids[:, cluster_intra_dist_avg_indices]
#    cluster_intra_dist = cluster_intra_dist[:, cluster_intra_dist_avg_indices]
#    cluster_indices_new = cluster_indices[:]
#
#    for j = 1:n_clus
#        cluster_indices_new[cluster_indices .== cluster_intra_dist_avg_indices[j]] .= j
#    end
#    cluster_indices = cluster_indices_new[:]
#
#    #### compute silhouette indices per cluster
#    stab_all = silhouettes(cluster_indices, pairwise(dist_func, W_all, dims=2))
#    cluster_stab_avg = zeros(Float64, n_clus)
#
#    for j = 1:n_clus
#        cluster_stab_avg[j] = mean(stab_all[cluster_indices .== j])
#    end
#
#    #### compute cluster standard deviations
#    cluster_W_avg = zeros(Float64, size(centroids))
#    cluster_W_std = zeros(Float64, size(cluster_W_avg))
#
#    for j = 1:n_clus
#        cluster_W_avg[:, j] = centroids[:, j]
#        cluster_W_std[:, j] = std(W_all[:, cluster_indices .== j], dims=2)
#    end
#
#    #### compute cluster H avg and std
#    cluster_H_avg = zeros(Float64, n_clus, size(H_all, 2))
#    cluster_H_std = zeros(Float64, n_clus, size(H_all, 2))
#
#    for j = 1:n_clus
#        cluster_H_avg[j, :] = mean(H_all[cluster_indices .== j, :], dims = 1)
#        cluster_H_std[j, :] = std(H_all[cluster_indices .== j, :], dims = 1)
#    end
#
#    RSCluNMFResults{T}(
#        cluster_W_avg,
#        cluster_W_std,
#        cluster_H_avg,
#        cluster_H_std,
#        cluster_indices,
#        cluster_stab_avg,
#        cluster_intra_dist
#    )
#end

#mutable struct RSCluNMFResults{T} <: AbstractRSCluNMFResults
#    cluster_W_avg::Matrix{T}
#    cluster_W_std::Matrix{T}
#    cluster_H_avg::Matrix{T}
#    cluster_H_std::Matrix{T}
#    cluster_indices::Vector{Int64}
#    cluster_stab_avg::Vector{Float64}
#    cluster_intra_dist::Matrix{Float64}
#
#    function RSCluNMFResults{T}(
#        cluster_W_avg::Matrix{T},
#        cluster_W_std::Matrix{T},
#        cluster_H_avg::Matrix{T},
#        cluster_H_std::Matrix{T},
#        cluster_indices::Vector{Int64},
#        cluster_stab_avg::Vector{Float64},
#        cluster_intra_dist::Matrix{Float64}) where T <: Real
#
#        new{T}(cluster_W_avg, cluster_W_std, cluster_H_avg, cluster_H_std, cluster_indices, cluster_stab_avg, cluster_intra_dist)
#    end
#end


text:
####
# Support for Fractional derivatives
#
# defined as
####

export LeftIntegral,LeftDerivative, RightDerivative, RightIntegral

@calculus_operator LeftIntegral
@calculus_operator LeftDerivative

@calculus_operator RightIntegral
@calculus_operator RightDerivative

# DLMF18.17.9 with μ=0.5 and α=β=0

function LeftIntegral(S::Jacobi,k)
    if S.b==0
        ConcreteLeftIntegral(S,k)
    else
        J=Jacobi(S.a,0.,domain(S))
        LeftIntegralWrapper(LeftIntegral(J,k)*Conversion(S,J),k)
    end
end

function LeftIntegral(S::Ultraspherical,k)
    if order(S) == 1/2
        LeftIntegralWrapper(SpaceOperator(LeftIntegral(Jacobi(S),k),S,S),0.5)
    else
        J = Jacobi(S)
        LeftIntegralWrapper(LeftIntegral(J,k)*Conversion(S,J),0.5)
    end
end

LeftIntegral(S::Chebyshev,k) = LeftIntegralWrapper(
    LeftIntegral(Ultraspherical(1//2,domain(S)),k)*Conversion(S,Ultraspherical(1//2,domain(S))),
    0.5)

function rangespace{T,DD<:Interval}(Q::ConcreteLeftIntegral{Jacobi{T,DD},Float64})
    μ=Q.order
    S=domainspace(Q)

    JacobiWeight(S.b+μ,0.,Jacobi(S.a-μ,S.b+μ,domain(S)))
end

function RightIntegral(S::Jacobi,k)
    if S.a==0
        ConcreteRightIntegral(S,k)
    else
        J=Jacobi(0.,S.b,domain(S))
        RightIntegralWrapper(RightIntegral(J,k)*Conversion(S,J),k)
    end
end

function rangespace{T,DD<:Interval}(Q::ConcreteRightIntegral{Jacobi{T,DD},Float64})
    μ=Q.order
    S=domainspace(Q)
    @assert S.a==0

    JacobiWeight(0.,S.a+μ,Jacobi(S.a+μ,S.b-μ,domain(S)))
end

for TYP in (:ConcreteLeftIntegral,:ConcreteRightIntegral)
    @eval begin
        bandinds{T,DD<:Interval}(Q::$TYP{Jacobi{T,DD},Float64})=(0,0)
        getindex{T,DD<:Interval}(Q::$TYP{Jacobi{T,DD},Float64},k::Integer,j::Integer) =
            jacobi_frac_getindex(domain(Q),0.,Q.order,k,j)
    end
end

jacobi_frac_getindex(d::Interval,α,μ,k::Integer,j::Integer) =
    jacobi_frac_getindex((arclength(d)/2)^μ,α,μ,k,j)
jacobi_frac_getindex(c::Number,α,μ,k::Integer,j::Integer) =
    k==j ? c*exp(lgamma(α+k)-lgamma(α+μ+k)) : zero(promote_type(typeof(c),typeof(α),typeof(μ)))

# jacobi_frac_addentries!(d::Interval,α,μ,A,kr::UnitRange)=
#     jacobi_frac_addentries!((arclength(d)/2)^μ,α,μ,A,kr)
# function jacobi_frac_addentries!(c::Number,α,μ,A,kr::UnitRange)
#     γ=c*gamma(α+1)/gamma(α+1+μ)
#     for k=1:first(kr)-1
#         γ*=(α+k)/(α+μ+k)
#     end
#     for k=kr
#        A[k,k]+=γ
#        γ*=(α+k)/(α+μ+k)
#         #should be gamma(S.α+k)/gamma(S.α+μ+k)=
#         #(S.α+k-1)/(S.α+μ+k-1)*gamma(S.α+k-1)/gamma(S.α+μ+k-1)
#     end
#     A
# end

function LeftIntegral{T,DD}(S::JacobiWeight{Jacobi{T,DD}},k)
    J=S.space
    @assert S.β==0
    @assert S.α==J.b
    ConcreteLeftIntegral(S,k)
end

function RightIntegral{T,DD}(S::JacobiWeight{Jacobi{T,DD}},k)
    J=S.space
    @assert S.β==J.a
    @assert S.α==0
    ConcreteRightIntegral(S,k)
end

function LeftIntegral{DD}(S::JacobiWeight{Chebyshev{DD}},k)
    # convert to Jacobi
    Q=LeftIntegral(JacobiWeight(S.α,S.β,Jacobi(S.space)),k)
    LeftIntegralWrapper(Q*Conversion(S,domainspace(Q)),k)
end

function RightIntegral{DD}(S::JacobiWeight{Chebyshev{DD}},k)
    # convert to Jacobi
    Q=RightIntegral(JacobiWeight(S.α,S.β,Jacobi(S.space)),k)
    RightIntegralWrapper(Q*Conversion(S,domainspace(Q)),k)
end

for (TYP,WRAP) in ((:LeftIntegral,:LeftIntegralWrapper),
                    (:RightIntegral,:RightIntegralWrapper))
    @eval function $TYP{PS<:PolynomialSpace}(S::JacobiWeight{PS},k)
        JS = JacobiWeight(S.α,S.β,Jacobi(S.space))
        $WRAP($TYP(JS,k)*Conversion(S,JS),k)
    end
end

#DLMF18.17.9
function rangespace{T,DD<:Interval}(Q::ConcreteLeftIntegral{JacobiWeight{Jacobi{T,DD},DD},Float64})
    μ=Q.order
    S=domainspace(Q)
    J=S.space
    if isapprox(S.α,-μ)
        Jacobi(J.a-μ,J.b+μ,domain(J))
    else
        JacobiWeight(S.α+μ,0.,Jacobi(J.a-μ,J.b+μ,domain(J)))
    end
end

function rangespace{T,DD<:Interval}(Q::ConcreteRightIntegral{JacobiWeight{Jacobi{T,DD},DD},Float64})
    μ=Q.order
    S=domainspace(Q)
    J=S.space
    @assert S.β==J.a
    @assert S.α==0
    if isapprox(S.β,-μ)
        Jacobi(J.a+μ,J.b-μ,domain(J))
    else
        JacobiWeight(0.,S.β+μ,Jacobi(J.a+μ,J.b-μ,domain(J)))
    end
end

for TYP in (:ConcreteLeftIntegral,:ConcreteRightIntegral)
    @eval bandinds{T,DD<:Interval}(Q::$TYP{JacobiWeight{Jacobi{T,DD},DD},Float64})=(0,0)
end

getindex{T,DD<:Interval}(Q::ConcreteLeftIntegral{JacobiWeight{Jacobi{T,DD},DD},Float64},k::Integer,j::Integer) =
    jacobi_frac_getindex(domain(Q),domainspace(Q).α,Q.order,k,j)

getindex{T,DD<:Interval}(Q::ConcreteRightIntegral{JacobiWeight{Jacobi{T,DD},DD},Float64},k::Integer,j::Integer) =
    jacobi_frac_getindex(domain(Q),domainspace(Q).β,Q.order,k,j)

function choosedomainspace{T<:Float64}(Q::LeftIntegral{UnsetSpace,T},sp::JacobiWeight)
    @assert sp.α>0 && isapproxinteger(sp.β)

    if isapprox(Q.order,sp.α)
        Legendre(domain(sp))
    else
        JacobiWeight(sp.α-Q.order,0.,Jacobi(0.,sp.α-Q.order,domain(sp)))
    end
end

choosedomainspace{T<:Float64}(Q::LeftIntegral{UnsetSpace,T},sp::PolynomialSpace)=
                JacobiWeight(-Q.order,0.,Jacobi(-Q.order,-Q.order,domain(sp)))

function choosedomainspace{T<:Float64}(Q::RightIntegral{UnsetSpace,T},sp::JacobiWeight)
    @assert sp.β>0  && isapproxinteger(sp.α)

    if isapprox(Q.order,sp.β)
        Legendre(domain(sp))
    else
        JacobiWeight(0.,sp.β-Q.order,Jacobi(sp.β-Q.order,0.))
    end
end
choosedomainspace{T<:Float64}(Q::RightIntegral{UnsetSpace,T},sp::PolynomialSpace)=
    JacobiWeight(0.,-Q.order,Jacobi(-Q.order,-Q.order,domain(sp)))

# Define Left/RightDerivative
for (DTYP,QTYP,DWRAP,QWRAP) in ((:LeftDerivative,:LeftIntegral,:LeftDerivativeWrapper,:LeftIntegralWrapper),
                            (:RightDerivative,:RightIntegral,:RightDerivativeWrapper,:RightIntegralWrapper))
    @eval begin
        function $DTYP(S::Space,k::Real)
            i=ceil(Int,k)
            r=i-k
            $DWRAP(i<0?$QTYP(S,-k):Derivative(i)*$QTYP(S,r),k)
        end
        $QTYP(S::SumSpace,k)=$QWRAP(DiagonalInterlaceOperator(map(s->$QTYP(s,k),S.spaces),SumSpace),k)
    end
end


text:
using Plots
rosenbrock(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

linspace(a,b,N)= (a+ (b-a)*(i)/(N-1) for i=0:(N-1))

x = y = sort(collect(Iterators.flatten([linspace(-1,-0.1,25), linspace(0,2,1000), linspace(2.1,3,25)])))

dpi = 600

l = Vector{Float64}()
for i in x
   for j in y
       v = [i j]
       append!(l, rosenbrock(v))
   end
end

surface(x,y,l, dpi=dpi)
png("./rosenbrock_surface.png")
contour(x,y,l, dpi=dpi)
png("./rosenbrock_contour.png")


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

#using Compat

const DNN_MAX_DIMENSION = 32
const DNN_QUERY_MAX_LENGTH = 128

typealias _uniPrimitive_s Void
typealias dnnPrimitive_t Ptr{_uniPrimitive_s}
typealias _dnnLayout_s Void
typealias dnnLayout_t Ptr{_dnnLayout_s}
typealias dnnPrimitiveAttributes_t Ptr{Cvoid}

# begin enum ANONYMOUS_1
typealias ANONYMOUS_1 Cint
const E_SUCCESS = (Int32)(0)
const E_INCORRECT_INPUT_PARAMETER = (Int32)(-1)
const E_UNEXPECTED_NULL_POINTER = (Int32)(-2)
const E_MEMORY_ERROR = (Int32)(-3)
const E_UNSUPPORTED_DIMENSION = (Int32)(-4)
const E_UNIMPLEMENTED = (Int32)(-127)
# end enum ANONYMOUS_1

# begin enum dnnError_t
typealias dnnError_t Cint
const E_SUCCESS = (Int32)(0)
const E_INCORRECT_INPUT_PARAMETER = (Int32)(-1)
const E_UNEXPECTED_NULL_POINTER = (Int32)(-2)
const E_MEMORY_ERROR = (Int32)(-3)
const E_UNSUPPORTED_DIMENSION = (Int32)(-4)
const E_UNIMPLEMENTED = (Int32)(-127)
# end enum dnnError_t

# begin enum ANONYMOUS_2
typealias ANONYMOUS_2 UInt32
const dnnAlgorithmConvolutionGemm = (UInt32)(0)
const dnnAlgorithmConvolutionDirect = (UInt32)(1)
const dnnAlgorithmConvolutionFFT = (UInt32)(2)
const dnnAlgorithmPoolingMax = (UInt32)(3)
const dnnAlgorithmPoolingMin = (UInt32)(4)
const dnnAlgorithmPoolingAvg = (UInt32)(5)
# end enum ANONYMOUS_2

# begin enum dnnAlgorithm_t
typealias dnnAlgorithm_t UInt32
const dnnAlgorithmConvolutionGemm = (UInt32)(0)
const dnnAlgorithmConvolutionDirect = (UInt32)(1)
const dnnAlgorithmConvolutionFFT = (UInt32)(2)
const dnnAlgorithmPoolingMax = (UInt32)(3)
const dnnAlgorithmPoolingMin = (UInt32)(4)
const dnnAlgorithmPoolingAvg = (UInt32)(5)
# end enum dnnAlgorithm_t

# begin enum ANONYMOUS_3
typealias ANONYMOUS_3 UInt32
const dnnResourceSrc = (UInt32)(0)
const dnnResourceFrom = (UInt32)(0)
const dnnResourceDst = (UInt32)(1)
const dnnResourceTo = (UInt32)(1)
const dnnResourceFilter = (UInt32)(2)
const dnnResourceScaleShift = (UInt32)(2)
const dnnResourceBias = (UInt32)(3)
const dnnResourceDiffSrc = (UInt32)(4)
const dnnResourceDiffFilter = (UInt32)(5)
const dnnResourceDiffScaleShift = (UInt32)(5)
const dnnResourceDiffBias = (UInt32)(6)
const dnnResourceDiffDst = (UInt32)(7)
const dnnResourceWorkspace = (UInt32)(8)
const dnnResourceMultipleSrc = (UInt32)(16)
const dnnResourceMultipleDst = (UInt32)(24)
const dnnResourceNumber = (UInt32)(32)
# end enum ANONYMOUS_3

# begin enum dnnResourceType_t
typealias dnnResourceType_t UInt32
const dnnResourceSrc = (UInt32)(0)
const dnnResourceFrom = (UInt32)(0)
const dnnResourceDst = (UInt32)(1)
const dnnResourceTo = (UInt32)(1)
const dnnResourceFilter = (UInt32)(2)
const dnnResourceScaleShift = (UInt32)(2)
const dnnResourceBias = (UInt32)(3)
const dnnResourceDiffSrc = (UInt32)(4)
const dnnResourceDiffFilter = (UInt32)(5)
const dnnResourceDiffScaleShift = (UInt32)(5)
const dnnResourceDiffBias = (UInt32)(6)
const dnnResourceDiffDst = (UInt32)(7)
const dnnResourceWorkspace = (UInt32)(8)
const dnnResourceMultipleSrc = (UInt32)(16)
const dnnResourceMultipleDst = (UInt32)(24)
const dnnResourceNumber = (UInt32)(32)
# end enum dnnResourceType_t

# begin enum ANONYMOUS_4
typealias ANONYMOUS_4 UInt32
const dnnBorderZeros = (UInt32)(0)
const dnnBorderExtrapolation = (UInt32)(3)
# end enum ANONYMOUS_4

# begin enum dnnBorder_t
typealias dnnBorder_t UInt32
const dnnBorderZeros = (UInt32)(0)
const dnnBorderExtrapolation = (UInt32)(3)
# end enum dnnBorder_t


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
using CSV, DataFrames, Pipe

languagesF = "../data/languages.csv"
formsF = "../data/forms.csv"

!(isfile(languagesF) && isfile(formsF)) && begin
      asjpZip = download(
            "https://zenodo.org/record/3843469/files/lexibank/asjp-v19.1.zip",
            "tmp/asjp-19.1.zip",
      )
      run(`unzip -o $asjpZip -d tmp/`)
      cp("tmp/lexibank-asjp-0c18d44/cldf/forms.csv", formsF)
      cp("tmp/lexibank-asjp-0c18d44/cldf/languages.csv", languagesF)
      for f in readdir("tmp")
            rm("tmp/" * f, recursive = true)
      end
end

forms = CSV.read(formsF)
languages = CSV.read(languagesF)

function cleanASJP(word)
      @pipe word |>
            replace(_, r"[ \*~\"]" => "") |>
            replace(_, r"(.)(.)(.)\$" => s"\2")
end

languages[!, :longname] = @pipe languages |>
      zip(_.classification_wals, _.Name) |>
      join.(_, ".") |>
      replace.(_, "-" => "_")

forms[!, :simplified] = cleanASJP.(forms.Value)

languages = languages[.!occursin.("Oth.", languages.classification_wals), :]

conceptCoverage = @pipe forms |>
      unique(_, [:Language_ID, :gloss_in_source]) |>
      groupby(_, :gloss_in_source) |>
      combine(nrow, _) |>
      sort(_, :nrow, rev = true)

concepts = conceptCoverage.gloss_in_source[1:40]

forms40 = forms[map(x -> x ∈ concepts, forms.gloss_in_source), :]

languageCoverage = @pipe forms40 |>
      unique(_, [:Language_ID, :gloss_in_source]) |>
      groupby(_, :Language_ID) |>
      combine(nrow, _) |>
      sort(_, :nrow, rev = true)

doculects = languageCoverage.Language_ID[languageCoverage.nrow.>=30]

forms40 = forms40[map(x -> x ∈ doculects, forms40.Language_ID), :]

asjpLong = innerjoin(
      forms40,
      languages[:, [:Name, :longname]],
      on = :Language_ID => :Name,
)

asjpWide = @pipe asjpLong |>
      groupby(_, [:longname, :gloss_in_source]) |>
      combine(:simplified => x -> join(x, "-"), _) |>
      unstack(_, :longname, :gloss_in_source, :simplified_function)

CSV.write("../data/asjp19wide.csv", asjpWide)


text:
#=display.jl - will replay the game
from start to finish, displaying the current game state.
Accepts 1 command line argument,<filename> => database
=#
include("board.jl")
#include("dParse.jl")
using BM

#include("start.jl")
using SQLite
pathToDatabase = ARGS[1]
db = SQLite.DB(pathToDatabase) #Opens the database gamefile
global board = BM.startGame(gameType)
#calculates all black piece positions
res = SQLite.query(db,"""SELECT "value" FROM meta WHERE "key" = 'type';""")
gameType = get(res[1,1])
#=----Gets last move number----=#
res = SQLite.query(db,"""SELECT MAX("move_number") FROM moves;""") #Finds the last played move (maximum move_number)
lastMoveIDNullable = res[1,1] #SQL query with max move_number (POSSIBLY "NULL" if no moves have been made)

if (!isnull(lastMoveIDNullable)) #Checks that lastMoveID was not NULL
lastMoveID = get(res[1,1]) #Parses the last move move_number as an Int

else #lastMoveID is NULL
lastMoveID = 0 #move_number "0" is unsused. This implies no moves have been made

end
currentMoveID=1

#=Helper functions=#
function displaymini()
  board = BM.startGame(gameType)
  global dboard = Array(Char,9,9)
  for (x,y) in keys(board)
    dboard[(x,y)]=board[(x,y)][1]
  end

  for y_index in (1:11)
    for x_index in (1:11)
      if x_index==1
        if y_index==1
          print("┌")
        elseif y_index==11
          print("└")
        elseif rem(y_index,2)==0
          print("|")
        else
          print("├")
        end
      elseif x_index==11
        if y_index==1
          print("┐")
        elseif y_index==11
          print("┘")
        elseif rem(y_index,2)==0
          print("|")
        else
          print("┤")
        end
      elseif rem(x_index,2)==1
        if y_index==1
          print("┬")
        elseif y_index==11
          print("┴")
        elseif rem(y_index,2)==1
          print("┼")
        else
          print("|")
        end
      else
        if rem(y_index,2)==1
          print("-")
        else
          if board[div(y_index,2),(6-div(x_index,2))]=="king"
            #print_with_color(:yellow, "$(board[div(y_index,2),(6-div(x_index,2))])")
            continue
          end
          #print(board[div(y_index,2),(6-div(x_index,2))])
        end
      end
    end
    print("\n")
  end
end

function displaystandard()

  board = BM.startGame(gameType)
  global dboard = Array(Char,9,9)
  for (x,y) in keys(board)
    dboard[(x,y)]=board[(x,y)][1]
  end

for y_index in (1:19)
  for x_index in (1:19)
    if x_index==1
      if y_index==1
        print("┌")
      elseif y_index==19
        print("└")
      elseif rem(y_index,2)==0
        print("|")
      else
        print("├")
      end
    elseif x_index==19
      if y_index==1
        print("┐")
      elseif y_index==19
        print("┘")
      elseif rem(y_index,2)==0
        print("|")
      else
        print("┤")
      end
    elseif rem(x_index,2)==1
      if y_index==1
        print("┬")
      elseif y_index==19
        print("┴")
      elseif rem(y_index,2)==1
        print("┼")
      else
        print("|")
      end
    else
      if rem(y_index,2)==1
        print("-")
      else
        if board[div(y_index,2),(10-div(x_index,2))]=='k' || board[div(y_index,2),(10-div(x_index,2))]=='K'
        #  print_with_color(:yellow, "$(dboard[div(y_index,2),(10-div(x_index,2))])")
          continue
        end
      #  print(dboard[div(y_index,2),(10-div(x_index,2))])
      end
    end
  end
  print("\n")
end

end#Func

function displaychu()
  board = BM.startGame(gameType)
  global dboard = Array(Char,9,9)
  for (x,y) in keys(board)
    dboard[(x,y)]=board[(x,y)][1]
  end

  for y_index in (1:25)
    for x_index in (1:25)
      if x_index==1
        if y_index==1
          print("┌")
        elseif y_index==25
          print("└")
        elseif rem(y_index,2)==0
          print("|")
        else
          print("├")
        end
      elseif x_index==25
        if y_index==1
          print("┐")
        elseif y_index==25
          print("┘")
        elseif rem(y_index,2)==0
          print("|")
        else
          print("┤")
        end
      elseif rem(x_index,2)==1
        if y_index==1
          print("┬")
        elseif y_index==25
          print("┴")
        elseif rem(y_index,2)==1
          print("┼")
        else
          print("|")
        end
      else
        if rem(y_index,2)==1
          print("-")
        else
          #if board[div(y_index,2),(13-div(x_index,2))]=='k' || board[div(y_index,2),(13-div(x_index,2))]=='K'
          #  print_with_color(:yellow, "$(dboard[div(y_index,2),(10-div(x_index,2))])")
          #  continue
        #  end
        #  print(dboard[div(y_index,2),(13-div(x_index,2))])
        end
      end
    end
    print("\n")
  end
end

function displaytenjiku()

  for y_index in (1:33)
    for x_index in (1:33)
      if x_index==1
        if y_index==1
          print("┌")
        elseif y_index==33
          print("└")
        elseif rem(y_index,2)==0
          print("|")
        else
          print("├")
        end
      elseif x_index==33
        if y_index==1
          print("┐")
        elseif y_index==33
          print("┘")
        elseif rem(y_index,2)==0
          print("|")
        else
          print("┤")
        end
      elseif rem(x_index,2)==1
        if y_index==1
          print("┬")
        elseif y_index==33
          print("┴")
        elseif rem(y_index,2)==1
          print("┼")
        else
          print("|")
        end
      else
        if rem(y_index,2)==1
          print("-")
        else
          #if board[div(y_index,2),(13-div(x_index,2))]=='k' || board[div(y_index,2),(13-div(x_index,2))]=='K'
          #  print_with_color(:yellow, "$(dboard[div(y_index,2),(10-div(x_index,2))])")
          #  continue
        #  end
        #  print(dboard[div(y_index,2),(13-div(x_index,2))])
        end
      end
    end
    print("\n")
  end
end

#=----=#
#=----Replays the game until move_id = lastMoveID----=#
while currentMoveID<=lastMoveID
  global board
  res = SQLite.query(db,"""SELECT sourcex,sourcey,targetx,targety,move_type,option FROM moves WHERE "move_number" = $(currentMoveID);""")
  sourcexNullable = res[1,1]
  sourceyNullable = res[1,2]
  targetxNullable = res[1,3]
  targetyNullable = res[1,4]
  move_type = get(res[1,5])
  optionNullable = res[1,6]

  if move_type == "move" #Regular Move
    targetx = get(targetxNullable)
    targety = get(targetyNullable)
    sourcex = get(sourcexNullable)
    sourcey = get(sourceyNullable)
    if !(ST.isEmpty(board[targetx,targety]))# capture
      #push(captures,board[targetx][targety])
    end
    board[targetx,targety].piece = board[sourcex,sourcey].piece
    board[targetx,targety].team = board[sourcex,sourcey].team
    ST.clear!(board[sourcex,sourcey])

  elseif move_type == "drop"
    option = get(optionNullable)
    try
      #deleteat!(captures,findfirst(x->x.piece==option))
    end
    targetx = get(targetxNullable)
    targety = get(targetyNullable)
    board[targetx,targety].piece = board[sourcex,sourcey].piece
    board[targetx,targety].team = board[sourcex,sourcey].team
  elseif move_type == "resign"
    #Do nothing
  end
  ST.saveBoard(board)

  if gameType == "standard"
    displaystandard()
  else
    displaymini()
  end
  sleep(1) #Waits half a second before running displaying a new board
  currentMoveID+=1
end

max_stars_repo_name: kinoroy/Shogi

max_issues_repo_name: kinoroy/Shulia

max_forks_repo_name: kinoroy/Shulia


text:
using BenchmarkTools

"""
    read_reverse_file(filename::String)

       Read an entire file and write its content in
       a new file in reverse order.
"""
function read_reverse_file(filename::String)
    open(filename) do file_id
        lines = readlines(file_id)
    end

    rev_filename = "rev_"+basename(filename)
    open(rev_filename, "w") do rev_file_id
        println(rev_file_id, [reverse(line) for line in reverse(lines)])
    end
end

# Get the text file name from the command line.
n, = size(ARGS)
if n < 1
    println("Usage: ")
    println("    julia " + ARG[0] + " filename")
    println("       --->      Please specify the filename.")
    exit()
end

filename = ARGS[1]
@btime read_reverse_file(filename)

println(" ")


text:
import Base.rand

"""
    rand(N::ProbabilisticNetwork)

Converts a probabilistic network into a deterministic one by performing random
draws. All interactions are treated as independent Bernoulli events. Note that
this network is *not* check for degeneracy, *i.e.* species can end up with no
interactions.
"""
function rand(N::ProbabilisticNetwork)
    # Get the correct network type
    newtype = typeof(N) <: AbstractUnipartiteNetwork ? UnipartiteNetwork : BipartiteNetwork
    return newtype(rand(Float64, size(N)).<=N.A, species_objects(N)...)
end

"""
    rand(N::ProbabilisticNetwork, n::T) where {T<:Integer}

Generates a number of random deterministic networks based on a probabilistic
network.
"""
function rand(N::ProbabilisticNetwork, n::T) where {T<:Integer}
    @assert n > 0
    return map(x -> rand(N), 1:n)
end


text:
#
# SoapySDR logger
# https://github.com/pothosware/SoapySDR/blob/1cf5a539a21414ff509ff7d0eedfc5fa8edb90c6/include/SoapySDR/Logger.h
#

@enum SoapySDRLogLevel begin
    SOAPY_SDR_FATAL    = 1 #!< A fatal error. The application will most likely terminate. This is the highest priority.
    SOAPY_SDR_CRITICAL = 2 #!< A critical error. The application might not be able to continue running successfully.
    SOAPY_SDR_ERROR    = 3 #!< An error. An operation did not complete successfully, but the application as a whole is not affected.
    SOAPY_SDR_WARNING  = 4 #!< A warning. An operation completed with an unexpected result.
    SOAPY_SDR_NOTICE   = 5 #!< A notice, which is an information with just a higher priority.
    SOAPY_SDR_INFO     = 6 #!< An informational message, usually denoting the successful completion of an operation.
    SOAPY_SDR_DEBUG    = 7 #!< A debugging message.
    SOAPY_SDR_TRACE    = 8 #!< A tracing message. This is the lowest priority.
    SOAPY_SDR_SSI      = 9 #!< Streaming status indicators such as "U" (underflow) and "O" (overflow).
end

#"""
#Send a message to the registered logger.
#\param logLevel a possible logging level
#\param message a logger message string
#"""
#SOAPY_SDR_API void SoapySDR_log(const SoapySDRLogLevel logLevel, const char *message);

#"""
#Send a message to the registered logger.
#\param logLevel a possible logging level
#\param format a printf style format string
#\param argList an argument list for the formatter
#"""
#SOAPY_SDR_API void SoapySDR_vlogf(const SoapySDRLogLevel logLevel, const char *format, va_list argList);

#"""
#Send a message to the registered logger.
#\param logLevel a possible logging level
#\param format a printf style format string
#"""
#static inline void SoapySDR_logf(const SoapySDRLogLevel logLevel, const char *format, ...)
#{
#    va_list argList;
#    va_start(argList, format);
#    SoapySDR_vlogf(logLevel, format, argList);
#    va_end(argList);
#}

#"""
#Typedef for the registered log handler function.
#"""
#typedef void (*SoapySDRLogHandler)(const SoapySDRLogLevel logLevel, const char *message);

"""
Register a new system log handler.
Platforms should call this to replace the default stdio handler.
Passing `NULL` restores the default.
"""
function SoapySDR_registerLogHandler(func)
    #SOAPY_SDR_API void SoapySDR_registerLogHandler(const SoapySDRLogHandler handler);
    ccall((:SoapySDR_registerLogHandler, lib), Cvoid, (Ptr{Cvoid},), func)
end

"""
Set the log level threshold.
Log messages with lower priority are dropped.
"""
function SoapySDR_setLogLevel(level::Cint)
    #SOAPY_SDR_API void SoapySDR_setLogLevel(const SoapySDRLogLevel logLevel);
end


text:
##### Multi-physics types

using CLIMAParameters.Planet: T_freeze, cv_l

using CloudMicrophysics.Microphysics_0M
using CloudMicrophysics.Microphysics_1M
import CloudMicrophysics

const CM1M = CloudMicrophysics.Microphysics_1M

export RainSnow_1M
export WarmRain_1M
export Diffusion
export Subsidence
export RemovePrecipitation

struct Subsidence{FT} <: TendencyDef{Source}
    D::FT
end

prognostic_vars(::Subsidence) = (Mass(), Energy(), TotalMoisture())

# Subsidence includes tendencies in Mass, Energy and TotalMoisture equations:

subsidence_velocity(subsidence::Subsidence{FT}, z::FT) where {FT} =
    -subsidence.D * z

struct PressureGradient <: TendencyDef{Flux{FirstOrder}} end
struct Pressure <: TendencyDef{Flux{FirstOrder}} end
struct Advect <: TendencyDef{Flux{FirstOrder}} end
struct Diffusion <: TendencyDef{Flux{SecondOrder}} end
struct MoistureDiffusion <: TendencyDef{Flux{SecondOrder}} end

"""
    RemovePrecipitation <: TendencyDef{Source}

A sink to `q_tot` when cloud condensate is exceeding a threshold.
The threshold is defined either in terms of condensate or supersaturation.
The removal rate is implemented as a relaxation term
in the CloudMicrophysics.jl Microphysics_0M module.
The default thresholds and timescale are defined in CLIMAParameters.jl.
"""
struct RemovePrecipitation <: TendencyDef{Source}
    " Set to true if using qc based threshold"
    use_qc_thr::Bool
end

prognostic_vars(::RemovePrecipitation) = (Mass(), Energy(), TotalMoisture())

"""
    PrecipitationFlux <: TendencyDef{Flux{FirstOrder}}

Computes the precipitation flux as a sum of air velocity and terminal velocity
multiplied by the advected variable.
"""
struct PrecipitationFlux <: TendencyDef{Flux{FirstOrder}} end

function remove_precipitation_sources(s::RemovePrecipitation, atmos, args)
    @unpack state, aux = args
    @unpack ts = args.precomputed

    FT = eltype(state)

    q = PhasePartition(ts)
    λ::FT = liquid_fraction(ts)
    I_l::FT = internal_energy_liquid(ts)
    I_i::FT = internal_energy_ice(ts)
    Φ::FT = gravitational_potential(atmos.orientation, aux)

    S_qt::FT = 0
    param_set = parameter_set(atmos)
    if s.use_qc_thr
        S_qt = remove_precipitation(param_set, q)
    else
        q_vap_sat = q_vap_saturation(ts)
        S_qt = remove_precipitation(param_set, q, q_vap_sat)
    end

    S_e::FT = (λ * I_l + (1 - λ) * I_i + Φ) * S_qt

    return (; S_ρ_qt = state.ρ * S_qt, S_ρ_e = state.ρ * S_e)
end

"""
    WarmRain_1M <: TendencyDef{Source}

A collection of source/sink terms related to 1-moment warm rain microphysics.
The microphysics process rates are implemented
in the CloudMicrophysics.jl Microphysics_1M module.
"""
struct WarmRain_1M <: TendencyDef{Source} end

prognostic_vars(::WarmRain_1M) =
    (Mass(), Energy(), TotalMoisture(), LiquidMoisture(), Rain())

function warm_rain_sources(atmos, args, ts)
    @unpack state, aux = args

    FT = eltype(state)

    q = PhasePartition(ts)
    T::FT = air_temperature(ts)
    I_l::FT = internal_energy_liquid(ts)
    q_rai::FT = state.precipitation.ρq_rai / state.ρ
    Φ::FT = gravitational_potential(atmos.orientation, aux)
    param_set = parameter_set(atmos)

    # autoconversion
    src_q_rai_acnv = conv_q_liq_to_q_rai(param_set, q.liq)
    # accretion
    src_q_rai_accr = accretion(
        param_set,
        CM1M.LiquidType(),
        CM1M.RainType(),
        q.liq,
        q_rai,
        state.ρ,
    )
    # rain evaporation
    src_q_rai_evap = evaporation_sublimation(
        param_set,
        CM1M.RainType(),
        q,
        q_rai,
        state.ρ,
        T,
    )

    S_qr::FT = src_q_rai_acnv + src_q_rai_accr + src_q_rai_evap
    S_ql::FT = -src_q_rai_acnv - src_q_rai_accr
    S_qt::FT = -S_qr
    S_e::FT = S_qt * (I_l + Φ)

    return (;
        S_ρ_qt = state.ρ * S_qt,
        S_ρ_ql = state.ρ * S_ql,
        S_ρ_qr = state.ρ * S_qr,
        S_ρ_e = state.ρ * S_e,
    )
end

"""
    RainSnow_1M <: TendencyDef{Source}

A collection of source/sink terms related to 1-moment rain and snow microphysics
The microphysics process rates are implemented
in the CloudMicrophysics.jl Microphysics_1M module.
"""
struct RainSnow_1M <: TendencyDef{Source} end

prognostic_vars(::RainSnow_1M) = (
    Mass(),
    Energy(),
    TotalMoisture(),
    LiquidMoisture(),
    IceMoisture(),
    Rain(),
    Snow(),
)

function rain_snow_sources(atmos, args, ts)
    @unpack state, aux = args

    FT = eltype(state)
    param_set = parameter_set(atmos)

    q_rai::FT = state.precipitation.ρq_rai / state.ρ
    q_sno::FT = state.precipitation.ρq_sno / state.ρ
    q = PhasePartition(ts)
    T::FT = air_temperature(ts)
    I_d::FT = internal_energy_dry(ts)
    I_v::FT = internal_energy_vapor(ts)
    I_l::FT = internal_energy_liquid(ts)
    I_i::FT = internal_energy_ice(ts)
    _T_freeze::FT = T_freeze(param_set)
    _L_f::FT = latent_heat_fusion(ts)
    _cv_l::FT = cv_l(param_set)
    Φ::FT = gravitational_potential(atmos.orientation, aux)

    # temporary vars for summming different source terms
    S_qr::FT = FT(0)
    S_qs::FT = FT(0)
    S_ql::FT = FT(0)
    S_qi::FT = FT(0)
    S_qt::FT = FT(0)
    S_e::FT = FT(0)

    # source of rain via autoconversion
    tmp = conv_q_liq_to_q_rai(param_set, q.liq)
    S_qr += tmp
    S_ql -= tmp
    S_e -= tmp * (I_l + Φ)

    # source of snow via autoconversion
    tmp = conv_q_ice_to_q_sno(param_set, q, state.ρ, T)
    S_qs += tmp
    S_qi -= tmp
    S_e -= tmp * (I_i + Φ)

    # source of rain water via accretion cloud water - rain
    tmp = accretion(
        param_set,
        CM1M.LiquidType(),
        CM1M.RainType(),
        q.liq,
        q_rai,
        state.ρ,
    )
    S_qr += tmp
    S_ql -= tmp
    S_e -= tmp * (I_l + Φ)

    # source of snow via accretion cloud ice - snow
    tmp = accretion(
        param_set,
        CM1M.IceType(),
        CM1M.SnowType(),
        q.ice,
        q_sno,
        state.ρ,
    )
    S_qs += tmp
    S_qi -= tmp
    S_e -= tmp * (I_i + Φ)

    # sink of cloud water via accretion cloud water - snow
    tmp = accretion(
        param_set,
        CM1M.LiquidType(),
        CM1M.SnowType(),
        q.liq,
        q_sno,
        state.ρ,
    )
    if T < _T_freeze # cloud droplets freeze to become snow)
        S_qs += tmp
        S_ql -= tmp
        S_e -= tmp * (I_i + Φ)
    else # snow melts, both cloud water and snow become rain
        α::FT = _cv_l / _L_f * (T - _T_freeze)
        S_ql -= tmp
        S_qs -= tmp * α
        S_qr += tmp * (1 + α)
        S_e -= tmp * ((1 + α) * I_l - α * I_i + Φ)
    end

    # sink of cloud ice via accretion cloud ice - rain
    tmp1 = accretion(
        param_set,
        CM1M.IceType(),
        CM1M.RainType(),
        q.ice,
        q_rai,
        state.ρ,
    )
    # sink of rain via accretion cloud ice - rain
    tmp2 = accretion_rain_sink(param_set, q.ice, q_rai, state.ρ)
    S_qi -= tmp1
    S_e -= tmp1 * (I_i + Φ)
    S_qr -= tmp2
    S_e += tmp2 * _L_f
    S_qs += tmp1 + tmp2

    # accretion rain - snow
    if T < _T_freeze
        tmp = accretion_snow_rain(
            param_set,
            CM1M.SnowType(),
            CM1M.RainType(),
            q_sno,
            q_rai,
            state.ρ,
        )
        S_qs += tmp
        S_qr -= tmp
        S_e += tmp * _L_f
    else
        tmp = accretion_snow_rain(
            param_set,
            CM1M.RainType(),
            CM1M.SnowType(),
            q_rai,
            q_sno,
            state.ρ,
        )
        S_qs -= tmp
        S_qr += tmp
        S_e -= tmp * _L_f
    end

    # rain evaporation sink (it already has negative sign for evaporation)
    tmp = evaporation_sublimation(
        param_set,
        CM1M.RainType(),
        q,
        q_rai,
        state.ρ,
        T,
    )
    S_qr += tmp
    S_e -= tmp * (I_l + Φ)

    # snow sublimation/deposition source/sink
    tmp = evaporation_sublimation(
        param_set,
        CM1M.SnowType(),
        q,
        q_sno,
        state.ρ,
        T,
    )
    S_qs += tmp
    S_e -= tmp * (I_i + Φ)

    # snow melt
    tmp = snow_melt(param_set, q_sno, state.ρ, T)
    S_qs -= tmp
    S_qr += tmp
    S_e -= tmp * _L_f

    # total qt sink is the sum of precip sources
    S_qt = -S_qr - S_qs

    return (;
        S_ρ_qt = state.ρ * S_qt,
        S_ρ_ql = state.ρ * S_ql,
        S_ρ_qi = state.ρ * S_qi,
        S_ρ_qr = state.ρ * S_qr,
        S_ρ_qs = state.ρ * S_qs,
        S_ρ_e = state.ρ * S_e,
    )
end


text:
# shortnames.jl -- containing the code for the Shortname submodule of MathGL.jl
using MathGL
export MathGL

# Export all types with shorter names
typealias Graph MathGL.Graph
typealias AbstractData MathGL.AbstractData
typealias Graph MathGL.Graph
typealias Data MathGL.Data

export Graph, AbstractData, Graph, Data

# Export all macros
export @mglplot

# Export the remaining functions or at least many of them. Should be more selective.
adjust       =  MathGL.adjust
alpha        =  MathGL.alpha
ambient      =  MathGL.ambient
arc          =  MathGL.arc
arrowsize    =  MathGL.arrowsize
aspect       =  MathGL.aspect
axis         =  MathGL.axis
axisstl      =  MathGL.axisstl
background   =  MathGL.background
ball         =  MathGL.ball
barwidth     =  MathGL.barwidth
box          =  MathGL.box
calcscr      =  MathGL.calcscr
calcxyz      =  MathGL.calcxyz
circle       =  MathGL.circle
clearframe   =  MathGL.clearframe
clf          =  MathGL.clf
closegif     =  MathGL.closegif
colorid      =  MathGL.colorid
columnplot   =  MathGL.columnplot
combine      =  MathGL.combine
cone         =  MathGL.cone
contfx       =  MathGL.contfx
contfy       =  MathGL.contfy
contfz       =  MathGL.contfz
contx        =  MathGL.contx
conty        =  MathGL.conty
contz        =  MathGL.contz
crange       =  MathGL.crange
crust        =  MathGL.crust
ctick        =  MathGL.ctick
curve        =  MathGL.curve
cut          =  MathGL.cut
defaultfont  =  MathGL.defaultfont
delframe     =  MathGL.delframe
densx        =  MathGL.densx
densy        =  MathGL.densy
densz        =  MathGL.densz
dew          =  MathGL.dew
diffuse      =  MathGL.diffuse
dots         =  MathGL.dots
drop         =  MathGL.drop
ellipse      =  MathGL.ellipse
endframe     =  MathGL.endframe
endgroup     =  MathGL.endgroup
errbox       =  MathGL.errbox
#=eval         =  MathGL.eval=#
exportmgld   =  MathGL.exportmgld
face         =  MathGL.face
facenum      =  MathGL.facenum
flow         =  MathGL.flow
fog          =  MathGL.fog
font         =  MathGL.font
fplot        =  MathGL.fplot
fsurf        =  MathGL.fsurf
getframe     =  MathGL.getframe
getheight    =  MathGL.getheight
getnumframe  =  MathGL.getnumframe
getrgb       =  MathGL.getrgb
getrgba      =  MathGL.getrgba
getwidth     =  MathGL.getwidth
grad         =  MathGL.grad
gridplot     =  MathGL.gridplot
height       =  MathGL.height
importmgld   =  MathGL.importmgld
inplot       =  MathGL.inplot
light        =  MathGL.light
line         =  MathGL.line
loadfont     =  MathGL.loadfont
mark         =  MathGL.mark
marksize     =  MathGL.marksize
mask         =  MathGL.mask
meshnum      =  MathGL.meshnum
mgl          =  MathGL.mgl
mpirecv      =  MathGL.mpirecv
mpisend      =  MathGL.mpisend
multiplot    =  MathGL.multiplot
newframe     =  MathGL.newframe
numthreads   =  MathGL.numthreads
origin       =  MathGL.origin
origintick   =  MathGL.origintick
palette      =  MathGL.palette
perspective  =  MathGL.perspective
pipe         =  MathGL.pipe
plotfactor   =  MathGL.plotfactor
plot         =  MathGL.plot
plotid       =  MathGL.plotid
polygon      =  MathGL.polygon
quadplot     =  MathGL.quadplot
quality      =  MathGL.quality
ranges       =  MathGL.ranges
rasterize    =  MathGL.rasterize
#=reset        =  MathGL.reset=#
resetframes  =  MathGL.resetframes
restorefont  =  MathGL.restorefont
rhomb        =  MathGL.rhomb
rotate       =  MathGL.rotate
rotatetext   =  MathGL.rotatetext
scheme       =  MathGL.scheme
setnumframe  =  MathGL.setnumframe
setsize      =  MathGL.setsize
showframe    =  MathGL.showframe
sphere       =  MathGL.sphere
startgif     =  MathGL.startgif
startgroup   =  MathGL.startgroup
stickplot    =  MathGL.stickplot
stop         =  MathGL.stop
subplot      =  MathGL.subplot
surf         =  MathGL.surf
ternary      =  MathGL.ternary
text         =  MathGL.text
ticklen      =  MathGL.ticklen
tickrotate   =  MathGL.tickrotate
tickshift    =  MathGL.tickshift
tickskip     =  MathGL.tickskip
title        =  MathGL.title
transptype   =  MathGL.transptype
traj         =  MathGL.traj
tricont      =  MathGL.tricont
triplot      =  MathGL.triplot
tuneticks    =  MathGL.tuneticks
vect         =  MathGL.vect
vect3        =  MathGL.vect3
view         =  MathGL.view
width        =  MathGL.width
writebmp     =  MathGL.writebmp
writebps     =  MathGL.writebps
writeeps     =  MathGL.writeeps
#=write        =  MathGL.write=#
writegif     =  MathGL.writegif
writejpg     =  MathGL.writejpg
writejson    =  MathGL.writejson
writeobj     =  MathGL.writeobj
writeoff     =  MathGL.writeoff
writepng     =  MathGL.writepng
writeprc     =  MathGL.writeprc
writestl     =  MathGL.writestl
writesvg     =  MathGL.writesvg
writetex     =  MathGL.writetex
writetga     =  MathGL.writetga
writexyz     =  MathGL.writexyz
xrange       =  MathGL.xrange
xrotate      =  MathGL.xrotate
xtick        =  MathGL.xtick
yrange       =  MathGL.yrange
yrotate      =  MathGL.yrotate
ytick        =  MathGL.ytick
zoom         =  MathGL.zoom
zoomaxis     =  MathGL.zoomaxis
zrange       =  MathGL.zrange
zrotate      =  MathGL.zrotate
ztick        =  MathGL.ztick

export adjust
export alpha
export ambient
export arc
export arrowsize
export aspect
export axis
export axisstl
export background
export ball
export barwidth
export box
export calcscr
export calcxyz
export circle
export clearframe
export clf
export closegif
export colorid
export columnplot
export combine
export cone
export contx
export conty
export contz
export contfx
export contfy
export contfz
export crange
export crust
export ctick
export curve
export cut
export defaultfont
export delframe
export densx
export densy
export densz
export dew
export diffuse
export dots
export drop
export ellipse
export endframe
export endgroup
export errbox
#=export eval=#
export exportmgld
export face
export facenum
export flow
export fog
export font
export fplot
export fsurf
export getframe
export getheight
export getnumframe
export getrgb
export getrgba
export getwidth
export grad
export gridplot
export height
export importmgld
export inplot
export light
export line
export loadfont
export mark
export marksize
export mask
export meshnum
export mgl
export mpirecv
export mpisend
export multiplot
export newframe
export numthr
export origin
export origintick
export palette
export perspective
export pipe
export plotfactor
export plot
export plotid
export polygon
export quadplot
export quality
export ranges
export rasterize
export reset
export resetframes
export restorefont
export rhomb
export rotate
export rotatetext
export scheme
export setnumframe
export setsize
export showframe
export sphere
export startgif
export startgroup
export stickplot
export stop
export subplot
export surf
export ternary
export text
export ticklen
export tickrotate
export tickshift
export tickskip
export title
export transptype
export traj
export tricont
export triplot
export tuneticks
export view
export vect
export vect3
export width
export writebmp
export writebps
export writeeps
#=export write=#
export writegif
export writejpg
export writejson
export writeobj
export writeoff
export writepng
export writeprc
export writestl
export writesvg
export writetex
export writetga
export writexyz
export xrange
export xrotate
export xtick
export yrange
export yrotate
export ytick
export zoom
export zoomaxis
export zrange
export zrotate
export ztick


text:
module KmerIndexes

export KmerIndex, count, locate

using Bio.Seq

"""
Dense k-mer index.
"""
immutable KmerIndex{K,T}
    data::Vector{Vector{T}}
end

function KmerIndex{K}(::Type{DNAKmer{K}}, seq::DNASequence, step::Integer=1)
    data = Vector{Vector{Int}}(4^K)
    for idx in 1:endof(data)
        data[idx] = Vector{Int}()
    end
    for (pos, kmer) in each(DNAKmer{K}, seq, step)
        idx = UInt64(kmer)
        push!(data[idx+1], pos)
    end
    return KmerIndex{K,Int}(data)
end

function Base.count{K}(kmer::DNAKmer{K}, index::KmerIndex{K})
    return length(locate(kmer, index))
end

function locate{K}(kmer::DNAKmer{K}, index::KmerIndex{K})
    idx = UInt64(kmer)
    return index.data[idx+1]
end

end # module


text:
function experiment_window()

    algorithms = [get_algorithm_desc("MOEA/D")]
    problems = [get_problem_desc("ZDT1")]

    function()
        CImGui.Begin("Experiment")

        if CImGui.Button("+ Algorithm")

        end

        if CImGui.Button("+ Problem")

        end

        CImGui.End()
    end
end