text:
# ---
# title: 1367. Linked List in Binary Tree
# id: problem1367
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List, Dynamic Programming, Tree
# link: <https://leetcode.com/problems/linked-list-in-binary-tree/description/>
# hidden: true
# ---
#
# Given a binary tree `root` and a linked list with `head` as the first node.
#
# Return True if all the elements in the linked list starting from the `head`
# correspond to some _downward path_ connected in the binary tree otherwise
# return False.
#
# In this context downward path means a path that starts at some node and goes
# downwards.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/02/12/sample_1_1720.png)**
#
#
#
#     Input: head = [4,2,8], root = [1,4,4,null,2,2,null,1,null,6,8,null,null,null,null,1,3]
#     Output: true
#     Explanation: Nodes in blue form a subpath in the binary Tree.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/02/12/sample_2_1720.png)**
#
#
#
#     Input: head = [1,4,2,6], root = [1,4,4,null,2,2,null,1,null,6,8,null,null,null,null,1,3]
#     Output: true
#
#
# **Example 3:**
#
#
#
#     Input: head = [1,4,2,6,8], root = [1,4,4,null,2,2,null,1,null,6,8,null,null,null,null,1,3]
#     Output: false
#     Explanation: There is no path in the binary tree that contains all the elements of the linked list from head.
#
#
#
#
# **Constraints:**
#
#   * `1 <= node.val <= 100` for each node in the linked list and binary tree.
#   * The given linked list will contain between `1` and `100` nodes.
#   * The given binary tree will contain between `1` and `2500` nodes.
#
#
## @lc code=start
using LeetCode

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


text:
"""
"""
function is_runtests(a::AbstractDict = ENV)
    result = lowercase(strip(get(a, "PREDICTMD_IS_RUNTESTS", ""))) ==
        "true"
    return result
end

"""
"""
function is_make_examples(a::AbstractDict = ENV)
    result = lowercase(strip(get(a, "PREDICTMD_IS_MAKE_EXAMPLES", ""))) ==
        "true"
    return result
end

"""
"""
function is_make_docs(a::AbstractDict = ENV)
    result = lowercase(strip(get(a, "PREDICTMD_IS_MAKE_DOCS", ""))) ==
        "true"
    return result
end

"""
"""
function is_deploy_docs(a::AbstractDict = ENV)
    result = lowercase(strip(get(a, "PREDICTMD_IS_DEPLOY_DOCS", ""))) ==
        "true"
    return result
end

"""
"""
is_docs_or_examples(a::AbstractDict = ENV) =
    is_make_examples(a) || is_make_docs(a) || is_deploy_docs(a)

ext: jl


text:
# psf-gauss.jl

using SPECTrecon: psf_gauss
using Test: @test, @testset, @test_throws, @inferred

@testset "psf" begin
    psf = @inferred psf_gauss()
    @test psf isa Array{Float32,3}

    ny = 4
    px = 5
    pz = 3
    psf = @inferred psf_gauss(; ny, px, pz, fwhm = zeros(ny))
    tmp = zeros(px,pz)
    tmp[(end+1)÷2,(end+1)÷2] = 1 # Kronecker impulse
    tmp = repeat(tmp, 1, 1, ny)
    @test psf == tmp

    ny = 4
    px = 5
    pz = 3
    psf = @inferred psf_gauss(; ny, px, pz,
        fwhm_x = fill(Inf, ny),
        fwhm_z = zeros(ny),
    )
    tmp = zeros(px,pz)
    tmp[:,(end+1)÷2] .= 1/px # wide in x
    tmp = repeat(tmp, 1, 1, ny)
    @test psf ≈ tmp
end


text:
Base.@kwdef struct PardisoJL <: SciMLLinearSolveAlgorithm
    nprocs::Union{Int, Nothing} = nothing
    solver_type::Union{Int, Pardiso.Solver, Nothing} = nothing
    matrix_type::Union{Int, Pardiso.MatrixType, Nothing} = nothing
    iparm::Union{Vector{Tuple{Int,Int}}, Nothing} = nothing
    dparm::Union{Vector{Tuple{Int,Int}}, Nothing} = nothing
end

MKLPardisoFactorize(;kwargs...) = PardisoJL(;solver_type = 0,kwargs...)
MKLPardisoIterate(;kwargs...) = PardisoJL(;solver_type = 1,kwargs...)
needs_concrete_A(alg::PardisoJL) = true

# TODO schur complement functionality

function init_cacheval(alg::PardisoJL, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
    @unpack nprocs, solver_type, matrix_type, iparm, dparm = alg
    A = convert(AbstractMatrix,A)

    solver =
    if Pardiso.PARDISO_LOADED[]
        solver = Pardiso.PardisoSolver()
        solver_type !== nothing && Pardiso.set_solver!(solver, solver_type)

        solver
    else
        solver = Pardiso.MKLPardisoSolver()
        nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs)

        solver
    end

    Pardiso.pardisoinit(solver) # default initialization

    if matrix_type !== nothing
        Pardiso.set_matrixtype!(solver, matrix_type)
    else
        if eltype(A) <: Real
            Pardiso.set_matrixtype!(solver, Pardiso.REAL_NONSYM)
        elseif eltype(A) <: Complex
            Pardiso.set_matrixtype!(solver, Pardiso.COMPLEX_NONSYM)
        else
            error("Number type not supported by Pardiso")
        end
    end
    verbose && Pardiso.set_msglvl!(solver, Pardiso.MESSAGE_LEVEL_ON)

    # pass in vector of tuples like [(iparm::Int, key::Int) ...]
    if iparm !== nothing
        for i in length(iparm)
            Pardiso.set_iparm!(solver, iparm[i]...)
        end
    end

    if dparm !== nothing
        for i in length(dparm)
            Pardiso.set_dparm!(solver, dparm[i]...)
        end
    end

    # Make sure to say it's transposed because its CSC not CSR
    Pardiso.set_iparm!(solver,12, 1)

    #=
    Note: It is recommended to use IPARM(11)=1 (scaling) and IPARM(13)=1 (matchings) for
    highly indefinite symmetric matrices e.g. from interior point optimizations or saddle point problems.
    It is also very important to note that the user must provide in the analysis phase (PHASE=11)
    the numerical values of the matrix A if IPARM(11)=1 (scaling) or PARM(13)=1 or 2 (matchings).

    The numerical values will be incorrect since the analysis is ran once and
    cached. If these two are not set, then Pardiso.NUM_FACT in the solve must
    be changed to Pardiso.ANALYSIS_NUM_FACT in the solver loop otherwise instabilities
    occur in the example https://github.com/SciML/OrdinaryDiffEq.jl/issues/1569
    =#
    Pardiso.set_iparm!(solver,11, 0)
    Pardiso.set_iparm!(solver,13, 0)

    Pardiso.set_phase!(solver, Pardiso.ANALYSIS)

    if alg.solver_type == 1
        # PARDISO uses a numerical factorization A = LU for the first system and
        # applies these exact factors L and U for the next steps in a
        # preconditioned Krylov-Subspace iteration. If the iteration does not
        # converge, the solver will automatically switch back to the numerical factorization.
        Pardiso.set_iparm!(solver,3,round(Int,abs(log10(reltol)),RoundDown) * 10 + 1)
    end

    Pardiso.pardiso(solver, u, A, b)

    return solver
end

function SciMLBase.solve(cache::LinearCache, alg::PardisoJL; kwargs...)
    @unpack A, b, u = cache
    A = convert(AbstractMatrix,A)

    if cache.isfresh
        Pardiso.set_phase!(cache.cacheval, Pardiso.NUM_FACT)
        Pardiso.pardiso(cache.cacheval, A, eltype(A)[])
    end

    Pardiso.set_phase!(cache.cacheval, Pardiso.SOLVE_ITERATIVE_REFINE)
    Pardiso.pardiso(cache.cacheval, u, A, b)
    return SciMLBase.build_linear_solution(alg,cache.u,nothing,cache)
end

# Add finalizer to release memory
# Pardiso.set_phase!(cache.cacheval, Pardiso.RELEASE_ALL)

export PardisoJL, MKLPardisoFactorize, MKLPardisoIterate


text:
__precompile__(true)

module DatasToolbox

using Reexport

@reexport using DataUtils
@reexport using DataHandlers
# for now we won't re-export the Python stuff

end

max_stars_repo_name: JuliaTagBot/DatasToolbox

max_issues_repo_name: JuliaTagBot/DatasToolbox

max_forks_repo_name: JuliaTagBot/DatasToolbox


text:
function EuropeanCallMC_DeltaBased(K, T, S, sig, r, div, N, M)

## setting parameters

dt = T/N
nudt = (r - div - 0.5 * sig^2) * dt
sigsdt = sig * sqrt(dt)
erddt = exp((r - div) * dt)

beta1 = -1

## variables to store the sums of the option values

sum_CT = 0
sum_CT2 = 0

function CND(X)
#   y::float64
#   Exponential::float64
#   SumA::float64
#   SumB::float64

  y = abs(X)

  if y > 37
    cnd = 0
  else
    Exponential = exp(-y^2/2)

    if y < 7.07106781186547
      SumA = 0.0352624965998911 * y + 0.700383064443688
      SumA = SumA * y + 6.37396220353165
      SumA = SumA * y + 33.912866078383
      SumA = SumA * y + 112.079291497871
      SumA = SumA * y + 221.213596169931
      SumA = SumA * y + 220.206867912376
      SumB = 0.0883883476483184 * y + 1.75566716318264
      SumB = SumB * y + 16.064177579207
      SumB = SumB * y + 86.7807322029461
      SumB = SumB * y + 296.564248779674
      SumB = SumB * y + 637.333633378831
      SumB = SumB * y + 793.826512519948
      SumB = SumB * y + 440.413735824752
      cnd = Exponential * SumA / SumB
    else
      SumA = y + 0.65
      SumA = y + 4 / SumA
      SumA = y + 3 / SumA
      SumA = y + 2 / SumA
      SumA = y + 1 / SumA
      cnd = Exponential / (SumA * 2.506628274631)

      end
  end

  if X > 0
    cnd = 1 - cnd
  end

cnd

end

function GDelta(CallPutFlag, S, X, T, r, b, v)

  #require("CND.jl")
  d1::Float64

  d1 = (log(S / X) + (b + v^2 / 2) * T) / (v * sqrt(T))

  if CallPutFlag == "c"
    delta = exp((b - r) * T) * CND(d1)

  elseif CallPutFlag == "p"
    delta = -exp((b - r) * T) * CND(-d1)
  end

    delta

end

## outer loop is the number of simulations bigger is better

for j = 1:M

  St = S
    cv = 0

    ## inner loop is simulating each price path of length N
  for i = 1:N

    t = (i-1) * dt
    delta = GDelta("c", St, K, t, r, div, sig)

    error = randn()
    Stn = St * exp(nudt + sigsdt * error) ## note, it only add the changes, does not store the path
    cv = cv + delta * (Stn-St*erddt)
    St = Stn

  end

    ## wraps up and stores all the values of each path

  CT = max(0, St - K) + beta1*cv
  sum_CT = sum_CT + CT
  sum_CT2 = sum_CT2 + CT^2

end

  ## discounting the average call value
call_value = sum_CT / M * exp(-r*T)

SD = (sqrt(sum_CT2 - 1/M * (sum_CT^2)) * exp(-2*r*T)) / (M-1)
SE = SD / sqrt(M)

println("The Call Value is $call_value")
println("The SD is $SD")
println("The SE is $SE")

MC = [call_value, SD, SE]

end


text:
module Xfoil

using xfoil_light_jll, Printf

export set_coordinates, pane, solve_alpha, get_xsep, bldump, alpha_sweep
export set_coordinates_cs, pane_cs, solve_alpha_cs, get_xsep_cs, bldump_cs, alpha_sweep_cs

# Fixed array limitations from Fortran shared libraries
const IBX=572
const IPX=5
const IQX=286
const ISX=2
const IZX=322
const IVX=229
const NPX=8
const IWK=36

function __init__()
    # wrap XFOIL global variables
    global xfoilglobals = get_globals()
    global xfoilglobals_cs = get_globals_cs()
    return nothing
end

# wrappers for XFOIL global variables
include("globals.jl")

# wrappers for XFOIL functions
include("libxfoil.jl")

# higher level functions
include("higherlevel.jl")

end #module


text:
# Paper notes
# Channel is emptyed with every setup to prevent unexpected behavior.

using Distributed

using Dates

# const BlockAndMetaOrNothing = Union{Tuple{Array{Int8,3},Array{Int64,1},Array{Int64,1},Int64}, Nothing}
const BlockAndMetaOrNothing = Union{Int64, Nothing}
const ArrayOrNothing = Union{Array, Nothing}
_single_boundary3 = nothing
# _b3_size = nothing # it is set in boundary_operator.jl
# _ch_block::RemoteChannel = nothing
# _ch_results::RemoteChannel = nothing
# _ch_block = nothing
# _ch_results = nothing
_ch_block = RemoteChannel(()->Channel{BlockAndMetaOrNothing}(32));
_ch_results = RemoteChannel(()->Channel{ArrayOrNothing}(32));
_data_size = nothing
_workers_running = false
_reference_time = nothing
# time data cane be used for measurement
# set it with _set_time_data()
# init on all workers is done in lsp_setup()
# data are computed in lsp_do_work_code_multiply_decode()
_time_data = nothing
_time_data_i = 1
_segmentation_data = nothing
_bgetter = nothing

# const _local_time_data = Dict()
# _setup_counter = 0

function set_single_boundary3(b3, block_size)
    global _single_boundary3, _b3_size
    _single_boundary3 = b3
    _b3_size = block_size
    # println("set boundary ")
end

function set_data_size(data_size)
    @debug "Set data size"
    global _data_size
    _data_size = data_size
end

function set_channels(ch_block, ch_results)
    global _ch_block, _ch_results
    _ch_block = ch_block
    _ch_results = ch_results
end

function set_time_data(tm_data=nothing)
    global _time_data, _time_data_i
    if tm_data == nothing
        tm_data = Dict()
    end
    _time_data = tm_data
    if _time_data != nothing
        # _time_data["total take time $(myid())"] = 0.
        # _time_data["total computation time $(myid())"] = 0.
        # _time_data["total put time $(myid())"] = 0.
        _time_data["cumulative take time $(myid())"] = 0.
        _time_data["cumulative computation time $(myid())"] = 0.
        _time_data["cumulative put time $(myid())"] = 0.
        _time_data["number of processed blocks $(myid())"] = 0
        _time_data_i = 1
    end
end

function _get_local_time_data()
    return _time_data
end

function get_time_data()
    global _time_data
    if _time_data != nothing
        # @info "_time_data before" _time_data
        for wid in workers()
            tmd = @spawnat wid LarSurf._get_local_time_data()
            tmdf = fetch(tmd)
            # @info "fetch tmd" tmdf
            merge!(_time_data, tmdf)

        end
        # @info "_time_data after" _time_data
    end
    return _time_data
end

"""
Lar Surf Parallel setup.
The reference time can be passed in to measure time.
"""
function lsp_setup(block_size; reference_time=nothing)
    global _workers_running, _reference_time
    global _ch_block, _ch_results
    global _b3_size
    global _single_boundary3
    if reference_time == nothing
        _reference_time = time()
    else
        _reference_time = reference_time
    end
    if nprocs() == 1
        addprocs(1)
        @info "Adding one worker"
    end
    # @everywhere using LarSurf
    # global _ch_block, _ch_results
    @info "lsp setup with block_size: $block_size"
    block_size = LarSurf.size_as_array(block_size)
    # println("block_size: $block_size")
    # Causes error = non existing SparseMatrixCSC
    LarSurf.set_param(boundary_allow_read_files=true, boundary_allow_write_files=true)
    # LarSurf.set_param(boundary_allow_read_files=false, boundary_allow_write_files=false)
    # tm = @elapsed b3, larmodel = LarSurf.get_boundary3(block_size, false)
    tm = @elapsed b3 = LarSurf.get_boundary3(block_size, false)
    @info "time for construction of b3 " tm
    # println("block_size: $block_size")
    # set on local
    # set_single_boundary3(b3, block_size)
    _b3_size = block_size
    _single_boundary3 = b3
    @debug "b3 calculated, _b3_size: $_b3_size"
    tim0 = time()
    @sync for wid in workers()
        @info "starting worker id: $wid"
        # ftch[wid] =
        @spawnat wid LarSurf.set_single_boundary3(b3, block_size)

        # @TODO this can be in next for and without sync
    end
    @info "time used for starting workers$(time() - tim0)"
    if _workers_running
        empty_channel(_ch_block)
        empty_channel(_ch_results)
        @debug "Workers are alredy initiated. No need to do it again."
    else
        for wid in workers()
                @debug "Initiating workers"
                @info "remote do for worker $(myid())"
                remote_do(
                    LarSurf.lsp_do_work_code_multiply_decode,
                    wid,
                    _ch_block,
                    _ch_results,
                    # LarSurf._ch_block,
                    # LarSurf._ch_results,
                )
        end
        _workers_running = true
    end
    @info "setup done in $(time() - _reference_time)"

end

# function lsp_deinit_workers(ch_block::RemoteChannel=nothing, ch_faces::RemoteChannel=nothing)

function lsp_deinit_workers()
    global _workers_running, _ch_block
    if _workers_running
        @debug "Sending 'nothing'"
        for wid in workers()
            put!(_ch_block, nothing)
        end
        _workers_running = false
    else
        @warn "No workers to deinit"
    end
end

"""
Extract surface after lsp_setup().

If set_time_data is set to Dict, the time measurements will be stored there.
"""
function lsp_get_surface(segmentation; voxelsize=[1,1,1])

    n, bgetter, data_size = LarSurf.lsp_setup_data(segmentation)
    for wid in workers()
        if _time_data != nothing
            @spawnat wid LarSurf.set_time_data(_time_data)
        end
    end

    @async LarSurf.lsp_job_enquing(n, bgetter)

    # Workers are started before in setup

    # print("===== Output Faces =====")
    @info "============== Sequential finished ======== in time $(time()-_reference_time) [s]"
    # cumulative_decoding_time = 0
    cumulative_doubled_filtration_time = 0
    tm_end = @elapsed begin
        numF = LarSurf.grid_number_of_faces(data_size)
        bigFchar = spzeros(Int8, numF)
        for i=1:n

            @debug "Collecting the information for block $i (first three messages)" maxlog=3
            faces_per_block = take!(_ch_results)
            if i == 1
                @info "First faces recived in time: $(time()-_reference_time) [s]"
                if _time_data != nothing
                    _time_data["first face recived"] = time()-_reference_time
                end
            end
            # TODO does it make sense to make an optimization? Copy of the faces into sparsearray have to be done anyway.
            # println(faces_per_block)
            # for  block_i=1:block_number
            tm_doubled_filtration = @elapsed for big_fid in faces_per_block
                # median time 31 ns
                # bigFchar[big_fid] = (bigFchar[big_fid] + 1) % 2
                # median time 14 ns
                bigFchar[big_fid] = if (bigFchar[big_fid] == 1) 0 else 1 end
            end
            cumulative_doubled_filtration_time = cumulative_doubled_filtration_time + tm_doubled_filtration
            # println(bigFchar)
            # end

        end
        @info "=== Parallel processing finished === Last faces recived in time: $(time()-_reference_time) [s]"
        # @info "Last faces recived in time: $(time()-_reference_time) [s] === end sequential ==="
        if _time_data != nothing
            _time_data["last face recived"] = time()-_reference_time
        end
        dropzeros!(bigFchar)
        # bigV, FVreduced = grid_Fchar_to_Vreduced_FVreduced(bigFchar, data_size; voxelsize=voxelsize)
        reduction_time = @elapsed bigV, FVreduced = grid_Fchar_to_Vreduced_FVreduced(bigFchar, data_size; voxelsize=voxelsize)
        # cumulative_decoding_time = cumulative_decoding_time + dec_time
    end
    @info "reciving and decoding time [s]: $tm_end"
    if _time_data != nothing
        _time_data["reciving and decoding time [s]"] = tm_end
        _time_data["V and FV reduction time [s]"] = reduction_time
        _time_data["cumulative double face filtration time [s]"] = cumulative_doubled_filtration_time
    end
    return bigV, FVreduced
    # return __make_return(V, filteredFV, Flin, larmodel, return_all)
    # return bigFchar
end

function set_data_local(data)
    global _segmentation_data = data
end

function set_bgetter(bgetter)
    global _bgetter = bgetter
end

function lsp_setup_data(segmentation)
    # global _b3_size
    # println("b3_size type $(typeof(_b3_size))")
    data_size = LarSurf.size_as_array(size(segmentation))
    n, bgetter = LarSurf.block_getter(segmentation, _b3_size, fixed_block_size=true)
    @sync for wid in workers()
        @info "set data size on worker: $wid"
        # ftch[wid] =
        @spawnat wid LarSurf.set_data_size(data_size)
        @spawnat wid LarSurf.set_bgetter(bgetter)
    end
    return n, bgetter, data_size

end

function lsp_job_enquing(n, bgetter)
    global _ch_block
    # global _b3_size
    # println("b3_size type $(typeof(_b3_size))")
    # n, bgetter = LarSurf.block_getter(segmentation, _b3_size, fixed_block_size=true)

    for block_id=1:n
        # @info "getting block $block_id"
        # tm_get_block = @elapsed block = LarSurf.get_block(block_id, bgetter...)
        # @info "putting block $block_id to channel"

        tm_put = @elapsed put!(_ch_block, block_id)
        # tm_put = @elapsed put!(_ch_block, (block..., block_id))
        if block_id == 1
            @info "First block sent in channel in time: $(time()-_reference_time) [s]"
            if _time_data != nothing
               _time_data["first block sent"] = time()-_reference_time
            end
        end

        # @info "== Job enquing get_block time: $tm_get_block, put time: $tm_put"
        # , channel size: $(lenght(_ch_block))"
    end
    @info "Last block sent in channel in time: $(time()-_reference_time) [s]"
    @info "Job enquing finished"
end

function lsp_do_work_code_multiply_decode(ch_block, ch_faces)
    @info "do work while initiated on worker $(myid())"
    global _data_size, _time_data_i
    # ref_time = 0.
    # after_take_time = 0
    # _local_time_data["first waiting time"] = 0.
    # _time_data[]
    while true
        @info "waiting for data on worker $(myid())" maxlog=2
        ref_time = time()

        fbl = take!(ch_block)
        after_take_time = time()
        # println("type of : $(typeof(ch_block))")
        if fbl == nothing
            @info "recived 'nothing' quit worker $(myid())"
            # put!(ch_block, nothing)
            return
        end
        # @info "code mul decode on block $(fbl[end]), worker $(myid()), " maxlog=3

         # $(time()-_reference_time)(first 3 messages per worker)" maxlog=3
        block_id = fbl
        data_size = _data_size
        block = LarSurf.get_block(block_id, _bgetter...)
        faces = LarSurf.code_multiply_decode(data_size, block..., block_id)
        # faces = LarSurf.code_multiply_decode(data_size, fbl...)
        after_work_time = time()
        put!(ch_faces, faces)
        after_put_time = time()
        # time measurement
        ttake = after_take_time - ref_time
        twork = after_work_time - after_take_time
        tput = after_put_time - after_work_time
        if _time_data != nothing
            _time_data["cumulative take time $(myid())"] += ttake
            _time_data["cumulative computation time $(myid())"] += twork
            _time_data["cumulative put time $(myid())"] += tput
            _time_data["number of processed blocks $(myid())"] += 1
            if _time_data_i == 1
                _time_data["first take time $(myid())"] = ttake
            end
            _time_data_i = 2
        end
        @debug "code mul decode on block $(fbl[end]), worker $(myid()), $ttake, $twork, $tput"
    end
end

function three_chain_coding(segmentation)
    return grid_to_linear(segmentation, 0)
end

"""
Decode 2-chain (Flin) to global IDs.
"""
function two_chain_decoding(Flin, data_size, offset, block_size, block_id)
    # TODO block_size and offset are in different order in two functions
    # sub_bgrid_fac..face_id() and  get_block()
# face from small to big

    # println("block_id: $block_id, data_size: $data_size, block_size: $block_size, offset: $offset")
    i, j, v = findnz(Flin)
    # fid_small = [
    #     fid
    #     for fid in j
    # ]
    ret = [
        sub_grid_face_id_to_orig_grid_face_id(
        data_size, block_size, offset, fid
        )[1]
        for fid in j
    ]
    # println("decode")
    # println(fid_small)
    # println(ret)
    return ret
end

"""
Calculate multiplication linearized volume cells with boundary matrix and
filter it for number 1.
"""
# function grid_get_surf_Fvec_larmodel(segmentation::AbstractArray)
function code_multiply_decode(
    data_size::Array,
    segmentation::Array, offset, block_size, block_id
    )
    # @debug "code multiply decode"
    # println("code multiply decode")

    segClin = three_chain_coding(segmentation)
    # println("segClin: $segClin")
    b3 = _single_boundary3

    # Multiplication step
    Flin = segClin' * b3
    # println("Flin: $Flin")

    ## Filtration
    sparse_filter!(Flin, 1, 1, 0)
    dropzeros!(Flin)
    # println("Flin filtered: $Flin")
    # println("_b3_size: $(typeof(_b3_size))")
    # _b3_size

    Flist = two_chain_decoding(Flin,
        # data_size, offset, block_size, block_id
        data_size, offset, _b3_size, block_id
    )
    return Flist
end

function empty_channel(ch::RemoteChannel)
    while isready(ch)
        take!(ch)
        @warn "Something was in channel. Now it is empty." maxlog=1
    end
end


text:
using Flux

l1 = Dense(32, 64)
l2 = Dense(64, 128)
l3 = Dense(128, 64)

topo = @functopo x => a => b => y

@testset "TopoChain" begin
    model = TopoChain(topo, l1, l2, l3)

    @test model[2] == l2
    @test length(model) == 3
end


text:
"""
Convert ASCII (RFC20) character `c` to lower case.
"""
ascii_lc(c::UInt8) = c in UInt8('A'):UInt8('Z') ? c + 0x20 : c

"""
Case insensitive ASCII character comparison.
"""
ascii_lc_isequal(a::UInt8, b::UInt8) = ascii_lc(a) == ascii_lc(b)

"""
Case insensitive ASCII string comparison.
"""
function ascii_lc_isequal(a, b)
    acu = codeunits(a)
    bcu = codeunits(b)
    len = length(acu)
    len != length(bcu) && return false
    for i = 1:len
        @inbounds !ascii_lc_isequal(acu[i], bcu[i]) && return false
    end
    return true
end


text:
using .SFEM
using BenchmarkTools

A = "G://SMESH//4.cen"
B = "G://SMESH//4.ele"
C = "G://SMESH//4.facein"
D = "G://SMESH//4.faceout"
E = "G://SMESH//4.node"

const EE, NU, FF = 3E7, 0.3, -1000.0
Node, Ele, Cen, Face_out, Face_in = IOsfem(A, B, C, D, E)
@time begin
    KE, Q = Ksfem(Node, Ele, Cen, Face_out, Face_in, EE, NU, FF)
end
#=
@time begin
    Disp = MKL(KE, Q)
end
=#

max_stars_repo_path: src/multi-core juSFEM/main.jl

max_stars_repo_name: Elliothuo/juSFEM

max_issues_repo_path: src/multi-core juSFEM/main.jl

max_issues_repo_name: Elliothuo/juSFEM

max_forks_repo_path: src/multi-core juSFEM/main.jl

max_forks_repo_name: Elliothuo/juSFEM


text:
database_kind(c::Type{Redis.RedisConnection}) = NonRelational

redis_id(::Type{T}, id) where T<:Model = "$T:$id"
redis_wildcard(::Type{T}) where T<:Model = "$T:*"

function insert!(mapper::DBMapper, ::Type{Redis.RedisConnection}, elem::T)  where T<:Model
    table = mapper.tables[T]
    conn = get_connection(mapper.pool)
    if table.primary_key.has_auto_value
        setid!(elem, mapper, generate_id(elem, idtype(T, mapper)))
    end
    id = redis_id(T, getid(elem, mapper))
    result = Redis.hmset(conn, id, marshal(mapper, elem))
    release_connection(mapper.pool, conn)
end

function select_one(mapper::DBMapper, ::Type{Redis.RedisConnection}, T::Type{<:Model}; kwargs...)
    params = Dict(kwargs...)
    id_field = idfield(mapper, T)
    id = params[id_field]
    conn = get_connection(mapper.pool)
    result = Redis.hgetall(conn, redis_id(T, id))
    release_connection(mapper.pool, conn)
    if isempty(result)
        return nothing
    else
        return T(;unmarshal(mapper, T, result)...)
    end
end

function clean_table!(mapper::DBMapper, ::Type{Redis.RedisConnection}, elem::Type{T})  where T<:Model
    conn = get_connection(mapper.pool)
    (cursor, results) = Redis.scan(conn, 0, "match", redis_wildcard(T))
    Redis.unlik(conn, keys...)
    while cursor != 0
        (cursor, results) = Redis.scan(conn, cursor, "match", redis_wildcard(T))
        Redis.unlik(conn, keys...)
    end
    release_connection(mapper.pool, conn)
end

function drop_table!(::DBMapper, ::Type{Redis.RedisConnection}, elem::Type{T})  where T<:Model

end
function update!(mapper::DBMapper, ::Type{Redis.RedisConnection}, elem::T; fields::Array{Symbol}=Symbol[])  where T<:Model
    id = getid(elem, mapper)
    if id === nothing
        throw("The element does not have id. Cannot update")
    end
    data = marshal(mapper, elem)
    if length(fields) > 0
        for f in keys(data)
            if !(f in fields)
                pop!(data, f)
            end
        end
    end
    id = redis_id(T, id)
    conn = get_connection(mapper.pool)
    result = Redis.hmset(conn, id, data)
    release_connection(mapper.pool, conn)
end

function extract_id_from_kwargs(mapper::DBMapper, T::Type{<:Model}; kwargs...)
    id_f = idfield(mapper, T)
    params = Dict(kwargs...)
    id = nothing
    if haskey(params, id_f)
        id = redis_id(T, pop!(params, id_f))
    end
    return (id, params)
end

function select_by_key_and_params(mapper::DBMapper, dbtype::Type{Redis.RedisConnection}, T::Type{<:Model}, pk, params)
    conn = get_connection(mapper.pool)
    params_key = collect(keys(params))
    elem_values = Redis.hmget(conn, pk, params_key...)
    release_connection(mapper.pool, conn)
    return unmarshal(mapper, T, Dict(zip(params_key, elem_values)); partial=true)
end

function iterate_all(mapper::DBMapper,  dbtype::Type{Redis.RedisConnection}, T::Type{<:Model},
                    conditions, apply_function, fields::Array{Symbol}=Symbol[])
    conn = get_connection(mapper.pool)
    cursor = -1
    conditions_fields = collect(keys(conditions))
    if isempty(fields)
        fields = column_names(mapper, T)
    end
    fields_to_obtain = union(conditions_fields, fields)
    while cursor != 0
        (cursor, results) = Redis.scan(conn, cursor == -1 ? 0 : cursor, "match", redis_wildcard(T))
        for elem_key in results
            elem_values = Redis.hmget(conn, elem_key, fields_to_obtain...)
            elem = unmarshal(mapper, T, Dict{Symbol, Any}(zip(fields_to_obtain, elem_values)); partial=true)
            ret = apply_function(elem_key, elem, conditions, conditions_fields)
            if (typeof(ret) <: Bool && ret == true) || (!(typeof(ret) <: Bool) && ret !== nothing)
                release_connection(mapper.pool, conn)
                return ret
            end
        end
    end
    release_connection(mapper.pool, conn)
end

function exists(mapper::DBMapper, dbtype::Type{Redis.RedisConnection}, T::Type{<:Model}; kwargs...)
    (id, params) = extract_id_from_kwargs(mapper, T; kwargs...)
    # Query with id
    if id !== nothing
        # Only id
        if length(params) == 0
            conn = get_connection(mapper.pool)
            result = Redis.exists(conn, id)
            release_connection(mapper.pool, conn)
            return result
        else
            # ID & fields
            elem = select_by_key_and_params(mapper, dbtype, T, id, params)
            return all([params[k] == elem[k] for k in keys(params)])
        end
    else
        found = (elem_id, elem, params, params_key)->all([params[k] == elem[k] for k in params_key])
        ret = iterate_all(mapper, dbtype, T, params, found)
        if ret !== nothing
            return ret
        end
    end
    return false
end
import Base:delete!
function delete!(mapper::DBMapper, dbtype::Type{Redis.RedisConnection}, T::Type{<:Model}; kwargs...)
    (id, params) = extract_id_from_kwargs(mapper, T; kwargs...)
    if id !== nothing
        if length(params) == 0
            conn = get_connection(mapper.pool)
            result = Redis.del(conn, id)
            release_connection(mapper.pool, conn)
        else
            elem = select_by_key_and_params(mapper, dbtype, T, id, params)
            if all([params[k] == elem[k] for k in keys(params)])
                Redis.del(conn, id)
            end
        end
    else
        all_keys_to_delete = []
        add_keys = (elem_key, elem, params, params_key)->begin
            if all([params[k] == elem[k] for k in params_key])
                push!(all_keys_to_delete, elem_key)
            end
            return nothing
        end
        ret = iterate_all(mapper, dbtype, T, params, add_keys)
        conn = get_connection(mapper.pool)
        for elem_key in all_keys_to_delete
            Redis.del(conn, elem_key)
        end
        release_connection(mapper.pool, conn)
    end
end

function select_all(mapper::DBMapper, dbtype::Type{Redis.RedisConnection}, T::Type{<:Model}; fields::Array{Symbol}=[], kwargs...) ::Array{T}
    conditions = Dict(kwargs...)
    all_elems = []
    found_function = (elem_key, elem, params, params_key)->begin
        if all([params[k] == elem[k] for k in params_key])
            push!(all_elems, elem)
        end
        return nothing
    end
    ret = iterate_all(mapper, dbtype, T, conditions, found_function, fields)
    if isempty(all_elems)
        return nothing
    else
        return [T(;unmarshal(mapper, T, elem, partial=!isempty(fields))...) for elem in all_elems]
    end
end

max_stars_repo_name: lucianolorenti/StructDatabaseMapping

max_issues_repo_name: lucianolorenti/StructDatabaseMapping

max_forks_repo_name: lucianolorenti/StructDatabaseMapping


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

import Base: values, convert
using Base.Threads
import Base.Threads: @threads, nthreads, threadid
using JuMP
using MathOptInterface
const MOI = MathOptInterface

function init(constructor)::MIP
    return MIP(constructor, Any[nothing for t = 1:nthreads()], Variable[], 1.0, 0)
end

function read!(mip::MIP, filename::AbstractString)::Nothing
    @threads for t = 1:nthreads()
        model = read_from_file(filename)
        mip.optimizers[t] = backend(model)
        _replace_zero_one!(mip.optimizers[t])
        if t == 1
            _assert_supported(mip.optimizers[t])
            mip.binary_variables = _get_binary_variables(mip.optimizers[t])
            mip.sense = _get_objective_sense(mip.optimizers[t])
        end
        _relax_integrality!(mip.optimizers[t])
        set_optimizer(model, mip.constructor)
        set_silent(model)
    end
    return
end

function _assert_supported(optimizer::MOI.AbstractOptimizer)::Nothing
    types = MOI.get(optimizer, MOI.ListOfConstraints())
    for (F, S) in types
        _assert_supported(F, S)
    end
end

function _assert_supported(F::DataType, S::DataType)::Nothing
    if F in [MOI.ScalarAffineFunction{Float64}, MOI.SingleVariable] && S in [
        MOI.LessThan{Float64},
        MOI.GreaterThan{Float64},
        MOI.EqualTo{Float64},
        MOI.Interval{Float64},
    ]
        return
    end
    if F in [MOI.SingleVariable] && S in [MOI.Integer, MOI.ZeroOne]
        return
    end
    error("MOI constraint not supported: $F in $S")
end

function _get_objective_sense(optimizer::MOI.AbstractOptimizer)::Float64
    sense = MOI.get(optimizer, MOI.ObjectiveSense())
    if sense == MOI.MIN_SENSE
        return 1.0
    elseif sense == MOI.MAX_SENSE
        return -1.0
    else
        error("objective sense not supported: $sense")
    end
end

_bounds_constraint(v::Variable) =
    MOI.ConstraintIndex{MOI.SingleVariable,MOI.Interval{Float64}}(v.index)

function _replace_zero_one!(optimizer::MOI.AbstractOptimizer)::Nothing
    constrs_to_delete = MOI.ConstraintIndex[]
    funcs = MOI.SingleVariable[]
    sets = Union{MOI.Interval,MOI.Integer}[]
    for ci in
        MOI.get(optimizer, MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.ZeroOne}())
        func = MOI.get(optimizer, MOI.ConstraintFunction(), ci)
        var = func.variable
        push!(constrs_to_delete, ci)
        push!(funcs, MOI.SingleVariable(var))
        push!(funcs, MOI.SingleVariable(var))
        push!(sets, MOI.Interval{Float64}(0.0, 1.0))
        push!(sets, MOI.Integer())
    end
    MOI.delete(optimizer, constrs_to_delete)
    MOI.add_constraints(optimizer, funcs, sets)
    return
end

function _get_binary_variables(optimizer::MOI.AbstractOptimizer)::Vector{Variable}
    vars = Variable[]
    for ci in
        MOI.get(optimizer, MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.Integer}())
        func = MOI.get(optimizer, MOI.ConstraintFunction(), ci)
        var = Variable(func.variable.value)

        MOI.is_valid(optimizer, _bounds_constraint(var)) ||
            error("$var is not interval-constrained")
        interval = MOI.get(optimizer, MOI.ConstraintSet(), _bounds_constraint(var))
        interval.lower == 0.0 || error("$var has lb != 0")
        interval.upper == 1.0 || error("$var has ub != 1")

        push!(vars, var)
    end
    return vars
end

function _relax_integrality!(optimizer::MOI.AbstractOptimizer)::Nothing
    indices =
        MOI.get(optimizer, MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.Integer}())
    MOI.delete(optimizer, indices)
end

"""
    solve_relaxation(mip::MIP)::Tuple{Symbol, Float64}

Solve the linear relaxation of `mip` and returns a tuple containing the
solution status (either :Optimal or :Infeasible) and the optimal objective
value. If the problem is infeasible, the optimal value is Inf for minimization
problems and -Inf for maximization problems..
"""
function solve_relaxation!(mip::MIP)::Tuple{Symbol,Float64}
    t = threadid()
    MOI.optimize!(mip.optimizers[t])
    try
        mip.lp_iterations += MOI.get(mip.optimizers[t], MOI.SimplexIterations())
    catch
        # ignore
    end
    status = MOI.get(mip.optimizers[t], MOI.TerminationStatus())
    if status == MOI.OPTIMAL
        obj = MOI.get(mip.optimizers[t], MOI.ObjectiveValue())
        return :Optimal, obj * mip.sense
    elseif status in [MOI.INFEASIBLE, MOI.INFEASIBLE_OR_UNBOUNDED]
        return :Infeasible, Inf * mip.sense
    end
    error("unknown status: $status")
end

"""
    values(mip::MIP, vars::Vector{Variable})::Array{Float64}

Returns a vector `vals` which describes the current primal values for the
decision variables `vars`. More specifically, `vals[j]` is the current
primal value of `vars[j]`.
"""
function values(mip::MIP, vars::Vector{Variable})::Array{Float64}
    return MOI.get(
        mip.optimizers[threadid()],
        MOI.VariablePrimal(),
        convert.(MOI.VariableIndex, vars),
    )
end

values(mip::MIP) =
    values(mip, MOI.get(mip.optimizers[threadid()], MOI.ListOfVariableIndices()))

"""
    set_bounds!(mip::MIP,
                vars::Vector{Variable},
                lb::Array{Float64},
                ub::Array{Float64})

Modify the bounds of the given variables. More specifically, sets
upper and lower bounds of `vars[j]` to `ub[j]` and `lb[j]`, respectively.
"""
function set_bounds!(
    mip::MIP,
    vars::Vector{Variable},
    lb::Array{Float64},
    ub::Array{Float64},
)::Nothing
    t = threadid()
    MOI.delete(mip.optimizers[t], _bounds_constraint.(vars))
    funcs = MOI.SingleVariable[]
    sets = MOI.Interval[]
    for j = 1:length(vars)
        push!(funcs, MOI.SingleVariable(vars[j]))
        push!(sets, MOI.Interval(lb[j], ub[j]))
    end
    MOI.add_constraints(mip.optimizers[t], funcs, sets)
    return
end

"""
    name(mip::MIP, var::Variable)::String

Return the name of the decision variable `var`.
"""
function name(mip::MIP, var::Variable)::String
    t = threadid()
    return MOI.get(mip.optimizers[t], MOI.VariableName(), convert(MOI.VariableIndex, var))
end

convert(::Type{MOI.VariableIndex}, v::Variable) = MOI.VariableIndex(v.index)

"""
    probe(mip::MIP, var)::Tuple{Float64, Float64}

Suppose that the LP relaxation of `mip` has been solved and that `var` holds
a fractional value `f`. This function returns two numbers corresponding,
respectively, to the the optimal values of the LP relaxations having the
constraints var=1 and var=0 enforced. If any branch is infeasible, the optimal
value for that branch is Inf for minimization problems and -Inf for maximization
problems.
"""
function probe(mip::MIP, var)::Tuple{Float64,Float64}
    set_bounds!(mip, [var], [1.0], [1.0])
    status_up, obj_up = solve_relaxation!(mip)
    set_bounds!(mip, [var], [0.0], [0.0])
    status_down, obj_down = solve_relaxation!(mip)
    set_bounds!(mip, [var], [0.0], [1.0])
    return obj_up * mip.sense, obj_down * mip.sense
end


text:
using LinearAlgebra
using TensorOperations

####################################################################
#
# construct \Omega matrix
#

function Omega_uu(N::Int)
    id = Matrix{Float64}(I, N, N)
    zs = zeros(N,N)
    return [zs -id; id zs]
end

function Omega_dd(N::Int)
    id = Matrix{Float64}(I, N, N)
    zs = zeros(N,N)
    return [zs id; -id zs]
end

function change_coordinate(X::Array{<:Number},S::AbstractArray{<:Number,2})
    size_X =  size(X)
    N = div(size_X[1], 2)
    rank_X = ndims(X)
    X_prime = zeros(Complex{Float64},size_X...)
    indeces_set = Iterators.product(ntuple(i->1:2N,rank_X)...)

    for inds_after in indeces_set
        for inds_before in indeces_set
            Xprime_element = X[inds_before...];
            for i in 1:rank_X
                Xprime_element = Xprime_element * S[inds_before[i],inds_after[i]]
            end;
            X_prime[inds_after...] += Xprime_element
        end;
    end;
    return X_prime
end

############################################################
#symmetrize tensors
function symmetrize_tensor(X::Array{<:Number})
    rank_X = ndims(X)
    if rank_X ==1; return X; end;

    #construct permutation_group
    permutation_group = Set([])
    inds_for_perms = Iterators.product(ntuple(i->1:(rank_X-i+1),rank_X)...)
    size_perm_set = length(inds_for_perms)
    for perm in inds_for_perms
        perm = collect(perm)
        vec_perm  = zeros(Int64,rank_X) # permutated [1:N]
        #generate an element of the permutation group from perm
        for i in 1:rank_X
            j = 1
            while true
                if perm[i] == 1 && vec_perm[j] == 0; vec_perm[j]=i; break;
                elseif vec_perm[j] != 0; j+=1;
                else j+=1; perm[i]-=1
                end;
            end;
        end;
        push!(permutation_group,vec_perm)
    end;

    indeces_set = Iterators.product(ntuple(i->1:2N,rank_X)...)
    X_sym = zeros(eltype(X),size(X)...)
    for inds in indeces_set
        for perm in permutation_group
            inds_perm = Tuple(permute!([inds...],perm)) #permute inds by perm
            X_sym[inds...] += 1.0 / size_perm_set * X[inds_perm...]
        end;
    end;
    return X_sym
end;

############################################################
# get tensors of the effecitive Hamiltonian including Lindblad terms
#
function get_effective_Hamiltonian_tensors!(A′::Array{Complex{Float64},1},B′::Array{Complex{Float64},2},C′::Array{Complex{Float64},3},D′::Array{Complex{Float64},4},A::Array{Float64,1},B::Array{Float64,2},C::Array{Float64,3},D::Array{Float64,4},vec_EF::Array{Array{Complex{Float64}},2},μ₁::Array{Float64,1},μ₂::Array{Float64,2})

    N = Int(length(A)/2)
    Ω_dd = Omega_dd(N)
    A′ .= A; B′ .= B; C′ .= C; D′ .= D;

    for i in 1:size(vec_EF)[1]
        E = vec_EF[i,1];
        F = vec_EF[i,2];
        real_E = real.(E);
        real_F = real.(F);
        exp_Re_L = dot(real_E,μ₁) + tr(real_F*μ₂)
        @tensor begin
            A′[i]       += im * exp_Re_L * E[i]
            B′[i,j]     += im * exp_Re_L * F[i,j]
            B′[i,j]     -= im *  real_E[i] * E[j]
            C′[i,j,k]   -= im *  real_E[i] * F[j,k]
            C′[i,j,k]   -= im *  real_F[i,j] *E[k]
            D′[i,j,k,l] -= im *  real_F[i,j] * F[k,l]
        end
    end

    return (A′,symmetrize_tensor(B′),symmetrize_tensor(C′),symmetrize_tensor(D′))
end

############################################################
# get raw 2rd-order moments from central moments
#
function get_2nd_raw_moment(μ::Array{<:Number,1},μ_dd::Array{<:Number,2})

    @tensor begin
        μ′_dd[i,j]  := μ_dd[i,j] + μ[i] * μ[j]
    end

    return symmetrize_tensor(μ′_dd)
end
############################################################
# get raw 3rd-order moments from central moments
#
function get_3rd_raw_moment(μ::Array{<:Number,1},μ_dd::Array{<:Number,2},μ_ddd::Array{<:Number,3})

    @tensor begin
        μ′_ddd[i,j,k]  := μ_ddd[i,j,k] + 3.0 * μ_dd[i,j] * μ[k] + μ[i] * μ[j] * μ[k]
    end

    return symmetrize_tensor(μ′_ddd)
end

############################################################
# get raw 4th-order moments from central moments
#
function get_4th_raw_moment(μ::Array{<:Number,1},μ_dd::Array{<:Number,2},μ_ddd::Array{<:Number,3},μ_dddd::Array{<:Number,4})

    @tensor begin
        μ′_dddd[i,j,k,l]  := μ_ddd[i,j,k,l] + 4. * μ_ddd[i,j,k] * μ[l] + 6.0 * μ_dd[i,j] * μ[k] * μ[l] + μ[i] * μ[j] * μ[k]
    end

    return symmetrize_tensor(μ′_dddd)
end

#################################################################################################
#  get symmetrized z coordinate tensors from a coordinate tensors
#
#

function get_sym_tensors_with_z_coord(array_of_tensors::Array{Array{Complex{Float64}},1})

    ret = Array{Array{Float64}}(undef,length(array_of_tensors))
    for i in 1:length(ret)
        ret[i] = real.(symmetrize_tensor(converta2z(array_of_tensors[i])))
    end
    return Tuple(ret)
end

max_stars_repo_name: ash-y/QuantumCOMFrame

max_issues_repo_name: ash-y/QuantumCOMFrame

max_forks_repo_name: ash-y/QuantumCOMFrame


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

"""
    describe_group(group_id, identity_store_id)
    describe_group(group_id, identity_store_id, params::Dict{String,<:Any})

Retrieves the group metadata and attributes from GroupId in an identity store.

# Arguments
- `group_id`: The identifier for a group in the identity store.
- `identity_store_id`: The globally unique identifier for the identity store, such as
  d-1234567890. In this example, d- is a fixed prefix, and 1234567890 is a randomly generated
  string that contains number and lower case letters. This value is generated at the time
  that a new identity store is created.

"""
function describe_group(
    GroupId, IdentityStoreId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return identitystore(
        "DescribeGroup",
        Dict{String,Any}("GroupId" => GroupId, "IdentityStoreId" => IdentityStoreId);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_group(
    GroupId,
    IdentityStoreId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return identitystore(
        "DescribeGroup",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "GroupId" => GroupId, "IdentityStoreId" => IdentityStoreId
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_user(identity_store_id, user_id)
    describe_user(identity_store_id, user_id, params::Dict{String,<:Any})

Retrieves the user metadata and attributes from UserId in an identity store.

# Arguments
- `identity_store_id`: The globally unique identifier for the identity store, such as
  d-1234567890. In this example, d- is a fixed prefix, and 1234567890 is a randomly generated
  string that contains number and lower case letters. This value is generated at the time
  that a new identity store is created.
- `user_id`: The identifier for a user in the identity store.

"""
function describe_user(
    IdentityStoreId, UserId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return identitystore(
        "DescribeUser",
        Dict{String,Any}("IdentityStoreId" => IdentityStoreId, "UserId" => UserId);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_user(
    IdentityStoreId,
    UserId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return identitystore(
        "DescribeUser",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("IdentityStoreId" => IdentityStoreId, "UserId" => UserId),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_groups(identity_store_id)
    list_groups(identity_store_id, params::Dict{String,<:Any})

Lists the attribute name and value of the group that you specified in the search. We only
support DisplayName as a valid filter attribute path currently, and filter is required.
This API returns minimum attributes, including GroupId and group DisplayName in the
response.

# Arguments
- `identity_store_id`: The globally unique identifier for the identity store, such as
  d-1234567890. In this example, d- is a fixed prefix, and 1234567890 is a randomly generated
  string that contains number and lower case letters. This value is generated at the time
  that a new identity store is created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filters"`: A list of Filter objects, which is used in the ListUsers and ListGroups
  request.
- `"MaxResults"`: The maximum number of results to be returned per request. This parameter
  is used in the ListUsers and ListGroups request to specify how many results to return in
  one page. The length limit is 50 characters.
- `"NextToken"`: The pagination token used for the ListUsers and ListGroups API operations.
  This value is generated by the identity store service. It is returned in the API response
  if the total results are more than the size of one page. This token is also returned when
  it is used in the API request to search for the next page.
"""
function list_groups(IdentityStoreId; aws_config::AbstractAWSConfig=global_aws_config())
    return identitystore(
        "ListGroups",
        Dict{String,Any}("IdentityStoreId" => IdentityStoreId);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_groups(
    IdentityStoreId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return identitystore(
        "ListGroups",
        Dict{String,Any}(
            mergewith(
                _merge, Dict{String,Any}("IdentityStoreId" => IdentityStoreId), params
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_users(identity_store_id)
    list_users(identity_store_id, params::Dict{String,<:Any})

Lists the attribute name and value of the user that you specified in the search. We only
support UserName as a valid filter attribute path currently, and filter is required. This
API returns minimum attributes, including UserId and UserName in the response.

# Arguments
- `identity_store_id`: The globally unique identifier for the identity store, such as
  d-1234567890. In this example, d- is a fixed prefix, and 1234567890 is a randomly generated
  string that contains number and lower case letters. This value is generated at the time
  that a new identity store is created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filters"`: A list of Filter objects, which is used in the ListUsers and ListGroups
  request.
- `"MaxResults"`: The maximum number of results to be returned per request. This parameter
  is used in the ListUsers and ListGroups request to specify how many results to return in
  one page. The length limit is 50 characters.
- `"NextToken"`: The pagination token used for the ListUsers and ListGroups API operations.
  This value is generated by the identity store service. It is returned in the API response
  if the total results are more than the size of one page. This token is also returned when
  it is used in the API request to search for the next page.
"""
function list_users(IdentityStoreId; aws_config::AbstractAWSConfig=global_aws_config())
    return identitystore(
        "ListUsers",
        Dict{String,Any}("IdentityStoreId" => IdentityStoreId);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_users(
    IdentityStoreId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return identitystore(
        "ListUsers",
        Dict{String,Any}(
            mergewith(
                _merge, Dict{String,Any}("IdentityStoreId" => IdentityStoreId), params
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
abstract type AbstractPoint end

struct PointI <: AbstractPoint
    x::Int
    y::Int
end

mutable struct PointM <: AbstractPoint
    x::Int
    y::Int
end

PointM(p::PointM) = PointM(p.x, p.y)

d(p::AbstractPoint) = abs(p.x) + abs(p.y)
move(p::PointI, d::PointI) = PointI(p.x+d.x, p.y+d.y)
move(p::PointM, d::PointM) = (p.x += d.x; p.y += d.y; p)

function simI()
    maxd = 0
    x = PointI(0, 0)
    @inbounds for i in 1:10^6
        x = move(x, PointI(2rand(Bool)-1, 2rand(Bool)-1))
        curd = d(x)
        maxd = max(maxd, curd)
    end
    maxd
end

function simM()
    maxd = 0
    x = PointM(0, 0)
    m = PointM(0, 0)
    @inbounds for i in 1:10^6
        m.x, m.y = 2rand(Bool)-1, 2rand(Bool)-1
        move(x, m)
        curd = d(x)
        maxd = max(maxd, curd)
    end
    maxd
end

function simI2()
    path = PointI[]
    x = PointI(0, 0)
    @inbounds for i in 1:10^6
        push!(path, x)
        x = move(x, PointI(2rand(Bool)-1, 2rand(Bool)-1))
    end
    path
end

function simM2()
    path = PointM[]
    x = PointM(0, 0)
    m = PointM(0, 0)
    @inbounds for i in 1:10^6
        push!(path, PointM(x))
        m.x, m.y = 2rand(Bool)-1, 2rand(Bool)-1
        move(x, m)
    end
    path
end

max_stars_repo_path: Chapter06/07. The efficiency of mutable versus immutable types/walk.jl


text:
using BenchmarkTools
using PolygonOps
using StaticArrays
using Test

const suite = BenchmarkGroup()
suite["inpolygon"] = BenchmarkGroup()

circle(t) = SVector(cos(t), sin(t))

circle_poly = map(circle, 0:(2pi/1000000):2pi)

#@show length(circle_poly)

# using Plots
# display(plot(map(x -> x[1], circle_poly),map(x -> x[2], circle_poly)))

# end == begin
push!(circle_poly, circle_poly[1])
for algo in (HaoSun(),HormannAgathos())
    @test inpolygon(SVector(0.1,0),circle_poly,algo) == 1
    @test inpolygon(SVector(0.,2.),circle_poly,algo) == 0
    @test inpolygon(SVector(2.,0.),circle_poly,algo) == 0
    @test inpolygon(SVector(1.,0.),circle_poly,algo) == -1
    @test inpolygon(SVector(0.,-2.),circle_poly,algo) == 0
end

suite["inpolygon"]["circle, HaoSun"] = @benchmarkable inpolygon(SVector(0,0),circle_poly, HaoSun())
suite["inpolygon"]["circle, HormanAgathos"] = @benchmarkable inpolygon(SVector(0,0),circle_poly, HormannAgathos())

results = run(suite)

for trial in results
    ctx = IOContext(stdout, :verbose => true, :compact => false)
    println(ctx)
    println(ctx, trial.first)
    println(ctx, trial.second)
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# Benchmark implementing the board logic for the game of go and
# exercising it by playing random games. Derived from
# http://www.lysator.liu.se/~gunnar/gtp/brown-1.0.tar.gz
import Base.getindex

const EMPTY = 0
const WHITE = 1
const BLACK = 2

# Function to find the opposite color.
other_color(color::Int) = WHITE + BLACK - color

# Used in the final_status[] array.
const DEAD = 0
const ALIVE = 1
const SEKI = 2
const WHITE_TERRITORY = 3
const BLACK_TERRITORY = 4
const UNKNOWN = 5

type XorRand
  state::UInt32
end

function xor_srand(rand::XorRand, seed::UInt32)
  rand.state = seed
end

function xor_randn(rand::XorRand, n::UInt32)
  rand.state $= rand.state << 13
  rand.state $= rand.state >> 17
  rand.state $= rand.state << 5
  rand.state % n
end

xor_randn(rand::XorRand, n::Int) = convert(Int, xor_randn(rand, convert(UInt32, n)))

# Offsets for the four directly adjacent neighbors. Used for looping.
const deltai = (-1, 1, 0, 0)
const deltaj = (0, 0, -1, 1)
neighbor(i::Int, j::Int, k::Int) = (i + deltai[k], j + deltaj[k])

type Board
  size::Int
  komi::Float64

  # Board represented by a 1D array. The first board_size*board_size
  # elements are used. Vertices are indexed row by row, starting with 0
  # in the upper left corner.
  board::Matrix{Int}

  # Stones are linked together in a circular list for each string.
  next_stone::Matrix{Int}

  # Storage for final status computations.
  final_status::Matrix{Int}

  # Point which would be an illegal ko recapture.
  ko_i::Int
  ko_j::Int

  # xor-shift random number generator.
  rand::XorRand

  function Board(n::Int)
    init(new(), n, convert(UInt32, 2463534242))
  end
end

function init(board::Board, n::Int, seed::UInt32)
  board.size = n
  board.komi = 0.0
  board.board = zeros(Int, n, n)
  board.next_stone = zeros(Int, n, n)
  board.final_status = zeros(Int, n, n)
  board.ko_i = 0
  board.ko_j = 0
  board.rand = XorRand(seed)
  board
end

getindex(board::Board, pos::Int) = board.board[pos]
getindex(board::Board, i::Int, j::Int) = board.board[i, j]

# Functions to convert between 1D and 2D coordinates. The 2D coordinate
# (i, j) points to row i and column j, starting with (1,1) in the
# upper left corner.
POS(board::Board, i::Int, j::Int) = (j - 1) * board.size + i
IJ(board::Board, pos::Int) = (1 + mod((pos - 1), board.size), 1 + fld(pos - 1, board.size))

function set_komi(board::Board, komi::Float64)
  board.komi = komi
end

function set_random_seed(board::Board, seed::UInt32)
  xor_srand(board.rand, seed)
end

function clear(board::Board)
  board.board[:] = EMPTY
end

is_pass_move(i::Int, j::Int) = i == 0 && j == 0

function on_board(board::Board, i::Int, j::Int)
  i >= 1 && i <= board.size && j >= 1 && j <= board.size
end

function legal_move(board::Board, i::Int, j::Int, color::Int)
  other = other_color(color)

  # Pass is always legal.
  if is_pass_move(i, j)
    return true
  end

  # Already occupied.
  if board[i, j] != EMPTY
    return false
  end

  # Illegal ko recapture. It is not illegal to fill the ko so we must
  # check the color of at least one neighbor.
  if i == board.ko_i && j == board.ko_j && ((on_board(board, i - 1, j) && board[i - 1, j] == other) || (on_board(board, i + 1, j) && board[i + 1, j] == other))
    return false
  end

  true
end

# Does the string at (i, j) have any more liberty than the one at (libi, libj)?
function has_additional_liberty(board::Board, i::Int, j::Int, libi::Int, libj::Int)
  start = POS(board, i, j)
  pos = start
  while true
    (ai, aj) = IJ(board, pos)
    for k = 1:4
      (bi, bj) = neighbor(ai, aj, k)
      if on_board(board, bi, bj) && board[bi, bj] == EMPTY && (bi != libi || bj != libj)
        return true
      end
    end

    pos = board.next_stone[pos]
    if pos == start
      break
    end
  end

  false
end

# Does (ai, aj) provide a liberty for a stone at (i, j)?
function provides_liberty(board::Board, ai::Int, aj::Int, i::Int, j::Int, color::Int)
  # A vertex off the board does not provide a liberty.
  if !on_board(board, ai, aj)
    return false
  end

  # An empty vertex IS a liberty.
  if board[ai, aj] == EMPTY
    return true
  end

  # A friendly string provides a liberty to (i, j) if it currently
  # has more liberties than the one at (i, j).
  if board[ai, aj] == color
    return has_additional_liberty(board, ai, aj, i, j)
  end

  # An unfriendly string provides a liberty if and only if it is
  # captured, i.e. if it currently only has the liberty at (i, j).
  !has_additional_liberty(board, ai, aj, i, j)
end

# Is a move at ij suicide for color?
function suicide(board::Board, i::Int, j::Int, color::Int)
  for k = 1:4
    if provides_liberty(board, neighbor(i, j, k)..., i, j, color)
      return false
    end
  end
  true
end

# Remove a string from the board array. There is no need to modify
# the next_stone array since this only matters where there are
# stones present and the entire string is removed.
function remove_string(board::Board, i::Int, j::Int)
  start = POS(board, i, j)
  pos = start
  removed = 0
  while true
    board.board[pos] = EMPTY
    removed += 1
    pos = board.next_stone[pos]
    if pos == start
      break
    end
  end
  removed
end

# Do two vertices belong to the same string? It is required that both
# pos1 and pos2 point to vertices with stones.
function same_string(board::Board, pos1::Int, pos2::Int)
  pos = pos1
  while true
    if pos == pos2
      return true
    end
    pos = board.next_stone[pos]
    if pos == pos1
      break
    end
  end
  false
end

# Play at (i, j) for color. No legality check is done here. We need
# to properly update the board array, the next_stone array, and the
# ko point.
function play_move(board::Board, i::Int, j::Int, color::Int)
  pos = POS(board, i, j)
  captured_stones = 0

  # Reset the ko point.
  board.ko_i = 0
  board.ko_j = 0

  # Nothing more happens if the move was a pass.
  if is_pass_move(i, j)
    return
  end

  # If the move is a suicide we only need to remove the adjacent
  # friendly stones.
  if suicide(board, i, j, color)
    for k = 1:4
      (ai, aj) = neighbor(i, j, k)
      if on_board(board, ai, aj) && board[ai, aj] == color
        remove_string(board, ai, aj)
      end
    end
    return
  end

  # Not suicide. Remove captured opponent strings.
  for k = 1:4
    (ai, aj) = neighbor(i, j, k)
    if on_board(board, ai, aj) && board[ai, aj] == other_color(color) && !has_additional_liberty(board, ai, aj, i, j)
      captured_stones += remove_string(board, ai, aj)
    end
  end

  # Put down the new stone. Initially build a single stone string by
  # setting next_stone[pos] pointing to itself.
  board.board[pos] = color
  board.next_stone[pos] = pos

  # If we have friendly neighbor strings we need to link the strings
  # together.
  for k = 1:4
    (ai, aj) = neighbor(i, j, k)
    pos2 = POS(board, ai, aj)
    # Make sure that the stones are not already linked together. This
    # may happen if the same string neighbors the new stone in more
    # than one direction.
    if on_board(board, ai, aj) && board[pos2] == color && !same_string(board, pos, pos2)
      # The strings are linked together simply by swapping the the
      # next_stone pointers.
      (board.next_stone[pos], board.next_stone[pos2]) = (board.next_stone[pos2], board.next_stone[pos])
    end
  end

  # If we have captured exactly one stone and the new string is a
  # single stone it may have been a ko capture.
  if captured_stones == 1 && board.next_stone[pos] == pos
    # Check whether the new string has exactly one liberty. If so it
    # would be an illegal ko capture to play there immediately. We
    # know that there must be a liberty immediately adjacent to the
    # new stone since we captured one stone.
    for k = 1:4
      (ai, aj) = neighbor(i, j, k)
      if on_board(board, ai, aj) && board[ai, aj] == EMPTY
        if !has_additional_liberty(board, i, j, ai, aj)
          board.ko_i = ai
          board.ko_j = aj
        end
        break
      end
    end
  end
end

# Generate a move.
function generate_move(board::Board, color::Int)
  moves = zeros(Int, 2, board.size * board.size)
  num_moves = 0

  for ai = 1:board.size, aj = 1:board.size
    # Consider moving at (ai, aj) if it is legal and not suicide.
    if legal_move(board, ai, aj, color) && !suicide(board, ai, aj, color)
      # Further require the move not to be suicide for the opponent...
      if !suicide(board, ai, aj, other_color(color))
        num_moves += 1
        moves[:,num_moves] = [ai, aj]
      else
        # ...however, if the move captures at least one stone,
        # consider it anyway.
        for k = 1:4
          (bi, bj) = neighbor(ai, aj, k)
          if on_board(board, bi, bj) && board[bi, bj] == other_color(color)
            num_moves += 1
            moves[:,num_moves] = [ai, aj]
            break
          end
        end
      end
    end
  end

  # Choose one of the considered moves randomly with uniform
  # distribution. (Strictly speaking the moves with smaller 1D
  # coordinates tend to have a very slightly higher probability to be
  # chosen, but for all practical purposes we get a uniform
  # distribution.)
  if num_moves > 0
    move = moves[:,1 + xor_randn(board.rand, num_moves)]
    return (move[1], move[2])
  else
    # But pass if no move was considered.
    return (0, 0)
  end
end

# Set a final status value for an entire string.
function set_final_status_string(board::Board, pos::Int, status::Int)
  start = pos
  while true
    board.final_status[pos] = status
    pos = board.next_stone[pos]
    if pos == start
      break
    end
  end
end

# Compute final status. This function is only valid to call in a
# position where generate_move() would return pass for at least one
# color.
#
# Due to the nature of the move generation algorithm, the final
# status of stones can be determined by a very simple algorithm:
#
# 1. Stones with two or more liberties are alive with territory.
# 2. Stones in atari are dead.
#
# Moreover alive stones are unconditionally alive even if the
# opponent is allowed an arbitrary number of consecutive moves.
# Similarly dead stones cannot be brought alive even by an arbitrary
# number of consecutive moves.
#
# Seki is not an option. The move generation algorithm would never
# leave a seki on the board.
#
# Comment: This algorithm doesn't work properly if the game ends with
#          an unfilled ko. If three passes are required for game end,
#          that will not happen.
#
function compute_final_status(board::Board)
  board.final_status[:] = UNKNOWN

  for i = 1:board.size, j = 1:board.size
    if board[i, j] == EMPTY
      for k = 1:4
        (ai, aj) = neighbor(i, j, k)
        if !on_board(board, ai, aj)
          continue
        end
        # When the game is finished, we know for sure that (ai, aj)
        # contains a stone. The move generation algorithm would
        # never leave two adjacent empty vertices. Check the number
        # of liberties to decide its status, unless it's known
        # already.
        #
        # If we should be called in a non-final position, just make
        # sure we don't call set_final_status_string() on an empty
        # vertex.
        pos = POS(board, ai, aj)
        if board.final_status[ai, aj] == UNKNOWN
          if board[ai, aj] != EMPTY
            if has_additional_liberty(board, ai, aj, i, j)
              set_final_status_string(board, pos, ALIVE)
            else
              set_final_status_string(board, pos, DEAD)
            end
          end
        end
        # Set the final status of the pos vertex to either black
        # or white territory.
        if board.final_status[i, j] == UNKNOWN
          if (board.final_status[ai, aj] == ALIVE) $ (board[ai, aj] == WHITE)
            board.final_status[i, j] = BLACK_TERRITORY
          else
            board.final_status[i, j] = WHITE_TERRITORY
          end
        end
      end
    end
  end
end

get_final_status(board::Board, i::Int, j::Int) = board.final_status[i, j]

function compute_score(board::Board)
  score = board.komi
  compute_final_status(board)
  for i = 1:board.size, j = 1:board.size
    status = get_final_status(board, i, j)
    if status == BLACK_TERRITORY
      score -= 1.0
    elseif status == WHITE_TERRITORY
      score += 1.0
    elseif (status == ALIVE) $ (board[i, j] == WHITE)
      score -= 1.0
    else
      score += 1.0
    end
  end
  score
end

function benchmark(num_games_per_point::Int)
  random_seed = convert(UInt32, 1)
  board_size = 9
  komi = 0.5

  board = Board(board_size)
  set_komi(board, komi)
  set_random_seed(board, random_seed)

  for i = 1:board.size, j = 1:board.size
    white_wins = 0
    black_wins = 0
    for k = 1:num_games_per_point
      passes = 0
      num_moves = 1
      color = WHITE
      clear(board)
      play_move(board, i, j, BLACK)
      while passes < 3 && num_moves < 600
        (movei, movej) = generate_move(board, color)
        play_move(board, movei, movej, color)
        if is_pass_move(movei, movej)
          passes += 1
        else
          passes = 0
        end
        num_moves += 1
        color = other_color(color)
      end
      if passes == 3
        if compute_score(board) > 0
          white_wins += 1
        else
          black_wins += 1
        end
      end
    end
#    @printf("%d %d %f\n", i - 1, j - 1, black_wins / (black_wins + white_wins))
  end
end

function main(args)
  n = 10
  if length(args) > 0
    n = parseint(args[1])
  end
  @time benchmark(n)
end


text:
module TestTools

using Test
using Try
using TryExperimental
using TryExperimental: @and_then, @or_else

if !@isdefined Returns
    Returns(y) = function constant(_args...; _kwrags...)
        return y
    end
end

function trygetnitems(xs)
    Try.@and_return trygetlength(xs)
    Ok(count(Returns(true), xs))
end

function nitems(xs)
    Try.or_else(trygetlength(xs)) do _
        Ok(count(Returns(true), xs))
    end |> Try.unwrap
end

function nitems2(xs)
    Try.@return trygetlength(xs)
    count(Returns(true), xs)
end

function test_and_return()
    @test Try.unwrap(trygetnitems(1:3)) == 3

    ch = foldl(push!, 1:3; init = Channel{Int}(3))
    close(ch)
    @test Try.unwrap(trygetnitems(ch)) == 3
end

function check_nitems(nitems)
    @test nitems(1:3) == 3

    ch = foldl(push!, 1:3; init = Channel{Int}(3))
    close(ch)
    @test nitems(ch) == 3
end

test_or_else() = check_nitems(nitems)
test_return() = check_nitems(nitems2)

try_map_prealloc(f, xs) =
    Try.and_then(trygetlength(xs)) do n
        Try.and_then(trygeteltype(xs)) do T
            ys = Vector{T}(undef, n)
            for (i, x) in zip(eachindex(ys), xs)
                ys[i] = f(x)
            end
            return Ok(ys)
        end
    end

function try_map_prealloc2(f, xs)
    T = @? trygeteltype(xs)  # macro-based short-circuiting
    n = @? trygetlength(xs)
    ys = Vector{T}(undef, n)
    for (i, x) in zip(eachindex(ys), xs)
        ys[i] = f(x)
    end
    return Ok(ys)
end

function test_and_then()
    @test Try.unwrap(try_map_prealloc(x -> x + 1, 1:3)) == 2:4
end

function test_or_return()
    @test Try.unwrap(try_map_prealloc2(x -> x + 1, 1:3)) == 2:4
end

function test_and()
    @test Try.and(Ok(1), Ok(2), Ok(3)) == Try.Ok(3)
    @test Try.and(Ok(1), Err(2), Ok(3)) == Try.Err(2)
    @test Try.and(Some(1), Some(2), Some(3)) == Some(3)
    @test Try.and(Some(1), nothing, Some(3)) === nothing
    @test Try.@and(Ok(1), Ok(2), Ok(3)) == Try.Ok(3)
    @test Try.@and(Ok(1), Err(2), Ok(3)) == Try.Err(2)
    @test Try.@and(Some(1), Some(2), Some(3)) == Some(3)
    @test Try.@and(Some(1), nothing, Some(3)) === nothing
end

function test_or()
    @test Try.or(Err(1), Ok(2), Err(3)) == Try.Ok(2)
    @test Try.or(Err(1), Err(2), Err(3)) == Try.Err(3)
    @test Try.or(nothing, Some(2), Some(3)) == Some(2)
    @test Try.or(nothing, nothing, nothing) === nothing
    @test Try.@or(Err(1), Ok(2), Err(3)) == Try.Ok(2)
    @test Try.@or(Err(1), Err(2), Err(3)) == Try.Err(3)
    @test Try.@or(nothing, Some(2), Some(3)) == Some(2)
    @test Try.@or(nothing, nothing, nothing) === nothing
end

function test_unwrap_or_else()
    @test Try.unwrap_or_else(length, Try.Ok(1)) == 1
    @test Try.unwrap_or_else(length, Try.Err("four")) == 4
end

function test_curry()
    value =
        tryconvert(String, 1) |>
        Try.or_else() do _
            Ok("123")
        end |>
        Try.and_then() do x
            y = tryparse(Int, x)
            if y === nothing
                Err(ErrorException(""))
            else
                Ok(y)
            end
        end |>
        Try.unwrap

    @test value == 123
end

function demo_macro(xs)
    i = firstindex(xs)
    y = nothing
    while true
        x = @or_else(trygetindex(xs, i)) do _
            return :oob
        end
        @and_then(trygetindex(xs, Try.unwrap(x))) do z
            if z > 1
                y = z
                break
            end
        end
        i += 1
    end
    return y
end

function test_macro()
    @test demo_macro(1:2:10) == 5
    @test demo_macro(1:0) === :oob
end

function test_map()
    @test Try.map(x -> x + 1, Ok(1)) == Try.Ok(2)
    @test Try.map(x -> x + 1, Err(KeyError(:a))) == Err(KeyError(:a))
    @test Try.map(x -> x + 1, Some(1)) == Some(2)
    @test Try.map(x -> x + 1, nothing) === nothing
end

end  # module


text:
using Ratios
using Compat.Test

r = SimpleRatio(1,2)
@test convert(Float64, r) == 0.5
@test convert(Float32, r) == 0.5f0
@test convert(BigFloat, r) == BigFloat(1)/2

r2 = SimpleRatio(2,3)
@test r*r2 == SimpleRatio(2,6) == SimpleRatio(1,3)
@test r2*3 == 3*r2 == 2
@test r*false == false*r == 0
@test r/r2 == SimpleRatio(3,4)
@test r/2 == SimpleRatio(1,4)
@test 2/r == 4
@test 4 == 2/r
@test r+1 == 1+r == SimpleRatio(3,2)
@test r-1 == SimpleRatio(-1,2)
@test 1-r == r
@test r+r2 == SimpleRatio(7,6)
@test r-r2 == SimpleRatio(-1,6)
@test r^2 == SimpleRatio(1,4)
@test -r == SimpleRatio(-1,2)
@test 0.2*r ≈ 0.1
@test r == 0.5
@test 0.5 == r

@test_throws OverflowError -SimpleRatio(0x02,0x03)

@test r + SimpleRatio(0x02,0x03) == SimpleRatio(7,6)

@test SimpleRatio(11, 10) == 11//10
@test 1//3 + SimpleRatio(1, 5) == 8//15

@test isfinite(SimpleRatio(0,0)) == false
@test isfinite(SimpleRatio(1,0)) == false
@test isfinite(SimpleRatio(2,1)) == true


text:
using DoubleExponentialFormulas
using LinearAlgebra: norm
using Test

quadts32 = QuadTS(Float32)
quadts64 = QuadTS(Float64)
quadtsBF = QuadTS(BigFloat)

# Test Base.show works without error
let
    io = IOBuffer()
    @test (Base.show(io, MIME("text/plain"), quadts32); true)
    @test (Base.show(io, MIME("text/plain"), quadts64); true)
    @test (Base.show(io, MIME("text/plain"), quadtsBF); true)
end

let
    f(x) = 1
    expect = 2

    I, E = quadts32(f)
    @test eltype(I) == Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadts64(f)
    @test eltype(I) == Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadtsBF(f)
    @test eltype(I) == BigFloat
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)
end

let
    f(x) = x
    expect = 0

    I, E = quadts32(f)
    @test eltype(I) == Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadts64(f)
    @test eltype(I) == Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadtsBF(f)
    @test eltype(I) == BigFloat
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)
end

let
    f(x) = 3*x^2
    expect = 2

    I, E = quadts32(f)
    @test eltype(I) == Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadts64(f)
    @test eltype(I) == Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadtsBF(f)
    @test eltype(I) == BigFloat
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)
end

let
    f(x) = 2/(1 + x^2)
    expect = π

    I, E = quadts32(f)
    @test eltype(I) == Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadts64(f)
    @test eltype(I) == Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadtsBF(f)
    @test eltype(I) == BigFloat
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)
end

# Singular points at x = ±1
let
    f(x) = (2 - x)/sqrt(1 - x^2)
    expect = 2π

    I, E = quadts32(f)
    @test eltype(I) == Float32
    @test isapprox(I, expect, rtol=1e-3)
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadts64(f)
    @test eltype(I) == Float64
    @test isapprox(I, expect, rtol=1e-7)
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadtsBF(f)
    @test eltype(I) == BigFloat
    @test isapprox(I, expect, rtol=1e-16)
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)
end

# Test non-scalar output with QuadTS
let
    f(x::AbstractFloat) = [1/(1 + x^2), 2/(1 + x^2)]
    expect = [BigFloat(π)/2, BigFloat(π)]

    I, E = quadts32(f)
    @test eltype(I) == Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadts64(f)
    @test eltype(I) == Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadtsBF(f)
    @test eltype(I) == BigFloat
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)
end


text:
using Distributions, SpecialFunctions, Random

# Functions for Categorical task simulations (Section 2.3.2)

# Run one learner on one Categorical environment
function run!(envtype::Multi, learnertype, learnerparam; ns = 5,
                    stochasticity = 0.1, changeprobability = 0.001,
                    nsteps = 10^3, iseedenvpolicy = 3, iseedlearner = 3)

    env, learner, callback = initialize(envtype, learnertype, learnerparam,
                                ns = ns, stochasticity = stochasticity,
                                changeprobability = changeprobability,
                                nsteps = nsteps,
                                iseedenvpolicy = iseedenvpolicy,
                                iseedlearner = iseedlearner)
    learn!(learner, env, callback)
    callback
 end

# Initialize one Categorical learner and one Categorical environment
function initialize(envtype::Multi, learnertype, learnerparameterset; ns = 5,
                      stochasticity = 0.1, changeprobability = 0.001,
                      nsteps = 10^3, iseedenvpolicy = 3,
                      iseedlearner = 3)

    env = setupenv(envtype, ns = ns, stochasticity = stochasticity,
                      changeprobability = changeprobability,
                      nsteps = nsteps,
                      iseedenvpolicy = iseedenvpolicy)

    learner = setuplearner(envtype, learnertype, learnerparameterset, ns = ns,
                      iseedlearner = iseedlearner)
    callback = CallbackMultinomial()
    env, learner, callback
end

# Set up a Categorical environment
function setupenv(envtype::Multi; ns = 5, stochasticity = 0.1, changeprobability = 0.001,
                    nsteps = 10^3, iseedenvpolicy = 3)

    Random.seed!(ENV_RNG, iseedenvpolicy)
    env = EnvMultinomial(ns = ns, stochasticity = stochasticity,
                    changeprobability = changeprobability,
                    nsteps = nsteps,
                    seed = iseedenvpolicy)
end

# Set up a Categorical learner
function setuplearner(envtype::Multi, learnertype, learnerparameterset;
                    ns = 5, iseedlearner = 3)
    if learnertype == TParticleFilter
        learner = TParticleFilter(ns = ns,
                                nparticles = learnerparameterset[1],
                                stochasticity = learnerparameterset[2],
                                changeprobability = learnerparameterset[3],
                                seed = iseedlearner)
    elseif learnertype == TLeakyIntegrator
        learner = TLeakyIntegrator(ns = ns, etaleak = learnerparameterset[1])
    elseif learnertype == TSmileOriginal
        learner = TSmileOriginal(ns = ns, m = learnerparameterset[1],
                        stochasticity = learnerparameterset[2])
    elseif learnertype == TVarSmile
        learner = TVarSmile(ns = ns, m = learnerparameterset[1],
                                stochasticity = learnerparameterset[2])
    elseif learnertype == TBayesFilter
        learner = TBayesFilter(ns = ns,
                            stochasticity = learnerparameterset[1],
                            changeprobability = learnerparameterset[2],
                            cutoff = 1e-16)
                            #cutoff = 0.)
    elseif learnertype == TMPN
        learner = TMPN(ns = ns,
                        N = learnerparameterset[1],
                        stochasticity = learnerparameterset[2],
                        changeprobability = learnerparameterset[3])
    elseif learnertype == TSOR
        learner = TSOR(ns = ns,
                        N = learnerparameterset[1],
                        stochasticity = learnerparameterset[2],
                        changeprobability = learnerparameterset[3],
                        seed = iseedlearner)
    end
    learner
end


text:
###################################################
## network/solvers/lp.jl
## Simple LP method for travel time inference
###################################################

"""
    `lpTimes`, optimize travel times to minimize MAPE error

    `continuityConstraint` can be
    - "none" : no constraint
    - "simple": simple ϵ-based continuity
    - "neighborhoods": continuity constraint for each neighborhood
"""
function lpTimes(s::IterativeState;
                 continuityConstraint::AbstractString="none",
                 velocityBound::Float64 = 0.1, args...) #args is solver args
    g = s.data.network.graph
    paths = s.paths
    tripData = s.trips
    roads = s.data.network.roads

    #Create the model (will be changed to avoid hard-coded parameters)
    # !BarConvTol needs to be changed
    m = Model(solver = GurobiSolver(TimeLimit=10000,
                                    Crossover = 0,
                                    Method=2,
                                    BarConvTol=1e-8; args...))

    # DECISION VARIABLES
    # Road times
    @variable(m, t[i=vertices(g), j=outneighbors(g,i)] >= s.data.minTimes[i,j])
    # Absolute difference between tripData times and computed times
    @variable(m, epsilon[d=eachindex(tripData)] >= 0)

    if continuityConstraint == "neighborhoods"
        edgesAndNeighbors = flatten(
        [
            [
                (src(edge), dst(edge), src(nearEdge), dst(nearEdge))
                 for nearEdge in findNearEdgesSameType(s.data.network, edge)
            ]
         for edge in edges(g)
        ])
        @variable(m, velocity[(i,j,p,q) = edgesAndNeighbors] >= 0)
    end

    # OBJECTIVE
    @objective(m, Min, sum(tripData[d].weight/tripData[d].time * epsilon[d] for d in eachindex(tripData)))

    # CONSTRAINTS
    # absolute values contraints (define epsilon), equal to time of first path
    @constraint(m, epsLower[d=eachindex(tripData)],
        sum(paths[d][1][edge] * t[src(edge), dst(edge)] for edge in keys(paths[d][1])) - tripData[d].time >=
        - epsilon[d])
    @constraint(m, epsUpper[d=eachindex(tripData)],
        sum(paths[d][1][edge] * t[src(edge), dst(edge)] for edge in keys(paths[d][1])) - tripData[d].time <=
        epsilon[d])

    # inequality constraints
    @constraint(m, inequalityPath[d=eachindex(tripData), p=1:(length(paths[d])-1)],
        sum(paths[d][p+1][edge] * t[src(edge), dst(edge)] for edge in keys(paths[d][p+1])) >=
        sum(paths[d][1][edge] * t[src(edge), dst(edge)] for edge in keys(paths[d][1]))
        )

    if continuityConstraint == "simple"
        for i in vertices(g), j in outneighbors(g,i)
            for edge in findNearEdgesSameType(s.data.network, Edge(i,j))
                p = src(edge)
                q = dst(edge)
                @constraint(m, t[i,j]/roads[i,j].distance - t[p,q]/roads[p,q].distance
                    <= velocityBound * (roads[i,j].distance + roads[p,q].distance) / 2)
                @constraint(m, t[i,j]/roads[i,j].distance - t[p,q]/roads[p,q].distance
                    >= - velocityBound * (roads[i,j].distance + roads[p,q].distance) / 2)
            end
        end
    elseif continuityConstraint == "neighborhoods"
        clusters = clusterEdges(s.data.network)
        for (k, cluster) in enumerate(clusters)
            println(k)
            for edge in cluster
                i=src(edge); j=dst(edge);
                for nearEdge in findNearEdgesSameType(s.data.network, edge)
                    p=src(nearEdge); q=dst(nearEdge);
                    @constraint(m, t[i,j]/roads[i,j].distance - t[p,q]/roads[p,q].distance
                        <= velocity[(i,j,p,q)] * (roads[i,j].distance + roads[p,q].distance) / 2)
                    @constraint(m, t[i,j]/roads[i,j].distance - t[p,q]/roads[p,q].distance
                        >= - velocity[(i,j,p,q)] * (roads[i,j].distance + roads[p,q].distance) / 2)
                end
            end
            vList = flatten(
            [
                [
                    (src(edge), dst(edge), src(nearEdge), dst(nearEdge))
                     for nearEdge in findNearEdgesSameType(s.data.network, edge)
                ]
                 for edge in cluster
            ])
            @constraint(m, sum(velocity[(i,j,p,q)] for (i,j,p,q) in vList) <= length(vList) * velocityBound)
        end
    end

    # SOLVE LP
    status = solve(m)
    times = getvalue(t)

    # Export result as sparse matrix
    result = spzeros(Float64, nv(g), nv(g))
    for i in vertices(g), j in outneighbors(g,i)
        result[i,j] = times[i,j]
    end

    return result
end


text:
module GPMethylation

using Statistics
using DataFrames, CSV, RData
using GaussianProcesses
using Distributions
using LinearAlgebra
using ThreadTools
using RCall
using ProgressMeter

export  loaddata, loadsamplemeta, loadprobemeta, loadall, loadmeta,
        gpreg, gpreg_const, gpreg_linear, gpreg_matern52, gpreg_all_threads, gpmodels,  gpmodels_prior, gpreg_matern52_prior, gpreg_linear_prior,
        run_gpregression, saverdata, gpstats, loadresults

include("loaddata.jl")
include("gp.jl")
include("analysis.jl")

end


text:
push!(LOAD_PATH,"../src/") # load module from local directory

using PolaronMobility
using Printf # used for some stdout

using Test

@testset "PolaronMobility" begin

include("FrohlichAlpha.jl") # Simple test of Frohlich alpha vs. literature values
include("FeynmanAthermal.jl") # Athermal Feynman tests
include("HellwarthEffectiveFrequency.jl") # Test Hellwarth et al. 1999 PRB 'B' multiple branch reduction scheme
include("FrostPolaronMobility2017.jl") # Reproduce values published in Frost 2017 PRB
include("MultipleBranches.jl") # Test explicit Oct 2019:-> explicit phonon branches
include("MultiplePhonons.jl") # Test 2021 work on multiple phonon modes
end

println("\nThat's me! If I finished without interupting, all tests have passed.")


text:
"""
Convert a vector `w` to a symmetric matrix `P` as
P = [
    w[1] ⋅ ⋅ ⋅ ⋅ ...;
    w[2] w[n+1]  ...;
    w[3] w[n+2]  ...;
    ...  ...        ;
    w[n] ...        ;
]
# Notes
It should be compatible with `convert_quadratic_to_linear_basis`.
"""
function convert_to_matrix(w::AbstractVector)::AbstractMatrix
    # n(n+1)/2 x 1 -> n x n
    _l = length(w)
    l = Int((-1 + sqrt(1+4*2*_l)) / 2)  # P is l x l lower triangular matrix with elements of w
    P = zeros(l, l)
    idx = 1
    for i in 1:l
        for j in i:l
            if i == j
                P[j, i] = w[idx]
            else
                P[j, i] = 0.5 * w[idx]
                P[i, j] = 0.5 * w[idx]
            end
            idx += 1
        end
    end
    P
    # # n^2 x 1 -> n x n for convenience...
    # l = Int(sqrt(length(w)))
    # P = zeros(l, l)
    # idx = 1
    # for i in 1:l
    #     for j in 1:l
    #         P[j, i] = w[idx]
    #         idx += 1
    #     end
    # end
    # P
end

"""
Basis transformation for linear ADP and IRL;
xᵀ * P * x = wᵀ * ϕ(x)
where ϕ(x) = [x1^2, x2*x1, x3*x1, ..., xn*x1, x2^2, x3*x2, x4*x2, ..., xn*x2, x3^2, x4*x3, ...]
# Notes
It should be compatible with `convert_to_matrix`.
"""
function convert_quadratic_to_linear_basis(x::AbstractVector)
    ## Notes
    ## ϕ(x) = (I_m ⊗ x) x;
    ## Px = (I_m ⊗ xᵀ) w
    ## where P ∈ ℝ^{m×n}, w ∈ ℝ^{nm×1}, I_m: m×m identity matrix, ⊗: Kronecker product
    # n = length(x)
    # kron(Matrix(I, n, n), x) * x
    ## for dimension reduction
    n = length(x)
    ϕ = zeros(Int(n*(n+1)/2))
    idx = 1
    for i in 1:n
        for j in i:n
            ϕ[idx] = x[j] * x[i]
            idx += 1
        end
    end
    ϕ
end


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

##==============================================================================
## GenericFunctionNodeData
##==============================================================================

## COMMON
# getSolveInProgress
# isSolveInProgress

#TODO  getFactorFunction = getFactorType
"""
    $SIGNATURES

Return reference to the user factor in `<:AbstractDFG` identified by `::Symbol`.
"""
getFactorFunction(fcd::GenericFunctionNodeData) = fcd.fnc.usrfnc!
getFactorFunction(fc::DFGFactor) = getFactorFunction(getSolverData(fc))
getFactorFunction(dfg::AbstractDFG, fsym::Symbol) = getFactorFunction(getFactor(dfg, fsym))

"""
    $SIGNATURES

Return user factor type from factor graph identified by label `::Symbol`.

Notes
- Replaces older `getfnctype`.
"""
getFactorType(data::GenericFunctionNodeData) = data.fnc.usrfnc!
getFactorType(fct::DFGFactor) = getFactorType(getSolverData(fct))
getFactorType(dfg::AbstractDFG, lbl::Symbol) = getFactorType(getFactor(dfg, lbl))

"""
    $SIGNATURES

If you know a variable is `::Type{<:Pose2}` but want to find its default prior `::Type{<:PriorPose2}`.

Assumptions
- The prior type will be defined in the same module as the variable type.
- Not exported per default, but can be used with knowledge of the caveats.

Example
```julia
using RoME
@assert RoME.PriorPose2 == DFG._getPriorType(Pose2)
```
"""
_getPriorType(_type::Type{<:InferenceVariable}) = getfield(_type.name.module, Symbol(:Prior, _type.name.name))

##==============================================================================
## Factors
##==============================================================================
# |                   | label | tags | timestamp | solvable | solverData |
# |-------------------|:-----:|:----:|:---------:|:--------:|:----------:|
# | SkeletonDFGFactor |   X   |   x  |           |          |            |
# | DFGFactorSummary  |   X   |   X  |     X     |          |            |
# | DFGFactor         |   X   |   X  |     X     |     X    |      X     |

##------------------------------------------------------------------------------
## label
##------------------------------------------------------------------------------

## COMMON
# getLabel

##------------------------------------------------------------------------------
## tags
##------------------------------------------------------------------------------

## COMMON
# getTags
# setTags!

##------------------------------------------------------------------------------
## timestamp
##------------------------------------------------------------------------------

## COMMON
# getTimestamp

setTimestamp(f::AbstractDFGFactor, ts::DateTime, timezone=localzone()) = setTimestamp(f, ZonedDateTime(ts,  timezone))
setTimestamp(f::DFGFactor, ts::ZonedDateTime) = DFGFactor(f.label, ts, f.nstime, f.tags, f.solverData, f.solvable, getfield(f,:_variableOrderSymbols))
setTimestamp(f::DFGFactorSummary, ts::ZonedDateTime) = DFGFactorSummary(f.label, ts, f.tags, f._variableOrderSymbols)

##------------------------------------------------------------------------------
## solvable
##------------------------------------------------------------------------------

## COMMON
# getSolvable
# setSolvable!
# isSolvable

##------------------------------------------------------------------------------
## solvable
##------------------------------------------------------------------------------

## COMMON

##------------------------------------------------------------------------------
## _variableOrderSymbols
##------------------------------------------------------------------------------

#TODO perhaps making _variableOrderSymbols imutable (NTuple) will be a save option
"""
$SIGNATURES

Get the variable ordering for this factor.
Should be equivalent to getNeighbors unless something was deleted in the graph.
"""
getVariableOrder(fct::DFGFactor) = fct._variableOrderSymbols::Vector{Symbol}
getVariableOrder(dfg::AbstractDFG, fct::Symbol) = getVariableOrder(getFactor(dfg, fct))

##------------------------------------------------------------------------------
## solverData
##------------------------------------------------------------------------------

"""
    $SIGNATURES

Retrieve solver data structure stored in a factor.
"""
function getSolverData(f::F) where F <: DFGFactor
  return f.solverData
end

setSolverData!(f::DFGFactor, data::GenericFunctionNodeData) = f.solverData = data

##------------------------------------------------------------------------------
## utility
##------------------------------------------------------------------------------

"""
    $SIGNATURES

Return `::Bool` on whether given factor `fc::Symbol` is a prior in factor graph `dfg`.
"""
function isPrior(dfg::AbstractDFG, fc::Symbol)
  fco = getFactor(dfg, fc)
  isPrior(getFactorType(fco))
end

function isPrior(::AbstractPrior)
  return true
end

function isPrior(::AbstractRelative)
  return false
end

##==============================================================================
## Layer 2 CRUD (none) and Sets
##==============================================================================

##==============================================================================
## TAGS - See CommonAccessors
##==============================================================================


text:
using NonLinearReactionAdvectionDiffusionWithFrontData;
using LinearAlgebra;

u_l(t) = -8 + 4*sin(2*π / T * t);# Прямая задача может быть с неоднородными ГУ
u_r(t) =  4 + sin(-2*π / T * t);# Но в дальнейшем, будем использовать только однородные.
qf(x) = 4*sin(3 * π * x);        # Коэффициент линейного усиления, который в обратной
# задаче необходимо определить, но при генерации априорной
# информации мы задаем некоторый коэффициент, который,
# собственно, после имея априорную информацию и будем определять.
ε = 0.2;                        # Малый параметр при старшей производной
a, b = 0, 1;                    # Область по X
t₀, T = 0, 0.28;                # Область по T
N, M = 50, 80;                  # Кол-во разбиений по X, T
h = (b-a)/N;                    # шаг по X
τ = (T-t₀)/M;                   # шаг по T
Xₙ = [a  + n*h for n in 0:N];   # Сетка по Х
Tₘ = [t₀ + m*τ for m in 0:M];   # Сетка по Т
qₙ =      qf.(Xₙ);               # Сеточные значения коэффициента линейного усиления
ulₘ=    u_l.(Tₘ);               # Сеточные значения левого  ГУ
urₘ=    u_r.(Tₘ);               # Сеточные значения правого ГУ
y₀ = u_init.(Xₙ);               # Начальные условия
nothing #hide

# Подготовим массивы, выбросив граничные точки
# ведь тестируемая функция — для внутреннего использования
qq = qₙ[2:N];
y = y₀[2:N];

dl, d, du = NonLinearReactionAdvectionDiffusionWithFrontData.DRP_y(y, 1, Xₙ, N, ε, ulₘ, urₘ, qq)

@test length(dl) == N - 2
@test length(d)  == N - 1
@test length(du) == N - 2

@test count(x -> x == 0, dl) == 0
@test count(x -> x == 0, d)  == 0
@test count(x -> x == 0, du) == 0

tdjac = NonLinearReactionAdvectionDiffusionWithFrontData.∂DRP_∂y(y, 1, Xₙ, N, ε, ulₘ, urₘ, qq)
@test typeof(tdjac) <: Tridiagonal
@test tdjac == Tridiagonal( dl, d, du )

# Сравним с якобианом автодифференцирования
jac = NonLinearReactionAdvectionDiffusionWithFrontData.∂directRP_∂y(y, 1, Xₙ, N, ε, ulₘ, urₘ, qq)
@test isapprox(tdjac, jac)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule AMReX_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("AMReX")
JLLWrappers.@generate_main_file("AMReX", UUID("9b6a1cb2-ee7e-5f22-a98b-ab382f85f133"))
end  # module AMReX_jll