text:
function SimpleBounds(arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return SimpleBounds((Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function get_lower(obj::SimpleBounds)
    return jcall(obj, "getLower", Vector{jdouble}, ())
end

function get_upper(obj::SimpleBounds)
    return jcall(obj, "getUpper", Vector{jdouble}, ())
end

function unbounded(::Type{SimpleBounds}, arg0::jint)
    return jcall(SimpleBounds, "unbounded", SimpleBounds, (jint,), arg0)
end


text:
module TensorOperations

using TupleTools
using Strided
using Strided: AbstractStridedView, UnsafeStridedView
using LinearAlgebra
using LinearAlgebra: mul!, BLAS.BlasFloat

# export macro API
export @tensor, @tensoropt, @optimalcontractiontree

export enable_blas, disable_blas, enable_cache, disable_cache, clear_cache, cachesize

# export function based API
export tensorcopy, tensoradd, tensortrace, tensorcontract, tensorproduct, scalar
export tensorcopy!, tensoradd!, tensortrace!, tensorcontract!, tensorproduct!

# Convenient type alias
const IndexTuple{N} = NTuple{N,Int}

# An exception type for reporting errors in the index specificatino
struct IndexError{S<:AbstractString} <: Exception
    msg::S
end

# Index notation
#----------------
@nospecialize
include("indexnotation/tensormacro.jl")
include("indexnotation/tensorexpressions.jl")
include("indexnotation/ncontree.jl")
@specialize
include("indexnotation/optimaltree.jl")
include("indexnotation/poly.jl")

# Implementations
#-----------------
include("implementation/indices.jl")
include("implementation/lrucache.jl")
include("implementation/tensorcache.jl")
include("implementation/stridedarray.jl")
include("implementation/diagonal.jl")

# Functions
#----------
include("functions/simple.jl")
include("functions/inplace.jl")

# Global package settings
#------------------------
# A switch for enabling/disabling the use of BLAS for tensor contractions
const _use_blas = Ref(true)
use_blas() = _use_blas[]
function disable_blas()
    _use_blas[] = false
    return
end
function enable_blas()
    _use_blas[] = true
    return
end

# A cache for temporaries of tensor contractions
const cache = LRU{Symbol,Any}()
const _use_cache = Ref(true)
use_cache() = _use_cache[]

"""
    disable_cache()

Disable the cache for further use but does not clear its current contents.
Also see [`clear_cache()`](@ref)
"""
function disable_cache()
    _use_cache[] = false
    return
end

"""
    enable_cache(; maxsize::Int = ..., maxrelsize::Real = 0.5)

(Re)-enable the cache for further use; set the maximal size `maxsize` (as number of bytes)
or relative size `maxrelsize`, as a fraction between 0 and 1, resulting in
`maxsize = floor(Int, maxrelsize * Sys.total_memory())`.
"""
function enable_cache(; kwargs...)
    _use_cache[] = true
    setsize!(cache; kwargs...)
    return
end

"""
    clear_cache()

Clear the current contents of the cache.
"""
function clear_cache()
    empty!(cache)
    return
end

"""
    cachesize()

Return the current memory size (in bytes) of all the objects in the cache.
"""
cachesize() = cache.currentsize

# Some precompile statements
#----------------------------
function _precompile_()
    AVector = Vector{Any}
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(Tuple{typeof(_findfirst), Base.Fix2{typeof(Base.isequal), Symbol}, AVector})
    @assert precompile(Tuple{typeof(_findnext), Base.Fix2{typeof(Base.isequal), Int64}, AVector, Int64})
    @assert precompile(Tuple{typeof(_intersect), Base.BitArray{1}, Base.BitArray{1}})
    @assert precompile(Tuple{typeof(_intersect), Base.BitSet, Base.BitSet})
    @assert precompile(Tuple{typeof(_intersect), UInt128, UInt128})
    @assert precompile(Tuple{typeof(_intersect), UInt32, UInt32})
    @assert precompile(Tuple{typeof(_intersect), UInt64, UInt64})
    @assert precompile(Tuple{typeof(_isemptyset), Base.BitArray{1}})
    @assert precompile(Tuple{typeof(_isemptyset), Base.BitSet})
    @assert precompile(Tuple{typeof(_isemptyset), UInt128})
    @assert precompile(Tuple{typeof(_isemptyset), UInt32})
    @assert precompile(Tuple{typeof(_isemptyset), UInt64})
    @assert precompile(Tuple{typeof(_ncontree!), AVector, Vector{Vector{Int64}}})
    @assert precompile(Tuple{typeof(_setdiff), Base.BitArray{1}, Base.BitArray{1}})
    @assert precompile(Tuple{typeof(_setdiff), Base.BitSet, Base.BitSet})
    @assert precompile(Tuple{typeof(_setdiff), UInt128, UInt128})
    @assert precompile(Tuple{typeof(_setdiff), UInt32, UInt32})
    @assert precompile(Tuple{typeof(_setdiff), UInt64, UInt64})
    @assert precompile(Tuple{typeof(_union), Base.BitArray{1}, Base.BitArray{1}})
    @assert precompile(Tuple{typeof(_union), Base.BitSet, Base.BitSet})
    @assert precompile(Tuple{typeof(_union), UInt128, UInt128})
    @assert precompile(Tuple{typeof(_union), UInt32, UInt32})
    @assert precompile(Tuple{typeof(_union), UInt64, UInt64})
    @assert precompile(Tuple{typeof(addcost), Power{:χ, Int64}, Power{:χ, Int64}})
    @assert precompile(Tuple{typeof(degree), Power{:x, Int64}})
    @assert precompile(Tuple{typeof(deindexify_contraction), Int, Int, Expr, Int, AVector, AVector, Int})
    @assert precompile(Tuple{typeof(deindexify_generaltensor), Int, Int, Expr, Int, AVector, AVector, Int})
    @assert precompile(Tuple{typeof(deindexify_linearcombination), Int, Int, Expr, Int, AVector, AVector, Int})
    @assert precompile(Tuple{typeof(deindexify), Int, Int, Expr, Int, AVector, AVector, Int})
    @assert precompile(Tuple{typeof(deindexify), Int, Int, Expr, Int, AVector, AVector})
    @assert precompile(Tuple{typeof(deindexify), Expr, Bool, Expr, Int64, AVector, AVector})
    @assert precompile(Tuple{typeof(deindexify), Nothing, Bool, Expr, Bool, AVector, AVector, Bool})
    @assert precompile(Tuple{typeof(deindexify), Symbol, Bool, Expr, Bool, AVector, AVector})
    @assert precompile(Tuple{typeof(disable_blas)})
    @assert precompile(Tuple{typeof(disable_cache)})
    @assert precompile(Tuple{typeof(enable_blas)})
    @assert precompile(Tuple{typeof(enable_cache)})
    @assert precompile(Tuple{typeof(expandconj), Expr})
    @assert precompile(Tuple{typeof(expandconj), Symbol})
    @assert precompile(Tuple{typeof(getallindices), Expr})
    @assert precompile(Tuple{typeof(getallindices), Int})
    @assert precompile(Tuple{typeof(getallindices), Symbol})
    @assert precompile(Tuple{typeof(geteltype), Expr})
    @assert precompile(Tuple{typeof(getindices), Symbol})
    @assert precompile(Tuple{typeof(getindices), Expr})
    @assert precompile(Tuple{typeof(getlhsrhs), Expr})
    @assert precompile(Tuple{typeof(hastraceindices), Expr})
    @assert precompile(Tuple{typeof(isassignment), Expr})
    @assert precompile(Tuple{typeof(isdefinition), Expr})
    @assert precompile(Tuple{typeof(isgeneraltensor), Expr})
    @assert precompile(Tuple{typeof(isindex), Symbol})
    @assert precompile(Tuple{typeof(isindex), Int})
    @assert precompile(Tuple{typeof(isnconstyle), Array{AVector, 1}})
    @assert precompile(Tuple{typeof(isscalarexpr), Expr})
    @assert precompile(Tuple{typeof(isscalarexpr), Float64})
    @assert precompile(Tuple{typeof(isscalarexpr), LineNumberNode})
    @assert precompile(Tuple{typeof(isscalarexpr), Symbol})
    @assert precompile(Tuple{typeof(istensor), Expr})
    @assert precompile(Tuple{typeof(istensorexpr), Expr})
    @assert precompile(Tuple{typeof(makegeneraltensor), Expr})
    @assert precompile(Tuple{typeof(makeindex), Int})
    @assert precompile(Tuple{typeof(makeindex), Symbol})
    @assert precompile(Tuple{typeof(makescalar), Expr})
    @assert precompile(Tuple{typeof(makescalar), Float64})
    @assert precompile(Tuple{typeof(maketensor), Expr})
    @assert precompile(Tuple{typeof(mulcost), Power{:χ, Int64}, Power{:χ, Int64}})
    @assert precompile(Tuple{typeof(ncontree), Vector{AVector}})
    @assert precompile(Tuple{typeof(optdata), Expr})
    @assert precompile(Tuple{typeof(optimaltree), Vector{AVector}, Base.Dict{Any, Power{:χ, Int64}}})
    @assert precompile(Tuple{typeof(parsecost), Expr})
    @assert precompile(Tuple{typeof(parsecost), Int64})
    @assert precompile(Tuple{typeof(parsecost), Symbol})
    @assert precompile(Tuple{typeof(processcontractorder), Expr, Nothing})
    @assert precompile(Tuple{typeof(processcontractorder), Expr, Int})
    @assert precompile(Tuple{typeof(processcontractorder), Int, Int})
    @assert precompile(Tuple{typeof(processcontractorder), Symbol, Nothing})
    @assert precompile(Tuple{typeof(storeset), Type{Base.BitArray{1}}, AVector, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{Base.BitArray{1}}, Array{Int64, 1}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{Base.BitArray{1}}, Base.Set{Int64}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{Base.BitSet}, AVector, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{Base.BitSet}, Array{Int64, 1}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{Base.BitSet}, Base.Set{Int64}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt128}, AVector, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt128}, Array{Int64, 1}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt128}, Base.Set{Int64}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt32}, AVector, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt32}, Array{Int64, 1}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt32}, Base.Set{Int64}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt64}, AVector, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt64}, Array{Int64, 1}, Int64})
    @assert precompile(Tuple{typeof(storeset), Type{UInt64}, Base.Set{Int64}, Int64})
    @assert precompile(Tuple{typeof(tensorify), Expr, Nothing})
    @assert precompile(Tuple{typeof(tensorify), Expr, Int})
    @assert precompile(Tuple{typeof(tensorify), Expr})
    @assert precompile(Tuple{typeof(tensorify), Int, Int})
    @assert precompile(Tuple{typeof(tree2expr), Int, Int})
    @assert precompile(Tuple{typeof(unique2), AVector})
    @assert precompile(Tuple{typeof(unique2), Array{Int64, 1}})
    @assert precompile(Tuple{typeof(use_blas)})
    @assert precompile(Tuple{typeof(use_cache)})
end
_precompile_()

end # module

ext: jl


text:
using NLopt
using Parameters

@with_kw type OptConfig
    LOWERXB = 1e-12
    XTOLREL = 1e-5
    XTOLABS = 0
    FTOLREL = 1e-5
    CTOLABS = 1e-5
    DEBUG   = false
    MAXEVAL = 0
    METHOD = :ineq
    OPTIMIZER = :LD_MMA
end

uniformweights(m::LikelihoodModel) = uniformweights(length(m.xs))
uniformweights(n) = fill(1/n, n)

function mple(m, reg, w0=uniformweights(m); config=OptConfig())
    @unpack_OptConfig config

    # generate closures
    f  = mple_obj(m,reg)
    df = dmple_obj(m,reg)

    function objective(x,g)
        # note the sign due to minimizing solver
        if length(g) == n
            g[:] = -df(x)
        end
        fx = f(x)
        DEBUG && @printf("f: f(x)=%f sum(x)=%f \n", fx, sum(x))
        -f(x)
    end

    local opt

    n = length(w0)
    opt = Opt(OPTIMIZER, n)

    min_objective!(opt, objective)
    lower_bounds!(opt, LOWERXB)
    upper_bounds!(opt, 1)

    xtol_rel!(opt, XTOLREL)
    xtol_abs!(opt, XTOLABS)
    ftol_rel!(opt, FTOLREL)

    maxeval!(opt, MAXEVAL)

    if METHOD == :ineq
        myineq(x, g)  = ((length(g) == n) && (g[:] = -1);
                         sum(x) - 1. - CTOLABS/2)
        myineq2(x, g) = ((length(g) == n) && (g[:] = -1);
                         1. - sum(x) - CTOLABS/2)

        inequality_constraint!(opt, myineq, CTOLABS/2)
        inequality_constraint!(opt, myineq2, CTOLABS/2)

    elseif METHOD == :auglag

        function eqconst(x, grad)
            fill!(grad, 1)
            sum(x) - 1
        end

        # store inner optimizer and overwrite with outer
        inneropt = opt
        opt  = Opt(:LD_AUGLAG, n)

        min_objective!(opt, objective)
        lower_bounds!(opt, LOWERXB)
        upper_bounds!(opt, 1)

        equality_constraint!(opt, eqconst, CTOLABS)

        xtol_rel!(opt, XTOLREL)
        xtol_abs!(opt, XTOLABS)
        ftol_rel!(opt, FTOLREL)

        local_optimizer!(opt, inneropt)
    end

    minf, minx, ret = optimize(opt, w0)

    if minx == w0
        warn("Optimizer returned initial proposal")
    end

    minx
end
#=
using JuMP
function jump(m,reg)
    n = length(m.xs)
    jm=Model()
    @variable(jm, 0 <= x[1:n] <= 1)

    f  = mple_obj(m,reg)
    df = dmple_obj(m,reg)

    JuMP.register(jm, :obj, n, (x...)->f(x), (g,x...)->(g[:] = df(x)))

    @constraint(jm, sum(x) == 1)
    @eval @NLobjective($jm, Max, $(Expr(:call, :obj, [Expr(:ref,:x, i) for i=1:n]...)))

    solve(jm)
end
=#

using MathProgBase

type StochasticProblem <: MathProgBase.AbstractNLPEvaluator
    n
    f
    df
end

SP  = StochasticProblem
MPB = MathProgBase

# initialization
MPB.initialize(d::SP, rf) = ()
MPB.features_available(d::SP) = [:Grad, :Jac]

# f
MPB.eval_f(d::SP, x) = d.f(x)
MPB.eval_grad_f(d::SP, g, x) = g[:]=d.df(x)

# constraints
MPB.eval_g(d::SP, g, x) = g[1]=sum(x)
MPB.jac_structure(d::SP) = ones(Int,d.n), collect(1:d.n)
MPB.eval_jac_g(d::SP, J, x) = J[:]=ones(Int, d.n)

MPB.isobjlinear(d::SP) = false
MPB.isobjquadratic(d::SP) = false
MPB.isconstrlinear(d::SP) = true

using Ipopt

function mpb(m, reg)
    nlm = MPB.NonlinearModel(IpoptSolver())
    n = length(m.xs)
    sp = SP(n, mple_obj(m,reg), dmple_obj(m,reg))
    MPB.loadproblem!(nlm, n, 1, fill(0, n), fill(1, n), [1], [1], :Max, sp)
    MPB.setwarmstart!(nlm, uniformweights(n))
    MPB.optimize!(nlm)
    @show MathProgBase.status(nlm)
end


text:
using My_iesh
using Test

# @testset "My_iesh.jl" begin
#     # Write your tests here.
#     @test f(2, 1) == 7
#     @test f(2, 3) == 13
# end
#
# @testset "oh hallo" begin
#     @test derivative_of_f(2, 1) == 2
end
#test NearestNeighbors
# plt3d= Plots.plot(testo3[:,1],testo3[:,2], testo3[:,3],
#      seriestype=:scatter, markersize = 7)
# display(plt3d)
#
#
# testo2 = similar(testo)
# testo2[:, 1] = testo[:, 1] .- testo[1,1]
# testo2[:, 2] = testo[:, 2] .- testo[1,2]
# testo2[:, 3] = testo[:, 3] .- testo[1,3]
#
# testo4 = similar(testo)
# for i in 1:13
#     testo3[i, :] = U' * testo2[i, :]
# end
# #
# # x = test_u2[:, 1]
# #
# # y = test_u2[:, 2]
# #
# # z = test_u2[:, 3]
# #
# # plt3d= Plots.plot(x,y, z,
# #      seriestype=:scatter, markersize = 7)
# # display(plt3d)
# #
# # theta = -π/4
# # rx = [[1 0 0]; [0 cos(theta) sin(theta)]; [0 -sin(theta) cos(theta)]]
# # ry = [[cos(theta) 0 -sin(theta)]; [0 1 0]; [sin(theta) 0 cos(theta)]]
# # r111 = rx * ry
# #
# #
# miu = zeros(Int32, 12)
# U = permutedims([[-1 0 1]/sqrt(2);[1 -2 1]/sqrt(6);[-1 -1 -1]/sqrt(3)])
# # for i in 1:13
#     for j in 2:13
#         temp = (testo[j, :] - testo[i, :])./abc
#         temp = temp .- floor.(temp .+ 1.0/2)
#         r = temp.*abc
#         rmi = floor(sum(r.*r)*10000)/10000
#         if rmi <= g^2
#             rb = U * r
#             # s = Integer.(sign.(rb) .* round.(abs.(rb)/g))
#             s = @. Integer(round(rb/g))
#             sint = Integer(s[1] * 9 + s[2] * 3 + s[3])
#             println(sint)
#             if sint == 4
#                 miu[1] = j
#             elseif sint == 2
#                 miu[2] = j
#             elseif sint == 10
#                 miu[3] = j
#             elseif sint == -8
#                 miu[4] = j
#             elseif sint == 12
#                 miu[5] = j
#             elseif sint == 6
#                 miu[6] = j
#             elseif sint == -4
#                 miu[7] = j
#             elseif sint == -2
#                 miu[8] = j
#             elseif sint == -10
#                 miu[9] = j
#             elseif sint == 8
#                 miu[10] = j
#             elseif sint == -12
#                 miu[11] = j
#             elseif sint == -6
#                 miu[12] = j
#             end
#         end
#     end
# # end
# #
# #         s = Integer.(sign.(rb) .* round.(rb/g))

#test for equilibrium: r0 =
# 1.476-1.476-0.000 2.952-1.476-1.476-1.476 1.476 0.000-2.952 1.476 1.476
# -0.852-2.557 1.704-0.000-0.852 2.557 0.852 2.557-1.704 0.000 0.852-2.557
# -2.410 0.000-2.410-0.000-2.410 0.000 2.410-0.000 2.410 0.000 2.410-0.000

# test_x = x0[100, :]
# test_r = r
#
# function test_potential(x)
#     V = 1.0/2.0*dot(x, d1_new * x)
#     return V
# end
#
#
#
# function test_force(x)
#     F = zeros(Float64, 3)
#     F = -1.0/2.0 * (d1_new + d1_new') * x
#     return F
# end
#
#
# a0 = test_force(test_r)
# q2 = grad(test_potential)
# a2 = -q2(r)
#
# q1 = x -> ForwardDiff.gradient(test_potential, x)
# a1 = -q1(r)
#
# a0./a2


text:
## External-viewer interface
function imshow(img, range)
    if ndims(img) == 2
        # only makes sense for gray scale images
        img = imadjustintensity(img, range)
    end
    tmp::String = "tmp.ppm"
    imwrite(img, tmp)
    cmd = `$imshow_cmd $tmp`
    spawn(cmd)
end

imshow(img) = imshow(img, [])

# 'illustrates' fourier transform
ftshow{T}(A::Array{T,2}) = imshow(log(1+abs(fftshift(A))),[])


text:
#=
Given a multiset of integers, return whether it can be partitioned into two subsets whose sums are the same.

For example, given the multiset {15, 5, 20, 10, 35, 15, 10}, it would return true, since we can split it up into {15, 5, 10, 15, 10} and {20, 35}, which both add up to 55.

Given the multiset {15, 5, 20, 10, 35}, it would return false, since we can't split it up into two subsets that add up to the same sum.
=#

using Test
include("Solutions/problem60_partition_multiset_in_two_equal_subset.jl")

@test ispartitionable([15, 5, 20, 10, 35, 15, 10])
@test !ispartitionable([15, 5, 20, 10, 35])

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
module SymSparseMatrixCSRTests

using SparseArrays
using SparseMatricesCSR
using LinearAlgebra
using Gridap.Algebra
using Test

maxnz=5
maxrows=5
maxcols=5
int_types=(Int32,Int)
float_types=(Float32,Float64)
Bi_types=(0,1)

for Ti in int_types
  for Tv in float_types
    for Bi in Bi_types
      I = Vector{Ti}()
      J = Vector{Ti}()
      V = Vector{Tv}()
      for (ik, jk, vk) in zip(rand(1:maxrows, maxnz), rand(1:maxcols, maxnz), rand(1:Tv(maxnz), maxnz-1) )
        push_coo!(SymSparseMatrixCSR{Bi,Tv,Ti},I,J,V,ik,jk,vk)
      end
      push_coo!(SymSparseMatrixCSR{Bi,Tv,Ti},I,J,V,maxrows,maxcols,maxnz)
      finalize_coo!(SymSparseMatrixCSR{Bi,Tv,Ti},I,J,V,maxrows, maxcols)
      SYMCSC = Symmetric(sparse(I, J, V, maxrows, maxcols),:U)
      SYMCSR = symsparsecsr(Val(Bi),I, J, V, maxrows, maxcols)

      @test is_entry_stored(SymSparseMatrixCSR{Bi,Tv,Ti},1,1)
      @test is_entry_stored(SymSparseMatrixCSR{Bi,Tv,Ti},1,2)
      @test !is_entry_stored(SymSparseMatrixCSR{Bi,Tv,Ti},2,1)

      _l = 10
      _I, _J, _V = allocate_coo_vectors(SymSparseMatrixCSR{Bi,Tv,Ti},_l)
      @test length(_I) == _l
      @test length(_J) == _l
      @test length(_V) == _l
      @test eltype(_I) == Ti
      @test eltype(_J) == Ti
      @test eltype(_V) == Tv

      V2 = copy(V)
      V2 .= 0
      SYMCSR2 = symsparsecsr(Val(Bi),I, J, V2, maxrows, maxcols)
      copy_entries!(SYMCSR2,SYMCSR)
      @test SYMCSR2 == SYMCSR
      copy_entries!(SYMCSR2,SYMCSR2)
      @test SYMCSR2 == SYMCSR

      @test size(SYMCSC)==size(SYMCSR)
      @test SYMCSC == SYMCSR

      @test convert(SymSparseMatrixCSR{Bi}, SYMCSR) === SYMCSR

      vold = getindex(SYMCSR,maxrows,maxcols)
      add_entry!(+,SYMCSR,1,maxrows,maxcols)
      @test getindex(SYMCSR,maxrows,maxcols) == vold+1

      LinearAlgebra.fillstored!(SYMCSR,0)
      @test all(x->x==0, nonzeros(SYMCSR))

    end
  end
end

end # module


text:
module MadelungConstant

import Plots

function Madelung(m)
	madelung = 0
	for i = -m:m
		for j = -m:m
			for k = -m:m
				if i == 0 && j == 0 && k == 0
					continue
				end
				madelung += ( (i + j + k) % 2 == 0 ? -1 : 1) / √(i^2 + j^2 + k^2)
			end
		end
	end
	madelung
end

function Madelung2(m)
	madelung = 0
	for i = -m:m
		for j = -m:m
			for k = -m:m
				if i == 0 && j == 0 && k == 0
					continue
				end
				r = √(i^2 + j^2 + k^2)
				if r > m
					continue
				end
				madelung += ((i + j + k) % 2 == 0 ? -1 : 1) / r
			end
		end
	end
	madelung
end

function plot(range=1:50)
    Plots.plot([Madelung, Madelung2], range, show=true)
end

end # module

max_stars_repo_name: YIsoda/MadelungConstant

max_issues_repo_name: YIsoda/MadelungConstant

max_forks_repo_name: YIsoda/MadelungConstant


text:
__precompile__()
module KaiserWindow_mod
#using Yeppp
export KaiserWindow

const XCUT=3.75
const NP=7
const NQ=9
const p=[1.0000000,
         3.5156229,
         3.0899424,
         1.2067492,
         0.2659732,
         0.03607680,
         0.004581300]

const q=[0.39894228,
         0.13285920e-1,
         0.22531900e-2,
        -0.15756500e-2,
         0.91628100e-2,
        -0.20577060e-1,
         0.26355370e-1,
        -0.16476330e-1,
         0.39237700e-2]

function bessel_i0(x::Float64)::Float64
ax::Float64=abs(x)
if ax < XCUT
   y::Float64=ax/XCUT
   y *= y
   poly::Float64=p[NP]
   for n=NP-1:-1:1
       poly=p[n]+y*poly
    end
else
        y = XCUT / ax;
        poly=q[NQ]
        for n=NQ-1:-1:1
            poly = q[n]+y*poly
        end
        poly *= exp(ax)/sqrt(ax)
end
return poly
end

function KaiserWindow(x::Float64,lx::Float64,b::Float64)::Float64
    a::Float64=x/lx
    if a <= 1.0
       return bessel_i0(b * sqrt(1.0 - a * a)) / bessel_i0(b);
    else
       return Float64(0.0)
    end
end

end


text:
module Essentials

using Reexport: @reexport

export dtype, kind, update!, reset!

function dtype end
function kind end
function update! end
function reset! end

include("QuantumOperators.jl")
include("QuantumNumbers.jl")
include("Spatials.jl")
include("DegreesOfFreedom.jl")
include("Frameworks.jl")
include("QuantumSystems.jl")

@reexport using .QuantumOperators
@reexport using .QuantumNumbers
@reexport using .Spatials
@reexport using .DegreesOfFreedom
@reexport using .Frameworks
@reexport using .QuantumSystems

end # module


text:
mesh_files = [("Abaqus", ".inp")]
for (file_type, ext) in mesh_files
    @testset "Import $file_type" begin
        for T in Floats
            @testset "$T" begin
                @testset "TriangleMesh" begin
                    ref_mesh = setup_TriangleMesh(T, UInt8)
                    mesh = import_mesh("./mesh/mesh_files/tri"*ext, T)
                    @test mesh.name == ref_mesh.name
                    for i in eachindex(ref_mesh.vertices)
                        @test mesh.vertices[i] ≈ ref_mesh.vertices[i]
                    end
                    @test mesh.polytopes == ref_mesh.polytopes
                    @test mesh.groups == ref_mesh.groups
                end

                @testset "QuadrilateralMesh" begin
                    ref_mesh = setup_QuadrilateralMesh(T, UInt8)
                    mesh = import_mesh("./mesh/mesh_files/quad"*ext, T)
                    @test mesh.name == ref_mesh.name
                    for i in eachindex(ref_mesh.vertices)
                        @test mesh.vertices[i] ≈ ref_mesh.vertices[i]
                    end
                    @test mesh.polytopes == ref_mesh.polytopes
                    @test mesh.groups == ref_mesh.groups
                end

                @testset "MixedPolygonMesh" begin
                    ref_mesh = setup_MixedPolygonMesh(T, UInt8)
                    mesh = import_mesh("./mesh/mesh_files/tri_quad"*ext, T)
                    @test mesh.name == ref_mesh.name
                    for i in eachindex(ref_mesh.vertices)
                        @test mesh.vertices[i] ≈ ref_mesh.vertices[i]
                    end
                    @test mesh.polytopes == ref_mesh.polytopes
                    @test mesh.groups == ref_mesh.groups
                end

                @testset "QuadraticTriangleMesh" begin
                    ref_mesh = setup_QuadraticTriangleMesh(T, UInt8)
                    mesh = import_mesh("./mesh/mesh_files/tri6"*ext, T)
                    @test mesh.name == ref_mesh.name
                    for i in eachindex(ref_mesh.vertices)
                        @test mesh.vertices[i] ≈ ref_mesh.vertices[i]
                    end
                    @test mesh.polytopes == ref_mesh.polytopes
                    @test mesh.groups == ref_mesh.groups
                end

                @testset "QuadraticQuadrilateralMesh" begin
                    ref_mesh = setup_QuadraticQuadrilateralMesh(T, UInt8)
                    mesh = import_mesh("./mesh/mesh_files/quad8"*ext, T)
                    @test mesh.name == ref_mesh.name
                    for i in eachindex(ref_mesh.vertices)
                        @test mesh.vertices[i] ≈ ref_mesh.vertices[i]
                    end
                    @test mesh.polytopes == ref_mesh.polytopes
                    @test mesh.groups == ref_mesh.groups
                end

                @testset "MixedQuadraticPolygonMesh" begin
                    ref_mesh = setup_MixedQuadraticPolygonMesh(T, UInt8)
                    mesh = import_mesh("./mesh/mesh_files/tri6_quad8"*ext, T)
                    @test mesh.name == ref_mesh.name
                    for i in eachindex(ref_mesh.vertices)
                        @test mesh.vertices[i] ≈ ref_mesh.vertices[i]
                    end
                    @test mesh.polytopes == ref_mesh.polytopes
                    @test mesh.groups == ref_mesh.groups
                end
            end
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using Falcons
using Healpix
#using Plots
using NPZ
#using StatsBase
using BenchmarkTools
using PyCall
using TickTock
using Base.Threads
using StaticArrays
using PyPlot
using Formatting
using ReferenceFrameRotations
using LinearAlgebra
using ProgressMeter
using HDF5
using Interpolations

hp = pyimport("healpy")
np = pyimport("numpy")
pd = pyimport("pandas")

include("../convolvor_function.jl")

data = npzread("/home/cmb/yusuket/program/beam_study/npz_beamdata/MFT_88.5GHz_000.0_166.7_xpol_v2.npz")
Nx = 360
Ny = 721
beam2d, x, y = get_2Ddata(data, Nx, Ny)

res = Resolution(512)
cmb_path = "/home/cmb/yusuket/program/MapData/Nside512/lensed_r0_512_non_smooth.fits"
skymap = hp.read_map(cmb_path, field=(0,1,2));
lmax = 3res.nside - 1
alm = hp.map2alm(skymap, lmax = lmax, datapath = "/home/cmb/yusuket/program/MapData/healpy-data/.", use_pixel_weights=true)
skyalm = AlmPair(alm, conj.(alm))

idx = parse(Int64, ARGS[1])
convolved_map = EffectiveBeamConvolution(res, skyalm, beam2d, idx)

println("Saving...")
np.savez_compressed("/group/cmb/litebird/usr/ytakase/beam_convolution/devided_maps/outputs/mapdata/idx=$idx", convolved_map = convolved_map)
println("=====Finish=====")


text:
using AVLTrees, BenchmarkTools
using Random
using Plots
using DataFrames

function batch_insert!(t::AVLTree{K,D}, v::Vector{K}) where {K,D}
    for i in v
        insert!(t, i, i)
    end
end

function batch_delete!(t::AVLTree{K,D}, v::Vector{K}) where {K,D}
    for i in v
        insert!(t, i, i)
    end
end

function batch_find(t::AVLTree{K,D}, v::Vector{K}) where {K,D}
    for i in v
        findkey(t, i)
    end
end

insertion_vec = []
deletion_vec = []
search_vec = []

d = DataFrame((op=[], time=[], n=[]))
x = [1000, 10000, 100000, 1000000]

for attempt in 1:20
for N in x
    global t = AVLTree{Int64,Int64}()
    rng = MersenneTwister(1111)
    global nums = rand(rng, Int64, N)
    global unique_nums = Vector{Int64}(undef, 1000)
    for i in pairs(unique_nums)
        r = rand(rng, Int64)
        while (r in nums)
            r = rand(rng, Int64)
        end
        unique_nums[i[1]] = r
    end
    for i in nums
        insert!(t, i, i)
    end

    search = @benchmark batch_find(t, nums[1:1000]) samples=1 evals=1
    insertion = @benchmark batch_insert!(t, unique_nums) samples=1 evals=1

    t = AVLTree{Int64,Int64}()
    rng = MersenneTwister(1111)
    nums = rand(rng, Int64, N)
    unique_nums = Vector{Int64}(undef, 1000)
    for i in pairs(unique_nums)
        r = rand(rng, Int64)
        while (r in nums)
            r = rand(rng, Int64)
        end
        unique_nums[i[1]] = r
    end
    for i in nums
        insert!(t, i, i)
    end

    deletion = @benchmark batch_delete!(t, nums[1:1000]) samples=1 evals=1

    push!(d, ("insert", minimum(insertion).time, N))
    push!(d, ("delete", minimum(deletion).time,N))
    push!(d, ("search", minimum(search).time,N))
    println("done $N")
end
end

c= combine(groupby(d, [:op,:n]), :time => minimum)

# plot(x, insertion_vec/1000, xscale=:log10, ylabel="us")
# plot(x, deletion_vec/1000, xscale=:log10, ylabel="us")
# plot(x, search_vec/1000, xscale=:log10, ylabel="us")

plot(
    x,
    [c[(c.op.=="insert"),:].time_minimum./1000,c[(c.op.=="delete"),:].time_minimum./1000, c[(c.op.=="search"),:].time_minimum./1000],
    xscale = :log10,
    ylabel = "operation time [us]",
    xlabel = "N",
    markershape =[:diamond :utriangle :dtriangle],
    labels= ["insert" "delete" "lookup"],
    legend=:topleft,
)

#savefig("result.svg")


text:
#!/usr/bin/env julia

# Reading from multiple sockets

# The ZMQ.jl wrapper implements ZMQ.recv as a blocking function. Nonblocking i/o
# in Julia is typically done using coroutines (Tasks).
# The @async macro puts its enclosed expression in a Task. When the macro is
# executed, its Task gets scheduled and execution continues immediately to
# whatever follows the macro.

# Note: the msreader example in the zguide is presented as a "dirty hack"
# using the ZMQ_DONTWAIT and EAGAIN codes. Since the ZMQ.jl wrapper API
# does not expose DONTWAIT directly, this example skips the hack and instead
# provides an efficient solution.

using ZMQ

# Prepare our context and sockets
context = ZMQ.Context()

# Connect to task ventilator
receiver = Socket(context, ZMQ.PULL)
ZMQ.connect(receiver, "tcp://localhost:5557")

# Connect to weather server
subscriber = Socket(context,ZMQ.SUB)
ZMQ.connect(subscriber,"tcp://localhost:5556")
ZMQ.set_subscribe(subscriber, "10001")

while true

    # Process any waiting tasks
    @async begin
        msg = unsafe_string(ZMQ.recv(receiver))
        println(msg)
    end

    # Process any waiting weather updates
    @async begin
        msg = unsafe_string(ZMQ.recv(subscriber))
        println(msg)
    end

    # Sleep for 1 msec
    sleep(0.001)
end

max_stars_repo_name: JailbreakFox/LightWeightRepository

max_issues_repo_name: JailbreakFox/LightWeightRepository

max_forks_repo_name: JailbreakFox/LightWeightRepository


text:
# ---
# title: 641. Design Circular Deque
# id: problem641
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Design, Queue
# link: <https://leetcode.com/problems/design-circular-deque/description/>
# hidden: true
# ---
#
# Design your implementation of the circular double-ended queue (deque).
#
# Your implementation should support following operations:
#
#   * `MyCircularDeque(k)`: Constructor, set the size of the deque to be k.
#   * `insertFront()`: Adds an item at the front of Deque. Return true if the operation is successful.
#   * `insertLast()`: Adds an item at the rear of Deque. Return true if the operation is successful.
#   * `deleteFront()`: Deletes an item from the front of Deque. Return true if the operation is successful.
#   * `deleteLast()`: Deletes an item from the rear of Deque. Return true if the operation is successful.
#   * `getFront()`: Gets the front item from the Deque. If the deque is empty, return -1.
#   * `getRear()`: Gets the last item from Deque. If the deque is empty, return -1.
#   * `isEmpty()`: Checks whether Deque is empty or not.
#   * `isFull()`: Checks whether Deque is full or not.
#
#
#
# **Example:**
#
#
#
#     MyCircularDeque circularDeque = new MycircularDeque(3); // set the size to be 3
#     circularDeque.insertLast(1);			// return true
#     circularDeque.insertLast(2);			// return true
#     circularDeque.insertFront(3);			// return true
#     circularDeque.insertFront(4);			// return false, the queue is full
#     circularDeque.getRear();  			// return 2
#     circularDeque.isFull();				// return true
#     circularDeque.deleteLast();			// return true
#     circularDeque.insertFront(4);			// return true
#     circularDeque.getFront();			// return 4
#
#
#
#
# **Note:**
#
#   * All values will be in the range of [0, 1000].
#   * The number of operations will be in the range of [1, 1000].
#   * Please do not use the built-in Deque library.
#
#
## @lc code=start
using LeetCode

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


text:
module test_julia_cartesian

using Test

start = CartesianIndex(1, 2, 5)
stop = CartesianIndex(1, 3, 7)
range = CartesianIndices(start:stop)
@test range isa CartesianIndices{3}
@test length(range) == 6
@test range[:, :, 1] == [CartesianIndex(1, 1, 1) CartesianIndex(1, 2, 1)]
@test range[:, :, 2] == [CartesianIndex(1, 1, 2) CartesianIndex(1, 2, 2)]
@test range[:, :, 3] == [CartesianIndex(1, 1, 3) CartesianIndex(1, 2, 3)]
@test range[:, 1, :] == [CartesianIndex(1, 1, 1) CartesianIndex(1, 1, 2) CartesianIndex(1, 1, 3)]
@test range[:, 2, :] == [CartesianIndex(1, 2, 1) CartesianIndex(1, 2, 2) CartesianIndex(1, 2, 3)]
@test range[1] == CartesianIndex(1, 1, 1)
@test range[2] == CartesianIndex(1, 2, 1)
@test range[3] == CartesianIndex(1, 1, 2)
@test range[4] == CartesianIndex(1, 2, 2)
@test range[5] == CartesianIndex(1, 1, 3)
@test range[6] == CartesianIndex(1, 2, 3)

end # module test_julia_cartesian


text:
module Uhlarm

using PortAudio
using LibSndFile

export @uhlarm

const SOUNDPATH = joinpath(dirname(@__FILE__()), "..", "sounds")

macro uhlarm(expr)
    quote
        output = try
            $(esc(expr))
        catch ex
            playrand("bad")
            rethrow(ex)
        end
        playrand("good")
        return output
    end
end

function playrand(dir="good")
    absdir = joinpath(SOUNDPATH, dir)
    files = readdir(absdir)
    path = joinpath(absdir, rand(files))
    write(PortAudioStream(), load(path))
end

end # module


text:
module AllenBrain

using Statistics
using JSON, LightGraphs, IndirectArrays, AxisArrays, OffsetArrays
using StaticArrays, CoordinateTransformations
using ImageMagick, ImageTransformations, ImageCore, Interpolations
using ColorVectorSpace
using IntervalSets, ProgressMeter, FileIO
using Unitful: μm
using Requires
import HTTP

export #
    # basic types
    BoundingBox,
    buffer,
    # ontology
    ontology,
    findvertices,
    structureids,
    boundingbox,
    # image loading, query, and download
    annotation,
    download_annotation,
    download_projection,
    query_reference2image,
    query_image2reference,
    sectionimage,
    download_insitu_images,
    splice_sectionimages,
    # genes
    query_insitu,
    # visualization
    colorize,
    colorize!,
    visualize_volume

const mousedir = joinpath(dirname(@__DIR__), "data", "mouse")

function dataset(species, category)
    if category == "ontology"
        species == "mouse" && return joinpath(mousedir, "structure_graph.json")
    elseif startswith(category, "annotation") && species == "mouse"
        if !endswith(category, ".nrrd")
            if !ismatch(r"[0-9]$", category)
                category = category*"_25"
            end
            category = category*".nrrd"
        end
        return joinpath(mousedir, category)
    end
    error(category, " for species ", species, " not found")
end

@inline mapfilter(f, t, ::Type{Axis{name}}, axs) where name =
    _mapfilter(f, t, Axis{name}, axs...)
@inline _mapfilter(f, t, ::Type{Axis{name}}, ::Axis{name}, axs...) where name =
    (f(t[1]), _mapfilter(f, Base.tail(t), Axis{name}, axs...)...)
@inline _mapfilter(f, t, ::Type{Axis{name}}, ::Axis, axs...) where name =
    (t[1], _mapfilter(f, Base.tail(t), Axis{name}, axs...)...)
_mapfilter(f, ::Tuple{}, ::Type{Axis{name}}) where name = ()

inmicrons(x::Real) = x
inmicrons(x) = x/(1μm)

include("types.jl")
include("ontology.jl")
include("images.jl")
include("projections.jl")
include("genes.jl")
include("visualize.jl")

function __init__()
    # NOTE: this was written against GLVisualize which is deprecated for
    # Makie. The code below likely needs updating.
    @require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" begin
        function visualize_volume(volumedata)
            window = Makie.glscreen()
            volume = Makie.visualize(volumedata, :absorption)
            Makie._view(volume, window)
            @schedule Makie.renderloop(window)
        end
    end
end

end # module


text:
module printerpolate

# Pkg dependencies
using DelimitedFiles, Gadfly

# Includes
include("functions.jl")

# exports
export interpolate
```
    interpolate(x, data): x is a floating number, data must be an array or a file of array of numbers. The function reads first two columns as x and y values.
```
function interpolate(x::Y, data::T) where{Y, T}
    if isa(data, String)
        xs = readdlm(data)[:,1]
        ys = readdlm(data)[:,2]
    elseif isa(data, Array)
        xs = data[:,1]
        ys = data[:,2]
    end
    return lagrangeinterpolation(x,xs,ys)
end

end # module


text:
export costfunction, LHD

# "Create tools to compute cost function"

# Type to hold cost functions
abstract type costfunction end

# Likelihood function
struct LHD <:costfunction
	A::Array{Float64,2}
	dψ::Array{Float64,2}
	b::Array{Float64,1}
	λ::Float64
	δ::Float64
	nb::Int64
	Ne::Int64

	# Cache variables
	dS::Array{Float64,1}
	logdS::Array{Float64,1}
	dψ_dS::Array{Float64,2}

	# Store gradient and hessian
	J::Array{Float64,1}
	G::Array{Float64,1}
	H::Array{Float64,2}
end

## Create constructor for likelihood function

function LHD(A::Array{Float64,2}, dψ::Array{Float64,2}, λ::Float64, δ::Float64)
	# Number of basis = p+2
	nb = size(A,1)
	@assert size(A,1)==size(A,2) "A must be a square matrix"

	Ne = size(dψ, 1)

	# Add L-2 regularization
	A += (λ/Ne)*I

	b = zeros(nb)
	b .= sum(A, dims=2)[:,1]
	rmul!(b, δ)

	return LHD(A, dψ, b, λ, δ, nb, Ne, zeros(Ne), zeros(Ne), zeros(Ne,nb),zeros(1), zeros(nb), zeros(nb,nb))
end

function Base.show(io::IO, L::LHD)
	if L.λ==0.0 && L.δ==0.0
	println(io,"Cost function with Ne = $(L.Ne) samples")
	elseif L.λ==0.0 && L.δ!=0.0
	println(io,"Cost function with δ = $(L.δ) - regularized log barrier and Ne = $(L.Ne) samples")
	elseif L.λ!=0.0 && L.δ==0.0
	println(io,"Cost function with λ = $(L.λ) L-2 regularization and Ne = $(L.Ne) samples")
	else
	println(io,"Cost function with λ = $(L.λ) L-2 regularization,  δ = $(L.δ) - regularized log barrier and Ne = $(L.Ne) samples")
	end
end

# Function to compute the cost, gradient and loss
function (lhd::LHD)(x::Array{Float64,1}, inplace::Bool; noutput::Int64=3)
	lhd.dS .= lhd.dψ*x
	lhd.dS .+= lhd.δ*sum(lhd.dψ, dims = 2)[:,1]

	@inbounds for i=1:lhd.Ne
		lhd.logdS[i] = log(lhd.dS[i])
	end

	if inplace
		# Compute gradient
		@inbounds for i=1:lhd.Ne
			lhd.dψ_dS[i,:] .= lhd.dψ[i,:]/lhd.dS[i]
		end

		# Return objective
		lhd.J .= 0.5*dot(x, lhd.A*x) - sum(lhd.logdS,dims = 1)[1]*(1/lhd.Ne) + dot(x,lhd.b)

		# Return gradient
		lhd.G .= lhd.A*x - sum(lhd.dψ_dS, dims = 1)[1,:]*(1/lhd.Ne) + lhd.b

		# Return hessian
		lhd.H .= lhd.A + BLAS.gemm('T', 'N', 1/lhd.Ne, lhd.dψ_dS, lhd.dψ_dS)
	else
		if noutput == 1
			# Return objective
			J = 0.5*dot(x, lhd.A*x) - sum(lhd.logdS,dims = 1)[1]*(1/lhd.Ne) + dot(x,lhd.b)
		return J

		elseif noutput == 2
			# Compute gradient
			@inbounds for i=1:lhd.Ne
				lhd.dψ_dS[i,:] .= lhd.dψ[i,:]/lhd.dS[i]
			end

			# Return objective
			J = 0.5*dot(x, lhd.A*x) - sum(lhd.logdS,dims = 1)[1]*(1/lhd.Ne) + dot(x,lhd.b)

			# Return gradient
			G = lhd.A*x - sum(lhd.dψ_dS, dims = 1)[1,:]*(1/lhd.Ne) + lhd.b
			return J,G

		elseif noutput == 3
			# Compute gradient
			@inbounds for i=1:lhd.Ne
				lhd.dψ_dS[i,:] .= lhd.dψ[i,:]/lhd.dS[i]
			end

			# Return objective
			J = 0.5*dot(x, lhd.A*x) - sum(lhd.logdS,dims = 1)[1]*(1/lhd.Ne) + dot(x,lhd.b)

			# Return gradient
			G = lhd.A*x - sum(lhd.dψ_dS, dims = 1)[1,:]*(1/lhd.Ne) + lhd.b

			# Return hessian
			H = lhd.A + BLAS.gemm('T', 'N', 1/lhd.Ne, lhd.dψ_dS, lhd.dψ_dS)

			return J,G,H
		end
	end
end


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.0
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

print("Hello world")

using Plots

# # Rapid Visualization on Jupyter
#
# Just `add IJulia.clear_output(true)` and display `p` plot object simultaneously will improve your visualizatoin life

ps=[]
for t in 1:0.1:10
    IJulia.clear_output(true)
    p = plot(x-> sin(x+t))
    p |> display
    push!(ps, p)
end

# # Dump result as tmp.gif

@gif for p in ps
    plot(p)
end


text:
map_to_extern(extern_func::Ptr{wasm_func_t}) = wasm_func_as_extern(extern_func)
map_to_extern(other) = error("Type $(typeof(other)) is not supported")

mutable struct WasmInstance
    wasm_instance_ptr::Ptr{wasm_instance_t}
    wasm_module::WasmModule
    wasm_externs::WasmVec{wasm_extern_vec_t,Ptr{wasm_extern_t}}

    WasmInstance(
        wasm_instance_ptr::Ptr{wasm_instance_t},
        wasm_module::WasmModule,
        wasm_externs = WasmVec{wasm_extern_vec_t,Ptr{wasm_extern_t}}(0, C_NULL),
    ) = finalizer(wasm_instance_delete, new(wasm_instance_ptr, wasm_module, wasm_externs))
end
function WasmInstance(store::WasmStore, wasm_module::WasmModule)
    module_imports = imports(wasm_module)
    n_expected_imports = length(module_imports.wasm_imports)
    @assert n_expected_imports == 0 "No imports provided, expected $n_expected_imports"

    empty_imports = WasmPtrVec(wasm_extern_t)
    # TODO: pass trap here to get better error messages
    # See `wasmtime/wasi.jl` for an example
    wasm_instance_ptr = wasm_instance_new(
        store.wasm_store_ptr,
        wasm_module.wasm_module_ptr,
        empty_imports,
        C_NULL,
    )
    @assert wasm_instance_ptr != C_NULL "Failed to create WASM instance"
    WasmInstance(wasm_instance_ptr, wasm_module)
end
function WasmInstance(store::WasmStore, wasm_module::WasmModule, host_imports)
    externs_vec = WasmPtrVec(map(map_to_extern, host_imports))
    WasmInstance(store, wasm_module, externs_vec)
end
function WasmInstance(
    store::WasmStore,
    wasm_module::WasmModule,
    externs_vec::WasmVec{wasm_extern_vec_t,Ptr{wasm_extern_t}},
)
    module_imports = imports(wasm_module)
    n_expected_imports = length(module_imports.wasm_imports)
    n_provided_imports = length(externs_vec)
    @assert n_expected_imports == n_provided_imports "$n_provided_imports imports provided, expected $n_expected_imports"

    wasm_trap_ptr = Ref(Ptr{wasm_trap_t}())
    wasm_instance_ptr = wasm_instance_new(
        store.wasm_store_ptr,
        wasm_module.wasm_module_ptr,
        externs_vec,
        wasm_trap_ptr,
    )
    if wasm_instance_ptr == C_NULL && wasm_trap_ptr[] != C_NULL
        trap_msg = WasmByteVec()
        wasm_trap_message(wasm_trap_ptr[], trap_msg)
        wasm_trap_delete(wasm_trap_ptr[])

        error_message = unsafe_string(trap_msg.data, trap_msg.size)
        error(error_message)
    end
    @assert wasm_instance_ptr != C_NULL "Failed to create WASM instance"
    WasmInstance(wasm_instance_ptr, wasm_module, externs_vec)
end

Base.show(io::IO, ::WasmInstance) = print(io, "WasmInstance()")
Base.unsafe_convert(::Type{Ptr{wasm_instance_t}}, instance::WasmInstance) =
    instance.wasm_instance_ptr

# TODO: One for each exporttype_type?
mutable struct WasmExport <: AbstractWasmExport
    # The wasm_exporttype_t refers to the export on the module side
    wasm_export_ptr::Ptr{wasm_exporttype_t}

    # The wasm_extern_t refers to the export on the instance side
    wasm_extern::WasmExtern
    wasm_instance::WasmInstance
    name::String

    function WasmExport(
        wasm_export_ptr::Ptr{wasm_exporttype_t},
        wasm_extern_ptr::Ptr{wasm_extern_t},
        wasm_instance::WasmInstance,
    )
        owned_wasm_export_ptr = wasm_exporttype_copy(wasm_export_ptr)
        @assert owned_wasm_export_ptr != C_NULL "Failed to copy WASM export"

        owned_wasm_extern_ptr = wasm_extern_copy(wasm_extern_ptr)
        @assert owned_wasm_extern_ptr != C_NULL "Failed to copy WASM extern"

        name_vec_ptr = wasm_exporttype_name(owned_wasm_export_ptr)
        name_vec = Base.unsafe_load(name_vec_ptr)
        name = unsafe_string(name_vec.data, name_vec.size)
        wasm_name_delete(name_vec_ptr)

        # TODO: Extract type here
        finalizer(
            new(
                owned_wasm_export_ptr,
                WasmExtern(owned_wasm_extern_ptr),
                wasm_instance,
                name,
            ),
        ) do wasm_export
            wasm_extern_delete(wasm_export.wasm_extern)
            wasm_exporttype_delete(wasm_export.wasm_export_ptr)
        end
    end
end

function (wasm_export::WasmExport)(args...)
    wasm_externtype_ptr = wasm_exporttype_type(wasm_export.wasm_export_ptr)
    @assert wasm_externtype_ptr != C_NULL "Failed to get export type for export $(wasm_export.name)"
    wasm_externkind = wasm_externtype_kind(wasm_externtype_ptr)
    @assert wasm_externkind == WASM_EXTERN_FUNC "Called export '$(wasm_export.name)' is not a function"

    wasm_export.wasm_extern(args...)
end

function exports(wasm_instance::WasmInstance)
    exports = WasmPtrVec(wasm_exporttype_t)
    wasm_module_exports(wasm_instance.wasm_module.wasm_module_ptr, exports)
    externs = WasmPtrVec(wasm_extern_t)
    wasm_instance_exports(wasm_instance.wasm_instance_ptr, externs)
    @assert length(exports) == length(externs)

    exports_vector = map(a -> WasmExport(a..., wasm_instance), zip(exports, externs))

    WasmExports{WasmInstance,WasmExport}(wasm_instance, exports_vector)
end


text:
using ChemostatPlots
using Test

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


text:
include("../benchmark.jl")
include("create_report.jl")

results = get_CUTEst_results()
#["one_phase"] = load("../results/one_phase/Dec10/summary.jld", "summary")

using DataFrames,CSV

df_results = Dict()
for (method_name,method_results) in results
  df_results[method_name] = DataFrame(name=String[], number_variables=Int[], number_constraints=Int[],it=Int[], time=Array{Float64,1}(),fval=Array{Float64,1}(), con=Array{Float64,1(), total_fval_evaluation=Int[], total_grad_evaluation=Int[], total_jac_evaluation=Int[], total_cons_evaluation=Int[], total_hess_evaluation=Int[], status=[])
  for (name,info) in method_results
    push!(df_results[method_name],[name,info.number_variables,info.number_constraints,info.it_count,info.total_time,info.fval,info.con_vio,info.total_fval_evaluation,info.total_grad_evaluation,info.total_jac_evaluation,info.total_cons_evaluation,info.total_hess_evaluation,string(info.status)])
  end
  CSV.write("$folder/table_$(folder)_$(method_name).csv", df_results[method_name])
end

###
### ITER/TIMING BREAK DOWN
###

df_summary = DataFrame(name=[],median_it=[],median_time=[],total_time=[])
for (method_name,method_results) in results
    push!(df_summary,[method_name,median(df_results[method_name][:it]),median(df_results[method_name][:time]),sum(df_results[method_name][:time])])
end

df_summary

CSV.write("$folder/runtime_summary.csv", df_summary)

using FreqTables

##
## FAILURE REASONS
##
freqtable(df_results["ipopt"][:status])
freqtable(df_results["one phase"][:status])

max_stars_repo_name: ohinder/OnePhase

max_issues_repo_name: ohinder/OnePhase

max_forks_repo_name: ohinder/OnePhase


text:
# using Revise
using Test, BifurcationKit, ForwardDiff, RecursiveArrayTools, LinearAlgebra
const BK = BifurcationKit
N = 10
M = 5
δ = 1e-8
####################################################################################################
# test the jacobian of the multiple shooting functional using Linear flow
# TODO do example with A matrix and exp(At)
vf(x, p) = x
flow(x, p, t) = (u=exp(t) .* x, t=t)
dflow(x, p, dx, t) = (flow(x, p, t)..., du = exp(t) .* dx)
section(x, T) = dot(x[1:N], ones(N))
section(x, T, dx, dT) = dot(dx[1:N], ones(N)) * T
section(x::BorderedArray, T) = section(vec(x.u[:,:]), T)
section(x::BorderedArray, T, dx, dT) = section(vec(x.u[:,:]), T, vec(dx.u[:,:]), dT)
par = nothing

fl = BK.Flow(vf, flow, dflow)

probSh = BK.ShootingProblem(M, fl,
	LinRange(0, 1, M+1) |> diff,
	section, false)

show(probSh)

poguess = VectorOfArray([rand(N) for ii=1:M])
	po = BorderedArray(poguess, 1.)

dpoguess = VectorOfArray([rand(N) for ii=1:M])
	dpo = BorderedArray(dpoguess, 2.)

# use of AbstractArray structure
pov = vcat(vec(po.u), po.p)
dpov = vcat(vec(dpo.u), dpo.p)
resv = probSh(pov, par)

dresv = probSh(pov, par, dpov; δ = δ)
resad = ForwardDiff.derivative(t -> probSh(pov .+ t .* dpov, par), 0.)
@test norm(resad[1:end-1] - dresv[1:end-1], Inf) < 1e-14

# use of BorderedArray structure
res = probSh(po, par)
@test norm(vec(res.u[:,:]) - resv[1:end-1] , Inf) ≈ 0
@test norm(res.p - resv[end], Inf) ≈ 0

dres = probSh(po, par, dpo; δ = δ)
@test norm(dres.p - dresv[end], Inf) ≈ 0
@test norm(vec(dres.u[:,:]) - dresv[1:end-1], Inf) ≈ 0
####################################################################################################
# test the jacobian of the multiple shooting functional using nonLinear flow
vf(x, p) = x.^2
flow(x, p, t) = (u=x ./ (1 .- t .* x),t=t)
dflow(x, p, dx, t) = (flow(x, p, t)..., du = dx ./ (1 .- t .* x).^2)

fl = BK.Flow(vf, flow, dflow)

probSh = BK.ShootingProblem(M, fl,
	LinRange(0,1,M+1) |> diff ,
	section, false)

poguess = VectorOfArray([rand(N) for ii=1:M])
	po = BorderedArray(poguess, 1.)

dpoguess = VectorOfArray([rand(N) for ii=1:M])
	dpo = BorderedArray(dpoguess, 2.)

# use of AbstractArray structure
pov = vcat(vec(po.u), po.p)
dpov = vcat(vec(dpo.u), dpo.p)
resv = probSh(pov, par)

dresv = probSh(pov, par, dpov; δ = δ)
resad = ForwardDiff.derivative(t -> probSh(pov .+ t .* dpov, par), 0.)
@test norm(resad[1:end-1] - dresv[1:end-1], Inf) < 1e-14

# use of BorderedArray structure
res = probSh(po, par)
@test norm(vec(res.u[:,:]) - resv[1:end-1] , Inf) ≈ 0
@test norm(res.p - resv[end], Inf) ≈ 0

dres = probSh(po, par, dpo; δ = δ)
@test norm(dres.p - dresv[end], Inf) ≈ 0
@test norm(vec(dres.u[:,:]) - dresv[1:end-1], Inf) ≈ 0
####################################################################################################
# test the hyperplane projections for Poincare Shooting
M = 1
normals = [rand(2) for ii=1:M]
for ii=1:M
	normals[ii] /= norm(normals[ii])
end
centers = [rand(2) for ii=1:M]

hyper = BK.SectionPS(normals, centers)
poguess = VectorOfArray([rand(2) for ii=1:M])

x = rand(2)
xb = BK.R(hyper, x, 1)
# test
for ii=1:M
	x2 = BK.E(hyper, xb, ii)
	# test that x2 in Sigma2
	@test dot(x2 - centers[ii], normals[ii]) < 1e-14
	# that we have Rk∘Ek = Id and Ek∘Rk = IdΣ
	@test BK.R(hyper, x2, ii) ≈ xb
end

# test of the derivatives of E and R
dx = rand(2)
_out1 = (BK.R(hyper, x .+ δ .* dx, 1) - BK.R(hyper, x, 1)) ./ δ
_out2 = zero(_out1)
_out2 = BK.dR!(hyper, _out2, dx, 1)

@test norm(_out2 - _out1, Inf) < 1e-5

dx = rand(2-1)
_out1 = ForwardDiff.derivative(t -> BK.E(hyper, xb .+ t .* dx,1), 0.)
_out2 = BK.dE(hyper, dx, 1)
@test norm(_out2 - _out1, Inf) < 1e-12

# flow for Poincare return map
# we consider the ODE dr = r⋅(1-r), dθ = 1 [2π]
# the flow is given by Φ(r0,θ0,t) = (r0 e^{t}/(1-r0+r0 e^{t}),θ+t)
vf(x, p) = [x[1]*(1-x[1]), 1.]
flow(x, p, t; k...) = (t = t, u = [exp(t) .* x[1] / (1-x[1]+x[1]*exp(t)),x[2]+t])
Π2(x, p, t = 0) = (t = 2pi -x[2], u = flow(x, p, 2pi-x[2]).u)
Π(x, p, t = 0) = flow(x, p, 2pi-x[2])	# return map
dflow(x, p, dx, t; k...) = (t = t, u = flow(x, p, t).u, du = ForwardDiff.derivative( z -> flow(x .+ z .* dx, p, t).u, 0),)
section(x, T) = dot(x[1:2], [1, 0])
# section(x::BorderedArray, T) = section(vec(x.u[:,:]), T)
par = nothing

fl = BK.Flow(F = vf, flow = Π, flowSerial = Π2, dfSerial = dflow)
sectionps = SectionPS(normals, centers)
probPSh = PoincareShootingProblem(flow = fl, M = M, section = sectionps)

ci = reduce(vcat, BK.projection(probPSh, poguess.u))
dci = rand(length(ci))

# we test that we have the analytical version of the flow
z0 = rand(2)
@test ForwardDiff.derivative(z -> flow(z0, par, z).u, 0.) ≈ vf(z0, par)

#test the show method
show(probPSh)

# test functional and its differential
probPSh(ci, par)
probPSh(ci, par, dci)

# test the analytical of the differential of the return Map
z0 = rand(2)
dz0 = rand(2)
_out0 = BK.diffPoincareMap(probPSh, z0, par, dz0, 1)
_out1 = ForwardDiff.derivative(z -> Π(z0 .+ z .* dz0, par).u, 0)
display(_out0)
display(_out1)

# test the analytical version of the functional
_out0 = probPSh(ci, par, dci)
	δ = 1e-6
	_out2 = (probPSh(ci .+ δ .* dci, par) .- probPSh(ci, par)) ./ δ
	_out1 = ForwardDiff.derivative(z -> probPSh(ci .+ z .* dci, par), 0)
	display(_out0)
	display(_out1)
	display(_out2)
	_out0 - _out1 |> display


text:
"""
    read(A::AbstractGeoArray)
    read(A::AbstractGeoStack)
    read(A::AbstractGeoSeries)

`read` will move a GeoData.jl object completely to memory.
"""
function Base.read(x::Union{AbstractGeoArray,AbstractGeoStack,AbstractGeoSeries})
    modify(x) do ds
        # Some backends don't implement `Array` properly.
        Array{eltype(ds),ndims(ds)}(undef, size(ds)) .= ds
    end
end

"""
    read!(src::Union{AbstractString,AbstractGeoArray}, dst::AbstractGeoArray)
    read!(src::Union{AbstractString,AbstractGeoStack}, dst::AbstractGeoStack)
    read!(scr::AbstractGeoSeries, dst::AbstractGeoSeries)

`read!` will copy the data from `src` to the object `dst`.

`src` can be an object or a file-path `String`.
"""
Base.read!(src::AbstractGeoArray, dst::AbstractArray) = dst .= src
function Base.read!(src::AbstractGeoStack, dst::AbstractGeoStack)
    map(k -> read!(src[k], dst[k]), keys(dst))
    return dst
end
function Base.read!(src::AbstractGeoSeries, dst::AbstractGeoSeries)
    map(read!, src, dst)
    return dst
end

# Filename methods
function Base.read!(filename::AbstractString, dst::AbstractGeoArray)
    src = GeoArray(filename;
        dims=dims(dst), refdims=refdims(dst), name=name(dst),
        metadata=metadata(dst), missingval=missingval(dst),
    )
    read!(src, dst)
end
function Base.read!(filenames::Union{NamedTuple,<:AbstractVector{<:AbstractString}}, dst::AbstractGeoStack)
    _readstack!(filenames, dst)
end
function Base.read!(filenames::AbstractString, dst::AbstractGeoStack)
    _readstack!(filenames, dst)
end
function Base.read!(filenames::AbstractVector{<:Union{AbstractString,NamedTuple}}, dst::AbstractGeoSeries)
    map((fn, d) -> read!(fn, d), filenames, dst)
    return dst
end

function _readstack!(filenames, dst)
    src = GeoStack(filenames;
        dims=dims(dst), refdims=refdims(dst), keys=keys(dst), metadata=metadata(dst),
        layermetadata=DD.layermetadata(dst), layermissingval=layermissingval(dst),
    )
    read!(src, dst)
end


text:
using DataDeps

register(DataDep(
  "soybean-large",
  "http://archive.ics.uci.edu/ml/datasets/Soybean+%28Large%29",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/soybean/soybean-large.data",
  "04b99f2728ded9f544022d2b4f6cce0ebe11ef5c9d44f10e21acd7093876507e",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:36,
        categoric_indices=2:36,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
@diffusion_process SIR begin
    :parameters
    (α, β, σ1, σ2) --> Float64

    :conjugate
    phi(t, u) --> (
        (0.0, 0.0),
        ((k - u[1] - u[2])*u[1], 0.0),
        (-u[1], u[1]),
        (0.0, 0.0),
        (0.0, 0.0),
        (0.0, 0.0),
    )
    nonhypo(x) --> x
    num_non_hypo --> 2

    :additional
    statespace -->BoundedStateSpace(((1,2), (0.0,0.0)), ((1,2), (1.0,1.0)))
end

function b(t, x, P::SIR)
    ℝ{2}(
        P.α*(1 - u[1] - u[2])*u[1] - P.β*u[1],
        P.β*u[1]
    )
end

function σ(t, x, P::SIR)
    @SMatrix Float64[
    -P.σ1*sq((1 - u[1] - u[2])*u[1])  -P.σ2*sq(u[1]);
    0.0   P.σ2*sq.(u[1])
]


text:
export AbstractTurnBuffer, SARD, SARDSA, SARDBuffer, SARDSABuffer,
       isfull, capacity, buffers, empty!

import Base: size, getindex, setindex!, length, push!, eltype, view,
             empty!, isempty, getproperty, lastindex

import DataStructures:isfull, capacity

"""
    AbstractTurnBuffer{names, types} <: AbstractArray{NamedTuple{names, types}, 1}

`AbstractTurnBuffer` is supertype of a collection of buffers to store the interactions between agents and environments.
It is a subtype of `AbstractArray{NamedTuple{names, types}, 1}` where `names` specifies which fields are to store
and `types` is the coresponding types of the `names`.

| Required Methods| Brief Description |
|:----------------|:------------------|
| `Base.push!(b::AbstractTurnBuffer{names, types}, s[, a, r, d, s′, a′])` | Push a turn info into the buffer. According to different `names` and `types` of the buffer `b`, it may accept different number of arguments |
| `isfull(b)` | Check whether the buffer is full or not |
| `capacity(b)` | The maximum length of buffer |
| `Base.length(b)` | Return the length of buffer |
| `Base.getindex(b::AbstractTurnBuffer{names, types})` | Return a turn of type `NamedTuple{names, types}` |
| `Base.empty!(b)` | Reset the buffer |
| **Optional Methods** | |
| `Base.size(b)` | Return `(length(b),)` by default |
| `Base.isempty(b)` | Check whether the buffer is empty or not. Return `length(b) == 0` by default |
| `Base.lastindex(b)` | Return `length(b)` by default |
"""
abstract type AbstractTurnBuffer{names, types} <: AbstractArray{NamedTuple{names, types}, 1} end

const SARD = (:state, :action, :reward, :isdone)
const SARDS = (:state, :action, :reward, :isdone, :nextstate)
const SARDSA = (:state, :action, :reward, :isdone, :nextstate, :nextaction)
const SARDBuffer = AbstractTurnBuffer{SARD}
const SARDSBuffer = AbstractTurnBuffer{SARDS}
const SARDSABuffer = AbstractTurnBuffer{SARDSA}

function isfull end
function capacity end
buffers(b::AbstractTurnBuffer) = getfield(b, :buffers)

size(b::AbstractTurnBuffer) = (length(b),)
isempty(b::AbstractTurnBuffer) = length(b) == 0
lastindex(b::AbstractTurnBuffer) = length(b)
getindex(b::AbstractTurnBuffer, i::Int) = eltype(b)(x[i] for x in buffers(b))
empty!(b::AbstractTurnBuffer) = for x in buffers(b) empty!(x) end
getproperty(b::AbstractTurnBuffer, f::Symbol) = getfield(buffers(b), f)


text:
import GZip
using StatsBase
using IterTools
using LightGraphs
using Formatting
#derrida_map calculates derrida map for a net up to a  H hamming distance
function derrida_map!(netstate::Vector{BNode},nrand_init::Integer=prod([length(sigma.state_range) for sigma in netstate]);tau::Integer=1,Hmax::Integer=sum([length(sigma.state_range)>1 for sigma in netstate]))
    index_nodes_to_flip = find(sigma->length(sigma.state_range)>1,netstate)
    N = length(index_nodes_to_flip)
    1<=Hmax<=N || error("Hmax must be less or equal than network size")
    Hmap = Tuple{Float64,Float64}[]
    pairs = [Dict{AbstractString,AbstractString}() for i in 1:Hmax]
    size_space_state = prod([length(sigma.state_range) for sigma in netstate])
    if nrand_init>size_space_state
        nrand_init = size_space_state
    end
    # tag_base = maximum([length(sigma.state_range) for sigma in netstate])
    nredundant = 0
    Xa = Vector{Integer}(length(netstate))
    Xb = Vector{Integer}(length(netstate))
    i = 0
    while i < nrand_init
        Xa = rand_init_net!(netstate)
        Xa_tag = join(map(string,round(Int,Xa)))
        for H in 1:Hmax
            already_visitted = true
            while already_visitted
                Xb = deepcopy(Xa)
                for j in sample(index_nodes_to_flip, H, replace = false, ordered = true)
                    Xb[j] = rand(filter(x -> x!=Xa[j],netstate[j].state_range))
                end
                Xb_tag = join(map(string,round(Int,Xb)))
                if in((Xa_tag => Xb_tag),pairs[H]) #|| in((Xb_tag => Xa_tag),pairs[H])
                    nredundant += 1
                else
                    pairs[H][Xa_tag] = Xb_tag
                    already_visitted = false
                end
            end
    		for t in 1:tau
                for k in eachindex(netstate)
                    init_node!(netstate[k],Xa[k])
                end
    			Xa = [evolve_node(netstate,nsigma) for nsigma in 1:length(netstate)]
    		end
            for t in 1:tau
                for k in eachindex(netstate)
                    init_node!(netstate[k],Xb[k])
                end
    			Xb = [evolve_node(netstate,nsigma) for nsigma in 1:length(netstate)]
    		end
            push!(Hmap,(H/N,sum(Xa.!=Xb)/N))
            if mod(i,10000)==0
                @printf("%0.2f%% progress, %d out of %d random initial conditions explored at distance H(0)=%d\n", 100*i/nrand_init, i, nrand_init,H)
            end
        end
        i += 1
    end
    println("$nredundant revisitted configuration pairs were discarded")
    return reshape(Hmap,Hmax,nrand_init)
end

function derrida_map2!(netstate::Vector{BNode},nrand_init::Integer=prod([length(sigma.state_range) for sigma in netstate]);tau::Integer=1,Hmax::Integer=sum([length(sigma.state_range)>1 for sigma in netstate]))
    index_nodes_to_flip = find(sigma->length(sigma.state_range)>1,netstate)
    N = length(index_nodes_to_flip)
    1<=Hmax<=N || error("Hmax must be less or equal than network size")
    Hmap = Tuple{Float64,Float64}[]
    pairs = [Dict{AbstractString,AbstractString}() for i in 1:Hmax]
    # tag_base = maximum([length(sigma.state_range) for sigma in netstate])
    nredundant = 0
    Xa = Vector{Integer}(length(netstate))
    Xb = Vector{Integer}(length(netstate))
    Xatp1 = Vector{Integer}(length(netstate))
    Xbtp1 = Vector{Integer}(length(netstate))
    range_lengths = [length(sigma.state_range) for sigma in netstate]
    state_space_size = prod(range_lengths)
    if nrand_init>state_space_size
        nrand_init = state_space_size
        println("Oversampling error: Random initial conditions bigger than possible number of configurations")
    end
    state_ranges = [shuffle(sigma.state_range) for sigma in netstate]
    state_space = product(state_ranges...)
    i = 0
    for Xa in cycle(state_space)
        if rand()>1/length(netstate[end].state_range)
            continue
        end
        for (k,sigma) in enumerate(netstate)
            init_node!(sigma,Xa[k])
        end
        Xa_tag = join(map(string,round(Int,[Xa...])))
        if haskey(pairs[1],Xa_tag)
            continue
        end
        for t in 1:tau
            Xatp1 = [evolve_node(netstate,nsigma) for nsigma in 1:length(netstate)]
            for k in eachindex(netstate)
                init_node!(netstate[k],Xatp1[k])
            end
        end
        for H in 1:Hmax
            Xb = deepcopy([Xa...])
            for j in shuffle(index_nodes_to_flip)[1:H]
                Xb[j] = rand(filter(x -> x!=Xa[j],netstate[j].state_range))
            end
            Xb_tag = join(map(string,round(Int,Xb)))
            if in((Xa_tag => Xb_tag),pairs[H]) || in((Xb_tag => Xa_tag),pairs[H])
                nredundant += 1
            end
            pairs[H][Xa_tag] = Xb_tag
            for k in eachindex(netstate)
                init_node!(netstate[k],Xb[k])
            end
            for t in 1:tau
                Xbtp1 = [evolve_node(netstate,nsigma) for nsigma in 1:length(netstate)]
                for k in eachindex(netstate)
                    init_node!(netstate[k],Xbtp1[k])
                end
            end
            push!(Hmap,(H/N,sum(Xatp1.!=Xbtp1)/N))
            if mod(i,10000)==0
                @printf("%0.2f%% progress, %d out of %d random initial conditions explored at distance H(0)=%d\n", 100*i/nrand_init, i, nrand_init,H)
            end
        end
        i += 1
        if i == nrand_init
            break
        end
    end
    println("$nredundant/$(Hmax*nrand_init) configuration pairs redundant")
    return reshape(Hmap,Hmax,nrand_init)
end

function derrida_map3!(netstate::Vector{BNode},nrand_init::Integer=prod([length(sigma.state_range) for sigma in netstate]);tau::Integer=1,Hmax::Integer=sum([length(sigma.state_range)>1 for sigma in netstate]))
    index_nodes_to_flip = find(sigma->length(sigma.state_range)>1,netstate)
    N = length(index_nodes_to_flip)
    1<=Hmax<=N || error("Hmax must be less or equal than network size")
    Hmap = Tuple{Float64,Float64}[]
    pairs = [Dict{AbstractString,AbstractString}() for i in 1:Hmax]
    # tag_base = maximum([length(sigma.state_range) for sigma in netstate])
    nredundant = 0
    Xa = Vector{Integer}(length(netstate))
    Xb = Vector{Integer}(length(netstate))
    Xatp1 = Vector{Integer}(length(netstate))
    Xbtp1 = Vector{Integer}(length(netstate))
    range_lengths = [length(sigma.state_range) for sigma in netstate]
    state_space_size = prod(range_lengths)
    if nrand_init>state_space_size
        nrand_init = state_space_size
        println("Oversampling error: Random initial conditions bigger than possible number of configurations")
    end
    pat_div = reverse(round(Int,cumprod(reverse(range_lengths))./reverse(range_lengths)))
    state_space_number = shuffle(collect(1:state_space_size))
    i = 0
    for conf_i in state_space_number
        conf_indices = mod(div((conf_i-1),pat_div),range_lengths)
        Xa = [sigma.state_range[conf_indices[k]+1] for (k,sigma) in enumerate(netstate)]
        Xa_tag = join(map(string,round(Int,Xa)))
        for k in eachindex(netstate)
            init_node!(netstate[k],Xa[k])
        end
        for t in 1:tau
            Xatp1 = [evolve_node(netstate,nsigma) for nsigma in 1:length(netstate)]
            for k in eachindex(netstate)
                init_node!(netstate[k],Xatp1[k])
            end
        end
        for H in 1:Hmax
            Xb = deepcopy(Xa)
            for j in shuffle(index_nodes_to_flip)[1:H]
                Xb[j] = rand(filter(x -> x!=Xa[j],netstate[j].state_range))
            end
            Xb_tag = join(map(string,round(Int,Xb)))
            if in((Xa_tag => Xb_tag),pairs[H]) || in((Xb_tag => Xa_tag),pairs[H])
                nredundant += 1
            end
            pairs[H][Xa_tag] = Xb_tag
            for k in eachindex(netstate)
                init_node!(netstate[k],Xb[k])
            end
            for t in 1:tau
                Xbtp1 = [evolve_node(netstate,nsigma) for nsigma in 1:length(netstate)]
                for k in eachindex(netstate)
                    init_node!(netstate[k],Xbtp1[k])
                end
            end
            push!(Hmap,(H/N,sum(Xatp1.!=Xbtp1)/N))
            if mod(i,10000)==0
                @printf("%0.2f%% progress, %d out of %d random initial conditions explored at distance H(0)=%d\n", 100*i/nrand_init, i, nrand_init,H)
            end
        end
        i += 1
        if i == nrand_init
            break
        end
    end
    println("$nredundant/$(Hmax*nrand_init) configuration pairs redundant")
    return reshape(Hmap,Hmax,nrand_init)
end

function derrida_map_att3!(netstate::Vector{BNode},state_space_tags::Set{Int},nrand_init::Integer=length(state_space_tags);tau::Integer=1,Hmax::Integer=sum([length(sigma.state_range)>1 for sigma in netstate]))
    index_nodes_to_flip = find(sigma->length(sigma.state_range)>1,netstate)
    N = length(index_nodes_to_flip)
    net_size = length(netstate)
    1<=Hmax<=N || error("Hmax must be less or equal than network size")
    Hmap = Tuple{Float64,Float64}[]
    pairs = [Dict{AbstractString,AbstractString}() for i in 1:Hmax]
    tag_base = maximum([length(sigma.state_range) for sigma in netstate])
    nredundant = 0
    Xa = Vector{Integer}(length(netstate))
    Xb = Vector{Integer}(length(netstate))
    Xatp1 = Vector{Integer}(length(netstate))
    Xbtp1 = Vector{Integer}(length(netstate))
    sampled_confs = shuffle(collect(state_space_tags))
    if nrand_init>length(state_space_tags)
        nrand_init = length(state_space_tags)
    end
    i = 0
    for conf_i in sampled_confs
        Xa = reverse(digits(conf_i,tag_base,net_size))
        Xa_tag = join(map(string,round(Int,Xa)))
        for k in eachindex(netstate)
            init_node!(netstate[k],Xa[k])
        end
        Xatp1 = evolve_net!(netstate,tau,keep=false)
        for H in 1:Hmax
            belongs_to_att = false
            already_visitted = true
            while !belongs_to_att && already_visitted
                Xb = deepcopy(Xa)
                for j in sample(index_nodes_to_flip, H, replace = false, ordered = true)
                    Xb[j] = rand(filter(x -> x!=Xa[j],netstate[j].state_range))
                end
                Xb_tag = parse(Int,join(map(string,Xa)),tag_base)
                if in(Xb_tag,state_space_tags)
                    belongs_to_att = true
                end
                if in((Xa_tag => Xb_tag),pairs[H]) || in((Xb_tag => Xa_tag),pairs[H])
                    nredundant += 1
                    belongs_to_att = false
                else
                    pairs[H][Xa_tag] = Xb_tag
                    already_visitted = false
                    for k in eachindex(netstate)
                        init_node!(netstate[k],Xb[k])
                    end
                    Xbtp1 = evolve_net!(netstate,tau,keep=false)
                    push!(Hmap,(H/N,sum(Xatp1.!=Xbtp1)/N))
                    push!(Hmap,(H/N,sum(Xatp1.!=Xbtp1)/N))
                    if mod(i,10000)==0
                        @printf("%0.2f%% progress, %d out of %d random initial conditions explored at distance H(0)=%d\n", 100*i/nrand_init, i, nrand_init,H)
                    end
                end
            end
        end
        i += 1
        if i == nrand_init
            break
        end
    end
    println("$nredundant revisitted configuration pairs were discarded")
    return reshape(Hmap,Hmax.nrand_init)
end

function full_derrida_map!(netstate::Vector{BNode},nrand_init::Integer;tau::Integer=1,Hmax::Integer=sum([length(sigma.state_range)>1 for sigma in netstate]))
    nconf = prod([length(sigma.state_range) for sigma in netstate])
    index_nodes_to_flip = find(sigma->length(sigma.state_range)>1,netstate)
    N = length(index_nodes_to_flip)
    1<=Hmax<=N || error("Hmax must be less or equal than network size")
    Hmap = Dict(i=>Tuple{AbstractFloat,AbstractFloat}[] for i in 1:Hmax)
    nredundant = 0
    Xatp1 = Vector{Integer}(length(netstate))
    Xbtp1 = Vector{Integer}(length(netstate))
    state_ranges = [shuffle(sigma.state_range) for sigma in netstate]
    state_space = collect(product(state_ranges...))
    i = 0
    for Xa in product([shuffle(sigma.state_range) for sigma in netstate]...)
        for (k,sigma) in enumerate(netstate)
            init_node!(sigma,Xa[k])
        end
        for t in 1:tau
            Xatp1 = [evolve_node(netstate,nsigma) for nsigma in 1:length(netstate)]
            for k in eachindex(netstate)
                init_node!(netstate[k],Xatp1[k])
            end
        end
        for Xb in product([shuffle(sigma.state_range) for sigma in netstate]...)
            if Xa==Xb
                continue
            end
            for (k,sigma) in enumerate(netstate)
                init_node!(sigma,Xb[k])
            end
            for t in 1:tau
                Xbtp1 = [evolve_node(netstate,nsigma) for nsigma in 1:length(netstate)]
                for k in eachindex(netstate)
                    init_node!(netstate[k],Xbtp1[k])
                end
            end
            H = sum(Xa.!=Xb)
            push!(Hmap[H],(H/N,sum(Xatp1.!=Xbtp1)/N))
            i+=1
            if mod(i,10000)==0
                @printf("%0.2f%% progress. H(0)=%d,%d\n", 100*i/nrand_init,H,length(Hmap[H]))
            end
            if i == nrand_init
                break
            end
        end
    end
    for H in sort(collect(keys(Hmap)))
        @show H,length(Hmap[H])
    end
    return Hmap
end

function derrida_map_mod!(netWTstate::Vector{BNode},netMutantstate::Vector{BNode},nrand_init::Integer=prod([length(sigma.state_range) for sigma in netMutantstate]);tau::Integer=1,fix_targets::Vector{Symbol}=Symbol[],Hmax::Integer=sum([length(sigma.state_range)>1 for sigma in netWTstate]))
    index_nodes_to_fix = Integer[]
    for id in fix_targets
        i = findfirst(x->x.id==string(id),netWTstate)
        j = findfirst(x->x.id==string(id),netMutantstate)
        (i != 0 && j != 0) || error("Trying to set unexisting $id node")
        i == j || error("Comparative derrida map requires networks to have the same size and same node order")
        push!(index_nodes_to_fix,i)
	end
    Hmin = length(index_nodes_to_fix)
    index_nodes_to_flip = find(x->length(x.state_range)>1, netWTstate)
    N = length(netWTstate)
    Hmin<=Hmax<=N || error("Hmax must be greater or equal than size of fixed targets, and equal or less than network size")
    deleteat!(index_nodes_to_flip,findin(index_nodes_to_flip,index_nodes_to_fix))
    Hmap = Tuple{Float64,Float64}[]
    size_Mutant_space_state = prod([length(sigma.state_range) for sigma in netMutantstate])
    if nrand_init>size_Mutant_space_state
        nrand_init = size_Mutant_space_state
    end
    pairs = Dict(i=>Dict{AbstractString,AbstractString}() for i in Hmin:Hmax)
    nredundant = 0
    XWTtp1 = Vector{Integer}(length(netWTstate))
    XMTtp1 = Vector{Integer}(length(netMutantstate))
    XWT = Vector{Integer}(length(netWTstate))
    XMT = Vector{Integer}(length(netMutantstate))
    i = 0
    while i < nrand_init
        XMT = rand_init_net!(netMutantstate)
        XMT_tag = join(map(string,round(Int,XMT)))
        for j in eachindex(netWTstate)
            init_node!(netWTstate[j],netMutantstate[j].state)
        end # the KO configuration is copied into the wild type (unperturbed) network
        for j in index_nodes_to_fix
            init_node!(netWTstate[j],rand(filter(x -> x!=netMutantstate[j].state,netWTstate[j].state_range)))
        end # those nodes to be shut off in the knock out network are set in the wild type network to any other state different than that
        tmpXWT = [sigma.state for sigma in netWTstate]
        for H in Hmin:Hmax
            already_visitted = true
            while already_visitted
                for j in eachindex(tmpXWT)
                    init_node!(netWTstate[j],tmpXWT[j])
                end # the KO configuration is copied into the wild type (unperturbed) network
                for j in sample(index_nodes_to_flip, H - Hmin, replace = false, ordered = true)
                    init_node!(netWTstate[j],rand(filter(x -> x!=netMutantstate[j].state,netWTstate[j].state_range)))
                end # the rest of the wildtype nodes are randomly set, up to H changes are done
                XWT = [sigma.state for sigma in netWTstate]
                XWT_tag = join(map(string,round(Int,XWT)))
                if (XMT_tag != XMT_tag) & (in((XWT_tag => XMT_tag),pairs[H]) || in((XMT_tag => XWT_tag),pairs[H]))
                    nredundant += 1
                else
                    pairs[H][XMT_tag] = XWT_tag
                    already_visitted = false
                end
            end
            for t in 1:tau
                XWTtp1 = [evolve_node(netWTstate,nsigma) for nsigma in 1:length(netWTstate)]
                for k in eachindex(netWTstate)
                    init_node!(netWTstate[k],XWTtp1[k])
                end
            end
            for t in 1:tau
                XMTtp1 = [evolve_node(netMutantstate,nsigma) for nsigma in 1:length(netMutantstate)]
                for k in eachindex(netMutantstate)
                    init_node!(netMutantstate[k],XMTtp1[k])
                end
            end
            # push!(Hmap,(H/N,sum(XWTtp1.!=XMTtp1)/N))
            push!(Hmap,(sum(XWT.!=XMT)/N,sum(XWTtp1.!=XMTtp1)/N))
            if mod(i,10000)==0
                @printf("%0.2f%% progress, %d out of %d random initial conditions explored at distance H(0)=%d\n", 100*i/nrand_init, i, nrand_init,H)
            end
        end
        i += 1
    end
    println("$nredundant revisitted configuration pairs were discarded")
    return reshape(Hmap,Hmax-Hmin+1,nrand_init)
end

function derrida_map_mod3!(netWTstate::Vector{BNode},netMutantstate::Vector{BNode},nrand_init::Integer=prod([length(sigma.state_range) for sigma in netMutantstate]);tau::Integer=1,fix_targets::Vector{Symbol}=Symbol[],Hmax::Integer=sum([length(sigma.state_range)>1 for sigma in netWTstate]))
    nconfMT = prod([length(sigma.state_range) for sigma in netMutantstate])
    if nrand_init>nconfMT
        nrand_init = nconfMT
        println("Oversampling error: Random initial conditions bigger than possible number of configurations")
    end
    index_nodes_to_fix = Integer[]
    for id in fix_targets
        i = findfirst(x->x.id==string(id),netWTstate)
        j = findfirst(x->x.id==string(id),netMutantstate)
        (i != 0 && j != 0) || error("Trying to set unexisting $id node")
        i == j || error("Comparative derrida map requires networks to have the same size and same node order")
        push!(index_nodes_to_fix,i)
    end
    Hmin = length(index_nodes_to_fix)
    index_nodes_to_flip = find(x->length(x.state_range)>1, netWTstate)
    N = length(netWTstate)
    Hmin<=Hmax<=N || error("Hmax must be greater or equal than size of fixed targets, and equal or less than network size")
    deleteat!(index_nodes_to_flip,findin(index_nodes_to_flip,index_nodes_to_fix))
    Hmap = Tuple{Float64,Float64}[]
    pairs = Dict(i=>Dict{AbstractString,AbstractString}() for i in Hmin:Hmax)
    nredundant = 0
    XWTtp1 = Vector{Integer}(length(netWTstate))
    XMTtp1 = Vector{Integer}(length(netMutantstate))
    XWT = Vector{Integer}(length(netWTstate))
    XMT = Vector{Integer}(length(netMutantstate))
    range_lengths = [length(sigma.state_range) for sigma in netMutantstate]
    pat_div = reverse(round(Int,cumprod(reverse(range_lengths))./reverse(range_lengths)))
    state_space_number = shuffle(collect(1:prod(range_lengths)))
    i = 0
    for conf_i in state_space_number
        conf_indices = mod(div((conf_i-1),pat_div),range_lengths)
        XMT = [sigma.state_range[conf_indices[k]+1] for (k,sigma) in enumerate(netMutantstate)]
        XMT_tag = join(map(string,round(Int,XMT)))
        for j in eachindex(netMutantstate)
            init_node!(netWTstate[j],XMT[j])
        end # the KO configuration is copied into the wild type (unperturbed) network
        for j in index_nodes_to_fix
            init_node!(netWTstate[j],rand(filter(x -> x!=XMT[j],netWTstate[j].state_range)))
        end # those nodes to be shut off in the knock out network are set in the wild  type network to any other state different than that
        tmpXWT = [sigma.state for sigma in netWTstate]
        for H in Hmin:Hmax
            already_visitted = true
            for j in eachindex(netMutantstate)
                init_node!(netMutantstate[j],XMT[j])
            end # the KO configuration is copied into the wild type (unperturbed) network
            while already_visitted
                for j in eachindex(tmpXWT)
                    init_node!(netWTstate[j],tmpXWT[j])
                end # the KO configuration is copied into the wild type (unperturbed) network
                for j in sample(index_nodes_to_flip, H - Hmin, replace = false, ordered = true)
                    init_node!(netWTstate[j],rand(filter(x -> x!=XMT[j],netWTstate[j].state_range)))
                end # the rest of the wildtype nodes are randomly set, up to H changes are done
                XWT = [sigma.state for sigma in netWTstate]
                XWT_tag = join(map(string,round(Int,XWT)))
                if (XMT_tag != XMT_tag) & (in((XWT_tag => XMT_tag),pairs[H]) || in((XMT_tag => XWT_tag),pairs[H]))
                    nredundant += 1
                else
                    pairs[H][XMT_tag] = XWT_tag
                    already_visitted = false
                end
            end
            for t in 1:tau
                XWTtp1 = [evolve_node(netWTstate,nsigma) for nsigma in 1:length(netWTstate)]
                for k in eachindex(netWTstate)
                    init_node!(netWTstate[k],XWTtp1[k])
                end
            end
            for t in 1:tau
                XMTtp1 = [evolve_node(netMutantstate,nsigma) for nsigma in 1:length(netMutantstate)]
                for k in eachindex(netMutantstate)
                    init_node!(netMutantstate[k],XMTtp1[k])
                end
            end
            # push!(Hmap,(H/N,sum(XWTtp1.!=XMTtp1)/N))
            push!(Hmap,(sum(XWT.!=XMT)/N,sum(XWTtp1.!=XMTtp1)/N))
            if mod(i,10000)==0
                @printf("%0.2f%% progress, %d out of %d random initial conditions explored at distance H(0)=%d\n", 100*i/nrand_init, i, nrand_init,H)
            end
        end
        i += 1
        if i==nrand_init
            break
        end
    end
    println("$nredundant revisitted configuration pairs were discarded")
    return reshape(Hmap,Hmax-Hmin+1,nrand_init)
end

# get_basins explores all possible configurations and returns a transition graph, with a dictionary storing labels for
# each configuration

function get_basins!(netstate::Vector{BNode},out_dir::AbstractString,root_name::AbstractString;trav::Symbol=:rw,update_mode::Symbol=:sy,max_traj::Integer=maximum(vcat([sigma.tau_in for sigma in netstate]...)))
    trav == :rw || trav == :l || error("Invalid option for traversal method")
    net_size = length(netstate)
    state_ranges = [sigma.state_range for sigma in netstate]
    itstspace = product(state_ranges...)
    size_space_state = length(itstspace)
    tags_nid = update_mode == :sy ? Dict{Integer,Integer}() : Dict{Vector{Integer},Integer}()
    nid_tags = update_mode == :sy ? Dict{Integer,Integer}() : Dict{Integer,Vector{Integer}}()
    tag_base = maximum([length(sigma.state_range) for sigma in netstate])
    transition_graph = DiGraph()
    println("Calculating attraction landscape for $root_name")
    tau_max = maximum(vcat([sigma.tau_in for sigma in netstate]...))
    net_hist_size = tau_max + 1
    if update_mode == :sy
        for leave in itstspace
            for k in eachindex(netstate)
                init_node!(netstate[k],leave[k])
            end
            leave_tag = parse(Int,string(map(x->round(Int,x),leave)...),tag_base)
            shoot_tag = parse(Int,string(round(Int,shoot)...),tag_base)
            if !haskey(tags_nid,leave_tag)
                add_vertex!(transition_graph)
                nvertex = nv(transition_graph)
                nid_tags[nvertex] = leave_tag
                tags_nid[leave_tag] = nvertex
            end
            shoot = evolve_net!(netstate,1,keep=false,update_mode=update_mode)
            if !haskey(tags_nid,shoot_tag)
                add_vertex!(transition_graph)
                nvertex = nv(transition_graph)
                nid_tags[nvertex] = shoot_tag
                tags_nid[shoot_tag] = nvertex
            end
            add_edge!(transition_graph,tags_nid[leave_tag],tags_nid[shoot_tag])
            if mod(nv(transition_graph),100000)==0
                @printf "%0.2f%% configurations explored\n" 100*nv(transition_graph)/size_space_state
            end
        end
    else
        hX = Array{Int8}(length(netstate),net_hist_size)
        hX_tag = Vector{Int}(net_hist_size)
        hXtp1 = Array{Int8}(length(netstate),net_hist_size)
        hXtp1_tag = Vector{Int}(net_hist_size)
        time_steps = 1
        sim_buffer = Array{Int8}(length(netstate),time_steps + net_hist_size)
        for (i,X) in enumerate(itstspace)
            hX[:] = init_net!(netstate;init_t0=:default,map((x,y)->Symbol(x)=>y,(sigma.id for sigma in netstate),X)...)
            # for (nsigma,sigma) in enumerate(netstate)
            #     init_node!(sigma,X[nsigma])
            #     hX[nsigma,:] = sigma.state
            # end
            for t in 1:net_hist_size
                hX_tag[t] = parse(Int,string(round(Int,hX[:,t])...),tag_base)
            end
            if !haskey(tags_nid,hX_tag)
                add_vertex!(transition_graph)
                nvertex = nv(transition_graph)
                nid_tags[nvertex] = hX_tag
                tags_nid[hX_tag] = nvertex
            end
            # hXtp1_tag[:] = hX_tag[:]
            # shift!(hXtp1_tag)
            for t in 0:max_traj-1
                hXtp1[:] = evolve_net!(netstate,1,keep=false,update_mode=update_mode,sim=sim_buffer,t0=t)
                for ti in 1:net_hist_size
                    hXtp1_tag[ti] = parse(Int,string(round(Int,hXtp1[:,ti])...),tag_base)
                end
                # Xtp1_tag = parse(Int,string(round(Int,Xtp1)...),tag_base)
                # push!(hXtp1_tag,Xtp1_tag)
                if !haskey(tags_nid,hXtp1_tag)
                    add_vertex!(transition_graph)
                    nvertex = nv(transition_graph)
                    nid_tags[nvertex] = hXtp1_tag
                    tags_nid[hXtp1_tag] = nvertex
                end
                add_edge!(transition_graph,tags_nid[hX_tag],tags_nid[hXtp1_tag]) || break
                hX_tag[:] = hXtp1_tag[:]
                # shift!(hXtp1_tag)
            end
            if mod(i,100000)==0
                @printf "%0.2f%% initial conditions explored\n" 100*i/size_space_state
            end
        end
    end
    size_space_state = nv(transition_graph)
    println("Decomposing transition graph into attraction basins")
    basin_node_list = weakly_connected_components(transition_graph)
    Natt = length(basin_node_list)
    println(Natt," attractors found")
    sort!(basin_node_list,by=length,rev=true)
    println("Saving individual basins")
    if !isdir(out_dir)
        mkpath(out_dir)
    end
    g = GZip.open(joinpath(out_dir,"$(root_name)_att_tagging.dat.gz"),"w")
    h = GZip.open(joinpath(out_dir,"$(root_name)_trans_dist.dat.gz"),"w")
    try
        basin_tree_list = Vector{LightGraphs.DiGraph}(Natt)
        basin_tags = Vector{Vector{Integer}}(Natt)
        att_list = Vector{Vector{Integer}}(Natt)
        periods = Vector{Integer}(Natt)
        branch_sizes = Vector{Vector{Integer}}(Natt)
        patt_basins = Vector{Vector{Integer}}(Natt)
        branch_tree_list = Vector{Vector{LightGraphs.DiGraph}}(Natt)
        branch_patt = Vector{Vector{Integer}}(Natt)
        for (bi,basin) in enumerate(basin_node_list)
            basin_tree_list[bi] = induced_subgraph(transition_graph,basin)[1]
            save(joinpath(out_dir,"$(root_name)_att_$(bi).gml.gz"),basin_tree_list[bi],join([root_name, bi],"_"),:gml,compress = true)
            println("$(root_name)_att_$(bi) basin saved")
            if update_mode == :sy
                writedlm(g,hcat(fill(bi,nv(basin_tree_list[bi])),[nid_tags[v] for v in basin]))
            else
                writedlm(g,hcat(fill(bi,nv(basin_tree_list[bi])),reshape([nid_tags[v][d] for d in 1:net_hist_size for v in basin],length(basin),net_hist_size)))
            end
            println("Computing attractor of basin $bi")
            probe = basin[1]
            path_to_att = saw(transition_graph,probe,length(basin))
            path_in_att = saw(transition_graph,path_to_att[end],length(basin))
            att_list[bi] = path_in_att
            periods[bi] = length(path_in_att)
            println("Attractor $bi has period ",periods[bi])
            patt_basin = Vector{Integer}()
            for point_in_att in path_in_att
                push!(patt_basin, findfirst(basin,point_in_att))
            end
            patt_basins[bi] = patt_basin
        end
        transition_graph = DiGraph()
        for (bi,basin_tree) in enumerate(basin_tree_list)
            println("Calculating branches associated to attractor $bi")
            edges_att = hcat(circshift(patt_basins[bi],1),patt_basins[bi])
            for pi in 1:length(patt_basins[bi])
                rem_edge!(basin_tree,edges_att[pi,:]...)
            end
            branch_node_list = weakly_connected_components(basin_tree)
            branch_tree_list[bi] = Vector{LightGraphs.DiGraph}(periods[bi])
            patt_to_look = deepcopy(patt_basins[bi])
            branch_patt[bi] = Vector{Integer}(periods[bi])
            for (bri,branch) in enumerate(branch_node_list)
                for point_in_att in patt_to_look
                    if point_in_att in branch
                        branch_patt[bi][bri] = findfirst(patt_basins[bi],point_in_att)
                        branch_tree_list[bi][bri] = induced_subgraph(basin_tree,branch)[1]
                        # println("Branch $bi:$bri has $(nv(branch_tree_list[bi][bri])) nodes")
                        deleteat!(patt_to_look,findfirst(patt_to_look,point_in_att))
                        break
                    end
                end
            end
        end
        basin_tree_list = Vector{LightGraphs.DiGraph}()
        for (bi,basin_tree) in enumerate(branch_tree_list)
            println("Computing transients distribution associated to attractor $bi")
            dists = Integer[]
            brsize = Vector{Integer}(periods[bi])
            for (bri,branch_tree) in enumerate(basin_tree)
                println("Measuring transient distances for branch $bi:$bri")
                dists_br = Integer[]
                if (trav==:rw)
                    maxbr_size = nv(branch_tree)
                    for pbi in vertices(branch_tree)
                      # push!(dists_br,length(a_star(basin_tree,pbranch,branch_patt[bri])))
                        push!(dists_br,length(non_backtracking_randomwalk(branch_tree,pbi,maxbr_size)-1))
                        if mod(pbi,100000)==0
                            @printf("%0.2f%% of branch %d:%d explored\n",100*pbi/maxbr_size, bri, bi)
                        end
                    end
                    filter!(x->x!=0,dists_br)
                else
                    m = mapreduce(x->in_neighbors(branch_tree,x),vcat,vertices(branch_tree))
                    d = sparsevec(m,ones(Integer,length(m)),nv(branch_tree))
                    l = 1
                    println("Level $l")
                    while !isempty(m)
                        m = mapreduce(x->in_neighbors(branch_tree,x),vcat,m)
                        d[m] += 1
                        l+=1
                        println("Level $l")
                    end
                    dists_br = nonzeros(d)
                end
                append!(dists,dists_br)
                brsize[branch_patt[bi][bri]] = length(dists_br)
            end
            branch_sizes[bi] = brsize
            writedlm(h,dists',',')
            println("Maximum transient distance:$(maximum(dists))")
        end
    finally
        close(g)
        close(h)
    end
    f = open(joinpath(out_dir,"$(root_name)_att_stats.txt"),"w")
    sum_att = sum(periods)
    sum_basin = size_space_state - sum_att
    try
        println("Writing statistics")
        for bi in eachindex(att_list)
            if bi > 1
                write(f,"\n")
            end
            write(f,"Attractor:$bi\t")
            period = periods[bi]
            write(f,"Period:$period\t")
            basin_size = mapreduce(x->nv(x),+,branch_tree_list[bi])-period
            write(f,"Basin size:$basin_size/$sum_basin ($(format(basin_size/sum_basin*100.0,precision=2,zeropadding=true))%)\n")
            if update_mode==:sy
                for (pi,conf_id) in enumerate(att_list[bi])
                    conf = reverse(digits(nid_tags[conf_id],tag_base,net_size))
                    write(f,"$(join(conf,","))\t$(branch_sizes[bi][pi])\n")
                end
            else
                for (pi,conf_id) in enumerate(att_list[bi])
                    conf = reverse(digits(nid_tags[conf_id][end],tag_base,net_size))
                    write(f,"$(join(conf,","))\t$(branch_sizes[bi][pi])\n")
                end
            end
        end
    finally
        close(f)
    end
    # return basin_tags,basin_tree_list,att_list
    return transition_graph
end

# get_attractor finds attractors exploring multiple random initial conditions (nrand_init) within a simulation
# time range limited by max_steps, after a transient given by transient_size
function find_attractor!(netstate::Vector{BNode},transient_size::Integer,max_steps::Integer,nrand_init::Integer;update_mode::Symbol=:sy)
    i = 0
    tag_base = maximum([length(sigma.state_range) for sigma in netstate])
    tau_max = update_mode==:sy ? 0 : maximum(vcat([sigma.tau_in for sigma in netstate]...))
    net_hist_size = tau_max + 1
    net_size = length(netstate)
    att_list = Vector{Vector{Tuple{Vararg{Tuple{Vararg{Int8,net_size}},net_hist_size}}}}()
    hX = Array{Int8}(net_size,net_hist_size)
    sim_trans_buffer = Array{Int8}(net_size,transient_size + net_hist_size)
    time_steps = 1
    sim_att_buffer = Array{Int8}(net_size,time_steps + net_hist_size)
    unsuccess = 0
    while i < nrand_init
        init_net!(netstate;init_t0=:rand,init_history=:rand)
        p_att = Vector{Tuple{Vararg{Tuple{Vararg{Int8,net_size}},net_hist_size}}}()
        hX[:] = evolve_net!(netstate,transient_size,keep=false,update_mode=update_mode,sim=sim_trans_buffer)
        probe = tuple(map(x->tuple(x...),(hX[:,ti] for ti in 1:net_hist_size))...)
        push!(p_att,probe)
        add_attractor = false
        already_found = false
        attractor_found = false
        tmax = 1
        for t in transient_size:transient_size+max_steps-1
            hX[:] = evolve_net!(netstate,1,keep=false,update_mode=update_mode,t0=t,sim=sim_att_buffer)
            ntp1 = tuple(map(x->tuple(x...),(hX[:,ti] for ti in 1:net_hist_size))...)
            tmax = t
            if probe[end] == ntp1[end]
                # println("Likely attractor found")
                for att in att_list
                    for p_att_i in p_att
                        if p_att_i in att
                            already_found = true
                            break
                        end
                    end
                    if already_found
                        break
                    end
                end
                attractor_found = true
                if !already_found
                    push!(att_list,p_att)
                    println("Attractor $(size(att_list,1)) has period $(size(p_att,1))")
                    for (j,att_i) in enumerate(p_att)
                        println(join(att_i[end],","))
                    end
                end
                break
            end
            push!(p_att,ntp1)
        end
        i += 1
        if !(already_found||attractor_found)
            println(tmax)
            print("probe:")
            println(join(probe[end],","))
            print("ntp1 :")
            println(join(hX[:,end],","))
            #  println("No attractors after reaching maxtimum time steps")
            unsuccess += 1
            # println("No attractor found for initial condition $i:")
            # for probe_i in probe
            #     conf = reverse(digits(probe_i,tag_base,net_size))
            #     println(join(conf,","))
            # end
        end
    end
    println(unsuccess," unsucessful tries after reaching maximmum time steps set")
    return att_list
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the conversion from Euler angle and axis to DCM.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export angleaxis_to_dcm

"""
    angleaxis_to_dcm(a::Number, v::AbstractVector)
    angleaxis_to_dcm(av::EulerAngleAxis)

Convert the Euler angle `a` [rad] and Euler axis `v` to a DCM.

Those values can also be passed inside the structure `ea` (see
[`EulerAngleAxis`](@ref)).

!!! warning
    It is expected that the vector `v` is unitary. However, no verification is
    performed inside the function. The user must handle this situation.

# Example

```jldoctest
julia> v = [1, 1, 1];

julia> v /= norm(v);

julia> angleaxis_to_dcm(pi / 2, v)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  0.333333   0.910684  -0.244017
 -0.244017   0.333333   0.910684
  0.910684  -0.244017   0.333333

julia> ea = EulerAngleAxis(pi / 2, v);

julia> angleaxis_to_dcm(ea)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  0.333333   0.910684  -0.244017
 -0.244017   0.333333   0.910684
  0.910684  -0.244017   0.333333
```
"""
@inline function angleaxis_to_dcm(a::Number, v::AbstractVector)
    # Check the arguments.
    (length(v) ≠ 3) && throw(ArgumentError("The provided vector for the Euler axis must have 3 elements."))

    sθ, cθ = sincos(a)
    aux    = 1 - cθ

    return DCM(
          cθ + v[1] * v[1] * aux,      v[1] * v[2] * aux + v[3] * sθ, v[1] * v[3] * aux - v[2] * sθ,
        v[1] * v[2] * aux - v[3] * sθ,   cθ + v[2] * v[2] * aux,      v[2] * v[3] * aux + v[1] * sθ,
        v[1] * v[3] * aux + v[2] * sθ, v[2] * v[3] * aux - v[1] * sθ,   cθ + v[3] * v[3] * aux
    )'
end

@inline angleaxis_to_dcm(av::EulerAngleAxis) = angleaxis_to_dcm(av.a, av.v)