text:
@testset "probabilityweightedmoment.jl" begin
    @testset "pwm(x, p, r, s)" begin

        # p < 0 throws error
        @test_throws AssertionError Extremes.pwm([1.0], -1, 1, 1)

        # r < 0 throws error
        @test_throws AssertionError Extremes.pwm([1.0], 1, -1, 1)

        #  s < 0 throws error
        @test_throws AssertionError Extremes.pwm([1.0], 1, 1, -1)

        # Probability weighted moment with p=1, r=0 and s=k
        μ = .5
        σ = 1.0
        pd = Logistic(μ,σ)
        y = rand(pd, 10000)
        # theoritical form of the pwm for the Logistic distribution (Greenwood and al., 1979)
        f(k::Int) = μ/(1+k) - σ/(1+k)*sum(1/i for i=1:k)

        for k=1:3
            @test f(k) ≈ Extremes.pwm(y,1,0,k) rtol=.1
        end

        # Probability weighted moment with p=1, r=j and s=0
        μ = .5
        σ = 1.0
        ξ = 0
        pd = GeneralizedExtremeValue(μ,σ,0)
        y = rand(pd, 10000)

        # theoritical form of the pwm for the Gumbel distribution (Greenwood and al., 1979)
        f2(j::Int) = μ/(1+j) + σ/(1+j) * (log(1+j) + MathConstants.eulergamma)

        for j=1:3
            @test f2(j) ≈ Extremes.pwm(y,1,j,0) rtol=.1
        end

    end

    include(joinpath("probabilityweightedmoment", "probabilityweightedmoment_gev_test.jl"))
    include(joinpath("probabilityweightedmoment", "probabilityweightedmoment_gp_test.jl"))
    include(joinpath("probabilityweightedmoment", "probabilityweightedmoment_gumbel_test.jl"))

end


text:
isdigit(x) = Int('0') <= Int(x) <= Int('9')

toint(x) = Int(x) - 48

function luhn(number)
    number = replace(number, " " => "")

    if length(number) < 2
        return false
    end

    sum = 0
    for (i, ch) in enumerate(reverse(number))
        x = Int(ch) - 48
        if x < 0 || x > 9
            return false
        end
        if i % 2 == 0
            x *= 2
            x = x > 9 ? x - 9 : x
        end
        sum += x
    end
    return rem(sum, 10) == 0
end

ext: jl

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    dfstree, dfstree_edges, cotree_edges, circuitrank

"""
    dfstree(adjfunc, graph, root; rootvalue = 0) -> Dict

Return tree node => predecessor dict of the connected component.
"""
function dfstree(adjfunc, graph, root; rootvalue = 0)
    stack = [root]
    pred = Dict(root => rootvalue)
    while !isempty(stack)
        i = pop!(stack)
        for adj in adjfunc(graph, i)
            if !haskey(pred, adj)
                pred[adj] = i
                push!(stack, adj)
            end
        end
    end
    return pred
end

"""
    edgedfstree(nbrfunc, graph, root; rootvalue = 0) -> Dict

Return tree node => edge to the predecessor dict of the connected component.
"""
function edgedfstree(nbrfunc, graph, root; rootvalue = 0)
    stack = [root]
    edges = Dict(root => rootvalue)
    while !isempty(stack)
        i = pop!(stack)
        for (inc, adj) in nbrfunc(graph, i)
            if !haskey(edges, adj)
                edges[adj] = inc
                push!(stack, adj)
            end
        end
    end
    return edges
end

"""
    cotree_edges(graph::UndirectedGraph, root) -> Set{Int}

Return a set of co-tree edges of the connected graph.
"""
function cotree_edges(graph::UndirectedGraph)
    root = pop!(nodeset(graph))
    treeedges = Set(values(edgedfstree(neighbors, graph, root)))
    return setdiff(edgeset(graph), treeedges)
end

"""
    circuitrank(graph::UndirectedGraph) -> Int

Return circuit rank (the number of cycles) of the undirected graph.
"""
function circuitrank(graph::UndirectedGraph)
    nodecount(graph) == 0 && return 0
    hascache(graph, :mincycles) && return length(graph.cache[:mincycles])
    conncount = length(connectedcomponents(graph))
    return edgecount(graph) - nodecount(graph) + conncount
end


text:
"""
    createrun(mlf::MLFlow, experiment_id; start_time=missing, tags=missing)

Creates a run associated to an experiment.

# Arguments
- `mlf`: [`MLFlow`](@ref) configuration.
- `experiment_id`: experiment identifier.

# Keywords
- `start_time`: if provided, must be a UNIX timestamp in milliseconds. By default, set to current time.
- `tags`: if provided, must be a key-value structure such as a dictionary.

# Returns
- an instance of type [`MLFlowRun`](@ref)
"""
function createrun(mlf::MLFlow, experiment_id; start_time=missing, tags=missing)
    endpoint = "runs/create"
    if ismissing(start_time)
        start_time = Int(trunc(datetime2unix(now()) * 1000))
    end
    result = mlfpost(mlf, endpoint; experiment_id=experiment_id, start_time=start_time, tags=tags)
    MLFlowRun(result["run"]["info"], result["run"]["data"])
end
"""
    createrun(mlf::MLFlow, experiment::MLFlowExperiment; start_time=missing, tags=missing)

Dispatches to `createrun(mlf::MLFlow, experiment_id; start_time=start_time, tags=tags)`
"""
createrun(mlf::MLFlow, experiment::MLFlowExperiment; start_time=missing, tags=missing) =
    createrun(mlf, experiment.experiment_id; start_time=start_time, tags=tags)

"""
    getrun(mlf::MLFlow, run_id)

Retrieves information about an MLFlow run.

# Arguments
- `mlf`: [`MLFlow`](@ref) configuration.
- `run_id::String`: run identifier.

# Returns
- an instance of type [`MLFlowRun`](@ref)
"""
function getrun(mlf::MLFlow, run_id)
    endpoint = "runs/get"
    result = mlfget(mlf, endpoint; run_id=run_id)
    MLFlowRun(result["run"]["info"], result["run"]["data"])
end

"""
    updaterun(mlf::MLFlow, run, status; end_time=missing)

Updates the status of an experiment's run.

# Arguments
- `mlf`: [`MLFlow`](@ref) configuration.
- `run`: one of [`MLFlowRun`](@ref), [`MLFlowRunInfo`](@ref), or `String`.
- `status`: either `String` and one of ["RUNNING", "SCHEDULED", "FINISHED", "FAILED", "KILLED"], or an instance of `MLFlowRunStatus`

# Keywords
- `end_time`: if provided, must be a UNIX timestamp in milliseconds. By default, set to current time.
"""
function updaterun(mlf::MLFlow, run_id::String, status::MLFlowRunStatus; end_time=missing)
    endpoint = "runs/update"
    kwargs = Dict(
        :run_id => run_id,
        :status => status.status,
        :end_time => end_time
    )
    if ismissing(end_time) && status.status == "FINISHED"
        end_time = Int(trunc(datetime2unix(now()) * 1000))
        kwargs[:end_time] = string(end_time)
    end
    result = mlfpost(mlf, endpoint; kwargs...)
    MLFlowRun(result["run_info"])
end
updaterun(mlf::MLFlow, run_id::String, status::String; end_time=missing) =
    updaterun(mlf, run_id, MLFlowRunStatus(status); end_time=end_time)
updaterun(mlf::MLFlow, run_info::MLFlowRunInfo, status::String; end_time=missing) =
    updaterun(mlf, run_info.run_id, MLFlowRunStatus(status); end_time=end_time)
updaterun(mlf::MLFlow, run::MLFlowRun, status::String; end_time=missing) =
    updaterun(mlf, run.info, MLFlowRunStatus(status), end_time=end_time)
updaterun(mlf::MLFlow, run_info::MLFlowRunInfo, status::MLFlowRunStatus; end_time=missing) =
    updaterun(mlf, run_info.run_id, status, end_time=end_time)
updaterun(mlf::MLFlow, run::MLFlowRun, status::MLFlowRunStatus; end_time=missing) =
    updaterun(mlf, run.info, status; end_time=end_time)

"""
    deleterun(mlf::MLFlow, run)

Deletes an experiment's run.

# Arguments
- `mlf`: [`MLFlow`](@ref) configuration.
- `run`: one of [`MLFlowRun`](@ref), [`MLFlowRunInfo`](@ref), or `String`.

# Returns
`true` if successful.

"""
function deleterun(mlf::MLFlow, run_id::String)
    endpoint = "runs/delete"
    mlfpost(mlf, endpoint; run_id=run_id)
    true
end
deleterun(mlf::MLFlow, run_info::MLFlowRunInfo) = deleterun(mlf, run_info.run_id)
deleterun(mlf::MLFlow, run::MLFlowRun) = deleterun(mlf, run.info)

"""
    searchruns(mlf::MLFlow, experiment_ids)

Searches for runs in an experiment.

# Arguments
- `mlf`: [`MLFlow`](@ref) configuration.
- `experiment_ids::AbstractVector{Integer}`: `experiment_id`s in which to search for runs. Can also be a single `Integer`.

# Keywords
- `filter::String`: filter as defined in [MLFlow documentation](https://mlflow.org/docs/latest/rest-api.html#search-runs)
- `filter_params::AbstractDict{K,V}`: if provided, `filter` is automatically generated based on `filter_params` using [`generatefilterfromparams`](@ref). One can only provide either `filter` or `filter_params`, but not both.
- `run_view_type::String`: one of `ACTIVE_ONLY`, `DELETED_ONLY`, or `ALL`.
- `max_results::Integer`: 50,000 by default.
- `order_by::String`: as defined in [MLFlow documentation](https://mlflow.org/docs/latest/rest-api.html#search-runs)
- `page_token::String`: paging functionality, handled automatically. Not meant to be passed by the user.

# Returns
- vector of [`MLFlowRun`](@ref) runs that were found in the list of experiments.

"""
function searchruns(mlf::MLFlow, experiment_ids::AbstractVector{<:Integer};
                    filter::String="",
                    filter_params::AbstractDict{K,V}=Dict{}(),
                    run_view_type::String="ACTIVE_ONLY",
                    max_results::Int64=50000,
                    order_by::AbstractVector{<:String}=["attribute.end_time"],
                    page_token::String=""
                    ) where {K,V}
    endpoint = "runs/search"
    run_view_type ∈ ["ACTIVE_ONLY", "DELETED_ONLY", "ALL"] || error("Unsupported run_view_type = $run_view_type")

    if length(filter_params) > 0 && length(filter) > 0
        error("Can only use either filter or filter_params, but not both at the same time.")
    end

    if length(filter_params) > 0
        filter = generatefilterfromparams(filter_params)
    end

    kwargs = (
        experiment_ids=experiment_ids,
        filter=filter,
        run_view_type=run_view_type,
        max_results=max_results,
        order_by=order_by
    )
    if !isempty(page_token)
        kwargs = (; kwargs..., page_token=page_token)
    end

    result = mlfpost(mlf, endpoint; kwargs...)
    haskey(result, "runs") || return MLFlowRun[]

    runs = map(x -> MLFlowRun(x["info"], x["data"]), result["runs"])

    # paging functionality using recursion
    if haskey(result, "next_page_token") && !isempty(result["next_page_token"])
        kwargs = (
            filter=filter,
            run_view_type=run_view_type,
            max_results=max_results,
            order_by=order_by,
            page_token=result["next_page_token"]
        )
        nextruns = searchruns(mlf, experiment_ids; kwargs...)
        return vcat(runs, nextruns)
    end

    runs
end
searchruns(mlf::MLFlow, experiment_id::Integer; kwargs...) =
    searchruns(mlf, [experiment_id]; kwargs...)
searchruns(mlf::MLFlow, exp::MLFlowExperiment; kwargs...) =
    searchruns(mlf, exp.experiment_id; kwargs...)
searchruns(mlf::MLFlow, exps::AbstractVector{MLFlowExperiment}; kwargs...) =
    searchruns(mlf, [getfield.(exps, :experiment_id)]; kwargs...)

"""
    logparam(mlf::MLFlow, run, key, value)
    logparam(mlf::MLFlow, run, kv)

Associates a key/value pair of parameters to the particular run.

# Arguments
- `mlf`: [`MLFlow`](@ref) configuration.
- `run`: one of [`MLFlowRun`](@ref), [`MLFlowRunInfo`](@ref), or `String`.
- `key`: parameter key (name). Automatically converted to string before sending to MLFlow because this is the only type that MLFlow supports.
- `value`: parameter value. Automatically converted to string before sending to MLFlow because this is the only type that MLFlow supports.

One could also specify `kv::Dict` instead of separate `key` and `value` arguments.
"""
function logparam(mlf::MLFlow, run_id::String, key, value)
    endpoint = "runs/log-parameter"
    mlfpost(mlf, endpoint; run_id=run_id, key=string(key), value=string(value))
end
logparam(mlf::MLFlow, run_info::MLFlowRunInfo, key, value) =
    logparam(mlf, run_info.run_id, key, value)
logparam(mlf::MLFlow, run::MLFlowRun, key, value) =
    logparam(mlf, run.info, key, value)
function logparam(mlf::MLFlow, run::Union{String,MLFlowRun,MLFlowRunInfo}, kv)
    for (k, v) in kv
        logparam(mlf, run, k, v)
    end
end

"""
    logmetric(mlf::MLFlow, run, key, value::T; timestamp, step) where T<:Real
    logmetric(mlf::MLFlow, run, key, values::AbstractArray{T}; timestamp, step) where T<:Real

Logs a metric value (or values) against a particular run.

# Arguments
- `mlf`: [`MLFlow`](@ref) configuration.
- `run`: one of [`MLFlowRun`](@ref), [`MLFlowRunInfo`](@ref), or `String`
- `key`: metric name.
- `value`: metric value, must be numeric.

# Keywords
- `timestamp`: if provided, must be a UNIX timestamp in milliseconds. By default, set to current time.
- `step`: step at which the metric value has been taken.
"""
function logmetric(mlf::MLFlow, run_id::String, key, value::T; timestamp=missing, step=missing) where T<:Real
    endpoint = "runs/log-metric"
    if ismissing(timestamp)
        timestamp = Int(trunc(datetime2unix(now()) * 1000))
    end
    mlfpost(mlf, endpoint; run_id=run_id, key=key, value=value, timestamp=timestamp, step=step)
end
logmetric(mlf::MLFlow, run_info::MLFlowRunInfo, key, value::T; timestamp=missing, step=missing) where T<:Real =
    logmetric(mlf::MLFlow, run_info.run_id, key, value; timestamp=timestamp, step=step)
logmetric(mlf::MLFlow, run::MLFlowRun, key, value::T; timestamp=missing, step=missing) where T<:Real =
    logmetric(mlf, run.info, key, value; timestamp=timestamp, step=step)

function logmetric(mlf::MLFlow, run::Union{String,MLFlowRun,MLFlowRunInfo}, key, values::AbstractArray{T}; timestamp=missing, step=missing) where T<:Real
    for v in values
        logmetric(mlf, run, key, v; timestamp=timestamp, step=step)
    end
end

"""
    logartifact(mlf::MLFlow, run, basefilename, data)

Stores an artifact (file) in the run's artifact location.

!!! note
    Assumes that artifact_uri is mapped to a local directory.
    At the moment, this only works if both MLFlow and the client are running on the same host or they map a directory that leads to the same location over NFS, for example.

# Arguments
- `mlf::MLFlow`: [`MLFlow`](@ref) onfiguration. Currently not used, but when this method is extended to support `S3`, information from `mlf` will be needed.
- `run`: one of [`MLFlowRun`](@ref), [`MLFlowRunInfo`](@ref) or `String`.
- `basefilename`: name of the file to be written.
- `data`: artifact content, an object that can be written directly to a file handle.

# Throws
- an `ErrorException` if an exception occurs during writing artifact.

# Returns
path of the artifact that was created.
"""
function logartifact(mlf::MLFlow, run_id::AbstractString, basefilename::AbstractString, data)
    mlflowrun = getrun(mlf, run_id)
    artifact_uri = mlflowrun.info.artifact_uri
    mkpath(artifact_uri)
    filepath = joinpath(artifact_uri, basefilename)
    try
        f = open(filepath, "w")
        write(f, data)
        close(f)
    catch e
        error("Unable to create artifact $(filepath): $e")
    end
    filepath
end
logartifact(mlf::MLFlow, run::MLFlowRun, basefilename::AbstractString, data) =
    logartifact(mlf, run.info, basefilename, data)
logartifact(mlf::MLFlow, run_info::MLFlowRunInfo, basefilename::AbstractString, data) =
    logartifact(mlf, run_info.run_id, basefilename, data)

"""
    logartifact(mlf::MLFlow, run, filepath)

Stores an artifact (file) in the run's artifact location.
The name of the artifact is calculated using `basename(filepath)`.

Dispatches on `logartifact(mlf::MLFlow, run, basefilename, data)` where `data` is the contents of `filepath`.

# Throws
- an `ErrorException` if `filepath` does not exist.
- an exception if such occurs while trying to read the contents of `filepath`.

"""
function logartifact(mlf::MLFlow, run_id::AbstractString, filepath::Union{AbstractPath,AbstractString})
    isfile(filepath) || error("File $filepath does not exist.")
    try
        f = open(filepath, "r")
        data = read(f)
        close(f)
        return logartifact(mlf, run_id, basename(filepath), data)
    catch e
        throw(e)
    finally
        if @isdefined f
            close(f)
        end
    end
end
logartifact(mlf::MLFlow, run::MLFlowRun, filepath::Union{AbstractPath,AbstractString}) =
    logartifact(mlf, run.info, filepath)
logartifact(mlf::MLFlow, run_info::MLFlowRunInfo, filepath::Union{AbstractPath,AbstractString}) =
    logartifact(mlf, run_info.run_id, filepath)

"""
    listartifacts(mlf::MLFlow, run)

Lists the artifacts associated with an experiment run.
According to [MLFlow documentation](https://mlflow.org/docs/latest/rest-api.html#list-artifacts), this API endpoint should return paged results, similar to [`searchruns`](@ref).
However, after some experimentation, this doesn't seem to be the case. Therefore, the paging functionality is not implemented here.

# Arguments
- `mlf::MLFlow`: [`MLFlow`](@ref) onfiguration. Currently not used, but when this method is extended to support `S3`, information from `mlf` will be needed.
- `run`: one of [`MLFlowRun`](@ref), [`MLFlowRunInfo`](@ref) or `String`.

# Keywords
- `path::String`: path of a directory within the artifact location. If set, returns the contents of the directory. By default, this is the root directory of the artifacts.
- `maxdepth::Int64`: depth of listing. Default is 1. This will only return the files/directories in the current `path`. To return all artifacts files and directories, use `maxdepth=-1`.

# Returns
A vector of `Union{MLFlowArtifactFileInfo,MLFlowArtifactDirInfo}`.
"""
function listartifacts(mlf::MLFlow, run_id::String; path::String="", maxdepth::Int64=1)
    endpoint = "artifacts/list"
    kwargs = (
        run_id=run_id,
    )
    kwargs = (; kwargs..., path=path)
    httpresult = mlfget(mlf, endpoint; kwargs...)
    "files" ∈ keys(httpresult) || return Vector{Union{MLFlowArtifactFileInfo,MLFlowArtifactDirInfo}}()
    "root_uri" ∈ keys(httpresult) || error("Malformed response from MLFlow REST API.")
    root_uri = httpresult["root_uri"]
    result = Vector{Union{MLFlowArtifactFileInfo,MLFlowArtifactDirInfo}}()
    maxdepth == 0 && return result

    for resultentry ∈ httpresult["files"]
        if resultentry["is_dir"] == false
            filepath = joinpath(root_uri, resultentry["path"])
            file_size = resultentry["file_size"]
            if typeof(file_size) <: Int
                filesize = file_size
            else
                filesize = parse(Int, file_size)
            end
            push!(result, MLFlowArtifactFileInfo(filepath, filesize))
        elseif resultentry["is_dir"] == true
            dirpath = joinpath(root_uri, resultentry["path"])
            push!(result, MLFlowArtifactDirInfo(dirpath))
            if maxdepth != 0
                nextdepthresult = listartifacts(mlf, run_id, path=resultentry["path"], maxdepth=maxdepth-1)
                result = vcat(result, nextdepthresult)
            end
        else
            isdirval = resultentry["is_dir"]
            @warn "Malformed response from MLFlow REST API is_dir=$isdirval - skipping"
            continue
        end
    end
    result
end
listartifacts(mlf::MLFlow, run::MLFlowRun; kwargs...) =
    listartifacts(mlf, run.info.run_id; kwargs...)
listartifacts(mlf::MLFlow, run_info::MLFlowRunInfo; kwargs...) =
    listartifacts(mlf, run_info.run_id; kwargs...)

max_stars_repo_name: JuliaAI/MLFlowClient

max_issues_repo_name: JuliaAI/MLFlowClient

max_forks_repo_name: JuliaAI/MLFlowClient


text:
function initial_linklist(p, data :: Data)

	@unpack N, cutoff, l = data

	first_atom = zeros(Int64, l, l)
	next_atom = zeros(Int64, N)

	for i in 1:N

		ic = trunc(Int64, p[i][1]/cutoff) + 1
		jc = trunc(Int64, p[i][2]/cutoff) + 1

		next_atom[i] = first_atom[ic, jc]
		first_atom[ic, jc] = i
	end

	return first_atom, next_atom
end


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

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

# ╔═╡ 0bf7ad90-e88f-11ea-147b-eb2f2eb606be
using PlutoUI, ShockTube, Unitful

# ╔═╡ 2f4c7d80-e893-11ea-3771-d3c8088a873a
md"## Calculating from flow rate settings"

# ╔═╡ 1f7b19b0-e88f-11ea-27b0-8db91c84aedb
Ace = Species("Acetone");

# ╔═╡ 27f9c870-e88f-11ea-0808-4ba3e41fc037
He = Species("Helium");

# ╔═╡ a6076a9e-e890-11ea-232c-61bbfd839410
begin
	Ace_rate_slider = @bind Ace_rate NumberField(LinRange(0, 20, 1001), default=10)
	He_rate_slider = @bind He_rate NumberField(LinRange(0, 50, 1001), default=30)

	md"""
	Acetone flow rate (sccm) $Ace_rate_slider

	Helium flow rate (slpm)  $He_rate_slider
	"""
end

# ╔═╡ 8cc5a110-e8b8-11ea-2cec-7f83d4d4f7ba
v̇Ace = Ace_rate*u"cm^3/minute"

# ╔═╡ 96e359ce-e8b8-11ea-29ae-d3c78b2e1dba
v̇He = He_rate*u"L/minute"

# ╔═╡ 6db8dae0-e88f-11ea-1355-bbce79b13063
ṅAce = ρm(Ace) * v̇Ace |> u"mol/s"

# ╔═╡ 70091700-e890-11ea-0307-31689d3e8e9f
ṅHe = ρm(He) * v̇He |> u"mol/s"

# ╔═╡ 85bf68b0-e890-11ea-04d6-e1cfff79278d
χAce = ṅAce / (ṅAce + ṅHe)

# ╔═╡ 3e9944d0-e893-11ea-30de-332f76096802
md"## Calculating from desired flow rate & acetone concentration"

# ╔═╡ b53e58a0-e893-11ea-111e-37b14cf3dcab
begin
	v̇_slider = @bind v̇ NumberField(LinRange(0, 50, 101))
	χAce_slider = @bind χAce_in NumberField(LinRange(0, 0.25, 101))
	md"""
	Net volume flow rate (slpm) $v̇_slider

	Acetone fraction $χAce_slider
	"""
end

# ╔═╡ 692684b0-e893-11ea-2af9-756ab48bcdd0
HeAce = Mixture(["Helium", "Acetone"], zs=[1 - χAce_in, χAce_in])

# ╔═╡ 3bcce8f0-e894-11ea-0ff0-25d9c9785a15
ṅ = ρm(HeAce) * v̇*u"L/minute" |> u"mol/s"

# ╔═╡ ab9018a0-e895-11ea-10ab-03e7704de77d
Ace_slpm = ṅ * χAce_in / ρm(Ace) |> u"cm^3 / minute"

# ╔═╡ 6b7a74e0-e895-11ea-015c-aba28e723af1
He_slpm = ṅ * (1 - χAce_in) / ρm(He) |> u"L/minute"

# ╔═╡ 209a1e20-eef9-11ea-2e0a-b3b5b43861b6
md""" ## Shock jump conditions for helium/acetone mixture"""

# ╔═╡ c238a120-eef9-11ea-35cd-f744665ba19a
begin
	Mach_slider = @bind Mach NumberField(range(1.1, 2.8, step=0.01))
	md"""
	Mach number $Mach_slider
	"""
end

# ╔═╡ 1da4d2f0-eef9-11ea-0752-cf0358151391
driver_gas = Species("N2", T=800)

# ╔═╡ 13cbc450-eef9-11ea-24e0-0fa6efd42922
driver_pressure = p(shockcalc(driver_gas, HeAce, Mach).driver) |> u"psi"

# ╔═╡ Cell order:
# ╟─2f4c7d80-e893-11ea-3771-d3c8088a873a
# ╠═0bf7ad90-e88f-11ea-147b-eb2f2eb606be
# ╠═1f7b19b0-e88f-11ea-27b0-8db91c84aedb
# ╠═27f9c870-e88f-11ea-0808-4ba3e41fc037
# ╟─a6076a9e-e890-11ea-232c-61bbfd839410
# ╟─8cc5a110-e8b8-11ea-2cec-7f83d4d4f7ba
# ╟─96e359ce-e8b8-11ea-29ae-d3c78b2e1dba
# ╟─6db8dae0-e88f-11ea-1355-bbce79b13063
# ╟─70091700-e890-11ea-0307-31689d3e8e9f
# ╟─85bf68b0-e890-11ea-04d6-e1cfff79278d
# ╟─3e9944d0-e893-11ea-30de-332f76096802
# ╟─b53e58a0-e893-11ea-111e-37b14cf3dcab
# ╠═692684b0-e893-11ea-2af9-756ab48bcdd0
# ╟─3bcce8f0-e894-11ea-0ff0-25d9c9785a15
# ╟─ab9018a0-e895-11ea-10ab-03e7704de77d
# ╟─6b7a74e0-e895-11ea-015c-aba28e723af1
# ╟─209a1e20-eef9-11ea-2e0a-b3b5b43861b6
# ╠═c238a120-eef9-11ea-35cd-f744665ba19a
# ╠═1da4d2f0-eef9-11ea-0752-cf0358151391
# ╠═13cbc450-eef9-11ea-24e0-0fa6efd42922


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -12.548880784274985 <= q <= 15.06837792907793)

                     add_NL_constraint(m, :(swish(0.4386422028907555 + 0.03049853723546203*swish(-0.33806544253006265 + -0.5419050015026476*swish(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.12420368558845318*swish(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.008939690927983968*swish(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + -0.7963936693468932*swish(0.35455126051631725 + -0.3141650319667155*swish(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + 0.26311214619974344*swish(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + -0.4178451532515153*swish(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + -0.19871150900740808*swish(-0.9712553762218641 + -0.9788234903921222*swish(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.6598676192535482*swish(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.7900118209668339*swish(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4])))) + swish(0.039843038969352484 + 0.2885956651060697*swish(-0.33806544253006265 + -0.5419050015026476*swish(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.12420368558845318*swish(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.008939690927983968*swish(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + 0.8477392827604406*swish(0.35455126051631725 + -0.3141650319667155*swish(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + 0.26311214619974344*swish(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + -0.4178451532515153*swish(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + -0.38122333260621444*swish(-0.9712553762218641 + -0.9788234903921222*swish(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.6598676192535482*swish(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.7900118209668339*swish(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4])))) + swish(0.5234531332153054 + -0.26529753551452195*swish(-0.33806544253006265 + -0.5419050015026476*swish(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.12420368558845318*swish(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.008939690927983968*swish(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + 0.29155449146294243*swish(0.35455126051631725 + -0.3141650319667155*swish(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + 0.26311214619974344*swish(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + -0.4178451532515153*swish(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + -0.357141816330512*swish(-0.9712553762218641 + -0.9788234903921222*swish(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.6598676192535482*swish(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.7900118209668339*swish(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module LightDistributions

using Random
using StatsFuns: poisinvcdf
using StaticArrays
using CUDAnative
using CuArrays
using MacroTools

# TODO extend Base.rand instead of defining own random
# import Base: rand

include("specfuns.jl")

############################
# Distribution generalities
############################

const _distributions = (
    :Categorical, :Poisson,
    :Uniform, :Normal, :Exponential, :Gamma,
    :Dirichlet
)
const _params = (
    Categorical = NamedTuple{(:p,), Tuple{Tuple{Vararg{Real}}}},
    Poisson     = NamedTuple{(:λ,), Tuple{Real}},
    Uniform     = NamedTuple{(:a,:b), Tuple{Real,Real}},
    Normal      = NamedTuple{(:μ,:σ), Tuple{Real,Real}},
    Exponential = NamedTuple{(:λ,), Tuple{Real}},
    Gamma       = NamedTuple{(:k,:θ), Tuple{Real,Real}},
    Dirichlet   = NamedTuple{(:α,), Tuple{Tuple{Vararg{Real}}}}
)

export Distribution
export sample, random, params, logpdf
for s in _distributions
    @eval export $s
end

abstract type Distribution{T,P} end

Base.eltype(::Distribution{T}) where T = T
Base.eltype(::Type{<:Distribution{T}}) where T = T

params(d::Distribution) = d.params

###############################
# Distribution implementations
###############################

# function homogenize(t::Tuple)
#     T = promote_type(typeof(t).parameters...)
#     Tuple{Vararg{T}}(t)
# end
#
# function homogenize(s::Tuple, t::Tuple)
#     T = promote_type(typeof(t).parameters..., typeof(t).parameters...)
#     Tuple{Vararg{T}}(s), Tuple{Vararg{T}}(t)
# end

include("scalars.jl")
include("arrays.jl")

for D in _distributions
    @eval params(::Type{<:$D}) = $(_params[D])
    @eval logpdf(::Type{<:$D}) = $(Symbol(:logpdf,D))
    @eval sample(::Type{<:$D}) = $(Symbol(:rand, D))
    @eval random(d::$D) = $(Symbol(:rand, D))(d.params...)
end

end # module


text:
function iscoherentmarket(f::Vector{Float64}, t::Vector{<:Real})
    length(f) == length(t) || throw(DimensionMismatch("`f` and `t` must have same length"))
    all(0 .< f .≤ 1) || throw(DomainError(f, "Must have `all(0 .< f .≤ 1)`"))
    all(0 .≤ t) || throw(DomainError(f, "Must have `all(0 .≤ t)`"))
end

function iscoherentmarket(f::Vector{Float64}, t::Vector{Int}, g::Vector{Int})
    length(f) == length(t) || throw(DimensionMismatch("`f` and `t` must have same length"))
    length(f) == length(g) || throw(DimensionMismatch("`f` and `g` must have same length"))
    all(0 .< f .≤ 1) || throw(DomainError(f, "Must have `all(0 .< f .≤ 1)`"))
    all(0 .≤ t) || throw(DomainError(f, "Must have `all(0 .≤ t)`"))
    all(0 .≤ g) || throw(DomainError(f, "Must have `all(0 .≤ g)`"))
end

"""
Contains information about a college application market.
"""
abstract type Market end

"""
    SameCostsMarket{U<:Real}

Contains information about a college application market with identical
application costs.

# Constructors

```julia
SameCostsMarket(f::Vector{Float64}, t::Vector{U}, h::Integer) where {U<:Real}
```
Construct the `SameCostsMarket` defined by admissions probabilities `f`,
utility values `t`, and application limit `h`.

```julia-repl
julia> SameCostsMarket([0.3, 0.2, 0.05], [2, 3, 4], 2)
SameCostsMarket{Int64}(3, [0.3, 0.2, 0.05], [2, 3, 4], 2, [0.6, 0.6000000000000001, 0.2], [0.7, 0.8, 0.95], [1, 2, 3])
```

```julia
SameCostsMarket(m)
```
Generate a random `SameCostsMarket{Int}` with `m` schools. `f` and `t` correlate negatively, mimicking
a realistic market.

# Internal API

This type contains the following fields:

 - `m`: Number of schools
 - `f`: Vector of admissions probabilities
 - `t`: Vector of utility values (must be sorted)
 - `h`: Number of schools student is allowed to apply to
 - `ft = f .* t`
 - `omf = 1 .- f`
 - `perm`: How the input data were permuted to sort by `t`

`U` is the eltype of `t`.

Internally, the schools must be indexed so that `t` is sorted ascending. `perm`
allows the original order to be recovered:

```julia-repl
julia> f = rand(3); t = rand(3); mkt = SameCostsMarket(f, t, 2);

julia> t[mkt.perm] == mkt.t
true

julia> t == mkt.t[invperm(mkt.perm)]
true
```
"""
struct SameCostsMarket{U<:Real} <: Market
    m::Int
    f::Vector{Float64}
    t::Vector{U}
    h::Int
    ft::Vector{Float64}      # = f .* t
    omf::Vector{Float64}     # = 1 .- f
    perm::Vector{Int}

    function SameCostsMarket(f::Vector{Float64}, t::Vector{U}, h::Integer) where {U<:Real}
        iscoherentmarket(f, t)
        perm = sortperm(t)

        return new{U}(
            length(f),
            f[perm],
            t[perm],
            min(h, length(f)),
            (f.*t)[perm],
            (1 .- f)[perm],
            perm)
    end

    function SameCostsMarket(m::Integer)
        t = ceil.(Int, -10 * log.(rand(m))) |> sort
        f = inv.(t .+ 10 * rand(m))
        perm = 1:m
        return new{Int}(
            m,
            f[perm],
            t[perm],
            m ÷ 2,
            (f.*t)[perm],
            (1 .- f)[perm],
            perm)
    end
end

"""
    VariedCostsMarket

Contains information about a college application market with varied application
costs.

# Constructors

```julia
VariedCostsMarket(f::Vector{Float64}, t::Vector{Int}, g::Vector{Int}, H::Int)
```
Construct the `VariedCostsMarket` defined by admissions probabilities `f`,
utility values `t`, application costs `g`, and application limit `h`.

```julia-repl
julia> mkt = VariedCostsMarket([0.1, 0.5, 0.3, 0.1], [12, 3, 4, 13], [1, 2, 1, 1], 3)
VariedCostsMarket(4, [0.5, 0.3, 0.1, 0.1], [3, 4, 12, 13], [2, 1, 1, 1], 3, [1.5, 1.2, 1.2000000000000002, 1.3], [0.5, 0.7, 0.9, 0.9], [2, 3, 1, 4])
```

```julia
VariedCostsMarket(m)
```
Generate a random `SameCostsMarket{Int}` with `m` schools. `f` and `t` correlate negatively, mimicking
a realistic market. Entries of `g` are random integers in `5:10`.

# Internal API

This type contains the following fields:

 - `m`: Number of schools
 - `f`: Vector of admissions probabilities
 - `t`: Vector of utility values (must be sorted)
 - `g`: Vector of application costs
 - `H`: Budget to spend on applications
 - `ft = f .* t`
 - `omf = 1 .- f`
 - `perm`: How the input data were permuted to sort by `t`

Internally, the schools must be indexed so that `t` is sorted ascending. `perm`
allows the original order to be recovered:

```julia-repl
julia> f = rand(3); t = rand(3); g = rand(5:10, 3); mkt = VariedCostsMarket(f, t, g, 10);

julia> t[mkt.perm] == mkt.t
true

julia> t == mkt.t[invperm(mkt.perm)]
true
```
"""
struct VariedCostsMarket <: Market
    m::Int
    f::Vector{Float64}
    t::Vector{Int}
    g::Vector{Int}
    H::Int
    ft::Vector{Float64}      # = f .* t
    omf::Vector{Float64}     # = 1 .- f
    perm::Vector{Int}

    function VariedCostsMarket(f::Vector{Float64}, t::Vector{Int}, g::Vector{Int}, H::Int)
        iscoherentmarket(f, t, g)
        perm = sortperm(t)

        return new(
            length(f),
            f[perm],
            t[perm],
            g[perm],
            min(H, sum(g)),
            (f.*t)[perm],
            (1 .- f)[perm],
            perm)
    end

    function VariedCostsMarket(m::Integer)
        t = ceil.(Int, -10 * log.(rand(m))) |> sort
        f = inv.(t .+ 10 * rand(m))
        g = rand(5:10, m)
        H = sum(g) ÷ 2
        perm = 1:m

        return new(
            m,
            f[perm],
            t[perm],
            g[perm],
            H,
            (f.*t)[perm],
            (1 .- f)[perm],
            perm)
    end
end

"""
    Market(f, t, h)

Equivalent to [`SameCostsMarket(f, t, h)`](@ref).
"""
function Market(f::Vector{Float64}, t::Vector{<:Real}, h::Int)
    return SameCostsMarket(f, t, h)
end

"""
    Market(f, t, g, H)

Equivalent to [`VariedCostsMarket(f, t, g, H)`](@ref).
"""
function Market(f::Vector{Float64}, t::Vector{Int}, g::Vector{Int}, H::Int)
    return VariedCostsMarket(f, t, g, H)
end

function valuation_nopermute_sorted(
    X::AbstractVector{<:Integer},
    mkt::Market;
)::Float64
    isempty(X) && return 0.0
    # @assert issorted(X)
    h = length(X)
    if h > 1
        res = mkt.ft[X[end]]
        cp = reverse(cumprod(reverse(mkt.omf[X[2:end]])))
        for j in 1:h-1
            res += mkt.ft[X[j]] * cp[j]
        end
        return res
    else
        return mkt.ft[X[1]]
    end
end

"""
    valuation(X, mkt)

Return the valuation of the portfolio `X` for the market `mkt`, which may be either
a `SameCostsMarket` or a `VariedCostsMarket`.

`SameCostsMarket` example (`h = 3` is irrelevant):

```julia-repl
julia> mkt = SameCostsMarket([0.1, 0.5, 0.3, 0.1], [12, 3, 4, 13], 3);

julia> round(valuation([1, 4], mkt), digits=2) # expected utility when applying to schools 1 and 4
2.38
```

`VariedCostsMarket` example (`g = [1, 2, 1, 1]` and `H = 4` are irrelevant):

```julia-repl
julia> mkt = VariedCostsMarket([0.1, 0.5, 0.3, 0.1], [12, 3, 4, 13], [1, 2, 1, 1], 4);

julia> round(valuation([1, 4], mkt), digits=2) # expected utility when applying to schools 1 and 4
2.38
```
"""
function valuation(
    X::AbstractVector{<:Integer},
    mkt::Market;
    invp::AbstractVector{<:Integer}=invperm(mkt.perm)
)::Float64
    return valuation_nopermute_sorted(sort(invp[X]), mkt)
end

max_stars_repo_name: maxkapur/OptimalApplication

max_issues_repo_name: maxkapur/OptimalApplication

max_forks_repo_name: maxkapur/OptimalApplication


text:
using PerceptualColourMaps
using Test, LinearAlgebra, Statistics

include("test_cmap.jl")
include("test_applycolourmaps.jl")
include("test_relief.jl")
include("test_utilities.jl")


text:
#===============================================================================

            The main macro of this package `@diffusion_process`
            used for defining diffusion processes

===============================================================================#
# define some names for the document
_DIMENSIONS = [:dimensions, :dim, :dims, :dimension]
_PARAMETERS = [:parameters, :param, :params]
_CONJUGATE = [:conjugate]
_EXTRA = [:additional, :extra]
_WIENER = [:wiener, :noise, :gaus, :gaussian]
_PROCESS = [:proc, :process, :state, :statespace]
_STATESPACE = [:statespace, :state, :domain, :domains]
_CONSTDIFF = [
    :constdiff,
    :constvola,
    :constdiffusivity,
    :constvolatility,
    :constantdiff,
    :constantvola,
]
_LINEAR = [:linear, :lineardiffusion]
_ELTYPE = [:eltype]
_NUMNONHYPO = [:num_non_hypo, :numnonhypo]
_PHI = [:phi, :ϕ, :φ]

"""
    Base.lowercase(s::Symbol)

Lowercase all letters in a symbol
"""
Base.lowercase(s::Symbol) = Symbol(lowercase(string(s)))

"""
    _symbol_in(s::Symbol, symbols)

Check if symbol `s` is listed in a list of symbols `symbols`
"""
_symbol_in(s::Symbol, symbols) = lowercase(s) in symbols

"""
    _symbol_in(s::QuoteNode, symbols)

Check if quote of a symbol `s` is listed in a list of symbols `symbols`
"""
_symbol_in(s::QuoteNode, symbols) = _symbol_in(eval(s), symbols)

"""
    _symbol_in(::Any, ::Any)

Return `false` by default i.e. if the first argument is not a symbol nor its
quote
"""
_symbol_in(::Any, ::Any) = false

"""
    diffusion_process(name, ex::Expr, p...)

Defines a diffusion process according to a template described in the
documentation of the github repository:
https://github.com/mmider/DiffusionDefinition.jl
"""
macro diffusion_process(name, ex::Expr, p...)
    parse_process(name, MacroTools.striplines(ex), p)
end

"""
    load_diffusion()

Displays available choices of predefined diffusion processes that can be loaded
"""
macro load_diffusion()
    println("The following diffusions have been defined inside the package ",
            "DiffusionDefinition.jl and can be loaded by calling a macro ",
            "@load_diffusion name-of-a-diffusion:")
    for name in keys(_NAMES_TO_ACCEPTEDNAMES)
        println("- $name")
    end
    nothing
end

"""
    load_diffusion(name)

Loads the predefined diffusion process.
"""
macro load_diffusion(name)
    if _symbol_in(name, _ADMISSIBLENAMES)
        if typeof(name) <: QuoteNode
            name = eval(name)
        end
        homedir = joinpath(@__DIR__, "..", "examples")
        importname = _ACCEPTEDNAMES_TO_NAMES[lowercase(name)]
        filepath = _NAMES_TO_PATH[importname]
        path = joinpath(homedir, filepath)
        isfile(path) && include(path)
        !isfile(path) && println("Error, diffusion $name is supposed to be ",
                                "defined but it seems the file does not exist")
        return Meta.parse("import DiffusionDefinition.$importname")
    else
        println("Diffusion $name does not seem to be defined...")
    end
    nothing
end

"""
    parse_process(name , ex::Expr, ::Any)

Parse a template defining a diffusion process, create a corresponding struct and
specified functions, evaluate them in the environment of a package and then
import the struct name to `Main` scope, in which the package has been imported
to.
"""
function parse_process(name , ex::Expr, ::Any)
    name, curly_name, template_args = process_name(name)
    p = (
        name = name,
        curly_name = curly_name,
        template_args = template_args,
        parameters = Vector{Tuple{Symbol,Union{Symbol,Expr},Symbol,Int64}}(
            undef,
            0
        ),
        dims = Dict{Symbol,Int64}(),
        extras = Dict{Symbol, Any}(),
        fns = Vector{Expr}(undef, 0),
    )

    # parse lines defining parameters to get all parameter names first
    parse_lines!(ex, p, x->(x == _PARAMETERS[1]))
    # parse all other lines
    parse_lines!(ex, p, x->(x != _PARAMETERS[1]))
    fill_unspecified_with_defaults(p)

    abstract_type = prepare_abstract_type(
        p.extras[_LINEAR[1]] ? :LinearDiffusion : :DiffusionProcess,
        p.dims,
        p.extras[_ELTYPE[1]],
        p.extras[_STATESPACE[1]],
    )
    struct_def = createstruct(abstract_type, p)
    add_constdiff_function!(p.fns, p)
    eval(struct_def)
    for fn in p.fns
        eval(fn)
    end
    Meta.parse("import DiffusionDefinition.$name")
end

"""
    process_name

Process the name of a struct by returning the pure name (without template
arguments), the name with template arguments and a list of template arguments.
"""
function process_name end
process_name(name::Symbol) = name, name, Any[]
process_name(name::Expr) = name.args[1], name, name.args[2:end]

"""
    parse_lines!(ex::Expr, p, condition)

Parse all lines of the expression `ex`, but process only those which satisfy
`condition`. `p` is a passed-around structure that accumulates processed
information.
"""
function parse_lines!(ex::Expr, p, condition)
    current_label = nothing
    for line in ex.args
        current_label = update_label(line, current_label)
        typeof(line) <: Union{Symbol, QuoteNode} && (continue)
        @assert line.head == :-->
        condition(current_label) && parse_line!(Val{current_label}(), line, p)
    end
end

"""
    update_label(line, current_label)

Update the label, which signifies what type of information a given line in a
template is supposed to be encoding.
"""
function update_label(line, current_label)
    _symbol_in(line, _DIMENSIONS) && return _DIMENSIONS[1]
    _symbol_in(line, _PARAMETERS) && return _PARAMETERS[1]
    _symbol_in(line, _CONJUGATE) && return _CONJUGATE[1]
    _symbol_in(line, _EXTRA) && return _EXTRA[1]
    current_label
end

#------------------------------------------------------------------------------#
#
#                          For parsing: PARAMETERS
#
#------------------------------------------------------------------------------#
"""
    parse_line!(::Val{:parameters}, line, p)

Parse a line that defines parameters of the diffusion. The line must be in a
format:
    name --> parameter-description
"""
function parse_line!(::Val{:parameters}, line, p)
    param_name = line.args[1]

    if typeof(param_name) <: Symbol
        parse_param_single_name(line, p)
    else
        @assert param_name.head == :tuple
        parse_param_multi_names(line, p)
    end
end

"""
    parse_param_single_name(line, p)

Parse a line that defines parameters of the diffusion. The line must be in one
of the formats:
    parameter_name --> (number_of_parameters, datatype)
    parameter_name --> datatype
In the former case defines `number_of_parameters`-many parameters, with names
`parameter_name`i and of `datatype` type. In the latter case defines a
single parameters with name `parameter_name` and of type `datatype`.
"""
function parse_param_single_name(line, p)
    name_stem, disambig_idx, generic_name = get_name_stem(
        line.args[1], p.parameters
    )

    data_type = line.args[2]

    if _is_datatype(data_type, p)
        num_params = 1
    else
        @assert data_type.head == :tuple
        num_params, data_type = data_type.args[1], data_type.args[2]
        generic_name = true
    end

    for i in 1:num_params
        name = ( generic_name ? Symbol(name_stem, disambig_idx+i) : name_stem )
        append!(p.parameters, [(name, data_type, name_stem, disambig_idx+i)])
    end
end

"""
    _is_datatype(sym, p)

Utility function that checks whether `sym` is a datatype. It returns true if
`sym` is either an in-built datatype (say float, StaticArray etc.) or if it is a
template argument.
"""
function _is_datatype(sym, p)
    (sym in p.template_args) && return true
    typeof(sym) <: Expr && sym.head == :tuple && return false
    typeof(eval(sym)) <: DataType
end

"""
    get_name_stem(name_stem::Symbol, parameters)

Get the stem of a name for a parameters and then add a disambiguation index.
Underscore `_` used in place of name is defaulted to `p`.
"""
function get_name_stem(name_stem::Symbol, parameters)
    if name_stem == :_
        name_stem = :p
    end
    disambig_idx = highest_idx_used(name_stem, parameters)
    name_stem, disambig_idx, (name_stem == :p)
end

"""
    highest_idx_used(name_stem, params)

Find the highest disambiguation index that has been used thus far for a given
`name_stem`.
"""
function highest_idx_used(name_stem, params)
    idxs = [p[4] for p in params if name_stem == p[3]]
    length(idxs) == 0 && return 0
    maximum(idxs)
end

"""
    parse_param_multi_names(line, p)

Parse a line that defines parameters of the diffusion. The line must be in one
of the formats:
    (p_name1, p_name2, ...) --> (number_of_parameters, datatype)
    (p_name1, p_name2, ...) --> datatype
    (p_name1, p_name2, ...) --> (datatype1, datatype2, ...)
In the former two cases defines `number_of_parameters`-many parameters, with
names `p_name1`, `p_name2`, ... and of `datatype` type. In the last case the
datatypes differ from parameter to parameter.
"""
function parse_param_multi_names(line, p)
    names = line.args[1].args
    data_types = line.args[2]

    if _is_datatype(data_types, p)
        data_types = repeat([data_types], length(names))
    else
        @assert data_types.head == :tuple
        if _is_datatype(data_types.args[1], p)
            @assert length(data_types.args) == length(names)
            @assert all( map(x->_is_datatype(x, p), data_types.args) )
            data_types = data_types.args
        else
            @assert typeof(eval(data_types.args[1])) <: Number
            @assert length(data_types.args) == 2
            @assert eval(data_types.args[1]) == length(names)
            @assert _is_datatype(data_types.args[2], p)
            data_types = repeat([data_types.args[2]], length(names))
        end
    end

    current_param_names = [param[1] for param in p.parameters]
    for (n,d) in zip(names, data_types)
        @assert !(n in current_param_names)
        append!(p.parameters, [(n,d,n,1)])
    end
end

#------------------------------------------------------------------------------#
#
#                       For parsing: CONJUGATE UPDATES
#
#------------------------------------------------------------------------------#

"""
    parse_line!(::Val{:conjugate}, line, p)

Parse a line that defines conjugate updates.
"""
function parse_line!(::Val{:conjugate}, line, p)
    if typeof(line.args[1]) <: Expr
        @assert line.args[1].head == :call
        function_name = line.args[1].args[1]
        fns = Vector{Expr}(undef, 0)
        if _symbol_in(function_name, _PHI)
            t_sym = line.args[1].args[2]
            x_sym = line.args[1].args[3]
            @assert line.args[2].head == :tuple
            num_funs = length(line.args[2].args)
            for expr in line.args[2].args
                add_phi_function!(expr, fns, t_sym, x_sym, p)
            end
        elseif function_name == :nonhypo
            add_nonhypo_function!(fns, line, p)
        end
        append!(p.fns, fns)
    elseif _symbol_in(line.args[1], _NUMNONHYPO)
        num = line.args[2]
        @assert typeof(num) <: Integer
        p.extras[_NUMNONHYPO[1]] = num
    end
end

"""
    add_nonhypo_function!(fns, line, p)

Add a definition of a function `nonhypo` that takes a current diffusion state
`x` and returns the elliptic coordinates.
"""
function add_nonhypo_function!(fns, line, p)
    fn_def = Expr(:call,
        :nonhypo,
        Expr(:(::),
            :P,
            p.name
        ),
        line.args[1].args[2]
    )
    push!(fns, Expr(:(=), fn_def, line.args[2]))
end

"""
    add_phi_function!(expr::Expr, fns::Vector{Expr}, t::Symbol, x::Symbol, p)

Construct a `phi` function used for conjugate updates.
"""
function add_phi_function!(
        expr::Expr,
        fns::Vector{Expr},
        t::Symbol,
        x::Symbol,
        p
    )
    idx = length(fns)
    name = p.name
    phi_def = Expr(:call,
        _PHI[1],
        Expr(:(::),
            Val{idx}
        ),
        :t,
        :x,
        Expr(:(::), :P, p.name)
    )

    @assert expr.head == :tuple
    @assert length(expr.args) == p.dims[_PROCESS[1]]
    phi_body = copy(expr)
    for i in 1:p.dims[_PROCESS[1]]
        cleanup_param_names!(phi_body.args[i], t, x, map(x->x[1], p.parameters))
    end

    phi = Expr(:(=), phi_def, phi_body)
    append!(fns, [phi])
end

"""
    cleanup_param_names!(expr, t::Symbol, x::Symbol, params)

A hepler function for cleaning up the names of variables in the definition of a
`phi` function.
"""
function cleanup_param_names!(expr, t::Symbol, x::Symbol, params)
    for (i,arg) in enumerate(expr.args)
        if typeof(arg) <: Symbol
            arg == t && (expr.args[i] = :t)
            arg == x && (expr.args[i] = :x)
            arg in params && (expr.args[i] = Expr(:., :P, QuoteNode(arg)))
        elseif typeof(arg) == Expr
            cleanup_param_names!(arg, t, x, params)
        end
    end
end

cleanup_param_names!(num::Number, t::Symbol, x::Symbol, params) = num

#------------------------------------------------------------------------------#
#
#                       For parsing: ADDITIONAL INFO
#
#------------------------------------------------------------------------------#

"""
    parse_line!(::Val{:additional}, line, p)

Parse a line that defines additional information about a diffusion process.
"""
function parse_line!(::Val{:additional}, line, p)
    name = line.args[1]
    _symbol_in(name, _STATESPACE) && (name = _STATESPACE[1])
    _symbol_in(name, _CONSTDIFF) && (name = _CONSTDIFF[1])
    _symbol_in(name, _LINEAR) && (name = _LINEAR[1])
    _symbol_in(name, _ELTYPE) && (name = _ELTYPE[1])

    p.extras[name] = line.args[2]
end

"""
    add_constdiff_function!(fns, p)

Add a definition of a function `consdiff` that indicates if the diffusion
coefficient is constant
"""
function add_constdiff_function!(fns, p)
    fn_def = Expr(:call,
        :constdiff,
        Expr(:(::),
            :P,
            p.name
        )
    )
    push!(fns, Expr(:(=), fn_def, p.extras[_CONSTDIFF[1]]))
end

#------------------------------------------------------------------------------#
#
#                   For parsing: DIFFUSION DIMENSION
#
#------------------------------------------------------------------------------#

"""
    parse_line!(::Val{:dimensions}, line, p)

Parse a line that defines the dimension of a diffusion process and the driving
Brownian motion.
"""
function parse_line!(::Val{:dimensions}, line, p)
    dim_of_what, dim = line.args[1], line.args[2]
    _symbol_in(dim_of_what, _WIENER) && (p.dims[_WIENER[1]] = dim)
    _symbol_in(dim_of_what, _PROCESS) && (p.dims[_PROCESS[1]] = dim)
end

#------------------------------------------------------------------------------#
#
#                            final routines
#
#------------------------------------------------------------------------------#

"""
    fill_unspecified_with_defaults(p)

Fill all unspecified variables with default values
"""
function fill_unspecified_with_defaults(p)
    fill_unspecified_with_defaults(Val{:dimensions}(), p)
    fill_unspecified_with_defaults(Val{:additional}(), p)
end

"""
    fill_unspecified_with_defaults(::Val{:dimensions}, p)

If unspecified, the dimension of the stochastic process and the driving Brownian
motion is set to 1.
"""
function fill_unspecified_with_defaults(::Val{:dimensions}, p)
    for key in [_WIENER[1], _PROCESS[1]]
        !haskey(p.dims, key) && (p.dims[key] = 1)
    end
end

"""
    fill_unspecified_with_defaults(::Val{:additional}, p)

If unspecified, there are no restriction on a state space, the volatility
coefficient is assumed constant, the diffusion is not linear and the datatype
of each coordinate is set to `Float64`.
"""
function fill_unspecified_with_defaults(::Val{:additional}, p)
    defaults = [
        (_STATESPACE[1], UnboundedStateSpace()),
        (_CONSTDIFF[1], false),
        (_LINEAR[1], false),
        (_ELTYPE[1], Float64),
        (_NUMNONHYPO[1], (
            haskey(p.extras, _STATESPACE[1]) ?
            p.extras[_STATESPACE[1]] :
            1
        )),
    ]
    for (key, default_val) in defaults
        !haskey(p.extras, key) && (p.extras[key] = default_val)
    end
end

"""
    prepare_abstract_type(stem, dims, data_type, state_restr)

Create a string defining a parent, abstract type from its `stem`, the dimensions
`dims` of the process and the driving Brownian motion, the datatype `data_type`
of each coordinate and the restrictions on the state space `state_restr`.
"""
function prepare_abstract_type(stem, dims, data_type, state_restr)
    Expr(:curly,
        stem,
        data_type,
        dims[_PROCESS[1]],
        dims[_WIENER[1]],
        state_restr,
    )
end

"""
    createstruct(abstract_type, name, params)

Create code that defines a struct defining a diffusion process.
"""
function createstruct(abstract_type, p)
    param_vec = map(p.parameters) do (name, data_type, _, _)
        Expr(:(::), name, data_type)
    end

    _new = (
        p.template_args == Any[] ?
        :new :
        Expr(:curly, :new, p.template_args...,)
    )

    constructor_def = Expr(:call, p.name, param_vec...,)

    constructor_def = (
        p.template_args == Any[] ?
        constructor_def :
        Expr(:where,
            constructor_def,
            p.template_args...,
        )
    )

    Expr(:struct,
        false, # immutable
        Expr(:<:,
            p.curly_name,
            abstract_type,
        ), # first line defining the struct
        Expr(:block,
            param_vec..., # parameters
            Expr(:function, # constructor
                constructor_def,
                Expr(:block,
                    Expr(:call,
                        _new,
                        map(x->x[1], p.parameters)...,
                    ),
                ), # body of the constructor
            ), # constructor
        ), # body defining the struct
    )
end


text:
using Distributions,StatsPlots, ForwardDiff, LinearAlgebra, PositiveFactorizations, Optim

n_samples = 2000

n = 100

x = randn(n,50)
theta = [repeat([1],10); repeat([0],40)]
y = x*theta + randn(n)

prior = Laplace(0,2)

ll(theta) = logpdf(MvNormal(I(n)),y-x*theta)

post(theta) = ll(theta) + sum(logpdf.(prior,theta))

init = rand(prior,50)

grad(theta) = ForwardDiff.gradient(post,theta)

ot = optimize(x->-1*ll(x),zeros(50),BFGS(), autodiff=:forward)

M = grad(ot.minimizer)*grad(ot.minimizer)'#mass matrix
M = inv(diagm(diag(M)))

mom_dist = MvNormal(M)

tol_psi = 1e18

#######
### start here
#

pars = zeros(n_samples,50)

uni = Uniform(0,1)

leap_step = 0.09

L = max(Int(ceil(1/leap_step)),20)

window_size = 10

window_dist = DiscreteUniform(1,window_size)

accept = 0

old_par = init#repeat([0],50)

leap_dist = Normal(leap_step,0.05)

for j = 1:500
    psi = rand(mom_dist)
    old_psi = psi
    new_par = old_par
    leap_step = rand(leap_dist)
    s = rand(window_dist)
    forward_leap = zeros(50,2)
    sum_p = 0
    for i = 1:(s*L) #forward leapfrog
        psi = psi + 1/2*leap_step*grad(new_par)
        new_par = new_par + inv(M)*leap_step*psi
        psi = psi + 1/2*leap_step*grad(new_par)
        pp = exp(logpdf(mom_dist,psi)+post(new_par))
        sum_p = sum_p + pp
        if pp/sum_p > rand(uni)
            forward_leap = [psi new_par]
        end
    end
    backward_leap = zeros(50,2)
    sum_p = 0
    psi = old_psi
    new_par = old_par
    for i = 1:((window_size-s)*L) #backward leapfrog
        psi = psi - 1/2*leap_step*grad(new_par)
        new_par = new_par - inv(M)*leap_step*psi
        psi = psi - 1/2*leap_step*grad(new_par)
        pp = exp(logpdf(mom_dist,psi)+post(new_par))
        sum_p = sum_p + pp
        if pp/sum_p > rand(uni)
            backward_leap = [psi new_par]
        end
    end
    #println(psi)
    if any(psi .> tol_psi)
        @info "Energy diverging. Quitting"
        break
    end
    r = post(forward_leap[:,2]) + logpdf(mom_dist,forward_leap[:,1]) - (post(backward_leap[:,2])+ logpdf(mom_dist,backward_leap[:,1]))
    a = rand(uni)
    r = min(1,exp(r))
    if r > a
        pars[j,:] = new_par
        global old_par = new_par
        global accept = accept + 1
    else
        pars[j,:] = old_par
    end
    println("Iteration ", j, " Accept ratio ", accept/j*100, "%")
end

leap_step = 0.15

L = max(Int(ceil(1/leap_step)),20)

accept = 0

old_par = reshape(mapslices(mode,pars[1:500,:],dims = 1),50)

for j = 1:2000
    psi = rand(mom_dist)
    old_psi = psi
    new_par = old_par
    leap_step = rand(leap_dist)
    s = rand(window_dist)
    forward_leap = zeros(50,2)
    sum_p = 0
    for i = 1:(s*L) #forward leapfrog
        psi = psi + 1/2*leap_step*grad(new_par)
        new_par = new_par + inv(M)*leap_step*psi
        psi = psi + 1/2*leap_step*grad(new_par)
        pp = exp(logpdf(mom_dist,psi)+post(new_par))
        sum_p = sum_p + pp
        if pp/sum_p > rand(uni)
            forward_leap = [psi new_par]
        end
    end
    backward_leap = zeros(50,2)
    sum_p = 0
    psi = old_psi
    new_par = old_par
    for i = 1:((window_size-s)*L) #backward leapfrog
        psi = psi - 1/2*leap_step*grad(new_par)
        new_par = new_par - inv(M)*leap_step*psi
        psi = psi - 1/2*leap_step*grad(new_par)
        pp = exp(logpdf(mom_dist,psi)+post(new_par))
        sum_p = sum_p + pp
        if pp/sum_p > rand(uni)
            backward_leap = [psi new_par]
        end
    end
    #println(psi)
    if any(psi .> tol_psi)
        @info "Energy diverging. Quitting"
        break
    end
    r = post(forward_leap[:,2]) + logpdf(mom_dist,forward_leap[:,1]) - (post(backward_leap[:,2])+ logpdf(mom_dist,backward_leap[:,1]))
    a = rand(uni)
    r = min(1,exp(r))
    if r > a
        pars[j,:] = new_par
        global old_par = new_par
        global accept = accept + 1
    else
        pars[j,:] = old_par
    end
    println("Iteration ", j, " Accept ratio ", accept/j*100, "%")
end

StatsPlots.density(pars[200:1700,1])

mm = mapslices(mode,pars[500:1126,:],dims = 1)

scatter(mm')

StatsPlots.density(pars[200:1700,1])
StatsPlots.density(pars[200:1700,2])
StatsPlots.density(pars[200:1700,3])
StatsPlots.density(pars[200:1700,4])
StatsPlots.density(pars[200:1700,5])
StatsPlots.density(pars[200:1700,6])
StatsPlots.density(pars[200:1700,7])
StatsPlots.density(pars[200:1700,8])
StatsPlots.density(pars[200:1700,9])
StatsPlots.density(pars[200:1700,10])
StatsPlots.density(pars[200:1700,11])
StatsPlots.density(pars[200:1700,12])
StatsPlots.density(pars[200:1700,13])
StatsPlots.density(pars[200:1700,14])
StatsPlots.density(pars[200:1700,15])
StatsPlots.density(pars[200:1700,16])

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
using Base.Test
using NGSIM
using NBInclude

# nbinclude(joinpath(dirname(@__FILE__), "..", "jnotebooks", "Demo.ipynb"))


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export l1_distance, l2_distance, squared_l2_distance, linf_distance, lp_distance

"""
    l1_distance(a, b)::Float64

Computes the Manhattan's distance between `a` and `b`
"""
function l1_distance(a, b)::Float64
    d::Float64 = 0.0 #zero(eltype(a))

    @inbounds @simd for i = 1:length(a)
	    m = a[i] - b[i]
        d += ifelse(m > 0, m, -m)
    end

    d
end

"""
    l2_distance(a, b)::Float64

Computes the Euclidean's distance betweem `a` and `b`
"""
function l2_distance(a, b)::Float64
    #d = zero(eltype(a))
    d::Float64 = 0.0

    @inbounds @simd for i = 1:length(a)
        m = a[i] - b[i]
        d += m * m
    end

    sqrt(d)
end

"""
    squared_l2_distance(a, b)::Float64

Computes the squared Euclidean's distance between `a` and `b`
"""
function squared_l2_distance(a, b)::Float64
    # d = zero(eltype(a))
    d::Float64 = 0.0

    @inbounds @simd for i = 1:length(a)
        m = a[i] - b[i]
        d += m * m
    end

    d
end

"""
    linf_distance(a, b)::Float64

Computes the max or Chebyshev'se distance
"""
function linf_distance(a, b)::Float64
   d::Float64 = 0.0 # d = zero(eltype(a))

    @inbounds @simd for i = 1:length(a)
        m = abs(a[i] - b[i])
        d = max(d, m)
    end

    d
end

"""
    lp_distance(p_::Real)

Creates a function that computes computes generic Minkowski's distance with the given `p_`
"""
function lp_distance(p::Real)
    p = convert(Float64, p)
    invp = 1.0 / p

    function _lp(a, b)::Float64
        d::Float64 = 0.0 # d = zero(eltype(a))

        @inbounds @simd for i = 1:length(a)
            m = abs(a[i] - b[i])
            d += m ^ p
        end

        d ^ invp
    end
end


text:
using PyPlot
pygui(true)

fig, ax = subplots()
xlim(0.0, 0.005)
ylim(0.0, 0.005)
grid()

circle1 = matplotlib[:patches][:Circle]((0.001, 0.001), radius = 0.001, facecolor="b")
circle2 = matplotlib[:patches][:Circle]((0.003, 0.003), radius = 0.001, facecolor="r")
p = matplotlib[:collections][:PatchCollection]([circle1, circle2])
ax[:add_collection](p)
show()


text:
# Particle: variant R2 (fully rational case; leja + repeated nodes + freeze)

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test

include("particle_test_utils.jl")

@bench @testset "NLEIGS: Particle variant R2" begin
    verbose = displaylevel

    nep, Σ, Ξ, v, nodes, xmin, xmax = particle_init(2)

    # solve nlep
    lambda, X, res, solution_info = nleigs(nep, Σ, Ξ=Ξ, logger=verbose > 0 ? 1 : 0, maxdgr=50, minit=30, maxit=100, v=v, nodes=nodes, return_details=verbose > 1)

    verify_lambdas(2, nep, lambda, X)

    if verbose > 1
        include("nleigs_residual_plot.jl")
        approx_Σ = [xmin-im*1e-10, xmin+im*1e-10, xmax+im*1e-10, xmax-im*1e-10]
        nleigs_residual_plot("Particle: variant R2", solution_info, approx_Σ)
    end
end