text:
using FileIO
using Images
using LinearAlgebra
using Statistics

"""
A Blob is a mask consisting of blocks of adjacent pixels meeting a threshold.
Only those pixels immediately above, below, left or right are considered to
be adjacent (not diagonals).  Blobs offer easy access to a mask consisting
of pixels in the blob and the perimeter either as a list of steps from
one perimeter point to the next or a list of CartesianIndex with the coordinates
of the perimeter points.
"""
struct Blob
    bounds::CartesianIndices
    mask::BitArray  # Mask of those pixels in the blob
    pstart::CartesianIndex{2} # Start of the perimeter
    psteps::Vector{Tuple{Int,Int}} # Steps around perimeter

    function Blob(bounds::CartesianIndices, mask::BitArray)
        function perimeter(m)  # Computes the perimeter steps
            stps = Tuple{Int,Int}[]
            if prod(size(m)) > 1
                msk(r, c) = checkbounds(Bool, m, r, c) && m[r, c]
                mod8(x) = (x + 7) % 8 + 1 # maintains 1...8
                next(cur, drs) = mod8(drs[findfirst(r -> msk((cur .+ steps[mod8(r)])...), drs)])
                steps = ((1, 0), (1, 1), (0, 1), (-1, 1), (-1, 0), (-1, -1), (0, -1), (1, -1))
                pre = (findfirst(r -> m[r, 1], 1:size(m, 1)), 1)
                @assert msk(pre...)
                ff = findfirst(r -> msk((pre .+ steps[mod8(r)])...), 1:5)
                if isnothing(ff)
                    FileIO.save(File(format"PNG", "dump.png"), mask)
                    @assert !isnothing(ff) "bounds=>$bounds, m[pre]=>$(m[pre...]) m[b]=>$([msk((pre .+ s)...) for s in steps])"
                end
                prevdir = next(pre, 1:5)
                start = pre .+ steps[prevdir]
                @assert msk(start...)
                curr = start
                while true
                    nextdir = next(curr, prevdir-2:prevdir+4)
                    push!(stps, steps[nextdir])
                    curr = curr .+ steps[nextdir]
                    @assert msk(curr...)
                    prevdir = nextdir
                    if curr == start && steps[next(curr, prevdir-2:prevdir+4)] == stps[1]
                        break # end at the start going the same direction
                    end
                end
                return (CartesianIndex{2}(start), stps)
            else
                return (CartesianIndex{2}(1, 1), stps)
            end
        end
        @assert ndims(mask) == 2
        @assert size(bounds) == size(mask)
        ci, stps = perimeter(mask)
        # cip = CartesianIndex([ci.I[i]+bounds.indices[i].start-1 for i in eachindex(ci.I)]...)
        return new(bounds, mask, ci, stps)
    end
end

Base.show(io::IO, b::Blob) = print(io, "Blob[$(b.bounds.indices), area=$(area(b)), perimeter=$(perimeterlength(b))]")

"""
    blob(img::AbstractArray, thresh::Function)::Vector{Blob}

Create a Vector of `Blob`s containing discontiguous regions meeting the threshold function.  The
result is sorted by `Blob` area.
"""
function blob(img::AbstractArray, thresh::Function)::Vector{Blob}
    function extractmask(res, i, minR, maxR, minC, maxC)
        mask = BitArray(undef, maxR - minR + 1, maxC - minC + 1)
        foreach(ci -> mask[ci] = (res[(ci.I .+ (minR - 1, minC - 1))...] == i), CartesianIndices(mask))
        return mask
    end
    res = zeros(UInt16, size(img))
    alias = Vector{Set{eltype(res)}}()
    prev, next = zero(eltype(res)), zero(eltype(res))
    for ci in CartesianIndices(img)
        ci.I[1] == 1 && (prev = 0) == 0  # Reset prev each column
        if thresh(img[ci])
            above = ci[2] > 1 ? res[(ci.I .- (0, 1))...] : zero(eltype(res))
            if above ≠ 0
                if (prev ≠ 0) && (above ≠ prev)
                    # Same blob => different indexes
                    ia = findfirst(al -> above in al, alias)
                    ip = findfirst(al -> prev in al, alias)
                    @assert !(isnothing(ia) || isnothing(ip)) "$prev or $above not in $alias"
                    if ia ≠ ip # merge regions
                        union!(alias[ia], alias[ip])
                        deleteat!(alias, ip)
                    end
                end
                prev = above
            elseif prev == 0
                prev = (next += one(next))
                push!(alias, Set{UInt16}(prev)) # new region
            end
            @assert prev ≠ 0
            res[ci] = prev
        else
            prev = zero(eltype(res))
        end
    end
    # Second pass combine the adjacent indices
    newidx, nblobs = zeros(eltype(res), next), length(alias)
    for i = 1:nblobs, id in alias[i]
        newidx[id] = i
    end
    rects = Dict{eltype(res),NTuple{4,Int}}(i => (100000, -1, 100000, -1) for i = 1:nblobs)
    for ci in CartesianIndices(img)
        if (bidx = res[ci]) ≠ zero(eltype(res)) # belongs to a blob
            ni = (res[ci] = newidx[bidx])
            rect = rects[ni]
            rects[ni] = (min(ci[1], rect[1]), max(ci[1], rect[2]), min(ci[2], rect[3]), max(ci[2], rect[4]))
        end
    end
    blobs = [
        Blob(CartesianIndices((rect[1]:rect[2], rect[3]:rect[4])), extractmask(res, i, rect...)) for (i, rect) in rects
    ]
    sort!(blobs, lt = (b1, b2) -> area(b1) > area(b2))
    return blobs
end

"""
    Base.CartesianIndices(b::Blob)

Bounds of the blob in the original image's coordinate system
"""
Base.CartesianIndices(b::Blob) = b.bounds

"""
    Base.getindex(b::Blob, ci::CartesianIndex)

Whether a pixel in the original image's coordinate system is in the blob.
"""
Base.getindex(b::Blob, ci::CartesianIndex{2}) = #
    (ci in b.bounds) && b.mask[map(i -> ci.I[i] - b.bounds.indices[i].start + 1, eachindex(ci.I))...]

"""
    perimeter(b::Blob)::Vector{CartesianIndex{2}}

Returns a vector of `CartesianIndex{2}` corresponding to the points around the
perimeter of the blob in the original image's coordinate system.
"""
function perimeter(b::Blob)::Vector{CartesianIndex{2}}
    pts, acc = CartesianIndex{2}[b.pstart], [b.pstart.I...]
    foreach(stp -> push!(pts, CartesianIndex{2}((acc .+= stp)...)), b.psteps[1:end-1])
    return pts
end

function perimeter(b::Blob, img::AxisArray, p::Real=2)
    perim = perimeter(b)
    sum( norm(img, perim[i], perim[i == length(p) ? 1 : i+1], p) for i in eachindex(p) )
end

"""
    perimeterlength(b::Blob)

Compute the length of the blob perimeter in pixels.  Diagonals are √2 and straights are 1.
"""
perimeterlength(b::Blob) = length(b.psteps) > 0 ? mapreduce(st -> sqrt(dot(st, st)), +, b.psteps) : 4.0

"""
    ecd(b::Blob, filled = true)

Computes the equivalent circular diameter (by default computes the ecd(..) based on the area including the area
of any interior holes.)  A = πr² => ecd = 2r = 2√(A/π)
"""
ecd(b::Blob, filled = true) = 2.0 * sqrt((filled ? filledarea(b) : area(b)) / π)

"""
    curvature(b::Blob, n::Int)

Compute an array that measures the angular difference between the pixel n before
and n after then current one for each point on the perimeter of the blob.
Negative indicates convex and positive indicates concave with large positive
values indicating a sharp concave angle.
"""
function curvature(b::Blob, n::Int)
    modn(i) = (i + length(b.psteps) - 1) % length(b.psteps) + 1
    stepsum(itr) = mapreduce(j -> b.psteps[modn(j)], (x, y) -> .+(x, y), itr, init = (0, 0))
    angles = Float64[]
    for i in eachindex(b.psteps)
        sm, sp = -1 .* stepsum(i-1:-1:i-n), stepsum(i:i+n-1)
        den² = dot(sm, sm) * dot(sp, sp)
        ac = den² > 0 ? dot(sm, sp) / sqrt(den²) : 1.0
        @assert (ac < 1.00001) && (ac > -1.00001) "ac=$ac"
        c = sign(sm[1] * sp[2] - sm[2] * sp[1]) / acos(min(1.0, max(-1.0, ac)))
        push!(angles, c)
    end
    return angles
end

"""
    splitblob(b::Blob, p1::CartesianIndex{2}, p2::CartesianIndex{2})

Split a Blob by drawing a line from p1 to p2 (assumed to be on the perimeter
or outside b) and reblobining.
"""
function splitblob(b::Blob, p1::CartesianIndex{2}, p2::CartesianIndex{2})
    mask = copy(b.mask)
    # Draw a line to divide the particles for reblobbing
    drawline(pt->mask[pt...]=false, p1, p2, true)
    res = blob(mask, p -> p)
    # Fix up the hoz and vert
    return map(b2 -> Blob(__offset(b2.bounds, b.bounds), b2.mask), res)
end

"""
    separate(b::Blob, concavity = 0.42, withinterior = true)::Vector{Blob}

Break b into many Blob(s) by looking for concave regions on the perimeter
and joining them with a short(ish) line.
"""
function separate(b::Blob, concavity = 0.42, withinterior = true)::Vector{Blob}
    modn(bb, i) = (i + length(bb.psteps) - 1) % length(bb.psteps) + 1
    # In the concatity vector, find peaks of positive regions
    function findmaxes(c::Vector{Float64}, minwidth::Int, concav::Float64, sym::Int)::Vector{Tuple{Int,Int}}
        first, st, maxi, res = c[1] > 0.0 ? -1 : 0, 0, 0, Vector{Tuple{Int,Int}}()
        for i in eachindex(c)
            if first < 0 # First keeps track of (possible) initial region above zero
                first = c[i] > c[-first] ? -i : (c[i] < 0.0 ? -first : first)
            end
            if st == 0 # In region below zero concavity
                if c[i] > 0.0
                    st, maxi = i, i
                end
            else # In region of positive concavity
                maxi = c[i] > c[maxi] ? i : maxi
                if c[i] < 0.0
                    if i - st >= minwidth && c[maxi] >= concav
                        push!(res, (maxi, sym))
                    end
                    st, maxi = 0, 0 # start new region of below zero concavity
                end
            end
        end
        if first > 0 # Deal with initial region
            lenfirst = (st ≠ 0 ? length(c) - st : 0) + findfirst(i -> c[i] < 0.0, eachindex(c))
            maxi = st ≠ 0 ? (c[maxi] > c[first] ? maxi : first) : first
            if lenfirst >= minwidth && c[maxi] >= concav
                push!(res, (maxi, sym))
            end
        end
        return res
    end
    besti = -1
    if length(b.psteps) > 20
        p, c = perimeter(b), curvature(b, 4)
        maxes = findmaxes(c, 4, concavity, 0) # maximum
        if withinterior
            for (i, ir) in enumerate(interiorregions(b))
                if area(ir) > 20
                    append!(p, perimeter(ir))
                    c = 0.0 .- curvature(ir, 4)
                    append!(maxes, findmaxes(c, 4, concavity, i))
                end
            end
        end
        # Find pairs of concavities to use as splitters
        bestj, bestlen = -1, 100_000_000::Int
        for (i, (ii, symi)) in enumerate(maxes), (ij, symj) in maxes[i+1:end]
            # Only break exterior to exterior, exterior to interior, interior to interior or from one interior to a different interior
            if symi == 0 || symj == 0 || symi != symj
                len = dot((p[ii].I .- p[ij].I),(p[ii].I .- p[ij].I))::Int
                # pick the shortest splitter
                if len < bestlen
                    besti, bestj, bestlen = ii, ij, len
                end
            end
        end
    end
    # There is a split point so split it and recursively separate the splits
    return besti == -1 ? [b] : mapreduce(separate, append!, splitblob(b, p[besti], p[bestj]), init = Blob[])
end

"""
    scorer(b::Blob)

A default function to score a blob as a candidate particle.  Smaller scores are more particle like.
"""
scorer(b::Blob, img::AbstractArray, minarea = 100) = # perimeter/π == ecd for a circle
    (area(b) < minarea ? 100.0 : minarea / area(b)) + perimeterlength(b) / (π * ecd(b, false))

"""
    commonarea(b1::Blob, b2::Blob)

Returns the number of pixels in common between `b1` and `b2`.
"""
commonarea(b1::Blob, b2::Blob) = count(ci -> b1[ci] && b2[ci], intersect(b1, b2))

"""
    multiseparate(img::Array, threshes, score; concavity=0.42, minarea = 10, dump = nothing)

Uses multiple thresholds to attempt to find the best separation of the distinct blobs in the image.  The best blob b
is defined as the one that produces particles that produce smaller `score(b)`.  So a blob will be split if splitting
the blob will produce multiple blobs of lower scores. The default function 'scorer(b::Blob)' looks for more circular
blobs.
"""
function multiseparate(img::AbstractArray, threshes; score = scorer, concavity = 0.42, minarea = 10, dump = nothing)::Vector{Blob}
    function segmenter(th)
        starters = blob(img, p -> p >= th)
        return length(starters) > 0 ? #
               filter(b -> area(b) > minarea, mapreduce(b -> separate(b, concavity), append!, starters)) : #
               Blob[]
    end
    best, cx = segmenter(threshes[1]), 0
    for th in threshes[2:end]
        blobs, newbest = segmenter(th), Blob[]
        # compare the new blobs to the best previous ones
        for bb in best
            # Find which `blobs` make up `bb`
            becomes = filter(b -> commonarea(b, bb) / area(b) > 0.8, blobs)
            # Should we split `bb` into `becomes`
            split = length(becomes) > 1 && mean(map(b->score(b, img), becomes)) < score(bb, img)
            if (!isnothing(dump)) && (length(becomes) > 1)
                open(dump * "[details].txt", "a") do io
                    write(
                        io,
                        "$(basename(dump))[$(cx+=1)][split=$split]: before = $(score(bb)), becomes = $(score.(becomes)), pl = $(perimeterlength(bb)), ecd=$(ecd(bb,false))\n",
                    )
                end
                FileIO.save(
                    File(format"PNG", dump * "[$cx, before, $split].png"),
                    NeXLParticle.colorizedimage([bb], img),
                )
                FileIO.save(
                    File(format"PNG", dump * "[$cx, after, $split].png"),
                    NeXLParticle.colorizedimage(becomes, img),
                )
            end
            if split
                # Use each `becomes` only once
                # deleteat!(blobs, indexin(becomes, blobs))
                append!(newbest, becomes) # split it up...
            else
                push!(newbest, bb)
            end
        end
        best = newbest
    end
    return best
end

"""
   area(b::Blob)

Area of the Blob in pixel count.

area(b::Blob, img::AxisArray)

Area of the Blob in the units associated with `img`.  Assumes that each pixel in the image
is the same size.
"""
area(b::Blob) = count(b.mask)
function area(b::Blob, img::AxisArray)
    pa = prod( abs(last(a)-first(a))/length(a) for a in AxisArrays.axes(img) ) # area of one pixel
    return pa*area(b)
end

"""
    maskedimage(b::Blob, img::AbstractMatrix, mark=missing, markvalue=0.5)

Extract the image data in `img` associate the the Blob `b`.
"""
function maskedimage(b::Blob, img::AbstractMatrix)
    trimmed = copy(img[b.bounds])[:, :]
    foreach(idx -> trimmed[idx] = 0, filter(i -> !b.mask[i], eachindex(trimmed)))
    return trimmed
end

"""
    colorizedimage(bs::Vector{Blob}, img::AbstractArray)
    colorizedimage(chords::Vector{Vector{CartesianIndex{2}}}, img::AbstractArray)

Create a colorized version of img and draw the blob or chords on it.
"""
function colorizedimage(bs::Vector{Blob}, img::AbstractArray)
    off(ci, bs) = [ci.I[i] + bs.bounds.indices[i].start - 1 for i in eachindex(ci.I)]
    colors =
        convert.(
            RGB,
            distinguishable_colors(
                length(bs) + 2,
                Color[RGB(253 / 255, 255 / 255, 255 / 255), RGB(0, 0, 0), RGB(0 / 255, 168 / 255, 45 / 255)],
                transform = deuteranopic,
            )[3:end],
        )
    res = RGB.(img)
    for (i, blob) in enumerate(bs)
        c = colors[i]
        foreach(ci -> res[ci] = blob[ci] ? 0.5 * c + 0.5 * img[ci] : img[ci], CartesianIndices(blob)) # draw interior
        foreach(ci -> res[off(ci, blob)...] = c, perimeter(blob)) # draw perimeter
        res[off(blob.pstart, blob)...] = RGB(1.0, 0.0, 0.0) # draw start of perimeter...
    end
    return res
end

"""
    intersect(b1::Blob, b2::Blob)

A CartesianIndices with the region in common between b1 and b2.
"""
Base.intersect(b1::Blob, b2::Blob) = CartesianIndices(tuple(map(
    i -> intersect(b1.bounds.indices[i], b2.bounds.indices[i]),
    eachindex(b1.bounds.indices),
)...))

function interiorregions(b::Blob)
    onedge(b1, b2) = any(map(i -> b1[i].start == 1 || b1[i].stop == b2[i].stop - b2[i].start + 1, eachindex(b2)))
    tmp = filter(ib -> !onedge(ib.bounds.indices, b.bounds.indices), blob(b.mask, p -> !p))
    return Blob[Blob(__offset(tb.bounds, b.bounds), tb.mask) for tb in tmp]
end

"""
    filledarea(b::Blob)
    filledarea(b::Blob, img::AxisArray)

Area of the blob plus interior regions.
"""
filledarea(b::Blob) = area(b) + sum(area.(interiorregions(b)))
function filledarea(b::Blob, img::AxisArray)
    pa = prod( abs(last(a)-first(a))/length(a) for a in AxisArrays.axes(img) ) # area of one pixel
    return pa*filledarea(b)
end

function __offset(ci::CartesianIndices{N,R}, base::CartesianIndices{N,R}) where { N, R<:Tuple}
    rs = map(
        i -> ci.indices[i].start+base.indices[i].start-1:ci.indices[i].stop+base.indices[i].start-1,
        eachindex(ci.indices),
    )
    return CartesianIndices{N,R}(tuple(rs...))
end

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
# Automatically generated using Clang.jl

@cenum mfb_mouse_button::UInt32 begin
    MOUSE_BTN_0 = 0
    MOUSE_BTN_1 = 1
    MOUSE_BTN_2 = 2
    MOUSE_BTN_3 = 3
    MOUSE_BTN_4 = 4
    MOUSE_BTN_5 = 5
    MOUSE_BTN_6 = 6
    MOUSE_BTN_7 = 7
end

const MOUSE_LEFT = MOUSE_BTN_1
const MOUSE_RIGHT = MOUSE_BTN_2
const MOUSE_MIDDLE = MOUSE_BTN_3

@cenum mfb_key::Int32 begin
    KB_KEY_UNKNOWN = -1
    KB_KEY_SPACE = 32
    KB_KEY_APOSTROPHE = 39
    KB_KEY_COMMA = 44
    KB_KEY_MINUS = 45
    KB_KEY_PERIOD = 46
    KB_KEY_SLASH = 47
    KB_KEY_0 = 48
    KB_KEY_1 = 49
    KB_KEY_2 = 50
    KB_KEY_3 = 51
    KB_KEY_4 = 52
    KB_KEY_5 = 53
    KB_KEY_6 = 54
    KB_KEY_7 = 55
    KB_KEY_8 = 56
    KB_KEY_9 = 57
    KB_KEY_SEMICOLON = 59
    KB_KEY_EQUAL = 61
    KB_KEY_A = 65
    KB_KEY_B = 66
    KB_KEY_C = 67
    KB_KEY_D = 68
    KB_KEY_E = 69
    KB_KEY_F = 70
    KB_KEY_G = 71
    KB_KEY_H = 72
    KB_KEY_I = 73
    KB_KEY_J = 74
    KB_KEY_K = 75
    KB_KEY_L = 76
    KB_KEY_M = 77
    KB_KEY_N = 78
    KB_KEY_O = 79
    KB_KEY_P = 80
    KB_KEY_Q = 81
    KB_KEY_R = 82
    KB_KEY_S = 83
    KB_KEY_T = 84
    KB_KEY_U = 85
    KB_KEY_V = 86
    KB_KEY_W = 87
    KB_KEY_X = 88
    KB_KEY_Y = 89
    KB_KEY_Z = 90
    KB_KEY_LEFT_BRACKET = 91
    KB_KEY_BACKSLASH = 92
    KB_KEY_RIGHT_BRACKET = 93
    KB_KEY_GRAVE_ACCENT = 96
    KB_KEY_WORLD_1 = 161
    KB_KEY_WORLD_2 = 162
    KB_KEY_ESCAPE = 256
    KB_KEY_ENTER = 257
    KB_KEY_TAB = 258
    KB_KEY_BACKSPACE = 259
    KB_KEY_INSERT = 260
    KB_KEY_DELETE = 261
    KB_KEY_RIGHT = 262
    KB_KEY_LEFT = 263
    KB_KEY_DOWN = 264
    KB_KEY_UP = 265
    KB_KEY_PAGE_UP = 266
    KB_KEY_PAGE_DOWN = 267
    KB_KEY_HOME = 268
    KB_KEY_END = 269
    KB_KEY_CAPS_LOCK = 280
    KB_KEY_SCROLL_LOCK = 281
    KB_KEY_NUM_LOCK = 282
    KB_KEY_PRINT_SCREEN = 283
    KB_KEY_PAUSE = 284
    KB_KEY_F1 = 290
    KB_KEY_F2 = 291
    KB_KEY_F3 = 292
    KB_KEY_F4 = 293
    KB_KEY_F5 = 294
    KB_KEY_F6 = 295
    KB_KEY_F7 = 296
    KB_KEY_F8 = 297
    KB_KEY_F9 = 298
    KB_KEY_F10 = 299
    KB_KEY_F11 = 300
    KB_KEY_F12 = 301
    KB_KEY_F13 = 302
    KB_KEY_F14 = 303
    KB_KEY_F15 = 304
    KB_KEY_F16 = 305
    KB_KEY_F17 = 306
    KB_KEY_F18 = 307
    KB_KEY_F19 = 308
    KB_KEY_F20 = 309
    KB_KEY_F21 = 310
    KB_KEY_F22 = 311
    KB_KEY_F23 = 312
    KB_KEY_F24 = 313
    KB_KEY_F25 = 314
    KB_KEY_KP_0 = 320
    KB_KEY_KP_1 = 321
    KB_KEY_KP_2 = 322
    KB_KEY_KP_3 = 323
    KB_KEY_KP_4 = 324
    KB_KEY_KP_5 = 325
    KB_KEY_KP_6 = 326
    KB_KEY_KP_7 = 327
    KB_KEY_KP_8 = 328
    KB_KEY_KP_9 = 329
    KB_KEY_KP_DECIMAL = 330
    KB_KEY_KP_DIVIDE = 331
    KB_KEY_KP_MULTIPLY = 332
    KB_KEY_KP_SUBTRACT = 333
    KB_KEY_KP_ADD = 334
    KB_KEY_KP_ENTER = 335
    KB_KEY_KP_EQUAL = 336
    KB_KEY_LEFT_SHIFT = 340
    KB_KEY_LEFT_CONTROL = 341
    KB_KEY_LEFT_ALT = 342
    KB_KEY_LEFT_SUPER = 343
    KB_KEY_RIGHT_SHIFT = 344
    KB_KEY_RIGHT_CONTROL = 345
    KB_KEY_RIGHT_ALT = 346
    KB_KEY_RIGHT_SUPER = 347
    KB_KEY_MENU = 348
end

const KB_KEY_LAST = KB_KEY_MENU

@cenum mfb_update_state::Int32 begin
    STATE_OK = 0
    STATE_EXIT = -1
    STATE_INVALID_WINDOW = -2
    STATE_INVALID_BUFFER = -3
    STATE_INTERNAL_ERROR = -4
end

@cenum mfb_key_mod::UInt32 begin
    KB_MOD_SHIFT = 1
    KB_MOD_CONTROL = 2
    KB_MOD_ALT = 4
    KB_MOD_SUPER = 8
    KB_MOD_CAPS_LOCK = 16
    KB_MOD_NUM_LOCK = 32
end

@cenum mfb_window_flags::UInt32 begin
    WF_RESIZABLE = 1
    WF_FULLSCREEN = 2
    WF_FULLSCREEN_DESKTOP = 4
    WF_BORDERLESS = 8
    WF_ALWAYS_ON_TOP = 16
end

const mfb_window = Cvoid
const mfb_timer = Cvoid
const mfb_active_func = Ptr{Cvoid}
const mfb_resize_func = Ptr{Cvoid}
const mfb_keyboard_func = Ptr{Cvoid}
const mfb_char_input_func = Ptr{Cvoid}
const mfb_mouse_button_func = Ptr{Cvoid}
const mfb_mouse_move_func = Ptr{Cvoid}
const mfb_mouse_scroll_func = Ptr{Cvoid}

# Skipping MacroDefinition: MFB_RGB ( r , g , b ) ( ( ( uint32_t ) r ) << 16 ) | ( ( ( uint32_t ) g ) << 8 ) | ( ( uint32_t ) b )

ext: jl


text:
using TensorOperations

#
#	sanity checks
#
# Not tested yet
function sanityLeftCan(M::AbstractArray{<:Number}; atol = 10e-10)

	@tensor tmp[β1, β2] := M[α, d, β1]*conj(M)[α, d, β2]

	return isapprox(tmp, Id(size(tmp)[1]), atol)
end

#
#	Can forms
#

"""
    leftCanSite(M; optSVD = false)

brings M into left canonical form. Either by QR-decomposition or iff optSVD = true
bs signular-value decomposition

# Arguments
-

"""
function leftCanSite(M::AbstractArray{<:Number}; optSVD = false)
	α_dim, d_dim, β_dim = size(M)
	M = reshape(M, *(size(M)[1:2]...), size(M)[3])

	if optSVD == true
		U,S,V = svd(M)
		AL = reshape(U, α_dim, d_dim, size(S)[1])

		return AL, diagm(0 => S), V'

	else
		Q, R = qr!(M)
		γ_dim = min(size(M)...)
		Q = Q*Matrix(I, size(Q)[2], γ_dim)
		AL = reshape(Q, α_dim, d_dim, γ_dim)

		return AL, R
	end
end

"""
    leftCanMPS(MPSvec; Cpre = nothing)

takes Cpre and brings the MPS chain Cpre - M.... M - stored in Mvec into left-canonical form.

# Arguments
- MPSvec: Vector cotaining the M's of the MPS which have to be transformed into the left-canonical form
- Cpre: singular values which come before the first element in MPSvec

return:
	Cnext, Remainder, Cvec

	perfoms SVD on last site so M = U, S, V were the last M in left can. form is formed from U
	S is returned as Cnext, and the V is returned as the Remainder.
	For all intermidiate steps S of the svd decomposition has been added to a vector Cvec.
	Note the last Cnext has not been added to Cvec by default! (Cvec just contains all intermidiate Cpre of each A)

"""
function leftCanMPS(MPSvec::SubArray{<:Any}; Cpre::Union{AbstractArray{<:Number}, Nothing} = nothing)
	if Cpre != nothing
		@tensor MPSvec[1][α, d, β] := Cpre[α, γ]*MPSvec[1][γ, d, β]
	end
	Cvec = []

	for site in 1:size(MPSvec[1:end-1])[1]
		F = leftCanSite(MPSvec[site]; optSVD = true)
		MPSvec[site] = F[1]
		@tensor newMPSnext[α, d, β] := F[2][α, γ]*F[3][γ, γ̃]*MPSvec[site+1][γ̃, d, β]
		MPSvec[site+1] = newMPSnext
		push!(Cvec, F[2])
	end

	F = leftCanSite(MPSvec[end]; optSVD = true)
	MPSvec[end] = F[1]
	Cnext = F[2]
	Remainder = F[3]

	return Cnext, Remainder, Cvec
end

"""
    rightCanSite(M; optSVD = false)

brings M into right canonical form. Either by QR-decomposition or iff optSVD = true
bs signular-value decomposition

# Arguments
-

"""
function rightCanSite(M::AbstractArray{<:Number}; optSVD = false)
	α_dim, d_dim, β_dim = size(M)
	M = reshape(M, size(M)[1], *(size(M)[2:3]...))

	if optSVD == true
		U, S, V = svd(M)
		AR = reshape(V', size(S)[1], d_dim, β_dim)

		return AR, diagm( 0 => S), U
	else
		L, Q = lq!(M)

		γ_dim = min(size(M)...)

		Q = Matrix(I, γ_dim, size(Q)[1])*Q

		AR = reshape(Q, γ_dim, d_dim, β_dim)

		return AR, L
	end
end

"""
    rightCanSite(M; optSVD = false)

brings M into right canonical form. Either by QR-decomposition or iff optSVD = true
bs signular-value decomposition

# Arguments
-

"""
function rightCanMPS(MPSvec::SubArray{<:Any}; Cnext::Union{AbstractArray{<:Number}, Nothing} = nothing)
	if Cnext != nothing

		@tensor MPSvec[end][α, d, β] = Cnext[γ, β] * MPSvec[end][α, d, γ]
	end

	Cvec = []
	for site in size(MPSvec)[1]:-1:2
		F = rightCanSite(MPSvec[site]; optSVD = true)
		MPSvec[site] = F[1]

		@tensor newMPSnext[α, d, β] := MPSvec[site-1][α, d, γ̃] * F[3][γ̃, γ]*F[2][γ, β]
		MPSvec[site-1] = newMPSnext
		push!(Cvec, F[2])
	end

	F = rightCanSite(MPSvec[1]; optSVD = true)

	MPSvec[1] = F[1]
	Cpre = F[2]
	Remainder = F[3]
	#push!(Cvec, Cpre)
	return Cpre, Remainder, Cvec[end:-1:1]
end

rightCanMPS(MPSvec::Vector{<:Any}; Cnext::Union{AbstractArray{<:Number}, Nothing} = nothing,
	    from::Int = 1, to::Int = size(MPSvec)[1]) = rightCanMPS(@view MPSvec[from:to]; Cnext = Cnext)
leftCanMPS(MPSvec::Vector{<:Any}; Cpre::Union{AbstractArray{<:Number}, Nothing} = nothing,
	   from::Int = 1, to::Int = size(MPSvec)[1]) = leftCanMPS(@view MPSvec[from:to]; Cpre = Cpre)
rightCanMPS(MPSvec::Vector{<:Any}, i::Int, j::Int) = rightCanMPS(MPSvec; Cnext = nothing, from = i, to = j)
leftCanMPS(MPSvec::Vector{<:Any}, i::Int, j::Int) = leftCanMPS(MPSvec; Cpre = nothing, from = i, to = j)

"""
    mixedCanMPS(MPSvec, site)

brings the MPSvec into mixed canonical from. That means all sites 1-site are in left can. form
while all sites from site+1-end are in right can. form.

# Arguments
-

return:

	Schmidt-Coefficients between left and right part.

"""
function mixedCanMPS(MPSvec::Vector{<:Any}, site::Int)

	MPS_size = size(MPSvec)[1]
	Cpre, Rr, Cvec = rightCanMPS(MPSvec, site+1, MPS_size)
	Cnext, Rl, Cvec = leftCanMPS(MPSvec, 1, site)

	F = svd(Cnext*Rl*Rr*Cpre)

	# multiply F.U and F.V into left respectively right MPS block
	@tensor res[α, d, β_new] := MPSvec[site][α, d, β]*F.U[β, β_new]
	MPSvec[site] = res

	@tensor res[α_new, d, β] := MPSvec[site+1][α, d, β] * F.V'[α_new, α]
	MPSvec[site+1] = res

	return F.S
end

"""
    getSchmidtVec(MPSvec, site)

# Arguments
-

return:

	Schmidt-Coefficients between left and right part.

"""
function getSchmidtVec(MPSvec::Vector{<:Any}, site::Int)
	S = mixedCanMPS(MPSvec, site)

	leftHalf = @view MPSvec[1:site]
	rightHalf = @view MPSvec[site+1:end]

	Svec_left = mergeMPS(leftHalf) #(phys_dim, ingoing, outgoing)
	Svec_left = Svec_left[:,1,:]
	Svec_right = mergeMPS(rightHalf)
	Svec_right = Svec_right[:,:,1]

	return S, Svec_left, Svec_right
end

max_stars_repo_name: Gawatz/finiteMPS

max_issues_repo_name: Gawatz/finiteMPS

max_forks_repo_name: Gawatz/finiteMPS


text:
import PyPlot
const plt = PyPlot

include("lagrange_interp.jl")

function myfunc(x)
    return log(x) + 1
end

function main()
    x = [1.0, 2.0, 3.0, 4.0, 5.0]
    y = myfunc.(x)

    NptsPlot = 100
    xo = collect(range(x[begin], stop=x[end], length=NptsPlot))
    yo = zeros(NptsPlot)
    for i in 1:NptsPlot
        @views yo[i] = lagrange_interp(x[1:4], y[1:4], xo[i])
    end

    plt.clf()
    plt.plot(x, y, marker="o", label="data")
    plt.plot(xo, yo, label="interp")
    plt.grid(true)
    plt.legend()
    plt.savefig("IMG_test_lagrange_interp.png", dpi=150)
end

@time main()
@time main()


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

abstract type AbstractQuasiGradientExecution end
# Execution API #
# ------------------------------------------------------------
restore_subproblems!(quasigradient::AbstractQuasiGradient) = restore_subproblems!(quasigradient, quasigradient.execution)
solve_master!(quasigradient::AbstractQuasiGradient) = solve_master!(quasigradient, quasigradient.execution)
iterate!(quasigradient::AbstractQuasiGradient) = iterate!(quasigradient, quasigradient.execution)
start_workers!(quasigradient::AbstractQuasiGradient) = start_workers!(quasigradient, quasigradient.execution)
close_workers!(quasigradient::AbstractQuasiGradient) = close_workers!(quasigradient, quasigradient.execution)
resolve_subproblems!(quasigradient::AbstractQuasiGradient) = resolve_subproblems!(quasigradient, quasigradient.execution)
current_decision(quasigradient::AbstractQuasiGradient) = current_decision(quasigradient, quasigradient.execution)
calculate_objective_value(quasigradient::AbstractQuasiGradient) = calculate_objective_value(quasigradient, quasigradient.execution)
timestamp(quasigradient::AbstractQuasiGradient) = timestamp(quasigradient, quasigradient.execution)
subobjectives(quasigradient::AbstractQuasiGradient) = subobjectives(quasigradient, quasigradient.execution)
set_subobjectives(quasigradient::AbstractQuasiGradient, Qs::AbstractVector) = set_subobjectives(quasigradient, Qs, quasigradient.execution)
# ------------------------------------------------------------
include("common.jl")
include("serial.jl")
include("synchronous.jl")
#include("asynchronous.jl")


text:
#=
    Elec grid modelling
 =#

mutable struct Grid <: AbstractGrid
    # Parameters
    powerMax::Float64
    # Variables
    carrier::EnergyCarrier
    cost_in::AbstractArray{Float64,3}
    cost_out::AbstractArray{Float64,3}
    # Inner constructor
    Grid(; powerMax = 36., carrier = Electricity()) = new(powerMax, carrier)
end

### Preallocation
function preallocate!(grid::Grid, nh::Int64, ny::Int64, ns::Int64)
     grid.carrier.power = convert(SharedArray,zeros(nh, ny, ns))
     grid.cost_in = convert(SharedArray,zeros(nh, ny, ns))
     grid.cost_out = convert(SharedArray,zeros(nh, ny, ns))
     return grid
end


text:
import Mongoc
using Dates
using NamedTupleTools
using DataStructures

struct TypeIsNotAMapping <: Exception
    msg::AbstractString
end

"""
Create a BSON document from an arbitrary value.
Recursively converts the internal structure of arrays, tuples, structs, etc.
"""
function Mongoc.BSON(value::Any)
    fnames = fieldnames(typeof(value))
    if isempty(fnames)
        throw(TypeIsNotAMapping("cannot create BSON from DataType $(typeof(value)), data must be a mapping with keys or fields accessible via the keys() or fieldnames() functions"))
    end

    doc = Mongoc.BSON()
    for (fname, ftype) in zip(fnames, typeof(value).types)
        doc[string(fname)] = getfield(value, fname)
    end
    return doc
end

# Mongoc only supports Float64, add other float types that fit into float64.
Base.setindex!(document::Mongoc.BSON, value::Union{Float16, Float32}, key::AbstractString) = setindex!(document, Float64(value), key)

# Mongoc only supports Int32 and Int64, and smaller signed int types.
Base.setindex!(document::Mongoc.BSON, value::Union{Int8, Int16}, key::AbstractString) = setindex!(document, Int64(value), key)

# Mongoc supports no unsigned integer types, add the ones which will fit into (signed) Int32/Int64.
Base.setindex!(document::Mongoc.BSON, value::Union{UInt8, UInt16}, key::AbstractString) = setindex!(document, Int32(value), key)
Base.setindex!(document::Mongoc.BSON, value::UInt32, key::AbstractString) = setindex!(document, Int64(value), key)

# Mongoc does not support using numbers as Dict keys.
Base.setindex!(document::Mongoc.BSON, value::Any, key::Number) = setindex!(document, value, string(key))

Base.setindex!(document::Mongoc.BSON, value::Tuple, key::AbstractString) = setindex!(document, collect(value), key)

function Base.setindex!(document::Mongoc.BSON, value::NamedTuple, key::AbstractString)
    tupleKeys = map(String, keys(value))
    tupleValues = values(value)
    data = Mongoc.BSON()
    for (k, v) in zip(tupleKeys, tupleValues)
        data[k] = v
    end
    setindex!(document, data, key)
end

Base.setindex!(document::Mongoc.BSON, value::DataType, key::AbstractString) = setindex!(document, string(value), key)

Base.setindex!(document::Mongoc.BSON, value::Any, key::AbstractString) = setindex!(document, Mongoc.BSON(value), key)

max_stars_repo_name: OliverColeman/BSONify

max_forks_repo_name: OliverColeman/BSONify


text:
using Revise
using RadiativeTransfer
using RadiativeTransfer.Absorption
using RadiativeTransfer.Scattering
using RadiativeTransfer.vSmartMOM
using Parameters
using LinearAlgebra

# Sets all the "specific" parameters
parameters = vSmartMOM.default_parameters();

# Generates all the derived attributes from above parameters
model = default_model(parameters);

pol_type = model.params.polarization_type;           # Polarization type (IQUV)
obs_geom = model.obs_geom::ObsGeometry;  #, # Solar Zenith, Viewing Zenith, Viewing Azimuthal
τRayl = model.τRayl;        # Rayleigh optical depth
    #nAer,                 # Number of aerosol species
τAer =  model.τAer ;                # Aerosol optical depth and single-scattering albedo
quadPoints = model.quad_points;
#wt_μ =  model.wt_μ  ;      # Quadrature points and weights
Ltrunc =  3  ;            # Trunction length for legendre terms
aerosol_optics = model.aerosol_optics;       # AerosolOptics (greek_coefs, ω̃, k, fᵗ)
GreekRayleigh =  model.greek_rayleigh ;       # Greek coefficients of Rayleigh Phase Function
τ_abs =    model.τ_abs  ;           # nSpec x N
architecture = model.params.architecture;

@unpack obs_alt, sza, vza, vaz = obs_geom   # Observational geometry properties

FT = eltype(sza)                    # Get the float-type to use
Nz = length(τRayl)                  # Number of vertical slices
nSpec = size(τ_abs, 1)              # Number of spectral points
##μ0 = cosd(sza)                      # μ0 defined as cos(θ); θ = sza
#iμ0 = vSmartMOM.nearest_point(qp_μ, μ0)       # Find the closest point to μ0 in qp_μ
arr_type = array_type(architecture)

# Output variables: Reflected and transmitted solar irradiation at TOA and BOA respectively
R = zeros(FT, length(vza), pol_type.n, nSpec)
T = zeros(FT, length(vza), pol_type.n, nSpec)
R_SFI = zeros(FT, length(vza), pol_type.n, nSpec)
T_SFI = zeros(FT, length(vza), pol_type.n, nSpec)

# Copy qp_μ "pol_type.n" times
#qp_μN = arr_type(reshape(transpose(repeat(qp_μ, 1, pol_type.n)),pol_type.n*size(qp_μ)[1],1))
#wt_μN = arr_type(reshape(transpose(repeat(wt_μ, 1, pol_type.n)),pol_type.n*size(wt_μ)[1],1))
#for i = 1:length(qp_μN)
#   @show(i,qp_μN[i])
#end
println("Processing on: ", architecture)
println("With FT: ", FT)

#=
Loop over number of truncation terms =#
SFI = true

m = 0

println("Fourier Moment: ", m)

# Azimuthal weighting
weight = m == 0 ? FT(0.5) : FT(1.0)
# Compute Z-moments of the Rayleigh phase matrix
# For m>=3, Rayleigh matrices will be 0, can catch with if statement if wanted
Rayl𝐙⁺⁺, Rayl𝐙⁻⁺ = Scattering.compute_Z_moments(pol_type, Array(quadPoints.qp_μ), GreekRayleigh, m, arr_type = arr_type);
# Number of aerosols
#@show size(aerosol_optics)
#nBand = length(aerosol_optics)
nAer  = length(aerosol_optics)
# Just for now:
iBand = 1
#nAer, nBand = size(aerosol_optics)
@show nAer#, nBand
dims = size(Rayl𝐙⁺⁺)

# Compute aerosol Z-matrices for all aerosols
Aer𝐙⁺⁺ = arr_type(zeros(FT, (dims[1], dims[2], nAer)))
Aer𝐙⁻⁺ = similar(Aer𝐙⁺⁺)

for i = 1:nAer
    @show aerosol_optics[i,1]
    Aer𝐙⁺⁺[:,:,i], Aer𝐙⁻⁺[:,:,i] = Scattering.compute_Z_moments(pol_type, Array(quadPoints.qp_μ), aerosol_optics[i].greek_coefs, m, arr_type = arr_type)
end

# R and T matrices for Added and Composite Layers for this m
added_layer = vSmartMOM.make_added_layer(FT, arr_type, dims, nSpec)
composite_layer = vSmartMOM.make_composite_layer(FT, arr_type, dims, nSpec)
I_static = Diagonal(arr_type(Diagonal{FT}(ones(dims[1]))));
scattering_interface = vSmartMOM.ScatteringInterface_00()
τ_sum = arr_type(zeros(nSpec)) #Suniti: declaring τ_sum to be of length nSpec
τ_λ = arr_type(zeros(nSpec))
iz = Nz
if iz==1
    τ_sum = τ_λ
else
    τ_sum = τ_sum + τ_λ
end
#@show(iz, Nz)
# Construct the atmospheric layer
# From Rayleigh and aerosol τ, ϖ, compute overall layer τ, ϖ
#@timeit "Constructing"

τ_λ, ϖ_λ, τ, ϖ, Z⁺⁺, Z⁻⁺ = vSmartMOM.construct_atm_layer(τRayl[iz], τAer[:,iz], aerosol_optics, Rayl𝐙⁺⁺, Rayl𝐙⁻⁺, Aer𝐙⁺⁺, Aer𝐙⁻⁺, τ_abs[:,iz], arr_type)
#@show(τ_λ)
#@show(ϖ_λ)
#@show(τ)
#@show(ϖ)
#sleep(5)
#for i=1:size(Z⁺⁺)[1]
#    @show(i,Z⁺⁺[i,:])
#end
# τ * ϖ should remain constant even though they individually change over wavelength
# @assert all(i -> (i ≈ τ * ϖ), τ_λ .* ϖ_λ)
# Compute doubling number
dτ_max = minimum([τ * ϖ, FT(0.01) * minimum(quadPoints.qp_μ)])
dτ, ndoubl = vSmartMOM.doubling_number(dτ_max, τ * ϖ) #Suniti
#@show(ndoubl, dτ_max, τ)
# Compute dτ vector
dτ_λ = arr_type(τ_λ ./ (FT(2)^ndoubl))
expk = exp.(-dτ_λ /quadPoints.μ₀) #Suniti
#@show(τ_λ, dτ_λ.*FT(2)^ndoubl)
#@show(τ, dτ*FT(2)^ndoubl,dτ, dτ_λ )
#@show(expk, exp.(-dτ /qp_μ[iμ0]))
#@show τ_sum
#@show dτ_λ, dτ
#scatter = true
vSmartMOM.elemental!(pol_type, SFI,  τ_sum, dτ_λ, dτ, ϖ_λ, ϖ, Z⁺⁺, Z⁻⁺, m, ndoubl, true, quadPoints, added_layer,  I_static, architecture)
vSmartMOM.doubling!(pol_type, SFI, expk, ndoubl, added_layer, I_static, architecture)

lambipoo = vSmartMOM.LambertianSurfaceScalar(0.1)
vSmartMOM.create_surface_layer!(lambipoo, added_layer, SFI, m, pol_type, quadPoints, τ_sum, architecture);
#added_layer_DNI = vSmartMOM.make_added_layer(FT, arr_type, dims, nSpec)
#composite_layer_DNI = vSmartMOM.make_composite_layer(FT, arr_type, dims, nSpec)
#vSmartMOM.elemental!(pol_type, false, iμ0, τ_sum, dτ_λ, dτ, ϖ_λ, ϖ, Z⁺⁺, Z⁻⁺, m, ndoubl, true, qp_μ, wt_μ, added_layer_DNI,  I_static, arr_type, architecture)


text:
# Julia wrapper for header: kinsol.h
# Automatically generated using Clang.jl

function KINCreate()
    ccall((:KINCreate, libsundials_kinsol), KINMemPtr, ())
end

function KINInit(kinmem, func::KINSysFn, tmpl::N_Vector)
    ccall((:KINInit, libsundials_kinsol), Cint, (KINMemPtr, KINSysFn, N_Vector), kinmem, func, tmpl)
end

function KINInit(kinmem, func, tmpl)
    __tmpl = convert(NVector, tmpl)
    KINInit(kinmem, KINSysFn_wrapper(func), convert(N_Vector, __tmpl))
end

function KINSol(kinmem, uu::N_Vector, strategy::Cint, u_scale::N_Vector, f_scale::N_Vector)
    ccall((:KINSol, libsundials_kinsol), Cint, (KINMemPtr, N_Vector, Cint, N_Vector, N_Vector), kinmem, uu, strategy, u_scale, f_scale)
end

function KINSol(kinmem, uu, strategy, u_scale, f_scale)
    __uu = convert(NVector, uu)
    __u_scale = convert(NVector, u_scale)
    __f_scale = convert(NVector, f_scale)
    KINSol(kinmem, convert(N_Vector, __uu), convert(Cint, strategy), convert(N_Vector, __u_scale), convert(N_Vector, __f_scale))
end

function KINSetErrHandlerFn(kinmem, ehfun::KINErrHandlerFn, eh_data)
    ccall((:KINSetErrHandlerFn, libsundials_kinsol), Cint, (KINMemPtr, KINErrHandlerFn, Ptr{Cvoid}), kinmem, ehfun, eh_data)
end

function KINSetErrFile(kinmem, errfp)
    ccall((:KINSetErrFile, libsundials_kinsol), Cint, (KINMemPtr, Ptr{FILE}), kinmem, errfp)
end

function KINSetInfoHandlerFn(kinmem, ihfun::KINInfoHandlerFn, ih_data)
    ccall((:KINSetInfoHandlerFn, libsundials_kinsol), Cint, (KINMemPtr, KINInfoHandlerFn, Ptr{Cvoid}), kinmem, ihfun, ih_data)
end

function KINSetInfoFile(kinmem, infofp)
    ccall((:KINSetInfoFile, libsundials_kinsol), Cint, (KINMemPtr, Ptr{FILE}), kinmem, infofp)
end

function KINSetUserData(kinmem, user_data)
    ccall((:KINSetUserData, libsundials_kinsol), Cint, (KINMemPtr, Any), kinmem, user_data)
end

function KINSetPrintLevel(kinmemm, printfl::Cint)
    ccall((:KINSetPrintLevel, libsundials_kinsol), Cint, (KINMemPtr, Cint), kinmemm, printfl)
end

function KINSetPrintLevel(kinmemm, printfl)
    KINSetPrintLevel(kinmemm, convert(Cint, printfl))
end

function KINSetMAA(kinmem, maa::Clong)
    ccall((:KINSetMAA, libsundials_kinsol), Cint, (KINMemPtr, Clong), kinmem, maa)
end

function KINSetMAA(kinmem, maa)
    KINSetMAA(kinmem, convert(Clong, maa))
end

function KINSetDampingAA(kinmem, beta::realtype)
    ccall((:KINSetDampingAA, libsundials_kinsol), Cint, (KINMemPtr, realtype), kinmem, beta)
end

function KINSetNumMaxIters(kinmem, mxiter::Clong)
    ccall((:KINSetNumMaxIters, libsundials_kinsol), Cint, (KINMemPtr, Clong), kinmem, mxiter)
end

function KINSetNumMaxIters(kinmem, mxiter)
    KINSetNumMaxIters(kinmem, convert(Clong, mxiter))
end

function KINSetNoInitSetup(kinmem, noInitSetup::Cint)
    ccall((:KINSetNoInitSetup, libsundials_kinsol), Cint, (KINMemPtr, Cint), kinmem, noInitSetup)
end

function KINSetNoInitSetup(kinmem, noInitSetup)
    KINSetNoInitSetup(kinmem, convert(Cint, noInitSetup))
end

function KINSetNoResMon(kinmem, noNNIResMon::Cint)
    ccall((:KINSetNoResMon, libsundials_kinsol), Cint, (KINMemPtr, Cint), kinmem, noNNIResMon)
end

function KINSetNoResMon(kinmem, noNNIResMon)
    KINSetNoResMon(kinmem, convert(Cint, noNNIResMon))
end

function KINSetMaxSetupCalls(kinmem, msbset::Clong)
    ccall((:KINSetMaxSetupCalls, libsundials_kinsol), Cint, (KINMemPtr, Clong), kinmem, msbset)
end

function KINSetMaxSetupCalls(kinmem, msbset)
    KINSetMaxSetupCalls(kinmem, convert(Clong, msbset))
end

function KINSetMaxSubSetupCalls(kinmem, msbsetsub::Clong)
    ccall((:KINSetMaxSubSetupCalls, libsundials_kinsol), Cint, (KINMemPtr, Clong), kinmem, msbsetsub)
end

function KINSetMaxSubSetupCalls(kinmem, msbsetsub)
    KINSetMaxSubSetupCalls(kinmem, convert(Clong, msbsetsub))
end

function KINSetEtaForm(kinmem, etachoice::Cint)
    ccall((:KINSetEtaForm, libsundials_kinsol), Cint, (KINMemPtr, Cint), kinmem, etachoice)
end

function KINSetEtaForm(kinmem, etachoice)
    KINSetEtaForm(kinmem, convert(Cint, etachoice))
end

function KINSetEtaConstValue(kinmem, eta::realtype)
    ccall((:KINSetEtaConstValue, libsundials_kinsol), Cint, (KINMemPtr, realtype), kinmem, eta)
end

function KINSetEtaParams(kinmem, egamma::realtype, ealpha::realtype)
    ccall((:KINSetEtaParams, libsundials_kinsol), Cint, (KINMemPtr, realtype, realtype), kinmem, egamma, ealpha)
end

function KINSetResMonParams(kinmem, omegamin::realtype, omegamax::realtype)
    ccall((:KINSetResMonParams, libsundials_kinsol), Cint, (KINMemPtr, realtype, realtype), kinmem, omegamin, omegamax)
end

function KINSetResMonConstValue(kinmem, omegaconst::realtype)
    ccall((:KINSetResMonConstValue, libsundials_kinsol), Cint, (KINMemPtr, realtype), kinmem, omegaconst)
end

function KINSetNoMinEps(kinmem, noMinEps::Cint)
    ccall((:KINSetNoMinEps, libsundials_kinsol), Cint, (KINMemPtr, Cint), kinmem, noMinEps)
end

function KINSetNoMinEps(kinmem, noMinEps)
    KINSetNoMinEps(kinmem, convert(Cint, noMinEps))
end

function KINSetMaxNewtonStep(kinmem, mxnewtstep::realtype)
    ccall((:KINSetMaxNewtonStep, libsundials_kinsol), Cint, (KINMemPtr, realtype), kinmem, mxnewtstep)
end

function KINSetMaxBetaFails(kinmem, mxnbcf::Clong)
    ccall((:KINSetMaxBetaFails, libsundials_kinsol), Cint, (KINMemPtr, Clong), kinmem, mxnbcf)
end

function KINSetMaxBetaFails(kinmem, mxnbcf)
    KINSetMaxBetaFails(kinmem, convert(Clong, mxnbcf))
end

function KINSetRelErrFunc(kinmem, relfunc::realtype)
    ccall((:KINSetRelErrFunc, libsundials_kinsol), Cint, (KINMemPtr, realtype), kinmem, relfunc)
end

function KINSetFuncNormTol(kinmem, fnormtol::realtype)
    ccall((:KINSetFuncNormTol, libsundials_kinsol), Cint, (KINMemPtr, realtype), kinmem, fnormtol)
end

function KINSetScaledStepTol(kinmem, scsteptol::realtype)
    ccall((:KINSetScaledStepTol, libsundials_kinsol), Cint, (KINMemPtr, realtype), kinmem, scsteptol)
end

function KINSetConstraints(kinmem, constraints::N_Vector)
    ccall((:KINSetConstraints, libsundials_kinsol), Cint, (KINMemPtr, N_Vector), kinmem, constraints)
end

function KINSetConstraints(kinmem, constraints)
    __constraints = convert(NVector, constraints)
    KINSetConstraints(kinmem, convert(N_Vector, __constraints))
end

function KINSetSysFunc(kinmem, func::KINSysFn)
    ccall((:KINSetSysFunc, libsundials_kinsol), Cint, (KINMemPtr, KINSysFn), kinmem, func)
end

function KINGetWorkSpace(kinmem, lenrw, leniw)
    ccall((:KINGetWorkSpace, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}, Ptr{Clong}), kinmem, lenrw, leniw)
end

function KINGetNumNonlinSolvIters(kinmem, nniters)
    ccall((:KINGetNumNonlinSolvIters, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, nniters)
end

function KINGetNumFuncEvals(kinmem, nfevals)
    ccall((:KINGetNumFuncEvals, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, nfevals)
end

function KINGetNumBetaCondFails(kinmem, nbcfails)
    ccall((:KINGetNumBetaCondFails, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, nbcfails)
end

function KINGetNumBacktrackOps(kinmem, nbacktr)
    ccall((:KINGetNumBacktrackOps, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, nbacktr)
end

function KINGetFuncNorm(kinmem, fnorm)
    ccall((:KINGetFuncNorm, libsundials_kinsol), Cint, (KINMemPtr, Ptr{realtype}), kinmem, fnorm)
end

function KINGetStepLength(kinmem, steplength)
    ccall((:KINGetStepLength, libsundials_kinsol), Cint, (KINMemPtr, Ptr{realtype}), kinmem, steplength)
end

function KINGetReturnFlagName(flag::Clong)
    ccall((:KINGetReturnFlagName, libsundials_kinsol), Cstring, (Clong,), flag)
end

function KINGetReturnFlagName(flag)
    KINGetReturnFlagName(convert(Clong, flag))
end

function KINFree(kinmem)
    ccall((:KINFree, libsundials_kinsol), Cvoid, (Ref{KINMemPtr},), kinmem)
end
# Julia wrapper for header: kinsol_bbdpre.h
# Automatically generated using Clang.jl

function KINBBDPrecInit(kinmem, Nlocal::sunindextype, mudq::sunindextype, mldq::sunindextype, mukeep::sunindextype, mlkeep::sunindextype, dq_rel_uu::realtype, gloc::KINBBDLocalFn, gcomm::KINBBDCommFn)
    ccall((:KINBBDPrecInit, libsundials_kinsol), Cint, (KINMemPtr, sunindextype, sunindextype, sunindextype, sunindextype, sunindextype, realtype, KINBBDLocalFn, KINBBDCommFn), kinmem, Nlocal, mudq, mldq, mukeep, mlkeep, dq_rel_uu, gloc, gcomm)
end

function KINBBDPrecGetWorkSpace(kinmem, lenrwBBDP, leniwBBDP)
    ccall((:KINBBDPrecGetWorkSpace, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}, Ptr{Clong}), kinmem, lenrwBBDP, leniwBBDP)
end

function KINBBDPrecGetNumGfnEvals(kinmem, ngevalsBBDP)
    ccall((:KINBBDPrecGetNumGfnEvals, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, ngevalsBBDP)
end
# Julia wrapper for header: kinsol_direct.h
# Automatically generated using Clang.jl

function KINDlsSetLinearSolver(kinmem, LS::SUNLinearSolver, A::SUNMatrix)
    ccall((:KINDlsSetLinearSolver, libsundials_kinsol), Cint, (KINMemPtr, SUNLinearSolver, SUNMatrix), kinmem, LS, A)
end

function KINDlsSetJacFn(kinmem, jac::KINDlsJacFn)
    ccall((:KINDlsSetJacFn, libsundials_kinsol), Cint, (KINMemPtr, KINDlsJacFn), kinmem, jac)
end

function KINDlsGetWorkSpace(kinmem, lenrw, leniw)
    ccall((:KINDlsGetWorkSpace, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}, Ptr{Clong}), kinmem, lenrw, leniw)
end

function KINDlsGetNumJacEvals(kinmem, njevals)
    ccall((:KINDlsGetNumJacEvals, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, njevals)
end

function KINDlsGetNumFuncEvals(kinmem, nfevals)
    ccall((:KINDlsGetNumFuncEvals, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, nfevals)
end

function KINDlsGetLastFlag(kinmem, flag)
    ccall((:KINDlsGetLastFlag, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, flag)
end

function KINDlsGetReturnFlagName(flag::Clong)
    ccall((:KINDlsGetReturnFlagName, libsundials_kinsol), Cstring, (Clong,), flag)
end

function KINDlsGetReturnFlagName(flag)
    KINDlsGetReturnFlagName(convert(Clong, flag))
end
# Julia wrapper for header: kinsol_ls.h
# Automatically generated using Clang.jl

function KINSetLinearSolver(kinmem, LS::SUNLinearSolver, A::SUNMatrix)
    ccall((:KINSetLinearSolver, libsundials_kinsol), Cint, (KINMemPtr, SUNLinearSolver, SUNMatrix), kinmem, LS, A)
end

function KINSetJacFn(kinmem, jac::KINLsJacFn)
    ccall((:KINSetJacFn, libsundials_kinsol), Cint, (KINMemPtr, KINLsJacFn), kinmem, jac)
end

function KINSetPreconditioner(kinmem, psetup::KINLsPrecSetupFn, psolve::KINLsPrecSolveFn)
    ccall((:KINSetPreconditioner, libsundials_kinsol), Cint, (KINMemPtr, KINLsPrecSetupFn, KINLsPrecSolveFn), kinmem, psetup, psolve)
end

function KINSetJacTimesVecFn(kinmem, jtv::KINLsJacTimesVecFn)
    ccall((:KINSetJacTimesVecFn, libsundials_kinsol), Cint, (KINMemPtr, KINLsJacTimesVecFn), kinmem, jtv)
end

function KINGetLinWorkSpace(kinmem, lenrwLS, leniwLS)
    ccall((:KINGetLinWorkSpace, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}, Ptr{Clong}), kinmem, lenrwLS, leniwLS)
end

function KINGetNumJacEvals(kinmem, njevals)
    ccall((:KINGetNumJacEvals, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, njevals)
end

function KINGetNumLinFuncEvals(kinmem, nfevals)
    ccall((:KINGetNumLinFuncEvals, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, nfevals)
end

function KINGetNumPrecEvals(kinmem, npevals)
    ccall((:KINGetNumPrecEvals, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, npevals)
end

function KINGetNumPrecSolves(kinmem, npsolves)
    ccall((:KINGetNumPrecSolves, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, npsolves)
end

function KINGetNumLinIters(kinmem, nliters)
    ccall((:KINGetNumLinIters, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, nliters)
end

function KINGetNumLinConvFails(kinmem, nlcfails)
    ccall((:KINGetNumLinConvFails, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, nlcfails)
end

function KINGetNumJtimesEvals(kinmem, njvevals)
    ccall((:KINGetNumJtimesEvals, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, njvevals)
end

function KINGetLastLinFlag(kinmem, flag)
    ccall((:KINGetLastLinFlag, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, flag)
end

function KINGetLinReturnFlagName(flag::Clong)
    ccall((:KINGetLinReturnFlagName, libsundials_kinsol), Cstring, (Clong,), flag)
end

function KINGetLinReturnFlagName(flag)
    KINGetLinReturnFlagName(convert(Clong, flag))
end
# Julia wrapper for header: kinsol_spils.h
# Automatically generated using Clang.jl

function KINSpilsSetLinearSolver(kinmem, LS::SUNLinearSolver)
    ccall((:KINSpilsSetLinearSolver, libsundials_kinsol), Cint, (KINMemPtr, SUNLinearSolver), kinmem, LS)
end

function KINSpilsSetPreconditioner(kinmem, psetup::KINSpilsPrecSetupFn, psolve::KINSpilsPrecSolveFn)
    ccall((:KINSpilsSetPreconditioner, libsundials_kinsol), Cint, (KINMemPtr, KINSpilsPrecSetupFn, KINSpilsPrecSolveFn), kinmem, psetup, psolve)
end

function KINSpilsSetJacTimesVecFn(kinmem, jtv::KINSpilsJacTimesVecFn)
    ccall((:KINSpilsSetJacTimesVecFn, libsundials_kinsol), Cint, (KINMemPtr, KINSpilsJacTimesVecFn), kinmem, jtv)
end

function KINSpilsGetWorkSpace(kinmem, lenrwLS, leniwLS)
    ccall((:KINSpilsGetWorkSpace, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}, Ptr{Clong}), kinmem, lenrwLS, leniwLS)
end

function KINSpilsGetNumPrecEvals(kinmem, npevals)
    ccall((:KINSpilsGetNumPrecEvals, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, npevals)
end

function KINSpilsGetNumPrecSolves(kinmem, npsolves)
    ccall((:KINSpilsGetNumPrecSolves, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, npsolves)
end

function KINSpilsGetNumLinIters(kinmem, nliters)
    ccall((:KINSpilsGetNumLinIters, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, nliters)
end

function KINSpilsGetNumConvFails(kinmem, nlcfails)
    ccall((:KINSpilsGetNumConvFails, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, nlcfails)
end

function KINSpilsGetNumJtimesEvals(kinmem, njvevals)
    ccall((:KINSpilsGetNumJtimesEvals, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, njvevals)
end

function KINSpilsGetNumFuncEvals(kinmem, nfevals)
    ccall((:KINSpilsGetNumFuncEvals, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, nfevals)
end

function KINSpilsGetLastFlag(kinmem, flag)
    ccall((:KINSpilsGetLastFlag, libsundials_kinsol), Cint, (KINMemPtr, Ptr{Clong}), kinmem, flag)
end

function KINSpilsGetReturnFlagName(flag::Clong)
    ccall((:KINSpilsGetReturnFlagName, libsundials_kinsol), Cstring, (Clong,), flag)
end

function KINSpilsGetReturnFlagName(flag)
    KINSpilsGetReturnFlagName(convert(Clong, flag))
end


text:
using .Plots
using .Plots: Plot

export makeplot, set_palette!, set_main_color!, default_palette, run_anim!

# color definitions
function default_palette()
    return palette([
        RGB(0.5586, 0.1211, 0.4688),
        RGB(0.0000, 0.4470, 0.7410),
        RGB(0.8500, 0.3250, 0.0980),
        RGB(0.9290, 0.6940, 0.1250),
        RGB(0.4660, 0.6740, 0.1880),
        RGB(0.3010, 0.7450, 0.9330),
        RGB(0.6350, 0.0780, 0.1840),
    ])
end

const PALETTE = Ref{ColorPalette}(default_palette())
const MAIN_COLOR = Ref{RGB}(default_palette()[1])

function set_main_color!(c)
    MAIN_COLOR[] = parse(Colorant, c)
end

function set_palette!(pal; set_main = true)
    pal = palette(pal)
    PALETTE[] = pal
    set_main && set_main_color!(pal[1])
    return
end

function get_color(id)
    pal = PALETTE[]
    return pal[mod1(id, length(pal))]
end
get_main_color() = MAIN_COLOR[]

# object shapes
function rectangle(pos::Point, w, h)
    x, y = pos
    return Shape(x .+ [0, w, w, 0], y .+ [0, 0, h, h])
end

function circle(pos::Point, r = 1; step = π/20)
    x, y = pos
    θ = range(-π, π; step)
    return Shape(x .+ r*sin.(θ), y .+ r*cos.(θ))
end

function circle_section(pos::Point, r = 1, θ0  = 0, θmax = 2π; step = π/20)
    x, y = pos
    coors = euclidian.(r, θ0 .+ range(-θmax/2, θmax/2; step), Ref(pos))

    return Shape(vcat(x, first.(coors)), vcat(y, last.(coors)))
end

# recipes
@recipe function f(s::RoomRectangle; px_meter = 100)
    w, h = s.width, s.height

    # set plot style
    seriestype  :=  :shape
    fillalpha := 0
    xlims := (- 0.01*w, 1.01*w)
    ylims := (- 0.01*h, 1.01*h)
    aspect_ratio := :equal
    label := ""
    linecolor --> get_main_color()
    linewidth --> 2
    legend --> false
    axis --> nothing
    framestyle --> :none
    size --> px_meter .* (w, h)

    # remove kwargs
    delete!(plotattributes, :px_meter)

    return rectangle((0, 0), w, h)
end

roomlims(s::RoomRectangle) = (s.width, s.height)

@recipe f(s::Rectangle) = [s]
@recipe function f(s::Vector{<:Rectangle})
    pos = getproperty.(s, :pos)
    w =getproperty.(s, :width)
    h = getproperty.(s, :height)

    # set plot style
    seriestype  := :shape
    label := ""
    linecolor --> get_main_color()
    fillcolor --> get_main_color()

    return rectangle.(pos, w, h)
end

@recipe function f(d::Door; add_checkpoints = false, radius = 0.1)
    pos = [d.pos, d.pos .+ (d.width, 0)]
    x, y = first.(pos), last.(pos)

    if add_checkpoints
        pos = [d.pos .+ (radius, 0), d.pos .+ (d.width - radius, 0)]
        x = hcat(x, first.(pos))
        y = hcat(y, last.(pos))
    end

    # set plot style
    seriestype  := :path
    label := ""
    marker := add_checkpoints ? [:square :vline] : :square
    linecolor := add_checkpoints ? [:white :gray] : :white
    linewidth := add_checkpoints ? [4 1] : 4
    markercolor := add_checkpoints ? [get_main_color() :gray] : get_main_color()
    markerstrokecolor := add_checkpoints ? [get_main_color() :gray] : get_main_color()
    markersize := add_checkpoints ? [4 6] : 4

    # remove kwargs
    delete!(plotattributes, :add_checkpoints)
    delete!(plotattributes, :radius)

    return x, y
end

@recipe function f(c::Checkpoint)

    # set plot style
    seriestype  := :scatter
    label --> ""
    marker --> :x
    markersize --> 6
    markercolor --> :gray

    return [c.pos[1]], [c.pos[2]]
end

@recipe function f(c::HCheckpoint)
    pos = [c.pos, c.pos .+ (c.width, 0)]

    # set plot style
    seriestype  := :path
    label := ""
    marker := :vline
    linecolor := :gray
    linewidth := 1
    markercolor := :gray
    markerstrokecolor := :gray
    markersize := 6

    return first.(pos), last.(pos)
end

@recipe f(p::Pedestrian) = [p]
@recipe f(d::Dict{Int, Pedestrian}) = collect(values(d))
@recipe function f(ps::Vector{<:Pedestrian}; add_view = false, add_personal = true)

    props = extract_shapes.(ps; add_view, add_personal)
    shps = reduce(vcat, getindex.(props, 1))
    cls = reduce(vcat, getindex.(props, 2))
    ops = reduce(vcat, getindex.(props, 3))
    lnstyle = reduce(vcat, getindex.(props, 4))
    lnalpha = reduce(vcat, getindex.(props, 5))

    # set plot style
    seriestype  := :shape
    label --> ""
    fillcolor --> permutedims(cls)
    linecolor --> permutedims(cls)
    linestyle --> permutedims(lnstyle)
    linealpha --> permutedims(lnalpha)
    fillalpha --> permutedims(ops)

    # remove kwargs
    delete!(plotattributes, :add_view)
    delete!(plotattributes, :add_personal)

    return shps
end

function extract_shapes(p::Pedestrian; add_view = false, add_personal = true)
    inds = [1]
    shps = [circle(p.pos, p.rlims[1])]
    if add_personal
        push!(shps, circle(p.pos, p.radius))
        push!(inds, 2)
    end
    if add_view
        push!(shps, circle_section(p.pos, norm(p.vel, 2), direction_angle(p.vel), p.φ))
        push!(inds, 3)
    end

    cls = fill(get_color(p.id), 3)[inds]
    ops = [1, 0, 0.2][inds]
    lnstyle = [:solid, :dash, :solid][inds]
    lnalpha = [1, 1, 0][inds]
    return shps, cls, ops, lnstyle, lnalpha
end

# room plot
function makeplot(
    room::Room;
    add_checkpoints = false,
    add_valid = false,
    k_valid = 300,
    radius = 0.1,
    kwargs...
)
    # basic setup
    plt = plot(; kwargs...)

    # plot walls
    plot!(plt, room.shape)

    # add doors
    for d in vcat(room.entrances, room.exits)
        plot!(plt, d; add_checkpoints, radius)
    end

    # add checkpoints
    if add_checkpoints && !isempty(room.checkpoints)
        for c in room.checkpoints
            plot!(plt, c)
        end
    end

    # add obstacles
    for o in room.obstacles
        plot!(plt, o)
    end

    # add valid positions
    if add_valid
        lims = roomlims(room.shape)
        xs = range(0, lims[1]; length = k_valid)
        ys = range(0, lims[2]; length = k_valid)
        heatmap!(plt, xs, ys, (x, y) -> isvalid(room, (x, y), radius);
            color = [get_main_color(), :white],
            opacity = 0.2,
            cbar = false,
        )
    end
    return plt
end

# model plot
function makeplot(
    model::Model;
    add_view = false,
    add_personal = true,
    kwargs...
)

    plt = makeplot(model.room; kwargs...)
    if !isempty(model.pedestrians)
        plot!(plt, model.pedestrians; add_view, add_personal)
    end
    return plt
end

function run_anim!(m::Model, iter, filename::String; at::Int = 5, fps::Int = 20)

    bar = Progress(iter; showspeed = true)
    anim = Animation()
    for i in 1:iter
        step!(m)

        # progress bar
        k = length(keys(m.pedestrians))
        next!(bar; showvalues = [
                (:iter, i),
                (:pedestrians, k)
        ])
        if mod(iter, at) == 0
            title = @sprintf "Pedestrians: %4.1f s" m.Δt * m.iter[]
            makeplot(m; title)
            frame(anim)
        end
        if m.stop[]
            finish!(bar)
            break
        end
    end
    return gif(anim, filename; fps)
end


text: using Test, JuLIP


text:
using FountainCodes, Test, Distributions

"""return an LT code object and a vector of source symbols"""
function init(K; Tv=GF256, M=K, δ=1e-9)
    dd = Soliton(K, M, δ)
    LT{Tv}(K, dd)
end

function test_decoding(K::Integer, N::Integer; Tv=GF256)
    rng = MersenneTwister(123)
    code = init(K; Tv)
    src = rand_nonzero(rng, Tv, K)
    G = generator_matrix(code, 1:N)
    enc = G'*src
    dec = decode(G, enc)
    @test dec == src
end
test_decoding(100, 120)

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
module ObjectDetectors

using GeometryBasics

export AbstractObjectDetection, DetectorType, HUMAN, MACHINE,
       class_name, confidence, image_width, image_height, detector_type, detector_name,
       centroid, bounding_box, bounding_box_detection, iou
export AbstractObjectDetector, detect
export BoundingBoxDetection, create_with_center
export OrientedBoundingBoxDetection, width, height, orientation, rotate_point
export PolygonDetection, vertices
export StaticObjectDetector, create_polygon

include("abstract_object_detection.jl")
include("bounding_box_detection.jl")
include("oriented_bounding_box_detection.jl")
include("polygon_detection.jl")

include("abstract_object_detector.jl")
include("static_object_detector.jl")

end # module


text:
# Default to nothing happening
reset_mask!(x) = nothing

"""
    DEQChain(pre_deq, deq, post_deq)
    DEQChain(layers...)

A Sequential Model containing a DEQ.

!!! note
    The Model should contain exactly 1 `AbstractDEQ` Layer
"""
struct DEQChain{P1,D<:AbstractDeepEquilibriumNetwork,P2}
    pre_deq::P1
    deq::D
    post_deq::P2
end

function DEQChain(layers...)
    pre_deq, post_deq, deq, encounter_deq = [], [], nothing, false
    for l in layers
        if typeof(l) <: AbstractDeepEquilibriumNetwork
            @assert !encounter_deq "Can have only 1 DEQ Layer in the Chain!!!"
            deq = l
            encounter_deq = true
            continue
        end
        push!(encounter_deq ? post_deq : pre_deq, l)
    end
    @assert encounter_deq "No DEQ Layer in the Chain!!! Maybe you wanted to use Chain"
    pre_deq = length(pre_deq) == 0 ? NoOpLayer() : (length(pre_deq) == 1 ? pre_deq[1] : FChain(pre_deq...))
    post_deq = length(post_deq) == 0 ? NoOpLayer() : (length(post_deq) == 1 ? post_deq[1] : FChain(post_deq...))
    return DEQChain(pre_deq, deq, post_deq)
end

Flux.@functor DEQChain

function (deq::DEQChain)(x; kwargs...)
    x1 = deq.pre_deq(x)
    x2, deq_soln = deq.deq(x1; kwargs...)
    x3 = deq.post_deq(x2)
    return x3, deq_soln
end

function get_and_clear_nfe!(model::DEQChain)
    nfe = model.deq.stats.nfe
    model.deq.stats.nfe = 0
    return nfe
end

function Base.show(io::IO, model::DEQChain)
    l1 = length(destructure_parameters(model)[1])
    println(io, "DEQChain(")
    print(io, "\t")
    show(io, model.pre_deq)
    print(io, "\n\t")
    show(io, model.deq)
    print(io, "\n\t")
    show(io, model.post_deq)
    return print(io, "\n) $l1 Trainable Parameters")
end


text:
mutable struct Variable
    data::Array
    grad::Array
    Variable(data) = new(data)
end

abstract type _Function end
function _function(f::_Function, forward::Function, input::Variable)::Variable
    x = input.data
    y::Array = forward(x)
    output = Variable(y)
    f.input = input
    return output
end

mutable struct Square <: _Function
    input::Variable
    Square() = new()
end
(s::Square)(input::Variable)::Variable = _function(s, x -> x .^ 2, input)
backward(s::Square, gy::Array) = 2 * s.input.data .* gy

mutable struct Exp <: _Function
    input::Variable
    Exp() = new()
end
(s::Exp)(input::Variable)::Variable = _function(s, x -> exp.(x), input)
backward(e::Exp, gy::Array) = exp.(e.input.data) .* gy


text:
using PyPlot
using UlamSpiral

# text spirals
text_spiral()
savefig("plots/ulam_spiral_0.pdf")
savefig("plots/ulam_spiral_0.png")

text_spiral( 100; fig_n=2, fontsize=16)
savefig("plots/ulam_spiral_1.pdf")
savefig("plots/ulam_spiral_1.png")

text_spiral(fig_n=3, color_comp="white", color_prime="black", width=20, fontsize=10)
savefig("plots/ulam_spiral_2.pdf")
savefig("plots/ulam_spiral_2.png")

# points
point_spiral()
savefig("plots/ulam_spiral_3.pdf")
savefig("plots/ulam_spiral_3.png")

point_spiral( 1000; fig_n=4, markersize=10)
savefig("plots/ulam_spiral_4.pdf")
savefig("plots/ulam_spiral_4.png")

point_spiral(fig_n=5, markersize=1, color="black", width=200)
plot_quadratic(; fig_n=5, width=200);
plot_quadratic(4, 2, 41; fig_n=5, width=200, color_prime="green");

savefig("plots/ulam_spiral_5.pdf")
savefig("plots/ulam_spiral_5.png")

text_spiral(fig_n=6, color_comp="white", color_prime="black", width=20, fontsize=10)
point_spiral(fig_n=6, markersize=3, marker="s", color="green", width=20, clear="off")
savefig("plots/ulam_spiral_6.pdf")
savefig("plots/ulam_spiral_6.png")


text:
"""
    Region(SPParams, TMParams; recurrent= true)
A Region is a set of neurons arranged in minicolumns with 2 input gates:
- _proximal_ (size: number of minicolumns Nc): feedforward input that activates some of the region's minicolumns
- _distal_ (size: number of neurons Nn): contextual input that causes individual neurons to become predictive and fire preferentially on the next time step

Each gate accepts 1 or more input SDRs of the correct size, which will be ORed together.
If `recurrent= true`, the distal gate implicitly receives the same layer's state from the last time step as input (in addition to any SDR given explicitly).

The activation of minicolumns resulting from the proximal synapses is defined by the [`SpatialPooler`](@ref) algorithm,
while the activation of individual neurons is defined by the [`TemporalMemory`](@ref).

# Example

```julia
r= Region(SPParams(), TMParams(), recurrent= true)
feedforward_input= bitrand(Nc(r))

# activate based on feedforward input and the region's previous activation
active, predictive= r(feedforward_input)
# learn (adapt synapses)
step!(r, feedforward_input)

# use explicit contextual input eg from another region, that will be ORed to the region's state
distal_input= sprand(Bool,Nn(r), 0.003)|> BitVector
active, predictive= r(feedforward_input, distal_input)
```
"""
struct Region
  sp::SpatialPooler
  tm::TemporalMemory
end
Region(spp::SPParams,tmp::TMParams; recurrent=true, distal_input_size=0)= begin
  @assert prod(spp.szₛₚ) == tmp.Nc "Number of minicolumns in Spatial Pooler and Temporal Memory in the same region must match"
  Region(SpatialPooler(spp), TemporalMemory(tmp, recurrent=recurrent, distal_input_size=distal_input_size))
end
Nc(r::Region)= r.tm.params.Nc
Nₙ(r::Region)= Nₙ(r.tm)
Base.size(r::Region)= (r.sp.params.szₛₚ..., r.tm.params.k)
Base.length(r::Region)= Nₙ(r)

"`distalSynapses(r::Region)` is the adjacency matrix of the Region's connected distal synapses in {pre- × post-}synaptic neuron format."
distalSynapses(r::Region)= Wd(r.tm)
"`proximalSynapses(r::Region)` is the adjacency matrix of the Region's connected proximal synapses in {pre- × post-}synaptic neuron format."
proximalSynapses(r::Region)= Wₚ(r.sp)

"""
    (r::Region)(proximal, distal=falses(0))

Stimulates the region with the given inputs, without adapting it (pure function).
If distal input isn't given explicitly, recurrent connections are still used (if present).

See also [`Region`](@ref)
"""
(r::Region)(proximal, distal=falses(0))= @chain gateCombine(proximal) begin
  r.sp
  r.tm(_, gateCombine(distal))
end
activate(r::Region,proximal,distal)= r(proximal,distal)

# Support easier chaining of regions
RegionOutput= NamedTuple{(:active, :predictive, :bursting)}
(r::Region)(proximal::RegionOutput, distal=falses(0))= r(proximal.active, distal)

"""
    step!(r::Region, proximal, distal=falses(0); learn= true)

Stimulates the region with the given inputs, adapting it following the learning rules.
If distal input isn't given explicitly, recurrent connections are still used (if present).
If `learn == false`, the synapses are not adapted.

See also [`Region`](@ref)
"""
step!(r::Region, proximal, distal=falses(0); learn= true)= @chain gateCombine(proximal) begin
  step!(r.sp, _, learn)
  step!(r.tm, _, gateCombine(distal), learn)
end

"`reset!(r::Region)` unlearns the region's synapses (and other internal state), restoring it to the state at construction"
reset!(r::Region)= begin
  reset!(r.sp)
  reset!(r.tm)
end

"""
    connectContext!(output::Region, input::Region)

Connect contextually 2 regions `output -> input` by wiring the activation of
`output` to the context gate of `input`. Only `input` changes.
"""
#connectContext!(input::Region, output::Region)

gateCombine(x::Vector{T} where T <: CellActivity)= length(x) > 0 ? reduce(.|, x) : falses(0)
gateCombine(x)= x


text:
#=
2D explicit BTE. This assumes MPI is being used
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("FVbte2d");

useLog("FVbte2dlog", level=3)

# constants and various functions are in another file
include("bte-parameters.jl")
# A set of callback functions for the boundary condition
include("bte-boundary.jl")

# Configuration setup
domain(2)
solverType(FV)
timeStepper(EULER_EXPLICIT)
dt = 2.5e-12;
nsteps = 10;
setSteps(dt, nsteps);

# direction and band numbers
ndirs = 16;
total_bands = 40;
np = Finch.config.num_procs;
rank = Finch.config.proc_rank;
cell_partitions = 1; # This will only be increased if at least 2*total_bands procs
# Divide bands among processes for primarily band based parallel
if np >= total_bands # If there are extra processes, use them for cell partitioning
    cell_partitions = Int(floor(np / total_bands));

    nbands = 1;
    low_band = Int(floor(rank/cell_partitions)) + 1;
    high_band = low_band;

else
    nbands = Int(floor(total_bands / np));
    low_band = rank*nbands + 1;
    high_band = (rank+1)*nbands;
    if mod(total_bands, np) > rank
        nbands += 1;
        low_band += rank;
        high_band += rank+1;
    end
end
MPI = Finch.MPI;
if cell_partitions > 1
    partition_num = Int(floor((rank+0.5) * cell_partitions / Finch.config.num_procs));
    cell_comm = MPI.Comm_split(MPI.COMM_WORLD, partition_num, rank); # communicator for this mesh partition
    # cell_rank = MPI.Comm_rank(cell_comm);
else
    cell_comm = MPI.Comm_dup(MPI.COMM_WORLD);
end
num_band_partitions = MPI.Comm_size(cell_comm);
if num_band_partitions > total_bands
    println("Error: number of band partitions = "*string(num_band_partitions)*", but there are only "*string(total_bands)*" bands.");
    exit(1);
end

# A simple mesh is internally generated for convenience
# This matches the mesh in model_setup_BTE.in
mesh(QUADMESH, # quad elements
    elsperdim=[20,5], # elements in each direction: 20 x 5 uniform grid
    interval=[0, 4e-6, 0, 1e-6],  # interval in each direction: a very small rectangle
    bids=3, # 3 boundary IDs for this mesh correspond to left, right, top/bottom
    partitions=cell_partitions) # If there are enough procs to also do cell partitioning

println(string(rank)*":"*string(Finch.config.partition_index)*", "*string(Finch.fv_grid.neighboring_partitions));

# Indices, Variables, etc.
direction = index("direction", range = [1,ndirs])
band = index("band", range = [1,nbands])

# These are all set as variables because they are unknown, but only I is solved for in the PDE.
I = variable("I", type=VAR_ARRAY, location=CELL, index = [direction, band]) # Intensity
Io = variable("Io", type=VAR_ARRAY, location=CELL, index = [band]) # Equilibrium intensity for each band
beta = variable("beta", type=VAR_ARRAY, location=CELL, index = [band]) # Relaxation time scale
temperature = variable("temperature", location=CELL) # temperature of each cell
G_last = variable("G_last", type=VAR_ARRAY, location=CELL, index = [band]) # integrated intensity from last step
G_next = variable("G_next", type=VAR_ARRAY, location=CELL, index = [band]) # integrated intensity for current step

# Coefficients and related numbers
(dir_x, dir_y, reflect) = get_directions_2d(ndirs)
(center_freq, delta_freq) = get_band_frequencies(total_bands, low_band, high_band);
group_v = get_group_speeds(center_freq);

# These are set as coefficients because they have known values.
Sx = coefficient("Sx", dir_x, type=VAR_ARRAY) # direction x component
Sy = coefficient("Sy", dir_y, type=VAR_ARRAY) # direction y component
vg = coefficient("vg", group_v, type=VAR_ARRAY) # group speed

# This extra storage is allocated due for parallel temperature update
num_cells = length(temperature.values);
uold = zeros(num_cells*num_band_partitions); # These strange names are taken from the fortran code
gnb = zeros(num_cells*num_band_partitions);
unew = zeros(num_cells*num_band_partitions);
gna = zeros(num_cells*num_band_partitions);
uprime = zeros(num_cells*num_band_partitions);

boundary(I, 1, FLUX, "isothermal_bdry(I, vg, Sx, Sy, band, direction, normal, 305)") # left (ID=1)
boundary(I, 2, FLUX, "isothermal_bdry(I, vg, Sx, Sy, band, direction, normal, 300)") # right (ID=2)
boundary(I, 3, FLUX, "symmetric_bdry(I, vg, Sx, Sy, band, direction, normal)") # top and bottom (ID=3)

init_temp = 300; # The initial equilibrium temperature everywhere
initial(I, [equilibrium_intensity(center_freq[b], delta_freq, init_temp) for d=1:ndirs, b=1:nbands])
initial(Io, [equilibrium_intensity(center_freq[b], delta_freq, init_temp) for b=1:nbands])
initial(beta, [get_time_scale(center_freq[b], init_temp) for b=1:nbands])
initial(temperature, init_temp);

# The flux and source terms of the conservation equation
# F and S in the following equation:
# Dt(int(u dx)/A) = int(S dx) - int(F.n ds)
# BTE:
# Dt(int(Iij dx)) = int((Io-Iij)*beta dx) ) - vg * int(Iij * Si.n ds)
flux(I, "vg[band] * upwind([Sx[direction];Sy[direction]] , I[direction,band])")
source(I, "(Io[band] - I[direction,band]) ./ beta[band]") # <- Shouldn't this be multiplied by beta, not divided?

assemblyLoops(I, [band, "elements", direction])

# Create an array to hold the values of I from the last step.
# To get initial values here we have to manually initialize.
evalInitialConditions();
get_integrated_intensity!(G_last.values, I.values, ndirs, nbands);

# After each time step the temperature, equilibrium I, and time scales are updated
function post_step()
    update_temperature(temperature.values, I.values, center_freq, delta_freq, true);
end
postStepFunction(post_step);

exportCode("bte2dcode") # uncomment to export generated code to a file
#importCode("bte2dcodein") # uncomment to import code from a file

if rank==0
    start = Base.Libc.time()
end

solve(I)

if rank==0
    time = Base.Libc.time() - start;
    println(time)
end

finalize_finch()

##### Uncomment below to plot ######
if rank == 0

xy = Finch.fv_info.cellCenters;

for i=1:1
    println(string(xy[1,i]) * ", \t" * string(xy[2,i]) * ", \t" * string(temperature.values[i]))
end

# using Plots
# pyplot();

# # p1 = plot(xy[1,:], xy[2,:], I.values[1,:], st=:surface)#, zlims=(0,Inf))
# # p2 = plot(xy[1,:], xy[2,:], I.values[2,:], st=:surface)#, zlims=(0,Inf))
# # p3 = plot(xy[1,:], xy[2,:], I.values[3,:], st=:surface)#, zlims=(0,Inf))
# # p4 = plot(xy[1,:], xy[2,:], I.values[4,:], st=:surface)#, zlims=(0,Inf))
# # p5 = plot(xy[1,:], xy[2,:], I.values[5,:], st=:surface)#, zlims=(0,Inf))
# # p6 = plot(xy[1,:], xy[2,:], I.values[6,:], st=:surface)#, zlims=(0,Inf))
# # p7 = plot(xy[1,:], xy[2,:], I.values[7,:], st=:surface)#, zlims=(0,Inf))
# # p8 = plot(xy[1,:], xy[2,:], I.values[8,:], st=:surface)#, zlims=(0,Inf))
# # p9 = plot(xy[1,:], xy[2,:], I.values[9,:], st=:surface)#, zlims=(0,Inf))
# # p10 = plot(xy[1,:], xy[2,:], I.values[10,:], st=:surface)#, zlims=(0,Inf))
# # p11 = plot(xy[1,:], xy[2,:], I.values[11,:], st=:surface)#, zlims=(0,Inf))
# # p12 = plot(xy[1,:], xy[2,:], I.values[12,:], st=:surface)#, zlims=(0,Inf))
# # p13 = plot(xy[1,:], xy[2,:], I.values[13,:], st=:surface)#, zlims=(0,Inf))
# # p14 = plot(xy[1,:], xy[2,:], I.values[14,:], st=:surface)#, zlims=(0,Inf))
# # p15 = plot(xy[1,:], xy[2,:], I.values[15,:], st=:surface)#, zlims=(0,Inf))
# # p16 = plot(xy[1,:], xy[2,:], I.values[16,:], st=:surface)#, zlims=(0,Inf))
# # display(plot(p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13, p14, p15, p16, layout=16))

# display(plot(xy[1,:], xy[2,:], temperature.values[:], st=:surface));

end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
export iauC2s
"""
P-vector to spherical coordinates.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   p      double[3]    p-vector

Returned:
   theta  double       longitude angle (radians)
   phi    double       latitude angle (radians)

Notes:

   1. The vector p can have any magnitude; only its direction is used.

   2. If p is null, zero theta and phi are returned.

   3. At either pole, zero theta is returned.

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauC2s(double p[3], double *theta, double *phi)

function iauC2s(p::AbstractVector{<:Real})

   # Allocate return value
   ref_theta = Ref{Float64}(0.0)
   ref_phi   = Ref{Float64}(0.0)

   ccall((:iauC2s, libsofa_c), Cvoid,
         (Ptr{Cdouble},
         Ref{Cdouble}, Ref{Cdouble}),
         convert(Array{Float64, 1}, p),
         ref_theta, ref_phi)

   return ref_theta[], ref_phi[]
end


text:
# ---
# title: 1439. Find the Kth Smallest Sum of a Matrix With Sorted Rows
# id: problem1439
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Heap
# link: <https://leetcode.com/problems/find-the-kth-smallest-sum-of-a-matrix-with-sorted-rows/description/>
# hidden: true
# ---
#
# You are given an `m * n` matrix, `mat`, and an integer `k`, which has its rows
# sorted in non-decreasing order.
#
# You are allowed to choose exactly 1 element from each row to form an array.
# Return the Kth **smallest** array sum among all possible arrays.
#
#
#
# **Example 1:**
#
#
#
#     Input: mat = [[1,3,11],[2,4,6]], k = 5
#     Output: 7
#     Explanation: Choosing one element from each row, the first k smallest sum are:
#     [1,2], [1,4], [3,2], [3,4], [1,6]. Where the 5th sum is 7.
#
# **Example 2:**
#
#
#
#     Input: mat = [[1,3,11],[2,4,6]], k = 9
#     Output: 17
#
#
# **Example 3:**
#
#
#
#     Input: mat = [[1,10,10],[1,4,5],[2,3,6]], k = 7
#     Output: 9
#     Explanation: Choosing one element from each row, the first k smallest sum are:
#     [1,1,2], [1,1,3], [1,4,2], [1,4,3], [1,1,6], [1,5,2], [1,5,3]. Where the 7th sum is 9.
#
#
# **Example 4:**
#
#
#
#     Input: mat = [[1,1,10],[2,2,9]], k = 7
#     Output: 12
#
#
#
#
# **Constraints:**
#
#   * `m == mat.length`
#   * `n == mat.length[i]`
#   * `1 <= m, n <= 40`
#   * `1 <= k <= min(200, n ^ m)`
#   * `1 <= mat[i][j] <= 5000`
#   * `mat[i]` is a non decreasing array.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Install dependencies

# install PkgFix if not present
if !isdir(joinpath(Pkg.dir(), "PkgFix"))
  Pkg.clone("https://github.com/OptimalDesignLab/PkgFix.jl.git")
end
start_dir=pwd()
cd(Pkg.dir("PkgFix"))
run(`git checkout upgrade_0.6`)
cd(start_dir)

using PkgFix  # from now on, use PkgFix instead of Pkg for everything

global const ARRAYVIEWS_URL = "https://github.com/JaredCrean2/ArrayViews.jl.git"
global const ARRAYVIEWS_VER = "work"

pkg_dict = PkgFix.installed()

if !haskey(pkg_dict, "ArrayViews")
  PkgFix.add(ARRAYVIEWS_URL, branch_ish=ARRAYVIEWS_VER)
else
  PkgFix.checkout("ArrayViews", ARRAYVIEWS_VER)
end


text:
module ASTMatching

struct UnionType{Base, Key, Arg} end

# Write your package code here.
include("patterns.jl")
include("treegen.jl")
include("codegen.jl")
include("macros.jl")

end


text:
#########
#File: c:\Users\digan\Dropbox\Dynamic_Networks\repos\ScoreDrivenExponentialRandomGraphs\_research\analysis_for_paper_revision\add_confidence_intervals\load_estimates_and_check_parametric_bootstrap_conf_intervals.jl
#Project: c:\Users\digan\Dropbox\Dynamic_Networks\repos\ScoreDrivenExponentialRandomGraphs\_research\analysis_for_paper_revision\add_confidence_intervals
#Created Date: Saturday April 24th 2021
#Author: Domenico Di Gangi,  <digangidomenico@gmail.com>
#-----
#Last Modified: Tuesday April 27th 2021 10:34:44 pm
#Modified By:  Domenico Di Gangi
#-----
#Description:
#-----
########

#region import and models
#%%
using DrWatson
@quickactivate "ScoreDrivenExponentialRandomGraphs"

using ProjUtilities
using Logging
using DataFrames
using ScoreDrivenERGM
import ScoreDrivenERGM:DynNets, Utilities, StaticNets, Scalings
using PyPlot
pygui(true)
#endregion

@elapsed df = collect_results( datadir("sims", "dgp&FIl_est"))
df["modelTag"] = string.(df["model"])

T = 300
N = 50
nSample = 50
indQuant = 1
modelTag = string(DynNets.SdErgmDirBin0Rec0_mle(scoreScalingType = "FISH_D"))

dgpSetting = DynNets.list_example_dgp_settings(DynNets.SdErgmDirBin0Rec0_mle()).dgpSetARlowlow

# # visual inspection of filter, dgp and conf bands
res = filter([:modelTag, :T, :N, :dgpSettings, :nSample, ] => (mTag,t,n, d, s) -> all((mTag==modelTag, t==T, n==N, d == dgpSetting, s == nSample)), df)[1,:]

using Statistics

model = res.model
quantilesVals = [0.975, 0.025]

n = 15
parDgp_T = res.allParDgpT[:,:,n]
parSdFilt_T = res.allfVecT_filt[:,:, n]
estSdResPar = res.allvEstSdResPar[:,n]
obsT = res.allObsT[n]
parSSFilt_T =  DynNets.estimate_single_snap_sequence(model, obsT)

# coverFlags = @sync @distributed  hcat
nQuant = length(quantilesVals)
confBandsT = zeros(model.staticModel.nErgmPar, T, 1, nQuant)
confBandsT_SS = zeros(model.staticModel.nErgmPar, T, 1, nQuant)
for t = 1:T
    parDgp_t = parDgp_T[:,t]
    parFiltSS_t = parSSFilt_T[:,t]
    parFiltSD_t = res.allfVecT_filt[:,t,n]
    confBandsT[:,t ,1, :] = mapslices((x -> quantile(x, quantilesVals)), StaticNets.get_par_boot_ergm_distrib(model.staticModel, parFiltSD_t, N; nSample = 100), dims=2)
    confBandsT_SS[:,t ,1, :] = mapslices((x -> quantile(x, quantilesVals)), StaticNets.get_par_boot_ergm_distrib(model.staticModel, parFiltSS_t, N; nSample = 100), dims=2)

end

DynNets.plot_filtered_and_conf_bands(model, N, parSdFilt_T, confBandsT; parDgpTIn=parDgp_T)
DynNets.plot_filtered_and_conf_bands(model, N, parSSFilt_T, confBandsT_SS; parDgpTIn=parDgp_T)

begin
obsT = res.allObsT[n]
ftot_0 = res.allftot_0[:,n]
modelIntegrated = deepcopy(model)
modelIntegrated.options["integrated"] = true

~, vEstSdResParIntegrated, fVecT_filt, ~, ~, conv_flag, ftot_0 = DynNets.estimate_and_filter(modelIntegrated, N, obsT; indTvPar = model.indTvPar, show_trace=true)

# vEstSdResParIntegrated[6] = 0.02
testIntFilt_T, _, sVecT, invScalT = DynNets.score_driven_filter( model, N, obsT, vEstSdResParIntegrated, model.indTvPar; ftot_0 = ftot_0 )
DynNets.plot_filtered(model, N, testIntFilt_T; parDgpTIn=parDgp_T)

end

invScalT
sVecT
#%%
## Sequence for rescaled Par Bootstrap confidence intervals for SDERGM
nParBootSample = 100
ftot_0 = res.allftot_0[:,n]

parBootSampMatsSdFilt_T_S = StaticNets.sample_ergm_sequence(model.staticModel, N, parSdFilt_T, nParBootSample)

#Estimate SD on data and filter
parBootSdFiltDistr = zeros(model.staticModel.nErgmPar, T, nParBootSample)

for pbs=1:nParBootSample
    @info "$pbs"
    obsT = DynNets.seq_of_obs_from_seq_of_mats(model, parBootSampMatsSdFilt_T_S[:,:,:,pbs] )

    # ~, vEstSdResPar, fVecT_filt, ~, ~, conv_flag, ftot_0 = DynNets.estimate_and_filter(model, N, obsT; indTvPar = model.indTvPar)
    # parBootSdFiltDistr[:,:, pbs] = fVecT_filt

    parBootSdFiltDistr[:,:, pbs] = DynNets.score_driven_filter( model, N, obsT, estSdResPar, model.indTvPar; ftot_0 = ftot_0 )[1]

end

nQuant = length(quantilesVals)
confBandsSdParBoot = mapslices((x -> quantile(x, quantilesVals)), parBootSdFiltDistr, dims=3)
confBandsSdParBoot = reshape(confBandsSdParBoot, (model.staticModel.nErgmPar, T, 1, nQuant))
DynNets.plot_filtered_and_conf_bands(model, N, parSdFilt_T, confBandsSdParBoot; parDgpTIn=parDgp_T)

# Obtain Distribution of observations from the filtered path, sampling each ergm_t using the filtered parameters nSample times

# filter each

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
function writecsvstep(t, x, ρ, u, v, p)
    file = open(joinpath("data", "data_$t.csv"), "w")
    write(file, "position,density,energy,velocity,pressure\n")
    writecsv(file, [x ρ u v p])
    close(file)
end


text:
using ExactCovers
using Base.Test

include("basic.jl")
include("queens.jl")
include("euler96.jl")


text:
using Flux
using Base.Iterators: repeated

##
# approximate this funcgtion on [0, 1]
# it is a small section of  the Lennard-Jones potential
f = x -> (1/(0.82+x))^12 -2 * (1/(.82+x))^6
# massage it so Flux can use it
trueFun(u) = [f(u[1])]

# Here we are taking measurements from the true model
# In real applications, we would just be given this
# data rather than generate it ourselves
nsamples = 200
xs = rand(nsamples)
samples = [ [x] for x in xs ]
y = trueFun.(samples)

xerr = range(0, 1, length=1000)
eval_nn(nn, x) = nn([x])[1].data
err_nn(nn) = maximum(abs, f.(xerr) - eval_nn.(Ref(nn), xerr))
maxfit_nn(nn) = maximum(abs, f.(xs) - eval_nn.(Ref(nn), xs))

function train_nn(nnodes)
    # This defines an ANN with an input layer, one hidden and an output layer
    nn = Chain(x -> x,             # input layer
               Dense(1, nnodes, Flux.sigmoid),   # first hidden layer
               Dense(nnodes, 1))        # output layer

    # extract the parameters:
    # the weights w^{[n]} and shifts b^{[n]}
    ps = Flux.params(nn)

    # errors at individual sample points
    sqerrors() = [ z[1]^2 for z in y - Flux.data(nn.(samples))]
    # loss functional => Least squares
    losssq() = sum(sqerrors())/length(samples)
    # for plotting the errors
    errcols() = [sqrt(e.data) for e in sqerrors()]

    cb() = nothing

    ## and now we can watch the training
    best_nn = deepcopy(nn)
    best_fit = 1e30

    for n = 1:nnodes, _ = 1:4
        Flux.train!(losssq, ps, repeated((), 2_000), ADAM(0.05), cb = cb)
        @show maxfit_nn(nn)
        if maxfit_nn(nn) < best_fit
            best_nn = deepcopy(nn)
        end
    end
    return nn
end

NN_nn = 3:5
errs_nn = Float64[]

for N in NN_nn
    @show N
    nn = train_nn(N)
    @show e = err_nn(nn)
    push!(errs_nn, e)
end
@show errs_nn

# ## Now produce plots of the two functions
# nn10 = train_nn(10)
# nn20 = train_nn(20)
# xp = range(0, 1, length = 300)
# y10 = eval_nn.(Ref(nn10), xp)
# y20 = eval_nn.(Ref(nn20), xp)
# plot(xp, f.(xp))
# plot!(xp, y10)
# plot!(xp, y20)

max_stars_repo_name: berkonat/ann-examples

max_issues_repo_name: berkonat/ann-examples

max_forks_repo_name: berkonat/ann-examples


text:
julia_version = v"1.4.2"
gcc_version = v"8"
include("../common.jl")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@testset "Visualization Diagnostics" begin
    air_passengers = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame
    log_air_passengers = log.(air_passengers.passengers)
    model = BasicStructural(log_air_passengers, 12)
    fit!(model)
    kf = kalman_filter(model)
    @test typeof(plotdiagnostics) <: Function
    @test length(methods(plotdiagnostics)) == 1
    @test typeof(plotdiagnostics!) <: Function
    @test length(methods(plotdiagnostics!)) == 2
    diagnostics = StateSpaceModels.PlotDiagnostics(kf)
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), diagnostics)
    @test length(rec) == 8
end


text:
using ADCME
using PyPlot

function residual_and_jacobian(θ, u)
    X = ae(u, config, θ) + 1.0     # (1)
    Xp = tf.gradients(X, u)[1]
    Xpp = tf.gradients(Xp, u)[1]
    up = [u[2:end];constant(zeros(1))]
    un = [constant(zeros(1)); u[1:end-1]]
    R = Xp .* ((up-un)/2h)^2 + X .* (up+un-2u)/h^2 - φ
    dRdu = Xpp .* ((up-un)/2h)^2 + Xp.*(up+un-2u)/h^2 - 2/h^2*X
    dRdun = -Xp[2:end]/h .* (up-un)[2:end]/2h + X[2:end]/h^2
    dRdup = Xp[1:end-1]/h .* (up-un)[1:end-1]/2h + X[1:end-1]/h^2
    J = spdiag(n-1,
        -1=>dRdun,
        0=>dRdu,
        1=>dRdup)     # (2)
    return R, J
end

config = [20,20,20,1]
n = 100
h = 1/n
x = collect(LinRange(0, 1.0, n+1))

φ = @. (1 - 2*x)*(-100*x^2*(2*x - 2) - 200*x*(1 - x)^2)/(100*x^2*(1 - x)^2 + 1)^2 - 2 - 2/(100*x^2*(1 - x)^2 + 1)
φ = φ[2:end-1]
θ = Variable(ae_init([1,config...]))
u0 = constant(zeros(n-1))
function L(u)    # (3)
  u_obs = (@. x * (1-x))[2:end-1]
  loss = mean((u - u_obs)^2)
end
loss, solution, grad = NonlinearConstrainedProblem(residual_and_jacobian, L, θ, u0)
X_pred = ae(collect(LinRange(0.0,0.25,100)), config, θ) + 1.0

sess = Session(); init(sess)
BFGS!(sess, loss, grad, θ)
x_pred, sol = run(sess, [X_pred, solution])

figure(figsize=(10,4))
subplot(121)
s = LinRange(0.0,0.25,100)
x_exact = @. 1/(1+100*s^2) + 1
plot(s, x_exact, "-", linewidth=3, label="Exact")
plot(s, x_pred, "o", markersize=2, label="Estimated")
legend()
xlabel("u")
ylabel("X(u)")

subplot(122)
s = LinRange(0.0,1.0,101)[2:end-1]
plot(s, (@. s * (1-s)), "-", linewidth=3, label="Exact")
plot(s, sol, "o", markersize=2, label="Estimated")
legend()
xlabel("x")
ylabel("u")
savefig("nn.png")


text:
"""
    function get_optimal(
        theta_init::Vector{Float64}, # initial point of parameters
        loss_func::Function; # lambda(theta)

        scale::Vector{Symbol} = fill(:direct, length(theta_init)),
        theta_bounds::Vector{Tuple{Float64,Float64}} = unscaling.(
            fill((-Inf, Inf), length(theta_init)),
            scale
            ),
        scan_tol::Float64 = 1e-3,
        loss_tol::Float64 = 1e-3,
        local_alg::Symbol = :LN_NELDERMEAD,
        silent::Bool = false,
        max_iter::Int = 10^5,
        loss_grad::Union{Function, Symbol} = :EMPTY
    )

Provides the optimization routine using the interface similar to [`get_endpoint`](@ref).
Currently it uses standard NLopt optimization but allows to use parameter scaling and autodiff method.

## Return
[`ProfilePoint`](@ref) object storing optimizations results.`
`value` and `loss` properties are equal to the optimal (minimal) `loss_func` value.

## Arguments
- `theta_init`: starting values of parameter vector ``\\theta``.
- `loss_func`: loss function ``\\Lambda\\left(\\theta\\right)`` for profile likelihood-based (PL) identification. Usually we use log-likelihood for PL analysis: ``\\Lambda( \\theta ) = - 2 ln\\left( L(\\theta) \\right)``.

## Keyword arguments
- `scale`: vector of scale transformations for each parameters' component. Possible values: `:direct` (`:lin`), `:log`, `:logit`. This option can speed up the optimization, especially for wide `theta_bounds`. The default value is `:direct` (no transformation) for all parameters.
- `theta_bounds`: vector of tuple `(lower_bound, upper_bound)` for each parameter. Bounds define the ranges for possible parameter values. Default bounds are `(-Inf,Inf)`.
- `scan_tol`: Absolute tolerance for all component of theta used as termination criterion.
- `loss_tol`: Absolute tolerance controlling `loss_func`.
- `local_alg`: algorithm of optimization. Derivative-free and gradient-based algorithms form NLopt package.
- `max_iter` : maximal number of optimizer iterations.
- `loss_grad` : This option declares the method for calcuting loss function gradient.
        This is required for gradient-based methods. The possible values
        - `:EMPTY` (default) not gradient is set. IT works only for gradient-free methods.
        - `:AUTODIFF` means autodifferentiation from `ForwardDiff` package is used.
        - `:FINITE` means finite difference method from `Calculus` is used.
        - It is also possible to set gradient function here `function(x::Vector{Float64})` which returns gradient vector.
"""
function get_optimal(
    theta_init::Vector{Float64}, # initial point of parameters
    loss_func::Function; # lambda(theta)

    scale::Vector{Symbol} = fill(:direct, length(theta_init)),
    theta_bounds::Vector{Tuple{Float64,Float64}} = unscaling.(
        fill((-Inf, Inf), length(theta_init)),
        scale
        ),
    scan_tol::Float64 = 1e-3,
    loss_tol::Float64 = 1e-3,
    local_alg::Symbol = :LN_NELDERMEAD,
    silent::Bool = false,
    max_iter::Int = 10^5,
    loss_grad::Union{Function, Symbol} = :EMPTY
)
    # dim of the theta vector
    n_theta = length(theta_init)

    # checking arguments
    # theta_bound[1] < theta_init < theta_bound[2]
    theta_init_outside_theta_bounds = .! [theta_bounds[i][1] < theta_init[i] < theta_bounds[i][2] for i in 1:n_theta]
    if any(theta_init_outside_theta_bounds)
        throw(ArgumentError("theta_init is outside theta_bound: $(findall(theta_init_outside_theta_bounds))"))
    end
    # 0 <= theta_bound[1] for :log
    less_than_zero_theta_bounds = (scale .== :log) .& [theta_bounds[i][1] < 0 for i in 1:n_theta]
    if any(less_than_zero_theta_bounds)
        throw(ArgumentError(":log scaled theta_bound min is negative: $(findall(less_than_zero_theta_bounds))"))
    end
    # 0 <= theta_bounds <= 1 for :logit
    less_than_zero_theta_bounds = (scale .== :logit) .& [theta_bounds[i][1] < 0 || theta_bounds[i][2] > 1 for i in 1:n_theta]
    if any(less_than_zero_theta_bounds)
        throw(ArgumentError(":logit scaled theta_bound min is outside range [0,1]: $(findall(less_than_zero_theta_bounds))"))
    end

    # checking arguments
    # when using :LN_NELDERMEAD initial parameters should not be zero
    if local_alg == :LN_NELDERMEAD
        zeroParameter = [isapprox(theta_init[i], 0., atol=1e-2) for i in 1:n_theta]
        if any(zeroParameter)
            @warn "Close-to-zero parameters found when using :LN_NELDERMEAD."
            show(findall(zeroParameter))
        end
    end

    # checking loss_grad
    is_gradient = occursin(r"^LD_", String(local_alg))
    if loss_grad == :EMPTY && is_gradient
        throw(ArgumentError("`loss_grad` must be set for gradient local fitter `$(local_alg)`"))
    end

    # progress info
    prog = ProgressUnknown("Fitter counter:"; spinner=false, enabled=!silent, showspeed=true)
    counter = 0
    supreme = nothing

    loss_func_rescaled = function(theta_g)
        theta = unscaling.(theta_g, scale)
        loss_value = loss_func(theta)

        counter += 1
        if (typeof(supreme) == Nothing || loss_value < supreme) && !isa(loss_value, ForwardDiff.Dual)
            supreme = loss_value
            ProgressMeter.update!(prog, counter, spinner="⠋⠙⠹⠸⠼⠴⠦⠧⠇⠏"; showvalues = [(:supreme,round(supreme; sigdigits=4))])
        end

        return loss_value
    end

    # transforming loss
    theta_init_g = scaling.(theta_init, scale)
    loss_func_g = if loss_grad == :EMPTY
        (theta_g, grad_g) -> loss_func_rescaled(theta_g) # no gradient, only for gradient-free methods
    elseif isa(loss_grad, Function)
        function(theta_g, grad_g) # gradient is directly declared
            # rescaling gradient function
            if length(grad_g) > 0
                theta = unscaling.(theta_g, scale)
                loss_grad_value = loss_grad(theta)
                for i in 1:n_theta
                    if scale[i] == :log
                        grad_g[i] = loss_grad_value[i] * theta[i] * log(10.)
                    elseif scale[i] == :logit
                        grad_g[i] = loss_grad_value[i] * theta[i] * (1. - theta[i]) * log(10.)
                    else # for :direct and :lin
                        grad_g[i] = loss_grad_value[i]
                    end
                end
            end

            return loss_func_rescaled(theta_g)
        end
    elseif loss_grad == :AUTODIFF
        function(theta_g, grad_g) # gradient is directly declared
            ForwardDiff.gradient!(grad_g, loss_func_rescaled, theta_g)
            return loss_func_rescaled(theta_g)
        end
    elseif loss_grad == :FINITE
        function(theta_g, grad_g) # gradient is directly declared
            Calculus.finite_difference!(loss_func_rescaled, theta_g, grad_g, :central)
            return loss_func_rescaled(theta_g)
        end
    end

    opt = Opt(local_alg, n_theta)
    opt.ftol_abs = loss_tol
    opt.xtol_abs = scan_tol
    opt.min_objective = loss_func_g
    opt.maxeval = max_iter

    # version 1: internal :LN_AUGLAG box constrains
    theta_bounds_g = scaling.(theta_bounds, scale)
    opt.lower_bounds = [tb[1] for tb in theta_bounds_g]
    opt.upper_bounds = [tb[2] for tb in theta_bounds_g]

    (optf, optx_g, ret) = optimize(opt, theta_init_g)
    optx = unscaling.(optx_g, scale)

    pp = ProfilePoint(optf, loss_func(optx), optx, ret, counter)

    return pp
end


text:
include("setup_1d.jl")

y = copy(x0)

println("dwt! by filtering (N=",N,"), ", L, " levels")
f(y, x0, L, wf) = for i = 1:tn; dwt!(y, x0, wf, L); end
f(y, x0, L, wf);
@time f(y, x0, L, wf);


text:
function _naive_history(block::Block{T}, window::Int) where {T}
    @assert window > 0
    times = DateTime[]
    values = Vector{T}[]
    buffer = CircularBuffer{T}(window)
    for (time, value) in block
        # Push a new value onto the buffer. (If the buffer is full, it will overwrite the
        # value that is falling out of the buffer, because it is circular.)
        push!(buffer, value)
        isfull(buffer) || continue
        push!(times, time)
        push!(values, copy(buffer))
    end
    return Block(times, values)
end

function _test_history(block::Block)
    n_in = block_node(block)
    for window in 1:length(block)
        expected = _naive_history(block, window)
        result = _eval(history(n_in, window))
        @test expected == result
    end
end

@testset "invalid" begin
    @test_throws ArgumentError history(n1, 0)
end

@testset "simple" begin
    _test_history(b3)
    _test_history(b4)
    _test_history(b_boolean)
end

@testset "vector" begin
    _test_history(_get_rand_svec_block(MersenneTwister(42), 3, 10))
    _test_history(_get_rand_vec_block(MersenneTwister(42), 3, 10))
end