text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Data: 05/02/2014

using Util

function validate(aem)

    validate_initial(aem)
    #validate_transition(aem)
end

function validate_initial(aem)

    params = aem.parameters

    marginal_dist = Array{Any}(params.n_initial)
    marginal_prob = Array{Any}(params.n_initial)

    for i = 1:params.n_initial
        marginal_dist[i] = vec(sum(params.N_initial[i], 2))
        marginal_prob[i] = marginal_dist[i] / sum(marginal_dist[i])
    end

    f = open(aem.initial_sample_filename, "r")

    readline(f)

    n_nmac = 0
    n_lines = 0

    sample_dist = Array{Any}(params.n_initial)
    sample_prob = Array{Any}(params.n_initial)

    for i = 1:params.n_initial
        sample_dist[i] = zeros(Int, params.r_initial[i])
    end

    for line = eachline(f) #(ix, line) in enumerate(eachline(f))
        values = float(split(chomp(line)))

        # 1       id
        # 2       A = [1:4]
        # 3       L = [1:5]
        # 4       chi = [90, 270]   : deg
        # 5       beta = [0:30:360] : deg
        # 6       C1 = [1, 2]
        # 7       C2 = [1, 2]
        # 8,9     v1, v2 = [50, 100, 200, 300, 400, 500, 600]   : kt
        # 10, 11  v1d, v2d = [-5, -2, -0.25, 0.25, 2, 5]        : kt/s
        # 12, 13  h1d, h2d = [-5000, -3000, -2000, -1000, -400, 400, 1000, 2000, 3000, 5000]    : ft/min
        # 14, 15  psi1d, psi2d = [-8, -6, -3.5, -1, -0.25, 0.25, 1, 3.5, 6, 8]  : deg/s
        # 16      hmd = [0, 0.0822896, 0.5, 1, 3]   : nm
        # 17      vmd = [0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 6000]  : ft

        for i = [2, 3, 4, 6, 7]
            if !(Int(values[i]) in 1:params.r_initial[i-1])
                println("sample: ", chomp(line))
                println("i: ", i, ", value: ", Int(values[i]), ", boundary: ", params.r_initial[i-1])
                error("The value violates the boundary.")
            end
            #println("Updating sample_dist[", i-1, "][", Int(values[i]), "] += 1")
            sample_dist[i - 1][Int(values[i])] += 1
        end

        for i = [5, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]
            if !isempty(params.boundaries[i-1]) && (values[i] < params.boundaries[i-1][1] || values[i] > params.boundaries[i-1][end])
                println("sample: ", chomp(line))
                println("i: ", i, ", value: ", values[i], ", boundary: ", params.boundaries[i-1][end])
                error("The value violates the boundary.")
            end

            j = findfirst(x -> (x > values[i]), params.boundaries[i - 1])

            if j == 0 && values[i] == params.boundaries[i - 1][end]
                sample_dist[i - 1][params.r_initial[i - 1]] += 1
            else
                sample_dist[i - 1][j - 1] += 1
            end
        end

        hmd_ft = 6076.12 * values[16]
        vmd = values[17]

        if hmd_ft <= 500 && vmd <= 200
            n_nmac += 1
        end

        n_lines += 1
    end

    close(f)

    for i = 1:params.n_initial
        sample_prob[i] = sample_dist[i] / sum(sample_dist[i])
    end

    for i = 1:params.n_initial
        println("variable #: ", i)

        for j = 1:params.r_initial[i]
            print(marginal_dist[i][j], " ")
            print(outputGFormatString(marginal_prob[i][j]), " ")
            print(sample_dist[i][j], " ")
            print(outputGFormatString(sample_prob[i][j]), " ")
            @printf("%.2f%%\n", abs((sample_prob[i][j] - marginal_prob[i][j]) / marginal_prob[i][j] * 100))
        end

        println()
    end

    println("n_nmac: ", n_nmac, " n_lines: ", n_lines)
    println()
    println("P(NMAC | enc, no TCAS) = ", n_nmac / n_lines, " ($n_nmac / $n_lines)")
    println()
end

function validate_transition(aem)

    # should validate conditional probabilities rather than marginal probability
    # check marginal probability here for convenience

    params = aem.parameters

    marginal_dist = Array{Any}(params.n_transition - params.n_initial)
    marginal_prob = Array{Any}(params.n_transition - params.n_initial)

    for i = 1:(params.n_transition - params.n_initial)
        marginal_dist[i] = vec(sum(params.N_transition[params.n_initial + i], 2))
        marginal_prob[i] = marginal_dist[i] / sum(marginal_dist[i])
    end

    f = open(aem.transition_sample_filename, "r")

    readline(f)

    sample_dist = Array{Any}(params.n_transition - params.n_initial)
    sample_prob = Array{Any}(params.n_transition - params.n_initial)

    for i = 1:(params.n_transition - params.n_initial)
        sample_dist[i] = zeros(Int, params.r_transition[params.n_initial + i])
    end

    bin = Array(Int, params.n_transition - params.n_initial)
    value = Array(AbstractFloat, params.n_transition - params.n_initial)

    n_sample = 0
    n_resample = zeros(Int, params.n_transition - params.n_initial)

    for line = eachline(f)
        values = float(split(chomp(line)))

        # 1     id
        # 2     time step
        # 3, 4  h1d, h2d        : ft/min
        # 5, 6  psi1d, psi2d    : deg/s

        # resampling_rate = 0 0 0 0 0 0 0 0 0 0 0.0487462 0.0505306 0.0794427 0.0827686 0 0

        id = Int(values[1])
        ts = Int(values[2])

        if ts != 0
            n_sample += 1
        end

        for i = [3 4 5 6]
            if !isempty(params.boundaries[i + 8]) && (values[i] < params.boundaries[i + 8][1] || values[i] > params.boundaries[i + 8][end])
                println("sample: ", chomp(line))
                println("i: ", i, ", value: ", values[i], ", boundary: ", params.boundaries[i + 8][end])
                error("The value violates the boundary.")
            end

            j = findfirst(x -> (x > values[i]), params.boundaries[i + 8])

            if ts != 0
                if j == 0 && values[i] == params.boundaries[i + 8][end]
                    sample_dist[i - 2][params.r_transition[params.n_initial + i - 2]] += 1
                else
                    sample_dist[i - 2][j - 1] += 1
                end

                if (j - 1) == bin[i - 2] && values[i] != value[i - 2]
                    n_resample[i - 2] += 1
                end
            end

            bin[i - 2] = j - 1
            value[i - 2] = values[i]
        end
    end

    close(f)

    for i = 1:(params.n_transition - params.n_initial)
        sample_prob[i] = sample_dist[i] / sum(sample_dist[i])
    end

    for i = 1:(params.n_transition - params.n_initial)
        println("variable #: ", params.n_initial + i)

        for j = 1:params.r_transition[params.n_initial + i]
            print(marginal_dist[i][j], " ")
            print(outputGFormatString(marginal_prob[i][j]), " ")
            print(sample_dist[i][j], " ")
            print(outputGFormatString(sample_prob[i][j]), " ")
            @printf("%.2f%%\n", abs((sample_prob[i][j] - marginal_prob[i][j]) / marginal_prob[i][j] * 100))
        end
        println()

        print(outputGFormatString(params.resample_rates[i + 10]), " ", outputGFormatString(n_resample[i] / n_sample), " (", n_resample[i], " / ", n_sample, ") ")
        @printf("%.2f%%\n", abs(((n_resample[i] / n_sample) - params.resample_rates[i + 10]) / params.resample_rates[i + 10] * 100))
        println()
    end
end


text:
@doc raw"""
    MinimumError <: AbstractImageBinarizationAlgorithm
    MinimumError()

    binarize([T,] img, f::MinimumError)
    binarize!([out,] img, f::MinimumError)

Under the assumption that the image histogram is a mixture of two Gaussian
distributions the binarization threshold is chosen such that the expected
misclassification error rate is minimised.

# Output

Return the binarized image as an `Array{Gray{T}}` of size `size(img)`. If
`T` is not specified, it is inferred from `out` and `img`.

# Details

Let ``f_i`` ``(i=1 \ldots I)`` denote the number of observations in the
``i``th bin of the histogram. Then the probability that an observation
belongs to the ``i``th bin is given by  ``p_i = \frac{f_i}{N}`` (``i = 1,
\ldots, I``), where ``N = \sum_{i=1}^{I}f_i``.

The minimum error thresholding method assumes that one can find a threshold
``T`` which partitions the data into two categories,  ``C_0`` and ``C_1``, such that
the data can be modelled by a mixture of two Gaussian distribution. Let
```math
P_0(T) = \sum_{i = 1}^T p_i \quad \text{and} \quad P_1(T) = \sum_{i = T+1}^I p_i
```
denote the cumulative probabilities,
```math
\mu_0(T) = \sum_{i = 1}^T i \frac{p_i}{P_0(T)} \quad \text{and} \quad \mu_1(T) = \sum_{i = T+1}^I i \frac{p_i}{P_1(T)}
```
denote the means, and
```math
\sigma_0^2(T) = \sum_{i = 1}^T (i-\mu_0(T))^2 \frac{p_i}{P_0(T)} \quad \text{and} \quad \sigma_1^2(T) = \sum_{i = T+1}^I (i-\mu_1(T))^2 \frac{p_i}{P_1(T)}
```
denote the variances of categories ``C_0`` and ``C_1``, respectively.

Kittler and Illingworth proposed to use the minimum error criterion function
```math
J(T) = 1 + 2 \left[ P_0(T) \ln \sigma_0(T) + P_1(T) \ln \sigma_1(T) \right] - 2 \left[P_0(T) \ln P_0(T) + P_1(T) \ln P_1(T) \right]
```
to assess the discreprancy between the mixture of Gaussians implied by a particular threshold ``T``,
and the piecewise-constant probability density function represented by the histogram.
The discrete value ``T`` which minimizes the function ``J(T)`` produces
the sought-after threshold value (i.e. the bin which determines the threshold).

# Arguments

The function argument is described in more detail below.

##  `img::AbstractArray`

The image that needs to be binarized. The image is automatically converted
to `Gray` in order to construct the requisite graylevel histogram.

# Example

Binarize the "cameraman" image in the `TestImages` package.

```julia
using TestImages, ImageBinarization

img = testimage("cameraman")
img_binary = binarize(img, MinimumError())
```

# References

1. J. Kittler and J. Illingworth, “Minimum error thresholding,” Pattern Recognition, vol. 19, no. 1, pp. 41–47, Jan. 1986. [doi:10.1016/0031-3203(86)90030-0](https://doi.org/10.1016/0031-3203%2886%2990030-0)
2. Q.-Z. Ye and P.-E. Danielsson, “On minimum error thresholding and its implementations,” Pattern Recognition Letters, vol. 7, no. 4, pp. 201–206, Apr. 1988. [doi:10.1016/0167-8655(88)90103-1](https://doi.org/10.1016/0167-8655%2888%2990103-1)
"""
struct MinimumError <: AbstractImageBinarizationAlgorithm end

function (f::MinimumError)(out::GenericGrayImage, img::GenericGrayImage)
    edges, counts = build_histogram(img,  256)
    t = find_threshold(HistogramThresholding.MinimumError(), counts[1:end], edges)
    @simd for i in CartesianIndices(img)
        out[i] = img[i] < t ? 0 : 1
    end
    out
end

(f::MinimumError)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))

ext: jl


text:
using Plots
using Printf

"""
    untildivergence(f, x, maxiter=50)

Iterate `f` on `x` until its module is greater than 2 or `maxiter` iterations have been done.

Returns the number of iterations.
"""
function untildivergence(f, x, maxiter = 50)
    i = 0
    while abs(x) < 2 && (i += 1) <= maxiter
        x = f(x)
    end
    return i
end

"""
    julia_set_wh(f; kwargs)

Approximate the Julia set of `f` in a grid using the escape time algorithm.

The grid is defined by its bounds `xmin`, `xmax`, `ymin` and `ymax` as well as its width `w` and height `h` in pixels. The optional parameter `maxiter` is passed on to `untildivergence`.
"""
function julia_set_wh(f; xmin::Real=-2, xmax::Real=2, ymin::Real=-2, ymax::Real=2, w::Int=1000, h::Int=1000, maxiter::Int=50)
    img = Array{Int,2}(undef, h, w)
    x = LinRange(xmin, xmax, w)
    y = LinRange(ymin, ymax, h)
    # creating labels for use with `heatmap`
    xticks = LinRange(1,w,11)
    yticks = LinRange(1,h,11)
    xlabels = map(z->@sprintf("%.2f",z), LinRange(xmin,xmax,11))
    ylabels = map(z->@sprintf("%.2f",z), LinRange(ymin,ymax,11))

    for i in 1:w, j in 1:h
        img[j,i] = untildivergence(f, x[i]+im*y[j], maxiter)
    end
    return img, xticks, xlabels, yticks, ylabels
end

"""
    julia_set(f, cellsize; kwargs)

Approximate the Julia of `f` set in a grid of cells of length `cellsize`.

The grid is defined by the same parameters as in `julia_set_wh`, except its width and height in pixels are computed with `cellsize`.
"""
function julia_set(f, cellsize::Real; xmin=-2, xmax=2, ymin=-2, ymax=2, maxiter=50)
    w = Integer(round(abs(xmax-xmin)/cellsize))
    h = Integer(round(abs(ymax-ymin)/cellsize))
    return julia_set_wh(f, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, w=w, h=h, maxiter=maxiter)
end

"""
    mandelbrot_set_wh(f; kwargs)

Approximate the Mandelbrot set of `f` in a grid using the escape time algorithm.

The grid is defined by its bounds `xmin`, `xmax`, `ymin` and `ymax` as well as its width `w` and height `h` in pixels. The optional parameter `maxiter` is passed on to `untildivergence`.
"""
function mandelbrot_set_wh(f; xmin::Real=-2, xmax::Real=2, ymin::Real=-2, ymax::Real=2, w::Integer=1000, h::Integer=1000, maxiter::Integer=50)
    img = Array{Int,2}(undef, h, w)
    x = LinRange(xmin, xmax, w)
    y = LinRange(ymin, ymax, h)
    # creating labels for use with `heatmap`
    xticks = LinRange(1,w,11)
    yticks = LinRange(1,h,11)
    xlabels = map(z->@sprintf("%.2f",z), LinRange(xmin,xmax,11))
    ylabels = map(z->@sprintf("%.2f",z), LinRange(ymin,ymax,11))

    for i in 1:w, j in 1:h
        img[j,i] = untildivergence(f(x[i]+im*y[j]), 0, maxiter)
    end
    return img, xticks, xlabels, yticks, ylabels
end

"""
    mandelbrot_set(f, cellsize; kwargs)

Approximate the Mandelbrot of `f` set in a grid of cells of length `cellsize`.

The grid is defined by the same parameters as in `julia_set_wh`, except its width and height in pixels are computed with `cellsize`.
"""
function mandelbrot_set(f, cellsize; xmin=-2, xmax=2, ymin=-2, ymax=2, maxiter=50)
    w = Integer(round(abs(xmax-xmin)/cellsize))
    h = Integer(round(abs(ymax-ymin)/cellsize))
    return mandelbrot_set_wh(f, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, w=w, h=h, maxiter=maxiter)
end

# the traditional julia/mandelbrot function
f(c) = z -> z^2 + c

img, xticks, xlabels, yticks, ylabels = mandelbrot_set(f, 0.00005, xmin=-2.01, xmax = -1.8, ymin = -0.1, ymax = 0.1)

img, xticks, xlabels, yticks, ylabels = julia_set(f(im), 0.0005, xmin=-2, xmax=2, ymin=-2,ymax=2)

heatmap(img,
        xticks = (xticks, xlabels),
        yticks = (yticks, ylabels),
        aspect_ratio=:equal
        )

# alternative, for better looking images
using Images, Colors
g = Gray.(img/51)

g = Gray.(.!(img .>= 50))

mset = img .>= 50

#save("julia_minus2.png", g)

include("FractalDimensions.jl")

B = boundary(mset)
Gray.(B)
points = array_to_points(B)

d, N, ε = box_counting_dimension(points, 10)

C, r = correlation_dimension(points)

loglog_regression(C,r)

max_stars_repo_name: csimal/chaos-and-fractals

max_issues_repo_name: csimal/chaos-and-fractals

max_forks_repo_name: csimal/chaos-and-fractals


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

export script_base, script_dasc, script_libcas098small, script_libcas0100star, script_libcas0100llcem
export process_jsons

using FilterNMACInfo

const DATADIR = joinpath(dirname(@__FILE__), "..", "..", "..", "data")

const DASC_JSON = joinpath(DATADIR, "dasc/json")
const DASC_CSV = joinpath(DATADIR, "dasc/csv")
const DASC_OUT = Pkg.dir("Datasets/data/dasc") #requires Datasets to be installed

const LIBCAS098SMALL_JSON = joinpath(DATADIR, "libcas098small/json")
const LIBCAS098SMALL_CSV = joinpath(DATADIR, "libcas098small/csv")
const LIBCAS098SMALL_OUT = Pkg.dir("Datasets/data/libcas098small")

const LIBCAS0100STAR_JSON = joinpath(DATADIR, "libcas0100star/json")
const LIBCAS0100STAR_CSV = joinpath(DATADIR, "libcas0100star/csv")
const LIBCAS0100STAR_OUT = Pkg.dir("Datasets/data/libcas0100star")

const LIBCAS0100LLCEM_JSON = joinpath(DATADIR, "libcas0100llcem/json")
const LIBCAS0100LLCEM_CSV = joinpath(DATADIR, "libcas0100llcem/csv")
const LIBCAS0100LLCEM_OUT = Pkg.dir("Datasets/data/libcas0100llcem")

#####################
#Add an entry here

#dasc set
function script_dasc(fromjson::Bool=true)
  script_base(DASC_JSON, DASC_CSV, DASC_OUT;
                 fromjson=fromjson, correct_coc=true)
end

#from APL 20151230, libcas0.9.8, MCTS iterations=500, testbatch
function script_libcas098small(fromjson::Bool=true)
  script_base(LIBCAS098SMALL_JSON, LIBCAS098SMALL_CSV, LIBCAS098SMALL_OUT;
                 fromjson=fromjson, correct_coc=true)
end

#Generated 20160413, libcas0.10.0, MCTS iterations=3000, 2ac, stardbn
function script_libcas0100star(fromjson::Bool=true)
  script_base(LIBCAS0100STAR_JSON, LIBCAS0100STAR_CSV, LIBCAS0100STAR_OUT;
                 fromjson=fromjson, correct_coc=false)
end

#Generated 20160422, libcas0.10.0, MCTS iterations=3000, 2ac, llcemdbn
function script_libcas0100llcem(fromjson::Bool=true)
  script_base(LIBCAS0100LLCEM_JSON, LIBCAS0100LLCEM_CSV, LIBCAS0100LLCEM_OUT;
                 fromjson=fromjson, correct_coc=false)
end

#####################
function script_base(jsondir::AbstractString, csvdir::AbstractString,
                        outdir::AbstractString;
                        fromjson::Bool=true, correct_coc::Bool=true,
                        verbose::Bool=true)
  if fromjson
    verbose && println("converting to csv...")
    mkpath(csvdir)
    convert2csvs(jsondir, csvdir)
  end
  tmpdir = mktempdir()
  verbose && println("tmp=$tmpdir")
  verbose && println("converting to dataframes...")
  csvs2dataframes(csvdir, tmpdir)
  if correct_coc
    verbose && println("correcting cocs...")
    correct_coc_stays!(tmpdir)
  end
  add_encounter_info!(tmpdir)

  mkpath(outdir)
  verbose && println("saving dataset...")
  make_dataset(tmpdir, jsondir, outdir)
end

"""
Convenience function to process a new dataset from directory of json files
dataname is the name of the dataset (no spaces)
jsonpath is the directory of input json files
adds 2 new datasets to Datasets/data (original and filtered NMAC)
"""
function process_jsons(dataname::AbstractString, jsonpath::AbstractString)
    csvpath = joinpath(jsonpath, "csv")
    outpath = Pkg.dir("Datasets", "data", dataname)
    script_base(jsonpath, csvpath, outpath; fromjson=true, correct_coc=false)

    filtname = dataname * "filt"
    filtpath = Pkg.dir("Datasets", "data", filtname)
    filter_nmac_info(FilterNMACInfo.isnmac, dataname, filtpath)
end


text:
abstract type SystemCategory <: PowerSystemCaseBuilderType end
struct PSYTestSystems <: SystemCategory end
struct PSITestSystems <: SystemCategory end
struct SIIPExampleSystems <: SystemCategory end
struct PSIDTestSystems <: SystemCategory end
struct PSSETestSystems <: SystemCategory end
struct MatPowerTestSystems <: SystemCategory end

const PACKAGE_DIR = joinpath(dirname(dirname(pathof(PowerSystemCaseBuilder))))

const SYSTEM_DESCRIPTORS_FILE = joinpath(PACKAGE_DIR, "src", "system_descriptor.jl")

const SERIALIZED_DIR = joinpath(PACKAGE_DIR, "data", "serialized_system")
const SERIALIZE_NOARGS_DIR = joinpath(PACKAGE_DIR, "data", "serialized_system", "NoArgs")
const SERIALIZE_FORECASTONLY_DIR =
    joinpath(PACKAGE_DIR, "data", "serialized_system", "ForecastOnly")
const SERIALIZE_RESERVEONLY_DIR =
    joinpath(PACKAGE_DIR, "data", "serialized_system", "ReserveOnly")
const SERIALIZE_FORECASTRESERVE_DIR =
    joinpath(PACKAGE_DIR, "data", "serialized_system", "ForecastReserve")

const SEARCH_DIRS = [
    SERIALIZE_NOARGS_DIR,
    SERIALIZE_FORECASTONLY_DIR,
    SERIALIZE_RESERVEONLY_DIR,
    SERIALIZE_FORECASTRESERVE_DIR,
]

const SERIALIZE_FILE_EXTENSIONS =
    [".json", "_validation_descriptors.json", "_time_series_storage.h5"]

function download_RTS(; kwargs...)
    PowerSystems.download("https://github.com/GridMod/RTS-GMLC", "master", "./data")
end


text:
"""
    cores(g)

Returns a vector `deg` such that if `deg[v]=k` then the vertex `v`
belongs to the `k`-core of `g` and not to the `k+1`-core.

See also [`kcore`](@ref).
"""
function cores(g::AGraph)
    n = nv(g)
    deg = degree(g)
    md = maximum(deg)
    bin = zeros(Int, md+1)
    for v=1:n
        bin[deg[v]+1] += 1
    end

    start = 1;
    for d=1:md+1
        num = bin[d]
        bin[d] = start
        start += num
    end
    #now bin[d+1] contains the position in vert of the first v with degree d

    #sort array vert according to degree and track their position in pos
    pos = zeros(Int, n)
    vert = zeros(Int, n)
    for v=1:n
        pos[v] = bin[deg[v]+1]
        vert[pos[v]] = v
        bin[deg[v]+1] += 1
    end

    # bring back bin to stating positions
    for d=md+1:-1:2
        bin[d] = bin[d-1]
    end
    bin[1] = 1

    #the main part of the algorithm
    for i=1:n
        v = vert[i]
        for u in neighbors(g, v)
            if deg[u] > deg[v] #move u one bin left
                du = deg[u]
                pu = pos[u]
                pw = bin[du+1]
                w = vert[pw]
                if u != w
                    pos[u], pos[w] = pw, pu
                    vert[pu], vert[pw] = w, u
                end
                bin[du+1] += 1
                deg[u] -= 1
            end
        end
    end
    return deg
end

"""
    kcore(g, k) -> (gnew, vmap)

Returns the `k`-core  of `g` along with a vertex map associating the mutated vertex
indexes to the old ones (as in [`rem_vertices!`](@ref)).

See also [`cores`](@ref)
"""
function kcore(g::AGraph, k::Integer)
    #TODO use subgraph
    gnew = copy(g)
    vmap = rem_vertices!(gnew, findall(d -> d < k, cores(g)))
    return gnew, vmap
end


text:
mutable struct example_t <: LCMType
    timestamp::Int64
    position::SVector{3, Float64}
    orientation::SVector{4, Float64}
    num_ranges::Int32
    ranges::Vector{Int16}
    name::String
    enabled::Bool
end

@lcmtypesetup(example_t,
  ranges => (num_ranges,),
)


text:
#############################################################################
#############################################################################
#
# This file contains units tests for polynomial factorization
#
#############################################################################
#############################################################################

"""
Test factorization of polynomials.
"""
function factor_test_poly(;N::Int = 25, seed::Int = 0, primes::Vector{Int} = [3,5,7,17,19])
    Random.seed!(seed)
    for prime in primes
        print("\ndoing prime = $prime \t")
        for _ in 1:N
            print(".")
            p = rand(Polynomial)
            factorization = factor(p, prime)
            pr = mod(expand_factorization(factorization),prime)
            @assert mod(p-pr,prime) == 0
        end
    end

    println("\nfactor_test_poly - PASSED")
end

"""
Test factorization of polynomials modulo p.
"""
function factor_test_poly_mod_p(;N::Int = 10, seed::Int = 0, primes::Vector{Int} = primes(3,25))
    Random.seed!(seed)
    x = x_poly()
    for prime in primes
        print("\ndoing prime = $prime \t")
        for _ in 1:N
            print(".")
            # Limiting to a polynomial of degree 6 since the ddfactor function scales poorly with an increase in the degree of the input polynomial
            p = rand(Polynomial, degree = rand(1:6))
            pmod = PolynomialModP(p, prime)
            factorization = factor(pmod)
            pr = mod(expand_factorization(factorization),prime)
            @assert mod(pmod.terms-pr,prime) == 0
        end
    end

    println("\nfactor_test_poly_mod_p - PASSED")
end


text:
function binforlength(len,bins,max_island_length)
    return floor(Int, ((len+eps())/max_island_length) * bins)
end

function count_islands_and_append_to_array!(genome_length, g::Vector
        ; max_no_of_reads::Int,read_length::Int)
        # returns arrays of length 100 , with maximum possible length of islands of a particular
        # read length mapped to be 100 and least possible length of a island ie. read length of a single island mapped to 1.

        data=[]
        number_of_reads = length(g)
        no_islands=zeros(max_no_of_reads)
        cum_length_islands=zeros(max_no_of_reads)
        hist_island_dist=[zeros(100) for i in 1:max_no_of_reads]
        for (r, len) in IslandIterator(g, read_length)
            if r<=max_no_of_reads
                no_islands[r]+=1
                cum_length_islands[r]+=len
                hist_island_dist[r][binforlength(len,100,r*read_length)]+=1.0
            end

        end

        for i in 1:max_no_of_reads
                push!(data, (genome_length,number_of_reads,i, no_islands[i], cum_length_islands[i],hist_island_dist[i]))
        end

        return data

end


text:
##############
# BlockSizes #
##############

abstract type AbstractBlockSizes{N} end

# Keeps track of the (cumulative) sizes of all the blocks in the `BlockArray`.
struct BlockSizes{N,VT<:NTuple{N,AbstractVector{Int}}} <: AbstractBlockSizes{N}
    cumul_sizes::VT
    # Takes a tuple of sizes, accumulates them and create a `BlockSizes`
    BlockSizes{N,VT}() where {N,VT<:NTuple{N,AbstractVector{Int}}} = new{N,VT}()
    BlockSizes{N,VT}(cs::VT) where {N,VT<:NTuple{N,AbstractVector{Int}}} = new{N,VT}(cs)
end

const DefaultBlockSizes{N} = BlockSizes{N,NTuple{N,Vector{Int}}}

BlockSizes{N}(cs::VT) where {N,VT<:NTuple{N,AbstractVector{Int}}} = BlockSizes{N,VT}(cs)
BlockSizes{N}() where N = BlockSizes{N,NTuple{N,Vector{Int}}}()
BlockSizes() = BlockSizes{0}()

BlockSizes(cs::NTuple{N,AbstractVector{Int}}) where N = BlockSizes{N}(cs)

function BlockSizes(sizes::Vararg{AbstractVector{Int}, N}) where {N}
    cumul_sizes = ntuple(k -> _cumul_vec(sizes[k]), Val(N))
    return BlockSizes(cumul_sizes)
end

Base.:(==)(a::BlockSizes, b::BlockSizes) = cumulsizes(a) == cumulsizes(b)

Base.dataids(b::BlockSizes) = _splatmap(dataids, b.cumul_sizes)
# _splatmap taken from Base:
_splatmap(f, ::Tuple{}) = ()
_splatmap(f, t::Tuple) = (f(t[1])..., _splatmap(f, tail(t))...)

function _cumul_vec(v::AbstractVector{T}) where {T}
    v_cumul = similar(v, length(v) + 1)
    z = one(T)
    v_cumul[1] = z
    for i in eachindex(v)
        z += v[i]
        v_cumul[i+1] = z
    end
    return v_cumul
end

@propagate_inbounds cumulsizes(block_sizes::BlockSizes) = block_sizes.cumul_sizes
@propagate_inbounds cumulsizes(block_sizes::AbstractBlockSizes, i) = cumulsizes(block_sizes)[i]
@propagate_inbounds cumulsizes(block_sizes::AbstractBlockSizes, i, j) = cumulsizes(block_sizes,i)[j]

@propagate_inbounds blocksize(block_sizes::AbstractBlockSizes, i, j) =
    cumulsizes(block_sizes, i, j+1) - cumulsizes(block_sizes, i, j)

# ntuple with Val was slow here. @generated it is!
@generated function blocksize(block_sizes::AbstractBlockSizes{N}, i::NTuple{N, Integer}) where {N}
    exp = Expr(:tuple, [:(blocksize(block_sizes, $k, i[$k])) for k in 1:N]...)
    return exp
end

# Gives the total sizes
@generated function Base.size(block_sizes::AbstractBlockSizes{N}) where {N}
    exp = Expr(:tuple, [:(cumulsizes(block_sizes, $i)[end] - 1) for i in 1:N]...)
    return quote
        @inbounds return $exp
    end
end

@inline function Base.size(block_sizes::AbstractBlockSizes{1})
    (cumulsizes(block_sizes,1)[end] - 1,)
end

@inline function Base.size(block_sizes::AbstractBlockSizes{2})
    (cumulsizes(block_sizes,1)[end] - 1,cumulsizes(block_sizes,2)[end] - 1)
end

function Base.show(io::IO, block_sizes::AbstractBlockSizes{N}) where {N}
    if N == 0
        print(io, "[]")
    else
        print(io, diff(cumulsizes(block_sizes,1)))
        for i in 2:N
            print(io, " × ", diff(cumulsizes(block_sizes,i)))
        end
    end
end

@inline function searchlinear(vec::AbstractVector, a)
    l = length(vec)
    @inbounds for i in 1:l
        vec[i] > a && return i - 1
    end
    return l
end

_find_block(bs::AbstractVector, i::Integer) = length(bs) > 10 ? last(searchsorted(bs, i)) : searchlinear(bs, i)

@inline function _find_block(block_sizes::AbstractBlockSizes, dim::Integer, i::Integer)
    bs = cumulsizes(block_sizes, dim)
    block = _find_block(bs, i)
    @inbounds cum_size = cumulsizes(block_sizes, dim, block) - 1
    return block, i - cum_size
end

@generated function nblocks(block_sizes::AbstractBlockSizes{N}) where {N}
    ex = Expr(:tuple, [:(nblocks(block_sizes, $i)) for i in 1:N]...)
    return quote
        @inbounds return $ex
    end
end

@inline @propagate_inbounds nblocks(block_sizes::AbstractBlockSizes, i::Integer) =
    length(cumulsizes(block_sizes,i)) - 1

function nblocks(block_sizes::AbstractBlockSizes, i::Vararg{Integer, N}) where {N}
    b = nblocks(block_sizes)
    return ntuple(k-> b[i[k]], Val(N))
end

# ntuple is yet again slower
@generated function Base.copy(block_sizes::BlockSizes{N}) where {N}
    exp = Expr(:tuple, [:(copy(cumulsizes(block_sizes, $k))) for k in 1:N]...)
    return quote
        BlockSizes($exp)
    end
end

# Computes the global range of an Array that corresponds to a given block_index
@generated function globalrange(block_sizes::AbstractBlockSizes{N}, block_index::NTuple{N, Integer}) where {N}
    indices_ex = Expr(:tuple, [:(cumulsizes(block_sizes, $i, block_index[$i]):cumulsizes(block_sizes, $i, block_index[$i] + 1) - 1) for i = 1:N]...)
    return quote
        $(Expr(:meta, :inline))
        @inbounds inds = $indices_ex
        return inds
    end
end

# I hate having these function definitions but the generated function above sometimes(!) generates bad code and starts to allocate
@inline function globalrange(block_sizes::AbstractBlockSizes{1}, block_index::NTuple{1, Integer})
    @inbounds v = (cumulsizes(block_sizes, 1, block_index[1]):cumulsizes(block_sizes, 1, block_index[1] + 1) - 1,)
    return v
end

@inline function globalrange(block_sizes::AbstractBlockSizes{2}, block_index::NTuple{2, Integer})
    @inbounds v = (cumulsizes(block_sizes, 1, block_index[1]):cumulsizes(block_sizes, 1, block_index[1] + 1) - 1,
                   cumulsizes(block_sizes, 2, block_index[2]):cumulsizes(block_sizes, 2, block_index[2] + 1) - 1)
    return v
end

@inline function globalrange(block_sizes::AbstractBlockSizes{3}, block_index::NTuple{3, Integer})
    @inbounds v = (cumulsizes(block_sizes, 1, block_index[1]):cumulsizes(block_sizes, 1, block_index[1] + 1) - 1,
                   cumulsizes(block_sizes, 2, block_index[2]):cumulsizes(block_sizes, 2, block_index[2] + 1) - 1,
                   cumulsizes(block_sizes, 3, block_index[3]):cumulsizes(block_sizes, 3, block_index[3] + 1) - 1)
    return v
end

"""
    blocksizes(A)

returns a subtype of `AbstractBlockSizes` that contains information about the
block sizes of `A`. Any subtype of AbstractBlockArrays must override this.
"""
blocksizes(A::AbstractBlockArray) = error("blocksizes for $(typeof(A)) is not implemented")

@inline nblocks(block_array::AbstractArray) = nblocks(blocksizes(block_array))

"""
    blocksize(A, inds)

Returns a tuple containing the size of the block at block index `inds`.

```jldoctest; setup = quote using BlockArrays end
julia> A = BlockArray(rand(5, 4, 6), [1, 4], [1, 2, 1], [1, 2, 2, 1]);

julia> blocksize(A, (1, 3, 2))
(1, 1, 2)

julia> blocksize(A, (2, 1, 3))
(4, 1, 2)
```
"""
@inline blocksize(block_array::AbstractArray, i::Integer...) =
    blocksize(blocksizes(block_array), i...)

@inline blocksize(block_array::AbstractArray{T,N}, i::NTuple{N, Integer}) where {T, N} =
    blocksize(blocksizes(block_array), i)

@inline Base.size(arr::AbstractBlockArray) = size(blocksizes(arr))

cumulsizes(A::AbstractArray) = cumulsizes(blocksizes(A))
@inline cumulsizes(A::AbstractArray, i) = cumulsizes(blocksizes(A), i)
@inline cumulsizes(A::AbstractArray, i, j) = cumulsizes(blocksizes(A), i, j)


text:
function initialize!(integrator, cache::DImplicitEulerConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.k[1] = integrator.du
end

function initialize!(integrator, cache::DImplicitEulerCache)
  integrator.kshortsize = 2
  @unpack k₁,k₂ = cache
  resize!(integrator.k, integrator.kshortsize)
  integrator.k .= [k₁,k₂]
  integrator.k[1] .= integrator.du
  nothing
end

@muladd function perform_step!(integrator, cache::DImplicitEulerConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack nlsolver = cache

  nlsolver.z = zero(u)
  nlsolver.tmp = zero(u)
  nlsolver.γ = 1
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  u = uprev + z

  if integrator.opts.adaptive && integrator.success_iter > 0
    # local truncation error (LTE) bound by dt^2/2*max|y''(t)|
    # use 2nd divided differences (DD) a la SPICE and Shampine

    # TODO: check numerical stability
    uprev2 = integrator.uprev2
    tprev = integrator.tprev

    dt1 = dt*(t+dt-tprev)
    dt2 = (t-tprev)*(t+dt-tprev)
    c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
    r = c*dt^2 # by mean value theorem 2nd DD equals y''(s)/2 for some s

    tmp = r*integrator.opts.internalnorm.((u - uprev)/dt1 - (uprev - uprev2)/dt2,t)
    atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  else
    integrator.EEst = 1
  end

  integrator.u = u
  integrator.du = (u-uprev)/dt

  if integrator.opts.calck
    integrator.k[2] = integrator.k[1]
    integrator.k[1] = integrator.du
  end
end

@muladd function perform_step!(integrator, cache::DImplicitEulerCache, repeat_step=false)
  @unpack t,dt,uprev,du,u,f,p = integrator
  @unpack atmp,nlsolver = cache
  @unpack tmp = nlsolver

  @. nlsolver.z = false
  @. nlsolver.tmp = false
  nlsolver.γ = 1
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  @.. u = uprev + z
  @.. du = z * inv(dt)

  if integrator.opts.adaptive && integrator.success_iter > 0
    # local truncation error (LTE) bound by dt^2/2*max|y''(t)|
    # use 2nd divided differences (DD) a la SPICE and Shampine

    # TODO: check numerical stability
    uprev2 = integrator.uprev2
    tprev = integrator.tprev

    dt1 = dt*(t+dt-tprev)
    dt2 = (t-tprev)*(t+dt-tprev)
    c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
    r = c*dt^2 # by mean value theorem 2nd DD equals y''(s)/2 for some s

    @.. tmp = r*integrator.opts.internalnorm((u - uprev)/dt1 - (uprev - uprev2)/dt2,t)
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  else
    integrator.EEst = 1
  end

  if integrator.opts.calck
    integrator.k[2] .= integrator.k[1]
    integrator.k[1] .= du
  end
end

function initialize!(integrator, cache::DABDF2ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.k[1] = integrator.du
end

@muladd function perform_step!(integrator, cache::DABDF2ConstantCache, repeat_step=false)
  @unpack t,f,p = integrator
  @unpack dtₙ₋₁,nlsolver = cache
  dtₙ, uₙ, uₙ₋₁, uₙ₋₂ = integrator.dt, integrator.u, integrator.uprev, integrator.uprev2

  if integrator.iter == 1 && !integrator.u_modified
    cache.dtₙ₋₁ = dtₙ
    perform_step!(integrator, cache.eulercache, repeat_step)
    integrator.fsalfirst = @.. (integrator.u - integrator.uprev) / dtₙ
    cache.fsalfirstprev = integrator.fsalfirst
    return
  end

  # precalculations
  ρ = dtₙ/dtₙ₋₁
  c1 = ρ^2/(1+2ρ)

  nlsolver.γ = (1+ρ)/(1+2ρ)
  nlsolver.α = 1//1

  nlsolver.z = zero(uₙ)

  nlsolver.tmp = -c1 * uₙ₋₁ + c1 * uₙ₋₂
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  uₙ = uₙ₋₁ + z
  integrator.fsallast = @.. z/dtₙ

  if integrator.opts.adaptive
    tmp = integrator.fsallast - (1+dtₙ/dtₙ₋₁)*integrator.fsalfirst + (dtₙ/dtₙ₋₁)*cache.fsalfirstprev
    est = (dtₙ₋₁+dtₙ)/6 * tmp
    atmp = calculate_residuals(est, uₙ₋₁, uₙ, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  ################################### Finalize

  if integrator.EEst < one(integrator.EEst)
    cache.fsalfirstprev = integrator.fsalfirst
    cache.dtₙ₋₁ = dtₙ
  end

  integrator.u = uₙ
  integrator.du = du = (nlsolver.α * z + nlsolver.tmp) * inv(nlsolver.γ * dtₙ)

  if integrator.opts.calck
    integrator.k[2] = integrator.k[1]
    integrator.k[1] = integrator.du
  end
  return
end

function initialize!(integrator, cache::DABDF2Cache)
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = du_alias_or_new(cache.nlsolver, integrator.fsalfirst)

  integrator.kshortsize = 2
  @unpack k₁,k₂ = cache.eulercache
  resize!(integrator.k, integrator.kshortsize)
  integrator.k .= [k₁,k₂]
  integrator.k[1] .= integrator.du
  nothing
end

@muladd function perform_step!(integrator, cache::DABDF2Cache, repeat_step=false)
  @unpack t,dt,du,f,p = integrator
  @unpack atmp,dtₙ₋₁,nlsolver = cache
  @unpack z,tmp = nlsolver
  uₙ,uₙ₋₁,uₙ₋₂,dtₙ = integrator.u,integrator.uprev,integrator.uprev2,integrator.dt

  if integrator.iter == 1 && !integrator.u_modified
    cache.dtₙ₋₁ = dtₙ
    perform_step!(integrator, cache.eulercache, repeat_step)
    @.. integrator.fsalfirst = (uₙ - uₙ₋₁) / dt
    cache.fsalfirstprev .= integrator.fsalfirst
    return
  end

  # precalculations
  ρ = dtₙ/dtₙ₋₁
  c1 = ρ^2/(1+2ρ)

  nlsolver.γ = (1+ρ)/(1+2ρ)
  nlsolver.α = 1//1
  @.. nlsolver.tmp = -c1 * uₙ₋₁ + c1 * uₙ₋₂
  nlsolver.z .= zero(eltype(z))
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  @.. uₙ = uₙ₋₁ + z
  @.. du = (nlsolver.α * z + nlsolver.tmp) * inv(nlsolver.γ * dt)

  @.. integrator.fsallast = du
  integrator.destats.nf += 1
  if integrator.opts.adaptive
    btilde0 = (dtₙ₋₁+dtₙ)*1//6
    btilde1 = 1+ρ
    btilde2 = ρ
    @.. tmp = btilde0*(integrator.fsallast - btilde1*integrator.fsalfirst + btilde2*cache.fsalfirstprev)
    calculate_residuals!(atmp, tmp, uₙ₋₁, uₙ, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  ################################### Finalize

  if integrator.EEst < one(integrator.EEst)
    @.. cache.fsalfirstprev = integrator.fsalfirst
    cache.dtₙ₋₁ = dtₙ
  end
  return
end

function initialize!(integrator, cache::DFBDFConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.du, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function perform_step!(integrator, cache::DFBDFConstantCache{max_order}, repeat_step=false) where max_order
  @unpack ts,u_history,order,u_corrector,bdf_coeffs,r,nlsolver,weights,ts_tmp, iters_from_event, nconsteps = cache
  @unpack t,dt,u,f,p,uprev = integrator

  k = order
  reinitFBDF!(integrator, cache)
  cache.u₀ = zero(u)
  if iters_from_event >= 1
    cache.u₀ = calc_Lagrange_interp(k,weights,t+dt,ts,u_history,cache.u₀)
  else
    cache.u₀ = u
  end
  markfirststage!(nlsolver)
  equi_ts = zeros(k-1)
  for i in 1:k-1
    equi_ts[i] = t - dt*i
  end

  fill!(u_corrector,zero(eltype(u)))
  if u isa Number
    for i in 1:k-1
      u_corrector[i] = calc_Lagrange_interp(k,weights,equi_ts[i],ts,u_history,u_corrector[i])
    end
    tmp = uprev * bdf_coeffs[k,2]
    for i in 1:k-1
      tmp += u_corrector[i] * bdf_coeffs[k,i+2]
    end
  else
    for i in 1:k-1
      @.. @views u_corrector[:,i] = $calc_Lagrange_interp(k,weights,equi_ts[i],ts,u_history,u_corrector[:,i])
    end
    tmp = uprev * bdf_coeffs[k,2]
    vc = _vec(tmp)
    for i in 1:k-1
      @.. @views vc += u_corrector[:,i] * bdf_coeffs[k,i+2]
    end
  end

  nlsolver.tmp = tmp + cache.u₀
  nlsolver.z = zero(nlsolver.z)
  nlsolver.γ = bdf_coeffs[k,1]
  nlsolver.α = 1//1
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) &&  return
  u = z + cache.u₀

  for j in 2:k
    r[j] = (1-j)
    for i in 2:k+1
      r[j] *= ((t+dt-j*dt)-ts[i])/(i*dt)
    end
  end

  terkp1 = z
  for j in 1:k+1
    terkp1 *= j*dt/(t+dt-ts[j])
  end

  lte = -1/(1+k)
  for j in 2:k
    lte -= (bdf_coeffs[k,j]//bdf_coeffs[k,1])*r[j]
  end
  lte *= terkp1

  if integrator.opts.adaptive
    for i in 1:k+1
      ts_tmp[i+1] = ts[i]
    end
    ts_tmp[1] = t+dt
    atmp = calculate_residuals(_vec(lte), _vec(uprev), _vec(u), integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)

    terk = estimate_terk(integrator, cache, k+1, Val(max_order), u)
    atmp = calculate_residuals(_vec(terk), _vec(uprev), _vec(u), integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
    cache.terk = integrator.opts.internalnorm(atmp,t)

    if k > 1
      terkm1 = estimate_terk(integrator, cache, k, Val(max_order), u)
      atmp = calculate_residuals(_vec(terkm1), _vec(uprev), _vec(u), integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      cache.terkm1 = integrator.opts.internalnorm(atmp,t)
    end
    if k > 2
      terkm2 = estimate_terk(integrator, cache, k-1, Val(max_order), u)
      atmp = calculate_residuals(_vec(terkm2), _vec(uprev), _vec(u), integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      cache.terkm2 = integrator.opts.internalnorm(atmp,t)
    end
    if nconsteps > k+1 && k < max_order
      atmp = calculate_residuals(_vec(terkp1), _vec(uprev), _vec(u), integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      cache.terkp1 = integrator.opts.internalnorm(atmp,t)
    else
      cache.terkp1 = zero(cache.terk)
    end
  end
  integrator.u = u
  integrator.fsallast = integrator.du = (nlsolver.α * z + nlsolver.tmp) * inv(nlsolver.γ * dt)
end

function initialize!(integrator, cache::DFBDFCache)
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = du_alias_or_new(cache.nlsolver, integrator.fsalfirst)
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  #integrator.f(integrator.fsalfirst, integrator.du, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point
  #integrator.destats.nf += 1
end

function perform_step!(integrator, cache::DFBDFCache{max_order}, repeat_step=false) where max_order
  @unpack ts,u_history,order,u_corrector,bdf_coeffs,r,nlsolver,weights,terk_tmp,terkp1_tmp,atmp,tmp,equi_ts,u₀,ts_tmp = cache
  @unpack t,dt,u,f,p,uprev = integrator

  reinitFBDF!(integrator, cache)
  k = order
  @.. u₀ = zero(u)
  if cache.iters_from_event >= 1
    calc_Lagrange_interp!(k,weights,t+dt,ts,u_history,u₀)
  else
    @.. u₀ = u
  end
  markfirststage!(nlsolver)

  for i in 1:k-1
    equi_ts[i] = t - dt*i
  end

  fill!(u_corrector,zero(eltype(u)))
  for i in 1:k-1
    @views calc_Lagrange_interp!(k,weights,equi_ts[i],ts,u_history,u_corrector[:,i])
  end

  @.. tmp = uprev * bdf_coeffs[k,2]
  vc = _vec(tmp)
  for i in 1:k-1
    @.. @views vc += u_corrector[:,i] * bdf_coeffs[k,i+2]
  end

  @.. nlsolver.tmp = tmp + u₀
  @.. nlsolver.z = zero(eltype(nlsolver.z))
  nlsolver.γ = bdf_coeffs[k,1]
  nlsolver.α = 1//1
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  @.. u = z + u₀

  for j in 2:k
    r[j] = (1-j)
    for i in 2:k+1
      r[j] *= ((t+dt-j*dt)-ts[i])/(i*dt)
    end
  end

  @.. terkp1_tmp = z
  for j in 1:k+1
    @.. terkp1_tmp *= j*dt/(t+dt-ts[j])
  end

  lte = -1/(1+k)
  for j in 2:k
    lte -= (bdf_coeffs[k,j]//bdf_coeffs[k,1])*r[j]
  end
  @.. terk_tmp = lte * terkp1_tmp
  if integrator.opts.adaptive
    @unpack abstol, reltol, internalnorm = integrator.opts
    for i in 1:k+1
      ts_tmp[i+1] = ts[i]
    end
    ts_tmp[1] = t+dt
    calculate_residuals!(atmp, _vec(terk_tmp), _vec(uprev), _vec(u), abstol, reltol, internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
    estimate_terk!(integrator, cache, k+1, Val(max_order))
    calculate_residuals!(atmp, _vec(terk_tmp), _vec(uprev), _vec(u), abstol, reltol, internalnorm, t)
    cache.terk = integrator.opts.internalnorm(atmp,t)

    if k > 1
      estimate_terk!(integrator, cache, k, Val(max_order))
      calculate_residuals!(atmp, _vec(terk_tmp), _vec(uprev), _vec(u), integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      cache.terkm1 = integrator.opts.internalnorm(atmp,t)
    end
    if k > 2
      estimate_terk!(integrator, cache, k-1, Val(max_order))
      calculate_residuals!(atmp, _vec(terk_tmp), _vec(uprev), _vec(u), integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      cache.terkm2 = integrator.opts.internalnorm(atmp,t)
    end
    if cache.nconsteps > k+1 && k < max_order
      calculate_residuals!(atmp, _vec(terkp1_tmp), _vec(uprev), _vec(u), integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      cache.terkp1 = integrator.opts.internalnorm(atmp,t)
    else
      cache.terkp1 = zero(cache.terkp1)
    end
  end
  @.. integrator.fsallast = integrator.du = (nlsolver.α * z + nlsolver.tmp) * inv(nlsolver.γ * dt) #TODO Lorenz plot seems not smooth
end


text:
a = [1 2  3  4
     5 6  7  8
     9 10 11 12]
da = DimensionalArray(a, (Y(10:30), Time((1:4)u"s")))

@testset "selectors with dim wrappers" begin
    @test da[Y<|At([10, 30]), Time<|At([1u"s", 4u"s"])] == [1 4; 9 12]
    @test_throws ArgumentError da[Y<|At([9, 30]), Time<|At([1u"s", 4u"s"])]
    @test view(da, Y<|At(20), Time<|At((3:4)u"s")) == [7, 8]
    @test view(da, Y<|Near(17), Time<|Near([1.3u"s", 3.3u"s"])) == [5, 7]
    @test view(da, Y<|Between(9, 31), Time<|At((3:4)u"s")) == [3 4; 7 8; 11 12]
end

@testset "selectors without dim wrappers" begin
    @test da[At(20:10:30), At(1u"s")] == [5, 9]
    @test view(da, Between(9, 31), Near((3:4)u"s")) == [3 4; 7 8; 11 12]
    @test view(da, Near(22), At(3u"s", 4u"s")) == [7, 8]
    @test view(da, At(20), At((2:3)u"s")) == [6, 7]
    @test view(da, Near<|13, Near<|[1.3u"s", 3.3u"s"]) == [1, 3]
    # Near works with a tuple input
    @test view(da, Near<|(13,), Near<|(1.3u"s", 3.3u"s")) == [1 3]
    @test view(da, Between(11, 20), At((2:3)u"s")) == [6 7]
    # Between also accepts a tuple input
    @test view(da, Between((11, 20)), Between((2u"s", 3u"s"))) == [6 7]
end

@testset "setindex! with selectors" begin
    da[Near(11), At(3u"s")] = 100
    @test a[1, 3] == 100
    da[Time<|Near(2.2u"s"), Y<|Between(10, 30)] = [200, 201, 202]
    @test a[1:3, 2] == [200, 201, 202]
end

a = [1 2  3  4
     5 6  7  8
     9 10 11 12]

@testset "more Unitful dims" begin
    dimz = Time<|1.0u"s":1.0u"s":3.0u"s", Y<|(1u"km", 4u"km")
    da = DimensionalArray(a, dimz)
    @test da[Y<|Between(2u"km", 3.9u"km"), Time<|At<|3.0u"s"] == [10, 11]
end

@testset "ad-hoc categorical indices" begin
    dimz = Time<|[:one, :two, :three], Y<|[:a, :b, :c, :d]
    da = DimensionalArray(a, dimz)
    @test da[Time<|At(:one, :two), Y<|At(:b)] == [2, 6]
    @test da[At([:one, :three]), At([:b, :c, :d])] == [2 3 4; 10 11 12]
end


text:
# fft_convolve.jl

export fft_conv!, fft_conv_adj!
export fft_conv, fft_conv_adj

"""
    imfilterz!(plan)
FFT-based convolution of `plan.img_compl`
and kernel `plan.ker_compl` (not centered),
storing result in `plan.workmat`.
"""
function imfilterz!(plan::PlanPSF)
    mul!(plan.img_compl, plan.fft_plan, plan.img_compl)
    mul!(plan.ker_compl, plan.fft_plan, plan.ker_compl)
    broadcast!(*, plan.img_compl, plan.img_compl, plan.ker_compl)
    mul!(plan.img_compl, plan.ifft_plan, plan.img_compl)
    fftshift2!(plan.ker_compl, plan.img_compl)
    plan.workmat .= real.(plan.ker_compl)
    return plan.workmat
end

"""
    fft_conv!(output, img, ker, plan)
Convolve 2D image `img` with 2D (symmetric!) kernel `ker` using FFT,
storing the result in `output`.
"""
function fft_conv!(
    output::AbstractMatrix{<:RealU},
    img::AbstractMatrix{<:RealU},
    ker::AbstractMatrix{<:RealU},
    plan::PlanPSF,
)
    @boundscheck size(output) == size(img) || throw("size output")
    @boundscheck size(img) == (plan.nx, plan.nz) || throw("size img")
    @boundscheck size(ker) == (plan.px, plan.pz) ||
        throw("size ker $(size(ker)) $(plan.px) $(plan.pz)")

    # filter image with a kernel, using replicate padding and fft convolution
    padrepl!(plan.img_compl, img, plan.padsize)

    pad2sizezero!(plan.ker_compl, ker, size(plan.ker_compl)) # zero pad kernel

    imfilterz!(plan)

    (M, N) = size(img)
    copyto!(output, (@view plan.workmat[plan.padsize[1] .+ (1:M),
                                        plan.padsize[3] .+ (1:N)]))
    return output
end

"""
    fft_conv(img, ker; T)
Convolve 2D image `img` with 2D (symmetric!) kernel `ker` using FFT.
"""
function fft_conv(
    img::AbstractMatrix{I},
    ker::AbstractMatrix{K};
    T::DataType = promote_type(I, K, Float32),
) where {I <: Number, K <: Number}

    ker ≈ reverse(ker, dims=:) || throw("asymmetric kernel")
    nx, nz = size(img)
    px, pz = size(ker)
    plan = plan_psf( ; nx, nz, px, pz, T, nthread = 1)[1]
    output = similar(Matrix{T}, size(img))
    fft_conv!(output, img, ker, plan)
    return output
end

"""
    fft_conv_adj!(output, img, ker, plan)
Adjoint of `fft_conv!`.
"""
function fft_conv_adj!(
    output::AbstractMatrix{<:RealU},
    img::AbstractMatrix{<:RealU},
    ker::AbstractMatrix{<:RealU},
    plan::PlanPSF{T},
) where {T}

    @boundscheck size(output) == size(img) || throw("size output")
    @boundscheck size(img) == (plan.nx, plan.nz) || throw("size img")
    @boundscheck size(ker) == (plan.px, plan.pz) ||
        throw("size ker $(size(ker)) $(plan.px) $(plan.pz)")

    padzero!(plan.img_compl, img, plan.padsize) # pad the image with zeros
    pad2sizezero!(plan.ker_compl, ker, size(plan.ker_compl)) # pad the kernel with zeros

    imfilterz!(plan)
    (M, N) = size(img)
    # adjoint of replicate padding
    plan.workvecz .= zero(T)
    for i = 1:plan.padsize[1]
        plus2di!(plan.workvecz, plan.workmat, i)
    end
    plus1di!(plan.workmat, plan.workvecz, 1+plan.padsize[1])

    plan.workvecz .= zero(T)
    for i = (plan.padsize[1]+M+1):size(plan.workmat, 1)
        plus2di!(plan.workvecz, plan.workmat, i)
    end
    plus1di!(plan.workmat, plan.workvecz, M+plan.padsize[1])

    plan.workvecx .= zero(T)
    for j = 1:plan.padsize[3]
        plus2dj!(plan.workvecx, plan.workmat, j)
    end
    plus1dj!(plan.workmat, plan.workvecx, 1+plan.padsize[3])

    plan.workvecx .= zero(T)
    for j = (plan.padsize[3]+N+1):size(plan.workmat, 2)
        plus2dj!(plan.workvecx, plan.workmat, j)
    end
    plus1dj!(plan.workmat, plan.workvecx, N+plan.padsize[3])

    copyto!(output,
        (@view plan.workmat[(plan.padsize[1]+1):(plan.padsize[1]+M),
                            (plan.padsize[3]+1):(plan.padsize[3]+N)]),
    )

    return output
end

"""
    fft_conv_adj(img, ker; T)
Adjoint of `fft_conv`.
"""
function fft_conv_adj(
    img::AbstractMatrix{I},
    ker::AbstractMatrix{K};
    T::DataType = promote_type(I, K, Float32),
) where {I <: Number, K <: Number}

    ker ≈ reverse(ker, dims=:) || throw("asymmetric kernel")
    nx, nz = size(img)
    px, pz = size(ker)
    plan = plan_psf( ; nx, nz, px, pz, T, nthread = 1)[1]
    output = similar(Matrix{T}, size(img))
    fft_conv_adj!(output, img, ker, plan)
    return output
end

"""
    fft_conv!(output, image3, ker3, plans)
Mutating version of convolving a 3D `image3`
with a set of 2D symmetric kernels
stored in 3D array `ker3`
using `foreach`.
"""
function fft_conv!(
    output::AbstractArray{<:RealU,3},
    image3::AbstractArray{<:RealU,3},
    ker3::AbstractArray{<:RealU,3},
    plans::Vector{<:PlanPSF},
)

    size(output) == size(image3) || throw(DimensionMismatch())

    fun = y -> fft_conv!(
            (@view output[:, y, :]),
            (@view image3[:, y, :]),
            (@view ker3[:, :, y]),
            plans[Threads.threadid()],
        )

    ntasks = length(plans)
    Threads.foreach(fun, foreach_setup(1:size(image3, 2)); ntasks)

    return output
end

"""
    fft_conv_adj!(output, image3, ker3, plans)
Adjoint of `fft_conv`.
"""
function fft_conv_adj!(
    output::AbstractArray{<:RealU,3},
    image3::AbstractArray{<:RealU,3},
    ker3::AbstractArray{<:RealU,3},
    plans::Vector{<:PlanPSF},
)

    size(output) == size(image3) || throw(DimensionMismatch())

    fun = y -> fft_conv_adj!(
            (@view output[:, y, :]),
            (@view image3[:, y, :]),
            (@view ker3[:, :, y]),
            plans[Threads.threadid()],
        )

    ntasks = length(plans)
    Threads.foreach(fun, foreach_setup(1:size(image3, 2)); ntasks)

    return output
end

"""
    fft_conv_adj2!(output, image2, ker3, plans)
Mutating version of adjoint of convolving a 2D `image2`
with each 2D kernel in the 3D array `ker3`.
"""
function fft_conv_adj2!(
    output::AbstractArray{<:RealU,3},
    image2::AbstractMatrix{<:RealU},
    ker3::AbstractArray{<:RealU,3},
    plans::Vector{<:PlanPSF},
)

    size(output, 1) == size(image2, 1) || throw("size 1")
    size(output, 3) == size(image2, 2) || throw("size 2")

    fun = y -> fft_conv_adj!(
            (@view output[:, y, :]),
            image2,
            (@view ker3[:, :, y]),
            plans[Threads.threadid()],
        )

    ntasks = length(plans)
    Threads.foreach(fun, foreach_setup(1:size(output, 2)); ntasks)

    return output
end


text:
import ..options
import ..@log

using .Schedulers

"""
Core algorithm of the library: the equality saturation step.
"""
function eqsat_step!(g::EGraph, theory::Vector{<:AbstractRule}, curr_iter,
        scheduler::AbstractScheduler, match_hist::MatchesBuf,
        params::SaturationParams, report)

    instcache = Dict{AbstractRule, Dict{Sub, EClassId}}()

    setiter!(scheduler, curr_iter)

    matches = @timeit report.to "Search" eqsat_search!(g, theory, scheduler, report; threaded=params.threaded)

    # matches = setdiff!(matches, match_hist)

    @timeit report.to "Apply" eqsat_apply!(g, matches, report, params)

    # union!(match_hist, matches)

    if report.reason === nothing && cansaturate(scheduler) && isempty(g.dirty)
        report.reason = Saturated()
    end
    @timeit report.to "Rebuild" rebuild!(g)

    return report, g
end

"""
Given an [`EGraph`](@ref) and a collection of rewrite rules,
execute the equality saturation algorithm.
"""
function saturate!(g::EGraph, theory::Vector{<:AbstractRule}, params=SaturationParams())
    curr_iter = 0

    sched = params.scheduler(g, theory, params.schedulerparams...)
    match_hist = MatchesBuf()
    report = Report(g)

    if !params.timer
        disable_timer!(report.to)
    end

    while true
        curr_iter+=1

        options.printiter && @info("iteration ", curr_iter)

        report, egraph = eqsat_step!(g, theory, curr_iter, sched, match_hist, params, report)

        # report.reason == :matchlimit && break
        if !(report.reason isa Nothing)
            break
        end

        if curr_iter >= params.timeout
            report.reason = Timeout()
            break
        end

        if params.eclasslimit > 0 && g.numclasses > params.eclasslimit
            report.reason = EClassLimit(params.eclasslimit)
            break
        end

        if reached(g, params.goal)
            report.reason = GoalReached()
            break
        end
    end
    report.iterations = curr_iter
    @log report

    return report
end


text:
# Determine MATLAB library path and provide facilities to load libraries with
# this path

function get_paths()
    global matlab_homepath = get(ENV, "MATLAB_HOME", "")
    if matlab_homepath == ""
        if is_linux()
            matlab_homepath = dirname(dirname(realpath(chomp(readstring(`which matlab`)))))
        elseif is_apple()
            apps = readdir("/Applications")
            filter!(app -> ismatch(r"^MATLAB_R[0-9]+[ab]\.app$", app), apps)
            if ~isempty(apps)
                matlab_homepath = joinpath("/Applications", maximum(apps))
            end
        elseif is_windows()
            default_dir = Int == Int32 ? "C:\\Program Files (x86)\\MATLAB" : "C:\\Program Files\\MATLAB"
            if isdir(default_dir)
                dirs = readdir(default_dir)
                filter!(dir -> ismatch(r"^R[0-9]+[ab]$", dir), dirs)
                if ~isempty(dirs)
                    matlab_homepath = joinpath(default_dir, maximum(dirs))
                end
            end
        end
    end

    if matlab_homepath == ""
        error("The MATLAB path could not be found. Set the MATLAB_HOME environmental variable to specify the MATLAB path.")
    end

    # Get path to MATLAB libraries
    global matlab_library_path = nothing
    if is_linux()
        matlab_library_path = joinpath(matlab_homepath, "bin", (Int == Int32 ? "glnx86" : "glnxa64"))
    elseif is_apple()
        matlab_library_path = joinpath(matlab_homepath, "bin", (Int == Int32 ? "maci" : "maci64"))
    elseif is_windows()
        matlab_library_path = joinpath(matlab_homepath, "bin", (Int == Int32 ? "win32" : "win64"))
    end

    if matlab_library_path != nothing && !isdir(matlab_library_path)
        matlab_library_path = nothing
    end
end
get_paths()

function open_matlab_library(lib::String)
    lib_path = matlab_library_path == nothing ? lib : joinpath(matlab_library_path, lib)
    ptr = Libdl.dlopen(lib_path, Libdl.RTLD_GLOBAL | Libdl.RTLD_LAZY)
    if ptr == C_NULL
        error("Failed to load $(lib)")
    end
    ptr
end


text:
using CSV, Distributions, HypothesisTests

data1 = CSV.read("../data/machine1.csv", header=false)[:,1]
data2 = CSV.read("../data/machine2.csv", header=false)[:,1]
xBar1, s1, n1 = mean(data1), std(data1), length(data1)
xBar2, s2, n2 = mean(data2), std(data2), length(data2)
delta0 = 0

v = ( s1^2/n1 + s2^2/n2 )^2 / ( (s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1)  )
testStatistic = ( xBar1-xBar2 - delta0 )  / sqrt( s1^2/n1 + s2^2/n2 )
pVal = 2*ccdf(TDist(v), abs(testStatistic))

println("Manually calculated degrees of freedom, v: ", v)
println("Manually calculated test statistic: ", testStatistic)
println("Manually calculated p-value: ", pVal, "\n")
println(UnequalVarianceTTest(data1, data2, delta0))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
mutable struct ADkernel{Kraw<:Kernel, Kdual<:Kernel, V<:AbstractVector, D<:Dual,CFG<:GradientConfig} <: Kernel
    raw::Kraw
    dual::Kdual
    hyp::V
    cfg::CFG
end
@inline raw(k::ADkernel) = k.raw
@inline dual(k::ADkernel) = k.dual

function to_autodiff(k::Kernel, duals::Vector{D}) where {D<:Dual}
    kerneltype = typeof(k)
    @assert !kerneltype.abstract
    @assert !kerneltype.hasfreetypevars

    typeparams = kerneltype.parameters

    fnames = fieldnames(kerneltype)
    values = []
    newtypes = []
    for field in fnames
        ftype = fieldtype(kerneltype, field)
        if ftype<:Float64
            push!(newtypes, D)
            val = getfield(k, field)
            dualval = D(val, zero(D).partials)
            push!(values, dualval)
        elseif ftype<:Vector{Float64}
            push!(newtypes, Vector{D})
            vals = getfield(k, field)
            dualvals = [D(val, zero(D).partials) for val in vals]
            push!(values, dualvals) # whatever
        elseif ftype<:Kernel
            subkernel = getfield(k, field)
            ad_subkernel = to_autodiff(subkernel, duals)
            push!(newtypes, typeof(ad_subkernel))
            push!(values, ad_subkernel)
        else
            push!(newtypes, ftype)
            push!(values, getfield(k, field))
        end
    end
    kerneltype.name.wrapper(values...)
end
function autodiff(k::Kernel)
    # TODO: smaller chunk sizes
    cfg = GradientConfig(nothing, zeros(num_params(k)), Chunk{num_params(k)}())
    D = eltype(cfg)
    CFG = typeof(cfg)
    hyp = get_params(k)
    duals = cfg.duals
    kdual = to_autodiff(k, duals)
    ad = ADkernel{typeof(k),typeof(kdual),typeof(hyp),eltype(cfg),typeof(cfg)}(k, kdual, hyp, cfg)
    return ad
end

# delegate (I think there are macros that would do this automatically)
num_params(ad::ADkernel) = num_params(raw(ad))
get_params(ad::ADkernel) = get_params(raw(ad))
get_param_names(ad::ADkernel) = get_param_names(raw(ad))
function set_params!(ad::ADkernel, hyp)
    set_params!(raw(ad), hyp)
    ad.hyp = get_params(raw(ad))
end
function seed!(ad::ADkernel)
    seed!(ad.cfg.duals, ad.hyp, ad.cfg.seeds)
    set_params!(dual(ad), ad.cfg.duals)
end
@inline cov(ad::ADkernel, r::Real) = cov(raw(ad), r)
@inline cov(ad::ADkernel, x1::AbstractVector, x2::AbstractVector) = cov(raw(ad), x1, x2)
@inline cov_ij(ad::ADkernel, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData, i::Int, j::Int, dim::Int) = cov_ij(raw(ad), X1, X2, data, i, j, dim)

KernelData(ad::ADkernel, X1::AbstractMatrix, X2::AbstractMatrix) = KernelData(raw(ad), X1, X2)
kernel_data_key(ad::ADkernel, X1::AbstractMatrix, X2::AbstractMatrix) = kernel_data_key(raw(ad), X1, X2)

@inline function dKij_dθ!(dK::AbstractVector, ad::ADkernel, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData,
                  i::Int, j::Int, dim::Int, npars::Int)
    seed!(ad)
    k_eval = cov_ij(dual(ad), X1, X2, data, i, j, dim)
    copyto!(dK, partials(k_eval))
    return dK
end
@inline function dKij_dθp(ad::ADkernel, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData, i::Int, j::Int, p::Int, dim::Int)
    seed!(ad)
    k_eval = cov_ij(dual(ad), X1, X2, data, i, j, dim)
    return partials(k_eval)[p]
end

# priors
get_priors(ad::ADkernel) = get_priors(raw(ad))
set_priors!(ad::ADkernel, priors) = set_priors!(raw(ad), priors)


text:
using Test

using CiphertextStealing
using Random: randstring

@testset "Smoke" begin
    for text_len in 1:100
        plaintext = Vector{UInt8}(randstring(text_len))
        iv = rand(UInt8, 16)
        key = rand(UInt8, 16)

        iv_backup = copy(iv)
        ciphertext = aes_128_cbc_encode_cts!(plaintext, key, iv)
        @test length(ciphertext) == text_len
        @test (text_len > 16) == (iv == iv_backup)
        iv_backup = copy(iv)
        deciphered = aes_128_cbc_decode_cts!(ciphertext, key, iv)
        @test (text_len > 16) == (iv == iv_backup)
        @test deciphered == plaintext
    end
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
using PowerModels

abstract type AbstractSparsePolyModel end

mutable struct SparsePolyModel <: AbstractSparsePolyModel
    n
    m
    numeq
    nbus
    ng
    nb
    supp
    coe
    dg
end

function normalize(coe)
    mc = maximum(abs.(coe))
    return coe./mc
end

function move_zero!(supp,coe)
    ind=[abs(coe[i])>=1e-8 for i=1:length(coe)]
    return supp[ind],coe[ind]
end

function fl_sum(vector)
    return mapreduce(x->x, +, vector, init = 0.0)
end

function bfind(A, l, a)
    if l==0
        return 0
    end
    low=1
    high=l
    while low<=high
        mid=Int(ceil(1/2*(low+high)))
        if ndims(A)==2
            temp=A[:, mid]
        else
            temp=A[mid]
        end
        if temp==a
           return mid
        elseif temp<a
           low=mid+1
        else
           high=mid-1
        end
    end
    return 0
end

# the real POP formulization (voltage-power)
function pop_opf_real(data::Dict{String, Any}; normal=true, AngleCons=false, LineLimit=false)
    silence()
    PowerModels.standardize_cost_terms!(data, order=2)
    ref = PowerModels.build_ref(data)[:it][pm_it_sym][:nw][0]
    nbus=length(ref[:bus])
    nb=length(ref[:branch])
    ng=length(ref[:gen])
    n=2*nbus+2*ng
    m=4*nbus+2*ng+length(keys(ref[:ref_buses]))
    if AngleCons==true
        m+=2*nb
    end
    if LineLimit==true||LineLimit=="relax"
        m+=2*nb
    end
    numeq=2*nbus+length(keys(ref[:ref_buses]))
    dg=2*ones(Int, m)
    startpoint=zeros(n)
    supp=Vector{Vector{Vector{UInt16}}}(undef, m+1)
    coe=Vector{Vector{Float64}}(undef, m+1)

    gens=collect(keys(ref[:gen]))
    sort!(gens)
    # objective function
    nc=2*ng+1
    coe[1]=Vector{Float64}(undef, nc)
    supp[1]=Vector{Vector{UInt16}}(undef, nc)
    coe[1][1]=sum(gen["cost"][3] for (i,gen) in ref[:gen])
    supp[1][1]=[]
    for i=1:ng
        gen=ref[:gen][gens[i]]
        coe[1][2*(i-1)+2:2*(i-1)+3]=[gen["cost"][2];gen["cost"][1]]
        supp[1][2*(i-1)+2:2*(i-1)+3]=[[2*nbus+i], [2*nbus+i;2*nbus+i]]
    end
    supp[1],coe[1]=move_zero!(supp[1],coe[1])
    k=2

    bus=collect(keys(ref[:bus]))
    sort!(bus)
    # voltage magnitude constraints
    for i=1:nbus
        supp[k]=[[], [i;i], [i+nbus;i+nbus]]
        coe[k]=[-ref[:bus][bus[i]]["vmin"]^2;1;1]
        supp[k+1]=[[], [i;i], [i+nbus;i+nbus]]
        coe[k+1]=[ref[:bus][bus[i]]["vmax"]^2;-1;-1]
        startpoint[i]=sqrt(ref[:bus][bus[i]]["vmin"]*ref[:bus][bus[i]]["vmax"]/2)
        startpoint[i+nbus]=sqrt(ref[:bus][bus[i]]["vmin"]*ref[:bus][bus[i]]["vmax"]/2)
        k+=2
    end

    if AngleCons==true||LineLimit==true||LineLimit=="relax"
        for (i, branch) in ref[:branch]
            g, b = PowerModels.calc_branch_y(branch)
            tr, ti = PowerModels.calc_branch_t(branch)
            g_fr = branch["g_fr"]
            b_fr = branch["b_fr"]
            g_to = branch["g_to"]
            b_to = branch["b_to"]
            tm = branch["tap"]
            vr = bfind(bus,nbus,branch["f_bus"])
            vt = bfind(bus,nbus,branch["t_bus"])
            srt=sort([vr;vt])
            dsrt=sort([vr;vr;vt;vt])
            ab1=(g+g_fr)^2+(b+b_fr)^2
            cd1=(-g*tr+b*ti)^2+(b*tr+g*ti)^2
            acbd1=(g+g_fr)*(-g*tr+b*ti)-(b+b_fr)*(b*tr+g*ti)
            bcad1=-(b+b_fr)*(-g*tr+b*ti)-(g+g_fr)*(b*tr+g*ti)
            ab2=(g+g_to)^2*tm^4+(b+b_to)^2*tm^4
            cd2=(g*tr+b*ti)^2+(-b*tr+g*ti)^2
            acbd2=-(g+g_to)*tm^2*(g*tr+b*ti)+(b+b_to)*tm^2*(-b*tr+g*ti)
            bcad2=(b+b_to)*tm^2*(g*tr+b*ti)+(g+g_to)*tm^2*(-b*tr+g*ti)

            # angle differences
            if AngleCons==true
                coe[k]=[tan(branch["angmax"]);tan(branch["angmax"]);-1;1]
                supp[k]=[srt, srt.+nbus, [vt;vr+nbus], [vr;vt+nbus]]
                coe[k+1]=[1;-1;-tan(branch["angmin"]);-tan(branch["angmin"])]
                supp[k+1]=[[vt;vr+nbus], [vr;vt+nbus], srt, srt.+nbus]
                if normal==true
                    coe[k]=normalize(coe[k])
                    coe[k+1]=normalize(coe[k+1])
                end
                k+=2
            end

            # thermal limits
            if LineLimit==true
                coe[k]=[branch["rate_a"]^2*tm^4;-ab1;-2*ab1;-ab1;-cd1;-cd1;-cd1;-cd1;-2*acbd1;-2*acbd1;-2*acbd1;-2*acbd1;-2*bcad1;2*bcad1;-2*bcad1;2*bcad1]
                supp[k]=[[], [vr;vr;vr;vr], [vr;vr;vr+nbus;vr+nbus], [vr+nbus;vr+nbus;vr+nbus;vr+nbus], [vt;vt;vr+nbus;vr+nbus], [vr;vr;vt+nbus;vt+nbus],
                dsrt, dsrt.+nbus, sort([vr;vr;vr;vt]), [vr;vr;sort([vr+nbus;vt+nbus])], [srt;vr+nbus;vr+nbus], sort([vr;vr;vr;vt]).+nbus, [sort([vr;vr;vt]);vr+nbus],
                [vr;vr;vr;vt+nbus], [vt;vr+nbus;vr+nbus;vr+nbus], [vr;sort([vr+nbus;vr+nbus;vt+nbus])]]
                supp[k],coe[k]=move_zero!(supp[k],coe[k])
                coe[k+1]=[branch["rate_a"]^2*tm^4;-ab2;-2*ab2;-ab2;-cd2;-cd2;-cd2;-cd2;-2*acbd2;-2*acbd2;-2*acbd2;-2*acbd2;-2*bcad2;2*bcad2;-2*bcad2;2*bcad2]
                supp[k+1]=[[], [vt;vt;vt;vt], [vt;vt;vt+nbus;vt+nbus], [vt+nbus;vt+nbus;vt+nbus;vt+nbus], [vr;vr;vt+nbus;vt+nbus], [vt;vt;vr+nbus;vr+nbus],
                dsrt, dsrt.+nbus, sort([vt;vt;vt;vr]), [vt;vt;sort([vt+nbus;vr+nbus])], [srt;vt+nbus;vt+nbus], sort([vt;vt;vt;vr]).+nbus, [sort([vt;vt;vr]);vt+nbus],
                [vt;vt;vt;vr+nbus], [vr;vt+nbus;vt+nbus;vt+nbus], [vt;sort([vt+nbus;vt+nbus;vr+nbus])]]
                supp[k+1],coe[k+1]=move_zero!(supp[k+1],coe[k+1])
                dg[k-1:k]=[4;4]
                if normal==true
                    coe[k]=normalize(coe[k])
                    coe[k+1]=normalize(coe[k+1])
                end
                k+=2
            elseif LineLimit=="relax"
                mvr=ref[:bus][bus[vr]]["vmin"]^2
                mvt=ref[:bus][bus[vt]]["vmin"]^2
                coe[k]=[branch["rate_a"]^2*tm^4/mvr;-ab1;-ab1;-cd1;-cd1;-2*acbd1;2*bcad1;-2*bcad1;-2*acbd1]
                supp[k]=[[], [vr;vr], [vr+nbus;vr+nbus], [vt;vt], [vt+nbus;vt+nbus], srt, [vr;vt+nbus], [vt;vr+nbus], srt.+nbus]
                supp[k],coe[k]=move_zero!(supp[k],coe[k])
                coe[k+1]=[branch["rate_a"]^2*tm^4/mvt;-ab2;-ab2;-cd2;-cd2;-2*acbd2;2*bcad2;-2*bcad2;-2*acbd2]
                supp[k+1]=[[], [vt;vt], [vt+nbus;vt+nbus], [vr;vr], [vr+nbus;vr+nbus], srt, [vt;vr+nbus], [vr;vt+nbus], srt.+nbus]
                supp[k+1],coe[k+1]=move_zero!(supp[k+1],coe[k+1])
                if normal==true
                    coe[k]=normalize(coe[k])
                    coe[k+1]=normalize(coe[k+1])
                end
                k+=2
            end
        end
    end

    # power generation bound
    zero_pgen=UInt16[]
    for i=1:ng
        gen=ref[:gen][gens[i]]
        if gen["pmax"]>=1e-6
            coe[k]=[-gen["pmin"]*gen["pmax"];gen["pmin"]+gen["pmax"];-1]
            if normal==true
                coe[k]=normalize(coe[k])
            end
            supp[k]=[[], [i+2*nbus], [i+2*nbus;i+2*nbus]]
            supp[k],coe[k]=move_zero!(supp[k],coe[k])
            startpoint[i+2*nbus]=(gen["pmin"]+gen["pmax"])/2
            k+=1
        else
            push!(zero_pgen, i)
            startpoint[i+2*nbus]=0
            numeq+=1
        end
        coe[k]=[-gen["qmin"]*gen["qmax"];gen["qmin"]+gen["qmax"];-1]
        if normal==true
            coe[k]=normalize(coe[k])
        end
        supp[k]=[[], [i+2*nbus+ng], [i+2*nbus+ng;i+2*nbus+ng]]
        supp[k],coe[k]=move_zero!(supp[k],coe[k])
        startpoint[i+2*nbus+ng]=0
        k+=1
    end

    # active/reactive power
    for (r, i) in enumerate(bus)
        bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]]
        bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]]
        coe[k]=zeros(4*length(ref[:bus_arcs][i])+3)
        supp[k]=Vector{Vector{UInt16}}(undef, 4*length(ref[:bus_arcs][i])+3)
        supp[k][1:3]=[[], [r;r], [r+nbus;r+nbus]]
        coe[k][1]=fl_sum(load["pd"] for load in bus_loads)
        coe[k+1]=zeros(4*length(ref[:bus_arcs][i])+3)
        supp[k+1]=Vector{Vector{UInt16}}(undef, 4*length(ref[:bus_arcs][i])+3)
        supp[k+1][1:3]=[[], [r;r], [r+nbus;r+nbus]]
        coe[k+1][1]=fl_sum(load["qd"] for load in bus_loads)
        sgs=fl_sum(shunt["gs"] for shunt in bus_shunts)
        sbs=fl_sum(shunt["bs"] for shunt in bus_shunts)
        coe[k][2:3]=[sgs;sgs]
        coe[k+1][2:3]=[-sbs;-sbs]
        j=1
        for flow in ref[:bus_arcs][i]
            branch=ref[:branch][flow[1]]
            vr = bfind(bus,nbus,branch["f_bus"])
            vt = bfind(bus,nbus,branch["t_bus"])
            srt=sort([vr;vt])
            g, b = PowerModels.calc_branch_y(branch)
            tr, ti = PowerModels.calc_branch_t(branch)
            g_fr = branch["g_fr"]
            b_fr = branch["b_fr"]
            g_to = branch["g_to"]
            b_to = branch["b_to"]
            tm = branch["tap"]
            a1=(g+g_fr)/tm^2
            b1=-(b+b_fr)/tm^2
            c1=(-g*tr+b*ti)/tm^2
            d1=(b*tr+g*ti)/tm^2
            a2=g+g_to
            b2=-(b+b_to)
            c2=-(g*tr+b*ti)/tm^2
            d2=-(-b*tr+g*ti)/tm^2

            supp[k][j+3:j+6]=[srt, [vt;vr+nbus], [vr;vt+nbus], srt.+nbus]
            supp[k+1][j+3:j+6]=[srt, [vt;vr+nbus], [vr;vt+nbus], srt.+nbus]
            if vr==r
                coe[k][2:3].+=a1
                coe[k][j+3:j+6]=[c1;-d1;d1;c1]
                coe[k+1][2:3].+=b1
                coe[k+1][j+3:j+6]=[d1;c1;-c1;d1]
            else
                coe[k][2:3].+=a2
                coe[k][j+3:j+6]=[c2;d2;-d2;c2]
                coe[k+1][2:3].+=b2
                coe[k+1][j+3:j+6]=[d2;-c2;c2;d2]
            end
            j+=4
        end
        if !isempty(ref[:bus_gens][i])
            for gen_id in ref[:bus_gens][i]
                gen=bfind(gens, ng, gen_id)
                push!(supp[k], [gen+2*nbus])
                push!(coe[k], -1)
                push!(supp[k+1], [gen+2*nbus+ng])
                push!(coe[k+1], -1)
            end
        end
        if normal==true
            coe[k]=normalize(coe[k])
            coe[k+1]=normalize(coe[k+1])
        end
        supp[k],coe[k]=move_zero!(supp[k],coe[k])
        supp[k+1],coe[k+1]=move_zero!(supp[k+1],coe[k+1])
        k+=2
    end

    # reference voltage
    for key in keys(ref[:ref_buses])
        i=bfind(bus,nbus,key)
        supp[k]=[[i+nbus;i+nbus]]
        coe[k]=[1]
        k+=1
    end

    # zero power generation
    for i in zero_pgen
        supp[k]=[[i+2*nbus]]
        coe[k]=[1]
        dg[k-1]=1
        k+=1
    end

    return SparsePolyModel(n,m,numeq,nbus,ng,nb,supp,coe,dg),startpoint
end

# the complex POP formulization (voltage-power)
function pop_opf_com(data::Dict{String, Any}; normal=true, AngleCons=false, LineLimit=false)
    silence()
    PowerModels.standardize_cost_terms!(data, order=2)
    ref = PowerModels.build_ref(data)[:it][pm_it_sym][:nw][0]
    nbus=length(ref[:bus])
    nb=length(ref[:branch])
    ng=length(ref[:gen])
    n=nbus+ng
    m=4*nbus+2*ng+length(keys(ref[:ref_buses]))
    if AngleCons==true
        m+=2*nb
    end
    if LineLimit==true||LineLimit=="relax"
        m+=2*nb
    end
    numeq=2*nbus+length(keys(ref[:ref_buses]))
    dg=2*ones(Int, m)
    supp=Vector{Vector{Vector{Vector{UInt16}}}}(undef, m+1)
    coe=Vector{Vector{complex(Float64)}}(undef, m+1)

    gens=collect(keys(ref[:gen]))
    sort!(gens)
    # objective function
    nc=5*ng+1
    coe[1]=Vector{complex(Float64)}(undef, nc)
    supp[1]=Vector{Vector{Vector{UInt16}}}(undef, nc)
    coe[1][1]=sum(gen["cost"][3] for (i,gen) in ref[:gen])
    supp[1][1]=[[], []]
    for i=1:ng
        gen=ref[:gen][gens[i]]
        coe[1][5*(i-1)+2:5*(i-1)+3]=[0.5*gen["cost"][2], 0.5*gen["cost"][2]]
        coe[1][5*(i-1)+4:5*(i-1)+6]=[0.25*gen["cost"][1], 0.5*gen["cost"][1], 0.25*gen["cost"][1]]
        supp[1][5*(i-1)+2:5*(i-1)+6]=[[[nbus+i], []], [[], [nbus+i]], [[nbus+i;nbus+i], []], [[nbus+i], [nbus+i]], [[], [nbus+i;nbus+i]]]
    end
    supp[1],coe[1]=move_zero!(supp[1],coe[1])
    k=2

    bus=collect(keys(ref[:bus]))
    sort!(bus)
    # voltage magnitude constraints
    for i=1:nbus
        supp[k]=[[[], []], [[i], [i]]]
        coe[k]=[-ref[:bus][bus[i]]["vmin"]^2;1]
        supp[k+1]=[[[], []], [[i], [i]]]
        coe[k+1]=[ref[:bus][bus[i]]["vmax"]^2;-1]
        k+=2
    end

    if AngleCons==true||LineLimit==true||LineLimit=="relax"
        for (i, branch) in ref[:branch]
            g, b = PowerModels.calc_branch_y(branch)
            tr, ti = PowerModels.calc_branch_t(branch)
            g_fr = branch["g_fr"]
            b_fr = branch["b_fr"]
            g_to = branch["g_to"]
            b_to = branch["b_to"]
            tm = branch["tap"]
            vr = bfind(bus,nbus,branch["f_bus"])
            vt = bfind(bus,nbus,branch["t_bus"])
            srt=sort([vr, vt])
            a1=g+g_fr
            b1=-(b+b_fr)
            c1=-g*tr+b*ti
            d1=b*tr+g*ti
            a2=(g+g_to)*tm^2
            b2=-(b+b_to)*tm^2
            c2=-(g*tr+b*ti)
            d2=-(-b*tr+g*ti)

            # angle differences
            if AngleCons==true
                coe[k]=[tan(branch["angmax"])+im;tan(branch["angmax"])-im]
                supp[k]=[[[vr], [vt]], [[vt], [vr]]]
                coe[k+1]=[-tan(branch["angmin"])-im;-tan(branch["angmin"])+im]
                supp[k+1]=[[[vr], [vt]], [[vt], [vr]]]
                if normal==true
                    coe[k]=normalize(coe[k])
                    coe[k+1]=normalize(coe[k+1])
                end
                k+=2
            end

            # thermal limits
            if LineLimit==true
                coe[k]=[branch["rate_a"]^2*tm^4;-(a1^2+b1^2);-(a1*c1+b1*d1)+(b1*c1-a1*d1)*im;-(a1*c1+b1*d1)+(a1*d1-b1*c1)*im;-(c1^2+d1^2)]
                supp[k]=[[[], []], [[vr;vr], [vr;vr]], [[vr;vr], srt], [srt, [vr;vr]], [srt, srt]]
                supp[k],coe[k]=move_zero!(supp[k],coe[k])
                coe[k+1]=[branch["rate_a"]^2*tm^4;-(a2^2+b2^2);-(a2*c2+b2*d2)+(a2*d2-b2*c2)*im;-(a2*c2+b2*d2)+(b2*c2-a2*d2)*im;-(c2^2+d2^2)]
                supp[k+1]=[[[], []], [[vt;vt], [vt;vt]], [srt, [vt;vt]], [[vt;vt], srt], [srt, srt]]
                supp[k+1],coe[k+1]=move_zero!(supp[k+1],coe[k+1])
                dg[k-1:k]=[4;4]
                if normal==true
                    coe[k]=normalize(coe[k])
                    coe[k+1]=normalize(coe[k+1])
                end
                k+=2
            elseif LineLimit=="relax"
                mvr=ref[:bus][bus[vr]]["vmin"]^2
                mvt=ref[:bus][bus[vt]]["vmin"]^2
                coe[k]=[branch["rate_a"]^2*tm^4/mvr;-(a1^2+b1^2);-(a1*c1+b1*d1)+(b1*c1-a1*d1)*im;-(a1*c1+b1*d1)+(a1*d1-b1*c1)*im;-(c1^2+d1^2)]
                supp[k]=[[[], []], [[vr], [vr]], [[vr], [vt]], [[vt], [vr]], [[vt], [vt]]]
                supp[k],coe[k]=move_zero!(supp[k],coe[k])
                coe[k+1]=[branch["rate_a"]^2*tm^4/mvt;-(a2^2+b2^2);-(a2*c2+b2*d2)+(a2*d2-b2*c2)*im;-(a2*c2+b2*d2)+(b2*c2-a2*d2)*im;-(c2^2+d2^2)]
                supp[k+1]=[[[], []], [[vt], [vt]], [[vr], [vt]], [[vt], [vr]], [[vr], [vr]]]
                supp[k+1],coe[k+1]=move_zero!(supp[k+1],coe[k+1])
                if normal==true
                    coe[k]=normalize(coe[k])
                    coe[k+1]=normalize(coe[k+1])
                end
                k+=2
            end
        end
    end

    # power generation bound
    zero_pgen=UInt16[]
    for i=1:ng
        gen=ref[:gen][gens[i]]
        if gen["pmax"]>=1e-6
            coe[k]=[-4*gen["pmin"]*gen["pmax"];2*gen["pmin"]+2*gen["pmax"];2*gen["pmin"]+2*gen["pmax"];-1;-2;-1]
            if normal==true
                coe[k]=normalize(coe[k])
            end
            supp[k]=[[[], []], [[i+nbus], []], [[], [i+nbus]], [[i+nbus;i+nbus], []], [[i+nbus], [i+nbus]], [[], [i+nbus;i+nbus]]]
            supp[k],coe[k]=move_zero!(supp[k],coe[k])
            k+=1
        else
            push!(zero_pgen, i)
            numeq+=1
        end
        coe[k]=[-4*gen["qmin"]*gen["qmax"];-2*gen["qmin"]*im-2*gen["qmax"]*im;2*gen["qmin"]*im+2*gen["qmax"]*im;1;-2;1]
        if normal==true
            coe[k]=normalize(coe[k])
        end
        supp[k]=[[[], []], [[i+nbus], []], [[], [i+nbus]], [[i+nbus;i+nbus], []], [[i+nbus], [i+nbus]], [[], [i+nbus;i+nbus]]]
        supp[k],coe[k]=move_zero!(supp[k],coe[k])
        k+=1
    end

    # active/reactive power
    for (r, i) in enumerate(bus)
        bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]]
        bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]]
        coe[k]=zeros(ComplexF64, 2*length(ref[:bus_arcs][i])+2)
        supp[k]=Vector{Vector{Vector{UInt16}}}(undef, 2*length(ref[:bus_arcs][i])+2)
        supp[k][1:2]=[[[], []], [[r], [r]]]
        coe[k][1]=fl_sum(load["pd"] for load in bus_loads)
        coe[k][2]=fl_sum(shunt["gs"] for shunt in bus_shunts)
        coe[k+1]=zeros(ComplexF64, 2*length(ref[:bus_arcs][i])+2)
        supp[k+1]=Vector{Vector{Vector{UInt16}}}(undef, 2*length(ref[:bus_arcs][i])+2)
        supp[k+1][1:2]=[[[], []], [[r], [r]]]
        coe[k+1][1]=fl_sum(load["qd"] for load in bus_loads)
        coe[k+1][2]=-fl_sum(shunt["bs"] for shunt in bus_shunts)
        j=1
        for flow in ref[:bus_arcs][i]
            branch=ref[:branch][flow[1]]
            vr = bfind(bus,nbus,branch["f_bus"])
            vt = bfind(bus,nbus,branch["t_bus"])
            g, b = PowerModels.calc_branch_y(branch)
            tr, ti = PowerModels.calc_branch_t(branch)
            g_fr = branch["g_fr"]
            b_fr = branch["b_fr"]
            g_to = branch["g_to"]
            b_to = branch["b_to"]
            tm = branch["tap"]
            a1=(g+g_fr)/tm^2
            b1=-(b+b_fr)/tm^2
            c1=(-g*tr+b*ti)/tm^2
            d1=(b*tr+g*ti)/tm^2
            a2=g+g_to
            b2=-(b+b_to)
            c2=-(g*tr+b*ti)/tm^2
            d2=-(-b*tr+g*ti)/tm^2

            supp[k][j+2:j+3]=[[[vr], [vt]], [[vt], [vr]]]
            supp[k+1][j+2:j+3]=[[[vr], [vt]], [[vt], [vr]]]
            if vr==r
                coe[k][2]+=a1
                coe[k][j+2:j+3]=[(c1+d1*im)/2, (c1-d1*im)/2]
                coe[k+1][2]+=b1
                coe[k+1][j+2:j+3]=[(-c1*im+d1)/2, (c1*im+d1)/2]
            else
                coe[k][2]+=a2
                coe[k][j+2:j+3]=[(c2-d2*im)/2, (c2+d2*im)/2]
                coe[k+1][2]+=b2
                coe[k+1][j+2:j+3]=[(c2*im+d2)/2, (-c2*im+d2)/2]
            end
            j+=2
        end
        if !isempty(ref[:bus_gens][i])
            for gen_id in ref[:bus_gens][i]
                gen=bfind(gens, ng, gen_id)
                push!(supp[k], [[gen+nbus], []], [[], [gen+nbus]])
                push!(coe[k], -0.5, -0.5)
                push!(supp[k+1], [[gen+nbus], []], [[], [gen+nbus]])
                push!(coe[k+1], 0.5*im, -0.5*im)
            end
        end
        if normal==true
            coe[k]=normalize(coe[k])
            coe[k+1]=normalize(coe[k+1])
        end
        supp[k],coe[k]=move_zero!(supp[k],coe[k])
        supp[k+1],coe[k+1]=move_zero!(supp[k+1],coe[k+1])
        k+=2
    end

    # reference voltage
    for key in keys(ref[:ref_buses])
        i=bfind(bus,nbus,key)
        supp[k]=[[[i;i], []], [[i], [i]], [[], [i;i]]]
        coe[k]=[1;-2;1]
        k+=1
    end

    # zero power generation
    for i in zero_pgen
        supp[k]=[[[i+nbus], []], [[], [i+nbus]]]
        coe[k]=[1;1]
        dg[k-1]=1
        k+=1
    end

    return SparsePolyModel(n,m,numeq,nbus,ng,nb,supp,coe,dg)
end


text:
struct VerySimpleRaster{T}
    mat::Matrix{T}
    nodata::T
    xs::StepRangeLen{Float64}
    ys::StepRangeLen{Float64}
    projection::String
    filename::String
end

VerySimpleRaster(mat::AbstractMatrix{T}, nodata, xll, yll, cellx, celly = cellx, proj = "", fname = "") where T =
    VerySimpleRaster{T}(mat, convert(T, nodata), range(xll, step = cellx, length = size(mat,1)), range(yll, step = celly, length = size(mat,2)), proj, abspath(fname))

bbox(vsr::VerySimpleRaster) = (
    xmin = minimum(vsr.xs),
    xmax = maximum(vsr.xs) + step(vsr.xs),
    ymin = minimum(vsr.ys),
    ymax = maximum(vsr.ys) + step(vsr.ys))

Base.show(io::IO, vsr::VerySimpleRaster{T}) where T = println(io,
"""
VerySimpleRaster{$T} with $(length(vsr.ys)) rows and $(length(vsr.xs)) columns
Resolution: $(round(step(vsr.xs), sigdigits = 5)), $(round(step(vsr.ys), sigdigits = 5))
Extent: $(round(bbox(vsr).xmin, sigdigits = 5)) to $(round(bbox(vsr).xmax, sigdigits = 5)) longitude, $(round(bbox(vsr).ymin, sigdigits = 5)) to $(round(bbox(vsr).ymax, sigdigits = 5)) latitude
Projection: $(vsr.projection)
Source file: $(vsr.filename)""")


text:
@recipe(Errorbars) do scene
    Theme(
        whiskerwidth = 10,
        color = :black,
        linewidth = 1,
        direction = :y,
        visible = theme(scene, :visible)
    )
end

function AbstractPlotting.plot!(plot::Errorbars{T}) where T <: Tuple{<:Any, <:Any, <:Any, <:Any}

    xs, ys, low, high = plot[1:4]
    xys = lift((x, y) -> Point2.(x, y), xs, ys)
    _plot_errorbars!(plot, xys, low, high)
end

function AbstractPlotting.plot!(plot::Errorbars{T}) where T <: Tuple{<:Any, <:Any, <:Any}
    xys, low, high = plot[1:3]
    _plot_errorbars!(plot, xys, low, high)
end

function AbstractPlotting.plot!(plot::Errorbars{T}) where T <: Tuple{<:Any, <:Any}
    xys, lowhigh = plot[1:2]
    _plot_errorbars!(plot, xys, lowhigh, lowhigh)
end

f_if(condition, f, arg) = condition ? f(arg) : arg

function _plot_errorbars!(plot, xys, low, high)

    @extract plot (whiskerwidth, color, linewidth, direction, visible)

    is_in_y_direction = lift(direction) do dir
        if dir == :y
            true
        elseif dir == :x
            false
        else
            error("Invalid direction $dir. Options are :x and :y.")
        end
    end

    linesegpairs = lift(xys, low, high, is_in_y_direction) do xys, low, high, is_in_y_direction

        [(xy .+ f_if(is_in_y_direction, reverse, Point2f0(-lo, 0)),
          xy .+ f_if(is_in_y_direction, reverse, Point2f0( hi, 0)),)
                for (xy, lo, hi) in zip(xys, low, high)]
    end

    scene = plot.parent

    whiskers = lift(linesegpairs, scene.camera.projectionview,
        scene.camera.pixel_space, whiskerwidth) do pairs, _, _, whiskerwidth

        endpoints = [p for pair in pairs for p in pair]

        screenendpoints = scene_to_screen(endpoints, scene)

        screenendpoints_shifted_pairs = map(screenendpoints) do sep
            (sep .+ f_if(is_in_y_direction[], reverse, Point(0, -whiskerwidth/2)),
             sep .+ f_if(is_in_y_direction[], reverse, Point(0,  whiskerwidth/2)))
        end

        screen_to_scene([p for pair in screenendpoints_shifted_pairs for p in pair], scene)
    end

    linesegments!(plot, linesegpairs, color = color, linewidth = linewidth, visible = visible)
    linesegments!(plot, whiskers, color = color, linewidth = linewidth, visible = visible)
    plot
end

function scene_to_screen(pts, scene)
    p4 = to_ndim.(Vec4f0, to_ndim.(Vec3f0, pts, 0.0), 1.0)
    p1m1 = Ref(scene.camera.projectionview[]) .* p4
    projected = Ref(inv(scene.camera.pixel_space[])) .* p1m1
    [Point2.(p[1:2]...) for p in projected]
end

function screen_to_scene(pts, scene)
    p4 = to_ndim.(Vec4f0, to_ndim.(Vec3f0, pts, 0.0), 1.0)
    p1m1 = Ref(scene.camera.pixel_space[]) .* p4
    projected = Ref(inv(scene.camera.projectionview[])) .* p1m1
    [Point2.(p[1:2]...) for p in projected]
end

function scene_to_screen(p::T, scene) where T <: Point
    p4 = to_ndim(Vec4f0, to_ndim(Vec3f0, p, 0.0), 1.0)
    p1m1 = scene.camera.projectionview[] * p4
    projected = inv(scene.camera.pixel_space[]) * p1m1
    T(projected[1:2]...)
end

function screen_to_scene(p::T, scene) where T <: Point
    p4 = to_ndim(Vec4f0, to_ndim(Vec3f0, p, 0.0), 1.0)
    p1m1 = scene.camera.pixel_space[] * p4
    projected = inv(scene.camera.projectionview[]) * p1m1
    T(projected[1:2]...)
end


text:
using TransformUtils
const TU = TransformUtils

# add AprilTag sightings from this pose
function addApriltags!(fg, pssym, posetags; bnoise=0.1, rnoise=0.1, lmtype=Point2, fcttype=Pose2Pose2, DAerrors=0.0, autoinit=true )
  @show currtags = ls(fg, r"l")
  for lmid in keys(posetags)
    @show lmsym = Symbol("l$lmid")
    if !(lmsym in currtags)
      @info "adding node $lmsym"
      addVariable!(fg, lmsym, lmtype, timestamp=getTimestamp(getVariable(fg, pssym)))
    end
    ppbr = nothing
    if lmtype == RoME.Point2
      ppbr = Pose2Point2BearingRange(Normal(posetags[lmid][:bearing][1],bnoise),
                                     Normal(posetags[lmid][:range][1],rnoise))
    elseif lmtype == Pose2
      dx, dy = posetags[lmid][:bP2t].translation[1], posetags[lmid][:bP2t].translation[2]
      dth = convert(RotXYZ, posetags[lmid][:bP2t].linear).theta3
      # ppbr = Pose2Pose2(
      #           MvNormal([ posetags[lmid][:pos][3];
      #                     -posetags[lmid][:pos][1];
      #                      posetags[lmid][:tRYc]],
      #                    diagm([0.1;0.1;0.05].^2)) )
      ppbr = fcttype(
                MvNormal([dx;
                          dy;
                          dth],
                         Matrix(Diagonal([0.05;0.05;0.03].^2))) )
    end
    if rand() > DAerrors
      # regular single hypothesis
      addFactor!(fg, [pssym; lmsym], ppbr, autoinit=autoinit)
    else
      # artificial errors to data association occur
      info("Forcing bad data association with $lmsym")
      xx,ll = ls(fg)
      @show ll2 = setdiff(ll, [lmsym])
      @show daidx = round(Int, (length(ll2)-1)*rand()+1)
      @show rda = ll2[daidx]
      addFactor!(fg, [pssym; lmsym; rda], ppbr, multihypo=[1.0;0.5;0.5], autoinit=autoinit)
    end
  end
  nothing
end

function addnextpose!(fg,
                      prev_psid,
                      new_psid,
                      pose_tag_bag,
                      poseTimes;
                      lmtype=Point2,
                      odotype=Pose2Pose2,
                      fcttype=Pose2Pose2,
                      DAerrors=0.0,
                      autoinit=true,
                      odopredfnc=nothing,
                      joyvel=nothing,
                      naiveFrac=0.6,
                      DXmvn = MvNormal(zeros(3),diagm([0.4;0.1;0.4].^2)) )
  #
  prev_pssym = Symbol("x$(prev_psid)")
  new_pssym = Symbol("x$(new_psid)")
  #naive odometry model

  donnpose = false
  if odopredfnc != nothing && joyvel != nothing
    donnpose = true
  end

  # let odoKDE,
    # # predict delta x y th odo if able
    #   nnpts = odopredfnc(joyvel)
    #   # replace theta points
    #   nnpts[:,3] .= rand(DXmvn, size(nnpts,1))[:,3]
    #   mvnpts = rand( DXmvn, round(Int, parametricOdoMix*size(nnpts, 1)) )
    #   odoKDE = manikde!([nnpts;mvnpts], Pose2)
    # else
    #   odoKDE = DXmvn

    # first pose with zero prior
    if odotype == Pose2Pose2
      addVariable!(fg, new_pssym, Pose2, timestamp=poseTimes[new_pssym])
      joyVals = zeros(25,4)
      for i in 1:25
        joyVals[i,:] .= joyvel[prev_pssym][i]
      end
      pp = !donnpose ? Pose2Pose2(DXmvn) : FluxModelsPose2Pose2(odopredfnc,joyVals,DXmvn,naiveFrac) # PyNeuralPose2Pose2(odopredfnc,joyvel[prev_pssym],DXmvn,naiveFrac)
      addFactor!(fg, [prev_pssym; new_pssym], pp, graphinit=autoinit)
    elseif odotype == VelPose2VelPose2
      donnpose ? error("Not implemented for VelPose2VelPose2 yet") : nothing
      addVariable!(fg, new_pssym, DynPose2(ut=round(Int, 200_000*(new_psid))), timestamp=poseTimes[new_pssym] )
      vpvp = VelPose2VelPose2(odoKDE, MvNormal(zeros(2),Matrix(Diagonal([0.2;0.1].^2))))
      addFactor!(fg, [prev_pssym; new_pssym], vpvp, graphinit=autoinit)
      #
    end
  # end

  addApriltags!(fg, new_pssym, pose_tag_bag, lmtype=lmtype, fcttype=fcttype, DAerrors=DAerrors)
  new_pssym
end

function prepCamLookup(imgseq; filenamelead="camera_image")
  # camcount = readdlm(datafolder*"$(imgfolder)/cam-count.csv",',')
  camlookup = Dict{Int, String}()
  count = -1
  for i in imgseq
    count += 1
    camlookup[count] = "$(filenamelead)$(i).jpeg"
    # camlookup[camcount[i,1]] = strip(camcount[i,2])
  end
  return camlookup
end

function detectTagsViaCamLookup(camlookup, imgfolder, imgsavedir)
  # AprilTag detector
  detector = AprilTagDetector()

  # extract tags from images
  IMGS = []
  TAGS = []
  psid = 1
  for psid in 0:(length(camlookup)-1)
    img = load("$(imgfolder)/$(camlookup[psid])")
    tags = detector(img)
    push!(TAGS, deepcopy(tags))
    push!(IMGS, deepcopy(img))
    foreach(tag->drawTagBox!(IMGS[psid+1],tag, width = 5, drawReticle = false), tags)
    Images.save(imgsavedir*"/tags/img_$(psid).jpg", IMGS[psid+1])
  end
  # psid = 1
  # img = load(datafolder*"$(imgfolder)/$(camlookup[psid])")

  # free the detector memory
  freeDetector!(detector)

  return IMGS, TAGS
end

function results2csv(fg; dir=".", filename="results.csv")

if !isdir(dir*"/results")
  mkdir(dir*"/results")
end
fid = open(dir*"/results/"*filename,"w")
prevpos = [0.0;0.0;0.0]
prevT = getTimestamp(getVariable(fg, :x0))
poses = ls(fg, r"x\d") |> sortDFG
for sym in poses # poses first
  # p = getKDE(fg, sym)
  # val = string(KDE.getKDEMax(p))
  val = getPPE(fg, sym).suggested
  # assuming Pose2 was used
  newT = getTimestamp(getVariable(fg, sym))
  DT = (newT-prevT).value*1e-3
  velx = (val[1] - prevpos[1])./DT
  vely = (val[2] - prevpos[2])./DT
  wVelx = isnan(velx) ? 0.0 : velx
  wVely = isnan(vely) ? 0.0 : vely
  # rotate velocity from world to body frame
  bVel = TU.R(-val[3])*[wVelx;wVely]
  velx = bVel[1]
  vely = bVel[2]
  if sym in [:x0;:x1]
    velx = 0.0
    vely = 0.0
  end
  prevT = newT
  prevpos = val
  println(fid, "$sym, $(val[1]), $(val[2]), $(val[3]), $velx, $vely")
end
close(fid)

end

## BUILD FACTOR GRAPH FOR SLAM SOLUTION,

function main(WP,
              resultsdir::String,
              camidxs,
              tag_bagl;
              maxlen = (length(tag_bagl)-1),
              BB=10,
              N=100,
              lagLength=100,
              dofixedlag=true,
              jldfile::String="",
              failsafe::Bool=false,
              show::Bool=false,
              odopredfnc=nothing,
              joyvel=nothing,
              poseTimes=nothing,
              multiproc=true,
              drawtree=true,
              batchResolve=true)
#

# Factor graph construction
fg = initfg()
getSolverParams(fg).logpath = resultsdir
getSolverParams(fg).multiproc=multiproc
getSolverParams(fg).drawtree=drawtree
getSolverParams(fg).showtree=false
getSolverParams(fg).maxincidence = 500
prev_psid = 0

# load from previous file
if jldfile != ""
  # fg, = loadjld( file=joinpath(resultsdir,jldfile) )
  fg = initfg()
  loadDFG(joinpath(resultsdir,jldfile), Main, fg)
  xx = ls(fg, r"x")
  prev_psid = parse(Int, string(xx[end])[2:end])
end

# fg.solverParams.isfixedlag = dofixedlag
# fg.solverParams.qfl = lagLength
defaultFixedLagOnTree!(fg, lagLength, limitfixeddown=false)

psid = 0
pssym = Symbol("x$psid")
# first pose with zero prior
addVariable!(fg, pssym, Pose2, timestamp=poseTimes[:x0])
# addFactor!(fg, [pssym], PriorPose2(MvNormal(zeros(3),diagm([0.01;0.01;0.001].^2))))
addFactor!(fg, [pssym], PriorPose2(MvNormal(zeros(3),Matrix(Diagonal([0.01;0.01;0.001].^2)))) )
#
# addVariable!(fg, pssym, DynPose2(ut=0))
# # addFactor!(fg, [pssym], PriorPose2(MvNormal(zeros(3),diagm([0.01;0.01;0.001].^2))))
# addFactor!(fg, [pssym], DynPose2VelocityPrior(MvNormal(zeros(3),Matrix(Diagonal([0.01;0.01;0.001].^2))),
#                                               MvNormal(zeros(2),Matrix(Diagonal([0.1;0.05].^2)))))
#

addApriltags!(fg, pssym, tag_bagl[psid], lmtype=Pose2, fcttype=Pose2Pose2) # DynPose2Pose2

show ? drawGraph(fg, show=true) : nothing

# quick solve as sanity check
failsafe ? @warn("recursive failsafe no longer as option") : nothing
# , N=N, drawpdf=true, show=show, recursive=failsafe
tree, smt, hist = solveTree!(fg)

# add other positions

# Gadfly.push_theme(:default)

# @sync begin

for psid in (prev_psid+1):1:maxlen
  # global prev_psid, maxlen
  @show psym = Symbol("x$psid")
  addnextpose!(fg, prev_psid, psid, tag_bagl[psid], poseTimes, lmtype=Pose2, odotype=Pose2Pose2, fcttype=Pose2Pose2, autoinit=true, odopredfnc=odopredfnc, joyvel=joyvel)
  # odotype=VelPose2VelPose2, fcttype=DynPose2Pose2
  # writeGraphPdf(fg)

  if psid % BB == 0 || psid == maxlen
    saveDFG(fg, resultsdir*"/racecar_fg_$(psym)_presolve")
    # , drawpdf=true, show=show, N=N, recursive=true
    tree, smt, hist = solveTree!(fg, tree)
  end

  # T1 = remotecall(saveDFG, WP, fg, resultsdir*"/racecar_fg_$(psym)")
  # @async fetch(T1)
  saveDFG(fg, resultsdir*"/racecar_fg_$(psym)")

  ## save factor graph for later testing and evaluation
  ensureAllInitialized!(fg)
  # T2 = remotecall(plotRacecarInterm, WP, fg, resultsdir, psym)
  # @async fetch(T2)
  if 1 < nprocs()
    remotecall(plotRacecarInterm, WP, fg, resultsdir, psym)
  else
    plotRacecarInterm(fg, resultsdir, psym)
  end

  # prepare for next iteration
  prev_psid = psid
end # for

# do a batch resolve first
if batchResolve
  dontMarginalizeVariablesAll!(fg)
  foreach(x->setSolvable!(fg, x, 1), setdiff(ls(fg),ls(fg,r"drt")))
  foreach(x->setSolvable!(fg, x, 1), setdiff(lsf(fg),lsf(fg,r"drt")))
  tree = solveTree!(fg)
end
saveDFG(fg, resultsdir*"/racecar_fg_final")

# extract results for later use as training data
results2csv(fg, dir=resultsdir, filename="results.csv")

# end #sync

return fg
end # main

#