text:
using MakieGallery: @substep

@block SimonDanisch ["tests"] begin
    ## Scatter
    @cell begin
        using GeometryBasics
        ## Some helpers
        data_2d = AbstractPlotting.peaks()
        args_2d = (-10..10, -10..10, data_2d)
        function n_times(f, n=10, interval=0.05)
            obs = Observable(f(1))
            @async for i in 2:n
                try
                    obs[] = f(i)
                    sleep(interval)
                catch e
                    @warn "Error!" exception=CapturedException(e, Base.catch_backtrace())
                end
            end
            return obs
        end
        cmap = cgrad(:RdYlBu; categorical=true);
        line_positions = Point2f0.([1:4; NaN; 1:4], [1:4; NaN; 2:5])
        cat = MakieGallery.loadasset("cat.obj")
        tex = MakieGallery.loadasset("diffusemap.png")
        meshes = GeometryBasics.normal_mesh.([Sphere(Point3f0(0.5), 1), Rect(Vec3f0(1, 0, 0), Vec3f0(1))])
        positions = Point2f0[(1, 1), (2, 2), (3, 2), (4, 4)]
        positions_connected = connect(positions, LineFace{Int}[(1, 2), (2, 3), (3, 4)])

        scatter(1:4; color=:red, markersize=0.3)
        @substep
        scatter(1:4; color=:red, markersize=10px)
        @substep
        scatter(1:4; color=:red, markersize=10, markerspace=Pixel)
        @substep
        scatter(1:4; color=:red, markersize=(1:4).*8, markerspace=Pixel)

        @substep
        scatter(1:4; marker='☼')
        @substep
        scatter(1:4; marker=['☼', '◒', '◑', '◐'])
        @substep
        scatter(1:4; marker="☼◒◑◐")
        @substep
        scatter(1:4; marker=rand(RGBf0, 10, 10), markersize=20px)
        # TODO rotation with markersize=px
        @substep
        scatter(1:4; marker='▲', markersize=0.3, rotations=LinRange(0, pi, 4))

        ## Meshscatter
        @substep
        meshscatter(1:4; color=1:4)

        @substep
        meshscatter(1:4; color=rand(RGBAf0, 4))
        @substep
        meshscatter(1:4; color=rand(RGBf0, 4))
        @substep
        meshscatter(1:4; color=:red)
        @substep
        meshscatter(rand(Point3f0, 10); color=rand(RGBf0, 10))
        @substep
        meshscatter(rand(Point3f0, 10); marker=Pyramid(Point3f0(0), 1f0, 1f0))

        ## Barplot
        barplot(sort(rand(10)); color=rand(10))
        @substep
        barplot(1:3; color=[:red, :green, :blue])
        @substep

        ## Lines
        @substep
        lines(line_positions)
        @substep
        lines(line_positions; linestyle=:dot)
        @substep
        lines(line_positions; linestyle=[0.0, 1.0, 2.0, 3.0, 4.0])
        @substep
        lines(line_positions; color=1:9)
        @substep
        lines(line_positions; color=rand(RGBf0, 9), linewidth=4)

        ## Linesegments
        @substep
        linesegments(1:4)
        @substep
        linesegments(1:4; linestyle=:dot)
        @substep
        linesegments(1:4; linestyle=[0.0, 1.0, 2.0, 3.0, 4.0])
        @substep
        linesegments(1:4; color=1:4)
        @substep
        linesegments(1:4; color=rand(RGBf0, 4), linewidth=4)

        ## Surface
        @substep
        surface(args_2d...)
        @substep
        surface(args_2d...; color=rand(size(data_2d)...))
        @substep
        surface(args_2d...; color=rand(RGBf0, size(data_2d)...))
        @substep
        surface(args_2d...; colormap=:magma, colorrange=(-3.0, 4.0))
        @substep
        surface(args_2d...; shading=false); wireframe!(args_2d..., linewidth=0.5)
        @substep
        surface(1:30, 1:31, rand(30, 31))
        @substep
        n = 20
        θ = [0;(0.5:n-0.5)/n;1]
        φ = [(0:2n-2)*2/(2n-1);2]
        x = [cospi(φ)*sinpi(θ) for θ in θ, φ in φ]
        y = [sinpi(φ)*sinpi(θ) for θ in θ, φ in φ]
        z = [cospi(θ) for θ in θ, φ in φ]
        surface(x, y, z)
        ## Polygons
        @substep
        poly(decompose(Point2f0, Circle(Point2f0(0), 1f0)))

        ## Image like!
        @substep
        image(rand(10, 10))
        @substep
        heatmap(rand(10, 10))

        ## Volumes
        @substep
        volume(rand(4, 4, 4), isovalue=0.5, isorange=0.01, algorithm=:iso)
        @substep
        volume(rand(4, 4, 4), algorithm=:mip)
        @substep
        volume(rand(4, 4, 4), algorithm=:absorption)
        @substep
        volume(rand(4, 4, 4), algorithm=Int32(5))

        @substep
        volume(rand(RGBAf0, 4, 4, 4), algorithm=:absorptionrgba)
        @substep
        contour(rand(4, 4, 4))

        ## Meshes

        @substep
        mesh(cat, color=tex)
        @substep
        mesh([(0.0, 0.0), (0.5, 1.0), (1.0, 0.0)], color = [:red, :green, :blue],
             shading = false)
        @substep
        mesh(meshes, color=[1, 2])
        @substep
        t = Tesselation(FRect2D(0, 0, 2, 2), (30, 30))
        m = GeometryBasics.mesh(t, pointtype=Point3f0, facetype=QuadFace{GLIndex})
        s = AbstractPlotting.mesh(m)
        wireframe!(s, m, color=("#ddd", 0.8), transparency=true, linewidth=2)
        scatter!(s, m, markersize=5px, color=:red)
        @substep

        ## Axis
        scene = lines(IRect(Vec2f0(0), Vec2f0(1)))
        axis = scene[Axis]
        axis.ticks.ranges = ([0.1, 0.2, 0.9], [0.1, 0.2, 0.9])
        axis.ticks.labels = (["😸", "♡", "𝕴"], ["β ÷ δ", "22", "≙"])
        scene
        @substep
        ## Text
        text("heyllo")
        @substep
        ## Colormaps
        heatmap(args_2d...; colormap=cgrad(:RdYlBu; categorical=true), interpolate=true)
        @substep
        image(args_2d...; colormap=cgrad(:RdYlBu; categorical=true), interpolate=true)
        @substep
        heatmap(args_2d...; colorrange=(-4.0, 3.0), colormap=cmap, highclip=:green,
                      lowclip=:pink, interpolate=true)
        @substep
        surface(args_2d...; colorrange=(-4.0, 3.0), highclip=:green, lowclip=:pink)
        @substep

        ## Animations
        annotations(n_times(i-> map(j-> ("$j", Point2f0(j*30, 0)), 1:i)), textsize=20,
                          limits=FRect2D(30, 0, 320, 50))
        @substep
        scatter(n_times(i-> Point2f0.((1:i).*30, 0)), markersize=20px,
                      limits=FRect2D(30, 0, 320, 50))
        @substep
        linesegments(n_times(i-> Point2f0.((2:2:2i).*30, 0)), limits=FRect2D(30, 0, 620, 50))
        @substep
        lines(n_times(i-> Point2f0.((2:2:2i).*30, 0)), limits=FRect2D(30, 0, 620, 50))
        @substep
        # Scaling
        scene = Scene(transform_func=(identity, log10))
        linesegments!(1:4, color=:black, linewidth=20, transparency=true)
        scatter!(1:4, color=rand(RGBf0, 4), markersize=20px)
        lines!(1:4, color=rand(RGBf0, 4))
        @substep

        # Views
        linesegments(positions)
        @substep
        scatter(positions_connected)
        @substep

        # Interpolation
        heatmap(data_2d, interpolate=true)
        @substep
        image(data_2d, interpolate=false)
        @substep

        # Categorical
        barplot(["hi", "ima", "string"], rand(3))
        @substep
        heatmap(["a", "b", "c"], ["α", "β", "γ"], rand(3, 3))

    end
end


text:
# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools
using Optim

x = -10:0.01:10
len_x = length(x)

options = Optim.Options(
    g_abstol=1e-8, g_reltol=1e-8, iterations=10^5, show_trace=false, store_trace=true,
);

let
    # Gaussian 1 peak
    seed1 = [[0.0], [15.0, 0.0, sqrt(0.2)]]
    # We generate some data using the Gaussian model
    y_data = Utils.gaussian(x, seed1) .+ 0.5.*randn(len_x)
    sol = fit_curve_model(:gaussian, x, y_data, seed1)
    plot(x, [y_data, sol.ymodel], lab=["Data" "Model"])
    gui()
    # Gaussian 2 peaks
    seed2 = [[1.0], [15.0, -1.0, sqrt(0.2)], [10.0, 1.5, sqrt(0.1)]]
    # We generate some data using the Gaussian model
    y_data = Utils.gaussian(x, seed2) .+ 0.5.*randn(len_x)
    sol = fit_curve_model(:gaussian, x, y_data, seed2; options=options)
    plot(x, [y_data, sol.ymodel], lab=["Data" "Model"])
    gui()
end

let
    # Lorentzian 1 peak
    seed1 = [[0.0], [15.0, 0.0, 0.5]]
    # We generate some data using the Lorentzian model
    y_data = Utils.lorentzian(x, seed1) .+ 0.5.*randn(len_x)
    sol = fit_curve_model(:lorentzian, x, y_data, seed1; options=options)
    plot(x, [y_data, sol.ymodel], lab=["Data" "Model"])
    gui()
    # Lorentzian 2 peaks
    seed2 = [[0.0], [15.0, -1.0, 0.5], [30.0, 1.0, 0.5]]
    # We generate some data using the Lorentzian model
    y_data = Utils.lorentzian(x, seed2) .+ 0.5.*randn(len_x)
    sol = fit_curve_model(:lorentzian, x, y_data, seed2; options=options)
    plot(x, [y_data, sol.ymodel], lab=["Data" "Model"])
    gui()
end

let
    # Voigtian 1 peak
    seed1 = [[0.0], [18.0, 0.0, 0.0, 1.53]]
    # We generate some data using the Voigtian model
    y_data = Utils.voigtian(x, seed1) .+ randn(len_x)
    sol = fit_curve_model(:voigtian, x, y_data, seed1; options=options)
    plot(x, [y_data, sol.ymodel], lab=["Data" "Model"])
    gui()
    # Voigtian 2 peaks
    seed2 = [[0.0], [10.0, -1.6, 0.0, 1.53], [24.0, 5.0, 0.0, 1.53]]
    # We generate some data using the Voigtian model
    y_data = Utils.voigtian(x, seed2) .+ randn(len_x)
    sol = fit_curve_model(:voigtian, x, y_data, seed2; options=options)
    plot(x, [y_data, sol.ymodel], lab=["Data" "Model"])
    gui()
end

ext: jl


text:
using Makie
AbstractPlotting.inline!(false)

O3d = [0,0,0]
A3d = [8,2,7]
B3d = [-2,-8,6]
C3d = [-3,6,-9]

limits = FRect3D(Vec3f0(-10, -10, -10), Vec3f0(20, 20, 20))

scene = Scene(resolution = (1200,1200))

AbstractPlotting.meshscatter!(scene, [O3d[1]],[O3d[2]],[O3d[3]], color = :dodgerblue,
    limits = limits, markersize = 0.5)
axis = scene[Axis]
axis.names.axisnames = ("posição x (m)", "posição y (m)", "posição z (m)")
axis.names.textsize = 5

AbstractPlotting.meshscatter!([A3d[1]],[A3d[2]],[A3d[3]], color = :orange, markersize = 0.5, transparency = true)
arrows!([0],[0],[0],[A3d[1]],[A3d[2]],[A3d[3]], arrowcolor = :orange, linecolor = :orange, linewidth = 5,
    arrowsize = 0.5)
arrows!([0,A3d[1],A3d[1]],[0,0,A3d[2]],[0,0,0],[A3d[1],0,0],[0,A3d[2],0],[0,0,A3d[3]],
    arrowcolor = :darkorange, linecolor = :orange, linewidth = 5, linestyle = :dash, arrowsize = 0.5)

AbstractPlotting.meshscatter!([B3d[1]],[B3d[2]],[B3d[3]], color = :red, markersize = 0.5, transparency = true)
arrows!([0],[0],[0],[B3d[1]],[B3d[2]],[B3d[3]], arrowcolor = :red, linecolor = :red, linewidth = 5,
    arrowsize = 0.5)
arrows!([0,B3d[1],B3d[1]],[0,0,B3d[2]],[0,0,0],[B3d[1],0,0],[0,B3d[2],0],[0,0,B3d[3]],
    arrowcolor = :darkred, linecolor = :red, linewidth = 5, linestyle = :dash, arrowsize = 0.5)

AbstractPlotting.meshscatter!([C3d[1]],[C3d[2]],[C3d[3]], color = :limegreen, markersize = 0.5, transparency = true)
arrows!([0],[0],[0],[C3d[1]],[C3d[2]],[C3d[3]], arrowcolor = :limegreen, linecolor = :limegreen, linewidth = 5,
    arrowsize = 0.5)
arrows!([0,0,C3d[1]],[0,C3d[2],C3d[2]],[0,0,0],[0,C3d[1],0],[C3d[2],0,0],[0,0,C3d[3]],
    arrowcolor = :forestgreen, linecolor = :limegreen, linewidth = 5, linestyle = :dash, arrowsize = 0.5)

scene

max_stars_repo_name: ghralima/mecanica

max_issues_repo_name: ghralima/mecanica

max_forks_repo_name: ghralima/mecanica


text:
CLASS_PARAMS = """    outlier_fraction::Float64
The fraction of outliers (number between 0 and 1) in the data is used to determine the score threshold to classify the
samples into inliers and outliers.

    classify::Union{Function, Nothing}
A function to transform a vector of scores to a vector of classes, where `"outlier"` represents an outlier and
`"normal"` represents a normal instance. See [`classify`](@ref) for a specific implementation."""

SCORE_PARAMS = """    combine::Function
A function to reduce a matrix, where each row represents an instance and each column represents a score of specific
detector, to a vector of scores for each instance. See `combine` for a specific implementation. *Note:* This function
is not called if the input to the evaluator consists of a single train/test scores tuple.

    normalize::Union{Function, Nothing}
A function to reduce a matrix, where each row represents an instance and each column a score of specific detector, to a
vector of scores for each instance. See [`normalize`](@ref) for a specific implementation."""

"""
    ScoreTransformer(combine = combine,
                     normalize = normalize)

Transform the results of a single or multiple outlier detection models to combined and normalized scores.

Parameters
----------
$SCORE_PARAMS
```
"""
@detector mutable struct ScoreTransformer <: MLJ.Static
    normalize::Function = scale_minmax
    combine::Function = combine_mean
end

"""
    ClassTransformer(normalize = scale_minmax,
                     combine = combine_mean,
                     classify = classify_percentile(0.9))

Transform the results of a single or multiple outlier detection models to binary classes, where `"normal"` represents
inliers and `"outlier"` represents outliers.

Parameters
----------
$CLASS_PARAMS

$SCORE_PARAMS
"""
default_percentile_threshold = classify_percentile(DEFAULT_THRESHOLD)
@detector mutable struct ClassTransformer <: MLJ.Static
    normalize::Function = scale_minmax
    combine::Function = combine_mean
    classify::Function = default_percentile_threshold
end

function MLJ.transform(ev::ScoreTransformer, _, scores::Tuple{Scores, Scores}...) # _ because there is no fitresult
    _, scores_test = to_scores(ev.normalize, ev.combine, scores...)
    to_univariate_finite(scores_test)
end

function MLJ.transform(ev::ClassTransformer, _, scores::Tuple{Scores, Scores}...) # _ because there is no fitresult
    _, classes_test = to_classes(ev.normalize, ev.combine, ev.classify, scores...)
    to_categorical(classes_test)
end

function to_scores(normalize::Function, combine::Function, scores::Tuple{Scores, Scores}...)::Tuple{Scores, Scores}
    # make sure that we have at least one tuple of train and test scores
    n_scores = length(scores)
    @assert n_scores > 0

    # normalize all scores (might be identity function)
    scores = map(score_tuple -> normalize(score_tuple...), scores)

    # [(train1, test1), (train2, test2)] -> matrix of scores for train [train1 train2]' and test [test1 test2]'
    # Note: The matrices contain one observation per row!
    reduce_cat = idx -> reduce(hcat, getfield.(scores, idx))
    scores_train, scores_test = combine(reduce_cat(1)), combine(reduce_cat(2))
    scores_train, scores_test
end

function to_classes(normalize::Function, combine::Function, classify::Function,
                    scores::Tuple{Scores, Scores}...)::Tuple{Labels, Labels}
    scores_train, scores_test = to_scores(normalize, combine, scores...)
    classify(scores_train, scores_test)
end


text:
# Functions for estimating probabilities using various estimators

# Estimators are described in:
# Hausser, Jean; Strimmer, Korbinian (2009-01-01).
# "Entropy Inference and the James-Stein Estimator, with Application to Nonlinear Gene Association Networks"
# https://arxiv.org/abs/0811.3579

# R implementation of estimators:
# https://cran.r-project.org/web/packages/entropy/

# Parameters:
# 	- frequencies, integer array
# 	- prior, number
function get_probabilities_dirichlet(frequencies, prior)
	prior_frequencies = fill(prior, size(frequencies))
	return (frequencies + prior_frequencies) / (sum(frequencies) + sum(prior_frequencies))
end

# Parameters:
# 	- frequencies, integer array
function get_probabilities_maximum_likelihood(frequencies)
	return frequencies / sum(frequencies)
end

# TODO: allow pre-calculated target to be passed in instead of uniform distribution
# Parameters:
# 	- frequencies, integer array
# 	- lambda, void
function get_probabilities_shrinkage(frequencies, lambda::Nothing)
	target = get_uniform_distribution(frequencies)
	n = sum(frequencies)
	normalized_frequencies = frequencies / n
	calculated_lambda = get_lambda(normalized_frequencies, target, n)
	return apply_shrinkage_formula(normalized_frequencies, target, calculated_lambda)
end
# Parameters:
# 	- frequencies, integer array
# 	- lambda, number
function get_probabilities_shrinkage(frequencies, lambda)
	target = get_uniform_distribution(frequencies)
	normalized_frequencies = get_normalized_frequencies(frequencies)
	return apply_shrinkage_formula(normalized_frequencies, target, lambda)
end

function apply_shrinkage_formula(normalized_frequencies, target, lambda)
	return lambda * target .+ (1 - lambda) * normalized_frequencies
end

function get_uniform_distribution(frequencies)
	return 1 / length(frequencies)
end

function get_normalized_frequencies(frequencies)
	return frequencies / sum(frequencies)
end

function get_lambda(normalized_frequencies, target, n)
	if n == 1 || n == 0
		return 1.0
	end
	# Add better comments about varu and msp
	# Unbiased estimator of variance of u
	varu = normalized_frequencies .* (1 .- normalized_frequencies) / (n - 1)
	msp = sum((normalized_frequencies .- target).^2) # misspecification ???

	# Estimate shrinkage intensity
	lambda = msp == 0 ? 1.0 : sum(varu) / msp

	# Make lambda be between 0 and 1 inclusive
	return lambda > 1 ? 1.0 : (lambda < 0 ? 0.0 : lambda)

end
function get_lambda(frequencies, get_target = get_uniform_distribution)
	n = sum(frequencies)
	if n == 1 || n == 0
		return 1
	end

	target = get_target(frequencies)
	normalized_frequencies = frequencies / n

	return get_lambda(normalized_frequencies, target, n)
end


text:
using Plots; pyplot()
using LinearAlgebra

################################################################################
################################################################################
################################################################################
# Optimization Algorithms
################################################################################
################################################################################
################################################################################

################################################################################
# A simple armijo line search with backtracking strategy
################################################################################
function armijo_bt(x::Array{Float64,1},p::Array{Float64,1},
                   f::Function,df::Function,
                   s::Float64,b::Float64,mu::Float64,maxit::Int64)
    a  = s
    it = 0

    f_0  = f(x)
    df_0 = df(x)
    dftp = dot(df_0,p)

    x_1 = x + a*p
    f_1 = f(x_1)

    armres = f_1 - f_0 - mu*a*dftp
    while armres > 0 && maxit > it
        it += 1
        a  *= b

        x_1 = x + a*p
        f_1 = f(x_1)

        armres = f_1 - f_0 - mu*a*dftp
    end
    return x_1, a
end

function armijo_bt2(x::Array{Float64,1},p::Array{Float64,1},
                   f::Function,df::Function,
                   s::Float64,b::Float64,mu::Float64,maxit::Int64)
    a  = s
    it = 0

    f_0  = f(x)
    df_0 = df(x)
    dftp = dot(df_0,p)

    x_1 = x + a*p
    f_1 = f(x_1)

    armres = f_1 - f_0 - mu*a*dftp
    while armres <= 0 && maxit > it
        it += 1
        a  *= 1/b

        x_1 = x + a*p
        f_1 = f(x_1)

        armres = f_1 - f_0 - mu*a*dftp
        if armres > 0
            a *= b
            x_1 = x + a*p
            return x_1, a
        end
    end

    while armres > 0 && maxit > it
        it += 1
        a  *= b

        x_1 = x + a*p
        f_1 = f(x_1)

        armres = f_1 - f_0 - mu*a*dftp
    end
    return x_1, a
end

################################################################################
# Steepest descent with Armijo (allowing step size increases)
################################################################################
function steepest_descent(x::Array{Float64,1},
                          f::Function,df::Function,
                          atol::Float64,rtol::Float64,mxto::Int64,
                          s::Float64,b::Float64,mu::Float64,mxti::Int64)
    res = norm(df(x))
    r_0 = copy(res)
    it  = 0

    while res > rtol*r_0 + atol && it < mxto
        it += 1
        p   = -df(x)
        # x,s = armijo_bt(x,p,f,df,s,b,mu,mxto)
        x,s = armijo_bt2(x,p,f,df,s,b,mu,mxto)
        res = norm(df(x))

        println("f(x) = ", f(x))
        println("||∇f(x)|| = ", res)
    end

    return x
end

################################################################################
# Newton's method with line search
################################################################################
function  newton_ls(x::Array{Float64,1},
                    f::Function,df::Function,d2f::Function,
                    atol::Float64,rtol::Float64,mxto::Int64,
                    s::Float64,b::Float64,mu::Float64,mxti::Int64)

    res = norm(df(x))
    r_0 = copy(res)
    it  = 0

    while res > rtol*r_0 + atol && it < mxto
        it += 1
        H   = d2f(x)
        p   = -H\df(x)
        # x,s = armijo_bt(x,p,f,df,s,b,mu,mxto)
        x,s = armijo_bt2(x,p,f,df,s,b,mu,mxto)
        res = norm(df(x))

        println("f(x) = ", f(x))
        println("||∇f(x)|| = ", res)
    end

    return x
end

################################################################################
# A Trust Region Solver
################################################################################
function quad_model(p::Array{Float64,1},g::Array{Float64,1},B::Array{Float64,2})
    phi_p  = dot(g,p)
    Bp     = B*p
    phi_p += 0.5*dot(p,Bp)
    return phi_p
end

mutable struct QuadMod
    m_k::Function
      g::Array{Float64,1}
      B::Array{Float64,2}
end

function (q::QuadMod)(p::Array{Float64,1})
    return q.m_k(p,q.g,q.B)
end

# A simple bisection method to find tau in steihaug_cg
function root_finder(p::Array{Float64,1},d::Array{Float64,1},D::Float64)
    func(tau::Float64) = norm(p + tau*d) - D
    a = 0.0
    b = 0.01
    while func(b) < 0
        b *= 2
    end

    if func(b) == 0
        # print("1")
        return b
    end

    it = 0
    while it < 2000
        c = (a + b)/2
        if abs(func(c)) < 1e-10 || (b-a)/2 < 1e-10
            # print("2", c)
            return c
        end
        if sign(func(c)) == sign(func(a))
            a = c
        else
            b = c
        end
        it += 1
    end
    # print("3")
    return 0.0
end

# p = rand(2)
# d = rand(2)
# D = norm(p+16*d)
# root_finder(p,d,D)

# The Steihaug CG method for the Trust Region Subproblem
function steihaug_cg(m_k::QuadMod,Delta::Float64,tol::Float64,maxit::Int64)
    n = length(m_k.g)

    p_i = zeros(n)

    r_i  = -copy(m_k.g)
    rt_i =  copy(r_i) # There is an option to use a preconditioner here.
    d_i  =  copy(rt_i)

    Bd  = (m_k.B)*d_i
    g_i = dot(d_i,Bd)

    for i in 1:maxit
        if g_i <= 0
            tau = root_finder(p_i,d_i,Delta)
            # println("Here", i)
            return p_i + tau*d_i
        else # g_i > 0
            c_i = dot(r_i,rt_i)
            a_i = c_i/g_i

            p_i = p_i + a_i*d_i

            if norm(p_i) >= Delta
                p_i = p_i - a_i*d_i
                tau = root_finder(p_i,d_i,Delta)
                # println("Here", i)
                return p_i + tau*d_i
            end

            r_i = r_i - a_i*Bd

            if norm(r_i)/norm(m_k.g) <= tol
                # println("Here", i)
                return p_i
            end

            rt_i  = copy(r_i) # There is an option to use a preconditioner here.

            b_i = dot(r_i,rt_i)/c_i
            d_i = rt_i + b_i*d_i

            Bd  = (m_k.B)*d_i
            g_i = dot(d_i,Bd)
        end
    end
    return p_i
end

function trust_region_steihaug(x::Array{Float64,1},atol::Float64,mxto::Int64,
                               mu::Float64,lam::Float64,
                               f::Function,df::Function,d2f::Function,
                               D::Float64,Dmax::Float64,cg_tol::Float64,
                               mxti::Int64,m_k::QuadMod)
    m_k.g = df(x)
    m_k.B = d2f(x)

    res  = norm(m_k.g)
    res0 = copy(res)
    it   = 0

    while res > atol && it < mxto
        p = steihaug_cg(m_k,D,cg_tol,mxti)

        f_old = f(x)
        f_new = f(x + p)
        rho_k = (f_old - f_new)/ -m_k(p)
        if rho_k > lam
            x = x + p

            D*=2
            D = min(D,Dmax)
        else
            while rho_k <= lam
                p     = 0.5*p
                f_new = f(x + p)
                rho_k = (f_old - f_new)/ -m_k(p)
                println("ρ = ", rho_k)
            end
            x = x + p
            D*= 0.5
        end

        m_k.g = df(x)
        m_k.B = d2f(x)
        res = norm(df(x))
        it += 1
        println("f(x) = ", f(x))
    end

    return x
end

################################################################################
# Tests for steepest_descent
################################################################################
test_x = 5*rand(2)
steepest_descent(test_x,rosenbrock,d_rb,1e-4,1e-4,20000,1.0,0.5,0.01,20)
steepest_descent(test_x,beale,d_ble,1e-4,1e-4,5000,1.0,0.5,0.01,20)
steepest_descent(test_x,rastrigin,d_ras,1e-4,1e-4,5000,1.0,0.5,0.01,20)
steepest_descent(test_x,goldstein_price,d_gp,1e-8,1e-8,5000,1.0,0.5,0.01,20)
steepest_descent(test_x,booth,d_bth,1e-4,1e-4,5000,1.0,0.5,0.01,20)

steepest_descent(rand(50),rosenbrock,d_rb,1e-4,1e-4,20000,1.0,0.5,0.01,20)

################################################################################
# Tests for Newton LS
################################################################################
test_x = 5*rand(2)
newton_ls(test_x,rosenbrock,d_rb,d2_rb,1e-8,1e-8,10,1.0,0.5,0.01,20)
newton_ls(test_x,beale,d_ble,d2_ble,1e-4,1e-4,200,1.0,0.5,0.01,20)
newton_ls(test_x,rastrigin,d_ras,d2_ras,1e-4,1e-4,200,1.0,0.5,0.01,20)
newton_ls(test_x,goldstein_price,d_gp,d2_gp,1e-8,1e-8,5000,1.0,0.5,0.01,20)
newton_ls(test_x,booth,d_bth,d2_bth,1e-4,1e-4,5000,1.0,0.5,0.01,20)

################################################################################
# Tests for Trust Region
################################################################################
qm = QuadMod(quad_model,rand(2),rand(2,2))
test_x = rand(2)
trust_region_steihaug(test_x,1e-8,20,0.25,0.75,rosenbrock,d_rb,d2_rb,
                      10.0,1000.0,1e-8,100,qm)
trust_region_steihaug(test_x,1e-8,50,0.25,0.75,beale,d_ble,d2_ble,
                      10.0,10000.0,1e-8,100,qm)
trust_region_steihaug(test_x,1e-4,30,0.25,0.75,rastrigin,d_ras,d2_ras,
                      10.0,1000.0,1e-8,100,qm)
trust_region_steihaug(test_x,1e-7,30,0.25,0.75,goldstein_price,d_gp,d2_gp,
                      10.0,1000.0,1e-8,100,qm)
trust_region_steihaug(test_x,1e-7,30,0.25,0.75,booth,d_bth,d2_bth,
                      10.0,1000.0,1e-8,100,qm)

trust_region_steihaug(rand(50),1e-4,200,0.25,0.75,rosenbrock,d_rb,d2_rb,
                      10.0,1000.0,1e-8,100,qm)
################################################################################
# Test Functions
################################################################################
# The Rosenbrock function. Global Minimizers at (1,1,1...,1)
function rosenbrock(x::Array{Float64,1})
    f = 0.0
    n = length(x)
    for i in 1:(n-1)
        f += 100*(x[i+1]-x[i]^2)^2 + (1-x[i])^2
    end
    return f
end

# Surface plot of Rosenbrock in 2D
x_grid = range(-2,stop=2,length=100)
y_grid = range(-1,stop=3,length=100)
rb_2D(x,y) = rosenbrock([x,y])

my_grad = cgrad(:roma, 1000, categorical = true, scale = :exp)
plot(x_grid,y_grid,rb_2D,st=:surface,c=my_grad,camera=(-30,30))
contour(x_grid,y_grid,rb_2D,levels=1000)

# Approximate gradient of Rosenbrock
d_rb(x::Array{Float64,1}) = CD_grad(x,1e-6,rosenbrock)

# Approximate Hessian of Rosenbrock
d2_rb(x::Array{Float64,1}) = CD_hess(x,1e-6,rosenbrock)
d2_rb(rand(2))
################################################################################
# The Beale function. Global Minimizer (3,0.5)
function beale(x::Array{Float64,1})
    return (1.5 - x[1] + x[1]*x[2])^2 + (2.25 - x[1] + x[1]*x[2]^2)^2 +
            (2.625 - x[1] + x[1]*x[2]^3)^2
end

# Surface plot of Beale
x_grid=range(-4,stop=4,length=100)
y_grid=range(-4,stop=4,length=100)
beale_2D(x,y) = (1.5 - x + x*y)^2 + (2.25 - x + x*y^2)^2 + (2.625 - x + x*y^3)^2

my_grad = cgrad(:roma, 1000, categorical = true, scale = :exp)
plot(x_grid,y_grid,beale_2D,st=:surface,c=my_grad,camera=(-30,30))
contour(x_grid,y_grid,beale_2D,levels=500)

# Approximate gradient of Beale
d_ble(x::Array{Float64,1}) = CD_grad(x,1e-6,beale)
d_ble(rand(2))

# Approximate Hessian of Beale
d2_ble(x::Array{Float64,1}) = CD_hess(x,1e-6,beale)
d2_ble(rand(2))
################################################################################

################################################################################
# The Rastrigin function. Global Minimizers at (0,0,0...0)
function rastrigin(x::Array{Float64,1})
    f = 0.0
    n = length(x)
    for i in 1:n
        f += x[i]^2 - 10*cos(2*pi*x[i])
    end
    return 10*n + f
end

# Surface plot of Rastrigin
x_grid=range(-5,stop=5,length=100)
y_grid=range(-5,stop=5,length=100)
ras_2D(x,y) = 10*2 + x^2 - 10*cos(2*pi*x) + y^2 - 10*cos(2*pi*y)

my_grad = cgrad(:roma, 1000, categorical = true, scale = :exp)
plot(x_grid,y_grid,ras_2D,st=:surface,c=my_grad,camera=(-30,15))
contour(x_grid,y_grid,ras_2D,levels=1000)

# Approximate gradient of Rastrigin
d_ras(x::Array{Float64,1}) = CD_grad(x,1e-6,rastrigin)
d_ras(rand(10))

# Approximate Hessian of Rastrigin
d2_ras(x::Array{Float64,1}) = CD_hess(x,1e-6,rastrigin)
println(d2_ras(rand(2)))
################################################################################

################################################################################
# The Goldstein-Price function. Global Minimizers at (0,-1)
function goldstein_price_2D(x::Float64,y::Float64)
    a = (1 + (x + y + 1)^2 * (19 - 14*x + 3*x^2 - 14*y + 6*x*y + 3*y^2))
    b = (30 + (2*x - 3*y)^2 *(18 - 32*x + 12*x^2 + 48*y - 36*x*y + 27*y^2))
    return a * b
end

function goldstein_price(x::Array{Float64,1})
    a = (1 + (x[1] + x[2] + 1)^2 * (19 - 14*x[1] + 3*x[1]^2 - 14*x[2] +
            6*x[1]*x[2] + 3*x[2]^2))
    b = (30 + (2*x[1] - 3*x[2])^2 *(18 - 32*x[1] + 12*x[1]^2 + 48*x[2] -
            36*x[1]*x[2] + 27*x[2]^2))
    return a * b
end

# Surface plot of Goldstein-Price
x_grid=range(-2,stop=2,length=100)
y_grid=range(-2,stop=2,length=100)

my_grad = cgrad(:roma, 1000, categorical = true, scale = :exp)
plot(x_grid,y_grid,goldstein_price_2D,st=:surface,c=my_grad,camera=(-30,30))
contour(x_grid,y_grid,goldstein_price_2D,levels=800)

# Approximate gradient of Goldstein-Price
d_gp(x::Array{Float64,1}) = CD_grad(x,1e-6,goldstein_price)
d_gp(rand(2))

# Approximate Hessian of Goldstein-Price
d2_gp(x::Array{Float64,1}) = CD_hess(x,1e-6,goldstein_price)
d2_gp(rand(2))
################################################################################
# The Booth function. Global Minimizers at (1,3)
function booth(x::Array{Float64,1})
    return (x[1] + 2*x[2] - 7)^2 + (2*x[1] + x[2] - 5)^2
end

# Surface plot of Booth
x_grid=range(-5,stop=5,length=100)
y_grid=range(-5,stop=5,length=100)

my_grad = cgrad(:roma, 1000, categorical = true, scale = :exp)
booth2d(x,y) = booth([x,y])
plot(x_grid,y_grid,booth2d,st=:surface,c=my_grad,camera=(-30,30))
contour(x_grid,y_grid,booth2d,levels=800)

# Approximate gradient of Booth
d_bth(x::Array{Float64,1}) = CD_grad(x,1e-6,booth)
d_bth(rand(2))

# Approximate Hessian of Booth
d2_bth(x::Array{Float64,1}) = CD_hess(x,1e-6,booth)
d2_bth(rand(2))

################################################################################
################################################################################
################################################################################
# Approximating Gradients and Hessians using Finite Differences
################################################################################
################################################################################
################################################################################
################################################################################
# A central difference function for finding partial derivatives
function CD(x::Array{Float64},d::Array{Float64},h::Float64,f::Function)
    return (f(x+d) - f(x-d))/(2*h)
end

# A function for approximating the gradient using finite differences
function CD_grad(x::Array{Float64},h::Float64,f::Function)
    n  = length(x)
    it = 1

    gradF = zeros(n)
    while it <= n
        gradF[it] = CD(x,h*(1:n .== it),h,f)
        it += 1
    end

    return gradF
end

function CD_hess(x::Array{Float64,1},h::Float64,f::Function)
     n = length(x)
     i = 1
     j = 1

    HessF = zeros(n,n)
    while i <= n
        while j <= n
            HessF[i,j] = (CD_grad(x+h*(1:n .== j),h,f)[i] -
                          CD_grad(x-h*(1:n .== j),h,f)[i] +
                          CD_grad(x+h*(1:n .== i),h,f)[j] -
                          CD_grad(x-h*(1:n .== i),h,f)[j])/(4.0*h)
            j+=1
        end
        j  = 1
        i += 1
    end

    return 0.5*(HessF + HessF')
end


text:
export GPC

abstract type AbstractGP end

# A collection of GPs (GPC == "GP Collection"). Used to keep track of GPs.
mutable struct GPC
    n::Int
    GPC() = new(0)
end

@nograd GPC

next_index(gpc::GPC) = gpc.n + 1

#
# Projecting AbstractGPs onto a finite dimensional marginal
#

import Base: rand, length
import Distributions: logpdf, ContinuousMultivariateDistribution

export mean, std, cov, marginals, rand, logpdf, elbo

"""
    FiniteGP{Tf<:AbstractGP, Tx<:AV, TΣy}

The finite-dimensional projection of the AbstractGP `f` at `x`.
"""
struct FiniteGP{Tf<:AbstractGP, Tx<:AV, TΣy} <: ContinuousMultivariateDistribution
    f::Tf
    x::Tx
    Σy::TΣy
end
FiniteGP(f::AbstractGP, x::AV, σ²::AV{<:Real}) = FiniteGP(f, x, Diagonal(σ²))
FiniteGP(f::AbstractGP, x::AV, σ²::Real) = FiniteGP(f, x, Fill(σ², length(x)))
FiniteGP(f::AbstractGP, x::AV) = FiniteGP(f, x, 1e-18)

length(f::FiniteGP) = length(f.x)

"""
    mean(f::FiniteGP)

The mean vector of `f`.
"""
mean(f::FiniteGP) = mean_vector(f.f, f.x)

"""
    cov(f::FiniteGP)

The covariance matrix of `f`.
"""
cov(f::FiniteGP) = cov(f.f, f.x) + f.Σy

"""
    cov(f::FiniteGP, g::FiniteGP)

The cross-covariance between `f` and `g`.
"""
cov(f::FiniteGP, g::FiniteGP) = cov(f.f, g.f, f.x, g.x)

"""
    marginals(f::FiniteGP)

Sugar, returns a vector of Normal distributions representing the marginals of `f`.
"""
marginals(f::FiniteGP) = Normal.(mean(f), sqrt.(cov_diag(f.f, f.x) .+ diag(f.Σy)))

"""
    rand(rng::AbstractRNG, f::FiniteGP, N::Int=1)

Obtain `N` independent samples from the marginals `f` using `rng`.
"""
function rand(rng::AbstractRNG, f::FiniteGP, N::Int)
    μ, C = mean(f), cholesky(Symmetric(cov(f)))
    return μ .+ C.U' * randn(rng, length(μ), N)
end
rand(f::FiniteGP, N::Int) = rand(Random.GLOBAL_RNG, f, N)
rand(rng::AbstractRNG, f::FiniteGP) = vec(rand(rng, f, 1))
rand(f::FiniteGP) = vec(rand(f, 1))

"""
    logpdf(f::FiniteGP, Y::AbstractMatrix{<:Real})

The log probability density of each column of `Y` under `f`.
"""
function logpdf(f::FiniteGP, y::AbstractVector{<:Real})
    μ, C = mean(f), cholesky(Symmetric(cov(f)))
    return -(length(y) * log(2π) + logdet(C) + Xt_invA_X(C, y - μ)) / 2
end

"""
    logpdf(f::FiniteGP, y::AbstractVector{<:Real})

The log probability density of `y` under `f`.
"""
function logpdf(f::FiniteGP, Y::AbstractMatrix{<:Real})
    μ, C = mean(f), cholesky(Symmetric(cov(f)))
    return -((size(Y, 1) * log(2π) + logdet(C)) .+ diag_Xt_invA_X(C, Y .- μ)) ./ 2
end

"""
   elbo(f::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)

The saturated Titsias-ELBO.
"""
function elbo(f::FiniteGP, y::AV{<:Real}, u::FiniteGP)
    consistency_check(f, y, u)
    chol_Σy = cholesky(f.Σy)

    A = cholesky(Symmetric(cov(u))).U' \ (chol_Σy.U' \ cov(f, u))'
    Λ_ε = cholesky(Symmetric(A * A' + I))
    δ = chol_Σy.U' \ (y - mean(f))

    return -(length(y) * log(2π) + logdet(chol_Σy) + logdet(Λ_ε) +
        sum(abs2, δ) - sum(abs2, Λ_ε.U' \ (A * δ)) +
        tr_Cf_invΣy(f, f.Σy, chol_Σy) - sum(abs2, A)) / 2
end

function consistency_check(f, y, u)
    @assert length(f) == size(y, 1)
end
Zygote.@nograd consistency_check

import Base: \
\(A::AbstractMatrix, B::Diagonal) = A \ Matrix(B)

\(A::Union{LowerTriangular, UpperTriangular}, B::Diagonal) = A \ Matrix(B)
\(A::Adjoint{<:Any, <:Union{LowerTriangular, UpperTriangular}}, B::Diagonal) = A \ Matrix(B)

# Compute tr(Cf / Σy) efficiently for different types of Σy. For dense Σy you obviously need
# to compute the entirety of Cf, which is bad, but for particular structured Σy one requires
# only a subset of the elements. Σy isa UniformScaling is version usually considered.
function tr_Cf_invΣy(f::FiniteGP, Σy::UniformScaling, chol_Σy::Cholesky)
    return sum(cov_diag(f.f, f.x)) / Σy.λ
end
function tr_Cf_invΣy(f::FiniteGP, Σy::Diagonal, chol_Σy::Cholesky)
    return sum(cov_diag(f.f, f.x) ./ diag(Σy))
end
tr_Cf_invΣy(f::FiniteGP, Σy::Matrix, chol_Σy::Cholesky) = tr(chol_Σy \ cov(f))
# function tr_Cf_invΣy(f::FiniteGP, Σy::BlockDiagonal, chol_Σy::Cholesky)
#     C = cholesky(Symmetric(_get_kernel_block_diag(f, cumulsizes(Σy, 1))))
#     return tr_At_A(chol_Σy.U' \ C.U')
# end

# function _get_kernel_block_diag(f::FiniteGP, cs)
#     k = kernel(f.f)
#     ids = map(n->cs[n]:cs[n+1]-1, 1:length(cs)-1)
#     xs = map(id->f.x[id], ids)
#     Σs = map(x->pw(k, x), xs)
#     return block_diagonal(Σs)
# end

# function _get_kernel_block_diag(f::FiniteGP{<:GP{<:BlockMean, <:BlockKernel}, <:BlockData}, cs)
#     k = kernel(f.f)
#     ids = map(n->cs[n]:cs[n+1]-1, 1:length(cs)-1)
#     @assert _test_block_consistency(ids, f)
#     xs = blocks(f.x)
#     Σs = map(n->pw(k.ks[n], xs[n]), 1:length(xs))
#     return block_diagonal(Σs)
# end

_test_block_consistency(ids, f) = length.(ids) == length.(blocks(f.x))
Zygote.@nograd _test_block_consistency

# """
#     elbo(f::FiniteGP, y::AV{<:Real}, u::FiniteGP, mε::AV{<:Real}, Λε::AM{<:Real})

# The unsaturated Titsias-ELBO.
# """
# function elbo(f::FiniteGP, y::AV{<:Real}, u::FiniteGP, mε::AV{<:Real}, Λε::AM{<:Real})
#     @assert length(u.x) == length(mε)
#     @assert size(Λε) == (length(mε), length(mε))
#     # do stuff.
# end

import Base: |, merge
export ←, |, Obs

"""
    Observation

Represents fixing a paricular (finite) GP to have a particular (vector) value.
"""
struct Observation{Tf<:FiniteGP, Ty<:Vector}
    f::Tf
    y::Ty
end

const Obs = Observation

←(f, y) = Observation(f, y)
get_f(c::Observation) = c.f
get_y(c::Observation) = c.y


text:
using SnoopCompile

### Log the compiles
# This only needs to be run once (to generate "/tmp/images_compiles.csv")

SnoopCompile.@snoop "/tmp/images_compiles.csv" begin
    include(Pkg.dir("Images", "test","runtests.jl"))
end

### Parse the compiles and generate precompilation scripts
# This can be run repeatedly

# IMPORTANT: we must have the module(s) defined for the parcelation
# step, otherwise we will get no precompiles for the Images module
using Images

data = SnoopCompile.read("/tmp/images_compiles.csv")

# The Images tests are run inside a module ImagesTest, so all
# the precompiles get credited to ImagesTest. Credit them to Images instead.
subst = Dict("ImagesTests"=>"Images")

# Blacklist helps fix problems:
# - MIME uses type-parameters with symbols like :image/png, which is
#   not parseable
blacklist = ["MIME"]

# Use these two lines if you want to create precompile functions for
# individual packages
pc, discards = SnoopCompile.parcel(data[end:-1:1,2], subst=subst, blacklist=blacklist)
SnoopCompile.write("/tmp/precompile", pc)

# Use these two lines if you want to add to your userimg.jl
pc = SnoopCompile.format_userimg(data[end:-1:1,2], subst=subst, blacklist=blacklist)
SnoopCompile.write("/tmp/userimg_Images.jl", pc)


text:
using MyExample
using Test

@testset "MyExample.jl" begin

    # 2x + 3y
  @test my_f(2,1) == 7
  @test my_f(2,3) == 13
  @test my_f(1,3) == 11

end

@testset "Derivative Tests" begin

  # dx(2x +3y) = 2
  @test derivative_of_my_f(2,1) == 2
  @test derivative_of_my_f(20,1) == 2
  @test derivative_of_my_f(20,10) == 2

end


text:
module matrix_io
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)

  module common
    const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
    include("entity_pb.jl")
  end
  module malos
    const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
    module v1
      const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
      module comm
        const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
        include("zwave_commands_pb.jl")
        include("comm_pb.jl")
      end
      module io
        const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
        include("io_pb.jl")
      end
      module maloseye
        const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
        include("maloseye_pb.jl")
      end
      module sense
        const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
        include("sense_pb.jl")
      end
      module driver
        const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
        include("driver_pb.jl")
      end
    end
  end
end


text:
module NMRData

# Write your package code here.
include("./get_path.jl")

export getBMRBexperimentpaths, getBMRBexperimentpath

end


text:
using StratiGraphics
using Meshes
using GeoStatsBase
using ReferenceTests, ImageIO
using Test, Plots, Random
gr(size=(600,400))

# workaround for GR warnings
ENV["GKSwstype"] = "100"

# environment settings
isCI = "CI" ∈ keys(ENV)
islinux = Sys.islinux()
visualtests = !isCI || (isCI && islinux)
datadir = joinpath(@__DIR__,"data")

include("dummysolver.jl")

@testset "StratiGraphics.jl" begin
  if visualtests
    rng    = MersenneTwister(2019)
    proc   = GeoStatsProcess(Dummy())
    env    = Environment(rng, [proc, proc], [0.5 0.5; 0.5 0.5], ExponentialDuration(rng, 1.0))
    record = simulate(env, LandState(zeros(50,50)), 10)
    strata = Strata(record)

    @test_reference "data/strata.png" plot(strata)

    for (i, fillbase, filltop) in [(1, NaN, NaN), (2, 0, NaN), (3, 0, 0)]
      rng  = MersenneTwister(2019)
      proc = GeoStatsProcess(Dummy())
      env  = Environment(rng, [proc, proc], [0.5 0.5; 0.5 0.5], ExponentialDuration(rng, 1.0))
      prob = SimulationProblem(CartesianGrid(50,50,20), :strata => Float64, 3)
      solv = StratSim(:strata => (environment=env, fillbase=fillbase, filltop=filltop))
      sol  = solve(prob, solv)
      plts = map(sol) do real
        R = asarray(real, :strata)
        heatmap(rotr90(R[1,:,:]))
      end
      @test_reference "data/voxel$i.png" plot(plts..., layout=(3,1))
    end
  end
end


text:
function draw_system(bd::CelestialObject, starttime, endtime=nothing; bodies=true, sois=true, objects=false, refline=false, kwargs...)
    traces = AbstractTrace[]
    push!(traces, draw_central_body(bd))
    if hasfield(typeof(bd), :orbit)
        push!(traces, draw_central_body_orbit(bd, starttime))
    end
    if bodies
        for body in bd.satellite_bodies
            push!(traces, draw_orbiting_body(body, starttime))
            append!(traces, draw_orbit(body, starttime, endtime; kwargs...))
        end
    end
    if sois
        for body in bd.satellite_bodies
            push!(traces, draw_soi(body, starttime))
        end
    end
    if objects
        for ob in bd.satellite_objects
            append!(traces, draw_orbit(ob, starttime, endtime; kwargs...))
        end
    end
    if refline
        push!(traces, draw_ref_line(orbit_plot_size(bd)))
    end
    return traces
end


text:
using WAV, Images, VideoIO, Colors
using Base.Iterators, Random, StatsBase

# attempt at non trash code lol
# simple objective: random isomorphism between wav file and an mp4 video
# given a wav file creates a music video that preserves the structure of the bits

function frames(fname; fps=60) #::Union{Nothing, Any}
	s, fs = wavread(fname)
	frames(s, fs, fps)
end

frames(s, fs, fps) = partition(eachrow(s), Int(fs / fps))

"dict from samples to grey shades of type RGB (idkwhylol)"
function gmap(s)
	us = unique(s)
	Dict(us .=> range(RGB(0), RGB(1), length=length(us)))
end

"dict colormap, but written dumb, how do i sample from multiple ranges"
function cmap(s; )
	us = unique(s)
	n = length(us)
	Dict(us .=> sample(all_colors_fast(), n, replace=false))
end

"horrible function for getting all colors, doesn't use rand. dont use this"
function all_colors()
	arr = RGB[]
	for i in 1:255
		for j in 1:255
			for k in 1:255
				f = [i, j, k] ./ 255
				push!(arr, RGB(f...))
			end
		end
	end
	arr
end

# thanks kimikage
mutable struct Xorshift24 <: AbstractRNG
    prev::UInt32
    Xorshift24(seed::Integer) = new(seed % 0xFFFFFF + 1)
end

function Random.rand(rng::Xorshift24, ::Random.SamplerType{UInt32}) # The parameters are just an example.  I don't know the randomness.
    x = rng.prev
    x ⊻= x >> 13
    x ⊻= (x << 19) >> 8
    x ⊻= x >> 2
    rng.prev = x
end

"better version of all_colors"
all_colors_fast() = reinterpret(RGB24, rand(Xorshift24(0), UInt32, 2^24 - 1))

function imgs_from_signal(s; fs=44_100, fps=60)
	d = cmap(s)
	map(x->reshape(hcat(x...), 30, 49), collect(frames(map(x->d[x], s), fs, fps))[1:end-1])
end

# todo: use AxisArrays for the times
# write sound and images into mp4
# save imgs, seems like FileIO was being a little jank.

# might use Plots mp4() to write then just add the sound on top? idk

# TEST_PATH="data/072520.wav"

function main(wavpath)
	rootname = split(split(wavpath, "/")[end], ".")[1]
	outname = "data/imgs"
	s, fs = wavread(wavpath)
	imgs = imgs_from_signal(s)
	if isdir(outname)
		error("folder exists brooo")
	else
		mkpath(outname)
		map(x->save("$outname/$(x[1]).png", x[2]), enumerate(imgs))
	end
end

main(ARGS[1])

max_stars_repo_name: anandijain/audio

max_issues_repo_name: anandijain/audio

max_forks_repo_name: anandijain/audio


text:
module SplitApplyCombine

using Base: @propagate_inbounds, @pure, promote_op, Callable
using Indexing
using Dictionaries
import Dictionaries: filterview

import Base: merge, merge!, size, IndexStyle, getindex, parent, axes, iterate

# collections -> scalar
if VERSION < v"1.4.0-DEV"
    export only
    include("only.jl")
end

# collections -> collections
export mapmany, mapview, MappedIterator, MappedArray, product, productview, ProductArray, filterview

# collections -> collections of collections
export group, groupfind, GroupDictionary, groupview, groupreduce, groupcount, groupsum, groupprod, groupunique, grouponly, groupfirst, grouplast
export splitdims, splitdimsview, SplitDimsArray

# colletions of collections -> collections
export flatten #, flattenview
export combinedims, combinedimsview, CombineDimsArray

# collections of collections -> collections of collections
export invert, invert! # a "transpose" for nested containers
export innerjoin, ⨝, leftgroupjoin

@inline former(a, b) = a
@inline latter(a, b) = b

include("map.jl")
include("group.jl")
include("product.jl")
include("innerjoin.jl")
include("leftgroupjoin.jl")
include("splitdims.jl")
include("combinedims.jl")
include("invert.jl")

# Silly definitions missing from Base
# (defining them here is piracy)
# ===================================

if !hasmethod(Base.haskey, (Any, Any))
    Base.haskey(a, i) = i ∈ keys(a)
end

if VERSION < v"1.2"
    keytype(a::AbstractArray) = eltype(keys(a))
    keytype(a) = Base.keytype(a)
end

end # module


text:
using Documenter, SparseUtils

makedocs(
    debug = true,
    strict = false,
    doctest = false,
    format = :html,
    sitename = "SparseUtils.jl",
    modules = [SparseUtils],
    pages = [
        "index.md"
    ]
)

deploydocs(
    repo = "github.com/jlapeyre/SparseUtils.jl.git",
    target = "build",
    julia  = "nightly",
    deps = nothing,
    make = nothing
)


text:
# Periodic Function

"""
    Periodic <: Isotropic{Euclidean}

Periodic kernel (covariance)
```math
k(x,x') = σ²\\exp(-2\\sin²(π|x-x'|/p)/ℓ²)
```
with length scale ``ℓ``, signal standard deviation ``σ``, and period ``p``.
"""
mutable struct Periodic <: Isotropic{Euclidean}
    "Squared length scale"
    ℓ2::Float64
    "Signal variance"
    σ2::Float64
    "Period"
    p::Float64
    "Priors for kernel parameters"
    priors::Array

    """
        Periodic(ll::Float64, lσ::Float64, lp::Float64)

    Create `Periodic` with length scale `exp(ll)`, signal standard deviation `exp(lσ)`, and
    period `exp(lp)`.
    """
    Periodic(ll::Float64, lσ::Float64, lp::Float64) =
        new(exp(2 * ll), exp(2 * lσ), exp(lp), [])
end

get_params(pe::Periodic) = Float64[log(pe.ℓ2) / 2, log(pe.σ2) / 2, log(pe.p)]
get_param_names(pe::Periodic) = [:ll, :lσ, :lp]
num_params(pe::Periodic) = 3

function set_params!(pe::Periodic, hyp::VecF64)
    length(hyp) == 3 || throw(ArgumentError("Periodic function has only three parameters"))
    pe.ℓ2, pe.σ2 = exp(2 * hyp[1]), exp(2 * hyp[2])
    pe.p = exp(hyp[3])
end

Statistics.cov(pe::Periodic, r::Float64) = pe.σ2 * exp(-2 / pe.ℓ2 * sin(π * r / pe.p)^2)

@inline dk_dll(pe::Periodic, r::Float64) =
    (s = 2 * sin(π * r / pe.p)^2 / pe.ℓ2; 2 * pe.σ2 * s * exp(-s)) # dK_dlogℓ
@inline dk_dlp(pe::Periodic, r::Float64) =
    (s = π * r / pe.p; t = 2 / pe.ℓ2; pe.σ2 * s * t * sin(2 * s) * exp(-t * sin(s)^2))  # dK_dlogp

@inline function dk_dθp(pe::Periodic, r::Float64, p::Int)
    if p==1
        dk_dll(pe, r)
    elseif p==2
        dk_dlσ(pe, r)
    elseif p==3
        dk_dlp(pe, r)
    else
        return NaN
    end
end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# 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.
##########################################################################

# using histogram from StatsBase

function Base.:+(h1::Histogram, h2::Histogram)::Histogram
	if h2 == nullHist return deepcopy(h1) end
	if h1 == nullHist return deepcopy(h2) end

	if (h1.closed == h2.closed) && (h1.isdensity == h2.isdensity) && (length(h1.edges) == 1 && length(h2.edges) == 1)
		e1, e2 = h1.edges[1], h2.edges[1] # shorthand
		if length(e1) < length(e2)
			h1, h2, e1, e2 = h2, h1, e2, e1 # swap h1 and h2, so that h1 is the 'longer' histogram
		end
		t1, t2 = typeof(e1), typeof(e2)
		# allow histograms to be added if they have the same bin size, with bins starting at the same point, but ending at different points
		if (t1 == t2) && (e1 != e2) && # note: if e1 == e2 here, then could just use merge(h1,h2)
			((t1 <: UnitRange && e1.start == e2.start) ||
			(t1 <: StepRange && e1.start == e2.start && e1.step == e2.step) ||
			(t1 <: StepRangeLen && e1.ref == e2.ref && e1.step == e2.step))

			h = deepcopy(h1)
			h.weights[1:length(h2.weights)] += h2.weights
			return h
		end
	end

	return merge(h1, h2) # requires that histograms have same field values, except for bin weights which just need to have equal length
end

function Base.:-(h::Histogram)::Histogram
	h2 = deepcopy(h)
	h2.weights = -h2.weights
	return h2
end

function Base.:-(h1::Histogram, h2::Histogram)::Histogram
	return h1 + (-h2)
end


text:
#list of added packages

#routines for adding the packages
add Atom Juno
add DataFrames Revise Distributions StatsBase GLM Dates NLopt ForwardDiff Formatting
add DataStructures LinearAlgebra Gadfly Measures Formatting CodecZlib CodecLz4
add CSV NLopt CategoricalArrays Serialization Distributed Random
add StaticArrays VegaLite CuArrays
add https://github.com/clintonTE/Finometrics

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# Mathematical & Numerical functions
############################################

export locate,
       interp,
       deriv,
       integ,
       cuminteg,
       smooth!,
       smooth,
       smooth_spline,
       smooth_plaw,
       gauss_laguerre_nw,
       gauss_legendre_nw,
       expi

############################################
# locate
############################################
# Binary search for vectors. Returns matching index.
# If value is outside of the array, return -1
############################################
function locate(arr::AbstractVector, #input vector
                x::Real) #input value

    n = length(arr)
    ascnd = (arr[n] >= arr[1])
    jl = 0
    ju = n+1
    while true
        if ju-jl <= 1
            break
        end
        #jm = int((ju+jl)/2)
        jm = round(Int,((ju+jl)/2))
        if ascnd == (x >= arr[jm])
            jl = jm
        else
            ju = jm
        end
    end

    if x == arr[1]
        return 1
    elseif x == arr[n]
        return n-1
    elseif ascnd && (x > arr[n] || x < arr[1])
        return -1
    elseif !ascnd && (x < arr[n] || x > arr[1])
        return -1
    else
        return jl
    end
end

############################################
# interp
############################################
# Interpolation function for tables
#
# Cubic interpolation & linear extrapolation (lerp)
############################################
function interp(xold::AbstractVector, #grid
                fold::AbstractVector, #values
                xnew::Real; #new value
                method=:cubic) #:cubic, :lin

    n=length(xold)

    #Determining if interpolation or extrapolation
    xmin=xold[1]
    xmax=xold[end]

    #Interpolation
    if xnew >= xmin && xnew <= xmax

        #locate the current index
        ind = locate(xold, xnew)
        ind = ind < 3 ? 3 : ind+1

        #TODO: method keyword to toggle linear interpolation
        #Linear interpolation
        if method == :lin
            ind=ind-1
            ind2=ind
            fnew=fold[ind2]+(fold[ind2]-fold[ind])*(xnew-x2)/(x2-x1)
        elseif method == :cubic
            #Cubic interpolation
            x1=xold[ind-2]
            x2=xold[ind-1]
            x3=xold[ind]
            y1=fold[ind-2]
            y2=fold[ind-1]
            y3=fold[ind]

            fnew = ((x3-xnew)*((x2-x3)*(x2-xnew)*y1+(x3-x1)*(x1-xnew)*y2)+(x1-x2)*(x1-xnew)*(x2-xnew)*y3)/((x1-x2)*(x1-x3)*(x2-x3))
        else
            error("Check method")
        end

    #Extrapolation (inner)
    elseif xnew < xmin
        x1=xold[1]
        x2=xold[2]
        fnew = fold[1]+(fold[2]-fold[1])*(xnew-x1)/(x2-x1)

    #Extrapolation (outer)
    elseif xnew > xmax
        x1=xold[n-1]
        x2=xold[n]
        fnew = fold[n]+(fold[n]-fold[n-1])*(xnew-x2)/(x2-x1)

    end

    return fnew
end

#vectorized interp
function interp(xold::AbstractVector,
                fold::AbstractVector,
                xnew::AbstractVector;
                method=:cubic)

    N = length(xnew)
    fnewarr = zeros(N)
    for i = 1:N
        fnewarr[i] = interp(xold, fold, xnew[i], method=method)
    end

    return fnewarr
end

############################################
# Deriv
############################################
# First derivative of vector
#
# Uses weighted average of two parabolas to compute the derivative
############################################
function deriv(x::AbstractVector, #grid
               f::AbstractVector) #values

    n=length(x)
    @assert n > 1
    @assert length(x) == length(f)

    dfdx=zeros(n)

    dfdx[1]=(f[2]-f[1])/(x[2]-x[1])
    n1=n-1
    dfdx[n]=(f[n]-f[n1])/(x[n]-x[n1])

    if(n==2)
        return dfdx
    end

    s=abs(x[2]-x[1])/(x[2]-x[1])
    for j in 2:n1
        scale=maximum([abs(f[j-1]),abs(f[j]),abs(f[j+1])])/abs(x[j])
        if(scale==0.0)
            scale=1.0
        end

        d1=(f[j+1]-f[j])/(x[j+1]-x[j])/scale
        d=(f[j]-f[j-1])/(x[j]-x[j-1])/scale
        tann1=d1/(s*sqrt(1.0+d1^2.0)+1.0)
        tann=d/(s*sqrt(1.0+d^2.0)+1.0)
        dfdx[j]=(tann1+tann)/(1.0-tann1*tann)*scale
    end

   return dfdx
end

############################################
# Integ
############################################
# Integration of arrays
#
# Uses weighted parabolas
############################################

#Weighted parabolas
function integ(x::AbstractVector,
               f::AbstractVector)

    function parcoe(f, x)
        a = similar(f)
        b = similar(f)
        c = similar(f)

        n = length(f)

        c[1] = 0.0
        b[1] = (f[2]-f[1])/(x[2]-x[1])
        a[1] = f[1]-x[1]*b[1]

        c[n] = 0.0
        b[n] = (f[n]-f[n-1])/(x[n]-x[n-1])
        a[n] = f[n]-x[n]*b[n]
        if n == 2; return a, b, c; end

        for j = 2:(n-1)
            d = (f[j]-f[j-1])/(x[j]-x[j-1])
            c[j] = f[j+1]/((x[j+1]-x[j])*(x[j+1]-x[j-1]))-f[j]/((x[j]-x[j-1])*(x[j+1]-x[j]))+f[j-1]/((x[j]-x[j-1])*(x[j+1]-x[j-1]))
            b[j] = d- (x[j]+x[j-1])*c[j]
            a[j] = f[j-1]-x[j-1]*d+x[j]*x[j-1]*c[j]
        end

        for j = 2:(n-2)
            if c[j] == 0.0; continue; end
            wt = abs(c[j+1])/(abs(c[j+1])+abs(c[j]))
            a[j] = a[j+1]+wt*(a[j]-a[j+1])
            b[j] = b[j+1]+wt*(b[j]-b[j+1])
            c[j] = c[j+1]+wt*(c[j]-c[j+1])
        end

        return a, b, c
    end

    n = length(x)
    if n == 1
        return x .* f
    end

    a, b, c = parcoe(f, x)
    fint = zeros(n)

    fint[1] = 0.0
    #finti = (a[1]+(b[1]/2.0+c[1]/3.0*x[2])*x[2])*x[2]-x[1]
    finti = (a[1]+b[1]/2.0*(x[2]+x[1])+c[1]/3.0*((x[2]+x[1])*x[2]+x[1]*x[1]))*(x[2]-x[1])

    #if unstable, use linear interpolation
    #TODO: check when this is needed
    finti = finti < 0.0 ? 0.5*(x[2]-x[1])*(f[1]+f[2]) : finti
    fint[2] = finti

    if n == 2; return fint; end

    for i = 2:(n-1)
        finti = (a[i]+b[i]/2.0*(x[i+1]+x[i])+c[i]/3.0*((x[i+1]+x[i])*x[i+1]+x[i]*x[i]))*(x[i+1]-x[i])

        #if unstable, use linear interpolation
        #TODO: check when this is needed
        finti = finti < 0.0 ? 0.5*(x[i+1]-x[i])*(f[i]+f[i+1]) : finti

        fint[i+1] = finti
    end
    return fint
end

#Cumulative integral
function cuminteg(x::AbstractVector,
                  f::AbstractVector;
                  extrapolate=:zlin) #:zlin, lin, :zplaw, :none

    fint = integ(x, f)

    #extrapolate from 0 to x1
    if extrapolate == :zlin #linear [0,x1]
        fint[1] = 0.5*x[1]*f[1]
#    elseif extrapolate_zero == :quad #quadratic
#        c = 0.0
#        b = (f[2]-f[1])/(x[2]-x[1])
#        a = f[1]-x[1]*b[1]
#        fint[1] = (a+(b/2.0+(c/3.0*x[1])*x[1])*x[1])
    elseif extrapolate == :lin #linear
        k = (f[2]-f[1])/(x[2]-x[1])
        b = f[1] - k*x[1]
        if b >= 0.0
            fint[1] = 0.5*k*x[1]^2+b*x[1]
        else
            x0 = -b/k
            fint[1] = 0.5*k*(x[1]^2.0-x0^2.0)+b*(x[1]-x0)
        end
    elseif extrapolate == :zplaw #powerlaw [0,x1]
        a=log(f[2]/f[1])/log(x[2]/x[1])
        c=f[1]/(x[1]^a)
        fint[1]=c*(x[1]^(a+1))/(a+1)
    elseif extrapolate == :none #no extrapolation
        fint[1] = 0.0
    else
        error("Unrecognized extrapolation method $(string(extrapolate))")
    end

    cfint = cumsum(fint)
    return cfint
end

#Simpsons rule (cumulative)
#may or may not work
function integ_Simpson(x,f)

    n = length(f)
    fint = Array(Float64,n)

    #Linear extrapolation between 0-x1
    fint[1] = (x[1]/(2.0*(x[1]-x[2])))*(-2.0*x[2]*f[1]+x[1]*(f[1]+f[2]))

    #Quadratic extrapolation between 0-x1
#    fint[1]=(x[1]*(6*x[2]*(x[2]-x[3])*x[3]*f[1]+x[1]^3.0*(f[2]-f[3])+3*x[1]*(x[3]^2*(f[1]+f[2])-x[2]^2.0*(f[1]+f[3]))+2*x[1]^2.0*(-x[3]*(f[1]+2*f[2])+x[2]*(f[1]+2*f[3]))))/(6*(x[1]-x[2])*(x[1]-x[3])*(x[2]-x[3]))
#    fint[1]=f[1]

    #Simpsons rule for everything else
    for i in 1:n-2

        d31=x[i+2]-x[i]
        d32=x[i+2]-x[i+1]
        d21=x[i+1]-x[i]

        x3=x[i+2]
        x2=x[i+1]
        x1=x[i]
        y3=f[i+2]
        y2=f[i+1]
        y1=f[i]

        integ=(d21/(6.0*d31*d32))*((3.0*x3^2.0)*(y1+y2)-2.0*x2*x3*(2.0*y1+y2)+(x2^2.0)*(y1-y3)+(x1^2.0)*(y2-y3)+2*x1*(x2*(y1+y2+y3)-x3*(y1+2.0*y2)))

#        println("i=$i $integ")

        fint[i+1]=fint[i]+abs(integ)

        if i==(n-2)
            integ=-(d32/(6.0*d21*d31))*((x3^2.0)*(y1-y2)+(x2^2.0)*(y1-y3)-(3.0*x1^2.0)*(y2+y3)+2.0*x1*x3*(2.0*y2+y3)-2.0*x2*(x3*(y1+y2+y3)-x1*(y2+2.0*y3)))
            fint[i+2]=fint[i+1]+integ
        end
    end

return fint
end

############################################
# smooth
############################################
#
# Smooth vector with Gaussian kernel
#
############################################
function smooth!(arr::AbstractVector, #input vector
                 N::Int=1; #number of smoothings
                 offs::Int=3) #offset of smoothing

    l = length(arr)
    offs = 3

    #Gaussian 7x1 kernel
    #XXX: accept other kernels too
    kernel = Float64[0.00443185,
                            0.05399100,
                            0.24197100,
                            0.39894200,
                            0.24197100,
                            0.05399100,
                            0.00443185]

    #XXX is this correct?
    logkernel = Float64[
                            0.00443185,
                            0.07342030,
                            0.27844800,
                            0.39894200,
                            0.16782600,
                            0.02927920,
                            0.00443185
                            ]

    #arr2 = zeros(l+2offs)
    for s = 1:N

        #Expand original array to avoid boundaries
        #XXX: extrapolate instead
        #arr2[1:offs] = ones(offs)*arr[1]
        #arr2[(offs+1):(end-offs)] = arr[:]
        #arr2[(end-offs+1):end] = ones(offs)*arr[l]
        #sarr = conv(arr2, kernel)
        #arr[:] = sarr[(1+offs+3):(end-offs-3)]

        #smooth only middle parts
        sarr = conv(arr, kernel)
        arr[(offs+1):(l-offs)] = sarr[(1+offs+3):(end-offs-3)]

    end

    return arr
end
smooth(arr, N=1; kvs...) = smooth!(deepcopy(arr), N; kvs...)

############################################
# smooth_spline
############################################
#
# Smooth vector with B-splines according to DeBoor's
# algorithm.
#
# Adapted from Fred Frigo's (Dec 8, 2001) Matlab source file
# that originates from http://www.psc.edu/~burkardt/src/splpak/splpak.f90
############################################
function smooth_spline(y::AbstractVector, dx::AbstractVector, smooth_factor::Real)

npoint = length(y)
p = smooth_factor

a = zeros(npoint, 4)
v = zeros(npoint, 7)

x = linspace(0.0, (npoint-1.0)/npoint, npoint)

# setupq

  v[1,4] = x[2]-x[1]

  for i = 2:npoint-1
    v[i,4] = x[i+1]-x[i]
    v[i,1] = dx[i-1]/v[i-1,4]
    v[i,2] = ((-1.0 .*dx[i])/v[i,4]) - (dx[i]/v[i-1,4])
    v[i,3] = dx[i+1]/v[i,4]
  end

  v[npoint,1] = 0.0
  for i = 2:npoint-1
    v[i,5] = (v[i,1]*v[i,1]) + (v[i,2]*v[i,2]) + (v[i,3]*v[i,3])
  end

  for i = 3:npoint-1
    v[i-1,6] = (v[i-1,2]*v[i,1]) + (v[i-1,3]*v[i,2])
  end

  v[npoint-1,6] = 0.0

  for i = 4: npoint-1
    v[i-2,7] = v[i-2,3]*v[i,1]
  end

  v[npoint-2,7] = 0.0
  v[npoint-1,7] = 0.0

#  Construct  q-transp. * y  in  qty.

  prev = (y[2]-y[1])/v[1,4]
  for i= 2:npoint-1
    diff = (y[i+1]-y[i])/v[i,4]
#    %qty(i) = diff-prev
    a[i,4] = diff - prev
    prev = diff
  end

# end setupq

#chol1d

#  Construct 6*(1-p)*q-transp .*(d**2)*q + p*r

  six1mp = 6.0 .*(1.0-p)
  twop = 2.0 .*p

  for i = 2: npoint-1
    v[i,1] = (six1mp .*v[i,5]) + (twop .*(v[i-1,4]) + v[i,4])
    v[i,2] = (six1mp .*v[i,6]) +( p .*v[i,4])
    v[i,3] = six1mp .*v[i,7]
  end

  if  npoint < 4
    u[1] = 0.0
    u[2] = a[2,4]/v[2,1]
    u[3] = 0.0

#  Factorization

  else
    for i = 2: npoint-2
      ratio = v[i,2]/v[i,1]
      v[i+1,1] = v[i+1,1]-(ratio .*v[i,2])
      v[i+1,2] = v[i+1,2]-(ratio .*v[i,3])
      v[i,2] = ratio
      ratio = v[i,3]./v[i,1]
      v[i+2,1] = v[i+2,1]-(ratio .*v[i,3])
      v[i,3] = ratio
    end

#  Forward substitution

    a[1,3] = 0.0
    v[1,3] = 0.0
    a[2,3] = a[2,4]
    for i = 2: npoint-2
      a[i+1,3] = a[i+1,4] - (v[i,2]*a[i,3]) - (v[i-1,3]*a[i-1,3])
    end

#%  Back substitution.

    a[npoint,3] = 0.0
    a[npoint-1,3] = a[npoint-1,3] / v[npoint-1,1]

    for i = npoint-2:-1:2
      a[i,3] = (a[i,3]/v[i,1]) - (a[i+1,3]*v[i,2]) - (a[i+2,3]*v[i,3])
    end

  end
#
#  Construct Q*U.
#
  prev = 0.0
  for i = 2: npoint
    a[i,1] = (a[i,3]-a[i-1,3])/v[i-1,4]
    a[i-1,1] = a[i,1]-prev
    prev = a[i,1]
  end

  a[npoint,1] = -1.0 .*a[npoint,1]

#end chol1d

  spline_sig = zeros(npoint)
  for i = 1: npoint
    spline_sig[i] = y[i]-(6.0 .*(1.0-p) .*dx[i] .*dx[i] .*a[i,1])
  end

    return spline_sig
end

#no weights
smooth_spline(y::AbstractVector, smooth_factor::Real) = smooth_spline(y, ones(length(y)), smooth_factor)

############################################
# smooth_plaw
############################################
#
# Smooth vector that has powerlaw behaviour by
# smoothing the ratio to some reference vector
# with either gaussian kernel or splines.
#
############################################
function smooth_plaw(x::AbstractVector,
                     ref::AbstractVector, #reference vector
                     N::Real=1; #number of smoothings OR smoothing factor
                     offs::Int=3, #offset of smoothing
                     method=:kernel) #:kernel, :spline

    @assert length(x) == length(ref)

    dr = x./ref

    if method == :kernel
        sdr = smooth(dr, int(N), offs=offs)
    elseif method == :spline
        sdr = smooth_spline(dr, N)
    else
        error("unrecognized method $(string(method))")
    end

    sx = ref .* sdr

    return sx
end
############################################
# gauss_laquerre_nw
############################################
#
# Gauss-Laguerra quadrature
#
# Calculates weights and nodes of Gaussian quadrature
# in interval [0,infty) with weight exp(-x)
#
# output: r (vector of nodes)
#         w (vector of weights)
############################################
function gauss_laguerre_nw(num::Real) #number of points
    #J=diagm([1:2:(2*num-1)])+diagm([1:num-1],1)+diagm([1:num-1],-1)
    J=diagm(0 => [1:2:(2*num-1)])+diagm(1 => [1:num-1])+diagm(-1 => [1:num-1])
    #J=spdiagm(sparsevec)+diagm(Pair{1,[1:num-1]})+diagm(Pair{-1,[1:num-1]})
    R,W=eig(J)
    r=sort(R)
    w=(W[1,:].^2.0)

    return r,w
end

############################################
# gauss_legendre_nw
############################################
#
# Gauss-Legendre quadrature
#
# calculates weights and nodes of Gaussian quadrature
# in interval [-1,1] with weight 1
#
# output: r (vector of nodes)
#         w (vector of weights)
############################################
function gauss_legendre_nw(num::Real)

    nlist=[1:num-1]
    J=diagm(nlist./sqrt(4 .*nlist.^2.0-1),1)+diagm(nlist./sqrt(4 .*nlist.^2.0-1),-1)
    R,W=eig(J)
    r=sort(R)
    w=2.0*(W[1,:].^2.0)

    return r,w
end

############################################
# expi
############################################
#
# Exponential integral for positive agruments
# after Cody and Thacher, Math. of Comp. 22,641 (1968)
#
############################################
function expi(n::Real,
              x::Real)

    @assert x >= 0.0

    #Table of predefined values
    X1=-1.0e20

    A0=-44178.5471728217
    A1=57721.7247139444
    A2=9938.31388962037
    A3=1842.11088668000
    A4=101.093806161906
    A5=5.03416184097568
    B0=76537.3323337614
    B1=32597.1881290275
    B2=6106.10794245759
    B3=635.419418378382
    B4=37.2298352833327

    C0=4.65627107975096e-07
    C1=0.999979577051595
    C2=9.04161556946329
    C3=24.3784088791317
    C4=23.0192559391333
    C5=6.90522522784444
    C6=0.430967839469389
    D1=10.0411643829054
    D2=32.4264210695138
    D3=41.2807841891424
    D4=20.4494785013794
    D5=3.31909213593302
    D6=0.103400130404874

    E0=-0.999999999998447
    E1=-26.6271060431811
    E2=-241.055827097015
    E3=-895.927957772937
    E4=-1298.85688746484
    E5=-545.374158883133
    E6=-5.66575206533869
    F1=28.6271060422192
    F2=292.310039388533
    F3=1332.78537748257
    F4=2777.61949509163
    F5=2404.01713225909
    F6=631.657483280800

    if x==X1
#        expint=ex1
        expint=0.0
        if n==1
            return expint
        end
        n1=n-1
        for i in 1:n1
            expint=(ex-x*expint)/(1.0*i)
        end
        return expint
    end

    ex=exp(-x)
    x1=x
    if x>4.0
        ex1=(ex+ex*(E0+(E1+(E2+(E3+(E4+(E5+E6/x)/x)/x)/x)/x)/x)/(x+F1+(F2+(F3+(F4+(F5+F6/x)/x)/x)/x)/x))/x
    elseif x>1.0
        ex1=ex*(C6+(C5+(C4+(C3+(C2+(C1+C0*x)*x)*x)*x)*x)*x)/(D6+(D5+(D4+(D3+(D2+(D1+x)*x)*x)*x)*x)*x)
    elseif x>0.0
        ex1=(A0+(A1+(A2+(A3+(A4+A5*x)*x)*x)*x)*x)/(B0+(B1+(B2+(B3+(B4+x)*x)*x)*x)*x)-log(x)
    else
        ex1=0.0
    end

    expint=ex1
    if n==1
        return expint
    end

    for i in 1:n-1
        expint=(ex-x*expint)/(1.0*i)
    end

    return expint
end

max_stars_repo_name: natj/toolbox

max_issues_repo_name: natj/toolbox

max_forks_repo_name: natj/toolbox


text:
import .Plots

Image(height, width, a::T) where {T<:Union{Plots.Plot,Plots.Subplot}} =
    Image(height, width, ImageFileData(height, width, a))

function ImageFileData(height, width, a::T) where {T<:Union{Plots.Plot,Plots.Subplot}}
    filename = tempname()*".pdf"
    Plots.pdf(a, filename)
    r = read(filename)
    rm(filename)
    ImageFileData(r, :pdf)
end


text:
using SuccessiveConvexProgrammings
const SCvxs = SuccessiveConvexifications
using Test

using LinearAlgebra
using Convex
const Cvx = Convex
# using Debugger
using Plots
ENV["GKSwstype"]="nul"  # do not show plot

function print_message(name::String)
    println(">"^6 * name * "<"^6)
end

function test_initial_guess()
    print_message("initial_guess")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u)
    X_0 = scvx.X_k
    U_0 = scvx.U_k
    X_k = zeros(N, n_x)
    U_k = ones(N-1, n_u)
    initial_guess!(scvx, X_k, U_k)
    @test norm(scvx.X_k - X_0) == 0.0
    @test norm(scvx.U_k - U_0) != 0.0
end

################## get_J ##################
function test_get_J()
    print_message("get_J")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                objs_path=[my_path_obj],
                objs_terminal=[my_terminal_obj],
                consts_path_ineq=[not_too_large_input],
                consts_initial_eq=[my_const_initial_eq],
                consts_terminal_eq=[my_const_terminal_eq],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    obj = SCvxs.get_J(scvx, X, U)
    @test obj == (
                  16.0
                  + (N-1) * scvx.λ * 2.0 * n_u
                  + 0.0 + scvx.λ * norm(my_const_terminal_eq(X[end, :]), 1)
                 )
end

## get_obj
function test_get_obj()
    print_message("get_obj")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                objs_path=[my_path_obj],
                objs_terminal=[my_terminal_obj],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    obj = SCvxs.get_obj(scvx, X, U)
    @test obj == 16.0
end

function test_calculate_objs_path()
    print_message("calculate_objs_path")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                objs_path=[my_path_obj],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    path_obj = SCvxs.calculate_objs_path(scvx, X, U)
    @test path_obj == 0.0
end

function test_calculate_objs_terminal()
    print_message("calculate_objs_terminal")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                objs_terminal=[my_terminal_obj],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    terminal_obj = SCvxs.calculate_objs_terminal(scvx, X)
    @test terminal_obj == 16.0
end

## get_obj_const_penalty
function test_get_obj_const_penalty()
    print_message("get_obj_const_penalty")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                consts_path_ineq=[not_too_large_input],
                consts_initial_eq=[my_const_initial_eq],
                consts_terminal_eq=[my_const_terminal_eq],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    obj = SCvxs.get_obj_const_penalty(scvx, X, U)
    @test obj == (N-1) * scvx.λ * 2.0 * n_u + 0.0 + scvx.λ * norm(my_const_terminal_eq(X[end, :]), 1)
end

function test_calculate_const_path_penalty()
    print_message("calculate_const_path_penalty")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                consts_path_ineq=[not_too_large_input],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    obj = SCvxs.calculate_const_path_penalty(scvx, X, U)
    @test obj == (N-1) * scvx.λ * 2.0 * n_u
end

function test_calculate_const_initial_penalty()
    print_message("calculate_const_initial_penalty")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                consts_initial_eq=[my_const_initial_eq],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    obj = SCvxs.calculate_const_initial_penalty(scvx, X, U)
    @test obj == 0.0
end

function test_calculate_const_terminal_penalty()
    print_message("calculate_const_terminal_penalty")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                consts_terminal_eq=[my_const_terminal_eq],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    obj = SCvxs.calculate_const_terminal_penalty(scvx, X)
    @test obj == scvx.λ * norm(my_const_terminal_eq(X[end, :]), 1)
end

################## get_L ##################
function test_get_L()
    print_message("get_L")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                objs_path=[my_path_obj],
                objs_terminal=[my_terminal_obj],
                consts_path_ineq=[not_too_large_input],
                consts_initial_eq=[my_const_initial_eq],
                consts_terminal_eq=[my_const_terminal_eq],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    D, W = X, U
    jacob_dict = Dict()
    obj = SCvxs.get_L(scvx, D, W, jacob_dict)
    @test obj == (
                  16.0
                  + (N-1) * scvx.λ * 2.0 * n_u
                  + 0.0 + scvx.λ * norm(my_const_terminal_eq(X[end, :]), 1)
                 )
end

## get_obj_linearised
function test_get_obj_linearised()
    print_message("get_obj_linearised")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                objs_path=[my_path_obj],
                objs_terminal=[my_terminal_obj],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    obj, jacob_dict = SCvxs.get_obj_linearised(scvx, X, U)
    @test obj == 0.0 + 16.0
end

function test_calculate_objs_path_linearised()
    print_message("calculate_objs_path_linearised")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                objs_path=[my_path_obj],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    D, W = X, U  # originally, it should be values from Convex.Variable
    jacob_dict = Dict()
    obj, _ = SCvxs.calculate_objs_path_linearised(scvx, D, W, jacob_dict)
    @test obj == 0.0
end
function test_calculate_objs_terminal_linearised()
    print_message("calculate_objs_terminal_linearised")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                objs_terminal=[my_terminal_obj],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    D, W = X, U  # originally, it should be values from Convex.Variable
    jacob_dict = Dict()
    obj, _ = SCvxs.calculate_objs_terminal_linearised(scvx, D, W, jacob_dict)
    @test obj == 16.0
end

## get_obj_const_penalty_linearised
function test_get_obj_const_penalty_linearised()
    print_message("get_obj_const_penalty_linearised")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                consts_path_ineq=[not_too_large_input],
                consts_path_eq=[my_const_path_eq],
                consts_initial_eq=[my_const_initial_eq],
                consts_terminal_eq=[my_const_terminal_eq],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    obj, jacob_dict = SCvxs.get_obj_const_penalty_linearised(scvx, X, U)
    @test obj == (N-1) * scvx.λ * 2.0 * n_u + 0.0 + scvx.λ * 2.0 * n_x
end

function test_calculate_const_path_penalty_linearised()
    print_message("calculate_const_path_penalty_linearised")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                consts_path_ineq=[not_too_large_input],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    D, W = X, U  # originally, it should be values from Convex.Variable
    jacob_dict = Dict()
    obj, _ = SCvxs.calculate_const_path_penalty_linearised(scvx, D, W, jacob_dict)
    @test obj == (N-1) * scvx.λ * 2.0 * n_u
end

function test_calculate_const_initial_penalty_linearised()
    print_message("calculate_const_initial_penalty_linearised")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                consts_initial_eq=[my_const_initial_eq],
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    D, W = X, U  # originally, it should be values from Convex.Variable
    jacob_dict = Dict()
    obj, _ = SCvxs.calculate_const_initial_penalty_linearised(scvx, D, W, jacob_dict)
    @test obj == 0.0
end

function test_calculate_const_terminal_penalty_linearised()
    print_message("calculate_const_terminal_penalty_linearised")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                consts_terminal_eq=[my_const_terminal_eq]
               )
    X, U = scvx.X_k, scvx.U_k  # zeros
    D, W = X, U  # originally, it should be values from Convex.Variable
    jacob_dict = Dict()
    obj, _ = SCvxs.calculate_const_terminal_penalty_linearised(scvx, D, W, jacob_dict)
    @test obj == 2.0 * n_x * scvx.λ
end

################## get_L ##################

function test_solve_cvx_subprob()
    print_message("solve_cvx_subprob")
    n_x, n_u, N = 4, 2, 11
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                objs_path=[my_path_obj],
                objs_terminal=[my_terminal_obj],
                consts_path_ineq=[not_too_large_input],
                consts_path_eq=[my_const_path_eq],
                consts_initial_eq=[my_const_initial_eq],
                consts_terminal_eq=[my_const_terminal_eq],
               )
    D, W, jacob_dict = SCvxs.solve_cvx_subprob(scvx)
    println("D: $(D.value)")
    println("W: $(W.value)")
end

function test_scvx_example(; verbose=false)
    n_x, n_u, N = 4, 2, 31
    scvx = SCvx(N=N, n_x=n_x, n_u=n_u,
                objs_path=[my_path_obj],
                objs_terminal=[my_terminal_obj],
                consts_path_ineq=[not_too_large_input2,
                                  not_too_small_input],
                consts_path_eq=[my_const_path_eq],
                consts_initial_eq=[my_const_initial_eq],
                # consts_terminal_eq=[my_const_terminal_eq],
               )
    append!(scvx, :consts_terminal_eq, [my_const_terminal_eq])
    X_0, U_0 = ones(N, n_x), ones(N-1, n_u)
    scvx = initial_guess!(scvx, X_0, U_0)
    @time solve!(scvx, verbose=verbose)
end

function test_solve(;verbose=false)
    print_message("solve")
    scvx = test_scvx_example(verbose=false)
    N = scvx.N
    @show scvx.obj_k, scvx.i, scvx.flag
    @show scvx.X_k
    @show scvx.U_k
    # initial condition
    tol = 1e-4
    @test norm(my_const_initial_eq(scvx.X_k[1, :], scvx.U_k[1, :])) < tol
    @test norm(my_const_terminal_eq(scvx.X_k[end, :])) < tol
    for n in 1:N-1
        @test all(not_too_large_input2(scvx.X_k[n, :],
                                       scvx.U_k[n, :]) .<= tol)
        @test norm(my_const_path_eq(scvx.X_k[n, :],
                                    scvx.U_k[n, :],
                                    scvx.X_k[n+1, :])) < tol
    end
    # plot
    data = Dict("time_state" => collect(1:N),
                "time_input" => collect(1:N-1),
                "state" => scvx.X_k_sol, "input" => scvx.U_k_sol)
    @show data["state"]
    @show data["input"]
    _plot(data)
    return scvx
end

function test_solve_and_flush(; verbose=false)
    print_message("solve and flush")
    scvx = test_solve(; verbose=verbose)
    X_k = deepcopy(scvx.X_k)
    U_k = deepcopy(scvx.U_k)
    flush!(scvx)
    @test norm(X_k[2:end, :] - scvx.X_k) == 0
    @test norm(U_k[2:end, :] - scvx.U_k) == 0
end

function _plot(data)
    log_dir = "data/test"
    mkpath(log_dir)
    for name in ["state", "input"]
        time = data["time_" * name]
        p = plot(1)
        for i in 1:length(data[name])
            if i == length(data[name])
                color = :purple
                lw = 3.0
            else
                color = :blue
                lw = 1.0
            end
            plot!(p, time, data[name][i],
                  linecolor=color, lw=lw, alpha=0.5*i/length(data[name]),
                  label="")
        end
        if name == "state"
            plot!(seriestype=:scatter,
                  [time[1]],
                  [zeros(size(data[name][1, :]))],
                  markercolor="red",
                  label="init cond",
                 )  # init cond
            plot!(seriestype=:scatter,
                  [time[end]],
                  [2*ones(size(data[name][end, :]))],
                  markercolor="red",
                  label="term cond",
                 )  # term cond
        elseif name == "input"
            plot!(time, 2*ones(size(time)),  # not too large
                  linecolor="red",
                  label="not too large",
                 )
            plot!(time, -2*ones(size(time)),  # not too small
                  linecolor="red",
                  label="not too small",
                 )
        end
        savefig(p, joinpath(log_dir, name*".pdf"))
    end
end

################## custom functions ##################
function my_path_obj(x::Array, u::Array)::Array
    return [0.5*(norm(x)^2 + norm(u)^2)]
end

function my_terminal_obj(x::Array)::Array
    return [norm(x .- 2.0)^2]
end

function not_too_large_input(x::Array, u::Array)::Array
    return u .+ 2.0  # <=0
end

function not_too_large_input2(x::Array, u::Array)::Array
    return u .- 2.0  # <=0
end

function not_too_small_input(x::Array, u::Array)::Array
    return -(u .+ 2.0)  # <=0
end

function my_const_path_eq(x::Array, u::Array, x_next::Array)::Array
    A = Matrix(I, 4, 4)
    B = [0 1; 0 1; 1 0; 1 0]
    dynamics(x, u) = A*x + B*u
    return x_next - dynamics(x, u)
end

function my_const_initial_eq(x::Array, u::Array)
    return x
end

function my_const_terminal_eq(x::Array)
    return (x .- 2.0)
end

################## test all ##################
function test_all()
    # test_initial_guess()

    # ################## get_J ##################
    # test_get_J()
    # ## get_obj
    # test_get_obj()
    # test_calculate_objs_path()
    # test_calculate_objs_terminal()

    # ## get_obj_const_penalty
    # test_get_obj_const_penalty()
    # test_calculate_const_path_penalty()
    # test_calculate_const_initial_penalty()
    # test_calculate_const_terminal_penalty()

    # ################## get_L ##################
    # test_get_L()
    # ## get_obj_linearised
    # test_get_obj_linearised()
    # test_calculate_objs_path_linearised()
    # test_calculate_objs_terminal_linearised()
    # ## get_obj_const_penalty_linearised
    # test_get_obj_const_penalty_linearised()
    # test_calculate_const_path_penalty_linearised()
    # test_calculate_const_initial_penalty_linearised()
    # test_calculate_const_terminal_penalty_linearised()

    # ################## solve ##################
    # test_solve_cvx_subprob()
    # test_solve(verbose=false)
    test_solve_and_flush(verbose=false)
    # ################## solve ##################
    # test_scvx_example()
end

test_all()


text:
"""
Given the original LP of the form
    min cx
    subject to
    -Inf <= A1 x <= b1
      b2 <= A2 x <= Inf
      bl <= A3 x <= bu
      b4 <= A4 x <= b4
      xl <= xu <= xu,
the canonical form is given by
    min cx + 0 s1 + 0 s2 + 0 s3
    subject to
    A1 x + s1 = b1
    A2 x - s2 = b2
    A3 x - s3 = bl
    A4 x      = b4
    xl <= x <= xu
    0 <= s1 <= Inf
    0 <= s2 <= Inf
    0 <= s3 <= bu - bl
"""

function canonical_form(standard::MatOI.LPForm{T, AT, VT}) where {T, AT, VT}
    # count the number of inequality constraints
    ineq = Int[]
    for i = 1:nrows(standard)
        if standard.c_lb[i] < standard.c_ub[i]
            push!(ineq, i)
        end
    end
    nineq = length(ineq)

    lp_nrows = nrows(standard)
    lp_ncols = ncols(standard) + nineq

    # objective coefficient
    c = append!(deepcopy(Array(standard.c)), zeros(nineq))
    @assert length(c) == lp_ncols

    # column bounds
    xl = append!(deepcopy(Array(standard.v_lb)), zeros(nineq))
    xu = append!(deepcopy(Array(standard.v_ub)), fill(Inf, nineq))
    @assert length(xl) == lp_ncols
    @assert length(xu) == lp_ncols
    for i = 1:nineq
        if standard.c_lb[ineq[i]] > -INF
            xu[ncols(standard) + i] = standard.c_ub[ineq[i]] - standard.c_lb[ineq[i]]
        end
    end

    # row bounds
    b = deepcopy(Array(standard.c_lb))

    # create the submatrix for slack
    I = Int64[]; J = Int64[]; V = Float64[]; j = 1;
    for i in ineq
        if standard.c_lb[i] > -INF
            push!(I, i); push!(J, j); push!(V, -1.0);
            j += 1
        elseif standard.c_ub[i] < INF
            b[i] = standard.c_ub[i]
            push!(I, i); push!(J, j); push!(V, 1.0);
            j += 1
        end
    end
    S = sparse(I, J, V, nrows(standard), nineq)

    # create the matrix in canonical form
    @assert size(standard.A,1) == size(S,1)
    A = [sparse(Matrix(standard.A)) S]

    # create a canonical form data
    # @show standard.TArray
    senses = fill(MatOI.EQUAL_TO, lp_nrows)
    canonical = MatOI.LPSolverForm{T, typeof(A), typeof(c)}(
        MOI.MIN_SENSE, c, A, b, senses, xl, xu
    )

    return canonical
end
canonical_form(standard::MatOI.LPSolverForm) = standard

function cpu2gpu(lp::MatOI.LPForm{T, AT, VT}) where {T, AT, VT}
    if VT <: CuArray
        return lp
    end
    c = convert(CuArray{T, 1}, lp.c)
    A = convert(CuArray{T, 2}, lp.A)
    c_lb = convert(CuArray{T, 1}, lp.c_lb)
    c_ub = convert(CuArray{T, 1}, lp.c_ub)
    v_lb = convert(CuArray{T, 1}, lp.v_lb)
    v_ub = convert(CuArray{T, 1}, lp.v_ub)

    return MatOI.LPForm{T, typeof(A), typeof(b)}(
        lp.direction, c, A, c_lb, c_ub, v_lb, v_ub
    )
end

function cpu2gpu(lp::MatOI.LPSolverForm{T, AT, VT}) where {T, AT, VT}
    if VT <: CuArray
        return lp
    end
    c = convert(CuArray{T, 1}, lp.c)
    A = convert(CuArray{T, 2}, lp.A)
    b = convert(CuArray{T, 1}, lp.b)
    v_lb = convert(CuArray{T, 1}, lp.v_lb)
    v_ub = convert(CuArray{T, 1}, lp.v_ub)

    return MatOI.LPSolverForm{T, typeof(A), typeof(b)}(
        lp.direction, c, A, b, lp.senses, v_lb, v_ub
    )
end


text:
module ExtensibleEffects
export effect, noeffect, NoEffect,
  runhandlers, @runhandlers,
  @insert_into_runhandlers,
  @syntax_eff, @syntax_eff_noautorun, noautorun,
  WriterHandler,
  ContextManagerHandler, @runcontextmanager, @runcontextmanager_, ContextManagerCombinedHandler,
  CallableHandler, @runcallable,  # Callable handler
  StateHandler, @runstate

using Compat
using Reexport
using TypeClasses

# re-export all DataTypes
@reexport using TypeClasses.DataTypes
@reexport using DataTypesBasic
export @pure  # users need @pure for the monadic syntax @syntax_eff

include("core.jl")
include("effecthandler.jl")
include("outereffecthandler.jl")
include("autorun.jl")
include("syntax.jl")
include("instances.jl")

end # module

max_stars_repo_name: schlichtanders/ExtensibleEffects

max_issues_repo_name: schlichtanders/ExtensibleEffects

max_forks_repo_name: schlichtanders/ExtensibleEffects


text:
using Fed
using Flux
using SHA

struct Server
    model::Chain
    dataset::DataLoader

    Server() = new(
        newmodel(),
        get_testdata(isflat=true)
    )
end

function start_server()
    server = Server()
    config = newconfig(9090, false)

    # config
    host = "127.0.0.1"
    port = 8080
    weights = flatten_model(server.model)
    strategy = Fed.Server.federated_averaging
    eval_hook = (weights::Vector{Float32}) -> evaluate(server, weights)

    central_node = Fed.Server.CentralNode(
        host,
        port,
        weights,
        strategy,
        eval_hook,
        config
    )

    # start the server
    @info "Server started on [http://$host:$port]"
    Fed.Server.start(central_node)
end


text:
using Documenter, ErrorAnalysis

makedocs(
    # options

)

makedocs(
    modules = [ErrorAnalysis],
    format = :html,
    sitename = "ErrorAnalysis.jl",
    pages = [
        "Home" => "index.md",
        "Functions" => "functions.md"
        # "Subsection" => [
        #     ...
        # ]
    ]
)

deploydocs(
    repo   = "github.com/crstnbr/ErrorAnalysis.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing,
    julia  = "release",
    osname = "linux"
)


text:
module Blosc2
include("Lib/lib.jl")
include("common/common.jl")
include("compression/compression.jl")
include("init.jl")
export Compressor, available_compressors, available_compressors_names, has_compressor, max_compression_overhead
export compressor_by_name, default_compressor_name, default_compressor
export CompressionParams, DecompressionParams, Context, CompressionContext, DecompressionContext
export compress!, compress, unsafe_compress!, decompress!, decompress, unsafe_decompress!, decompress_items!
export sizes, uncompressed_sizeof, uncompressed_length, compressed_sizeof, max_compressed_size, make_compress_buffer
export filter_pipeline, has_filter, available_filter_names, filter_description, max_filters_count
export Lib
end # module


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # VV10 Non-local correlation kernel
# One of the largest deficiencies of semilocal functionals is the lack of long-range correlation effects. The most notable effect is the lack of dispersion in the interactions between molecules. VV10 was expressly created to bridge the gap between the expensive true non-local correlation and a computational tractable form. We will begin by writing the overall expression:
#
# $$E_c^{\rm{nl}} = \frac{1}{2}\int \int d{\bf r}d{\bf r'}\rho({\bf r})\Phi({\bf r},{\bf r'})\rho({\bf r'}),$$
#
# where the two densities are tied together through the $\Phi$ operator.
#
# For VV10 we have:
# $$
# \begin{align}
# \Phi &= -\frac{3}{2gg'(g + g')},\\
#  g &= \omega_0({\rm r}) R^2 + \kappa({\rm r)},\\
#  g' &= \omega_0({\rm r}) R^2 + \kappa({\rm r')},
# \end{align}
# $$
#
# where $\omega_{0}$:
#
# $$
# \begin{align}
# \omega_{0}(r) &= \sqrt{\omega_{g}^2(r) + \frac{\omega_p^2(r)}{3}} \\
# \omega_g^2(r) &= C \left | \frac{\nabla \rho({\bf r})}{\rho({\bf r})} \right |^4 \\
# \omega_p^2(r) &= 4 \pi \rho({\bf r}),
# \end{align}
# $$
#
# and finally:
#
# $$\kappa({\bf r}) = b * \frac{3 \pi}{2} \left [ \frac{\rho({\bf r})}{9\pi} \right ]^\frac{1}{6}.$$
#
# While there are several expressions, this is actually relatively easy to compute. First let us examine how the VV10 energy is reintegrated:
#
# $$E_c^{\rm{VV10}} = \int d{\bf r} \rho{\bf r} \left [ \beta + \frac{1}{2}\int d{\bf r'} \rho{\bf r'} \Phi({\bf r},{\bf r'}) \right].$$
#
#

# +
using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy") # used only to cast to Psi4 arrays
using LinearAlgebra: dot
import Formatting: printfmt
include("ks_helper.jl")

mol = psi4.geometry("""
He 0 0 -5
He 0 0  5
symmetry c1
""")
options = (Dict("BASIS" => "aug-cc-pvdz",
                "DFT_SPHERICAL_POINTS" => 110,
                "DFT_RADIAL_POINTS" => 20))
# -

# ## VV10 coefficients
# First let us build set and build a few static coefficients:

coef_C = 0.0093
coef_B = 5.9
β = 1.0 / 32.0 * (3.0 / (coef_B^2.0))^(3.0 / 4.0)

# ## VV10 kernel
# First let us construct a function that computes $\omega_0$, and $\kappa$ quantities. We simplify the following piece contained in $\omega_g$:
# $$\left |\frac{\nabla \rho({\bf r})}{\rho({\bf r})} \right|^4$$
#
# by recalling that
#
# $$\gamma({\bf r}) = \nabla\rho({\bf r})\cdot\nabla\rho({\bf r}),$$
#
# therefore, we can simplify the above to:
#
# $$\left |\frac{\nabla \rho({\bf r})}{\rho({\bf r})} \right |^4 = \left | \frac{\gamma({\bf r})}{\rho({\bf r})\cdot \rho({\bf r})} \right | ^2 $$

function compute_vv10_kernel(ρ, γ)
    κ_pref = coef_B * (1.5 * π) / ((9.0 * π)^(1.0 / 6.0))

    # Compute R quantities
    Wp = (4.0 / 3.0) * π * ρ
    Wg = @. coef_C * ((γ / (ρ * ρ))^2.0)
    W0 = @. sqrt(Wg + Wp)

    κ = @. ρ^(1.0 / 6.0) * κ_pref

    W0, κ
end

# ## VV10 energy and gradient evaluation
#
# The next block of code computes the VV10 energy and its gradient. In the very end we plug this function in
# the Khon-Sham solver `ks_solver` that is in a separate file `ks_helper.jl`. We conveniently separated the solver to focus in the non-local kernel.

# +
function compute_vv10(D, Vpot)

    nbf = D.shape[1]
    Varr = zeros(nbf, nbf)

    total_e = 0.0
    tD = 2np.array(D)

    points_func = Vpot.properties()[1]
    superfunc = Vpot.functional()

    xc_e = 0.0
    vv10_e = 0.0

    # First loop over the outer set of blocks
    for l_block in 1:Vpot.nblocks()

        # Obtain general grid information
        l_grid = Vpot.get_block(l_block-1)
        l_w = np.array(l_grid.w())
        l_x = np.array(l_grid.x())
        l_y = np.array(l_grid.y())
        l_z = np.array(l_grid.z())
        l_npoints = size(l_w,1)

        points_func.compute_points(l_grid)

        # Compute the functional itself
        ret = superfunc.compute_functional(points_func.point_values(), -1)

        xc_e += dot(l_w, np.array(ret["V"])[1:l_npoints])
        v_ρ = np.array(ret["V_RHO_A"])[1:l_npoints]
        v_γ = np.array(ret["V_GAMMA_AA"])[1:l_npoints]

        # Begin VV10 information
        l_ρ = np.array(points_func.point_values()["RHO_A"])[1:l_npoints]
        l_γ = np.array(points_func.point_values()["GAMMA_AA"])[1:l_npoints]

        l_W0, l_κ = compute_vv10_kernel(l_ρ, l_γ)

        phi_kernel = zero(l_ρ)
        phi_U = zero(l_ρ)
        phi_W = zero(l_ρ)

        # Loop over the inner set of blocks
        for r_block in 1:Vpot.nblocks()

            # Repeat as for the left blocks
            r_grid = Vpot.get_block(r_block-1)
            r_w = np.array(r_grid.w())
            r_x = np.array(r_grid.x())
            r_y = np.array(r_grid.y())
            r_z = np.array(r_grid.z())
            r_npoints = size(r_w,1)

            points_func.compute_points(r_grid)

            r_ρ = np.array(points_func.point_values()["RHO_A"])[1:r_npoints]
            r_γ = np.array(points_func.point_values()["GAMMA_AA"])[1:r_npoints]

            r_W0, r_κ = compute_vv10_kernel(r_ρ, r_γ)

            newaxis = [CartesianIndex()]

            # Build the distance matrix
            R2  = (l_x[:, newaxis] .- r_x').^2
            R2 += (l_y[:, newaxis] .- r_y').^2
            R2 += (l_z[:, newaxis] .- r_z').^2

            # Build g
            g  = @. l_W0[:, newaxis] * R2 + l_κ[:, newaxis]
            gp = @. r_W0' * R2 + r_κ'

            F_kernel = @. -1.5(r_w * r_ρ)' / (g * gp * (g + gp))
            F_U = @. F_kernel * ((1.0 / g) + (1.0 / (g + gp)))
            F_W = F_U .* R2

            phi_kernel += sum(F_kernel, dims=2)
            phi_U += -sum(F_U, dims=2)
            phi_W += -sum(F_W, dims=2)
        end

        # Compute those derivatives
        κ_dn = l_κ ./ (6l_ρ)
        w0_dγ = @. coef_C * l_γ / (l_W0 * l_ρ^4.0)
        w0_dρ = @. 2.0 / l_W0 * (π/3.0 - coef_C * l_γ^2.0 / (l_ρ^5.0))

        # Sum up the energy
        vv10_e += sum(@. l_w * l_ρ * (β + 0.5phi_kernel))

        # Perturb the derivative quantities
        v_ρ += @. β + phi_kernel + l_ρ * (κ_dn * phi_U + w0_dρ * phi_W)
        v_ρ *= 0.5

        v_γ += l_ρ .* w0_dγ .* phi_W

        # Recompute to l_grid
        lpos = np.array(l_grid.functions_local_to_global()) .+ 1
        points_func.compute_points(l_grid)
        nfunctions = size(lpos,1)

        # Integrate the LDA and GGA quantities
        phi = np.array(points_func.basis_values()["PHI"])[1:l_npoints, 1:nfunctions]
        phi_x = np.array(points_func.basis_values()["PHI_X"])[1:l_npoints, 1:nfunctions]
        phi_y = np.array(points_func.basis_values()["PHI_Y"])[1:l_npoints, 1:nfunctions]
        phi_z = np.array(points_func.basis_values()["PHI_Z"])[1:l_npoints, 1:nfunctions]

        # LDA
        # Vtmp[ab] = ϕ[pb] Vρ[p] lw[p] ϕ[pa]
        Vtmp = phi' * (v_ρ .* l_w .* phi)

        # GGA
        l_ρ_x = np.array(points_func.point_values()["RHO_AX"])[1:l_npoints]
        l_ρ_y = np.array(points_func.point_values()["RHO_AY"])[1:l_npoints]
        l_ρ_z = np.array(points_func.point_values()["RHO_AZ"])[1:l_npoints]

        tmp_grid = 2l_w .* v_γ
        # Vtmp[ab] += ϕ[pb] tmp_grid[p] ∇ₓρ[p] ∇ₓϕ[pa] ∀ x,y,z
        Vtmp += phi' * (tmp_grid .* l_ρ_x .* phi_x)
        Vtmp += phi' * (tmp_grid .* l_ρ_y .* phi_y)
        Vtmp += phi' * (tmp_grid .* l_ρ_z .* phi_z)

        # Sum back to the correct place
        Varr[lpos, lpos] += Vtmp + Vtmp'
    end

    printfmt("   VV10 NL energy: {:16.8f}\n", vv10_e)

    xc_e += vv10_e
    return xc_e, Varr
end

ks_solver("VV10", mol, options, compute_vv10)
# -

# Refs:
#  - Vydrov O. A.; Van Voorhis T., *J. Chem. Phys.*, **2010**, *133*, 244103


text:
using HilbertSpaceFillingCurve
using Test

@testset "hilbert" begin

d = 10
for ndims in 2:3, nbits in [8,16]
    p = hilbert(d, ndims, nbits)
    @test d == hilbert(p, ndims, nbits)
end

@test_throws AssertionError hilbert(d, 2, 64)

end


text:
@doc Markdown.doc"""
    lift(a::padic) -> fmpz

Returns the positive canonical representative in $\mathbb{Z}$. $a$ needs
to be integral.
"""
function lift(a::padic)
  b = fmpz()
  R = parent(a)
  ccall((:padic_get_fmpz, libflint), Nothing, (Ref{fmpz}, Ref{padic}, Ref{FlintPadicField}), b, a, R)
  return b
end

function _lift(a::padic)
  R = parent(a)
  v = valuation(a)
  if v >= 0
    return fmpq(lift(a))
  else
    m = prime(R)^-v
    return fmpq(lift(m * a))//m
  end
end

function Base.setprecision(f::Generic.Poly{padic}, N::Int)
  f = deepcopy(f)
  for i=1:length(f)
    f.coeffs[i].N = N
  end
  return f
end

function setprecision!(f::Generic.Poly{padic}, N::Int)
  for i=1:length(f)
    f.coeffs[i].N = N
  end
  return f
end