text:
#This file defines methods for using the EM algorithm to converge on a joint inference of basic parameters and the daily contact rate

#Define the initial 'guess' for the contact rate based on Google mobility data
cts₀ = [KenyaCoVSD.ct_kenya(t,0.55) for (t,d) in enumerate(Date(2020,2,21):Day(1):Date(2021,4,15))]

#Define variable transformation

trans_one_group = as((R₀ = as(Real, 0.0, 7.0),
        χ=as(Real, 0.1, 30.0),
        p_test=as(Real, 0.0, 50.0),
        extra_transmissibility = as(Real, 0.5, 2.5),
        influx_exposed_new_variant = as(Real, 0.0, 10e3),
        p_test_boost = as(Real, 1.0, 3.0),
        χ_boost= as(Real, 1.0, 3.0),
        E₀=as(Real, 0.0, 10e3)))
D = Diagonal(0.1*ones(TransformVariables.dimension(trans_one_group)))

#Modified MCMC methods for working with daily fitted ct value

function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal,cts,q₀;serowaningrate = 1/365,num_chains = 1)
        println("Starting MCMC parameter inference for $(model.areaname)")
        n = length(q₀)
        l_area(x) = model.log_likelihood(x,model,serowaningrate,cts) + model.log_priors(x)
        ℓ = TransformedLogDensity(trans, l_area)#transformed log-likelihood
        ∇ℓ = LogDensityProblems.ADgradient(:ForwardDiff, ℓ)#transformed log-likelihood gradient wrt the parameters

        results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇ℓ, samples,
                                initialization = (q = q₀,κ=GaussianKineticEnergy(D),ϵ = stepsize),
                                warmup_stages = fixed_stepsize_warmup_stages(M=Symmetric))
        # results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇ℓ, samples,
        #                         initialization = (q = q₀,κ=GaussianKineticEnergy(D)))

        transformed_results = TransformVariables.transform.(trans,results.chain)
        val = zeros(length(transformed_results),length(transformed_results[1]),1)

        for i = 1:size(val,1),j = 1:size(val,2)
                val[i,j,1] = transformed_results[i][j]
        end

        for chain_reps = 2:num_chains
                println("Running additional chain $(chain_reps) for $(model.areaname)")
                results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇ℓ, samples,
                                        initialization = (q = q₀,κ=GaussianKineticEnergy(D),ϵ = stepsize),
                                        warmup_stages = fixed_stepsize_warmup_stages(M=Symmetric,stepsize_search = nothing),
                                        reporter = NoProgressReport())
                transformed_results = TransformVariables.transform.(trans,results.chain)
                for i = 1:size(val,1),j = 1:size(val,2)
                        val[i,j,chain_reps] = transformed_results[i][j]
                end
        end

        chn = Chains(val,[String(k) for k in keys(transformed_results[1])])

        model.MCMC_results = KenyaCoVSD.MCMCResults(chn,
                                        [l_area(transformed_results[i]) - model.log_priors(transformed_results[i]) for i = 1:length(transformed_results)],
                                        results.tree_statistics)

        return nothing
end

function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal,cts;serowaningrate = 1/365,num_chains = 1)
        println("Searching for a good initial condition for $(model.areaname)")
        l_area(x) = model.log_likelihood(x,model,serowaningrate,cts) + model.log_priors(x)
        f(x) = -transform_logdensity(trans, l_area,x)
        searchrange = fill((-3.,3.),TransformVariables.dimension(trans))
        res = bboptimize(f; SearchRange = searchrange,PopulationSize=1000,MaxSteps=30000)
        q₀ = best_candidate(res)
        inferparameters!(model,samples,trans,stepsize,D,cts,q₀;serowaningrate=serowaningrate,num_chains=num_chains)
        return nothing
end


text:
# run_logit_model_1.jl

# USER MUST SET THESE:
collatedPath = "/Volumes/Neurobio/Assad Lab/allison/Collated Analyses/JULIA_CSVs/All_SNc_nosplitsesh"
compositesavepath="" # the save directory will be formed automagically!
runID = 1 # set runID nonZero to specify which run it is.

# Use the latest versions of our files and dependencies!
cd("/Users/lilis/Documents/GitHub/HSOManalysisPackages/JuliaLanguage_Hierarchical_Models")
include("/Users/lilis/Documents/GitHub/HSOManalysisPackages/JuliaLanguage_Hierarchical_Models/function_library/file_modules.jl")
refresh_tools("/Users/lilis/Documents/GitHub/HSOManalysisPackages/JuliaLanguage_Hierarchical_Models/function_library", exact=true)
println("	")
# CANNED FOR LOGIT:
modelpackagefunction = bootlogit_modelpackage1
postprocessingfunction = bootlogit_postprocessingfunction1

warning("suppressing figures from individual sessions, but they will still be saved to host folders for that session.")
println("	")
(fails, results) = run_collated_model(collatedPath, modelpackagefunction; pathIDx = [],runFails=false, failDirs=[], postprocessingfunction=postprocessingfunction,
    compositesavepath=compositesavepath, runID=runID, suppressFigures=true);

ext: jl


text:
import RegressionTables: columns

align_table(tab) = "l" * ("r" ^ (columns(tab) - 1))

function column_widths(tab, align)
    n_columns = columns(tab)

    if length(align) != n_columns
        error("align string has invalid length.")
    end

    # construct column width, first from the maximum of the bodies' column widths.
    colWidths = zeros(Int64, n_columns)
    for colIndex = 1:n_columns
        if (align[colIndex] == 'l') || (align[colIndex] == 'r') || (align[colIndex] == 'c')
            colWidths[colIndex] = maximum([length(b[r,colIndex]) for b in tab.bodies for r=1:size(b,1) ])
        else
            error("Invalid character in align string. Only 'l', 'r', 'c' are allowed.")
        end
    end

    colWidths
end

function header(tab, header, settings, colWidths)
  # header
  headerLabels = Vector{String}(undef,0)
  headerWidths = Vector{Int64}(undef,0)
  headerCellStartEnd = Vector{Vector{Int64}}(undef,0)
  # first column is empty (top left)
  push!(headerLabels, "")
  push!(headerWidths, colWidths[1])
  push!(headerCellStartEnd, [1,1])
  # first regression result
  push!(headerLabels, header[1,2])
  push!(headerWidths, colWidths[2])
  push!(headerCellStartEnd, [2,2])
  if columns(tab)>2
      for rIndex = 3:size(header,2)
          if header[1,rIndex] == header[1,rIndex-1]
              headerWidths[end] += length(settings.headercolsep) + colWidths[rIndex]
              headerCellStartEnd[end][2] += 1
          else
              push!(headerLabels, header[1,rIndex])
              push!(headerWidths, colWidths[rIndex])
              push!(headerCellStartEnd, [rIndex, rIndex])
          end
      end
  end

  headerArray = Array{String}(undef,1,length(headerLabels))
  headerArray[1,1] = ""
  for i = 2:size(headerArray,2)
      headerArray[1,i] = settings.encapsulateRegressand(headerLabels[i],headerCellStartEnd[i][1],headerCellStartEnd[i][2] )
  end

  (headerArray=headerArray, headerLabels=headerLabels, headerWidths=headerWidths, headerCellStartEnd=headerCellStartEnd)
end

"""
    now it could be that a columns of the header is wider than the colWidth of the bodies.
    in that case, make the last column of the bodies wider.
"""
function adjust_widths!(colWidths, header, settings)
    @unpack headerArray, headerCellStartEnd, headerLabels, headerWidths = header

    adjusted = false

    for i = 2:size(headerArray,2)
        totalWidth = sum([colWidths[cind] for cind in headerCellStartEnd[i][1] : headerCellStartEnd[i][2]]) + length(settings.colsep)*(headerCellStartEnd[i][2] - headerCellStartEnd[i][1])
        if length(headerArray[1,i])>totalWidth
            adjusted = true
            # extend width of cells
            colWidths[headerCellStartEnd[i][2]] +=  (length(headerArray[1,i])-totalWidth)
            headerWidths[i] += (length(headerArray[1,i])-totalWidth)
        end
    end
    (adjusted = adjusted,
     h = (headerArray=headerArray, headerLabels=headerLabels, headerWidths=headerWidths, headerCellStartEnd=headerCellStartEnd))
end

# distinguish two cases:
#   if headerrule gives a string of length one, put into table, and repeat the string
function headerrule(header, settings)
    @unpack headerArray, headerCellStartEnd, headerLabels, headerWidths = header
    headerArray = copy(headerArray)

    hr = settings.headerrule(headerCellStartEnd)
    if length(hr) == 1
        secondRow = Array{String}(undef,1,length(headerLabels))
        secondRow[1,1] = ""
        for i = 2:size(secondRow,2)
            secondRow[1,i] = (length(settings.midrule) == 1 ? settings.midrule ^ headerWidths[i] : settings.midrule)
        end
        headerArray = [headerArray; secondRow]
        print_headerrule_separately = false
    else
        print_headerrule_separately = true
    end
    (headerArray = headerArray,
     hr = (headerrule=hr, print_headerrule_separately=print_headerrule_separately))
end


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

"""
# Measure (abbr: msr)
"""
# sets
const MsrSet   = Set{Symbol}([:flow,:power])
const UnitDict = Dict{Symbol,_UF.FreeUnits}(:flow   => u"m^3/hr",
                                            :power  => u"MW")
const MaxVal   = Dict{Symbol,_UF.Quantity}( :flow   => (Inf)u"m^3/hr",
                                            :power  => (Inf)u"MW")
# functions
get_max(msr::Symbol) = MaxVal[msr]
get_unit(msr::Symbol) = UnitDict[msr]
check_unit(msr::Symbol,unit::_UF.FreeUnits) =
    _UF.dimension(UnitDict[msr]) == _UF.dimension(unit)
has_msr(kwargs::Iterators.Pairs) =
    !isempty(collect(intersect(MsrSet,keys(kwargs))))
get_msr(ugf::AbstractUGF) = [ugf.msr]
get_msr(kwargs::Iterators.Pairs) = collect(intersect(MsrSet,keys(kwargs)))[1]
get_msr(std::AbstractSTD) = haskey(std.props,:msr) ? [std.props[:msr]] : [] ;
get_msr(ntw::AbstractNetwork) = haskey(ntw.props,:msr) ? [ntw.props[:msr]] : [] ;
function set_msr!(ntw::AbstractNetwork)
    msr = []
    for ne in elements(ntw)
        if haskey(ne,:ugf) union!(msr,get_msr(ne[:ugf])) end
        if haskey(ne,:std) union!(msr,get_msr(ne[:std])) end
        if haskey(ne,:ntw) union!(msr,get_msr(ne[:ntw][1])) end
    end
    ntw.props[:msr] = msr[1]                                                    # currently, only one msr is supported.
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Base.Markdown
import Base.Markdown: MD, Paragraph, Header, Italic, Bold, LineBreak, plain, term, html, rst, Table, Code, LaTeX, Footnote
import Base: show

# Basics
# Equality is checked by making sure the HTML output is
# the same – the structure itself may be different.

@test md"foo" == MD(Paragraph("foo"))
@test md"foo *bar* baz" == MD(Paragraph(["foo ", Italic("bar"), " baz"]))
@test md"""foo
bar""" == MD(Paragraph(["foo bar"]))
@test md"""foo\
bar""" == MD(Paragraph(["foo", LineBreak(), "bar"]))

@test md"#no title" == MD(Paragraph(["#no title"]))
@test md"# title" == MD(Header{1}("title"))
@test md"""
  #
  empty
  """ == MD(Header{1}(""), Paragraph("empty"))
@test md"## section" == MD(Header{2}("section"))
@test md"# title *foo* `bar` **baz**" ==
    MD(Header{1}(["title ", Italic("foo")," ",Code("bar")," ",Bold("baz")]))
@test md"""
h1
===""" == md"# h1"
@test md"""
h2
   ---""" == md"## h2"

@test md"**foo *bar* baz**" == MD(Paragraph(Bold(["foo ", Italic("bar"), " baz"])))
@test md"*foo **bar** baz*" == MD(Paragraph(Italic(["foo ", Bold("bar"), " baz"])))

@test md"""```julia
foo
```
""" == MD(Code("julia", "foo"))

@test md"foo ``bar`` baz" == MD(Paragraph(["foo ", LaTeX("bar"), " baz"]))

@test md"""
```math
...
```
""" == MD(LaTeX("..."))

code_in_code = md"""
````
```
````
"""
@test code_in_code == MD(Code("```"))
@test plain(code_in_code) == "````\n```\n````\n"

let text = "Foo ```bar` ``baz`` ```\n",
    md = Markdown.parse(text)
    @test text == Markdown.plain(md)
end

@test isempty(Markdown.parse("\r"))
@test Markdown.parse("hello\r") == MD(Paragraph(["hello"]))
@test Markdown.parse("hello\r*julia*") == MD(Paragraph(Any["hello ", Italic(Any["julia"])]))

@test md"A footnote [^foo]." == MD(Paragraph(["A footnote ", Footnote("foo", nothing), "."]))

@test md"[^foo]: footnote" == MD([Footnote("foo", Any[Paragraph(Any["footnote"])])])

let text =
    """
    A paragraph with some footnotes,[^1] and another.[^note]

    [^1]: Footnote text for the first.

    [^note]: A longer footnote:

        Indented paragraphs are part of the footnote.

            some.code

        And *another* paragraph.

    This isn't part of the footnote.
    """,
    md = Markdown.parse(text)
    @test length(md.content) == 4
    @test isa(md.content[1], Markdown.Paragraph)
    @test isa(md.content[2], Markdown.Footnote)
    @test isa(md.content[3], Markdown.Footnote)
    @test isa(md.content[4], Markdown.Paragraph)

    @test md.content[2].id == "1"
    @test md.content[3].id == "note"

    @test length(md.content[3].text) == 4

    let expected =
            """
            A paragraph with some footnotes,[^1] and another.[^note]

            [^1]: Footnote text for the first.

            [^note]:
                A longer footnote:

                Indented paragraphs are part of the footnote.

                ```
                some.code
                ```

                And *another* paragraph.

            This isn't part of the footnote.
            """
        @test Markdown.plain(md) == expected
    end
    let expected =
            """
            A paragraph with some footnotes,[1]_ and another.[note]_

            .. [1] Footnote text for the first.

            .. [note]
               A longer footnote:

               Indented paragraphs are part of the footnote.

               .. code-block:: julia

                   some.code

               And *another* paragraph.

            This isn't part of the footnote.
            """
        @test Markdown.rst(md) == expected
    end
    let html = Markdown.html(md)
        @test contains(html, ",<a href=\"#footnote-1\" class=\"footnote\">[1]</a>")
        @test contains(html, ".<a href=\"#footnote-note\" class=\"footnote\">[note]</a>")
        @test contains(html, "<div class=\"footnote\" id=\"footnote-1\"><p class=\"footnote-title\">1</p>")
        @test contains(html, "<div class=\"footnote\" id=\"footnote-note\"><p class=\"footnote-title\">note</p>")
    end
    let latex = Markdown.latex(md)
        @test contains(latex, ",\\footnotemark[1]")
        @test contains(latex, ".\\footnotemark[note]")
        @test contains(latex, "\n\\footnotetext[1]{Footnote text for")
        @test contains(latex, "\n\\footnotetext[note]{A longer footnote:\n")
    end
end

let doc = md"""
* one
* two

1. pirate
2. ninja
3. zombie"""
    @test isa(doc.content[1], Markdown.List)
    @test isa(doc.content[2], Markdown.List)
    @test doc.content[1].items[1][1].content[1] == "one"
    @test doc.content[1].items[2][1].content[1] == "two"
    @test doc.content[2].items[1][1].content[1] == "pirate"
    @test doc.content[2].items[2][1].content[1] == "ninja"
    @test doc.content[2].items[3][1].content[1] == "zombie"
end

let doc = Markdown.parse(
        """
        A paragraph...
        - one
        - two
           * three
           * four
        ... another paragraph.
        """
    )
    @test length(doc.content) === 3
    @test isa(doc.content[1], Markdown.Paragraph)
    @test isa(doc.content[2], Markdown.List)
    @test isa(doc.content[3], Markdown.Paragraph)

    @test length(doc.content[2].items) === 2
    @test doc.content[2].items[1][1].content[1] == "one"
    @test length(doc.content[2].items[2]) == 2
    @test doc.content[2].items[2][1].content[1] == "two"

    @test isa(doc.content[2].items[2][2], Markdown.List)
    @test length(doc.content[2].items[2][2].items) === 2
    @test doc.content[2].items[2][2].items[1][1].content[1] == "three"
    @test doc.content[2].items[2][2].items[2][1].content[1] == "four"
end

@test md"Foo [bar]" == MD(Paragraph("Foo [bar]"))
@test md"Foo [bar](baz)" != MD(Paragraph("Foo [bar](baz)"))
@test md"Foo \[bar](baz)" == MD(Paragraph("Foo [bar](baz)"))

# Basic plain (markdown) output
@test md"foo" |> plain == "foo\n"
@test md"foo *bar* baz" |> plain == "foo *bar* baz\n"
@test md"# title" |> plain == "# title\n"
@test md"## section" |> plain == "## section\n"
@test md"## section `foo`" |> plain == "## section `foo`\n"
@test md"""Hello

---
World""" |> plain == "Hello\n\n---\n\nWorld\n"
@test md"[*a*](b)" |> plain == "[*a*](b)\n"
@test md"""
> foo
>
>   * bar
>
> ```
> baz
> ```""" |> plain == """> foo\n>\n>   * bar\n>\n> ```\n> baz\n> ```\n\n"""

# Terminal (markdown) output

# multiple whitespace is ignored
@test sprint(term, md"a  b") == "  a b\n"
@test sprint(term, md"[x](https://julialang.org)") == "  x (https://julialang.org)\n"
@test sprint(term, md"[x](@ref)") == "  x\n"
@test sprint(term, md"[x](@ref something)") == "  x\n"
@test sprint(term, md"![x](https://julialang.org)") == "  (Image: x)\n"

# enumeration is normalized
let doc = Markdown.parse(
        """
        1. a
        3. b
        """
    )
    @test contains(sprint(term, doc), "1. ")
    @test contains(sprint(term, doc), "2. ")
    @test !contains(sprint(term, doc), "3. ")
end

# HTML output
@test md"foo *bar* baz" |> html == "<p>foo <em>bar</em> baz</p>\n"
@test md"something ***" |> html == "<p>something ***</p>\n"
@test md"# h1## " |> html == "<h1>h1##</h1>\n"
@test md"## h2 ### " |> html == "<h2>h2</h2>\n"
@test md"###### h6" |> html == "<h6>h6</h6>\n"
@test md"####### h7" |> html == "<p>####### h7</p>\n"
@test md"   >" |> html == "<blockquote>\n</blockquote>\n"
@test md"1. Hello" |> html == "<ol>\n<li><p>Hello</p>\n</li>\n</ol>\n"
@test md"* World" |> html == "<ul>\n<li><p>World</p>\n</li>\n</ul>\n"
@test md"# title *blah*" |> html == "<h1>title <em>blah</em></h1>\n"
@test md"## title *blah*" |> html == "<h2>title <em>blah</em></h2>\n"
@test md"<https://julialang.org>" |> html == """<p><a href="https://julialang.org">https://julialang.org</a></p>\n"""
@test md"<mailto://a@example.com>" |> html == """<p><a href="mailto://a@example.com">mailto://a@example.com</a></p>\n"""
@test md"<https://julialang.org/not a link>" |> html == "<p><https://julialang.org/not a link></p>\n"
@test md"""<https://julialang.org/nota
link>""" |> html == "<p><https://julialang.org/nota link></p>\n"
@test md"""Hello

---
World""" |> html == "<p>Hello</p>\n<hr />\n<p>World</p>\n"
@test md"`escape</code>`" |> html == "<p><code>escape</code></code></p>\n"

@test md"""
    code1

    code2
""" |> html == "<pre><code>code1\n\ncode2</code></pre>\n" # single code block

# @test md"""
# - Foo
#  ---
# - Bar""" |> html == "<ul>\n<li>Foo</li>\n</ul>\n<hr />\n<ul>\n<li>Bar</li>\n</ul>\n"
@test md"""
h1
===
h2
---
not
== =""" |> html == "<h1>h1</h1>\n<h2>h2</h2>\n<p>not == =</p>\n"

# Latex output
book = md"""
# Title

Some discussion

> A quote

## Section *important*

Some **bolded**

- list1
- list2
"""
@test latex(book) == "\\section{Title}\nSome discussion\n\n\\begin{quote}\nA quote\n\n\\end{quote}\n\\subsection{Section \\emph{important}}\nSome \\textbf{bolded}\n\n\\begin{itemize}\n\\item list1\n\n\n\\item list2\n\n\\end{itemize}\n"
table = md"""
 a | b
---|---
 1 | 2
"""
@test latex(table) ==
    "\\begin{tabular}\n{r | r}\na & b \\\\\n\\hline\n1 & 2 \\\\\n\\end{tabular}\n"

# mime output
let out =
    @test sprint(show, "text/plain", book) ==
        "  Title\n  ≡≡≡≡≡≡≡\n\n  Some discussion\n\n  |  A quote\n\n  Section important\n  ===================\n\n  Some bolded\n\n    •    list1\n      \n    •    list2\n      \n"
    @test sprint(show, "text/markdown", book) ==
        """
        # Title

        Some discussion

        > A quote

        ## Section *important*

        Some **bolded**

          * list1
          * list2
        """
end
let out =
    """
    <div class="markdown"><h1>Title</h1>
    <p>Some discussion</p>
    <blockquote>
    <p>A quote</p>
    </blockquote>
    <h2>Section <em>important</em></h2>
    <p>Some <strong>bolded</strong></p>
    <ul>
    <li><p>list1</p>
    </li>
    <li><p>list2</p>
    </li>
    </ul>
    </div>"""
    @test sprint(show, "text/html", book) == out
end
let out =
    """
    \\section{Title}
    Some discussion

    \\begin{quote}
    A quote

    \\end{quote}
    \\subsection{Section \\emph{important}}
    Some \\textbf{bolded}

    \\begin{itemize}
    \\item list1

    \\item list2

    \\end{itemize}
    """
    @test sprint(show, "text/latex", book) == out
end
let out =
    """
    Title
    *****

    Some discussion

        A quote

    Section *important*
    ===================

    Some **bolded**

    * list1
    * list2
    """
    @test sprint(show, "text/rst", book) == out
end

# rst rendering
for (input, output) in (
        md"foo *bar* baz"     => "foo *bar* baz\n",
        md"something ***"     => "something ***\n",
        md"# h1## "           => "h1##\n****\n\n",
        md"## h2 ### "        => "h2\n==\n\n",
        md"###### h6"         => "h6\n..\n\n",
        md"####### h7"        => "####### h7\n",
        md"   >"              => "    \n\n",
        md"1. Hello"          => "1. Hello\n",
        md"* World"           => "* World\n",
        md"``x + y``"         => ":math:`x + y`\n",
        md"# title *blah*"    => "title *blah*\n************\n\n",
        md"## title *blah*"   => "title *blah*\n============\n\n",
        md"[`x`](:func:`x`)"  => ":func:`x`\n",
        md"[`x`](:obj:`x`)"   => ":obj:`x`\n",
        md"[`x`](:ref:`x`)"   => ":ref:`x`\n",
        md"[`x`](:exc:`x`)"   => ":exc:`x`\n",
        md"[`x`](:class:`x`)" => ":class:`x`\n",
        md"[`x`](:const:`x`)" => ":const:`x`\n",
        md"[`x`](:data:`x`)"  => ":data:`x`\n",
        md"[`x`](:???:`x`)"   => "```x`` <:???:`x`>`_\n",
        md"[x](y)"            => "`x <y>`_\n",
    )
    @test rst(input) == output
end

# Interpolation / Custom types
mutable struct Reference
    ref
end

ref(x) = Reference(x)

ref(mean)

show(io::IO, m::MIME"text/plain", r::Reference) =
    print(io, "$(r.ref) (see Julia docs)")

mean_ref = md"Behaves like $(ref(mean))"
@test plain(mean_ref) == "Behaves like mean (see Julia docs)\n"
@test html(mean_ref) == "
Behaves like mean (see Julia docs)
\n"

show(io::IO, m::MIME"text/html", r::Reference) =
    Markdown.withtag(io, :a, :href=>"test") do
        Markdown.htmlesc(io, Markdown.plaininline(r))
    end
@test html(mean_ref) == "
Behaves like mean (see Julia docs)
\n"

@test md"""
````julia
foo()
````""" == md"""
```julia
foo()
```"""

# GH tables
@test md"""
    a  | b
    ---|---
    1  | 2""" == MD(Table(Any[["a","b"],
                              ["1","2"]], [:r, :r]))

@test md"""
    | a  |  b | c |
    | :-- | --: | --- |
    | d`gh`hg | hgh**jhj**ge | f |""" == MD(Table(Any[["a","b","c"],
                                                      Any[["d",Code("gh"),"hg"],
                                                          ["hgh",Bold("jhj"),"ge"],
                                                          "f"]],
                                                  [:l, :r, :r]))
@test md"""
no|table
no error
""" == MD([Paragraph(Any["no|table no error"])])

let t = """a   |   b
    :-- | --:
    1   |   2
    """
    @test Markdown.parse(t) == MD(Table(Any[Any["a", "b"], Any["1", "2"]], [:l, :r]))
end

let text =
    """
    | a   |   b |
    |:--- | ---:|
    | 1   |   2 |
    """,
    table = Markdown.parse(text)
    @test text == Markdown.plain(table)
end
let text =
    """
    | Markdown | Table |  Test |
    |:-------- |:-----:| -----:|
    | foo      | `bar` | *baz* |
    | `bar`    |  baz  | *foo* |
    """,
    table = Markdown.parse(text)
    @test text == Markdown.plain(table)
end
let text =
    """
    | a        |   b |
    |:-------- | ---:|
    | `x \\| y` |   2 |
    """,
    table = Markdown.parse(text)
    @test text == Markdown.plain(table)
end

# LaTeX extension
let in_dollars =
    """
    We have \$x^2 < x\$ whenever:

    \$|x| < 1\$

    etc.
    """,
    in_backticks =
    """
    We have ``x^2 < x`` whenever:

    ```math
    |x| < 1
    ```

    etc.
    """,
    out_plain =
    """
    We have \$x^2 < x\$ whenever:

    \$\$
    |x| < 1
    \$\$

    etc.
    """,
    out_rst =
    """
    We have :math:`x^2 < x` whenever:

    .. math::

        |x| < 1

    etc.
    """,
    out_latex =
    """
    We have \$x^2 < x\$ whenever:

    \$\$|x| < 1\$\$
    etc.

    """,
    dollars   = Markdown.parse(in_dollars),
    backticks = Markdown.parse(in_backticks),
    latex_doc = MD(
        Any[Paragraph(Any["We have ", LaTeX("x^2 < x"), " whenever:"]),
            LaTeX("|x| < 1"),
            Paragraph(Any["etc."])
    ])

    @test out_plain == Markdown.plain(dollars)
    @test out_plain == Markdown.plain(backticks)

    @test out_rst   == Markdown.rst(dollars)
    @test out_rst   == Markdown.rst(backticks)

    @test out_latex == Markdown.latex(dollars)
    @test out_latex == Markdown.latex(backticks)

    @test latex_doc == dollars
    @test latex_doc == backticks
end

# Nested backticks for inline code and math.
let t_1 = "`code` ``math`` ```code``` ````math```` `````code`````",
    t_2 = "`` `math` `` ``` `code` ``code`` ``` ```` `math` ``math`` ```math``` ````",
    t_3 = "`` ` `` ``` `` ` `` ` ` ```",
    t_4 = """`code
    over several
    lines` ``math
    over several
    lines`` ``math with
    ` some extra ` ` backticks`
    ``""",
    t_5 = "``code at end of string`",
    t_6 = "```math at end of string``"
    @test Markdown.parse(t_1) == MD(Paragraph([
        Code("code"),
        " ",
        LaTeX("math"),
        " ",
        Code("code"),
        " ",
        LaTeX("math"),
        " ",
        Code("code"),
    ]))
    @test Markdown.parse(t_2) == MD(Paragraph([
        LaTeX("`math`"),
        " ",
        Code("`code` ``code``"),
        " ",
        LaTeX("`math` ``math`` ```math```"),
    ]))
    @test Markdown.parse(t_3) == MD(Paragraph([
        LaTeX("`"),
        " ",
        Code("`` ` `` ` `"),
    ]))
    @test Markdown.parse(t_4) == MD(Paragraph([
        Code("code over several lines"),
        " ",
        LaTeX("math over several lines"),
        " ",
        LaTeX("math with ` some extra ` ` backticks`")
    ]))
    @test Markdown.parse(t_5) == MD(Paragraph([
        "`",
        Code("code at end of string"),
    ]))
    @test Markdown.parse(t_6) == MD(Paragraph([
        "`",
        LaTeX("math at end of string"),
    ]))
end

# Admonitions.
let t_1 =
        """
        # Foo

        !!! note

        !!! warning "custom title"

        ## Bar

        !!! danger ""

        !!!
        """,
    t_2 =
        """
        !!! note
            foo bar baz

        !!! warning "custom title"
            - foo
            - bar
            - baz

            foo bar baz

        !!! danger ""

            ```
            foo
            ```

                bar

            # baz
        """,
    m_1 = Markdown.parse(t_1),
    m_2 = Markdown.parse(t_2)

    # Content Tests.

    @test isa(m_1.content[2], Markdown.Admonition)
    @test m_1.content[2].category == "note"
    @test m_1.content[2].title == "Note"
    @test m_1.content[2].content == []

    @test isa(m_1.content[3], Markdown.Admonition)
    @test m_1.content[3].category == "warning"
    @test m_1.content[3].title == "custom title"
    @test m_1.content[3].content == []

    @test isa(m_1.content[5], Markdown.Admonition)
    @test m_1.content[5].category == "danger"
    @test m_1.content[5].title == ""
    @test m_1.content[5].content == []

    @test isa(m_1.content[6], Markdown.Paragraph)

    @test isa(m_2.content[1], Markdown.Admonition)
    @test m_2.content[1].category == "note"
    @test m_2.content[1].title == "Note"
    @test isa(m_2.content[1].content[1], Markdown.Paragraph)

    @test isa(m_2.content[2], Markdown.Admonition)
    @test m_2.content[2].category == "warning"
    @test m_2.content[2].title == "custom title"
    @test isa(m_2.content[2].content[1], Markdown.List)
    @test isa(m_2.content[2].content[2], Markdown.Paragraph)

    @test isa(m_2.content[3], Markdown.Admonition)
    @test m_2.content[3].category == "danger"
    @test m_2.content[3].title == ""
    @test isa(m_2.content[3].content[1], Markdown.Code)
    @test isa(m_2.content[3].content[2], Markdown.Code)
    @test isa(m_2.content[3].content[3], Markdown.Header{1})

    # Rendering Tests.
    let out = Markdown.plain(m_1),
        expected =
            """
            # Foo

            !!! note
            \n\n
            !!! warning "custom title"
            \n\n
            ## Bar

            !!! danger ""
            \n\n
            !!!
            """
        @test out == expected
    end
    let out = Markdown.rst(m_1),
        expected =
            """
            Foo
            ***
            \n
            .. note::
            \n\n
            .. warning:: custom title
            \n\n
            Bar
            ===
            \n
            .. danger::
            \n\n
            !!!
            """
        @test out == expected
    end
    let out = Markdown.latex(m_1),
        expected =
            """
            \\section{Foo}
            \\begin{quote}
            \\textbf{note}

            Note

            \\end{quote}
            \\begin{quote}
            \\textbf{warning}

            custom title

            \\end{quote}
            \\subsection{Bar}
            \\begin{quote}
            \\textbf{danger}
            \n\n
            \\end{quote}
            !!!

            """
        @test out == expected
    end
    let out = Markdown.html(m_1),
        expected =
            """
            <h1>Foo</h1>
            <div class="admonition note"><p class="admonition-title">Note</p></div>
            <div class="admonition warning"><p class="admonition-title">custom title</p></div>
            <h2>Bar</h2>
            <div class="admonition danger"><p class="admonition-title"></p></div>
            <p>!!!</p>
            """
        @test out == expected
    end

    let out = Markdown.plain(m_2),
        expected =
            """
            !!! note
                foo bar baz

            !!! warning "custom title"
                  * foo
                  * bar
                  * baz

                foo bar baz

            !!! danger ""
                ```
                foo
                ```

                ```
                bar
                ```

                # baz

            """
        @test out == expected
    end
    let out = Markdown.rst(m_2),
        expected =
            """
            .. note::
               foo bar baz

            .. warning:: custom title
               * foo
               * bar
               * baz

               foo bar baz

            .. danger::
               .. code-block:: julia

                   foo

               .. code-block:: julia

                   bar

               baz
               ***

            """
        @test out == expected
    end
end

# Nested Lists.
let text =
        """
        1. A paragraph
           with two lines.

               indented code

           > A block quote.

        - one

         two

        - one

          two

        - baz

        + ```
          foo
          ```

        1. foo
        2. bar
        3. baz
        """,
    md = Markdown.parse(text)

    # Content and structure tests.

    @test length(md.content) == 6
    @test length(md.content[1].items) == 1
    @test length(md.content[1].items[1]) == 3
    @test isa(md.content[1].items[1][1], Markdown.Paragraph)
    @test isa(md.content[1].items[1][2], Markdown.Code)
    @test isa(md.content[1].items[1][3], Markdown.BlockQuote)
    @test length(md.content[2].items) == 1
    @test isa(md.content[2].items[1][1], Markdown.Paragraph)
    @test isa(md.content[3], Markdown.Paragraph)
    @test length(md.content[4].items) == 1
    @test isa(md.content[4].items[1][1], Paragraph)
    @test isa(md.content[4].items[1][2], Paragraph)
    @test length(md.content[5].items) == 2
    @test isa(md.content[5].items[1][1], Markdown.Paragraph)
    @test isa(md.content[5].items[2][1], Markdown.Code)
    @test length(md.content[6].items) == 3
    @test md.content[6].items[1][1].content[1] == "foo"
    @test md.content[6].items[2][1].content[1] == "bar"
    @test md.content[6].items[3][1].content[1] == "baz"

    # Rendering tests.
    let expected =
            """
            1. A paragraph with two lines.

                ```
                indented code
                ```

                > A block quote.

              * one

            two

              * one

                two

              * baz
              * ```
                foo
                ```

            1. foo
            2. bar
            3. baz
            """
        @test expected == Markdown.plain(md)
    end

    let expected =
            """

A paragraph with two lines.

indented code

A block quote.

one

two

one

two

baz

foo

foo

bar

baz

            """
        @test expected == Markdown.html(md)
    end

    let expected =
            """
            1. A paragraph with two lines.

               .. code-block:: julia

                   indented code

                   A block quote.

            * one

            two

            * one

              two

            * baz
            * .. code-block:: julia

                  foo

            1. foo
            2. bar
            3. baz
            """
        @test expected == Markdown.rst(md)
    end
end

# Ordered list starting number.
let text =
        """
        42. foo
        43. bar

        1. foo
        2. bar

        - foo
        - bar
        """,
    md = Markdown.parse(text)

    @test md.content[1].ordered == 42
    @test md.content[2].ordered == 1
    @test md.content[3].ordered == -1

    let expected =
            """

foo

bar

foo

bar

foo

bar

            """
        @test expected == Markdown.html(md)
    end

    let expected =
            """
            \\begin{itemize}
            \\item[42. ] foo

            \\item[43. ] bar

            \\end{itemize}
            \\begin{itemize}
            \\item[1. ] foo

            \\item[2. ] bar

            \\end{itemize}
            \\begin{itemize}
            \\item foo

            \\item bar

            \\end{itemize}
            """
        @test expected == Markdown.latex(md)
    end
end

# different output depending on whether color is requested:
let buf = IOBuffer()
    show(buf, "text/plain", md"*emph*")
    @test String(take!(buf)) == "  emph\n"
    show(buf, "text/markdown", md"*emph*")
    @test String(take!(buf)) == "*emph*\n"
    show(IOContext(buf, :color=>true), "text/plain", md"*emph*")
    @test String(take!(buf)) == "  \e[4memph\e[24m\n"
end

# table rendering with term #25213
t = """
    a   |   b
    :-- | --:
    1   |   2"""
@test sprint(Markdown.term, Markdown.parse(t), 0) == "a b\n– –\n1 2\n"

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
# Autogenerated wrapper script for Telescope_jll for x86_64-w64-mingw32-cxx11
export libtelescope, telescope_h

using SDL2_jll
using SDL2_image_jll
using SDL2_mixer_jll
using SDL2_net_jll
using SDL2_ttf_jll
using Shaderc_jll
JLLWrappers.@generate_wrapper_header("Telescope")
JLLWrappers.@declare_library_product(libtelescope, "libtelescope.dll")
JLLWrappers.@declare_file_product(telescope_h)
function __init__()
    JLLWrappers.@generate_init_header(SDL2_jll, SDL2_image_jll, SDL2_mixer_jll, SDL2_net_jll, SDL2_ttf_jll, Shaderc_jll)
    JLLWrappers.@init_library_product(
        libtelescope,
        "bin\\libtelescope.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        telescope_h,
        "include\\telescope.h",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Base: copy, copyto!, ==
import Random: rand, seed!
import RandomNumbers: AbstractRNG, gen_seed, split_uint, seed_type

"""
```julia
AbstractXoroshiro64 <: AbstractRNG{UInt32}
```

The base abstract type for `Xoroshiro64Star` and `Xoroshiro64StarStar`.
"""
abstract type AbstractXoroshiro64 <: AbstractRNG{UInt32} end

for (star, starstar) in (
        (true, false),
        (false, true),
    )
    rng_name = Symbol(string("Xoroshiro64", star ? "Star" : starstar ? "StarStar" : ""))
    @eval begin
        mutable struct $rng_name <: AbstractXoroshiro64
            x::UInt32
            y::UInt32
            function $rng_name(seed::NTuple{2, UInt32}=gen_seed(UInt32, 2))
                r = new(0, 0)
                seed!(r, seed)
                r
            end
        end

        $rng_name(seed::Integer) = $rng_name(split_uint(seed % UInt64))

        @inline function xorshift_next(r::$rng_name)
            p = r.x * 0x9E3779BB
            $(starstar ? :(p = xorshift_rotl(p, 5) * (5 % UInt32)) : nothing)
            s1 = r.y ⊻ r.x
            r.x = xorshift_rotl(r.x, 26) ⊻ s1 ⊻ (s1 << 9)
            r.y = xorshift_rotl(s1, 13)
            p
        end
    end
end

@inline seed_type(::Type{T}) where T <: AbstractXoroshiro64 = NTuple{2, UInt32}

function copyto!(dest::T, src::T) where T <: AbstractXoroshiro64
    dest.x = src.x
    dest.y = src.y
    dest
end

copy(src::T) where T <: AbstractXoroshiro64 = copyto!(T(), src)

==(r1::T, r2::T) where T <: AbstractXoroshiro64 = r1.x == r2.x && r1.y == r2.y

seed!(r::AbstractXoroshiro64, seed::Integer) = seed!(r, split_uint(seed % UInt64))
function seed!(r::AbstractXoroshiro64, seed::NTuple{2, UInt32}=gen_seed(UInt32, 2))
    all(==(0), seed) && error("0 cannot be the seed")
    r.x = seed[1]
    r.y = seed[2]
    xorshift_next(r)
    xorshift_next(r)
    r
end

@inline rand(r::AbstractXoroshiro64, ::Type{UInt32}) = xorshift_next(r)


text:
module JANAF
using CSV # used for reading in JANAF tables

# Export useful functions
export JanafDict, Cp, S, Hs, Hf, Gf, Kf

"""
    JanafDict(path,specie_array)

Read in each specie in `specie_array` located in
JANAF table folder `path`, and place tables into
a dictionary in memory.

Tables are in `.txt` format, downloaded each from:
https://janaf.nist.gov/

Standard reference states are at 298.15K and 0.1MPa.

# Examples
Read in N2 and O2 JANAF JANAF tables:
```julia-repl
julia> path = "/my/absolute/path/";
julia> specie_array = ["N2","O2"];
julia> jd = JanafDict(path,specie_array);
```
"""
function JanafDict(path::String,specie_array::Array{String,1})
	# Initialize dictionary
	jdict = Dict()

	# TODO
	# Need to find a way to convert user-written inputs to
	# the odd JANAF title format, i.e. N-023 vs N2.

	# Loop through each specie
	for specie in specie_array
		# Read in data, convert to array
		data = Array(CSV.read(path*specie*".txt",delim='\t',ignorerepeated=true,header=2))

		# Loop through thermo quantities
		for i=2:8
			# Convert "INFINITE" string and create numeric array
			if data[1,i] == "INFINITE"
				data[1,i] = "1.0e12" # create "infinity" for interpolation
				data[:,i] = parse.(Float64,data[:,i]) # convert strings to floats for vector

			end # end comparison if
		end # end loop through thermo quantities

		# Typecast array to float
		data = Float64.(data)

		# Place data into dictionary
		jdict[specie] = data

	end # end specie loop

	# Return back parsed dictionary
	return jdict

end # end JanafDict()

"""
	interp(table,T)

If `table` is a n x 2 table with `table[:,1]` containing
temperatures and `table[:,2]` containing the property of
interest, `interp(table,T)` will interpolate `table[:,2]`
using the input `T` via a binary search and linear
interpolation between the nearest bins in the table.
"""
function interp(table::Array{Float64,2}, T::Number)
    low::Int64 = 1 # initialize low index
    high::Int64 = size(table,1) # initialize high index
    binlow::Int64 = low # Initialize low bin
    binhigh::Int64 = high # Initilize high bin

    # Perform binary search:
    while low <= high
        mid = low + div((high-low),2) # find mid, integer divide

        # Check if T is at mid, then return exact value:
        if table[mid,1] == T
            return table[mid,2]
        # Else if T>midpoint T, shift up low bound:
        elseif table[mid,1] < T
            low = mid + 1
        # Else if T<midpoint T, shift down the upper bound:
        elseif table[mid,1] > T
            high = mid - 1
        end

        # Test for out of bounds
        if table[end,1] < T
            println("WARNING: Temperature $T is larger than bound $(table[end,1])!")
            return table[end,2]
        end

        # Test for negative temperatures
        T < 0 && error("Temperature error: cannot input negative temperatures!")

        # For last iteration where location is found, the
        #   high and low bounds cross and stop the loop.
        binlow = high # set high bound to the lower bin
        binhigh = low # set low bound to the higher bin
    end

	# Return linear interpolation between bins at desired temperature
    return table[binlow,2] + (table[binhigh,2] - table[binlow,2])/(table[binhigh,1] - table[binlow,1])*(T - table[binlow,1])

end # end interp()

"""
    Cp(jd,specie,T)

Specific heat of `specie` at constant pressure
at `T` temperature in kelvin. Interpolates JanafDict
given as `jd`. Returns specific heat in J/K/mol.

# Example
```julia-repl
julia> Cp(jd,"N2",405)
29.2628
```
"""
function Cp(jd,specie::String,T::Number)
	return interp(jd[specie][:,1:2],T)
end # end Cp()

"""
   S(jd,specie,T)

Entropy of `specie` at `T` temperature in kelvin.
Interpolates JanafDict given as `jd`. Returns entropy
in J/K/mol.

# Example
```julia-repl
julia> S(jd,"N2",802)
221.0918
```
"""
function S(jd,specie::String,T::Number)
	return interp(jd[specie][:,1:2:3],T)
end # end S()

"""
   Hs(jd,specie,T)

Sensible enthalpy of `specie` at `T` temperature
in kelvin. Interpolates JanafDict given as `jd`.
Returns sensible enthalpy in kJ/mol. In JANAF
table as H-H0(Tr).

# Example
```julia-repl
julia> Hs(jd,"N2",506)
6.08998
```
"""
function Hs(jd,specie::String,T::Number)
	return interp(jd[specie][:,1:4:5],T)
end # end Hs()

"""
   Hf(jd,specie,T)

Standard enthalpy of formation of `specie` at
`T` temperature in kelvin. Interpolates
JanafDict given as `jd`. Returns formation
enthalpy in kJ/mol. In JANAF table as ΔfH0.

# Example
```julia-repl
julia> Hf(jd,"H2O",802)
-245.9564
```
"""
function Hf(jd,specie::String,T::Number)
	return interp(jd[specie][:,1:5:6],T)
end # end Hf()

"""
   Gf(jd,specie,T)

Standard Gibbs free energy of formation of
`specie` at `T` temperature in kelvin.
Interpolates JanafDict given as `jd`. Returns
formation enthalpy in kJ/mol. In JANAF table
as ΔfG0.

# Example
```julia-repl
julia> Gf(jd,"H2O",370)
-225.2807
```
"""
function Gf(jd,specie::String,T::Number)
	return interp(jd[specie][:,1:6:7],T)
end # end Gf()

"""
   Kf(jd,specie,T)

Equilibrium constant of `specie` at `T`
temperature in kelvin. Interpolates JanafDict
given as `jd`. Returns equilibrium constant in
kJ/mol. In JANAF table as log(Kf), and this
function returns Kf.

# Example
```julia-repl
julia> Kf(jd,"H2O",623)
8.506090518593942e17
```
"""
function Kf(jd,specie::String,T::Number)
	tmp = interp(jd[specie][:,1:7:8],T)
	return 10^tmp
end # end Kf()

end # module


text:
# convertion of time units
using Test
using Dates
using NCDatasets

filename = tempname()

for (timeunit,factor) in [("days",1),("hours",24),("minutes",24*60),("seconds",24*60*60)]

    NCDatasets.NCDataset(filename,"c") do ds
        NCDatasets.defDim(ds,"time",3)
        v = NCDatasets.defVar(ds,"time",Float64,("time",), attrib = [
            "units" => "$(timeunit) since 2000-01-01 00:00:00"])
        v[:] = [DateTime(2000,1,2), DateTime(2000,1,3), DateTime(2000,1,4)]
        #v.var[:] = [1.,2.,3.]

        # write "scalar" value
        v[3] = DateTime(2000,1,5)
        @test v[3] == DateTime(2000,1,5)

        # time origin
        v[3] = 0
        @test v[3] == DateTime(2000,1,1)
    end

    NCDatasets.NCDataset(filename,"r") do ds
        v2 = ds["time"].var[:]
        @test v2[1] == 1. * factor

        v2 = ds["time"][:]
        @test v2[1] == DateTime(2000,1,2)
    end

    rm(filename)
end

NCDatasets.NCDataset(filename,"c") do ds
    NCDatasets.defDim(ds,"time",3)

    v2 = NCDatasets.defVar(ds,"time2",
                           [DateTime(2000,1,2), DateTime(2000,1,3), DateTime(2000,1,4)],("time",), attrib = [
                               "units" => NCDatasets.CFTime.DEFAULT_TIME_UNITS
                           ])

    @test v2[:] == [DateTime(2000,1,2), DateTime(2000,1,3), DateTime(2000,1,4)]
    @test v2.attrib["units"] == NCDatasets.CFTime.DEFAULT_TIME_UNITS
end

# test fill-value in time axis
filename = tempname()
NCDatasets.NCDataset(filename,"c") do ds
    NCDatasets.defDim(ds,"time",3)
    v = NCDatasets.defVar(ds,"time",Float64,("time",), attrib = [
        "units" => "days since 2000-01-01 00:00:00",
        "_FillValue" => -99999.])
    v[:] = [DateTime(2000,1,2), DateTime(2000,1,3), missing]
    # load a "scalar" value
    @test v[1] == DateTime(2000,1,2)
end
rm(filename)

# test fill-value in time axis
filename = tempname()
NCDatasets.NCDataset(filename,"c") do ds
    NCDatasets.defDim(ds,"time",3)
    v = NCDatasets.defVar(ds,"time",Float64,("time",), attrib = [
        "units" => "days since 2000-01-01 00:00:00",
        "_FillValue" => -99999.])
    v[:] = [1.,2.,3.]
    # load a "scalar" value
    @test v[1] == DateTime(2000,1,2)
end
rm(filename)

# test non-standard calendars
filename = tempname()
NCDatasets.NCDataset(filename,"c") do ds
    NCDatasets.defDim(ds,"time",3)
    v = @test_logs (:warn,r".*unsupported.*") NCDatasets.defVar(ds,"time",Float64,("time",), attrib = [
        "units" => "days since 2000-01-01 00:00:00",
        "calendar" => "I_made_this_up"])
    v.var[:] = [1.,2.,3.]
    # load a "scalar" value
    @test v[1] == 1.
end
rm(filename)


text:
using ..Turing
using ..Bijectors
using LinearAlgebra

import ..AbstractMCMC: AbstractSampler
import ..DynamicPPL
import ..DynamicPPL: Model, AbstractContext, VarInfo, AbstractContext, VarName,
    _getindex, getsym, getfield, settrans!,  setorder!,
    get_and_set_val!, istrans, tilde, dot_tilde, get_vns_and_dist
import .Optim
import .Optim: optimize
import ..ForwardDiff
import NamedArrays
import StatsBase
import Printf

struct MLE end
struct MAP end

"""
    OptimizationContext{C<:AbstractContext} <: AbstractContext

The `OptimizationContext` transforms variables to their constrained space, but
does not use the density with respect to the transformation. This context is
intended to allow an optimizer to sample in R^n freely.
"""
struct OptimizationContext{C<:AbstractContext} <: AbstractContext
    context::C
end

# assume
function DynamicPPL.tilde(rng, ctx::OptimizationContext, spl, dist, vn::VarName, inds, vi)
    return DynamicPPL.tilde(ctx, spl, dist, vn, inds, vi)
end

function DynamicPPL.tilde(ctx::OptimizationContext{<:LikelihoodContext}, spl, dist, vn::VarName, inds, vi)
    r = vi[vn]
    return r, 0
end

function DynamicPPL.tilde(ctx::OptimizationContext, spl, dist, vn::VarName, inds, vi)
    r = vi[vn]
    return r, Distributions.logpdf(dist, r)
end

# observe
function DynamicPPL.tilde(rng, ctx::OptimizationContext, sampler, right, left, vi)
    return DynamicPPL.tilde(ctx, sampler, right, left, vi)
end

function DynamicPPL.tilde(ctx::OptimizationContext{<:PriorContext}, sampler, right, left, vi)
    return 0
end

function DynamicPPL.tilde(ctx::OptimizationContext, sampler, dist, value, vi)
    return Distributions.logpdf(dist, value)
end

# dot assume
function DynamicPPL.dot_tilde(rng, ctx::OptimizationContext, sampler, right, left, vn::VarName, inds, vi)
    return DynamicPPL.dot_tilde(ctx, sampler, right, left, vn, inds, vi)
end

function DynamicPPL.dot_tilde(ctx::OptimizationContext{<:LikelihoodContext}, sampler, right, left, vn::VarName, _, vi)
    vns, dist = get_vns_and_dist(right, left, vn)
    r = getval(vi, vns)
    return r, 0
end

function DynamicPPL.dot_tilde(ctx::OptimizationContext, sampler, right, left, vn::VarName, _, vi)
    vns, dist = get_vns_and_dist(right, left, vn)
    r = getval(vi, vns)
    return r, loglikelihood(dist, r)
end

# dot observe
function DynamicPPL.dot_tilde(ctx::OptimizationContext{<:PriorContext}, sampler, right, left, vn, _, vi)
    return 0
end

function DynamicPPL.dot_tilde(ctx::OptimizationContext{<:PriorContext}, sampler, right, left, vi)
    return 0
end

function DynamicPPL.dot_tilde(ctx::OptimizationContext, sampler, right, left, vn, _, vi)
    vns, dist = get_vns_and_dist(right, left, vn)
    r = getval(vi, vns)
    return loglikelihood(dist, r)
end

function DynamicPPL.dot_tilde(ctx::OptimizationContext, sampler, dists, value, vi)
    return sum(Distributions.logpdf.(dists, value))
end

function getval(
    vi,
    vns::AbstractVector{<:VarName},
)
    r = vi[vns]
    return r
end

function getval(
    vi,
    vns::AbstractArray{<:VarName},
)
    r = reshape(vi[vec(vns)], size(vns))
    return r
end

"""
    OptimLogDensity{M<:Model,C<:Context,V<:VarInfo}

A struct that stores the log density function of a `DynamicPPL` model.
"""
struct OptimLogDensity{M<:Model,C<:AbstractContext,V<:VarInfo}
    "A `DynamicPPL.Model` constructed either with the `@model` macro or manually."
    model::M
    "A `DynamicPPL.AbstractContext` used to evaluate the model. `LikelihoodContext` or `DefaultContext` are typical for MAP/MLE."
    context::C
    "A `DynamicPPL.VarInfo` struct that will be used to update model parameters."
    vi::V
end

"""
    OptimLogDensity(model::Model, context::AbstractContext)

Create a callable `OptimLogDensity` struct that evaluates a model using the given `context`.
"""
function OptimLogDensity(model::Model, context::AbstractContext)
    init = VarInfo(model)
    return OptimLogDensity(model, context, init)
end

"""
    (f::OptimLogDensity)(z)

Evaluate the log joint (with `DefaultContext`) or log likelihood (with `LikelihoodContext`)
at the array `z`.
"""
function (f::OptimLogDensity)(z)
    spl = DynamicPPL.SampleFromPrior()

    varinfo = DynamicPPL.VarInfo(f.vi, spl, z)
    f.model(varinfo, spl, f.context)
    return -DynamicPPL.getlogp(varinfo)
end

function (f::OptimLogDensity)(F, G, H, z)
    # Throw an error if a second order method was used.
    if H !== nothing
        error("Second order optimization is not yet supported.")
    end

    spl = DynamicPPL.SampleFromPrior()

    if G !== nothing
        # Calculate log joint and the gradient
        l, g = gradient_logp(
            z,
            DynamicPPL.VarInfo(f.vi, spl, z),
            f.model,
            spl,
            f.context
        )

        # Use the negative gradient because we are minimizing.
        G[:] = -g

        # If F is something, return that since we already have the
        # log joint.
        if F !== nothing
            F = -l
            return F
        end
    end

    # No gradient necessary, just return the log joint.
    if F !== nothing
        F = f(z)
        return F
    end

    return nothing
end

"""
    ModeResult{
        V<:NamedArrays.NamedArray,
        M<:NamedArrays.NamedArray,
        O<:Optim.MultivariateOptimizationResults,
        S<:NamedArrays.NamedArray
    }

A wrapper struct to store various results from a MAP or MLE estimation.
"""
struct ModeResult{
    V<:NamedArrays.NamedArray,
    O<:Optim.MultivariateOptimizationResults,
    M<:OptimLogDensity
} <: StatsBase.StatisticalModel
    "A vector with the resulting point estimates."
    values :: V
    "The stored Optim.jl results."
    optim_result :: O
    "The final log likelihood or log joint, depending on whether `MAP` or `MLE` was run."
    lp :: Float64
    "The evaluation function used to calculate the output."
    f :: M
end
#############################
# Various StatsBase methods #
#############################

function Base.show(io::IO, ::MIME"text/plain", m::ModeResult)
    print(io, "ModeResult with maximized lp of ")
    Printf.@printf(io, "%.2f", m.lp)
    println(io)
    show(io, m.values)
end

function Base.show(io::IO, m::ModeResult)
    show(io, m.values.array)
end

function StatsBase.coeftable(m::ModeResult)
    # Get columns for coeftable.
    terms = StatsBase.coefnames(m)
    estimates = m.values.array[:,1]
    stderrors = StatsBase.stderror(m)
    tstats = estimates ./ stderrors

    StatsBase.CoefTable([estimates, stderrors, tstats], ["estimate", "stderror", "tstat"], terms)
end

function StatsBase.informationmatrix(m::ModeResult; hessian_function=ForwardDiff.hessian, kwargs...)
    # Calculate Hessian and information matrix.

    # Convert the values to their unconstrained states to make sure the
    # Hessian is computed with respect to the untransformed parameters.
    spl = DynamicPPL.SampleFromPrior()

    # NOTE: This should be converted to islinked(vi, spl) after
    # https://github.com/TuringLang/DynamicPPL.jl/pull/124 goes through.
    vns = DynamicPPL._getvns(m.f.vi, spl)

    linked = DynamicPPL._islinked(m.f.vi, vns)
    linked && invlink!(m.f.vi, spl)

    # Calculate the Hessian.
    varnames = StatsBase.coefnames(m)
    H = hessian_function(m.f, m.values.array[:, 1])
    info = inv(H)

    # Link it back if we invlinked it.
    linked && link!(m.f.vi, spl)

    return NamedArrays.NamedArray(info, (varnames, varnames))
end

StatsBase.coef(m::ModeResult) = m.values
StatsBase.coefnames(m::ModeResult) = names(m.values)[1]
StatsBase.params(m::ModeResult) = StatsBase.coefnames(m)
StatsBase.vcov(m::ModeResult) = StatsBase.informationmatrix(m)
StatsBase.loglikelihood(m::ModeResult) = m.lp

####################
# Optim.jl methods #
####################

"""
    Optim.optimize(model::Model, ::MLE, args...; kwargs...)

Compute a maximum likelihood estimate of the `model`.

# Examples

```julia-repl
@model function f(x)
    m ~ Normal(0, 1)
    x ~ Normal(m, 1)
end

model = f(1.5)
mle = optimize(model, MLE())

# Use a different optimizer
mle = optimize(model, MLE(), NelderMead())
```
"""
function Optim.optimize(model::Model, ::MLE, options::Optim.Options=Optim.Options(); kwargs...)
    return _mle_optimize(model, options; kwargs...)
end
function Optim.optimize(model::Model, ::MLE, init_vals::AbstractArray, options::Optim.Options=Optim.Options(); kwargs...)
    return _mle_optimize(model, init_vals, options; kwargs...)
end
function Optim.optimize(model::Model, ::MLE, optimizer::Optim.AbstractOptimizer, options::Optim.Options=Optim.Options(); kwargs...)
    return _mle_optimize(model, optimizer, options; kwargs...)
end
function Optim.optimize(
    model::Model,
    ::MLE,
    init_vals::AbstractArray,
    optimizer::Optim.AbstractOptimizer,
    options::Optim.Options=Optim.Options();
    kwargs...
)
    return _mle_optimize(model, init_vals, optimizer, options; kwargs...)
end

function _mle_optimize(model::Model, args...; kwargs...)
    ctx = OptimizationContext(DynamicPPL.LikelihoodContext())
    return _optimize(model, OptimLogDensity(model, ctx), args...; kwargs...)
end

"""
    Optim.optimize(model::Model, ::MAP, args...; kwargs...)

Compute a maximum a posterior estimate of the `model`.

# Examples

```julia-repl
@model function f(x)
    m ~ Normal(0, 1)
    x ~ Normal(m, 1)
end

model = f(1.5)
map_est = optimize(model, MAP())

# Use a different optimizer
map_est = optimize(model, MAP(), NelderMead())
```
"""

function Optim.optimize(model::Model, ::MAP, options::Optim.Options=Optim.Options(); kwargs...)
    return _map_optimize(model, options; kwargs...)
end
function Optim.optimize(model::Model, ::MAP, init_vals::AbstractArray, options::Optim.Options=Optim.Options(); kwargs...)
    return _map_optimize(model, init_vals, options; kwargs...)
end
function Optim.optimize(model::Model, ::MAP, optimizer::Optim.AbstractOptimizer, options::Optim.Options=Optim.Options(); kwargs...)
    return _map_optimize(model, optimizer, options; kwargs...)
end
function Optim.optimize(
    model::Model,
    ::MAP,
    init_vals::AbstractArray,
    optimizer::Optim.AbstractOptimizer,
    options::Optim.Options=Optim.Options();
    kwargs...
)
    return _map_optimize(model, init_vals, optimizer, options; kwargs...)
end

function _map_optimize(model::Model, args...; kwargs...)
    ctx = OptimizationContext(DynamicPPL.DefaultContext())
    return _optimize(model, OptimLogDensity(model, ctx), args...; kwargs...)
end

"""
    _optimize(model::Model, f::OptimLogDensity, optimizer=Optim.LBFGS(), args...; kwargs...)

Estimate a mode, i.e., compute a MLE or MAP estimate.
"""
function _optimize(
    model::Model,
    f::OptimLogDensity,
    optimizer::Optim.AbstractOptimizer = Optim.LBFGS(),
    args...;
    kwargs...
)
    return _optimize(model, f, f.vi[DynamicPPL.SampleFromPrior()], optimizer, args...; kwargs...)
end

function _optimize(
    model::Model,
    f::OptimLogDensity,
    options::Optim.Options = Optim.Options(),
    args...;
    kwargs...
)
    return _optimize(model, f, f.vi[DynamicPPL.SampleFromPrior()], Optim.LBFGS(), args...; kwargs...)
end

function _optimize(
    model::Model,
    f::OptimLogDensity,
    init_vals::AbstractArray = f.vi[DynamicPPL.SampleFromPrior()],
    options::Optim.Options = Optim.Options(),
    args...;
    kwargs...
)
    return _optimize(model, f,init_vals, Optim.LBFGS(), options, args...; kwargs...)
end

function _optimize(
    model::Model,
    f::OptimLogDensity,
    init_vals::AbstractArray = f.vi[DynamicPPL.SampleFromPrior()],
    optimizer::Optim.AbstractOptimizer = Optim.LBFGS(),
    options::Optim.Options = Optim.Options(),
    args...;
    kwargs...
)
    # Do some initialization.
    spl = DynamicPPL.SampleFromPrior()

    # Convert the initial values, since it is assumed that users provide them
    # in the constrained space.
    f.vi[spl] = init_vals
    link!(f.vi, spl)
    init_vals = f.vi[spl]

    # Optimize!
    M = Optim.optimize(Optim.only_fgh!(f), init_vals, optimizer, options, args...; kwargs...)

    # Warn the user if the optimization did not converge.
    if !Optim.converged(M)
        @warn "Optimization did not converge! You may need to correct your model or adjust the Optim parameters."
    end

    # Get the VarInfo at the MLE/MAP point, and run the model to ensure
    # correct dimensionality.
    f.vi[spl] = M.minimizer
    invlink!(f.vi, spl)
    vals = f.vi[spl]
    link!(f.vi, spl)

    # Make one transition to get the parameter names.
    ts = [Turing.Inference.Transition(DynamicPPL.tonamedtuple(f.vi), DynamicPPL.getlogp(f.vi))]
    varnames, _ = Turing.Inference._params_to_array(ts)

    # Store the parameters and their names in an array.
    vmat = NamedArrays.NamedArray(vals, varnames)

    return ModeResult(vmat, M, -M.minimum, f)
end


text:
export quantum_ising_hamiltonian, lattice_average_operator

"""
    sigmax(hilbert::AbstractHilbert, site::Int)

Builds a σₓ operator acting on the `site`-th site of the Many body Hilbert
space `hilbert`.

Note: M-dimensional Hilbert spaces are treated as spin (M-1)//2 spins.
"""
function QuantumOpticsBase.sigmax(h::AbstractHilbert, i::Int)
    N = shape(h)[i]
    S = (N-1)//2 # map Fock spaces to Spins

    # Build the local operator matrix
    D   = ComplexF64[complex(sqrt(real((S + 1)*2*a - a*(a+1)))) for a=1:N-1]
    mat = diagm(1=>D, -1=>D)

    return KLocalOperatorRow(h, [i], mat)
end

"""
    sigmay(hilbert::AbstractHilbert, site::Int)

Builds a σ_y operator acting on the `site`-th site of the Many body Hilbert
space `hilbert`.

Note: Hilbert must have local space dimension 2 on that site.
"""
function QuantumOpticsBase.sigmay(h::AbstractHilbert, i::Int)
    N = shape(h)[i]
    S = (N-1)//2 # map Fock spaces to Spins

    # Build the local operator matrix
    D   = ComplexF64[1im*complex(sqrt(real((S + 1)*2*a - a*(a+1)))) for a=1:N-1]
    mat = diagm(-1=>D, 1=>-D)

    return KLocalOperatorRow(h, [i], mat)
end

"""
    sigmaz(hilbert::AbstractHilbert, site::Int)

Builds a σ_z operator acting on the `site`-th site of the Many body Hilbert
space `hilbert`.

Note: M-dimensional Hilbert spaces are treated as spin (M-1)//2 spins.
"""
function QuantumOpticsBase.sigmaz(h::AbstractHilbert, i::Int)
    N = shape(h)[i]
    S = (N-1)//2 # map Fock spaces to Spins

    # Build the local operator matrix
    D   = ComplexF64[complex(2*m) for m=S:-1:-S]
    mat = diagm(0 => D)

    return KLocalOperatorRow(h, [i], mat)
end

"""
    sigmam(hilbert::AbstractHilbert, site::Int)

Builds a σ₋ operator acting on the `site`-th site of the Many body Hilbert
space `hilbert`.

Note: M-dimensional Hilbert spaces are treated as spin (M-1)//2 spins.
"""
function QuantumOpticsBase.sigmam(h::AbstractHilbert, i::Int)
    N = shape(h)[i]
    S = (N-1)//2 # map Fock spaces to Spins

    S2 = (S+1)*S
    D = [complex(sqrt(float(S2 - m*(m-1)))) for m=S:-1:-S+1]
    mat = diagm(-1 => D)

    return KLocalOperatorRow(h, [i], mat)
end

"""
    sigmap(hilbert::AbstractHilbert, site::Int)

Builds a σ₊ operator acting on the `site`-th site of the Many body Hilbert
space `hilbert`.

Note: M-dimensional Hilbert spaces are treated as spin (M-1)//2 spins.
"""
function QuantumOpticsBase.sigmap(h::AbstractHilbert, i::Int)
    N = shape(h)[i]
    S = (N-1)//2 # map Fock spaces to Spins

    S2 = (S+1)*S
    D = ComplexF64[complex(sqrt(float(S2 - m*(m+1)))) for m=S-1:-1:-S]
    mat = diagm(1 => D)

    return KLocalOperatorRow(h, [i], mat)
end

"""
    destroy(hilbert::AbstractHilbert, site::Int)

Builds a bosonic destruction operator acting on the `site`-th site of the Many
body Hilbert space `hilbert`.

Note: spin-`m//2` spaces are treated as fock spaces with dimension `m+1`
"""
function QuantumOpticsBase.destroy(h::AbstractHilbert, i::Int)
    N = shape(h)[i] - 1

    D = complex.(sqrt.(1.:N))
    mat = diagm(1 => D)

    return KLocalOperatorRow(h, [i], mat)
end

"""
    create(hilbert::AbstractHilbert, site::Int)

Builds a bosonic creation operator acting on the `site`-th site of the Many
body Hilbert space `hilbert`.

Note: spin-`m//2` spaces are treated as fock spaces with dimension `m+1`
"""
function QuantumOpticsBase.create(h::AbstractHilbert, i::Int)
    N = shape(h)[i] - 1

    D = complex.(sqrt.(1.:N))
    mat = diagm(-1 => D)

    return KLocalOperatorRow(h, [i], mat)
end

"""
    number(hilbert::AbstractHilbert, site::Int)

Builds a bosonic creation operator acting on the `site`-th site of the Many
body Hilbert space `hilbert`.

Note: spin-`m//2` spaces are treated as fock spaces with dimension `m+1`
"""
function QuantumOpticsBase.number(h::AbstractHilbert, i::Int)
    N = shape(h)[i] - 1

    D = complex.(0.:N)
    mat = diagm(0 => D)

    return KLocalOperatorRow(h, [i], mat)
end

quantum_ising_hamiltonian(lattice, hilb; kwargs...) =
    quantum_ising_hamiltonian(STD_REAL_PREC, lattice, hilb; kwargs...)
"""
    quantum_ising_hamiltonian(lattice, g, V)

Constructs the Quantum Ising Hamiltonian defined as

```math
H = \\sum_i g/2 σ_i^{(x)} + \\V/4 sum_{i,j} σ_i^{z}σ_j^{z}
```
"""
function quantum_ising_hamiltonian(T::Type, lattice, hilb::AbstractHilbert; g, V)
    H = LocalOperator(T, hilb)

    for i=vertices(lattice)
        H += g/2.0 * sigmax(hilb, i)
    end
    for e=edges(lattice)
        i = src(e)
        j = dst(e)
        H += V/4.0 * sigmaz(hilb, i) * sigmaz(hilb, j)
    end

    return H
end

lattice_average_operator(lattice, hilb, op) =
    lattice_average_operator(STD_REAL_PREC, lattice, hilb, op)
function lattice_average_operator(T::Type, lattice, hilb::AbstractHilbert, op::Function)
    O = LocalOperator(T, hilb)
    for i=vertices(lattice)
        O += op(hilb, i)
    end
    O = O/nv(lattice)

    return O
end


text:
@testset "lanczos" begin
	@testset "givens" begin
		symOrth = KrylovMethods.symOrtho

		v     = ([0;randn()],[randn();0],randn(2))
		for k=1:3
			println("Givens rotation for v=[$(v[k][1]),$(v[k][2])]")
			c,s,r = symOrth(v[k][1],v[k][2])
			@test abs(c.^2 + s.^2 - 1) < 1e-14
			Q = [c s; -s c]
			w = [norm(v[k]);0]
			@test norm(Q*v[k]-w) < 1e-13
		end
	end
	@testset "bidiag" begin
		println("=== Testing Lanczos Bidiagonalization ===")
		A = sprandn(10,10,.4) + I
		p = randn(10)
		k = 9

		U,B,V = lanczosBidiag(A,p,k)

		@test norm(Matrix(U'*A*V-B))/norm(Matrix(A)) < 1e-15
	end

	@testset "tridiag" begin
		# test stopping for low-rank matrix
		A = randn(100,10)
		H = A*A'
		T,V = lanczosTridiag(H,ones(100),100,doReorth=true)
		@test size(T,1)==11
		@test size(V,2)==11

		# A   = matrixdepot("rosser",4)
		# Aff  = x -> A*x
		# T1,V1 = lanczosTridiag(A,ones(size(A,1)),size(A,1),tol=1e-10,doReorth=true)
		# T2,V2 = lanczosTridiag(Aff,ones(size(A,1)),size(A,1),tol=1e-10,doReorth=false)
		# sT1  = svd(Matrix(T1)).S
		# sT2  = svd(Matrix(T2)).S
		# sA   = svd(Matrix(A)).S
		#
		# @test abs(sT1[1]-sA[1])/sA[1] < 1e-4
		# @test abs(sT2[1]-sA[1])/sA[1] < 1e-4
		# B = V1'*V1
		# @test norm(B - Diagonal(diag(B))) < 1e-10
		# @test norm(diag(B)- ones(size(B,1))) < 1e-10
	end
end


text:
"""
    struct Box{T, CO} <: AbstractVolumePrimitive{T}

Volume primitive describing a three-dimensional [Box](@ref) with its surfaces
being parallel to the `xy`, `xy` and `yz` plane.

## Parametric types
* `T`: Precision type.
* `CO`: Describes whether the surface belongs to the primitive.
    It can be `ClosedPrimitive`, i.e. the surface points belong to the primitive,
    or `OpenPrimitive`, i.e. the surface points do not belong to the primitive.

## Fields
* `hX::T`: Half of the width in `x` dimension (in m).
* `hY::T`: Half of the width in `y` dimension (in m).
* `hZ::T`: Half of the width in `z` dimension (in m).
* `origin::CartesianPoint{T}`: The position of the center of the `Box`.
* `rotation::SMatrix{3,3,T,9}`: Matrix that describes a rotation of the `Box` around its `origin`.

## Definition in Configuration File

A `Box` is defined in the configuration file as part of the `geometry` field
of an object through the field `box`.

Example definitions of a `Box` looks like this:
```yaml
box:
  widths: [2, 4, 6] # => hX = 1; hY = 2; hZ = 3;
  origin: [0, 0, 0] # [x, y, z] - Optional; Default: [0, 0, 0]
  rotate: # Optional; Default: no rotation
    Z: 0
```

The halfwidths `hX`, `hY` and `hZ` can also be defined directly in the configuration file:
```yaml
box:
  halfwidths: [1, 2, 3] # => hX = 1; hY = 2; hZ = 3;
```
or
```yaml
box:
  hX: 1
  hY: 2
  hZ: 3
```

See also [Constructive Solid Geometry (CSG)](@ref).
"""
struct Box{T,CO} <: AbstractVolumePrimitive{T,CO}
    hX::T
    hY::T
    hZ::T

    origin::CartesianPoint{T}
    rotation::SMatrix{3,3,T,9}
end

#Type conversion happens here
function Box{T}(CO,hX, hY, hZ, origin, rotation) where {T}
    _hX = _csg_convert_args(T, hX)
    _hY = _csg_convert_args(T, hY)
    _hZ = _csg_convert_args(T, hZ)
    Box{T,CO}(_hX, _hY, _hZ, origin, rotation)
end

#Type promotion happens here
function Box(CO, hX::TX, hY::TY, hZ::TZ, origin::PT, rotation::ROT) where {TX, TY, TZ, PT, ROT}
    eltypes = _csg_get_promoted_eltype.((TX, TY, TZ, PT, ROT))
    T = float(promote_type(eltypes...))
    Box{T}(CO,T(hX), T(hY), T(hZ), origin, rotation)
end

function Box(::Type{CO}=ClosedPrimitive;
    hX = 1,
    hY = 1,
    hZ = 1,
    origin = zero(CartesianPoint{Int}),
    rotation = one(SMatrix{3, 3, Int, 9})
) where {CO}
    Box(CO, hX, hY, hZ, origin, rotation)
end

function Box{T}(::Type{CO}=ClosedPrimitive;
    hX = 1.0,
    hY = 1.0,
    hZ = 1.0,
    origin = zero(CartesianPoint{Float64}),
    rotation = one(SMatrix{3, 3, Float64, 9})
) where {T, CO}
    Box{T}(CO, hX, hY, hZ, origin, rotation)
end

Box{T, CO}( b::Box{T, CO}; COT = CO,
            origin::CartesianPoint{T} = b.origin,
            rotation::SMatrix{3,3,T,9} = b.rotation) where {T, CO<:Union{ClosedPrimitive, OpenPrimitive}} =
    Box{T, COT}(b.hX, b.hY, b.hZ, origin, rotation)

function _in(pt::CartesianPoint{T}, b::Box{<:Any, ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.x) <= b.hX + csgtol &&
    abs(pt.y) <= b.hY + csgtol &&
    abs(pt.z) <= b.hZ + csgtol
end

_in(pt::CartesianPoint{T}, b::Box{<:Any, OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T} =
    abs(pt.x) < b.hX - csgtol && abs(pt.y) < b.hY - csgtol && abs(pt.z) < b.hZ - csgtol

function Geometry(::Type{T}, ::Type{Box}, dict::AbstractDict, input_units::NamedTuple, transformations::Transformations{T}) where {T}
    length_unit = input_units.length
    origin = get_origin(T, dict, length_unit)
    rotation = get_rotation(T, dict, input_units.angle)
    hX, hY, hZ = if haskey(dict, "widths")
        @assert length(dict["widths"]) == 3
        _parse_value(T, dict["widths"], length_unit)/2
    elseif haskey(dict, "halfwidths")
        @assert length(dict["halfwidths"]) == 3
        _parse_value(T, dict["halfwidths"], length_unit)
    elseif haskey(dict, "hX") && haskey(dict, "hZ") && haskey(dict, "hZ")
        _parse_value(T, dict["hX"], length_unit),
        _parse_value(T, dict["hY"], length_unit),
        _parse_value(T, dict["hZ"], length_unit)
    elseif haskey(dict, "x") && haskey(dict, "y") && haskey(dict, "z")
        @warn "Deprecation warning: Detected old primitive definition for `Box`.
            Please update your configuration file to the new format
            via `widths`, `halfwidths` or `hX`, `hY` and `hZ`
            in combination with possible fields `origin` and `rotate`.
            The old definition overwrites the optional field `origin`."
        x = parse_interval_of_primitive(T, "x", dict, length_unit)
        y = parse_interval_of_primitive(T, "y", dict, length_unit)
        z = parse_interval_of_primitive(T, "z", dict, length_unit)
        μx = typeof(x) <: Real ? zero(T) : mean(x)
        μy = typeof(y) <: Real ? zero(T) : mean(y)
        μz = typeof(z) <: Real ? zero(T) : mean(z)
        origin = CartesianPoint{T}(μx, μy, μz)
        hX = typeof(x) <: Real ? x : (x[2] - x[1])/2
        hY = typeof(y) <: Real ? y : (y[2] - y[1])/2
        hZ = typeof(z) <: Real ? z : (z[2] - z[1])/2
        hX, hY, hZ
    end
    box = Box{T}(ClosedPrimitive,
        hX = hX,
        hY = hY,
        hZ = hZ,
        origin = origin,
        rotation = rotation
    )
    transform(box, transformations)
end

function Dictionary(b::Box{T})::OrderedDict{String, Any} where {T}
    dict = OrderedDict{String,Any}()
    # dict["widths"] = [2*b.hX, 2*b.hY, 2*b.hZ]
    dict["hX"] = b.hX
    dict["hY"] = b.hY
    dict["hZ"] = b.hZ
    if b.origin != zero(CartesianVector{T}) dict["origin"] = b.origin end
    if b.rotation != one(SMatrix{3,3,T,9}) dict["rotation"] = Dictionary(b.rotation) end
    OrderedDict{String,Any}("box" => dict)
end

function vertices(b::Box{T}) where {T}
    return (
        b.rotation * SVector{3,T}(-b.hX, -b.hY, -b.hZ) .+ b.origin,
        b.rotation * SVector{3,T}(+b.hX, -b.hY, -b.hZ) .+ b.origin,
        b.rotation * SVector{3,T}(+b.hX, +b.hY, -b.hZ) .+ b.origin,
        b.rotation * SVector{3,T}(-b.hX, +b.hY, -b.hZ) .+ b.origin,
        b.rotation * SVector{3,T}(-b.hX, -b.hY, +b.hZ) .+ b.origin,
        b.rotation * SVector{3,T}(+b.hX, -b.hY, +b.hZ) .+ b.origin,
        b.rotation * SVector{3,T}(+b.hX, +b.hY, +b.hZ) .+ b.origin,
        b.rotation * SVector{3,T}(-b.hX, +b.hY, +b.hZ) .+ b.origin,
    )
end

function sample(b::Box{T})::Vector{CartesianPoint{T}} where {T}
    [vertices(b)...]
end

function surfaces(b::Box{T}) where {T}
    vs = vertices(b)
    return (
        Quadrangle{T}((vs[1], vs[2], vs[3], vs[4])),
        Quadrangle{T}((vs[5], vs[6], vs[2], vs[1])),
        Quadrangle{T}((vs[8], vs[7], vs[6], vs[5])),
        Quadrangle{T}((vs[6], vs[7], vs[3], vs[2])),
        Quadrangle{T}((vs[7], vs[8], vs[4], vs[3])),
        Quadrangle{T}((vs[8], vs[5], vs[1], vs[4])),
    )
end

extremum(b::Box{T}) where {T} = hypot(b.hX, b.hY, b.hZ)


text:
"""
    volumerhs!(::Val{dim}, ::Val{N}, ::Val{nstate}, ::Val{nviscstate},
               ::Val{nauxstate}, flux!, source!, rhs::Array, Q, Qvisc, auxstate,
               vgeo, t, D, elems) where {dim, N, nstate, nviscstate,

Computational kernel: Evaluate the volume integrals on right-hand side of a
`DGBalanceLaw` semi-discretization.

See [`odefun!`](@ref) for usage.
"""
function volumerhs!(::Val{dim}, ::Val{N},
                    ::Val{nstate}, ::Val{nviscstate},
                    ::Val{nauxstate},
                    flux!, source!,
                    rhs::Array,
                    Q, Qvisc, auxstate, vgeo, t,
                    D, elems) where {dim, N, nstate, nviscstate,
                                     nauxstate}
  DFloat = eltype(Q)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  nelem = size(Q)[end]

  Q = reshape(Q, Nq, Nq, Nqk, nstate, nelem)
  Qvisc = reshape(Qvisc, Nq, Nq, Nqk, nviscstate, nelem)
  rhs = reshape(rhs, Nq, Nq, Nqk, nstate, nelem)
  vgeo = reshape(vgeo, Nq, Nq, Nqk, _nvgeo, nelem)
  auxstate = reshape(auxstate, Nq, Nq, Nqk, nauxstate, nelem)

  s_F = MArray{Tuple{3, Nq, Nq, Nqk, nstate}, DFloat}(undef)

  source! !== nothing && (l_S = MArray{Tuple{nstate}, DFloat}(undef))
  l_Q = MArray{Tuple{nstate}, DFloat}(undef)
  l_Qvisc = MArray{Tuple{nviscstate}, DFloat}(undef)
  l_aux = MArray{Tuple{nauxstate}, DFloat}(undef)

  l_F = MArray{Tuple{3, nstate}, DFloat}(undef)

  @inbounds for e in elems
    for k = 1:Nqk, j = 1:Nq, i = 1:Nq
      MJ = vgeo[i, j, k, _MJ, e]
      # MJI = vgeo[i, j, k, _MJI, e]
      ξx, ξy, ξz = vgeo[i,j,k,_ξx,e], vgeo[i,j,k,_ξy,e], vgeo[i,j,k,_ξz,e]
      ηx, ηy, ηz = vgeo[i,j,k,_ηx,e], vgeo[i,j,k,_ηy,e], vgeo[i,j,k,_ηz,e]
      ζx, ζy, ζz = vgeo[i,j,k,_ζx,e], vgeo[i,j,k,_ζy,e], vgeo[i,j,k,_ζz,e]

      for s = 1:nstate
        l_Q[s] = Q[i, j, k, s, e]
      end

      for s = 1:nviscstate
        l_Qvisc[s] = Qvisc[i, j, k, s, e]
      end

      for s = 1:nauxstate
        l_aux[s] = auxstate[i, j, k, s, e]
      end

      flux!(l_F, l_Q, l_Qvisc, l_aux, t)

      for s = 1:nstate
        s_F[1,i,j,k,s] = MJ * (ξx * l_F[1, s] + ξy * l_F[2, s] + ξz * l_F[3, s])
        s_F[2,i,j,k,s] = MJ * (ηx * l_F[1, s] + ηy * l_F[2, s] + ηz * l_F[3, s])
        s_F[3,i,j,k,s] = MJ * (ζx * l_F[1, s] + ζy * l_F[2, s] + ζz * l_F[3, s])
      end

      if source! !== nothing
        source!(l_S, l_Q, l_aux, t)

        for s = 1:nstate
          rhs[i, j, k, s, e] += l_S[s]
        end
      end
    end

    # loop of ξ-grid lines
    for s = 1:nstate, k = 1:Nqk, j = 1:Nq, i = 1:Nq
      MJI = vgeo[i, j, k, _MJI, e]
      for n = 1:Nq
        rhs[i, j, k, s, e] += MJI * D[n, i] * s_F[1, n, j, k, s]
      end
    end
    # loop of η-grid lines
    for s = 1:nstate, k = 1:Nqk, j = 1:Nq, i = 1:Nq
      MJI = vgeo[i, j, k, _MJI, e]
      for n = 1:Nq
        rhs[i, j, k, s, e] += MJI * D[n, j] * s_F[2, i, n, k, s]
      end
    end
    # loop of ζ-grid lines
    if Nqk > 1
      for s = 1:nstate, k = 1:Nqk, j = 1:Nq, i = 1:Nq
        MJI = vgeo[i, j, k, _MJI, e]
        for n = 1:Nqk
          rhs[i, j, k, s, e] += MJI * D[n, k] * s_F[3, i, j, n, s]
        end
      end
    end
  end
end

"""
    facerhs!(::Val{dim}, ::Val{N}, ::Val{nstate}, ::Val{nviscstate},
             ::Val{nauxstate}, numerical_flux!,
             numerical_boundary_flux!, rhs::Array, Q, Qvisc, auxstate,
             vgeo, sgeo, t, vmapM, vmapP, elemtobndy,
             elems) where {dim, N, nstate, nviscstate, nauxstate}

Computational kernel: Evaluate the surface integrals on right-hand side of a
`DGBalanceLaw` semi-discretization.

See [`odefun!`](@ref) for usage.
"""
function facerhs!(::Val{dim}, ::Val{N},
                  ::Val{nstate}, ::Val{nviscstate},
                  ::Val{nauxstate},
                  numerical_flux!,
                  numerical_boundary_flux!,
                  rhs::Array, Q, Qvisc, auxstate,
                  vgeo, sgeo,
                  t, vmapM, vmapP, elemtobndy,
                  elems) where {dim, N, nstate, nviscstate, nauxstate}
  DFloat = eltype(Q)

  if dim == 1
    Np = (N+1)
    Nfp = 1
    nface = 2
  elseif dim == 2
    Np = (N+1) * (N+1)
    Nfp = (N+1)
    nface = 4
  elseif dim == 3
    Np = (N+1) * (N+1) * (N+1)
    Nfp = (N+1) * (N+1)
    nface = 6
  end

  l_QM = MArray{Tuple{nstate}, DFloat}(undef)
  l_QviscM = MArray{Tuple{nviscstate}, DFloat}(undef)
  l_auxM = MArray{Tuple{nauxstate}, DFloat}(undef)

  l_QP = MArray{Tuple{nstate}, DFloat}(undef)
  l_QviscP = MArray{Tuple{nviscstate}, DFloat}(undef)
  l_auxP = MArray{Tuple{nauxstate}, DFloat}(undef)

  l_F = MArray{Tuple{nstate}, DFloat}(undef)

  @inbounds for e in elems
    for f = 1:nface
      for n = 1:Nfp
        nM = (sgeo[_nx, n, f, e], sgeo[_ny, n, f, e], sgeo[_nz, n, f, e])
        sMJ, vMJI = sgeo[_sMJ, n, f, e], sgeo[_vMJI, n, f, e]
        idM, idP = vmapM[n, f, e], vmapP[n, f, e]

        eM, eP = e, ((idP - 1) ÷ Np) + 1
        vidM, vidP = ((idM - 1) % Np) + 1,  ((idP - 1) % Np) + 1

        # Load minus side data
        for s = 1:nstate
          l_QM[s] = Q[vidM, s, eM]
        end

        for s = 1:nviscstate
          l_QviscM[s] = Qvisc[vidM, s, eM]
        end

        for s = 1:nauxstate
          l_auxM[s] = auxstate[vidM, s, eM]
        end

        # Load plus side data
        for s = 1:nstate
          l_QP[s] = Q[vidP, s, eP]
        end

        for s = 1:nviscstate
          l_QviscP[s] = Qvisc[vidP, s, eP]
        end

        for s = 1:nauxstate
          l_auxP[s] = auxstate[vidP, s, eP]
        end

        bctype =
            numerical_boundary_flux! === nothing ? 0 : elemtobndy[f, e]
        if bctype == 0
          numerical_flux!(l_F, nM, l_QM, l_QviscM, l_auxM, l_QP, l_QviscP,
                          l_auxP, t)
        else numerical_boundary_flux!(l_F, nM, l_QM, l_QviscM, l_auxM, l_QP,
                                      l_QviscP, l_auxP, bctype, t)
        end

        #Update RHS
        for s = 1:nstate
          rhs[vidM, s, eM] -= vMJI * sMJ * l_F[s]
        end
      end
    end
  end
end

function volumeviscterms!(::Val{dim}, ::Val{N}, ::Val{nstate},
                          ::Val{states_grad}, ::Val{ngradstate},
                          ::Val{nviscstate}, ::Val{nauxstate},
                          viscous_transform!, gradient_transform!, Q::Array,
                          Qvisc, auxstate, vgeo, t, D,
                          elems) where {dim, N, states_grad, ngradstate,
                                        nviscstate, nstate, nauxstate}
  DFloat = eltype(Q)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  nelem = size(Q)[end]
  ngradtransformstate = length(states_grad)

  Q = reshape(Q, Nq, Nq, Nqk, nstate, nelem)
  Qvisc = reshape(Qvisc, Nq, Nq, Nqk, nviscstate, nelem)
  auxstate = reshape(auxstate, Nq, Nq, Nqk, nauxstate, nelem)
  vgeo = reshape(vgeo, Nq, Nq, Nqk, _nvgeo, nelem)

  s_G = MArray{Tuple{Nq, Nq, Nqk, ngradstate}, DFloat}(undef)

  l_Q = MArray{Tuple{ngradtransformstate}, DFloat}(undef)
  l_aux = MArray{Tuple{nauxstate}, DFloat}(undef)
  l_G = MArray{Tuple{ngradstate}, DFloat}(undef)
  l_Qvisc = MArray{Tuple{nviscstate}, DFloat}(undef)
  l_gradG = MArray{Tuple{3, ngradstate}, DFloat}(undef)

  @inbounds for e in elems
    for k = 1:Nqk, j = 1:Nq, i = 1:Nq
      for s = 1:ngradtransformstate
        l_Q[s] = Q[i, j, k, states_grad[s], e]
      end

      for s = 1:nauxstate
        l_aux[s] = auxstate[i, j, k, s, e]
      end

      gradient_transform!(l_G, l_Q, l_aux, t)
      for s = 1:ngradstate
        s_G[i, j, k, s] = l_G[s]
      end
    end

    # Compute gradient of each state
    for k = 1:Nqk, j = 1:Nq, i = 1:Nq
      ξx, ξy, ξz = vgeo[i,j,k,_ξx,e], vgeo[i,j,k,_ξy,e], vgeo[i,j,k,_ξz,e]
      ηx, ηy, ηz = vgeo[i,j,k,_ηx,e], vgeo[i,j,k,_ηy,e], vgeo[i,j,k,_ηz,e]
      ζx, ζy, ζz = vgeo[i,j,k,_ζx,e], vgeo[i,j,k,_ζy,e], vgeo[i,j,k,_ζz,e]

      for s = 1:ngradtransformstate
        l_Q[s] = Q[i, j, k, states_grad[s], e]
      end

      for s = 1:ngradstate
        Gξ = Gη = Gζ = zero(DFloat)
        for n = 1:Nq
          Gξ += D[i, n] * s_G[n, j, k, s]
          Gη += D[j, n] * s_G[i, n, k, s]
          dim == 3 && (Gζ += D[k, n] * s_G[i, j, n, s])
        end
        l_gradG[1, s] = ξx * Gξ + ηx * Gη + ζx * Gζ
        l_gradG[2, s] = ξy * Gξ + ηy * Gη + ζy * Gζ
        l_gradG[3, s] = ξz * Gξ + ηz * Gη + ζz * Gζ
      end

      viscous_transform!(l_Qvisc, l_gradG, l_Q, l_aux, t)

      for s = 1:nviscstate
        Qvisc[i, j, k, s, e] = l_Qvisc[s]
      end
    end
  end
end

function faceviscterms!(::Val{dim}, ::Val{N}, ::Val{nstate}, ::Val{states_grad},
                        ::Val{ngradstate}, ::Val{nviscstate},
                        ::Val{nauxstate}, viscous_penalty!,
                        viscous_boundary_penalty!, gradient_transform!,
                        Q::Array, Qvisc, auxstate, vgeo, sgeo, t, vmapM, vmapP,
                        elemtobndy, elems) where {dim, N, states_grad,
                                                  ngradstate, nviscstate,
                                                  nstate, nauxstate}
  DFloat = eltype(Q)

  if dim == 1
    Np = (N+1)
    Nfp = 1
    nface = 2
  elseif dim == 2
    Np = (N+1) * (N+1)
    Nfp = (N+1)
    nface = 4
  elseif dim == 3
    Np = (N+1) * (N+1) * (N+1)
    Nfp = (N+1) * (N+1)
    nface = 6
  end

  ngradtransformstate = length(states_grad)

  l_QM = MArray{Tuple{ngradtransformstate}, DFloat}(undef)
  l_auxM = MArray{Tuple{nauxstate}, DFloat}(undef)
  l_GM = MArray{Tuple{ngradstate}, DFloat}(undef)

  l_QP = MArray{Tuple{ngradtransformstate}, DFloat}(undef)
  l_auxP = MArray{Tuple{nauxstate}, DFloat}(undef)
  l_GP = MArray{Tuple{ngradstate}, DFloat}(undef)

  l_Qvisc = MArray{Tuple{nviscstate}, DFloat}(undef)

  @inbounds for e in elems
    for f = 1:nface
      for n = 1:Nfp
        nM = (sgeo[_nx, n, f, e], sgeo[_ny, n, f, e], sgeo[_nz, n, f, e])
        sMJ, vMJI = sgeo[_sMJ, n, f, e], sgeo[_vMJI, n, f, e]
        idM, idP = vmapM[n, f, e], vmapP[n, f, e]

        eM, eP = e, ((idP - 1) ÷ Np) + 1
        vidM, vidP = ((idM - 1) % Np) + 1,  ((idP - 1) % Np) + 1

        # Load minus side data
        for s = 1:ngradtransformstate
          l_QM[s] = Q[vidM, states_grad[s], eM]
        end

        for s = 1:nauxstate
          l_auxM[s] = auxstate[vidM, s, eM]
        end

        gradient_transform!(l_GM, l_QM, l_auxM, t)

        # Load plus side data
        for s = 1:ngradtransformstate
          l_QP[s] = Q[vidP, states_grad[s], eP]
        end

        for s = 1:nauxstate
          l_auxP[s] = auxstate[vidP, s, eP]
        end

        gradient_transform!(l_GP, l_QP, l_auxP, t)

        bctype =
            viscous_boundary_penalty! === nothing ? 0 : elemtobndy[f, e]
        if bctype == 0
          viscous_penalty!(l_Qvisc, nM, l_GM, l_QM, l_auxM, l_GP,
                                  l_QP, l_auxP, t)
        else
          viscous_boundary_penalty!(l_Qvisc, nM, l_GM, l_QM, l_auxM,
                                           l_GP, l_QP, l_auxP, bctype, t)
        end

        for s = 1:nviscstate
          Qvisc[vidM, s, eM] += vMJI * sMJ * l_Qvisc[s]
        end

      end
    end
  end

end

"""
    initauxstate!(::Val{dim}, ::Val{N}, ::Val{nauxstate}, auxstatefun!,
                  auxstate, vgeo, elems) where {dim, N, nauxstate}

Computational kernel: Initialize the auxiliary state

See [`DGBalanceLaw`](@ref) for usage.
"""
function initauxstate!(::Val{dim}, ::Val{N}, ::Val{nauxstate}, auxstatefun!,
                       auxstate, vgeo, elems) where {dim, N, nauxstate}

  # Should only be called in this case I think?
  @assert nauxstate > 0

  DFloat = eltype(auxstate)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  nelem = size(auxstate)[end]

  vgeo = reshape(vgeo, Nq, Nq, Nqk, _nvgeo, nelem)
  auxstate = reshape(auxstate, Nq, Nq, Nqk, nauxstate, nelem)

  l_aux = MArray{Tuple{nauxstate}, DFloat}(undef)

  @inbounds for e in elems
    for k = 1:Nqk, j = 1:Nq, i = 1:Nq
      x, y, z = vgeo[i,j,k,_x,e], vgeo[i,j,k,_y,e], vgeo[i,j,k,_z,e]
      for s = 1:nauxstate
        l_aux[s] = auxstate[i, j, k, s, e]
      end

      auxstatefun!(l_aux, x, y, z)

      for s = 1:nauxstate
        auxstate[i, j, k, s, e] = l_aux[s]
      end
    end
  end
end

"""
    elem_grad_field!(::Val{dim}, ::Val{N}, ::Val{nstate}, Q, vgeo, D, elems, s,
                     sx, sy, sz) where {dim, N, nstate}

Computational kernel: Compute the element gradient of state `s` of `Q` and store
it in `sx`, `sy`, and `sz` of `Q`.

!!! warning

    This does not compute a DG gradient, but only over the element. If ``Q_s``
    is discontinuous you may want to consider another approach.

"""
function elem_grad_field!(::Val{dim}, ::Val{N}, ::Val{nstate}, Q, vgeo,
                          D, elems, s, sx, sy, sz) where {dim, N, nstate}

  DFloat = eltype(vgeo)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  nelem = size(vgeo)[end]

  vgeo = reshape(vgeo, Nq, Nq, Nqk, _nvgeo, nelem)
  Q = reshape(Q, Nq, Nq, Nqk, nstate, nelem)

  s_f = MArray{Tuple{Nq, Nq, Nqk}, DFloat}(undef)
  l_fξ = MArray{Tuple{Nq, Nq, Nqk}, DFloat}(undef)
  l_fη = MArray{Tuple{Nq, Nq, Nqk}, DFloat}(undef)
  l_fζ = MArray{Tuple{Nq, Nq, Nqk}, DFloat}(undef)

  @inbounds for e in elems
    for k = 1:Nqk, j = 1:Nq, i = 1:Nq
      s_f[i,j,k] = Q[i,j,k,s,e]
    end

    # loop of ξ-grid lines
    l_fξ .= 0
    for k = 1:Nqk, j = 1:Nq, i = 1:Nq
      for n = 1:Nq
        l_fξ[i, j, k] += D[i, n] * s_f[n, j, k]
      end
    end
    # loop of η-grid lines
    l_fη .= 0
    for k = 1:Nqk, j = 1:Nq, i = 1:Nq
      for n = 1:Nq
        l_fη[i, j, k] += D[j, n] * s_f[i, n, k]
      end
    end
    # loop of ζ-grid lines
    l_fζ .= 0
    if Nqk > 1
      for k = 1:Nqk, j = 1:Nq, i = 1:Nq
        for n = 1:Nq
          l_fζ[i, j, k] += D[k, n] * s_f[i, j, n]
        end
      end
    end

    for k = 1:Nqk, j = 1:Nq, i = 1:Nq
      ξx, ξy, ξz = vgeo[i,j,k,_ξx,e], vgeo[i,j,k,_ξy,e], vgeo[i,j,k,_ξz,e]
      ηx, ηy, ηz = vgeo[i,j,k,_ηx,e], vgeo[i,j,k,_ηy,e], vgeo[i,j,k,_ηz,e]
      ζx, ζy, ζz = vgeo[i,j,k,_ζx,e], vgeo[i,j,k,_ζy,e], vgeo[i,j,k,_ζz,e]

      Q[i,j,k,sx,e] = ξx * l_fξ[i,j,k] + ηx * l_fη[i,j,k] + ζx * l_fζ[i,j,k]
      Q[i,j,k,sy,e] = ξy * l_fξ[i,j,k] + ηy * l_fη[i,j,k] + ζy * l_fζ[i,j,k]
      Q[i,j,k,sz,e] = ξz * l_fξ[i,j,k] + ηz * l_fη[i,j,k] + ζz * l_fζ[i,j,k]
    end
  end
end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
const _pkg_root = dirname(dirname(@__FILE__))
const _pkg_deps = joinpath(_pkg_root,"deps")
const _pkg_assets = joinpath(_pkg_root,"assets")

!isdir(_pkg_assets) && mkdir(_pkg_assets)

deps = [
    "https://cdnjs.cloudflare.com/ajax/libs/uikit/3.0.0-beta.42/css/uikit.min.css",
    "https://cdnjs.cloudflare.com/ajax/libs/uikit/3.0.0-beta.42/js/uikit.min.js",
    "https://cdnjs.cloudflare.com/ajax/libs/uikit/3.0.0-beta.42/js/uikit-icons.min.js",
    "https://cdn.jsdelivr.net/npm/css-toggle-switch@latest/dist/toggle-switch.css",
]

for dep in deps
    download(dep, joinpath(_pkg_assets, splitdir(dep)[2]))
end


text:
using Pluto

# phases
# 1) create/discover/load causal model  (causal estimand)
# 2) identify statistical model from causal model (statistical estimand)
# 3) estimate effects given data with the statistical model (estimate)
# 4) post-estimate sanity checks (i.e. refute as MS calls it )

# similar projects/sites
# 1) Colliders in Epidemiology: https://watzilei.com/shiny/collider/
# 2) Dagitty: http://www.dagitty.net/dags.html

UI_TYPE = "dash"

if UI_TYPE == "dash"
    println("running dash")
    include("app.jl")
else
    println("running pluto")
    Pluto.run()
end

max_stars_repo_name: Magier/Aetia

max_issues_repo_name: Magier/Aetia

max_forks_repo_name: Magier/Aetia


text:
# defining an empty array
l = []

# pushing values into an empty container
push!(l, 1)             # push a value
push!(l, 2, 3)          # or push a couple values
append!(l, [4, 5, 6])   # or append a container

# defining an undefined array
arr1 = Array{Float64}(undef, 1)
arr2 = Array{Int}(undef, 1)

# defining an undefined array of different sizes
arr3 = Array{Real}(undef, 2, 3)         # 2D
arr4 = Array{Number}(undef, 5, 6, 7)    # 3D

# range sets
r = range(1, 10, length=100)    # step range
u = 1:.005:2                    # unit range

# Float64 of zeros and ones
z = zeros(10)
o = ones(5)

# list comprehension
l1 = [i for i in 1:20]
l2 = [2*i for i in 1:20]

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
abstract type HamiltonianModel end

@doc raw"
    hamiltonian(model<:HamiltonianModel)

return the hamiltonian of the `model` as a two-site tensor operator.
"
function hamiltonian end
struct Ising{T<:Real} <: HamiltonianModel
    β::T
end

const isingβc = log(1+sqrt(2))/2


text:
println("Parsing deployments")
deps = parse_deps(creek)


text:
using QuantLib

const numRows = 5
const numCols = 5
const swaptionVols = [0.1490, 0.1340, 0.1228, 0.1189, 0.1148, 0.1290, 0.1201, 0.1146, 0.1108, 0.1040, 0.1149, 0.1112, 0.1070, 0.1010, 0.0957, 0.1047, 0.1021, 0.0980, 0.0951, 0.1270, 0.1000, 0.0950, 0.0900, 0.1230, 0.1160]
const swaptionLengths = [Dates.Year(1), Dates.Year(2), Dates.Year(3), Dates.Year(4), Dates.Year(5)]

function generate_flatforward_ts{C <: QuantLib.Time.BusinessCalendar}(cal::C, settlementDate::Date)
  flat_rate = Quote(0.04875825)

  ffts = FlatForwardTermStructure(settlementDate, cal, flat_rate, QuantLib.Time.Actual365())

  return ffts
end

function calibrate_model{M <: ShortRateModel}(model::M, helpers::Vector{SwaptionHelper})
  om = QuantLib.Math.LevenbergMarquardt()
  calibrate!(model, helpers, om, QuantLib.Math.EndCriteria(400, 100, 1.0e-8, 1.0e-8, 1.0e-8))

  for i=1:numRows
    j = numCols - (i - 1)
    k = (i - 1) * numCols + j

    npv = model_value!(helpers[i])
    implied = implied_volatility!(helpers[i], npv, 1e-4, 1000, 0.05, 0.50)
    diff = implied - swaptionVols[k]

    println(@sprintf("%i x %i: model %.5f%%, market: %.5f%% (%.5f%%)", i, Int(swaptionLengths[j]), implied * 100, swaptionVols[k] * 100, diff * 100))
  end
end

function main()
  cal = QuantLib.Time.TargetCalendar()
  settlementDate = Date(2002, 2, 19)
  todays_date = Date(2002, 2, 15)
  set_eval_date!(settings, todays_date)

  const swaptionMats = [Dates.Year(1), Dates.Year(2), Dates.Year(3), Dates.Year(4), Dates.Year(5)]

  # flat yield term strucutre implying 1x5 swap at 5%
  rhTermStructure = generate_flatforward_ts(cal, settlementDate)

  # Define the ATM/OTM/ITM swaps
  fixedLegFrequency = QuantLib.Time.Annual()
  fixedLegConvention = QuantLib.Time.Unadjusted()
  floatingLegConvention = QuantLib.Time.ModifiedFollowing()
  fixedLegDayCounter = QuantLib.Time.EuroThirty360()
  floatingLegFrequency = QuantLib.Time.Semiannual()

  swapType = Payer()
  dummyFixedRate = 0.03
  indexSixMonths = euribor_index(QuantLib.Time.TenorPeriod(Dates.Month(6)), rhTermStructure)

  startDate = QuantLib.Time.advance(Dates.Year(1), cal, settlementDate, floatingLegConvention)
  maturity = QuantLib.Time.advance(Dates.Year(5), cal, startDate, floatingLegConvention)

  fixedSchedule = QuantLib.Time.Schedule(startDate, maturity, QuantLib.Time.TenorPeriod(fixedLegFrequency), fixedLegConvention, fixedLegConvention, QuantLib.Time.DateGenerationForwards(), false, cal)
  floatSchedule = QuantLib.Time.Schedule(startDate, maturity, QuantLib.Time.TenorPeriod(floatingLegFrequency), floatingLegConvention, floatingLegConvention, QuantLib.Time.DateGenerationForwards(), false, cal)

  swap = VanillaSwap(swapType, 1000.0, fixedSchedule, dummyFixedRate, fixedLegDayCounter, indexSixMonths, 0.0, floatSchedule, indexSixMonths.dc, DiscountingSwapEngine(rhTermStructure))

  fixedATMRate = fair_rate(swap)
  fixedOTMRate = fixedATMRate * 1.2
  fixedITMRate = fixedATMRate * 0.8

  atmSwap = VanillaSwap(swapType, 1000.0, fixedSchedule, fixedATMRate, fixedLegDayCounter, indexSixMonths, 0.0, floatSchedule, indexSixMonths.dc, DiscountingSwapEngine(rhTermStructure))
  otmSwap = VanillaSwap(swapType, 1000.0, fixedSchedule, fixedOTMRate, fixedLegDayCounter, indexSixMonths, 0.0, floatSchedule, indexSixMonths.dc, DiscountingSwapEngine(rhTermStructure))
  itmSwap = VanillaSwap(swapType, 1000.0, fixedSchedule, fixedITMRate, fixedLegDayCounter, indexSixMonths, 0.0, floatSchedule, indexSixMonths.dc, DiscountingSwapEngine(rhTermStructure))

  times = zeros(0)
  swaptions = Vector{SwaptionHelper}(numRows)

  # First Model
  modelG2 = G2(rhTermStructure)
  hullWhiteModel = HullWhite(rhTermStructure)
  hullWhiteModel2 = HullWhite(rhTermStructure)
  blackKarasinski = BlackKarasinski(rhTermStructure)

  for i = 1:numRows
    j = numCols - (i - 1)
    k = (i - 1) * numCols + j

    sh = SwaptionHelper(swaptionMats[i], swaptionLengths[j], Quote(swaptionVols[k]), indexSixMonths, indexSixMonths.tenor, indexSixMonths.dc, indexSixMonths.dc,
                        rhTermStructure, G2SwaptionEngine(modelG2, 6.0, 16))

    times = add_times_to!(sh, times)
    swaptions[i] = sh
  end

  tg = QuantLib.Time.TimeGrid(times, 30)

  # models
  # modelG2 = G2(rhTermStructure)
  # hullWhiteModel = HullWhite(rhTermStructure)
  # hullWhiteModel2 = HullWhite(rhTermStructure)
  # blackKarasinski = BlackKarasinski(rhTermStructure)
  #
  # for swaptionHelper in swaptions
  #   swaptionHelper.pricingEngine = G2SwaptionEngine(modelG2, 6.0, 16)
  # end

  println("G2 (analytic formulae) calibration")
  calibrate_model(modelG2, swaptions)
  println("calibrated to: ")
  println(@sprintf("a = %.6f, sigma = %.6f", get_params(modelG2)[1], get_params(modelG2)[2]))
  println(@sprintf("b = %.6f, eta = %.6f", get_params(modelG2)[3], get_params(modelG2)[4]))
  println(@sprintf("rho = %.6f", get_params(modelG2)[5]))

  # Hull White
  for i in eachindex(swaptions)
    swaptions[i] = update_pricing_engine(swaptions[i], JamshidianSwaptionEngine(hullWhiteModel))
  end

  println("")
  println("Hull-White (analytic formulae) calibration")
  calibrate_model(hullWhiteModel, swaptions)
  println("calibrated to: ")
  println(@sprintf("a = %.6f, sigma = %.6f", get_params(hullWhiteModel)[1], get_params(hullWhiteModel)[2]))

  # Hull White (numeric)
  for i in eachindex(swaptions)
    swaptions[i] = update_pricing_engine(swaptions[i], TreeSwaptionEngine(hullWhiteModel2, tg))
  end

  println("")
  println("Hull-White (numerical) calibration")
  calibrate_model(hullWhiteModel2, swaptions)
  println("calibrated to: ")
  println(@sprintf("a = %.6f, sigma = %.6f", get_params(hullWhiteModel2)[1], get_params(hullWhiteModel2)[2]))

  # Black Karasinski
  for i in eachindex(swaptions)
    swaptions[i] = update_pricing_engine(swaptions[i], TreeSwaptionEngine(blackKarasinski, tg))
  end

  println("")
  println("Black-Karasinski (numerical) calibration")
  calibrate_model(blackKarasinski, swaptions)
  println("calibrated to: ")
  println(@sprintf("a = %.6f, sigma = %.6f", get_params(blackKarasinski)[1], get_params(blackKarasinski)[2]))

  # ATM Bermudan swaption pricing
  println("")
  println(@sprintf("Payer bermudan swaption struck at %.6f %% (ATM)", fixedATMRate * 100.0))

  swapLeg = swap.legs[1] # Fixed Leg

  bermudanDates = Vector{Date}(length(swapLeg.coupons))
  for i=1:length(swapLeg.coupons)
    bermudanDates[i]  = accrual_start_date(swapLeg.coupons[i])
  end

  bermudanExercise = BermudanExercise(bermudanDates)

  bermudanSwaption = Swaption(atmSwap, bermudanExercise)

  bermudanSwaption = update_pricing_engine(bermudanSwaption, TreeSwaptionEngine(modelG2, 50))

  println(@sprintf("G2 (tree):       %.6f", npv(bermudanSwaption)))

  bermudanSwaption = update_pricing_engine(bermudanSwaption, FdG2SwaptionEngine(modelG2))

  println(@sprintf("G2 (fdm):       %.6f", npv(bermudanSwaption)))

  bermudanSwaption = update_pricing_engine(bermudanSwaption, TreeSwaptionEngine(hullWhiteModel, 50))

  println(@sprintf("HW (tree):       %.6f", npv(bermudanSwaption)))

  bermudanSwaption = update_pricing_engine(bermudanSwaption, FdHullWhiteSwaptionEngine(hullWhiteModel))

  println(@sprintf("HW (fdm):       %.6f", npv(bermudanSwaption)))

  bermudanSwaption = update_pricing_engine(bermudanSwaption, TreeSwaptionEngine(hullWhiteModel2, 50))

  println(@sprintf("HW (num, tree):       %.6f", npv(bermudanSwaption)))

  bermudanSwaption = update_pricing_engine(bermudanSwaption, FdHullWhiteSwaptionEngine(hullWhiteModel2))

  println(@sprintf("HW (num, fdm):       %.6f", npv(bermudanSwaption)))

  bermudanSwaption = update_pricing_engine(bermudanSwaption, TreeSwaptionEngine(blackKarasinski, 50))

  println(@sprintf("BK:       %.6f", npv(bermudanSwaption)))

  # # OTM Bermudan swaption pricing
  # println("")
  # println(@sprintf("Payer bermudan swaption struck at %.6f %% (OTM)", fixedOTMRate * 100.0))
  #
  # otmBermudanSwaption = Swaption(otmSwap, bermudanExercise)
  #
  # update_pricing_engine!(otmBermudanSwaption, TreeSwaptionEngine(modelG2, 300))
  #
  # println(@sprintf("G2 (tree):       %.6f", npv(otmBermudanSwaption)))
  #
  # update_pricing_engine!(otmBermudanSwaption, FdG2SwaptionEngine(modelG2))
  #
  # println(@sprintf("G2 (fdm):       %.6f", npv(otmBermudanSwaption)))
  #
  # update_pricing_engine!(otmBermudanSwaption, TreeSwaptionEngine(hullWhiteModel, 50))
  #
  # println(@sprintf("HW (tree):       %.6f", npv(otmBermudanSwaption)))
  #
  # update_pricing_engine!(otmBermudanSwaption, FdHullWhiteSwaptionEngine(hullWhiteModel))
  #
  # println(@sprintf("HW (fdm):       %.6f", npv(otmBermudanSwaption)))
  #
  # update_pricing_engine!(otmBermudanSwaption, TreeSwaptionEngine(hullWhiteModel2, 50))
  #
  # println(@sprintf("HW (num, tree):       %.6f", npv(otmBermudanSwaption)))
  #
  # update_pricing_engine!(otmBermudanSwaption, FdHullWhiteSwaptionEngine(hullWhiteModel2))
  #
  # println(@sprintf("HW (num, fdm):       %.6f", npv(otmBermudanSwaption)))
  #
  # update_pricing_engine!(otmBermudanSwaption, TreeSwaptionEngine(blackKarasinski, 50))
  #
  # println(@sprintf("BK:       %.6f", npv(otmBermudanSwaption)))
  #
  # # ITM Bermudan swaption pricing
  # println("")
  # println(@sprintf("Payer bermudan swaption struck at %.6f %% (OTM)", fixedITMRate * 100.0))
  #
  # itmBermudanSwaption = Swaption(itmSwap, bermudanExercise)
  #
  # update_pricing_engine!(itmBermudanSwaption, TreeSwaptionEngine(modelG2, 50))
  #
  # println(@sprintf("G2 (tree):       %.6f", npv(itmBermudanSwaption)))
  #
  # update_pricing_engine!(itmBermudanSwaption, FdG2SwaptionEngine(modelG2))
  #
  # println(@sprintf("G2 (fdm):       %.6f", npv(itmBermudanSwaption)))
  #
  # update_pricing_engine!(itmBermudanSwaption, TreeSwaptionEngine(hullWhiteModel, 50))
  #
  # println(@sprintf("HW (tree):       %.6f", npv(itmBermudanSwaption)))
  #
  # update_pricing_engine!(itmBermudanSwaption, FdHullWhiteSwaptionEngine(hullWhiteModel))
  #
  # println(@sprintf("HW (fdm):       %.6f", npv(itmBermudanSwaption)))
  #
  # update_pricing_engine!(itmBermudanSwaption, TreeSwaptionEngine(hullWhiteModel2, 50))
  #
  # println(@sprintf("HW (num, tree):       %.6f", npv(itmBermudanSwaption)))
  #
  # update_pricing_engine!(itmBermudanSwaption, FdHullWhiteSwaptionEngine(hullWhiteModel2))
  #
  # println(@sprintf("HW (num, fdm):       %.6f", npv(itmBermudanSwaption)))
  #
  # update_pricing_engine!(itmBermudanSwaption, TreeSwaptionEngine(blackKarasinski, 50))
  #
  # println(@sprintf("BK:       %.6f", npv(itmBermudanSwaption)))
end