text:
## Handling of the topological archive internally used to recognize topologies.
using Pkg.Artifacts

const CRYSTAL_NETS_VERSION = v"0.1.0"
const arc_location = artifact"archives"

"""
    const CRYSTAL_NETS_ARCHIVE::Dict{String,String}

The archive used to recognize known topologies.

You probably don't need to access it directly: rely on [`recognize_topology`](@ref) to read
and the various archive functions like [`add_to_current_archive!`](@ref) to write.
"""
const CRYSTAL_NETS_ARCHIVE = if isdir(arc_location) && !isempty(readdir(arc_location))
    flag, parsed = parse_arcs(arc_location)
    if !flag
        error("""CrystalNets.jl appears to have a broken installation (the archive version is older than that package's).
        Please rebuild CrystalNets.jl with `import Pkg; Pkg.build("CrystalNets")`.
        """)
    end
    parsed
else
    error("""CrystalNets.jl appears to have a broken installation (missing default archive).
          Please rebuild CrystalNets.jl with `import Pkg; Pkg.build("CrystalNets")`.

          This issue may be due to an error while using `CrystalNets.clean_default_archive!(...)`.
          If you have not encountered any such error and you did not modify or erase the archive located at $arc_location please open an issue at https://github.com/coudertlab/CrystalNets.jl/issues/new
          """)
end

"""
    const REVERSE_CRYSTAL_NETS_ARCHIVE::Dict{String,String}

Reverse of [`CRYSTAL_NETS_ARCHIVE`](@ref).

Can be used to query the topological genome of known nets, as in:
```jldoctest
julia> REVERSE_CRYSTAL_NETS_ARCHIVE["dia"]
"3 1 2 0 0 0 1 2 0 0 1 1 2 0 1 0 1 2 1 0 0"

julia> topological_genome(CrystalNet(PeriodicGraph(ans)))
dia
```

!!! note
    It is also possible to directly access the topological genome as a `PeriodicGraph`
    by parsing the name as a [`TopologicalGenome`](@ref):
    ```jldoctest
    julia> parse(TopologicalGenome, "pcu").genome
    PeriodicGraph3D(1, PeriodicEdge3D[(1, 1, (0,0,1)), (1, 1, (0,1,0)), (1, 1, (1,0,0))])

    julia> string(parse(TopologicalGenome, "nbo").genome) == REVERSE_CRYSTAL_NETS_ARCHIVE["nbo"]
    true
    ```
"""
const REVERSE_CRYSTAL_NETS_ARCHIVE = Dict{String,String}(id => (startswith(key, "unstable") ? key[10:end] : key) for (key, id) in CRYSTAL_NETS_ARCHIVE)

export clean_default_archive!,
       set_default_archive!,
       empty_default_archive!,
       change_current_archive!,
       refresh_current_archive!,
       add_to_current_archive!,
       make_archive,
       REVERSE_CRYSTAL_NETS_ARCHIVE

function _reset_archive!()
    global CRYSTAL_NETS_ARCHIVE
    global REVERSE_CRYSTAL_NETS_ARCHIVE
    global arc_location
    empty!(CRYSTAL_NETS_ARCHIVE)
    merge!(CRYSTAL_NETS_ARCHIVE, last(parse_arcs(arc_location)))
    empty!(REVERSE_CRYSTAL_NETS_ARCHIVE)
    merge!(REVERSE_CRYSTAL_NETS_ARCHIVE, Dict{String,String}(last(x) => first(x) for x in CRYSTAL_NETS_ARCHIVE))
    nothing
end

function validate_archive(custom_arc)::Dict{String,String}
    arc = try
        Serialization.deserialize(custom_arc)
    catch
        nothing
    end
    if arc isa Dict{String,String} && !isempty(arc) && isnumeric(first(first(keys(arc))))
        @ifwarn @info "Processing input as a serialized CrystalNets.jl archive"
    else
        try
            flag, parsed = parse_arc(custom_arc)
            if flag
                arc = parsed
            else
                @ifwarn begin
                    @info "Processing input as a topological .arc file"
                    if occursin("Made by CrystalNets.jl", readline(custom_arc))
                        @info "This archive was generated by an older version of CrystalNets."
                    end
                    @info "Keys will be converted to the topological genome used by CrystalNets. This may take a while."
                end
                arc_per_thread = [Pair{String,String}[] for _ in 1:nthreads()]
                Threads.@threads for (key, id) in collect(parsed)
                    _g = topological_genome(CrystalNet(PeriodicGraph(key)))
                    if _g.unstable || !isempty(_g.error)
                        if _g.unstable
                            println(stderr, "Net ", id, " is unstable.")
                        else
                            println(stderr, "Failed for net ", id, " with error:", g.error)
                        end
                        continue
                    end
                    genome = string(_g.genome)
                    if !isempty(genome)
                        push!(arc_per_thread[threadid()], (genome => id))
                    end
                end
                arc = Dict{String,String}(pop!(arc_per_thread))
                for _arc in arc_per_thread
                    merge!(arc, Dict{String,String}(_arc))
                end
            end
        catch e
            print(stderr, """
            Impossible to parse input as a topological archive. Please use a format
            similar to that of the RCSR Systre .arc file, with for each entry at
            least the "id" and "key" fields.

            This error may also occur if the given archive was empty. If you
            wish to set an empty archive, use `CrystalNets.empty_default_archive!()`
            """, '\n', "Encountered error: ")
            rethrow()
            # showerror(stderr, e)
            # Base.display_error(stderr, e, catch_backtrace())
        end
    end
    return arc
end

"""
    clean_default_archive!(custom_arc; validate=true, refresh=true)

Erase the default archive used by CrystalNets.jl to recognize known topologies
and replace it with a new one from the file located at `custom_arc`.

The `validate` parameter controls whether the new file is checked and converted
to a format usable by CrystalNets.jl. If unsure, leave it set.

The `refresh` optional parameter controls whether the current archive should be
replaced by the new default one.

!!! warning
    This archive will be kept and used for subsequent runs of CrystalNets.jl, even
    if you restart your Julia session.

    To only change the archive for the current session, use [`change_current_archive!(custom_arc)`](@ref change_current_archive!).

    See also [`refresh_current_archive!`](@ref) for similar uses.

!!! warning
    The previous default archive cannot be recovered afterwards, so make sure to
    keep a copy if necessary. The default archive is the set of ".arc" files located
    at `joinpath(dirname(dirname(pathof(CrystalNets))), "archives")`.
"""
function clean_default_archive!(custom_arc=nothing; validate=true, refresh=true, name="new")
    rm(arc_location; recursive=true)
    mkdir(arc_location)
    if validate
        arc = validate_archive(custom_arc)
        export_arc(joinpath(arc_location, name*".arc"), false, arc)
    else
        cp(custom_arc, arc_location)
    end
    if refresh
        refresh_current_archive!()
    end
    nothing
end

"""
    set_default_archive!()

Set the current archive as the new default archive.

!!! warning
    This archive will be kept and used for subsequent runs of CrystalNets.jl, even
    if you restart your Julia session.
"""
function set_default_archive!(name="new")
    global CRYSTAL_NETS_ARCHIVE
    export_arc(joinpath(arc_location, name*".arc"))
end

"""
    empty_default_archive!(; refresh=true)

Empty the default archive. This will prevent CrystalNets from recognizing any
topology before they are explicitly added.

The `refresh` optional parameter controls whether the current archive should also
be emptied.

!!! warning
    This empty archive will be kept and used for subsequent runs of CrystalNets.jl, even
    if you restart your Julia session. If you only want to empty the current archive,
    do `empty!(CrystalNets.CRYSTAL_NETS_ARCHIVE)`.
"""
function empty_default_archive!(; refresh=true)
    global CRYSTAL_NETS_ARCHIVE
    global REVERSE_CRYSTAL_NETS_ARCHIVE
    export_arc(arc_location, true)
    if refresh
        empty!(CRYSTAL_NETS_ARCHIVE)
        empty!(REVERSE_CRYSTAL_NETS_ARCHIVE)
    end
    nothing
end

function _change_current_archive!(newarc)
    global CRYSTAL_NETS_ARCHIVE
    global REVERSE_CRYSTAL_NETS_ARCHIVE
    empty!(CRYSTAL_NETS_ARCHIVE)
    empty!(REVERSE_CRYSTAL_NETS_ARCHIVE)
    merge!(CRYSTAL_NETS_ARCHIVE, newarc)
    merge!(REVERSE_CRYSTAL_NETS_ARCHIVE,
           Dict{String,String}(last(x) => first(x) for x in CRYSTAL_NETS_ARCHIVE))
    nothing
end

"""
    change_current_archive!(custom_arc; validate=true)

Erase the current archive used by CrystalNets.jl to recognize known topologies and
replace it with the archive stored in the file located at `custom_arc`.

The `validate` optional parameter controls whether the new file is checked and converted
to a format usable by CrystalNets.jl. If unsure, leave it set.

!!! note
    This modification will only last for the duration of this Julia session.

    If you wish to change the default archive and use it for subsequent runs, use
    [`clean_default_archive!`](@ref).

!!! warning
    Using an invalid archive will make CrystalNets.jl unusable. If this happens,
    simply run [`refresh_current_archive!()`](@ref) to revert to the
    default archive.
"""
function change_current_archive!(custom_arc; validate=true)
    arc::Dict{String,String} = if validate
         validate_archive(custom_arc)
    else
        last(parse_arc(custom_arc))
    end
    _change_current_archive!(arc)
end

"""
    refresh_current_archive!()

Revert the current topological archive to the default one.
"""
function refresh_current_archive!()
    _change_current_archive!(last(parse_arcs(arc_location)))
end

function _update_archive!(id, genome)
    global CRYSTAL_NETS_ARCHIVE
    global REVERSE_CRYSTAL_NETS_ARCHIVE
    CRYSTAL_NETS_ARCHIVE[genome] = id
    REVERSE_CRYSTAL_NETS_ARCHIVE[id] = genome
    nothing
end

"""
    add_to_current_archive!(id, genome)

Mark `genome` as the topological genome associated with the name `id` in the
current archive.

The input `id` and `genome` are not modified by this operation.

!!! note
    This modification will only last for the duration of this Julia session.

    If you wish to save the archive and use it for subsequent runs, use
    [`set_default_archive!`](@ref) after calling this function.
"""
function add_to_current_archive!(id::AbstractString, genome::AbstractString)
    if !isnumeric(first(genome))
        throw(ArgumentError(lazy"""
            This genome ("$genome") does not look like a genome. Are you sure you did not mix `id` and `genome`?

            If you really want to associate this id with this genome, use `CrystalNets._update_archive!(id, genome)`
            """))
    end
    global CRYSTAL_NETS_ARCHIVE
    for (x,y) in CRYSTAL_NETS_ARCHIVE
        if x == genome
            y == id && return
            throw(ArgumentError(lazy"""
                This genome is already registered under the name "$y".

                If you really want to change the name associated with it, use `CrystalNets._update_archive!(id, genome)`
                """))
        end
        if y == id
            throw(ArgumentError(lazy"""
                The name $id already corresponds to a different genome: "$x"

                If you really want to store another genome with the same name, use `CrystalNets._update_archive!(id, genome)`
                """))
        end
    end
    _update_archive!(id, genome)
end

"""
    make_archive(path, destination=nothing)

Make an archive from the files located in the directory given by `path` and export
it to `destination`, if specified. Each file of the directory should correspond
to a unique topology: if a topology is encountered multiple times, it will be assigned
the name of the latest file that bore it.

The archive can then be used with [`change_current_archive!(destination; validate=false)`](@ref change_current_archive!)
for instance.
"""
function make_archive(path, destination, verbose=false)
    arc = Dict{String,String}()
    Threads.@threads for f in readdir(path; join=false)
        name = splitext(f)[1]
        verbose && print("Handling "*name*"... ")
        flag = false
        flagerror = Ref{Any}(Tuple{Vector{Int},String}[])
        genomes::Vector{Tuple{Vector{Int},String}} = try
            x = topological_genome(UnderlyingNets(parse_chemfile(path*f)))
            verbose && println(name*" done.")
            x
        catch e
            flag = true
            flagerror[] = e
            Tuple{Vector{Int},String}[]
        end
        for (i, (vmap, genome)) in enumerate(genomes)
            if startswith(genome, "unstable") || genome == "non-periodic"
                flag = true
                push!(flagerror[]::Vector{Vector{Int}}, (vmap, genome))
                continue
            end
            arc[genome] = length(genomes) == 1 ? name : (name * '_' * string(i))
        end
        if flag
            e = flagerror[]
            isinterrupt(e) && rethrow()
            if e isa Vector{Tuple{Vector{Int},String}}
                for (vmap, instability) in e
                    println(stderr, "The component of file ", f, " containing atoms ",
                                    vmap, " was found to be ", instability, '.')
                end
            else
                println(stderr, "Failed for file ", f, " with error:")
                showerror(stderr, e)
                Base.display_error(stderr, e, catch_backtrace())
                println(stderr)
            end
        end
    end
    if !(destination isa Nothing)
        export_arc(destination, false, arc)
    end
    return arc
end


text:
##### Beginning of file

import PredictMDExtra

function is_travis_ci(a::AbstractDict = ENV)
    result = ( lowercase(strip(get(a, "CI", ""))) == "true" ) &&
        ( lowercase(strip(get(a, "TRAVIS", ""))) == "true" ) &&
        ( lowercase(strip(get(a, "CONTINUOUS_INTEGRATION", "")) ) == "true")
    return result
end

if is_travis_ci()
    Test.@testset "Testing _print_list_of_package_imports()" begin
        PredictMDExtra._print_list_of_package_imports()
    end
end

##### End of file

ext: jl


text:
const agenda_url_prefix = "http://somervillecityma.iqm2.com/Citizens/Detail_Meeting.aspx?ID="

const agenda_version = 1
const agenda_schema = Legolas.Schema("agenda", agenda_version)
const Agenda = Legolas.@row("agenda@1", item::Union{Int,Nothing,Missing}, content::String,
                            is_heading::Bool = isnothing(item))

function agenda_cache_path(cache_dir, id)
    return joinpath(cache_dir, "v$(agenda_version)", "agenda-$(id).arrow")
end

"""
    get_agenda_items(meeting_link; cache_dir=nothing)

Return `DataFrame` of agenda items from meeting `meeting_link`. `meeting_link` can be either
the full link (`"http://somervillecityma.iqm2.com/Citizens/Detail_Meeting.aspx?ID=2570"`) or
the meeting ID (`2570`).

When `cache_dir` is not empty, agendas are read from the cache if they've been previously
downloaded. If they do not exist in the cache, they'll be added to it.
"""
function get_agenda_items(meeting_link; cache_dir=nothing)
    isnothing(cache_dir) && return request_agenda_items(meeting_link)

    # Figure out caching info
    id = replace(meeting_link, agenda_url_prefix => "")
    cache_path = agenda_cache_path(cache_dir, id)
    if isnothing(tryparse(Int, id))
        @warn "Unexpected agenda meeting agenda id; may yield wonky cache path!" meeting_link id cache_path
    end

    # Has it been cached before?
    if !isfile(cache_path)
        items = request_agenda_items(meeting_link)
        mkpath(dirname(cache_path))
        Legolas.write(cache_path, items, agenda_schema)
    end
    return DataFrame(Legolas.read(cache_path))
end

"""
    request_agenda_items(meeting_id; verbose=false)

Return `DataFrame` of agenda items from meeting `meeting_link`. `meeting_link` can be either
the full link (`"http://somervillecityma.iqm2.com/Citizens/Detail_Meeting.aspx?ID=2570"`) or
the meeting ID (`2570`).
"""
function request_agenda_items(meeting_id; verbose=false)
    return request_agenda_items(agenda_url_prefix * string(meeting_id); verbose)
end

function request_agenda_items(meeting_link::AbstractString; verbose=false)
    if !startswith(meeting_link, agenda_url_prefix)
        meeting_link = agenda_url_prefix * meeting_link
        @warn "Inserting missing agenda prefix" meeting_link
    end

    r = HTTP.get(meeting_link)
    html = with_logger(NullLogger()) do
        return root(parsehtml(String(r.body)))
    end

    rows = findall("//td[@class='Title']", html)
    items = map(rows) do r
        length(elements(r)) == 0 && return nothing
        return only(elements(r))
    end
    filter!(!isnothing, items)

    # Headings have `name` "strong"; main items are numbered ("number : description").
    # Remove items that are not headings or top-level items.
    filter!(i -> i.name == "strong" || contains(i.content, ":"), items)

    items = map(items) do i
        i.name == "strong" && return Agenda(; item=nothing, i.content)
        x = split(i.content, " : "; limit=2)
        item = tryparse(Int, x[1])
        content = length(x) == 2 && !isnothing(item) ? x[2] : x[1]
        return Agenda(; item, content)
    end
    return DataFrame(items)
end

"""
    filter_agenda(agenda::DataFrame, search_terms; case_invariant=true)

Return copy of `agenda` filtered such that only rows containing an item from
`search_terms` in the description of the item will be preserved. When `case_invariant`
is true, ignores uppercase vs lowercase letters.
"""
function filter_agenda(agenda::DataFrame, search_terms; case_invariant=true)
    nrow(agenda) == 0 && return agenda
    case_invariant && (search_terms = lowercase.(search_terms))
    @debug search_terms

    check_item = (is_heading, content) -> begin
        is_heading && return false
        case_invariant && (content = lowercase(content))
        return any(occursin(p, content) for p in search_terms)
    end
    return filter([:is_heading, :content] => check_item, agenda)
end

#####
## All together now
#####

"""
    search_agendas_for_content(start_date, stop_date, search_terms;
                               cache_dir=nothing, case_invariant=true)

Return `NamedTuple` containing `meetings` (a `DataFrame` of meetings found within the given `start_date`-`stop_date`
range) and `items` (a `DataFrame` of all items containing at least one of the `search_terms`).
"""
function search_agendas_for_content(start_date, stop_date, search_terms; cache_dir=nothing,
                                    case_invariant=true)
    meetings = request_meetings(start_date, stop_date)
    @info """Found agendas for $(nrow(meetings)) meetings between $(start_date) and $(stop_date)!
             Searching their agendas for $(search_terms)..."""

    relevant_items = DataFrame()
    p = Progress(nrow(meetings) + 1)
    next!(p)
    _get_relevant_items = (link, date, name) -> begin
        try
            next!(p)
            total_items = 0
            items = get_agenda_items(link; cache_dir)
            total_items = nrow(items)
            items = filter_agenda(items, search_terms; case_invariant)
            nrow(items) == 0 && return (0, total_items, nothing)

            insertcols!(items, :meeting => name, :date => date, :meeting_link => link)
            relevant_items = vcat(relevant_items, items)
            return (nrow(items), total_items, nothing)
        catch e
            return (0, total_items, e)
        end

        return true
    end
    transform!(meetings,
               [:link, :date, :name] => ByRow(_get_relevant_items) => [:num_items,
                                                                       :total_items,
                                                                       :failed_parsing])

    # Summarize
    num_rel_m = nrow(relevant_items) == 0 ? 0 : length(unique(relevant_items.meeting_link))
    num_failed = nrow(meetings) - count(isnothing.(meetings.failed_parsing))
    num_irrel = nrow(meetings) - num_failed - num_rel_m
    failed = num_failed == 0 ? "" :
             "\n  -> $(num_failed) meetings that failed parsing (may be relevant)"
    @info """For the $(nrow(meetings)) meetings between $(start_date) and $(stop_date):
            -> $(length(unique(relevant_items.meeting_link))) meeting(s) with a total of $(nrow(relevant_items)) relevant item(s)
            -> $(num_irrel) meeting(s) with no relevant items$failed
          """
    return (; items=relevant_items, meetings)
end

"""
    display_items_by_meeting(items::DataFrame)

Display `items` grouped by meeting, along with the meeting's date, name, and url.
"""
function display_items_by_meeting(items::DataFrame)
    gdf = groupby(items, :meeting_link)
    for g in gdf
        meeting = only(unique(g.meeting))
        date = only(unique(g.date))
        link = only(unique(g.meeting_link))
        println(" * $date - $meeting: $link")
        for item in g.content
            println("   -> ", item)
        end
        print("\n")
    end
    return nothing
end


text:
module MCTS

using POMDPs
using POMDPModelTools
using POMDPPolicies
using POMDPSimulators
#using Blink
using CPUTime
using Random
using Printf

export
    MCTSSolver,
    MCTSPlanner,
    DPWSolver,
    DPWPlanner,
    BeliefMCTSSolver,
    AbstractMCTSPlanner,
    AbstractMCTSSolver,
    solve,
    action,
    action_info,
    rollout,
    StateNode,
    RandomActionGenerator,
    RolloutEstimator,
    next_action,
    clear_tree!,
    estimate_value,
    init_N,
    init_Q,
    children,
    n_children,
    isroot,
    default_action,
    get_state_node

export
    AbstractStateNode,
    StateActionStateNode,
    DPWStateActionNode,
    DPWStateNode,

    ExceptionRethrow,
    ReportWhenUsed

abstract type AbstractMCTSPlanner{P<:Union{MDP,POMDP}} <: Policy end
abstract type AbstractMCTSSolver <: Solver end
abstract type AbstractStateNode end

include("requirements_info.jl")
include("domain_knowledge.jl")
include("vanilla.jl")
include("dpw_types.jl")
include("dpw.jl")
include("action_gen.jl")
include("util.jl")
include("default_action.jl")
include("belief_mcts.jl")

include("visualization.jl")

end # module


text:
__precompile__(true)

module Microeconometrics

using DataFrames
using Dates
using Format: format
using LinearAlgebra
using Optim
using SparseArrays
using SpecialFunctions: lgamma
using StatsFuns
using StatsModels
    const SM = StatsModels
using SuiteSparse

using StatsBase: AbstractWeights, CoefTable, RegressionModel
using StatsBase: AnalyticWeights, FrequencyWeights, ProbabilityWeights, Weights
using StatsBase: mean, sum, pweights

import Base:        copy, isequal, show, sum
import Statistics:  mean
import StatsBase:   fit, coef, coefnames, coeftable, confint, stderror, vcov
import StatsBase:   deviance, loglikelihood, nulldeviance, nullloglikelihood
import StatsBase:   adjr2, dof, dof_residual, nobs, r2
import StatsBase:   fitted, predict, residuals, response
import StatsModels: terms, termvars, schema, apply_schema, has_schema, modelcols

include("./inference/corr.jl")
include("./data/types.jl")
include("./data/weights.jl")
include("./general/types.jl")
include("./data/utils.jl")
include("./general/etable.jl")
include("./general/statsmodels.jl")
include("./general/utils.jl")
include("./inference/adjfactor.jl")
include("./inference/hausman.jl")
include("./inference/utils.jl")
include("./inference/vcov.jl")
include("./estimation/ols.jl")
include("./estimation/iv.jl")
include("./estimation/logit.jl")
include("./estimation/probit.jl")
include("./estimation/cloglog.jl")
include("./estimation/gompit.jl")
include("./estimation/poisson.jl")
include("./estimation/ivpoisson.jl")
include("./estimation/mullahy.jl")
include("./estimation/ipw.jl")
include("./estimation/abadie.jl")
include("./estimation/frölichmelly.jl")
include("./estimation/tan.jl")

export
    @micromodel, Microdata,
    Homoscedastic, Heteroscedastic, Clustered, CrossCorrelated,
    OLS, Logit, Probit, Cloglog, Gompit, Poisson,
    IV, IVPoisson, Mullahy,
    IPW, Abadie, FrölichMelly, Tan

export
    fit, coef, confint, stderror, vcov,
    deviance, loglikelihood, nulldeviance, nullloglikelihood,
    adjr2, dof, dof_residual, nobs, r2,
    fitted, predict, residuals, response,
    coefnames, coeftable, etable,
    hausman_1s, hausman_2s, pval, tstat

end # module


text:
type TPQRCompactWY{S}
  R::StridedMatrix{S}
  V::StridedMatrix{S}
  T::Array{Float64,2}
end

function tpqrt!( A, B, l, T, work )
  m, n = size(B)
  nb = size(T,1)
  info = 0
  ccall((:dtpqrt_,Base.libblas_name), Void, (Ptr{Int64}, Ptr{Int64},
    Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Float64},
    Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}), &m,
    &n, &l, &nb, A, &max(1, stride(A,2)), B, &max(1, stride(B,2)), T, &max(1,
    stride(T,2)), work, &info)
  TPQRCompactWY(A, B, T)
end

tpqrt!(A,B,l,T) = tpqrt!(A,B,l,T, Array(Float64,size(T,1)*size(A,2)))
tpqrt!(A,B,l) = tpqrt!(A,B,l, Array(Float64,min(36,size(A,2)),size(A,2)))
tpqrt!(A,B) = tpqrt!(A,B,0.0)

function tpmqrt!( Q::TPQRCompactWY, A )
  m, n = size(A)
  k = l = n
  nb = min( n, 36 )
  B = zeros(m, n)
  work = Array(Float64, nb*n)
  info = Array(Int64, 1)
  trans = 'N'
  side = 'L'
  ccall((:dtpmqrt_,Base.libblas_name), Void, (Ptr{Uint8}, Ptr{Uint8},Ptr{Int64},
  Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64},
  Ptr{Float64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64},
  Ptr{Float64}, Ptr{Int64}),
  &side, &trans, &m, &n, &k, &l, &nb, Q.V, &max(1, stride(Q.V,2)), Q.T,
  &max(1, stride(Q.T,2)), A, &max(1, stride(A,2)), B, &max(1, stride(B,2)),
  work, info)
  # Returns Q1 and Q2, where QA = [Q1; Q2]
  A, B
end

function mylarfg!( n::Int64, alpha::Ptr{Float64}, X::Ptr{Float64}, incx::Int64, tau::Ptr{Float64} )
  ccall((:dlarfg_, Base.libblas_name), Void, (Ptr{Int64}, Ptr{Float64},
    Ptr{Float64}, Ptr{Int64}, Ptr{Float64}), &n, alpha, X, &incx, tau)
end

function mygemv!(trans::Char, m::Int64, n::Int64, alpha::Float64,
  A::Ptr{Float64}, lda::Int64,
  X::Ptr{Float64}, incx::Int64, beta::Float64, Y::Ptr{Float64}, incy::Int64)
  ccall((:dgemv_,Base.libblas_name), Void,
    (Ptr{UInt8}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64},
     Ptr{Float64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64},
     Ptr{Float64}, Ptr{Float64}, Ptr{Int64}),
    &trans, &m, &n, &alpha,
    A, &lda, X, &incx,
    &beta, Y, &incy)
end

function myger!(m::Int64, n::Int64, α::Float64, x::Ptr{Float64}, incx::Int64,
              y::Ptr{Float64}, incy::Int64, A::Ptr{Float64}, lda::Int64)
  ccall((:dger_, Base.libblas_name), Void,
    (Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Float64},
     Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Float64},
     Ptr{Int64}),
    &m, &n, &α, x,
    &incx, y, &incy, A,
    &lda)
end

function mytrmm!(side::Char, uplo::Char, trans::Char, diag::Char,
                m::Int64,n::Int64, alpha::Float64, A::Ptr{Float64}, lda::Int64,
                B::Ptr{Float64}, ldb::Int64)
  ccall((:dtrmm_, Base.libblas_name), Void,
    (Ptr{Char}, Ptr{Char}, Ptr{Char}, Ptr{Char}, Ptr{Int64}, Ptr{Int64},
      Ptr{Float64}, Ptr{Float64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}),
    &side, &uplo, &trans, &diag,
    &m, &n, &alpha, A, &lda, B, &ldb)
end

@inline function loc(A::Ptr{Float64}, lda::Int64, i::Int64, j::Int64)
  return A + sizeof(eltype(A))*(lda*(j-1)+i-1)
end

@inline function loc(A::Ptr{Float64}, i::Int64)
  return A + sizeof(eltype(A))*(i-1)
end

# QR factorization of triangle on top of rectangle
function tpqrf!( m::Int64, n::Int64, A::Ptr{Float64}, lda::Int64,
                 B::Ptr{Float64}, ldb::Int64, tau::Ptr{Float64}, work::Ptr{Float64} )
  for j = 1:n
    n1 = n-j
    # Construct Householder reflector
    mylarfg!(m+1, loc(A, lda, j, j), loc(B, ldb, 1, j), 1, tau)
    # Apply reflector to remainder
    BLAS.blascopy!(n1, loc(A, lda, j, j+1), lda, loc(work, j+1), 1)
    mygemv!('T', m, n1, 1.0, loc(B, ldb, 1, j+1), ldb, loc(B, ldb, 1, j), 1, 1.0,
      loc(work, j+1), 1)
    BLAS.axpy!(n1, -unsafe_load(tau), loc(work, j+1), 1,
      loc(A, lda, j, j+1), lda)
    myger!(m, n1, -unsafe_load(tau), loc(B, ldb, 1, j), 1, loc(work, j+1), 1, loc(B, ldb, 1, j+1), ldb)
  end
end

function tpmqrf!( m::Int64, n::Int64, A::Ptr{Float64}, lda::Int64,
                  B::Ptr{Float64}, ldb::Int64, work::Ptr{Float64})
  for j = n:-1:1
    n1 = n-j
    # Compute tau from the stored Householder vector
    t = 2/(1.0+BLAS.dot(m, loc(B, ldb, 1, j), 1, loc(B, ldb, 1, j), 1))
    a = unsafe_load(loc(A, lda, j, j))

    # Apply reflector to submatrix. Confusing pointers!
    BLAS.blascopy!(n1+1, loc(A, lda, j, j), lda, loc(work, j), 1)
    mygemv!('T', m, n1, 1.0, loc(B, ldb, 1, j+1), ldb, loc(B, ldb, 1, j), 1, 1.0,
      loc(work, j+1), 1)
    BLAS.axpy!(n1+1, -t, loc(work, j), 1, loc(A, lda, j, j), lda)
    # unsafe_store!(loc(A, lda, j, j), a-t*a)
    myger!(m, n1, -t, loc(B, ldb, 1, j), 1, loc(work, j+1), 1, loc(B, ldb, 1, j+1), ldb)
    w = unsafe_load(loc(work,j))
    BLAS.scal!(m, -t*w, loc(B,ldb,1,j), 1)
  end
end

# Compute Q'A for A split into nxk A and mxk B
function tpmrqrf!( m::Int64, n::Int64, k::Int64, A::Ptr{Float64}, lda::Int64,
                   B::Ptr{Float64}, ldb::Int64,
                   V::Ptr{Float64}, ldv::Int64, work::Ptr{Float64} )
  for j = 1:n
    # Compute tau
    t = 2/(1.0+BLAS.dot(m, loc(V, ldv, 1, j), 1, loc(V, ldv, 1, j), 1))
    BLAS.blascopy!( k, loc(A,j), lda, work, 1 )
    mygemv!('T', m, k, 1.0, loc(B,1), ldb, loc(V,ldv,1,j), 1, 1.0, work, 1)
    BLAS.axpy!(k, -t, work, 1, loc(A,j), lda)
    myger!(m,k,-t, loc(V,ldv,1,j), 1, work, 1, loc(B,1), ldb)
  end
end

A_mul_B!(A::TPQRCompactWY, B::StridedVecOrMat) = tpmqrt!(A,B)

# Sparse matrix multiple vector multiplication
function mkl_cscmm!( Y::Array, A::SparseMatrixCSC, X::Array )
  m,n = size(X)
  trans = 'N'
  alpha = 1.0
  beta = 0.0
  ccall((:mkl_dcscmm_, Base.libblas_name), Void,
  (Ptr{UInt8}, Ptr{Int64}, Ptr{Int64},
  Ptr{Int64}, Ptr{Float64}, Ptr{UInt8}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64},
  Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Float64}, Ptr{Float64},
  Ptr{Int64}),
  &trans, &m, &n, &m, &alpha, "GUNF", A.nzval, A.rowval, A.colptr,
  pointer(A.colptr,2), X, &stride(X,2), &beta, Y,
  &stride(Y,2))
end

# Sparse symmetric matrix multiple vector multiplication (only store upper
# triangle of A!)
function mkl_cscsymm!( Y::Array, A::SparseMatrixCSC, X::Array )
  m,n = size(X)
  trans = 'N'
  alpha = 1.0
  beta = 0.0
  ccall((:mkl_dcscmm_, Base.libblas_name), Void,
  (Ptr{UInt8}, Ptr{Int64}, Ptr{Int64},
  Ptr{Int64}, Ptr{Float64}, Ptr{UInt8}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64},
  Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Float64}, Ptr{Float64},
  Ptr{Int64}),
  &trans, &m, &n, &m, &alpha, "SUNF", A.nzval, A.rowval, A.colptr,
  pointer(A.colptr,2), X, &stride(X,2), &beta, Y,
  &stride(Y,2))
end

function dlarfbtr!( V::StridedMatrix, T::StridedMatrix, R::StridedMatrix,
  work::StridedMatrix )
  m,n = size(R)
  k = size(T,1)
  for j = 1:k
    @simd for i = 1:k
      @inbounds work[i,j] = R[i,j]
    end
  end
  for j = 1:k
    @simd for i = k+1:m
      @inbounds R[i,j] = V[i,j]
    end
  end

  BLAS.trmm!( 'L', 'L', 'T', 'U', 1.0, sub(V, 1:k, 1:k), work )
  BLAS.trmm!( 'L', 'U', 'N', 'N', 1.0, T, work )
  BLAS.trmm!( 'R', 'U', 'N', 'N', -1.0, work, sub(R,k+1:m,1:k))
  BLAS.trmm!( 'L', 'L', 'N', 'U', 1.0, sub(V,1:k,1:k), work )
  for j = 1:k
    @simd for i = 1:k
      @inbounds R[i,j] = R[i,j] - work[i,j]
    end
  end
end

function dlarfb!( V::StridedMatrix, T::StridedMatrix, C::StridedMatrix,
  work::StridedMatrix )
  m,n = size(C)
  k = size(T,1)
  side = 'L'
  trans = 'N'
  direct = 'F'
  storev = 'C'
  ccall((:dlarfb_, Base.libblas_name), Void,
  (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8},
   Ptr{Int64}, Ptr{Int64}, Ptr{Int64},
   Ptr{Float64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64},
   Ptr{Float64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}),
  &side, &trans, &direct, &storev, &m, &n, &k, V, &stride(V,2),
  T, &stride(T,2), C, &stride(C,2), work, &stride(work,2))
end

function gemqrtr!( F::Base.LinAlg.QRCompactWY, R::Array{Float64,2} )
  m,n = size(R)
  nb = 36
  work = zeros(n,nb)
  for i = (div(n-1,nb)*nb+1):-nb:1
    ib = min(n-i+1,nb)
    V = sub( F.factors, i:m, i:i+ib-1 )
    T = sub( F.T, 1:ib, i:i+ib-1 )
    dlarfbtr!( V, T, sub( R, i:m, i:i+ib-1 ), sub(work, 1:ib, 1:ib) )
    if i+ib-1 < n
      dlarfb!( V, T, sub( R, i:m, i+ib:n ), work )
    end
  end
end


text:
# to be called before the REPL is initialized

using VimBindings
# VimBindings.init()


text:
struct RectilinearGrid{FT, TX, TY, TZ, FX, FY, FZ, VX, VY, VZ, Arch} <: AbstractRectilinearGrid{FT, TX, TY, TZ, Arch}
    architecture :: Arch
    Nx :: Int
    Ny :: Int
    Nz :: Int
    Hx :: Int
    Hy :: Int
    Hz :: Int
    Lx :: FT
    Ly :: FT
    Lz :: FT
    # All directions can be either regular (FX, FY, FZ) <: Number
    # or stretched (FX, FY, FZ) <: AbstractVector
    Δxᶠᵃᵃ :: FX
    Δxᶜᵃᵃ :: FX
    xᶠᵃᵃ  :: VX
    xᶜᵃᵃ  :: VX
    Δyᵃᶠᵃ :: FY
    Δyᵃᶜᵃ :: FY
    yᵃᶠᵃ  :: VY
    yᵃᶜᵃ  :: VY
    Δzᵃᵃᶠ :: FZ
    Δzᵃᵃᶜ :: FZ
    zᵃᵃᶠ  :: VZ
    zᵃᵃᶜ  :: VZ

    function RectilinearGrid{TX, TY, TZ}(arch::Arch,
                                         Nx, Ny, Nz,
                                         Hx, Hy, Hz,
                                         Lx::FT, Ly::FT, Lz::FT,
                                         Δxᶠᵃᵃ :: FX, Δxᶜᵃᵃ :: FX,
                                          xᶠᵃᵃ :: VX,  xᶜᵃᵃ :: VX,
                                         Δyᵃᶠᵃ :: FY, Δyᵃᶜᵃ :: FY,
                                          yᵃᶠᵃ :: VY,  yᵃᶜᵃ :: VY,
                                         Δzᵃᵃᶠ :: FZ, Δzᵃᵃᶜ :: FZ,
                                          zᵃᵃᶠ :: VZ,  zᵃᵃᶜ :: VZ) where {Arch, FT,
                                                                          TX, TY, TZ,
                                                                          FX, VX, FY,
                                                                          VY, FZ, VZ}

        return new{FT, TX, TY, TZ, FX, FY, FZ, VX, VY, VZ, Arch}(arch, Nx, Ny, Nz,
                                                                 Hx, Hy, Hz, Lx, Ly, Lz,
                                                                 Δxᶠᵃᵃ, Δxᶜᵃᵃ, xᶠᵃᵃ, xᶜᵃᵃ,
                                                                 Δyᵃᶠᵃ, Δyᵃᶜᵃ, yᵃᶠᵃ, yᵃᶜᵃ,
                                                                 Δzᵃᵃᶠ, Δzᵃᵃᶜ, zᵃᵃᶠ, zᵃᵃᶜ)
    end
end

const XRegRectilinearGrid  = RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Number}
const YRegRectilinearGrid  = RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Any,    <:Number}
const ZRegRectilinearGrid  = RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Any,    <:Any,    <:Number}
const HRegRectilinearGrid  = RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Number, <:Number}
const XYRegRectilinearGrid = RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Number, <:Number}
const XZRegRectilinearGrid = RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Number, <:Any,    <:Number}
const YZRegRectilinearGrid = RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Any,    <:Number, <:Number}
const  RegRectilinearGrid  = RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Number, <:Number, <:Number}

regular_dimensions(::XRegRectilinearGrid)  = tuple(1)
regular_dimensions(::YRegRectilinearGrid)  = tuple(2)
regular_dimensions(::ZRegRectilinearGrid)  = tuple(3)
regular_dimensions(::XYRegRectilinearGrid) = (1, 2)
regular_dimensions(::XZRegRectilinearGrid) = (1, 3)
regular_dimensions(::YZRegRectilinearGrid) = (2, 3)
regular_dimensions(::RegRectilinearGrid)   = (1, 2, 3)

"""
    RectilinearGrid([architecture = CPU(), FT = Float64];
                    size,
                    x = nothing,
                    y = nothing,
                    z = nothing,
                    halo = nothing,
                    extent = nothing,
                    topology = (Periodic, Periodic, Bounded))

Creates a `RectilinearGrid` with `size = (Nx, Ny, Nz)` grid points.

Positional arguments
=================

- `architecture`: Specifies whether arrays of coordinates and spacings are stored
                  on the CPU or GPU. Default: `architecture = CPU()`.

- `FT` : Floating point data type. Default: `FT = Float64`.

Keyword arguments
=================

- `size` (required): A tuple prescribing the number of grid points in non-`Flat` directions.
                     `size` is a 3-tuple for 3D models, a 2-tuple for 2D models, and either a
                     scalar or 1-tuple for 1D models.

- `topology`: A 3-tuple `(TX, TY, TZ)` specifying the topology of the domain.
              `TX`, `TY`, and `TZ` specify whether the `x`-, `y`-, and `z` directions are
              `Periodic`, `Bounded`, or `Flat`. The topology `Flat` indicates that a model does
              not vary in those directions so that derivatives and interpolation are zero.
              The default is `topology = (Periodic, Periodic, Bounded)`.

- `extent`: A tuple prescribing the physical extent of the grid in non-`Flat` directions.
            All directions are contructed with regular grid spacing and the domain (in the
            case that no direction is `Flat`) is x ∈ (0, Lx), y ∈ (0, Ly), and z ∈ (-Lz, 0),
            which is most appropriate for oceanic applications with z = 0 usually being the
            ocean's surface.

- `x`, `y`, and `z`: Each of `x, y, z` are either (i) 2-tuples that specify the end points of the domain
                     in their respect directions (in which case scalar values may be used in `Flat`
                     directions), or (ii) arrays or functions of the corresponding indices `i`, `j`, or `k`
                     that specify the locations of cell faces in the `x`-, `y`-, or `z`-direction, respectively.
                     For example, to prescribe the cell faces in `z` we need to provide a function that takes
                     `k` as argument and retuns the location of the faces for indices `k = 1` through `k = Nz + 1`,
                     where `Nz` is the `size` of the stretched `z` dimension.

**Note**: _Either_ `extent`, or _all_ of `x`, `y`, and `z` must be specified.

- `halo`: A tuple of integers that specifies the size of the halo region of cells surrounding
          the physical interior for each non-`Flat` direction. The default is 3 halo cells in every direction.

The physical extent of the domain can be specified via `x`, `y`, and `z` keyword arguments
indicating the left and right endpoints of each dimensions, e.g. `x = (-π, π)` or via
the `extent` argument, e.g. `extent = (Lx, Ly, Lz)`, which specifies the extent of each dimension
in which case 0 ≤ x ≤ Lx, 0 ≤ y ≤ Ly, and -Lz ≤ z ≤ 0.

A grid topology may be specified via a tuple assigning one of `Periodic`, `Bounded`, and `Flat`
to each dimension. By default, a horizontally periodic grid topology `(Periodic, Periodic, Bounded)`
is assumed.

Constants are stored using floating point values of type `FT`. By default this is `Float64`.
Make sure to specify the desired `FT` if not using `Float64`.

Grid properties
===============

- `(Nx, Ny, Nz) :: Int`: Number of physical points in the ``(x, y, z)``-direction.

- `(Hx, Hy, Hz) :: Int`: Number of halo points in the ``(x, y, z)``-direction.

- `(Lx, Ly, Lz) :: FT`: Physical extent of the grid in the ``(x, y, z)``-direction.

- `(Δxᶜᵃᵃ, Δyᵃᶜᵃ, Δzᵃᵃᶜ)`: Grid spacing in the ``(x, y, z)``-direction between cell centers.
                           Defined at cell centers in ``x``, ``y``, and ``z``.

- `(Δxᶠᵃᵃ, Δyᵃᶠᵃ, Δzᵃᵃᶠ)`: Grid spacing in the ``(x, y, z)``-direction between cell faces.
                           Defined at cell faces in ``x``, ``y``, and ``z``.

- `(xᶜᵃᵃ, yᵃᶜᵃ, zᵃᵃᶜ)`: ``(x, y, z)`` coordinates of cell centers.

- `(xᶠᵃᵃ, yᵃᶠᵃ, zᵃᵃᶠ)`: ``(x, y, z)`` coordinates of cell faces.

Examples
========

* A default grid with `Float64` type:

```jldoctest
julia> using Oceananigans

julia> grid = RectilinearGrid(size=(32, 32, 32), extent=(1, 2, 3))
32×32×32 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0)  regularly spaced with Δx=0.03125
├── Periodic y ∈ [0.0, 2.0)  regularly spaced with Δy=0.0625
└── Bounded  z ∈ [-3.0, 0.0] regularly spaced with Δz=0.09375
```

* A default grid with `Float32` type:

```jldoctest
julia> using Oceananigans

julia> grid = RectilinearGrid(Float32; size=(32, 32, 16), x=(0, 8), y=(-10, 10), z=(-π, π))
32×32×16 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 8.0)          regularly spaced with Δx=0.25
├── Periodic y ∈ [-10.0, 10.0)       regularly spaced with Δy=0.625
└── Bounded  z ∈ [-3.14159, 3.14159] regularly spaced with Δz=0.392699
```

* A two-dimenisional, horizontally-periodic grid:

```jldoctest
julia> using Oceananigans

julia> grid = RectilinearGrid(size=(32, 32), extent=(2π, 4π), topology=(Periodic, Periodic, Flat))
32×32×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── Periodic x ∈ [3.60072e-17, 6.28319) regularly spaced with Δx=0.19635
├── Periodic y ∈ [7.20145e-17, 12.5664) regularly spaced with Δy=0.392699
└── Flat z
```

* A one-dimensional "column" grid:

```jldoctest
julia> using Oceananigans

julia> grid = RectilinearGrid(size=256, z=(-128, 0), topology=(Flat, Flat, Bounded))
1×1×256 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── Flat x
├── Flat y
└── Bounded  z ∈ [-128.0, 0.0]    regularly spaced with Δz=0.5
```

* A horizontally-periodic regular grid with cell interfaces stretched hyperbolically near the top:

```jldoctest
julia> using Oceananigans

julia> σ = 1.1; # stretching factor

julia> Nz = 24; # vertical resolution

julia> Lz = 32; # depth (m)

julia> hyperbolically_spaced_faces(k) = - Lz * (1 - tanh(σ * (k - 1) / Nz) / tanh(σ));

julia> grid = RectilinearGrid(size = (32, 32, Nz), x = (0, 64),
                              y = (0, 64), z = hyperbolically_spaced_faces)
32×32×24 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 64.0)   regularly spaced with Δx=2.0
├── Periodic y ∈ [0.0, 64.0)   regularly spaced with Δy=2.0
└── Bounded  z ∈ [-32.0, -0.0] variably spaced with min(Δz)=0.682695, max(Δz)=1.83091
```

* A three-dimensional grid with regular spacing in x,
  cell interfaces at Chebyshev nodes in y, and cell interfaces
  stretched in z hyperbolically near the top:

```jldoctest
julia> using Oceananigans

julia> Nx, Ny, Nz = 32, 30, 24;

julia> Lx, Ly, Lz = 200, 100, 32; # (m)

julia> chebychev_nodes(j) = - Ly/2 * cos(π * (j - 1) / Ny);

julia> σ = 1.1; # stretching factor

julia> hyperbolically_spaced_faces(k) = - Lz * (1 - tanh(σ * (k - 1) / Nz) / tanh(σ));

julia> grid = RectilinearGrid(size = (Nx, Ny, Nz),
                              topology=(Periodic, Bounded, Bounded),
                              x = (0, Lx),
                              y = chebychev_nodes,
                              z = hyperbolically_spaced_faces)
32×30×24 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 200.0)  regularly spaced with Δx=6.25
├── Bounded  y ∈ [-50.0, 50.0] variably spaced with min(Δy)=0.273905, max(Δy)=5.22642
└── Bounded  z ∈ [-32.0, -0.0] variably spaced with min(Δz)=0.682695, max(Δz)=1.83091
```
"""
function RectilinearGrid(architecture::AbstractArchitecture = CPU(),
                         FT = Float64;
                         size,
                         x = nothing,
                         y = nothing,
                         z = nothing,
                         halo = nothing,
                         extent = nothing,
                         topology = (Periodic, Periodic, Bounded))

    if architecture == GPU() && !has_cuda()
        throw(ArgumentError("Cannot create a GPU grid. No CUDA-enabled GPU was detected!"))
    end

    TX, TY, TZ, size, halo, x, y, z = validate_rectilinear_grid_args(topology, size, halo, FT, extent, x, y, z)

    Nx, Ny, Nz = size
    Hx, Hy, Hz = halo

    Lx, xᶠᵃᵃ, xᶜᵃᵃ, Δxᶠᵃᵃ, Δxᶜᵃᵃ = generate_coordinate(FT, topology[1], Nx, Hx, x, architecture)
    Ly, yᵃᶠᵃ, yᵃᶜᵃ, Δyᵃᶠᵃ, Δyᵃᶜᵃ = generate_coordinate(FT, topology[2], Ny, Hy, y, architecture)
    Lz, zᵃᵃᶠ, zᵃᵃᶜ, Δzᵃᵃᶠ, Δzᵃᵃᶜ = generate_coordinate(FT, topology[3], Nz, Hz, z, architecture)

    return RectilinearGrid{TX, TY, TZ}(architecture,
                                       Nx, Ny, Nz,
                                       Hx, Hy, Hz,
                                       FT(Lx), FT(Ly), FT(Lz),
                                       Δxᶠᵃᵃ, Δxᶜᵃᵃ, xᶠᵃᵃ, xᶜᵃᵃ,
                                       Δyᵃᶠᵃ, Δyᵃᶜᵃ, yᵃᶠᵃ, yᵃᶜᵃ,
                                       Δzᵃᵃᶠ, Δzᵃᵃᶜ, zᵃᵃᶠ, zᵃᵃᶜ)
end

""" Validate user input arguments to the `RectilinearGrid` constructor. """
function validate_rectilinear_grid_args(topology, size, halo, FT, extent, x, y, z)
    TX, TY, TZ = validate_topology(topology)
    size = validate_size(TX, TY, TZ, size)
    halo = validate_halo(TX, TY, TZ, halo)

    # Validate the rectilinear domain
    x, y, z = validate_rectilinear_domain(TX, TY, TZ, FT, size, extent, x, y, z)

    return TX, TY, TZ, size, halo, x, y, z
end

#####
##### Showing grids
#####

x_domain(grid::RectilinearGrid) = domain(topology(grid, 1), grid.Nx, grid.xᶠᵃᵃ)
y_domain(grid::RectilinearGrid) = domain(topology(grid, 2), grid.Ny, grid.yᵃᶠᵃ)
z_domain(grid::RectilinearGrid) = domain(topology(grid, 3), grid.Nz, grid.zᵃᵃᶠ)

# architecture = CPU() default, assuming that a DataType positional arg
# is specifying the floating point type.
RectilinearGrid(FT::DataType; kwargs...) = RectilinearGrid(CPU(), FT; kwargs...)

function Base.summary(grid::RectilinearGrid)
    FT = eltype(grid)
    TX, TY, TZ = topology(grid)

    return string(size_summary(size(grid)),
                  " RectilinearGrid{$FT, $TX, $TY, $TZ} on ", summary(architecture(grid)),
                  " with ", size_summary(halo_size(grid)), " halo")
end

function Base.show(io::IO, grid::RectilinearGrid, withsummary=true)
    TX, TY, TZ = topology(grid)

    x₁, x₂ = domain(topology(grid, 1), grid.Nx, grid.xᶠᵃᵃ)
    y₁, y₂ = domain(topology(grid, 2), grid.Ny, grid.yᵃᶠᵃ)
    z₁, z₂ = domain(topology(grid, 3), grid.Nz, grid.zᵃᵃᶠ)

    x_summary = domain_summary(TX(), "x", x₁, x₂)
    y_summary = domain_summary(TY(), "y", y₁, y₂)
    z_summary = domain_summary(TZ(), "z", z₁, z₂)

    longest = max(length(x_summary), length(y_summary), length(z_summary))

    x_summary = dimension_summary(TX(), "x", x₁, x₂, grid.Δxᶜᵃᵃ, longest - length(x_summary))
    y_summary = dimension_summary(TY(), "y", y₁, y₂, grid.Δyᵃᶜᵃ, longest - length(y_summary))
    z_summary = dimension_summary(TZ(), "z", z₁, z₂, grid.Δzᵃᵃᶜ, longest - length(z_summary))

    if withsummary
        print(io, summary(grid), '\n')
    end

    return print(io, "├── ", x_summary, '\n',
                     "├── ", y_summary, '\n',
                     "└── ", z_summary)
end

#####
##### Utilities
#####

function Adapt.adapt_structure(to, grid::RectilinearGrid)
    TX, TY, TZ = topology(grid)
    return RectilinearGrid{TX, TY, TZ}(nothing,
                                       grid.Nx, grid.Ny, grid.Nz,
                                       grid.Hx, grid.Hy, grid.Hz,
                                       grid.Lx, grid.Ly, grid.Lz,
                                       Adapt.adapt(to, grid.Δxᶠᵃᵃ),
                                       Adapt.adapt(to, grid.Δxᶜᵃᵃ),
                                       Adapt.adapt(to, grid.xᶠᵃᵃ),
                                       Adapt.adapt(to, grid.xᶜᵃᵃ),
                                       Adapt.adapt(to, grid.Δyᵃᶠᵃ),
                                       Adapt.adapt(to, grid.Δyᵃᶜᵃ),
                                       Adapt.adapt(to, grid.yᵃᶠᵃ),
                                       Adapt.adapt(to, grid.yᵃᶜᵃ),
                                       Adapt.adapt(to, grid.Δzᵃᵃᶠ),
                                       Adapt.adapt(to, grid.Δzᵃᵃᶜ),
                                       Adapt.adapt(to, grid.zᵃᵃᶠ),
                                       Adapt.adapt(to, grid.zᵃᵃᶜ))
end

@inline xnode(::Face  , i, grid::RectilinearGrid) = @inbounds grid.xᶠᵃᵃ[i]
@inline xnode(::Center, i, grid::RectilinearGrid) = @inbounds grid.xᶜᵃᵃ[i]
@inline ynode(::Face  , j, grid::RectilinearGrid) = @inbounds grid.yᵃᶠᵃ[j]
@inline ynode(::Center, j, grid::RectilinearGrid) = @inbounds grid.yᵃᶜᵃ[j]
@inline znode(::Face  , k, grid::RectilinearGrid) = @inbounds grid.zᵃᵃᶠ[k]
@inline znode(::Center, k, grid::RectilinearGrid) = @inbounds grid.zᵃᵃᶜ[k]

all_x_nodes(::Type{Face}  , grid::RectilinearGrid) = grid.xᶠᵃᵃ
all_x_nodes(::Type{Center}, grid::RectilinearGrid) = grid.xᶜᵃᵃ
all_y_nodes(::Type{Face}  , grid::RectilinearGrid) = grid.yᵃᶠᵃ
all_y_nodes(::Type{Center}, grid::RectilinearGrid) = grid.yᵃᶜᵃ
all_z_nodes(::Type{Face}  , grid::RectilinearGrid) = grid.zᵃᵃᶠ
all_z_nodes(::Type{Center}, grid::RectilinearGrid) = grid.zᵃᵃᶜ

cpu_face_constructor_x(grid::XRegRectilinearGrid) = x_domain(grid)
cpu_face_constructor_y(grid::YRegRectilinearGrid) = y_domain(grid)
cpu_face_constructor_z(grid::ZRegRectilinearGrid) = z_domain(grid)

function with_halo(new_halo, old_grid::RectilinearGrid)

    size = (old_grid.Nx, old_grid.Ny, old_grid.Nz)
    topo = topology(old_grid)

    x = cpu_face_constructor_x(old_grid)
    y = cpu_face_constructor_y(old_grid)
    z = cpu_face_constructor_z(old_grid)

    # Remove elements of size and new_halo in Flat directions as expected by grid
    # constructor
    size     = pop_flat_elements(size, topo)
    new_halo = pop_flat_elements(new_halo, topo)

    new_grid = RectilinearGrid(architecture(old_grid), eltype(old_grid);
                               size = size,
                               x = x, y = y, z = z,
                               topology = topo,
                               halo = new_halo)

    return new_grid
end

function on_architecture(new_arch::AbstractArchitecture, old_grid::RectilinearGrid)
    old_properties = (old_grid.Δxᶠᵃᵃ, old_grid.Δxᶜᵃᵃ, old_grid.xᶠᵃᵃ, old_grid.xᶜᵃᵃ,
                      old_grid.Δyᵃᶠᵃ, old_grid.Δyᵃᶜᵃ, old_grid.yᵃᶠᵃ, old_grid.yᵃᶜᵃ,
                      old_grid.Δzᵃᵃᶠ, old_grid.Δzᵃᵃᶜ, old_grid.zᵃᵃᶠ, old_grid.zᵃᵃᶜ)

    new_properties = Tuple(arch_array(new_arch, p) for p in old_properties)

    TX, TY, TZ = topology(old_grid)

    return RectilinearGrid{TX, TY, TZ}(new_arch,
                                       old_grid.Nx, old_grid.Ny, old_grid.Nz,
                                       old_grid.Hx, old_grid.Hy, old_grid.Hz,
                                       old_grid.Lx, old_grid.Ly, old_grid.Lz,
                                       new_properties...)
end

return_metrics(::RectilinearGrid) = (:xᶠᵃᵃ, :xᶜᵃᵃ, :yᵃᶠᵃ, :yᵃᶜᵃ, :zᵃᵃᶠ, :zᵃᵃᶜ)

#####
##### Get minima of grid
#####

function min_Δx(grid::RectilinearGrid)
    topo = topology(grid)
    if topo[1] == Flat
        return Inf
    else
        return min_number_or_array(grid.Δxᶜᵃᵃ)
    end
end

function min_Δy(grid::RectilinearGrid)
    topo = topology(grid)
    if topo[2] == Flat
        return Inf
    else
        return min_number_or_array(grid.Δyᵃᶜᵃ)
    end
end

function min_Δz(grid::RectilinearGrid)
    topo = topology(grid)
    if topo[3] == Flat
        return Inf
    else
        return min_number_or_array(grid.Δzᵃᵃᶜ)
    end
end

@inline min_number_or_array(var) = var
@inline min_number_or_array(var::AbstractVector) = minimum(parent(var))

isrectilinear(::RectilinearGrid) = true


text:
using GenomicTiles
using BioFeatures
using DataFrames
using Gadfly
using Clustering #need for assignments()
using Distances

function get_sub_regions(start, stop)
    subregions = DataFrame( start = [start, start + 50], stop = [start+60, start + 100],rtype=["A","B"])
    return Features( subregions )
end

#
# GET SOME COORDINATES TO SERVE AS A BASE THAT WILL BE SPLIT INTO TILES
#
region_lengths = int( rand(100)*1000 )
starts = 1:1000:100000
stops = starts + region_lengths

# Read these into a BioFeature - FeatureRegion
regions = DataFrame(start = starts, stop= stops)
FeaturesToBeSplitIntoTiles = Features( regions )

tiles = nothing
tiles_matrix = nothing
for features = values( FeaturesToBeSplitIntoTiles.features )
    for feature in values(features)
         subregions = get_sub_regions( feature[:start][1],feature[:stop][1] )
         tiles = make_tiles(feature, num_tiles=30, default_value = 0)
         #println(tiles)
         for tile in tiles
            dmap = @Compat.Dict("A"=> 3.0,"B"=-3.0,"MIXED"=0.0)
            tile = set_tile_value( tile, subregions, :rtype,dmap)
         end
         tiles_matrix = add_tiles_to_tiles_matrix(tiles_matrix, tiles)
    end
end

tiles_matrix[1,:]

(tiles_matrix_clustered,kmeans_cluster_object,sil_df) = cluster_tiles(tiles_matrix',
                                                                      method="kmeans",
                                                                      num_clusters_range=5:10,
                                                                      kmeans_repeats=10)
path = "/home/nic/Desktop/test-kmeans.png"
colour_scale =Scale.discrete_color_manual("red","green")
plot_matrix2(tiles_matrix_clustered,path,colour_scale=colour_scale)
path = "/home/nic/Desktop/test-kmeans-bar.png"
plot_cluster_bar_chart(tiles_matrix',assignments(kmeans_cluster_object),path)

sil_df
Gadfly.plot(sil_df,x="r",y="ave_sil",color="k",Geom.line)

#tiles_matrix
#distance_matrix = pairwise(SqEuclidean(),tiles_matrix')
#length(assignments(kmeans_cluster_object))
#sil_val = silhouettes(kmeans_cluster_object,distance_matrix)
#mean(sil_val)
#=
(tiles_matrix_clustered,cluster_object,sil_val) = cluster_tiles(tiles_matrix, method="dbscan", eps=1,minpts=2)
path = "/home/nic/Desktop/test-dbscan.png"
colour_scale =Scale.discrete_color_manual("red","green")
plot_matrix2(tiles_matrix_clustered,path,colour_scale=colour_scale)
path = "/home/nic/Desktop/test-dbscan-bar.png"
plot_cluster_bar_chart(tiles_matrix, assignments( cluster_object), path)
=#

#cluster_object.assignments =cluster_object.assignments + 5
#(tiles_matrix_clustered,cluster_assignments) = cluster_tiles(tiles_matrix, method="affprop")
#path = "/home/nic/Desktop/test-affprop.png"
#plot_matrix2(tiles_matrix_clustered,path)
#path = "/home/nic/Desktop/test-affprop-bar.png"
#plot_cluster_bar_chart(tiles_matrix,cluster_assignments,path)
#cluster_assignments
#tiles_matrix'[1,:]
#tiles_matrix'[100,:]
#cluster_assignments


text:
function export_traffic_matrix(tm::ConsideredTrafficMatrix, basename="$(@__DIR__)")
    export_traffic_matrix(tm.matrix, "$(basename)_$(tm.subproblem_edge.src)_to_$(tm.subproblem_edge.dst)")
end

function export_traffic_matrix(tm::Dict{Edge{Int}, Float64}, filename::String)
    open(filename, "w") do f
        write(f, "Source, Destination, Traffic\n")
        for (edge, value) in tm
            write(f, "$(edge.src), $(edge.dst), $(value)\n")
        end
    end
end

"""
Summarises as text the differences between two different routings.
"""
function summary(rd::RoutingData, base_flows, oblivious_flows; filename="$(@__DIR__)/summary.txt", epsilon=1.e-2)
    # Helper functions.
    stringNode(n) = "$(get_prop(rd.g, n, :name))"
    stringEdge(rd::Tuple) = "($(stringNode(rd[1])), $(stringNode(rd[2])))"
    stringEdge(rd::Edge) = "($(stringNode(src(rd))), $(stringNode(dst(rd))))"
    stringPath(p) = join(map(stringEdge, p), ", ")

    # Get rid of JuMP's data structures.
    base_routing_matrix = zeros(n_demands(rd), n_paths(rd))
    for k in eachindex(base_flows)
        base_routing_matrix[k[1], k[2]] = base_flows[k...]
    end
    oblivious_routing_matrix = zeros(n_demands(rd), n_paths(rd))
    for k in eachindex(oblivious_flows)
        oblivious_routing_matrix[k[1], k[2]] = oblivious_flows[k...]
    end

    # Write the summary.
    list_edges = collect(edges(rd.g))
    open(filename, "w") do f
        write(f, "Total load\n")
        write(f, "==========\n\n")

        base_path_loads = sum(base_routing_matrix[d, :] for d in 1:n_demands(rd))
        oblivious_path_loads = sum(oblivious_routing_matrix[d, :] for d in 1:n_demands(rd))

        base_edge_throughputs = [any(in(e, rd.paths_edges[p]) for p in findall(base_path_loads .> 0)) ? sum(base_path_loads[p] * in(e, rd.paths_edges[p]) for p in findall(base_path_loads .> 0)) : 0.0 for e in 1:n_edges(rd)]
        oblivious_edge_throughputs = [any(in(e, rd.paths_edges[p]) for p in findall(oblivious_path_loads .> 0)) ? sum(oblivious_path_loads[p] * in(e, rd.paths_edges[p]) for p in findall(oblivious_path_loads .> 0)) : 0.0 for e in 1:n_edges(rd)]

        base_edge_loads = [base_edge_throughputs[i] / get_prop(rd.g, list_edges[i], :capacity) for i in 1:n_edges(rd)]
        oblivious_edge_loads = [oblivious_edge_throughputs[i] / get_prop(rd.g, list_edges[i], :capacity) for i in 1:n_edges(rd)]

        write(f, "Load    | Base    | Oblivious\n")
        write(f, "Maximum |  $(@sprintf("% 2.2f", 100 * maximum(base_edge_loads)))  | $(@sprintf("% 2.2f", 100 * maximum(oblivious_edge_loads)))\n")
        write(f, "Average |  $(@sprintf("% 2.2f", 100 * mean(base_edge_loads)))  | $(@sprintf("% 2.2f", 100 * mean(oblivious_edge_loads)))\n")
        write(f, "Minimum |  $(@sprintf("% 2.2f", 100 * minimum(base_edge_loads)))  | $(@sprintf("% 2.2f", 100 * minimum(oblivious_edge_loads)))\n")
        write(f, "\n")

        for i in 1:n_edges(rd)
            s = src(list_edges[i])
            t = dst(list_edges[i])

            write(f, "  - Edge #$(i): $(stringEdge((s, t)))\n") # TODO: Write a specialisation of stringEdge that directly works with edges?
            write(f, "    Base throughput: $(round(base_edge_throughputs[i], digits=2)) Mbps\n")
            write(f, "    Base load: $(round(100 * base_edge_loads[i], digits=2))%\n")
            write(f, "    Oblivious throughput: $(round(oblivious_edge_throughputs[i], digits=2)) Mbps\n")
            write(f, "    Oblivious load: $(round(100 * oblivious_edge_loads[i], digits=2))%\n")
        end

        write(f, "\n\n\n")

        for d in 1:n_demands(rd)
            k_d = rd.demand_id_to_demand[rd.path_id_to_demand_id[d]]

            write(f, "Demand #$(d): from $(src(k_d)) to $(dst(k_d))\n")
            write(f, "======\n\n")

            # When a route is used with a weight higher than epsilon, consider it is used (boolean).
            base_routes = vec(copy(base_routing_matrix[d, :]))
            base_routes[base_routes .>= epsilon] .= 1
            base_routes[base_routes .< epsilon] .= 0
            oblivious_routes = vec(copy(oblivious_routing_matrix[d, :]))
            oblivious_routes[oblivious_routes .>= epsilon] .= 1
            oblivious_routes[oblivious_routes .< epsilon] .= 0

            base_n_routes = sum(base_routes)
            oblivious_n_routes = sum(oblivious_routes)

            # Statistics on the number of possible paths for this demand (be they actually used or not).
            write(f, "Number of possible paths for this demand: $(length(rd.demand_id_to_path_ids[d]))\n")
            write(f, "Minimum path length (number of hops): $(minimum(map(p -> length(rd.paths[p]), rd.demand_id_to_path_ids[d])))\n")
            write(f, "Average path length (number of hops): $(round(mean(map(p -> length(rd.paths[p]), rd.demand_id_to_path_ids[d])), digits=2))\n")
            write(f, "Maximum path length (number of hops): $(maximum(map(p -> length(rd.paths[p]), rd.demand_id_to_path_ids[d])))\n")
            write(f, "\n")

            # Compare the number of routes.
            if base_n_routes != oblivious_n_routes
                write(f, "Base and oblivious solutions do not have the same number of routes: \n")
                write(f, "  - base: $(convert(Int, base_n_routes)) route$((base_n_routes != 1.0) ? "s" : "")\n")
                write(f, "  - oblivious: $(convert(Int, oblivious_n_routes)) route$((oblivious_n_routes != 1.0) ? "s" : "")\n")
            else
                write(f, "Base and oblivious solutions have the same number of routes. \n")
            end

            # Compare the actual routes.
            if any(base_routes .!= oblivious_routes)
                write(f, "Base and oblivious solutions do not have the same routes: \n")

                for r in 1:length(base_routes)
                    if base_routes[r] == 0 && oblivious_routes[r] == 1
                        write(f, "  - difference for route #$(r): used only by oblivious\n")
                        write(f, "    $(stringPath(rd.paths_edges[r]))\n")
                    elseif base_routes[r] == 1 && oblivious_routes[r] == 0
                        write(f, "  - difference for route #$(r): used only by base\n")
                        write(f, "    $(stringPath(rd.paths_edges[r]))\n")
                    end
                end
            else
                write(f, "Base and oblivious solutions have the same routes. \n")
            end

            # Compare the route weights.
            if norm(vec(base_routing_matrix[d, :]) - vec(oblivious_routing_matrix[d, :])) >= epsilon
                write(f, "Base and oblivious solutions do not have the same route weights: \n")

                for r in 1:nPaths
                    if abs(base_routing_matrix[d, r] - oblivious_routing_matrix[d, r]) >= epsilon
                        write(f, "  - significant difference for route #$(r): base $(round(base_routing_matrix[d, r], digits=2)), oblivious $(round(oblivious_routing_matrix[d, r], digits=2))\n")
                        write(f, "    $(stringPath(rd.paths_edges[r]))\n")
                    end
                end
            else
                write(f, "Base and oblivious solutions have the same route weights. \n")
            end

            write(f, "\n")

            # Loads of the links for each used path.
            write(f, "Loads for the base routing: \n")
            for p in rd.demand_id_to_path_ids[d]
                if base_routing_matrix[d, p] <= epsilon
                    continue
                end

                write(f, "  - Path #$(p): \n")

                for (i, e) in enumerate(rd.paths_edges[p])
                    write(f, "      - $(stringEdge(e)): $(base_edge_loads[i])\n")
                end
            end

            write(f, "Loads for the oblivious routing: \n")
            for p in rd.demand_id_to_path_ids[d]
                if oblivious_routing_matrix[d, p] <= epsilon
                    continue
                end

                write(f, "  - Path #$(p): \n")

                for (i, e) in enumerate(rd.paths_edges[p])
                    write(f, "      - $(stringEdge(e)): $(oblivious_edge_loads[i])\n")
                end
            end

            write(f, "\n\n\n")
        end
    end
end

function export_routing(routing::Routing, filename="$(@__DIR__)/oblivious_routing.txt"; epsilon=1.e-2)
    rd = routing.data
    nodeName(n) = get_prop(rd.g, n, :name)

    routing_matrix = routing_to_matrix(routing)
    open(filename, "w") do f
        write(f, "Demand number, Path proportion, Vertices (names), Edges (name origin - name destination)\n")
        for d in 1:n_demands(rd)
            for p in findall(x -> x > epsilon, routing_matrix[d, :])
                write(f, "$(d), $(routing_matrix[d, p]), ")

                # Path as a list of vertices.
                for edge in routing.paths[p]
                    write(f, "$(nodeName(src(edge))) ")
                end
                write(f, "$(nodeName(dst(routing.paths[p][end]))), ")

                # Path as a list of edges.
                for edge in routing.paths[p]
                    write(f, "($(nodeName(src(edge))) - $(nodeName(dst(edge))))")
                end

                write(f, "\n")
            end
        end
    end
end


text:
function compute()
	SIZE=1001
	ans=1
	ans+=sum(4*i*i-6*(i-1) for i in 3:2:SIZE)

	return ans
end

println(compute())


text:
# ---
# title: 1261. Find Elements in a Contaminated Binary Tree
# id: problem1261
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, Tree
# link: <https://leetcode.com/problems/find-elements-in-a-contaminated-binary-tree/description/>
# hidden: true
# ---
#
# Given a binary tree with the following rules:
#
#   1. `root.val == 0`
#   2. If `treeNode.val == x` and `treeNode.left != null`, then `treeNode.left.val == 2 * x + 1`
#   3. If `treeNode.val == x` and `treeNode.right != null`, then `treeNode.right.val == 2 * x + 2`
#
# Now the binary tree is contaminated, which means all `treeNode.val` have been
# changed to `-1`.
#
# You need to first recover the binary tree and then implement the
# `FindElements` class:
#
#   * `FindElements(TreeNode* root)` Initializes the object with a contamined binary tree, you need to recover it first.
#   * `bool find(int target)` Return if the `target` value exists in the recovered binary tree.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/11/06/untitled-
# diagram-4-1.jpg)**
#
#
#
#     **Input**
#     ["FindElements","find","find"]
#     [[[-1,null,-1]],[1],[2]]
#     **Output**
#     [null,false,true]
#     **Explanation**
#     FindElements findElements = new FindElements([-1,null,-1]);
#     findElements.find(1); // return False
#     findElements.find(2); // return True
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2019/11/06/untitled-diagram-4.jpg)**
#
#
#
#     **Input**
#     ["FindElements","find","find","find"]
#     [[[-1,-1,-1,-1,-1]],[1],[3],[5]]
#     **Output**
#     [null,true,true,false]
#     **Explanation**
#     FindElements findElements = new FindElements([-1,-1,-1,-1,-1]);
#     findElements.find(1); // return True
#     findElements.find(3); // return True
#     findElements.find(5); // return False
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2019/11/07/untitled-
# diagram-4-1-1.jpg)**
#
#
#
#     **Input**
#     ["FindElements","find","find","find","find"]
#     [[[-1,null,-1,-1,null,-1]],[2],[3],[4],[5]]
#     **Output**
#     [null,true,false,false,true]
#     **Explanation**
#     FindElements findElements = new FindElements([-1,null,-1,-1,null,-1]);
#     findElements.find(2); // return True
#     findElements.find(3); // return False
#     findElements.find(4); // return False
#     findElements.find(5); // return True
#
#
#
#
# **Constraints:**
#
#   * `TreeNode.val == -1`
#   * The height of the binary tree is less than or equal to `20`
#   * The total number of nodes is between `[1, 10^4]`
#   * Total calls of `find()` is between `[1, 10^4]`
#   * `0 <= target <= 10^6`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module InvisibleBarrier

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Invisible Barrier" => Ahorn.EntityPlacement(
        Maple.Barrier,
        "rectangle"
    ),
)

Ahorn.minimumSize(entity::Maple.Barrier) = 8, 8
Ahorn.resizable(entity::Maple.Barrier) =  true, true

Ahorn.selection(entity::Maple.Barrier) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Barrier, room::Maple.Room)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.4, 0.4, 0.4, 0.8), (0.0, 0.0, 0.0, 0.0))
end

end


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

abstract type IPAddr end

Base.isless(a::T, b::T) where {T<:IPAddr} = isless(a.host, b.host)
(dt::Type{<:Integer})(ip::IPAddr) = dt(ip.host)::dt

struct IPv4 <: IPAddr
    host::UInt32
    IPv4(host::UInt32) = new(host)
    IPv4(a::UInt8,b::UInt8,c::UInt8,d::UInt8) = new(UInt32(a)<<24|
                                                    UInt32(b)<<16|
                                                    UInt32(c)<<8|
                                                    d)
    function IPv4(a::Integer,b::Integer,c::Integer,d::Integer)
        if !(0<=a<=255 && 0<=b<=255 && 0<=c<=255 && 0<=d<=255)
            throw(ArgumentError("IPv4 field out of range (must be 0-255)"))
        end
        IPv4(UInt8(a),UInt8(b),UInt8(c),UInt8(d))
    end
end

"""
    IPv4(host::Integer) -> IPv4

Returns an IPv4 object from ip address `host` formatted as an [`Integer`](@ref).

```jldoctest
julia> IPv4(3223256218)
ip"192.30.252.154"
```
"""
function IPv4(host::Integer)
    if host < 0
        throw(ArgumentError("IPv4 address must be positive"))
    elseif typemax(typeof(host)) > typemax(UInt32) && host > typemax(UInt32)
        throw(ArgumentError("IPv4 address must fit within 32 bits"))
    else
        return IPv4(UInt32(host))
    end
end

# constructor: ("1.2.3.4")
IPv4(str::AbstractString) = parse(IPv4, str)

show(io::IO,ip::IPv4) = print(io,"ip\"",ip,"\"")
print(io::IO,ip::IPv4) = print(io,string((ip.host&(0xFF000000))>>24),".",
                                  string((ip.host&(0xFF0000))>>16),".",
                                  string((ip.host&(0xFF00))>>8),".",
                                  string(ip.host&0xFF))

struct IPv6 <: IPAddr
    host::UInt128
    IPv6(host::UInt128) = new(host)
    IPv6(a::UInt16,b::UInt16,c::UInt16,d::UInt16,
     e::UInt16,f::UInt16,g::UInt16,h::UInt16) = new(UInt128(a)<<(7*16)|
                            UInt128(b)<<(6*16)|
                            UInt128(c)<<(5*16)|
                            UInt128(d)<<(4*16)|
                            UInt128(e)<<(3*16)|
                            UInt128(f)<<(2*16)|
                            UInt128(g)<<(1*16)|
                            h)
    function IPv6(a::Integer,b::Integer,c::Integer,d::Integer,
          e::Integer,f::Integer,g::Integer,h::Integer)
        if !(0<=a<=0xFFFF && 0<=b<=0xFFFF && 0<=c<=0xFFFF && 0<=d<=0xFFFF &&
             0<=e<=0xFFFF && 0<=f<=0xFFFF && 0<=g<=0xFFFF && 0<=h<=0xFFFF)
            throw(ArgumentError("IPv6 field out of range (must be 0-65535)"))
        end
        IPv6(UInt16(a),UInt16(b),UInt16(c),UInt16(d),
             UInt16(e),UInt16(f),UInt16(g),UInt16(h))
    end
end

"""
    IPv6(host::Integer) -> IPv6

Returns an IPv6 object from ip address `host` formatted as an [`Integer`](@ref).

```jldoctest
julia> IPv6(3223256218)
ip"::c01e:fc9a"
```
"""
function IPv6(host::Integer)
    if host < 0
        throw(ArgumentError("IPv6 address must be positive"))
        # We allow passing bigger integer types, but need to be careful to avoid overflow
        # Let's hope promotion rules are sensible
    elseif typemax(typeof(host)) > typemax(UInt128) && host > typemax(UInt128)
        throw(ArgumentError("IPv6 address must fit within 128 bits"))
    else
        return IPv6(UInt128(host))
    end
end

IPv6(str::AbstractString) = parse(IPv6, str)

# Suppress leading '0's and "0x"
print_ipv6_field(io,field::UInt16) = print(io,string(field, base = 16))

print_ipv6_field(io,ip,i) = print_ipv6_field(io,ipv6_field(ip,i))
function ipv6_field(ip::IPv6,i)
    if i < 0 || i > 7
        throw(BoundsError())
    end
    UInt16(ip.host&(UInt128(0xFFFF)<<(i*16))>>(i*16))
end

show(io::IO, ip::IPv6) = print(io,"ip\"",ip,"\"")
# RFC 5952 compliant show function
# http://tools.ietf.org/html/rfc5952
function print(io::IO,ip::IPv6)
    i = 8
    m = 0
    longest_sub_i = -1
    while i!=0
        i-=1
        field = ipv6_field(ip,i)
        if field == 0 && longest_sub_i == -1
            # Find longest subsequence of 0
            longest_sub_i,j,m,c = i,i,1,1
            while j != 0
                j-=1
                if ipv6_field(ip,j) == 0
                    c += 1
                else
                    c = 0
                end
                if c > m
                    if j+c != longest_sub_i+1
                        longest_sub_i = j+c-1
                    end
                    m = c
                end
            end
            # Prevent single 0 from contracting to :: as required
            if m == 1
                longest_sub_i = 9
            end
        end
        if i == longest_sub_i
            print(io,":")
            i -= m-1
            if i == 0
                print(io,":")
                break
            end
        else
            if i != 7
                print(io,":")
            end
            print_ipv6_field(io,field)
        end
    end
end

# Parsing

const ipv4_leading_zero_error = """
Leading zeros in IPv4 addresses are disallowed due to ambiguity.
If the address is in octal or hexadecimal, convert it to decimal, otherwise remove the leading zero.
"""

function parse(::Type{IPv4}, str::AbstractString)
    fields = split(str,'.')
    i = 1
    ret = 0
    for f in fields
        if isempty(f)
            throw(ArgumentError("empty field in IPv4 address"))
        end
        if length(f) > 1 && f[1] == '0'
            throw(ArgumentError(ipv4_leading_zero_error))
        else
            r = parse(Int, f, base = 10)
        end
        if i != length(fields)
            if r < 0 || r > 255
                throw(ArgumentError("IPv4 field out of range (must be 0-255)"))
            end
            ret |= UInt32(r) << ((4-i)*8)
        else
            if r > ((UInt64(1)<<((5-length(fields))*8))-1)
                throw(ArgumentError("IPv4 field too large"))
            end
            ret |= r
        end
        i+=1
    end
    IPv4(ret)
end

function parseipv6fields(fields,num_fields)
    if length(fields) > num_fields
        throw(ArgumentError("too many fields in IPv6 address"))
    end
    cf = 7
    ret = UInt128(0)
    for f in fields
        if isempty(f)
            # ::abc:... and ..:abc::
            if cf != 7 && cf != 0
                cf -= num_fields-length(fields)
            end
            cf -= 1
            continue
        end
        ret |= UInt128(parse(Int, f, base = 16))<<(cf*16)
        cf -= 1
    end
    ret
end
parseipv6fields(fields) = parseipv6fields(fields,8)

function parse(::Type{IPv6}, str::AbstractString)
    fields = split(str,':')
    if length(fields) > 8
        throw(ArgumentError("too many fields in IPv6 address"))
    elseif length(fields) == 8
        return IPv6(parseipv6fields(fields))
    elseif in('.',fields[end])
        return IPv6((parseipv6fields(fields[1:(end-1)],6))
            | parse(IPv4, fields[end]).host )
    else
        return IPv6(parseipv6fields(fields))
    end
end

#
# This support IPv4 addresses in the common dot (IPv4) or colon (IPv6)
# separated formats. Most other common formats use a standard integer encoding
# of the appropriate size and should use the appropriate constructor
#

function parse(::Type{IPAddr}, str::AbstractString)
    if ':' in str
        return parse(IPv6, str)
    else
        return parse(IPv4, str)
    end
end

macro ip_str(str)
    return parse(IPAddr, str)
end

struct InetAddr{T<:IPAddr}
    host::T
    port::UInt16
end

InetAddr(ip::IPAddr, port) = InetAddr{typeof(ip)}(ip, port)

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
## context

## A context stores objects where named values are looked up.
mutable struct Context
    view ## of what? a Dict, Module, CompositeKind, DataFrame.parent.
    parent ## a context or nothing
    _cache::Dict
end
function Context(view, parent=nothing)
    d = Dict()
    d["."] = isa(view, AnIndex) ? view.value : view
    Context(view, parent, d)
end

## used by renderTokens
function ctx_push(ctx::Context, view)
    Context(view, ctx) ## add context as a parent
end

function ctx_pop(ctx::Context)
    ctx.parent
end

# we have some rules here
# * Each part of a dotted name should resolve only against its parent.
# * Any falsey value prior to the last part of the name should yield ''.
# * The first part of a dotted name should resolve as any other name.
function lookup_dotted(ctx::Context, dotted)
    for key in split(dotted, ".")
        nctx = lookup(Context(ctx.view), key)
        falsy(nctx) && return nothing
        ctx = Context(nctx, ctx)
    end
    ctx.view
end

## Lookup value by key in the context
function lookup(ctx::Context, key)
    if haskey(ctx._cache, key)
        value = ctx._cache[key]
    else
        context = ctx
        value = nothing
        while value == nothing && context != nothing
            ## does name have a .?
            if occursin(r"\.", key)
                value = lookup_dotted(context, key)
                value != nothing && break

                ## do something with "."
                ## we use .[ind] to refer to value in parent of given index;
                m = match(r"^\.\[(.*)\]$", key)
                m == nothing && break
#                m == nothing && error("Not implemented. Can use Composite Kinds in the view.")

                idx = m.captures[1]
                vals = context.parent.view

                # this has limited support for indices: "end", or a number, but no
                # arithmetic, such as `end-1`.
                if isa(vals, Vector) # supports getindex(v, i)?
                   if idx == "end"
                       value = AnIndex(-1, vals[end])
                   else
                       ind = Base.parse(Int, idx)
                       value = AnIndex(ind, vals[ind])
                   end
                    break
                end
            else
                ## strip leading, trailing whitespace in key
                value = lookup_in_view(context.view, stripWhitespace(key))
            end
            context = context.parent

        end

        ## cache
        ctx._cache[key] = value
    end

    return(value)
end

## Lookup value in an object by key
## This of course varies based on the view.
## we special case dataframes here, so that we don't have to assume package is loaded
function lookup_in_view(view, key)
    if Tables.istable(view)
        if isempty(Tables.rows(view))
            return nothing
        else
            rows = Tables.rows(view)
            sch = Tables.schema(rows)
            if sch == nothing
                ## schema is unknown or non inferrable
                ## What to do?
                r = getfield(first(rows), 1)
                k = occursin(r"^:", key)  ? Symbol(key[2:end])  : key
                if isa(r, Pair)
                    return k == r.first  ? r.second : nothing
                else
                    return k in propertynames(r) ? getproperty(r, k) : nothing
                end
            else
                # work with a dictionary from the IteratorRow interface
                # follows "Sinks (transferring data from one table to another)"
                out = Any[]
                for row in rows
                    rD = Dict()
                    Tables.eachcolumn(sch, row) do val, col, name
                        rD[name] = val
                    end
                    k = occursin(r"^:", key)  ? Symbol(key[2:end])  : key
                    push!(out, get(rD,k, nothing))
                end
                return out
            end
        end
    elseif  is_dataframe(view)

        if occursin(r"^:", key)  key = key[2:end] end
        key = Symbol(key)
        out = nothing
        if haskey(view, key)
                out = view[1, key] ## first element only
        end
        out
    else
        _lookup_in_view(view, key)
    end
end

function _lookup_in_view(view::AbstractDict, key)
    ## is it a symbol?

    if occursin(r"^:", key)
        key = Symbol(key[2:end])
    end

    get(view, key, nothing)

end

function _lookup_in_view(view::NamedTuple, key)
    ## is it a symbol?
    if occursin(r"^:", key)
        key = Symbol(key[2:end])
    end
     if haskey(view, key)
        getindex(view, key)
    else
        nothing
    end
end

function _lookup_in_view(view::Module, key)

    hasmatch = false
    re = Regex("^$key\$")
    for i in names(view, all=true)
        if occursin(re, string(i))
            hasmatch = true
            break
        end
    end

    out = nothing
    if hasmatch
        out = getfield(view, Symbol(key))  ## view.key
    end
    out

end

## Default is likely not great,
function _lookup_in_view(view, key)

    if occursin(r"^:", key)
        k = Symbol(key[2:end])
    else
        k = key
    end

    # check propertyname, then fieldnames
    if k in propertynames(view)
        getproperty(view, k)
    else

        nms = fieldnames(typeof(view))

        # we match on symbol or string for fieldname
        if isa(k, Symbol)
            has_match = k in nms
        else
            has_match = Symbol(k) in nms
        end

        out = nothing
        if has_match
            out = getfield(view, Symbol(k))  ## view.key
        end

        out
    end
end


text:
# This script shows how one might generate a daylight plot
# for a custom time series time using the Quantity type found
# in my TidalFluxQuantities package
using TidalFluxQuantities, ADCPDataProcessing, PIEMetData, TidalFluxExampleData, DaylightPlots, TidalFluxPlots, Plots

pyplot()

# This stuff is all just loading and processing some
# example tidal data

setADCPdatadir!(Pkg.dir("TidalFluxExampleData","data","adcp"))
setmetdatadir!(Pkg.dir("TidalFluxExampleData","data","met"))

creek = Creek{:sweeney}()
deps = parse_deps(creek)
adata = load_data.(deps)

cs = parse_cs(creek)
csd = load_data(cs)
M = ADCPDataProcessing.depth_mask(adata[1])
H = atmoscorrect(adata[1])

# Here we define tlims. We use the argument `x...` to catch
# any additional arguments that might be passed to `daylight`,
# such as the depth masks that we use in ADCPDataProcessing
# to get rid of invalid data

DaylightPlots.tlims(q::Quantity,x...) = DaylightPlots.tlims(TidalFluxQuantities.times(q),quantity(q))

# Here we make our plot. `H::Stage` and `M::Mask` have
# recipes defined in TidalFluxPlots for plotting. Along
# with our definition of `tlims` above, this allows us
# to make daylight plots with our quantities.

daylight(H,M,loc=Point(42.35,-71.05),tz=tz"America/New_York",color=:black,leg=false,grid=false)

savefig("tides.png")


text:
function get(link::Link, ::Type{Any})
    # TODO: handle big integers/floats
    t = getnextraw(link)
    if t == TK_INT
        get(link, WInteger)
    elseif t == TK_UINT8
        get(link, UInt8)
    elseif t == TK_INT32_LE
        get(link, Int)
    elseif t == TK_INT64_LE
        get(link, Int64)
    elseif t == TK_FUNC || t == TK_ARRAY
        get(link, WExpr)
    elseif t == TK_STR
        get(link, String)
    elseif t == TK_REAL
        get(link, WReal)
    elseif t == TK_FLOAT64_LE
        get(link, Float64)
    elseif t == TK_SYM
        get(link, WSymbol)
    elseif t == TK_ERROR
        throw(MathLinkError(link))
    else
        # TODO: if we end up here then need to fix link
        error("Unsupported data type $t ($(Int(t)))")
    end
end

function get(link::Link, ::Type{WExpr})
    nargs = getargcount(link)
    head = get(link, Any)
    args = [get(link, Any) for i = 1:nargs]
    WExpr(head, args)
end
function put(link::Link, w::WExpr)
    puttype(link, TK_FUNC)
    putargcount(link, length(w.args))
    put(link, w.head)
    for arg in w.args
        put(link, arg)
    end
    nothing
end

function put(link::Link, list::Union{AbstractVector, Tuple})
    putfunction(link, WSymbol("List"), length(list))
    for x in list
        put(link, x)
    end
end

function get(link::Link, ::Type{T}) where {T<:Tuple}
    F = fieldtypes(T)
    nargs = getargcount(link)
    _ = get(link, Any)
    @assert nargs == length(F)
    map(S -> get(link, S), F)
end
function get(link::Link, ::Type{Vector{T}}) where {T}
    nargs = getargcount(link)
    _ = get(link, Any)
    T[get(link, T) for i = 1:nargs]
end
function get(link::Link, ::Type{Vector})
    nargs = getargcount(link)
    _ = get(link, Any)
    [get(link, Any) for i = 1:nargs]
end

get(link, ::Type{UInt8}) = get(link, Int8)

function put(link::Link, w::WInteger)
    puttype(link, TK_INT)
    put(link, w.value)
end
function put(link::Link, x::BigInt)
    put(link, WInteger(string(x)))
end
function get(link::Link, ::Type{BigInt})
    parse(BigInt, get(link, WInteger).value)
end

const log10_2 = log10(2)
function put(link::Link, x::BigFloat)
    str = Base.MPFR.string_mpfr(x,"%.Re")
    str = replace(str, "e" => "`$(log10_2*precision(x))*^")
    put(link, WReal(str))
end
function get(link::Link, ::Type{BigFloat})
    str = get(link, WReal).value
    BigFloat(replace(replace(str, r"`[0-9\.]*"=>""), "*^"=>"e"))
end

put(link::Link, ::typeof(pi)) = put(link, W"Pi")
put(link::Link, ::typeof(MathConstants.e)) = put(link, W"E")
put(link::Link, ::typeof(MathConstants.golden)) = put(link, W"GoldenRatio")
put(link::Link, ::typeof(MathConstants.eulergamma)) = put(link, W"EulerGamma ")
put(link::Link, ::typeof(MathConstants.catalan)) = put(link, W"Catalan ")


text:
function assign_scenario(s::BasePriceScenario, config::VanillaOptions, params::VanillaParameters)
    params.cpricebase
end

function assign_scenario(s::OptimalPriceScenario, config::VanillaOptions, params::VanillaParameters)
    fill(Inf, config.N)
end

function assign_scenario(s::Scenario, config::VanillaOptions, params::VanillaParameters)
    error("$(s) is not a valid scenario for v2013R (Vanilla flavour)");
end


text:
"""
# Exponential Search
It works in O(Log n) time
Exponential search involves two steps:
- Find range where element is present
- Do Binary Search in above found range.
### Time Complexity :
O(Log n)
### Auxiliary Space :
The above implementation of Binary Search is recursive and requires O(Log n) space. With iterative Binary Search, we need only O(1) space.
Applications of Exponential Search:
Exponential Binary Search is particularly useful for unbounded searches, where size of array is infinite. Please refer Unbounded Binary Search for an example.
It works better than Binary Search for bounded arrays, and also when the element to be searched is closer to the first element.
"""

"""
	Binary Search
"""
function binary_search(arr::AbstractArray{T,1}, l::T, r::T, x::T) where T <: Real
	n = size(arr)[1]
	if(r>=l)
		mid = Int(ceil(l+(r-l)/2));
		println(mid)
		if(arr[mid] == x)
			return "Element present at index $mid"
		elseif(arr[mid] > x)
			binary_search(arr, l, mid-1, x)
		else
			binary_search(arr, mid+1, r, x)
		end
	else
		return "Element not present in array"
	end
end

"""
	 exponential_search(arr::AbstractArray{T,1}, x::T) where {T <: Real}

Exponential Search in 1-D array
Time Complexity:  O(Log n)
"""
function exponential_search(arr::AbstractArray{T,1}, x::T) where {T <: Real}
	n = size(arr)[1]
	if(arr[1] == x)
		return "Elemenet present at index 1"
	end

	i = 1
	while( i<n && arr[i]<=x)
		i = i * 2
	end
	return binary_search( arr, Int(ceil(i / 2)), min(i, n), x)
end

# Arguments
arr = [1, 2, 3, 4, 13, 15, 20];
x = 4;
n = size(arr)[1]
l = 1;
r = n;

exponential_search(arr, x)

max_stars_repo_name: madisonestabrook/Julia

max_issues_repo_name: madisonestabrook/Julia

max_forks_repo_name: madisonestabrook/Julia


text:
using Revise
using Basins
using Plots
using DynamicalSystems
using DifferentialEquations

#
#
# ω=1.
# F = 0.2
# ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1)
# integ_df  = integrator(ds; alg=Tsit5(),  reltol=1e-8, save_everystep=false)
# xg = range(-2.2,2.2,length=150)
# yg = range(-2.2,2.2,length=150)
#
# @time bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω)
#
# sa,sb = compute_saddle(integ_df, bsn, [1], [2]; N=1000)

# plot(xg,yg,bsn.basin', seriestype=:heatmap)
# s = Dataset(sa)
# plot!(s[:,1],s[:,2],seriestype=:scatter)

# io = open("myfile.txt", "w");
# for v in s
#     write(io, string(v[1], " ", v[2], "; \n"));
# end
# close(io);

# Forced pendulum

# Equations of motion:
function forced_pendulum!(du, u, p, t)
    d = p[1]; F = p[2]; omega = p[3]
    du[1] = u[2]
    du[2] = -d*u[2] - sin(u[1])+ F*cos(omega*t)
end

# We have to define a callback to wrap the phase in [-π,π]
function affect!(integrator)
    if integrator.u[1] < 0
        integrator.u[1] += 2*π
    else
        integrator.u[1] -= 2*π
    end
end

condition(u,t,integrator) = (integrator.u[1] < -π  || integrator.u[1] > π)

cb = DiscreteCallback(condition,affect!)

#d, F ,w
F = 1.66
ω = 1.
d=0.2
p=[d, F, ω]
#p=[0.15, 0.2, 0.1]
df = ODEProblem(forced_pendulum!,rand(2),(0.0,20.0), p)
integ_df  = init(df, alg=AutoTsit5(Rosenbrock23()); reltol=1e-9, abstol=1e-9, save_everystep=false, callback=cb)

xres=200
yres=200

# range for forced pend
xg = range(-pi,pi,length=xres)
yg = range(-2.,4.,length=yres)

# compute basin
@time bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω)

Na = length(unique(bsn.basin))

sa,sb = compute_saddle(integ_df, bsn, [1], [2,3]; N=1000)

plot(xg,yg,bsn.basin', seriestype=:heatmap)
s = Dataset(sa)
plot!(s[:,1],s[:,2],seriestype=:scatter, markercolor=:blue)


text:
"""
This module demonstrates how to write an abstract type hierarchy and explore
it's shape and behaviours.
"""
module PersonalAssets

using InteractiveUtils

export Asset
export Property
export Investment
export Cash
export House
export Apartment
export FixedIncome
export Equity

abstract type Asset end

abstract type Property <: Asset end
abstract type Investment <: Asset end
abstract type Cash <: Asset end

abstract type House <: Property end
abstract type Apartment <: Property end

abstract type FixedIncome <: Investment end
abstract type Equity <: Investment end

export subtypetree

function subtypetree(roottype, level=1, indent=2)
  level == 1 && println(roottype)
  for type in subtypes(roottype)
    join(fill(" ", level*indent)) * string(type) |> println
    subtypetree(type, level+1, indent)
  end
end

end # PersonalAssets

max_stars_repo_name: BradLyman/learn-julia

max_issues_repo_name: BradLyman/learn-julia

max_forks_repo_name: BradLyman/learn-julia


text:
@symbol_func function beta_N(cur_reactor::AbstractReactor)
  cur_beta = 4.027e-2

  cur_beta *= cur_reactor.epsilon

  cur_beta *= cur_reactor.R_0

  cur_beta *= cur_reactor.n_bar

  cur_beta *= cur_reactor.T_bar

  cur_beta /= cur_reactor.B_0

  cur_beta /= cur_reactor.I_P

  cur_beta *= ( 1 + cur_reactor.f_D )

  cur_beta *= ( 1 + cur_reactor.nu_n )

  cur_beta *= ( 1 + cur_reactor.nu_T )

  cur_beta /= ( 1 + cur_reactor.nu_n + cur_reactor.nu_T )

  cur_beta
end


text:
function spherical_to_cartesian(ρ, θ, φ)
    x = ρ * sin(θ) * cos(φ)
    y = ρ * sin(θ) * sin(φ)
    z = ρ * cos(θ)
    Vector3(x, y, z)
end
spherical_to_cartesian(spherical::Vector3) = spherical_to_cartesian(spherical[1], spherical[2], spherical[3])

function cartesian_to_spherical(x, y, z)
    ρ = √(x ^ 2 + y ^ 2 + z ^ 2)
    θ = acos(z / ρ)
    φ = atan(y, x)
    Vector3(ρ, θ, φ)
end
cartesian_to_spherical(xyz::Vector3) = cartesian_to_spherical(xyz[1], xyz[2], xyz[3])

ρ_hat(θ, φ) = Vector3(sin(θ) * cos(φ), sin(θ) * sin(φ), cos(θ))
ρ_hat(Ω::Vector3) = normalize(Ω)
θ_hat(θ, φ) = Vector3(cos(θ) * cos(φ), cos(θ) * sin(φ), -sin(θ))
φ_hat(θ, φ) = Vector3(-sin(φ), cos(φ), 0)
∂Ω_∂ρ(θ, φ) = zero(Vector3)
∂Ω_∂θ(θ, φ) = Vector3(cos(θ) * cos(φ), cos(θ) * sin(φ), -sin(θ))
∂Ω_∂φ(θ, φ) = Vector3(-sin(θ) * sin(φ), sin(θ) * cos(φ), 0)

for func in (:θ_hat, :φ_hat, :∂Ω_∂ρ, :∂Ω_∂θ, :∂Ω_∂φ)
    @eval $func(Ω::Vector3) = let (_, θ, φ) = cartesian_to_spherical(Ω)
        $func(θ, φ)
    end
end


text:
function derivative{T<:FloatingPoint}(f::Function,x::Array{T,1})

    n = length(f(x))
    m = length(x)

	e  = eps(T)^(1/3)*Base.maxabs([x,one(T)])
    dh = eye(m)*e
	deriv = Array(T,n,m)

    for i = 1:m
	    f1 = f(x+2*dh[:,i])
	    f2 = f(x+dh[:,i])
	    f3 = f(x-dh[:,i])
	    f4 = f(x-2*dh[:,i])
        deriv[:,i] = (-f1+8*f2-8*f3+f4)/(12*e)
    end

    return deriv

end

max_stars_repo_name: dingliumath/SolveDSGE

max_issues_repo_name: dingliumath/SolveDSGE

max_forks_repo_name: dingliumath/SolveDSGE