text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
module Q

using InteractiveUtils: subtypes

abstract type Creature end

Base.@kwdef mutable struct _Common
    hp::Int = 200
    atk::Int = 100
end

Base.@kwdef mutable struct Goblin <: Creature
    _common::_Common = _Common(88, 77)
    weapon::String = "dagger"
end

Base.@kwdef mutable struct Slime <: Creature
    _common::_Common = _Common(55, 44)
    sp::String = "acid"
end

for T in subtypes(Creature)
    n = length(fieldnames(_Common))
    @eval $(nameof(T))(a...) = $T(_Common(a[1:$n]...), a[$n+1:end]...)
    @eval Base.getproperty(x::$T, p::Symbol) =
        p ∈ fieldnames(_Common) ? getfield(x._common, p) : getfield(x, p)
    @eval Base.setproperty!(x::$T, p::Symbol, v) =
        p ∈ fieldnames(_Common) ? setfield!(x._common, p, v) : setfield!(x, p, v)
    @eval Base.propertynames(x::$T) = (fieldnames(_Common)...,
        (p for p in fieldnames($T) if first(string(p)) != '_')...)
    @eval function Base.show(io::IO, x::$T)
        props = getproperty.(Ref(x), propertynames(x))
        print(io, nameof($T), '(', repr(first(props)))
        for p in props[2:end] print(io, ", ", repr(p)) end
        print(io, ')')
    end
end

end

@show Q.Goblin()
@show propertynames(Q.Goblin())
@show Q.Slime()
@show propertynames(Q.Slime())
println()
@show g = Q.Goblin(123, 45, "dirty dagger")
@show g.hp
@show g.atk
@show g.weapon
@show g.hp = 99
@show g
@show g.atk = 11
@show g
@show g.weapon = "nothing"
@show g
println()
@show s = Q.Slime(123, 45, "weak acid")
@show s.hp
@show s.atk
@show s.sp
@show s.hp = 66
@show s
@show s.atk = 99
@show s
@show s.sp = "strong acid"
@show s;

# %%
T = Q.Goblin
n = 2
:($(nameof(T))(a...) = $T(_Common(a[1:$n]...), a[$n+1:end]...)) |> Base.remove_linenums! |> print

# %%
:(Base.getproperty(x::$T, p::Symbol) =
    p ∈ fieldnames(_Common) ? getfield(x._common, p) : getfield(x, p)) |> Base.remove_linenums! |> print

# %%
:(Base.setproperty!(x::$T, p::Symbol, v) =
    p ∈ fieldnames(_Common) ? setfield!(x._common, p, v) : setfield!(x, p, v)) |> Base.remove_linenums! |> print

# %%
:(Base.propertynames(x::$T) = (fieldnames(_Common)...,
    (p for p in fieldnames($T) if first(string(p)) != '_')...)) |> Base.remove_linenums! |> print

# %%
:(function Base.show(io::IO, x::$T)
    props = getproperty.(Ref(x), propertynames(x))
    print(io, nameof($T), '(', repr(first(props)))
    for p in props[2:end] print(io, ", ", repr(p)) end
    print(io, ')')
end) |> Base.remove_linenums! |> print

# %%
module O

using InteractiveUtils: subtypes

abstract type Creature end

Base.@kwdef mutable struct _Common
    hp::Int = 200
    atk::Int = 100
end
for p in fieldnames(_Common)
    @eval $p(x::Creature) = getfield(x._common, $(QuoteNode(p)))
    @eval $(Symbol(:set_, p, :!))(x::Creature, v) = setfield!(x._common, $(QuoteNode(p)), v)
end

Base.@kwdef mutable struct Goblin <: Creature
    _common::_Common = _Common(88, 77)
    weapon::String = "dagger"
end

Base.@kwdef mutable struct Slime <: Creature
    _common::_Common = _Common(55, 44)
    sp::String = "acid"
end

for T in subtypes(Creature)
    n = length(fieldnames(_Common))
    @eval $(nameof(T))(a...) = $T(_Common(a[1:$n]...), a[$n+1:end]...)
    for p in fieldnames(T)
        first(string(p)) == '_' && continue
        @eval $p(x::$T) = getfield(x, $(QuoteNode(p)))
        @eval $(Symbol(:set_, p, :!))(x::$T, v) = setfield!(x, $(QuoteNode(p)), v)
    end
end

end

@show O.Goblin()
@show O.Slime()
println()
@show g = O.Goblin(123, 45, "dirty dagger")
@show O.hp(g)
@show O.atk(g)
@show O.weapon(g)
@show O.set_hp!(g, 99)
@show g
@show O.set_atk!(g, 11)
@show g
@show O.set_weapon!(g, "nothing")
@show g
println()
@show s = O.Slime(123, 45, "weak acid")
@show O.hp(s)
@show O.atk(s)
@show O.sp(s)
@show O.set_hp!(s, 66)
@show s
@show O.set_atk!(s, 99)
@show s
@show O.set_sp!(s, "strong acid")
@show s;

# %%
p = :hp
:($p(x::Creature) = getfield(x._common, $(QuoteNode(p)))) |> Base.remove_linenums! |> print

# %%
:($(Symbol(:set_, p, :!))(x::Creature, v) = setfield!(x._common, $(QuoteNode(p)), v)) |> Base.remove_linenums! |> print

# %%
T = O.Goblin
n = 2
:($(nameof(T))(a...) = $T(_Common(a[1:$n]...), a[$n+1:end]...)) |> Base.remove_linenums! |> print

# %%
:($p(x::$T) = getfield(x, $(QuoteNode(p)))) |> Base.remove_linenums! |> print

# %%
:($(Symbol(:set_, p, :!))(x::$T, v) = setfield!(x, $(QuoteNode(p)), v)) |> Base.remove_linenums! |> print

# %%
p = :hp
:(x = $p) |> Base.remove_linenums! |> print

# %%
p = :hp
:(x = $(QuoteNode(p))) |> Base.remove_linenums! |> print

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
# # Spectrahedron
#
# This example shows an optimization problem over the spectraplex:
# ```math
# S = \{X \in \mathbb{S}_+^n, Tr(X) = 1\}
# ```
# with $\mathbb{S}_+^n$ the set of positive semidefinite matrices.
# Linear optimization with symmetric objective $D$ over the spetraplex consists in computing the
# leading eigenvector of $D$.
#
# The package also exposes `UnitSpectrahedronLMO` which corresponds to the feasible set:
# ```math
# S_u = \{X \in \mathbb{S}_+^n, Tr(X) \leq 1\}
# ```

using FrankWolfe
using LinearAlgebra
using Random
using SparseArrays

# The objective function will be the symmetric squared distance to a set of known or observed entries $Y_{ij}$ of the matrix.
# ```math
# f(X) = \sum_{(i,j) \in L} 1/2 (X_{ij} - Y_{ij})^2
# ```

# ## Setting up the input data, objective, and gradient

# Dimension, number of iterations and number of known entries:
n = 500
k = 10000
n_entries = 50

Random.seed!(41)

const entry_indices = unique!([minmax(rand(1:n, 2)...) for _ in 1:n_entries])
const entry_values = randn(length(entry_indices))

function f(X)
    r = zero(eltype(X))
    for (idx, (i, j)) in enumerate(entry_indices)
        r += 1/2 * (X[i,j] - entry_values[idx])^2
        r += 1/2 * (X[j,i] - entry_values[idx])^2
    end
    return r
end

function grad!(storage, X)
    storage .= 0
    for (idx, (i, j)) in enumerate(entry_indices)
        storage[i,j] += (X[i,j] - entry_values[idx])
        storage[j,i] += (X[j,i] - entry_values[idx])
    end
end

# Note that the `ensure_symmetry = false` argument to `SpectraplexLMO`.
# It skips an additional step making the used direction symmetric.
# It is not necessary when the gradient is a `LinearAlgebra.Symmetric` (or more rarely a `LinearAlgebra.Diagonal` or `LinearAlgebra.UniformScaling`).

const lmo = FrankWolfe.SpectraplexLMO(1.0, n, false)
const x0 = FrankWolfe.compute_extreme_point(lmo, spzeros(n, n))

target_tolerance = 1e-6;

#src the following two lines are used only to precompile the functions
FrankWolfe.frank_wolfe(f, grad!, lmo, x0, max_iteration=2, line_search=FrankWolfe.MonotonousStepSize()) #src
FrankWolfe.lazified_conditional_gradient(f, grad!, lmo, x0, max_iteration=2, line_search=FrankWolfe.MonotonousStepSize()) #src

# ## Running standard and lazified Frank-Wolfe

Xfinal, Vfinal, primal, dual_gap, trajectory = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.MonotonousStepSize(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
    epsilon=target_tolerance,
)

Xfinal, Vfinal, primal, dual_gap, trajectory_lazy = FrankWolfe.lazified_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.MonotonousStepSize(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
    epsilon=target_tolerance,
)

# ## Plotting the resulting trajectories

data = [trajectory, trajectory_lazy]
label = ["FW", "LCG"]
FrankWolfe.plot_trajectories(data, label, xscalelog=true)

ext: jl


text:
"""
DNRM2 returns the euclidean norm of a vector via the function
name, so that

   DNRM2 := sqrt( x'*x )

-- This version written on 25-October-1982.
   Modified on 14-October-1993 to inline the call to DLASSQ.
   Sven Hammarling, Nag Ltd.
"""
function tron_dnrm2(n,x,incx)
    one = 1.0
    zero = 0.0

    if n < 1 || incx < 1
        xnorm = zero
    elseif n == 1
        xnorm = abs(x[1])
    else
        scale = zero
        ssq = one

        # The following loop is equivalent to this call to the LAPACK
        # auxiliary routine:
        # CALL DLASSQ( N, X, INCX, SCALE, SSQ )

        for ix=1:incx:1+(n-1)*incx
            if x[ix] != zero
                absxi = abs(x[ix])
                if scale < absxi
                    ssq = one + ssq*(scale/absxi)^2
                    scale = absxi
                else
                    ssq = ssq + (absxi/scale)^2
                end
            end
        end
        xnorm = scale*sqrt(ssq)
    end

    return xnorm
end

dnrm2(n,x,incx) = tron_dnrm2(n, x, incx)


text:
@defcomp neteconomy begin

	PBACKTIME	= Variable(index=[time])	#Backstop price
	COST1		= Variable(index=[time])	#Adjusted cost for backstop
    ABATECOST   = Variable(index=[time])    #Cost of emissions reductions (trillions 2010 USD per year)
    C           = Variable(index=[time])    #Consumption (trillions 2010 US dollars per year)
    CPC         = Variable(index=[time])    #Per capita consumption (thousands 2010 USD per year)
    CPRICE      = Variable(index=[time])    #Carbon price (2010$ per ton of CO2)
    I           = Variable(index=[time])    #Investment (trillions 2010 USD per year)
    MCABATE     = Variable(index=[time])    #Marginal cost of abatement (2010$ per ton CO2)
    Y           = Variable(index=[time])    #Gross world product net of abatement and damages (trillions 2010 USD per year)
    YNET        = Variable(index=[time])    #Output net of damages equation (trillions 2010 USD per year)

	SIGMA		= Parameter(index=[time])	#CO2-equivalent-emissions output ratio
    DAMAGES     = Parameter(index=[time])   #Damages (Trillion $)
    l           = Parameter(index=[time])   #Level of population and labor
    MIU         = Parameter(index=[time])   #Emission control rate GHGs
    S           = Parameter(index=[time])   #Gross savings rate as fraction of gross world product
    YGROSS      = Parameter(index=[time])   #Gross world product GROSS of abatement and damages (trillions 2010 USD per year)
    expcost2    = Parameter()               #Exponent of control cost function
	pback		= Parameter()				#Cost of backstop 2010$ per tCO2 2015
	gback		= Parameter()				#Initial cost decline backstop cost per period

    function run_timestep(p, v, d, t)
		#Define function for PBACKTIME
        if is_first(t)
            v.PBACKTIME[t] = p.pback
        else
            v.PBACKTIME[t] = v.PBACKTIME[t - 1] * (1 - p.gback)
        end

		#Define function for COSTL
		v.COST1[t] = v.PBACKTIME[t] * p.SIGMA[t] / p.expcost2 / 1000

        #Define function for YNET
        v.YNET[t] = p.YGROSS[t] - p.DAMAGES[t]

        #Define function for ABATECOST
        v.ABATECOST[t] = p.YGROSS[t] * v.COST1[t] * (p.MIU[t]^p.expcost2)

        #Define function for MCABATE (equation from GAMS version)
        v.MCABATE[t] = v.PBACKTIME[t] * p.MIU[t]^(p.expcost2 - 1)

        #Define function for Y
        v.Y[t] = v.YNET[t] - v.ABATECOST[t]

        #Define function for I
        v.I[t] = p.S[t] * v.Y[t]

        #Define function for C
        v.C[t] = v.Y[t] - v.I[t]

        #Define function for CPC
        v.CPC[t] = 1000 * v.C[t] / p.l[t]

        #Define function for CPRICE (equation from GAMS version of DICE2016)
        v.CPRICE[t] = v.PBACKTIME[t] * (p.MIU[t]^(p.expcost2 - 1))
    end
end


text:
using Distributed
@everywhere using Pkg
@everywhere Pkg.activate(@__DIR__)
@everywhere begin
    using CalibrateEDMF
    ce = pkgdir(CalibrateEDMF)
    gs = joinpath(ce, "experiments", "grid_search")
    include(joinpath(gs, "grid_search.jl"))
    include(joinpath(gs, "loss_map.jl"))
    include(joinpath(gs, "quick_plots.jl"))
end
using Test

config_path = joinpath(@__DIR__, "grid_search_integration_test_config.jl")
include(config_path)
config = get_config()  #
out_dir = joinpath(@__DIR__, "output", "grid_search")
rm(out_dir, recursive = true, force = true)  # delete old test output
@time begin
    grid_search(config, config_path, out_dir)
    @info "Grid search complete"
end
y_root = joinpath(out_dir, "general_stochastic_ent_params_{1}.general_stochastic_ent_params_{2}/0.3_0.2")
y_dirs = [
    joinpath(y_root, "TRMM_LBA.1/Output.TRMM_LBA.1"),
    joinpath(y_root, "Rico.1/Output.Rico.1"),
    joinpath(y_root, "TRMM_LBA.2/Output.TRMM_LBA.1"),
]
nc_files = joinpath.(y_dirs, ["stats/Stats.TRMM_LBA.nc", "stats/Stats.Rico.nc", "stats/Stats.TRMM_LBA.nc"])
@test all(isfile.(nc_files))

# update reference data paths
config["reference"]["y_dir"] = y_dirs

@time begin
    compute_loss_map(config, out_dir)
    @info "Loss map complete"
end

quick_plots(joinpath(out_dir, "loss.nc"), out_dir)


text:
function RealVectorFormat()
    return RealVectorFormat(())
end

function RealVectorFormat(arg0::JString, arg1::JString, arg2::JString)
    return RealVectorFormat((JString, JString, JString), arg0, arg1, arg2)
end

function RealVectorFormat(arg0::JString, arg1::JString, arg2::JString, arg3::NumberFormat)
    return RealVectorFormat((JString, JString, JString, NumberFormat), arg0, arg1, arg2, arg3)
end

function RealVectorFormat(arg0::NumberFormat)
    return RealVectorFormat((NumberFormat,), arg0)
end

function format(obj::RealVectorFormat, arg0::RealVector)
    return jcall(obj, "format", JString, (RealVector,), arg0)
end

function format(obj::RealVectorFormat, arg0::RealVector, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (RealVector, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function get_available_locales(::Type{RealVectorFormat})
    return jcall(RealVectorFormat, "getAvailableLocales", Vector{Locale}, ())
end

function get_format(obj::RealVectorFormat)
    return jcall(obj, "getFormat", NumberFormat, ())
end

function get_prefix(obj::RealVectorFormat)
    return jcall(obj, "getPrefix", JString, ())
end

function get_real_vector_format(::Type{RealVectorFormat})
    return jcall(RealVectorFormat, "getRealVectorFormat", RealVectorFormat, ())
end

function get_real_vector_format(::Type{RealVectorFormat}, arg0::Locale)
    return jcall(RealVectorFormat, "getRealVectorFormat", RealVectorFormat, (Locale,), arg0)
end

function get_separator(obj::RealVectorFormat)
    return jcall(obj, "getSeparator", JString, ())
end

function get_suffix(obj::RealVectorFormat)
    return jcall(obj, "getSuffix", JString, ())
end

function parse(obj::RealVectorFormat, arg0::JString)
    return jcall(obj, "parse", ArrayRealVector, (JString,), arg0)
end

function parse(obj::RealVectorFormat, arg0::JString, arg1::ParsePosition)
    return jcall(obj, "parse", ArrayRealVector, (JString, ParsePosition), arg0, arg1)
end


text:
#
include("mesh.jl")

mutable struct Connectivity
    ndof::Int64
    e2dof::Array{Array{Int64,1},1}
    bnddofs::Array{Int64, 1}
end

function connectivity_k0(msh::mesh)
    #= Compute the connectivity table for the cas k=0 =#
    e2edg = msh.edg
    e2nod = msh.elt
    nedg  = msh.nedg
    nelt = length(e2nod)
    conn = [zeros(Int64, length(e2edg[i])) for i=1:nelt]
    edgflagged = zeros(Bool, nedg)
    for e=1:nelt
        edgloc  = e2edg[e]
        nedgloc = length(edgloc)
        for ee=1:nedgloc
            if edgflagged[edgloc[ee]] != true
                conn[e][ee] = edgloc[ee]
                edgflagged[edgloc[ee]] = true
            else
                conn[e][ee] = -edgloc[ee]
            end
        end
    end
    # dofs on the boundary
    nbnddof = sum(msh.freeedg)
    bnddofs = zeros(Int64, nbnddof)
    it = 0
    for e in eachindex(msh.freeedg)
        if msh.freeedg[e]
            it += 1
            bnddofs[it] = e
        end
    end
    #
    return return Connectivity(nedg, conn, bnddofs)
end


text:
# Implement NumberToPattern

const NUCLEOTIDES = ('A', 'C', 'G', 'T')

function number_to_pattern(encoding, k)
    pattern = Char[]
    for _ in 1:k
        encoding, remainder = divrem(encoding, 4)
        push!(pattern, NUCLEOTIDES[remainder+1])
    end
    join(reverse(pattern))
end

function main()
    encoding = parse(Int64, readline())
    k = parse(Int64, readline())
    result = number_to_pattern(encoding, k)
    println(result)
end

main()


text:
export iauD2tf
"""
Decompose days to hours, minutes, seconds, fraction.

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

Status:  vector/matrix support function.

Given:
   ndp     int     resolution (Note 1)
   days    double  interval in days

Returned:
   sign    char    '+' or '-'
   ihmsf   int[4]  hours, minutes, seconds, fraction

Notes:

   1. The argument ndp is interpreted as follows:

   ndp         resolution
   :      ...0000 00 00
   -7         1000 00 00
   -6          100 00 00
   -5           10 00 00
   -4            1 00 00
   -3            0 10 00
   -2            0 01 00
   -1            0 00 10
   0            0 00 01
   1            0 00 00.1
   2            0 00 00.01
   3            0 00 00.001
   :            0 00 00.000...

   2. The largest positive useful value for ndp is determined by the
   size of days, the format of double on the target platform, and
   the risk of overflowing ihmsf[3].  On a typical platform, for
   days up to 1.0, the available floating-point precision might
   correspond to ndp=12.  However, the practical limit is typically
   ndp=9, set by the capacity of a 32-bit int, or ndp=4 if int is
   only 16 bits.

   3. The absolute value of days may exceed 1.0.  In cases where it
   does not, it is up to the caller to test for and handle the
   case where days is very nearly 1.0 and rounds up to 24 hours,
   by testing for ihmsf[0]=24 and setting ihmsf[0-3] to zero.

This revision:  2013 June 18

SOFA release 2018-01-30

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

# void iauD2tf(int ndp, double days, char *sign, int ihmsf[4])

function iauD2tf(ndp::Int, days::Real)

   # Allocate output
   ref_sign = Ref{UInt8}(0)
   ihmsf  = zeros(Int32, 4)

   ccall((:iauD2tf, libsofa_c), Cvoid,
         (Cint, Cdouble, Ref{Cuchar}, Ptr{Cint}),
         convert(Int32, ndp), convert(Float64, days),
         ref_sign, ihmsf)

   return convert(Char, ref_sign[]), SVector{4}(ihmsf)
end


text:
Base.show(io::IO, uri::URI) = print(io, "URI(\"", uri.uri, "\")")
Base.show(io::IO, client::Client) = print(io, "Client(URI(\"", client.uri, "\"))")
Base.getindex(client::Client, database::String) = Database(client, database)

"""
    Client(host, port)
    Client(uri)
    Client()
    Client(pool; [try_pop=false])

Creates a `Client`, which represents a connection to a MongoDB database.

See also [`Mongoc.ClientPool`](@ref).

# Examples:

These lines are equivalent.

```julia
c = Mongoc.Client()
c = Mongoc.Client("localhost", 27017)
c = Mongoc.Client("mongodb://localhost:27017")
```
"""
Client(host::String, port::Int) = Client(URI("mongodb://$host:$port"))
Client(uri::String) = Client(URI(uri))
Client() = Client("localhost", 27017)

"""
    set_appname!(client::Client, appname::String)

Sets the application name for this client.

This string, along with other internal driver details,
is sent to the server as part of the initial connection handshake.

# C API

* [`mongoc_client_set_appname`](http://mongoc.org/libmongoc/current/mongoc_client_set_appname.html).
"""
function set_appname!(client::Client, appname::String)
    ok = mongoc_client_set_appname(client.handle, appname)
    if !ok
        error("Couldn't set appname=$appname for client $client.")
    end
    nothing
end

"""
    ping(client::Client) :: BSON

Pings the server, testing wether it is reachable.

One thing to keep in mind is that operations on MongoDB are *lazy*,
which means that a client reaches a server only when it needs to
transfer documents.

# Example

```julia
julia> client = Mongoc.Client() # nothing happens here between client and server
Client(URI("mongodb://localhost:27017"))

julia> Mongoc.ping(client) # connection to server happens here
BSON("{ "ok" : 1.0 }")
```
"""
function ping(client::Client) :: BSON
    return command_simple(client["admin"], BSON("""{ "ping" : 1 }"""))
end

"""
    get_server_mongodb_version(client::Client) :: VersionNumber

Queries the version for the MongoDB server instance.
"""
function get_server_mongodb_version(client::Client) :: VersionNumber
    bson_server_status = command_simple(client["admin"], BSON("""{ "serverStatus" : 1 }"""))
    return VersionNumber(bson_server_status["version"])
end

"""
    find_databases(client::Client; options::Union{Nothing, BSON}=nothing) :: Cursor

Queries for databases.
"""
function find_databases(client::Client; options::Union{Nothing, BSON}=nothing) :: Cursor
    options_handle = options == nothing ? C_NULL : options.handle
    cursor_handle = mongoc_client_find_databases_with_opts(client.handle, options_handle)
    if cursor_handle == C_NULL
        error("Couldn't execute query.")
    end
    return Cursor(client, cursor_handle)
end

"""
    get_database_names(client::Client; options::Union{Nothing, BSON}=nothing) :: Vector{String}

Helper method to get a list of names for all databases.

See also [`Mongoc.find_databases`](@ref).
"""
function get_database_names(client::Client; options::Union{Nothing, BSON}=nothing) :: Vector{String}
    result = Vector{String}()
    for bson_database in find_databases(client, options=options)
        push!(result, bson_database["name"])
    end
    return result
end

"""
    has_database(client::Client, database_name::String;
        options::Union{Nothing, BSON}=nothing) :: Bool

Helper method to check if there is a database named `database_name`.

See also [`Mongoc.find_databases`](@ref).
"""
function has_database(client::Client, database_name::String;
                      options::Union{Nothing, BSON}=nothing) :: Bool
    for bson_database in find_databases(client, options=options)
        if bson_database["name"] == database_name
            return true
        end
    end
    return false
end


text:
"A non-allocating versin of norm(x::SVector, Inf)"
@inline infnorm(x::StaticVector) = maximum(abs.(x))


text:
function f(x)
  println("Calling f($x)")
end

function f(x, y::Number)
  println("Calling y::Number f($x, $y)")
end

function f(x, y::Array)
  println("Calling y::Array f($x, $y)")
end

function f(x, y::Matrix)
  println("Calling y::Matrix f($x, $y)")
end

function f(x, y::Dict)
  println("Calling y::Dict f($x, $y)")
end

function f(x, y::Vector{Dict{Symbol, Any}})
  println("Calling y::Vector{Dict{Symbol, Any}} f($x, $y)")
end

function f(x, y::Vector{Dict{String, Any}})
  println("Calling y::Vector{Dict{String, Any}} f($x, $y)")
end

f(5)

f(5, 3)

f(5, [2])

f(5, [2.0])

f(5, [1 2;3 4])

dct1 = Dict{Symbol, Any}(:y => 1)
f(5, dct1)

dct2 = Dict{Symbol, Any}[
  Dict(:theta => 0.1),
  Dict(:theta => [0.4]),
  Dict(:theta => [0.1, 0.2, 0.4]),
  Dict(:theta => [0.4 0.5; 0.6 0.7]),
]
f(5, dct2)

dct3 = Dict[
  Dict("theta" => 0.1),
  Dict("theta" => [0.4]),
  Dict("theta" => [0.1, 0.2, 0.4]),
  Dict("theta" => [0.4 0.5; 0.6 0.7]),
]
f(5, dct3)


text:
# First script, mostly exploration
using Pkg
Pkg.activate(".")
using MAT, CSV, Makie, CairoMakie, DataFrames

max_stars_repo_name: fieldofnodes/MLPPProblemsInJulia

max_issues_repo_name: fieldofnodes/MLPPProblemsInJulia

max_forks_repo_name: fieldofnodes/MLPPProblemsInJulia


text:
module _Sqrt

using Brainstorm._Math._GMP: set!, add!, mul!, neg!, lsh!, rsh!

const bits_per_digit = log2(10)
const guard_rounds = 1

#
# Compute int(sqrt(x) * 10^digits)
#
# This is done using Newton-Raphson method.
#
sqrt(digits::Integer, x::Integer) = begin
    r, e = _reciprocalsqrt(digits, x)

    # (r / 2^e) * n * 10^digits
    s = big(5)^digits
    mul!(r, x)
    mul!(r, s)
    rsh!(r, e - digits)

    r
end

rsqrt(digits::Integer, x::Integer) = begin
    r, e = _reciprocalsqrt(digits, x)

    # (r / 2^e) * 10^digits
    s = big(5)^digits
    mul!(r, s)
    rsh!(r, e - digits)

    r
end

_reciprocalsqrt(digits, x) = begin
    # initial guess is accurate to eᵢ bits
    eᵢ = Int64(50)
    rᵢ = trunc(BigInt, (Int64(1) << eᵢ) / Base.sqrt(x))

    target_precision = ceil(Int64, digits * bits_per_digit) << guard_rounds

    oᵢ, wᵢ, dᵢ = big(0), big(0), big(0)
    while eᵢ <= target_precision
        peᵢ = eᵢ
        eᵢ = eᵢ << 1

        #
        # http://www.numberworld.org/y-cruncher/algorithms/invsqrt.html
        # http://cs.stackexchange.com/a/37645
        #
        #      w_i+1 = r_i^2
        #            = r_i'^2 / 2^(2*e_i)
        #  => w_i+1' = r_i'^2
        #
        #      d_i+1 = 1 - w_i+1 * x
        #            = 1 - r_i'^2 / 2^(2*e_i) * x
        #            = (2^(2*e_i) - r_i'^2 * x) / 2^(2*e_i)
        #  => d_i+1' = 2^(2*e_i) - w_i+1' * x
        #
        #      r_i+1 = r_i * (1 + d_i+1 / 2)
        #            = r_i' / 2^e_i * (1 + (d_i+1' / 2) / 2^(2*e_i))
        #            = r_i' / 2^e_i * (2^(2*e_i) + d_i+1' / 2) / 2^(2*e_i)
        #            = r_i' * (2^(2*e_i) + d_i+1' / 2) / 2^e_i / 2^(2*e_i)
        #  => r_i+1' = r_i' * (2^(2*e_i) + d_i+1' / 2) / 2^e_i
        #

        # oᵢ = 1 << eᵢ
        set!(oᵢ, 1)
        lsh!(oᵢ, eᵢ)

        # wᵢ = rᵢ * rᵢ
        set!(wᵢ, rᵢ)
        mul!(wᵢ, rᵢ)
        # dᵢ = oᵢ - wᵢ * x
        set!(dᵢ, wᵢ)
        mul!(dᵢ, x)
        neg!(dᵢ)
        add!(dᵢ, oᵢ)
        # rᵢ = (rᵢ * (oᵢ + dᵢ >> 1)) >> peᵢ
        rsh!(dᵢ, 1)
        add!(dᵢ, oᵢ)
        mul!(rᵢ, dᵢ)
        rsh!(rᵢ, peᵢ)
    end

    rᵢ, eᵢ
end

end


text:
#=
        project1.jl -- This is where the magic happens!

    All of your code must either live in this file, or be `include`d here.
=#

#=
    If you want to use packages, please do so up here.
    Note that you may use any packages in the julia standard library
    (i.e. ones that ship with the julia language) as well as Statistics
    (since we use it in the backend already anyway)
=#

# Example:
# using LinearAlgebra

#=
    If you're going to include files, please do so up here. Note that they
    must be saved in project1_jl and you must use the relative path
    (not the absolute path) of the file in the include statement.

    [Good]  include("somefile.jl")
    [Bad]   include("/pathto/project1_jl/somefile.jl")
=#

# Example
# include("myfile.jl")
include("algo_zeroth_order.jl")
include("algo_first_order.jl")
include("algo_util.jl")

"""
    optimize(f, g, x0, n, prob)

Arguments:
    - `f`: Function to be optimized
    - `∇f`: Gradient function for `f`
    - `x0`: (Vector) Initial position to start from
    - `n`: (Int) Number of evaluations allowed. Remember `g` costs twice of `f`
    - `prob`: (String) Name of the problem. So you can use a different strategy for each problem. E.g. "simple1", "secret2", etc.

Returns:
    - The location of the minimum
"""
function optimize(f, ∇f, x0, n, prob_name)

    if prob_name == "simple1"
        method = HookeJeevesDynamic(α=0.3, γ=0.3)
        x, _ = solve(method, f, x0, n)

    elseif prob_name == "simple2"
        method = HookeJeevesDynamic(α=0.3)
        x, _ = solve(method, f, x0, n)

    elseif prob_name == "simple3"
        method = HookeJeevesDynamic(α=0.3)
        x, _ = solve(method, f, x0, n)

    elseif prob_name == "secret1"
        method = HookeJeevesDynamic(α=0.5, γ=0.5)
        x, _ = solve(method, f, x0, n)

    elseif prob_name == "secret2"
        method = HookeJeevesDynamic(α=1.0, γ = 0.3)
        x, _ = solve(method, f, x0, n)

    else
        # method = GradientDescent(3e-4)
        # method = GDMomentum(α=3e-4, β=0.95)
        # method = GDMomentumNesterov(α=3e-4, β=0.6)
        method = Adam(α=3e-4)
        # method = GDApproxLineSearch(α=3e-4)

        x, _ = solve(method, f, ∇f, x0, n)

        # method = HookeJeeves(α=1.0)
        # x, _ = solve(method, f, x0, n)

        # method = HookeJeevesDynamic(α=0.3)
        # x, _ = solve(method, f, x0, n)

        # method = CoordinateDescentAcceleration()
        # x, x_hist = solve(method, f, ∇f, x0, 5; num_eval_termination=false)
    end
        # println("Pre: $(round.(x0; digits=3)) -> $(round(f(x0); digits=3))")
        # println("Pos: $(round.(x; digits=3)) -> $(round(f(x); digits=3)) \n")
    return x
end


text:
@enum frequencyUnit begin
   Hz = 1
   kHz = 2
   MHz = 3
   GHz = 4
end

@enum parameterType begin
    S = 1 # Scattering parameters
    Y = 2 # Admittance parameters
    Z = 3 # Impedance parameters
    H = 4 # Hybrid-h parameters
    G = 5 # Hybrid-g parameters
end
@enum parameterFormat begin
    DB = 1 # decibel angle (decibel = 20log10|magnitude|) [deg]
    MA = 2 # magnitude angle [deg]
    RI = 3 #real-imaginary
end

function skip_comments(sentences::Array, i::Integer)
    while match(r"^(\s*)$", sentences[i]) != nothing || match(r"^(\s*\!)", sentences[i]) != nothing || sentences[i] == ""
        i += 1
    end
    return i
end

function transform_vector_to_matrix_mode(arr::Array, matrix_mode::String, freq_n::Int)
    freq_params = Dict() # separate matrices based on frequencies
    i = 1
    k = Int(length(arr)/freq_n)
    for m=1:freq_n
        varr = arr[(k*(m-1) + 2):k*m]
        n = Int(round(0.5*(sqrt(4*k - 3) - 1)))
        # TODO: Use LinearAlgebra to construct the diagonal matrices
        if matrix_mode == "UPPER"
            n = round(0.5*(sqrt(4*k - 3) - 1))
            freq_params[arr[k*(m-1) + 1]] = [2i <= j ? (arr[(n*(n+1)-(n-i+1+1)*(n-i+1))>>>1 + 1 + (j-1)%n - (i - 1)]) : 0 for i=1:2n, j=1:n]

        elseif matrix_mode == "LOWER"
            n = Int(round(0.5*(sqrt(4*k - 3) - 1)))
            freq_params[arr[k*(m-1) + 1]] = [2i >= j ? varr[(((i+1)*i)>>>1 + (n-1) - ((n - j )%n)- (i-1))] : 0 for i=1:n, j=1:2n]
        elseif matrix_mode == "FULL"
            n = Int(round(sqrt((k - 1)/2)))
            freq_params[arr[k*(m-1) + 1]] = [varr[(i-1)*2n + (j-1)%2n + 1] for i=1:n, j=1:2n]
        end
        i += 1
    end
    return freq_params
 end

 function scatter_params_to_admittance(s_params::Array, ref_resistance::Int)
    n, m = size(s_params)
    if n != m
        throw("Provided matrix should be square")
    end

    return inv(I + s_params)*(I - s_params)/ref_resistance
 end

 function admittance_params_to_scatter(y_params::Array, ref_resistance::Int)
    n, m = size(y_params)
    if n != m
        throw("Provided matrix should be square")
    end

    return inv(I + y_params*ref_resistance)*(I - y_params*ref_resistance)
 end

 function impedance_params_to_scatter(z_params::Array, ref_resistance::Int)
    n, m = size(z_params)
    if n != m
        throw("Provided matrix should be square")
    end

    return inv(I + z_params/ref_resistance)*(z_params/ref_resistance - I)
 end

 function scatter_params_to_impedance(s_params::Array, ref_resistance::Int)
    n, m = size(s_params)
    if n != m
        throw("Provided matrix should be square")
    end

    return inv(I - s_params)*(I + s_params)*ref_resistance
 end

 function to_complex(carr)
    return carr[1] + carr[2]im
 end

 """
 Converts parameters in Magnitude-Angle to Real-Imaginary format
...
# Arguments
- `ma_params` : Parameters in the form (magnitude[linear], angle[rad])
...
"""
function ma_to_ri(ma_params::Array)
    return [ma_params[1]*cos(ma_params[2]), ma_params[1]*sin(ma_params[2])]
end

"""
Converts parameters in Real-Imaginary to Magnitude-Angle format
...
# Arguments
- `ri_params` : Parameters in the form (real, imaginary)
...
"""
function ri_to_ma(ri_params::Array)
   return [sqrt(ri_params[1]*ri_params[1] + ri_params[2]*ri_params[2]), atan(ri_params[2]/ri_params[1])]
end

"""
Converts parameters in Decibel-Angle to Magnitude-Angle format
...
# Arguments
- `db_params` : Parameters in the form (magnitude[db], angle[rad])

Note: for Touchstoner ver2.0 files decibel = 20log10
...
"""
function db_to_ma(db_params::Array)
   return [10^(db_params[1]/20), db_params[2]]
end

"""
Converts parameters in Magnitude-Angle to Decibel-Angle format
...
# Arguments
- `ma_params` : Parameters in the form (magnitude[linear], angle[rad])

Note: for Touchstoner ver2.0 files decibel = 20log10
...
"""
function ma_to_db(ma_params::Array)
   return [20log10(ma_params[1]), ma_params[2]]
end


text:
module Shapes

include("PartialSphere.jl")
include("Cylinder.jl")
include("Grid.jl")
include("PartialTorus.jl")

function choose(shape::Symbol)
    apps = Dict(
            :PartialSphere =>PartialSphere,
            :Cylinder => Cylinder,
            :Grid => Grid,
            :PartialTorus => PartialTorus,
        )
    if shape in keys(apps)
        return apps[shape].assemble
    else
        error("Error hint: shape \"$(string(shape))\" not recognized")
    end
end

end

max_stars_repo_name: yuhangwang/SimShape

max_issues_repo_name: yuhangwang/SimShape

max_forks_repo_name: yuhangwang/SimShape


text:
using BinaryBuilder

name = "Blosc"
version = v"1.14.3"

# Collection of sources required to build Blosc
sources = [
    ArchiveSource("https://github.com/Blosc/c-blosc/archive/v$(version).tar.gz", "7217659d8ef383999d90207a98c9a2555f7b46e10fa7d21ab5a1f92c861d18f7"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/c-blosc-*
mkdir build
cd build
CMAKE_FLAGS=(-DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release)
CMAKE_FLAGS+=(-DBUILD_TESTS=OFF -DBUILD_BENCHMARKS=OFF)
CMAKE_FLAGS+=(-DBUILD_STATIC=OFF)
CMAKE_FLAGS+=(-DCMAKE_SHARED_LIBRARY_LINK_C_FLAGS="")
CMAKE_FLAGS+=(-DCMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS="")
CMAKE_FLAGS+=(-DPREFER_EXTERNAL_ZLIB=ON)
CMAKE_FLAGS+=(-DPREFER_EXTERNAL_ZSTD=ON)
CMAKE_FLAGS+=(-DPREFER_EXTERNAL_LZ4=ON)
#CMAKE_FLAGS+=(-DPREFER_EXTERNAL_SNAPPY=ON)
cmake ${CMAKE_FLAGS[@]} ..
make -j${nproc}
make install

if [[ "${target}" == *-mingw* ]]; then
    # Manually move dlls from lib/ to bin/
    mv ${prefix}/lib/*.${dlext} ${libdir}
fi

install_license ../LICENSES/*.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libblosc", :libblosc)
]

# Dependencies that must be installed before this package can be built
dependencies = [
	Dependency("Zlib_jll"),
	Dependency("Zstd_jll"),
	Dependency("Lz4_jll"),
	# Dependency("Snappy_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
type VarStats
    n_obs::Int64
    min::Float64
    max::Float64
    var::Float64
    sum::Float64

    VarStats() = new(0, +Inf, -Inf, 0, 0)
end

count(vs::VarStats) = vs.n_obs
sum(vs::VarStats) = vs.sum
minimum(vs::VarStats) = vs.min
maximum(vs::VarStats) = vs.max
mean(vs::VarStats) = vs.sum / vs.n_obs
var(vs::VarStats) = vs.var
std(vs::VarStats) = sqrt(vs.var)

function updatestats!(vs::VarStats, value::Real)
    vs.n_obs += 1
    vs.sum += value
    if value < vs.min
        vs.min = value
    end
    if value > vs.max
        vs.max = value
    end
    if vs.n_obs > 2
        n = vs.n_obs
        vs.var = (n - 2) / (n - 1) * vs.n_obs + 1 / n * (value - mean(vs))^2
    end
end

function updatestats!{T<:Real}(varstats::Vector{VarStats}, values::Vector{T})
    for i=1:length(varstats)
        updatestats!(varstats[i], values[i])
    end
end


text:
# Chipmunk's axis-aligned 2D bounding box type along with a few handy routines.

# Constructs a cpBB centered on a point with the given extents (half sizes).
function BB_for_extents(c::Vect, hw::Real, hh::Real)
    BB(ccall(dlsym(libchipmunk, :cpBBNewForExtents), Ptr{Void}, (Vect, Cdouble, Cdouble,), c, hw, hh))
end

# Constructs a cpBB for a circle with the given position and radius.
function BB_for_circle(p::Vect, r::Real)
    BB(ccall(dlsym(libchipmunk, :cpBBNewForCircle), Ptr{Void}, (Vect, Cdouble,), p, r))
end

# Returns true if a and b intersect.
function intersects(a::BB, b::BB)
    a.l <= b.r && b.l <= a.r && a.b <= b.t && b.b <= a.t
end

# Returns true if other lies completely within bb.
function containsBB(bb::BB, other::BB)
    bb.l <= other.l && bb.r >= other.r && bb.b <= other.b && bb.t >= other.t
end

# Returns true if bb contains v.
function contains_vect(bb::BB, v::Vect)
    bb.l <= v.x && bb.r >= v.x && bb.b <= v.y && bb.t >= v.y
end

# Returns a bounding box that holds both bounding boxes.
function merge(a::BB, b::BB)
    BB(min(a.l, b.l), min(a.b, b.b), max(a.r, b.r), max(a.t, b.t))
end

# Returns a bounding box that holds both bb and v.
function expand(bb::BB, v::Vect)
    BB(min(bb.l, v.x), min(bb.b, v.y), max(bb.r, v.x), max(bb.t, v.y))
end

# Returns the center of a bounding box.
function center(bb::BB)
    ccall(dlsym(libchipmunk, :cpBBCenter), Vect, (BB,), bb)
end

# Returns the area of the bounding box.
function area(bb::BB)
    (bb.r - bb.l)*(bb.t - bb.b)
end

# Merges a and b and returns the area of the merged bounding box.
function merged_area(a::BB, b::BB)
    (max(a.r, b.r) - min(a.l, b.l))*(max(a.t, b.t) - min(a.b, b.b))
end

# Returns the fraction along the segment query the cpBB is hit. Returns INFINITY if it doesn't hit.
function segment_query(bb::BB, a::Vect, b::Vect)
    ccall(dlsym(libchipmunk, :cpBBQuery), Cdouble, (BB, Vect, Vect,), bb, a, b)
end

# Return true if the bounding box intersects the line segment with ends a and b.
function intersects_segment(bb::BB, a::Vect, b::Vect)
    ccall(dlsym(libchipmunk, :cpBBIntersectsSegment), Bool, (BB, Vect, Vect,), bb, a, b)
end

# Clamp a vector to a bounding box.
function clamp_vect(bb::BB, v::Vect)
    Vect(clamp(v.x, bb.l, bb.r), clamp(v.y, bb.b, bb.t))
end

# Wrap a vector to a bounding box.
function wrap_vect(bb::BB, v::Vect)
    ccall(dlsym(libchipmunk, :cpBBWrapVect), Vect, (BB, Vect,), bb, v)
end

# Returns a bounding box offseted by v.
function offset(bb::BB, v::Vect)
    ccall(dlsym(libchipmunk, :cpBBOffset), Vect, (BB, Vect,), bb, v)
end

export BB_for_extents, BB_for_circle, intersects, containsBB, contains_vect, merge, expand,
center, area, merged_area, segment_query, intersects_segment, clamp_vect, wrap_vect, offset


text:
# Utility functions for motif finding algorithms

DNA = Dict("A" => 1, "T" => 2, "C" => 3, "G" => 4)

function align(sequences, s, l)
    """
    Return the length(*s*) x *l* alignment matrix for *sequences*
    """
    align = []
    for i = 1:length(s)
        # add the motif substring to align array
        push!(align, sequences[i][s[i]:s[i]+l-1])
    end
    return align
end

function getprofile(alignment)
    """
    Return the profile matrix of an alignment
    """
    profile = zeros(Int32, 4, length(alignment[1]))
    for i = 1:length(alignment)
        for j = 1:length(alignment[1])
            # get the character at position [i,j] in alignment
            nucl = alignment[i][[j]]
            profile[DNA[nucl], j] += 1
        end
    end
    return profile
end

score(profile) = sum(maximum(profile, 1))

function score(sequences, s, l)
    """
    Find the alignment and calculate the score of the corresponding profile
    matrix.
    """
    alignment = align(sequences, s, l)
    profile = getprofile(alignment)
    return score(profile)
end

function probscore(motif, profile)
    """
    Calculate the probability of seeing *motif* given the matrix *profile*.
    Return a probability in [0, 1]
    """

    # normalize profile
    P = profile ./ sum(profile, 1)

    return prod(P[DNA[motif[[i]]], i] for i = 1:length(motif))
end

max_stars_repo_name: dillondaudert/JuliaMotifs

max_issues_repo_name: dillondaudert/JuliaMotifs

max_forks_repo_name: dillondaudert/JuliaMotifs


text:
"""
	set_interaction!(itx, args...)

Set interaction parameters for the whole interaction or specific monomers.
"""
set_interaction!(itx::AbstractInteraction, args...) = nothing
set_interaction!(itx::AbstractInteraction, mon::Monomer, args...) = set_interaction!(itx, mon.id, args...)
set_interaction!(itx::AbstractInteraction, mon1::Monomer, mon2::Monomer, args...) = set_interaction!(itx, mon1.id, mon2.id, args...)

"""
    energy(itx)

Return the energy of the `AbstractInteraction` due to the current system configuration.
"""
energy(itx::AbstractInteraction) = nothing

"""
    bulk_energy(itx)

Return the energy of the `AbstractInteraction` due to the homogeneous bulk composition.
"""
bulk_energy(itx::AbstractInteraction) = nothing

"""
    potential!(itx, alpha, potential)

Add the potential contribution for monomer ID `alpha` to the grid.
"""
potential!(itx::AbstractInteraction, alpha::Integer, pot::FieldGrid) = nothing

setup!(itx::AbstractInteraction, sys::AbstractSystem)  = nothing


text:
module index.reserved
using base: @kwdf
using pest
using pest.parse

function pest()
    arr keywords    = ["compute","compute_euler", "compute_mcarlo", "compute_algol", "compute_circle"]
    arr methods     = ["sequences", "verbose", "silent"]
    arr datatype    = [(bigfloat), (bigint)]
    char comment    = "//"
    char mcomment   = "#|  |#"
    char jlmcomment = "/\###"
    char jlcomment  = '/\#'
    comment  = jlcomment
    mcomment = jlmcomment
end

###
# sample twister program:
#
 function compute( n )
   @compute_mcarlo( n )
       sequence = n
   end
 end

 compute( 50 )
#
# prints the monte-carlo algorithm
# for the first 50 sequences
#
###

max_stars_repo_name: timo-cmd/DANCE

max_issues_repo_name: timo-cmd/Twister

max_forks_repo_name: timo-cmd/DANCE


text:
module MoodleQuestions

include("read_write.jl")

# Read from SWAD file
export read_swad
# Read from TXT file
export read_txt
# Save to Moodle
export save_to_moodle

# Export serve version
include("serve.jl")
export serve_quiz

end # module

max_stars_repo_name: dmolina/MoodleQuestions

max_issues_repo_name: dmolina/MoodleQuestions

max_forks_repo_name: dmolina/MoodleQuestions


text:
(rd,wr) = Base.redirect_stdin()
write(wr,"println(Base.roottask === Base.current_task());\n")
write(wr,"@assert false\n")

using Test
using RootTaskREPL

@testset "check_root_task" begin
        @test Base.roottask === Base.current_task()
end

close(wr)


text:
mutable struct TestDev <: AbstractDaqDevice
    "Handling of background data acquisition"
    task::DAQTask
    nchans::Int

end


text:
"""
    $(TYPEDEF)
Kerr spacetime with spin `spin`. We assume that we are in unitless scales where
the black hole mass defines the length, so `M=1`.
"""
struct Kerr{S}
    spin::S
end

"""
    $(TYPEDEF)
For a spacetime `g` construct metric tensor at position `p`.
There is an optional cache that can be used to cache intermediate results for increased
speed.
"""

@inline function _Δ(g::Kerr, r)
    return r^2 - 2*r + g.spin^2
end

@inline function _Σ(g::Kerr, r, θ)
    return r^2 + (g.spin*cos(θ))^2
end

@inline function _Ξ(g::Kerr, r, θ, Δ)
    return (r^2 + g.spin^2)^2 - Δ*(g.spin*sin(θ))^2
end

@inline function _ωZamo(g::Kerr, r, Ξ)
    return 2*g.spin*r/Ξ
end
struct MetricTensor{M, P, C}
    g::M
    p::P
    cache::C
end

@inline function _ℛ(g, r, λ, η, Δ)
    a = g.spin
    return (r^2 + a^2 - a*λ)^2 - Δ*(η + (a-λ)^2)
end

"""
    $(SIGNATURES)
Construct the metric tensor for spacetime `g` at the location `t`, `r`, `θ`, `ϕ`.

Returns a `MetricTensor` object. To access the components of the tensor please
the components function.
"""
function metric(g::Kerr, t, r, θ, ϕ)
    Δ = _Δ(g, r)
    Σ = _Σ(g, r, θ)
    Ξ = _Ξ(g, r, θ,  Δ)
    Ξsin2 = Ξ*sin(θ)^2
    ω = _ωZamo(g, r, Ξ,)
    cache = (;Δ, Σ, Ξ, Ξsin2, ω)
    return MetricTensor(g, (t,r,θ,ϕ), cache)
end

"""
    $(SIGNATURES)
Returns the components of the metric tensor `g` in SMatrix.
"""
function components(g::MetricTensor{K, C}) where {K<:Kerr, C}
    Δ, Σ, Ξ, Ξsin2, ω = g.cache
    gtt = -Δ*Σ/Ξ + ω^2*Ξsin2/Σ
    gtϕ = -ω*Ξsin2/Σ
    grr = Σ/Δ
    gθθ = Σ
    gϕϕ = Ξsin2/Σ
    T = eltype(gtt)
    return @SMatrix [ gtt     zero(T) zero(T) gtϕ    ;
                      zero(T) grr     zero(T) zero(T);
                      zero(T) zero(T) gθθ     zero(T);
                      gtϕ     zero(T) zero(T) gϕϕ
                    ]
end

"""
    $(SIGNATURES)
Computes the radial potential of the metric tensor assuming the energy scaled angular momentum
`λ`` and carter constant `η``
"""
function ℛ(g::MetricTensor, λ, η)
    r = g.p[2]
    Δ = g.cache.Δ
    return _ℛ(g.g, r, λ, η, Δ)
end

struct Minkowski{T} end
Minkowski() = Minkowski{Float64}()
metric(m::Minkowski, t, r, θ, ϕ) = MetricTensor(m, (t,r,θ,ϕ), (;))
components(::MetricTensor{<:Minkowski{T}, P,C}) where {T,P,C} = Diagonal(SVector(-one(T), one(T), one(T), one(T)))


text:
"""
    RationalQuadraticKernel(; α = 2.0)

The rational-quadratic kernel is a Mercer kernel given by the formula:
```
    κ(x,y)=(1+||x−y||²/α)^(-α)
```
where `α` is a shape parameter of the Euclidean distance. Check [`GammaRationalQuadraticKernel`](@ref) for a generalization.
"""
struct RationalQuadraticKernel{Tα<:Real} <: SimpleKernel
    α::Vector{Tα}
    function RationalQuadraticKernel(;alpha::T=2.0, α::T=alpha) where {T}
        @check_args(RationalQuadraticKernel, α, α > zero(T), "α > 1")
        return new{T}([α])
    end
end

kappa(κ::RationalQuadraticKernel, d²::T) where {T<:Real} = (one(T)+d²/first(κ.α))^(-first(κ.α))
metric(::RationalQuadraticKernel) = SqEuclidean()

Base.show(io::IO, κ::RationalQuadraticKernel) = print(io, "Rational Quadratic Kernel (α = ", first(κ.α), ")")

"""
`GammaRationalQuadraticKernel([ρ=1.0[,α=2.0[,γ=2.0]]])`
The Gamma-rational-quadratic kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y)=(1+ρ^(2γ)||x−y||^(2γ)/α)^(-α)
```
where `α` is a shape parameter of the Euclidean distance and `γ` is another shape parameter.
"""
struct GammaRationalQuadraticKernel{Tα<:Real, Tγ<:Real} <: SimpleKernel
    α::Vector{Tα}
    γ::Vector{Tγ}
    function GammaRationalQuadraticKernel(;alpha::Tα=2.0, gamma::Tγ=2.0, α::Tα=alpha, γ::Tγ=gamma) where {Tα<:Real, Tγ<:Real}
        @check_args(GammaRationalQuadraticKernel, α, α > one(Tα), "α > 1")
        @check_args(GammaRationalQuadraticKernel, γ, γ >= one(Tγ), "γ >= 1")
        return new{Tα, Tγ}([α], [γ])
    end
end

kappa(κ::GammaRationalQuadraticKernel, d²::T) where {T<:Real} = (one(T)+d²^first(κ.γ)/first(κ.α))^(-first(κ.α))
metric(::GammaRationalQuadraticKernel) = SqEuclidean()

Base.show(io::IO, κ::GammaRationalQuadraticKernel) = print(io, "Gamma Rational Quadratic Kernel (α = ", first(κ.α), ", γ = ", first(κ.γ), ")")


text:
abstract type SimpleSDMSource end

latitudes(::Type{T}) where {T <: SimpleSDMSource} = (-90.0, 90.0)
longitudes(::Type{T}) where {T <: SimpleSDMSource} = (-180.0, 180.0)

struct WorldClim <: SimpleSDMSource
    resolution::AbstractFloat
    function WorldClim(resolution::AbstractFloat)
        resolution ∈ [2.5, 5.0, 10.0] || throw(ArgumentError("The resolution argument ($(resolution)) must be 2.5, 5, or 10"))
        return new(resolution)
    end
end

WorldClim() = WorldClim(10.0)

struct BioClim <: SimpleSDMSource end
longitudes(::Type{BioClim}) = (-180.0001388888, 179.9998611111)
latitudes(::Type{BioClim}) = (-90.0001388888, 83.9998611111)

struct EarthEnv <: SimpleSDMSource
    full::Bool
end

EarthEnv() = EarthEnv(false)
latitudes(::Type{EarthEnv}) = (-56.0, 90.0)
longitudes(::Type{EarthEnv}) = (-180.0, 180.0)


text:
module DistributedJets

using Distributed, DistributedArrays, DistributedOperations, JSON, Jets, LinearAlgebra, Statistics

#
# DArray extensions
#
function DistributedArrays.DArray(init::Function, _pids::AbstractArray{Int}, rngs::Vector{UnitRange{Int}}...)
    pids = vec(_pids)
    np = mapreduce(length, *, rngs)
    @assert length(pids) >= np
    pids = reshape(pids[1:np], map(rng->length(rng), rngs))
    idxs,cuts = DistributedArrays.chunk_idxs(rngs...)
    dims = ntuple(i->idxs[end][i][end], length(rngs))
    id = DistributedArrays.next_did()
    DistributedArrays.DArray(id, init, dims, pids, idxs, cuts)
end

function DistributedArrays.chunk_idxs(rngs::Vector{UnitRange{Int}}...)
    chunks = map(rng->length(rng), rngs)
    cuts = [prescribedist(rng) for rng in rngs]
    ndim = length(rngs)
    idxs = Array{NTuple{ndim,UnitRange{Int}}}(undef, chunks...)
    for cidx in CartesianIndices(chunks)
        idxs[cidx] = ntuple(i->(cuts[i][cidx[i]]:cuts[i][cidx[i]+1]-1), ndim)
    end
    idxs, cuts
end

function prescribedist(rng::Vector{UnitRange{Int}})
    n = length(rng)
    cuts = Vector{Int}(undef, n+1)
    for j=1:n
        cuts[j] = rng[j][1]
    end
    cuts[n+1] = rng[n][end]+1
    cuts
end

Jets.indices(cuts::AbstractVector{Int}) = [cuts[i]:(cuts[i+1]-1) for i=1:(length(cuts)-1)]
Jets.indices(A::DArray{T}) where {T} = A.indices
Jets.indices(A::DArray{T,1}) where {T} = indices(A.cuts[1])
Jets.indices(A::DArray{T}, i::Integer) where {T} = indices(A.cuts[i])

#
# Distributed block operators
#
struct JetDSpace{T,S<:Jets.JetAbstractSpace{T,1}} <: Jets.JetAbstractSpace{T,1}
    blkspaces::DArray{S,1,Array{S,1}}
    blkindices::Vector{UnitRange{Int}}
    indices::Vector{UnitRange{Int}}
end

function JetDSpace(blkspaces::DArray{S,1}) where {S<:JetBSpace}
    blkₒ,iₒ = 1,1
    blkindices = Vector{UnitRange{Int}}(undef, length(blkspaces))
    indices = Vector{UnitRange{Int}}(undef, length(blkspaces))
    for i = 1:length(blkspaces)
        blkspacesᵢ = blkspaces[i]
        blkindices[i] = blkₒ:(blkₒ+nblocks(blkspacesᵢ)-1)
        indices[i] = iₒ:(iₒ+length(blkspacesᵢ)-1)
        blkₒ = blkindices[i][end] + 1
        iₒ = indices[i][end] + 1
    end

    JetDSpace(blkspaces, blkindices, indices)
end

Base.:(==)(R::JetDSpace, S::JetDSpace) = R.blkspaces == S.blkspaces && R.blkindices == S.blkindices && R.indices == S.indices

"""
    size(R)
returns the size of R
"""
Base.size(R::JetDSpace) = (R.indices[end][end],)
"""
    eltype(R)
returns the element type of R
"""
Base.eltype(R::Type{JetDSpace{T,S}}) where {T,S} = T
Base.eltype(R::Type{JetDSpace{T}}) where {T} = T
Base.vec(R::JetDSpace) = R

Jets.indices(R::JetDSpace) = R.indices
Jets.space(R::JetDSpace, iblock::Integer) where {T,S} = R.blkspaces[iblock]
Jets.nblocks(R::JetDSpace) = R.blkindices[end][end]

DistributedArrays.localindices(R::JetDSpace) = R.indices[findfirst(pid->pid==myid(), procs(R))]
DistributedArrays.localpart(R::JetDSpace) = localpart(R.blkspaces)[1]
"""
    localblockindices(R)
returns the block indicies that are local to the calling process
"""
localblockindices(R::JetDSpace) = R.blkindices[findfirst(pid->pid==myid(), procs(R))]
localblockindices(R::Jets.JetBSpace) = 1:length(R.indices)
localblockindices(R::Jets.JetAbstractSpace) = 1:1
"""
    blockmap(A)
returns map between process id and block indicies
"""
blockmap(R::JetDSpace) = R.blkindices

Distributed.procs(R::JetDSpace) = procs(R.blkspaces)
Distributed.nprocs(R::JetDSpace) = length(procs(R.blkspaces))

struct DBArray{T,A<:Jets.BlockArray{T},B} <: AbstractArray{T,1}
    darray::DArray{A,1,B}
    indices::Vector{UnitRange{Int}}
    blkindices::Vector{UnitRange{Int}}
end

function DBArray_localpart(f, blkidxs)
    arrays = [f(blkidx) for blkidx in blkidxs]
    indices = Vector{UnitRange{Int}}(undef, length(blkidxs))

    ibeg = 1
    for i = 1:length(arrays)
        iend = ibeg + length(arrays[i]) - 1
        indices[i] = ibeg:iend
        ibeg = iend + 1
    end

    Jets.BlockArray(arrays, indices)
end

function DBArray_local_length(future)
    x = fetch(future)
    isa(x, RemoteException) && throw(x)
    length(x)
end

function DBArray(f::Function, blkidxs::Vector{UnitRange{Int}}, pids=workers())
    futures = Dict{Int,Future}()
    for (ipid,pid) in enumerate(pids)
        futures[pid] = remotecall(DBArray_localpart, pid, f, blkidxs[ipid])
    end

    idxs = Vector{UnitRange{Int}}(undef, length(pids))
    ibeg = 1
    for (ipid,pid) in enumerate(pids)
        iend = ibeg + remotecall_fetch(DBArray_local_length, pid, futures[pid]) - 1
        idxs[ipid] = ibeg:iend
        ibeg = iend + 1
    end

    darray = DArray(i->[fetch(futures[myid()])], pids, idxs)
    DBArray(darray, idxs, blkidxs)
end

function DBArray(f::Function, nblks::Tuple, pids::AbstractArray, dist::AbstractArray)
    blkidxs, cuts = DistributedArrays.chunk_idxs([nblks...], dist) # cuts is not used
    DBArray(f, [blkidx[1] for blkidx in blkidxs], pids)
end

DBArray(f::Function, nblks::Tuple, pids::AbstractArray) = DBArray(f, nblks, pids, DistributedArrays.defaultdist(nblks, pids))
DBArray(f::Function, nblks::Tuple) = DBArray(f, nblks, workers()[1:min(nworkers(),nblks[1])], DistributedArrays.defaultdist(nblks, workers()[1:min(nworkers(),nblks[1])]))

function Jets.space(x::DBArray)
    blkspaces = DArray(I->[space(localpart(x)) for i in I[1]], (nprocs(x),), procs(x))
    JetDSpace(blkspaces)
end

# DBArray array interface implementation <--
Base.IndexStyle(::Type{T}) where {T<:DBArray} = IndexLinear()
Base.size(x::DBArray) = (x.indices[end][end],)
Jets.indices(x::DBArray) = indices(x.darray)
Jets.indices(x::DBArray, i::Integer) = indices(x.darray, i)

DBArray_local_getindex(x, j) = getindex(localpart(x), j)

function Base.getindex(x::DBArray, i::Int)
    ipid = findfirst(rng->i∈rng, x.indices)::Int
    pid = procs(x)[ipid]
    iₒ = x.indices[ipid][1]
    remotecall_fetch(DBArray_local_getindex, pid, x, i-iₒ+1)
end

function Base.similar(x::DBArray, ::Type{T}) where {T}
    darray = DArray(I->[similar(localpart(x), T)], procs(x)[:], indices(x, 1))
    DBArray(darray, x.indices, x.blkindices)
end
Base.similar(x::DBArray{T}) where {T} = similar(x, T)
Base.similar(x::DBArray, ::Type{T}, n::Integer) where {T} = length(x) == n ? similar(x, T) : Array{T}(undef, n)
Base.similar(x::DBArray, ::Type{T}, dims::Tuple{Int}) where {T} = similar(x, T, dims[1])

DBArray_local_norm(x::DBArray, p) = norm(localpart(x), p)

function LinearAlgebra.norm(x::DBArray{T}, p::Real=2) where {T}
    _T = float(real(T))
    pids = procs(x)
    z = zeros(_T, length(pids))
    @sync for (ipid,pid) in enumerate(pids)
        @async begin
            z[ipid] = remotecall_fetch(DBArray_local_norm, pid, x, p)::_T
        end
    end
    if p == Inf
        maximum(z)
    elseif p == -Inf
        minimum(z)
    elseif p == 0 || p == 1
        sum(z)
    else
        _p = _T(p)
        mapreduce(_z->_z^_p, +, z)^(one(_T)/_p)
    end
end

DBArray_local_dot(x::DBArray, y::DBArray) = dot(localpart(x), localpart(y))

function LinearAlgebra.dot(x::DBArray{T}, y::DBArray{T}) where {T}
    pids = procs(x)
    z = zeros(T, length(pids))
    @sync for (ipid,pid) in enumerate(pids)
        @async begin
            z[ipid] = remotecall_fetch(DBArray_local_dot, pid, x, y)::T
        end
    end
    sum(z)
end

DBArray_local_extrema(x::DBArray) = extrema(localpart(x))

function Base.extrema(x::DBArray{T}) where {T}
    pids = procs(x)
    mn = zeros(T, length(pids))
    mx = zeros(T, length(pids))
    @sync for (ipid,pid) in enumerate(pids)
        @async begin
            mn[ipid],mx[ipid] = remotecall_fetch(DBArray_local_extrema, pid, x)::Tuple{T,T}
        end
    end
    minimum(mn),maximum(mx)
end

DBArray_local_mapreduce(f, op, x; kw...) = mapreduce(f, op, localpart(x); kw...)

function Base.mapreduce(f, op, x::DBArray{T}; kw...) where {T}
    pids = procs(x)
    mx = zeros(T, length(pids))
    @sync for (ipid,pid) in enumerate(pids)
        @async begin
            mx[ipid] = remotecall_fetch(DBArray_local_mapreduce, pid, f, op, x; kw...)::T
        end
    end
    mapreduce(f, op, mx)
end

DBArray_local_mean(x::DBArray) = sum(localpart(x))

function Statistics.mean(x::DBArray{T}) where {T}
    pids = procs(x)
    y = zeros(T, length(pids))
    @sync for (ipid,pid) in enumerate(pids)
        @async begin
            y[ipid] = remotecall_fetch(DBArray_local_mean, pid, x)
        end
    end
    sum(y) / length(x)
end

DBArray_local_var(x::DBArray, μ) = sum((localpart(x) .- μ).^2)

function Statistics.var(x::DBArray{T}; corrected=true, mean=nothing) where {T}
    pids = procs(x)
    μ = mean === nothing ? Statistics.mean(x) : mean
    y = zeros(T, length(pids))
    @sync for (ipid,pid) in enumerate(pids)
        @async begin
            y[ipid] = remotecall_fetch(DBArray_local_var, pid, x, μ)
        end
    end
    σ² = sum(y)
    if corrected
        σ² /= (length(x) - 1)
    else
        σ² /= length(x)
    end
    σ²
end

_localpart(x::DBArray) = localpart(x.darray)
DistributedArrays.localpart(x::DBArray) = _localpart(x)[1]::Jets.BlockArray
"""
    procs(R)
returns Julia processes associated with R
"""
Distributed.procs(x::DBArray) = procs(x.darray)
"""
    nprocs(R)
returns the number of Julia processes associated with R
"""
Distributed.nprocs(x::DBArray) = length(procs(x))
Jets.nblocks(x::DBArray) = x.blkindices[end][end]
"""
    localindices(R)
returns the indicies of R that are local to the calling process
"""
DistributedArrays.localindices(x::DBArray) = x.indices[findfirst(ipid->ipid==myid(), procs(x))]
localblockindices(x::DBArray) = x.blkindices[findfirst(ipid->ipid==myid(), procs(x))]
blockmap(x::DBArray) = x.blkindices

"""
    collect(d)
collect the block array `d` to the calling process, returning a block array
"""
function Base.collect(x::DBArray{T,A}) where {T,B,A<:Jets.BlockArray{T,B}}
    _x = B[]
    _indices = UnitRange{Int}[]
    n = 0
    for pid in procs(x)
        y = remotecall_fetch(localpart, pid, x)::Jets.BlockArray{T,B}
        _x = [_x; y.arrays]
        for i = 1:length(y.indices)
            push!(_indices, ((n+y.indices[i][1]):(n+y.indices[i][end])))
        end
        n = _indices[end][end]
    end
    Jets.BlockArray(_x, _indices)
end

"""
    convert(Array,d)
collect the block array `d` and convert it to a Julia array
"""
Base.convert(::Jets.BlockArray, x::DBArray) = collect(x)
Base.convert(::Array, x::DBArray) = convert(Array, collect(x))

DBArray_local_fill!(x::DBArray, a) = fill!(localpart(x), a)

function Base.fill!(x::DBArray, a)
    @sync for pid in procs(x)
        @async remotecall_fetch(DBArray_local_fill!, pid, x, a)
    end
    x
end

# DBArrays are vectors
function Base.reshape(x::DBArray, R::JetDSpace)
    length(x) == length(R) || error("dimension mismatch, unable to reshap distributed block arrays")
    x
end
# -->

# DBArray broadcasting implementation --<
struct DBArrayStyle <: Broadcast.AbstractArrayStyle{1} end
Base.BroadcastStyle(::Type{<:DBArray}) = DBArrayStyle()
DBArrayStyle(::Val{1}) = DBArrayStyle()

function Base.similar(bc::Broadcast.Broadcasted{DBArrayStyle}, ::Type{T}) where {S,T}
    A = find_dbarray(bc)
    similar(A)
end
find_dbarray(bc::Broadcast.Broadcasted) = find_dbarray(bc.args)
find_dbarray(args::Tuple) = find_dbarray(find_dbarray(args[1]), Base.tail(args))
find_dbarray(x) = x
find_dbarray(a::DBArray, rest) = a
find_dbarray(::Any, rest) = find_dbarray(rest)

DistributedArrays.localpart(bc::Broadcast.Broadcasted{DBArrayStyle}) = Broadcast.Broadcasted{Jets.BlockArrayStyle}(bc.f, map(arg->localpart(arg), bc.args))

function DBArray_local_copyto!(dest, bc)
    copyto!(localpart(dest), localpart(bc))
    nothing
end

function Base.copyto!(dest::DBArray, bc::Broadcast.Broadcasted{DBArrayStyle})
    @sync for pid in procs(dest)
        @async remotecall_fetch(DBArray_local_copyto!, pid, dest, bc)
    end
    dest
end
# -->

DistributedArrays.empty_localpart(T,N,::Type{A}) where {B<:Jets.BlockArray,A<:Vector{B}} = Jets.BlockArray([Array{T}(undef, ntuple(zero, N))], [0:0])
for f in (:Array, :ones, :rand, :zeros)
    @eval (Base.$f)(R::JetDSpace) = DBArray(DArray(i->[($f)(localpart(R.blkspaces)[1])], procs(R), indices(R)), R.indices, R.blkindices)
end

getblocklocal(x::DBArray, δblock) = getblock(localpart(x), δblock)
"""
    getblock(d,i)
retrive block i from a distributed block operator
"""
function Jets.getblock(x::DBArray{T,A}, iblock::Integer) where {T,B,A<:Jets.BlockArray{T,B}}
    ipid = findfirst(rng->iblock∈rng, x.blkindices)
    remotecall_fetch(getblocklocal, procs(x)[ipid], x, iblock - x.blkindices[ipid][1] + 1)::B
end
"""
    getblock!(d, i, x)
get block `i::Int` and put it into `x`
"""
function Jets.getblock!(x::DBArray{T,A}, iblock::Integer, xblock::AbstractArray) where {T,B,A<:Jets.BlockArray{T,B}}
    ipid = findfirst(rng->iblock∈rng, x.blkindices)::Int
    xblock .= remotecall_fetch(getblocklocal, procs(x)[ipid], x, iblock - x.blkindices[ipid][1] + 1)::B
end

function setblocklocal!(x::DBArray, δblock::Integer, xblock)
    setblock!(localpart(x), δblock, xblock)
    nothing
end
"""
    setblock!(d,i,x)
put `x` into block `i::Int`
"""
function Jets.setblock!(x::DBArray, iblock::Integer, xblock)
    ipid = findfirst(rng->iblock∈rng, x.blkindices)::Int
    remotecall_fetch(setblocklocal!, procs(x)[ipid], x, iblock - x.blkindices[ipid][1] + 1, xblock)
end

Jets.JopBlock(A::DArray{T,2}; kwargs...) where {T<:Jop} = JopNl(Jets.JetBlock(A; kwargs...))
Jets.JopBlock(A::DArray{T,2}; kwargs...) where {T<:JopLn} = JopLn(Jets.JetBlock(A; kwargs...))

function JetDBlock_build(I, ops::DArray{T,2}, dadom) where {T<:Jop}
    irng = indices(ops, 1)[I[1][1]]
    jrng = indices(ops, 2)[I[2][1]]
    _ops = T[ops[i,j] for i in irng, j in jrng]
    [Jets.JopBlock(_ops; dadom=dadom) for k=1:1, l=1:1]
end

JetDBlock_countidxs(spaces) = length(localpart(spaces)[1])
JetDBlock_countblks(spaces) = nblocks(localpart(spaces)[1])

function JetDBlock_count_idxs_blks(blkspaces)
    pids = procs(blkspaces)
    n = length(blkspaces)
    countidxs = zeros(Int, n)
    countblks = zeros(Int, n)
    @sync for (i,pid) in enumerate(pids)
        @async begin
            countidxs[i] = remotecall_fetch(JetDBlock_countidxs, pid, blkspaces)::Int
            countblks[i] = remotecall_fetch(JetDBlock_countblks, pid, blkspaces)::Int
        end
    end
    countidxs,countblks
end

function JetDBlock_collect_idxs(countidxs, countblks, pids)
    idxs = Vector{UnitRange{Int}}(undef, length(pids))
    lst = 0
    for i = 1:length(idxs)
        frst = lst + 1
        lst = frst + countidxs[i] - 1
        idxs[i] = frst:lst
    end
    blkidxs = Vector{UnitRange{Int}}(undef, length(pids))
    lst = 0
    for i = 1:length(idxs)
        frst = lst + 1
        lst = frst + countblks[i] - 1
        blkidxs[i] = frst:lst
    end
    idxs,blkidxs
end

function JetDBlock_space(blkspaces, blkidxs, idxs)
    length(blkidxs) == 1 && (return blkspaces[1])
    JetDSpace(blkspaces, blkidxs, idxs)
end

function localpart_diag_darray_spaces(_ops)
    ops = state(localpart(_ops)[1]).ops
    for j = 1:size(ops,1)
        opsj = @view ops[j,:]
        s = mapreduce(!iszero, +, opsj)
        s == 1 || error("Expecting one non-zero column block per row in block diagonal construction, got $s.")
    end
    ops = state(localpart(_ops)[1]).ops[:,1:1]
    A = @blockop adjoint.(JopLn.(jet.(ops))) dadom=true
    range(A)
end

function Jets.JetBlock(ops::DArray{T,2}; perfstatfile="", isdiag=false) where {T<:Jop}
    pids = procs(ops)

    n1,n2 = length(indices(ops, 1)),length(indices(ops, 2))

    (isdiag && n2 != 1) && error("block diagonal jets must define distribution along rows.")

    _ops = DArray(I->JetDBlock_build(I, ops, n2 > 1), (n1,n2), pids[:], [n1,n2])

    blkspaces_rng = DArray(I->[range(localpart(_ops)[1]) for i in I[1]], (n1,), pids[:,1], [n1])

    local blkspaces_dom
    if isdiag
        # TODO: fragile
        # assumption that there is only one column block per block-diagnol block
        blkspaces_dom = DArray(I->[localpart_diag_darray_spaces(_ops) for i in I[1]], (n1,), pids[:,1], [n1])
    else
        blkspaces_dom = DArray(I->[domain(localpart(_ops)[1]) for i in I[1]], (n2,), pids[1,:], [n2])
    end

    countidxs_dom,countblks_dom = JetDBlock_count_idxs_blks(blkspaces_dom)
    countidxs_rng,countblks_rng = JetDBlock_count_idxs_blks(blkspaces_rng)

    idxs_rng,blkidxs_rng = JetDBlock_collect_idxs(countidxs_rng, countblks_rng, pids[:,1])
    idxs_dom,blkidxs_dom = JetDBlock_collect_idxs(countidxs_dom, countblks_dom, isdiag ? pids[:,1] : pids[1,:])

    dom = JetDBlock_space(blkspaces_dom, blkidxs_dom, idxs_dom)
    rng = JetDBlock_space(blkspaces_rng, blkidxs_rng, idxs_rng)

    perfstatfile == "" || rm(perfstatfile, force=true)

    Jet(dom = dom, rng = rng, f! = JetDBlock_f!, df! = JetDBlock_df!, df′! = JetDBlock_df′!, s = (ops=_ops, dom=dom, rng=rng, blockmap=indices(ops), isdiag=isdiag, perfstatfile=perfstatfile))
end

function addmasterpid(pids)
    if myid() ∉ pids
        return [myid();pids]
    end
    pids
end

function JetDBlock_local_f!(d, _m, ops)
    op = localpart(ops)[1]
    mul!(localpart(d), op, localpart(_m))
    nothing
end

function _JetDBlock_f!(d::DBArray, m::AbstractArray, ops)
    pids = procs(ops)
    _m = bcast(m, addmasterpid(pids))
    @sync for pid in pids
        @async remotecall_fetch(JetDBlock_local_f!, pid, d, _m, ops)
    end
    nothing
end

function JetDBlock_f!(d::DBArray, m::AbstractArray; ops, perfstatfile, kwargs...)
    _JetDBlock_f!(d, m, ops)
    perfstat(ops, "f", perfstatfile)
    d
end

function JetDBlock_local_df!(d, _m, ops)
    op = localpart(ops)[1]
    mul!(localpart(d), JopLn(op), localpart(_m))
    nothing
end

function _JetDBlock_df!(d::DBArray, m::AbstractArray, ops)
    pids = procs(ops)
    _m = bcast(m, addmasterpid(pids))
    @sync for pid in pids
        @async remotecall_fetch(JetDBlock_local_df!, pid, d, _m, ops)
    end
    nothing
end

function JetDBlock_df!(d::DBArray, m::AbstractArray; ops, perfstatfile, kwargs...)
    _JetDBlock_df!(d, m, ops)
    perfstat(ops, "df", perfstatfile)
    d
end

function JetDBlock_local_df′!(_m, d, ops)
    op = localpart(ops)[1]
    mul!(localpart(_m), (JopLn(op))', localpart(d))
    nothing
end

function _JetDBlock_df′!(m::AbstractArray, d::DBArray, ops, dom)
    m .= 0
    pids = procs(ops)
    _m = TypeFutures(m, zeros, addmasterpid(pids), dom)
    @sync for pid in pids
        @async remotecall_fetch(JetDBlock_local_df′!, pid, _m, d, ops)
    end
    reduce!(_m)
    nothing
end

function JetDBlock_df′!(m::AbstractArray, d::DBArray; ops, dom, perfstatfile, kwargs...)
    _JetDBlock_df′!(m, d, ops, dom)
    perfstat(ops, "df′", perfstatfile)
    m
end

function JetDBlock_local_row_cols_f!(_d, m, ops, blockmap, irow_pid, icol_pid, irow)
    __d = localpart(_d)
    dtmp = similar(__d)
    for icol in blockmap[irow_pid,icol_pid][2]
        op = getblock(ops, blockmap, irow, icol)
        if !iszero(op)
            __d .+= mul!(dtmp, op, getblock(m, icol))
        end
    end
    nothing
end

function JetDBlock_local_rows_f!(d, m, ops, irow_pid, blockmap)
    pids = procs(ops)[irow_pid,:]
    _pids = addmasterpid(pids)
    for irow in blockmap[irow_pid,1][1]
        dblk = getblock(d, irow)
        dblk .= 0
        _d = ArrayFutures(dblk, _pids)
        @sync for (icol_pid,col_pid) in enumerate(pids)
            @async remotecall_fetch(JetDBlock_local_row_cols_f!, col_pid, _d, m, ops, blockmap, irow_pid, icol_pid, irow)
        end
        reduce!(_d)
    end
    nothing
end

function JetDBlock_local_rows_diag_f!(d, m, ops, irow_pid, blockmap)
    for irow in blockmap[irow_pid,1][1]
        mul!(getblock(d, irow), getblock(ops, blockmap, irow, irow), getblock(m, irow))
    end
    nothing
end

function JetDBlock_f!(d::DBArray, m::DBArray; ops, blockmap, isdiag, perfstatfile, kwargs...)
    pids = procs(ops)
    @sync for (irow_pid,row_pid) in enumerate(pids[:,1])
        if isdiag
            @async remotecall_fetch(JetDBlock_local_rows_diag_f!, row_pid, d, m, ops, irow_pid, blockmap)
        else
            @async remotecall_fetch(JetDBlock_local_rows_f!, row_pid, d, m, ops, irow_pid, blockmap)
        end
    end
    perfstat(ops, "f", perfstatfile)
    d
end

function JetDBlock_local_row_cols_df!(_d, m, ops, blockmap, irow_pid, icol_pid, irow)
    __d = localpart(_d)
    dtmp = similar(__d)
    for icol in blockmap[irow_pid,icol_pid][2]
        op = getblock(ops, blockmap, irow, icol)
        if !iszero(op)
            __d .+= mul!(dtmp, JopLn(op), getblock(m, icol))
        end
    end
    nothing
end

function JetDBlock_local_rows_df!(d, m, ops, irow_pid, blockmap)
    pids = procs(ops)[irow_pid,:]
    _pids = addmasterpid(pids)
    for irow in blockmap[irow_pid,1][1]
        dblk = getblock(d, irow)
        dblk .= 0
        _d = ArrayFutures(dblk, _pids)
        @sync for (icol_pid,col_pid) in enumerate(pids)
            @async remotecall_fetch(JetDBlock_local_row_cols_df!, col_pid, _d, m, ops, blockmap, irow_pid, icol_pid, irow)
        end
        reduce!(_d)
    end
    nothing
end

function JetDBlock_local_rows_diag_df!(d, m, ops, irow_pid, blockmap)
    for irow in blockmap[irow_pid,1][1]
        mul!(getblock(d, irow), JopLn(getblock(ops, blockmap, irow, irow)), getblock(m, irow))
    end
    nothing
end

function JetDBlock_df!(d::DBArray, m::DBArray; ops, blockmap, isdiag, perfstatfile, kwargs...)
    pids = procs(ops)
    @sync for (irow_pid,row_pid) in enumerate(pids[:,1])
        if isdiag
            @async remotecall_fetch(JetDBlock_local_rows_diag_df!, row_pid, d, m, ops, irow_pid, blockmap)
        else
            @async remotecall_fetch(JetDBlock_local_rows_df!, row_pid, d, m, ops, irow_pid, blockmap)
        end
    end
    perfstat(ops, "df", perfstatfile)
    d
end

function JetDBlock_local_col_rows_df′!(_m, d, ops, blockmap, icol_pid, irow_pid, icol)
    __m = localpart(_m)
    mtmp = similar(__m)
    for irow in blockmap[irow_pid,icol_pid][1]
        op = getblock(ops, blockmap, irow, icol)
        if !iszero(op)
            __m .+= mul!(mtmp, (JopLn(op))', getblock(d, irow))
        end
    end
    nothing
end

function JetDBlock_local_cols_df′!(m, d, ops, icol_pid, blockmap)
    pids = procs(ops)[:,icol_pid]
    _pids = addmasterpid(pids)
    for icol in blockmap[1,icol_pid][2]
        mblk = getblock(m, icol)
        mblk .= 0
        _m = ArrayFutures(mblk, _pids)
        @sync for (irow_pid,row_pid) in enumerate(pids)
            @async remotecall_fetch(JetDBlock_local_col_rows_df′!, row_pid, _m, d, ops, blockmap, icol_pid, irow_pid, icol)
        end
        reduce!(_m)
    end
    nothing
end

function JetDBlock_local_cols_diag_df′!(m, d, ops, icol_pid, blockmap)
    for icol in blockmap[icol_pid,1][1]
        mul!(getblock(m, icol), (JopLn(getblock(ops, blockmap, icol, icol)))', getblock(d, icol))
    end
    nothing
end

function JetDBlock_df′!(m::DBArray, d::DBArray; ops, blockmap, isdiag, perfstatfile, kwargs...)
    pids = isdiag ? (procs(ops))' : procs(ops)
    @sync for (icol_pid,col_pid) in enumerate(pids[1,:])
        if isdiag
            @async remotecall_fetch(JetDBlock_local_cols_diag_df′!, col_pid, m, d, ops, icol_pid, blockmap)
        else
            @async remotecall_fetch(JetDBlock_local_cols_df′!, col_pid, m, d, ops, icol_pid, blockmap)
        end
    end
    perfstat(ops, "df′", perfstatfile)
    m
end

function _perfstat_local(ops, operation)
    _ops = state(localpart(ops)[1]).ops
    _s = [perfstat(_ops[iblkrow,iblkcol]) for iblkrow=1:size(_ops,1), iblkcol=1:size(_ops,2)][:]
    Dict("id"=>myid(), "hostname"=>gethostname(), "operation"=>operation, "localblock"=>_s)
end

function Jets.perfstat(ops, operation, perfstatfile)
    if perfstatfile != ""
        pids = procs(ops)
        s = Vector{Dict{String,Any}}(undef, length(pids))
        @sync for (i, pid) in enumerate(procs(ops))
            @async s[i] = remotecall_fetch(_perfstat_local, pid, ops, operation)
        end
        if isfile(perfstatfile)
            _s = JSON.parse(read(perfstatfile, String))
            push!(_s["step"], Dict("pid"=>s))
            write(perfstatfile, json(_s, 1))
        else
            write(perfstatfile, json(Dict("step"=>[Dict("pid"=>s)]), 1))
        end
    end
    nothing
end

function JetDBlock_local_point!(ops, _mₒ)
    op = localpart(ops)[1]
    Jets.point!(jet(op), localpart(_mₒ))
    nothing
end

function Jets.point!(j::Jet{D,R,typeof(JetDBlock_f!)}, mₒ::AbstractArray) where {D<:Jets.JetAbstractSpace,R<:Jets.JetAbstractSpace}
    ops = state(j).ops
    pids = procs(ops)
    _mₒ = bcast(mₒ, addmasterpid(pids))
    @sync for pid in pids
        @async remotecall_fetch(JetDBlock_local_point!, pid, ops, _mₒ)
    end
    j
end

function JetDBlock_local_row_cols_point!(mₒ, ops, blockmap, irow_pid, icol_pid, irow)
    icols = blockmap[irow_pid,icol_pid][2]
    for icol in icols
        op = getblock(ops, blockmap, irow, icol)
        if !iszero(op)
            Jets.point!(jet(op), getblock(mₒ, icol))
        end
    end
    nothing
end

function JetDBlock_local_rows_point!(mₒ, ops, irow_pid, blockmap)
    pids = procs(ops)[irow_pid,:]
    _pids = addmasterpid(pids)
    for irow in blockmap[irow_pid,1][1]
        @sync for (icol_pid,col_pid) in enumerate(pids)
            @async remotecall_fetch(JetDBlock_local_row_cols_point!, col_pid, mₒ, ops, blockmap, irow_pid, icol_pid, irow)
        end
    end
    nothing
end

function JetDBlock_local_rows_diag_point!(mₒ, ops, irow_pid, blockmap)
    for irow in blockmap[irow_pid,1][1]
        op = getblock(ops, blockmap, irow, irow)
        Jets.point!(jet(op), getblock(mₒ, irow))
    end
    nothing
end

function Jets.point!(j::Jet{D,R,typeof(JetDBlock_f!)}, mₒ::DBArray) where {D<:Jets.JetAbstractSpace,R<:Jets.JetAbstractSpace}
    pids = procs(j)
    @sync for (irow_pid,row_pid) in enumerate(pids[:,1])
        if state(j).isdiag
            @async remotecall_fetch(JetDBlock_local_rows_diag_point!, row_pid, mₒ, state(j).ops, irow_pid, state(j).blockmap)
        else
            @async remotecall_fetch(JetDBlock_local_rows_point!, row_pid, mₒ, state(j).ops, irow_pid, state(j).blockmap)
        end
    end
    j
end

Distributed.procs(jet::J) where {D,R,J<:Jet{D,R,typeof(JetDBlock_f!)}} = procs(state(jet).ops)
Distributed.procs(A::T) where {D,R,J<:Jet{D,R,typeof(JetDBlock_f!)},T<:Jop{J}} = procs(jet(A))
Distributed.nprocs(A::T) where {D,R,J<:Jet{D,R,typeof(JetDBlock_f!)},T<:Jop{J}} = length(procs(A))

DistributedArrays.localpart(j::Jet{D,R,typeof(JetDBlock_f!)}) where {D,R} = jet(localpart(state(j).ops)[1])
DistributedArrays.localpart(A::T) where {D,R,J<:Jet{D,R,typeof(JetDBlock_f!)},T<:JopLn{J}} = JopLn(localpart(state(A).ops)[1])
DistributedArrays.localpart(A::T) where {D,R,J<:Jet{D,R,typeof(JetDBlock_f!)},T<:JopAdjoint{J}} = JopAdjoint(localpart(state(A).ops)[1])
DistributedArrays.localpart(A::T) where {D,R,J<:Jet{D,R,typeof(JetDBlock_f!)},T<:Jop{J}} = localpart(state(A).ops)[1]

function localblockindices(j::Jet{D,R,typeof(JetDBlock_f!)}) where {D,R}
    blkmap = blockmap(j)
    pids = procs(j)
    mypid = myid()
    for i2 = 1:size(pids,2), i1 = 1:size(pids,1)
        if mypid == pids[i1,i2]
            return blkmap[i1,i2]
        end
    end
    (0:0,0:0)
end
localblockindices(j::Jet{D,R,typeof(JetDBlock_f!)}, i) where {D,R} = localblockindices(j)[i]
localblockindices(A::T) where {D,R,J<:Jet{D,R,typeof(JetDBlock_f!)},T<:Jop{J}} = localblockindices(jet(A))
localblockindices(A::T, i) where {D,R,J<:Jet{D,R,typeof(JetDBlock_f!)},T<:Jop{J}} = localblockindices(jet(A), i)

localblockindices(j::Jet{D,R,typeof(Jets.JetComposite_f!)}) where {D,R} = error("localblockindices does not support composite operators; consider using `localblockindices(range(F))` or `localblockindices(domain(F))`")
localblockindices(A::T) where {D,R,J<:Jet{D,R,typeof(Jets.JetComposite_f!)},T<:Jop{J}} = localblockindices(jet(A))
localblockindices(A::T, i) where {D,R,J<:Jet{D,R,typeof(Jets.JetComposite_f!)},T<:Jop{J}} = localblockindices(jet(A))

blockmap(jet::T) where {D,R,T<:Jet{D,R,typeof(JetDBlock_f!)}} = state(jet).blockmap
blockmap(op::Jop) = blockmap(jet(op))

JetDBlock_local_getblock(ops, δi, δj) = getblock(localpart(ops)[1], δi, δj)

function Jets.getblock(ops::DArray, blkmap, i::Integer, j::Integer)
    local ipid, jpid
    for _ipid in 1:size(blkmap,1), _jpid in 1:size(blkmap,2)
        if i ∈ blkmap[_ipid,_jpid][1] && j ∈ blkmap[_ipid,_jpid][2]
            ipid,jpid = _ipid,_jpid
            break
        end
    end
    pid = procs(ops)[ipid,jpid]
    remotecall_fetch(JetDBlock_local_getblock, pid, ops, i - blkmap[ipid,jpid][1][1] + 1, j - blkmap[ipid,jpid][2][1] + 1)
end

Jets.isblockop(A::Jop{<:Jet{<:JetAbstractSpace,<:JetAbstractSpace,typeof(JetDBlock_f!)}}) = true

Jets.getblock(jet::Jet{D,R,typeof(JetDBlock_f!)}, i, j) where {D,R} = getblock(state(jet).ops, blockmap(jet), i, j)

close_localpart(ops::DArray) = close(localpart(ops)[1])

function Base.close(j::Jet{D,R,typeof(JetDBlock_f!)}) where {D,R}
    ops = state(j).ops
    @sync for pid in procs(ops)
        @async remotecall_fetch(close_localpart, pid, ops)
    end
end

export DBArray, JetDSpace, blockmap, localblockindices

end