text:
using Mustache, HttpCommon
using StatsBase: RealVector, fweights
import Statistics: quantile, mean

# used tableau 10 colormap
HTML_COLORS = ["#4E79A7","#F28E2C","#E15759","#76B7B2","#59A14F", "#EDC949","#AF7AA1","#FF9DA7","#9C755F","#BAB0AB"]

quantile(v::Dict{Int, Int}, p::Number) = quantile(v |> keys |> collect, v |> values |> collect |> fweights, [p])[1]
quantile(v::Dict{Int, Int}, p::RealVector) = quantile(v |> keys |> collect, v |> values |> collect |> fweights, p)
mean(v::Dict{Int, Int}) = mean(v |> keys |> collect, v |> values |> collect |> fweights)

function schema2html(e::Entry; pad = "", max_vals=100, parent_updated=nothing, parent_key="")
	c = HTML_COLORS[((length(pad)÷2)%length(COLORS))+1]
	filled_percent = isnothing(parent_updated) ? "" : ", filled = $(10000 * e.updated ÷ parent_updated / 100)%"
	sorted_counts = sort(collect(e.counts), by=x->x[2], rev=true)
	ret_str = """
$pad
$pad[Scalar - $(join(types(e)))], $(length(keys(e.counts))) unique values,
(updated = $(e.updated)$filled_percent, min=$(Tuple(sorted_counts[end])), max=$(Tuple(sorted_counts[1])))
"""
	i = 0
    for (key, val) in sorted_counts
		key_repr = key isa String ? escapeHTML(key) : "$key"
		ret_str *= pad*" "^2 * "
$key_repr: $val
\n"
		i += 1
		if i == max_vals
			ret_str *= pad*" "^2 * "
and other $(length(e.counts) - i) values
\n"
			break
		end
    end
	ret_str * "$pad
\n"
end

# because sometimes there is empty list in all jsons, this helps to determine the pruning of such element
function schema2html(e::Nothing; pad = "", max_vals=100, parent_updated=nothing, parent_key="")
	c = HTML_COLORS[((length(pad)÷2)%length(COLORS))+1]
	ret_str = """
$pad[Empty list element], this list is empty in all JSONs, can not infer schema, suggesting to delete key $parent_key
"""
	ret_str
end

function schema2html(e::ArrayEntry; pad = "", max_vals=100, parent_updated=nothing, parent_key="")
 	c = HTML_COLORS[((length(pad)÷2)%length(COLORS))+1]
	# todo: fix it all so it is different method for array of entries and the rest so only nested things are truly nested
	filled_percent = isnothing(parent_updated) ? "" : ", filled=$(10000 * e.updated ÷ parent_updated / 100)%"
	quantiles = quantile(e.l, [0.1, 0.5, 0.9])
	min_val = minimum(keys(e.l))
	max_val = maximum(keys(e.l))
	mean_val = round(mean(e.l), digits=2)
	q10_val = round(quantiles[1], digits=2)
	median_val = round(quantiles[2], digits=2)
	q90_val = round(quantiles[3], digits=2)
	ret_str = """
$pad
[List] (updated=$(e.updated)$filled_percent, mean=$mean_val,
min=$min_val, max=$max_val, 10th percentile=$q10_val, median=$median_val, 90th percentile=$q90_val)
$pad
with following frequencies
$pad

"""
 	i = 0
 	for (key, val) in sort(collect(e.l))
		ret_str *= pad*" "^2 * "
$key: $val
\n"
		i += 1
		if i == max_vals
			ret_str *= pad*" "^2 * "
and other $(length(e.l) - i) values
\n"
			break
		end
 	end
	ret_str *= """
$pad

$pad

$pad
and data
"""
	ret_str *= schema2html(e.items, pad=pad*" "^2, max_vals=max_vals, parent_key="$parent_key.items")
	ret_str * """
$pad

$pad

"""
end

function schema2html(e::DictEntry; pad = "", max_vals=100, parent_updated=nothing, parent_key="")
	c = HTML_COLORS[((length(pad)÷2)%length(COLORS))+1]
    if isempty(e.childs)
    	return pad * """
Empty Dict
\n"""
    end
	class = isempty(pad) ? "top_dict" : "nested"
	filled_percent = isnothing(parent_updated) ? "" : ", filled=$(10000 * e.updated ÷ parent_updated / 100)%"
	ret_str = pad * """
[Dict] (updated=$(e.updated)$filled_percent)\n"""
	i = 0
    for (key, val) in sort(e.childs)
		child_key = """$parent_key["$key"]"""
		ret_str *= pad*" "^2 * """
$key - $child_key<input type="checkbox" name="$child_key" value="$child_key">\n"""
		ret_str *= schema2html(val, pad=pad*" "^4, max_vals=max_vals, parent_updated=e.updated, parent_key=child_key)
		ret_str *= pad*" "^2 * "
\n"
		i += 1
		if i == max_vals
			ret_str *= pad*" "^2 * "
and other $(length(e.childs) - i) values
\n"
			break
		end
    end
	ret_str * pad * "
\n"
end

function generate_html(sch::DictEntry; max_vals=100)
	tpl = mt"""
	<!DOCTYPE html>

    	<meta charset="UTF-8">
    	<title>Json schema dump</title>

		    <button type="button" id="copy_clipboard">Copy selectors to clipboard</button>

		{{{list_dump}}}

"""

	d = Dict(
		"list_dump" => schema2html(sch, max_vals=max_vals),
	)
	return(render(tpl, d))
end

function generate_html(sch::DictEntry, file_name ; max_vals=100)
	s = generate_html(sch; max_vals = max_vals)
	open(file_name, "w") do f
 		write(f, s)
	end
end

# Base.show(io, m::MIME{Symbol("text/html")}, sch::JsonGrinder.DictEntry) = print(io, repr(m,sch))

Base.repr(::MIME"text/html", sch::DictEntry; max_vals = 100, context = nothing) = generate_html(sch; max_vals = max_vals)
Base.showable(::MIME"text/html", ::DictEntry) = true


text:
using Quantiles
using Test

@testset "Quantiles.jl" begin
    # Write your own tests here.
end

ext: jl


text:
function initclimate(p::ClimateParams, firstyear, usecache)
	if usecache && firstyear > YEARS[1]
		coeffs = cached_coeff_state[:,:,div(firstyear-1800, 10) + 1]
		s = ClimateState([interpolatespline(p.lambda, coeffs[i,:]) for i=1:75])
	else
		s = ClimateState()
		init_temperatures!(s)
		init_oceancarbon!(s, p)
		init_biocarbon!(s, p)
		init_concentration_and_forcing!(s, p)
	end
	return s
end

# firstyear & lastyear: years to start & end the model run (max 1765-2500)
# usecache: whether to use the precalculated cache of results 1765-2010 to speed up run times
# timestep [years]: time step for solving the climate differential equations
# lambda [°C/(W/m^2)]: climate sensitivity, related to CS per doubling of CO2 by CS [°C/(2xCO2)] = 3.7 W/m^2 * lambda
# oceantempfeedback: set to 0 to disable temperature feedback in ocean carbon uptake (default is 1)
# bioQ10factor: set to 1 to disable temperature feedback in the terrestrial biosphere (default is 2)
# equilibriumCO2 [ppm]: use 278 to start model in preindustrial times
function getparams(annualemissions;  firstyear::Int=2010, lastyear::Int=2100, usecache::Bool=true,
						timestep::Float64=0.01, lambda::Float64=0.8, rcp::String="RCP45",
						oceantempfeedback::Float64=1.0, bioQ10factor::Float64=2.0, equilibriumCO2::Float64=278.0)

	p = ClimateParams(timestep, lambda, 0.0, 0.0, oceantempfeedback, bioQ10factor, equilibriumCO2)
	if usecache && oceantempfeedback == 1.0 && bioQ10factor == 2.0 && equilibriumCO2 == 278.0
		p.aerosolforcingfactor = interpolatespline(lambda, cached_coeff_forcing)
		p.fertilization = interpolatespline(lambda, cached_coeff_fertilization)
	else
		calibrateforcing!(p, rcp)
		calibratefertilization!(annualemissions, p, rcp)
	end
	return p, firstyear, lastyear, usecache, rcp
end

# scen: BAU, 3DEG, 2DEG, 1DEG, 15DEG, 2DEGB, 15DEGB, 1PCT, 1PCTD, 3GTC, 2C-3CS, 2C-45CS, 2C-15CS
function runscen(scen::String;  kwargs...)
	results, p, annualemissions, rcp = solveclimate(scen;  kwargs...)
	printresults(2010:10:2100, results, p, annualemissions, rcp)
end

solveclimate(; kwargs...) = solveclimate("2DEG";  kwargs...)
function solveclimate(scen::String;  kwargs...)
	annualemissions = getscenario(scen)
	p, firstyear, lastyear, usecache, rcp = getparams(annualemissions; kwargs...)
	rcp = scen[1:3] == "RCP" ? scen : rcp
	results = solveclimate(annualemissions, p, firstyear, lastyear, usecache, rcp)
	return results, p, annualemissions, rcp
end

function solveclimate(annualemissions;  kwargs...)
	p, firstyear, lastyear, usecache, rcp = getparams(annualemissions; kwargs...)
	results = solveclimate(annualemissions, p, firstyear, lastyear, usecache, rcp)
	return results, p, rcp
end

function solveclimate(annualemissions, p::ClimateParams,
						firstyear::Int=2010, lastyear::Int=2100, usecache::Bool=true, rcp::String="RCP45")
	stabilitycheck(p.timestep)

	firstyear = usecache ? firstyear : YEARS[1]
	roundedyear = firstyear >= 1800 ?  min(2010, 10*floor(Int, firstyear/10)) :  1765
	s = initclimate(p, roundedyear, usecache)

	results = Vector{ClimateState}(undef, iyear(lastyear))

	# don't run the first time step if using the cache (cached results are taken after the iteration)
	startyear = float(roundedyear)
	if usecache
		startyear = roundedyear + p.timestep
		results[iyear(roundedyear)] = copy(s)
	end

	# main time loop of the model
	for t::Float64 in startyear:p.timestep:lastyear
		# the shorter version is type unstable because of the closure bug
		#emissions = Dict(g => interpolate(t, annualemissions[g]) for g in GAS3)
		emissions = Dict(	:CO2 => interpolate(t, annualemissions[:CO2]),
							:CH4 => interpolate(t, annualemissions[:CH4]),
							:N2O => interpolate(t, annualemissions[:N2O]))

		temperatures!(s, p)		# accuracy is better if temperatures are calculated first
		carbonbalance!(emissions, s, p)
		oceancarbon!(s, p)
		biocarbon!(s, p)
		concentrations!(emissions, s, p)
		radiativeforcing!(t, s, p, rcp)

		if t ≈ round(t)
			results[iyear(t)] = copy(s)
		end
	end
	return results
end

function printresults(yearrange, res, p::ClimateParams, annualemissions, rcp::String="RCP45")
	ly = length(yearrange)
	out_carbon = Matrix(undef, ly,3)
	out_emissions = Matrix(undef, ly,3)
	out_concentration = Matrix(undef, ly,3)
	out_temp_forcing = Matrix(undef, ly,6)
	out_temp_ocean = Matrix(undef, ly,5)

	row = 0
	for i in iyear.(yearrange)
		row += 1
		RF_nonCO2 = otherforcing[rcp][i] + p.aerosolforcingfactor * aerosolforcing[rcp][i] +
			sum(res[i].RadiativeForcing[g] for g in GAS if g != :CO2)
		out_carbon[row,:] =
			[res[i].Concentration[:CO2]/ppm_per_GtC  sum(res[i].BioReservoir)  res[i].DeltaDIC/micromol_per_kg_per_GtC]
		out_emissions[row,:] =
			[annualemissions[:CO2][i]  annualemissions[:CH4][i]  annualemissions[:N2O][i]]
		out_concentration[row,:] =
			[res[i].Concentration[:CO2]  res[i].Concentration[:CH4]  res[i].Concentration[:N2O]]
		out_temp_forcing[row,:] =
			[res[i].Temp_global  res[i].TotalRadiativeForcing*p.lambda  res[i].Temp_land  res[i].TotalRadiativeForcing  RF_nonCO2  res[i].Temp_land./res[i].Temp[1]]
		out_temp_ocean[row,:] =
			[res[i].Temp[1]  res[i].Temp[2]  res[i].Temp[5]  res[i].Temp[10]  res[i].Temp[maxlayers]]
	end

	#println("\nAtmosphereCarbon BioCarbon OceanCarbon")
	#display([yearrange out_carbon])
	println("\nEmissions_CO2 Emissions_CH4 Emissions_N2O")
	display([yearrange out_emissions])
	println("\nConcentration_CO2 Concentration_CH4 Concentration_N2O")
	display([yearrange out_concentration])
	println("\nTemp_global EquilibriumTemperature Temp_land TotalRadiativeForcing TotalNonCO2forcing R_lambda")
	display([yearrange out_temp_forcing])
	println("\nTemp_oceanlayer (layers 1,2,5,10,max)")
	display([yearrange out_temp_ocean])

	println("\naerosolforcingfactor & fertilization: ", p.aerosolforcingfactor, "  ", p.fertilization)
end

max_stars_repo_name: niclasmattsson/ClimateCalculator

max_issues_repo_name: niclasmattsson/ClimateCalculator

max_forks_repo_name: niclasmattsson/ClimateCalculator


text:
function XC_c_pw( Rhoe )
    # J.P. Perdew and Y. Wang, PRB 45, 13244 (1992)

    third = 1.0/3.0
    pi34 = 0.6203504908994
    rs = pi34/Rhoe^third

    a = 0.031091
    a1 = 0.21370
    b1 = 7.5957
    b2 = 3.5876
    b3 = 1.6382
    b4 = 0.49294

    # interpolation formula
    rs12 = sqrt(rs)
    rs32 = rs * rs12
    rs2 = rs^2

    om = 2.0 * a * (b1 * rs12 + b2 * rs + b3 * rs32 + b4 * rs2)
    dom = 2.0 * a * (0.5 * b1 * rs12 + b2 * rs + 1.5 * b3 * rs32 + 2.0 * b4 * rs2)
    olog = log(1.0 + 1.0/om)
    ec = -2.0 * a * (1.0 + a1 * rs) * olog
    vc = -2.0*a*(1.0 + 2.0/3.0 * a1 * rs) * olog - 2.0/3.0 * a * (1.0 + a1*rs) * dom/ (om * (om + 1.0) )

    return ec, vc
end


text:
include("ReadDataFile.jl")

module SetProb

using Main.ReadDataFile

export
    CParam, read_prob_data

mutable struct CParam
    rho   ::Float64
    bulk  ::Float64
    cc    ::Float64
    zz    ::Float64

    beta  ::Float64

    # default constructor
    function CParam()
        new(0.0, 0.0, 0.0, 0.0, 0.0)
    end
end

# Set the material parameters for the acoustic equations

function read_prob_data(parms::CParam)
    r = Reader("setprob.data")

    # density:
    parms.rho  = read_flt(r)

    # bulk modulus:
    parms.bulk = read_flt(r)

    # sound speed:
    parms.cc = sqrt(parms.bulk/parms.rho)

    # impedance:
    parms.zz = parms.cc*parms.rho

    # beta for initial conditions:
    parms.beta = read_flt(r)

    close_reader(r)
end # read_prob_data

end # SetProb

# using .SetProb
# using Base.Test.@test
#
# function test_read_data()
#     const parms = CParam()
#     read_prob_data(parms)
#     println("CParam = $parms")
#     return true
# end

# @test test_read_data()


text:
using Knet,CUDArt

include("../src/op/conv.jl")
include("../src/op/pool.jl")

x = CudaArray(rand(6,5,4,3))
w = CudaArray(rand(3,3,4,2))
c = conv()
s = Knet.infersize(c, size(w), size(x))
y = CudaArray(Float64, s[end])
Knet.forw(c, w, x, y)
dy = similar(y); rand!(dy)
dw = similar(w)
dx = similar(x)
Knet.back(c, dy, dw, dx; x=(w,x))

p = pool()
s = Knet.infersize(p, size(x))
y = CudaArray(Float64, s[end])
Knet.forw(p, x, y)
dy = similar(y); rand!(dy)
dx = similar(x)
Knet.back(p, dy, dx; x=x, y=y)


text:
using FlightGNC
using UnPack, ComponentArrays, Transducers
using LinearAlgebra, DifferentialEquations  # for callbacks

function main(p_M_0::Vector, v_M_0::Vector, p_T_0::Vector, v_T_0::Vector, s_guidance::BPNG; A_M_max = Inf)
    @unpack dim = s_guidance
    # Environments
    pursuer = PointMassMissile(dim, A_M_max)
    evador = PointMassMissile(dim, Inf)
    env = PursuerEvadorMissile(pursuer, evador)

    # Initial condition
    x0_pursuer = State(pursuer)(p_M_0, v_M_0)
    x0_evador = State(evador)(p_T_0, v_T_0)
    x0 = State(env)(x0_pursuer, x0_evador)

    # Simulation parameters
    Δt      = 0.01
    t_sim_f = 50

    # callbacks
    function condition_stop(u, t, integrator)
        (p_M, v_M, p_T, v_T) = (u.pursuer.p, u.pursuer.v, u.evador.p , u.evador.v)
        r = norm(p_T-p_M)
        ṙ = dot(p_T-p_M, v_T-v_M) / r
        r < 0.1  || (r < 10 && ṙ >= 0)
    end
    affect!(integrator) = terminate!(integrator)  # See DiffEq.jl documentation
    cb_stop    = DiscreteCallback(condition_stop, affect!)
    cb = CallbackSet(cb_stop)  # useful for multiple callbacks

    # Execute Simulation
    # prob: DE problem, df: DataFrame
    @time prob, df = FSimBase._sim(
                         x0,  # initial condition
                         apply_inputs(Dynamics!(env);
                                      u_pursuer=BPNG_cmd(s_guidance),
                                      u_evador=(x, params, t) -> zeros(evador.dim));  # dynamics!; apply_inputs is exported from FS and is so useful for systems with inputs
                         tf=t_sim_f,
                         savestep=Δt,  # savestep is NOT simulation step
                         solver=Tsit5(),
                         callback=cb,
                         reltol=1e-6
                        )  # sim is exported from FS

	return df
end


text:
# Here are leftovers from  my previous attempts to run the EWS from Julia, including programming the gaussian filter rolling windows.

#Program the gaussian filter as a function and run it before calculating the EWS. Note here that the \sigma is the standard deviation of complete time series. The following code doesn't work!
# import Statistics: var, mean, std
# function gaussian_filter(xin, xout)
#     sd = std(skipmissing(xin))

#     if sd > 0 # std non-zero
#         xout[:].=e^(.-xin^2 ./(2*sd^2))./ (sd*sqrt(2*pi))
#     else # time series is probably constant / sd == 0, so avoiding division by zero
#         xout[:].=0
#     end
# end

## Define input and output dimensions
indims = InDims("Time")
outdims = OutDims("Time")

# test the function
test = mapCube(
    gaussian_filter,
    subsetcube(d2, Scale = "Fast Oscillations"),
    indims=indims, outdims=outdims)

    # Not sure this is needed. The TS are normalized to 0-mean and 1-variance
# import Statistics: var, mean, std
# std_dev = mapslices(
#     std ∘ skipmissing,
#     subsetcube(d2, Scale = "Fast Oscillations"),
#     dims = "Time")

max_stars_repo_name: juanrocha/ESDL

max_issues_repo_name: juanrocha/ESDL

max_forks_repo_name: juanrocha/ESDL


text:
# Copyright 2019 Digital Domain 3.0
#
# Licensed under the Apache License, Version 2.0 (the "Apache License")
# with the following modification; you may not use this file except in
# compliance with the Apache License and the following modification to it:
# Section 6. Trademarks. is deleted and replaced with:
#
# 6. Trademarks. This License does not grant permission to use the trade
#    names, trademarks, service marks, or product names of the Licensor
#    and its affiliates, except as required to comply with Section 4(c) of
#    the License and to reproduce the content of the NOTICE file.
#
# You may obtain a copy of the Apache License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the Apache License with the above modification is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the Apache License for the specific
# language governing permissions and limitations under the Apache License.

using Test
using Handles

@HandleType Foo
@HandleType Bar
@HandleType Baz

hey(x::Integer) = x+1
hey(x::Foo) = string("foo", x)
hey(x::Bar) = string("bar", x)

@testset "handles" begin
  @test Foo(1) == 1

  a = sort(map( f->f+1 + Foo(1), (2 .* Foo.([10:-1:-2...]))))
  @test eltype(a) == Foo
  @test a[1] == -2 == Foo(-2)
  @test a[2] < a[3]
  @test a[end] == 22

  @test hey(Foo(2)) == "foo2"
  @test hey(Bar(3)) == "bar3"
  @test hey(Baz(4)) == 5

end


text:
# Early tests -- keep around for now.
module OldTests

using Docile
@docstrings

@doc {
    :file    => "methoddoc.md",
    :section => "Tests"
    } ->
function g(x)
end

@doc tex"""
The `tex` multiline string allows LaTeX input without escaping
characters.

$$
\frac{a + b}{a - b}
$$

Some more text.
""" {
    :key => :value
    } ->
h(x) = x

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut tristique
vitae, sagittis vel odio. Maecenas convallis ullamcorper ultricies.

```julia
a = 2
b = 3
a + b
```
""" {
    :key => :value
    } ->
function f(x::Int)
end

@doc """
Nam dictum, odio nec pretium volutpat, arcu ante placerat erat, non
tristique elit urna et turpis. Quisque mi metus, ornare sit amet
fermentum et, tincidunt et orci. Fusce eget orci a orci congue
vestibulum.

A code block that runs:

```julia
a = rand(10,10)
b = rand(10)
a * b
```

A code block that fails:

```julia
sqrt(-1)
```

""" {
    :parameters => {
        (:x, "the string argument")
        }
    } ->
f(x::String) = x

@doc """
Suspendisse lectus leo, consectetur in tempor sit amet, placerat quis
neque. Etiam luctus porttitor lorem, sed suscipit est rutrum non.
Curabitur lobortis nisl a enim congue semper. Aenean commodo ultrices
imperdiet. Vestibulum ut justo vel sapien venenatis tincidunt. Phasellus
eget dolor sit amet ipsum dapibus condimentum vitae quis lectus. Aliquam
ut massa in turpis dapibus convallis.

* one
* two
* three
* four

""" ->
function f(x::Float64)
end

@doc """
Nam dictum, odio nec pretium volutpat, arcu ante placerat erat, non
tristique elit urna et turpis.

```julia
a = rand(5,5)
@show @which f(a)
```
""" ->
f(x::Matrix) = x^2

## types ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––

@doc """
Suspendisse lectus leo, consectetur in tempor sit amet, placerat quis
neque. Etiam luctus porttitor lorem, sed suscipit est rutrum non.
Curabitur lobortis nisl a enim congue semper. Aenean commodo ultrices
imperdiet.
""" {
    :tags => {"one", "two", "three"}
    } ->
abstract FooAbs1

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut tristique
vitae, sagittis vel odio. Maecenas convallis ullamcorper ultricies.
Curabitur ornare, ligula semper consectetur sagittis, nisi diam iaculis
velit, id fringilla sem nunc vel mi.
""" ->
abstract FooAbs2

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
abstract FooAbs3{S}

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
abstract FooAbs4{S}

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
abstract Foo1 <: FooAbs1

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
abstract Foo2 <: FooAbs2

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
typealias Foo3 Int

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
typealias Foo4 Int

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
immutable Foo5
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
immutable Foo6
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
immutable Foo7 <: FooAbs1
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
immutable Foo8 <: FooAbs2
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
immutable Foo9{T}
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
immutable Foo10{T}
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
immutable Foo11{T} <: FooAbs1
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
immutable Foo12{T} <: FooAbs2
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
immutable Foo13{T,S} <: FooAbs3{S}
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
immutable Foo14{T,S} <: FooAbs4{S}
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
type Foo15
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
type Foo16
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
type Foo17 <: FooAbs1
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
type Foo18 <: FooAbs2
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
type Foo19{T}
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
type Foo20{T}
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
type Foo21{T} <: FooAbs1
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
type Foo22{T} <: FooAbs2
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
type Foo23{T,S} <: FooAbs3{S}
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
type Foo24{T,S} <: FooAbs4{S}
end

## macros –––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
macro foo(args)
end

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
macro foo(args)
end

## constants ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––

@doc """ Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
const FOO = "foo"

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
const FOO = "foo"

## globals ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
global FOO = "foo"

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
global FOO = "foo"

## global constants –––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
global const FOO = "foo"

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
global const FOO = "foo"

## non constants ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
FOO = "foo"

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
FOO = "foo"

## symbols (functions and modules) ––––––––––––––––––––––––––––––––––––––––––––––––––––––

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" {
    :key => :value
    } ->
f

@doc """
Vivamus fermentum semper porta. Nunc diam velit, adipiscing ut
tristique vitae, sagittis vel odio. Maecenas convallis ullamcorper
ultricies. Curabitur ornare, ligula semper consectetur sagittis, nisi
diam iaculis velit, id fringilla sem nunc vel mi.
""" ->
g

@doc """
Documentation for the Module itself.
""" {
    :tags => {"One", "Two", "Three"}
    } ->
OldTests

# fix for `@doc` with no docstring
@doc {
    :file => "methoddoc.md"
    } ->
function h(x, y, z)

end
@doc {
    :file => "methoddoc.md"
    } ->
type FooBar30

end

@doc {
    :file => "methoddoc.md"
    } ->
abstract FooBar31

@doc {
    :file => "methoddoc.md"
    } ->
global const FOOBAR1 = 4

@doc {
    :file => "methoddoc.md"
    } ->
const FOOBAR2 = 4

@doc {
    :file => "methoddoc.md"
    } ->
FOOBAR3 = 4

@doc {
    :file => "methoddoc.md"
    } ->
h(x, y, z) = x

doctest(OldTests)
doctest(OldTests; verbose = true)

end


text:
export getSparseSF

##### Sparse getSF #########
function getSF(bSF::MPIFile, frequencies, sparseTrafo::AbstractString;
               redFactor=0.1,  kargs...)

	return transformAndGetSparseSF(bSF, frequencies, sparseTrafo;
     redFactor=redFactor, kargs...)
end

# transforms selected frequencies, do compression end return sparse matrix
# parameter globalComp: if true, each frequency gets the same number of non-zero entries
#						if false, the number of non-zero entries is determined by a changing threshold
function transformAndGetSparseSF(bSF::MPIFile,frequencies,sparseTrafo::String;
    thresh=0.0, redFactor=0.1, bgcorrection=false, bgCorrection=bgcorrection,
    loadasreal=false, compAna=nothing, combine=true, useCOM=false, depth=4, useDFFoV=false, kargs...)

    if useDFFoV
      #z = findCenterOfDfFov(bSF;dispRes=false,combine=combine,useCOM=useCOM,depth=depth)

      dfFovSize=dfFov(bSF)
      fovSize=calibFov(bSF)
      gridSF=calibSize(bSF)
      # adding one and floor afterwards to avoid zeros in the grid
      dfGridSize = floor.(Int, gridSF .* dfFovSize ./ fovSize .+ 1)

      origin = RegularGridPositions(calibSize(bSF),calibFov(bSF),[0.0,0.0,0.0])
      grid = RegularGridPositions(dfGridSize, dfFovSize, [0.0,0.0,0.0])
    else
      grid = RegularGridPositions(calibSize(bSF),calibFov(bSF),[0.0,0.0,0.0])
    end

    basisTrafo = linearOperator(sparseTrafo,shape(grid))

    nFreq = rxNumFrequencies(bSF)*rxNumChannels(bSF)
    N = length(grid)
    NRed = max(1,floor(Int,redFactor*N))
    l = length(frequencies)

    buffer = zeros(ComplexF32, N)
    indices = [zeros(Int32,1) for i=1:l]
    data = [zeros(ComplexF32,1) for i=1:l]
    numCoeff = zeros(Int32,l)

    for k = 1:l
        SF = vec(map(ComplexF32, systemMatrix(bSF, frequencies[k], bgCorrection) ))

        if useDFFoV
          A = MPIFiles.interpolate(reshape(SF,shape(origin)...), origin, grid)
          buffer = vec(A)
        else
          buffer = SF
        end

        buffer[:] = transpose(basisTrafo) * buffer
        # compression
        if thresh == 0.0
            indices[k] = round.(Int32,reverse(sortperm(abs.(buffer)),dims=1)[1:NRed])
        else
            max, = findmax(abs.(buffer))
            t = thresh*max
            indices[k] = findall(x->x>t,abs.(buffer))
        end
        numCoeff[k] = length(indices[k])
        data[k] = buffer[indices[k]]
        #calcSigmaStep(compAna,buffer,indices[k],k)
    end

    @info "compression factor = " sum(numCoeff)/(N*l)

    sparse = loadsparsedata(bSF,data,indices,l,N,numCoeff,loadasreal)

    return sparse, grid
end

function loadsparsedata(f,data,indices,l,nPos,numCoeff,loadasreal::Bool)
  N = nPos
  M = l
  S = ComplexF32[]
  I = Int64[]
  for k=1:l
    append!(S, data[k])
    append!(I, indices[k])
  end
  if loadasreal
    S = reshape(S,sum(numCoeff),1)
    S = converttoreal(S)
    M *= 2
    I = vcat(I,I)
    indptr = collect(0:l*2) .* [0;numCoeff;numCoeff] .+ 1
  else
    #indptr = collect(0:l) .* [0;numCoeff] .+ 1
    indptr = cumsum([0;numCoeff]) .+ 1
  end
  return SparseMatrixCSC(N,M,indptr,I[:],S[:])
end


text:
# In this example, we create piece-wise linear functions that has a single line segment that has a small slope. The interpretation of using such a function as a transport map is to move mass from a large coverage to a concentrated coverage. The intervals in the domain and range that correspond to this segment is the *focus interval*, and is controlled by `range_percentage` and `domain_percentage`. We then fit a logistic-logit function to each of the piece-wise linear functions. The fitted parameters can then be interpolated so that we get a logistic-logit function for any given focus interval.

include("../src/MonotoneMaps.jl")
import .MonotoneMaps
#import MonotoneMaps

import PyPlot # works better in REPL.

#using Plots; pyplot() # easier to work with in .jmd / Weave.jl.

import Random
using LinearAlgebra

PyPlot.close("all")
fig_num = 1

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

Random.seed!(25)

# Get the set of piece-wise linear functions.
p_lb = 0.0
p_ub = 1.0
range_percentage = 0.1
N_itp_samples = 10
domain_percentage = 0.7
infos, zs, p_range = MonotoneMaps.getendomorphismpiecewiselinear(p_lb, p_ub, range_percentage; N_itp_samples = N_itp_samples, domain_percentage = domain_percentage)

# construct the set of piece-wise functions.
fs = collect( xx->MonotoneMaps.evalpiecewise2Dlinearfunc(xx, infos[i]) for i = 1:length(infos) )

# visualize.
display_t = LinRange(0.0, 1.0, 5000)

PyPlot.figure(fig_num)
fig_num += 1

for i = length(fs):length(fs)
    PyPlot.plot(display_t, fs[i].(display_t), label = "f[$(i)]")
end

PyPlot.legend()
PyPlot.xlabel("x")
PyPlot.ylabel("y")
PyPlot.title("Piece-wise linear functions, each with a single focus subinterval")

# display_mat = [ fs[c](display_t[r]) for r = 1:length(display_t), c = 1:length(fs) ]
# display_labels = ["$(i)-th function" for j = 1:1, i = 1:length(fs)]
#
# default(titlefont = (20, "times"), legendfontsize = 15, guidefont = (18, :black),
# tickfont = (12, :black))
#
# plot_handle = plot(display_t, display_mat,
# label = display_labels,
# title = "Piece-wise linear functions, each with a single focus subinterval",
# xlabel = "x",
# ylabel = "y",
# linewidth = 2, legend = :outerright, aspect_ratio=:equal, size = (800,800))
# display(plot_handle)

# fit the compact sigmoids (composite function of applying probit then logistic functions).
p0 = [0.5; 0.0]
p_lb = [0.1; -5.0]
p_ub = [0.6; 5.0]
optim_algorithm = :LN_BOBYQA # a local derivative-free algorithm. For other algorithms in NLopt, see https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/
costfuncs, minxs, rets = MonotoneMaps.getcompactsigmoidparameters(infos; p0 = p0, p_lb = p_lb, p_ub = p_ub, optim_algorithm = optim_algorithm)
qs = collect( tt->MonotoneMaps.evalcompositelogisticprobit(tt, minxs[i][1], minxs[i][2]) for i = 1:length(minxs) )

## visualize oracle vs. fitted sigmoids.
PyPlot.figure(fig_num)
fig_num += 1

for l = 1:length(qs)
    PyPlot.plot(display_t, fs[l].(display_t), label = "f[$(l)]")
    PyPlot.plot(display_t, qs[l].(display_t), "--", label = "q[$(l)]")
end

PyPlot.legend()
PyPlot.xlabel("x")
PyPlot.ylabel("")
PyPlot.title("target vs fit")

# display_mat = [ fs[c](display_t[r]) for r = 1:length(display_t), c = 1:length(fs) ]
# display_labels = ["$(i)-th piece-wise linear" for j = 1:1, i = 1:length(fs)]
#
# display_mat2 = [ qs[c](display_t[r]) for r = 1:length(display_t), c = 1:length(fs) ]
# display_labels2 = ["$(i)-th logistic-logit" for j = 1:1, i = 1:length(qs)]
#
# plot_handle2 = plot(display_t, display_mat,
# label = display_labels,
# title = "Piece-wise linear and logistic-logit functions",
# xlabel = "x",
# ylabel = "y",
# linewidth = 2, legend = :outerright, aspect_ratio=:equal, size = (800,800))
#
# plot!(display_t, display_mat2,
# label = display_labels2,
# style = :dash)
#
# display(plot_handle2)

# print fit cost.
for l = 1:length(costfuncs)
    println("cost of p_star $(l) is ", costfuncs[l](minxs[l]))
end


text:
# dirichlet distribution sample

module DirichletDistribution
  using Distributions

  function main()
    tau = [1, 2, 3, 4, 5]

    for i in 1:3
      pi_ = rand(Dirichlet(tau))
      println("π = $(pi_), Sum: $(sum(pi_))")
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
"""
    poly(vertices, indices; kwargs...)
    poly(points; kwargs...)
    poly(shape; kwargs...)
    poly(mesh; kwargs...)

Plots a polygon based on the arguments given.
When vertices and indices are given, it functions similarly to `mesh`.
When points are given, it draws one polygon that connects all the points in order.
When a shape is given (essentially anything decomposable by `GeometryBasics`), it will plot `decompose(shape)`.

    poly(coordinates, connectivity; kwargs...)

Plots polygons, which are defined by
`coordinates` (the coordinates of the vertices) and
`connectivity` (the edges between the vertices).

## Attributes
$(ATTRIBUTES)
"""
@recipe(Poly) do scene
    Attributes(;
        color = theme(scene, :color),
        visible = theme(scene, :visible),
        strokecolor = RGBAf0(0,0,0,0),
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        strokewidth = 1.0,
        shading = false,
        # we turn this false for now, since otherwise shapes look transparent
        # since we use meshes, which are drawn into a different framebuffer because of fxaa
        # if we use fxaa=false, they're drawn into the same
        # TODO, I still think this is a bug, since they should still use the same depth buffer!
        fxaa = false,
        linestyle = nothing,
        overdraw = false,
        transparency = false,
    )
end

const PolyElements = Union{Polygon, MultiPolygon, Circle, Rect, AbstractMesh, VecTypes, AbstractVector{<:VecTypes}}

convert_arguments(::Type{<: Poly}, v::AbstractVector{<: PolyElements}) = (v,)
convert_arguments(::Type{<: Poly}, v::Union{Polygon, MultiPolygon}) = (v,)

convert_arguments(::Type{<: Poly}, args...) = ([convert_arguments(Scatter, args...)[1]],)
convert_arguments(::Type{<: Poly}, vertices::AbstractArray, indices::AbstractArray) = convert_arguments(Mesh, vertices, indices)
convert_arguments(::Type{<: Poly}, m::GeometryBasics.Mesh) = (m,)

function plot!(plot::Poly{<: Tuple{Union{GeometryBasics.Mesh, GeometryPrimitive}}})
    mesh!(
        plot, plot[1],
        color = plot[:color], colormap = plot[:colormap], colorrange = plot[:colorrange],
        shading = plot[:shading], visible = plot[:visible], overdraw = plot[:overdraw]
    )
    wireframe!(
        plot, plot[1],
        color = plot[:strokecolor], linestyle = plot[:linestyle],
        linewidth = plot[:strokewidth], visible = plot[:visible], overdraw = plot[:overdraw]
    )
end

# Poly conversion
poly_convert(geometries) = triangle_mesh.(geometries)
poly_convert(meshes::AbstractVector{<:AbstractMesh}) = meshes
poly_convert(polys::AbstractVector{<:Polygon}) = triangle_mesh.(polys)
function poly_convert(multipolygons::AbstractVector{<:MultiPolygon})
    return [merge(triangle_mesh.(multipoly.polygons)) for multipoly in multipolygons]
end

poly_convert(mesh::GeometryBasics.Mesh) = mesh

poly_convert(polygon::Polygon) = triangle_mesh(polygon)

function poly_convert(polygon::AbstractVector{<: VecTypes})
    return poly_convert([convert_arguments(Scatter, polygon)[1]])
end

function poly_convert(polygons::AbstractVector{<: AbstractVector{<: VecTypes}})
    return map(polygons) do poly
        s = GeometryBasics.split_intersections(poly)
        merge(triangle_mesh.(Polygon.(s)))
    end
end

to_line_segments(polygon) = convert_arguments(PointBased(), polygon)[1]
to_line_segments(polygon::GeometryBasics.Mesh) = convert_arguments(PointBased(), polygon)[1]

function to_line_segments(meshes::AbstractVector)
    line = Point2f0[]
    for (i, mesh) in enumerate(meshes)
        points = to_line_segments(mesh)
        append!(line, points)
        # push!(line, points[1])
        # dont need to separate the last line segment
        if i != length(meshes)
            push!(line, Point2f0(NaN))
        end
    end
    return line
end

function to_line_segments(polygon::AbstractVector{<: VecTypes})
    result = Point2f0.(polygon)
    push!(result, polygon[1])
    return result
end

function plot!(plot::Poly{<: Tuple{<: Union{Polygon, AbstractVector{<: PolyElements}}}})
    geometries = plot[1]
    meshes = lift(poly_convert, geometries)
    mesh!(plot, meshes;
        visible = plot.visible,
        shading = plot.shading,
        color = plot.color,
        colormap = plot.colormap,
        colorrange = plot.colorrange,
        overdraw = plot.overdraw,
        fxaa = plot.fxaa,
        transparency = plot.transparency
    )
    outline = lift(to_line_segments, geometries)
    lines!(
        plot, outline, visible = plot.visible,
        color = plot.strokecolor, linestyle = plot.linestyle,
        linewidth = plot.strokewidth,
        overdraw = plot.overdraw, transparency = plot.transparency
    )
end

function plot!(plot::Mesh{<: Tuple{<: AbstractVector{P}}}) where P <: Union{AbstractMesh, Polygon}
    meshes = plot[1]
    color_node = plot.color
    attributes = Attributes(visible = plot.visible, shading = plot.shading, fxaa=plot.fxaa)

    attributes[:colormap] = get(plot, :colormap, nothing)
    attributes[:colorrange] = get(plot, :colorrange, nothing)

    bigmesh = if color_node[] isa AbstractVector && length(color_node[]) == length(meshes[])
        # One color per mesh
        lift(meshes, color_node, attributes.colormap, attributes.colorrange) do meshes, colors, cmap, crange
            # Color are reals, so we need to transform it to colors first
            single_colors = if colors isa AbstractVector{<:Number}
                interpolated_getindex.((to_colormap(cmap),), colors, (crange,))
            else
                to_color.(colors)
            end
            real_colors = RGBAf0[]
            # Map one single color per mesh to each vertex
            for (mesh, color) in zip(meshes, single_colors)
                append!(real_colors, Iterators.repeated(RGBAf0(color), length(coordinates(mesh))))
            end
            # real_colors[] = real_colors[]
            if P <: AbstractPolygon
                meshes = triangle_mesh.(meshes)
            end
            return pointmeta(merge(meshes), color=real_colors)
        end
    else
        attributes[:color] = color_node
        lift(meshes) do meshes
            return merge(GeometryBasics.mesh.(meshes))
        end
    end
    mesh!(plot, attributes, bigmesh)
end


text:
macro def(name, definition)
  esc(quote
    macro $name()
      esc($(Expr(:quote, definition)))
    end
  end)
end

# TODO decide if this is wanted and/or necessary
@def add_generic_fields begin
    method_string::String
    n::Int64
    x::Array{T}
    f_x::T
    f_calls::Int64
    g_calls::Int64
    h_calls::Int64
end

@def add_linesearch_fields begin
    x_ls::Array{T}
    g_ls::Array{T}
    alpha::T
    mayterminate::Bool
    lsr::LineSearches.LineSearchResults
end

@def initial_linesearch begin
    (similar(initial_x), # Buffer of x for line search in state.x_ls
    similar(initial_x), # Buffer of g for line search in state.g_ls
    LineSearches.alphainit(one(T), initial_x, g, f_x), # Keep track of step size in state.alpha
    false, # state.mayterminate
    LineSearches.LineSearchResults(T))
end


text:
part1(input) = simulate(input, 80)
part2(input) = simulate(input, 256)

function simulate(input, steps)
    x0 = x1 = x2 = x3 = x4 = x5 = x6 = x7 = x8 = 0
    for c in read(input)[1:2:end]
        if c == UInt8('1')
            x1 += 1
        elseif c == UInt8('2')
            x2 += 1
        elseif c == UInt8('3')
            x3 += 1
        elseif c == UInt8('4')
            x4 += 1
        elseif c == UInt8('5')
            x5 += 1
        end
    end
    while true
        x7 += x0; (steps -= 1) > 0 || break
        x8 += x1; (steps -= 1) > 0 || break
        x0 += x2; (steps -= 1) > 0 || break
        x1 += x3; (steps -= 1) > 0 || break
        x2 += x4; (steps -= 1) > 0 || break
        x3 += x5; (steps -= 1) > 0 || break
        x4 += x6; (steps -= 1) > 0 || break
        x5 += x7; (steps -= 1) > 0 || break
        x6 += x8; (steps -= 1) > 0 || break
    end
    return x0 + x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8
end


text:
# outputfile for 4 components of the system - useful for electrolyte solutions with more than one salt
function writefile(il::String,c::String,time, prob_an, prob_an2, prob_cat, prob_wat)
  file = open("timecorr-$il-$c.dat","w")
  @printf(file,"# Data for survivor probability calculation - %s \n", il)
  @printf(file,"# 1 - Simulation time (ns)\n")
  @printf(file,"# 2 - ANION1  timecorrelation function\n")
  @printf(file,"# 3 - ANION2  timecorrelation function\n")
  @printf(file,"# 4 - CATION  timecorrelation function\n")
  @printf(file,"# 5 - WATER   timecorrelation function\n")
  for i in 1:length(time)
    @printf(file,"%f %f %f %f %f\n",time[i],prob_an[i],prob_an2[i], prob_cat[i], prob_wat[i])
  end
end

# outputfile for 3 components of the system - useful for electrolyte solutions
function writefile(il::String,c::String,time, prob_an, prob_cat, prob_wat)
  file = open("timecorr-$il-$c.dat","w")
  @printf(file,"# Data for survivor probability calculationi - %s \n", il)
  @printf(file,"# 1 - Simulation time (ns)\n")
  @printf(file,"# 2 - ANION  timecorrelation function\n")
  @printf(file,"# 3 - CATION timecorrelation function\n")
  @printf(file,"# 4 - WATER  timecorrelation function\n")
  for i in 1:length(time)
    @printf(file,"%f %f %f %f\n",time[i],prob_an[i], prob_cat[i], prob_wat[i])
  end
end

# output file for two components of the system (solvent)
function writefile(il::String,c::String,time, prob_an, prob_cat, name1::String, name2::String)
  file = open("timecorr-$il-$c-ions.dat","w")
  @printf(file,"# Data for survivor probability calculation\n")
  @printf(file,"# 1 - Simulation time (ns)\n")
  @printf(file,"# 2 - %s probability survivor function\n", name1)
  @printf(file,"# 3 - %s probability survivor function\n", name2)
  for i in 1:length(time)
    @printf(file,"%f %f %f\n",time[i],prob_an[i], prob_cat[i])
  end
end

# output file for one component of the system
function writefile(il::String,c::String,time, prob_wat, name1::String)
  file = open("timecorr-$il-$c-water.dat","w")
  @printf(file,"# Data for survivor probability calculation\n")
  @printf(file,"# 1 - Simulation time (ns)\n")
  @printf(file,"# 2 - %s  timecorrelation function\n", name1)
  for i in 1:length(time)
    @printf(file,"%f %f\n", time[i], prob_wat[i])
  end
end


text:
using FactCheck, GeoJSON, Turf, Compat

# small tests
polygon = Vector[Vector{Float64}[[-1.,-1.],[1.,-1.],[1.,1.],[-1.,1.]]]
point = [0.9,0.9]
@fact inside(Point(point), Polygon(polygon)) --> true
point = [1.0,1.1]
@fact inside(Point(point), Polygon(polygon)) --> false

polygon = Vector[Vector{Float64}[[-1.,-1.],[1.,-1.],[1.,1.],[-1.,1.]],
                 Vector{Float64}[[-0.5,-0.5],[0.5,-0.5],[0.5,0.5],[-0.5,0.5]]]

point = [-0.9,-0.9]
@fact inside(Point(point), Polygon(polygon)) --> true
point = [-1.0,-1.1]
@fact inside(Point(point), Polygon(polygon)) --> false
point = [0.0,0.1]
@fact inside(Point(point), Polygon(polygon)) --> false

# test simple polygon
polygon = GeoJSON.Polygon(Vector{Vector{Float64}}[Vector{Float64}[[0,0],[0,100],[100,0],[0,0]]])
point_in = GeoJSON.Point(Float64[50, 50])
point_out = GeoJSON.Point(Float64[140, 150])
@fact Turf.inside(point_in, polygon) --> true
@fact Turf.inside(point_out, polygon) --> false

# test concave polygon
polygon = GeoJSON.Polygon(Vector{Vector{Float64}}[Vector{Float64}[[0,0],[50,50],[0,100],[100,100],[100,0],[0,0]]])
point_in = GeoJSON.Point(Float64[75, 75])
point_out = GeoJSON.Point(Float64[25, 50])
@fact Turf.inside(point_in, polygon) --> true
@fact Turf.inside(point_out, polygon) --> false

# test polygon with hole
feature = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/poly-with-hole.geojson"))
polygon = feature.geometry
in_hole = GeoJSON.Point(Float64[-86.69208526611328, 36.20373274711739])
point_in = GeoJSON.Point(Float64[-86.72229766845702, 36.20258997094334])
point_out = GeoJSON.Point(Float64[-86.75079345703125, 36.18527313913089])
@fact Turf.inside(in_hole, polygon) --> false
@fact Turf.inside(point_in, polygon) --> true
@fact Turf.inside(point_out, polygon) --> false

# test multipolygon with hole
feature = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/multipoly-with-hole.geojson"))
polygon = feature.geometry
in_hole = GeoJSON.Point(Float64[-86.69208526611328, 36.20373274711739])
point_in1 = GeoJSON.Point(Float64[-86.72229766845702, 36.20258997094334])
point_in2 = GeoJSON.Point(Float64[-86.75079345703125, 36.18527313913089])
point_out = GeoJSON.Point(Float64[-86.75302505493164, 36.23015046460186])
@fact Turf.inside(in_hole, polygon) --> false
@fact Turf.inside(point_in1, polygon) --> true
@fact Turf.inside(point_in2, polygon) --> true
@fact Turf.inside(point_out, polygon) --> false

# test with a single point
polygon = GeoJSON.Polygon(Vector{Vector{Float64}}[Vector{Float64}[[0,0],[0,100],[100,0],[0,0]]])
point_in = GeoJSON.Point(Float64[50, 50])
polygonFC = GeoJSON.FeatureCollection(GeoJSON.Feature[GeoJSON.Feature(polygon)])
pointFC = GeoJSON.FeatureCollection(GeoJSON.Feature[GeoJSON.Feature(point_in)])
counted = Turf.within(pointFC, polygonFC)
@fact typeof(counted) --> GeoJSON.FeatureCollection
@fact length(counted.features) --> 1

# test with multiple points and multiple polygons
polygon1 = GeoJSON.Polygon(Vector{Vector{Float64}}[Vector{Float64}[[0,0],[10,0],[10,10],[0,10]]])
polygon2 = GeoJSON.Polygon(Vector{Vector{Float64}}[Vector{Float64}[[10,0],[20,10],[20,20],[20,0]]])
polygonFC = GeoJSON.FeatureCollection(GeoJSON.Feature[GeoJSON.Feature(polygon1),GeoJSON.Feature(polygon2)])
point1 = GeoJSON.Feature(GeoJSON.Point(Float64[1, 1]), @compat Dict{AbstractString,Any}("population" => 500))
point2 = GeoJSON.Feature(GeoJSON.Point(Float64[1, 3]), @compat Dict{AbstractString,Any}("population" => 400))
point3 = GeoJSON.Feature(GeoJSON.Point(Float64[14, 2]), @compat Dict{AbstractString,Any}("population" => 600))
point4 = GeoJSON.Feature(GeoJSON.Point(Float64[13, 1]), @compat Dict{AbstractString,Any}("population" => 500))
point5 = GeoJSON.Feature(GeoJSON.Point(Float64[19, 7]), @compat Dict{AbstractString,Any}("population" => 200))
point6 = GeoJSON.Feature(GeoJSON.Point(Float64[100, 7]), @compat Dict{AbstractString,Any}("population" => 200))
pointFC = GeoJSON.FeatureCollection(GeoJSON.Feature[point1, point2, point3, point4, point5, point6])
counted = Turf.within(pointFC, polygonFC)
@fact typeof(counted) --> GeoJSON.FeatureCollection
@fact length(counted.features) --> 5

# test tag (which is dependent on filter)
pointFC = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/tagPoints.geojson"))
polygonFC = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/tagPolygons.geojson"))
taggedFC = Turf.tag(pointFC, polygonFC, "polyID", "containingPolyID")
@fact length(taggedFC.features) --> length(pointFC.features)
countFC = Turf.filter(taggedFC, "containingPolyID", 4)
@fact length(countFC.features) --> 6


text:
# Container data structure (structure-of-arrays style) for DG elements
# TODO: Taal refactor, remove NVARS, POLYDEG?
mutable struct ElementContainer1D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer
  inverse_jacobian::Vector{RealT}        # [elements]
  node_coordinates::Array{RealT, 3}      # [orientation, i, elements]
  surface_flux_values::Array{uEltype, 3} # [variables, direction, elements]
  cell_ids::Vector{Int}                  # [elements]
  # internal `resize!`able storage
  _node_coordinates::Vector{RealT}
  _surface_flux_values::Vector{uEltype}
end

nvariables(::ElementContainer1D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = NVARS
polydeg(::ElementContainer1D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = POLYDEG
Base.eltype(::ElementContainer1D{RealT, uEltype}) where {RealT, uEltype} = uEltype

# Only one-dimensional `Array`s are `resize!`able in Julia.
# Hence, we use `Vector`s as internal storage and `resize!`
# them whenever needed. Then, we reuse the same memory by
# `unsafe_wrap`ping multi-dimensional `Array`s around the
# internal storage.
function Base.resize!(elements::ElementContainer1D, capacity)
  n_nodes = polydeg(elements) + 1
  nvars = nvariables(elements)
  @unpack _node_coordinates, _surface_flux_values,
          inverse_jacobian, cell_ids = elements

  resize!(inverse_jacobian, capacity)

  resize!(_node_coordinates, 1 * n_nodes * capacity)
  elements.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
                                          (1, n_nodes, capacity))

  resize!(_surface_flux_values, nvars * 2 * 1 * capacity)
  elements.surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values),
                                             (nvars, 2 * 1, capacity))

  resize!(cell_ids, capacity)

  return nothing
end

function ElementContainer1D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  nan_RealT = convert(RealT, NaN)
  nan_uEltype = convert(uEltype, NaN)

  # Initialize fields with defaults
  inverse_jacobian = fill(nan_RealT, capacity)

  _node_coordinates = fill(nan_RealT, 1 * n_nodes * capacity)
  node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
                                 (1, n_nodes, capacity))

  _surface_flux_values = fill(nan_uEltype, NVARS * 2 * 1 * capacity)
  surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values),
                                    (NVARS, 2 * 1, capacity))

  cell_ids = fill(typemin(Int), capacity)

  return ElementContainer1D{RealT, uEltype, NVARS, POLYDEG}(
    inverse_jacobian, node_coordinates, surface_flux_values, cell_ids,
    _node_coordinates, _surface_flux_values)
end

# Return number of elements
@inline nelements(elements::ElementContainer1D) = length(elements.cell_ids)
# TODO: Taal performance, 1:nelements(elements) vs. Base.OneTo(nelements(elements))
@inline eachelement(elements::ElementContainer1D) = Base.OneTo(nelements(elements))

# Create element container and initialize element data
function init_elements(cell_ids, mesh::TreeMesh1D,
                       equations::AbstractEquations{1, NVARS},
                       basis::LobattoLegendreBasis{T, NNODES}, ::Type{RealT}, ::Type{uEltype}) where {RealT<:Real, uEltype<:Real,NVARS, T, NNODES}
  # Initialize container
  n_elements = length(cell_ids)
  elements = ElementContainer1D{RealT, uEltype, NVARS, NNODES-1}(n_elements)

  init_elements!(elements, cell_ids, mesh, basis.nodes)
  return elements
end

function init_elements!(elements, cell_ids, mesh::TreeMesh1D, nodes)
  # Store cell ids
  elements.cell_ids .= cell_ids

  # Calculate inverse Jacobian and node coordinates
  for element in eachelement(elements)
    # Get cell id
    cell_id = cell_ids[element]

    # Get cell length
    dx = length_at_cell(mesh.tree, cell_id)

    # Calculate inverse Jacobian as 1/(h/2)
    elements.inverse_jacobian[element] = 2/dx

    # Calculate node coordinates
      for i in eachindex(nodes)
        elements.node_coordinates[1, i, element] = (
            mesh.tree.coordinates[1, cell_id] + dx/2 * nodes[i])
      end
  end

  return elements
end

# Container data structure (structure-of-arrays style) for DG interfaces
# TODO: Taal refactor, remove NVARS, POLYDEG?
mutable struct InterfaceContainer1D{uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer
  u::Array{uEltype, 3}      # [leftright, variables, interfaces]
  neighbor_ids::Matrix{Int} # [leftright, interfaces]
  orientations::Vector{Int} # [interfaces]
  # internal `resize!`able storage
  _u::Vector{uEltype}
  _neighbor_ids::Vector{Int}
end

nvariables(::InterfaceContainer1D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = NVARS
polydeg(::InterfaceContainer1D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = POLYDEG
Base.eltype(::InterfaceContainer1D{uEltype}) where {uEltype} = uEltype

# See explanation of Base.resize! for the element container
function Base.resize!(interfaces::InterfaceContainer1D, capacity)
  n_nodes = polydeg(interfaces) + 1
  nvars = nvariables(interfaces)
  @unpack _u, _neighbor_ids, orientations = interfaces

  resize!(_u, 2 * nvars * capacity)
  interfaces.u = unsafe_wrap(Array, pointer(_u),
                             (2, nvars, capacity))

  resize!(_neighbor_ids, 2 * capacity)
  interfaces.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ),
                                        (2, capacity))

  resize!(orientations, capacity)

  return nothing
end

function InterfaceContainer1D{uEltype, NVARS, POLYDEG}(capacity::Integer) where {uEltype<:Real, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  nan = convert(uEltype, NaN)

  # Initialize fields with defaults
  _u = fill(nan, 2 * NVARS * capacity)
  u = unsafe_wrap(Array, pointer(_u),
                  (2, NVARS, capacity))

  _neighbor_ids = fill(typemin(Int), 2 * capacity)
  neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ),
                             (2, capacity))

  orientations = fill(typemin(Int), capacity)

  return InterfaceContainer1D{uEltype, NVARS, POLYDEG}(
    u, neighbor_ids, orientations,
    _u, _neighbor_ids)
end

# Return number of interfaces
@inline ninterfaces(interfaces::InterfaceContainer1D) = length(interfaces.orientations)

# Create interface container and initialize interface data in `elements`.
function init_interfaces(cell_ids, mesh::TreeMesh1D,
                         elements::ElementContainer1D{RealT, uEltype, NVARS, POLYDEG}) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG}
  # Initialize container
  n_interfaces = count_required_interfaces(mesh, cell_ids)
  interfaces = InterfaceContainer1D{uEltype, NVARS, POLYDEG}(n_interfaces)

  # Connect elements with interfaces
  init_interfaces!(interfaces, elements, mesh)
  return interfaces
end

# Count the number of interfaces that need to be created
function count_required_interfaces(mesh::TreeMesh1D, cell_ids)
  count = 0

  # Iterate over all cells
  for cell_id in cell_ids
    for direction in eachdirection(mesh.tree)
      # Only count interfaces in positive direction to avoid double counting
      if direction == 1
        continue
      end

      # Skip if no neighbor exists
      if !has_any_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      count += 1
    end
  end

  return count
end

# Initialize connectivity between elements and interfaces
function init_interfaces!(interfaces, elements, mesh::TreeMesh1D)
  # Construct cell -> element mapping for easier algorithm implementation
  tree = mesh.tree
  c2e = zeros(Int, length(tree))
  for element in eachelement(elements)
    c2e[elements.cell_ids[element]] = element
  end

  # Reset interface count
  count = 0

  # Iterate over all elements to find neighbors and to connect via interfaces
  for element in eachelement(elements)
    # Get cell id
    cell_id = elements.cell_ids[element]

    # Loop over directions
    for direction in eachdirection(mesh.tree)
      # Only create interfaces in positive direction
      if direction == 1
        continue
      end

      # Skip if no neighbor exists and current cell is not small
      if !has_any_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      count += 1

      if has_neighbor(mesh.tree, cell_id, direction)
        neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id]
        if has_children(mesh.tree, neighbor_cell_id) # Cell has small neighbor
          interfaces.neighbor_ids[2, count] = c2e[mesh.tree.child_ids[1, neighbor_cell_id]]
        else # Cell has same refinement level neighbor
          interfaces.neighbor_ids[2, count] = c2e[neighbor_cell_id]
        end
      else # Cell is small and has large neighbor
        parent_id = mesh.tree.parent_ids[cell_id]
        neighbor_cell_id = mesh.tree.neighbor_ids[direction, parent_id]
        interfaces.neighbor_ids[2, count] = c2e[neighbor_cell_id]
      end

      interfaces.neighbor_ids[1, count] = element
      # Set orientation (x -> 1)
      interfaces.orientations[count] = 1
    end
  end

  @assert count == ninterfaces(interfaces) ("Actual interface count ($count) does not match " *
                                            "expectations $(ninterfaces(interfaces))")
end

# Container data structure (structure-of-arrays style) for DG boundaries
# TODO: Taal refactor, remove NVARS, POLYDEG?
mutable struct BoundaryContainer1D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer
  u::Array{uEltype, 3}              # [leftright, variables, boundaries]
  neighbor_ids::Vector{Int}         # [boundaries]
  orientations::Vector{Int}         # [boundaries]
  neighbor_sides::Vector{Int}       # [boundaries]
  node_coordinates::Array{RealT, 2} # [orientation, elements]
  n_boundaries_per_direction::SVector{2, Int} # [direction]
  # internal `resize!`able storage
  _u::Vector{uEltype}
  _node_coordinates::Vector{RealT}
end

nvariables(::BoundaryContainer1D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = NVARS
polydeg(::BoundaryContainer1D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = POLYDEG
Base.eltype(::BoundaryContainer1D{RealT, uEltype}) where {RealT, uEltype} = uEltype

# See explanation of Base.resize! for the element container
function Base.resize!(boundaries::BoundaryContainer1D, capacity)
  n_nodes = polydeg(boundaries) + 1
  nvars = nvariables(boundaries)
  @unpack _u, _node_coordinates,
          neighbor_ids, orientations, neighbor_sides = boundaries

  resize!(_u, 2 * nvars * capacity)
  boundaries.u = unsafe_wrap(Array, pointer(_u),
                             (2, nvars, capacity))

  resize!(_node_coordinates, 1 * capacity)
  boundaries.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates ),
                                            (1, capacity))

  resize!(neighbor_ids, capacity)

  resize!(orientations, capacity)

  resize!(neighbor_sides, capacity)

  return nothing
end

function BoundaryContainer1D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  nan_RealT = convert(RealT, NaN)
  nan_uEltype = convert(uEltype, NaN)

  # Initialize fields with defaults
  _u = fill(nan_uEltype, 2 * NVARS * capacity)
  u = unsafe_wrap(Array, pointer(_u),
                  (2, NVARS, capacity))

  neighbor_ids = fill(typemin(Int), capacity)

  orientations = fill(typemin(Int), capacity)

  neighbor_sides = fill(typemin(Int), capacity)

  _node_coordinates = fill(nan_RealT, 1 * capacity)
  node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
                                 (1, capacity))

  n_boundaries_per_direction = SVector(0, 0)

  return BoundaryContainer1D{RealT, uEltype, NVARS, POLYDEG}(
    u, neighbor_ids, orientations, neighbor_sides,
    node_coordinates, n_boundaries_per_direction,
    _u, _node_coordinates)
end

# Return number of boundaries
nboundaries(boundaries::BoundaryContainer1D) = length(boundaries.orientations)

# Create boundaries container and initialize boundary data in `elements`.
function init_boundaries(cell_ids, mesh::TreeMesh1D,
                         elements::ElementContainer1D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG}
  # Initialize container
  n_boundaries = count_required_boundaries(mesh, cell_ids)
  boundaries = BoundaryContainer1D{RealT, uEltype, NVARS, POLYDEG}(n_boundaries)

  # Connect elements with boundaries
  init_boundaries!(boundaries, elements, mesh)
  return boundaries
end

# Count the number of boundaries that need to be created
function count_required_boundaries(mesh::TreeMesh1D, cell_ids)
  count = 0

  # Iterate over all cells
  for cell_id in cell_ids
    for direction in eachdirection(mesh.tree)
      # If neighbor exists, current cell is not at a boundary
      if has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # If coarse neighbor exists, current cell is not at a boundary
      if has_coarse_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # No neighbor exists in this direction -> must be a boundary
      count += 1
    end
  end

  return count
end

# Initialize connectivity between elements and boundaries
function init_boundaries!(boundaries, elements, mesh::TreeMesh1D)
  # Reset boundaries count
  count = 0

  # Initialize boundary counts
  counts_per_direction = MVector(0, 0)

  # OBS! Iterate over directions first, then over elements, and count boundaries in each direction
  # Rationale: This way the boundaries are internally sorted by the directions -x, +x, -y etc.,
  #            obviating the need to store the boundary condition to be applied explicitly.
  # Loop over directions
  for direction in eachdirection(mesh.tree)
    # Iterate over all elements to find missing neighbors and to connect to boundaries
    for element in eachelement(elements)
      # Get cell id
      cell_id = elements.cell_ids[element]

      # If neighbor exists, current cell is not at a boundary
      if has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # If coarse neighbor exists, current cell is not at a boundary
      if has_coarse_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # Create boundary
      count += 1
      counts_per_direction[direction] += 1

      # Set neighbor element id
      boundaries.neighbor_ids[count] = element

      # Set neighbor side, which denotes the direction (1 -> negative, 2 -> positive) of the element
      if direction == 2
        boundaries.neighbor_sides[count] = 1
      else
        boundaries.neighbor_sides[count] = 2
      end

      # Set orientation (x -> 1)
      boundaries.orientations[count] = 1

      # Store node coordinates
      enc = elements.node_coordinates
      if direction == 1 # -x direction
        boundaries.node_coordinates[:, count] .= enc[:, 1,  element]
      elseif direction == 2 # +x direction
        boundaries.node_coordinates[:, count] .= enc[:, end, element]
      else
        error("should not happen")
      end
    end
  end

  @assert count == nboundaries(boundaries) ("Actual boundaries count ($count) does not match " *
                                            "expectations $(nboundaries(boundaries))")
  @assert sum(counts_per_direction) == count

  boundaries.n_boundaries_per_direction = SVector(counts_per_direction)

  return boundaries.n_boundaries_per_direction
end


text:
# Add packages

using PyPlot
using DataFrames
using Vann
using DataAssim

# Model choices

snow_choice = TinBasic
hydro_choice = Gr4j

# Filter choices

filter_choice = particle_filter

nens = 50

# Load data

path_inputs = Pkg.dir("Vann", "data/atnasjo")

date, tair, prec, q_obs, frac = load_data(path_inputs)

# Compute potential evapotranspiration

epot = epot_zero(date)

# Initilize model

tstep = 24.0

time = date[1]

st_snow = eval(Expr(:call, snow_choice, tstep, time, frac))
st_hydro = eval(Expr(:call, hydro_choice, tstep, time))

# Run calibration

param_snow, param_hydro = run_model_calib(st_snow, st_hydro, tair, prec, epot, q_obs)

# Run model and filter

st_snow = eval(Expr(:call, snow_choice, tstep, time, param_snow, frac))
st_hydro = eval(Expr(:call, hydro_choice, tstep, time, param_hydro))

q_res = eval(Expr(:call, filter_choice, st_snow, st_hydro, prec, tair, epot, q_obs, nens))

# Add results to dataframe

q_sim  = q_res[:, 1]
q_min  = q_res[:, 2]
q_max  = q_res[:, 3]

q_obs = round(q_obs, 2)
q_sim = round(q_sim, 2)
q_min = round(q_min, 2)
q_max = round(q_max, 2)

df_res = DataFrame(date = date, q_obs = q_obs, q_sim = q_sim, q_min = q_min, q_max = q_max)

# Plot results

days_warmup = 3*365

df_res = df_res[days_warmup:end, :]

fig = plt[:figure](figsize = (12,7))

plt[:style][:use]("ggplot")

plt[:plot](df_res[:date], df_res[:q_obs], linewidth = 1.2, color = "k", label = "Observed", zorder = 1)
plt[:fill_between](df_res[:date], df_res[:q_max], df_res[:q_min], facecolor = "r", edgecolor = "r", label = "Simulated", alpha = 0.55, zorder = 2)
plt[:legend]()

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
## circulant_embedding.jl : Gaussian random field generator using fft; only for uniformly spaced GRFs

"""
CirculantEmbedding <: GaussianRandomFieldGenerator

A [`GaussiandRandomFieldGenerator`](@ref) that uses FFT to compute samples of the Gaussian random field. Circulant embedding can only be applied if the points are specified on a structured grid.

# Examples
```jldoctest
julia> m = Matern(0.1,1.0)
Matérn (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> c = CovarianceFunction(2,m)
2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> pts1 = 0:0.02:1; pts2 = 0:0.02:1
0.0:0.02:1.0

julia> grf = GaussianRandomField(c,CirculantEmbedding(),pts1,pts2)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a circulant embedding

julia> contourf(grf)
[...]

julia> plot_eigenvalues(grf)
[...]

```
With appropriate ordering, the covariance matrix of a Gaussian random field is a (nested block) Toeplitz matrix. This matrix can be embedded into a larger (nested block) circulant matrix, whose eigenvalues can be rapidly computed using FFT. A difficulty here is that although the covariance matrix is positive semi-definite, its circulant extension in general is not. As a remedy, one can add so-called *ghost points* outside the domain of interest using the optional flag `padding`.
```jldoctest
julia> m = Matern(1,1)
Matérn (λ=1.0, ν=1.0, σ=1.0, p=2.0)

julia> c = CovarianceFunction(1,m)
1d Matérn covariance function (λ=1.0, ν=1.0, σ=1.0, p=2.0)

julia> pts = range(0,stop = 1,length = 256)
0.0:0.00392156862745098:1.0

julia> grf = GaussianRandomField(c,CirculantEmbedding(),pts)
WARNING: negative eigenvalue -0.001465931113950698 detected, Gaussian random field will be approximated (ignoring all negative eigenvalues)
WARNING: increase padding if possible
Gaussian random field with 1d Matérn covariance function (λ=1.0, ν=1.0, σ=1.0, p=2.0) on a 256-point structured grid, using a circulant embedding

julia> grf = GaussianRandomField(c,CirculantEmbedding(),pts,padding=5)
Gaussian random field with 1d Matérn covariance function (λ=1.0, ν=1.0, σ=1.0, p=2.0) on a 256-point structured grid, using a circulant embedding

```
For faster sampling (but slower initialization), use the optional argument `measure` (default=true).

See also: [`Cholesky`](@ref), [`Spectral`](@ref), [`KarhunenLoeve`](@ref)
"""
struct CirculantEmbedding <: GaussianRandomFieldGenerator end

# whether a covariance function is even in every dimension
# this simplifies the circulant embedding algorithm
Base.iseven(::CovarianceStructure) = false
Base.iseven(::IsotropicCovarianceStructure) = true

function _GaussianRandomField(mean, cov::CovarianceFunction{d}, method::CirculantEmbedding,
                              pts::Vararg{AbstractRange,d};
                              minpadding::Union{Int,Dims{d}} = 0,
                              measure::Bool = true,
                              primes::Bool = false) where {d}
    # normalize points
    normedpts = map(x -> x .- first(x), pts)

    # compute size of minimum circulant embedding
    ρ = cov.cov
    factors = primes ? [2, 3, 5, 7] : [2]
    dims = circulant_minsize.(Ref(ρ), normedpts, minpadding, Ref{Vector{Int}}(factors))

    # compute eigenvalues of the circulant embedding
    Λ = circulant_eigvals(ρ, normedpts, dims)

    # store sqrt of eigenvalues for more efficient sampling
    n₋, λ₋ = 0, 0.0
    M = length(Λ)
    @inbounds for i in eachindex(Λ)
        λ = Λ[i]
        if λ < 0
            # reset negative eigenvalues to zero
            Λ[i] = 0

            # update statistics about negative eigenvalues
            n₋ += 1
            λ < λ₋ && (λ₋ = λ)
        else
            Λ[i] = sqrt(λ / M)
        end
    end
    λ₋ < 0 && @warn begin
        "$n₋ negative eigenvalues ≥ $λ₋ detected, Gaussian random field will be " *
        "approximated (ignoring all negative eigenvalues); increase padding if possible"
    end

    # optimize
    P = measure ? plan_fft(Λ, flags=FFTW.MEASURE) : plan_fft(Λ)
    data = (Λ, P)

    GaussianRandomField{CirculantEmbedding,typeof(cov),typeof(pts),typeof(mean),typeof(data)}(mean, cov, pts, data)
end

"""
    circulant_minsize(cov::CovarianceStructure, pts::AbstractRange,
                      minpadding::Int, factors::Vector{Int})

Return size of minimum circulant embedding of covariance matrix of points `pts` with
covariance function `cov` and minimum padding `minpadding` that can be written as product of
integers in `factors`.

Typically `factors` is chosen to be `[2]` or `[2, 3, 5, 7]` to speed up FFT computations.
"""
function circulant_minsize(cov::CovarianceStructure, pts::AbstractRange,
                           minpadding::Int, factors::Vector{Int})
    if iseven(cov)
        2 * nextprod(factors, length(pts) + minpadding - 1)
    else
        2 * nextprod(factors, length(pts) + minpadding)
    end
end

"""
    circulant_eigvals(cov::CovarianceStructure, pts::NTuple{N,AbstractRange},
                      dims::Dims{N}) where {N}

Compute eigenvalues of circulant embedding with dimensions `dims` of covariance matrix of
points `pts` with covariance function `cov`.
"""
@generated function circulant_eigvals(cov::CovarianceStructure{T},
                                      pts::NTuple{N,AbstractRange},
                                      dims::Dims{N}) where {T,N}
    quote
        if iseven(cov)
            # compute first row of circulant embedding
            # since the covariance function is even in every dimension,
            # instead of 2n only n+1 values have to be computed in every dimension
            C = Array{$T}(undef, div.(dims, 2) .+ 1)
            x = Array{$T}(undef, $N)

            @nloops $N i C d -> begin
                @inbounds x[d] = extrapolate(pts[d], i_d)
            end begin
                @inbounds (@nref $N C i) = apply(cov, x)
            end

            # compute eigenvalues of circulant embedding
            # since the covariance function is even in every dimension,
            # we can perform a real even FFT and extend the resulting array of eigenvalues
            FFTW.r2r!(C, FFTW.REDFT00)
            Λ = Array{eltype(C)}(undef, dims)
            Imax = CartesianIndex(size(Λ) .+ 2)
            @inbounds for i in CartesianIndices(Λ)
                Λ[i] = C[min(i, Imax - i)]
            end
            Λ
        else
            # compute first row of circulant embedding
            C = zeros($T, dims)
            x = Array{$T}(undef, $N)

            dims2 = dims .+ 2
            mids = div.(dims2, 2)
            @nloops $N i C d -> begin
                @inbounds begin
                    m = mids[d]
                    i_d == m && continue

                    if i_d < m
                        x[d] = extrapolate(pts[d], i_d)
                    else
                        x[d] = -extrapolate(pts[d], dims2[d] - i_d)
                    end
                end
            end begin
                @inbounds (@nref $N C i) = apply(cov, x)
            end

            # compute eigenvalues of circulant embedding
            real.(fft(C))
        end
    end
end

# extrapolate ranges
extrapolate(r::Union{StepRangeLen,LinRange}, i::Integer) = Base.unsafe_getindex(r, i)
extrapolate(r::UnitRange{T}, i::Integer) where T = convert(T, r.start + (i - 1))
extrapolate(r::Base.OneTo{T}, i::Integer) where T = convert(T, i)
extrapolate(r::AbstractRange{T}, i::Integer) where T =
    convert(T, first(r) + (i - 1)*step(r))

# returns the required dimension of the random points
randdim(grf::GaussianRandomField{CirculantEmbedding}) = length(grf.data[1])

# sample function
function _sample(grf::GaussianRandomField{CirculantEmbedding}, xi)
    v, P = grf.data

    # compute multiplication with square root of circulant embedding via FFT
    y = v .* reshape(xi, size(v))
    w = P * y

    # extract realization of random field
    μ, σ = grf.mean, std(grf.cov)
    z = Array{eltype(grf.cov)}(undef, length.(grf.pts))
    @inbounds for i in CartesianIndices(z)
        wi = w[i]
        z[i] = μ[i] + σ * (real(wi) + imag(wi))
    end
    z
end

Base.show(io::IO, ::CirculantEmbedding) = print(io, "circulant embedding")


text:
@doc raw"""
First Poincaré invariant for a guiding center particle in an θ-pinch magnetic
field of the form ``B(x,y,z) = B_0 \, e_z``.

The loop for the first Poincaré invariant is initialized by
```math
q (\tau) = \begin{pmatrix}
r_x \cos (2\pi \tau) \\
y_0 \\
r_z \sin (2\pi \tau) \\
\end{pmatrix}
```
with parameters
```math
B_0 = 1, \quad
r_x = 0.5, \quad
r_z = 0.3, \quad
y_0 = 0.0, \quad
u_0 = 0.5, \quad
\mu = 2.5 \times 10^{-6} .
```
"""
module ThetaPinchField

    import ElectromagneticFields.ThetaPinch

    ThetaPinch.@code() # inject magnetic field code

    μ_loop() = 2.5E-6

    function f_loop(t)
        μ  = 2.5E-6
        Y0 = 0.0
        u0 = 4E-4
        r0 = 0.5
        r1 = 0.3

        Xt = r0*cos(2π*t)
        Zt = r1*sin(2π*t)

        qt = [Xt, Y0, Zt, u0]

        return qt
    end

    include("guiding_center_4d_common.jl")
    include("guiding_center_4d_equations.jl")
    include("guiding_center_4d_loop.jl")
    include("guiding_center_4d_diagnostics.jl")

end


text:
module Infiltrator

using REPL.LineEdit: getproperty
using REPL
using REPL.LineEdit

export @infiltrate, @exfiltrate, @withstore, safehouse, exfiltrated

REPL_HOOKED = Ref{Bool}(false)
function __init__()
  clear_store!(store)
  if VERSION >= v"1.5.0-DEV.282"
    if isdefined(Base, :active_repl_backend)
        pushfirst!(Base.active_repl_backend.ast_transforms, exit_transform)
        REPL_HOOKED[] = true
    else
      atreplinit() do repl
        @async begin
          iter = 0
          # wait for active_repl_backend to exist
          while !isdefined(Base, :active_repl_backend) && iter < 20
            sleep(0.05)
            iter += 1
          end
          if isdefined(Base, :active_repl_backend)
            pushfirst!(Base.active_repl_backend.ast_transforms, exit_transform)
            REPL_HOOKED[] = true
          end
        end
      end
    end
  end
end

"""
    @infiltrate cond = true

`@infiltrate` sets an infiltration point (or breakpoint).

When the infiltration point is hit, it will drop you into an interactive REPL session that
lets you inspect local variables and the call stack as well as execute aribtrary statements
in the context of the current functions module.

This macro also accepts an optional argument `cond` that must evaluate to a boolean,
and then this macro will serve as a "conditinal breakpoint", which starts inspections only
when its condition is `true`.
```
"""
macro infiltrate(cond = true)
  quote
    if $(esc(cond))
      $(start_prompt)(@__MODULE__, Base.@locals, $(String(__source__.file)), $(__source__.line))
    end
  end
end

"""
    @exfiltrate

Assigns all local variables into the global storage.
"""
macro exfiltrate()
  quote
    for (k, v) in Base.@locals
      try
        Core.eval(getfield(getfield($(@__MODULE__), :store), :store), Expr(:(=), k, QuoteNode(v)))
      catch err
        println(stderr, "Assignment to store variable failed.")
        Base.display_error(stderr, err, catch_backtrace())
      end
    end
  end
end

const prompt_color = if Sys.iswindows()
  "\e[33m"
else
  "\e[38;5;166m"
end

# these are used for easier testing
const TEST_TERMINAL_REF = Ref{Any}(nothing)
const TEST_REPL_REF = Ref{Any}(nothing)
const TEST_NOSTACK = Ref{Any}(false)

mutable struct Session
  store::Module
  exiting::Bool
  disabled::Set
end
function Base.show(io::IO, s::Session)
  n = length(get_store_names(s))
  print(io, "Safehouse ($(n) variable$(n == 1 ? "" : "s"))")
end
function Base.show(io::IO, ::MIME"text/plain", s::Session)
  names = get_store_names(s)
  n = length(names)
  println(io, "Safehouse with $(n) stored variable$(n == 1 ? "" : "s"):")
  for (i, (var, val)) in enumerate(names)
    if i > 6
      println(io, "...")
      break
    end
    one_line_show(io, var, val)
  end
end
function Base.getproperty(sp::Session, s::Symbol)
  m = getfield(sp, :store)
  if isdefined(m, s)
    getproperty(m, s)
  else
    throw(UndefVarError(s))
  end
end
Base.propertynames(s::Session) = keys(get_store_names(s))

"""
    store
    safehouse
    exfiltrated

Global storage for storing values while `@infiltrate`ing or `@exfiltrate`ing.
"""
const store = Session(Module(), false, Set())
const safehouse = store
const exfiltrated = store

"""
    @withstore ex

Evaluates the expression `ex` in the context of the global store.

Mainly intended for interactive use, as changes to the stores state will not
propagate into the returned expression.
"""
macro withstore(ex)
  withstore(ex)
end

function withstore(ex)
  ns = get_store_names(store)
  return Expr(:let,
    Expr(:block, map(x->Expr(:(=), x...), [(k, maybe_quote(v)) for (k, v) in ns])...),
    Expr(:block, ex)
  )
end

"""
    clear_disabled!(s = safehouse)

Clear all disabled infiltration points.
"""
function clear_disabled!(s = store)
  empty!(getfield(s, :disabled))
  return nothing
end

"""
    end_session!(s = safehouse)

End this infiltration session (reverts the effect of `@exit` in the `debug>` REPL).

Only needs to be manually called on Julia versions prior to 1.5.
"""
function end_session!(s::Session = store)
  setfield!(s, :exiting, false)
  return nothing
end

"""
    clear_store!(s = safehouse)

Reset the store used for global symbols.
"""
clear_store!(s::Session = store) = set_store!(s, Module())

"""
    set_store!(s = safehouse, m::Module)

Set the module backing the store `s`.
"""
set_store!(m::Module) = set_store!(store, m)
function set_store!(s::Session, m::Module)
  setfield!(s, :store, m)
  nothing
end

function start_prompt(mod, locals, file, fileline;
                        terminal = TEST_TERMINAL_REF[],
                        repl = TEST_REPL_REF[],
                        nostack = TEST_NOSTACK[]
                      )
  getfield(store, :exiting) && return
  (file, fileline) in getfield(store, :disabled) && return

  if terminal === nothing || repl === nothing
    if isdefined(Base, :active_repl) && isdefined(Base.active_repl, :t)
      repl = Base.active_repl
      terminal = Base.active_repl.t
    else
      println("Infiltrator.jl needs a proper Julia REPL.")
      return
    end
  end
  io = Base.pipe_writer(terminal)

  trace = stacktrace()
  start = something(findlast(x -> x.func === Symbol("start_prompt"), trace), 0) + 2
  last = something(findlast(x -> x.func === Symbol("eval"), trace),
                   length(trace)) - 2
  trace = trace[start:last]
  current = trace[1]

  println(io, "Infiltrating $current:")
  println(io)
  debugprompt(mod, locals, trace, terminal, repl, nostack, file=file, fileline=fileline)
  println(io)
end

function show_help(io)
  println(io, """
    Code entered is evaluated in the current functions module. Note that you cannot change local
    variables, but can assign to globals in a permanent store module.

    The following commands are special cased:
      - `?`: Print this help text.
      - `@trace`: Print the current stack trace.
      - `@locals`: Print local variables.
      - `@exfiltrate`: Save all local variables into the store.
      - `@toggle`: Toggle infiltrating at this `@infiltrate` spot (clear all with `Infiltrator.clear_disabled!()`).
      - `@continue`: Continue to the next infiltration point or exit (shortcut: Ctrl-D).
      - `@doc symbol`: Get help for `symbol` (same as in the normal Julia REPL).
      - `@exit`: Stop infiltrating for the remainder of this session and exit (on Julia versions prior to
        1.5 this needs to be manually cleared with `Infiltrator.end_session!()`).
  """)
end

function show_trace(io, trace)
  for (i, frame) in enumerate(trace)
    println(io, "[", i, "] ", frame)
  end
  println(io)
end

function strlimit(str, limit = 30)
  strwidth = 0
  for (i, c) in enumerate(str)
    strwidth += textwidth(c)
    strwidth >= limit && return string(first(str, i - 1), "…")
  end
  return str
end

function show_locals(io, locals)
  for (var, val) in locals
    one_line_show(io, var, val)
  end
  println(io)
end

function one_line_show(io::IO, var, val)
  width = max(displaysize(io)[2], 40)
  prefix = string("- ", strlimit(string(var), width ÷ 3), "::", strlimit(string(typeof(val)), width ÷ 3), " = ")
  println(io, prefix, strlimit(repr(val), max(width - textwidth(prefix), 10)))
end

function is_complete(s)
  ast = Base.parse_input_line(String(take!(copy(LineEdit.buffer(s)))), depwarn=false)
  return !(isa(ast, Expr) && ast.head === :incomplete)
end

const PROMPT = Ref{Any}()

function debugprompt(mod, locals, trace, terminal, repl, nostack = false; file, fileline)
  io = Base.pipe_writer(terminal)

  try
    if isassigned(PROMPT)
      panel = PROMPT[]
      panel.complete = InfiltratorCompletionProvider(mod, locals)
    else
      panel = PROMPT[] = REPL.LineEdit.Prompt("infil> ";
                prompt_prefix = prompt_color,
                prompt_suffix = Base.text_colors[:normal],
                complete = InfiltratorCompletionProvider(mod, locals),
                on_enter = is_complete)
      panel.hist = REPL.REPLHistoryProvider(Dict{Symbol,Any}(:Infiltrator => panel))
      REPL.history_reset_state(panel.hist)
    end
    search_prompt, skeymap = LineEdit.setup_search_keymap(panel.hist)
    search_prompt.complete = REPL.LatexCompletions()

    panel.on_done = (s, buf, ok) -> begin
      if !ok
        LineEdit.transition(s, :abort)
        REPL.LineEdit.reset_state(s)
        return false
      end

      line = String(take!(buf))

      if isempty(line)
        println(io)
        LineEdit.reset_state(s)
        return true
      end

      sline = strip(line)
      if sline == "?"
        show_help(io)
        LineEdit.reset_state(s)
        return true
      elseif sline == "@trace"
        show_trace(io, trace)
        LineEdit.reset_state(s)
        return true
      elseif sline == "@locals"
        show_locals(io, locals)
        LineEdit.reset_state(s)
        return true
      elseif sline == "@exfiltrate"
        n = length(locals)
        println(io, "Exfiltrating $(n) local variable$(n == 1 ? "" : "s") into the store.\n")

        for (k, v) in locals
          try
            Core.eval(getfield(store, :store), Expr(:(=), k, QuoteNode(v)))
          catch err
            println(io, "Assignment to store variable failed.")
            Base.display_error(io, err, catch_backtrace())
          end
        end
        LineEdit.reset_state(s)
        return true
      elseif sline == "@toggle"
        spot = (file, fileline)
        ds = getfield(store, :disabled)
        if spot in ds
          delete!(ds, spot)
          println(io, "Enabled infiltration at this infiltration point.\n")
        else
          push!(ds, spot)
          println(io, "Disabled infiltration at this infiltration point.\n")
        end
        LineEdit.reset_state(s)
        return true
      elseif sline == "@exit"
        setfield!(store, :exiting, true)
        if !REPL_HOOKED[]
          println(io, "Revert the effect of this with `Infiltrator.end_session!()` or you will not be able to enter a new session!")
        end
        LineEdit.transition(s, :abort)
        LineEdit.reset_state(s)
        return true
      elseif sline == "@continue"
        LineEdit.transition(s, :abort)
        LineEdit.reset_state(s)
        return true
      end
      ok = true
      result = nothing

      expr = try
        Base.parse_input_line(line)
      catch err
        if get(io, :color, false)
          println(io, "Parsing error.")
        else
          printstyled(io, "Parsing error.", color = Base.error_color())
        end

        LineEdit.reset_state(s)
        return false
      end

      @static if VERSION >= v"1.2.0-DEV.253"
        try
          result = interpret(io, expr, mod, locals)
        catch err
          ok = false
          result = Base.catch_stack()
          if nostack
            result = map(r -> Any[first(r), []], result)
          else
            result = map(result) do (err, bt)
              return err, crop_backtrace(bt)
            end
          end
        end
        REPL.print_response(repl, (result, !ok), true, true)
      else
        try
          result = interpret(io, expr, mod, locals)
        catch err
          ok = false
          result = (err, nostack ? Any[] : crop_backtrace(catch_backtrace()))
        end
        REPL.print_response(repl, ok ? result : result[1], ok ? nothing : result[2], true, true)
      end
      println(io)
      LineEdit.reset_state(s)

      return true
    end

    panel.keymap_dict = LineEdit.keymap(Dict{Any,Any}[skeymap, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults])

    REPL.run_interface(terminal, REPL.LineEdit.ModalInterface([panel, search_prompt]))
  catch e
    e isa InterruptException || rethrow(e)
  end
end

function get_store_names(s = store)
  m = getfield(s, :store)
  ns = names(m, all = true)
  out = Dict()
  for n in ns
    if isdefined(m, n) && n !== :eval && n !== :include && n !== :anonymous
      out[n] = getfield(m, n)
    end
  end
  out
end

maybe_quote(x) = (isa(x, Expr) || isa(x, Symbol)) ? QuoteNode(x) : x

function interpret(io, expr, mod, locals)
  symbols = merge(get_store_names(), locals)
  Meta.isexpr(expr, :toplevel) && (expr = expr.args[end])
  res = gensym()
  eval_expr = Expr(:let,
    Expr(:block, map(x->Expr(:(=), x...), [(k, maybe_quote(v)) for (k, v) in symbols])...),
    Expr(:block,
      Expr(:(=), res, expr),
      Expr(:tuple, res, Expr(:tuple, [k for (k, v) in locals]...))
    )
  )
  eval_res, res = Core.eval(mod, eval_expr)

  assignment = nothing
  if Meta.isexpr(expr, :(=))
    if expr.args[1] isa Symbol
      assignment = expr.args[1]
    elseif Meta.isexpr(expr.args[1], :call) && expr.args[1].args[1] isa Symbol
      assignment = expr.args[1].args[1]
    end
  elseif Meta.isexpr(expr, :function)
    assignment = expr.args[1].args[1]
  end
  if assignment !== nothing
    if haskey(locals, assignment)
      str = "Cannot assign a new value to local variable `$(assignment)`."
      if get(io, :color, false)
        println(io, str)
      else
        printstyled(io, str, color = Base.error_color())
      end
    else
      try
        Core.eval(getfield(store, :store), Expr(:(=), assignment, QuoteNode(eval_res)))
      catch err
        println(io, "Assignment to store variable failed.")
        Base.display_error(io, err, catch_backtrace())
      end
    end
  end

  eval_res
end

function exit_transform(ex)
  return quote
    out = $(ex)
    $(@__MODULE__).end_session!()
    out
  end
end

# backtraces

function find_first_topelevel_scope(bt::Vector{<:Union{Base.InterpreterIP,Ptr{Cvoid}}})
  for (i, ip) in enumerate(bt)
      st = Base.StackTraces.lookup(ip)
      ind = findfirst(st) do frame
          linfo = frame.linfo
          if linfo isa Core.CodeInfo
              linetable = linfo.linetable
              if isa(linetable, Vector) && length(linetable) ≥ 1
                  lin = first(linetable)
                  if isa(lin, Core.LineInfoNode) && lin.method === Symbol("top-level scope")
                      return true
                  end
              end
          else
              return frame.func === Symbol("top-level scope")
          end
      end
      ind === nothing || return i
  end
  return
end

function crop_backtrace(bt)
  i = find_first_topelevel_scope(bt)
  return bt[1:(i === nothing ? end : i)]
end

# completions

mutable struct InfiltratorCompletionProvider <: REPL.CompletionProvider
  mod::Module
  locals::Dict{Symbol, Any}
end

function LineEdit.complete_line(c::InfiltratorCompletionProvider, s)
  partial = REPL.beforecursor(s.input_buffer)
  full = LineEdit.input_string(s)

  ret, range, should_complete = completions(c, full, partial)
  return ret, partial[range], should_complete
end

function completions(c::InfiltratorCompletionProvider, full, partial)
  # repl backend completions
  comps, range, should_complete = REPL.REPLCompletions.completions(full, lastindex(partial), c.mod)
  ret = map(REPL.REPLCompletions.completion_text, comps)

  # local completions
  prepend!(ret, filter!(v -> startswith(v, partial), vcat(string.(keys(c.locals)), string.(keys(get_store_names())))))

  # Infiltrator commands completions
  commands = ["?", "@trace", "@locals", "@toggle", "@exit", "@continue", "@exfiltrate"]
  prepend!(ret, filter!(c -> startswith(c, partial), commands))

  unique!(ret), range, should_complete
end

end # module


text:
using CatmullRom, Plots

x(t) = cos(t)
y(t,m=4) = sin(t) * (sin(t/2)^m) # m=0..7

n_between = 24
given = collect(range(0.0, length=15, stop=2.0*pi));
xs = x.(given);
ys = y.(given, 5);

points = collect(zip(xs,ys)

crpoints = catmullrom(points), n_between); #
cxs, cys = crpoints;                       # and
cxs, cys = catmullrom(points), n_between); # or

plot(xs, ys, linecolor=:lightgreen, linewidth=9, size=(500,500), xaxis=nothing, yaxis=nothing, legend=nothing)
plot!(cxs, cys, linecolor=:black, linewidth=2, size=(500,500), xaxis=nothing, yaxis=nothing, legend=nothing)

crpoints = catmullrom_by_arclength(collect(zip(xs,ys)), arcpoints_min=6, arcpoints_max=16);
cxs, cys = crpoints;

plot(xs, ys, linecolor=:lightgreen, linewidth=10, size=(500,500), xaxis=nothing, yaxis=nothing, legend=nothing)
plot!(cxs, cys, linecolor=:black, linewidth=2, size=(500,500), xaxis=nothing, yaxis=nothing, legend=nothing)


text:
using FixedEffectModels

tests = ["reg.jl",
		 "RegressionResult.jl",
		 "partial_out.jl"
		 ]

println("Running tests:")

for test in tests
	try
		include(test)
		println("\t\033[1m\033[32mPASSED\033[0m: $(test)")
	 catch e
	 	println("\t\033[1m\033[31mFAILED\033[0m: $(test)")
	 	showerror(stdout, e, backtrace())
	 	rethrow(e)
	 end
end

using Test
@test FixedEffectModels.pinvertible(Symmetric([1.0 1.0; 1.0 1.0])) ≈ [1.0 1.0; 1.0 1.0]