text:
module TestStationary

import LinearAlgebra.SingularException
import IterativeSolvers: DiagonalIndices, FastLowerTriangular,
       FastUpperTriangular, forward_sub!, backward_sub!, OffDiagonal,
       gauss_seidel_multiply!, StrictlyUpperTriangular, StrictlyLowerTriangular

using IterativeSolvers
using Test
using Random
using LinearAlgebra
using SparseArrays

@testset "Stationary solvers" begin

n = 10
m = 6
ω = 1.2
Random.seed!(1234322)

@testset "SparseMatrix{$T, $Ti}" for T in (Float32, Float64, ComplexF32, ComplexF64), Ti in (Int64, Int32)
    @testset "Sparse? $sparse" for sparse = (true, false)
        # Diagonally dominant
        if sparse
            A = sprand(T, n, n, 4 / n) + 2n * I
        else
            A = rand(T, n, n) + 2n * I
        end

        b = rand(T, n)
        x0 = rand(T, n)
        tol = √eps(real(T))

        for solver in (jacobi, gauss_seidel)
            xi = solver(A, b, maxiter=2n)
            @test norm(b - A * xi) / norm(b) ≤ tol
        end

        for solver in (jacobi!, gauss_seidel!)
            xi = solver(copy(x0), A, b, maxiter=2n)
            @test norm(b - A * xi) / norm(b) ≤ tol
        end

        for solver in (sor, ssor)
            xi = solver(A, b, ω, maxiter=2n)
            @test norm(b - A * xi) / norm(b) ≤ tol
        end

        for solver in (sor!, ssor!)
            xi = solver(copy(x0), A, b, ω, maxiter=2n)
            @test norm(b - A * xi) / norm(b) ≤ tol
        end
    end
end

@testset "SOR and Gauss-Seidel" begin
    A = sprand(10, 10, 4/10) + 4I
    x = ones(10)
    b = A * x

    # Gauss-Seidel and SOR should coincide when ω = 1
    gauss_seidel_it = IterativeSolvers.gauss_seidel_iterable(zeros(10), A, b, maxiter = 5)
    sor_it = IterativeSolvers.sor_iterable(zeros(10), A, b, 1.0, maxiter = 5)

    for _ = zip(gauss_seidel_it, sor_it)
        @test gauss_seidel_it.x ≈ sor_it.x
    end
end

@testset "Should throw for singular diagonals" begin
    A = [0.0 1.0; 1.0 0.0]
    B = sparse(A)
    b = rand(2)

    for matrix in (A, B)
        for solver in (jacobi, gauss_seidel)
            @test_throws LinearAlgebra.SingularException solver(matrix, b)
        end

        for solver in (sor, ssor)
            @test_throws LinearAlgebra.SingularException solver(A, b, 0.5)
        end
    end
end

@testset "Diagonal" begin
    A = sprand(10, 10, .3) + 10I
    B = sparse([0. 0. 0.; 0. 1. 0.; 0. 0. 1.])
    D = DiagonalIndices(A)
    @test A.nzval[D.diag] == Diagonal(A).diag
    @test_throws SingularException DiagonalIndices(B)

    x = ones(10)
    y = zeros(10)
    ldiv!(y, D, x)
    @test y ≈ Diagonal(A) \ ones(10)
end

@testset "Forward substitution" begin
    A = sprand(10, 10, .3) + 10I
    D = DiagonalIndices(A)
    L = FastLowerTriangular(A, D)

    x = rand(10)
    y = copy(x)

    forward_sub!(L, x)

    @test x ≈ LowerTriangular(A) \ y
end

@testset "Forward substitution with update" begin
    B = rand(3, 3) + 10I
    A = sparse(B)
    D = DiagonalIndices(A)
    L = FastLowerTriangular(A, D)
    x = rand(3)
    y = ones(3)
    z = copy(x)
    α = 2.0
    β = 3.0

    forward_sub!(α, L, x, β, y)

    z[1] = β * y[1] + α * z[1] / B[1,1]
    z[2] = β * y[2] + α * (z[2] - B[2,1] * z[1]) / B[2, 2]
    z[3] = β * y[3] + α * (z[3] - B[3,1] * z[1] - B[3,2] * z[2]) / B[3, 3]

    @test x ≈ z
end

@testset "Backward substitution" begin
    A = sprand(10, 10, .3) + 10I
    D = DiagonalIndices(A)
    L = FastUpperTriangular(A, D)

    x = rand(10)
    y = copy(x)

    backward_sub!(L, x)

    @test x ≈ UpperTriangular(A) \ y
end

@testset "Backward substitution with update" begin
    B = rand(3, 3) + 10I
    A = sparse(B)
    D = DiagonalIndices(A)
    L = FastUpperTriangular(A, D)
    x = rand(3)
    y = ones(3)
    z = copy(x)
    α = 2.0
    β = 3.0

    backward_sub!(α, L, x, β, y)

    z[3] = β * y[3] + α * z[3] / B[3,3]
    z[2] = β * y[2] + α * (z[2] - B[2,3] * z[3]) / B[2, 2]
    z[1] = β * y[1] + α * (z[1] - B[1,2] * z[2] - B[1,3] * z[3]) / B[1, 1]

    @test x ≈ z
end

@testset "mul! with OffDiagonal type" begin
    A = sprand(10, 10, .3) + 2I
    D = DiagonalIndices(A)
    O = OffDiagonal(A, D)

    x = rand(10)
    y = ones(10)
    mul!(1.0, O, x, 0.0, y)
    @test y ≈ A * x - Diagonal(A) * x

    x = rand(10)
    y = ones(10)
    mul!(1.0, O, x, 1.0, y)
    @test y ≈ A * x - Diagonal(A) * x + ones(10)

    x = rand(10)
    y = ones(10)
    mul!(2.0, O, x, 3.0, y)
    @test y ≈ 2 * (A * x - Diagonal(A) * x) + 3 * ones(10)
end

@testset "Gauss-Seidel multiplication upper" begin
    A = sprand(10, 10, .3) + 2I
    D = DiagonalIndices(A)
    U = StrictlyUpperTriangular(A, D)

    x = rand(10)
    y = copy(x)
    b = rand(10)

    # x ← b - U * x
    gauss_seidel_multiply!(-1.0, U, x, 1.0, b, x)

    @test x ≈ b - triu(A, 1) * y
end

@testset "Gauss-Seidel multiplication lower" begin
    A = sprand(10, 10, .3) + 2I
    D = DiagonalIndices(A)
    L = StrictlyLowerTriangular(A, D)

    x = rand(10)
    y = copy(x)
    b = rand(10)

    # x ← b - L * x
    gauss_seidel_multiply!(-1.0, L, x, 1.0, b, x)

    @test x ≈ b - tril(A, -1) * y
end
end

end # module


text:
import Base.read, Base.eof, Base.close

# the conversion never returns 0x00 as a continuation character
const EMPTY = 0x00

type Latin1Buffer <: IO
    io::IO
    padding::Uint8

    Latin1Buffer(io::IO) = new(io, EMPTY)
end

function Base.read(l1::Latin1Buffer, ::Type{Uint8})
    if l1.padding == EMPTY
        c = read(l1.io, Uint8)
        if c >= 0x80
            c, l1.padding = c < 0xc0 ? (0xc2, c) : (0xc3, uint8(c - 0x40))
        end
        c
    else
        tmp = l1.padding
        l1.padding = EMPTY
        tmp
    end
end

Base.eof(l1::Latin1Buffer) = l1.padding == EMPTY && eof(l1.io)

Base.close(l1::Latin1Buffer) = close(l1.io)

ext: jl


text:
"""This notebook is space delimited, as you can see.

In the newer versions of this script, a flag is not necessary to have it parsed that way.

Just don't include any code or md tags!"""

a = rand(5,5)

	b = a + 5 # This is tabbed in the script (but not in the notebook!)

"""This is markdown squished on top of a code block"""
    c = a+250

"""#Markdown header!"""


text:
export CartesianTuple

struct CartesianTuple{N, T<:Number, U<:ComponentType}
    coordinates::SVector{N,T}

    function CartesianTuple{N,T,U}(args::SVector{N,T}) where {N,T,U}
        @assert !any(isnan.(args))
        new{N,T,U}(args)
    end

    function CartesianTuple{N,T,U}(args::Vararg{T,N}) where {N,T,U}
        @assert !any(isnan.(args))
        new{N,T,U}(SVector(args))
    end

    function CartesianTuple{N,T,U}(arg::T) where {N,T,U}
        @assert !isnan(arg)
        new{N,T,U}(@SVector fill(arg, N))
    end

    CartesianTuple{N,T,U}() where {N,T,U} = new{N,T,U}(zeros(SVector{N,T}))
end

CartesianTuple{N,T,U}(args::Vararg{Any,N}) where {N,T,U} = CartesianTuple{N,T,U}(convert.(T, args)...)
CartesianTuple{N,T,U}(arg) where {N,T,U} = CartesianTuple{N,T,U}(convert(T, arg))
CartesianTuple{N,T,U}(args::AbstractArray{T,N}) where {N,T,U} = CartesianTuple{N,T,U}(args...)


text:
export all_fuel, all_fuel_recursive

fuel(mass) = mass ÷ 3 - 2
all_fuel(modules) = sum(fuel, modules)

function fuel_recursive(mass)
    next_fuel = fuel(mass)
    return next_fuel < 1 ? 0 : next_fuel + fuel_recursive(next_fuel)
end

all_fuel_recursive(modules) = sum(fuel_recursive, modules)


text:
module EE263

# Write your package code here.

export name
export greeting

function greeting()
    println("Hello, World!")
end

function name()
    println("Larry")
end

end


text:
gensym_unhashed(s::Symbol) = Symbol(Base.replace(string(gensym(s)), "#"=>""))

constraint_walk(f, expr) = postwalk(expr) do x
    @capture(x, p_ == 0) ? f(p) : x
end

function_walk(f, expr) = postwalk(expr) do x
    @capture(x, g_(a__)) && issymbol(g) ? f(g, a) : x
end

symbol_walk(f, ex) = postwalk(x -> issymbol(x) ? f(x) : x, ex)

atom_walk(f, x) = walk(x, x -> (@capture(x, y_(ys__)) && issymbol(y)) ? f(x) : atom_walk(f, x) , f)

issymbol(x) = x isa Symbol && Base.isidentifier(x)
isfunction(x) = @capture(x, s_(xs__)) && s isa Symbol && Base.isidentifier(s)

function symbols(x::Expr)
    ls = Symbol[]
    symbol_walk(x) do s
        push!(ls, s)
        s
    end
    unique(ls)
end
symbols(x::Symbol) = [x]
symbols(::Number) = []


text:
"""Convenience function to add an id attribute to a html element."""
attr(name::Symbol, val::AS) = ifelse(isempty(val), "", " $name=\"$val\"")

"""Convenience function for a sup item."""
html_sup(id::String, in::AS) =  "$in"

"""Convenience function for a header."""
html_hk(hk::String, t::AS; id::String="") = "<$hk$(attr(:id, id))>$t</$hk>"

"""Convenience function to introduce a hyper reference."""
function html_ahref(link::AS, name::Union{Int,AS};
                    title::AS="", class::AS="")
    a  = "$name"
    return a
end

"""Convenience function to introduce a hyper reference relative to a key."""
html_ahref_key(k::AS, n::Union{Int,AS}) = html_ahref("#$k", n)

"""Convenience function to introduce a div block."""
html_div(name::AS, in::AS) = "
$in
"

"""Convenience function to introduce a named span block."""
html_span(name::AS, in::AS) = "$in"

"""
Convenience function to introduce an image. The alt is escaped just in case the
user adds quotation marks in the alt string.
"""
html_img(src::AS, alt::AS="") = ""

"""
Convenience function to introduce a code block. With indented blocks, the
language will be resolved at build time; as a result no lines will be hidden
for such blocks.
"""
function html_code(c::AS, lang::AS="")::String
    isempty(c)    && return ""
    # if it's plaintext
    isempty(lang) && return "
$c
"
    # if it's html code escape it if it isn't
    lang == "html" && !is_html_escaped(c) && (c = htmlesc(c))
    # remove hidden lines if any
    c = html_skip_hidden(c, lang)
    # if empty (e.g. via #hideall) return ""
    c == "" && return ""
    return "
$c
"
end

const REGEX_CODE_HIDE = Regex(raw"(?:^|[^\S\r\n]*?)#(\s)*?(?i)hide(all)?")
const REGEX_LIT_HIDE  = Regex(raw"(?:^|[^\S\r\n]*?)#src")

"""Convenience function to process a code string and hide some lines."""
function html_skip_hidden(c::AS, lang::AS)::String
    # if the language is not one of CODE_LANG, just return
    # without hiding anything
    lang in keys(CODE_LANG) || return c
    # otherwise retrive the symbol marking a comment
    _, comsym = CODE_LANG[lang]
    # read code line by line and write to buffer
    buf = IOBuffer()
    for line in split(c, '\n')
        m  = match(REGEX_CODE_HIDE, line)
        ml = match(REGEX_LIT_HIDE, line)
        if m === ml === nothing
            println(buf, line)
        elseif m !== nothing
            # does it have a "all" or not?
            isnothing(m.captures[2]) && continue
            # if it doesn return an empty string
            return ""
        end
        # in other cases: skip the line
    end
    # strip as there may be a stray `\n`
    return strip(String(take!(buf)))
end

"""Convenience function to introduce inline code."""
html_code_inline(c::AS) = "$c"

"""Insertion of a visible red message in HTML to show there was a problem."""
html_err(mess::String="") =
    "
// $mess //
"

"""Helper function to get the relative url of the current page."""
function url_curpage()
    if FD_ENV[:STRUCTURE] < v"0.2"
        return _url_curpage()
    end
    return _url_curpage2()
end

function _url_curpage()
    # go from /pages/.../something.md to /pub/.../something.html note that if
    # on windows then it would be \\ whence the PATH_SEP
    rp = replace(locvar("fd_rpath"),
            Regex("^pages$(escape_string(PATH_SEP))") => "pub$(PATH_SEP)")
    rp = unixify(rp)
    rp = splitext(rp)[1] * ".html"
    startswith(rp, "/") || (rp = "/" * rp)
    return rp
end

function _url_curpage2()
    # get the relative path to current page and split extension (.md)
    fn = splitext(locvar("fd_rpath"))[1]
    # if it's not `index` then add `index`:
    if splitdir(fn)[2] != "index"
        fn = joinpath(fn, "index")
    end
    fn *= ".html"
    # unixify
    fn = unixify(fn)
    # if it does not start with "/", add a "/" in front
    startswith(fn, "/") || (fn = "/" * fn)
    # append .html and return
    return fn
end

# Copied from https://github.com/JuliaLang/julia/blob/acb7bd93fb2d5adbbabeaed9f39ab3c85495b02f/stdlib/Markdown/src/render/html.jl#L25-L31
const _htmlescape_chars = LittleDict(
            '<' => "<",
            '>' => ">",
            '"' => """,
            '&' => "&"
            )
for ch in "'`!\$%()=+{}[]"
    _htmlescape_chars[ch] = "&#$(Int(ch));"
end

const _htmlesc_to = values(_htmlescape_chars) |> collect

"""Internal function to check if some html code has been escaped."""
is_html_escaped(cs::AS) =
    !isnothing(findfirst(ss -> occursin(ss, cs), _htmlesc_to))

"""
Internal function to reverse the escaping of some html code (in order to avoid
double escaping when pre-rendering with highlight, see issue 326).
"""
function html_unescape(cs::String)
    # this is a bit inefficient but whatever, `cs` shouldn't  be very long.
    for (ssfrom, ssto) in _htmlescape_chars
        cs = replace(cs, ssto => ssfrom)
    end
    return cs
end


text:
#!/usr/bin/env julia

using LaTeXStrings
import PyPlot.plt
using NPZ

img = npzread("img.npy")
tau = npzread("tau.npy")
mask = npzread("mask.npy")
mask_height = npzread("mask_height.npy")

tau8 = tau .> 8.0

nx, ny, nv = size(img)

# Estimate how many pixels show emission above some threshold
max_pix = maximum(img)
println("Max pix value ", max_pix)
above = (img .> (1e-5 * max_pix))
println("Total pixels traced ", sum(above))
println("Fraction of pixes unmasked ", sum(mask)/(nx * ny * nv))
println("Fraction of pixes unmasked-height ", sum(mask_height)/(nx * ny * nv))
println("Fraction of pixels brighter than max * 1e-3, ", sum(above)/(nx * ny * nv))

println("Total pixels mask ", sum(mask))
println("Total pixels mask height ", sum(mask_height))

fig, ax = plt[:subplots](nrows=5, ncols=nv, figsize=(1.5 * nv, 8.0))
for k=1:nv
    ax[1,k][:imshow](img[:,:,k], origin="lower", interpolation="none")
    ax[2,k][:imshow](tau[:,:,k], origin="lower", interpolation="none")
    ax[3,k][:imshow](tau8[:,:,k], origin="lower", interpolation="none")
    ax[4,k][:imshow](mask[:,:,k], origin="lower", interpolation="none")
    ax[5,k][:imshow](mask_height[:,:,k], origin="lower", interpolation="none")
end
fig[:subplots_adjust](left=0.03, right=0.97, wspace=0.01)
plt[:savefig]("chmaps.png")


text:
const _DEFAULT_ENGINE_KEY = "texengine"
const _DEFAULT_PREAMBLE_KEY = "preamble"

function _set_default_texengine(::T; export_prefs::Bool=false) where T <: LaTeXEngine
    set_preferences!(LatexSVG, _DEFAULT_ENGINE_KEY => string(nameof(T)); export_prefs=export_prefs)
    @info "The default LaTeX engine is set to $(string(nameof(T))). This will take effect after you restart your Julia session."
    return nothing
end

function _set_default_preamble(pre::Vector{<:AbstractString}; export_prefs::Bool=false)
    set_preferences!(LatexSVG, _DEFAULT_PREAMBLE_KEY => String.(pre); export_prefs=export_prefs)
    info = "The default preamble is set to:\n\n  " * join((pre), "\n  ") * "\n\nThis will take effect after you restart your Julia session."
    @info info
    return nothing
end
_set_default_preamble(pre::AbstractString; export_prefs::Bool=false) = _set_default_preamble(split(strip(pre), "\n"); export_prefs=export_prefs)

"""
    config(key::String)

Obtain the persistent settings for the LaTeX engine and the preamble. `key` can be `"texengine"` or `"preamble"`. Use [`config!`](@ref) to configure these settings.
"""
function config(key::String; info=true)
    if @has_preference(key)
        return @load_preference(key)
    else
        if info
            @info "Preference \"$(key)\" is unset."
        end
        return nothing
    end
end

"""
    config!(; texengine=nothing, preamble=nothing, export_prefs=false)

Configure persistent settings for the LaTeX engine and the preamble:
- `texengine` can be [`XeLaTeX`](@ref), [`PDFLaTeX`](@ref), or [`LuaLaTeX`](@ref)
- `preamble` can be an `AbstractString` or a `Vector` of `AbstractString`s.

Two examples:
```julia-repl
julia> config!(texengine=PDFLaTeX)
```
This sets the default LaTeX engine to be `PDFLaTeX` *for all future sessions*.

```julia-repl
julia> config!(preamble=["\\usepackage{mathtools}", "\\usepackage{xcolor}"])
```
This sets the default preamble to be
```latex
\\usepackage{mathtools}
\\usepackage{xcolor}
```
*for all future sessions*, replacing the original default of
```latex
\\usepackage{amsmath,amsthm,amssymb}
\\usepackage{color}
```

You can also call `config!()` without any keyword argument. In this case your current LaTeX engine and preamble will be stored as the default.

These preferences are stored in a `LocalPreferences.toml` file in your active project. If `export_prefs=true`, they will instead be written into your `Project.toml`. Either way, if you have multiple projects using this package you need to set the preference for them individually.
"""
function config!(;
    texengine::Union{Nothing, T, Type{T}}=nothing,
    preamble::Union{Nothing, <:AbstractString, Vector{<:AbstractString}}=nothing,
    export_prefs::Bool=false
) where T <: LaTeXEngine
    if texengine !== nothing
        _set_default_texengine(T(); export_prefs=export_prefs)
    end
    if preamble !== nothing
        _set_default_preamble(preamble; export_prefs=export_prefs)
    end
    if texengine === nothing && preamble === nothing
        _set_default_texengine(LatexSVG.texengine(); export_prefs=export_prefs)
        _set_default_preamble(LatexSVG.preamble(); export_prefs=export_prefs)
    end
    return nothing
end


text:
combinations_upto(x, n) = Iterators.flatten(combinations(x, i) for i in 1:n)

"""
    (term1, term2, ...) ^ n

Generate all interactions of terms up to order ``n``.

!!! warning
    If any term is an `InteractionTerm`, then nonsensical interactions may
    arise, e.g. `a & a & b`.
"""
function Base.:(^)(args::TermTuple, deg::ConstantTerm{<:Integer})
    deg.n > 0 || throw(ArgumentError("power should be greater than zero (got $deg)"))
    tuple(((&)(terms...) for terms in combinations_upto(args, deg.n))...)
end

function Base.:(^)(::TermTuple, deg::AbstractTerm)
    throw(ArgumentError("power should be an integer constant (got $deg)"))
end

function StatsModels.apply_schema(
    t::FunctionTerm{typeof(^)},
    sch::StatsModels.FullRank,
    ctx::Type{<:RegressionModel}
)
    length(t.args_parsed) == 2 ||
        throw(ArgumentError("invalid term $t: should have exactly two arguments"))
    first, second = t.args_parsed
    second isa ConstantTerm{<:Integer} ||
        throw(ArgumentError("invalid term $t: power should be an integer (got $second)"))
    apply_schema.(first^second, Ref(sch), ctx)
end


text:
#-----------------------------------------------------------

msg("Test fit SVM as a Statsmodel on DataFrame")

using RDatasets
iris = dataset("datasets", "iris")
iris = iris[iris[:Species] .!= "setosa", :]

lm=labelmap(iris[:Species])
length(lm.vs) == 2
(labelencode(lm, iris[:Species]) - 1.5) * 2


text:
returnAssignmentName = :result

# Return setup functions

function findCustomReturnName(expression)
    returnNameExpr = getAgreementsFromExpressions(:returnName, :(=), expression)
    if length(returnNameExpr) == 0
        return returnAssignmentName
    elseif length(returnNameExpr) != 1 || typeof(returnNameExpr[1]) != Symbol
        throw(ArgumentError(
"""
Custom Return names should be added in the format 'resultName = \$customName'
"""
                            ))
    else
        return returnNameExpr[1]
    end
end

# Expression list setup by Agreement

function getAgreementsFromExpressions(agreementType::Symbol, head::Symbol,
                                      expressions)
    conditions = []
    for expr in expressions
        if typeof(expr) == Expr &&
            expr.head == head && expr.args[1] == agreementType
            append!(conditions, expr.args[2:end])
        end
    end
    return conditions
end

function lookAndFillType(agreement::Agreement,
                         topExpr)
    typeSym = Symbol(agreement.type)
    conditions = getAgreementsFromExpressions(typeSym, :call, topExpr.args)
    append!(agreement.expressions, conditions)
end

# Agreement setups on function Body

function addRequirements!(functionBody::Expr, agreement::Agreement)
    start = createCheckExpressions(agreement)
    pushfirst!(functionBody.args, start...)
end

function addEnsures!(functionBody::Expr, agreement::Agreement,
                     returnName::Symbol)
    finish = createCheckExpressions(agreement)
    push!(finish, :(return $returnName))
    # Look for "return" keywords and replace them
    newFunctionBody = MacroTools.postwalk(functionBody) do x
        @capture(x, return a_) || return x
        newReturn = Expr(:block)
        newReturn.args = vcat([:($returnName = $a)], finish)
        return newReturn
    end
    functionBody.args = newFunctionBody.args
    functionBody.head = newFunctionBody.head
end

function generateEnsureFunction(returnName::Symbol)
    return (x, y) -> addEnsures!(x, y, returnName)
end

# Finding function

function seekFunctionDefinition(expressions)
    for expr ∈ expressions
        if typeof(expr) == Expr && expr.head == :function
            return expr
        end
    end
end

"""
@contract

Take in a :block expression and define a function wraped inside its contract. A
contract is defined by "require" expressions and "ensure" expressions. It will
constrain function calls and returns to have specific user defined conditions
(It'll wrap the code with "if" statements).

```julia-repl

julia> @contract begin
           require(a>=-1)
           ensure(result>0)
           function foo(a)
               return a+1
           end
       end
foo (generic function with 1 method)

julia> foo(-2)
ERROR: ContractBreachException: Breach on Requirement Expression 'a >= -1' in function 'foo'
...

julia> foo(-1)
ERROR: ContractBreachException: Breach on Ensure Expression 'result > 0' in function 'foo'
...

julia> foo(0)
1
```
"""
macro contract(expr)
    @assert expr.head == :block
    functionExpr = seekFunctionDefinition(expr.args)
    @assert typeof(functionExpr) == Expr """
Check the arguments previous to the function declaration!
"""
    functionName = string(functionExpr.args[1].args[1])
    functionBody = functionExpr.args[2]

    returnName = findCustomReturnName(expr.args)

    agreements = [
        Agreement(require, addRequirements!, functionName),
        Agreement(ensure, generateEnsureFunction(returnName), functionName),
    ]

    for agreement in values(agreements)
        lookAndFillType(agreement, expr)
        agreement(functionBody)
    end

    return functionExpr |> esc
end


text:
# ----------- Empty constrained selection ------------ #

struct ConstrainedEmptySelection <: ConstrainedSelection end

has_top(cas::ConstrainedEmptySelection, addr) = false
dump_queries(cas::ConstrainedEmptySelection) = Set()
get_top(cas::ConstrainedEmptySelection, addr) = error("ConstrainedEmptySelection has no queries!")
get_sub(cas::ConstrainedEmptySelection, addr) = cas
isempty(cas::ConstrainedEmptySelection) = true
function merge(cl::T, sel::ConstrainedEmptySelection) where T <: CallSite
    cl_selection = get_selection(cl)
    return cl_selection
end

# Pretty printing.
function Base.display(chs::ConstrainedEmptySelection; show_values = false)
    println("  __________________________________\n")
    println("            Empty Selection\n")
    println("  __________________________________\n")
end

# ----------- Empty unconstrained selection ------------ #

struct UnconstrainedEmptySelection <: UnconstrainedSelection end

has_top(uas::UnconstrainedEmptySelection, addr) = false
dump_queries(uas::UnconstrainedEmptySelection) = Set()
get_sub(cas::UnconstrainedEmptySelection, addr) = cas
isempty(cas::UnconstrainedEmptySelection) = true


text:
export BinomialHomotopy

"""
    BinomialHomotopy

Construct the homotopy ``x^A - t γ + (1-t) b``
where each column in A denotes a monomial.
"""
struct BinomialHomotopy <: AbstractHomotopy
    A::Matrix{Int}
    b::Vector{ComplexF64}
    γ::Vector{Float64}

end

function BinomialHomotopy(A::Matrix{Int}, b::Vector{ComplexF64})
    γ = sign.(real.(b))
    BinomialHomotopy(A, b, γ)
end

cache(H::BinomialHomotopy, x, t) = HomotopyNullCache()
Base.size(H::BinomialHomotopy) = size(H.A)

function evaluate!(u, H::BinomialHomotopy, x, t, c::HomotopyNullCache)
    @unpack A, b, γ = H
    n = length(b)
    @inbounds for i = 1:n
        xAᵢ = one(x[1])
        for j = 1:n
            if A[j, i] != 0
                xAᵢ *= x[j]^A[j, i]
            end
        end
        u[i] = xAᵢ - (t * γ[i] + (1.0 - t) * b[i])
    end
    u
end

function dt!(u, H::BinomialHomotopy, x, t, c::HomotopyNullCache)
    @unpack b, γ = H
    n = length(b)
    @inbounds for i = 1:n
        u[i] = b[i] - γ[i]
    end
    u
end

function jacobian!(U, H::BinomialHomotopy, x, t, c::HomotopyNullCache)
    @unpack A, b, γ = H
    n = length(b)
    @inbounds for i = 1:n
        xAᵢ = one(x[1])
        for j = 1:n
            if A[j, i] != 0
                xAᵢ *= x[j]^A[j, i]
            end
        end
        for j = 1:n
            if A[j, i] == 0
                U[i, j] = 0
            else
                U[i, j] = A[j, i] * (@fastmath xAᵢ / x[j])
            end
        end
    end
    U
end

function evaluate_and_jacobian!(u, U, H::BinomialHomotopy, x, t, c::HomotopyNullCache)
    @unpack A, b, γ = H
    n = length(b)
    @inbounds for i = 1:n
        xAᵢ = one(x[1])
        for j = 1:n
            if A[j, i] != 0
                xAᵢ *= x[j]^A[j, i]
            end
        end
        u[i] = xAᵢ - (t * γ[i] + (1.0 - t) * b[i])

        for j = 1:n
            if A[j, i] == 0
                U[i, j] = 0
            else
                U[i, j] = A[j, i] * (@fastmath xAᵢ / x[j])
            end
        end
    end
    nothing
end

function evaluate(H::BinomialHomotopy, x, t, c::HomotopyNullCache)
    T = promote_type(eltype(x), ComplexF64)
    evaluate!(zeros(T, size(H.A, 1)), H, x, t, c)
end
(H::BinomialHomotopy)(x, t, c = cache(H, x, t)) = evaluate(H, x, t, c)

function jacobian(H::BinomialHomotopy, x, t, c::HomotopyNullCache)
    T = promote_type(eltype(x), ComplexF64)
    jacobian!(zeros(T, size(H.A)), H, x, t, c)
end
function dt(H::BinomialHomotopy, x, t, c::HomotopyNullCache)
    dt!(zeros(ComplexF64, size(H.A, 1)), H, x, t, c)
end


text:
struct MPSDescription{T<:AbstractFloat,M<:AbstractMatrix}
  Q::M
  q₁::Vector{T}
  q₂::T
  c₁::Vector{T}
  C::M
  c₂::Vector{T}
  A::M
  b::Vector{T}
  lb::Vector{T}
  ub::Vector{T}
  vars::Vector{Symbol}
  name::String
end

function (::Type{MPSDescription})(::Type{T}, ::Type{Matrix}, n::Int, m₁::Int,
  m₂::Int, name::AbstractString = "QP") where {T<:AbstractFloat}
  MPSDescription{T,Matrix{T}}(zeros(T, n, n), zeros(T, n), zero(T), fill(convert(T, -Inf), m₁),
    zeros(T, m₁, n), fill(convert(T, Inf), m₁), zeros(T, m₂, n), zeros(T, m₂),
    zeros(T, n), fill(convert(T, Inf), n), Symbol[Symbol(:x_,k) for k in 1:n],
    name)
end

function (::Type{MPSDescription})(::Type{T}, ::Type{SparseMatrixCSC}, n::Int, m₁::Int,
  m₂::Int, name::AbstractString = "QP") where {T<:AbstractFloat}
  MPSDescription{T,SparseMatrixCSC{T,Int}}(spzeros(T, n, n), zeros(T, n), zero(T), fill(convert(T, -Inf), m₁),
    spzeros(T, m₁, n), fill(convert(T, Inf), m₁), spzeros(T, m₂, n), zeros(T, m₂),
    zeros(T, n), fill(convert(T, Inf), n), Symbol[Symbol(:x_,k) for k in 1:n],
    name)
end

MPSDescription(n::Int, m₁::Int, m₂::Int, name::AbstractString = "QP")     =
  MPSDescription(Float64, Matrix, n, m₁, m₂, name)
MPSDescription(t::Type{T}, name::AbstractString = "QP") where {T<:AbstractFloat} =
  MPSDescription(t, 0, 0, 0, name)
MPSDescription(name::AbstractString = "QP")                               =
  MPSDescription(Float64, name)

SparseMPSDescription(n::Int, m₁::Int, m₂::Int, name::AbstractString = "QP")     =
  MPSDescription(Float64, SparseMatrixCSC, n, m₁, m₂, name)
SparseMPSDescription(t::Type{T}, name::AbstractString = "QP") where {T<:AbstractFloat} =
  SparseMPSDescription(t, 0, 0, 0, name)
SparseMPSDescription(name::AbstractString = "QP")                               =
  SparseMPSDescription(Float64, name)

function MPSDescription(::Type{T}, Q::AbstractMatrix, q₁::AbstractVector, q₂::Real,
  c₁::AbstractVector, C::AbstractMatrix, c₂::AbstractVector, A::AbstractMatrix, b::AbstractVector,
  lb::AbstractVector, ub::AbstractVector, vars::AbstractVector{S},
  name::AbstractString) where {T<:AbstractFloat,S<:Union{Symbol,Char,AbstractString}}

  Q, C, A = convert.(Matrix{T}, (Q, C, A))

  _mpscheck(Q, q₁, q₂, c₁, C, c₂, A, b, lb, ub, vars)

  MPSDescription{T,Matrix{T}}(Q, q₁, q₂, c₁, C, c₂, A, b, lb, ub, Symbol.(vars), name)
end

function MPSDescription(::Type{T}, Q::AbstractSparseMatrix, q₁::AbstractVector, q₂::Real,
  c₁::AbstractVector, C::AbstractSparseMatrix, c₂::AbstractVector, A::AbstractSparseMatrix, b::AbstractVector,
  lb::AbstractVector, ub::AbstractVector, vars::AbstractVector{S},
  name::AbstractString) where {T<:AbstractFloat,S<:Union{Symbol,Char,AbstractString}}

  Q, C, A = convert.(SparseMatrixCSC{T,Int}, (Q, C, A))

  _mpscheck(Q, q₁, q₂, c₁, C, c₂, A, b, lb, ub, vars)

  MPSDescription{T,SparseMatrixCSC{T,Int}}(Q, q₁, q₂, c₁, C, c₂, A, b, lb, ub, Symbol.(vars), name)
end

function _mpscheck(Q::AbstractMatrix, q₁::AbstractVector, q₂::Real, c₁::AbstractVector,
  C::AbstractMatrix, c₂::AbstractVector, A::AbstractMatrix, b::AbstractVector,
  lb::AbstractVector, ub::AbstractVector, vars::AbstractVector{S}) where {S<:Union{Symbol,Char,AbstractString}}

  n₁, n₂  = size(Q)
  m₁, n₃  = size(C)
  m₂, n₄  = size(A)
  n₅      = length(vars)

  if n₁ ≠ n₂
    throw(DomainError(Q ,"MPSDescription: `Q` matrix must be symmetric"))
  elseif n₁ ≠ length(q₁)
    throw(DomainError(q₁ ,"MPSDescription: `q₁` must have $(n₁) elements"))
  elseif n₁ ≠ n₃
    throw(DomainError(C ,"MPSDescription: `C` must have $(n₁) columns"))
  elseif m₁ ≠ length(c₁) || m₁ ≠ length(c₂)
    throw(DomainError((c₁, c₂) ,"MPSDescription: Both `c₁` and `c₂` must have $(m₁) elements"))
  elseif n₁ ≠ n₄
    throw(DomainError(A ,"MPSDescription: `A` must have $(n₁) columns"))
  elseif m₂ ≠ length(b)
    throw(DomainError(b ,"MPSDescription: `b` must have $(m₂) elements"))
  elseif n₁ ≠ length(lb) || n₁ ≠ length(ub)
    throw(DomainError((lb, ub) ,"MPSDescription: Both `lb` and `ub` must have $(n₁) elements"))
  elseif n₁ ≠ length(vars)
    throw(DomainError(vars ,"MPSDescription: `vars` must have $(n₁) elements"))
  end
end

MPSDescription(Q::AbstractMatrix, q₁::AbstractVector, q₂::Real,
  c₁::AbstractVector, C::AbstractMatrix, c₂::AbstractVector,
  A::AbstractMatrix, b::AbstractVector, lb::AbstractVector,
  ub::AbstractVector, vars::AbstractVector{S}, name::AbstractString) where {S<:Union{Symbol,Char,AbstractString}} =
  MPSDescription(Float64, Q, q₁, q₂, c₁, C, c₂, A, b, lb, ub, vars, name)

"""
    name(qp) -> qpname

Return the problem instance `qp`'s name. `qp` can be either `MPSDescription` or
`CanonicalDescription`. `qpname` is a `String`.

See also: `variables`, `objective`, `equalities`, `inequalities`, `bounds`,
  `parseqps` and `canonical`.
"""
name(qp::MPSDescription)          = qp.name

"""
    variables(qp) -> vars

Return the problem instance `qp`'s variable names. `qp` can be either
`MPSDescription` or `CanonicalDescription`. `vars` is a `Vector{Symbol}`.

See also: `name`, `objective`, `equalities`, `inequalities`, `bounds`, `parseqps`
  and `canonical`.
"""
variables(qp::MPSDescription)     = qp.vars

"""
    objective(qp) -> (Q, q₁, q₂)

Return the problem instance `qp`'s objective part. `qp` can be either
`MPSDescription` or `CanonicalDescription`. `Q` is a `Matrix{T}`, `q₁` is a
`Vector{T}` and `q₂` is a `T`, where `T` is the first argument in `parseqps`.

See also: `name`, `variables`, `equalities`, `inequalities`, `bounds`, `parseqps`
  and `canonical`.
"""
objective(qp::MPSDescription)     = (qp.Q, qp.q₁, qp.q₂ )

"""
    equalities(qp) -> (A, b)

Return the problem instance `qp`'s `A, b` pair that defines the equalities among
the problem's variables. `qp` can be either `MPSDescription` or `CanonicalDescription`.
`A` is a `Matrix{T}` and `b` is a `Vector{T}`, where `T` is the first argument
in `parseqps`.

See also: `name`, `variables`, `objective`, `inequalities`, `bounds`, `parseqps`
  and `canonical`.
"""
equalities(qp::MPSDescription)    = (qp.A, qp.b         )

"""
    inequalities(qp::MPSDescription)        -> (C, c₁, c₂)
    inequalities(qp::CanonicalDescription)  -> (C̃, c̃)

Return the problem instance `qp`'s part that defines the inequalities among
the problem's variables.

When `qp` is `MPSDescription`, `C` is a `Matrix{T}`, and `c₁` and `c₂` are
`Vector{T}`, where `T` is the first argument in `parseqps`. When `qp` is
`CanonicalDescription`, `C̃` is a `Matrix{T}` and `c̃` is a `Vector{T}`.

See also: `name`, `variables`, `objective`, `equalities`, `bounds`, `parseqps`
and `canonical`.
"""
inequalities(qp::MPSDescription)  = (qp.C, qp.c₁, qp.c₂ )

"""
    bounds(qp) -> (lb, ub)

Return the lower and upper bounds, `lb` and `ub`, defined on the variables of the
problem instance `qp`. `qp` can be either `MPSDescription` or `CanonicalDescription`.
Both `lb` and `ub` are `Vector{T}`, where `T` is the first argument in `parseqps`.

See also: `name`, `variables`, `objective`, `inequalities`, `inequalities`,
  `parseqps` and `canonical`.
"""
bounds(qp::MPSDescription)        = (qp.lb, qp.ub       )

function _compact(stream, ::MIME"text/plain", qp::MPSDescription)
  print(stream, "QP")
end

function _full(stream, ::MIME"text/plain", qp::MPSDescription)
  println(stream, _header(qp))
  println(stream, "---")
  println(stream, "  minimize    0.5*x'*Q*x + q₁'*x + q₂")
  println(stream, "  subject to       A*x = b")
  println(stream, "              c₁ ≤ C*x ≤ c₂")
  println(stream, "              lb ≤  x  ≤ ub")
  print(stream, "---")
end

_header(qp::MPSDescription) = string("QP instance: ", name(qp), ".")
_header(qp::MPSDescription{T,M}) where {T<:AbstractFloat,M<:AbstractSparseMatrix} =
  string("Sparse QP instance: ", name(qp), ".")

show(stream::IO, qp::MPSDescription)                          =
  show(stream, MIME("text/plain"), qp)
show(stream::IO, mime::MIME"text/plain", qp::MPSDescription)  =
  get(stream, :compact, false) ? _compact(stream, mime, qp) : _full(stream, mime, qp)


text:
#include("../ReadWrite/Header.jl")

"""
**SeisGeometry**

*Update headers with Geometry information beginning with shot and receiver coordinates (calculation of midpoints, absolute and vector offsets, azimuth, coordinate rotation, and binning of headers values)*

**IN**
* in: input filename (the .seish file is updated)
* ang=90 inline direction measured in degrees CC from East
* gamma=1 (vp/vs ratio for PS Asymptotic Conversion Point gathers (use gamma=1 for PP data))
* osx=0
* osy=0
* ogx=0
* ogy=0
* omx=0
* omy=0
* ohx=0
* ohy=0
* oh=0
* oaz=0
* dsx=1
* dsy=1
* dgx=1
* dgy=1
* dmx=1
* dmy=1
* dhx=1
* dhy=1
* dh=1
* daz=1

**OUT**

*Credits: A. Stanton, 2015*

"""

function SeisGeometry(in;ang=90,gamma=1,osx=0,osy=0,ogx=0,ogy=0,omx=0,omy=0,ohx=0,ohy=0,oh=0,oaz=0,dsx=1,dsy=1,dgx=1,dgy=1,dmx=1,dmy=1,dhx=1,dhy=1,dh=1,daz=1)

	rad2deg = 180/pi
	deg2rad = pi/180
	gammainv = 1/gamma
	if (ang > 90)
		ang2=-deg2rad*(ang-90)
	else
		ang2=deg2rad*(90-ang)
	end

	filename = ParseHeaderName(in)
	stream = open(filename,"r+")
	nhead = 27
	@compat nx = round(Int,filesize(stream)/(4*length(fieldnames(Header))))
	for j=1:nx
		h = GrabHeader(stream,j)
		h.hx = h.gx - h.sx
		h.hy = h.gy - h.sy
		h.h = sqrt((h.hx^2) + (h.hy^2))
		h.az = rad2deg*atan2((h.gy-h.sy),(h.gx-h.sx))
		if (h.az < 0)
			h.az += 360.0
		end
		h.mx = h.sx + h.hx/(1 + gammainv);
		h.my = h.sy + h.hy/(1 + gammainv);
		sx_rot = (h.sx-osx)*cos(ang2) - (h.sy-osy)*sin(ang2) + osx;
		sy_rot = (h.sx-osx)*sin(ang2) + (h.sy-osy)*cos(ang2) + osy;
		gx_rot = (h.gx-ogx)*cos(ang2) - (h.gy-ogy)*sin(ang2) + ogx;
		gy_rot = (h.gx-ogx)*sin(ang2) + (h.gy-ogy)*cos(ang2) + ogy;
		mx_rot = (h.mx-omx)*cos(ang2) - (h.my-omy)*sin(ang2) + omx;
		my_rot = (h.mx-omx)*sin(ang2) + (h.my-omy)*cos(ang2) + omy;
		hx_rot = (h.hx-ohx)*cos(ang2) - (h.hy-ohy)*sin(ang2) + ohx;
		hy_rot = (h.hx-ohx)*sin(ang2) + (h.hy-ohy)*cos(ang2) + ohy;
		h.isx = convert(Int32,round((sx_rot-osx)/dsx))
		h.isy = convert(Int32,round((sy_rot-osy)/dsy))
		h.igx = convert(Int32,round((gx_rot-ogx)/dgx))
		h.igy = convert(Int32,round((gy_rot-ogy)/dgy))
		h.imx = convert(Int32,round((mx_rot-omx)/dmx))
		h.imy = convert(Int32,round((my_rot-omy)/dmy))
		h.ihx = convert(Int32,round((hx_rot-ohx)/dhx))
		h.ihy = convert(Int32,round((hy_rot-ohy)/dhy))
		h.ih = convert(Int32,round((h.h-oh)/dh))
		h.iaz = convert(Int32,round((h.az-oaz)/daz))
		PutHeader(stream,h,j)
	end
	close(stream)

end


text:
"""
    LinearLayeredModel <: SeisModel1D

A `LinearLayeredModel` contains `n` points at which velocities are defined, with
linear interpolation between them.  Hence there are `n - 1` layers.

Discontinuities are represented by two layers with the same radii.
"""
struct LinearLayeredModel <: SeisModel1D
    "Earth radius in km"
    a :: Float64
    "Number of layers"
    n :: Int
    "Radii in km of top of each layer"
    r :: Vector{Float64}
    "Physical parameters"
    vp :: Vector{Float64}
    vs :: Vector{Float64}
    density :: Vector{Float64}
    "Optional anisotropic parameters"
    aniso :: Bool
    vph :: Vector{Float64}
    vpv :: Vector{Float64}
    vsh :: Vector{Float64}
    vsv :: Vector{Float64}
    eta :: Vector{Float64}
    "Optional attenuation parameters"
    attenuation :: Bool
    Qμ :: Vector{Float64}
    Qκ :: Vector{Float64}
end

"""
    LinearLayeredModel(; kwargs...) -> m

Construct a `LinearLayeredModel` by specifying keyword arguments.  The length of
arrays defines the number of layers and the presence of arguments determines
whether the model has density, anisotropy and attenuation.

Note that radii are given in km, velocities in km/s, densities in g/cm³
and frequency in Hz.

### Keyword arguments
#### Radii
- `r`: The radius of a set of nodes ranging from the centre of the planet to
  the surface, in km.  When the same value of `r` is repeated twice, a
  discontinuity is introduced at that radius.  The first value of `r` must
  be 0 km.  This is the only mandatory argument.
#### Isotropic properties
These must be given each as a `Vector` of values at each radial node.
Each parameter (`vp`, `density`, `Qμ`, etc.) must have the correct number of layers.
- `vp`: P-wave velocity
- `vs`: S-wave velocity
- `density`: Material density
#### Anisotropic parameters
If the model is anisotropic and the radial anisotropy parameters are all
provided, then the Voigt average velocities are constructed automatically
unless `vp` and/or `vs` are also given separately.
- `vpv`: Vertical P-wave velocity
- `vph`: Horizontal P-wave velocity
- `vsv`: Vertical S-wave velocity
- `vsh`: Horizontally-polarised S-wave velocity
- `eta`: The ratio C₁₃/(Vph - 2Vsv), where cᵢⱼ is the Voigt elasticity
  matrix, the 1-direction is horizontal and the 3-direction is radial
#### Attenuation parameters
ASCII versions (`Qmu`, `Qkappa`) of these exist for ease of use as well as
the Unicode versions (`Qμ`, `Qκ`).  If both are supplied, an error is thrown.
- `Qμ` or `Qmu`: Shear modulus quality factor
- `Qκ` or `Qkappa`: Bulk modulus quality factor
"""
function LinearLayeredModel(; r, vp=nothing, vs=nothing, density=nothing,
                            vph=nothing, vpv=nothing, vsh=nothing, vsv=nothing,
                            eta=nothing, Qμ=nothing, Qmu=nothing,
                            Qκ=nothing, Qkappa=nothing, fref=nothing)
    n = length(r)
    first(r) == 0 || throw(ArgumentError("first value of r must be 0 km"))
    all(x->x>=0, diff(r)) || throw(ArgumentError("radii must increase monotonically"))
    R = last(r)
    aniso = attenuation = false
    # Anisotropy
    if all(!isnothing, (vph, vpv, vsh, vsv, eta))
        aniso = true
        # Set isotropic velocities automatically if not provided
        if any(isnothing, (vp, vs))
            voigt_vels = voigt_velocities.(vpv, vsv, vph, vsh, eta)
            vp === nothing && (vp = first.(voigt_vels))
            vs === nothing && (vs = last.(voigt_vels))
        end
    end
    # Attenuation
    # Default to ASCII versions if given
    Qmu !== nothing && (Qμ = Qmu)
    Qkappa !== nothing && (Qκ = Qkappa)
    if all(!isnothing, (Qμ, Qκ))
        attenuation = true
    end
    # Convert to correct shape, and check for length
    vp = _set_coeffs_size(LinearLayeredModel, n, vp)
    vs = _set_coeffs_size(LinearLayeredModel, n, vs)
    density = _set_coeffs_size(LinearLayeredModel, n, density)
    vph = _set_coeffs_size(LinearLayeredModel, n, vph)
    vsh = _set_coeffs_size(LinearLayeredModel, n, vsh)
    vpv = _set_coeffs_size(LinearLayeredModel, n, vpv)
    vsv = _set_coeffs_size(LinearLayeredModel, n, vsv)
    eta = _set_coeffs_size(LinearLayeredModel, n, eta)
    Qμ = _set_coeffs_size(LinearLayeredModel, n, Qμ)
    Qκ = _set_coeffs_size(LinearLayeredModel, n, Qκ)
    LinearLayeredModel(R, n, r, vp, vs, density, aniso, vph, vpv,
        vsh, vsv, eta, attenuation, Qμ, Qκ)
end

_set_coeffs_size(::Type{LinearLayeredModel}, n, v::Nothing) = []
function _set_coeffs_size(::Type{LinearLayeredModel}, n, v::Vector)
    length(v) == n || throw(ArgumentError("number of values must be $n"))
    v
end

function findlayer(m::LinearLayeredModel, r::Real)
    r < 0. && error("Radius cannot be negative")
    r > m.a && error("Radius $r km is greater than Earth radius for model ($(m.a) km)")
    for i in 1:(length(m.r) - 1)
        r < m.r[i+1] && return i
    end
    return length(m.r) - 1
end

for sym in model_variables_SeisModel1D
    @eval function ($sym)(m::LinearLayeredModel, r::Real; depth::Bool=false)
        length(m.$sym) > 0 ||
            throw(ArgumentError("'$(split(string($sym), ".")[end])' not defined for model"))
        depth && (r = radius(m, r))
        ir = findlayer(m, r)
        r0 = m.r[ir]
        r1 = m.r[ir+1]
        (m.$sym[ir]*(r1 - r) + m.$sym[ir+1]*(r - r0))/(r1 - r0)
    end
end

function evaluate(m::LinearLayeredModel, field::Symbol, r; depth::Bool=false)
    y = getfield(m, field)
    length(y) > 0 ||
        throw(ArgumentError("'$field' not defined for model"))
    depth && (r = radius(m, r))
    ir = findlayer(m, r)
    r0 = m.r[ir]
    r1 = m.r[ir+1]
    (y[ir]*(r1 - r) + y[ir+1]*(r - r0))/(r1 - r0)
end

function mass(m::LinearLayeredModel, r; depth::Bool=false)
    hasdensity(m) || throw(ArgumentError("model does not contain density"))
    depth && (r = radius(m, r))
    l = findlayer(m, r)
    r *= 1.e3
    M = 0.
    for i in 1:l
        m.r[i] == m.r[i+1] && continue
        R0 = m.r[i]*1.e3
        R = i == l ? r : m.r[i+1]*1.e3
        dρ_dr = 1.e3*(m.density[i+1] - m.density[i])/(m.r[i+1]*1.e3 - R0)
        ρ0 = 1.e3*m.density[i] - dρ_dr*R0
        M += ρ0*(R^3 - R0^3)/3 + dρ_dr*(R^4 - R0^4)/4
    end
    4*π*M
end


text:
# Single line comments start with a hash (pound) symbol.
#= Multiline comments can be written
   by putting '#=' before the text  and '=#'
   after the text. They can also be nested.
=#

####################################################
## 1. Primitive Datatypes and Operators
####################################################

# Everything in Julia is an expression.

# There are several basic types of numbers.
typeof(3)       # => Int64
typeof(3.2)     # => Float64
typeof(2 + 1im) # => Complex{Int64}
typeof(2 // 3)  # => Rational{Int64}

# All of the normal infix operators are available.
1 + 1      # => 2
8 - 1      # => 7
10 * 2     # => 20
35 / 5     # => 7.0
10 / 2     # => 5.0  # dividing integers always results in a Float64
div(5, 2)  # => 2    # for a truncated result, use div
5 \ 35     # => 7.0
2^2        # => 4    # power, not bitwise xor
12 % 10    # => 2

# Enforce precedence with parentheses
(1 + 3) * 2  # => 8

# Bitwise Operators
~2         # => -3 # bitwise not
3 & 5      # => 1  # bitwise and
2 | 4      # => 6  # bitwise or
xor(2, 4)  # => 6  # bitwise xor
2 >>> 1    # => 1  # logical shift right
2 >> 1     # => 1  # arithmetic shift right
2 << 1     # => 4  # logical/arithmetic shift left

# Use the bitstring function to see the binary representation of a number.
bitstring(12345)
# => "0000000000000000000000000000000000000000000000000011000000111001"
bitstring(12345.0)
# => "0100000011001000000111001000000000000000000000000000000000000000"

# Boolean values are primitives
true
false

# Boolean operators
!true   # => false
!false  # => true
1 == 1  # => true
2 == 1  # => false
1 != 1  # => false
2 != 1  # => true
1 < 10  # => true
1 > 10  # => false
2 <= 2  # => true
2 >= 2  # => true
# Comparisons can be chained
1 < 2 < 3  # => true
2 < 3 < 2  # => false

# Strings are created with "
"This is a string."

# Character literals are written with '
'a'

# Strings are UTF8 encoded. Only if they contain only ASCII characters can
# they be safely indexed.
ascii("This is a string")[1]
# => 'T': ASCII/Unicode U+0054 (category Lu: Letter, uppercase)
# Julia indexes from 1
# Otherwise, iterating over strings is recommended (map, for loops, etc).

# $ can be used for string interpolation:
"2 + 2 = $(2 + 2)" # => "2 + 2 = 4"
# You can put any Julia expression inside the parentheses.

# Another way to format strings is the printf macro from the stdlib Printf.
using Printf
@printf "%d is less than %f\n" 4.5 5.3  # => 5 is less than 5.300000

# Printing is easy
println("I'm Julia. Nice to meet you!") # => I'm Julia. Nice to meet you!

# String can be compared lexicographically
"good" > "bye" # => true
"good" == "good" # => true
"1 + 2 = 3" == "1 + 2 = $(1 + 2)" # => true

####################################################
## 2. Variables and Collections
####################################################

# You don't declare variables before assigning to them.
someVar = 5  # => 5
someVar  # => 5

# Accessing a previously unassigned variable is an error
try
    someOtherVar  # => ERROR: UndefVarError: someOtherVar not defined
catch e
    println(e)
end

# Variable names start with a letter or underscore.
# After that, you can use letters, digits, underscores, and exclamation points.
SomeOtherVar123! = 6  # => 6

# You can also use certain unicode characters
☃ = 8  # => 8
# These are especially handy for mathematical notation
2 * π # => 6.283185307179586

# A note on naming conventions in Julia:
#
# * Word separation can be indicated by underscores ('_'), but use of
#   underscores is discouraged unless the name would be hard to read
#   otherwise.
#
# * Names of Types begin with a capital letter and word separation is shown
#   with CamelCase instead of underscores.
#
# * Names of functions and macros are in lower case, without underscores.
#
# * Functions that modify their inputs have names that end in !. These
#   functions are sometimes called mutating functions or in-place functions.

# Arrays store a sequence of values indexed by integers 1 through n:
a = Int64[] # => 0-element Array{Int64,1}

# 1-dimensional array literals can be written with comma-separated values.
b = [4, 5, 6] # => 3-element Array{Int64,1}: [4, 5, 6]
b = [4; 5; 6] # => 3-element Array{Int64,1}: [4, 5, 6]
b[1]    # => 4
b[end]  # => 6

# 2-dimensional arrays use space-separated values and semicolon-separated rows.
matrix = [1 2; 3 4] # => 2×2 Array{Int64,2}: [1 2; 3 4]

# Arrays of a particular type
b = Int8[4, 5, 6] # => 3-element Array{Int8,1}: [4, 5, 6]

# Add stuff to the end of a list with push! and append!
# By convention, the exclamation mark '!'' is appended to names of functions
# that modify their arguments
push!(a, 1)    # => [1]
push!(a, 2)    # => [1,2]
push!(a, 4)    # => [1,2,4]
push!(a, 3)    # => [1,2,4,3]
append!(a, b)  # => [1,2,4,3,4,5,6]

# Remove from the end with pop
pop!(b)  # => 6
b # => [4,5]

# Let's put it back
push!(b, 6)  # => [4,5,6]
b # => [4,5,6]

a[1]  # => 1  # remember that Julia indexes from 1, not 0!

# end is a shorthand for the last index. It can be used in any
# indexing expression
a[end]  # => 6

# we also have popfirst! and pushfirst!
popfirst!(a)  # => 1
a # => [2,4,3,4,5,6]
pushfirst!(a, 7)  # => [7,2,4,3,4,5,6]
a # => [7,2,4,3,4,5,6]

# Function names that end in exclamations points indicate that they modify
# their argument.
arr = [5,4,6]  # => 3-element Array{Int64,1}: [5,4,6]
sort(arr)   # => [4,5,6]
arr         # => [5,4,6]
sort!(arr)  # => [4,5,6]
arr         # => [4,5,6]

# Looking out of bounds is a BoundsError
try
    a[0]
    # => ERROR: BoundsError: attempt to access 7-element Array{Int64,1} at
    # index [0]
    # => Stacktrace:
    # =>  [1] getindex(::Array{Int64,1}, ::Int64) at .\array.jl:731
    # =>  [2] top-level scope at none:0
    # =>  [3] ...
    # => in expression starting at ...\LearnJulia.jl:180
    a[end + 1]
    # => ERROR: BoundsError: attempt to access 7-element Array{Int64,1} at
    # index [8]
    # => Stacktrace:
    # =>  [1] getindex(::Array{Int64,1}, ::Int64) at .\array.jl:731
    # =>  [2] top-level scope at none:0
    # =>  [3] ...
    # => in expression starting at ...\LearnJulia.jl:188
catch e
    println(e)
end

# Errors list the line and file they came from, even if it's in the standard
# library. You can look in the folder share/julia inside the julia folder to
# find these files.

# You can initialize arrays from ranges
a = [1:5;]  # => 5-element Array{Int64,1}: [1,2,3,4,5]
a2 = [1:5]  # => 1-element Array{UnitRange{Int64},1}: [1:5]

# You can look at ranges with slice syntax.
a[1:3]  # => [1, 2, 3]
a[2:end]  # => [2, 3, 4, 5]

# Remove elements from an array by index with splice!
arr = [3,4,5]
splice!(arr, 2) # => 4
arr # => [3,5]

# Concatenate lists with append!
b = [1,2,3]
append!(a, b) # => [1, 2, 3, 4, 5, 1, 2, 3]
a # => [1, 2, 3, 4, 5, 1, 2, 3]

# Check for existence in a list with in
in(1, a)  # => true

# Examine the length with length
length(a)  # => 8

# Tuples are immutable.
tup = (1, 2, 3)  # => (1,2,3)
typeof(tup) # => Tuple{Int64,Int64,Int64}
tup[1] # => 1
try
    tup[1] = 3
    # => ERROR: MethodError: no method matching
    # setindex!(::Tuple{Int64,Int64,Int64}, ::Int64, ::Int64)
catch e
    println(e)
end

# Many array functions also work on tuples
length(tup) # => 3
tup[1:2]    # => (1,2)
in(2, tup)  # => true

# You can unpack tuples into variables
a, b, c = (1, 2, 3)  # => (1,2,3)
a # => 1
b # => 2
c # => 3

# Tuples are created even if you leave out the parentheses
d, e, f = 4, 5, 6  # => (4,5,6)
d # => 4
e # => 5
f # => 6

# A 1-element tuple is distinct from the value it contains
(1,) == 1 # => false
(1) == 1  # => true

# Look how easy it is to swap two values
e, d = d, e  # => (5,4)
d # => 5
e # => 4

# Dictionaries store mappings
emptyDict = Dict()  # => Dict{Any,Any} with 0 entries

# You can create a dictionary using a literal
filledDict = Dict("one" => 1, "two" => 2, "three" => 3)
# => Dict{String,Int64} with 3 entries:
# =>  "two" => 2, "one" => 1, "three" => 3

# Look up values with []
filledDict["one"]  # => 1

# Get all keys
keys(filledDict)
# => Base.KeySet for a Dict{String,Int64} with 3 entries. Keys:
# =>  "two", "one", "three"
# Note - dictionary keys are not sorted or in the order you inserted them.

# Get all values
values(filledDict)
# => Base.ValueIterator for a Dict{String,Int64} with 3 entries. Values:
# =>  2, 1, 3
# Note - Same as above regarding key ordering.

# Check for existence of keys in a dictionary with in, haskey
in(("one" => 1), filledDict)  # => true
in(("two" => 3), filledDict)  # => false
haskey(filledDict, "one")     # => true
haskey(filledDict, 1)         # => false

# Trying to look up a non-existent key will raise an error
try
    filledDict["four"]  # => ERROR: KeyError: key "four" not found
catch e
    println(e)
end

# Use the get method to avoid that error by providing a default value
# get(dictionary, key, defaultValue)
get(filledDict, "one", 4)   # => 1
get(filledDict, "four", 4)  # => 4

# Use Sets to represent collections of unordered, unique values
emptySet = Set()  # => Set(Any[])
# Initialize a set with values
filledSet = Set([1, 2, 2, 3, 4])  # => Set([4, 2, 3, 1])

# Add more values to a set
push!(filledSet, 5)  # => Set([4, 2, 3, 5, 1])

# Check if the values are in the set
in(2, filledSet)   # => true
in(10, filledSet)  # => false

# There are functions for set intersection, union, and difference.
otherSet = Set([3, 4, 5, 6])         # => Set([4, 3, 5, 6])
intersect(filledSet, otherSet)      # => Set([4, 3, 5])
union(filledSet, otherSet)          # => Set([4, 2, 3, 5, 6, 1])
setdiff(Set([1,2,3,4]), Set([2,3,5])) # => Set([4, 1])

####################################################
## 3. Control Flow
####################################################

# Let's make a variable
someVar = 5

# Here is an if statement. Indentation is not meaningful in Julia.
if someVar > 10
    println("someVar is totally bigger than 10.")
elseif someVar < 10    # This elseif clause is optional.
    println("someVar is smaller than 10.")
else                    # The else clause is optional too.
    println("someVar is indeed 10.")
end
# => prints "some var is smaller than 10"

# For loops iterate over iterables.
# Iterable types include Range, Array, Set, Dict, and AbstractString.
for animal = ["dog", "cat", "mouse"]
    println("$animal is a mammal")
    # You can use $ to interpolate variables or expression into strings
end
# => dog is a mammal
# => cat is a mammal
# => mouse is a mammal

# You can use 'in' instead of '='.
for animal in ["dog", "cat", "mouse"]
    println("$animal is a mammal")
end
# => dog is a mammal
# => cat is a mammal
# => mouse is a mammal

for pair in Dict("dog" => "mammal", "cat" => "mammal", "mouse" => "mammal")
    from, to = pair
    println("$from is a $to")
end
# => mouse is a mammal
# => cat is a mammal
# => dog is a mammal

for (k, v) in Dict("dog" => "mammal", "cat" => "mammal", "mouse" => "mammal")
    println("$k is a $v")
end
# => mouse is a mammal
# => cat is a mammal
# => dog is a mammal

# While loops loop while a condition is true
let x = 0
    while x < 4
        println(x)
        x += 1  # Shorthand for x = x + 1
    end
end
# => 0
# => 1
# => 2
# => 3

# Handle exceptions with a try/catch block
try
    error("help")
catch e
    println("caught it $e")
end
# => caught it ErrorException("help")

####################################################
## 4. Functions
####################################################

# The keyword 'function' creates new functions
# function name(arglist)
#   body...
# end
function add(x, y)
    println("x is $x and y is $y")

    # Functions return the value of their last statement
    x + y
end

add(5, 6)
# => x is 5 and y is 6
# => 11

# Compact assignment of functions
f_add(x, y) = x + y  # => f_add (generic function with 1 method)
f_add(3, 4)  # => 7

# Function can also return multiple values as tuple
fn(x, y) = x + y, x - y # => fn (generic function with 1 method)
fn(3, 4)  # => (7, -1)

# You can define functions that take a variable number of
# positional arguments
function varargs(args...)
    return args
    # use the keyword return to return anywhere in the function
end
# => varargs (generic function with 1 method)

varargs(1, 2, 3)  # => (1,2,3)

# The ... is called a splat.
# We just used it in a function definition.
# It can also be used in a function call,
# where it will splat an Array or Tuple's contents into the argument list.
add([5,6]...)  # this is equivalent to add(5,6)

x = (5, 6)  # => (5,6)
add(x...)  # this is equivalent to add(5,6)

# You can define functions with optional positional arguments
function defaults(a, b, x=5, y=6)
    return "$a $b and $x $y"
end
# => defaults (generic function with 3 methods)

defaults('h', 'g')  # => "h g and 5 6"
defaults('h', 'g', 'j')  # => "h g and j 6"
defaults('h', 'g', 'j', 'k')  # => "h g and j k"
try
    defaults('h')  # => ERROR: MethodError: no method matching defaults(::Char)
    defaults()  # => ERROR: MethodError: no method matching defaults()
catch e
    println(e)
end

# You can define functions that take keyword arguments
function keyword_args(;k1=4, name2="hello")  # note the ;
    return Dict("k1" => k1, "name2" => name2)
end
# => keyword_args (generic function with 1 method)

keyword_args(name2="ness")  # => ["name2"=>"ness", "k1"=>4]
keyword_args(k1="mine")     # => ["name2"=>"hello", "k1"=>"mine"]
keyword_args()              # => ["name2"=>"hello", "k1"=>4]

# You can combine all kinds of arguments in the same function
function all_the_args(normalArg, optionalPositionalArg=2; keywordArg="foo")
    println("normal arg: $normalArg")
    println("optional arg: $optionalPositionalArg")
    println("keyword arg: $keywordArg")
end
# => all_the_args (generic function with 2 methods)

all_the_args(1, 3, keywordArg=4)
# => normal arg: 1
# => optional arg: 3
# => keyword arg: 4

# Julia has first class functions
function create_adder(x)
    adder = function (y)
        return x + y
    end
    return adder
end
# => create_adder (generic function with 1 method)

# This is "stabby lambda syntax" for creating anonymous functions
(x -> x > 2)(3)  # => true

# This function is identical to create_adder implementation above.
function create_adder(x)
    y -> x + y
end
# => create_adder (generic function with 1 method)

# You can also name the internal function, if you want
function create_adder(x)
    function adder(y)
        x + y
    end
    adder
end
# => create_adder (generic function with 1 method)

add_10 = create_adder(10) # => (::getfield(Main, Symbol("#adder#11")){Int64})
                          # (generic function with 1 method)
add_10(3) # => 13

# There are built-in higher order functions
map(add_10, [1,2,3])  # => [11, 12, 13]
filter(x -> x > 5, [3, 4, 5, 6, 7])  # => [6, 7]

# We can use list comprehensions
[add_10(i) for i = [1, 2, 3]]   # => [11, 12, 13]
[add_10(i) for i in [1, 2, 3]]  # => [11, 12, 13]
[x for x in [3, 4, 5, 6, 7] if x > 5] # => [6, 7]

####################################################
## 5. Types
####################################################

# Julia has a type system.
# Every value has a type; variables do not have types themselves.
# You can use the `typeof` function to get the type of a value.
typeof(5)  # => Int64

# Types are first-class values
typeof(Int64)     # => DataType
typeof(DataType)  # => DataType
# DataType is the type that represents types, including itself.

# Types are used for documentation, optimizations, and dispatch.
# They are not statically checked.

# Users can define types
# They are like records or structs in other languages.
# New types are defined using the `struct` keyword.

# struct Name
#   field::OptionalType
#   ...
# end
struct Tiger
    taillength::Float64
    coatcolor  # not including a type annotation is the same as `::Any`
end

# The default constructor's arguments are the properties
# of the type, in the order they are listed in the definition
tigger = Tiger(3.5, "orange")  # => Tiger(3.5,"orange")

# The type doubles as the constructor function for values of that type
sherekhan = typeof(tigger)(5.6, "fire")  # => Tiger(5.6,"fire")

# These struct-style types are called concrete types
# They can be instantiated, but cannot have subtypes.
# The other kind of types is abstract types.

# abstract Name
abstract type Cat end  # just a name and point in the type hierarchy

# Abstract types cannot be instantiated, but can have subtypes.
# For example, Number is an abstract type
subtypes(Number)  # => 2-element Array{Any,1}:
                  # =>  Complex
                  # =>  Real
subtypes(Cat)  # => 0-element Array{Any,1}

# AbstractString, as the name implies, is also an abstract type
subtypes(AbstractString)  # => 4-element Array{Any,1}:
                          # =>  String
                          # =>  SubString
                          # =>  SubstitutionString
                          # =>  Test.GenericString

# Every type has a super type; use the `supertype` function to get it.
typeof(5) # => Int64
supertype(Int64)    # => Signed
supertype(Signed)   # => Integer
supertype(Integer)  # => Real
supertype(Real)     # => Number
supertype(Number)   # => Any
supertype(supertype(Signed))  # => Real
supertype(Any)      # => Any
# All of these type, except for Int64, are abstract.
typeof("fire")      # => String
supertype(String)   # => AbstractString
# Likewise here with String
supertype(SubString)  # => AbstractString

# <: is the subtyping operator
struct Lion <: Cat  # Lion is a subtype of Cat
    maneColor
    roar::AbstractString
end

# You can define more constructors for your type
# Just define a function of the same name as the type
# and call an existing constructor to get a value of the correct type
Lion(roar::AbstractString) = Lion("green", roar)
# This is an outer constructor because it's outside the type definition

struct Panther <: Cat  # Panther is also a subtype of Cat
    eyeColor
    Panther() = new("green")
    # Panthers will only have this constructor, and no default constructor.
end
# Using inner constructors, like Panther does, gives you control
# over how values of the type can be created.
# When possible, you should use outer constructors rather than inner ones.

####################################################
## 6. Multiple-Dispatch
####################################################

# In Julia, all named functions are generic functions
# This means that they are built up from many small methods
# Each constructor for Lion is a method of the generic function Lion.

# For a non-constructor example, let's make a function meow:

# Definitions for Lion, Panther, Tiger
function meow(animal::Lion)
    animal.roar  # access type properties using dot notation
end

function meow(animal::Panther)
    "grrr"
end

function meow(animal::Tiger)
    "rawwwr"
end

# Testing the meow function
meow(tigger)  # => "rawwwr"
meow(Lion("brown", "ROAAR"))  # => "ROAAR"
meow(Panther()) # => "grrr"

# Review the local type hierarchy
Tiger   <: Cat  # => false
Lion    <: Cat  # => true
Panther <: Cat  # => true

# Defining a function that takes Cats
function pet_cat(cat::Cat)
    println("The cat says $(meow(cat))")
end
# => pet_cat (generic function with 1 method)

pet_cat(Lion("42")) # => The cat says 42
try
    pet_cat(tigger) # => ERROR: MethodError: no method matching pet_cat(::Tiger)
catch e
    println(e)
end

# In OO languages, single dispatch is common;
# this means that the method is picked based on the type of the first argument.
# In Julia, all of the argument types contribute to selecting the best method.

# Let's define a function with more arguments, so we can see the difference
function fight(t::Tiger, c::Cat)
    println("The $(t.coatcolor) tiger wins!")
end
# => fight (generic function with 1 method)

fight(tigger, Panther())  # => The orange tiger wins!
fight(tigger, Lion("ROAR")) # => The orange tiger wins!

# Let's change the behavior when the Cat is specifically a Lion
fight(t::Tiger, l::Lion) = println("The $(l.maneColor)-maned lion wins!")
# => fight (generic function with 2 methods)

fight(tigger, Panther())  # => The orange tiger wins!
fight(tigger, Lion("ROAR")) # => The green-maned lion wins!

# We don't need a Tiger in order to fight
fight(l::Lion, c::Cat) = println("The victorious cat says $(meow(c))")
# => fight (generic function with 3 methods)

fight(Lion("balooga!"), Panther())  # => The victorious cat says grrr
try
    fight(Panther(), Lion("RAWR"))
    # => ERROR: MethodError: no method matching fight(::Panther, ::Lion)
    # => Closest candidates are:
    # =>   fight(::Tiger, ::Lion) at ...
    # =>   fight(::Tiger, ::Cat) at ...
    # =>   fight(::Lion, ::Cat) at ...
    # => ...
catch e
    println(e)
end

# Also let the cat go first
fight(c::Cat, l::Lion) = println("The cat beats the Lion")
# => fight (generic function with 4 methods)

# This warning is because it's unclear which fight will be called in:
try
    fight(Lion("RAR"), Lion("brown", "rarrr"))
    # => ERROR: MethodError: fight(::Lion, ::Lion) is ambiguous. Candidates:
    # =>   fight(c::Cat, l::Lion) in Main at ...
    # =>   fight(l::Lion, c::Cat) in Main at ...
    # => Possible fix, define
    # =>   fight(::Lion, ::Lion)
    # => ...
catch e
    println(e)
end
# The result may be different in other versions of Julia

fight(l::Lion, l2::Lion) = println("The lions come to a tie")
# => fight (generic function with 5 methods)
fight(Lion("RAR"), Lion("brown", "rarrr"))  # => The lions come to a tie

# Under the hood
# You can take a look at the llvm  and the assembly code generated.

square_area(l) = l * l  # square_area (generic function with 1 method)

square_area(5)  # => 25

# What happens when we feed square_area an integer?
code_native(square_area, (Int32,), syntax = :intel)
    #         .text
    # ; Function square_area {
    # ; Location: REPL[116]:1       # Prologue
    #         push    rbp
    #         mov     rbp, rsp
    # ; Function *; {
    # ; Location: int.jl:54
    #         imul    ecx, ecx      # Square l and store the result in ECX
    # ;}
    #         mov     eax, ecx
    #         pop     rbp           # Restore old base pointer
    #         ret                   # Result will still be in EAX
    #         nop     dword ptr [rax + rax]
    # ;}

code_native(square_area, (Float32,), syntax = :intel)
    #         .text
    # ; Function square_area {
    # ; Location: REPL[116]:1
    #         push    rbp
    #         mov     rbp, rsp
    # ; Function *; {
    # ; Location: float.jl:398
    #         vmulss  xmm0, xmm0, xmm0  # Scalar single precision multiply (AVX)
    # ;}
    #         pop     rbp
    #         ret
    #         nop     word ptr [rax + rax]
    # ;}

code_native(square_area, (Float64,), syntax = :intel)
    #         .text
    # ; Function square_area {
    # ; Location: REPL[116]:1
    #         push    rbp
    #         mov     rbp, rsp
    # ; Function *; {
    # ; Location: float.jl:399
    #         vmulsd  xmm0, xmm0, xmm0  # Scalar double precision multiply (AVX)
    # ;}
    #         pop     rbp
    #         ret
    #         nop     word ptr [rax + rax]
    # ;}

# Note that julia will use floating point instructions if any of the
# arguments are floats.
# Let's calculate the area of a circle
circle_area(r) = pi * r * r     # circle_area (generic function with 1 method)
circle_area(5)  # 78.53981633974483

code_native(circle_area, (Int32,), syntax = :intel)
    #         .text
    # ; Function circle_area {
    # ; Location: REPL[121]:1
    #         push    rbp
    #         mov     rbp, rsp
    # ; Function *; {
    # ; Location: operators.jl:502
    # ; Function *; {
    # ; Location: promotion.jl:314
    # ; Function promote; {
    # ; Location: promotion.jl:284
    # ; Function _promote; {
    # ; Location: promotion.jl:261
    # ; Function convert; {
    # ; Location: number.jl:7
    # ; Function Type; {
    # ; Location: float.jl:60
    #         vcvtsi2sd       xmm0, xmm0, ecx     # Load integer (r) from memory
    #         movabs  rax, 497710928              # Load pi
    # ;}}}}}
    # ; Function *; {
    # ; Location: float.jl:399
    #         vmulsd  xmm1, xmm0, qword ptr [rax] # pi * r
    #         vmulsd  xmm0, xmm1, xmm0            # (pi * r) * r
    # ;}}
    #         pop     rbp
    #         ret
    #         nop     dword ptr [rax]
    # ;}

code_native(circle_area, (Float64,), syntax = :intel)
    #         .text
    # ; Function circle_area {
    # ; Location: REPL[121]:1
    #         push    rbp
    #         mov     rbp, rsp
    #         movabs  rax, 497711048
    # ; Function *; {
    # ; Location: operators.jl:502
    # ; Function *; {
    # ; Location: promotion.jl:314
    # ; Function *; {
    # ; Location: float.jl:399
    #         vmulsd  xmm1, xmm0, qword ptr [rax]
    # ;}}}
    # ; Function *; {
    # ; Location: float.jl:399
    #         vmulsd  xmm0, xmm1, xmm0
    # ;}
    #         pop     rbp
    #         ret
    #         nop     dword ptr [rax + rax]
    # ;}

max_stars_repo_name: jameszhan/notes-ml

max_issues_repo_name: jameszhan/notes-ml

max_forks_repo_name: jameszhan/notes-ml


text:
println("Welcome to rate2.jl, where rates for CIRIS are plotted")
using DataFrames
using Gadfly

path = "/scratch/cns7ae/losalamos9/"
outfile = "rates-plot.pdf"
xmax = 1e15
xmin = 1e13
speciesName = "O3"
speciesNum  = 7

plottheme = Gadfly.Theme(
    highlight_width=0pt,
    default_point_size=2pt,
    line_width=2pt,
    minor_label_font_size=13pt,
    minor_label_color=colorant"black",
    major_label_font_size=15pt,
    major_label_color=colorant"black",
    key_label_color=colorant"black",
    key_title_color=colorant"black",
    key_title_font_size=8pt,
    key_label_font_size=6pt,
    grid_color=colorant"black")

o_prod_filename  = "o_prod.wsv"
o2_prod_filename = "o2_prod.wsv"
o3_prod_filename = "o3_prod.wsv"
o_dest_filename  = "o_dest.wsv"
o2_dest_filename = "o2_dest.wsv"
o3_dest_filename = "o3_dest.wsv"

ostar_prod_filename  = "ostar_prod.wsv"
o2star_prod_filename = "o2star_prod.wsv"
o3star_prod_filename = "o3star_prod.wsv"
ostar_dest_filename  = "ostar_dest.wsv"
o2star_dest_filename = "o2star_dest.wsv"
o3star_dest_filename = "o3star_dest.wsv"

o_prod_output = path*o_prod_filename
o2_prod_output = path*o2_prod_filename
o3_prod_output = path*o3_prod_filename
o_dest_output = path*o_dest_filename
o2_dest_output = path*o2_dest_filename
o3_dest_output = path*o3_dest_filename

ostar_prod_output = path*ostar_prod_filename
#o2star_prod_output = path*o2star_prod_filename
o3star_prod_output = path*o3star_prod_filename
ostar_dest_output = path*ostar_dest_filename
o2star_dest_output = path*o2star_dest_filename
o3star_dest_output = path*o3star_dest_filename

#=
Generate the dataframes containing the input and output dataframes
=#
println("Generating initial dataframes")

o_prod_reactions = readtable(
    o_prod_output,
    header = true,
    names = [:Count,:Time, :Fluence, :ID, :Label],
    eltypes = [Int64,Float64,Float64,Int64,String]
)

println("Generated o prod")

println("Now reading $o2_prod_output")
o2_prod_reactions = readtable(
    o2_prod_output,
    header = true,
    names = [:Count,:Time, :Fluence, :ID, :Label],
    eltypes = [Int64,Float64,Float64,Int64,String]
)

println("Generated o2 prod")

o3_prod_reactions = readtable(
    o3_prod_output,
    header = true,
    names = [:Count,:Time, :Fluence, :ID, :Label],
    eltypes = [Int64,Float64,Float64,Int64,String]
)

println("Generated o3 prod")

o_dest_reactions = readtable(
    o_dest_output,
    header = true,
    names = [:Count,:Time, :Fluence, :ID, :Label],
    eltypes = [Int64,Float64,Float64,Int64,String]
)

println("Generated o dest")

o2_dest_reactions = readtable(
    o2_dest_output,
    header = true,
    names = [:Count,:Time, :Fluence, :ID, :Label],
    eltypes = [Int64,Float64,Float64,Int64,String]
)

println("Generated o2 dest")

o3_dest_reactions = readtable(
    o3_dest_output,
    header = true,
    names = [:Count,:Time, :Fluence, :ID, :Label],
    eltypes = [Int64,Float64,Float64,Int64,String]
)

println("Generated o3 dest")

println("Now plotting p1")

p1 = plot(
#    destdf[xmax.>destdf[:Fluence].>xmin,:],
    o_prod_reactions,
    x="Fluence",
    y="Count",
    color="Label",
    Geom.line,
#    Geom.smooth(method=:lm,smoothing=0.9),
    Guide.title("O Prod. Reactions"),
    Guide.xlabel("Fluence"),
    Guide.colorkey("Reaction #"),
    Guide.ylabel("Count"),
#    Coord.Cartesian(xmin=fluenceArr[1],xmax=fluenceArr[size(fluenceArr,1)])
    Scale.x_log10(minvalue=xmin,maxvalue=xmax),
#    Scale.y_log10
)

push!(p1,plottheme)
println("Plotted p1")

p2 = plot(
#    destdf[xmax.>destdf[:Fluence].>xmin,:],
    o2_prod_reactions,
    x="Fluence",
    y="Count",
    color="Label",
    Geom.line,
#    Geom.smooth(method=:lm,smoothing=0.9),
    Guide.title("O2 Prod. Reactions"),
    Guide.xlabel("Fluence"),
    Guide.colorkey("Reaction #"),
    Guide.ylabel("Count"),
#    Coord.Cartesian(xmin=fluenceArr[1],xmax=fluenceArr[size(fluenceArr,1)])
    Scale.x_log10(minvalue=xmin,maxvalue=xmax),
#    Scale.y_log10
)

push!(p2,plottheme)
println("Plotted p2")

p3 = plot(
#    destdf[xmax.>destdf[:Fluence].>xmin,:],
    o3_prod_reactions,
    x="Fluence",
    y="Count",
    color="Label",
    Geom.line,
#    Geom.smooth(method=:loess,smoothing=0.9),
    Guide.title("O3 Prod. Reactions"),
    Guide.xlabel("Fluence"),
    Guide.colorkey("Reaction #"),
    Guide.ylabel("Count"),
#    Coord.Cartesian(xmin=fluenceArr[1],xmax=fluenceArr[size(fluenceArr,1)])
    Scale.x_log10(minvalue=xmin,maxvalue=xmax),
#    Scale.y_log10
)

push!(p3,plottheme)
println("Plotted p3")

r1 = plot(
#    destdf[xmax.>destdf[:Fluence].>xmin,:],
    o_dest_reactions,
    x="Fluence",
    y="Count",
    color="Label",
    Geom.line,
#    Geom.smooth(method=:lm,smoothing=0.9),
    Guide.title("O Dest. Reactions"),
    Guide.xlabel("Fluence"),
    Guide.colorkey("Reaction #"),
    Guide.ylabel("Count"),
#    Coord.Cartesian(xmin=fluenceArr[1],xmax=fluenceArr[size(fluenceArr,1)])
    Scale.x_log10(minvalue=xmin,maxvalue=xmax),
#    Scale.y_log10
)

push!(r1,plottheme)
println("Plotted r1")

r2 = plot(
#    destdf[xmax.>destdf[:Fluence].>xmin,:],
    o2_dest_reactions,
    x="Fluence",
    y="Count",
    color="Label",
    Geom.line,
#    Geom.smooth(method=:lm,smoothing=0.9),
    Guide.title("O2 Dest. Reactions"),
    Guide.xlabel("Fluence"),
    Guide.colorkey("Reaction #"),
    Guide.ylabel("Count"),
#    Coord.Cartesian(xmin=fluenceArr[1],xmax=fluenceArr[size(fluenceArr,1)])
    Scale.x_log10(minvalue=xmin,maxvalue=xmax),
#    Scale.y_log10
)

push!(r2,plottheme)
println("Plotted r2")

r3 = plot(
#    destdf[xmax.>destdf[:Fluence].>xmin,:],
    o3_dest_reactions,
    x="Fluence",
    y="Count",
    color="Label",
    Geom.line,
#    Geom.smooth(method=:lm,smoothing=0.9),
    Guide.title("O3 Dest. Reactions"),
    Guide.xlabel("Fluence"),
    Guide.colorkey("Reaction #"),
    Guide.ylabel("Count"),
#    Coord.Cartesian(xmin=fluenceArr[1],xmax=fluenceArr[size(fluenceArr,1)])
    Scale.x_log10(minvalue=xmin,maxvalue=xmax),
#    Scale.y_log10
)

push!(r3,plottheme)
println("Plotted r3")

# Now plotting excited species reactions

ostar_prod_reactions = readtable(
    ostar_prod_output,
    header = true,
    names = [:Count,:Time, :Fluence, :ID, :Label],
    eltypes = [Int64,Float64,Float64,Int64,String]
)

println("Generated o* prod")

ostar_dest_reactions = readtable(
    ostar_dest_output,
    header = true,
    names = [:Count,:Time, :Fluence, :ID, :Label],
    eltypes = [Int64,Float64,Float64,Int64,String]
)

println("Generated o* dest")

println("Generated dataframe of reactions output")

s1 = plot(
#    destdf[xmax.>destdf[:Fluence].>xmin,:],
    ostar_prod_reactions,
    x="Fluence",
    y="Count",
    color="Label",
    Geom.line,
#    Geom.smooth(method=:lm,smoothing=0.9),
    Guide.title("O* Prod. Reactions"),
    Guide.xlabel("Fluence"),
    Guide.colorkey("Reaction #"),
    Guide.ylabel("Count"),
#    Coord.Cartesian(xmin=fluenceArr[1],xmax=fluenceArr[size(fluenceArr,1)])
    Scale.x_log10(minvalue=xmin,maxvalue=xmax),
#    Scale.y_log10
)

push!(s1,plottheme)
println("Plotted s1")

t1 = plot(
#    destdf[xmax.>destdf[:Fluence].>xmin,:],
    ostar_dest_reactions,
    x="Fluence",
    y="Count",
    color="Label",
    Geom.line,
#    Geom.smooth(method=:lm,smoothing=0.9),
    Guide.title("O* Dest. Reactions"),
    Guide.xlabel("Fluence"),
    Guide.colorkey("Reaction #"),
    Guide.ylabel("Count"),
#    Coord.Cartesian(xmin=fluenceArr[1],xmax=fluenceArr[size(fluenceArr,1)])
    Scale.x_log10(minvalue=xmin,maxvalue=xmax),
#    Scale.y_log10
)

push!(t1,plottheme)
println("Plotted t1")

#println("Generating final output pdf")
#draw(PDF(path*"o_reactions.pdf",6inch,9inch),vstack(p1,r1))
#draw(PDF(path*"o2_reactions.pdf",6inch,9inch),vstack(p2,r2))
#draw(PDF(path*"o3_reactions.pdf",6inch,9inch),vstack(p3,r3))
#draw(PDF(path*"o*_reactions.pdf",6inch,9inch),vstack(s1,t1))
#draw(PDF(path*"o2*_reactions.pdf",6inch,9inch),vstack(s2,t2))
#draw(PDF(path*"o3*_reactions.pdf",6inch,9inch),vstack(s3,t3))

println("Ending rate2.jl!")
#quit()