text:
# Copyright 2020 Adam Liddell
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test
using CSV
using Printf
using DataStructures

using S2Cell

@info "Low-level internals"
@test_throws ArgumentError("Cannot convert UV to XYZ with invalid face: 6") S2Cell._s2_face_uv_to_xyz(6, (0, 0))

encode_tbl = CSV.File("./s2_encode_corpus.csv.gz"; types=[Float64, Float64, Int32, UInt64, String])
decode_tbl = CSV.File("./s2_decode_corpus.csv.gz"; types=[UInt64, String, Float64, Float64, Int32])

@info "Cell to token"
@test S2Cell.cell_id_to_token(S2Cell.Cell(0)) == "X"

for row in encode_tbl
    try
        @test S2Cell.cell_id_to_token(S2Cell.Cell(UInt64(row[:cell_id]))) == row[:token]
    catch e
        println(row)
        break
    end
end

@info "Token to cell"
@test S2Cell.token_to_cell_id("x").id == 0
@test S2Cell.token_to_cell_id("X").id == 0

@test_throws S2Cell.InvalidToken("Cannot convert S2 token with length > 16 characters") S2Cell.token_to_cell_id('a' ^ 17)

for row in encode_tbl
    @test S2Cell.token_to_cell_id(row[:token]).id == row[:cell_id]
end

@info "LL to *"
@test_throws ArgumentError S2Cell.lat_lon_to_cell_id(0., 0., -1)
@test_throws ArgumentError S2Cell.lat_lon_to_cell_id(0., 0., 31)

for row in encode_tbl
    @test S2Cell.lat_lon_to_cell_id(row[:lat], row[:lon], row[:level]).id == row[:cell_id]
end

@test_throws ArgumentError S2Cell.lat_lon_to_token(0, 0, -1)
@test_throws ArgumentError S2Cell.lat_lon_to_token(0, 0, 31)

for row in encode_tbl
    @test S2Cell.lat_lon_to_token(row[:lat], row[:lon], row[:level]) == row[:token]
end

@info "* to LL"
@test_throws S2Cell.InvalidCellID S2Cell.cell_id_to_lat_lon((UInt64(0b110) << S2Cell._S2_POS_BITS))

for row in decode_tbl
    @test all(S2Cell.cell_id_to_lat_lon(row[:cell_id]) .≈ (row[:lat], row[:lon]))
end

@test_throws S2Cell.InvalidToken S2Cell.token_to_lat_lon('a' ^ 17)

@test_throws S2Cell.InvalidToken S2Cell.token_to_lat_lon(@sprintf "%016x" (UInt64(0b110) << S2Cell._S2_POS_BITS))

for row in decode_tbl
    @test all(S2Cell.token_to_lat_lon(row[:token]) .≈ (row[:lat], row[:lon]))
end

@info "Token miscellany"
for (token, expected) in [
    ("3", "3"), ("2ef59bd352b93ac3", "2ef59bd352b93ac3"), ("  2ef59bd352b93ac3", "2ef59bd352b93ac3"),
    ("2ef59bd352b93ac3  ", "2ef59bd352b93ac3"), (" 2ef59bd352b93ac3 ", "2ef59bd352b93ac3"),
    ("2EF", "2ef"), ("2eF", "2ef"), ("2ef000", "2ef"), ("", "X"), ("x", "X"),
    ]
    @test S2Cell.token_to_canonical_token(token) == expected
end

evens = [ (UInt64(1) << even_number, true) for even_number in 0:2:S2Cell._S2_POS_BITS]
odds = [ (UInt64(1) << even_number, false) for even_number in 1:2:S2Cell._S2_POS_BITS]
trivial = [
    (UInt64(0b1111010101010101010101010101010101010101010101010101010101010101), false),  # Invalid face
    (UInt64(0), false),
]

for (id, is_valid) in cat(trivial, evens, odds, dims=1)
    @test S2Cell.cell_id_is_valid(id) == is_valid
end

for row in decode_tbl
    S2Cell.cell_id_is_valid(row[:cell_id])
end

for (token, is_valid) in [
    ("", false),  # Invalid cell ID
    ("x", false),  # Invalid cell ID
    ("X", false),  # Invalid cell ID
    ("2ef", true),
    ("2EF", true),
    ("2ef0000000000000", true),
    ("2ef00000000000000", false),  # Too long
    ("2efinvalid", false),  # Invalid characters
    ("2efx", false),  # Incorrect use of X
    ]
    @test S2Cell.token_is_valid(token) == is_valid
end

for row in decode_tbl
    @test S2Cell.token_is_valid(row[:token])
end

@info "Level shifting"
@test_throws S2Cell.InvalidCellID S2Cell.cell_id_to_level(0)

for row in encode_tbl
    @test S2Cell.cell_id_to_level(row[:cell_id]) == row[:level]
end

@test_throws S2Cell.InvalidToken S2Cell.token_to_level('a' ^ 17)
@test_throws S2Cell.InvalidToken S2Cell.token_to_level("")

for row in encode_tbl
    @test S2Cell.token_to_level(row[:token]) == row[:level]
end

@test_throws S2Cell.InvalidCellID S2Cell.cell_id_to_parent_cell_id(0)
@test_throws ArgumentError S2Cell.cell_id_to_parent_cell_id(3458764513820540928)
@test_throws ArgumentError S2Cell.cell_id_to_parent_cell_id(3383782026652942336, -1)
@test_throws ArgumentError S2Cell.cell_id_to_parent_cell_id(3383782026652942336, 31)
@test_throws ArgumentError S2Cell.cell_id_to_parent_cell_id(3383782026652942336, 16)

points = DefaultDict(Dict)
for row in encode_tbl
    points[(row[:lat], row[:lon])][row[:level]] = row[:cell_id]
end
for levels_dict in rand(collect(values(points)), 100)
    for level in keys(levels_dict)
        if level > 0
            @test S2Cell.cell_id_to_parent_cell_id(levels_dict[level]) == levels_dict[level-1]
        else
            @test_throws ArgumentError S2Cell.cell_id_to_parent_cell_id(levels_dict[level])
        end

        for other in keys(levels_dict)
            if other <= level
                @test S2Cell.cell_id_to_parent_cell_id(levels_dict[level], other) == levels_dict[other]
            else
                @test_throws ArgumentError S2Cell.cell_id_to_parent_cell_id(levels_dict[level], other)
            end
        end
    end
end

@test_throws S2Cell.InvalidToken S2Cell.token_to_parent_token('a' ^ 17)
@test_throws ArgumentError S2Cell.token_to_parent_token("3")
@test_throws ArgumentError S2Cell.token_to_parent_token("2ef59bd34", -1)
@test_throws ArgumentError S2Cell.token_to_parent_token("2ef59bd34", 31)
@test_throws ArgumentError S2Cell.token_to_parent_token("2ef59bd34", 16)

points = DefaultDict(Dict)
for row in encode_tbl
    points[(row[:lat], row[:lon])][row[:level]] = row[:token]
end
for levels_dict in rand(collect(values(points)), 100)
    for level in keys(levels_dict)
        if level > 0
            @test S2Cell.token_to_parent_token(levels_dict[level]) == levels_dict[level-1]
        else
            @test_throws ArgumentError S2Cell.token_to_parent_token(levels_dict[level])
        end

        for other in keys(levels_dict)
            if other <= level
                @test S2Cell.token_to_parent_token(levels_dict[level], other) == levels_dict[other]
            else
                @test_throws ArgumentError S2Cell.token_to_parent_token(levels_dict[level], other)
            end
        end
    end
end


text:
"""
    mutable struct StandardModel

`StandardModel` is used to store a constraint based metabolic model with
meta-information.  Meta-information is defined as annotation details, which
include gene-reaction-rules, formulas, etc.

This model type seeks to keep as much meta-information as possible, as opposed
to `CoreModel` and `CoreModelCoupled`, which keep the bare neccessities only.
When merging models and keeping meta-information is important, use this as the
model type.  If meta-information is not important, use the more efficient core
model types.  See [`CoreModel`](@ref) and [`CoreModelCoupled`](@ref) for
comparison.

In this model, reactions, metabolites, and genes are stored in ordered
dictionaries indexed by each struct's `id` field.  For example,
`model.reactions["rxn1_id"]` returns a `Reaction` where the field `id` equals
`"rxn1_id"`.  This makes adding and removing reactions efficient.

Note that the stoichiometric matrix (or any other core data, e.g. flux bounds)
is not stored directly as in `CoreModel`.  When this model type is used in
analysis functions, these core data structures are built from scratch each time
an analysis function is called.  This can cause performance issues if you run
many small analysis functions sequentially.  Consider using the core model
types if performance is critical.

See also: [`Reaction`](@ref), [`Metabolite`](@ref), [`Gene`](@ref)

# Fields
```
id :: String
reactions :: OrderedDict{String, Reaction}
metabolites :: OrderedDict{String, Metabolite}
genes :: OrderedDict{String, Gene}
```

# Example
```
model = load_model(StandardModel, "my_model.json")
keys(model.reactions)
```
"""
mutable struct StandardModel <: MetabolicModel
    id::String
    reactions::OrderedDict{String,Reaction}
    metabolites::OrderedDict{String,Metabolite}
    genes::OrderedDict{String,Gene}

    StandardModel(
        id = "";
        reactions = OrderedDict{String,Reaction}(),
        metabolites = OrderedDict{String,Metabolite}(),
        genes = OrderedDict{String,Gene}(),
    ) = new(id, reactions, metabolites, genes)
end

# MetabolicModel interface follows
"""
    reactions(model::StandardModel)

Return a vector of reaction id strings contained in `model`.
The order of reaction ids returned here matches the order used to construct the
stoichiometric matrix.
"""
reactions(model::StandardModel)::StringVecType = collect(keys(model.reactions))

"""
    n_reactions(model::StandardModel)

Return the number of reactions contained in `model`.
"""
n_reactions(model::StandardModel)::Int = length(model.reactions)

"""
    metabolites(model::StandardModel)

Return a vector of metabolite id strings contained in `model`.
The order of metabolite strings returned here matches the order used to construct
the stoichiometric matrix.
"""
metabolites(model::StandardModel)::StringVecType = collect(keys(model.metabolites))

"""
n_metabolites(model::StandardModel)

Return the number of metabolites in `model`.
"""
n_metabolites(model::StandardModel)::Int = length(model.metabolites)

"""
    genes(model::StandardModel)

Return a vector of gene id strings in `model`.
"""
genes(model::StandardModel)::StringVecType = collect(keys(model.genes))

"""
    n_genes(model::StandardModel)

Return the number of genes in `model`.
"""
n_genes(model::StandardModel)::Int = length(model.genes)

"""
    stoichiometry(model::StandardModel)

Return the stoichiometric matrix associated with `model` in sparse format.
"""
function stoichiometry(model::StandardModel)::SparseMat
    n_entries = 0
    for (_, r) in model.reactions
        for _ in r.metabolites
            n_entries += 1
        end
    end

    MI = Vector{Int}()
    RI = Vector{Int}()
    SV = Vector{Float64}()
    sizehint!(MI, n_entries)
    sizehint!(RI, n_entries)
    sizehint!(SV, n_entries)

    # establish the ordering
    rxns = reactions(model)
    met_idx = Dict(mid => i for (i, mid) in enumerate(metabolites(model)))

    # fill the matrix entries
    for (ridx, rid) in enumerate(rxns)
        for (mid, coeff) in model.reactions[rid].metabolites
            haskey(met_idx, mid) || throw(
                DomainError(
                    mid,
                    "Metabolite $(mid) not found in model but occurs in stoichiometry of $(rid)",
                ),
            )
            push!(MI, met_idx[mid])
            push!(RI, ridx)
            push!(SV, coeff)
        end
    end
    return SparseArrays.sparse(MI, RI, SV, n_metabolites(model), n_reactions(model))
end

"""
    lower_bounds(model::StandardModel)::Vector{Float64}

Return the lower bounds for all reactions in `model` in sparse format.
"""
lower_bounds(model::StandardModel)::Vector{Float64} =
    sparse([model.reactions[rxn].lb for rxn in reactions(model)])

"""
    upper_bounds(model::StandardModel)::Vector{Float64}

Return the upper bounds for all reactions in `model` in sparse format.
Order matches that of the reaction ids returned in `reactions()`.
"""
upper_bounds(model::StandardModel)::Vector{Float64} =
    sparse([model.reactions[rxn].ub for rxn in reactions(model)])

"""
    bounds(model::StandardModel)::Tuple{Vector{Float64},Vector{Float64}}

Return the lower and upper bounds, respectively, for reactions in `model`.
Order matches that of the reaction ids returned in `reactions()`.
"""
bounds(model::StandardModel)::Tuple{Vector{Float64},Vector{Float64}} =
    (lower_bounds(model), upper_bounds(model))

"""
    balance(model::StandardModel)

Return the balance of the linear problem, i.e. b in Sv = 0 where S is the stoichiometric matrix
and v is the flux vector.
"""
balance(model::StandardModel)::SparseVec = spzeros(length(model.metabolites))

"""
    objective(model::StandardModel)

Return sparse objective vector for `model`.
"""
objective(model::StandardModel)::SparseVec =
    sparse([model.reactions[rid].objective_coefficient for rid in keys(model.reactions)])

"""
    reaction_gene_association(model::StandardModel, id::String)

Return the gene reaction rule in string format for reaction with `id` in `model`.
Return `nothing` if not available.
"""
reaction_gene_association(model::StandardModel, id::String)::Maybe{GeneAssociation} =
    _maybemap(identity, model.reactions[id].grr)

"""
    metabolite_formula(model::StandardModel, id::String)

Return the formula of reaction `id` in `model`.
Return `nothing` if not present.
"""
metabolite_formula(model::StandardModel, id::String)::Maybe{MetaboliteFormula} =
    _maybemap(_parse_formula, model.metabolites[id].formula)

"""
    metabolite_charge(model::StandardModel, id::String)

Return the charge associated with metabolite `id` in `model`.
Return nothing if not present.
"""
metabolite_charge(model::StandardModel, id::String)::Maybe{Int} =
    model.metabolites[id].charge

"""
    metabolite_compartment(model::StandardModel, id::String)

Return compartment associated with metabolite `id` in `model`.
Return `nothing` if not present.
"""
metabolite_compartment(model::StandardModel, id::String)::Maybe{String} =
    model.metabolites[id].compartment

"""
    reaction_subsystem(id::String, model::StandardModel)

Return the subsystem associated with reaction `id` in `model`.
Return `nothing` if not present.
"""
reaction_subsystem(model::StandardModel, id::String)::Maybe{String} =
    model.reactions[id].subsystem

"""
    metabolite_notes(model::StandardModel, id::String)::Notes

Return the notes associated with metabolite `id` in `model`.
Return an empty Dict if not present.
"""
metabolite_notes(model::StandardModel, id::String)::Maybe{Notes} =
    model.metabolites[id].notes

"""
    metabolite_annotations(model::StandardModel, id::String)::Annotations

Return the annotation associated with metabolite `id` in `model`.
Return an empty Dict if not present.
"""
metabolite_annotations(model::StandardModel, id::String)::Maybe{Annotations} =
    model.metabolites[id].annotations

"""
    gene_notes(model::StandardModel, id::String)::Notes

Return the notes associated with gene `id` in `model`.
Return an empty Dict if not present.
"""
gene_notes(model::StandardModel, id::String)::Maybe{Notes} = model.genes[id].notes

"""
    gene_annotations(model::StandardModel, id::String)::Annotations

Return the annotation associated with gene `id` in `model`.
Return an empty Dict if not present.
"""
gene_annotations(model::StandardModel, id::String)::Maybe{Annotations} =
    model.genes[id].annotations

"""
    reaction_notes(model::StandardModel, id::String)::Notes

Return the notes associated with reaction `id` in `model`.
Return an empty Dict if not present.
"""
reaction_notes(model::StandardModel, id::String)::Maybe{Notes} = model.reactions[id].notes

"""
    reaction_annotations(model::StandardModel, id::String)::Annotations

Return the annotation associated with reaction `id` in `model`.
Return an empty Dict if not present.
"""
reaction_annotations(model::StandardModel, id::String)::Maybe{Annotations} =
    model.reactions[id].annotations

"""
    reaction_stoichiometry(model::StandardModel, rid::String)::Dict{String, Float64}

Return the stoichiometry of reaction with ID `rid`.
"""
reaction_stoichiometry(m::StandardModel, rid::String)::Dict{String,Float64} =
    m.reactions[rid].metabolites

"""
Base.convert(::Type{StandardModel}, model::MetabolicModel)

Convert any `MetabolicModel` into a `StandardModel`.
Note, some data loss may occur since only the generic interface is used during
the conversion process.
"""
function Base.convert(::Type{StandardModel}, model::MetabolicModel)
    if typeof(model) == StandardModel
        return model
    end

    id = "" # TODO: add accessor to get model ID
    modelreactions = OrderedDict{String,Reaction}()
    modelmetabolites = OrderedDict{String,Metabolite}()
    modelgenes = OrderedDict{String,Gene}()

    gids = genes(model)
    metids = metabolites(model)
    rxnids = reactions(model)

    for gid in gids
        modelgenes[gid] = Gene(
            gid;
            notes = gene_notes(model, gid),
            annotations = gene_annotations(model, gid),
        )
    end

    for mid in metids
        modelmetabolites[mid] = Metabolite(
            mid;
            charge = metabolite_charge(model, mid),
            formula = _maybemap(_unparse_formula, metabolite_formula(model, mid)),
            compartment = metabolite_compartment(model, mid),
            notes = metabolite_notes(model, mid),
            annotations = metabolite_annotations(model, mid),
        )
    end

    S = stoichiometry(model)
    lbs, ubs = bounds(model)
    ocs = objective(model)
    for (i, rid) in enumerate(rxnids)
        rmets = Dict{String,Float64}()
        for (j, stoich) in zip(findnz(S[:, i])...)
            rmets[metids[j]] = stoich
        end
        modelreactions[rid] = Reaction(
            rid;
            metabolites = rmets,
            lb = lbs[i],
            ub = ubs[i],
            grr = reaction_gene_association(model, rid),
            objective_coefficient = ocs[i],
            notes = reaction_notes(model, rid),
            annotations = reaction_annotations(model, rid),
            subsystem = reaction_subsystem(model, rid),
        )
    end

    return StandardModel(
        id;
        reactions = modelreactions,
        metabolites = modelmetabolites,
        genes = modelgenes,
    )
end

ext: jl


text:
# function OLDminDistGauss(bd::BallTreeDensity, dRoot::Int64, atTree::BallTreeDensity, aRoot::Int64)
#   #@show "minDistGauss",  dRoot, aRoot
#   atCenter = center(atTree, aRoot)
#   densCenter = center(bd, dRoot)
#   bw = bwMax(bd, dRoot)
#   result = 0.0
#   tmp = 0.0
#   for (k=1:Ndim(atTree))
#     tmp = abs( atCenter[k] - densCenter[k] )
#     tmp -= (rangeB(atTree, aRoot))[k] + (rangeB(bd, dRoot))[k];
#     tmp = (tmp > 0) ? tmp : 0.0
#     if ( bwUniform(bd) )
#       result -= (tmp*tmp)/bw[k]
#     else
#       result -= (tmp*tmp)/bw[k] + log(bwMin(bd, dRoot)[k]) end
#   end
#   result = exp(result/2.0)
#   return result
# end
#
# function OLDmaxDistGauss(bd::BallTreeDensity, dRoot::Int64, atTree::BallTreeDensity, aRoot::Int64)
#   atCenter = center(atTree, aRoot)
#   densCenter = center(bd, dRoot)
#   bw = bwMin(bd, dRoot)
#   tmp = 0.0
#   result = 0.0
#   for k in 1:Ndim(atTree.bt)
#     tmp = abs( atCenter[k] - densCenter[k] )
#     tmp+= rangeB(atTree, aRoot)[k] + rangeB(bd, dRoot)[k]
#     if ( bwUniform(bd) )
#         result -= (tmp*tmp)/bw[k]
#     else
#         result -= (tmp*tmp)/bw[k] + log(bwMax(bd, Root)[k])
#     end
#   end
#   result = exp(result/2.0);
#   return result
# end

function minDistGauss!(restmp::Array{Float64, 1}, bd::BallTreeDensity, dRoot::Int64, atTree::BallTreeDensity, aRoot::Int64)
  #@show "minDistGauss",  dRoot, aRoot
  #atCenter = center(atTree, aRoot)
  #densCenter = center(bd, dRoot)
  #bw = bwMax(bd, dRoot)
  @fastmath @inbounds begin
    restmp[1] = 0.0
    #tmp = 0.0
    for (k=1:Ndim(atTree.bt))
      restmp[2] = abs( center(atTree.bt, aRoot, k) - center(bd.bt, dRoot, k) )
      restmp[2] -= rangeB(atTree.bt, aRoot, k) + rangeB(bd.bt, dRoot, k)
      restmp[2] = (restmp[2] > 0) ? restmp[2] : 0.0
      if ( bwUniform(bd) )
        restmp[1] -= (restmp[2]*restmp[2])/bwMax(bd, dRoot, k)
      else
        restmp[1] -= (restmp[2]*restmp[2])/bwMax(bd, dRoot, k) + log(bwMin(bd, dRoot, k))
      end
    end
    restmp[1] = exp(restmp[1]/2.0)
  end
  nothing
end

function maxDistGauss!(rettmp::Array{Float64, 1}, bd::BallTreeDensity, dRoot::Int64, atTree::BallTreeDensity, aRoot::Int64)
  #atCenter = center(atTree, aRoot)
  #densCenter = center(bd, dRoot)
  #bw = bwMin(bd, dRoot)

    #tmp = 0.0
    rettmp[1] = 0.0
    for k in 1:Ndim(atTree.bt)
      @fastmath begin
        #@inbounds tmp = abs( center(atTree.bt, aRoot)[k] - center(bd.bt, dRoot)[k] )
        #@inbounds tmp += rangeB(atTree.bt, aRoot)[k] + rangeB(bd.bt, dRoot)[k]
        #@inbounds tmp = abs( atTree.bt.centers[((aRoot-1)*atTree.bt.dims+k)] - bd.bt.centers[((dRoot-1)*bd.bt.dims+k)] )
        #@inbounds tmp += atTree.bt.ranges[((aRoot-1)*atTree.bt.dims+k)] + bd.bt.ranges[((dRoot-1)*bd.bt.dims+k)]
        @inbounds rettmp[2] = abs( center(atTree.bt, aRoot, k) - center(bd.bt, dRoot, k) )
        @inbounds rettmp[2] += rangeB(atTree.bt, aRoot,k) + rangeB(bd.bt, dRoot, k)
        if ( bwUniform(bd) )
            #@inbounds result = result - (tmp*tmp)/(bwMin(bd, dRoot)[k])
            @inbounds rettmp[1] -= (rettmp[2]*rettmp[2])/bwMin(bd, dRoot, k)
        else
            #@inbounds result = result - (tmp*tmp)/(bwMin(bd, dRoot)[k]) + (log(bwMax(bd.bt, Root))[k])
            @inbounds rettmp[1] -= (rettmp[2]*rettmp[2])/(bwMin(bd, dRoot, k)) + (log(bwMax(bd.bt, Root, k)))
        end
      end
    end
    rettmp[1] = exp(rettmp[1]/2.0)
    nothing
end

#   Bounds on kernel values between points in this subtree & another
function maxDistKer!(rettmp, bd::BallTreeDensity, dRoot::Int64, atTree::BallTreeDensity, aRoot::Int64)
#switch(getType(atTree))
#  { case Gaussian:
   maxDistGauss!(rettmp, bd, dRoot, atTree, aRoot)
#    case Laplacian:      return maxDistLaplace(dRoot,atTree,aRoot);
#    case Epanetchnikov:  return maxDistEpanetch(dRoot,atTree,aRoot);
#  }
  nothing
end

function minDistKer!(rettmp, bd::BallTreeDensity, dRoot::Int64, atTree::BallTreeDensity, aRoot::Int64)
#switch(getType())
#  { case Gaussian:
   minDistGauss!(rettmp, bd, dRoot, atTree, aRoot)
#    case Laplacian:      return minDistLaplace(dRoot,atTree,aRoot);
#    case Epanetchnikov:  return minDistEpanetch(dRoot,atTree,aRoot);
#  }
  nothing
end

# need to declare these here, for kernel
type pArrHdls
    pMin::Array{Float64,1}
    pMax::Array{Float64,1}
    pAdd::Array{Float64,2}
    pErr::Array{Float64,1}
    mini::Array{Float64,1}
    maxi::Array{Float64,1}
end

function pushDownLocal(atTree::BallTreeDensity, aRoot::Int64, hdl::pArrHdls)
    if !(isLeaf(atTree, aRoot))
      close = atTree.left(aRoot);
      if (close != NO_CHILD)
        hdl.pAdd[1,close] += hdl.pAdd[1,aRoot]
      end
      close = right(atTree, aRoot);
      if (close != NO_CHILD)
        hdl.pAdd[1,close] += hdl.pAdd[1,aRoot]
      end
      hdl.pAdd[1,aRoot] = 0
    end
end

function pushDownAll(locations::BallTreeDensity, hdl::pArrHdls)
  for j in root():(leafFirst(locations,root())-1)
      hdl.pAdd[1,  left(locations, j)] += hdl.pAdd[1,j]
      hdl.pAdd[1, right(locations, j)] += hdl.pAdd[1,j]
      hdl.pAdd[1,j] = 0
    end
    for j in leafFirst(locations, root()):(leafLast(locations, root())+1)
      hdl.pMin[j] += hdl.pAdd[1,j] - hdl.pErr[j]
      hdl.pMax[j] += hdl.pAdd[1,j] + hdl.pErr[j]
      hdl.pAdd[1,j] = 0; hdl.pErr[j] = 0;
    end
end

function recurseMinMax(atTree::BallTreeDensity, aRoot::Int64, hdl::pArrHdls)
  l = left(atTree, aRoot); r = right(atTree, aRoot);
  if !(isLeaf(atTree, l)) recurseMinMax(atTree, l, hdl) end
  if !(isLeaf(atTree, r)) recurseMinMax(atTree, r, hdl) end
  hdl.pMin[aRoot] = hdl.pMin[l]; hdl.pMax[aRoot] = hdl.pMax[l];
  if (hdl.pMin[aRoot] > hdl.pMin[r]) hdl.pMin[aRoot] = hdl.pMin[r] end
  if (hdl.pMax[aRoot] < hdl.pMax[r]) hdl.pMax[aRoot] = hdl.pMax[r] end
end

function evalDirect(bd::BallTreeDensity, dRoot::Int64, atTree::BallTreeDensity, aRoot::Int64, hdl::pArrHdls)
  firstFlag = true;
  minVal=2e22;
  maxVal=0.0;
  restmp = Array{Float64,1}(2)
  d = 0.0
  for j in leafFirst(atTree.bt, aRoot):leafLast(atTree.bt, aRoot)
    for i in leafFirst(bd.bt,dRoot):leafLast(bd.bt, dRoot)
      if (bd != atTree || i != j)                              # Check leave-one-out condition;
        #d = weight(bd.bt, i) * maxDistKer(bd, i, atTree, j)   #  Do direct N^2 kernel evaluation
        maxDistKer!(restmp, bd, i, atTree, j)
        d = bd.bt.weights[i] * restmp[1]                       #  Do direct N^2 kernel evaluation
        hdl.pMin[j] = hdl.pMin[j] + d
        hdl.pMax[j] = hdl.pMax[j] + d
      end
    end
    #D = 0.0
    #D = @parallel (+) for i=leafFirst(bd,dRoot):leafLast(bd, dRoot)
    #  (bd != atTree || i != j) ? weight(bd, i) * maxDistKer(bd, i, atTree, j) : 0.0               # Check leave-one-out condition;
    #end
    #hdl.pMin[j] += D
    #hdl.pMax[j] += D

    @inbounds if (hdl.pMin[j] < minVal) minVal = hdl.pMin[j]; end  # determine min & max value in this block
    @inbounds if (hdl.pMax[j] > maxVal) maxVal = hdl.pMax[j]; end
  end
  @inbounds hdl.pMin[aRoot] = minVal; hdl.pMax[aRoot] = maxVal;
  nothing
end

# Recursively evaluate the density implied by the samples of the
# subtree (rooted at dRoot) of densTree at the locations given by
# the subtree (rooted at aRoot) of *this, to within the error
# percentage "maxErr"
function evaluate(bd::BallTreeDensity, dRoot::Int64, atTree::BallTreeDensity, aRoot::Int64, maxErr::Float64, hdl::pArrHdls)

  #result = 0.0
  #tmp = 0.0
  restmp = Array{Float64,1}(2)

  # find the minimum and maximum effect of these two balls on each other
  minDistKer!(restmp, bd, dRoot, atTree, aRoot)
  Kmax = restmp[1]
  maxDistKer!(restmp, bd, dRoot, atTree, aRoot)
  Kmin = restmp[1]

  total = hdl.pMin[ aRoot ];		   	        # take pmin of data below this level
  #total += hdl.pAdd[1,aRoot] - hdl.pErr[aRoot]; # add lower bound from local expansion
  total += weight(bd, dRoot)*Kmin;              # also add minimum for this block

  # if the weighted contribution of this multiply is below the
  #    threshold, no need to recurse; just treat as constant
  if ( Kmax - Kmin <= maxErr * total)                     # APPROXIMATE: PERCENT
    Kmin *= weight(bd, dRoot); Kmax *= weight(bd, dRoot);

    if (bd == atTree && aRoot==dRoot)                  # LEAVE-ONE-OUT (and same subtree)
      for k in leafFirst(atTree, aRoot):(leafLast(atTree, aRoot))
        hdl.pMin[k] += Kmin * (1.0 - weight(bd, k)/weight(bd, dRoot))   # leave our weight out of it
        hdl.pMax[k] += Kmax * (1.0 - weight(bd, k)/weight(bd, dRoot))
      end
      recurseMinMax(atTree, aRoot, hdl)
    else                                                 #     NO L-O-O => just add away
      # hdl.pAdd[1,aRoot] += (Kmin + Kmax)/2.0; hdl.pErr[aRoot] = (Kmax-Kmin)/2.0;
      # !!! Should *not* do this -- instead add to local expansion (constant term)
      for k in leafFirst(atTree, aRoot):(leafLast(atTree, aRoot))
        hdl.pMin[k] += Kmin;
        hdl.pMax[k] += Kmax;
      end

      if !(isLeaf(atTree, aRoot))
          hdl.pMin[aRoot] += Kmin; hdl.pMax[aRoot] += Kmax
      end
    end

  else
    if (Npts(bd, dRoot)*Npts(atTree, aRoot)<=DirectSize)  # DIRECT EVALUATION
        evalDirect(bd, dRoot, atTree, aRoot, hdl)
    else
        # RECURSE ON SUBTREES

        # Find the subtree in closest to the other tree's left child and do
        # that first so that the values are higher and there is a better
        # chance of being able to skip a recursion.
        close_ = closer!(atTree, left(atTree, aRoot), right(atTree, aRoot), bd, left(bd, dRoot))
        if (left(bd, dRoot) != NO_CHILD && close_ != NO_CHILD)
          evaluate(bd, left(bd, dRoot), atTree, close_, maxErr, hdl);
        end
        far   = (close_ == left(atTree, aRoot)) ? right(atTree, aRoot) : left(atTree, aRoot);
        if (left(bd, dRoot) != NO_CHILD && far != NO_CHILD)
          evaluate(bd, left(bd, dRoot), atTree, far, maxErr, hdl);
        end

        # Now the same thing for the density's right child
        close_ = closer!(atTree, left(atTree, aRoot), right(atTree, aRoot), bd, right(bd, dRoot))
        if (right(bd, dRoot) != NO_CHILD && close_ != NO_CHILD)
          evaluate(bd, right(bd, dRoot) , atTree, close_, maxErr, hdl)
        end
        far   = (close_ == left(atTree, aRoot)) ? right(atTree, aRoot) : left(atTree, aRoot);
        if (right(bd, dRoot) != NO_CHILD && far != NO_CHILD)
          evaluate(bd, right(bd, dRoot), atTree, far, maxErr, hdl);
        end

        # Propogate additions in children's minimum value to this node
        if !(isLeaf(atTree, aRoot))
          hdl.pMin[aRoot] = hdl.pMin[ left(atTree, aRoot) ]
          hdl.pMax[aRoot] = hdl.pMax[ left(atTree, aRoot) ]
          if (right(atTree, aRoot) != NO_CHILD)
            if (hdl.pMin[aRoot] > hdl.pMin[ right(atTree, aRoot) ])
              hdl.pMin[aRoot] = hdl.pMin[ right(atTree, aRoot) ]
            end
            if (hdl.pMax[aRoot] < hdl.pMax[ right(atTree, aRoot) ])
              hdl.pMax[aRoot] = hdl.pMax[ right(atTree, aRoot) ]
            end
          end
        end
    end
  end#?? this one is extra

  return Union{}
end

# Dual Tree evaluation: estimate the values at this ball tree's
# points given the other tree as the samples from a distribution.
function evaluate(bd::BallTreeDensity, locations::BallTreeDensity, p::Array{Float64,1}, maxErr::Float64)
  #println("BallTreeDensity::evaluate(bd, locations, double*, maxerr) -- starting")

  if bd.bt.dims != locations.bt.dims
    error("evaluate -- dimensions of two BallTreeDensities must match")
  end
  if (p == Union{})
    error("evaluate -- density p must exist")
  end

  # TODO -- this internal declaration is making it slow!!
  hdl = pArrHdls(zeros(2*locations.bt.num_points),
                 zeros(2*locations.bt.num_points),
                 zeros(1,0),#zeros(1,2*locations.bt.dims),
                 zeros(0),#zeros(2*locations.bt.num_points),
                 [0.], [0.])

  evaluate(bd, root(), locations, root(), 2.0*maxErr, hdl)

  # Compute & account for the kernel f'ns normalization constant
  #norm = 1.0
  #switch(getType()) {
    #case Gaussian:
    norm = (2.0*pi)^((bd.bt.dims)/2.0)
    if (bwUniform(bd))
      for i in 1:bd.bt.dims
          norm *= sqrt(bd.bandwidthMax[i])
      end
    end

    #case Laplacian: norm = pow(2, ((double)Ndim()) );
    #                if (bwUniform())
    #                  for (unsigned int i=0;i<Ndim();i++) norm *= bandwidthMax[i];
    #                break;
    #case Epanetchnikov: norm = pow(4.0/3, ((double)Ndim()) );
    #                if (bwUniform())
    #                  for (unsigned int i=0;i<Ndim();i++) norm *= bandwidthMax[i];
    #                break;
  #}
  lRoot = root()
  #@show size(hdl.pMin), size(hdl.pMax)
  if (bd == locations)    # if we need to do leave-one-out
    for j in leafFirst(locations, lRoot):(leafLast(locations, lRoot))
      #println("if bd==locations j=$(j)")
      p[getIndexOf(locations, j)] = .5*(hdl.pMin[j]+hdl.pMax[j])/norm/(1-weight(bd, j))
    end
  else
    for j in leafFirst(locations, lRoot):(leafLast(locations, lRoot))
      #println("else bd==locations j=$(j)")
      p[getIndexOf(locations, j)] = .5*(hdl.pMin[j]+hdl.pMax[j])/norm;
    end
  end
  hdl.pMin = zeros(0)
  hdl.pMax = zeros(0)
  Union{}
end

function makeDualTree(bd1::BallTreeDensity, bd2::BallTreeDensity, errTol::Float64)
    #println("makeDualTree(::BTD,::BTD) -- UNTESTED!")
    #densTree = bd1
    #atTree = bd2
    pRes = zeros(Npts(bd2))
    evaluate(bd1, bd2, pRes, errTol)
    return pRes
end

function evaluateDualTree(bd::BallTreeDensity, pos::Array{Float64,2}, lvFlag::Bool=false, errTol::Float64=1e-3)
    #dim = size(pos,1)
    if (bd.bt.dims != size(pos,1)) error("bd and pos must have the same dimension") end
    if (lvFlag)
        p = makeDualTree(bd, errTol)
    else
      posKDE = makeBallTreeDensity(pos, ones(size(pos,2))./size(pos,2));
      p = makeDualTree(bd,posKDE,errTol)
    end
    return p
end
# should make this be the Union again TODO ??
function evaluateDualTree(bd::BallTreeDensity, pos::Union{Array{Float64,1},LinSpace{Float64}}, lvFlag::Bool=false, errTol::Float64=1e-3)
    pos2 = zeros(1,length(pos))
    pos2[1,:] = pos[:]
    return evaluateDualTree(bd, pos2, lvFlag, errTol)
end

function makeDualTree(bd::BallTreeDensity, errTol::Float64)
    #densTree = bd ## just a new pointer to the same data...
    pRes = zeros(Npts(bd))
    #println("makeDualTree -- leave one out, going to call evaluate(densTree")
    evaluate(bd, bd, pRes, errTol) # this seems excessive

    return pRes
end

function evaluateDualTree(bd::BallTreeDensity, pos::BallTreeDensity, lvFlag::Bool=false, errTol::Float64=1e-3)
    if (bd.bt.dims != pos.bt.dims) error("bd and pos must have the same dimension") end
    if (lvFlag)
        p = makeDualTree(bd, errTol)
    else
        posKDE = pos
        p = makeDualTree(dens,posKDE,errTol)
    end
    return p
end

function evalAvgLogL(bd1::BallTreeDensity, bd2::BallTreeDensity)
  L = evaluateDualTree(bd1, bd2, true)
  #printBallTree(bd1)
  W = getWeights(bd2)
  ind = find(L==0)
  if sum(find(W[ind])) > 0
    #println("evalAvgLogL -- in if")
    ll=-Inf
  else
    #println("evalAvgLogL -- in else")
    L[ind] = 1
    ll = (log(L)')*W
  end
  return ll
end

function evalAvgLogL(bd1::BallTreeDensity, at::Array{Float64,1})
    error("evalAvgLogL(bd1::BallTreeDensity, at::Array{Float64,1}) -- not implemented yet")
end

function entropy(bd::BallTreeDensity)
    H = -evalAvgLogL(bd,bd)
    return H[1]
end

function updateBandwidth!(bd::BallTreeDensity, bw::Array{Float64, 1})
    if (bd.multibandwidth==0)
        bd.bandwidth = bw
        bd.bandwidthMax = bd.bandwidthMin = bd.bandwidth[(bd.bt.num_points*bd.bt.dims+1):end]
    else
        error("updateBandwidth! -- multibandwidth==0 ELSE not implemented yet")
    end
end

function nLOO_LL(alpha::Float64, bd::BallTreeDensity)
  alpha = alpha.^2 # assume always a Gaussian
  updateBandwidth!(bd,bd.bandwidth*alpha)
  H = entropy(bd)
  updateBandwidth!(bd,bd.bandwidth/alpha)
  return H
end

function golden(bd::BallTreeDensity, ax::Float64, bx::Float64, cx::Float64, tol::Float64=1e-2)
#GOLDEN   Minimize function of one variable using golden section search
#
#   xmin, fmin = golden(npd, f, ax, bx, cx, tol) computes a local minimum
#   of f. xmin is the computed local minimizer of f and fmin is
#   f(xmin). xmin is computed to an relative accuracy of TOL.
#
#   The parameters ax, bx and cx must satisfy the following conditions:
#   ax < bx < cx, f(bx) < f(ax) and f(bx) < f(cx).
#
#   xmin satisfies ax < xmin < cx. golden is guaranteed to succeed if f
#   is continuous between ax and cx
#
#   Roman Geus, ETH Zuerich, 9.12.97

    C = (3.0-sqrt(5.0))/2.0;  R = 1.0-C;

    x0 = ax;  x3 = cx;
    if (abs(cx-bx) > abs(bx-ax))
      x1 = bx;
      x2 = bx + C*(cx-bx)
    else
      x2 = bx;
      x1 = bx - C*(bx-ax)
    end

    #@show x1, x2
    f1 = nLOO_LL(x1,bd)
    f2 = nLOO_LL(x2,bd)

    k = 1;
    #tic()
    while abs(x3-x0) > tol*(abs(x1)+abs(x2))
    #  fprintf(1,'k=%4d, |a-b|=%e\n', k, abs(x3-x0));
      if f2 < f1
        x0 = x1;  x1 = x2;
        x2 = R*x1 + C*x3   # x2 = x1+c*(x3-x1)
        f1 = f2;
        f2 = nLOO_LL(x2,bd)
      else
        x3 = x2; x2 = x1;
        x1 = R*x2 + C*x0   # x1 = x2+c*(x0-x2)
        f2 = f1
        f1 = nLOO_LL(x1,bd)
      end
      k += 1
    #  [x0,x1,x2,x3,f1,f2]
    end
    #println("Time for while loop $(toc())")

    if f1 < f2
      xmin = x1
      fmin = f1
    else
      xmin = x2
      fmin = f2
    end
    return xmin, fmin
end

function neighborMinMax(bd::BallTreeDensity)
    tmp = (2*bd.bt.ranges).^2
    rang = reshape(bd.bt.ranges[1:(floor(Int64,end/2.0))],bd.bt.dims,bd.bt.num_points)
    maxm = sqrt(sum( (2.0*rang[:,1]).^2 ))
    ssumt = sqrt(sum( (2.0*rang[:,1:(bd.bt.num_points-1)]).^2 ,1))
    minm = minimum(ssumt[:])
    minm = maximum([minm,1e-6]);
    return minm, maxm
end

function ksize(bd::BallTreeDensity, t::ASCIIString="lcv")
    # sticking with lcv to start
    Nd = bd.bt.dims; Np = bd.bt.num_points;
    ##if (t=="lcv" || 'unif','lcvp','unifp')
    minm,maxm = neighborMinMax(bd)
    p = kde!(getPoints(bd),[(minm+maxm)/2.0],getWeights(bd))#,getType(bd));
    ks, dummy =  golden(p,2.0*minm/(minm+maxm),1.0,2.0*maxm/(minm+maxm)) #(p,nLOO_LL
    ks = ks * ( minm + maxm )/2.0
    ##end
    npd = kde!(getPoints(p),[ks],getWeights(p))
    return npd
end

function kde!(points::Array{Float64,2}, autoselect::ASCIIString="lcv")
  p = kde!(points, [1.0])
  #BEFORE
  # p = ksize(p, autoselect)
  #AFTER with independent weights on each dimension
  dims = size(points,1)
  bwds = zeros(dims)
  for i in 1:dims
    pp = ksize(marginal(p,[i]), autoselect)
    # pp2 = kde!(samplePts[i,:], "lcv")
    bwds[i] = getBW(pp)[1]
  end
  p = kde!(points, bwds)

  return p
end

#function kde!(points::Array{Float64,2}, autoselect::ASCIIString="lcv")
#  p = kde!(points, [retrieveTCP( points)])
#  return p
#end

function kde!(points::Array{Float64,1}, autoselect::ASCIIString="lcv")
  return kde!(points',autoselect)
end

function getKDERange(bd::BallTreeDensity; extend::Float64=0.1)
  if (bd.bt.dims == 1)
    pts = getPoints(bd)
    rangeV = [minimum(pts),maximum(pts)]
  else
      return error("getKDERange(::BTD) -- multidimensional not implemented yet")
  end

  dr = extend*(rangeV[2]-rangeV[1])
  rangeV[1] = rangeV[1] - dr;
  rangeV[2] = rangeV[2] + dr;
  return rangeV
end

function getKDERangeLinspace(bd::BallTreeDensity; extend::Float64=0.1, N::Int=201)
  v = getKDERange(bd,extend=extend)
  return linspace(v[1],v[2],N)
end

function getKDEMax(p::BallTreeDensity;N=200)
  m = zeros(p.bt.dims)
  for i in 1:p.bt.dims
    mm = marginal(p,[i])
    rangeV = getKDERange(mm)
    X = linspace(rangeV[1],rangeV[2],N)
    yV = evaluateDualTree(mm,X)
    m[i] = X[findfirst(yV,maximum(yV))]
  end
  return m
end

function getKDEMean(p::BallTreeDensity)
  return vec(Base.mean(getPoints(p),2))
end

function intersIntgAppxIS(p::BallTreeDensity, q::BallTreeDensity;N=201)
  ndims = Ndim(p)
  xx = zeros(ndims, N)
  dx = zeros(ndims)
  LD = Array{LinSpace,1}(ndims)
  for d in 1:ndims
    LD[d] = getKDERangeLinspace(marginal(p,[d]), N=N, extend=0.3)
    dx[d] = LD[d][2]-LD[d][1]
  end
  xx[1,:] = collect(getKDERangeLinspace(marginal(p,[1]), N=N, extend=0.3))
  acc = 0.0
  if ndims == 1
    yy=evaluateDualTree(p,xx[:])
    yy= yy .* evaluateDualTree(q,xx[:])
    acc = sum(yy)*dx[1]
  elseif ndims == 2
    for i in 1:N
      for j in 1:N
        @inbounds xx[2,j] = LD[2][i]
      end
      yy=evaluateDualTree(p,xx)
      yy = yy .* evaluateDualTree(q,xx)
      acc += (dx[1]*sum(yy))*dx[2]
    end
  else
    error("Can't do higher dimensions yet")
  end
  return acc
end

function test03()
    spls = zeros(1,3)
    spls[1,:] = [0.5172, 0.7169, 0.4049]
    p1000 = kde!(spls,"lcv");
    printBallTree(p1000);

    # should get
    #D= 1
    #N= 3
    #centers= 0.5609     0.46105           0      0.4049      0.5172      0.7169
    #weights= 1     0.66667           0     0.33333     0.33333     0.33333
    #ranges= 0.156     0.05615           0           0           0           0
    #leftch= 1  3  0  3  4  5
    #rightch= 5           4           0  4294967295  4294967295  4294967295
    #highest= 5  4  0  3  4  5
    #lowest= 3  3  0  3  4  5
    #perm= 0  0  0  2  0  1
    #means= 0.54633     0.46105           0      0.4049      0.5172      0.7169
    #bw= 0.05517    0.041674           0    0.038521    0.038521    0.038521
end
function test04()
    spls = zeros(2,3)
    spls[1,:] = [0.5172, 0.7169, 0.4049]
    spls[2,:] = [0.0312, 1.0094, 2.0204]
    p1000 = kde!(spls,"lcv");
    printBallTree(p1000);
    # should get
    #BallTreeDensity
    #D= 2
    #N= 3
    #centers= 0.5609      1.0258     0.61705      0.5203           0           0      0.5172      0.0312      0.7169      1.0094      0.4049      2.0204
    #weights= 1     0.66667           0     0.33333     0.33333     0.33333
    #ranges= 0.156      0.9946     0.09985      0.4891           0           0           0           0           0           0           0           0
    #leftch= 1  3  0  3  4  5
    #rightch= 5           4           0  4294967295  4294967295  4294967295
    #highest= 5  4  0  3  4  5
    #lowest= 3  3  0  3  4  5
    #perm= 0  0  0  0  1  2
    #means= 0.54633      1.0203     0.61705      0.5203           0           0      0.5172      0.0312      0.7169      1.0094      0.4049      2.0204
    #bw= 1.0268      1.6697      1.0201      1.2494           0           0      1.0101      1.0101      1.0101      1.0101      1.0101      1.0101
end


text:
include("setup.jl")

@testset "ca_roots" begin
    withenv(
        "JULIA_SSL_CA_ROOTS_PATH" => nothing,
    ) do
        @test ca_roots_path() isa String
        @test ispath(ca_roots_path())
        if Sys.iswindows() || Sys.isapple()
            @test ca_roots_path() == BUNDLED_CA_ROOTS
            @test ca_roots() === nothing
        else
            @test ca_roots_path() != BUNDLED_CA_ROOTS
            @test ca_roots() == ca_roots_path()
        end
        unset = ca_roots(), ca_roots_path()
        value = "Why hello!"
        ENV["JULIA_SSL_CA_ROOTS_PATH"] = value
        @test ca_roots() == value
        @test ca_roots_path() == value
        ENV["JULIA_SSL_CA_ROOTS_PATH"] = ""
        @test ca_roots() == unset[1]
        @test ca_roots_path() == unset[2]
    end
end

@testset "verify_host" begin
    withenv(
        "JULIA_NO_VERIFY_HOSTS" => nothing,
        "JULIA_SSL_NO_VERIFY_HOSTS" => nothing,
        "JULIA_SSH_NO_VERIFY_HOSTS" => nothing,
    ) do
        @testset "verify everything" begin
            for url in TEST_URLS
                @test verify_host(url) # cover this API once
                for transport in (nothing, "ssl", "ssh", "xyz")
                    @test verify_host(url, transport)
                end
            end
            clear_vars!(ENV)
        end

        @testset "bad patterns fail safely" begin
            patterns = [
                "~", "* *", "*~*", "***", "∀", "~, ***",
                ".com", "*com", ".*com", ".example.com", "*example.com",
            ]
            for pattern in patterns
                ENV["JULIA_NO_VERIFY_HOSTS"] = pattern
                for url in TEST_URLS, transport in TRANSPORTS
                    verify = without_warnings() do
                        verify_host(url, transport)
                    end
                    @test verify
                end
            end
            clear_vars!(ENV)
        end

        @testset "only ignore bad patterns in list" begin
            patterns = ["ok.com,~", "^, ok.com ,, !"]
            for pattern in patterns
                ENV["JULIA_NO_VERIFY_HOSTS"] = pattern
                for url in TEST_URLS
                    verify = without_warnings() do
                        verify_host(url)
                    end
                    @test verify
                end
                verify = without_warnings() do
                    verify_host("ok.com")
                end
                @test !verify
        end
            clear_vars!(ENV)
        end

        @testset "verify nothing" begin
            for pattern in ["**", "example.com,**", "**,, blah"]
                ENV["JULIA_NO_VERIFY_HOSTS"] = pattern
                for url in TEST_URLS, transport in TRANSPORTS
                    @test !verify_host(url, transport)
                end
            end
            clear_vars!(ENV)
        end

        @testset "SSL no verify" begin
            for pattern in ["**", "example.com,**", "**, blah"]
                ENV["JULIA_SSL_NO_VERIFY_HOSTS"] = pattern
                for url in TEST_URLS, transport in TRANSPORTS
                    no_verify = transport in ["ssl", "tls"]
                    @test verify_host(url, transport) == !no_verify
                end
            end
            clear_vars!(ENV)
        end

        @testset "SSH no verify" begin
            for pattern in ["**", "example.com,**", "**, blah"]
                ENV["JULIA_SSH_NO_VERIFY_HOSTS"] = pattern
                for url in TEST_URLS, transport in TRANSPORTS
                    no_verify = transport == "ssh"
                    @test verify_host(url, transport) == !no_verify
                end
            end
            clear_vars!(ENV)
        end

        @testset "complex scenario" begin
            ENV["JULIA_NO_VERIFY_HOSTS"] = "**.safe.example.com"
            ENV["JULIA_SSL_NO_VERIFY_HOSTS"] = "ssl.example.com"
            ENV["JULIA_SSH_NO_VERIFY_HOSTS"] = "ssh.example.com,*.ssh.example.com"
            for transport in TRANSPORTS
                for url in TEST_URLS
                    @test verify_host(url, transport)
                end
                hosts = [
                    "safe.example.com",
                    "api.SAFE.example.COM",
                    "v1.API.safe.eXample.com",
                ]
                for host in hosts, (url, valid) in host_variants(host)
                    @test verify_host(url, transport) == !valid
                end
                hosts = [
                    "ssl.example.com",
                    "SSL.example.com",
                    "ssl.Example.COM",
                ]
                for host in hosts, (url, valid) in host_variants(host)
                    no_verify = valid && transport in ["ssl", "tls"]
                    @test verify_host(url, transport) == !no_verify
                end
                hosts = [
                    "sub.ssl.example.com",
                    "sub.SSL.example.com",
                    "ssl..example.com",
                ]
                for host in hosts, (url, valid) in host_variants(host)
                    @test verify_host(url, transport)
                end
                hosts = [
                    "ssh.example.com",
                    "ssh.EXAMPLE.com",
                    "sub.ssh.example.com",
                    "sub.ssh.example.COM",
                ]
                for host in hosts, (url, valid) in host_variants(host)
                    no_verify = valid && transport == "ssh"
                    @test verify_host(url, transport) == !no_verify
                end
                hosts = [
                    "v1.api.ssh.example.com",
                    "123.api.SSH.example.COM",
                ]
                for host in hosts, (url, valid) in host_variants(host)
                    @test verify_host(url, transport)
                end
            end
            clear_vars!(ENV)
        end
    end
end


text:
@testset "Ellipse" begin

    img = zeros(Gray{N0f8}, 10, 10)
    expected = copy(img)
    res = @inferred draw(img, Ellipse(Point(5,5), 5, 1))
    expected[5, 1:9] .= 1
    @test all(expected .== res) == true

    res = @inferred draw(img, Ellipse(5, 5, 5, 5), oneunit(Gray{N0f8}))
    expected = Gray{N0f8}[  0 0 1 1 1 1 1 0 0 0
                            0 1 1 1 1 1 1 1 0 0
                            1 1 1 1 1 1 1 1 1 0
                            1 1 1 1 1 1 1 1 1 0
                            1 1 1 1 1 1 1 1 1 0
                            1 1 1 1 1 1 1 1 1 0
                            1 1 1 1 1 1 1 1 1 0
                            0 1 1 1 1 1 1 1 0 0
                            0 0 1 1 1 1 1 0 0 0
                            0 0 0 0 0 0 0 0 0 0 ]
    @test all(expected .== res) == true

    res = @inferred draw(img, Ellipse(CirclePointRadius(5, 5, 5)))
    @test all(expected .== res) == true

    res = @inferred draw(img, Ellipse(5, 5, 5, 5), Gray{N0f8}(0.5))
    expected = map(i->(i==1 ? Gray{N0f8}(0.5) : Gray{N0f8}(0)), expected)
    @test all(expected .== res) == true

    res = @inferred draw(img, Ellipse(CartesianIndex(5,5), 5, 3))
    expected = Gray{N0f8}[  0 0 0 0 0 0 0 0 0 0
                            0 0 0 0 0 0 0 0 0 0
                            0 1 1 1 1 1 1 1 0 0
                            1 1 1 1 1 1 1 1 1 0
                            1 1 1 1 1 1 1 1 1 0
                            1 1 1 1 1 1 1 1 1 0
                            0 1 1 1 1 1 1 1 0 0
                            0 0 0 0 0 0 0 0 0 0
                            0 0 0 0 0 0 0 0 0 0
                            0 0 0 0 0 0 0 0 0 0 ]
    @test all(expected .== res) == true

    img = zeros(Gray{N0f8}, 10, 10)
    expected = Gray{N0f8}[ 0 0 0 0 0 0 0 0 0 0
                           0 0 0 1 1 1 1 1 0 0
                           0 0 1 1 1 1 1 1 1 0
                           0 1 1 1 1 1 1 1 1 1
                           0 1 1 1 1 0 1 1 1 1
                           0 1 1 1 0 0 0 1 1 1
                           0 1 1 1 1 0 1 1 1 1
                           0 1 1 1 1 1 1 1 1 1
                           0 0 1 1 1 1 1 1 1 0
                           0 0 0 1 1 1 1 1 0 0 ]
    res = @inferred draw(img, Ellipse(CirclePointRadius(Point(6, 6), 5; thickness = 4, fill = false)))
    @test all(expected .== res)

    err = ArgumentError("Thickness 7 should be smaller than 5.")
    @test_throws err Ellipse(CartesianIndex(6, 6), 5, 5; thickness = 7, fill = false)
end


text:
module Car1DExample
export CarX, CarZ, CarU, car_system, WallObsDynamics, LeaderFollowerControl

using OneVision
using OneVision: ℝ, 𝕋, ℕ, @kwdef, @unpack, @asserteq
using OneVision.Examples
using Random
using StaticArrays
using Plots
using Statistics: mean

import OneVision

@kwdef struct CarX{R} <: FieldVector{3,R}
    pos::R
    velocity::R
    "Acceleration from the last time step"
    acc::R
end

StaticArrays.similar_type(::Type{CarX{A}}, ::Type{B}) where {A, B} = CarX{B}

struct CarU{R} <: FieldVector{1,R}
    acc::R
end

StaticArrays.similar_type(::Type{CarU{A}}, ::Type{B}) where {A, B} = CarU{B}

CarU{X}(x::Tuple{X}) where X = CarU(x[1])

struct CarZ{R}
    detected::Bool
    distance::R
end

struct CarSystem{Noise} <: SysDynamics
    dt::ℝ
    mass::ℝ
    v_noise::Noise
end

function car_system(dt::ℝ, v_noise::Function = _ -> 0.0; mass::ℝ = 1.0)
    CarSystem(dt, mass, v_noise)
end

function OneVision.sys_forward(dy::CarSystem, x::CarX{R}, u, t::𝕋)::CarX{R} where R
    @unpack dt, mass, v_noise = dy
    CarX(
        pos = x.pos + x.velocity * dt + 0.5 * u.acc / mass * dt^2,
        velocity = x.velocity + u.acc / mass * dt + v_noise(t),
        acc = u.acc,
    )
end

"The wall observation model used for simulation.\n"
@kwdef struct WallObsDynamics <: ObsDynamics
    wall_position::Union{ℝ,Nothing}
    detector_range::ℝ
end

function OneVision.obs_forward(
    dy::WallObsDynamics, x::CarX, z::CarZ{R}, t::𝕋
)::CarZ{R} where R
    if z.detected
        z
    elseif (isnothing(dy.wall_position)
            || x.pos + dy.detector_range < dy.wall_position)
        CarZ(false, 0.0)
    else
       CarZ(true, dy.wall_position)
    end
end

"""
The simple wall dynamics model used by OneVision.
"""
struct WallObsModel <: ObsDynamics end

OneVision.obs_forward(dy::WallObsModel, x::CarX, z::CarZ, t::𝕋) = z

"""
A controller struct - holds all necessary information
"""
@kwdef struct LeaderFollowerControl <: CentralControlStateless{CarU{ℝ}}
    "Will output u=0 before this time"
    warm_up_time::𝕋
    "Will only use bang bang control when difference is larger than this"
    bang_bang_tol::ℝ = Inf
    k_v::ℝ = 3.0
    k_x::ℝ = 2.0
    stop_distance::ℝ = 4.0
    target_v::ℝ = 2.0
end

OneVision.control_one(
    lf::LeaderFollowerControl, xs, zs, t::𝕋, id::ℕ
)::CarU{ℝ} = begin
    tol = lf.bang_bang_tol
    function bang_bang(x̂, x, k, tol)
        if abs(x̂ - x) ≤ tol
            (x̂ - x) * k
        else
            sign(x̂ - x) * tol * k
        end
    end

    x, z = xs[id], zs[id]
    if t ≤ lf.warm_up_time
        acc = 0.0
    elseif z.detected && z.distance - x.pos ≤ lf.stop_distance
        acc = bang_bang(0.0, x.velocity, lf.k_v, tol)
    elseif id == 1
        # the leader
        acc = bang_bang(lf.target_v, x.velocity, lf.k_v, tol)
    else
        # a follower
        leader = xs[1]
        acc = (leader.acc # mimic the leader
                + bang_bang(leader.velocity, x.velocity, lf.k_v, Inf)
                + bang_bang(leader.pos, x.pos, lf.k_x, Inf))
    end
    CarU(acc)
end

"""
Compute the time-averaged distance between the leader and the follower
using a quadratic norm, i.e., √(1/T ∫ (x_1 - x_2)^2 dt) .
"""
function avg_distance(result::TrajectoryData)::ℝ
    pos_data = result["pos"]
    @asserteq size(pos_data, 2) 2 "there should only be 2 agents"
    dis = (view(pos_data, :, 1) .- view(pos_data, :, 2)).^2
    √(mean(dis))
end

function run_example(;
        setting::ExampleSetting,
        CF::CFName = onevision_cf,
        has_obstacle = true, use_bang_bang = true,
        seed = 1234,
        plot_result = true, log_prediction = false)
    @unpack time_end, freq, noise, sensor_noise, delays, H, model_error = setting
    t_end = 𝕋(ceil(time_end * freq))
    delta_t = 1 / freq
    times::Vector{𝕋} = collect(1:t_end)
    idx_to_time(xs) = (xs .- 1) .* delta_t
    rng = MersenneTwister(seed)

    agent_info(id) = begin
        v_noise = randn(rng, ℝ, t_end) * noise
        sys_dy = car_system(delta_t, t -> v_noise[t])
        obs_dy =
            if id == 1 && has_obstacle; WallObsDynamics(wall_position = 30.0, detector_range = 8.0)
            else WallObsModel() end

        sys_dy, obs_dy
    end

    N = 2
    world_dynamics = WorldDynamics([agent_info(i) for i in 1:N])
    init_states = fill((CarX(0.0, 0.0, 0.0), CarZ(false, 0.0), CarU(0.0)), N)
    ΔT = delays.ΔT
    world_model = WorldDynamics(
        fill((car_system(delta_t; mass = 1.0 / (1 + model_error)), WallObsModel()), N))
    function xs_observer(xs, t)
        (x -> x + randn(rng, CarX{ℝ}) * sensor_noise).(xs)
    end

    comps = ["pos", "velocity", "acceleration", "obstacle"]
    function record_f(xs, zs, us)
        [(x -> x.pos).(xs) (x -> x.velocity).(xs) (u -> u.acc).(us) (z -> z.distance).(zs)]
    end

    central = LeaderFollowerControl(
        warm_up_time = delays.act,
        bang_bang_tol = use_bang_bang ? 1.0 : Inf)
    x_weights = SVector{N}(fill(CarX(100.0, 10.0, 0.0), N))
    u_weights = SVector{N}(fill(CarU(1.0), N))
    loss_model = RegretLossModel(central, world_model, x_weights, u_weights)
    framework = mk_cf(
        CF, world_model, central, delays, loss_model;
        X = CarX{ℝ}, Z = CarZ{ℝ}, H)

    result, (logs, loss) = simulate(
        world_dynamics,
        delays,
        framework,
        init_states,
        (comps, record_f),
        times;
        loss_model,
        xs_observer,
    )

    if log_prediction
        # this is out-dated
        t0, t1 = (0.0, 11.0)
        id = 1
        x_indices = [t for t in result.times if t0 ≤ (t - 1) * delta_t ≤ t1]
        sorted_log = sort!([x for x in logs[id]], by = x -> x[1])
        δv_x = idx_to_time((x -> x[1]).(sorted_log))
        anim = @gif for (t, log) in sorted_log
            ps = []
            for (c, label) in enumerate(["δx", "δv"])
                δv_y = [log.δxz[1][c] for (_, log) in sorted_log]
                p = plot(idx_to_time(x_indices),
                        result.values[x_indices,:,c], label = ["agent 1" "agent 2"])
                plot!(p, δv_x, δv_y, label = label)
                t -= 1
                local ts = idx_to_time(t:t + size(log.x̃, 1) - 1)
                plot!(p, ts, (x -> x[c]).(log.x̃); label = ["central 1" "central 2"])
                push!(ps, p)
            end
            plot(ps...; layout = (2, 1))
        end
        anim |> display
    end
    plot_result && display(visualize(result; delta_t, loss))

    loss, Dict("avg distance" => avg_distance(result))
end

end # Car1DExampleLabeled


text:
using GNSSTools
using SatelliteToolbox
using ProgressMeter
using Makie
using GLMakie
using FileIO
# using GLMakie: Scene, surface!, scatter!, lines!, update_cam!, save
using Colors, GeometryBasics
# using AbstractPlotting, AbstractPlotting.MakieLayout

"""
Constellation info from
[GPS Performance standard](https://www.gps.gov/technical/ps/2020-SPS-performance-standard.pdf)
[USCG NAVCEN Constellation Status](https://www.navcen.uscg.gov/?Do=constellationStatus)
"""

GM = 3.986004418e14  # m³s⁻²
eop = get_eop();
scale_a = Rₑ + 35786*1000

a = Rₑ+20000*1000;
plane_num = 5;
sat_per_plane = 6;
incl = 56;
ΔΩ = 360/plane_num;
orbital_period = 2π*sqrt(a^3 / GM)  # seconds

Δt = 1 # seconds
t_range = Array(0:Δt:orbital_period);

constellation = define_constellation(a, plane_num,
                                     sat_per_plane, incl,
                                     t_range; obs_lla=missing,
                                     print_steps=true,
                                     eop=eop,
                                     ΔΩ=ΔΩ,
                                     show_plot=false,
                                     Ω₀=-30);

# Earth sphere
u = Array(range(0, 2π, length=200)) .- (55*pi/180);
v = Array(range(0, 1π, length=200));
u, v = meshgrid(u, v);
E_x = Rₑ.*cos.(u).*sin.(v);
E_y = Rₑ.*sin.(u).*sin.(v);
E_z = Rₑ.*cos.(v);

# Earth imagery
# earth = load(download("https://svs.gsfc.nasa.gov/vis/a000000/a002900/a002915/bluemarble-2048.png"));
earth = load("bluemarble-2048.png");

# View angle
ϕ = deg2rad(-80)
θ = deg2rad(80)
scaler = 3
camera_pos = (scaler*scale_a).*[cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)]

# Parameters
light_ϕ = deg2rad(150+90)
light_θ = deg2rad(90)
light_scaler = 149597870700/Rₑ
light_camera_pos = (light_scaler*scale_a).*[cos(light_ϕ)*sin(light_θ), sin(light_ϕ)*sin(light_θ), cos(light_θ)]
lightpos = Vec3f0(light_camera_pos...)
resolution = (1000, 1000)
save_resolution_scale = 2
pt_scale = 2
pt_size = (save_resolution_scale*pt_scale) * 9e8
strokewidth = (save_resolution_scale*pt_scale) * 3
linewidth = (save_resolution_scale*pt_scale) * 2
linestyle = nothing
pt_color = :lightgrey
line_color = :lightgrey
background_color = :white

function rotate_earth(i)
    to_eci = rECEFtoECI(ITRF(), TEME(), constellation.epoch+constellation.t_range[i], eop)
    # Transform Earth and CO from ECEF to ECI frames
    E_xt = zeros(size(E_x))
    E_yt = zeros(size(E_y))
    E_zt = zeros(size(E_z))
    for j in 1:size(E_x)[1]
        for k in 1:size(E_x)[2]
            x, y, z = to_eci*[E_x[j,k], E_y[j,k], E_z[j,k]]
            E_xt[j,k] = x
            E_yt[j,k] = y
            E_zt[j,k] = z
        end
    end
    return (E_xt, E_yt, E_zt)
end

function get_obj_positions(constellation)
    i = 1
    sat_num = length(constellation.satellites)
    to_eci = rECEFtoECI(ITRF(), TEME(), constellation.epoch+constellation.t_range[i], eop)
    sat_pos = zeros(length(constellation.satellites), 3)
    planes = []
    for j in 1:length(constellation.satellites)
        t = constellation.satellites[j].t[i]
        to_eci = rECEFtoECI(ITRF(), TEME(), t, eop)
        r_ecef = constellation.satellites[j].r_ecef[i,:]
        # r_eci = kepler_to_sv(constellation.satellites[j].orbit[i]).r
        r_eci = to_eci*r_ecef
        sat_pos[j,:] = r_eci
        if (constellation.satellites[j].id%sat_per_plane) == 1
            r_ecis = zeros(size(constellation.satellites[j].r_ecef))
            for k in 1:size(r_ecis)[1]
                t = constellation.satellites[j].t[k]
                to_eci = rECEFtoECI(ITRF(), TEME(), t, eop)
                # r_ecis[k,:] =  kepler_to_sv(constellation.satellites[j].orbit[k]).r
                r_ecis[k,:] = to_eci*constellation.satellites[j].r_ecef[k,:]
            end
            push!(planes, r_ecis)
        end
    end
    return (sat_pos, planes)
end

# Create scene
scene = Scene(resolution=resolution, backgroundcolor=background_color)
E_xt, E_yt, E_zt = rotate_earth(1)
surf = surface!(scene, E_xt, E_yt, E_zt, shading=true, color=earth,
                show_axis=false, #lightposition=lightpos,
                colormap=(:white, :white),
                ambient = Vec3f0(0.5, 0.5, 0.5),
                diffuse = Vec3f0(0.4, 0.4, 0.2),
                specular = Vec3f0(0.5, 0.5, 0.5),
                shininess = 0.5f0)
sat_pos, planes = get_obj_positions(constellation)
scatter_sat = scatter!(scene, sat_pos[:,1], sat_pos[:,2], sat_pos[:,3],
                       strokewidth=strokewidth, strokecolor=pt_color)
for j in 1:length(planes)
    line_planes = lines!(scene, planes[j][:,1], planes[j][:,2],
                         planes[j][:,3],
                         color=line_color, linewidth=linewidth,
                         linestyle=linestyle)
end

# GEO satellites
a = Rₑ + 35785*1000
Ωs = [-117.0, -107.3, -129.0]  # degrees
eccentricity = 0.0
inclination = 0.0
true_anomaly = 0.0
arg_of_perigee = 0.0
orbital_period_geo = 2π*sqrt(a^3 / GM)  # seconds
t_range_geo = Array(0:Δt:orbital_period_geo);  # seconds
geo_pt_color = :lightblue
geo_line_color = geo_pt_color
scatter_geo = zeros(length(Ωs), 3)
for i in 1:length(Ωs)
    Ω = deg2rad(Ωs[i])
    init_orbit = init_orbit_propagator(Val(:twobody), 0., a, eccentricity,
                                       inclination, Ω, arg_of_perigee,
                                       true_anomaly)
    r, v = propagate_to_epoch!(init_orbit, t_range_geo)
    scatter_geo[i,:] = r[1]  # ECI TEME
    if i == 1
        plane_geo = Array{Float64}(undef, length(r), 3)
        for j in 1:length(r)
            plane_geo[j,:] = r[j]
        end
        line_planes_geo = lines!(scene, plane_geo[:,1], plane_geo[:,2],
                                 plane_geo[:,3],
                                 color=geo_line_color, linewidth=linewidth,
                                 linestyle=linestyle)
    end
end
scatter_sat = scatter!(scene, scatter_geo[:,1], scatter_geo[:,2], scatter_geo[:,3],
                       strokewidth=strokewidth, strokecolor=geo_pt_color)

display(scene)
update_cam!(scene, cameracontrols(scene), Vec3f0(camera_pos), Vec3f0(0), Vec3f0(0, 0, 1))

scene.center=false;
save("gps.png", scene, resolution=(save_resolution_scale*resolution[1],
     save_resolution_scale*resolution[2]));


text:
module VarianceComponentTest

import Distributions
if VERSION < v"0.4-"
       using Docile
end

export
       readgeno!,
       readAnnotate!,
       kinshipcoef,
       vctestnullsim,
       vctest,
       loopwinFixsize,
       loopwinAnnot,
       gwasvctest

include("readgeno!.jl")
include("readAnnotate!.jl")
include("kinshipcoef.jl")
include("vctestnullsim.jl")
include("vctest.jl")
include("loopwinFixsize.jl")
include("loopwinAnnot.jl")
include("gwasvctest.jl")

end


text:
benchmark_suite = true

global N = 0

for n = 2:16
    N = n
    include("benchmark2.jl")
end

benchmark_suite = false;


text:
using Distributions

# params (computed from empirical data)
LAMBDA = 0.6882627
ALPHA = 1.3441739239406807
SIGMA = 0.4203008

mutable struct Model
    post_list
    top_list
    tick
    n_users
end

mutable struct Post
    id::Int64
    quality::Float64
    votes::Int64
    views::Int64
    timestamp::Int64
    score::Float64
end

function generate_attention(n)
    pois = Poisson(5)
    return rand(pois, n)
end

# Wie legen wir n fest?
# Aus den Daten holen (wie viele Votes werde pro Minute abgegeben?)

# refactor
function user_interaction!(post_list, attention_list)
    for att in attention_list
        for i in 1:att
            post_list[i].views += 1
            if post_list[i].quality > 0.5
                post_list[i].votes += 1
            end
        end
    end
    return post_list
end

function score_single(post, tick)
    age = (tick - post.timestamp) / 60
    votes = post.votes
    return ((votes - 1)^0.8)/((age + 2)^1.8)
end

function hnscore!(post_list, tick)
    for p in post_list
        p.score = score_single(p, tick)
    end
    return post_list
end

function sort_toplist(post_list)
    toplist = sort(post_list, rev=true, by=(post -> post.score))
    return toplist
end

function bias_lambda(slice, alpha=ALPHA, sd=SIGMA)
    sd * (- alpha * sin(2 * pi  * slice))
end

function draw_arrival_rate(tick, ticks_per_day, arrival_rate_mean=LAMBDA; alpha=ALPHA, arrival_rate_sd=SIGMA)
    rel_slice = mod(tick, ticks_per_day)
    lambda = arrival_rate_mean + bias_lambda(rel_slice, alpha, arrival_rate_sd)
    dist = Poisson(lambda)
    return rand(dist)
end

function tick()

end

l = [draw_arrival_rate(100, 1440) for i in 1:100]

histogram(l)

using Plots

post_list = [Post(rand(Uniform(0, 1)), 1, 0, 6, 0) for i in 1:100]
attention_list = generate_attention(100)
user_interaction!(post_list, attention_list)
hnscore!(post_list, 8)
toplist = sort_toplist(post_list)

post_list

# votes + age -> wie verhalten sie sich in Bezug auf score?


text:
coll = [1:50]
for e in coll
   # process(e) executed for every element e in coll
end

for n = 1:10  println(n^3)  end
# prints:
#1
#8
#27
#64
#125
#216
#343
#512
#729
#1000

# iterate over index:
arr = [x^2 for x in 1:10]
#10-element Array{Int64,1}:
#   1
#   4
#   9
#  16
#  25
#  36
#  49
#  64
#  81
# 100

for i = 1:length(arr)
   println("the $i-th element is $(arr[i])")
end
# the 1-th element is 1
# the 2-th element is 4
# the 3-th element is 9
# the 4-th element is 16
# the 5-th element is 25
# the 6-th element is 36
# the 7-th element is 49
# the 8-th element is 64
# the 9-th element is 81
# the 10-th element is 100

for n = 1:5, m = 1:5
    println("$n * $m = $(n * m)")
end
#
#1 * 1 = 1
#1 * 2 = 2
#1 * 3 = 3
#1 * 4 = 4
#1 * 5 = 5
#2 * 1 = 2
#2 * 2 = 4
#2 * 3 = 6
#2 * 4 = 8
#2 * 5 = 10
#3 * 1 = 3
#3 * 2 = 6
#3 * 3 = 9
#3 * 4 = 12
#3 * 5 = 15
#4 * 1 = 4
#4 * 2 = 8
#4 * 3 = 12
#4 * 4 = 16
#4 * 5 = 20
#5 * 1 = 5
#5 * 2 = 10
#5 * 3 = 15
#5 * 4 = 20
#5 * 5 = 25

# while loop:
a = 10
b = 15
while a < b
    # process(a)
    println(a)
    a += 1
end
# prints:
#10
#11
#12
#13
#14

# loop over array that is changing in the loop:
arr = [1,2,3,4]
while !isempty(arr)
    print(pop!(arr), ", ")
end
# => 4, 3, 2, 1,
# julia> arr
# 0-element Array{Int64,1}

# break:
a = 10
b = 150
while a < b
    # process(a)
    println(a)
    a += 1
    if a >= 50
        break
    end
end

# continue:
for n in 1:10
    if 3 <= n <= 6
        continue # skip one iteration
    end
    println(x)
end

max_stars_repo_path: Module 1/Chapter04/Old/repetitions.jl

max_issues_repo_path: Module 1/Chapter04/Old/repetitions.jl

max_forks_repo_path: Module 1/Chapter04/Old/repetitions.jl


text:
"""
    latitudes(layer::T) where {T <: SimpleSDMLayer}

Returns an iterator with the latitudes of the SDM layer passed as its argument.
This returns the latitude at the center of each cell in the grid.
"""
latitudes(layer::T) where {T <: SimpleSDMLayer} = LinRange(layer.bottom+stride(layer, 2), layer.top-stride(layer, 2), size(layer,1))

"""
    longitudes(layer::T) where {T <: SimpleSDMLayer}

Returns an iterator with the longitudes of the SDM layer passed as its argument.
This returns the longitudes at the center of each cell in the grid.
"""
longitudes(layer::T) where {T <: SimpleSDMLayer} = LinRange(layer.left+stride(layer, 1), layer.right-stride(layer, 1), size(layer,2))

"""
    _layers_are_compatible(l1::X, l2::Y) where {X <: SimpleSDMLayer, Y <: SimpleSDMLayer}

    Internal function to verify if layers are compatible, i.e. have the same
    size and bounding coordinates.

"""
function _layers_are_compatible(l1::X, l2::Y) where {X <: SimpleSDMLayer, Y <: SimpleSDMLayer}
    size(l1) == size(l2) || throw(ArgumentError("The layers have different sizes"))
    l1.top == l2.top || throw(ArgumentError("The layers have different top coordinates"))
    l1.left == l2.left || throw(ArgumentError("The layers have different left coordinates"))
    l1.bottom == l2.bottom || throw(ArgumentError("The layers have different bottom coordinates"))
    l1.right == l2.right || throw(ArgumentError("The layers have different right coordinates"))
end

function _layers_are_compatible(layers::Array{T}) where {T <: SimpleSDMLayer}
    all(layer -> _layers_are_compatible(layer, layers[1]), layers)
end


text:
"""
$(TYPEDSIGNATURES)

Generate `ODESystem`, dependent variables, and parameters from an `ODEProblem`.
"""
function modelingtoolkitize(prob::DiffEqBase.ODEProblem)
    prob.f isa DiffEqBase.AbstractParameterizedFunction &&
                            return prob.f.sys
    @parameters t

    if prob.p isa Tuple || prob.p isa NamedTuple
        p = [x for x in prob.p]
    else
        p = prob.p
    end

    has_p = !(p isa Union{DiffEqBase.NullParameters,Nothing})

    var(x, i) = Num(Sym{FnType{Tuple{symtype(t)}, Real}}(nameof(Variable(x, i))))
    vars = ArrayInterface.restructure(prob.u0,[var(:x, i)(ModelingToolkit.value(t)) for i in eachindex(prob.u0)])
    params = has_p ? reshape([Num(Sym{Real}(nameof(Variable(:α, i)))) for i in eachindex(p)],size(p)) : []
    var_set = Set(vars)

    D = Differential(t)
    mm = prob.f.mass_matrix

    if mm === I
        lhs = map(v->D(v), vars)
    else
        lhs = map(mm * vars) do v
            if iszero(v)
                0
            elseif v in var_set
                D(v)
            else
                error("Non-permuation mass matrix is not supported.")
            end
        end
    end

    if DiffEqBase.isinplace(prob)
        rhs = similar(vars, Num)
        prob.f(rhs, vars, params, t)
    else
        rhs = prob.f(vars, params, t)
    end

    eqs = vcat([lhs[i] ~ rhs[i] for i in eachindex(prob.u0)]...)

    sts = vec(collect(vars))
    params = if ndims(params) == 0
        [params[1]]
    else
        vec(collect(params))
    end
    default_u0 = Dict(sts .=> vec(collect(prob.u0)))
    default_p = has_p ? Dict(params .=> vec(collect(prob.p))) : Dict()

    de = ODESystem(
        eqs, t, sts, params,
        defaults=merge(default_u0, default_p),
    )

    de
end

"""
$(TYPEDSIGNATURES)

Generate `SDESystem`, dependent variables, and parameters from an `SDEProblem`.
"""
function modelingtoolkitize(prob::DiffEqBase.SDEProblem)
    prob.f isa DiffEqBase.AbstractParameterizedFunction &&
                            return (prob.f.sys, prob.f.sys.states, prob.f.sys.ps)
    @parameters t
    if prob.p isa Tuple || prob.p isa NamedTuple
        p = [x for x in prob.p]
    else
        p = prob.p
    end
    var(x, i) = Num(Sym{FnType{Tuple{symtype(t)}, Real}}(nameof(Variable(x, i))))
    vars = ArrayInterface.restructure(prob.u0,[var(:x, i)(ModelingToolkit.value(t)) for i in eachindex(prob.u0)])
    params = p isa DiffEqBase.NullParameters ? [] :
             reshape([Num(Sym{Real}(nameof(Variable(:α, i)))) for i in eachindex(p)],size(p))

    D = Differential(t)

    rhs = [D(var) for var in vars]

    if DiffEqBase.isinplace(prob)
        lhs = similar(vars, Any)

        prob.f(lhs, vars, params, t)

        if DiffEqBase.is_diagonal_noise(prob)
            neqs = similar(vars, Any)
            prob.g(neqs, vars, params, t)
        else
            neqs = similar(vars, Any, size(prob.noise_rate_prototype))
            prob.g(neqs, vars, params, t)
        end
    else
        lhs = prob.f(vars, params, t)
        if DiffEqBase.is_diagonal_noise(prob)
            neqs = prob.g(vars, params, t)
        else
            neqs = prob.g(vars, params, t)
        end
    end
    deqs = vcat([rhs[i] ~ lhs[i] for i in eachindex(prob.u0)]...)

    params = if ndims(params) == 0
        [params[1]]
    else
        Vector(vec(params))
    end

    de = SDESystem(deqs,neqs,t,Vector(vec(vars)),params)

    de
end

"""
$(TYPEDSIGNATURES)

Generate `OptimizationSystem`, dependent variables, and parameters from an `OptimizationProblem`.
"""
function modelingtoolkitize(prob::DiffEqBase.OptimizationProblem)

    if prob.p isa Tuple || prob.p isa NamedTuple
        p = [x for x in prob.p]
    else
        p = prob.p
    end

    vars = reshape([Num(Sym{Real}(nameof(Variable(:x, i)))) for i in eachindex(prob.u0)],size(prob.u0))
    params = p isa DiffEqBase.NullParameters ? [] :
             reshape([Num(Sym{Real}(nameof(Variable(:α, i)))) for i in eachindex(p)],size(Array(p)))

    eqs = prob.f(vars, params)
    de = OptimizationSystem(eqs,vec(vars),vec(params))
    de
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Stack <: Edit
        col::Array{Symbol,1}
        on::Array{Symbol,1}
        var::Symbol
    end

Concatenate side-by-side DataFrames into one normal-form DataFrame.

# Arguments
- `col::Array{Symbol,1}`: final column names
- `on::Array{Symbol,1}`: column name indicator specifying where to stack
- `var::Symbol`: column name for storing indicator
"""
mutable struct Stack <: Edit
    "final column names"
    col::Array{Symbol,1}
    "column name indicator specifying where to stack"
    on::Array{Symbol,1}
    "column name for storing indicator"
    var::Symbol
end

function Stack(; col, on, var, )
    Stack(col, on, var, )
end

"""Get [`Stack`](@ref) `col`."""
get_col(value::Stack) = value.col
"""Get [`Stack`](@ref) `on`."""
get_on(value::Stack) = value.on
"""Get [`Stack`](@ref) `var`."""
get_var(value::Stack) = value.var

"""Set [`Stack`](@ref) `col`."""
set_col!(value::Stack, val) = value.col = val
"""Set [`Stack`](@ref) `on`."""
set_on!(value::Stack, val) = value.on = val
"""Set [`Stack`](@ref) `var`."""
set_var!(value::Stack, val) = value.var = val

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
"""
Baseclass for all penalties.
"""
abstract type Penalty end

"""
Penalties that are applied element-wise.
"""
abstract type ElementPenalty <: Penalty end
abstract type ProxableElementPenalty <: ElementPenalty end

"""
Penalties that are applied to the entire parameter array.
"""
abstract type ArrayPenalty <: Penalty end


text:
############################################################
# joStack ##############################################
############################################################

export joStack, joStackException

struct joStackException <: Exception
    msg :: String
end

############################################################
## outer constructors

"""
    julia> op = joStack(op1[,op2][,...];[weights=...,][name=...])

Stack operator (single block column) composed from different JOLI operators

# Signature

    joStack(ops::joAbstractLinearOperator...;
        weights::LocalVector{WDT}=zeros(0),name::String="joStack")
            where {WDT<:Number}

# Arguments

- `op#`: JOLI operators (subtypes of joAbstractLinearOperator)
- keywords
    - `weights`: vector of waights for each operator
    - `name`: custom name

# Notes
- all operators must have the same # of columns (N)
- all given operators must have same domain/range types
- the domain/range types of joStack are equal to domain/range types of the given operators

# Example

define operators

    a=rand(ComplexF64,4,4);
    A=joMatrix(a;DDT=ComplexF32,RDT=ComplexF64,name="A")
    b=rand(ComplexF64,8,4);
    B=joMatrix(b;DDT=ComplexF32,RDT=ComplexF64,name="B")
    c=rand(ComplexF64,6,4);
    C=joMatrix(c;DDT=ComplexF32,RDT=ComplexF64,name="C")

define weights if needed

    w=rand(ComplexF64,3)

basic stack in function syntax

    S=joStack(A,B,C)

basic stack in [] syntax

    S=[A; B; C]

weighted stack

    S=joStack(A,B,C;weights=w)

"""
function joStack(ops::joAbstractLinearOperator...;
           weights::LocalVector{WDT}=zeros(0),name::String="joStack") where {WDT<:Number}

    isempty(ops) && throw(joStackException("empty argument list"))
    l=length(ops)
    for i=1:l
        ops[i].n==ops[1].n || throw(joStackException("size mismatch for $i operator"))
        deltype(ops[i])==deltype(ops[1]) || throw(joStackException("domain type mismatch for $i operator"))
        reltype(ops[i])==reltype(ops[1]) || throw(joStackException("range type mismatch for $i operator"))
    end
    (length(weights)==l || length(weights)==0) || throw(joStackException("lenght of weights vector does not match number of operators"))
    ws=Base.deepcopy(weights)
    ms=zeros(Int,l)
    ns=zeros(Int,l)
    mo=zeros(Int,l)
    no=zeros(Int,l)
    for i=1:l
        ms[i]=ops[i].m
        ns[i]=ops[1].n
    end
    for i=2:l
        mo[i]=mo[i-1]+ms[i-1]
    end
    m=sum(ms)
    n=ops[1].n
    weighted=(length(ws)==l)
    fops=Vector{joAbstractLinearOperator}(undef,0)
    fops_T=Vector{joAbstractLinearOperator}(undef,0)
    fops_A=Vector{joAbstractLinearOperator}(undef,0)
    fops_C=Vector{joAbstractLinearOperator}(undef,0)
    for i=1:l
        if weighted
            push!(fops,ws[i]*ops[i])
            push!(fops_T,ws[i]*transpose(ops[i]))
            push!(fops_A,conj(ws[i])*adjoint(ops[i]))
            push!(fops_C,conj(ws[i])*conj(ops[i]))
        else
            push!(fops,ops[i])
            push!(fops_T,transpose(ops[i]))
            push!(fops_A,adjoint(ops[i]))
            push!(fops_C,conj(ops[i]))
        end
    end
    return joCoreBlock{deltype(fops[1]),reltype(fops[1])}(name*"($l)",m,n,l,ms,ns,mo,no,ws,
                      fops,fops_T,fops_A,fops_C,@joNF,@joNF,@joNF,@joNF)
end
"""
    julia> op = joStack(l,op;[weights=...,][name=...])

Stack operator composed from l-times replicated square JOLI operator

# Signature

    joStack(l::Integer,op::joAbstractLinearOperator;
        weights::LocalVector{WDT}=zeros(0),name::String="joStack")
            where {WDT<:Number}

# Arguments

- `l`: # of replcated blocks
- `op`: JOLI operators (subtypes of joAbstractLinearOperator)
- keywords
    - `weights`: vector of waights for each operator
    - `name`: custom name

# Notes
- all operators must have the same # of columns (N)
- all given operators must have same domain/range types
- the domain/range types of joStack are equal to domain/range types of the given operators

# Example

define operator

    a=rand(ComplexF64,4,4);
    A=joMatrix(a;DDT=ComplexF32,RDT=ComplexF64,name="A")

define weights if needed

    w=rand(ComplexF64,3)

basic stack

    S=joStack(3,A)

weighted stack

    S=joStack(3,A;weights=w)

"""
function joStack(l::Integer,op::joAbstractLinearOperator;
           weights::LocalVector{WDT}=zeros(0),name::String="joStack") where {WDT<:Number}

    (length(weights)==l || length(weights)==0) || throw(joStackException("lenght of weights vector does not match number of operators"))
    ws=Base.deepcopy(weights)
    ms=zeros(Int,l)
    ns=zeros(Int,l)
    mo=zeros(Int,l)
    no=zeros(Int,l)
    for i=1:l
        ms[i]=op.m
        ns[i]=op.n
    end
    for i=2:l
        mo[i]=mo[i-1]+ms[i-1]
    end
    m=l*op.m
    n=op.n
    weighted=(length(weights)==l)
    fops=Vector{joAbstractLinearOperator}(undef,0)
    fops_T=Vector{joAbstractLinearOperator}(undef,0)
    fops_A=Vector{joAbstractLinearOperator}(undef,0)
    fops_C=Vector{joAbstractLinearOperator}(undef,0)
    for i=1:l
        if weighted
            push!(fops,ws[i]*op)
            push!(fops_T,ws[i]*transpose(op))
            push!(fops_A,conj(ws[i])*adjoint(op))
            push!(fops_C,conj(ws[i])*conj(op))
        else
            push!(fops,op)
            push!(fops_T,transpose(op))
            push!(fops_A,adjoint(op))
            push!(fops_C,conj(op))
        end
    end
    return joCoreBlock{deltype(op),reltype(op)}(name*"($l)",m,n,l,ms,ns,mo,no,ws,
                      fops,fops_T,fops_A,fops_C,@joNF,@joNF,@joNF,@joNF)
end


text:
"""
Package: Forecast

    co2(full = false)

Return dataset with atmospheric Carbon Dioxide Dry Air Mole Fractions from quasi-continuous measurements at Mauna Loa, Hawaii.

K.W. Thoning, A.M. Crotwell, and J.W. Mund (2020), Atmospheric Carbon Dioxide Dry Air Mole Fractions from continuous measurements at Mauna Loa, Hawaii, Barrow, Alaska, American Samoa and South Pole. 1973-2019, Version 2020-08 National Oceanic and Atmospheric Administration (NOAA), Global Monitoring Laboratory (GML), Boulder, Colorado, USA https://doi.org/10.15138/yaf1-bk21 FTP path: ftp://aftp.cmdl.noaa.gov/data/greenhouse_gases/co2/in-situ/surface/

# Arguments
- `full`: if `true` Returns the full original dataset from 1973 to 2020 in a DataFrame, otherwise returns the subset used in "STL: A Seasonal-Trend Decomposition Procedure Based on Loess" from Cleveland et. al. Its default value is `false`.

# Returns
Dataframe containing the descrived dataset.

# Examples
```julia-repl
julia> co2()
[ Info: Dataset used in Cleveland et al. paper
4612×2 DataFrame
  Row │ date        co2
      │ Date        Float64?
──────┼────────────────────────
    1 │ 1974-05-17      333.38
    2 │ 1974-05-18      333.11
    3 │ 1974-05-19      333.46
   [...]
```
"""
function co2(full::Bool = false)

    data = "data/co2.csv.gz"
    path = joinpath(splitdir(@__DIR__)[1], data)

    co2_df = GZip.open(path, "r") do io
        CSV.read(io,DataFrame)
    end

    if full
        @info "Full dataset from 1973 to 2020"
        return co2_df
    else
        @info "Dataset used in Cleveland et al. paper"
    end

    # recode missing values
    co2_df.value = replace(co2_df.value, -999.99 => missing)

    dates_co2 = Dates.Date(1973,1,1):Dates.Day(1):Dates.Date(2019,12,31)
    co2_sdf = DataFrame([dates_co2, co2_df.value],[:date,:co2])

    # paper dates
    # using 17 may instead 17 april (as the paper claims)
    # since there seems to be no data for april in 1974
    dates_co2_stl = Dates.Date(1974,5,17):Dates.Day(1):Dates.Date(1986,12,31)

    # revome leap year day
    dates_co2_stl = filter(dates_co2_stl) do x
        (Dates.isleapyear(x) & (Dates.month(x) == 2)) ? Dates.day(x) != 29 : true
    end

    @where(co2_sdf, in.(:date,  Ref(dates_co2_stl)))

end

"""
Package: Forecast

    seaborne(full = false)

Return estimates of world seaborne trade from AIS data collected by Marine Traffic.

By default a DataFrame containing deadweight imports for France, Germany and the United Kingdom from 2015-04-01 to 2021-05-02 is returned, otherwise a DataFrame is returned for the same countries with import and exports for the below fields:

num_pc:     number of port calls
mtc:        metric tons of cargo
dwt:        deadweight tonnage
suffix_ma:  30-day moving averages

Data available at UN COMTRADE Monitor.Cerdeiro, Komaromi, Liu and Saeed (2020).

# Returns
 Dataframe containing the seaborne dataset.

# Examples
```julia-repl
julia> seaborne()
[ Info: Seaborne deadweight trade imports from AIS
2199×4 DataFrame
  Row │ Date        France   Germany  UK
      │ Date        Int64    Int64    Int64
──────┼──────────────────────────────────────
    1 │ 2015-04-01   507946   878377  599573
    2 │ 2015-04-02   332043  1501614  772714
    3 │ 2015-04-03   810077   941663  262994
   [...]
```
"""
function seaborne(full::Bool = false)

    data = "data/seaborne.csv.gz"
    path = joinpath(splitdir(@__DIR__)[1], data)

    sb_df = GZip.open(path, "r") do io
        CSV.read(io, DataFrame, dateformat = "mm/dd/yyyy HH:MM:SS pp", types = Dict(:date => Date))
    end

    if full
        @info "Full dataset with estimates of world seaborne trade from AIS"
        return sb_df
    else
        @info "Seaborne deadweight trade imports from AIS"
    end

    # Group by Country and Flow to select Imports
    gbt = groupby(sb_df, [:country_name, :flow], sort=true)
    DataFrame(Dict(:Date => gbt[1].date ,
                   :UK => gbt[1].dwt, :Germany => gbt[3].dwt, :France => gbt[5].dwt))
end

"""
Package: Forecast

    quakes()

Return the number of earthquakes per year on earth with a magnitude higher or equal to six from 1950 to 2020. The data has been collected from https://earthquake.usgs.gov/ and aggregated.

# Examples
```julia-repl
julia> quakes()
71×2 DataFrame
 Row │ year        quakes
     │ Date        Int64
─────┼────────────────────
   1 │ 1950-01-01     138
   2 │ 1951-01-01     151
   3 │ 1952-01-01     181
   [...]
```
"""
function quakes()

    data = "data/quakes.csv.gz"
    path = joinpath(splitdir(@__DIR__)[1], data)

    qk_df = GZip.open(path, "r") do io
        CSV.read(io,DataFrame, dateformat = "yyyy", types = Dict(:year => Date))
    end

    qk_df

end

"""
Package: Forecast

    air()

Return the classic Box & Jenkins airline data. Monthly totals of international airline passengers from 1949 to 1960.

Box, G. E. P., Jenkins, G. M. and Reinsel, G. C. (1976) _Time Series Analysis, Forecasting and Control._ Third Edition. Holden-Day. Series G.

# Returns
Dataframe containing the descrived dataset.

# Examples
```julia-repl
julia> air()
71×2 DataFrame
 Row │ year        quakes
     │ Date        Int64
─────┼────────────────────
   1 │ 1950-01-01     138
   2 │ 1951-01-01     151
   3 │ 1952-01-01     181
   [...]
```
"""
function air()

    data = "data/air.csv.gz"
    path = joinpath(splitdir(@__DIR__)[1], data)

    df = GZip.open(path, "r") do io
        CSV.read(io,DataFrame)
    end

    t =  DataFrame([Date.(df.year,df.month)],[:Date])
    x = DataFrame([df.passengers],[:Passengers])

    [t x]
end

"""
Package: Forecast

    london()

Return ten years of monthly data about weather and crime in Greater London from 2008 to 2018.

Data has been collected and joined from london.gov.uk and metoffice.gov.uk (Heathrow Station).

# Weather Variables

- `MaxTemp`:  Mean daily maximum temperature in C°
- `MinTemp`:  Mean daily minimum temperature in C°
- `AirFrost`: Days of air frost
- `Rain`:     Total rainfall in mm
- `Sun`:      Total sunshine durationin hours

# Crime Variables and its aggregated categories
```
┌─────────────────────────────┬────────────────────────────────────────┐
│                       Crime │                               Category │
├─────────────────────────────┼────────────────────────────────────────┤
│                    Burglary │            Burglary in Other Buildings │
│                             │                 Burglary in a Dwelling │
│                      Damage │            Criminal Damage To Dwelling │
│                             │       Criminal Damage To Motor Vehicle │
│                             │      Criminal Damage To Other Building │
│                             │                  Other Criminal Damage │
│                       Drugs │                       Drug Trafficking │
│                             │                            Other Drugs │
│                             │                    Possession Of Drugs │
│                       Fraud │                     Counted per Victim │
│                             │                  Other Fraud & Forgery │
│                       Other │                         Going Equipped │
│                             │                       Other Notifiable │
│                     Robbery │                      Business Property │
│                             │                      Personal Property │
│                      Sexual │                           Other Sexual │
│                             │                                   Rape │
│                       Theft │                  Handling Stolen Goods │
│                             │ Motor Vehicle Interference & Tampering │
│                             │                            Other Theft │
│                             │                     Other Theft Person │
│                             │               Theft From Motor Vehicle │
│                             │                       Theft From Shops │
│                             │          Theft/Taking Of Motor Vehicle │
│                             │            Theft/Taking of Pedal Cycle │
│                    Violence │                    Assault with Injury │
│                             │                         Common Assault │
│                             │                   Grievous Bodily Harm │
│                             │                             Harassment │
│                             │                                 Murder │
│                             │                       Offensive Weapon │
│                             │                         Other violence │
│                             │                           Wounding/GBH │
└─────────────────────────────┴────────────────────────────────────────┘
```

# Returns
Dataframe containing the descrived dataset.

# Examples
```julia-repl
julia> london()
132×15 DataFrame
 Row │ Date        MaxTemp  MinTemp  AirFrost
     │ Date        Float64  Float64  Int64
─────┼─────────────────────────────────────────[...]
   1 │ 2008-01-01     10.4      4.7         0
   2 │ 2008-02-01     11.0      2.0         7
   3 │ 2008-03-01     10.6      3.7         2
   [...]
```
"""
function london()

    data = "data/london.csv.gz"
    path = joinpath(splitdir(@__DIR__)[1], data)

    df = GZip.open(path, "r") do io
        CSV.read(io,DataFrame)
    end

    df = GZip.open(path, "r") do io
        CSV.read(io, DataFrame, dateformat = "yyyy-mm-dd", types = Dict(:Date => Date))
    end

    df

end


text:
@testset "inverse.jl" begin
    InverseFunctions.test_inverse(log1pexp, randn())
    InverseFunctions.test_inverse(logexpm1, randexp())

    InverseFunctions.test_inverse(log1mexp, -randexp())

    InverseFunctions.test_inverse(log2mexp, log(2) - randexp())

    InverseFunctions.test_inverse(logistic, randn())
    InverseFunctions.test_inverse(logit, rand())
end


text:
using DataFrames: DataFrames

@testset "stats" begin
    idata = load_arviz_data("centered_eight")

    @testset "compare" begin
        idata2 = load_arviz_data("non_centered_eight")
        model_dict = Dict("a" => idata, "b" => idata2)
        loo_dict = Dict("a" => loo(idata), "b" => loo(idata2))
        df = compare(model_dict)
        @test df isa DataFrames.DataFrame
        df2 = compare(loo_dict)
        @test df2 isa DataFrames.DataFrame
        @test df == df2
    end

    @testset "hdi" begin
        rng = Random.MersenneTwister(42)
        x = randn(rng, 100)
        @test hdi(x) == ArviZ.arviz.hdi(x)
    end

    @testset "r2_score" begin
        rng = Random.MersenneTwister(42)
        ytrue = randn(rng, 100)
        ypred = randn(rng, 100)
        df = r2_score(ytrue, ypred)
        @test df isa DataFrames.DataFrame
        @test all(df == ArviZ.todataframes(ArviZ.arviz.r2_score(ytrue, ypred)))
    end

    @testset "loo" begin
        df = loo(idata; pointwise=false)
        @test df isa DataFrames.DataFrame
        @test all(df == ArviZ.todataframes(ArviZ.arviz.loo(idata; pointwise=false)))
    end

    @testset "psislw" begin
        @testset for sz in ((1000,), (10, 1000))
            logw = randn(sz)
            logw_smoothed, k = psislw(copy(logw), 0.9)

            # check against PSIS.jl
            result = psis(copy(logw), 0.9)
            @test exp.(logw_smoothed) ≈ result.weights
            @test k ≈ result.pareto_shape

            # check against Python ArviZ
            logw_smoothed2, k2 = ArviZ.arviz.psislw(copy(logw), 0.9)
            # NOTE: sometimes the smoothed weights disagree, while the shapes agree
            @test_skip logw_smoothed ≈ logw_smoothed2
            if length(sz) == 1
                @test k ≈ k2[] # k2 is a 0-dimensional array
            else
                @test k ≈ k2
            end
        end
    end

    @testset "waic" begin
        df = waic(idata; pointwise=false)
        @test df isa DataFrames.DataFrame
        @test all(df == ArviZ.todataframes(ArviZ.arviz.waic(idata; pointwise=false)))
    end

    @testset "loo_pit" begin
        ret = loo_pit(idata; y="obs")
        @test ret == ArviZ.arviz.loo_pit(idata; y="obs")
    end

    @testset "summarystats" begin
        rng = MersenneTwister(42)
        nchains, ndraws = 4, 10
        idata = convert_to_inference_data(
            Dict(
                "a" => randn(rng, nchains, ndraws), "b" => randn(rng, nchains, ndraws, 3, 4)
            ),
        )

        s = summarystats(idata)
        @test s isa DataFrames.DataFrame
        @test string(first(names(summarystats(idata)))) == "variable"
        @test string(first(names(summarystats(idata; fmt="wide")))) == "variable"
        @test :variable in propertynames(summarystats(idata; fmt="wide"))
        @test "a" ∈ s.variable
        @test "b" ∉ s.variable
        @test "b[0, 0]" ∉ s.variable
        @test "b[1, 1]" ∈ s.variable

        s2 = summarystats(idata; fmt="long")
        @test s2 isa DataFrames.DataFrame
        @test string(first(names(s2))) == "statistic"
        @test "mean" ∈ s2.statistic

        s3 = summarystats(idata; fmt="xarray")
        @test s3 isa ArviZ.Dataset
    end

    @testset "ArviZ.summary" begin
        rng = MersenneTwister(42)
        nchains, ndraws = 4, 10
        data = Dict(
            "a" => randn(rng, nchains, ndraws), "b" => randn(rng, nchains, ndraws, 3, 4)
        )

        @test ArviZ.summary(data) isa DataFrames.DataFrame
    end
end