text:
"""
	experiment(score_fun, parameters, data, savepath; save_entries...)
Eval score function on test/val/train data and save.
"""
function experiment(score_fun, parameters, data, savepath; verb=true, save_result=true, save_entries...)
	tr_data, val_data, tst_data = data

	# extract scores
	tr_scores, tr_eval_t, _, _, _ = @timed score_fun(tr_data[1])
	val_scores, val_eval_t, _, _, _ = @timed score_fun(val_data[1])
	tst_scores, tst_eval_t, _, _, _ = @timed score_fun(tst_data[1])

	# now save the stuff
	savef = joinpath(savepath, savename(parameters, "bson", digits=5))
	result = (
		parameters = parameters,
		tr_scores = tr_scores,
		tr_labels = tr_data[2],
		tr_eval_t = tr_eval_t,
		val_scores = val_scores,
		val_labels = val_data[2],
		val_eval_t = val_eval_t,
		tst_scores = tst_scores,
		tst_labels = tst_data[2],
		tst_eval_t = tst_eval_t
		)
	result = Dict{Symbol, Any}([sym=>val for (sym,val) in pairs(merge(result, save_entries))]) # this has to be a Dict
	if save_result
		tagsave(savef, result, safe = true)
		verb ? (@info "Results saved to $savef") : nothing
	end
	result
end

"""
	edit_params(data, parameters)
This modifies parameters according to data. Default version only returns the input arg.
Overload for models where this is needed.
"""
function edit_params(data, parameters)
	parameters
end

"""
	check_params(savepath, parameters)

Returns `true` if the model with given parameters wasn't already trained and saved.
"""
function check_params(savepath, parameters)
	if ~isdir(savepath)
		return true
	end
	# filter out duplicates created by tagsave
	fs = filter(x->!(occursin("_#", x)), readdir(savepath))
	# filter out model files
	fs = filter(x->!(startswith(x, "model")), fs)
	# if the first argument name contains a "_", than the savename is parsed wrongly
	saved_params = map(x -> DrWatson.parse_savename("_"*x)[2], fs)
	# now filter out saved models where parameter names are different or missing
	pkeys = collect(keys(parameters))
	filter!(ps->intersect(pkeys, Symbol.(collect(keys(ps))))==pkeys, saved_params)
	for params in saved_params
		all(map(k->params[String(k)] == parameters[k], pkeys)) ? (return false) : nothing
	end
	true
end


text:
using Images, ImageView, ImageSegmentation

# convert image segments to mean color value
segment_to_image(segments) = map(i->segment_mean(segments, i), labels_map(segments))

# one object
img = load("sample-images/cat-3352842_640.jpg");
imshow(img);

# find segments
segments_1 = felzenszwalb(img, 75)
segments_2 = felzenszwalb(img, 75, 150)
segments_3 = felzenszwalb(img, 75, 350)

# define helper values
img_width = size(img, 2)

new_img = hcat(
    segment_to_image(segments_1),
    segment_to_image(segments_2),
    segment_to_image(segments_3)
)

new_img[:, img_width] = new_img[:, img_width*2] = colorant"black"

imshow(new_img)

# multiple objects
img = load("sample-images/bird-3183441_640.jpg");
#imshow(img);

# find segments
segments_1 = felzenszwalb(img, 10)
segments_2 = felzenszwalb(img, 30, 50)
segments_3 = felzenszwalb(img, 35, 300)

# create a new image which will be a stack
# of 3 different results of felzenszwalb algorithm
img_width = size(img, 2)

new_img = hcat(
    segment_to_image(segments_1),
    segment_to_image(segments_2),
    segment_to_image(segments_3)
)

new_img[:, img_width] = new_img[:, img_width*2] = colorant"black"

imshow(new_img)

ext: jl


text:
module RocketSyncActorTest

using Test
using Rocket

@testset "SyncActor" begin

    println("Testing: actor SyncActor")

    @testset begin
        actor  = KeepActor{Int}()
        synced = SyncActor{Int, KeepActor{Int}}(actor)

        source = interval(1) |> take(5)

        subscribe!(source, synced)

        wait(synced)

        @test actor.values == [ 0, 1, 2, 3, 4 ]
    end

    @testset begin
        @test sync(void(Int)) isa SyncActor{Int, VoidActor{Int}}
    end

    @testset begin
        values = Int[]

        factory  = lambda(on_next = (d) -> push!(values, d))
        synced   = sync(factory)

        subscribe!(interval(1) |> take(5), synced)

        wait(synced)

        @test values == [ 0, 1, 2, 3, 4 ]
    end

    @testset begin
        completions = []

        factory  = lambda(on_complete = () -> push!(completions, 1))
        synced   = sync(factory)

        subscribe!(completed(), synced)

        wait(synced)

        @test completions == [ 1 ]
    end

    @testset begin
        errors = []

        factory  = lambda(on_error = (d) -> push!(errors, d))
        synced   = sync(factory)

        subscribe!(faulted("e"), synced)

        wait(synced)

        @test errors == [ "e" ]
    end

    @testset begin
        source = never(Int)
        actor  = sync(void(Int), timeout = 100)

        subscribe!(source, actor)

        @test_throws SyncActorTimedOutException wait(actor)
    end

    struct DummyActor end

    @testset begin
        @test_throws InvalidActorTraitUsageError sync(DummyActor())
    end
end

end


text: print("Julia >>")

max_stars_repo_name: Akuramu/HelloWorld

max_issues_repo_name: Akuramu/HelloWorld

max_forks_repo_name: Akuramu/HelloWorld


text:
"""
`Data` groups the instance-parsing functions and the types involved.

The parsing is extensible, as users can define any new types and extend
the parsing functions to the types. It is advised to always return the
same type for the same format type given as input, to avoid type-unstability
and have a clear mapping.
However, "auto" formats can also be created, to simplify the user's life
(in special during testing) paying for a little overhead (and risk
of misdetection).

The parsing function is extended for parametric types (e.g., `Classic_G2KP{D,
S, P}`) so the data type returned (e.g., `G2KP{D, S, P}`) has the demand/N/M
(`D`), length/width (`S`), and profit/area (`P`) fields with the right types.
The function is also extend to `Val` types for convenience.  They often map to
a parametric type with sensible safe defaults (e.g.  `Val{:Classic_G2KP}` ->
`Classic_G2KP{Int, Int, Int}`). But a finer-grained control is also possible,
as you can define something like (`Val{:Classic_G2KP_16}` ->
`Classic_G2KP_16{Int16, Int16, Int16}`). `Val` parameters are also important
because module `CommandLine` currently takes `String`s that are passed to `(Val
∘ Symbol)` in order to define the instance format to be used.

This module already extends the instance reading for many instance formats.
From the 2DCPackGen (further abbreviated here as `CPG`, Cutting and Packing
Generator, because `struct`s cannot start with a digit), the module supports:
`SLOPP`, `MHLOPPW`, `ODPW`, and `SSSCSP`. They all follow the same pattern:
`SLOPP` (the `struct` type with the data), `CPG_SLOPP{D, S, P}` (for the
clearly typed format type), and `Val{:CPG_SLOPP}` for the safe defaults (i.e.,
the same as `CPG{Int, Int, Int}`, note this will be `32` in some machined and
`64` in others). For details on these formats look at the 2DCPackGen tool.

Some file formats aggregate multiple instances in a single file. To allow
generic functions to deal with this possibility, it is recommended to
specialize [`GuillotineModels.Data.is_collection`](@ref). This function should
be specialized for the format type (not the return of the parsing).
A format that does not specialize it is assumed to return a single instance
upon parsing (i.e., [`GuillotineModels.Data.read_from_string`](@ref)).

The `Classic_G2KP`/`G2KP` format/data is also provided in the same fashion.
The format is described in the method documentation.

The function may `@warn` the users of "strange" (but not incorrect)
structure of the given instance, and `error` if the format is incorrect.

The current methods dispatch on the type of the first object, even if all
current format types have empty objects. Consequently, instantiations like
`CPG_SLOPP{Int8, Int16, Int32}()` and `Val(:CPG_SLOPP)` are expected,
instead of just the type itself. In the future may be possible to
add parsing options/configuration inside the format objects.
"""
module Data

import TimerOutputs.@timeit
import ..TIMER
import ..solution_value
using AutoHashEquals: @auto_hash_equals
using UnPack: @unpack

export read_from_string, read_from_file, is_collection
export GenericParseError

# From classic_g2kp
export Classic_G2KP, G2KP
include("classic_g2kp.jl")

# From cpg
export CPG_SLOPP, CPG_MHLOPPW, CPG_ODPW, CPG_SSSCSP
export SLOPP, MHLOPPW, ODPW, SSSCSP
include("cpg.jl")

# TODO: consider the assymetry that `read_from_file`/`read_from_string`
# has the real code inside `read_from_string` but `write_to_file`
# has a `filepath` and a `io` versions and it is the io version that do
# the real work.
"""
    read_from_file(format, filepath :: AbstractString)

Reads the given file and apply `read_from_string` to its contents.

The (already provided) generic body of this function should be enough for most
formats without the need of explicit specialization.

Note: this function remove all '\\r' characters from the obtained `String` to
avoid 'dos' vs 'unix' format inconsistencies.

See also: [`read_from_string`](@ref)
"""
@timeit TIMER function read_from_file(
	format, filepath :: AbstractString
)
	return open(filepath) do f
		read_from_string(format, replace(read(f, String), '\r' => ""))
	end
end

"""
    write_to_file(format, instance, filepath :: AbstractString[, mode])
    write_to_file(format, instance, io)

Wite the `instance` in `filepath`/`io` using `format`.

NOTE: this function is not yet implemented for most format and instance
types.

`mode` is passed to `open`, by default it is: "w+".

New `format`/`instance` combinations should implement a method taking an
`io::IO` third parameter (otherwise ambiguity problems may arise).  The generic
method taking a `filename` just open the file and call the `io`-like version
over it.

Note: the `write_to_string` methods implemented by this package does not
use '\\r\\n' (carriage return line endings) just '\\n'.

To write to a `String`, pass an empty `IOBuffer` object and after the call
execute `read(seekstart(iob), String)` to extract the `String` (this is the
most efficient way).

See also: [`read_from_file`, `read_from_string`](@ref)
"""
@timeit TIMER function write_to_file(
	format, instance, filepath :: AbstractString, mode :: AbstractString = "w+"
)
	return open(filepath, mode) do io
		write_to_file(format, instance, io)
	end
end

"""
    is_collection(format) :: Bool

Indicate if a format is parsed to a single instance or to a collection.

If the format may store a variable number of instances, then passing the
format object to this function should return `true`; if the format always
store a single instance, then it should return false.

If this function is not specialized it assumes the format parses to a single
instance.

If the function returns `true` for some format, it is assumed that the
value returned by [`read_from_string`](@ref) (and its file counterpart)
respect the `Base.iterate` interface; if a file in this format has just a
single instance, this single instance should be wrapped in the same container
used when there are multiple instances (any overhead caused by wrapping will
probably not be larger than the overhead caused by type unstability).

Ideally, the value returned should not depend on the value of the
argument, but only the type of the argument.
"""
is_collection(format) = false

"""
    abstract type ParseError <: Exception end

Supertype for parse errors.
"""
abstract type ParseError <: Exception end

"""
    GenericParseError(format, error_message)

Most general error thrown by `parse_from_string` when it fails.

Ideally, methods extending `parse_from_string` should use a subtype of
`ParseError`, either custom or this one.
"""
struct GenericParseError <: ParseError
	format :: Any
	error_message :: String
end

function Base.showerror(io :: IO, e :: GenericParseError)
	print(io,
		"Trying to parse $(e.format) has resulted in the following error: " *
		e.error_message
	)
end

end # module


text:
immutable FLSynapseParameter
end

@withkw type FLSynapse
    param::FLSynapseParameter = FLSynapseParameter()
    W::Matrix{Float}  # synaptic weight
    rI::Vector{Float} # postsynaptic rate
    rJ::Vector{Float} # presynaptic rate
    g::Vector{Float}  # postsynaptic conductance
    P::Matrix{Float}  # <rᵢrⱼ>⁻¹
    q::Vector{Float}  # P * r
    u::Vector{Float} # force weight
    w::Vector{Float} # output weight
    f::Float = 0 # postsynaptic traget
    z::Float = 0.5randn()  # output z ≈ f
    records::Dict = Dict()
end

function FLSynapse(pre, post; σ = 1.5, p = 0.0, α = 1)
    rI, rJ, g = post.r, pre.r, post.g
    W = σ * 1 / √pre.N * randn(post.N, pre.N) # normalized recurrent weight
    w = 1 / √post.N * (2rand(post.N) - 1) # initial output weight
    u = 2rand(post.N) - 1 # initial force weight
    P = α * eye(post.N) # initial inverse of C = <rr'>
    q = zeros(post.N)
    FLSynapse(;@symdict(W, rI, rJ, g, P, q, u, w)...)
end

@replace function forward!(c::FLSynapse, param::FLSynapseParameter)
    z = dot(w, rI)
    BLAS.A_mul_B!(q, P, rJ)
    BLAS.A_mul_B!(g, W, rJ)
    BLAS.axpy!(z, u, g)
end

@replace function plasticity!(c::FLSynapse, param::FLSynapseParameter, dt::Float, t::Float)
    C = 1 / (1 + dot(q, rI))
    BLAS.axpy!(C * (f - z), q, w)
    BLAS.ger!(-C, q, q, P)
end


text:
# See benchmark.jl for definitions

# This file collects data and tests functions.
# The intention is preparing input for a benchmark suite,
# but the output from this file actually is sufficient for some purposes.
#
# Both log file, results plots and tables are printed to the same file.
# Viewing the plots in a text editor is probably possible, see UnicodePlots.jl.
# Pull request welcome if someone can figure out how to do it.
# Running this file on a Windows laptop with all browsers takes 5-10 minutes.

if !isdefined(:SRCPATH)
    const SRCPATH = Base.source_dir() == nothing ? Pkg.dir("WebSockets", "benchmark") :Base.source_dir()
    const LOGGINGPATH = realpath(joinpath(SRCPATH, "../logutils/"))
    SRCPATH ∉ LOAD_PATH && push!(LOAD_PATH, SRCPATH)
    LOGGINGPATH ∉ LOAD_PATH && push!(LOAD_PATH, LOGGINGPATH)
    include(joinpath(SRCPATH, "functions_open_browsers.jl"))
end

"A vector of message sizes"
const VSIZE = reverse([i^3 * 1020 for i = 1:12])
include(joinpath(SRCPATH, "functions_open_browsers.jl"))
include(joinpath(SRCPATH, "functions_benchmark.jl"))
#
prepareworker()
# Load modules on both processes
import HTTP
using WebSockets
using ws_jce
using UnicodePlots
import IndexedTables.table
import ws_hts: listen_hts, getws_hts
#
remotecall_fetch(ws_jce.clog, 2, "ws_jce ", :green, " is ready")
# Start async HTS server on this process and check that it is up and running
const TIMEOUT = Base.Dates.Second(20)
hts_task = start_hts(TIMEOUT)

"""
Prepare logging for this process

    logs are written to console as well as to a file, which
    is different for the other process.
    Logs in other processes appear in console with a delay, hence the timestamp
    is interesting.
    To write log file buffer to disk immediately, call zflush().
    To drop duplicate console log, use zlog(..) instead of clog(..)
"""
const ID = "Benchmark"
const LOGFILE = "benchmark_prepare.log"
import logutils_ws: logto, clog, zlog, zflush, clog_notime
fbm = open(joinpath(SRCPATH, "logs", LOGFILE), "w")
logto(fbm)
clog(ID, "Started async HTS and prepared parallel worker")
zflush()

#
#  Do an initial test run, HTS-JCE, brief text output
#

const INITSIZE = 130560
const INITN = 200

# Measured time interval vectors [ns]
# Time measurements are taken directly both at server and client
(testid, serverlatencies, clientlatencies) = HTS_JCE(INITN, INITSIZE)
# Calculate speeds [ns/b] and averaged speeds (bandwidths)
(serverspeeds, clientspeeds,
    serverbandwidth, clientbandwidth) = serverandclientspeeds(INITSIZE, serverlatencies, clientlatencies)
# Store plots and a table in dictionaries
vars= [:serverspeeds, :clientspeeds];
init_plots = Dict(testid => lp(vars, testid));
init_tables = Dict(testid => table(eval.(vars)..., names = vars));
init_serverbandwidths = Dict(testid => serverbandwidth);
init_clientbandwidths = Dict(testid => clientbandwidth);
# Sleep to avoid interspersing with worker output to REPL
sleep(2)
# Brief output to file and console
clog(testid, " Initial test run with messagesize ", INITSIZE, " bytes \n\t",
    "serverbandwidth = ", :yellow,  round(serverbandwidth, 4), :normal, " [ns/b] = [s/GB]\n\t",
    :normal, "clientbandwidth = ", :yellow, round(clientbandwidth, 4), :normal, " [ns/b] = [s/GB]")

#
#  Continue initial test run with HTS-BCE, brief text output for each browser
#

COUNTBROWSER.value = 0
serverbandwidth = 0.
clientbandwidth = 0.
serverbandwidths = Vector{Float64}()
clientbandwidths = Vector{Float64}()
t1 = Vector{Int}()
t2 = Vector{Int}()
t3 = Vector{Int}()
t4 = Vector{Int}()
browser = ""
success = true
while success
    # Measured time interval vectors [ns] for the next browser in line
    # Time measurements are taken only at server; a message pattern
    # is used to distinguish server and client performance
    (browser, t1, t2, t3, t4) = HTS_BCE(INITN, INITSIZE)
    if browser != ""
        # Calculate speeds [ns/b] and averaged speeds (bandwidths)
        (serverspeeds, clientspeeds,
            serverbandwidth, clientbandwidth) =
            serverandclientspeeds_indirect(INITSIZE, t1, t2, t3, t4)
        # Store plots and a table in dictionaries
        testid = "HTS_BCE " * browser
        push!(init_plots, testid => lp(vars, testid));
        push!(init_tables, testid => table(eval.(vars)..., names = vars));
        push!(init_serverbandwidths, testid => serverbandwidth);
        push!(init_clientbandwidths, testid => clientbandwidth);
        # Brief output to file and console
        clog(testid, " Initial test run with messagesize ", INITSIZE, " bytes \n\t",
            "serverbandwidth = ", :yellow,  round(serverbandwidth, 4), :normal, " [ns/b] = [s/GB]\n\t",
            :normal, "clientbandwidth = ", :yellow, round(clientbandwidth, 4), :normal, " [ns/b] = [s/GB]")
    else
        success = false
    end
end

#
#  Collect HTS_JCE bandwidths vs message size
#

const SAMPLES = 200
serverbandwidths = Dict{String, Vector{Float64}}()
clientbandwidths = Dict{String, Vector{Float64}}()
for msgsiz in VSIZE
    # Measured time interval vectors [ns]
    # Time measurements are taken directly both at server and client
    (testid, serverlatencies, clientlatencies) = HTS_JCE(SAMPLES, msgsiz)
    # Find averaged speed (bandwidth) scalars
    (_, _,
        sbw, cbw) = serverandclientspeeds(msgsiz, serverlatencies, clientlatencies)
    # Assign storage
    !haskey(serverbandwidths, testid) && push!(serverbandwidths, testid => Vector{Float64}())
    !haskey(clientbandwidths, testid) && push!(clientbandwidths, testid => Vector{Float64}())
    # Store bandwidths from this size
    push!(serverbandwidths[testid], sbw)
    push!(clientbandwidths[testid], cbw)
end

#
#  Collect HTS_BCE bandwidths
#  For individual and available browsers
#  This opens a lot of browser tabs; there may
#  be no easy way of closing them except
#  by closing Julia. internet explorer even
#  opens a separate window every time.
#
for msgsiz in VSIZE
    COUNTBROWSER.value = 0
    t1 = Vector{Int}()
    t2 = Vector{Int}()
    t3 = Vector{Int}()
    t4 = Vector{Int}()
    browser = ""
    success = true
    while success
        # Measured time interval vectors [ns] for the next browser in line
        # Time measurements are taken only at server; a message pattern
        # is used to distinguish server and client performance
        (browser, t1, t2, t3, t4) = HTS_BCE(SAMPLES, msgsiz)
        testid = "HTS_BCE " * browser
        if browser != ""
            # Find averaged speed (bandwidth) scalars
            (_, _,
                sbw, cbw) =
                    serverandclientspeeds_indirect(msgsiz, t1, t2, t3, t4)
            # Assign storage
            !haskey(serverbandwidths, testid) && push!(serverbandwidths, testid => Vector{Float64}())
            !haskey(clientbandwidths, testid) && push!(clientbandwidths, testid => Vector{Float64}())
            # Store bandwidths from this size
            push!(serverbandwidths[testid], sbw)
            push!(clientbandwidths[testid], cbw)
        else
            success = false
        end
    end
end

#
#   Measurements are done. Close server and log file, open results log file.
#
ws_hts.close_hts()
clog(ID, "Closing HTS server")
const RESULTFILE = "benchmark_results.log"
clog(ID, "Results are summarized in ", joinpath(SRCPATH, "logs", RESULTFILE))
fbmr = open(joinpath(SRCPATH, "logs", RESULTFILE), "w")
logto(fbmr)
close(fbm)

#
# Find optimum message size and nominal 100% bandwidths
# Make and store plots and tables
#
test_bestserverbandwidths = Dict{String, Float64}()
test_bestclientbandwidths = Dict{String, Float64}()
test_bestserverlatencies = Dict{String, Float64}()
test_bestclientlatencies = Dict{String, Float64}()
test_plots = Dict()
test_tables = Dict()
test_latency_plots = Dict()
test_latency_tables = Dict()
serverbandwidth = Vector{Float64}()
clientbandwidth = Vector{Float64}()
serverlatency = Vector{Float64}()
clientlatency = Vector{Float64}()
bestserverbandwidth = 0.
bestclientbandwidth = 0.
bestserverlatency = 0.
bestclientlatency = 0.

for testid in keys(serverbandwidths)
    serverbandwidth = serverbandwidths[testid]
    clientbandwidth = clientbandwidths[testid]
    # Store the optimal bandwidths in a dictionary
    bestserverbandwidth = minimum(serverbandwidth)
    bestclientbandwidth = minimum(clientbandwidth)
    push!(test_bestserverbandwidths, testid => bestserverbandwidth);
    push!(test_bestclientbandwidths,  testid => bestclientbandwidth);
    # Store msgsiz-bandwidth line plots and tables in dictionaries
    vars = [:serverbandwidth, :clientbandwidth]
    tvars = vcat([:VSIZE], vars)
    push!(test_plots,  testid => lp([:VSIZE, :VSIZE], vars, testid));
    push!(test_tables,  testid => table(eval.(tvars)..., names = tvars));
    # Store msgsiz-latency line plots and tables in dictionaries
    serverlatency = serverbandwidth .* VSIZE
    clientlatency = clientbandwidth .* VSIZE
    bestserverlatency = minimum(serverlatency)
    bestclientlatency = minimum(clientlatency)
    push!(test_bestserverlatencies,  testid => bestserverlatency);
    push!(test_bestclientlatencies,  testid => bestclientlatency);
    vars = [:serverlatency, :clientlatency]
    tvars = vcat([:VSIZE], vars)
    push!(test_latency_plots,  testid => lp([:VSIZE, :VSIZE], vars, testid));
    push!(test_latency_tables,  testid => table(eval.(tvars)..., names = tvars));

    # Brief output to file and console
    clog_notime(testid, :normal, " Varying message size: \n\t",
        "bestserverbandwidth = ", :yellow,  round(bestserverbandwidth, 4), :normal, " [ns/b] = [s/GB]",
        " @ size = ", VSIZE[findfirst(serverbandwidth, bestserverbandwidth)], " b\n\t",
        :normal, "bestclientbandwidth = ", :yellow, round(bestclientbandwidth, 4), :normal, " [ns/b] = [s/GB]",
        " @ size = ", VSIZE[findfirst(clientbandwidth, bestclientbandwidth)], " b\n\t",
        "bestserverlatency = ", :yellow,  Int(round(bestserverlatency)), :normal, " [ns] ",
        " @ size = ", VSIZE[findfirst(serverlatency, bestserverlatency)], " b\n\t",
        :normal, "bestclientlatency = ", :yellow, Int(round(bestclientlatency)), :normal, " [ns]",
        " @ size = ", VSIZE[findfirst(clientlatency, bestclientlatency)], " b\n\t"
        )
end

#
#   Full text output
#   The plots are not currently readably encoded in the text file
#

clog_notime(ID, :bold, :yellow, " Plots of all samples :init_plots [ns/b], message size ", INITSIZE, " b ", SAMPLES, " samples" )
foreach(values(init_plots)) do pls
    foreach(values(pls)) do pl
        clog_notime(pl)
    end
end
clog_notime(ID, :bold, :yellow, " Tables of all samples, :init_tables, message size ", INITSIZE, " b ", SAMPLES, " samples" )
for (ke, ta) in  init_tables
    clog_notime(ke, "\n=> ", ta, "\n")
end

clog_notime(ID, :bold, :yellow, " Plots of varying size messages :test_plots [ns/b],\n\t VSIZE = ", VSIZE)
foreach(values(test_plots)) do pls
    foreach(values(pls)) do pl
        clog_notime(pl)
    end
end

clog_notime(ID, :bold, :yellow, " Tables of varying size messages :test_tables [ns/b]")
for (ke, ta) in  test_tables
    clog_notime(ke, "\n=> ", ta, "\n")
end

clog_notime(ID, :bold, :yellow, " Plots of varying size messages :test_latency_plots [ns],\n\t VSIZE = ", VSIZE)
foreach(values(test_latency_plots)) do pls
    foreach(values(pls)) do pl
        clog_notime(pl)
    end
end

clog_notime(ID, :bold, :yellow, " Tables of varying size messages :test_latency_tables [ns]")
for (ke, ta) in  test_latency_tables
    clog_notime(ke, "\n=> ", ta, "\n")
end

clog_notime(ID, :bold, :yellow, " Dictionary  :test_bestserverlatencies [ns]")
for (ke, va) in  test_bestserverlatencies
    clog_notime(ke, " => \t", Int(round(va)))
end

clog_notime(ID, :bold, :yellow, " Dictionary  :test_bestclientlatencies [ns]")
for (ke, va) in  test_bestclientlatencies
    clog_notime(ke, " => \t", Int(round(va)))
end

clog_notime(ID, :bold, :yellow, " Dictionary  :test_bestserverbandwidths [ns/b]")
for (ke, va) in  test_bestserverbandwidths
    clog_notime(ke, " => \t", round(va, 4))
end

clog_notime(ID, :bold, :yellow, " Dictionary  :test_bestclientbandwidths [ns/b]")
for (ke, va) in  test_bestclientbandwidths
    clog_notime(ke, " => \t", round(va, 4))
end

zflush()


text:
@symbol_func function plasma_parameter(cur_reactor::AbstractReactor)
  cur_lambda = 5.4e6

  cur_lambda *= sqrt( cur_reactor.T_bar ^ 3 )

  cur_lambda /= sqrt( cur_reactor.n_bar )

  cur_lambda
end


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module HandData
using DelimitedFiles

export HandModel, HandInput, HandOutput, empty_hand_output, load_hand_input

struct HandModel
    bone_names::Vector{String}
    parents::Vector{Int}
    base_relatives::Vector{Matrix{Float64}}
    inverse_base_absolutes::Vector{Matrix{Float64}}
    base_positions::Matrix{Float64}
    weights::Matrix{Float64}
    triangles::Vector{Vector{Int}}
    is_mirrored::Bool
end

struct HandInput
    model::HandModel
    correspondences::Vector{Int}
    points::Matrix{Float64}
    theta::Vector{Float64}
    "Is present only for 'complicated' kind of problems."
    us::Union{Vector{Vector{Float64}}, Nothing}
end

mutable struct HandOutput
    objective::Vector{Float64}
    jacobian::Matrix{Float64}
end

empty_hand_output() = HandOutput([], Matrix{Float64}(undef, 0, 0))

function load_hand_model(model_dir::AbstractString)::HandModel
    delimeter = ':'
    bones_fn = joinpath(model_dir, "bones.txt")
    A = readdlm(bones_fn, delimeter)

    n_bones = size(A, 1)
    bone_names = A[:, 1]
    parents = A[:, 2] .+ 1 #julia indexing

    transforms = A[:, 3:18]
    transforms = permutedims(reshape(transforms, (n_bones, 4, 4)), (3, 2, 1))
    base_relatives = [ transforms[:, :, i] for i ∈ 1:size(transforms, 3) ]

    inverse_absolute_transforms  = A[:, 19:34]
    inverse_absolute_transforms = permutedims(reshape(inverse_absolute_transforms, (n_bones, 4, 4)), (3, 2, 1));
    inverse_base_absolutes = [ inverse_absolute_transforms[:, :, i] for i ∈ 1:size(inverse_absolute_transforms, 3) ]

    vertices_fn = joinpath(model_dir, "vertices.txt")
    A = readdlm(vertices_fn, delimeter);

    base_positions = A[:, 1:3]';
    n_vertices = size(A, 1);

    weights = zeros(Float64, n_bones, n_vertices);
    for i_vert = 1:n_vertices
        for i=0:A[i_vert, 9]-1
            i_bone = A[i_vert, 9 + i * 2 + 1] + 1; #julia indexing
            weights[i_bone, i_vert] = A[i_vert, 9 + i * 2 + 2];
        end
    end

    triangles_fn = joinpath(model_dir, "triangles.txt")
    T = readdlm(triangles_fn, delimeter, Int)
    triangles = [ T[i, :] .+ 1 for i ∈ 1:size(T, 1) ] #julia indexing

    HandModel(
        bone_names,
        parents,
        base_relatives,
        inverse_base_absolutes,
        [base_positions; ones(Float64, 1, size(base_positions, 2))],
        weights,
        triangles,
        false
    )
end

function load_hand_input(fn::AbstractString, iscomplicated::Bool)::HandInput
    modeldir = joinpath(dirname(fn), "model")
    model = load_hand_model(modeldir)
    open(fn) do io
        line = split(readline(io), " ")
        n_pts = parse(Int, line[1])
        n_theta = parse(Int, line[2])

        corrs = Vector{Int}(undef, n_pts)
        pts = Matrix{Float64}(undef, 3, n_pts)
        for i in 1:n_pts
            line = split(readline(io), " ")
            corrs[i] = parse(Int, line[1]) + 1 #julia indexing
            for j in 1:3
                pts[j, i] = parse(Float64, line[j + 1])
            end
        end

        if iscomplicated
            us = [parse.(Float64, split(readline(io), " ")) for i=1:n_pts]
        else
            us = nothing
        end

        theta = [parse(Float64, readline(io)) for i=1:n_theta]

        return HandInput(model, corrs, pts, theta, us)
    end
end

end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
"""
This module provides `EquationOfState` types and `apply` methods to calculate
energy, pressure, or bulk modulus of an `EquationOfState` on
a (an) volume (array of volumes).
"""
module Collections

using Unitful
using Unitful: AbstractQuantity

using EquationsOfState: EnergyForm, PressureForm, BulkModulusForm

export apply,
       EquationOfState,
       FiniteStrainEquationOfState,
       Murnaghan,
       BirchMurnaghan2nd,
       BirchMurnaghan3rd,
       BirchMurnaghan4th,
       PoirierTarantola2nd,
       PoirierTarantola3rd,
       PoirierTarantola4th,
       Vinet,
       AntonSchmidt,
       BreenanStacey

# ============================================================================ #
#                                     Types                                    #
# ============================================================================ #
"""
    EquationOfState{T}

An abstraction of equations of state, where `T` specifies the elements' common type.
"""
abstract type EquationOfState{T} end

"""
    FiniteStrainEquationOfState{T} <: EquationOfState{T}

An abstraction of finite strain equations of state, where `T` specifies the elements' common type.
"""
abstract type FiniteStrainEquationOfState{T} <: EquationOfState{T} end

"""
    Murnaghan(v0, b0, bp0, e0)

Create a Murnaghan equation of state. The elements' type will be handled automatically.

This equation of state can have units. The units are specified in [`Unitful.jl`](https://github.com/PainterQubits/Unitful.jl)'s
`@u_str` style.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `bp0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.

# Examples
```jldoctest
julia> Murnaghan(1, 2, 3.0)
Murnaghan{Float64}(1.0, 2.0, 3.0, 0.0)

julia> Murnaghan(Int8(1), 2//1, 3.0, 4)
Murnaghan{Float64}(1.0, 2.0, 3.0, 4.0)

julia> Murnaghan(1u"nm^3", 2u"GPa", 3, 3.0u"eV")
Murnaghan{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 3.0, 3.0 eV)
```
"""
struct Murnaghan{T} <: EquationOfState{T}
    v0::T
    b0::T
    bp0::T
    e0::T
end
function Murnaghan(v0, b0, bp0, e0)
    T = Base.promote_typeof(v0, b0, bp0, e0)
    return Murnaghan{T}(convert.(T, [v0, b0, bp0, e0])...)
end
Murnaghan(v0::Real, b0::Real, bp0::Real) = Murnaghan(v0, b0, bp0, 0)
Murnaghan(v0::AbstractQuantity, b0::AbstractQuantity, bp0) =
    Murnaghan(v0, b0, bp0, 0 * upreferred(Unitful.J))

"""
    BirchMurnaghan2nd(v0, b0, e0)

Create a Birch–Murnaghan 2nd order equation of state. The elements' type will be handled automatically.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.

See also: [`BirchMurnaghan3rd`](@ref), [`BirchMurnaghan4th`](@ref)

# Examples
```jldoctest
julia> BirchMurnaghan2nd(1, 2.0)
BirchMurnaghan2nd{Float64}(1.0, 2.0, 0.0)

julia> BirchMurnaghan2nd(Int8(1), 2//1, 0.0)
BirchMurnaghan2nd{Float64}(1.0, 2.0, 0.0)

julia> BirchMurnaghan2nd(1u"nm^3", 2u"GPa", 3.0u"eV")
BirchMurnaghan2nd{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 3.0 eV)
```
"""
struct BirchMurnaghan2nd{T} <: FiniteStrainEquationOfState{T}
    v0::T
    b0::T
    e0::T
end
function BirchMurnaghan2nd(v0, b0, e0)
    T = Base.promote_typeof(v0, b0, e0)
    return BirchMurnaghan2nd{T}(convert.(T, [v0, b0, e0])...)
end
BirchMurnaghan2nd(v0::Real, b0::Real) = BirchMurnaghan2nd(v0, b0, 0)
BirchMurnaghan2nd(v0::AbstractQuantity, b0::AbstractQuantity) =
    BirchMurnaghan2nd(v0, b0, 0 * upreferred(Unitful.J))

"""
    BirchMurnaghan3rd(v0, b0, bp0, e0)

Create a Birch–Murnaghan 3rd order equation of state. The elements' type will be handled automatically.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `bp0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.

See also: [`BirchMurnaghan2nd`](@ref), [`BirchMurnaghan4th`](@ref)

# Examples
```jldoctest
julia> BirchMurnaghan3rd(1, 2.0, 3)
BirchMurnaghan3rd{Float64}(1.0, 2.0, 3.0, 0.0)

julia> BirchMurnaghan3rd(Int8(1), 2//1, 4, 0.0)
BirchMurnaghan3rd{Float64}(1.0, 2.0, 4.0, 0.0)

julia> BirchMurnaghan3rd(1u"nm^3", 2u"GPa", 4.0, 3u"eV")
BirchMurnaghan3rd{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 4.0, 3.0 eV)
```
"""
struct BirchMurnaghan3rd{T} <: FiniteStrainEquationOfState{T}
    v0::T
    b0::T
    bp0::T
    e0::T
end
function BirchMurnaghan3rd(v0, b0, bp0, e0)
    T = Base.promote_typeof(v0, b0, bp0, e0)
    return BirchMurnaghan3rd{T}(convert.(T, [v0, b0, bp0, e0])...)
end
BirchMurnaghan3rd(v0::Real, b0::Real, bp0::Real) = BirchMurnaghan3rd(v0, b0, bp0, 0)
BirchMurnaghan3rd(v0::AbstractQuantity, b0::AbstractQuantity, bp0) =
    BirchMurnaghan3rd(v0, b0, bp0, 0 * upreferred(Unitful.J))

"""
    BirchMurnaghan4th(v0, b0, bp0, bpp0, e0)

Create a Birch–Murnaghan 4th order equation of state. The elements' type will be handled automatically.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `bp0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
- `bpp0`: the second-order pressure-derivative bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.

See also: [`BirchMurnaghan2nd`](@ref), [`BirchMurnaghan4th`](@ref)

# Examples
```jldoctest
julia> BirchMurnaghan4th(1, 2.0, 3, 4)
BirchMurnaghan4th{Float64}(1.0, 2.0, 3.0, 4.0, 0.0)

julia> BirchMurnaghan4th(Int8(1), 2//1, 4, 5.0, Float16(6))
BirchMurnaghan4th{Float64}(1.0, 2.0, 4.0, 5.0, 6.0)

julia> BirchMurnaghan4th(1u"nm^3", 2u"GPa", 3.0, 4u"1/GPa", 5u"eV")
BirchMurnaghan4th{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 3.0, 4.0 GPa^-1, 5.0 eV)
```
"""
struct BirchMurnaghan4th{T} <: FiniteStrainEquationOfState{T}
    v0::T
    b0::T
    bp0::T
    bpp0::T
    e0::T
end
function BirchMurnaghan4th(v0, b0, bp0, bpp0, e0)
    T = Base.promote_typeof(v0, b0, bp0, bpp0, e0)
    return BirchMurnaghan4th{T}(convert.(T, [v0, b0, bp0, bpp0, e0])...)
end
BirchMurnaghan4th(v0::Real, b0::Real, bp0::Real, bpp0::Real) =
    BirchMurnaghan4th(v0, b0, bp0, bpp0, 0)
BirchMurnaghan4th(
    v0::AbstractQuantity,
    b0::AbstractQuantity,
    bp0,
    bpp0::AbstractQuantity,
) = BirchMurnaghan4th(v0, b0, bp0, bpp0, 0 * upreferred(Unitful.J))

"""
    PoirierTarantola2nd(v0, b0, e0)

Create a Poirier–Tarantola order equation of state. The elements' type will be handled automatically.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.

See also: [`PoirierTarantola3rd`](@ref), [`PoirierTarantola4th`](@ref)

# Examples
```jldoctest
julia> PoirierTarantola2nd(1, 2.0)
PoirierTarantola2nd{Float64}(1.0, 2.0, 0.0)

julia> PoirierTarantola2nd(Int8(1), 2//1, 3.0)
PoirierTarantola2nd{Float64}(1.0, 2.0, 3.0)

julia> PoirierTarantola2nd(1u"nm^3", 2u"GPa", 3.0u"eV")
PoirierTarantola2nd{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 3.0 eV)
```
"""
struct PoirierTarantola2nd{T} <: FiniteStrainEquationOfState{T}
    v0::T
    b0::T
    e0::T
end
function PoirierTarantola2nd(v0, b0, e0)
    T = Base.promote_typeof(v0, b0, e0)
    return PoirierTarantola2nd{T}(convert.(T, [v0, b0, e0])...)
end
PoirierTarantola2nd(v0::Real, b0::Real) = PoirierTarantola2nd(v0, b0, 0)
PoirierTarantola2nd(v0::AbstractQuantity, b0::AbstractQuantity) =
    PoirierTarantola2nd(v0, b0, 0 * upreferred(Unitful.J))

"""
    PoirierTarantola3rd(v0, b0, bp0, e0)

Create a Poirier–Tarantola 3rd order equation of state. The elements' type will be handled automatically.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `bp0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.

See also: [`PoirierTarantola2nd`](@ref), [`PoirierTarantola4th`](@ref)

# Examples
```jldoctest
julia> PoirierTarantola3rd(1, 2.0, 3)
PoirierTarantola3rd{Float64}(1.0, 2.0, 3.0, 0.0)

julia> PoirierTarantola3rd(Int8(1), 2//1, 3.0, Float16(4))
PoirierTarantola3rd{Float64}(1.0, 2.0, 3.0, 4.0)

julia> PoirierTarantola3rd(1u"nm^3", 2u"GPa", 3, 4.0u"eV")
PoirierTarantola3rd{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 3.0, 4.0 eV)
```
"""
struct PoirierTarantola3rd{T} <: FiniteStrainEquationOfState{T}
    v0::T
    b0::T
    bp0::T
    e0::T
end
function PoirierTarantola3rd(v0, b0, bp0, e0)
    T = Base.promote_typeof(v0, b0, bp0, e0)
    return PoirierTarantola3rd{T}(convert.(T, [v0, b0, bp0, e0])...)
end
PoirierTarantola3rd(v0::Real, b0::Real, bp0::Real) = PoirierTarantola3rd(v0, b0, bp0, 0)
PoirierTarantola3rd(v0::AbstractQuantity, b0::AbstractQuantity, bp0) =
    PoirierTarantola3rd(v0, b0, bp0, 0 * upreferred(Unitful.J))

"""
    PoirierTarantola4th(v0, b0, bp0, bpp0, e0)

Create a Birch–Murnaghan 4th order equation of state. The elements' type will be handled automatically.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `bp0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
- `bpp0`: the second-order pressure-derivative bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.

See also: [`PoirierTarantola2nd`](@ref), [`PoirierTarantola3rd`](@ref)

# Examples
```jldoctest
julia> PoirierTarantola4th(1, 2.0, 3, 4)
PoirierTarantola4th{Float64}(1.0, 2.0, 3.0, 4.0, 0.0)

julia> PoirierTarantola4th(Int8(1), 2//1, 3.0, Float16(4), 5)
PoirierTarantola4th{Float64}(1.0, 2.0, 3.0, 4.0, 5.0)

julia> PoirierTarantola4th(1u"nm^3", 2u"GPa", 3, 4u"1/GPa", 5.0u"eV")
PoirierTarantola4th{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 3.0, 4.0 GPa^-1, 5.0 eV)
```
"""
struct PoirierTarantola4th{T} <: FiniteStrainEquationOfState{T}
    v0::T
    b0::T
    bp0::T
    bpp0::T
    e0::T
end
function PoirierTarantola4th(v0, b0, bp0, bpp0, e0)
    T = Base.promote_typeof(v0, b0, bp0, bpp0, e0)
    return PoirierTarantola4th{T}(convert.(T, [v0, b0, bp0, bpp0, e0])...)
end
PoirierTarantola4th(v0::Real, b0::Real, bp0::Real, bpp0::Real) =
    PoirierTarantola4th(v0, b0, bp0, bpp0, 0)
PoirierTarantola4th(
    v0::AbstractQuantity,
    b0::AbstractQuantity,
    bp0,
    bpp0::AbstractQuantity,
) = PoirierTarantola4th(v0, b0, bp0, bpp0, 0 * upreferred(Unitful.J))

"""
    Vinet(v0, b0, bp0, e0)

Create a Vinet equation of state. The elements' type will be handled automatically.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `bp0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure. Its default value is `0u"eV"` (`0`), if other parameters have (no) units.

# Examples
```jldoctest
julia> Vinet(1, 2.0, 3)
Vinet{Float64}(1.0, 2.0, 3.0, 0.0)

julia> Vinet(Int8(1), 2//1, 3.0, Float16(4))
Vinet{Float64}(1.0, 2.0, 3.0, 4.0)

julia> Vinet(1u"nm^3", 2u"GPa", 3, 4.0u"eV")
Vinet{Quantity{Float64,D,U} where U where D}(1.0 nm^3, 2.0 GPa, 3.0, 4.0 eV)
```
"""
struct Vinet{T} <: EquationOfState{T}
    v0::T
    b0::T
    bp0::T
    e0::T
end
function Vinet(v0, b0, bp0, e0)
    T = Base.promote_typeof(v0, b0, bp0, e0)
    return Vinet{T}(convert.(T, [v0, b0, bp0, e0])...)
end
Vinet(v0::Real, b0::Real, bp0::Real) = Vinet(v0, b0, bp0, 0)
Vinet(v0::AbstractQuantity, b0::AbstractQuantity, bp0) =
    Vinet(v0, b0, bp0, 0 * upreferred(Unitful.J))

struct AntonSchmidt{T} <: EquationOfState{T}
    v0::T
    β::T
    n::T
    e∞::T
end
function AntonSchmidt(v0, β, n, e∞)
    T = Base.promote_typeof(v0, β, n, e∞)
    return AntonSchmidt{T}(convert.(T, [v0, β, n, e∞])...)
end
AntonSchmidt(v0::Real, β::Real, n::Real) = AntonSchmidt(v0, β, n, 0)

struct BreenanStacey{T} <: EquationOfState{T}
    v0::T
    b0::T
    γ0::T
    e0::T
end
function BreenanStacey(v0, b0, γ0, e0)
    T = Base.promote_typeof(v0, b0, γ0, e0)
    return BreenanStacey{T}(convert.(T, [v0, b0, γ0, e0])...)
end
BreenanStacey(v0::Real, b0::Real, γ0::Real) = BreenanStacey(v0, b0, γ0, 0)
# =================================== Types ================================== #

# ============================================================================ #
#                               Energy evaluation                              #
# ============================================================================ #
"""
    apply(EnergyForm(), eos::EquationOfState)
    apply(PressureForm(), eos::EquationOfState)
    apply(BulkModulusForm(), eos::EquationOfState)

Return a function that takes a volume as a variable, suitable for mapping onto an array.

# Examples
```jldoctest
julia> using EquationsOfState, EquationsOfState.Collections

julia> f = apply(EnergyForm(), Vinet(1, 2, 3));

julia> map(f, 1:1:10)
10-element Array{Float64,1}:
 0.0
 0.367905230584308
 0.7652477289745814
 1.0516459435179233
 1.2560420090256408
 1.405149833626178
 1.5165867441792136
 1.6017034530570884
 1.6679539823686644
 1.7203642945516917

julia> g = apply(PressureForm(), Vinet(1, 2, 3));

julia> map(g, 1:1:10)
10-element Array{Float64,1}:
  0.0
 -0.45046308428750254
 -0.3384840350043251
 -0.24010297221667418
 -0.17314062272722755
 -0.12795492664586872
 -0.09677154467733216
 -0.07468060255179591
 -0.05864401631176751
 -0.04674768462396211

julia> h = apply(BulkModulusForm(), BirchMurnaghan3rd(1, 2, 3));

julia> map(h, 1:1:10)
10-element Array{Float64,1}:
 2.0
 0.9216086833346415
 0.444903691617472
 0.2540009203153288
 0.16193296566524193
 0.11130584492987289
 0.08076305569984538
 0.06103515625
 0.047609811583958425
 0.03808959181078831
```
"""
apply(form::EnergyForm, eos::EquationOfState) = v -> apply(form, eos, v)
"""
    apply(EnergyForm(), eos::EquationOfState, v)

Return the energy of an `EquationOfState` on volume `v`. If `eos` has units,
`v` must also has.
"""
function apply(::EnergyForm, eos::Murnaghan, v)
    v0, b0, bp0, e0 = fieldvalues(eos)

    x = bp0 - 1
    y = (v0 / v)^bp0
    return e0 + b0 / bp0 * v * (y / x + 1) - v0 * b0 / x
end
function apply(::EnergyForm, eos::BirchMurnaghan2nd, v)
    v0, b0, e0 = fieldvalues(eos)

    f = (cbrt(v0 / v)^2 - 1) / 2
    return e0 + 9 / 2 * b0 * v0 * f^2
end
function apply(::EnergyForm, eos::BirchMurnaghan3rd, v)
    v0, b0, bp0, e0 = fieldvalues(eos)

    eta = cbrt(v0 / v)
    xi = eta^2 - 1
    return e0 + 9 / 16 * b0 * v0 * xi^2 * (6 + bp0 * xi - 4 * eta^2)
end
function apply(::EnergyForm, eos::BirchMurnaghan4th, v)
    v0, b0, bp0, bpp0, e0 = fieldvalues(eos)

    f = (cbrt(v0 / v)^2 - 1) / 2
    h = b0 * bpp0 + bp0^2
    return e0 + 3 / 8 * v0 * b0 * f^2 * ((9h - 63bp0 + 143) * f^2 + 12 * (bp0 - 4) * f + 12)
end
function apply(::EnergyForm, eos::PoirierTarantola2nd, v)
    v0, b0, e0 = fieldvalues(eos)

    return e0 + b0 / 2 * v0 * log(v / v0)^(2 / 3)
end
function apply(::EnergyForm, eos::PoirierTarantola3rd, v)
    v0, b0, bp0, e0 = fieldvalues(eos)

    x = cbrt(v / v0)
    xi = -3 * log(x)
    return e0 + b0 / 6 * v0 * xi^2 * ((bp0 - 2) * xi + 3)
end
function apply(::EnergyForm, eos::PoirierTarantola4th, v)
    v0, b0, bp0, bpp0, e0 = fieldvalues(eos)

    x = cbrt(v / v0)
    xi = log(x)
    h = b0 * bpp0 + bp0^2
    return e0 + b0 / 24v0 * xi^2 * ((h + 3bp0 + 3) * xi^2 + 4 * (bp0 + 2) * xi + 12)
end
function apply(::EnergyForm, eos::Vinet, v)
    v0, b0, bp0, e0 = fieldvalues(eos)

    x = cbrt(v / v0)
    xi = 3 / 2 * (bp0 - 1)
    return e0 + 9b0 * v0 / xi^2 * (1 + (xi * (1 - x) - 1) * exp(xi * (1 - x)))
end
function apply(::EnergyForm, eos::AntonSchmidt, v)
    v0, β, n, e∞ = fieldvalues(eos)

    x = v / v0
    η = n + 1
    return e∞ + β * v0 / η * x^η * (log(x) - 1 / η)
end
# ============================= Energy evaluation ============================ #

# ============================================================================ #
#                              Pressure evaluation                             #
# ============================================================================ #
apply(::PressureForm, eos::EquationOfState) = v -> apply(PressureForm(), eos, v)
"""
    apply(PressureForm(), eos::EquationOfState, v)

Return the pressure of an `EquationOfState` on volume `v`. If `eos` has units,
`v` must also has.
"""
function apply(::PressureForm, eos::Murnaghan, v)
    v0, b0, bp0 = fieldvalues(eos)

    return b0 / bp0 * ((v0 / v)^bp0 - 1)
end
function apply(::PressureForm, eos::BirchMurnaghan2nd, v)
    v0, b0 = fieldvalues(eos)

    f = ((v0 / v)^(2 / 3) - 1) / 2
    return 3b0 * f * (1 + 2f)^(5 / 2)
end
function apply(::PressureForm, eos::BirchMurnaghan3rd, v)
    v0, b0, bp0 = fieldvalues(eos)

    eta = (v0 / v)^(1 / 3)
    return 3 / 2 * b0 * (eta^7 - eta^5) * (1 + 3 / 4 * (bp0 - 4) * (eta^2 - 1))
end
function apply(::PressureForm, eos::BirchMurnaghan4th, v)
    v0, b0, bp0, bpp0 = fieldvalues(eos)

    f = ((v0 / v)^(2 / 3) - 1) / 2
    h = b0 * bpp0 + bp0^2
    return b0 / 2 * (2f + 1)^(5 / 2) * ((9h - 63bp0 + 143) * f^2 + 9 * (bp0 - 4) * f + 6)
end
function apply(::PressureForm, eos::PoirierTarantola2nd, v)
    v0, b0 = fieldvalues(eos)

    x = (v / v0)^(1 / 3)
    return -b0 / x * log(x)
end
function apply(::PressureForm, eos::PoirierTarantola3rd, v)
    v0, b0, bp0 = fieldvalues(eos)

    x = v / v0
    xi = log(x)
    return -b0 * xi / 2x * ((bp0 - 2) * xi - 2)
end
function apply(::PressureForm, eos::PoirierTarantola4th, v)
    v0, b0, bp0, bpp0 = fieldvalues(eos)

    x = (v / v0)^(1 / 3)
    xi = log(x)
    h = b0 * bpp0 + bp0^2
    return -b0 * xi / 6 / x * ((h + 3bp0 + 3) * xi^2 + 3 * (bp0 + 6) * xi + 6)
end
function apply(::PressureForm, eos::Vinet, v)
    v0, b0, bp0 = fieldvalues(eos)

    x = (v / v0)^(1 / 3)
    xi = 3 / 2 * (bp0 - 1)
    return 3b0 / x^2 * (1 - x) * exp(xi * (1 - x))
end
function apply(::PressureForm, eos::AntonSchmidt, v)
    v0, β, n = fieldvalues(eos)

    x = v / v0
    return -β * x^n * log(x)
end
function apply(::PressureForm, eos::BreenanStacey, v)
    v0, b0, γ0 = fieldvalues(eos)

    x = v0 / v
    return b0 / 2 / γ0 * x^(4 / 3) * (exp(2γ0 * (1 - x)) - 1)
end
# ============================ Pressure evaluation =========================== #

# ============================================================================ #
#                            Bulk modulus evaluation                           #
# ============================================================================ #
apply(::BulkModulusForm, eos::EquationOfState) = v -> apply(BulkModulusForm(), eos, v)
"""
    apply(BulkModulusForm(), eos::EquationOfState, v)

Return the bulk modulus of an `EquationOfState` on volume `v`. If `eos` has units,
`v` must also has.
"""
function apply(::BulkModulusForm, eos::BirchMurnaghan2nd, v)
    v0, b0 = fieldvalues(eos)

    f = ((v0 / v)^(2 / 3) - 1) / 2
    return b0 * (7f + 1) * (2f + 1)^(5 / 2)
end
function apply(::BulkModulusForm, eos::BirchMurnaghan3rd, v)
    v0, b0, bp0 = fieldvalues(eos)

    f = ((v0 / v)^(2 / 3) - 1) / 2
    return b0 / 2 * (2f + 1)^(5 / 2) * ((27 * f^2 + 6f) * (bp0 - 4) - 4f + 2)
end
function apply(::BulkModulusForm, eos::BirchMurnaghan4th, v)
    v0, b0, bp0, bpp0 = fieldvalues(eos)

    f = ((v0 / v)^(2 / 3) - 1) / 2
    h = b0 * bpp0 + bp0^2
    return b0 / 6 * (2f + 1)^(5 / 2) *
           ((99h - 693bp0 + 1573) * f^3 + (27h - 108bp0 + 105) * f^2 + 6f * (3bp0 - 5) + 6)
end
function apply(::BulkModulusForm, eos::PoirierTarantola2nd, v)
    v0, b0 = fieldvalues(eos)

    x = (v / v0)^(1 / 3)
    return b0 / x * (1 - log(x))
end
function apply(::BulkModulusForm, eos::PoirierTarantola3rd, v)
    v0, b0, bp0 = fieldvalues(eos)

    x = v / v0
    xi = log(x)
    return -b0 / 2x * (((bp0 - 2) * xi + 2 - 2bp0) * xi + 2)
end
function apply(::BulkModulusForm, eos::PoirierTarantola4th, v)
    v0, b0, bp0, bpp0 = fieldvalues(eos)

    x = (v / v0)^(1 / 3)
    xi = log(x)
    h = b0 * bpp0 + bp0^2
    return -b0 / (6x) *
           ((h + 3bp0 + 3) * xi^3 - 3 * xi^2 * (h + 2bp0 + 1) - 6xi * (bp0 + 1) - 6)
end
function apply(::BulkModulusForm, eos::Vinet, v)
    v0, b0, bp0 = fieldvalues(eos)

    x = (v / v0)^(1 / 3)
    xi = 3 / 2 * (bp0 - 1)
    return -b0 / (2 * x^2) * (3x * (x - 1) * (bp0 - 1) + 2 * (x - 2)) * exp(-xi * (x - 1))
end
function apply(::BulkModulusForm, eos::AntonSchmidt, v)
    v0, β, n = fieldvalues(eos)

    x = v / v0
    return β * x^n * (1 + n * log(x))
end
# ========================== Bulk modulus evaluation ========================= #

# ============================================================================ #
#                                 Miscellaneous                                #
# ============================================================================ #
# This is a helper function and should not be exported.
fieldvalues(eos::EquationOfState) = [getfield(eos, i) for i in 1:nfields(eos)]

function Base.:(==)(x::T, y::T) where {T<:EquationOfState}
    return all(getfield(x, i) == getfield(y, i) for i in 1:fieldcount(T))
end

Unitful.upreferred(::typeof(dimension(u"J"))) = u"eV"
Unitful.upreferred(::typeof(dimension(u"m^3"))) = u"angstrom^3"
Unitful.upreferred(::typeof(dimension(u"Pa"))) = u"eV/angstrom^3"
# =============================== Miscellaneous ============================== #

end


text:
# Copyright 2021 Gravifer
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Pkg
Pkg.activate(normpath(joinpath(@__DIR__, "../..")))
# using Revise
using BenchmarkTools
using Plots
using OffsetArrays

myblue = RGB(0.368417, 0.506779, 0.709798) # <==> 94, 129, 181 (the default style in Wolfram Mathematica)

"""
    yt0012(ξ)

Compute ``yₜ`` (i.e., the half-thickness relative to the chord)
of NACA0012 airfoils from the normalized chord position ``ξ=x/c``.

Additionally defined upon ``[1, 1.005]`` for a 2-order smooth closure.

See also: https://en.wikipedia.org/wiki/NACA_airfoil

TODO: generalize for 00xx airfoils; furthermore, for all 4-digit airfoils.

# Examples
```jldoctest
julia> yt0012(0.13)
0.051193437032115506
```
"""
Base.@pure function yt0012(ξ::T where {T<:Real})
    ε = 0.0050544107511712
    return 5 * 0.12 *  (
        if    0 <= ξ <= 1
        0.2969 * sqrt(ξ)  -
        0.1260 *      ξ   -
        0.3516 *      ξ^2 +
        0.2843 *      ξ^3 -
        0.1015 *      ξ^4
        elseif 1 < ξ <= 1+ε # 2-order smooth joining of the tail edge.
        0.03423616658680458 * sqrt(1+ε - ξ)   +
        0.12491972188290168 *     (1+ε - ξ)   +
        11.579398976610927   *     (1+ε - ξ)^2 +
        12.21202728162769    *     (1+ε - ξ)^3
        else
            0
        end)
end

const Mx, My = 100+1, 20+1

include(normpath(joinpath(@__DIR__, "../../src/misc-util.jl")))
using .__CFD2021__misc_util__: tuplejoin
"""
- usage of ``tuplejoin``:
```jldoctest
  julia> temp_v = ([0,0.5,1], [0,0.5,1])
  julia> temp_t = (axes(v) for v in temp_v)
  julia> @btime tuplejoin(temp_t...)
```
"""

include(normpath(joinpath(@__DIR__, "../src/transfinite-interpolate.jl")))
import .transfinite_interpolate: transfinite_interpolate_2d!, transfinite_interpolate_2d

"""
- Try to use the in-place version when applicable;
  there is a tremendous performance boost.
```jldoctest
  julia> v_lo = ([LinRange(0, 1, 10^4)...], [LinRange(0, 1, 10^4)...]); v_hi = ([LinRange(1, 2, 10^4)...], [LinRange(1, 2, 10^4)...])

  julia> interpolated_array = @btime Array{Float64}(undef, (length(v) for v in v_lo)...)
  20.500 μs (7 allocations: 762.94 MiB)

  julia> @btime transfinite_interpolate_2d!(interpolated_array, v_lo, v_hi)
  1.535 s (1 allocation: 96 bytes)

  julia> res = @btime transfinite_interpolate_2d(v_lo, v_hi)
  7.581 s (100000005 allocations: 2.98 GiB)
```
"""

"""
- example of a cartesian grid;
```jldoctest
  x_lo = ([LinRange(0, 1, Mx)...], [LinRange(0, 0, My)...]); x_hi = ([LinRange(0, 1, Mx)...], [LinRange(1, 1, My)...]);
  y_lo = ([LinRange(0, 0, Mx)...], [LinRange(0, 1, My)...]); y_hi = ([LinRange(1, 1, Mx)...], [LinRange(0, 1, My)...]);
```
"""

"""
- example of a half-O grid aroundNACA0012;
"""
x_lo = ([sinpi(ξ)^2 for ξ in LinRange(0, 0.5, Mx)], -exp.(LinRange(0, log(6), My)) .+ 1);
y_lo = (yt0012.(x_lo[1]), [LinRange(0, 0, My)...]);
x_hi = (5cospi.(LinRange(1, 0, Mx)), +exp.(LinRange(0, log(5), My)));
y_hi = (5sinpi.(LinRange(1, 0, Mx)), [LinRange(0, 0, My)...]);

"""
- example of a whole-O grid aroundNACA0012;
```jldoctest
  x_lo = ([         sinpi.(LinRange(-.5,0, Mx)).^2;          sinpi.(LinRange(0, .5, Mx)).^2],
             exp.(LinRange(0, log(5), My)                  )   );
  y_lo = ([-yt0012.(sinpi.(LinRange(-.5,0, Mx)).^2); yt0012.(sinpi.(LinRange(0, .5, Mx)).^2)],
          [LinRange(0, 0, My)...]);
  x_hi = (5cospi.(LinRange( 2, 0, 2Mx)),
             exp.(LinRange(0, log(5), My)                  )   );
  y_hi = (5sinpi.(LinRange( 2, 0, 2Mx)),
          [LinRange(0, 0, My)...]);
```
"""

p = plot(; xlims=(-5, 5), ylims=(-5, 5), aspect_ratio=:equal)
# p = plot(xlims = (-.6, .6), ylims = (0, .5), aspect_ratio = :equal)
interpolated_x = Array{Float64}(undef, (length(v) for v in x_lo)...);
interpolated_y = Array{Float64}(undef, (length(v) for v in y_lo)...);
@btime transfinite_interpolate_2d!(interpolated_x, x_lo, x_hi)
@btime transfinite_interpolate_2d!(interpolated_y, y_lo, y_hi)
for u in axes(interpolated_x)[begin]  plot!(p,interpolated_x[u,:], interpolated_y[u,:], label = nothing, color = myblue); end
for u in axes(interpolated_x)[ end ]  plot!(p,interpolated_x[:,u], interpolated_y[:,u], label = nothing, color = myblue); end
# for u in axes(interpolated_x)[begin]  plot!(p,interpolated_x[u,:],-interpolated_y[u,:], label = nothing, color = myblue); end
# for u in axes(interpolated_x)[ end ]  plot!(p,interpolated_x[:,u],-interpolated_y[:,u], label = nothing, color = myblue); end
plot!(p, yt0012, 0, 1; label=nothing, color=:red)
display(p) # savefig(normpath(joinpath(@__DIR__, "img/tf-hO.svg")))


text:
module OptEvapCool

export
kB, g, c, ε₀, h̄, # Constants
AtomSpecies, Rb87, kappa, # Atoms
uniform_cloud, boltzmann_velocities, # Atom initialisation
harmonic_boltzmann_positions,
GaussianBeam, HarmonicField, UniformField, # Field initialisation
gravity, acceleration, potential,
no_evap, energy_evap, ellipsoid_evap, # Evaporation
exponential_ramp, linear_ramp,
GlobalSensor, measurer, measure, # Measurement and plotting
plot_energy, plot_temperature, plot_speed,
plot_collrate, plot_number, plot_density,
plot_psd, rollmean, rolling,
equilibrium_collrate, equilibrium_speeds, # Theoretical predictions
harmonic_theory, purdue_theory,
savefig, savecsv, filetime, loadsensor, # File saving/loading
SimulationConditions, # Setup
evolve # Perform the simulation

include("./conditions.jl")
include("./measure.jl")
include("./dsmc.jl")
include("./theory.jl")

end


text:
fname = :ARCLDLt
c = Combi(hessian_dense,PDataLDLt,solve_modelARCDiag,preprocessLDLt,decreaseFact,Tparam())
include("Template.jl")


text:
#=
    Sorting networks that process 1,2,..16 values are given.
    Sorting networks that process 17..25 values are given in "swapsort_17to25.jl".

    These sorting networks should minimize conditional swaps.
    The first sixteen are known to minimize conditional swaps.

    Functions are written with groups of internally parallelizable
    statements given together and sequential steps are separated.

    `swapsort(x1::T, .., xn::T)::NTuple{n,T}`
    `swapsort(tup::NTuple{n,T})::NTuple{n,T}`
    `swapsort(vec::Vector{T})::NTuple{N,T}`

    networks were selected using software by John Gamble
        http://pages.ripco.net/~jgamble/nw.html
        http://search.cpan.org/dist/Algorithm-Networksort/
=#

const VALS = ([Val{i} for i=1:16]...,)

@inline function swapsort(vec::Vector{T}) where {T}
    n = length(vec)
    n <= 16 && return swapsort(VALS[length(vec)], vec)
    throw(ErrorException("swapsort(vector) not defined where length(vector) > 16"))
end

#    sort 1 value with 0 minmaxs in 1 stage

function swapsort(a::T) where {T}
    return (a,)
end

@inline swapsort(x::NTuple{1,T}) where {T} = swapsort(x[1])
@inline swapsort(::Type{Val{1}}, x::Vector{T}) where {T} = swapsort(x[1])

#    sort 2 values with 1 minmaxs in 1 stage

function swapsort(a::T, b::T) where {T}

    a, b = minmax(a, b)

    return a,b
end

@inline swapsort(x::NTuple{2,T}) where {T} = swapsort(x[1], x[2])
@inline swapsort(::Type{Val{2}}, x::Vector{T}) where {T} = swapsort(x[1], x[2])

#    sort 3 values with 3 minmaxs in 3 parallel stages

function swapsort(a::T, b::T, c::T) where {T}

    b, c = minmax(b, c)

    a, c = minmax(a, c)

    a, b = minmax(a, b)

    return a, b, c
end

@inline swapsort(x::NTuple{3,T}) where {T} = swapsort(x[1], x[2], x[3])
@inline swapsort(::Type{Val{3}}, x::Vector{T}) where {T} = swapsort(x[1], x[2], x[3])

#    sort 4 values with 5 minmaxs in 3 parallel stages

function swapsort(a::T, b::T, c::T, d::T) where {T}

    a, b = minmax(a, b)
    c, d = minmax(c, d)

    a, c = minmax(a, c)
    b, d = minmax(b, d)

    b, c = minmax(b, c)

    return a, b, c, d
end

@inline swapsort(x::NTuple{4,T}) where {T} = swapsort(x[1], x[2], x[3], x[4])
@inline swapsort(::Type{Val{4}}, x::Vector{T}) where {T} = swapsort(x[1], x[2], x[3], x[4])

#    sort 5 values with 9 minmaxs in 5 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T) where {T}

    b, c = minmax(b, c)
    d, e = minmax(d, e)

    b, d = minmax(b, d)
    a, c = minmax(a, c)

    c, e = minmax(c, e)
    a, d = minmax(a, d)

    a, b = minmax(a, b)
    c, d = minmax(c, d)

    b, c = minmax(b, c)

    return a, b, c, d, e
end

@inline swapsort(x::NTuple{5,T}) where {T} = swapsort(x[1], x[2], x[3], x[4], x[5])
@inline swapsort(::Type{Val{5}}, x::Vector{T}) where {T} = swapsort(x[1], x[2], x[3], x[4], x[5])

#    sort 6 values with 12 minmaxs in 6 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T) where {T}

    b, c = minmax(b, c)
    e, f = minmax(e, f)

    a, c = minmax(a, c)
    d, f = minmax(d, f)

    a, b = minmax(a, b)
    d, e = minmax(d, e)
    c, f = minmax(c, f)

    a, d = minmax(a, d)
    b, e = minmax(b, e)

    c, e = minmax(c, e)
    b, d = minmax(b, d)

    c, d = minmax(c, d)

    return a, b, c, d, e, f
end

@inline swapsort(x::NTuple{6,T}) where {T} = swapsort(x[1], x[2], x[3], x[4], x[5], x[6])
@inline swapsort(::Type{Val{6}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6])

#    sort 7 values with 16 minmaxs in 6 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T) where {T}

    a, e = minmax(a, e)
    b, f = minmax(b, f)
    c, g = minmax(c, g)

    a, c = minmax(a, c)
    b, d = minmax(b, d)
    e, g = minmax(e, g)

    c, e = minmax(c, e)
    d, f = minmax(d, f)
    a, b = minmax(a, b)

    c, d = minmax(c, d)
    e, f = minmax(e, f)

    b, e = minmax(b, e)
    d, g = minmax(d, g)

    b, c = minmax(b, c)
    d, e = minmax(d, e)
    f, g = minmax(f, g)

    return a, b, c, d, e, f, g
end

@inline swapsort(x::NTuple{7,T}) where {T} = swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7])
@inline swapsort(::Type{Val{7}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7])

#    sort 8 values with 19 minmaxs in 7 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T) where {T}

    a, b = minmax(a, b)
    c, d = minmax(c, d)
    e, f = minmax(e, f)
    g, h = minmax(g, h)

    a, c = minmax(a, c)
    b, d = minmax(b, d)
    e, g = minmax(e, g)
    f, h = minmax(f, h)

    b, c = minmax(b, c)
    f, g = minmax(f, g)
    a, e = minmax(a, e)
    d, h = minmax(d, h)

    b, f = minmax(b, f)
    c, g = minmax(c, g)

    b, e = minmax(b, e)
    d, g = minmax(d, g)

    c, e = minmax(c, e)
    d, f = minmax(d, f)

    d, e = minmax(d, e)

    return a, b, c, d, e, f, g, h
end

@inline swapsort(x::NTuple{8,T}) where {T} = swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8])
@inline swapsort(::Type{Val{8}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8])

#    sort 9 values with 25 minmaxs in 9 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T) where {T}

    a, b = minmax(a, b)
    d, e = minmax(d, e)
    g, h = minmax(g, h)

    b, c = minmax(b, c)
    e, f = minmax(e, f)
    h, i = minmax(h, i)

    a, b = minmax(a, b)
    d, e = minmax(d, e)
    g, h = minmax(g, h)
    c, f = minmax(c, f)

    a, d = minmax(a, d)
    b, e = minmax(b, e)
    f, i = minmax(f, i)

    d, g = minmax(d, g)
    e, h = minmax(e, h)
    c, f = minmax(c, f)

    a, d = minmax(a, d)
    b, e = minmax(b, e)
    f, h = minmax(f, h)
    c, g = minmax(c, g)

    b, d = minmax(b, d)
    e, g = minmax(e, g)

    c, e = minmax(c, e)
    f, g = minmax(f, g)

    c, d = minmax(c, d)

    return a, b, c, d, e, f, g, h, i
end

@inline swapsort(x::NTuple{9,T}) where {T} = swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9])
@inline swapsort(::Type{Val{9}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9])

#    sort 10 values with 29 minmaxs in 9 parallel stages

function swapsort(a::T,b::T,c::T,d::T,e::T,
              f::T,g::T,h::T,i::T,j::T) where {T}

    e, j = minmax(e, j)
    d, i = minmax(d, i)
    c, h = minmax(c, h)
    b, g = minmax(b, g)
    a, f = minmax(a, f)

    b, e = minmax(b, e)
    g, j = minmax(g, j)
    a, d = minmax(a, d)
    f, i = minmax(f, i)

    a, c = minmax(a, c)
    d, g = minmax(d, g)
    h, j = minmax(h, j)

    a, b = minmax(a, b)
    c, e = minmax(c, e)
    f, h = minmax(f, h)
    i, j = minmax(i, j)

    b, c = minmax(b, c)
    e, g = minmax(e, g)
    h, i = minmax(h, i)
    d, f = minmax(d, f)

    c, f = minmax(c, f)
    g, i = minmax(g, i)
    b, d = minmax(b, d)
    e, h = minmax(e, h)

    c, d = minmax(c, d)
    g, h = minmax(g, h)

    d, e = minmax(d, e)
    f, g = minmax(f, g)

    e, f = minmax(e, f)

    return a,b,c,d,e,f,g,h,i,j
end

@inline swapsort(x::NTuple{10,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10])
@inline swapsort(::Type{Val{10}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10])

#    sort 11 values with 35 minmaxs in 9 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T) where {T}

    a, b = minmax(a, b)
    c, d = minmax(c, d)
    e, f = minmax(e, f)
    g, h = minmax(g, h)
    i, j = minmax(i, j)

    b, d = minmax(b, d)
    f, h = minmax(f, h)
    a, c = minmax(a, c)
    e, g = minmax(e, g)
    i, k = minmax(i, k)

    b, c = minmax(b, c)
    f, g = minmax(f, g)
    j, k = minmax(j, k)
    a, e = minmax(a, e)
    d, h = minmax(d, h)

    b, f = minmax(b, f)
    g, k = minmax(g, k)
    e, i = minmax(e, i)

    f, j = minmax(f, j)
    c, g = minmax(c, g)
    a, e = minmax(a, e)
    d, i = minmax(d, i)

    b, f = minmax(b, f)
    g, k = minmax(g, k)
    c, d = minmax(c, d)
    i, j = minmax(i, j)

    b, e = minmax(b, e)
    h, k = minmax(h, k)
    d, f = minmax(d, f)
    g, i = minmax(g, i)

    c, e = minmax(c, e)
    h, j = minmax(h, j)
    f, g = minmax(f, g)

    d, e = minmax(d, e)
    h, i = minmax(h, i)

    return a,b,c,d,e,f,g,h,i,j,k
end

@inline swapsort(x::NTuple{11,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11])
@inline swapsort(::Type{Val{11}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11])

#    sort 12 values with 39 minmaxs in 9 parallel stages

function swapsort(a::T,b::T,c::T,d::T,e::T,
              f::T,g::T,h::T,i::T,j::T,
              k::T,l::T) where {T}

    a, b = minmax(a, b)
    c, d = minmax(c, d)
    e, f = minmax(e, f)
    g, h = minmax(g, h)
    i, j = minmax(i, j)
    k, l = minmax(k, l)

    b, d = minmax(b, d)
    f, h = minmax(f, h)
    j, l = minmax(j, l)
    a, c = minmax(a, c)
    e, g = minmax(e, g)
    i, k = minmax(i, k)

    b, c = minmax(b, c)
    f, g = minmax(f, g)
    j, k = minmax(j, k)
    a, e = minmax(a, e)
    h, l = minmax(h, l)

    b, f = minmax(b, f)
    g, k = minmax(g, k)
    d, h = minmax(d, h)
    e, i = minmax(e, i)

    f, j = minmax(f, j)
    c, g = minmax(c, g)
    a, e = minmax(a, e)
    h, l = minmax(h, l)
    d, i = minmax(d, i)

    b, f = minmax(b, f)
    g, k = minmax(g, k)
    c, d = minmax(c, d)
    i, j = minmax(i, j)

    b, e = minmax(b, e)
    h, k = minmax(h, k)
    d, f = minmax(d, f)
    g, i = minmax(g, i)

    c, e = minmax(c, e)
    h, j = minmax(h, j)
    f, g = minmax(f, g)

    d, e = minmax(d, e)
    h, i = minmax(h, i)

    return a,b,c,d,e,f,g,h,i,j,k,l
end

@inline swapsort(x::NTuple{12,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12])
@inline swapsort(::Type{Val{12}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12])

#    sort 13 values with 45 minmaxs in 10 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T) where {T}

    b, h = minmax(b, h)
    j, l = minmax(j, l)
    d, e = minmax(d, e)
    f, i = minmax(f, i)
    a, m = minmax(a, m)
    c, g = minmax(c, g)

    a, b = minmax(a, b)
    c, d = minmax(c, d)
    e, g = minmax(e, g)
    i, l = minmax(i, l)
    h, m = minmax(h, m)
    f, j = minmax(f, j)

    a, c = minmax(a, c)
    d, h = minmax(d, h)
    k, l = minmax(k, l)
    b, e = minmax(b, e)
    g, m = minmax(g, m)

    h, i = minmax(h, i)
    l, m = minmax(l, m)
    e, j = minmax(e, j)
    g, k = minmax(g, k)

    d, e = minmax(d, e)
    f, g = minmax(f, g)
    i, j = minmax(i, j)
    k, l = minmax(k, l)
    b, h = minmax(b, h)

    c, g = minmax(c, g)
    j, l = minmax(j, l)
    b, d = minmax(b, d)
    e, h = minmax(e, h)
    i, k = minmax(i, k)
    a, f = minmax(a, f)

    c, f = minmax(c, f)
    g, i = minmax(g, i)
    j, k = minmax(j, k)

    b, c = minmax(b, c)
    d, f = minmax(d, f)
    h, i = minmax(h, i)
    e, g = minmax(e, g)

    c, d = minmax(c, d)
    e, f = minmax(e, f)
    g, h = minmax(g, h)
    i, j = minmax(i, j)

    d, e = minmax(d, e)
    f, g = minmax(f, g)

    return a,b,c,d,e,f,g,h,i,j,k,l,m
end

@inline swapsort(x::NTuple{13,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13])
@inline swapsort(::Type{Val{13}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13])

#    sort 14 values with 51 minmaxs in 10 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T) where {T}

    a, b = minmax(a, b)
    c, d = minmax(c, d)
    e, f = minmax(e, f)
    g, h = minmax(g, h)
    i, j = minmax(i, j)
    k, l = minmax(k, l)
    m, n = minmax(m, n)

    a, c = minmax(a, c)
    e, g = minmax(e, g)
    i, k = minmax(i, k)
    b, d = minmax(b, d)
    f, h = minmax(f, h)
    j, l = minmax(j, l)

    a, e = minmax(a, e)
    i, m = minmax(i, m)
    b, f = minmax(b, f)
    j, n = minmax(j, n)
    c, g = minmax(c, g)
    d, h = minmax(d, h)

    a, i = minmax(a, i)
    b, j = minmax(b, j)
    c, k = minmax(c, k)
    d, l = minmax(d, l)
    e, m = minmax(e, m)
    f, n = minmax(f, n)

    f, k = minmax(f, k)
    g, j = minmax(g, j)
    d, m = minmax(d, m)
    h, l = minmax(h, l)
    b, c = minmax(b, c)
    e, i = minmax(e, i)

    b, e = minmax(b, e)
    h, n = minmax(h, n)
    c, i = minmax(c, i)
    f, g = minmax(f, g)
    j, k = minmax(j, k)

    c, e = minmax(c, e)
    l, n = minmax(l, n)
    d, i = minmax(d, i)
    h, m = minmax(h, m)

    g, i = minmax(g, i)
    k, m = minmax(k, m)
    d, f = minmax(d, f)
    h, j = minmax(h, j)

    d, e = minmax(d, e)
    f, g = minmax(f, g)
    h, i = minmax(h, i)
    j, k = minmax(j, k)
    l, m = minmax(l, m)

    g, h = minmax(g, h)
    i, j = minmax(i, j)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n
end

@inline swapsort(x::NTuple{14,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14])
@inline swapsort(::Type{Val{14}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14])

#    sort 15 values with 56 minmaxs in 10 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T) where {T}

    a, b = minmax(a, b)
    c, d = minmax(c, d)
    e, f = minmax(e, f)
    g, h = minmax(g, h)
    i, j = minmax(i, j)
    k, l = minmax(k, l)
    m, n = minmax(m, n)

    a, c = minmax(a, c)
    e, g = minmax(e, g)
    i, k = minmax(i, k)
    m, o = minmax(m, o)
    b, d = minmax(b, d)
    f, h = minmax(f, h)
    j, l = minmax(j, l)

    a, e = minmax(a, e)
    i, m = minmax(i, m)
    b, f = minmax(b, f)
    j, n = minmax(j, n)
    c, g = minmax(c, g)
    k, o = minmax(k, o)
    d, h = minmax(d, h)

    a, i = minmax(a, i)
    b, j = minmax(b, j)
    c, k = minmax(c, k)
    d, l = minmax(d, l)
    e, m = minmax(e, m)
    f, n = minmax(f, n)
    g, o = minmax(g, o)

    f, k = minmax(f, k)
    g, j = minmax(g, j)
    d, m = minmax(d, m)
    n, o = minmax(n, o)
    h, l = minmax(h, l)
    b, c = minmax(b, c)
    e, i = minmax(e, i)

    b, e = minmax(b, e)
    h, n = minmax(h, n)
    c, i = minmax(c, i)
    l, o = minmax(l, o)
    f, g = minmax(f, g)
    j, k = minmax(j, k)

    c, e = minmax(c, e)
    l, n = minmax(l, n)
    d, i = minmax(d, i)
    h, m = minmax(h, m)

    g, i = minmax(g, i)
    k, m = minmax(k, m)
    d, f = minmax(d, f)
    h, j = minmax(h, j)

    d, e = minmax(d, e)
    f, g = minmax(f, g)
    h, i = minmax(h, i)
    j, k = minmax(j, k)
    l, m = minmax(l, m)

    g, h = minmax(g, h)
    i, j = minmax(i, j)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o
end

@inline swapsort(x::NTuple{15,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14], x[15])
@inline swapsort(::Type{Val{15}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14], x[15])

#    sort 16 values with 60 minmaxs in 10 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T) where {T}

    a, b = minmax(a, b)
    c, d = minmax(c, d)
    e, f = minmax(e, f)
    g, h = minmax(g, h)
    i, j = minmax(i, j)
    k, l = minmax(k, l)
    m, n = minmax(m, n)
    o, p = minmax(o, p)

    a, c = minmax(a, c)
    e, g = minmax(e, g)
    i, k = minmax(i, k)
    m, o = minmax(m, o)
    b, d = minmax(b, d)
    f, h = minmax(f, h)
    j, l = minmax(j, l)
    n, p = minmax(n, p)

    a, e = minmax(a, e)
    i, m = minmax(i, m)
    b, f = minmax(b, f)
    j, n = minmax(j, n)
    c, g = minmax(c, g)
    k, o = minmax(k, o)
    d, h = minmax(d, h)
    l, p = minmax(l, p)

    a, i = minmax(a, i)
    b, j = minmax(b, j)
    c, k = minmax(c, k)
    d, l = minmax(d, l)
    e, m = minmax(e, m)
    f, n = minmax(f, n)
    g, o = minmax(g, o)
    h, p = minmax(h, p)

    f, k = minmax(f, k)
    g, j = minmax(g, j)
    d, m = minmax(d, m)
    n, o = minmax(n, o)
    h, l = minmax(h, l)
    b, c = minmax(b, c)
    e, i = minmax(e, i)

    b, e = minmax(b, e)
    h, n = minmax(h, n)
    c, i = minmax(c, i)
    l, o = minmax(l, o)
    f, g = minmax(f, g)
    j, k = minmax(j, k)

    c, e = minmax(c, e)
    l, n = minmax(l, n)
    d, i = minmax(d, i)
    h, m = minmax(h, m)

    g, i = minmax(g, i)
    k, m = minmax(k, m)
    d, f = minmax(d, f)
    h, j = minmax(h, j)

    d, e = minmax(d, e)
    f, g = minmax(f, g)
    h, i = minmax(h, i)
    j, k = minmax(j, k)
    l, m = minmax(l, m)

    g, h = minmax(g, h)
    i, j = minmax(i, j)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p
end

@inline swapsort(x::NTuple{16,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14], x[15], x[16])
@inline swapsort(::Type{Val{16}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14], x[15], x[16])


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

    ThermalConsistencyOutput.jl

    Author: Tom Price
    Date:   July 2019

===============================================================================#

struct ThermalConsistencyOutput
    tm_check            ::Vector{TmCheck1w}
    delta_Tm            ::Tuple{Float_T, Bool}
    valid               ::Bool
end

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
path = @__DIR__

loadbutton = filepicker()
columnbuttons = Observable{Any}(dom"div"())
data = Observable{Any}(Phantom)
map!(f->begin
        global recon = JLD2.load(FileIO.File{FileIO.DataFormat{:JLD2}}(f),"recon")
        print("Loaded! ... $f\n")
    end
    , data, loadbutton)

# Ploting recon
global recon = ifftc(kdata)

hh, ww = 420,550
# l = PlotlyJS.Layout(;title="Reconstruction", yaxis_title="y", yaxis=attr(scaleanchor="x"),
#     xaxis_title="x",height=hh,width=ww,
#     modebar=attr(orientation="v"),scene=attr(aspectratio=attr(x=1,y=1,z=1)))
l = PlotlyJS.Layout(;title="Reconstruction",yaxis_title="y",
    xaxis_title="x",height=hh,width=ww,
    yaxis=attr(scaleanchor="x"),
    modebar=attr(orientation="v"),xaxis=attr(constrain="domain"),hovermode="closest")
p = PlotlyJS.plot(PlotlyJS.heatmap(z=abs.(recon),showscale=false,colorscale="Greys",transpose = false),l)
plt = Observable{Any}(p)
# PlotlyJS.savefig(p, path*"/assets/phantom.png", width=320, height=300)

recon = dom"div"(plt)
content!(w, "div#content", recon)


text:
module BenchCollisionModels

using BenchmarkTools
using LatticeBoltzmann

suite = BenchmarkGroup()

function initialize_benchmark(q = D2Q9(), τ = 1.0, scale = 2)
    benchmark_problem = LatticeBoltzmann.TGV(q, τ, scale)

    LatticeBoltzmann.initialize(
        LatticeBoltzmann.ZeroVelocityInitialCondition(),
        q,
        benchmark_problem,
    )
end

for q in [D2Q9(), D2Q13(), D2Q17(), D2Q21(), D2Q37()]
    suite[string(q)] = BenchmarkGroup([string(q), "collision_models"])

    srt_equilibrium = SRT(1.0)
    trt_equilibrium = TRT(1.0, 1.0)
    mrt_equilibrium = MRT(q, 1.0)

    τ = 1.0
    f = initialize_benchmark(q, τ)

    suite[string(q)]["srt", 1.0] = @benchmarkable(
        collide!(cm, $q, f_in, f_out),
        setup = (cm = $srt_equilibrium; f_in = copy($f); f_out = copy($f))
    )

    suite[string(q)]["srt", 1.3] = @benchmarkable(
        collide!(cm, $q, f_in, f_out),
        setup = (cm = SRT(1.3); f_in = copy($f); f_out = copy($f))
    )

    suite[string(q)]["srt", 0.8] = @benchmarkable(
        collide!(cm, $q, f_in, f_out),
        setup = (cm = SRT(0.8); f_in = copy($f); f_out = copy($f))
    )

    suite[string(q)]["trt", (1.0, 1.0)] = @benchmarkable(
        collide!(cm, $q, f_in, f_out),
        setup = (cm = $trt_equilibrium; f_in = copy($f); f_out = copy($f))
    )

    suite[string(q)]["trt", (0.9, 1.1)] = @benchmarkable(
        collide!(cm, $q, f_in, f_out),
        setup = (cm = TRT(0.9, 1.1); f_in = copy($f); f_out = copy($f))
    )

    suite[string(q)]["mrt-equilibrium", 1.0] = @benchmarkable(
        collide!(cm, $q, f_in, f_out),
        setup = (cm = $mrt_equilibrium; f_in = copy($f); f_out = copy($f))
    )

    scale = 2
    benchmark_problem = LatticeBoltzmann.DecayingShearFlow(
        1.0 / (2 * q.speed_of_sound_squared),
        scale,
        static = true,
    )
    f_force = LatticeBoltzmann.initialize(
        LatticeBoltzmann.ZeroVelocityInitialCondition(),
        q,
        benchmark_problem,
    )
    suite[string(q)]["srt", "force"] = @benchmarkable(
        collide!(cm, $q, f_in, f_out),
        setup = (cm = CollisionModel(SRT, $q, $benchmark_problem);
        f_in = copy($f_force);
        f_out = copy($f_force))
    )

    suite[string(q)]["trt", "force"] = @benchmarkable(
        collide!(cm, $q, f_in, f_out),
        setup = (cm = CollisionModel(TRT, $q, $benchmark_problem);
        f_in = copy($f_force);
        f_out = copy($f_force))
    )
end

end  # module
BenchCollisionModels.suite


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Models to simulate a power plant.
#*--------------------------------------------------------------------
#* Author: Argimiro R. Secchi
#* $Id: power_plant.mso 195 2007-03-07 20:30:12Z arge $
#*-------------------------------------------------------------------
# Declaracao de tipos
export Temperatura
typealias DanaTemperatura DanaRealParametric
type _Temperatura
	function _Temperatura(_::Dict{Symbol,Any})
		fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
		fields[:Default]=600
		fields[:Lower]=273.16
		fields[:Upper]=900
		fields[:Unit]="K"
		drive!(fields,_)
		new(fields)
	end
	value::Dict{Symbol,Any}
end
typealias Temperatura DanaRealParametric{_Temperatura}


text:
export Coulomb

"""
    Coulomb(; cutoff, nl_only, weight_14, coulomb_const, force_units, energy_units)

The Coulomb electrostatic interaction.
"""
struct Coulomb{C, W, T, F, E} <: GeneralInteraction
    cutoff::C
    nl_only::Bool
    weight_14::W
    coulomb_const::T
    force_units::F
    energy_units::E
end

const coulombconst = 138.93545764u"kJ * mol^-1 * nm" # 1 / 4πϵ0

function Coulomb(;
                    cutoff=NoCutoff(),
                    nl_only=false,
                    weight_14=1,
                    coulomb_const=coulombconst,
                    force_units=u"kJ * mol^-1 * nm^-1",
                    energy_units=u"kJ * mol^-1")
    return Coulomb{typeof(cutoff), typeof(weight_14), typeof(coulomb_const), typeof(force_units), typeof(energy_units)}(
        cutoff, nl_only, weight_14, coulomb_const, force_units, energy_units)
end

@inline @inbounds function force(inter::Coulomb{C},
                                    coord_i,
                                    coord_j,
                                    atom_i,
                                    atom_j,
                                    box_size,
                                    weight_14::Bool=false) where C
    dr = vector(coord_i, coord_j, box_size)
    r2 = sum(abs2, dr)

    cutoff = inter.cutoff
    coulomb_const = inter.coulomb_const
    qi, qj = atom_i.charge, atom_j.charge

    params = (coulomb_const, qi, qj)

    if cutoff_points(C) == 0
        f = force_divr_nocutoff(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        r2 > cutoff.sqdist_cutoff && return ustrip.(zero(coord_i)) * inter.force_units

        f = force_divr_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        r2 > cutoff.sqdist_cutoff && return ustrip.(zero(coord_i)) * inter.force_units

        if r2 < cutoff.activation_dist
            f = force_divr_nocutoff(inter, r2, inv(r2), params)
        else
            f = force_divr_cutoff(cutoff, r2, inter, params)
        end
    end

    if weight_14
        return f * dr * inter.weight_14
    else
        return f * dr
    end
end

@fastmath function force_divr_nocutoff(::Coulomb, r2, invr2, (coulomb_const, qi, qj))
    (coulomb_const * qi * qj) / √(r2 ^ 3)
end

@inline @inbounds function potential_energy(inter::Coulomb{C},
                                            coord_i,
                                            coord_j,
                                            atom_i,
                                            atom_j,
                                            box_size,
                                            weight_14::Bool=false) where C
    dr = vector(coord_i, coord_j, box_size)
    r2 = sum(abs2, dr)

    cutoff = inter.cutoff
    coulomb_const = inter.coulomb_const
    qi, qj = atom_i.charge, atom_j.charge
    params = (coulomb_const, qi, qj)

    if cutoff_points(C) == 0
        pe = potential(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        r2 > cutoff.sqdist_cutoff && return ustrip(zero(box_size[1])) * inter.energy_units

        pe = potential_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        r2 > cutoff.sqdist_cutoff && return ustrip(zero(box_size[1])) * inter.energy_units

        if r2 < cutoff.activation_dist
            pe = potential(inter, r2, inv(r2), params)
        else
            pe = potential_cutoff(cutoff, r2, inter, params)
        end
    end

    if weight_14
        return pe * inter.weight_14
    else
        return pe
    end
end

@fastmath function potential(::Coulomb, r2, invr2, (coulomb_const, qi, qj))
    (coulomb_const * qi * qj) * √invr2
end


text:
#######################################################################################################################################################################################################
#
# Changes made to this function
# General
#     2020-Mar-30: migrate the function from SCOPE
#     2021-Oct-21: rename the function to average_transmittance
#     2022-Feb-02: fix documentation
# To do
#     TODO: figure out where does this equation comes from
#
#######################################################################################################################################################################################################
"""

    average_transmittance(α::FT, nr::FT) where {FT<:AbstractFloat}

Return the average transmittance of isotropic radiation across an interface between two dielectrics, given
- `α` angle of incidence
- `nr` Index of refraction

# References
- Stern (1964) Transmission of isotropic radiation across an interface between two dielectrics. Applied Optics 3(1): 111-113.
- Allen (1973) Transmission of isotropic light across a dielectric surface in two and three dimensions. Journal of the Optical Society of America 63(6): 664-666.
"""
function average_transmittance(α::FT, nr::FT) where {FT<:AbstractFloat}
    @assert 0 < α <= 90;

    # some shortcuts to avoid overly comlicated equation
    _a     = (nr + 1) ^ 2 / 2;
    _a³    = _a ^ 3;
    _n²    = nr ^ 2;
    _n⁴    = nr ^ 4;
    _n⁶    = nr ^ 6;
    _n²p   = _n² + 1;
    _n²p²  = _n²p ^ 2;
    _n²p³  = _n²p ^ 3;
    _n²m²  = (_n² - 1) ^ 2;
    _k     = -1 * _n²m² / 4;
    _k²    = _k ^ 2;
    _sin²α = sind(α) ^ 2;
    _b₂    = _sin²α - _n²p/2;
    _b₁    = (α==90 ? 0 : sqrt(_b₂ ^ 2 + _k));
    _b     = _b₁ - _b₂;
    _b³    = _b ^ 3;
    _npanm = 2 * _n²p * _a - _n²m²;
    _npbnm = 2 * _n²p * _b - _n²m²;

    # S polarization
    _ts  = ( _k² / (6*_b³) + _k/_b - _b/2 ) - ( _k² / (6*_a³) + _k/_a - _a/2 );

    # P polarization
    _tp₁ = -2 * _n² * (_b - _a) / _n²p²;
    _tp₂ = -2 * _n² * _n²p * log(_b / _a) / _n²m²;
    _tp₃ = _n² * (1/_b - 1/_a) / 2;
    _tp₄ = 16 * _n⁴ * (_n⁴ + 1) * log(_npbnm / _npanm) / (_n²p³ * _n²m²);
    _tp₅ = 16 * _n⁶ * (1/_npbnm - 1/_npanm) / _n²p³;
    _tp  = _tp₁ + _tp₂ + _tp₃ + _tp₄ + _tp₅;

    return  (_ts + _tp) / (2 * _sin²α)
end


text:
# External packages
using Revise, DelimitedFiles, SymPy

# Add Exasim to Julia search path
cdir = pwd(); ii = findlast("Exasim", cdir);
include(cdir[1:ii[end]] * "/Installation/setpath.jl");

# Exasim packages
using Preprocessing, Mesh, Gencode, Postprocessing

pde = Array{Any, 1}(undef, 2);
mesh = Array{Any, 1}(undef, 2);

# create pde and mesh for each PDE model
include("pdeapp1.jl");
include("pdeapp2.jl");

# call exasim to generate and run C++ code to solve the PDE models
sol,pde,mesh,master,dmd,compilerstr,runstr = exasim(pde,mesh);


text:
# basisfunctions.jl

"""
    fourier(k::Int)

fourier implements the Fourier basis of functions ϕ_k, defined by ϕ_0≡1, and if k≥1 and k is odd, by
ϕ_k(x)=sqrt(2)sin((k+1)*π*x) and if k≥1 and even, by ϕ_k(x)=sqrt(2)cos(π*k*x).
This function is defined for k≥0.

# Examples

```julia
x = -1.0:0.01:2.0
y = fourier(3).(x)
```
"""
function fourier(k::Int)
  k ≥ 0 || throw(AssertionError("k should be non-negative"))
  if k == 0
    return x -> 1.0
  elseif isodd(k)
    kplusonefloat = float(k+1)
    return x -> 1.4142135623730951 * sinpi(kplusonefloat*x)
  else
    kfloat = float(k)
    return x -> 1.4142135623730951 * cospi(kfloat*x)
  end
end

"""
    faberschauderone(x)

Implements the first Faber-Schauder function defined by 1-2x for 0≤x≤1/2 and
-1+2x for 0.5≤x≤1, and is 1-periodically extended to all x∈R.

#See also: faberschauder

#Warning

Note the difference between faberschauderone and faberschauder. The first is a
function that takes a Float64 and returns a Float64, the second takes (j,k) and
returns an anonymous function that takes a Float64 and returns a Float64.

#Examples

```julia
x=-2.0:0.001:2.0
y=faberschauderone.(x)
```
"""
function faberschauderone(x)
    y = mod(x, 1.0)
    if 0≤y≤0.5
        return 1.0 - 2*y
    else
        return -1.0 + 2*y
    end
end

"""
    faberschauder(j::Int, k::Int)

faberschauder implements the k-th Faber-Schauder function of level j. Here, j≥0 and 1≤k≤2^j.
It is a one-periodic function and defined on [0,1] by 2^(j+1)(x-(k-1)2^(-j)) on
(k-1)2^(-j)≤x≤(k-1/2)2^(-j) and 1 - 2^(j+1)(x-(k-1/2)2^(-j)) on
[(k-1/2)2^(-j), k2^(-j)] and zero outside these intervals.

# See also: faberschauderone.

#Warning

Note the difference between faberschauderone and faberschauder. The first is a
function that takes a Float64 and returns a Float64, the second takes (j,k) and
returns an anonymous function that takes a Float64 and returns a Float64.

#Example with Plots
```julia
using Plots
J=2
x=0.0:0.001:1.0
p=plot()
for j in 0:J
  for k in 1:2^j
      y = faberschauder(j,k).(x)
      plot!(p,x,y)
  end
end
display(p)
```
"""
function faberschauder(j::Int,k::Int)
  j ≥ 0 || throw(AssertionError("j should be a nonnegative integer."))
  1≤k≤2^j || throw(AssertionError("k should be an integer between 1 and 2^j."))
  return function(x::Float64)
    y = mod(x,1.0)
    if y ≤ ldexp(float(k-1), -j) || y ≥ ldexp(float(k), -j)
      return 0.0
    elseif ldexp(float(k-1), -j) < y ≤ ldexp(k-0.5, -j)
      return ldexp(y, j+1)-2*(k-1)
    else
      return 2*k - ldexp(y, j+1)
    end
  end
end


text:
module TestTesting

import SplittablesTesting

SplittablesTesting.test_ordered(Any[1:10, 1:11])
SplittablesTesting.test_unordered(Any[Set(1:10), Set(1:11)])

end  # module


text:
using JuMP
using GLPK
using GraphPlot
using LightGraphs
using Colors
using LinearAlgebra

G = [
    0 1 0 0 0 0;
    1 0 1 1 0 0;
    0 1 0 0 1 1;
    0 1 0 0 1 0;
    0 0 1 1 0 0;
    0 0 1 0 0 0;
]

g = SimpleGraph(G)

gplot(g)

vertex_cover = Model(GLPK.Optimizer)

@variable(vertex_cover, y[1:nv(g)], Bin)
@constraint(vertex_cover, [i = 1:nv(g), j = 1:nv(g); G[i,j] == 1], y[i] + y[j] >= 1)
@objective(vertex_cover, Min, sum(y))

optimize!(vertex_cover)
@show value.(y);

membership = round.(Int, value.(y)) # Change to Int
membership = membership + ones(Int, nv(g)) # Make the color groups one indexed
nodecolor = [colorant"red", colorant"blue"] # Blue to represent vertices in the cover
nodefillc = nodecolor[membership]
gplot(g, nodefillc = nodefillc)

G = [
    0 1 0 0 0 0 0 0 0 1 0 ;
    1 0 1 0 0 0 0 0 0 0 1;
    0 1 0 1 0 1 0 0 0 0 0;
    0 0 1 0 1 0 0 0 0 0 0;
    0 0 0 1 0 1 0 0 0 0 0;
    0 0 1 0 1 0 1 0 0 0 0;
    0 0 0 0 0 1 0 1 0 0 0;
    0 0 0 0 0 0 1 0 1 0 1;
    0 0 0 0 0 0 0 1 0 1 1;
    1 0 0 0 0 0 0 0 1 0 1;
    0 1 0 0 0 0 0 1 1 1 0;
]

g = SimpleGraph(G)

gplot(g)

dominating_set = Model(GLPK.Optimizer)

@variable(dominating_set, x[1:nv(g)], Bin)
@constraint(dominating_set, [i = 1:nv(g)], dot(G[i,:], x) >= 1)
@objective(dominating_set, Min, sum(x))

optimize!(dominating_set)
@show value.(x);

membership = convert(Array{Int},value.(x)) # Change to Int
membership = membership + ones(Int, nv(g)) # Make the color groups one indexed
nodecolor = [colorant"red", colorant"blue"] # Blue to represent vertices in the set
nodefillc = nodecolor[membership]
gplot(g, nodefillc = nodefillc)

G = [
    0 0 0 0 1 0 0 0;
    0 0 0 0 0 1 0 0;
    0 0 0 0 0 0 1 0;
    0 0 0 0 0 0 0 1;
    1 0 0 0 0 1 0 1;
    0 1 0 0 1 0 1 0;
    0 0 1 0 0 1 0 1;
    0 0 0 1 1 0 1 0;
]

g = SimpleGraph(G)

gplot(g)

matching = Model(GLPK.Optimizer)

@variable(matching, m[i = 1:nv(g), j = 1:nv(g)], Bin)
@constraint(matching, [i = 1:nv(g)], sum(m[i,:]) <= 1)
@constraint(matching, [i = 1:nv(g), j = 1:nv(g); G[i,j] == 0], m[i,j] == 0)
@constraint(matching, [i = 1:nv(g), j = 1:nv(g)], m[i,j] == m[j,i])
@objective(matching, Max, sum(m))

optimize!(matching)
@show value.(m);

G = [
    0 1 0 0 1 1 0 0 0 0;
    1 0 1 0 0 0 1 0 0 0;
    0 1 0 1 0 0 0 1 0 0;
    0 0 1 0 1 0 0 0 1 0;
    1 0 0 1 0 0 0 0 0 1;
    1 0 0 0 0 0 1 0 0 1;
    0 1 0 0 0 1 0 1 0 0;
    0 0 1 0 0 0 1 0 1 0;
    0 0 0 1 0 0 0 1 0 1;
    0 0 0 0 1 1 0 0 1 0;
]

g = SimpleGraph(G)

gplot(g)

k = nv(g)

k_colouring = Model(GLPK.Optimizer)

@variable(k_colouring, z[1:k], Bin)
@variable(k_colouring, c[1:nv(g),1:k], Bin)
@constraint(k_colouring, [i = 1:nv(g)], sum(c[i,:]) == 1)
@constraint(k_colouring, [i = 1:nv(g), j = 1:nv(g), l = 1:k; G[i,j] == 1], c[i,l] + c[j,l] <= 1)
@constraint(k_colouring, [i = 1:nv(g), l = 1:k], c[i,l] <= z[l])

@objective(k_colouring, Min, sum(z))

optimize!(k_colouring)
@show value.(z);
@show value.(c);

c = value.(c)
membership = zeros(nv(g))
for i in 1:nv(g)
    for j in 1:k
        if c[i,j] == 1
            membership[i] = j
            break
        end
    end
end
membership = convert(Array{Int},membership)

nodecolor = distinguishable_colors(nv(g), colorant"green")
nodefillc = nodecolor[membership]
gplot(g, nodefillc = nodefillc)


text:
# This file illustrates the calculation of energy spectral density and
# total energy of a complex exponential pulse.

using FFTW
using Plots

# Construct the signal
f = 20
fs = 100
l = 2^10
ts = 1 / fs
t = (0 : l -1) * ts
x = exp.(1im * 2π * f * t)

# Compute fft
xf = fft(x) * ts
ff = 0 : fs / (l - 1) : fs

# Compute energy
esd = abs.(xf).^2
Δf = fs / (l - 1)
Ef = sum(esd) * Δf              # Energy calculated in time
Et = sum(abs.(x).^2) * ts       # Energy calculated in frequency
@show Ef, Et

# Plots
plt = plot(layout=1)
plot!(ff, esd, subplot=1)
xlabel!("Frequency [Hz]")
ylabel!("Ψ(f)")
title!("Energy Spectal Density")
display(plt)


text:
module ConstantCellFields

using Gridap
import Gridap: gradient

const ConstantCellFieldLike{D,T,N,R<:FieldLike{D,T,N}} = ConstantCellValue{R}

function gradient(f::ConstantCellFieldLike)
  v = f.value
  g = gradient(v)
  ConstantCellValue(g,f.length)
end

end # module


text:
module PartitionedPoissonTests

using PartitionedPoisson
using Test

include("../compile/gridap_warmup.jl")

@time @testset "seq 2d" begin poisson(mode=:seq,nc=(4,4),np=(2,2)) end

end # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "gdk_pixbuf"
version = v"2.38.2" # we are actually on master

# Collection of sources required to build gdk-pixbuf
sources = [
    ArchiveSource("https://gitlab.gnome.org/GNOME/gdk-pixbuf/-/archive/3c7740498fd31b6746dd7e04601886766a6644b7/gdk-pixbuf-3c7740498fd31b6746dd7e04601886766a6644b7.tar.bz2",
                  "9fad057e8c51bc4373948a02c8ee7d8afe254b361fb4abc43767fce43982dd25"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gdk-pixbuf-*/
mkdir build && cd build

# Get a local gettext for msgfmt cross-building
apk add gettext

FLAGS=()
if [[ "${target}" == *-apple-* ]] || [[ "${target}" == *-mingw* ]]; then
    FLAGS+=(-Dx11=false)
fi

meson .. \
    -Dgir=false \
    -Dman=false \
    -Dinstalled_tests=false \
    -Dgio_sniffing=false \
    "${FLAGS[@]}" \
    --cross-file="${MESON_TARGET_TOOLCHAIN}"
ninja -j${nproc}
ninja install

# Cleanup `loaders.cache` file, we're going to generate a new one on the user's machine
rm -f ${prefix}/lib/gdk-pixbuf-2.0/2.10.0/loaders/loaders.cache
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libgdk_pixbuf-2", "libgdk_pixbuf-2.0"], :libgdkpixbuf),
    ExecutableProduct("gdk-pixbuf-query-loaders", :gdk_pixbuf_query_loaders),
    FileProduct("lib/gdk-pixbuf-2.0/2.10.0/loaders", :gdk_pixbuf_loaders_dir),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Glib_jll", v"2.59.0"; compat="2.59.0"),
    Dependency("JpegTurbo_jll"),
    Dependency("libpng_jll"),
    # TODO: v4.3.0 is available, use that next time
    Dependency("Libtiff_jll"; compat="4.1.0"),
    Dependency("Xorg_libX11_jll"),
    BuildDependency("Xorg_xproto_jll"),
    BuildDependency("Xorg_kbproto_jll"),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: remi-garcia/Yggdrasil

max_issues_repo_name: remi-garcia/Yggdrasil

max_forks_repo_name: remi-garcia/Yggdrasil


text:
function stratification(r,Nmax,rbottom,iopt)

# defines stratification on radial grid.

# r       - radial grid points
# Nmax    - stratification at CMB (r=1)
# Rbottom - radius of base of stratified layer
# iopt =  - constant stratification (=0); linear stratification (=1)

#   constructs density stratification
   n = length(r);
   N2 = zeros(n);

   N2_peak = Nmax * Nmax;

#  stratified layer
   if (iopt == 1)
       N2[1:n] .= (r[1:n] .- r[1]) * N2_peak / (r[n] - r[1]);
   end

   if (iopt == 0 )
       N2[1:n] .= N2_peak;
   end

   return N2

end

max_stars_repo_name: bbuffett/MacWaves

max_issues_repo_name: bbuffett/MacWaves

max_forks_repo_name: bbuffett/MacWaves


text:
using GPUArrays
import CuArrays.CUBLAS

function batched_scal!(A::CuArray{T, 1}, B::CuArray{T, 3}) where T
    gpu_call(B, (A, B)) do state, A, B
        batch = @linearidx(A)
        s = zero(eltype(B))
        for i in axes(B, 1), j in axes(B, 2)
            B[i, j, batch] *= A[batch]
        end
        return
    end
    return B
end

for (gemm, elty) in
        ((:dgemm_,:Float64),
         (:sgemm_,:Float32),
         (:zgemm_,:ComplexF64),
         (:cgemm_,:ComplexF32))

    @eval begin
        function batched_gemm!(transA::AbstractChar, transB::AbstractChar, alpha::($elty), A::CuArray{$elty, 3}, B::CuArray{$elty, 3}, beta::($elty), C::CuArray{$elty, 3})
            CUBLAS.gemm_strided_batched!(transA, transB, alpha, A, B, beta, C)
        end

        function batched_gemm(transA::AbstractChar, transB::AbstractChar, alpha::($elty), A::CuArray{$elty, 3}, B::CuArray{$elty, 3})
            CUBLAS.gemm_strided_batched(transA, transB, alpha, A, B)
        end

        function batched_gemm(transA::AbstractChar, transB::AbstractChar, A::CuArray{$elty, 3}, B::CuArray{$elty, 3})
            CUBLAS.gemm_strided_batched(transA, transB, A, B)
        end
    end
end