text:
using GLMakie, VideoIO
using Images, ImageDraw, BoundingSphere
using Printf

"""
function objecttracker(img)

Return image which tracks green colored objects
from the input image

## Details
Steps

- Converts image to HSV space
- Set HSV value range to create binary mask
  which can identify green objects
- Find contours from that binary mask
- Pick the largest contour area and find its
  bounding box
- Draw the bounding box on the input image
- return the resultant image

"""
function objecttracker(
    img,
    h = 150,
    s = 150,
    v = 255,
    lh = 120,
    ls = 80,
    lv = 90,
    boundingboxvar = 10,
)
    hsv_img = HSV.(img)
    channels = channelview(float.(hsv_img))
    hue_img = channels[1, :, :]
    val_img = channels[3, :, :]
    satur_img = channels[2, :, :]
    mask = zeros(size(hue_img))
    # h, s, v = 150,  150, 255
    # h1, s1, v1 = 120,  80, 90
    h, s, v = h, s, v
    h1, s1, v1 = lh, ls, lv
    ex = boundingboxvar
    for ind in eachindex(hue_img)
        if hue_img[ind] <= h && satur_img[ind] <= s / 255 && val_img[ind] <= v / 255
            if hue_img[ind] >= h1 && satur_img[ind] >= s1 / 255 && val_img[ind] >= v1 / 255
                mask[ind] = 1
            end
        end
    end
    cnts1 = find_contours(mask)
    if length(cnts1) != 0
        nums = findmax(map(x -> length(x), cnts1))[2]
        maxcontourarea = cnts1[nums]
        pts = map(x -> [float(x[1]), float(x[2])], maxcontourarea)
        center, radius = boundingsphere(pts)
        center = [round(Int, center[1]), round(Int, center[2])]
        radius = round(Int, radius)
        # draw_contours(img, RGB(1,0,0), cnts1)
        final_image = draw(
            RGB{N0f8}.(img),
            Polygon(
                RectanglePoints(
                    ImageDraw.Point(center[2] - radius - ex, center[1] - radius - ex),
                    ImageDraw.Point(center[2] + radius + ex, center[1] + radius + ex),
                ),
            ),
            RGB(0, 1, 0),
        )
        # final_image = image
    else
        final_image = img
    end
    return final_image
end;

# rotate direction clocwise
function clockwise(dir)
    return (dir) % 8 + 1
end

# rotate direction counterclocwise
function counterclockwise(dir)
    return (dir + 6) % 8 + 1
end

# move from current pixel to next in given direction
function move(pixel, image, dir, dir_delta)
    newp = pixel + dir_delta[dir]
    height, width = size(image)
    if (0 < newp[1] <= height) && (0 < newp[2] <= width)
        if image[newp] != 0
            return newp
        end
    end
    return CartesianIndex(0, 0)
end

# finds direction between two given pixels
function from_to(from, to, dir_delta)
    delta = to - from
    return findall(x -> x == delta, dir_delta)[1]
end

function detect_move(image, p0, p2, nbd, border, done, dir_delta)
    dir = from_to(p0, p2, dir_delta)
    moved = clockwise(dir)
    p1 = CartesianIndex(0, 0)
    while moved != dir ## 3.1
        newp = move(p0, image, moved, dir_delta)
        if newp[1] != 0
            p1 = newp
            break
        end
        moved = clockwise(moved)
    end

    if p1 == CartesianIndex(0, 0)
        return
    end

    p2 = p1 ## 3.2
    p3 = p0 ## 3.2
    done .= false
    while true
        dir = from_to(p3, p2, dir_delta)
        moved = counterclockwise(dir)
        p4 = CartesianIndex(0, 0)
        done .= false
        while true ## 3.3
            p4 = move(p3, image, moved, dir_delta)
            if p4[1] != 0
                break
            end
            done[moved] = true
            moved = counterclockwise(moved)
        end
        push!(border, p3) ## 3.4
        if p3[1] == size(image, 1) || done[3]
            image[p3] = -nbd
        elseif image[p3] == 1
            image[p3] = nbd
        end

        if (p4 == p0 && p3 == p1) ## 3.5
            break
        end
        p2 = p3
        p3 = p4
    end
end

function find_contours(image)
    nbd = 1
    lnbd = 1
    image = Float64.(image)
    contour_list = Vector{typeof(CartesianIndex[])}()
    done = [false, false, false, false, false, false, false, false]

    # Clockwise Moore neighborhood.
    dir_delta = [
        CartesianIndex(-1, 0),
        CartesianIndex(-1, 1),
        CartesianIndex(0, 1),
        CartesianIndex(1, 1),
        CartesianIndex(1, 0),
        CartesianIndex(1, -1),
        CartesianIndex(0, -1),
        CartesianIndex(-1, -1),
    ]

    height, width = size(image)

    for i = 1:height
        lnbd = 1
        for j = 1:width
            fji = image[i, j]
            is_outer = (image[i, j] == 1 && (j == 1 || image[i, j-1] == 0)) ## 1 (a)
            is_hole = (image[i, j] >= 1 && (j == width || image[i, j+1] == 0))

            if is_outer || is_hole
                # 2
                border = CartesianIndex[]

                from = CartesianIndex(i, j)

                if is_outer
                    nbd += 1
                    from -= CartesianIndex(0, 1)

                else
                    nbd += 1
                    if fji > 1
                        lnbd = fji
                    end
                    from += CartesianIndex(0, 1)
                end

                p0 = CartesianIndex(i, j)
                detect_move(image, p0, from, nbd, border, done, dir_delta) ## 3
                if isempty(border) ##TODO
                    push!(border, p0)
                    image[p0] = -nbd
                end
                push!(contour_list, border)
            end
            if fji != 0 && fji != 1
                lnbd = abs(fji)
            end

        end
    end
    return contour_list
end

# a contour is a vector of 2 int arrays
function draw_contour(image, color, contour)
    for ind in contour
        image[ind] = color
    end
end
function draw_contours(image, color, contours)
    for cnt in contours
        draw_contour(image, color, cnt)
    end
end

println("Step 1: open camera")
f = VideoIO.opencamera()

println("Step 2: Set original image shower")

img = read(f)
fig = Figure(size = (1000, 700), title = "Object Tracker")
ax = GLMakie.Axis(
    fig[1, 1],
    aspect = DataAspect(),
    xlabel = "x",
    xlabelcolor = :black,
    ylabel = "y label",
    ylabelcolor = :white,
    title = "Input Image",
    backgroundcolor = :black,
    labelcolor = :white,
)
ax2 = GLMakie.Axis(
    fig[2, 1],
    aspect = DataAspect(),
    xlabel = "x label",
    ylabel = "y label",
    title = "Resultant Image",
    backgroundcolor = :gray80,
)
Textbox(
    fig[1, 2],
    placeholder = "Steps \n - Converts image to HSV space \n - Set HSV value range to create \n binary mask which can identify green objects \n - Find contours from that binary mask \n - Pick the largest contour area and \n find its bounding box \n - Draw the bounding box on the input image \n - return the resultant image \n ",
    bordercolor = RGB(0),
    textcolor = :black,
    textcolor_placeholder = :black,
)
lsgrid = labelslidergrid!(
    fig,
    [
        "Upper Hue Limit",
        "Upper Saturation Limit",
        "Upper Value Limit",
        "Lower Hue Limit",
        "Lower Saturation Limit",
        "Lower Value Limit",
        "BoundingBox Size",
    ],
    [0:1:255, 0:1:255, 0:1:255, 0:1:255, 0:1:255, 0:1:255, 10:1:30];
    formats = [
        x -> "$(round(x, digits = 1))$s" for
        s in ["UpH", "UpS", "UpV", "LowH", "LowS", "LowV", "Pix"]
    ],
    width = 350,
    tellheight = false,
)
fig[2, 2] = lsgrid.layout
sliderobservables = [s.value for s in lsgrid.sliders]
hidedecorations!(ax)
hidedecorations!(ax2)
node = Node(rotr90(img))
node2 = Node(rotr90(img))
makieimg = image!(ax, node)
makieimg2 = image!(ax2, node2)

running = Node(false)

function toggle_running()
    running[] = !running[] # or more complex logic
    set_close_to!(lsgrid.sliders[1], 150)
    set_close_to!(lsgrid.sliders[2], 150)
    set_close_to!(lsgrid.sliders[3], 255)
    set_close_to!(lsgrid.sliders[4], 120)
    set_close_to!(lsgrid.sliders[5], 80)
    set_close_to!(lsgrid.sliders[6], 90)
    set_close_to!(lsgrid.sliders[7], 10)
end
Textbox(
    fig[3, 1],
    placeholder = "Object Tracker using VideoIO.jl,\n Makie.jl and Images.jl",
    bordercolor = RGB(0),
    textcolor = :black,
    textcolor_placeholder = :black,
)
fig[3, 2] = buttongrid = GridLayout()
btn = Button(fig, label = @lift($running ? "Usedefaults" : "Usedefaults "))
buttongrid[1, 1:1] = [btn]

on(btn.clicks) do _
    toggle_running()
end

# h, s, v = 150,  150, 255
# h1, s1, v1 = 120,  80, 90
set_close_to!(lsgrid.sliders[1], 150)
set_close_to!(lsgrid.sliders[2], 150)
set_close_to!(lsgrid.sliders[3], 255)
set_close_to!(lsgrid.sliders[4], 120)
set_close_to!(lsgrid.sliders[5], 80)
set_close_to!(lsgrid.sliders[6], 90)
set_close_to!(lsgrid.sliders[7], 10)
zebra = 0
display(fig)
while isopen(f)
    global zebra = zebra + 1
    read!(f, img)
    data = [150, 150, 255, 120, 80, 90, 10]
    node[] = rotr90(img)
    lift(sliderobservables...) do slvalues...
        data = [slvalues...]
    end

    node2[] = rotr90(
        objecttracker(img, data[1], data[2], data[3], data[4], data[5], data[6], data[7]),
    )
    # save(@sprintf("test_%04d.png", i) , fig.scene)
    sleep(1 / VideoIO.framerate(f))
end

max_stars_repo_name: ashwani-rathee/Sigil

max_issues_repo_name: ashwani-rathee/Sigil

max_forks_repo_name: ashwani-rathee/Sigil


text:
#####
##### unit tests for building blocks
#####

using TrustRegionMethods: ResidualModel, MinimizationModel, cauchy_point, dogleg_boundary,
    ges_kernel, solve_model, unconstrained_optimum

"Return a closure that evaluates to the objective function of a model."
function model_objective(model::ResidualModel)
    @unpack r, J = model
    function(p)
        0.5 * sum(abs2, r .+ J * p)
    end
end

"Basic sanity check for solver results (splat into tail arguments)."
function is_consistent_solver_results(Δ, p, p_norm, on_boundary)
    all(isfinite, p) &&
        (0 ≤ p_norm ≤ Δ) &&
        (norm(p, 2) ≈ p_norm) &&
        ((p_norm ≈ Δ) == on_boundary)
end

@testset "Nonlinear model constructor sanity checks" begin
    r = ones(2)
    J = ones(2, 2)
    @test_throws ArgumentError ResidualModel(fill(Inf, 2), J)
    @test_throws ArgumentError ResidualModel(fill(NaN, 2), J)
    @test_throws ArgumentError ResidualModel(r, fill(-Inf, 2, 2))
    @test_throws ArgumentError ResidualModel(r, fill(NaN, 2, 2))
    @test_throws ArgumentError ResidualModel(r, ones(2, 3))
    @test_throws ArgumentError ResidualModel(r, ones(3, 3))
end

@testset "Cauchy point" begin
    for _ in 1:100
        # random parameters
        n = rand(2:10)
        r = randn(n)
        J = randn(n, n)
        Δ = abs(randn())
        model = ResidualModel(r, J)
        m_obj = model_objective(model)

        # brute force calculation
        pS = - Δ .* normalize(model.J' * model.r, 2)
        opt = Optim.optimize(τ -> m_obj(pS .* τ), 0, 1, Optim.Brent())

        # calculate and compare
        pC, _, _= pC_results = @inferred cauchy_point(Δ, model)
        @test is_consistent_solver_results(Δ, pC_results...)
        @test pC ≈ (Optim.minimizer(opt) .* pS) atol = eps()^0.25 * n

        # estimated decrease invariant
        g_norm = norm(J' * r, 2)
        @test m_obj(zeros(n)) - m_obj(pC) ≥ 0.5 * g_norm * min(Δ, g_norm / norm(J' * J, 2))
    end
end

@testset "dogleg quadratic" begin
    for _ in 1:100
        n = rand(2:10)
        pC = randn(n)
        D = randn(n)
        while dot(D, pC) ≤ 0
            D = randn(n)
        end
        τ = rand()
        Δ = norm(pC .+ τ .* D, 2)
        @test @inferred(dogleg_boundary(Δ, D, pC)) ≈ τ
    end
end

@testset "eigensolver type stability" begin
    # we test this separately because it makes it easier to debug type inference failures
    n = 10
    model = ResidualModel(rand(n), rand(n, n))
    λ, gap, y1, y2 = @inferred ges_kernel(1.0, MinimizationModel(model), I)
    @test λ isa Float64
    @test gap::Float64 ≥ 0
    @test y1 isa Vector{Float64} && y2 isa Vector{Float64}
end

@testset "trust region solver tests" begin
    for _ in 1:100

        # random problem
        n = rand(2:10)
        model = ResidualModel(rand(n), rand(n, n))
        Δ = abs(randn())
        m_obj = model_objective(model)

        # unconstrained optimum and Cauchy point
        pU, pU_norm = @inferred unconstrained_optimum(model)
        @test pU_norm ≥ 0
        pC, _, _ = pC_results = cauchy_point(Δ, model)
        @test is_consistent_solver_results(Δ, pC_results...)

        # Dogleg
        pD, _, _, = pD_results = @inferred solve_model(Dogleg(), Δ, model)
        @test is_consistent_solver_results(Δ, pD_results...)
        @test m_obj(pD) ≤ m_obj(pC) # improve on Cauchy point

        # GES
        pG, _, _ = pG_results = @inferred solve_model(GeneralizedEigenSolver(), Δ, model)
        @test is_consistent_solver_results(Δ, pG_results...)
        @test m_obj(pG) ≤ m_obj(pD) # improve on dogleg
    end
end

@testset "singularities" begin
    # just some basic sanity check to see if the solver can deal with these
    singular1 = ResidualModel(ones(2), ones(2, 2))
    Δ = 1.0
    @test is_consistent_solver_results(Δ, cauchy_point(1.0, singular1)...)
    @test is_consistent_solver_results(Δ, solve_model(Dogleg(), 1.0, singular1)...)
    # FIXME: solver below could do better, just that hard case is not implemented.
    # check for that when it is.
    @test is_consistent_solver_results(Δ, solve_model(GeneralizedEigenSolver(), 1.0,
                                                      singular1)...)
end

@testset "ForwardDiff wrapper" begin
    @testset "different input types" begin
        r = 0:3
        f(x) = x .+ r
        ff = ForwardDiff_wrapper(f, 4)
        J = Diagonal(ones(4))
        for x ∈ (randn(4), zeros(Int, 4), range(1//5, 3//7; length = 4))
            @test ff(x) == (residual = Float64.(x .+ r), Jacobian = J)
        end
    end

    @testset "handle infeasible" begin
        ff = ForwardDiff_wrapper(x -> all(x .> 0) ? x : x .+ NaN, 4)
        @test !all(isfinite, ff(-ones(4)).residual)
        @test ff(ones(4)) == (residual = ones(4), Jacobian = Diagonal(ones(4)))
    end
end

@testset "printing" begin       # just test that printing is defined
    ff = ForwardDiff_wrapper(x -> Diagonal(ones(2)) * x, 2)
    @test repr(ff) isa AbstractString
    res = trust_region_solver(ff, ones(2))
    @test repr(res) isa AbstractString
    @test repr(TrustRegionResult(1, [1.0], nothing, 1.0, false, 99)) isa AbstractString
end

ext: jl


text:
using PredicateComposition
using Test

@testset "PredicateComposition.jl" begin
    tstdata = [(1:i).^2 for i=1:10]
    @test filter(length ≺ 3, tstdata) == [[1],[1,4]]
    @test filter((length ⪰ 3) ⩓ (length ⪯ 5), tstdata) == [[1,4,9],[1,4,9,16],[1,4,9,16,25]]

    counteven(A) = count(iseven,A)
    countodd(A) = count(isodd,A)
    @test filter(MIN(length,counteven) ≺ 3, tstdata) == [[1],[1,4],[1,4,9],[1,4,9,16],[1,4,9,16,25]]
    @test filter(MAX(countodd,counteven) ≺ 3, tstdata) == [[1],[1,4],[1,4,9],[1,4,9,16]]

    @test map(SUM(length,counteven), tstdata) == [1, 3, 4, 6, 7, 9, 10, 12, 13, 15]
    @test SUM(length,counteven).(tstdata) == [1, 3, 4, 6, 7, 9, 10, 12, 13, 15]

    @test filter(length ≣ counteven, tstdata) == []
    @test filter(length ≺ counteven, tstdata) == []
    @test filter(length ⪯ counteven, tstdata) == []
    @test filter(length ≻ counteven, tstdata) == tstdata
    @test filter(length ⪰ counteven, tstdata) == tstdata

    @test filter(length ≣ countodd, tstdata) == [[1]]
    @test filter(length ≺ countodd, tstdata) == []
    @test filter(length ⪯ countodd, tstdata) == [[1]]
    @test filter(length ≻ countodd, tstdata) == tstdata[2:end]
    @test filter(length ⪰ countodd, tstdata) == tstdata

    @test count(iseven ⩔ isodd, 1:100) == 100
    @test count(!iseven ⩓ !isodd, 1:100) == 0
    @test count(iscntrl ⩔ isnumeric, Char.(0:100)) == 42
    @test count(iscntrl ⩔ isnumeric, Char.(0:10000)) == 598
    @test count(isletter ⩔ !isascii, Char.(0:100)) == 30
    @test count(isletter ⩔ !isascii, Char.(0:10000)) == 9925
end


text:
module TestEnvironmentSampler

include("preamble.jl")
include("toyenv.jl")
#using LyceumBase: asvec, ntimesteps, collate

function extract_tid(batch)
    u = unique(batch.actions)
    @assert length(u) == 1
    @assert length(u[1]) == 1
    tid = u[1][1]
    return tid
end

unwrap(x) = (@assert length(x) == 1; x[])
function unwraptraj(τ::Trajectory)
    Trajectory(unwrap.(τ.S), unwrap.(τ.O), unwrap.(τ.A), τ.R, τ.sT, τ.oT, τ.done)
end

@testset "single threaded" begin
    ntimesteps = 30
    env_kwargs = (max_length = 10, reward_scale = 5)
    sampler = EnvironmentSampler(n -> ntuple(i -> ToyEnv(; env_kwargs...), n))

    i = 1
    B = rollout(sampler, ntimesteps, reset! = reset!, nthreads = 1) do a, o
        a .= i
        i += 1
    end
    A = StructArray(B)

    @test length(A) == ntimesteps / env_kwargs.max_length
    @test all(τ -> length(τ) == env_kwargs.max_length, A)

    U = map(unwraptraj, A)

    @test all(U) do τ
        τ.S == [0, cumsum(τ.A[1:end-1])...]
    end
    @test all(U) do τ
        τ.O == 0:length(τ)-1
    end
    @test let start = 1
        all(U) do τ
            pass = τ.A == start:(start+length(τ)-1)
            start += length(τ)
            pass
        end
    end
    @test all(U) do τ
        τ.R == env_kwargs.reward_scale .* τ.S
    end
    @test all(U) do τ
        τ.sT[] == sum(τ.A)
    end
    @test all(U) do τ
        τ.oT[] == length(τ)
    end
    @test all(τ -> τ.done, U)
end

@testset "multi-threaded (nthreads = $nthreads)" for nthreads in (
    Threads.nthreads(),
    div(Threads.nthreads(), 2),
)
    if Threads.nthreads() < 2
        @warn "Cannot test EnvironmentSampler multi-threading with Threads.nthreads() < 2"
    end

    env_kwargs = (
        max_length = 50,
        reward_scale = 5,
        step_hook = _ -> isodd(Threads.threadid()) && busyloop(0.001),
    )
    ntimesteps = 5 * Threads.nthreads() * env_kwargs.max_length
    sampler = EnvironmentSampler(n -> ntuple(i -> ToyEnv(; env_kwargs...), n))

    tcounts = zeros(Int, Threads.nthreads())
    B = rollout(sampler, ntimesteps, reset! = reset!, nthreads = nthreads) do a, o
        a .= Threads.threadid()
        tcounts[Threads.threadid()] += 1
    end
    A = StructArray(B)

    @test length(A) == ntimesteps / env_kwargs.max_length
    @test all(τ -> length(τ) == env_kwargs.max_length, A)

    U = map(unwraptraj, A)

    @test all(U) do τ
        τ.S == [0, cumsum(τ.A[1:end-1])...]
    end
    @test all(U) do τ
        τ.O == 0:length(τ)-1
    end
    @test all(τ -> length(unique(τ.A)) == 1, U)
    @test all(U) do τ
        τ.R == env_kwargs.reward_scale .* τ.S
    end
    @test all(U) do τ
        τ.sT[] == sum(τ.A)
    end
    @test all(U) do τ
        τ.oT[] == length(τ)
    end
    @test all(τ -> τ.done, U)

    nonempty = filter(B -> ntrajectories(B) > 0, sampler.buffers)
    # test that all threads were utilized
    @test length(nonempty) == nthreads
    @test begin
        minn, maxx = extrema(map(ntrajectories, nonempty))
        # test that the number of trajectories collected by each thread differs by at most 2
        maxx - minn <= 2
    end
end

end # module


text:
module MPM

using Distributions
using OBC
using SAMC
import SAMC: energy, propose!, save!, reject!, record

export MPMPrior, MPMParams, MPMSampler, calc_g, gen_points,
    gen_posterior_points, calc_pvals, error_points, predict, error_moments_cube,
    var_error_hists, e_error_hists, e_error_eff, mpm_classifier, sample!

include("mpmcls.jl")

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

type MPMPrior
    D :: Int
    mu :: Vector{Float64}
    mu_sigmas :: Vector{Float64} # Unvariate variances on mu
    kappa :: Float64
    S :: Matrix{Float64}
end

function MPMPrior(;D=2, args...)
    kappa = 6.0
    p = MPMPrior(D, zeros(D), ones(D)*30, kappa, (kappa-1-D)*eye(D,D))
    for (sym, val) in args
        if sym in fieldnames(MPMPrior)
            setfield!(p, sym, val)
        else
            error("Quantity $sym not recognized in MPMPrior constructor")
        end
    end
    p
end

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

type MPMParams
    mu :: Vector{Float64}
    sigma :: Matrix{Float64}
    sigpre :: Matrix{Float64}
    lam :: Matrix{Float64}
    energy :: Float64 # I dunno... does this belong here?
    # d :: Vector{Float64}
    function MPMParams(mu, sigpre, lam)
        @assert size(sigpre) == (size(mu,1),size(mu,1))
        @assert size(lam,1) == size(mu,1)
        new(mu, sigpre'*sigpre, sigpre, lam, -Inf)
    end
end

function copy!(to::MPMParams, from::MPMParams)
    to.mu[:] = from.mu
    to.sigma[:] = from.sigma
    to.sigpre[:] = from.sigpre
    to.lam[:] = from.lam
    to.energy = from.energy
end

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

type MPMSampler <: Sampler
    curr :: MPMParams
    old :: MPMParams
    prior :: MPMPrior
    data :: Matrix{Int}
    usepriors :: Bool
    d :: Vector{Float64}
end

function MPMSampler(prior::MPMPrior, data::Matrix, obj::MPMParams, d::Float64)
    MPMSampler(prior, data, obj, d*ones(size(data,1)))
end

function MPMSampler(prior::MPMPrior, data::Matrix, obj::MPMParams, d::Vector{Float64})
    @assert size(data,1) > size(data,2)
    if !(eltype(data) <: Integer)
        warn("Data is being truncated from type $(eltype(data)) to Int, dataloss might occur")
    end
    MPMSampler(deepcopy(obj), deepcopy(obj), prior, round(Int, data), true, d)
end

const blocks = 3

function propose!(obj::MPMSampler, block::Int, sigma::Float64)
    curr = obj.curr
    prior = obj.prior

    if block == 1 #Modify means
        curr.mu[:] += randn(prior.D) * sigma*0.1

    elseif block == 2

        #Modify covariances
        #temp = curr.sigma .* (0.5 .+ rand(prior.D, prior.D)) #.* curr.sigma ./ mean(curr.sigma)
        #temp = curr.sigma .+ diagm(abs(clamp(diag(curr.sigma),0.001,Inf)) .* (rand(prior.D).-0.5))
        #temp = (temp'*temp).^(1/2)
        #if isposdef(temp) && false #rand(1:100) > 2
            #curr.sigma = temp
        #else
        for i=1:prior.D*prior.D
            curr.sigpre[i] += rand(Normal(0.0,sigma*0.1))
        end
        At_mul_B!(curr.sigma, curr.sigpre, curr.sigpre) # sigma = sigpre'*sigpre

        #curr.sigma = rand(InverseWishart(propmoves.priorkappa,
                #curr.sigma*(propmoves.priorkappa-1-prior.D)))
        #end

    elseif block == 3
        # Modify lambdas
        for i in 1:length(curr.lam)
            if rand(Bool)
                curr.lam[i] += randn() * sigma*0.1
            end
        end

    #elseif block == 3
        #Modify di's
        #curr.d = clamp(curr.d + randn(self.n) * 0.2, 8,12) #FIXME

    end
    nothing
end

function energy(obj::MPMSampler, block::Int=0) #block currently unused
    accum = 0.0
    #likelihoods
    for i in 1:size(obj.data,1), j in 1:size(obj.data,2)
        accum -= logpdf(Poisson(obj.d[i] * exp(obj.curr.lam[j,i])),
                    obj.data[i,j])
    end
    #lambdas
    accum -= sum(logpdf(MultivariateNormal(obj.curr.mu, obj.curr.sigma), obj.curr.lam))

    if obj.usepriors
        # Sigma
        #accum -= sum(logpdf(Normal(0.0, obj.propmoves.premove), obj.sigpre))
        accum -= logpdf(InverseWishart(obj.prior.kappa, obj.prior.S), obj.curr.sigma)
        # mu
        accum -= sum(logpdf(
            MvNormal(obj.prior.mu, obj.prior.mu_sigmas), obj.curr.mu))
    end
    obj.curr.energy = accum
    return accum
end

record(obj::MPMSampler) = deepcopy(obj.curr)

reject!(obj::MPMSampler) = copy!(obj.curr, obj.old)
save!(obj::MPMSampler) = copy!(obj.old, obj.curr)

function gen_points(n, dmean, pt)
    D = size(pt.lam,1)
    newpts = Array(Float64,D,n)
    for i=1:n
        tempmv = rand(MultivariateNormal(pt.mu, pt.sigma))
        for j=1:D
            rate = dmean*exp(tempmv[j])
            if rate > 1000
                newpts[j,i] = rand(Normal(rate,sqrt(rate)))
            else
                try
                    newpts[j,i] = rand(Poisson(rate))
                catch x
                    @show newpts[j,i]
                    @show rate
                    throw(x)
                end
            end
        end
    end

    # Or using lams below is another way to go, but it is not ideal for getting a true
    # 'feel' for the posterior fit
    #
    #for (i,lamind)=enumerate(sample(1:size(pt.lam,2), n))
        #newpts[:,i] = map(x->rand(Poisson(d*exp(x))), pt.lam[:,lamind])
    #end

    return newpts
end

function gen_posterior_points(n, dmean, posterior)
    npost = length(posterior)
    nper = iceil(n/npost)
    pts = Array(Float64, size(posterior[1].lam, 1), nper*npost)
    for i in 1:npost
        pts[:, (i-1)*nper+1:i*nper] = gen_points(nper, dmean, posterior[i])
    end
    return pts
end

function calc_pvals(Ts, points, db, d=10.0)
    numpoints = size(points, 2) # FIXME Which direction are we going again?
    ndims = size(points,1)
    tvals = Array(Float64, ndims, length(Ts), length(db))
    tru_ts = Array(Float64, ndims, length(Ts))
    pvals = similar(tru_ts)
    for i in 1:length(db)
        newpts = gen_points(numpoints, d, db[i])
        for dim=1:ndims
            for (j,T) in enumerate(Ts)
                tvals[dim,j,i] = T(vec(newpts[dim,:]))
            end
        end
    end
    for dim=1:ndims
        for (j,T) in enumerate(Ts)
            tru_ts[dim,j] = T(vec(points[dim,:]))
            pvals[dim,j] = mean(tru_ts[dim,j] .< tvals[dim,j,:]) + 0.5*mean(tru_ts[dim,j] .== tvals[dim,j,:])
        end
    end
    return pvals, tru_ts, tvals
end

function calc_g(points, db, numlam; dmean=10.0)
    # Note: dmean here could actually be d if we knew valid values of d for
    # each point we are calculating the effective density for
    numpts = size(points, 2)
    dims = size(points, 1)
    dblen = length(db)
    res = zeros(numpts)
    lams = Array(Float64, dims, numlam)
    llams = Array(Float64, dims, numlam)
    rlams = Array(Float64, dims, numlam)
    points = floor(points)
    lgpoints = lgamma(points .+ 1)
    for i in 1:dblen # each draw of theta
        curr = db[i]
        rand!(MultivariateNormal(curr.mu, curr.sigma), lams)
        for k=1:numlam*dims
            rlams[k] = dmean * exp(lams[k])
            llams[k] = log(dmean)+lams[k]
        end
        for j in 1:numpts # each point
            accumlam = 0.0
            for s in 1:numlam # each lambda
                accumD = 0.0
                for d in 1:dims # each dimension
                    accumD += points[d,j]*llams[d,s] - rlams[d,s] - lgpoints[d,j]
                end
                accumlam += exp(accumD)
            end
            res[j] += (accumlam/numlam)
        end
    end
    return log(res / length(db))
end

# Utility convenience function
function mpm_model(data; burn=1000, thin=50, d=100.0, usepriors=true)
    D = size(data, 2)
    cov = eye(D,D) .* 0.1
    mu = ones(D)
    prior = MPM.MPMPrior(D=D, kappa=10., S=eye(D)*0.05*10)

    start = MPM.MPMParams(mu, cov,
        clamp(log(data'/d),-8.0,Inf)) #lam
    obj = MPM.MPMSampler(prior, data, deepcopy(start), d)
    obj.usepriors = usepriors
    #mymh_a = MPM.MHRecord(obj_a,burn=burn,thin=thin)
    return MPM.AMWGRecord(obj, blocks, burn=burn,thin=thin)
end

end


text:
function instantiate_linear_equation_of_state(T, α, β)
    eos = LinearEquationOfState(T, α=α, β=β)
    return eos.α == T(α) && eos.β == T(β)
end

function instantiate_roquet_equations_of_state(T, flavor; coeffs=nothing)
    eos = (coeffs == nothing ? RoquetIdealizedNonlinearEquationOfState(T, flavor) :
                               RoquetIdealizedNonlinearEquationOfState(T, flavor, polynomial_coeffs=coeffs))
    return typeof(eos.polynomial_coeffs.R₁₀₀) == T
end

function instantiate_seawater_buoyancy(T, EquationOfState)
    buoyancy = SeawaterBuoyancy(T, equation_of_state=EquationOfState(T))
    return typeof(buoyancy.gravitational_acceleration) == T
end

function density_perturbation_works(arch, T, eos)
    grid = RegularCartesianGrid(T, N=(3, 3, 3), L=(1, 1, 1))
    C = datatuple(TracerFields(arch, grid))
    density_anomaly = ρ′(2, 2, 2, grid, eos, C)
    return true
end

function buoyancy_frequency_squared_works(arch, T, buoyancy)
    grid = RegularCartesianGrid(N=(3, 3, 3), L=(1, 1, 1))
    C = datatuple(TracerFields(arch, grid))
    N² = buoyancy_frequency_squared(2, 2, 2, grid, buoyancy, C)
    return true
end

function thermal_expansion_works(arch, T, eos)
    grid = RegularCartesianGrid(T, N=(3, 3, 3), L=(1, 1, 1))
    C = datatuple(TracerFields(arch, grid))
    α = thermal_expansion(2, 2, 2, grid, eos, C)
    return true
end

function haline_contraction_works(arch, T, eos)
    grid = RegularCartesianGrid(N=(3, 3, 3), L=(1, 1, 1))
    C = datatuple(TracerFields(arch, grid))
    β = haline_contraction(2, 2, 2, grid, eos, C)
    return true
end

@testset "Buoyancy" begin
    println("Testing buoyancy...")

    @testset "Equations of State" begin
        for T in float_types
            @test instantiate_linear_equation_of_state(T, 0.1, 0.3)

            testcoeffs = (R₀₁₀ = π, R₁₀₀ = ℯ, R₀₂₀ = 2π, R₀₁₁ = 2ℯ, R₂₀₀ = 3π, R₁₀₁ = 3ℯ, R₁₁₀ = 4π)
            for flavor in (:linear, :cabbeling, :cabbeling_thermobaricity, :freezing, :second_order)
                @test instantiate_roquet_equations_of_state(T, flavor)
                @test instantiate_roquet_equations_of_state(T, flavor, coeffs=testcoeffs)
            end

            for EOS in EquationsOfState
                @test instantiate_seawater_buoyancy(T, EOS)
            end

            for arch in archs
                @test density_perturbation_works(arch, T, RoquetIdealizedNonlinearEquationOfState())
            end

            for arch in archs
                for EOS in EquationsOfState
                    buoyancy = SeawaterBuoyancy(T, equation_of_state=EOS(T))
                    @test buoyancy_frequency_squared_works(arch, T, buoyancy)
                end

                for buoyancy in (BuoyancyTracer(), nothing)
                    @test buoyancy_frequency_squared_works(arch, T, buoyancy)
                end
            end

            for arch in archs
                for EOS in EquationsOfState
                    @test thermal_expansion_works(arch, T, EOS())
                    @test haline_contraction_works(arch, T, EOS())
                end
            end
        end
    end
end


text:
# This file contains functions for the Radial Big Data Model, following Khezrimotlagh et al. (2019) algorithm (we use KZCT as acronym here)
# D. Khezrimotlagh, J. Zhu and W.D. Cook et al. / European Journal of Operational Research 274 (2019) 1047–1054

"""
    initialsubset(X,Y)

Creates the initial subset as in KZCT algorithm:
1. Determines p the size of the subset and create an empty vector or matrix to fill out
2. Determines the minimum value for each input m and the maximum value for each output s
3. Select a subset with one DMU per input, with DMU's input equal to the minimum value for this input m and one DMU per output
with DMU's output equal to the maximum value for this output s: m + s DMU's selected
4. To fill out the subset up to p DMUs, implement the following algorithm:
    For each unselected DMU_j,
    {u = 0,
        {For each i and k_i,
            If x_ij <= k_ith percentile of x^i then
                u = u + 1}
            {For each r and k'_r
                If y_rj >= k'th percentile of y^r then
                    u = u + 1}
    pre_score_j = u }

Where x^i and y^r represent the ith inputs and the rth outputs of all DMUS.
The indexes k_i and k'_r are changed from 0 to 100 and refer to the percentiles of x^i and y^r respectively.
Sorting DMUs in descending order by the assigned pre-scores, the remaining DMUs (to construct the subsample with size p)
are selected as those having the greatest pre-scores.

"""
function initialsubset(X::Union{Vector,Matrix},Y::Union{Vector,Matrix}, n::Int64, s::Int64, m::Int64)
    # Determines the size of the initial subset Dˢ
    p = convert(Int64, round(sqrt(n)))

    # Select DMU's with the minimum input or the maximum output
    dmus_selected = Vector{Int64}()
    for x in 1:m
        push!(dmus_selected, argmin(X[:,x]))
    end

    for y in 1:s
        push!(dmus_selected, argmax(Y[:,y]))
    end
    unique!(dmus_selected)

    initial_index = collect(1:n)
    dmus_unselected = initial_index[Not(dmus_selected)]

    # Algorithm implementation to fill out the subsample up to the size p

    # select unselected DMUs
    X_unselected = @view X[Not(dmus_selected), :]
    Y_unselected = @view Y[Not(dmus_selected), :]

    # create the pre-scores matrix (one column with the pre-score and the other to store the initial index)
    pre_scores = zeros(length(dmus_unselected))

    @inbounds for i in 1:size(X_unselected, 1)
        for x in 1:m
            xquant = quantile!(X[:,x], 0.01:0.01:1)
            for k in 1:100
                if X_unselected[i,x] <= xquant[k]
                    pre_scores[i] = pre_scores[i] + 1
                end
            end
        end
        for y in 1:s
            yquant = quantile!(Y[:,y], 0.01:0.01:1)
            for k in 1:100
                if Y_unselected[i,y] >= yquant[k]
                    pre_scores[i] = pre_scores[i] + 1
                end
            end
        end
    end

    # sort pre_scores to get the best DMUs in unselected DMUs
    score_perm = sortperm(pre_scores, rev = true)
    dmus_unselected = dmus_unselected[score_perm]
    pre_scores = pre_scores[score_perm, :]

    # add the number of DMUs needed to fill out the subsample
    z = p - length(dmus_selected)
    new_dmus = dmus_unselected[1:z]
    dmus_selected = vcat(dmus_selected, new_dmus)
    dmus_unselected = initial_index[Not(dmus_selected)]

    # Creates the subsets Dˢ and D_excluding_Dˢ
    return dmus_selected, dmus_unselected
end

"""
    bestpracticesfinder(Subset)

Identify best-practices in the subsample selected
"""
function bestpracticesfinder(scores::Vector, orient::Symbol, atol::Float64)
    index_bestpractices = Vector{Int64}()
    for i in 1:length(scores)
        if orient == :Input
            if scores[i] >= 1 - atol
                push!(index_bestpractices, i)
            end
        elseif orient == :Output
            if scores[i] <= 1 + atol
                push!(index_bestpractices, i)
            end
        end
    end
    return index_bestpractices
end

"""
    deabigdata(X, Y)

Compute the big data radial model using data envelopment analysis for inputs X and outputs Y.

# Optional Arguments
- `orient=:Input`: chooses the radially oriented input mode. For the radially oriented output model choose `:Output`.
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `atol=1e-6`: tolerance for DMU to be considered efficient.
- `names`: a vector of strings with the names of the decision making units.
"""
function deabigdata(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
        orient::Symbol = :Input, rts::Symbol = :CRS, slack::Bool = true, atol::Float64 = 1e-6,
        names::Union{Vector{String},Nothing} = nothing,
        optimizer::Union{DEAOptimizer,Nothing} = nothing, progress::Bool = true)

    #=
    The algorithm is such as:
    1. Start
    2. D <- get a sample of DMUs
    3. D^S <- select a subsample of D
    4. B^S <- find best-practices in D^S
    5. E <- find exterior DMUs in D excluding D^S respect to the hull of B^S
    6. If E = {}, then F = B^S and all DMUs are altready evaluated (go to step 7)
    Otherwise
        6.1. F <- find best-practice DMUs in D^S Union E
        6.2. Evaluate DMUs in D excluding (D^S Union E) respoect to the F's hull
    7. End
    =#

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    if nx != ny
        throw(DimensionMismatch("number of@rows in X and Y ($nx, $ny) are not equal"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    n = nx
    disposX = :Strong
    disposY = :Strong

    # Select initial subset Dˢ using Khezrimotlagh et al. (2019) algorithm
    Dˢ, D_excluding_Dˢ = initialsubset(X, Y, n, s, m)

    # Find the best-practices Bˢ in Dˢ
    Dˢ_evaluation = dea(X[Dˢ, :], Y[Dˢ, :], orient = orient, rts = rts,slack = slack,
                    disposX = disposX, disposY = disposY, optimizer = optimizer, progress = progress)

    index_bestpractices = bestpracticesfinder(efficiency(Dˢ_evaluation), orient, atol)

    length(index_bestpractices) > 0 || throw(ErrorException("No efficient DMUs found in the initial subset. Consider increasing the tolerance `atol`"))

    Bˢ = Dˢ[index_bestpractices]

    # Find exterior DMUs in D excluding Dˢ respect to the hull of Bˢ
    D_excluding_Dˢ_evaluation = dea(X[D_excluding_Dˢ, :], Y[D_excluding_Dˢ, :], orient = orient,
                        rts = rts, slack = slack, Xref = X[Bˢ, :], Yref = Y[Bˢ, :],
                        disposX = disposX, disposY = disposY, optimizer = optimizer, progress = progress)

    index_exteriors = bestpracticesfinder(efficiency(D_excluding_Dˢ_evaluation), orient, atol)

    E = D_excluding_Dˢ[index_exteriors]

    # If E is not empty, then we have to find best practice DMUs in D^S Union E, otherwise, F = B_S and all DMUs are already evaluated
    if size(E,1) > 0
        # Create the subset Dˢ_union_E
        Dˢ_union_E = vcat(Dˢ, E)

        sort!(Dˢ_union_E, rev = false)

        Dˢ_union_E_evaluation = dea(X[Dˢ_union_E, :], Y[Dˢ_union_E, :], orient = orient, rts = rts,slack = slack,
                        disposX = disposX, disposY = disposY,optimizer = optimizer, progress = progress)

        # Find the index of best practices DMUs
        index_bestpractices = bestpracticesfinder(efficiency(Dˢ_union_E_evaluation), orient, atol)

        # Best practices DMUs are the efficient DMUs F of the sample D
        F = Dˢ_union_E[index_bestpractices]

        # # We then evaluate the rest of DMUs
        initial_index = collect(1:n)

        D_excluding_Dˢ_union_E = initial_index[Not(Dˢ_union_E)]

        D_excluding_Dˢ_union_E_evaluation = dea(X[D_excluding_Dˢ_union_E, :], Y[D_excluding_Dˢ_union_E, :], orient = orient,
                            rts = rts, slack = slack, Xref = X[F, :], Yref = Y[F, :],
                            disposX = disposX, disposY = disposY, optimizer = optimizer, progress = progress)

        subset1  = Dˢ_union_E
        results1 = Dˢ_union_E_evaluation
        subset2  = D_excluding_Dˢ_union_E
        results2 = D_excluding_Dˢ_union_E_evaluation
    else
        F = Bˢ

        subset1  = Dˢ
        results1 = Dˢ_evaluation
        subset2  = D_excluding_Dˢ
        results2 = D_excluding_Dˢ_evaluation
    end

    # Get results
    resultsperm = sortperm(vcat(subset1, subset2), rev = false)
    effi = vcat(efficiency(results1), efficiency(results2))
    effi = effi[resultsperm]

    if slack
        slackX = vcat(slacks(results1, :X), slacks(results2, :X))
        slackY = vcat(slacks(results1, :Y), slacks(results2, :Y))

        slackX = slackX[resultsperm, :]
        slackY = slackY[resultsperm, :]
    else
        slackX = Array{Float64}(undef, 0, 0)
        slackY = Array{Float64}(undef, 0, 0)
    end

    Xtarget = vcat(targets(results1, :X), targets(results2, :X))
    Ytarget = vcat(targets(results1, :Y), targets(results2, :Y))

    Xtarget = Xtarget[resultsperm, :]
    Ytarget = Ytarget[resultsperm, :]

    lambdas1 = peersmatrix(results1)
    lambdas2 = peersmatrix(results2)

    index_lambda_to_keep = findall(x -> x in F, subset1)
    lambdaeff = vcat(lambdas1[:,index_lambda_to_keep],lambdas2)
    lambdaindex = vcat(subset1, subset2)
    lambda_matrix = lambdaeff[resultsperm, :]

    # Create the sparse matrix for lambdas
    lambdaeff = spzeros(n, n)
    lambdaeff[:, convert.(Int, lambdaindex[index_lambda_to_keep])] = lambda_matrix

    # return results
    return RadialDEAModel(n, m, s, orient, rts, disposX, disposY, names, effi, slackX, slackY, lambdaeff, Xtarget, Ytarget)
end


text:
addprocs(4)
@everywhere include("src/SearchBall.jl")
include("src/other/MyUtils.jl")

using GeometryTypes, Plots, MyUtils, Statistics
plotly()
sb = SearchBall

const dims = 2

# ws = sb.WorldState([[4.5, 3.0]], [[1.0, 1.0], [5.0, 8.0]], [-1.0, -1.0])
ws(myPlayers) = sb.WorldState([[-1.0, -1.0]], myPlayers, [-1.0, -1.0])
res = [90, 60]
n_opps = 5
runs = 10000
fields = SharedMatrix{Float64}[]

# TODO: parallelisieren
@time for r in 1:runs

    new_pos = rand(n_opps, dims)
    @inbounds new_pos[:, 1] .= project.(new_pos[:, 1], 0.0, 1.0, 0.0, 9.0)
    @inbounds new_pos[:, 2] .= project.(new_pos[:, 2], 0.0, 1.0, 0.0, 6.0)

    # TODO: move player away, if they are to close

    field_corners = [Point2(0.0, 0.0), Point2(9.0, 0.0), Point2(9.0, 6.0), Point2(0.0, 6.0)]
    scalar_field, insides = sb.sample(ws([@view(new_pos[i, :]) for i in 1:n_opps]), 1, field_corners, resolution=res) # FIXME: needs config instead of ws
    @inbounds mean_val = median(scalar_field[insides])
    @inbounds scalar_field[.!insides] = mean_val

    push!(fields, scalar_field)
end

scalar_field = mean(fields)'

using JLD
save(string("r", runs,"opps", n_opps,"_", res[1], "x", res[2], ".jld"), "scalar_field", scalar_field)
heatmap(scalar_field, color=:inferno_r, aspect_ratio=:equal)


text:
module OptimalSpacecraftTrajectories

include("./Derivatives.jl")

end


text:
println("**********Testing Cluster**********")

facts("Creating a Cluster") do
  cl1 = Cluster()
  @fact typeof(cl1) => Cluster
  @fact isempty(cl1) => true
end

facts("Indexing Feature Vectors") do
  cl1 = Cluster()
  fv1 = FeatureVector(["word1" => 4, "word2" => 3])
  fv2 = FeatureVector()
  cl1["fv1"] = fv1
  cl1["fv2"] = fv2

  @fact cl1["fv1"] => fv1
  @fact cl1["fv2"] => fv2
  @fact cl1.vector_sum.map => (fv1+fv2).map
  @fact centroid(cl1).map => ((fv1+fv2)/2).map
end

facts("Getting Keys and Values") do
  cl1 = Cluster()
  fv1 = FeatureVector()
  fv2 = FeatureVector()
  cl1["fv1"] = fv1
  cl1["fv2"] = fv2

  @fact [key for key in  keys(cl1)] => ["fv1","fv2"]
  @fact [value for value in  values(cl1)] => [fv1, fv2]
end


text:
module SYNJULIA

mutable struct Position
    x::Int
    y::Int
end

struct Size
    width::Int
    height::Int
end

struct Widget
    width::String
    position::Position
end

end # module


text:
"""
    IdentityMultiple{T} < AbstractMatrix{T} where T

A scalar multiple of the identity matrix of given order and numeric type.

### Fields

- `M` -- uniform scaling operator of type `T`
- `n` -- size of the identity matrix

### Notes

This type can be used to create multiples of the identity of given size. Since
only the multiple and the order are stored, the allocations are minimal.

Internally, the type wraps Julia's lazy multiple of the identity operator,
`UniformScaling`. `IdentityMultiple` subtypes `AbstractMatix`, hence it can be
used in usual matrix arithmetic and for dispatch on `AbstractMatrix`.

The difference between `UniformScaling` and `IdentityMultiple` is that while the
size of the former is generic, the size of the latter is fixed.

### Examples

The easiest way to create an identity multiple is to use the callable version
of `LinearAlgebra.I`:

```jldoctest identitymultiple
julia> using MathematicalSystems: IdentityMultiple

julia> I2 = I(2)
IdentityMultiple{Float64} of value 1.0 and order 2

julia> I2 + I2
IdentityMultiple{Float64} of value 2.0 and order 2

julia> 4*I2
IdentityMultiple{Float64} of value 4.0 and order 2
```

The numeric type (default `Float64`) can be passed as a second argument:

```jldoctest identitymultiple
julia> I2r = I(2, Rational{Int})
IdentityMultiple{Rational{Int64}} of value 1//1 and order 2

julia> I2r + I2r
IdentityMultiple{Rational{Int64}} of value 2//1 and order 2

julia> 4*I2r
IdentityMultiple{Rational{Int64}} of value 4//1 and order 2
```

To create the matrix with a value different from the default (`1.0`), there are
two ways. Either pass the value through the callable `I`, as in

```jldoctest identitymultiple
julia> I2 = I(2.0, 2)
IdentityMultiple{Float64} of value 2.0 and order 2

julia> I2r = I(2//1, 2)
IdentityMultiple{Rational{Int64}} of value 2//1 and order 2
```

Or use the lower level constructor passing the `UniformScaling` (`I`):

```jldoctest identitymultiple
julia> I2 = IdentityMultiple(2.0*I, 2)
IdentityMultiple{Float64} of value 2.0 and order 2

julia> I2r = IdentityMultiple(2//1*I, 2)
IdentityMultiple{Rational{Int64}} of value 2//1 and order 2
```
"""
struct IdentityMultiple{T} <: AbstractMatrix{T}
    M::UniformScaling{T}
    n::Int
end

Base.IndexStyle(::Type{<:IdentityMultiple}) = IndexLinear()
Base.size(𝐼::IdentityMultiple) = (𝐼.n, 𝐼.n)
Base.getindex(𝐼::IdentityMultiple, inds...) = getindex(𝐼.M, inds...)
Base.getindex(𝐼::IdentityMultiple{T}, ind) where {T} =
    rem(ind-1, 𝐼.n+1) == 0 ? 𝐼.M.λ : zero(T)
Base.setindex!(𝐼::IdentityMultiple, X, inds...) = error("cannot store a value in an `Identity`")

Base.:(*)(x::Number, 𝐼::IdentityMultiple) = IdentityMultiple(x * 𝐼.M, 𝐼.n)

function Base.:(*)(𝐼::IdentityMultiple, v::AbstractVector)
    @assert 𝐼.n == length(v)
    return 𝐼.M.λ * v
end

function Base.:(*)(𝐼::IdentityMultiple, A::AbstractMatrix)
    @assert 𝐼.n == size(A, 1)
    return 𝐼.M.λ * A
end

function Base.:(*)(A::AbstractMatrix, 𝐼::IdentityMultiple)
    @assert size(A, 2) == 𝐼.n
    return A * 𝐼.M.λ
end

function Base.:(+)(𝐼1::IdentityMultiple, 𝐼2::IdentityMultiple)
    @assert 𝐼1.n == 𝐼2.n
    return IdentityMultiple(𝐼1.M + 𝐼2.M, 𝐼1.n)
end

function Base.:(*)(𝐼1::IdentityMultiple, 𝐼2::IdentityMultiple)
    @assert 𝐼1.n == 𝐼2.n
    return IdentityMultiple(𝐼1.M * 𝐼2.M, 𝐼1.n)
end

function Base.show(io::IO, ::MIME"text/plain", 𝐼::IdentityMultiple{T}) where T
    print(io, "$(typeof(𝐼)) of value $(𝐼.M.λ) and order $(𝐼.n)")
end

# callable identity matrix
LinearAlgebra.I(n::Int, N=Float64) = IdentityMultiple(one(N)*I, n)

LinearAlgebra.I(λ::Number, n::Int) = IdentityMultiple(λ*I, n)


text:
# time.jl
# some timing tests

# times reported below are all with julia 1.5
# on 2017 iMac (4.2GHz Intel Core i7) with Mojave 10.14.6

using ConstantArrays
#using FillArrays: Ones
using FillArrays: Ones, Trues # for 0.9
using BenchmarkTools

dim = (2^8,2^8,2^4)
x = randn(dim)
y = similar(x)

mask0 = trues(dim)
maskc = ConstantArray(dim)
maskf = Ones{Bool}(dim)
#const Trues = Ones{Bool, N, Axes} where {N, Axes} # lazy trues()
maskt = Trues(dim)
@assert mask0 == maskt

@assert mask0 == maskc == maskf == maskt

# function wrappers for timing tests

function getter(mask, x) # getindex
	return x[mask]
end

function masker!(y, mask, x) # copyto! and getindex, ala getindex!
	copyto!(y, x[mask])
end

function get_check()
	@assert vec(x) ==
		 getter(mask0, x) ==
		 getter(maskc, x) ==
		 getter(maskf, x) ==
		 getter(maskt, x)
end

# x[mask]
# FillArrays is even slower than base here; it needs an appropriate getindex
if false
	@info "getter"
	get_check()
	@btime getter($mask0, $x) # 1 ms 4 alloc 8MiB
	@btime getter($maskc, $x) # 28 ns 2 alloc 80
	@btime getter($maskf, $x) # 1.6 ms 2 alloc 8MiB -> 28 ns 2 alloc 80
	@btime getter($maskt, $x) # 1.6 ms 2 alloc 8MiB -> 28 ns 2 alloc 80
end

function masker_check()
	@assert masker!(y, maskc, x) == masker!(y, maskf, x) ==
		masker!(y, mask0, x) == x
end

if true
	@info "masker!"
	masker_check()
	@btime masker!($y, $mask0, $x) # 2 ms 4 alloc 8MiB
	@btime masker!($y, $maskc, $x) # 3 ms 4 alloc 8MiB => was 1.5× slower :(
								# before i overloaded getindex; then it became:
								# 0.6 ms 2 alloc 80 byte => 3× faster!
	@btime masker!($y, $maskf, $x) # 2.6 ms (2 allocations: 8MiB) -> 0.6 ms 2 alloc 80
	@btime masker!($y, $maskt, $x) # 2.6 ms (2 allocations: 8MiB) -> 0.6 ms 2 alloc 80
end

#=
# here is the key fix (now in module):
Base.getindex(x::AbstractArray{T,D}, ::Trues{D}) where {T,D} = vec(x)
=#

# much faster after providing that "getindex" method:
if false
	get_check()
	@info "updated"
	@btime getter($maskf, $x) # 28 ns 2 alloc 80 :)
	@btime masker!($y, $maskf, $x) # 66 ms 2 alloc 80 :)
end

#include("time-set.jl") # setindex!
#include("time-dot.jl") # .*

nothing


text:
using OsmoticShockModel

exp_bounds = [ 0 3 ; -1 3 ; -1 2.5 ]'

fparam = (pp) -> OsmoticShockModel.SimData( [0,1.0,2.0,2.5,3.0,6.0,7.0,8.0,12.0,14.0] , exp(pp[1]) , exp(pp[2]),exp(pp[3]),exp(pp[4]),exp(pp[5]),exp(pp[6]) )
fpcomp =  (pp) -> fparam( [-1;0.;pp[1];pp[2];3;pp[3]] )

# simulate reference trajectory:
xr_ref_compact = [.5,2.5,1]
sd_ref = fpcomp(xr_ref_compact)
#ed_ref = OsmoticShockModel.create_oshock_exp_data( sd_ref , collect(linspace(0,10,100)) , 0.005 , 0.02 )
ed_ref = OsmoticShockModel.create_oshock_exp_data( sd_ref , collect(linspace(0,10,100)) , 0.001 , 0.01 )
plot_exp(ed_ref)

ssr_a = OsmoticShockModel.eval_param2ssr( ed_ref , fpcomp([2.;2.;1.]) , 0.001 )

xi = [2,2,1]
for zi=1:10000

end


text:
###############################################################################
### Main proximal gradient function

"""
    x, history = solve_proxgrad(pb, x0, extrapolation; tol, itmax, x1, logstep, kwargs...)

    Run the proximal gradient algorithm solving `pb`, starting from `x0`, using extrapolation function `extrapolation` and return
    optimal poit `x` and indicators history `history`.

    kwargs holds arguments that will be passed to `extrapolation`.
"""
function solve_proxgrad(pb::AbstractProblem{T},
                        x0,
                        extrapolation::Function;
                        tol = 1e-9,
                        itmax = 5000,
                        x1 = nothing,
                        printstep = nothing,
                        logstep = 1,
                        linesearch = false,
                        αuser = nothing,
                        saveiter = true,
                        proxgrad_per_it = 1,
                        kwargs...) where {T}
    n = pb.n
    L = get_gradlips(pb)
    μ = get_μ_cvx(pb)
    @show L, μ
    α = 1/L
    !isnothing(αuser) && (α = αuser)

    @assert length(x0) == n
    x = copy(x0)
    x_old = zeros(n).+1e9

    y = copy(x)
    y_old = copy(y)
    memory = Dict{Symbol, Any}(
        :t => 1.,
        :t_old => 1.,
        :μ => μ,
        :L => L,
        kwargs...
    )

    it=1

    Lk = linesearch ? 1e-10 : L

    ## Secondary stuff: history and kpis
    isnothing(printstep) && (printstep = floor(itmax / 30))

    history = Dict{Symbol, Any}()
    history[:iter_f] = SortedDict{Int, Float64}(0 => F(pb, x))
    history[:iter_x] = SortedDict{Int, Any}(0 => x)
    history[:iter_y] = SortedDict{Int, Any}(0 => y)
    history[:iter_step_length] = SortedDict{Int, Any}()
    history[:iter_support_x] = SortedDict{Int, Any}()
    history[:iter_support_y] = SortedDict{Int, Any}()
    cur_support = SortedSet{Int64}()

    # Number of identifications, desidentifications
    history[:num_id] = 0
    history[:num_desid] = 0
    history[:desidtime_to_num] = SortedDict{Int, Int}()
    history[:coord_to_idit] = Dict{Int, Int}()

    println("==== Proximal gradient solve")
    println("Problem type:          ", typeof(pb))
    println("Lipschitz constant:    ", L)
    !linesearch && println("Step size:             ", α)
    linesearch && println("Step size:             ", "backtracking")

    print("\nit       obj.                       step               struct x          struct y")
    while norm(x-x_old)>tol && it<itmax
        printtest = length(setdiff(get_support(T, y_old), get_support(T, y))) > 0
        if mod(it, printstep) == 0 || printtest
            logstring = @sprintf "\n%6d   %.19e  %.5e" it F(pb, x) norm(x-x_old)
            logstring2 = @sprintf "%16s  %16s" sprint(show, collect(get_support(T, x))) sprint(show, collect(get_support(T, y)))

            col = printtest ? :red : :green
            printstyled(logstring*logstring2, color=col)
        end

        x_old = copy(x)
        y_old = copy(y)
        memory[:t_old] = memory[:t]
        it == 1 && !isnothing(x1) && (y_old = copy(x1))

        ## Foward Backward
        if linesearch
            Lk = fwbw_backtracking(pb, x, ∇f(pb, x), Lk)
            α = 1/Lk
        end

        u = x - α * ∇f(pb, x)
        y = prox_αg(pb, u, α)

        ## Extrapolation step
        p = 1/20
        q = (p^2 + (2-p)^2)/2
        r = 4.0
        memory[:t] = (p + sqrt(q+r*memory[:t]^2))/2
        memory[:u] = u
        memory[:α] = α

        x = extrapolation(pb, y, y_old, it, memory)

        cur_supp_x = get_support(T, x)
        cur_supp_y = get_support(T, y)

        ## history
        if mod(it-1, logstep) == 0
            ## relevant iterate
            relevant_iterate = x
            if haskey(memory, :evaluate_funcval)
                relevant_iterate = memory[:z_old]
                @printf "%16s" collect(get_support(T, relevant_iterate))
            end
            history[:iter_f][proxgrad_per_it*it] = F(pb, relevant_iterate)
            saveiter && (history[:iter_x][proxgrad_per_it*it] = relevant_iterate)
            saveiter && (history[:iter_y][proxgrad_per_it*it] = y)
            history[:iter_step_length][proxgrad_per_it*it] = norm(x-x_old)
            history[:iter_support_x][proxgrad_per_it*it] = get_support(T, relevant_iterate)
            history[:iter_support_y][proxgrad_per_it*it] = get_support(T, y)
        end

        if cur_supp_x != last(history[:iter_support_x]).second
            history[:iter_support_x][proxgrad_per_it*it] = cur_supp_x
        end

        if cur_supp_y != last(history[:iter_support_y]).second
            history[:iter_support_y][proxgrad_per_it*it] = cur_supp_y
        end

        ## Support
        old_support = cur_support
        cur_support = get_support(T, x)
        newelems = setdiff(cur_support, old_support)
        # !isempty(newelems) && printstyled("$it \tSupport length: $(length(cur_support))\t\tNew elements    : $(collect(newelems))\n", color=:green)

        delelems = setdiff(old_support, cur_support)
        # !isempty(delelems) && printstyled("$it \tSupport length: $(length(cur_support))\t\tRemoved elements: $(collect(delelems))\n", color=:red)

        history[:num_id] += length(newelems)
        history[:num_desid] += length(delelems)

        ## Count number of ins and outs
        for elem in newelems
            @assert !haskey(history[:coord_to_idit], elem)
            history[:coord_to_idit][elem] = it
        end
        for elem in delelems
            desidtime = it-history[:coord_to_idit][elem]

            history[:desidtime_to_num][desidtime] = get(history[:desidtime_to_num], desidtime, 0) + 1
            delete!(history[:coord_to_idit], elem)
        end

        it += 1
    end
    @printf "\n%6d   %.19e  %.5e\n" it F(pb, x) norm(x-x_old)

    support = get_support(T, x)

    println("")
    println("Final solution nzeros: ", length(support), " \t(thresh: 1e-13)")
    print("Null coordinates: "); print(collect(support))
    println()
    println("# identifications      : ", history[:num_id])
    println("#     - faulty id      : ", history[:num_desid])
    println()
    @assert length(history[:coord_to_idit]) == length(support)
    println("Time of false identification to number of occurences:")
    display(collect(history[:desidtime_to_num]))
    # println("Time of identification if final coordinates:")
    # display(collect(history[:coord_to_idit]))

    return x, history
end

function fwbw_backtracking(pb, x, gradf, L)
    Lk = L
    η = 2

    # println("Backtracking -- Lk = $Lk")
    u = x - (1/Lk) * gradf
    y_ls = prox_αg(pb, u, 1/Lk)

    # @assert f(pb, y_ls) ≤ f(pb, x) + dot(gradf, y_ls-x) + Lk/2 * norm(y_ls - x)^2

    k = 0
    while (f(pb, y_ls) > f(pb, x) + dot(gradf, y_ls-x) + Lk/2 * norm(y_ls - x)^2) && (k < 50)
        Lk *= η

        u = x - (1/Lk) * gradf
        y_ls = prox_αg(pb, u, 1/Lk)

        # @show -f(pb, y_ls) + f(pb, x) + dot(gradf, y_ls-x) + Lk/2 * norm(y_ls - x)^2, Lk

        k += 1
    end

    # printstyled("its, final Lk = $k, $Lk\n", color=:red)

    return Lk
end


text:
@lazylmul UpperTriangular{T, DefaultBandedBlockBandedMatrix{T}} where T
@lazylmul UnitUpperTriangular{T, DefaultBandedBlockBandedMatrix{T}} where T
@lazylmul LowerTriangular{T, DefaultBandedBlockBandedMatrix{T}} where T
@lazylmul UnitLowerTriangular{T, DefaultBandedBlockBandedMatrix{T}} where T

@lazyldiv UpperTriangular{T, DefaultBandedBlockBandedMatrix{T}} where T
@lazyldiv UnitUpperTriangular{T, DefaultBandedBlockBandedMatrix{T}} where T
@lazyldiv LowerTriangular{T, DefaultBandedBlockBandedMatrix{T}} where T
@lazyldiv UnitLowerTriangular{T, DefaultBandedBlockBandedMatrix{T}} where T

@inline hasmatchingblocks(A) =
    cumulsizes(blocksizes(A),1) == cumulsizes(blocksizes(A),2)

#TODO: non-matchin g blocks
isblockbanded(A::AbstractTriangular) =
    isblockbanded(parent(A))
isbandedblockbanded(A::AbstractTriangular) =
    isbandedblockbanded(parent(A))
function blockbandwidths(A::Union{UpperTriangular,UnitUpperTriangular})
    P = parent(A)
    if hasmatchingblocks(P)
        (min(0,blockbandwidths(P,1)), blockbandwidth(P,2))
    else
        blockbandwidths(P)
    end
end
function blockbandwidths(A::Union{LowerTriangular,UnitLowerTriangular})
    P = parent(A)
    if hasmatchingblocks(P)
        (blockbandwidth(P,1), min(0,blockbandwidth(P,2)))
    else
        blockbandwidths(P)
    end
end
subblockbandwidths(A::AbstractTriangular) = subblockbandwidths(parent(A))

triangularlayout(::Type{Tri}, ::ML) where {Tri,ML<:AbstractBlockBandedLayout} = Tri{ML}()

_triangular_matrix(::Val{'U'}, ::Val{'N'}, A) = UpperTriangular(A)
_triangular_matrix(::Val{'L'}, ::Val{'N'}, A) = LowerTriangular(A)
_triangular_matrix(::Val{'U'}, ::Val{'U'}, A) = UnitUpperTriangular(A)
_triangular_matrix(::Val{'L'}, ::Val{'U'}, A) = UnitLowerTriangular(A)

function _matchingblocks_triangular_mul!(::Val{'U'}, UNIT, A::AbstractMatrix{T}, dest) where T
    # impose block structure
    b = PseudoBlockArray(dest, BlockSizes((cumulsizes(blocksizes(A),1),)))

    Bs = blocksizes(A)
    N = nblocks(Bs,1)

    for K = 1:N
        b_2 = view(b, Block(K))
        Ũ = _triangular_matrix(Val('U'), UNIT, view(A, Block(K,K)))
        materialize!(Lmul(Ũ, b_2))
        JR = Block(K+1):blockrowstop(A,K)
        if !isempty(JR)
            muladd!(one(T), view(A, Block(K), JR), view(b,JR), one(T), b_2)
        end
    end
    dest
end

function _matchingblocks_triangular_mul!(::Val{'L'}, UNIT, A::AbstractMatrix{T}, dest) where T
    # impose block structure
    b = PseudoBlockArray(dest, BlockSizes((cumulsizes(blocksizes(A),1),)))

    Bs = blocksizes(A)
    N = nblocks(Bs,1)

    for K = N:-1:1
        b_2 = view(b, Block(K))
        L̃ = _triangular_matrix(Val('L'), UNIT, view(A, Block(K,K)))
        materialize!(Lmul(L̃, b_2))
        JR = blockrowstart(A,K):Block(K-1)
        if !isempty(JR)
            muladd!(one(T), view(A, Block(K), JR), view(b,JR), one(T), b_2)
        end
    end

    dest
end

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{UPLO,UNIT,<:AbstractBlockBandedLayout},
                                   <:AbstractStridedLayout}) where {UPLO,UNIT}
    U,x = M.A,M.B
    @boundscheck size(U,1) == size(x,1) || throw(BoundsError())
    if hasmatchingblocks(U)
        _matchingblocks_triangular_mul!(Val(UPLO), Val(UNIT), triangulardata(U), x)
    else # use default
        materialize!(MulAdd{BandedBlockBandedColumnMajor,
                            typeof(MemoryLayout(typeof(x))),
                            typeof(MemoryLayout(typeof(dest)))}(one(T), U, copy(x), zero(T), dest))
    end
end

for UNIT in ('U', 'N')
    @eval begin
        @inline function materialize!(M::MatLdivVec{<:TriangularLayout{'U',$UNIT,<:AbstractBlockBandedLayout},
                                        <:AbstractStridedLayout})
            U,dest = M.A,M.B
            T = eltype(dest)

            A = triangulardata(U)
            @assert hasmatchingblocks(A)

            @boundscheck size(A,1) == size(dest,1) || throw(BoundsError())

            # impose block structure
            b = PseudoBlockArray(dest, BlockSizes((cumulsizes(blocksizes(A),1),)))

            Bs = blocksizes(A)
            N = nblocks(Bs,1)

            for K = N:-1:1
                b_2 = view(b, Block(K))
                Ũ = _triangular_matrix(Val('U'), Val($UNIT), view(A, Block(K,K)))
                materialize!(Ldiv(Ũ, b_2))

                if K ≥ 2
                    KR = blockcolstart(A, K):Block(K-1)
                    V_12 = view(A, KR, Block(K))
                    b̃_1 = view(b, KR)
                    materialize!(MulAdd(-one(T), V_12, b_2, one(T), b̃_1))
                end
            end

            dest
        end

        @inline function materialize!(M::MatLdivVec{<:TriangularLayout{'L',$UNIT,<:AbstractBlockBandedLayout},
                                        <:AbstractStridedLayout})
            L,dest = M.A, M.B
            T = eltype(dest)
            A = triangulardata(L)
            @assert hasmatchingblocks(A)

            @boundscheck size(A,1) == size(dest,1) || throw(BoundsError())

            # impose block structure
            b = PseudoBlockArray(dest, BlockSizes((cumulsizes(blocksizes(A),1),)))

            Bs = blocksizes(A)
            N = nblocks(Bs,1)

            for K = 1:N
                b_2 = view(b, Block(K))
                L̃ = _triangular_matrix(Val('L'), Val($UNIT), view(A, Block(K,K)))
                materialize!(Ldiv(L̃, b_2))

                if K < N
                    KR = Block(K+1):blockcolstop(A, K)
                    V_12 = view(A, KR, Block(K))
                    b̃_1 = view(b, KR)
                    materialize!(MulAdd(-one(T), V_12, b_2, one(T), b̃_1))
                end
            end

            dest
        end
    end
end

##################
# UpperBandedBlockBandedBlock
#
#   views of the blocks satisfy the BandedMatrix interface, and are memory-compatible
#   with BLASBandedMatrix.
##################

# const UpperBandedBlockBandedBlock{T} = SubArray{T,2,UpperTriangular{T,BandedBlockBandedMatrix{T}},Tuple{BlockSlice1,BlockSlice1},false}
#
#
# isbanded(::UpperBandedBlockBandedBlock) = true
# # Not type stable, but infers Union type so should be "fast"
# MemoryLayout(B::Type{<:UpperBandedBlockBandedBlock}) =
#     ==(parentindices(B)...) ? UpperTriangularLayout(BandedColumnMajor()) : BandedColumnMajor()
#
# function inblockbands(V::UpperBandedBlockBandedBlock)
#     A = parent(V)
#     K_sl, J_sl = parentindices(V)
#     K, J = K_sl.block, J_sl.block
#     -A.l ≤ Int(J-K) ≤ A.u
# end
#
#
# ######################################
# # BandedMatrix interface  for Blocks #
# ######################################
# @inline bandwidths(V::UpperBandedBlockBandedBlock) = subblockbandwidths(parent(V))
#
#
#
# # gives the columns of parent(V).data that encode the block
# blocks(V::UpperBandedBlockBandedBlock)::Tuple{Int,Int} = Int(first(parentindices(V)).block),
#                                                     Int(last(parentindices(V)).block)
#
#
# function tribandeddata(V::UpperBandedBlockBandedBlock)
#     A = parent(parent(V))
#     K_sl, J_sl = parentindices(V)
#     K, J = K_sl.block, J_sl.block
#     tribandeddata(UpperTriangular(view(A, K, J)))
# end
#
#
# @inline inbands_getindex(V::UpperBandedBlockBandedBlock, k::Int, j::Int) =
#     banded_getindex(bandeddata(V), bandwidths(V)..., k, j)
#
# @inline inbands_setindex!(V::UpperBandedBlockBandedBlock, v, k::Int, j::Int) =
#     banded_setindex!(bandeddata(V), bandwidths(V)..., v, k, j)


text:
using Colors, ImageCore, Base.Test

immutable ArrayLF{T,N} <: AbstractArray{T,N}
    A::Array{T,N}
end
Base.linearindexing{A<:ArrayLF}(::Type{A}) = Base.LinearFast()
Base.size(A::ArrayLF) = size(A.A)
Base.getindex(A::ArrayLF, i::Int) = A.A[i]
Base.setindex!(A::ArrayLF, val, i::Int) = A.A[i] = val

immutable ArrayLS{T,N} <: AbstractArray{T,N}
    A::Array{T,N}
end
Base.linearindexing{A<:ArrayLS}(::Type{A}) = Base.LinearSlow()
Base.size(A::ArrayLS) = size(A.A)
Base.getindex{T,N}(A::ArrayLS{T,N}, i::Vararg{Int,N}) = A.A[i...]
Base.setindex!{T,N}(A::ArrayLS{T,N}, val, i::Vararg{Int,N}) = A.A[i...] = val

@testset "ChannelView" begin

@testset "grayscale" begin
    a0 = [Gray(U8(0.2)), Gray(U8(0.4))]
    for (a, VT, LI) in ((copy(a0), Array, Base.LinearFast()),
                       (ArrayLF(copy(a0)), ChannelView, Base.LinearFast()),
                       (ArrayLS(copy(a0)), ChannelView, Base.LinearSlow()))
        v = ChannelView(a)
        @test isa(channelview(a), VT)
        @test Base.linearindexing(v) == LI
        @test isa(colorview(Gray, v), typeof(a))
        @test ndims(v) == 2 - ImageCore.squeeze1
        @test size(v) == (ImageCore.squeeze1 ? (2,) : (1, 2))
        @test eltype(v) == U8
        @test parent(v) === a
        @test v[1] == U8(0.2)
        @test v[2] == U8(0.4)
        @test_throws BoundsError v[0]
        @test_throws BoundsError v[3]
        v[1] = 0.8
        @test a[1] === Gray(U8(0.8))
        @test_throws BoundsError (v[0] = 0.6)
        @test_throws BoundsError (v[3] = 0.6)
        c = similar(v)
        @test isa(c, ChannelView{U8,1,Array{Gray{U8},1}})
        @test length(c) == 2
        c = similar(v, ImageCore.squeeze1 ? 3 : (1,3))
        @test isa(c, ChannelView{U8,1,Array{Gray{U8},1}})
        @test length(c) == 3
        c = similar(v, Float32)
        @test isa(c, ChannelView{Float32,1,Array{Gray{Float32},1}})
        @test length(c) == 2
        c = similar(v, Float16, ImageCore.squeeze1 ? (5,5) : (1,5,5))
        @test isa(c, ChannelView{Float16,2,Array{Gray{Float16},2}})
        @test size(c) == (ImageCore.squeeze1 ? (5,5) : (1,5,5))
    end
end

@testset "RGB, HSV, etc" begin
    for T in (RGB, BGR, RGB1, RGB4, HSV, Lab, XYZ)
        a0 = [T(0.1,0.2,0.3), T(0.4, 0.5, 0.6)]
        for (a, VT) in ((copy(a0), T<:Union{BGR,RGB1,RGB4} ? ChannelView : Array),
                        (ArrayLS(copy(a0)), ChannelView))
            v = ChannelView(a)
            @test isa(channelview(a), VT)
            @test isa(colorview(T, v), typeof(a))
            @test ndims(v) == 2
            @test size(v) == (3,2)
            @test eltype(v) == Float64
            @test parent(v) === a
            @test v[1] == v[1,1] == 0.1
            @test v[2] == v[2,1] == 0.2
            @test v[3] == v[3,1] == 0.3
            @test v[4] == v[1,2] == 0.4
            @test v[5] == v[2,2] == 0.5
            @test v[6] == v[3,2] == 0.6
            @test_throws BoundsError v[0,1]
            @test_throws BoundsError v[4,1]
            @test_throws BoundsError v[2,0]
            @test_throws BoundsError v[2,3]
            v[2] = 0.8
            @test a[1] == T(0.1,0.8,0.3)
            v[2,1] = 0.7
            @test a[1] == T(0.1,0.7,0.3)
            @test_throws BoundsError (v[0,1] = 0.7)
            @test_throws BoundsError (v[4,1] = 0.7)
            @test_throws BoundsError (v[2,0] = 0.7)
            @test_throws BoundsError (v[2,3] = 0.7)
            c = similar(v)
            @test isa(c, ChannelView{Float64,2,Array{T{Float64},1}})
            @test size(c) == (3,2)
            c = similar(v, (3,4))
            @test isa(c, ChannelView{Float64,2,Array{T{Float64},1}})
            @test size(c) == (3,4)
            @test_throws DimensionMismatch similar(v, (5,4))
            c = similar(v, Float32)
            @test isa(c, ChannelView{Float32,2,Array{T{Float32},1}})
            @test size(c) == (3,2)
            c = similar(v, Float16, (3,5,5))
            @test isa(c, ChannelView{Float16,3,Array{T{Float16},2}})
            @test size(c) == (3,5,5)
            @test_throws DimensionMismatch similar(v, Float16, (2,5,5))
        end
    end
    a = reshape([RGB(1,0,0)])  # 0-dimensional
    v = channelview(a)
    @test indices(v) === (Base.OneTo(3),)
end

@testset "Gray+Alpha" begin
    for T in (AGray,GrayA)
        a = [T(0.1f0,0.2f0), T(0.3f0,0.4f0), T(0.5f0,0.6f0)]
        v = ChannelView(a)
        @test isa(channelview(a), ChannelView)
        @test isa(colorview(T, v), Array)
        @test ndims(v) == 2
        @test size(v) == (2,3)
        @test eltype(v) == Float32
        @test parent(v) === a
        @test v[1] == v[1,1] == 0.1f0
        @test v[2] == v[2,1] == 0.2f0
        @test v[3] == v[1,2] == 0.3f0
        @test v[4] == v[2,2] == 0.4f0
        @test v[5] == v[1,3] == 0.5f0
        @test v[6] == v[2,3] == 0.6f0
        @test_throws BoundsError v[0,1]
        @test_throws BoundsError v[3,1]
        @test_throws BoundsError v[2,0]
        @test_throws BoundsError v[2,4]
        v[2] = 0.8
        @test a[1] == T(0.1f0,0.8f0)
        v[2,1] = 0.7
        @test a[1] == T(0.1f0,0.7f0)
        @test_throws BoundsError (v[0,1] = 0.7)
        @test_throws BoundsError (v[3,1] = 0.7)
        @test_throws BoundsError (v[2,0] = 0.7)
        @test_throws BoundsError (v[2,4] = 0.7)
        c = similar(v)
        @test isa(c, ChannelView{Float32,2,Array{T{Float32},1}})
        @test size(c) == (2,3)
        c = similar(v, (2,4))
        @test isa(c, ChannelView{Float32,2,Array{T{Float32},1}})
        @test size(c) == (2,4)
        @test_throws DimensionMismatch similar(v, (3,4))
        c = similar(v, Float64)
        @test isa(c, ChannelView{Float64,2,Array{T{Float64},1}})
        @test size(c) == (2,3)
        c = similar(v, Float16, (2,5,5))
        @test isa(c, ChannelView{Float16,3,Array{T{Float16},2}})
        @test size(c) == (2,5,5)
        @test_throws DimensionMismatch similar(v, Float16, (3,5,5))
    end
end

@testset "Alpha+RGB, HSV, etc" begin
    for T in (ARGB, ABGR, AHSV, ALab, AXYZ,
              RGBA, BGRA, HSVA, LabA, XYZA)
        a = [T(0.1,0.2,0.3,0.4), T(0.5,0.6,0.7,0.8)]
        v = ChannelView(a)
        @test isa(channelview(a), ChannelView)
        @test isa(colorview(T, v), Array)
        @test ndims(v) == 2
        @test size(v) == (4,2)
        @test eltype(v) == Float64
        @test parent(v) === a
        @test v[1] == v[1,1] == 0.1
        @test v[2] == v[2,1] == 0.2
        @test v[3] == v[3,1] == 0.3
        @test v[4] == v[4,1] == 0.4
        @test v[5] == v[1,2] == 0.5
        @test v[6] == v[2,2] == 0.6
        @test v[7] == v[3,2] == 0.7
        @test v[8] == v[4,2] == 0.8
        @test_throws BoundsError v[0,1]
        @test_throws BoundsError v[5,1]
        @test_throws BoundsError v[2,0]
        @test_throws BoundsError v[2,3]
        v[2] = 0.9
        @test a[1] == T(0.1,0.9,0.3,0.4)
        v[2,1] = 0.7
        @test a[1] == T(0.1,0.7,0.3,0.4)
        @test_throws BoundsError (v[0,1] = 0.7)
        @test_throws BoundsError (v[5,1] = 0.7)
        @test_throws BoundsError (v[2,0] = 0.7)
        @test_throws BoundsError (v[2,3] = 0.7)
        c = similar(v)
        @test isa(c, ChannelView{Float64,2,Array{T{Float64},1}})
        @test size(c) == (4,2)
        c = similar(v, (4,4))
        @test isa(c, ChannelView{Float64,2,Array{T{Float64},1}})
        @test size(c) == (4,4)
        @test_throws DimensionMismatch similar(v, (5,4))
        c = similar(v, Float32)
        @test isa(c, ChannelView{Float32,2,Array{T{Float32},1}})
        @test size(c) == (4,2)
        c = similar(v, Float16, (4,5,5))
        @test isa(c, ChannelView{Float16,3,Array{T{Float16},2}})
        @test size(c) == (4,5,5)
        @test_throws DimensionMismatch similar(v, Float16, (3,5,5))
    end
end

end

@testset "ColorView" begin

@testset "grayscale" begin
    _a0 = [U8(0.2), U8(0.4)]
    a0 = ImageCore.squeeze1 ? _a0 : reshape(_a0, (1, 2))
    for (a, VT, LI) in ((copy(a0), Array{Gray{U8}}, Base.LinearFast()),
                        (ArrayLF(copy(a0)), ColorView{Gray{U8}}, Base.LinearFast()),
                        (ArrayLS(copy(a0)), ColorView{Gray{U8}}, Base.LinearSlow()))
        @test_throws ErrorException ColorView(a)
        v = ColorView{Gray}(a)
        @test isa(colorview(Gray,a), VT)
        @test Base.linearindexing(v) == LI
        @test isa(channelview(v), typeof(a))
        @test ndims(v) == 1
        @test size(v) == (2,)
        @test eltype(v) == Gray{U8}
        @test parent(v) === a
        @test v[1] == Gray(U8(0.2))
        @test v[2] == Gray(U8(0.4))
        @test_throws BoundsError v[0]
        @test_throws BoundsError v[3]
        v[1] = 0.8
        @test a[1] === U8(0.8)
        @test_throws BoundsError (v[0] = 0.6)
        @test_throws BoundsError (v[3] = 0.6)
        c = similar(v)
        @test isa(c, ColorView{Gray{U8},1,Array{U8,1}})
        @test length(c) == 2
        c = similar(v, ImageCore.squeeze1 ? 3 : (1,3))
        @test isa(c, ColorView{Gray{U8},1,Array{U8,1}})
        @test length(c) == 3
        c = similar(v, Gray{Float32})
        @test isa(c, ColorView{Gray{Float32},1,Array{Float32,1}})
        @test length(c) == 2
        c = similar(v, Gray{Float16}, ImageCore.squeeze1 ? (5,5) : (1,5,5))
        @test isa(c, ColorView{Gray{Float16},2,Array{Float16,2}})
        @test size(c) == (ImageCore.squeeze1 ? (5,5) : (1,5,5))
        c = similar(v, Float32)
        @test isa(c, Array{Float32, 1})
    end
    # two dimensional images and linear indexing
    _a0 = U8[0.2 0.4; 0.6 0.8]
    a0 = ImageCore.squeeze1 ? _a0 : reshape(_a0, (1, 2, 2))
    for (a, VT, LI) in ((copy(a0), Array{Gray{U8}}, Base.LinearFast()),
                        (ArrayLF(copy(a0)), ColorView{Gray{U8}}, Base.LinearFast()),
                        (ArrayLS(copy(a0)), ColorView{Gray{U8}}, Base.LinearSlow()))
        @test_throws ErrorException ColorView(a)
        v = ColorView{Gray}(a)
        @test isa(colorview(Gray,a), VT)
        @test Base.linearindexing(v) == LI
        @test isa(channelview(v), typeof(a))
        @test ndims(v) == 2
        @test size(v) == (2,2)
        @test eltype(v) == Gray{U8}
        @test parent(v) === a
        @test v[1] == Gray(U8(0.2))
        @test v[2] == Gray(U8(0.6))
        @test_throws BoundsError v[0]
        @test_throws BoundsError v[5]
        v[1] = 0.9
        @test a[1] === U8(0.9)
        @test_throws BoundsError (v[0] = 0.6)
        @test_throws BoundsError (v[5] = 0.6)
    end
end

@testset "RGB, HSV, etc" begin
    for T in (RGB, BGR, RGB1, RGB4, HSV, Lab, XYZ)
        a0 = [0.1 0.2 0.3; 0.4 0.5 0.6]'
        for (a, VT) in ((copy(a0), T<:Union{BGR,RGB1,RGB4} ? ColorView : Array),
                        (ArrayLS(copy(a0)), ColorView))
            @test_throws ErrorException ColorView(a)
            v = ColorView{T}(a)
            @test isa(colorview(T,a), VT)
            @test isa(channelview(v), typeof(a))
            @test ndims(v) == 1
            @test size(v) == (2,)
            @test eltype(v) == T{Float64}
            @test parent(v) === a
            @test v[1] == T(0.1,0.2,0.3)
            @test v[2] == T(0.4,0.5,0.6)
            @test_throws BoundsError v[0]
            @test_throws BoundsError v[3]
            v[2] = T(0.8, 0.7, 0.6)
            @test a == [0.1 0.2 0.3; 0.8 0.7 0.6]'
            @test_throws BoundsError (v[0] = T(0.8, 0.7, 0.6))
            @test_throws BoundsError (v[3] = T(0.8, 0.7, 0.6))
            c = similar(v)
            @test isa(c, ColorView{T{Float64},1,Array{Float64,2}})
            @test size(c) == (2,)
            c = similar(v, 4)
            @test isa(c, ColorView{T{Float64},1,Array{Float64,2}})
            @test size(c) == (4,)
            c = similar(v, T{Float32})
            @test isa(c, ColorView{T{Float32},1,Array{Float32,2}})
            @test size(c) == (2,)
            c = similar(v, T)
            @test isa(c, ColorView{T{Float64},1,Array{Float64,2}})
            @test size(c) == (2,)
            c = similar(v, T{Float16}, (5,5))
            @test isa(c, ColorView{T{Float16},2,Array{Float16,3}})
            @test size(c) == (5,5)
        end
    end
    a = rand(RGBA{U8}, 5, 5)
    vc = channelview(a)
    @test isa(colorview(RGBA, vc), Array{RGBA{U8},2})
    @test_throws ArgumentError colorview(ARGB, vc)
end

@testset "Gray+Alpha" begin
    for T in (AGray,GrayA)
        a = [0.1f0 0.2f0; 0.3f0 0.4f0; 0.5f0 0.6f0]'
        v = ColorView{T}(a)
        @test isa(colorview(T,a), ColorView{T{Float32}})
        @test isa(channelview(v), Array)
        @test ndims(v) == 1
        @test size(v) == (3,)
        @test eltype(v) == T{Float32}
        @test parent(v) === a
        @test v[1] == T(0.1f0, 0.2f0)
        @test v[2] == T(0.3f0, 0.4f0)
        @test v[3] == T(0.5f0, 0.6f0)
        @test_throws BoundsError v[0]
        @test_throws BoundsError v[4]
        v[2] = T(0.8, 0.7)
        @test a[1,2] === 0.8f0
        @test a[2,2] === 0.7f0
        @test_throws BoundsError (v[0] = T(0.8,0.7))
        @test_throws BoundsError (v[4] = T(0.8,0.7))
        c = similar(v)
        @test isa(c, ColorView{T{Float32},1,Array{Float32,2}})
        @test size(c) == (3,)
        c = similar(v, (4,))
        @test isa(c, ColorView{T{Float32},1,Array{Float32,2}})
        @test size(c) == (4,)
        c = similar(v, T{Float64})
        @test isa(c, ColorView{T{Float64},1,Array{Float64,2}})
        @test size(c) == (3,)
        c = similar(v, T{Float16}, (5,5))
        @test isa(c, ColorView{T{Float16},2,Array{Float16,3}})
        @test size(c) == (5,5)
    end
end

@testset "Alpha+RGB, HSV, etc" begin
    for T in (ARGB, ABGR, AHSV, ALab, AXYZ,
              RGBA, BGRA, HSVA, LabA, XYZA)
        a = [0.1 0.2 0.3 0.4; 0.5 0.6 0.7 0.8]'
        v = ColorView{T}(a)
        @test isa(colorview(T,a), ColorView{T{Float64}})
        @test isa(channelview(v), Array)
        @test ndims(v) == 1
        @test size(v) == (2,)
        @test eltype(v) == T{Float64}
        @test parent(v) === a
        @test v[1] == T(0.1,0.2,0.3,0.4)
        @test v[2] == T(0.5,0.6,0.7,0.8)
        @test_throws BoundsError v[0]
        @test_throws BoundsError v[3]
        v[2] = T(0.9,0.8,0.7,0.6)
        @test a[1,2] == 0.9
        @test a[2,2] == 0.8
        @test a[3,2] == 0.7
        @test a[4,2] == 0.6
        @test_throws BoundsError (v[0] = T(0.9,0.8,0.7,0.6))
        @test_throws BoundsError (v[3] = T(0.9,0.8,0.7,0.6))
        c = similar(v)
        @test isa(c, ColorView{T{Float64},1,Array{Float64,2}})
        @test size(c) == (2,)
        c = similar(v, 4)
        @test isa(c, ColorView{T{Float64},1,Array{Float64,2}})
        @test size(c) == (4,)
        c = similar(v, T{Float32})
        @test isa(c, ColorView{T{Float32},1,Array{Float32,2}})
        @test size(c) == (2,)
        c = similar(v, T{Float16}, (5,5))
        @test isa(c, ColorView{T{Float16},2,Array{Float16,3}})
        @test size(c) == (5,5)
    end
end

end

nothing


text:
"""
An implementation of simple method to calculate Ewald energy
"""
function calc_ewald( pw::PWGrid, Sf, Xpos, Nspecies::Int, atm2species,
                     Zv::Array{Float64}; sigma=nothing )
    #
    Ω  = pw.Ω
    r  = pw.r
    Ns = pw.Ns
    G2 = pw.gvec.G2
    Npoints = prod(Ns)
    LatVecs = pw.LatVecs
    #
    # Generate array of distances
    dr = gen_dr_center( r, LatVecs )
    #
    # Generate charge density
    #
    if sigma==nothing
      sigma = 0.25*ones(Nspecies)
    end
    #
    g1  = zeros( Float64, Npoints )
    rho_is = zeros( Float64, Npoints, Nspecies )
    Rho = zeros(Float64, Npoints)
    #
    for isp = 1:Nspecies
        c1 = 2*sigma[isp]^2
        cc1 = sqrt(2*pi*sigma[isp]^2)^3
        for ip=1:Npoints
            g1[ip] = exp(-dr[ip]^2/c1)/cc1
        end
        #
        g1 = Zv[isp] * g1[:]
        #
        ctmp = R_to_G( Ns, g1 )
        for ip = 1:Npoints
            ctmp[ip] = ctmp[ip]*Sf[ip,isp]
        end
        rho_is[:,isp] = real( G_to_R(Ns, ctmp) )
        intrho = sum(rho_is[:,isp])*Ω/Npoints
        @printf("Species %d, intrho: %18.10f\n", isp, intrho)
        Rho[:] = Rho[:] + rho_is[:,isp]
    end
    #
    intrho = sum(Rho)*Ω/Npoints
    @printf("Integrated total Gaussian rho: %18.10f\n", intrho)
    #
    # Solve Poisson equation and calculate Hartree energy
    ctmp = 4.0*pi*R_to_G( Ns, Rho )
    ctmp[1] = 0.0
    for ip = 2:Npoints
        ctmp[ip] = ctmp[ip] / G2[ip]
    end
    phi = real( G_to_R( Ns, ctmp ) )
    Ehartree = 0.5*dot( phi, Rho ) * Ω/Npoints
    #
    Eself = 0.0
    Natoms = size(Xpos)[2]
    for ia = 1:Natoms
        isp = atm2species[ia]
        Eself = Eself + Zv[isp]^2/(2*sqrt(pi))*(1.0/sigma[isp])
    end
    E_nn = Ehartree - Eself
    @printf("Ehartree, Eself, E_nn = %18.10f %18.10f %18.10f\n", Ehartree, Eself, E_nn)
    return E_nn
end


text:
module NipponColors

using Colors
using PyCall

export plot_color, dict_color

const dirc = @__DIR__

function plot_color()
    c = (
        ro = colorant"rgb(12,12,12)",
        sora = colorant"rgb(88,178,220)",
        ukon = colorant"rgb(239,187,36)",
        nae = colorant"rgb(134,192,102)",
        kohbai = colorant"rgb(225,107,140)",
        fujimurasaki = colorant"rgb(138,107,190)",
    )

    return c
end

function dict_color(file::AbstractString = "$(dirc)/../assets/nippon.json"; mode = :hex)
    dn = dict_name(file; mode = mode)
    fn = eval(mode)

    return fn(dn)
end

function hex(D::AbstractDict)
    dict = Dict{String,RGB}()
    for iter in eachindex(D.keys)
        try
            k = D.keys[iter]
            v = D.vals[iter]
            dict[k] = parse(Colorant, v)
        catch
            continue
        end
    end

    return dict
end

function rgb(D::AbstractDict)
    dict = Dict{String,RGB}()
    for iter in eachindex(D.keys)
        try
            k = D.keys[iter]
            v = D.vals[iter]
            dict[k] = parse(Colorant, "rgb($(v[1]),$(v[2]),$(v[3]))")
        catch
            continue
        end
    end

    return dict
end

function dict_name(file = "$(dirc)/../assets/nippon.json"; mode = :hex)
    py"""
    from json import load

    def convert_json(in_fname):
        with open(in_fname, 'r') as fp:
            c1 = load(fp)
            romanji = [d['romanji'].lower() for d in c1]
            kanji = [d['kanji'] for d in c1]
            color = [d['hex'] for d in c1]
            rgb = [d['rgb'] for d in c1]
            hls = [d['hls'] for d in c1]

            dict_list = []
            for i in range(len(color)):
                dict_list.append({
                    'romanji': romanji[i],
                    'kanji': kanji[i],
                    'hex': color[i],
                    'rgb': rgb[i],
                    'hls': hls[i]
                })

        return dict_list
    """
    list = py"convert_json"(file)
    fn = eval(mode)

    return fn(list)
end

function hex(list)
    dict = Dict{String,String}()
    for ele in list
        dict[ele["romanji"]] = ele["hex"]
        dict[ele["kanji"]] = ele["hex"]
    end

    return dict
end

function rgb(list)
    dict = Dict{String,Vector}()
    for ele in list
        dict[ele["romanji"]] = ele["rgb"]
        dict[ele["kanji"]] = ele["rgb"]
    end

    return dict
end

end


text:
using RayTracingWeekend
using Test

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


text:
################
# This file includes all the functions
# for computing PSF using CPU
# Developed by Xuanwen
# xwghua@gmail.com
###############

function fmax_complx(ma,mb)
    return ((ma+mb)+sign(abs(ma)-abs(mb))*(ma-mb))/2
end

function calcPSF(p1, p2, p3, fobj, NA, xspace, yspace, λ, M, n,boundary,calcMode=0)
    println("==Start PSF calculation==")
    ### param initial
    k = 2*pi*n/λ
    alpha = asin(NA/n)
    xlength = length(xspace)
    ylength = length(yspace)
    zeroline = zeros(Complex64,(1, ylength))

    pattern = zeros(Complex64,(self.xlength,self.ylength))
    centerPT = ceil(xlength/2)

    println("Calculation mode is",calcMode)
    if calcMode == 0
        calcOCTA(noundary)
        patternQUAT = fmax_complx()
        patternHALF = fmax_complx()
        PSF_field =

    psf = 0;
    return psf
end

max_stars_repo_name: xwghua/DRiP-public

max_issues_repo_name: xwghua/DRiP-public

max_forks_repo_name: xwghua/DRiP-public


text:
module Utils

include("utils/factor.jl")
include("utils/cutils.jl")
include("utils/simple_graph.jl")

# Helper functions for normalize; may be useful elsewhere
function _entries(xs::Array)
    ys = map(_entries, xs)
    result = []
    for y in ys
        append!(result, y)
    end
    return result
end

_entries(d::Dict) = _entries(collect(map(_entries, values(d))))

_entries(xs::Tuple) = _entries([y for y in xs])

_entries(x::Number) = [x]

_entries(::Nothing) = []

_transform(f::Function, xs::Array) = [_transform(f,x) for x in xs]

_transform(f::Function, d::Dict) = Dict([(k,_transform(f,x)) for (k,x) in d])

_transform(f::Function, xs::Tuple) = tuple([_transform(f,x) for x in xs]...)

_transform(f::Function, x::Number) = f(x)

_transform(::Function, ::Nothing) = nothing

"""
    normalize(xs)

Normalize an array of non-negative reals to sum to 1
"""
function normalize(xs)
    tot = 0.0
    ys = _entries(xs)
    for y in ys
        if y < 0
            error("Negative probability in $(ys)")
        end
        tot += y
    end
    # In the case of learning, it is legitimately possible that a case never
    # happens, so all statistics are zero. Therefore we accept this case.
    if tot == 0.0
        f = x -> 1.0 / length(ys)
    else
        f = x -> x / tot
    end
    return _transform(f, xs)
end

end


text:
# ---
# title: 557. Reverse Words in a String III
# id: problem557
# author: zhwang
# date: 2022-02-15
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/reverse-words-in-a-string-iii/description/>
# hidden: true
# ---
#
# Given a string, you need to reverse the order of characters in each word
# within a sentence while still preserving whitespace and initial word order.
#
# **Example 1:**
#
#
#
#     Input: "Let's take LeetCode contest"
#     Output: "s'teL ekat edoCteeL tsetnoc"
#
#
# **Note:** In the string, each word is separated by single space and there will
# not be any extra space in the string.
#
#
## @lc code=start
using LeetCode
reverse_words_557(s::String) = join([reverse(i) for i in split(s)], " ")

## @lc code=end


text:
# Solver
function check(model::Z3Model)
    result = model.ptr[:check]()
    if result == z3[:sat]
        model.status = :Feasible
    elseif result == z3[:unsat]
        model.status = :Infeasible
    else
        model.status = :Unknown
    end
    return model.status
end

function get_unsat_core(model::Z3Model)
    IIS = Int64[]
    if model.status != :Infeasible
        warn("IIS available for infeasible models only. Make sure to run check(model) first")
    else
        unsat_core = model.ptr[:unsat_core]()
        for (i,c) in enumerate(model.constr)
            if unsat_core[:__contains__](c)
                push!(IIS,i)
            end
        end
    end
    return IIS
end

function checkbound(model::Z3Model,bound)
    push!(model)
    if model.sense == :Min
        untrackedconstraint(model, model.objvar <= bound)
    else
        untrackedconstraint(model, model.objvar >= bound)
    end
    status = check(model)
    pop!(model)
    return status
end

function solve(m::Z3Model,lb=-1e12,ub=1e12; AbsGap=1,Precision=1)
    println("lb = $lb, ub = $ub")
    if m.status == :Unchecked
        status = check(m)
        m.status = status
        if status == :Infeasible
            warn("Infeasible Model")
            return :Infeasible
        end
    end
    rlb = rlog(lb)
    rub = rlog(ub)
    if abs(ub-lb) <= AbsGap
        m.status = :Optimal
        return :Optimal
    elseif rub - rlb > 1
        mid = 10^(0.5*(rub + rlb))
    else
        mid = 0.5*(ub + lb)
    end
    mid = round(mid,Precision)
    status = checkbound(m,mid)
    if status == :Feasible
        solve(m,lb,mid,AbsGap=AbsGap)
    else
        solve(m,mid,ub,AbsGap=AbsGap)
    end
end

function rlog(n)
    n = round(n, 0)
    if n == 0
        return 0
    elseif n < 0
        return round(-log10(abs(n)),0)
    else
        return round(log10(n),0)
    end
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function ps_m5(A)
    return ps_m5!(copy(A))
end

@inline function ps_m5!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    max_memslots=6
    memslots=Vector{Matrix{T}}(undef,max_memslots)
    n=size(A,1)
    for j=1:max_memslots
        memslots[j]=Matrix{T}(undef,n,n)
    end
    # The first slots are precomputed nodes [:A]
    memslots[1]=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: B2 B3 B4 Bb5 B5 Bb6 B6 T2k8
    # Computing B2 with operation: mult
    mul!(memslots[2],memslots[1],memslots[1])
    # Computing B3 with operation: mult
    mul!(memslots[3],memslots[2],memslots[1])
    # Computing B4 with operation: mult
    mul!(memslots[4],memslots[3],memslots[1])
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.013092041015625
    coeff2=0.0109100341796875
    coeff3=-0.009273529052734375
    coeff4=0.008008956909179688
    coeff5=-0.0070078372955322266
    memslots[5] .= coeff2.*memslots[1] .+ coeff3.*memslots[2] .+ coeff4.*memslots[3] .+ coeff5.*memslots[4]
    mul!(memslots[5],true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots[6],memslots[4],memslots[5])
    # Deallocating Bb5 in slot 5
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B5
    coeff1=-0.0390625
    coeff2=0.02734375
    coeff3=-0.0205078125
    coeff4=0.01611328125
    coeff5=1.0
    # Smart lincomb recycle B5
    memslots[6] .= coeff2.*memslots[1] .+ coeff3.*memslots[2] .+ coeff4.*memslots[3] .+ coeff5.*memslots[6]
    mul!(memslots[6],true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots[5],memslots[4],memslots[6])
    # Deallocating B4 in slot 4
    # Deallocating Bb6 in slot 6
    # Computing T2k8 = x*I+x*A+x*B2+x*B3+x*B6
    coeff1=1.0
    coeff2=0.5
    coeff3=-0.125
    coeff4=0.0625
    coeff5=1.0
    # Smart lincomb recycle A
    memslots[1] .= coeff2.*memslots[1] .+ coeff3.*memslots[2] .+ coeff4.*memslots[3] .+ coeff5.*memslots[5]
    mul!(memslots[1],true,I*coeff1,true,true)
    # Deallocating B2 in slot 2
    # Deallocating B3 in slot 3
    # Deallocating B6 in slot 5
    return memslots[1] # Returning T2k8
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
module Taste

export Bitter, Salty, Sour, Sweet, うま味,
    isbitter, issalty, issour, issweet, issavory, measureumami

include("bitter.jl")
include("salty.jl")
include("sour.jl")
include("sweet.jl")
include("うま味.jl")

import Base.svd
function svd(A::うま味)
    println("Bitterness: ", A.concentration * √eps())
    println("Saltiness: ", A.concentration * rand()/1000)
    println("Savoriness: ", 1-rand()*eps()^0.2)
    println("Sourness: ", 0.0)
    println("Sweetness: ", 0.0)
end

end # module


text:
# calls to prepare artificial julia functions accepting C_NULL

arbcall"int acb_mat_approx_eig_qr(acb_ptr E, void * L, acb_mat_t R, const acb_mat_t A, const mag_t tol, slong maxiter, slong prec)"
arbcall"int acb_mat_approx_eig_qr(acb_ptr E, void * L, acb_mat_t R, const acb_mat_t A, void * tol, slong maxiter, slong prec)"

arbcall"int acb_mat_approx_eig_qr(acb_ptr E, acb_mat_t L, void * R, const acb_mat_t A, const mag_t tol, slong maxiter, slong prec)"
arbcall"int acb_mat_approx_eig_qr(acb_ptr E, acb_mat_t L, void * R, const acb_mat_t A, void * tol, slong maxiter, slong prec)"

arbcall"int acb_mat_approx_eig_qr(acb_ptr E, void * L, void * R, const acb_mat_t A, const mag_t tol, slong maxiter, slong prec)"
arbcall"int acb_mat_approx_eig_qr(acb_ptr E, void * L, void * R, const acb_mat_t A, void * tol, slong maxiter, slong prec)"

arbcall"void acb_mat_eig_enclosure_rump(acb_t lambda, void * J, acb_mat_t R, const acb_mat_t A, const acb_t lambda_approx, const acb_mat_t R_approx, slong prec)"

arbcall"int acb_mat_eig_simple_rump(acb_ptr E, void * L, acb_mat_t R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"
arbcall"int acb_mat_eig_simple_rump(acb_ptr E, acb_mat_t L, void * R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"
arbcall"int acb_mat_eig_simple_rump(acb_ptr E, void * L, void * R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"

arbcall"int acb_mat_eig_simple_vdhoeven_mourrain(acb_ptr E, void * L, acb_mat_t R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"
arbcall"int acb_mat_eig_simple_vdhoeven_mourrain(acb_ptr E, acb_mat_t L, void * R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"
arbcall"int acb_mat_eig_simple_vdhoeven_mourrain(acb_ptr E, void * L, void * R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"

arbcall"int acb_mat_eig_simple(acb_ptr E, void * L, acb_mat_t R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"
arbcall"int acb_mat_eig_simple(acb_ptr E, acb_mat_t L, void * R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"
arbcall"int acb_mat_eig_simple(acb_ptr E, void * L, void * R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"


text:
module SpringCollab2020SpeedBasedMusicParamTrigger

using ..Ahorn, Maple

@mapdef Trigger "SpringCollab2020/SpeedBasedMusicParamTrigger" SpeedBasedMusicParamTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    paramName::String="fade", minSpeed::Number=0.0, maxSpeed::Number=1.0, activate::Bool=true)

const placements = Ahorn.PlacementDict(
    "Speed-Based Music Param (Spring Collab 2020)" => Ahorn.EntityPlacement(
        SpeedBasedMusicParamTrigger,
        "rectangle"
    )
)

end