text:
import Base.iterate
using IPFitting.Tools: tfor
using Base.Threads: SpinLock, nthreads

struct ObservationsIterator
   configs::Vector{Dat}   # remember the database
   Icfg::Vector{Int}      # vector of indices over which we iterate
end

"""
```
observations(db::LsqDB, Icfg=:)
observations(configs::Vector{Dat}, Icfg=:)
```

Creates an iterator over all observations stored in the
list of configurations. I.e., the resulting loop will iterate
over all configurations and within each configuration over all
observations stored in that configuration. E.g.,
```
for (okey, cfg, icfg) in observations(configs)
   # okey -> String key
   # cfg :: Dat
   # icfg :: Int, the index of cfg in configs
end
```
"""
observations(db::LsqDB, args...) = observations(db.configs, args...)
observations(configs::Vector{Dat}) = observations(configs, 1:length(configs))
observations(configs::Vector{Dat}, Icfg) =
     ObservationsIterator(configs, collect(Icfg))
observations(configs::Vector{Dat}, ::Colon) =
     observations(configs, 1:length(configs))

iterate(iter::ObservationsIterator) = iterate(iter, 1)

function iterate(iter::ObservationsIterator, i::Integer)
   if i > length(iter.Icfg)
      return nothing
   end
   # the observation keys for the current config
   okeys = collect(keys(iter.configs[iter.Icfg[i]].D))
   # get the next observation
   return iterate(iter, (i=i, okeys=okeys, ikey=1))
end

function iterate(iter::ObservationsIterator, state::NamedTuple)
   # if we got through all observations, then advance to the next config
   if state.ikey > length(state.okeys)
      return iterate(iter, state.i+1)
   end
   # if there is another observation left...
   obskey = state.okeys[state.ikey]
   icfg = iter.Icfg[state.i]
   d = iter.configs[icfg]
   return (obskey, d, icfg),  # returned to user
          (i=state.i, okeys=state.okeys, ikey=state.ikey+1)   # next state
end

tfor_observations(db::LsqDB, callback; kwargs...) =
      tfor_observations(db.configs, callback; kwargs...)

"""
```
tfor_observations(configs::Vector{Dat}, callback; kwargs...)
tfor_observations(db::LsqDB, callback; kwargs...)
```

Create a multi-threaded for loop over all observations in the db of
list of configs. The callback is expected to be of the form
```
callback(n, obskey::AbstractString, cfg::Dat, lock::SpinLock)
```

### Kwargs
* `verbose=true`
* `msg = "Loop over observations"`
"""
function tfor_observations(configs::Vector{Dat}, callback;
                           verbose=true,
                           msg = "Loop over observations",
                           maxnthreads=nthreads())

   # collect a complete list of observations
   idats = Int[]; sizehint!(idats, 3*length(configs))
   okeys = String[]; sizehint!(okeys, 3*length(configs))
   for (okey, _, idat) in observations(configs)
      push!(idats, idat)
      push!(okeys, okey)
   end
   # a rough cost estimate
   costs = [ length(configs[i]) for i in idats ]

   # start a threaded for loop
   lck = SpinLock()
   tfor( n -> callback(n, okeys[n], configs[idats[n]], lck),
         1:length(idats),
         verbose=verbose, msg = msg,
         costs = costs,
         maxnthreads=maxnthreads )

   return nothing
end


text:
function chebyshev_polynomial(order::S,x::Array{T,1}) where {T<:AbstractFloat,S<:Integer}

  polynomial      = Array{T}(undef,length(x),order+1)
  polynomial[:,1] = ones(T,length(x))

  for i = 2:order+1
      for j = 1:length(x)
      if i == 2
        polynomial[j,i] = x[j]
      else
        polynomial[j,i] = 2*x[j]*polynomial[j,i-1]-polynomial[j,i-2]
      end
    end
  end

  return polynomial

end

function chebyshev_polynomial(order::S,x::T) where {T<:AbstractFloat,S<:Integer}

  polynomial = ones(T,1,order+1)

  for i = 2:order+1
    if i == 2
      polynomial[1,i] = x
    else
      polynomial[1,i] = 2*x*polynomial[1,i-1]-polynomial[1,i-2]
    end
  end

  return polynomial

end

ext: jl


text:
module ScenTrees
include("TreeStructure.jl")
include("TreeApprox.jl")
include("StochPaths.jl")
include("LatticeApprox.jl")
include("KernelDensityEstimation.jl")
include("bushinessNesDistance.jl")
export tree_approximation!,lattice_approximation,Tree,Lattice,stage,height,leaves,nodes,root,
        part_tree,build_probabilities!,tree_plot,plot_hd,plot_lattice,gaussian_path1D,gaussian_path2D,
        running_maximum1D,running_maximum2D,path,bushiness_nesdistance,kernel_scenarios
end


text:
using ImageNoise
using ImageCore, ImageTransformations, ImageQualityIndexes
using Test, ReferenceTests, TestImages, Random

include("testutils.jl")

@testset "ImageNoise" begin
# ApplyNoise
@info "Test: ApplyNoise"
include("ApplyNoise/AdditiveWhiteGaussianNoise.jl")

# ReduceNoise
@info "Test: ReduceNoise"
include("ReduceNoise/NonlocalMean.jl")
include("ReduceNoise/BM3DDenoise.jl")
end

nothing


text:
@testset "parse single statement" begin
    free = GTP.Free()
    sf = GTP.SemiFree(1)
    dep = GTP.Dependent(1)

    AB = Segment(:A, :B)
    CD = Segment(:C, :D)

    hp = @hp points(A, B)
    @test isempty(hp.constraints)
    @test hp.points == OrderedDict(:A => free, :B => free)

    hp = @hp Points(A, B, C)
    @test isempty(hp.constraints)
    @test hp.points == OrderedDict(:A => free, :B => free, :C => free)

    hp = @hp P := P ∈ Segment(A, B)
    @test hp.points == OrderedDict(:A => free, :B => free, :P => sf)
    @test hp.constraints == [Appartenence(:P, AB)]

    hp = @hp C := Segment(A, B) ⟂ Segment(C, D)
    @test hp.points == Dict(:A => free, :B => free, :C => sf, :D => free)
    @test hp.constraints == [Perpendicularity(AB, CD)]

    hp = @hp C := Segment(A, B) ∥ Segment(C, D)
    @test hp.points == Dict(:A => free, :B => free, :C => sf, :D => free)
    @test hp.constraints == [Parallelism(AB, CD)]

    hp = @hp C := Segment(A, B) ≅ Segment(C, D)
    @test hp.points == Dict(:A => free, :B => free, :C => sf, :D => free)
    @test hp.constraints == [Congruence(AB, CD)]

    hp = @hp P := Segment(A, B) ∩ Segment(C, D)
    @test hp.points == Dict(:A => free, :B => free, :C => free, :P => dep, :D => free)
    @test hp.constraints == [Intersection(:P, AB, CD)]

    hp = @hp H := A ↓ Segment(C, D)
    @test hp.points == Dict(:A => free, :D => free, :C => free, :H => dep)
    @test hp.constraints == [Projection(:H, :A, CD)]

    hp1 = @hp M := Midpoint(A, B)
    hp2 = @hp M := Midpoint(Segment(A, B))
    hp3 = @hp M := midpoint(A, B)
    hp4 = @hp M := Midpoint(Segment(A, B))
    for hp in (hp1, hp2, hp3, hp4)
        @test hp.points == Dict(:A => free, :B => free, :M => dep)
        @test hp.constraints == [Midpoint(:M, AB)]
    end

    hp1 = @hp D := Parallelogram(A, B, C, D)
    hp2 = @hp C := parallelogram(A, B, C, D)

    @test hp1.points == Dict(:A => free, :B => free, :C => free, :D => dep)
    @test hp2.points == Dict(:A => free, :B => free, :C => dep, :D => free)
    @test hp1.constraints == [Parallelogram(:A, :B, :C, :D)]
    @test hp2.constraints == [Parallelogram(:A, :B, :C, :D)]

    # THESIS
    th = @th Segment(A, B) ∥ Segment(C, D)
    @test th.constraints == [Parallelism(AB, CD)]

    th = @th P ∈ Segment(C, D)
    @test th.constraints == [Appartenence(:P, CD)]
    @test string(th) == "THESIS:\n------------\nP ∈ CD\n"
end

@testset begin "test printing"
    hp = @hp begin
        points(A, B, C)
        C := Segment(A, B) ⟂ Segment(A, C)
        M₁ := Midpoint(A, B)
        M₂ := Midpoint(A, C)
        M₃ := Midpoint(B, C)
        H := A ↓ Segment(B, C)
        O := Circle(M₁, M₂, M₃)
    end

    th = @th H ∈ circle(O, M₁)
    @test string(Theorem(hp, th)) == "POINTS:\n------------\nA : free\nB : free\nC : semifree by (1)\nM₁ : dependent by (2)\nM₂ : dependent by (3)\nM₃ : dependent by (4)\nH : dependent by (5)\nO : dependent by (6)\n\nHYPOTHESIS:\n------------\n(1) AB ⟂ AC\n(2) M₁ ∈ AB ∧ AM₁ ≅ M₁B\n(3) M₂ ∈ AC ∧ AM₂ ≅ M₂C\n(4) M₃ ∈ BC ∧ BM₃ ≅ M₃C\n(5) H ∈ BC ∧ AH ⟂ BC\n(6) OM₁ ≅ OM₂ ≅ OM₃\n\nTHESIS:\n------------\nOH ≅ OM₁\n\nPROOF:\n------\nProof missing\n"

end

@testset "parse multiple statements" begin
    free = GTP.Free()
    dep12 = GTP.Dependent(1, 2)
    dep3 = GTP.Dependent(3)
    hp = @hp begin
        points(A, B, C)
        B := Segment(A, B) ⟂ Segment(B, C)
        B := Segment(A, B) ≅ Segment(B, C)
        M := Midpoint(A, B)
    end

    AB = Segment(:A, :B)
    BC = Segment(:B, :C)

    @test hp.points == Dict(:A => free, :B => dep12, :C => free, :M => dep3)
    @test hp.constraints == [Perpendicularity(AB, BC), Congruence(AB, BC), Midpoint(:M, AB)]

    th = @th begin
        Segment(A, B) ⟂ Segment(B, C)
        Segment(A, B) ≅ Segment(B, C)
    end
    @test th.constraints == [Perpendicularity(AB, BC), Congruence(AB, BC)]
end

@testset "point status algebra" begin
    free = GTP.Free()
    sf1 = GTP.SemiFree(1)
    sf2 = GTP.SemiFree(2)
    dep = GTP.Dependent(3)

    for s in (free, sf1, dep)
        @test GTP._update_status(free, s) == s
        @test GTP._update_status(s, free) == s
    end

    @test GTP._update_status(sf1, sf2) == GTP.Dependent(1, 2)

    @test_throws ArgumentError GTP._update_status(sf1, dep)
    @test_throws ArgumentError GTP._update_status(dep, sf1)
    @test_throws ArgumentError GTP._update_status(dep, dep)
end


text:
# Using user recipes from Plots.

using RecipesBase, Printf, UnPack

# Build a recipe which acts on a custom type.
@recipe function f(meta::MetaVLSV, var::AbstractString; op=:mag, axisunit=RE, normal=:y)
   @unpack ncells, coordmin, coordmax = meta
   if ndims(meta) == 1
      data = readvariable(meta, var)

      x = LinRange(coordmin[1], coordmax[1], ncells[1])

      @series begin
         seriestype --> :line  # use := if you want to force it
         x, data
      end
   elseif ndims(meta) == 2
      # Check if ecliptic or polar run
      if ncells[2] == 1 && ncells[3] != 1
         plotrange = [coordmin[1], coordmax[1], coordmin[3], coordmax[3]]
         sizes = [ncells[1], ncells[2]]
         PLANE = "XZ"
         axislabels = ['X', 'Z']
      elseif ncells[3] == 1 && ncells[2] != 1
         plotrange = [coordmin[1], coordmax[1], coordmin[2], coordmax[2]]
         sizes = [ncells[1], ncells[2]]
         PLANE = "XY"
         axislabels = ['X', 'Y']
      end

      x, y = Vlasiator.get_axis(axisunit, plotrange, sizes)
      data = Vlasiator.prep2d(meta, var, op)'
      unitstr = axisunit == RE ? "R_E" : "m"

      strx = L"\textrm{%$(axislabels[1])}[%$unitstr]"
      stry = L"\textrm{%$(axislabels[2])}[%$unitstr]"

      @series begin
         seriestype --> :heatmap  # use := if you want to force it
         seriescolor --> :turbo
         xguide --> strx
         yguide --> stry
         title --> @sprintf "t= %4.1fs" meta.time
         x, y, data
      end
   elseif ndims(meta) == 3
      pArgs = set_args(meta, var, axisunit; normal, origin)
      data = prep2dslice(meta, var, normal, pArgs)'
      x, y = Vlasiator.get_axis(axisunit, plotrange, sizes)

      @series begin
         seriestype --> :heatmap  # use := if you want to force it
         seriescolor --> :turbo
         title --> @sprintf "t= %4.1fs" meta.time
         x, y, data
      end
   end
end

@userplot VDFSlice

@recipe function f(h::VDFSlice; species="proton", unit=SI, unitv="km/s", slicetype=:nothing,
   vslicethick=0.0, center=:nothing, fmin=-Inf, fmax=Inf, weight=:particle, flimit=-1.0,
   verbose=false)

   if length(h.args) != 2 || !(typeof(h.args[1]) == MetaVLSV) ||
      !(length(h.args[2]) == 3)
      error("Input argument error for VDF slices.  Got: $(typeof.(h.args))")
   end
   meta, location = h.args

   v1, v2, r1, r2, fweight, strx, stry, str_title =
      Vlasiator.prep_vdf(meta, location;
         species, unit, unitv, slicetype, vslicethick, center, weight, flimit, verbose)

   isinf(fmin) && (fmin = minimum(fweight))
   isinf(fmax) && (fmax = maximum(fweight))

   verbose && @info "Active f range is $fmin, $fmax"

   # setups
   legend := false
   grid := true

   #TODO: Plots 1.22.6 needs improvements on colorbar!
   @series begin
      seriestype := :histogram2d
      seriescolor --> :turbo
      colorbar := true
      colorbar_scale --> :log10
      clims --> (fmin, fmax)
      colorbar_title := "f(v)"
      weights --> fweight
      bins --> r1, r2
      v1, v2
   end

end


text:
###### function that reads a fasta file ################
function reading(filepath::String, max_entries=max_num_read_fasta)
    # read the file
    reads = read(filepath, String);
    # process the fasta file to get the DNA part
    dna_reads = [join(split(i, "\n")[2:end]) for i in split(reads, '>') if !isempty(i)];
    dna_reads = length(dna_reads) > max_entries ? dna_reads[1:max_entries] : dna_reads;
    return dna_reads
end

function read_fasta(filepath::String, max_entries=max_num_read_fasta)
    #= read a fasta file =#
    dna_reads = reading(filepath, max_entries);
    # convert all DNA seqeunce to uppercase
    return [uppercase(i) for i in dna_reads]
end

function read_fasta_that_has_ground_truth(filepath::String, max_entries=max_num_read_fasta)
    #= read the fasta file that contains the ground truth
    (ground truths are in capital letters, as fasta files
    from public databases typically contains it), return
    1. dna reads and 2. ground truth

    note: ok maybe it's not really ground truth since
    by defn ground truth is unknown but we still keep that
    labeled info
    =#

    #= read a fasta file =#
    dna_reads = reading(filepath, max_entries);
    # ground truth motifs
    ground_truth_motif_strings = [join([j for j in i if isuppercase(j)]) for i in dna_reads];
    # convert all DNA seqeunce to uppercase
    dna_reads = [uppercase(i) for i in dna_reads];
    return (str = dna_reads,
            labeled_ground_truth = ground_truth_motif_strings)
end

function get_JASPAR_PFM(jaspar_path::String)
    read_jaspar = read(jaspar_path, String);
    s = split(read_jaspar, "\n");
    parsed = [[parse(Float64, i.match) for i in eachmatch(nmbr_regex, s[j])] for j = 2:length(s)-1];
    mat = reduce(hcat, parsed);
    return mat ./ sum(mat, dims=2)
end
########################################################

###### reverse complement of a string ##################
function reverse_complement(s::String)
    join(islowercase(s[si]) ? s[si] : DNA_complement[s[si]] for si = length(s):-1:1)
end
########################################################

###### create data-set for as specified in DNAdataset###
function dna2dummy(dna_string::String; F=Float32)
    v = Array{F,1}(undef, 4*length(dna_string));
    for (index, alphabet) in enumerate(dna_string)
        start = (index-1)*4+1;
        v[start:start+3] = dummy[uppercase(alphabet)];
    end
    return v
end

function sim_data_2_dummy(dna_sim_data_vec::Vector{sim_dna_str_w_motif}; F=Float32)
    how_many_strings = length(dna_sim_data_vec);
    @assert F <: Real "input F must be subtype of Real"
    @assert how_many_strings != 0 "There aren't DNA strings found in the input";
    len = length(dna_sim_data_vec[1].str); # length of each dna string in data
    # note that for the simulated data here, each string always has
    # the same length
    S = Array{F, 2}(undef, (4*len, how_many_strings));
    for i = 1:how_many_strings
        @inbounds S[:, i] = dna2dummy(dna_sim_data_vec[i].str);
    end
    return S
end

"""
get the set of dummy-vectors from a set of dna-strings
the dummy-vectors are all of same length (for now)
"""
function DNA_strings_2_dummy_vectors_same_length(dna_strings::Vector{String}; F=Float32)
    how_many_strings = length(dna_strings);
    @assert how_many_strings != 0 "There aren't DNA strings found in the input";
    len = length(dna_strings[1]); # length of each dna string in data
    S = Array{F, 2}(undef, (4*len, how_many_strings));
    for i = 1:how_many_strings
        @inbounds S[:, i] = dna2dummy(dna_strings[i]);
    end
    return S
end
########################################################


text:
struct Header
    metainfo::Vector{MetaInfo}
    sampleID::Vector{String}

    function Header(metainfo::Vector{MetaInfo}, sampleID::Vector{String})
        return new(metainfo, sampleID)
    end
end

"""
    VCF.Header()
Create an empty VCF header.
"""
function Header()
    return Header(MetaInfo[], String[])
end

"""
    VCF.Header(metainfo::Vector, sampleID::Vector)
Create a VCF header with `metainfo` and `sampleID`.
"""
function Header(metainfo::Vector, sampleID::Vector)
    return Header(convert(Vector{MetaInfo}, metainfo), convert(Vector{String}, sampleID))
end

function Base.eltype(::Type{Header})
    return MetaInfo
end

function Base.length(header::Header)
    return length(header.metainfo)
end

function Base.iterate(header::Header, i = 1)
    return i > lastindex(header.metainfo) ? nothing : (header.metainfo[i], i + 1)
end

function Base.findall(header::Header, tag::AbstractString)
    return Base.filter(m -> isequaltag(m, tag), header.metainfo)
end

function Base.pushfirst!(header::Header, metainfo)
    pushfirst!(header.metainfo, convert(MetaInfo, metainfo))
    return header
end

function Base.push!(header::Header, metainfo)
    push!(header.metainfo, convert(MetaInfo, metainfo))
    return header
end

function Base.show(io::IO, header::Header)
    println(io, summary(header), ':')
    tags = BioCore.metainfotag.(header.metainfo)
    println(io, "  metainfo tags: ", join(unique(tags), ' '))
      print(io, "     sample IDs: ", join(header.sampleID, ' '))
end

function Base.write(io::IO, header::Header)
    n = 0
    for metainfo in header.metainfo
        n += write(io, metainfo, '\n')
    end
    n += write(io, "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO")
    if !isempty(header.sampleID)
        n += write(io, '\t', "FORMAT")
        for id in header.sampleID
            n += write(io, '\t', id)
        end
    end
    n += write(io, '\n')
    return n
end


text:
### press-rybicki.jl ---  Perform fast but approximate Lomb–Scargle periodogram
#
# Copyright (C) 2016 Astropy Developers
# Copyright (C) 2017 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: periodogram, lomb scargle
#
# This file is a part of LombScargle.jl.
#
# License is BSD 3-clause "New" or "Revised".
#
### Code:

struct FastGLSPlan{T,A,B<:AbstractVector{T},C,D,E,F,G,H,I,J,K,L} <: PeriodogramPlan
    times::A
    signal::B
    freq::C
    sumw::D
    w::E
    y::B
    YY::T
    fftgrid::F
    bfft_vect::G
    bfft_grid::H
    bfft_plan::I
    Mfft::J
    noise::K
    norm::Symbol
    P::L
end

struct FastGLSPlan_fit_mean{T,A,B<:AbstractVector{T},C,D,E,F,G,H,I,J,K,L} <: PeriodogramPlan
    times::A
    signal::B
    freq::C
    sumw::D
    w::E
    y::B
    YY::T
    fftgrid::F
    bfft_vect::G
    bfft_grid::H
    bfft_plan::I
    Mfft::J
    noise::K
    norm::Symbol
    P::L
end

include("extirpolation.jl")

# Fast, but approximate, method to compute the Lomb–Scargle periodogram for
# evenly spaced frequency grid.  See
# * Press, W. H., Rybicki, G. B. 1989, ApJ, 338, 277 (URL:
#   http://dx.doi.org/10.1086/167197,
#   Bibcode: http://adsabs.harvard.edu/abs/1989ApJ...338..277P)
# This is adapted from Astropy implementation of the method.  See
# `lombscargle_fast' function.

function _press_rybicki!(P, times::AbstractVector{<:Real}, y::AbstractVector{<:Real},
                         w::AbstractVector{<:Real}, t0, df, N , f0, YY::Real,
                         fftgrid, bfft_vec::AbstractVector{Complex{T}},
                         grid, plan, Nfft::Integer, Mfft::Integer) where {T<:AbstractFloat}
    #---------------------------------------------------------------------------
    # 1. compute functions of the time-shift tau at each frequency
    Ch, Sh = trig_sum!(grid, fftgrid, bfft_vec, plan, times, w .* y, df, N, Nfft, t0, f0, 1, Mfft)
    C2, S2 = trig_sum!(grid, fftgrid, bfft_vec, plan, times, w,      df, N, Nfft, t0, f0, 2, Mfft)
    #---------------------------------------------------------------------------
    # 2. Compute the periodogram, following Zechmeister & Kurster
    #    and using tricks from Press & Rybicki.
    tan_2ωτ = S2 ./ C2
    C2w = 1 ./ (sqrt.(1 .+ tan_2ωτ .* tan_2ωτ)) # = cos(2 * ωτ)
    S2w = tan_2ωτ .* C2w # = sin(2 * ωτ)
    Cw  = sqrt.((1 .+ C2w) ./ 2) # = cos(ωτ)
    Sw  = sign.(S2w) .* sqrt.((1 .- C2w) ./ 2) # = sin(ωτ)
    return P .= ((Ch .* Cw .+ Sh .* Sw) .^ 2 ./ ((1 .+ C2 .* C2w .+ S2 .* S2w) ./ 2) .+
                 (Sh .* Cw .- Ch .* Sw) .^ 2 ./ ((1 .- C2 .* C2w .- S2 .* S2w) ./ 2)) ./ YY
end

function _press_rybicki_fit_mean!(P, times::AbstractVector{<:Real},
                                  y::AbstractVector{<:Real}, w::AbstractVector{<:Real},
                                  t0, df, N , f0, YY::Real, fftgrid,
                                  bfft_vec::AbstractVector{Complex{T}}, grid, plan,
                                  Nfft::Integer, Mfft::Integer) where {T<:AbstractFloat}
    #---------------------------------------------------------------------------
    # 1. compute functions of the time-shift tau at each frequency
    Ch, Sh = trig_sum!(grid, fftgrid, bfft_vec, plan, times, w .* y, df, N, Nfft, t0, f0, 1, Mfft)
    C2, S2 = trig_sum!(grid, fftgrid, bfft_vec, plan, times, w,      df, N, Nfft, t0, f0, 2, Mfft)
    C, S   = trig_sum!(grid, fftgrid, bfft_vec, plan, times, w,      df, N, Nfft, t0, f0, 1, Mfft)
    #---------------------------------------------------------------------------
    # 2. Compute the periodogram, following Zechmeister & Kurster
    #    and using tricks from Press & Rybicki.
    tan_2ωτ = (S2 .- 2 .* S .* C) ./ (C2 .- (C .* C .- S .* S))
    C2w = 1 ./ (sqrt.(1 .+ tan_2ωτ .* tan_2ωτ)) # = cos(2 * ωτ)
    S2w = tan_2ωτ .* C2w # = sin(2 * ωτ)
    Cw  = sqrt.((1 .+ C2w) ./ 2) # = cos(ωτ)
    Sw  = sign.(S2w) .* sqrt.((1 .- C2w) ./ 2) # = sin(ωτ)
    return P .= ((Ch .* Cw .+ Sh .* Sw) .^ 2 ./
                 ((1 .+ C2 .* C2w .+ S2 .* S2w) ./ 2 .- (C .* Cw .+ S .* Sw) .^ 2) .+
                 (Sh .* Cw .- Ch .* Sw) .^ 2 ./
                 ((1 .- C2 .* C2w .- S2 .* S2w) ./ 2 .- (S .* Cw .- C .* Sw) .^ 2)) ./ YY
end

_periodogram!(P::AbstractVector, times, p::FastGLSPlan) =
    _press_rybicki!(P, times, p.y, p.w, minimum(p.times), step(p.freq),
                    length(p.freq), minimum(p.freq), p.YY, p.fftgrid, p.bfft_vect,
                    p.bfft_grid, p.bfft_plan, length(p.bfft_vect), p.Mfft)
_periodogram!(P::AbstractVector, times, p::FastGLSPlan_fit_mean) =
    _press_rybicki_fit_mean!(P, times, p.y, p.w, minimum(p.times), step(p.freq),
                             length(p.freq), minimum(p.freq), p.YY, p.fftgrid, p.bfft_vect,
                             p.bfft_grid, p.bfft_plan, length(p.bfft_vect), p.Mfft)


text:
module SuperLU

    import Base: \, lufact, size, nfilled

    export
        DMat,
        NCMat,
        superlu_options_t,
        SuperLUStat_t,
        SuperMatrix,
        nv,
        sp_ienv

    include("superlu_h.jl")
    include("superlu_types.jl")
    include("SuperLUCode.jl")

end # module


text:
"""
    Poisson Likelihood(λ=1.0)

## Arguments
- `λ::Real` : Poisson rate

---

[Poisson Likelihood](https://en.wikipedia.org/wiki/Poisson_distribution) where a Poisson distribution is defined at every point in space (careful, it's different from continous Poisson processes)
```math
    p(y|f) = Poisson(y|\\lambda \\sigma(f))
```
Where `σ` is the logistic function.
Augmentation details will be released at some point (open an issue if you want to see them)
"""
mutable struct PoissonLikelihood{T<:Real,A<:AbstractVector{T}} <: EventLikelihood{T}
    λ::T
    c::A
    γ::A
    θ::A
    function PoissonLikelihood{T}(λ::T) where {T<:Real}
        return new{T,Vector{T}}(λ)
    end
    function PoissonLikelihood{T}(
        λ::T, c::A, γ::A, θ::A
    ) where {T<:Real,A<:AbstractVector{T}}
        return new{T,A}(λ, c, γ, θ)
    end
end

function PoissonLikelihood(λ::T=1.0) where {T<:Real}
    return PoissonLikelihood{T}(λ)
end

implemented(::PoissonLikelihood, ::Union{<:AnalyticVI,<:GibbsSampling}) = true

function init_likelihood(
    likelihood::PoissonLikelihood{T}, ::AbstractInference{T}, ::Integer, nSamplesUsed::Int
) where {T}
    return PoissonLikelihood{T}(
        likelihood.λ, zeros(T, nSamplesUsed), zeros(T, nSamplesUsed), zeros(T, nSamplesUsed)
    )
end

function (l::PoissonLikelihood)(y::Real, f::Real)
    return pdf(Poisson(get_p(l, l.λ, f)), y)
end

function Distributions.loglikelihood(l::PoissonLikelihood, y::Real, f::Real)
    return logpdf(Poisson(expec_count(l, f)), y)
end

function expec_count(l::PoissonLikelihood, f)
    return get_p(l, l.λ, f)
end

function get_p(::PoissonLikelihood, λ::Real, f)
    return λ * logistic.(f)
end

function Base.show(io::IO, l::PoissonLikelihood{T}) where {T}
    return print(io, "Poisson Likelihood (λ = $(l.λ))")
end

function compute_proba(
    l::PoissonLikelihood{T}, μ::AbstractVector{<:Real}, σ²::AbstractVector{<:Real}
) where {T<:Real}
    N = length(μ)
    pred = zeros(T, N)
    sig_pred = zeros(T, N)
    for i in 1:N
        x = pred_nodes .* sqrt.(max(σ²[i], zero(T))) .+ μ[i]
        pred[i] = dot(pred_weights, get_p(l, l.λ, x))
        sig_pred[i] = dot(pred_weights, get_p(l, l.λ, x) .^ 2) - pred[i]^2
    end
    return pred, sig_pred
end

### Local Updates ###

function local_updates!(
    l::PoissonLikelihood{T}, y::AbstractVector, μ::AbstractVector, diag_cov::AbstractVector
) where {T}
    @. l.c = sqrt(abs2(μ) + diag_cov)
    @. l.γ = 0.5 * l.λ * safe_expcosh(-0.5 * μ, 0.5 * l.c)
    @. l.θ = (y + l.γ) / l.c * tanh(0.5 * l.c)
    return l.λ = sum(y) / sum(expectation.(logistic, μ, diag_cov))
end

function sample_local!(l::PoissonLikelihood, y::AbstractVector, f::AbstractVector)
    @. l.γ = rand(Poisson(l.λ * logistic(f))) # Sample n
    return set_ω!(l, rand.(PolyaGamma.(y + Int.(l.γ), abs.(f)))) # Sample ω
end

### Global Updates ###

@inline ∇E_μ(l::PoissonLikelihood, ::AOptimizer, y::AbstractVector) = (0.5 * (y - l.γ),)
@inline ∇E_Σ(l::PoissonLikelihood, ::AOptimizer, y::AbstractVector) = (0.5 * l.θ,)

## ELBO Section ##
function expec_loglikelihood(
    l::PoissonLikelihood{T}, ::AnalyticVI, y, μ::AbstractVector, Σ::AbstractVector
) where {T}
    tot = 0.5 * (dot(μ, (y - l.γ)) - dot(l.θ, abs2.(μ)) - dot(l.θ, Σ))
    tot += Zygote.@ignore(
        sum(y * log(l.λ)) - sum(logfactorial, y) - logtwo * sum((y + l.γ))
    )
    return tot
end

AugmentedKL(l::PoissonLikelihood, y::AbstractVector) = PoissonKL(l) + PolyaGammaKL(l, y)

PoissonKL(l::PoissonLikelihood) = PoissonKL(l.γ, l.λ)

PolyaGammaKL(l::PoissonLikelihood, y::AbstractVector) = PolyaGammaKL(y + l.γ, l.c, l.θ)


text:
module Struct
export @struct
import Compat.String

function struct_string(bytes::Vector{UInt8})
    last_byte = findfirst(bytes, 0)
    String(last_byte == 0 ? bytes : bytes[1:last_byte-1])
end

macro struct(typename, contents)
    if contents.head != :block
        error("Invalid struct declaration")
    end

    read_blk = Expr(:call, esc(typename))
    write_blk = Expr(:block)

    for typedecl in contents.args
        if typedecl.head == :line
            continue
        elseif typedecl.head != :(::)
            error("invalid struct declaration")
        end
        fieldname = typedecl.args[1]
        fieldtype = typedecl.args[2]
        field = :(x.$(fieldname))
        if isa(fieldtype, Expr) && fieldtype.head != :curly
            if fieldtype.head != :call
                error("invalid struct declaration")
            end
            # Type has size parameters

            if fieldtype.args[1] in (:ASCIIString, :UTF8String, :String)
                typedecl.args[2] = fieldtype.args[1]
                push!(read_blk.args, :(struct_string(read(io, UInt8, ($(fieldtype.args[2:end]...))))))
                push!(write_blk.args, quote
                    if sizeof($field) > $(fieldtype.args[2])
                        throw(ArgumentError($("size of $(fieldtype.args[1]) must be less than $(fieldtype.args[2])")))
                    end
                    write(io, $field)
                    for i = sizeof($field)+1:$(fieldtype.args[2])
                        write(io, UInt8(0))
                    end
                end)
            else
                typedecl.args[2] = :(Array{$(fieldtype.args[1]), $(length(fieldtype.args)-1)})
                push!(read_blk.args, :(read(io, $(esc(fieldtype.args[1])), ($(fieldtype.args[2:end]...)))))
                push!(write_blk.args, quote
                    if size($field) != $(tuple(fieldtype.args[2:end]...))
                        throw(ArgumentError($("size of $(fieldname) must be $((fieldtype.args[2:end]...,)), but was ")*string(size($field))))
                    end
                    write(io, $field)
                end)
            end
        else
            push!(read_blk.args, :(read(io, $fieldtype)))
            push!(write_blk.args, :(write(io, $field)))
        end
    end
    quote
        type $(esc(typename))
            $(esc(contents))
        end
        Base.read(io::IO, ::Type{$(esc(typename))}) = $read_blk
        Base.write(io::IO, x::$(esc(typename))) = $write_blk
    end
end
end


text:
include("Core.jl")
include("LocalHamiltonian.jl")
include("RegionIter.jl")

# Fused Hamiltonian
include("FusedHamiltonian.jl")


text:
#temp

function very_crude_c_solver(y::Float64,i::Int,j::Int,k::Int,maph::MAPHDist)
    quadgk(u -> (maph.α*exp(maph.T*u))[i]*(exp(maph.T*(y-u))*maph.T0[:,j])[k] , 0, y, rtol=1e-8) |> first
end

function poi_bound(λ, ε)
    K = 0
    pdf = (λ^K/factorial(K))*exp(-λ)

    while pdf < 1-ε
        K += 1
        pdf += (λ^K/factorial(K))*exp(-λ)
    end

    return K

end

function poi_prob(k, λ)

    return (λ^k/factorial(k))*exp(-λ)

end

function uniformizaition_integral(y::Float64, row_vec, col_vec, maph::MAPHDist)
    r = maximum(abs.(diag(maph.T)))
    P = I + maph.T./r
    R = poi_bound(r*y,0.1)

    column_vecs = zeros(length(col_vec),R)

    column_vecs[:,1] = col_vec
    for i = 2:R
        column_vecs[:,i] = P*column_vecs[:,i-1]
    end

    row_vecs = zeros(R,length(row_vec))
    row_vecs[1,:] = (poi_prob(R+1,r*y)*row_vec)

    for i = 2:R
        row_vecs[i,:] = row_vecs[i-1,:]'*P + poi_prob(i+1,r*y)*row_vec
    end

    row_vecs = row_vecs[end:-1:1,1:1:end]
    integral = 0
    for i = 1:R
        integral += row_vecs[i,:]'*column_vecs[:,i]/r
    end

    return integral
end

function sufficient_stats(observation::SingleObs, maph::MAPHDist; c_solver = very_crude_c_solver)::MAPHSufficientStats
    stats = MAPHSufficientStats(maph)

    p, q = model_size(maph)

    a(y::Float64) = maph.α*exp(maph.T*y)
    b(y::Float64,j::Int) = exp(maph.T*y)*maph.T0[:,j]
    c(y::Float64,i::Int,j::Int,k::Int) = very_crude_c_solver(y,i,j,k,maph)

    D = Diagonal(maph.T)
    PT = I-inv(Diagonal(maph.T))*maph.T
    PT0 = -inv(Diagonal(maph.T))*maph.T0
    A = inv(I-PT)*PT0
    PA = maph.α*A

    EB(y::Float64, i::Int, j::Int) = maph.α[i] * b(y,j)[i]*PA[j]/ (maph.α*b(y,j))
    EZ(y::Float64, i::Int, j::Int,k::Int) = c(y,i,j,k)*PA[j]/(maph.α*b(y,j))
    ENT(y::Float64,i::Int,j::Int,k::Int) = i !=k ? maph.T[i,:].*c(y,i,j,k)*PA[j]/(maph.α*b(y,j)) : zeros(p)
    ENA(y::Float64,i::Int,j::Int) = PA[j]*a(y)[i]*maph.T0[i,j]/(maph.α*b(y,j))

    stats.B = [sum([EB(observation.y, i, j) for j = 1:q]) for i =1:p]
    stats.Z = [sum([EZ(observation.y,i,j,i) for j =1:q]) for i = 1:p]

    for i = 1:p
        for k = (q+1):(q+p)
            V = sum([ENT(observation.y,i,j,k-q) for j in 1:q])
            stats.N[i,k] = V[k-q]
        end

        for j = 1:q
            stats.N[i,j] = ENA(observation.y,i,j)
        end
    end

    return stats
end

sufficient_stats(data::MAPHObsData, maph::MAPHDist; c_solver = very_crude_c_solver)::MAPHSufficientStats = mean([sufficient_stats(d, maph) for d in data])


text:
### A Pluto.jl notebook ###
# v0.17.2

using Markdown
using InteractiveUtils

# ╔═╡ 227de887-4284-46fb-83e2-22d33830c9cf
using LaTeXStrings

# ╔═╡ 1eb8d032-dbeb-4544-8bcd-e512f19d6cd5
θ = round(1/3, digits=2)

# ╔═╡ ec80fb74-5aa2-11ec-0e31-bb7b0d8fe4f7
md"""
``\theta`` = $(round (θ, digits = 2))
"""

# ╔═╡ 76f42d5d-7b83-4d11-a52b-261ad8b55567
md"""
``\theta`` = $θ
"""

# ╔═╡ 18fbc708-96ea-4a2a-b631-442295d59b86
md"""
θ = $(round (θ, digits = 2))
"""

# ╔═╡ 147eb4f3-bbf7-49d4-b7db-2c017878ff21
θ

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"

[compat]
LaTeXStrings = "~1.3.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0"
manifest_format = "2.0"

[[deps.LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"
"""

# ╔═╡ Cell order:
# ╠═227de887-4284-46fb-83e2-22d33830c9cf
# ╠═1eb8d032-dbeb-4544-8bcd-e512f19d6cd5
# ╠═ec80fb74-5aa2-11ec-0e31-bb7b0d8fe4f7
# ╠═76f42d5d-7b83-4d11-a52b-261ad8b55567
# ╟─18fbc708-96ea-4a2a-b631-442295d59b86
# ╠═147eb4f3-bbf7-49d4-b7db-2c017878ff21
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
export DGCNN, EdgeConv

function CreateSingleKNNGraph(X::AbstractArray, K::Int)
    F, N = size(X)
    kdtree = KDTree(X)
    cat([X[:, (knn(kdtree, X[:, i], K + 1, true)[1][2:K+1])] for i = 1:N]..., dims = 3)
end

@nograd CreateSingleKNNGraph

struct EdgeConv
    layers::AbstractArray
    K::Int
    mlp::Any
    maxpool_K::Any
end

function EdgeConv(layers::AbstractArray, K::Int)
    m = []

    for i = 1:(length(layers)-1)
        if i == 1
            push!(m, conv_bn_block(2 * layers[i], layers[i+1], (1,)))
        else
            push!(m, conv_bn_block(layers[i], layers[i+1], (1,)))
        end
    end

    EdgeConv(layers, K, Chain(m...), MaxPool((K,)))
end

function (m::EdgeConv)(X)
    F, N, B = size(X)
    # X: [F, N, B]

    KNNGraph = cat([CreateSingleKNNGraph(X[:, :, i], m.K) for i = 1:B]..., dims = 4)
    # KNNGraph: [F, K, N, B]

    X = reshape(X, F, 1, N, B)
    # X: [F, 1, N, B]

    X = cat([X for i = 1:m.K]..., dims = 2)
    # X: [F, K, N, B]

    X = cat(X, KNNGraph - X, dims = 1)
    # X: [2*F, K, N, B]

    X = PermutedDimsArray(X, (2, 3, 1, 4))
    # X: [K, N, 2*F, B]

    X = reshape(X, N * m.K, 2 * F, B)
    # X: [K*N, 2*F, B]

    X = m.mlp(X)
    # X: [K*N, an, B]

    an = size(X)[2]
    X = reshape(X, m.K, an * N, B)
    # X: [K, N*an, B]

    X = m.maxpool_K(X)
    # X: [1, N*an, B]

    X = reshape(X, N, an, B)
    # X: [N, an, B]

    X = PermutedDimsArray(X, (2, 1, 3))
    # X: [an, N, B]

    return X
end

@functor EdgeConv

"""
    DGCNN(num_classes::Int=10, K::Int=10, npoints::Int=1024)

Flux implementation of Dynamic-Graph CNN classification model.

### Fields:

- `num_classes` - Number of classes in dataset.
- `K`           - k nearest neighbour to be used EdgeConv.
- `npoints`     - Number of points in input PointCloud.

"""
struct DGCNN
    EdgeConv1::EdgeConv
    EdgeConv2::EdgeConv
    conv_3::Any
    maxpool_3::Any
    fc_4::Any
    drop_4::Any
    fc_5::Any
    drop_5::Any
    fc_6::Any
end

function DGCNN(num_classes::Int = 10, K::Int = 10, npoints::Int = 1024)
    DGCNN(
        EdgeConv([3, 32, 64, 64], K),
        EdgeConv([64, 128, 256], K),
        conv_bn_block(256, 1024, (1,)),
        MaxPool((npoints,)),
        fc_bn_block(1024, 512),
        Dropout(0.5),
        fc_bn_block(512, 256),
        Dropout(0.5),
        Dense(256, num_classes),
    )
end

function (m::DGCNN)(X)

    # X: [3, N, B]

    C, N, B = size(X)
    X = m.EdgeConv1(X)
    # X: [64, N, B]

    X = m.EdgeConv2(X)
    # X: [128, N, B]

    X = PermutedDimsArray(X, (2, 1, 3))
    # X: [N, 128, B]

    X = m.conv_3(X)
    # X: [N, 1024, B]

    X = m.maxpool_3(X)
    # X: [1, 1024, B]

    X = reshape(X, :, B)
    # X: [1024, B]

    X = m.drop_4(m.fc_4(X))
    # X: [512, B]

    X = m.drop_5(m.fc_5(X))
    # X: [256, B]

    X = m.fc_6(X)
    # X: [num_classes, B]

    return softmax(X, dims = 1)

end

@functor DGCNN


text:
"""
    JarrowRudd()

[Jarrow-Rudd model](http://simulations.lpsm.paris/binomial_trees/jr).
"""
struct JarrowRudd <: Model end

"""
    evaluate(O, JarrowRudd(), risk_neutral = true, N = 1000)

Evaluate option `O` using `JarrowRudd` binomial model (defaults to the risk-neutral version).

# Arguments
- `O::Option`: option
- `risk_neutral`: `true` if risk neutral, `false` if equal probability.
- `N`: number of paths to simulate
"""
function evaluate(O::Option, m::JarrowRudd, risk_neutral::Bool = true, N::Int64 = 1000)
    Δt = O.t / N
    R = exp(O.r * Δt)

    if risk_neutral
         U = exp((O.r - O.σ^2 / 2) * Δt + O.σ * √Δt)
         D = exp((O.r - O.σ^2 / 2) * Δt - O.σ * √Δt)
         p = (R - D) / (U - D)
         q = (U - R) / (U - D)
    else
        p = q = 0.5
    end

    if O.call == -1
        Z = [max(0, O.k - O.s * exp((2 * i - N) * O.σ * √Δt)) for i = 0:N]
    elseif O.call == 1
        Z = [max(0, O.s * exp((2 * i - N) * O.σ * √Δt) - O.k) for i = 0:N]
    end

    for n = N-1:-1:0, i = 0:n
        if O.call == -1
            x = O.k - O.s * exp((2 * i - n) * O.σ * √Δt)
        elseif O.call == 1
            x = O.s * exp((2 * i - n) * O.σ * √Δt) - O.k
        end
        y = (q * Z[i+1] + p * Z[i+2]) / exp(O.r * Δt)
        Z[i+1] = max(x, y)
    end

    return Z[1]
end


text:
using OrdinaryDiffEq, StochasticDelayDiffEq, ParameterizedFunctions, Test, Random
using ParameterizedFunctions.ModelingToolkit # macro hygiene
f = @ode_def LotkaVolterra begin
  dx = 1.5x - x*y
  dy = -3y + x*y
end

prob = ODEProblem(f,big.([1.0;1.0]),(big(0.0),big(10.0)))

using DiffEqDevTools

dts = big(1/2).^(8:-1:4)
test_setup = Dict(:alg=>Vern9(),:reltol=>1e-25,:abstol=>1e-25)
sim1 = analyticless_test_convergence(dts,prob,Tsit5(),test_setup)
sim2 = analyticless_test_convergence(dts,prob,Vern9(),test_setup)

@test abs(sim1.𝒪est[:final]-5) < 0.2
@test abs(sim2.𝒪est[:final]-9) < 0.2

function f2(du,u,p,t)
  du .= 1.01u
end
function g2(du,u,p,t)
  du .= 1.01u
end
prob = SDEProblem(f2,g2,[1.0;1.0],(0.0,1.0))

using StochasticDiffEq

dts = (1/2).^(7:-1:3)
test_dt = 1/2^8
Random.seed!(100)
sim1 = analyticless_test_convergence(dts,prob,SRIW1(),test_dt,trajectories=100)
@test abs(sim1.𝒪est[:final]-1.5) < 0.4
@show sim1.𝒪est[:final]

dts = (1/2).^(7:-1:4)
test_dt = 1/2^8
sim2 = analyticless_test_convergence(dts,prob,SRIW1(),test_dt,trajectories=100, use_noise_grid=false)
@test abs(sim2.𝒪est[:final]-1.5) < 0.3
@show sim2.𝒪est[:final]

# EnsembleProblem

function prob_func(prob, i, repeat)
    remake(prob,seed=seeds[i])
end

u₀ = [1.0,1.0]
function f2!(du,u,p,t)
  du[1] = -273//512*u[1]
  du[2] = -1//160*u[1]-(-785//512+sqrt(2)/8)*u[2]
end
function g2!(du,u,p,t)
  du[1,1] = 1//4*u[1]
  du[1,2] = 1//16*u[1]
  du[2,1] = (1-2*sqrt(2))/4*u[1]
  du[2,2] = 1//10*u[1]+1//16*u[2]
end
dts = 1 .//2 .^(3:-1:0)
tspan = (0.0,3.0)

h2(z) = z^2 # but apply it only to u[1]

prob = SDEProblem(f2!,g2!,u₀,tspan,noise_rate_prototype=zeros(2,2))

numtraj = Int(1e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h2(sol[end][1]),false),
        prob_func = prob_func
        )
sim = test_convergence(dts,ensemble_prob,DRI1(),save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,weak_timeseries_errors=false,weak_dense_errors=false,expected_value=exp(-3.0))

@test abs(sim.𝒪est[:weak_final]-2.0) < 0.3
@show sim.𝒪est[:weak_final]

### SDDE

function hayes_modelf(du,u,h,p,t)
    τ,a,b,c,α,β,γ = p
    du .= a.*u .+ b .* h(p,t-τ) .+ c
end
function hayes_modelg(du,u,h,p,t)
    τ,a,b,c,α,β,γ = p
    du .= α.*u .+ β.*h(p,t-τ) .+ γ
end
h(p,t) = (ones(1) .+ t);
tspan = (0.,10.)

pmul = [1.0,-4.,-2.,10.,-1.3,-1.2, 1.1]
padd = [1.0,-4.,-2.,10.,-0.0,-0.0, 0.1]

prob = SDDEProblem(hayes_modelf, hayes_modelg, [1.], h, tspan, pmul; constant_lags = (pmul[1],));
dts = (1/2).^(7:-1:3)
test_dt = 1/2^8
sim2 = analyticless_test_convergence(dts,prob,RKMil(),test_dt,trajectories=100, use_noise_grid=false)
@test abs(sim2.𝒪est[:final]-1.0) < 0.3


text:
using TopicModels

##################################################################################################################################
# Fit and Visualize Real-World Text Data

exdir = joinpath(dirname(pathof(TopicModels)), "..", "examples")

testDocuments = readDocs(open(joinpath(exdir, "cora.documents")))
testLexicon = readLexicon(open(joinpath(exdir, "cora.lexicon")))

corpus = Corpus(testDocuments,testLexicon)
model = Model(fill(0.1, 10), fill(0.01,length(testLexicon)), corpus)
state = State(model,corpus)

#@time Juno.@run trainModel(model, state, 30)
@time trainModel(model, state, 30)
topWords = topTopicWords(model, state, 10)

# visualize the fit
# using Plots, UMAP
# @time embedding = umap(state.topics, 2, n_neighbors=10)
# maxlabels = vec(map(i->i[1], findmax(state.topics,dims=1)[2]))
# scatter(embedding[1,:], embedding[2,:], zcolor=maxlabels, title="UMAP: Max-Likelihood Doc Topics on Learned", marker=(2, 2, :auto, stroke(0)))

##################################################################################################################################
# Fit, Validate, and Visualize Synthetic Data Derived from a Fully-Generative Simulation (Poisson-distributed document-length)

k = 10
lexLength = 1000
corpLambda = 1000 # poisson parameter for random doc length
corpLength = 100
scaleK = 0.01
scaleL = 0.01
testCorpus = LdaCorpus(k, lexLength, corpLambda, corpLength, scaleK, scaleL)

testModel = Model(testCorpus.alpha, testCorpus.beta, testCorpus)
testState = State(testModel, testCorpus)
    @time trainModel(testModel, testState, 100)

# compute validation metrics on a single fit
CorpusARI(testState,testModel,testCorpus) # ARI for max. likelihood. document topics
DocsARI(testState,testCorpus) # ARI for actual word topics


text:
#Written By Johnathan Bizzano
module Regression
    using ...Functions
    using ...Data
    using ..BuiltInFunctions

    include("RegressableSummation.jl")
    include("MultiVarPolynomial.jl")

end


text:
module HaltonSequences

using Primes
export haltonvalue,Halton,HaltonPoint

const MAXLEN = typemax(UInt32)

"""
    haltonvalue(n, base[, T])

Return the `n`th element of the Halton sequence with prime base `base` (defaults to 2). If
a `<:Real` type `T` is given, compute the value in that type.

# Examples
```jldoctest
julia> haltonvalue(4,3,Rational)
4//9
```
"""
function haltonvalue(n::Integer,base::Integer=2,T::Type=float(Int))
    haltonvalue!(zero(T),n,base)
end

function haltonvalue!(x::Real,n::Integer,base::Integer=2)
    x = zero(x)
    s = one(x)
    while n > 0
        s /= base
        n,r = divrem(n,base)
        x += s*r
    end
    return x
end

const val = haltonvalue  # for terser qualified uses

########################################################
# Vector of Halton sequence values
########################################################
"""
    Halton(base=2; <keyword arguments>)
    Halton{T<:Real}(base=2; <keyword arguments>)

Create a `Vector`-like object for a Halton sequence. Optionally, use a `<:Real` type
`T` for the elements.

# Arguments
- `base::Integer`: prime number base for the sequence
- `start=1` (keyword): starting/offset index
- `length=typemax(UInt32)` (keyword): length of the sequence vector

# Examples
```jldoctest
julia> Halton(3,length=4)
4-element Halton{Float64}:
 0.3333333333333333
 0.6666666666666666
 0.1111111111111111
 0.4444444444444444
```
```jldoctest
julia> Halton{Rational}(2,length=7)
7-element Halton{Rational}:
 1//2
 1//4
 3//4
 1//8
 5//8
 3//8
 7//8
```
```jldoctest
julia> sum( Halton()[1:10000] )/10000
0.4998324768066406
```

 See also: [`HaltonPoint`](@ref), [`haltonvalue`](@ref)
"""
struct Halton{T<:Real} <: AbstractArray{T, 1}
    base::Int
    start::Int
    length::Int
    # Constructor for when the type is explicit.
    function Halton{T}(base::Integer=2;start=1,length=MAXLEN) where {T<:Real}
        new(base,start,length)
    end
end

# Constructor when no type given. Also includes a check for prime base.
function Halton(base::Integer=2;kw...)
    isprime(base) || error("Base must be a prime number.")
    Halton{float(Int)}(base;kw...)
end

# Serve the AbstractArray interface
Base.size(H::Halton) = (H.length,)
Base.IndexStyle(::Type{<:Halton}) = IndexLinear()
function Base.getindex(H::Halton{T}, n::Integer) where T
    n <= H.length || throw(BoundsError(H,n))
    haltonvalue(n+H.start-1,H.base,T)
end

########################################################
# Vector of multidimensional points based on different bases
########################################################
"""
    HaltonPoint(ndim::Integer; <keyword arguments>)
    HaltonPoint(base::Vector{<:Integer}; <keyword arguments>)
    HaltonPoint{T<:Real}(base::Vector{<:Integer}; <keyword arguments>)

Create a `Vector`-like object for a sequence of `ndim`-dimensional vectors constructed
from Halton sequences. If given, `base` should be a vector of primes and then
`ndim=length(base)`. If `T` is given, it is the element type for the vectors.

# Keyword arguments
- `start=1`: starting/offset index
- `length=typemax(UInt32)`: length of the sequence vector

# Examples
```jldoctest
julia> HaltonPoint(4,length=5)
5-element HaltonPoint{Float64}:
 [0.5, 0.333333, 0.2, 0.142857]
 [0.25, 0.666667, 0.4, 0.285714]
 [0.75, 0.111111, 0.6, 0.428571]
 [0.125, 0.444444, 0.8, 0.571429]
 [0.625, 0.777778, 0.04, 0.714286]
```
```jldoctest
julia> h = HaltonPoint([5,11,17]);  sum(h[10001:20000])/10000
3-element Array{Float64,1}:
 0.49996162048000015
 0.5000650725546567
 0.49992980926952507
```

 See also: [`Halton`](@ref), [`haltonvalue`](@ref)
"""
struct HaltonPoint{T<:Real} <: AbstractArray{T, 1}
    base::Vector{Int}
    start::Int
    length::Int
    # Constructor when the type and base are explicit.
    function HaltonPoint{T}(base::Vector{S};start=1,length=MAXLEN) where {T<:Real} where {S<:Integer}
        new(base,start,length)
    end
    # Constructor when the type is explicit, but bases are not given.
    function HaltonPoint{T}(ndim=Integer;kw...) where {T<:Real}
        function firstprimes(k::Integer)
            n = zeros(UInt16,k)
            n[1] = 2
            for i = 2:k
                n[i] = nextprime(n[i-1]+1)
            end
            return n
        end
        HaltonPoint{T}(firstprimes(ndim);kw...)
    end
end

# Constructor when the type is not given.
HaltonPoint(b_or_n;kw...) = HaltonPoint{float(Int)}(b_or_n;kw...)

# Serve the AbstractArray interface
Base.size(P::HaltonPoint) = (P.length,)
Base.IndexStyle(::Type{<:HaltonPoint}) = IndexLinear()
Base.eltype(::HaltonPoint{T}) where T = Vector{T}
function Base.getindex(P::HaltonPoint{T}, n::Integer) where T
    [ haltonvalue(n+P.start-1,P.base[i],T) for i in eachindex(P.base) ]
end
# This seems to be needed for range/vector indexing, since the eltype is not scalar.
function Base.getindex(P::HaltonPoint{T}, n::Union{AbstractVector,AbstractRange}) where T
    [ P[n[i]] for i in eachindex(n) ]
end

end # module


text:
using RecipesBase
using Printf

"""
    plot(image::IntensityMap)

where `image` is templated off of EHTImage struct.

# Details
This was created to be close to the ehtim display object. It takes an
EHTImage object and plots it according to EHT conventions.

Note that is does not save the figure.
"""
@recipe function f(image::IntensityMap; uvscale=rad2μas)

    #Define some constants
    #Construct the image grid in μas
    fovx, fovy = uvscale.(fov(image))
    xitr, yitr = imagepixels(image)

    tickfontsize --> 11
    guidefontsize --> 14
    size --> (500,400)
    xaxis --> "ΔRA  (μas)"
    yaxis --> "ΔDEC (μas)"
    seriescolor --> :afmhot
    aspect_ratio --> 1
    bar_width --> 0
    xlims --> (-fovx/2,fovx/2)
    ylims --> (-fovy/2,fovy/2)
    #left_margin --> -2mm
    #right_margin --> 5mm
    z = image/uvscale.(1)^2
    seriestype := :heatmap
    #fontfamily --> "sans serif"
    colorbar_title --> "Jy/μas²"
    xflip --> true
    widen := false
    #framestyle --> :box
    title --> @sprintf("flux = %.2f Jy", flux(image))
    linecolor-->:black
    tick_direction --> :out
    uvscale.(collect(xitr)),uvscale.(collect(yitr)),z
end

"""
    plot(image::AbstractModel)

where `image` is templated off of EHTImage struct.

# Details
This was created to be close to the ehtim display object. It takes an
EHTImage object and plots it according to EHT conventions.

Note that is does not save the figure.
"""
@recipe function f(m::AbstractModel; uvscale=rad2μas,
                   fovx = 2*radialextent(m), fovy=2*radialextent(m), dim=(512, 512))

    ny, nx = dim
    image = intensitymap(m, fovx, fovy, nx, ny)
    psizex,psizey = pixelsizes(image)
    #Define some constants
    #Construct the image grid in μas
    fovx, fovy = uvscale.(fov(image))
    xitr, yitr = uvscale.(imagepixels(image))
    tickfontsize --> 11
    guidefontsize --> 14
    size --> (500,400)
    xaxis --> "ΔRA  (μas)"
    yaxis --> "ΔDEC (μas)"
    seriescolor --> :afmhot
    aspect_ratio --> 1
    bar_width --> 0
    xlims --> (-fovx/2,fovx/2)
    ylims --> (-fovy/2,fovy/2)
    #left_margin --> -2mm
    #right_margin --> 5mm
    z = image*psizex*psizey
    seriestype := :heatmap
    #fontfamily --> "sans serif"
    colorbar_title --> "Jy/px²"
    xflip --> true
    widen := false
    #framestyle --> :box
    title --> @sprintf("flux = %.2f Jy", flux(image))
    linecolor-->:black
    tick_direction --> :out
    collect(xitr),collect(yitr),z
end


text:
module GeneralFunctions
using Dierckx
using TwoFAST

export BiasDict13, BiasDictlhd, BiasDict22
export extend_func, eval_extend_func
export Plhd, Pllhd, Plhdmultipole
export ξ
export best_q, henry_q, show_best_q

include("BiasDictConstructors.jl")
include("ExtendedFunctions.jl")
include("Plhd.jl")
include("xi.jl")
include("Bestq.jl")
end

max_stars_repo_name: JosephTomlinson/GeneralBiasPk

max_issues_repo_name: JosephTomlinson/GeneralBiasPk

max_forks_repo_name: JosephTomlinson/GeneralBiasPk


text:
## Getting the ID and metadata of all food webs archived on mangal.io

## Read data from the mangal.io database

# query the ID number and other metadata for all ecological networks archived on mangal.io
# id: network ID
# S: number of species in the network
# L: number of interactions in the network
# P: number of interactions of predation in the network
# H: number of interactions of herbivory in the network

number_of_networks = count(MangalNetwork)
count_per_page = 100
number_of_pages = convert(Int, ceil(number_of_networks/count_per_page))

mangal_networks = DataFrame(fill(0, (number_of_networks, 5)),
                [:id, :S, :L, :P, :H])

global cursor = 1
@progress "Paging networks" for page in 1:number_of_pages
   global cursor
   networks_in_page = Mangal.networks("count" => count_per_page, "page" => page-1)
   @progress "Counting items" for current_network in networks_in_page
       S = count(MangalNode, current_network)
       L = count(MangalInteraction, current_network)
       P = count(MangalInteraction, current_network, "type" => "predation")
       H = count(MangalInteraction, current_network, "type" => "herbivory")
       mangal_networks[cursor,:] .= (current_network.id, S, L, P, H)
       cursor = cursor + 1
   end
end

# Filter for food webs (i.e. networks mainly having trophic links)
fw = (mangal_networks.P .+ mangal_networks.H) ./ mangal_networks.L .>= 0.5

mangal_foodwebs = mangal_networks[fw, :]

# Write file
CSV.write(joinpath("data", "mangal_foodwebs.csv"), mangal_foodwebs)

max_stars_repo_name: PoisotLab/worldwide-links-prediction

max_issues_repo_name: PoisotLab/worldwide-links-prediction

max_forks_repo_name: PoisotLab/worldwide-links-prediction


text:
module mpiFFT4jl

module slab
include("slab.jl")
end

# module pencil
# include("pencil.jl")
# end

end


text:
# This file is a part of SIS3316Digitizers.jl, licensed under the MIT License (MIT).

"""
    AbstractVMEGateway

Abstract type for VME Gateways.

Required functionality for subtypes:

* [`read_registers`](@ref)
* [`write_registers!`](@ref)
* [`read_fifo_space!`](@ref)
* [`Base.close`](@ref)
"""
abstract type AbstractVMEGateway end

export AbstractVMEGateway

@inline Base.Broadcast.broadcastable(gw::AbstractVMEGateway) = Ref(gw)

"""
    read_registers!(
        gw::AbstractVMEGateway,
        addrs::AbstractVector{<:Unsigned}, values::AbstractVector{<:Unsigned};
        timeout::Float64 = ...
    )::typeof(values)

Read VME register values via gateway `gw`.

Vectors `addrs` and `values` must have the same length. Returns `values`.

There is no guarantee on the order in which the read operations are executed.
Depending on the hardware, read operations may be executed in batches.
The `timeout` value is applied to each internal read operation.
"""
function read_registers! end
export read_registers!

"""
    write_registers!(
        gw::AbstractVMEGateway,
        addrs::AbstractVector{<:Unsigned}, values::AbstractVector{<:Unsigned};
        timeout::Float64 = ...
    )::Nothing

Write VME register values via gateway `gw`.

Vectors `addrs` and `values` must have the same length.

There is no guarantee on the order in which the write operations are executed.
Depending on the hardware, write operations may be executed in batches.
The `timeout` value is applied to each internal write operation.
"""
function write_registers! end
export write_registers!

"""
    read_bulk!(
        gw::AbstractVMEGateway,
        addr::Unsigned,
        data::AbstractVector{UInt32};
        timeout::Float64 = gw.default_timeout
    )

Fast bulk data read transfer.
"""
function read_bulk! end
export read_bulk!

"""
    struct DummyVMEGateway{Addr<:Unsigned,Value<:Unsigned} <: AbstractVMEGateway

Dummy VME Gateway.
"""
struct DummyVMEGateway{Addr<:Unsigned,Value<:Unsigned} <: AbstractVMEGateway
    state::IdDict{Addr,Value}
end

export DummyVMEGateway

function DummyVMEGateway{Addr,Value}() where {Addr<:Unsigned,Value<:Unsigned}
    DummyVMEGateway{Addr,Value}(IdDict{Addr,Value}())
end

Base.deepcopy(gw::DummyVMEGateway{Addr,Value}) where {Addr,Value} = DummyVMEGateway{Addr,Value}()

function read_registers!(
    gw::DummyVMEGateway{Addr,Value},
    addrs::AbstractVector{<:Addr},
    values::AbstractVector{<:Value};
    timeout::Float64 = 0.005
) where {Addr,Value}
    @debug "read_registers!(::DummyVMEGateway, $addrs, $values)"
    length(addrs) == length(values) || throw(ArgumentError("Number of addresses does not match number of values"))
    values .= [get(gw.state, a, zero(Value)) for a in addrs]
    values
end

function write_registers!(
    gw::DummyVMEGateway{Addr,Value},
    addrs::AbstractVector{<:Addr},
    values::AbstractVector{<:Value};
    timeout::Float64 = 0.005
) where {Addr,Value}
    @debug "write_registers!(::DummyVMEGateway, $addrs, $values)"
    length(addrs) == length(values) || throw(ArgumentError("Number of addresses does not match number of values"))
    for (a, v) in zip(addrs, values)
        gw.state[a] = v
    end
    nothing
end

Base.isopen(gw::DummyVMEGateway) = true

function Base.close(gw::DummyVMEGateway)
    empty!(gw.state)
    # @info "Closed $gw" # logging causes trouble in finalizers
end

function read_bulk!(gw::DummyVMEGateway, addr::UInt32, data::AbstractVector{UInt32})
    # @debug "read_bulk!(::DummyVMEGateway, $(repr(addr)), data)"
    for i in eachindex(data)
        data[i] = addr + i - firstindex(data)
    end
    data
end


text:
@testset "Britter-McQuaid plume tests" begin
    # Britter-McQuaid example, *Guidelines for Consequence Analysis of Chemical
    # Releases* CCPS, 1999, pg 122
    ex = Scenario( Dict([
    :mass_emission_rate => (0.23*425.6),
    :release_height => 10.0,
    :jet_density => 1.76,
    :release_temperature => (273.15-162),
    :windspeed => 10.9,
    :ambient_density => 1.224,
    :ambient_temperature => 298,
    :pasquill_gifford => "F"
    ]))
    # known answers
    # set 1 is covers the short-distance correlation
    # set 2 is the answer from the above example, which uses the interpolations
    x₁, c₁ = 2.26*20, 1.1827612798486666
    x₂, c₂ = 367.0, 0.0872919843565787

    # missing model parameters
    @test_throws MissingException plume(ambient, model=:brittermcquaid)

    @testset "Britter-McQuaid plume tests for class $class" for class in ["A", "B", "C", "D", "E", "F"]
        # because the windspeed is at 10m, the class should not impact the
        # calculations but this is a check that getting the corresponding
        # correlation doesn't throw any errors or do anything deeply strange
        s = Scenario( ex; pasquill_gifford=class )
        britter_mcquaid = plume(s, model=:brittermcquaid)
        @test britter_mcquaid(x₁,0,0) ≈ c₁
        @test britter_mcquaid(x₂,0,0) ≈ c₂

    end
end


text:
using DashBootstrapComponents, DashHtmlComponents

badges = html_span([
    dbc_badge("Primary", color="primary", className="mr-1"),
    dbc_badge("Secondary", color="secondary", className="mr-1"),
    dbc_badge("Success", color="success", className="mr-1"),
    dbc_badge("Warning", color="warning", className="mr-1"),
    dbc_badge("Danger", color="danger", className="mr-1"),
    dbc_badge("Info", color="info", className="mr-1"),
    dbc_badge("Light", color="light", className="mr-1"),
    dbc_badge("Dark", color="dark"),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using StaticArrays
using ClimateMachine.VariableTemplates
using ClimateMachine.DGMethods: nodal_update_auxiliary_state!
import ClimateMachine.DGMethods:
    BalanceLaw,
    vars_state_auxiliary,
    vars_state_conservative,
    vars_state_gradient,
    vars_state_gradient_flux,
    flux_first_order!,
    flux_second_order!,
    source!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    init_state_auxiliary!,
    update_auxiliary_state!,
    init_state_conservative!,
    boundary_state!,
    wavespeed,
    LocalGeometry,
    number_state_conservative,
    number_state_gradient
using ClimateMachine.DGMethods.NumericalFluxes:
    NumericalFluxFirstOrder, NumericalFluxSecondOrder, NumericalFluxGradient
import ClimateMachine.DGMethods.NumericalFluxes:
    numerical_flux_first_order!, boundary_flux_second_order!

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

abstract type AdvectionDiffusionProblem end
struct AdvectionDiffusion{dim, P, fluxBC, no_diffusion} <: BalanceLaw
    problem::P
    function AdvectionDiffusion{dim}(
        problem::P,
    ) where {dim, P <: AdvectionDiffusionProblem}
        new{dim, P, false, false}(problem)
    end
    function AdvectionDiffusion{dim, fluxBC}(
        problem::P,
    ) where {dim, P <: AdvectionDiffusionProblem, fluxBC}
        new{dim, P, fluxBC, false}(problem)
    end
    function AdvectionDiffusion{dim, fluxBC, no_diffusion}(
        problem::P,
    ) where {dim, P <: AdvectionDiffusionProblem, fluxBC, no_diffusion}
        new{dim, P, fluxBC, no_diffusion}(problem)
    end
end

# Stored in the aux state are:
#   `coord` coordinate points (needed for BCs)
#   `u` advection velocity
#   `D` Diffusion tensor
vars_state_auxiliary(::AdvectionDiffusion, FT) =
    @vars(coord::SVector{3, FT}, u::SVector{3, FT}, D::SMatrix{3, 3, FT, 9})
function vars_state_auxiliary(
    ::AdvectionDiffusion{dim, P, fluxBC, true},
    FT,
) where {dim, P, fluxBC}
    @vars begin
        coord::SVector{3, FT}
        u::SVector{3, FT}
    end
end

# Density is only state
vars_state_conservative(::AdvectionDiffusion, FT) = @vars(ρ::FT)

# Take the gradient of density
vars_state_gradient(::AdvectionDiffusion, FT) = @vars(ρ::FT)
vars_state_gradient(
    ::AdvectionDiffusion{dim, P, fluxBC, true},
    FT,
) where {dim, P, fluxBC} = @vars()

# The DG auxiliary variable: D ∇ρ
vars_state_gradient_flux(::AdvectionDiffusion, FT) = @vars(σ::SVector{3, FT})
vars_state_gradient_flux(
    ::AdvectionDiffusion{dim, P, fluxBC, true},
    FT,
) where {dim, P, fluxBC} = @vars()

"""
    flux_first_order!(m::AdvectionDiffusion, flux::Grad, state::Vars,
                       aux::Vars, t::Real)

Computes non-diffusive flux `F` in:

```
∂ρ
-- = - ∇ • (u ρ - σ) = - ∇ • F
∂t
```
Where

 - `u` is the advection velocity
 - `ρ` is the advected quantity
 - `σ` is DG auxiliary variable (`σ = D ∇ ρ` with D being the diffusion tensor)
"""
function flux_first_order!(
    ::AdvectionDiffusion,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρ = state.ρ
    u = aux.u
    flux.ρ += u * ρ
end

"""
flux_second_order!(m::AdvectionDiffusion, flux::Grad, auxDG::Vars)

Computes diffusive flux `F` in:

```
∂ρ
-- = - ∇ • (u ρ - σ) = - ∇ • F
∂t
```
Where

 - `u` is the advection velocity
 - `ρ` is the advected quantity
 - `σ` is DG auxiliary variable (`σ = D ∇ ρ` with D being the diffusion tensor)
"""
function flux_second_order!(::AdvectionDiffusion, flux::Grad, auxDG::Vars)
    σ = auxDG.σ
    flux.ρ += -σ
end
flux_second_order!(
    ::AdvectionDiffusion{dim, P, fluxBC, true},
    flux::Grad,
    auxDG::Vars,
) where {dim, P, fluxBC} = nothing
flux_second_order!(
    m::AdvectionDiffusion,
    flux::Grad,
    state::Vars,
    auxDG::Vars,
    auxHDG::Vars,
    aux::Vars,
    t::Real,
) = flux_second_order!(m, flux, auxDG)

"""
    compute_gradient_argument!(m::AdvectionDiffusion, transform::Vars, state::Vars,
                   aux::Vars, t::Real)

Set the variable to take the gradient of (`ρ` in this case)
"""
function compute_gradient_argument!(
    m::AdvectionDiffusion,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.ρ = state.ρ
end

"""
    compute_gradient_flux!(m::AdvectionDiffusion, transform::Vars, gradvars::Vars,
               aux::Vars)

Set the variable to take the gradient of (`ρ` in this case)
"""
function compute_gradient_flux!(
    m::AdvectionDiffusion,
    auxDG::Vars,
    gradvars::Grad,
    aux::Vars,
)
    ∇ρ = gradvars.ρ
    D = aux.D
    auxDG.σ = D * ∇ρ
end
compute_gradient_flux!(
    m::AdvectionDiffusion,
    auxDG::Vars,
    gradvars::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) = compute_gradient_flux!(m, auxDG, gradvars, aux)

"""
    source!(m::AdvectionDiffusion, _...)

There is no source in the advection-diffusion model
"""
source!(m::AdvectionDiffusion, _...) = nothing

"""
    wavespeed(m::AdvectionDiffusion, nM, state::Vars, aux::Vars, t::Real)

Wavespeed with respect to vector `nM`
"""
function wavespeed(m::AdvectionDiffusion, nM, state::Vars, aux::Vars, t::Real)
    u = aux.u
    abs(dot(nM, u))
end

"""
    init_state_auxiliary!(m::AdvectionDiffusion, aux::Vars, geom::LocalGeometry)

initialize the auxiliary state
"""
function init_state_auxiliary!(
    m::AdvectionDiffusion,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.coord = geom.coord
    init_velocity_diffusion!(m.problem, aux, geom)
end

has_variable_coefficients(::AdvectionDiffusionProblem) = false
function update_auxiliary_state!(
    dg::DGModel,
    m::AdvectionDiffusion,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    if has_variable_coefficients(m.problem)
        nodal_update_auxiliary_state!(dg, m, Q, t, elems) do m, state, aux, t
            update_velocity_diffusion!(m.problem, m, state, aux, t)
        end
        return true
    end
    return false
end

function init_state_conservative!(
    m::AdvectionDiffusion,
    state::Vars,
    aux::Vars,
    coords,
    t::Real,
)
    initial_condition!(m.problem, state, aux, coords, t)
end

Neumann_data!(problem, ∇state, aux, x, t) = nothing
Dirichlet_data!(problem, state, aux, x, t) = nothing

function boundary_state!(
    nf,
    m::AdvectionDiffusion,
    stateP::Vars,
    auxP::Vars,
    nM,
    stateM::Vars,
    auxM::Vars,
    bctype,
    t,
    _...,
)
    if bctype == 1 # Dirichlet
        Dirichlet_data!(m.problem, stateP, auxP, auxP.coord, t)
    elseif bctype ∈ (2, 4) # Neumann
        stateP.ρ = stateM.ρ
    elseif bctype == 3 # zero Dirichlet
        stateP.ρ = 0
    end
end

function boundary_state!(
    nf::CentralNumericalFluxSecondOrder,
    m::AdvectionDiffusion,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n⁻::SVector,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
    _...,
)

    if bctype ∈ (1, 3) # Dirchlet
        # Just use the minus side values since Dirchlet
        diff⁺.σ = diff⁻.σ
    elseif bctype == 2 # Neumann with data
        FT = eltype(diff⁺)
        ngrad = number_state_gradient(m, FT)
        ∇state = Grad{vars_state_gradient(m, FT)}(similar(
            parent(diff⁺),
            Size(3, ngrad),
        ))
        # Get analytic gradient
        Neumann_data!(m.problem, ∇state, aux⁻, aux⁻.coord, t)
        compute_gradient_flux!(m, diff⁺, ∇state, aux⁻)
        # compute the diffusive flux using the boundary state
    elseif bctype == 4 # zero Neumann
        FT = eltype(diff⁺)
        ngrad = number_state_gradient(m, FT)
        ∇state = Grad{vars_state_gradient(m, FT)}(similar(
            parent(diff⁺),
            Size(3, ngrad),
        ))
        # Get analytic gradient
        ∇state.ρ = SVector{3, FT}(0, 0, 0)
        # convert to auxDG variables
        compute_gradient_flux!(m, diff⁺, ∇state, aux⁻)
    end
    nothing
end
boundary_state!(
    nf::CentralNumericalFluxSecondOrder,
    m::AdvectionDiffusion{dim, P, fluxBC, true},
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n⁻::SVector,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
    _...,
) where {dim, P, fluxBC} = nothing

function boundary_flux_second_order!(
    nf::CentralNumericalFluxSecondOrder,
    m::AdvectionDiffusion{dim, P, true},
    F,
    state⁺,
    diff⁺,
    hyperdiff⁺,
    aux⁺,
    n⁻,
    state⁻,
    diff⁻,
    hyperdiff⁻,
    aux⁻,
    bctype,
    t,
    _...,
) where {dim, P}

    # Default initialize flux to minus side
    if bctype ∈ (1, 3) # Dirchlet
        # Just use the minus side values since Dirchlet
        flux_second_order!(m, F, state⁻, diff⁻, hyperdiff⁻, aux⁻, t)
    elseif bctype == 2 # Neumann data
        FT = eltype(diff⁺)
        ngrad = number_state_gradient(m, FT)
        ∇state = Grad{vars_state_gradient(m, FT)}(similar(
            parent(diff⁺),
            Size(3, ngrad),
        ))
        # Get analytic gradient
        Neumann_data!(m.problem, ∇state, aux⁻, aux⁻.coord, t)
        # get the diffusion coefficient
        D = aux⁻.D
        # exact the exact data
        ∇ρ = ∇state.ρ
        # set the flux
        F.ρ = -D * ∇ρ
    elseif bctype == 4 # Zero Neumann
        FT = eltype(diff⁺)
        F.ρ = SVector{3, FT}(0, 0, 0)
    end
    nothing
end
boundary_flux_second_order!(
    ::CentralNumericalFluxSecondOrder,
    ::AdvectionDiffusion{dim, P, true, true},
    _...,
) where {dim, P} = nothing

struct UpwindNumericalFlux <: NumericalFluxFirstOrder end
function numerical_flux_first_order!(
    ::UpwindNumericalFlux,
    ::AdvectionDiffusion,
    fluxᵀn::Vars{S},
    n::SVector,
    state⁻::Vars{S},
    aux⁻::Vars{A},
    state⁺::Vars{S},
    aux⁺::Vars{A},
    t,
) where {S, A}
    un⁻ = dot(n, aux⁻.u)
    un⁺ = dot(n, aux⁺.u)
    un = (un⁺ + un⁻) / 2

    fluxᵀn.ρ = un ≥ 0 ? un * state⁻.ρ : un * state⁺.ρ
end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
global anyerrors = false

tests = [
   "mock test" => "00_allgood.jl",
   "Linear gradient" => "01_gradient.jl",
   "Boundaries detection" => "02_boundaries.jl",
   "Triangulation wombling" => "03_delaunay.jl",
   "Direction of change" => "04_direction.jl",
   "REQUIRE: SimpleSDMLayers" => "R1_simplesdmlayers.jl"
]

for test in tests
   try
      include(joinpath("units", test.second))
      println("\033[1m\033[32m✓\033[0m\t$(test.first)")
   catch e
      global anyerrors = true
      println("\033[1m\033[31m×\033[0m\t$(test.first)")
      println("\033[1m\033[38m→\033[0m\ttest/$(test.second)")
      showerror(stdout, e, backtrace())
      println()
      break
   end
end

if anyerrors
   throw("Tests failed")
end


text:
using LinearAlgebra: lmul!, norm
using Parameters: Parameters
using NQCDistributions: ElectronicDistribution

"""
    eCMM{T} <: DynamicsMethods.Method

# References

- [HeGong2021](@cite)
- [HeWu2021](@cite)
"""
struct eCMM{T} <: DynamicsMethods.Method
    γ::T
    temp_q::Vector{T}
    temp_p::Vector{T}
    function eCMM{T}(n_states::Integer, γ) where {T}
        new{T}(γ, zeros(n_states), zeros(n_states))
    end
end

function NQCDynamics.Simulation{eCMM}(atoms::Atoms{T}, model::Model; γ=0, kwargs...) where {T}
    NQCDynamics.Simulation(atoms, model, eCMM{T}(NQCModels.nstates(model), γ); kwargs...)
end

function DynamicsMethods.DynamicsVariables(sim::Simulation{<:eCMM}, v, r, ::ElectronicDistribution)

    F = NQCModels.nstates(sim)
    radius = sqrt(2 + 2F*sim.method.γ)
    points = generate_random_points_on_nsphere(2F, radius)
    qmap = points[begin:div(F,2)+1]
    pmap = points[div(F,2)+2:end]

    ComponentVector(v=v, r=r, pmap=pmap, qmap=qmap)
end

function generate_random_points_on_nsphere(n, radius)
    randomnormal = randn(n)
    return randomnormal ./ norm(randomnormal) .* radius
end

function DynamicsMethods.motion!(du, u, sim::AbstractSimulation{<:eCMM}, t)
    dr = DynamicsUtils.get_positions(du)
    dv = DynamicsUtils.get_velocities(du)
    r = DynamicsUtils.get_positions(u)
    v = DynamicsUtils.get_velocities(u)

    DynamicsUtils.velocity!(dr, v, r, sim, t)
    acceleration!(dv, u, sim)
    set_mapping_force!(du, u, sim)
end

function acceleration!(dv, u, sim::Simulation{<:eCMM})

    Parameters.@unpack γ, temp_q, temp_p = sim.method

    Calculators.evaluate_derivative!(sim.calculator, DynamicsUtils.get_positions(u))
    qmap = get_mapping_positions(u)
	pmap = get_mapping_momenta(u)
    for I in eachindex(dv)
        D = sim.calculator.derivative[I]
        mul!(temp_q, D, qmap)
        mul!(temp_p, D, pmap)
        dv[I] = dot(qmap, temp_q)
        dv[I] += dot(pmap, temp_p)
        dv[I] -= tr(D) * 2γ
    end
    lmul!(-1/2, dv)
    DynamicsUtils.divide_by_mass!(dv, sim.atoms.masses)
end

function set_mapping_force!(du, u, sim::Simulation{<:eCMM})
    Calculators.evaluate_potential!(sim.calculator, DynamicsUtils.get_positions(u))
    V = sim.calculator.potential
    mul!(get_mapping_positions(du), V, get_mapping_momenta(u))
    mul!(get_mapping_momenta(du), V, get_mapping_positions(u))
    lmul!(-1, get_mapping_momenta(du))
end

function DynamicsUtils.classical_hamiltonian(sim::Simulation{<:eCMM}, u)
    r = DynamicsUtils.get_positions(u)
    v = DynamicsUtils.get_velocities(u)

    kinetic = DynamicsUtils.classical_kinetic_energy(sim, v)

    Calculators.evaluate_potential!(sim.calculator, r)
    V = sim.calculator.potential
    qmap = get_mapping_positions(u)
    pmap = get_mapping_momenta(u)
    potential = (pmap'V*pmap + qmap'V*qmap - 2sim.method.γ*tr(V)) / 2

    H = kinetic + potential
    return H
end

mapping_kernel(qmap, pmap, γ) = (qmap.^2 .+ pmap.^2)./2 .- γ

function inverse_mapping_kernel(qmap, pmap, γ)
    F = length(qmap)
    prefactor = (1 + F) / (1 + F*γ)^2
    subtractor = (1 - γ) / (1 + F*γ)
    return prefactor .* (qmap.^2 .+ pmap.^2)./2 .- subtractor
end

function Estimators.diabatic_population(sim::Simulation{<:eCMM}, u)
    qmap = get_mapping_positions(u)
    pmap = get_mapping_momenta(u)
    inverse_mapping_kernel(qmap, pmap, sim.method.γ)
end

function Estimators.initial_diabatic_population(sim::Simulation{<:eCMM}, u)
    qmap = get_mapping_positions(u)
    pmap = get_mapping_momenta(u)
    mapping_kernel(qmap, pmap, sim.method.γ)
end

function TimeCorrelationFunctions.evaluate_normalisation(
    sim::AbstractSimulation{<:eCMM},
    ::TimeCorrelationFunctions.PopulationCorrelationFunction
)
    NQCModels.nstates(sim)
end


text:
using MAGMA
using Test, LinearAlgebra, Random

@testset "singular_value_decomposition" begin
    include("svds/runtests.jl")
end


text:
# GUI event loops and toolkit integration for Python, most importantly
# to support plotting in a window without blocking.  Currently, we
# support wxWidgets, Qt4, and GTK+.

############################################################################

# global variable to specify default GUI toolkit to use
gui = :default # one of :default, :wx, :qt, or :gtk

pyexists(mod) = try
    pyimport(mod)
    true
  catch
    false
  end

pygui_works(gui::Symbol) = gui == :default ||
    ((gui == :wx && pyexists("wx")) ||
     (gui == :gtk && pyexists("gtk")) ||
     (gui == :qt && (pyexists("PyQt4") || pyexists("PySide"))))

# get or set the default GUI; doesn't affect running GUI
function pygui()
    global gui
    if gui::Symbol != :default && pygui_works(gui::Symbol)
        return gui::Symbol
    else
        for g in (:qt, :wx, :gtk)
            if pygui_works(g)
                gui::Symbol = g
                return gui::Symbol
            end
        end
        error("No supported Python GUI toolkit is installed.")
    end
end
function pygui(g::Symbol)
    global gui
    if g != gui::Symbol
        if g != :wx && g != :gtk && g != :qt && g != :default
            throw(ArgumentError("invalid gui $g"))
        elseif !pygui_works(g)
            error("Python GUI toolkit for $g is not installed.")
        end
        gui::Symbol = g
    end
    return g
end

############################################################################
# Event loops for various toolkits.

# call doevent(status) every sec seconds
function install_doevent(doevent::Function, sec::Real)
    timeout = Base.TimeoutAsyncWork(doevent)
    Base.start_timer(timeout,sec,sec)
    return timeout
end

# GTK:
function gtk_eventloop(sec::Real=50e-3)
    gtk = pyimport("gtk")
    events_pending = gtk["events_pending"]
    main_iteration = gtk["main_iteration"]
    function doevent(async, status::Int32) # handle all pending
        while pycall(events_pending, Bool)
            pycall(main_iteration, PyObject)
        end
    end
    install_doevent(doevent, sec)
end

# Qt4: (PyQt4 or PySide module)
function qt_eventloop(QtModule="PyQt4", sec::Real=50e-3)
    QtCore = pyimport("$QtModule.QtCore")
    instance = QtCore["QCoreApplication"]["instance"]
    AllEvents = QtCore["QEventLoop"]["AllEvents"]
    processEvents = QtCore["QCoreApplication"]["processEvents"]
    maxtime = PyObject(50)
    function doevent(async, status::Int32)
        app = pycall(instance, PyObject)
        if app.o != pynothing::PyPtr
            app["_in_event_loop"] = true
            pycall(processEvents, PyObject, AllEvents, maxtime)
        end
    end
    install_doevent(doevent, sec)
end

# wx:  (based on IPython/lib/inputhookwx.py, which is 3-clause BSD-licensed)
function wx_eventloop(sec::Real=50e-3)
    wx = pyimport("wx")
    GetApp = wx["GetApp"]
    EventLoop = wx["EventLoop"]
    EventLoopActivator = wx["EventLoopActivator"]
    function doevent(async, status::Int32)
        app = pycall(GetApp, PyObject)
        if app.o != pynothing::PyPtr
            app["_in_event_loop"] = true
            evtloop = pycall(EventLoop, PyObject)
            ea = pycall(EventLoopActivator, PyObject, evtloop)
            Pending = evtloop["Pending"]
            Dispatch = evtloop["Dispatch"]
            while pycall(Pending, Bool)
                pycall(Dispatch, PyObject)
            end
            pydecref(ea) # deactivate event loop
            pycall(app["ProcessIdle"], PyObject)
        end
    end
    install_doevent(doevent, sec)
end

# cache running event loops (so that we don't start any more than once)
const eventloops = (Symbol=>TimeoutAsyncWork)[]

function pygui_start(gui::Symbol=pygui(), sec::Real=50e-3)
    pygui(gui)
    if !haskey(eventloops, gui)
        if gui == :wx
            eventloops[gui] = wx_eventloop(sec)
        elseif gui == :gtk
            eventloops[gui] = gtk_eventloop(sec)
        elseif gui == :qt
            try
                eventloops[gui] = qt_eventloop("PyQt4", sec)
            catch
                eventloops[gui] = qt_eventloop("PySide", sec)
            end
        else
            throw(ArgumentError("unsupported GUI type $gui"))
        end
    end
    gui
end

function pygui_stop(gui::Symbol=pygui())
    if haskey(eventloops, gui)
        Base.stop_timer(pop!(eventloops, gui))
        true
    else
        false
    end
end

pygui_stop_all() = for gui in keys(eventloops); pygui_stop(gui); end

############################################################################


text:
type Split
    account::Account
    netdebit::Monetary
end
debit(split::Split) = split.netdebit
credit(split::Split) = -split.netdebit
currency(split::Split) = currency(split.netdebit)

type Entry
    id::Int64
    date::Date
    description::Compat.UTF8String
    splits::Vector{Split}

    function Entry(
        date::Date,
        description::Compat.UTF8String,
        splits::Vector{Split})

        dd = DefaultDict(Symbol, Int64, 0)
        for s in splits
            dd[currency(s)] += debit(s).val
        end
        if all(kv -> last(kv) == 0, dd)
            new(-1, date, description, splits)
        else
            throw(UnbalancedException())
        end
    end
end

id(e::Entry) = e.id
id!(e::Entry, o::Int64) = (e.id = o)
date(e::Entry) = e.date
isless(e::Entry, f::Entry) = (e.date, e.id) < (f.date, f.id)


text:
import RRSTExperiments: InputEdge

function opposite_node(edge::Edge, node)
  @assert(Graphs.src(edge) == node || Graphs.dst(edge) == node)
  return Graphs.src(edge) == node ? Graphs.dst(edge) : Graphs.src(edge)
end

"""
    get_edge_path(tree, edge)

Finds unique path in spanning tree `tree` between endpoints of `edge`and returns
list edges on this path (if `edge ∈ tree`, then `[edge]` is returned).
"""
function get_edge_path(tree::Array{Edge{Int64}}, edge::Edge{Int64})
  source = Graphs.src(edge)
  target = Graphs.dst(edge)

  queue = Tuple{Int, Array{Int}}[] # Breadth-first search
  push!(queue, (source, Int[]))
  while !isempty(queue)
    node, path = pop!(queue)
    indices = findall(e -> Graphs.src(e) == node || Graphs.dst(e) == node, tree)
    for index in indices
      if index ∉ path
        new_source = opposite_node(tree[index], node)
        path_copy = [copy(path); index]
        if (new_source == target)
          return [tree[i] for i in path_copy]
        else
          push!(queue, (new_source, path_copy))
        end
      end
    end
  end
end

"""
    get_admissible_nodes(y, nodes, edges)

Returns a list of addmissible nodes, which are:

* each `v ∈ y` is admissible node
* each `v ∈ nodes` reachible from any node in `y` by a directed path formed
of edges `e ∈ edges` is admissible
"""
function get_admissible_nodes(m::Int64, y::Array{Int64}, nodes::Array{Int64}, edges::Array{Edge{Int64}})
  queue = copy(y) # every e ∈ y is admissible

  visited = falses(m)
  for i in y
    visited[i] = true
  end

  while !isempty(queue)
    v_1 = pop!(queue)
    indices = findall(e -> Graphs.src(e) == v_1, edges)
    for i in indices
         v_2 = opposite_node(edges[i], v_1)
         if !visited[v_2]
           push!(queue, v_2)
           visited[v_2] = true
         end
    end
  end

  return filter(v -> visited[v], nodes)
end

"""
    build_admissible_graph(g, t_x, w_1, t_y, w_2)

Builds admissible graph for graph `g`, using minimum spanning trees `t_x`
(with cost vector `w_1`) and `t_y` (with cost vector `w_2`)
"""
function build_admissible_graph(edge_indices::Dict{Edge{Int64}, Int64}, t_x::Array{Edge{Int64}}, w_1::Array{Float64}, t_y::Array{Edge{Int64}}, w_2::Array{Float64})
  # First step
  v_0 = values(edge_indices) |> collect
  e_0 = Edge{Int64}[]
  indices = [(e, w_1[edge_indices[e]], w_2[edge_indices[e]]) for e in keys(edge_indices)]
  for (e, w1_e, w2_e) in indices

    if e ∉ t_x
      for (f, w1_f, _) in indices
        if f ∈ get_edge_path(t_x, e) && w1_e == w1_f # && C_e^* == C_f^*
          push!(e_0, Edge(edge_indices[e], edge_indices[f]))
        end
      end
    end

    if e ∉ t_y
      for (f, _, w2_f) in indices
        if f ∈ get_edge_path(t_y, e) && w2_e == w2_f # && \overline{c}_e^* == \overline{c}_f^*
          push!(e_0, Edge(edge_indices[f], edge_indices[e]))
        end
      end
    end
  end

  # Second step
  y = [edge_indices[e] for e in setdiff(t_y, t_x)]
  admissible_nodes = get_admissible_nodes(length(edge_indices), y, v_0, e_0)
  v_0 = filter(v -> v ∈ admissible_nodes, v_0)
  e_0 = filter(e -> Graphs.src(e) ∈ admissible_nodes && Graphs.dst(e) ∈ admissible_nodes, e_0)

  return (v_0, e_0)
end

"""
    find_δ_star(g, t_x, w_1, t_y, w_2)

Finds δ*, which is the smallest value of δ for which inequality originaly not
tight becomes tight
"""
function find_δ_star(edge_indices::Dict{Edge{Int64}, Int64}, t_x, w_1, t_y, w_2)
  δ = Inf;

  edges = keys(edge_indices)
  for e₁ in edges
    for e₂ in edges

      if e₁ ∉ t_x && e₂ ∈ get_edge_path(t_x, e₁) && 0 < w_1[edge_indices[e₁]] - w_1[edge_indices[e₂]] < δ
        δ = w_1[edge_indices[e₁]] - w_1[edge_indices[e₂]]
      end

      if e₂ ∉ t_y && e₁ ∈ get_edge_path(t_y, e₂) && 0 < w_2[edge_indices[e₂]] - w_2[edge_indices[e₁]] < δ
        δ = w_2[edge_indices[e₂]] - w_2[edge_indices[e₁]]
      end
    end
  end

  @assert(δ > 0)
  return δ
end

"""
    modify_costs_with_δ!(edge_indices, v_0, δ_star, w1_star, w2_star)

Modifies vectors of costs `w1_star` and `w2_star` using value of `δ_star`according to rules:

* if some edge is in `v_0`, then w1_star(δ) = w1_star - δ, w2_star(δ) = w2_star
* if some edge is not in `v_0`, then w1_star(δ) = w1_star, w2_star(δ) = w2_star-δ
"""
function modify_costs_with_δ!(edge_indices::Dict{Edge{Int64}, Int64}, v_0, δ_star, w1_star, w2_star)
  for e in keys(edge_indices)
    i = edge_indices[e]
    if i ∈ v_0
      w1_star[i] -= δ_star
    else
      w2_star[i] -= δ_star
    end
  end

  return (w1_star, w2_star)
end

"""
    get_path(x, y, admissible_graph)

Finds a path from one of the edges in `y` to one of the edges in `x` using
an admissible graph `admissible_graph`
"""
function get_path(x::Array, y::Array, v_a::Array{Int64}, e_a::Array{Edge{Int64}})
  queue = Tuple{Int, Array{Int}}[]

  # Put all edges of Y to queue
  for e in y
    push!(queue, (e, [e]))
  end

  shortest_path = nothing
  while !isempty(queue)
    (v, path) = pop!(queue)
    if v ∈ v_a
      next_edges = filter(i -> Graphs.dst(e_a[i]) ∉ path,
          findall(x -> Graphs.src(x) == v, e_a))
      for e_n in next_edges
        v_n = opposite_node(e_a[e_n], v)
        path_copy = [copy(path); v_n]

        # path found
        if v_n ∈ x && (shortest_path === nothing || length(path_copy) < length(shortest_path))
          shortest_path = path_copy
        else
          push!(queue, (v_n, path_copy))
        end
      end
    end
  end

  return shortest_path
end

function get_adjacent_edges(g, v; directed = false)
  edges = Graphs.Edge[]
  for e_x in Graphs.edges(g) #visit neighbours
    if Graphs.src(e_x) == v || (Graphs.dst(e_x) == v && !directed)
      push!(edges, e_x)
    end
  end
  return edges
end

function is_acyclic(g, t; directed = false)
  # Workaround for https://github.com/JuliaLang/Graphs.jl/issues/220
  edges = Graphs.Edge[]
  for e in t
    push!(edges, Graphs.Edge(Graphs.src(e), Graphs.dst(e)))
  end
  tree_graph = Graphs.graph(collect(Graphs.vertices(g)), edges)

  if directed # if graph is directed use Graphs library
    return !Graphs.test_cyclic_by_dfs(tree_graph)
  else # custom implementation of DFS otherwise
    visited_vertices = falses(Graphs.nv(tree_graph))
    visited_edges = falses(Graphs.ne(tree_graph))

    for v in Graphs.vertices(tree_graph)
      if visited_vertices[v]
        continue
      end

      edges_to_visit = filter(e -> !visited_edges[Graphs.edge_index(e)], get_adjacent_edges(tree_graph, v, directed=directed))
      stack = [(v, edges_to_visit)]

      while !isempty(stack)
        (v_c, edges) = pop!(stack)
        if visited_vertices[v_c]
          return false
        end

        visited_vertices[v_c] = true
        for e_c in edges
          visited_edges[Graphs.edge_index(e_c)] = true
          v_n = opposite_node(e_c, v_c)
          edges_to_visit = filter(e -> !visited_edges[Graphs.edge_index(e)], get_adjacent_edges(tree_graph, v_n, directed=directed))
          push!(stack, (v_n, edges_to_visit))
        end
      end
    end

    return true
  end
end

function update_trees(edge_indices::Dict{Edge{Int64}, Int64},
                      t_x::Array{Edge{Int64}},
                      w_1::Array{Float64},
                      t_y::Array{Edge{Int64}},
                      w_2::Array{Float64},
                      path::Array{Edge{Int64}},
                      z::Array{Int64})
  path_len = length(path)
  if path_len == 2
    e = path[1] # e ∈ Y
    f = path[2] # f ∈ X
    i_e = edge_indices[e]
    i_f = edge_indices[f]

    can_modify_x = e ∉ t_x && f ∈ get_edge_path(t_x, e)
    can_modify_y = f ∉ t_y && e ∈ get_edge_path(t_y, f)

    if can_modify_x && (!can_modify_y || w_1[i_e] - w_1[i_f] < w_2[i_f] - w_2[i_e])
      t_x = filter(edge -> edge != f, t_x)
      push!(t_x, e)
    else
      t_y = filter(edge -> edge != e, t_y)
      push!(t_y, f)
    end

  elseif edge_indices[path[2]] ∈ z # case 2
    e_x = path[1:2:path_len-1]
    g_x = path[2:2:path_len]
    e_y = Array{Graphs.Edge}[]
    g_y = Array{Graphs.Edge}[]
    if isodd(path_len) # case (a)
      e_y = path[3:2:path_len]
      g_y = path[2:2:path_len-1]
    else # case (b)
      e_y = path[3:2:path_len-1]
      g_y = path[2:2:path_len-1]
    end

    t_x = filter(e -> e ∉ g_x, t_x)
    t_x = t_x ∪ e_x

    t_y = filter(e -> e ∉ g_y, t_y)
    t_y = t_y ∪ e_y
  else #case 3
    e_x = Array{Graphs.Edge}[]
    g_x = Array{Graphs.Edge}[]
    e_y = path[1:2:path_len-1]
    g_y = path[2:2:path_len]
    if isodd(path_len) # case (a)
      e_x = path[3:2:path_len]
      g_x = path[2:2:path_len-1]
    else # case (b)
      e_x = path[3:2:path_len-1]
      g_x = path[2:2:path_len-1]
    end

    t_x = filter(e -> e ∉ e_x, t_x)
    t_x = t_x ∪ g_x

    t_y = filter(e -> e ∉ e_y, t_y)
    t_y = t_y ∪ g_y
  end

  return t_x, t_y
end

function get_objective_value(edge_indices::Dict{Edge{Int64}, Int64},
                            t_x::Array{Edge{Int64}},
                            w_1::Array{Float64},
                            t_y::Array{Edge{Int64}},
                            w_2::Array{Float64})
  sum = 0;
  for (e, i) in edge_indices
    if e ∈ t_x
      sum += w_1[i]
    end
    if e ∈ t_y
      sum += w_2[i]
    end
  end
  return sum
end

function get_initial_trees(n::Int, A::Array{InputEdge})
  g = SimpleGraph(n)
  for e in A
    add_edge!(g, e.i, e.j)
  end

  w1_mat = spzeros(Float64, n, n)
  w2_mat = spzeros(Float64, n, n)
  for e in A
    w1_mat[e.i, e.j] = e.C
    w1_mat[e.j, e.i] = e.C

    w2_mat[e.i, e.j] = e.c + e.d
    w2_mat[e.j, e.i] = e.c + e.d
  end

  t_x = kruskal_mst(g, w1_mat)
  t_y = kruskal_mst(g, w2_mat)

  return t_x, t_y
end

function get_first_stage_solution(t_x::Array{Edge{Int64}})::Vector{Tuple{Int, Int}}
  return [(e.src, e.dst) for e ∈ t_x]
end

# Exported function
function solve_rec_st_with_algorithm(n::Int, A::Array{InputEdge}, k::Int)
  L = n - k - 1 # |V| - K - 1, K - recovery parameter

  w1 = [e.C for e in A]
  w2 = [e.c + e.d for e in A]

  t_x, t_y = get_initial_trees(n, A)

  edges = [e.i < e.j ? Edge(e.i, e.j) : Edge(e.j, e.i) for e in A]
  edge_indices = Dict([e => i for (i, e) in enumerate(edges)])

  w1_star, w2_star = copy(w1), copy(w2)

  while length(t_x ∩ t_y) < L
    y = map(e -> edge_indices[e], setdiff(t_y, t_x))
    x = map(e -> edge_indices[e], setdiff(t_x, t_y))
    z = map(e -> edge_indices[e], t_x ∩ t_y)
    w = map(e -> edge_indices[e], setdiff(setdiff(edges, t_x), t_y))

    if isempty(x) && isempty(y) # no further improvements could be done
      break
    end

    # Check if X ∩ V^0 ≠ ∅
    v_0, e_0 = build_admissible_graph(edge_indices, t_x, w1_star, t_y, w2_star)
    while isempty(filter(e -> e ∈ v_0, x)) # no path from Y to X
      # Find δ*
      δ_star = find_δ_star(edge_indices, t_x, w1_star, t_y, w2_star)

      # Modify costs
      w1_star, w2_star = modify_costs_with_δ!(edge_indices, v_0, δ_star, w1_star, w2_star)

      # Add some new nodes to admission graph
      v_0, e_0 = build_admissible_graph(edge_indices, t_x, w1_star, t_y, w2_star)
      # repeat until there is a path from y to x
    end

    # Theorem 4: ∃(T_X', T_Y') satysfying SSOC for θ and |Z'|=|Z|+1
    # Find path from Y to X first
    path_indices = get_path(x, y, v_0, e_0)
    @assert(path_indices !== nothing)
    path = [findfirst(x -> x == i, edge_indices) for i ∈ path_indices]

    # Next modify t_x and t_y as described in the proof of Theorem 4
    t_x, t_y = update_trees(edge_indices, t_x, w1_star, t_y, w2_star, path, z)
  end

  objective_value = get_objective_value(edge_indices, t_x, w1, t_y, w2)
  tree = get_first_stage_solution(t_x)
  return objective_value, tree
end


text:
import GaussQuadrature
using OffsetArrays
using PyPlot
import FractionalTimeDG.graded_mesh

include("FPDE_problem.jl")

Nx = 20
cutoff_t = 1.0
N_outer = 6
mg = 1 + rt / α # mesh grading
Nquad = 15
u(x, t) = refsoln(x, t, α, Nquad)
rx = 4
rt = 3
Mt = 2rt
Mx = 2rx
ppImax = 20

store = Store(α, rt+1, Mt, ppImax)
x_ = collect(range(0, L, length=Nx+1))
x = OffsetVector(x_, 0:Nx)
grid = SpatialGrid(x, rx, ("essential", "essential"), Mx)

U0 = spatial_L2_projection(u0, grid)
t = graded_mesh(max_t, cutoff_t, N_outer, mg)
Nt = length(t) - 1
U = FPDEDG(1.0, f, U0, grid, t, rt, Mt, store)
τ = range(-1, 1, length=ppImax)
pcwise_t, Uerr = L2error(U, u, t, τ, grid, Mx)

jump = evaluate_jumps(U, U0, t, grid, Mx)
Uhat, jumpU = reconstruction(U, U0, store)
pcwise_t, Uhaterr = L2error(Uhat, u, t, τ, grid, Mx)

figure(6)
line1 = semilogy(t[0:Nt-1], abs.(jump[0:Nt-1]), "C0o")
line2 = semilogy(pcwise_t, Uerr, "C1-")
line3 = semilogy(pcwise_t, Uhaterr, "k:")
legend((line1[1], line2[2], line3[1]),
       (L"jump in $U$", L"error in $U$", L"error in $\widehat U$"))
xlabel(L"$t$")
PyPlot.grid(true)
savefig("fig6b.eps")


text:
##
##################################################

"""
 rotations; find values
 Real Givens
 This subroutine computes c and s such that,

 [c -s] * [a, b] = [r,0]; c^2 + s^2 = 1

 and

 r = sqrt(|a|^2 + |b|^2).

"""
@inline function givensrot(a::T, b::T) where {T <: Real}
    u, v = givens(a,b,2,1)
    s = sign(v)
    s*u.c,s*u.s,s*v
end

# @inline function givensrota(a::T,b::T) where {T <: Real}

#     iszero(b) && return (sign(a) * one(T), zero(T), abs(a))
#     iszero(a) && return(zero(T), -sign(b) * one(T), abs(b))

#     r = hypot(a,b)
#     return(a/r,-b/r,r)

#     u, v = givens(a,b,2,1)
#     sign(v)*u.c,u.s,abs(v)
# end

## givens rotation
##################################################
# Compute Givens rotation zeroing b
#
# G1 [ ar + i*ai ] = [ nrm ]
# G1 [    b      ] = [     ]
#
# all variables real (nrm complex)
# returns (copmlex(cr, ci), s) with
# u=complex(cr,ci), v=s; then [u -v; v conj(u)] * [complex(ar, ai), s] has 0
#
# XXX: Could hand write this if needed, here we just use `givens` with a flip
# to get what we want, a real s part, not s part
function givensrot(a::Complex{T},b::Complex{T}) where {T <: Real}
    G, r = givens(b, a, 1, 2)
    G.s, -real(G.c), r
end

function givensrot(a::Complex{T},b::T) where {T <: Real}
    u,v,r = givensrot(a, complex(b, zero(T)))
    u,real(v),r
end

####   Operations on [,[ terms

## The zero_index and stop_index+1 point at "D" matrices
##
## Let a D matrix be one of [1 0; 0 1] or [-1 0; 0 1] (D^2 = I). Then we have this move
## D    --->   D  (we update the rotator)
##   [       [
##
## this is `dflip`
function dflip(a::RealRotator{T}, d=one(T)) where {T}
    return RealRotator(a.c, sign(d)*a.s, a.i)
#    a.s = sign(d)*a.s
end

# get d from rotator which is RR(1,0) or RR(-1, 0)
function getd(a::RealRotator{T}) where {T}
    c, s = vals(a)
    norm(s) <= 4eps(T) || error("a is not a diagonal rotator")
    sign(c)
end

## This is main case
#  Q           D Q
#     D --> D

"""
   D  --> D
U           V
"""
function Dflip(r::ComplexComplexRotator{T}, d::ComplexComplexRotator{T}) where {T}
    !is_diagonal(d) && error("d must be diagonal rotator")

    # D is fixed,
    alpha = d.c
    r.s = r.s * conj(alpha)
end

##   U --> U Da
## D           Da
## (not the reverse!)
function Dflip(d::ComplexRealRotator{T}, r::ComplexRealRotator{T}) where {T}
#    !is_diagonal(d) && error("d must be diagonal rotator")

    alpha = d.c
    c,s = vals(r)
    vals!(r, c*conj(alpha), s)
end

## We have this for left fuse and for deflation
#
#  Di                         Di
#    Qi+1             Si+1     Di+1       Si+1
#        Qi+2    -->    Si+2    Di+2    =    Si+2  * diagm([alpha, I, conj(alpha)])
#           ...           ...     ...          ...
#             Qj            Sj      Dj          Sj
function cascade(Qs, D, alpha, i, j)
    # Q = CR(c,s) -> S = CR(c*conj(alpha), s)

    for k in (i+1):j
        c,s = vals(Qs[k])
        Qs[k] = Rotator(c*conj(alpha), s, idx(Qs[k]))
    end

    D[i] *= alpha
    D[j+1] *= conj(alpha)
end

## Fuse
## fuse combines two rotations, a and b, into one,

## For ComplexRealRotator, the result of a*b will not have a real sign
## we output by rotating by alpha.
## return alpha so a*b = f(ab) * [alpha 0; 0 conj(alpha)]
## for left with have uv -> (u') Di
@inline function fuse(a::ComplexRealRotator{T}, b::ComplexRealRotator{T},::Type{Val{:left}}) where {T}
    #    idx(a) == idx(b) || error("can't fuse")
    u = a.c * b.c - conj(a.s) * b.s
    v = conj(a.c) * b.s + a.s * b.c
    s = norm(v)

    alpha =  iszero(v) ? one(Complex{T}) : conj(v)/s

    c = u * alpha

    Rotator(c,s,idx(a)), conj(alpha)
end

# for right we have uv -> (v') Di
@inline function fuse(a::ComplexRealRotator{T}, b::ComplexRealRotator{T}, ::Type{Val{:right}}) where {T}
#    idx(a) == idx(b) || error("can't fuse")
    u = a.c * b.c - conj(a.s) * b.s
    v = conj(a.c) * b.s + a.s * b.c
    s = norm(v)

    alpha =  iszero(v) ? one(Complex{T}) : conj(v)/s

    c = u * alpha

    Rotator(c, s, idx(b)), conj(alpha)
end

## Fuse for general rotation
## We have two functions as it seems a bit faster
fuse(a::AbstractRotator{T}, b::AbstractRotator{T}, dir, d) where {T} = fuse(a,b,dir)

@inline function fuse(a::AbstractRotator{T}, b::AbstractRotator{T},::Type{Val{:left}}) where {T}
    #    idx(a) == idx(b) || error("can't fuse")

    ac,as = vals(a)
    bc,bs = vals(b)
    u = ac * bc - conj(as) * bs
    v = conj(ac)*bs + as * bc
    return Rotator(u,v,idx(a)), one(T)

end

@inline function fuse(a::AbstractRotator{T}, b::AbstractRotator{T}, ::Type{Val{:right}}) where {T}
    #    idx(a) == idx(b) || error("can't fuse")

    ac,as = vals(a)
    bc,bs = vals(b)

    u = ac * bc - conj(as)*bs
    v = conj(ac)*bs + as*bc
    return Rotator(u, v, idx(b)), one(T)

end

# Turnover: Q1    Q3   | x x x |      Q1
#              Q2    = | x x x | = Q3    Q2  <-- misfit=3 Q1, Q2 shift;
#                      | x x x |
#
# misfit is Val{:right} for <-- (right to left turnover), Val{:left} for -->
#
# This is the key computation once matrices are written as rotators
# We wrote this for complex rotators where sine part may be complex
# so we make use of alpha and beta, which isn't otherwise needed
# could streamline, has small expense, so we break out

@inline function _turnover(Q1::RealRotator{T}, Q2::RealRotator{T}, Q3::RealRotator{T}) where {T}
#    i,j,k = idx(Q1), idx(Q2), idx(Q3)
#    (i == k) || error("Need to have a turnover up down up or down up down: have i=$i, j=$j, k=$k")
#    abs(j-i) == 1 || error("Need to have |i-j| == 1")

    c1,s1 = vals(Q1)
    c2,s2 = vals(Q2)
    c3,s3 = vals(Q3)

    # key is to find U1,U2,U3 with U2'*U1'*U3' * (Q1*Q2*Q3) = I
    # do so by three Givens rotations to make (Q1*Q2*Q3) upper triangular

    # initialize c4 and s4
    a = c1*c2*s3 + s1*c3
    b = s2*s3
    # check norm([a,b]) \approx 1
    c4, s4, temp = givensrot(a,b)#, Val{true})

    # initialize c5 and s5

    a = c1*c3 - s1*c2*s3
    b = temp
    # check norm([a,b]) \approx 1
    c5, s5, alpha = givensrot(a, b)

    # second column
    u = -c1*s3 - s1*c2*c3
    v = c1*c2*c3 - s1*s3
    w = s2 * c3

    a = c4*c5*v - s4*c5*w + s5*u
    b = c4*w + s4*v

    c6, s6, beta = givensrot(a,b)

    (c4, s4, c5, s5, c6, s6)
end

@inline function _turnover(Q1::AbstractRotator{T}, Q2::AbstractRotator{T}, Q3::AbstractRotator{T}) where {T}
#    i,j,k = idx(Q1), idx(Q2), idx(Q3)
#    (i == k) || error("Need to have a turnover up down up or down up down: have i=$i, j=$j, k=$k")
#    abs(j-i) == 1 || error("Need to have |i-j| == 1")

    c1,s1 = vals(Q1)
    c2,s2 = vals(Q2)
    c3,s3 = vals(Q3)

    # key is to find U1,U2,U3 with U2'*U1'*U3' * (Q1*Q2*Q3) = I
    # do so by three Givens rotations to make (Q1*Q2*Q3) upper triangular

    # initialize c4 and s4
    a = conj(c1)*c2*s3 + s1*c3
    b = s2*s3
    # check norm([a,b]) \approx 1
    c4, s4, temp = givensrot(a,b)#, Val{true})

    # initialize c5 and s5

    a = c1*c3 - conj(s1)*c2*s3
    b = temp
    # check norm([a,b]) \approx 1
    c5, s5, alpha = givensrot(a, b)

    alpha = alpha/norm(alpha)
    c5 *= conj(alpha) # make diagonal elements 1
    s5 *= alpha

    # second column
    u = -c1*conj(s3) - conj(s1)*c2*conj(c3)
    v = conj(c1)*c2*conj(c3) - s1*conj(s3)
    w = s2 * conj(c3)

    a = c4*conj(c5)*v - conj(s4)*conj(c5)*w + s5*u
    b = conj(c4)*w + s4*v

    c6, s6, beta = givensrot(a,b)

    beta = beta/norm(beta)
    c6 *= conj(beta) # make diagonal elements 1
    s6 *= beta

    (c4, s4, c5, s5, c6, s6)
end

@inline function turnover(Q1::AbstractRotator{T}, Q2::AbstractRotator{T}, Q3::AbstractRotator{T}, ::Type{Val{:right}}) where {T}

    c4,s4,c5,s5,c6,s6 = _turnover(Q1,Q2,Q3)

    return (Rotator(conj(c5), -s5, idx(Q1)),
            Rotator(conj(c6), -s6, idx(Q2)),
            Rotator(conj(c4), -s4, idx(Q2)) # misfit is right one
            )

end

@inline turnover(Q1::AbstractRotator{T}, Q2::AbstractRotator{T}, Q3::AbstractRotator{T}) where {T} = turnover(Q1, Q2, Q3, Val{:right})

@inline function turnover(Q1::AbstractRotator{T}, Q2::AbstractRotator{T}, Q3::AbstractRotator{T}, ::Type{Val{:left}}) where {T}

    c4,s4,c5,s5,c6,s6 = _turnover(Q1,Q2,Q3)

    return (Rotator(conj(c6), -s6, idx(Q2)), # misfit is left one
            Rotator(conj(c4), -s4, idx(Q2)),
            Rotator(conj(c5), -s5, idx(Q3)))

end

## passthrough
## Pass a rotator through a diagonal matrix with phase shifts
## D U -> U' D'
## Here D[i] = D[1]
## usually call with view(state.d, idx(U):idx(U)+1)
@inline function passthrough(D, U::ComplexRealRotator{T}) where {T}
    i = idx(U)
    alpha, beta = D[i], D[i+1]

    c, s = vals(U)
    u = c * alpha * conj(beta)
    v = s
    U = Rotator(u, v, idx(U))

    D[i], D[i+1] = beta, alpha
    return (D, U)
end

@inline function passthrough(D::ComplexRealRotator{T}, U::ComplexRealRotator{T}) where {T}
    norm(D.s) <= 1e2*eps(T) || error("D not diagonal")
    alpha, ds = vals(D)
    c,s = vals(U)

    U = RealRotator(c * alpha * alpha, s, idx(U))
    D = RealRotator(conj(alpha), ds, idx(D))

    return (D, U)

end