text:
"""
    Probability
Represents a probability that something (typically mutation) will happen.

Possible to specify RNG implementation. If not specified, `rng_default` will be used.
"""
struct Probability{P<:Real, R<:AbstractRNG}
    p::P
    rng::R
    function Probability(p::P, rng::R) where {P, R}
         @assert 0 <= p <= 1
         return new{P, R}(p, rng)
    end
end
Probability(p::Real) = Probability(p, rng_default)
Probability(p::Integer) = Probability(p, rng_default)
Probability(p::Integer, rng) = p == 1 ? Probability(1.0, rng) : Probability(p / 100.0, rng)

"""
    apply(p::Probability)
    apply(p::Number)

Return true with a probability of ´p.p´ (subject to `p.rng` behaviour).
"""
apply(p::Probability) = rand(p.rng) < p.p
apply(p::Real) = apply(Probability(p))

"""
    apply(f, p::Probability, or = () -> nothing)
    apply(f, p::Real, or = () -> nothing)

Call `f` with probability `p.p` (subject to `p.rng` behaviour). If `f` is not called then call `or`.
"""
apply(f, p::Probability, or = () -> nothing) =  apply(p) ? f() : or()
apply(f, p::Real, or = () -> nothing) = apply(p) ? f() : or()

"""
    MutationShield <: DecoratingTrait
    MutationShield(t, allowed...)

Shields its associated vertex from being selected for mutation.

Any types in `allowed` will be allowed to mutate the vertex if supplied when calling `allow_mutation`.

Note that vertex might still be modified if an adjacent vertex is mutated in a way which propagates to a shielded vertex.
"""
struct MutationShield{T<:MutationTrait, S} <:DecoratingTrait
    t::T
    allowed::S
    MutationShield(t::T, allowed::Tuple) where T = new{T, typeof(allowed)}(t, allowed)
end
MutationShield(t, allowed...) = MutationShield(t, allowed)

NaiveNASlib.base(t::MutationShield) = t.t
allow_mutation(v::AbstractVertex, ms...) = allow_mutation(trait(v), ms...)
allow_mutation(t::DecoratingTrait, ms...) = allow_mutation(base(t), ms...)
allow_mutation(::MutationTrait, ms...) = true
allow_mutation(::Immutable, ms...) = false
allow_mutation(t::MutationShield, ms...) = !isempty(ms) && all(mt -> any(amt -> mt <: amt, t.allowed), typeof.(mutationleaves(ms)))

@functor MutationShield

"""
    AbstractVertexSelection

Abstract type for determining how to select vertices from a `CompGraph` or an array of vertices.
"""
abstract type AbstractVertexSelection end

"""
    AllVertices

Select all vertices in `g`.
"""
struct AllVertices <:AbstractVertexSelection end
select(::AllVertices, g::CompGraph, ms...) = vertices(g)
select(::AllVertices, vs::AbstractArray, ms...) = vs

"""
    FilterMutationAllowed

Filters out only the vertices for which mutation is allowed from another selection.
"""
struct FilterMutationAllowed{S} <:AbstractVertexSelection
    s::S
end
FilterMutationAllowed() = FilterMutationAllowed(AllVertices())
select(s::FilterMutationAllowed, x, ms...) = filter(v -> allow_mutation(v, ms...), select(s.s, x, ms...))

"""
    SelectWithMutation{S, M} <: AbstractVertexSelection
    SelectWithMutation(m::M)
    SelectWithMutation(s::S, m::M)

Adds `m` to the list of `ms` when selecting vertices using `s`.

Useful when calling `select` from inside a function without the context of an `AbstractMutation`.
"""
struct SelectWithMutation{S, M} <: AbstractVertexSelection
    s::S
    m::M
end
SelectWithMutation(m) = SelectWithMutation(FilterMutationAllowed(), m)
select(s::SelectWithMutation, x, ms...) = select(s.s, x, (s.m, ms...))

"""
    ApplyIf <: DecoratingTrait
    ApplyIf(predicate::Function, apply::Function, base::MutationTrait)

Enables calling `apply(v)` for an `AbstractVertex v` which has this trait if 'predicate(v) == true'.

Motivating use case is to have a way to remove vertices which have ended up as noops, e.g. element wise and concatenation vertices with a single input or identity activation functions.
"""
struct ApplyIf{P, A, M<:MutationTrait} <: DecoratingTrait
    predicate::P
    apply::A
    base::M
end
RemoveIfSingleInput(t) = ApplyIf(v -> length(inputs(v)) == 1, remove!, t)
NaiveNASlib.base(t::ApplyIf) = t.base

check_apply(g::CompGraph) = foreach(check_apply, vertices(g))
check_apply(v::AbstractVertex) = check_apply(trait(v), v)
check_apply(t::DecoratingTrait, v) = check_apply(base(t), v)
check_apply(t::ApplyIf, v) = t.predicate(v) && t.apply(v)
function check_apply(t, v) end

@functor ApplyIf

"""
    PersistentArray{T, N} <: AbstractArray{T, N}
    PersistentArray(savedir::String, nr::Integer, generator;suffix=".jls")
    PersistentArray(savedir::String, suffix::String, data::Array)

Simple persistent array. Can be created from serialized data and can be asked to persist its elements using [`persist`](@ref).

Note that once initialized, the array is not backed by the serialized data. Adding/deleting files is not reflected in data and vice versa.
"""
struct PersistentArray{T, N} <: AbstractArray{T, N}
    savedir::String
    suffix::String
    data::Array{T,N}
end
function PersistentArray(savedir::String, nr::Integer, generator;suffix=".jls")
    data = map(1:nr) do i
        filename = joinpath(savedir, "$i$suffix")
        isfile(filename) && return deserialize(filename)
        return generator(i)
    end
    return PersistentArray(savedir, suffix, data)
end
"""
    persist(a::PersistentArray)

Serializes the elements of `a`, one file per element.
"""
function persist(a::PersistentArray)
    mkpath(a.savedir)
    for (i, v) in enumerate(a)
        serialize(filename(a, i), v)
    end
end
filename(a::PersistentArray, i::Int) = joinpath(a.savedir, "$i$(a.suffix)")
function Base.rm(a::PersistentArray; force=true, recursive=true)
    foreach(i -> rm(a, i; force=force, recursive=recursive), 1:length(a))
    if readdir(a.savedir) |> isempty
        rm(a.savedir; force=force, recursive=recursive)
    end
end
Base.rm(a::PersistentArray, i::Int; force=false, recursive=true) = rm(filename(a,i), force=force, recursive=recursive)

Base.size(a::PersistentArray) = size(a.data)
Base.getindex(a::PersistentArray, i::Int) = getindex(a.data, i)
Base.getindex(a::PersistentArray, I::Vararg{Int, N}) where N = getindex(a.data, I...)
Base.setindex!(a::PersistentArray, v, i::Int) = setindex!(a.data, v, i)
Base.setindex!(a::PersistentArray, v, I::Vararg{Int, N}) where N = setindex!(a.data, v, I...)
Base.similar(a::PersistentArray, t::Type{S}, dims::Dims) where S = PersistentArray(a.savedir, a.suffix, similar(a.data,t, dims))

"""
    BoundedRandomWalk{T <: Real, R <: Function}
    BoundedRandomWalk(lb, ub, rfun =  (x...) -> 0.2randn(rng_default))

Generates steps for a random walk with bounds `[lb, ub]`.

Main use case is with learning rate mutation to prevent it from accidentally drifting to some catastropically high value.

# Examples
```julia-repl
julia> import NaiveGAflux: BoundedRandomWalk

julia> using Random

julia> rng = MersenneTwister(0);

julia> brw = BoundedRandomWalk(-2.34, 3.45, () -> 10randn(rng));

julia> extrema(cumsum([brw() for i in 1:10000]))
(-2.3400000000000007, 3.4500000000000153)
```
"""
struct BoundedRandomWalk{T <: Real, R <: Function}
    lb::T
    ub::T
    state::Ref{T}
    rfun::R
end
BoundedRandomWalk(lb::T,ub::T, rfun = (x...) -> 0.2randn(rng_default)) where T = BoundedRandomWalk(lb,ub, Ref(zero(ub)), rfun)

function(r::BoundedRandomWalk)(x...)
    y =  ((r.ub - r.lb) * clamp(r.rfun(x...), -1, 1) + (r.ub + r.lb)) / 2 - r.state[]
    r.state[] += y
    return y
 end

 """
    FluxOptimizer

Alias for Flux.Optimise.AbstractOptimiser
"""
const FluxOptimizer = Flux.Optimise.AbstractOptimiser

"""
    ShieldedOpt{O} <: Flux.Optimise.AbstractOptimiser
    ShieldedOpt(o)

Shields `o` from mutation by `OptimizerMutation`.
"""
struct ShieldedOpt{O<:FluxOptimizer} <: FluxOptimizer
    opt::O
end
Flux.Optimise.apply!(o::ShieldedOpt, args...) = Flux.Optimise.apply!(o.opt, args...)

"""
    mergeopts(t::Type{T}, os...) where T
    mergeopts(os::T...)

Merge all optimizers of type `T` in `os` into one optimizer of type `T`.

Defaults to `T(prod(learningrate.(os)))`.
"""
function mergeopts(t::Type{T}, os...) where T
    merged = mergeopts(filter(o -> isa(o, T), os)...)
    return vcat(filter(o -> !isa(o, T), os)..., merged)
end
mergeopts() = []
mergeopts(t::Type{T}, os::T...) where T = [mergeopts(os...)]
mergeopts(os...) = first(@set os[1].eta = (prod(learningrate.(os))))
mergeopts(os::ShieldedOpt{T}...) where T = ShieldedOpt(mergeopts(map(o -> o.opt, os)...))
mergeopts(os::WeightDecay...) = WeightDecay(mapreduce(o -> o.wd, *, os))

"""
    optmap(fopt, x, felse=identity)
    optmap(fopt, felse=identity)

Return `fopt(x)` if `x` is an optimizer, else return `felse(x)`.

Call without x to return `x -> optmap(fopt, x, felse)`
"""
optmap(fopt, felse=identity) = x -> optmap(fopt, x, felse)
optmap(fopt, x, felse) = felse(x)
optmap(fopt, o::FluxOptimizer, felse) = fopt(o)

function clearstate!(s) end
clearstate!(s::AbstractDict) = empty!(s)

cleanopt!(o::T) where T = (foreach(fn -> clearstate!(getfield(o, fn)), fieldnames(T)); return o)
cleanopt!(o::ShieldedOpt) = (cleanopt!(o.opt); return o)
cleanopt!(o::Flux.Optimiser) = (foreach(cleanopt!, o.os); return o)

"""
    Singleton{T}
    Singleton(val::T)

Wrapper for `val` which prevents it from being copied if the wrapping singleton is copied using `copy` and `deepcopy`.

Also makes sure that only one unique copy of `val` is created when deserializing a serialized `Singleton`.
"""
struct Singleton{T}
    val::T
    function Singleton(val::T) where T
        s = new{T}(val)
        singletons[val] = s
        return s
    end
end
val(s::Singleton) = s.val

const singletons = WeakKeyDict()
function Serialization.deserialize(s::AbstractSerializer, ::Type{Singleton{T}}) where T
    val = deserialize(s)
    return get!(()-> Singleton(val), singletons, val)
end

Base.deepcopy_internal(s::Singleton, stackdict::IdDict) = s
Base.copy(s::Singleton) = s

Base.iterate(s::Singleton, state...) = iterate(val(s), state)

Base.length(s::Singleton) = length(val(s))
Base.size(s::Singleton) = size(val(s))

Base.IteratorEltype(s::Singleton) = Base.IteratorEltype(val(s))
Base.IteratorSize(s::Singleton) = Base.IteratorSize(val(s))

"""
     PrefixLogger{L<:AbstractLogger, S} <: AbstractLogger
     PrefixLogger(wrappedLogger::L, prefix::S)

Add `prefix` to messages and forwards to `wrappedLogger`.
"""
struct PrefixLogger{L<:AbstractLogger, S} <: AbstractLogger
      wrapped::L
      prefix::S
end
PrefixLogger(prefix::String) = PrefixLogger(current_logger(), prefix)

Logging.min_enabled_level(l::PrefixLogger) = Logging.min_enabled_level(l.wrapped)
Logging.shouldlog(l::PrefixLogger, args...) =  Logging.shouldlog(l.wrapped, args...)
Logging.handle_message(l::PrefixLogger, level, message, args...; kwargs...) = Logging.handle_message(l.wrapped, level, l.prefix * message, args...; kwargs...)

"""
    GlobalPool{PT}
    GlobalPool(PT)

Global pool combined with `Flux.flatten`. Pool type is determined by `PT`, typically `MaxPool` or `MeanPool`.

Useful because in the context of mutations, one most likely want to keep them as one single vertex to avoid mutations happening between them.
"""
struct GlobalPool{PT} end
GlobalPool(PT) = GlobalPool{PT}()

# About 50% faster on GPU to create a MeanPool and use it compared to dropdims(mean(x, dims=[1:2]), dims=(1,2)). CBA to figure out why...
(::GlobalPool{PT})(x::AbstractArray{<:Any, N}) where {N, PT} = Flux.flatten(PT(size(x)[1:N-2])(x))

NaiveNASflux.layertype(gp::GlobalPool) = gp
NaiveNASflux.layer(gp::GlobalPool) = gp

"""
    ninputs(model)

Return the number of model inputs.
"""
ninputs(cg::CompGraph) = length(cg.inputs)


text:
# PS2  Dayuan Xie
# using Pkg
# Pkg.add("Optim")
# Pkg.add("HTTP")
# Pkg.add("GLM")

using Optim
using HTTP
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using CSV
using FreqTables

function q6()

# Question 1
f(x) = -x[1]^4-10x[1]^3-2x[1]^2-3x[1]-2
minusf(x) = x[1]^4+10x[1]^3+2x[1]^2+3x[1]+2
startval = rand(1)
result = optimize(minusf, startval, BFGS())
println(result)

# Question 2
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS1-julia-intro/nlsw88.csv"
df = CSV.read(HTTP.get(url).body)
X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.married.==1

function ols(beta, X, y)
    ssr = (y.-X*beta)'*(y.-X*beta)
    return ssr
end

beta_hat_ols = optimize(b -> ols(b, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
println(beta_hat_ols.minimizer)

bols = inv(X'*X)*X'*y
df.white = df.race.==1
bols_lm = lm(@formula(married ~ age + white + collgrad), df)
println(bols_lm)

# Question 3
function logit(alpha, X, d)
    loglike=d'*X*alpha-sum(log.(ones(size(X,1),1)+exp.(X*alpha)))
    return loglike
end
alpha_hat_logit = optimize(a -> -(logit(a, X, y)), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
println(alpha_hat_logit.minimizer)

# question 4
alpha_glm= glm(@formula(married ~ age + white + collgrad), df, Binomial(), LogitLink())
println(alpha_glm)

# Question 5
freqtable(df, :occupation) # note small number of obs in some occupations
df = dropmissing(df, :occupation)
df[df.occupation.==8 ,:occupation] .= 7
df[df.occupation.==9 ,:occupation] .= 7
df[df.occupation.==10,:occupation] .= 7
df[df.occupation.==11,:occupation] .= 7
df[df.occupation.==12,:occupation] .= 7
df[df.occupation.==13,:occupation] .= 7
freqtable(df, :occupation) # problem solved

X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.occupation

y1=y.==1
y2=y.==2
y3=y.==3
y4=y.==4
y5=y.==5
y6=y.==6
y7=y.==7

Y = cat(y1,y2,y3,y4,y5,y6,dims=2)

function csum(A)
    Colsum = []
    for i=1:size(A,1)
        sum(A[i,:])
        push!(Colsum,sum(A[i,:]))
    end
        return Colsum
end

function mlogit(alpha, X, d)
    loglike=tr(d*(X*alpha)')-sum(log.(csum(exp.(X*alpha))))
    return loglike
end

beta_hat_mlogit1 = optimize(b -> -mlogit(b, X, Y), zeros(size(X,2),6), LBFGS(), Optim.Options(g_tol=10-5, iterations=100_000, show_trace=true))
println(beta_hat_mlogit1.minimizer)

beta_hat_mlogit2 = optimize(b -> -mlogit(b, X, Y), rand(size(X,2),6), LBFGS(), Optim.Options(g_tol=10-5, iterations=100_000, show_trace=true))
println(beta_hat_mlogit2.minimizer)

beta_hat_mlogit3 = optimize(b -> -mlogit(b, X, Y), 2*rand(size(X,2),6).-1, LBFGS(), Optim.Options(g_tol=10-5, iterations=100_000, show_trace=true))
println(beta_hat_mlogit3.minimizer)

end

# Question 6
q6()

ext: jl


text:
#
# Defines AstrodynamicalCalculations methods for
# OrbitalStates types
#

using OrbitalStates
import OrbitalStates: CartesianStateVector

export massparameter, massparameters, normalized_massparameter
export primary_massparameter, secondary_massparameter

include("States.jl")
include("Systems.jl")


text:
"""
An implementation of simple method to calculate Ewald energy
"""
function calc_ewald( pw::PWGrid, Xpos, Sf; sigma=0.25 )
    #
    Npoints = pw.Npoints
    Ω  = pw.Ω
    r  = pw.r
    Ns = pw.Ns
    G2 = pw.G2
    #
    # Generate array of distances
    center = sum(pw.LatVecs,dims=2)/2
    dr = gen_dr( r, center )
    #
    # Generate charge density
    rho = gen_rho( Ns, dr, sigma, Sf )
    intrho = sum(rho)*Ω/Npoints
    @printf("sigma, int_rho: %10.5f %18.10f\n", sigma, intrho)

    #
    # Solve Poisson equation and calculate Hartree energy
    #
    ctmp = 4.0*pi*R_to_G( Ns, rho )
    ctmp[1] = 0.0
    for ip = 2:Npoints
      ctmp[ip] = ctmp[ip] / G2[ip]
    end
    phi = real( G_to_R( Ns, ctmp ) )
    Ehartree = 0.5*dot( phi, rho ) * Ω/Npoints
    #
    Eself = 1.0/(2*sqrt(pi))*(1.0/sigma)*size(Xpos,2)
    E_nn = Ehartree - Eself
    #@printf("Ehartree, Eself, E_nn = %20.16f %20.16f %20.16f\n", Ehartree, Eself, E_nn)
    return E_nn
end


text:
# NOTE: this file keeps patches for julia itself to make JET.jl test keep passing

@static if VERSION < v"1.7.0-DEV.0"
    # take in https://github.com/JuliaLang/julia/pull/38832
    @eval Base function _split(str::AbstractString, splitter::F, limit::Integer, keepempty::Bool, strs::Vector) where F
        # Forcing specialization on `splitter` improves performance (roughly 30% decrease in runtime)
        # and prevents a major invalidation risk (1550 MethodInstances)
        i = 1 # firstindex(str)
        n = lastindex(str)::Int
        r = findfirst(splitter,str)::Union{Nothing,Int,UnitRange{Int}}
        if r !== nothing
            j, k = first(r), nextind(str,last(r))::Int
            while 0 < j <= n && length(strs) != limit-1
                if i < k
                    if keepempty || i < j
                        push!(strs, @inbounds SubString(str,i,prevind(str,j)::Int))
                    end
                    i = k
                end
                (k <= j) && (k = nextind(str,j)::Int)
                r = findnext(splitter,str,k)::Union{Nothing,Int,UnitRange{Int}}
                r === nothing && break
                j, k = first(r), nextind(str,last(r))::Int
            end
        end
        if keepempty || i <= ncodeunits(str)::Int
            push!(strs, @inbounds SubString(str,i))
        end
        return strs
    end
end


text:
module GridapGeosciencesTests

include("sequential/runtests.jl")

include("mpi/runtests.jl")

end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
struct Chart{R, T}
    num_tokens::Int
    active::Dict{T, Vector{Vector{ActiveArc{R}}}} # organized by next needed constituent then by stop
    passive::Dict{T, Vector{Vector{PassiveArc{R}}}} # organized by head then by start
end

Chart{R, T}(num_tokens::Integer) where {R, T} =
    Chart(num_tokens, Dict{T, Vector{Vector{ActiveArc{R}}}}(),
                Dict{T, Vector{Vector{PassiveArc{R}}}}())

num_tokens(chart::Chart) = chart.num_tokens

function _active_storage(chart::Chart{R, T}, next_needed::T, stop::Integer) where {R, T}
    v = get!(chart.active, next_needed) do
        [Vector{ActiveArc{R}}() for _ in 0:num_tokens(chart)]
    end
    v[stop + 1]
end

function _passive_storage(chart::Chart{R, T}, head::T, start::Integer) where {R, T}
    v = get!(chart.passive, head) do
        [Vector{PassiveArc{R}}() for _ in 0:num_tokens(chart)]
    end
    v[start + 1]
end

storage(chart::Chart, arc::ActiveArc) = _active_storage(chart, next_needed(arc), stop(arc))
storage(chart::Chart, arc::PassiveArc) = _passive_storage(chart, head(arc), start(arc))

mates(chart::Chart, candidate::ActiveArc) = _passive_storage(chart, next_needed(candidate), stop(candidate))
mates(chart::Chart, candidate::PassiveArc) = _active_storage(chart, head(candidate), start(candidate))

function Base.push!(chart::Chart, arc::AbstractArc)
    push!(storage(chart, arc), arc)
end

Base.in(arc::AbstractArc, chart::Chart) = arc ∈ storage(chart, arc)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dashbio_nglmoleculeviewer

"""
    dashbio_nglmoleculeviewer(;kwargs...)

A NglMoleculeViewer component.
The NglMoleculeViewer is used to render schematic diagrams
of biomolecules in ribbon-structure representations.
Read more about the component here:
https://github.com/IvoLeist/dash_ngl
Read more about the used WebGL protein viewer here:
https://github.com/arose/ngl
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks.
The ID needs to be unique across all of the components in an app.
- `data` (optional): The data (in JSON format) that will be used to display the molecule
filename: name of the used pdb/cif file
ext: file extensions (pdb or cif)
selectedValue: pdbString
chain: ALL if the whole molecule shoud be displayed, e.g. A for showing only chain A
aaRange: ALL if the whole molecule should be displayed, e.g. 1:50 for showing only 50 atoms
color: chain color
chosen.atoms: string of the chosen Atoms, e.g. 50,100,150
              --> chosen eatoms changed to colored 'ball'
chosen.residues: string of the chosen residues, e.g. 50,100,150
                 --> C alpha of chosen residue changed to colored 'ball'
config.input: content of the pdb file
config.type: format of config.input
uploaded: bool if file from local storage (false) or uploaded by user (true)
resetView: bool if the selection did not change but the view should be resettet (true). data has the following type: Array of lists containing elements 'filename', 'ext', 'selectedValue', 'chain', 'aaRange', 'color', 'chosen', 'config', 'uploaded', 'resetView'.
Those elements have the following types:
  - `filename` (String; required)
  - `ext` (String; optional)
  - `selectedValue` (String; required)
  - `chain` (String; required)
  - `aaRange` (String; required)
  - `color` (String; required)
  - `chosen` (optional): . chosen has the following type: lists containing elements 'residues', 'atoms'.
Those elements have the following types:
  - `residues` (String; required)
  - `atoms` (String; required)
  - `config` (optional): . config has the following type: lists containing elements 'input', 'type'.
Those elements have the following types:
  - `input` (String; required)
  - `type` (String; required)
  - `uploaded` (Bool; required)
  - `resetView` (Bool; required)s
- `downloadImage` (Bool; optional): flag if download image was selected
- `height` (String | Real; optional): The height (in px or as a number) of the container
in which the molecules will be displayed.
- `imageParameters` (optional): Parameters (in JSON format) for exporting the image. imageParameters has the following type: lists containing elements 'antialias', 'transparent', 'trim', 'defaultFilename'.
Those elements have the following types:
  - `antialias` (Bool; optional)
  - `transparent` (Bool; optional)
  - `trim` (Bool; optional)
  - `defaultFilename` (String; optional)
- `molStyles` (optional): The data (in JSON format) that will be used to style the displayed molecule
representations: one or multiple selected molecule representation
 - Possible molecule styles:
   'backbone,'ball+stick','cartoon', 'hyperball','licorice','line',
   'ribbon',''rope','spacefill','surface','trace','tube'
 - Possible additional representations:
   'axes','axes+box','helixorient','unitcell'
chosenAtomsColor: color of the 'ball+stick' representation of the chosen atoms
chosenAtomsRadius: radius of the 'ball+stick' representation of the chosen atoms
molSpacingXaxis: distance on the xAxis between each molecule. molStyles has the following type: lists containing elements 'representations', 'chosenAtomsColor', 'chosenAtomsRadius', 'molSpacingXaxis', 'sideByside'.
Those elements have the following types:
  - `representations` (Array of Strings; optional)
  - `chosenAtomsColor` (String; required)
  - `chosenAtomsRadius` (Real; required)
  - `molSpacingXaxis` (Real; required)
  - `sideByside` (Bool; required)
- `pdbString` (String; optional): Variable which defines how many molecules should be shown and/or which chain
The following format needs to be used:
pdbID1.chain:start-end@atom1,atom2_pdbID2.chain:start-end
. indicates that only one chain should be shown
: indicates that a specific amino acids range should be shown (e.g. 1-50)
@ indicates that chosen atoms should be highlighted (e.g. @50,100,150)
 _ indicates that more than one protein should be shown
- `stageParameters` (optional): Parameters (in JSON format) for the stage object of ngl.
Currently implemented are render quality, background color and camera type
quality: auto, low, medium, high (default: auto)
backgroundColor: white / black (default: white)
cameraType: perspective / orthographic (default: perspective). stageParameters has the following type: lists containing elements 'quality', 'backgroundColor', 'cameraType'.
Those elements have the following types:
  - `quality` (String; optional)
  - `backgroundColor` (String; optional)
  - `cameraType` (String; optional)
- `width` (String | Real; optional): The width (in px or as a number) of the container
in which the molecules will be displayed.
"""
function dashbio_nglmoleculeviewer(; kwargs...)
        available_props = Symbol[:id, :data, :downloadImage, :height, :imageParameters, :molStyles, :pdbString, :stageParameters, :width]
        wild_props = Symbol[]
        return Component("dashbio_nglmoleculeviewer", "NglMoleculeViewer", "dash_bio", available_props, wild_props; kwargs...)
end


text:
using AbstractPlotting
using MakieGallery
using Test
using GLMakie

# Download reference images from master
MakieGallery.current_ref_version[] = "master"

const MINIMAL = get(ENV, "ABSTRACTPLOTTING_MINIMAL", "false")

# does this machine have a OPENGL?
const OPENGL = GLMakie.WORKING_OPENGL
OPENGL || (MINIMAL = "true")

@show OPENGL MINIMAL

OPENGL && begin @info "OpenGL detected"; using GLMakie end
OPENGL || @warn "No OpenGL detected!  Software tests only."

@info "Running conversion tests"
include("conversions.jl")

@info "Running quaternion tests"
include("quaternions.jl")

@info "Running projection tests"
include("projection_math.jl")

@info "Running shorthand tests"
include("shorthands.jl")

@info "Running basic Scene test"
@testset "basic functionality" begin
    scene = scatter(rand(4))
    @test scene[Axis].ticks.title_gap[] == 3
    scene[Axis].ticks.title_gap = 4
    @test scene[Axis].ticks.title_gap[] == 4
    @test scene[Axis].tickmarks.length[] == (3, 3)
end

if MINIMAL == "false"

    @info "Minimal was false; running full tests!"

    # Load all entries in the database

    database = MakieGallery.load_database()

    ## Exceptions are made on basis of title and not index,
    # because index may change as MakieGallery changes.

    # we found which ones are the slowest, and kicked those out!
    slow_examples = [
        "Animated time series",
        "Animation",
        "Lots of Heatmaps",
        "Chess Game",
        "Line changing colour",
        "Line changing colour with Observables",
        "Colormap collection",
        "Record Video",
        "Animated surface and wireframe",
        "Moire",
        "Line GIF",
        "Electrostatic repulsion",
        "pong",
        "pulsing marker",
        "Travelling wave",
        "Axis theming",
        "Legend",
        "Color Legend",
        "DifferentialEquations path animation",
        "Interactive Differential Equation",
        "Spacecraft from a galaxy far, far away",
        "WorldClim visualization",
        "Image on Geometry (Moon)",
        "Image on Geometry (Earth)",
        "Interaction with Mouse",
        "Air Particulates",
        "Normals of a Cat",
    ]

    # diffeq is also slow, as are analysis heatmaps.  Colormap collection likes to
    # fail a lot.
    filter!(MakieGallery.database) do entry
        !("diffeq" in entry.tags) &&
        !(entry.unique_name in (:analysis, :colormap_collection, :lots_of_heatmaps)) &&
        !(entry.title in slow_examples)
     end

    # Download is broken on CI
    if get(ENV, "CI", "false") == "true"
        printstyled("CI detected\n"; bold = true, color = :yellow)
        println("Filtering out examples which download")
        filter!(entry-> !("download" in entry.tags), database)
    end

else
    # Load only short tests
    database = MakieGallery.load_tests()
end

# one last time to make sure
filter!(database) do entry
    !("diffeq" in entry.tags) &&
    !(entry.unique_name in (:analysis, :colormap_collection, :lots_of_heatmaps, :window_resizing, :aspect_ratios_stretching_circles))
end

# Here, we specialize on two cases.
# If there is no opengl, then we have to run software-only tests, i.e., eval the
# examples and test that nothing errors.
# If there is opengl, then we use MakieGallery's test protocol.

if !OPENGL # run software only tests...
    # Make sure we don't include Makie in the usings
    empty!(MakieGallery.plotting_backends)
    push!(MakieGallery.plotting_backends, "AbstractPlotting")

    @test AbstractPlotting.current_backend[] isa Missing # should we change this, so e.g. CairoMakie works?

    @info "Starting minimal software tests"
    filter!(database) do example
        !("record" in example.tags)
    end

    @testset "Gallery short tests" begin
        # iterate over database
        @testset "$(database[i].title) (#$i)" for i in 1:length(database)

            printstyled("Running "; bold = true, color = :blue)
            print(database[i].title * "\n(index $i)\n")

            # evaluate the entry
            try
                MakieGallery.eval_example(database[i])
                @test true
            catch e
                # THis is ok, since we try to record something, which we can't
                # without backend
                if e isa LoadError && e.error isa MethodError && (
                        e.error.f == AbstractPlotting.backend_display ||
                        e.error.f == AbstractPlotting.format2mime
                    )
                    @test true
                else
                    @test rethrow(e)
                end
            end
            @debug AbstractPlotting.current_backend[] isa Missing
       end
    end

else # full MakieGallery comparisons here
    @info("Running full tests - artifacts will be stored!")

    tested_diff_path = joinpath(@__DIR__, "tested_different")
    test_record_path = joinpath(@__DIR__, "test_recordings")

    rm(tested_diff_path, force = true, recursive = true)
    mkpath(tested_diff_path)

    rm(test_record_path, force = true, recursive = true)
    mkpath(test_record_path)

    examples = MakieGallery.record_examples(test_record_path)

    @test length(examples) == length(database)

    printstyled("Running ", color = :green, bold = true)
    println("visual regression tests")
    MakieGallery.run_comparison(test_record_path, tested_diff_path)
end


text:
module TSXPlayground

using Random

include("utils.jl")
include("TSX.jl")
include("UnsafeAtomics.jl")

using .TSX
using .UnsafeAtomics: unordered, monotonic, acquire, release, acq_rel, seq_cst

include("random_increments.jl")

end


text:
using FixedPointSinCosApproximations
using Test

include("second_order.jl")
include("third_order.jl")
include("fourth_order.jl")
include("fifth_order.jl")
include("sixth_order.jl")
include("approximations.jl")


text:
using PlutoSliderServer
ENV["JULIA_DEBUG"] = PlutoSliderServer

test_dir = tempname(cleanup=false)
cp(@__DIR__, test_dir)

try
    # open the folder on macos:
    run(`open $(test_dir)`)
catch
end

port = 2341

# We want to use our local development version of Pluto for the frontend assets. This allows us to change and test the Pluto assets without having to restart the slider server.
cdn = "http://127.0.0.1:$(port)/pluto_asset/"
# cdn = nothing

PlutoSliderServer.run_directory(
    test_dir;
    SliderServer_port=port,
    SliderServer_host="127.0.0.1",
    SliderServer_watch_dir=true,
    Export_pluto_cdn_root=cdn,
)


text:
using .ForestFire

using GLMakie

# PARAMETERS
nx = 1000;
ny = 1000;
density = 0.5;
transmission_prob = 0.8;

# Initialize forest
forest = Observable(createforest(nx, ny, density));
forest[][div(nx, 2), div(ny, 2)] = 2;

# Initialize Makie scene
# fig = Figure()
# ax = Axis(fig[1, 1])
fig, ax = plotforest(forest, nx, ny)

# Interactive run button
run = Button(fig[2, 1]; label = "run", tellwidth=false)

isrunning = Observable(false)
on(run.clicks) do clicks; isrunning[] = !isrunning[]; end
on(run.clicks) do clicks
    @async while isrunning[]
        isopen(fig.scene) || break
        forest[] = burn(forest[], transmission_prob)
        sleep(0.1)
    end
end

# Start fires where the user clicks
ax1 = content(fig[1, 1])
Makie.deactivate_interaction!(ax1, :rectanglezoom)

spoint = select_point(ax1.scene)

on(spoint) do coord
    x, y = coord
    i = trunc(Int, x)
    j = trunc(Int, y)
    forest[][i, j] = 2
    forest[] = forest[]
end

max_stars_repo_name: TheSquareRoot/ForestFire

max_issues_repo_name: TheSquareRoot/ForestFire

max_forks_repo_name: TheSquareRoot/ForestFire


text:
struct Scale{T, N} <: Bijector{N}
    a::T
end

Scale(a::T; dim::Val{D} = Val(0)) where {T<:Real, D} = Scale{T, D}(a)
Scale(a::A; dim::Val{D} = Val(N)) where {T, D, N, A<:AbstractArray{T, N}} = Scale{A, D}(a)

(b::Scale)(x) = b.a .* x
(b::Scale{<:Real})(x::AbstractArray) = b.a .* x
(b::Scale{<:AbstractMatrix})(x::AbstractArray) = b.a * x
(b::Scale{<:AbstractVector{<:Real}, 2})(x::AbstractMatrix{<:Real}) = b.a .* x

inv(b::Scale{T, D}) where {T, D} = Scale(inv(b.a); dim = Val(D))
inv(b::Scale{<:AbstractVector, D}) where {D} = Scale(inv.(b.a); dim = Val(D))

# We're going to implement custom adjoint for this
logabsdetjac(b::Scale{T, N}, x) where {T, N} = _logabsdetjac_scale(b.a, x, Val(N))

_logabsdetjac_scale(a::Real, x::Real, ::Val{0}) = log(abs(a))
_logabsdetjac_scale(a::Real, x::AbstractVector, ::Val{0}) = fill(log(abs(a)), length(x))
_logabsdetjac_scale(a::Real, x::AbstractVector, ::Val{1}) = log(abs(a)) * length(x)
_logabsdetjac_scale(a::Real, x::AbstractMatrix, ::Val{1}) = fill(log(abs(a)) * size(x, 1), size(x, 2))
_logabsdetjac_scale(a::AbstractVector, x::AbstractVector, ::Val{1}) = sum(log.(abs.(a)))
_logabsdetjac_scale(a::AbstractVector, x::AbstractMatrix, ::Val{1}) = fill(sum(log.(abs.(a))), size(x, 2))
_logabsdetjac_scale(a::AbstractMatrix, x::AbstractVector, ::Val{1}) = log(abs(det(a)))
_logabsdetjac_scale(a::AbstractMatrix, x::AbstractMatrix{T}, ::Val{1}) where {T} = log(abs(det(a))) * ones(T, size(x, 2))

# Adjoints for 0-dim and 1-dim `Scale` using `Real`
function _logabsdetjac_scale(a::TrackedReal, x::Real, ::Val{0})
    return track(_logabsdetjac_scale, a, data(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::Real, ::Val{0})
    return _logabsdetjac_scale(data(a), data(x), Val(0)), Δ -> (inv(data(a)) .* Δ, nothing, nothing)
end

# Need to treat `AbstractVector` and `AbstractMatrix` separately due to ambiguity errors
function _logabsdetjac_scale(a::TrackedReal, x::AbstractVector, ::Val{0})
    return track(_logabsdetjac_scale, a, data(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::AbstractVector, ::Val{0})
    da = data(a)
    J = fill(inv.(da), length(x))
    return _logabsdetjac_scale(da, data(x), Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end

function _logabsdetjac_scale(a::TrackedReal, x::AbstractMatrix, ::Val{0})
    return track(_logabsdetjac_scale, a, data(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::AbstractMatrix, ::Val{0})
    da = data(a)
    J = fill(size(x, 1) / da, size(x, 2))
    return _logabsdetjac_scale(da, data(x), Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end

# adjoints for 1-dim and 2-dim `Scale` using `AbstractVector`
function _logabsdetjac_scale(a::TrackedVector, x::AbstractVector, ::Val{1})
    return track(_logabsdetjac_scale, a, data(x), Val(1))
end
@grad function _logabsdetjac_scale(a::TrackedVector, x::AbstractVector, ::Val{1})
    # ∂ᵢ (∑ⱼ log|aⱼ|) = ∑ⱼ δᵢⱼ ∂ᵢ log|aⱼ|
    #                 = ∂ᵢ log |aᵢ|
    #                 = (1 / aᵢ) ∂ᵢ aᵢ
    #                 = (1 / aᵢ)
    da = data(a)
    J = inv.(da)
    return _logabsdetjac_scale(da, data(x), Val(1)), Δ -> (J .* Δ, nothing, nothing)
end

function _logabsdetjac_scale(a::TrackedVector, x::AbstractMatrix, ::Val{1})
    return track(_logabsdetjac_scale, a, data(x), Val(1))
end
@grad function _logabsdetjac_scale(a::TrackedVector, x::AbstractMatrix, ::Val{1})
    da = data(a)
    Jᵀ = repeat(inv.(da), 1, size(x, 2))
    return _logabsdetjac_scale(da, data(x), Val(1)), Δ -> (Jᵀ * Δ, nothing, nothing)
end

# TODO: implement analytical gradient for scaling a vector using a matrix
# function _logabsdetjac_scale(a::TrackedMatrix, x::AbstractVector, ::Val{1})
#     track(_logabsdetjac_scale, a, data(x), Val{1})
# end
# @grad function _logabsdetjac_scale(a::TrackedMatrix, x::AbstractVector, ::Val{1})
#     throw
# end


text:
"""Transport layer over Message Queues"""
type AMQPTransportClient
    queue::String
    correlation_id::Int
end

immutable AMQPTransport <: AbstractTransport
    queue::String
    mode::Symbol # :client / :server
    mqconn::Any
    mqchan::Any
    consumer_tag::String
    msgchan::Channel{AMQPClient.Message}
    client::AMQPTransportClient

    function AMQPTransport(queue::String, mode::Symbol; virtualhost="/", host="localhost", port=AMQPClient.AMQP_DEFAULT_PORT, auth_params=AMQPClient.DEFAULT_AUTH_PARAMS)
        Logging.info("opening connection...")
        conn = connection(;virtualhost=virtualhost, host=host, port=port, auth_params=auth_params)
        Logging.info("opening channel...")
        chan1 = channel(conn, AMQPClient.UNUSED_CHANNEL, true)
        Logging.info("opening message channel...")
        msgchan = Channel{AMQPClient.Message}(1) # should use unbuffered channels in Julia v0.6

        rpc_fn = (rcvd_msg) -> begin
            put!(msgchan, rcvd_msg)
        end
        Logging.info("initializing as ", mode)

        if mode === :client
            # create a reply queue for the client
            client_queue = string(Base.Random.uuid1())
            success, message_count, consumer_count = queue_declare(chan1, client_queue; exclusive=true)
            @assert success

            # start a consumer task
            success, consumer_tag = basic_consume(chan1, client_queue, rpc_fn)
            @assert success
            client = AMQPTransportClient(client_queue, 0)

            success, message_count, consumer_count = queue_declare(chan1, queue; durable=true)
            @assert success
            # wait till server queue is set up (if queue is not durable)
            while consumer_count == 0
                Logging.info("waiting for consumers...")
                sleep(5)
                success, message_count, consumer_count = queue_declare(chan1, queue; durable=true)
            end
            new(queue, mode, conn, chan1, consumer_tag, msgchan, client)
        elseif mode === :server
            # create a server queue
            Logging.info("creating server queue...")
            success, message_count, consumer_count = queue_declare(chan1, queue; durable=true)
            Logging.info("created server queue...")
            @assert success

            # start a consumer task
            success, consumer_tag = basic_consume(chan1, queue, rpc_fn)
            @assert success
            client = AMQPTransportClient("", 0)
            new(queue, mode, conn, chan1, consumer_tag, msgchan, client)
        else
            error("Unknown mode $mode")
        end
    end
end

function sendrecv(conn::AMQPTransport, msgstr)
    Logging.debug("sending request: ", msgstr)

    if conn.client.correlation_id == typemax(Int)
        conn.client.correlation_id = 0
    end
    conn.client.correlation_id += 1

    M = AMQPClient.Message(convert(Vector{UInt8}, msgstr), content_type="application/octet-stream", delivery_mode=PERSISTENT, reply_to=conn.client.queue, correlation_id=string(conn.client.correlation_id))
    basic_publish(conn.mqchan, M; exchange=default_exchange_name(), routing_key=conn.queue)
    respstr = ""

    while true
        rcvd_msg = take!(conn.msgchan)
        correlation_id = convert(String, rcvd_msg.properties[:correlation_id])
        basic_ack(conn.mqchan, rcvd_msg.delivery_tag)
        # ignore messages if not matching correlation id (could be a resend of something already processed)
        if correlation_id == string(conn.client.correlation_id)
            respstr = String(rcvd_msg.data)
            break
        end
    end

    Logging.debug("received response: ", respstr)
    respstr
end

function sendresp(conn::AMQPTransport, msgstr)
    Logging.debug("sending response: ", msgstr)
    M = AMQPClient.Message(convert(Vector{UInt8}, msgstr), content_type="application/octet-stream", delivery_mode=PERSISTENT, correlation_id=string(conn.client.correlation_id))
    basic_publish(conn.mqchan, M; exchange=default_exchange_name(), routing_key=conn.client.queue)
    nothing
end

function recvreq(conn::AMQPTransport)
    rcvd_msg = take!(conn.msgchan)
    conn.client.queue = convert(String, rcvd_msg.properties[:reply_to])
    conn.client.correlation_id = parse(Int, convert(String, rcvd_msg.properties[:correlation_id]))
    reqstr = String(rcvd_msg.data)
    Logging.debug("received request: ", reqstr)
    basic_ack(conn.mqchan, rcvd_msg.delivery_tag)
    reqstr
end

function close(conn::AMQPTransport)
    if conn.mode === :client
        queue = conn.client.queue
    elseif conn.mode === :server
        queue = conn.queue
    else
        error("Unknown mode $mode")
    end

    queue_purge(conn.mqchan, queue)
    basic_cancel(conn.mqchan, conn.consumer_tag)
    queue_delete(conn.mqchan, queue)
    close(conn.mqchan)
    close(conn.mqconn)
    nothing
end

max_stars_repo_name: tanmaykm/JuliaWebAPIPlugins

max_issues_repo_name: tanmaykm/JuliaWebAPIPlugins

max_forks_repo_name: tanmaykm/JuliaWebAPIPlugins


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown]
# https://discourse.julialang.org/t/how-to-correctly-define-and-use-global-variables-in-the-module-in-julia/65720/8

# %%
# The `::Float64`'s below are redundant.
module A
global a = 100.0::Float64
function f(x)::Float64
    a = x^2
    return a
end
println("f(6.0) = ", f(6.0), ",  a = ", a)
end;

# %%
# A is equivalent to the following.
module B
a = 100.0
function f(x)::Float64
    a = x^2
    return a
end
println("f(6.0) = ", f(6.0), ",  a = ", a)
end;

# %%
# Don't write return-type in general.
# https://github.com/JuliaLang/julia/blob/master/doc/src/manual/functions.md#return-type
module B′
a = 100.0
function f(x)
    a = x^2
    return a
end
println("f(6.0) = ", f(6.0), ",  a = ", a)
end;

# %%
# We need `global` for the immutable global variable `a` in for-loop.
module C
a = 100.0
function f(x)
    global a = x^2
    return a
end
println("f(6.0) = ", f(6.0), ",  a = ", a)
end;

# %%
# We don't need `global` for the content `a[]` of the mutable global variable `a`.
module D
a = Ref(100.0)
function f(x)
    a[] = x^2
    return a[]
end
println("f(6.0) = ", f(6.0), ",  a[] = ", a[])
end;

# %%
# D is type-unstable.  We need `const` for type stability.
module E
const a = Ref(100.0)
function f(x)
    a[] = x^2
    return a[]
end
println("f(6.0) = ", f(6.0), ",  a[] = ", a[])
end;

# %%
@code_warntype D.f(6.0)

# %%
@code_warntype E.f(6.0)

# %%

max_stars_repo_path: 0014/global and const.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0014/global and const.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0014/global and const.jl

max_forks_repo_name: genkuroki/public


text:
using WarpedRegressors: BoxCox

@testset "extras" begin
    rng = MersenneTwister(123456)

    b = BoxCox(0.1)
    binv = inv(b)

    x = randn(rng)
    y = b(x)

    @test binv(b(x)) ≈ x
    @test b(binv(x)) ≈ x

    @test log(abs(first(Zygote.gradient(b, x)))) ≈ Bijectors.logabsdetjac(b, x)
end


text:
using LinearAlgebra
using FFTW

const m = 5
const sigma = 2.0
const K = 200000

@testset "NFFT in multiple dimensions" begin
    for (u,N) in enumerate([(256,), (32,32), (12,12,12), (6,6,6,6)])
      for pre in [NFFT.LUT, NFFT.FULL]
        eps = [1e-7, 1e-3, 1e-6, 1e-4]
        for (l,window) in enumerate([:kaiser_bessel, :gauss, :kaiser_bessel_rev, :spline])
            D = length(N)
            @info "Testing in $D dimensions using $window window"

            M = prod(N)
            x = rand(Float64,D,M) .- 0.5
            p = NFFTPlan(x, N, m, sigma, window, K, precompute = pre,
                         flags = FFTW.ESTIMATE)

            fHat = rand(Float64,M) + rand(Float64,M)*im
            f = ndft_adjoint(p, fHat)
            fApprox = nfft_adjoint(p, fHat)
            e = norm(f[:] - fApprox[:]) / norm(f[:])
            @debug "error adjoint nfft "  e
            @test e < eps[l]

            gHat = ndft(p, f)
            gHatApprox = nfft(p, f)
            e = norm(gHat[:] - gHatApprox[:]) / norm(gHat[:])
            @debug "error nfft "  e
            @test e < eps[l]
        end
      end
    end
end

@testset "High level NFFT " begin
    N = (32,32)
    eps =  1e-3
    D = length(N)

    M = prod(N)
    x = rand(Float64,D,M) .- 0.5

    fHat = rand(Float64,M) + rand(Float64,M)*im
    f = ndft_adjoint(x, N, fHat)
    fApprox = nfft_adjoint(x, N, fHat)
    e = norm(f[:] - fApprox[:]) / norm(f[:])
    @debug "error adjoint nfft "  e
    @test e < eps

    gHat = ndft(x, f)
    gHatApprox = nfft(x, f)
    e = norm(gHat[:] - gHatApprox[:]) / norm(gHat[:])
    @debug "error nfft "  e
    @test e < eps

    f = ndft_adjoint(x, N, ComplexF32.(fHat))
    fApprox = nfft_adjoint(x, N, ComplexF32.(fHat))
    e = norm(f[:] - fApprox[:]) / norm(f[:])
    @debug "error adjoint nfft "  e
    @test e < eps

    gHat = ndft(x, ComplexF32.(f))
    gHatApprox = nfft(x, ComplexF32.(f))
    e = norm(gHat[:] - gHatApprox[:]) / norm(gHat[:])
    @debug "error nfft "  e
    @test e < eps
end

@testset "Abstract sampling points" begin
    M, N = rand(100:200, 2)
    x = range(-0.4, stop=0.4, length=M)
    p = NFFTPlan(x, N, flags = FFTW.MEASURE)
end

@testset "Directional NFFT $D dim" for D in 2:3 begin
    # NFFT along a specified dimension should give the same result as
    # running a 1D NFFT on every slice along that dimension
        eps = 1e-4
        N = tuple( 2*rand(4:8,D)... )
        M = prod(N)
        for d in 1:D
            x = rand(M) .- 0.5

            f = rand(ComplexF64,N)
            p_dir = NFFTPlan(x, d, N)
            fHat_dir = nfft(p_dir, f)
            g_dir = nfft_adjoint(p_dir, fHat_dir)

            p = NFFTPlan(x, N[d])
            fHat = similar(fHat_dir)
            g = similar(g_dir)

            sz = size(fHat)
            Rpre = CartesianIndices( sz[1:d-1] )
            Rpost = CartesianIndices( sz[d+1:end] )
            for Ipost in Rpost, Ipre in Rpre
                idx = [Ipre, :, Ipost]
                fview = f[idx...]
                fHat[idx...] = nfft(p, vec(fview))

                fHat_view = fHat_dir[idx...]
                g[idx...] = nfft_adjoint(p, vec(fHat_view))
            end

            e = norm( fHat_dir[:] - fHat[:] )
            @test e < eps

            e = norm( g_dir[:] - g[:] ) / norm(g[:])
            @test e < eps
        end
    end
end


text:
using Printf
using JLD2
using Oceananigans.Utils
using Oceananigans.Models
using Oceananigans.Utils: TimeInterval, pretty_filesize, prettykeys

using Oceananigans.Fields: boundary_conditions

default_included_properties(::NonhydrostaticModel) = [:grid, :coriolis, :buoyancy, :closure]
default_included_properties(::ShallowWaterModel) = [:grid, :coriolis, :closure]
default_included_properties(::HydrostaticFreeSurfaceModel) = [:grid, :coriolis, :buoyancy, :closure]

"""
    JLD2OutputWriter{I, T, O, IF, IN, KW} <: AbstractOutputWriter

An output writer for writing to JLD2 files.
"""
mutable struct JLD2OutputWriter{O, T, FS, D, IF, IN, KW} <: AbstractOutputWriter
              filepath :: String
               outputs :: O
              schedule :: T
          field_slicer :: FS
            array_type :: D
                  init :: IF
             including :: IN
                  part :: Int
          max_filesize :: Float64
                 force :: Bool
               verbose :: Bool
               jld2_kw :: KW
end

noinit(args...) = nothing

"""
    JLD2OutputWriter(model, outputs; prefix, schedule,
                              dir = ".",
                     field_slicer = FieldSlicer(),
                       array_type = Array{Float32},
                     max_filesize = Inf,
                            force = false,
                             init = noinit,
                        including = [:grid, :coriolis, :buoyancy, :closure],
                          verbose = false,
                             part = 1,
                          jld2_kw = Dict{Symbol, Any}())

Construct a `JLD2OutputWriter` for an Oceananigans `model` that writes `label, output` pairs
in `outputs` to a JLD2 file.

The argument `outputs` may be a `Dict` or `NamedTuple`. The keys of `outputs` are symbols or
strings that "name" output data. The values of `outputs` are either `AbstractField`s, objects that
are called with the signature `output(model)`, or `WindowedTimeAverage`s of `AbstractFields`s,
functions, or callable objects.

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

  ## Filenaming

  - `prefix` (required): Descriptive filename prefixed to all output files.

  - `dir`: Directory to save output to.
           Default: "." (current working directory).

  ## Output frequency and time-averaging

  - `schedule` (required): `AbstractSchedule` that determines when output is saved.

  ## Slicing and type conversion prior to output

  - `field_slicer`: An object for slicing field output in ``(x, y, z)``, including omitting halos.
                    Has no effect on output that is not a field. `field_slicer = nothing` means
                    no slicing occurs, so that all field data, including halo regions, is saved.
                    Default: `FieldSlicer()`, which slices halo regions.

  - `array_type`: The array type to which output arrays are converted to prior to saving.
                  Default: `Array{Float32}`.

  ## File management

  - `max_filesize`: The writer will stop writing to the output file once the file size exceeds `max_filesize`,
                    and write to a new one with a consistent naming scheme ending in `part1`, `part2`, etc.
                    Defaults to `Inf`.

  - `force`: Remove existing files if their filenames conflict.
             Default: `false`.

  ## Output file metadata management

  - `init`: A function of the form `init(file, model)` that runs when a JLD2 output file is initialized.
            Default: `noinit(args...) = nothing`.

  - `including`: List of model properties to save with every file.
                 Default: `[:grid, :coriolis, :buoyancy, :closure]`

  ## Miscellaneous keywords

  - `verbose`: Log what the output writer is doing with statistics on compute/write times and file sizes.
               Default: `false`.

  - `part`: The starting part number used if `max_filesize` is finite.
            Default: 1.

  - `jld2_kw`: Dict of kwargs to be passed to `jldopen` when data is written.

Example
=======

Write out 3D fields for u, v, w, and a tracer c, along with a horizontal average:

```jldoctest jld2_output_writer
using Oceananigans
using Oceananigans.Utils: hour, minute

model = NonhydrostaticModel(grid=RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)), tracers=(:c,))
simulation = Simulation(model, Δt=12, stop_time=1hour)

function init_save_some_metadata!(file, model)
    file["author"] = "Chim Riggles"
    file["parameters/coriolis_parameter"] = 1e-4
    file["parameters/density"] = 1027
    return nothing
end

c_avg =  Field(Average(model.tracers.c, dims=(1, 2)))

# Note that model.velocities is NamedTuple
simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities,
                                                          prefix = "some_data",
                                                          schedule = TimeInterval(20minute),
                                                          init = init_save_some_metadata!)

# output
JLD2OutputWriter scheduled on TimeInterval(20 minutes):
├── filepath: ./some_data.jld2
├── 3 outputs: (u, v, w)
├── field slicer: FieldSlicer(:, :, :, with_halos=false)
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB
```

and a time- and horizontal-average of tracer `c` every 20 minutes of simulation time
to a file called `some_averaged_data.jld2`

```jldoctest jld2_output_writer
simulation.output_writers[:avg_c] = JLD2OutputWriter(model, (; c=c_avg),
                                                     prefix = "some_averaged_data",
                                                     schedule = AveragedTimeInterval(20minute, window=5minute))

# output
JLD2OutputWriter scheduled on TimeInterval(20 minutes):
├── filepath: ./some_averaged_data.jld2
├── 1 outputs: c averaged on AveragedTimeInterval(window=5 minutes, stride=1, interval=20 minutes)
├── field slicer: FieldSlicer(:, :, :, with_halos=false)
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB
```
"""
function JLD2OutputWriter(model, outputs; prefix, schedule,
                                   dir = ".",
                          field_slicer = FieldSlicer(),
                            array_type = Array{Float32},
                          max_filesize = Inf,
                                 force = false,
                                  init = noinit,
                             including = default_included_properties(model),
                               verbose = false,
                                  part = 1,
                               jld2_kw = Dict{Symbol, Any}())

    # Convert each output to WindowedTimeAverage if schedule::AveragedTimeWindow is specified
    schedule, outputs = time_average_outputs(schedule, outputs, model, field_slicer)

    mkpath(dir)
    filepath = joinpath(dir, prefix * ".jld2")
    force && isfile(filepath) && rm(filepath, force=true)

    initialize_jld2_file!(filepath, init, jld2_kw, including, outputs, model)

    return JLD2OutputWriter(filepath, outputs, schedule, field_slicer,
                            array_type, init, including, part, max_filesize,
                            force, verbose, jld2_kw)
end

function initialize_jld2_file!(filepath, init, jld2_kw, including, outputs, model)
    try
        jldopen(filepath, "a+"; jld2_kw...) do file
            init(file, model)
            saveproperties!(file, model, including)

            # Serialize properties in `including`.
            for property in including
                serializeproperty!(file, "serialized/$property", getproperty(model, property))
            end

            # Serialize the location and boundary conditions of each output.
            for (i, (field_name, field)) in enumerate(pairs(outputs))
                file["timeseries/$field_name/serialized/location"] = location(field)
                serializeproperty!(file, "timeseries/$field_name/serialized/boundary_conditions", boundary_conditions(field))
            end
        end
    catch err
        @warn """Initialization of $filepath failed because
                 $(typeof(err)): $(sprint(showerror, err))"""
    end

    return nothing
end

initialize_jld2_file!(writer::JLD2OutputWriter, model) =
    initialize_jld2_file!(writer.filepath, writer.init, writer.jld2_kw, writer.including, writer.outputs, model)

function iteration_zero_exists(filepath)
    file = jldopen(filepath, "r")

    zero_exists = try
        t₀ = file["timeseries/t/0"]
        true
    catch
        false
    finally
        close(file)
    end

    return zero_exists
end

function write_output!(writer::JLD2OutputWriter, model)

    # Catch an error that occurs when a simulation is initialized but not time-stepped:
    if model.clock.iteration == 0 && iteration_zero_exists(writer.filepath)
        if writer.force
            # Re-initialize file:
            rm(writer.filepath, force=true)
            initialize_jld2_file!(writer, model)
        else
            error("Attempting to overwrite data at iteration 0, possibly because a simulation is being
                  re-initialized. Use `force=true` when constructing JLD2OutputWriter to
                  replace any existing files and avoid this error.")
        end
    end

    verbose = writer.verbose

    # Fetch JLD2 output and store in dictionary `data`
    verbose && @info @sprintf("Fetching JLD2 output %s...", keys(writer.outputs))

    tc = Base.@elapsed data = Dict((name, fetch_and_convert_output(output, model, writer)) for (name, output)
                                   in zip(keys(writer.outputs), values(writer.outputs)))

    verbose && @info "Fetching time: $(prettytime(tc))"

    # Start a new file if the filesize exceeds max_filesize
    filesize(writer.filepath) >= writer.max_filesize && start_next_file(model, writer)

    # Write output from `data`
    path = writer.filepath
    verbose && @info "Writing JLD2 output $(keys(writer.outputs)) to $path..."

    start_time, old_filesize = time_ns(), filesize(path)
    jld2output!(path, model.clock.iteration, model.clock.time, data, writer.jld2_kw)
    end_time, new_filesize = time_ns(), filesize(path)

    verbose && @info @sprintf("Writing done: time=%s, size=%s, Δsize=%s",
                              prettytime((end_time - start_time) / 1e9),
                              pretty_filesize(new_filesize),
                              pretty_filesize(new_filesize - old_filesize))

    return nothing
end

"""
    jld2output!(path, iter, time, data, kwargs)

Write the (name, value) pairs in `data`, including the simulation
`time`, to the JLD2 file at `path` in the `timeseries` group,
stamping them with `iter` and using `kwargs` when opening
the JLD2 file.
"""
function jld2output!(path, iter, time, data, kwargs)
    jldopen(path, "r+"; kwargs...) do file
        file["timeseries/t/$iter"] = time
        for (name, datum) in data
            file["timeseries/$name/$iter"] = datum
        end
    end
    return nothing
end

function start_next_file(model, writer::JLD2OutputWriter)
    verbose = writer.verbose
    sz = filesize(writer.filepath)
    verbose && @info begin
        "Filesize $(pretty_filesize(sz)) has exceeded maximum file size $(pretty_filesize(writer.max_filesize))."
    end

    if writer.part == 1
        part1_path = replace(writer.filepath, r".jld2$" => "_part1.jld2")
        verbose && @info "Renaming first part: $(writer.filepath) -> $part1_path"
        mv(writer.filepath, part1_path, force=writer.force)
        writer.filepath = part1_path
    end

    writer.part += 1
    writer.filepath = replace(writer.filepath, r"part\d+.jld2$" => "part" * string(writer.part) * ".jld2")
    writer.force && isfile(writer.filepath) && rm(writer.filepath, force=true)
    verbose && @info "Now writing to: $(writer.filepath)"

    initialize_jld2_file!(writer, model)

    return nothing
end

Base.summary(ow::JLD2OutputWriter) =
    string("JLD2OutputWriter writing ", prettykeys(ow.outputs), " to ", ow.filepath, " on ", summary(ow.schedule))

function Base.show(io::IO, ow::JLD2OutputWriter)

    averaging_schedule = output_averaging_schedule(ow)
    Noutputs = length(ow.outputs)

    print(io, "JLD2OutputWriter scheduled on ", summary(ow.schedule), ":", '\n',
              "├── filepath: ", ow.filepath, '\n',
              "├── $Noutputs outputs: ", prettykeys(ow.outputs), show_averaging_schedule(averaging_schedule), '\n',
              "├── field slicer: ", summary(ow.field_slicer), '\n',
              "├── array type: ", show_array_type(ow.array_type), '\n',
              "├── including: ", ow.including, '\n',
              "└── max filesize: ", pretty_filesize(ow.max_filesize))
end


text:
module Referenceables

# Use markdown files as the docstring:
for (name, path) in [
    :Referenceables => joinpath(dirname(@__DIR__), "README.md"),
    :referenceable => joinpath(@__DIR__, "referenceable.md"),
]
    include_dependency(path)
    str = replace(read(path, String), "```julia" => "```jldoctest $name")
    @eval @doc $str $name
end

export referenceable

struct ReferenceableArray{T, N, A <: AbstractArray{T, N}} <: AbstractArray{Ref{T}, N}
    x::A
end

parenttype(::Type{<:ReferenceableArray{<:Any, <:Any, A}}) where A = A

struct ReferenceableDict{K, V, A <: AbstractDict{K, V}} <: AbstractDict{K, Ref{V}}
    x::A
end

parenttype(::Type{<:ReferenceableDict{<:Any, <:Any, A}}) where A = A

"""
    RefIndexable{inbounds}(x, i::Tuple) where {inbounds isa Bool} <: Ref{valtype(x)}

Like `Base.RefArray`, but works with arbitrary containers and indices.
"""
struct RefIndexable{inbounds, T, A, I} <: Ref{T}
    x::A
    i::I
end

RefIndexable{inbounds}(x, i) where {inbounds} =
    RefIndexable{inbounds, eltype(x), typeof(x), typeof(i)}(x, i)

RefIndexable(x, i) = RefIndexable{false}(x, i)

Base.show(io::IO, x::RefIndexable) =
    if get(io, :limit, false)
        print(io, '↪', x[])
    else
        Base.show_default(io, x)
    end

@inline Base.getindex(x::RefIndexable{inbounds}) where {inbounds} =
    if inbounds
        @inbounds x.x[x.i...]
    else
        x.x[x.i...]
    end

@inline Base.setindex!(x::RefIndexable{inbounds}, value) where {inbounds} =
    if inbounds
        @inbounds x.x[x.i...] = value
    else
        x.x[x.i...] = value
    end

const Referenceable = Union{ReferenceableArray, ReferenceableDict}

@inline function Base.getindex(x::ReferenceableArray, i::Int...)
    @boundscheck checkbounds(x, i...)
    return RefIndexable{true}(x.x, i)
end

@inline function Base.getindex(x::ReferenceableDict, i...)
    return RefIndexable(x.x, i...)
end

referenceable(x::AbstractArray) = ReferenceableArray(x)
referenceable(x::AbstractDict) = ReferenceableDict(x)

Base.size(A::ReferenceableArray) = size(A.x)
Base.axes(A::ReferenceableArray) = axes(A.x)
Base.IndexStyle(::Type{A}) where {A <: ReferenceableArray} =
    Base.IndexStyle(parenttype(A))

Base.length(A::ReferenceableDict) = length(A.x)
Base.iterate(A::ReferenceableDict) = iterate(A.x)
Base.iterate(A::ReferenceableDict, state) = iterate(A.x, state)

end # module


text:
export @experiment_cmd, Experiment

using BSON
using Markdown

Base.@kwdef mutable struct Experiment
    agent::Any
    env::Any
    stop_condition::Any
    hook::Any
    description::String
end

macro experiment_cmd(s)
    Experiment(s)
end

function Experiment(s::String)
    m = match(r"(?<source>\w+)_(?<method>\w+)_(?<env>\w+)(\((?<game>\w*)\))?", s)
    isnothing(m) &&
        throw(ArgumentError("invalid format, got $s, expected format is a local dir or a predefined experiment like dopamine_dqn_atari(pong)`"))
    Experiment(
        Val(Symbol(m[:source])),
        Val(Symbol(m[:method])),
        Val(Symbol(m[:env])),
        m[:game],
    )
end

function Base.run(x::Experiment)
    display(Markdown.parse(x.description))
    run(x.agent, x.env, x.stop_condition, x.hook)
    x
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export grid

"""
    grid(;kwargs...)
    grid(children::Any;kwargs...)
    grid(children_maker::Function;kwargs...)

A Grid component.
Flexbox grid system with variable amount of columns. For more information, see: https://mantine.dev/core/grid/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Content
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `align` (a value equal to: "stretch", "center", "flex-end", "flex-start"; optional): Set grid align-content property
- `className` (String; optional): Often used with CSS to style elements with common properties
- `columns` (Real; optional): Amount of columns in each row
- `grow` (Bool; optional): Should columns in the last row take 100% of grid width
- `gutter` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Spacing between columns predefined value from theme.spacing or number for gutter in px
- `justify` (a value equal to: "space-between", "space-around", "center", "flex-end", "flex-start"; optional): Set grid justify-content property
- `style` (Dict; optional): Inline style override
"""
function grid(; kwargs...)
        available_props = Symbol[:children, :id, :align, :className, :columns, :grow, :gutter, :justify, :style]
        wild_props = Symbol[]
        return Component("grid", "Grid", "dash_mantine_components", available_props, wild_props; kwargs...)
end

grid(children::Any; kwargs...) = grid(;kwargs..., children = children)
grid(children_maker::Function; kwargs...) = grid(children_maker(); kwargs...)


text:
using LinearAlgebra
export TFTDRandomFourier, TFTD, TFTDAAFT, TFTDIAAFT

# Efficient linear regression formula from dmbates julia discourse post (nov 2019)
# https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/27?page=2
function linreg(x, y)
    (N = length(x)) == length(y) || throw(DimensionMismatch())
    ldiv!(cholesky!(Symmetric([float(N) sum(x); 0.0 sum(abs2, x)], :U)), [sum(y), dot(x, y)])
end

function linear_trend(x)
    l = linreg(0.0:1.0:length(x)-1.0 |> collect, x)
    return [l[1] + l[2]*a for a in x]
end

"""
    TFTD(phases::Bool = true, fϵ = 0.05)

The `TFTDRandomFourier` (or just `TFTD` for short) surrogate was proposed by Lucio et al. (2012)[^Lucio2012] as
a combination of truncated Fourier surrogates[^Nakamura2006] ([`TFTS`](@ref)) and
detrend-retrend surrogates.

The `TFTD` part of the name comes from the fact that it uses a combination of truncated
Fourier transforms (TFT) and de-trending and re-trending (D) the time series
before and after surrogate generation. Hence, it can be used to generate surrogates also from
(strongly) nonstationary time series.

## Implementation details

Here, a best-fit linear trend is removed/added from the signal prior to and after generating
the random Fourier signal. In principle, any trend can be removed, but so far, we only provide
the linear option.

See also: [`TFTDAAFT`](@ref), [`TFTDIAAFT`](@ref).

[^Nakamura2006]: Nakamura, Tomomichi, Michael Small, and Yoshito Hirata. "Testing for nonlinearity in irregular fluctuations with long-term trends." Physical Review E 74.2 (2006): 026205.
[^Lucio2012]: Lucio, J. H., Valdés, R., & Rodríguez, L. R. (2012). Improvements to surrogate data methods for nonstationary time series. Physical Review E, 85(5), 056202.
"""
struct TFTDRandomFourier <: Surrogate
    phases::Bool
    fϵ

    function TFTDRandomFourier(phases::Bool, fϵ = 0.05)
        if !(0 < fϵ ≤ 1)
            throw(ArgumentError("`fϵ` must be on the interval  (0, 1] (indicates fraction of lowest frequencies to be preserved)"))
        end
        new(phases, fϵ)
    end
end

const TFTD = TFTDRandomFourier
TFTD() = TFTD(true)
TFTD(fϵ::Real) = TFTDRandomFourier(true, fϵ)

function surrogenerator(x::AbstractVector, rf::TFTD, rng = Random.default_rng())
    # Detrended time series
    m = mean(x)
    trend = linear_trend(x)
    x̂ = x .- m .- trend

    # Pre-plan and allocate Fourier transform
    forward = plan_rfft(x̂)
    inverse = plan_irfft(forward * x̂, length(x̂))
    𝓕 = forward * x̂
    n = length(𝓕)

    # Polar coordinate representation of the Fourier transform
    rx = abs.(𝓕)
    ϕx = angle.(𝓕)
    ϕs = similar(ϕx)

    permutation = zeros(Int, length(x))
    idxs = collect(1:length(x))

    # Initialize surrogate
    s = similar(x)

    init = (forward = forward, inverse = inverse,
        rx = rx, ϕx = ϕx, n = n, m = m,
        𝓕 = 𝓕, ϕs = ϕs,
        trend = trend, x̂ = x̂,
        permutation, idxs)

    return SurrogateGenerator(rf, x, s, init, rng)
end

function (sg::SurrogateGenerator{<:TFTD})()
    fϵ = sg.method.fϵ
    s = sg.s

    init_fields = (:forward, :inverse,
        :rx, :ϕx, :n, :m,
        :𝓕, :ϕs, :trend, :x̂, :permutation, :idxs)

    forward, inverse,
        rx, ϕx, n, m,
        𝓕, ϕs, trend, x̂,
        permutation, idxs = getfield.(Ref(sg.init), init_fields)

    # Surrogate starts out as a random permutation of x̂
    sample!(sg.rng, idxs, permutation; replace = false)
    s .= @view x̂[permutation]

    mul!(𝓕, forward, s)
    ϕs .= angle.(𝓕)

    # Frequencies are ordered from lowest when taking the Fourier
    # transform, so by keeping the 1:n_preserve first phases intact,
    # we are only randomizing the high-frequency components of the
    # signal.
    n_preserve = ceil(Int, abs(fϵ * n))
    ϕs[1:n_preserve] .= @view ϕx[1:n_preserve]

    # Updated spectrum is the old amplitudes with the mixed phases.
    𝓕 .= rx .* exp.(ϕs .* 1im)

    # TODO: Unfortunately, we can't do inverse transform in-place yet, but
    # this is an open PR in FFTW.
    s .= inverse*𝓕 .+ m .+ trend

    return s
end

"""
    TFTDAAFT(fϵ = 0.05)

[`TFTDAAFT`](@ref)[^Lucio2012] are similar to [`TFTD`](@ref) surrogates, but also re-scales
back to the original values of the time series. `fϵ ∈ (0, 1]` is the fraction of the powerspectrum
corresponding to the lowermost frequencies to be preserved.

See also: [`TFTD`](@ref), [`TFTDIAAFT`](@ref).

[^Lucio2012]: Lucio, J. H., Valdés, R., & Rodríguez, L. R. (2012). Improvements to surrogate data methods for nonstationary time series. Physical Review E, 85(5), 056202.
"""
struct TFTDAAFT <: Surrogate
    fϵ

    function TFTDAAFT(fϵ = 0.05)
        if !(0 < fϵ ≤ 1)
            throw(ArgumentError("`fϵ` must be on the interval  (0, 1] (indicates fraction of lowest frequencies to be preserved)"))
        end
        new(fϵ)
    end
end

function surrogenerator(x::AbstractVector, method::TFTDAAFT, rng = Random.default_rng())
    init = (
        gen = surrogenerator(x, TFTS(method.fϵ), rng),
        ix = zeros(Int, length(x)),
        x_sorted = sort(x),
    )
    s = similar(x)

    return SurrogateGenerator(method, x, s, init, rng)
end

function (sg::SurrogateGenerator{<:TFTDAAFT})()
    s = sg.s
    tfts_gen, ix, x_sorted = sg.init.gen, sg.init.ix, sg.init.x_sorted
    s .= tfts_gen()
    sortperm!(ix, s)
    s[ix] .= x_sorted

    return s
end

"""
    TFTDIAAFT(fϵ = 0.05; M::Int = 100, tol::Real = 1e-6, W::Int = 75)

[`TFTDIAAFT`](@ref)[^Lucio2012] are similar to [`TFTDAAFT`](@ref), but adds an iterative
procedure to better match the periodograms of the surrogate and the original time series,
analogously to how [`IAAFT`](@ref) improves upon [`AAFT`](@ref).

`fϵ ∈ (0, 1]` is the fraction of the powerspectrum corresponding to the lowermost
frequencies to be preserved. `M` is the maximum number of iterations. `tol` is the
desired maximum relative tolerance between power spectra. `W` is the number of
bins into which the periodograms are binned when comparing across iterations.

See also: [`TFTD`](@ref), [`TFTDAAFT`](@ref).

[^Lucio2012]: Lucio, J. H., Valdés, R., & Rodríguez, L. R. (2012). Improvements to surrogate data methods for nonstationary time series. Physical Review E, 85(5), 056202.
"""
struct TFTDIAAFT <: Surrogate
    fϵ
    M::Int
    tol::Real
    W::Int

    function TFTDIAAFT(fϵ = 0.05; M::Int = 100, tol::Real = 1e-6, W::Int = 75)
        if !(0 < fϵ ≤ 1)
            throw(ArgumentError("`fϵ` must be on the interval  (0, 1] (indicates fraction of lowest frequencies to be preserved)"))
        end
        new(fϵ, M, tol, W)
    end
end

function surrogenerator(x::AbstractVector, method::TFTDIAAFT, rng = Random.default_rng())
    # Surrogate starts out as a TFTDRandomFourier surrogate
    gen = surrogenerator(x, TFTDRandomFourier(true, method.fϵ), rng)

    # Pre-allocate forward transform for periodogram; can be re-used.
    𝓕, x̂, forward = gen.init.𝓕, gen.init.x̂, gen.init.forward
    𝓕p = prepare_spectrum(x̂, forward)

    # Initial power spectra and their interpolated versions.
    xpower = zeros(length(𝓕));
    powerspectrum!(𝓕p, xpower, x̂, forward)
    spower = copy(xpower)
    xpowerᵦ = interpolated_spectrum(xpower, method.W)
    spowerᵦ = interpolated_spectrum(spower, method.W)

    init = (
        gen = gen,
        ix = zeros(Int, length(x)),
        x_sorted = sort(x),
        𝓕p = 𝓕p,
        xpower = xpower,
        spower = spower,
        xpowerᵦ = xpowerᵦ,
        spowerᵦ = spowerᵦ,
    )

    s = similar(x)
    return SurrogateGenerator(method, x, s, init, rng)
end

function (sg::SurrogateGenerator{<:TFTDIAAFT})()
    x, s, rng = sg.x, sg.s, sg.rng
    fϵ, M, W, tol = sg.method.fϵ, sg.method.M, sg.method.W, sg.method.tol
    n_preserve = ceil(Int, abs(fϵ * length(x)))

    tftd_gen = sg.init.gen
    𝓕, forward, ϕs, ϕx, rx, trend, x̂ = getfield.(Ref(tftd_gen.init),
        (:𝓕, :forward, :ϕx, :ϕs, :rx, :trend, :x̂)
    )
    x_sorted, ix, 𝓕p, xpower, spower, xpowerᵦ, spowerᵦ = getfield.(Ref(sg.init),
        (:x_sorted, :ix, :𝓕p, :xpower, :spower, :xpowerᵦ, :spowerᵦ)
    )

    sum_old, sum_new = 0.0, 0.0
    iter = 1

    # Surrogate starts out as a TFTDRandomFourier realization of `x`.
    s .= tftd_gen()

    while iter <= M
        # Detrend and take transform (steps (vii-viii) in Lucio et al.)
        s .= s .- trend
        mul!(𝓕, forward, s)

        # Rescaling the power spectrum, keeping some percentage of lowermost
        # frequencies, then re-trending (steps vii-x in Lucio et al.)
        ϕs .= angle.(𝓕)
        ϕs[1:n_preserve] .= @view ϕx[1:n_preserve]
        𝓕 .= rx .* exp.(ϕs .* 1im)
        s .= s .+ trend

        # Adjusting amplitudes
        sortperm!(ix, s)
        s[ix] .= x_sorted

        # Compare power spectra
        powerspectrum!(𝓕p, spower, s, forward)
        interpolated_spectrum!(spowerᵦ, spower, W)
        if iter == 1
            sum_old = sum((xpowerᵦ .- xpowerᵦ) .^ 2) / sum(xpowerᵦ .^ 2)
        else
            sum_new = sum((xpowerᵦ .- spowerᵦ) .^ 2) / sum(xpowerᵦ .^ 2)
            if abs(sum_old - sum_new) < tol
                iter = M + 1
            else
                sum_old = sum_new
            end
        end

        iter += 1
    end

    return s
end


text:
using FixedSizeArrays
# I only import it and do not use "using" so that Datastructures.status does not conflict with MathProgBase.status
import DataStructures
import GeometryTypes.decompose, GeometryTypes.isdecomposable

function fulldecompose{T}(poly::Polyhedron{3,T})
  ine = SimpleHRepresentation(getinequalities(poly))
  ext = SimpleVRepresentation(getgenerators(poly))

  # I need to do division so if T is e.g. Integer, I need to use another type
  RT = typeof(one(T)/2)

  A = ine.A
  #rayinface{T<:Real}(r::Vector{T}, i::Integer) = myeqzero(dot(r, A[i,:])) && !myeqzero(r)
  #vertinface{T<:Real}(r::Vector{T}, i::Integer) = myeqzero(dot(r, A[i,:])) && !myeqzero(r)
  R = ext.R
  V = isempty(ext.V) ? [0 0 0] : ext.V

  # Intersection of rays with the limits of the scene
  (xmin, xmax) = extrema(map((j)->V[j,1], 1:size(V,1)))
  (ymin, ymax) = extrema(map((j)->V[j,2], 1:size(V,1)))
  (zmin, zmax) = extrema(map((j)->V[j,3], 1:size(V,1)))
  width = max(xmax-xmin, ymax-ymin, zmax-zmin)
  if width == zero(T)
    width = 2
  end
  scene = HyperRectangle{3,RT}([(xmin+xmax)/2-width, (ymin+ymax)/2-width, (zmin+zmax)/2-width], 2*width*ones(RT,3))
  function exit_point(start, ray)
    times = max((Vector(minimum(scene))-start) ./ ray, (Vector(maximum(scene))-start) ./ ray)
    times[ray .== 0] = Inf # To avoid -Inf with .../(-0)
    time = minimum(times)
    start + time * ray
  end

  triangles = DataStructures.Stack(Tuple{Tuple{Vector{Float64},Vector{Float64},Vector{Float64}},Int64})
  for i in 1:size(A, 1)
    xray = nothing
    yray = nothing
    zray = A[i,:]
    if myeqzero(zray)
      continue
    end
    newface = true
    for j in 1:i-1
      if myeqzero(cross(zray, A[j,:])) && (i in ine.linset || dot(zray, A[j,:]) > 0) # If A[j,:] is almost 0, it is always true...
        # parallel and equality or inequality and same sense
        # TODO is it possible that A[i,:] is stronger than A[j,:] ?
        newface = false
      end
    end
    if !newface
      continue
    end
    if i > 1 && A[i,:] == A[i-1,:] # Same row, only need to check i-1 since the rows are sorted
      continue
    end

    # Checking rays
    counterclockwise(a, b) = dot(cross(a, b), zray)
    line = nothing
    lineleft = false
    lineright = false
    function checkleftright(r::Vector, lin::Bool)
      cc = counterclockwise(r, line)
      if !myeqzero(cc)
        if cc < 0 || lin
          lineleft = true
        end
        if cc > 0 || lin
          lineright = true
        end
      end
    end
    for j in 1:size(R, 1)
      if myeqzero(dot(R[j,:], zray)) && !myeqzero(R[j,:])
        if line != nothing
          checkleftright(R[j,:], j in ext.Rlinset)
        else
          if j in ext.Rlinset
            line = vec(R[j,:])
            if xray != nothing
              checkleftright(xray, false) # false otherwise line wouldn't be nothing
            end
            if yray != nothing
              checkleftright(yray, false)
            end
          end
          if xray == nothing || counterclockwise(R[j,:], xray) > 0
            xray = vec(R[j,:])
          end
          if yray == nothing || counterclockwise(R[j,:], yray)  < 0
            yray = vec(R[j,:])
          end
        end
      end
    end

    # Checking vertices
    face_vert = IntSet([])
    for j in 1:size(V, 1)
      if myeq(dot(V[j,:], zray), ine.b[i])
        push!(face_vert, j)
      end
    end

    if line != nothing
      if isempty(face_vert)
        center = zeros(typeof(V), size(V, 2))
      else
        center = vec(V[first(face_vert), :])
      end
      hull = DataStructures.Stack(Any)
      push!(hull, exit_point(center, line))
      if lineleft
        push!(hull, exit_point(center, cross(zray, line)))
      end
      push!(hull, exit_point(center, -line))
      if lineright
        push!(hull, exit_point(center, cross(line, zray)))
      end
      hulls = (hull,)
    else
      #if length(face_vert) < 3 # Wrong, they are also the rays
      #  error("Not enough vertices and rays to form a face, it may be because of numerical rounding. Otherwise, please report this bug.")
      #end
      if length(face_vert) < 3 && (xray == nothing || (length(face_vert) < 2 && (yray == xray || length(face_vert) < 1)))
        continue
      end
      face_verts = Vector{Int}(length(face_vert))
      idx = 1
      for v in face_vert
        face_verts[idx] = v
        idx += 1
      end
      if xray == nothing
        sweep_norm = cross(zray, [1,0,0])
        if sum(abs(sweep_norm)) == 0
          sweep_norm = cross(zray, [0,1,0])
        end
      else
        sweep_norm = cross(zray, xray)
      end
      sort!(face_verts, by = j -> dot(V[j,:], sweep_norm))
      function getsemihull(sign_sense)
        hull = DataStructures.Stack(Int)
        prev = sign_sense == 1 ? face_verts[1] : face_verts[length(face_verts)]
        cur = prev
        for j in (sign_sense == 1 ? (2:length(face_verts)) : ((length(face_verts)-1):-1:1))
          while prev != cur && counterclockwise(V[cur,:] - V[prev,:], V[face_verts[j],:] - V[prev,:]) >= 0
            cur = prev
            pop!(hull)
            if !isempty(hull)
              prev = DataStructures.top(hull)
            end
          end
          if yray != nothing && counterclockwise(V[face_verts[j],:] - V[cur,:], yray) >= 0
            break
          else
            push!(hull, cur)
            prev = cur
            cur = face_verts[j]
          end
        end
        push!(hull, cur)
        hull
      end
      xtoy_hull = getsemihull(1)
      if yray == nothing
        ytox_hull = getsemihull(-1)
      else
        ytox_hull = DataStructures.Stack(Any)
        push!(ytox_hull, face_verts[1])
        if DataStructures.top(xtoy_hull) != face_verts[1]
          push!(ytox_hull, DataStructures.top(xtoy_hull))
        end
        push!(ytox_hull, exit_point(V[DataStructures.top(xtoy_hull),:], yray))
        push!(ytox_hull, exit_point(V[face_verts[1],:], xray))
      end
      hulls = (xtoy_hull, ytox_hull)
    end
    for hull in hulls
      if length(hull) >= 3
        a = pop!(hull)
        if isa(a, Int)
          a = V[a,:]
        end
        b = pop!(hull)
        if isa(b, Int)
          b = V[b,:]
        end
        while !isempty(hull)
          c = pop!(hull)
          if isa(c, Int)
            c = V[c,:]
          end
          push!(triangles, ((a,b,c), i))
          b = c
        end
      end
    end

  end
  ntri = length(triangles)
  points  = Vector{FixedSizeArrays.Point{3,RT}}(3*ntri)
  faces   = Vector{GeometryTypes.Face{3,Int,0}}(ntri)
  ns = Vector{GeometryTypes.Normal{3,RT}}(3*ntri)
  for i in 1:ntri
    tri = pop!(triangles)
    normal = vec(A[tri[2],:])
    for j = 1:3
      idx = 3*(i-1)+j
      #ns[idx] = -normal
      ns[idx] = normal
    end
    faces[i] = Array(3*(i-1)+(1:3))
    k = 1
    for k = 1:3
      # reverse order of the 3 vertices so that if I compute the
      # normals with the `normals` function, they are in the good
      # sense.
      # I know I don't use the `normals` function but I don't know
      # what is the OpenGL convention so I don't know if it cares
      # about the order of the vertices.
      points[3*i-k+1] = tri[1][k]
    end
  end
  # If the type of ns is Rational, it also works.
  # The normalized array in in float but then it it recast into Rational
  map!(normalize, ns)
  (points, faces, ns)
end

isdecomposable{T<:Point, S<:Polyhedron}(::Type{T}, ::Type{S}) = true
isdecomposable{T<:Face, S<:Polyhedron}(::Type{T}, ::Type{S}) = true
isdecomposable{T<:Normal, S<:Polyhedron}(::Type{T}, ::Type{S}) = true
function decompose{N, T1, T2}(PT::Type{Point{N, T1}}, poly::Polyhedron{N, T2})
  points = fulldecompose(poly)[1]
  map(PT, points)
end
function decompose{N, T, O, T2}(FT::Type{Face{N, T, O}}, poly::Polyhedron{3, T2})
  faces = fulldecompose(poly)[2]
  decompose(FT, faces)
end
function decompose{NT<:Normal, T}(::Type{NT}, poly::Polyhedron{3,T})
  ns = fulldecompose(poly)[3]
  map(NT, ns)
end


text:
export TimeSpec

# Most constants are defined here, except:
#
# BUF           src/Types/SeisIOBuf.jl
# KW            src/Types/KWDefs.jl
# PhaseCat      src/Types/Quake/PhaseCat.jl
# flat_resp     src/Types/InstResp.jl
# RespStage     src/Types/InstResp.jl
# type_codes    src/Types/Methods/0_type_codes.jl
#
# ...and submodule-specific constants, which are in the submodule declaration files (e.g. SEED.jl)

# Type aliases
const ChanSpec    = Union{Integer, UnitRange, Array{Int64, 1}}
const FloatArray  = Union{AbstractArray{Float64, 1}, AbstractArray{Float32, 1}}

@doc """
    TimeSpec = Union{Real, DateTime, String}

# Time Specification
Most functions that allow time specification use two reserved keywords to track
time: `s` (start/begin) time and `t` (termination/end) time. Exact behavior
of each is given in the table below.

* Real numbers are interpreted as seconds
* DateTime values are as in the Dates package
* Strings should use ISO 8601 *without* time zone (`YYYY-MM-DDThh:mm:ss.s`); UTC is assumed. Equivalent Unix `strftime` format codes are `%Y-%m-%dT%H:%M:%S` or `%FT%T`.

## **parsetimewin Behavior**
In all cases, parsetimewin outputs a pair of strings, sorted so that the first string corresponds to the earlier start time.

| typeof(s) | typeof(t) | Behavior                                          |
|:------    |:------    |:-------------------------------------             |
| DateTime  | DateTime  | sort                                              |
| DateTime  | Real      | add *t* seconds to *s*, then sort                 |
| DateTime  | String    | convert *t* => DateTime, then sort                |
| DateTime  | String    | convert *t* => DateTime, then sort                |
| Real      | DateTime  | add *s* seconds to *t*, then sort                 |
| Real      | Real      | treat *s*, *t* as seconds from current time; sort |
| String    | DateTime  | convert *s* => DateTime, then sort                |
| String    | Real      | convert *s* => DateTime, then sort                |

Special behavior with (Real, Real): *s* and *t* are converted to seconds from
the start of the current minute. Thus, for `s=0` (the default), the data request
begins (or ends) at the start of the minute in which the request is submitted.

See also: `Dates.DateTime`
""" TimeSpec
const TimeSpec    = Union{Real, DateTime, String}
const ChanOpts    = Union{String, Array{String, 1}, Array{String, 2}}
const bad_chars = Dict{String, Any}(
  "File" => (0x22, 0x24, 0x2a, 0x2f, 0x3a, 0x3c, 0x3e, 0x3f, 0x40, 0x5c, 0x5e, 0x7c, 0x7e, 0x7f),
  "HTML" => (0x22, 0x26, 0x27, 0x3b, 0x3c, 0x3e, 0xa9, 0x7f),
  "Julia" => (0x24, 0x5c, 0x7f),
  "Markdown" => (0x21, 0x23, 0x28, 0x29, 0x2a, 0x2b, 0x2d, 0x2e, 0x5b, 0x5c, 0x5d, 0x5f, 0x60, 0x7b, 0x7d),
  "SEED" => (0x2e, 0x7f),
  "Strict" => (0x20, 0x21, 0x22, 0x23, 0x24, 0x25, 0x26, 0x27, 0x28, 0x29, 0x2a,
               0x2b, 0x2c, 0x2d, 0x2e, 0x2f, 0x3a, 0x3b, 0x3c, 0x3d, 0x3e, 0x3f,
               0x40, 0x5b, 0x5c, 0x5d, 0x5e, 0x60, 0x7b, 0x7c, 0x7d, 0x7e, 0x7f) )
const datafields = (:id, :name, :loc, :fs, :gain, :resp, :units, :src, :notes, :misc, :t, :x)
const days_per_month = Int32[31,28,31,30,31,30,31,31,30,31,30,31]
const dtchars = (0x2d, 0x2e, 0x30, 0x31, 0x32, 0x33, 0x34, 0x35, 0x36, 0x37, 0x38, 0x39, 0x3a, 0x54)
const dtconst = 62135683200000000
const regex_chars = String[Sys.iswindows() ? "/" : "\\", "\$", "(", ")", "+", "?", "[", "\\0",
"\\A", "\\B", "\\D", "\\E", "\\G", "\\N", "\\P", "\\Q", "\\S", "\\U", "\\U",
"\\W", "\\X", "\\Z", "\\a", "\\b", "\\c", "\\d", "\\e", "\\f", "\\n", "\\n",
"\\p", "\\r", "\\s", "\\t", "\\w", "\\x", "\\x", "\\z", "]", "^", "{", "|", "}"]
const show_os = 8
const sac_float_k = String[ "delta", "depmin", "depmax", "scale", "odelta",
                            "b", "e", "o", "a", "internal1",
                            "t0", "t1", "t2", "t3", "t4",
                            "t5", "t6", "t7", "t8", "t9",
                            "f", "resp0", "resp1", "resp2", "resp3",
                            "resp4", "resp5", "resp6", "resp7", "resp8",
                            "resp9", "stla", "stlo", "stel", "stdp",
                            "evla", "evlo", "evel", "evdp", "mag",
                            "user0", "user1", "user2", "user3", "user4",
                            "user5", "user6", "user7", "user8", "user9",
                            "dist", "az", "baz", "gcarc", "internal2",
                            "internal3", "depmen", "cmpaz", "cmpinc", "xminimum",
                            "xmaximum", "yminimum", "ymaximum", "unused1", "unused2",
                            "unused3", "unused4", "unused5", "unused6", "unused7" ]
const sac_int_k = String[   "nzyear", "nzjday", "nzhour", "nzmin", "nzsec",
                            "nzmsec", "nvhdr", "norid", "nevid", "npts",
                            "internal4", "nwfid", "nxsize", "nysize", "unused8",
                            "iftype", "idep", "iztype", "unused9", "iinst",
                            "istreg", "ievreg", "ievtyp", "iqual", "isynth",
                            "imagtyp", "imagsrc", "unused10", "unused11", "unused12",
                            "unused13", "unused14", "unused15", "unused16", "unused17",
                            "leven", "lpspol", "lovrok", "lcalda", "unused18" ]
const sac_string_k = String["kstnm", "kevnm", "khole", "ko", "ka", "kt0", "kt1", "kt2",
                            "kt3", "kt4", "kt5", "kt6", "kt7", "kt8", "kt9", "kf", "kuser0",
                            "kuser1", "kuser2", "kcmpnm", "knetwk", "kdatrd", "kinst" ]
const sac_nul_c = UInt8[0x2d, 0x31, 0x32, 0x33, 0x34, 0x35, 0x20, 0x20]
const sac_nul_f = -12345.0f0
const sac_nul_i = Int32(-12345)
const sac_nul_start = 0x2d
const sac_nul_Int8 = UInt8[0x31, 0x32, 0x33, 0x34, 0x35]
const segy_ftypes  = Array{DataType, 1}([UInt32, Int32, Int16, Any, Float32, Any, Any, Int8]) # Note: type 1 is IBM Float32
const segy_units = Dict{Int16, String}(0 => "unknown", 1 => "Pa", 2 => "V", 3 => "mV", 4 => "A", 5 => "m", 6 => "m/s", 7 => "m/s2", 8 => "N", 9 => "W")
const seis_inst_codes = ('H', 'J', 'L', 'M', 'N', 'P', 'Z')
const seisio_file_begin = UInt8[0x53, 0x45, 0x49, 0x53, 0x49, 0x4f]
const sμ = 1000000.0
const vSeisIO = Float32(0.54)
const unindexed_fields = (:c, :n)
const webhdr = Dict("User-Agent" => "Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/86.0.4240.111 Safari/537.36") # lol
const xml_endtime = 19880899199000000
const μs = 1.0e-6


text:
function find_hashtags(line)
  matched_phrases = SubString[]
  references = SubString[]
  reg = r"(#\w+)|(#\[[\w|\s]+\])|(@\w+)"i  # eigher a series of characters without white-space or anything inside [] or references identified by @somename
  offset = 1
  m = match(reg, line, offset)
  while m != nothing
    captures = m.captures
    if captures[1] != nothing
      push!(matched_phrases, captures[1][2:end])
    elseif captures[2] != nothing
      push!(matched_phrases, captures[2][3:end-1])
    elseif captures[3] != nothing
      push!(references, captures[3][2:end])
    end
    offset = m.offset + 1
    m = match(reg, line, offset)
  end
  return matched_phrases, references
end

function word_stems(words::Array{SubString})
  td = TokenDocument(words)
  stem!(td)
  return td.tokens
end

"""
Update your knowledge graph.

"""
function update!(maindb::SQLite.DB, inputfile::AbstractString)
  for line in eachline(inputfile)
    line = lstrip(line)
    if startswith(line, "*")  # if this is a proposition
      # Check for hash-tagged words
      nodes, refs = find_hashtags(line)
      if length(nodes) != nothing
        stems = word_stems(nodes)
        # populate tables
        refs = join(refs, ";")
        add_proposition!(maindb, line, nodes, stems, refs)
      end
    end
  end
end

"Copy the bibtex entry from the original file to the programs central file"
function copy_bibentry()
end

## ArgParse
function parse_commandline()
  settings = ap.ArgParseSettings(
                version="0.1",
                prog="Notes2graph",
                description="Create a graph from your notes written in Pandoc format. Each note to be considered should have a bullet point (*). Moreover, keywords should be tagged with a hash (#). If you want to tag compound words that have spaces between them, wrap them in square brackets. E.g. #[compound word].",
                add_version=true,
                add_help=true,
                autofix_names=true,
                )

  argtable = ap.add_arg_table(
                settings,
                "savelocation",
                Dict(
                  :help => "Path to a directory that keeps the database file of Notes2Graph. If the first time running the program, the database file will be saved in this directory.",
                  :required => true,
                  :arg_type => String,
                ),
                ["--notefile", "-n"],
                Dict(
                  :help => "Add your notes to the graph. Supply notes file.",
                  :arg_type => String
                ),
                ["--reffile", "-r"],
                Dict(
                  :help => "If a notefile is mentioned, add a reference file in the form of a bibtex file. This option is not implemented yet.",
                  :arg_type => String
                ),
                ["--descriptions", "-d"],
                Dict(
                  :help => "Shows all the available descriptions for the provided concept. To write concepts with multiple words and spaces, wrap them in double quotations (\")."
                ),
                ["--related", "-t"],
                Dict(
                  :help => "Shows all related concepts. If you want to know more about the relation between each pair of concepts, use -r1 and -r2.",
                  :action => :store_true
                ),
                "--r1",
                Dict(
                  :help => "The first concept of two concepts, to show their relation."
                ),
                "--r2",
                Dict(
                  :help => "The second concept of two concepts, to show their relation."
                )
                )

  parsed_args = ap.parse_args(ARGS, settings, as_symbols=true)
  return parsed_args
end

"Print a list in a readable format"
function pretty_print(yourlist)
  for (index, item) in enumerate(yourlist)
    print("$index  ")
    println(item)
  end
end

function main()
  parsed_args = parse_commandline()

  savelocation = parsed_args[:savelocation]

  # load/create databases
  dbfile = joinpath(savelocation, "notes2graphdb.sqlite")
  bibfile = joinpath(savelocation, "internalbibtex.bib")
  if !isfile(dbfile) # create tables to be filled
    bibfile = create_bibfile(savelocation)
    maindb = initialize_database(savelocation);
  else # load the tables from file
    maindb = load_database(savelocation)
  end

  if parsed_args[:notefile] != nothing
    update!(maindb, parsed_args[:notefile])
  else
    if parsed_args[:descriptions] != nothing
      nodeid = find_nodeid(maindb, parsed_args[:descriptions])
      descrs, refs = descriptions(maindb, nodeid)
      println("# Descriptions for $(parsed_args[:descriptions]):")
      pretty_print(descrs)
      println("# References:")
      pretty_print(refs)
      if parsed_args[:related] == true
        relatedconcepts = related_concepts(maindb, nodeid)
        println("# Related concepts:")
        println(relatedconcepts)
      end
    end
    if parsed_args[:r1] != nothing && parsed_args[:r2] != nothing
      descrs = relation_descr(maindb,  parsed_args[:r1], parsed_args[:r2])
      println("# Relations between $(parsed_args[:r1]) and $(parsed_args[:r2]):")
      pretty_print([i[1] for i in descrs])
      println("# References for the raltions:")
      pretty_print([i[2] for i in descrs])
    end
  end
  # parse text files
  # parse bibtex files and add entries to the databases internal bibtex
  # get stems of words
  # fill tables
  # call information
  # edit entries in tables
  # pretty print results
end