text:
# This file is a part of ParallelProcessingTools.jl, licensed under the MIT License (MIT).

function _current_thread_selected(threadsel::Union{Integer,AbstractVector{<:Integer}})
    tid = threadid()
    checkindex(Bool, tid:tid, threadsel)
end

@static if VERSION >= v"1.3.0-alpha.0"

# From Julia PR 32477:
function _run_on(t::Task, tid)
    @assert !istaskstarted(t)
    t.sticky = true
    ccall(:jl_set_task_tid, Cvoid, (Any, Cint), t, tid-1)
    schedule(t)
    return t
end

# Adapted from Julia PR 32477:
function _threading_run(func, threadsel::AbstractVector{<:Integer})
    tasks = Vector{Task}(undef, length(eachindex(threadsel)))
    for tid in threadsel
        i = firstindex(tasks) + (tid - first(threadsel))
        tasks[i] = _run_on(Task(func), tid)
    end
    foreach(wait, tasks)
    return nothing
end

_threading_run(func, threadsel::Integer) = _threading_run(func, threadsel:threadsel)

function _thread_exec_func(threadsel, expr)
    quote
        local thread_body_wrapper_fun
        let threadsel_eval = $(esc(threadsel))
            function thread_body_wrapper_fun()
                $(esc(expr))
            end
            if _current_thread_selected(threadsel_eval)
                thread_body_wrapper_fun()
            else
                _threading_run(thread_body_wrapper_fun, threadsel_eval)
            end
            nothing
        end
    end
end

else #VERSION < v"1.3.0-alpha.0"

const _thread_local_error_err = ThreadLocal{Any}(undef)
const _thread_local_error_set = ThreadLocal{Bool}(undef)

_clear_thread_local_errors() = fill!(getallvalues(_thread_local_error_set), false)

function _check_thread_local_errors()
    i = something(findfirst(isequal(true), getallvalues(_thread_local_error_set)), 0)
    (i > 0) && throw(getallvalues(_thread_local_error_err)[i])
    nothing
end

function _set_thread_local_error(err)
    _thread_local_error_err[] = err
    _thread_local_error_set[] = true
end

function _check_threadsel(threadsel::Union{Integer,AbstractVector{<:Integer}})
    if !checkindex(Bool, allthreads(), threadsel)
        throw(ArgumentError("Thread selection not within available threads"))
    end
    threadsel
end

function _run_on_threads(f)
    try
        @assert(!Base.Threads.in_threaded_loop[], "Can't nest threaded execution")
        _clear_thread_local_errors()
        Base.Threads.in_threaded_loop[] = true
        ccall(:jl_threading_run, Ref{Cvoid}, (Any,), f)
    finally
        Base.Threads.in_threaded_loop[] = false
        _check_thread_local_errors()
    end
end

function _thread_exec_func(threadsel, expr)
    quote
        local thread_body_wrapper_fun
        let threadsel_eval = $(esc(threadsel))
            function thread_body_wrapper_fun()
                try
                    if Base.Threads.threadid() in threadsel_eval
                        $(esc(expr))
                    end
                catch err
                    _set_thread_local_error(err)
                    rethrow()
                end
            end
            if _current_thread_selected(threadsel_eval)
                thread_body_wrapper_fun()
            else
                _run_on_threads(thread_body_wrapper_fun)
            end
            nothing
        end
    end
end

end # Julia version-dependent code

"""
    allthreads()

Convencience function, returns `1:Base.Threads.nthreads()`.
"""
allthreads() = 1:nthreads()
export allthreads

"""
    @onthreads threadsel expr

Execute code in `expr` in parallel on the threads in `threadsel`.

`threadsel` should be a single thread-ID or a range (or array) of thread-ids.
If `threadsel == Base.Threads.threadid()`, `expr` is run on the current
tread with only minimal overhead.

Note: Currently, multiple `@onthreads` sections will not run in parallel
to each other, even if they use disjunct sets of threads, due to limitations
of the Julia multithreading implementation. This restriction is likely to
disappear in future Julia versions.

In contrast to `Base.Threads.@threads`, `@onthreads` does forward
exceptions to the caller.

Example 1:

```juliaexpr
tlsum = ThreadLocal(0.0)
data = rand(100)
@onthreads allthreads() begin
    tlsum[] = sum(workpart(data, allthreads(), Base.Threads.threadid()))
end
sum(getallvalues(tlsum)) ≈ sum(data)
```

Example 2:

```julia
# Assuming 4 threads:
tl = ThreadLocal(42)
threadsel = 2:3
@onthreads threadsel begin
    tl[] = Base.Threads.threadid()
end
getallvalues(tl)[threadsel] == [2, 3]
getallvalues(tl)[[1,4]] == [42, 42]
```
"""
macro onthreads(threadsel, expr)
    _thread_exec_func(threadsel, expr)
end
export @onthreads

macro onallthreads(expr)
    quote
        Base.depwarn("`@onallthreads expr` is deprecated, use `@onthreads allthreads() expr` instead.", nothing)
        $(_thread_exec_func(:(ParallelProcessingTools.allthreads()), expr))
    end
end
export @onallthreads

function ThreadLocal{T}(f::Base.Callable) where {T}
    result = ThreadLocal{T}(undef)
    result.value
    @onthreads allthreads() result.value[threadid()] = f()
    result
end

"""
    @mt_async expr

Spawn a Julia task running `expr` asynchronously.

Compatible with `@sync`. Uses a multi-threaded task scheduler if available (on
Julia >= v1.3).

Equivalent to `Base.@async` on Julia <= v1.2, equivalent to
`Base.Threads.@spawn` on Julia >= v1.3.
"""
macro mt_async(expr)
    # Code taken from Base.@async and Base.Threads.@spawn:
    thunk = esc(:(()->($expr)))
    var = esc(Base.sync_varname)
    quote
        local task = Task($thunk)
        @static if VERSION >= v"1.3.0-alpha.0"
            task.sticky = false
        end
        if $(Expr(:isdefined, var))
            push!($var, task)
        end
        schedule(task)
        task
    end
end
export @mt_async

"""
    @mt_out_of_order begin expr... end

Runs all top-level expressions in `begin expr... end` on parallel tasks.
On Julia >= v1.3, the tasks will run multi-threaded.

Example:

```
@mt_out_of_order begin
    a = foo()
    bar()
    c = baz()
end

will run `a = foo()`, `bar()` and `c = baz()` in parallel and in arbitrary
order, results of assignments will appear in the outside scope.
"""
macro mt_out_of_order(ex)
    if !(ex isa Expr && ex.head == :block)
        throw(ErrorException("@mt_out_of_order expects a code block as it's argument"))
    end

    exprs = ex.args
    idxs = eachindex(ex.args)
    tasks = gensym(:tasks)
    handle_results = Vector{Expr}()
    for i in idxs
        if exprs[i] isa Expr && exprs[i].head == :(=)
            trg = exprs[i].args[1]
            val = exprs[i].args[2]
            if val isa Expr
                exprs[i] = :(push!($tasks, @mt_async($(esc(val)))))
                push!(handle_results, :($(esc(trg)) = fetch(popfirst!($tasks))))
            else
                exprs[i] = esc(exprs[i])
            end
        elseif exprs[i] isa Expr
            ftvar = gensym()
            exprs[i] = :(push!($tasks, @mt_async($(esc(exprs[i])))))
            push!(handle_results, :(wait(popfirst!($tasks))))
        else
            exprs[i] = esc(exprs[i])
        end
    end
    pushfirst!(exprs, :($tasks = Vector{Task}()))
    append!(exprs, handle_results)
    push!(exprs, :(@assert isempty($tasks)))
    push!(exprs, :nothing)
    ex
end
export @mt_out_of_order


text:
import HTTP, JSON2
using Test
using Dash

@testset "callback! single output" begin
    app = dash()
    app.layout = html_div() do
            dcc_input(id = "my-id", value="initial value", type = "text"),
            html_div(id = "my-div")
        end

    callback!(app, Output("my-div", "children"), Input("my-id", "value")) do value
        return value
    end
    @test length(app.callbacks) == 1
    @test haskey(app.callbacks, Symbol("my-div.children"))
    @test app.callbacks[Symbol("my-div.children")].func("test") == "test"

    handler = make_handler(app)
    request = HTTP.Request("GET", "/_dash-dependencies")
    resp = HTTP.handle(handler, request)
    deps = JSON2.read(String(resp.body))

    @test length(deps) == 1
    cb = deps[1]
    @test cb.output == "my-div.children"
    @test cb.inputs[1].id == "my-id"
    @test cb.inputs[1].property == "value"
    @test :clientside_function in keys(cb)
    @test isnothing(cb.clientside_function)

    handler = Dash.make_handler(app)
    test_json = """{"output":"my-div.children","changedPropIds":["my-id.value"],"inputs":[{"id":"my-id","property":"value","value":"test"}]}"""
    request = HTTP.Request("POST", "/_dash-update-component", [], Vector{UInt8}(test_json))
    response = HTTP.handle(handler, request)
    @test response.status == 200
    resp_obj = JSON2.read(String(response.body))
    @test in(:multi, keys(resp_obj))
    @test resp_obj.response.var"my-div".children == "test"

end

@testset "callback! multi output" begin
    app = dash()
    app.layout = html_div() do
            dcc_input(id = "my-id", value="initial value", type = "text"),
            html_div(id = "my-div"),
            html_div(id = "my-div2")
        end
    callback!(app, [Output("my-div","children"), Output("my-div2","children")], Input("my-id","value"), State("my-id","type")) do value, state
        return value, state
    end
    @test length(app.callbacks) == 1
    @test haskey(app.callbacks, Symbol("..my-div.children...my-div2.children.."))
    @test app.callbacks[Symbol("..my-div.children...my-div2.children..")].func("value", "state") == ("value", "state")

    handler = Dash.make_handler(app)
    test_json = """{"output":"..my-div.children...my-div2.children..","changedPropIds":["my-id.value"],"inputs":[{"id":"my-id","property":"value","value":"test"}], "state":[{"id":"my-id","property":"type","value":"state"}]}"""
    request = HTTP.Request("POST", "/_dash-update-component", [], Vector{UInt8}(test_json))
    response = HTTP.handle(handler, request)
    @test response.status == 200
    resp_obj = JSON2.read(String(response.body))
    @test in(:multi, keys(resp_obj))
    @test resp_obj.response[Symbol("my-div")].children == "test"
    @test resp_obj.response[Symbol("my-div2")].children == "state"

    app = dash()
    app.layout = html_div() do
            dcc_input(id = "my-id", value="initial value", type = "text"),
            html_div(id = "my-div"),
            html_div(id = "my-div2")
        end
    callback!(app, [Output("my-div","children")], Input("my-id","value"), State("my-id","type")) do value, state
        return (value,)
    end
    @test length(app.callbacks) == 1
    @test haskey(app.callbacks, Symbol("..my-div.children.."))
    @test app.callbacks[Symbol("..my-div.children..")].func("value", "state") == ("value", )

    handler = Dash.make_handler(app)
    test_json = """{"output":"..my-div.children..","changedPropIds":["my-id.value"],"inputs":[{"id":"my-id","property":"value","value":"test"}], "state":[{"id":"my-id","property":"type","value":"state"}]}"""
    request = HTTP.Request("POST", "/_dash-update-component", [], Vector{UInt8}(test_json))
    response = HTTP.handle(handler, request)
    @test response.status == 200
    resp_obj = JSON2.read(String(response.body))
    @test in(:multi, keys(resp_obj))
    @test resp_obj.response[Symbol("my-div")].children == "test"

end
@testset "callback! multi output flat" begin
    app = dash()
    app.layout = html_div() do
            dcc_input(id = "my-id", value="initial value", type = "text"),
            dcc_input(id = "my-id2", value="initial value", type = "text"),
            html_div(id = "my-div"),
            html_div(id = "my-div2")
        end
    callback!(app, Output("my-div","children"), Output("my-div2","children"),
                   Input("my-id","value"), Input("my-id", "value"), State("my-id","type")) do value, value2, state
        return value * value2, state
    end
    @test length(app.callbacks) == 1
    @test haskey(app.callbacks, Symbol("..my-div.children...my-div2.children.."))
    @test app.callbacks[Symbol("..my-div.children...my-div2.children..")].func("value", " value2", "state") == ("value value2", "state")

    app = dash()
    app.layout = html_div() do
            dcc_input(id = "my-id", value="initial value", type = "text"),
            dcc_input(id = "my-id2", value="initial value", type = "text"),
            html_div(id = "my-div"),
            html_div(id = "my-div2")
        end
    callback!(app, Output("my-div","children"),
                   Input("my-id","value"), Input("my-id", "value")) do value, value2
        return value * value2
    end
    @test length(app.callbacks) == 1
    @test haskey(app.callbacks, Symbol("my-div.children"))
    @test app.callbacks[Symbol("my-div.children")].func("value", " value2") == "value value2"

end
@testset "callback! checks" begin

    app = dash()

    app.layout = html_div() do
            dcc_input(id = "my-id", value="initial value", type = "text"),
            html_div(id = "my-div"),
            html_div(id = "my-div2")
        end

    callback!(app, Output("my-div","children"), Input("my-id","value")) do value
        return value
    end
    callback!(app, Output("my-div2","children"), Input("my-id","value")) do value
        return "v_$(value)"
    end

    @test length(app.callbacks) == 2
    @test haskey(app.callbacks, Symbol("my-div.children"))
    @test haskey(app.callbacks, Symbol("my-div2.children"))
    @test app.callbacks[Symbol("my-div.children")].func("value") == "value"
    @test app.callbacks[Symbol("my-div2.children")].func("value") == "v_value"

    #output already exists
    @test_throws ErrorException callback!(app, Output("my-div2","children"), Input("my-id","value")) do value
        return "v_$(value)"
    end

    #output same as input
    @test_throws ErrorException callback!(app,
        Output("my-id","value"),
        Input("my-id","value"),
        State("my-id","value")) do value
        return "v_$(value)"
    end

    #output same as input
    @test_throws ErrorException callback!(app,
        Output("my-id","value"),
        [Input("my-id","value"), Input("my-div","children")]) do value
        return "v_$(value)"
    end

    #output same as input
    @test_throws ErrorException callback!(app,
        [Output("my-id","value"), Output("my-div","children")],
        Input("my-id","value")) do value
        return "v_$(value)"
    end

    #empty input
    @test_throws ErrorException callback!(app,
        [Output("my-id","value"), Output("my-div","children")],
        Input[]) do value
        return "v_$(value)"
    end

    #duplicated output
    @test_throws ErrorException callback!(app,
        [Output("my-div","value"), Output("my-div","value")],
        Input("my-id","value")) do value
        return "v_$(value)"
    end

    app = dash()

    app.layout = html_div() do
            dcc_input(id = "my-id", value="initial value", type = "text"),
            dcc_input(id = "my-id2", value="initial value", type = "text"),
            html_div(id = "my-div"),
            html_div(id = "my-div2")
        end
    #empty output
    @test_throws ErrorException callback!(app,
        Input("my-id","value")) do value
        return "v_$(value)"
    end
    #empty input
    @test_throws ErrorException callback!(app,
        Output("my-div2", "children")) do value
        return "v_$(value)"
    end

    #wrong args order
    @test_throws ErrorException callback!(app,
        Input("my-id","value"), Output("my-div", "children")) do value
        return "v_$(value)"
    end

    @test_throws ErrorException callback!(app,
        Output("my-div2", "children"), Input("my-id","value"), Output("my-div", "children")) do value
        return "v_$(value)"
    end

    @test_throws ErrorException callback!(app,
        Output("my-div2", "children"), Input("my-id","value"), State("my-div", "children"), Input("my-id2", "value")) do value, value2
        return "v_$(value)"
    end

end

@testset "pattern-match" begin
    app = dash()
    app.layout = html_div() do
            dcc_input(id = (type="test", index = 1), value="initial value", type = "text"),
            dcc_input(id = (type="test", index = 2), value="initial value", type = "text"),
            html_div(id = (type = "test-out", index = 1)),
            html_div(id = (type = "test-out", index = 2))
        end

    changed_key = """{"index":2,"type":"test"}.value"""
    callback!(app, Output((type = "test-out", index = MATCH), "children"), Input((type="test", index = MATCH), "value")) do value
        context = callback_context()
        @test haskey(context.inputs, changed_key)
        @test context.inputs[changed_key] == "test"
        @test length(context.triggered) == 1
        @test context.triggered[1].prop_id == changed_key
        @test context.triggered[1].value == "test"
        return string(value)
    end
    @test length(app.callbacks) == 1
    out_key = Symbol("""{"index":["MATCH"],"type":"test-out"}.children""")
    @test haskey(app.callbacks, out_key)

    handler = Dash.make_handler(app)
    request = (
        output = string(out_key),
        outputs = (id = (index = 1, type="test_out"), property = "children"),
        changedPropIds = [changed_key],
        inputs = [
            (id = (index=2, type="test"), property = "value", value = "test")
        ]
    )
    test_json = JSON2.write(request)
    request = HTTP.Request("POST", "/_dash-update-component", [], Vector{UInt8}(test_json))
    response = HTTP.handle(handler, request)

    @test response.status == 200
    resp_obj = JSON2.read(String(response.body))
    @test in(:multi, keys(resp_obj))
    @test resp_obj.response.var"{\"index\":1,\"type\":\"test_out\"}".children == "test"

end

@testset "pattern-match ALL" begin
    app = dash()
    app.layout = html_div() do
            dcc_input(id = (type="test", index = 1), value="initial value", type = "text"),
            dcc_input(id = (type="test", index = 2), value="initial value", type = "text"),
            html_div(id = (type = "test-out", index = 1)),
            html_div(id = (type = "test-out", index = 2))
        end

    first_key = """{"index":1,"type":"test"}.value"""
    changed_key = """{"index":2,"type":"test"}.value"""
    callback!(app, Output((type = "test-out", index = ALL), "children"), Input((type="test", index = ALL), "value")) do value
        context = callback_context()
        @test haskey(context.inputs, first_key)
        @test context.inputs[first_key] == "test 1"
        @test haskey(context.inputs, changed_key)
        @test context.inputs[changed_key] == "test"
        @test length(context.triggered) == 1
        @test context.triggered[1].prop_id == changed_key
        @test context.triggered[1].value == "test"
        return value
    end
    @test length(app.callbacks) == 1
    out_key = Symbol("""{"index":["ALL"],"type":"test-out"}.children""")
    @test haskey(app.callbacks, out_key)

    handler = Dash.make_handler(app)
    request = (
        output = string(out_key),
        outputs = [
           [
                (id = (index=1, type="test-out"), property = "children"),
                (id = (index=2, type="test-out"), property = "children")
            ]
        ],
        changedPropIds = [changed_key],
        inputs = [
            [
                (id = (index=1, type="test"), property = "value", value = "test 1"),
                (id = (index=2, type="test"), property = "value", value = "test")
            ]
        ]
    )
    test_json = JSON2.write(request)
    request = HTTP.Request("POST", "/_dash-update-component", [], Vector{UInt8}(test_json))
    response = HTTP.handle(handler, request)
    @test response.status == 200
    resp_obj = JSON2.read(String(response.body))
    @test in(:multi, keys(resp_obj))

    @test resp_obj.response.var"{\"index\":1,\"type\":\"test-out\"}".children =="test 1"
    @test resp_obj.response.var"{\"index\":2,\"type\":\"test-out\"}".children =="test"

end

@testset "invalid multi out" begin
    app = dash()
    app.layout = html_div() do
        html_div(id="a"),
        html_div(id="b"),
        html_div(id="c"),
        html_div(id="d"),
        html_div(id="e"),
        html_div(id="f")
    end

    callback!(app, [Output("b", "children")], Input("a", "children")) do a
        return (1,2)
    end

    #invalid number of outputs
    @test_throws Dash.InvalidCallbackReturnValue Dash.process_callback_call(app, :var"..b.children..", [(id = "b", property = "children")], [(value = "aaa",)], [])

    callback!(app, Output("c", "children"), Output("d", "children"), Input("a", "children")) do a
        return 1
    end
    #result not array
    @test_throws Dash.InvalidCallbackReturnValue Dash.process_callback_call(
            app, Symbol("..c.children...d.children.."),
            [(id = "b", property = "children")], [(value = "aaa",)], []
            )

    app = dash()
    app.layout = html_div() do
        html_div(id=(index = 1,)),
        html_div(id=(index = 2,)),
        html_div(id=(index = 3,))
    end

    callback!(app, Output((index = ALL,), "children"), Input((index = ALL,), "css")) do inputs
        return [1,2]
    end

    #pattern-match output elements length not match to specs
    @test_throws Dash.InvalidCallbackReturnValue Dash.process_callback_call(
            app, Symbol("""{"index":["ALL"]}.children"""),
            [[
                (id = """{"index":1}""", property = "children"),
                (id = """{"index":2}""", property = "children"),
                (id = """{"index":3}""", property = "children"),
            ]],
            [[(value = "aaa",), (value = "bbb",), (value = "ccc",)]], []
            )

    app = dash()
    app.layout = html_div() do
        html_div(id=(index = 1,)),
        html_div(id=(index = 2,)),
        html_div(id=(index = 3,))
    end

    callback!(app, Output((index = ALL,), "children"), Input((index = ALL,), "css")) do inputs
        return 1
    end

    #pattern-match output element not array
    @test_throws Dash.InvalidCallbackReturnValue Dash.process_callback_call(
            app, Symbol("""{"index":["ALL"]}.children"""),
            [[
                (id = """{"index":1}""", property = "children"),
                (id = """{"index":2}""", property = "children"),
                (id = """{"index":3}""", property = "children"),
            ]],
            [[(value = "aaa",), (value = "bbb",), (value = "ccc",)]], []
            )
end

@testset "clientside callbacks function" begin
    app = dash()
    app.layout = html_div() do
            dcc_input(id = "my-id", value="initial value", type = "text"),
            html_div(id = "my-div")
        end

    callback!(ClientsideFunction("namespace", "func_name"),app, Output("my-div","children"), Input("my-id","value"))

    @test length(app.callbacks) == 1
    @test haskey(app.callbacks, Symbol("my-div.children"))
    @test app.callbacks[Symbol("my-div.children")].func isa ClientsideFunction
    @test app.callbacks[Symbol("my-div.children")].func.namespace == "namespace"
    @test app.callbacks[Symbol("my-div.children")].func.function_name == "func_name"

    handler = make_handler(app)
    request = HTTP.Request("GET", "/_dash-dependencies")
    resp = HTTP.handle(handler, request)
    deps = JSON2.read(String(resp.body))

    @test length(deps) == 1
    cb = deps[1]
    @test cb.output == "my-div.children"
    @test cb.inputs[1].id == "my-id"
    @test cb.inputs[1].property == "value"
    @test :clientside_function in keys(cb)
    @test cb.clientside_function.namespace == "namespace"
    @test cb.clientside_function.function_name == "func_name"
end
@testset "clientside callbacks string" begin
    app = dash()
    app.layout = html_div() do
            dcc_input(id = "my-id", value="initial value", type = "text"),
            html_div(id = "my-div")
        end

    callback!(
        """
        function(input_value) {
            return (
                parseFloat(input_value_1, 10)
            );
        }
        """
        , app,
        Output("my-div", "children"),
        Input("my-id", "value")
        )

    @test length(app.callbacks) == 1
    @test haskey(app.callbacks, Symbol("my-div.children"))
    @test app.callbacks[Symbol("my-div.children")].func isa ClientsideFunction
    @test app.callbacks[Symbol("my-div.children")].func.namespace == "_dashprivate_my-div"
    @test app.callbacks[Symbol("my-div.children")].func.function_name == "children"
    @test length(app.inline_scripts) == 1
    @test occursin("clientside[\"_dashprivate_my-div\"]", app.inline_scripts[1])
    @test occursin("ns[\"children\"]", app.inline_scripts[1])

    handler = make_handler(app)
    request = HTTP.Request("GET", "/_dash-dependencies")
    resp = HTTP.handle(handler, request)
    deps = JSON2.read(String(resp.body))

    @test length(deps) == 1
    cb = deps[1]
    @test cb.output == "my-div.children"
    @test cb.inputs[1].id == "my-id"
    @test cb.inputs[1].property == "value"
    @test :clientside_function in keys(cb)
    @test cb.clientside_function.namespace == "_dashprivate_my-div"
    @test cb.clientside_function.function_name == "children"
    request = HTTP.Request("GET", "/")
    resp = HTTP.handle(handler, request)
    body = String(resp.body)
    @test occursin("clientside[\"_dashprivate_my-div\"]", body)
    @test occursin("ns[\"children\"]", body)
end

ext: jl


text:
module 2022_01_comext

# Write your package code here.

end


text:
"""
we need to give back number of false positive and false negatives and min,max x,y,x of block containing all data
IMPORTANT - we require at least 7 in y dim and 32 in x dim of block thread
adapted from https://github.com/JuliaGPU/CUDA.jl/blob/afe81794038dddbda49639c8c26469496543d831/src/mapreduce.jl
goldBoolGPU3d - array holding 3 dimensional data of gold standard bollean array
segmBoolGPU3d - array with 3 dimensional the data we want to compare with gold standard
reducedGold - the smallest boolean block (3 dim array) that contains all positive entris from both masks
reducedSegm - the smallest boolean block (3 dim array) that contains all positive entris from both masks
numberToLooFor - number we will analyze whether is the same between two sets
loopNumbYdim - number of times the single lane needs to loop in order to get all needed data - in this kernel it will be exactly a y dimension of a slice
xdim - length in x direction of source array
loopNumbXdim - in case the x dim will be bigger than number of threads we will create second inner loop
cuda arrays holding just single value wit atomically reduced result
,fn,fp
,minxRes,maxxRes
,minyRes,maxyRes
,minZres,maxZres

"""
function getBoolCubeKernel(goldBoolGPU3d
        ,segmBoolGPU3d
        ,reducedGoldA
        ,reducedSegmA
        ,reducedGoldB
        ,reducedSegmB
        ,loopNumbYdim::UInt16
        ,xdim::UInt16
        ,loopNumbXdim::UInt16
        ,numberToLooFor::T
        ,fn::CuDeviceVector{UInt32, 1}
        ,fp::CuDeviceVector{UInt32, 1}
        ,minxRes::CuDeviceVector{UInt32, 1}
        ,maxxRes::CuDeviceVector{UInt32, 1}
        ,minyRes::CuDeviceVector{UInt32, 1}
        ,maxyRes::CuDeviceVector{UInt32, 1}
        ,minzRes::CuDeviceVector{UInt32, 1}
        ,maxzRes::CuDeviceVector{UInt32, 1}
        ,dataBdim
) where T

   anyPositive = false # true If any bit will bge positive in this array - we are not afraid of data race as we can set it multiple time to true
   #creates shared memory and initializes it to 0
   shmemSum = @cuStaticSharedMem(Float32,(32,2))
   interMediateFPFN=@cuStaticSharedMem(UInt32,(2))
   #incrementing appropriate number of times

    #1 - false negative; 2- false positive
    locArr= (Float32(0.0), Float32(0.0))
    offsetIter= UInt8(1)
    #needed to get the borders in metadata terms - so blocks that have sth of our intrest
    minX =@cuStaticSharedMem(Float32, 1)
    maxX= @cuStaticSharedMem(Float32, 1)
    minY = @cuStaticSharedMem(Float32, 1)
    maxY= @cuStaticSharedMem(Float32, 1)
    minZ = @cuStaticSharedMem(Float32, 1)
    maxZ= @cuStaticSharedMem(Float32, 1)
    #those will be used in order to establish how much results we have in each border queue
    #will be accesses=d atomically from threads - to be experimented upon wheather we have better way
    topFP =@cuStaticSharedMem(UInt16, 1)
    topFN =@cuStaticSharedMem(UInt16, 1)

    bottomFP =@cuStaticSharedMem(UInt16, 1)
    bottomFN =@cuStaticSharedMem(UInt16, 1)

    leftFP =@cuStaticSharedMem(UInt16, 1)
    leftFN =@cuStaticSharedMem(UInt16, 1)

    rightFP =@cuStaticSharedMem(UInt16, 1)
    rightFN =@cuStaticSharedMem(UInt16, 1)

    anteriorFP =@cuStaticSharedMem(UInt16, 1)
    anteriorFN =@cuStaticSharedMem(UInt16, 1)

    posteriorFP =@cuStaticSharedMem(UInt16, 1)
    posteriorFN =@cuStaticSharedMem(UInt16, 1)

    isAnyPositive = @cuStaticSharedMem(Bool, 1)
    #resetting
    minX[1]= Float32(1110.0)
    maxX[1]= Float32(0.0)
    minY[1]= Float32(1110.0)
    maxY[1]= Float32(0.0)
    minZ[1]= Float32(1110.0)
    maxZ[1]= Float32(0.0)

    topFP = UInt16(0)
    topFN = UInt16(0)
    bottomFP= UInt16(0)
    bottomFN= UInt16(0)
    leftFP= UInt16(0)
    leftFN= UInt16(0)
    rightFP= UInt16(0)
    rightFN= UInt16(0)
    anteriorFP= UInt16(0)
    anteriorFN= UInt16(0)
    posteriorFP= UInt16(0)
    posteriorFN= UInt16(0)

    isAnyPositive[1]= false

    sync_threads()
    #we need nested x,y,z iterations so we will iterate over the matadata and on its basis over the  data in the main arrays
    #first loop over the metadata
    #dataBdim - indicats dimensions of data blocks
    @iter3dOuter(xname,yName , zName, loopXMeta,loopYMeta,loopZmeta,
         begin
         #inner loop is over the data indicated by metadata
         @iter3dInMeta(xOuter*dataBdim[1] ,yOuter*dataBdim[2], zzOuter*dataBdim[3],true,zdim ,inBlockLoopX,inBlockLoopY,inBlockLoopZ
                         ,begin
                                boolGold=    goldBoolGPU3d[x,y,z]==numberToLooFor
                                boolSegm=    segmBoolGPU3d[x,y,z]==numberToLooFor
                                fpXOrFn = (boolGold  ⊻ boolSegm)
                                #we need to also collect data about how many fp and fn we have in main part and borders
                                #important in case of corners we will first analyze z and y dims and z dim on last resort only !
                                if(fpXOrFn)
                                    #here we store total count first false nehative, second false positive
                                    @inbounds locArr[boolGold+ boolSegm*2]+=1

                                    if(xdim ==1) #left
                                        incrementAtomFPifTrue(boolSegm,leftFP,leftFN)
                                    elseif(xdim == inBlockLoopDims[1] )  # right
                                        incrementAtomFPifTrue(boolSegm,rightFP, rightFN)
                                    elseif(ydim == 0 )  # posterior
                                        incrementAtomFPifTrue(boolSegm,posteriorFP , posteriorFN)
                                    elseif(ydim == inBlockLoopDims[2] )  # anterior
                                        incrementAtomFPifTrue(boolSegm,anteriorFP, anteriorFN)
                                    elseif(zdim == 0 )  # top
                                        incrementAtomFPifTrue(boolSegm,topFP,topFN)
                                    elseif(zdim == inBlockLoopDims[3] )  # bottom
                                        incrementAtomFPifTrue(boolSegm,bottomFP, bottomFN)
                                    end
                                    isAnyPositive[1]= true #- we just mark that there was some fp or fn in this block
                                end#if fpXOrFn

                                #in case some is positive we can go futher with looking for max,min in dims and add to the new reduced boolean arrays waht we are intrested in
                                if(boolGold  || boolSegm)
                                    #passing data to new arrays needed for running final algorithm
                                    reducedGoldA[x,y,z]=boolGold
                                    reducedSegmA[x,y,z]=boolSegm
                                    reducedGoldB[x,y,z]=boolGold
                                    reducedSegmB[x,y,z]=boolSegm
                                end#if boolGold  || boolSegm
                            end)#ex

                #now we are just after we iterated over a single data block  we need to update the data about block metadata

                #we save data about border data blocks
                sync_threads()
                @redWitAct(offsetIter,shmemSum,  locArr[1],+,     locArr[2],+ )
                sync_threads()

                #now we need to reduce locArr - to get total amount of fp and fn

               #save the data about number of fp and fn of this block and accumulate also this sum for global sum
                @ifXY 1 1 if(isAnyPositive[1]) minX[1]= min(minX[1],xOuter) end
                @ifXY 2 1 if(isAnyPositive[1]) maxX[1]= max(maxX[1],xOuter) end
                @ifXY 3 1 if(isAnyPositive[1]) minY[1]= min(minY[1],yOuter) end
                @ifXY 4 1 if(isAnyPositive[1]) maxY[1]= max(maxY[1],yOuter) end
                @ifXY 5 1 if(isAnyPositive[1]) minZ[1]= min(minZ[1],zOuter) end
                @ifXY 6 1 if(isAnyPositive[1]) maxZ[1]= max(maxZ[1],zOuter) end
                @ifXY 7 1 if(isAnyPositive[1]) isAnyPositive[1]= false end #reset

                @ifXY 1 2 if(isAnyPositive[1]) setMetaLeftFP(metaData,leftFP[1]) end
                @ifXY 2 2 if(isAnyPositive[1]) setMetaLeftFN(metaData,leftFN[1]) end

                @ifXY 3 2 if(isAnyPositive[1]) setMetaRightFP(metaData,rightFP[1]) end
                @ifXY 4 2 if(isAnyPositive[1]) setMetaRightFN(metaData,rightFN[1]) end

                @ifXY 5 2 if(isAnyPositive[1]) setMetaPosteriorFP(metaData,posteriorFP[1]) end
                @ifXY 6 2 if(isAnyPositive[1]) setMetaPosteriorFN(metaData,posteriorFN[1]) end

                @ifXY 7 2 if(isAnyPositive[1]) setMetaAnteriorFP(metaData,anteriorFP[1]) end
                @ifXY 8 2 if(isAnyPositive[1]) setMetaAnteriorFN(metaData,anteriorFN[1]) end

                @ifXY 9 2 if(isAnyPositive[1]) setMetaTopFP(metaData,topFP[1]) end
                @ifXY 10 2 if(isAnyPositive[1]) setMetaTopFN(metaData,topFN[1]) end

                @ifXY 11 2 if(isAnyPositive[1]) setMetaBottomFP(metaData,bottomFP[1]) end
                @ifXY 12 2 if(isAnyPositive[1])  setMetaBottomFN(metaData,bottomFN[1]) end

                @ifXY 1 3 if(isAnyPositive[1]) setMetaDataTotalFpCount(metaData,shmemSum[1,2], xOuter,yOuter,zOuter) end
                @ifXY 2 3 if(isAnyPositive[1]) setMetaDataTotalFnCount(metaData,shmemSum[1,1], xOuter,yOuter,zOuter) end
                #now in order to get total count  of fp and fn per block  we will do this subtraction

                @ifXY 3 3 if(isAnyPositive[1]) setMetaDataMainFpCount(metaData,shmemSum[1,2] -leftFP[1] - rightFP[1]- posteriorFP[1] -anteriorFP[1]- topFP[1] -bottomFP[1] , xOuter,yOuter,zOuter) end
                @ifXY 4 3 if(isAnyPositive[1]) setMetaDataMainFnCount(metaData,shmemSum[1,1] -leftFN[1] - rightFN[1] -posteriorFN[1]- anteriorFN[1]- topFN[1] -bottomFN[1]   , xOuter,yOuter,zOuter) end
                @ifXY 4 4 setMetaDataXYZ(metaData, xOuter,yOuter,zOuter  )

                #adding the total fp and fn of a block to total result per thread block - later it will be pushed as total amiunt of fp and fn
                @ifXY 5 3 interMediateFPFN[1] +=shmemSum[1,1]
                @ifXY 6 3 interMediateFPFN[2] +=shmemSum[1,2]
                sync_threads()

                #cleaning
                clearSharedMemWarpLong(shmemSum, UInt8(6), Float32(0.0))
                locArr= (Float32(0.0), Float32(0.0))

                #set the x,y,z coordinates - so we will able to query it efficiently also with linear index
                #what is important later as we will use only part of meta data this indicies will need to be updated

     end) #outer loop
                #consider ceating tuple structure where we will have  number of outer tuples the same as z dim then inner tuples the same as y dim and most inner tuples will have only the entries that are fp or fn - this would make us forced to put results always in correct spots

        # outer loop expession  )

    #in order to have global data

    @ifXY 1 1 atomicMinSet(minxRes[],minX[1])
    @ifXY 2 1 atomicMaxSet(maxxRes[],maxX[1])

    @ifXY 3 1 atomicMinSet(minyRes[],minY[1])
    @ifXY 4 1 atomicMaxSet(maxyRes[],maxY[1])

    @ifXY 5 1 atomicMinSet(minzRes[],minY[1])
    @ifXY 6 1 atomicMaxSet(maxzRes[],maxZ[1])

    @ifXY 7 1 atomicAdd(fn[],interMediateFPFN[2])
    @ifXY 8 1 atomicAdd(fp[],interMediateFPFN[1])
   return
   end

   """
   atomically increment by one fp if true otherwise increment fn by 1
   """
   function incrementAtomFPifTrue(bool,fp,fn)
        if(bool)
            CUDA.@atomic @inbounds fp[]+=1
        end
        CUDA.@atomic @inbounds fn[]+=1

   end


text:
###################################
# variables shared across scripts #
###################################
ALL_PROBLEM_NAMES = ["qap10", "qap15", "nug08-3rd", "nug20"]

# restart lengths to try for eahc problem
RESTART_LENGTHS_DICT = Dict(
  "qap10" => 4 .^ collect(1:9),
  "qap15" => 4 .^ collect(1:9),
  "nug08-3rd" => 4 .^ collect(1:9),
  "nug20" => 4 .^ collect(1:9),
)

# other parameter values
ITERATION_LIMIT = 500000

max_stars_repo_name: pedersor/google-research

max_issues_repo_name: pedersor/google-research

max_forks_repo_name: pedersor/google-research


text:
using BlsData
using Test

let
    b = Bls()
    data = get_data(b, "PRS85006092")
    data = get_data(b, ["PRS85006092", "LNS11000000"])
    data = get_data(b, "PRS85006092"; startyear=1990)
    data = get_data(b, ["PRS85006092", "LNS11000000"]; endyear=2014)
    data = get_data(b, ["PRS85006092", "LNS11000000"]; startyear=1980, endyear=2014)
    data = get_data(b, "PRS85006092"; catalog=true)

    nothing
end

let
    b = Bls("test")

    nothing
end

let
    # Test empty BlsSeries
    empty_response = BlsData.EMPTY_RESPONSE()
    @test isempty(empty_response)

    nothing
end


text:
import HomotopyContinuation;	const HC = HomotopyContinuation
import DynamicPolynomials;		const DP = DynamicPolynomials
using FastClosures
"""
	W0(t, tf, χ, args)
Calculate the correction operator, on the n subspace, associated to the a pulse
on the j subspace.
#Arguments
- t: time.
- tf: duration of the snap gate.
- χ: dispersive coupling constant
- n: subspace on which we calculate the correction operator.
- args: tuple with specification about the control
args[1] is an int that tells to which Pauli matrix the control couples to.
args[2] is the Fourier order of envelope.
args[3] is the envelope (fs or fc)
"""
function W0(t, tf, χ, args)
	control = args[1]
	order = args[2]
	envelope = args[3]

	W0t = zeros(4)
	if control == 1
		W0t[1] = envelope(t, tf, order)
	elseif control == 2
		W0t[2] = -envelope(t, tf, order)
	else
		W0t[3] = envelope(t, tf, order)/2
	end
	return W0t
end
"""
Calculate the correction operator W0 in the interaction picture.
#Arguments
- t: time.
- tf: duration of the snap gate.
- χ: dispersive coupling constant
- levels: the energy levels of the cavity that are resonantly driven by the ideal
Hamiltonian H0.
- α: the angle of the driving pulse in the second half of the time evolution.
- n: subspace on which we calculate the correction operator.
- W0_args: array with specification about the control; cf. W0.
- N: the size of the Hilbert space of the cavity.
"""
function W0I(t, tf, χ, levels, α, n, W0_args, N)
	W0t = W0(t, tf, χ, W0_args)
	indice = findall(k->k==n, levels)
	if length(indice) == 0
		W0It = W0t
	else
		T = tf/2.0
		if t < T + 10.0*eps()
			intfx, intfy = intf(t, T, zero(t))
			vn = (intfx, intfy, zero(t), zero(t))
			unitary_tranformation!(W0t, vn)
		else
			αn = α[indice[1]]
			intfx, intfy = intf(t - T, T, αn)
			vn = (intfx, intfy, zero(t), zero(t))
			unitary_tranformation!(W0t, vn)
			intfx, intfy = intf(T, T, zero(t))
			vn = (intfx, intfy, zero(t), zero(t))
			unitary_tranformation!(W0t, vn)
		end
	end
	return W0t
end
function delta_coeffs(system_params; eps=0.0)
	N = system_params["N"]
	M = zeros(3*N, 1)
	control = 3
	order = 0
	W0_args = (control, order, (t, tf, order)-> 1.0)
	indice = 1
	V(t) = zeros(4, N)
	for n=0:N-1
		Mn = view(M, n*3 + 1: (n + 1)*3)#zeros(3)
		linear_coeffs_col!(Mn, n, W0_args, indice, V, system_params)
		#M[n*3+1:(n+1)*3] .= Mn
	end

	#set elements smaller than eps to zero
	M[abs.(M) .< eps] .= 0.0
	return M
end
"""
Calculate the matrix with the coefficients of the free parameters in the
subspace n of the bosonic mode. The product of the coefficients with the
respective free parameters is the 1st order term of the Magnus expansion of the
correction operator. Each column is associated to one free parameter.
#Arguments
- n: the subspace of the bosonic mode.
- max_order: the maximum value of the parameter 'order'.
- V: The Hamiltonian in the interaction picture. It includes lower order
correction operators.
- system_params: dictionary with the system parameters.
#Keyword arguments
- eps: if the absolute value of any coefficients is smaller than eps, it will
be set to zero. The default value is 0.
"""
function linear_coeffs(n, max_order, V, system_params; eps::Real=0.0)
	M = zeros(3, 4*max_order)
	indice = 1
	for control = 1: 2
		for order = 1: max_order
			W0_args = (control, order, fc)
			linear_coeffs_col!(M, n, W0_args, indice, V, system_params)
			W0_args = (control, order, fs)
			linear_coeffs_col!(M, n, W0_args, indice + 1, V, system_params)
			indice += 2
		end
	end
	#set elements smaller than eps to zero
	M[abs.(M) .< eps] .= 0.0
	return M
end
"""
Update the elements of the matrix with the coefficients of the free parameters.
#Arguments
- M: the matrix with the coefficients of the free parameters.
- n: the subspace of the bosonic mode.
- W0_args: array with specifications about the control; cf. W0.
- indice: indice of the column being updated.
- V: The Hamiltonian in the interaction picture. It includes lower order
correction operators.
- system_params: dictionary with the system parameters.
"""
function linear_coeffs_col!(M, n, W0_args, indice, V, system_params)
	tf = system_params["tf"]
	χ = system_params["chi"]
	levels = system_params["levels"]
	α = system_params["alpha"]
	N = system_params["N"]

	funcV = @closure t-> V(t)[:, n + 1]
	funcW = @closure t-> W0I(t, tf, χ, levels, α, n, W0_args, N)

	f0 = zeros(4, 3)
	tspan = (zero(tf), tf)
	args = (funcV, funcW, zeros(4), zeros(4))
	prob = DE.ODEProblem(df!, f0, tspan, args)
	sol = DE.solve(prob, DE.Vern9(), saveat=[tf], reltol=1e-10, abstol=1e-12)
	sol = sol(tf)

	M[:, indice] .= -(sol[1: 3, 2] .+ sol[1: 3, 3])
end
"""
Calculate the tensor with the coefficients of the quadratic products of the free
parameters in the subspace n of the bosonic mode. The coefficients (times the
quadratic products of the free parameters) appear in the 2nd term of the Magnus
expansion of the modified Hamiltonian (original Hamiltonian + correction
Hamiltonian).
#Arguments
- n: the subspace of the bosonic mode.
- max_order: the maximum value of the parameter 'order'.
- system_params: dictionary with the system parameters.
#Keyword arguments
- eps: if the absolute value of any coefficients is smaller than eps, it will
be set to zero. The default value is 0.
"""
function nonlinear_coeffs(n, max_order, system_params; eps::Real = 0.0)
	M = zeros(3, 4*max_order, 4*max_order)
	indice = 1
	for control = 1: 2
		for order = 1: max_order
			W0_args1 = (control, order, fc)
			nonlinear_coeffs_row!(M, n, max_order, W0_args1, indice, system_params)
			W0_args1 = (control, order, fs)
			nonlinear_coeffs_row!(M, n, max_order, W0_args1, indice+1, system_params)
			indice += 2
		end
	end
	#set elements smaller than eps to zero
	M[abs.(M) .< eps] .= 0.0
	return M
end
"""
Update the elements of the tensor with the coefficients of the quadratic products
of the free parameters.
#Arguments
- M: the tensor to be updated.
- n: the subspace of the bosonic mode.
- max_order: the maximum value of the parameter 'order'.
- W0_args1: array with specifications about control 1; cf. W0.
- indice1: indice of the elements being updated (2nd indice of the tensor).
- system_params: dictionary with the system parameters.
"""
function nonlinear_coeffs_row!(M, n, max_order, W0_args1, indice1, system_params)
	indice2 = 1
	for control = 1: 2
		for order = 1: max_order
			W0_args2 = (control, order, fc)
			index = (indice1, indice2)
			nonlinear_coeffs_col!(M, n, W0_args1, W0_args2, index, system_params)
			W0_args2 = (control, order, fs)
			index = (indice1, indice2 + 1)
			nonlinear_coeffs_col!(M, n, W0_args1, W0_args2, index, system_params)
			indice2 += 2
		end
	end
end
"""
Update the elements of the tensor with the coefficients of the quadratic products
of the free parameters.
#Arguments
- M: the tensor to be updated.
- n: the subspace of the bosonic mode.
- max_order: the maximum value of the parameter 'order'.
- W0_args1 and W0_args2: array with specifications about controls 1 and 2; cf. W0.
- index: index of the elements being updated (2nd and 3rd index of the tensor).
- system_params: dictionary with the system parameters.
"""
function nonlinear_coeffs_col!(M, n, W0_args1, W0_args2, index, system_params)
	tf = system_params["tf"]
	χ = system_params["chi"]
	levels = system_params["levels"]
	α = system_params["alpha"]
	N = system_params["N"]

	funcW1 = @closure t-> W0I(t, tf, χ, levels, α, n, W0_args1, N)
	funcW2 = @closure t-> W0I(t, tf, χ, levels, α, n, W0_args2, N)

	f0 = zeros(4, 3)
	tspan = (zero(tf), tf)
	args = (funcW1, funcW2, zeros(4))
	prob = DE.ODEProblem(df2!, f0, tspan, args)
	sol = DE.solve(prob, DE.Vern9(), saveat=[tf], reltol=1e-10, abstol=1e-12)
	sol = sol(tf)

	M[:, index[1], index[2]] .= -sol[1: 3, 3]
end
"""
	df!(dOmega, Omega, args, t)
The integral of this function is related to the coefficients of the linear terms in the
polynomial system of equations.
"""
function df!(dOmega, Omega, args, t)
	V = args[1]
	W0 = args[2]
	temp1 = args[3]
	temp2 = args[4]
	#k=1
	dOmega[:, 1] .= -V(t)
	dOmega[:, 2] .= -W0(t)
	#k=2
	comm!(temp1, Omega[:, 2], dOmega[:, 1])
	comm!(temp2, Omega[:, 1], dOmega[:, 2])
	dOmega[:, 3] .= (temp1 .+ temp2)./2
	return nothing
end
"""
	df2!(dOmega, Omega, args, t)
The integral of this function is related to the coefficients of the nonlinear terms in the
polynomial system of equations.
"""
function df2!(dOmega, Omega, args, t)
	W01 = args[1]
	W02 = args[2]
	temp1 = args[3]::Array{Float64,1}
	#k=1
	dOmega[:,1] .= -W01(t)
	dOmega[:,2] .= -W02(t)
	#k=2
	comm!(temp1, Omega[:,2], dOmega[:,1])
	dOmega[:,3] .= temp1./2
	return nothing
end
"""
	correction_coeffs(Mlinear, Mnonlinear, Omega, max_order)
Generate the system of polynomial equations in a given number subspace and find its
solutions using homotopy continuation. Lagrange multipliers are used to find the solutions
with smallest norm.
"""
function correction_coeffs(Mlinear, Mnonlinear, Omega, max_order)
	HC.@polyvar x[1: 4*max_order]

	p1 = sum(Mlinear[1, :].*x) + sum(x.*(Mnonlinear[1, :, :]*x)) - Omega[1]
	p2 = sum(Mlinear[2, :].*x) + sum(x.*(Mnonlinear[2, :, :]*x)) - Omega[2]
	p3 = sum(Mlinear[3, :].*x) + sum(x.*(Mnonlinear[3, :, :]*x)) - Omega[3]
	norm2 = sum(x.*x)
	# Jacobians
	Jp1 = HC.differentiate(p1, x[1: end])
	Jp2 = HC.differentiate(p2, x[1: end])
	Jp3 = HC.differentiate(p3, x[1: end])
	Jnorm = HC.differentiate(norm2, x[1: end])
	# define Langrange equations
	HC.@polyvar l[1: 3]
	a = Jnorm - l[1]*Jp1 - l[2]*Jp2 - l[3]*Jp3
	C = vcat(a, [p1, p2, p3])

	res = HC.solve(C, show_progress=false)
	real_sols = HC.real_solutions(res)
	if length(real_sols) == 0
		res = HC.solve(C, show_progress=true)
		real_sols = HC.real_solutions(res)
	end

	real_sols = map(p -> p[1: 4*max_order], real_sols)

	coeffs = real_sols[1]
	for i = 2: length(real_sols)
		if norm(coeffs) > norm(real_sols[i])
			coeffs = real_sols[i]
		end
	end
	return coeffs
end
"""
	W(t, tf, χ, coeffs, max_order, Ncorr, N)
Correction Hamiltonian in the qubit frame.
# Arguments
- `t::Real`: time.
- `tf::Real`: gate time,
- `χ::Tuple`: Tuple with first, second, and third order dispersive couplings.
- `coeffs::Array`: corection coefficients.
- `max_order::Integer`: maximum number of harmonics of the correction pulse.
- `Ncorr::Integer`: bosonic subspace that should be corrected.
- `N::Integer`: number of dimensions of the truncated bosonic space.
"""
function W(t, tf, χ, coeffs, max_order, Ncorr, N)
	Wt = zeros(4, N)
	χef = (χ[1], χ[2]/2, χ[3]/6)
	for j = 0: Ncorr - 1
		coeffsj = view(coeffs, 4*max_order*j + 1:4*max_order*(j+1))
		Wj!(Wt, t, tf, χef, j, coeffsj, max_order, N)
	end
	Wt[3, :] .= coeffs[end]/2
	return Wt
end
"""
	Wj!(Wt, t, tf, χef, j, coeffsj, max_order, N)
Update `Wt` in the bosonic subspace associated with the `j`th number state.
"""
function Wj!(Wt, t, tf, χef, j, coeffsj, max_order, N)
	gx, gy = g(t, tf, coeffsj, max_order)
	for n = 0: N - 1
		dφn = χef[1]*n + χef[2]*n*(n - 1) + χef[3]*n*(n - 1)*(n - 2)
		dφj = χef[1]*j + χef[2]*j*(j - 1) + χef[3]*j*(j - 1)*(j - 2)
		φ = (dφn - dφj)*t
		cosφ = cos(φ)
		sinφ = sin(φ)
		Wt[1, n + 1] += + gx*cosφ - gy*sinφ
		Wt[2, n + 1] += - gx*sinφ - gy*cosφ
	end
	return nothing
end
"""
	g(t, tf, coeffs, max_order)
Correction pulse envelope.
# Arguments
- `t::Real`: time.
- `tf::Real`: gate time,
- `coeffs::Array`: corection coefficients.
- `max_order::Integer`: maximum number of harmonics of the correction pulse.
"""
function g(t, tf, coeffs, max_order)
	gx, gy = zero(t), zero(t)
	for order = 1: max_order
		cos_order = fc(t, tf, order)
		sin_order = fs(t, tf, order)

		l = 2*(order - 1) + 1
		gx += coeffs[l]*cos_order + coeffs[l + 1]*sin_order
		l += 2*max_order
		gy += coeffs[l]*cos_order + coeffs[l + 1]*sin_order
	end
	return gx, gy
end
"""
	WI(t, tf, χ, coeffs, max_order, levels, α, Ncorr, N)
Correction Hamiltonian in the interaction picture.
# Arguments
- `t::Real`: time.
- `tf::Real`: gate time,
- `χ::Tuple`: Tuple with first, second, and third order dispersive couplings.
- `coeffs::Array`: corection coefficients.
- `max_order::Integer`: maximum number of harmonics of the correction pulse.
- `levels::StaticArray`: levels in which one wants to imprint a phase.
- `α::Tuple`: the ''drinving angles'', which corresponds to π .- φ, where φ is the array with
angles that one wishes to imprint.
- `Ncorr::Integer`: bosonic subspace that should be corrected.
- `N::Integer`: number of dimensions of the truncated bosonic space.
"""
function WI(t, tf, χ, coeffs, max_order, levels, α, Ncorr, N)
	T = tf/2
	Wt = W(t, tf, χ, coeffs, max_order, Ncorr, N)
	v = zeros(4, N)
	if t < T + 10.0*eps()
		v[1: 2, :] .= intf(t, T, levels, α.*0, N)
		unitary_tranformation!(Wt, v, N)
	else
		v[1: 2, :] .= intf(t - T, T, levels, α, N)
		unitary_tranformation!(Wt, v, N)
		v[1: 2, :] .= intf(T, T, levels, α.*0, N)
		unitary_tranformation!(Wt, v, N)
	end
	return Wt
end
"""
	correction_coeffs2(Mlinear, Omega, max_order)
Calculate the correction coefficients using the fully linear method.
"""
function correction_coeffs2(Mlinear, Omega, max_order)
	if max_order > 1
		println("Danger!!! Correct correction_coeffs()")
	end
	coeffs = pinv(Mlinear)*Omega[1: 3]
	return coeffs
end
"""
	correction_matrix(max_order, system_params)
Correction matrix for the (linear) harmonic Magnus method.
"""
function correction_matrix(max_order, system_params)
	N = system_params["N"]
	χ = system_params["chi"]
	tf = system_params["tf"]
	levels = system_params["levels"]
	α = system_params["alpha"]
	Ncorr = system_params["Ncorr"]

	M = zeros(3*N, 4*max_order*N + 1)
	for i = 1: 4*max_order*N + 1
		coeffs = zeros(4*max_order*N + 1)
		coeffs[i] = 1.0
		func = @closure t-> WI(t, tf, χ, coeffs, max_order, levels, α, Ncorr, N)

		f0 = zeros(4, N)
		tspan = (zero(tf), tf)
		args = (func,)
		prob = DE.ODEProblem(df3!, f0, tspan, args)
		sol = DE.solve(prob, DE.Vern9(), saveat=[tf], reltol=1e-10, abstol=1e-12)
		sol = sol(tf)

		M[:, i] = reshape(sol[1: 3, :], 3*N)
	end
	return M
end
"""
	df3!(dOmega, Omega, args, t)
Auxiliary function to compute the integral of a function W0(t).
"""
function df3!(dOmega, Omega, args, t)
	W0 = args[1]
	dOmega .= W0(t)
	return nothing
end
"""
	U(system_params, σ, coeffs, max_order; frame::String = "qubit", matrix::Bool = false)
Calculate the evolution operator.
# Keyword arguments
- `frame::String`: The frame in which the evolution must be calculated. The frame can be the
`qubit` frame, the `lab` frame, or the `interaction` frame. Default value is `qubit`.
- `matrix::Bool`: if `true`, U is returned as a matrix. If `false`, U is returned as a
tensor. Default value is `false`.
"""
function U(system_params, σ, coeffs, max_order; frame::String = "qubit",
		   matrix::Bool = false)
	tf = system_params["tf"]
	χ = system_params["chi"]
	levels = system_params["levels"]
	α = system_params["alpha"]
	Ncorr = system_params["Ncorr"]
	N = system_params["N"]

	if frame == "qubit"
		Faux = @closure t-> (F(t, tf, χ, levels, α, N)
							 .+ W(t, tf, χ, coeffs, max_order, Ncorr, N))
	elseif frame == "interaction"
		Faux = @closure t-> (VI(t, tf, χ, levels, α, N)
							 .+ WI(t, tf, χ, coeffs, max_order, levels, α, Ncorr, N))
	elseif frame == "lab"
		Faux = @closure t-> Flab_corrected(t, tf, χ, levels, α, coeffs, max_order, Ncorr, N)
	end

	sol = solveU(H, χ, Faux, tf, N, σ)
	if matrix == false
		return sol
	else
		return tensor2matrix(sol)
	end
end
"""
	get_coeffs(magnus_order, max_order, Ncorr, system_params)
Calculate the correction coefficients.
"""
function get_coeffs(magnus_order, max_order, Ncorr, system_params)
	N = system_params["N"]
	tf = system_params["tf"]
	χ = system_params["chi"]
	levels = system_params["levels"]
	α = system_params["alpha"]

	Faux = @closure t-> VI(t, tf, χ, levels, α, N)
	coeffs = get_coeffs_aux(Faux, magnus_order, max_order, system_params)
	return coeffs
end
"""
	get_coeffs(coeffs, magnus_order, max_order, Ncorr, system_params)
Calculate the correction coefficients.
"""
function get_coeffs(coeffs, magnus_order, max_order, Ncorr, system_params)
	N = system_params["N"]
	tf = system_params["tf"]
	χ = system_params["chi"]
	levels = system_params["levels"]
	α = system_params["alpha"]

	Faux = @closure t-> (VI(t, tf, χ, levels, α, N)
						 .+ WI(t, tf, χ, coeffs, max_order, levels, α, Ncorr, N))

	coeffs = get_coeffs_aux(Faux, magnus_order, max_order, system_params)
	return coeffs
end
"""
	get_coeffs_aux(Faux, magnus_order, max_order, system_params)
Calculate the correction coefficients.
"""
function get_coeffs_aux(Faux, magnus_order, max_order, system_params)
	tf = system_params["tf"]
	Ncorr = system_params["Ncorr"]
	N = system_params["N"]
	Omega = solveOmega_snap(Faux, tf, 2*magnus_order, N)
	Omega_sum = sum(Omega, dims=1)[1, :, :]

	coeffs = zeros(N*4*max_order + 1)
	for n = 0: Ncorr - 1
		Omega_sum_n = Omega_sum[:, n + 1]
		Mlinear = linear_coeffs(n, max_order, Faux, system_params; eps=1e-10)
		Mnonlinear = nonlinear_coeffs(n, max_order, system_params; eps=1e-10)
		coeffsn = correction_coeffs(Mlinear, Mnonlinear, Omega_sum_n, max_order)

		index = 4*max_order*n + 1: 4*max_order*(n+1)
		coeffs[index] .= coeffsn
	end
	return coeffs
end
"""
	fcorrected(t, tf, α, coeffs)
Calculate corrected envelope function in number state subspace.
"""
function fcorrected(t, tf, α, coeffs)
	gx, gy = g(t, tf, coeffs, max_order)
	theta = pi/2
	T = tf/2
	if t < T
		fx = theta*fc(t, T, 1) + gx
		fy = gy
	else
		fct = theta*fc(t - T, T, 1)
		fx = fct*cos(α) + gx
		fy = fct*sin(α) + gy
	end
	return fx, fy
end
"""
	Fcorrected(t, tf, χ, levels, α, coeffs, Ncorr, N)
Calculate the corrected driving Hamiltonian in the qubit frame.
"""
function Fcorrected(t, tf, χ, levels, α, coeffs, Ncorr, N)
	χef = (χ[1], χ[2]/2, χ[3]/6)
	Ft = zeros(4, N)
	for i = 1: length(levels)
		Fj!(Ft, t, tf, χef, levels[i], α[i], N, false)
	end
	for i = 1: Ncorr
		coeffsj = view(coeffs, 4*max_order*(i - 1) + 1: 4*max_order*i)
		Wj!(Ft, t, tf, χef, i - 1, coeffsj, max_order, N)
	end
	return Ft
end
"""
	correction_pulse_envelope(t, tf, χ, coeffs, max_order, Ncorr)
Calculate the correction pulse envelope in the lab frame.
"""
function correction_pulse_envelope(t, tf, χ, coeffs, max_order, Ncorr)
	χef = (χ[1], χ[2]/2, χ[3]/6)
	gx_lab, gy_lab = zero(t), zero(t)
	for n = 0: Ncorr - 1
		coeffsn = view(coeffs, 4*max_order*n + 1: 4*max_order*(n+1))
		gx, gy = g(t, tf, coeffsn, max_order)
		dφn = χef[1]*n + χef[2]*n*(n - 1) + χef[3]*n*(n - 1)*(n - 2)
		φ = dφn*t
		cosφ, sinφ = cos(φ), sin(φ)
		gx_lab += gx*cosφ + gy*sinφ
		gy_lab += gx*sinφ - gy*cosφ
	end
	return gx_lab, gy_lab
end
"""
	pulse_envelope(t, tf, χ, levels, α)
Calculate the (uncorrected) pulse envelope in the lab frame
"""
function pulse_envelope(t, tf, χ, levels, α)
	χef = (χ[1], χ[2]/2, χ[3]/6)
	fx_lab, fy_lab = zero(t), zero(t)
	for i = 1: length(levels)
		fx, fy = f(t, tf, α[i])
		li = levels[i]
		dφ = χef[1]*li + χef[2]*li*(li - 1) + χef[3]*li*(li - 1)*(li - 2)
		φ = dφ*t
		cosφ, sinφ = cos(φ), sin(φ)
		fx_lab += fx*cosφ + fy*sinφ
		fy_lab += fx*sinφ - fy*cosφ
	end
	return fx_lab, fy_lab
end
"""
	Flab(t, tf, χ, levels, α, N; ideal_interaction::Bool=false)
Calculate the driving Hamiltonian in the lab frame.
"""
function Flab(t, tf, χ, levels, α, N; ideal_interaction::Bool=false)
	χef = (χ[1], χ[2]/2, χ[3]/6)
	Ft = zeros(4, N)
	if ideal_interaction == true
		for i = 1: length(levels)
			indice = levels[i] + 1
			fx, fy = f(t, tf, α[i])
			li = levels[i]
			dφ = χef[1]*li + χef[2]*li*(li - 1) + χef[3]*li*(li - 1)*(li - 2)
			φ = dφ*t
			cosφ, sinφ = cos(φ), sin(φ)
			Ft[1, indice] = fx*cosφ + fy*sinφ
			Ft[2, indice] = fx*sinφ - fy*cosφ
		end
		for i = 1: N
			n = i - 1
			ω = χef[1]*n + χef[2]*n*(n - 1) + χef[3]*n*(n - 1)*(n - 2)
			Ft[3, i] = ω/2
			#Ft[3, i] = χ*(i-1)/2.0
		end
	else
		fx, fy = pulse_envelope(t, tf, χ, levels, α)
		for i = 1: N
			Ft[1, i] = fx
			Ft[2, i] = fy
			n = i - 1
			ω = χef[1]*n + χef[2]*n*(n - 1) + χef[3]*n*(n - 1)*(n - 2)
			Ft[3, i] = ω/2
			#Ft[3, i] = χ*(i-1)/2.0
		end
	end
	return Ft
end
"""
	Flab_corrected(t, tf, χ, levels, α, coeffs, max_order, Ncorr, N)
Calculate the corrected driving Hamiltonian in the lab frame.
"""
function Flab_corrected(t, tf, χ, levels, α, coeffs, max_order, Ncorr, N)
	χef = (χ[1], χ[2]/2, χ[3]/6)
	Ft = zeros(4, N)
	original_pulse = pulse_envelope(t, tf, χ, levels, α)
	correction_pulse = correction_pulse_envelope(t, tf, χ, coeffs, max_order, Ncorr)
	fx, fy =  original_pulse .+ correction_pulse
	for i = 1: N
		Ft[1, i] = fx
		Ft[2, i] = fy
		n = i - 1
		Ft[3, i] = (χef[1]*n + χef[2]*n*(n - 1) + χef[3]*n*(n - 1)*(n - 2))/2
		#Ft[3, i] = χ*(i-1)/2.0
	end
	return Ft
end

max_stars_repo_name: thalesfr/MagnusCorrection

max_issues_repo_name: thalesfr/MagnusCorrection

max_forks_repo_name: thalesfr/MagnusCorrection


text:
"""
we need to give back number of false positive and false negatives and min,max x,y,x of block containing all data
IMPORTANT - in order to avoid bound checking on every iteration we need to keep the dimension of the resulting block be divided by data block cube size for example 32
IMPORTANT - we assume that x dim can not be bigger than 1024
adapted from https://github.com/JuliaGPU/CUDA.jl/blob/afe81794038dddbda49639c8c26469496543d831/src/mapreduce.jl
goldBoolGPU3d - array holding 3 dimensional data of gold standard bollean array
segmBoolGPU3d - array with 3 dimensional the data we want to compare with gold standard
reducedGold - the smallest boolean block (3 dim array) that contains all positive entris from both masks
reducedSegm - the smallest boolean block (3 dim array) that contains all positive entris from both masks
numberToLooFor - number we will analyze whether is the same between two sets
loopNumbYdim - number of times the single lane needs to loop in order to get all needed data - in this kernel it will be exactly a y dimension of a slice
xdim - length in x direction of source array
loopNumbXdim - in case the x dim will be bigger than number of threads we will create second inner loop
cuda arrays holding just single value wit atomically reduced result
,fn,fp
,minxRes,maxxRes
,minyRes,maxyRes
,minZres,maxZres
"""
function getBoolCubeKernel(goldBoolGPU3d
        ,segmBoolGPU3d
        ,reducedGoldA
        ,reducedSegmA
        ,reducedGoldB
        ,reducedSegmB
        ,loopNumbYdim::UInt16
        ,xdim::UInt16
        ,loopNumbXdim::UInt16
        ,numberToLooFor::T
        ,IndexesArray
        ,fn::CuDeviceVector{UInt32, 1}
        ,fp::CuDeviceVector{UInt32, 1}
        ,minxRes::CuDeviceVector{UInt32, 1}
        ,maxxRes::CuDeviceVector{UInt32, 1}
        ,minyRes::CuDeviceVector{UInt32, 1}
        ,maxyRes::CuDeviceVector{UInt32, 1}
        ,minZres::CuDeviceVector{UInt32, 1}
        ,maxZres::CuDeviceVector{UInt32, 1}
        ,warpNumber
) where T

   anyPositive = false # true If any bit will bge positive in this array - we are not afraid of data race as we can set it multiple time to true
#creates shared memory and initializes it to 0
   shmemSum = createAndInitializeShmem(wid,threadIdxX(),lane)
# incrementing appropriate number of times

  #0 - false negative; 1- false positive; 2 -minx; 3 max x; 4 miny; 5 maxy
  locArr= zeros(MVector{6,UInt16})

  @iter3dAdditionalzActs(arrDims,loopXdim,loopYdim,loopZdim,
    #inner expression
    if(  @inbounds($arrAnalyzed[x,y,z])  ==numberToLooFor)
        #updating variables needed to calculate means
        sumX+=Float32(x) ;  sumY+=Float32(y)  ; sumZ+=Float32(z)   ; count+=Float32(1)
    end,
    #after z expression - we get slice wise true counts from it
    begin
        sync_threads()
        #reducing count only
        if(z<=arrDims[3])
            countTemp = count
            @redWitAct(offsetIter,shmemSum, count,+)
            #saving to global memory count of this slice
            @ifXY 1 1 begin
                 $countPerZ[z]=(shmemSum[1,1] - oldZVal[1] )
                oldZVal[1]=shmemSum[1,1]
            end
            #clear shared memory only first row was used and sync threads
            clearSharedMemWarpLong(shmemSum, UInt8(1), Float32(0.0))
            count=countTemp#to preserve proper value for total count
        end#if ar dims
    end )#if bool in arr

   return
   end