text:
export compute_fem_bdm_div_matrix1, compute_fem_bdm_div_matrix,
        compute_fem_bdm_mass_matrix, compute_fem_bdm_mass_matrix1,
        compute_fem_bdm_skew_matrix, impose_bdm_traction_boundary_condition1,
        impose_bdm_traction_boundary_condition

@doc raw"""
    compute_fem_bdm_div_matrix1(mmesh::Mesh)

Computes the coefficient matrix for

$$\int_\Omega \text{div} \tau \delta u dx$$

Here $\tau \in \mathbb{R}^{2}$ is a second-order tensor (not necessarily symmetric). `mmesh`
uses the BDM1 finite element. The output is a `mmesh.nelem × 2mmesh.nedge` matrix.
"""
function compute_fem_bdm_div_matrix1(mmesh::Mesh)
    @assert mmesh.elem_type == BDM1
    ii = zeros(Int64, get_ngauss(mmesh) * mmesh.elem_ndof)
    jj = zeros(Int64, get_ngauss(mmesh) * mmesh.elem_ndof)
    vv = zeros(get_ngauss(mmesh) * mmesh.elem_ndof)
    @eval ccall((:BDMDiveMatrixMfem, $LIBMFEM), Cvoid, (Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Cdouble}),
        $ii, $jj, $vv)
    sparse(ii, jj, vv, mmesh.nelem, 2mmesh.nedge)
end

@doc raw"""
    compute_fem_bdm_div_matrix(mmesh::Mesh)

Computes the coefficient matrix for

$$\int_\Omega \text{div} \tau \delta u dx$$

Here $\tau \in \mathbb{R}^{2\times 2}$ is a fourth-order tensor (not necessarily symmetric). `mmesh`
uses the BDM1 finite element. The output is a `2mmesh.nelem × 4mmesh.nedge` matrix.
"""
function compute_fem_bdm_div_matrix(mmesh::Mesh)
    C = compute_fem_bdm_div_matrix1(mmesh)
    [C spzeros(mmesh.nelem, 2mmesh.nedge)
    spzeros(mmesh.nelem, 2mmesh.nedge) C]
end

@doc raw"""
    compute_fem_bdm_mass_matrix(alpha::Union{Array{Float64,1}, PyObject},beta::Union{Array{Float64,1}, PyObject}, mmesh::Mesh)

Computes

$$\int_\Omega A\sigma : \delta \tau dx$$

Here

$$A\sigma = \alpha \sigma + \beta \text{tr} \sigma I$$

Here $\sigma$ and $\tau$ are both fourth-order tensors. The output is a  `4mmesh.nedge × 4mmesh.nedge` matrix.
"""
function compute_fem_bdm_mass_matrix(alpha::Union{Array{Float64,1}, PyObject},beta::Union{Array{Float64,1}, PyObject}, mmesh::Mesh)
    @assert mmesh.elem_type == BDM1
    bdm_inner_product_matrix_mfem_ = load_op_and_grad(libmfem,"bdm_inner_product_matrix_mfem", multiple=true)
    alpha,beta = convert_to_tensor(Any[alpha,beta], [Float64,Float64])
    indices, vv = bdm_inner_product_matrix_mfem_(alpha,beta)
    RawSparseTensor(indices, vv, 2mmesh.ndof, 2mmesh.ndof)
end

@doc raw"""
    compute_fem_bdm_mass_matrix(mmesh::Mesh)

Same as [`compute_fem_bdm_mass_matrix`](@ref)
"""
function compute_fem_bdm_mass_matrix(mmesh::Mesh)
    C = compute_fem_bdm_mass_matrix1(mmesh)
    [C spzeros(mmesh.ndof, mmesh.ndof)
    spzeros(mmesh.ndof, mmesh.ndof) C]
end

"""
    compute_fem_bdm_mass_matrix(alpha::Array{Float64,1},beta::Array{Float64,1}, mmesh::Mesh)
"""
function compute_fem_bdm_mass_matrix(alpha::Array{Float64,1},beta::Array{Float64,1}, mmesh::Mesh)
    @assert mmesh.elem_type == BDM1
    @assert length(alpha)==length(beta)==get_ngauss(mmesh)
    N = mmesh.elem_ndof^2 * get_ngauss(mmesh) * 6;
    indices = zeros(Int64, 2N)
    vv = zeros(N)
    @eval ccall((:BDMInnerProductMatrixMfem_forward_Julia, $LIBMFEM),
        Cvoid, (Ptr{Clonglong}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), $indices, $vv, $alpha, $beta)
    # @info minimum(indices[1:2:end].+1), maximum(indices[1:2:end].+1)
    # @info minimum(indices[2:2:end].+1), maximum(indices[2:2:end].+1)
    sparse(indices[1:2:end].+1, indices[2:2:end].+1, vv, 4mmesh.nedge, 4mmesh.nedge)
end

@doc raw"""
    compute_fem_bdm_mass_matrix1(alpha::Array{Float64,1}, mmesh::Mesh)

Computes

$$\int_\Omega \alpha\sigma \cdot \delta \tau dx$$

Here $\alpha$ is a scalar, and $\sigma$ and $\delta \tau$ are second order tensors.

The returned value is a `2mmesh.nedge × 2mmesh.nedge` matrix.
"""
function compute_fem_bdm_mass_matrix1(alpha::Array{Float64,1}, mmesh::Mesh)
    @assert mmesh.elem_type == BDM1
    @assert length(alpha)==get_ngauss(mmesh)
    N = mmesh.elem_ndof^2 * get_ngauss(mmesh);
    indices = zeros(Int64, 2N)
    vv = zeros(N)
    @eval ccall((:BDMInnerProductMatrixMfem1_forward_Julia, $LIBMFEM),
        Cvoid, (Ptr{Clonglong}, Ptr{Cdouble},  Ptr{Cdouble}), $indices, $vv, $alpha)
    sparse(indices[1:2:end].+1, indices[2:2:end].+1, vv, 2mmesh.nedge, 2mmesh.nedge)
end

@doc raw"""
    compute_fem_bdm_mass_matrix1(mmesh::Mesh)

Same as [`compute_fem_bdm_mass_matrix1`](@ref), except that $\alpha\equiv 1$
"""
function compute_fem_bdm_mass_matrix1(mmesh::Mesh)
    compute_fem_bdm_mass_matrix1(ones(get_ngauss(mmesh)), mmesh)
end

@doc raw"""
    compute_fem_bdm_skew_matrix(mmesh::Mesh)

Computes
$$\int_\Omega \sigma : v dx$$
where

$$v = \begin{bmatrix}0 & \rho \\-\rho & 0 \end{bmatrix}$$

Here $\sigma$ is a fourth-order tensor.

The returned value is a `mmesh.nelem × 4mmesh.nedge` matrix.
"""
function compute_fem_bdm_skew_matrix(mmesh::Mesh)
    @assert mmesh.elem_type == BDM1
    N = mmesh.elem_ndof * get_ngauss(mmesh) * 2;
    ii = zeros(Int64, N)
    jj = zeros(Int64, N)
    vv = zeros(N)
    @eval ccall((:BDMSkewSymmetricMatrixMfem, $LIBMFEM),
        Cvoid, (Ptr{Clonglong}, Ptr{Clonglong},  Ptr{Cdouble}), $ii, $jj, $vv)
    sparse(ii, jj, vv, mmesh.nelem, 4mmesh.nedge)
end

@doc raw"""
    impose_bdm_traction_boundary_condition1(gN::Array{Float64, 1}, bdedge::Array{Int64, 2}, mesh::Mesh)

Imposes the BDM traction boundary condition

$$\int_{\Gamma} \sigma \mathbf{n} g_N ds$$

Here $\sigma$ is a second-order tensor. `gN` is defined on the Gauss points, e.g.

```julia
gN = eval_f_on_boundary_edge(func, bdedge, mesh)
```
"""
function impose_bdm_traction_boundary_condition1(gN::Array{Float64, 1}, bdedge::Array{Int64, 2}, mesh::Mesh)
    @assert mesh.elem_type == BDM1
    # sort bdedge so that bdedge[i,1] < bdedge[i,2]
    for i = 1:size(bdedge, 1)
        if bdedge[i,1]>bdedge[i,2]
            bdedge[i,:] = [bdedge[i,2]; bdedge[i,1]]
        end
    end
    order = mesh.lorder
    D = _edge_dict(mesh)
    node_x = zeros(size(bdedge, 1))
    node_y = zeros(size(bdedge, 1))
    ngauss = Int64(@eval ccall((:ComputeFemTractionTermMfem_forward_getNGauss, $LIBMFEM), Cint,
            (Cint,), Int32($order)))
    bdN = size(bdedge, 1)
    @assert length(gN) == ngauss * bdN

    dofs = zeros(Int64, 2bdN)
    for i = 1:bdN
        e = D[(bdedge[i,1], bdedge[i,2])]
        dofs[2*i-1] = e
        dofs[2*i] = e + mesh.nedge
    end

    out = zeros(2bdN)
    bdnode_x, bdnode_y = _traction_get_nodes(bdedge, mesh)
    bdnode_x = bdnode_x'[:]
    bdnode_y = bdnode_y'[:]

    sn = get_boundary_edge_orientation(bdedge, mesh)
    @eval ccall((:ComputeBDMTractionBoundaryMfem_forward_Julia, $LIBMFEM), Cvoid,
        (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint),
        $out, $sn, $gN, $bdnode_x, $bdnode_y, Int32($bdN), Int32($order))
    dofs, out
end

@doc raw"""
    impose_bdm_traction_boundary_condition(g1:Array{Float64, 1}, g2:Array{Float64, 1},
    bdedge::Array{Int64, 2}, mesh::Mesh)

Imposes the BDM traction boundary condition

$$\int_{\Gamma} \sigma \mathbf{n} \cdot \mathbf{g}_N ds$$

Here $\sigma$ is a fourth-order tensor. $\mathbf{g}_N = \begin{bmatrix}g_{N1}\\ g_{N2}\end{bmatrix}$
See [`impose_bdm_traction_boundary_condition1`](@ref).

Returns a `dof` vector and a `val` vector.
"""
function impose_bdm_traction_boundary_condition(g1::Array{Float64, 1}, g2::Array{Float64, 1},
                        bdedge::Array{Int64, 2}, mesh::Mesh)
    d1, v1 = impose_bdm_traction_boundary_condition1(g1, bdedge, mesh)
    d2, v2 = impose_bdm_traction_boundary_condition1(g2, bdedge, mesh)
    [d1;d2 .+ mesh.ndof], [v1;v2]
end


text:
using ReferenceImplementations
using Test

@testset "ReferenceImplementations.jl" begin

    # test long form
    @refimpl function f(x)
        sleep(1)
        return sin(x)
    end
    @test @refimpl(f(1.0)) == sin(1.0)
    f(x) = sin(x)
    @test f(1.0) == @refimpl f(1.0)

    # test short form
    @refimpl g(x) = isnothing(sleep(1)) && return cos(x)
    @test @refimpl(g(1.0)) == cos(1.0)
    g(x) = cos(x)
    @test g(1.0) == @refimpl g(1.0)

end

##
@testset "kwargs" begin
    @refimpl f(x; y=0) = return y
    f(x; y=0) = 0
    @test @refimpl(f(1; y=1)) == 1
    @test @refimpl(f, f(1; y=1)) == 1

    @refimpl kwf(y; x=0.) = cos(x)
    kwf(y; x=0.) = sin(x)
    @test @refimpl kwf(0.; x=0.) == 1.
end

##
@testset "refimpl a specific function" begin
    # fake naive implementation
    @refimpl function f(x)
        return 0
    end
    f(x) = x

    @refimpl function g(x)
        return f(x) * 2
    end
    g(x) = f(x) * 3
    @test @refimpl(f, g(1)) == 0
    @test @refimpl(g, g(1)) == 2
    @test @refimpl(g(1)) == 0
end

##
@testset "nesting" begin

    # fake naive implementation
    @refimpl function f(x)
        print("ref f\n")
        return x
    end

    @refimpl function g(x)
        print("ref g\n")
        return f(x) + 1
    end

    @test 2 == @refimpl g(1)
end

## shouldn't replace methods that don't have reference implementations
@testset "refimpl only specific type signatures" begin
    @refimpl f(x::Int) = x
    f(x::Float64) = x
    @test @refimpl f(0.) == 0.
    @test @refimpl f f(0.) == 0.
end

## no composite arguments (i.e. function definition and other expressions)
@testset "no composites" begin
    @test_throws ArgumentError("To define a new reference implementation, the argument"*
        " of @refimpl must contain only a function definition.") try @eval (@refimpl begin
                    f(x) = 1
                    1 + 1
                end) catch err; throw(err.error) end
end

## escaping local
@testset "escaping locals" begin
    x = 0.0
    @refimpl f(x) = sin(x)
    @test @refimpl f(x) == 0.0
end

## type parameters
@testset "type params" begin
    x = 0.0
    @refimpl f(x::T) where T = sin(x)
    @test @refimpl f(x) == 0.0
end

ext: jl


text:
using ZulipFinBot
using Test

@testset "ZulipFinBot.jl" begin
    # Write your tests here.
end


text:
"""
    mutate(p::StaticParticle{T}, x::T) where {T<:Any}

Mutate a particle - update the location.

- ```p``` Particle to update
- ```x``` New location of particle

"""
function mutate(p::StaticParticle{T}, x::T) where {T<:Any}

    StaticParticle(x, p.ℓw)

end

"""
    mutate(p::DynamicParticle{T}, x::T, τ::Real = p.τ) where {T<:Any}

    Mutate a particle - update the location.

    - ```p``` Particle to update
    - ```x``` New location of particle
    - ```τ``` New temperature (or index) for ℓdens (defaults to old temperature)

"""
function mutate(p::DynamicParticle{T}, x::T, τ::Real = p.τ) where {T<:Any}
    # recalculates cached density
    # logweight = true
    Particle(x, p.ℓw, p.ℓdens, τ, true)

end


text:
import MPI

# min (x1+x2)^2+(x1+x2)*x3 + (x1+x2)*x4
# st.
#     x1 * x2 = 10
#     x2^2 + x3*x1 < 5
#     x2^2 + x4*x1 < 6
# x1, x2 , x3, x4 free variables

include("../../JuliaInterface/ParPipsNlp.jl")

using ParPipsNlp

function str_init_x0(nodeid, x0)
    if(nodeid == 0)
        assert(length(x0) ==2)
        x0[1] = 1.0
        x0[2] = 1.0
    elseif nodeid == 1
        assert(length(x0) ==1)
        x0[1] = 1.0
    elseif nodeid == 2
        assert(length(x0) ==1)
        x0[1] = 1.0
    end
end

function str_prob_info(nodeid,mode,col_lb,col_ub,row_lb,row_ub)

    if(mode == :Values)
        #setting col lb ub, row lb ub
        if(nodeid==0)
            fill!(col_lb,-Inf)
            fill!(col_ub, Inf)
            fill!(row_lb, 10.0)
            fill!(row_ub, 10.0)
        elseif (nodeid == 1 )
            fill!(col_lb, -Inf)
            fill!(col_ub, Inf)
            fill!(row_lb, -Inf)
            fill!(row_ub, 5)
        elseif (nodeid == 2)
            fill!(col_lb, -Inf)
            fill!(col_ub, Inf)
            fill!(row_lb, -Inf)
            fill!(row_ub, 6)
        else
            assert(false)
        end

        @show col_lb
        @show col_ub
        @show row_lb
        @show row_ub
    end

    if(nodeid==0)
        return (2,1)
    elseif(nodeid == 1 || nodeid == 2)
        return (1,1)
    else
        assert(false)
    end
end

function str_eval_f(nodeid,x0,x1)
    @show nodeid
    @show x0
    @show x1
    fval = 0.0
    if nodeid ==  0
        x01 = x0[1]
        x02 = x0[2]
        fval = (x01+x02)^2
    elseif nodeid == 1
        x01 = x0[1]
        x02 = x0[2]
        x13 = x1[1]
        fval = (x01+x02)*x13
    elseif nodeid == 2
        x01 = x0[1]
        x02 = x0[2]
        x24 = x1[1]
        fval = (x01+x02)*x24
    else
        assert(false)
    end
    return fval
end

function str_eval_g(nodeid,x0,x1,new_eq_g, new_inq_g)
    if nodeid == 0
        assert(x0 == x1)
        assert(length(new_inq_g) == 0 )
        assert(length(new_eq_g) == 1)
        x01 = x0[1]
        x02 = x0[2]
        new_eq_g[1] = x01 * x02
    elseif nodeid == 1
        x01 = x0[1]
        x02 = x0[2]
        x13 = x1[1]
        assert(length(new_inq_g) == 1 )
        assert(length(new_eq_g) == 0)
        new_inq_g[1] = x02^2 + x13*x01
    elseif nodeid == 2
        x01 = x0[1]
        x02 = x0[2]
        x24 = x1[1]
        assert(length(new_inq_g) == 1 )
        assert(length(new_eq_g) == 0)
        new_inq_g[1] = x02^2 + x24*x01
    else
        assert(false)
    end
end

function str_eval_grad_f(rowid,colid,x0,x1,new_grad_f)
    @show rowid
    @show colid
    @show x0
    @show x1
    if (0,0) == (rowid, colid)
        assert(length(new_grad_f) == 2)
        assert(x1 == x0)
        x01 = x0[1]
        x02 = x0[2]
        new_grad_f[1] = 2*(x01 + x02)
        new_grad_f[2] = 2*(x01 + x02)
    elseif (1,0) == (rowid, colid)
        assert(length(new_grad_f) == 2)
        x01 = x0[1]
        x02 = x0[2]
        x13 = x1[1]
        new_grad_f[1] = x13
        new_grad_f[2] = x13
    elseif (1,1) == (rowid, colid)
        assert(length(new_grad_f) == 1)
        x01 = x0[1]
        x02 = x0[2]
        x13 = x1[1]
        new_grad_f[1] = x01+x02
    elseif (2,0) == (rowid, colid)
        assert(length(new_grad_f) == 2)
        x01 = x0[1]
        x02 = x0[2]
        x24 = x1[1]
        new_grad_f[1] = x24
        new_grad_f[2] = x24
    elseif (2,2) == (rowid, colid)
        assert(length(new_grad_f) == 1)
        x01 = x0[1]
        x02 = x0[2]
        x24 = x1[1]
        new_grad_f[1] = x01+x02
    else
        assert(false)
    end
end

function str_eval_jac_g(rowid,colid,x0,x1,mode,e_rowidx,e_colptr,e_values,i_rowidx,i_colptr,i_values)
    if(mode == :Structure)
        if (rowid,colid) == (0,0)
            return (2,0)
        elseif (rowid, colid) == (1,1)
            return (0,1)
        elseif (rowid, colid) == (2,2)
            return (0,1)
        elseif (rowid, colid) == (1,0)
            return (0,2)
        elseif (rowid, colid) == (2,0)
            return (0,2)
        else
            assert(false)
        end
    else
        x01 = x0[1]
        x02 = x0[2]
        if (rowid,colid) == (0,0)
            assert(x1==x0)
            assert(length(e_rowidx) == length(e_values) == 2)
            assert(length(i_rowidx) == length(i_values) == 0)
            e_rowidx[1] = 1
            e_rowidx[2] = 1
            e_colptr[1] = 1
            e_colptr[2] = 2
            e_colptr[3] = 3
            e_values[1] = x02
            e_values[2] = x01
        elseif (rowid, colid) == (1,1)
            assert(length(e_rowidx) == length(e_values) == 0)
            assert(length(i_rowidx) == length(i_values) == 1)
            x13 = x1[1]
            i_rowidx[1] = 1
            i_colptr[1] = 1
            i_colptr[2] = 2
            i_values[1] = x01
        elseif (rowid, colid) == (2,2)
            assert(length(e_rowidx) == length(e_values) == 0)
            assert(length(i_rowidx) == length(i_values) == 1)
            x24 = x1[1]
            i_rowidx[1] = 1
            i_colptr[1] = 1
            i_colptr[2] = 2
            i_values[1] = x01
        elseif (rowid, colid) == (1,0)
            assert(length(e_rowidx) == length(e_values) == 0)
            assert(length(i_rowidx) == length(i_values) == 2)
            x13 = x1[1]
            i_rowidx[1] = 1
            i_rowidx[2] = 1
            i_colptr[1] = 1
            i_colptr[2] = 2
            i_colptr[3] = 3
            i_values[1] = x13
            i_values[2] = 2.0*x02
        elseif (rowid, colid) == (2,0)
            assert(length(e_rowidx) == length(e_values) == 0)
            assert(length(i_rowidx) == length(i_values) == 2)
            x24 = x1[1]
            i_rowidx[1] = 1
            i_rowidx[2] = 1
            i_colptr[1] = 1
            i_colptr[2] = 2
            i_colptr[3] = 3
            i_values[1] = x24
            i_values[2] = 2.0*x02
        else
            assert(false)
        end
        convert_to_c_idx(i_rowidx)
        convert_to_c_idx(i_colptr)
        convert_to_c_idx(e_rowidx)
        convert_to_c_idx(e_colptr)
    end
end

function str_eval_h(rowid,colid,x0,x1,obj_factor,lambda,mode,rowidx,colptr,values)
    if(mode == :Structure)
        if(rowid,colid) == (0,0)
            return 3
        elseif (rowid,colid) == (1,1)
            return 0
        elseif (rowid,colid) == (2,2)
            return 0
        elseif (rowid,colid) == (1,0)
            return 3
        elseif (rowid,colid) == (2,0)
            return 3
        elseif (rowid,colid) == (0,1)
            return 2
        elseif (rowid,colid) == (0,2)
            return 2
        else
            assert(false)
        end
    else
        assert(length(lambda) == 1)
        @show (rowid,colid)
        @show colptr
        @show rowidx
        @show values
        if(rowid,colid) == (0,0)
            @show "diag 0,0 "
            rowidx[1] = 1
            rowidx[2] = 2
            rowidx[3] = 2
            colptr[1] = 1
            colptr[2] = 3
            colptr[3] = 4
            values[1] = 2.0 * obj_factor
            values[2] = 2.0 * obj_factor + 1.0 * lambda[1]
            values[3] = 2.0 * obj_factor
        elseif (rowid,colid) == (1,1)
            @show "1,1"
        elseif (rowid,colid) == (2,2)
            @show "2,2"
        elseif (rowid,colid) == (1,0) #the diagnal contribution to 1st stage
            @show "diag  contr - 1, 0 "
            rowidx[1] = 1
            rowidx[2] = 2
            rowidx[3] = 2
            colptr[1] = 1
            colptr[2] = 3
            colptr[3] = 4
            values[1] = 0.0
            values[2] = 0.0
            values[3] = 2.0 * lambda[1]
        elseif (rowid,colid) == (2,0) #the diagnal contribution  to 1st stage
            @show "diag contr - 2, 0 "
            rowidx[1] = 1
            rowidx[2] = 2
            rowidx[3] = 2
            colptr[1] = 1
            colptr[2] = 3
            colptr[3] = 4
            values[1] = 0.0
            values[2] = 0.0
            values[3] = 2.0 * lambda[1]
        elseif (rowid,colid) == (0,1) #liking border
            @show "linking border - 0 ,1 "
            colptr[1] = 1
            colptr[2] = 2
            colptr[3] = 3
            rowidx[1] = 1
            rowidx[2] = 1
            values[1] = 1.0 * obj_factor + 1.0 * lambda[1]
            values[2] = 1.0 * obj_factor
        elseif (rowid,colid) == (0,2) #liking border
            @show "linking border - ", (rowid,colid)
            colptr[1] = 1
            colptr[2] = 2
            colptr[3] = 3
            rowidx[1] = 1
            rowidx[2] = 1
            values[1] = 1.0 * obj_factor + 1.0 * lambda[1]
            values[2] = 1.0 * obj_factor
        else
            assert(false)
        end
        convert_to_c_idx(colptr)
        convert_to_c_idx(rowidx)
    end
end

function convert_to_c_idx(indicies)
    for i in 1:length(indicies)
        indicies[i] = indicies[i] - 1
    end
end

function create()
    comm = MPI.COMM_WORLD
    println("[$(MPI.Comm_rank(comm))/$(MPI.Comm_size(comm))] create problem ")
    @show comm
    model = FakeModel(:Min,0, 2,
        str_init_x0, str_prob_info, str_eval_f, str_eval_g, str_eval_grad_f, str_eval_jac_g, str_eval_h)
    prob = createProblemStruct(comm, model)
    println("end create problem ")
    return prob
end

function solve(prob)
    println("Solve problem")
    @show prob
    ret = solveProblemStruct(prob)
    println("end solve problem")
end

function main()
    MPI.Init()

    prob = create()

    solve(prob)

    MPI.Finalize()
end

main()

max_stars_repo_name: jalving/PIPS

max_issues_repo_name: jalving/PIPS

max_forks_repo_name: jalving/PIPS


text:
function BandstructureGUI(
        uc::UC,
        path = getReciprocalPath(uc, :Gamma, :K, :M, :Gamma),
        slider_range = -2:0.1:2
    ) where {
        LS, D, S <: AbstractSite{LS,D}, UC <: AbstractUnitcell{S}
    }
    labels = unique(uc |> bonds .|> label)
    H = getHoppingHamiltonianDict(uc)
    trigger = Node(nothing)

    fig = Figure()
    sublayout = fig[1, 1] = GridLayout()
    bs_axis = sublayout[1, 1] = Axis(fig)
    uc_scene = sublayout[1, 2] = LScene(fig, scenekw = (raw=false,))
    sublayout[1, 2] = LRect(fig, color=:transparent, strokecolor = :gray24, strokewidth=1)
    colsize!(sublayout, 2, Aspect(1, 1))

    N_rows = div(length(labels), 4)
    sliders = fig[2, 1] = GridLayout()
    for (i, label) in enumerate(labels)
        s = sliders[2 + 2div(i-1, 4), mod1(i, 4)] = LSlider(
            fig, range = slider_range, width = Relative(1)
        )
        sliders[1 + 2div(i-1, 4), mod1(i, 4)] = LText(
            fig,
            map(v -> @sprintf("Coupling %s = %0.2f", label, v), s.value),
            tellwidth=false
        )
        on(v -> (H.couplings[label] = v; trigger[] = trigger[]), s.value)
    end

    path_signal = Node(path)
    bs = Node{AbstractBandstructure}(getBandstructure(uc, path_signal[], H))
    on(_ -> bs[] = getBandstructure(uc, path_signal[], H), trigger)
    plot!(bs_axis, bs, strokewidth=0, markersize=5)
    setup_axis!(bs_axis, bs)

    bz = getBrillouinZone(uc)
    plot!(uc_scene, bz)
    plot!(uc_scene, path_signal, color=:red, overdraw=true)
    setup_axis!(uc_scene, bz)

    fig
end


text:
using Base.Test
using Base.@__FILE__

include("../src/mlst_download_functions.jl")
include("../src/build_db_functions.jl")
include("../src/calling_functions.jl")
include("../src/db_graph.jl")

TEST_DIR = (dirname(@__FILE__))
TMPDIR = mktempdir()
VERSION = "testing"

# encapsulate in a big test to get a summary of all in the end; if a testset fails, it will show the specific testset results,
# otherwise will show just the total # of test passes.
@testset "mentalist_tests" begin

  c_jejuni_pubmlst_dir = string(TMPDIR, "/", "c_jejuni_pubmlst")
  c_jejuni_pubmlst_db_file = string(TMPDIR, "/", "c_jejuni_pubmlst_31.db")
  c_jejuni_pubmlst_copy_dir = string(TMPDIR, "/", "c_jejuni_pubmlst_copy")
  c_jejuni_fastq_file = string("../data/SRR5824107_small.fastq.gz")
  l_pneumophila_cgmlst_dir = string(TMPDIR, "/", "l_pneumophila_cgmlst")
  l_pneumophila_cgmlst_db_file = string(TMPDIR,"/", "l_pneumophila_cgmlst_31.db")
  l_pneumophila_fastq_file = string("../data/SRR801793.fastq.gz")

  c_jejuni_pubmlst_loci_files = ""
  @testset "download_pubmlst_scheme" begin
    mkdir(c_jejuni_pubmlst_dir)
    c_jejuni_pubmlst_loci_files, c_jejuni_profile_file = download_pubmlst_scheme("Campylobacter jejuni", c_jejuni_pubmlst_dir)
    c_jejuni_allele_filenames = open(string(TEST_DIR, "/", "test_data", "/", "c_jejuni_pubmlst_allele_filenames.txt"))
    for c_jejuni_allele_filename in eachline(c_jejuni_allele_filenames)
       @test isfile(string(c_jejuni_pubmlst_dir, "/", chomp(c_jejuni_allele_filename)))
    end
  end

  l_pneumophila_cgmlst_loci_files = ""
  @testset "download_cgmlst_scheme" begin
    mkdir(l_pneumophila_cgmlst_dir)
    l_pneumophila_cgmlst_loci_files = download_cgmlst_scheme("Legionella pneumophila", l_pneumophila_cgmlst_dir)
    l_pneumophila_allele_filenames = open(string(TEST_DIR, "/", "test_data", "/", "l_pneumophila_cgmlst_allele_filenames.txt"))
    for l_pneumophila_allele_filename in eachline(l_pneumophila_allele_filenames)
      @test isfile(string(l_pneumophila_cgmlst_dir, "/", chomp(l_pneumophila_allele_filename)))
    end
  end

  @testset "build_l_pneumophila_db" begin
    K::Int8 = 31
    results, loci = kmer_class_for_each_locus(K, l_pneumophila_cgmlst_loci_files)
    @test typeof(results) == Vector{Tuple{Dict{UInt64,Vector{Int16}},Vector{Int16},Dict{UInt64,Int64}}}
    @test typeof(loci) == Vector{String}

    kmer_classification = combine_loci_classification(K, results, loci)
    loci_list, weight_list, alleles_list, kmer_list, allele_ids_per_locus = kmer_classification
    @test typeof(loci_list) == Vector{Int32}
    @test typeof(weight_list) == Vector{Int16}
    @test typeof(alleles_list) == Vector{Int16}
    @test typeof(kmer_list) == Vector{String}
    @test typeof(allele_ids_per_locus) == Vector{Int64}
    @test length(weight_list) == length(kmer_list)
    @test all([x != 0 for x in weight_list])
    @test all([length(kmer) == K for kmer in kmer_list])

    profile = nothing
    args = Dict("k" => K, "fasta_files" => l_pneumophila_cgmlst_loci_files)
    save_db(K, kmer_classification, loci, l_pneumophila_cgmlst_db_file, profile, args, VERSION)
    @test isfile(l_pneumophila_cgmlst_db_file)
  end

  @testset "build_c_jejuni_db" begin
    K::Int8 = 31
    results, loci = kmer_class_for_each_locus(K, c_jejuni_pubmlst_loci_files)
    @test typeof(results) == Vector{Tuple{Dict{UInt64,Vector{Int16}},Vector{Int16},Dict{UInt64,Int64}}}
    @test typeof(loci) == Vector{String}

    kmer_classification = combine_loci_classification(K, results, loci)
    loci_list, weight_list, alleles_list, kmer_list, allele_ids_per_locus = kmer_classification
    @test typeof(loci_list) == Vector{Int32}
    @test typeof(weight_list) == Vector{Int16}
    @test typeof(alleles_list) == Vector{Int16}
    @test typeof(kmer_list) == Vector{String}
    @test typeof(allele_ids_per_locus) == Vector{Int64}
    @test length(weight_list) == length(kmer_list)
    @test all([x != 0 for x in weight_list])
    @test all([length(kmer) == K for kmer in kmer_list])

    profile = nothing
    args = Dict("k" => K, "fasta_files" => c_jejuni_pubmlst_loci_files)
    save_db(K, kmer_classification, loci, c_jejuni_pubmlst_db_file, profile, args, VERSION)
    @test isfile(c_jejuni_pubmlst_db_file)
  end

  @testset "call_l_pneumophila" begin
    # open the db:
    kmer_db, loci, loci2alleles, k, profile, build_args = open_db(l_pneumophila_cgmlst_db_file)
    @test typeof(kmer_db) == Dict{Bio.Seq.Kmer{Bio.Seq.DNANucleotide,31},Vector{Tuple{Int16,Int16,Vector{Int16}}}}
    @test typeof(loci) == Vector{String}
    @test typeof(loci2alleles) == Dict{Int16,Vector{Int16}}

    # count kmers:
    kmer_count = count_kmers(DNAKmer{k}, [l_pneumophila_fastq_file])
    # test kmer counts:
    @test kmer_count[DNAKmer{k}("AACTAAAACACTGGCTTTACTGCCCCAAATG")] == 19
    @test kmer_count[DNAKmer{k}("TTATGCCTATGAAAATATGATTTTATCTGAA")] == 24
    @test kmer_count[DNAKmer{k}("AAAGGAATTACAATGAATATCTAAATAATTC")] == 14
    @test kmer_count[DNAKmer{k}("AGGAAGTTTTGGTTATATGCTATATGATGCG")] == 21
    @test kmer_count[DNAKmer{k}("AAAAGAGCCAATTTGATGAGCCACCACATCA")] == 20

    # call:
    votes, loci_votes = count_votes(kmer_count, kmer_db, loci2alleles)
    @test typeof(votes) == Dict{Int16,Dict{Int16,Int64}}
    @test typeof(loci_votes) == DataStructures.DefaultDict{Int16,Int64,Int64}
    @test loci_votes[29] == 29212
    @test loci_votes[306] == 10971
    @test loci_votes[1090] == 20743
    @test loci_votes[1316] == 47172
    @test loci_votes[1333] == 16469
    @test votes[29][1] == 3815
    @test votes[29][2] == -640
    @test votes[29][3] == -24757
    @test votes[306][2] == -385
    @test votes[306][18] == -4224
    @test votes[1090][11] == -1917
    @test votes[1090][16] == 439

    # some parameters:
    kmer_thr, max_mutations, output_votes = 1, 5, true
    # call:
    allele_calls, voting_result = call_alleles(k, kmer_count, votes, loci_votes, loci, loci2alleles, l_pneumophila_cgmlst_loci_files, kmer_thr, max_mutations, output_votes)
    @test typeof(allele_calls) == Vector{AlleleCall}
    @test allele_calls[1].allele == "1"
    @test allele_calls[1].depth == 6
    @test allele_calls[2].allele == "1"
    @test allele_calls[2].depth == 9
    @test allele_calls[132].allele == "N" # novel_allele
    @test allele_calls[662].allele == "N" # novel_allele
    @test allele_calls[536].allele == "1"
    @test allele_calls[536].flag == "+" # multiple alleles possible
    @test allele_calls[652].allele == "1"
    @test allele_calls[652].flag == "+" # multiple alleles possible
  end

  @testset "call_c_jejuni" begin
    # open the db:
    kmer_db, loci, loci2alleles, k, profile, build_args = open_db(c_jejuni_pubmlst_db_file)
    @test typeof(kmer_db) == Dict{Bio.Seq.Kmer{Bio.Seq.DNANucleotide,31},Vector{Tuple{Int16,Int16,Vector{Int16}}}}
    @test typeof(loci) == Vector{String}
    @test typeof(loci2alleles) == Dict{Int16,Vector{Int16}}

    # count kmers:
    kmer_count = count_kmers(DNAKmer{k}, [c_jejuni_fastq_file])

    # test kmer counts:
    @test kmer_count[DNAKmer{k}("TCATTTAAGGACTTTTCAGTGATTAAAATCA")] == 2
    @test kmer_count[DNAKmer{k}("CACTCCAATTTTTTCAAATAAAGTAGCTAAG")] == 0
    @test kmer_count[DNAKmer{k}("ATTCTTTTACTCCTATTATCGGTTATACTAA")] == 0
    @test kmer_count[DNAKmer{k}("GAAAAAAGTAATCCAAGGTGCGCAAAAAGCA")] == 4
    @test kmer_count[DNAKmer{k}("AAATATAGTCAATAAATTATAAAAAAAACTT")] == 0

    # call:
    votes, loci_votes = count_votes(kmer_count, kmer_db, loci2alleles)
    @test typeof(votes) == Dict{Int16,Dict{Int16,Int64}}
    @test typeof(loci_votes) == DataStructures.DefaultDict{Int16,Int64,Int64}

    @test loci_votes[1] == 3554
    @test loci_votes[2] == 2065
    @test loci_votes[3] == 1434
    @test loci_votes[4] == 2448
    @test loci_votes[5] == 1683
    @test loci_votes[6] == 3073
    @test loci_votes[7] == 1246

    @test votes[7][288] == 485
    @test votes[7][306] == 440
    @test votes[7][520] == -551
    @test votes[4][306] == 1539
    @test votes[4][29] == 484
    @test votes[4][413] == -769

    # some parameters:
    kmer_thr, max_mutations, output_votes = 2, 5, true
    # call:
    allele_calls, voting_result = call_alleles(k, kmer_count, votes, loci_votes, loci, loci2alleles, c_jejuni_pubmlst_loci_files, kmer_thr, max_mutations, output_votes)
    @test typeof(allele_calls) == Vector{AlleleCall}
    @test [ac.allele for ac in allele_calls] == ["2", "17", "2", "3", "2", "1", "5"]
  end

  @testset "call_c_jejuni_after_moving_db" begin
    mkdir(c_jejuni_pubmlst_copy_dir)
    cp(c_jejuni_pubmlst_db_file, string(c_jejuni_pubmlst_copy_dir, "/", basename(c_jejuni_pubmlst_db_file)))
    cp(c_jejuni_pubmlst_dir, string(c_jejuni_pubmlst_copy_dir, "/", basename(c_jejuni_pubmlst_dir)))
    # open the db:
    kmer_db, loci, loci2alleles, k, profile, build_args = open_db(string(dirname(c_jejuni_pubmlst_copy_dir), "/", basename(c_jejuni_pubmlst_db_file)))
    @test typeof(kmer_db) == Dict{Bio.Seq.Kmer{Bio.Seq.DNANucleotide,31},Vector{Tuple{Int16,Int16,Vector{Int16}}}}
    @test typeof(loci) == Vector{String}
    @test typeof(loci2alleles) == Dict{Int16,Vector{Int16}}

    # count kmers:
    kmer_count = count_kmers(DNAKmer{k}, [c_jejuni_fastq_file])

    # test kmer counts:
    @test kmer_count[DNAKmer{k}("TCATTTAAGGACTTTTCAGTGATTAAAATCA")] == 2
    @test kmer_count[DNAKmer{k}("CACTCCAATTTTTTCAAATAAAGTAGCTAAG")] == 0
    @test kmer_count[DNAKmer{k}("ATTCTTTTACTCCTATTATCGGTTATACTAA")] == 0
    @test kmer_count[DNAKmer{k}("GAAAAAAGTAATCCAAGGTGCGCAAAAAGCA")] == 4
    @test kmer_count[DNAKmer{k}("AAATATAGTCAATAAATTATAAAAAAAACTT")] == 0

    # call:
    votes, loci_votes = count_votes(kmer_count, kmer_db, loci2alleles)
    @test typeof(votes) == Dict{Int16,Dict{Int16,Int64}}
    @test typeof(loci_votes) == DataStructures.DefaultDict{Int16,Int64,Int64}

    @test loci_votes[1] == 3554
    @test loci_votes[2] == 2065
    @test loci_votes[3] == 1434
    @test loci_votes[4] == 2448
    @test loci_votes[5] == 1683
    @test loci_votes[6] == 3073
    @test loci_votes[7] == 1246

    @test votes[7][288] == 485
    @test votes[7][306] == 440
    @test votes[7][520] == -551
    @test votes[4][306] == 1539
    @test votes[4][29] == 484
    @test votes[4][413] == -769

    # some parameters:
    kmer_thr, max_mutations, output_votes = 2, 5, true
    # call:
    allele_calls, voting_result = call_alleles(k, kmer_count, votes, loci_votes, loci, loci2alleles, c_jejuni_pubmlst_loci_files, kmer_thr, max_mutations, output_votes)
    @test typeof(allele_calls) == Vector{AlleleCall}
    @test [ac.allele for ac in allele_calls] == ["2", "17", "2", "3", "2", "1", "5"]
  end

  @testset "complement_alleles" begin
    @test complement_alleles([], 0) == Int16[]
    @test complement_alleles([], 1) == [Int16(1)]
    @test complement_alleles([-1], 0) == Int16[]
    @test complement_alleles([-1], 1) == [Int16(1)]
    @test complement_alleles([0], 0) == Int16[]
    @test complement_alleles([0], 1) == [Int16(1)]
    @test complement_alleles([1], 0) == Int16[]
    @test complement_alleles([1], 1) == Int16[]
    @test complement_alleles([0,1,3,5,7,9], 10) == Int16[2,4,6,8,10]
    @test complement_alleles([0,2,4,6,8,10], 10) == Int16[1,3,5,7,9]
  end

  @testset "twin" begin
    @test twin(DNAKmer{0}, DNAKmer("")) == DNAKmer("")
    @test twin(DNAKmer{1}, DNAKmer("A")) == DNAKmer("T")
    @test twin(DNAKmer{1}, DNAKmer("C")) == DNAKmer("G")
    @test twin(DNAKmer{1}, DNAKmer("G")) == DNAKmer("C")
    @test twin(DNAKmer{1}, DNAKmer("T")) == DNAKmer("A")
    @test twin(DNAKmer{4}, DNAKmer("ACTG")) == DNAKmer("CAGT")
  end

  # rm(TMPDIR, recursive=true)
end

max_stars_repo_name: matnguyen/MentaLiST

max_issues_repo_name: matnguyen/MentaLiST

max_forks_repo_name: matnguyen/MentaLiST


text:
using DifferentialEquations, DataFrames
using Plots
using CSV
include("odeFileExample.jl")
include("variableNames.jl")
f=DDEFunction(ddeFile,Symbol.(syms))
#f=ODEFunction(odeFile,Symbol.(syms))
y0=ones(length(syms))
p=1
h(p,t)=ones(length(syms))
params=1
prob=DDEProblem(f,ones(length(syms)),h,(0.0,10.0))
#prob=ODEProblem(f,ones(length(syms)),(0.0,10.0))
print("Solving Equations \n")
@timev sol=solve(prob)
plot(sol)
png("testOut.png")
df=DataFrame(sol)
CSV.write("out.csv",df)


text:
"""
    TfidfTransformer()

The following is taken largely from scikit-learn's documentation:
https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/feature_extraction/text.py

Convert a collection of raw documents to a matrix of TF-IDF features.

"TF" means term-frequency while "TF-IDF" means term-frequency times
inverse document-frequency.  This is a common term weighting scheme in
information retrieval, that has also found good use in document
classification.

The goal of using TF-IDF instead of the raw frequencies of occurrence
of a token in a given document is to scale down the impact of tokens
that occur very frequently in a given corpus and that are hence
empirically less informative than features that occur in a small
fraction of the training corpus.

The formula that is used to compute the TF-IDF for a term `t` of a
document `d` in a document set is `tf_idf(t, d) = tf(t, d) *
idf(t)`. Assuming `smooth_idf=false`, `idf(t) = log [ n / df(t) ] + 1`
where `n` is the total number of documents in the document set and
`df(t)` is the document frequency of `t`. The document frequency is
the number of documents in the document set that contain the term
`t`. The effect of adding “1” to the idf in the equation above is that
terms with zero idf, i.e., terms that occur in all documents in a
training set, will not be entirely ignored. (Note that the idf formula
above differs from that appearing in standard texts, `idf(t) = log [ n
/ (df(t) + 1) ])`.

If `smooth_idf=true` (the default), the constant “1” is added to the
numerator and denominator of the idf as if an extra document was seen
containing every term in the collection exactly once, which prevents
zero divisions: `idf(t) = log [ (1 + n) / (1 + df(t)) ] + 1`.

The parameters `max_doc_freq` and `min_doc_freq` restrict the vocabulary
that the transformer will consider. `max_doc_freq` indicates that terms in only
up to the specified percentage of documents will be considered. For example, if
`max_doc_freq` is set to 0.9, terms that are in more than 90% of documents
will be removed. Similarly, the `min_doc_freq` parameter restricts terms in the
other direction. A value of 0.01 means that only terms that are at least in 1% of
documents will be included.
"""
mutable struct TfidfTransformer <: AbstractTextTransformer
    max_doc_freq::Float64
    min_doc_freq::Float64
    smooth_idf::Bool
end

function TfidfTransformer(; max_doc_freq::Float64 = 1.0, min_doc_freq::Float64 = 0.0, smooth_idf::Bool = true)
    transformer = TfidfTransformer(max_doc_freq, min_doc_freq, smooth_idf)
    message = MMI.clean!(transformer)
    isempty(message) || @warn message
    return transformer
end

struct TfidfTransformerResult
    vocab::Vector{String}
    idf_vector::Vector{Float64}
end

get_result(::TfidfTransformer, idf::Vector{<:AbstractFloat}, vocab::Vector{String}, ::SparseMatrixCSC) =
    TfidfTransformerResult(vocab, idf)

function build_tfidf!(doc_term_mat::SparseMatrixCSC{T},
                      tfidf::SparseMatrixCSC{F},
                      idf_vector::Vector{F}) where {T <: Real, F <: AbstractFloat}
    rows = rowvals(doc_term_mat)
    dtmvals = nonzeros(doc_term_mat)
    tfidfvals = nonzeros(tfidf)
    @assert size(dtmvals) == size(tfidfvals)

    p, n = size(doc_term_mat)

    # TF tells us what proportion of a document is defined by a term
    words_in_documents = F.(sum(doc_term_mat; dims=1))
    oneval = one(F)

    @inbounds for i = 1:n
        for j in nzrange(doc_term_mat, i)
            row = rows[j]
            tfidfvals[j] = dtmvals[j] / max(words_in_documents[i], oneval) * idf_vector[row]
        end
    end

    return tfidf
end

function _transform(::TfidfTransformer,
                    result::TfidfTransformerResult,
                    v::Corpus)
    doc_terms = build_dtm(v, result.vocab)
    tfidf = similar(doc_terms.dtm, eltype(result.idf_vector))
    build_tfidf!(doc_terms.dtm, tfidf, result.idf_vector)

    # here we return the `adjoint` of our sparse matrix to conform to
    # the `n x p` dimensions throughout MLJ
    return adjoint(tfidf)
end

# for returning user-friendly form of the learned parameters:
function MMI.fitted_params(::TfidfTransformer, fitresult)
    vocab = fitresult.vocab
    idf_vector = fitresult.idf_vector
    return (vocab = vocab, idf_vector = idf_vector)
end

## META DATA

MMI.metadata_pkg(TfidfTransformer,
             name="$PKG",
             uuid="7876af07-990d-54b4-ab0e-23690620f79a",
             url="https://github.com/JuliaAI/MLJText.jl",
             is_pure_julia=true,
             license="MIT",
             is_wrapper=false
)

MMI.metadata_model(TfidfTransformer,
               input_scitype = Union{
                   AbstractVector{<:AbstractVector{STB.Textual}},
                   AbstractVector{<:STB.Multiset{<:ScientificNGram}},
                   AbstractVector{<:STB.Multiset{STB.Textual}}
                   },
               output_scitype = AbstractMatrix{STB.Continuous},
               docstring = "Build TF-IDF matrix from raw documents",
               path = "MLJText.TfidfTransformer"
               )


text:
import Base.GC.gc
import Base.summarysize

mutable struct PPO{T,TP,TV,OP,OV} <: AbstractPolicyAgent where{T<:Real,TP<:AbstractPolicy,TV<:AbstractFuncApprox,OP<:AbstractOptimizer,OV<:AbstractOptimizer}
    π::TP
    vf::TV
    popt::OP
    vopt::OV
    γ::T
    λ::T
    ϵ::T  # clip param
    α::T  # entropy
    steps_per_batch::Int
    num_epochs::Int
    batch_size::Int
    counter::Int
    states::Vector
    actions::Vector
    blogps::Array{T, 1}
    rewards::Array{T, 1}
    values::Array{T, 1}
    terminals::Array{Bool,1}
    gp::Array{T}
    ψ::Array{T}
    gv::Array{T}
    ϕ::Array{T}

    function PPO(π::AbstractPolicy, vf::AbstractFuncApprox, popt::AbstractOptimizer, vopt::AbstractOptimizer, γ::T, λ::T, ϵ::T, α::T, steps_per_batch::Int, num_epochs::Int, batch_size::Int, ::Type{TS}, ::Type{TA}) where {T<:Real,TS,TA}
        π2 = clone(π)
        vf2 = clone(vf)
        popt2 = clone(popt, π2)
        vopt2 = clone(vopt, vf2)
        states = Array{TS, 1}(undef, steps_per_batch)
        actions = Array{TA, 1}(undef, steps_per_batch)
        blogps = Array{T, 1}(undef, steps_per_batch)
        rewards = Array{T, 1}(undef, steps_per_batch)
        values = Array{T, 1}(undef, steps_per_batch+1)
        terminals = Array{Bool, 1}(undef, steps_per_batch)
        gp = zeros(size(get_params(π)))
        ψ = zeros(size(get_params(π)))
        gv = zeros(size(get_params(vf)))
        ϕ = zeros(size(get_params(vf)))
        new{T,typeof(π),typeof(vf),typeof(popt),typeof(vopt)}(π2, vf2, popt2, vopt2, γ, λ, ϵ, α, steps_per_batch, num_epochs, batch_size, 0, states, actions, blogps, rewards, values, terminals, gp, ψ, gv, ϕ)
    end
end

function act!(agent::PPO, env::AbstractEnvironment, rng::AbstractRNG)
    agent.counter += 1
    counter = agent.counter
    # push!(agent.states, env.state)
    agent.states[counter] = deepcopy(env.state)
    logp = get_action!(agent.π, env.state, rng)
    # push!(agent.actions, agent.π.action)
    # push!(agent.blogps, logp)
    agent.actions[counter] = deepcopy(agent.π.action)
    agent.blogps[counter] = logp

    v = call(agent.vf, env.state, 1)
    # push!(agent.values, v)
    agent.values[counter] = v

    reward = step!(env, agent.π.action, rng)
    # push!(agent.rewards, reward)
    agent.rewards[counter] = reward
    terminal = is_terminal(env)
    # push!(agent.terminals, terminal)
    agent.terminals[counter] = terminal
    if counter ≥ agent.steps_per_batch #length(agent.states) ≥ agent.steps_per_batch
        if !terminal
            v′ = call(agent.vf, env.state, 1)
            # push!(agent.values, v′)
            agent.values[counter+1] = v′
        else
            # push!(agent.values, 0.0)
            agent.values[counter+1] = 0.0
        end

        ppo_update!(agent, rng)

        # empty!(agent.states)
        # empty!(agent.actions)
        # empty!(agent.blogps)
        # empty!(agent.rewards)
        # empty!(agent.values)
        # empty!(agent.terminals)
        agent.counter = 0
    end
end

function ppo_make_targets!(LamAdvs, LamRets, γ, λ, rewards, values, terminals)
    prevAdv = 0.0
    T = length(rewards)
    for t in T:-1:1
        r = rewards[t]
        v = values[t]
        v′ =values[t+1]
        terminal = terminals[t]
        δ = r + γ * v′ * (1-terminal) - v
        LamAdvs[t] = δ + γ * λ * (1-terminal) * prevAdv
        LamRets[t] = LamAdvs[t] + v
        prevAdv = LamAdvs[t]
    end
end

function ppo_policy_grad!(g, ψ, π, state, action, blogp, Adv, ϵ, α)
    logp = gradient_logp!(ψ, π, state, action)
    ratio = exp(logp - blogp)
    unclipped = ratio * Adv
    clipped = clamp(ratio, 1.0 - ϵ, 1.0 + ϵ) * Adv
    if (unclipped ≤ clipped)  # if clipped is lower then the ratio has been clipped and gradient would be 0
        @. g += unclipped * ψ
    end
    if α > 0
        gradient_entropy!(ψ, π, state, reuse=true)
        @. g += α * ψ
    end
end

function ppo_value_grad!(g, ϕ, vf, state, G)
    v = call_gradient!(ϕ, vf, state, 1)
    δ = G - v
    @. g += δ * ϕ
end

function ppo_update!(agent::PPO, rng)
    N = length(agent.blogps)
    LamAdvs = zeros(N)
    LamRets = zeros(N)

    ppo_make_targets!(LamAdvs, LamRets, agent.γ, agent.λ, agent.rewards, agent.values, agent.terminals)

    # gp = zeros(size(get_params(agent.π)))
    # ψ = zeros(size(get_params(agent.π)))
    # gv = zeros(size(get_params(agent.vf)))
    # ϕ = zeros(size(get_params(agent.vf)))
    # gp = agent.gp
    # ψ = agent.ψ
    # gv = agent.gv
    # ϕ = agent.ϕ
    for epoch in 1:agent.num_epochs
        idxs = randperm(rng, N)
        i = 1
        count = 0
        while i <= N
            t = idxs[i]
            ppo_policy_grad!(agent.gp, agent.ψ, agent.π, agent.states[t], agent.actions[t], agent.blogps[t], LamAdvs[t], agent.ϵ, agent.α)
            ppo_value_grad!(agent.gv, agent.ϕ, agent.vf, agent.states[t], LamRets[t])
            count += 1
            if count >= agent.batch_size || i == N
                @. agent.gp /= count
                @. agent.gv /= count
                update!(agent.popt, agent.gp)
                update!(agent.vopt, agent.gv)
                count = 0
                fill!(agent.gp, 0.0)
                fill!(agent.gv, 0.0)
            end
            i += 1
        end
    end
    # gc(true)
    # gc(true)
end

function new_episode!(agent::PPO, rng::AbstractRNG)
    nothing
end

function clone(agent::PPO)::PPO
    a = PPO(agent.π, agent.vf, agent.popt, agent.vopt, agent.γ, agent.λ, agent.ϵ, agent.α, agent.steps_per_batch, agent.num_epochs, agent.batch_size, eltype(agent.states), eltype(agent.actions))
    a.states = deepcopy(agent.states)
    a.actions = deepcopy(agent.actions)
    a.blogps = deepcopy(agent.blogps)
    a.rewards = deepcopy(agent.rewards)
    a.values = deepcopy(agent.values)
    a.terminals = deepcopy(agent.terminals)
    a.counter = agent.counter
    a.gp .= agent.gp
    a.gv .= agent.gv
    a.ψ .= agent.ψ
    a.ϕ .= agent.ϕ
    return a
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
using PyPlot
using JLD

D = load("/Users/milan/cams/entropy/gridded_linlog16_nozeroes.jld")

varnames = D["varnames"]
Hlin = D["Hlin"]
Hlog = D["Hlog"]

# vargroups
aervars = [startswith(name,"aer") for name in varnames]
qvars = [name in ["clwc","crwc","cswc","ciwc"] for name in varnames]
hvars = [name in ["h2o2","hcho","hno3","ho2"] for name in varnames]
humid = [name in ["q"] for name in varnames]
nitro = [startswith(name,"n") for name in varnames]
temp = [name in ["t"] for name in varnames]
lnsp = [name in ["lnsp"] for name in varnames]
co2 = [name in ["co2"] for name in varnames]
ozone = [name in ["o3","o3s","go3"] for name in varnames]
velo = [name in ["w","vo","d","z","etadot"] for name in varnames]
alc = [endswith(name,"oh") for name in varnames]
alk = [name in ["c2h4","c2h6","c3h8","c5h8","ch4"] for name in varnames]
sul = [name in ["so2"] for name in varnames]
all_groups = Array{Bool}(aervars+qvars+hvars+humid+nitro+temp+lnsp+co2+ozone+velo+alc+alk+sul)
others = Array{Bool}(.~(all_groups))

# rs = 10.0 .^ Array(-1:0.05:1)
# Hlintheo = [bitentropy(LinQuant16Array(exp.(r*randn(Float32,1000000)))) for r in rs]
# Hlogtheo = [bitentropy(LogQuant16Array(exp.(r*randn(Float32,1000000)))) for r in rs]

## PLOT
pygui(true)
fig,ax = subplots(1,1,figsize=(6,6))

ms = 35
al = 0.7

ax.scatter(Hlin[aervars],Hlog[aervars],ms,alpha=al,"k",label="Aerosols",marker="^")
ax.scatter(Hlin[humid],Hlog[humid],ms,alpha=al,"C4",label="Humidity",marker="o")
ax.scatter(Hlin[qvars],Hlog[qvars],ms,alpha=al,"C0",label="Cloud water species",marker="o")
ax.scatter(Hlin[hvars],Hlog[hvars],ms,alpha=al,"C2",label="Hydrogens",marker="p")
ax.scatter(Hlin[nitro],Hlog[nitro],ms,alpha=al,"C1",label="Nitrogens",marker="v")
ax.scatter(Hlin[temp],Hlog[temp],ms,alpha=al,"C3",label="Temperature",marker="<")
ax.scatter(Hlin[co2],Hlog[co2],ms,alpha=al,"C9",label="Carbon dioxide",marker=">")
ax.scatter(Hlin[lnsp],Hlog[lnsp],ms,alpha=al,"C2",label="Log surface pressure",marker="o")
ax.scatter(Hlin[ozone],Hlog[ozone],ms,alpha=al,"C5",label="Ozone",marker="d")
ax.scatter(Hlin[velo],Hlog[velo],ms,alpha=al,"C6",label="Dynamics",marker="s")
ax.scatter(Hlin[alc],Hlog[alc],ms,alpha=al,"C7",label="Alcohols",marker="<")
ax.scatter(Hlin[alk],Hlog[alk],ms,alpha=al,"C8",label="Alkanes",marker="h")
ax.scatter(Hlin[sul],Hlog[sul],ms,alpha=al,"C4",label="Sulphates",marker=">")
ax.scatter(Hlin[others],Hlog[others],ms,alpha=al,"grey",label="Others",marker="o")

# ax.plot(Hlintheo,Hlogtheo)
#
# for (i,varname) in enumerate(varnames)
#     ax.text(Hlin[i],Hlog[i],varname,fontsize=8)
# end
ax.legend(loc=4,ncol=2,fontsize=9)

n=16

ax.plot([0,n],[0,n],"grey",lw=0.5,zorder=-1)

ax.set_xlim(0,n+1)
ax.set_ylim(0,n+1)
ax.plot([0,n],[n,n],"k")
ax.plot([n,n],[0,n],"k")
ax.text(1,n+0.1,"maximum entropy")
ax.fill_between([0,n,n,n+2],[n,n,0,0],[n+2,n+2,n+2,n+2],alpha=0.3)

ax.set_xlabel("linear quantization entropy [bit]")
ax.set_ylabel("logarithmic quantization entropy [bit]")
ax.set_title("$n-bit quantization entropy",loc="left")

tight_layout()


text:
module Avatars

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export avatar

register_normal_element("q__avatar", context = @__MODULE__)

function avatar(args...;
                wrap::Function = StippleUI.DEFAULT_WRAPPER,
                kwargs...)

  wrap() do
    q__avatar(args...; attributes([kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

end


text:
Makie.plottype(::MultiScaleTreeGraph.Node) = Viz{<:Tuple{MultiScaleTreeGraph.Node}}

"""
    viz(opf::MultiScaleTreeGraph.Node; kwargs...)
    viz!(opf::MultiScaleTreeGraph.Node; kwargs...)

Vizualise the 3D geometry of an MTG (usually read from an OPF). This function search for
the `:geometry` attribute in each node of the MTG, and build the vizualisation using the
`mesh` field, or the reference meshes and the associated transformation matrix if missing.

This function needs 3D information first.

# Examples

```julia
using MultiScaleTreeGraph, PlantGeom, GLMakie

file = joinpath(dirname(dirname(pathof(PlantGeom))),"test","files","simple_plant.opf")
# file = joinpath(dirname(dirname(pathof(PlantGeom))),"test","files","coffee.opf")

opf = read_opf(file)
viz(opf)

# If you need to plot the opf several times, you better cache the mesh in the node geometry
# like so:
transform!(opf, refmesh_to_mesh!)

# Then plot it again like before, and it will be faster:
viz(opf)

# We can also color the 3d plot with several options:
# With one shared color:
viz(opf, color = :red)
# One color per reference mesh:
viz(opf, color = Dict(1 => :burlywood4, 2 => :springgreen4, 3 => :burlywood4))

# Or just changing the color of some:
viz(opf, color = Dict(1 => :burlywood4))
# One color for each vertex of the refmesh 1:
viz(opf, color = Dict(1 => 1:nvertices(get_ref_meshes(opf))[1]))

# Or coloring by opf attribute, e.g. using the mesh max Z coordinates (NB: need to use
# `refmesh_to_mesh!` before, see above):
transform!(opf, :geometry => (x -> zmax(x.mesh)) => :z_max, ignore_nothing = true)
viz(opf, color = :z_max)

# Or even coloring by the value of the Z coordinates of each vertex:
transform!(opf, :geometry => (x -> [i.coords[3] for i in x.mesh.points]) => :z, ignore_nothing = true)
viz(opf, color = :z, showfacets = true)
```

    viz!(ref_meshes; kwargs...)

Plot all reference meshes in a single 3d plot using Makie.

# Examples

```julia
using PlantGeom, GLMakie

file = joinpath(dirname(dirname(pathof(PlantGeom))),"test","files","simple_plant.opf")
opf = read_opf(file)
meshes = get_ref_meshes(opf)

viz(meshes)
# With one shared color:
viz(meshes, color = :green)
# One color per reference mesh:
viz(meshes, color = Dict(1 => :burlywood4, 2 => :springgreen4, 3 => :burlywood4))
# Or just changing the color of some:
viz(meshes, color = Dict(1 => :burlywood4, 3 => :burlywood4))
# One color for each vertex of the refmesh 0:
viz(meshes, color = Dict(2 => 1:nvertices(meshes)[2]))
# Colors as a vector (no missing values allowed here):
viz(meshes, color = [:burlywood4, :springgreen4, :burlywood4])
```
"""
viz, viz!

function Makie.plot!(plot::Viz{<:Tuple{MultiScaleTreeGraph.Node}})
    # function Makie.plot!(plot::Viz{<:Tuple{MultiScaleTreeGraph.Node}}, fig::Makie.GridPosition; axis = NamedTuple(), kwargs...)

    # ax = Makie.Axis(fig[1, 2]; axis...)

    # Mesh list:
    opf = plot[:object][]

    # Plot options:
    color = plot[:color][]
    colormap = get_colormap(plot[:colormap][])
    color_missing = RGBA(0, 0, 0, 0.3)
    ref_meshes = get_ref_meshes(opf)
    colorbar = false

    # Color the meshes:
    if isa(color, Symbol) || typeof(color) <: Colorant
        if color == :slategray3
            # No color is given, we inherit from the default of MeshViz `Viz` (:slategray3) .
            # What we do in this case is use the color given by each reference mesh in the MTG
            # See here for default value in MeshViz:
            # https://github.com/JuliaGeometry/MeshViz.jl/blob/6e37908c78c06212f09229e3e8d92483535ffa16/src/MeshViz.jl#L50
            attr_color = false
            color = get_ref_meshes_color(ref_meshes)
        elseif color in get_attributes(opf)
            # The user provides the name of an attribute from the MTG, coloring using the attribute:
            attr_color = true
            colorbar = true
            # Get the attribute values without nothing values:
            range_val = attribute_range(opf, color)
            # Make a temporary name for our color to use as attribute:
            key_cache = MultiScaleTreeGraph.cache_name(color)

            # Compute the color of each mesh based on the attribute value (use color_missing
            # if no value).
            transform!(
                opf,
                color =>
                    (x -> if x === nothing
                        color_missing
                    else
                        # get the color based on a colormap and the normalized attribute value
                        get_color(colormap, x, range_val)
                    end
                    ) => key_cache
            )
        else
            # User-input single color-value. Here we give the same color to all reference meshes:
            attr_color = false
            color = Dict(zip(keys(ref_meshes.meshes), repeat([color], length(ref_meshes.meshes))))
        end
    elseif length(color) != length(ref_meshes.meshes) && !isa(color, Dict)
        error(
            "color argument should be of type Colorant ",
            "(see [Colors.jl](https://juliagraphics.github.io/Colors.jl/stable/)), or ",
            "an MTG attribute, or a Dict{Int,Colorant} mapping reference meshes to a color."
        )
    else
        attr_color = false
    end

    # If color is not an attribute, it should have the same length as number of RefMeshes, or
    # if not, we provide the reference mesh color:
    if attr_color == false && length(color) != length(ref_meshes.meshes)
        new_color = Dict{Int,Any}(color)
        ref_cols = get_ref_meshes_color(ref_meshes)
        missing_mesh_input = setdiff(collect(keys(ref_cols)), collect(keys(color)))
        for i in missing_mesh_input
            push!(new_color, i => ref_cols[i])
        end
        color = new_color
    end

    facetcolor = plot[:facetcolor][]
    showfacets = plot[:showfacets][]
    colormap = plot[:colormap][]

    if attr_color == false
        # Make the plot, case where the color is a color for each reference mesh:
        traverse!(
            opf,
            node -> viz!(
                # ax,
                plot,
                node[:geometry].mesh === nothing ? refmesh_to_mesh(node) : node[:geometry].mesh,
                color = color[get_ref_mesh_index!(node, ref_meshes)],
                facetcolor = facetcolor,
                showfacets = showfacets,
                colormap = colormap
            )
            ;
            # scale = scale,
            # symbol = symbol,
            # link = link,
            filter_fun = node -> node[:geometry] !== nothing
        )
        #? NB: implement scale / symbol / link / filter_fun filtering to be able to plot only
        #? a subset of the plant/scene. This will be especially usefull when we have different
        #? kind of geometries at different scales of representation.
    else
        # Make the plot, case where the color is an attribute from the MTG:
        traverse!(
            opf,
            function (node)
                viz!(
                    # ax,
                    plot,
                    node[:geometry].mesh === nothing ? refmesh_to_mesh(node) : node[:geometry].mesh,
                    color = node[key_cache],
                    facetcolor = facetcolor,
                    showfacets = showfacets,
                    colormap = colormap
                )
                pop!(node, key_cache) # Remove the cached variable
            end;
            # scale = scale,
            # symbol = symbol,
            # link = link,
            filter_fun = node -> node[:geometry] !== nothing
        )
    end

    # if colorbar
    # cb = Makie.Colorbar(plot.parent, label = string(color), colormap = colormap, colorrange = range_val)
    # end

    # return Makie.AxisPlot(ax, plot)
end


text:
module IO

export outMCMC

using DelimitedFiles

macro name(arg)
    x = string(arg)
    quote
        $x
    end
end

function outMCMC(folder::String,thisVar,output)
        out0 = open(pwd()*"/$(thisVar)Out", "a")
        writedlm(out0, output)
        close(out0)
end

end


text:
using NonNegLeastSquaresMLJInterface: NonNegativeLeastSquareRegressor
using MLJBase
using Random
using Test

@testset "Testing NonNegativeLeastSquareRegressor MLJ Interface" begin
    Random.seed!(123)
    n = 100
    X = randn(n, 3)
    y = X * [2, -10, 7] + randn(n)
    model = NonNegativeLeastSquareRegressor(;alg=:nnls)
    mach = machine(model, X, y)

    @testset " With fit_intercept" begin
        # Testing fit and fitted_params method
        fit!(mach)
        fp = fitted_params(mach)

        @test all(fp.coefs .>= 0)
        @test fp.intercept >= 0

        # testing predict method
        ypred = predict(mach)

        mse = sum((ypred - y).^2) / n
        avg_mse = sum((y .- mean(y)).^2) / n

        @test mse < avg_mse
    end

    @testset " Without fit_intercept" begin
        model.fit_intercept = false
        # Testing fit and fitted_params method
        fit!(mach)
        fp = fitted_params(mach)

        @test all(fp.coefs .>= 0)
        @test fp.intercept === nothing

        # testing predict method
        ypred = predict(mach)

        mse = sum((ypred - y).^2) / n
        avg_mse = sum((y .- mean(y)).^2) / n

        @test mse < avg_mse
    end

end


text:
# This file was generated, do not modify it.

using MLJ, RDatasets, Random

data = dataset("datasets", "USArrests")
names(data)

describe(data, :mean, :std)

X = select(data, Not(:State))
X = coerce(X, :UrbanPop=>Continuous);

@load PCA pkg=MultivariateStats

pca_mdl = PCA(pratio=1)
pca = machine(pca_mdl, X)
fit!(pca)

W = transform(pca, X);

schema(W).names

r = report(pca)
cumsum(r.principalvars ./ r.tvar)

data = dataset("ISLR", "OJ")

X = select(data, [:PriceCH, :PriceMM, :DiscCH, :DiscMM, :SalePriceMM,
                  :SalePriceCH, :PriceDiff, :PctDiscMM, :PctDiscCH]);

Random.seed!(1515)

@pipeline SPCA(std = Standardizer(),
               pca = PCA(pratio=1-1e-4))
spca_mdl = SPCA()
spca = machine(spca_mdl, X)
fit!(spca)
W = transform(spca, X)
names(W)

r  = report(spca).reports[1]
cs = cumsum(r.principalvars ./ r.tvar)

using PyPlot

figure(figsize=(8,6))

bar(1:length(cs), cs)
plot(1:length(cs), cs, color="red", marker="o")

xlabel("Number of PCA features", fontsize=14)
ylabel("Ratio of explained variance", fontsize=14)

savefig("assets/literate/ISL-lab-10-g1.svg") # hide

Random.seed!(1515)

@load KMeans
@pipeline SPCA2(std = Standardizer(),
                pca = PCA(),
                km = KMeans(k=3))

spca2_mdl = SPCA2()
spca2 = machine(spca2_mdl, X)
fit!(spca2)

assignments = report(spca2).reports[1].assignments
mask1 = assignments .== 1
mask2 = assignments .== 2
mask3 = assignments .== 3;

using PyPlot

figure(figsize=(8, 6))
for (m, c) in zip((mask1, mask2, mask3), ("red", "green", "blue"))
    plot(W[m, 1], W[m, 2], ls="none", marker=".", markersize=10, color=c)
end

xlabel("PCA-1", fontsize=13)
ylabel("PCA-2", fontsize=13)
legend(["Group 1", "Group 2", "Group 3"], fontsize=13)

savefig("assets/literate/ISL-lab-10-cluster.svg") # hide


text:
target_quantity(m) =
    target_scitype_union(m) <: Tuple ? true : false
is_probabilistic(::Type{<:Model}) = false
is_probabilistic(::Type{<:Deterministic}) = false
is_probabilistic(::Type{<:Probabilistic}) = true

# TODO: depreciate? currently used by ensembles.jl
# output_is(modeltype::Type{<:Model}) =
#     [is_probabilistic(modeltype), output_kind(modeltype), target_quantity(modeltype)]
# output_is(model::Model) = output_is(typeof(model))

function coretype(M)
    if isdefined(M, :name)
        return M.name
    else
        return coretype(M.body)
    end
end

name(M::Type{<:Model}) = split(string(coretype(M)), '.')[end] |> String

if VERSION < v"1.0.0"
    import Base.info
end

function info(M::Type{<:Supervised})

    message = "$M has a bad trait declaration. "

    load_path != "unknown" || error(message*"MLJBase.load_path($M) should be defined so that "*
                                    "using MLJ; import MLJ.load_path($M) loads $M into current namespace.")

    # check target_scitype_union:
    T = target_scitype_union(M)
    T <: Union{Found,NTuple{<:Found}} ||
        error(message*"target_scitype_union($M) (defining upper bound of target scitype) "*
              "is not a subtype of Found. ")

    input_scitype_union(M) <: Union{Missing, Found} ||
        error(message*"input_scitype_union($M) (defining upper bound of input scitypes) not a subtype of Union{Missing,Found}. ")

    input_is_multivariate(M) in [false, true] ||
        error(message*"input_is_multivariate($M) must return true or false. ")

    is_pure_julia(M) in [true, false, :unknown] ||
        error(message*"is_pure_julia($M) must return true or false. ")

    d = Dict{Symbol,Any}()
    d[:is_supervised] = true
    d[:load_path] = load_path(M)
    d[:name] = name(M)
    d[:target_scitype_union] = target_scitype_union(M)
    d[:input_scitype_union] = input_scitype_union(M)
    d[:input_is_multivariate] = input_is_multivariate(M)
    d[:is_pure_julia] = is_pure_julia(M)
    d[:package_name] = package_name(M)
    d[:package_uuid] = package_uuid(M)
    d[:package_url] = package_url(M)
    d[:is_probabilistic] = is_probabilistic(M)
    d[:is_wrapper] = is_wrapper(M)
    return d
end

function info(M::Type{<:Unsupervised})

    message = "$M has a bad trait declaration. "

    load_path != "unknown" || error(message*"MLJBase.load_path($M) should be defined so that "*
                                    "using MLJ; import MLJ.load_path($M) loads $M into current namespace.")

    output_scitype_union(M) <: Union{Missing,Found} ||
        error(message*"output_scitype_union($M) (defining upper bound of output scitypes) is not a subtype of Union{Missing,Found}. ")

    output_is_multivariate(M) in [false, true] ||
        error(message*"output_is_multivariate($M) must return true or false.")

    input_scitype_union(M) <: Union{Missing, Found} ||
        error(message*"input_scitype_union($M) (defining upper bound of input scitypes) not a subtype of Union{Missing,Found}. ")

    input_is_multivariate(M) in [false, true] ||
        error(message*"input_is_multivariate($M) must return true or false.")

    is_pure_julia(M) in [true, false, :unknown] ||
        error(message*"is_pure_julia($M) must return true or false. ")

    d = Dict{Symbol,Any}()
    d[:is_supervised] = false
    d[:load_path] = load_path(M)
    d[:name] = name(M)
    d[:output_scitype_union] = output_scitype_union(M)
    d[:output_is_multivariate] = output_is_multivariate(M)
    d[:input_scitype_union] = input_scitype_union(M)
    d[:input_is_multivariate] = input_is_multivariate(M)
    d[:is_pure_julia] = is_pure_julia(M)
    d[:package_name] = package_name(M)
    d[:package_uuid] = package_uuid(M)
    d[:package_url] = package_url(M)
    d[:is_wrapper] = is_wrapper(M)

    return d
end

info(model::Model) = info(typeof(model))


text:
# This file was generated, do not modify it. # hide
ŷ_tuned = predict(mach_restored, Xtest);
ŷ_tuned[1]

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
module Broadcast

using ..Cartesian
import Base.promote_eltype
import Base.@get!
import Base.num_bit_chunks, Base.@_msk_end, Base.unsafe_bitgetindex
import Base.(.+), Base.(.-), Base.(.*), Base.(./), Base.(.\)
import Base.(.==), Base.(.<), Base.(.!=), Base.(.<=)
export broadcast, broadcast!, broadcast_function, broadcast!_function, bitbroadcast
export broadcast_getindex, broadcast_setindex!

## Broadcasting utilities ##

droparg1(a, args...) = args

longer_tuple(x::(), retx::Tuple, y::(), rety::Tuple) = retx
longer_tuple(x::(), retx::Tuple, y::Tuple, rety::Tuple) = rety
longer_tuple(x::Tuple, retx::Tuple, y::(), rety::Tuple) = retx
longer_tuple(x::Tuple, retx::Tuple, y::Tuple, rety::Tuple) =
    longer_tuple(droparg1(x...), retx, droparg1(y...), rety)
longer_tuple(x::Tuple, y::Tuple) = longer_tuple(x, x, y, y)

longer_size(x::Union(AbstractArray,Number)) = size(x)
longer_size(x::Union(AbstractArray,Number), y::Union(AbstractArray,Number)...) =
    longer_tuple(size(x), longer_size(y...))

# Calculate the broadcast shape of the arguments, or error if incompatible
broadcast_shape() = ()
function broadcast_shape(As::Union(AbstractArray,Number)...)
    sz = longer_size(As...)
    nd = length(sz)
    bshape = ones(Int, nd)
    for A in As
        for d = 1:ndims(A)
            n = size(A, d)
            if n != 1
                if bshape[d] == 1
                    bshape[d] = n
                elseif bshape[d] != n
                    error("arrays could not be broadcast to a common size")
                end
            end
        end
    end
    return tuple(bshape...)::typeof(sz)
end

# Check that all arguments are broadcast compatible with shape
function check_broadcast_shape(shape::Dims, As::Union(AbstractArray,Number)...)
    for A in As
        if ndims(A) > length(shape)
            error("cannot broadcast array to have fewer dimensions")
        end
        for k in 1:ndims(A)
            n, nA = shape[k], size(A, k)
            if n != nA != 1
                error("array could not be broadcast to match destination")
            end
        end
    end
end

## Broadcasting core
# Generate the body for a broadcasting function f_broadcast!(B, A1, A2, ..., A$narrays),
# using function f, output B, and inputs As...
# B must have already been set to the appropriate size.

# version using cartesian indexing
function gen_broadcast_body_cartesian(nd::Int, narrays::Int, f::Function)
    F = Expr(:quote, f)
    quote
        @assert ndims(B) == $nd
        @ncall $narrays check_broadcast_shape size(B) k->A_k
        @nloops($nd, i, B,
            d->(@nexprs $narrays k->(j_d_k = size(A_k, d) == 1 ? 1 : i_d)), # pre
            begin # body
                @nexprs $narrays k->(@inbounds v_k = @nref $nd A_k d->j_d_k)
                @inbounds (@nref $nd B i) = (@ncall $narrays $F v)
            end)
    end
end

# version using start/next for iterating over the arguments
function gen_broadcast_body_iter(nd::Int, narrays::Int, f::Function)
    F = Expr(:quote, f)
    quote
        @assert ndims(B) == $nd
        @ncall $narrays check_broadcast_shape size(B) k->A_k
        @nexprs 1 d->(@nexprs $narrays k->(state_k_0 = state_k_{$nd} = start(A_k)))
        @nexprs $nd d->(@nexprs $narrays k->(skip_k_d = size(A_k, d) == 1))
        @nloops($nd, i, B,
            d->(@nexprs $narrays k->(state_k_{d-1} = state_k_d)),           # pre
            d->(@nexprs $narrays k->(skip_k_d || (state_k_d = state_k_0))), # post
            begin # body
                @nexprs $narrays k->(@inbounds (v_k, state_k_0) = next(A_k, state_k_0))
                @inbounds (@nref $nd B i) = (@ncall $narrays $F v)
            end)
    end
end

## Broadcasting cores specialized for returning a BitArray

const bitcache_chunks = 64 # this can be changed
const bitcache_size = 64 * bitcache_chunks # do not change this

function dumpbitcache(Bc::Vector{Uint64}, bind::Int, C::Vector{Bool})
    ind = 1
    nc = min(bitcache_chunks, length(Bc)-bind+1)
    for i = 1:nc
        u = uint64(1)
        c = uint64(0)
        for j = 1:64
            C[ind] && (c |= u)
            ind += 1
            u <<= 1
        end
        Bc[bind] = c
        bind += 1
    end
end

# using cartesian indexing
function gen_broadcast_body_cartesian_tobitarray(nd::Int, narrays::Int, f::Function)
    F = Expr(:quote, f)
    quote
        @assert ndims(B) == $nd
        @ncall $narrays check_broadcast_shape size(B) k->A_k
        C = Array(Bool, bitcache_size)
        Bc = B.chunks
        ind = 1
        cind = 1
        @nloops($nd, i, B,
            d->(@nexprs $narrays k->(j_d_k = size(A_k, d) == 1 ? 1 : i_d)), # pre
            begin # body
                @nexprs $narrays k->(@inbounds v_k = @nref $nd A_k d->j_d_k)
                @inbounds C[ind] = (@ncall $narrays $F v)
                ind += 1
                if ind > bitcache_size
                    dumpbitcache(Bc, cind, C)
                    cind += bitcache_chunks
                    ind = 1
                end
            end)
        if ind > 1
            @inbounds C[ind:bitcache_size] = false
            dumpbitcache(Bc, cind, C)
        end
    end
end

# using start/next
function gen_broadcast_body_iter_tobitarray(nd::Int, narrays::Int, f::Function)
    F = Expr(:quote, f)
    quote
        @assert ndims(B) == $nd
        @ncall $narrays check_broadcast_shape size(B) k->A_k
        C = Array(Bool, bitcache_size)
        Bc = B.chunks
        ind = 1
        cind = 1
        @nexprs 1 d->(@nexprs $narrays k->(state_k_0 = state_k_{$nd} = start(A_k)))
        @nexprs $nd d->(@nexprs $narrays k->(skip_k_d = size(A_k, d) == 1))
        @nloops($nd, i, B,
            d->(@nexprs $narrays k->(state_k_{d-1} = state_k_d)),           # pre
            d->(@nexprs $narrays k->(skip_k_d || (state_k_d = state_k_0))), # post
            begin # body
                @nexprs $narrays k->(@inbounds (v_k, state_k_0) = next(A_k, state_k_0))
                @inbounds C[ind] = (@ncall $narrays $F v)
                ind += 1
                if ind > bitcache_size
                    dumpbitcache(Bc, cind, C)
                    cind += bitcache_chunks
                    ind = 1
                end
            end)
        if ind > 1
            @inbounds C[ind:bitcache_size] = false
            dumpbitcache(Bc, cind, C)
        end
    end
end

function gen_broadcast_function(genbody::Function, nd::Int, narrays::Int, f::Function)
    As = [symbol("A_"*string(i)) for i = 1:narrays]
    body = genbody(nd, narrays, f)
    @eval begin
        local _F_
        function _F_(B, $(As...))
            $body
        end
        _F_
    end
end

function gen_broadcast_function_tobitarray(genbody::Function, nd::Int, narrays::Int, f::Function)
    As = [symbol("A_"*string(i)) for i = 1:narrays]
    body = genbody(nd, narrays, f)
    @eval begin
        local _F_
        function _F_(B::BitArray, $(As...))
            $body
        end
        _F_
    end
end

for (Bsig, Asig, gbf, gbb) in
    ((BitArray                          , Union(Array,BitArray,Number)                   ,
      :gen_broadcast_function_tobitarray, :gen_broadcast_body_iter_tobitarray     ),
     (Any                               , Union(Array,BitArray,Number)                   ,
      :gen_broadcast_function           , :gen_broadcast_body_iter                ),
     (BitArray                          , Any                                     ,
      :gen_broadcast_function_tobitarray, :gen_broadcast_body_cartesian_tobitarray),
     (Any                               , Any                                     ,
      :gen_broadcast_function           , :gen_broadcast_body_cartesian           ))

    @eval let cache = Dict{Function,Dict{Int,Dict{Int,Function}}}()
        global broadcast!
        function broadcast!(f::Function, B::$Bsig, As::$Asig...)
            nd = ndims(B)
            narrays = length(As)

            cache_f    = @get! cache      f       Dict{Int,Dict{Int,Function}}()
            cache_f_na = @get! cache_f    narrays Dict{Int,Function}()
            func       = @get! cache_f_na nd      $gbf($gbb, nd, narrays, f)

            func(B, As...)
            B
        end
    end  # let broadcast_cache
end

broadcast(f::Function, As...) = broadcast!(f, Array(promote_eltype(As...), broadcast_shape(As...)), As...)

bitbroadcast(f::Function, As...) = broadcast!(f, BitArray(broadcast_shape(As...)), As...)

broadcast!_function(f::Function) = (B, As...) -> broadcast!(f, B, As...)
broadcast_function(f::Function) = (As...) -> broadcast(f, As...)

broadcast_getindex(src::AbstractArray, I::AbstractArray...) = broadcast_getindex!(Array(eltype(src), broadcast_shape(I...)), src, I...)
@ngenerate N typeof(dest) function broadcast_getindex!(dest::AbstractArray, src::AbstractArray, I::NTuple{N, AbstractArray}...)
    check_broadcast_shape(size(dest), I...)  # unnecessary if this function is never called directly
    checkbounds(src, I...)
    @nloops N i dest d->(@nexprs N k->(j_d_k = size(I_k, d) == 1 ? 1 : i_d)) begin
        @nexprs N k->(@inbounds J_k = @nref N I_k d->j_d_k)
        @inbounds (@nref N dest i) = (@nref N src J)
    end
    dest
end

@ngenerate N typeof(A) function broadcast_setindex!(A::AbstractArray, x, I::NTuple{N, AbstractArray}...)
    checkbounds(A, I...)
    shape = broadcast_shape(I...)
    @nextract N shape d->(length(shape) < d ? 1 : shape[d])
    if !isa(x, AbstractArray)
        @nloops N i d->(1:shape_d) d->(@nexprs N k->(j_d_k = size(I_k, d) == 1 ? 1 : i_d)) begin
            @nexprs N k->(@inbounds J_k = @nref N I_k d->j_d_k)
            @inbounds (@nref N A J) = x
        end
    else
        X = x
        # To call setindex_shape_check, we need to create fake 1-d indexes of the proper size
        @nexprs N d->(fakeI_d = 1:shape_d)
        Base.setindex_shape_check(X, (@ntuple N fakeI)...)
        k = 1
        @nloops N i d->(1:shape_d) d->(@nexprs N k->(j_d_k = size(I_k, d) == 1 ? 1 : i_d)) begin
            @nexprs N k->(@inbounds J_k = @nref N I_k d->j_d_k)
            @inbounds (@nref N A J) = X[k]
            k += 1
        end
    end
    A
end

## elementwise operators ##

.*(As::AbstractArray...) = broadcast(*, As...)
.%(A::AbstractArray, B::AbstractArray) = broadcast(%, A, B)

eltype_plus(As::AbstractArray...) = promote_eltype(As...)
eltype_plus(As::AbstractArray{Bool}...) = typeof(true+true)

.+(As::AbstractArray...) = broadcast!(+, Array(eltype_plus(As...), broadcast_shape(As...)), As...)

type_minus(T, S) = promote_type(T, S)
type_minus(::Type{Bool}, ::Type{Bool}) = typeof(true-true)

function .-(A::AbstractArray, B::AbstractArray)
    broadcast!(-, Array(type_minus(eltype(A), eltype(B)), broadcast_shape(A,B)), A, B)
end

type_div(T,S) = promote_type(T,S)
type_div{T<:Integer,S<:Integer}(::Type{T},::Type{S}) = typeof(one(T)/one(S))
type_div{T,S}(::Type{Complex{T}},::Type{Complex{S}}) = Complex{type_div(T,S)}
type_div{T,S}(::Type{Complex{T}},::Type{S})          = Complex{type_div(T,S)}
type_div{T,S}(::Type{T},::Type{Complex{S}})          = Complex{type_div(T,S)}

function ./(A::AbstractArray, B::AbstractArray)
    broadcast!(/, Array(type_div(eltype(A), eltype(B)), broadcast_shape(A, B)), A, B)
end

function .\(A::AbstractArray, B::AbstractArray)
    broadcast!(\, Array(type_div(eltype(A), eltype(B)), broadcast_shape(A, B)), A, B)
end

type_pow(T,S) = promote_type(T,S)
type_pow{S<:Integer}(::Type{Bool},::Type{S}) = Bool
type_pow{S}(T,::Type{Rational{S}}) = type_pow(T, type_div(S, S))

function .^(A::AbstractArray, B::AbstractArray)
    broadcast!(^, Array(type_pow(eltype(A), eltype(B)), broadcast_shape(A, B)), A, B)
end

## element-wise comparison operators returning BitArray ##

for (f, scalarf, bitf, bitfbody) in ((:.==, :(==), :biteq , :(~a $ b)),
                                     (:.< , :<   , :bitlt , :(~a & b)),
                                     (:.!=, :!=  , :bitneq, :( a $ b)),
                                     (:.<=, :<=  , :bitle , :(~a | b)))
    @eval begin
        ($f)(A::AbstractArray, B::AbstractArray) = bitbroadcast($scalarf, A, B)
        ($bitf)(a::Uint64, b::Uint64) = $bitfbody
        function ($f)(A::AbstractArray{Bool}, B::AbstractArray{Bool})
            local shape
            try
                shape = promote_shape(size(A), size(B))
            catch
                return bitbroadcast($scalarf, A, B)
            end
            F = BitArray(shape)
            Fc = F.chunks
            Ac = bitpack(A).chunks
            Bc = bitpack(B).chunks
            if !isempty(Ac) && !isempty(Bc)
                for i = 1:length(Fc) - 1
                    Fc[i] = ($bitf)(Ac[i], Bc[i])
                end
                msk = @_msk_end length(F)
                Fc[end] = msk & ($bitf)(Ac[end], Bc[end])
            end
            return F
        end
    end
end

# note: the following are not broadcasting, but need to be defined here to avoid
# ambiguity warnings

for (f, cachef, scalarf) in ((:.==, :bitcache_eq , :(==)),
                             (:.< , :bitcache_lt , :<   ),
                             (:.!=, :bitcache_neq, :!=  ),
                             (:.<=, :bitcache_le , :<=  ))
    for (sigA, sigB, expA, expB, shape) in ((:Any, :AbstractArray,
                                             :A, :(B[ind]),
                                             :(size(B))),
                                            (:AbstractArray, :Any,
                                             :(A[ind]), :B,
                                             :(size(A))))
        @eval begin
            ($cachef)(A::AbstractArray, B::AbstractArray, l::Int, ind::Int, C::Vector{Bool}) = 0
            function ($cachef)(A::$sigA, B::$sigB, l::Int, ind::Int, C::Vector{Bool})
                left = l - ind + 1
                @inbounds begin
                    for j = 1:min(bitcache_size, left)
                        C[j] = ($scalarf)($expA, $expB)
                        ind += 1
                    end
                    C[left+1:bitcache_size] = false
                end
                return ind
            end
            function ($f)(A::$sigA, B::$sigB)
                F = BitArray($shape)
                l = length(F)
                l == 0 && return F
                Fc = F.chunks
                C = Array(Bool, bitcache_size)
                ind = 1
                cind = 1
                for i = 1:div(l + bitcache_size - 1, bitcache_size)
                    ind = ($cachef)(A, B, l, ind, C)
                    dumpbitcache(Fc, cind, C)
                    cind += bitcache_chunks
                end
                return F
            end
        end
    end
end

## specialized element-wise operators for BitArray

(.^)(A::BitArray, B::AbstractArray{Bool}) = (B .<= A)
(.^)(A::AbstractArray{Bool}, B::AbstractArray{Bool}) = (B .<= A)

function bitcache_pow{T}(Ac::Vector{Uint64}, B::Array{T}, l::Int, ind::Int, C::Vector{Bool})
    left = l - ind + 1
    @inbounds begin
        for j = 1:min(bitcache_size, left)
            C[j] = unsafe_bitgetindex(Ac, ind) ^ B[ind]
            ind += 1
        end
        C[left+1:bitcache_size] = false
    end
    return ind
end
function (.^){T<:Integer}(A::BitArray, B::Array{T})
    local shape
    try
        shape = promote_shape(size(A), size(B))
    catch
        return bitbroadcast(^, A, B)
    end
    F = BitArray(shape)
    l = length(F)
    l == 0 && return F
    Ac = A.chunks
    Fc = F.chunks
    C = Array(Bool, bitcache_size)
    ind = 1
    cind = 1
    for i = 1:div(l + bitcache_size - 1, bitcache_size)
        ind = bitcache_pow(Ac, B, l, ind, C)
        dumpbitcache(Fc, cind, C)
        cind += bitcache_chunks
    end
    return F
end

for (sigA, sigB) in ((BitArray, BitArray),
                     (AbstractArray{Bool}, BitArray),
                     (BitArray, AbstractArray{Bool}))
    @eval function (.*)(A::$sigA, B::$sigB)
        try
            return bitpack(A) & bitpack(B)
        catch
            return bitbroadcast(&, A, B)
        end
    end
end

end # module

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
module Denumpyfy

# Write your package code here.

replace_exponents(str) = replace(str, "**" => "^")

# range(0, self.n_grid):
function replace_range(str)
    isempty(str) && return str
    !occursin(r"range\(", str) && return str
    caps = match(r"range\((.*)\):$", str).captures[1]
    #caps = filter(!isspace, caps)

    # Handle single argument case by checking for existence of a comma
    # TODO Assume there is a `len` call?
    !occursin(r",", caps) && return string("1:", caps)
    for m in eachmatch(r"(range\((.+?)\):)", str, overlap = false)
        #cap = (m.match)[2:end-1] # get the stuff inside the [...]
        splits = split(m.captures[2], ',')
        args = handle_indexing_arg.(splits, false)
        #idx = "range(" * join(idx, ',') * ')'
        args = join(args, ':')
        str = replace(str, m.captures[1] => args)
    end
    str
end

    # Try and be smart about 1-based indexing
    # first, second = split(caps, ',')

    # # Handle either one being a positive digit and make it 1-based index
    # if all(isdigit, first)
    #     first = parse(Int, first) + 1
    # end
    # if all(isdigit, second)
    #     second = parse(Int, second) + 1
    # end
    # if second == "-1"
    #     second = "end"
    # end
    # if second == "-2"
    #     second = "end-1"
    # end
    # if second == "-3"
    #     second = "end-2"
    # end
    # # If you need more may the lord have mercy on your code
    # return string(first,':', second)
#end

replace_eye(str) = replace(str, "eye(" => "I(")

# INPUT:
# A possible indexing like
# 1, 2 => 1,2
# -1, ngrid - 1 => end,ngrid-1
# 0, 1 => 1,2
function handle_indexing_arg(str, ignore_range = true)
        str = filter(!isspace, str)
        # Special case "0"
        if ignore_range
            str == "0" && return "begin"
        else
            str == "0" && return "1"
        end
        # If is all digits, bump
        all(isdigit, str) && return string(parse(Int, str) + 1)
        # if uses "-1" or some other negative integer

        # TODO add range(1, n_grid - 1): cases
        if ignore_range
            str == "-1" && return "end"
            str == "-2" && return "end-1"
            str == "-3" && return "end-2"
            str == "-4" && return "end-3"
            str == "-5" && return "end-4"
            str == "-6" && return "end-5"
            str == "-7" && return "end-6"
        else
            # lop off the -1 at the end
            endswith(str, "-1") && return str[begin:end-2]
        end
        str
end

replace_len(str) = replace(str, "len(" => "length(")

# linspace(half_delta - x_out, x_out - half_delta, n_grid)
replace_linspace(str) = replace(str, "linspace(" => "LinRange(")

# def construct(self, tol, it_max):
replace_def(str) = replace(str, r"def (.*):$" => s"function \1")

# if condition:
replace_if(str) = replace(str, r"if (.*):$" => s"if \1")
# while cond:
replace_while(str) = replace(str, r"while (.*):$" => s"while \1")
# else:
replace_else(str) = replace(str, r"else:$" => "else")

# zeros((n_grid, n_grid, n_grid))
replace_zeros(str) = replace(str, r"zeros\(\((.*)\)\)" => s"zeros(\1)")
# empty(()) == what in Julia???
# ones((n_grid, n_grid))
replace_ones(str) = replace(str, r"ones\(\((.*)\)\)" => s"ones(\1)")

# TODO handle slices

# [\d] increase?
function replace_indexing(str)
    isempty(str) && return str
    !occursin(r"\[|\]", str) && return str

    for m in eachmatch(r"\[(.+?)\]", str, overlap = false)
        cap = (m.match)[2:end-1] # get the stuff inside the [...]
        splits = split(cap, ',')
        idx = handle_indexing_arg.(splits)
        idx = '[' * join(idx, ',') * ']'
        str = replace(str, m.match => idx)
    end
    return str
end

# More/other letters possible but maybe are problematic with collisions :/
function replace_greek(str)
    for (greek, letter) in zip("αβργΔϵθκστω", ["alpha", "beta", "rho", "gamma", "delta", "epsilon", "theta", "kappa", "sigma", "tau", "omega"])
        str = replace(str, letter => greek)
    end
    str
end

#l_x px_x

const replacements! = [
    replace_exponents,
    replace_range,
    replace_len,
    replace_linspace,
    replace_def,
    replace_if,
    replace_else,
    replace_while,
    replace_zeros,
    replace_ones,
    replace_greek,
    replace_eye,
    replace_indexing
]

function denumpyfy(file)
    jlfilename = replace(file, r"^(\w+).py$" => s"\1.jl")
    contents = readlines(file, keep = false) .|> String
    for f in replacements!
        contents .= f.(contents)
    end
    open(jlfilename, "w") do x
        for s in contents
            println(x, s)
        end
    end
end

export replace_exponents
export replace_range
export replace_len
export replace_linspace
export replace_def
export replace_if
export replace_else
export replace_while
export replace_zeros
export replace_ones
export replace_greek
export handle_indexing_arg
export replace_eye
export replace_indexing
export denumpyfy

end


text:
using ParameterHandling: Positive, Bounded

@testset "parameters_scalar.jl" begin
    @testset "postive" begin
        @testset "$val" for val in [5.0, 0.001f0, 1.0e-7]
            p = positive(val)
            test_parameter_interface(p)
            @test value(p) ≈ val
            @test typeof(value(p)) === typeof(val)
        end

        # Test edge cases around the size of the value relative to the error tol.
        @test_throws ArgumentError positive(-0.1)
        @test_throws ArgumentError positive(1e-11)
        @test value(positive(1e-11, exp, 1e-12)) ≈ 1e-11
    end

    @testset "bounded" begin
        @testset "$val" for val in [-0.05, -0.1 + 1e-12, 2.0 - 1e-11, 2.0 - 1e-12]
            p = bounded(val, -0.1, 2.0)
            test_parameter_interface(p)
            @test value(p) ≈ val
        end

        @test_throws ArgumentError bounded(-0.05, 0.0, 1.0)
    end
end


text:
dt = .01
t_end = 1e3 # ms
in_current = 2.

ns_of_neurons = 10 .^ [1,2,3,4]

V_REST = -65. # * b2.mV
V_RESET = -65. # * b2.mV
FIRING_THRESHOLD = -50. # * b2.mV
MEMBRANE_RESISTANCE = 10. # * b2.Mohm
MEMBRANE_TIME_SCALE = 10. # * b2.ms
ABSOLUTE_REFRACTORY_PERIOD = 0. # * b2.ms

###############################################################################
## EventBasedIntegrator
###############################################################################

using Pkg; Pkg.activate("./../../core/"); Pkg.instantiate()
push!(LOAD_PATH, string(pwd(),"./../../core/src/lifintegrator/core/"))
using EventBasedIntegrator, LinearAlgebra, ProgressMeter
using TimerOutputs
using CPUTime, BenchmarkTools, HDF5

include("./eventbasedlifintegrator_defs/eventbasedlifintegrator_defs.jl")

function balanced_julia(n_of_neurons; doEuler = false)
    network_julia, spike_monitor_julia, profiling_time_julia =
        simulate_balanced_network(10 * in_current; # should actually be voltage (R = 1 Ω)
                            n_of_neurons = n_of_neurons,
                            simulation_time = t_end,
                            doEuler = doEuler, dt = dt,
                            v_rest=V_REST,
                            v_reset=V_RESET,
                            firing_threshold=FIRING_THRESHOLD,
                            membrane_resistance=MEMBRANE_RESISTANCE,
                            membrane_time_scale=MEMBRANE_TIME_SCALE,
                            abs_refractory_period=ABSOLUTE_REFRACTORY_PERIOD)

    spike_times_julia = [spike[2] for spike in spike_monitor_julia.spikes]
    return profiling_time_julia, spike_times_julia
end

# event based

fid = h5open("./data/balancednet_julia.h5", "w")
write(fid, "ns_of_neurons", ns_of_neurons)
profiling_times_julia = zeros(length(ns_of_neurons))
for i in 1:length(ns_of_neurons)
    t, spts = balanced_julia(ns_of_neurons[i]; doEuler = false)
    write(fid, string("spike_train_",ns_of_neurons[i]), spts)
    profiling_times_julia[i] = t
end
write(fid, "profiling_times", profiling_times_julia)
close(fid)

# euler

fid = h5open("./data/balancednet_julia_euler.h5", "w")
write(fid, "ns_of_neurons", ns_of_neurons)
profiling_times_julia_euler = zeros(length(ns_of_neurons))
for i in 1:length(ns_of_neurons)
    t, spts = balanced_julia(ns_of_neurons[i]; doEuler = true)
    write(fid, string("spike_train_",ns_of_neurons[i]), spts)
    profiling_times_julia_euler[i] = t
end
write(fid, "profiling_times", profiling_times_julia_euler)
close(fid)

# get the (2) weights for init of python network?

###############################################################################
## brian2 using C++ (cpp_standalone) and exact integration
###############################################################################

using PyCall

py"""
import sys
import brian2 as b2

sys.path.append('/Users/Bernd/Documents/PhD/Drafts/draft_mnistbenchmarks/BioPlausibleShallowDeepLearning/scripts/LIF_integrator_benchmarks/brian2_defs/')
from brian2_defs import simulate_balanced_network
"""

function balanced_brian2(n_of_neurons; first_init = false)
    py"""
    import sys
    import brian2 as b2

    if $first_init:
        b2.set_device('cpp_standalone', build_on_run=False)
    else:
        #brian2.__init__.clear_cache()
        b2.device.reinit()
        b2.device.activate()
        b2.set_device('cpp_standalone', build_on_run=False)
    spike_monitor_brian2, profiling_summary = simulate_balanced_network($in_current * b2.namp,
                            n_of_neurons = $n_of_neurons,
                            simulation_time = $t_end * b2.ms,
                            dt = $dt,
                            v_rest=$V_REST * b2.mV,
                            v_reset=$V_RESET * b2.mV,
                            firing_threshold=$FIRING_THRESHOLD * b2.mV,
                            membrane_resistance=$MEMBRANE_RESISTANCE * b2.Mohm,
                            membrane_time_scale=$MEMBRANE_TIME_SCALE * b2.ms,
                            abs_refractory_period=$ABSOLUTE_REFRACTORY_PERIOD * b2.ms)
    """
    spike_monitor_brian2, profiling_summary = py"spike_monitor_brian2", py"profiling_summary"
    profiling_times = []
    for x in profiling_summary push!(profiling_times, x[2][1]) end
    profiling_time_brian2 = sum(profiling_times)

    spike_monitor_brian2 = py"spike_monitor_brian2"
    spike_times_brian2 = 1e3 .* spike_monitor_brian2.spike_trains()[0]
    return profiling_time_brian2, spike_times_brian2
end

fid = h5open("./data/balancednet_brian2.h5", "w")
write(fid, "ns_of_neurons", ns_of_neurons)
profiling_times_brian2 = zeros(length(ns_of_neurons))
for i in 1:length(ns_of_neurons)
    t, spts = balanced_brian2(ns_of_neurons[i]; first_init = (i == 1))
    write(fid, string("spike_train_",ns_of_neurons[i]), spts)
    profiling_times_brian2[i] = t
end
write(fid, "profiling_times", profiling_times_brian2)
close(fid)

###############################################################################

#print(string("Time elapsed brian2 (profiling)[ms]: ", 1e3 .* profiling_time_brian2,"\n"))
#print(string("Time elapsed julia (median_of_profiling) [ms]: ", 1e3 .* profiling_time_julia,"\n"))