text:
import Base: *, isless

# Auxiliary functions
islastsibling(node::Node) = node.sibling === node

struct LCS{T}
    v::Vector{Tuple{T, T}}
    l::Int
end

LCS(t::Tuple{T, T}, l::Int) where T = LCS{T}([t], l)
isless(seq1::LCS, seq2::LCS) = seq1.l < seq2.l
*(seq1::LCS{T}, seq2::LCS{T}) where T = LCS{T}(vcat(seq1.v, seq2.v), seq1.l + seq2.l)

# No one expects the dynamic programming :-)
# Generally LCS actually can return more than one sequence, but in practice greedy solution should be enough
function lcs(node1::T, node2::T,
    eq = (x, y) -> x == y, is_wildcard = x -> false,
    memo = Dict{Tuple{T, T}, LCS{T}}()) where T

    init = (node1, node2)
    init in keys(memo) && return memo[init]

    if (!is_wildcard(node1) && !is_wildcard(node2) && !eq(node1, node2)) ||
            (!is_wildcard(node1) && is_wildcard(node2)) ||
            (is_wildcard(node1) && !is_wildcard(node2))

        seq1 = islastsibling(node2) ? LCS{T}([], 0) : lcs(node1, node2.sibling, eq, is_wildcard, memo)
        seq2 = islastsibling(node1) ? LCS{T}([], 0) : lcs(node1.sibling, node2, eq, is_wildcard, memo)
        memo[init] = max(seq1, seq2)
   elseif !is_wildcard(node1) && !is_wildcard(node2) && eq(node1, node2)
       if islastsibling(node1) || islastsibling(node2)
           memo[init] = LCS((node1, node2), 1)
       else
           memo[init] = LCS((node1, node2), 1) * lcs(node1.sibling, node2.sibling, eq, is_wildcard, memo)
       end
   else
       seq1 = islastsibling(node2) ? LCS{T}([], 0) : lcs(node1, node2.sibling, eq, is_wildcard, memo)
       seq2 = islastsibling(node1) ? LCS{T}([], 0) : lcs(node1.sibling, node2, eq, is_wildcard, memo)

       seq3 =  if isleaf(node1) || isleaf(node2)
                   LCS((node1, node2), 1)
               else
                   subseq3 = lcs(node1.child, node2.child, eq, is_wildcard, memo)
                   LCS((node1, node2), subseq3.l + 1)
               end
       seq4 = islastsibling(node1) || islastsibling(node2) ? LCS{T}([], 0) : lcs(node1.sibling, node2.sibling, eq, is_wildcard, memo)

       memo[init] = max(seq3*seq4, seq1, seq2)
   end

   return memo[init]
end

# sometimes root node get random narrow branches, in order to compare
# profiles it's easier to remove this spurious branches
function simplify!(node)
    children = collect(node)
    sort!(children, by = x -> length(x.data.span), rev = true)
    for i in 2:length(children)
        prunebranch!(children[i])
    end

    node
end

span(node) = length(node.data.span)
delta(node1, node2) = span(node1) - span(node2)
sf(node) = string(node.data.sf)

function _add_diff_node_to_data(target::T, baseline::T, indent = 0, data = Vector{NamedTuple{(:level, :delta, :x1, :x2, :sf), Tuple{Int,Int,Float64,Float64,String}}}(undef, 0)) where T
    push!(data,
            (level=indent,
            delta=delta(target, baseline),
            x1=target.data.span.start,
            x2=target.data.span.stop+.9,
            sf="$(sf(target))"))
    if target == baseline
        for child in target
            _add_diff_node_to_data(child, child, indent + 1, data)
        end
    else
        seq = lcs(target.child, baseline.child, (x, y) -> sf(x) == sf(y), x -> occursin(r"\[\d+\]:\d+$", sf(x)))
        seq = Dict(seq.v)
        for child in target
            _add_diff_node_to_data(child, get(seq, child, child), indent + 1, data)
        end
    end

    data
end

function _plotdiffflamegraph(target, baseline; kwargs...)
    data = _add_diff_node_to_data(target, baseline)
    domain = sort([extrema(map(x -> x.delta, data))..., 0])
    if domain[1] == 0
        range = ["#d0d0d0", "#d73027"]
        domain = domain[[1, 3]]
    elseif domain[3] == 0
        range = ["#528fe0", "#d0d0d0"]
        domain = domain[[1, 3]]
    else
        range = ["#528fe0", "#d0d0d0", "#d73027"]
    end

    return data |> VegaLite.@vlplot(
        mark={type=:rect, stroke="#505050"},
        transform=[{calculate="datum.level+0.9", as=:level2}],
        selection={grid={type=:interval, bind=:scales},
            highlight={type=:single, on=:mouseover, empty=:none},
            select={type=:multi, empty=:none}},
        x={:x1, axis=nothing},
        x2=:x2,
        y={"level:q", axis=nothing},
        y2={"level2:q"},
        fillOpacity={condition=[
                {selection=:highlight, value=0.6},
                {selection=:select, value=0.5}
            ],
            value=1},
        strokeWidth={condition=[
            {selection=:highlight, value=0.5},
            {selection=:select, value=0.5}
            ],
            value=0},
        width=get(kwargs, :width, width(_options)),
        height=get(kwargs, :height, height(_options)),
        color={field="delta", type="quantitative",
            scale={
                domain=domain,
                range=range
            },
            legend={title=nothing, orient=:bottom}},
        tooltip=:sf,
        title="Diff Profile Results"
    )
end

"""
    ProfileVega.diffview(baseline, target=Profile.fetch(); kwargs...)

View differential flamegraph which allow one to compare performance before and after a code change. Colors corresponds to
`target` - `baseline` delta,
If a frame appeared more times in `target`, it is red, less times, it is blue. The saturation is relative to the delta.

Typically it should be used together with `@profbase`. `@profdiffview` accepts all [FlameGraphs](https://github.com/timholy/FlameGraphs.jl) `kwargs` as well as
the following options:

    - `width` - width of the chart in pixels
    - `height`- height of the chart in pixels
    - `negate` - default: false. If `false` shows `target` flamechart with colorization relative to `baseline`.
        If `true` shows `baseline` flamechart with colorization relative to `target`.
    - `simplify` - default: true. Removes spurious branches (all small span branches originated from root).

"""
function diffview(baseline, target::Vector{UInt64}=Profile.fetch(); kwargs...)
    target = FlameGraphs.flamegraph(target; filter(x -> !(x[1] in fieldnames(Options)), kwargs)...)
    baseline = FlameGraphs.flamegraph(baseline; filter(x -> !(x[1] in fieldnames(Options)), kwargs)...)

    if get(kwargs, :simplify, _options.simplify)
        target = simplify!(target)
        baseline = simplify!(baseline)
    end

    if get(kwargs, :negate, _options.negate)
        return _plotdiffflamegraph(baseline, target; kwargs...)
    else
        return _plotdiffflamegraph(target, baseline; kwargs...)
    end
end

"""
    @profbase f(args...)

Returns the stack frame which can be used later in `diffview` function. Typical usage of this macro

```
baseline = @profbase myfunction()
```
"""
macro profbase(ex)
    return quote
        Profile.clear()
        Profile.@profile $(esc(ex))
        deepcopy(Profile.fetch())
    end
end

"""
    @profdiffview baseline f(args...)

Clear the Profile buffer, profile `f(args...)`, and view differential flamegraph which allow one to
compare performance before and after a code change. Colors corresponds to `target` - `baseline` delta,
If a frame appeared more times in `target`, it is red, less times, it is blue. The saturation is relative to the delta.

Typically it should be used together with `@profbase`. `@profdiffview` accepts all [FlameGraphs](https://github.com/timholy/FlameGraphs.jl) `kwargs` as well as
the following options:

    - `width` - width of the chart in pixels
    - `height`- height of the chart in pixels
    - `negate` - default: false. If `false` shows `target` flamechart with colorization relative to `baseline`.
        If `true` shows `baseline` flamechart with colorization relative to `target`.
    - `simplify` - default: true. Removes spurious branches (all small span branches originated from root).

#### Example
```
baseline = @profbase myfun()
@profdiffview baseline mymodifiedfun() combine=false width=1200 height=800
```
"""
macro profdiffview(baseline, ex, args...)
    return quote
        Profile.clear()
        Profile.@profile $(esc(ex))
        diffview($(esc(baseline)); $(esc.(args)...))
    end
end

max_stars_repo_name: davidanthoff/ProfileVega

max_issues_repo_name: davidanthoff/ProfileVega

max_forks_repo_name: davidanthoff/ProfileVega


text:
# ==============================================================================
# Parseur lisant une instance de probleme de partitionnement d'ensembles (SPA) bi-objectif

function loadInstance2SPA(fname::String)
    f = open(fname)
    nbcontraintes, nbvar = parse.(Int, split(readline(f))) # nombre de contraintes , nombre de variables
    L = zeros(Int, nbcontraintes, nbvar)    # matrice des contraintes
    c1 = zeros(Int, nbvar)                  # vecteur des couts
    c2 = zeros(Int, nbvar)                  # deuxième vecteur des couts
    nb = zeros(Int, nbvar)
    for i in 1:nbvar
        flag = 1
        for valeur in split(readline(f))
            if flag == 1
                c1[i] = parse(Int, valeur)
                flag +=1
            elseif flag == 2
                c2[i] = parse(Int, valeur)
                flag +=1
            elseif flag == 3
                nb[i] = parse(Int, valeur)
                flag +=1
            else
                j = parse(Int, valeur)
                L[j,i] = 1
            end
        end
    end
    close(f)
    return nbvar, nbcontraintes, L, c1, c2
end

ext: jl

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine


text:
function compute_angle(origin,point)

    degree = atan((origin[2] - point[2]),(point[1] - origin[1])) * (180/pi)
    if degree > 90
        degree = 450 - degree
    else
        degree = 90 - degree
    end
    return degree
end

function compute_distance(origin,point)
    distance = sqrt((origin[1] - point[1])^2 + (origin[2] - point[2])^2)
    return distance
end

f = open("day10.txt");
ast_list = []
for (y, line) in enumerate(readlines(f))
    for (x, value) in enumerate(line)
        if value == '#'
            push!(ast_list,[x-1,y-1])
        end

    end
end
detected = -1
station = -1
for ast in ast_list
    detected_angles = Set()
    for other in ast_list
        if ast != other
            push!(detected_angles,compute_angle(ast,other))
        end
    end
    if length(detected_angles) > detected
        global detected
        global station
        detected = length(detected_angles)
        station = ast
    end
end

println(detected)
println(station)

filter!(x->x!=station,ast_list)

angle_dict = Dict()
for ast in ast_list
    angle = compute_angle(station,ast)

    if haskey(angle_dict, angle) == false
        push!(angle_dict,angle=>[ast])
    else
        value = angle_dict[angle]
        append!(value,[ast])
    end
end
iter = 0
for angle in sort([i for i in keys(angle_dict)],by=abs)
    remove = typemax(Float64)
    vaporized = nothing
    for i in angle_dict[angle]
        dist = compute_distance(station,i)
        if dist < remove
            remove = dist
            vaporized = i
        end
    end
    global iter
    iter += 1
    if iter == 200
        solution = vaporized[1]*100 + vaporized[2]
        println("Vaporize $vaporized, solution is $solution")
        break
    end
end


text:
"""
    Broadcaster

A broadcaster that calls all of the registered functions.

- `register`: add a function to the list of functions to be called.
- `broadcast`: call `f(adapter, msg_or_event)` for all of registered functions if `applicable`.
- `launch`: launch the adapter and start broadcasting messages.
"""
struct Broadcaster{Adapter<:ProtocolAdapter}
    adapter::Adapter
    methods::Vector{Any}
end
Broadcaster(adp::ProtocolAdapter) = Broadcaster(adp, [])

function register(f, b::Broadcaster)
    push!(b.methods, f)
end

struct ShutdownBroadcaster <: Exception end

"""
    broadcast(b::Broadcast, msg)

Call `f(adapter, msg_or_event)` for all of registered functions if `applicable`, return if `ShutdownBroadcaster` is thrown from any of the callback functions.
"""
function broadcast(b::Broadcaster, msg)
    args = b.adapter, msg
    for f in b.methods
        applicable(f, args...) || continue
        try
            f(args...)
        catch e
            if e isa ShutdownBroadcaster
                return false
            else
                @error "Error occurred executing broadcaster callbacks" exception = (e, catch_backtrace())
            end
        end
    end
    return true
end

function launch(b::Broadcaster, server, qq, key; check_version = true, kwargs...)
    @sync begin
        @async with_catch_backtrace(get_output_channel(b.adapter)) do
            MiraiBots.loop(b.adapter, server, qq, key; kwargs...)
        end
        receive_or_throw(b.adapter)
        check_version && @assert check_adapter_compatibility(b)
        @info "Adapter is connected to mirai."
        handle_message(b, Events.ProtocolAdapterConnected())
        for msg in get_output_channel(b.adapter)
            handle_message(b, msg)
        end
    end
end

function handle_message(bot, msg)
    if msg isa MiraiBots.ExceptionAndBacktrace
        @error "Error occurred in adapter" exception = (msg.exception, msg.backtrace)
    else
        broadcast(bot, msg) || close(bot.adapter)
    end
end

function check_adapter_compatibility(b)
    resp = send(b.adapter, Commands.about())
    is_adapter_compatibile(b.adapter, resp.data.version)
end

is_adapter_compatibile(::ProtocolAdapter, version) = version.major == 2


text:
abstract type Question end

mutable struct Stringq <: Question
    re::Regex
    label
    hint
end

"""
    stringq(re::Regex; label="", hint="")

Match string answer with regular expression

Arguments:

* `re`: a regular expression for grading
* `label`: optional label for the form element
* `hint`: optional plain-text hint that can be seen on hover

"""
stringq(re::Regex; label="", hint="") = Stringq(re, label, hint)

##
mutable struct Numericq <: Question
    val
    tol
    units
    label
    hint
end

"""
    numericq(value, tol=1e-3; label="", hint="", units="")

Match a numeric answer

Arguments:

* `value`: the numeric answer
* `tol`: ``|answer - value| \\le tol`` is used to determine correctness
* `label`: optional label for the form element
* `hint`: optional plain-text hint that can be seen on hover
* `units`: a string indicating expected units.

"""
function numericq(val, tol=1e-3, args...;
                  label="",
                  hint="", units::AbstractString="")

    Numericq(val, tol,  units, label, hint)
end

numericq(val::Int; kwargs...) = numericq(val, 0; kwargs...)

##
mutable struct Radioq <: Question
    choices
    answer
    values
    labels
    label
    hint
    inline
end

"""
    radioq(choices, answer; label="", hint="", keep_order=false)

Multiple choice question (one of several)

Arguments:

* `choices`: indexable collection of choices. As seen in the example, choices can be formatted with markdown.

* `answer::Int`: index of correct choice

* `keep_order::Boolean`: if `true` keeps display order of choices, otherwise they are shuffled.

* `inline::Bool`: hint to render inline (or not) if supported

* `label`: optional label for the form element

* `hint`: optional plain-text hint that can be seen on hover

Example:

```
choices = ["beta", raw"``\\beta``", "`beta`"]
answer = 2
radioq(choices, answer; hint="Which is the Greek symbol?")
```

"""
function radioq(choices, answer::Integer;
                label="", hint="", inline::Bool=(hint!=""),
                keep_order::Bool=false)
    inds = collect(1:length(choices))
    values = copy(inds)
    labels = choices
    !keep_order && shuffle!(inds)

    Radioq(choices[inds], findfirst(isequal(answer), inds),
           values, labels[inds], label, hint, inline)
end

##
mutable struct Multiq <: Question
    choices
    answer
    values
    labels
    label
    hint
    inline
end

"""
    multiq(choices, answers; label="", hint="", keep_order=false)

Multiple choice question (one *or more* of several)

Arguments:

* `choices`: indexable collection of choices. As seen in the example, choices can be formatted with markdown.

* `answers::Vector{Int}`: index of correct choice(s)

* `keep_order::Boolean`: if `true` keeps display order of choices, otherwise they are shuffled.

* `inline::Bool`: hint to render inline (or not) if supported

* `label`: optional label for the form element

* `hint`: optional plain-text hint that can be seen on hover

Example:

```
choices = ["pear", "tomato", "banana"]
answers = [1,3]
multiplecq(choices, answers; hint="not the red one!")
```

"""
function multiq(choices, answers;
                label="", hint="", inline::Bool=(hint!=""),
                keep_order::Bool=false)
    inds = collect(1:length(choices))
    values = copy(inds)
    labels = choices
    !keep_order && shuffle!(inds)

    Multiq(choices[inds], findall(in(answers), inds),
           values, labels[inds], label, hint, inline)
end

"""
    booleanq(ans; [label, hint])
True of false questions:

Example:

```
booleanq(true; label="Does it hurt...")
```

"""
function booleanq(ans::Bool;
                  labels::Vector=["true", "false"],
                  label="", hint::AbstractString="",
                  inline::Bool=true)
    choices = labels[1:2]
    ans = 2 - ans
    radioq(choices, ans;
           label=label, hint=hint, inline=inline, keep_order=true)
end

"""
    yesnoq(ans)

Boolean question with `yes` or `no` labels.

Examples:

```
yesnoq("yes")
yesnoq(true)
```

"""
yesnoq(ans::AbstractString, args...; kwargs...) = radioq(["Yes", "No"], ans == "yes" ? 1 : 2, args...; keep_order=true, kwargs...)
yesnoq(ans::Bool, args...; kwargs...) = yesnoq(ans ? "yes" : "no", args...;kwargs...)


text:
using BenchmarkTools

function me1(elementsInLeftChildBV, numerator, weight, current_meanl, current_meanr)
    dl = 0.0
    dr = 0.0
    f = weight  # simplified version
    for count in 1:length(f)
        @inbounds idx = count # f.parent.refs[count] + ooo
        @inbounds ni = numerator[count]
        @inbounds wi = weight[count]
        @inbounds eli = elementsInLeftChildBV[idx]
        if eli
            dl += (-log(ni / (wi * current_meanl)) - (ni - wi * current_meanl) / (wi * current_meanl))
        else
            dr += (-log(ni / (wi * current_meanr)) - (ni - wi * current_meanr) / (wi * current_meanr))
        end
    end
    return dl, dr
end

function me2(elementsInLeftChildBV, numerator, weight, current_meanl, current_meanr)
    dl = 0.0
    dr = 0.0
    f = weight  # simplified version
    for count in 1:length(f)
        @inbounds idx = count # f.parent.refs[count] + ooo
        @inbounds ni = numerator[count]
        @inbounds wi = weight[count]
        @inbounds eli = elementsInLeftChildBV[idx]
        wiTimesMean = ifelse(eli, wi * current_meanl, wi * current_meanr)
        addition = -log(ni / wiTimesMean) - (ni - wiTimesMean) / wiTimesMean
        if eli
            dl += addition
        else
            dr += addition
        end
    end
    return dl, dr
end

n = 2000000
eLR = (rand(n) .< .3)
nu = rand(n)
weight = rand(n)
mr = rand()
ml = rand()
me1(eLR,nu,weight,ml,mr)
me2(eLR,nu,weight,ml,mr)

@assert me1(eLR, nu, weight, ml, mr) == me2(eLR, nu, weight, ml, mr)
@benchmark me1($eLR, $nu, $weight, $ml, $mr)
@benchmark me2($eLR, $nu, $weight, $ml, $mr)
@btime me1($eLR, $nu, $weight, $ml, $mr)
@btime me2($eLR, $nu, $weight, $ml, $mr)
0
0
@code_warntype me2(eLR, nu, weight, ml, mr)


text:
@assert begin
    local_interp=InterpolationFunction(
        :(a0+a1*xi+a2*eta),
        [:xi, :eta],
        [:a0, :a1, :a2],
        [
            [0, 0]
        ]
    )

    abs_interp=InterpolationFunction(
        :(a0+a1*x+a2*y),
        [:x, :y],
        [:a0, :a1, :a2],
        [
            [0, 0]
        ]
    )

    dmn=Domain(
        local_interp,
        [
            (0.0, 1.0),
            (0.0, :(1.0-xi))
        ],
        [
            [0.0, 0.0],
            [1.0, 0.0],
            [0.0, 1.0]
        ]
    )

    fin=Finel(
        [
            [1.0, 1.0],
            [2.0, 1.0],
            [2.0, 2.0]
        ],
        dmn
    )

    vals=[
        0.0, 0.0, 1.0
    ]

    regr=finel_get_regr(fin, abs_interp)

    coefs=regress(regr, vals)

    isright=isapprox(coefs, [-1.0, 0.0, 1.0])

    Jac=Matrix{Float64}(undef, 2, 2)
    from_finel_jac!(Jac, fin, 0.0, 0.0)

    isright && isapprox(
        Jac,
        [1.0 1.0; 0.0 1.0]
    )
end


text:
"""
   dependenton(command::String, logicset::LogicalCombo; verbose=true)

Checks whether variables on the left of the infix operator (⊥) depend for their
distribution on variables on the right of the operator.

#### Operators

* x ⊥ y: checks if the distribution of potential values of x does not change when y changes.
**Alternatively:** x independentof y

* x !⊥ y: checks if the distribution of potential values of x changes when y changes.
**Alternatively:** x dependenton y

#### Examples - API

```
julia> logicset = logicalparse("a, b, c ∈ 1:4; a > b", verbose=false)

julia> dependenton("a ⊥ c", logicset)
True: a is independent of c
true

julia> dependenton("a ⊥ b", logicset)
False: a is dependent on b
range( a | b == b ) =   {}
––––––––––––––––––– – –––––––
range( a | b == 1 ) = {2,3,4}
range( a | b == 2 ) =  {3,4}
range( a | b == 3 ) =   {4}
false

julia> dependenton("a !⊥ b", logicset)
a is dependent on b
true
```

#### Examples - REPL

```
abstractlogic> a, b, c ∈ 1:4; a > b [clear]
Activeset Already Empty

a, b, c ∈ 1:4            Feasible Outcomes: 64   Perceived Outcomes: 64 ✓        :1 4 2
a > b                    Feasible Outcomes: 24   Perceived Outcomes: 36 ✓        :4 2 1

abstractlogic> a ⊥ c
True: a is independent of c

abstractlogic> a ⊥ b
False: a is dependent on b
range( a | b == b ) =   {}
––––––––––––––––––– – –––––––
range( a | b == 1 ) = {2,3,4}
range( a | b == 2 ) =  {3,4}
range( a | b == 3 ) =   {4}

abstractlogic> a !⊥ b
True: a is dependent on b
```
"""
function dependenton(command::String, logicset::LogicalCombo; verbose=true)

    m = match(r"^(.*?)([⊥!]{1,2}|(dependenton|independentof))(.*?)$",command)
    left, operator, right  = m.captures[[1,2,4]] .|> strip

    ranges = range(logicset)

    rightranges = ranges[Symbol(right)]
    (length(rightranges) == 1) && return println("$left cannot vary in $right if right is constant!")

    leftranges = ranges[Symbol(left)]
    (length(rightranges) == 1) && return println("$left cannot vary in $right if left is constant!")

    leftsets =
      [range(logicalparse("$right = $v", logicset=logicset, verbose=false))[Symbol(left)]
      for v in rightranges ]

    independence = all(length.(leftsets) .== length(leftranges))
    indeptext = (independence ? "independent of" : "dependent on")

    (operator == "⊥")  && independence &&
      return (verbose && println("True: $left is $indeptext $right"); true)
    (operator == "!⊥") && !independence &&
      return (verbose && println("True: $left is $indeptext $right"); true)

   if verbose
       println("False: $left is $indeptext $right")
       txtout = "| " * join(["range( $left \\| $right == $right )", " = ", "{}"], " | ") * " | \n"
       txtout *= "| " * join(fill(":---:", 3), " | ") * " | \n"
       for i in 1:length(rightranges)
         txtout *= "| " * "range( $left \\| $right == $(rightranges[i]) ) " * " | = | {" *
           join(leftsets[i], ",") * "} | \n"
       end
       printmarkdown(txtout)
       println()
   end

   return false
end


text:
using HDF5

TRAIN_DATA_PATH="datasets/train_catvnoncat.h5"
TEST_DATA_PATH="datasets/test_catvnoncat.h5"

train_data_x_orig = h5read(TRAIN_DATA_PATH, "train_set_x")
train_data_y_orig = h5read(TRAIN_DATA_PATH, "train_set_y")

test_data_x_orig = h5read(TEST_DATA_PATH, "test_set_x")
test_data_y_orig = h5read(TEST_DATA_PATH, "test_set_y")

classes = h5read(TEST_DATA_PATH, "list_classes")

train_data_y_orig = reshape(train_data_y_orig, (1, size(train_data_y_orig)[1]))

test_data_y_orig = reshape(test_data_y_orig, (1, size(test_data_y_orig)[1]))

function reshape_image_array(image_array::Array{UInt8,4})
    return reshape(image_array, :, size(image_array)[4])
end
function reshape_image_array_reverse(image_array::Array{UInt8,2})
    return reshape(image_array, 64, 64, 3, :)
end

test_data_x = reshape_image_array(test_data_x_orig) ./ 255
test_data_x = Array{Float32, 2}(test_data_x)

train_data_x = reshape_image_array(train_data_x_orig) ./ 255
train_data_x = Array{Float32, 2}(train_data_x)

using Statistics, BenchmarkTools

include("../../NeuralNetwork.jl")

layer_dims = [12288, 7, 1]

# @btime neural_network_dense($train_data_x, $train_data_y_orig, $layer_dims, 100, 0.0075)

parameters, activations = neural_network_dense(train_data_x, train_data_y_orig, layer_dims, 1000, 0.0075)

predicts, accuracy = predict(train_data_x, train_data_y_orig, parameters, activations)

predicts, accuracy = predict(test_data_x, test_data_y_orig, parameters, activations)


text:
#=
Draw arrow at point (x_b,y_b) with length m at angle a
=#

function draw_arrow(ctx,x_b,y_b,m,a,mycolor=(1,1,1))

    set_source_rgb(ctx,mycolor...)

    move_to(ctx,x_b,y_b)

    x_f = m * cos(a)
    y_f = m * sin(a)
    line_to(ctx,x_f,y_f)

    stroke(ctx)
end


text:
push!(LOAD_PATH, "../LF_common/")
using m_LF3d
include("../LF_common/sparse_LF3d.jl")
include("build_nabla2_x.jl")
include("build_nabla2_y.jl")
include("build_nabla2_z.jl")

const ⊗ = kron

function test_main( NN::Array{Int64} )

  println("NN = ", NN)

  Nx = NN[1]
  Ny = NN[2]
  Nz = NN[3]
  Npoints = Nx*Ny*Nz

  AA = [0.0, 0.0, 0.0]
  BB = [16.0, 16.0, 16.0]

  # Initialize LF
  LF = init_LF3d_c( NN, AA, BB, verbose=true )

  @time nabla2_x = build_nabla2_x(LF)
  @time nabla2_y = build_nabla2_y(LF)
  @time nabla2_z = build_nabla2_z(LF)

  @time combine_nabla2( NN, nabla2_x, nabla2_y, nabla2_z )

end

function combine_nabla2( NN, nabla2_x, nabla2_y, nabla2_z )

  Nx = NN[1]
  Ny = NN[2]
  Nz = NN[3]
  Npoints = Nx*Ny*Nz

  # number of non-zero per-column
  nnz = Nx + Ny + Nz - 2

  # combined results
  rowval = zeros(Int64,nnz*Npoints)
  colptr = zeros(Int64,Npoints+1)
  colptr[1] = 1
  nzval = zeros(Float64,nnz*Npoints)

  # shortcuts
  x_rowval = nabla2_x.rowval
  y_rowval = nabla2_y.rowval
  z_rowval = nabla2_z.rowval
  #
  x_nzval = nabla2_x.nzval
  y_nzval = nabla2_y.nzval
  z_nzval = nabla2_z.nzval

  c_rowval = zeros(Int64, Nx+Ny+Nz)  # temporary combined indices of row
  c_nzval  = zeros(Float64, Nx+Ny+Nz)

  idxSorted = zeros(Int64, Nx+Ny+Nz)

  idxGbl = 0

  println(size(x_rowval))

  compxyz = Array{ASCIIString}(Nx+Ny+Nz)
  compxyz[1:Nx] = "x"
  compxyz[Nx+1:Nx+Ny] = "y"
  compxyz[Nx+Ny+1:Nx+Ny+Nz] = "z"

  idxDiag = 0
  for colGbl = 1:Npoints
    #
    idxDiag = idxDiag + 1
    #println("\nidxDiag = ", idxDiag)
    #
    idx = 0
    for i in range(1,Nx) + Nx*(colGbl-1)
      idx = idx + 1
      c_rowval[idx] = x_rowval[i]
      c_nzval[idx] = x_nzval[i]
    end
    for j in range(1,Ny) + Ny*(colGbl-1)
      idx = idx + 1
      c_rowval[idx] = y_rowval[j]
      c_nzval[idx] = y_nzval[j]
    end
    for k in range(1,Nz) + Nz*(colGbl-1)
      idx = idx + 1
      c_rowval[idx] = z_rowval[k]
      c_nzval[idx] = z_nzval[k]
    end
    # Sort the indices
    idxSorted = sortperm(c_rowval)
    println(compxyz[idxSorted])
    #println(c_rowval[idxSorted])
    #
  end

end

test_main( [3, 2, 4] )
#test_main( [40, 40, 40] )
#test_main( [80, 80, 80] )


text:
export SU2

"""
    Represents SU(2).

fields: u and v.
"""
struct SU2
    a::Real
    v::Complex
end


text:
include("common.jl")
@testset "Lazy interpolants" begin
    ## simple problems
    @testset "simple problems" begin
        prob_inplace = prob_dde_1delay
        prob_notinplace = prob_dde_1delay_scalar_notinplace

        # Vern6
        @testset "Vern6" begin
            sol = solve(prob_inplace, MethodOfSteps(Vern6()))

            @test sol.errors[:l∞] < 8.7e-4
            @test sol.errors[:final] < 7.5e-6
            @test sol.errors[:l2] < 5.3e-4

            sol2 = solve(prob_notinplace, MethodOfSteps(Vern6()))

            @test sol.t ≈ sol2.t && sol[1, :] ≈ sol2.u
        end

        # Vern7
        @testset "Vern7" begin
            sol = solve(prob_inplace, MethodOfSteps(Vern7()))

            @test sol.errors[:l∞] < 4.0e-4
            @test sol.errors[:final] < 3.5e-7
            @test sol.errors[:l2] < 1.9e-4

            sol2 = solve(prob_notinplace, MethodOfSteps(Vern7()))

            @test sol.t ≈ sol2.t && sol[1, :] ≈ sol2.u
        end

        # Vern8
        @testset "Vern8" begin
            sol = solve(prob_inplace, MethodOfSteps(Vern8()))

            @test sol.errors[:l∞] < 2.0e-3
            @test sol.errors[:final] < 1.8e-5
            @test sol.errors[:l2] < 9.0e-4

            sol2 = solve(prob_notinplace, MethodOfSteps(Vern8()))

            @test sol.t ≈ sol2.t && sol[1, :] ≈ sol2.u
        end

        # Vern9
        @testset "Vern9" begin
            sol = solve(prob_inplace, MethodOfSteps(Vern9()))

            @test sol.errors[:l∞] < 1.5e-3
            @test sol.errors[:final] < 3.8e-6
            @test sol.errors[:l2] < 6.2e-4

            sol2 = solve(prob_notinplace, MethodOfSteps(Vern9()))

            @test sol.t ≈ sol2.t && sol[1, :] ≈ sol2.u
        end
    end

    # model of Mackey and Glass
    @testset "Mackey and Glass" begin
        prob = prob_dde_mackey

        # Vern6
        solve(prob, MethodOfSteps(Vern6()))

        # Vern7
        solve(prob, MethodOfSteps(Vern7()))

        # Vern8
        solve(prob, MethodOfSteps(Vern8()))

        # Vern9
        solve(prob, MethodOfSteps(Vern9()))
    end
end


text:
using Test
using LinearAlgebra
using CartesianMesh
using PolynomialBasis
using ImplicitDomainQuadrature
using HDGElasticity

function allapprox(v1,v2)
    return all(v1 .≈ v2)
end

function allapprox(v1,v2,atol)
    @assert length(v1) == length(v2)
    return all([isapprox(v1[i],v2[i],atol=atol) for i = 1:length(v1)])
end

mesh = UniformMesh([0.,0.],[2.,1.],[1,1])
poly = InterpolatingPolynomial(1,2,1)
coeffs = reshape(ones(4),4,1)
dgmesh = HDGElasticity.DGMesh(mesh,coeffs,poly)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,1,2,coeffs,poly)

cellmap = HDGElasticity.CellMap(dgmesh.domain[1])
Dhalf = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(1.,2.))
stabilization = 0.1
normals = HDGElasticity.reference_normals()

lop = HDGElasticity.local_operator(ufs.vbasis,ufs.vquads[1,1],
    ufs.fquads[1,1],dgmesh.facemaps,Dhalf,stabilization,cellmap)
facelhops = HDGElasticity.local_hybrid_operator(ufs.vbasis,ufs.sbasis,
    ufs.fquads[1,1],dgmesh.facemaps,normals,Dhalf,stabilization,cellmap)
lhop = hcat(facelhops...)

function bc_displacement(coords;alpha=0.1,beta=0.1)
    disp = copy(coords)
    disp[1,:] .*= alpha
    disp[2,:] .*= beta
    return disp
end

function compute_rhs(lhop,disp,ndofs)
    H = vcat(disp...)
    return lhop*H
end

H1c = [0. 2.
       0. 0.]
H2c = [2. 2.
       0. 1.]
H3c = [0. 2.
       1. 1.]
H4c = [0. 0.
       0. 1.]
Hcoords = [H1c,H2c,H3c,H4c]
Hdisp = vec.(bc_displacement.(Hcoords))

rhs = compute_rhs(lhop,Hdisp,20)
sol = lop\rhs

L = -Dhalf*reshape(sol[1:12],3,:)
U = reshape(sol[13:20],2,:)

testL = zeros(3,4)
testL[1:2,:] .= 0.6
testU = zeros(2,4)
testU[:,2] .= [0.,0.1]
testU[:,3] .= [0.2,0.]
testU[:,4] .= [0.2,0.1]

@test allapprox(L,testL,1e-14)
@test allapprox(U,testU,1e-14)

function plane_distance_function(coords,n,x0)
    return [n'*(coords[:,idx]-x0) for idx in 1:size(coords)[2]]
end

coords = HDGElasticity.nodal_coordinates(mesh,poly.basis)
coeffs = reshape(plane_distance_function(coords,[1.,1.]/sqrt(2.),[0.5,0.]),4,1)
dgmesh = HDGElasticity.DGMesh(mesh,coeffs,poly)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,1,4,coeffs,poly)

update!(poly,coeffs[:,1])
update!(ufs.imap,ufs.icoeffs[1])

lop = HDGElasticity.local_operator(ufs.vbasis,ufs.vquads[1,1],
    ufs.fquads[1,1],dgmesh.facemaps,ufs.iquad,ufs.imap,-ufs.inormals[1],
    Dhalf,stabilization,cellmap)
facelhops = HDGElasticity.local_hybrid_operator(ufs.vbasis,ufs.sbasis,
    ufs.fquads[1,1],dgmesh.facemaps,normals,Dhalf,
    stabilization,cellmap)
lhop = hcat(facelhops...)
ilhop = HDGElasticity.local_hybrid_operator_on_interface(ufs.vbasis,ufs.sbasis,
    ufs.iquad,ufs.imap,-ufs.inormals[1],Dhalf,stabilization,cellmap)

HIrefcoords = reshape(ufs.icoeffs[1],2,:)
HIcoords = hcat([cellmap(HIrefcoords[:,i]) for i in 1:size(HIrefcoords)[2]]...)
HIdisp = vec(bc_displacement(HIcoords))

rhs = compute_rhs(lhop,Hdisp,20)
rI = ilhop*HIdisp
rhs2 = rhs+rI

sol = lop\rhs2
L = -Dhalf*reshape(sol[1:12],3,:)
U = reshape(sol[13:20],2,:)

testL = zeros(3,4)
testL[1:2,:] .= 0.6
testU = zeros(2,4)
testU[:,2] .= [0.,0.1]
testU[:,3] .= [0.2,0.]
testU[:,4] .= [0.2,0.1]

@test allapprox(L,testL,1e-14)
@test allapprox(U,testU,1e-14)

lop = HDGElasticity.local_operator(ufs.vbasis,ufs.vquads[2,1],
    ufs.fquads[2,1],dgmesh.facemaps,ufs.iquad,ufs.imap,ufs.inormals[1],
    Dhalf,stabilization,cellmap)
facelhops = HDGElasticity.local_hybrid_operator(ufs.vbasis,ufs.sbasis,
    ufs.fquads[2,1],dgmesh.facemaps,normals,Dhalf,
    stabilization,cellmap)
lhop = hcat(facelhops...)
ilhop = HDGElasticity.local_hybrid_operator_on_interface(ufs.vbasis,ufs.sbasis,
    ufs.iquad,ufs.imap,ufs.inormals[1],Dhalf,stabilization,cellmap)

lochyb = hcat(lhop,ilhop)

isactiveface = [length(fq) > 0 ? true : false for fq in ufs.fquads[2,1]]
faceids = findall(isactiveface)
Hfacedisp = vcat(vec.(bc_displacement.(Hcoords[faceids]))...)
Hdisp = vcat(Hfacedisp,HIdisp)

rhs = lochyb*Hdisp

sol = lop\rhs
L = -Dhalf*reshape(sol[1:12],3,:)
U = reshape(sol[13:20],2,:)

@test allapprox(L,testL,1e-12)
@test allapprox(U,testU,1e-12)

Hdisp = vcat(vec.(bc_displacement.(Hcoords[faceids],beta=0.2))...)
HIdisp = vec(bc_displacement(HIcoords,beta=0.2))
H = vcat(Hdisp,HIdisp)

rhs = lochyb*H

sol = lop\rhs
L = -Dhalf*reshape(sol[1:12],3,:)
U = reshape(sol[13:20],2,:)

testL = zeros(3,4)
testL[1,:] .= 0.7
testL[2,:] .= 1.1
testU = zeros(2,4)
testU[:,2] .= [0.,0.2]
testU[:,3] .= [0.2,0.]
testU[:,4] .= [0.2,0.2]

@test allapprox(L,testL,1e-12)
@test allapprox(U,testU,1e-12)

function shear_bc_displacement(coords;alpha=0.1)
    disp = similar(coords)
    disp[1,:] = alpha*coords[2,:]
    disp[2,:] .= 0.0
    return disp
end

Hdisp = vcat(vec.(shear_bc_displacement.(Hcoords[faceids]))...)
HIdisp = vec(shear_bc_displacement(HIcoords))

H = vcat(Hdisp,HIdisp)

rhs = lochyb*H

sol = lop\rhs
L = -Dhalf*reshape(sol[1:12],3,:)
U = reshape(sol[13:20],2,:)

testL = zeros(3,4)
testL[3,:] .= 0.2
testU = zeros(2,4)
testU[:,2] .= [0.1,0.]
testU[:,3] .= [0.,0.]
testU[:,4] .= [0.1,0.]

@test allapprox(L,testL,1e-12)
@test allapprox(U,testU,1e-12)


text:
import Multirate
import Multirate: NaiveResamplers
using Base.Test

Tx           = Float32
Th           = Float32
numFilters   = 32
x            = Array{Tx}(1:101)
resampleRate = 33/32

cutoffFreq      = 0.45
transitionWidth = 0.05
(hLen, β)       = Multirate.kaiserlength( transitionWidth, samplerate = numFilters )
hLen            = ceil( Int, hLen/numFilters ) .* numFilters
h               = Multirate.firdes( hLen, cutoffFreq, DSP.kaiser, samplerate = numFilters, beta = β ) .* numFilters
h               = convert( Vector{Th}, h)

@time yNaive = NaiveResamplers.naivefilt( h, x, resampleRate, numFilters )
@time yArb   = Multirate.filt( h, x, resampleRate )

commonLen = min( length(yNaive), length(yArb) )

yNaiveCommon = yNaive[1:commonLen]
yArbCommon = yArb[1:commonLen]

isapprox( yNaiveCommon, yArbCommon, atol=eps() ) ||  display( [ (1:commonLen) yNaiveCommon yArbCommon (yNaiveCommon.-yArbCommon)] )


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: textract
using AWS.Compat
using AWS.UUIDs

"""
    AnalyzeDocument()

Analyzes an input document for relationships between detected items.  The types of information returned are as follows:    Form data (key-value pairs). The related information is returned in two Block objects, each of type KEY_VALUE_SET: a KEY Block object and a VALUE Block object. For example, Name: Ana Silva Carolina contains a key and value. Name: is the key. Ana Silva Carolina is the value.   Table and table cell data. A TABLE Block object contains information about a detected table. A CELL Block object is returned for each cell in a table.   Lines and words of text. A LINE Block object contains one or more WORD Block objects. All lines and words that are detected in the document are returned (including text that doesn't have a relationship with the value of FeatureTypes).    Selection elements such as check boxes and option buttons (radio buttons) can be detected in form data and in tables. A SELECTION_ELEMENT Block object contains information about a selection element, including the selection status. You can choose which type of analysis to perform by specifying the FeatureTypes list.  The output is returned in a list of Block objects.  AnalyzeDocument is a synchronous operation. To analyze documents asynchronously, use StartDocumentAnalysis. For more information, see Document Text Analysis.

# Required Parameters
- `Document`: The input document as base64-encoded bytes or an Amazon S3 object. If you use the AWS CLI to call Amazon Textract operations, you can't pass image bytes. The document must be an image in JPEG or PNG format. If you're using an AWS SDK to call Amazon Textract, you might not need to base64-encode image bytes that are passed using the Bytes field.
- `FeatureTypes`: A list of the types of analysis to perform. Add TABLES to the list to return information about the tables that are detected in the input document. Add FORMS to return detected form data. To perform both types of analysis, add TABLES and FORMS to FeatureTypes. All lines and words detected in the document are included in the response (including text that isn't related to the value of FeatureTypes).

# Optional Parameters
- `HumanLoopConfig`: Sets the configuration for the human in the loop workflow for analyzing documents.
"""
analyze_document(Document, FeatureTypes; aws_config::AWSConfig=global_aws_config()) = textract("AnalyzeDocument", Dict{String, Any}("Document"=>Document, "FeatureTypes"=>FeatureTypes); aws_config=aws_config)
analyze_document(Document, FeatureTypes, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = textract("AnalyzeDocument", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Document"=>Document, "FeatureTypes"=>FeatureTypes), args)); aws_config=aws_config)

"""
    DetectDocumentText()

Detects text in the input document. Amazon Textract can detect lines of text and the words that make up a line of text. The input document must be an image in JPEG or PNG format. DetectDocumentText returns the detected text in an array of Block objects.  Each document page has as an associated Block of type PAGE. Each PAGE Block object is the parent of LINE Block objects that represent the lines of detected text on a page. A LINE Block object is a parent for each word that makes up the line. Words are represented by Block objects of type WORD.  DetectDocumentText is a synchronous operation. To analyze documents asynchronously, use StartDocumentTextDetection. For more information, see Document Text Detection.

# Required Parameters
- `Document`: The input document as base64-encoded bytes or an Amazon S3 object. If you use the AWS CLI to call Amazon Textract operations, you can't pass image bytes. The document must be an image in JPEG or PNG format. If you're using an AWS SDK to call Amazon Textract, you might not need to base64-encode image bytes that are passed using the Bytes field.

"""
detect_document_text(Document; aws_config::AWSConfig=global_aws_config()) = textract("DetectDocumentText", Dict{String, Any}("Document"=>Document); aws_config=aws_config)
detect_document_text(Document, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = textract("DetectDocumentText", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Document"=>Document), args)); aws_config=aws_config)

"""
    GetDocumentAnalysis()

Gets the results for an Amazon Textract asynchronous operation that analyzes text in a document. You start asynchronous text analysis by calling StartDocumentAnalysis, which returns a job identifier (JobId). When the text analysis operation finishes, Amazon Textract publishes a completion status to the Amazon Simple Notification Service (Amazon SNS) topic that's registered in the initial call to StartDocumentAnalysis. To get the results of the text-detection operation, first check that the status value published to the Amazon SNS topic is SUCCEEDED. If so, call GetDocumentAnalysis, and pass the job identifier (JobId) from the initial call to StartDocumentAnalysis.  GetDocumentAnalysis returns an array of Block objects. The following types of information are returned:    Form data (key-value pairs). The related information is returned in two Block objects, each of type KEY_VALUE_SET: a KEY Block object and a VALUE Block object. For example, Name: Ana Silva Carolina contains a key and value. Name: is the key. Ana Silva Carolina is the value.   Table and table cell data. A TABLE Block object contains information about a detected table. A CELL Block object is returned for each cell in a table.   Lines and words of text. A LINE Block object contains one or more WORD Block objects. All lines and words that are detected in the document are returned (including text that doesn't have a relationship with the value of the StartDocumentAnalysis FeatureTypes input parameter).    Selection elements such as check boxes and option buttons (radio buttons) can be detected in form data and in tables. A SELECTION_ELEMENT Block object contains information about a selection element, including the selection status. Use the MaxResults parameter to limit the number of blocks that are returned. If there are more results than specified in MaxResults, the value of NextToken in the operation response contains a pagination token for getting the next set of results. To get the next page of results, call GetDocumentAnalysis, and populate the NextToken request parameter with the token value that's returned from the previous call to GetDocumentAnalysis. For more information, see Document Text Analysis.

# Required Parameters
- `JobId`: A unique identifier for the text-detection job. The JobId is returned from StartDocumentAnalysis. A JobId value is only valid for 7 days.

# Optional Parameters
- `MaxResults`: The maximum number of results to return per paginated call. The largest value that you can specify is 1,000. If you specify a value greater than 1,000, a maximum of 1,000 results is returned. The default value is 1,000.
- `NextToken`: If the previous response was incomplete (because there are more blocks to retrieve), Amazon Textract returns a pagination token in the response. You can use this pagination token to retrieve the next set of blocks.
"""
get_document_analysis(JobId; aws_config::AWSConfig=global_aws_config()) = textract("GetDocumentAnalysis", Dict{String, Any}("JobId"=>JobId); aws_config=aws_config)
get_document_analysis(JobId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = textract("GetDocumentAnalysis", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("JobId"=>JobId), args)); aws_config=aws_config)

"""
    GetDocumentTextDetection()

Gets the results for an Amazon Textract asynchronous operation that detects text in a document. Amazon Textract can detect lines of text and the words that make up a line of text. You start asynchronous text detection by calling StartDocumentTextDetection, which returns a job identifier (JobId). When the text detection operation finishes, Amazon Textract publishes a completion status to the Amazon Simple Notification Service (Amazon SNS) topic that's registered in the initial call to StartDocumentTextDetection. To get the results of the text-detection operation, first check that the status value published to the Amazon SNS topic is SUCCEEDED. If so, call GetDocumentTextDetection, and pass the job identifier (JobId) from the initial call to StartDocumentTextDetection.  GetDocumentTextDetection returns an array of Block objects.  Each document page has as an associated Block of type PAGE. Each PAGE Block object is the parent of LINE Block objects that represent the lines of detected text on a page. A LINE Block object is a parent for each word that makes up the line. Words are represented by Block objects of type WORD. Use the MaxResults parameter to limit the number of blocks that are returned. If there are more results than specified in MaxResults, the value of NextToken in the operation response contains a pagination token for getting the next set of results. To get the next page of results, call GetDocumentTextDetection, and populate the NextToken request parameter with the token value that's returned from the previous call to GetDocumentTextDetection. For more information, see Document Text Detection.

# Required Parameters
- `JobId`: A unique identifier for the text detection job. The JobId is returned from StartDocumentTextDetection. A JobId value is only valid for 7 days.

# Optional Parameters
- `MaxResults`: The maximum number of results to return per paginated call. The largest value you can specify is 1,000. If you specify a value greater than 1,000, a maximum of 1,000 results is returned. The default value is 1,000.
- `NextToken`: If the previous response was incomplete (because there are more blocks to retrieve), Amazon Textract returns a pagination token in the response. You can use this pagination token to retrieve the next set of blocks.
"""
get_document_text_detection(JobId; aws_config::AWSConfig=global_aws_config()) = textract("GetDocumentTextDetection", Dict{String, Any}("JobId"=>JobId); aws_config=aws_config)
get_document_text_detection(JobId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = textract("GetDocumentTextDetection", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("JobId"=>JobId), args)); aws_config=aws_config)

"""
    StartDocumentAnalysis()

Starts the asynchronous analysis of an input document for relationships between detected items such as key-value pairs, tables, and selection elements.  StartDocumentAnalysis can analyze text in documents that are in JPEG, PNG, and PDF format. The documents are stored in an Amazon S3 bucket. Use DocumentLocation to specify the bucket name and file name of the document.   StartDocumentAnalysis returns a job identifier (JobId) that you use to get the results of the operation. When text analysis is finished, Amazon Textract publishes a completion status to the Amazon Simple Notification Service (Amazon SNS) topic that you specify in NotificationChannel. To get the results of the text analysis operation, first check that the status value published to the Amazon SNS topic is SUCCEEDED. If so, call GetDocumentAnalysis, and pass the job identifier (JobId) from the initial call to StartDocumentAnalysis. For more information, see Document Text Analysis.

# Required Parameters
- `DocumentLocation`: The location of the document to be processed.
- `FeatureTypes`: A list of the types of analysis to perform. Add TABLES to the list to return information about the tables that are detected in the input document. Add FORMS to return detected form data. To perform both types of analysis, add TABLES and FORMS to FeatureTypes. All lines and words detected in the document are included in the response (including text that isn't related to the value of FeatureTypes).

# Optional Parameters
- `ClientRequestToken`: The idempotent token that you use to identify the start request. If you use the same token with multiple StartDocumentAnalysis requests, the same JobId is returned. Use ClientRequestToken to prevent the same job from being accidentally started more than once. For more information, see Calling Amazon Textract Asynchronous Operations.
- `JobTag`: An identifier that you specify that's included in the completion notification published to the Amazon SNS topic. For example, you can use JobTag to identify the type of document that the completion notification corresponds to (such as a tax form or a receipt).
- `KMSKeyId`: The KMS key used to encrypt the inference results. This can be in either Key ID or Key Alias format. When a KMS key is provided, the KMS key will be used for server-side encryption of the objects in the customer bucket. When this parameter is not enabled, the result will be encrypted server side,using SSE-S3.
- `NotificationChannel`: The Amazon SNS topic ARN that you want Amazon Textract to publish the completion status of the operation to.
- `OutputConfig`: Sets if the output will go to a customer defined bucket. By default, Amazon Textract will save the results internally to be accessed by the GetDocumentAnalysis operation.
"""
start_document_analysis(DocumentLocation, FeatureTypes; aws_config::AWSConfig=global_aws_config()) = textract("StartDocumentAnalysis", Dict{String, Any}("DocumentLocation"=>DocumentLocation, "FeatureTypes"=>FeatureTypes); aws_config=aws_config)
start_document_analysis(DocumentLocation, FeatureTypes, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = textract("StartDocumentAnalysis", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("DocumentLocation"=>DocumentLocation, "FeatureTypes"=>FeatureTypes), args)); aws_config=aws_config)

"""
    StartDocumentTextDetection()

Starts the asynchronous detection of text in a document. Amazon Textract can detect lines of text and the words that make up a line of text.  StartDocumentTextDetection can analyze text in documents that are in JPEG, PNG, and PDF format. The documents are stored in an Amazon S3 bucket. Use DocumentLocation to specify the bucket name and file name of the document.   StartTextDetection returns a job identifier (JobId) that you use to get the results of the operation. When text detection is finished, Amazon Textract publishes a completion status to the Amazon Simple Notification Service (Amazon SNS) topic that you specify in NotificationChannel. To get the results of the text detection operation, first check that the status value published to the Amazon SNS topic is SUCCEEDED. If so, call GetDocumentTextDetection, and pass the job identifier (JobId) from the initial call to StartDocumentTextDetection. For more information, see Document Text Detection.

# Required Parameters
- `DocumentLocation`: The location of the document to be processed.

# Optional Parameters
- `ClientRequestToken`: The idempotent token that's used to identify the start request. If you use the same token with multiple StartDocumentTextDetection requests, the same JobId is returned. Use ClientRequestToken to prevent the same job from being accidentally started more than once. For more information, see Calling Amazon Textract Asynchronous Operations.
- `JobTag`: An identifier that you specify that's included in the completion notification published to the Amazon SNS topic. For example, you can use JobTag to identify the type of document that the completion notification corresponds to (such as a tax form or a receipt).
- `KMSKeyId`: The KMS key used to encrypt the inference results. This can be in either Key ID or Key Alias format. When a KMS key is provided, the KMS key will be used for server-side encryption of the objects in the customer bucket. When this parameter is not enabled, the result will be encrypted server side,using SSE-S3.
- `NotificationChannel`: The Amazon SNS topic ARN that you want Amazon Textract to publish the completion status of the operation to.
- `OutputConfig`: Sets if the output will go to a customer defined bucket. By default Amazon Textract will save the results internally to be accessed with the GetDocumentTextDetection operation.
"""
start_document_text_detection(DocumentLocation; aws_config::AWSConfig=global_aws_config()) = textract("StartDocumentTextDetection", Dict{String, Any}("DocumentLocation"=>DocumentLocation); aws_config=aws_config)
start_document_text_detection(DocumentLocation, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = textract("StartDocumentTextDetection", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("DocumentLocation"=>DocumentLocation), args)); aws_config=aws_config)


text:
module System
    export libmpi, mpiexec
    using Preferences, Libdl
    const libmpi = @load_preference("libmpi")
    const mpiexec_path = @load_preference("mpiexec")
    mpiexec(;adjust_PATH=true, adjust_LIBPATH=true) = `$mpiexec_path`
    mpiexec(f;adjust_PATH=true, adjust_LIBPATH=true) = f(`$mpiexec_path`)

    libmpi_handle = C_NULL
    function __init__()
        global libmpi_handle = Libdl.dlopen(libmpi, Libdl.RTLD_LAZY | Libdl.RTLD_GLOBAL)
    end
end


text:
using Sophus
using Test
using StaticArrays

@testset "Sophus.jl" begin
    # Write your own tests here.
    so = SO2(Float16(3.2))
    @test typeof(so.tangent) == Float16
    @test so.DoF == Int16(1)
    @test generator(so) == SMatrix{2,2,Float16}(0, -1, 1, 0)
    @test matrix(so) == SMatrix{2,2,Float16}(0, -3.2, 3.2, 0)
    @test transformation(so) == exp.(SMatrix{2,2,Float16}(0, -3.2, 3.2, 0))
end


text:
function _pad_spaces(slen::Int, maxlen::Int)
    " "^max(0, maxlen - slen + 1)
end

_pad_spaces(str::String, maxlen::Int) = _pad_spaces(length(str), maxlen)

function Base.show(io::IO, ::MIME"text/plain", flux_res::FluxSummary)
    longest_biomass_len =
        maximum(length(k) for k in keys(flux_res.biomass_fluxes); init = 0)
    longest_import_len = maximum(length(k) for k in keys(flux_res.import_fluxes); init = 0)
    longest_export_len = maximum(length(k) for k in keys(flux_res.export_fluxes); init = 0)

    if !isempty(flux_res.unbounded_fluxes)
        longest_unbounded_len =
            maximum([length(k) for k in keys(flux_res.unbounded_fluxes)])
        word_pad = max(
            longest_biomass_len,
            longest_export_len,
            longest_import_len,
            longest_unbounded_len,
        )
    else
        word_pad = max(longest_biomass_len, longest_export_len, longest_import_len)
    end

    println(io, "Biomass")
    for (k, v) in flux_res.biomass_fluxes
        println(io, "  ", k, ":", _pad_spaces(k, word_pad), round(v, digits = 4))
    end

    println(io, "Import")
    for (k, v) in flux_res.import_fluxes
        println(io, "  ", k, ":", _pad_spaces(k, word_pad), round(v, digits = 4))
    end

    println(io, "Export")
    for (k, v) in flux_res.export_fluxes
        println(io, "  ", k, ":", _pad_spaces(k, word_pad), round(v, digits = 4))
    end

    if !isempty(flux_res.unbounded_fluxes)
        println(io, "Unbounded")
        for (k, v) in flux_res.unbounded_fluxes
            println(io, "  ", k, ":", _pad_spaces(k, word_pad), round(v, digits = 4))
        end
    end
    return nothing
end


text:
fsabs(x) = sum(abs(x))
fsumsq(x) = sum(x.^2)

facts("FixedDimProblem") do
  context("1-dimensional, single-objective sumabs") do
    ss = symmetric_search_space(1)
    p = BlackBoxOptim.FixedDimProblem("sumabs", [fsabs], ss, [0.0])

    @fact is_fixed_dimensional(p)         => true
    @fact is_any_dimensional(p)           => false
    @fact is_single_objective_problem(p)  => true
    @fact is_multi_objective_problem(p)   => false
    @fact numdims(p)                      => 1
    @fact search_space(p)                 => ss

    @fact eval1([0.0], p)                 => 0.0
    @fact eval1([1.2], p)                 => 1.2
    @fact eval1([-1.9], p)                => 1.9
  end

  context("3-dimensional, single-objective sumabs") do
    ss = symmetric_search_space(3)
    p = BlackBoxOptim.FixedDimProblem("sumabs", [fsabs], ss, [0.0])

    @fact is_fixed_dimensional(p)         => true
    @fact is_any_dimensional(p)           => false
    @fact is_single_objective_problem(p)  => true
    @fact is_multi_objective_problem(p)   => false
    @fact numdims(p)                      => 3
    @fact search_space(p)                 => ss

    @fact eval1([0.0, 1.0, 2.0], p)       => 3.0
    @fact eval1([-1.0, 1.0, 2.0], p)      => 4.0
  end

  context("1-dimensional, multi-objective sumabs_sumsq") do
    ss = symmetric_search_space(1)
    p = BlackBoxOptim.FixedDimProblem("sumabs_sumsq", [fsabs, fsumsq], ss, [0.0, 0.0])

    @fact is_fixed_dimensional(p)         => true
    @fact is_any_dimensional(p)           => false
    @fact is_single_objective_problem(p)  => false
    @fact is_multi_objective_problem(p)   => true
    @fact numdims(p)                      => 1
    @fact search_space(p)                 => ss

    @fact eval1([0.0], p)                 => 0.0
    @fact eval1([1.2], p)                 => 1.2
    @fact eval1([-1.9], p)                => 1.9

    @fact evalall([0.0], p)               => [0.0, 0.0]
    @fact evalall([1.2], p)               => [1.2, 1.2^2]
    @fact evalall([-1.9], p)              => [1.9, 1.9^2]
  end

  context("Fixed dim problem from an anydim one") do
    ap = anydim_problem("sumabs", fsabs, (0.0, 1.0), 0.0)
    p1 = as_fixed_dim_problem(ap, 1)

    @fact is_fixed_dimensional(p1)         => true
    @fact is_any_dimensional(p1)           => false
    @fact is_single_objective_problem(p1)  => true
    @fact is_multi_objective_problem(p1)   => false
    @fact numdims(p1)                      => 1

    @fact eval1([0.0], p1)                 => 0.0
    @fact eval1([1.2], p1)                 => 1.2
    @fact eval1([-1.9], p1)                => 1.9

    p3 = as_fixed_dim_problem(ap, 3)

    @fact is_fixed_dimensional(p3)         => true
    @fact is_any_dimensional(p3)           => false
    @fact is_single_objective_problem(p3)  => true
    @fact is_multi_objective_problem(p3)   => false
    @fact numdims(p3)                      => 3

    @fact eval1([0.0, 1.0, 2.0], p3)       => 3.0
    @fact eval1([-1.0, 1.0, 2.0], p3)      => 4.0
  end
end

facts("ShiftedAndBiasedProblem") do
  context("Only x shifted 1-dim") do
    ss = symmetric_search_space(1)
    subp = BlackBoxOptim.FixedDimProblem("sumabs", [fsabs], ss, [0.0])
    sp = BlackBoxOptim.shifted(subp)

    @fact is_fixed_dimensional(sp)         => true
    @fact is_any_dimensional(sp)           => false
    @fact is_single_objective_problem(sp)  => true
    @fact is_multi_objective_problem(sp)   => false
    @fact numdims(sp)                      => 1

    xs = sp.xshift
    @fact eval1(xs + [0.0], sp)                 => 0.0
    @fact eval1(xs + [1.2], sp)                 => 1.2
    @fact eval1(xs + [-1.9], sp)                => 1.9
  end

  context("Shifted and biased 2-dim") do
    ss = symmetric_search_space(2)
    subp = BlackBoxOptim.FixedDimProblem("sumabs", [fsabs], ss, [0.0])
    sp = BlackBoxOptim.shifted(subp; funcshift = 1.3)

    @fact is_fixed_dimensional(sp)         => true
    @fact is_any_dimensional(sp)           => false
    @fact is_single_objective_problem(sp)  => true
    @fact is_multi_objective_problem(sp)   => false
    @fact numdims(sp)                      => 2

    xs = sp.xshift
    @fact eval1(xs + [0.0, 1.0], sp)       => 1.0 + 1.3
    @fact eval1(xs + [1.2, -1.3], sp)      => 2.5 + 1.3
    @fact eval1(xs + [-1.9, 0.0], sp)      => 1.9 + 1.3
  end

  context("Within ftol") do
    ss = symmetric_search_space(2)
    subp = BlackBoxOptim.FixedDimProblem("sumabs", [fsabs], ss, [0.0])

    @fact fitness_is_within_ftol(subp, 0.1, 0.2)    => false
    @fact fitness_is_within_ftol(subp, 0.1, 0.09)   => true
  end
end


text:
_append(x::Tuple, y) = (x..., y)
_append(x::Tuple, y::Tuple) = (x..., y...)
flatten(x::Tuple) = _flatten((), x)
_flatten(result::Tuple, x::Tuple) = _flatten(_append(result, first(x)), Base.tail(x))
_flatten(result::Tuple, x::Tuple{}) = result

function rfftfreq(n::Integer,s::Real)::Vector{Float64}
    Float64[(n/2 - i)/s for i = n/2:-1:0]
end

@inline fsqrt(x) = @fastmath(sqrt(x))

function fftfreq(n::Integer,s::Real)::Vector{Float64}
    if iseven(n)
        return vcat(Float64[(n/2 - i)/s for i = n/2:-1:1],Float64[-i/s for i = n/2:-1:1])
    else return vcat(Float64[(n/2 - i)/s for i = n/2:-1:0],Float64[-i/s for i = (n-1)/2:-1:1])
    end
end

function mycopy!(out::Array{<:Real,N},inp::Array{<:Real,N}) where N
    @mthreads for j in TRANGE
        l = REAL_RANGES[j]
        for i in l
            @inbounds out[i] = inp[i]
        end
    end
end

mycopy!(o::ScalarField,i::ScalarField) = mycopy!(o.field.data,i.field.data)

mycopy!(o::VectorField,i::VectorField) = (mycopy!(o.rr.x,i.rr.x);
                                          mycopy!(o.rr.y,i.rr.y);
                                          mycopy!(o.rr.z,i.rr.z))

mycopy!(o::SymTrTenField,i::SymTrTenField) = (mycopy!(o.rr.xx,i.rr.xx);
                                              mycopy!(o.rr.xy,i.rr.xy);
                                              mycopy!(o.rr.xz,i.rr.xz);
                                              mycopy!(o.rr.yy,i.rr.yy);
                                              mycopy!(o.rr.yz,i.rr.yz))

@par function myscale!(field::AbstractArray{<:Real,N}) where N
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
        x = 1/(NRX*NY*NZ)
            @msimd for i in XRANGE2X
                @inbounds field[i,j,k] = x*field[i,j,k]
            end
    end
end

@par function fullmyscale!(field::AbstractArray{<:Real,N}) where N
    @mthreads for l in TRANGE
        x = 1/(NRX*NY*NZ)
        @msimd for i in REAL_RANGES[l]
            @inbounds field[i] = x*field[i]
        end
    end
end

myscale!(o::ScalarField) = myscale!(o.field.data)

myscale!(o::VectorField) = (myscale!(o.rr.x);
                            myscale!(o.rr.y);
                            myscale!(o.rr.z))

myscale!(o::SymTrTenField) = (myscale!(o.rr.xx);
                              myscale!(o.rr.xy);
                              myscale!(o.rr.xz);
                              myscale!(o.rr.yy);
                              myscale!(o.rr.yz))

myscale!(o::SymTenField) = (myscale!(o.rr.xx);
                            myscale!(o.rr.xy);
                            myscale!(o.rr.xz);
                            myscale!(o.rr.yy);
                            myscale!(o.rr.yz);
                            myscale!(o.rr.zz))

fullmyscale!(o::ScalarField) = fullmyscale!(o.field.data)

fullmyscale!(o::VectorField) = (fullmyscale!(o.rr.x);
                            fullmyscale!(o.rr.y);
                            fullmyscale!(o.rr.z))

fullmyscale!(o::SymTrTenField) = (fullmyscale!(o.rr.xx);
                              fullmyscale!(o.rr.xy);
                              fullmyscale!(o.rr.xz);
                              fullmyscale!(o.rr.yy);
                              fullmyscale!(o.rr.yz))

fullmyscale!(o::SymTenField) = (fullmyscale!(o.rr.xx);
                            fullmyscale!(o.rr.xy);
                            fullmyscale!(o.rr.xz);
                            fullmyscale!(o.rr.yy);
                            fullmyscale!(o.rr.yz);
                            fullmyscale!(o.rr.zz))

function myfourier!(field::A) where {T,N,N2,L,A<:Union{<:ScalarField{T,N,N2,L},<:VectorField{T,N,N2,L},<:SymTrTenField{T,N,N2,L},<:SymTenField{T,N,N2,L}}}
    rfft!(field)
    dealias!(field)
    myscale!(field)
    return nothing
end

function fullfourier!(field::A) where {T,N,N2,L,A<:Union{<:ScalarField{T,N,N2,L},<:VectorField{T,N,N2,L},<:SymTrTenField{T,N,N2,L},<:SymTenField{T,N,N2,L}}}
    rfft!(field)
    fullmyscale!(field)
    return nothing
end

function dealias!(f::AbstractArray{T,3}) where {T<:Complex}
    @mthreads for i in RANGEC
        @inbounds begin
            DEALIAS[i] && (f[i] = zero(T))
        end
   end
end

dealias!(o::VectorField) = (dealias!(o.c.x);
                            dealias!(o.c.y);
                            dealias!(o.c.z))

dealias!(o::SymTrTenField) = (dealias!(o.c.xx);
                              dealias!(o.c.xy);
                              dealias!(o.c.xz);
                              dealias!(o.c.yy);
                              dealias!(o.c.yz))

dealias!(o::SymTenField) = (dealias!(o.c.xx);
                              dealias!(o.c.xy);
                              dealias!(o.c.xz);
                              dealias!(o.c.yy);
                              dealias!(o.c.yz);
                              dealias!(o.c.zz))

function splitrange(lr,nt)
    a = UnitRange{Int}[]
    sizehint!(a,nt)
    n = lr÷nt
    r = lr%nt
    stop = 0
    init = 1
    for i=1:r
        stop=init+n
        push!(a,init:stop)
        init = stop+1
    end
    for i=1:(nt-r)
        stop=init+n-1
        push!(a,init:stop)
        init = stop+1
    end
    return (a...,)
end

function compute_shells2D(kx,ky,Nx,Ny)
    nShells2D = min(Nx,Ny÷2)
    maxdk2D = max(kx[2],ky[2])
    kh = zeros(nShells2D)
    numPtsInShell2D = zeros(Int,nShells2D)

    @inbounds for j=1:Ny
        for i=1:Nx
            K = fsqrt(kx[i]^2 + ky[j]^2)
            ii = round(Int,K/maxdk2D)+1
            ii > nShells2D && break
            kh[ii] += K
            numPtsInShell2D[ii] += 1
        end
    end

    @inbounds @simd for i in 1:length(kh)
        kh[i] = kh[i]/numPtsInShell2D[i]
    end

    return nShells2D, maxdk2D, numPtsInShell2D, kh
end

function calculate_Zf(d,kf,kh)
    Zf = zeros(size(kh))
    zft = haskey(d,:Zf) ? Symbol(d[:Zf]) : :cutoff
    if zft == :cutoff
        for (i,k) in enumerate(kh)
            Zf[i] = k <= kf
        end
    else
        for (i,k) in enumerate(kh)
            Zf[i] = tanh((kf-k) / (0.25*kf)) * (((kf-k) <= 0.0) ? 0.0 : 1.0)
        end
    end
    return Zf
end

function calculate_Zf(d,kf,kx,ky)
    Zf = zeros(length(kx),length(ky))
    zft = haskey(d,:Zf) ? Symbol(d[:Zf]) : :cutoff
    if zft == :tanh
        for j in eachindex(ky)
            for i in eachindex(kx)
                k = sqrt(kx[i]^2 + ky[j]^2)
                Zf[i,j] = tanh((kf-k) / (0.25*kf)) * (((kf-k) <= 0.0) ? 0.0 : 1.0)
            end
        end
    else
        for j in eachindex(ky)
            for i in eachindex(kx)
                k = sqrt(kx[i]^2 + ky[j]^2)
                Zf[i,j] = k <= kf
            end
        end
    end
    return Zf
end

@inline function Gaussfilter(Δ²::Real, k2::Real)
    aux = -Δ²/24
    return exp(aux*k2)
end

@inline Gaussfilter(Δ²::Real, i::Integer,j::Integer,k::Integer) = Gaussfilter(Δ²,K[i,j,k]⋅K[i,j,k])

@inline function boxfilter(Δ²::Real, k2::Real)
    aux = 0.5*fsqrt(Δ²*k2)
    return k2 == zero(k2) ? oneunit(k2) : sin(aux)/aux
end

@inline boxfilter(Δ::Real,i::Integer,j::Integer,k::Integer) = boxfilter(Δ,K[i,j,k]⋅K[i,j,k])

@inline function cutofffilter(Δ²::Real, k2::Real)
    aux = inv(Δ²)*π^2
    return k2 < aux
end

@inline cutofffilter(Δ²::Real,i::Integer,j::Integer,k::Integer) = boxfilter(Δ²,K[i,j,k]⋅K[i,j,k])

mymax(x::Float64,y::Float64) = max(abs(x),abs(y))
mymax(x::Float64,v::Vec{Float64}) = mymax(x,max(abs(v.x),abs(v.y),abs(v.z)))
mymax(v::Vec{Float64},x::Float64) = mymax(x,v)
mymax(v::Vec{Float64}) = mymax(v.x,mymax(v.y,v.z))
mymax(x::Float64) = abs(x)

function find_max(reduction,v)
    @mthreads for k in RZRANGE
        ti = Threads.threadid()
        umax = mymax(v[1,1,k])
        @inbounds for j in RYRANGE
            for i in RXRANGE
                umax = mymax(umax,mymax(v[i,j,k]))
            end
        end
        reduction[ti] = max(umax,reduction[ti])
    end
end

mysetfourier!(f) = (dealias!(f); setfourier!(f))


text:
#= ###################################################

* This pre-loads and saves the radial profiles of the bessel trasformations used in "likewhite.jl".
* These are used in `YXcov(x1, x2)`.
* This also saves the cls.
* The variables can be uploaded with
```
using JLD
@load "src/saved_cls/vKXnatxsq.jld"
@load "src/saved_cls/cls.jld"
```

!!!!!!!!!!
this needs updating
!!!!!!!!!

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

using PrimordialB
using JLD

cls  = PrimordialB.class(
    r           = 0.2,
    r0          = 1.0,
    lmax        = 20_000,
    omega_b     = 0.0224567,
    omega_cdm   = 0.118489,
    tau_reio    = 0.128312,
    theta_s     = 0.0104098,
    logA_s_1010 = 3.29056,
    n_s         = 0.968602
)

# --- compute the radial profile
vKE2atxsq, vKB2atxsq, ugrid = PrimordialB.vecKn_atxsq(2, cls)
vKE3atxsq, vKB3atxsq, ugrid = PrimordialB.vecKn_atxsq(3, cls)
vKE4atxsq, vKB4atxsq, ugrid = PrimordialB.vecKn_atxsq(4, cls)
vKE5atxsq, vKB5atxsq, ugrid = PrimordialB.vecKn_atxsq(5, cls)
vKE6atxsq, vKB6atxsq, ugrid = PrimordialB.vecKn_atxsq(6, cls)

# --- save
@save "src/saved_cls/vKXnatxsq.jld" ugrid vKE2atxsq vKB2atxsq vKE3atxsq vKB3atxsq vKE4atxsq vKB4atxsq vKE5atxsq vKB5atxsq vKE6atxsq vKB6atxsq

# --- also save the cls
@save "src/saved_cls/cls.jld" cls


text:
function steadystate_iterative!(rho0::Trho, H::AbstractMatrix, rates::DecayRates, J::Vector{TJ}, Jdagger::TJd, method!::Function, args...;kwargs...) where {Trho<:AbstractMatrix,TJ<:AbstractMatrix,TJd<:Union{Vector{TJ},Nothing}}
    # Size of the Hilbert space
    M = size(H,1)
    # Non-Hermitian Hamiltonian
    iHnh = nh_hamiltonian(H,rates,J,Jdagger)
    Jrho_cache = similar(rho0)

    # In-place update of y = Lx where L and x are respectively the vectorized
    # Liouvillian and the vectorized density matrix. y[1] is set to the trace
    # of the density matrix so as to enforce a trace one non-trivial solution.
    mvecmul!(y::AbstractVector, x::AbstractVector) = residual!(y,x,iHnh,rates,J,Jdagger,Jrho_cache)

    # Solution x must satisfy L.x = y with y[1] = tr(x) = 1 and y[j≠1] = 0.
    x0 = similar(rho0, M^2+1)
    x0 .= zero(eltype(rho0))
    x0[1] = zero(x0[1])
    x0[2:end] .= reshape(rho0, M^2)

    y = similar(rho0, M^2+1)
    y .= zero(y)
    y[1] = one(y[1])

    # Define the linear map lm: rho ↦ L(rho)
    lm = LinearMap{eltype(iHnh)}(mvecmul!, length(y)::Int, length(y)::Int; ismutating=true, issymmetric=false, ishermitian=false, isposdef=false)

    log = haskey(kwargs, :log) ? kwargs[:log] : false

    # Perform the stabilized biconjugate gradient procedure and devectorize rho
    if !log
        rho0 .= @views reshape(method!(x0,lm,y,args...;kwargs...)[2:end],(M,M))
        return rho0
    else
        R, history = method!(x0,lm,y,args...;kwargs...)
        rho0 .= @views reshape(R[2:end],(M,M))
        return rho0, history
    end
end

function nh_hamiltonian(H::AbstractMatrix, rates::Nothing, J::Vector{T}, Jdagger::Nothing) where {T<:AbstractMatrix}
    iHnh = -im*H
    for i=1:length(J)
        iHnh .+= -0.5 .* J[i]'J[i]
    end
    return iHnh
end
function nh_hamiltonian(H::AbstractMatrix, rates::Vector{Q}, J::Vector{T}, Jdagger::Nothing) where {Q<:Number,T<:AbstractMatrix}
    iHnh = -im*H
    for i=1:length(J)
        iHnh .+= -0.5rates[i] .* J[i]'J[i]
    end
    return iHnh
end
function nh_hamiltonian(H::AbstractMatrix, rates::Matrix{Q}, J::Vector{T}, Jdagger::Nothing) where {Q<:Number,T<:AbstractMatrix}
    iHnh = -im*H
    for i=1:length(J), j=1:length(J)
        iHnh -= 0.5rates[i,j] .* J[i]'J[j]
    end
    return iHnh
end
function nh_hamiltonian(H::AbstractMatrix, rates::Nothing, J::Vector{T}, Jdagger::Vector{T}) where {T<:AbstractMatrix}
    iHnh = -im*H
    for i=1:length(J)
        iHnh .+= -0.5 .* Jdagger[i]*J[i]
    end
    return iHnh
end
function nh_hamiltonian(H::AbstractMatrix, rates::Vector{Q}, J::Vector{T}, Jdagger::Vector{T}) where {Q<:Number,T<:AbstractMatrix}
    iHnh = -im*H
    for i=1:length(J)
        iHnh .+= -0.5rates[i] .* Jdagger[i]*J[i]
    end
    return iHnh
end
function nh_hamiltonian(H::AbstractMatrix, rates::Matrix{Q}, J::Vector{T}, Jdagger::Vector{T}) where {Q<:Number,T<:AbstractMatrix}
    iHnh = -im*H
    for i=1:length(J), j=1:length(J)
        iHnh -= 0.5rates[i,j] .* Jdagger[i]*J[j]
    end
    return iHnh
end

# With Jdagger
function residual!(y::AbstractVector, x::AbstractVector, iHnh::Tm, rates::Nothing, J::Vector{Tm}, Jdagger::Vector{Tm}, Jrho_cache::Tm) where {T<:T_blas,Tm<:DenseMatrix{T}}
    M = size(iHnh,1)
    y  .= zero(eltype(y));
    drho= @views reshape(y[2:end], M, M)
    rho = @views reshape(x[2:end], M, M)

    BLAS.gemm!('N', 'N', one(eltype(y)), iHnh, rho, one(eltype(y)), drho)
    BLAS.gemm!('N', 'C', one(eltype(y)), rho, iHnh, one(eltype(y)), drho)
    for i=1:length(J)
        BLAS.gemm!('N','N', one(eltype(y)), J[i], rho, zero(eltype(y)), Jrho_cache)
        BLAS.gemm!('N','N', one(eltype(y)), Jrho_cache, Jdagger[i], one(eltype(y)), drho)
    end

    y[1] = tr(rho)

    return y
end

function residual!(y::AbstractVector, x::AbstractVector, iHnh::Tm, rates::Vector{Q}, J::Vector{Tm}, Jdagger::Vector{Tm}, Jrho_cache::Tm) where {Q<:Number,T<:T_blas,Tm<:DenseMatrix{T}}
    M = size(iHnh,1)
    y  .= zero(eltype(y));
    drho= @views reshape(y[2:end], M, M)
    rho = @views reshape(x[2:end], M, M)

    BLAS.gemm!('N', 'N', one(eltype(y)), iHnh, rho, one(eltype(y)), drho)
    BLAS.gemm!('N', 'C', one(eltype(y)), rho, iHnh, one(eltype(y)), drho)
    for i=1:length(J)
        BLAS.gemm!('N','N', eltype(y)(rates[i]), J[i], rho, zero(eltype(y)), Jrho_cache)
        BLAS.gemm!('N','N', one(eltype(y)), Jrho_cache, Jdagger[i], one(eltype(y)), drho)
    end

    y[1] = tr(rho)

    return y
end

function residual!(y::AbstractVector, x::AbstractVector, iHnh::Tm, rates::Matrix{Q}, J::Vector{Tm}, Jdagger::Vector{Tm}, Jrho_cache::Tm) where {Q<:Number,T<:T_blas,Tm<:DenseMatrix{T}}
    M = size(iHnh,1)
    y  .= zero(eltype(y));
    drho= @views reshape(y[2:end], M, M)
    rho = @views reshape(x[2:end], M, M)

    BLAS.gemm!('N', 'N', one(eltype(y)), iHnh, rho, one(eltype(y)), drho)
    BLAS.gemm!('N', 'C', one(eltype(y)), rho, iHnh, one(eltype(y)), drho)
    for i=1:length(J), j=1:length(J)
        BLAS.gemm!('N','N', eltype(y)(rates[i,j]), J[i], rho, zero(eltype(y)), Jrho_cache)
        BLAS.gemm!('N','N', one(eltype(y)), Jrho_cache, Jdagger[j], one(eltype(y)), drho)

        BLAS.gemm!('N','N', -0.5one(eltype(y)), Jdagger[j], Jrho_cache, one(eltype(y)), drho)

        BLAS.gemm!('N','N', eltype(y)(rates[i,j]), rho, Jdagger[j], zero(eltype(y)), Jrho_cache)
        BLAS.gemm!('N','N', -0.5one(eltype(y)), Jrho_cache, J[i], one(eltype(y)), drho)
    end

    y[1] = tr(rho)

    return y
end

# Without
function residual!(y::AbstractVector, x::AbstractVector, iHnh::Tm, rates::Nothing, J::Vector{Tm}, Jdagger::Nothing, Jrho_cache::Tm) where {T<:T_blas,Tm<:DenseMatrix{T}}
    M = size(iHnh,1)
    y  .= zero(eltype(y));
    drho= @views reshape(y[2:end], M, M)
    rho = @views reshape(x[2:end], M, M)

    BLAS.gemm!('N', 'N', one(eltype(y)), iHnh, rho, one(eltype(y)), drho)
    BLAS.gemm!('N', 'C', one(eltype(y)), rho, iHnh, one(eltype(y)), drho)
    for Ji=J
        BLAS.gemm!('N','N', one(eltype(y)), Ji, rho, zero(eltype(y)), Jrho_cache)
        BLAS.gemm!('N','C', one(eltype(y)), Jrho_cache, Ji, one(eltype(y)), drho)
    end

    y[1] = tr(rho)

    return y
end

function residual!(y::AbstractVector, x::AbstractVector, iHnh::Tm, rates::Vector{Q}, J::Vector{Tm}, Jdagger::Nothing, Jrho_cache::Tm) where {Q<:Number,T<:T_blas,Tm<:DenseMatrix{T}}
    M = size(iHnh,1)
    y  .= zero(eltype(y));
    drho= @views reshape(y[2:end], M, M)
    rho = @views reshape(x[2:end], M, M)

    BLAS.gemm!('N', 'N', one(eltype(y)), iHnh, rho, one(eltype(y)), drho)
    BLAS.gemm!('N', 'C', one(eltype(y)), rho, iHnh, one(eltype(y)), drho)
    for i=1:length(J)
        BLAS.gemm!('N','N', eltype(y)(rates[i]), J[i], rho, zero(eltype(y)), Jrho_cache)
        BLAS.gemm!('N','C', one(eltype(y)), Jrho_cache, J[i], one(eltype(y)), drho)
    end

    y[1] = tr(rho)

    return y
end

function residual!(y::AbstractVector, x::AbstractVector, iHnh::Tm, rates::Matrix{Q}, J::Vector{Tm}, Jdagger::Nothing, Jrho_cache::Tm) where {Q<:Number,T<:T_blas,Tm<:DenseMatrix{T}}
    M = size(iHnh,1)
    y  .= zero(eltype(y));
    drho= @views reshape(y[2:end], M, M)
    rho = @views reshape(x[2:end], M, M)

    BLAS.gemm!('N', 'N', one(eltype(y)), iHnh, rho, one(eltype(y)), drho)
    BLAS.gemm!('N', 'C', one(eltype(y)), rho, iHnh, one(eltype(y)), drho)
    for i=1:length(J), j=1:length(J)
        BLAS.gemm!('N','N', eltype(y)(rates[i,j]), J[i], rho, zero(eltype(y)), Jrho_cache)
        BLAS.gemm!('N','C', one(eltype(y)), Jrho_cache, J[j], one(eltype(y)), drho)

        BLAS.gemm!('C','N', -0.5one(eltype(y)), J[j], Jrho_cache, one(eltype(y)), drho)

        BLAS.gemm!('N','C', eltype(y)(rates[i,j]), rho, J[j], zero(eltype(y)), Jrho_cache)
        BLAS.gemm!('N','N', -0.5one(eltype(y)), Jrho_cache, J[i], one(eltype(y)), drho)
    end

    y[1] = tr(rho)

    return y
end

# No-BLAS fallback

# With Jdagger
function residual!(y::AbstractVector, x::AbstractVector, iHnh::AbstractMatrix, rates::Nothing, J::Vector{TJ}, Jdagger::Vector{TJ}, Jrho_cache::Trho) where {Trho<:AbstractMatrix,TJ<:AbstractMatrix}
    M = size(iHnh,1)
    y  .= zero(eltype(y));
    drho  = Trho(@views reshape(y[2:end], M, M))
    rho = Trho(@views reshape(x[2:end], M, M))

    drho .= iHnh * rho .+ rho *iHnh'
    for i=1:length(J)
        Jrho_cache .= J[i] * rho
        drho .+= Jrho_cache * Jdagger[i]
    end
    @views y[2:end] .= reshape(drho, M^2)

    y[1] = tr(rho)

    return y
end

function residual!(y::AbstractVector, x::AbstractVector, iHnh::AbstractMatrix, rates::Vector{Q}, J::Vector{TJ}, Jdagger::Vector{TJ}, Jrho_cache::Trho) where {Q<:Number,Trho<:AbstractMatrix,TJ<:AbstractMatrix}
    M = size(iHnh,1)
    y  .= zero(eltype(y));
    drho= @views reshape(y[2:end], M, M)
    rho = @views reshape(x[2:end], M, M)

    drho .= iHnh * rho .+ rho *iHnh'
    for i=1:length(J)
        Jrho_cache .= rates[i] .* J[i] * rho
        drho .+= Jrho_cache * Jdagger[i]
    end

    y[1] = tr(rho)

    return y
end

function residual!(y::AbstractVector, x::AbstractVector, iHnh::AbstractMatrix, rates::Matrix{Q}, J::Vector{TJ}, Jdagger::Vector{TJ}, Jrho_cache::Trho) where {Q<:Number,Trho<:AbstractMatrix,TJ<:AbstractMatrix}
    M = size(iHnh,1)
    y  .= zero(eltype(y));
    drho= @views reshape(y[2:end], M, M)
    rho = @views reshape(x[2:end], M, M)

    drho .= iHnh * rho .+ rho *iHnh'
    for i=1:length(J), j=1:length(J)
        Jrho_cache .= rates[i,j] .* J[i] * rho
        drho .+= Jrho_cache * Jdagger[j]

        drho .+= -0.5 .* Jdagger[j] * Jrho_cache

        Jrho_cache .= rates[i,j] .* rho * Jdagger[j]
        drho .+= -0.5 .* Jrho_cache * J[i]
    end

    y[1] = tr(rho)

    return y
end

# Without
function residual!(y::AbstractVector, x::AbstractVector, iHnh::AbstractMatrix, rates::Nothing, J::Vector{TJ}, Jdagger::Nothing, Jrho_cache::Trho) where {Trho<:AbstractMatrix,TJ<:AbstractMatrix}
    M = size(iHnh,1)
    y  .= zero(eltype(y));
    drho  = Trho(@views reshape(y[2:end], M, M))
    rho = Trho(@views reshape(x[2:end], M, M))

    drho .= iHnh * rho .+ rho *iHnh'
    for Ji=J
        Jrho_cache .= Ji * rho
        drho .+= Jrho_cache * Ji'
    end
    @views y[2:end] .= reshape(drho, M^2)

    y[1] = tr(rho)

    return y
end

function residual!(y::AbstractVector, x::AbstractVector, iHnh::AbstractMatrix, rates::Vector{Q}, J::Vector{TJ}, Jdagger::Nothing, Jrho_cache::Trho) where {Q<:Number,Trho<:AbstractMatrix,TJ<:AbstractMatrix}
    M = size(iHnh,1)
    y  .= zero(eltype(y));
    drho= @views reshape(y[2:end], M, M)
    rho = @views reshape(x[2:end], M, M)

    drho .= iHnh * rho .+ rho *iHnh'
    for i=1:length(J)
        Jrho_cache .= rates[i] .* J[i] * rho
        drho .+= Jrho_cache * J[i]'
    end
    @views y[2:end] .= reshape(drho, M^2)

    y[1] = tr(rho)

    return y
end

function residual!(y::AbstractVector, x::AbstractVector, iHnh::AbstractMatrix, rates::Matrix{Q}, J::Vector{TJ}, Jdagger::Nothing, Jrho_cache::Trho) where {Q<:Number,Trho<:AbstractMatrix,TJ<:AbstractMatrix}
    M = size(iHnh,1)
    y  .= zero(eltype(y));
    drho= @views reshape(y[2:end], M, M)
    rho = @views reshape(x[2:end], M, M)

    drho .= iHnh * rho .+ rho *iHnh'
    for i=1:length(J), j=1:length(J)
        Jrho_cache .= rates[i,j] .* J[i] * rho
        drho .+= Jrho_cache * J[j]'

        drho .+= -0.5 .* J[j]' * Jrho_cache

        Jrho_cache .= rates[i,j] .* rho * J[j]'
        drho .+= -0.5 .* Jrho_cache * J[i]
    end

    y[1] = tr(rho)

    return y
end


text:
using Gadfly
using ANOVADDPTest
using DataFrames
using Distributions
using Random

function generate_plot_poisson(N, G, distinct_groups)
    rng = MersenneTwister(1)
    m = PoissonDDP(rng, N, G)
    x = rand(rng, 1:G, N)
    y = rand(rng, Poisson(2), N)
    for i in 1:N
        x[i] in distinct_groups && (y[i] = rand(rng, Poisson(4)))
    end
    data = PoissonData(x, y)
    pred = PoissonData(expandgrid(1:G, 0:8)...)
    ch = train(rng, m, data, pred)
    df = DataFrame(x = pred.x, y = pred.y, f = mean(ch.f))
    plot(df, x = :y, y = :f, color = :x, Geom.line, Scale.color_discrete())
end

generate_plot_poisson(1000, 3, [])
generate_plot_poisson(1000, 3, [2])
generate_plot_poisson(1000, 3, [3])
generate_plot_poisson(1000, 3, [2, 3])


text:
# Analyze convergence of Monte Carlo tree search for various horizon lenghts

using SubmodularMaximization
using PyPlot
using Statistics
using Base.Threads
using Base.Iterators
using JLD2
using Printf

close("all")

println("Loading rollout library")
library_file = "./data/rollout_library.jld2"
@time @load library_file data num_robots trials
rollouts = data

# Create a map from number of robots to a vector target filters
filters = Dict(map(num_robots) do num_robots
                 num_robots => vcat(map(trials) do trial
                                      rollout = rollouts[num_robots, trial]

                                      vcat([x.target_filters for x in
                                            rollout.samples]...)
                                    end...)
               end)

# Concatenate and sort filter data (increasing) by robot
probabilities = Dict(map(num_robots) do num_robots
                       my_filters = filters[num_robots]
                       num_robots => sort(vcat(map(my_filters) do filter
                                            @views filter.data[:]
                                          end...))
                     end)

#
# Plot cumulative distribution functions
#

# The cdf value at p is the fraction of the sum of probabilities less than p
xs = 10 .^ (-1 .* reverse(0:.1:10))
map(num_robots) do num_robots
  let probabilities = probabilities[num_robots]

    cdf_vals = map(xs) do p
      index = findfirst(x->x > p, probabilities)

      # If there is no value greater than p, then all is less than p
      total_value = sum(probabilities)
      if isnothing(index)
        1.0
      else
        sum(probabilities[1:index-1]) / total_value
      end
    end

    label = string(num_robots, " robots")
    semilogx(xs, cdf_vals, label=label)
  end
end

ylabel("Cumulative value frac. less than")
xlabel("Filter probability")
legend()
grid()

save_fig("fig", "cell_value_cdf")

#
# Plot the cell value fraction
#

figure()

# We are interested in the fraction of cell values that are less than a given
# value
xs = 10 .^ (-1 .* reverse(0:.1:10))
map(num_robots) do num_robots
  let probabilities = probabilities[num_robots]

    cdf_vals = map(xs) do p
      index = findfirst(x->x > p, probabilities)

      # If there is no value greater than p, then all are less than p
      if isnothing(index)
        1.0
      else
        (index - 1) / length(probabilities)
      end
    end

    label = string(num_robots, " robots")
    semilogx(xs, cdf_vals, label=label)
  end
end

ylabel("Cumulative frac. cells less than")
xlabel("Filter probability")
legend()
grid()

save_fig("fig", "cell_value_fraction")

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 5c3096b0-01b6-11eb-2c01-5ba0aac27b70
using CardioModels, PlutoUI, Plots, Distributions

# ╔═╡ 0c720f5f-a8be-4d8e-a803-461118c3b2be
md"
# DeBoer Model

The following equations manifest the model:

#### Effective pressure:
The systolic pressure is transformed to be in accordance with the sigmoid-shaped activation curve of the baroreflex. The constant $S_0$ defines the working point.

$S_{n}^{\star} = f(S_n)\hspace{5mm}\text{with} \hspace{5mm} f(S_n) = S_0 + 18\text{atan}\left(\frac{S_{n} -S_0 }{18}\right)$

#### Baroreflex on Heart Rate
The transfer function $G_a(z)$ relates the effective systolic pressure to the RR interval. It models fast vagally mediated and slower, sympathetic, baroreflex effects.

$I_{n} = G_a (z) S_{n}^{\star}$

#### Baroreflex on Peripheral Resistance
The peripheral resistance is influenced by the baroreflex mediated through the sympathetic nervous system reflected in the transfer function $G_b(z)$. $T_0$ is an empirical offset constant (around 3000 ms).

$T_{ n } = R_n C = T_{ 0 } + G_b(z)S_{n}^{\star}$

#### Properties of Myocardium & Respiration
Starling's law governs the pulse pressure through the coupling constant $\gamma$. Respiration also enters through the pulse pressure, modeled as a sine.

$P_{ n } = \gamma I_{n}  + A_{\rho}\sin(2\pi f_{\rho} t) \hspace{5mm} \text{with} \hspace{5mm} t = \sum I_n$

#### Windkessel
Diastolic pressure is derived from the two-element Windkessel model.

$D_{ n } = S_{n} \exp \left( - \frac{I_{n}}{ T_{n} } \right)$

#### Time Propagation
The next beat is calculated through the standard relation:

$S_{n+1} = P_{n} + D_{n}$
"

# ╔═╡ 2e8472f6-01b4-11eb-1f49-b35e9ee7c70b
model = DeBoerModel();

# ╔═╡ 373eadc4-01b6-11eb-0ef3-73d0e4518a24
md"""
N:
$(@bind N NumberField(1:1000, default = 100))
"""

# ╔═╡ 232d91d8-01b6-11eb-08d6-3343be95b77f
S, D, I, T, ρ = predict!(model, N, burnIn = 100);

# ╔═╡ b19bf544-b153-41a0-a013-992da6cdd0cd
md"""
Show Spectrum:
$(@bind psd CheckBox())
"""

# ╔═╡ 579f32aa-01b6-11eb-3a79-8f2207aeb03c
begin
	if psd
		psdplot(S)
	else
		s = plot(S, lab = "", ylabel = "mmHg", color = :red)
		plot!(D, lab = "", title = "BP", ylabel = "mmHg", color = :green)
		t = plot(T, lab = "", title = "T", ylabel = "a.u.", color = :orange)
		i = plot(I, lab = "", title = "I", ylabel = "ms", color = :blue)
		r = plot(ρ, lab = "", title = "ρ", ylabel = "a.u.", xlab = "beats", color = :black)
		plot(s,t,i,r, layout = (4,1), dpi = 500)
	end
end

# ╔═╡ ac3a7eac-6bba-11eb-37ec-49854b22053f
md"
# Simulating a phenylephrine injection

"

# ╔═╡ 53b158a2-6bbb-11eb-0fa7-87390c33cb64
begin
	model_phe = DeBoerModel(a0 = 9, A = 1,  noiseI = Normal(0, 25), noiseS = Normal(0, 2))
	# pre-injection
	Npre = 50
	S1,D1,I1,T1,R1 = predict!(model_phe, Npre, burnIn = 100)

	# injection
	duration = 15
	S2,D2,I2,T2,R2 = phenylephrine(model_phe, 1000, duration)

	# post-injection
	Npost = 10 # length of plateau
	S3,D3,I3,T3,R3 = predict!(model_phe, Npost, burnIn = 0)

	# back to default
	S4,D4,I4,T4,R4 = phenylephrine(model_phe, -1000, 10)

	plot([S1; S2; S3; S4], lab = "", title = "BP", ylab = "mmHg", color = :red)
	pBP = plot!([D1; D2; D3; D4], lab = "", color = :green)
	pIBI = plot([I1; I2; I3; I4], lab = "", ylab = "ms", title = "I", color = :blue)
	plot(pBP, pIBI, layout = (2,1))
	vline!([Npre Npre], color = :black, lab = "", dpi = 500)

end

# ╔═╡ b312bc92-6bbf-11eb-131e-c946cf3eb23f
begin
	xbp = 0:30
	yms = 9 .* xbp
	scatter(S2 .- minimum(S2), I2 .- minimum(I2), lab = "", title = "Simulated response")
	plot!(xbp, yms, color = :black, lab = "BRS 9 ms/mmHg", ylab = "\\Delta I [ms]", xlab = "\\Delta S [mmHg]", dpi = 500)

end

# ╔═╡ Cell order:
# ╠═5c3096b0-01b6-11eb-2c01-5ba0aac27b70
# ╟─0c720f5f-a8be-4d8e-a803-461118c3b2be
# ╠═2e8472f6-01b4-11eb-1f49-b35e9ee7c70b
# ╟─373eadc4-01b6-11eb-0ef3-73d0e4518a24
# ╠═232d91d8-01b6-11eb-08d6-3343be95b77f
# ╟─b19bf544-b153-41a0-a013-992da6cdd0cd
# ╟─579f32aa-01b6-11eb-3a79-8f2207aeb03c
# ╟─ac3a7eac-6bba-11eb-37ec-49854b22053f
# ╠═53b158a2-6bbb-11eb-0fa7-87390c33cb64
# ╠═b312bc92-6bbf-11eb-131e-c946cf3eb23f


text:
using Printf
using Statistics
using LinearAlgebra
using MDBM
using PyPlot
pygui(true);
using Interpolations
using Plots

dim=2 #DoF of the system (after Cauchy transcription)
mult=3 #Ns=dim*mult
n=40 #timestep number
#choosing numerical simulation method
#RK4: 4th order Runge-Kutta
method="RK4"
i=20 #number of iteration
iT=8

#parameters of the delayed Mathieu equation
tau=2*pi
omega=1
kappa=0.1
delta=1.5
epsilon=1
b=0.25
T=2*pi/omega

#timestep
dt=tau/(n-1)

#matrix coefficients
function A(t)
 [-kappa -(delta+epsilon*cos(omega*t)); 1 0]
end
B=[0 b; 0 0]

#checking the relation of tau and T
function checkT(tau,T)
    if tau <= T
        println("OK")
    else
        println("Error")
    end
end

checkT(tau,T)
k=Int16(floor(T/tau)) #k factor

#initial function ()
vect1=rand(n,mult*dim)*(-10)+ones(n,mult*dim)*5

#interpolate
function it(A)
    matrdim=size(A,2)-1
    scaleitp=real(A[:,1])
    knots=(scaleitp,)
    ARe=real(A)
    AIm=imag(A)

    imRe=[interpolate(knots, ARe[:,2],Gridded(Linear()))]
    for j = 3:matrdim+1
        imRetemp=[interpolate(knots, real(ARe[:,j]), Gridded(Linear()))]
        imRe=vcat(imRe,imRetemp)
    end
    imIm=[interpolate(knots, AIm[:,2],Gridded(Linear()))]
    for j = 3:matrdim+1
        imImtemp=[interpolate(knots, real(AIm[:,j]), Gridded(Linear()))]
        imIm=vcat(imIm,imImtemp)
    end
    return(hcat(imRe,imIm))
end

#substitution in interpolating function
function sub(it,t)
    subdim=size(it,1)
    out=zeros(ComplexF64,subdim)
    for j = 1:subdim
        out[j]=it[j,1](t)+it[j,2](t)*im
    end
    return(out)
end
collect(-tau:dt:iT*tau)
function simulation(init)
        solreturn=zeros(ComplexF64,1,1+mult*dim) #matrix for time integration solution (test)
        eval=zeros(ComplexF64,mult*dim,1) #eigenvalue matrix

        sol00=init #S0 inital random array
        #iteration loop
        for g = 1:i
        sol=zeros(iT*n-iT,mult*dim) #empty solution vector
        tvect0=collect(-tau:dt:iT*tau+dt/2) #discretized time vector
        sol0=hcat(tvect0,vcat(sol00,sol)) #constructing solution array
            #time integration
            for m = 2:dim:(dim*mult-(dim-2)) #choosing "one physical system"
            sol0m=hcat(sol0[:,1],sol0[:,m:m+dim-1])
                #time integration steps
                for j = 0:iT*n-iT-1
                    sol0mit=it(sol0m)

                    Y1=sub(sol0mit,j*dt)
                    Y2=Y1+(dt/2)*(A(j*dt)*Y1+B*sub(sol0mit,-tau+j*dt))
                    Y3=Y1+(dt/2)*(A((j+0.5)*dt)*Y2+B*sub(sol0mit,-tau+(j+0.5)*dt))
                    Y4=Y1+dt*(A((j+0.5)*dt)*Y3+B*sub(sol0mit,-tau+(j+0.5)*dt))
                    Y=Y1+(dt/6)*((A(j*dt)*Y1+B*sub(sol0mit,-tau+j*dt))+2*(A((j+0.5)*dt)*Y2+B*sub(sol0mit,-tau+(j+0.5)*dt))+2*(A((j+0.5)*dt)*Y3+B*sub(sol0mit,-tau+(j+0.5)*dt))+(A((j+1)*dt)*Y4+B*sub(sol0mit,-tau+(j+1)*dt)))

                    sol0m[n+j+1,2:2+(dim-1)]=transpose(Y) #filling solution array of the physical system
                end
                sol0[:,m:m+dim-1]=sol0m[:,2:2+(dim-1)] #filling iterational solution array
            end
            solreturn=vcat(solreturn,sol0) #filling global solution array (test)
            S=sol0[(iT-1)*n-((iT-1)-1):(iT)*n-(iT-1),2:mult*dim+1] #taking array Sj (n x Ns)
            V=sol0[iT*n-(iT-1):(iT+1)*n-iT,2:mult*dim+1] #taking array Vj (n x Ns)
            #H=pinv(S)*V #calculating matrix Hj (Ns x Ns)
            H=pinv(S)*V
            Am=eigvecs(H) #calculation eigenvectors of Hj, creating matrix Am (Ns x Ns)
            sol00=zeros(ComplexF64,n,mult*dim) #calculating matrix Sj+1 (n x Ns)
            for o = 1:n
                sol00[o,:]=Am*V[o,:]
            end
            sol00=sol00
            evm=eigvals(H) #obtaining dominant eigenvalues in the iteration step
            eval=hcat(eval,evm) #filling up the eigenvalue array
        end
        return(eval[:,2:i+1]) #returning the global eigenvalue array
        #return(solreturn)
end

aa=simulation(vect1)

aa[6000:6100,:]

function evplot(aaa,k)
U=aaa[:,k]
eigRe1=[real(U[n,1]) for n in 1:size(U,1)]
eigIm1=[imag(U[n,1]) for n in 1:size(U,1)]
plot1=Plots.plot(eigRe1,eigIm1,seriestype=:scatter)
Re_circle=[cos(n) for n in 0:0.05:2*pi]
Im_circle=[sin(n) for n in 0:0.05:2*pi]
return(Plots.plot!(plot1,Re_circle,Im_circle))
end

evplot(aa,1)
evplot(aa,2)
evplot(aa,3)
evplot(aa,4)
evplot(aa,5)
evplot(aa,6)
evplot(aa,7)
evplot(aa,8)
evplot(aa,9)
evplot(aa,10)
evplot(aa,12)
evplot(aa,11)
evplot(aa,13)
evplot(aa,14)
evplot(aa,15)
evplot(aa,16)
evplot(aa,17)
evplot(aa,18)
evplot(aa,19)
evplot(aa,20)
evplot(aa,21)
evplot(aa,22)
evplot(aa,23)
evplot(aa,24)
evplot(aa,25)
evplot(aa,26)
evplot(aa,27)
evplot(aa,28)
evplot(aa,29)
evplot(aa,30)
evplot(aa,31)
evplot(aa,32)

solt1=real(aa)
v1plot=[solt1[j,2] for j in 1:size(solt1,1)]
v2plot=[solt1[j,3] for j in 1:size(solt1,1)]
v3plot=[solt1[j,4] for j in 1:size(solt1,1)]
v4plot=[solt1[j,5] for j in 1:size(solt1,1)]
v5plot=[solt1[j,6] for j in 1:size(solt1,1)]
v6plot=[solt1[j,7] for j in 1:size(solt1,1)]
Plots.plot(v1plot)
Plots.plot(v2plot)
Plots.plot(v3plot)
Plots.plot(v4plot)
Plots.plot(v5plot)
Plots.plot(v6plot)

function avarr(A)
    dim1=size(A,1)
    dim2=size(A,2)
    sum=0
    for i=1:dim1
            for j=1:dim2
                sum=sum+norm(A[i,j])
            end
    end
    return(sum/(dim1*dim2))
end

avarr([1 2; 3 4])
(1/2.5)*[1 2; 3 4]

norm(-1)
avarr([1 2+4*im; 3 4])


text:
module SparseMatrix

using LinearAlgebra #: Adjoint
import SparseArrays #: SparseMatrixCSC, nnz

export AbstractSparseM
export SparseCSR, SparseCSC, SparseCOO, BlockCOO
export spdiag
export block

include("_sparse.jl")
include("_mul.jl")
include("_spdiag.jl")
include("_block.jl")

end