text:
#
# see: https://github.com/JuliaWeb/Mux.jl
# and:
#

@debug "server starting up"

include( "uses.jl")
include( "types.jl")
include( "logger.jl")
include( "examples.jl")
include( "display_constants.jl")
include( "static_texts.jl")
include( "table_libs.jl")
include( "text_html_libs.jl")
include( "runner_libs.jl" )
include( "base_results.jl")
include( "queue_libs.jl")

# example for json3 StructTypes.@Struct T

function get_thing( thing::AbstractArray, key :: AbstractString, default :: AbstractString )
   for i in thing
      if i[1] == key
         return i[2]
      end
   end # loop
   default
end #get

function addqstrdict( app, req  :: Dict )
   req[:parsed_querystring] = qstrtodict(req[:query])
   @debug "req=$req"
   return app(req)
end

# Better error handling
function errorCatch( app, req  :: Dict )
   try
      app(req)
   catch e
      @error "Error occured!"
      io = IOBuffer()
      showerror(io, e)
      err_text = takebuf_string(io)
      @error err_text
      resp = withHeaders(JSON3.write(Dict("message" => err_text, "error" => true)), req)
      resp[:status] = 500
      return resp
   end
end

function d100( v :: Number ) :: Number
   v/100.0
end

function web_map_params( req  :: Dict ) :: TaxBenefitSystem
   querydict = req[:parsed_querystring]
   sys = deepcopy( BASE_PARAMS )
   d = req[:parsed_querystring]

   sys.ubi.abolished = false
   sys.ubi.adult_amount = d["bi_adult"]/WEEKS_PER_YEAR
   sys.ubi.child_amount = d["bi_child"]/WEEKS_PER_YEAR
   sys.ubi.universal_pension  = d["bi_pensioner"]/WEEKS_PER_YEAR
   sys.ubi.adult_age = d["bi_adult_age"]
   sys.ubi.retirement_age = d["bi_pens_age"]
   sys.ubi.mt_bens_treatment =
      if d["ubi_mtbens_abolish"]
         ub_abolish
      elseif d["ubi_mtbens_keep_as_is"]
         ub_as_is
      elseif d["ubi_mtbens_keep_housing"]
         ub_keep_housing
      else
         @assert false "no assignment for ubi.mt_bens_treatment"
      end
   sys.ubi.abolish_sickness_bens = d["ubi_abolish_sick"]
   sys.ubi.abolish_pensions = d["ubi_abolish_pensions"]
   sys.ubi.abolish_jsa_esa = d["ubi_abolish_esa"]
   sys.ubi.abolish_others = d["ubi_abolish_others"]
   sys.ubi.ub_as_mt_income = d["ubi_as_mt_income"]
   sys.ubi.ub_taxable = d["ubi_taxable"]
   sys.it.personal_allowance = d["it_pa"]/WEEKS_PER_YEAR
   br = d["it_basic_rate"] /=100.0
   if br == 0
      sys.it.non_savings_rates[1:3] .= 0.0
   else
      bincr = br-sys.it.non_savings_rates[2]
      sys.it.non_savings_rates[1:3] .+= bincr
      sys.it.non_savings_rates[1] = max(0, sys.it.non_savings_rates[1])
   end
   sys.it.non_savings_rates[4] = d["it_higher_rate"] / 100.0
   sys.it.non_savings_rates[5] = d["it_top_rate"] / 100.0
   make_ubi_pre_adjustments!( sys )
   return sys
end

function web_map_settings( req  :: Dict ) :: Settings
   querydict = req[:parsed_querystring]
   settings = deepcopy( BASE_SETTINGS )
   # ...
   return settings
end

# Headers -- set Access-Control-Allow-Origin for either dev or prod
# this is from https://github.com/JuliaDiffEq/DiffEqOnlineServer
#
function add_headers( json :: AbstractString ) :: Dict
    headers  = HttpCommon.headers()
    headers["Content-Type"] = "application/json; charset=utf-8"
    headers["Access-Control-Allow-Origin"] = "*"
    Dict(
       :headers => headers,
       :body=> json
    )
end

#
# This is my attempt at starting a task using the 1.3 @spawn macro
# 1st parameter is a function that returns String (probably a json string) and accepts the req Dict
#
function do_in_thread( the_func::Function, req :: Dict ) :: Dict
   response = @spawn the_func( req )
   # note that the func returns a string but response is a Future type
   # line below converts response to a string
   @debug "do_in_thread response is $response"
   json = fetch( response )
   add_headers( json )
end

function get_results( u :: AbstractString ):: Dict
   uuid = UUID(u)
   if haskey( CACHED_RESULTS[res.cache_key])

   end
end

function get_progress( u :: AbstractString ) :: Dict
   uuid = UUID(u)
   state = ( uuid=u, phase="missing", count=0, total=0 )
   if haskey( PROGRESS, uuid )
      @debug "get_progress entered uuid=|$u| key is in progress;"
      p = PROGRESS[uuid]
      @debug "phase is $(p.progress.phase)"
      if p.progress.phase == "end"
         @debug "get_progress: phase end"
         state = STASHED_RESULTS[uuid]
         @debug "got state for STASHED_RESULTS; uuid=$uuid"
         delete!( PROGRESS, uuid )
         @debug( "deleted $uuid progress")
      else
         state = ( uuid=p.progress.uuid, phase=p.progress.phase, count=p.progress.count, total=11_000 )
      end
   end
   @debug "PROGRESS now $PROGRESS"
   json = JSON3.write( state )
   return add_headers( json )
end

function get_base( req :: Dict )
   @debug "returning base results"
   return add_headers( JSON3.write(BASE_TEXT_OUTPUT))
 end

function submit_model( req :: Dict )
   query = req[:query]
   @debug "submit model entered with query $query"
   if haskey( CACHED_RESULTS, query )
      @debug "returning results from cache"
      return add_headers( JSON3.write(CACHED_RESULTS[query]))
   end
   sys = web_map_params( req )
   @debug "submit_model  - starting mapping params"
   settings = web_map_settings( req )
   uuid = submit_job( req[:query], sys, settings )
   @debug "submit_model uuid=$uuid"
   json = add_headers( JSON3.write( uuid ))
   return json
end

#
# from diffeq thingy instead of Mux.defaults
#
# ourstack = stack(Mux.todict, errorCatch, Mux.splitquery, Mux.toresponse) # from DiffEq
#
@app stbdemo = (
   Mux.defaults,
   addqstrdict,
   page( respond("
A Scottish Tax-Benefit Model
")),
   # handle main stb run a bit differently so we can
   # check & cache the results.
   # page("/stb", req -> web_do_one_run_cached( req )),
   page( "/bi/progress/:uuid", req -> get_progress((req[:params][:uuid]))), # note no Headers
   page("/bi/run/", req -> submit_model( req )),
   page("/bi/base/", req -> get_base( req )),
   Mux.notfound(),
)


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")

polyorder = 1
numqp = 2

elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
levelsetbasis = HermiteTensorProductBasis(2)

x0 = [0.0, 0.0]
meshwidths = [3.0, 1.0]
nelements = [3, 1]
points = interpolation_points(elasticitybasis)
dim, nf = size(interpolation_points(levelsetbasis))
dgmesh = CutCellDG.DGMesh(x0, meshwidths, nelements, points)
cgmesh = CutCellDG.CGMesh(x0, meshwidths, nelements, nf)

normal = [1.0, 0.0]
xI = [0.5, 0.0]
levelset = CutCellDG.LevelSet(
    x -> plane_distance_function(x, normal, xI),
    cgmesh,
    levelsetbasis,
)
cutmesh = CutCellDG.CutMesh(dgmesh, levelset)

interfacequads = CutCellDG.InterfaceQuadratures(cutmesh, levelset, numqp)
@test size(interfacequads.quads) == (2, 1)
testcelltoquad = [1, 0, 0]
@test allequal(interfacequads.celltoquad, testcelltoquad)

testnormals = repeat(normal, inner = (1, numqp))
@test allapprox(
    testnormals,
    CutCellDG.interface_normals(interfacequads, 1),
    1e3eps(),
)
@test interfacequads[1, 1] ≈ interfacequads[-1, 1]

testtangents = repeat([0.0, 1.0], inner = (1, numqp))
testscaleareas = repeat([0.5], numqp)
@test allapprox(
    testtangents,
    CutCellDG.interface_tangents(interfacequads, 1),
    1e3eps(),
)
@test allapprox(
    testscaleareas,
    CutCellDG.interface_scale_areas(interfacequads, 1),
)

ext: jl


text:
@symbol_func function K_WM(cur_reactor::AbstractReactor)
  cur_K = volume(cur_reactor)

  cur_K /= 2

  cur_K /= 4 * cur_reactor.pi * 1e-7

  cur_K /= safe_symbol(cur_reactor, :R_0) ^ 3

  cur_K *= 1e-9

  cur_K
end


text:
"""
RDA (R data archive) reading context.

* Stores flags that define how R objects are read and converted into Julia objects.
* Maintains the list of R objects that could be referenced later in the RDA stream.
"""
struct RDAContext{T<:RDAIO}
    io::T                      # RDA input stream

    # RDA format properties
    fmtver::UInt32             # format version
    Rver::VersionNumber        # R version used to write the file
    Rmin::VersionNumber        # minimal R version to read the file
    natenc::Union{String, Nothing} # native encoding at serialization time

    kwdict::Dict{Symbol,Any}   # options defining RDA deserialization behaviour

    # intermediate data
    ref_tab::Vector{RSEXPREC}  # SEXP array for references
end

int2ver(v::Integer) = VersionNumber(v >> 16, (v >> 8) & 0xff, v & 0xff)

function RDAContext(io::RDAIO; kwoptions...)
    fmtver = readuint32(io)
    rver = int2ver(readint32(io))
    rminver = int2ver(readint32(io))
    natenc = fmtver == 3 ? readnchars(io, readint32(io)) : nothing
    kwdict = Dict{Symbol,Any}(kwoptions)
    RDAContext(io, fmtver, rver, rminver, natenc, kwdict, RSEXPREC[])
end

"""
Register R object, so that it could be referenced later
(by its index in the reference table).
"""
function registerref!(ctx::RDAContext, obj::RSEXPREC)
    push!(ctx.ref_tab, obj)
    obj
end


text:
# https://www.edfplus.info/specs/edftexts.html

# This is a mapping from Onda-appropriate unit names to EDF standard physical dimensions.
# This doesn't cover everything, just units that we've seen often. Note that casing
# matters for the EDF standard physical dimensions, so we can't just always normalize
# to a single case when processing; this results in some redundant-looking values below.
const STANDARD_UNITS = Dict("nanovolt" => ["nV"],
                            "microvolt" => ["uV", "\xb5V", "\xb5V\xb2"],
                            "millivolt" => ["mV"],
                            "centivolt" => ["cV"],
                            "decivolt" => ["dV"],
                            "volt" => ["V", "v", "volts"],
                            "millimiter" => ["mm"],
                            "centimeter" => ["cm"],
                            "milliliter" => ["mL", "ml"],
                            "degrees_celsius" => ["degC", "degc"],
                            "degrees_fahrenheit" => ["degF", "degf"],
                            "kelvin" => ["K"],
                            "percent" => ["%"],
                            "liter_per_minute" => ["L/m", "l/m", "LPM", "Lpm", "lpm", "LpM", "L/min", "l/min"],
                            "millimeter_of_mercury" => ["mmHg", "mmhg", "MMHG"],
                            "beat_per_minute" => ["B/m", "b/m", "bpm", "BPM", "BpM", "Bpm"],
                            "centimeter_of_water" => ["cmH2O", "cmh2o", "cmH20"],
                            "ohm" => ["Ohm", "ohms", "Ohms", "ohm"],
                            "unknown" => ["", "\"\"", "#", "u", "none", "---"],
                            "relative" => ["rel."])

# The case-sensitivity of EDF physical dimension names means you can't/shouldn't
# naively convert/lowercase them to compliant Onda unit names, so we have to be
# very conservative here and error if we don't recognize the input.
function edf_to_onda_unit(edf_physical_dimension::AbstractString, unit_alternatives=STANDARD_UNITS)
    edf_physical_dimension = replace(edf_physical_dimension, r"\s"=>"")
    for (onda_unit, potential_edf_matches) in unit_alternatives
        any(==(edf_physical_dimension), potential_edf_matches) && return onda_unit
    end
    error("""
          Failed to convert EDF physical dimension label `$(edf_physical_dimension)`
          to known Onda unit; please either open a PR to add this unknown unit
          to `OndaEDF.STANDARD_UNITS` (if the unit obeys the EDF standard),
          or otherwise preprocess your EDF such that physical dimension labels
          contain known/standard values.
          """)
end

function onda_to_edf_unit(onda_sample_unit::String)
    haskey(STANDARD_UNITS, onda_sample_unit) && return first(STANDARD_UNITS[onda_sample_unit])
    error("""
          Failed to convert Onda unit `$(onda_sample_unit)` to EDF physical dimension
          label; please either open a PR to add this unknown unit to `OndaEDF.STANDARD_UNITS`
          (if the unit obeys the EDF standard), or otherwise preprocess your input data such
          that their unit names are known/standard values.
          """)
end

const STANDARD_LABELS = Dict(# This EEG channel name list is a combined 10/20 and 10/10
                             # physical montage; channels are ordered from left-to-right,
                             # front-to-back w.r.t a top-down, nose-up view of the head.
                             # 10/20 channel names that aren't in the 10/10 system (and
                             # vice versa) are interleaved into their appropriate location
                             # on the head relative to other channels.
                             #
                             # Note that other commonly used (but non-EDF-standard) labels
                             # used as references include `LE` ("Linked Ear") and `AR`
                             # ("Average Reference"). See the following paper for details:
                             # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5479869/
                             ["eeg"] => ["pg1", "nz", "pg2",
                                         "fp1", "fpz", "fp2",
                                         "af7", "af3", "afz", "af4", "af8",
                                         "f9", "f7", "f5", "f3", "f1", "fz", "f2", "f4", "f6", "f8", "f10",
                                         "ft9", "ft7", "fc5", "fc3", "fc1", "fcz", "fc2", "fc4", "fc6", "ft8", "ft10",
                                         "a1", "m1", "t9", "t7", "t3", "c5", "c3", "c1", "cz", "c2", "c4", "c6", "t4", "t8", "t10", "a2", "m2",
                                         "tp9", "tp7", "cp5", "cp3", "cp1", "cpz", "cp2", "cp4", "cp6", "tp8", "tp10",
                                         "t5", "p9", "p7", "p5", "p3", "p1", "pz", "p2", "p4", "p6", "p8", "p10", "t6",
                                         "po7", "po3", "poz", "po4", "po8",
                                         "o1" => ["01"], "oz", "o2" => ["02"],
                                         "iz"],
                             # It is very common in the wild to see "EKG1", "EKG2", etc., and it's not possible
                             # by label alone to tell whether such channels refer to I, II, etc. or aVL, aVR,
                             # etc., so there's a burden on users to preprocess their EKG labels
                             ["ecg", "ekg"] => ["i" => ["1"], "ii" => ["2"], "iii" => ["3"],
                                                "avl"=> ["ecgl", "ekgl", "ecg", "ekg", "l"], "avr"=> ["ekgr", "ecgr", "r"], "avf",
                                                "v1", "v2", "v3", "v4", "v5", "v6", "v7", "v8", "v9",
                                                "v1r", "v2r", "v3r", "v4r", "v5r", "v6r", "v7r", "v8r", "v9r",
                                                "x", "y", "z"],
                             # EOG should not have any channel names overlapping with EEG channel names
                             ["eog", "eeg"] => ["left"=> ["eogl", "loc", "lefteye", "leye", "e1", "eog1", "l", "left eye", "leog", "log", "li"],
                                                "right"=> ["eogr", "roc", "righteye", "reye", "e2", "eog2", "r", "right eye", "reog", "rog", "re"]],
                             ["emg"] => ["chin1" => ["chn", "chin_1", "chn1", "kinn", "menton", "submental", "submentalis", "submental1", "subm1", "chin", "mentalis", "chinl", "chinli", "chinleft", "subm_1", "subment"],
                                         "chin2" => ["chn2", "chin_2", "submental2", "subm2", "chinr", "chinre", "chinright", "subm_2"],
                                         "chin3" => ["chn3", "submental3", "subm3", "chincenter"],
                                         "intercostal"=> ["ic"],
                                         "left_anterior_tibialis"=> ["lat", "lat1", "l", "left", "leftlimb", "tibl", "tibli", "plml", "leg1", "lleg", "lleg1", "legl", "jambe_l", "leftleg"],
                                         "right_anterior_tibialis"=> ["rat", "rat1", "r", "right", "rightlimb", "tibr", "tibre", "plmr", "leg2", "leg3", "rleg", "rleg1", "legr", "jambe_r", "rightleg"]],
                             # it is common to see ambiguous channels, which could be leg or face EMG
                             # if leg EMG is present in separate channels,
                             # post-processing might map "emg_ambiguous"
                             # to chin channels (for example)
                             ["emg_ambiguous", "emg"] => ["1" => ["aux1", "l"],
                                                          "2" => ["aux2", "r"],
                                                          "3" => ["aux3"]],
                             ["heart_rate"] => ["heart_rate"=> ["hr", "pulse", "pulso", "pr", "pulserate"]],
                             ["snore"] => ["snore" => ["ronquido", "ronquido derivad", "schnarchen", "ronfl", "schnarchmikro"]],
                             ["positive_airway_pressure", "pap"] => ["ipap", "epap", "cpap"],
                             ["pap_device_cflow"] => ["pap_device_cflow"=> ["cflow", "airflow", "flow"]],
                             ["pap_device_cpres"] => ["pap_device_cpres"=> ["cpres"]],
                             ["pap_device_leak"] => ["pap_device_leak"=> ["leak", "airleak"]],
                             ["ptaf"] => ["ptaf"],
                             ["respiratory_effort"] => ["chest" => ["thorax", "torax", "brust", "thor"], "abdomen"=> ["abd", "abdo", "bauch"]],
                             ["tidal_volume"] => ["tidal_volume"=> ["tvol", "tidal"]],
                             ["spo2"] => ["spo2"],
                             ["sao2"] => ["sao2", "osat"],
                             ["etco2"] => ["etco2 "=> ["capno"]])

const STANDARD_EXTRACTORS = [edf -> extract_channels_by_label(edf, signal_names, channel_names)
                             for (signal_names, channel_names) in STANDARD_LABELS]


text:
#=
A Region is a (possibly infinite) set of locations in space.
Regions include points, curves, surfaces, and volumes
A ParametricRegion and an ImplicitRegion are examples of a Region.
=#

abstract type AbstractRegion end
show(io::IO, s::AbstractRegion) =
  print(io, "$(typeof(s))(...)")

AbstractRegions = Vector{<:AbstractRegion}

abstract type Region0D <: AbstractRegion end
abstract type Region1D <: AbstractRegion end
abstract type Region2D <: AbstractRegion end
abstract type Region3D <: AbstractRegion end

macro defregion(name_typename, parent, fields...)
  (name, typename) = name_typename isa Symbol ?
    (name_typename, Symbol(string(map(uppercasefirst,split(string(name_typename),'_'))...))) :
    name_typename.args
  name_str = string(name)
  struct_name = esc(typename)
  field_names = map(field -> field.args[1].args[1], fields)
  field_types = map(field -> esc(field.args[1].args[2]), fields)
  field_inits = map(field -> field.args[2], fields)
  field_renames = map(Symbol ∘ string, field_names)
  field_replacements = Dict(zip(field_names, field_renames))
  struct_fields = map((name,typ) -> :($(name) :: $(typ)), field_names, field_types)
  mk_param(name,typ,init) = Expr(:kw, name, init)
  opt_params = map(mk_param, field_renames, field_types, map(init -> replace_in(init, field_replacements), field_inits))
  key_params = map(mk_param, field_names, field_types, field_renames)
  constructor_name = esc(name)
  predicate_name = esc(Symbol("is_", name_str))
  mk_convert(name,typ) = :($(esc(name)))
  field_converts = map(mk_convert, field_names, field_types)
  selector_names = map(field_name -> esc(Symbol(name_str, "_", string(field_name))), field_names)
  quote
    export $(constructor_name), $(struct_name), $(predicate_name) #, $(selector_names...)
    struct $struct_name <: $parent
      $(struct_fields...)
    end
    @noinline $(constructor_name)($(opt_params...); $(key_params...)) =
      $(struct_name)($(field_converts...))
    $(predicate_name)(v::$(struct_name)) = true
    $(predicate_name)(v::Any) = false
    #$(map((selector_name, field_name) -> :($(selector_name)(v::$(struct_name)) = v.$(field_name)),
    #      selector_names, field_names)...)
    KhepriBase.meta_program(v::$(struct_name)) =
        Expr(:call, $(Expr(:quote, name)), $(map(field_name -> :(meta_program(v.$(field_name))), field_names)...))
  end
end

@defregion(annulus, Region2D, center::Loc=u0(), outer_radius::Real=1, inner_radius::Real=1/2)

area(r::Annulus) = annulus_area(r.outer_radius, r.inner_radius)


text:
# # Multivariate maps

using SimpleSDMLayers
using StatsPlots
using StatsPlots.PlotMeasures
using Statistics

# **Justification for this use case:** We can show more than one (specifically,
# two or three) variables on a single map, using a bivariate or trivariate color
# scale. In order to illustrate these mappings, we will look at the joint
# distribution of three measures: eveness of land use, terrain roughness, and
# proportion of urbanized land.

boundaries = (left=-12.0, right=30.0, bottom=36.0, top=72.0)
layer1 = convert(
    Float16,
    SimpleSDMPredictor(EarthEnv, HabitatHeterogeneity, 2; resolution=5, boundaries...),
)
layer2 = convert(
    Float16, SimpleSDMPredictor(EarthEnv, Topography, 7; resolution=5, boundaries...)
)

# The landcover data are finer than the other layers, so we will coarsen them.
# Because of rounding issues, the left and right cordinates need to be rounded.

layer3 = coarsen(
    convert(Float16, SimpleSDMPredictor(EarthEnv, LandCover, 9; boundaries...)),
    mean,
    (5, 5),
)
layer3.left = round(layer3.left; digits=0)
layer3.right = round(layer3.right; digits=0)

# We finally mask everything according to the first layer

layer2 = mask(layer1, layer2);
layer3 = mask(layer1, layer3);

# Note that bivariate maps usually work best when used with 9 classes in total
# (so 3 for each side). The next decision is to take a bivaraite color palette,
# and the combinations below are [commonly
# used](https://www.joshuastevens.net/cartography/make-a-bivariate-choropleth-map/).
# Note that you can definitely use [diverging
# colors](https://www.personal.psu.edu/cab38/ColorSch/Schemes.html) if you want.
# If you use colors in the RGBA format (*e.g.* `colorant"#ef0ce8c4"`), the color
# map will account for transparency.

p0 = colorant"#e8e8e8"
bv_pal_1 = (p0=p0, p1=colorant"#64acbe", p2=colorant"#c85a5a")
bv_pal_2 = (p0=p0, p1=colorant"#73ae80", p2=colorant"#6c83b5")
bv_pal_3 = (p0=p0, p1=colorant"#9972af", p2=colorant"#c8b35a")
bv_pal_4 = (p0=p0, p1=colorant"#be64ac", p2=colorant"#5ac8c8")

# The bivariate map itself is a call to plot. Internally, this will transform
# the layers into quantiles (determined by the `classes` keyword, defaults to
# 3):

plot(layer1, layer3; st=:bivariate, bv_pal_3...)

# Note that you can use the `bivariate` shorthand as well:

pl1 = bivariate(layer1, layer3; classes=3, frame=:box, bv_pal_4...)
xaxis!(pl1, "Longitude")
yaxis!(pl1, "Latitude")

# We can repeat essentially the same process for the legend:

pl2 = bivariatelegend(layer1, layer3; classes=3, bv_pal_4...)
xaxis!(pl2, layernames(EarthEnv, HabitatHeterogeneity, 2))
yaxis!(pl2, layernames(EarthEnv, LandCover, 9))

# And now, we can plot the legend next to the map - future releases of the
# package will hopefully offer this in a far more user friendly way.

plot(pl1, pl2; layout=@layout [a{0.75w} b])

# Using the `subplot` and `inset` arguments of Plots.jl, we can have the legend
# within the figure. Note how in this example we expand the limits on the x axis
# to make the legend fit, but also use more classes in the map to have a
# smoother result.

p1 = bivariate(
    layer1,
    layer3;
    classes=6,
    bv_pal_2...,
    frame=:box,
    xlim=(-24, maximum(longitudes(layer1))),
)
xaxis!(p1, "Longitude")
yaxis!(p1, "Latitude")
p2 = bivariatelegend!(
    layer1,
    layer3;
    bv_pal_2...,
    inset=(1, bbox(0.04, 0.05, 0.28, 0.28, :top, :left)),
    subplot=2,
    xlab=layernames(EarthEnv, HabitatHeterogeneity, 2),
    ylab=layernames(EarthEnv, LandCover, 9),
    guidefontsize=7,
)

# Using a trivariate mapping follows the same process, with layers representing
# the red, green, and blue channel respectively.

plot(layer1, layer2, layer3; st=:trivariate, frame=:grid)

# There are two options for this type of plots. The first is `quantiles=true`
# (which maps quantiles rather than raw values), and the second is
# `simplex=false`, which makes all values sum to 1 within a pixel. For example:

trivariate(layer1, layer2, layer3; quantiles=true, simplex=true, frame=:grid)

# It is a good idea to question whether using `simplex` is appropriate. The
# legend can also be plotted using `trivariatelegend`:

trivariatelegend(layer1, layer2, layer3; quantiles=true, simplex=true)

# The legend function admits three additional arguments for the names of the
# `red`, `green`, and `blue` channels:

trivariatelegend(
    layer1,
    layer2,
    layer3;
    quantiles=true,
    simplex=true,
    red="Heterogeneous",
    green="Rough",
    blue="Urbanized",
)

# We can also combine the two elements:

tri1 = trivariate(
    layer1,
    layer2,
    layer3;
    xlim=(-24, maximum(longitudes(layer1))),
    frame=:grid,
    grid=false,
    simplex=false,
)
xaxis!(tri1, "Longitude")
yaxis!(tri1, "Latitude")
tri2 = trivariatelegend!(
    layer1,
    layer2,
    layer3;
    inset=(1, bbox(0.01, 0.03, 0.35, 0.35, :top, :left)),
    subplot=2,
    red="Heterogeneity",
    green="Roughness",
    blue="Urban",
    annotationfontsize=6,
    simplex=false,
)


text:
"""
Methods to find a low-rank subset
"""

"""
Return (error, indices) for the subset of size n, whose k-rank approximation has lowest error, as measured by the 2-norm.
Checks all subsets of size n. Has complexity Ω(|data_set|^n).
"""
function find_low_rank_subset_checkall(data_set,n,k)
    f(S) = let M = hcat([data_set[i] for i in S]...);
        truncated_svd(M, k)[1], S end
    minimum(f , powerset(1:length(data_set), n, n) )
end

find_low_rank_subset_checkall(A :: Array{T, 2}, n, k) where T <: Number = find_low_rank_subset_checkall(collect(eachcol(A)), n, k)

"""
Return (error, indices) for the subset of size n, whose k-rank approximation with vectors from data set has lowest error, as measured by the 2-norm.
Checks all subsets of size k. Has complexity Ω(|data_set|^k).
"""
function find_low_rank_subset_sample_rep(data_set,n,k)
    f(E) = let M = hcat([data_set[i] for i in E]...);
    M_times_pinv = M * pinv(M);
    distances_squared = [let d = v - M_times_pinv*v; dot(d,d) end for v in data_set]
    indices, distances = findmin(distances_squared, 1:n);
    sum(distances), indices
    end
    minimum(f, powerset(1:length(data_set), k, k))
end

find_low_rank_subset_sample_rep(A :: Array{T, 2}, n, k) where T <: Number = find_low_rank_subset_sample_rep(collect(eachcol(A)), n, k)

"""
Return an array of current_trajectories for each iteration step.
"""
function FastLowRank(data_set, n, k; initial_indices = [], step = L2_step, num_trajectories = k*trunc(Int(length(data_set)/n)), num_iterations = 50, verbose = false)
    current_trajectories = Vector{Trajectory}(undef, num_trajectories)
    # all_trajectories = Vector{Vector{Trajectory}}(undef, num_iterations)

    # initialize
    for t in 1:length(initial_indices)
        current_trajectories[t] = LowRank.Trajectory(data_set, n, k, initial_indices[t])
    end
    for t in length(initial_indices)+1:num_trajectories
        current_trajectories[t] = LowRank.Trajectory(data_set, n, k)
    end

    # do iterations
    is_converged = [false for i in 1:num_trajectories]
    best_converged_error = Inf
    best_error = Inf
    for i in 1:num_iterations
        for (ind, traj) in enumerate(current_trajectories)
            is_converged[ind] && continue
            new_traj = step(traj)
            best_error = min(best_error, new_traj.error)
            if Set(traj.indices) == Set(new_traj.indices)
                is_converged[ind] = true
                best_converged_error = min(best_converged_error, new_traj.error)
            end
            current_trajectories[ind] = new_traj
        end

        # @show best_error, best_converged_error
        if best_converged_error - best_error < eps(best_error)
            verbose && println("converged in $i steps")
            break
        elseif i == num_iterations
            verbose && println("did not converge in $i steps")
        end
    end
    best_error, best_trajectory_index = findmin((x -> x.error).(current_trajectories))
    best_trajectory = current_trajectories[best_trajectory_index]
    verbose && println("error = $best_error")
    return best_trajectory.error, best_trajectory.indices
end

FastLowRank(A :: Array{T, 2}, n, k; kwargs...) where T <: Number = FastLowRank(collect(eachcol(A)), n, k; kwargs...)


text:
using Econometrics, Test
using Econometrics: ConvergenceException, Hermitian
using CSV, RDatasets
# Exceptions
@testset "Exceptions" begin
	data = dataset("Ecdat", "Crime")
	data.AvgSen2 = 2 * data.AvgSen
	# Rank-deficient model
	@test isa(fit(EconometricModel,
			      @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + AvgSen2),
				  data),
			  EconometricModel)
	@test isa(fit(RandomEffectsEstimator,
				  @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris),
				  data,
				  panel = :County,
				  time = :Year),
			  EconometricModel)
	# Model specifics
	@test_throws(ArgumentError("The between estimator requires the panel identifier."),
				 fit(BetweenEstimator,
				 	 @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(Year)),
					 data))
	@test_throws(ArgumentError("Absorbing features is only implemented for least squares."),
				 fit(BetweenEstimator,
				 	 @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(Year)),
					 data,
					 panel = :County))
	@test_throws(ArgumentError("The random effects estimator requires the panel and temporal identifiers."),
				 fit(RandomEffectsEstimator,
				 	 @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(Year)),
					 data,
					 panel = :County))
	@test_throws(ArgumentError("Absorbing features is only implemented for least squares."),
				 fit(RandomEffectsEstimator,
				 	 @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(Year)),
					 data,
					 panel = :County,
					 time = :Year))
	data = dataset("datasets", "iris")
	@test_throws(ConvergenceException,
				 fit(EconometricModel,
				 	 @formula(Species ~ SepalLength + SepalWidth),
				 	 data))
	# Robust variance covariance estimators
	@test_throws(ArgumentError("Robust variance covariance estimators only available for continous response models."),
				 fit(EconometricModel,
				 	 @formula(Species ~ SepalLength + SepalWidth),
				 	 data,
					 vce = HC1))
	@test_throws(ArgumentError("When absorbing features, only HC2-HC4 are unavailable."),
				 fit(EconometricModel,
				 	 @formula(SepalLength ~ SepalLength + absorb(Species)),
				 	 data,
					 vce = HC2))
end
# Linear Models
@testset "Linear Models" begin
	@testset "Balanced Panel Data" begin
    	# Balanced Panel Data
    	data = dataset("Ecdat", "Crime")
    	# Pooling
    	model = fit(EconometricModel,
                	@formula(CRMRTE ~ PrbConv + AvgSen + PrbPris),
                	data)
		@test sprint(show, model) == "Continuous Response Model\nNumber of observations: 630\nNull Loglikelihood: 1633.30\nLoglikelihood: 1643.30\nR-squared: 0.0312\nLR Test: 19.99 ∼ χ²(3) ⟹  Pr > χ² = 0.0002\nFormula: CRMRTE ~ 1 + PrbConv + AvgSen + PrbPris\nVariance Covariance Estimator: OIM\n──────────────────────────────────────────────────────────────────────────────────────\n                    PE           SE      t-value  Pr > |t|         2.50%        97.50%\n──────────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.0186413    0.00431066    4.32447    <1e-4    0.0101762     0.0271064  \nPrbConv      -0.00116473   0.000422062  -2.75961    0.0060  -0.00199356   -0.0003359  \nAvgSen        0.000236185  0.000268216   0.88058    0.3789  -0.000290526   0.000762897\nPrbPris       0.0273394    0.00817642    3.34369    0.0009   0.0112829     0.0433959  \n──────────────────────────────────────────────────────────────────────────────────────"
    	@test nullloglikelihood(model) ≈ 1633.305 rtol = 1e-3
    	@test loglikelihood(model) ≈ 1643.304 rtol = 1e-3
    	@test dof(model) == 5
    	# This corresponds to BIC from R's lm, Stata uses rank(vcov(model))
    	@test bic(model) ≈ -3254.379 rtol = 1e-3
    	@test coef(model) ≈ [0.0186413, -0.0011647, 0.0002362, 0.0273394] rtol = 1e-5
    	@test vcov(model) ≈ [ 1.858177e-05 -1.620480e-07 -6.464872e-07 -2.860941e-05
                          	 -1.620480e-07  1.781365e-07 -1.715589e-09  1.286486e-07
                          	 -6.464872e-07 -1.715589e-09  7.193971e-08  8.182148e-09
                             -2.860941e-05  1.286486e-07  8.182148e-09  6.685387e-05
                        	] rtol = 1e-6
    	@test confint(model) ≈ [ 0.0101761950  0.0271063941
                            	-0.0019935581 -0.0003358998
                            	-0.0002905263  0.0007628970
                             	0.0112828681  0.0433959412
                           	   ] rtol = 1e-6
    	# Between
    	model = fit(BetweenEstimator,
                	@formula(CRMRTE ~ PrbConv + AvgSen + PrbPris),
                	data,
					panel = :County)
		@test sprint(show, model) == "Between Estimator\nCounty with 630 groups\nBalanced groups with size 7\nNumber of observations: 90\nNull Loglikelihood: 239.56\nLoglikelihood: 244.40\nR-squared: 0.1029\nWald: 3.29 ∼ F(3, 86) ⟹ Pr > F = 0.0246\nFormula: CRMRTE ~ 1 + PrbConv + AvgSen + PrbPris\nVariance Covariance Estimator: OIM\n───────────────────────────────────────────────────────────────────────────────────\n                    PE          SE       t-value  Pr > |t|        2.50%      97.50%\n───────────────────────────────────────────────────────────────────────────────────\n(Intercept)  -0.00464339   0.0186571   -0.24888     0.8040  -0.0417325   0.0324457 \nPrbConv      -0.00354113   0.0018904   -1.87322     0.0644  -0.00729911  0.00021685\nAvgSen        0.000848049  0.00116477   0.728086    0.4685  -0.00146743  0.00316353\nPrbPris       0.0730299    0.0332057    2.19932     0.0305   0.0070191   0.139041  \n───────────────────────────────────────────────────────────────────────────────────"
    	@test nullloglikelihood(model) ≈ 239.5638 rtol = 1e-3
    	@test loglikelihood(model) ≈ 244.4479 rtol = 1e-3
    	@test dof(model) == 5
    	@test coef(model) ≈ [-0.00464339, -0.00354113, 0.00084805, 0.07302994] rtol = 1e-6
    	@test vcov(model) ≈ [ 3.480883e-04 -5.985740e-06 -1.382250e-05 -5.104706e-04
                         	 -5.985740e-06  3.573593e-06 -1.212431e-08  8.538922e-06
                         	 -1.382250e-05 -1.212431e-08  1.356678e-06  3.953835e-06
                         	 -5.104706e-04  8.538922e-06  3.953835e-06  1.102622e-03
                        	] rtol = 1e-6
    	# Using Stata's value, R's plm has confint incorrect
    	@test confint(model) ≈ [-0.0417325 0.0324457
                            	-0.0072991 0.0002169
                            	-0.0014674 0.0031635
                             	 0.0070191 0.1390408
                           	   ] rtol = 1e-5
    	# Within
    	model = fit(EconometricModel,
                	@formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(County)),
                	data)
		@test sprint(show, model) == "Continuous Response Model\nNumber of observations: 630\nNull Loglikelihood: 1633.30\nLoglikelihood: 2273.95\nR-squared: 0.8707\nWald: 0.16 ∼ F(3, 537) ⟹ Pr > F = 0.9258\nFormula: CRMRTE ~ 1 + PrbConv + AvgSen + PrbPris + absorb(County)\nVariance Covariance Estimator: OIM\n──────────────────────────────────────────────────────────────────────────────────────\n                   PE           SE        t-value  Pr > |t|         2.50%       97.50%\n──────────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.0314527   0.00204638   15.3699       <1e-43   0.0274328    0.0354726  \nPrbConv       6.65981e-6  0.0001985     0.0335507    0.9732  -0.000383272  0.000396592\nAvgSen        7.83181e-5  0.000127904   0.612318     0.5406  -0.000172936  0.000329572\nPrbPris      -0.0013419   0.00405182   -0.331185     0.7406  -0.00930126   0.00661746 \n──────────────────────────────────────────────────────────────────────────────────────"
    	model = fit(EconometricModel,
                	@formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(County + Year)),
                	data)
		@test sprint(show, model) == "Continuous Response Model\nNumber of observations: 630\nNull Loglikelihood: 1633.30\nLoglikelihood: 2290.25\nR-squared: 0.8775\nWald: 0.23 ∼ F(3, 531) ⟹ Pr > F = 0.8767\nFormula: CRMRTE ~ 1 + PrbConv + AvgSen + PrbPris + absorb(County + Year)\nVariance Covariance Estimator: OIM\n──────────────────────────────────────────────────────────────────────────────────────\n                    PE           SE       t-value  Pr > |t|         2.50%       97.50%\n──────────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.0319651    0.00206378   15.4886      <1e-44   0.027911     0.0360193  \nPrbConv       7.90362e-5   0.000196089   0.403062    0.6871  -0.00030617   0.000464242\nAvgSen       -9.4788e-5    0.000134924  -0.702531    0.4827  -0.000359838  0.000170262\nPrbPris       0.000979533  0.0040432     0.242267    0.8087  -0.00696309   0.00892216 \n──────────────────────────────────────────────────────────────────────────────────────"
    	# Random Effects
    	model = fit(RandomEffectsEstimator,
                	@formula(CRMRTE ~ PrbConv + AvgSen + PrbPris),
                	data,
					panel = :County,
					time = :Year)
		@test sprint(show, model) == "One-way Random Effect Model\nLongitudinal dataset: County, Year\nBalanced dataset with 90 panels of length 7\nindividual error component: 0.0162\nidiosyncratic error component: 0.0071\nρ: 0.8399\nNumber of observations: 630\nNull Loglikelihood: 2225.79\nLoglikelihood: 2226.01\nR-squared: 0.0007\nWald: 0.15 ∼ F(3, 626) ⟹ Pr > F = 0.9291\nFormula: CRMRTE ~ PrbConv + AvgSen + PrbPris\nVariance Covariance Estimator: OIM\n──────────────────────────────────────────────────────────────────────────────────────\n                    PE           SE       t-value  Pr > |t|         2.50%       97.50%\n──────────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.0309293    0.00266706   11.5968      <1e-27   0.0256918    0.0361668  \nPrbConv      -3.96634e-5   0.000198471  -0.199845    0.8417  -0.000429413  0.000350086\nAvgSen        8.2428e-5    0.000127818   0.644887    0.5192  -0.000168575  0.000333431\nPrbPris      -0.000123327  0.00404272   -0.030506    0.9757  -0.00806226   0.00781561 \n──────────────────────────────────────────────────────────────────────────────────────"
		# IV Pooling
    	model = fit(EconometricModel,
                	@formula(CRMRTE ~ PrbConv + (AvgSen ~ PrbPris)),
                	data)
		@test sprint(show, model) == "Continuous Response Model\nNumber of observations: 630\nNull Loglikelihood: 1633.30\nLoglikelihood: -611.40\nR-squared: NaN\nLR Test: -4489.40 ∼ χ²(2) ⟹  Pr > χ² = 1.0000\nFormula: CRMRTE ~ 1 + PrbConv + (AvgSen ~ PrbPris)\nVariance Covariance Estimator: OIM\n─────────────────────────────────────────────────────────────────────────────────\n                   PE          SE        t-value  Pr > |t|       2.50%     97.50%\n─────────────────────────────────────────────────────────────────────────────────\n(Intercept)   2.17878     23.0244      0.0946291    0.9246  -43.0357    47.3932  \nPrbConv       0.00456765   0.0638118   0.0715801    0.9430   -0.120743   0.129879\nAvgSen       -0.240139     2.57602    -0.093221     0.9258   -5.29883    4.81855 \n─────────────────────────────────────────────────────────────────────────────────"
    	# IV Between
    	model = fit(BetweenEstimator,
                	@formula(CRMRTE ~ PrbConv + (AvgSen ~ PrbPris)),
                	data,
					panel = :County)
		@test sprint(show, model) == "Between Estimator\nCounty with 630 groups\nBalanced groups with size 7\nNumber of observations: 90\nNull Loglikelihood: 239.56\nLoglikelihood: 161.00\nR-squared: NaN\nWald: 0.73 ∼ F(2, 86) ⟹ Pr > F = 0.4844\nFormula: CRMRTE ~ 1 + PrbConv + (AvgSen ~ PrbPris)\nVariance Covariance Estimator: OIM\n─────────────────────────────────────────────────────────────────────────────────\n                   PE          SE       t-value  Pr > |t|       2.50%      97.50%\n─────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.250667    0.255826     0.979834    0.3299  -0.257898   0.759233  \nPrbConv      -0.00331719  0.00481737  -0.688589    0.4929  -0.0128938  0.00625942\nAvgSen       -0.0242107   0.0286331   -0.845549    0.4002  -0.0811314  0.03271   \n─────────────────────────────────────────────────────────────────────────────────"
    	# IV Within
    	model = fit(EconometricModel,
	 				@formula(CRMRTE ~ PrbConv + (AvgSen ~ PrbPris) + absorb(County)),
                	data)
		@test sprint(show, model) == "Continuous Response Model\nNumber of observations: 630\nNull Loglikelihood: 1633.30\nLoglikelihood: 2245.89\nR-squared: NaN\nWald: 0.04 ∼ F(2, 537) ⟹ Pr > F = 0.9583\nFormula: CRMRTE ~ 1 + PrbConv + (AvgSen ~ PrbPris) + absorb(County)\nVariance Covariance Estimator: OIM\n──────────────────────────────────────────────────────────────────────────────────────\n                    PE           SE       t-value  Pr > |t|         2.50%       97.50%\n──────────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.0389703    0.025508      1.52777     0.1272  -0.0111374    0.089078   \nPrbConv       2.46646e-5   0.000215805   0.114291    0.9090  -0.000399261  0.000448591\nAvgSen       -0.000826364  0.00285293   -0.289654    0.7722  -0.00643064   0.00477791 \n──────────────────────────────────────────────────────────────────────────────────────"
    	model = fit(EconometricModel,
                	@formula(CRMRTE ~ PrbConv + (AvgSen ~ PrbPris) + absorb(County + Year)),
					data)
		@test sprint(show, model) == "Continuous Response Model\nNumber of observations: 630\nNull Loglikelihood: 1633.30\nLoglikelihood: 2280.31\nR-squared: NaN\nWald: 0.09 ∼ F(2, 531) ⟹ Pr > F = 0.9116\nFormula: CRMRTE ~ 1 + PrbConv + (AvgSen ~ PrbPris) + absorb(County + Year)\nVariance Covariance Estimator: OIM\n────────────────────────────────────────────────────────────────────────────────────\n                   PE           SE      t-value  Pr > |t|         2.50%       97.50%\n────────────────────────────────────────────────────────────────────────────────────\n(Intercept)  0.027409     0.0208187    1.31656     0.1886  -0.0134881    0.0683061  \nPrbConv      6.00404e-5   0.000214966  0.279301    0.7801  -0.000362248  0.000482329\nAvgSen       0.000462028  0.0023309    0.198219    0.8429  -0.00411688   0.00504094 \n────────────────────────────────────────────────────────────────────────────────────"
    	# IV Random
    	model = fit(RandomEffectsEstimator,
                	@formula(CRMRTE ~ PrbConv + (AvgSen ~ PrbPris)),
                	data,
					panel = :County,
					time = :Year)
		@test sprint(show, model) == "One-way Random Effect Model\nLongitudinal dataset: County, Year\nBalanced dataset with 90 panels of length 7\nindividual error component: 0.0413\nidiosyncratic error component: 0.0074\nρ: 0.9691\nNumber of observations: 630\nNull Loglikelihood: 2268.01\nLoglikelihood: 2248.34\nR-squared: NaN\nWald: 0.03 ∼ F(2, 626) ⟹ Pr > F = 0.9671\nFormula: CRMRTE ~ PrbConv + (AvgSen ~ PrbPris)\nVariance Covariance Estimator: OIM\n───────────────────────────────────────────────────────────────────────────────────────\n                    PE           SE        t-value  Pr > |t|         2.50%       97.50%\n───────────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.037747     0.0241684     1.56183      0.1188  -0.00971405   0.085208   \nPrbConv       1.39581e-5   0.000200244   0.0697054    0.9445  -0.000379273  0.000407189\nAvgSen       -0.000688924  0.00266514   -0.258494     0.7961  -0.00592262   0.00454478 \n───────────────────────────────────────────────────────────────────────────────────────"
  	end
	@testset "Heteroscedasticity Consistent Variance Covariance Estimators" begin
		data = dataset("Ecdat", "Crime")
    	# Pooling
    	model = fit(EconometricModel,
                	@formula(CRMRTE ~ PrbConv + AvgSen + PrbPris),
                	data)
		@test vcov(model) ≈ vcov(model, OIM)
		@test vcov(model, HC0) ≈ [ 1.51708E-05  3.58518E-07 -7.30083E-07 -2.03240E-05
									 3.58518E-07  2.85571E-07 -3.58783E-08 -4.13434E-07
									-7.30083E-07 -3.58783E-08  9.36232E-08 -1.46277E-07
									-2.03240E-05 -4.13434E-07 -1.46277E-07  5.17390E-05
								   ] rtol = 1e-6
		@test vcov(model, HC1) ≈ [ 1.52677E-05  3.60808E-07 -7.34748E-07 -2.04539E-05
		   							 3.60808E-07  2.87396E-07 -3.61076E-08 -4.16076E-07
									-7.34748E-07 -3.61076E-08  9.42214E-08 -1.47212E-07
									-2.04539E-05 -4.16076E-07 -1.47212E-07  5.20696E-05
								   ] rtol = 1e-6
		@test vcov(model, HC2) ≈ [ 1.56208E-05  9.05631E-08 -7.36000E-07 -2.08761E-05
									 9.05631E-08  6.78782E-07 -5.03354E-08 -3.93296E-08
									-7.36000E-07 -5.03354E-08  9.59715E-08 -1.59383E-07
									-2.08761E-05 -3.93296E-08 -1.59383E-07  5.27586E-05
								   ] rtol = 1e-6
		@test vcov(model, HC3) ≈ [ 1.66839E-05 -9.88963E-07 -7.13420E-07 -2.23048E-05
									-9.88963E-07  2.16675E-06 -1.04165E-07  1.50077E-06
									-7.13420E-07 -1.04165E-07  9.98145E-08 -2.14382E-07
									-2.23048E-05  1.50077E-06 -2.14382E-07  5.50463E-05
								   ] rtol = 1e-6
		@test vcov(model, HC4) ≈ [ 3.21597E-05 -2.14119E-05 -4.18897E-09 -4.43245E-05
									-2.14119E-05  2.97197E-05 -1.09486E-06  3.08355E-05
									-4.18897E-09 -1.09486E-06  1.38192E-07 -1.26498E-06
									-4.43245E-05  3.08355E-05 -1.26498E-06  8.67924E-05
								   ] rtol = 1e-6
		dv = fit(EconometricModel,
				 @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + County),
				 data,
				 contrasts = Dict(:County => DummyCoding()))
		fe = fit(EconometricModel,
				 @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(County)),
				 data)
		@test vcov(fe)[2:end,2:end] ≈ vcov(dv)[2:4, 2:4]
		@test vcov(fe, HC1)[2:end,2:end] ≈ vcov(dv, HC1)[2:4, 2:4]
	end
end
# Nominal Models
@testset "Nominal Models" begin
	data = joinpath(dirname(pathof(Econometrics)), "..", "data", "insure.csv") |>
		CSV.read |>
		(data -> select(data, [:insure, :age, :male, :nonwhite, :site])) |>
		dropmissing |>
		(data -> categorical!(data, [:insure, :site]))
	model = fit(EconometricModel,
				@formula(insure ~ age + male + nonwhite + site),
				data,
				contrasts = Dict(:insure => DummyCoding(base = "Uninsure")))
	@test sprint(show, model) == "Probability Model for Nominal Response\nCategories: Uninsure, Indemnity, Prepaid\nNumber of observations: 615\nNull Loglikelihood: -555.85\nLoglikelihood: -534.36\nR-squared: 0.0387\nLR Test: 42.99 ∼ χ²(10) ⟹ Pr > χ² = 0.0000\nFormula: insure ~ 1 + age + male + nonwhite + site\n───────────────────────────────────────────────────────────────────────────────────────────────────\n                                       PE         SE       t-value  Pr > |t|       2.50%     97.50%\n───────────────────────────────────────────────────────────────────────────────────────────────────\ninsure: Indemnity ~ (Intercept)   1.28694     0.59232     2.17271     0.0302   0.123689   2.4502   \ninsure: Indemnity ~ age           0.00779612  0.0114418   0.681372    0.4959  -0.0146743  0.0302666\ninsure: Indemnity ~ male         -0.451848    0.367486   -1.22957     0.2193  -1.17355    0.269855 \ninsure: Indemnity ~ nonwhite     -0.217059    0.425636   -0.509965    0.6103  -1.05296    0.618843 \ninsure: Indemnity ~ site: 2       1.21152     0.470506    2.57493     0.0103   0.287497   2.13554  \ninsure: Indemnity ~ site: 3       0.207813    0.366293    0.56734     0.5707  -0.511547   0.927172 \ninsure: Prepaid ~ (Intercept)     1.55666     0.596327    2.61041     0.0093   0.385533   2.72778  \ninsure: Prepaid ~ age            -0.00394887  0.0115993  -0.340439    0.7336  -0.0267287  0.018831 \ninsure: Prepaid ~ male            0.109846    0.365187    0.300793    0.7637  -0.607343   0.827035 \ninsure: Prepaid ~ nonwhite        0.757718    0.419575    1.80592     0.0714  -0.0662835  1.58172  \ninsure: Prepaid ~ site: 2         1.32456     0.469789    2.81947     0.0050   0.401941   2.24717  \ninsure: Prepaid ~ site: 3        -0.380175    0.372819   -1.01973     0.3083  -1.11235    0.352001 \n───────────────────────────────────────────────────────────────────────────────────────────────────"
	β, V, σ = coef(model), vcov(model), stderror(model)
	@test β ≈ [1.286943, 0.0077961, -0.4518496, -0.2170589, 1.211563, 0.2078123,
               1.556656, -0.0039489, 0.1098438, 0.7577178, 1.324599, -0.3801756
              ] rtol = 1e-3
  	@test V ≈ [ 0.35084518 -0.00590456 -0.02519209 -0.01837512 -0.08995794 -0.08264268  0.29928937 -0.00510107 -0.02257270 -0.01585492 -0.07392512 -0.06808314
               -0.00590456  0.00013092 -0.00046225 -0.00028066  0.00041004  0.00040470 -0.00509454  0.00011354 -0.00039387 -0.00024122  0.00033847  0.00035236
               -0.02519209 -0.00046225  0.13504644  0.01565749  0.01009559  0.00933967 -0.02288847 -0.00039497  0.11365141  0.01369872  0.00945927  0.00802593
               -0.01837512 -0.00590456 -0.02519209  0.18116609  0.01886597 -0.02821348 -0.01638629 -0.00023179  0.01314666  0.15068111  0.01664999 -0.02345786
               -0.08995794 -0.00590456 -0.02519209 -0.01837512  0.22138224 0.064050660 -0.07438947  0.00035751  0.00845363  0.01594691  0.19895500  0.05189742
               -0.08264268  0.00040470  0.00933967 -0.02821348  0.06405066 0.134170290 -0.06796097  0.00034828  0.00826890 -0.02353016  0.05182727  0.11060487
                0.29928937 -0.00509454 -0.02288847 -0.01638629 -0.07438947 -0.06796097  0.35560782 -0.00601966 -0.02683161 -0.02105510 -0.08991231 -0.08165509
               -0.00510107  0.00011354 -0.00039497 -0.00023179  0.00035751  0.00034828 -0.00601966  0.00013455 -0.00047223 -0.00028680  0.00039990  0.0004179
               -0.02257270 -0.00039387  0.11365141  0.01314666  0.00845363  0.00826890 -0.02683161 -0.00047223  0.13336253  0.01625685  0.01101870  0.00922333
               -0.01585492 -0.00024122  0.01369872  0.15068111  0.01594691 -0.02353016 -0.02105510 -0.00028680  0.01625685  0.17604391  0.02099458 -0.03008391
               -0.07392512  0.00033847  0.00945927  0.01664999  0.19895500  0.05182727 -0.08991231  0.00039990  0.01101870  0.02099458  0.22070774  0.06273376
               -0.06808314  0.00035236  0.00802593 -0.02345786  0.05189742  0.11060487 -0.08165509  0.00041790  0.00922333 -0.03008391  0.06273376  0.13899387
              ] |> Hermitian rtol = 1e-3
	@test σ ≈ [0.5923219, 0.0114418, 0.3674867, 0.4256361, 0.4705127, 0.3662926,
  			   0.5963286, 0.0115994, 0.3651883, 0.4195759, 0.4697954, 0.3728188
			  ] rtol = 1e-4
end
# Ordinal Models
@testset "Ordinal Models" begin
	data = dataset("Ecdat", "Kakadu") |>
		(data -> select(data, [:RecParks, :Sex, :Age, :Schooling]))
	data.RecParks = convert(Vector{Int}, data.RecParks)
	data.RecParks = levels!(categorical(data.RecParks, ordered = true), collect(1:5))
	model = fit(EconometricModel,
              	@formula(RecParks ~ Age + Sex + Schooling),
              	data)
  	@test sprint(show, model) == "Probability Model for Ordinal Response\nCategories: 1 < 2 < 3 < 4 < 5\nNumber of observations: 1827\nNull Loglikelihood: -2677.60\nLoglikelihood: -2657.40\nR-squared: 0.0075\nLR Test: 40.42 ∼ χ²(3) ⟹ Pr > χ² = 0.0000\nFormula: RecParks ~ Age + Sex + Schooling\n──────────────────────────────────────────────────────────────────────────────────────────\n                          PE          SE        t-value  Pr > |t|        2.50%      97.50%\n──────────────────────────────────────────────────────────────────────────────────────────\nAge                  0.00943647  0.00254031    3.71469     0.0002   0.00445424   0.0144187\nSex: male           -0.0151659   0.0846365    -0.179188    0.8578  -0.181161     0.150829 \nSchooling           -0.103902    0.0248742    -4.17711     <1e-4   -0.152687    -0.0551174\n(Intercept): 1 | 2  -2.92405     0.191927    -15.2352      <1e-48  -3.30047     -2.54763  \n(Intercept): 2 | 3  -1.54922     0.171444     -9.03632     <1e-18  -1.88547     -1.21297  \n(Intercept): 3 | 4  -0.298938    0.166904     -1.79108     0.0734  -0.626281     0.0284051\n(Intercept): 4 | 5   0.669835    0.167622      3.9961      <1e-4    0.341083     0.998587 \n──────────────────────────────────────────────────────────────────────────────────────────"
  	β, V, σ = coef(model), vcov(model), stderror(model)
  	@test β ≈ [ 0.009437926, -0.015143049, -0.103911316,
               -2.92391240,  -1.549266200, -0.298963900, 0.6698249] rtol = 1e-4
	@test V ≈ [ 0.00000646  -0.0000007000  0.000013250  0.00031473  0.0003161358  0.000321276  0.0003270734
               -0.00000070   0.0071633800 -0.000184930  0.00294759  0.0029647442  0.002995366  0.0029904485
                0.00001325  -0.0001849300  0.000618730  0.00288365  0.0028548209  0.002791361  0.0027318928
                0.00031473   0.0029475900  0.002883650  0.03684071  0.0289239290  0.026986237  0.026498765
                0.00031614   0.0029647400  0.002854820  0.02892393  0.0293971651  0.027160433  0.0266006675
                0.000321276  0.0029953660  0.002791361  0.02698623  0.0271604331  0.027860143  0.0270590982
                0.0003270734 0.002990448   0.002731893  0.02649876  0.0266006675  0.027059098  0.0281002038
              ] rtol = 1e-3
    @test σ ≈ [0.002540789, 0.084636770, 0.024874381,
               0.191939349, 0.171456015, 0.166913580, 0.167631154] rtol = 1e-4
	data = data[data.RecParks .≠ 5,:]
	model = fit(EconometricModel,
              	@formula(RecParks ~ Age + Sex + Schooling),
              	data)
	@test sprint(show, model) == "Probability Model for Ordinal Response\nCategories: 1 < 2 < 3 < 4\nNumber of observations: 1200\nNull Loglikelihood: -1502.61\nLoglikelihood: -1500.73\nR-squared: 0.0013\nLR Test: 3.76 ∼ χ²(3) ⟹ Pr > χ² = 0.2886\nFormula: RecParks ~ Age + Sex + Schooling\n───────────────────────────────────────────────────────────────────────────────────────────\n                          PE          SE        t-value  Pr > |t|        2.50%       97.50%\n───────────────────────────────────────────────────────────────────────────────────────────\nAge                  0.00063565  0.00329173    0.193105    0.8469  -0.00582255   0.00709385\nSex: male            0.0896101   0.106068      0.84484     0.3984  -0.118489     0.297709  \nSchooling           -0.0509934   0.0300878    -1.69482     0.0904  -0.110024     0.00803731\n(Intercept): 1 | 2  -2.58123     0.228989    -11.2723      <1e-27  -3.0305      -2.13197   \n(Intercept): 2 | 3  -1.12053     0.212192     -5.28071     <1e-6   -1.53684     -0.704217  \n(Intercept): 3 | 4   0.486702    0.210541      2.31167     0.0210   0.0736316    0.899772  \n───────────────────────────────────────────────────────────────────────────────────────────"
  	β, V, σ = coef(model), vcov(model), stderror(model)
  	@test β ≈ [ 0.00063565,  0.08961012, -0.05099339,
               -2.5812334,  -1.1205278,   0.48670172] rtol = 1e-4
	@test V ≈ [ 0.00001084 -5.386e-06   0.00002186 0.00052578 0.0005259  0.00052967
  	 		   -5.386e-06   0.01125032 -0.00021416 0.00465658 0.00470423 0.00484966
  	  			0.00002186 -0.00021416  0.00090528 0.00435314 0.00431928 0.00424373
  	  			0.00052578  0.00465658  0.00435314 0.05243606 0.04424607 0.04221716
  	  			0.0005259   0.00470423  0.00431928 0.04424607 0.04502565 0.0424849
  	  			0.00052967  0.00484966  0.00424373 0.04221716 0.0424849  0.04432745
              ] rtol = 1e-3
    @test σ ≈ [0.0032917, 0.1060675, 0.0300878,
			   0.2289892, 0.2121925, 0.2105408] rtol = 1e-4
  	data = joinpath(dirname(pathof(Econometrics)), "..", "data", "auto.csv") |>
		CSV.read |>
    	(data -> select(data, [:rep77, :foreign, :length, :mpg])) |>
    	dropmissing |>
		categorical!
  	data.rep77 = levels!(categorical(data.rep77; ordered = true),
                       	 ["Poor", "Fair", "Average", "Good", "Excellent"])
  	model = fit(EconometricModel,
              	@formula(rep77 ~ foreign + length + mpg),
              	data)
  	@test sprint(show, model) == "Probability Model for Ordinal Response\nCategories: Poor < Fair < Average < Good < Excellent\nNumber of observations: 66\nNull Loglikelihood: -89.90\nLoglikelihood: -78.25\nR-squared: 0.1295\nLR Test: 23.29 ∼ χ²(3) ⟹ Pr > χ² = 0.0000\nFormula: rep77 ~ foreign + length + mpg\n─────────────────────────────────────────────────────────────────────────────────────────────\n                                    PE         SE     t-value  Pr > |t|       2.50%    97.50%\n─────────────────────────────────────────────────────────────────────────────────────────────\nforeign: Foreign                2.89681    0.790641   3.66387    0.0005   1.31684     4.47678\nlength                          0.0828275  0.02272    3.64558    0.0005   0.0374253   0.12823\nmpg                             0.230768   0.0704548  3.2754     0.0017   0.0899749   0.37156\n(Intercept): Poor | Fair       17.9275     5.55119    3.22948    0.0020   6.83431    29.0206 \n(Intercept): Fair | Average    19.8651     5.59648    3.54956    0.0007   8.68139    31.0487 \n(Intercept): Average | Good    22.1033     5.70894    3.8717     0.0003  10.6949     33.5117 \n(Intercept): Good | Excellent  24.6921     5.89075    4.19168    <1e-4   12.9204     36.4639 \n─────────────────────────────────────────────────────────────────────────────────────────────"
  	β, V, σ = coef(model), vcov(model), stderror(model)
  	@test β ≈ [2.89679875, 0.08282676, 0.23076532, 17.92728, 19.86486, 22.10311, 24.69193] atol = 1e-3
  	@test V ≈ [ 0.62578335	0.0111044038 0.012321411 2.4421745	2.4836507	2.5750058	2.6957587
             	0.0111044	0.0005186744 0.001222385 0.1243671	0.1258979	0.1284727	0.1319772
                0.01232141	0.0012223847 0.004974678 0.3328169	0.3368057	0.3430257	0.3549711
                2.44217451	0.1243670884 0.332816913 30.9544666	31.0279194	31.6005528	32.5099187
                2.48365067	0.1258979463 0.336805731 31.0279194	31.4396158	31.9986912	32.9170873
                2.57500582	0.1284726666 0.343025653 31.6005528	31.9986912	32.7044311	33.6189243
                2.69575871	0.1319772302 0.354971081 32.5099187	32.9170873	33.6189243	34.8533724
              ] rtol = 1e-2
    @test σ ≈ [0.79106469, 0.02277442, 0.07053140,
               5.56367384, 5.60710405, 5.71877881, 5.90367448] rtol = 1e-2
end
# Formula displays
@testset "Formula Displays" begin
	data = dataset("Ecdat", "Crime")
	# Continuous Response
	f = @formula(CRMRTE ~ 1 + PrbConv)
	model = fit(EconometricModel, f, data)
	@test Econometrics.clean_fm(model) == "Formula: CRMRTE ~ 1 + PrbConv"
	f = @formula(CRMRTE ~ 0 + PrbConv)
	model = fit(EconometricModel, f, data)
	@test Econometrics.clean_fm(model) == "Formula: CRMRTE ~ 0 + PrbConv"
	f = @formula(CRMRTE ~ -1 + PrbConv)
	model = fit(EconometricModel, f, data)
	@test Econometrics.clean_fm(model) == "Formula: CRMRTE ~ -1 + PrbConv"
	f = @formula(CRMRTE ~ -1 + PrbConv + PrbConv & PrbPris^2 + absorb(County + Year))
	model = fit(EconometricModel, f, data)
	@test Econometrics.clean_fm(model) ==
		"Formula: CRMRTE ~ -1 + PrbConv + absorb(County + Year) + PrbConv & (PrbPris ^ 2)"
	f = @formula(log(CRMRTE) ~ 1)
	model = fit(EconometricModel, f, data)
	@test Econometrics.clean_fm(model) ==
		"Formula: log(CRMRTE) ~ 1"
	# Random Effects Model
	f = @formula(CRMRTE ~ 1 + PrbConv)
	model = fit(RandomEffectsEstimator, f, data, panel = :County, time = :Year)
	@test Econometrics.clean_fm(model) == "Formula: CRMRTE ~ PrbConv"
	# Nominal Response Model
	data = joinpath(dirname(pathof(Econometrics)), "..", "data", "insure.csv") |>
		CSV.read |>
		(data -> select(data, [:insure, :age, :male, :nonwhite, :site])) |>
		dropmissing |>
		(data -> categorical!(data, [:insure, :site]))
	f = @formula(insure ~ 0 + age^2 + male + nonwhite + site)
	model = fit(EconometricModel, f, data,
				contrasts = Dict(:insure => DummyCoding(base = "Uninsure")))
	@test Econometrics.clean_fm(model) == "Formula: insure ~ 0 + (age ^ 2) + male + nonwhite + site"
	# Ordinal Response Model
	data = dataset("Ecdat", "Kakadu") |>
		(data -> select(data, [:RecParks, :Sex, :Age, :Schooling]))
	data.RecParks = convert(Vector{Int}, data.RecParks)
	data.RecParks = levels!(categorical(data.RecParks, ordered = true), collect(1:5))
	model = fit(EconometricModel,
				@formula(RecParks ~ Age * Age + Sex),
              	data)
	@test Econometrics.clean_fm(model) == "Formula: RecParks ~ Age + Sex + Age & Age"
end


text:
using Revise
using Mitosis
using Random, Test, LinearAlgebra, Statistics
using Distributions, StatsBase
using Mitosis: outer
const d = 4

# Forward model and data generation: multivariate Wiener
# X = x0 + μ*t + σ*W(t)

μ = randn(d)
σ = randn(d,d)
Σ = outer(σ)
n = 500
t = [0.0; sort(1000*rand(n-1))]
n = length(t)
function bm(x, s, μ, σ)
    y = float(x)
    ys = [y]
    for i in 1:length(s)-1
        Δ = s[i+1] - s[i]
        k = kernel(Gaussian; μ=AffineMap(I(d), Δ*μ), Σ=ConstantMap(Σ*Δ))
        y = rand(k(y))
        push!(ys, y)
    end
    ys
end
x0 = zeros(d)
x = bm(x0, t, μ, σ)

# Conjugate analysis of increments

using Distributions
E(i, d) = Vector(I(d)[:,i])

function gibbs(t, x, (μ0, Σ0), πμ, πΣ, iterations)
    n = length(t)
    d = length(μ0)
    μ = μ0
    Σ = Σ0
    out = [(μ, σ)]
    for iter in 1:iterations
        # conjugate μ
        p = πμ
        for i in n-1:-1:1
            Δ = t[i+1] - t[i]
            r = Mitosis.Regression(μ, [ConstantMap(Δ*E(i, d)) for i in 1:d], identity)
            k = kernel(Gaussian; μ=r, Σ=ConstantMap(Σ*Δ))
            con = Mitosis.conjugate(k, x[i])
            _, pi = backward(BF(), con, x[i+1])
            _, p = fuse(p, pi)
        end
        μ = rand(convert(Gaussian{(:μ,:Σ)}, p))

        α = params(πΣ)[1]
        β = Matrix(params(πΣ)[2])

        for i in n-1:-1:1
            Δ = t[i+1] - t[i]
            α = α + 1/2
            k = kernel(Gaussian; μ=AffineMap(I(d), Δ*μ), Σ=ConstantMap(Σ*Δ))
            β = β + Mitosis.outer(x[i+1] .- mean(k(x[i])))/(2Δ)
        end
        Σ = rand(InverseWishart(α, β))

        push!(out, (μ, Σ))
    end
    out
end

F0, Γ0 = zeros(d), zeros(d,d)
πμ = Gaussian{(:F,:Γ)}(F0, Γ0)
πΣ = InverseWishart(d+4, zeros(d,d) + 20I)
var(πΣ)
Σ0 = mean(πΣ)
chain = gibbs(t, x, (F0, Σ0), πμ, πΣ, 2000)
μs = first.(chain)
Σs = last.(chain)

display([μ mean(μs)])
display([vec(Σ) vec(mean(Σs))])


text:
"""
    Vector{UInt8} -> String

Returns a String representing a Bitcoin address
"""
function h160_2_address(h160::Vector{UInt8}, testnet::Bool=false, type::String="P2SH")
    testnet ? i = 1 : i = 2
    result = copy(h160)
    result = pushfirst!(result, SCRIPT_TYPES[type][i])
    return String(base58checkencode(result))
end

"""
    point2address(P::Secp256k1.Point, compressed::Bool, testnet::Bool) -> String

Returns the Base58 Bitcoin address given an Secp256k1.Point
Compressed is set to true if not provided.
Testnet is set to false by default.
"""
function point2address(P::T, compressed::Bool=true, testnet::Bool=false, type::String="P2PKH") where {T<:Secp256k1.Point}
    s = serialize(P, compressed=compressed)
    h160 = ripemd160(sha256(s))
    return h160_2_address(h160, testnet, type)
end

"""
    wif(kp::KeyPair, compressed::Bool=true, testnet::Bool=false) -> String

Returns a KeyPair in Wallet Import Format
Compressed is set to true if not provided.
Testnet is set to false by default.
"""
function wif(kp::KeyPair, compressed::Bool=true, testnet::Bool=false)
    secret_bytes = bytes(kp.𝑑)
    if testnet
        prefix = 0xef
    else
        prefix = 0x80
    end
    result = pushfirst!(secret_bytes, prefix)
    if compressed
        return String(base58checkencode(push!(result, 0x01)))
    else
        return String(base58checkencode(result))
    end
end

@deprecate address(P::T, compressed::Bool, testnet::Bool)  where {T<:Secp256k1.Point} point2address(P::T, compressed::Bool, testnet::Bool, type::String)  where {T<:Secp256k1.Point}


text:
export iauRefco
"""
Determine the constants A and B in the atmospheric refraction model
dZ = A tan Z + B tan^3 Z.

Z is the "observed" zenith distance (i.e. affected by refraction)
and dZ is what to add to Z to give the "topocentric" (i.e. in vacuo)
zenith distance.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   phpa   double    pressure at the observer (hPa = millibar)
   tc     double    ambient temperature at the observer (deg C)
   rh     double    relative humidity at the observer (range 0-1)
   wl     double    wavelength (micrometers)

Returned:
   refa   double*   tan Z coefficient (radians)
   refb   double*   tan^3 Z coefficient (radians)

Notes:

   1. The model balances speed and accuracy to give good results in
   applications where performance at low altitudes is not paramount.
   Performance is maintained across a range of conditions, and
   applies to both optical/IR and radio.

   2. The model omits the effects of (i) height above sea level (apart
   from the reduced pressure itself), (ii) latitude (i.e. the
   flattening of the Earth), (iii) variations in tropospheric lapse
   rate and (iv) dispersive effects in the radio.

   The model was tested using the following range of conditions:

      lapse rates 0.0055, 0.0065, 0.0075 deg/meter
      latitudes 0, 25, 50, 75 degrees
      heights 0, 2500, 5000 meters ASL
      pressures mean for height -10% to +5% in steps of 5%
      temperatures -10 deg to +20 deg with respect to 280 deg at SL
      relative humidity 0, 0.5, 1
      wavelengths 0.4, 0.6, ... 2 micron, + radio
      zenith distances 15, 45, 75 degrees

   The accuracy with respect to raytracing through a model
   atmosphere was as follows:

                           worst         RMS

      optical/IR           62 mas       8 mas
      radio               319 mas      49 mas

   For this particular set of conditions:

      lapse rate 0.0065 K/meter
      latitude 50 degrees
      sea level
      pressure 1005 mb
      temperature 280.15 K
      humidity 80%
      wavelength 5740 Angstroms

   the results were as follows:

      ZD       raytrace     iauRefco   Saastamoinen

      10         10.27        10.27        10.27
      20         21.19        21.20        21.19
      30         33.61        33.61        33.60
      40         48.82        48.83        48.81
      45         58.16        58.18        58.16
      50         69.28        69.30        69.27
      55         82.97        82.99        82.95
      60        100.51       100.54       100.50
      65        124.23       124.26       124.20
      70        158.63       158.68       158.61
      72        177.32       177.37       177.31
      74        200.35       200.38       200.32
      76        229.45       229.43       229.42
      78        267.44       267.29       267.41
      80        319.13       318.55       319.10

   deg        arcsec       arcsec       arcsec

   The values for Saastamoinen's formula (which includes terms
   up to tan^5) are taken from Hohenkerk and Sinclair (1985).

   3. A wl value in the range 0-100 selects the optical/IR case and is
   wavelength in micrometers.  Any value outside this range selects
   the radio case.

   4. Outlandish input parameters are silently limited to
   mathematically safe values.  Zero pressure is permissible, and
   causes zeroes to be returned.

   5. The algorithm draws on several sources, as follows:

      a. The formula for the saturation vapour pressure of water as
      a function of temperature and temperature is taken from
      Equations (A4.5-A4.7) of Gill (1982).

      b. The formula for the water vapour pressure, given the
      saturation pressure and the relative humidity, is from
      Crane (1976), Equation (2.5.5).

      c. The refractivity of air is a function of temperature,
      total pressure, water-vapour pressure and, in the case
      of optical/IR, wavelength.  The formulae for the two cases are
      developed from Hohenkerk & Sinclair (1985) and Rueger (2002).

      d. The formula for beta, the ratio of the scale height of the
      atmosphere to the geocentric distance of the observer, is
      an adaption of Equation (9) from Stone (1996).  The
      adaptations, arrived at empirically, consist of (i) a small
      adjustment to the coefficient and (ii) a humidity term for the
      radio case only.

      e. The formulae for the refraction constants as a function of
      n-1 and beta are from Green (1987), Equation (4.31).

References:

   Crane, R.K., Meeks, M.L. (ed), "Refraction Effects in the Neutral
   Atmosphere", Methods of Experimental Physics: Astrophysics 12B,
   Academic Press, 1976.

   Gill, Adrian E., "Atmosphere-Ocean Dynamics", Academic Press,
   1982.

   Green, R.M., "Spherical Astronomy", Cambridge University Press,
   1987.

   Hohenkerk, C.Y., & Sinclair, A.T., NAO Technical Note No. 63,
   1985.

   Rueger, J.M., "Refractive Index Formulae for Electronic Distance
   Measurement with Radio and Millimetre Waves", in Unisurv Report
   S-68, School of Surveying and Spatial Information Systems,
   University of New South Wales, Sydney, Australia, 2002.

   Stone, Ronald C., P.A.S.P. 108, 1051-1058, 1996.

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauRefco(double phpa, double tc, double rh, double wl,
#               double *refa, double *refb)

function iauRefco(phpa::Real, tc::Real, rh::Real, wl::Real)
   # Preallocate return values
   ref_refa = Ref{Float64}(0.0)
   ref_refb = Ref{Float64}(0.0)

   ccall((:iauRefco, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Cdouble, Cdouble,
         Ref{Cdouble}, Ref{Cdouble}),
          convert(Float64, phpa),
          convert(Float64, tc),
          convert(Float64, rh),
          convert(Float64, wl),
          ref_refa, ref_refb)

   return ref_refa[], ref_refb[]
end


text:
module YYZRecipes

# Write your package code here.

using Requires
using RecipesBase

function __init__()
    @require Distributions="31c24e10-a181-5473-b8eb-7969acd0382f" include("Distributions.jl")
end

end


text:
# Model permittivity of doped silicon, reproducing the model used in Fu, C. J.,
# & Zhang, Z. M. (2006). Nanoscale radiation heat transfer for silicon at different doping levels.
# (International Journal of Heat and Mass Transfer, 49(9-10)), and 1703-1718 American Journal of
# Physics 44, 353 (1976); doi: 10.1119/1.10194.
## Dependencies
using Roots
# Conversion factors are loaded from call through filmUtilities.jl
"""

	blkFD(arg::Float64)::Float64

Blakemore approximation of the Fermi-Dirac integral, accurate to with two percent of true values.
"""
function blkFD(arg::Float64)::Float64

	if arg < -5.0

		return exp(arg)
	elseif arg < 0.7

		return exp(arg) - 0.3 * exp(2.0 * arg) + 0.06 * exp(3.0 * arg)
	elseif arg < 1.1

		return exp(arg) * /(1.0, 1.0 + 0.27 * exp(arg))
	elseif arg < 1.3

		return /(4.0 * ^(arg, 1.5), 3.0 * sqrt(pi)) * (1.0 + /(1.15, arg^2))
	else

		return /(4.0 * ^(arg, 1.5), 3.0 * sqrt(pi)) * (1.0 + /(pi^2, 8.0 * arg^2))
	end
end
precompile(blkFD, (Float64, ))
"""

	siRspBck(enr::Float64, tmp::Float64)::ComplexF64

Temperature dependent background permittivity of intrinsic silicon. Loss is assumed to be
negligible in comparison to free carrier contributions.
"""
# Reproduction of the model reported in Lee, B. J., Zhang, Z. M., Early, E. A., DeWitt, D. P., &
# Tsai, B. K. (2005). Modeling radiative properties of silicon with coatings and comparison with
# reflectance measurements. Journal of thermophysics and heat transfer, 19(4), 558-565.
@inline function siRspBck(enr::Float64, tmp::Float64)::ComplexF64

	return ^(sqrt(4.565 + /(97.3, ^(3.648,2) - ^(/(enr * 1.24, muEv), 2))) + (-1.864 * ^(10, -4) + 5.394 * /(^(10, -3), ^(3.648,2) - ^(/(enr * 1.24, muEv), 2))) * tmp, 2) - /(^(hBEv * 1.0e13, 2), enr * (enr + im * 1.0e10 * hBEv))
end
precompile(siRspBck, (Float64, Float64))
"""

	siEnrGap(tmp::Float64)::Float64

Temperature dependent energy gap for intrinsic silicon in electron volts.
"""
@inline function siEnrGap(tmp::Float64)::Float64

	return 1.17 - /(0.000473 * tmp^2.0, tmp + 636.0)
end
precompile(siEnrGap, (Float64, ))
"""

	itrnV(tmp::Float64)::Float64

Number of intrinsic valence band carriers in units of cm^{-3} as a function of temperature in
Kelvin.
"""
@inline function itrnV(tmp::Float64)::Float64

	return 2.66 * ^(10.0, 19.0) * ^(/(tmp, 300), 1.5)
end
precompile(itrnV, (Float64, ))
"""

	itrnC(tmp::Float64)::Float64

Number of intrinsic conduction band carriers in units of cm^{-3} as a function of temperature in
Kelvin.
"""
@inline function itrnC(tmp::Float64)::Float64

	return 2.86 * ^(10.0, 19.0) * ^(/(tmp, 300), 1.5)
end
precompile(itrnC, (Float64, ))
"""

	function dnrCcn(enrFrm::Float64, tmp::Float64, acp::dptDsc)::Float64

Acceptor concentration function.
"""

@inline function acpCcn(enrFrm::Float64, tmp::Float64, acp::dptDsc)::Float64

	return acp.ccn * (1.0 - /(4.0, 4.0  + exp(/(enrFrm - acp.enr, blzK * tmp))))
end
precompile(acpCcn, (Float64, Float64, dptDsc))
"""

	function dnrCcn(enrFrm::Float64, tmp::Float64, dnr::dptDsc)::Float64

Donor concentration function.
"""
@inline function dnrCcn(enrFrm::Float64, tmp::Float64, dnr::dptDsc)::Float64

	return dnr.ccn * (1.0 - /(2.0, 2.0  + exp(/(siEnrGap(tmp) - dnr.enr - enrFrm, blzK * tmp))))
end
precompile(dnrCcn, (Float64, Float64, dptDsc))
"""

	frmZero(enrFrm::Float64, tmp::Float64, dnr::dptVal, acp::dptVal)::Float64

Net charge concentration as a function of supposed Fermi energy and temperature.
Used as an implicit equation for determining the correct Fermi energy.
"""
@inline function frmZero(enrFrm::Float64, tmp::Float64, dnr::dptDsc, acp::dptDsc)::Float64

	return itrnV(tmp) * blkFD(-/(enrFrm, blzK * tmp)) + /(dnr.ccn, 1.0 + 2.0 * exp(/(enrFrm - siEnrGap(tmp) + dnr.enr, blzK * tmp))) - itrnC(tmp) * blkFD(/(enrFrm - siEnrGap(tmp), blzK * tmp)) - /(acp.ccn, 1.0 + 4.0 * exp(/(acp.enr - enrFrm, blzK * tmp)))
end
precompile(frmZero, (Float64, Float64, dptDsc, dptDsc))
"""

	prmMSi(tmp::Float64, dnr::dptDsc, acp::dptDsc)::siDsc

Computation of model parameters for doped silicon, see preamble of source file for details.
"""
function prmMSi(tmp::Float64, dnr::dptDsc, acp::dptDsc)::siDsc

	# Self-consistent Fermi level.
	frmFnc(flvl) = frmZero(flvl, tmp, dnr, acp)
	enrFrm = find_zero(frmFnc, /(siEnrGap(tmp), 2.0))

	# Number of ionized acceptor and donor dopants (carriers), units cm^-3.
	cncDnr = dnrCcn(enrFrm, tmp, dnr)
	cncAcp = acpCcn(enrFrm, tmp, acp)

	# Room temperature total electron scattering time.
	stETR = (141.0 + /(19.5, 1.0 + ^(/(dnr.ccn, 1.3 * ^(10.0, 17.0)), 0.91))) * 1.0e-15
	# Room temperature total hole scattering time.
	stHTR = (10.0 + /(94.0, 1.0 + ^(/(acp.ccn, 1.9 * ^(10.0, 17.0)), 0.76))) * 1.0e-15

	# Room temperature electron-lattice scattering time.
	stELR = 2.23 * 1.0e-13
	# Room temperature hole-lattice scattering time.
	stHLR = 1.06 * 1.0e-13
	# Room temperature electron-impurity scattering time.
	stEIR = /(stETR * stELR, stELR - stETR)
	# Room temperature hole-impurity scattering time.
	stHIR = /(stHTR * stHLR, stHLR - stHTR)

	# Include temperature dependence scattering times.
	# Electrons
	stEI = stEIR * ^(/(tmp, 300), 1.5)
	stEL = stELR * ^(/(tmp, 300), -3.8)
	# Holes
	stHI = stHIR * ^(/(tmp, 300), 1.5)
	stHL = stHLR * ^(/(tmp, 300), -3.6)

	# Total donor and acceptor decay rates.
	dnrDcy = hBEv * (/(1.0, stEI) + /(1.0, stEL))
	acpDcy = hBEv * (/(1.0, stHI) + /(1.0, stHL))

	# Electron and hole effective masses.
	mEffE = 0.27
	mEffH = 0.37

	# Conversion factors for donor and acceptor contributions.
	dnrFac = /(cncDnr * hBEv^2 * ^(10.0, 9.0), mEffE * 0.5199895 * 55.26349406 * 0.01112265)
	acpFac = /(cncAcp * hBEv^2 * ^(10.0, 9.0), mEffH * 0.5199895 * 55.26349406 * 0.01112265)

	return siDsc(tmp, dnrFac, acpFac, dnrDcy, acpDcy)
end
precompile(prmMSi, (Float64, dptDsc, dptDsc))


text:
module ChoroplethMaps

using DataFrames
using Proj4
using Gadfly
using Compose

export mapify, choroplethmap, choroplethlayer, graticule, Provider, projstr

# TODO: To be able to layer the plot over a background tile, use an invocation like this:
# compose(plot_context(), bitmap("image/png", Array{UInt8}(readall("price.png")), 0, 0, 1, 1))
# Hack this in with underguide perhaps?

include("providers/provider.jl")
include("graticule.jl")
include("underguide.jl")

"""
`mapify(df::DataFrame, provider::Provider.AbstractProvider;
 key::Symbol=:GEOID, plotgroup::Symbol=:blank_symbol,
 projection::AbstractString="epsg:3857", keepcols=[])`

Takes in a DataFrame and returns a new DataFrame containing polygons pulled
from the `provider`'s geodata, suitable for plotting with Gadfly.

`df` is the input DataFrame which will be joined with the `provider` polygons,
using field `key` as the join key. `key` MUST exist in both `df` and
`provider`.

`plotgroup` is the `provider` field that differentiates unique shapes,
and need only be specfied if it differs from `key`. `projection` is
an EPSG numeric code specifying the desired projection of the choropleth map;
the default is EPSG:3857 "Web Mercator". `keepcols` specifies any additional
columns (as symbols) from the provider that should be kept in the returned
DataFrame

Example:

    df = DataFrame(NAME=["Alabama","Mississippi"], POP=[4849000, 2994000])
    mapframe = mapify(df, Provider.statesummary(), key=:NAME)
    # Then to plot with Gadfly:
    plot(mapframe, x=:CM_X, y=:CM_Y, group=:NAME, color=:POP,
      Geom.polygon(preserve_order=true, fill=true),
      Coord.cartesian(fixed=true))
"""
function mapify(df::DataFrame, provider::Provider.AbstractProvider;
                  key::Symbol=:GEOID,
                  plotgroup::Symbol=:blank_symbol,
                  projection::AbstractString="epsg:3857",
                  keepcols=[])

  polygons = DataFrame(CM_POLYGONS=Provider.polygons(provider))
  records = DataFrame(Provider.records(provider))
  if plotgroup == :blank_symbol
    plotgroup = key
  end
  dropallbut!(records, vcat(unique([key,plotgroup]), keepcols))
  pdata = hcat(records, polygons)

  joined = join(df, pdata, on=key, kind=:inner)
  joined = splitmultipolygons(joined, plotgroup)

  newdf = add_xyp_cols!(similar(joined, 0))
  delete!(newdf, :CM_POLYGONS)

  projfrom = Projection(projstr(Provider.projection(provider)))
  projto = Projection(projstr(projection))

  for row in eachrow(joined)
    for pt in row[:CM_POLYGONS]
      copyrow!(newdf, row, exclude=[:CM_POLYGONS])

      # TODO: check that transform short-circuits if projfrom == projto
      tpt = transform(projfrom, projto, pt)

      push!(newdf[:CM_X], tpt[1])
      push!(newdf[:CM_Y], tpt[2])
      push!(newdf[:CM_P], projection)
    end
  end

  return newdf
end

"""
`choroplethmap(df::DataFrame, args...; group::Symbol=:NAME, color::Symbol=:feature, namedargs...)`

Produces a Gadfly plot using a mapify'd DataFrame. `group` refers to the
column in `df` that identifies the shapes to be plotted (it "groups" the
polygon points), and `color` refers to the column containing the statistic
of interest for the choropleth map. Set `graticule=false` to turn off the
default graticule. `args` and `namedargs` are optional, and
are passed through to Gadfly.plot untouched. This function sets up a number of
defaults for Gadfly. If you intend to heavily style your maps using
Gadfly.Theme, you may find it easier to copy the code from this function and
modify the underlying call to Gadfly.plot directly.
"""
function choroplethmap(df::DataFrame, args...; group::Symbol=:NAME, color::Symbol=:feature, graticule::Union{Compose.Context, Bool}=true ,namedargs...)
  if isa(graticule, Bool)
    graticule = graticule == false ? compose(context()) : ChoroplethMaps.graticule(df)
  end
  Gadfly.plot(df, x=:CM_X, y=:CM_Y, group=group, color=color,
              Geom.polygon(preserve_order=true, fill=true),
              Coord.cartesian(fixed=true),
              underguide(graticule),
              Guide.xticks( label=false), Guide.yticks( label=false),
              Guide.XLabel(""), Guide.YLabel(""),
              Theme(grid_color=colorant"rgba(0,0,0,0.0)"),
              args...; namedargs...)
end

"""
`choroplethlayer(df::DataFrame, args...; fill::Bool=true, group::Symbol=:NAME, color::Symbol=:feature, namedargs...)`

Similar to `choroplethmap`, but produces a layer suitable for passing into
Gadfly.plot. `fill` controls whether the polygons in the layer are filled
or merely outlines. You likely will want to add
`Coord.cartesian(fixed=true)` to your Gadfly.plot call.
"""
function choroplethlayer(df::DataFrame, args...; fill::Bool=true, group::Symbol=:NAME, color::Symbol=:feature, namedargs...)
  if fill
    Gadfly.layer(df, x=:CM_X, y=:CM_Y, group=group, color=color,
                 Geom.polygon(preserve_order=true, fill=fill), args...; namedargs...)
  else
    # Omit the color aesthetic, since it makes no sense here and may not have
    # been passed in at all
    Gadfly.layer(df, x=:CM_X, y=:CM_Y, group=group,
                 Geom.polygon(preserve_order=true, fill=fill), args...)
  end
end

"""
`projstr(str::AbstractString)`

Returns a Proj4 projection string based on the contents of `str`, which
can have any of the following forms:

* "EPSG:<wkid>"
* "ESRI:<wkid>"
* "<wkid>"
* "+proj=..."

In each of the first three forms, `<wkid>` refers to an integer well-known id.
The third form is provided as a convenience, but be aware that it may not
return the exact projection you expect -- it's better to specify the authority
explicitly as in the first two forms, if using a wkid.
The last form represents a valid Proj4 string, resulting in the contents
of `str` simply being returned as the result.
"""
function projstr(str::AbstractString)
  # Plain Proj4 string?
  m = match(r"\+proj"i, str)
  isa(m, RegexMatch) && return str

  # EPSG wkid?
  m = match(r"epsg:(\d+)"i, str)
  isa(m, RegexMatch) && return Proj4.epsg[parse(Int64, m.captures[1])]

  # ESRI wkid?
  m = match(r"esri:(\d+)"i, str)
  isa(m, RegexMatch) && return Proj4.esri[parse(Int64, m.captures[1])]

  try
    return projstr(parse(Int64, str))
  catch
    error("Unknown crs/srs: $str")
  end
end

"""
`projstr(wkid::Integer)`

Returns a Proj4 projection string based on the contents of `wkid`. First
tries to interpret the wkid as an EPSG code, and, if that fails, tries to
interpret it as an ESRI code.
"""
function projstr(wkid::Integer)
  try
    return Proj4.epsg[wkid]
  catch
    return Proj4.esri[wkid]
  end
end

function dropallbut!(df, keep)
  for name in names(df)
    !(name in keep) && delete!(df, name)
  end
end

function add_xyp_cols!(df)
  lat = DataArray(Float64[], Bool[])
  lon = DataArray(Float64[], Bool[])
  prj = DataArray(AbstractString[], Bool[])
  insert!(df, 1, prj, :CM_P)
  insert!(df, 1, lat, :CM_Y)
  insert!(df, 1, lon, :CM_X)
end

# Political boundaries often contain multiple polygons buried inside a single
# "polygon" that thematically represent pieces of the same unit; eg. a territory
# with islands. To ensure these are rendered correctly, these have to be split
# up while still maintaining all other associated data.
function splitmultipolygons(df, key)

  indices = findindices(df)

  # "Convert" key column into a string so we can safely append an ordered tag
  # that will be used for render grouping
  newkeys = AbstractString[ string(k) for k in df[key] ]
  rename!(df, key, :CM_ORIG_KEY)
  insert!(df, 1, DataArray(newkeys), key)
  newdf = similar(df, 0)

  for (row, polyindices) in zip(eachrow(df), indices)
    polyset = row[:CM_POLYGONS]
    for (pidxs, idx) in zip(polyindices, collect(1:length(polyindices)))
      copyrow!(newdf, row, exclude=[:CM_POLYGONS, key])
      push!(newdf[key], "$(row[key])_$idx")
      start, finish = pidxs
      push!(newdf[:CM_POLYGONS], polyset[start:finish])
    end
  end
  return newdf
end

function findindices(df)
  indices = Array{Array{Int64}}[]
  for row in eachrow(df)
    polygon = row[:CM_POLYGONS]
    rowindices = Array{Int64}[]
    idx = 1
    while idx < length(polygon)
      start = idx
      opener = polygon[idx]
      idx += 1 # begin search at next index
      while polygon[idx] != opener && idx != length(polygon)
        idx += 1
      end
      push!(rowindices, Int64[start, idx])
      idx += 1
    end
    push!(indices, rowindices)
  end
  return indices
end

function copyrow!(df, row; exclude=Symbol[])
  for name in names(row)
    !(name in exclude) && push!(df[name], row[name])
  end
end

function objdiff(a, b)
  typeof(a) != typeof(b) && error("Objects must be of same type to diff")
  return filter((f)->getfield(a, f) != getfield(b, f), fieldnames(a))
end

end # module


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

module Random

using Base.dSFMT
using Base.GMP: Limb, MPZ
import Base: copymutable, copy, copy!, ==, hash

export srand,
       rand, rand!,
       randn, randn!,
       randexp, randexp!,
       bitrand,
       randstring,
       randsubseq, randsubseq!,
       shuffle, shuffle!,
       randperm, randperm!,
       randcycle, randcycle!,
       AbstractRNG, MersenneTwister, RandomDevice,
       GLOBAL_RNG, randjump

abstract type AbstractRNG end

abstract type FloatInterval{T<:AbstractFloat} end

struct CloseOpen{  T<:AbstractFloat} <: FloatInterval{T} end # interval [0,1)
struct Close1Open2{T<:AbstractFloat} <: FloatInterval{T} end # interval [1,2)

const FloatInterval_64 = FloatInterval{Float64}
const CloseOpen_64     = CloseOpen{Float64}
const Close1Open2_64   = Close1Open2{Float64}

CloseOpen(  ::Type{T}=Float64) where {T<:AbstractFloat} = CloseOpen{T}()
Close1Open2(::Type{T}=Float64) where {T<:AbstractFloat} = Close1Open2{T}()

const BitFloatType = Union{Type{Float16},Type{Float32},Type{Float64}}

function __init__()
    try
        srand()
    catch ex
        Base.showerror_nostdio(ex,
            "WARNING: Error during initialization of module Random")
    end
end

include("RNGs.jl")
include("generation.jl")
include("normal.jl")
include("misc.jl")

## rand & rand! & srand docstrings

"""
    rand([rng=GLOBAL_RNG], [S], [dims...])

Pick a random element or array of random elements from the set of values specified by `S`;
`S` can be

* an indexable collection (for example `1:n` or `['x','y','z']`),
* an `Associative` or `AbstractSet` object,
* a string (considered as a collection of characters), or
* a type: the set of values to pick from is then equivalent to `typemin(S):typemax(S)` for
  integers (this is not applicable to [`BigInt`](@ref)), and to ``[0, 1)`` for floating
  point numbers;

`S` defaults to [`Float64`](@ref).

# Examples
```julia-repl
julia> rand(Int, 2)
2-element Array{Int64,1}:
 1339893410598768192
 1575814717733606317

julia> rand(MersenneTwister(0), Dict(1=>2, 3=>4))
1=>2
```

!!! note
    The complexity of `rand(rng, s::Union{Associative,AbstractSet})`
    is linear in the length of `s`, unless an optimized method with
    constant complexity is available, which is the case for `Dict`,
    `Set` and `IntSet`. For more than a few calls, use `rand(rng,
    collect(s))` instead, or either `rand(rng, Dict(s))` or `rand(rng,
    Set(s))` as appropriate.
"""
rand

"""
    rand!([rng=GLOBAL_RNG], A, [S=eltype(A)])

Populate the array `A` with random values. If `S` is specified
(`S` can be a type or a collection, cf. [`rand`](@ref) for details),
the values are picked randomly from `S`.
This is equivalent to `copy!(A, rand(rng, S, size(A)))`
but without allocating a new array.

# Examples
```jldoctest
julia> rng = MersenneTwister(1234);

julia> rand!(rng, zeros(5))
5-element Array{Float64,1}:
 0.590845
 0.766797
 0.566237
 0.460085
 0.794026
```
"""
rand!

"""
    srand([rng=GLOBAL_RNG], seed) -> rng
    srand([rng=GLOBAL_RNG]) -> rng

Reseed the random number generator. If a `seed` is provided, the RNG will give a
reproducible sequence of numbers, otherwise Julia will get entropy from the system. For
`MersenneTwister`, the `seed` may be a non-negative integer or a vector of [`UInt32`](@ref)
integers. `RandomDevice` does not support seeding.

# Examples
```jldoctest
julia> srand(1234);

julia> x1 = rand(2)
2-element Array{Float64,1}:
 0.590845
 0.766797

julia> srand(1234);

julia> x2 = rand(2)
2-element Array{Float64,1}:
 0.590845
 0.766797

julia> x1 == x2
true
```
"""
srand

end # module

max_stars_repo_name: baggepinnen/julia

max_issues_repo_name: baggepinnen/julia

max_forks_repo_name: afniedermayer/julia


text:
@testset "context" begin
    @test isnothing(get(current_context(), :say, nothing))
    with_context(:say=>"foo") do
        @test current_context()[:say] == "foo"
        with_context(:say => "bar") do
            @test current_context()[:say] == "bar"
        end
        @test current_context()[:say] == "foo"
    end
    @test isnothing(get(current_context(), :say, nothing))
end


text:
# Autogenerated wrapper script for ConnectFourSolver_jll for aarch64-linux-musl-cxx11
export libconnect4jl

using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("ConnectFourSolver")
JLLWrappers.@declare_library_product(libconnect4jl, "libconnect4jl.so")
function __init__()
    JLLWrappers.@generate_init_header(libcxxwrap_julia_jll)
    JLLWrappers.@init_library_product(
        libconnect4jl,
        "lib/libconnect4jl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using QuantEcon
using Plots
pyplot()

# == Plot functions == #

function plot_spectral_density(arma::ARMA)
    (w, spect) = spectral_density(arma, two_pi=false)
    p = plot(w, spect, color=:blue, linewidth=2, alpha=0.7,
             xlims=(0, pi), xlabel="frequency", ylabel="spectrum",
             title="Spectral density", yscale=:log, legend=:none, grid=false)
    return p
end

function plot_autocovariance(arma::ARMA)
    acov = autocovariance(arma)
    n = length(acov)
    N = repmat(0:(n - 1), 1, 2)'
    heights = [zeros(1,n); acov']
    p = scatter(0:(n - 1), acov, title="Autocovariance", xlims=(-0.5, n - 0.5),
                xlabel="time", ylabel="autocovariance", legend=:none, color=:blue)
    plot!(-1:(n + 1), zeros(1, n + 3), color=:red, linewidth=0.5)
    plot!(N, heights, color=:blue, grid=false)
    return p
end

function plot_impulse_response(arma::ARMA)
    psi = impulse_response(arma)
    n = length(psi)
    N = repmat(0:(n - 1), 1, 2)'
    heights = [zeros(1,n); psi']
    p = scatter(0:(n - 1), psi, title="Impulse response", xlims=(-0.5, n - 0.5),
                xlabel="time", ylabel="response", legend=:none, color=:blue)
    plot!(-1:(n + 1), zeros(1, n + 3), color=:red, linewidth=0.5)
    plot!(N, heights, color=:blue, grid=false)
    return p
end

function plot_simulation(arma::ARMA)
    X = simulation(arma)
    n = length(X)
    p = plot(0:(n - 1), X, color=:blue, linewidth=2, alpha=0.7,
             xlims=(0.0, n), xlabel="time", ylabel="state space",
            title="Sample path", legend=:none, grid=:false)
    return p
end

function quad_plot(arma::ARMA)
    p = plot(plot_impulse_response(arma), plot_autocovariance(arma),
         plot_spectral_density(arma), plot_simulation(arma))
    return p
end

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
abstract type TreeOptions end

## Box tree
# Implements quad- (2D) and oct-tree (3D) cluster strategies
include("boundingbox.jl")
include("boxtree.jl")

max_stars_repo_name: JoshuaTetzner/FastBEAST

max_issues_repo_name: JoshuaTetzner/FastBEAST

max_forks_repo_name: JoshuaTetzner/FastBEAST


text:
@testset "argmin_x" begin

#test argmin_x.jl for CSC sparse matrix format
  A = sprandn(100,100,0.01)+sparse(I, 100, 100);
  A = A'*A;
  while rank(Array(A))<100
    A = A + sparse(I, 100, 100)
  end
  xt = randn(100)
  b  = A*xt;

  maxit = 100;
  pp    = 2
  p     = 2
  TF    = Float64
  x_solve_tol_ref = 1e-5

  log_PARSDMM = log_type_PARSDMM(zeros(maxit,pp),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),zeros(maxit),
  zeros(maxit),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),TimerOutput());

  #test for zero initial guess
    (x,iter,relres)=argmin_x(A,b,zeros(100),x_solve_tol_ref,5,log_PARSDMM)
    @test norm(A*x-b)/norm(b) <= x_solve_tol_ref
    @test isapprox(norm(A*x-b)/norm(b),relres,rtol=x_solve_tol_ref)

  #test for good initial guess
    (x2,iter2,relres2)=argmin_x(A,b,xt+randn(100).*1e-4,1e-5,5,log_PARSDMM)
    @test iter2<iter
    @test norm(A*x2-b)/norm(b) <= 1e-5
    @test isapprox(norm(A*x2-b)/norm(b),relres2,rtol=1e-5)

  #test if 10*eps(Float32/Float64) is achieved
    (x,iter,relres)=argmin_x(A,b,zeros(100),10*eps(),15,log_PARSDMM)
    @test norm(A*x-b)/norm(b) <= 20*eps()

#test argmin_x.jl for CDS sparse matrix format
  A=sprandn(100,100,0.01)+sparse(I, 100, 100);
  A=A'*A;
  while rank(Array(A))<100
    A=A+sparse(I, 100, 100)
  end
  xt=randn(100)
  b=A*xt;

  (R,offset) = mat2CDS(A)

  y=zeros(Float64,length(b))
  (x,iter,relres)=argmin_x(A,b,zeros(100),1e-5,5,log_PARSDMM,offset,y)

  @test norm(A*x-b)/norm(b) <= 1e-5
  @test norm(A*x-b)/norm(b) <= 2.0*relres

#test for more accurate solutions with CDS
  y=zeros(Float64,length(b))
  (x,iter,relres)=argmin_x(A,b,zeros(100),1e-10,5,log_PARSDMM,offset,y)

  @test norm(A*x-b)/norm(b) <= 1e-10
  @test norm(A*x-b)/norm(b) <= 2.0*relres

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
# BFS demo using cvdlab SuiteSparseGraphBLAS Julia wrapper
# (can be found at https://github.com/cvdlab/SuiteSparseGraphBLAS.jl)

#=
    Since Julia does not allow to have two direct dependencies with the same name
    in a single project, in order to execute the following code, you have to go inside
    the directory "sparse/examples/GraphBLAS_Julia_cvdlab/"; the corresponding environment
    will be activated by the code below. Additionally, the necessary src files are
    manually added below. Before executing the other demos, head back to the main
    project environment. Still, it is suggested to end the ongoing Julia session and
    to execute the other demos in another Julia session.
    Suggested flow of execution:

    $ cd /sparse/examples/GraphBLAS_Julia_cvdlab
    $ julia

    julia>[paste the code below...]

    [end Julia session...]

    $ cd ../../
    $ julia

    julia> [paste the code from another demo...]
=#

# Activate "sparse/examples/GraphBLAS_Julia_cvdlab/" environment
using Pkg
Pkg.activate(".")

using SuiteSparseGraphBLAS, BenchmarkTools

# Manually add the necessary src files
include("../../src/GraphBLAS_sparse_square_matrix_cvd.jl")
include("../../src/GraphBLAS_bfs_cvd.jl")

nrows = 100
ncols = 100

A = from_type(Int64, nrows, ncols)

GraphBLAS_sparse_square_matrix_cvd!(A)

bench = @benchmark GraphBLAS_bfs_cvd(A, 1)

display(bench)


text:
using Optim, Random
obj = build_loss_objective(prob1,Tsit5(),L2Loss(t,data),
                           maxiters=10000,verbose=false)

### Optim Method

println("Use Optim Brent to fit the parameter")
result = Optim.optimize(obj, 1.0, 10.0)
@test_broken result.minimizer[1] ≈ 1.5 atol=3e-1

println("Use Optim BFGS to fit the parameter")
result = Optim.optimize(obj, [1.0], Optim.BFGS())
@test result.minimizer[1] ≈ 1.5 atol=3e-1
#sol_optimized2 = solve(prob)
#plot!(sol_optimized2,leg=false)

cost_function2 = build_loss_objective(prob2,Tsit5(),L2Loss(t,data),
                                      maxiters=10000,verbose=false)
result_bfgs = Optim.optimize(cost_function2, [1.0,2.5], Optim.BFGS())
@test_broken result_bfgs.minimizer ≈ [1.5;3.0] atol=3e-1

Random.seed!(200)
cost_function3 = build_loss_objective(prob3,Tsit5(),L2Loss(t,data),
                                      maxiters=10000,verbose=false)
result_bfgs = Optim.optimize(cost_function3, [1.3,0.8,2.8,1.2], Optim.BFGS())
@test result_bfgs.minimizer ≈ [1.5;1.0;3.0;1.0] atol=5e-1

#=
println("Use Optim NelderMead to fit the parameter")
result_neldermead = Optim.optimize(cost_function2, [1.0,2.5], Optim.NelderMead())
@test result_neldermead.minimizer ≈ [1.5;3.0] atol=3e-1

println("Use Optim SimulatedAnnealing to fit the parameter")
options = Optim.Options(iterations=100_000)
result_simulatedannealing = Optim.optimize(obj, [1.0], Optim.SimulatedAnnealing(),options)
@test result_simulatedannealing.minimizer ≈ [1.5] atol=3e-1

result_simulatedannealing = Optim.optimize(cost_function2, [1.0,2.5], Optim.SimulatedAnnealing(),options)
@test result_simulatedannealing.minimizer ≈ [1.5;3.0] atol=3e-1

println("Use Optim Particle Swarm Optimization to fit the parameter")
initial_x = [1.0]
lower = [0.0]
upper = [5.0]
n_particles = 4
options = Optim.Options(iterations=100)
result_pso = Optim.optimize(obj, initial_x, Optim.ParticleSwarm(lower, upper, n_particles), options)
@test result_pso.minimizer[1] ≈ 1.5 atol=3e-1

initial_x = [1.0,2.5]
lower = [0.0,2.0]
upper = [5.0,10.0]
n_particles = 4
options = Optim.Options(iterations=100)
result_pso = Optim.optimize(cost_function2, initial_x, Optim.ParticleSwarm(lower, upper, n_particles), options)
@test result_pso.minimizer ≈ [1.5;3.0] atol=3e-1
=#


text:
# Local Tangent Space Alignment (LTSA)
# ---------------------------
# Principal Manifolds and Nonlinear Dimension Reduction via Local Tangent Space Alignment,
# Zhang, Zhenyue; Hongyuan Zha (2004),  SIAM Journal on Scientific Computing 26 (1): 313–338.
# doi:10.1137/s1064827502419154.

#### LTSA type
struct LTSA{T <: AbstractFloat} <: SpectralResult
    k::Int
    λ::AbstractVector{T}
    proj::Projection{T}

    LTSA{T}(k::Int, λ::AbstractVector{T}, proj::Projection{T}) where T = new(k, λ, proj)
end

## properties
outdim(M::LTSA) = size(M.proj, 1)
projection(M::LTSA) = M.proj

eigvals(M::LTSA) = M.λ
neighbors(M::LTSA) = M.k

## show & dump
function show(io::IO, M::LTSA)
    print(io, "LTSA(outdim = $(outdim(M)), neighbors = $(neighbors(M)))")
end

function dump(io::IO, M::LTSA)
    show(io, M)
    println(io, "eigenvalues: ")
    Base.showarray(io, transpose(M.λ), header=false, repr=false)
    println(io)
    println(io, "projection:")
    Base.showarray(io, M.proj, header=false, repr=false)
end

## interface functions
function transform(::Type{LTSA}, X::DenseMatrix{T}; d::Int=2, k::Int=12) where T<:AbstractFloat
    n = size(X, 2)

    # Construct NN graph
    D, I = find_nn(X, k)

    B = spzeros(n,n)
    for i=1:n
        # re-center points in neighborhood
        μ = mean(X[:,I[:,i]], dims=2)
        δ_x = X[:,I[:,i]] .- μ

        # Compute orthogonal basis H of θ'
        θ_t = svd(δ_x).V[:,1:d]

        # Construct alignment matrix
        G = hcat(ones(k)./sqrt(k), θ_t)
        B[I[:,i], I[:,i]] =  B[I[:,i], I[:,i]] + Matrix{Float64}(LinearAlgebra.I, k, k) - G*transpose(G)
    end

    # Align global coordinates
    λ, V = decompose(B, d)
    return LTSA{T}(k, λ, transpose(V))
end