text:
@mix @columns struct LC50{LC}
    # Field           | Default  | Flat  | Bounds      | Description
    LC50::LC          | 10.0     | true  | (0.0, 1e9)  | "Lethal concentration for 50% of individuals."
end

@mix @columns struct HillCoefficient{HC}
    # Field                 | Default  | Flat  | Bounds      | Description
    hillcoefficient::HC     | 0.1      | true  | (0.0, 1e9) | "Hill coefficient, or shape of a log logistic function"
end

@mix @columns struct Exposure{XP}
    # Field           | Default  | Flat  | Bounds      | Description
    exposure::XP      | 100.0      | true  | (0.0, 1e9) | "Exposure: environmental concentration to which the population is exposed to"
end

"""
Extends CellRule for rules of survival effect

For best performance these should be chained with other
CellRule or following an NeighborhoodRule.
"""
abstract type SurvRule{R,W} <: CellRule{R,W} end

"""
Extends SurvRule for survival rules using a heterogenous
concentration layer.
"""

abstract type SurvMapRule{R,W} <: SurvRule{R,W} end

# Base.key(rule::SurvMapRule{R,W,K}) where {R,W,K} = K
# layer(rule::SurvMapRule, data) = axkey(rule)

DynamicGrids.precalcrules(rule::SurvMapRule, data) = begin
    precalclayer(layer(rule, data), rule, data)
end

DynamicGrids.precalcrules(rule::SurvRule, data) = rule

# Exact survival solutions

"""
Simple fixed survival loglogistic using exact solution.

Cumulative function of loglogistic: F(x; α, β) =  1/(1 + (x/ α)^{-β} ),
    where α>0 is the scale and β>0 is the shape.

Dose-response log-logistic survival model is defined as:
    Surv(x; α, β) = 1-F(x; α, β),
    where α is the LC50 (lethal concentration for 50% of individuals)
    and β is known as the Hill's coefficient, a measure of ultrasensitivity (i.e. how steep is the response curve)
    See [wikipedia on Hill equation](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry))

    Surv(x; α, β) = 1-F(x; α, β) = 1/(1 + (x/ α)^β )

$(FIELDDOCTABLE)
"""
@Exposure @LC50 @HillCoefficient struct SurvLogLogistic{R,W} <: SurvRule{R,W} end

@inline applyrule(data, rule::SurvLogLogistic, population, args...) = begin
    population > zero(population) || return zero(population)
    @fastmath population * ( 1/ (1 + (rule.exposure / rule.LC50)^rule.hillcoefficient) )
end

"""
LogLogistic survival based on an exposure layer
$(FIELDDOCTABLE)
"""
@Layers @LC50 @HillCoefficient struct SurvLogLogisticMap{R,W} <: SurvMapRule{R,W} end

@inline applyrule(data, rule::SurvLogLogisticMap, population, index, args...) = begin
    population > zero(population) || return zero(population)
    exposure = layer(rule, data, index)
    @fastmath population * ( 1/ (1 + (exposure / rule.LC50)^rule.hillcoefficient) )
end

@mix @columns struct LC503{LC1,LC2,LC3}
    # Field           | Default  | Flat  | Bounds      | Description
    LC501::LC1        | 10.0     | true  | (0.0, 1e9)  | "Lethal concentration for 50% of individuals."
    LC502::LC2        | 10.0     | true  | (0.0, 1e9)  | "Lethal concentration for 50% of individuals."
    LC503::LC3        | 10.0     | true  | (0.0, 1e9)  | "Lethal concentration for 50% of individuals."
end

@Layers @LC503 @HillCoefficient @CarryCap @InstrinsicGrowthRate struct GrowthSurvLogLogisticMap3{R,W} <: SurvMapRule{R,W} end

@inline applyrule(data, rule::GrowthSurvLogLogisticMap3, (population1,population2,population3), index, args...) = begin
    growthCarrycapEffect =  rule.intrinsicrate * (1-(population1+population2+population3)/rule.carrycap)
    exposure = layer(rule, data, index)
    @fastmath (
        population1 * ( 1/ (1 + (exposure / rule.LC501)^rule.hillcoefficient) ) * growthCarrycapEffect,
        population2 * ( 1/ (1 + (exposure / rule.LC502)^rule.hillcoefficient) ) * growthCarrycapEffect,
        population3 * ( 1/ (1 + (exposure / rule.LC503)^rule.hillcoefficient) ) * growthCarrycapEffect
    )
end

"""
LogLogistic survival based on an exposure layer
$(FIELDDOCTABLE)
"""
@Layers @HillCoefficient struct SurvLogLogisticMap2{R,W} <: SurvMapRule{R,W} end

@inline applyrule(data, rule::SurvLogLogisticMap2, population, index, args...) = begin
    population > zero(population) || return zero(population)
    exposure = layer(rule, data, index,1)
    LC50 = layer(rule, data, index,2)
    @fastmath population * ( 1/ (1 + (exposure / LC50)^rule.hillcoefficient) )
end

"""
Simple layer mask. Values below a certain threshold are replaced with zero.
$(FIELDDOCTABLE)
"""

#? Should we gather MaskSurvMap with MaskGrowthMap?
@Layers @Timestep struct MaskSurvMap{R,W,ST} <: SurvMapRule{R,W}
    threshold::ST | 0.5 | true  | (0.0, 1.0) | "Minimum viability index."
end

@inline applyrule(data, rule::MaskSurvMap, population, index, args...) =
    layer(rule, data, index) >= rule.threshold ? population : zero(population)


text:
using Dolo

# # model = Model("examples/models/rbc.yaml")
# model = Model("examples/models/consumption_savings_iid.yaml")

model = Model("experiments/rbc.yaml")

F = Dolo.Euler(model, ignore_constraints=false);

sol = improved_time_iteration(model)

z0 = SVector(model.calibration[:exogenous][1:2]...)
x0 = F.x0

F(x0,x0,z0,z0) - F(x0,x0)

r = F(F.x0, F.x0; set_future=false)

##### Compare with dolo.py
import Dolo: norm

import Dolo: MSM
import Dolo: df_A

f = u->F(u, F.x0; set_future=false, ignore_constraints=false)
df = u->df_A(F, u, F.x0; set_future=false)
v1 , j_num = Dolo.DiffFun(f, F.x0)

f0 = u->F(u, F.x0; set_future=false, ignore_constraints=true)
v1 , j_0 = Dolo.DiffFun(f0, F.x0)

j = Dolo.df_A(F, F.x0, F.x0)

println(v1.data[1])

println("j_0")
println(j_0.data[1])

println("j num")
println(j_num.data[1])

println("j phi")
println(j.data[1])

function PhiTest(u,v)
    sq = sqrt(u^2+v^2)
    p = u+v-sq
    return p
end

function PhiPhiTest(x, a, b)
    f = x
    y = PhiTest(f,x-a)
    z = PhiTest(-y, b-x)
    return z
end

fun(u) = PhiPhiTest(u,0, 1)
dfun(u) = (fun(u+0.000001) - fun(u))/0.000001

zvec = range(-1,1;length=100)
fvec = fun.(zvec)
dfvec = dfun.(zvec)

plot(zvec, fvec)

plot(zvec, dfvec)

l = [Dolo.PhiPhi(SVector(x), SMatrix{1,1}(1.0), SVector(x), SVector(0.0), SVector(1.0)) for x in zvec]
fvec2= [e[1][1] for e in l]
dfvec2= [e[2][1,1] for e in l]

plot(zvec, fvec)
plot!(zvec, fvec2)

plot(zvec, dfvec)
plot!(zvec, dfvec2)

using SimplePlots

(ddf - j).data[1]

function newton(fun, dfun, x0::MSM; maxit=50, dampen=1.0, verbose=false, bcksteps = 5)

    r0 = fun(x0)

    err = norm(r0)

    local x1

    x1data = deepcopy(x0.data)
    x1 = MSM(x1data, x0.sizes)

    for n=1:maxit

        r0 = fun(x0)
        err_0 = norm(r0)

        j = dfun(x0)
        δ = (j\r0)
        for i=0:(bcksteps-1)
            u = 0.5^i
            x1.data[:] .= x0.data - δ.data*u
            r1 = fun(x1)
            err_1 = norm(r1)
            if verbose
                println( "-    $i: ", err_1)
            end
            if err_1<err_0
                break
            end
        end

        if verbose
            println(n, " | ", norm(r0), " | ",  norm(δ))
        end

        x0 = x1

    end

    return x0

end

sol = newton(f, df, F.x0; maxit=15, verbose=true, bcksteps = 20)

norm(F(sol, F.x0))

F(sol, sol; set_future=true);

function loop(F, K)

    x0 = F.x0
    x1 = x0

    for k=1:K

        res = F(x0, x0; set_future=true)
        ε = norm(res)

        f = u->F(u, x0; set_future=false)
        df = u->df_A(F, u, x0; set_future=false)

        x1 = newton(f, df, x0; verbose=false)

        δ = x1 - x0
        η = norm(δ)

        println(ε, " : ", η)

        x0 = x1
    end

    return x1

end

loop(F, 100);

# tab = Dolo.tabulate(model, F.dr.dr, :k)
# plot(tab[:k], tab[:k])
# plot!(tab[:j], tab[:c])

using SimplePlots

tab = Dolo.tabulate(model, F.dr.dr, :w)
plot(tab[:w], tab[:w])
plot!(tab[:w], tab[:c])

# l = Dolo.df_B(F, F.x0, F.x0)

# Dolo.newton(u->F(u,F.x0), F.x0)

sol = Dolo.time_iteration(model, maxit=1000);
tab = tabulate(model, sol.dr, :m)
pl = plot(tab[:m], tab[:m])
plot!(pl, tab[:m], tab[:c])

Dolo.time_iteration(model; ignore_constraints=true);

Dolo.improved_time_iteration(model);

# grid, exo = Dolo.discretize(model)

println("Endo Grid size: $(Dolo.n_nodes(grid.endo))")
println("Exo Grid size: $(Dolo.n_nodes(grid.exo))")
println("Exo Grid I size: $(Dolo.n_inodes(exo, 1))")

grid, exo = Dolo.discretize(model;)

println("Endo Grid size: $(Dolo.n_nodes(grid.endo))")
println("Exo Grid size: $(Dolo.n_nodes(grid.exo))")
println("Exo Grid I size: $(Dolo.n_inodes(exo, 1))")

sol0 = time_iteration(model)

F = Dolo.Euler(model)

@time Dolo.loop_iti(F, F.x0)

@time Dolo.loop_ti(F, F.x0)

# μ = Dolo.ergodic_distribution(model, sol)

# using Plots

ext: jl


text:
import Projections, Test, Distributions, LinearAlgebra

include("tests.jl")

divergences = [Projections.Burg(),
               Projections.Hellinger(),
               Projections.ChiSquare(),
               Projections.ModifiedChiSquare(),
               Projections.KullbackLeibler()]

norms = [Projections.Linf(),
         Projections.Lone(),
         Projections.Ltwo()]

constraints = vcat(divergences, norms)

n  = 100
εs = vcat(1e-5:0.01:0.1, 0.2, 0.5, 1)

Test.@testset "All tests" begin

    Test.@testset "Generating functions" begin
        Test.@testset "$(typeof(d))" for d in divergences
            test_generate(d)
        end
    end

    Test.@testset "DRO" begin
        Test.@testset "ε = $ε" for ε in εs
            Test.@testset "Constraint = $(Projections.name(d))" for d in constraints
                q   = LinearAlgebra.normalize(rand(n), 1)
                c1  = rand(n)
                c2  = max.(min.(rand(n), 0.8), 0.2)

                Test.@testset "standard" begin
                    m = Projections.DRO(d, q, c1, ε)
                    test_model(m)
                end
                Test.@testset "sorted" begin
                    m = Projections.DRO(d, sort(q; rev = true), sort(c1), ε)
                    test_model(m)
                end
                Test.@testset "bounds" begin
                    m = Projections.DRO(d, q, c2, ε)
                    test_model(m)
                end
            end
        end
    end

    Test.@testset "Simplex" begin
        q  = rand(n)
        lb = rand(n)./n
        ub = 1 .+ rand(n)

        m = Projections.Simplex(q, lb, ub)
        test_model(m)
    end
end

max_stars_repo_name: VaclavMacha/Projections

max_issues_repo_name: VaclavMacha/Projections

max_forks_repo_name: VaclavMacha/Projections


text:
# ---
# title: 316. Remove Duplicate Letters
# id: problem316
# author: Indigo
# date: 2020-12-20
# difficulty: Medium
# categories: String, Stack, Greedy
# link: <https://leetcode.com/problems/remove-duplicate-letters/description/>
# hidden: true
# ---
#
# Given a string `s`, remove duplicate letters so that every letter appears once
# and only once. You must make sure your result is **the smallest in
# lexicographical order** among all possible results.
#
# **Note:** This question is the same as 1081:
# <https://leetcode.comhttps://leetcode.com/problems/smallest-subsequence-of-distinct-characters/>
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "bcabc"
#     Output: "abc"
#
#
# **Example 2:**
#
#
#
#     Input: s = "cbacdcbc"
#     Output: "acdb"
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 104`
#   * `s` consists of lowercase English letters.
#
#
## @lc code=start
using LeetCode

function remove_duplicate_letters(s::String)::String
    stk = Char[]
    added = Set{Char}()
    remain = fill(0, 26)
    for ch in s
        remain[ch - 'a' + 1] += 1
    end
    for ch in s
        if !(ch in added)
            while !isempty(stk) && stk[end] > ch && remain[stk[end] - 'a' + 1] > 0
                pop!(added, pop!(stk))
            end
            push!(added, ch)
            push!(stk, ch)
        end
        remain[ch - 'a' + 1] -= 1
    end
    return join(stk)
end
## @lc code=end


text:
"""
A simple common interface to diverse graphics packages backends useful for
drawings in SDD algorithms.
"""
module SDDGraphics

using Reexport

include("canvas.jl")

export
  canvassize,
  width,
  height,
  rectregion,
  xlims,
  ylims,
  xmin,
  ymin,
  xmax,
  ymax,
  square,
  insiderectregion,
  pointtopixel,
  pixeltopoint,
  pointtocanvas,
  canvastopoint

include("colors.jl")
include("coloringfunctions.jl")

export
  colortype,
  color,
  bgcolor,
  fgcolor,
  colorscheme,
  colormap,
  colorarray,
  updatecolorarray,
  colorinterpolation,
  colorinterpolationbg,
  colorinterpolationfg,
  coloringfunction,
  RadialColoringFunction, RadialCF,
  AngleColoringFunction, AngleCF

include("drawing.jl")

#ToDo!: Documentation.
"""
    configure(; kwargs...)

Shortcut function to configure graphics.

#### Arguments
- `canvassize=(w,h)` or `canvas=(w,h)`: Set the canvas size.
"""
  function configure(; kwargs...)
    for (k,v) in kwargs
      if k == :canvassize || k == :canvas
        canvassize(v...)
      elseif k == :width || k == :w
        width(v)
      elseif k == :height || k == :h
        height(v)
      elseif k == :xmin
        xmin(v)
      elseif k == :xmax
        xmax(v)
      elseif k == :ymin
        ymin(v)
      elseif k == :ymax
        ymax(v)
      elseif k == :xlims
        xlims(v...)
      elseif k == :ylims
        width(v...)
      elseif k == :rectregion || k == :rect
        rectregion(v...)
      elseif k == :color
        color(v)
      elseif k == :bgcolor || k == :bgc
        bgcolor(v)
      elseif k == :fgcolor || k == :fgc
        fgcolor(v)
      elseif k == :colormap || k == :cm || k == :colorscheme || k == :cs
        colorscheme(v)
      elseif k == :colormapinv || k == :cminv || k == :colorschemeinv || k == :csinv
        colorscheme(v,true)
      elseif k == :coloringfunction || k == :cf
        coloringfunction(v)
      elseif k == :axes
        axes(v)
      elseif k == :pointsize || k == :ps
        pointsize(v)
      elseif k == :linewidth || k == :lw
        linewidth(v)
      elseif k == :style
        style(v)
      #elseif k == :dash
      #  dash(v)
      end
    end
  end

export
  pointsize,
  linewidth,
  style, # :stroke, :fill, :fillstroke, :fillgrad
  #dash
  configure,
  drawingkind,
  newdrawing,
  drawpixel,
  drawpoint,
  drawline,
  #drawray,
  drawlinesegment,
  drawcircle,
  drawarc,
  #drawpath, #?
  #drawbox, #?
  #drawrect, #?
  axes,
  drawing

include("backends/images.jl") # Default backend
include("backends.jl")

export
  backend,
  backends,
  supported

end


text:
using MultiQuad, Test, Unitful

@testset "quad" begin
    rtol = 1e-10

    xmin = 0
    xmax = 4
    func1(x) = x^2 * exp(-x)
    res = 2 - 26 / exp(4)
    int, err = quad(func1, xmin, xmax, rtol = rtol)
    @test isapprox(int, res, atol = err)
    @test typeof(quad(
        func1,
        xmin,
        xmax,
        rtol = rtol,
        method = :vegas,
    )[1]) == Float64
    @test typeof(quad(
        func1,
        xmin,
        xmax,
        rtol = rtol,
        method = :suave,
    )[1]) == Float64

    xmin2 = 0 * u"m"
    xmax2 = 4 * u"m"
    func2(x) = x^2
    res2 = 64 / 3
    int2, err2 = quad(func2, xmin2, xmax2, rtol = rtol)
    @test isapprox(int2, res2 * u"m^3", atol = err2)
    @test typeof(quad(
        func2,
        xmin2,
        xmax2,
        rtol = rtol,
        method = :vegas,
    )[1]) == typeof(1.0 * u"m^3")
    @test typeof(quad(
        func2,
        xmin2,
        xmax2,
        rtol = rtol,
        method = :suave,
    )[1]) == typeof(1.0 * u"m^3")

    f(x) = x^4

    @test quad(f, -1, 1, method=:gausslegendre, order=100000)[1] ≈ 2/5

    @test quad(f, method=:gausshermite, order=10000)[1] ≈ 3(√π)/4

    @test quad(f, method=:gausslaguerre, order=10000)[1] ≈ 24
    @test quad(f, method=:gausslaguerre, order=3, α=1.0)[1] ≈ 120

    @test quad(f, -1, 1, method=:gausschebyshev, order=3)[1] ≈ 3π/8
    @test quad(f, -1, 1, method=:gausschebyshev, order=3, kind=1)[1] ≈ 3π/8
    @test quad(f, -1, 1, method=:gausschebyshev, order=3, kind=2)[1] ≈ π/16
    @test quad(f, -1, 1, method=:gausschebyshev, order=3, kind=3)[1] ≈ 3π/8
    @test quad(f, -1, 1, method=:gausschebyshev, order=3, kind=4)[1] ≈ 3π/8

    @test quad(f, -1, 1, method=:gaussradau, order=3)[1] ≈ 2/5
    @test quad(f, -1, 1, method=:gausslobatto, order=4)[1] ≈ 2/5

    @test_throws ErrorException quad(func1, 0, 4, method = :none)
end

@testset "dblquad" begin
    rtol = 1e-10

    xmin = 0
    xmax = 4
    ymin(x) = x^2
    ymax(x) = x^3

    func1(y, x) = y * x^3
    res1 = 241664 / 5
    int1, err1 = dblquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        rtol = rtol,
        method = :cuhre,
    )
    int1b, err1b = dblquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        rtol = rtol,
        method = :hcubature,
    )

    @test dblquad(
        (y,x)->x*y,
        0,
        2,
        0,
        3,
        rtol = rtol,
        method = :hcubature,
    )[1] ≈ 9

    @test isapprox(int1, res1, atol = err1)
    @test isapprox(int1b, res1, atol = err1b)
    @test typeof(dblquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        rtol = rtol,
        method = :vegas,
    )[1]) == Float64
    @test typeof(dblquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        rtol = rtol,
        method = :suave,
    )[1]) == Float64
    @test typeof(dblquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        rtol = rtol,
        method = :divonne,
    )[1]) == Float64

    xmin2 = 0 * u"m"
    xmax2 = 4 * u"m"
    ymin2(x) = x^2 * u"m"
    ymax2(x) = x^3
    func2(y, x) = y * x^3

    int2, err2 = dblquad(func2, xmin2, xmax2, ymin2, ymax2, rtol = rtol)
    int2b, err2b = dblquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        rtol = rtol,
        method = :hcubature,
    )

    @test isapprox(int2, res1 * u"m^10", atol = err2)
    @test isapprox(int2b, res1 * u"m^10", atol = err2b)
    @test typeof(dblquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        rtol = rtol,
        method = :vegas,
    )[1]) == typeof(1.0 * u"m^10")
    @test typeof(dblquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        rtol = rtol,
        method = :suave,
    )[1]) == typeof(1.0 * u"m^10")
    @test typeof(dblquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        rtol = rtol,
        method = :divonne,
    )[1]) == typeof(1.0 * u"m^10")

    @test_throws ErrorException dblquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        rtol = rtol,
        method = :none,
    )
end

@testset "tblquad" begin
    rtol = 1e-10

    xmin = 0
    xmax = 4
    ymin(x) = x^2
    ymax(x) = x^3
    zmin(x, y) = 3
    zmax(x, y) = 4

    func1(z, y, x) = y * x^3 * sin(z)
    res = 241664 / 5 * (cos(3) - cos(4))
    int1, err1 = tplquad(func1, xmin, xmax, ymin, ymax, zmin, zmax, rtol = rtol)
    int1b, err1b = tplquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        zmin,
        zmax,
        rtol = rtol,
        method = :hcubature,
    )
    @test isapprox(int1, res, atol = err1)
    @test isapprox(int1b, res, atol = err1b)
    @test typeof(tplquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        zmin,
        zmax,
        rtol = rtol,
        method = :vegas,
    )[1]) == Float64
    @test typeof(tplquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        zmin,
        zmax,
        rtol = rtol,
        method = :suave,
    )[1]) == Float64

    @test tplquad(
        (z,y,x)->x*y*z,
        0,
        2,
        0,
        2,
        0,
        3,
        rtol = rtol,
        method = :hcubature,
    )[1] ≈ 18

    @test typeof(tplquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        zmin,
        zmax,
        rtol = rtol,
        method = :divonne,
    )[1]) == Float64

    xmin2 = 0 * u"m"
    xmax2 = 4 * u"m"
    ymin2(x) = x^2 * u"m"
    ymax2(x) = x^3
    zmin2(x, y) = 3
    zmax2(x, y) = 4

    func2(z, y, x) = y * x^3 * sin(z)
    int2, err2 = tplquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        zmin2,
        zmax2,
        rtol = rtol,
    )
    int2b, err2b = tplquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        zmin2,
        zmax2,
        rtol = rtol,
        method = :hcubature,
    )
    @test isapprox(int2, res * u"m^10", atol = err2)
    @test isapprox(int2b, res * u"m^10", atol = err2b)
    @test typeof(tplquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        zmin2,
        zmax2,
        rtol = rtol,
        method = :vegas,
    )[1]) == typeof(1.0 * u"m^10")
    @test typeof(tplquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        zmin2,
        zmax2,
        rtol = rtol,
        method = :suave,
    )[1]) == typeof(1.0 * u"m^10")
    @test typeof(tplquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        zmin2,
        zmax2,
        rtol = rtol,
        method = :divonne,
    )[1]) == typeof(1.0 * u"m^10")

    @test_throws ErrorException tplquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        zmin2,
        zmax2,
        rtol = rtol,
        method = :none,
    )[1]
end


text:
module LED

import ..ioctl

export led_matrix, RGB565, led_clear, @colorant_str

using Colors, ColorTypes, FixedPointNumbers, Dates, Mmap

function _led_fb_device()
    try
        for devname in readdir("/sys/class/graphics")
            sysfname = joinpath("/sys/class/graphics",devname,"name")
            if startswith(devname, "fb") && isfile(sysfname)
                if startswith(readline(sysfname),"RPi-Sense FB")
                    return joinpath("/dev",devname)
                end
            end
        end
    catch e
        @warn("LED: Sense Hat not found.")
        return nothing
    end
end

const LED_FB_DEVICE_PATH = _led_fb_device()
const LED_FB_DEVICE = Ref{IOStream}()

function __init__()
    if LED_FB_DEVICE_PATH != nothing
        LED_FB_DEVICE[] = open(LED_FB_DEVICE_PATH,"w+")
    end
end

const U5 = Normed{UInt8,5}
const U6 = Normed{UInt8,6}

struct RGB565 <: AbstractRGB{Normed{UInt8, 8}}
    data::UInt16
end

function RGB565(r::U5, g::U6, b::U5)
    RGB565( (UInt16(reinterpret(r)) << 11) |
            (UInt16(reinterpret(g)) << 5) |
            (UInt16(reinterpret(b))) )
end

RGB565(r::Real, g::Real, b::Real) =
    RGB565(convert(U5,r), convert(U6,g), convert(U5,b))
RGB565(c::Union{Color,Colorant}) = RGB565(red(c), green(c), blue(c))

ColorTypes.red(c::RGB565) = U5(c.data >> 11, Val{true})
ColorTypes.green(c::RGB565) = U6((c.data >> 5) & 0x3f, Val{true})
ColorTypes.blue(c::RGB565) = U5(c.data & 0x1f, Val{true})

ColorTypes.ccolor(::Type{RGB565}, ::Type{Csrc}) where {Csrc <: Colorant} = RGB565
ColorTypes.base_color_type(::Type{RGB565}) = RGB565

"""
    led_matrix()

Returns an 8x8 matrix of `RGB565` elements that is memory-mapped to the Sense HAT LED Matrix.

While it is possible to invoke the function multiple times (each returning different
arrays), it is generally preferable to assign it once into a `const` variable so as to
minimise the number of open file handlers.

# Example

```
using SenseHat
using ColorTypes

const LED = led_matrix()

LED[:] = SenseHat.JULIA_LOGO
sleep(3)
LED[:] = RGB(0,0,0)
```

"""
led_matrix() = Mmap.mmap(LED_FB_DEVICE[], Array{RGB565,2}, (8,8); grow=false)

const FBIO_GET_GAMMA = 61696
const FBIO_SET_GAMMA = 61697
const FBIO_RESET_GAMMA = 61698

@enum(GammaScale, GAMMA_DEFAULT, GAMMA_LOW, GAMMA_USER)

"""
    getgamma()

Return the current gamma correction table for the Sense HAT LED matrix, as a vector of 32
`UFixed{UInt8,5}` numbers.

"""
function getgamma()
    buffer = zeros(UInt8, 32)
    ioctl(LED_FB_DEVICE[], FBIO_GET_GAMMA, buffer)
    return buffer
end

"""
    setgamma(v::GammaScale)
    setgamma(table::Vector)

Set the gamma correction for the Sense HAT LED matrix. This can either be one of the
predefined settings (`GAMMA_DEFAULT`, `GAMMA_LOW`, `GAMMA_USER`) or a table of values
between 0 and 1 (the actual values will be converted to `UFixed{UInt8, 5}`).

"""
function setgamma end

function setgamma(v::GammaScale)
    ioctl(LED_FB_DEVICE[], FBIO_RESET_GAMMA, Cint(v))
    return nothing
end

function setgamma(table::Vector{U5})
    length(table) == 32 || error("Table must contain 32 values.")
    ioctl(LED_FB_DEVICE[], FBIO_SET_GAMMA, table)
    return nothing
end
setgamma(table::AbstractVector) = setgamma(convert(Vector{U5}, table))

function led_display(X)
    Base.depwarn("`led_display(X)` has been deprecated, use `led_matrix()[:] = X` instead.", :led_display)
    led_matrix()[:] = X
end

"""
    led_clear()

Sets the SenseHat LED matrix to all black.
"""
function led_clear()
    led_matrix()[:] .= RGB565(0,0,0)
    return
end

end # module


text:
# mask.jl

using MIRT: embed, embed!, mask_or, mask_outline, maskit, getindex!

using Test: @test, @testset, @test_throws, @inferred
#using SparseArrays: sparse

function embed_test()
	mask = [false true true; true false false]
	a = [0 2 3; 1 0 0]; v = a[mask]
	b = similar(a)
	@test (@inferred embed(v,mask)) == a
	@test (@inferred embed([v 2v], mask)) == cat(dims=3, a, 2a)
	@test embed!(b, v, mask ; filler=-1) == a + (-1) * .!mask
#	@test embed(sparse(1:3),mask) == sparse([0 2 3; 1 0 0]) # later
	true
end

function getindex!_test()
	N = (2^2,2^3)
	mask = rand(N...) .< 0.5
	K = sum(mask)
	T = ComplexF32
	x = randn(T, N...)
	y0 = x[mask]
	y1 = Array{T}(undef, K)
	@test 0 == @allocated getindex!(y1, x, mask)
	@test y0 == y1
	true
end

function mask_or_test()
	mask2 = trues(3,4)
	@test (@inferred mask_or(mask2)) === mask2
	@test sum(mask_or(mask2)) == length(mask2)
	mask3 = trues(3,4,5)
	@test (@inferred mask_or(mask3)) == trues(3,4)
	@test_throws String mask_or(trues(1,))
	true
end

function mask_outline_test()
	mask2 = trues(3,4)
	@test (@inferred mask_outline(mask2)) == falses(3,4)
	mask3 = trues(3,4,5)
	@test (@inferred mask_outline(mask3)) == falses(3,4)
	true
end

function maskit_test()
	mask = [false true true; true false false]
	@inferred maskit([7 2 3; 1 7 7], mask)
	data = [7 2 3; 1 7 7]
	@test maskit(data, mask) == [1,2,3]
	data = repeat(data, 1, 1, 2)
	@test maskit(data, mask) == repeat(1:3, 1, 2)
	@test_throws DimensionMismatch maskit([0], mask)
	true
end

@testset "embed" begin
	@test embed_test()
end
@testset "getindex!" begin
	@test getindex!_test()
end
@testset "mask_or" begin
	@test mask_or_test()
end
@testset "mask_outline" begin
	@test mask_outline_test()
end
@testset "maskit" begin
	@test maskit_test()
end


text:
export ImPlotPoint, ImPlotRange, ImPlotLimits, ImPlotStyle, ImPlotInputMap, ImPlotContext
export CreateContext, DestroyContext, GetCurrentContext, SetCurrentContext, SetImGuiContext
export EndPlot
export SetNextPlotLimits, SetNextPlotLimitsX, SetNextPlotLimitsY, LinkNextPlotLimits, FitNextPlotAxes
export SetPlotYAxis, HideNextItem, IsPlotHovered, IsPlotXAxisHovered, IsPlotYAxisHoevered, IsPlotQueried
export DragLineX, DragLineY, DragPoint
export SetLegendLocation, SetMousePosLocation, IsLegendEntryHovered
export BeginLegendDragDropSource, EndLegendDragDropSource, BeginLegendPopup, EndLegendPopup
export PopStyleColor, PopStyleVar
export SetNextLineStyle, SetNextFillStyle, SetNextMarkerStyle, SetNextErrorBarStyle
export GetStyleColorName, GetMarkerName, PopColormap, GetColormapSize, ShowColormapScale
export GetColormapName, PushPlotClipRect, PopPlotClipRect
export ShowStyleSelector, ShowColormapSelector, ShowUserGuide, ShowDemoWindow


text:
export jacobi, jacobi!, gauss_seidel, gauss_seidel!, sor, sor!, ssor, ssor!

import LinearAlgebra.SingularException
import Base: getindex, iterate

function check_diag(A::AbstractMatrix)
    for i = 1 : size(A, 1)
        if iszero(A[i,i])
            throw(SingularException(i))
        end
    end
end

"""
    jacobi(A, b) -> x

Same as [`jacobi!`](@ref), but allocates a solution vector `x` initialized with zeros.
"""
jacobi(A::AbstractMatrix, b; kwargs...) = jacobi!(zerox(A, b), A, b; kwargs...)

"""
    jacobi!(x, A::AbstractMatrix, b; maxiter=10) -> x

Performs exactly `maxiter` Jacobi iterations.

Allocates a single temporary vector and traverses `A` columnwise.

Throws `LinearAlgebra.SingularException` when the diagonal has a zero. This check
is performed once beforehand.
"""
function jacobi!(x, A::AbstractMatrix, b; maxiter::Int=10)
    check_diag(A)
    iterable = DenseJacobiIterable(A, x, similar(x), b, maxiter)
    for _ = iterable end
    x
end

mutable struct DenseJacobiIterable{matT,vecT,solT,rhsT}
    A::matT
    x::solT
    next::vecT
    b::rhsT
    maxiter::Int
end

start(::DenseJacobiIterable) = 1
done(it::DenseJacobiIterable, iteration::Int) = iteration > it.maxiter
function iterate(j::DenseJacobiIterable, iteration::Int=start(j))
    if done(j, iteration) return nothing end

    n = size(j.A, 1)

    copyto!(j.next, j.b)

    # Computes next = b - (A - D)x
    for col = 1 : n
        @simd for row = 1 : col - 1
            @inbounds j.next[row] -= j.A[row, col] * j.x[col]
        end

        @simd for row = col + 1 : n
            @inbounds j.next[row] -= j.A[row, col] * j.x[col]
        end
    end

    # Computes x = D \ next
    for col = 1 : n
        @inbounds j.x[col] = j.next[col] / j.A[col, col]
    end

    nothing, iteration + 1
end

"""
    gauss_seidel(A, b) -> x

Same as [`gauss_seidel!`](@ref), but allocates a solution vector `x` initialized with zeros.
"""
gauss_seidel(A::AbstractMatrix, b; kwargs...) = gauss_seidel!(zerox(A, b), A, b; kwargs...)

"""
    gauss_seidel!(x, A::AbstractMatrix, b; maxiter=10) -> x

Performs exactly `maxiter` Gauss-Seidel iterations.

Works fully in-place and traverses `A` columnwise.

Throws `LinearAlgebra.SingularException` when the diagonal has a zero. This check
is performed once beforehand.
"""
function gauss_seidel!(x, A::AbstractMatrix, b; maxiter::Int=10)
    check_diag(A)
    iterable = DenseGaussSeidelIterable(A, x, b, maxiter)
    for _ = iterable end
    x
end

mutable struct DenseGaussSeidelIterable{matT,solT,rhsT}
    A::matT
    x::solT
    b::rhsT
    maxiter::Int
end

start(::DenseGaussSeidelIterable) = 1
done(it::DenseGaussSeidelIterable, iteration::Int) = iteration > it.maxiter

function iterate(s::DenseGaussSeidelIterable, iteration::Int=start(s))
    if done(s, iteration) return nothing end

    n = size(s.A, 1)

    for col = 1 : n
        @simd for row = 1 : col - 1
            @inbounds s.x[row] -= s.A[row, col] * s.x[col]
        end

        s.x[col] = s.b[col]
    end

    for col = 1 : n
        @inbounds s.x[col] /= s.A[col, col]
        @simd for row = col + 1 : n
            @inbounds s.x[row] -= s.A[row, col] * s.x[col]
        end
    end

    nothing, iteration + 1
end

"""
    sor(A, b, ω::Real) -> x

Same as [`sor!`](@ref), but allocates a solution vector `x` initialized with zeros.
"""
sor(A::AbstractMatrix, b, ω::Real; kwargs...) =
    sor!(zerox(A, b), A, b, ω; kwargs...)

"""
    sor!(x, A::AbstractMatrix, b, ω::Real; maxiter=10) -> x

Performs exactly `maxiter` SOR iterations with relaxation parameter `ω`.

Allocates a single temporary vector and traverses `A` columnwise.

Throws `LinearAlgebra.SingularException` when the diagonal has a zero. This check
is performed once beforehand.
"""
function sor!(x, A::AbstractMatrix, b, ω::Real; maxiter::Int=10)
    check_diag(A)
    iterable = DenseSORIterable(A, x, similar(x), b, ω, maxiter)
    for _ = iterable end
    x
end

mutable struct DenseSORIterable{matT,solT,vecT,rhsT,numT}
    A::matT
    x::solT
    tmp::vecT
    b::rhsT
    ω::numT
    maxiter::Int
end

start(::DenseSORIterable) = 1
done(it::DenseSORIterable, iteration::Int) = iteration > it.maxiter
function iterate(s::DenseSORIterable, iteration::Int=start(s))
    if done(s, iteration) return nothing end

    n = size(s.A, 1)

    for col = 1 : n
        @simd for row = 1 : col - 1
            @inbounds s.tmp[row] -= s.A[row, col] * s.x[col]
        end

        s.tmp[col] = s.b[col]
    end

    for col = 1 : n
        @inbounds s.x[col] += s.ω * (s.tmp[col] / s.A[col, col] - s.x[col])
        @simd for row = col + 1 : n
            @inbounds s.tmp[row] -= s.A[row, col] * s.x[col]
        end
    end

    nothing, iteration + 1
end

"""
    ssor(A, b, ω::Real) -> x

Same as [`ssor!`](@ref), but allocates a solution vector `x` initialized with zeros.
"""
ssor(A::AbstractMatrix, b, ω::Real; kwargs...) =
    ssor!(zerox(A, b), A, b, ω; kwargs...)

"""
    ssor!(x, A::AbstractMatrix, b, ω::Real; maxiter=10) -> x

Performs exactly `maxiter` SSOR iterations with relaxation parameter `ω`. Each iteration
is basically a forward *and* backward sweep of SOR.

Allocates a single temporary vector and traverses `A` columnwise.

Throws `LinearAlgebra.SingularException` when the diagonal has a zero. This check
is performed once beforehand.
"""
function ssor!(x, A::AbstractMatrix, b, ω::Real; maxiter::Int=10)
    check_diag(A)
    iterable = DenseSSORIterable(A, x, similar(x), b, ω, maxiter)
    for _ = iterable end
    x
end

mutable struct DenseSSORIterable{matT,solT,vecT,rhsT,numT}
    A::matT
    x::solT
    tmp::vecT
    b::rhsT
    ω::numT
    maxiter::Int
end

start(::DenseSSORIterable) = 1
done(it::DenseSSORIterable, iteration::Int) = iteration > it.maxiter
function iterate(s::DenseSSORIterable, iteration::Int=start(s))
    if done(s, iteration) return nothing end

    n = size(s.A, 1)

    for col = 1 : n
        @simd for row = 1 : col - 1
            @inbounds s.tmp[row] -= s.A[row, col] * s.x[col]
        end

        s.tmp[col] = s.b[col]
    end

    for col = 1 : n
        @inbounds s.x[col] += s.ω * (s.tmp[col] / s.A[col, col] - s.x[col])
        @simd for row = col + 1 : n
            @inbounds s.tmp[row] -= s.A[row, col] * s.x[col]
        end
    end

    for col = n : -1 : 1
        s.tmp[col] = s.b[col]
        @simd for row = col + 1 : n
            @inbounds s.tmp[row] -= s.A[row, col] * s.x[col]
        end
    end

    for col = n : -1 : 1
        @simd for row = 1 : col - 1
            @inbounds s.tmp[row] -= s.A[row, col] * s.x[col]
        end

        @inbounds s.x[col] += s.ω * (s.tmp[col] / s.A[col, col] - s.x[col])
    end

    nothing, iteration + 1
end


text:
# Simple examples to illustrate how to obtain a Bayes (junction) tree using
# with beautiful Tex labels in `.dot` and `.tex` format.

using Test

using IncrementalInference

# EXPERIMENTAL FEATURE, 4Q19: need `sudo apt install dot2tex`
import IncrementalInference: generateTexTree

@testset "testing generateTexTree" begin

# Create a dummy factor graph, with variables and constraints.
fg = initfg()

# Add four pose variables, with 'x' symbol.
addVariable!(fg, :x1, ContinuousScalar)
addVariable!(fg, :x2, ContinuousScalar)
addVariable!(fg, :x3, ContinuousScalar)
addVariable!(fg, :x4, ContinuousScalar)

# Add two landmark variables, with 'l' symbol.
addVariable!(fg, :l1, ContinuousScalar)
addVariable!(fg, :l2, ContinuousScalar)

# Add the pose chain constraints (odometry and priors).
addFactor!(fg, [:x1], Prior(Normal()))
addFactor!(fg, [:x1;:x2], LinearRelative(Normal()))
addFactor!(fg, [:x2;:x3], LinearRelative(Normal()))
addFactor!(fg, [:x3;:x4], LinearRelative(Normal()))

# Add the pose-landmark constraints (range measurements)
addFactor!(fg, [:x1;:l1], LinearRelative(Normal()))
addFactor!(fg, [:x2;:l1], LinearRelative(Normal()))
addFactor!(fg, [:x3;:l1], LinearRelative(Normal()))
addFactor!(fg, [:x2;:l2], LinearRelative(Normal()))
addFactor!(fg, [:x3;:l2], LinearRelative(Normal()))
addFactor!(fg, [:x4;:l2], LinearRelative(Normal()))

# Let's take a peek to see how our factor graph looks like.
# drawGraph(fg, show=true)
# As well as our tree (AMD ordering)
tree = buildTreeReset!(fg)
# drawTree(tree, show=true, imgs=false)

# Now, let's generate the corresponding `.dot` and `.tex`.
texTree = generateTexTree(tree, filepath=joinpath(@__DIR__,"tmp","bt"))

# All you have to do now is compile your newly created `.tex` file, probably
# include the `bm` package (`\usepackage{bm}`), and enjoy!

# fake a count of 1, since we are actually testing generateTexTree
@test true

end  # testset


text:
# This code generates the 2l+1 dimensional matrix Fourier coefficients
# for a function defined over the rotation and rotation/reflection groups
#
# The routines in this file are based upon:
#
#   FFTs on the Rotation Group
#   Peter J. Kostelec and Daniel N. Rockmore
#   J. Fourier. Anal. Appl., Vol. 14, Issue. 2, p. 145-179, 2008

import Base.fft

function w(B::Int, k::Int)
    s = 0.0

    for j in 0:B-1
        s += sin((2j + 1) * (2k + 1) * (π / (4B))) / (2j + 1)
    end

    s *= (2/B) * sin(π * (2k + 1) / (4B))
end

function S₁(f::Function, B::Int, k::Int, j₂::Int)
    β = π * (2k + 1) / (4B)
    γ = (2π * j₂) / (2B)

    A = [exp(-im * ((2π * j₁) / (2B)) * (B - 1)) * f(O3((2π * j₁) / (2B), β, γ)) for j₁ in 0:2B-1]
    A′ = ifft(A)
end

function S₂(f::Function, B::Int, k::Int)
    S = [S₁(f, B, k, j₂) for j₂ in 0:2B-1]
    A = Vector{Vector{Complex128}}(2B)
    A′ = Vector{Vector{Complex128}}(2B)

    for M′ in -(B-1):B
        A[M′ + B] = [exp(-im * ((2π * j₂) / (2B)) * (B - 1)) * S[j₂ + 1][M′ + B] for j₂ in 0:2B-1]
        A′[M′ + B] = ifft(A[M′ + B])
    end

    return hcat(A′...)
end

function fast_fourier_coeffs(f::Function, B::Int)
    S = [S₂(f, B, k) for k in 0:2B-1]
    f̂ = [zeros(Complex128, 2l + 1, 2l + 1) for l in 0:B-1]

    for l in 0:B-1
        for M in -l:l, M′ in -l:l
            for k in 0:2B-1
                β = π * (2k + 1) / (4B)
                f̂[l + 1][M′ + l + 1, M + l + 1] += wigner_d̃(l, M, M′, β) * w(B, k) * S[k + 1][M′ + B, M + B] / sqrt(2 / (2l + 1))
            end
        end
    end

    return f̂
end

fft(f::Function, ::Type{O3}, L::Int) = fast_fourier_coeffs(f, L)


text:
@testset "Common" begin
    @test get_code_center_frequency_ratio(GPSL1) ≈ 1/1540
end


text:
# norm, det, inv, eig, trace, dev
"""
    norm(::Vec)
    norm(::SecondOrderTensor)
    norm(::FourthOrderTensor)

Computes the norm of a tensor.

# Examples
```jldoctest
julia> A = rand(Tensor{2,3})
3×3 Tensor{2,3,Float64,9}:
 0.590845  0.460085  0.200586
 0.766797  0.794026  0.298614
 0.566237  0.854147  0.246837

julia> norm(A)
1.7377443667834924
```
"""
@inline LinearAlgebra.norm(v::Vec) = sqrt(dot(v, v))
@inline LinearAlgebra.norm(S::SecondOrderTensor) = sqrt(dcontract(S, S))

# special case for Tensor{4, 3} since it is faster than unrolling
@inline LinearAlgebra.norm(S::Tensor{4, 3}) = sqrt(mapreduce(abs2, +, S))

@generated function LinearAlgebra.norm(S::FourthOrderTensor{dim}) where {dim}
    idx(i,j,k,l) = compute_index(get_base(S), i, j, k, l)
    ex = Expr[]
    for l in 1:dim, k in 1:dim, j in 1:dim, i in 1:dim
        push!(ex, :(get_data(S)[$(idx(i,j,k,l))]))
    end
    exp = reducer(ex, ex)
    return quote
      $(Expr(:meta, :inline))
      @inbounds return sqrt($exp)
    end
end

"""
    det(::SecondOrderTensor)

Computes the determinant of a second order tensor.

# Examples
```jldoctest
julia> A = rand(SymmetricTensor{2,3})
3×3 SymmetricTensor{2,3,Float64,6}:
 0.590845  0.766797  0.566237
 0.766797  0.460085  0.794026
 0.566237  0.794026  0.854147

julia> det(A)
-0.1005427219925894
```
"""
@inline LinearAlgebra.det(t::SecondOrderTensor{1}) = @inbounds t[1,1]
@inline LinearAlgebra.det(t::SecondOrderTensor{2}) = @inbounds (t[1,1] * t[2,2] - t[1,2] * t[2,1])
@inline function LinearAlgebra.det(t::SecondOrderTensor{3})
    @inbounds (t[1,1] * (t[2,2]*t[3,3] - t[2,3]*t[3,2]) -
                  t[1,2] * (t[2,1]*t[3,3] - t[2,3]*t[3,1]) +
                  t[1,3] * (t[2,1]*t[3,2] - t[2,2]*t[3,1]))
end

"""
    inv(::SecondOrderTensor)

Computes the inverse of a second order tensor.

# Examples
```jldoctest
julia> A = rand(Tensor{2,3})
3×3 Tensor{2,3,Float64,9}:
 0.590845  0.460085  0.200586
 0.766797  0.794026  0.298614
 0.566237  0.854147  0.246837

julia> inv(A)
3×3 Tensor{2,3,Float64,9}:
  19.7146   -19.2802    7.30384
   6.73809  -10.7687    7.55198
 -68.541     81.4917  -38.8361
```
"""
@generated function Base.inv(t::Tensor{2, dim}) where {dim}
    Tt = get_base(t)
    idx(i,j) = compute_index(Tt, i, j)
    if dim == 1
        ex = :($Tt((dinv, )))
    elseif dim == 2
        ex = quote
            v = get_data(t)
            $Tt((v[$(idx(2,2))] * dinv, -v[$(idx(2,1))] * dinv,
                -v[$(idx(1,2))] * dinv,  v[$(idx(1,1))] * dinv))
        end
    else # dim == 3
        ex = quote
            v = get_data(t)
            $Tt(((v[$(idx(2,2))]*v[$(idx(3,3))] - v[$(idx(2,3))]*v[$(idx(3,2))]) * dinv,
                -(v[$(idx(2,1))]*v[$(idx(3,3))] - v[$(idx(2,3))]*v[$(idx(3,1))]) * dinv,
                 (v[$(idx(2,1))]*v[$(idx(3,2))] - v[$(idx(2,2))]*v[$(idx(3,1))]) * dinv,

                -(v[$(idx(1,2))]*v[$(idx(3,3))] - v[$(idx(1,3))]*v[$(idx(3,2))]) * dinv,
                 (v[$(idx(1,1))]*v[$(idx(3,3))] - v[$(idx(1,3))]*v[$(idx(3,1))]) * dinv,
                -(v[$(idx(1,1))]*v[$(idx(3,2))] - v[$(idx(1,2))]*v[$(idx(3,1))]) * dinv,

                 (v[$(idx(1,2))]*v[$(idx(2,3))] - v[$(idx(1,3))]*v[$(idx(2,2))]) * dinv,
                -(v[$(idx(1,1))]*v[$(idx(2,3))] - v[$(idx(1,3))]*v[$(idx(2,1))]) * dinv,
                 (v[$(idx(1,1))]*v[$(idx(2,2))] - v[$(idx(1,2))]*v[$(idx(2,1))]) * dinv))
        end
    end
    return quote
        $(Expr(:meta, :inline))
        dinv = 1 / det(t)
        @inbounds return $ex
    end
end

@generated function Base.inv(t::SymmetricTensor{2, dim}) where {dim}
    Tt = get_base(t)
    idx(i,j) = compute_index(Tt, i, j)
    if dim == 1
        ex = :($Tt((dinv, )))
    elseif dim == 2
        ex = quote
            v = get_data(t)
            $Tt((v[$(idx(2,2))] * dinv, -v[$(idx(2,1))] * dinv,
                 v[$(idx(1,1))] * dinv))
        end
    else # dim == 3
        ex = quote
            v = get_data(t)
            $Tt(((v[$(idx(2,2))]*v[$(idx(3,3))] - v[$(idx(2,3))]*v[$(idx(3,2))]) * dinv,
                -(v[$(idx(2,1))]*v[$(idx(3,3))] - v[$(idx(2,3))]*v[$(idx(3,1))]) * dinv,
                 (v[$(idx(2,1))]*v[$(idx(3,2))] - v[$(idx(2,2))]*v[$(idx(3,1))]) * dinv,

                 (v[$(idx(1,1))]*v[$(idx(3,3))] - v[$(idx(1,3))]*v[$(idx(3,1))]) * dinv,
                -(v[$(idx(1,1))]*v[$(idx(3,2))] - v[$(idx(1,2))]*v[$(idx(3,1))]) * dinv,

                 (v[$(idx(1,1))]*v[$(idx(2,2))] - v[$(idx(1,2))]*v[$(idx(2,1))]) * dinv))
        end
    end
    return quote
        $(Expr(:meta, :inline))
        dinv = 1 / det(t)
        @inbounds return $ex
    end
end

function Base.inv(t::Tensor{4, dim}) where {dim}
    fromvoigt(Tensor{4, dim}, inv(tovoigt(t)))
end

function Base.inv(t::SymmetricTensor{4, dim, T}) where {dim, T}
    frommandel(SymmetricTensor{4, dim}, inv(tomandel(t)))
end

Base.:\(S1::SecondOrderTensor, S2::AbstractTensor) = inv(S1) ⋅ S2

"""
    eigvals(::SymmetricTensor)

Compute the eigenvalues of a symmetric tensor.
"""
@inline LinearAlgebra.eigvals(S::SymmetricTensor) = eigvals(eigen(S))

"""
    eigvecs(::SymmetricTensor)

Compute the eigenvectors of a symmetric tensor.
"""
@inline LinearAlgebra.eigvecs(S::SymmetricTensor) = eigvecs(eigen(S))

struct Eigen{T, dim, M}
    values::Vec{dim, T}
    vectors::Tensor{2, dim, T, M}
end

struct FourthOrderEigen{dim,T,M}
    values::Vector{T}
    vectors::Vector{SymmetricTensor{2,dim,T,M}}
end

# destructure via iteration
function Base.iterate(E::Union{Eigen,FourthOrderEigen}, state::Int=1)
    return iterate((eigvals(E), eigvecs(E)), state)
end

"""
    eigen(A::SymmetricTensor{2})

Compute the eigenvalues and eigenvectors of a symmetric second order tensor
and return an `Eigen` object. The eigenvalues are stored in a `Vec`,
sorted in ascending order. The corresponding eigenvectors are stored
as the columns of a `Tensor`.

See [`eigvals`](@ref) and [`eigvecs`](@ref).

# Examples
```jldoctest
julia> A = rand(SymmetricTensor{2, 2});

julia> E = eigen(A);

julia> E.values
2-element Tensor{1,2,Float64,2}:
 -0.1883547111127678
  1.345436766284664

julia> E.vectors
2×2 Tensor{2,2,Float64,4}:
 -0.701412  0.712756
  0.712756  0.701412
```
"""
LinearAlgebra.eigen(::SymmetricTensor{2})

"""
    eigvals(::Union{Eigen,FourthOrderEigen})

Extract eigenvalues from an `Eigen` or `FourthOrderEigen` object,
returned by [`eigen`](@ref).
"""
@inline LinearAlgebra.eigvals(E::Union{Eigen,FourthOrderEigen}) = E.values
"""
    eigvecs(::Union{Eigen,FourthOrderEigen})

Extract eigenvectors from an `Eigen` or `FourthOrderEigen` object,
returned by [`eigen`](@ref).
"""
@inline LinearAlgebra.eigvecs(E::Union{Eigen,FourthOrderEigen}) = E.vectors

"""
    eigen(A::SymmetricTensor{4})

Compute the eigenvalues and second order eigentensors of a symmetric fourth
order tensor and return an `FourthOrderEigen` object. The eigenvalues and
eigentensors are sorted in ascending order of the eigenvalues.

See also [`eigvals`](@ref) and [`eigvecs`](@ref).
"""
function LinearAlgebra.eigen(S::SymmetricTensor{4,dim,T}) where {dim,T}
    # Scale from the right only; tovoigt(S; offdiagscale=2) scales from left and right.
    S′ = SymmetricTensor{4,dim,T}((i,j,k,l) -> k == l ? S[i,j,k,l] : 2*S[i,j,k,l])
    v = tovoigt(S′)
    E = eigen(v)
    perm = sortperm(E.values)
    values = T[E.values[i] for i in perm]
    vectors = [fromvoigt(SymmetricTensor{2,dim,T}, view(E.vectors, :, i)) for i in perm]
    return FourthOrderEigen(values, vectors)
end

"""
    sqrt(S::SymmetricTensor{2})

Calculate the square root of the positive definite symmetric
second order tensor `S`, such that `√S ⋅ √S == S`.

# Examples
```jldoctest
julia> S = rand(SymmetricTensor{2,2}); S = tdot(S)
2×2 SymmetricTensor{2,2,Float64,3}:
 0.937075  0.887247
 0.887247  0.908603

julia> sqrt(S)
2×2 SymmetricTensor{2,2,Float64,3}:
 0.776178  0.578467
 0.578467  0.757614

julia> √S ⋅ √S ≈ S
true
```
"""
Base.sqrt(::SymmetricTensor{2})

Base.sqrt(S::SymmetricTensor{2,1}) = SymmetricTensor{2,1}((sqrt(S[1,1]),))

# https://en.m.wikipedia.org/wiki/Square_root_of_a_2_by_2_matrix
function Base.sqrt(S::SymmetricTensor{2,2})
    s = √(det(S))
    t = √(tr(S)+2s)
    return SymmetricTensor{2,2}((S[1,1]+s, S[2,1], S[2,2]+s)) / t
end

function Base.sqrt(S::SymmetricTensor{2,3,T}) where T
    E = eigen(S)
    λ = E.values
    Φ = E.vectors
    z = zero(T)
    Λ = Tensor{2,3}((√(λ[1]), z, z, z, √(λ[2]), z, z, z, √(λ[3])))
    return symmetric(Φ⋅Λ⋅Φ')
end

"""
    tr(::SecondOrderTensor)

Computes the trace of a second order tensor.

# Examples
```jldoctest
julia> A = rand(SymmetricTensor{2,3})
3×3 SymmetricTensor{2,3,Float64,6}:
 0.590845  0.766797  0.566237
 0.766797  0.460085  0.794026
 0.566237  0.794026  0.854147

julia> tr(A)
1.9050765715072775
```
"""
@generated function LinearAlgebra.tr(S::SecondOrderTensor{dim}) where {dim}
    idx(i,j) = compute_index(get_base(S), i, j)
    ex = Expr[:(get_data(S)[$(idx(i,i))]) for i in 1:dim]
    exp = reduce((ex1, ex2) -> :(+($ex1, $ex2)), ex)
    @inbounds return exp
end

Statistics.mean(S::SecondOrderTensor) = tr(S) / 3

"""
    vol(::SecondOrderTensor)

Computes the volumetric part of a second order tensor
based on the additive decomposition.

# Examples
```jldoctest
julia> A = rand(SymmetricTensor{2,3})
3×3 SymmetricTensor{2,3,Float64,6}:
 0.590845  0.766797  0.566237
 0.766797  0.460085  0.794026
 0.566237  0.794026  0.854147

julia> vol(A)
3×3 SymmetricTensor{2,3,Float64,6}:
 0.635026  0.0       0.0
 0.0       0.635026  0.0
 0.0       0.0       0.635026

julia> vol(A) + dev(A) ≈ A
true
```
"""
vol(S::SecondOrderTensor) = mean(S) * one(S)

"""
    dev(::SecondOrderTensor)

Computes the deviatoric part of a second order tensor.

# Examples
```jldoctest
julia> A = rand(Tensor{2, 3});

julia> dev(A)
3×3 Tensor{2,3,Float64,9}:
 0.0469421  0.460085   0.200586
 0.766797   0.250123   0.298614
 0.566237   0.854147  -0.297065

julia> tr(dev(A))
0.0
```
"""
@inline function dev(S::SecondOrderTensor)
    Tt = get_base(typeof(S))
    trace = tr(S) / 3
    Tt(
        @inline function(i, j)
            @inbounds  v = i == j ? S[i,j] - trace : S[i,j]
            v
        end
    )
end

# http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/
"""
    rotate(x::Vec{3}, u::Vec{3}, θ::Number)

Rotate a three dimensional vector `x` around another vector `u` a total of `θ` radians.

# Examples
```jldoctest
julia> x = Vec{3}((0.0, 0.0, 1.0))
3-element Tensor{1,3,Float64,3}:
 0.0
 0.0
 1.0

julia> u = Vec{3}((0.0, 1.0, 0.0))
3-element Tensor{1,3,Float64,3}:
 0.0
 1.0
 0.0

julia> rotate(x, u, π/2)
3-element Tensor{1,3,Float64,3}:
 1.0
 0.0
 6.123233995736766e-17
```
"""
function rotate(x::Vec{3}, u::Vec{3}, θ::Number)
    ux = u ⋅ x
    u² = u ⋅ u
    c = cos(θ)
    s = sin(θ)
    (u * ux * (1 - c) + u² * x * c + sqrt(u²) * (u × x) * s) / u²
end


text:
@testset "serve" begin
    # q: request
    # p: reply
    # |: batch_wait_timeout_s reached

    @testset "batch_wait_timeout_s = 0" begin
        # case 1: qpqpqp
        struct DummyModel end

        function OL.handle(m::DummyModel, msg::OL.RequestMsg)
            xs = msg.msg
            @debug "value received in model" xs typeof(xs)
            sleep(0.1)
            res = vec(sum(xs; dims=1)) .+ 1
            OL.async_rep(msg.from, res)
        end
        model = @actor DummyModel() name="model"

        server = @actor OL.BatchStrategy(model;max_batch_size=2) name="server"

        worker = OL.Mailbox()
        t = @elapsed for i in 1:10
            put!(server,OL.RequestMsg([i], worker))
            sleep(0.11)
        end

        for i in 1:10
            msg = take!(worker)
            @test msg.msg == i+1
        end

        # case 2: qqqqqpqqqqp
        t = @elapsed begin
            for i in 1:5
                put!(server,OL.RequestMsg([i], worker))
            end
            for i in 1:5
                @test take!(worker).msg == i+1
            end
        end

        t = @elapsed begin
            for i in 1:64
                put!(server,OL.RequestMsg([i], worker))
            end
            for i in 1:64
                @test take!(worker).msg == i+1
            end
        end
    end

end


text:
ALL_TASK_OFFSETS[("linux-x86", v"1.5.4")] = Dict(
    :END => 360,
    :backtrace => 28,
    :bufsz => 204,
    :copy_stack => 208,
    :ctx => 44,
    :donenotify => 16,
    :eh => 212,
    :exception => 24,
    :excstack => 220,
    :gcstack => 216,
    :locks => 228,
    :logstate => 32,
    :next => 0,
    :prio => 226,
    :queue => 4,
    :result => 20,
    :sizeof_ctx => 156,
    :start => 36,
    :state => 12,
    :sticky => 40,
    :stkbuf => 200,
    :tid => 224,
    :timing_stack => 356,
    :tls => 8,
    :tls_base_context => 3408,
)


text:
# For proximal methods: https://github.com/kul-forbes/ProximalOperators.jl

module NonsmoothOptim

using LinearAlgebra

import JuMP

export OptimisationSense, Minimise, Maximise,
       UnconstrainedNonSmoothProblem, ConstrainedNonSmoothProblem, ProjectedConstrainedNonSmoothProblem,
       NonSmoothSolver, UnconstrainedNonSmoothSolver, ConstrainedNonSmoothSolver, ProjectedConstrainedNonSmoothSolver,
       solve,
       StepSizeRule, step, ConstantStepSize, ConstantStepLength, InverseStepSize,
       SubgradientMethod, ProjectedSubgradientMethod, BundleMethod

## Define an optimisation problem.

@enum OptimisationSense Minimise Maximise

"""
    abstract type NonSmoothProblem

Represents a nonsmooth optimisation problem. Any implementation must provide the following properties:

  * `f`: the function to optimise. It takes a a vector of `dimension` Float64 and generates a single Float64.
  * `g`: a function that returns any subgradient of `f`. It takes a vector of `dimension` Float64 and generates
    a vector of `dimension` Float64.
  * `dimension`: the dimension of the problem, i.e. the number of variables (an `Integer`).
  * `sense`: the `OptimisationSense` when processing `f` (i.e. minimise or maximise).

No automatic differentiation is performed in this package. Indeed, these techniques are usually not ensured to be
correct for nonsmooth functions. https://arxiv.org/pdf/1809.08530.pdf
"""
abstract type NonSmoothProblem end

"""
    struct UnconstrainedNonSmoothProblem <: NonSmoothProblem

Represents an unconstrained nonsmooth optimisation problem. This type provides a barebones implementation of the
`NonSmoothProblem` interface.
"""
struct UnconstrainedNonSmoothProblem <: NonSmoothProblem
  f::Function
  g::Function
  dimension::Int
  sense::OptimisationSense
end

"""
    abstract type ConstrainedNonSmoothProblem <: NonSmoothProblem

Represents a constrained nonsmooth optimisation problem. It does not imply any new properties with respect to a
`NonSmoothProblem`, because they might heavily depend on the kind of optimisation solver that is used
(indicator function, polytope, projection operator, proximale operator, etc.)
"""
abstract type ConstrainedNonSmoothProblem <: NonSmoothProblem end

"""
    struct ProjectedConstrainedNonSmoothProblem

Represents a constrained nonsmooth optimisation problem whose constraints are imposed through a projection operator.
This structure provides one more member to the `NonSmoothProblem` interface:

  * `project`: the projection operator. It takes a a vector of `dimension` Float64 and generates a vector of
    `dimension` Float64. The returned point must be feasible according to the definition of the optimisation problem.

TODO: a stricter interface for projection operators? I.e. closed-form expressions or iterative solvers (they need parameters; maybe use a smooth optimisation package for this, with just constraints to be added, like `NLPModels`?).
"""
struct ProjectedConstrainedNonSmoothProblem <: ConstrainedNonSmoothProblem
  f::Function
  g::Function
  project::Function
  dimension::Int
  sense::OptimisationSense
end

## Interface for a solver.

abstract type NonSmoothSolver end # Also contains the solver parameters, like tolerance, number of iterations, etc.
abstract type UnconstrainedNonSmoothSolver <: NonSmoothSolver end
abstract type ConstrainedNonSmoothSolver <: NonSmoothSolver end
abstract type ProjectedConstrainedNonSmoothSolver <: ConstrainedNonSmoothSolver end

function solve
  # Input:
  #   - A problem to solve.
  #   - An NSO method.
  #   - An initial iterate.
  # Keyword arguments:
  #   - info_callback: (k, f, x, g, f_best, x_best, t_iter) -> nothing
  # Output: best solution found (Vector{Float64})
end

## Rules for step sizes, likely to be used by several solvers.

abstract type StepSizeRule end

function step
  # Input:
  #   - Rule to apply.
  #   - Number of iterations.
  #   - Norm of the subgradient.
  # Output: step (Float64)
end

struct ConstantStepSize <: StepSizeRule
  step::Float64
end

function step(rule::ConstantStepSize, ::Int, ::Float64)
  return rule.step
end

struct ConstantStepLength <: StepSizeRule
  step::Float64
end

function step(rule::ConstantStepLength, ::Int, sg_norm::Float64)
  return rule.step / sg_norm
end

struct InverseStepSize <: StepSizeRule # Square summable, but not summable
  step::Float64
end

function step(rule::InverseStepSize, n_iter::Int, ::Float64)
  if n_iter <= 0
    return rule.step
  else
    return rule.step / n_iter
  end
end

# Actual code.
include("subgradient.jl")
include("bundle.jl")

end


text:
module RelocatableFolders

import Scratch, SHA

export @path

function safe_isfile(file)
    return try
        isfile(file)
    catch err
        err isa Base.IOError || rethrow()
        false
    end
end

function safe_ispath(file)
    return try
        ispath(file)
    catch err
        err isa Base.IOError || rethrow()
        false
    end
end

"""
    @path expr

Define a "relocatable" path. Its contents will be available regardless of
whether the original path still exists or not. The contents of the path is
stored within the returned `Path` object and the folder structure is recreated
as a scratchspace if the original does not exist anymore. Calling any path
manipulation functions, such as `joinpath`, on the `Path` will return a path
to a valid folder.
"""
macro path(expr)
    file = string(__source__.file)
    dir = safe_isfile(file) ? dirname(file) : pwd()
    return :($(Path)($__module__, $dir, $(esc(expr))))
end

struct Path <: AbstractString
    is_dir::Bool
    mod::Module
    path::String
    hash::String
    files::Dict{String,Vector{UInt8}}

    function Path(mod::Module, dir, path::AbstractString)
        path = isabspath(path) ? path : joinpath(dir, path)
        path = normpath(path)
        safe_ispath(path) || throw(ArgumentError("not a path: `$path`"))
        is_dir = isdir(path)
        dir = is_dir ? path : dirname(path)
        files = Dict{String,Vector{UInt8}}()
        ctx = SHA.SHA1_CTX()
        for (root, _, fs) in walkdir(dir), f in fs
            fullpath = joinpath(root, f)
            if is_dir || path == fullpath
                include_dependency(fullpath)
                SHA.update!(ctx, codeunits(fullpath))
                content = read(fullpath)
                SHA.update!(ctx, content)
                files[relpath(fullpath, dir)] = content
            end
        end
        return new(is_dir, mod, dir, string(Base.SHA1(SHA.digest!(ctx))), files)
    end
end

Base.show(io::IO, path::Path) = print(io, repr(getroot(path)))
Base.ncodeunits(f::Path) = ncodeunits(getpath(f))
Base.isvalid(f::Path, index::Integer) = isvalid(getpath(f), index)
Base.iterate(f::Path) = iterate(getpath(f))
Base.iterate(f::Path, state::Integer) = iterate(getpath(f), state)
Base.String(f::Path) = String(getpath(f))

function getpath(f::Path)
    safe_ispath(f.path) && return getroot(f)
    dir = Scratch.get_scratch!(f.mod, f.hash)
    if !isempty(f.files) && !safe_ispath(joinpath(dir, first(keys(f.files))))
        cd(dir) do
            for (file, blob) in f.files
                mkpath(dirname(file))
                write(file, blob)
            end
        end
    end
    return getroot(f, dir)
end

getroot(p::Path, root = p.path) = p.is_dir ? root : joinpath(root, first(keys(p.files)))

end # module


text:
include("dependencies_for_runtests.jl")

include("regression_tests/shallow_water_bickley_jet_regression.jl")

@testset "Shallow Water Regression" begin
    @info "Running shallow water regression tests..."

    for arch in archs
        for formulation in (VectorInvariantFormulation(), ConservativeFormulation())
            @testset "Shallow Water Bickley jet simulation [$(typeof(arch)), $(typeof(formulation))]" begin
                @info "  Testing shallow water Bickley jet simulation regression [$(typeof(arch)), $(typeof(formulation))]"
                run_shallow_water_regression(arch, formulation; regenerate_data = false)
            end
        end
    end
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 623b2240-c520-11eb-0b1e-4b03cb67ac24
begin
	using PlutoUI,Plots, ForwardDiff
	plotlyjs()
end

# ╔═╡ c4ec78ff-b9c7-4867-81df-9384b965881e
md"## Polinomio de Taylor"

# ╔═╡ bc296b17-b6ca-476d-b9f4-1743559f0a6d
f(x) = sin(x)

# p₀(x=0) = f(x=0) = cos(0) = 1
# f'(x) = -sen(x) // f'(x=0) = -sen(0) = 0

# ╔═╡ 916c7595-d8ee-4ce1-9083-f403613131f7
md"""
a₀ = $(@bind a₀ Slider(-2:0.1:2,default=1,show_value=true))

a₁ = $(@bind a₁ Slider(-2:0.1:2,default=0,show_value=true))

a₂ = $(@bind a₂ Slider(-2:0.1:2,default=-0.5,show_value=true))

a₃ = $(@bind a₃ Slider(-2:0.1:2,default=0,show_value=true))

a₄ = $(@bind a₄ Slider(-2:0.1:2,default=0,show_value=true))
"""

# ╔═╡ 0831c272-bd4a-459a-9d8b-75bc4f926022
md"""
p₀(x) = $(a₀)

p₁(x) = $(a₀) + $(a₁) x

p₂(x) = $(a₀) + $(a₁) x + $(a₂) x²

p₃(x) = $(a₀) + $(a₁) x + $(a₂) x² + $(a₃) x³

p₄(x) = $(a₀) + $(a₁) x + $(a₂) x² + $(a₃) x³ + + $(a₄) x⁴
"""

# ╔═╡ 21815a23-2d33-4813-baca-21adcfa4a85f
begin
	p₀(x) = a₀
	p₁(x) = a₀ + a₁ * x
	p₂(x) = a₀ + a₁ * x + a₂ * x^2
	p₃(x) = a₀ + a₁ * x + a₂ * x^2 +  a₃ * x^3
	p₄(x) = a₀ + a₁ * x + a₂ * x^2 +  a₃ * x^3 + a₄ * x^4
	# p₂(x) = a₀ + a₁ * x + a₂ * x^2
end

# ╔═╡ 9f161f71-89db-44d3-84f3-be08176ff802
begin
	plot(f,xlim=[-2,2],label=false,width=2)
	plot!(p₀,label="p₀",width=2)
	plot!(p₁,label="p₁",width=2)
	plot!(p₂,label="p₂",width=2)
	plot!(p₃,label="p₃",width=2)
	plot!(p₄,label="p₄",width=2)
end

# ╔═╡ Cell order:
# ╟─623b2240-c520-11eb-0b1e-4b03cb67ac24
# ╟─c4ec78ff-b9c7-4867-81df-9384b965881e
# ╠═bc296b17-b6ca-476d-b9f4-1743559f0a6d
# ╟─916c7595-d8ee-4ce1-9083-f403613131f7
# ╟─0831c272-bd4a-459a-9d8b-75bc4f926022
# ╟─21815a23-2d33-4813-baca-21adcfa4a85f
# ╟─9f161f71-89db-44d3-84f3-be08176ff802

max_stars_repo_name: akielbowicz/cbc-analisis-ingenieria

max_issues_repo_name: akielbowicz/cbc-analisis-ingenieria

max_forks_repo_name: akielbowicz/cbc-analisis-ingenieria


text:
"""
Calculate both the spatial and the operational entanglement entropies of a
region A, using the SVD. The latter is the "entanglement of particles"
introduced by Wiseman and Vaccaro in 2003.
"""
function operational_entanglement(L::Int, N::Int, l::Int, precision::Int, alpha::Float64, nu::Vector{Float64}, eigenvalues::Bool, probabilities::Bool)

    setprecision(precision)
#set_bigfloat_precision(precision)
# =======================================================
#Calculating the probability distribution (PN) of the number of particles in the targeted partition.
    scaling_factor= BigFloat
    scaling_factor=10.0^300
    normsDim=l+1
    #P1 = zeros(BigFloat, 2)
    PN_old = zeros(BigFloat, normsDim)
    PN_new = zeros(BigFloat, normsDim)
    if nu[1]>=0 && nu[1]<=1
       nuVal=nu[1]
    elseif nu[1]>1
       nuVal=1.0
    else
       nuVal=0.0
    end
    PN_old[1]=scaling_factor*(1-nuVal)
    PN_old[2]=scaling_factor*nuVal
    NN=2
    for i= 2:normsDim-1
       if nu[i]>=0 && nu[i]<=1
          nuVal=nu[i]
       elseif nu[i]>1
          nuVal=1.0
       else
          nuVal=0.0
       end
       PN_new[1]=(1-nuVal)*PN_old[1]
       PN_new[NN+1]=nuVal*PN_old[NN]
       for k= 1:NN-1
          PN_new[k+1]=(1-nuVal)*PN_old[k+1]+nuVal*PN_old[k]
       end
       NN=NN+1
       for k=1: NN
          PN_old[k]=PN_new[k]
       end
    end
    if l==1
       for k=1: normsDim
          PN_new[k]=PN_old[k]
       end
    end
    PN_new/= scaling_factor
    err_PN = abs(sum(PN_new) - 1.0)
    if err_PN > 1e-12
        @warn("probability distribution error: $(err_PN)")

    end
# =======================================================
# =======================================================
#Calculating the second Renyi entanglement entropy.
    Salpha = zeros(BigFloat, 3)
    Salpha_old = zeros(BigFloat, normsDim)
    Salpha_new = zeros(BigFloat, normsDim)
    #alpha=2.0
    if nu[1]>=0 && nu[1]<=1
       nuVal=nu[1]
    elseif nu[1]>1
       nuVal=1.0
    else
       nuVal=0.0
    end
    Salpha_old[1]=scaling_factor*(1-nuVal)^alpha
    Salpha_old[2]=scaling_factor*nuVal^alpha
    NN=2
    for i= 2:normsDim-1
       if nu[i]>=0 && nu[i]<=1
          nuVal=nu[i]
       elseif nu[i]>1
          nuVal=1.0
       else
          nuVal=0.0
       end
       Salpha_new[1]=(1-nuVal)^alpha*Salpha_old[1]
       Salpha_new[NN+1]=nuVal^alpha*Salpha_old[NN]
       for k= 1:NN-1
          Salpha_new[k+1]=(1-nuVal)^alpha*Salpha_old[k+1]+nuVal^alpha*Salpha_old[k]
       end
       NN=NN+1
       for k=1: NN
          Salpha_old[k]=Salpha_new[k]
       end
    end

    if l==1
       for k=1: normsDim
          Salpha_new[k]=Salpha_old[k]
       end
    end
    Salpha_new/= scaling_factor
# =======================================================
    LogSumPNalpha= BigFloat
    SalphalogSalpha= BigFloat
    PnlogPN= BigFloat
    sigma2_n= BigFloat
    sigma2_n_PN= BigFloat
    n= BigFloat
    n2= BigFloat

    H4= BigFloat
    H5= BigFloat
    H6= BigFloat

    S1_sp= BigFloat
    Salpha_sp= BigFloat
    S1_op= BigFloat
    Salpha_op= BigFloat
# =======================================================

    S1_sp=0.0
    Salpha_sp=0.0
    for i= 1:normsDim-1
       if nu[i]>0 && nu[i]<1
             S1_sp-=(nu[i]*log(nu[i])+(1-nu[i])*log((1-nu[i])))
             Salpha_sp-=log(nu[i]^alpha+(1-nu[i])^alpha)
       end
    end
    err_Salpha = abs(sum(Salpha_new) - exp(-1.0*Salpha_sp))/exp(-1.0*Salpha_sp)
    if err_Salpha > 1e-12
        @warn("Salpha_op error: $(err_Salpha)")

    end
    Salpha_sp/=(alpha-1)
# =======================================================

    sigma2_n=0.0
    for i= 1:normsDim-1
       if nu[i]>0 && nu[i]<1
          sigma2_n+=(nu[i]*(1-nu[i]))
       end
    end

    sigma2_n_PN=0.0
    n=0
    n2=0
    for i= 1: normsDim
       if PN_new[i]>0
          n+=PN_new[i]*(i-1)
          n2+=PN_new[i]*(i-1)^2
       end
    end
    sigma2_n_PN=n2-n^2

    err_sigma2 = abs(sigma2_n_PN - sigma2_n)/sigma2_n
    if err_sigma2 > 1e-12
        @warn("sigma2_n error: $(err_sigma2)")

    end
# =======================================================

    PnlogPN=0.0
    for i= 1: normsDim
       if PN_new[i]>0
          PnlogPN+=PN_new[i]*log(PN_new[i])
       end
    end
# =======================================================

    S1_op=0.0
    S1_op=S1_sp+PnlogPN

    PNlogSalpha=0.0
    for i= 1: normsDim
       if Salpha_new[i]>0
          PNlogSalpha-=PN_new[i]*log(Salpha_new[i])
       end
    end

    Salpha_op= PNlogSalpha+alpha*PnlogPN
    Salpha_op/=(alpha-1)

    LogSumPNalpha=sum(PN_new.^alpha)
    LogSumPNalpha=-log(LogSumPNalpha)/(alpha-1)

# =======================================================

    H4=0.0
    H5=0.0
    H6=0.0
    H = zeros(BigFloat, normsDim)
    for i= 1: normsDim
       if PN_new[i]>0
          H[i]+=Salpha_new[i]/PN_new[i]^alpha
          H4+=PN_new[i]*H[i]
          H5+=PN_new[i]*H[i]^(1/(alpha-1))
          H6+=Salpha_new[i]^(1/alpha)
       end

    end

    H4=log(H4)/(1-alpha)
    H5=-log(H5)
    H6=-log(H6)*(alpha/(alpha-1))

# =======================New=============================
# =======================================================

     if l==100
        if probabilities
           Output="L$(L)N$(N)l$(l)alpha$(alpha)Pn.dat"
           open(Output, "w") do fp
              write(fp, "# L=$(L), N=$(N), l=$(l), alpha=$(alpha) \n")
              write(fp, "# ==============================================================================\n")
              write(fp, @sprintf "#%7s%24s%24s%24s\n" "n"  "P(n)" "P(n)^alpha"  "P(n)_alpha" )

              sumPnalpha=0.0
              sumPn_alpha=0.0
              sumSalpha_new =sum(Salpha_new)
              sumPN_newalpha =sum(PN_new.^alpha)
              for i=1: normsDim
                 write(fp, @sprintf "%8s%24.12E%24.12E%24.12E\n" "$(i-1)"  PN_new[i]  PN_new[i]^alpha/sumPN_newalpha Salpha_new[i]/sumSalpha_new)
                 sumPnalpha= sumPnalpha+PN_new[i]^alpha/sumPN_newalpha
                 sumPn_alpha= sumPn_alpha+Salpha_new[i]/sumSalpha_new
              end
               write(fp, "#$(sumPnalpha)  \n")
               write(fp, "#$(sumPn_alpha)  \n")
               flush(fp)
           end
        end

#___________________________________________eigvals_______________________________________________
        if eigenvalues
           Outputn="L$(L)N$(N)l$(l)alpha$(alpha)eigvals.dat"
           open(Outputn, "w") do fn
              write(fn, "# L=$(L), N=$(N), l=$(l) \n")
              write(fn, "# ===============================\n")
              write(fn, @sprintf "#%7s%24s\n" "n"  "nu" )
              for i=1: normsDim-1
                 if nu[i]>=0 && nu[i]<=1
                    nuVal=nu[i]
                 elseif nu[i]>1
                    nuVal=1.0
                 else
                    nuVal=0.0
                 end
                 write(fn, @sprintf "%8s%24.12E\n" "$(i)"  nuVal)
              end
              flush(fn)
           end
        end
#__________________________________________________________________________________________

     end

    S1_sp, S1_op, Salpha_sp, Salpha_op, sigma2_n,LogSumPNalpha,H4,H5,H6
end


text:
module LIFs
using LightGraphs
using SpikingNeuralNets: voltage, spike, spiketime
import SpikingNeuralNets: AbstractSNN, potential, thresh
export LIF, memory, channels, haschannel, currentFn, varFn, vars, add_channel!, remove_channel!, I_syn, potential, thresh

"""
A data-type for LIF spiking neural networks.

# Arguments
- `graph::G`: the graph structure of neurons (vertices) and weights (edges)
- `I::Vector{Bool}`: specifies whether each neuron is an input neuron
- `O::Vector{Bool}`: specifies whether each neuron is an output neuron
- `S::Vector{Dict{ST, ST}}`: a memory S[neuron][time] of spike counts
- `V::Vector{Vector{VT}}`: a length m voltage memory V[neuron][time],
    where V[1] is the least recent voltage and V[end] is the current recent voltage
- `m::Integer`: the length of configuration memory
- `dt::VT`: the integration time-step of the network
- `Vrest::Vector{VT}`: the resting voltage of each neuron in the network
- `Vθ::Vector{VT}`: the voltage threshold of each neuron in the network
- `Vreset::Vector{VT}`: the reset voltage of each neuron in the network
- `τ::Vector{VT}`: the membrane time constant of each neuron in the network
- `R::Vector{VT}`: the membrane resistance of each neuron in the network
- `τref::Vector{VT}`: the refractory period of each neuron in the network
- `currentFns::Dict{String, Function}`: the mapping of channel names to the current
    functions I(lif::LIF, n::Integer, s::Vector{ST})::Real for each channel
- `varFns::Dict{String, Function}`: the mapping of channel names to the variable
    update functions s(lif::LIF, n::Integer, s::Vector{ST})::Vector{ST} for each channel
- `vars::Dict{Tuple{String, Integer}, Array{VT}}`: the mapping of channel names
    and neuron index to the (gating) parameters Array{VT}
"""
struct LIF{VT, ST, G} <: AbstractSNN{VT, ST, G}
    graph::G
    I::Vector{Bool}
    O::Vector{Bool}
    S::Vector{Dict{ST, ST}}
    V::Vector{Vector{VT}}
    m::ST
    dt::VT
    Vrest::Vector{VT}
    Vθ::Vector{VT}
    Vreset::Vector{VT}
    τ::Vector{VT}
    R::Vector{VT}
    τref::Vector{VT}
    currentFns::Dict{String, Function}
    varFns::Dict{String, Function}
    vars::Dict{Tuple{String, Integer}, Array{VT}}

    """
        LIF{VT,ST,G}(graph::G; m::ST=ST(1), dt::VT=VT(1), Vrest::Vector{<:VT}=zeros(VT, nv(graph)),
                     Vθ::Vector{<:VT}=zeros(VT, nv(graph)),
                     Vreset::Vector{VT}=zeros(VT, nv(graph)),
                     τ::Vector{<:VT}=ones(VT, nv(graph)),
                     R::Vector{<:VT}=ones(VT, nv(graph)),
                     τref::Vector{<:VT}=zeros(VT, nv(graph)),
                     currentFns::Dict{String, Function}=Dict{String, Function}(),
                     varFns::Dict{String, Function}=Dict{String, Function}(),
                     vars::Dict{<:Tuple{String, <:Integer}, <:Array{<:VT}}=Dict{Tuple{String, Int}, Vector{VT}}(),
                     S::Vector{<:Dict{<:ST, <:ST}}=[Dict{ST, ST}() for _ in vertices(graph)],
                     V::Vector{<:Vector{<:VT}}=[rand(nv(graph)).*(Vθ.-Vrest).+Vrest]) where {VT<:Real, ST<:Integer, G<:AbstractGraph{ST}}

    Create a new LIF{VT,ST,G} with the given graph, parameters, and initial configuration.
    """
    function LIF{VT,ST,G}(graph::G; m::ST=ST(1), dt::VT=VT(1), Vrest::Vector{<:VT}=zeros(VT, nv(graph)),
                          Vθ::Vector{<:VT}=zeros(VT, nv(graph)),
                          Vreset::Vector{VT}=zeros(VT, nv(graph)),
                          τ::Vector{<:VT}=ones(VT, nv(graph)),
                          R::Vector{<:VT}=ones(VT, nv(graph)),
                          τref::Vector{<:VT}=zeros(VT, nv(graph)),
                          currentFns::Dict{String, Function}=Dict{String, Function}(),
                          varFns::Dict{String, Function}=Dict{String, Function}(),
                          vars::Dict{<:Tuple{String, <:Integer}, <:Array{<:VT}}=Dict{Tuple{String, Int}, Vector{VT}}(),
                          S::Vector{<:Dict{<:ST, <:ST}}=[Dict{ST, ST}() for _ in vertices(graph)],
                          V::Vector{<:Vector{<:VT}}=[rand(nv(graph)).*(Vθ.-Vrest).+Vrest]) where {VT<:Real, ST<:Integer, G<:AbstractGraph{ST}}
        # Input neurons are those without incoming connections
        I = map(n -> length(inneighbors(graph, n)) == 0, vertices(graph))
        # Output neurons are those without outgoing connections
        O = map(n -> length(outneighbors(graph, n)) == 0, vertices(graph))
        return new(graph, I, O, S, V, m, dt, Vrest, Vθ, Vreset, τ, R, τref, currentFns, varFns, vars)
    end
end

"""
    LIF{VT,ST}(graph::G; m::ST=ST(1), dt::VT=VT(1), Vrest::Vector{<:VT}=zeros(VT, nv(graph)),
               Vθ::Vector{<:VT}=zeros(VT, nv(graph)),
               Vreset::Vector{VT}=zeros(VT, nv(graph)),
               τ::Vector{<:VT}=ones(VT, nv(graph)),
               R::Vector{<:VT}=ones(VT, nv(graph)),
               τref::Vector{<:VT}=zeros(VT, nv(graph)),
               currentFns::Dict{String, Function}=Dict{String, Function}(),
               varFns::Dict{String, Function}=Dict{String, Function}(),
               vars::Dict{<:Tuple{String, <:Integer}, <:Array{<:VT}}=Dict{Tuple{String, Int}, Vector{VT}}(),
               S::Vector{<:Dict{<:ST, <:ST}}=[Dict{ST, ST}() for _ in vertices(graph)],
               V::Vector{<:Vector{<:VT}}=[rand(nv(graph)).*(Vθ.-Vrest).+Vrest]) where {VT<:Real, ST<:Integer, G<:AbstractGraph{ST}}

Create a new LIF{VT,ST,G} with the given graph and starting configuration.
Shorthand for the more explicit call LIF{VT,ST,G}(graph...)
"""

function LIF{VT,ST}(graph::G; m::ST=ST(1), dt::VT=VT(1), Vrest::Vector{<:VT}=zeros(VT, nv(graph)),
                    Vθ::Vector{<:VT}=zeros(VT, nv(graph)),
                    Vreset::Vector{VT}=zeros(VT, nv(graph)),
                    τ::Vector{<:VT}=ones(VT, nv(graph)),
                    R::Vector{<:VT}=ones(VT, nv(graph)),
                    τref::Vector{<:VT}=zeros(VT, nv(graph)),
                    currentFns::Dict{String, Function}=Dict{String, Function}(),
                    varFns::Dict{String, Function}=Dict{String, Function}(),
                    vars::Dict{<:Tuple{String, <:Integer}, <:Array{<:VT}}=Dict{Tuple{String, Int}, Vector{VT}}(),
                    S::Vector{<:Dict{<:ST, <:ST}}=[Dict{ST, ST}() for _ in vertices(graph)],
                    V::Vector{<:Vector{<:VT}}=[rand(nv(graph)).*(Vreset.-Vrest).+Vrest]) where {VT<:Real, ST<:Integer, G<:AbstractGraph{ST}}
    return LIF{VT,ST,G}(graph; S=S, V=V, m=m, dt=dt, Vrest=Vrest, Vθ=Vθ, Vreset=Vreset,
                        τ=τ, R=R, τref=τref, currentFns=currentFns, varFns=varFns, vars=vars)
end

"""
    memory(lif::LIF)::Integer

Return the size of the network `lif`'s voltage memory.
"""
function memory(lif::LIF{VT, ST})::ST where {VT<:Real, ST<:Integer}
    return lif.m
end

"""
    channels(lif::LIF)::Base.KeySet{String, Dict{String, Function}}

Return a set of channel labels for all channels in the network `lif`.
"""
channels(lif::LIF)::Base.KeySet{String, Dict{String, Function}} = keys(lif.currentFns)

"""
    haschannel(lif::LIF, channel::String)::Bool

Return true if the network `lif` has a channel with the label `channel`.
"""
haschannel(lif::LIF, channel::String)::Bool = haskey(lif.currents, channel)

"""
    haschannel(lif::LIF, channel::String, n::Integer)::Bool

Return true if the neuron `n` in the network `lif` has a channel with the label `channel`.
"""
haschannel(lif::LIF, channel::String, n::Integer)::Bool = haskey(lif.vars, (channel, n))

"""
    currentFn(lif::LIF, channel::String)::Union{Function, Nothing}

Return the current function for the channel `channel` in the network `lif`.
"""
currentFn(lif::LIF, channel::String)::Union{Function, Nothing} = get(lif.currentFns, channel, nothing)

"""
    varFn(lif::LIF, channel::String)::Union{Function, Nothing}

Return the variable update function for the channel `channel` in the network `lif`.
"""
varFn(lif::LIF, channel::String)::Union{Function, Nothing} = get(lif.varFns, channel, nothing)

"""
    vars(lif::LIF{VT}, channel::String, n::Integer)::Union{Array{VT}, Nothing} where {VT<:Real}

Return the (gating) variables for the channel `channel` of neuron `n` in the network `lif`.
"""
function vars(lif::LIF{VT}, channel::String, n::Integer)::Union{Array{VT}, Nothing} where {VT<:Real}
    return get(lif.vars, (channel, n), nothing)
end

"""
    vars(lif::LIF{VT}, channel::String)::Vector{Union{Array{VT}, Nothing}} where {VT<:Real}

Return the (gating) variables for the channel `channel` of all neurons in the network `lif`.
"""
function vars(lif::LIF{VT}, channel::String)::Vector{Array{VT}} where {VT<:Real}
    v = Vector{Array{VT}}(undef, size(lif))
    @inbounds for n in neurons(lif)
        v[n] = vars(lif, channel, n)
    end
    return v
 end

"""
    add_channel!(lif::LIF{VT,ST,G}, name::String, currentFn::Function, varFn::Function,
                 vars::Vector{<:Array{<:VT}}=fill(zeros(VT, 0), size(lif)),
                 neurons::AbstractVector{Bool}=trues(size(lif)))::LIF{VT,ST,G} where {VT<:Real,ST<:Integer,G<:AbstractGraph}

Add a new channel type to the network `lif` with the label `channel`, the current
function `currentFn`, the variable update function `varFn`, the initial variables
`vars`, for every neuron in `neurons`.
"""
function add_channel!(lif::LIF{VT,ST,G}, name::String, currentFn::Function, varFn::Function,
                      vars::Vector{<:Array{<:VT}}=fill(zeros(VT, 0), size(lif)),
                      neurons::AbstractVector{<:Integer}=findall(.!lif.I))::LIF{VT,ST,G} where {VT<:Real,ST<:Integer,G<:AbstractGraph}
    lif.currentFns[name] = currentFn
    lif.varFns[name] = varFn
    for n in neurons
        lif.vars[name, n] = vars[n]
    end
    return lif
end

"""
    remove_channel!(lif::LIF{VT,ST,G}, name::String)::LIF{VT,ST,G} where {VT<:Real,ST<:Integer,G<:AbstractGraph}

Remove all channels from the network `lif` with the label `channel`.
"""
function remove_channel!(lif::LIF{VT,ST,G}, name::String)::LIF{VT,ST,G} where {VT<:Real,ST<:Integer,G<:AbstractGraph}
    delete!(lif.currentFns, name)
    delete!(lif.varFns, name)
    for n in neurons(lif)
        delete!(lif.vars, (name, n))
    end
    return lif
end

"""
    Isyn!(lif::LIF, n::Integer)::Real

Calculate the the total synaptic input to neuron `n` in the network `lif`,
calculating and updating the (gating) variables `lif.vars[channel, n]` for
each `channel` of the neuron.
"""
function Isyn!(lif::LIF, n::Integer)::Real
    I = 0

    @inbounds for channel in channels(lif)
        if haskey(lif.vars, (channel, n))
            # calculate the gating variables
            lif.vars[channel, n] .+= lif.varFns[channel](lif, n, lif.vars[channel, n])
            # calculate the synaptic current
            I += lif.currentFns[channel](lif, n, lif.vars[channel, n])
        end
    end
    return I
end

"""
    potential(lif::LIF, n::Integer)::Real

Calculate the potential of neuron `n` in the network `lif` according
to a leaky integrate-and-fire rule. Mutates the (gating) variables for neuron `n`
while calculating the synaptic input (see `Isyn!`).
"""
function potential(lif::LIF, n::Integer)::Real
    v = voltage(lif, n)
    I = Isyn!(lif, n)  # make sure to update the gating variables
    if !iszero(spike(lif, n)) || spiketime(lif, n) * lif.dt < lif.τref[n]
        return lif.Vreset[n]
    else
        return v + lif.dt * (lif.Vrest[n] - v + (lif.R[n] * I)) / lif.τ[n]
    end
end

"""
    thresh(lif::LIF, n::Integer, x::Real)::Real

The threshold/step transfer function.
"""
thresh(lif::LIF, n::Integer, x::Real)::Real = thresh(x, lif.Vθ[n])
end


text:
@testset "202.happy-number.jl" begin
    @test all(ishappy(i) for i in (1, 7, 10, 13, 19, 23, 28, 31, 123456))
    @test !any(ishappy(i) for i in 2:6)
    @test !any(ishappy(i) for i in (8, 9, 11, 12, 12345))
end


text:
using Printf

weights = [3 6 8 6 3; 7 6 5 7 3; 4 10 4 1 6; 10 4 3 1 2; 6  1  7 3 9]
A = [0 68 50 99 64 2 7 84 0 84;
12 58 97 80 80 78 47 19 65 69;
82 52 113 156 131 75 86 103 70 84;
126 84 71 126 170 104 149 101 70 142;
84 143 120 81 120 198 136 128 149 121;
112 175 178 89 87 170 165 181 210 218]

function cumulative(weights)
    rows, cols = size(weights) # dimensjonene til en matrisen
    ret_matrix = copy(weights) # initialize the return matrix, assuming input matrix is of type Int
    for i = 2:rows # loop through rows, first row is the same
        for j = 1:cols # loop through columns
            if (j==1) # if at first column
                w_min = min(ret_matrix[i-1,j], ret_matrix[i-1,j+1]) # find smalles value from previous row
            elseif (j == cols) # if at last column
                w_min = min(ret_matrix[i-1,j-1], ret_matrix[i-1,j]) # find smalles value from previous row
            else
                w_min = min(ret_matrix[i-1,j-1],ret_matrix[i-1,j], ret_matrix[i-1,j+1])
            end
            ret_matrix[i,j] = w_min + ret_matrix[i,j] # add smallest weight to current element
        end
    end
    return ret_matrix
end

# returns tuple of column and value of biggest weight
# -1 input at either side indicate edge columns
# always pick the "leftest" column
function col_min(col,left,ctr,right) # col is the column at which we are checking which way to go..
    if left < 0
        if ctr <= right
            return col
        else
            return col+1
        end
    elseif right < 0
        if left <= ctr
            return col-1
        else
            return col
        end
    else
        if left <= ctr
            if left <= right
                return col-1
            else
                return col+1
            end
        else
            if ctr <= right
                return col
            else
                return col+1
            end
        end
    end
end

function back_track(path_weights)
    w = copy(path_weights)
    rows, cols = size(path_weights)
    col = argmin(w[rows,:])
    println(col)
    path = []
    push!(path, (rows, col))
    for row in (rows-1):-1:1
        if col == 1 # 1 and cols indicate edge columns, therefore special cases
            # set path column to the next minimum weight column
            col = col_min(col, -1, w[row,col], w[row,col+1])
            # add tuple of the row and column indexes of the weight in current path
            push!(path, (row, col))
        elseif col == cols
            col = col_min(col, w[row,col-1], w[row,col], -1)
            push!(path, (row, col))
        else
            col = col_min(col, w[row,col-1], w[row,col],  w[row,col+1])
            push!(path, (row, col))
        end
    end
    return path
end

# returns tuple of column and value of biggest weight
# -1 input at either side indicate edge columns
# always pick the "leftest" column
function col_min(col,left,ctr,right) # col is the column at which we are checking which way to go..
    if left < 0
        if ctr <= right
            return col
        else
            return col+1
        end
    elseif right < 0
        if left <= ctr
            return col-1
        else
            return col
        end
    else
        if left <= ctr
            if left <= right
                return col-1
            else
                return col+1
            end
        else
            if ctr <= right
                return col
            else
                return col+1
            end
        end
    end
end

function seam_carving(weights)
    w = copy(weights)
    rows, cols = size(weights) # get dimensions
    #energy = Array{Int64}(0, cols) # Total energy per path
    paths = [] # keep track of the path
    path_energy = [] # keep track of total energy per path
    return_path = []
    min_energy = Inf # initialize smallest total energy count to infinity
    # "c" every path starting from column 1, while "col" is at every row of path
    for c = 1:cols
        push!(paths,[])
        col = c
        push!(paths[c], (rows, col)) # first row and column index is already set
        path_energy = w[rows,col] # same for energy
        for row = (rows-1):-1:1
            if col == 1 # 1 and cols indicate edge columns, therefore special cases
                # set path column to the next minimum weight column
                col = col_min(col, -1, w[row,col], w[row,col+1])
                # cumulate weights to form total path energy
                path_energy = path_energy + w[row,col]
                # add tuple of the row and column indexes of the weight in current path
                push!(paths[c], (row, col))
            elseif col == cols
                col = col_min(col, w[row,col-1], w[row,col], -1)
                path_energy = path_energy + w[row,col]
                push!(paths[c], (row, col))
            else
                col = col_min(col, w[row,col-1], w[row,col],  w[row,col+1])
                path_energy = path_energy + w[row,col]
                push!(paths[c], (row, col))
            end
        end
        if path_energy < min_energy
            return_path = paths[c]
            min_energy = path_energy # set current smalles energy path
        end
    end
    return return_path
end

max_stars_repo_name: evenlwanvik-student/algdat

max_issues_repo_name: evenlwanvik-student/algdat

max_forks_repo_name: evenlwanvik-student/algdat


text:
const Vec = Vector{<:Real}

const VECT_LETTER_INDICIES = Dict{Symbol, Int}(
    :x => 1,
    :y => 2,
    :z => 3,
    :w => 4
)

function set_vector_lindex!(index::Pair{Symbol, Int})
    VECT_LETTER_INDICIES[first(index)] = last(index)
end

function Base.getproperty(v::Vec, key::Symbol)
    i = get(VECT_LETTER_INDICIES, key, missing)
    @assert !ismissing(i) "Vector letter index '$key' not found.
    You may want to define others letter indicies for higher dimension vectors.\n"
    try
        return v[i]
    catch e
        throw(error("attempt to access the $(i)th dimension of $(length(v))D Vector."))
    end
end

function Base.setproperty!(v::Vec, key::Symbol, value::Real)
    i = get(VECT_LETTER_INDICIES, key, missing)

    @assert !ismissing(i) "Vector letter index '$key' not found.
    You may want to define other letter indicies for vectors with dimension higher than 4.\n"

    try
        v[i] = value
    catch e
        throw(error("attempt to set component '$key' of $(length(v))D Vector."))
    end
end

extend!(v::Vec, d = 1) = push!(v, [0 for i in 0:d]...)

"""
        mag(v::Vector{<:Real})

Returns the `magnitude` of `v`.
"""
mag(v::Vec) = sqrt(sum(v .^ 2))
"""
        mag2(v::Vector{<:Real})

Returns the `magnitude` squared of `v`.
"""
mag2(v::Vec) = sum(v .^ 2)
"""
        mag²(v::Vector{<:Real})

Returns the `magnitude` squared of `v`.
"""
mag²(v::Vec) = mag2(v)

Base.abs(v::Vec) = mag(v)
"""
        angle(v::Vector{<:Real})

Returns the `angle` (in radians) that a 2D vector `v` makes
with the positive x-axis.
"""
Base.angle(v::Vec) = atan(v[2], v[1])

"""
        rot(v::Vector{<:Real}, φ)

Returns a new vector where its `magnitude` is the same as `v`,
but its angle is the sum of `φ` and the angle of `v`.
"""
rot(v::Vec, φ::Real) = Float64[cos(angle(v) + φ), sin(angle(v) + φ)] * mag(v)

function limit(v::Vec, max::Real)
    m = mag2(v)
    nv = copy(v)
    if m > max ^ 2
        nv = (nv / sqrt(m)) * max
    end
    return nv
end

unitv(v::Vec) = v / abs(v)

"""
        randv()

Returns a vector of magnitude `1` pointing in a random direction.
"""
randv() = withangle(rand(Float64) * 2π)
"""
        withangle(angle [, abs = 1])

Returns a vector of magnitude `abs` and angle `angle`.
"""
withangle(angle::Real, m::Real = 1) = m * [cos(angle), sin(angle)]
"""
        withmag(v, m)

Returns a vector with the same direction as `v`, but with magnitude `m`.
"""
withmag(v::Vec, m::Real) = m * unitv(v)


text:
@testset "Computing distributions" begin
    ts::Vector{Float64} = [-1.29166, -0.0138715, -0.166332, -1.10423, -0.123669, -0.333283, -1.63198, 0.104893, -0.525271, -1.02037, -0.217897, -0.411236, -1.56985, 0.936447, 0.00540328, -0.817672, -0.66606, -0.697607, -0.848608, 0.377567, -1.11741, 1.73203, 1.36521, 0.0153917, -0.14448, -0.939471, -0.397516, -0.372424, -1.29114, 1.09698, -0.176772, 0.075443, -0.750092, -0.718689, -0.890431, -0.255283, -1.82462, 0.838895, 0.465424, -0.59368, -0.324783, -1.10423, 0.596902, -1.11463, 1.2048, 1.05063, -0.109606, -0.264874, -1.60936, -0.0519972, -0.114344, -1.08595, -0.183041, -0.424608, -1.55924, 0.751244, -0.385355, -0.199314, -1.69702, 1.39298, 1.46684, 0.00402413, -0.957294, -0.392595, -0.980378, -0.0983017, -0.65489, -1.06448, 0.360916, -1.6192, 1.50238, 1.76332, 1.29788, 0.0212719, -0.0279064, -0.396213, -1.38781, 0.0035086, -0.109723, -0.850478, -0.609453, -0.206752, -1.23392, 0.993249, -0.342366, 0.176762, -0.554082, -0.865173, -0.393957, -0.464267, -1.08186, 0.688451, -0.846321, 0.706613, -0.475809, -0.0211409, -1.70142, 1.6549, 1.78347, 0.396218, -1.10942, 1.68192, 1.43142, -0.00429295, -0.322581, -1.32556, 0.0219429, -0.0295917, -0.405318, -1.38712, 0.00763867, -0.139915, -0.948333, -0.384921, -0.386374, -1.29076, 1.1008, -0.171684, 0.0695431, -0.756431, -0.715106, -0.88513, -0.334492, -1.71579, 0.564329, -0.481763, -0.407606, -1.03805, 0.465979, -1.32223, 1.46151, 1.80423, -0.0718727, -1.85739, 1.34291, 1.99842, 0.901415, -0.51084, 0.34916, -0.688936, -0.337897, -1.23121, 0.958487, -0.403893, 0.234622, -0.620584, -0.676877, -0.742962, 0.0569321, -1.3654, 1.81456, 0.502851, -0.81959, 1.1458, 0.977908, -0.122892, -0.405326, -1.55857, 0.301308, -0.772338, -0.25191, -1.50797, 1.38587, 0.569416, -1.23556, 1.86485, 1.31952, 0.0751918, -0.774186, -0.672818, -0.980621, -0.165931, -1.87667, 0.260215, -0.831583, -0.213941, -1.70374, 1.35865, 1.42889, 0.0350163, -0.938797, -0.402145, -1.17729, -0.00212732, -0.18619, -1.11735, -0.092832, -0.303358, -1.56958, -0.028755, -0.11704, -0.996922, -0.315869, -0.452005, -1.30119, 1.12457, -0.135106, 0.0117999, -0.924947, -0.447779, -0.992671]

    embedding = embed(ts, 3, 2)

    if invariantset(embedding)
      println("The embedding forms an invariant set")
    else
      println("The embedding does not form an invariant set")
      return
    end

    # Embed using all points except the last (to allow projection of all vertices
    # in the triangulation).
    points, simplex_inds = triangulate(embedding[1:end-1, :])
    image_points = embedding[2:end, :]

    @testset "Simplex vertices matches embedding[1:end-1, :]" begin
      @test points == embedding[1:end-1, :]
    end

    P = markovmatrix(points, image_points, simplex_inds)
    P = mm_parse(points, image_points, simplex_inds)

    P2 = markovmatrixp(points, image_points, simplex_inds)
    @show sum(P2, 2)
    @testset "Markov matrix" begin
        @test all(sum(P, 2) .≈ 1.0)
        @test all(sum(P2, 2) .≈ 1.0)
    end

    invmeasure, inds_nonzero_simplices = invariantdist(P)
    @show invmeasure, inds_nonzero_simplices

    @testset "Invariant measure" begin
        # There must be at least one simplex with nonzero measure; otherwise, the invariant
        # distribution cannot sum to 1.
        @test size(inds_nonzero_simplices, 1) >= 1

        # The invariant distribution must sum to 1 within reasonable tolerance.
        @test sum(invmeasure) ≈ 1.0 atol = 1/10^6
    end

    nonzero_simplices = simplex_inds[inds_nonzero_simplices, :]
    vertices = points
    centroids, radii = SimplexSplitting.centroids_radii2(vertices, nonzero_simplices)
    #@show size(centroids), size(invmeasure)
    n = 3
    nonempty_bins, measure = get_nonempty_bins(centroids, invmeasure, [n, n, n])

    unique_nonempty_bins = unique(nonempty_bins, 1)

    @testset "Bin measure is 1" begin
        @test sum(measure) ≈ 1.0 atol = 1/10^6
    end

    @testset "All centroids are assigned to bins" begin
        @test !any(nonempty_bins .== 0)
    end

    @testset "Joint distribution" begin
        joint = jointdist(nonempty_bins, measure)
        @show joint
    end

    @testset "Marginal distribution" begin
        Py, Pxy, Pxz = marginaldists(unique_nonempty_bins, measure)
        @show Py
        @show Pxy
        @show Pxz
    end

    joint = jointdist(nonempty_bins, measure)
    Py, Pxy, Pxz = marginaldists(unique_nonempty_bins, measure)

    #@show size(joint), size(Py), size(Pxy), size(Pxz)

    @testset "Joint and marginals are the same size" begin
        sizes = [size(distribution, 1) for distribution in (joint, Py, Pxy, Pxz)]
        @test length(unique(sizes)) == 1
    end
end