text:
function format(doc::WeaveDoc)
    formatted = AbstractString[]
    docformat = doc.format

    #Complete format dictionaries with defaults
    formatdict = docformat.formatdict
    get!(formatdict, :termstart, formatdict[:codestart])
    get!(formatdict, :termend, formatdict[:codeend])
    get!(formatdict, :out_width, nothing)
    get!(formatdict, :out_height, nothing)
    get!(formatdict, :fig_pos, nothing)
    get!(formatdict, :fig_env, nothing)

    for chunk in copy(doc.chunks)
        result = format_chunk(chunk, formatdict, docformat)
        push!(formatted, result)
    end

    return formatted
end

function format_chunk(chunk::DocChunk, formatdict, docformat)
    return chunk.content
end

function format_chunk(chunk::CodeChunk, formatdict, docformat)
    #Fill undefined options with format specific defaults
    chunk.options[:out_width] == nothing &&
        (chunk.options[:out_width] =  formatdict[:out_width])
    chunk.options[:fig_pos] == nothing &&
        (chunk.options[:fig_pos] =  formatdict[:fig_pos])

    #Only use floats if chunk has caption or sets fig_env
    if chunk.options[:fig_cap] != nothing && chunk.options[:fig_env] == nothing
        (chunk.options[:fig_env] =  formatdict[:fig_env])
    end

    if haskey(formatdict, :indent)
        chunk.content = indent(chunk.content, formatdict[:indent])
    end

    if !chunk.options[:eval]
        if chunk.options[:echo]
            result = "$(formatdict[:codestart])$(chunk.content)\n$(formatdict[:codeend])"
            return result
        else
            r = ""
            return r
        end
    end

    if chunk.options[:term]
        result = format_termchunk(chunk, formatdict)
    else

    if chunk.options[:echo]
        result = "$(formatdict[:codestart])$(chunk.content)\n$(formatdict[:codeend])\n"
    else
        result = ""
    end

    if (strip(chunk.output)!= "" || strip(chunk.rich_output) != "") && (chunk.options[:results] != "hidden")
        if chunk.options[:results] != "markup" && chunk.options[:results] != "hold"
            strip(chunk.output) ≠ "" && (result *= "$(chunk.output)\n")
            strip(chunk.rich_output) ≠ "" && (result *= "$(chunk.rich_output)\n")
        else
            if chunk.options[:wrap]
                chunk.output = "\n" * wraplines(chunk.output,
                                        chunk.options[:line_width])
            end

            if haskey(formatdict, :indent)
                chunk.output = indent(chunk.output, formatdict[:indent])
            end
            strip(chunk.output) ≠ "" &&
                (result *= "$(formatdict[:outputstart])$(chunk.output)\n$(formatdict[:outputend])\n")
            strip(chunk.rich_output) ≠ "" && (result *= chunk.rich_output * "\n")
        end
    end

    end

    #Handle figures
    if chunk.options[:fig] && length(chunk.figures) > 0
        if chunk.options[:include]
            result *= formatfigures(chunk, docformat)
        end
    end

    return result
end

function format_termchunk(chunk, formatdict)
    if chunk.options[:echo] && chunk.options[:results] != "hidden"
        result = "$(formatdict[:termstart])$(chunk.output)\n" * "$(formatdict[:termend])\n"
        #chunk.options[:term_state] == :text && (result*= "$(formatdict[:termend])\n")
    else
        result = ""
    end
    return result
end

function indent(text, nindent)
    return join(map(x->
                    string(repeat(" ", nindent), x), split(text, "\n")), "\n")
end

function wraplines(text, line_width=75)
    result = AbstractString[]
    lines = split(text, "\n")
    for line in lines
        if length(line) > line_width
            push!(result, wrapline(line, line_width))
        else
            push!(result, line)
        end
    end

    #return result
    return strip(join(result, "\n"))
end

function wrapline(text, line_width=75)
result = ""
    while length(text) > line_width
        result*= text[1:line_width] * "\n"
        text = text[(line_width+1):end]
    end
result *= text
end

type Tex
    description::AbstractString
    formatdict::Dict{Symbol,Any}
end

const tex = Tex("Latex with custom code environments",
                @compat Dict{Symbol,Any}(:codestart => "\\begin{juliacode}",
                                         :codeend => "\\end{juliacode}",
                                         :outputstart => "\\begin{juliaout}",
                                         :outputend => "\\end{juliaout}",
                                         :termstart => "\\begin{juliaterm}",
                                         :termend => "\\end{juliaterm}",
                                         :fig_ext => ".pdf",
                                         :extension =>"tex",
                                         :out_width=> "\\linewidth",
                                         :fig_env=> "figure",
                                         :fig_pos => "htpb",
                                         :doctype => "tex",
                                         :mimetypes => ["application/pdf", "image/png", "text/latex", "text/plain"]
                                         ))

const texminted = Tex("Latex using minted for highlighting",
                      @compat Dict{Symbol,Any}(
                                         :codestart => "\\begin{minted}[mathescape, fontsize=\\small, xleftmargin=0.5em]{julia}",
                                         :codeend => "\\end{minted}",
                                         :outputstart => "\\begin{minted}[fontsize=\\small, xleftmargin=0.5em, mathescape, frame = leftline]{text}",
                                         :outputend => "\\end{minted}",
                                         :termstart=> "\\begin{minted}[fontsize=\\footnotesize, xleftmargin=0.5em, mathescape]{julia}",
                                         :termend => "\\end{minted}",
                                         :fig_ext => ".pdf",
                                         :extension =>"tex",
                                         :out_width => "\\linewidth",
                                         :fig_env=> "figure",
                                         :fig_pos => "htpb",
                                         :doctype => "texminted",
                                         :mimetypes => ["application/pdf", "image/png", "text/latex", "text/plain"]
                                         ))

type Pandoc
  description::AbstractString
  formatdict::Dict{Symbol,Any}
end

const pandoc = Pandoc("Pandoc markdown",
                        @compat Dict{Symbol,Any}(
                                :codestart => "~~~~{.julia}",
                                :codeend=>"~~~~~~~~~~~~~\n\n",
                                :outputstart=>"~~~~",
                                :outputend=>"~~~~\n\n",
                                :fig_ext=>".png",
                                :out_width=>nothing,
                                :extension=>"md",
                                #Prefer png figures for markdown conversion, svg doesn't work with latex
                                :mimetypes => ["image/png", "image/jpg", "image/svg+xml", "text/markdown", "text/plain"],
                                :doctype=>"pandoc"
                                               ))

const md2html = Pandoc("Markdown to HTML (requires Pandoc)",
                      @compat Dict{Symbol,Any}(
                              :codestart => "````julia",
                              :codeend=> "````\n\n",
                              :outputstart=> "````",
                              :outputend=> "````\n\n",
                              :fig_ext=> ".svg",
                              :extension=> "md",
                              :mimetypes => ["image/svg+xml", "image/png", "image/jpg",
                                  "text/html", "text/markdown",  "text/plain"],
                              :doctype=> "md2html"))

const md2pdf = Pandoc("Markdown to pdf (requires Pandoc and xelatex)",
                      @compat Dict{Symbol,Any}(
                              :codestart => "````julia",
                              :codeend=> "````\n\n",
                              :outputstart=> "````",
                              :outputend=> "````\n\n",
                              :fig_ext=> ".pdf",
                              :extension=> "md",
                              :mimetypes => ["application/pdf", "image/png", "image/jpg",
                                  "text/latex", "text/plain"],
                              :doctype=> "md2pdf"))

type Markdown
   description::AbstractString
   formatdict::Dict{Symbol,Any}
end

const github = Markdown("Github markdown",
                        @compat Dict{Symbol,Any}(
                                :codestart => "````julia",
                                :codeend=> "````\n\n",
                                :outputstart=> "````",
                                :outputend=> "````\n\n",
                                :fig_ext=> ".png",
                                :extension=> "md",
                                :doctype=> "github"
                                               ))

type MultiMarkdown
  description::AbstractString
  formatdict::Dict{Symbol,Any}
end

const multimarkdown = MultiMarkdown("MultiMarkdown",
                        @compat Dict{Symbol,Any}(
                                :codestart => "````julia",
                                :codeend=> "````\n\n",
                                :outputstart=> "````",
                                :outputend=> "````\n\n",
                                :fig_ext=> ".png",
                                :extension=> "md",
                                :doctype=> "github"
                                               ))

type Rest
    description::AbstractString
    formatdict::Dict{Symbol,Any}
end

const rst = Rest("reStructuredText and Sphinx",
                 @compat Dict{Symbol,Any}(
                                :codestart => ".. code-block:: julia\n",
                                :codeend => "\n\n",
                                :outputstart => "::\n",
                                :outputend => "\n\n",
                                :indent=> 4,
                                :fig_ext => ".png",
                                :extension => "rst",
                                :out_width => "15 cm",
                                :doctype => "rst"
                                ))

type AsciiDoc
    description::AbstractString
    formatdict::Dict{Symbol,Any}
end

#asciidoc -b html5 -a source-highlighter=pygments ...
const adoc = AsciiDoc("AsciiDoc",
        @compat Dict{Symbol,Any}(
        :codestart => "[source,julia]\n--------------------------------------",
        :codeend => "--------------------------------------\n\n",
        :outputstart => "--------------------------------------",
        :outputend => "--------------------------------------\n\n",
        :fig_ext => ".png",
        :extension => "txt",
        :out_width => "600",
        :doctype => "asciidoc"
))

function formatfigures(chunk, docformat::Tex)
    fignames = chunk.figures
    caption = chunk.options[:fig_cap]
    width = chunk.options[:out_width]
    height = chunk.options[:out_height]
    f_pos = chunk.options[:fig_pos]
    f_env = chunk.options[:fig_env]
    result = ""
    figstring = ""

    #Set size
    attribs = ""
    width == nothing || (attribs = "width=$width")
    (attribs != "" && height != nothing ) && (attribs *= ",")
    height == nothing    || (attribs *= "height=$height")

    if f_env != nothing
        result *= """\\begin{$f_env}[$f_pos]\n"""
    end

    for fig = fignames

        if splitext(fig)[2] == ".tex" #Tikz figures
            figstring *= "\\resizebox{$width}{!}{\\input{$fig}}\n"
        else
            figstring *= "\\includegraphics[$attribs]{$fig}\n"
        end
    end

    # Figure environment
    if caption != nothing
        result *= string("\\center\n",
                         "$figstring",
                         "\\caption{$caption}\n")
    else
        result *= figstring
    end

    if chunk.options[:name] != nothing && f_env !=nothing
        label = chunk.options[:name]
        result *= "\\label{fig:$label}\n"
    end

    if f_env != nothing
        result *= "\\end{$f_env}\n"
    end

   return result
end

function formatfigures(chunk, docformat::Pandoc)
    fignames = chunk.figures
    caption = chunk.options[:fig_cap]
    result = ""
    figstring = ""
    attribs = ""
    width = chunk.options[:out_width]
    height = chunk.options[:out_height]

    #Build figure attibutes
    width == nothing || (attribs = "width=$width")
    (attribs ≠ "" && height ≠ nothing ) && (attribs *= " ")
    height == nothing   || (attribs *= "height=$height")
    attribs == ""    || (attribs = "{$attribs}")
    length(fignames) > 0 || (return "")

    if caption != nothing
        result *= "![$caption]($(fignames[1]))$attribs\n"
        for fig = fignames[2:end]
            result *= "![]($fig)$attribs\n"
            println("Warning, only the first figure gets a caption\n")
        end
    else
        for fig in fignames
            result *= "![]($fig)$attribs\\ \n\n"
        end
    end
    return result
end

function formatfigures(chunk, docformat::Markdown)
    fignames = chunk.figures
    caption = chunk.options[:fig_cap]
    result = ""
    figstring = ""

    length(fignames) > 0 || (return "")

    if caption != nothing
        result *= "![$caption]($(fignames[1]))\n"
        for fig = fignames[2:end]
            result *= "![]($fig)\n"
            println("Warning, only the first figure gets a caption\n")
        end
    else
        for fig in fignames
            result *= "![]($fig)\n"
        end
    end
    return result
end

function formatfigures(chunk, docformat::MultiMarkdown)
    fignames = chunk.figures
    caption = chunk.options[:fig_cap]
    result = ""
    figstring = ""

    if chunk.options[:out_width] == nothing
      width = ""
    else
      width = "width=$(chunk.options[:out_width])"
    end

    length(fignames) > 0 || (return "")

    if caption != nothing
       result *= "![$caption][$(fignames[1])]\n\n"
       result *= "[$(fignames[1])]: $(fignames[1]) $width\n"
        for fig = fignames[2:end]
          result *= "![][$fig]\n\n"
          result *= "[$fig]: $fig $width\n"
          println("Warning, only the first figure gets a caption\n")
        end
    else
        for fig in fignames
          result *= "![][$fig]\n\n"
          result *= "[$fig]: $fig $width\n"
        end
    end
    return result
end

function formatfigures(chunk, docformat::Rest)
    fignames = chunk.figures
    caption = chunk.options[:fig_cap]
    width = chunk.options[:out_width]
    result = ""
    figstring = ""

    for fig=fignames
        figstring *= @sprintf(".. image:: %s\n   :width: %s\n\n", fig, width)
    end

    if caption != nothing
        result *= string(".. figure:: $(fignames[1])\n",
                         "   :width: $width\n\n",
                         "   $caption\n\n")
    else
        result *= figstring
        return result
    end
end

function formatfigures(chunk, docformat::AsciiDoc)
    fignames = chunk.figures
    caption = chunk.options[:fig_cap]
    width = chunk.options[:out_width]
    result = ""
    figstring = ""

    for fig=fignames
        figstring *= @sprintf("image::%s[width=%s]\n", fig, width)
    end

    if caption != nothing
        result *= string("image::$(fignames[1])",
        "[width=$width,",
        "title=\"$caption\"]")
    else
        result *= figstring
        return result
    end
end

#Add new supported formats here
const formats = @compat Dict{AbstractString, Any}("tex" => tex,
                                          "texminted" => texminted,
                                          "pandoc" => pandoc,
                                          "md2html" => md2html,
                                          "md2pdf" => md2pdf,
                                          "github" => github,
                                          "multimarkdown" => multimarkdown,
                                          "rst" => rst,
                                          "asciidoc" => adoc
                                          )


text:
module Layers
export Layer, next, top_layer, insert

include("exceptions.jl")

"""
## Request Execution Stack

The Request Execution Stack is separated into composable layers.

Each layer is defined by a nested type `Layer{Next}` where the `Next`
parameter defines the next layer in the stack.
The `request` method for each layer takes a `Layer{Next}` type as
its first argument and dispatches the request to the next layer
using `request(Next, ...)`.

The example below defines three layers and three stacks each with
a different combination of layers.

```julia
abstract type Layer end
abstract type Layer1{Next <: Layer} <: Layer end
abstract type Layer2{Next <: Layer} <: Layer end
abstract type Layer3 <: Layer end

request(::Type{Layer1{Next}}, data) where Next = "L1", request(Next, data)
request(::Type{Layer2{Next}}, data) where Next = "L2", request(Next, data)
request(::Type{Layer3}, data) = "L3", data

const stack1 = Layer1{Layer2{Layer3}}
const stack2 = Layer2{Layer1{Layer3}}
const stack3 = Layer1{Layer3}
```

```julia
julia> request(stack1, "foo")
("L1", ("L2", ("L3", "foo")))

julia> request(stack2, "bar")
("L2", ("L1", ("L3", "bar")))

julia> request(stack3, "boo")
("L1", ("L3", "boo"))
```

This stack definition pattern gives the user flexibility in how layers are
combined but still allows Julia to do whole-stack compile time optimisations.

e.g. the `request(stack1, "foo")` call above is optimised down to a single
function:
```julia
julia> code_typed(request, (Type{stack1}, String))[1].first
CodeInfo(:(begin
    return (Core.tuple)("L1", (Core.tuple)("L2", (Core.tuple)("L3", data)))
end))
```
"""
abstract type Layer{Next} end

"""
    next(::Type{S}) where {T, S<:Layer{T}}

Return the next `Layer` in the stack

Example:
stack = MessageLayer{ConnectionPoolLayer{StreamLayer{Union{}}}
next(stack)  # ConnectionPoolLayer{StreamLayer{Union{}}}
"""
next(::Type{S}) where {T, S<:Layer{T}} = T

"""
    top_layer(::Type{T}) where T <: Layer

Return the parametric type of the top most `Layer` in the stack

Example:
stack = MessageLayer{ConnectionPoolLayer{StreamLayer{Union{}}}
top_layer(stack)  # MessageLayer
"""
top_layer(::Type{T}) where T <: Layer = T.name.wrapper
top_layer(::Type{Union{}}) = Union{}

"""
    insert(stack::Type{<:Layer}, layer_before::Type{<:Layer}, custom_layer::Type{<:Layer})

Insert your `custom_layer` in-front of the `layer_before`

Example:
stack = MessageLayer{ConnectionPoolLayer{StreamLayer{Union{}}}
result = insert(stack, MessageLayer, TestLayer)  # TestLayer{MessageLayer{ConnectionPoolLayer{StreamLayer{Union{}}}}}
"""
function insert(stack::Type{<:Layer}, layer_before::Type{<:Layer}, custom_layer::Type{<:Layer})
    new_stack = Union
    head_layer = top_layer(stack)
    rest_stack = stack

    while true
        if head_layer === layer_before
            return new_stack{custom_layer{rest_stack}}
        else
            head_layer === Union{} && break
            new_stack = new_stack{head_layer{T}} where T
            rest_stack = next(rest_stack)
            head_layer = top_layer(rest_stack)
        end
    end
    throw(LayerNotFoundException("$layer_before not found in $stack"))
end

end

ext: jl


text:
using Base.Test
using PropertyGraph
using UUID
using Compat

# run the test within a function scope to avoid interference with other tests

# define custom types required for the test
type TestCustomPersonType
	firstname::String
	lastname::String
end

type TestKnowsType
	someproperty::String
end

function graphpopulationtest()
	# create a graph
	g = Graph()

	# create a vertex and add it to the graph
	fred = add!(g, Vertex("Person", @compat Dict{String,Any}("name"=>"Fred","age"=>45)))
	# and test that it was initialised as expected
	# it should belong to the graph
	@test fred.graph == g
	# and have the property values passed in on initialisation
	@test get(fred,"name") == "Fred"
	@test get(fred,"age") == 45

	setpropertyvalue!(fred, "other", 100)
	@test get(fred,"other") == 100

	setpropertyvalues!(fred, @compat Dict{String,Any}("other"=>101,"new"=>"new value"))
	@test get(fred,"other") == 101
	@test get(fred,"new") == "new value"

	# retrieval of a property not set should return UnspecifiedValue
	@test get(fred,"otherpropertynotset") == UnspecifiedValue
	# unless a default is passed in to the get call in which case the default should be returned
	@test get(fred,"otherpropertynotset",0) == 0

	# attempting to add the same vertex again should result in error
	caughtexception = UnspecifiedValue
	try
		add!(g, fred)
	catch ex
		caughtexception = ex
	end

	@test isa(caughtexception, VertexAlreadyBelongsToSpecifiedGraphException)

	# make a new graph
	g2 = Graph()

	caughtexception = UnspecifiedValue
	try
		add!(g2, fred)
	catch ex
		caughtexception = ex
	end
	@test isa(caughtexception, VertexAlreadyBelongsToAnotherGraphException)

	# make some more vertices
	sally = add!(g, Vertex("Person", @compat Dict{String,Any}("name"=>"Sally","age"=>25)))

	johanobject = TestCustomPersonType("Johan", "Smith")
	johan = vertexforobject("Person", johanobject,  @compat Dict{String,Any}("name"=>"Johan","age"=>28))
	@test johan.containedobject == johanobject
	add!(g, johan)

	# the graph should now have 3 vertices
	@test length(g.vertices) == 3

	# create some edges
	e = Edge(fred,sally)
	e2 = Edge("knows", fred, sally, @compat Dict{String,Any}("since"=>2004))
	@test get(e2,"since") == 2004

	#including one with an associated object
	fredknowsjohanobject = TestKnowsType("SomeValue")
	e3 = edgeforobject("knows", fredknowsjohanobject, fred,johan)
	@test e3.containedobject == fredknowsjohanobject

	# add the edges to the graph
	add!(g, e)
	add!(g, e2)
	add!(g, e3)

	#the graph should now have 3 edges
	@test length(g.edges) == 3

	# test the edges are in the expected sets
	@test in(e, fred.outgoingedges)
	@test in(e2, fred.outgoingedges)
	@test in(e3, fred.outgoingedges)
	@test length(fred.outgoingedges) == 3

	@test in(e, sally.incomingedges)
	@test in(e2, sally.incomingedges)
	@test length(sally.incomingedges) == 2

	@test in(e3, johan.incomingedges)
	@test length(johan.incomingedges) == 1

	# test adding an edge already in the graph
	# this should result in exception
	caughtexception = UnspecifiedValue
	try
		add!(g, e)
	catch ex
		caughtexception = ex
	end
	@test isa(caughtexception, EdgeAlreadyBelongsToSpecifiedGraphException)

	# test adding an edge to another graph
	# this should result in exception
	caughtexception = UnspecifiedValue
	try
		add!(g2, e)
	catch ex
		caughtexception = ex
	end
	@test isa(caughtexception, EdgeAlreadyBelongsToAnotherGraphException)

	# create another vertex but don't add it to the graph
	unaddedVertex = Vertex("Person")

	# create another edge involving the new vertex as tail
	e4 = Edge(unaddedVertex,sally)

	# attempt adding the edge to the graph involving a vertex with tail not in the graph
	caughtexception = UnspecifiedValue
	try
		add!(g, e4)
	catch ex
		caughtexception = ex
	end

	@test isa(caughtexception, EdgeTailDoesNotBelongToGraphException)

	# create another edge involving the new vertex as head
	e5 = Edge(sally, unaddedVertex)

	# attempt adding the edge to the graph involving a vertex with head not in the graph
	caughtexception = UnspecifiedValue
	try
		add!(g, e5)
	catch ex
		caughtexception = ex
	end
	@test isa(caughtexception, EdgeHeadDoesNotBelongToGraphException)
end
graphpopulationtest()


text:
using Test
using TuringTutorials

@testset "TuringTutorials" begin
    include("build.jl")
end

max_stars_repo_name: staticfloat/TuringTutorials

max_issues_repo_name: staticfloat/TuringTutorials

max_forks_repo_name: staticfloat/TuringTutorials


text:
# In this script, we will design a PID controller by optimization. The system model has uncertain parameters, and we pay a price not only for poor performance of the closed-loop system, but also for a high variance in the performance. In addition to this, we place a constraint on the 90:th percentile of the maximum of the sensitivity function. This way, we will get a doubly robust controller as a result :)
# We start by defining the system and some initial controller parameters
using MonteCarloMeasurements, Optim, ControlSystems, Plots
using MonteCarloMeasurements: ∓
p = 1 ∓ 0.1
ζ = 0.3 ∓ 0.05
ω = 1 ∓ 0.1
P = tf([p*ω], [1, 2ζ*ω, ω^2])
w = exp10.(LinRange(-2,3,100))
params = [1,0.1,0.1]
Msc = 1.5 # Constraint on Ms

# We now define the cost function, which includes the constraint on the maximum sensitivity function
function cost(params)
    kp,ki,kd = params
    C = pid(kp=kp,ki=ki,kd=kd)
    G = feedback(P*C) # Closed-loop system
    S = 1/(1 + P*C)   # Sensitivity function
    local Gd
    try
        @unsafe Gd = c2d(G,0.1) # Discretize the system. This might fail for some parameters, so we catch these cases and return a high value
    catch
        return 1000
    end
    y  = step(Gd,15)[1][:] # This is the time-domain simulation
    Ms = maximum(bode(S, w)[1]) # Maximum of the sensitivity function
    q  = quantile(Ms, 0.9)
    performance = mean(abs, 1 .- y)   # This is our performance measure
    robustness = (q > Msc ? 10000(q-Msc) : 0) # This is our robustness constraint
    variance = 10std(performance)     # This is the price we pay for high variance in the performance
    mean(performance) + robustness + variance
end

# We are now ready to test the cost function. This will take very long time to compile the first time it's called since we use StaticParticles (~60s on my machine), but should be very fast after that (~6ms)
@time cost(params)
#
res = Optim.optimize(cost, params, SimulatedAnnealing(), Optim.Options(iterations=200, show_trace=false));
println("Final cost: ", res.minimum)
# We can now perform the same computations as above to visualize the found controller
fig = plot(layout=2)
for params = (params, res.minimizer)
    kp,ki,kd = params
    C  = pid(kp =kp,ki =ki,kd =kd)
    G  = feedback(P*C)
    S  = 1/(1 + P*C)
    @unsafe Gd = c2d(G,0.1)
    y,t,_ = step(Gd,15)
    y     = y[:]
    mag   = bode(S, w)[1][:]
    plot!(t,y, title="Time response", subplot=1, legend=false)
    plot!(w, mag, title="Sensitivity function", xscale=:log10, yscale=:log10, subplot=2, legend=false)
end
hline!([Msc], l=(:black, :dash), subplot=2)
fig

# Other things that could potentially be relevant is adding a probabilistic constraint on the time-domain output, such as the probability of having the step response go above 1.5 must be < 0.05 etc.


text:
push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/movie/")
# push!(LOAD_PATH, "/n/home07/iczekala/DiskJockey/")
# push!(LOAD_PATH, "/pool/scout0/DiskJockey/")

using consonance

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    # "--opt1"
    # help = "an option with an argument"
    "--run_index", "-r"
    help = "Output run index"
    default = 0
    arg_type = Int
end

parsed_args = parse_args(ARGS, s)
run_index = parsed_args["run_index"]

imgdir = moviedir * @sprintf("g%02d/", run_index)

using DiskJockey.constants
using DiskJockey.image

import YAML
config = YAML.load(open("config.yaml"))
dpc = config["parameters"]["dpc"][1]

import PyPlot.plt
using LaTeXStrings

global norm = plt[:Normalize](vmin=log10(vmax_g) - 4, vmax=log10(vmax_g), clip=false)

# Plot the channel maps using sky convention
function plot_chmaps(img::image.SkyImage, id::Int)

    fig, ax = plt[:subplots](ncols=nchan, figsize=(6.4, 1.6))
    # Image needs to be flipped along RA dimension
    ext = (img.ra[end], img.ra[1], img.dec[1], img.dec[end])

    for col=1:nchan

        ax[col][:xaxis][:set_ticklabels]([])
        ax[col][:yaxis][:set_ticklabels]([])

        # Set all labels but the leftmost to blank
        # if col != 1
        #     ax[col][:xaxis][:set_ticklabels]([])
        #     ax[col][:yaxis][:set_ticklabels]([])
        # else
        #     ax[col][:set_xlabel](L"$\Delta \alpha$ ('')", size=12)
        #     ax[col][:set_ylabel](L"$\Delta \delta$ ('')", size=12)
        #     ax[col][:tick_params](axis="both", which="major", labelsize=10)
        # end

        #Flip the frame for Sky convention
        frame = flipdim(img.data[:,:,col], 2)
        frame += 1e-10 #Add a tiny bit so that we don't have log10(0)
        lframe = log10(frame)
        max = maximum(lframe)
        ax[col][:imshow](lframe, extent=ext, norm=norm, interpolation="none", origin="lower", cmap=plt[:get_cmap]("PuBu"))
        if col==1
            ax[col][:annotate](@sprintf("%.1f km/s", vels[col]), (0.1, 0.78), xycoords="axes fraction", size=12)
        else
            ax[col][:annotate](@sprintf("%.1f", vels[col]), (0.1, 0.78), xycoords="axes fraction", size=12)
        end

    end

    # Annotate the parameters
    p = pars[id]
    mass = p.M_star
    r_c = p.r_c
    incl = p.incl
    label_incl = L"$i$: " * @sprintf("%.0f", incl) * L"${}^\circ$"
    label_mass = L"$M_\ast$: " * @sprintf("%.2f", mass) * L" $M_\odot$"
    fig[:text](0.01, 0.87, label_mass, size=16)
    fig[:text](0.25, 0.87, label_incl, size=16)

    fig[:subplots_adjust](wspace=0.08, top=0.9, bottom=0.01, left=0.02, right=0.98)

    #plt.savefig(outdir * @sprintf("g%04d.png", id), dpi=150)
    plt[:savefig](imgdir * @sprintf("g%04d.svg", id), transparent=true)
    plt[:close]("all")

end

# How many frames per process?
start = run_index * nframes_per_proc + 1

ids = Int[i for i=start:(start + nframes_per_proc)]

for id in ids
    if id <= nframes
        fname = imgdir * @sprintf("gimage%04d.out", id)
        im = imread(fname)
        skim = imToSky(im, dpc)
        plot_chmaps(skim, id)
    else
        break
    end
end

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
#################
### User Info ###
#################

"""
Gitea User
"""
FieldTags.@tag immutable User
	id::Int64
	userName::String => json:"login"
	fullName::String => json:"full_name"
	email::String => json:"email"
	avatarURL::String => json:"avatar_url"
end

##########################
### Application Tokens ###
##########################

"""
AccessToken represents a API access token.
"""
immutable AccessToken
	name::String
	sha1::String
end

####################
### User Email ###
####################

immutable Email
	email::String
	verified::Bool
	primary::Bool
end

####################
### User GPG Key ###
####################

immutable GPGKeyEmail
	email::String
	verified::Bool
end

FieldTags.@tag immutable GPGKey
	id::Int64
	primaryKeyID::String => json:"primary_key_id"
	keyID::String => json:"key_id"
	publicKey::String => json:"public_key"
	emails::Vector{GPGKeyEmail}
	subkeys::Vector{GPGKey}
	canSign::Bool => json:"can_sign"
	canEncryptComms::Bool => json:"can_encrypt_comms"
	canEncryptStorage::Bool => json:"can_encrypt_storage"
	canCertify::Bool => json:"can_certify"
	created::Nullable{DateTime} => json:"created_at,format:y-m-dTH:M:SZ"
	expires::Nullable{DateTime} => json:"expires_at,format:y-m-dTH:M:SZ"
end

####################
### User SSH Key ###
####################

FieldTags.@tag immutable PublicKey
	id::Int64
	key::String
	url::Nullable{String}
	title::Nullable{String}
	create::Nullable{DateTime} => json:"created_at,format:y-m-dTH:M:SZ"
end

##############
### Issues ###
##############

const StateType = String
const StateOpen = "open"
const StateClosed = "closed"

FieldTags.@tag immutable PullRequestMeta
	hasMerged::Bool => json:"merged"
	merged::DateTime => json:"merged_at,format:y-m-dTH:M:SZ"
end

FieldTags.@tag immutable Label
	id::Int64
	name::String
	color::String
	url::String
end

FieldTags.@tag immutable Milestone
	id::Int64
	title::String
	description::String
	state::StateType
	openIssues::Int => json:"open_issues"
	closedIssues::Int => json:"closed_issues"
	closed::Nullable{DateTime} => json:"closed_at,format:y-m-dTH:M:SZ"
	deadline::Nullable{DateTime} => json:"due_on,format:y-m-dTH:M:SZ"
end

FieldTags.@tag immutable Issue
	id::Int64
	url::String
	index::Int64 => json:"number"
	poster::User => json:"user"
	title::String
	body::String
	labels::Vector{Label}
	milestone::Nullable{Milestone}
	assignee::Nullable{User}
	state::StateType
	comments::Int
	created::DateTime => json:"created_at,format:y-m-dTH:M:SZ"
	updated::DateTime => json:"updated_at,format:y-m-dTH:M:SZ"

	pullRequest::Nullable{PullRequestMeta} => json:"pull_request"
end

FieldTags.@tag immutable Comment
	id::Int64
	htmlURL::String => json:"html_url"
	prURL::String => json:"pull_request_url"
	issueURL::String => json:"issue_url"
	poster::Nullable{User}
	body::String
	created::DateTime => json:"created_at,format:y-m-dTH:M:SZ"
	updated::DateTime => json:"updated_at,format:y-m-dTH:M:SZ"
end

####################
### Repositories ###
####################

immutable Permission
	admin::Bool
	push::Bool
	pull::Bool
end

FieldTags.@tag immutable Repository
	id::Int64
	owner::User
	name::String
	fullName::String => json:"full_name"
	description::String
	private::Bool
	fork::Bool
	mirror::Bool
	htmlURL::String => json:"html_url"
	sshURL::String => json:"ssh_url"
	cloneURL::String => json:"clone_url"
	website::String
	stars::Int => json:"stars_count"
	forks::Int => json:"forks_count"
	watchers::Int => json:"watchers_count"
	openIssues::Int => json:"open_issues_count"
	defaultBranch::String => json:"default_branch"
	created::DateTime => json:"created_at,format:y-m-dTH:M:SZ"
	updated::DateTime => json:"updated_at,format:y-m-dTH:M:SZ"
	permissions::Nullable{Permission}
end

##################
### Hook Types ###
##################

FieldTags.@tag immutable Hook
	id::Int64
	hookType::String => json:"type"
	url::Nullable{String}
	config::Dict{String,String}
	events::Vector{String}
	active::Bool
	created::DateTime => json:"created_at,format:y-m-dTH:M:SZ"
	updated::DateTime => json:"updated_at,format:y-m-dTH:M:SZ"
end

#####################
### Hook Payloads ###
#####################

immutable UserPayload
	name::String
	email::String
	username::String
end

immutable CommitVerificationPayload
	verified::Bool
	reason::String
	signature::String
	payload::String
end

FieldTags.@tag immutable CommitPayload
	id::String
	message::String
	url::String
	author::UserPayload
	committer::UserPayload
	verification::Nullable{CommitVerificationPayload}
	timestamp::DateTime => json:",format:y-m-dTH:M:SZ"
end

#########################
### Repository Branch ###
#########################

immutable Branch
	name::String
	commit::CommitPayload
end

#############################
### Repository Deploy Key ###
#############################

FieldTags.@tag immutable DeployKey
	id::Int64
	key::String
	url::String
	title::String
	created::DateTime => json:"created_at,format:y-m-dTH:M:SZ"
	readonly::Bool => json:"read_only"
end

####################
### Organization ###
####################

FieldTags.@tag immutable Organization
	id::Int64
	username::String
	fullName::String => json:"full_name"
	avatarURL::String => json:"avatar_url"
	description::String
	website::String
	location::String
end

immutable Team
	id::Int64
	name::String
	description::String
	permission::String
end

####################
### Pull Request ###
####################

FieldTags.@tag immutable PRBranchInfo
	name::String
	ref::String
	sha::String
	repoID::Int64 => json:"repo_id"
	repo::Nullable{Repository}
end

FieldTags.@tag immutable PullRequest
	id::Int64
	url::String
	index::Int64
	poster::Nullable{User} => json:"user"
	title::String
	body::String
	labels::Vector{Label}
	milestone::Nullable{Milestone}
	assignee::Nullable{User}
	state::StateType
	comments::Int

	htmlURL::String => json:"html_url"
	diffURL::String => json:"diff_url"
	patchURL::String => json:"patch_url"

	mergeable::Bool
	hasMerged::Bool
	merged::Nullable{DateTime} => json:"merged_at,format:y-m-dTH:M:SZ"
	mergedCommitID::Nullable{String} => json:"merge_commit_sha"
	mergedBy::Nullable{User} => json:"merged_by"

	base::Nullable{PRBranchInfo}
	head::Nullable{PRBranchInfo}
	mergeBase::String => json:"merge_base"
end

###############
### Release ###
###############

FieldTags.@tag immutable Release
	id::Int64
	tagName::String => json:"name"
	target::String => json:"target_commitish"
	title::String => json:"name"
	note::String => json:"body"
	url::String
	tarURL::String => json:"tarball_url"
	zipURL::String => json:"zipball_url"
	isDraft::Bool => json:"draft"
	isPrerelease::Bool => json:"prelease"
	createdAt::DateTime => json:"created_at,format:y-m-dTH:M:SZ"
	publishedAt::DateTime => json:"published_at,format:y-m-dTH:M:SZ"
	publisher::User => json:"author"
end

##############
### Status ###
##############

const StatusState = String
const StatusPending = "pending"
const StatusSuccess = "success"
const StatusError   = "error"
const StatusFailure = "failure"
const StatusWarning = "warning"

FieldTags.@tag immutable Status
	id::Int64
	state::StatusState
	targetURL::String => json:"target_url"
	description::String
	url::String
	context::String
	creator::User
	created::DateTime => json:"created_at,format:y-m-dTH:M:SZ"
	updated::DateTime => json:"updated_at,format:y-m-dTH:M:SZ"
end

FieldTags.@tag immutable CombinedStatus
	state::StatusState
	sha::String
	totalCount::Int => json:"total_count"
	statuses::Vector{Status}
	repository::Repository
	commitURL::String => json:"commit_url"
	url::String
end


text:
using Makie

meshscatter(rand(100), rand(100), rand(100))

AbstractPlotting.current_scene()

scatter(
    Point3f0[(1,0,0), (0,1,0), (0,0,1)],
    marker = [:x, :circle, :cross]
)

s = heatmap(rand(100, 100))
Makie.center!(s)
s
scene = contour(-1..1, -1..1, rand(100, 100))
scene

surface(rand(100, 100))
surface(-1..1, -1..1, rand(100, 100))
surface(LinRange(-1, 1, 100), range(-1, stop = 1, length = 100) .+ (rand(100) .* 0.1), rand(100, 100))
N = 40
r = range(-1, stop = 1, length = N)
x = map(hcat(repeated(r, N)...)) do x
   x + (rand() * 0.01)
end
y = vcat(repeated(r', N)...)
surface(x, y, rand(N, N))
surface(x, y, rand(N, N) * 0.1, color = rand(RGBAf0, N, N))

wf = wireframe(LinRange(-1, 1, 100), range(-1, stop = 1, length = 100) .+ (rand(100) .* 0.1), rand(100, 100))
wf = wireframe(-1..1, -1..1, rand(100, 100))
wf = wireframe(Sphere(Point3f0(0), 1f0))
wf = wireframe(x, y, rand(N, N))

scene = Scene()
x = map([:dot, :dash, :dashdot], [2, 3, 4]) do ls, lw
    linesegments!(scene, range(1, stop = 5, length = 100), rand(100), rand(100), linestyle = ls, linewidth = lw)
end
push!(x, scatter!(scene, range(1, stop = 5, length = 100), rand(100), rand(100)))
scene

l = Makie.legend(x, ["attribute $i" for i in 1:4])
io = VideoStream(scene, homedir()*"/Desktop")
record(io) = (for i = 1:35; recordframe!(io); sleep(1/30); end);
l[:position] = (0, 1)
record(io)
l[:backgroundcolor] = RGBA(0.95, 0.95, 0.95)
record(io)
l[:strokecolor] = RGB(0.8, 0.8, 0.8)
record(io)
l[:gap] = 30
record(io)
l[:textsize] = 19
record(io)
l[:linepattern] = Point2f0[(0,-0.2), (0.5, 0.2), (0.5, 0.2), (1.0, -0.2)]
record(io)
l[:scatterpattern] = decompose(Point2f0, Circle(Point2f0(0.5, 0), 0.3f0), 9)
record(io)
l[:markersize] = 2f0
record(io)
finish(io, "mp4")

using Makie, Colors, GeometryTypes
GLAbstraction.DummyCamera <: GLAbstraction.Camera
scene = Scene()
cmap = collect(LinRange(to_color(:red), to_color(:blue), 20))
l = Makie.legend(cmap, 1:4)
l[:position] = (1.0,1.0)

l[:textcolor] = :blue
l[:strokecolor] = :black
l[:strokewidth] = 1
l[:textsize] = 15
l[:textgap] = 5

using Makie
scene = Scene(resolution = (500, 500))
x = [0, 1, 2, 0]
y = [0, 0, 1, 2]
z = [0, 2, 0, 1]
color = [:red, :green, :blue, :yellow]
i = [0, 0, 0, 1]
j = [1, 2, 3, 2]
k = [2, 3, 1, 3]
indices = [1, 2, 3, 1, 3, 4, 1, 4, 2, 2, 3, 4]
m = mesh(x, y, z, indices, color = color)
m[:color] = [:blue, :red, :blue, :red]

using Makie, GeometryTypes

scene = Scene()
sub = Scene(scene, offset = Vec3f0(1, 2, 0))
scatter(sub, rand(10), rand(10), camera = :orthographic)
sub[:camera]
lines(sub, rand(10), rand(10), camera = :orthographic)
axis(LinRange(0, 2, 4), range(0, stop = 2, length = 4))
center!(scene)
sub[:offset] = Vec3f0(0, 0, 0)

using Makie, GeometryTypes

scene = Scene()

# create a subscene from which the next scenes you will plot
# `inherit` attributes. Need to add the attribute you care about,
# in this case the rotation
sub = Scene(scene, rotation = Vec4f0(0, 0, 0, 1))

meshscatter(sub, rand(30) + 1.0, rand(30), rand(30))
meshscatter(sub, rand(30), rand(30), rand(30) .+ 1.0)
r = range(-2, stop = 2, length = 4)
Makie.axis(r, r, r)
center!(scene)

axis = Vec3f0(0, 0, 1)
io = VideoStream(scene, homedir()*"/Desktop/", "rotation")
for angle = range(0, stop = 2pi, length = 100)
    sub[:rotation] = Makie.qrotation(axis, angle)
    recordframe!(io)
    sleep(1/15)
end
finish(io, "gif")

keys = (
    :lol, :pos, :test
)
positions, pos, test = get.(attributes, keys, scene)

function myvisual(scene, args, attributes)
    keys = (
        :positions, :color, :blah, shared...,
    )
    Scene(zip(keys, getindex.(attributes, keys))

end

using Makie, GeometryTypes

function plot(scene::S, A::AbstractMatrix{T}) where {T <: AbstractFloat, S <: Scene}
    N, M = size(A)
    sub = Scene(scene, scale = Vec3f0(1))
    attributes = Dict{Symbol, Any}()

    plots = map(1:M) do i
        lines(sub, 1:N, A[:, i])
    end
    labels = get(attributes, :labels) do
        map(i-> "y $i", 1:M)
    end

    lift_node(to_node(A), to_node(Makie.getscreen(scene).area)) do a, area
        xlims = (1, size(A, 1))
        ylims = extrema(A)
        stretch = Makie.to_nd(fit_ratio(area, (xlims, ylims)), Val{3}, 1)
        sub[:scale] = stretch
    end
    l = legend(scene, plots, labels)
    # Only create one axis per scene
    xlims = range(1, stop = size(A, 1), length = min(N, 10))
    a, b = extrema(A)
    ylims = range(a, stop = b, length = 5)
    a = get(scene, :axis) do
        xlims = (1, size(A, 1))
        ylims = extrema(A)
        area = Reactive.value(Makie.getscreen(scene).area)
        stretch = Makie.to_nd(fit_ratio(area, (xlims, ylims)), Val{3}, 1)
        axis(LinRange((xlims .* stretch[1])..., 4), LinRange((ylims .* stretch[2])..., 4))
    end
    center!(scene)
end

using Makie, GeometryTypes

scene = Scene()

x = (0.5rel, 0.5rel)

x .* Point2f0(0.5, 0.5)
StaticArrays.similar_type(NTuple{2, Float32}, Int)

Makie.VecTypes{3, Float32}
Makie.Units.to_absolute(scene, x)
to_positions(scene, (rand(10) .* rel, rand(10) .* rel))

convert(Relative{Float64}, 1)

convert(NTuple{2, Relative{Float64}}, Tuple(Relative.(widths(scene))))

plot(scene, rand(11, 2))


text:
export PARSDMM_precompute_distribute

"""
Precomputes and distributes some quantities that serve as input for PARSDMM.jl
"""
function PARSDMM_precompute_distribute(
                                      TD_OP    ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}},
                                      set_Prop,
                                      comp_grid,
                                      options
                                      ) where {TF<:Real,TI<:Integer}

N = prod(comp_grid.n)

#add the identity matrix as the operator for the squared distance from the point we want to project
if options.feasibility_only==false
  push!(TD_OP,SparseMatrixCSC{TF}(LinearAlgebra.I,N,N));
  push!(set_Prop.TD_n,comp_grid.n)
  push!(set_Prop.AtA_offsets,[0])
  push!(set_Prop.banded,true)
  push!(set_Prop.AtA_diag,true)
  push!(set_Prop.ncvx,false)
  push!(set_Prop.dense,false)
  push!(set_Prop.tag,("distance squared","identity","matrix",""))
end

p     = length(TD_OP);

AtA=Vector{SparseMatrixCSC{TF,TI}}(undef,p)
for i=1:p
  if set_Prop.dense[i]==true
    if set_Prop.AtA_diag[i]==true
      AtA[i]=SparseMatrixCSC{TF}(LinearAlgebra.I,N,N)
    else
      error("provided a dense non orthogoal transform-domain operator")
    end
  else
    AtA[i] = TD_OP[i]'*TD_OP[i]
  end
end

#if all AtA are banded -> convert to CDS (DIA) format
if sum(set_Prop.banded[1:p].=true)==p
  for i=1:p
    (AtA[i],set_Prop.AtA_offsets[i])  = mat2CDS(AtA[i])
    set_Prop.AtA_offsets[i]=convert(Vector{TI},set_Prop.AtA_offsets[i])
  end
   AtA=convert(Vector{Array{TF,2}},AtA);
  set_Prop.AtA_offsets=set_Prop.AtA_offsets[1:p]
end

#allocate arrays of vectors
y       = Vector{Vector{TF}}(undef,p);
l       = Vector{Vector{TF}}(undef,p);

for ii=1:p #initialize all rho's, gamma's, y's and l's
    y[ii]       = zeros(TF,size(TD_OP[ii],1))
    l[ii]       = zeros(TF,length(y[ii]));
end

if options.parallel==true
  TD_OP = distribute(TD_OP)
  y     = distribute(y)
  l     = distribute(l)
end

return TD_OP,AtA,l,y
end


text:
@testset "8.string-to-integer-atoi.jl" begin
    @test my_atoi("   -42") == -42
    @test my_atoi("") == 0
    @test my_atoi("4193 with words") == 4193
    @test my_atoi("words and 987") == 0
    @test my_atoi("-91283472332") == -2147483648
    @test my_atoi("+-0123") == 0
    @test my_atoi("91283472332") == typemax(Int32)
end


text:
@assert ARGS == ["a", "b"]
@assert Base.JLOptions().quiet == 1 # --quiet
@assert Base.JLOptions().check_bounds == 2 # overriden default when testing


text:
using IDRsSolver

println("FDM discretisation of a 3D convection-diffusion-reaction problem on a unit cube")
println("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")

# Define system

# Defaults:
h = 0.025
eps1 = 1.
beta = [0/sqrt(5), 250/sqrt(5), 500/sqrt(5)]
r = 400

# Generate matrix
m = round(Int, 1/h)-1
n = m*m*m
Sx = sparse(Tridiagonal((-eps1/h^2-beta[1]/(2*h))*ones(m-1),2*eps1/h^2*ones(m),(-eps1/h^2+beta[1]/(2*h))*ones(m-1)))
Sy = sparse(Tridiagonal((-eps1/h^2-beta[2]/(2*h))*ones(m-1),2*eps1/h^2*ones(m),-eps1/h^2+beta[2]/(2*h)*ones(m-1)))
Sz = sparse(Tridiagonal((-eps1/h^2-beta[3]/(2*h))*ones(m-1),2*eps1/h^2*ones(m),-eps1/h^2+beta[3]/(2*h)*ones(m-1)))
Is = speye(m)
A = kron(kron(Is,Is),Sx) + kron(kron(Is,Sy),Is)+ kron(kron(Sz,Is),Is) -r*speye(n)

x = linspace(h,1-h,m)
sol = kron(kron(x.*(1-x),x.*(1-x)),x.*(1-x))
b = A*sol

println(" ");
println("The parameters of the problem are :")
println("Gridsize h = ", (h),";")
println("Number of equations = ", (n),";")
println("Diffusion parameter = ", (eps1),";")
println("Convection parameters = (",(beta[1]),",",(beta[2]),",",(beta[3]),");")
println("Reaction parameter = ",(r)," (Note: positive reaction parameter gives negative shift to matrix);")
println(" ")

# Defaults for the iterative solvers:

tol = 1e-8
maxiter = 1000

s = 4
println("IDR(4) iteration...")
time = @elapsed x, _ = idrs( A, b; s=s, tol=tol, maxiter=maxiter)

println("Final accuracy: ", norm(b-A*x)/norm(b))
println("CPU time: ", time, "s.")


text:
struct IterateOutput
    index::Int
    dates::ClosedInterval{Date}
end

iter_index(x::IterateOutput) = x.index
iter_dates(x::IterateOutput) = x.dates

struct IterateTimelineTable{T, MNames, FNames, N1, N2}
    parent::MarketData{T, MNames, FNames, N1, N2}
    index_dict::Dict{T, Vector{IterateOutput}}
    key_vec::Vector{T}
    function IterateTimelineTable(data::MarketData{T, MNames, FNames, N1, N2}, index_dict, key_vec) where {T, MNames, FNames, N1, N2}
        @assert Set(key_vec) ⊆ Set(keys(index_dict)) "Some Keys are Missing"
        @assert Set(keys(index_dict)) ⊆ Set(keys(data.firmdata)) "Some Keys are not in the Parent Data"
        new{T, MNames, FNames, N1, N2}(data, index_dict, key_vec)
    end
end

parent(data::IterateTimelineTable) = data.parent

function construct_id_dict(
    ids::Vector{T},
    date_starts::Vector{Date},
    date_ends::Vector{Date},
) where {T}
    @assert length(ids) == length(date_starts) == length(date_ends) "All Vectors must be of same length"
    out = Dict{T, Vector{IterateOutput}}()
    sizehint!(out, length(ids))
    for (i, id) in enumerate(ids)
        if !haskey(out, id)
            out[id] = IterateOutput[]
        end
        push!(out[id], IterateOutput(i, date_starts[i] .. date_ends[i]))
    end
    out
end

construct_id_dict(ids::Vector{T}, dates::ClosedInterval{Vector{Date}}) where {T} = construct_id_dict(ids, dates.left, dates.right)

function Base.iterate(iter::IterateTimelineTable{T, MNames, FNames}, state=1) where {T, MNames, FNames}
    if state > length(iter.key_vec)
        return nothing
    end
    (iter[state], state+1)
end

function Base.eltype(::IterateTimelineTable{T, MNames, FNames, N1, N2}) where {T, MNames, FNames, N1, N2}
    Tuple{T, Vector{IterateOutput}}
end
function Base.length(iter::IterateTimelineTable)
    length(iter.key_vec)
end

total_length(data::IterateTimelineTable) = sum(length.(values(data.index_dict)))

function Base.getindex(data::IterateTimelineTable, i::Int)
    1 <= i <= length(data) || throw(BoundsError(data, i))
    id = data.key_vec[i]
    (id, data.index_dict[id])
end

Base.firstindex(data::IterateTimelineTable) = 1
Base.lastindex(data::IterateTimelineTable) = length(data)

function validate_iterator(
    data::IterateTimelineTable,
    out_vector::Vector
)
    all_ids = iter_index.(vcat(values(data.index_dict)...)) |> sort
    all_ids == 1:length(out_vector)
end

function Base.getindex(data::MarketData{T}, ids::Vector{T}, dates::ClosedInterval{Vector{Date}}) where {T}
    id_dict = construct_id_dict(ids, dates)
    IterateTimelineTable(
        data,
        id_dict,
        unique(ids)
    )
end

function Base.show(io::IO, data::IterateTimelineTable)
    println(io, "Iterable set of TimelineTable with $(length(data)) unique firms")
    println(io, "and a total number of iterations of $(total_length(data)) and a parent of")
    show(io, parent(data))
end


text:
import Statistics

@inline function marginal_median(a::AbstractVector{V}) where V <: AbstractVector{T} where T
    return Statistics.median.(a)
end


text:
function test_dropout_layer(backend::Backend, T, eps)
  println("-- Testing Dropout on $(typeof(backend)){$T}...")

  tensor_dim = abs(rand(Int)) % 6 + 1
  dims = tuple((abs(rand(Int, tensor_dim)) % 8 + 1)...)
  println("    > $dims")

  ratio = convert(T, 0.6)
  input = rand(T, dims)
  input_blob = make_blob(backend, input)
  diff_blob  = make_blob(backend, input)

  println("    > Setup")
  layer = DropoutLayer(bottoms=[:input], ratio=ratio)
  state = setup(backend, layer, Blob[input_blob], Blob[diff_blob])

  println("    > Forward")
  forward(backend, state, Blob[input_blob])
  got_output = zeros(T, size(input))
  copy!(got_output, input_blob)
  rand_vals = zeros(T, size(input))
  copy!(rand_vals, state.rand_vals)
  expected_output = input .* (rand_vals .> ratio) / (1-ratio)
  @test all(abs(got_output - expected_output) .< eps)

  println("    > Backward")
  top_diff = rand(T, size(input))
  copy!(diff_blob, top_diff)
  backward(backend, state, Blob[input_blob], Blob[diff_blob])
  expected_grad = top_diff .* (rand_vals .> ratio) / (1-ratio)
  got_grad = zeros(T, size(expected_grad))
  copy!(got_grad, diff_blob)
  @test all(abs(got_grad - expected_grad) .< eps)

  shutdown(backend, state)
end
function test_dropout_layer(backend::Backend)
  test_dropout_layer(backend, Float64, 1e-10)
  test_dropout_layer(backend, Float32, 1e-4)
end

if test_cpu
  test_dropout_layer(backend_cpu)
end
if test_gpu
  test_dropout_layer(backend_gpu)
end


text:
#!/usr/bin/env julia --project
using ClimateMachine
ClimateMachine.init()
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Grids: polynomialorder
using ClimateMachine.DGmethods: vars_state_conservative
using ClimateMachine.HydrostaticBoussinesq

using Test
using CLIMAParameters
using CLIMAParameters.Planet: grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

function config_simple_box(FT, N, resolution, dimensions; BC = nothing)
    if BC == nothing
        problem = HomogeneousBox{FT}(dimensions...)
    else
        problem = HomogeneousBox{FT}(dimensions...; BC = BC)
    end

    _grav::FT = grav(param_set)
    cʰ = sqrt(_grav * problem.H) # m/s
    model = HydrostaticBoussinesqModel{FT}(param_set, problem, cʰ = cʰ)

    config = ClimateMachine.OceanBoxGCMConfiguration(
        "homogeneous_box",
        N,
        resolution,
        model,
    )

    return config
end

function run_homogeneous_box(; imex::Bool = false, BC = nothing)
    FT = Float64

    # DG polynomial order
    N = Int(4)

    # Domain resolution and size
    Nˣ = Int(20)
    Nʸ = Int(20)
    Nᶻ = Int(20)
    resolution = (Nˣ, Nʸ, Nᶻ)

    Lˣ = 4e6   # m
    Lʸ = 4e6   # m
    H = 400   # m
    dimensions = (Lˣ, Lʸ, H)

    timestart = FT(0)    # s
    timeend = FT(6 * 3600) # s

    if imex
        solver_type = ClimateMachine.IMEXSolverType(
            linear_model = LinearHBModel,
            linear_solver = ClimateMachine.ColumnwiseLUSolver.SingleColumnLU,
        )
        Nᶻ = Int(400)
        Courant_number = 0.1
    else
        solver_type = ClimateMachine.ExplicitSolverType(
            solver_method = LSRK144NiegemannDiehlBusch,
        )
        Courant_number = 0.4
    end

    driver_config = config_simple_box(FT, N, resolution, dimensions; BC = BC)

    grid = driver_config.grid
    vert_filter = CutoffFilter(grid, polynomialorder(grid) - 1)
    exp_filter = ExponentialFilter(grid, 1, 8)
    modeldata = (vert_filter = vert_filter, exp_filter = exp_filter)

    solver_config = ClimateMachine.SolverConfiguration(
        timestart,
        timeend,
        driver_config,
        init_on_cpu = true,
        ode_solver_type = solver_type,
        modeldata = modeldata,
        Courant_number = Courant_number,
    )

    result = ClimateMachine.invoke!(solver_config)

    maxQ = Vars{vars_state_conservative(driver_config.bl, FT)}(maximum(
        solver_config.Q,
        dims = (1, 3),
    ))
    minQ = Vars{vars_state_conservative(driver_config.bl, FT)}(minimum(
        solver_config.Q,
        dims = (1, 3),
    ))

    @test maxQ.θ ≈ minQ.θ
end

@testset "$(@__FILE__)" begin
    boundary_conditions = [
        (
            OceanBC(Impenetrable(NoSlip()), Insulating()),
            OceanBC(Impenetrable(NoSlip()), Insulating()),
            OceanBC(Penetrable(KinematicStress()), Insulating()),
        ),
        (
            OceanBC(Impenetrable(FreeSlip()), Insulating()),
            OceanBC(Impenetrable(NoSlip()), Insulating()),
            OceanBC(Penetrable(KinematicStress()), Insulating()),
        ),
        (
            OceanBC(Impenetrable(NoSlip()), Insulating()),
            OceanBC(Impenetrable(FreeSlip()), Insulating()),
            OceanBC(Penetrable(KinematicStress()), Insulating()),
        ),
        (
            OceanBC(Impenetrable(FreeSlip()), Insulating()),
            OceanBC(Impenetrable(FreeSlip()), Insulating()),
            OceanBC(Penetrable(KinematicStress()), Insulating()),
        ),
    ]

    for BC in boundary_conditions
        run_homogeneous_box(imex = true, BC = BC)
        run_homogeneous_box(imex = false, BC = BC)
    end
end

max_stars_repo_name: sandreza/CLIMA

max_issues_repo_name: sandreza/CLIMA

max_forks_repo_name: sandreza/CLIMA


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -6.019960436071491 <= q <= 11.696817216382863)

                     add_NL_constraint(m, :((0.9042935565441375 + 0.31236426273908524*(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4]))*(1 + erf((0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4]))/sqrt(2)))/2 + -0.4128166127064836*(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4]))*(1 + erf((0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4]))/sqrt(2)))/2 + -0.2476169775927386*(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))*(1 + erf((0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))/sqrt(2)))/2)*(1 + erf((0.9042935565441375 + 0.31236426273908524*(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4]))*(1 + erf((0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4]))/sqrt(2)))/2 + -0.4128166127064836*(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4]))*(1 + erf((0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4]))/sqrt(2)))/2 + -0.2476169775927386*(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))*(1 + erf((0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + (0.8600467319105785 + -0.4006275013725942*(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4]))*(1 + erf((0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4]))/sqrt(2)))/2 + 0.22106887260345687*(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4]))*(1 + erf((0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4]))/sqrt(2)))/2 + 0.7210213948433983*(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))*(1 + erf((0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))/sqrt(2)))/2)*(1 + erf((0.8600467319105785 + -0.4006275013725942*(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4]))*(1 + erf((0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4]))/sqrt(2)))/2 + 0.22106887260345687*(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4]))*(1 + erf((0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4]))/sqrt(2)))/2 + 0.7210213948433983*(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))*(1 + erf((0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + (0.982756228863964 + 0.34792310275142935*(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4]))*(1 + erf((0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4]))/sqrt(2)))/2 + 0.7924515917531365*(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4]))*(1 + erf((0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4]))/sqrt(2)))/2 + -0.9399735455074207*(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))*(1 + erf((0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))/sqrt(2)))/2)*(1 + erf((0.982756228863964 + 0.34792310275142935*(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4]))*(1 + erf((0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4]))/sqrt(2)))/2 + 0.7924515917531365*(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4]))*(1 + erf((0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4]))/sqrt(2)))/2 + -0.9399735455074207*(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))*(1 + erf((0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

#
# remote node failure detection (RNFD) actor
# ----------------------------------------------
# if an actor on a (relative) remote node is added to a
# supervisor an actor is started scanning the remote link
# in regular intervals. If a ProcessExitedException is
# detected, it notifies the supervisor about it.
#
const scan_interval = 1.0

struct RNFD{L,S,T}
    sv::L      # supervisor link
    lks::S     # links to remote actors
    pids::T    # supervised pids
end
struct Add{L}
    lk::L
end
struct Remove{L}
    lk::L
end
struct Scan end

#
# RNFD behavior
#
function (rfd::RNFD)(msg::Add)
    if msg.lk ∉ rfd.lks
        push!(rfd.lks, msg.lk)
        msg.lk.pid ∈ rfd.pids || push!(rfd.pids, msg.lk.pid)
    end
end
function (rfd::RNFD)(msg::Remove)
    if msg.lk ∈ rfd.lks
        filter!(≠(msg.lk), rfd.lks)
        pids = (lk.pid for lk ∈ rfd.lks)
        filter!(∈(pids), rfd.pids)
    end
end
function (rfd::RNFD)(::Scan)
    isempty(rfd.lks) && return nothing
    excs = Int[]
    for lk in rfd.lks
        try
            isready(lk.chn)
        catch exc
            filter!(≠(lk), rfd.lks)
            exc isa ProcessExitedException && push!(excs, exc.worker_id)
        end
    end
    if !isempty(excs)
        filter!(lk->lk.pid ∉ excs, rfd.lks)
        filter!(p->p ∉ excs, rfd.pids)
        send(rfd.sv, NodeFailure(unique(excs)))
    end
end

#
# Supervisor behavior for NodeFailure
#

#
# remove temporary childs from both s and fchilds
#
function remove_temporary!(s, fchilds)
    act = task_local_storage("_ACT")
    filter!(fchilds) do child
        if child.info.restart == :temporary
            log_warn("temporary actor $(isnothing(child.name) ? :noname : child.name) failed, $(ProcessExitedException(child.lk.pid))")
            filter!(c->c.lk!=child.lk, act.conn)
            filter!(c->c.lk!=child.lk, s.childs)
            return false
        else
            return true
        end
    end
end
function restart_child!(c::Child, pid::Int)
    log_warn("supervisor: restarting child $(isnothing(c.name) ? :noname : c.name) on pid $pid")
    if c.lk isa Link
        lk = !isnothing(c.start) ? c.start(pid) : spawn(c.init; pid)
        c.lk.chn = lk.chn
        c.lk.pid = lk.pid
        isnothing(c.name) || update!(lk, c.name, s=:name)
    end
end
function restart!(s::Supervisor, cs::Vector{Child}, pids::Vector{Int})
    if s.option[:strategy] == :one_for_one
        for (i, c) in enumerate(cs)
            restart_child!(c, pids[i])
            rnfd_add(s, c.lk)
        end
    elseif s.option[:strategy] == :one_for_all
        log_warn("supervisor: restarting all")
        for child in s.childs
            child ∈ cs ?
                begin
                    restart_child!(child, pids[findfirst(==(child),cs)])
                    rnfd_add(s, child.lk)
                end :
                child.lk.mode ≠ :rnfd && shutdown_restart_child!(child)
        end
    else
        log_warn("supervisor: restarting rest")
        ix = findfirst(c->c ∈ cs, s.childs)
        for child in s.childs[ix:end]
            child ∈ cs ?
                begin
                    restart_child!(child, pids[findfirst(==(child),cs)])
                    rnfd_add(s, child.lk)
                end :
                child.lk.mode ≠ :rnfd && shutdown_restart_child!(child)
        end
    end
end
#
# return spare pids for failed childs cs and delete them from
# the s.option[:spares] dict entry
#
function spare_pids!(s::Supervisor, cs)
    spares = if haskey(s.option, :spares) && !isempty(s.option[:spares])
        s.option[:spares]
    else
        used = unique(map(c->c.lk.pid, s.childs))
        filter(p->p ∉ used, reverse(procs()))
    end
    filter!(p->p ∈ procs(), spares)
    pids = map(c->c.lk.pid, cs)
    p_old = sort(unique(pids))
    if length(p_old) ≤ length(spares)
        p_new = spares[1:length(p_old)]
        rp = [p_old[i]=>p_new[i] for i ∈ 1:length(p_old)]
        replace!(pids, rp...)
    elseif !isempty(spares)
        pids = rand(spares, length(pids))
    else
        pids = rand(procs(), length(pids))
    end
    haskey(s.option, :spares) && filter!(p->p ∉ pids, s.option[:spares])
    return pids
end
function (s::Supervisor)(msg::NodeFailure)
    foreach(msg.pids) do pid
        log_warn("supervisor: Process $pid exited!")
    end
    failed_childs = filter(c->c.lk.pid ∈ msg.pids, s.childs)
    remove_temporary!(s, failed_childs)
    if !isempty(failed_childs)
        if restart_limit!(s)
            log_warn("supervisor: restart limit $(s.option[:max_restarts]) exceeded!")
            send(self(), Exit(:shutdown, fill(nothing, 3)...))
        else
            restart!(s, failed_childs, spare_pids!(s, failed_childs))
        end
    end
end

#
# RNFD API
#
"""
    rnfd_start(sv::Link;; interval=1, kwargs...)

Start a RNFD actor and return a link to it.

# Arguments
- `sv::Link`: supervisor to the actor,
- `interval=1`: interval in seconds for checking remote nodes.
"""
function rnfd_start(sv::Link; interval=1, kwargs...)
    lk = spawn(RNFD(sv, Link[], Int[]); mode = :rnfd, kwargs...)
    exec(lk, supervise, sv)
    timer = Timer(interval; interval) do t
        send(lk, Scan())
    end
    term!(lk, (exp)->close(timer))
    return lk
end

#
# get a link to the RNFD actor, create it if it doesn't exist.
#
function rnfd(s::Supervisor)
    i = findfirst(c->c.lk.mode == :rnfd, s.childs)
    return !isnothing(i) ? s.childs[i].lk : rnfd_start(self())
end
rnfd_exists(s::Supervisor) = !isnothing(findfirst(c->c.lk.mode == :rnfd, s.childs))
#
# add a remote child to an RNFD actor,
# create it first if it doesn't exist.
#
rnfd_add(s::Supervisor, child::Link) = send(rnfd(s), Add(child))


text:
macro register(name)
    poly = esc(name)
    quote
        Base.convert(::Type{P}, p::P) where {P<:$poly} = p
        #Base.convert(P::Type{<:$poly}, p::$poly) where {T} = P(coeffs(p), p.var) # needs replacing
        Base.promote_rule(::Type{$poly{T}}, ::Type{$poly{S}}) where {T,S} =
            $poly{promote_type(T, S)}
        Base.promote_rule(::Type{$poly{T}}, ::Type{S}) where {T,S<:Number} =
            $poly{promote_type(T, S)}

        function (p::$poly)(x::AbstractVector)
            Base.depwarn(
                "Calling p(x::AbstractVector is deprecated. Use p.(x) instead.",
                Symbol("(p::AbstractPolynomial)"),
            )
            return p.(x)
        end

        #$poly(coeffs::AbstractVector{T}, var::SymbolLike = :x) where {T} =
        #    $poly{T}(coeffs, Symbol(var))
        #$poly(n::Number, var = :x) = $poly([n], var)
        #$poly{T}(n::S, var = :x) where {T,S<:Number} = $poly(T(n), var)
        #$poly{T}(x::AbstractVector{S}, var = :x) where {T,S<:Number} = $poly(T.(x), var)
    end
end


text:
# Load Julia packages (libraries) needed

using StatisticalRethinking, Optim
#gr(size=(600,600));

# ### snippet 2.6

# Grid of 1001 steps

p_grid = range(0, step=0.001, stop=1);

# all priors = 1.0

prior = ones(length(p_grid));

# Binomial pdf

likelihood = [pdf(Binomial(9, p), 6) for p in p_grid];

# As Uniform prior has been used, unstandardized posterior is equal to likelihood

posterior = likelihood .* prior;

# Scale posterior such that they become probabilities

posterior = posterior / sum(posterior);

# ### snippet 3.3

# Sample using the computed posterior values as weights

N = 10000
samples = sample(p_grid, Weights(posterior), N);

# In StatisticalRethinkingJulia samples will always be stored
# in an MCMCChains.Chains object.

chn = MCMCChains.Chains(reshape(samples, N, 1, 1), ["toss"]);

# Describe the chain

MCMCChains.describe(chn)

# Plot the chain

plot(chn)

# Compute the MAP (maximum_a_posteriori) estimate

x0 = [0.5]
lower = [0.0]
upper = [1.0]

function loglik(x)
  ll = 0.0
  ll += log.(pdf.(Beta(1, 1), x[1]))
  ll += sum(log.(pdf.(Binomial(9, x[1]), repeat([6], 1))))
  -ll
end

(qmap, opt) = quap(loglik, x0, lower, upper)

# Show optimization results

opt

# Fit quadratic approcimation

quapfit = [qmap[1], std(samples, mean=qmap[1])]

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 4)
p[1] = scatter(1:length(p_grid), samples, markersize = 2, ylim=(0.0, 1.3), lab="Draws")

# analytical calculation

w = 6
n = 9
x = 0:0.01:1
p[2] = plot( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")
density!(p[2], samples, lab="Sample density")

# quadratic approximation

p[3] = plot( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")
plot!( p[3], x, pdf.(Normal( quapfit[1], quapfit[2] ) , x ), lab="Quap approximation")

# ### snippet 2.7

w = 6; n = 9; x = 0:0.01:1
p[4] = plot( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")
plot!(p[4], x, pdf.(Normal( 0.67 , 0.16 ) , x ), lab="Normal approximation")
plot(p..., layout=(2, 2))

# End of `02/clip-06-07.jl`


text:
dispersion(obj::EconometricModel{<:LinearModelEstimators}) = true
dispersion(obj::EconometricModel{<:Union{NominalResponse,OrdinalResponse}}) = false
isiv(obj::EconometricModel{<:LinearModelEstimators}) = obj.iv > 0
coef(obj::EconometricModel) = obj.β
coefnames(obj::EconometricModel) = obj.vars
deviance(obj::EconometricModel) = -2 * loglikelihood(obj)
deviance(obj::EconometricModel{<:LinearModelEstimators}) =
    weights(obj) |> (
        wts -> isa(wts, FrequencyWeights) ?
            sum(w * (y - ŷ)^2 for (w, y, ŷ) in zip(wts, response(obj), fitted(obj))) :
            sum((y - ŷ)^2 for (y, ŷ) in zip(response(obj), fitted(obj)))
    )
"""
	hasintercept(obj::EconometricModel)::Bool

Return whether the model has an intercept.
"""
hasintercept(obj::EconometricModel) = true
hasintercept(obj::EconometricModel{<:LinearModelEstimators}) =
    !any(t -> isa(t, InterceptTerm{false}), terms(obj.f.rhs))
islinear(obj::EconometricModel{<:LinearModelEstimators}) =
    !(isa(obj.estimator, NominalResponse) || isa(obj.estimator, OrdinalResponse))
nulldeviance(obj::EconometricModel{<:LinearModelEstimators}) =
    meanresponse(obj) |> (ȳ -> sum((y - ȳ)^2 for y in response(obj)))
loglikelihood(obj::EconometricModel{<:LinearModelEstimators}) =
    √(deviance(obj) / dof_residual(obj)) |> (
        ϕ -> (sum(
            w * logpdf(Normal(μ, ϕ), y)
            for (w, μ, y) in zip(weights(obj), response(obj), fitted(obj))
        ))
    )
function loglikelihood(obj::EconometricModel{<:Union{NominalResponse,OrdinalResponse}})
    y = response(obj)
    categories = levels(y)
    b = mapreduce(elem -> (eachindex(categories) .== elem)', vcat, y)
    μ = predict(obj)
    wts = weights(obj)
    sum(
        wts[idx[1]] * logpdf(Categorical(collect(μ[idx[1], :])), idx[2])
        for idx in findall(b)
    )
end
function nullloglikelihood(obj::EconometricModel{<:LinearModelEstimators})
    ϕ = √(nulldeviance(obj) / (dof_residual(obj) + (dof(obj) - hasintercept(obj))))
    μ = meanresponse(obj)
    sum(w * logpdf(Normal(μ, ϕ), y) for (w, y) in zip(weights(obj), response(obj)))
end
function nullloglikelihood(obj::EconometricModel{<:Union{NominalResponse,OrdinalResponse}})
    y = response(obj)
    categories = levels(y)
    b = mapreduce(elem -> (eachindex(categories) .== elem)', vcat, y)
    μ = collect(vec(mean(b, dims = 1)))
    wts = weights(obj)
    sum(wᵢ * logpdf(Categorical(μ), yᵢ) for (yᵢ, wᵢ) in zip(y, wts))
end
# StatsBase.score(obj::StatisticalModel) = error("score is not defined for $(typeof(obj)).")
nobs(obj::EconometricModel) = sum(weights(obj))
dof(obj::EconometricModel{<:LinearModelEstimators}) =
    length(coef(obj)) + dispersion(obj) + obj.iv
dof(obj::EconometricModel{<:ContinuousResponse}) =
    length(coef(obj)) +
    dispersion(obj) +
    obj.iv +
    (obj.estimator.groups |> (g -> isempty(g) ? 0 : sum(length(group) - 1 for group in g)))
dof(obj::EconometricModel{<:Union{NominalResponse,OrdinalResponse}}) =
    length(coef(obj)) + dispersion(obj)
dof_residual(obj::EconometricModel{<:LinearModelEstimators}) =
    nobs(obj) - dof(obj) + dispersion(obj) + hasintercept(obj.f)
# Verify the dof_residual(obj::EconometricModel{<:Union{NominalResponse,OrdinalResponse}})
dof_residual(obj::EconometricModel{<:Union{NominalResponse,OrdinalResponse}}) =
    nobs(obj) - dof(obj) + dispersion(obj) + (length(obj.estimator.categories) - 1)
mss(obj::EconometricModel{<:LinearModelEstimators}) =
    meanresponse(obj) |>
    (ȳ -> sum(w * abs2(ŷ - ȳ) for (ŷ, w) in zip(fitted(obj), weights(obj))))
rss(obj::EconometricModel{<:LinearModelEstimators}) =
    residuals(obj) |> (û -> û' * Diagonal(weights(obj)) * û)
r2(obj::EconometricModel) = 1 - loglikelihood(obj) / nullloglikelihood(obj)
r2(obj::EconometricModel{<:LinearModelEstimators}) =
    (isiv(obj) || !hasintercept(obj)) ? NaN : rss(obj) |> (r -> 1 - r / (r + mss(obj)))
function adjr2(obj::EconometricModel{<:LinearModelEstimators})
    (isiv(obj) || !hasintercept(obj)) && return (NaN)
    ℓℓ = loglikelihood(obj)
    ℓℓ₀ = nullloglikelihood(obj)
    k = dof(obj)
    if variant == :McFadden
        1 - (ll - k) / ll0
    else
        error(":McFadden is the only currently supported variant")
    end
end
informationmatrix(obj::EconometricModel; expected::Bool = true) = obj.Ψ
"""
	vcov(obj::EconometricModel)
	vcov(obj::EconometricModel{<:LinearModelEstimators}, vce::VCE = obj.vce)

Return the variance-covariance matrix for the coefficients of the model.
The `vce` argument allows to request variance estimators.
"""
vcov(obj::EconometricModel) = informationmatrix(obj)
function vcov(obj::EconometricModel{<:LinearModelEstimators})
    vcov(obj::EconometricModel{<:LinearModelEstimators}, obj.vce)
end
function vcov(obj::EconometricModel{<:LinearModelEstimators}, vce::VCE)
    vce == OIM && return deviance(obj) / dof_residual(obj) * informationmatrix(obj)
    X = modelmatrix(obj)
    Ψ = informationmatrix(obj)
    û = residuals(obj)
    m = nobs(obj)
    k = dof(obj) - hasintercept(obj)
    p = dof_residual(obj)
    λ = vce == HC1 ? m / p : 1
    if vce == HC0 || vce == HC1
        ũ = û .^ 2
    else
        h = leverage(obj)
        if vce == HC2
            ũ = û .^ 2 ./ (1 .- h)
        elseif vce == HC3
            ũ = û .^ 2 ./ (1 .- h) .^ 2
        elseif vce == HC4
            ũ = û .^ 2 ./ (1 .- h) .^ min.(4, m / k * h)
        end
    end
    Ω = X' * Diagonal(ũ) * X
    Hermitian(λ * Ψ * Ω * Ψ)
end
"""
	stderror(obj::EconometricModel)
	stderror(obj::EconometricModel{<:LinearModelEstimators}, vce::VCE = obj.vce)

Return the standard errors for the coefficients of the model.
The `vce` argument allows to request variance estimators.
"""
stderror(obj::EconometricModel, vce::VCE) = sqrt.(diag(vcov(obj, vce)))
"""
	confint(obj::EconometricModel; se::AbstractVector{<:Real} = stderror(obj), level::Real = 0.95)

Compute the confidence intervals for coefficients, with confidence level `level` (by default, 95%).
`se` can be provided as a precomputed value.
"""
function confint(
    obj::EconometricModel;
    se::AbstractVector{<:Real} = stderror(obj),
    level::Real = 0.95,
)
    @assert zero(level) ≤ level ≤ one(level)
    α = 1 - level
    se * quantile(TDist(dof_residual(obj)), 1 - α / 2) |>
    (σ -> coef(obj) |> (β -> hcat(β .- σ, β .+ σ)))
end
weights(obj::EconometricModel) = obj.wts
isfitted(obj::EconometricModel) = !isempty(obj.Ψ)
fitted(obj::EconometricModel) = obj.ŷ
response(obj::EconometricModel) = obj.y
meanresponse(obj::EconometricModel{<:LinearModelEstimators}) =
    mean(response(obj), weights(obj))
modelmatrix(obj::EconometricModel) = obj.X
function leverage(obj::EconometricModel)
    X = modelmatrix(obj)
    Ψ = informationmatrix(obj)
    ω = weights(obj)
    diag(X * Ψ * X' * Diagonal(ω))
end
residuals(obj::EconometricModel{<:LinearModelEstimators}) = response(obj) - fitted(obj)
function residuals(obj::EconometricModel{<:NominalResponse})
    @unpack y = obj
    @unpack categories = obj.estimator
    # Change from levels to the categories after StatsModels patch
    l = levels(y)
    y = [findfirst(isequal(x), l) for x in y]
    b = mapreduce(elem -> (eachindex(categories) .== elem)', vcat, y)
    b - predict(obj)
end
predict(obj::EconometricModel) = fitted(obj)
predict(obj::EconometricModel{<:NominalResponse}) =
    mapslices(softmax, fitted(obj), dims = 2)
function predict(obj::EconometricModel{<:OrdinalResponse})
    y = response(obj)
    ŷ = fitted(obj)
    β = coef(obj)
    outcomes = obj.estimator.categories
    ζ = vcat(-Inf, β[end-length(outcomes)+2:end], Inf)
    mapreduce(
        row ->
            (
                cdf.(Logistic(), ζ[2:end] .- ŷ[row]) .-
                cdf.(Logistic(), ζ[1:end-1] .- ŷ[row])
            )',
        vcat,
        eachindex(y),
    )
end
"""
	coeftable(obj::EconometricModel;
		  level::Real = 0.95)
	coeftable(obj::EconometricModel{<:LinearModelEstimators};
		  level::Real = 0.95,
		  vce::VCE = obj.vce)

Return a table of class `CoefTable` with coefficients and related statistics.
`level` determines the level for confidence intervals (by default, 95%).
`vce` determines the variance-covariance estimator (by default, `OIM`).
"""
function coeftable(obj::EconometricModel; level::Real = 0.95, vce::VCE = obj.vce)
    β = coef(obj)
    σ = stderror(obj, vce)
    t = β ./ σ
    p = 2 * ccdf.(TDist(dof_residual(obj)), abs.(t))
    mat = hcat(β, σ, t, p, confint(obj, se = σ, level = level))
    lims = (100 * (1 - level) / 2, 100 * (1 - (1 - level) / 2))
    colnms = [
        "PE   ",
        "SE   ",
        "t-value",
        "Pr > |t|",
        string(@sprintf("%.2f", lims[1]), "%"),
        string(@sprintf("%.2f", lims[2]), "%"),
    ]
    rownms = coefnames(obj)[2]
    CoefTable(mat, colnms, rownms, 4)
end
function coeftable(obj::EconometricModel{<:NominalResponse}; level::Real = 0.95)
    @assert zero(level) ≤ level ≤ one(level)
    β = coef(obj)
    σ = stderror(obj)
    t = β ./ σ
    p = 2 * ccdf.(TDist(dof_residual(obj)), abs.(t))
    mat = hcat(β, σ, t, p, confint(obj, se = σ, level = level))
    lims = (100 * (1 - level) / 2, 100 * (1 - (1 - level) / 2))
    colnms = [
        "PE   ",
        "SE   ",
        "t-value",
        "Pr > |t|",
        string(@sprintf("%.2f", lims[1]), "%"),
        string(@sprintf("%.2f", lims[2]), "%"),
    ]
    vars = coefnames(obj)
    rownms = [string(lhs, " ~ ", rhs) for lhs in vars[1] for rhs in vars[2]]
    CoefTable(mat, colnms, rownms, 4)
end
function coeftable(obj::EconometricModel{<:OrdinalResponse}; level::Real = 0.95)
    @assert zero(level) ≤ level ≤ one(level)
    β = coef(obj)
    σ = stderror(obj)
    t = β ./ σ
    p = 2 * ccdf.(TDist(dof_residual(obj)), abs.(t))
    mat = hcat(β, σ, t, p, confint(obj, se = σ, level = level))
    lims = (100 * (1 - level) / 2, 100 * (1 - (1 - level) / 2))
    colnms = [
        "PE   ",
        "SE   ",
        "t-value",
        "Pr > |t|",
        string(@sprintf("%.2f", lims[1]), "%"),
        string(@sprintf("%.2f", lims[2]), "%"),
    ]
    vars = coefnames(obj)
    outcomes = obj.estimator.categories
    rownms = vcat(
        vars[2],
        [
            string("(Intercept): ", outcomes[i], " | ", outcomes[i+1])
            for i in 1:length(outcomes)-1
        ],
    )
    CoefTable(mat, colnms, rownms, 4)
end


text:
using Clapeyron, Test

@testset "SAFT models" begin
    T = 298.15
    V = 1e-4

    @testset "ogSAFT" begin
        system = ogSAFT(["water","ethylene glycol"])
        z = [0.5, 0.5]
        @test Clapeyron.a_seg(system, V, T, z) ≈ -2.0332062924093366 rtol = 1e-6
        @test Clapeyron.a_chain(system, V, T, z) ≈ -0.006317441684202759 rtol = 1e-6
        @test Clapeyron.a_assoc(system, V, T, z) ≈ -4.034042081699316 rtol = 1e-6
    end

    @testset "CKSAFT" begin
        system = CKSAFT(["carbon dioxide", "2-propanol"])
        z = [0.5, 0.5]
        @test Clapeyron.a_seg(system, V, T, z) ≈ -1.24586302917188 rtol = 1e-6
        @test Clapeyron.a_chain(system, V, T, z) ≈ -0.774758615408493 rtol = 1e-6
        @test Clapeyron.a_assoc(system, V, T, z) ≈ -1.2937079004096872 rtol = 1e-6
    end

    @testset "sCKSAFT" begin
        system = sCKSAFT(["benzene","acetic acid"])
        z = [0.5, 0.5]
        @test Clapeyron.a_seg(system, V, T, z) ≈ -3.1809330810925256 rtol = 1e-6
        @test Clapeyron.a_assoc(system, V, T, z) ≈ -3.3017434376105514 rtol = 1e-6
    end

    @testset "BACKSAFT" begin
        system = BACKSAFT(["carbon dioxide"])
        z = [1.]
        @test Clapeyron.a_hcb(system, V, T, z) ≈ 1.0118842111801198 rtol = 1e-6
        @test Clapeyron.a_chain(system, V, T, z) ≈ -0.14177009317268635 rtol = 1e-6
        @test Clapeyron.a_disp(system, V, T, z) ≈ -2.4492518566426296 rtol = 1e-6
    end

    @testset "LJSAFT" begin
        system = LJSAFT(["ethane","1-propanol"])
        z = [0.5, 0.5]
        @test Clapeyron.a_seg(system, V, T, z) ≈ -2.207632433058473 rtol = 1e-6
        @test Clapeyron.a_chain(system, V, T, z) ≈ -0.04577483379871112 rtol = 1e-6
        @test Clapeyron.a_assoc(system, V, T, z) ≈ -1.3009761155167205 rtol = 1e-6
    end

    @testset "CPA" begin
        system = CPA(["ethanol","benzene"])
        z = [0.5, 0.5]
        @test Clapeyron.a_SRK(system, V, T, z) ≈ 4.510022402195623 rtol = 1e-6
        @test Clapeyron.a_assoc(system, V, T, z) ≈ -1.1575210505284332 rtol = 1e-6
    end

    @testset "SAFTVRSW" begin
        system = SAFTVRSW(["water", "ethane"])
        z = [0.5, 0.5]
        @test Clapeyron.a_mono(system, V, T, z) ≈ -1.4367205951569462 rtol = 1e-6
        @test Clapeyron.a_chain(system, V, T, z) ≈ 0.024000058201261557 rtol = 1e-6
        @test Clapeyron.a_assoc(system, V, T, z) ≈ -0.5238154638538838 rtol = 1e-6
    end

    @testset "softSAFT" begin
        system = softSAFT(["hexane","1-propanol"])
        z = [0.5,0.5]
        @test Clapeyron.a_LJ(system, V, T, z) ≈ -3.960728242264164 rtol = 1e-6
        @test Clapeyron.a_chain(system, V, T, z) ≈ 0.3736728407455211 rtol = 1e-6
        @test Clapeyron.a_assoc(system, V, T, z) ≈ -2.0461376618069034 rtol = 1e-6
    end

    @testset "PCSAFT" begin
        system = PCSAFT(["butane", "ethanol"])
        z = [0.5, 0.5]
        @test Clapeyron.a_hc(system, V, T, z) ≈ 3.1148229872928654 rtol = 1e-6
        @test Clapeyron.a_disp(system, V, T, z) ≈ -6.090736508783152 rtol = 1e-6
        @test Clapeyron.a_assoc(system, V, T, z) ≈ -1.6216064387201956 rtol = 1e-6
    end

    @testset "sPCSAFT" begin
        system = sPCSAFT(["propane", "methanol"])
        z = [0.5, 0.5]
        @test Clapeyron.a_hc(system, V, T, z) ≈ 2.024250583187793 rtol = 1e-6
        @test Clapeyron.a_disp(system, V, T, z) ≈ -4.138653131750594 rtol = 1e-6
        @test Clapeyron.a_assoc(system, V, T, z) ≈ -1.1459701721909195 rtol = 1e-6
        #difference in this error is almost exactly 0.5, suspicious
    end

    @testset "SAFTVRMie" begin
        system = SAFTVRMie(["methanol", "water"])
        z = [0.5, 0.5]
        @test Clapeyron.a_mono(system, V, T, z) ≈ -0.9729134860869052 rtol = 1e-6
        @test Clapeyron.a_chain(system, V, T, z) ≈ -0.02834738013535014 rtol = 1e-6
        @test Clapeyron.a_assoc(system, V, T, z) ≈ -4.180807072390184 rtol = 1e-6
    end

    @testset "SAFTVRQMie" begin
        system = SAFTVRQMie(["helium"])
        z = [1.]
        @test Clapeyron.a_mono(system, V, T, z) ≈ 0.12286776703976324 rtol = 1e-6
    end

    @testset "SAFTgammaMie" begin
        system = SAFTgammaMie(["methanol","butane"])
        V_γMie = exp10(-3.5)
        z = [0.5,0.5]
        @test Clapeyron.a_mono(system, V_γMie, T, z) ≈ -1.0400249396482548 rtol = 1e-6
        @test Clapeyron.a_chain(system, V_γMie, T, z) ≈ -0.07550931466871749 rtol = 1e-6
        @test Clapeyron.a_assoc(system, V_γMie, T, z) ≈ -0.8205840455850311 rtol = 1e-6
    end
end

@testset "Cubic models" begin
    T = 333.15
    V = 1e-3
    p = 1e5
    z = [0.5, 0.5]

    @testset "vdW Models" begin
            system = vdW(["ethane","undecane"])
            @test Clapeyron.a_res(system, V, T, z) ≈ -0.7088380780265725 rtol = 1e-6
            @test Clapeyron.cubic_poly(system, p, T, z)[1][1] ≈ -0.0002475728429728521 rtol = 1e-6
            @test Clapeyron.cubic_abp(system, V, T, z)[1] ≈ 2.591962221523441 rtol = 1e-6
    end

    @testset "RK Models" begin
        @testset "Default RK" begin
            system = RK(["ethane","undecane"])
            @test Clapeyron.a_res(system, V, T, z) ≈ -0.9825375012134132 rtol = 1e-6
            @test Clapeyron.cubic_poly(system, p, T, z)[1][1] ≈ -0.00022230043592123767 rtol = 1e-6
            @test Clapeyron.cubic_abp(system, V, T, z)[1] ≈ 3.357807485319651 rtol = 1e-6
        end

        @testset "SRK" begin
            system = SRK(["ethane","undecane"])
            @test Clapeyron.a_res(system, V, T, z) ≈ -1.2572506872856557 rtol = 1e-6
        end

        @testset "PSRK" begin
            system = PSRK(["ethane","undecane"])
            @test Clapeyron.a_res(system, V, T, z) ≈ -1.1994085074805483 rtol = 1e-6
        end

        @testset "RK w/ BMAlpha" begin
            system = RK(["ethane","undecane"];alpha = BMAlpha)
            @test Clapeyron.a_res(system, V, T, z) ≈ -1.2569334957019538 rtol = 1e-6
        end

        @testset "RK w/ PenelouxTranslation" begin
            system = RK(["ethane","undecane"];translation = PenelouxTranslation)
            @test Clapeyron.a_res(system, V, T, z) ≈ -0.9800472116681871 rtol = 1e-6
        end

        @testset "RK w/ KayRule" begin
            system = RK(["ethane","undecane"];mixing = KayRule)
            @test Clapeyron.a_res(system, V, T, z) ≈ -0.8176850121211936 rtol = 1e-6
        end

        @testset "RK w/ HVRule" begin
            system = RK(["methanol","benzene"];mixing = HVRule, activity=Wilson)
            @test Clapeyron.a_res(system, V, T, z) ≈ -0.5209112693371991 rtol = 1e-6
        end

        @testset "RK w/ MHV1Rule" begin
            system = RK(["methanol","benzene"];mixing = MHV1Rule, activity=Wilson)
            @test Clapeyron.a_res(system, V, T, z) ≈ -0.5091065987876959 rtol = 1e-6
        end

        @testset "RK w/ MHV2Rule" begin
            system = RK(["methanol","benzene"];mixing = MHV2Rule, activity=Wilson)
            @test Clapeyron.a_res(system, V, T, z) ≈ -0.5071048453490448 rtol = 1e-6
        end

        @testset "RK w/ WSRule" begin
            system = RK(["methanol","benzene"];mixing = WSRule, activity=Wilson)
            @test Clapeyron.a_res(system, V, T, z) ≈ -0.572912690389026 rtol = 1e-6
        end
    end

    @testset "PR Models" begin
        @testset "Default PR" begin
            system = PR(["ethane","undecane"])
            @test Clapeyron.a_res(system, V, T, z) ≈ -1.244772730766631 rtol = 1e-6
            @test Clapeyron.cubic_poly(system, p, T, z)[1][1] ≈ -0.00023285390449318037 rtol = 1e-6
            @test Clapeyron.cubic_abp(system, V, T, z)[1] ≈ 4.268630968024985 rtol = 1e-6
        end

        @testset "PR78" begin
            system = PR78(["ethane","undecane"])
            @test Clapeyron.a_res(system, V, T, z) ≈ -1.246269686941749 rtol = 1e-6
        end

        @testset "VTPR" begin
            system = VTPR(["ethane","undecane"])
            @test Clapeyron.a_res(system, V, T, z) ≈ -1.23401133834732 rtol = 1e-6
        end

        @testset "UMRPR" begin
            system = UMRPR(["ethane","undecane"])
            @test Clapeyron.a_res(system, V, T, z) ≈ -1.1447318247939071 rtol = 1e-6
        end

        @testset "PR w/ BMAlpha" begin
            system = PR(["ethane","undecane"];alpha = BMAlpha)
            @test Clapeyron.a_res(system, V, T, z) ≈ -1.244507550417118 rtol = 1e-6
        end

        @testset "PR w/ TwuAlpha" begin
            system = PR(["ethane","undecane"];alpha = TwuAlpha)
            @test Clapeyron.a_res(system, V, T, z) ≈ -1.2650743158660063 rtol = 1e-6
        end

        @testset "PR w/ MTAlpha" begin
            system = PR(["ethane","undecane"];alpha = MTAlpha)
            @test Clapeyron.a_res(system, V, T, z) ≈ -1.2542333213442207 rtol = 1e-6
        end

        @testset "PR w/ RackettTranslation" begin
            system = PR(["ethane","undecane"];translation = RackettTranslation)
            @test Clapeyron.a_res(system, V, T, z) ≈ -1.2453840474643165 rtol = 1e-6
        end

        @testset "PR w/ MTTranslation" begin
            system = PR(["ethane","undecane"];translation = MTTranslation)
            @test Clapeyron.a_res(system, V, T, z) ≈ -1.24391715980424 rtol = 1e-6
        end

        @testset "PR w/ HVRule" begin
            system = PR(["methanol","benzene"];mixing = HVRule, activity=Wilson)
            @test Clapeyron.a_res(system, V, T, z) ≈ -0.632982061564318 rtol = 1e-6
        end

        @testset "PR w/ MHV1Rule" begin
            system = PR(["methanol","benzene"];mixing = MHV1Rule, activity=Wilson)
            @test Clapeyron.a_res(system, V, T, z) ≈ -0.6211434867194446 rtol = 1e-6
        end

        @testset "PR w/ MHV2Rule" begin
            system = PR(["methanol","benzene"];mixing = MHV2Rule, activity=Wilson)
            @test Clapeyron.a_res(system, V, T, z) ≈ -0.6210836570941939 rtol = 1e-6
        end

        @testset "PR w/ LCVMRule" begin
            system = PR(["methanol","benzene"];mixing = LCVMRule, activity=Wilson)
            @test Clapeyron.a_res(system, V, T, z) ≈ -0.6286234264772204 rtol = 1e-6
        end

        @testset "PR w/ WSRule" begin
            system = PR(["methanol","benzene"];mixing = WSRule, activity=Wilson)
            @test Clapeyron.a_res(system, V, T, z) ≈ -0.669085674824878 rtol = 1e-6
        end
    end
end

@testset "Activity models" begin
    T = 333.15
    p = 1e5
    z = [0.5,0.5]

    @testset "Wilson" begin
        system = Wilson(["methanol","benzene"])
        @test Clapeyron.activity_coefficient(system,p,T,z)[1] ≈ 1.530046633499114 rtol = 1e-6
    end

    @testset "NRTL" begin
        system = NRTL(["methanol","benzene"])
        @test Clapeyron.activity_coefficient(system,p,T,z)[1] ≈ 1.5309354738922405 rtol = 1e-6
    end

    @testset "UNIQUAC" begin
        system = UNIQUAC(["methanol","benzene"])
        @test Clapeyron.activity_coefficient(system,p,T,z)[1] ≈ 1.3630421218486388 rtol = 1e-6
    end

    @testset "UNIFAC" begin
        system = UNIFAC(["methanol","benzene"])
        @test Clapeyron.activity_coefficient(system,p,T,z)[1] ≈ 1.5322232657797463 rtol = 1e-6
    end

    @testset "COSMOSAC02" begin
        system = COSMOSAC02(["water","ethanol"])
        @test Clapeyron.activity_coefficient(system,p,T,z)[1] ≈ 1.3871817962565904 rtol = 1e-6
    end

    @testset "COSMOSAC10" begin
        system = COSMOSAC10(["water","ethanol"])
        @test Clapeyron.activity_coefficient(system,p,T,z)[1] ≈ 1.4015660588643404 rtol = 1e-6
    end

    @testset "COSMOSACdsp" begin
        system = COSMOSACdsp(["water","ethanol"])
        @test Clapeyron.activity_coefficient(system,p,T,z)[1] ≈ 1.4398951117248127 rtol = 1e-6
    end
end

@testset "Ideal models" begin
    T = 298.15
    V = 1e-4
    z = [1.]

    @testset "Joback" begin
        system = JobackIdeal(["hexane"])
        @test Clapeyron.C_p(system,298.15) ≈ 143.22076150138616 rtol = 1e-6
        @test Clapeyron.T_b(system) ≈ 336.88 rtol = 1e-6
        @test Clapeyron.crit_pure(system)[1] ≈ 500.2728274871347 rtol = 1e-6
        @test Clapeyron.a_ideal(system,V,T,z) ≈ 9.210841420941021 rtol = 1e-6
    end

    @testset "Reid" begin
        system = ReidIdeal(["butane"])
        @test Clapeyron.a_ideal(system,V,T,z) ≈ 9.210842104089576 rtol = 1e-6
    end

    @testset "Walker" begin
        system = WalkerIdeal(["hexane"])
        @test Clapeyron.a_ideal(system,V,T,z) ≈ 179.51502015696653 rtol = 1e-6
    end

    @testset "Monomer" begin
        system = MonomerIdeal(["hexane"])
        @test Clapeyron.a_ideal(system,V,T,z) ≈ -20.368750666236373 rtol = 1e-6
    end
end

@testset "Multi-parameter models" begin
    T = 298.15
    V = 1e-4

    @testset "IAPWS95" begin
        z = [1.]
        system = IAPWS95()
        @test Clapeyron.a_ideal(system, V, T, z) ≈ 7.932118505056652 rtol = 1e-6
        @test Clapeyron.a_res(system, V, T, z) ≈ -2.1152657050144347e14 rtol = 1e-6
    end

    @testset "PropaneRef" begin
        z   = [1.]
        system = PropaneRef()
        @test Clapeyron.a_ideal(system, V, T, z) ≈ 0.6426994942361217 rtol = 1e-6
        @test Clapeyron.a_res(system, V, T, z) ≈ -2.436280448227229 rtol = 1e-6
    end

    @testset "GERG2008" begin
        z   = [1.]
        system = GERG2008(["water"])
        @test Clapeyron.a_ideal(system, V, T, z) ≈ 4.500099760879548 rtol = 1e-6
        @test Clapeyron.a_res(system, V, T, z) ≈ -10.122119572808764 rtol = 1e-6
        z   = [0.25,0.25,0.25,0.25]
        system = GERG2008(["water","carbon dioxide","hydrogen sulfide","argon"])
        @test Clapeyron.a_ideal(system, V, T, z) ≈ 3.1135835641766594 rtol = 1e-6
        @test Clapeyron.a_res(system, V, T, z) ≈ -1.1706377677539772 rtol = 1e-6
    end
end

@testset "SPUNG models" begin
    T = 298.15
    V = 1e-4
    z = [1.]

    @testset "SRK" begin
        system = SPUNG(["ethane"])
        @test Clapeyron.shape_factors(system, V, T, z)[1] ≈ 0.7289071312821193 rtol = 1e-6
    end

    @testset "PCSAFT" begin
        system = SPUNG(["ethane"],PropaneRef(),PCSAFT(["ethane"]),PCSAFT(["propane"]))
        @test Clapeyron.shape_factors(system, V, T, z)[1] ≈ 1.3499576779924594 rtol = 1e-6
    end
end