text:
"""
Load and plot simulations results on confidence bands' coverages
"""

#region import and models
using Pkg
Pkg.activate(".")
Pkg.instantiate()
using DrWatson
using DataFrames
using PyPlot
pygui(true)
using JLD2
using SharedArrays
using Statistics
using ScoreDrivenERGM
using TableView
using Blink
viewtab(df) = body!(Window(), showtable(df))
#endregion

# #region load and plot coverage simulations

@elapsed df = collect_results( datadir("sims", "dgp&FIl_conf"))
df["modelTag"] = string.(df["model"])

df = df[getfield.(df.model, :scoreScalingType) .== "FISH_D",:]
df = df[getfield.(df.dgp, :type) .== "AR",:]
viewtab(df[[:modelTag, :N, :T]])

df.confBndCoverFiltPar = [zeros(1, 50) for i = 1: nrow(df)]
df.confBndCoverPar = [zeros(1, 50) for i = 1: nrow(df)]
df.avgConfBndCoverFiltPar = zeros(nrow(df))
df.avgConfBndCoverPar = zeros(nrow(df))

obsShift = 1
for res in eachrow(df)
    T = res.T
    limitSample = sum(res.allConfBandsFiltPar[1,1,1,1,:] .!= 0 )
    nBands = size(res.allConfBandsFiltPar)[3]
    avgCoverFiltPar = zeros(nBands, limitSample)
    avgCoverPar = zeros(nBands, limitSample)

    constInds = any(res.allvEstSdResPar[3:3:6, 1:limitSample] .< 0.00005, dims=1)

    # Threads.@threads
    for n = 1:limitSample
        coverFiltParUnc = DynNets.conf_bands_coverage(res.allParDgpT[:,1:T-obsShift,n],   res.allConfBandsFiltPar[:,1+obsShift:end,:,:,n])

        coverParUnc = DynNets.conf_bands_coverage(res.allParDgpT[:,1:T-obsShift,n],   res.allConfBandsPar[:,1+obsShift:end,:,:,n])

        avgCoverFiltPar[1, n] =  mean(coverFiltParUnc[:,:,1])
        avgCoverPar[1, n] =  mean(coverFiltParUnc[:,:,1])

    end

    res.avgConfBndCoverFiltPar = mean(avgCoverFiltPar[.!constInds])
    res.avgConfBndCoverPar = mean(avgCoverPar[.!constInds])

    res.confBndCoverFiltPar = avgCoverFiltPar
    res.confBndCoverPar = avgCoverPar
end

viewtab(df[df.N.==100,:][[:modelTag, :m,  :dgp, :N, :T, :avgConfBndCoverFiltPar, :confBndCoverFiltPar]])

#region obtain box plots to compare subsets of the simulations run
begin

parUncMethod = "WHITE-MLE"
parUncMethod = "NPB-MVN" #
limitSample =50
indB = 1
tVals = [300]
nVals = [100  ]
model = DynNets.SdErgmDirBin0Rec0_mle(scoreScalingType="FISH_D")
model = DynNets.SdErgmDirBin0Rec0_pmle(scoreScalingType="FISH_D")
dgpSetting = DynNets.list_example_dgp_settings(DynNets.SdErgmDirBin0Rec0_mle()).dgpSetARlowlow
dgpSetting.opt.B[1] = 1

modelTags = [DynNets.name(model)] # ["SdErgmDirBin0Rec0_mle(Bool[1, 1], scal = HESS_D)"]# unique(df["modelTag"])
# modelTags =unique(df["modelTag"])
nNVals = length(nVals)
nTVals = length(tVals)
nModels = length(modelTags)
nSample = 50 #length(unique(df["nSample"]))>1 ? missing : (unique(df["nSample"])[1] )
nErgmPar = 2

indQuant = 1
nBands = 2
allMeanDgpStd =zeros(nNVals, nTVals, nModels)
allMeanDgpStdConf =zeros(nNVals, nTVals, nModels)
allAvgCover =zeros(nErgmPar ,nNVals, nTVals, nModels, nBands, limitSample)
allConstInds = falses(nErgmPar, nNVals, nTVals, nModels, limitSample)
allErrInds = trues(nErgmPar, nNVals, nTVals, nModels, limitSample)

avgCover = SharedArray(zeros(nErgmPar, nBands, limitSample))
avgCoverJoint = SharedArray(zeros(nErgmPar, nBands, limitSample))
indM=1
obsShift = 1

for (indT, T) in Iterators.enumerate(tVals)
    for (indN, N) in Iterators.enumerate(nVals)
        for (indM, modelTag) in Iterators.enumerate(modelTags)

            res = filter([:modelTag, :T, :N, :dgp, :S, :m] => (mtag,t,n, d, s, m) -> all((mtag==modelTag, t==T, n==N, d == dgpSetting, s == nSample, m==parUncMethod)), df)

            nrow(res) != 1 ? error("$N, $T,  $(size(res))") : res = res[1,:]

            # Threads.@threads
            for n = 1:limitSample
                coverFiltParUnc = DynNets.conf_bands_coverage(res.allParDgpT[:,1:T-obsShift,n],   res.allConfBandsFiltPar[:,1+obsShift:end,:,:,n])

                coverParUnc = DynNets.conf_bands_coverage(res.allParDgpT[:,1:T-obsShift,n],   res.allConfBandsPar[:,1+obsShift:end,:,:,n])

                for indPar = 1:2
                    avgCover[indPar, 1, n] =  mean(coverFiltParUnc[indPar,:,1])

                    avgCover[indPar, 2, n] =  mean(coverParUnc[indPar,:,1])
                end

            end

            variability = res.allfVecT_filt |>  x -> dropdims(std(x, dims=2), dims=2) |> x-> replace(e -> (isnan(e) ? 0 : e), x ) |> x ->((mean = mean(x) , conf =  1.96.*std(x) ))

            allMeanDgpStd[indN, indT, indM] = round.(variability.mean, sigdigits=2)
            allMeanDgpStdConf[indN, indT, indM] = round.(variability.conf, sigdigits=2)
            allAvgCover[:, indN, indT, indM, :, 1:limitSample] = avgCover
            allErrInds[:, indN, indT, indM, 1:limitSample] = res.errInds[:, 1:limitSample]
            allConstInds[:, indN, indT, indM, 1:limitSample] .= any(res.allvEstSdResPar[3:3:6, 1:limitSample] .< 0.00005, dims=1)

        end
    end
end

nominalLevel = 0.95
parNames = ["θ", "η", "mean θ η"]
BandNames = ["Parameters + Filtering Uncertainty", "Parameters Uncertainty"]

fig, ax1 = plt.subplots(3, length(tVals),figsize=(12, 6), sharey =true)
fig.canvas.set_window_title("Confidence Bands' Coverages $(BandNames[indB])")
fig.subplots_adjust(left=0.075, right=0.95, top=0.9, bottom=0.25)
fig.suptitle("Confidence Bands' Coverages $(BandNames[indB]) DGP = $(dgpSetting.type),\n filter = $(modelTags[indM]), Cov-Estimate : $(parUncMethod)")

for (indT, T) in Iterators.enumerate(tVals)
    for indPar = 1:3

        if indPar == 3
            data = [c[(.!indC).&(.!indE)] for (c, indC, indE) in Iterators.zip(eachrow(dropdims(mean(allAvgCover[:,:,indT,indM,indB,:], dims=1), dims=1)), eachrow(allConstInds[1, :, indT, indM, :]), eachrow(allErrInds[1, :, indT, indM, :]))]
        else
            data = [c[(.!indC).&(.!indE)] for (c, indC, indE) in Iterators.zip(eachrow(allAvgCover[indPar,:,indT,indM,indB,:]), eachrow(allConstInds[indPar, :, indT, indM, :]), eachrow(allErrInds[indPar, :, indT, indM, :]))]
        end
        # data = [c for (c, ind) in Iterators.zip(eachrow(allAvgCover[indPar,:,indT,indM,indB,:]), eachrow(allConstInds[indPar, :, indT, indM, :]))]

        bp = ax1[indPar, indT].boxplot(data, notch=0, sym="+", vert=1, whis=1.5, showfliers =true, showmeans=true)

        ax1[indPar, indT].yaxis.grid(true, linestyle="-", which="major", color="lightgrey", alpha=0.5)

        # Hide these grid behind plot objects
        xlims = ax1[indPar, indT].get_xlim()
        ax1[indPar, indT].hlines(nominalLevel, xlims[1], xlims[2], linestyle=":" , colors = "r")
        ax1[indPar, indT].set_ylim([0.70, 1])
        ax1[indPar, indT].set_axisbelow(true)
        ax1[indPar, indT].set_title("T = $T")
        # ax1[indPar, indT].set_xlabel("Network Size")
        ax1[indPar, indT].set_ylabel("$(parNames[indPar])")
        if indPar ==3
            xlab = ["N = $n \n <std($(parNames[indPar]) )> =  $(allMeanDgpStd[indN, indT, indM])  " for (indN,n) in Iterators.enumerate(nVals) ] #\n ($(allMeanDgpStdConf[indPar, indN, indT, indM])
            # ax1[indPar, indT].set_xticklabels(xlab, rotation=0, fontsize=8)
        end
    end
end
tight_layout()

end
#endregion

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
"""
    pop_index(x)

Returns the index `x` that corresponds to `pop`. New subtypes of `AbstractIndex`
must implement a unique `pop_index` method.
"""
pop_index(x::SimpleIndex{name}) where {name} = SimpleIndex{name}(pop(values(x)))
pop_index(x::Index{name}) where {name} = Index{name}(pop(keys(x)), pop(values(x)))

"""
    popfirst_index(x)

Returns the index `x` that corresponds to `pop`. New subtypes of `AbstractIndex`
must implement a unique `popfirst_index` method.
"""
popfirst_index(x::SimpleIndex) = pop_index(x)
popfirst_index(x::Index{name}) where {name} = Index{name}(popfirst(keys(x)), pop(values(x)))

"""
    pop_index!(x)

Returns the index `x` that corresponds to `pop!`.
"""
function pop_index!(x::SimpleIndex)
    pop!(values(x))
    return x
end
function pop_index!(x::AbstractIndex)
    pop!(keys(x))
    pop!(values(x))
    return x
end

"""
    popfirst_index!(x)

Returns the index `x` that corresponds to `popfirst!`.
"""
popfirst_index!(x::SimpleIndex) = pop_index!(x)
function popfirst_index!(x::Index)
    popfirst!(keys(x))
    pop!(values(x))
    return x
end

"""
    push_index!(x, items...)

Returns the index `x` appropriate for `push(vector, items...)`. If `items` are
`Pairs` then the first element of each pair is treated as a key to be appended to
those of of `x`.
"""
function push_index!(x::SimpleIndex, items...)
    set_length!(values(x), length(x) + length(items))
    return x
end
function push_index!(x::AbstractIndex, items...)
    len = length(x) + length(items)
    set_length!(keys(x), len)
    set_length!(values(x), len)
    return x
end
function push_index!(x::AbstractIndex, items::Pair...)
    len = length(x) + length(items)
    set_length!(values(x), len)
    ks = keys(x)
    f = unique_offset(ks)
    for i in items
        k, v = i
        if k in ks
            push!(ks, add_offset(ks, k, f))
        else
            push!(ks, k)
        end
    end
    return x
end

"""
    pushfirst_index!(x, items...)
"""
pushfirst_index!(x::SimpleIndex, items...) = push_index!(x, items...)
pushfirst_index!(x::AbstractIndex, items...) = push_index!(x, items...)
function pushfirst_index!(x::AbstractIndex, items::Pair...)
    len = length(x) + length(items)
    set_length!(values(x), len)
    ks = keys(x)
    f = unique_offset(ks)
    for i in items
        k, v = i
        if k in ks
            pushfirst!(ks, add_offset(ks, k, f))
        else
            pushfirst!(ks, k)
        end
    end
    return x
end

ext: jl


text:
# Define core types and core functions
# Base functions are provided as core featuers of OpenGene, so it is not inside a sub-module
# ===========

include("Types/Types.jl")


text:
function GeometricMean()
    return GeometricMean(())
end

function GeometricMean(arg0::GeometricMean)
    return GeometricMean((GeometricMean,), arg0)
end

function GeometricMean(arg0::SumOfLogs)
    return GeometricMean((SumOfLogs,), arg0)
end

function aggregate(obj::GeometricMean, arg0::GeometricMean)
    return jcall(obj, "aggregate", void, (GeometricMean,), arg0)
end

function clear(obj::GeometricMean)
    return jcall(obj, "clear", void, ())
end

function copy(obj::GeometricMean)
    return jcall(obj, "copy", GeometricMean, ())
end

function evaluate(obj::GeometricMean, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function get_n(obj::GeometricMean)
    return jcall(obj, "getN", jlong, ())
end

function get_result(obj::GeometricMean)
    return jcall(obj, "getResult", jdouble, ())
end

function increment(obj::GeometricMean, arg0::jdouble)
    return jcall(obj, "increment", void, (jdouble,), arg0)
end


text:
function graft!(a, n, clr::Int, colors::Vararg{Int})
	graft!(a, n, clr)
	for c in colors
		graft!(a, n, c)
	end

	return nothing
end
function graft!(a, n, colors::Vector{Bool})
	for (i,c) in enumerate(colors)
		if c
			graft!(a, n, i)
		end
	end

	return nothing
end
function graft!(a::ARGNode, n::ARGNode, c::Int)
	# Check that `n` does not have an ancestor
	@assert !hasancestor(n, c)
	# Proceed: set `a` as ancestor of `n`
	set_ancestor!(n, a, c)
	# If `n` is not already in children of `a`, put it there
	if !haschild(a, n)
		add_child!(a, n)
	end

	return nothing
end

function graft!(a::RootNode, n::ARGNode, c::Int)
	@assert !hasancestor(n, c)
	set_ancestor!(n, a, c)

	return nothing
end

function set_root!(n::ARGNode, c::Int)
	graft!(RootNode(), n, c)
end

function unset_child!(a::ARGNode, n, clr)
	@assert ancestor(n, clr) == a
	@assert haschild(a, n)
	# Is a the ancestor of n for another color?
	# If yes, we're not doing anything
	for (i,c) in color(n)
		if c && (i != clr)
			return nothing
		end
	end
	# Else, remove n from the children of a
	deleteat!(findfirst(==(n), a.children))
end
unset_child!(a, n, clr) = nothing

function prune!(n::ARGNode, c::Int)
	a = ancestor(n, c)
	unset_child!(a, n, c)
	set_ancestor!(n, nothing, c)
end

# function set_global_root!(arg::ARG)
# 	gr = RootNode([])
# 	for (clr, r) in arg.roots
# 		@assert isroot(r, clr)
# 		prune!(r, clr)
# 		graft!(gr, r, clr)
# 	end
# end

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
using MLStyle

function refblock(expr, refs::Dict{Symbol, Symbol}, scopesafe :: Bool)
    rec(x) = refblock(x, refs, true)
    @match expr begin
        a && if scopesafe end && (
                # 函数定义，开启新的scope
                Expr(:function, _...)   ||
                Expr(:->, _...)         ||
                Expr(:(=), Expr(:call, _...), _...)
        ) => refblock(a, Dict{Symbol, Symbol}(), false)

        :($a = &$b) =>
            begin
                refs[a] = b
                nothing
            end

        a::Symbol && if haskey(refs, a) end => refs[a]

        Expr(hd, tl...) => Expr(hd, map(rec, tl)...)

        a => a
    end
end
macro refblock(expr)
    refblock(expr, Dict{Symbol, Symbol}(), true) |> esc
end

function myf(x)
    @refblock begin
        a = &x
        a = a + 1
        println(x)
    end
end
myf(1)

max_stars_repo_name: thautwarm/MLStyle-Playground

max_issues_repo_name: thautwarm/MLStyle-Playground

max_forks_repo_name: thautwarm/MLStyle-Playground


text:
export gauss_quadrature

# Dispatch on Val, so the return Vec length can be inferred.
# Should implement a non-value based version which returns a Vector at some point.
include("legendre_line.jl")      # Accurate to degree 2p-1
include("legendre_hypercube.jl") # Accurate to degree 2p-1
include("legendre_triangle.jl")  # Accurate to degree p

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# Diffusion maps
# --------------
# Diffusion maps,
# Coifman, R. & Lafon, S., Applied and Computational Harmonic Analysis, Elsevier, 2006, 21, 5-30

"""

    DiffMap{T <: Real} <: AbstractDimensionalityReduction

The `DiffMap` type represents diffusion maps model constructed for `T` type data.
"""
struct DiffMap{T <: Real} <: AbstractDimensionalityReduction
    t::Int
    α::Real
    ɛ::Real
    λ::AbstractVector{T}
    K::AbstractMatrix{T}
    proj::Projection{T}
end

## properties
outdim(R::DiffMap) = size(R.proj, 1)
eigvals(R::DiffMap) = R.λ

## custom
"""Returns the kernel matrix of the diffusion maps model `R`"""
kernel(R::DiffMap) = R.K

## show
summary(io::IO, R::DiffMap) = print(io, "Diffusion Maps(outdim = $(outdim(R)), t = $(R.t), α = $(R.α), ɛ = $(R.ɛ))")
function show(io::IO, R::DiffMap)
    summary(io, R)
    io = IOContext(io, :limit=>true)
    println(io)
    println(io, "Kernel: ")
    Base.print_matrix(io, R.K, "[", ",","]")
    println(io)
    println(io, "Embedding:")
    Base.print_matrix(io, transform(R), "[", ",","]")
end

## interface functions
"""
    fit(DiffMap, data; maxoutdim=2, t=1, α=1.0, ɛ=1.0)

Fit a isometric mapping model to `data`.

# Arguments
* `data`: a matrix of observations. Each column of `data` is an observation.

# Keyword arguments
* `maxoutdim`: a dimension of the reduced space.
* `t`: a number of transitions
* `α`: a normalization parameter
* `ɛ`: a Gaussian kernel variance (the scale parameter)

# Examples
```julia
M = fit(DiffMap, rand(3,100)) # construct diffusion map model
R = transform(M)              # perform dimensionality reduction
```
"""
function fit(::Type{DiffMap}, X::AbstractMatrix{T}; maxoutdim::Int=2, t::Int=1, α::Real=0.0, ɛ::Real=1.0) where {T<:Real}
    # compute kernel matrix
    sumX = sum(X.^ 2, dims=1)
    L = exp.(-( transpose(sumX) .+ sumX .- 2*transpose(X) * X ) ./ convert(T, ɛ))
    # L = pairwise((x,y)-> exp(-norm(x-y,2)^2/ε), Xtr)

    # Calculate Laplacian & normalize it
    if α > 0
        D = transpose(sum(L, dims=1))
        L ./= (D * transpose(D)) .^ convert(T, α)
    end
    D = Diagonal(vec(sum(L, dims=1)))
    M = inv(D)*L

    # D = Diagonal(vec(sum(L, dims=1)))
    # D⁻ᵅ = inv(D^α)
    # Lᵅ = D⁻ᵅ*L*D⁻ᵅ
    # Dᵅ = Diagonal(vec(sum(Lᵅ, dims=1)))
    # M = inv(Dᵅ)*Lᵅ

    # Eigendecomposition & reduction
    F = eigen(M, permute=false, scale=false)
    λ = real.(F.values)
    idx = sortperm(λ, rev=true)[2:maxoutdim+1]
    λ = λ[idx]
    V = real.(F.vectors[:,idx])
    Y = (λ.^t) .* V'

    return DiffMap{T}(t, α, ɛ, λ, L, Y)
end

"""
    transform(R::DiffMap)

Transforms the data fitted to the diffusion map model `R` into a reduced space representation.
"""
transform(R::DiffMap) = R.proj


text:
using HetaSimulator, Plots

################################## Model Upload ###########################################
p = load_platform("./boehm");
m = p.models[:nameless]

################################## Simple test simulation ######################################

sim(m, tspan = (0, 600)) |> plot

################################## Loading conditions and measurements ##########################

cond_df = read_conditions("./boehm/conditions.csv");
measurement_df = read_measurements("./boehm/measurements.csv");

add_conditions!(p, cond_df)
add_measurements!(p, measurement_df)

sim(p) |> plot

################################ Fitting #######################

fit_cons = [
  :Epo_degradation_BaF3=>0.026982514033029,
  :k_exp_hetero=>0.0000100067973851508,
  :k_exp_homo=>0.006170228086381,
  :k_imp_hetero=>0.0163679184468,
  :k_imp_homo=>97749.3794024716,
  :k_phos=>15766.5070195731,
  :sd_pSTAT5A_rel=>3.85261197844677,
  :sd_pSTAT5B_rel=>6.59147818673419,
  :sd_rSTAT5A_rel=>3.15271275648527
]

fit1 = fit(p, fit_cons; alg = CVODE_BDF()) # default algorithm doesn't work

################################## plot fitted #################

sim(p) |> plot # default
sim(p, parameters_upd = fit1.optim) |> plot # best fit
sim(p, parameters_upd = [:k_phos=>1000.]) |> plot # modified plot

max_stars_repo_name: hetalang/qsp-benchmark

max_issues_repo_name: hetalang/qsp-benchmark

max_forks_repo_name: hetalang/qsp-benchmark


text:
"""
    black_body(T, K₀, σ)
    black_body(T)

Thermal infrared, *i.e.* longwave radiation emitted from a black body at temperature T.

- `T`: temperature of the object in Celsius degree
- `K₀`: absolute zero (°C)
- `σ` (``W\\ m^{-2}\\ K^{-4}``) [Stefan-Boltzmann constant](https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_law)

# Note

`K₀` and `σ` are taken from [`Constants`](@ref) if not provided.

"""
function black_body(T, K₀, σ)
  Tₖ = T - K₀
  σ * (Tₖ^4.0)
end

function black_body(T)
    constants = Constants()
    black_body(T, constants.K₀, constants.σ)
end

"""
Thermal infrared, *i.e.* longwave radiation emitted from an object at temperature T.

- `T`: temperature of the object in Celsius degree
- `ε` object [emissivity](https://en.wikipedia.org/wiki/Emissivity) (not to confuse with ε the
ratio of molecular weights from [`Constants`](@ref)). A typical value for a leaf is 0.955.
- `K₀`: absolute zero (°C)
- `σ` (``W\\ m^{-2}\\ K^{-4}``) [Stefan-Boltzmann constant](https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_law)

# Note

`K₀` and `σ` are taken from [`Constants`](@ref) if not provided.

# Examples

```julia
# Thermal infrared radiation of water at 25 °C:
grey_body(25.0, 0.96)
```
"""
function grey_body(T, ε, K₀, σ)
  ε * black_body(T, K₀, σ)
end

function grey_body(T, ε)
    constants = Constants()
    grey_body(T, ε, constants.K₀, constants.σ)
end

"""
    net_longwave_radiation(T₁,T₂,ε₁,ε₂,F₁,K₀,σ)
    net_longwave_radiation(T₁,T₂,ε₁,ε₂,F₁)

Net longwave radiation fluxes (*i.e.* thermal radiation, W m-2) between an object and another.
The object of interest is at temperature T₁ and has an emissivity ε₁, and the object with
wich it exchanges energy is at temperature T₂ and has an emissivity ε₂.

If the result is positive, then the object of interest gain energy.

# Arguments

- `T₁` (Celsius degree): temperature of the target object (object 1)
- `T₂` (Celsius degree): temperature of the object with which there is potential exchange (object 2)
- `ε₁`: object 1 emissivity
- `ε₂`: object 2 emissivity
- `F₁`: view factor (0-1), *i.e.* visible fraction of object 2 from object 1 (see note)
- `K₀`: absolute zero (°C)
- `σ` (``W\\ m^{-2}\\ K^{-4}``) [Stefan-Boltzmann constant](https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_law)

# Note

`F₁`, the view factor (also called shape factor) is a coefficient applied to the semi-hemisphere
field of view of object 1 that "sees" object 2. E.g. a leaf can be viewed as a plane. If one side
of the leaf sees only object 2 in its field of view (e.g. the sky), then `F₁ = 1`.
Then the net longwave radiation flux for this part of the leaf is multiplied by its actual
surface to get the exchange. Note that we apply reciprocity between the two objects for
the view factor (they have the same value), *i.e.*: A₁F₁₂ = A₂F₂₁.

Then, if we take a leaf as object 1, and the sky as object 2, the visible fraction of
sky viewed by the leaf would be:

- `0.5` if the leaf is on top of the canopy, *i.e.* the upper side of the leaf sees the sky,
the side bellow sees other leaves and the soil.
- between 0 and 0.5 if it is within the canopy and partly shaded by other objects.

Note that `A₁` for a leaf is twice its common used leaf area, because `A₁` is the **total**
leaf area of the object that exchange energy.

```julia
# Net thermal radiation fluxes between a leaf and the sky considering the leaf at the top of
# the canopy:
Tₗ = 25.0 ; Tₐ = 20.0
ε₁ = 0.955 ; ε₂ = 1.0
Rₗₗ = net_longwave_radiation(Tₗ,Tₐ,ε₁,ε₂,1.0)
Rₗₗ

# Rₗₗ is the net longwave radiation flux between the leaf and the atmosphere per surface area.
# To get the actual net longwave radiation flux we need to multiply by the surface of the
# leaf, e.g. for a leaf of 2cm²:
leaf_area = 2e-4 # in m²
Rₗₗ * leaf_area

# The leaf lose ~0.0055 W towards the atmosphere.
```

# References

Cengel, Y, et Transfer Mass Heat. 2003. A practical approach. New York, NY, USA: McGraw-Hill.
"""
function net_longwave_radiation(T₁,T₂,ε₁,ε₂,F₁,K₀,σ)
    (black_body(T₂,K₀,σ) - black_body(T₁,K₀,σ)) / (1.0 / ε₁ + 1.0 / ε₂ - 1.0) * F₁
end

function net_longwave_radiation(T₁,T₂,ε₁,ε₂,F₁)
    constants = Constants()
    net_longwave_radiation(T₁,T₂,ε₁,ε₂,F₁,constants.K₀,constants.σ)
end

"""
    atmosphere_emissivity(Tₐ,eₐ)

Emissivity of the atmoshpere at a given temperature and vapor pressure.

# Arguments

- `Tₐ` (°C): air temperature
- `eₐ` (kPa): air vapor pressure
- `K₀` (°C): absolute zero

# Examples

```julia
Tₐ = 20.0
VPD = 1.5
atmosphere_emissivity(Tₐ, vapor_pressure(Tₐ,VPD))
```

# References

Leuning, R., F. M. Kelliher, DGG de Pury, et E.-D. SCHULZE. 1995. Leaf nitrogen,
photosynthesis, conductance and transpiration: scaling from leaves to canopies ». Plant,
Cell & Environment 18 (10): 1183‑1200.
"""
function atmosphere_emissivity(Tₐ,eₐ,K₀)
    0.642 * (eₐ * 100 / (Tₐ - K₀))^(1 / 7)
end

function atmosphere_emissivity(Tₐ,eₐ)
    atmosphere_emissivity(Tₐ,eₐ,Constants().K₀)
end


text:
function FormationsTrack(Depth,WellUnit;DepthFrom=minimum(Depth),DepthTo=maximum(Depth),
                        FmC1=RGB(0.0,0.0,1.0), FmC2=RGB(1.0,0.0,0.0), FmAlpha=0.7)

    ax=gca()                       #If Units are Provided
    ntop=size(WellUnit,1)
    if ntop==1
            zcolor=range(FmC1, stop=FmC2, length=2)
        else
            zcolor=range(FmC1, stop=FmC2, length=ntop)
    end
    zcols = map(col -> (red(col), green(col), blue(col)), zcolor)

    for i=1:ntop
        ax.axhspan(WellUnit[i,1],WellUnit[i,2],xmin=0,xmax=1,alpha=FmAlpha,color=zcols[i])
        ax.hlines([WellUnit[i,1],WellUnit[i,2]],0,1, linewidth=0.5)
    end
    setp(ax.get_yticklabels(),visible=false)
    setp(ax.get_xticklabels(),visible=false)
    ax.set_ylim([DepthTo,DepthFrom])
end


text:
"""
Tests for R2BP dynamics.
"""
module R2BPTests

    using ModelingToolkit, AstrodynamicalModels, Test

    @testset "R2BP Model Constructors" begin
        model = R2BP(; stm=false, structural_simplify=false)
        model = R2BP(; stm=true, structural_simplify=true)
        @test true
    end

    @testset "R2BP Model Calculations" begin
        vectorfield = R2BPFunction()

        r = [-11e3, 0, 5e3]
        v = [0.0, 0.0, 0.0]
        μ = 398600.4354360959

        @test isapprox(
                vectorfield(vcat(r,v), [μ], NaN),
                [0.0, 0.0, 0.0, 0.00248543, -0.0, -0.00112974];
                atol = 1e-8
        )
    end

end


text:
# Basic model stuff
# model types
type AffineModelType <: ModelType end
type TermStructureConsistentModelType <: ModelType end
type CalibratedModelType <: ModelType end # general model type

# type aliases
typealias AffineModel Model{AffineModelType}
typealias TermStructureConsistentModel Model{TermStructureConsistentModelType}
typealias CalibratedModel Model{CalibratedModelType}

type ModelCommon
  observe::ObserverMixin
end

ModelCommon() = ModelCommon(ObserverMixin())

get_observers(m::Model) = m.common.observe.observers


text:
# backward algorithms
######################### Backward Regression Algorithm ########################
struct BackwardRegression{T, AT<:AbstractMatOrFac{T}, B<:AbstractVector{T},
                        V<:AbstractVector{T}, FT, N} <: AbstractMatchingPursuit{T}
    A::AT # matrix
    b::B # target
    r::V # residual
    AiQR::FT # updatable QR factorization of Ai
    δ²::V # increase in SQUARED residual norm
    isfast::N # toggles fast implementation
end
const BR = BackwardRegression
const BackwardGreedy = BackwardRegression
const BOOMP = BR # a.k.a. Backward Optimized OMP

function BR(A::AbstractMatOrFac, b::AbstractVector; isfast::Bool = true)
    n, m = size(A)
    r = zeros(eltype(A), n)
    AiQR = UpdatableQR(A)
    rescaling = colnorms(A)
    δ² = fill(-Inf, m)
    return BR(A, b, r, AiQR, δ², Val(isfast))
end

# calculates a solution to a full column rank (not underdetermined)
# linear system Ax = b with the backward greedy algorithm
# max_ε is the residual norm tolerance
# max_δ is the largest marginal increase in residual norm before the algorithm terminates
# k is the desired sparsity of the solution
# whichever criterion is hit first
function br(A::AbstractMatOrFac, b::AbstractVector, max_ε::Real, max_δ::Real, k::Int; isfast::Bool = true)
    m = size(A, 2)
    P = BR(A, b, isfast = isfast)
    x = sparsevec(1:m, P.AiQR \ b)
    for _ in m:-1:k+1
        backward_step!(P, x, max_ε, max_δ, P.isfast) || break
    end
    return x
end

# keyword version
function br(A::AbstractMatOrFac, b::AbstractVector; max_residual::Real = Inf,
            max_increase::Real = Inf, sparsity::Int = 0, isfast::Bool = true)
     br(A, b, max_residual, max_increase, sparsity, isfast = isfast)
end

function update!(P::BR, x::SparseVector)
    backward_step!(P, x, Inf, Inf)
end

function backward_step!(P::Union{FR, BR}, x::SparseVector, max_ε::Real, max_δ::Real,
                                                        isfast::Val = Val(true))
    nnz(x) > 0 || return false
    residual!(P.r, P.A, x, P.b)
    normr = norm(P.r)
    δ² = isfast isa Val{true} ? backward_δ!(P, x) : naive_backward_δ!(P, x)
    min_δ², i = findmin(δ²) # drop the atom that leads to the minimum increase of the residual norm
    new_norm = sqrt(min_δ² + normr^2) # since min_δ is the squared marginal increase in norm
    if new_norm < max_ε && min_δ² < max_δ^2
        _dropindex!(x, P.AiQR, i) # i is index into x.nzval, NOT into x
        ldiv!!(x.nzval, P.AiQR, P.b, P.r)
        return true
    else
        ldiv!!(x.nzval, P.AiQR, P.b, P.r)
        return false
    end
end

function get_gamma(P::Union{FR, BR})
    F = P.AiQR
    AA = F.R \ Matrix(F.R' \ I)
    ip = invperm(F.perm)
    return γ = diag(AA)[ip]
end

# updates vector of square root of increase in SQUARED norm δ²
# this is done for compatability with RelevanceMatchingPursuit,
# since |r_{-i}|^2 - |r|^2 = |<φ, r>|^2 / |ψ|^2
function backward_δ!(P::Union{FR, BR}, x::SparseVector)
    γ = get_gamma(P)
    @. P.δ²[x.nzind] = x.nzval^2  / γ
    # return P.δ²[x.nzind]
end

# WARNING: assumes P.r == P.b - P.A*x when function is called
# or normr = norm(P.b - P.A*x)
function naive_backward_δ!(P::Union{FR, BR}, x::SparseVector, normr = norm(P.r))
    # reduce all variables to support of x
    A = @view P.A[:, x.nzind]
    δ² = @view P.δ²[x.nzind]
    n = length(x.nzind)
    # y = @view x.nzval[1:n-1] # temporary storage for coefficients of smaller problem
    y = similar(x.nzval, n-1)
    for i in 1:n # could be parallelized if y and P.r are separate for each thread
        a = @view A[:, i]
        Ai = @view A[:, filter(!=(i), 1:length(x.nzind))]
        remove_column!(P.AiQR, i)
        ldiv!!(y, P.AiQR, P.b, P.r)
        P.r .= P.b
        mul!(P.r, Ai, y, -1, 1)
        δ²[i] = norm(P.r)^2 - normr^2 # squared residual norm increase
        add_column!(P.AiQR, a, i)
    end
    return δ²
end

################################################################################
# see "An Efficient Implementation of the Backward Greedy Algorithm for Sparse Signal Reconstruction"
# WARNING: potentially more susceptible to ill-conditioned systems
# except for research purposes, the isfast = true option of BackwardRegression should be used
# IDEA: could implement corresponding forward regression which keeps track of AA⁻¹ instead of QR
struct FastBackwardRegression{T, AT<:AbstractMatOrFac{T}, B<:AbstractVector{T},
            V<:AbstractVector{T}, AAT<:AbstractMatOrFac{T}} <: AbstractMatchingPursuit{T}
    A::AT # matrix
    b::B # target
    r::V # residual
    AA⁻¹::AAT # stores AA⁻¹ corresonding to active set in upper left block
    Ab::V # stores A' * b
    δ²::V # increase in SQUARED residual norm
end
const FBR = FastBackwardRegression
function FBR(A::AbstractMatOrFac, b::AbstractVector)
    r = similar(b)
    return FBR(A, b, r, qr(A))
end
# utilizing the already existing QR-factorization
function FBR(A::AbstractMatOrFac, b::AbstractVector, r::AbstractVector, F::Factorization)
    n, m = size(A)
    AA = F.R \ Matrix(F.R' \ I(m))
    Ab = A'b
    δ² = zeros(m)
    return FBR(A, b, r, AA, Ab, δ²)
end
function FBR(A::AbstractMatOrFac, b::AbstractVector, r::AbstractVector, F::UpdatableQR)
    n, m = size(F)
    E = F.uqr
    AA = E.R1 \ Matrix(E.R1' \ I(m))
    ip = invperm(F.perm)
    AA = AA[ip, ip]
    Ab = A'b
    δ² = zeros(m)
    return FBR(A, b, r, AA, Ab, δ²)
end
# build FBR object from existing ForwardRegression object
function FBR(P::FR)
    FBR(P.A, P.b, P.r, P.AiQR)
end
# TODO: has some overlap with BackwardRegression, which could be consolidated
function fbr(A::AbstractMatOrFac, b::AbstractVector;
     max_residual::Real = Inf, max_increase::Real = Inf, sparsity::Int = 0)
     fbr(A, b, max_residual, max_increase, sparsity)
end

function fbr(A::AbstractMatOrFac, b::AbstractVector, max_ε::Real, max_δ::Real, k::Int)
    m = size(A, 2)
    P = FBR(A, b)
    x = sparsevec(1:m, P.AA⁻¹ * P.Ab)
    for i in m:-1:k+1
        backward_step!(P, x, max_ε, max_δ) || break
    end
    return x
end

function backward_step!(P::FBR, x::SparseVector, max_ε::Real, max_δ::Real)
    nnz(x) > 0 || return false
    normr = norm(residual!(P.r, P.A, x, P.b))
    δ² = backward_δ!(P, x)
    min_δ², i = findmin(δ²) # drop the atom that leads to the minimum increase of the residual norm
    if min_δ² + normr^2 < 0
        println(min_δ² + normr^2)
        println([min_δ², normr^2])
        throw("numerical instability encountered in backward step")
    end
    new_norm = sqrt(min_δ² + normr^2) # since δ is the marginal increase of squared norm
    if new_norm < max_ε && min_δ² < max_δ^2
        _dropindex!(x, P, i) # i is index into x.nzval, NOT into x
        _solve!(P, x)
        return true
    else
        _solve!(P, x)
        return false
    end
end

# solves the least squares problem constrained to non-zero elements of x
# WARNING: assumes P.AA⁻¹ is updated correctly
function _solve!(P::FBR, x::SparseVector)
    AA⁻¹, Ab = @views P.AA⁻¹[1:nnz(x), 1:nnz(x)], P.Ab[x.nzind]
    mul!(x.nzval, AA⁻¹, Ab)
    return x
end

function backward_δ!(P::FBR, x::SparseVector)
    m = nnz(x)
    AA⁻¹ = @views P.AA⁻¹[1:m, 1:m]
    δ² = @view P.δ²[1:m]
    γ = @view AA⁻¹[diagind(AA⁻¹)]
    @. δ² = x.nzval^2 / γ # since x = AA⁻¹ * A' * P.b
end

# drops the ith index of x and updates the matrices AA⁻¹, AA⁻¹A accordingly
function _dropindex!(x::SparseVector, P::FBR, i::Int)
    j = x.nzind[i]
    m = nnz(x)
    deleteat!(x.nzind, i)
    deleteat!(x.nzval, i)
    A, AA⁻¹ = @views P.A[:, x.nzind], P.AA⁻¹[1:m, 1:m]
    active = 1:m .!= i
    G, g, γ = AA⁻¹[active, active], AA⁻¹[i, active], AA⁻¹[i, i]
    @. P.AA⁻¹[1:m-1, 1:m-1] = G - g * g' / γ # inv(A'A)
    return x, P
end

################################# LACE #########################################
# Least Absolute Coefficient Elimination
struct LACE{T, AT<:AbstractMatOrFac{T}, BT<:AbstractVector{T}, V<:AbstractVector{T}, FT}
    A::AT # dictionary
    b::BT # target
    r::V # residual
    AiQR::FT # UpdatableQR
end
function LACE(A::AbstractMatOrFac, b::AbstractVector)
    r = similar(b)
    n, m = size(A)
    n ≥ m || throw("A needs to be overdetermined but is of size ($n, $m)")
    AiQR = UpdatableQR(A)
    LACE(A, b, r, AiQR)
end

function lace(A::AbstractMatOrFac, b::AbstractVector;
        max_residual::Real = Inf, max_increase::Real = Inf, sparsity::Int = 0)
    return lace(A, b, max_residual, max_increase, sparsity)
end

# A is overdetermined linear system, b is target, ε is tolerable residual norm
function lace(A::AbstractMatOrFac, b::AbstractVector, ε::Real, δ::Real, k::Int)
    n, m = size(A)
    L = LACE(A, b)
    x = sparse(ones(eltype(A), size(A, 2)))
    ldiv!!(x.nzval, L.AiQR, L.b, L.r)
    for _ in m:-1:k+1
        backward_step!(L, x, ε, δ) || break
    end
    return x
end
# input x is assumed to be ls-solution with current active set
function update!(P::LACE, x::SparseVector)
    i = argmin(abs, x.nzval) # choose least absolute coefficient magnitude from current support
    _dropindex!(x, P.AiQR, i) # drops ith atom in active set
    ldiv!!(x.nzval, P.AiQR, P.b, P.r)
    return x
end

function backward_step!(P::LACE, x::SparseVector, max_ε::Real, max_δ::Real)
    nnz(x) > 0 || return false
    residual!(P.r, P.A, x, P.b)
    normr = norm(P.r)

    i = argmin(abs, x.nzval) # choose least absolute coefficient magnitude from current support
    j = x.nzind[i] # remember which atom we are deleting
    _dropindex!(x, P.AiQR, i) # i is index into x.nzval, NOT into x
    ldiv!!(x.nzval, P.AiQR, P.b, P.r)

    residual!(P.r, P.A, x, P.b) # new residual
    δ² = norm(P.r)^2 - normr^2 # change in residual magnitude
    newnorm = sqrt(normr^2 + δ²)
    if newnorm < max_ε && δ² < max_δ^2
        return true
    else
        Aj = @view P.A[:,j] # if we can't accept the deletion,
        addindex!(x, P.AiQR, Aj, j) # add back the atom we deleted
        ldiv!!(x.nzval, P.AiQR, P.b, P.r)
        return false
    end
end


text:
#!/usr/bin/env julia
doc""" Script to build the documentation of KullbackLeibler.jl.

- MIT Licensed, (C) 2018 Lilian Besson (Naereen)
  https://GitHub.com/Naereen/KullbackLeibler.jl
"""

__author__ = "Lilian Besson (Naereen)"
__version__ = "0.1"

# https://juliadocs.github.io/Documenter.jl/latest/
using Documenter

include("src/KullbackLeibler.jl")
# include("KullbackLeibler.jl")
using KullbackLeibler

makedocs(
    format = :html,
    sitename = "KullbackLeibler.jl",
    pages = [
        "index.md",
        # "About" => "README.md"
    ],
    # WARNING maybe I need more options right here?
)


text:
using SparseArrays, LightXML, ViewerGL
GL = ViewerGL
include("./../src/Lar.jl")

#=== EXTRACT NODEs AND STREETs FROM OSM FILE ===#
xdoc  = parse_file("examples/map.osm")
xroot = root(xdoc)

node_list = get_elements_by_tagname(xroot, "node")
way_list  = get_elements_by_tagname(xroot, "way")

#=== EXTRACT VERTEX SET AND VERTEX POSITION MAP ===#

Vidxs = [parse(Int, attribute(v, "id")) for v in node_list]
Vmax  = max(Vidxs...)
Vmap  = SparseArrays.spzeros(Int, Vmax)
Vnum  = length(Vidxs)
for i = 1 : Vnum   Vmap[Vidxs[i]] = i   end

lats  = [parse(Float64, attribute(v, "lat")) for v in node_list]
lmin  = min(lats...)
lmax  = max(lats...)
lats  = (lats .- lmin) ./ (lmax - lmin)

lons  = [parse(Float64, attribute(v, "lon")) for v in node_list]
lmin  = min(lons...)
lmax  = max(lons...)
lons  = (lons .- lmin) ./ (lmax - lmin)

#=== BUILD LAR MODEL ===#

model = Lar.Model([lats'; lons'])

for street in way_list
    nodes = [
             Vmap[parse(Int, attribute(n, "ref"))]
             for n in get_elements_by_tagname(street, "nd")
            ]
    for i = 1 : length(nodes) - 1
        newcell = SparseArrays.spzeros(Int8, Vnum)
        newcell[nodes[  i]] = 1
        newcell[nodes[i+1]] = 1
        Lar.addModelCells!(model, 1, convert(Lar.ChainOp, newcell'))
    end
end

#===
    extraction of a box of
     0.31<lats<0.55
     0.65<lons<0.75
===#

toDelV = convert(Array{Int64,1},
    [i  for i = 1 : size(model, 0, 2)
        if model.G[1, i] < 0.31
        || model.G[1, i] > 0.35
        || model.G[2, i] < 0.65
        || model.G[2,i] > 0.75
    ]
)
Lar.deleteModelVertices!(model, toDelV)

GL.VIEW([GL.GLGrid(model.G, Lar.cop2lar(model.T[1]))])

#=== NOW IT LOOPS ===#
LarA = Lar.Arrangement

sigma           = spzeros(Int8, 0)
return_edge_map = false
multiproc       = false

model, edge_map = LarA.planar_arrangement_1(model, sigma, true, multiproc)

bicon_comps = LarA.biconnected_components(model.T[1])
V           = convert(Lar.Points, model.G')
copEV       = model.T[1]

# arrangement of isolated components
n = size(bicon_comps, 1)
shells = Array{Lar.Chain, 1}(undef, n)
boundaries = Array{Lar.ChainOp, 1}(undef, n)
EVs = Array{Lar.ChainOp, 1}(undef, n)
# for each component
for p in 2 : n
    ev = copEV[sort(bicon_comps[p]), :]
    # computation of 2-cells
    fe = LarA.minimal_2cycles(V, ev)
    # exterior cycle
    shell_num = LarA.get_external_cycle(V, ev, fe)
    # decompose each fe (co-boundary local to component)
    EVs[p] = ev
    tokeep = setdiff(1:fe.m, shell_num)
    boundaries[p] = fe[tokeep, :]
    shells[p] = fe[shell_num, :]
end

p = 1
ev = copEV[sort(bicon_comps[p]), :]
# computation of 2-cells
fe = LarA.minimal_2cycles(V, ev)                         # <---| Here it loops
# exterior cycle
shell_num = LarA.get_external_cycle(V, ev, fe)
# decompose each fe (co-boundary local to component)
EVs[p] = ev
tokeep = setdiff(1:fe.m, shell_num)
boundaries[p] = fe[tokeep, :]
shells[p] = fe[shell_num, :]


text:
# 2D plots
"""
    plot([x,] y, [z,] [axes,] [supp,] [dims,] [h,] [curve]) -> Gaston.Figure

Plots data `x`, `y`, and optionally `z` and `supp`, in two or three dimensions
as specified by `dims`, with axes configuration `axes`, and curve configuration
`curve`, on the figure with handle `h`.

If `x` is omitted, then `x=1:length(t)` is assumed.

    plot(x, f::Function, args...) -> Gaston.Figure

Assume `y = f.(x)`.

    plot(c, args...) -> Gaston.Figure

If `c` is a complex vector, plot its real part vs its imaginary part.

    plot(P::Matrix{Union{Gaston.Figure}, Nothing}) -> Gaston.Figure

Create a 'muliplot' with the layout in `P`. Returns a new figure.
"""
function plot(x::Coord, y::Coord, z::Coord = nothing, axes::Axes = Axes() ;
              supp::Coord = nothing,
              dims::Int   = 2,
              handle      = gnuplot_state.current,
              args...)

    ## Process optional keyword arguments.
    axesconf = parse(axes)
    curveconf = parse(args)

    # create curve
    c = Curve(x, y, z, supp, curveconf)

    # push curve we just created to current figure
    # create new figure
    fig = newfigure(handle)
    push!(fig, Plot(dims=dims, axesconf = axesconf, curves = [c]))
    #fig.subplots[1] = Plot(dims=dims, axesconf = axesconf, curves = [c])

    # write gnuplot data to file
    write_data(c, fig.subplots[1].dims, fig.subplots[1].datafile)

    return fig
end

### Alternative `plot` methods
plot(y, a::Axes = Axes() ; args...) = plot(1:length(y), y, nothing, a ; args...)
plot(x, y, a::Axes = Axes() ; args...) = plot(x, y, nothing, a ; args...)
## Complex vector
plot(c::ComplexCoord, a::Axes = Axes() ; args...) = plot(real(c), imag(c), nothing, a ; args...)
## Function
plot(x, f::Function, a::Axes = Axes() ; args...) = plot(x, f.(x), nothing, a ; args...)

# Add a curve to an existing figure
"""
plot!(x, y [, z] [,supp] [, dims] [, h] [, curve...]) -> Gaston.Figure

Adds the curve specified by data `x`, `y`, and optionally `z` and `supp`, in
two or three dimensions as specified by `dims`, and curve configuration
`curve`, to the figure with handle `h`.
"""
function plot!(x::Coord, y::Coord, z::Coord = nothing;
               supp::Coord = nothing,
               dims::Int   = 2,
               handle      = gnuplot_state.current,
               args...)

    # Process optional keyword arguments.
    curveconf = parse(args)

    # determine handle
    handles = gethandles()
    handle ∈ handles || error("Cannot append curve to non-existing handle")

    # create curve
    c = Curve(x, y, z, supp, curveconf)

    # push new curve to current figure
    fig = gnuplot_state.figs[findfigure(handle)]
    push!(fig, c)

    # write gnuplot data to a file
    write_data(c, fig.subplots[1].dims, fig.subplots[1].datafile, append=true)

    return fig
end

### Alternative `plot!` methods
# Omit `x` coordinate
plot!(y ; args...) = plot!(1:length(y), y ; args...)
# Complex vector
plot!(c::ComplexCoord ; args...) = plot!(real(c), imag(c) ; args...)
# Function
plot!(x, f::Function ; args...) = plot!(x, f.(x) ; args...)

### Multiplot
function plot(P::FigArray ; handle::Int = 0)
    # Create new figure
    fig = newfigure(handle)
    if P isa Vector
        fig.layout = (size(P)[1], 1)
    else
        fig.layout = size(P)
    end
    for p in P
        if p isa Figure
            push!(fig, p.subplots[1])
        else
            push!(fig, nothing)
        end
    end
    return fig
end


text:
mutable struct Author <: Model
    id::DBId{Integer}
    name::String
    age::Integer
    country::String
    date::DateTime
end
function Author(;id::Union{Integer, Nothing} = nothing,
                name::String="",
                age::Integer=0,
                country::String="",
                date::DateTime=now())
    return Author(id, name, age, country, date)
end
mutable struct Book <: Model
    id::DBId{String}
    title::String
    author::ForeignKey{Author}
    data::Dict{String, Integer}
end
function Book(;id::Union{String, Nothing}=nothing,
               title::String="",
               author::Foreign{Author}=Author(),
               data::Dict{String, Integer}=Dict())
    return Book(id, title, author, data)
end

function _test_on_delete()
    mapper = DBMapper(creator)
    register!(mapper, Author)
    register!(mapper, Book)
    configure_relation(mapper, Book, :author, on_delete=Cascade())

    author = Author(name="Author 1", age=50)
    insert!(mapper, Author(name="Author 1", age=50))
    insert!(mapper, Author(name="Author 2", age=50))
    insert!(mapper, Author(name="Author 3", age=50))
    insert!(mapper, Author(name="Author 4", age=50))

    author1 = select_one(mapper, Author, name="Author 1")
    author3 = select_one(mapper, Author, name="Author 3")

    insert!(mapper, Book(author=author1, title="Book 1", data=Dict("key"=>"Some 1")))
    insert!(mapper, Book(author=author1, title="Book 2", data=Dict("key"=>"Some 2")))
    insert!(mapper, Book(author=author1, title="Book 3", data=Dict("key"=>"Some 3")))

    insert!(mapper, Book(author=author3, data=Dict("key"=>"Some 3")))
    delete!(mapper, author3)

    #select_all(mapper, Book, author=author1)

end
function _test_basic_functionalities(creator)
    mapper = DBMapper(creator)

    register!(mapper, Author)
    register!(mapper, Book)
    configure_relation(mapper, Book, :author, on_delete=Cascade())
    @test haskey(mapper.tables, Author)
    @test haskey(mapper.tables, Book)

    create_table(mapper, Author)
    create_table(mapper, Book)

    author = Author(name="pirulo", age=50, country="Argentina")
    insert!(mapper, author)
    @test !isnothing(author.id.x)
    id = author.id.x

    author = Author(name="Author 1", age=3, country="Argentina")
    insert!(mapper, author)

    author = Author(name="Author 2", age=3, country="Brasil")
    insert!(mapper, author)

    author = Author(name="Author 5", age=25, country="Italia")
    insert!(mapper, author)

    author = Author(name="Author 3", age=3, country="Uruguay")
    insert!(mapper, author)

    @test StructDatabaseMapping.exists(mapper, Author, name="Enrique Banch") == false
    @test StructDatabaseMapping.exists(mapper, Author, name="pirulo") == true
    @test StructDatabaseMapping.exists(mapper, author) == true

    @test StructDatabaseMapping.exists(mapper, Author, name="Author 3", age=3) == true
    @test StructDatabaseMapping.exists(mapper, Author, name="Author 3", age=4) == false
    @test StructDatabaseMapping.exists(mapper, Author, pk=author.id.x, age=3) == true
    @test StructDatabaseMapping.exists(mapper, Author, pk=author.id.x, age=4) == false

    authors = select_all(mapper, Author, age=3)
    @test isa(authors[1], Author)
    @test length(authors) == 3

    authors = select_all(mapper, Author)
    @test length(authors) == 5

    authors = select_all(mapper, Author, age=3, fields=[:name, :country])
    @test Set([author.country for author in authors]) == Set(["Brasil", "Argentina", "Uruguay"])

    a = select_one(mapper, Author, id=999)
    @test isnothing(a)
    a = select_one(mapper, Author, id=id)
    @test a.name == "pirulo"

    a.name = "otro_pirulo"
    a.age = 5
    update!(mapper, a; fields=[:name])
    author = select_one(mapper, Author, id=id)
    @test author.name == "otro_pirulo"
    @test author.age == 50

    author.name = "some_other_name"
    author.age = 5
    update!(mapper, author)
    author = select_one(mapper, Author, id=id)
    @test author.name == "some_other_name"
    @test author.age == 5

    book = Book(id="super_string_id", author=author,
                data=Dict{String, Integer}("some_data"=>5))
    insert!(mapper, book)

    a = select_one(mapper, Book, id="bbb")
    @test isnothing(a)
    a = select_one(mapper, Book, id="super_string_id")

    @test a.id.x == "super_string_id"
    @test get(a.author, mapper).name == "some_other_name"
    @test a.data["some_data"] == 5

    book = select_one(mapper, Book, id="super_string_id")
    delete!(mapper, book)

    author = select_one(mapper, Author, id=id)
    @test author.name == "some_other_name"
    @test author.age == 5
    delete!(mapper, Author, name="some_other_name")

    author = select_one(mapper, Author, id=id)
    @test isnothing(author)

end

max_stars_repo_name: lucianolorenti/StructDatabaseMapping

max_issues_repo_name: lucianolorenti/StructDatabaseMapping

max_forks_repo_name: lucianolorenti/StructDatabaseMapping


text:
__precompile__()

"""
A Julia package for dynamical billiard systems in two dimensions.

The goals of the package is to provide a flexible, easy-to-use
and intuitive framework for
fast implementation of billiard systems of arbitrary construction.
"""
module DynamicalBilliards

using LinearAlgebra
using StaticArrays
using Elliptic
import Base: show, eltype, getindex

const SV = SVector{2}
export SVector

cossin(a) = ((x, y) = sincos(a); (y, x))

##########################################
# Core                                   #
##########################################
include("billiards/particles.jl")
include("billiards/obstacles.jl")
include("billiards/billiardtable.jl")
include("billiards/standard_billiards.jl")

include("timeevolution/collisions.jl")
include("timeevolution/propagation.jl")
include("timeevolution/timeseries.jl")
include("timeevolution/highleveltimes.jl")

##########################################
# Advanced                               #
##########################################
include("boundary/boundarymap.jl")
include("boundary/phasespacetools.jl")

include("poincare.jl")
include("lyapunov_spectrum.jl")

include("mushroomtools.jl")
export MushroomTools

include("raysplitting.jl")
include("parallel.jl")
include("testing.jl")

####################################################
# Plotting Routines (loaded when `Using PyPlot`)   #
####################################################
using Requires
function __init__()
    @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" begin
        import .PyPlot: plot
        dir = joinpath(@__DIR__, "plotting")
        for f in readdir(dir)
            include(joinpath(dir, f))
        end
    end
end

###################
# Update messages #
###################
# if !isfile(joinpath(@__DIR__, "update_v3.0.0"))
# printstyled(stdout,
# """
# \nUpdate message: DynamicalBilliards v3.0
#
# The new version v3.0 of DynamicalBilliards has
# a reworked (and better!) propagation
# algorithm, a new sexy way to animate billiards
# (multiple particle support!)
# the Ellipse obstacle and many other things!
# Please see the changelog,
# because there have been a small
# amount of breaking changes!\n
# """; color = :light_magenta)
# touch(joinpath(@__DIR__, "update_v3.0.0"))
# end

end#module


text:
# lambdify an expression

## Mapping of SymPy Values into julia values
val_map = Dict(
               "Zero"             => :(0),
               "One"              => :(1),
               "NegativeOne"      => :(-1),
               "Half"             => :(1/2),
               "Pi"               => :pi,
               "Exp1"             => :ℯ,
               "Infinity"         => :Inf,
               "NegativeInfinity" => :(-Inf),
               "ComplexInfinity"  => :Inf, # error?
               "ImaginaryUnit"    => :im,
               "BooleanTrue"      => :true,
               "BooleanFalse"     => :false
               )

## Mapping of Julia function names into julia ones
## most are handled by Symbol(fnname), this catches exceptions
_heaviside(x) = 1//2 * (1 + sign(x))
function _piecewise(args...)
    as = copy([args...])
    val, cond = pop!(as)
    ex = Expr(:call, :ifelse, cond, val, :nothing)
    while length(as) > 0
        val, cond = pop!(as)
        ex = Expr(:call, :ifelse, cond, val, ex)
    end
    ex
end

## Hack to avoid Expr(:call,  :*,2, x)  being  2x and  not  2*x
__prod__(args...) =  prod(args)
export __prod__

fn_map = Dict(
              "Add" => :+,
              "Sub" => :-,
              "Mul" => :__prod__,
              "Div" => :/,
              "Pow" => :^,
              "re"  => :real,
              "im"  => :imag,
              "Abs" => :abs,
              "Min" => :min,
              "Max" => :max,
              "Poly" => :identity,
              "Heaviside" => :(_heaviside),
              "Piecewise" => :(_piecewise),
              "And" => :(&),
              "Or" => :(|),
              "Less" => :(<),
              "LessThan" => :(<=),
              "StrictLessThan" => :(<),
              "Equal" => :(==),
              "Equality" => :(==),
              "Unequality" => :(!==),
              "StrictGreaterThan" => :(>),
              "GreaterThan" => :(>=),
"Greater" => :(>),
"conjugate" => :conj
              )

map_fn(key, fn_map) = haskey(fn_map, key) ? fn_map[key] : Symbol(key)

Base.convert(::Type{Expr}, x::SymbolicObject) = walk_expression(x)

function walk_expression(ex; values=Dict(), fns=Dict())

    fns_map = merge(fn_map, fns)
    vals_map = merge(val_map, values)

    fn = Introspection.funcname(ex)

    # special case `F(t) = ...` output from ODE
    # this may be removed if it proves a bad idea....
    if fn == "Equality" && lhs(ex).is_Function
        return walk_expression(rhs(ex), values=values, fns=fns)
    end

    if fn == "Symbol" || fn == "Dummy"
        str_ex = string(ex)
        return get(vals_map, str_ex, Symbol(str_ex))
    elseif fn in ["Integer" , "Float"]
        return N(ex)
    elseif fn == "Rational"
        return convert(Int, numer(ex))//convert(Int, denom(ex))
        ## piecewise requires special treatment
    elseif fn == "Piecewise"
        return _piecewise([walk_expression(cond, values=values, fns=fns) for cond in Introspection.args(ex)]...)
    elseif fn == "ExprCondPair"
        val, cond = Introspection.args(ex)
        return (val, walk_expression(cond, values=values, fns=fns))
    elseif fn == "Tuple"
        return walk_expression.(Introspection.args(ex), values=values, fns=fns)
    elseif haskey(vals_map, fn)
        return vals_map[fn]
    end

    as = Introspection.args(ex)

    Expr(:call, map_fn(fn, fns_map), [walk_expression(a, values=values, fns=fns) for a in as]...)
end

"""
     lambdify(ex, vars; typ, fns, values, use_julia_code=false, invoke_latest=true)

Take a symbolic expression and returns an anonymous `Julia` function.

Converts from a SymPy object to a `Julia` expression by walking the SymPy expression tree and converting each step, (e.g. essentially calls `convert(Expr,  ex)`).

Then creates a function. The function arguments are based on `free_symbols` and its ordering, unless `vars` is specified directly.

* `use_julia_code=true` will use SymPy's conversion to an expression, the default is `false`

* `invoke_latest=true` calls `Base.invokelatest` to work around world age issues. This is the safe default, but setting to `false` will result in faster-executing functions.

Example:

```jldoctest
julia> using SymPy

julia> @vars x y z
(x, y, z)

julia> ex = x^2 * sin(x)
 2
x ⋅sin(x)

julia> fn = lambdify(ex)
#87 (generic function with 1 method)

julia> fn(pi)
1.2086779438644711e-15

julia> ex = x + 2y + 3z
x + 2⋅y + 3⋅z

julia> fn = lambdify(ex)
#87 (generic function with 1 method)

julia> fn(1,2,3) # order is by free_symbols
14

julia> ex(x=>1, y=>2, z=>3)
14

julia> fn = lambdify(ex, (y,x,z))
#87 (generic function with 1 method)

julia> fn(1,2,3)
13
```

"""
function lambdify(ex::Sym, vars=free_symbols(ex);
                  typ=Any, fns=Dict(), values=Dict(),
                  use_julia_code=false,
                  invoke_latest=true
                  )
    # if :julia_code printer is there, use it
    if use_julia_code
        body = Meta.parse(sympy.julia_code(ex)) # issue here with 2.*...
    else
        body = walk_expression(ex, fns=fns, values=values)
    end

    try
        syms = typ == Any ? map(Symbol,vars) : map(s->Expr(:(::),s,typ), Symbol.(vars))
        fn = eval(Expr(:function, Expr(:call, gensym(), syms...), body))
        if invoke_latest
            (args...) -> Base.invokelatest(fn, args...) # https://github.com/JuliaLang/julia/pull/19784
        else
            fn
        end
    catch err
        throw(ArgumentError("Expression does not lambdify"))
    end
end

function _lambdify(ex::Sym, vars=free_symbols(ex))
    body = convert(Expr, ex)
    syms = Symbol.(vars)
    fn = eval(Expr(:function, Expr(:call, gensym(), syms...), body))
    fn
end

# from @mistguy cf. https://github.com/JuliaPy/SymPy.jl/issues/218
# T a data type to convert to, when specified
function lambdify(exs::Array{S, N}, vars = union(free_symbols.(exs)...); T::DataType=Nothing, kwargs...) where {S <: Sym, N}
    f = lambdify.(exs, (vars,)) # prevent broadcast in vars
    if T == Nothing
        (args...) -> map.(f, args...)
    else
        (args...) -> convert(Array{T,N}, map.(f, args...))
    end
end

export(lambdify)

Base.convert(::Type{Function}, ex::Sym) = lambdify(ex)


text:
export JDFFile, @jdf_str, path, getindex

import Base: getindex, view

"""
    jdf"path/to/JDFfile.jdf"

    JDFFile("path/to/JDFfile.jdf")

Define a JDF file, which you use with methods like  `names` and `size`.

## Example
using JDF, DataFrames
df = DataFrame(a = 1:3, b = 1:3)
savejdf(df, "plsdel.jdf")

names(jdf"plsdel.jdf") # [:a, :b]
size(jdf"plsdel.jdf") # (2, 3)

size(jdf"plsdel.jdf", 1) # (2, 3)

size(jdf"plsdel.jdf", 1) # (2, 3)

# clean up
rm("plsdel.jdf", force = true, recursive = true)
"""
struct JDFFile{T<:AbstractString}
    path::T
end

"""
    jdf"path/to/JDFfile.jdf"

    JDFFile("path/to/JDFfile.jdf")

Define a JDF file, which you can apply `names` and `size`.

## Example
using JDF, DataFrames
df = DataFrame(a = 1:3, b = 1:3)
savejdf(df, "plsdel.jdf")

names(jdf"plsdel.jdf") # [:a, :b]
ncol(jdf"plsdel.jdf") # 2
size(jdf"plsdel.jdf") # (2, 3)

size(jdf"plsdel.jdf", 1) # (2, 3)

size(jdf"plsdel.jdf", 1) # (2, 3)

# clean up
rm("plsdel.jdf", force = true, recursive = true)
"""
macro jdf_str(path)
    return :(JDFFile($path))
end

"""
    path(jdf::JDFFile)

Return the path of the JDF
"""
path(jdf) = getfield(jdf, :path)

function Base.getindex(file::JDFFile, col::String)
   getindex(file, Symbol(col))
end

function Base.getindex(file::JDFFile, col::Symbol)
    # TODO make it load from column loader for faster access
    JDF.load(file; cols = [col])[col]
end

# function Base.getindex(file::JDFFile, rows, col::String)
#     # TODO make it load from column loader for faster access
#     getfield(JDF.load(file; cols = [col]), Symbol(col))[rows]
# end

# function Base.getindex(file::JDFFile, rows, cols::AbstractVector{String})
#     JDF.load(file; cols = cols)[rows, :]
# end

# Base.view(file::JDFFile, rows, cols) = getindex(file, rows, cols)

# getindex(file::JDFFile, rows, cols) = JDF.load(file)[rows, cols]

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
# This is an attempt to include the required for n-body simulations fields into structures
# Indeed, tools of ConcreteAbstractions.jl seem to be more suitable for the fields inhereting
DiffEqBase.@def position_velocity_mass begin
    r::SVector{3,cType}
    v::SVector{3,cType}
    m::mType
end

abstract type Body
end

struct MassBody{cType<:Real, mType<:Real} <: Body
    @position_velocity_mass
end

struct ChargedParticle{cType<:Real, mType<:Real, qType<:Real} <: Body
    @position_velocity_mass
    q::qType
end

struct MagneticParticle{cType<:Real, mType<:Real, mmType<:Real} <: Body
    @position_velocity_mass
    mm::SVector{3,mmType}
end

struct WaterMolecule <: Body
    O::MassBody
    H1::MassBody
    H2::MassBody
end

function generate_bodies_in_cell_nodes(n::Integer, m::Real, v_dev::Real, L::Real; T=Float64)
    rng = MersenneTwister(n)
    velocities = v_dev * randn(rng, T, (3,n))
    bodies = MassBody[]
    count = 1
    dL = T(L / (ceil(n^(1/3))))
    for x = dL/2:dL:L, y = dL/2:dL:L, z = dL/2:dL:L
        if count > n
            break
        end
        r = SVector(x, y, z)
        v = SVector{3,T}(velocities[:,count])
        body = MassBody{T,T}(r, v, m)
        push!(bodies, body)
        count += 1
    end
    return bodies
end


text:
using Plots
using GDAL
using Shapefile
using StatsBase
using CSV
using DataFrames
using Statistics
using GBIF

cd("/Users/kiristern/Documents/GitHub/SDM")
include("../BioClim/src/lib/SDMLayer.jl")
include("../BioClim/src/lib/gdal.jl")
include("../BioClim/src/lib/worldclim.jl")
include("../BioClim/src/lib/bioclim.jl")
include("../BioClim/src/lib/shapefiles.jl")
include("explo_fnc.jl")

# Import csv file
df = CSV.read("data/prey_class/prey_class_onedf.csv", header=true)

#remove rows with missing values
df_drop = dropmissing(df, [:decimalLatitude, :decimalLongitude, :year])

#new df of just genus, occurence, year
df_new = select(df_drop, [:class, :decimalLatitude, :decimalLongitude, :year])

# Rename coordinate names
rename!(df_new, :decimalLatitude => :latitude)
rename!(df_new, :decimalLongitude => :longitude)
#for simplicity sake of running Tim's predictive model
rename!(df_new, :class => :species)

# convert longitude "string" into "Float64"
#df_new.longitude = parse.(Float64, df_new.longitude)

# prepare data according to explo_fnc function
# df = prepare_csvdata(df)

df_new = dropmissing(df_new, [:species])

#select for one genus at a time
prey = [df_new[df_new.species .== u,:] for u in unique(df_new.species)]
occ = prey[2]

# Get the worldclim data by their layer number
cd("/Users/kiristern/Documents/GitHub/BioClim/")

@info "Extract and crop bioclim variables"
@time wc_vars = [clip(worldclim(i), occ) for i in 1:19];
# Make the prediction for each layer
@info "Predictions for each layer"
@time predictions = [bioclim(wc_vars[i], occ) for i in 1:length(wc_vars)];
# Make the final prediction by taking the minimum
@info "Minimum-consensus aggregation"
@time prediction = reduce(minimum, predictions);
# Get the threshold for NaN given a percentile
@info "Threshold estimation"
@time threshold = first(quantile(prediction[occ], [0.05]))
@info "5% threshold:\t$(round(threshold; digits=3))"
# Filter the predictions based on the threshold
@info "Final prediction filtering"
@time for i in eachindex(prediction.grid)
    prediction.grid[i] < threshold && (prediction.grid[i] = NaN)
end

worldmap = clip(worldshape(50), prediction)

sdm_plot = plot([0.0], lab = " ", msw = 0.0, ms = 0.0, size = (900, 450), frame=:box,
                title = first(unique(occ.species)))
xaxis!(sdm_plot, (prediction.left, prediction.right), "Longitude")
yaxis!(sdm_plot, (prediction.bottom, prediction.top), "Latitude")

for p in worldmap
    sh = Shape([pp.x for pp in p.points], [pp.y for pp in p.points])
    plot!(sdm_plot, sh, c = :lightgrey, lab = " ")
end

heatmap!(
    sdm_plot,
    longitudes(prediction), latitudes(prediction), prediction.grid,
    aspectratio = 92.60/60.75, c = :BuPu,
    clim = (0.0, maximum(filter(!isnan, prediction.grid)))
    )

for p in worldmap
    xy = map(x -> (x.x, x.y), p.points)
    plot!(sdm_plot, xy, c=:grey, lab = " ", leg=false)
end
sdm_plot

savefig("../SDM/figures/sdm_map-$(first(unique(occ.species))).pdf")

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
using HTTP
using Random

"HTTPUpgradeResponse is a collection of data returned by a HTTP upgrade, regardless of package."
struct HTTPUpgradeResponse
    io::IO
    statuscode::Int
    headers::HeaderList
    excess::AbstractVector{UInt8}
end

"""
HTTPAdapter is an abstract type for an HTTP package that can be used for a handshake GET.

Any subtype of HTTPAdapter should implement the `dohandshake` method.
"""
abstract type HTTPAdapter end

"""
    dohandshake(http::HTTPAdapter, uri::String, headers::HeaderList)

Do an HTTP GET request to `uri` including headers `headers`.

The `headers` list will contain all WebSocket upgrade specific headers, such as `Connection`,
`Upgrade`, and `Sec-WebSocket-Key`.
"""
dohandshake(::HTTPAdapter, uri::String, headers::HeaderList) :: HTTPUpgradeResponse = error("Implement this in your subtype")

"AbstractHandshakeResult is a  supertype for either a successful or an unsuccessful handshake."
abstract type AbstractHandshakeResult end

"GoodHandshake is returned when a successful handshake has been made."
struct GoodHandshake <: AbstractHandshakeResult
    # The socket to use for all WebSocket communication.
    io::IO

    # Any excess bytes read from the socket during the handshake, that should actually be part of
    # the WebSocket communication.
    excess::AbstractVector{UInt8}
end

"BadHandshake is returned when an unsuccessful handshake has been made."
struct BadHandshake <: AbstractHandshakeResult end

"""
    issuccessful(::AbstractHandshakeResult)

True if a handshake was succcessful, false otherwise.
"""
issuccessful(::GoodHandshake) = true
issuccessful(::BadHandshake) = false

"""WebSocketHandshake represents a way to make a WebSocket connection handshake.

The only handshake detailed in the specification is an HTTP handshake, represented by the
`HTTPHandshake` type (regardless of HTTP package used to implement it). The specification does
mention the possibility of other types of handshakes, even though these would be outside of the
specification.
"""
abstract type WebSocketHandshake end

"""HTTPHandshake implements the HTTP handshake from the WebSocket specification.

It requires a random number generator to generate the WebSocket random key used to verify that the
other side is actually a WebSocket server. The `HTTPAdapter` supplied implements the actual HTTP
GET request, which returns an `HTTPUpgradeResponse`.
"""
struct HTTPHandshake <: WebSocketHandshake
    handshakelogic::HTTPHandshakeLogic
    http::HTTPAdapter

    HTTPHandshake(rng::Random.AbstractRNG, http::HTTPAdapter) = new(HTTPHandshakeLogic(rng), http)
end

"""
    performhandshake(h::HTTPHandshake, uri::String)

Do a handshake with the server at `uri`, with parameters supplied by `h`. Validate the handshake,
and return a good or bad result, depending on the validation.
"""
function performhandshake(h::HTTPHandshake, uri::String) :: AbstractHandshakeResult
    try
        upgraderesponse = dohandshake(h.http, uri, getrequestheaders(h.handshakelogic))
        validation = validateresponse(h.handshakelogic, upgraderesponse.statuscode, upgraderesponse.headers)
        if issuccessful(validation)
            GoodHandshake(upgraderesponse.io, upgraderesponse.excess)
        else
            BadHandshake()
        end
    catch ex
        @error "Exception during WebSocket handshake:", ex
        BadHandshake()
    end
end

#
# Integration with the WebSocketUpgrade module.
#

struct WebSocketUpgradeAdapter <: HTTPAdapter end

function dohandshake(::WebSocketUpgradeAdapter, uri::String, headers::HeaderList) :: HTTPUpgradeResponse
    response, socket = WebSocketUpgrade.websocketupgrade(uri, headers)
    # TODO: Handle excess when available
    HTTPUpgradeResponse(socket, response.status, response.headers, b"")
end

#
# Integration with HTTP.jl
#
# This implements the handshake HTTP request using the HTTP.jl package.
#

# Side note: The name `HTTPjlAdapter` isn't great, but not worse than `HTTPHTTPAdapter`.
struct HTTPjlAdapter <: HTTPAdapter end

function dohandshake(::HTTPjlAdapter, uri::String, headers::HeaderList) :: HTTPUpgradeResponse
    socket, response, excess = HTTP.openraw("GET", uri, headers)
    responseheaders = [String(a) => String(b) for (a, b) in response.headers]

    HTTPUpgradeResponse(socket, response.status, response.headers, excess)
end


text:
module DrawSimpleGraphs

using Plots, SimpleGraphs, SimpleDrawing

import SimpleDrawing: newdraw, finish
export newdraw, finish, draw, draw_nodes, draw_edges
export draw_labels

function draw_one_edge(a, b, c, d, hue = "black")
    draw_segment(a, b, c, d, color = hue)
end
"""
    draw_edges(G)
Draw the edges of `G`.
"""
function draw_edges(G::SimpleGraph)
    xy = getxy(G)
    hue = get_line_color(G)
    for ee in elist(G)
        u, v = ee
        draw_one_edge(xy[u][1], xy[u][2], xy[v][1], xy[v][2], hue)
    end
    finish()
end

function draw_one_node(x, y, hue = "black", fill = "white", node_size = 6)
    draw_point(
        x,
        y,
        markerstrokecolor = hue,
        markercolor = fill,
        marker = node_size,
        markerstrokewidth = 1,
    )
end

"""
    draw_nodes(G)
Draw the vertices of `G`.
"""
function draw_nodes(G::SimpleGraph)
    hue = get_line_color(G)
    node_size = get_vertex_size(G)
    xy = getxy(G)
    for v in vlist(G)
        x, y = xy[v]
        fill = get_vertex_color(G, v)
        draw_one_node(x, y, hue, fill, node_size)
    end
    finish()
end

"""
    draw(G)
Draws `G` in its current embedding.
(If `G` does not have an embedding, then it is given a circular
embedding.)

Optional second argument, `clear_first`. If set to `false` the drawing
window is not erased before drawing the graph.
"""
function SimpleDrawing.draw(G::SimpleGraph, clear_first::Bool = true)
    ensure_embed(G)
    if clear_first
        newdraw()
    end
    draw_edges(G)
    draw_nodes(G)
    finish()
end

function draw_labels(G::SimpleGraph, fontsize = 8)
    ensure_embed(G)
    xy = getxy(G)
    for v in G.V
        x, y = xy[v]
        annotate!((x, y, string(v), fontsize))
    end
    finish()
end

include("KnightTourDrawing.jl")

end  # end of module DrawSimpleGraphs


text:
using Plots
using LinearAlgebra
using Statistics
using Random
using Distributions
using PDMats
using PDMatsExtras
using ShiftedArrays
using Parameters
using Polynomials

using DelimitedFiles
using DataFrames

# Import data
data = readdlm("test/otrok_whiteman/real_data/data_deseasoned.csv", ',')
regstatematch = readdlm("test/otrok_whiteman/real_data/regstatematch.csv", ',')

# Data as matrix
datamat = Matrix(data[2:end, :])     # remove first row containing col names
datamat = Float64.(datamat)         # transform entires to Float64
regstatematch = Matrix(regstatematch[2:end, :])
regassign = Int64.(regstatematch[:, 3])

# Estimate HDFM

nlevels = 2

nvar = size(datamat)[2]

nfactors = [1, max(regassign...)]

fassign = ones(Int, nvar, nlevels)
fassign[:, 2] = regassign

flags = [3, 3]

varlags = 3 * ones(Int, nvar)

hdfmpriors = HDFMStruct(nlevels = nlevels,
    nfactors = nfactors,
    factorassign = fassign,
    factorlags = flags,
    errorlags = varlags,
    ndraws = 1000,
    burnin = 50)

results = OW2LevelEstimator(datamat, hdfmpriors)

medians = Any[]
quant33 = Any[]
quant66 = Any[]
stds = Any[]

j = 1
for i in 1:size(results.F)[1]
    push!(stds, std(results.F[i, j, :]))
    push!(quant33, quantile(results.F[i, j, :], 0.05))
    push!(quant66, quantile(results.F[i, j, :], 0.95))
    push!(medians, median(results.F[i, j, :]))
end

plot(results.means.F[:, j])
plot(medians)
plot!(quant33)
plot!(quant66)
plot!(results.means.F[:, j] - stds)
plot!(results.means.F[:, j] + stds)

j = 2
plot(results.means.F[:, 1+j])
plot(results.means.F[:, 1])
plot!(results.means.F[:, 1+j])


text:
#####
##### 'Tupled closure' implementation: 1-tuple, 2-tuple, and then n-tuple by induction
#####

function closure_summary(closures::Tuple, padchar="│")
    Nclosures = length(closures)
    if Nclosures == 1
        return string("Tuple with 1 closure:", '\n',
                      "$padchar   └── ", summary(closures[1]))
    else
        return string("Tuple with $Nclosures closures:", '\n',
         Tuple(string("$padchar   ├── ", summary(c), '\n') for c in closures[1:end-1])...,
                      "$padchar   └── ", summary(closures[end]))
    end
end

#####
##### Kernel functions
#####

funcs     = [:∂ⱼ_τ₁ⱼ, :∂ⱼ_τ₂ⱼ, :∂ⱼ_τ₃ⱼ, :∇_dot_qᶜ, :maybe_tupled_ivd_upper_diagonal, :maybe_tupled_ivd_lower_diagonal, :maybe_tupled_implicit_linear_term]
alt_funcs = [:∂ⱼ_τ₁ⱼ, :∂ⱼ_τ₂ⱼ, :∂ⱼ_τ₃ⱼ, :∇_dot_qᶜ, :ivd_upper_diagonal, :ivd_lower_diagonal, :implicit_linear_term]

for (f, alt_f) in zip(funcs, alt_funcs)
    @eval begin
        @inline $f(i, j, k, grid, closures::Tuple{<:Any}, Ks, args...) =
                    $alt_f(i, j, k, grid, closures[1], Ks[1], args...)

        @inline $f(i, j, k, grid, closures::Tuple{<:Any, <:Any}, Ks, args...) = (
                    $alt_f(i, j, k, grid, closures[1], Ks[1], args...)
                  + $alt_f(i, j, k, grid, closures[2], Ks[2], args...))

        @inline $f(i, j, k, grid, closures::Tuple{<:Any, <:Any, <:Any}, Ks, args...) = (
                    $alt_f(i, j, k, grid, closures[1], Ks[1], args...)
                  + $alt_f(i, j, k, grid, closures[2], Ks[2], args...)
                  + $alt_f(i, j, k, grid, closures[3], Ks[3], args...))

        @inline $f(i, j, k, grid, closures::Tuple{<:Any, <:Any, <:Any, <:Any}, Ks, args...) = (
                    $alt_f(i, j, k, grid, closures[1], Ks[1], args...)
                  + $alt_f(i, j, k, grid, closures[2], Ks[2], args...)
                  + $alt_f(i, j, k, grid, closures[3], Ks[3], args...)
                  + $alt_f(i, j, k, grid, closures[4], Ks[4], args...))

        @inline $f(i, j, k, grid, closures::Tuple, Ks, args...) = (
                    $alt_f(i, j, k, grid, closures[1], Ks[1], args...)
                  + $f(i, j, k, grid, closures[2:end], Ks[2:end], args...))
    end
end

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

with_tracers(tracers, closure_tuple::Tuple) = Tuple(with_tracers(tracers, closure) for closure in closure_tuple)

function calculate_diffusivities!(diffusivity_fields_tuple, closure_tuple::Tuple, args...)
    for (α, closure) in enumerate(closure_tuple)
        @inbounds diffusivity_fields = diffusivity_fields_tuple[α]
        calculate_diffusivities!(diffusivity_fields, closure, args...)
    end
    return nothing
end

function add_closure_specific_boundary_conditions(closure_tuple::Tuple, bcs, args...)
    # So the last closure in the tuple has the say...
    for closure in closure_tuple
        bcs = add_closure_specific_boundary_conditions(closure, bcs, args...)
    end
    return bcs
end

#####
##### Compiler-inferrable time_discretization for tuples
#####

const ETD = ExplicitTimeDiscretization
const VITD = VerticallyImplicitTimeDiscretization

@inline combine_time_discretizations(disc) = disc
@inline combine_time_discretizations(::ETD, ::VITD)  = VerticallyImplicitTimeDiscretization()
@inline combine_time_discretizations(::VITD, ::ETD)  = VerticallyImplicitTimeDiscretization()
@inline combine_time_discretizations(::VITD, ::VITD) = VerticallyImplicitTimeDiscretization()
@inline combine_time_discretizations(::ETD, ::ETD)   = ExplicitTimeDiscretization()

@inline combine_time_discretizations(d1, d2, other_discs...) =
    combine_time_discretizations(combine_time_discretizations(d1, d2), other_discs...)

@inline time_discretization(closures::Tuple) = combine_time_discretizations(time_discretization.(closures)...)


text:
include("DislocationUtil.jl")
include("FEMUtil.jl")

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
export  excelSubpopReport

function excelSubpopReport( mpSim::MPsim, timeGrid::Vector{T},
    subpops::Subpopulation...; filename::AbstractString = "compReport",
    overwrite::Bool = true, timeFactor::Real = 12.0 ) where T <: Real

    filename, overwrite = setupFile( filename, overwrite )
    tStart = now()
    subpopReport = subpopulationPopReport( mpSim, timeGrid, subpops... )
    reportGenerationTime = ( now() - tStart ).value / 1000.0

    if isempty( subpopReport )
        return
    end  # if isempty( subpopReport )

    tStart = now()

    XLSX.openxlsx( filename, mode = overwrite ? "w" : "rw" ) do xf
        dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
        ws = xf[ 1 ]
        nReports = 1

        # Prep the new sheet.
        if overwrite
            XLSX.rename!( ws, string( "Report 1  ", dateStr ) )
        else
            nReports = length( filter( name -> startswith( name, "Report" ),
                XLSX.sheetnames( xf ) ) ) + 1
            ws = XLSX.addsheet!( xf, string( "Report ", nReports, "  ",
                dateStr ) )
        end  # if overwrite

        # Sheet header
        ws[ "A1" ] = "Simulation length"
        ws[ "B1" ] = mpSim.simLength / timeFactor
        ws[ "C1" ] = "years"
        ws[ "A2" ] = "Data generation time"
        ws[ "B2" ] = reportGenerationTime
        ws[ "C2" ] = "seconds"
        ws[ "A3" ] = "Excel generation time"
        ws[ "C3" ] = "seconds"

        subpopReport[ :, :timePoint ] /= timeFactor
        nrows, ncols = size( subpopReport )
        ii = 0

        for name in names( subpopReport )
            ii += 1
            ws[ 5, ii ] = string( name )

            for jj in 1:nrows
                ws[ jj + 5, ii ] = subpopReport[ jj, ii ]
            end  # for jj in 1:nrows
        end  # for ii in 1:ncols

        excelGenerationTime = ( now() - tStart ).value / 1000.0
        ws[ "B3" ] = excelGenerationTime
    end  # XLSX.openxlsx( ... ) do xf

end  # excelSubpopReport( mpSim, timeGrid, nodes, filename, overwrite,
     #   timeFactor )

excelSubpopReport( mpSim::MPsim, timeRes::Real,
    subpops::Subpopulation...; filename::AbstractString = "compReport",
    overwrite::Bool = true, timeFactor::Real = 12.0 ) =
    excelSubpopReport( mpSim, MP.generateTimeGrid( mpSim, timeRes ),
    subpops...; filename = filename, overwrite = overwrite,
    timeFactor = timeFactor )


text:
abstract type MeanPrior{T} end

import Base: +, -, *, convert

include("constantmean.jl")
include("zeromean.jl")
include("empiricalmean.jl")


text:
@doc raw"""
    PoissonEquation(h, ensemble) ::GainEquation

Returns a Poisson equation struct  representing the equation ``\nabla\cdot(p\nabla \phi) = - \tilde h``, where ``p`` is a probability density and ``\tilde h = h-\int h p dx``.
The container contains the following fields
* `:h': ``h`` itself
* `:positions': an i.i.d. sample from ``p``, represented as a matrix
* `:H': the evaluation of ``h`` at the sample points
* `:mean_H': the sample average of `H'
* `:potential': the evaluation of ``\phi`` at the sample points
* `:gain': the evaluation of ``K=\nabla \phi`` at the sample points

    solve!(eq::PoissonEquation, method::GainEstimationMethod)

Fills the field `:gain' with appropriate values.
The fields `:H', `:mean_H', and `:potential' are stored to be re-used.

    update!(eq::PoissonEquation, ensemble)

Fills the fields `:positions', `:H', and `:mean_H' according to the new samples from `ensemble'.

    update!(eq::PoissonEquation)

Updates fields `:H', and `:mean_H' to be consistent with `:positions'.
"""
struct PoissonEquation{T, TF} <: GainEquation
    h::TF
    positions::Matrix{T}
    H::Matrix{T}
    mean_H::Matrix{T}
    potential::Matrix{T}
    gain::Array{T,3}

    PoissonEquation(h, positions, H, mean_H, potential, gain) =
        if  size(positions, 2) == size(H, 2) == size(potential, 2) == size(gain, 2) &&
            size(H, 1) == size(mean_H, 1) == size(potential, 1) == size(gain, 3) &&
            size(positions, 1) == size(gain, 1) &&
            size(mean_H,2) == 1
            new{eltype(positions), typeof(h)}(h, positions, H, mean_H, potential, gain)
        else
            throw(DimensionMismatch("Inconsistent dimensions of positions, H, mean_H, potential, or gain."))
        end
end

function PoissonEquation(h::Function, pos::AbstractMatrix)
    H = mapslices(h, pos, dims=1)
    mean_H = Statistics.mean(H, dims=2)
    return PoissonEquation(h, pos, H, mean_H, zeros(eltype(pos), size(H, 1), size(pos, 2)), zeros(eltype(pos), size(pos, 1), size(pos, 2), size(H, 1)))
end

function PoissonEquation(h::Function, ens::UnweightedParticleEnsemble)
    pos = copy(ens.positions)
    return PoissonEquation(h, pos)
end

state_dim(eq::PoissonEquation)       = size(eq.positions, 1)
obs_dim(eq::PoissonEquation)         = size(eq.H, 1)
no_of_particles(eq::PoissonEquation) = size(eq.positions, 2)

Base.show(io::IO, eq::PoissonEquation) = print(io, "Poisson equation for the gain
    # of particles:        ", size(eq.positions, 2),"
    hidden dimension:      ", state_dim(eq),"
    observed dimension:    ", obs_dim(eq))

function Htilde(eq::PoissonEquation)
    return eq.H .- eq.mean_H
end

function update!(eq::PoissonEquation)
    eq.H         .= mapslices(eq.h, eq.positions, dims=1)
    eq.mean_H    .= Statistics.mean(eq.H, dims=2)
    return eq
end

function update!(eq::PoissonEquation, pos::AbstractMatrix)
    eq.positions .= pos
    update!(eq)
end

function update!(eq::PoissonEquation, ens::UnweightedParticleEnsemble)
    update!(eq, ens.positions)
end

function GainEquation(state_model::DiffusionStateModel, obs_model::DiffusionObservationModel, ens::UnweightedParticleEnsemble)
    return PoissonEquation(obs_model.observation_function, ens)
end

function GainEquation(state_model::DiffusionStateModel, obs_model::DiffusionObservationModel, N::Int)
    ens = UnweightedParticleEnsemble(state_model, N)
    return GainEquation(state_model, obs_model, ens)
end

function GainEquation(filt_prob::AbstractFilteringProblem, ens::UnweightedParticleEnsemble)
    return GainEquation(state_model(filt_prob), obs_model(filt_prob), ens)
end

function GainEquation(filt_prob::AbstractFilteringProblem, N::Int)
    return GainEquation(state_model(filt_prob), obs_model(filt_prob), N)
end


text:
"PDDL state description."
mutable struct State
    types::Set{Term} # Object type declarations
    facts::Set{Term} # Boolean-valued fluents
    fluents::Dict{Symbol,Any} # All other fluents
end

"PDDL action description."
struct Action
    name::Symbol # Name of action
    args::Vector{Var} # Action parameters
    types::Vector{Symbol} # Parameter types
    precond::Term # Precondition of action
    effect::Term # Effect of action
end

Action(term::Term, precond::Term, effect::Term) =
    Action(term.name, get_args(term), Symbol[], precond, effect)

Base.:(==)(a1::Action, a2::Action) = (a1.name == a2.name &&
    Set(a1.args) == Set(a2.args) && Set(a1.types) == Set(a2.types) &&
    a1.precond == a2.precond && a1.effect == a2.effect)

"PDDL event description."
struct Event
    name::Symbol # Name of event
    precond::Term # Precondition / trigger of event
    effect::Term # Effect of event
end

Base.:(==)(e1::Event, e2::Event) = (e1.name == e2.name &&
    e1.precond == e2.precond && e1.effect == e2.effect)

"PDDL planning domain with events and axioms."
@kwdef mutable struct Domain
    name::Symbol # Name of domain
    requirements::Dict{Symbol,Bool} = Dict() # PDDL requirements used
    types::Dict{Symbol,Vector{Symbol}} = Dict() # Types and their subtypes
    constants::Vector{Const} = [] # List of constants
    constypes::Dict{Const,Symbol} = Dict() # Types of constants
    predicates::Dict{Symbol,Term} = Dict() # Dictionary of predicates
    predtypes::Dict{Symbol,Vector{Symbol}} = Dict() # Predicate type signatures
    functions::Dict{Symbol,Term} = Dict() # Dictionary of function declarations
    functypes::Dict{Symbol,Vector{Symbol}} = Dict() # Function type signatures
    funcdefs::Dict{Symbol,Any} = Dict() # Dictionary of function definitions
    axioms::Vector{Clause} = [] # Axioms / derived predicates
    actions::Dict{Symbol,Action} = Dict() # Action definitions
    events::Vector{Event} = [] # Event definitions
    _extras::Dict{Symbol,Any} # Extra fields
end

function Domain(name::Symbol, header::Dict{Symbol,Any}, body::Dict{Symbol,Any})
    h_extras = filter(item -> !(first(item) in fieldnames(Domain)), header)
    b_extras = filter(item -> !(first(item) in fieldnames(Domain)), body)
    extras = merge!(h_extras, b_extras)
    header = filter(item -> first(item) in fieldnames(Domain), header)
    axioms = Clause[get(body, :axioms, []); get(body, :deriveds, [])]
    body = filter(item -> first(item) in fieldnames(Domain), body)
    body[:axioms] = axioms
    body[:actions] = Dict(act.name => act for act in body[:actions])
    return Domain(;name=name, _extras=extras, header..., body...)
end

Base.copy(d::Domain) =
    Domain(; Dict(fn => getfield(d, fn) for fn in fieldnames(typeof(d)))...)

Base.getproperty(d::Domain, s::Symbol) =
    hasfield(Domain, s) ? getfield(d, s) : d._extras[s]

Base.setproperty!(d::Domain, s::Symbol, val) =
    hasfield(Domain, s) && s != :_extras ?
        setfield!(d, s, val) : setindex!(d._extras, val, s)

Base.propertynames(d::Domain, private=false) = private ?
    tuple(fieldnames(Domain)..., keys(d._extras)...) :
    tuple(filter(f->f != :_extras, fieldnames(Domain))..., keys(d._extras)...)

"Get domain constant type declarations as a set of facts."
function get_const_facts(domain::Domain)
  return Set([@julog($ty(:o)) for (o, ty) in domain.constypes])
end

"Get domain constant type declarations as a list of clauses."
function get_const_clauses(domain::Domain)
   return [@julog($ty(:o) <<= true) for (o, ty) in domain.constypes]
end

"Get domain type hierarchy as a list of clauses."
function get_type_clauses(domain::Domain)
    clauses = [[Clause(@julog($ty(X)), Term[@julog($s(X))]) for s in subtys]
               for (ty, subtys) in domain.types if length(subtys) > 0]
    return length(clauses) > 0 ? reduce(vcat, clauses) : Clause[]
end

"Get all proof-relevant Horn clauses for PDDL domain."
function get_clauses(domain::Domain)
   return [domain.axioms; get_const_clauses(domain); get_type_clauses(domain)]
end

"Get list of predicates that are never modified by actions in the domain."
function get_static_predicates(domain::Domain)
    ground = t ->
        substitute(t, Subst(v => Const(gensym()) for v in Julog.get_vars(t)))
    diffs = [effect_diff(ground(act.effect)) for act in values(domain.actions)]
    derived = p -> any(unify(p, ax.head) !== nothing for ax in domain.axioms)
    modified = p -> any(contains_term(d, p) for d in diffs)
    return Term[p for p in values(domain.predicates)
                if !derived(p) && !modified(p)]
end

"Get list of functions that are never modified by actions in the domain."
function get_static_functions(domain::Domain)
    ground = t ->
        substitute(t, Subst(v => Const(gensym()) for v in Julog.get_vars(t)))
    diffs = [effect_diff(ground(act.effect)) for act in values(domain.actions)]
    modified = p -> any(contains_term(d, p) for d in diffs)
    return Term[p for p in values(domain.functions) if !modified(p)]
end

"PDDL planning problem."
@kwdef mutable struct Problem
    name::Symbol # Name of problem
    domain::Symbol # Name of associated domain
    objects::Vector{Const} # List of objects
    objtypes::Dict{Const,Symbol} # Types of objects
    init::Vector{Term} # Predicates that hold in initial state
    goal::Term # Goal formula
    metric::Union{Tuple{Int,Term},Nothing} # Metric direction (+/-) and formula
end

function Problem(name::Symbol, header::Dict{Symbol,Any}, body::Dict{Symbol,Any})
    header = filter(item -> first(item) in fieldnames(Problem), header)
    body = filter(item -> first(item) in fieldnames(Problem), body)
    return Problem(;name=name, header..., body...)
end

function Problem(state::State, goal::Term=@julog(and()),
                 metric::Union{Tuple{Int,Term},Nothing}=nothing;
                 name=:problem, domain=:domain)
    objtypes = Dict{Const,Symbol}(get_args(t)[1] => t.name for t in state.types)
    objects = collect(keys(objtypes))
    init = Term[get_facts(state); get_assignments(state)]
    return Problem(Symbol(name), Symbol(domain),
                   objects, objtypes, init, goal, metric)
end

Base.copy(p::Problem) =
    Problem(; Dict(fn => getfield(p, fn) for fn in fieldnames(typeof(p)))...)

"Get object type declarations as a list of clauses."
function get_obj_clauses(problem::Problem)
    return [@julog($ty(:o) <<= true) for (o, ty) in problem.objtypes]
end


text:
module FindDefinition

export finddef, finddefs

import MacroTools

ismacrocall(_) = false
ismacrocall(expr::Expr) = expr.head == :macrocall

ismoduledef(_) = false
ismoduledef(e::Expr) = e.head == :module

isinclude(_) = false
isinclude(e::Expr) = e.head == :call && e.args[1] == :include

isquote(_) = false
isquote(::QuoteNode) = true
isquote(e::Expr) = e.head == :quote

args(_) = []
args(e::Expr) = e.args

function collect_nodes(f, _module, ast; skip_quoted = true)
    skip(arg) = skip_quoted && isquote(arg)
    init = f(ast) ? [(from_module = _module, ex = ast)] : []
    if ismoduledef(ast)
        modulename = ast.args[2]
        _module = getproperty(_module, modulename)
    end
    reduce( ∪, collect_nodes(f, _module, arg) for arg in args(ast) if !skip(arg); init )
end

collect_includefiles(_module, ast) = [(; from_module = mod, file = ex.args[2])
                                            for (mod, ex) in collect_nodes(isinclude, _module, ast)]

function parsefile(path, fromline = 1; greedy=true)
    str = read(path, String)
    newlines = only.(findall("\n", str))
    res = Expr[]
    pos = fromline == 1 ? 1 : newlines[fromline-1]+1
    while pos < length(str)
        lineshift = searchsortedfirst(newlines, pos) - 1
        fix_lnn(e) = e
        fix_lnn(lnn::LineNumberNode) = LineNumberNode(lnn.line + lineshift, path)

        expr, pos = Meta.parse(str, pos)
        isnothing(expr) && continue ## why can this happen?
        expr = MacroTools.postwalk(fix_lnn, expr)
        push!(res, expr)
        greedy || break
    end
    res
end

function rec_find_expr(f, _module, ast, includepath)
    res = collect_nodes(f, _module, ast)
    includes = collect_includefiles(_module, ast)
    for (mod, inc) in includes
        if !(inc isa AbstractString)
            @warn "Cannot resolve non-literal include `include($inc)` from module $mod. Skipped."
            continue
        end
        inc = joinpath(includepath, inc)
        if !isfile(inc)
            @warn "Included file not found: $inc."
            continue
        end
        for expr in parsefile(inc)
            append!(res, rec_find_expr(f, mod, expr, dirname(inc)))
        end
    end
    res
end

trysplitdef(ex) = try
    MacroTools.splitdef(ex)
catch
    nothing
end
# MacroTools.isdef is broken. See MacroTools issue #154
isfunctiondef(ex) = trysplitdef(ex) != nothing

istypeassertion(ex) = false
istypeassertion(ex::Expr) = ex.head == :(::)

isinterpolated(ex) = false
isinterpolated(ex::Expr) = ex.head == :$

dequalify(name) = name
dequalify(name::Expr) = (name.head == :(.) || error("Not a qualified name: $name");
                         name.args[2].value)
dequalify(name::GlobalRef) = name.name
argtypes(method) = Base.tail(tuple(method.sig.types...))
function defines_this_method(_module, ex, method)
    d = trysplitdef(ex)
    isnothing(d) && return false
    qualified_name = get(d,:name,nothing)
    qualified_name === nothing && return false
    istypeassertion(qualified_name) && return false # a type method definition, e.g. (foo::Foo)(x,y) = ...
    if isinterpolated(qualified_name)
        @warn "Cannot resolve interpolated function name $qualified_name in function definition $ex"
        return false
    end
    dequalify(qualified_name) != method.name && return false
    # TODO should be possible without `eval`, using `Meta.lower`
    newname = gensym()
    d[:name] = newname
    newexp = MacroTools.combinedef(d)
    f = try
        _module.eval(newexp)
    catch e
        @warn("Failed to evaluate symgen-ed function definition in module $_module.",
                definition = newexp)
        showerror(stderr, e, catch_backtrace())
        println(stderr)
        return false
    end
    exmethod = only(methods(f))
    return argtypes(exmethod) == argtypes(method) &&
           Base.kwarg_decl(exmethod) == Base.kwarg_decl(method)
end

macrocall_linenumnode(macrocall) = macrocall.args[2]
function find_definitions(method::Method)
    _module = method.module
    eval_method = only(methods(_module.eval))
    file, line = string(eval_method.file), eval_method.line
    module_def = only(parsefile(file, line; greedy=false))
    res = []
    mcls = rec_find_expr(ismacrocall, _module, module_def, dirname(file))
    for (_module, macrocall) in mcls
        lnn = macrocall_linenumnode(macrocall)
        expanded = try
            macroexpand(_module, macrocall)
        catch e
            @warn("Failed to expand macrocall in $_module:$lnn. Skipping.",
                    macrocall=macrocall)
            showerror(stderr, e, catch_backtrace())
            println(stderr)
            continue
        end
        funcdefs = collect_nodes(isfunctiondef, _module, expanded)
        for (_module, funcdef) in funcdefs
            if defines_this_method(_module, funcdef, method)
                push!(res, (; lnn, funcdef))
            end
        end
    end
    res
end

find_definitions(f::Function) = reduce( ∪, find_definitions(m) for m in methods(f) )

finddef(m::Method) = last(find_definitions(m)).lnn
finddefs(f::Function) = map(finddef, methods(f))

end


text:
module CovSel

using ProximalBase
using DataStructures: BinaryMaxHeap
using CoordinateDescent, Distributions
using Printf, SparseArrays, LinearAlgebra

export
  ADMMOptions,
  covsel!, covselpath,
  #covsel_refit!, covselpath_refit

  # difference estimation
  #
  fusedGraphicalLasso, fusedGraphicalLasso!,   # Danaher et al.
  fusedNeighborhoodSelection,

  # Direct Difference Estimation and Inference
  differencePrecisionActiveShooting, differencePrecisionActiveShooting!, differencePrecisionRefit,
  CDDirectDifferenceLoss,

  differencePrecisionIHT,
  differenceLatentPrecisionIHT

include("utils.jl")

include("admmCovSel.jl")
include("tvCovSel.jl")
include("LatentCovSel_IHT.jl")

include("directDifference.jl")
include("fusedDifference.jl")

include("generate.jl")
include("evaluation.jl")

# alternative implementations
include("altImplementation.jl")

end


text:
using Documenter, Shapes

makedocs(;
    modules=[Shapes],
    format=Documenter.HTML(
        prettyurls = get(ENV, "GITHUB_ACTIONS", nothing) == "true",
        canonical = "https://lyceum.github.io/Shapes.jl/stable/"
    ),
    pages=[
        "Home" => "index.md",
        "API" => "api.md",
    ],
    sitename="Shapes.jl",
    authors = "Colin Summers",
    clean = true,
    doctest = true,
    checkdocs = :exports,
    linkcheck = :true,
)

deploydocs(
    repo = "github.com/Lyceum/Shapes.jl.git",
    push_preview = true
)


text:
using Distributed
addprocs(2)

@everywhere using Pkg
@everywhere Pkg.activate("./julia-multiprocessing/TestModule1")

@everywhere using TestModule1

#%%
print_nprocs()

func1()

func2()

test_interp(3.1415)
sin(3.1415)

TestModule1.create_interp()

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples