text:
"""
	NEOSServer(;email)

Construct a `NEOSServer` object. The `email` argument should take the users valid
email address (required for solvers like CPLEX).
"""
mutable struct NEOSServer
	useragent::String
	host::String
	contenttype::String
	email::String
	NEOSServer(;email="") = new("NEOS.jl", "https://neos-server.org:3333", "text/xml", email)
end

function addemail!(s::NEOSServer, email::String)
	s.email = email
end

type NEOSJob
	number::Int
	password::String
end

_add_text(value, arg::String) = add_text(value, arg)
_add_text(value, arg) = add_text(new_child(value, "int"), string(arg))

function buildxml(s::NEOSServer, name::String, args...)
	xml = XMLDocument()
	mthd = 	create_root(xml, "methodCall")
	mname = new_child(mthd, "methodName")
	add_text(mname, name)
	if length(args) > 0
		params = new_child(mthd, "params")
		for a in args
			param = new_child(params, "param")
			value = new_child(param, "value")
			_add_text(value, a)
		end
	end
	return string(xml)
end

function apimethod(s::NEOSServer, name::String, args...)
	xml = buildxml(s, name, args...)
	res = send(s, xml)
	parseresponse(res)
end
function send(s::NEOSServer, xml::String)
	hdrs = Dict{String, String}("user-agent" => s.useragent, "host" => s.host,
	"content-type" => s.contenttype, "content-length" => string(length(xml)))
	res = post(s.host; headers=hdrs, data=xml)
end

function getvalues!(values, node)
	types = ["int", "i4", "string", "double", "base64", "dateTime.iso8601"]
	for child in child_nodes(node)
		if name(child) in types
			push!(values, content(child))
		else
			getvalues!(values, child)
		end
	end
end

function parseresponse(res)
	if res.status != 200
		error("XML-RPC failed with code: $(res.status)")
	end
	parameters = Any[]
	xml = parse_string(convert(String, res.data))

	xroot = root(xml)
	getvalues!(parameters, xroot)
	return parameters
end

# ===========================================
#
# 	NEOS API Methods
#
function neosHelp(s::NEOSServer)
	return apimethod(s, "help")[1]
end

messages = [:emailHelp, :welcome, :version, :ping, :printQueue]
for m in messages
	m_str = string(m)
	@eval 	function ($m)(s::NEOSServer)
				return apimethod(s, $(m_str))[1]
			end
end

lists = [:listAllSolvers, :listCategories]
for m in lists
	m_str = string(m)
	@eval 	function ($m)(s::NEOSServer)
				return apimethod(s, $(m_str))
			end
end

function getSolverTemplate(s::NEOSServer, category::String, solvername::String, inputMethod::String)
	apimethod(s, "getSolverTemplate", category, solvername, inputMethod)[1]
end

function listSolversInCategory(s::NEOSServer, category::String)
	apimethod(s, "listSolversInCategory", category)
end

function submitJob(s::NEOSServer, xmlstring::String, print_level::Int = 1)
	res = apimethod(s, "submitJob", xmlstring)
	if print_level >= 1
		println("==================")
		println("NEOS Job submitted")
		println("number:\t$(res[1])")
		println("pwd:\t$(res[2])")
		println("==================")
	end
	return NEOSJob(parse(Int, res[1]), res[2])
end

jobapimethods = [:getJobStatus, :killJob]
for m in jobapimethods
	m_str = string(m)
	@eval 	function ($m)(s::NEOSServer, j::NEOSJob)
				return apimethod(s, $(m_str), j.number, j.password)[1]
			end
end

function getJobInfo(s::NEOSServer, j::NEOSJob)
	apimethod(s, "getJobInfo", j.number, j.password)
end

for s in ["", "NonBlocking"]
	_final_str = "getFinalResults$s"
	_final = Symbol(_final_str)
	_intermediate_str = "getIntermediateResults$s"
	_intermediate = Symbol(_intermediate_str)
	@eval function ($_final)(s::NEOSServer, j::NEOSJob)
		decode_to_string(apimethod(s, $(_final_str), j.number, j.password)[1])
	end

	@eval function ($(_intermediate))(s::NEOSServer, j::NEOSJob; offset=0)
		(results, offset) = apimethod(s, $(_intermediate_str), j.number, j.password, offset)
		decode_to_string(results), parse(Int, offset)
	end
end

function decode_to_string(s)
	String(decode(Base64, replace(s, "\n", "")))
end


text:
"""
		Implementation of Hyperspherical Variational Auto-Encoders

		Original paper: https://arxiv.org/abs/1804.00891

		SVAEbase(q,g,zdim,hue,μzfromhidden,κzfromhidden)

		q --- encoder - in this case encoder only encodes from input to a hidden
						layer which is then transformed into parameters for the latent
						layer by `μzfromhidden` and `κzfromhidden` functions
		g --- decoder
		zdim --- dimension of the latent space
		hue --- Hyperspherical Uniform Entropy that is part of the KL divergence but depends only on dimensionality so can be computed in constructor
		μzfromhidden --- function that transforms the hidden layer to μ parameter of the latent layer by normalization
		κzfromhidden --- transforms hidden layer to κ parameter of the latent layer using softplus since κ is a positive scalar
"""
struct SVAEbase <: SVAE
	q
	g
	zdim
	hue
	μzfromhidden
	κzfromhidden

	"""
	SVAEbase(q, g, hdim, zdim, T) Constructor of the S-VAE where `zdim > 3` and T determines the floating point type (default Float32)
	"""
end

SVAEbase(q, g, hdim::Int, zdim::Int, T = Float32) = SVAEbase(q, g, zdim, convert(T, huentropy(zdim)), Adapt.adapt(T, Chain(Dense(hdim, zdim), x -> normalizecolumns(x))), Adapt.adapt(T, Dense(hdim, 1, x -> σ.(x) .* 100)))
function SVAEbase(inputDim::Int, hiddenDim::Int, latentDim::Int, numLayers::Int, nonlinearity::String, layerType::String, T = Float32)
	encoder = Adapt.adapt(T, FluxExtensions.layerbuilder(inputDim, hiddenDim, hiddenDim, numLayers, nonlinearity, "", layerType))
    decoder = Adapt.adapt(T, FluxExtensions.layerbuilder(latentDim, hiddenDim, inputDim, numLayers + 1, nonlinearity, "linear", layerType))
	return SVAEbase(encoder, decoder, hiddenDim, latentDim, T)
end

Flux.@treelike(SVAEbase)

"""
	loss(m::SVAEbase, x)

	Loss function of the S-VAE combining reconstruction error and the KL divergence
"""
function loss(m::SVAEbase, x, β)
	(μz, κz) = zparams(m, x)
	z = samplez(m, μz, κz)
	xgivenz = m.g(z)
	return Flux.mse(x, xgivenz) + β * mean(kldiv(m, κz))
end

function wloss(m::SVAEbase, x, β, d)
	(μz, κz) = zparams(m, x)
	z = samplez(m, μz, κz)
	zp = samplehsuniform(size(z))
	Ω = d(z, zp)
	xgivenz = m.g(z)
	return Flux.mse(x, xgivenz) + β * Ω
end

function mem_wloss(svae::SVAEbase, mem::KNNmemory, x, y, β, d, α)
	(μz, κz) = zparams(svae, x)
	z = samplez(svae, μz, κz)
	xgivenz = svae.g(z)

	repetitions = 10
	(priorμ, priorκ) = zparams(svae, repeat(mem.M', 1, repetitions))
	priorlabels = repeat(mem.V, repetitions)
	priorsamples = samplez(svae, priorμ, priorκ)

	loss = 0

	if count(y .== 1) > 0
		anom_ids = findall(y .== 1)
		anom_ids = anom_ids[rand(1:length(anom_ids), length(y) * 2)]
		μzanom = μz[:, anom_ids]
		κzanom = κz[anom_ids]
		zanom = samplez(svae, μzanom, collect(κzanom'))

		norm_ids = findall(y .== 0)
		norm_ids = norm_ids[rand(1:length(norm_ids), length(y) * repetitions)]
		μznorm = μz[:, norm_ids]
		κznorm = κz[norm_ids]
		znorm = samplez(svae, μznorm, collect(κznorm'))

		Ωnorm = d(znorm, priorsamples[:, priorlabels .== 0])
		Ωanom = d(zanom, priorsamples[:, priorlabels .== 1])

		return Flux.mse(x, xgivenz) + β * (α .* Ωnorm .+ (1 - α) .* Ωanom)
	else
		Ωnorm = d(z, priorsamples[:, priorlabels .== 0])
		return Flux.mse(x, xgivenz) + β * Ωnorm
	end

	# if isnan(loss) | (loss == Inf) | (loss == -Inf)
	# 	println("Watch out, loss is $loss")
	# 	println("There were $(count(y .== 1)) anomalies.")
	# end

	return loss
end

ext: jl

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
reponame = "arnavs/InstantiationTest"
tag1 = "master"
tag2 = "v0.1.0"
sha = "181f673" # 7 chars
pdir = pwd()

# Clean environment.
isdir(".projects") == false || rm(".projects", recursive = true)

prefix, path = activate_github(reponame)
@test prefix == "master"
@test path == joinpath(pdir, ".projects", "InstantiationTest-master")
@test Base.active_project() == joinpath(path, "Project.toml")

@warn "Shouldn't take too long on new precompilations, etc."
prefix, path = activate_github(reponame, tag = tag1)
@test prefix == "master"
@test path == joinpath(pdir, ".projects", "InstantiationTest-master")
@test Base.active_project() == joinpath(path, "Project.toml")

prefix, path = activate_github(reponame, tag = tag2)
@test prefix == "v0.1.0"
@test path == joinpath(pdir, ".projects", "InstantiationTest-v0.1.0")
@test Base.active_project() == joinpath(path, "Project.toml")

InstantiateFromURL.copy_env(reponame, "v0.1.0", "v0.1.0-edit")
prefix, path = activate_github(reponame, tag = "v0.1.0-edit")
@test prefix == "v0.1.0-edit"
@test path == joinpath(pdir, ".projects", "InstantiationTest-v0.1.0-edit")
@test Base.active_project() == joinpath(path, "Project.toml")

prefix, path = activate_github(reponame, sha = sha)
@test prefix == sha
@test path == joinpath(pdir, ".projects", "InstantiationTest-$sha")
@test Base.active_project() == joinpath(path, "Project.toml")


text:
using DataDeps

register(DataDep(
  "indian-liver-patient",
  "http://archive.ics.uci.edu/ml/datasets/ILPD+%28Indian+Liver+Patient+Dataset%29",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/00225/Indian%20Liver%20Patient%20Dataset%20(ILPD).csv",
  "84feac16488de5cf89bd22bd802c77f25841fe93e9ddd32355683e94d46c3425",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=11,
        feature_indices=1:10,
        categoric_indices=2,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 6b473b2d-4326-46b4-af38-07b61de287fc
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
		Pkg.PackageSpec(name="ImageIO", version="0.5"),
		Pkg.PackageSpec(name="ImageShow", version="0.2"),
		Pkg.PackageSpec(name="FileIO", version="1.6"),
		Pkg.PackageSpec(name="PNGFiles", version="0.3.6"),
		Pkg.PackageSpec(name="Colors", version="0.12"),
		Pkg.PackageSpec(name="ColorVectorSpace", version="0.8"),
		Pkg.PackageSpec(name="PlutoUI", version="0.7"),
		Pkg.PackageSpec(name="HypertextLiteral", version="0.5")
	])

	using Colors, ColorVectorSpace, ImageShow, FileIO
	using PlutoUI
	using HypertextLiteral
	using LinearAlgebra
end

# ╔═╡ 2e8c4a48-d535-44ac-a1f1-4cb26c4aece6
filter!(LOAD_PATH) do path
	path != "@v#.#"
end;

# ╔═╡ 60532aa0-740c-11eb-0402-af8ff117f042
md"Show grid lines $(@bind show_grid CheckBox(default=true))"

# ╔═╡ 35904b8e-7a28-4dbc-bbf9-b45da448452c
let

	range = -2:.1:2
	md"""

	a $(@bind a Slider(range; default=1.0, show_value=true))

	b $(@bind b Slider(range; default=0.0, show_value=true))

	c $(@bind c Slider(range; default=0.0, show_value=true))

	d $(@bind d Slider(range; default=1.0, show_value=true))

	**Re-run this cell to reset to identity transformation**
	"""
end

# ╔═╡ f085296d-48b1-4db6-bb87-db863bb54049
A = [
	a b
	c d
	]

# ╔═╡ d1757b2c-7400-11eb-1406-d937294d5388
md"**_Det(A)_ = $a * $d - $c * $b =  $(det(A))**"

# ╔═╡ c536dafb-4206-4689-ad6d-6935385d8fdf
md"""
## Appendix
"""

# ╔═╡ fb509fb4-9608-421d-9c40-a4375f459b3f
det_A = det(A)

# ╔═╡ 0ec60f6e-7627-11eb-17e5-f965a81f3545
if det_A != 0
	invA = inv(A)
end

# ╔═╡ 40655bcc-6d1e-4d1e-9726-41eab98d8472
img_sources = [
	"https://user-images.githubusercontent.com/6933510/108605549-fb28e180-73b4-11eb-8520-7e29db0cc965.png" => "Corgis",
	"https://user-images.githubusercontent.com/6933510/108883855-39690f80-7606-11eb-8eb1-e595c6c8d829.png" => "Arrows",
	"https://images.squarespace-cdn.com/content/v1/5cb62a904d546e33119fa495/1589302981165-HHQ2A4JI07C43294HVPD/ke17ZwdGBToddI8pDm48kA7bHnZXCqgRu4g0_U7hbNpZw-zPPgdn4jUwVcJE1ZvWQUxwkmyExglNqGp0IvTJZamWLI2zvYWH8K3-s_4yszcp2ryTI0HqTOaaUohrI8PISCdr-3EAHMyS8K84wLA7X0UZoBreocI4zSJRMe1GOxcKMshLAGzx4R3EDFOm1kBS/fluffy+corgi?format=2500w" => "Alan"
]

# ╔═╡ c0c90fec-0e55-4be3-8ea2-88b8705ee258
md"""
### Choose an image:

$(@bind img_source Select(img_sources))
"""

# ╔═╡ 4fcb4ac1-1ad1-406e-8776-4675c0fdbb43
img_original = load(download(img_source));

# ╔═╡ 52a8009e-761c-11eb-2dc9-dbccdc5e7886
typeof(img_original)

# ╔═╡ 7d0096ad-d89a-4ade-9679-6ee95f7d2044
function trygetpixel(img::AbstractMatrix, x::Float64, y::Float64)
	rows, cols = size(img)

	"The linear map [-1,1] ↦ [0,1]"
	f = t -> (t - -1.0)/(1.0 - -1.0)

	i = floor(Int, rows * f(-y))
	j = floor(Int, cols * f(x * (rows / cols)))

	if 1 < i ≤ rows && 1 < j ≤ cols
		img[i,j]
	else
		zero(eltype(img))
	end
end

# ╔═╡ 83d45d42-7406-11eb-2a9c-e75efe62b12c
function with_gridlines(img::Array{<:Any,2}; n=8)

	sep_i = size(img, 1) ÷ n
	sep_j = size(img, 2) ÷ n

	result = copy(img)
	# stroke = zero(eltype(img))#RGBA(RGB(1,1,1), 0.75)

	stroke = RGBA(1, 1, 1, 0.75)

	result[1:sep_i:end, :] .= stroke
	result[:, 1:sep_j:end] .= stroke

	# a second time, to create a line 2 pixels wide
	result[2:sep_i:end, :] .= stroke
	result[:, 2:sep_j:end] .= stroke

	return result
end

# ╔═╡ 55898e88-36a0-4f49-897f-e0850bd2b0df
img = if show_grid
	with_gridlines(img_original)
else
	img_original
end;

# ╔═╡ 8e0505be-359b-4459-9de3-f87ec7b60c23
[
	if det_A == 0
		RGB(1.0, 1.0, 1.0)
	else
		in_x, in_y = invA*[out_x, out_y]
		trygetpixel(img, in_x, in_y)
	end

	for out_y in LinRange(2, -2, 300),
		out_x in LinRange(-2, 2, 300)
]

# ╔═╡ Cell order:
# ╠═6b473b2d-4326-46b4-af38-07b61de287fc
# ╟─2e8c4a48-d535-44ac-a1f1-4cb26c4aece6
# ╟─c0c90fec-0e55-4be3-8ea2-88b8705ee258
# ╟─60532aa0-740c-11eb-0402-af8ff117f042
# ╠═8e0505be-359b-4459-9de3-f87ec7b60c23
# ╠═35904b8e-7a28-4dbc-bbf9-b45da448452c
# ╟─f085296d-48b1-4db6-bb87-db863bb54049
# ╟─d1757b2c-7400-11eb-1406-d937294d5388
# ╟─c536dafb-4206-4689-ad6d-6935385d8fdf
# ╟─fb509fb4-9608-421d-9c40-a4375f459b3f
# ╟─0ec60f6e-7627-11eb-17e5-f965a81f3545
# ╟─40655bcc-6d1e-4d1e-9726-41eab98d8472
# ╠═4fcb4ac1-1ad1-406e-8776-4675c0fdbb43
# ╠═52a8009e-761c-11eb-2dc9-dbccdc5e7886
# ╠═55898e88-36a0-4f49-897f-e0850bd2b0df
# ╠═7d0096ad-d89a-4ade-9679-6ee95f7d2044
# ╠═83d45d42-7406-11eb-2a9c-e75efe62b12c


text:
export Avac, vacA, vacExpVal

function Avac(A::QuTerm,fac)
    ii = length(A.bares.v)
    while ii >= 1
        O = A.bares.v[ii]
        if O.t in (BosonCreate_,FermionCreate_,TLSCreate_)
            ii==length(A.bares.v) && return (A,fac)
            break
        elseif O.t in (BosonDestroy_,FermionDestroy_,TLSDestroy_)
            return (A,zero(fac))
        elseif O.t == TLSz_
            # vacuum is eigenvalue of σz
            fac = -fac
        elseif O.t in (TLSx_,TLSy_)
            break
        else
            error("should not be reached")
        end
        ii -= 1
    end
    Anew = QuTerm(A.nsuminds,A.δs,A.params,A.expvals,A.corrs,BaseOpProduct(A.bares.v[1:ii]))
    return (Anew, fac)
end

function vacA(A::QuTerm,fac)
    ii = 1
    while ii <= length(A.bares.v)
        O = A.bares.v[ii]
        if O.t in (BosonDestroy_,FermionDestroy_,TLSDestroy_)
            ii==1 && return (A,fac)
            break
        elseif O.t in (BosonCreate_,FermionCreate_,TLSCreate_)
            return (A,zero(fac))
        elseif O.t == TLSz_
            # vacuum is eigenvalue of σz
            fac = -fac
        elseif O.t in (TLSx_,TLSy_)
            break
        else
            error("should not be reached")
        end
        ii += 1
    end
    Anew = QuTerm(A.nsuminds,A.δs,A.params,A.expvals,A.corrs,BaseOpProduct(A.bares.v[ii:end]))
    return (Anew, fac)
end

Avac(A::QuExpr) = QuExpr(Avac(t,s) for (t,s) in normal_form(A).terms)
vacA(A::QuExpr) = QuExpr(vacA(t,s) for (t,s) in normal_form(A).terms)

"""
    Avac(A::QuExpr), vacA(A::QuExpr)

Simplify operator by assuming it is applied to the vacuum from the left or
right, respectively. To be precise, `Avac(A)` returns ``A'`` such that ``A'|0⟩ =
A|0⟩``, while `vacA(A)` does the same for ``⟨0|A``."""
Avac, vacA

function _TLS_to_pm_normal(A::QuExpr)
    An = QuExpr()
    terms_to_clean = collect(A.terms)
    while !isempty(terms_to_clean)
        (t,s) = pop!(terms_to_clean)
        v = t.bares.v
        iTLS = findfirst(A->A.t ∈ (TLSx_,TLSy_,TLSz_), v)
        if iTLS === nothing
            _add_with_normal_order!(An,t,s)
        else
            O = v[iTLS]
            if O.t == TLSx_
                # σˣ = σ⁺ + σ⁻
                for On in (TLSCreate(O.name,O.inds),TLSDestroy(O.name,O.inds))
                    vn = [v[1:iTLS-1]; On; v[iTLS+1:end]]
                    tn = QuTerm(t.nsuminds,t.δs,t.params,t.expvals,t.corrs,BaseOpProduct(vn))
                    push!(terms_to_clean,tn=>s)
                end
            elseif O.t == TLSy_
                # σʸ = -i σ⁺ + i σ⁻
                for (On,sn) in ((TLSCreate(O.name,O.inds),-1im),(TLSDestroy(O.name,O.inds),1im))
                    vn = [v[1:iTLS-1]; On; v[iTLS+1:end]]
                    tn = QuTerm(t.nsuminds,t.δs,t.params,t.expvals,t.corrs,BaseOpProduct(vn))
                    push!(terms_to_clean,tn=>s*sn)
                end
            elseif O.t == TLSz_
                # σᶻ = 2σ⁺σ⁻ - 1
                for (On,sn) in (([TLSCreate(O.name,O.inds),TLSDestroy(O.name,O.inds)],2),(BaseOperator[],-1))
                    vn = [v[1:iTLS-1]; On; v[iTLS+1:end]]
                    tn = QuTerm(t.nsuminds,t.δs,t.params,t.expvals,t.corrs,BaseOpProduct(vn))
                    push!(terms_to_clean,tn=>s*sn)
                end
            end
        end
    end
    An
end

"""
    vacExpVal(A::QuExpr,S::QuExpr=1)

Calculate the vacuum expectation value ``⟨0|S^\\dagger A S|0⟩``, i.e., the
expectation value ``⟨ψ|A|ψ⟩`` for the state defined by ``|ψ⟩= S|0⟩```.
"""
function vacExpVal(A::QuExpr,stateop::QuExpr=QuExpr(QuTerm()))
    # simplify down as much as possible by applying vacuum from left and right
    # convert TLSx/y/z operators to TLSCreate/TLSDestroy to ensure that no bare operators survive
    vsAsv = vacA(Avac(_TLS_to_pm_normal(stateop' * A * stateop)))
    _map_quexpr_ops(t -> (@assert isempty(t.bares); t), vsAsv)
end


text:
@testset "Clique cliquetree utilities" begin
    g = cycle_graph(4)
    add_edge!(g, 1, 3)
    clique = maximal_cliques(g)
    @test get_nv(clique) == nv(g)
    @test get_ne(clique) == ne(g)
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
function srk_optimize(alg,dx,pop_size,imin,imax,jmin,jmax;
                    NLoptRandSeed = 0,parameter_minmax=5,max_eval=Int(1e8),
                    x_L = -parameter_minmax*ones(M),
                    x_U =  parameter_minmax*ones(M),
                    rand_minmax = 3,
                    initCon = rand_minmax*(2rand(44)-1),
                    tol = 1e-2,ftol = 1e-15,tol2 = 1e-5,
                    counterSteps=Int(1e5),counterSteps2=Int(1e6),
                    initStepSize=[],gpuEnabled=true,ptx_str  = "integration.ptx",
                    cudaCores = 1664,initStepSize2=1e-6,outfile="",constrain_c = true)
  ##Parameters
  N = 28
  N2= 20
  N = EQ_C ? 30 : END_C ? 28 : 26
  N2= INEQ_C ? 20 : 16
  if EQ_C && INEQ_C
    error("Can only have EQ_C or INEQ_C!")
  end
  M = 44
  count = [0]
  timeNow = now()
  timeNow = Dates.format(timeNow, "y-m-d-HH-MM-SS")
  outfile = open(joinpath(Pkg.dir("SRKGenerator"),"output","Opti$timeNow.txt"), "w")

  ## Script Start
  coef_ans = Vector{Float32}(36)
  powz = Vector{Int8}(36)
  poww = Vector{Int8}(36)
  g_L = -tol*ones(N)
  g_U =  tol*ones(N)
  sizei = length(imin:dx:imax)
  sizej = length(jmin:dx:jmax)
  totArea = (imax-imin)*(jmax-jmin)/(sizei*sizej)
  if gpuEnabled
    iarr = convert(Vector{Float32},collect(imin:dx:imax))
    jarr = convert(Vector{Float32},collect(jmin:dx:jmax))
    numCards = length(cudaCores)
    totCores = sum(cudaCores)
    equalDiv = sizei*sizej÷totCores + 1
    portions = equalDiv*cudaCores
    if numCards > 1
      startIdx = [0 cumsum(portions)[1:end-1]]
    else
      startIdx = [0]
    end
    g_iarr = Vector{CUDArt.CudaArray{Float32,1}}(numCards)
    g_jarr = Vector{CUDArt.CudaArray{Float32,1}}(numCards)
    g_tmp = Vector{CUDArt.CudaArray{Int32,1}}(numCards)
    g_coefs = Vector{CUDArt.CudaArray{Float32,1}}(numCards)
    ans = Vector{Int32}(numCards)
    integrationFuncs = Vector{CUDArt.CuFunction}(numCards)
    CUDArt.init(devices(dev->true))
    devlist = devices(dev->true)
    streams = [(device(dev); Stream()) for dev in devlist]
    for i = 1:numCards
      device(i-1)
      g_iarr[i]  = CudaArray(iarr)
      g_jarr[i]  = CudaArray(jarr)
      g_tmp[i]   = CudaArray(Int32,cudaCores[i])
      md = CuModule(ptx_str,false)
      integrationFuncs[i] = CuFunction(md,"integration")
    end
  end

  eV = [1;1;1;1]

  eval_f = (x,grad) -> f_maker(x,ans,coef_ans,powz,poww,integrationFuncs,cudaCores,numCards,g_coefs,g_iarr,g_jarr,sizei,sizej,equalDiv,startIdx,g_tmp,totArea,counterSteps,counterSteps2,outfile,gpuEnabled,count)
  eval_g = (tmp,x,grad) -> g_maker(x,tmp,eV,counterSteps,counterSteps2,outfile,count)
  eval_g_ineq = (tmp,x,grad) -> g_ineq_maker(x,tmp,eV,counterSteps,counterSteps2,outfile,count)
  opt = Opt(alg,M) #:LD_SLSQP, :LN_COBYLA (semi), :GN_ISRES, :LN_AUGLAG_EQ support equality constraints
  lower_bounds!(opt,x_L)
  upper_bounds!(opt,x_U)
  max_objective!(opt, eval_f)
  maxeval!(opt::Opt, max_eval::Integer)
  equality_constraint!(opt,eval_g,tol*ones(N))
  if constrain_c
    inequality_constraint!(opt,eval_g_ineq,tol*ones(N2))
  end
  ftol_abs!(opt,ftol)
  xtol_abs!(opt,ftol)
  if initStepSize != []
    initial_step!(opt,initStepSize)
  end
  population!(opt,pop_size)
  NLopt.srand(NLoptRandSeed)

  setupString = """

  -----------------Setup------------------
  Options: alg=$alg,dx=$dx,max_eval=$max_eval,popSize=$pop_size
  imin=$imin,jmin=$jmin,imax=$imax,jmax=$jmax,
  parameter_minmax=$parameter_minmax,
  randSeed=$NLoptRandSeed, rand_minmax=$rand_minmax,
  END_C=$END_C, INEQ_C=$INEQ_C

  Initial Condition:

  $initCon
  """
  println(setupString)
  write(outfile,setupString)
  flush(outfile)

  maxf,maxx,ret = optimize(opt,initCon)

  println(ret)
  tmp2 = Vector{Float64}(N)
  eval_g(tmp2,maxx,0)
  maxErr = maximum(tmp2)
  mathx = translateToMathematica(maxx)
  println(mathx)
  flush(STDOUT)
  julString = printForJulia(maxx)
  resString = """

  -----------------Final Result------------------
  Options: alg=$alg,dx=$dx,max_eval=$max_eval,popSize=$pop_size
  imin=$imin,jmin=$jmin,imax=$imax,jmax=$jmax,
  parameter_minmax=$parameter_minmax,
  randSeed=$NLoptRandSeed,rand_minmax=$rand_minmax,
  END_C=$END_C, INEQ_C=$INEQ_C

  $(count[1]) steps

  Completed with maxf : $maxf

  Constraints at $maxErr

  $julString

  $mathx

  """
  println(resString)
  if gpuEnabled
    CUDArt.close(devices(dev->true))
  end
  write(outfile,resString)
  flush(outfile)

  return resString
end


text:
# This file includes modified source code from https://github.com/odow/SDDP.jl
# as at d90faae19c90f1fa03636ebe1cee92b083c355c2

#  Copyright 2017, Oscar Dowson and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################

"""simulate function"""
function simulate(hydromodel::HydroPowerModel, number_replications::Int = 1;
    sampling_scheme = SDDP.InSampleMonteCarlo(  max_depth = hydromodel.params["stages"],
                                                terminate_on_dummy_leaf = false),
    kwargs...)
    solution = Dict{Symbol, Any}()

    start_time = time()
    solution[:simulations] = SDDP.simulate( hydromodel.policygraph,
                                            number_replications,
                                            sampling_scheme = sampling_scheme,
                                            custom_recorders = Dict{Symbol, Function}(
                                                :powersystem => build_sol_powermodels,
                                                :reservoirs => build_sol_reservoirs,
                                                :objective => objective_value
                                            ))
    solution[:solve_time] = end_time = time() - start_time

    solution[:params] = hydromodel.params
    solution[:machine] = Dict(
        :cpu => Sys.cpu_info()[1].model,
        :memory => string(Sys.total_memory()/2^30, " Gb")
    )

    # add original data dict
    solution[:data] = hydromodel.alldata

    return solution

end

"""PowerModels solution build"""
function build_sol_powermodels(sp::JuMP.Model)
    solve_time = 0.0
    try
        solve_time = MOI.get(sp, MOI.SolveTime())
    catch
        solve_time = 0.0
    end
    status = JuMP.termination_status(sp)
    built_sol = PowerModels.build_solution(sp.ext[:pm],solve_time,
        solution_builder = get_solution)
end

"""Reservoir solution build"""
function build_sol_reservoirs(sp::JuMP.Model)
    store = Dict{Symbol, Any}(
            :reservoir => JuMP.value.(sp[:reservoir]),
            :outflow => JuMP.value.(sp[:outflow]),
            :spill => JuMP.value.(sp[:spill])
        )
    return store
end

""
function get_solution(pm::AbstractPowerModel, sol::Dict{String,<:Any})
    PowerModels.solution_opf!(pm, sol)
    add_kcl_deficit(sol, pm)
end

""
function add_kcl_deficit(sol, pm::AbstractPowerModel)
    PowerModels.add_setpoint!(sol, pm, "bus", "deficit", :deficit; status_name="bus_type", inactive_status_value = 4)
end


text:
module kmean
using Distributions
using Random
using Printf

function k_mean(X::AbstractMatrix{<:Real},
     k::Integer,
     w::Array{<:Real};
     maxiter::Integer=200, # maximun iteration
     tol::Real=1e-6 # tolearance of change at converage
     )
     d, n = size(X)
     (2 <= k < n) || throw(ArgumentError("k must be 2 <= k < n, k=$k given."))
     centers =  init_centers(X, k)
     # println(centers)
     centers_new = Array{Float64, 2}(undef, d, k)
     # println(centers)
     clustDist = zeros(n, 2)
     delta = 1.0
     iter_num = 1
     inertia = 0
     while delta > tol
         if iter_num > maxiter
             break
         end

         for i in range(1, length=n, step =1)
             distList = Float64[]
             for j in range(1, length=k, step=1)
                 append!(distList, distance(centers[:,j], X[:, i]))
             end
             # println(distList)
             dist, ind = findmin(distList)
             # print(ind)
             # * w[i]
             clustDist[i, :] = [dist, ind]
         end
         inertia = 0
         for j in range(1, length = k, step = 1)
              dInx = clustDist[:, 2].==j
              dist_t = clustDist[:, 1].==j
              inertia += mean(dist_t)
              # println(dInx)
              ptInCluster = X[:, dInx]
              # ptWeight = w[dInx]
              # ptWeight = ptWeight ./ sum(ptWeight)
              #* ptWeight
              centers_new[:,j] = mean(ptInCluster, dims=2)
              # println(centers_new[:,j])
         end
         # println(centers_new)
         # println(centers_new - centers)
         delta = sum((centers_new - centers).^2)
         centers = copy(centers_new)
         @printf("iter:\t%d\tdelta:\t%0.4f\n", iter_num, delta)
         # println("delta: $delta")
         iter_num += 1
     end
     # println(centers)
     cent = deepcopy(centers)
     inertia -= sum(var(cent, dims=2))/k
     # println(centers)
     return centers, inertia
end

function init_centers(X::AbstractMatrix{<:Real},
     k::Integer)
     d, n = size(X)

     # init_index = rand(1:n, k)
     # init_centers = X[:,init_index]
     init_index = Inf64[]
     index = rand(1:n)
     init_cluster = X[:,index]
     init_centers = Array{Float64, 2}(undef, d, k)
     distListTotal = Array{Float64, 1}(undef, n)
     init_centers[:,1] = init_cluster

     for i in 2:k
         sum_all = 0
         for j in 1:n
             distList = Float64[]
             for j in 1:i
                 append!(distList, distance(init_centers[:,i], X[:, j]))
             end
             dist, ind = findmin(distList)
             distListTotal[j] = dist
             sum_all += dist
         end
         sum_all *= rand()
         for j in 1:n
             di = distListTotal[k]
             sum_all -= di
             if sum_all > 0
                 continue
             end
             init_centers[:,i] = X[:, j]
             break
         end
     end
     return init_centers
end

function distance(v1, v2)
    sqrt(sum((v1 - v2).^2))
end

function k_mean_predict(X::AbstractMatrix{<:Real},
    centers::AbstractMatrix{<:Real})
    d, n = size(X)
    d2, k = size(centers)
    (d == d2) || throw(ArgumentError(
    "dims of data $d and center $d2 should be equal"))
    clustInd = Int[]
    for i in range(1, length=n, step =1)
        distList = Float64[]
        for j in range(1, length=k, step=1)
            append!(distList, distance(centers[:,j], X[:, i]))
        end
        # println(distList)
        dist, ind = findmin(distList)
        # print(ind)
        append!(clustInd, ind)
    end
    return clustInd
end

function K_mean(X::AbstractMatrix{<:Real},
     k::Integer,
     w::Array{<:Real};
     maxiter::Integer=100, # maximun iteration
     tol::Real=1e-4, # tolearance of change at converage
     random_state::Integer=10
     )
     d, n = size(X)
     centers_ls = Array{Float64, 3}(undef, random_state, d, k)
     inertia_ls = []
    for iter in 1:random_state
        centers, inertia = k_mean(X, k, w)
        if sum(isnan.(centers)) > 0
            continue
        end
        # println(centers_ls)
        # println(typeof(centers_ls))
        println(size(centers))
        centers_ls[iter,:,:] = centers
        append!(inertia_ls, inertia)
    end
    inertia_min, min_ind = findmin(inertia_ls)
    return centers_ls[min_ind,:,:]
end

end  # module kmean


text:
# This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license

@deprecate airy(z::Number) airyai(z)
@deprecate airyx(z::Number) airyaix(z)
@deprecate airyprime(z::Number) airyaiprime(z)

function _airy(k::Integer, z::Complex{Float64})
    Base.depwarn("`airy(k,x)` is deprecated, use `airyai(x)`, `airyaiprime(x)`, `airybi(x)` or `airybiprime(x)` instead.",:airy)
    id = Int32(k==1 || k==3)
    if k == 0 || k == 1
        return _airy(z, id, Int32(1))
    elseif k == 2 || k == 3
        return _biry(z, id, Int32(1))
    else
        throw(DomainError(k, "`k` must be between 0 and 3."))
    end
end
function _airyx(k::Integer, z::Complex{Float64})
    Base.depwarn("`airyx(k,x)` is deprecated, use `airyaix(x)`, `airyaiprimex(x)`, `airybix(x)` or `airybiprimex(x)` instead.",:airyx)
    id = Int32(k==1 || k==3)
    if k == 0 || k == 1
        return _airy(z, id, Int32(2))
    elseif k == 2 || k == 3
        return _biry(z, id, Int32(2))
    else
        throw(DomainError(k, "`k` must be between 0 and 3."))
    end
end

for afn in (:airy,:airyx)
    _afn = Symbol("_"*string(afn))
    suf  = string(afn)[5:end]
    @eval begin
        function $afn(k::Integer, z::Complex{Float64})
            afn = $(QuoteNode(afn))
            suf = $(QuoteNode(suf))
            Base.depwarn("`$afn(k,x)` is deprecated, use `airyai$suf(x)`, `airyaiprime$suf(x)`, `airybi$suf(x)` or `airybiprime$suf(x)` instead.",$(QuoteNode(afn)))
            $_afn(k,z)
        end

        $afn(k::Integer, z::Complex) = $afn(k, float(z))
        $afn(k::Integer, z::Complex{<:AbstractFloat}) = throw(MethodError($afn,(k,z)))
        $afn(k::Integer, z::Complex{Float32}) = Complex{Float32}($afn(k, Complex{Float64}(z)))
        $afn(k::Integer, x::Real) = $afn(k, float(x))
        $afn(k::Integer, x::AbstractFloat) = real($afn(k, complex(x)))
    end
end


text:
module LinearVortexSource

using LinearAlgebra
using StaticArrays
using Base.Iterators

using ..AeroMDAO: AbstractPanel, AbstractPanel2D, Panel2D, WakePanel2D, Point2D, collocation_point, p1, p2, transform_panel, affine_2D, panel_length, panel_angle, panel_tangent, panel_normal, panel_dist, rotation, inverse_rotation, midpair_map, panel_velocity, pressure_coefficient, wake_panel, wake_panels, panel_points

include("singularities.jl")

export constant_source_velocity, linear_source_velocity_a, linear_source_velocity_b, linear_vortex_velocity_a, linear_vortex_velocity_b

include("matrix.jl")

export total_velocity, source_velocity, vortex_velocity, vortex_influence_matrix, source_influence_matrix, neumann_boundary_condition, kutta, two_point_matrix, source_matrix, vortex_matrix, constant_source_matrix, constant_source_boundary_condition

end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
# add Merly    # pkg mode

using Merly

@page "/" HTTP.Response(200,"Hello World from Merly!")    # 1
@page "/hello/:user" HTTP.Response(200,string("<b>Hello ",request.params["user"],"!</b>")) # 2

@route POST "/post" HTTP.Response(200,"I did something!")  # 3

start(host = "127.0.0.1", port = 8086, verbose = true)    # 9


text:
using Documenter, DocumenterMarkdown, ObservableCollector

makedocs(sitename="ObservableCollector.jl",
    modules=[ObservableCollector],
    format = Markdown(),
    pages = [
        "Overview" => "index.md",
    ],
    clean=true)

if !haskey(ENV, "TRAVIS") || ENV["TRAVIS"]==false
    run(`mkdocs build`)
end

deploydocs(
    repo   = "github.com/skleinbo/ObservableCollector.jl.git",
    deps   = Deps.pip("mkdocs", "pygments", "python-markdown-math"),
    make   = () -> run(`mkdocs build`),
    target = "site"
)


text:
#=##############################################################################################
Filename: system.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: this file defines the `Parameters` type.
=###############################################################################################

"""
Objects inheriting from the `Parameters` type contain dot indexed members and are passed into action functions. A single `<: Parameters` struct is defined for the entire simulation and must contain members required by all Action functions. Action functions are overloaded to accept an `Aircraft` struct and a `step_range::AbstractArray` object and return default values.

For example:

struct MyParameters{ta, tb} <: Parameters
    a::ta
    b::tb
end

"""
abstract type Parameters end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
# Hamming Distance
# ================
#
# The Hamming distance algorithm.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioAlignments.jl/blob/master/LICENSE.md

function hamming_distance(::Type{T}, a, b) where T
    m = length(a)
    @assert m == length(b)
    anchors = [AlignmentAnchor(0, 0, 0, OP_START)]
    if m == 0
        # empty string
        return 0, anchors
    end
    was_match = a[1] == b[1]
    d = ifelse(was_match, T(0), T(1))
    for i in 2:m
        if a[i] == b[i]
            if !was_match
                push!(anchors, AlignmentAnchor(i - 1, i - 1, i - 1, OP_SEQ_MISMATCH))
            end
            was_match = true
        else
            if was_match
                push!(anchors, AlignmentAnchor(i - 1, i - 1, i - 1, OP_SEQ_MATCH))
            end
            d += T(1)
            was_match = false
        end
    end
    push!(anchors, AlignmentAnchor(m, m, m, was_match ? OP_SEQ_MATCH : OP_SEQ_MISMATCH))
    return d, anchors
end


text:
using CollegeStratBase, CollegeStratData
using Test

@testset "All" begin
    include("wage_regr_test.jl")
    include("datasettings_test.jl")
    include("helper_test.jl")
    include("data_files_test.jl")
    include("dataframe_xy_test.jl")
    include("load_data_test.jl");
    include("worker_moments_test.jl");
end

# -----------


text:
# rect_im.jl

using MIRT: rect_im, image_geom
using MIRTjim: jim

using Plots: plot
using Test: @test, @test_throws, @inferred

"""
    rect_im_show()
"""
function rect_im_show()
	ig = image_geom(nx=2^8, ny=2^8, fov=100)

	x1 = rect_im(ig, [[[0.5, 0, 3, 20]*ig.dx..., 0, 1]';], oversample=3)
	p1 = jim(x1)

	x2 = rect_im(ig, :default ; oversample=3, chat=true)
	p2 = jim(x2, title="default rects")

	x3 = rect_im(ig, :my_rect ; oversample=3, chat=true)
	p3 = jim(x3, title="my rect")

	x4 =rect_im(ig, :smiley ; oversample=3, chat=true)
	p4 = jim(x4, title="smiley")

	plot(p1, p2, p3, p4)
end

ig = image_geom(nx=2^8, dx=3)
@test_throws String rect_im(ig, :bad)
rect_im(ig, :smiley, how=:fast, replace=true)
rect_im(ig, how=:slow, replace=true)
rect_im(32, ny=30, dx=3, params=:default, how=:slow, return_params=true)
rect_im(32, 30)
@test_throws String rect_im(32, :default, how=:bad)

rect_im_show()


text:
"""
    simplex_basis(a, b, i, j)

Evaluate 2D orthonormal polynomial at simplex (a, b) of order (i, j)

Translated from Simplex2DP.m

"""
function simplex_basis(a::T, b::T, i, j) where {T<:Real}
    # x, a, b, n
    h1 = JacobiP(a, 0, 0, i)
    h2 = JacobiP(b, 2 * i + 1, 0, j)

    return sqrt(2.0) * h1 * h2 * (1 - b)^i
end

simplex_basis(a::AbstractVector{T}, b::AbstractVector{T}, i, j) where {T<:Real} =
    [simplex_basis(a[k], b[k], i, j) for k in eachindex(a)]

function ∂simplex_basis(a::T, b::T, id, jd) where {T<:Real}
    #fa = jacobi(a, id, 0, 0)
    #dfa = djacobi(a, id, 0, 0)
    #gb = jacobi(b, jd, 2*id+1, 0)
    #dgb = djacobi(b, jd, 2*id+1, 0)
    fa = JacobiP(a, 0, 0, id)
    dfa = ∂JacobiP(a, 0, 0, id)
    gb = JacobiP(b, 2 * id + 1, 0, jd)
    dgb = ∂JacobiP(b, 2 * id + 1, 0, jd)

    # r-derivative
    # d/dr = da/dr d/da + db/dr d/db = (2/(1-s)) d/da = (2/(1-b)) d/da
    dmodedr = dfa * gb
    if id > 0
        dmodedr *= (0.5 * (1.0 - b))^(id - 1)
    end

    # s-derivative
    # d/ds = ((1+a)/2)/((1-b)/2) d/da + d/db
    dmodeds = dfa * (gb * (0.5 * (1.0 + a)))
    if id > 0
        dmodeds *= (0.5 * (1.0 - b))^(id - 1)
    end

    tmp = dgb * (0.5 * (1.0 - b))^id
    if id > 0
        tmp -= 0.5 * id * gb * ((0.5 * (1.0 - b))^(id - 1))
    end
    dmodeds += fa * tmp

    # normalization
    dmodedr *= 2^(id + 0.5)
    dmodeds *= 2^(id + 0.5)

    return dmodedr, dmodeds
end

function ∂simplex_basis(a::AbstractVector{T}, b::AbstractVector{T}, id, jd) where {T<:Real}
    dmodedr = zero(a)
    dmodeds = zero(b)

    for i in eachindex(a)
        dmodedr[i], dmodeds[i] = ∂simplex_basis(a[i], b[i], id, jd)
    end

    return dmodedr, dmodeds
end


text:
module HMM

using Distributions

import Distributions.rand
import Distributions.fit

immutable HiddenMarkovModel{TP, K}
    theta::Vector{TP}
    A::Matrix{Float64}
    pi::Vector{Float64}

    function HiddenMarkovModel(theta, A, pi)
        @assert length(theta) == K
        @assert size(A) == (K, K)
        @assert length(pi) == K

        tol = 1e-6

        #@assert sum(pi) - 1.0 < tol, sum(pi)

        for k=1:K
            #@assert(sum(A[k, :]) - 1.0 < tol, string(k, " ", sum(A[k, :])))
        end

        return new(theta, A, pi)
    end
end

function HiddenMarkovModel{TP}(::Type{TP}, K::Int64)
    pi = rand(K)
    pi /= sum(pi)

    theta = [TP() for k=1:K]

    A = rand(K, K)
    for k=1:K
        A[k, :] /= sum(A[k, :])
    end

    return HiddenMarkovModel{TP, K}(theta, A, pi)
end

function HiddenMarkovModel{TP}(theta::Vector{TP})
    K = length(theta)
    pi = rand(K)
    pi /= sum(pi)

    A = rand(K, K)
    for k=1:K
        A[k, :] /= sum(A[k, :])
    end

    return HiddenMarkovModel{TP, K}(theta, A, pi)
end

function rand{TP, K}(hmm::HiddenMarkovModel{TP, K}, len::Int64)
    z = zeros(Int, len)

    _pi = prior_distribution(hmm)
    z[1] = rand(_pi)
    _A = transition_distributions(hmm)
    for i=2:len
        z[i] = rand(_A[z[i-1]])
    end

    x = Array(typeof(rand(hmm.theta[1])), len)
    for i=1:len
        x[i] = rand(hmm.theta[z[i]])
    end

    return (x, z)
end

function transition_distributions{TP, K}(hmm::HiddenMarkovModel{TP, K})
    A = Array(Categorical, K)
    for k=1:K
        A[k, :] = Categorical(vec(hmm.A[k, :]))
    end

    return A
end

function viterbi{TP, K, TO}(hmm::HiddenMarkovModel{TP, K}, x::Array{TO})
    len = length(x)
    z = zeros(Int64, len)

    backptr = zeros(Int64, len-1, K)
    pr = zeros(K)

    for k=1:K
        pr[k] = logpdf(hmm.theta[k], x[1]) + log(hmm.pi[k])
    end

    for i=2:len
        next_pr = zeros(K)
        for k=1:K
            maxk, maxp = -1, -Inf
            for t=1:K
                prob = pr[t] + log(hmm.A[t, k]) + logpdf(hmm.theta[k], x[i])
                if prob > maxp
                    maxk, maxp = t, prob
                end
            end
            backptr[i-1, k] = maxk
            next_pr[k] = maxp
        end
        pr = next_pr
    end

    z[len] = indmax(pr)
    for i=len-1:-1:1
        z[i] = backptr[i, z[i+1]]
    end

    return z
end

numstates{TP, K}(::HiddenMarkovModel{TP, K}) = K
numstates{TP, K}(::Type{HiddenMarkovModel{TP, K}}) = K
partype{TP, K}(::HiddenMarkovModel{TP, K}) = TP
partype{TP, K}(::Type{HiddenMarkovModel{TP, K}}) = TP

function fit{HMM <: HiddenMarkovModel}(hmm_type::Type{HMM}, x; conv_eps=1e-3, init_params=None, fit_param=fit, print_iteration=true)
    len = size(x, 1)
    K = numstates(hmm_type)
    TP = partype(hmm_type)

    hmm = if init_params == None
        HiddenMarkovModel(TP, K)
    else
        HiddenMarkovModel(init_params)
    end

    old_likelihood = -Inf
    while true
        alpha = zeros(len, K)
        beta = zeros(len, K)
        c = zeros(len)

        alpha[1, :] = hmm.pi .* [pdf(hmm.theta[k], slicedim(x, 1, 1))[1] for k=1:K]
        c[1] = sum(alpha[1, :])
        alpha[1, :] /= K

        for i=2:len
            for k=1:K
                alpha[i, k] = pdf(hmm.theta[k], slicedim(x, 1, i))[1]
                a = 0.
                for t=1:K
                    a += hmm.A[t, k] * alpha[i-1, t]
                end
                alpha[i, k] *= a
            end
            c[i] = sum(alpha[i, :])
            alpha[i, :] /= c[i]
        end

        beta[len, :] = 1.
        for i=len-1:-1:1
            for k=1:K
                b = 0.
                for t=1:K
                    b += beta[i+1, t] * pdf(hmm.theta[t], slicedim(x, 1, i+1))[1] * hmm.A[k, t]
                end
                beta[i, k] = b / c[i+1]
            end
        end

        likelihood = sum(log(c))
        if print_iteration == true; println("EM iteration. log-likelihood=", likelihood); end

        gamma = alpha .* beta
        ksi = zeros(len, K, K)

        for i=2:len
            for k=1:K
                for t=1:K
                    ksi[i, t, k] = alpha[i-1, t] * beta[i, k] * pdf(hmm.theta[k], slicedim(x, 1, i))[1] * hmm.A[t, k] / c[i]
                end
            end
        end

        A = zeros(K, K)
        for k=1:K
            for t=1:K
                for i=2:len
                    A[k, t] += ksi[i-1, k, t]
                end
            end
            A[k, :] /= sum(A[k, :])
        end

        theta = Array(TP, 0)
        for k=1:K
            push!(theta, fit_param(TP, x, vec(gamma[:, k])))
        end

        hmm = HiddenMarkovModel{TP, K}(theta, A, vec(gamma[1, :] / sum(gamma[1, :])))

#assert(old_likelihood < likelihood, "likelihood should monotonically increase")
        if abs(old_likelihood - likelihood) < conv_eps
            return hmm
        end

        old_likelihood = likelihood
    end
end

function fit(::Type{Normal}, x::Vector{Float64}, weights::Vector{Float64})
    mu = sum(x .* weights) / sum(weights)
    y = x - mu
    var = sum((y .^ 2) .* weights) / sum(weights)

    return Normal(mu, sqrt(var))
end

prior_distribution{TP, K}(hmm::HiddenMarkovModel{TP, K}) = Categorical(hmm.pi)

export HiddenMarkovModel, transition_distributions, prior_distribution, rand
export numstates, partype
export viterbi, fit

end

max_stars_repo_name: laisun/koan

max_issues_repo_name: laisun/koan

max_forks_repo_name: laisun/koan


text:
module Walker2014Sampling

import Base:
    rand

using Random:
    AbstractRNG

export
    rand,
    Walker2014Sampler

"""
    Walker2014Sampler(k::Int)

Walker's sampler with parameter `k`.
"""
struct Walker2014Sampler
    k::Int
    p::Vector{Float64}
    q::Vector{Float64} # internal parameter
    r::Vector{Float64} # internal parameter
end

function Walker2014Sampler(k::Int)
    Walker2014Sampler(k, zeros(2k - 1), zeros(k), zeros(2k - 1))
end

function get_q!(q, p)
    k = length(q)
    for j = 1:k
        q[j] = 0.0
        for l = j:(j + k - 1)
            q[j] += p[l]
        end
    end
    q .^= -1
    return nothing
end

function get_r!(r, q, p)
    k = length(q)
    for j = 1:(2k - 1)
        r[j] = 0.0
        for l in max(1, 1 + j - k):min(j, k)
            r[j] += q[l]
        end
    end
    @inbounds for j in 1:(2k - 1)
        r[j] = iszero(p[j]) ? 0.0 : p[j] * r[j]
    end
    sumr = sum(r)
    r ./= sumr
    return nothing
end

function rcat(rng, p)
    k = 0
    u = rand(rng)
    F = zero(eltype(p))
    while F < u
        k += 1
        F += p[k]
    end
    return k
end

"""
    rand(rng::AbstractRNG, p::Function, s::Walker2014Sampler, x0::Int)

Draw the next state in a MCMC with target unnormalized pdf `p` using the sampler `s`, assuming that the current state is `x0`. Both `rng` and `s` are modified in the process.
"""
function rand(rng::AbstractRNG, p::Function, s::Walker2014Sampler, x0::Int)
    lmin = max(1 + s.k - x0, 1)
    for l = 1:(lmin - 1)
        s.p[l] = 0.0
    end
    for l = lmin:(2 * s.k - 1)
        j = x0 - s.k + l
        s.p[l] = p(j)
    end
    get_q!(s.q, s.p)
    get_r!(s.r, s.q, s.p)
    return x0 += rcat(rng, s.r) - s.k
end

end # module


text:
module DataStructure101

    file_name = [
        "stack",
        "queue",
        "deque",
        "circular_queue",
        "sparse_array",
        "linked_list",
        "binary_tree",
        "union_find",
        "heap",
        "min_max_heap",
        "priority_queue",
        "graph",
    ]

    for f in file_name
        include("$f.jl")
    end

end  # module DataStructure101


text:
export instanceResolve

function instanceResolve(
        globalTC::GlobalTC,
        localImplicits::InstResolCtx,
        target::HMT,
        ln::LineNumberNode,
    )::IR.Expr
    tcstate = globalTC.tcstate
    target = tcstate.prune(target)

    candicates = Pair{InstRec, Vector{HMT}}[]
    for rec in localImplicits
        if !rec.isPruned
            rec.t = tcstate.prune(rec.t)
        end
        is_less, evidences = less_than_under_evidences(target, rec.t)
        if is_less
            push!(candicates, rec => evidences)
        end
    end
    !isempty(candicates) && begin
        if length(candicates) === 1
            rec, evidences = candicates[1]
            explicits = IR.Expr[instanceResolve(globalTC, localImplicits, evi, ln) for evi in evidences]
            return IR.applyExplicits(rec.gensym, explicits, target, ln)
        else
            throw(MLError(
                ln,
                DuplicateInstanceError([a.first for a in candicates])))
        end
    end

    tHead = getTConsHead(target)
    globalImplicits = get!(globalTC.globalImplicits, tHead) do
        InstRec[]
    end
    for i in eachindex(globalImplicits)
        rec = globalImplicits[i]

        if !rec.isPruned
            rec.t = tcstate.prune(rec.t)
        end
        is_less, evidences = less_than_under_evidences(target, rec.t)
        if is_less
            push!(candicates, rec => evidences)
        end
    end
    !isempty(candicates) && begin
        if length(candicates) === 1
            rec, evidences = candicates[1]
            explicits = IR.Expr[instanceResolve(globalTC, localImplicits, evi, ln) for evi in evidences]
            return IR.applyExplicits(rec.gensym, explicits, target, ln)
        else
            throw(MLError(
                ln,
                DuplicateInstanceError([a.first for a in candicates])))
        end
    end

    throw(MLError(
        ln,
        InstanceNotFound(target)))
end

function instanceResolveAndType(
    globalTC::GlobalTC,
    localImplicits::InstResolCtx,
    target::HMT,
    ln::LineNumberNode,
)
tcstate = globalTC.tcstate
target = tcstate.prune(target)

candicates = Pair{InstRec, Vector{HMT}}[]
for rec in localImplicits
    if !rec.isPruned
        rec.t = tcstate.prune(rec.t)
    end
    is_less, evidences = less_than_under_evidences(target, rec.t)
    if is_less
        push!(candicates, rec => evidences)
    end
end
!isempty(candicates) && begin
    if length(candicates) === 1
        rec, evidences = candicates[1]
        explicits = IR.Expr[instanceResolve(globalTC, localImplicits, evi, ln) for evi in evidences]
        return rec.t, IR.applyExplicits(rec.gensym, explicits, target, ln)
    else
        throw(MLError(
            ln,
            DuplicateInstanceError([a.first for a in candicates])))
    end
end

tHead = getTConsHead(target)
globalImplicits = get!(globalTC.globalImplicits, tHead) do
    InstRec[]
end
for i in eachindex(globalImplicits)
    rec = globalImplicits[i]

    if !rec.isPruned
        rec.t = tcstate.prune(rec.t)
    end
    is_less, evidences = less_than_under_evidences(target, rec.t)
    if is_less
        push!(candicates, rec => evidences)
    end
end

!isempty(candicates) && begin
    if length(candicates) === 1
        rec, evidences = candicates[1]
        explicits = IR.Expr[instanceResolve(globalTC, localImplicits, evi, ln) for evi in evidences]
        return rec.t, IR.applyExplicits(rec.gensym, explicits, target, ln)
    else
        throw(MLError(
            ln,
            DuplicateInstanceError([a.first for a in candicates])))
    end
end

throw(MLError(
    ln,
    InstanceNotFound(target)))
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

function randombubbles()
	store = []
	mycircle(r,n) = Lar.circle(r)([n])
	for k=1:50
		r = rand()
		n = 8
		while true
			n = abs(rand(Int8)+1)
			if n>2 break end
		end
		scale = Lar.s(0.25,0.25)
		transl = Lar.t(rand(2)...)
		s = Lar.Struct([ transl, scale, mycircle(r,n) ])
		push!(store, Lar.struct2lar(s))
	end
	s = Lar.Struct(store)
	V,EV = Lar.struct2lar(s)
	V = GL.normalize2(V)
	GL.VIEW([ GL.GLLines(V,EV) ])

	W,EW = Lar.fragmentlines((V,EV))
	V,FVs,EVs = Lar.arrange2D(W,EW)
end

# ////////////////////////////////////////////////////////////
# generation of 2D arrangement
V,FVs,EVs = randombubbles()

# native OpenGL visualization
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,1,1));
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,3,1));
GL.VIEW(GL.GLExplode(V,FVs,1.,1.,1.,99,1));
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,99,1));

GL.VIEW(GL.GLExplode(V,EVs,1.2,1.2,1.2,1,1));
GL.VIEW(GL.GLExplode(V,EVs,1,1,1,1,1));


text:
"""
App level functionality -- loading and managing app-wide components like configs, models, initializers, etc.
"""
module App

using Revise
using Genie

const ASSET_FINGERPRINT = ""

### PRIVATE ###

"""
    bootstrap(context::Union{Module,Nothing} = nothing) :: Nothing

Kickstarts the loading of a Genie app by loading the environment settings.
"""
function bootstrap(context::Union{Module,Nothing} = Genie.default_context(context)) :: Nothing
  if haskey(ENV, "GENIE_ENV") && isfile(joinpath(Genie.ENV_PATH, ENV["GENIE_ENV"] * ".jl"))
    isfile(joinpath(Genie.ENV_PATH, Genie.GLOBAL_ENV_FILE_NAME)) && Base.include(context, joinpath(Genie.ENV_PATH, Genie.GLOBAL_ENV_FILE_NAME))
    isfile(joinpath(Genie.ENV_PATH, ENV["GENIE_ENV"] * ".jl")) && Base.include(context, joinpath(Genie.ENV_PATH, ENV["GENIE_ENV"] * ".jl"))
  else
    ENV["GENIE_ENV"] = Genie.Configuration.DEV
    Core.eval(context, Meta.parse("const config = Configuration.Settings(app_env = Configuration.DEV)"))
  end

  for f in fieldnames(typeof(context.config))
    setfield!(Genie.config, f, getfield(context.config, f))
  end

  printstyled("""

   _____         _
  |   __|___ ___|_|___
  |  |  | -_|   | | -_|
  |_____|___|_|_|_|___|

  """, color = :red, bold = true)

  printstyled("| Web: https://genieframework.com\n", color = :light_black, bold = true)
  printstyled("| GitHub: https://github.com/genieframework/Genie.jl\n", color = :light_black, bold = true)
  printstyled("| Docs: https://genieframework.github.io/Genie.jl\n", color = :light_black, bold = true)
  printstyled("| Gitter: https://gitter.im/essenciary/Genie.jl\n", color = :light_black, bold = true)
  printstyled("| Twitter: https://twitter.com/GenieMVC\n\n", color = :light_black, bold = true)

  printstyled("Genie v$(Genie.Configuration.GENIE_VERSION)\n", color = :green, bold = true)
  printstyled("Active env: $(ENV["GENIE_ENV"] |> uppercase)\n\n", color = :light_blue, bold = true)

  nothing
end

end

### THIS IS LOADED INTO the Genie module !!!

using Logging
using .Configuration

### PUBLIC ###

"""
    newmodel(model_name::String; context::Union{Module,Nothing} = nothing) :: Nothing

Creates a new SearchLight `model` file.
"""
function newmodel(model_name::String; context::Union{Module,Nothing} = nothing) :: Nothing
  context = default_context(context)

  try
    Core.eval(context, :(SearchLight.Generator.newmodel($model_name)))

    load_resources()
  catch ex
    @error ex
  end

  nothing
end

"""
    newcontroller(controller_name::String) :: Nothing

Creates a new `controller` file. If `pluralize` is `false`, the name of the controller is not automatically pluralized.
"""
function newcontroller(controller_name::String; pluralize::Bool = true) :: Nothing
  Generator.newcontroller(Dict{String,Any}("controller:new" => controller_name), pluralize = pluralize)
  load_resources()

  nothing
end

"""
    newresource(resource_name::String; pluralize::Bool = true, context::Union{Module,Nothing} = nothing) :: Nothing

Creates all the files associated with a new resource.
If `pluralize` is `false`, the name of the resource is not automatically pluralized.
"""
function newresource(resource_name::String; pluralize::Bool = true, context::Union{Module,Nothing} = nothing) :: Nothing
  context = default_context(context)

  Generator.newresource(Dict{String,Any}("resource:new" => resource_name), pluralize = pluralize)

  try
    pluralize || error("SearchLight resources need to be pluralized")
    Core.eval(context, :(SearchLight.Generator.newresource(uppercasefirst($resource_name)))) # SearchLight resources don't work on singular
  catch ex
    @error ex
    @warn "Skipping SearchLight"
  end

  load_resources()

  nothing
end

"""
    newmigration(migration_name::String, context::Union{Module,Nothing} = nothing) :: Nothing

Creates a new SearchLight migration file.
"""
function newmigration(migration_name::String; context::Union{Module,Nothing} = nothing) :: Nothing
  context = default_context(context)

  try
    Core.eval(context, :(SearchLight.Generator.new_migration(Dict{String,Any}("migration:new" => $migration_name))))
  catch ex
    @error ex
  end

  nothing
end

"""
    newtablemigration(migration_name::String, context::Union{Module,Nothing} = nothing) :: Nothing

Creates a new migration prefilled with code for creating a new table.
"""
function newtablemigration(migration_name::String; context::Union{Module,Nothing} = nothing) :: Nothing
  context = default_context(context)

  try
    Core.eval(context, :(SearchLight.Generator.new_table_migration(Dict{String,Any}("migration:new" => $migration_name))))
  catch ex
    @error ex
  end

  nothing
end

"""
    newtask(task_name::String) :: Nothing

Creates a new Genie `Task` file.
"""
function newtask(task_name::String) :: Nothing
  endswith(task_name, "Task") || (task_name = task_name * "Task")
  Genie.Toolbox.new(Dict{String,Any}("task:new" => task_name), Genie.config)

  nothing
end

### PRIVATE ###

"""
    load_libs(root_dir::String = Genie.LIB_PATH) :: Nothing

Recursively adds subfolders of `lib/` to LOAD_PATH.
The `lib/` folder, if present, is designed to host user code in the form of .jl files.
This function loads user code into the Genie app.
"""
function load_libs(root_dir::String = Genie.LIB_PATH) :: Nothing
  isdir(root_dir) || return nothing

  push!(LOAD_PATH, root_dir)

  for (root, dirs, files) in walkdir(root_dir)
    for dir in dirs
      p = joinpath(root, dir)
      in(p, LOAD_PATH) || push!(LOAD_PATH, p)
    end
  end

  nothing
end

"""
    load_resources(root_dir::String = RESOURCES_PATH) :: Nothing

Recursively adds subfolders of `resources/` to LOAD_PATH.
"""
function load_resources(root_dir::String = Genie.RESOURCES_PATH) :: Nothing
  isdir(root_dir) || return nothing

  push!(LOAD_PATH, root_dir)

  for (root, dirs, files) in walkdir(root_dir)
    for dir in dirs
      p = joinpath(root, dir)
      in(p, LOAD_PATH) || push!(LOAD_PATH, joinpath(root, dir))
    end
  end

  nothing
end

"""
    load_helpers(root_dir::String = HELPERS_PATH) :: Nothing

Recursively adds subfolders of `helpers/` to LOAD_PATH.
"""
function load_helpers(root_dir::String = Genie.HELPERS_PATH) :: Nothing
  isdir(root_dir) || return nothing

  push!(LOAD_PATH, root_dir)

  for (root, dirs, files) in walkdir(root_dir)
    for dir in dirs
      p = joinpath(root, dir)
      in(p, LOAD_PATH) || push!(LOAD_PATH, joinpath(root, dir))
    end
  end

  nothing
end

"""
    load_configurations(root_dir::String = CONFIG_PATH; context::Union{Module,Nothing} = nothing) :: Nothing

Loads (includes) the framework's configuration files into the app's module `context`.
The files are set up with `Revise` to be automatically reloaded.
"""
function load_configurations(root_dir::String = Genie.CONFIG_PATH; context::Union{Module,Nothing} = nothing) :: Nothing
  context = default_context(context)

  secrets_path = joinpath(root_dir, Genie.SECRETS_FILE_NAME)
  isfile(secrets_path) && Revise.track(context, secrets_path, define = true)

  nothing
end

"""
    load_initializers(root_dir::String = CONFIG_PATH; context::Union{Module,Nothing} = nothing) :: Nothing

Loads (includes) the framework's initializers.
The files are set up with `Revise` to be automatically reloaded.
"""
function load_initializers(root_dir::String = Genie.CONFIG_PATH; context::Union{Module,Nothing} = nothing) :: Nothing
  context = default_context(context)

  dir = joinpath(root_dir, Genie.INITIALIZERS_FOLDER)

  isdir(dir) || return nothing

  f = readdir(dir)
  for i in f
    fi = joinpath(dir, i)
    endswith(fi, ".jl") && Revise.track(context, fi, define = true)
  end

  nothing
end

"""
    load_plugins(root_dir::String = PLUGINS_PATH; context::Union{Module,Nothing} = nothing) :: Nothing

Loads (includes) the framework's plugins initializers.
"""
function load_plugins(root_dir::String = Genie.PLUGINS_PATH; context::Union{Module,Nothing} = nothing) :: Nothing
  context = default_context(context)

  isdir(root_dir) || return nothing

  for i in readdir(root_dir)
    fi = joinpath(root_dir, i)
    endswith(fi, ".jl") && Revise.track(context, fi, define = true)
  end

  nothing
end

"""
    load_routes_definitions(routes_file::String = Genie.ROUTES_FILE_NAME; context::Union{Module,Nothing} = nothing) :: Nothing

Loads the routes file.
"""
function load_routes_definitions(routes_file::String = Genie.ROUTES_FILE_NAME; context::Union{Module,Nothing} = nothing) :: Nothing
  context = default_context(context)

  isfile(routes_file) && Revise.track(context, routes_file, define = true)

  nothing
end

"""
    secret_token(; context::Union{Module,Nothing} = nothing) :: String

Wrapper around /config/secrets.jl SECRET_TOKEN `const`.
Sets up the secret token used in the app for encryption and salting.
If there isn't a valid secrets file, a temporary secret token is generated for the current session only.
"""
function secret_token(; context::Union{Module,Nothing} = nothing) :: String
  context = default_context(context)

  if isdefined(context, :SECRET_TOKEN)
    context.SECRET_TOKEN
  else
    @warn "SECRET_TOKEN not configured - please make sure that you have a valid secrets.jl file.
          You can generate a new secrets.jl file with a random SECRET_TOKEN using Genie.Generator.write_secrets_file()
          or use the included /app/config/secrets.jl.example file as a model."

    st = Generator.secret_token()
    Core.eval(@__MODULE__, Meta.parse("""const SECRET_TOKEN = "$st" """))

    st
  end
end

"""
    default_context(context::Union{Module,Nothing})

Sets the module in which the code is loaded (the app's module)
"""
function default_context(context::Union{Module,Nothing})
  try
    context === nothing ? Main.UserApp : context
  catch ex
    @error ex
    @__MODULE__
  end
end

"""
    load(; context::Union{Module,Nothing} = nothing) :: Nothing

Main entry point to loading a Genie app.
"""
function load(; context::Union{Module,Nothing} = nothing) :: Nothing
  context = default_context(context)

  App.bootstrap(context)

  load_configurations(context = context)

  Core.eval(Genie, Meta.parse("""const SECRET_TOKEN = "$(secret_token(context = context))" """))
  Core.eval(Genie, Meta.parse("""const ASSET_FINGERPRINT = "$(App.ASSET_FINGERPRINT)" """))

  load_initializers(context = context)
  load_helpers()

  load_libs()
  load_resources()

  load_routes_definitions(context = context)

  load_plugins(context = context)

  nothing
end


text:
"""
    TriangularDist(a,b,c)

The *triangular distribution* with lower limit `a`, upper limit `b` and mode `c` has probability density function

```math
f(x; a, b, c)= \\begin{cases}
        0 & \\mathrm{for\\ } x < a, \\\\
        \\frac{2(x-a)}{(b-a)(c-a)} & \\mathrm{for\\ } a \\le x \\leq c, \\\\[4pt]
        \\frac{2(b-x)}{(b-a)(b-c)} & \\mathrm{for\\ } c < x \\le b, \\\\[4pt]
        0 & \\mathrm{for\\ } b < x,
        \\end{cases}
```

```julia
TriangularDist(a, b)        # Triangular distribution with lower limit a, upper limit b, and mode (a+b)/2
TriangularDist(a, b, c)     # Triangular distribution with lower limit a, upper limit b, and mode c

params(d)       # Get the parameters, i.e. (a, b, c)
minimum(d)      # Get the lower bound, i.e. a
maximum(d)      # Get the upper bound, i.e. b
mode(d)         # Get the mode, i.e. c
```

External links

* [Triangular distribution on Wikipedia](http://en.wikipedia.org/wiki/Triangular_distribution)

"""
struct TriangularDist{T<:Real} <: ContinuousUnivariateDistribution
    a::T
    b::T
    c::T
    TriangularDist{T}(a::T, b::T, c::T) where {T <: Real} = new{T}(a, b, c)
end

function TriangularDist(a::T, b::T, c::T; check_args=true) where {T <: Real}
    check_args && @check_args(TriangularDist, a <= c <= b)
    return TriangularDist{T}(a, b, c)
end

TriangularDist(a::T, b::T) where {T <: Real} = TriangularDist(a, b, middle(a, b))

TriangularDist(a::Real, b::Real, c::Real) = TriangularDist(promote(a, b, c)...)
TriangularDist(a::Integer, b::Integer, c::Integer) = TriangularDist(float(a), float(b), float(c))
TriangularDist(a::Real, b::Real) = TriangularDist(promote(a, b)...)
TriangularDist(a::Integer, b::Integer) = TriangularDist(float(a), float(b))

@distr_support TriangularDist d.a d.b

#### Conversions
convert(::Type{TriangularDist{T}}, a::Real, b::Real, c::Real) where {T<:Real} = TriangularDist(T(a), T(b), T(c))
convert(::Type{TriangularDist{T}}, d::TriangularDist{S}) where {T<:Real, S<:Real} = TriangularDist(T(d.a), T(d.b), T(d.c), check_args=false)

#### Parameters

params(d::TriangularDist) = (d.a, d.b, d.c)
partype(::TriangularDist{T}) where {T<:Real} = T

#### Statistics

mode(d::TriangularDist) = d.c

mean(d::TriangularDist) = (d.a + d.b + d.c) / 3

function median(d::TriangularDist)
    (a, b, c) = params(d)
    m = middle(a, b)
    c >= m ? a + sqrt((b - a) * (c - a)/2) :
             b - sqrt((b - a) * (b - c)/2)
end

_pretvar(a::Real, b::Real, c::Real) = a*a + b*b + c*c - a*b - a*c - b*c

function var(d::TriangularDist)
    (a, b, c) = params(d)
    _pretvar(a, b, c) / 18
end

function skewness(d::TriangularDist{T}) where T<:Real
    (a, b, c) = params(d)
    sqrt2 * (a + b - 2c) * (2a - b - c) * (a - 2b + c) / ( 5 * _pretvar(a, b, c)^(T(3)/2) )
end

kurtosis(d::TriangularDist{T}) where {T<:Real} = T(-3)/5

entropy(d::TriangularDist{T}) where {T<:Real} = one(T)/2 + log((d.b - d.a) / 2)

#### Evaluation

function pdf(d::TriangularDist{T}, x::Real) where T<:Real
    (a, b, c) = params(d)
    x <= a ? zero(T) :
    x <  c ? 2 * (x - a) / ((b - a) * (c - a)) :
    x == c ? 2 / (b - a) :
    x <= b ? 2 * (b - x) / ((b - a) * (b - c)) : zero(T)
end

function cdf(d::TriangularDist{T}, x::Real) where T<:Real
    (a, b, c) = params(d)
    x <= a ? zero(T) :
    x <  c ? (x - a)^2 / ((b - a) * (c - a)) :
    x == c ? (c - a) / (b - a) :
    x <= b ? 1 - (b - x)^2 / ((b - a) * (b - c)) : one(T)
end

function quantile(d::TriangularDist, p::Real)
    (a, b, c) = params(d)
    c_m_a = c - a
    b_m_a = b - a
    rl = c_m_a / b_m_a
    p <= rl ? a + sqrt(b_m_a * c_m_a * p) :
              b - sqrt(b_m_a * (b - c) * (1 - p))
end

function mgf(d::TriangularDist{T}, t::Real) where T<:Real
    if t == zero(t)
        return one(T)
    else
        (a, b, c) = params(d)
        u = (b - c) * exp(a * t) - (b - a) * exp(c * t) + (c - a) * exp(b * t)
        v = (b - a) * (c - a) * (b - c) * t^2
        return 2u / v
    end
end

function cf(d::TriangularDist{T}, t::Real) where T<:Real
    # Is this correct?
    if t == zero(t)
        return one(Complex{T})
    else
        (a, b, c) = params(d)
        u = (b - c) * cis(a * t) - (b - a) * cis(c * t) + (c - a) * cis(b * t)
        v = (b - a) * (c - a) * (b - c) * t^2
        return -2u / v
    end
end

#### Sampling

function rand(rng::AbstractRNG, d::TriangularDist)
    (a, b, c) = params(d)
    b_m_a = b - a
    u = rand(rng)
    b_m_a * u < (c - a) ? d.a + sqrt(u * b_m_a * (c - a)) :
                          d.b - sqrt((1 - u) * b_m_a * (b - c))
end


text:
"""
    refine_triangulation(points, simplex_inds, maxsize, k, factor)

Refine a triangulation recursively. Continue splitting some fraction of the simplices
with size reduction factor of 2 until until the maximum simplex size is below some
reference size.

Returns a tuple containing the refined triangulation and some information about it, as
follows: (points, simplex_inds, centroids, radii)
"""
function refine_recursive(points, simplex_inds, maxsize, k; niter = 1)
    if niter == 1
        println(size(points, 2), "D triangulation containing ", size(simplex_inds, 1), " simplices.")
        print("Target simplex radius = ", round(maxsize, 4), ". ")
        println("Split simplices into subsimplices with radii reduced by  factor of ", k, ".")
    end

    centroids, radii = centroids_radii2(points, simplex_inds)

    # If convergence is reached, return the triangulation and the centroids and radii of
    # the simplices furnishing the triangulation
    maxradius = maximum(radii)
    print("Refinement #", niter, "\tMaximum simplex radius = ", round(maxradius, 4),  ". ")

    if maxradius < maxsize || niter > 100
        if niter > 100
            println("\n\t\tMaximum number of iterations (100) reached. Ending refinement.")
        else
            println("\n\t\tRefinement finished.\n")
        end
        simplexvolumes = simplex_volumes(points, simplex_inds)

        return points, simplex_inds, centroids, radii, simplexvolumes
    end

    # If no convergence, continue splitting some fraction of the largest simplices.
    if length(radii) == 1
        percentile = 0.0
        split_indices = [1]
    else
        # What portion of the largest simplices to split (1 - percentile)?
        percentile = 0.999
        split_indices = []

        # If the simplices in the triangulation are very regular, there might not be
        # any simplices larger than a given percentile. If so, reduce percentile iteratively
        # until there are simplices larger than the percentile.
        while length(split_indices) == 0
            # Check if all elements are approximately the same size. If so, we must split all
            # of them.
            if all(y-> isapprox(radii[1], y), radii)
                split_indices =  find(radii .> 0)
            else
                split_indices = find(radii .> quantile(radii, percentile))
            end
            percentile = percentile - 0.02
        end
    end

    # The number of simplices to split
    n_split_simplices = length(split_indices)

    # The dimension of the space
    E = size(points, 2)

    # Generic rules for splitting a simplex
    splitting_rules = simplicial_subdivision(k, E)

    # Rules for forming the strictly new vertices of the subtriangulation
    rules = splitting_rules[1]

    # Array where each row represents one of the new simplices in the splitted simplex
    subtriangulation = splitting_rules[2]

    # How many new vertices are created each split?
    n_newvertices_eachsplit = size(rules, 1)

    println(" Splitting ", n_split_simplices, "/", size(simplex_inds, 1), " (", @sprintf("%.4f", 1.0 - percentile) ," %)", " simplices into ", n_split_simplices * n_newvertices_eachsplit - 1, " subsimplices.")

    # We need an array that can accomodate all of them. Each row in this
    # array will be a new vertex. Stacks of n_newvertices_eachsplit * E arrays.
    # We have as many as we have simplices to split.
    new_vertices = zeros(n_newvertices_eachsplit * n_split_simplices, E)

    # Fill the array by looping over all simplices that we want to split
    for i = 1:n_split_simplices
        # Figure out what the row indices corresponding to the ith simplex
        # must be. Marks the beginning of each of the simplex stacks in new_vertices
        ind = n_newvertices_eachsplit * (i -1)

        # Index of the simplex we need to split
        simplex_idx = split_indices[i]

        # Get the vertices of the simplex currently being splitted. Each of the
        # n_newvertices_eachsplit new vertices will be a linear combination
        # of these vertices. Each row is a vertex.
        vertices = points[simplex_inds[simplex_idx, :], :]

        # Generate the strictly new vertices for each sub
        for j = 1:n_newvertices_eachsplit
            # The index for a particular new vertex of the ith new simplex
            ind_newvertex = ind + j

            # Compute the jth vertex of the subtriangulation of the ith simplex

            # Go the jth new subsimplex of the ith splitted simplex. The entries of this
            # vector
            jth_subsimplex = rules[j, :]

            # Pick the corresponding original vertices with indices contained in rules[j, :]
            original_vertices = vertices[rules[j, :], :]

            new_vertices[ind_newvertex, :] = sum(original_vertices, dims=1) ./ k

        end
    end

    # Find the unique new vertices
    new_vertices_noreps = unique(new_vertices, 1)
    Ind = Vector{Int}(size(new_vertices, 1))
    for i = 1:size(new_vertices_noreps, 1)
        for j = 1:size(new_vertices, 1)
            if new_vertices_noreps[i, :] == new_vertices[j, :]
                Ind[j] = i
            end
        end
    end

    # Combine old and newly introduced vertices
    num_vertices_beforesplit = size(points, 1)
    allpoints = vcat(points, new_vertices_noreps)

    # Update the Ind array, so that we start at the new vertices
    Ind = Ind + num_vertices_beforesplit
    num_simplices_each_split = size(subtriangulation, 1)

    # The subsimplices formed by the splitting. Each row contains E + 1 indices referencing
    # the vertices furnishing that particular subsimplex (now found in the updated
    # points array).
    newtriangulation = Array{Float64}(num_simplices_each_split * n_split_simplices, E + 1)

    # For each simplex that we need to split,
    for i = 1:n_split_simplices
        # The beginning of the stack we need to fill.
        index = num_simplices_each_split * (i - 1)

        # Figure out what the row indices corresponding to the ith simplex
        # must be. Marks the beginning of each of the simplex stacks in new_vertices
        ind = n_newvertices_eachsplit * (i -1)

        # Index of the simplex we need to split
        simplex_idx = split_indices[i]

        # Pick the indices of the original vertices. Should be a column vector.
        inds_original_vertices = simplex_inds[simplex_idx, :]

        # Indices of the new vertices. Should be a column vector
        inds_new_vertices = Ind[(ind + 1):(ind + n_newvertices_eachsplit)]
        inds_all_vertices = vcat(inds_original_vertices, inds_new_vertices)

        # Go through each of the subsimplices formed by splitting that simplex
        for j = 1:num_simplices_each_split
            newtriangulation[index + j, :] = inds_all_vertices[subtriangulation[j, :]]
        end
    end

    # Indices of the simplices that are not split
    untouched_indices = setdiff(1:size(simplex_inds, 1), split_indices)

    all_simplex_inds = round.(Int, vcat(simplex_inds[untouched_indices, :],
                                  newtriangulation))

    refine_recursive(allpoints, all_simplex_inds, maxsize, k, niter = niter + 1)
end


text:
# This file is part of Intersections.jl
export svs

"""
    svs(postinglists, intersect2=baezayates)

Computes the intersection of the ordered lists in `postinglists` using a
small vs small strategy. Accepts an intersection algorithm of two sets.
"""
function svs(postinglists, intersect2=baezayates)
    sort!(postinglists, by=length, rev=true)
    curr = intersect2(pop!(postinglists), pop!(postinglists))

    if length(postinglists) > 0
        let prev = similar(curr)

            while length(postinglists) > 0
                empty!(prev)
                prev, curr = curr, intersect2(curr, pop!(postinglists), prev)
            end
        end
    end

    curr
end


text:
@testset "VoronoiTree" begin
    data = [randn(ComplexF64, 12) for i = 1:10_000]

    tree = HC.VoronoiTree(first(data), 1)
    for (i, d) in enumerate(data)
        insert!(tree, d, i)
    end

    @test length(tree) == length(data)
    @test sort!(collect(tree)) == 1:10_000
    @test all(i -> HC.search_in_radius(tree, data[i], 1e-12) == i, 1:length(data))
    let tree = HC.VoronoiTree(first(data), 1)
        @test all(enumerate(data)) do (i, d)
            HC.add!(tree, d, i, 1e-12) == (i, true) &&
                HC.search_in_radius(tree, d, 1e-12) == i
        end
    end

    data = [randn(ComplexF64, 12) for i = 1:10_000]
    tree = HC.VoronoiTree(first(data), 1)
    for (i, d) in enumerate(data)
        insert!(tree, d, i)
    end
    d = data[9342] .+ 1e-5
    @test HC.search_in_radius(tree, d, 1e-4) == 9342
    @test HC.search_in_radius(tree, d, 1e-6) === nothing

    # Test many points with nearly indentical distance to the inserted point
    p = shuffle!([[cis(k / 100 * 2π)] for k = 0:99])
    tree = HC.VoronoiTree(p[1], 1)
    for (i, pi) in enumerate(p)
        insert!(tree, pi, i)
    end
    @test isnothing(HC.search_in_radius(tree, [0.0], 1e-5))
    insert!(tree, [1e-5], 101)
    @test isnothing(HC.search_in_radius(tree, [0.0], 1e-6))
    @test HC.search_in_radius(tree, [0.0], 1e-4) == 101
end


text:
type PLRMat{T<:Number}
    # type for Partitioned Low Rank matrix
    blockType
    childPLRMat::Array{PLRMat{T},2}
    # if the matrix is compressible we use a low rank matrix
    Ulowrank::LowRankMatrix{T}
    Udense::Array{T,2}
    # size of the matrix
    n::Int64
    m::Int64
    kMax::Int64 # maximum rank
    epsilon::Float64

    # default constructor
    PLRMat() = new()
end

function dense2PLR{T<:Number}(A::Array{T,2}, kMax::Int64 ; epsilon=1e-6)
    # creating a new node in the tree
    node = PLRMat{T}()
    # computing the size of A
    (n,m) = size(A)
    # filling some fields
    node.n    = n;
    node.m    = m;
    node.kMax = kMax;
    node.epsilon = epsilon;
    if n <= kMax || m <= kMax
        node.Udense = A
        node.blockType = 'd'
    else
        # computing the randomized SVD
        (U,S,V) = randSVD(A, kMax+1)
        # testing the last eigenvalue
        if (S[end] < epsilon)
            # if it's true it means that the matrix is
            # compressible
            node.Ulowrank = LowRankMatrix(U[:,1:kMax], spdiagm(S[1:kMax])*V[:,1:kMax]',n,m,kMax)
            node.blockType = 'c'
        else
            # if the matrix is not compressible we need to build the
            # hierarchical structure
            node.blockType = 'h'
            # creating the array with the PlR matrices
            node.childPLRMat = Array(PLRMat{T}, 4,4)
            # size of the blocks of the partition
            nlocal = floor(Integer, n/4);
            mlocal = floor(Integer, m/4);
            # indices for the new partition (in 16 blocks)
            indi = [0, nlocal, 2*nlocal, 3*nlocal, n];
            indj = [0, mlocal, 2*mlocal, 3*mlocal, m];
            for ii = 1:4
                for jj = 1:4
                    # matrix will be divide in 16 equally spaced
                    indicesi = indi[ii]+1:indi[ii+1]
                    indicesj = indj[jj]+1:indj[jj+1]
                    node.childPLRMat[ii,jj] =  dense2PLR(A[indicesi,indicesj], kMax; epsilon=epsilon)
                end
            end
        end
    end

    return node
end

function matvecPLR{T}(H::PLRMat{T},b::Array{T,1})
    # check that they have the correct size
    @assert H.m == size(b)[1]
    # if the matrix is just a dense matrix
    if H.blockType == 'd'
        return H.Udense*b
    # if the matrix is a compressed matrix
    elseif H.blockType == 'c'
        return H.Ulowrank*b
    else
        # otherwise just call it in a recursive fashion
        # defining the index parition
        nlocal = floor(Integer, H.n/4);
        mlocal = floor(Integer, H.m/4);
        indj = [0, mlocal, 2*mlocal, 3*mlocal, H.m];
        indi = [0, nlocal, 2*nlocal, 3*nlocal, H.n];
        y = zeros(H.n)[:]
        for ii = 1:4
            for jj = 1:4
                # the indices for each one of the applications
                indicesi = indi[ii]+1:indi[ii+1]
                indicesj = indj[jj]+1:indj[jj+1]
                y[indicesi] += matvecPLR(H.childPLRMat[ii,jj], b[indicesj]);
            end
        end
        return y
    end
end

import Base.*
function *{T}(M::PLRMat{T}, v::Array{T,1})
    @assert M.m == length(v) # make sure vector is right size for matrix
    return matvecPLR(M,v)
end
;


text:
function solve(s::String)
    maximum(regex -> length(regex.match), eachmatch(r"[aeiou]+"i, s))
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
const CARRYCAP_PARAM      = Param(100000.0; bounds=(0.0, 10.0))
const INTRINSICRATE_PARAM = Param(0.1,      bounds=(0.0, 10.0))
const THRESHOLD_PARAM     = Param(0.5;      bounds=(0.0, 1.0))

"""
    GrowthRule <: CellRule

Abstract supertype for growth dynamics rules.
"""
abstract type GrowthRule{R,W} <: CellRule{R,W} end

"""
    ExponentialGrowth <: CellRule

    ExponentialGrowth(; rate, timestep, [nsteps_type])
    ExponentialGrowth{R}(; rate, timestep, [nsteps_type])
    ExponentialGrowth{R,W}(; rate, timestep, [nsteps_type])

Exponential growth of population size N based on an intrinsic growth rate ``r``, using the
exact solution between timesteps ``t`` and ``t-1``:

```math
N_t = N_{t-1}e^{r t}
```

# Keywords

- `rate`: Intrinsic growth rate. May be a `Number`, an `Aux` array or another `Grid`.
- `timestep`: Time step for the growth rate calculation, in a type compatible with the
    simulation `tspan`.
- `nsteps_type`: Specify the floating point type to use when `nsteps` is generated from the
    timestep, if it is required for type-stability or performance. The default is `Float64`.

Pass grid `Symbol`s to `R` or both `R` and `W` type parameters to use to specific grids.
"""
struct ExponentialGrowth{R,W,GR,TS,S} <: GrowthRule{R,W}
    rate::GR
    timestep::TS
    nsteps::S
end
function ExponentialGrowth{R,W}(;
    rate=INTRINSICRATE_PARAM,
    timestep,
    nsteps_type=Float64,
) where {R,W}
    ExponentialGrowth{R,W}(rate, timestep, zero(nsteps_type))
end

modifyrule(rule::ExponentialGrowth, data) = precalc_nsteps(rule, data)

@inline function applyrule(data, rule::ExponentialGrowth, N, I)
    N > zero(N) || return zero(N)
    rt = get(data, rule.rate, I...) * rule.nsteps

    return @fastmath N * exp(rt)
end

"""
    LogisticGrowth <: GrowthRule

    LogisticGrowth(; rate, carrycap, timestep, [nsteps_type])
    LogisticGrowth{R}(; rate, carrycap, timestep, [nsteps_type])
    LogisticGrowth{R,W}(; rate, carrycap, timestep, [nsteps_type])

Logistic growth rate of population size ``N`` based on an intrinsic growth rate ``r`` and
carry capacity ``K``, using the exact solution between timesteps ``t+1`` and ``t``:

```math
N_{t+1} = (N_t K) / (N_t + (K - N_t) e^{-rt})
```
Saturation only applies with positive growth.

# Keywords

These may be a `Number`, an `Aux` array or another `Grid`:

- `rate`: Intrinsic growth rate.
- `carrycap`: Carrying capacity.
- `timestep`: Time step for the growth rate, in a type compatible with the simulation `tspan`.
- `nsteps_type`: Specify the floating point type to use when `nsteps` is generated from the
    timestep, if it is required for type-stability or performance. The default is `Float64`.

Pass grid `Symbol`s to `R` or both `R` and `W` type parameters to use to specific grids.
"""
struct LogisticGrowth{R,W,GR,CC,TS,S} <: GrowthRule{R,W}
    rate::GR
    carrycap::CC
    timestep::TS
    nsteps::S
end
function LogisticGrowth{R,W}(;
    rate=INTRINSICRATE_PARAM,
    carrycap=CARRYCAP_PARAM,
    timestep,
    nsteps_type=Float64,
) where {R,W}
    LogisticGrowth{R,W}(rate, carrycap, timestep, zero(nsteps_type))
end

modifyrule(rule::LogisticGrowth, data) = precalc_nsteps(rule, data)

@inline function applyrule(data, rule::LogisticGrowth, N, I)
    N > zero(N) || return zero(N)

    rt = get(data, rule.rate, I...) * rule.nsteps
    k = get(data, rule.carrycap, I...)

    if rt > zero(rt)
        return @fastmath (N * k) / (N + (k - N) * exp(-rt))
    else
        return @fastmath N * exp(rt)
    end
end

"""
    ThresholdGrowth <: CellRule

    ThresholdGrowth(; rate, threshold)
    ThresholdGrowth{R}(; rate, threshold)
    ThresholdGrowth{R,W}(; rate, threshold)

Simple threshold mask. Values below a certain threshold are replaced with zero.

# Keywords

These may be a `Number`, an `Aux` array or another `Grid`.

- `rate`: Intrinsic growth rate.
- `threshold`: Minimum viability threshold below which population falls to zero.

Pass grid `Symbol`s to `R` or both `R` and `W` type parameters to use to specific grids.
"""
struct ThresholdGrowth{R,W,GR,Th} <: GrowthRule{R,W}
    rate::GR
    threshold::Th
end
function ThresholdGrowth{R,W}(;
    rate=INTRINSICRATE_PARAM,
    threshold=THRESHOLD_PARAM,
) where {R,W}
    ThresholdGrowth{R,W}(rate, threshold)
end

@inline function applyrule(data, rule::ThresholdGrowth, N, I)
    r = get(data, rule.rate, I...)
    t = get(data, rule.threshold, I...)

    return r >= t ? N : zero(N)
end

precalc_nsteps(rule, data) = precalc_nsteps(rule.timestep, rule, data)
precalc_nsteps(rulestep::DatePeriod, rule, data) =
    @set rule.nsteps = typeof(rule.nsteps)(currenttimestep(data) / Millisecond(rulestep))
precalc_nsteps(rulestep, rule, data) =
    @set rule.nsteps = typeof(rule.nsteps)(timestep(data) / rulestep)


text: function println() end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
@doc raw"""
    HyperbolicDiffusionEquations2D

The linear hyperbolic diffusion equations in two space dimensions.
A description of this system can be found in Sec. 2.5 of the book "I Do Like CFD, Too: Vol 1".
The book is freely available at http://www.cfdbooks.com/ and further analysis can be found in
the paper by Nishikawa [DOI: 10.1016/j.jcp.2007.07.029](https://doi.org/10.1016/j.jcp.2007.07.029)
"""
struct HyperbolicDiffusionEquations2D <: AbstractHyperbolicDiffusionEquations{2, 3}
  Lr::Float64
  Tr::Float64
  nu::Float64
  resid_tol::Float64
end

function HyperbolicDiffusionEquations2D()
  # diffusion coefficient
  nu = parameter("nu", 1.0)
  # relaxation length scale
  Lr = parameter("Lr", 1.0/(2.0*pi))
  # relaxation time
  Tr = Lr*Lr/nu
  # stopping tolerance for the pseudotime "steady-state"
  resid_tol = parameter("resid_tol", 1e-12)
  HyperbolicDiffusionEquations2D(Lr, Tr, nu, resid_tol)
end

get_name(::HyperbolicDiffusionEquations2D) = "HyperbolicDiffusionEquations2D"
varnames_cons(::HyperbolicDiffusionEquations2D) = @SVector ["phi", "q1", "q2"]
varnames_prim(::HyperbolicDiffusionEquations2D) = @SVector ["phi", "q1", "q2"]
default_analysis_quantities(::HyperbolicDiffusionEquations2D) = (:l2_error, :linf_error, :residual)

# Set initial conditions at physical location `x` for pseudo-time `t`
function initial_conditions_poisson_periodic(x, t, equation::HyperbolicDiffusionEquations2D)
  # elliptic equation: -νΔϕ = f
  # depending on initial constant state, c, for phi this converges to the solution ϕ + c
  if iszero(t)
    phi = 0.0
    q1  = 0.0
    q2  = 0.0
  else
    phi = sin(2.0*pi*x[1])*sin(2.0*pi*x[2])
    q1  = 2*pi*cos(2.0*pi*x[1])*sin(2.0*pi*x[2])
    q2  = 2*pi*sin(2.0*pi*x[1])*cos(2.0*pi*x[2])
  end
  return @SVector [phi, q1, q2]
end

function initial_conditions_poisson_nonperiodic(x, t, equation::HyperbolicDiffusionEquations2D)
  # elliptic equation: -νΔϕ = f
  if t == 0.0
    phi = 1.0
    q1  = 1.0
    q2  = 1.0
  else
    phi = 2.0*cos(pi*x[1])*sin(2.0*pi*x[2]) + 2.0 # ϕ
    q1  = -2.0*pi*sin(pi*x[1])*sin(2.0*pi*x[2])   # ϕ_x
    q2  = 4.0*pi*cos(pi*x[1])*cos(2.0*pi*x[2])    # ϕ_y
  end
  return @SVector [phi, q1, q2]
end

function initial_conditions_harmonic_nonperiodic(x, t, equation::HyperbolicDiffusionEquations2D)
  # elliptic equation: -νΔϕ = f
  if t == 0.0
    phi = 1.0
    q1  = 1.0
    q2  = 1.0
  else
    C   = 1.0/sinh(pi)
    phi = C*(sinh(pi*x[1])*sin(pi*x[2]) + sinh(pi*x[2])*sin(pi*x[1]))
    q1  = C*pi*(cosh(pi*x[1])*sin(pi*x[2]) + sinh(pi*x[2])*cos(pi*x[1]))
    q2  = C*pi*(sinh(pi*x[1])*cos(pi*x[2]) + cosh(pi*x[2])*sin(pi*x[1]))
  end
  return @SVector [phi, q1, q2]
end

function initial_conditions_jeans_instability(x, t, equation::HyperbolicDiffusionEquations2D)
  # gravity equation: -Δϕ = -4πGρ
  # Constants taken from the FLASH manual
  # https://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel.pdf
  rho0 = 1.5e7
  delta0 = 1e-3
  #
  phi = rho0*delta0 # constant background pertubation magnitude
  q1  = 0.0
  q2  = 0.0
  return @SVector [phi, q1, q2]
end

function initial_conditions_eoc_test_coupled_euler_gravity(x, t, equation::HyperbolicDiffusionEquations2D)

  # Determine phi_x, phi_y
  G = 1.0 # gravitational constant
  C = -2.0*G/pi
  A = 0.1 # perturbation coefficient must match Euler setup
  rho1 = A * sin(pi * (x[1] + x[2] - t))
  # intialize with ansatz of gravity potential
  phi = C * rho1
  q1  = C * A * pi * cos(pi*(x[1] + x[2] - t)) # = gravity acceleration in x-direction
  q2  = q1                                     # = gravity acceleration in y-direction

  return @SVector [phi, q1, q2]
end

function initial_conditions_sedov_self_gravity(x, t, equation::HyperbolicDiffusionEquations2D)
  # for now just use constant initial condition for sedov blast wave (can likely be improved)
  phi = 0.0
  q1  = 0.0
  q2  = 0.0
  return @SVector [phi, q1, q2]
end

# Apply boundary conditions
function boundary_conditions_harmonic_nonperiodic(u_inner, orientation, direction, x, t,
                                                  surface_flux_function,
                                                  equation::HyperbolicDiffusionEquations2D)
  # elliptic equation: -νΔϕ = f
  C   = 1.0/sinh(pi)
  phi = C*(sinh(pi*x[1])*sin(pi*x[2]) + sinh(pi*x[2])*sin(pi*x[1]))
  q1  = C*pi*(cosh(pi*x[1])*sin(pi*x[2]) + sinh(pi*x[2])*cos(pi*x[1]))
  q2  = C*pi*(sinh(pi*x[1])*cos(pi*x[2]) + cosh(pi*x[2])*sin(pi*x[1]))
  u_boundary = @SVector [phi, q1, q2]

  # Calculate boundary flux
  if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end

function boundary_conditions_poisson_nonperiodic(u_inner, orientation, direction, x, t,
                                                 surface_flux_function,
                                                 equation::HyperbolicDiffusionEquations2D)
  # elliptic equation: -νΔϕ = f
  phi = 2.0*cos(pi*x[1])*sin(2.0*pi*x[2]) + 2.0 # ϕ
  q1  = -2.0*pi*sin(pi*x[1])*sin(2.0*pi*x[2])   # ϕ_x
  q2  = 4.0*pi*cos(pi*x[1])*cos(2.0*pi*x[2])    # ϕ_y
  u_boundary = @SVector [phi, q1, q2]

  # Calculate boundary flux
  if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end

function boundary_conditions_sedov_self_gravity(u_inner, orientation, direction, x, t,
                                                surface_flux_function,
                                                equation::HyperbolicDiffusionEquations2D)
  u_boundary = initial_conditions_sedov_self_gravity(x, t, equation)

  # Calculate boundary flux
  if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end

# Apply source terms
function source_terms_poisson_periodic(ut, u, x, element_id, t, n_nodes, equation::HyperbolicDiffusionEquations2D)
  # elliptic equation: -νΔϕ = f
  # analytical solution: phi = sin(2πx)*sin(2πy) and f = -8νπ^2 sin(2πx)*sin(2πy)
  inv_Tr = inv(equation.Tr)
  C = -8.0*equation.nu*pi*pi

  for j in 1:n_nodes, i in 1:n_nodes
    x1 = x[1, i, j, element_id]
    x2 = x[2, i, j, element_id]
    tmp1 = sin(2.0*pi*x1)
    tmp2 = sin(2.0*pi*x2)
    ut[1, i, j, element_id] -= C*tmp1*tmp2
    ut[2, i, j, element_id] -= inv_Tr * u[2, i, j, element_id]
    ut[3, i, j, element_id] -= inv_Tr * u[3, i, j, element_id]
  end

  return nothing
end

function source_terms_poisson_nonperiodic(ut, u, x, element_id, t, n_nodes, equation::HyperbolicDiffusionEquations2D)
  # elliptic equation: -νΔϕ = f
  # analytical solution: ϕ = 2cos(πx)sin(2πy) + 2 and f = 10π^2cos(πx)sin(2πy)
  inv_Tr = inv(equation.Tr)

  for j in 1:n_nodes, i in 1:n_nodes
    x1 = x[1, i, j, element_id]
    x2 = x[2, i, j, element_id]
    ut[1, i, j, element_id] += 10 * pi^2 * cos(pi*x1) * sin(2.0*pi*x2)
    ut[2, i, j, element_id] -= inv_Tr * u[2, i, j, element_id]
    ut[3, i, j, element_id] -= inv_Tr * u[3, i, j, element_id]
  end

  return nothing
end

function source_terms_harmonic(ut, u, x, element_id, t, n_nodes, equation::HyperbolicDiffusionEquations2D)
  # harmonic solution ϕ = (sinh(πx)sin(πy) + sinh(πy)sin(πx))/sinh(π), so f = 0
  inv_Tr = inv(equation.Tr)

  for j in 1:n_nodes, i in 1:n_nodes
    ut[2, i, j, element_id] -= inv_Tr * u[2, i, j, element_id]
    ut[3, i, j, element_id] -= inv_Tr * u[3, i, j, element_id]
  end

  return nothing
end

# The coupled EOC test does not require additional sources
function source_terms_eoc_test_coupled_euler_gravity(ut, u, x, element_id, t, n_nodes, equation::HyperbolicDiffusionEquations2D)
  return source_terms_harmonic(ut, u, x, element_id, t, n_nodes, equation)
end

# Calculate 1D flux in for a single point
@inline function calcflux(u, orientation, equation::HyperbolicDiffusionEquations2D)
  phi, q1, q2 = u

  if orientation == 1
    f1 = -equation.nu*q1
    f2 = -phi/equation.Tr
    f3 = zero(phi)
  else
    f1 = -equation.nu*q2
    f2 = zero(phi)
    f3 = -phi/equation.Tr
  end

  return SVector(f1, f2, f3)
end

@inline function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::HyperbolicDiffusionEquations2D)
  # Obtain left and right fluxes
  f_ll = calcflux(u_ll, orientation, equation)
  f_rr = calcflux(u_rr, orientation, equation)

  λ_max = sqrt(equation.nu / equation.Tr)

  return 0.5 * (f_ll + f_rr - λ_max * (u_rr - u_ll))
end

@inline function flux_upwind(u_ll, u_rr, orientation, equation::HyperbolicDiffusionEquations2D)
  # Obtain left and right fluxes
  phi_ll, p_ll, q_ll = u_ll
  phi_rr, p_rr, q_rr = u_rr
  f_ll = calcflux(u_ll, orientation, equation)
  f_rr = calcflux(u_rr, orientation, equation)

  # this is an optimized version of the application of the upwind dissipation matrix:
  #   dissipation = 0.5*R_n*|Λ|*inv(R_n)[[u]]
  λ_max = sqrt(equation.nu/equation.Tr)
  f1 = 1/2 * (f_ll[1] + f_rr[1]) - 1/2 * λ_max * (phi_rr - phi_ll)
  if orientation == 1 # x-direction
    f2 = 1/2 * (f_ll[2] + f_rr[2]) - 1/2 * λ_max * (p_rr - p_ll)
    f3 = 1/2 * (f_ll[3] + f_rr[3])
  else # y-direction
    f2 = 1/2 * (f_ll[2] + f_rr[2])
    f3 = 1/2 * (f_ll[3] + f_rr[3]) - 1/2 * λ_max * (q_rr - q_ll)
  end

  return SVector(f1, f2, f3)
end

# Determine maximum stable time step based on polynomial degree and CFL number
function calc_max_dt(u, element_id, invjacobian, cfl,
                     equation::HyperbolicDiffusionEquations2D, dg)
  λ_max = sqrt(equation.nu / equation.Tr)
  dt = cfl * 2 / (nnodes(dg) * invjacobian * λ_max)

  return dt
end

# Convert conservative variables to primitive
@inline cons2prim(u, equation::HyperbolicDiffusionEquations2D) = u

# Convert conservative variables to entropy found in I Do Like CFD, Too, Vol. 1
@inline function cons2entropy(u, equation::HyperbolicDiffusionEquations2D)
  phi, q1, q2 = u
  w1 = phi
  w2 = equation.Lr^2 * q1
  w3 = equation.Lr^2 * q2

  return SVector(w1, w2, w3)
end

# Calculate entropy for a conservative state `u` (here: same as total energy)
@inline entropy(u, equation::HyperbolicDiffusionEquations2D) = energy_total(u, equation)

# Calculate total energy for a conservative state `u`
@inline function energy_total(u, equation::HyperbolicDiffusionEquations2D)
  # energy function as found in equation (2.5.12) in the book "I Do Like CFD, Vol. 1"
  phi, q1, q2 = u
  return 0.5 * (phi^2 + equation.Lr^2 * (q1^2 + q2^2))
end


text:
"""
    Pendulum{R}

The most simple canonical nonlinear system, consisting of a simple pendulum controlled by
torque at the base. Has 2 state and 1 control.

# Constructor
    Pendulum(; kwargs...)

with keyword arguments
* `mass` - mass of the pendulum, in kg (default = 1.0)
* `length` - length of the pendulum, in m (default = 0.5)
* `b` - friction coefficient (default=0.1)
* `I` - inertia, in kg⋅m² (default=0.25)
* `g` - gravity, in kg/m² (default=9.81)
"""

@autodiff mutable struct Pendulum{T} <: ContinuousDynamics
    mass::T
    len::T
    b::T
    lc::T
    I::T
    g::T
end
function Pendulum(;mass=1., len=0.5, b=0.1, lc=0.5, I=0.25, g=9.81)
    T = eltype(promote(mass, len, b, lc, I, g))
    Pendulum{T}(mass, len, b, lc, I, g)
end

RobotDynamics.state_dim(::Pendulum) = 2
RobotDynamics.control_dim(::Pendulum) = 1

function dynamics(p::Pendulum, x, u)
    m = p.mass * p.lc * p.lc
    @SVector [x[2],
              u[1]/m - p.g*sin(x[1])/p.lc - p.b*x[2]/m]
end
function dynamics!(p::Pendulum, ẋ, x, u)
    ẋ .= dynamics(p, SA[x[1], x[2]], SA[u[1]])
end
@inline Base.position(::Pendulum, x) = @SVector zeros(3)
orientation(::Pendulum, x) = expm((pi-x[1])*@SVector [1,0,0.])


text:
# function parallelmuffin(wlt,taut,t,rhot,x,psf,psfadj,p,taup,fty,rhop,taus,s,rhos,mu,spatialwlt,μt,nspat)
#
#     wlt = myidwt(wlt, nspat, taut[:,:,1,:], rhot, t[:,:,1,:], spatialwlt)
#     b = fty + taup + rhop*p + taus + rhos*s
#     wlt_b = wlt + b
#
#     nxy = (size(x))[1]
#     nxypsf = (size(psf))[1]
#     psfcbe = zeros(Complex64,nxy,nxy)
#     psfpad = zeros(Float64,nxy,nxy)
#     psfpad[1:nxypsf,1:nxypsf] = psf[:,:]
#     psfcbe = 1./(abs(fft(psfpad)).^2+mu)
#     x = real(ifft(psfcbe.*fft(wlt_b)))
#
#     # tmp1 = 0.0
#     # tmp2 = zeros(Float64,nxy,nxy)
#     for b in 1:nspat
#                 hx = dwt(x,wavelet(spatialwlt[b]))
#                 tmp = hx - taut[:,:,1,b]/rhot
#                 t[:,:,1,b] = prox_u(tmp,μt/rhot)
#                 taut[:,:,1,b] = taut[:,:,1,b] + rhot*(t[:,:,1,b]-hx)
#                 # tmp1 = vecnorm([tmp2 (hx-(admmst.t)[:,:,z,b])],2)
#                 # tmp2 = (hx-(admmst.t)[:,:,z,b])
#     end
#     # tmp2[:] = 0
#
#
#     tmp = x-taup/rhop
#     p = max(0,tmp)
#     taup = taup + rhop*(p-x)
#
#     return wlt,x,t,taut,p,taup
#
# end

function parallelmuffin(wlt,taut,t,rhot,x,psf,psfadj,p,taup,fty,rhop,taus,s,rhos,mu,spatialwlt,μt,nspat)

    wlt = convert(Array,wlt)
    x = convert(Array,x)
    p = convert(Array,p)
    taup = convert(Array,taup)
    t = convert(Array,taut)
    taut = convert(Array,taut)

    wlt = myidwt(wlt, nspat, taut[:,:,1,:], rhot, t[:,:,1,:], spatialwlt)
    b = fty + taup + rhop*p + taus + rhos*s
    wlt_b = wlt + b

    nxy = (size(x))[1]
    nxypsf = (size(psf))[1]
    psfcbe = zeros(Complex64,nxy,nxy)
    psfpad = zeros(Float64,nxy,nxy)
    psfpad[1:nxypsf,1:nxypsf] = psf[:,:]
    psfcbe = 1./(abs(fft(psfpad)).^2+mu)
    x = real(ifft(psfcbe.*fft(wlt_b)))

    # tmp1 = 0.0
    # tmp2 = zeros(Float64,nxy,nxy)
    for b in 1:nspat
                hx = dwt(x,wavelet(spatialwlt[b]))
                tmp = hx - taut[:,:,1,b]/rhot
                t[:,:,1,b] = prox_u(tmp,μt/rhot)
                taut[:,:,1,b] = taut[:,:,1,b] + rhot*(t[:,:,1,b]-hx)
                # tmp1 = vecnorm([tmp2 (hx-(admmst.t)[:,:,z,b])],2)
                # tmp2 = (hx-(admmst.t)[:,:,z,b])
    end
    # tmp2[:] = 0

    tmp = x-taup/rhop
    p = max(0,tmp)
    taup = taup + rhop*(p-x)

    return wlt,x,t,taut,p,taup

end

@parallel for i in 2:nfreq+1
    z = i-1
    admmst.wlt[:,:,z],admmst.x[:,:,z],admmst.t[:,:,z,:],admmst.taut[:,:,z,:],admmst.p[:,:,z],admmst.taup[:,:,z] =
                                @fetchfrom(i,parallelmuffin(admmst.wlt[:,:,z], admmst.taut[:,:,z,:], admmst.t[:,:,z,:], rhot, admmst.x[:,:,z],
                                psfst.mypsf[:,:,z], psfst.mypsfadj[:,:,z], admmst.p[:,:,z], admmst.taup[:,:,z],
                                fty[:,:,z], rhop, admmst.taus[:,:,z], admmst.s[:,:,z], rhos, admmst.mu, spatialwlt,
                                μt, nspat))
end

image = randn(2048,2048,16);
psf = randn(2048,2048,16);
resultats = SharedArray(Float64,2048,2048,16);

@time @parallel for z in 1:16
    resultats[:,:,z] = (image[:,:,z].*psf[:,:,z])
    psf[:,:,z] = psf[:,:,z]*2
    println(sum(resultats[:,:,z]))
    # rr = RemoteRef(myid())
    # put!(rr,psf[:,:,z])
end

@time @parallel for z in 1:16
    resultats = @fetchfrom(z+1, gc())
    image = @fetchfrom(z+1, gc())
    psf = @fetchfrom(z+1, gc())
end