text:
export codegen, JuliaASTCodegenCtx

abstract type CompileCtx end
abstract type AbstractJuliaASTCtx <: CompileCtx end

const codegen_passes = Dict{Symbol,Function}()

# abstract type AbstractQASMCtx <: CompileCtx end
# struct QASMCodegenCtx <: AbstractQASMCtx
#     codegen_pass::Vector{Any}
#     registers::Vector{Any}
# end

# function codegen(ctx::QASMCodegenCtx, ir::YaoIR)
#     lines = String[]
#     for pass in ctx.codegen_pass
#         push!(lines, codegen_passes[pass](ctx, ir))
#     end
#     return join(["// code generated by Yao",
#         "OPENQASM 2.0;", join(lines, ";\n"), "\n"], "\n")
# end

struct JuliaASTCodegenCtx <: AbstractJuliaASTCtx
    stub_name::Any
    circuit::Symbol
    registers::Vector{Any}
    locations::Symbol
    ctrl_locations::Symbol
    codegen_pass::Vector{Any}
end

# generate to JuliaAST by default
codegen(ir::YaoIR) = codegen(JuliaASTCodegenCtx(ir), ir)

function codegen(ctx::JuliaASTCodegenCtx, ir)
    ex = Expr(:block)
    for pass in ctx.codegen_pass
        push!(ex.args, codegen_passes[pass](ctx, ir))
    end

    if ir.name isa Symbol
        push!(ex.args, :($(ir.name)))
    else
        push!(ex.args, nothing)
    end
    return ex
end

macro codegen(ex)
    defs = splitdef(ex)
    defs[:name] isa Symbol ||
        throw(ParseError("@codegen expect function not callable or lambda function"))
    name = Symbol(:codegen_, defs[:name])
    quoted_name = QuoteNode(defs[:name])
    defs[:name] = name
    defs[:body] = quote
        @timeit_debug to $(string(name)) begin
            $(defs[:body])
        end
    end

    quote
        codegen_passes[$(quoted_name)] = $(esc(combinedef(defs)))
    end
end

flatten_locations(parent, x) = IRTools.xcall(Base, :getindex, parent, x)

function flatten_locations(pr, v, parent, x)
    loc = insert!(pr, v, Statement(to_locations(x)))
    return insert!(pr, v, Statement(flatten_locations(parent, loc)))
end

# merge location in runtime
merge_location_ex(l1, l2) = :(merge_locations($l1, $l2))
# merge literal location in compile time
merge_location_ex(l1::AbstractLocations, l2::AbstractLocations) = merge_locations(l1, l2)

function extract_closure_captured_variables!(pr::IRTools.Pipe, circ, ir::YaoIR)
    v = push!(pr, IRTools.xcall(Base, :getproperty, circ, QuoteNode(:free)))
    for (k, each) in enumerate(arguements(ir))
        x = push!(pr, IRTools.xcall(Base, :getfield, v, k))
        push!(pr, Statement(Expr(:(=), IRTools.Slot(each), x)))
    end
end

function scan_registers(register, stack::Vector, pr::IRTools.Pipe, v, st::Statement)
    if st.expr.args[2] === :new
        push!(stack, st.expr.args[3])
    elseif st.expr.args[2] === :prev
        pop!(stack)
    end

    if length(stack) == 1
        delete!(pr, v)
        return register
    else
        pr[v] = Statement(IRTools.Slot(first(stack)))
        return v
    end
end

@codegen function call(ctx::JuliaASTCodegenCtx, ir::YaoIR)
    @timeit_debug to "create defs" defs = signature(ir)
    if ir.name isa Symbol
        @timeit_debug to "create name" defs[:name] = :(::$(generic_circuit(ir.name)))
    end
    @timeit_debug to "create body" defs[:body] =
        Expr(:call, circuit(ir.name), ctx.stub_name, Expr(:tuple, arguements(ir)...))
    @timeit_debug to "combinedef" code = combinedef(defs)
    return code
end

@codegen function evaluate(ctx::JuliaASTCodegenCtx, ir::YaoIR)
    @timeit_debug to "create defs" defs = signature(ir)
    @timeit_debug to "create name" defs[:name] = GlobalRef(YaoLang, :evaluate)

    if ir.name isa Symbol
        @timeit_debug to "create args" defs[:args] = Any[:(::$(generic_circuit(ir.name))), ir.args...]
    else
        @timeit_debug to "create args" defs[:args] = Any[ir.name, ir.args...]
    end
    @timeit_debug to "create body" defs[:body] =
        Expr(:call, circuit(ir.name), ctx.stub_name, Expr(:tuple, arguements(ir)...))
    return combinedef(defs)
end

@codegen function circuit(ctx::JuliaASTCodegenCtx, ir::YaoIR)
    empty!(ctx.registers)
    pr = IRTools.Pipe(ir.body)
    circ = IRTools.argument!(pr)
    register = IRTools.argument!(pr)
    locations = IRTools.argument!(pr)

    # extract arguements from closure
    @timeit_debug to "extract variables" extract_closure_captured_variables!(pr, circ, ir)

    @timeit_debug to "IRTools.Pipe loop" for (v, st) in pr
        if is_quantum(st)
            head = st.expr.args[1]
            if head === :register
                register = scan_registers(register, ctx.registers, pr, v, st)
            elseif head in [:gate, :ctrl]
                locs = map(x -> flatten_locations(pr, v, locations, x), st.expr.args[3:end])
                pr[v] = Statement(st; expr = Expr(:call, st.expr.args[2], register, locs...))
            elseif head === :measure
                measure_ex = Expr(:call, GlobalRef(YaoAPI, :measure!))

                if (st.expr.args[2] isa Expr) && (st.expr.args[2].head === :parameters)
                    kwargs = first(st.expr.args[2])
                    if kwargs.args[1] === :reset_to
                        push!(measure_ex.args, Expr(:call, ResetTo, kwargs.args[2]))
                    elseif (kwargs.args[1] === :remove) && (kwargs.args[2] == true)
                        push!(measure_ex.args, Expr(:call, RemoveMeasured))
                    end
                    args = st.expr.args[3:end]
                else
                    args = st.expr.args[2:end]
                end

                # contains operator
                if length(args) == 2
                    push!(measure_ex.args, args[1])
                    push!(measure_ex.args, register)
                end

                push!(measure_ex.args, register)
                loc = flatten_locations(pr, v, locations, args[end])
                push!(measure_ex.args, loc)

                pr[v] = Statement(st; expr = measure_ex)
            else # reserved for extending keywords
                throw(ParseError("Invalid keyword: $head"))
            end
        end
    end

    def = Dict{Symbol,Any}()
    def[:name] = ctx.stub_name
    def[:args] = Any[
        :($(ctx.circuit)::$(YaoLang.Circuit)),
        :($(last(ctx.registers))::$AbstractRegister),
        :($(ctx.locations)::Locations),
    ]

    @timeit_debug to "IRTools.finish" ssa = IRTools.finish(pr)
    @timeit_debug to "build_codeinfo" code = build_codeinfo(ir.mod, def, ssa)
    return code
end

@codegen function ctrl_circuit(ctx::JuliaASTCodegenCtx, ir::YaoIR)
    if hasmeasure(ir)
        def = Dict{Symbol,Any}()
        def[:name] = ctx.stub_name
        def[:args] = Any[
            :($(ctx.circuit)::$(YaoLang.Circuit)),
            :($(gensym(:register))::$AbstractRegister),
            :($(ctx.locations)::$Locations),
            :($(ctx.ctrl_locations)::$CtrlLocations),
        ]
        def[:body] = quote
            error("cannot control quantum program contains measurement")
        end
        return combinedef(def)
    end

    empty!(ctx.registers)

    pr = IRTools.Pipe(ir.body)
    circ = IRTools.argument!(pr)
    register = IRTools.argument!(pr)
    locations = IRTools.argument!(pr)
    ctrl_locations = IRTools.argument!(pr)

    # extract arguements from closure
    @timeit_debug to "extract variables" extract_closure_captured_variables!(pr, circ, ir)

    @timeit_debug to "IRTools.Pipe loop" for (v, st) in pr
        if is_quantum(st)
            head = st.expr.args[1]
            if head === :register
                register = scan_registers(register, ctx.registers, pr, v, st)
            elseif head === :gate
                locs = flatten_locations(pr, v, locations, st.expr.args[3])
                pr[v] =
                    Statement(st; expr = Expr(:call, st.expr.args[2], register, locs, ctrl_locations))
            elseif head === :ctrl
                locs = flatten_locations(pr, v, locations, st.expr.args[3])
                ctrl_locs = flatten_locations(pr, v, locations, st.expr.args[4])
                ctrl_locs = insert!(pr, v, merge_location_ex(ctrl_locations, ctrl_locs))
                pr[v] = Statement(st; expr = Expr(:call, st.expr.args[2], register, locs, ctrl_locs))
            else # reserved for extending keywords
                throw(ParseError("Invalid keyword: $head"))
            end
        end
    end

    def = Dict{Symbol,Any}()
    def[:name] = ctx.stub_name
    def[:args] = Any[
        :($(ctx.circuit)::$(YaoLang.Circuit)),
        :($(last(ctx.registers))::$AbstractRegister),
        :($(ctx.locations)::$Locations),
        :($(ctx.ctrl_locations)::$CtrlLocations),
    ]

    @timeit_debug to "IRTools.finish" ssa = IRTools.finish(pr)
    @timeit_debug to "build_codeinfo" code = build_codeinfo(ir.mod, def, ssa)
    return code
end

@codegen function create_symbol(ctx::JuliaASTCodegenCtx, ir::YaoIR)
    # only create symbol when its a function declaration
    if ir.name isa Symbol
        :(Core.@__doc__ const $(ir.name) = $(generic_circuit(ir.name))())
    end
end

@codegen function code_yao_runtime_stub(ctx::JuliaASTCodegenCtx, ir::YaoIR)
    def = Dict{Symbol,Any}()
    def[:name] = GlobalRef(YaoLang, :code_yao)
    def[:args] = Any[:(::$(generic_circuit(ir.name))), ir.args...]
    def[:body] = ir
    if !isempty(ir.whereparams)
        def[:whereparams] = ir.whereparams
    end
    return combinedef(def)
end

function JuliaASTCodegenCtx(ir::YaoIR, pass = collect(Any, keys(codegen_passes)))
    stub_name = gensym(rm_annotations(ir.name))
    JuliaASTCodegenCtx(
        stub_name,
        gensym(:circ),
        Any[gensym(:register)],
        gensym(:locations),
        gensym(:ctrl_locations),
        pass,
    )
end


text:
include("item_expected_information.jl")
include("item_observed_information.jl")
include("latent_information.jl")

##########################################
#   Latent information                   #
##########################################
"""
```julia
latent_information(examinee::AbstractExaminee, item::AbstractItem)
```

# Description

An abstraction of `_latent_information(latent::AbstractLatent, parameters::AbstractParametersBinary)` on an examinee and an item.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`examinee::AbstractExaminee`** : Required.
- **`item::AbstractItem`** : Required.

# Output
A `Float64` scalar.
"""
function latent_information(examinee::AbstractExaminee, item::AbstractItem)
    _latent_information(examinee.latent, item.parameters)::Float64
end

"""
```julia
latent_information(examinee::AbstractExaminee, items::Vector{<:AbstractItem})
```

# Description

An abstraction of `_latent_information(latent::AbstractLatent, parameters::AbstractParametersBinary)` on an examinee and items.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`examinee::AbstractExaminee`** : Required.
- **`items::Vector{<:AbstractItem}`** : Required.

# Output
A `Float64` vector.

# Example

Compute the Fisher information for the latent/s of the examinees (second derivatives of the likelihood with respect to ``\theta``) and each item.

````E_\theta_n \\[ I(\theta_n | b_i) \\] ````

```julia
examinee = Examinee(); #default examinee factory
items = [Item() for i = 1 : 30]; #default item factory
latent_information(examinee, items) #compute the information wrt θ
```
"""
function latent_information(examinee::AbstractExaminee, items::Vector{<:AbstractItem})
    mapreduce( i -> _latent_information(examinee.latent, i.parameters), vcat, items)::Vector{Float64}
end

"""
```julia
latent_information(
    examinees::Vector{<:AbstractExaminee},
    items::Vector{<:AbstractItem},
)
```

# Description

An abstraction of `_latent_information(latent::AbstractLatent, parameters::AbstractParametersBinary)` on examinees and items.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`examinees::Vector{<:AbstractExaminee}`** : Required.
- **`items::Vector{<:AbstractItem}`** : Required.

# Output
A ``N \times I`` matrix of latent informations (matrices or scalars).

```@example
examinees = [Examinee(); for n = 1 : 100]; #default examinee factory
items = [Item() for i = 1 : 30]; #default item factory
latent_information(examinees, items) #compute the information wrt θs
```
"""
function latent_information(
    examinees::Vector{<:AbstractExaminee},
    items::Vector{<:AbstractItem},
)
    return [_latent_information(e.latent, i.parameters) for i in items, e in examinees]
end

##########################################
#   Item expected information            #
##########################################
"""
    item_expected_information(
        item::AbstractItem,
        examinee::AbstractExaminee,
    )

# Description

Abstraction of _item_expected_information(latent, parameters) on Vector{<:AbstractExaminee} and items::Vector{<:AbstractItem}.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`item::AbstractItem`** : Required.
- **`examinee::AbstractExaminee`** : Required.

# Output
A matrix of expected informations (matrices for 2PL and 3PL or scalars 1PL).
"""
function item_expected_information(
    item::AbstractItem,
    examinee::AbstractExaminee,
)
    return _item_expected_information(item.parameters, examinee.latent)
end

"""
```julia
item_expected_information(
    items::Vector{<:AbstractItem},
    examinees::Vector{<:AbstractExaminee},
)
```

# Description

Abstraction of _item_expected_information(latent, parameters) on Vector{<:AbstractExaminee} and items::Vector{<:AbstractItem}.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`items::Vector{<:AbstractItem}`** : Required. Size I x 1.
- **`examinees::Vector{<:AbstractExaminee}`** : Required. Size N x 1.

# Output
A ``N \times I`` matrix of expected informations (matrices or scalars).

```@example
examinee = Examinee(1); # default examinee random factory (1-D latent)
item_1PL = Item(1, Parameters1PL()); # default 1PL item random factory
response_1PL = answer(examinee, item_1PL) # generate responses
item_exp_info_1PL = item_expected_information(item_1PL, examinee)

item_2PL = Item(1, Parameters2PL()); # default 2PL item random factory
response_2PL = answer(examinee, item_2PL) # generate responses
item_exp_info_2PL = item_expected_information(item_2PL, examinee)
"""
function item_expected_information(
    items::Vector{<:AbstractItem},
    examinees::Vector{<:AbstractExaminee},
)
    return [
        _item_expected_information(i.parameters, e.latent) for i in items, e in examinees
    ]
end

##########################################
#   Item observed information            #
##########################################

"""
```julia
item_observed_information(
    item::AbstractItem,
    examinee::AbstractExaminee,
    response::AbstractResponse,
)
```

# Description

It computes the item observed information for an `examinee` who answered `response` .

```@example
examinee = Examinee(1); # default examinee random factory (1-D latent)
item_2PL = Item(1, Parameters2PL()); # default 2PL item random factory
response_2PL = answer(examinee, item_2PL) # generate responses
item_obs_info_2PL = item_observed_information(item_2PL, examinee, response_2PL)
item_exp_info_2PL = item_expected_information(item_2PL, examinee)
# for 2PL and 1PL item parameters, the observed information is equal to the expected information.
@assert item_obs_info_2PL == item_exp_info_2PL

examinee = Examinee(1); # default examinee random factory (1-D latent)
item_3PL = Item(1, Parameters3PL()); # default 3PL item random factory
response_3PL = answer(examinee, item_3PL) # generate responses
# for 3PL item parameters, the observed information is different from the expected information.
item_obs_info_3PL = item_observed_information(item_3PL, examinee, response_3PL)
item_exp_info_3PL = item_expected_information(item_3PL, examinee)
```
"""
function item_observed_information(
    item::AbstractItem,
    examinee::AbstractExaminee,
    response::AbstractResponse,
)
    return  _item_observed_information(item.parameters, examinee.latent, response.val)
end

"""
```julia
item_observed_information(
    items::Vector{<:AbstractItem},
    examinees::Vector{<:AbstractExaminee},
    responses::Vector{<:AbstractResponse},
)
```

# Description

It computes the sum of the item observed informations across the responses.
Items (latents) must be of the same type.

# Example 1

```@example
examinees = [Examinee(n) for n = 1 : 100]; #default examinee factory (1-D latent)
items_2PL = [Item(i, Parameters2PL()) for i = 1 : 30]; #2PL item factory
items_3PL = [Item(i + 30, string(i + 30), Parameters3PL()) for i = 1 : 30]; #3PL item factory, id and idx must be different
items = vcat([items_2PL, items_3PL]...) # create vector of items

responses = answer(examinees, items) # generate responses
#obs_item_info = item_observed_information(items, examinees, responses) # do not run
#^ This returns an error since Julia cannot sum ``2 \times 2`` matrices with ``3 \times 3`` matrices.
# Instead do this ⌄.
obs_items_info_2PL = item_observed_information(items_2PL, examinees, responses)
obs_items_info_3PL = item_observed_information(items_3PL, examinees, responses)
```

# Example 2. Compute the standard errors of item parameter estimates using the observed responses and estimated abilities.

```@example
# The package LinearAlgebra is required to compute the inverse of the matrix.
using LinearAlgebra

# Randomly generate examinees, items and responses.
examinees = [Examinee(n) for n = 1 : 10]
items = [Item(i) for i = 1:40];
responses = answer(examinees, items);
responses_per_item = map(i -> get_responses_by_item_id(i.id, responses), items);

items_obs_info = [sum([item_observed_information(items[i], examinees[r.examinee_idx], r) for r in responses_per_item[i]] for i in 1:40];

#Using the inverse of the matrix:
inv_items_obs_info = inv.(items_obs_info);
item_parameters_standard_errors_inv = [[sqrt(i[1,1]), sqrt(i[2,2])] for i in inv_items_obs_info]

#Using the inverse of the diagonal:
item_parameters_standard_errors_diag = [[sqrt(1/i[1,1]), sqrt(1/i[2,2])] for i in items_obs_info]
```
"""
function item_observed_information(
    items::Vector{<:AbstractItem},
    examinees::Vector{<:AbstractExaminee},
    responses::Vector{<:AbstractResponse},
)
    return mapreduce( r ->
        item_observed_information(
            r,
            get_examinee_by_id(r.examinee_id, examinees),
            get_item_by_id(r.item_id, items),
        ),
        +,
        responses)
end

##########################################
#   Externals for matrix data            #
##########################################
"""
```julia
latent_information(
    latents_matrix::Matrix{Float64},
    parameters_matrix::Matrix{Float64},
)
```

# Description

It computes the information function (IIF) for item parameters at latents values provided in matrix form.
Not suitable for 3PL models, for such a kind of model use `latent_information_3PL()`.
It follows the parametrization ``a \theta - b ``.
See the docs of [`latent_information(examinee::AbstractExaminee, item::AbstractItem)`](@ref) for details.

# Arguments

- **`latents_matrix::Matrix{Float64}`** : Required. A ``\text{n_latents} \times N`` matrix with latents values.
- **`parameters_matrix::Matrix{Float64}`** : Required. A ``(\text{n_latents} + 1) \times I`` matrix with item parameters. intercept (b) must be in first row, latents coefficients ``(a_j)`` in next rows ``(2, \\ldots, n_latents + 1)``.

# Output

A ``I \times N`` `Float64` matrix.
"""
function latent_information(
    latents_matrix::Matrix{Float64},
    parameters_matrix::Matrix{Float64},
)
    p = probability(parameters_matrix, latents_matrix)
    if size(parameters_matrix, 2) > 1
        return _matrix_cols_vec(p, parameters_matrix[:, end], (x, y) -> _p1p(x) * y^2)
    else
        return _p1p.(p)
    end
end

"""
```julia
latent_information_3PL(
    latents_matrix::Matrix{Float64},
    parameters_matrix::Matrix{Float64},
)
```

# Description

Only for models which has guessing parameter (c) in last row of parameters_matrix.
It computes the information function (IIF) for item parameters at latents values provided in matrix form.
It follows the parametrization ``a \theta - b ``.

# Arguments

- **`latents_matrix::Matrix{Float64}`** : Required. A ``\text{n_latents} \times N`` matrix with latents values.
- **`parameters_matrix::Matrix{Float64}`** : Required. A ``(\text{n_latents} + 1) \times I`` matrix with item parameters. intercept (b) must be in first row, latents coefficients (a_j) in next rows (2, ..., n_latents + 1).

# Output

A ``I \times N`` `Float64` matrix.
"""
function latent_information_3PL(
    latents_matrix::Matrix{Float64},
    parameters_matrix::Matrix{Float64},
)
    p = probability_3PL(parameters_matrix, latents_matrix)
    return _matrix_cols_vec(p, parameters_matrix[:, 2], (x, y) -> y^2 * (1 - x) * x) .*
           _matrix_cols_vec(p, parameters_matrix[:, end], (x, y) -> ((x - y) / (1 - y))^2)
end

ext: jl


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5

# Short description of this file: Simple Example Neural PDE

using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
using Plots
#GalacticOptim
# 3D PDE
@parameters x y t θ
@variables u(..)
@derivatives Dxx''~x
@derivatives Dyy''~y
@derivatives Dt'~t

# 3D PDE
eq  = Dt(u(x,y,t,θ)) ~ Dxx(u(x,y,t,θ)) + Dyy(u(x,y,t,θ))
# Initial and boundary conditions
bcs = [u(x,y,0,θ) ~ exp(x+y)*cos(x+y),
       u(0,y,t,θ) ~ exp(y)*cos(y+4t),
       u(2,y,t,θ) ~ exp(2+y)*cos(2+y+4t) ,
       u(x,0,t,θ) ~ exp(x)*cos(x+4t),
       u(x,2,t,θ) ~ exp(x+2)*cos(x+2+4t)]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,2.0),
           y ∈ IntervalDomain(0.0,2.0),
           t ∈ IntervalDomain(0.0,2.0)]

# Discretization
dx = 0.25; dy= 0.25; dt = 0.25
# Neural network
chain = FastChain(FastDense(3,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

discretization = NeuralPDE.PhysicsInformedNN([dx,dy,dt],
                                             chain,
                                             strategy = NeuralPDE.StochasticTraining(include_frac=0.9))

pde_system = PDESystem(eq,bcs,domains,[x,y,t],[u])
prob = NeuralPDE.discretize(pde_system,discretization)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

res = GalacticOptim.solve(prob, ADAM(0.1), progress = false; cb = cb, maxiters=3000)
phi = discretization.phi

ts,xs = [domain.domain.lower:dx:domain.domain.upper for domain in domains]

analytic_sol_func(t,x) = [exp(-t)*sin(pi*x), exp(-t)*cos(pi*x), (1+pi^2)*exp(-t)]
u_real  = [[analytic_sol_func(t,x)[i] for t in ts for x in xs] for i in 1:3]

u_predict  = [[phi([t,x],res.minimizer)[i] for t in ts for x in xs] for i in 1:3]
diff_u = [abs.(u_real[i] .- u_predict[i] ) for i in 1:3]

for i in 1:3
    p1 = plot(xs, ts, u_real[i], st=:surface,title = "u$i, analytic");
    p2 = plot(xs, ts, u_predict[i], st=:surface,title = "predict");
    p3 = plot(xs, ts, diff_u[i],linetype=:contourf,title = "error");
    plot(p1,p2,p3)
    savefig("sol_u$i")
end


text:
@deprecate calibrationerror(estimator, data...) estimator(predictions_targets(data...)...)

@deprecate BiasedSKCE(k1::Kernel, k2::Kernel) BiasedSKCE(k1 ⊗ k2)
@deprecate UnbiasedSKCE(k1::Kernel, k2::Kernel) UnbiasedSKCE(k1 ⊗ k2)
@deprecate BlockUnbiasedSKCE(k1::Kernel, k2::Kernel) BlockUnbiasedSKCE(k1 ⊗ k2)
@deprecate BlockUnbiasedSKCE(k1::Kernel, k2::Kernel, blocksize::Int) BlockUnbiasedSKCE(
    k1 ⊗ k2, blocksize
)
@deprecate UCME(k1::Kernel, k2::Kernel, data...) UCME(k1 ⊗ k2, data...)
@deprecate UCME(kernel::Kernel, testdata...) UCME(
    kernel, predictions_targets(testdata...)...
)
@deprecate UCME(kernel::Kernel, testdata::Tuple{<:Any,<:Any}) UCME(
    kernel, predictions_targets(testdata...)...
)
@deprecate UCME(kernel::Kernel, testdata::AbstractVector{<:Tuple{<:Any,<:Any}}) UCME(
    kernel, predictions_targets(map(first, testdata), map(last, testdata))...
)
@deprecate UCME(
    kernel::Kernel, predictions::AbstractMatrix{<:Real}, targets::AbstractVector
) UCME(kernel, predictions_targets(predictions, targets)...)

@deprecate predictions_targets((predictions, targets)::Tuple{<:Any,<:Any}) predictions_targets(
    predictions, targets
) false
@deprecate predictions_targets(data::AbstractVector{<:Tuple{<:Any,<:Any}}) predictions_targets(
    map(first, data), map(last, data)
) false
@deprecate predictions_targets(predictions::AbstractMatrix{<:Real}, targets::AbstractVector) predictions_targets(
    [predictions[:, i] for i in axes(predictions, 2)], targets
) false
@deprecate predictions_targets(predictions::AbstractVector, targets::AbstractVector) tuple(
    predictions, targets
) false

@deprecate unbiasedskce(kernel, predictions, targets) UnbiasedSKCE(kernel)(
    predictions, targets
) false


text:
localIP =  Sockets.localhost
DG = dg_new(localIP,8888,"QRS")
DG.data["QPoint"].data[2]("0","500","500")
Int32.(DG.data["WidthQRS"].data[1]("0","10",""))

info1 = helloIBox.findnode(tree, "ClassQRS") #запрос группы
mask =helloIBox.getMask(helloIBox.parseType(info1["attrs"]))

a,b,c,d = helloIBox.getSegBegsType(DG)

using StructArrays
newSeg = StructArray(ibeg = [10], iend = [20], type = [1])
helloIBox.addSeg!(DG,newSeg,"simpleAdd")

DG = dg_new(localIP,8888,"QRS")
allObj = Dict{String,Any}()
allObj["history"] = []
allObj["state"] = 0
allObj["dataStorage"] = Dict{String,Any}()
allObj["dataStorage"]["QRS"] = DG
command1Path = joinpath(Base.@__DIR__, "files","oxy115829.002","command_FT.json")
command_FT= String(read(command1Path))
helloIBox.parseCommand(localIP, 8888, allObj,command_FT)
@test allObj["dataStorage"]["QRS"].result.QPoint[1:2]==[10,30]
@test allObj["dataStorage"]["QRS"].result.WidthQRS[1:2]==[10,15]

#добавление с расширением границ, если тип одинаковый
command1Path = joinpath(Base.@__DIR__, "files","oxy115829.002","command_FT_addRaw.json")
command_FT= String(read(command1Path))
helloIBox.parseCommand(localIP, 8888, allObj,command_FT)
@test allObj["dataStorage"]["QRS"].result.QPoint[1:2]==[10,50]
@test allObj["dataStorage"]["QRS"].result.WidthQRS[1:2]==[35,5]
ibeg =  DG.data[DG.ibegdata].data[1]("0","","")
#удаление сегента
command1Path = joinpath(Base.@__DIR__, "files","oxy115829.002","command_FT_delete.json")
command_FT= String(read(command1Path))
helloIBox.parseCommand(localIP, 8888, allObj,command_FT)
@test allObj["dataStorage"]["QRS"].result.QPoint[1:2]==[101,221]
@test allObj["dataStorage"]["QRS"].result.WidthQRS[1:2]==[23,21]


text:
using Revise
using AdFem
using LinearAlgebra
using SparseArrays

A = rand(10, 10)
bd = rand(1:10, 5)
bdval = rand(5)
rhs = rand(10)
A1, rhs1 = impose_Dirichlet_boundary_conditions(A, rhs, bd, bdval)
A2, rhs2 = impose_Dirichlet_boundary_conditions(constant(sparse(A)), rhs, bd, bdval)

sess = Session(); init(sess)
@info maximum(abs.(Array(run(sess, A2)) - A1))
@info maximum(abs.(run(sess, rhs2) - rhs1))


text:
const _club_base = 0x1F0D1
const _diamond_base = 0x1F0C1
const _heart_base = 0x1F0B1
const _spade_base = 0x1F0A1

function _make_chars(base::Integer)
    vals = collect(base:base+10)
    append!(vals, base + 12)
    append!(vals, base + 13)
    return Char.(vals)

end

_cards = _make_chars(_club_base)
append!(_cards, _make_chars(_diamond_base))
append!(_cards, _make_chars(_heart_base))
append!(_cards, _make_chars(_spade_base))

import Base.Char

"""
    Char(c::Card)

Return the unicode character for that playing card.

```
julia> Char(4♣)
'🃔': Unicode U+1F0D4 (category So: Symbol, other)
```
"""
function Char(c::Card)
    return @inbounds _cards[index(c)]
end


text:
using GeoClustering
using Documenter

DocMeta.setdocmeta!(GeoClustering, :DocTestSetup, :(using GeoClustering); recursive=true)

makedocs(;
    modules=[GeoClustering],
    authors="RohitRathore1",
    repo="https://github.com/RohitRathore1/GeoClustering.jl/blob/{commit}{path}#{line}",
    sitename="GeoClustering.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://RohitRathore1.github.io/GeoClustering.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/RohitRathore1/GeoClustering.jl",
)


text:
using IterativeSolvers
using FactCheck

facts("svdvals_tr") do
#Thick restart methods
for method in (:thick, :harmonic) context("Thick restart with method=$method") do
  for elty in (Float32, Float64)
    context("Diagonal Matrix{$elty}") do
        n = 30
        A = full(Diagonal(elty[1.0:n;]))
        q = convert(Vector{elty}, ones(n)/√n)
        σ, L = svdvals_tr(A, q, 5, 10, tol=1e-5, maxiter=30, method=method)
        @fact norm(σ - [n:-1.0:n-4;]) --> less_than(5^2*1e-5)
    end

    context("Rectangular Matrix{$elty}") do
        srand(1)
        m = 300
        n = 200
        k = 5
        l = 10

        A = convert(Matrix{elty}, randn(m,n))
        q = convert(Vector{elty}, randn(n))|>x->x/norm(x)
        σ, L = svdvals_tr(A, q, k, l, tol=1e-5, maxiter=30, method=method)
        @fact norm(σ - svdvals(A)[1:k]) --> less_than(k^2*1e-5)
    end
  end
end end
end #svdvals_tr

facts("BrokenArrowBidiagonal") do
    B = IterativeSolvers.BrokenArrowBidiagonal([1, 2, 3], [1, 2], Int[])
    @fact full(B) --> [1 0 1; 0 2 2; 0 0 3]
end


text:
# Temperature Modulation
const GLOBALKERNEL_ = Function[kse]

"Global Kernel"
function globalkernel!(k)
  global GLOBALKERNEL_
  GLOBALKERNEL_[1] = k
end

"Retrieve global kernel"
function globalkernel()
  global GLOBALKERNEL_
  GLOBALKERNEL_[1]
end

"Temporarily set global kernel"
function withkernel(thunk, k)
  globalkernel!(k)
  res = thunk()
  globalkernel!(kse)
  res
end

# Ctx Based Temperature Modulation

"Context for particular kernel "
struct KernelContext{K}
  kernel::K
end

"Evaluate f(args...) under kernels kc"
function ctxapl(kc::KernelContext, f, args...)
  old = getctx(:kernel)  # FIXME: What if kernel doesn't exist
  setctx!(:kernel, kc.kernel)
  res = f(args...)
  setctx!(:kernel, old)
  res
end

# function testkernelctx()
#   softgt(0.3, 0.2, kernel = globalkernel())
#   softgt(0.3, 0.2, kernel = getctx(:kernel))
# end

# (Cassette-based) Temperature Modulation
Cassette.@context AlphaCtx

"""
`f(args)` where temperature controlled with temperature `α`

```julia
x = normal(0.0, 1.0)
atα(10, rand, x ==ₛ 0.3)
```
"""
@inline function atα(α, f, args...)
  ctx = Cassette.disablehooks(AlphaCtx(metadata = α))
  Cassette.overdub(ctx, f, args...)
end

# @inline Cassette.overdub(ctx::AlphaCtx, ::typeof(kse), x, α) = kse(x, ctx.metadata)
@inline Cassette.overdub(ctx::AlphaCtx, ::typeof(globalkernel)) = kseα(ctx.metadata)

"""

```julia
x = normal(0.0, 1.0)
@atα 100 rand(y ==ₛ 0.0)
```
"""
macro atα(a, fexpr)
  :(atα($(esc(a)), $(esc.(fexpr.args)...)))
end
@pre fexpr.head == :call


text:
using .util_convexopt

"project a vector 'u' to l1 unit ball "
function proj_l1_ball!(u::Array{T, 1}) where {T<:AbstractFloat}
    sum = Inf
    shrink = 0.0
    while (sum > 1.0)
        sum = 0.0
        cnt = 0

        for i=1:length(u)
            u[i] = max(u[i] - shrink, 0.0)
            sum += abs(u[i])
            if u[i] != 0.0
                cnt += 1
            end
        end
        if cnt > 0
            shrink = (sum - 1.0) / cnt
        else
            break
        end
    end
end

"""
    function prox_nuclear!(y, x, stepsize)

Proximal operator of nuclear norm to x. y=prox_{s f}(x)

is_S1=true : (S1, l1) norm for TNV
is_S1=false: (S∞, l1)
"""
function prox_nuclear_wo_svd!(y, x, stepsize, is_S1=true)
    H, W, M, C = size(x) # Dim : 2 for 2D, 3 for 3D (C >= 2)
    x_HxWxCx2 = PermutedDimsArray(x, [1,2,4,3])
    y_HxWxCx2 = PermutedDimsArray(y, [1,2,4,3])
    x_ = reshape(x_HxWxCx2, H*W, C, 2)
    y_ = reshape(y_HxWxCx2, H*W, C, 2)

    EPS = 1e-8

    # x_svd = mapslices(svd!, x_, dims=[2,3])
    # x_svd = dropdims(x_svd_, dims=4)
    Threads.@threads for hw=1:H*W
        x_hw = view(x_, hw, :, :) # [C x 2]
        y_hw = view(y_, hw, :, :)

        BB1 = 0.0
        BB2 = 0.0
        BB4 = 0.0

        for k=1:C
            BB1 += x_hw[k, 1] * x_hw[k, 1]
            BB2 += x_hw[k, 1] * x_hw[k, 2]
            BB4 += x_hw[k, 2] * x_hw[k, 2]
        end

        # # eigenvalues of BB
        T = BB1 + BB4
        D = BB1 * BB4 - BB2 * BB2

        det = sqrt( max( ( T*T / 4.0) - D, 0.0) )

        eig1 = max(  (T / 2.0) + det , 0.0 )
        eig2 = max(  (T / 2.0) - det , 0.0 )
        σ1 = sqrt(eig1)
        σ2 = sqrt(eig2)

        if abs(BB2) > EPS
            v0 = BB2
            V1 = eig1 - BB4
            V2 = eig2 - BB4

            len1 = sqrt(v0*v0 + V1*V1)
            len2 = sqrt(v0*v0 + V2*V2)

            if len1 > EPS
                V1 /= len1
                V3 = v0 / len1
            end
            if len2 > EPS
                V2 /= len2
                V4 = v0 / len2
            end
        else
            if BB1 > BB4
                V1 = V4 = 1.0
                V2 = V3 = 0.0
            else
                V1 = V4 = 0.0
                V2 = V3 = 1.0
            end
        end

        if is_S1
            # (S1, l1) Eq. (11) in IPOL paper
            ss1 = max(σ1 - stepsize, 0.0)
            ss2 = max(σ2 - stepsize, 0.0)
        else
            # (S∞, l1) prox to l∞
            # l∞(x) = x - λ prox_(*/λ) (x / λ)
            # proj_dual_l1!(σ1/sigma, 1.0/sigma)
            # proj_dual_l1: x ./= max.(1.0, abs.(x) / weight)
            pp = [σ1 / stepsize, σ2 / stepsize]
            proj_l1_ball!(pp)

            ss1 = σ1 - stepsize * pp[1]
            ss2 = σ2 - stepsize * pp[2]
        end

        if σ1 > EPS
            ss1 /= σ1
        end
        if σ2 > EPS
            ss2 /= σ2
        end

        t1 = ss1 * V1 * V1 + ss2 * V2 * V2
        t2 = ss1 * V1 * V3 + ss2 * V2 * V4
        t3 = ss1 * V3 * V3 + ss2 * V4 * V4

        # B * t
        for c=1:C
            y_hw[c,1] = x_hw[c,1]*t1 + x_hw[c,2]*t2
            y_hw[c,2] = x_hw[c,1]*t2 + x_hw[c,2]*t3
        end
    end
    return y
end

"""
Example:
x = rand(20, 30, 2, 10);
y = zeros(20, 30, 2, 10);
prox_l∞11!(y, x, 0.01);
"""
function prox_l∞11!(y, x, stepsize)
    H, W, M, C = size(x) # Dim : 2 for 2D, 3 for 3D (C >= 2)
    x_HxWxCx2 = PermutedDimsArray(x, [1,2,4,3])
    y_HxWxCx2 = PermutedDimsArray(y, [1,2,4,3])
    x_ = reshape(x_HxWxCx2, H*W, C, 2)
    y_ = reshape(y_HxWxCx2, H*W, C, 2)

    # x_svd = mapslices(svd!, x_, dims=[2,3])
    # x_svd = dropdims(x_svd_, dims=4)
    Threads.@threads for hw=1:H*W
        x_hw = view(x_, hw, :, :) # [C x 2]
        y_hw = view(y_, hw, :, :)

        pp1 = zeros(C)
        pp2 = zeros(C)
        for c=1:C
            pp1[c] = abs(x_hw[c,1]) / stepsize
            pp2[c] = abs(x_hw[c,2]) / stepsize
        end

        proj_l1_ball!(pp1)
        proj_l1_ball!(pp2)

        for c=1:C
            x1 = x_hw[c,1]
            x2 = x_hw[c,2]

            y_hw[c,1] = x1 - stepsize * sign(x1) * pp1[c]
            y_hw[c,2] = x2 - stepsize * sign(x2) * pp2[c]
        end
    end
end

function grad2d!(du, u) where {T<:AbstractFloat}
    du1, du2 = view(du, :, :, 1), view(du, :, :, 2)
    du1 .= circshift!(du1, u, [0, -1]) .- u
    du2 .= circshift!(du2, u, [-1, 0]) .- u
    du1[:, end] .= 0.0
    du2[end, :] .= 0.0
end

function div2d!(divp, p) where {T<:AbstractFloat}
    p1, p2 = view(p,:,:,1), view(p,:,:,2)
    p1_x, p2_y = view(divp,:,:,1), view(divp,:,:,2)

    # p2_y is temp variable to save memory
    p1_x         .=  p1 - circshift!(p2_y, p1, [0, 1])
    p1_x[:, end] .= -p1[:, end-1]
    p1_x[:,   1] .=  p1[:, 1]

    circshift!(p2_y, p2, [1, 0])
    p2_y         .= p2 - circshift!(p2_y, p2, [1, 0])
    p2_y[end, :] .= -p2[end-1, :]
    p2_y[1,   :] .=  p2[1, :]

    p1_x .+= p2_y
    return p1_x
end

function _recon2d_ctv_primaldual!(u::Array{T, 3}, A, b0::Array{T, 3}, niter, w_data, sigmas, tau, type, ϵ, nverbose) where {T<:AbstractFloat}
    At = A'
    H, W, C = size(u)

    # normalize the data fidelity with respect.
    # we consider mean(residual) instead of sum
    w_data = w_data / C
    #w_data = w_data / (C  )

    b = reshape(b0, (size(b0, 1)*size(b0, 2), C))
    # b = vec(b0)

    ubar = deepcopy(u)
    u_prev = similar(u)

    p1 = zeros(size(b))
    p2 = zeros(H, W, 2, C) # q in Duran 2016
    g  = similar(p2)
    p2_temp = similar(p2)

    p_adjoint = zeros(H, W, C)

    du = similar(p2)
    divp2 = similar(p2)

    data1 = copy(b)

    invσ11 = 1. / (sigmas[1] / w_data + 1.0) # l2 data, (Handa Sec 8.1)
    invσ2 = 1. / sigmas[2]

    res_primals = zeros(niter)
    res_duals = zeros(niter)

    if nverbose > 0
        p_adjoint_prev = similar(p_adjoint)
    #    if flag_dual_residual == true
            # variables for computinig residuals

            p1_prev = similar(p1)
            p2_prev = similar(p2)
            du_prev = similar(du)
            data1_prev = similar(data1)
        # end
    end

    # init p1
    # Threads.@threads for c=1:C
    #     @views mul!(data1[:,c], A, vec(ubar[:,:,c]))
    #     p1[:,c] .= data1[:,c]
    #     @views data1[:,c] .-= b[:,c]
    # end

    for it=1:niter
        if nverbose > 0 && it % nverbose == 0
            copy!(p_adjoint_prev, p_adjoint)
            # if flag_dual_residual == true
                copy!(data1_prev, data1)
                copy!(p1_prev, p1)
                copy!(p2_prev, p2)
                copy!(du_prev, du)
            # end
        end

        Threads.@threads for c=1:C
            ubar_c = view(ubar, :, :, c)
            data1_c = view(data1, :, c)
            p_adjoint_c = vec(view(p_adjoint, :, :, c))
            p1_c = view(p1, :, c)

            # dual update: data fidelity
            data1_c .= mul!(data1_c, A, vec(ubar_c)) .- view(b, :, c)
            p1_c .= (p1_c .+ sigmas[1] .* data1_c) .* invσ11
            mul!(p_adjoint_c, At, p1_c)

            # compute gradient for TNV later
            du_c = view(du, :, :, :, c)
            grad2d!(du_c, ubar_c)
        end

        u_prev .= u

        # projection onto dual of (S1,l1)
        p2_temp .= du .+ invσ2 .* p2

        if type == "tnv" || type == "S1l1"
            prox_nuclear_wo_svd!(g, p2_temp, invσ2)
        elseif type == "S∞l1"
            prox_nuclear_wo_svd!(g, p2_temp, invσ2, false)
        elseif type == "l∞11" # strong coupling
            prox_l∞11!(g, p2_temp, invσ2)
        else
            error("not supported type")
        end

        # prox_nuclear!(g, p2_temp, invσ2)
        p2 .+= sigmas[2] .* ( du .- g )

        Threads.@threads for c=1:C
            p_adjoint_c = view(p_adjoint, :, :, c)
            p_adjoint_c .-= div2d!(view(p2_temp,:,:,:,c), view(p2,:,:,:,c))
        end

        # primal update
        # u .= u .- tau .* p_adjoint
        u .= max.(u .- tau .* p_adjoint, 0.0) # positivity constraint

        # acceleration
        ubar .= 2 .* u .- u_prev

        # primal residual
        # if it > 1
        # end

        if nverbose > 0 && it % nverbose == 0
            # compute primal energy (optional)
            res_primal = sum(abs.((u_prev - u) / tau .- (p_adjoint_prev .- p_adjoint))) / length(u)
            res_primals[it] = res_primal

            residual = res_primal
            res_dual = 0.0

            # if flag_dual_residual == true
                res_dual1 = sum(abs.((p1_prev-p1) / sigmas[1] .- (data1_prev .- data1)))
                res_dual2 = sum(abs.((p2_prev-p2) / sigmas[2] .- (du_prev .- du)))
                res_dual = (res_dual1 + res_dual2) / (length(p1) + length(p2))
                res_duals[it] = res_dual
            # end
            residual += res_dual

            if residual < ϵ && it > 1
                @info "$it Stopping condition is met. $res_primal $res_dual"
                return (u, it, res_primals, res_duals)
            end

            energy = sum(data1.^2) /  length(data1)
            println("$it, data term: $energy, primal_res: $res_primal, dual_res: $res_dual")
        end
    end

    return (u, niter, res_primals, res_duals)
end

"""
    recon2d_ctv_primaldual!(u::Array{T, 2}, A, b::Array{T, 2}, niter::Int, w_tv::T, type="tnv", c=1.0)

Reconstruct a 2d image by Collaborative TV using Primal Dual optimization optimizer

# Args
u : Initial guess of images
A : Forward opeartor
b : Projection data
niter : number of iterations
w_tv : weight for TV term
type (string) : tnv, S∞l1, l∞11
c : See 61 page in 2016_Chambolle,Pock_An_introduction_to_continuous_optimization_for_imagingActa_Numerica

For Collaborative TV, refer to:
Duran,Moeller,Sbert,Cremers_On_the_Implementation_of_Collaborative_TV_Regularization_-_Application_toImage_Processing_On_Line
"""
function recon2d_ctv_primaldual!(u::Array{T, 3}, A, b::Array{T, 3}, niter::Int, w_data, type="tnv", ϵ=1e-6, nverbose=10, c=1.0) where {T <: AbstractFloat}
    if size(u, 3) != size(b, 3)
        error("The channel size of u and b should match.")
    end

    op_A_norm = compute_opnorm(A, 6) # for safety
    println("@ opnorm of forward projection operator: $op_A_norm")
    ops_norm = [op_A_norm, sqrt(8)]

    sigmas = zeros(length(ops_norm))
    for i=1:length(ops_norm)
        sigmas[i] = 1.0 / (ops_norm[i] * c)
    end

    tau = c / sum(ops_norm)

    sigmas .*= 0.99
    tau *= 0.99

    println("@ step sizes sigmas: ", sigmas, ", tau: $tau")

    return _recon2d_ctv_primaldual!(u, A, b, niter, w_data, sigmas, tau, type, ϵ, nverbose)
end

max_stars_repo_name: QPImaging/xfromProjections

max_issues_repo_name: QPImaging/xfromProjections

max_forks_repo_name: QPImaging/xfromProjections


text:
# Return ODESystem pretty printing for our wrapper types
#Base.show(io::IO, ::MIME"text/plain", x::IonChannel) = Base.display(isbuilt(x) ? x.sys : x)


text:
# QuEST_jl/src/base/QASM_logging.jl
#

function clearRecordedQASM(qureg ::QuEST_Types.Qureg) ::Nothing
    ccall(:clearRecordedQASM, Cvoid, (QuEST_Types.Qureg,), qureg)
    return nothing
end

function printRecordedQASM(qureg ::QuEST_Types.Qureg) ::Nothing
    ccall(:printRecordedQASM, Cvoid, (QuEST_Types.Qureg,), qureg)
    return nothing
end

function startRecordingQASM(qureg ::QuEST_Types.Qureg) ::Nothing
    ccall(:startRecordingQASM, Cvoid, (QuEST_Types.Qureg,), qureg)
    return nothing
end

function stopRecordingQASM(qureg ::QuEST_Types.Qureg) ::Nothing
    ccall(:stopRecordingQASM, Cvoid, (QuEST_Types.Qureg,), qureg)
    return nothing
end

function writeRecordedQASMToFile(qureg    ::QuEST_Types.Qureg,
                                 filename ::String) ::Nothing
    @assert begin
        ios = open(filename, "w")
        close(ios) === nothing
    end

    ccall(:writeRecordedQASMToFile, Cvoid, (QuEST_Types.Qureg, Cstring), qureg, filename)
    return nothing
end


text:
using Hecke

boundexp = ARGS[1]
c4fields = ARGS[2]
startfield = ARGS[3]
number = ARGS[4]

sprint_formatted(fmt, args...) = @eval @sprintf($fmt, $(args...))

if length(ARGS) == 5
  basename = ARGS[5]
else
  basename = ""
end

boundexp = parse(Int, boundexp)
startfield = parse(Int, startfield)
number = parse(Int, number)

@assert mod(boundexp, 3) == 0

boundquad = div(boundexp, 3)

lall = Hecke._read_fields(c4fields)

l = lall[startfield:(startfield + number - 1)]

width = length(string(length(lall)))

for (i, f) in enumerate(l)
  K, _ = NumberField(f[1], "a")
  @show K
  z = Hecke.Dic3_extensions(fmpz(10)^boundexp, K)
  @show length(z)
  if length(z) != 0
    fname = basename * sprint_formatted("%0$(width)d", startfield + i - 1)
    @show fname
    file = open(fname, "a")
    write(file, "# $(K.pol)\n")
    for (p, D) in z
      write(file, "($p, $D)\n")
    end
    close(file)
  end
end


text:
# vim: set et ts=2 sw=2:
@@with_kw mutable struct Table
  ys=[]; xs=[]; rows=[]; cols=[] end

@with_kw mutable struct Row
  cells=[]; dom=0; y=no end

goalp(s,c=it.char) = c.less in s || c.more in s || c.klass in s

max_stars_repo_name: timm/jut

max_issues_repo_name: timm/jut

max_forks_repo_name: timm/jut


text:
using Test
using Distributed
using LinearAlgebra
using NPZ

@everywhere using Random
@everywhere using DPMMSubClusters
include("multinomial_tests.jl")
include("niw_tests.jl")
include("unitests.jl")
include("module_tests.jl")

addprocs(2)
@everywhere using Random
@everywhere using DPMMSubClusters

include("multinomial_tests.jl")
include("niw_tests.jl")
include("unitests.jl")
include("module_tests.jl")


text:
#  Functions and Constants required for traversal

const DNAPAIRS = [(DNA_A, DNA_T),   (DNA_T, DNA_A),
                  (DNA_G, DNA_C),   (DNA_C, DNA_G),
                  (DNA_G, DNA_T),   (DNA_T, DNA_G)]

const DNAGAPS  = [(DNA_A, DNA_Gap), (DNA_Gap, DNA_A),
                  (DNA_C, DNA_Gap), (DNA_Gap, DNA_C),
                  (DNA_G, DNA_Gap), (DNA_Gap, DNA_G),
                  (DNA_T, DNA_Gap), (DNA_Gap, DNA_T)]

const DNAMISMATCH = [(DNA_A, DNA_A), (DNA_A, DNA_C),
                     (DNA_A, DNA_G), (DNA_C, DNA_A),
                     (DNA_C, DNA_C), (DNA_C, DNA_T),
                     (DNA_G, DNA_A), (DNA_G, DNA_G),
                     (DNA_T, DNA_C), (DNA_T, DNA_T)]

const RNAPAIRS = [(RNA_A, RNA_U), (RNA_U, RNA_A),
                  (RNA_G, RNA_C), (RNA_C, RNA_G),
                  (RNA_G, RNA_U), (RNA_U, RNA_G)]

const RNAGAPS  = [(RNA_A, RNA_Gap), (RNA_Gap, RNA_A),
                  (RNA_C, RNA_Gap), (RNA_Gap, RNA_C),
                  (RNA_G, RNA_Gap), (RNA_Gap, RNA_G),
                  (RNA_U, RNA_Gap), (RNA_Gap, RNA_U)]

const RNAMISMATCH = [(RNA_A, RNA_A), (RNA_A, RNA_C),
                     (RNA_A, RNA_G), (RNA_C, RNA_A),
                     (RNA_C, RNA_C), (RNA_C, RNA_U),
                     (RNA_G, RNA_A), (RNA_G, RNA_G),
                     (RNA_U, RNA_C), (RNA_U, RNA_U)]

const PairsType = Vector{Tuple{DNA, DNA}}

function twobitalpha{S <: BioSequences.Sequence}( seq::S )
   t = eltype(seq)
   if t == DNANucleotide
      return DNAAlphabet{2}
   elseif t == RNANucleotide
      return RNAAlphabet{2}
   else
      error("RNASuffixArrays are only built for DNA/RNA input sequences!")
   end
end

pairs{n}(::Type{BioSequences.DNAAlphabet{n}})      = DNAPAIRS
pairs{n}(::Type{BioSequences.RNAAlphabet{n}})      = RNAPAIRS
gaps{n}(::Type{BioSequences.DNAAlphabet{n}})       = DNAGAPS
gaps{n}(::Type{BioSequences.RNAAlphabet{n}})       = RNAGAPS
mismatches{n}(::Type{BioSequences.DNAAlphabet{n}}) = DNAMISMATCH
mismatches{n}(::Type{BioSequences.RNAAlphabet{n}}) = RNAMISMATCH

gap{n}(::Type{DNAAlphabet{n}})                = DNA_Gap
gap{n}(::Type{RNAAlphabet{n}})                = RNA_Gap

onehot{I <: Integer}(x::I) = 0x01 << (x-1)

encodeindex{A <: Alphabet, N <: NucleicAcid}(::Type{A}, x::N) = (x == DNA_Gap || x == RNA_Gap) ? 0 : BioSequences.encode(A, x)+1
index{A <: Alphabet}(::Type{A}, func=pairs) = map( x->map(y->encodeindex(A, y), x), func(A) )

revoffset( x::Int, seq::BioSequence ) = revoffset( x, length(seq) )
revoffset( x, len ) = len - x + 1

Base.reverse( seq::BioSequences.ReferenceSequence ) = BioSequences.ReferenceSequence( reverse( String( seq ) ) )


text:
"""
    A, B, C, D = partitionABCD(ABCD, m)

Partition `ABCD` into `A`, `B`, `C`, `D` for an `m`-input state-space system.
"""
function partitionABCD(ABCD, m=NaN)
    if isnan(m)
        n = minimum(size(ABCD)) - 1
        m = size(ABCD, 2) - n
    else
        n = size(ABCD, 2) - m
    end
    r = size(ABCD, 1) - n

    A = ABCD[1:n, 1:n]
    B = ABCD[1:n, n+1:n+m]
    C = ABCD[n+1:n+r, 1:n]
    D = ABCD[n+1:n+r, n+1:n+m]

    return A, B, C, D
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
using Test
using NiLang, NiLang.AD
using LinearAlgebra
using Random

@testset "inv" begin
    Random.seed!(2)

    @test check_inv(i_inv!, (randn(3, 3), randn(3, 3)))
    @test check_inv(⊕(det), (0.3, randn(3, 3)))
    @test check_inv(⊕(logdet), (0.3, randn(3, 3)))
    @test check_grad(⊕(det), (0.3, randn(3, 3)), iloss=1)
    @test check_grad(⊕(logdet), (0.3, randn(3, 3)), iloss=1)

    @i function loss(out!, y, A)
        i_inv!(y, A)
        out! += y[1,1]
    end
    @test check_grad(loss, (0.0, randn(3, 3), randn(3, 3)); iloss=1)
end

@testset "affine" begin
    Random.seed!(2)
    A = randn(5, 5)
    b = randn(5)
    x = randn(5)
    y! = zeros(5)
    @test i_affine!(y!, A, b, x)[1] ≈ A*x + b
end


text:
export  PlaneStrainViscoelasticityProny, PlaneStressViscoelasticityProny, ViscoelasticityProny

@doc raw"""
    mutable struct ViscoelasticityProny
        ρ::Float64 # density
        E::Float64 # Young's Modulus
        ν::Float64 # Poisson's ratio
        τ::Float64 # relaxation time
        c::Float64 # coefficient
        planestress::Bool
        g0::Float64
        g1::Float64
        J::Float64
        σ0::Array{Float64} # stress at last time step
        σ0_::Array{Float64} # σ0 to be updated in `commitHistory`
        ε0::Array{Float64} # strain at last time step
        ε0_::Array{Float64} # ε0 to be updated in `commitHistory`
    end

Prony series viscoelasticity model.

$$E(t) = E_0 - cE_0(1-\exp(-t/\tau))$$

In principle, we can have different Prony series for both shear modulus $G$ and bulk modulus $K$,

$$\sigma_{ij} = 3K \left( \frac{1}{3} \epsilon_{kk}\delta_{ij}\right) + 2G \left( \epsilon_{ij} - \frac{1}{3}\epsilon_{kk} \delta_{ij} \right)$$

We only implemented $E(t)$ case for simplicity.
"""
mutable struct ViscoelasticityProny
    ρ::Float64 # density
    E::Float64 # Young's Modulus
    ν::Float64 # Poisson's ratio
    τ::Float64 # relaxation time
    c::Float64 # coefficient
    planestress::Bool
    g0::Float64
    g1::Float64
    J::Array{Float64,1}
    J_::Array{Float64,1}
    H::Array{Float64, 2}
    σ0::Array{Float64} # stress at last time step
    σ0_::Array{Float64} # σ0 to be updated in `commitHistory`
    ε0::Array{Float64} # strain at last time step
    ε0_::Array{Float64} # ε0 to be updated in `commitHistory`
end

function ViscoelasticityProny(prop::Dict{String, Any})
    ρ = Float64(prop["rho"])
    E = Float64(prop["E"])
    ν = Float64(prop["nu"])
    τ = Float64(prop["tau"])
    c = Float64(prop["c"])
    planestress = Bool(prop["planestress"])
    H = zeros(3,3)
    if planestress
        H[1,1] = E/(1. -ν*ν)
        H[1,2] = H[1,1]*ν
        H[2,1] = H[1,2]
        H[2,2] = H[1,1]
        H[3,3] = E/(2.0*(1.0+ν))
    else
        H[1,1] = E*(1. -ν)/((1+ν)*(1. -2. *ν));
        H[1,2] = H[1,1]*ν/(1-ν);
        H[2,1] = H[1,2];
        H[2,2] = H[1,1];
        H[3,3] = H[1,1]*0.5*(1. -2. *ν)/(1. -ν);
    end

    planestress = Bool(prop["planestress"])
    σ0 = zeros(3); σ0_ = zeros(3)
    ε0 = zeros(3); ε0_ = zeros(3)
    J = zeros(3);  J_ = zeros(3)
    G0 = E/(2(1+ν))
    g0 = G0*(1-c)
    g1 = G0*c
    ViscoelasticityProny(ρ, E, ν, τ, c, planestress, g0, g1, J, J_, H, σ0, σ0_, ε0, ε0_)
end

function getStress(mat::ViscoelasticityProny,  strain::Array{Float64},  Dstrain::Array{Float64}, Δt::Float64)
    τ = mat.τ
    g0, g1, J = mat.g0, mat.g1, mat.J
    σ = mat.H * strain - 2*g0*Dstrain + 2*g1*(- Dstrain) + exp(-Δt/τ)*J
    J = exp(-Δt/τ) * J + 2*g1*(strain - Dstrain)
    dΔσdΔε = mat.H
    mat.σ0_ = σ
    mat.ε0_ = strain
    mat.J_ = J
    return σ, dΔσdΔε
end

function getTangent(mat::ViscoelasticityProny)
    error("Not implemented")
end

function commitHistory(mat::ViscoelasticityProny)
    mat.σ0 = mat.σ0_
    mat.ε0 = mat.ε0_
    mat.J = mat.J_
end

"""
    PlaneStrainViscoelasticityProny(prop::Dict{String, Any})

Returns the plane strain viscoelasticity Prony shear series material model.
"""
function PlaneStrainViscoelasticityProny(prop::Dict{String, Any})
    prop["planestress"] = false
    ViscoelasticityProny(prop)
end

"""
    PlaneStrainViscoelasticityProny(prop::Dict{String, Any})

Returns the plane stress viscoelasticity Prony shear series material model.
"""
function PlaneStressViscoelasticityProny(prop::Dict{String, Any})
    prop["planestress"] = true
    ViscoelasticityProny(prop)
end


text:
module Currencies

# Data obtained from https://datahub.io/core/country-codes
# Consider downloading data in build.jl?

using DelimitedFiles

export Currency

const (data,headers) = readdlm(joinpath(@__DIR__,"data","country-codes.csv"),',',header=true)

struct Currency{T}
    function Currency{T}() where T
        c = new()
        list[T] = c
        return c
    end
end
const list = Dict{Symbol,Currency}()

const (nrow,ncol) = size(data)
for i in 1:nrow
    currencies = split(data[i,10],",")
    currency_units = split(string(data[i,12]),",")
    currency_names = split(data[i,13],",")
    currency_codes = split(string(data[i,14]),",")
    for (currency,currency_unit,currency_name,currency_code) in zip(currencies,currency_units,currency_names,currency_codes)
        if (length(currency) == 3) & !isdefined(Currencies,Symbol(currency))
            @eval Currencies begin
                $(Symbol(currency)) = Currency{Symbol($(currency))}()
                unit(::Currency{Symbol($(currency))}) = parse(Int,$(currency_unit))
                name(::Currency{Symbol($(currency))}) = $(currency_name)
                code(::Currency{Symbol($(currency))}) = parse(Int,$(currency_code))
                Base.show(io::Base.IO,::Currency{Symbol($(currency))}) = print(io,$(currency))
            end
        end
    end
end

end


text:
#### Land sources
export PhaseChange

function heaviside(x::FT) where {FT}
    if x >= FT(0)
        output = FT(1)
    else
        output = FT(0)
    end
    return output
end

abstract type SoilSource{FT <: AbstractFloat} end

"""
    PhaseChange <: SoilSource
The function which computes the freeze/thaw source term for Richard's equation,
assuming the timescale is the maximum of the thermal timescale and the timestep.
"""
Base.@kwdef struct PhaseChange{FT} <: SoilSource{FT}
    "Timestep"
    Δt::FT = FT(NaN)
    "Timescale for temperature changes"
    τLTE::FT = FT(NaN)
end

function land_source!(
    f::Function,
    land::LandModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    f(land, source, state, diffusive, aux, t, direction)
end

function land_source!(
    source_type::PhaseChange,
    land::LandModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    FT = eltype(state)

    _ρliq = FT(ρ_cloud_liq(land.param_set))
    _ρice = FT(ρ_cloud_ice(land.param_set))
    _Tfreeze = FT(T_freeze(land.param_set))
    _LH_f0 = FT(LH_f0(land.param_set))
    _g = FT(grav(land.param_set))

    ϑ_l, θ_i = get_water_content(land.soil.water, aux, state, t)
    eff_porosity = land.soil.param_functions.porosity - θ_i
    θ_l = volumetric_liquid_fraction(ϑ_l, eff_porosity)
    T = get_temperature(land.soil.heat, aux, t)

    # The inverse matric potential is only defined for negative arguments.
    # This is calculated even when T > _Tfreeze, so to be safe,
    # take absolute value and then pick the sign to be negative.
    # But below, use heaviside function to only allow freezing in the
    # appropriate temperature range (T < _Tfreeze)
    ψ = -abs(_LH_f0 / _g / _Tfreeze * (T - _Tfreeze))
    hydraulics = land.soil.water.hydraulics
    θstar =
        land.soil.param_functions.porosity *
        inverse_matric_potential(hydraulics, ψ)

    τft = max(source_type.Δt, source_type.τLTE)
    freeze_thaw =
        FT(1) / τft * (
            _ρliq *
            (θ_l - θstar) *
            heaviside(_Tfreeze - T) *
            heaviside(θ_l - θstar) - _ρice * θ_i * heaviside(T - _Tfreeze)
        )
    source.soil.water.ϑ_l -= freeze_thaw / _ρliq
    source.soil.water.θ_i += freeze_thaw / _ρice
end

# sources are applied additively

@generated function land_source!(
    stuple::Tuple,
    land::LandModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    N = fieldcount(stuple)
    return quote
        Base.Cartesian.@nexprs $N i -> land_source!(
            stuple[i],
            land,
            source,
            state,
            diffusive,
            aux,
            t,
            direction,
        )
        return nothing
    end
end

max_forks_repo_name: jm-c/CLIMA


text:
# This file is a part of Julia, but is derived from
# https://github.com/floitsch/double-conversion which has the following license
#
# Copyright 2006-2014, the V8 project authors. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
#     * Redistributions of source code must retain the above copyright
#       notice, this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above
#       copyright notice, this list of conditions and the following
#       disclaimer in the documentation and/or other materials provided
#       with the distribution.
#     * Neither the name of Google Inc. nor the names of its
#       contributors may be used to endorse or promote products derived
#       from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

const kMinExp = -60
const kMaxExp = -32

function roundweed(buffer,len,rest,tk,unit,kappa,too_high::UInt64,unsafe_interval::UInt64)
    small = too_high - unit
    big = too_high + unit
    while rest < small &&
            unsafe_interval - rest >= tk &&
            (rest + tk < small ||
            small - rest >= rest + tk - small)
        buffer[len-1] -= 1
        rest += tk
    end
    if rest < big &&
        unsafe_interval - rest >= tk &&
        (rest + tk < big ||
        big - rest > rest + tk - big)
        return false, kappa
    end
    return (2 * unit <= rest) && (rest <= unsafe_interval - 4 * unit), kappa
end

const SmallPowersOfTen = [
        0, 1, 10, 100, 1000, 10000, 100000,
        1000000, 10000000, 100000000, 1000000000]

function bigpowten(n,n_bits)
    guess = (n_bits + 1) * 1233 >> 12
    guess += 1
    i = SmallPowersOfTen[guess+1]
    return n < i ? (SmallPowersOfTen[guess], guess-1) : (i,guess)
end

function digitgen(low,w,high,buffer)
    unit::UInt64 = 1
    one = Float(unit << -w.e, w.e)
    too_high = Float(high.s+unit,high.e)
    unsafe_interval = too_high - Float(low.s-unit,low.e)
    integrals = too_high.s >> -one.e
    fractionals = too_high.s & (one.s-1)
    divisor, kappa = bigpowten(integrals, 64 + one.e)
    len = 1
    rest = UInt64(0)
    while kappa > 0
        digit = div(integrals,divisor)
        buffer[len] = 0x30 + digit
        len += 1
        integrals %= divisor
        kappa -= 1
        rest = (UInt64(integrals) << -one.e) + fractionals
        if rest < unsafe_interval.s
            r, kappa = roundweed(buffer, len, rest, UInt64(divisor) << -one.e,
                        unit,kappa,(too_high - w).s,unsafe_interval.s)
            return r, kappa, len
        end
        divisor = div(divisor,10)
    end
    while true
        fractionals *= 10
        unit *= 10
        unsafe_interval = Float(unsafe_interval.s*10,unsafe_interval.e)
        digit = fractionals >> -one.e
        buffer[len] = 0x30 + digit
        len += 1
        fractionals &= one.s - 1
        kappa -= 1
        if fractionals < unsafe_interval.s
            r, kappa = roundweed(buffer,len,fractionals,one.s,
                        unit,kappa,(too_high - w).s*unit,unsafe_interval.s)
            return r, kappa, len
        end
    end
end

function fastshortest(v, buffer = Array{UInt8}(17))
    f = normalize(Float64(v))
    bound_minus, bound_plus = normalizedbound(v)
    ten_mk_min_exp = kMinExp - (f.e + FloatSignificandSize)
    ten_mk_max_exp = kMaxExp - (f.e + FloatSignificandSize)
    cp = binexp_cache(ten_mk_min_exp,ten_mk_max_exp)
    scaled_w = f * cp
    scaled_bound_minus = bound_minus * cp
    scaled_bound_plus = bound_plus * cp
    r, kappa, len = digitgen(scaled_bound_minus,scaled_w,
                             scaled_bound_plus,buffer)
    decimal_exponent = -cp.de + kappa
    return r, len, decimal_exponent+len-1
end


text:
using JuMP, Ipopt
using Plots, LinearAlgebra
default(dpi=100)
using BenchmarkTools
using Random
using StatsBase, StatsModels

clearconsole()

N = 15 #data points
x = collect(range(0,1,length=N)) #training points and data
t = sin.(2*pi.*x) + 0.2 .* randn(size(x))
X = collect(range(-0.5,1.5,length=100)) #test points and data
T = sin.(2*pi.*X)

a = 5e-3
b = 11.1
D = 8 #poly degree
M = D+1 #model param num
W = 0.5.*[1 -1]

s = zeros(size(X))
m = s

for j = 1:N

end

t2 = t
x2 = x
ϕ = broadcast(^, x2',0:D)
Sinv = a.*I + b.*(ϕ*ϕ')
Φ = broadcast(^, X', 0:D)
for k = 1:length(X)
    m[k] = b * Φ[:,k]'*(inv(Sinv)*ϕ*t2)
    s[k] = 1/b + Φ[:,k]'*(inv(Sinv)*Φ[:,k])
end #end for

scatter(x2,t2);
plot!(X,T)
plot!(X,m)

max_stars_repo_name: gbarsih/Safe-Optimal-Rendezvous

max_issues_repo_name: gbarsih/Safe-Optimal-Rendezvous

max_forks_repo_name: gbarsih/Safe-Optimal-Rendezvous


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "QPALM"
version = v"1.0.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/kul-optec/QPALM.git","9bd94ef60ca5d82cd8d5f05ac7efd5a449e501c0"),
    GitSource("https://github.com/kul-optec/LADEL.git", "70ad69a621756e72b164dee302a826f9cc111c55"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd QPALM
rmdir LADEL
mv ../LADEL .
cmake -Bbuild-jl -SQPALM \
    -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
    -DCMAKE_INSTALL_PREFIX:PATH=${prefix} \
    -DQPALM_WITH_JULIA:BOOL=On \
    -DCMAKE_BUILD_TYPE:STRING=Release \
    -DCMAKE_POSITION_INDEPENDENT_CODE:BOOL=On
cmake --build build-jl -j
cmake --install build-jl --component julia_modules
mv LICENSE LICENSE.QPALM
install_license LICENSE.QPALM
mv LADEL/LICENSE LICENSE.LADEL
install_license LICENSE.LADEL
mv LADEL/thirdparty/SuiteSparse/AMD/Doc/License.txt LICENSE.AMD
install_license LICENSE.AMD
"""

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

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

# Dependencies that must be installed before this package can be built
dependencies = Dependency[]

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

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
using Test;

using AdalmPluto;

@testset "libIIO/device.jl" begin
    toggleNoAssertions(true);

    try
        global uri = AdalmPluto.scan("usb", 1, false);
        global ctx = AdalmPluto.createContext(uri);
        global txd,rxd = AdalmPluto.findTRXDevices(ctx);
        global phy = C_iio_context_find_device(ctx, AdalmPluto.PHY_DEVICE_NAME);
    catch
        @error "Could not get a valid uri or context or devices, skipping the remaining tests";
        C_iio_context_destroy(ctx);
        return;
    end

    @testset "Device properties" begin
        #  C_iio_device_get_context (device::Ptr, …)
        @test C_iio_device_get_context(rxd) == ctx;

        saved_data_ptr = C_iio_device_get_data(rxd);
        test_data = Cuchar.(collect("Test data"));
        #  C_iio_device_set_data (device::Ptr, …)
        @test C_iio_device_set_data(rxd, Ptr{Cvoid}(pointer(test_data))) === nothing;
        #  C_iio_device_get_data (device::Ptr, …)
        @test C_iio_device_get_data(rxd) == Ptr{Cvoid}(pointer(test_data));
        C_iio_device_set_data(rxd, saved_data_ptr);

        # NOTE : I'm not actually sure if that ID is constant to the model or the radio or the session
        #  C_iio_device_get_id (device::Ptr, …)
        # @test C_iio_device_get_id(rxd) == "iio:device4";

        #  C_iio_device_get_name (device::Ptr, …)
        @test C_iio_device_get_name(rxd) == AdalmPluto.RX_DEVICE_NAME;

        # Find what that function does, and implement its test
        #  C_iio_device_set_kernel_buffers_count (device::Ptr, …)
        @test_skip C_iio_device_set_kernel_buffers_count(device, UInt32(2));
    end

    # Are there triggers in the AdalmPluto ? Or are those functions meant for other radios ?
    # From libiio doc: Some devices, mostly low speed ADCs and DACs, require a trigger
    @testset "Triggers" begin
        #  C_iio_device_is_trigger (device::Ptr, …)
        @test C_iio_device_is_trigger(rxd) == false;

        for devIdx = 0:C_iio_context_get_devices_count(ctx)-1
            global device = C_iio_context_get_device(ctx, UInt32(devIdx));
            if C_iio_device_is_trigger(device); break; end;
        end

        if !C_iio_device_is_trigger(device)
            @warn "Could not find a trigger device, skipping trigger tests";
        else
            #  C_iio_device_get_trigger (device::Ptr, …)
            @test_skip C_iio_device_get_trigger(rxd) == (-2, Ptr{iio_device}(0));
            #  C_iio_device_set_trigger (device::Ptr, …)
            @test_skip C_iio_device_set_trigger(rxd, device) == 0;
        end

    end

    @testset "Device attributes" begin
        #  C_iio_device_get_attrs_count (device::Ptr, …)
        @test (global nb_attrs = C_iio_device_get_attrs_count(phy)) > 0;

        # uncomment to display all the available attributes
        #  for devIdx = 0:C_iio_context_get_devices_count(ctx)-1
            #  device = C_iio_context_get_device(ctx, UInt32(devIdx));
            #  nb_attrs = C_iio_device_get_attrs_count(device);
            #  for attrIdx = 0:nb_attrs-1
                #  @show attr = C_iio_device_get_attr(device, UInt32(attrIdx));
                #  @show C_iio_device_attr_read(device, attr);
                #  @show C_iio_device_attr_read_bool(device, attr);
                #  @show C_iio_device_attr_read_double(device, attr);
                #  @show C_iio_device_attr_read_longlong(device, attr);
            #  end
        #  end

        if nb_attrs < 1
            @warn "No device attributes, skipping the related tests";
        else
            #  C_iio_device_get_attr (device::Ptr, …)
            @test (global first_attr = C_iio_device_get_attr(phy, UInt32(0))) != "";
            #  C_iio_device_find_attr (device::Ptr, …)
            @test C_iio_device_find_attr(phy, first_attr) == first_attr;

            attr = "calib_mode";
            #  C_iio_device_attr_read (device::Ptr, …)
            @test C_iio_device_attr_read(phy, attr)[1] > 0;
            @test C_iio_device_attr_read(phy, "dummy")[1] < 0;
            saved_attr = C_iio_device_attr_read(phy, attr)[2];
            #  C_iio_device_attr_write (device::Ptr, …)
            @test C_iio_device_attr_write(phy, attr, "manual") > 0 && C_iio_device_attr_read(phy, attr)[2] == "manual";
            @test C_iio_device_attr_write(phy, attr, "auto") > 0 && C_iio_device_attr_read(phy, attr)[2] == "auto";
            #  C_iio_device_attr_write_raw (device::Ptr, …)
            new_calib_mode = b"manual\0";
            @test (C_iio_device_attr_write_raw(phy, attr, Ptr{Cvoid}(pointer(new_calib_mode)), UInt64(sizeof(new_calib_mode))) > 0
                && C_iio_device_attr_read(phy, attr)[2] == "manual");
            # put back the original value
            C_iio_device_attr_write(phy, attr, saved_attr);

            #  C_iio_device_attr_read_bool (device::Ptr, …)
            @test C_iio_device_attr_read_bool(phy, "xo_correction") == (0, true);   # no boolean attributes, so we read a number
            #  C_iio_device_attr_write_bool (device::Ptr, …)
            @test_skip C_iio_device_attr_write_bool(phy, nothing) == 0;             # there are no boolean attributes to write to

            #  C_iio_device_attr_read_double (device::Ptr, …)
            @test C_iio_device_attr_read_double(phy, "xo_correction")[2] > 0;       # same here, real value is an integer
            #  C_iio_device_attr_write_double (device::Ptr, …)
            @test_skip C_iio_device_attr_write_double(phy, nothing) == 0;           # there are no double attributes to write to

            #  C_iio_device_attr_read_longlong (device::Ptr, …)
            @test (global saved_longlong = C_iio_device_attr_read_longlong(phy, "xo_correction")[2]) > 0;
            #  C_iio_device_attr_write_longlong (device::Ptr, …)
            # Could be more robust as we suppose the current value was strictly inferior to the maximum available value
            @test (C_iio_device_attr_write_longlong(phy, "xo_correction", saved_longlong + 1) == 0
                && C_iio_device_attr_read_longlong(phy, "xo_correction") == (0, saved_longlong + 1));
            C_iio_device_attr_write_longlong(phy, "xo_correction", saved_longlong);
        end
    end

    @testset "Buffer attributes" begin
        #  C_iio_device_get_buffer_attrs_count (device::Ptr, …)
        @test (global nb_attrs = C_iio_device_get_buffer_attrs_count(rxd)) > 0;

        # uncomment to display all the available attributes
        #  for devIdx = 0:C_iio_context_get_devices_count(ctx)-1
            #  device = C_iio_context_get_device(ctx, UInt32(devIdx));
            #  nb_attrs = C_iio_device_get_buffer_attrs_count(device);
            #  for attrIdx = 0:nb_attrs-1
                #  @show attr = C_iio_device_get_buffer_attr(device, UInt32(attrIdx));
                #  @show C_iio_device_buffer_attr_read(device, attr);
                #  @show C_iio_device_buffer_attr_read_bool(device, attr);
                #  @show C_iio_device_buffer_attr_read_double(device, attr);
                #  @show C_iio_device_buffer_attr_read_longlong(device, attr);
            #  end
        #  end
        # uncomment to display all the available attributes of a specific device
        #  @show C_iio_device_buffer_attr_read(rxd, "");

        if nb_attrs < 1
            @warn "No device attributes, skipping the related tests";
        else
            #  C_iio_device_get_buffer_attr (device::Ptr, …)
            @test (global first_attr = C_iio_device_get_buffer_attr(rxd, UInt32(0))) != "";
            #  C_iio_device_find_buffer_attr (device::Ptr, …)
            @test C_iio_device_find_buffer_attr(rxd, first_attr) == first_attr;

            attr = "watermark";
            #  C_iio_device_buffer_attr_read (device::Ptr, …)
            @test C_iio_device_buffer_attr_read(rxd, attr)[1] > 0;
            @test C_iio_device_buffer_attr_read(rxd, "dummy")[1] < 0;
            #  saved_attr = C_iio_device_attr_read(rxd, attr)[2];
            #  C_iio_device_buffer_attr_write (device::Ptr, …)
            @test_skip C_iio_device_buffer_attr_write(rxd, attr, nothing);                  # there are no attributes available for user writing
            #  C_iio_device_buffer_attr_write_raw (device::Ptr, …)
            @test_skip C_iio_device_buffer_attr_write_raw(rxd, attr, nothing);              # there are no attributes available for user writing

            #  C_iio_device_buffer_attr_read_bool (device::Ptr, …)
            @test C_iio_device_buffer_attr_read_bool(rxd, attr)[1] == 0;
            #  C_iio_device_buffer_attr_write_bool (device::Ptr, …)
            @test_skip C_iio_device_buffer_attr_write_bool(rxd, attr, nothing);             # there are no boolean attributes to write to

            #  C_iio_device_buffer_attr_read_double (device::Ptr, …)
            @test C_iio_device_buffer_attr_read_double(rxd, attr)[2] > 0;
            #  C_iio_device_buffer_attr_write_double (device::Ptr, …)
            @test_skip C_iio_device_buffer_attr_write_double(rxd, attr, nothing);           # there are no double attributes to write to

            #  C_iio_device_buffer_attr_read_longlong (device::Ptr, …)
            @test C_iio_device_buffer_attr_read_longlong(rxd, attr)[2] > 0;
            #  C_iio_device_buffer_attr_write_longlong (device::Ptr, …)
            @test_skip C_iio_device_buffer_attr_write_longlong(rxd, attr, nothing);         # there are no double attributes available for user writing
        end
    end

    @testset "Channels" begin
        #  C_iio_device_get_channels_count (device::Ptr, …)
        @test (global nb_channels = C_iio_device_get_channels_count(rxd)) > 0;

        if nb_channels < 1
            @warn "No channels found for the tested device, skipping related tests"
        else
            #  C_iio_device_get_channel (device::Ptr, …)
            @test C_iio_device_get_channel(rxd, UInt32(0)) != Ptr{iio_channel}(0);

            #  C_iio_device_find_channel (device::Ptr, …)
            @test C_iio_device_find_channel(phy, "voltage0", false) != Ptr{iio_channel}(0); # RX channel
        end
    end

    # Those functions take a function pointer as an argument, hence their non-implementation
    @testset "Placeholders" begin
        # C_iio_device_attr_read_all
        @test_skip C_iio_device_attr_read_all() == Any;
        # C_iio_device_attr_write_all
        @test_skip C_iio_device_attr_write_all() == Any;
        # C_iio_device_buffer_attr_read_all
        @test_skip C_iio_device_buffer_attr_read_all() == Any;
        # C_iio_device_buffer_attr_write_all
        @test_skip C_iio_device_buffer_attr_write_all() == Any;
    end

    C_iio_context_destroy(ctx);
end


text:
using Test
using MPI

MPI.Init()

comm = MPI.COMM_WORLD
comm_rank = MPI.Comm_rank(comm)
comm_size = MPI.Comm_size(comm)

@test MPI.Comm_compare(comm, comm) == MPI.IDENT
@test comm == comm
@test MPI.Comm_compare(comm, MPI.COMM_SELF) == (MPI.Comm_size(comm) == 1 ? MPI.CONGRUENT : MPI.UNEQUAL)
MPI.Barrier(comm)
comm2 = MPI.Comm_dup(comm)
@test MPI.Comm_compare(comm, comm2) == MPI.CONGRUENT
MPI.Barrier(comm2)
comm3 = MPI.Comm_dup(comm2)
@test MPI.Comm_compare(comm, comm3) == MPI.CONGRUENT
MPI.Barrier(comm3)
MPI.free(comm2)
MPI.Barrier(comm3)
# Don't free comm2
MPI.free(comm3)

splitcomm1 = MPI.Comm_split(comm, comm_rank % 2, 0)
@test div(comm_size, 2) <= MPI.Comm_size(splitcomm1) <= div(comm_size+1, 2)
@test MPI.Comm_rank(splitcomm1) == div(comm_rank, 2)

splitcomm2 = MPI.Comm_split(comm, comm_rank < 2 ? 0 : nothing, 0)
if comm_rank < 2
    @test MPI.Comm_size(splitcomm2) <= 2
    @test MPI.Comm_rank(splitcomm2) == comm_rank
else
    @test splitcomm2 == MPI.COMM_NULL
end

splitcomm3 = MPI.Comm_split_type(comm, comm_rank < 2 ? MPI.COMM_TYPE_SHARED : nothing, 0)
if comm_rank < 2
    @test MPI.Comm_size(splitcomm3) <= 2
    @test MPI.Comm_rank(splitcomm3) < 2
else
    @test splitcomm3 == MPI.COMM_NULL
end

MPI.Finalize()
@test MPI.Finalized()


text:
struct Dataset
    dataset::String
    description::String
    version::String
end

Dataset(dataset::String) = Dataset(dataset,"","")

@memoize function datasets(database::String)

    url = string(BioMart.BIOMART_URL, "?type=datasets&requestid=biomaRt&mart=$(database)")
    r = HTTP.get(url, cookies=true)

    out = Dataset[]

    for line in split(String(r.body), '\n')
        els = split(line, '\t')
        length(els) < 5 && continue

        push!(out, Dataset(
            els[2], els[3], els[5]
        ))
    end
    out
end

"""
    BioMart.datasets(database::Database)
    BioMart.datasets(database::String)

List the datasets available for a given database.

Examples :

    BioMart.datasets("ENSEMBL_MART_ENSEMBL")

    dbs = BioMart.databases()
    BioMart.datasets(dbs[1])

"""
datasets(database::Database) = datasets(database.name)

"""
    BioMart.list_attributes(dataset::Dataset)
    BioMart.list_attributes(dataset::String)

List the attributes available for a given dataset.

Examples :

    BioMart.list_attributes("hsapiens_gene_ensembl")
"""
@memoize function list_attributes(dataset::String)

    url = string(BioMart.BIOMART_URL, "?type=configuration&dataset=$(dataset)")
    r = HTTP.get(url, cookies=true)
    data = String(r.body)

    xdoc = parse_string(data)
    xroot = root(xdoc)

    out = collect_elems(xroot)
    DataFrame(out)
end

list_attributes(dataset::Dataset) = list_attributes(dataset.dataset)

"""go recursively over the XML and collect `AttributeDescription`s"""
function collect_elems(node, out)

    if name(node) == "AttributeDescription"
        el = XMLElement(node)
        push!(out, (
                internalName = attribute(el, "internalName"),
                displayName  = attribute(el, "displayName"),
            )
        )
    else
        for c in child_nodes(node)
            collect_elems(c, out)
        end
    end
    out
end
collect_elems(node) = collect_elems(node, [])


text:
using Plots

function plot_star(star::Star)
    p1 = plot(star.r / Rsun, star.grid / Msun, ylabel="\$(M/M_{sun})\$", xlabel="\$(R/R_{sun})\$", color="blue")
    p2 = plot(star.r / Rsun, star.l / Lsun, xlabel="\$(R/R_{sun})\$", ylabel="\$(L/L_{sun})\$", color="blue")
    p3 = plot(star.r / Rsun, star.P, xlabel="\$(R/R_{sun})\$", ylabel="\$(dyne/cm^2)\$", color="blue")
    p4 = plot(star.r / Rsun, star.T, xlabel="\$(R/R_{sun})\$", ylabel="\$T\$", color="blue")
    plot(p1, p2, p3, p4,layout=(2,2), legend=false)
end


text:
using DependencyParser.DependencyParsing

function build_model(settings::Settings)
  word_ids = Dict{String, Integer}(
    "here"=>1,
    "some"=>2,
    "words"=>3,
    DependencyParser.DependencyParsing.UNKNOWN_TOKEN => 4,
    DependencyParser.DependencyParsing.NULL_TOKEN => 5,
    DependencyParser.DependencyParsing.ROOT_TOKEN => 6
  )
  tag_ids = Dict{String, Integer}(
    "TAG1"=>7,
    "TAG2"=>8,
    "TAG3"=>9,
    DependencyParser.DependencyParsing.UNKNOWN_TOKEN => 10,
    DependencyParser.DependencyParsing.NULL_TOKEN => 11,
    DependencyParser.DependencyParsing.ROOT_TOKEN => 12,
  )
  label_ids = Dict{String, Integer}(
    "LABEL1"=>13,
    "LABEL2"=>14,
    DependencyParser.DependencyParsing.NULL_TOKEN => 15
  )

  embeddings = rand(Float64, 15, settings.embeddings_size)
  hidden_bias = rand(Float64, settings.hidden_size)
  hidden_layer = Dense(settings.embeddings_size * settings.batch_size, settings.hidden_size, bias=hidden_bias)
  output_layer = Dense(settings.hidden_size, 10)

  Model(embeddings, hidden_layer, output_layer, word_ids, tag_ids, label_ids)
end

function build_correct_batch(model::Model, config::Configuration)
  stack = config.stack |> collect

  Integer[
    stack[begin] |> word_id -> model.word_ids[config.sentence.tokens[word_id].name],
    stack[begin + 1] |> word_id -> model.word_ids[config.sentence.tokens[word_id].name],
    model.word_ids[ROOT_TOKEN],
    config.buffer[begin] |> word_id -> model.word_ids[config.sentence.tokens[word_id].name],
    config.buffer[begin + 1] |> word_id -> model.word_ids[config.sentence.tokens[word_id].name],
    config.buffer[begin + 2] |> word_id -> model.word_ids[config.sentence.tokens[word_id].name],
    model.word_ids[NULL_TOKEN],
    model.word_ids[NULL_TOKEN],
    model.word_ids[NULL_TOKEN],
    model.word_ids[NULL_TOKEN],
    model.word_ids[NULL_TOKEN],
    model.word_ids[NULL_TOKEN],
    model.word_ids[config.sentence.tokens[1].name],
    model.word_ids[config.sentence.tokens[3].name],
    model.word_ids[NULL_TOKEN],
    model.word_ids[NULL_TOKEN],
    model.word_ids[NULL_TOKEN],
    model.word_ids[NULL_TOKEN],
    stack[begin] |> word_id -> model.tag_ids[config.sentence.pos_tags[word_id].name],
    stack[begin + 1] |> word_id -> model.tag_ids[config.sentence.pos_tags[word_id].name],
    model.tag_ids[ROOT_TOKEN],
    config.buffer[begin] |> word_id -> model.tag_ids[config.sentence.pos_tags[word_id].name],
    config.buffer[begin + 1] |> word_id -> model.tag_ids[config.sentence.pos_tags[word_id].name],
    config.buffer[begin + 2] |> word_id -> model.tag_ids[config.sentence.pos_tags[word_id].name],
    model.tag_ids[NULL_TOKEN],
    model.tag_ids[NULL_TOKEN],
    model.tag_ids[NULL_TOKEN],
    model.tag_ids[NULL_TOKEN],
    model.tag_ids[NULL_TOKEN],
    model.tag_ids[NULL_TOKEN],
    model.tag_ids[config.sentence.pos_tags[1].name],
    model.tag_ids[config.sentence.pos_tags[3].name],
    model.tag_ids[NULL_TOKEN],
    model.tag_ids[NULL_TOKEN],
    model.tag_ids[NULL_TOKEN],
    model.tag_ids[NULL_TOKEN],
    model.label_ids[NULL_TOKEN],
    model.label_ids[NULL_TOKEN],
    model.label_ids[NULL_TOKEN],
    model.label_ids[NULL_TOKEN],
    model.label_ids[NULL_TOKEN],
    model.label_ids[NULL_TOKEN],
    model.label_ids["SBJ"],
    model.label_ids["IOBJ"],
    model.label_ids[NULL_TOKEN],
    model.label_ids[NULL_TOKEN],
    model.label_ids[NULL_TOKEN],
    model.label_ids[NULL_TOKEN],
  ]
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
module ProcessLibrariesIO

export
dependencies_by_month,
monthly_adjacency_matrices,
get_depversions_and_adj_mats,
make_depversions_and_adj_mats

using Serialization
using JuliaDB

import Dates
import LightGraphs

"""Add Project_Node and Dependency_Project_Node columns to `deps`

Node IDs are natural numbers between 1 and {number of unique IDs}
"""
function _IDs_to_nodes(deps)
    project_ids = vcat(select(deps, :Project_ID), select(deps, :Dependency_Project_ID)) |> unique
    node_numbers = Dict(project_ids[i] => i for i in 1:length(project_ids))

    deps = setcol(deps, :Project_Node, :Project_ID => pid->node_numbers[pid])
    setcol(deps, :Dependency_Project_Node, :Dependency_Project_ID => pid->node_numbers[pid])
end

"""All dependencies of the latest version of each project in each month"""
function dependencies_by_month(versions, dependencies)
    parsed_dates = select(versions, :Published_Timestamp) .|>
        x -> Dates.DateTime(x, Dates.dateformat"y-m-d H:M:S \U\T\C")
    versions = setcol(versions, :Timestamp, parsed_dates)
    versions = setcol(versions, :Month, :Timestamp => ts -> floor(ts, Dates.Month))

    # Latest ID for each project for each month:
    # Sort by timestamp, then number
    latest_version_each_month = groupby(
        (Version_ID=x->sort(x)[end].ID,),
        versions, (:Project_ID, :Month),
        select=(:Timestamp, :Number, :ID))

    # Inner join to get only the deps that match a selected versionID
    # We lose packages here if they don't have any dependents and aren't depended on by anything.
    lonely_packages = setdiff(
        versions.columns.Project_ID,
        dependencies.columns.Project_ID,
        dependencies.columns.Dependency_Project_ID)

    if length(lonely_packages) > 0
        @warn "$(length(lonely_packages)) unconnected packages omitted"
    end

    depversions = join(
        select(latest_version_each_month, (:Version_ID, :Month)),
        dependencies, lkey=:Version_ID, rkey=:Version_ID)

    # Just need to rekey all project_ids to natural numbers and reindex
    depversions |> _IDs_to_nodes |> t -> reindex(t, :Project_Node)
end

"Number of projects in the table"
function number_of_nodes(depversion)
    select(depversion, (:Project_Node, :Dependency_Project_Node)) |>
        maximum ∘ maximum
end

"DiGraph from dependency table"
function graph_from_table(depversion,
                          nnodes = number_of_nodes(depversion))
    g = LightGraphs.SimpleDiGraph(nnodes)
    for row in depversion
        # Create edge from dependency to dependent
        LightGraphs.add_edge!(g, row.Dependency_Project_Node, row.Project_Node)
    end
    return g
end

"""
Make an adjacency matrix from each month.

The nth column contains the dependents of package n. The nth row contains the dependees.
"""
function latest_at_month(deps, month)
    deps = filter(row->row.Month <= month, deps)
    # This allocates memory. Could compute row indexes if we need to.
    groupby(
        (x = grp -> filter(row -> row.Month == maximum(grp.Month), grp),),
        deps, :Project_Node) |> flatten
end

"One project adjacency matrix for each month in deps"
function monthly_adjacency_matrices(deps)
    months = deps.columns.Month
    months = minimum(months):Dates.Month(1):maximum(months)
    nnodes = number_of_nodes(deps)

    # This is slow - do it in parallel

    # optimize grouping in latest_at_month
    if colnames(deps)[1] !== :Project_Node
        deps = reindex(deps, :Project_Node)
    end

    # Drop columns we don't need to avoid allocating excess memory
    deps = select(deps, (:Project_Node, :Month, :Dependency_Project_Node))

    # pre-allocate
    result = Array{Any}(undef, length(months))

    Threads.@threads for i in 1:length(months)
        # Can't use enumerate with @threads
        month = months[i]
        result[i] = latest_at_month(deps, month) |>
            t -> graph_from_table(t, nnodes) |>
            LightGraphs.adjacency_matrix
    end
    result
end

function tensor(adj_mats)
    # Guessing that float float calculations are faster, but haven't checked.
    X = Array{Float64}(undef, (size(adj_mats[1])..., length(adj_mats)))
    [X[:,:,i]=f for (i,f) in adj_mats|>enumerate]
    X
end

# These should obviously be a Make/Drake/DAG thing, but I haven't bothered to make that.

function filter_to_platform(table, platform)
    filter(row -> row.Platform == platform, table)
end

function platpath(platform)
    path = "../data/processed/pacman/$platform/"
    run(`mkdir -p $path`)
    path
end

function make_depversions_and_adj_mats(platform)
    versions = loadtable("../data/sample-1.4/$(platform)_versions-1.4.0-2018-12-22.csv")
    dependencies = loadtable("../data/sample-1.4/$(platform)_dependencies-1.4.0-2018-12-22.csv", type_detect_rows=4000)

    versions = filter_to_platform(versions, platform)
    dependences = filter_to_platform(dependencies, platform)

    depversions = ProcessLibrariesIO.dependencies_by_month(versions, dependencies)
    adj_mats = ProcessLibrariesIO.monthly_adjacency_matrices(depversions);

    ppath = platpath(platform)
    serialize(ppath * "meta.jls", depversions)
    serialize(ppath * "adj-mats.jls", adj_mats)

    depversions, adj_mats
end

function get_meta(platform)
    ppath = platpath(platform)
    if isfile(ppath * "meta.jls")
        meta = deserialize(ppath * "meta.jls")
    else
        meta = make_depversions_and_adj_mats(platform)[1]
    end
    meta
end

function get_adj_mats(platform)
    ppath = platpath(platform)
    if isfile(ppath * "adj-mats.jls")
        adj_mats = deserialize(ppath * "adj-mats.jls")
    else
        adj_mats = make_depversions_and_adj_mats(platform)[2]
    end
    adj_mats
end

"Check ../data/processed and make only if missing"
function get_depversions_and_adj_mats(platform)
    ppath = platpath(platform)
    if isfile(ppath * "meta.jls")
        adj_mats = deserialize(ppath * "adj-mats.jls")
        meta = deserialize(ppath * "meta.jls")
    else
        meta, adj_mats = make_depversions_and_adj_mats(platform)
    end
    meta, adj_mats
end

PLATFORMS = [
   "Elm",
   "Cargo",
   "Pypi",
   "CRAN",
   "Maven",
   "NPM",
  ]

function _make_all()
    # Printing causes segfaults with @threads.
    #= oldlogger = Logging.global_logger() =#
    #= Logging.global_logger(Logging.NullLogger()) =#
    for p in PLATFORMS
        get_depversions_and_adj_mats(p)
        @info "$p done!"
    end
    #= Logging.global_logger(oldlogger) =#
end

end


text:
# utilities.jl

using LinearAlgebra
using CSV
using DataFrames

"""
    setup(path::String)

Setup a first-passage time problem by specifying a path to a .csv file containing the
transition.

The .csv file should have the column names "condition", "from", and "to". Currently, the
transition rates are set to be equal to the number of states in each condition. This keeps
the processing times comparable between conditions, and when fitting the scaling parameter
τ, that τ represents the transition rate per jump.
"""
function setup(path::String)
    df = CSV.read(path, DataFrame, stringtype=String)
    @assert names(df) == ["from", "to"] "Incorrect column names in .csv file."

    statenames = []
    ctr = 1
    nametoidx = Dict()
    rows = []
    columns = []
    # Getting a mapping from state name to its index in the W matrix
    for row = 1:nrow(df)
        if df[row, :from] ∉ statenames
            push!(statenames, df[row, :from])
            nametoidx[df[row, :from]] = ctr
            ctr += 1
        end

        if df[row, :to] ∉ statenames
            push!(statenames, df[row, :to])
            nametoidx[df[row, :to]] = ctr
            ctr += 1
        end
        push!(rows, nametoidx[df[row, :to]])
        push!(columns, nametoidx[df[row, :from]])
    end

    # Setting transition rates
    nstates = length(statenames)
    W = zeros(nstates, nstates)
    for (col, row) in zip(columns, rows)
        W[row, col] = nstates
    end
    W = setdiagonal!(W)

    @assert 0 in W[diagind(W)] "No absorbing states provided."

    # Separating transient and absorbing submatrices
    (transient, absorbing) = getdims(W)
    T = W[transient, transient]
    A = W[absorbing, transient]

    return T, A
end

"""
    setdiagonal!(W::Matrix)

Set the diagonal correctly for W-matrices (van Kampen, 2007, Ch. V.2).
"""
function setdiagonal!(W::Matrix)
    W[diagind(W)] .= 0.0
    W[diagind(W)] = -sum(W, dims=1)
    return W
end

"""
    rescale!(W::Matrix, τ)

Rescale a matrix a matrix by a constant τ.
"""
function rescale!(W::Matrix{Float64}, τ)
    W *= τ
end

"""
    getdims(W::Matrix)

Takes a W-matrix and returns the indices of the transient and absorbing states.
"""
function getdims(W::Matrix)
    transient = []
    absorbing =[]
    for i = 1:size(W, 1)
        if W[i, i] == 0
            push!(absorbing, i)
        else
            push!(transient, i)
        end
    end
    return transient, absorbing
end


text:
using POMDPs, Crux, Flux, POMDPGym, Random, Distributions, POMDPPolicies, PyCall
init_mujoco_render() # Required for visualization
pyimport("safety_gym")
Crux.set_function("isfailure", POMDPGym.isfailure)

# Define the disturbance distribution based on a normal distribution
px = MvNormal([0f0, 0f0], [0.2f0, 0.2f0])

# Construct the MDP
mdp = AdditiveAdversarialPOMDP(EpisodicSafetyGym(pomdp=GymPOMDP(Symbol("Safexp-PointGoal1"), version=:v0), λ=0.01), px)
mdp_nom = AdditiveAdversarialPOMDP(EpisodicSafetyGym(pomdp=GymPOMDP(Symbol("Safexp-PointGoal1"), version=:v0), λ=0.01), px)
# mdp = AdditiveAdversarialPOMDP(GymPOMDP(Symbol("Safexp-PointGoal1"), version=:v0), px)
S = state_space(mdp)

# construct the models for the protagonist
QSA() = ContinuousNetwork(Chain(Dense(62, 256, relu), Dense(256, 256, relu), Dense(256, 1)))
A() = ContinuousNetwork(Chain(Dense(60, 256, relu), Dense(256, 256, relu), Dense(256, 2, tanh)))
Protag() = ActorCritic(A(), DoubleNetwork(QSA(), QSA()))

# And the models for the antagonist
Pf() = ContinuousNetwork(Chain(Dense(62, 256, relu), Dense(256, 256, relu), Dense(256, 1, (x)-> -softplus(-(x-2)))))
function SquashedG()
    base = Chain(Dense(60, 256, relu), Dense(256, 256, relu))
    mu = ContinuousNetwork(Chain(base..., Dense(256, 2)))
    logΣ = ContinuousNetwork(Chain(base..., Dense(256, 2)))
    SquashedGaussianPolicy(mu, logΣ)
end
Antag() = ActorCritic(SquashedG(), Pf())

# Solve with a traditional RL approach
TD3_params() = (N=100000, S=S, buffer_size=Int(1e5), buffer_init=1000, max_steps=1000)

𝒮_td3 = TD3(;π=ActorCritic(A(), DoubleNetwork(QSA(), QSA())), TD3_params()...)
π_td3 = solve(𝒮_td3, mdp)

# Solver with adversarial approach
ARL_params() = (TD3_params()..., desired_AP_ratio=1.0)

# solve with IS
𝒮_isarl = ISARL_Continuous(;π=AdversarialPolicy(Protag(), Antag()), ARL_params()..., px)
π_isarl = solve(𝒮_isarl, mdp, mdp_nom)

# solve with RARL
𝒮_rarl = Crux.RARL_Continuous(;π=AdversarialPolicy(Protag(), Protag(), GaussianNoiseExplorationPolicy(0.2f0)), ARL_params()...)
π_rarl = solve(𝒮_rarl, mdp, mdp_nom)

D_vanilla = steps!(Sampler(mdp, π_td3, S=S, required_columns=[:cost], max_steps=1000), Nsteps=5000)
sum(D_vanilla[:r])/5
sum(D_vanilla[:cost])

D_isarl = steps!(Sampler(mdp, protagonist(π_isarl), S=S, required_columns=[:cost], max_steps=1000), Nsteps=5000)
sum(D_isarl[:r])/5
sum(D_isarl[:cost])

s = tovec(𝒮_isarl.buffer[:s][:, 100], S)
antagonist(π_isarl).A.μ(s)
exp.(antagonist(π_isarl).A.logΣ(s))

using Plots
histogram(𝒮_isarl.buffer[:x][:])


text:
# rotate quaternion
quaternion_rotate(q1::Quaternion,q2::Quaternion) = q2 * q1 / q2

# rotate vector
vector_rotate(v::AbstractVector,q::Quaternion) = Vmat(quaternion_rotate(Quaternion(v), q))
∂vector_rotate∂q(p::AbstractVector, q::Quaternion) = VLmat(q) * Lmat(Quaternion(p)) * Tmat() + VRᵀmat(q) * Rmat(Quaternion(p))

# rotate matrix
function matrix_rotate(A::AbstractMatrix,q::Quaternion)
    c1 = Vmat(quaternion_rotate(Quaternion(A[SVector{3}(1,2,3)]), q))
    c2 = Vmat(quaternion_rotate(Quaternion(A[SVector{3}(4,5,6)]), q))
    c3 = Vmat(quaternion_rotate(Quaternion(A[SVector{3}(7,8,9)]), q))

    return [c1 c2 c3]
end

# rotation matrix
rotation_matrix(q::Quaternion) = VRᵀmat(q) * LVᵀmat(q)
function ∂rotation_matrix∂q(q::Quaternion, p::AbstractVector{T}) where T #; attjac::Bool=false) where T
    M = ∂VRᵀmat∂q(LVᵀmat(q) * p) + VRᵀmat(q) * ∂LVᵀmat∂q(p)
    # if attjac
    #     return (M * LVᵀmat(q))::SMatrix{3,3,T,9}
    # else
    #     return M::SMatrix{3,4,T,12}
    # end
end
function ∂rotation_matrix_inv∂q(q::Quaternion, p::AbstractVector{T}) where T #; attjac::Bool=false) where T
    M = ∂rotation_matrix∂q(inv(q), p) * Tmat()
    # if attjac
    #     return (M * LVᵀmat(q))::SMatrix{3,3,T,9}
    # else
    #     return M::SMatrix{3,4,T,12}
    # end
end


text:
using KahanSummation, Documenter

makedocs(
    modules = [KahanSummation],
    clean = false,
    format = :html,
    sitename = "KahanSummation.jl",
    authors = "Jeff Bezanson, Jeffrey Sarnoff, and other contributors",
    pages = [
        "Home" => "index.md",
    ],
)

deploydocs(
    julia = "1.0",
    repo = "github.com/JuliaMath/KahanSummation.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
module DirectedAcyclicGraphs

include("dags.jl")
include("trees.jl")

end # module