text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

module MemoryModule

using LinearAlgebra, ForwardDiff, Tensors, NLsolve, Parameters

import ..AbstractMaterial, ..AbstractMaterialState
import ..Utilities: Symm2, Symm4, isotropic_elasticity_tensor, lame, debang
import ..integrate_material!  # for method extension

# parametrically polymorphic for any type representing ℝ
export GenericMemory, GenericMemoryDriverState, GenericMemoryParameterState, GenericMemoryVariableState

# specialization for Float64
export Memory, MemoryDriverState, MemoryParameterState, MemoryVariableState

@with_kw mutable struct GenericMemoryDriverState{T <: Real} <: AbstractMaterialState
    time::T = zero(T)
    strain::Symm2{T} = zero(Symm2{T})
end

"""Parameter state for Memory material.

- `E`: Young's modulus
- `nu`: Poisson's ratio
- `R0`: initial yield strength
- `Kn`: plasticity multiplier divisor (drag stress)
- `nn`: plasticity multiplier exponent
- `C1`, `D1`: parameters governing behavior of backstress X1
- `C2`, `D2`: parameters governing behavior of backstress X2
- `Q0`: The initial isotropic hardening saturation value. Has the units of stress.
- `QM`: The asymptotic isotropic hardening saturation value reached with high strain amplitude.

  Has the units of stress.

- `mu`: Controls the rate of evolution of the strain-amplitude dependent isotropic hardening saturation value.
- `b`: Controls the rate of evolution for isotropic hardening.
- `eta`: Controls the balance between memory surface kinematic and isotropic hardening.

  Dimensionless, support `[0,1]`.
  At `0`, the memory surface hardens kinematically.
  At `1`, the memory surface hardens isotropically.

  Initially, the value `1/2` was used by several authors. Later, values `< 1/2` have been suggested
  to capture the progressive process of memorization.

- `m`: memory evanescence exponent. Controls the non-linearity of the memory evanescence.
- `pt`: threshold of equivalent plastic strain, after which the memory evanescence starts.
- `xi`: memory evanescence strength multiplier.
"""
@with_kw struct GenericMemoryParameterState{T <: Real} <: AbstractMaterialState
    E::T = 0.0
    nu::T = 0.0
    R0::T = 0.0
    Kn::T = 0.0
    nn::T = 0.0
    C1::T = 0.0
    D1::T = 0.0
    C2::T = 0.0
    D2::T = 0.0
    Q0::T = 0.0
    QM::T = 0.0
    mu::T = 0.0
    b::T = 0.0
    eta::T = 0.0
    m::T = 0.0
    pt::T = 0.0
    xi::T = 0.0
end

"""Problem state for Memory material.

- `stress`: stress tensor
- `X1`: backstress 1
- `X2`: backstress 2
- `plastic_strain`: plastic part of strain tensor
- `cumeq`: cumulative equivalent plastic strain (scalar, ≥ 0)
- `R`: yield strength
- `q`: size of the strain memory surface (~plastic strain amplitude)
- `zeta`: strain memory surface kinematic hardening variable
- `jacobian`: ∂σij/∂εkl
"""
@with_kw struct GenericMemoryVariableState{T <: Real} <: AbstractMaterialState
    stress::Symm2{T} = zero(Symm2{T})
    X1::Symm2{T} = zero(Symm2{T})
    X2::Symm2{T} = zero(Symm2{T})
    plastic_strain::Symm2{T} = zero(Symm2{T})
    cumeq::T = zero(T)
    R::T = zero(T)
    q::T = zero(T)
    zeta::Symm2{T} = zero(Symm2{T})
    jacobian::Symm4{T} = zero(Symm4{T})
end

# TODO: Does this eventually need a {T}?
@with_kw struct MemoryOptions <: AbstractMaterialState
    nlsolve_method::Symbol = :trust_region
end

@with_kw mutable struct GenericMemory{T <: Real} <: AbstractMaterial
    drivers::GenericMemoryDriverState{T} = GenericMemoryDriverState{T}()
    ddrivers::GenericMemoryDriverState{T} = GenericMemoryDriverState{T}()
    variables::GenericMemoryVariableState{T} = GenericMemoryVariableState{T}()
    variables_new::GenericMemoryVariableState{T} = GenericMemoryVariableState{T}()
    parameters::GenericMemoryParameterState{T} = GenericMemoryParameterState{T}()
    dparameters::GenericMemoryParameterState{T} = GenericMemoryParameterState{T}()
    options::MemoryOptions = MemoryOptions()
end

MemoryDriverState = GenericMemoryDriverState{Float64}
MemoryParameterState = GenericMemoryParameterState{Float64}
MemoryVariableState = GenericMemoryVariableState{Float64}
Memory = GenericMemory{Float64}

"""
    state_to_vector(sigma::U, R::T, X1::U, X2::U) where U <: Symm2{T} where T <: Real

Adaptor for `nlsolve`. Marshal the problem state into a `Vector`.
"""
function state_to_vector(sigma::U, R::T, X1::U, X2::U) where U <: Symm2{T} where T <: Real
    return [tovoigt(sigma); R; tovoigt(X1); tovoigt(X2)]::Vector{T}
end

"""
    state_from_vector(x::AbstractVector{<:Real})

Adaptor for `nlsolve`. Unmarshal the problem state from a `Vector`.
"""
function state_from_vector(x::AbstractVector{T}) where T <: Real
    sigma::Symm2{T} = fromvoigt(Symm2{T}, @view x[1:6])
    R::T = x[7]
    X1::Symm2{T} = fromvoigt(Symm2{T}, @view x[8:13])
    X2::Symm2{T} = fromvoigt(Symm2{T}, @view x[14:19])
    return sigma, R, X1, X2
end

"""
    strain_memory_explicit_update(q, zeta, plastic_strain, dp, cumeq, pt, n, eta, xi, m)

Internal helper function for what it says on the tin.

Return `(dq, dzeta)`, the computed increments for `q` and `zeta`
for the given input.
"""
function strain_memory_explicit_update(q, zeta, plastic_strain, dp, cumeq, pt, n, eta, xi, m)
    dq = zero(q)
    dzeta = zero(zeta)
    JF = sqrt(1.5)*norm(dev(plastic_strain - zeta))
    FF = 2.0/3.0*JF - q
    if FF > 0.0
        nF = 1.5*dev(plastic_strain - zeta)/JF
        nnF = dcontract(n, nF)
        if nnF > 0
            dq = 2.0/3.0*eta*nnF*dp
            dzeta = 2.0/3.0*(1.0 - eta)*nnF*nF*dp
        end
    else
        # Memory evanescence term
        if cumeq >= pt
            dq = -xi*q^m*dp
        end
    end
    return dq, dzeta
end

"""
    integrate_material!(material::GenericMemory{T}) where T <: Real

Material model with a strain memory effect.

This is similar to the Chaboche material with two backstresses, with both
kinematic and isotropic hardening, but this model also features a strain
memory term.

Strain memory is used to be able to model strain amplitude-dependent isotropic
hardening. In practice, the transition from a tensile test curve to cyclic
behavior can be captured with this model.

See:

    D. Nouailhas, J.-L. Chaboche, S. Savalle, G. Cailletaud. On the constitutive
    equations for cyclic plasticity under nonproportional loading. International
    Journal of Plasticity 1(4) (1985), 317--330.
    https://doi.org/10.1016/0749-6419(85)90018-X
"""
function integrate_material!(material::GenericMemory{T}) where T <: Real
    p = material.parameters
    v = material.variables
    dd = material.ddrivers
    d = material.drivers
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q0, QM, mu, b, eta, m, pt, xi = p
    lambda, elastic_mu = lame(E, nu)

    @unpack strain, time = d
    dstrain = dd.strain
    dtime = dd.time
    @unpack stress, X1, X2, plastic_strain, cumeq, R, q, zeta, jacobian = v

    # elastic part
    jacobian = isotropic_elasticity_tensor(lambda, elastic_mu)
    stress += dcontract(jacobian, dstrain)

    # resulting deviatoric plastic stress (accounting for backstresses Xm)
    seff_dev = dev(stress - X1 - X2)
    # von Mises yield function
    f = sqrt(1.5)*norm(seff_dev) - (R0 + R)  # using elastic trial problem state
    if f > 0.0
        # Explicit update to memory-surface
        g! = create_nonlinear_system_of_equations(material)
        x0 = state_to_vector(stress, R, X1, X2)
        res = nlsolve(g!, x0; method=material.options.nlsolve_method, autodiff = :forward)
        converged(res) || error("Nonlinear system of equations did not converge!")
        x = res.zero
        stress, R, X1, X2 = state_from_vector(x)

        # using the new problem state
        seff_dev = dev(stress - X1 - X2)
        f = sqrt(1.5)*norm(seff_dev) - (R0 + R)

        dotp = ((f >= 0.0 ? f : 0.0)/Kn)^nn
        dp = dotp*dtime
        n = sqrt(1.5)*seff_dev/norm(seff_dev)

        plastic_strain += dp*n
        cumeq += dp

        dq, dzeta = strain_memory_explicit_update(q, zeta, plastic_strain, dp, cumeq, pt, n, eta, xi, m)
        q += dq
        zeta += dzeta

        # Compute the new Jacobian, accounting for the plastic contribution.
        drdx = ForwardDiff.jacobian(debang(g!), x)
        drde = zeros((length(x),6))
        drde[1:6, 1:6] = -tovoigt(jacobian)
        jacobian = fromvoigt(Symm4, (drdx\drde)[1:6, 1:6])
    end
    variables_new = GenericMemoryVariableState{T}(stress = stress,
                                                  X1 = X1,
                                                  X2 = X2,
                                                  R = R,
                                                  plastic_strain = plastic_strain,
                                                  cumeq = cumeq,
                                                  q = q,
                                                  zeta = zeta,
                                                  jacobian = jacobian)
    material.variables_new = variables_new
    return nothing
end

"""
    create_nonlinear_system_of_equations(material::GenericMemory{T}) where T <: Real

Create and return an instance of the equation system for the incremental form of
the evolution equations of the Memory material.

Used internally for computing the plastic contribution in `integrate_material!`.

The input `material` represents the problem state at the end of the previous
timestep. The created equation system will hold its own copy of that state.

The equation system is represented as a mutating function `g!` that computes the
residual:

```julia
    g!(F::V, x::V) where V <: AbstractVector{<:Real}
```

Both `F` (output) and `x` (input) are length-19 vectors containing
[sigma, R, X1, X2], in that order. The tensor quantities sigma, X1,
X2 are encoded in Voigt format.

The function `g!` is intended to be handed over to `nlsolve`.
"""
function create_nonlinear_system_of_equations(material::GenericMemory{T}) where T <: Real
    p = material.parameters
    v = material.variables
    dd = material.ddrivers
    d = material.drivers
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q0, QM, mu, b, eta, m, pt, xi = p
    lambda, elastic_mu = lame(E, nu)

    # Old problem state (i.e. the problem state at the time when this equation
    # system instance was created).
    #
    # Note this does not include the elastic trial; this is the state at the
    # end of the previous timestep.
    @unpack strain, time = d
    dstrain = dd.strain
    dtime = dd.time
    @unpack stress, X1, X2, plastic_strain, cumeq, R, q, zeta, jacobian = v

    jacobian = isotropic_elasticity_tensor(lambda, elastic_mu)

    # Explicit update of memory surface.
    #
    # Compute the residual. F is output, x is filled by NLsolve.
    # The solution is x = x* such that g(x*) = 0.
    function g!(F::V, x::V) where V <: AbstractVector{<:Real}
        stress_new, R_new, X1_new, X2_new = state_from_vector(x)  # tentative new values from nlsolve

        seff_dev = dev(stress_new - X1_new - X2_new)
        f = sqrt(1.5)*norm(seff_dev) - (R0 + R_new)

        dotp = ((f >= 0.0 ? f : 0.0)/Kn)^nn
        dp = dotp*dtime
        n = sqrt(1.5)*seff_dev/norm(seff_dev)

        dstrain_plastic = dp*n

        # Strain memory - explicit update
        plastic_strain_new = plastic_strain + dstrain_plastic
        dq, dzeta = strain_memory_explicit_update(q, zeta, plastic_strain_new, dp, cumeq, pt, n, eta, xi, m)
        q_new = q + dq
        zeta_new = zeta + dzeta

        # The equations are written in an incremental form.
        # TODO: multiply the equations by -1 to make them easier to understand in the context of the rest of the model.

        dstrain_elastic = dstrain - dstrain_plastic
        tovoigt!(view(F, 1:6), stress - stress_new + dcontract(jacobian, dstrain_elastic))

        F[7] = R - R_new + b*((QM + (Q0 - QM)*exp(-2.0*mu*q_new)) - R_new)*dp

        tovoigt!(view(F, 8:13), X1 - X1_new + dp*(2.0/3.0*C1*n - D1*X1_new))
        tovoigt!(view(F, 14:19), X2 - X2_new + dp*(2.0/3.0*C2*n - D2*X2_new))
        return nothing
    end
    return g!
end

end


text:
using SparseDiffTools
include("test_data.jl")

for i = 1:5
    g = graphs[i]
    out_colors = greedy_d1(g)
    for v = 1:length(g.vertices)
        color = out_colors[v]
        for j in neighbors(g, v)
            if out_colors[j] == color
                 @test false
            end
        end
    end
    @test true
end

ext: jl


text:
include("utils.jl")

root = dirname(dirname(@__FILE__))
config = Dict{String, Vector{Float64}}()
checkpoint = loadcheckpoint("params_final")

for (id, subdict) in checkpoint
    params = subdict["params"]
    error = subdict["kstep_mse"]
    params = params[error.!==nothing]
    error = error[error.!==nothing]
    bestparams = params[argmin(error)]
    config[id] = bestparams
end

open(joinpath(root, "config", "config.json"),"w") do f
    JSON.print(f, config)
end

println("Config saved!")


text:
using NURBS.Bases

export Curve, BCurve, NCurve, evalpoint, evaltangent, evalderiv1

struct Curve{B<:Basis1D, T<:Real} <: Object1D
    basis::B
    d1b::B
    points::Array{T, 2}

    function Curve(basis::B, points::Array{T, 2}) where {B<:Basis1D, T<:Real}
        d1b = deriv(basis)
        new{B, T}(basis, d1b, points)
    end
end

const BCurve = Curve{BSplineBasis}
const NCurve = Curve{NURBSBasis}

function evalpoint(crv::Curve, u)
    return crv.basis(u, crv.points)
end

function evaltangent(crv::Curve, u)
    t = crv.d1b(u, crv.points)'
    t ./ norm(t)
end

function evalderiv1(crv::Curve, u)
    return crv.d1b(u, crv.points)'
end


text:
using Documenter, Pkg3

makedocs(
    modules = [Pkg3],
    format = :html,
    sitename = "Pkg3.jl",
    pages = Any[
        "Pkg3" => "index.md",
    ]
)

deploydocs(
    repo = "github.com/JuliaLang/Pkg3.jl",
    target = "build",
    julia = "nightly",
    deps = nothing,
    make = nothing,
)

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
# constructs f which is the transition matrix between different skill levels
# also contstructs the sgrid (agin this uses uniform weights where the weights sum to the difference between shi and slo, not 1)
function persistent_skill_process(sH_over_sL::AbstractFloat, pLH::AbstractFloat,
                                  pHL::AbstractFloat, ns::Int)
    f1 = [[1-pLH pLH];[pHL 1-pHL]] # f1[i,j] is prob of going from i to j
    ss_skill_distr = [pHL/(pLH+pHL); pLH/(pLH+pHL)]
    slo    = 1.0 / (ss_skill_distr'*[1;sH_over_sL])
    sgrid  = slo*[1;sH_over_sL]
    sscale = sgrid[2] - sgrid[1]
    swts   = (sscale/ns)*ones(ns) # Quadrature weights
    f      = f1 ./ repeat(swts',ns,1)
    return f, sgrid, swts, sscale
end

# constructs the cash on hand grid (known as both agrid and xgrid and kinda used interchangebly throughout the code
function cash_grid(sgrid::AbstractArray, ω::AbstractFloat, H::AbstractFloat,
                   r::AbstractFloat, η::AbstractFloat, γ::AbstractFloat,
                   T::AbstractFloat, zlo::AbstractFloat, na::Int)
    smin = minimum(sgrid)*zlo                                   # lowest possible skill
    xlo_ss = ω*smin*H - (1+r)*η*exp(-γ) + T + sgrid[1]*ω*H*0.05 # lowest SS possible cash on hand

    xlo = xlo_ss                        # lower bound on cash on hand - could be < xlo_ss
    xhi = max(xlo*2, xlo + 12.0)         # upper bound on cash on hand
    xscale = (xhi-xlo)                  # size of w grids

    # Make grids
    xgrid = collect(range(xlo,stop = xhi, length = na)) # Evenly spaced grid
    xwts  = (xscale/na)*ones(na)          # Quadrature weights
    return xgrid, xwts, xlo, xhi, xscale
end

"""
```
generate_us_and_zs(ni, nz)
```

There is no need to recall this function, unless one wants to undo the seeding
in all past saved output.
"""
function generate_us_and_zs(ni, nz)
    us = rand(ni, 8)
    uz = rand(ni, 8)

    zgrid  = collect(range(0., stop = 2., length = nz))
    zprob  = [2*mollifier_hetdsgegovdebt(zgrid[i], 2., 0.) / nz for i=1:nz]
    zprob /= sum(zprob)

    zcdf = cumsum(zprob)
    zave = 0.5 * zgrid[1:nz-1] + 0.5 * zgrid[2:nz]
    zs   = zsample(uz, zgrid, zcdf, ni, nz)

    return us, zs
    #=
    JLD2.jldopen("$HETDSGEGOVDEBT/reference/us_zs.jld2", true, true, true, IOStream) do file
        file["us"] = us
        file["zs"] = zs
    end
    =#
end


text:
module Rankings

"""
    @ranking_str -> Dict

Parse human description of the ranking into dictionary
of alternatives' names as keys and their priorities
in the ranking as values.

# Examples
```jldoctest
julia> print(ranking"a₃,a₂,a₁∼a₄")
Dict("a₃"=>0,"a₁"=>-2,"a₄"=>-2,"a₂"=>-1)

julia> print(ranking"B (E A)(C D) F")
Dict("B"=>0,"A"=>-1,"C"=>-2,"D"=>-2,"E"=>-1,"F"=>-3)

julia> print(ranking"(O₂,O₃),O₄,(O₁,O₅)")
Dict("O₄"=>-1,"O₃"=>0,"O₁"=>-2,"O₅"=>-2,"O₂"=>0)

julia> print(ranking"a≻b≻d≻c")
Dict("c"=>-3,"b"=>-1,"a"=>0,"d"=>-2)

julia> print(ranking"AIFDHCBEJG")
Dict("B"=>-6,"A"=>0,"I"=>-1,"J"=>-8,"C"=>-5,"D"=>-3,"G"=>-9,"E"=>-7,"F"=>-2,"H"=>-4)

julia> print(ranking"54132")
Dict("4"=>-1,"1"=>-2,"5"=>0,"2"=>-4,"3"=>-3)
```
"""
:(@ranking_str)
:(@Ранжирование_str)

using Tokenize
using LinearAlgebra: ⋅

export @ranking_str, ⦸, @Ранжирование_str, consensus, ⊖

include("op_obslash.jl")
include("parse_data.jl")
include("negs_utils.jl")
include("consensus.jl")

macro ranking_str(s::String) ; parse_ranking(s) ; end
macro Ранжирование_str(s::String) ; parse_ranking(s) ; end

end # module


text:
# Hello, World with Julia
using DataFrames
print("Hello Ajay!")

max_stars_repo_name: ajk-lab/WsJulia

max_issues_repo_name: ajk-lab/WsJulia

max_forks_repo_name: ajk-lab/WsJulia


text:
@safetestset "Escaping quotes" begin
  @safetestset "Double quoted arguments" begin
    using Genie, Genie.Renderer
    using Genie.Renderer.Html

    @test Html.html("""
Good morning
""").body |> String ==
      """
Good morning
"""

    @test Html.html("""
Good morning
""").body |> String ==
      """
Good morning
"""

    @test Html.html("""
Good morning
""").body |> String ==
      """
Good morning
"""

    @test Html.html("""
Good morning
""").body |> String ==
      """
Good morning
"""
  end

  @safetestset "Single quoted arguments" begin
    using Genie, Genie.Renderer
    using Genie.Renderer.Html

    @test Html.html("""
Good morning
""").body |> String ==
      """
Good morning
"""

    @test Html.html("""
Good morning
""").body |> String ==
      """
Good morning
"""

    @test Html.html("""
Good morning
""").body |> String ==
      """
Good morning
"""
  end

  @safetestset "Arguments in templates" begin
    using Genie, Genie.Renderer
    using Genie.Renderer.Html

    @test Html.html(filepath("views/argsencoded1.jl.html")).body |> String ==
      """<!DOCTYPE html>
Greetings
"""

    @test_throws LoadError Html.html(filepath("views/argsencoded2.jl.html")).body |> String ==
      """<!DOCTYPE html>
Greetings
"""

    @test_throws LoadError Html.html(filepath("views/argsencoded3.jl.html")).body |> String ==
      """<!DOCTYPE html>
Greetings
"""

    @test Html.html(filepath("views/argsencoded1.jl.html"), layout=filepath("views/layoutargsencoding.jl.html")).body |> String ==
      """<!DOCTYPE html>

Layout header

Greetings

Layout footer

"""
  end
end


text:
"""
    edge_index(g::GNNGraph)

Return a tuple containing two vectors, respectively storing
the source and target nodes for each edges in `g`.

```julia
s, t = edge_index(g)
```
"""
edge_index(g::GNNGraph{<:COO_T}) = g.graph[1:2]

edge_index(g::GNNGraph{<:ADJMAT_T}) = to_coo(g.graph, num_nodes=g.num_nodes)[1][1:2]

get_edge_weight(g::GNNGraph{<:COO_T}) = g.graph[3]

get_edge_weight(g::GNNGraph{<:ADJMAT_T}) = to_coo(g.graph, num_nodes=g.num_nodes)[1][3]

Graphs.edges(g::GNNGraph) = zip(edge_index(g)...)

Graphs.edgetype(g::GNNGraph) = Tuple{Int, Int}

# """
#     eltype(g::GNNGraph)
#
# Type of nodes in `g`,
# an integer type like `Int`, `Int32`, `Uint16`, ....
# """
function Base.eltype(g::GNNGraph{<:COO_T})
    s, t = edge_index(g)
    w = get_edge_weight(g)
    return w !== nothing ? eltype(w) : eltype(s)
end

Base.eltype(g::GNNGraph{<:ADJMAT_T}) = eltype(g.graph)

function Graphs.has_edge(g::GNNGraph{<:COO_T}, i::Integer, j::Integer)
    s, t = edge_index(g)
    return any((s .== i) .& (t .== j))
end

Graphs.has_edge(g::GNNGraph{<:ADJMAT_T}, i::Integer, j::Integer) = g.graph[i,j] != 0

graph_type_symbol(g::GNNGraph{<:COO_T}) = :coo
graph_type_symbol(g::GNNGraph{<:SPARSE_T}) = :sparse
graph_type_symbol(g::GNNGraph{<:ADJMAT_T}) = :dense

Graphs.nv(g::GNNGraph) = g.num_nodes
Graphs.ne(g::GNNGraph) = g.num_edges
Graphs.has_vertex(g::GNNGraph, i::Int) = 1 <= i <= g.num_nodes
Graphs.vertices(g::GNNGraph) = 1:g.num_nodes

function Graphs.neighbors(g::GNNGraph, i; dir=:out)
    @assert dir ∈ (:in, :out)
    if dir == :out
        outneighbors(g, i)
    else
        inneighbors(g, i)
    end
end

function Graphs.outneighbors(g::GNNGraph{<:COO_T}, i::Integer)
    s, t = edge_index(g)
    return t[s .== i]
end

function Graphs.outneighbors(g::GNNGraph{<:ADJMAT_T}, i::Integer)
    A = g.graph
    return findall(!=(0), A[i,:])
end

function Graphs.inneighbors(g::GNNGraph{<:COO_T}, i::Integer)
    s, t = edge_index(g)
    return s[t .== i]
end

function Graphs.inneighbors(g::GNNGraph{<:ADJMAT_T}, i::Integer)
    A = g.graph
    return findall(!=(0), A[:,i])
end

Graphs.is_directed(::GNNGraph) = true
Graphs.is_directed(::Type{<:GNNGraph}) = true

"""
    adjacency_list(g; dir=:out)
    adjacency_list(g, nodes; dir=:out)

Return the adjacency list representation (a vector of vectors)
of the graph `g`.

Calling `a` the adjacency list, if `dir=:out` than
`a[i]` will contain the neighbors of node `i` through
outgoing edges. If `dir=:in`, it will contain neighbors from
incoming edges instead.

If `nodes` is given, return the neighborhood of the nodes in `nodes` only.
"""
function adjacency_list(g::GNNGraph, nodes; dir=:out, with_eid=false)
    @assert dir ∈ [:out, :in]
    s, t = edge_index(g)
    if dir == :in
        s, t = t, s
    end
    T = eltype(s)
    idict = 0
    dmap = Dict(n => (idict += 1) for n in nodes)
    adjlist = [T[] for _=1:length(dmap)]
    eidlist = [T[] for _=1:length(dmap)]
    for (eid, (i, j)) in enumerate(zip(s, t))
        inew = get(dmap, i, 0)
        inew == 0 && continue
        push!(adjlist[inew], j)
        push!(eidlist[inew], eid)
    end
    if with_eid
        return adjlist, eidlist
    else
        return adjlist
    end
end

# function adjacency_list(g::GNNGraph, nodes; dir=:out)
#     @assert dir ∈ [:out, :in]
#     fneighs = dir == :out ? outneighbors : inneighbors
#     return [fneighs(g, i) for i in nodes]
# end

adjacency_list(g::GNNGraph; dir=:out) = adjacency_list(g, 1:g.num_nodes; dir)

"""
    adjacency_matrix(g::GNNGraph, T=eltype(g); dir=:out, weighted=true)

Return the adjacency matrix `A` for the graph `g`.

If `dir=:out`, `A[i,j] > 0` denotes the presence of an edge from node `i` to node `j`.
If `dir=:in` instead, `A[i,j] > 0` denotes the presence of an edge from node `j` to node `i`.

User may specify the eltype `T` of the returned matrix.

If `weighted=true`, the `A` will contain the edge weigths if any, otherwise the elements of `A` will be either 0 or 1.
"""
function Graphs.adjacency_matrix(g::GNNGraph{<:COO_T}, T::DataType=eltype(g); dir=:out, weighted=true)
    if g.graph[1] isa CuVector
        # TODO revisit after https://github.com/JuliaGPU/CUDA.jl/pull/1152
        A, n, m = to_dense(g.graph, T; num_nodes=g.num_nodes, weighted)
    else
        A, n, m = to_sparse(g.graph, T; num_nodes=g.num_nodes, weighted)
    end
    @assert size(A) == (n, n)
    return dir == :out ? A : A'
end

function Graphs.adjacency_matrix(g::GNNGraph{<:ADJMAT_T}, T::DataType=eltype(g); dir=:out, weighted=true)
    @assert dir ∈ [:in, :out]
    A = g.graph
    if !weighted
        A = binarize(A)
    end
    A = T != eltype(A) ? T.(A) : A
    return dir == :out ? A : A'
end

function _get_edge_weight(g, edge_weight)
    if edge_weight === true || edge_weight === nothing
        ew = get_edge_weight(g)
    elseif edge_weight === false
        ew = nothing
    elseif edge_weight isa AbstractVector
        ew = edge_weight
    else
        error("Invalid edge_weight argument.")
    end
    return ew
end

"""
    degree(g::GNNGraph, T=nothing; dir=:out, edge_weight=true)

Return a vector containing the degrees of the nodes in `g`.

# Arguments
- `g`: A graph.
- `T`: Element type of the returned vector. If `nothing`, is
       chosen based on the graph type and will be an integer
       if `edge_weight=false`.
- `dir`: For `dir=:out` the degree of a node is counted based on the outgoing edges.
         For `dir=:in`, the ingoing edges are used. If `dir=:both` we have the sum of the two.
- `edge_weight`: If `true` and the graph contains weighted edges, the degree will
                be weighted. Set to `false` instead to just count the number of
                outgoing/ingoing edges.
                In alternative, you can also pass a vector of weights to be used
                instead of the graph's own weights.
"""
function Graphs.degree(g::GNNGraph{<:COO_T}, T::TT=nothing; dir=:out, edge_weight=true) where
        TT<:Union{Nothing, Type{<:Number}}
    s, t = edge_index(g)

    edge_weight = _get_edge_weight(g, edge_weight)
    edge_weight = edge_weight === nothing ? ones_like(s) : edge_weight

    T = isnothing(T) ? eltype(edge_weight) : T
    degs = fill!(similar(s, T, g.num_nodes), 0)

    if dir ∈ [:out, :both]
        degs = degs .+ NNlib.scatter(+, edge_weight, s, dstsize=(g.num_nodes,))
    end
    if dir ∈ [:in, :both]
        degs = degs .+ NNlib.scatter(+, edge_weight, t, dstsize=(g.num_nodes,))
    end
    return degs
end

# TODO:: Make efficient
Graphs.degree(g::GNNGraph, i::Union{Int, AbstractVector}; dir=:out) = degree(g; dir)[i]

function Graphs.degree(g::GNNGraph{<:ADJMAT_T}, T::TT=nothing; dir=:out, edge_weight=true) where TT
        TT<:Union{Nothing, Type{<:Number}}
    # edge_weight=true or edge_weight=nothing act the same here
    @assert !(edge_weight isa AbstractArray) "passing the edge weights is not support by adjacency matrix representations"
    @assert dir ∈ (:in, :out, :both)
    if T === nothing
        Nt = eltype(g)
        if edge_weight === false && !(Nt <: Integer)
            T = Nt == Float32 ? Int32 :
                Nt == Float16 ? Int16 : Int
        else
            T = Nt
        end
    end
    A = adjacency_matrix(g)
    if edge_weight === false
        A = binarize(A)
    end
    A = eltype(A) != T ? T.(A) : A
    return dir == :out ? vec(sum(A, dims=2)) :
           dir == :in  ? vec(sum(A, dims=1)) :
                  vec(sum(A, dims=1)) .+ vec(sum(A, dims=2))
end

function Graphs.laplacian_matrix(g::GNNGraph, T::DataType=eltype(g); dir::Symbol=:out)
    A = adjacency_matrix(g, T; dir=dir)
    D = Diagonal(vec(sum(A; dims=2)))
    return D - A
end

"""
    normalized_laplacian(g, T=Float32; add_self_loops=false, dir=:out)

Normalized Laplacian matrix of graph `g`.

# Arguments

- `g`: A `GNNGraph`.
- `T`: result element type.
- `add_self_loops`: add self-loops while calculating the matrix.
- `dir`: the edge directionality considered (:out, :in, :both).
"""
function normalized_laplacian(g::GNNGraph, T::DataType=Float32;
                        add_self_loops::Bool=false, dir::Symbol=:out)
    Ã = normalized_adjacency(g, T; dir, add_self_loops)
    return I - Ã
end

function normalized_adjacency(g::GNNGraph, T::DataType=Float32;
                        add_self_loops::Bool=false, dir::Symbol=:out)
    A = adjacency_matrix(g, T; dir=dir)
    if add_self_loops
        A = A + I
    end
    degs = vec(sum(A; dims=2))
    inv_sqrtD = Diagonal(inv.(sqrt.(degs)))
    return inv_sqrtD * A * inv_sqrtD
end

@doc raw"""
    scaled_laplacian(g, T=Float32; dir=:out)

Scaled Laplacian matrix of graph `g`,
defined as ``\hat{L} = \frac{2}{\lambda_{max}} L - I`` where ``L`` is the normalized Laplacian matrix.

# Arguments

- `g`: A `GNNGraph`.
- `T`: result element type.
- `dir`: the edge directionality considered (:out, :in, :both).
"""
function scaled_laplacian(g::GNNGraph, T::DataType=Float32; dir=:out)
    L = normalized_laplacian(g, T)
    # @assert issymmetric(L) "scaled_laplacian only works with symmetric matrices"
    λmax = _eigmax(L)
    return  2 / λmax * L - I
end

# _eigmax(A) = eigmax(Symmetric(A)) # Doesn't work on sparse arrays
function _eigmax(A)
    x0 = _rand_dense_vector(A)
    KrylovKit.eigsolve(Symmetric(A), x0, 1, :LR)[1][1] # also eigs(A, x0, nev, mode) available
end

_rand_dense_vector(A::AbstractMatrix{T}) where T = randn(float(T), size(A, 1))
_rand_dense_vector(A::CUMAT_T)= CUDA.randn(size(A, 1))

# Eigenvalues for cuarray don't seem to be well supported.
# https://github.com/JuliaGPU/CUDA.jl/issues/154
# https://discourse.julialang.org/t/cuda-eigenvalues-of-a-sparse-matrix/46851/5

"""
    graph_indicator(g)

Return a vector containing the graph membership
(an integer from `1` to `g.num_graphs`) of each node in the graph.
"""
function graph_indicator(g; edges=false)
    if isnothing(g.graph_indicator)
        gi = ones_like(edge_index(g)[1], Int, g.num_nodes)
    else
        gi = g.graph_indicator
    end
    if edges
        s, t = edge_index(g)
        return gi[s]
    else
        return gi
    end
end

function node_features(g::GNNGraph)
    if isempty(g.ndata)
        return nothing
    elseif length(g.ndata) > 1
        @error "Multiple feature arrays, access directly through `g.ndata`"
    else
        return g.ndata[1]
    end
end

function edge_features(g::GNNGraph)
    if isempty(g.edata)
        return nothing
    elseif length(g.edata) > 1
        @error "Multiple feature arrays, access directly through `g.edata`"
    else
        return g.edata[1]
    end
end

function graph_features(g::GNNGraph)
    if isempty(g.gdata)
        return nothing
    elseif length(g.gdata) > 1
        @error "Multiple feature arrays, access directly through `g.gdata`"
    else
        return g.gdata[1]
    end
end

"""
    is_bidirected(g::GNNGraph)

Check if the directed graph `g` essentially corresponds
to an undirected graph, i.e. if for each edge it also contains the
reverse edge.
"""
function is_bidirected(g::GNNGraph)
    s, t = edge_index(g)
    s1, t1 = sort_edge_index(s, t)
    s2, t2 = sort_edge_index(t, s)
    all((s1 .== s2) .& (t1 .== t2))
end

"""
    has_self_loops(g::GNNGraph)

Return `true` if `g` has any self loops.
"""
function Graphs.has_self_loops(g::GNNGraph)
    s, t = edge_index(g)
    any(s .== t)
end

"""
    has_multi_edges(g::GNNGraph)

Return `true` if `g` has any multiple edges.
"""
function has_multi_edges(g::GNNGraph)
    s, t = edge_index(g)
    idxs, _ = edge_encoding(s, t, g.num_nodes)
    length(union(idxs)) < length(idxs)
end

@non_differentiable edge_index(x...)
@non_differentiable adjacency_list(x...)
@non_differentiable graph_indicator(x...)
@non_differentiable has_multi_edges(x...)
@non_differentiable Graphs.has_self_loops(x...)
@non_differentiable is_bidirected(x...)
@non_differentiable normalized_adjacency(x...) # TODO remove this in the future
@non_differentiable normalized_laplacian(x...) # TODO remove this in the future
@non_differentiable scaled_laplacian(x...) # TODO remove this in the future


text:
totoro_coords=[484 800;
441 796;
485 794;
480 793;
522 791;
465 789;
419 787;
406 785;
474 785;
520 785;
540 784;
466 783;
519 782;
429 781;
464 781;
459 780;
480 779;
462 777;
583 774;
426 771;
456 771;
394 770;
579 770;
384 768;
511 768;
601 768;
464 767;
446 766;
493 766;
509 766;
571 766;
609 765;
587 764;
589 764;
410 763;
521 763;
575 763;
464 761;
580 761;
543 760;
594 760;
530 759;
414 758;
448 757;
430 756;
617 756;
626 756;
497 755;
507 755;
410 754;
497 754;
504 754;
506 754;
507 754;
533 754;
566 754;
437 753;
543 752;
563 752;
601 752;
610 752;
454 751;
513 750;
559 750;
566 750;
603 750;
465 748;
546 748;
413 747;
488 747;
491 747;
583 747;
500 746;
401 745;
528 745;
646 744;
454 743;
600 743;
444 742;
631 742;
486 739;
532 739;
490 738;
494 738;
546 738;
652 738;
545 737;
600 737;
437 736;
440 736;
631 736;
433 735;
508 735;
583 734;
613 734;
490 733;
562 733;
443 732;
503 731;
666 731;
566 730;
629 730;
451 728;
437 727;
567 727;
618 727;
541 726;
671 726;
454 725;
553 725;
600 725;
449 724;
477 724;
556 724;
574 724;
474 723;
499 723;
584 723;
584 722;
672 722;
588 721;
470 720;
471 720;
510 719;
523 719;
562 719;
665 717;
671 717;
520 716;
580 716;
620 715;
626 715;
667 715;
683 715;
585 714;
567 713;
497 712;
554 712;
559 712;
544 711;
549 711;
622 709;
674 709;
644 708;
595 707;
620 707;
519 706;
542 706;
603 706;
653 705;
593 703;
533 701;
547 701;
528 700;
518 699;
625 699;
609 697;
638 697;
696 697;
509 695;
604 695;
636 694;
666 694;
684 694;
526 693;
565 693;
619 693;
634 693;
669 692;
658 691;
677 690;
697 690;
570 688;
576 688;
607 688;
649 687;
678 687;
694 685;
534 684;
567 683;
671 683;
622 682;
665 682;
671 682;
682 682;
631 681;
573 680;
681 679;
650 675;
577 674;
593 673;
631 672;
602 670;
533 669;
619 669;
645 667;
620 665;
691 665;
659 662;
530 660;
701 659;
705 658;
664 657;
327 654;
331 648;
170 639;
328 639;
155 636;
158 636;
328 634;
181 632;
336 632;
519 630;
346 628;
174 627;
331 627;
149 626;
324 626;
178 625;
327 625;
334 625;
307 624;
167 622;
309 620;
185 618;
180 615;
154 614;
171 611;
319 609;
192 607;
325 607;
142 606;
145 606;
302 606;
335 606;
310 605;
178 601;
187 601;
345 600;
180 598;
146 597;
174 597;
180 597;
297 594;
152 591;
295 591;
330 591;
188 590;
330 590;
505 588;
318 587;
502 587;
157 584;
168 579;
159 577;
163 574;
312 574;
332 573;
187 572;
307 572;
151 571;
324 570;
142 569;
151 569;
180 568;
184 568;
180 566;
190 566;
176 565;
324 565;
188 563;
145 560;
157 560;
186 560;
192 560;
227 560;
182 559;
192 559;
300 559;
308 559;
159 558;
226 558;
266 558;
175 557;
206 557;
237 557;
159 555;
280 554;
159 553;
187 551;
240 551;
278 551;
311 551;
226 550;
490 550;
251 549;
273 549;
159 547;
165 547;
180 547;
188 546;
245 546;
310 546;
206 544;
236 544;
271 544;
159 542;
217 542;
161 541;
178 541;
192 541;
229 540;
273 540;
293 540;
307 540;
157 539;
239 538;
261 538;
296 537;
297 537;
287 534;
485 534;
173 533;
203 532;
227 532;
321 532;
248 531;
216 529;
245 529;
286 529;
154 527;
197 527;
477 527;
314 526;
136 523;
218 523;
282 523;
288 523;
472 522;
312 521;
248 519;
256 519;
194 518;
217 518;
221 518;
160 517;
206 517;
240 516;
162 514;
245 514;
209 513;
236 513;
237 513;
143 512;
131 511;
135 511;
153 510;
491 510;
454 509;
242 508;
284 508;
142 507;
249 507;
331 506;
340 506;
271 505;
293 505;
347 505;
466 505;
121 504;
332 504;
346 504;
213 503;
486 503;
277 501;
489 501;
487 500;
200 499;
240 498;
275 498;
489 498;
438 497;
45 496;
141 496;
280 496;
121 495;
344 495;
352 495;
135 494;
145 493;
294 493;
458 493;
64 492;
346 492;
255 491;
429 491;
481 491;
90 490;
121 489;
171 489;
335 489;
421 489;
486 488;
218 487;
442 487;
150 486;
422 486;
488 486;
256 485;
337 484;
106 483;
216 483;
123 482;
139 482;
186 482;
413 482;
478 482;
110 481;
112 481;
254 480;
352 480;
183 479;
304 479;
467 479;
217 478;
428 478;
272 477;
209 476;
218 476;
324 476;
227 475;
236 475;
336 475;
355 475;
266 474;
297 474;
222 473;
225 473;
248 473;
253 473;
281 473;
439 473;
102 472;
149 472;
236 472;
284 472;
310 472;
494 472;
405 471;
207 470;
470 470;
36 469;
54 469;
84 469;
102 469;
103 469;
126 469;
166 469;
213 469;
313 469;
466 469;
476 469;
124 468;
185 468;
273 468;
363 468;
270 467;
108 465;
285 465;
359 465;
488 465;
233 464;
252 464;
270 464;
432 464;
279 463;
313 463;
104 462;
175 462;
201 462;
212 462;
391 462;
474 462;
396 461;
146 460;
148 460;
397 460;
332 459;
364 459;
101 458;
370 458;
342 457;
369 453;
439 453;
109 452;
402 452;
438 451;
471 451;
478 451;
371 450;
365 449;
477 449;
479 449;
353 447;
381 447;
395 447;
27 446;
32 445;
129 445;
373 445;
154 443;
380 443;
93 442;
130 441;
94 440;
424 439;
335 438;
97 437;
180 437;
354 437;
128 436;
367 436;
145 435;
148 435;
366 434;
380 434;
428 434;
173 433;
87 432;
117 432;
440 432;
163 431;
191 431;
406 431;
209 430;
315 430;
292 429;
126 428;
363 428;
467 428;
301 427;
361 427;
274 426;
246 425;
324 425;
329 425;
445 425;
123 422;
452 422;
164 421;
384 421;
439 421;
263 420;
291 420;
329 420;
368 420;
73 419;
310 419;
154 418;
369 418;
376 418;
434 418;
360 417;
456 417;
288 416;
440 416;
373 415;
438 415;
174 414;
254 414;
331 414;
338 414;
347 414;
63 413;
246 412;
384 412;
65 411;
395 411;
128 410;
373 410;
387 410;
80 409;
340 409;
91 408;
80 407;
82 407;
367 407;
394 406;
431 406;
444 406;
61 405;
382 405;
123 403;
320 403;
326 403;
103 402;
81 401;
58 400;
149 400;
376 399;
115 398;
380 398;
332 397;
359 397;
99 396;
361 396;
449 396;
115 395;
122 395;
400 395;
451 395;
138 394;
116 393;
377 392;
382 392;
423 392;
82 391;
52 389;
90 389;
133 389;
413 387;
61 385;
360 385;
106 384;
357 384;
55 383;
78 382;
379 382;
385 382;
60 380;
121 380;
390 380;
391 380;
67 379;
379 378;
76 377;
99 376;
405 376;
406 376;
386 375;
99 373;
231 373;
387 373;
59 372;
427 372;
427 371;
95 370;
100 370;
403 369;
439 369;
370 368;
401 368;
245 367;
378 366;
61 364;
104 364;
215 363;
46 362;
90 362;
416 362;
422 360;
382 359;
414 359;
41 358;
28 357;
386 357;
424 357;
64 356;
426 356;
428 356;
39 354;
42 354;
59 354;
155 353;
300 353;
396 353;
430 353;
72 352;
298 352;
327 351;
139 350;
422 350;
95 349;
62 348;
41 347;
44 346;
137 344;
76 342;
57 341;
59 341;
67 341;
350 341;
388 341;
405 341;
101 339;
400 336;
46 333;
89 332;
45 331;
98 331;
395 331;
407 330;
60 329;
63 328;
269 326;
40 325;
50 324;
86 324;
262 320;
32 319;
32 318;
75 317;
290 317;
291 316;
191 315;
353 315;
99 314;
252 313;
81 311;
259 311;
71 310;
75 310;
372 310;
96 309;
288 309;
44 308;
359 308;
338 307;
355 305;
21 303;
40 303;
57 303;
25 302;
111 302;
112 302;
346 302;
105 301;
129 301;
53 300;
31 299;
103 299;
173 299;
112 298;
114 298;
126 298;
97 297;
114 297;
5 296;
59 295;
56 294;
97 294;
133 294;
6 293;
70 293;
118 292;
129 292;
135 292;
410 292;
384 291;
11 287;
58 287;
59 287;
66 287;
93 287;
19 285;
56 282;
60 282;
74 279;
415 277;
82 276;
52 274;
28 269;
72 268;
27 267;
52 267;
77 266;
63 265;
41 264;
21 263;
410 259;
18 257;
3 255;
415 255;
1 249;
37 248;
47 248;
417 247;
47 243;
24 242;
87 241;
18 240;
36 240;
72 235;
33 233;
40 233;
55 233;
8 232;
11 228;
36 228;
88 228;
49 225;
83 223;
31 220;
412 220;
29 218;
91 217;
84 212;
80 208;
84 208;
36 206;
22 205;
68 204;
78 204;
20 203;
32 201;
42 201;
62 201;
77 201;
50 199;
92 199;
78 198;
39 196;
406 196;
48 191;
46 189;
32 188;
60 187;
82 187;
26 186;
49 186;
61 186;
56 183;
31 182;
31 180;
45 177;
35 175;
393 175;
104 173;
68 172;
40 171;
405 170;
407 170;
79 169;
82 169;
38 165;
50 164;
39 161;
108 160;
57 159;
132 158;
47 156;
87 154;
101 154;
55 152;
81 152;
122 150;
125 150;
94 147;
125 147;
55 146;
58 145;
67 145;
372 142;
78 141;
149 141;
368 140;
121 137;
149 137;
366 137;
102 136;
143 136;
90 135;
99 134;
136 134;
142 133;
67 132;
149 131;
169 130;
390 130;
82 129;
160 129;
383 127;
104 126;
118 125;
99 123;
360 123;
146 120;
97 119;
112 118;
361 117;
376 117;
96 116;
61 115;
185 115;
69 114;
135 114;
150 112;
371 112;
117 110;
62 109;
165 109;
203 107;
366 107;
93 106;
190 106;
215 105;
341 105;
354 105;
326 104;
132 103;
146 103;
345 103;
115 102;
168 102;
183 102;
369 102;
229 101;
334 101;
111 100;
231 100;
353 100;
162 98;
262 98;
104 96;
273 96;
352 96;
260 94;
267 94;
308 94;
225 93;
117 91;
153 91;
253 91;
298 91;
344 91;
135 90;
159 90;
149 87;
317 87;
90 85;
217 85;
286 85;
99 84;
201 83;
259 83;
299 83;
322 82;
239 81;
257 81;
265 81;
159 79;
281 79;
140 78;
163 78;
170 77;
99 75;
140 75;
212 75;
295 75;
202 74;
317 74;
162 73;
167 73;
140 72;
161 72;
203 72;
311 71;
203 70;
293 70;
185 68;
266 68;
132 67;
195 67;
239 66;
263 65;
292 65;
311 65;
311 64;
296 62;
124 61;
312 61;
154 58;
149 57;
281 57;
307 57;
179 56;
223 55;
309 55;
278 54;
133 50;
294 50;
191 49;
246 49;
265 49;
283 46;
147 45;
260 45;
189 44;
254 43;
268 43;
285 42;
184 40;
164 39;
182 36;
167 35;
199 34;
258 34;
165 33;
256 33;
167 32;
202 31;
172 30;
225 30;
264 30;
269 29;
274 29;
280 27;
166 25;
175 25;
147 22;
252 22;
168 20;
283 17;
274 15;
270 12]


text:
# Script to replicate table 4 and 5 (regression from equation 15 and 16)

using DataFrames, FixedEffectModels, RegressionTables, StatsBase, Statistics, LinearAlgebra, CSV

# -------------- Previous scripts -----------------------------------------------------------------------------------------------------------------------

include("src/data_wrangling.jl") # formats the raw input data obtained from WIOD

# -------------- Country-sector level I/O-Table ---------------------------------------------------------------------------------------------------------

# function to format matrix to DataFrame in main function "country_sector_level_head_reis_index"
function head_reis_df(result::Matrix, year::Int64, output::String)
    namestable = ["year", "row_country", "row_sector"]
    namestable = [namestable; repeat(ctrys, inner=n_sectors) .* "_" .* string.(repeat(1:35, outer=n_ctrys))]

    if (output == "ID")
        data = UpperTriangular(result)
    else
        data = UpperTriangular(result)
        for diagonal in 1:size(data, 1)
            data[diagonal, diagonal] = 0.0 # set diagonal entries = 0.0 as we filter them out later
        end
    end

    table = hcat(fill(year, n_ctrys*n_sectors), repeat(ctrys, inner=n_sectors), repeat(1:35, outer=n_ctrys), UpperTriangular(data))
    table = DataFrame(table, namestable)

    table_long = stack(table, Not([:year, :row_country, :row_sector]))
    table_long = subset(table_long, :value => ByRow(x -> x != 0.0))

    table_long = hcat(table_long, DataFrame(reduce(vcat, permutedims.(split.(table_long.variable, "_"))), [:col_country, :col_sector]))
    table_long = table_long[:, Not(:variable)]

    transform!(table_long, [:row_country, :col_country] .=> ByRow(string) .=> [:row_country, :col_country], renamecols=false) # change type to String
    transform!(table_long, [:year, :row_sector] .=> ByRow(Int) .=> [:year, :row_sector], renamecols=false) # change type to Int64
    transform!(table_long, :value .=> ByRow(Float64) .=> :value, renamecols=false) # change type to Float64
    transform!(table_long, :col_sector .=> ByRow(x -> parse(Int64,x)) .=> :col_sector, renamecols=false) # change type to Int64

    return table_long
end

function country_sector_level_head_reis_index(df::DataFrame, year::Int64, output::String)

    ID, FD, GO, VA, IV, other = obtain_matrices(df, year) # N*S×N*S, N*S×N, N*S×1, N*S×1, N*S×N, N*S×1

    # net inventory correction
    iv_correction = [GO[i]/(GO[i]-IV[i, j]) for i in 1:n_ctrys*n_sectors, j in 1:n_ctrys] # N*S×N, to match country wise inventory changes
    # notice difference in column size, need spread inventory correction of country column to country-sector columns
    # i.e. ID[1,1:35]/iv_correction[1,1], however problem with indexing => ceil(Int, 34/35)=1, ceil(Int, 36/35)=2
    ID_corrected = [ID[i, j]*iv_correction[i, ceil(Int, j/n_sectors)] for i in 1:n_ctrys*n_sectors, j in 1:n_ctrys*n_sectors] # N*S×N*S
    FD_corrected = [FD[i, j]*iv_correction[i, j] for i in 1:n_ctrys*n_sectors, j in 1:n_ctrys] # N*S×N

    replace!(ID_corrected, NaN => 0.0) # sometimes zeros are interpreted as NaN (particular for sector 35, private households with employed persons)
    replace!(FD_corrected, NaN => 0.0)

    # footnote 22 on page 21: substitute 1e-18 for any zero entries otherwise trade costs explode to infinity!
    ID = replace!(x -> x <= 0.0 ? 1e-18 : x, ID_corrected) # N*S×N*S
    FD = replace!(x -> x <= 0.0 ? 1e-18 : x, FD_corrected) # N*S×N

    # claculate Head-Reis index
    θ = 5 # assumption on trade elasticity
    τ_ID = [((ID[i, j] * ID[j, i]) / (ID[i, i] * ID[j, j]))^(-1/(2θ)) for i in 1:n_ctrys*n_sectors, j in 1:n_ctrys*n_sectors] # bilater country-sector trade costs
    τ_FD = [((FD[i, ceil(Int, j/n_sectors)] * FD[j, ceil(Int, i/n_sectors)]) / (FD[i, ceil(Int, i/n_sectors)] * FD[j, ceil(Int, j/n_sectors)]))^(-1/(2θ)) for i in 1:n_ctrys*n_sectors, j in 1:n_ctrys*n_sectors]

    result = output == "ID" ? τ_ID : τ_FD

    head_reis_df(result::Matrix, year::Int64, output::String)

end

t_tab4 = DataFrame(year=Int[], row_country=String[], row_sector=Int64[], col_country=String[], col_sector=Int64[], value=Float64[])
t_tab5 = DataFrame(year=Int[], row_country=String[], row_sector=Int64[], col_country=String[], col_sector=Int64[], value=Float64[])

for i in years
    append!(t_tab4, country_sector_level_head_reis_index(df, i, "ID"))
    append!(t_tab5, country_sector_level_head_reis_index(df, i, "FD"))
end

CSV.write("clean/t_tab4.csv", t_tab4) # use CSV since XLSX cannot store that many rows
CSV.write("clean/t_tab5.csv", t_tab5)

# -------------- Table 4 -----------------------------------------------------------------------------------------------------------------------------

t_tab4 = CSV.read("clean/t_tab4.csv", DataFrame) # also imports types correctly!
t_tab5 = CSV.read("clean/t_tab4.csv", DataFrame) # also imports types correctly!

function reg_tab_4_5(df::DataFrame)

    df[:, :row_origin_destination] .= df[:, :row_country] .* "_" .* string.(df[:, :row_sector]) # for standard error clustering
    df[:, :col_origin_destination] .= df[:, :row_country] .* "_" .* string.(df[:, :row_sector]) # for standard error clustering
    df[:, :origin_destination] .= df[:, :row_origin_destination] .* "___" .* df[:, :col_origin_destination] # for fixed effects
    transform!(df, :row_sector => ByRow(x -> x in goods_sectors ? 1 : 0) => :row_sector_dummy) # dummy for *reporting* sector, "goods" = 1

    df_goods = subset(df, :row_sector_dummy => ByRow(x -> x==1)) # only goods sectors
    df_services = subset(df, :row_sector_dummy => ByRow(x -> x!=1)) # only service sectors

    regression = Dict{String,Any}() # initializing dictionary to store regression results (Any could be replaced by type FixedEffectModel)
    df_names = ["all", "goods", "services"] # names of dictionary "keys"
    df_array = [df, df_goods, df_services] # names of DataFrames used

    for i in 1:length(df_names)
        push!(regression, "reg_1_" * df_names[i] =>
        FixedEffectModels.reg(df_array[i], @formula(log(value) ~ year + fe(origin_destination)), Vcov.cluster(:row_origin_destination, :col_origin_destination, :year)))

        push!(regression, "reg_2_" * df_names[i] =>
        FixedEffectModels.reg(df_array[i], @formula(log(value) ~ year + fe(origin_destination)), Vcov.cluster(:row_origin_destination, :col_origin_destination, :year),
            contrasts = Dict(:year => DummyCoding(base=1995)))) # treat years as dummies
    end

    reg = RegressionTables.regtable(regression["reg_1_all"], regression["reg_2_all"], regression["reg_1_goods"], regression["reg_2_goods"],
        regression["reg_1_services"], regression["reg_2_services"];
        renderSettings = asciiOutput(), regression_statistics=[:nobs, :r2], print_fe_section=true, estimformat="%0.4f") # output on the console

    return reg

end

reg_tab_4_5(t_tab4)

# -------------- Table 5 -----------------------------------------------------------------------------------------------------------------------------

reg_tab_4_5(t_tab5) # is not correct!
# should input should be a N*S×N by matrix not N*S×N*S
# results are still suprisingly close


text:
# Autogenerated wrapper script for WaveFD_jll for i686-linux-musl
export libillumination, libprop2DAcoIsoDenQ_DEO2_FDTD, libprop2DAcoTTIDenQ_DEO2_FDTD, libprop2DAcoVTIDenQ_DEO2_FDTD, libprop3DAcoIsoDenQ_DEO2_FDTD, libprop3DAcoTTIDenQ_DEO2_FDTD, libprop3DAcoVTIDenQ_DEO2_FDTD, libspacetime

using CompilerSupportLibraries_jll
using FFTW_jll
JLLWrappers.@generate_wrapper_header("WaveFD")
JLLWrappers.@declare_library_product(libillumination, "libillumination.so")
JLLWrappers.@declare_library_product(libprop2DAcoIsoDenQ_DEO2_FDTD, "libprop2DAcoIsoDenQ_DEO2_FDTD.so")
JLLWrappers.@declare_library_product(libprop2DAcoTTIDenQ_DEO2_FDTD, "libprop2DAcoTTIDenQ_DEO2_FDTD.so")
JLLWrappers.@declare_library_product(libprop2DAcoVTIDenQ_DEO2_FDTD, "libprop2DAcoVTIDenQ_DEO2_FDTD.so")
JLLWrappers.@declare_library_product(libprop3DAcoIsoDenQ_DEO2_FDTD, "libprop3DAcoIsoDenQ_DEO2_FDTD.so")
JLLWrappers.@declare_library_product(libprop3DAcoTTIDenQ_DEO2_FDTD, "libprop3DAcoTTIDenQ_DEO2_FDTD.so")
JLLWrappers.@declare_library_product(libprop3DAcoVTIDenQ_DEO2_FDTD, "libprop3DAcoVTIDenQ_DEO2_FDTD.so")
JLLWrappers.@declare_library_product(libspacetime, "libspacetime.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, FFTW_jll)
    JLLWrappers.@init_library_product(
        libillumination,
        "lib/libillumination.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprop2DAcoIsoDenQ_DEO2_FDTD,
        "lib/libprop2DAcoIsoDenQ_DEO2_FDTD.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprop2DAcoTTIDenQ_DEO2_FDTD,
        "lib/libprop2DAcoTTIDenQ_DEO2_FDTD.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprop2DAcoVTIDenQ_DEO2_FDTD,
        "lib/libprop2DAcoVTIDenQ_DEO2_FDTD.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprop3DAcoIsoDenQ_DEO2_FDTD,
        "lib/libprop3DAcoIsoDenQ_DEO2_FDTD.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprop3DAcoTTIDenQ_DEO2_FDTD,
        "lib/libprop3DAcoTTIDenQ_DEO2_FDTD.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprop3DAcoVTIDenQ_DEO2_FDTD,
        "lib/libprop3DAcoVTIDenQ_DEO2_FDTD.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libspacetime,
        "lib/libspacetime.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module LSHFunctions

using Distributions, LinearAlgebra, SparseArrays

#========================
Common types/utilities used throughout the module
========================#

include(joinpath("utils", "vecops.jl"))
include("LSHBase.jl")
include("intervals.jl")
include("similarities.jl")

#========================
Hash functions
========================#

include(joinpath("hashes", "simhash.jl"))
include(joinpath("hashes", "minhash.jl"))
include(joinpath("hashes", "lphash.jl"))
include(joinpath("hashes", "sign_alsh.jl"))

# Must be placed last, since it uses the definitions of LSHFunction subtypes
# defined in the other files.
include(joinpath("hashes", "lshfunction.jl"))

#========================
Function hashing
========================#

include(joinpath("function_hashing", "monte_carlo.jl"))

#========================
Exports
========================#

# Similarity functions, norms, inner products
export cossim, inner_prod, ℓ1, ℓ2, ℓp, L1, L2, Lp, ℓ1_norm, ℓ2_norm,
       ℓp_norm, L1_norm, L2_norm, Lp_norm, jaccard

# Hash functions
export SimHash, L1Hash, L2Hash, MIPSHash, SignALSH, MinHash,
       LSHFunction, MonteCarloHash, SymmetricLSHFunction,
       AsymmetricLSHFunction

# Helper / utility functions for LSHFunctions
export index_hash, query_hash, n_hashes, hashtype, similarity, lsh_family,
       embedded_similarity, collision_probability, @interval, HashCompressor

end # module


text:
# Implementation of BasicCounting algorithm from the paper
# Datar et.al.
# "Maintaining stream statistics over sliding windows"
# http://ilpubs.stanford.edu:8090/504/1/2001-34.pdf
#
using Base.Test

"""
 Count the number of 1's in the binary vector.
"""
function numOfOnes(w)
  h::Array{Int,1} = zeros(2)
  for e in w
    h[e+1] += 1
  end
  h[2]
end

"""
 Naive counting for the benchmark.
"""
function naiveCounting(seq_length = 100, w_length = 10)
    n::Array{Int64,1} = []
    y = rand([0,1], seq_length)
    arr = zeros(2)
    for t = w_length : length(y)
        window = y[ (t-w_length + 1) : t]
        n_of_ones = numOfOnes(window)
        push!(n, n_of_ones)
    end
    n
end

function hist_int(v::Array{Int,1})
    a=unique(v)
    hist::Array{Int,1} = zeros(maximum(a))
    for e in a
        hist[e]+=1
    end
    hist
end

"""
 ..each data element has an arrival time which increments by one at each arrival ,
 with the leftmost element considered to have arrived at time 1.

 timestamp corresponds to the position of an active data element in the current
 window.. we timestamp the active data elements from right to left, with the most
 recent element being at position 1.

 histograms are constructed by counting 1's in buckets by timestamps (?)

 IMPORTANT: Histogram must be updates so that Invariants 1,2 satisfied.

"""
function basic_counting(x::Array{Int, 1}, N=10)
    n = length(x)
    TOTAL::Int=0
    LAST::Int = 0
    k=2 # k=[1/ϵ]
    # p.4: "The estimate is TOTAL minus half of LAST"
    # timestamps are naturally given by indexing buckets[] array
    buckets::Array{Int,1} = zeros(N)
    arrival_times::Array{Int,1} = zeros(N)
    # expiration_times = zeros(N)
    max_num_of_buckets = (k/2.0+1.0)*(log2(2.0*N/k+1.0)+1)
    println("max_num_of_buckets=$max_num_of_buckets")
    for t = 1:n
        if x[t]==1
            TOTAL += 1
            # 1.. If the timestamp of the last bucket indicates expiry, delete
            # that bucket and update the counter LAST containing the size of the
            # last bycket and the counter TOTAL (total size of buckets)
            buckets = circshift(buckets,1)
            if buckets[1]!=0
                TOTAL += buckets[1]
            end
            buckets[1] = x[t]
            println(buckets)
            println(TOTAL)
            arrival_times[1] = t
            arrival_times = circshift(arrival_times,1)
            # println(arrival_times,"\n")
        end
    end
end

"""
 Merge just two buckets with given indices.
 ([1,1,2,4], 1, 2) -> [2,2,4,0]
 ([1,1,2,4,0,0], 1, 2) -> [2,2,4,0,0,0]
"""
function merge_buckets!(buckets::Array{Int,1}, idx1::Int, idx2::Int)
    val = buckets[idx1]
    @assert val == buckets[idx2]
    @assert idx1 <= idx2
    buckets[idx1] = 2 * val
    buckets[idx2] = 0
    buckets[idx2:end-1] = buckets[idx2+1:end]
    buckets[end]=0
    buckets
end

"""
 Collapse/merge buckets to keep the Invariant 2.
 Example:
 ([32,16,8,8,8,2,2,1,1], val=8, num=2) -> [32,16,16,8,2,2,1,1,0]
"""
function collapse_buckets(EH_histogram::Array{Int,1}; val=1, num=1)
    @assert num >= 2
    n = length(EH_histogram)
    out::Array{Int,1} = zeros(n)
    count::Int = 0
    first_val_index::Int = 0
    last_val_index::Int = 0
    last_non_zero::Int=0
    for i = 1:n
        if EH_histogram[i] == val
            count += 1
            last_val_index = i
        end
        if EH_histogram[i]==0
            last_non_zero = i-1
        end
    end
    @assert count >= num
    if last_non_zero == 0
        last_non_zero = n
    end
    first_val_index = last_val_index - count + 1
    # shift_left = count - num + 1
    out[1:first_val_index-1] = EH_histogram[1:first_val_index-1]
    out[first_val_index] = val*num
    out[first_val_index + 1 : first_val_index + 1 + count - num - 1] = val #todo
    # @assert sum(EH_histogram) == sum(out)
    ## Untouched tail
    width_of_untouched_tail = last_non_zero - last_val_index
    ## Index of the first element of the tail
    tail_first_el = first_val_index + count - num + 1
    out[tail_first_el : tail_first_el + width_of_untouched_tail-1] = EH_histogram[last_val_index+1:last_non_zero]
    out
end

function test_collapse_buckets()
    println("... Run test_collapse_buckets()")
    @test collapse_buckets([32,16,8,8,8,2,2,1,1], val=8, num=2) == [32,16,16,8,2,2,1,1,0]
    @test collapse_buckets([32,16,16,8,2,2,1,1,0], val=2, num=2) == [32,16,16,8,4,1,1,0,0]
    @test collapse_buckets([32,16,16,8,4,1,1,0,0], val=16, num=2) == [32,32,8,4,1,1,0,0,0]
    @test collapse_buckets([1,1,1,1,0], val=1, num=2) == [2,1,1,0,0]
    @test collapse_buckets([2,1,1,0,0], val=1, num=2) == [2,2,0,0,0]
    println(".. Stop")
end

"""
 Check the Invariant 2: if there are (k/2 + 2) buckets with the number of
 elements of the same size.
 EH_histogram example : [32,16,8,8,4,4,2,1,1]
 Return: [] if there no such elements, otherwise - the element's value.
"""
function check_condition(EH_histogram::Array{Int,1}, k::Int)
    uniq_vals = unique(EH_histogram)
    h::Array{Int,1} = zeros(maximum(uniq_vals))
    out_of_condition_values::Array{Int,1} = []
    for el in EH_histogram
        if el > 0
            h[el] += 1
            if h[el] >= convert(Int, k/2 + 2)
                push!(out_of_condition_values, el)
            end
        end
    end
    return sort(unique(out_of_condition_values))
end

# convert(Int, 4/2 + 2)

"""
 The state of the exponential histogram.
 Example:
 [3, 2, 1, 1, 1, 0, 0]
 head = 1
 tail = 5
 total = 8
 maxLen = 7
"""
mutable struct EHState
    head::Int
    tail::Int
    total::Int
    maxLen::Int
    function EHState(maxLen::Int)
        new(0,0,0,maxLen)
    end
end

function checkEHState()
    # for i = 1:n
    #     e = new_buckets[i]
    #     if e != 0
    #         # @assert e <= 2^(n-i) # todo check boundaries
    #     else
    #         break
    #     end
    # end
end

"Find tail of the buckets array"
function find_tail(buckets::Array{Int,1})
    tail::Int=-1
    flag::Bool=false
    n::Int = length(buckets)
    for k=1:n
        if flag
            break
        else
            if buckets[k] == 0
                tail = k-1
                flag = true
            end
            if k==n
                tail=n
                flag=true
            end
        end
    end
    tail
end
# find_tail([1,2,3,0,0,1])
# find_tail([1,2,3,5,5,1])

"Check buckets array structure - if it consists of consecutive
positive numbers and after that only consecutive zeros.
Or
- only zeros
- or all positive.
- exponentiality = MERGE TO ALREADY IMPLEMENTED FUNC
"
function check_buckets(buckets::Array{Int,1}) end

"Find first 2 same elements"
function find2same(buckets::Array{Int,1}, vl::Int)
    flag::Bool=false
    cand_idx::Int=0
    for k=2:length(buckets)
        if flag
            break
        else
            if buckets[k-1]==vl && buckets[k]==vl
                cand_idx = k
                flag=true
            end
        end
    end
    cand_idx
end

"""
 Fill from right to left.
 2: Translate one state to another.
 Test cases:
 [0,0,0,0,0,0,0,0,0,0], 1 -> [0,0,0,0,0,0,0,0,0,1]
 [0,0,0,0,0,0,0,0,0,1], 1 -> [0,0,0,0,0,0,0,0,1,1]
 [0,0,0,0,0,0,0,0,1,1], 1 -> [0,0,0,0,0,0,0,1,1,1]
 [0,1,1,1,1,1,1,1,1,1], 1 -> [1,1,1,1,1,1,1,1,1,1]
 [32,16,8,8,4,4,2,1,1], 1 -> [32,16,8,8,4,4,2,2,1]

 STORE 1'S IN THE BUCKETS OF SIZES OF POWER 2.
"""
function basic_counting2!(buckets::Array{Int,1},
                          ehstate::EHState,
                          x_new::Int;
                          k::Int=2)
    n = length(buckets)
    ## initialization
    if ehstate.head == 0 && ehstate.tail == 0
        buckets[end] = x_new
        ehstate.total += x_new
        ehstate.head = n
        ehstate.tail = n
    else
        if buckets[1] == 0 ## filling up
            buckets = circshift(buckets, -1)
            buckets[end] = x_new
            ehstate.total += x_new
            ehstate.head -= 1
            ehstate.tail = n
        else ## All buckets are not empty
            ehstate.total += x_new
            ##=== Start add new data
            if ehstate.tail < n
                @assert buckets[ehstate.tail+1] == 0
                @assert buckets[ehstate.tail] != 0
                buckets[ehstate.tail+1] = x_new
            else
                buckets[end] = x_new # add new data
            end
            ehstate.tail = find_tail(buckets)
            ##=== End add new data
            ch = check_condition(buckets, k)
            if ch == [] ## -> don't do collapse
                if ehstate.head==1 && ehstate.tail == n
                    buckets[1] = 0 # release memory
                    buckets = circshift(buckets, -1)
                end
                println("Don't do collapse:", buckets, " tail=", ehstate.tail)
            else # -> do collapse
                println("TOTAL (before collapse)=", ehstate.total)
                ##=== Start Find same succesive elements to collapse
                println("candidates indices (ch)=$ch")
                v :: Int = ch[1] ##??? fix - joining only 1-s
                cand_idx = find2same(buckets, v)
                ## === End Find same succesive elements to join
                if cand_idx != 0
                    println("Join ", cand_idx-1, " and ", cand_idx, " in ", buckets)
                    # println("Join ", buckets[cand_idx-1], " and ", buckets[cand_idx], " in ", buckets)
                    @assert buckets[cand_idx-1] == buckets[cand_idx]
                    buckets=collapse_buckets(buckets, val=buckets[cand_idx], num=2)
                    # -- sync circshift and collapse
                    ehstate.tail = find_tail(buckets)
                    println("After collapse: ", buckets)
                end
                ##=== End Find same succesive elements to collapse
            end
        end
    end
    buckets, ehstate
end

function test_check_condition()
    println("... Run test_check_condition()")
    @test check_condition([32,16,1,1,1], 2) == [1]
    @test check_condition([32,16,1,1,1], 4) == []
    @test check_condition([32,16,16,16,16,1,1,1], 4) == [16]
    println(".. Stop")
end

""" Sliding window demo. """
function sliding_window_demo()
    x = [1:100;]
    window=zeros(10)
    for i=1:length(x)
        window=circshift(window,1)
        window[1]=x[i]
        println(window)
    end
end

function test_basic_counting()
    println("... Run test_basic_counting()")
    r,s=basic_counting2!(zeros(Int,10), EHState(10), 1)
    for k=1:550
        println("k=$k")
        r, s = basic_counting2!(r, s, 1,k=4)
        # println(r, " total:", s.total, " head:", s.head, " tail:", s.tail)
    end
    # r = basic_counting!([32,16,8,8,4,4,4,2,1,1], 1)

    @test merge_buckets!([1,1,2,3], 1, 2) == [2,2,3,0]
    @test merge_buckets!([1,1], 1, 2) == [2,0]
    # @test basic_counting([1,1,2,4,4,8,8,16,32],1) == [32,16,8,8,4,4,2,2,1]
    # @test basic_counting([32,16,8,8,4,4,2,2,1],1) == [32,16,16,8,4,2,1]
    # @test check_condition([32,16,8,8,4,4,2,1,1] == true)
    # @test check_condition([32,16,8,8,4,4,2,2,1] == true)
    println(".. Stop ")
end

test_check_condition()
test_collapse_buckets()
test_basic_counting()

# tic()
# r = naiveCounting(1000000, 10000)
# toc()

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
"""FM Signal"""
function fmWave(P::Process)
    # Only parameter is the FM signal
    T = P.transient_t0:P.savedt:P.tmax
    p = parameter_functions(P)
    Δ𝑓 = 1.0
    sol = zeros(size(T))
    pint = 0.0
    for i ∈ 2:lastindex(T)
        t = T[i]
        pint += (p(t-P.savedt) + p(t))*P.savedt/2 # Crude integration, should be fine
        sol[i] = cos(2π*(t + Δ𝑓*pint))
    end
    return sol
end

fmWaveSim = Process(
    process = fmWave,
    parameter_profile = (stepNoise),
    parameter_profile_parameters = ((0.0, 100.0), 10.0, 0.5, 0.0),
    transient_t0 = 0.0,
    t0 = 0.0,
    savedt = 0.01,
    tmax = 100.0)
export fmWaveSim

"""Noisy Sine"""
function noisySine(P::Process)
    seed(P.solver_rng)
    sol = [sin(t) + parameter_function(P)(t)*randn() for t in P.transient_t0:P.savedt:P.tmax]
end

noisySineSim = Process(
    process = noisySine,
    X0 = [0.0],
    parameter_profile = unitStep,
    parameter_profile_parameters = (100.0, 0.0, 0.1), # (threshold, baseline, stepHeight)
    transient_t0 = -10.0,
    dt = 0.01,
    savedt = 0.01,
    tmax = 200.0,
    alg = nothing,
    solver_opts=Dict())
export noisySineSim

"""Noisy Shifty Scaly Sine"""
function noisyShiftyScalySine(P::Process)
    # Parameters (η, C, A)
    seed(P.solver_rng)
    (η, C, A) = parameter_functions(P)
    sol = [A(t)*sin(t) + η(t)*randn() + C(t) for t in P.transient_t0:P.savedt:P.tmax]
end

function noisyShiftyShcalySine(P::Process)
    # Parameters (η, C, A)
    seed(P.solver_rng)
    (η, C, A) = parameter_functions(P)
    sol = [A(t)*(sin(t) + η(t)*randn() + C(t)) for t in P.transient_t0:P.savedt:P.tmax]
end

"""Shifty Scaly Sine"""
function shcalySine(P::Process)
    seed(P.solver_rng)
    A = parameter_function(P)
    sol = [A(t)*(sin(t) + 1.0*randn() + 1.0) for t in P.transient_t0:P.savedt:P.tmax]
end

noisyShiftyScalySineSim = Process(
    process = noisyShiftyScalySine,
    X0 = [0.0],
    parameter_profile = (ramp, constantParameter, constantParameter), # (η, C, A)
    parameter_profile_parameters = ((0.001, 2.0, 0.0), (1.0,), (1.0,)),
    transient_t0 = 0.0,
    t0 = 0.0,
    dt = 0.1,
    savedt = 0.1,
    tmax = 1000.0,
    alg = nothing,
    solver_opts=Dict())
export noisyShiftyScalySineSim

noisyShiftyShcalySineSim = Process(
    process = noisyShiftyShcalySine,
    X0 = [0.0],
    parameter_profile = (constantParameter, constantParameter, constantParameter), # (η, C, A)
    parameter_profile_parameters = ((1.0,), (1.0,), (1.0,)),
    transient_t0 = 0.0,
    t0 = 0.0,
    dt = 0.05,
    savedt = 0.05,
    tmax = 500.0,
    alg = nothing,
    solver_opts=Dict())
export noisyShiftyShcalySineSim

shcalySineSim = Process(
    process = shcalySine,
    X0 = [0.0],
    parameter_profile = ramp, # (η, C, A)
    parameter_profile_parameters = (0.001, 1.0, 0.0),
    transient_t0 = 0.0,
    t0 = 0.0,
    dt = 0.01,
    savedt = 0.01,
    tmax = 1000.0,
    alg = nothing,
    solver_opts=Dict())
export shcalySineSim


text:
using FieldTracer, Meshes
using Test

@testset "FieldTracer.jl" begin
   @testset "Seeding" begin
      x = 0.0:0.1:1.0
      y = 1.0:0.1:2.0
      seeds = select_seeds(x, y; nsegment=(2,2))
      @test seeds[:,4] == [0.75, 1.75]
      z = 2.0:0.1:3.0
      seeds = select_seeds(x, y, z; nsegment=(2,2,2))
      @test seeds[:,2] == [0.75, 1.25, 2.25]
   end

   @testset "2D structured mesh" begin
      # bilinear interpolation function in a normalized rectangle
      x, y = 0.1, 0.2
      Q00, Q10, Q01, Q11 = 3.0, 40.0, 5.0, 60.0
      out = FieldTracer.bilin_reg(x, y, Q00, Q10, Q01, Q11)
      @test out ≈ 7.460 atol=1e-5

      # Euler2 and RK4 functions
      # https://se.mathworks.com/help/matlab/ref/streamline.html
      ds, maxstep = 0.1, 150
      x = range(0, 1, step=0.1)
      y = range(0, 1, step=0.1)
      # ndgrid
      xgrid = [i for i in x, _ in y]
      ygrid = [j for _ in x, j in y]
      u = copy(xgrid)
      v = -copy(ygrid)
      xstart = 0.1
      ystart = 0.9

      xt, yt = FieldTracer.euler(maxstep, ds, xstart, ystart, x, y, u, v)

      @test length(xt) == 141

      xt, yt = FieldTracer.rk4(maxstep, ds, xstart, ystart, x, y, u, v)

      @test length(xt) == 140

      grid = CartesianGrid((length(x)-1,length(y)-1),(0.,0.),(0.1,0.1))

      # default direction is both
      xt, yt = trace(u, v, xstart, ystart, grid; alg=Euler(), maxstep, ds)

      @test length(xt) == 148

      include("test_structured.jl")
      # asymptotic field
      @test test_trace_asymptote()     # single precision
      @test test_trace_asymptote(true) # double precision

      # dipole field
      @test test_trace_dipole()        # dipole tracing in 2D
   end

   @testset "3D structured mesh" begin
      # trilinear interpolation function in a normalized box
      x, y, z = 0.1, 0.2, 0.3
      Q = [3.0, 40.0, 5.0, 60.0, 3.0, 40.0, 5.0, 60.0]
      out = FieldTracer.trilin_reg(x, y, z, Q)
      @test out ≈ 7.460 atol=1e-5

      x = range(0, 10, length=15)
      y = range(0, 10, length=20)
      z = range(0, 10, length=25)
      bx = fill(1.0, length(x), length(y), length(z))
      by = fill(1.0, length(x), length(y), length(z))
      bz = fill(1.0, length(x), length(y), length(z))
      xs, ys, zs = 1.0, 1.0, 1.0

      # Euler 2nd order
      x1, y1, z1 = trace(bx, bz, bz, xs, ys, zs, x, y, z;
         alg=Euler(), ds=0.2, maxstep=200)

      @test length(x1) == 170 && x1[end] ≈ y1[end] ≈ z1[end]

      # RK4 by default
      x1, y1, z1 = trace(bx, bz, bz, xs, ys, zs, x, y, z;
         ds=0.2, maxstep=200, direction="forward")

      @test length(x1) == 152 && x1[end] ≈ y1[end] ≈ z1[end]

      Δx = x[2] - x[1]
      Δy = y[2] - y[1]
      Δz = z[2] - z[1]

      grid = CartesianGrid((length(x)-1, length(y)-1, length(z)-1),
         (0., 0., 0.),
         (Δx, Δy, Δz))

      # default direction is both
      x1, y1, z1 = trace(bx, bz, bz, xs, ys, zs, grid;
         alg=Euler(), ds=0.2, maxstep=200)

      @test length(x1) == 170 && x1[end] ≈ y1[end] ≈ z1[end]
   end

   @testset "2D unstructured mesh" begin
      include("test_unstructured.jl")
      @test test_trace_unstructured2D()
   end

end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdtreeselect

"""
    'feffery'_antdtreeselect(;kwargs...)

An AntdTreeSelect component.

Keyword arguments:
- `id` (String; optional)
- `allowClear` (Bool; optional)
- `bordered` (Bool; optional)
- `className` (String; optional)
- `defaultValue` (String | Real | Array of String | Reals; optional)
- `disabled` (Bool; optional)
- `listHeight` (Real; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `maxTagCount` (Real; optional)
- `multiple` (Bool; optional)
- `placeholder` (String; optional)
- `size` (String; optional)
- `style` (Dict; optional)
- `treeCheckStrictly` (Bool; optional)
- `treeCheckable` (Bool; optional)
- `treeData` (optional)
- `treeDefaultExpandAll` (Bool; optional)
- `treeDefaultExpandedKeys` (Array of Strings; optional)
- `treeExpandedKeys` (Array of Strings; optional)
- `value` (String | Real | Array of String | Reals; optional)
- `virtual` (Bool; optional)
"""
function 'feffery'_antdtreeselect(; kwargs...)
        available_props = Symbol[:id, :allowClear, :bordered, :className, :defaultValue, :disabled, :listHeight, :loading_state, :maxTagCount, :multiple, :placeholder, :size, :style, :treeCheckStrictly, :treeCheckable, :treeData, :treeDefaultExpandAll, :treeDefaultExpandedKeys, :treeExpandedKeys, :value, :virtual]
        wild_props = Symbol[]
        return Component("'feffery'_antdtreeselect", "AntdTreeSelect", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: Beddow/feffery-antd-components

max_issues_repo_name: Beddow/feffery-antd-components

max_forks_repo_name: Beddow/feffery-antd-components


text:
using Pkg
Pkg.activate("./GraphFusedElasticNet.jl/")
Pkg.instantiate()
using DelimitedFiles
using CSV
using DataFrames

workdir() = pwd()
splits_relpath() = "processed_data/splitlevels.csv"
vertexinfo_relpath() = "processed_data/vertex_data.csv"
betas_relpath() = "best_betas/"
betafiles = readdir(joinpath(workdir(), betas_relpath()))

# output settings
# optionally we can set a numer of files for the output
# so that file size is smaller, in this case, they can best_betas
# read in order and reassambled by horizontal concatenation
num_files = 20
smoothprobs_relpath() = "output_smooth_probs/"

# read the data from the graph
# actualy we only need the number of nodes
types=Dict(
    "nodes" => String,
    "vertex" => Int,
    "taz" => Int,
    "hours" => Int)
graphdata = CSV.read(vertexinfo_relpath(), types=types)
num_nodes = size(graphdata, 1)

function getsplit(s::String)
    x = split(s, '_')[3]
    x = split(x, '.')[1]
    return x
end

# read files with betas
betas = Dict{Tuple{Int, Int}, Vector{Float64}}()
for file in betafiles
    bin = getsplit(file)
    a, b = parse.(Int, String.(split(bin, '-')))
    β = vec(readdlm(joinpath(betas_relpath(), file)))
    betas[a, b] = β
end

splitlevels = vec(readdlm(splits_relpath()))
nbins = length(splitlevels) - 1

# descend in each bin and multiplyprobabilities
probs = ones(nbins, num_nodes)
for (key, value) in betas
    a = key[1] + 1
    b = key[2]
    m = (a + b) ÷ 2
    p = 1.0 ./ (1.0 .+ exp.(-value))
    p = reshape(p, 1, :)
    probs[a:m, :] .*= p
    probs[m:b, :] .*= 1.0 .- p
end
probs ./= sum(probs, dims=1)

# write output
if num_files == 1
    fname = joinpath(smoothprobs_relpath(), "probs.csv")
    writedlm(fname, probs, ',')
else
    chunksize = if num_nodes % num_files == 0
        num_nodes ÷ num_files
    else
        (num_nodes ÷ num_files) + 1
    end
    for i in 1:num_files
        from = (i - 1) * chunksize + 1
        to = min(i * chunksize, num_nodes)
        chunk = probs[:, from:to]
        fname = joinpath(smoothprobs_relpath(), "probs-chunk_$i.csv")
        writedlm(fname, chunk, ',')
    end
end

# probs to densities
# Δ = splitlevels[2:end] .- splitlevels[1:end-1]
# dens = probs ./ Δ

max_stars_repo_name: mauriciogtec/gfen-reproduce

max_issues_repo_name: mauriciogtec/gfen-reproduce

max_forks_repo_name: mauriciogtec/gfen-reproduce


text:
@testset "Format real numbers" begin
    for x in [0.1111111111111111111111111, 0.111111111f0, π, 111112.2234]
        for w in 5:15
            s_x = TableLogger.format_table_value(w, x)
            @test s_x isa String
            @test length(s_x) <= w
        end
    end

    x = 0.12345
    w = 3
    s_x = TableLogger.format_table_value(w, x)
    @test s_x == "0.1"

    x = 0.12345
    w = 5
    s_x = TableLogger.format_table_value(w, x)
    @test s_x == "0.12"

    x = 0.12345
    w = 5
    s_x = TableLogger.format_table_value(w, x; default_precision=5)
    @test s_x == "0.123"
end


text:
module YaoBlocksQobj

export convert_to_qobj, convert_to_qbir
export U1, U2, U3

using Configurations
using IBMQClient.Schema
using YaoBlocks

include("qobj.jl")
include("qbir.jl")

end


text:
################################################################################
## MOI wrapper
#################################################################################
## This file is largely inspired from MOI wrappers existing in other
## JuliaOpt packages:
## https://www.juliaopt.org/
## The authors are indebted to the developers of JuliaOpt for
## the current MOI wrapper.
##
## In particular, KNITRO.jl and IpOPT.jl where used as an inspiration

RESHOP_OPTIONS = ["solver"]

# TODO(Xhub) is this useful?
const MOIU = MathOptInterface.Utilities

# Scalar Functions
const SF = Union{MOI.ScalarAffineFunction{Float64},
				 MOI.ScalarQuadraticFunction{Float64}}

# Vector Types
const VAF = MOI.VectorAffineFunction{Float64}
const VOV = MOI.VectorOfVariables

const ALLV = Union{MOI.SingleVariable,MOI.VectorOfVariables}

# Scalar Sets
const SS = Union{MOI.EqualTo{Float64},
					  MOI.GreaterThan{Float64},
					  MOI.LessThan{Float64},
					  MOI.Interval{Float64}}
# LinSets
const LS = Union{MOI.EqualTo{Float64},
					  MOI.GreaterThan{Float64},
					  MOI.LessThan{Float64}}
# VecLinSets
const VLS = Union{MOI.Nonnegatives,
				  MOI.Nonpositives,
				  MOI.Zeros}

# NonConvex Sets
const NCS = Union{MOI.ZeroOne,
									MOI.Integer,
									MOI.Semicontinuous,
									MOI.Semiinteger,
									MOI.SOS1,
									MOI.SOS2}

# NonPolyhedral Cones
const NPC = Union{MOI.SecondOrderCone,
MOI.RotatedSecondOrderCone,
MOI.ExponentialCone,
MOI.DualExponentialCone,
MOI.PowerCone,
MOI.DualPowerCone}

mutable struct Optimizer <: MOI.AbstractOptimizer
    # Get number of solve for restart.
    number_solved::Int
    sense::MOI.OptimizationSense
    status::Cint

    vov_mapping::Dict{MOI.ConstraintIndex{VOV, <: Union{VLS, MOI.SOS1{Float64}, MOI.SOS2{Float64}}}, Vector{Cuint}}
    vaf_mapping::Dict{MOI.ConstraintIndex{VAF, <: VLS}, Vector{Cuint}}
    quadfn_mapping::Dict{MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, <: SS}, Cuint}
    sos_sets::Dict{MOI.ConstraintIndex{VOV, <: Union{MOI.SOS1{Float64}, MOI.SOS2{Float64}}}, Union{MOI.SOS1{Float64}, MOI.SOS2{Float64}}}
    fake_cons_name::Dict{MOI.ConstraintIndex, String}
    delvars::Set{MOI.VariableIndex}

    ctx::Union{Ptr{context}, Nothing}
    ctx_dest::Ptr{context}
    mdl::Ptr{reshop_model}
    mdl_solver::Ptr{reshop_model}
    options::Ptr{reshop_options}
    rhp_options::Dict{String, Union{Cdouble, Cint, Bool, Cstring, String}}
    gams_dir::String
    avar_cache::Union{Ptr{abstract_var}, Nothing}
    solver_name::String
    start_nl_cons::Int
    len_nl_cons::Cuint
end

function helper_options(ctx, options, reshop_opts::Ptr{reshop_options})
    solver_name = ""
    rhp_options = Dict{String, Union{Cdouble, Cint, Bool, Cstring, String}}()
    for (name, value) in options
        if string(name) == "solver"
            solver_name = value
        else
            res = reshop_option_set(reshop_opts, string(name), value)
            if (res != 0)
                 rhp_set_option(ctx, string(name), value)
            end
            rhp_options[string(name)] = value
        end
    end
    return (solver_name, rhp_options)
end

function Optimizer(;options...)
    reshop_set_printops(stdout)
    # Create ReSHOP context.
    ctx = ctx_alloc()

    # TODO this is quite a hack just for the "output" option.
    # Refactoring option in ReSHOP will enable us to move on
    reshop_opts = reshop_options_alloc()
    solver_name, rhp_options = helper_options(ctx, options, reshop_opts)

    model = Optimizer(0, MOI.FEASIBILITY_SENSE, 0,
                      Dict{MOI.ConstraintIndex{VOV, <: VLS}, Cuint}(), Dict{MOI.ConstraintIndex{VAF, <: VLS}, Cuint}(),
                      Dict{MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, <: SS}, Cuint}(),
                      Dict{MOI.ConstraintIndex{VOV, <: Union{MOI.SOS1{Float64}, MOI.SOS2{Float64}}}, Union{MOI.SOS1{Float64}, MOI.SOS2{Float64}}}(),
                      Dict{MOI.ConstraintIndex,String}(), Set{MOI.VariableIndex}(),
                      ctx, Ptr{context}(C_NULL), Ptr{reshop_model}(C_NULL), Ptr{reshop_model}(C_NULL),
                      reshop_opts, rhp_options, "", nothing, solver_name, -1, 0)

    finalizer(MOI.empty!, model)
    return model
end

# Print Optimizer.
function Base.show(io::IO, model::Optimizer)
    println(io, "A MathOptInterface model with backend:")
    println(io, model.ctx)
    return
end

# copy
MOIU.supports_default_copy_to(model::Optimizer, copy_names::Bool) = true
function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; kws...)
    return MOI.Utilities.automatic_copy_to(model, src; kws...)
end

# Some utilities.
number_variables(model::Optimizer) = ctx_numvar(model.ctx) - length(model.delvars)
number_constraints(model::Optimizer) = ctx_m(model.ctx)

# Getter for solver's name.
function MOI.get(model::Optimizer, ::MOI.SolverName)
    if model.solver_name == ""
        return "ReSHOP"
    else
        return "ReSHOP ($(model.solver_name))"
    end
end

# MOI.Silent.
MOI.supports(model::Optimizer, ::MOI.Silent) = true
function MOI.get(model::Optimizer, ::MOI.Silent)
    return reshop_option_get(model.options, "output") == 0xF
end

function MOI.set(model::Optimizer, ::MOI.Silent, value)
    val = value ? 0xF : 0x8
    reshop_option_set(model.options, "output", val)
    return
end

# MOI.TimeLimitSec.
MOI.supports(model::Optimizer, ::MOI.TimeLimitSec) = false

# MOI.RawParameters
function MOI.supports(model::Optimizer, param::MOI.RawParameter)
    name = param.name
    if name in RESHOP_OPTIONS
        return true
    end
    return false
end

function MOI.set(model::Optimizer, p::MOI.RawParameter, value)
    if !MOI.supports(model, p)
        throw(MOI.UnsupportedAttribute)
    end
    if p.name == "solver"
         model.solver_name = value
    end
    #reshop_option_set(model.options, p.name, value)
    return
end

function MOI.get(model::Optimizer, p::MOI.RawParameter)
    if haskey(model.options, p.name)
        return model.options[p.name]
    end
    error("RawParameter with name $(p.name) is not set.")
end

##################################################
# Optimize
##################################################
function MOI.optimize!(model::Optimizer)
    model.number_solved += 1
    # If gams_dir already existed (from a past solve), get rid of it
    # TODO(xhub) model.gams_dir could be a list of directory to remove
    if (!isempty(model.gams_dir))
        try
            rm(model.gams_dir, recursive=true, force=true)
        catch
            iswin && run(`cmd /C RMDIR /s /q $(model.gams_dir)`)
        end
    end
    # TODO check if gams_dir and ctx_dest already exists, do not reallocate then.
    model.ctx_dest, model.gams_dir = reshop_setup_gams()

    # Calling from emp, we already have a mdl object
    if model.mdl == C_NULL
        model.mdl = reshop_alloc(model.ctx)
    end
    model.mdl_solver = reshop_alloc(model.ctx_dest)
    model.status = reshop_solve(model.mdl, model.mdl_solver, model.ctx_dest, model.solver_name)
    reshop_postprocess(model.mdl_solver)

    return
end

function MOI.empty!(model::Optimizer)
    ctx_dealloc(model.ctx_dest)
    model.ctx_dest = C_NULL
    reshop_free(model.mdl)
    model.mdl = C_NULL
    reshop_free(model.mdl_solver)
    model.mdl_solver = C_NULL
    if model.ctx != nothing
        ctx_dealloc(model.ctx)
        model.ctx = ctx_alloc()
        helper_options(model.ctx, model.rhp_options, model.options)
    end

    model.number_solved = 0
    model.vov_mapping = Dict()
    model.vaf_mapping = Dict()
    model.quadfn_mapping = Dict()
    model.sos_sets       = Dict()
    model.fake_cons_name = Dict()
    model.delvars    = Set()
    # TODO(xhub) cleanup avar?
    # TODO(xhub) check that it is not necessary to do anything here
    #set_options(model, model.options)
    if (!isempty(model.gams_dir))
        try
            rm(model.gams_dir, recursive=true, force=true)
        catch
            iswin && run(`cmd /C RMDIR /s /q $(model.gams_dir)`)
        end
    end
    model.start_nl_cons = -1
    model.len_nl_cons = 0
    return
end

function MOI.is_empty(model::Optimizer)
    return model.number_solved == 0 &&
           length(model.vov_mapping) == 0 &&
           length(model.vaf_mapping) == 0 &&
           length(model.quadfn_mapping) == 0 &&
           length(model.sos_sets) == 0 &&
           ctx_m(model.ctx) == 0 &&
           ctx_numvar(model.ctx) == 0
end

include(joinpath("MOI_wrapper", "getters.jl"))
include(joinpath("MOI_wrapper", "supports.jl"))
include(joinpath("MOI_wrapper", "variables.jl"))
include(joinpath("MOI_wrapper", "constraints.jl"))
include(joinpath("MOI_wrapper", "objective.jl"))
include(joinpath("MOI_wrapper", "results.jl"))
include(joinpath("MOI_wrapper", "utils.jl"))
include(joinpath("MOI_wrapper", "changes.jl"))
include(joinpath("MOI_wrapper", "nlp.jl"))


text:
module FunctionTest1
   using FactCheck
   include(Pkg.dir("MethylUtils", "src","Functions.jl"))
   include(Pkg.dir("MethylUtils", "src","FeatureServices.jl"))
   file1 = Pkg.dir("MethylUtils", "testdata", "moabs-format-1000lines.txt")
   file2 = Pkg.dir("MethylUtils", "testdata", "ensGene1000lines.txt.gz")

   ##### READ TEST DATA #####

   ## cytosine counts ##
   cytosine_counts         =   load_moabs_format_as_sequence_features(file1)
   println("Cytosine Counts")
   println(cytosine_counts)

   ## annotation file ##
   ensgene_gene_regions_ft =   make_ensgene_gene_feature_table(file2, "mm10")
   ensgene_gene_regions    =   FeatureRegions(ensgene_gene_regions_ft)
   println(ensgene_gene_regions_ft)

   ####====####

   ###### METHYLATION ESTIMATION FUNCTONS #################

   cytosine_methylation    =   run_mc_hmc_calc( cytosine_counts )
   println("Cytosines Methylation")
   println(cytosine_methylation)
   println(typeof(cytosine_methylation))

   ####====####

   ####  TEST FILTER ####
   for key in keys(cytosine_methylation.features)
       vec = cytosine_methylation.features[key]
       filter!( x -> x.hmc != 2, vec)
   end

   println("Cytosine Methylation after filtering")
   println(cytosine_methylation)

   ###### TILING RELATED  ##############
   by_gene = cytosine_methylation_by_gene( cytosine_methylation, ensgene_gene_regions )
   println(by_gene)

   ####====####

   ###### UTILS ############

   # TEST INTERSECT BED
   println("Intersect $file1 and $file2")
   mm10_ensgene_path = Pkg.dir("MethylUtils", "testdata", "ensGene1000lines.bed")
   feature_table_to_bed(ensgene_gene_regions_ft, mm10_ensgene_path )
   path = intersect_files(file1, mm10_ensgene_path, overwrite=true)

   println("intersected files path: $path")

   ####====####

end

  # TEST MLML
  # test that we have order cytosines

  # include(Pkg.dir("MethylUtils", "src","FeatureServices.jl"))
  # println("start")
  # local_store_path =  Pkg.dir("MethylUtils", "testdata","store", "annotation")
  # SET_LOCAL_STORE_PATH(local_store_path)
  # moabs_bs_path = Pkg.dir("MethylUtils", "testdata","moabs-format-1000lines.txt")
  # organism_id = "mm10"
  # ensgene_bodies =  ensgene_regions( organism_id, "ensGene1000lines.txt.gz")
  # println(ensgene_bodies)
  # println(moabs_bs_path)
  # intersect_bed_path = intersect_bed_with_gene_regions( moabs_bs_path, ensgene_bodies)
  # println(intersect_bed_path)
  # # rm the file for testing
  # SET_LOCAL_STORE_PATH("nothing")

  #rm files for testing next time
  # Pkg.dir("MethylUtils", "testdata", "ensGene1000lines.bed")
  # Pkg.dir("MethylUtils", "testdata", "moabs-format-1000lines.txt.intersect.bed")