text:
using LinearAlgebra, Statistics

abstract type AbstractStateSpaceFilter end
abstract type AbstractKalmanFilter  <: AbstractStateSpaceFilter end
abstract type NonlinearFilter       <: AbstractStateSpaceFilter end

abstract type LinearKalmanFilter    <: AbstractKalmanFilter end
abstract type NonlinearKalmanFilter <: AbstractKalmanFilter end

# LinearGaussianSSM     LinearKalmanFilter
#
# NonlinearGaussianSSM  NonlinearKalmanFilter
#                       NonlinearFilter
#
# NonlinearSSM          NonlinearFilter

######################################################################
# Kalman filter
######################################################################

struct KalmanFilter <: LinearKalmanFilter end
const  KF = KalmanFilter

## methods

# basic Kalman update, once we have the predicted state, an observation,
# and the observation matrix/covariance.
function update_kalman(pred, y, G, W)

    innovation     = y - G * mean(pred)
    innovation_cov = G * cov(pred) * G' + W
    K = cov(pred) * G' * inv(innovation_cov)

    @show K
    mean_update = mean(pred) + K * innovation
    cov_update  = Symmetric((Matrix(I, size(cov(pred))) - K * G) * cov(pred))

    MvNormal(mean_update, cov_update)

end

function update( m      :: LinearGaussianSSM,
                 pred   :: AbstractMvNormal,
                 y      :: Vector,
                 filter :: KalmanFilter,
                 t      :: Real = 0.0)

    update_kalman(pred, y, m.G(t), m.W(t))

end

function update( m    :: LinearGaussianSSM,
                 pred :: AbstractMvNormal,
                 y    :: Vector,
                 t    :: Real = 0.0)

    update(m, pred, y, KalmanFilter(), t)

end

function update!( m      :: LinearGaussianSSM,
                  fs     :: FilteredState,
                  y      :: Vector;
                  u      :: Vector = zeros(m.nu),
                  filter :: KalmanFilter=KalmanFilter(),
                  t      :: Float64=0.0)

    x_pred = predict(m, fs.state[end], u=u, t=t)
    x_filt = update_kalman(x_pred, y, m.G(t), m.W(t))
    push!(fs.state, x_filt)
    fs.observations = [fs.observations y]

end

######################################################################
# Extended Kalman filter
######################################################################
struct ExtendedKalmanFilter <: NonlinearKalmanFilter end

const EKF = ExtendedKalmanFilter

## methods

function update( m      :: NonlinearGaussianSSM,
                 pred   :: AbstractMvNormal,
                 y      :: Vector,
                 filter :: NonlinearKalmanFilter=EKF(),
                 t      :: Float64=0.0)

    G = observation_matrix(m, pred, t)
    W = m.W(t)
    return update_kalman(pred, y, G, W)

end

function update!( m      :: NonlinearGaussianSSM,
                  fs     :: FilteredState,
                  y      :: Vector;
                  u      :: Vector = zeros(m.nu),
                  filter :: NonlinearKalmanFilter=EKF(),
                  t      :: Float64 = 0.0)

    x_pred = predict(m, fs.state[end], u=u, t=t-1)
    G = observation_matrix(m, x_pred, t)
    W = m.W(t)
    x_filt = update_kalman(x_pred, y, G, W)
    push!(fs.state, x_filt)
    fs.observations = [fs.observations y]

end

######################################################################
# Unscented Kalman filter
######################################################################

struct UnscentedKalmanFilter <: NonlinearKalmanFilter

    α :: Float64
    β :: Float64
    κ :: Float64

    function UnscentedKalmanFilter()

       α = 1e-3
       β = 2.0
       κ = 0.0

       new(α, β, κ)

    end

end

const UKF = UnscentedKalmanFilter

"""
Data structure containing sigma points and their respective weights for the Unscented Kalman Filter

#### Fields
- `χ` : The matrix containing the 2L+1 sigma vectors, where L is the length of the system state vector.
- `wm` : Vector containing the weights for the corresponding sigma vectors required to reconsruct the predicted mean.
- `wc` : Vector containing the weights for the corresponding sigma vectors required to reconsruct the predicted covariance.
"""
struct SigmaPoints{T<:Real}
    χ  :: Matrix{T}
    wm :: Vector{T}
    wc :: Vector{T}
end

"""
# calcSigmaPoints
Function to calculate the sigma points for the Unscented Transform

`calcSigmaPoints(state, α, β, κ)`
#### Parameters
- `state` : AbstractMvNormal type representing the current state estimate (mean vector with covariance matrix).
- `α` : Determines the spread of the sigma points around the mean of the system state. It is usually set to a small positive value. The default is set to 1e-3.
- `β` : Is used to incorporate prior knowledge of the distribution of the state. For Gaussian distributions a value of 2.0 is optimal. The default is set to 2.0.
- `κ` : Is a secondary scaling parameter which also determines the spread of sigma points around the mean of the system. This parameter allows for additional 'fine tuning'. It is usually set to a small value. The default is set to 0.0

`calcSigmaPoints(state, filter)`
#### Parameters
- `state` : AbstractMvNormal type representing the current state estimate (mean vector with covariance matrix).
- `filter` : UKFParameters type containing the α, β and κ parameters for the Unscented Kalman Filter.
"""
function calcSigmaPoints( state :: AbstractMvNormal,
                          α     :: T,
                          β     :: T,
                          κ     :: T) where T

    x  = mean(state)
    p  = cov(state)
    L  = length(x)
    χ  = zeros(L, 2L+1)
    wm = zeros(2L+1)
    wc = zeros(2L+1)
    λ  = α^2 * (L+κ) - L
    χ[:,1] = x
    wm[1] = λ/(L+λ)
    wc[1] = wm[1] + (1 - α^2 + β)
    wm[2:end] = 1/(2*(L+λ))
    wc[2:end] = wm[2:end]
    sp = sqrt(L+λ)*full(chol(p))
    for i = 2:L+1
        χ[:,i] = x + sp[:,i-1]
        χ[:,i+L] = x - sp[:,i-1]
    end
    return SigmaPoints(χ, wm, wc)
end

function calcSigmaPoints( state :: AbstractMvNormal, filter :: UKF)

    calcSigmaPoints(state, filter.α, filter.β, filter.κ)

end

"""
# timeUpdate
Function to calculate the predicted mean and the predicted covariance given a UKF model and the sigma points with their corresponding weights

`timeUpdate(m, sp)`
#### Parameters
- `m` : AbstractGaussianSSM type containing the parameters of the Unscented Kalman Filter model.
- `sp` : SigmaPoints type containing the matrix of sigma vectors and their corresponding weights
"""
function timeUpdate(m :: AbstractGaussianSSM, sp :: SigmaPoints, t::Real=0.0)

    L, M  = size(sp.χ)
    χ_x = zeros(L, M)
    x_pred = zeros(L)
    p_pred = zeros(L, L)
    for i in 1:2L+1
        χ_x[:,i] = m.f(sp.χ[:,i])
        x_pred += sp.wm[i] * χ_x[:,i]
    end
    for i in 1:2L+1
        p_pred += sp.wc[i] * (χ_x[:,i] - x_pred)*(χ_x[:,i] - x_pred)'
    end
    p_pred += m.V(t)
    return MvNormal(x_pred, p_pred), SigmaPoints(χ_x, sp.wm, sp.wc)

end

function predict(m::AbstractGaussianSSM, x::AbstractMvNormal, filter::UKF)

    sigPoints = calcSigmaPoints(x, filter)
    pred_state, new_sigPoints = timeUpdate(m, sigPoints)
    return pred_state, new_sigPoints

end

function observe(m::AbstractGaussianSSM, x::AbstractMvNormal, sp::SigmaPoints, y, t::Real=0.0)

    obsLength = length(y)
    L, M = size(sp.χ)
    y_trans = zeros(obsLength, M)
    y_pred = zeros(obsLength)
    for i in 1:2L+1
        y_trans[:,i] = m.g(sp.χ[:,i])
        y_pred += sp.wm[i] * y_trans[:,i]
    end

    P_xy = zeros(L, obsLength)
    P_yy = zeros(obsLength, obsLength)
    for i in 1:2L+1
        resy = (y_trans[:,i] - y_pred)
        P_xy += sp.wc[i] * (sp.χ[:,i] - mean(x)) * resy'
        P_yy += sp.wc[i] * resy * resy'
    end
    P_yy += m.W(t)
    return MvNormal(y_pred, P_yy), P_xy

end

function innovate( m     :: AbstractGaussianSSM,
                   x     :: AbstractMvNormal,
                   yPred :: AbstractMvNormal,
                   P_xy  :: Matrix,
                   sp    :: SigmaPoints,
                   y     :: Vector)

    kalmanGain = P_xy * inv(cov(yPred))
    new_x = mean(x) + kalmanGain * (y - mean(yPred))
    new_cov = cov(x) - kalmanGain * cov(yPred) * kalmanGain'
    return MvNormal(new_x, new_cov)
end

function update(m::AbstractGaussianSSM, x::AbstractMvNormal, sp::SigmaPoints, y::Vector)
    yPred, P_xy = observe(m, x, sp, y)
    return innovate(m, x, yPred, P_xy, sp, y)
end

######################################################################
# Ensemble Kalman filter
######################################################################

struct EnsembleKalmanFilter <: NonlinearFilter

    nparticles :: Int

    EnsembleKalmanFilter() = new(100)

end

const EnKF = EnsembleKalmanFilter

function predict( m        :: AbstractGaussianSSM,
                  ensemble :: Matrix,
                  filter   :: EnKF=EnKF();
                  u        :: Vector=zeros(m.nu),
                  t        :: Real=0.0)

    ensemble_new = similar(ensemble)
    CI = control_input(m, u, t)
    for i in 1:size(ensemble, 2)
        F = process_matrix(m, ensemble[:, i], t)
        ensemble_new[:, i] = F * ensemble[:, i] .+ CI
    end
    return ensemble_new

end

function update( m        :: AbstractGaussianSSM,
                 ensemble :: Matrix,
                 y        :: Vector,
                 kf       :: EnKF = EnsembleKalmanFilter(),
                 t        :: Real=0.0)

    P = cov(ensemble')
    ensemble_updated = similar(ensemble)
    for i in 1:kf.nparticles
        G = observation_matrix(m, ensemble[:, i], t)
        W = observation_matrix(m, ensemble[:, i], t)
        innovation = y - G * ensemble[:, i]
        innovation_cov = G * P * G' + W
        K = P * G' * inv(innovation_cov)
        ensemble_updated[:, i] = ensemble[:, i] + K * innovation
        cov_update = Symmetric((Matrix(I,size(P)) - K * G) * P)
    end
    ensemble_updated
end

import Base.filter

function filter( m      :: AbstractGaussianSSM,
                 y      :: Array{T},
                 x0     :: AbstractMvNormal,
                 kf     :: UnscentedKalmanFilter ) where T

    x_filtered = Array(AbstractMvNormal, size(y, 2))
    loglik = 0.0
    x_pred, sigma_points = predict(m, x0, kf)
    x_filtered[1] = update(m, x_pred, sigma_points, y[:, 1])
    for i in 2:size(y, 2)
        x_pred, sigma_points = predict(m, x_filtered[i-1], kf)
        # Check for missing values in observation
        if any(isnan(y[:, i]))
            x_filtered[i] = x_pred
        else
            x_filtered[i] = update(m, x_pred, sigma_points, y[:, i])
            loglik += logpdf(observe(m, x_filtered[i],
                calcSigmaPoints(x_filtered[i], filter.α, filter.β,
                                filter.κ), y[:, 1])[1], y[:, 1])
        end
        loglik += logpdf(x_pred, mean(x_filtered[i]))
    end
    FilteredState(y, x_filtered, loglik, false)
end

function filter( m              :: AbstractGaussianSSM,
                 y              :: Array{T},
                 x0             :: AbstractMvNormal,
                 kf             :: EnsembleKalmanFilter = EnsembleKalmanFilter();
                 u              :: Matrix{T} = zeros(m.nu, size(y, 2)),
                 times          :: Vector{T} = zeros(size(y, 2)),
                 return_ensemble = false) where T

    nt = size(y, 2)
    ensemble = rand(x0, kf.nparticles)
    ensemble = predict(m, ensemble, kf, u=u[:,1], t=times[1])
    loglik = 0.0
    if return_ensemble
        x_filtered = zeros(m.nx, kf.nparticles, nt)
        x_filtered[:, :, 1] = update(m, ensemble, y[:, 1], kf, times[1])
        for i in 2:nt
            ensemble = predict(m, x_filtered[:, :, i-1], kf, u=u[:, i], t=times[i])
            x_filtered[:, :, i] = update(m, ensemble, y[:, i], kf, times[i])
        end
        return x_filtered
    else
        x_filtered = Array(AbstractMvNormal, size(y, 2))
        ensemble = update(m, ensemble, y[:, 1], kf, times[1])
        x_filtered[1] = MvNormal(vec(mean(ensemble, 2)), cov(ensemble'))
        for j in 1:kf.nparticles
            loglik += logpdf(x_filtered[1], ensemble[:, j])
        end
        if !any(isnan(y[:, 1]))
            for j in 1:kf.nparticles
                loglik += logpdf(observe(m, x_filtered[1], times[1]), y[:,1])
            end
        end
        for i in 2:nt
            ensemble = predict(m, ensemble, kf, u=u[:, i], t=times[i])
            x_filtered[i] = MvNormal(vec(mean(ensemble, 2)), cov(ensemble'))
            for j in 1:kf.nparticles
                loglik += logpdf(x_filtered[i], ensemble[:, j])
            end
            if ! any(isnan(y[:, i]))
                for j in 1:kf.nparticles
                    loglik += logpdf(observe(m, x_filtered[1], times[1]), y[:,i])
                end
            end
        end
        FilteredState(y, x_filtered, u, times, loglik, false)
    end
end


text:
module utils

include("metric.jl")

end

ext: jl


text:
function spec_clustering(data,k)

    w = computing_similarity(data)
    d = diagm(vec(sum(w,1)))
    l = d-w
    temp = eig(l)
    temp = temp[2]
    e_map = temp[:,1:(k-1)] #seems that the largest k eigenvalue works?
    #e_map = temp[:,(end-k+1):end]
    model = Kmeans(k=k)
    train!(model,e_map)
    predict!(model)
    return model
end

function computing_similarity(data)
    n_sample = size(data,1)
    w = zeros(n_sample,n_sample)
    for i = 1:n_sample
        for j = 1:n_sample
            w[i,j] = count_sim(data[i,:],data[j,:])
        end
    end
    return w
end

function count_sim(x::Vector,y::Vector;
                   types="Gaussian",
                    gamma = 1)
    if types == "Gaussian"
        return exp(-gamma*norm(x-y)^2)
    end
end

function test_spec_cluster()

    X, y = make_blo()
    clu = length(unique(y))
    model = spec_clustering(X,clu)
    predictions = model.clusters
    plot_in_2d(model)

end


text:
using BenchmarkTools
const SUITE = BenchmarkGroup()
SUITE["post_evolve"] = include("post_evolve.jl")


text:
using IncrementalInference
using Graphs
# using KernelDensityEstimate, Gadfly # for vstack

using Base: Test

fg = emptyFactorGraph()

N=100

doors = [-100.0;0.0;100.0;300.0]'
cov = [3.0]

v1 = addNode!(fg,:x1, ContinuousScalar, N=N)
f1  = addFactor!(fg, [:x1;], Obsv2(doors, cov', [1.0]))

tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addNode!(fg,:x2, ContinuousScalar, N=N)
addFactor!(fg,[:x1, :x2],Odo([50.0]',[2.0]',[1.0]))

v3=addNode!(fg, :x3, ContinuousScalar, N=N) # 4.0*randn(1,N)+getVal(v2)+50.0
addFactor!(fg,[:x2,:x3],Odo([50.0]',[4.0]',[1.0]))

l1=addNode!(fg, :l1, ContinuousScalar, N=N) # 0.5*randn(1,N)+getVal(v3)+64.0
addFactor!(fg, [:x1,:l1], Ranged([64.0],[0.5],[1.0]))
addFactor!(fg, [:x2,:l1], Ranged([16.0],[0.5],[1.0]))

l2=addNode!(fg, :l2, ContinuousScalar, N=N) # 0.5*randn(1,N)+getVal(v3)+64.0
addFactor!(fg, [:x3,:l2], Ranged([16.0],[0.5],[1.0]))

# for thesis
# v4=addNode!(fg,:x4,4.0*randn(1,N)+getVal(v2)+50.0, N=N)
# addFactor!(fg,[:x4,:l2],Odo([50.0]',[4.0]',[1.0]))

# Graphs.plot(fg.g)
# writeGraphPdf(fg);
# run(`evince fg.pdf`)

# p = IncrementalInference.getEliminationOrder(fg, ordering=:qr)
p = [7,10,1,3,5,12];
p = [7,10,1,3,5];

println()
fge = deepcopy(fg)
println("Building Bayes net...")
buildBayesNet!(fge, p)

tree = emptyBayesTree()
buildTree!(tree, fge, p)

# println("Bayes Net")
# sleep(0.1)
#fid = open("bn.dot","w+")
#write(fid,to_dot(fge.bn))
#close(fid)

@test num_vertices(tree.bt) == 3

# Michael reference -- x2->x1, x2->x3, x2->x4, x2->l1, x4->x3, l1->x3, l1->x4
if false
  println("Bayes Tree")
  # Graphs.plot(tree.bt)
  fid = open("bt.dot","w+")
  write(fid,Graphs.to_dot(tree.bt))
  close(fid)
  run(`dot bt.dot -Tpdf -o bt.pdf`)
end

# GraphViz.Graph(to_dot(tree.bt))
#Michael reference 3sig -- x2l1x4x3    x1|x2

import IncrementalInference: buildCliquePotentials

println("Find potential functions for each clique")
cliq = tree.cliques[1] # start at the root
buildCliquePotentials(fg, tree, cliq); # fg does not have the marginals as fge does

# now update all factor graph vertices used for this tree
for v in Graphs.vertices(fg.g)
  IncrementalInference.dlapi.updatevertex!(fg, v)
end

# TODO -- add testing to ensure this is the correct tree!

# run(`evince bt.pdf`)


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
using Test

DDEProblemLibrary.importddeproblems()

const alg = MethodOfSteps(RK4(); constrained = false)
const prob = DDEProblemLibrary.prob_dde_constant_1delay_scalar

# reference solution with delays specified
@testset "reference" begin
  sol = solve(prob, alg)

  @test sol.errors[:l∞] < 5.6e-5
  @test sol.errors[:final] < 1.8e-6
  @test sol.errors[:l2] < 2.0e-5
end

# problem without delays specified
const prob_wo = remake(prob; constant_lags = nothing)

# solutions with residual control
@testset "residual control" begin
  sol = solve(prob_wo, alg)

  @test sol.errors[:l∞] < 1.4e-4
  @test sol.errors[:final] < 3.5e-6
  @test sol.errors[:l2] < 6.5e-5

  sol = solve(prob_wo, alg; abstol = 1e-9, reltol = 1e-6)

  @test sol.errors[:l∞] < 5.0e-8
  @test sol.errors[:final] < 3.4e-9
  @test sol.errors[:l2] < 2.2e-8

  sol = solve(prob_wo, alg; abstol = 1e-13, reltol = 1e-13)

  # relaxed tests to prevent floating point issues
  @test sol.errors[:l∞] < 5e-11
  @test sol.errors[:final] < 4.5e-12
  @test sol.errors[:l2] < 7.7e-11 # 7.7e-12
end

######## Now show that non-residual control is worse
# solutions without residual control
@testset "non-residual control" begin
  sol = solve(prob_wo, MethodOfSteps(OwrenZen5(); constrained = false))

  @test sol.errors[:l∞] > 1.1e-1
  @test sol.errors[:final] > 1.2e-3
  @test sol.errors[:l2] > 4.5e-2

  sol = solve(prob_wo, MethodOfSteps(OwrenZen5(); constrained = false); abstol=1e-13, reltol=1e-13)

  @test sol.errors[:l∞] > 1.1e-1
  @test sol.errors[:final] > 1.2e-3
  @test sol.errors[:l2] > 5.5e-2
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

@testset "Repeat" begin
    function hex3RGB()
        clusterelements = SVector((0,0),(-1,0),(-1,1))
        colors = [colorant"red",colorant"green",colorant"blue"]
        names = ["R","G","B"]
        eltlattice = Repeat.HexBasis1()
        clusterbasis = Repeat.LatticeBasis(( -1,2),(2,-1))
        lattice = Repeat.LatticeCluster(clusterbasis,eltlattice,clusterelements)
        properties =  DataFrame(Color = colors, Name = names)
        return Repeat.ClusterWithProperties(lattice,properties)
    end

    #spherepoint tests
    @test isapprox(Repeat.Lenslets.spherepoint(1,π/2,0.0), [0.0,1.0,0.0])
    @test isapprox(Repeat.Lenslets.spherepoint(1,0.0,π/2), [1.0,0.0,0.0])
    @test isapprox(Repeat.Lenslets.spherepoint(1,0,0.0), [0.0,0.0,1.0])
    @test isapprox(Repeat.Lenslets.spherepoint(1,0.0,π/4), [sqrt(2)/2,0.0,sqrt(2)/2])

    """ Create a LatticeCluser with three elements at (0,0),(-1,0),(-1,1) coordinates in the HexBasis1 lattice"""
    function hex3cluster()
        clusterelts = SVector((0,0),(-1,0),(-1,1))
        eltlattice = Repeat.HexBasis1()
        clusterbasis = Repeat.LatticeBasis(( -1,2),(2,-1))
        return Repeat.LatticeCluster(clusterbasis,eltlattice,clusterelts)
    end

    @test [-1 2;2 -1] == Repeat.basismatrix(Repeat.clusterbasis(hex3RGB()))

    function basistest(a::Repeat.AbstractLatticeCluster)
        return Repeat.clusterbasis(a)
    end

    @test basistest(hex3cluster()) == basistest(hex3RGB())

    #LatticeCluster testset
    cluster = Repeat.Lenslets.hex9()

    for iter in 1:100
        (i,j) = rand.((1:1000,1:1000))
        coords,tileindex = Repeat.cluster_coordinates_from_tile_coordinates(cluster,i,j)
        reconstructed = Repeat.tilecoordinates(cluster,coords...,tileindex)
        @test all((i,j) .== reconstructed)
    end

    #verify that the 0,0 cluster is correct
    for (index,element) in pairs(Repeat.clusterelements(cluster))
        coords, tileindex = Repeat.cluster_coordinates_from_tile_coordinates(cluster, element...)
        @test all(coords .== 0)
        @test tileindex == index
    end
end


text:
using LinearAlgebra
import Base: size, permutedims, permutedims!, getindex, setindex!, ndims, show, convert
import Base: iterate, length, similar, adjoint, eltype, +, -, *, /, ^, BroadcastStyle
import Strided: UnsafeStridedView

export FermiMDArray, FermiMDrand, FermiMDzeros, FermiSparse, diagonalize

abstract type AbstractFermiArray{T,N} <: AbstractArray{T,N} end
"""

    FermiMDArray{T,N}

Fermi array object held entirely in memory. Thin wrap around a standard Julia array, representing a dense array of type T and rank N.

# Struct tree

**FermiMDArray** <: AbstractArray
"""
struct FermiMDArray{T,N} <: AbstractFermiArray{T,N}
    data::Array{T,N}
end

struct Fermi4SymArray{T} <: AbstractFermiArray{T,1}
    data::Array{T,1}
end

struct FermiSparse{Td,Ti,N} <: AbstractFermiArray{Td,N}
    indexes::Vector{NTuple{N, Ti}}
    data::Vector{Td}
end

function ndims(A::FermiSparse{Td,Ti,N}) where {Td, Ti, N}
    return N
end

show(io::IO, ::MIME"text/plain", A::FermiSparse{Td, Ti, N}) where {Td, Ti, N} = show(io,A)

function show(io::IO, A::FermiSparse{Td, Ti, N}) where {Td, Ti, N}
    println(io, "FermiSparse{$Td, $Ti, $N}:")
    print(io, A.data)
    #if length(A.data) > 10
    #    for (x,idx) in zip(A.data[1:5], A.indexes[1:5])
    #        idx_str = "("
    #        for i in idx
    #            idx_str *= "$i,"
    #        end
    #        idx_str = idx_str[1:end-1]*") => $x"
    #        println(io, idx_str)
    #    end
    #    println(io," ⋮")
    #    for (x,idx) in zip(A.data[end-5:end], A.indexes[end-5:end])
    #        idx_str = "("
    #        for i in idx
    #            idx_str *= "$i,"
    #        end
    #        idx_str = idx_str[1:end-1]*") => $x"
    #        println(io, idx_str)
    #    end
    #else
    #    for (x,idx) in zip(A.data, A.indexes)
    #        idx_str = "("
    #        for i in idx
    #            idx_str *= "$i,"
    #        end
    #        idx_str = idx_str[1:end-1]*") => $x"
    #        println(io, idx_str)
    #    end
    #end
end

function index2(i::Signed, j::Signed)::Signed
    if i < j
        return (j * (j + 1)) >> 1 + i
    else
        return (i * (i + 1)) >> 1 + j
    end
end

# This function is not useful currently, but it will be for direct computations
#function index4(i::Signed , j::Signed, k::Signed, l::Signed)::Signed
#    return index2(index2(i,j), index2(k,l))
#end
#
#function getindex(A::Fermi4SymArray, I::Vararg{Signed,4})
#    idx = index4((I .- 1)...) + 1
#    return A.data[idx]
#end
#
#function getindex(A::Fermi4SymArray, i::Signed)
#    return A.data[i]
#end

# Creates a FermiMDArray from a native Julia Array
function FermiMDArray(A::AbstractArray)
    return FermiMDArray(Array(A))
end

# If trying to create a FermiMDArray from another FermiMDArray, do nothing
function FermiMDArray(A::FermiMDArray{T,N}) where {T,N}
    return A
end

# If trying to create a FermiMDArray from a number, returns the number
function FermiMDArray(num::Number)
    return num
end

"""
    FermiMDzeros(x...)

Create a Julia Array as `zeros(x...)` wrapped in a FermiMDArray.
"""
function FermiMDzeros(x...)
    data = zeros(x...)
    return FermiMDArray(data)
end

"""
    FermiMDrand(x...)

Create a Julia Array as `rand(x...)` wrapped in a FermiMDArray.
"""
function FermiMDrand(x...)
    data = rand(x...)
    return FermiMDArray(data)
end

function similar(A::FermiMDArray)
    data = similar(A.data)
    return FermiMDArray(data)
end

function similar(A::FermiMDArray, dims::Dims)
    data = similar(A.data, dims)
    return FermiMDArray(data)
end

# Basic methods for AbstractArrays in Julia
function size(A::T, i...) where T <: AbstractFermiArray
    return size(A.data, i...)
end

function getindex(A::AbstractFermiArray, I...)
    return FermiMDArray(A.data[I...])
end

function setindex!(A::AbstractFermiArray, val, I...)
    A.data[I...] = val
end

function ndims(A::FermiMDArray)
    return ndims(A.data)
end

function length(A::AbstractFermiArray)
    return length(A.data)
end

function eltype(A::AbstractFermiArray)
    return eltype(A.data)
end

# Special methods for LinearAlgebra
function permutedims(A::FermiMDArray,tup)
    FermiMDArray(permutedims(A.data,tup))
end

function permutedims!(A::FermiMDArray,tup)
    A.data .= permutedims(A.data,tup)
end

"""
    Fermi.diagonalize(A::FermiMDArray; sortby=x->x)

Diagonalize a NxN FermiMDArray, returning eigenvalues and eigenvectors.

# Example
```
julia> A = FermiMDArray([3 6; 6 4])
julia> ϵ,ν = diagonalize(A)
julia> ϵ
Fermi Memory-held Dense Array - 2-element Array{Float64,1}:
 -2.5207972893961488
  9.520797289396146
```
"""
function diagonalize(A::FermiMDArray; sortby=x->x, hermitian=true)
    if hermitian
        vals, vecs = LinearAlgebra.eigen(Hermitian(A.data), sortby=sortby)
        return FermiMDArray(vals), FermiMDArray(vecs)
    else
        vals, vecs = LinearAlgebra.eigen(A.data, sortby=sortby)
        return FermiMDArray(vals), FermiMDArray(vecs)
    end
end

function LinearAlgebra.Hermitian(A::FermiMDArray)
    return FermiMDArray(LinearAlgebra.Hermitian(A.data))
end

function LinearAlgebra.diag(A::FermiMDArray)
    return FermiMDArray(LinearAlgebra.diag(A.data))
end

function LinearAlgebra.factorize(A::FermiMDArray)
    return factorize(A.data)
end

function adjoint(A::FermiMDArray)
    return FermiMDArray(adjoint(A.data))
end

# Basic mathematical methods
function Base.:+(A::FermiMDArray, B::FermiMDArray)
    return FermiMDArray(A.data + B.data)
end

function Base.:+(A::FermiMDArray, B::AbstractArray)
    return FermiMDArray(A.data + B)
end

function Base.:+(A::AbstractArray, B::FermiMDArray)
    return FermiMDArray(A + B.data)
end

function Base.:-(A::FermiMDArray, B::FermiMDArray)
    return FermiMDArray(A.data - B.data)
end

function Base.:-(A::FermiMDArray, B::AbstractArray)
    return FermiMDArray(A.data - B)
end

function Base.:-(A::AbstractArray, B::FermiMDArray)
    return FermiMDArray(A - B.data)
end

function Base.:*(A::FermiMDArray, B::Number)
    return FermiMDArray(B*A.data)
end

function Base.:*(A::Number, B::FermiMDArray)
    return FermiMDArray(A*B.data)
end

function Base.:*(A::FermiMDArray, B::FermiMDArray)
    return FermiMDArray(A.data*B.data)
end

function Base.:^(A::FermiMDArray, B::Integer)
    return FermiMDArray(A.data^B)
end

function Base.:^(A::FermiMDArray, B::AbstractFloat)
    return FermiMDArray(A.data^B)
end

function UnsafeStridedView(A::FermiMDArray)
    UnsafeStridedView(A.data)
end

# Methods to allow Broadcasting
function Base.BroadcastStyle(::Type{<:FermiMDArray})
    Broadcast.ArrayStyle{FermiMDArray}()
end

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{FermiMDArray}}, ::Type{ElType}) where ElType
    FermiMDArray(similar(Array{ElType}, axes(bc)))
end

function Base.convert(::Type{T}, A::FermiMDArray) where T<:FermiMDArray{Float64}
    newdata = Base.convert(Array{Float64}, A.data)
    return FermiMDArray(newdata)
end


text:
using QuBase
using Base.Test

include("multest.jl")
include("tensortest.jl")
include("labelbasistest.jl")
include("operatortest.jl")


text:
function loadNonStopWords(fileName::String)::Tuple{Vector{String},Dict{String,Int}}
    file = open(fileName)
    nonStopWords = Vector{String}()
    nonStopWordIDs = Dict{String,Int}()
    for phrase in readlines(file)
        if phrase != ""
            push!(nonStopWords, phrase)
            nonStopWordIDs[phrase] = length(nonStopWords)
        end
    end
    close(file)
    (nonStopWords, nonStopWordIDs)
end

function loadDualAspectPairFreqs(fileName::String)::Vector{Tuple{Tuple{Int,Int},Int}}
    file = open(fileName)
    dualAspectPairFreqs = Vector{Tuple{Tuple{Int,Int},Int}}()
    for line in readlines(file)
        if line == ""
            continue
        end
        fields = split(line)
        @assert length(fields) == 3
        i = parse(Int, fields[1])
        j = parse(Int, fields[2])
        freq = parse(Int, fields[3])
        push!(dualAspectPairFreqs, ((i, j), freq))
    end
    close(file)
    dualAspectPairFreqs
end

function loadConns(fileName::String, n::Int)::Vector{Dict{Int,Float64}}
    sparseConns = [Dict{Int,Float64}() for i in 1:n]
    file = open(fileName)
    for line in readlines(file)
        if line == ""
            continue
        end
        fields = split(line)
        @assert length(fields) == 3
        i = parse(Int, fields[1])
        j = parse(Int, fields[2])
        conn = parse(Float64, fields[3])
        sparseConns[i][j] = conn
    end
    close(file)
    sparseConns
end

function deltaModularity(sparseConns::Vector{Dict{Int,Float64}}, sumConn::Vector{Float64}, totalConn::Float64, u::Int, oldCu::Set{Int}, newCu::Set{Int})
    result = 0.0
    for v in oldCu
        if v != u
            wuv = sumConn[u] * sumConn[v] / totalConn
            result -= get(sparseConns[u], v, 0.0) - wuv
        end
    end
    for v in newCu
        wuv = sumConn[u] * sumConn[v] / totalConn
        result += get(sparseConns[u], v, 0.0) - wuv
    end
    result
end

function sumDualAspectConn(sparseConns::Vector{Dict{Int,Float64}})::Vector{Float64}
    numPoints = length(sparseConns)
    sumConn = Vector{Float64}(undef, numPoints)
    for i in 1:numPoints
        sumI = 0.0
        for (j, conn) in sparseConns[i]
            sumI += conn
        end
        sumConn[i] = sumI
    end
    sumConn
end

function findDualAspectTopics(sparseConns::Vector{Dict{Int,Float64}}, maxIter::Int)::Tuple{Vector{Int},Vector{Set{Int}}}
    sumConn = sumDualAspectConn(sparseConns)
    totalConn = sum(sumConn)
    numPoints = length(sumConn)
    clusters = collect(Set{Int}(i) for i in 1:numPoints)
    clusterIDs = collect(1:numPoints)
    numClusters = numPoints

    numCPUs = length(Sys.cpu_info())
    for iter in 1:maxIter
        ###########################################################################################
        # step 1: randomly pick a destination with positive gain for each point
        jobs = Channel{Tuple{Int,Int}}(numCPUs)
        Segments = Vector{Tuple{Int,Int}}
        WeightedDests = Vector{Tuple{Float64,Int}}
        jobOutputs = Channel{Tuple{Segments,Vector{WeightedDests}}}(numCPUs)

        function makeJobs(batchSize::Int)
            for i in 1:batchSize:numPoints
                put!(jobs, (i, min(i + batchSize - 1, numPoints)))
            end
        end

        function runJob(iCPU::Int)
            mySegments = Segments()
            myWeightedDests = Vector{WeightedDests}()
            for job in jobs
                push!(mySegments, job)
                jobWeightedDests = WeightedDests(undef, job[2] - job[1] + 1)
                # println("cpu $iCPU: computing weighted dests for $job")
                for u in job[1]:job[2]
                    # println("cpu $iCPU: computing for $u")
                    oldCuID = clusterIDs[u]
                    oldCu = clusters[oldCuID]
                    neighborClusterIDs = Set{Int}()
                    for (neighbor, conn) in sparseConns[u]
                        push!(neighborClusterIDs, clusterIDs[neighbor])
                    end
                    neighborClusterIDs = [clusterID for clusterID in neighborClusterIDs]
                    numNeighborClusters = length(neighborClusterIDs)
                    attractions = fill(0.0, numNeighborClusters)

                    for (i, newCuID) in enumerate(neighborClusterIDs)
                        if newCuID != oldCuID
                            newCu = clusters[newCuID]
                            attractions[i] = max(0.0, deltaModularity(sparseConns, sumConn, totalConn, u, oldCu, newCu))
                        end
                    end
                    # println("cpu $iCPU: attractions of $u is $attractions")
                    sumAttractions = sum(attractions)
                    thres = rand() * sumAttractions
                    partialSumAttr = 0.0
                    weightedDest = (0.0, oldCuID)
                    for dest in 1:numNeighborClusters
                        attr = attractions[dest]
                        partialSumAttr += attr
                        if attr > 0.0 && partialSumAttr >= thres
                            weightedDest = (attr, neighborClusterIDs[dest])
                            break
                        end
                    end
                    jobWeightedDests[u - job[1] + 1] = weightedDest
                    # if u % 10000 == 0
                    #     println("cpu $iCPU: $u of $numPoints have been processed")
                    # end
                end # u
                push!(myWeightedDests, jobWeightedDests)
                # println("cpu $iCPU: finish computing weighted dests for $job")
            end # job
            put!(jobOutputs, (mySegments, myWeightedDests))
        end # runJob

        bind(jobs, @async makeJobs(100))
        for iCPU in 1:numCPUs
            Threads.@spawn runJob(iCPU)
        end
        weightedDests = WeightedDests(undef, numPoints)
        for iCPU in 1:numCPUs
            segmentedWeightedDests = take!(jobOutputs)
            for (segID, segment) in enumerate(segmentedWeightedDests[1])
                weightedDests[segment[1]:segment[2]] = segmentedWeightedDests[2][segID]
            end
        end

        ###########################################################################################
        # step 2: merge point-movement requests
        sortedPoints = sortperm(weightedDests, rev=true)
        clusterChanged = fill(false, numClusters)
        mergedRequests = Vector{Tuple{Int,Int}}()
        mergedGain = 0.0
        for u in sortedPoints
            gain = weightedDests[u][1]
            if gain <= 0.0
                break
            end
            oldCuID = clusterIDs[u]
            newCuID = weightedDests[u][2]
            if clusterChanged[oldCuID] || clusterChanged[newCuID]
                continue
            end
            push!(mergedRequests, (u, newCuID))
            mergedGain += gain
            clusterChanged[oldCuID] = true
            clusterChanged[newCuID] = true
        end

        numMerged = length(mergedRequests)
        if numMerged == 0
            break
        end

        ###########################################################################################
        # step 3: move points according to the merged requests
        for request in mergedRequests
            u = request[1]
            newCuID = request[2]
            oldCuID = clusterIDs[u]
            delete!(clusters[oldCuID], u)
            push!(clusters[newCuID], u)
            clusterIDs[u] = newCuID
        end

        ###########################################################################################
        # step 4: remove empty clusters
        iLeft = 1
        iRight = numClusters
        while length(clusters[iRight]) == 0 && iRight > iLeft
            iRight -= 1
        end
        while iLeft < iRight
            if length(clusters[iLeft]) == 0
                temp = clusters[iLeft]
                clusters[iLeft] = clusters[iRight]
                clusters[iRight] = temp
                iRight -= 1
                while length(clusters[iRight]) == 0 && iRight > iLeft
                    iRight -= 1
                end
                for u in clusters[iLeft]
                    clusterIDs[u] = iLeft
                end
            end
            iLeft += 1
        end
        if length(clusters[iLeft]) == 0
            numClusters = iLeft - 1
        else
            numClusters = iLeft
        end

        clusters = clusters[1:numClusters]
        #println("Topic Iter $iter: moved $numMerged of $numPoints points to gain $mergedGain, $numClusters clusters remain")
    end # iter

    (clusterIDs, clusters)
end

function loadAspects(fileName::String, nonStopWordIDs::Dict{String,Int})::Tuple{Vector{Set{Int}},Vector{Int}}
    file = open(fileName)
    aspects = Vector{Set{Int}}()
    numPhrases = length(nonStopWordIDs)
    aspectIDs = Vector{Int}(undef, numPhrases)
    for line in readlines(file)
        fields = split(line)
        if length(fields) != 3
            continue
        end
        aspectID = parse(Int, fields[1])
        card = parse(Int, fields[2])
        phrase = fields[3]
        phraseID = nonStopWordIDs[phrase]
        if aspectID > length(aspects)
            @assert aspectID == length(aspects) + 1
            push!(aspects, Set{Int}())
        end
        push!(aspects[aspectID], phraseID)
        aspectIDs[phraseID] = aspectID
    end
    close(file)
    (aspects, aspectIDs)
end

function saveTopics(fileName::String, topics::Vector{Set{Int}}, nonStopWords::Vector{String}, dualAspectPairFreqs::Vector{Tuple{Tuple{Int,Int},Int}})
    file = open(fileName, "w")
    for (topicID, topic) in enumerate(topics)
        for pairID in topic
            pairFreq = dualAspectPairFreqs[pairID]
            pair = pairFreq[1]
            freq = pairFreq[2]
            phrase1 = nonStopWords[pair[1]]
            phrase2 = nonStopWords[pair[2]]
            println(file, "$topicID $freq $phrase1 $phrase2")
        end
    end
    close(file)
end

function divideGroups(dualAspectPairFreqs::Vector{Tuple{Tuple{Int,Int},Int}}, regions::Vector{Vector{Int}}, aspectIDs::Vector{Int}, maxPairsInGroup::Int)::Vector{Vector{Int}}
    groups = Vector{Vector{Int}}()
    for region in regions
        if length(region) < maxPairsInGroup
            push!(groups, region)
            continue
        end

        subRegions = Dict{Tuple{Int,Int},Vector{Int}}()
        for pairID in region
            pairFreq = dualAspectPairFreqs[pairID]
            pair = pairFreq[1]
            i = pair[1]
            j = pair[2]
            aspectI = aspectIDs[i]
            aspectJ = aspectIDs[j]
            if haskey(subRegions, (aspectI, aspectJ))
                push!(subRegions[(aspectI, aspectJ)], pairID)
            else
                subRegions[(aspectI, aspectJ)] = [pairID]
            end
        end

        for (aspectPair, group) in subRegions
            push!(groups, group)
        end
    end
    groups
end

function getRegionalConns(sparseConns::Vector{Dict{Int,Float64}}, region::Vector{Int})::Vector{Dict{Int,Float64}}
    regionalConns = [Dict{Int,Float64}() for i in region]
    idMap = Dict{Int,Int}()
    for (localID, globalID) in enumerate(region)
        idMap[globalID] = localID
    end
    for (localID, globalID) in enumerate(region)
        for (globalNeighbor, conn) in sparseConns[globalID]
            if haskey(idMap, globalNeighbor)
                regionalConns[localID][idMap[globalNeighbor]] = conn
            end
        end
    end
    regionalConns
end

maxIter = 1000
maxPairsInGroup = 10000
nonStopWords, nonStopWordIDs = loadNonStopWords("DualAspectTopicRepresentation/workspace/NIPS.nonStopWords.txt")
dualAspectPairFreqs = loadDualAspectPairFreqs("DualAspectTopicRepresentation/workspace/NIPS.DAPairFreq.txt")
sparseConns = loadConns("DualAspectTopicRepresentation/workspace/NIPS.DAPairConns.txt", length(dualAspectPairFreqs))
aspect1, aspectIDs1 = loadAspects("DualAspectTopicRepresentation/workspace/NIPS.aspects.1.txt", nonStopWordIDs)
aspect2, aspectIDs2 = loadAspects("DualAspectTopicRepresentation/workspace/NIPS.aspects.2.txt", nonStopWordIDs)
aspect3, aspectIDs3 = loadAspects("DualAspectTopicRepresentation/workspace/NIPS.aspects.3.txt", nonStopWordIDs)
regions0 = [[i for i in 1:length(dualAspectPairFreqs)]]
regions1 = divideGroups(dualAspectPairFreqs, regions0, aspectIDs3, maxPairsInGroup)
regions2 = divideGroups(dualAspectPairFreqs, regions1, aspectIDs2, maxPairsInGroup)
regions3 = divideGroups(dualAspectPairFreqs, regions2, aspectIDs1, maxPairsInGroup)
topics = Vector{Set{Int}}()
numRegions = length(regions3)
for (iRegion, region) in enumerate(regions3)
    println("finding topics in region $iRegion of $numRegions ...")
    regionalConns = getRegionalConns(sparseConns, region)
    regionalTopicIDs, regionalTopics = findDualAspectTopics(regionalConns, maxIter)
    for regionalTopic in regionalTopics
        topic = Set{Int}()
        for regionalID in regionalTopic
            push!(topic, region[regionalID])
        end
        push!(topics, topic)
    end
end
saveTopics("DualAspectTopicRepresentation/workspace/NIPS.topics.txt", topics, nonStopWords, dualAspectPairFreqs)

max_stars_repo_path: src/05 find_topics.jl

max_issues_repo_path: src/05 find_topics.jl

max_forks_repo_path: src/05 find_topics.jl


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

# camera/acquisition.jl: Camera object acquistion features

"""
  isrunning(::Camera) -> Bool

  Determine if the camera is currently acquiring images.
"""
function isrunning(cam::Camera)
  pbIsStreaming = Ref(bool8_t(false))
  spinCameraIsStreaming(cam, pbIsStreaming)
  return (pbIsStreaming[] == 0x01)
end

"""
  start!(::Camera)

  Start acquistion on specified camera.
"""
function start!(cam::Camera)
  spinCameraBeginAcquisition(cam)
  return cam
end

"""
  stop!(::Camera)

  Stop acquistion on specified camera.
"""
function stop!(cam::Camera)
  spinCameraEndAcquisition(cam)
  return cam
end

"""
  triggermode(::Camera) -> String

  Camera trigger mode.
"""
triggermode(cam::Camera) = get(SpinEnumNode(cam, "TriggerMode"))

"""
  triggermode!(::Camera, ::String) -> String

  Set camera trigger mode, returns set mode.
"""
triggermode!(cam::Camera, mode) = set!(SpinEnumNode(cam, "TriggerMode"), mode)

"""
  triggersource!(::Camera, ::AbstractString) -> String

  Set camera to use specified trigger source. Trigger mode is disabled during
  update and restored after new source is set. Method returns set trigger source.
"""
function triggersource!(cam::Camera, src)
  initmode = triggermode(cam)   # Save current trigger mode
  setsrc = set!(SpinEnumNode(cam, "TriggerSource"), src)
  triggermode!(cam, initmode)   # Restore initial trigger mode
  return setsrc
end

"""
  triggersource(::Camera) -> String

  Return current camera trigger source.
"""
triggersource(cam::Camera) = get(SpinEnumNode(cam, "TriggerSource"))

"""
  trigger!(::Camera)

  Emit software trigger to specified camera.
"""
function trigger!(cam::Camera)

  if triggersource(cam) != "Software"
    @error "Camera is not set to software trigger source"
  end

  if triggermode(cam) != "On"
    @error "Camera is not set to trigger mode"
  end

  hNodeMap = Ref(spinNodeMapHandle(C_NULL))
  spinCameraGetNodeMap(cam, hNodeMap)

  hTriggerSoftware = Ref(spinNodeHandle(C_NULL))
  spinNodeMapGetNode(hNodeMap[], "TriggerSoftware", hTriggerSoftware);
  spinCommandExecute(hTriggerSoftware[])

end

"""
  exposure(::Camera) -> Float, String

  Camera exposure time and mode.
"""
exposure(cam::Camera) = (get(SpinFloatNode(cam, "ExposureTime")), get(SpinEnumNode(cam, "ExposureAuto")))

"""
  exposure!(::Camera)

  Activate (continuous) automatic exposure control on specified camera.
"""
exposure!(cam::Camera) = set!(SpinEnumNode(cam, "ExposureAuto"), "Continuous")

"""
  exposure!(::Camera, ::Number) -> Float

  Set exposure time on camera to specified number of microseconds. The requested
  value is clamped to range supported by the camera, and the actual set value
  is returned. This function disables automatic exposure.
"""
function exposure!(cam::Camera, t)
  set!(SpinEnumNode(cam, "ExposureAuto"), "Off")
  set!(SpinFloatNode(cam, "ExposureTime"), t)
end

"""
  exposure_range(::Camera) -> (Float, Float)

  Exposure time limits in microseconds.
"""
exposure_limits(cam::Camera) = range(SpinFloatNode(cam, "ExposureTime"))

"""
  autoexposure_limits!(::Camera, (lower, upper); clampwarn=true) -> (Float, Float)

  Write lower and upper limits of the Auto Exposure Time (us) value.
  Values exceeding range are clamped to the allowable extrema and returned, and
  a warning is given, which can be disabled with `clampwarn=false`.
"""
function autoexposure_limits!(cam::Camera, lims; clampwarn=true)
  set!(SpinFloatNode(cam, cam.names["AutoExposureTimeLowerLimit"]), lims[1], clampwarn=clampwarn)
  set!(SpinFloatNode(cam, cam.names["AutoExposureTimeUpperLimit"]), lims[2], clampwarn=clampwarn)
  autoexposure_limits(cam)
end

"""
  autoexposure_limits(::Camera) -> (Float, Float)

  Lower and upper limits of the Auto Exposure Time (us) value.
"""
function autoexposure_limits(cam::Camera)
  (get(SpinFloatNode(cam, cam.names["AutoExposureTimeLowerLimit"])),
   get(SpinFloatNode(cam, cam.names["AutoExposureTimeUpperLimit"])))
end

"""
  autoexposure_control_priority!(cam::Camera, priority)

  Select whether to adjust gain or exposure first. `priority` can be one of `"Gain"`, `"ExposureTime"`.
  See Spinnaker SDK docs section 12.6.2.6 `enum AutoExposureControlPriorityEnums` for more information.
"""
function autoexposure_control_priority!(cam::Camera, priority)
  set!(SpinEnumNode(cam, "AutoExposureControlPriority"), priority)
  autoexposure_control_priority(cam)
end

"""
  autoexposure_control_priority(cam::Camera)

  Get the auto-exposure control priority. One of `"Gain"`, `"ExposureTime"`.
  See Spinnaker SDK docs section 12.6.2.6 `enum AutoExposureControlPriorityEnums` for more information.
"""
function autoexposure_control_priority(cam::Camera)
  get(SpinEnumNode(cam, "AutoExposureControlPriority"))
end

"""
  framerate(::Camera) -> Float

  Camera frame rate.
"""
function framerate(cam::Camera)
  get(SpinFloatNode(cam, "AcquisitionFrameRate"))
end

"""
  framerate!(::Camera)

  Activate (continuous) automatic framerate control on specified camera.
"""
framerate!(cam::Camera) = set!(SpinEnumNode(cam, "AcquisitionFrameRateAuto"), "Continuous")

"""
  framerate!(::Camera, ::Number) -> Float

  Set framerate on camera to specified number of frames per second. The requested
  value is clamped to range supported by the camera, and the actual set value
  is returned.
"""
function framerate!(cam::Camera, fps)
  try
    set!(SpinEnumNode(cam, "AcquisitionFrameRateAuto"), "Off")
  catch e
  end
  set!(SpinBooleanNode(cam, cam.names["AcquisitionFrameRateEnabled"]), true)
  set!(SpinFloatNode(cam, "AcquisitionFrameRate"), fps)
end

"""
  framerate_limits(::Camera) -> (Float, Float)

  Framerate limits in microseconds.
"""
framerate_limits(cam::Camera) = range(SpinFloatNode(cam, "AcquisitionFrameRate"))


text:
using DataFrames
using GLM
using RCall
using Test

@testset "Cubic spline matches at points" begin
    using Glissa
    x = collect(0:0.1:pi)
    f = cos.(x)
    cs = cubic_spline(x, f, [0, 0]; slope=FreeSlope())
    for i1 in 1:length(x)
        @test abs(cs(x[i1]) - f[i1]) < 1e-9
    end
end

function points_slope(logx, logy)
    df = DataFrame(logx = logx, logy = logy)
    model = @formula(logy ~ logx)
    regress = lm(model, df)
    slope = coef(regress)[2]
    slope
end

@testset "Cubic spline error reduces as dx reduces" begin
    using Glissa
    Δ = (1e-5) * 2 .^ (0:14)
    maxerr = zeros(Float64, length(Δ))
    for j in 1:length(Δ)
        x = collect(0:(Δ[j]):(pi/2))
        f = cos.(x)
        cs = cubic_spline(x, f, [0, 0]; slope=FreeSlope())

        xx = collect(0:(Δ[j]/10):(pi/2))
        mm = zero(Float64)
        for x1 in xx
            mm = max(mm, abs(cs(x1) - cos(x1)))
        end
        maxerr[j] = mm
    end
    slope = points_slope(log.(Δ), log.(maxerr))
    @test 0.7 < slope
    @test slope < 1.0
end

@testset "Cubic spline is continuous" begin
    using Glissa
    Δ = 0.1
    x = collect(0:Δ:pi)
    f = sin.(x)
    cs = cubic_spline(x, f, [0, 0]; slope=FreeSlope())
    for i1 in 1:(length(x) - 1)
        d = Δ / 10
        prev = 100
        # Assert values approach the expected value.
        for xl in range(x[i1] + Δ / 10, length = 5, stop = x[i1] + 1e-9)
            v = abs(sin(x[i1]) - cs(xl))
            @test v <= prev
            prev = v
        end
        # Assert they get close-enough for the final value.
        @test prev < 1e-5
    end
    for i1 in 2:length(x)
        d = Δ / 10
        prev = 100
        # Look at approach from below.
        for xl in range(x[i1] - Δ / 10, length = 5, stop = x[i1] - 1e-9)
            v = abs(sin(x[i1]) - cs(xl))
            @test v <= prev
            prev = v
        end
        @test prev < 1e-5
    end
end

@testset "Cubic spline is continuous for uneven x" begin
    using Glissa
    Δ = 0.1
    x = vcat([0.0, 0.1, 0.2, 0.4, 0.5, 1], 2:5)
    f = sin.(x)
    cs = cubic_spline(x, f, [0, 0]; slope=FreeSlope())
    for i1 in 1:(length(x) - 1)
        d = Δ / 10
        prev = 100
        # Assert values approach the expected value.
        for xl in range(x[i1] + Δ / 10, length = 5, stop = x[i1] + 1e-9)
            v = abs(sin(x[i1]) - cs(xl))
            @test v <= prev
            prev = v
        end
        # Assert they get close-enough for the final value.
        @test prev < 1e-5
    end
    for i1 in 2:length(x)
        d = Δ / 10
        prev = 100
        # Look at approach from below.
        for xl in range(x[i1] - Δ / 10, length = 5, stop = x[i1] - 1e-9)
            v = abs(sin(x[i1]) - cs(xl))
            @test v <= prev
            prev = v
        end
        @test prev < 1e-5
    end
end

@testset "Cubic spline overshoot" begin
    using Glissa
    # RPN 15A data from Fritsch-Carlson
    x = [7.99, 8.09, 8.19, 8.7, 9.2, 10, 12, 15, 20]
    f = [0, 2.76429e-5, 4.37498e-2, 0.169183, 0.469428, 0.943740, 0.998636, 0.999919, 0.999994]
    cs = cubic_spline(x, f, [0.0, 0.0]; slope = Monotonic())
    xx = 8.08:0.0001:8.1
    # plot(xx, cs.(xx))
    ϵ = 1e-8
    for i in 2:(length(x) - 1)
        left = (cs(x[i] - ϵ) - cs(x[i] - 2ϵ)) / ϵ
        right = (cs(x[i] + 2ϵ) - cs(x[i] + ϵ))/ϵ
        d = abs(left - right)
        @test d < 1e-6
    end
    ϵ = 1e-5
    second_derivative = zeros(Float64, length(x) - 2)
    for i in 2:(length(x) - 1)
        left = (cs(x[i] - 3ϵ) - 2*cs(x[i] - 2ϵ) + cs(x[i] - ϵ)) / ϵ^2
        right = (cs(x[i] + ϵ) - 2*cs(x[i] + 2ϵ) + cs(x[i] + 3ϵ)) / ϵ^2
        second_derivative[i - 1] = abs(left - right)
    end
    @test second_derivative[1] > 5 * second_derivative[5]
end

@testset "Hyman criterion matches splinefun in R" begin
    using Glissa
    using RCall
    # this hyman filter comes from the R splinefun source code.
    # x is the knots. y is values. b is gradient.
    R"""hyman_filter <- function(x, y, b)
    {
        n <- length(x)
        ss <- diff(y) / diff(x)
        S0 <- c(ss[1L], ss)
        S1 <- c(ss, ss[n-1L])
        t1 <- pmin(abs(S0), abs(S1))
        sig <- b
        ind <- S0*S1 > 0
        sig[ind] <- S1[ind]
        ind <- sig >= 0
        if(sum(ind)) b[ind] <- pmin(pmax(0, b[ind]), 3*t1[ind])
        ind <- !ind
        if(sum(ind)) b[ind] <- pmax(pmin(0, b[ind]), -3*t1[ind])
        b
    }
    """
    x = collect(0:0.1:π)
    y = sin.(x)
    yp = cos.(x)
    hf = collect(R"""hyman_filter($x, $y, $yp)""")
    Glissa.project_to_monotonicity!(x, y, yp)
    @test maximum(abs.(yp .- hf)) < 1e-12

    # This example is constructed by hand to need monotonic smoothing.
    x = [7.99, 8.09, 8.19, 8.7, 9.2, 10, 12, 15, 20]
    f = [0, 2.76429e-5, 4.37498e-2, 0.169183, 0.469428, 0.943740, 0.998636, 0.999919, 0.999994]
    fp = similar(f)
    Glissa.global_derivatives!(x, f, fp)
    hf = collect(R"""hyman_filter($x, $f, $fp)""")
    # Assert that the uncorrected slopes are different from corrected ones.
    @test maximum(abs.(hf - fp)) > 1e-10
    Glissa.project_to_monotonicity!(x, f, fp)
    # After correction, they agree with the R version.
    @test maximum(abs.(fp .- hf)) < 1e-12

    R"""spl_coef_conv <- function(x, y, b)
    {
        n <- length(x)
        h <- diff(x); y <- -diff(y)
        b0 <- b[-n]; b1 <- b[-1L]
        cc <- -(3*y + (2*b0 + b1)*h) / h^2
        c1 <- (3*y[n-1L] + (b0[n-1L] + 2*b1[n-1L])*h[n-1L]) / h[n-1L]^2
        c <- c(cc, c1)
        dd <- (2*y/h + b0 + b1) / h^2
        d <- c(dd, dd[n-1L])
        list(c, d)
    }
    """
    cd = R"""spl_coef_conv($x, $f, $fp)"""
    c = collect(cd[1])
    d = collect(cd[2])
    C = zeros(Float64, 4, length(x) - 1)
    Glissa.hyman_coefficients!(x, f, fp, C)
    # The R functions define vectors that are 1 longer.
    @test maximum(abs.(C[3, :] - c[1:(end-1)])) < 1e-12
    @test maximum(abs.(C[4, :] - d[1:(end-1)])) < 1e-12
end

@testset "Cubic spline is compatible with forward diff" begin
    using Glissa
    using ForwardDiff
    x = [7.99, 8.09, 8.19, 8.7, 9.2, 10, 12, 15, 20]
    x2 = collect(8:0.01:10)
    y = [0, 2.76429e-5, 4.37498e-2, 0.169183, 0.469428, 0.943740, 0.998636, 0.999919, 0.999994]

    function g(f::AbstractVector{T}) where {T}
        cs = Glissa.cubic_spline(x, f, zeros(T, 2); slope = Monotonic())
        cs.(x2)
    end

    interp = g(y)
    jacob = ForwardDiff.jacobian(g, y)
    @test size(jacob) == (201, length(x))
    @test jacob[1, 1] > 0
    @test abs(jacob[1, 9]) < 1e-7
end

@testset "Cubic spline is compatible with reverse diff" begin
    using Glissa
    using ReverseDiff
    x = [7.99, 8.09, 8.19, 8.7, 9.2, 10, 12, 15, 20]
    x2 = collect(8:0.01:10)
    f = [0, 2.76429e-5, 4.37498e-2, 0.169183, 0.469428, 0.943740, 0.998636, 0.999919, 0.999994]
    function g(f::AbstractVector{T}, x1, x2) where {T}
        cs = Glissa.cubic_spline(x1, f, zeros(T, 2); slope = Monotonic())
        cs.(x2)
    end
    interp = g(f, x, x2)
    f_tape = ReverseDiff.GradientTape(g, (f, x, x2))
    compiled_f_tape = ReverseDiff.compile(f_tape)
    inputs = (f, x, x2)
    cfg = ReverseDiff.GradientConfig(inputs)
    results = similar(interp)
    # This doesn't work, and I'd like it to work.
    # ReverseDiff.gradient!(results, compiled_f_tape, inputs)
end

@testset "Cubic spline integrates a flat line" begin
    using Glissa
    using Random
    using Distributions
    x = Float64[0, 1, 3, 5]
    y = Float64[1, 1, 1, 1]
    cs = Glissa.cubic_spline(x, y, [0.0, 0.0]; slope = Monotonic())
    @test abs(cs(2.5) - 1) < 1e-6
    exact(z) = z
    rng = MersenneTwister(9234724)
    for j in 1:100
        a = rand(rng, Uniform(0, 5))
        b = rand(rng, Uniform(a, 5))
        found = Glissa.integrate(cs, a, b)
        @test abs(found - exact(b) + exact(a)) < 1e-9
    end
end

@testset "Cubic spline integrates a sloped line" begin
    using Glissa
    using Random
    using Distributions
    x = Float64[0, 1, 3, 5]
    y = Float64[0, 1, 3, 5]
    cs = Glissa.cubic_spline(x, y, [1.0, 1.0]; slope = Monotonic())
    @test abs(cs(2.5) - 2.5) < 1e-6
    exact(z) = 0.5 * z^2
    rng = MersenneTwister(9292734)
    for j in 1:100
        a = rand(rng, Uniform(0, 5))
        b = rand(rng, Uniform(a, 5))
        found = Glissa.integrate(cs, a, b)
        @test abs(found - exact(b) + exact(a)) < 1e-9
    end
end

@testset "Cubic spline integrates a quadratic line" begin
    using Glissa
    using Random
    using Distributions
    Δ = 0.1
    x = collect(0:Δ:5)
    y = x.^2
    cs = Glissa.cubic_spline(x, y, [0.0, 10.0]; slope = Monotonic())
    for q in [0.3, 1.2, 4.7, 4.99]
        @test abs(cs(q) - q^2) < Δ^3
    end
    exact(z) = z^3 / 3
    rng = MersenneTwister(89192734)
    for j in 1:100
        a = rand(rng, Uniform(0, 5))
        b = rand(rng, Uniform(a, 5))
        found = Glissa.integrate(cs, a, b)
        @test abs(found - exact(b) + exact(a)) < 8*Δ^3
    end
end

@testset "Cubic spline integrates a cosine" begin
    using Glissa
    using Random
    using Distributions
    Δ = π / 30
    x = collect(0:Δ:π)
    y = cos.(x)
    cs = Glissa.cubic_spline(x, y, [0.0, 0.0]; slope = Monotonic())
    for q in [0.3, 1.2, 1.7, 2.99]
        @test abs(cs(q) - cos(q)) < Δ^3
    end
    exact(z) = sin(z)
    rng = MersenneTwister(89192734)
    for j in 1:100
        a = rand(rng, Uniform(0, 3))
        b = rand(rng, Uniform(a, 3))
        found = Glissa.integrate(cs, a, b)
        @test abs(found - exact(b) + exact(a)) < Δ^3
    end
end

@testset "Cubic spline derivative" begin
    using Glissa
    cs = CubicSpline{Int,Float64}([0, 1, 2],
        Float64[
            1 2;
            3 7;
            4 0;
            1 -1
        ]
    )
    csprime = Glissa.derivative(cs)
    @test size(csprime.c) == (3, 2)
    @test csprime.c ≈ Float64[
        3 7;
        8 0;
        3 -3
    ]
end

@testset "Cubic spline multiplication" begin
    using Glissa
    cs1 = CubicSpline{Int,Float64}([0, 1, 2],
        Float64[
            1 2;
            3 7;
            4 0
        ]
    )
    using Glissa
    cs2 = CubicSpline{Int,Float64}([0, 1, 2],
        Float64[
            2 2;
            1 -2
        ]
    )
    cs3 = cs1 * cs2
    @test cs3.c ≈ Float64[
        2 4;
        7 10;
        11 -14;
        4 0
    ]
end


text:
using Pkg

# We need these for the doctests. We install them before we load any
# package so that we don't run into precompile problems
Pkg.add(PackageSpec(url="https://github.com/fund-model/MimiFUND.jl", rev="master"))
Pkg.add(PackageSpec(url="https://github.com/anthofflab/MimiDICE2010.jl", rev="master"))

using Mimi
import Electron
using Test
using Documenter

Electron.prep_test_env()

@testset "Mimi" begin

    @info("test_main.jl")
    @time include("test_main.jl")

    @info("test_composite.jl")
    @time include("test_composite.jl")

    @info("test_composite_parameters.jl")
    @time include("test_composite_parameters.jl")

    @info("test_main_variabletimestep.jl")
    @time include("test_main_variabletimestep.jl")

    @info("test_broadcast.jl")
    @time include("test_broadcast.jl")

    @info("test_metainfo.jl")
    @time include("test_metainfo.jl")

    @info("test_metainfo_variabletimestep.jl")
    @time include("test_metainfo_variabletimestep.jl")

    @info("test_references.jl")
    @time include("test_references.jl")

    @info("test_units.jl")
    @time include("test_units.jl")

    @info("test_model_structure.jl")
    @time include("test_model_structure.jl")

    @info("test_model_structure_variabletimestep.jl")
    @time include("test_model_structure_variabletimestep.jl")

    @info("test_delete.jl")
    @time include("test_delete.jl")

    @info("test_replace_comp.jl")
    @time include("test_replace_comp.jl")

    @info("test_tools.jl")
    @time include("test_tools.jl")

    @info("test_parameter_labels.jl")
    @time include("test_parameter_labels.jl")

    @info("test_parametertypes.jl")
    @time include("test_parametertypes.jl")

    @info("test_defaults.jl")
    @time include("test_defaults.jl")

    @info("test_marginal_models.jl")
    @time include("test_marginal_models.jl")

    @info("test_adder.jl")
    @time include("test_adder.jl")

    @info("test_getindex.jl")
    @time include("test_getindex.jl")

    @info("test_getindex_variabletimestep.jl")
    @time include("test_getindex_variabletimestep.jl")

    @info("test_components.jl")
    @time include("test_components.jl")

    @info("test_variables_model_instance.jl")
    @time include("test_variables_model_instance.jl")

    @info("test_getdataframe.jl")
    @time include("test_getdataframe.jl")

    @info("test_mult_getdataframe.jl")
    @time include("test_mult_getdataframe.jl")

    @info("test_clock.jl")
    @time include("test_clock.jl")

    @info("test_timesteps.jl")
    @time include("test_timesteps.jl")

    @info("test_timesteparrays.jl")
    @time include("test_timesteparrays.jl")

    @info("test_dimensions.jl")
    @time include("test_dimensions.jl")

    @info("test_datum_storage.jl")
    @time include("test_datum_storage.jl")

    @info("test_connectorcomp.jl")
    @time include("test_connectorcomp.jl")

    @info("test_firstlast.jl")
    @time include("test_firstlast.jl")

    @info("test_explorer_model.jl") # BROKEN
    @time include("test_explorer_model.jl")

    @info("test_explorer_sim.jl")
    @time include("test_explorer_sim.jl")

    @info("test_explorer_compositecomp.jl")
    @time include("test_explorer_compositecomp.jl")

    @info("mcs/runtests.jl")
    @time include("mcs/runtests.jl")

    @info("doctests")
    @time doctest(Mimi)

    for app in Electron.applications()
        close(app)
    end

end


text:
using NLopt

function myfunc(x::Vector, grad::Vector)
    if length(grad) > 0
        grad[1] = 0
        grad[2] = 0.5/sqrt(x[2])
    end
    return sqrt(x[2])
end

function myconstraint(x::Vector, grad::Vector, a, b)
    if length(grad) > 0
        grad[1] = 3a * (a*x[1] + b)^2
        grad[2] = -1
    end
    (a*x[1] + b)^3 - x[2]
end

opt = Opt(:LD_SLSQP, 2)
opt.lower_bounds = [-Inf, 0.]
opt.xtol_rel = 1e-4

opt.min_objective = myfunc
equality_constraint!(opt, (x,g) -> sum(x)==3, 1e-8)
inequality_constraint!(opt, (x,g) -> myconstraint(x,g,-1,1), 1e-8)

(minf,minx,ret) = optimize(opt, [1.234, 5.678])
numevals = opt.numevals # the number of function evaluations
println("got $minf at $minx after $numevals iterations (returned $ret)")


text:
using ComsolUtils
using Test

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


text:
function get_pull_requests(
    api::GitHub.GitHubAPI,
    repo::GitHub.Repo,
    state::String;
    per_page::Integer=100,
    page_limit::Integer=100,
)
    repo_name = repo.name
    repo_owner = repo.owner.login
    paginated_prs = GitForge.paginate(
        GitHub.get_pull_requests,
        api,
        repo_owner,
        repo_name;
        state=state,
        per_page=per_page,
        page_limit=page_limit,
    )

    return @mock unique(paginated_prs)
end

function get_pull_requests(
    api::GitLab.GitLabAPI,
    project::GitLab.Project,
    state::String;
    per_page::Integer=100,
    page::Integer=1,
)
    project_id = project.id
    paginated_prs = GitForge.paginate(
        GitLab.get_pull_requests, api, project_id; state=state, per_page=per_page, page=page
    )

    return @mock unique(paginated_prs)
end

not_pr_fork(repo::GitHub.Repo, pr::GitHub.PullRequest) = !pr.head.repo.fork
not_pr_fork(repo::GitLab.Project, pr::GitLab.MergeRequest) = repo.id == pr.project_id

my_prs(username::String, pr::GitHub.PullRequest) = lower(pr.user.login) == lower(username)
function my_prs(username::String, pr::GitLab.MergeRequest)
    return lower(pr.author.username) == lower(username)
end

function get_pr_titles(
    forge::Forge, repo::Union{GitHub.Repo,GitLab.Project}, username::AbstractString
)
    state = repo isa GitHub.Repo ? "open" : "opened"
    all_open_prs = @mock get_pull_requests(forge, repo, state)

    return [
        convert(String, strip(pr.title)) for
        (pr, _) in all_open_prs if not_pr_fork(repo, pr) && my_prs(username, pr)
    ]
end


text:
using BenchmarkTools

function callMethod(a::A)
    sum = 0
    sum += a.benchmark()
    sum += a.benchmark2(11, 1.f0)
    sum += a.benchmark2(22)
    sum += a.x
    sum += a.y
    sum += a.z
    return sum
end

function callFunction(a::A)
    sum = 0
    sum += benchmark(a)
    sum += benchmark2(a, 11, 1.f0)
    sum += benchmark2(a, 22)
    sum += a.x
    sum += a.y
    sum += a.z
    return sum
end

function benchmarkClassMethod()
    a = A(1.2f0, 3, false)

    println("----Call method:")
    b = @benchmarkable callMethod($a)
    tune!(b)
    display(run(b))

    println("\n----Call function:")
    b = @benchmarkable callFunction($a)
    tune!(b)
    display(run(b))
    println()

    return 0
end


text:
import DCELs.squeeze_edge

@testset "Squeezing - empty" begin
    u, v = Vertex(0, 0), Vertex(1, 0)
    a1, a2 = Edge(), Edge()
    a1.twin, a2.twin = a2, a1
    a1.dest, a2.dest = u, v

    @test squeeze_edge(u, a1) == a1
end

@testset "Squeezing - line" begin
    u, v = Vertex(0, 0), Vertex(1, 0)
    a1, a2 = Edge(), Edge()

    a1.twin, a2.twin = a2, a1
    a1.prev, a2.prev = a2, a1
    a1.next, a2.next = a2, a1
    a1.dest, a2.dest = u, v

    u.edge, v.edge = a1, a2

    z = Vertex(2, 0)
    h1, h2 = Edge(), Edge()
    h1.twin, h2.twin = h2, h1
    h1.dest, h2.dest = v, z
    @test squeeze_edge(v, h1) == a2
end

@testset "Squeezing - cross" begin
    u, v, w, x, c = Vertex(1, 0), Vertex(0, 1), Vertex(-1, 0), Vertex(0, -1), Vertex(0, 0)
    a1, a2, b1, b2, c1, c2, d1, d2 = [Edge() for _ in 1:8]

    a1.twin, a2.twin = a2, a1
    b1.twin, b2.twin = b2, b1
    c1.twin, c2.twin = c2, c1
    d1.twin, d2.twin = d2, d1

    a1.dest, a2.dest = u, c
    b1.dest, b2.dest = v, c
    c1.dest, c2.dest = w, c
    d1.dest, d2.dest = x, c

    a1.prev, a2.prev = b2, a1
    b1.prev, b2.prev = c2, b1
    c1.prev, c2.prev = d2, c1
    d1.prev, d2.prev = a2, d1

    a1.next, a2.next = a2, d1
    b1.next, b2.next = b2, a1
    c1.next, c2.next = c2, b1
    d1.next, d2.next = d2, c1

    u.edge, v.edge, w.edge, x.edge = a1, b1, c1, d1
    c.edge = a2

    z = Vertex(1, -1)
    h1, h2 = Edge(), Edge()
    h1.twin, h2.twin = h2, h1
    h1.dest, h2.dest = c, z
    @test squeeze_edge(c, h1) == a2
end

@testset "Squeezing - triangle" begin
    u, v, w = Vertex(1, 0), Vertex(0, 1), Vertex(-1, 0)
    a1, a2, b1, b2, c1, c2 = [Edge() for _ in 1:6]

    u.edge, v.edge, w.edge = a1, b1, c1

    a1.twin, a2.twin = a2, a1
    b1.twin, b2.twin = b2, b1
    c1.twin, c2.twin = c2, c1

    a1.dest, b1.dest, c1.dest = u, v, w
    a2.dest, b2.dest, c2.dest = w, u, v

    a1.next, b1.next, c1.next = b1, c1, a1
    a1.prev, b1.prev, c1.prev = c1, a1, b1

    a2.next, b2.next, c2.next = c2, a2, b2
    a2.prev, b2.prev, c2.prev = b2, c2, a2

    x = Vertex(2, 0)
    h1, h2 = Edge(), Edge()
    h1.twin, h2.twin = h2, h1
    h1.dest, h2.dest = u, x
    @test squeeze_edge(u, h1) == b2

    x = Vertex(0, 0.5)
    h1, h2 = Edge(), Edge()
    h1.twin, h2.twin = h2, h1
    h1.dest, h2.dest = u, x
    @test squeeze_edge(u, h1) == a1
end


text:
using Documenter, BLASFEO

makedocs(;
    modules=[BLASFEO],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/JuliaEmbedded/BLASFEO.jl/blob/{commit}{path}#L{line}",
    sitename="BLASFEO.jl",
    authors="Ian McInerney, Imperial College London",
    assets=String[],
)

deploydocs(;
    repo="github.com/JuliaEmbedded/BLASFEO.jl",
)


text:
# ---
# title: 365. Water and Jug Problem
# id: problem365
# author: Indigo
# date: 2021-06-14
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/water-and-jug-problem/description/>
# hidden: true
# ---
#
# You are given two jugs with capacities _x_ and _y_ litres. There is an
# infinite amount of water supply available. You need to determine whether it is
# possible to measure exactly _z_ litres using these two jugs.
#
# If _z_ liters of water is measurable, you must have _z_ liters of water
# contained within **one or both buckets** by the end.
#
# Operations allowed:
#
#   * Fill any of the jugs completely with water.
#   * Empty any of the jugs.
#   * Pour water from one jug into another till the other jug is completely full or the first jug itself is empty.
#
# **Example 1:** (From the famous [_" Die Hard"_
# example](https://www.youtube.com/watch?v=BVtQNK_ZUJg))
#
#
#
#     Input: x = 3, y = 5, z = 4
#     Output: True
#
#
# **Example 2:**
#
#
#
#     Input: x = 2, y = 6, z = 5
#     Output: False
#
#
#
#
# **Constraints:**
#
#   * `0 <= x <= 10^6`
#   * `0 <= y <= 10^6`
#   * `0 <= z <= 10^6`
#
#
## @lc code=start
using LeetCode

function can_measure_water(jug1_cap::Int, jug2_cap::Int, target::Int)
    target > jug1_cap + jug2_cap && return false
    return target % gcd(jug2_cap, jug1_cap) == 0 ? true : false
end
## @lc code=end


text:
# This file is a part of JuliaFEM
# License is MIT: see https://github.com/JuliaFEM/FEModels.jl/blob/master/LICENSE

module FEModels

using GZip

"""
    fe_download(model, mesh)

Download finite element mesh from FEModels repository.
"""
function fe_download(model, mesh, version)
#   url = "https://github.com/JuliaFEM/FEModels/blob/$version/$model/$mesh.gz?raw=true"
    url = "https://raw.githubusercontent.com/JuliaFEM/FEModels/$version/$model/$mesh.gz"
    isdir(model) || mkdir(model)
    local_file = joinpath(model, mesh)
    local_file_packed = "$local_file.gz"
    download(url, local_file_packed)
    fh_packed = GZip.open(local_file_packed)
    data = read(fh_packed)
    close(fh_packed)
    fh_unpacked = open(local_file, "w")
    write(fh_unpacked, data)
    close(fh_unpacked)
    return local_file
end

"""
    fe_list()

List all finite element models / meshes in FEModels repository.
"""
function fe_list()
    error("Not implemented.")
end

export fe_download

end


text:
module DensityScatter

include("densitycalc.jl")

export performBinning, calc_densities,fa, smoothe;

end # module


text:
"""
    jacobi(x, p::Integer, α, β)

Evaluate the Legendre polynomial with parameters `α`, `β` of degree `p` at `x`
using the three term recursion [Karniadakis and Sherwin, Spectral/hp Element
Methods for CFD, Appendix A].
"""
function jacobi(x, p::Integer, α, β)
    T = typeof( (2+α+β)*x / 2 )
    a = one(T)
    b = ((2+α+β)*x + α - β) / 2

    if p <= 0
        return a
    elseif p == 1
        return b
    end

    for n in 2:p
        a1 = 2n*(n+α+β)*(2n-2+α+β)
        a2 = (2n-1+α+β)*(α+β)*(α-β)
        a3 = (2n-2+α+β)*(2n-1+α+β)*(2n+α+β)
        a4 = 2*(n-1+α)*(n-1+β)*(2n+α+β)
        a, b = b, ( (a2+a3*x)*b - a4*a ) / a1
    end

    b
end

"""
    jacobi_and_derivative(x, p::Integer, α, β)

Evaluate the Jacobi polynomial with parameters `α`, `β` of degree `p` and its
derivative at `x` using the three term recursion [Karniadakis and Sherwin,
Spectral/hp Element Methods for CFD, Appendix A].
"""
function jacobi_and_derivative(x, p::Integer, α, β)
    T = typeof( (2+α+β)*x / 2 )
    # Coefficients for the polynomial...
    a = one(T)
    b = ((2+α+β)*x + α - β) / 2
    # ... and for the derivative. The derivative is computed via eq. (A.1.8)
    # of Karniadakis and Sherwin as Jacobi polynomial with parameters α+1, β+1.
    aa = one(T)
    bb = ((4+α+β)*x + α - β) / 2

    if p <= 0
        return a, zero(T)
    elseif p == 1
        return b, (2+α+β)/2
    end

    for n in 2:p
        a1 = 2n*(n+α+β)*(2n-2+α+β)
        a2 = (2n-1+α+β)*(α+β)*(α-β)
        a3 = (2n-2+α+β)*(2n-1+α+β)*(2n+α+β)
        a4 = 2*(n-1+α)*(n-1+β)*(2n+α+β)
        a, b = b, ( (a2+a3*x)*b - a4*a ) / a1

        b1 = 2n*(n+2+α+β)*(2n+α+β)
        b2 = (2n+1+α+β)*(α+β+2)*(α-β)
        b3 = (2n+α+β)*(2n+1+α+β)*(2n+2+α+β)
        b4 = 2*(n+α)*(n+β)*(2n+2+α+β)
        aa, bb = bb, ( (b2+b3*x)*bb - b4*aa ) / b1
    end

    b, (α+β+p+1)*aa/2
end

"""
    gauss_jacobi_nodes_and_weights(p, α, β, T=Float64::Type, tol=4*eps(T), maxit=100)

Compute the Gauss-Jacobi nodes and weights for polynomials of degree `p` with
parameters `α`, `β` using the scalar type `T`, tolerance `tol` and maximal number
of Newton iterations `maxit` [Karniadakis and Sherwin, Spectral/hp Element
Methods for CFD, Appendix B].
"""
function gauss_jacobi_nodes_and_weights(p, α, β, T=Float64::Type, tol=4*eps(T), maxit=1000)
    T = promote_type(typeof(α), typeof(β), T)
    nodes = Vector{T}(undef, p+1)
    weights = Vector{T}(undef, p+1)
    for j in 0:p
        x = -cospi(T(2j+1)/(2p+2))
        if j > 0
            x = (x + nodes[j]) / 2
        end
        for k in 1:maxit
            s = zero(T)
            for l in 1:j
                s += 1 / (x - nodes[l])
            end
            pol, der = jacobi_and_derivative(x, p+1, α, β)
            Δ = -pol / (der - s*pol)
            x = x + Δ
            abs(Δ) <= tol*abs(x) && break
        end
        pol, der = jacobi_and_derivative(x, p+1, α, β)
        nodes[j+1] = x
        weights[j+1] = 2^(α+β+1) * gamma(α+p+2) * gamma(β+p+2) /
                        (gamma(p+2) * gamma(α+β+p+2) * (1-x^2) * der^2)
    end

    return nodes, weights
end

"""
    jacobi_vandermonde(nodes, α, β)

Computes the Vandermonde matrix with respect to the Jacobi polynomials with
parameters `α`, `β` and the nodal basis on `nodes`.
The Vandermonde matrix `V` is the transformation matrix from the modal Jacobi
basis to the nodal Lagrange basis associated with `nodes`.
"""
function jacobi_vandermonde(nodes::AbstractVector, α, β)
    T = eltype(nodes)
    pp1 = length(nodes)
    V = Array{T}(undef, pp1, pp1)
    for j in 1:pp1, (i,x) in enumerate(nodes)
      V[i, j] = jacobi(x, j-1, α, β)
    end
    V
end

jacobi_vandermonde(basis::NodalBasis{Line}, α, β) = jacobi_vandermonde(basis.nodes, α, β)


text:
pages = [
        "DiffEqParamEstim.jl: Parameter Estimation for Differential Equations" => "index.md",
        "Tutorials" => Any[
          "tutorials/ODE_inference.md",
          "tutorials/global_optimization.md",
          "tutorials/jump.md",
          "tutorials/generalized_likelihood.md",
          "tutorials/stochastic_evaluations.md"
        ],
        "Methods" => Any[
          "methods/recommended_methods.md",
          "methods/optimization_based_methods.md",
          "methods/collocation_loss.md",
          "methods/alternative_objectives.md"
        ],
    ]


text:
using Nystrom
using Documenter

DocMeta.setdocmeta!(Nystrom, :DocTestSetup, :(using Nystrom); recursive=true)

makedocs(;
    modules=[Nystrom],
    authors="Luiz M. Faria <maltezfaria@gmail.com> and contributors",
    repo="https://github.com/WaveProp/Nystrom.jl/blob/{commit}{path}#{line}",
    sitename="Nystrom.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://WaveProp.github.io/Nystrom.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/WaveProp/Nystrom.jl",
    devbranch="main",
)

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
include("../src/Rand.jl")
include("../src/Padding.jl")
include("../src/Iteration.jl")
include("../src/BlockCiphers.jl")
include("../src/CipherModes.jl")
include("testing.jl")

importall Padding
importall CipherModes
importall BlockCiphers

println("testing encryption modes:")

for key_len = [128, 192, 256]
    for plain = {zeros(Uint8, 20), zeros(Uint8, 0), zeros(Uint8, 1)}
        for key = [AES(zeros(Uint8, int32(key_len/8))), DES(zeros(Uint8, 8)), DES3(zeros(Uint8, 8 * 3)), DES3ede(zeros(Uint8, 8 * 3))]
            print("testing ", typeof(key), " ")

            base_len = length(plain)
            padded_len = block_size(key) * int64(ceil(length(plain) / block_size(key)))
            if(length(plain) % block_size(key) == 0)
                padded_len += block_size(key)
            end
            iv_len = block_size(key)

            out = ecb_encrypt(Pkcs7, key, plain)
            @assert_equal(padded_len + 0, length(out))
            @assert_not_equal(plain, out)
            @assert_equal(plain, ecb_decrypt(Pkcs7, key, out))

            out = cbc_encrypt(Pkcs7, key, plain)
            @assert_equal(padded_len + iv_len, length(out))
            @assert_not_equal(plain, out)
            @assert_equal(plain, cbc_decrypt(Pkcs7, key, out))

            out = pcbc_encrypt(Pkcs7, key, plain)
            @assert_equal(padded_len + iv_len, length(out))
            @assert_not_equal(plain, out)
            @assert_equal(plain, pcbc_decrypt(Pkcs7, key, out))

            out = ctr_encrypt(key, 31337, plain)
            @assert_equal(base_len, length(out)) # no padding :D
            if(length(plain) != 0)
                @assert_not_equal(plain, out)
            end
            @assert_equal(plain, ctr_decrypt(key, 31337, out))

            out = cfb_encrypt(key, plain)
            @assert_equal(base_len + iv_len, length(out)) # no padding + iv though :D
            @assert_not_equal(plain, out)
            @assert_equal(plain, cfb_decrypt(key, out))

            out = ofb_encrypt(key, plain)
            @assert_equal(base_len + iv_len, length(out)) # no padding + iv though :D
            @assert_not_equal(plain, out)
            @assert_equal(plain, ofb_decrypt(key, out))
            println()
        end
    end
end

println()


text:
module XML2011

using Dates
using EzXML
using Printf

include("types.jl")
include("moeda.jl")
include("posicao.jl")
include("pais.jl")
include("write.jl")
include("read.jl")
include("validate.jl")

end # module


text:
abstract type AbstractMDIntegrator end

struct LeapFrog <: AbstractMDIntegrator
    Δt::Float64
end

mutable struct VelocityVerlet <: AbstractMDIntegrator
    Δt::Float64
end

LeapFrog(Δt::Unitful.Time) = LeapFrog(austrip(Δt))

function propagate!(coordinates, velocities, potential, mass, integrator::LeapFrog)
    velocities .+= force(potential, coordinates).*integrator.Δt ./mass
    coordinates .+= velocities.*integrator.Δt
    return  coordinates, velocities
end

function propagate!(r1, v1, f1, r0, v0, f0, potential, mass, integrator::VelocityVerlet)
    r1 .= r0 .+ v0.*integrator.Δt .+ integrator.Δt^2*f1./mass
    v1 .= v0 .+ integrator.Δt.*(f0 .+ force!(f1,r1))./(2*mass)
    return r1, v1
end


text:
using DataStructures

import Base: eof, close, position
export FeatureMap, close, value, eof, advance!

type FeatureMap
    reader::BamReader
    feature::Array{Float64,1}
    featureWidth::Int64
    position::Int64
    posQueue::Deque{Int64}
    value::Float64
end

function FeatureMap(reader::BamReader, feature)
    fm = FeatureMap(reader, feature, (length(feature)-1)/2, 0, Deque{Int64}(), 0.0)

    advance!(fm)
    fm
end
close(fm::FeatureMap) = close(fm.reader)
value(fm::FeatureMap) = fm.value
position(fm::FeatureMap) = fm.position
eof(fm::FeatureMap) = fm.position < 0

function advance!(fm::FeatureMap)
    fm.position += 1

    # remove old positions that fell outside the sliding window
    while length(fm.posQueue) > 0 && front(fm.posQueue) < fm.position-fm.featureWidth
        shift!(fm.posQueue)
    end

    # skip ahead to the next spot we have data if our queue is empty
    if length(fm.posQueue) == 0
        fm.position = fm.reader.position-fm.featureWidth
    end

    # add new positions in the sliding window
    while fm.reader.position != -1 && fm.reader.position <= fm.position+fm.featureWidth
        push!(fm.posQueue, fm.reader.position)
        advance!(fm.reader)
    end

    # compute the feature value for this position
    fm.value = 0
    for pos in fm.posQueue
        fm.value += fm.feature[pos - fm.position + fm.featureWidth + 1]
    end
end


text:
import MatrixDepot: publish_user_generators, include_generator, Group, FunctionName
using Logging

"random symmetric matrix"
function randsym(::Type{T}, n) where T
 A = zeros(T, n, n)
  for j = 1:n
      for i = j:n
          A[i,j] = randn()
      end
  end
  A = A + tril(A, -1)'
  return A
end
randsym(n) = randsym(Float64, n)
include_generator(FunctionName, "randsym", randsym)
include_generator(Group, :random, randsym)
include_generator(Group, :symmetric, randsym)

# update the database
MatrixDepot.publish_user_generators()

n = rand(1:8)
@test matrixdepot("randsym", n) !== nothing
@test mdinfo("randsym") !== nothing
@test "randsym" in MatrixDepot.mdlist(:random)
@test "randsym" in MatrixDepot.mdlist(:symmetric)

@test_logs min_level=Logging.Warn MatrixDepot.init()
n = rand(1:8)
@test matrixdepot("randsym", n) !== nothing
@test mdinfo("randsym") !== nothing
@test mdinfo("randsym") == Base.Docs.doc(randsym)
@test "randsym" in MatrixDepot.mdlist(:random)
@test "randsym" in MatrixDepot.mdlist(:symmetric)

setgroup!(:testgroup, "randsy*")
@test mdlist(:testgroup) == ["randsym"]
setgroup!(:testgroup, ["*/1138_bus"])
@test mdlist(:testgroup) == ["HB/1138_bus"]
deletegroup!(:testgroup)
@test_throws ArgumentError mdlist(:testgroup)

@test_throws ArgumentError include_generator(Group, :lkjasj, sin)


text:
# Part 1 - parse passports
passes = open(f -> strip.(readlines(f)), "day05_input.txt")
seatids = Set(map(p -> parse(Int32, join([c ∈ ['F', 'L'] ? '0' : '1' for c in p]), base=2), passes))
println("Part 1: Highest seat ID is $(maximum(seatids))")

# Part 2
println("Part 2: seat ID is $(first([s+1 for s in seatids if (s+1) ∉ seatids && (s+2) ∈ seatids]))")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
"""
This function builds a BezierCurve which leaves point p vertically upwards and
arrives point q vertically upwards. It may create a loop if necessary.
It assumes the view is [0,1]. That can be modified using the `xview` and
`yview` keyword arguments (default: `0:1`).
"""
function directed_curve(x1, x2, y1, y2; xview = 0:1, yview = 0:1, root::Symbol = :bottom)
    if root in (:left, :right)
        # flip x/y to simplify
        x1, x2, y1, y2, xview, yview = y1, y2, x1, x2, yview, xview
    end
    x = Float64[x1, x1]
    y = Float64[y1]

    minx, maxx = extrema(xview)
    miny, maxy = extrema(yview)
    dist = sqrt((x2 - x1)^2 + (y2 - y1)^2)
    flip = root in (:top, :right)
    need_loop = (flip && y1 <= y2) || (!flip && y1 >= y2)

    # these points give the initial/final "rise"
    # note: this is a function of distance between points and axis scale
    y_offset = if need_loop
        0.3dist
    else
        min(0.3dist, 0.5*abs(y2-y1))
    end
    y_offset = max(0.02*(maxy-miny), y_offset)

    if flip
        # got the other direction
        y_offset *= -1
    end
    push!(y, y1 + y_offset)

    # try to figure out when to loop around vs just connecting straight
    if need_loop
        if abs(x2-x1) > 0.1 * (maxx - minx)
            # go between
            sgn = x2 > x1 ? 1 : -1
            x_offset = 0.5 * abs(x2-x1)
            append!(x, [x1 + sgn * x_offset, x2 - sgn * x_offset])
        else
            # add curve points which will create a loop
            x_offset = 0.3 * (maxx - minx) * (rand(Bool) ? 1 : -1)
            append!(x, [x1 + x_offset, x2 + x_offset])
        end
        append!(y, [y1 + y_offset, y2 - y_offset])
    end

    append!(x, [x2, x2])
    append!(y, [y2 - y_offset, y2])
    if root in (:left,:right)
        # flip x/y to simplify
        x,y = y,x
    end
    x, y
end

function shorten_segment(x1, y1, x2, y2, shorten)
    xshort = shorten * (x2-x1)
    yshort = shorten * (y2-y1)
    x1+xshort, y1+yshort, x2-xshort, y2-yshort
end

# """
#     shorten_segment_absolute(x1, y1, x2, y2, shorten)
#
# Remove an amount `shorten` from the end of the line [x1,y1] -> [x2,y2].
# """
# function shorten_segment_absolute(x1, y1, x2, y2, shorten)
#     if x1 == x2 && y1 == y2
#         return x1, y1, x2, y2
#     end
#     t = shorten/sqrt(x1*(x1-2x2) + x2^2 + y1*(y1-2y2) + y2^2)
#     x1, y1, (1.0-t)*x2 + t*x1, (1.0-t)*y2 + t*y1
# end

"""
    nearest_intersection(xs, ys, xd, yd, vec_xy_d)

Find where the line defined by [xs,ys] -> [xd,yd] intersects with the closed shape who's
vertices are stored in `vec_xy_d`. Return the intersection that is closest to the point
[xs,ys] (the source node).
"""
function nearest_intersection(xs, ys, xd, yd, vec_xy_d)
    if xs == xd && ys == yd
        return xs, ys, xd, yd
    end
    t = Vector{Float64}(undef, 2)
    xvec = Vector{Float64}(undef, 2)
    yvec = Vector{Float64}(undef, 2)
    xy_d_edge = Vector{Float64}(undef, 2)
    ret = Vector{Float64}(undef, 2)
    A = Array{Float64}(undef, 2, 2)
    nearest = Inf
    for i in 1:length(vec_xy_d)-1
        xvec .= [vec_xy_d[i][1], vec_xy_d[i+1][1]]
        yvec .= [vec_xy_d[i][2], vec_xy_d[i+1][2]]
        A .= [-xs+xd -xvec[1]+xvec[2] ; -ys+yd -yvec[1]+yvec[2]]
        t .= (A + eps()*I)\[xs-xvec[1] ; ys-yvec[1]]
        xy_d_edge .= [(1-t[2])*xvec[1] + t[2]*xvec[2], (1-t[2])*yvec[1] + t[2]*yvec[2]]
        if 0 <= t[2] <= 1
            tmp = abs2(xy_d_edge[1] - xs) + abs2(xy_d_edge[2] - ys)
            if tmp < nearest
                ret .= xy_d_edge
                nearest = tmp
            end
        end
    end
    xs, ys, ret[1], ret[2]
end

function nearest_intersection(xs, ys, xd, yd, vec_xy_d::GeometryTypes.Circle)
    if xs == xd && ys == yd
        return xs, ys, xd, yd
    end

    α = atan(ys - yd, xs - xd)
    xd = xd + vec_xy_d.r*cos(α)
    yd = yd + vec_xy_d.r*sin(α)

    xs, ys, xd, yd
end

function nearest_intersection(xs, ys, zs, xd, yd, zd, vec_xyz_d)
    # TODO make 3d work.
end

"""
Randomly pick a point to be the center control point of a bezier curve,
which is both equidistant between the endpoints and normally distributed
around the midpoint.
"""
function random_control_point(xi, xj, yi, yj, curvature_scalar)
    xmid = 0.5 * (xi+xj)
    ymid = 0.5 * (yi+yj)

    # get the angle of y relative to x
    theta = atan((yj-yi) / (xj-xi)) + 0.5pi

    # calc random shift relative to dist between x and y
    dist = sqrt((xj-xi)^2 + (yj-yi)^2)
    dist_from_mid = curvature_scalar * (rand()-0.5) * dist

    # now we have polar coords, we can compute the position, adding to the midpoint
    (xmid + dist_from_mid * cos(theta),
     ymid + dist_from_mid * sin(theta))
end

function control_point(xi, xj, yi, yj, dist_from_mid)
    xmid = 0.5 * (xi+xj)
    ymid = 0.5 * (yi+yj)

    # get the angle of y relative to x
    theta = atan((yj-yi) / (xj-xi)) + 0.5pi

    # dist = sqrt((xj-xi)^2 + (yj-yi)^2)
    # dist_from_mid = curvature_scalar * 0.5dist

    # now we have polar coords, we can compute the position, adding to the midpoint
    (xmid + dist_from_mid * cos(theta),
     ymid + dist_from_mid * sin(theta))
end

function annotation_extent(p, annotation; width_scalar=0.06, height_scalar=0.096)
    str = string(annotation[3])
    position = annotation[1:2]
    plot_size = get(p, :size, (600, 400))
    fontsize = annotation[4]
    xextent_length = width_scalar*(600/plot_size[1])*fontsize*length(str)^0.8
    xextent = [position[1]-xextent_length,position[1]+xextent_length]
    yextent_length = height_scalar*(400/plot_size[2])*fontsize
    yextent = [position[2] - yextent_length, position[2] + yextent_length]

    [xextent, yextent]
end

function clockwise_difference(angle1, angle2)
    pi - abs(abs(angle1 - angle2) - pi)
end

function clockwise_mean(angles)
    if clockwise_difference(angles[2], angles[1]) > angles[2] - angles[1]
        return mean(angles) + pi
    else
        return mean(angles)
    end
end

"""
    unoccupied_angle(x1, y1, x, y)

Starting from the point [x1,y1], find the angle theta such that a line leaving at an angle
theta will have maximum distance from the points [x[i],y[i]]
"""
function unoccupied_angle(x1, y1, x, y)
    @assert length(x) == length(y)

    if length(x) == 1
        return atan(y[1] - y1, x[1] - x1) + pi
    end

    max_range = zeros(2)
    # Calculate all angles between the point [x1,y1] and all points [x[i],y[i]], make sure
    # that all of the angles are between 0 and 2pi
    angles = [atan(y[i] - y1, x[i] - x1) for i in 1:length(x)]
    for i in 1:length(angles)
        if angles[i] < 0
            angles[i] += 2pi
        end
    end
    # Sort all of the angles and calculate which two angles subtend the largest gap.
    sort!(angles)
    max_range .= [angles[end], angles[1]]
    for i in 2:length(x)
        if (clockwise_difference(angles[i], angles[i-1])
            > clockwise_difference(max_range[2], max_range[1]))
            max_range .= [angles[i-1], angles[i]]
        end
    end
    # Return the angle that is in the middle of the two angles subtending the largest
    # empty angle.
    clockwise_mean(max_range)
end

function process_edge_attribute(attr, source, destiny, weights)
    if isnothing(attr) || (attr isa Symbol)
        return attr
    elseif attr isa LightGraphs.AbstractGraph
        mat = incidence_matrix(attr)
        attr = [mat[si, di] for (si, di) in zip(source, destiny)][:] |> permutedims
    elseif attr isa Function
        attr = [attr(si, di, wi) for (i, (si, di, wi)) in
                enumerate(zip(source, destiny, weights))][:] |> permutedims
    elseif attr isa Dict
        attr = [attr[(si, di)] for (si, di) in zip(source, destiny)][:] |> permutedims
    elseif all(size(attr) .!= 1)
        attr = [attr[si, di] for (si, di) in zip(source, destiny)][:] |> permutedims
    end
    attr
end
# Function from Plots/src/components.jl
"get an array of tuples of points on a circle with radius `r`"
function partialcircle(start_θ, end_θ, n = 20, r = 1)
    Tuple{Float64,Float64}[(r*cos(u), r*sin(u)) for u in
                           range(start_θ, stop = end_θ, length = n)]
end

function partialcircle(start_θ, end_θ, circle_center::Array{T,1}, n = 20, r = 1) where T
    Tuple{Float64,Float64}[(r*cos(u) + circle_center[1], r*sin(u) + circle_center[2]) for u in
                           range(start_θ, stop = end_θ, length = n)]
end

function partialellipse(start_θ, end_θ, n = 20, major_axis = 2, minor_axis = 1)
    Tuple{Float64,Float64}[(major_axis*cos(u), minor_axis*sin(u)) for u in
                           range(start_θ, stop = end_θ, length = n)]
end

function partialellipse(start_θ, end_θ, ellipse_center::Array{T,1}, n = 20, major_axis = 2,
                        minor_axis = 1) where T
    Tuple{Float64,Float64}[(major_axis*cos(u) + ellipse_center[1],
                            minor_axis*sin(u) + ellipse_center[2]) for u in
                           range(start_θ, stop = end_θ, length = n)]
end

# for chord diagrams:
function arcshape(θ1, θ2)
    vcat(
        partialcircle(θ1, θ2, 15, 1.05),
        reverse(partialcircle(θ1, θ2, 15, 0.95))
    )
end

# x and y limits for arc diagram ()
function arcdiagram_limits(x, source, destiny)
    @assert length(x) >= 2
    margin = abs(0.1*(x[2] - x[1]))
    xmin, xmax = extrema(x)
    r = abs(0.5 * (xmax - xmin))
    mean_upside = mean(source .< destiny)
    ylims = if mean_upside == 1.0
        (-margin, r+margin)
    elseif mean_upside == 0.0
        (-r-margin, margin)
    else
        (-r-margin, r+margin)
    end
    (xmin-margin, xmax+margin), ylims
end

function islabel(item)
    ismissing(item) && return false
    ((item isa AbstractFloat) && isnan(item)) && return false
    !in(item, (nothing, false, ""))
end

function replacement_kwarg(sym, name, plotattributes, graph_aliases)
    replacement = name
    for alias in graph_aliases[sym]
        if haskey(plotattributes, alias)
            replacement = plotattributes[alias]
        end
    end
    replacement
end

macro process_aliases(plotattributes, graph_aliases)
    ex = Expr(:block)
    attributes = getfield(__module__, graph_aliases) |> keys
    ex.args = [Expr(:(=), esc(sym), :($(esc(replacement_kwarg))($(QuoteNode(sym)), $(esc(sym)), $(esc(plotattributes)), $(esc(graph_aliases))))) for sym in attributes]
    ex
end

function remove_aliases!(sym, plotattributes, graph_aliases)
    for alias in graph_aliases[sym]
        if haskey(plotattributes, alias)
            delete!(plotattributes, alias)
        end
    end
end

# From Plots/src/utils.jl
isnothing(x::Nothing) = true
isnothing(x) = false

# From Plots/src/Plots.jl
ignorenan_extrema(x) = Base.extrema(x)
# From Plots/src/utils.jl
ignorenan_extrema(x::AbstractArray{F}) where {F<:AbstractFloat} = NaNMath.extrema(x)
# From Plots/src/components.jl
function extrema_plus_buffer(v, buffmult = 0.2)
    vmin,vmax = extrema(v)
    vdiff = vmax-vmin
    zero_buffer = vdiff == 0 ? 1.0 : 0.0
    buffer = (vdiff+zero_buffer) * buffmult
    vmin - buffer, vmax + buffer
end


text:
abstract type A

end

# include("temp2.jl")

type B <: A
    x
    y
end

type C <: A
    x
    z
end

prop1(b::B) = b.x
prop2(b::B) = b.y
prop1(c::C) = c.w
prop2(c::C) = c.z  # changed from prop2(c::C)=c.w

mysum(a::A) = prop1(a) + prop2(a)
function getx(a::A)
  b = a.x
  return b
end


text:
# This file was generated, do not modify it. # hide
lgb = LGBMRegressor() #initialised a model with default params
lgbm = machine(lgb, features[train, :], targets[train, 1])
boostrange = range(lgb, :num_iterations, lower=2, upper=500)
curve = learning_curve!(lgbm, resampling=CV(nfolds=5),
                        range=boostrange, resolution=100,
                        measure=rms)

figure(figsize=(8,6))
plot(curve.parameter_values, curve.measurements)
xlabel("Number of rounds", fontsize=14)
ylabel("RMSE", fontsize=14)

plt.savefig(joinpath(@OUTPUT, "lgbm_hp1.svg")) # hide


text:
# More generic implementation.
# https://github.com/ChrisRackauckas/Lists.jl/blob/master/src/Lists.jl
#
abstract type AbstractData end
abstract type AbstractNode end
struct NullNode <: AbstractNode end
const null = NullNode()

mutable struct Data <: AbstractData
    t::Int
    v::Float64
end

mutable struct Node{P, N} <: AbstractNode
    prev:: P
    next:: N
    data:: AbstractData
    function Node{P,N}(p,n,d) where {P,N}
        new(p,n,d)
    end
end
# Node(d::AbstractData) = Node{AbstractNode, AbstractNode}(null, null, d)
Node(d::AbstractData) = Node{AbstractNode, AbstractNode}(null, null, d)
Node(p::AbstractNode, n::AbstractNode, d::AbstractData) = Node{AbstractNode, AbstractNode}(p,n,d)

mutable struct List{T<:AbstractNode}
    # https://github.com/JuliaLang/julia/issues/20877
    # https://docs.julialang.org/en/latest/manual/constructors.html#Case-Study:-Rational-1
    head::T
    tail::T
    # function List{T}() where T <: AbstractNode
    #     new(null, null)
    # end
    function List{T}(nditem::T) where T <: AbstractNode
        new(nditem, nditem)
    end
end
### Start outer constructors
List() = List{NullNode}()
List(nd::AbstractNode) = List{AbstractNode}(nd)
### End outer constructors

function add!(list::List, data::Data)
    node = Node{NullNode, AbstractNode}(null, list.head, data)
    list.head = node
    node.next.prev = node
    node.prev = null # --todo
    println("Added!")
end
# list = List(Node(Data(1,0.1)))
# add!(list, Data(1, 1.0))
# add!(list, Data(1, 1.0))
nd1=Node(Data(1,0.1))
nd2=Node(Data(2,0.2))
println(typeof(nd1.prev) <: NullNode)
nd1.prev = nd2
println(typeof(nd1.prev) <: NullNode)

## tests
# function test_node()
#     init_node = Node(Data(1, 0.1))
#     node1 = Node(init_node, null, Data(2, 0.1))
#     node2 = Node(node1, null, Data(3, 0.1))
#     node3 = Node(node2, null, Data(4, 0.1))
#     typeof(node2.prev) <: NullNode
#     typeof(node3.next) <: NullNode
# end
# test_node()

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
# Autogenerated wrapper script for CyrusSASL_jll for powerpc64le-linux-gnu
export libanonymous, libcrammd5, libdigestmd5, libntlm, libotp, libplan, libsasl2, libscram

using OpenSSL_jll
using SQLite_jll
using Gdbm_jll
JLLWrappers.@generate_wrapper_header("CyrusSASL")
JLLWrappers.@declare_library_product(libanonymous, "libanonymous.so.3")
JLLWrappers.@declare_library_product(libcrammd5, "libcrammd5.so.3")
JLLWrappers.@declare_library_product(libdigestmd5, "libdigestmd5.so.3")
JLLWrappers.@declare_library_product(libntlm, "libntlm.so.3")
JLLWrappers.@declare_library_product(libotp, "libotp.so.3")
JLLWrappers.@declare_library_product(libplan, "libplain.so.3")
JLLWrappers.@declare_library_product(libsasl2, "libsasl2.so.3")
JLLWrappers.@declare_library_product(libscram, "libscram.so.3")
function __init__()
    JLLWrappers.@generate_init_header(OpenSSL_jll, SQLite_jll, Gdbm_jll)
    JLLWrappers.@init_library_product(
        libanonymous,
        "lib/sasl2/libanonymous.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcrammd5,
        "lib/sasl2/libcrammd5.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdigestmd5,
        "lib/sasl2/libdigestmd5.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libntlm,
        "lib/sasl2/libntlm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libotp,
        "lib/sasl2/libotp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libplan,
        "lib/sasl2/libplain.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsasl2,
        "lib/libsasl2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libscram,
        "lib/sasl2/libscram.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct WatchState
    filename::String
    mtime ::Float64
    WatchState(filename) = new(filename, mtime(filename))
end
function poll_until_changed(files::Set{String}; interval = 1.)
    watched = WatchState[]
    for f in files
        !isfile(f) && return
        push!(watched, WatchState(f))
    end
    active = true
    while active
        for w in watched
            if !isfile(w.filename)
                active = false
                break;
            end
            if mtime(w.filename) != w.mtime
                active = false
                break;
            end
        end
        sleep(interval)
    end
end

function poll_folders(on_change, folders; interval = 1.)
    watched = Dict{String, Float64}()
    while true
        walked = Set{String}()
        for folder in folders
            !isdir(folder) && continue
            for (base, _, files) in walkdir(folder)
                for f in files
                    path = joinpath(base, f)
                    new_time = mtime(path)
                    if new_time > get(watched, path, -1.)
                        on_change(path, new_time, false)
                    end
                    watched[path] = new_time
                    push!(walked, path)
                end
            end
        end

        for path in keys(watched)
            if !(path in walked)
                on_change(path, -1., true)
                delete!(watched, path)
            end
        end
        sleep(interval)
    end
end