text:
abstract type AbstractFiniteDifferenceSpace <: AbstractSpace end

abstract type Staggering end

""" Cell center location """
struct CellCenter <: Staggering end

""" Cell face location """
struct CellFace <: Staggering end

struct FiniteDifferenceSpace{S <: Staggering, M <: Meshes.IntervalMesh, G} <:
       AbstractFiniteDifferenceSpace
    staggering::S
    mesh::M
    center_local_geometry::G
    face_local_geometry::G
end

function FiniteDifferenceSpace{S}(
    mesh::Meshes.IntervalMesh,
) where {S <: Staggering}
    CT = Meshes.coordinate_type(mesh)
    AIdx = Geometry.coordinate_axis(CT)
    # TODO: FD operators  hardcoded to work over the 3-axis, need to generalize
    # similar to spectral operators
    @assert AIdx == (3,) "FiniteDifference operations only work over the 3-axis (ZPoint) domain"
    FT = eltype(CT)
    face_coordinates = collect(mesh.faces)
    Mxξ = Geometry.Axis2Tensor{
        FT,
        Tuple{Geometry.CartesianAxis{AIdx}, Geometry.CovariantAxis{AIdx}},
        SMatrix{1, 1, FT, 1},
    }
    Mξx = Geometry.Axis2Tensor{
        FT,
        Tuple{Geometry.ContravariantAxis{AIdx}, Geometry.CartesianAxis{AIdx}},
        SMatrix{1, 1, FT, 1},
    }
    LG = Geometry.LocalGeometry{CT, FT, Mxξ, Mξx}
    nface = length(face_coordinates)
    ncent = nface - 1
    center_local_geometry = DataLayouts.VF{LG}(Array{FT}, ncent)
    face_local_geometry = DataLayouts.VF{LG}(Array{FT}, nface)
    for i in 1:ncent
        # centers
        coord⁻ = Geometry.component(face_coordinates[i], 1)
        coord⁺ = Geometry.component(face_coordinates[i + 1], 1)
        # at the moment we use a "discrete Jacobian"
        # ideally we should use the continuous quantity via the derivative of the warp function
        # could we just define this then as deriv on the mesh element coordinates?
        coord = (coord⁺ + coord⁻) / 2
        Δcoord = coord⁺ - coord⁻
        J = Δcoord
        WJ = Δcoord
        ∂x∂ξ = SMatrix{1, 1}(J)
        ∂ξ∂x = SMatrix{1, 1}(inv(J))
        center_local_geometry[i] = Geometry.LocalGeometry(
            CT(coord),
            J,
            WJ,
            Geometry.AxisTensor(
                (
                    Geometry.CartesianAxis{AIdx}(),
                    Geometry.CovariantAxis{AIdx}(),
                ),
                ∂x∂ξ,
            ),
            Geometry.AxisTensor(
                (
                    Geometry.ContravariantAxis{AIdx}(),
                    Geometry.CartesianAxis{AIdx}(),
                ),
                ∂ξ∂x,
            ),
        )
    end
    for i in 1:nface
        coord = Geometry.component(face_coordinates[i], 1)
        if i == 1
            # bottom face
            coord⁺ = Geometry.component(face_coordinates[2], 1)
            J = coord⁺ - coord
            WJ = J / 2
        elseif i == nface
            # top face
            coord⁻ = Geometry.component(face_coordinates[i - 1], 1)
            J = coord - coord⁻
            WJ = J / 2
        else
            coord⁺ = Geometry.component(face_coordinates[i + 1], 1)
            coord⁻ = Geometry.component(face_coordinates[i - 1], 1)
            J = (coord⁺ - coord⁻) / 2
            WJ = J
        end
        ∂x∂ξ = SMatrix{1, 1}(J)
        ∂ξ∂x = SMatrix{1, 1}(inv(J))
        face_local_geometry[i] = Geometry.LocalGeometry(
            CT(coord),
            J,
            WJ,
            Geometry.AxisTensor(
                (
                    Geometry.CartesianAxis{AIdx}(),
                    Geometry.CovariantAxis{AIdx}(),
                ),
                ∂x∂ξ,
            ),
            Geometry.AxisTensor(
                (
                    Geometry.ContravariantAxis{AIdx}(),
                    Geometry.CartesianAxis{AIdx}(),
                ),
                ∂ξ∂x,
            ),
        )
    end
    return FiniteDifferenceSpace(
        S(),
        mesh,
        center_local_geometry,
        face_local_geometry,
    )
end

const CenterFiniteDifferenceSpace = FiniteDifferenceSpace{CellCenter}
const FaceFiniteDifferenceSpace = FiniteDifferenceSpace{CellFace}

function FiniteDifferenceSpace{S}(
    space::FiniteDifferenceSpace,
) where {S <: Staggering}
    FiniteDifferenceSpace(
        S(),
        space.mesh,
        space.center_local_geometry,
        space.face_local_geometry,
    )
end

Base.length(space::FiniteDifferenceSpace) = length(coordinates_data(space))

nlevels(space::FiniteDifferenceSpace) = length(space)

local_geometry_data(space::CenterFiniteDifferenceSpace) =
    space.center_local_geometry

local_geometry_data(space::FaceFiniteDifferenceSpace) =
    space.face_local_geometry

left_boundary_name(space::FiniteDifferenceSpace) =
    propertynames(space.mesh.boundaries)[1]

right_boundary_name(space::FiniteDifferenceSpace) =
    propertynames(space.mesh.boundaries)[2]


text:
# import log functions
@ms include("$APPDIR/utils/log.jl")

# read command line arguments
@ms include("$APPDIR/src/args.jl")
@ms global args = parse_commandline()

# import all routines
@ms include("$APPDIR/src/jump_model.jl")
@ms include("$APPDIR/src/instance.jl")
@ms include("$APPDIR/src/solution.jl")
@ms include("$APPDIR/src/solver.jl")

"""
Main function
"""
function main()

    # print active parameters
    display_active_args()

    ln0("\n=================== Start MAIN() =================\n")
    # clock time at the beginning (in s)
    t_begin = time()

    inst = Instance()
    sv = Solver(inst)
    solve!(sv)
    pprint(sv.sol)
    export_to_file(sv.sol)

    # total elapsed time of the MAIN
    sec = round((time() - t_begin), digits=3)
    ln0("\nTotal elapsed time MAIN():\t$(sec)s")
    ln0("\n==================== End MAIN() ==================\n")

end

ext: jl

max_stars_repo_name: antcrepin/flight-route-optimizer

max_issues_repo_name: antcrepin/flight-route-optimizer

max_forks_repo_name: antcrepin/flight-route-optimizer


text:
using LightGraphs, MetaGraphs, Xpress,  LinearAlgebra, SparseArrays, SuiteSparse

function load_matpower(filename)
  instance_name = split(filename, '.')[1]

  touch(instance_name*".temp")
  f = open(filename)
  out = open(instance_name*".temp", "w")

  # removing all ';' at end of lines
  while !eof(f)
    line = readline(f)
    if length(line) > 0 && line[1] != '%' && line[1] != 'f'
      s = split(line, ";")
      println(out, s[1])
    end
  end
  close(f)
  close(out)

  data = DelimitedFiles.readdlm(instance_name*".temp")
  rm(instance_name*".temp")
  return data
end

function find_numarray(i_start, data)
  i_debut = i_start
  while !isa(data[i_debut, 1], Int)
    i_debut+=1
  end
  i_fin=i_debut
  while !isa(data[i_fin,1], SubString)
    i_fin += 1
  end
  i_debut, i_fin-1
end

checkfor(data, line_ind, name) = (data[line_ind, 1] == name) || error("Expected ", name, " at line ", line_ind, ", got ", data[line_ind,1], " instead.")

function read_network_graph(instance_path::String)

    data = load_matpower(instance_path)

    ## Bus load and shunt information
    i_debut, i_fin = find_numarray(1, data)
    checkfor(data, i_debut-1, "mpc.bus")
    nb_bus = i_fin-i_debut+1
    index_bus = Dict( data[i+i_debut-1,1] => i for i=1:nb_bus)
    if index_bus != Dict( i => i for i=1:nb_bus)
        println("!!! Specific numerotation of buses !!! \n")
    end
    ## Bus generator information
    i_debut, i_fin = find_numarray(i_fin+1, data)
    checkfor(data, i_debut-1, "mpc.gen")

    #initialize network graph G
    G = MetaGraph(nb_bus)
    for i in 1:nb_bus
        set_props!(G, i, Dict(:name => "node$i"))
    end

    ## Link information
    i_debut, i_fin = find_numarray(i_fin+1, data)
    checkfor(data, i_debut-1, "mpc.branch")
    for i=i_debut:i_fin
      if data[i, 11] == 0
        #@warn("link $(data[i,1])⟶$(data[i,2]) breaker out of service !")
      else
          orig = index_bus[data[i,1]]
          dest = index_bus[data[i,2]]
          LightGraphs.add_edge!(G, orig, dest)
      end
    end
    return G
end

function read_sparsity_pattern(instance_path::String)

    data = load_matpower(instance_path)

    ## Bus load and shunt information
    i_debut, i_fin = find_numarray(1, data)
    checkfor(data, i_debut-1, "mpc.bus")
    nb_bus = i_fin-i_debut+1
    index_bus = Dict( data[i+i_debut-1,1] => i for i=1:nb_bus)

    ## Bus generator information
    i_debut, i_fin = find_numarray(i_fin+1, data)
    checkfor(data, i_debut-1, "mpc.gen")

    #initialize network graph G
    sp = spzeros(nb_bus,nb_bus)

    ## Link information
    i_debut, i_fin = find_numarray(i_fin+1, data)
    checkfor(data, i_debut-1, "mpc.branch")
    for i=i_debut:i_fin
      if data[i, 11] == 0
        #@warn("link $(data[i,1])⟶$(data[i,2]) breaker out of service !")
      else
          orig = index_bus[data[i,1]]
          dest = index_bus[data[i,2]]
          sp[orig,dest] = 1
      end
    end

    sp_sym = sp + sp'

    diag = zeros(nb_bus)

    for i=1:nb_bus
        sum_col = sum(sp_sym[i,:])
        diag[i] = Int(sum_col + 1)
    end
    return sp_sym + Diagonal(diag)
    # return sp+sp'+nb_bus*sparse(I, nb_bus, nb_bus)
end

function read_sparsity_pattern_real(instance_path::String)
    data = load_matpower(instance_path)
    ## Bus load and shunt information
    i_debut, i_fin = find_numarray(1, data)
    checkfor(data, i_debut-1, "mpc.bus")
    nb_bus = i_fin-i_debut+1
    index_bus = Dict( data[i+i_debut-1,1] => i for i=1:nb_bus)

    ## Bus generator information
    i_debut, i_fin = find_numarray(i_fin+1, data)
    checkfor(data, i_debut-1, "mpc.gen")

    sp = spzeros(2*nb_bus,2*nb_bus)

    ## Link information
    i_debut, i_fin = find_numarray(i_fin+1, data)
    checkfor(data, i_debut-1, "mpc.branch")
    for i=i_debut:i_fin
      if data[i, 11] == 0
        #@warn("link $(data[i,1])⟶$(data[i,2]) breaker out of service !")
      else
          orig = index_bus[data[i,1]]
          dest = index_bus[data[i,2]]
          node_orig_Re = orig
          node_orig_Im = nb_bus + orig
          node_dest_Re = dest
          node_dest_Im = nb_bus + dest               #bus i => node i for iRe, node nb_bus+i for iIm
          sp[node_orig_Re, node_dest_Re] = 1
          sp[node_orig_Re, node_dest_Im] = 1
          sp[node_orig_Im, node_dest_Re] = 1
          sp[node_orig_Im, node_dest_Im] = 1
      end
    end
    sp_sym = sp + sp'
    diag = zeros(2*nb_bus)
    for i=1:2*nb_bus
        sum_col = sum(sp_sym[i,:])
        diag[i] = Int(sum_col + 1)
    end
    return sp_sym + Diagonal(diag)
end

function chordal_ext_cholesky(sparsity_pattern)
    A = sparsity_pattern
    nb_edges_A = (nnz(A) - size(A,1))/2
    #computing cholesky factorisation of A
    F = cholesky(A) #NOTE: order = F.p
    # println(F.p)
    # computing L + LT
    L = sparse(F.L)
    nb_edges_L = nnz(L) - size(A,1)
    nb_added_edges = nb_edges_L - nb_edges_A
    SP = L + L'
    #inverting permutation to get chordal extension of sparsity_pattern
    H = SP[invperm(F.p), invperm(F.p)]
    return H, F.p, nb_added_edges
end

function construct_graph_from_matrix(L)
    n = size(L,1)
    H = MetaGraph(n)
    for i in 1:n
        set_props!(H, i, Dict(:name => "node$i"))
    end
    for i in 1:n
        for j in 1:i
            if L[i,j] != 0
                MetaGraphs.add_edge!(H,i,j)
            end
        end
    end
    return H
end

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

function MD_chordal_ext(G)
    n = nv(G)
    G_α = copy(G)
    H = copy(G)
    order = []
    nb_added_edges = 0
     for i=1:n
         degrees = [length(neighbors(G_α,s)) for s=1:(n-i+1)]
         perm = sortperm(degrees)
         v = perm[1]
         push!(order, Meta.parse(get_prop(G_α, v,:name)[5:end]))
         # println("v = ",v)
         # println(get_prop(G_α, v,:name))
         Nv = neighbors(G_α,v)
         for w in Nv
             for u in Nv
                 if w < u
                     if !has_edge(G_α, w, u)
                         add_edge!(G_α, w, u)
                         nodew = Meta.parse(get_prop(G_α, w,:name)[5:end])
                         nodeu = Meta.parse(get_prop(G_α, u,:name)[5:end])
                         add_edge!(H, nodew, nodeu)
                         nb_added_edges += 1
                    end
                     # println("edge = ",get_prop(G_α, w,:name),"--",get_prop(G_α, u,:name))
                end
            end
        end
        rem_vertex!(G_α,v)
    end
    return H, order, nb_added_edges
end

##########################################################################################################
function cliques_max(G, order)
    n = nv(G)
    C0 = []
    C = C0
    cliques_dict = Dict{Int64,Array{Any,1}}()
    k = 0
    for i=1:n
        v = order[i]
        Nv = neighbors(G,v)
        # println(v)
        C = [Nv[j] for j=1:length(Nv) if Nv[j] ∉ order[1:i]]
        C = [C ; v]
        # println(C)
        if C ⊈ C0
            k +=1
            cliques_dict[k] = C
            C0 = C
        end
    end
    return cliques_dict
end

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

function weighted_graph(cliques_dict)
    nb_cl = length(cliques_dict)
    G = MetaGraph(nb_cl)
    for i in 1:nb_cl
        for j in 1:i
            inter = length(intersect(cliques_dict[i], cliques_dict[j]))
            if inter > 0
                LightGraphs.add_edge!(G,i,j)
                set_prop!(G, i, j, :weight, inter)
            end
        end
    end
    return G
end

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

function Prim_algo(G)
    A = MetaGraph(nv(G))
    for i = 1:nv(G)
        set_prop!(A,i, :name, "C$i")
    end
    V = Set([ v for v in LightGraphs.vertices(G)])
    Vprim = Set([1])
    while Vprim != V
        max_weight = 0
        edge_to_add = LightGraphs.Edge(0,0)
        node_to_add = 0
        for edge in LightGraphs.edges(G)
            u = src(edge)
            v = dst(edge)
            if (u ∈ Vprim && v ∉ Vprim)
                weight_uv = get_prop(G,LightGraphs.Edge(u,v), :weight)
                if weight_uv > max_weight
                    max_weight = weight_uv
                    edge_to_add = LightGraphs.Edge(u,v)
                    node_to_add = v
                end
            elseif (u ∉ Vprim && v ∈ Vprim)
                weight_uv = get_prop(G,LightGraphs.Edge(u,v), :weight)
                if weight_uv > max_weight
                    max_weight = weight_uv
                    edge_to_add = LightGraphs.Edge(u,v)
                    node_to_add = u
                end
            end
        end
        LightGraphs.add_edge!(A, edge_to_add)
        set_prop!(A, edge_to_add, :weight, max_weight)
        push!(Vprim, node_to_add)
    end
    return A
end

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

function mergingIP(cliques_dict, clique_tree,size_max)
    m = Model(with_optimizer(Xpress.Optimizer))
    x = Dict((Meta.parse(get_prop(clique_tree,src(e), :name)[2:end]),Meta.parse(get_prop(clique_tree,dst(e), :name)[2:end])) => @variable(m, base_name = "x_$(Meta.parse(get_prop(clique_tree,src(e), :name)[2:end]))_$(Meta.parse(get_prop(clique_tree,dst(e), :name)[2:end]))", binary=true) for e in edges(clique_tree))
    edge = collect(edges(clique_tree))[1]
    exp = Dict(i => @expression(m, 0*x[(Meta.parse(get_prop(clique_tree,src(edge), :name)[2:end]),Meta.parse(get_prop(clique_tree,dst(edge), :name)[2:end]))]) for i in keys(cliques_dict))
    for e in edges(clique_tree)
        C1 = Meta.parse(get_prop(clique_tree,src(e), :name)[2:end])
        C2 = Meta.parse(get_prop(clique_tree,dst(e), :name)[2:end])
        #x[(C1,C2)] = @variable(m, "x_$(C1)_$(C2)", Bin)
        size_merge_clique = length(union(cliques_dict[C1], cliques_dict[C2]))
        @constraint(m, size_merge_clique*x[(C1,C2)] <= size_max)
        add_to_expression!(exp[C1], 1, x[(C1,C2)])
        add_to_expression!(exp[C2], 1, x[(C1,C2)])
    end

    for (Ci, ex_ctr) in exp
        @constraint(m, exp[Ci]<=1)
    end
    @objective(m, Max, sum(length(intersect(cliques_dict[Meta.parse(get_prop(clique_tree,src(e), :name)[2:end])],cliques_dict[Meta.parse(get_prop(clique_tree,dst(e), :name)[2:end])]))*(length(intersect(cliques_dict[Meta.parse(get_prop(clique_tree,src(e), :name)[2:end])],cliques_dict[Meta.parse(get_prop(clique_tree,dst(e), :name)[2:end])]))+1)/2*x[(Meta.parse(get_prop(clique_tree,src(e), :name)[2:end]),Meta.parse(get_prop(clique_tree,dst(e), :name)[2:end]))] for e in edges(clique_tree)))
    optimize!(m)
    new_cliques_dict = copy(cliques_dict)
    new_clique_tree = copy(clique_tree)
    for (cliques,var) in x
        if value.(x[cliques]) == 1
            list_nodes1 = cliques_dict[cliques[1]]
            list_nodes2 = cliques_dict[cliques[2]]
            list_nodes_merge = union(list_nodes1, list_nodes2)
            new_cliques_dict[cliques[1]] = list_nodes_merge
            delete!(new_cliques_dict, cliques[2])
            v_iter = MetaGraphs.filter_vertices(new_clique_tree, :name, "C$(cliques[1])")
            v = 0
            for i in v_iter
                v = i
            end
            w_iter = MetaGraphs.filter_vertices(new_clique_tree, :name, "C$(cliques[2])")
            w = 0
            for j in w_iter
                w = j
            end
            N = neighbors(new_clique_tree, w)
            for n in N
                if n != v
                    LightGraphs.add_edge!(new_clique_tree, v, n)
                    weight = get_prop(new_clique_tree, LightGraphs.Edge(w,n), :weight)
                    set_prop!(new_clique_tree, v, n, :weight, weight)
                    # println("add_edge", get_prop(new_clique_tree,v, :name), get_prop(new_clique_tree,n, :name))
                end
            end
            rem_vertex!(new_clique_tree, w)
        end
    end
    new_list_maximal_cliques_clean = Dict{Int64, Array{Any,1}}()
    compteur = 1
    Dict_new_numbers = Dict{Int64, Int64}()
    for (clique, list_nodes) in new_cliques_dict
            new_list_maximal_cliques_clean[compteur] = list_nodes
            Dict_new_numbers[clique] = compteur
            compteur += 1
    end
     for i in 1:nv(new_clique_tree)
         old_name = get_prop(new_clique_tree, i, :name)
         new_number = Dict_new_numbers[parse(Int64,old_name[2:end])]
         set_prop!(new_clique_tree, i, :name, "C$new_number")
     end
     return new_list_maximal_cliques_clean, new_clique_tree
end

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

function Molzahn_merging(cliques_dict, clique_tree, number_max_of_cliques)
    new_cliques_dict = copy(cliques_dict)
    new_clique_tree = copy(clique_tree)
    costs_fusion = Dict( e => cost_fusion(e, clique_tree, cliques_dict) for e in edges(clique_tree))
    while length(new_cliques_dict) > number_max_of_cliques
        (min_cost, edge_to_merge) = findmin(costs_fusion)
        clique1 = Meta.parse(get_prop(new_clique_tree, src(edge_to_merge), :name)[2:end])
        clique2 = Meta.parse(get_prop(new_clique_tree, dst(edge_to_merge), :name)[2:end])
        # println("fusion cliques $clique1 $clique2")
        list_nodes1 = new_cliques_dict[clique1]
        list_nodes2 = new_cliques_dict[clique2]
        list_nodes_fusion = union(list_nodes1,list_nodes2)
        # println("merged nodes : $list_nodes_fusion ")
        new_cliques_dict[clique1] = list_nodes_fusion
        delete!(new_cliques_dict, clique2)
        v_iter = MetaGraphs.filter_vertices(new_clique_tree, :name, "C$(clique1)")
        v = 0
        for i in v_iter
            v = i
        end
        w_iter = MetaGraphs.filter_vertices(new_clique_tree, :name, "C$(clique2)")
        w = 0
        for j in w_iter
            w = j
        end
        # println(w)
        N = neighbors(new_clique_tree, w)
        for n in N
            if n != v
                LightGraphs.add_edge!(new_clique_tree, v, n)
                weight = get_prop(new_clique_tree, LightGraphs.Edge(w,n), :weight)
                set_prop!(new_clique_tree, v, n, :weight, weight)
                # println("add_edge", get_prop(new_clique_tree,v, :name), get_prop(new_clique_tree,n, :name))
            end
        end
        rem_vertex!(new_clique_tree, w)
        costs_fusion = Dict(e => cost_fusion(e, new_clique_tree, new_cliques_dict) for e in edges(new_clique_tree))
    end
    new_list_maximal_cliques_clean = Dict{Int64, Array{Any,1}}()
    compteur = 1
    for (clique, list_nodes) in new_cliques_dict
            new_list_maximal_cliques_clean[compteur] = list_nodes
            compteur += 1
    end
    return new_list_maximal_cliques_clean
end

function cost_fusion(e, clique_tree, cliques_dict)
    clique1 = Meta.parse(get_prop(clique_tree, src(e), :name)[2:end])
    clique2 = Meta.parse(get_prop(clique_tree, dst(e), :name)[2:end])
    list_nodes1 = cliques_dict[clique1]
    list_nodes2 = cliques_dict[clique2]
    nintersect = length(intersect(list_nodes1,list_nodes2))
    elimmaxcliques1 = length(list_nodes1);
    elimmaxcliques2 = length(list_nodes2);
    lnewmaxcliques = elimmaxcliques1 + elimmaxcliques2 - nintersect
    nvarafter = (lnewmaxcliques)*(2*lnewmaxcliques+1) - (elimmaxcliques1*(2*elimmaxcliques1+1)+elimmaxcliques2*(2*elimmaxcliques2+1))
    ocostbefore = (nintersect)*(2*nintersect+1)
    cost = nvarafter - ocostbefore
    return cost
end

##########################################################################################################
#OUR DECOMPOSITION
##########################################################################################################
function cost_IPiteration(i, array_edges, cliques_dict, alpha, beta, m, nlc_before_merge)
    clique1 = array_edges[i,1]
    clique2 = array_edges[i,2]
    list_nodes1 = cliques_dict[clique1]
    # println(list_nodes1)
    list_nodes2 = cliques_dict[clique2]
    # println(list_nodes2)
    size_intersect = length(intersect(list_nodes1,list_nodes2))
    size_clique1 = length(list_nodes1)
    size_clique2 = length(list_nodes2)
    size_fusion_clique = size_clique1 + size_clique2 - size_intersect
    # println("Size clique 1 : $size_clique1, Size clique 2 : $size_clique2")
    #q1before_merge = sum(length(list_nodes)^3 for (clique, list_nodes) in cliques_dict) NOTE: inutile d'utiliser la somme totale car c'est la même pour tout le monde
    #q1 = q1before_merge - size_clique1^3 - size_clique2^3 + size_fusion_clique^3
    nlc = nlc_before_merge - size_intersect*(2*size_intersect+1)
    # println("nlc : $nlc")
    q2 = (m+nlc)^3
    cost = alpha*(-size_clique1^3-size_clique2^3+size_fusion_clique^3) + beta*q2
    return cost
end

function cost_Molzahn(i, array_edges, cliques_dict)
    clique1 = array_edges[i,1]
    clique2 = array_edges[i,2]
    list_nodes1 = cliques_dict[clique1]
    list_nodes2 = cliques_dict[clique2]
    nintersect = length(intersect(list_nodes1,list_nodes2))
    elimmaxcliques1 = length(list_nodes1);
    elimmaxcliques2 = length(list_nodes2);
    lnewmaxcliques = elimmaxcliques1 + elimmaxcliques2 - nintersect
    nvarafter = (lnewmaxcliques)*(2*lnewmaxcliques+1) - (elimmaxcliques1*(2*elimmaxcliques1+1)+elimmaxcliques2*(2*elimmaxcliques2+1))
    ocostbefore = (nintersect)*(2*nintersect+1)
    cost = nvarafter - ocostbefore
    return cost
end

function convert_clique_tree_to_array(A)
    array_edges = Array{Int64}(undef, ne(A),2)
    i = 0
    for e in edges(A)
        i += 1
        clique1 = Meta.parse(get_prop(A, src(e), :name)[2:end])
        clique2 = Meta.parse(get_prop(A, dst(e), :name)[2:end])
        array_edges[i,1] = clique1
        array_edges[i,2] = clique2
    end
    return array_edges
end

function fusion_Mosek_kmax_iterations(cliques_dict, clique_tree_array, number_max_of_cliques, alpha, beta, m, nlc, kmax)
    new_cliques_dict = copy(cliques_dict)
    array_edges = copy(clique_tree_array)
    nlc_before_merge = nlc
    costs_fusion = [cost_IPiteration(i, array_edges, new_cliques_dict, alpha, beta, m, nlc_before_merge) for i in 1:size(array_edges,1)]
    # println(costs_fusion)
    println("begin while...")
    it=0
    while length(new_cliques_dict) > number_max_of_cliques
        it+=1
        (min_cost, edge_to_merge) = findmin(costs_fusion)
        # println("Min cost : ", min_cost)
        # println("Edge to merge : ", edge_to_merge)
        clique1 = array_edges[edge_to_merge,1]
        clique2 = array_edges[edge_to_merge,2]
        # println("fusion cliques $clique1 $clique2")
        list_nodes1 = new_cliques_dict[clique1]
        list_nodes2 = new_cliques_dict[clique2]
        list_nodes_fusion = union(list_nodes1,list_nodes2)
        nintersect = length(intersect(list_nodes1,list_nodes2))
        # println("merged nodes : $list_nodes_fusion ")
        new_cliques_dict[clique1] = list_nodes_fusion
        delete!(new_cliques_dict, clique2)
        #delete edge in the clique tree array
        array_edges = array_edges[setdiff(1:end, edge_to_merge), :]
        #rename C2 into C1
        for i in 1:size(array_edges,1)
            if array_edges[i,1] == clique2
                array_edges[i,1] = clique1
            elseif array_edges[i,2] == clique2
                array_edges[i,2] = clique1
            end
        end
        # println(array_edges)
        nlc_before_merge -= (nintersect)*(2*nintersect+1)
        if it <= kmax
        # println("calculation new costs...")
            costs_fusion = [cost_IPiteration(i, array_edges, new_cliques_dict, alpha, beta, m, nlc_before_merge) for i in 1:size(array_edges,1)]
        # println(costs_fusion)
        else
            costs_fusion = [cost_Molzahn(i, array_edges, new_cliques_dict) for i in 1:size(array_edges,1)]
        end
    end
    new_list_maximal_cliques_clean = Dict{Int64, Array{Any,1}}()
    compteur = 1
    for (clique, list_nodes) in new_cliques_dict
            new_list_maximal_cliques_clean[compteur] = list_nodes
            compteur += 1
    end
    return new_list_maximal_cliques_clean
end

##########################################################################################################
#WRITE DECOMPOSITIONS
##########################################################################################################
function get_nb_linking_constraints(cliquetree)
    nb = 0
    for e in LightGraphs.edges(cliquetree)
        weight = get_prop(cliquetree, e, :weight)
        nb += weight*(2*weight+1)
    end
    return nb
end

function write_file_decompositions(data_path, instances, algo, repo_name, FUSION, fusion_type = "", size_max_IPmerging=0, nb_times=0, kmax=0)
    for instance in instances
        println(instance)
        instance_path = joinpath(data_path, instance)
        if algo == "MD_real"
            G = read_network_graph_real_numbers(instance_path)
        elseif algo == "cholesky" || algo == "Metis"
            sp = read_sparsity_pattern(instance_path)
        elseif algo == "cholesky_real"
            sp = read_sparsity_pattern_real(instance_path)
        else
            G = read_network_graph(instance_path)
        end
        if algo == "MD"
            H, order, nb_added_edges = MD_chordal_ext(G)
        elseif algo == "cholesky" || algo == "cholesky_real"
            L, order, nb_added_edges = chordal_ext_cholesky(sp)
            H = construct_graph_from_matrix(L)
        else
            println("error : algo $algo does not exist")
            exit()
        end
        nb_nodes = 0
        if algo == "cholesky" || algo == "cholesky_real"
            nb_nodes = size(sp,1)
        else
            nb_nodes = nv(G)
        end
        cliques_dict = cliques_max(H,order)
        GC = weighted_graph(cliques_dict)
        A = Prim_algo(GC)

        clique_tree = A
        if FUSION == "YES"
            if fusion_type == "IP"
                for i in 1:nb_times
                    new_cliques_dict, new_clique_tree = mergingIP(cliques_dict, clique_tree, size_max_IPmerging)
                    cliques_dict = new_cliques_dict
                    clique_tree = new_clique_tree
                end
            elseif fusion_type == "Molzahn"
                nb_max = 0.1*nb_nodes
                new_cliques_dict = Molzahn_merging(cliques_dict,A, nb_max)
                GC2 = weighted_graph(new_cliques_dict)
                clique_tree = Prim_algo(GC2)
                cliques_dict = new_cliques_dict
            elseif fusion_type == "our_merging_heuristic"
                nb_max = 0.1*nb_nodes
                nlc = get_nb_linking_constraints(A)
                infos_data_linear_regression = readdlm("..\\data\\parameters_linear_regression\\data_objscale10-4_$(instance[1:end-2]).csv", ';')
                m = infos_data_linear_regression[2,2]
                println(m)
                alpha = infos_data_linear_regression[13,1]
                println(alpha)
                beta = infos_data_linear_regression[13,2]
                println(beta)
                array_edges = convert_clique_tree_to_array(A)
                new_cliques_dict = fusion_Mosek_kmax_iterations(cliques_dict,array_edges, nb_max, alpha, beta, m, nlc, kmax)
                GC2 = weighted_graph(new_cliques_dict)
                clique_tree = Prim_algo(GC2)
                cliques_dict = new_cliques_dict
            end
        end

        isdir("..\\data\\decompositions\\blocks_"*repo_name) || mkpath("..\\data\\decompositions\\blocks_"*repo_name)
        f = open("..\\data\\decompositions\\blocks_"*repo_name*"\\$(instance[1:end-2])_sdp_blocks.txt", "w")
        if algo == "MD_real" || algo == "cholesky_real"
            nb_bus = Int(nb_nodes[instance]/2)
            for (clique, nodes_list) in cliques_dict
                for node in nodes_list
                    if node <= nb_bus
                        @printf(f, "%10s %20s\n", "B$clique", "VOLT_$(node)_Re")
                    else
                        bus = node - nb_bus
                        @printf(f, "%10s %20s\n", "B$clique", "VOLT_$(bus)_Im")
                    end
                end
            end
        else
            for (clique, nodes_list) in cliques_dict
                for node in nodes_list
                    @printf(f, "%10s %20s\n", "B$clique", "VOLT_$(node)_Re")
                    @printf(f, "%10s %20s\n", "B$clique", "VOLT_$(node)_Im")
                end
            end
        end
        close(f)

        isdir("..\\data\\decompositions\\cliquetree_"*repo_name) || mkpath("..\\data\\decompositions\\cliquetree_"*repo_name)
        f = open("..\\data\\decompositions\\cliquetree_"*repo_name*"\\$(instance[1:end-2])_sdp_cliquetree.txt", "w")

        if FUSION == "YES"
            for edge in LightGraphs.edges(clique_tree)
                clique1 = get_prop(clique_tree, src(edge), :name)[2:end]
                clique2 = get_prop(clique_tree, dst(edge), :name)[2:end]
                @printf(f, "%10s %20s\n", "B$clique1", "B$clique2")
            end
        else
            for edge in LightGraphs.edges(A)
                clique1 = src(edge)
                clique2 = dst(edge)
                @printf(f, "%10s %20s\n", "B$clique1", "B$clique2")
            end
        end
        close(f)
    end

    return
end

max_stars_repo_name: JulieSliwak/CliqueMerging

max_issues_repo_name: JulieSliwak/CliqueMerging

max_forks_repo_name: JulieSliwak/CliqueMerging


text:
ENV["JULIA_CXX_RTTI"]=1
using Cxx
using Test
#  include("../src/ROS.jl")
using ROS

ROS.@genNew
ROS.@updateAll

include("publisher.jl")
include("subscriber.jl")
include("serviceServer.jl")
include("serviceClient.jl")
include("actionlib_client.jl")
include("actionlib_server.jl")

as = Threads.@spawn testActionServer()
ac = Threads.@spawn testActionClient()
#  wait(as)
wait(ac)
s = Threads.@spawn testSubscriber()
p = Threads.@spawn testPublisher()
#  wait(s)
wait(p)
ss = Threads.@spawn testServiceServer()
sc = Threads.@spawn testServiceClient()
#  wait(ss)
wait(sc)

include("tf.jl")
include("time.jl")
include("roscpp.jl")
include("package.jl")
include("nodehandle.jl")


text:
export SymNumber
struct SymNumber <: Number
    ex::Expression
end

SymNumber(name::Symbol) = SymNumber(Variable(name))
SymNumber(x::Number) = SymNumber(Numeric(x))
SymNumber(x::Irrational{sym}) where sym = SymNumber(Constant(sym))
SymNumber(x::SymReal) = SymNumber(term(x))
SymNumber(x::SymComplex) = SymNumber(term(x))

term(x::SymNumber) = x.ex
Base.show(io::IO, t::SymNumber) = print(io, t.ex)

for (M, f, arity) in DiffRules.diffrules()
    op = :($M.$f)

    if arity == 1
        @eval $op(x::SymNumber) = track($op, x)
    elseif arity == 2
        @eval $op(x::SymNumber, y::SymNumber) = track($op, x, y)
    end
end

Base.promote_rule(::Type{<:SymNumber}, ::Type{<:Number}) = SymNumber
Base.promote_rule(::Type{<:SymNumber}, ::Type{<:SymComplex}) = SymNumber

Base.convert(::Type{<:SymNumber}, x::SymNumber) = x
Base.convert(::Type{<:SymNumber}, x::Number) = SymNumber(x)
Base.convert(::Type{<:SymNumber}, x::Irrational) = SymNumber(x)

Base.convert(::Type{<:SymNumber}, x::SymReal) = SymNumber(x)
Base.convert(::Type{<:SymNumber}, x::SymComplex) = SymNumber(x)


text:
"""
    const PointType = UInt8

Stores certain information about a grid point via bit-flags.

Right now there are:

    `const update_bit      = 0x01`
    `const undepleted_bit  = 0x02`
    `const pn_junction_bit = 0x04`

How to get information out of a PointType variable `pt`:
1. `pt & update_bit == 0` -> do not update this point (for fixed points)
2. `pt & update_bit >  0` -> do update this point
3. `pt & undepleted_bit > 0` -> this point is undepleted
4. `pt & pn_junction_bit > 0` -> this point belongs to the solid state detector. So it is in the volume of the n-type or p-type material.
"""
const PointType       = UInt8

# Point types for electric potential calculation
const update_bit      = 0x01 # parse(UInt8, "00000001", base=2) # 1 -> do update; 0 -> do not update
const undepleted_bit  = 0x02 # parse(UInt8, "00000010", base=2) # 0 -> depleted point; 1 -> undepleteded point
const pn_junction_bit = 0x04 # parse(UInt8, "00001000", base=2) # 0 -> point is not inside a bubble; 1 -> point is inside a bubble
# const bubble_bit      = 0x08 # parse(UInt8, "00000100", base=2) # 0 -> point is not inside a bubble; 1 -> point is inside a bubble

const max_pointtype_value = update_bit + undepleted_bit + pn_junction_bit #+ bubble_bit

"""
    PointTypes{T, N, S} <: AbstractArray{T, N}

PointTypes stores the point type of each grid point.
"""
struct PointTypes{T, N, S} <: AbstractArray{T, N}
    data::Array{PointType, N}
    grid::Grid{T, N, S}
end

@inline size(pts::PointTypes{T, N, S}) where {T, N, S} = size(pts.data)
@inline length(pts::PointTypes{T, N, S}) where {T, N, S} = length(pts.data)
@inline getindex(pts::PointTypes{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(pts.data, I...)
@inline getindex(pts::PointTypes{T, N, S}, i::Int) where {T, N, S} = getindex(pts.data, i)
@inline getindex(pts::PointTypes{T, N, S}, s::Symbol) where {T, N, S} = getindex(pts.grid, s)

function in(pt::AbstractCoordinatePoint{T}, pts::PointTypes{T, 3, S})::Bool where {T <: SSDFloat, S}
    i1::Int, i2::Int, i3::Int = searchsortednearest(pt, pts.grid)
    return pts.data[i1, i2, i3] & pn_junction_bit > 0
end

"""
    get_active_volume(pts::PointTypes{T}) where {T}
Returns an approximation of the active volume of the detector by summing up the cell volumes of
all depleted cells.
"""
function get_active_volume(pts::PointTypes{T, 3, :cylindrical}) where {T}
    active_volume::T = 0

    only_2d::Bool = length(pts.grid.φ) == 1
    cyclic::T = pts.grid.φ.interval.right

    r_ext::Vector{T} = get_extended_ticks(pts.grid.r)
    φ_ext::Vector{T} = get_extended_ticks(pts.grid.φ)
    z_ext::Vector{T} = get_extended_ticks(pts.grid.z)
    Δr_ext::Vector{T} = diff(r_ext)
    Δφ_ext::Vector{T} = diff(φ_ext)
    Δz_ext::Vector{T} = diff(z_ext)

    mpr::Vector{T} = midpoints(r_ext)
    mpφ::Vector{T} = midpoints(φ_ext)
    mpz::Vector{T} = midpoints(z_ext)
    Δmpφ::Vector{T} = diff(mpφ)
    Δmpz::Vector{T} = diff(mpz)
    Δmpr_squared::Vector{T} = T(0.5) .* ((mpr[2:end].^2) .- (mpr[1:end-1].^2))
    if r_ext[2] == 0
        Δmpr_squared[1] = T(0.5) * (mpr[2]^2)
    end

    isclosed::Bool = typeof(pts.grid.φ.interval).parameters[2] == :closed
    for iz in eachindex(pts.grid.z)
        if !isclosed || only_2d
            for iφ in eachindex(pts.grid.φ)
                for ir in eachindex(pts.grid.r)
                    pt::PointType = pts[ir, iφ, iz]
                    if (pt & pn_junction_bit > 0) && (pt & undepleted_bit == 0) && (pt & update_bit > 0)
                        dV::T = Δmpz[iz] * Δmpφ[iφ] * Δmpr_squared[ir]
                        active_volume += dV
                    end
                end
            end
        elseif isclosed && !only_2d
            for iφ in eachindex(pts.grid.φ)
                for ir in eachindex(pts.grid.r)
                    pt::PointType = pts[ir, iφ, iz]
                    if (pt & pn_junction_bit > 0) && (pt & undepleted_bit == 0) && (pt & update_bit > 0)
                        dV = Δmpz[iz] * Δmpφ[iφ] * Δmpr_squared[ir]
                        active_volume += if iφ == length(pts.φ) || iφ == 1
                            dV / 2
                        else
                            dV
                        end
                    end
                end
            end
        end
    end
    if cyclic > 0
        active_volume *= 2π / cyclic
    end

    f::T = 10^6
    return active_volume * f * Unitful.cm * Unitful.cm * Unitful.cm
end

function PointTypes(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = typeof(grid[1].ticks[1])
    S = get_coordinate_system(grid)
    N = get_number_of_dimensions(grid)
    PointTypes{T, N, S}( nt.values, grid )
end

Base.convert(T::Type{PointTypes}, x::NamedTuple) = T(x)

function NamedTuple(pts::PointTypes{T, 3}) where {T}
    return (
        grid = NamedTuple(pts.grid),
        values = pts.data,
    )
end

Base.convert(T::Type{NamedTuple}, x::PointTypes) = T(x)

@recipe function f( pts::PointTypes{T, 3, :cylindrical};
                    r = missing,
                    φ = missing,
                    z = missing ) where {T}
    g::Grid{T, 3, :cylindrical} = pts.grid

    seriescolor --> :viridis
    st --> :heatmap
    aspect_ratio --> 1

    cross_section::Symbol, idx::Int = if ismissing(φ) && ismissing(r) && ismissing(z)
        :φ, 1
    elseif !ismissing(φ) && ismissing(r) && ismissing(z)
        φ_rad::T = T(deg2rad(φ))
        while !(g.φ.interval.left <= φ_rad <= g.φ.interval.right)
            if φ_rad > g.φ.interval.right
                φ_rad -= g.φ.interval.right - g.φ.interval.left
            elseif φ_rad < g.φ.interval.left
                φ_rad += g.φ.interval.right - g.φ.interval.left
            end
        end
        :φ, searchsortednearest(g.φ, φ_rad)
    elseif ismissing(φ) && !ismissing(r) && ismissing(z)
        :r, searchsortednearest(g.r, T(r))
    elseif ismissing(φ) && ismissing(r) && !ismissing(z)
        :z, searchsortednearest(g.z, T(z))
    else
        error(ArgumentError, ": Only one of the keywords `r, φ, z` is allowed.")
    end
    value::T = if cross_section == :φ
        g.φ[idx]
    elseif cross_section == :r
        g.r[idx]
    elseif cross_section == :z
        g.z[idx]
    end

    @series begin
        clims --> (0, max_pointtype_value)
        if cross_section == :φ
            title --> "Point Type Map @$(cross_section) = $(round(rad2deg(value), sigdigits = 2))"
            xlabel --> "r / m"
            ylabel --> "z / m"
            size --> ( 400, 350 / (g.r[end] - g.r[1]) * (g.z[end] - g.z[1]) )
            g.r, g.z, pts.data[:, idx,:]'
        elseif cross_section == :r
            title --> "Point Type Map @$(cross_section) = $(round(value, sigdigits = 2))"
            g.φ, g.z, pts.data[idx,:,:]'
        elseif cross_section == :z
            title --> "Point Type Map @$(cross_section) = $(round(value, sigdigits = 2))"
            proj --> :polar
            g.φ, g.r, pts.data[:,:,idx]
        end
    end
end

@recipe function f( pts::PointTypes{T, 3, :cartesian};
                    x = missing,
                    y = missing,
                    z = missing ) where {T}
    g::Grid{T, 3, :cartesian} = pts.grid

    seriescolor --> :viridis
    st --> :heatmap
    aspect_ratio --> 1

    cross_section::Symbol, idx::Int = if ismissing(x) && ismissing(y) && ismissing(z)
        :x, 1
    elseif !ismissing(x) && ismissing(y) && ismissing(z)
        :x, searchsortednearest(g[:x], T(x))
    elseif ismissing(x) && !ismissing(y) && ismissing(z)
        :y, searchsortednearest(g[:y], T(y))
    elseif ismissing(x) && ismissing(y) && !ismissing(z)
        :z, searchsortednearest(g.z, T(z))
    else
        error(ArgumentError, ": Only one of the keywords `r, y, z` is allowed.")
    end
    value::T = if cross_section == :x
        g[:x][idx]
    elseif cross_section == :y
        g[:y][idx]
    elseif cross_section == :z
        g.z[idx]
    end

    @series begin
        clims --> (0, max_pointtype_value)
        title --> "Point Type Map @$(cross_section) = $(round(value, sigdigits = 2))"
        if cross_section == :x
            xlabel --> "y / m"
            ylabel --> "z / m"
            g[:y], g.z, pts.data[idx, :, :]'
        elseif cross_section == :y
            xlabel --> "x / m"
            ylabel --> "z / m"
            g[:x], g.z, pts.data[:, idx, :]'
        elseif cross_section == :z
            xlabel --> "x / m"
            ylabel --> "y / m"
            g[:x], g[:y], pts.data[:,:,idx]'
        end
    end
end


text:
## Optimization =======================================

function perturb(rng::AbstractRNG, p, Cp)
    d = MvNormal(pmean(p), 1.1Cp + 1e-12I)
    Particles(rng, nparticles(p[1]), d)
end

"""
    res = optimize([rng::AbstractRNG,] f,p,τ=1,iters=10000)

Find the minimum of Function `f`, starting with initial distribution described by `p::Vector{Particles}`. `τ` is the initial temperature.
"""
function optimize(rng::AbstractRNG,f,p,τ=1; τi=1.005, iters=10000, tol=1e-8)
    p = deepcopy(p)
    N = nparticles(p[1])
    we = zeros(N)
    for i = 1:iters
        y  = -(f(p).particles)
        we .= exp.(τ.*y)
        j  = sample(rng,1:N, ProbabilityWeights(we), N, ordered=true)
        foreach(x->(x.particles .= x.particles[j]), p); # @test length(unique(p[1].particles)) == length(unique(j))
        Cp = pcov(p)
        tr(Cp) < tol && (@info "Converged at iteration $i"; return p)
        p = perturb(rng, p, Cp)
        τ *= τi
    end
    p
end
optimize(f,p,τ=1; kwargs...) = optimize(Random.GLOBAL_RNG,f,p,τ; kwargs...)


text:
export FfOrdFracIdl

################################################################################
#
#  Fractional Ideals
#
################################################################################

mutable struct FfOrdFracIdl
  order::GenericRound2.Order
  num::FfOrdIdl
  den::RingElem
  norm::RingElem
  basis_matrix::FakeFracFldMat
  basis_mat_inv::FakeFracFldMat

  function FfOrdFracIdl(O::GenericRound2.Order)
    z = new()
    z.order = O
    return z
  end

  function FfOrdFracIdl(a::FfOrdIdl, b::RingElem)
    z = new()
    z.order = order(a)
    b = Hecke.AbstractAlgebra.MPolyFactor.make_monic(b)
    z.num = a
    z.den = b
    return z
  end

   function FfOrdFracIdl(a::FfOrdIdl)
    z = new()
    O = order(a)
    z.order = order(a)
    z.num = a
    z.den = O.R(1)
    return z
  end

  function FfOrdFracIdl(O::GenericRound2.Order, a::FakeFracFldMat)
    z = new()
    z.order = O
    z.basis_matrix = a
    return z
  end
end

Hecke.order(a::FfOrdFracIdl) = a.order

function AbstractAlgebra.iszero(x::FfOrdFracIdl)
  return iszero(numerator(x))
end

################################################################################
#
#  Basis matrix
#
################################################################################

function assure_has_basis_matrix(a::FfOrdFracIdl)
  if isdefined(a, :basis_matrix)
    return nothing
  end
  if !isdefined(a, :num)
    error("Not a valid fractional ideal")
  end

  a.basis_matrix = FakeFracFldMat(basis_matrix(numerator(a, copy = false)), denominator(a))
  return nothing
end

function Hecke.basis_matrix(x::FfOrdFracIdl; copy::Bool = true)
  assure_has_basis_matrix(x)
  if copy
    return deepcopy(x.basis_matrix)
  else
    return x.basis_matrix
  end
end

################################################################################
#
#  Numerator and denominator
#
################################################################################

function assure_has_numerator_and_denominator(a::FfOrdFracIdl)
  if isdefined(a, :num) && isdefined(a, :den)
    return nothing
  end
  if !isdefined(a, :basis_matrix)
    error("Not a valid fractional ideal")
  end

  a.num = FfOrdIdl(order(a), numerator(basis_matrix(a, copy = false)))
  a.den = denominator(basis_matrix(a, copy = false))
  return nothing
end

function Base.numerator(x::FfOrdFracIdl; copy::Bool = true)
  assure_has_numerator_and_denominator(x)
  if copy
    return deepcopy(x.num)
  else
    return x.num
  end
end

function Base.denominator(x::FfOrdFracIdl; copy::Bool = true)
  assure_has_numerator_and_denominator(x)
  if copy
    return deepcopy(x.den)
  else
    return x.den
  end
end

################################################################################
#
#  Binary operations
#
################################################################################

function Base.prod(a::T, b::T) where T <: FfOrdFracIdl
  A = a.num*b.num
  return FfOrdFracIdl(A, a.den*b.den)
end

Base.:*(A::T, B::T) where T <: FfOrdFracIdl = prod(A, B)

################################################################################
#
#  Powering
#
################################################################################

function Base.:^(A::FfOrdFracIdl, a::Int)

  O = order(A)
  if a == 0
    B = FfOrdFracIdl(ideal(order(A), 1), O.R(1))
    return B
  end

  if a == 1
    return A
  end

  if a < 0
    return inv(A^(-a))
  end

  if a == 2
    return A*A
  end

  if mod(a, 2) == 0
    return (A^div(a, 2))^2
  else
    return A * A^(a - 1)
  end
end

################################################################################
#
#  Simplification
#
################################################################################

function Hecke.simplify(A::FfOrdFracIdl)
  assure_has_numerator_and_denominator(A)
  if isone(A.den)
    return A
  end

  b = basis_matrix(A.num)
  g = gcd(denominator(A), content(b))

  if g != 1
    A.num = divexact(A.num, g)
    A.den = divexact(A.den, g)
  end

  return A
end

################################################################################
#
#   Is integral
#
################################################################################

Hecke.is_integral(I::FfOrdIdl) = true

function Hecke.is_integral(I::FfOrdFracIdl)
  simplify(I)
  return denominator(I) == 1
end

################################################################################
#
#  Ad hoc binary operations
#
################################################################################

function Base.:*(A::FfOrdIdl, B::FfOrdFracIdl)
  z = FfOrdFracIdl(A*numerator(B, copy = false), denominator(B))
  return z
end

Base.:*(A::FfOrdFracIdl, B::FfOrdIdl) = FfOrdFracIdl(numerator(A, copy = false)*B, denominator(A))

function Base.:*(x::GenericRound2.OrderElem, y::FfOrdFracIdl)
  #parent(x) !== order(y) && error("Orders of element and ideal must be equal")
  return FfOrdIdl(parent(x), x) * y
end

Base.:*(x::FfOrdFracIdl, y::GenericRound2.OrderElem) = y * x

################################################################################
#
#  Colon
#
################################################################################

function Hecke.colon(I::FfOrdFracIdl, J::FfOrdFracIdl)
  # Let I = a/c and J = b/d, a and b integral ideals, c, d \in Z, then
  # \{ x \in K | xJ \subseteq I \} = \{ x \in K | xcb \subseteq da \}

  O = order(I)
  II = numerator(I, copy = false)*O(denominator(J, copy = false))
  JJ = numerator(J, copy = false)*O(denominator(I, copy = false))
  return Hecke.colon(II, JJ)
end


text:
using Distributed
@everywhere using StatsBase, DataStructures, IterTools, SharedArrays

@everywhere @views correlate(ts, i; max_correlation_length=300) = takewhile(<=(max_correlation_length), imap(t -> t - ts[i], ts[i:end]));
@everywhere @views correlate(ts; max_correlation_length=300) = idx -> correlate(ts, idx; max_correlation_length=max_correlation_length);
# ts = SharedArray{UInt}(parse.(UInt, readlines("./channels/ch0.txt")));
@everywhere ⊎ = merge!

y(ts) = Iterators.accumulate(⊎ , imap(counter ∘ correlate(ts), 1:length(ts)))

G₂(times, n::Int; max_correlation_length=300)::Accumulator{UInt, UInt} = @distributed (⊎) for i ∈ 1:n
                @views (counter ∘ correlate(times, max_correlation_length=max_correlation_length))(i)
                end

g₂(times::SharedArray{UInt}, n::Int)::Dict{UInt, Float64} = begin
    G2 = G₂(times, n);
    norm(G2)
end

norm(x) = Dict{Any, Float64}(Pair.(keys(x), values(x) ./ mean(values(x))))

bin(dict, size_of_bin) = Dict(Pair(i, sum(get(dict, j, (mean∘ values)(dict)) for j in i:i+size_of_bin)) for i in 3:size_of_bin:(length(dict)))


text:
# using DataFrames
using ..ChemistryFeaturization.Utils.SpeciesFeatureUtils
using ..ChemistryFeaturization.Codec

include("abstractfeatures.jl")

"""
    SpeciesFeatureDescriptor{A}

A feature object that encodes features associated with individual atoms (collectively depicted by a representation of type `A`) that depend upon their local environment in some way (if your feature is defined only by elemental identity, you should use ElementFeatureDescriptor!)

Type parameter represents the structure representation(s) from which this feature descriptor is able to compute features.

## Fields
- `name::String`: the name of the feature
- `compute_f`: function that takes in a structure of type `A` and returns values of this feature for every atom in the structure
- `encoder_decoder::AbstractCodec`: codec that encodes/decodes values of this feature
- `categorical::Bool`: flag for whether the feature is categorical or continuous-valued
- `encodable_elements::Vector{String}`: list of elements (by symbol) that can be encoded by this feature
"""
struct SpeciesFeatureDescriptor{A,C<:AbstractCodec} <: AbstractAtomFeatureDescriptor
    name::String
    compute_f::Function
    encoder_decoder::C
    categorical::Bool
    encodable_elements::Vector{String}
end

"""
    SpeciesFeatureDescriptor(name::String)

Construct a `SpeciesFeatureDescriptor` from one of the built-in options (currently supports a variety of MolecularGraph functions, see Utils.SpeciesFeatureUtils.sfd_names_props for details).
"""
function SpeciesFeatureDescriptor(name::String)
    @assert name in keys(sfd_names_props) "That species feature isn't one of the built-in ones, you'll have to construct it directly. Consult ChemistryFeaturiztion.Utils.SpeciesFeatureUtils.sfd_names_props for built-in species features."

    info = sfd_names_props[name]
    categorical = info[:categorical]
    local codec
    if categorical
        possible_vals = info[:possible_vals]
        if length(possible_vals) == 2
            @assert possible_vals == [true, false]
            codec = DirectCodec(1)
        else
            codec = OneHotOneCold(true, possible_vals)
        end
    else
        # TODO: figure out default binning situation for continuous-valued SFD's
        #codec = OneHotOneCold(false, )
    end
    SpeciesFeatureDescriptor{info[:A],typeof(codec)}(
        name,
        info[:compute_f],
        codec,
        categorical,
        info[:encodable_elements],
    )
end

# pretty printing, long version
function Base.show(io::IO, ::MIME"text/plain", fd::SpeciesFeatureDescriptor{A}) where {A}
    st = "$(typeof(fd)) $(fd.name):\n   categorical: $(fd.categorical)\n   works on: $(A)"
    print(io, st)
end

encodable_elements(fd::SpeciesFeatureDescriptor) = fd.encodable_elements

function get_value(sfd::SpeciesFeatureDescriptor{A}, a::AbstractAtoms{<:A}) where {A}
    @assert all([el in encodable_elements(sfd) for el in elements(a)]) "Feature $(sfd.name) cannot encode some element(s) in this structure!"
    sfd.compute_f(a.structure)
end


text:
using Parameters

@with_kw struct GeneralParams

    # input parameters (Nusselt scaling)
    Re::Float64 = 40.8
    Ka::Float64 = 3923.0
    Ct::Float64 = 0.0
    k_x::Float64 = 0.0729
    k_z::Float64 = 0.0729
    h0::Float64 = 0.9147

    # derived (Shkadov scaling)
    eta::Float64 = (3.0 * Re)^(4 / 9) / Ka^(2 / 3)
    delta::Float64 = (3.0 * Re)^(11 / 9) / Ka^(1 / 3)
    zeta::Float64 = Ct * (3.0 * Re)^(2 / 9) / Ka^(1 / 3)
    k0x::Float64 = k_x / sqrt(eta)
    k0z::Float64 = k_z / sqrt(eta)

    Lx::Float64 = 2 * pi / k0x
    Lz::Float64 = 2 * pi / k0z

    # initial conditions
    stSine::Float64 = 0.2
    spSine::Float64 = 0.2
    stCosine::Int64 = 0
    spCosine::Int64 = 0
    stNwaves::Int64 = 1
    spNwaves::Int64 = 1
    stNoise::Int64 = 0
    spNoise::Int64 = 0

    # mesh parameters
    # 2D simulation
    # NN1::Int64 = 5
    # NN2::Int64 = 1
    # N1::Int64 = 32
    # N2::Int64 = 1
    # M1::Int64 = 64
    # M2::Int64 = 2

    # 3D simulation
    NN1::Int64 = 5
    NN2::Int64 = 5
    N1::Int64 = 32
    N2::Int64 = 32
    M1::Int64 = 64
    M2::Int64 = 64
end

parameters = GeneralParams()


text:
hexIDs = UInt16[0x0001, 0x0002, 0x0003]
@test findhex(0x0000, hexIDs) == -1
@test findhex(0x0001, hexIDs) == 1
@test findhex(0x0002, hexIDs) == 2
@test findhex(0x0003, hexIDs) == 3

# win32 with gaps
if has_restricted
  cfile = path*"/SampleFiles/Restricted/03_02_27_20140927.euc.ch"
  printstyled("  win32\n", color=:light_green)
  @info(string(timestamp(), ": win32 tests use an SSL-encrypted tarball."))

  redirect_stdout(out) do
    fname = path*"/SampleFiles/Restricted/2014092709*.cnt"
    cfile = path*"/SampleFiles/Restricted/03_02_27_20140927.euc.ch"
    S = read_data("win32", fname, cf=cfile, v=3)

    # There should be 8 channels
    @test S.n==8

    # There should be exactly 360000 points per channel (all are 100 Hz)
    nx = [length(S.x[i]) for i=1:1:S.n]
    @test minimum(nx)==360000
    @test maximum(nx)==360000

    # Check against SAC files
    printstyled("    checking read integrity against HiNet SAC files\n", color=:light_green)
    testfiles = path*"/SampleFiles/Restricted/" .* ["20140927000000.V.ONTA.E.SAC",
                                                    "20140927000000.V.ONTA.H.SAC",
                                                    "20140927000000.V.ONTA.N.SAC",
                                                    "20140927000000.V.ONTA.U.SAC",
                                                    "20140927000000.V.ONTN.E.SAC",
                                                    "20140927000000.V.ONTN.H.SAC",
                                                    "20140927000000.V.ONTN.N.SAC",
                                                    "20140927000000.V.ONTN.U.SAC"]

    # SAC files prepared in SAC with these commands from day-long Ontake files
    # beginning at midnight Japan time converted to SAC with win32 precompiled
    # utilities
    #
    # SAC commands that generate these from day-long SAC files:
    # cut b 32400 35999.99
    # r /home/josh/SAC/20140927000000.V.*.SAC"
    # ch b 0 nzjday 270 nzhour 09
    # w over
    # q

    U = SeisData()
    for f in testfiles
      T = read_data("sac", f)[1]
      push!(U, T)
    end

    # ID correspondence
    j = Array{Int64,1}(undef, S.n)
    for (n,i) in enumerate(S.id)
      id = split(i, '.')
      id[3] = ""
      id[2] = "V"*id[2]
      id[1] = "JP"
      c = id[4][3:3]
      if c == "Z"
        id[4] = "U"
      else
        id[4] = c
      end
      id = join(id, '.')
      j[n] = findid(id, U)
    end

    inds = Array{Array{Int64,1}}(undef, S.n)
    for (n,i) in enumerate(j)
      inds[n] = findall(abs.(S.x[n]-U.x[i]).>eps())
    end

    # The only time gaps should be what are in the logs:
    # ┌ Warning: Time gap detected! (15.0 s at V.ONTA.H, beginning 2014-09-27T09:58:00)
    # └ @ SeisIO ~/.julia/dev/SeisIO/src/Formats/Win32.jl:137
    # ┌ Warning: Time gap detected! (15.0 s at V.ONTA.U, beginning 2014-09-27T09:58:00)
    # └ @ SeisIO ~/.julia/dev/SeisIO/src/Formats/Win32.jl:137
    # ┌ Warning: Time gap detected! (15.0 s at V.ONTA.N, beginning 2014-09-27T09:58:00)
    # └ @ SeisIO ~/.julia/dev/SeisIO/src/Formats/Win32.jl:137
    # ┌ Warning: Time gap detected! (15.0 s at V.ONTA.E, beginning 2014-09-27T09:58:00)
    # └ @ SeisIO ~/.julia/dev/SeisIO/src/Formats/Win32.jl:137

    for k = 1:S.n
      i = j[k]
      if !isempty(inds[i])
        @test length(inds[i]) == 1500
        @test div(first(inds[i]), 6000) == 58
        @test div(last(inds[i]), 6000) == 58
        r₀ = rem(first(inds[i]), 6000)
        r₁ = rem(last(inds[i]), 6000)
        @test round((r₁ - r₀)/100, digits=1) == 15.0
      end
    end
  end

  # Converting the mean to single-point precision gives exactly the same
  # result as SAC conversion from win32; however, the average computed
  # over an hour will be slightly different.
  #
  # This can be verified using the script "ontake_test.jl" in ../../internal_tests/

  # Now test the other two bits types, 4-bit Int ...
  printstyled("    testing Int4 and Int24 handling\n", color=:light_green)
  fname = path*"/SampleFiles/Restricted/2014092700000302.cnt"
  cfile = path*"/SampleFiles/Restricted/03_02_27_20140927*ch"
  S = SeisData()
  redirect_stdout(out) do
    read_data!(S, "win32", fname, cf=cfile, v=3)
  end
  i = findid("V.ONTA.23.EHH", S)
  @test length(S.x[i]) == 60*S.fs[i]
  @test maximum(S.x[i]) == 11075.0
  @test minimum(S.x[i]) == -5026.0

  # ...and 24-bit bigendian Int...
  fname = path*"/SampleFiles/Restricted/2014092712000302.cnt"
  read_data!(S, "win32", fname, cf=cfile)
  @test length(S.x[1]) == round(Int64, 60*S.fs[1]) == S.t[1][end,1]
  ii = findlast(S.id.=="V.ONTA.23.EHH")
  @test maximum(S.x[ii]) == 14896.0
  @test minimum(S.x[ii]) == -12651.0
  for id in unique(S.id)
    @test length(findall(S.id.==id))==2
  end

  # ...and 32-bit bigendian Int ...
  cfile = path*"/SampleFiles/Restricted/chandata_20140927.txt"
  fname = path*"/SampleFiles/Restricted/2014092712370207VM.cnt"
  read_data!(S, "win32", fname, cf=cfile)
else
  printstyled("  win32 data format not tested (files not found)\n", color=:red)
end


text:
module TestBioServices

using BioServices.EUtils
using BioServices.UMLS
using BioServices.GGGenome
using BioServices.EBIProteins
using XMLDict
using Test

all_tests = [
    ("eutils.jl",      "      Testing: EUtils"),
    ("umls.jl",        "      Testing: UMLS"),
    ("gggenome.jl",    "      Testing: GGGenome"),
    ("ebiproteins.jl", "      Testing: EBIProteins")
    ]

println("Running tests:")

for (t, test_string) in all_tests
    println("-----------------------------------------")
    println("-----------------------------------------")
    println(test_string)
    println("-----------------------------------------")
    println("-----------------------------------------")
    include(t)
end

end # TestBioServices

max_stars_repo_name: BioJulia/BioServices

max_issues_repo_name: BioJulia/BioServices

max_forks_repo_name: BioJulia/BioServices


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAuthorizationV1beta1NonResourceRule <: SwaggerModel
    nonResourceURLs::Any # spec type: Union{ Nothing, Vector{String} } # spec name: nonResourceURLs
    verbs::Any # spec type: Union{ Nothing, Vector{String} } # spec name: verbs

    function IoK8sApiAuthorizationV1beta1NonResourceRule(;nonResourceURLs=nothing, verbs=nothing)
        o = new()
        validate_property(IoK8sApiAuthorizationV1beta1NonResourceRule, Symbol("nonResourceURLs"), nonResourceURLs)
        setfield!(o, Symbol("nonResourceURLs"), nonResourceURLs)
        validate_property(IoK8sApiAuthorizationV1beta1NonResourceRule, Symbol("verbs"), verbs)
        setfield!(o, Symbol("verbs"), verbs)
        o
    end
end # type IoK8sApiAuthorizationV1beta1NonResourceRule

const _property_map_IoK8sApiAuthorizationV1beta1NonResourceRule = Dict{Symbol,Symbol}(Symbol("nonResourceURLs")=>Symbol("nonResourceURLs"), Symbol("verbs")=>Symbol("verbs"))
const _property_types_IoK8sApiAuthorizationV1beta1NonResourceRule = Dict{Symbol,String}(Symbol("nonResourceURLs")=>"Vector{String}", Symbol("verbs")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiAuthorizationV1beta1NonResourceRule }) = collect(keys(_property_map_IoK8sApiAuthorizationV1beta1NonResourceRule))
Swagger.property_type(::Type{ IoK8sApiAuthorizationV1beta1NonResourceRule }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAuthorizationV1beta1NonResourceRule[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthorizationV1beta1NonResourceRule }, property_name::Symbol) =  _property_map_IoK8sApiAuthorizationV1beta1NonResourceRule[property_name]

function check_required(o::IoK8sApiAuthorizationV1beta1NonResourceRule)
    (getproperty(o, Symbol("verbs")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAuthorizationV1beta1NonResourceRule }, name::Symbol, val)
end


text:
"
Do I really need this? This is my hack to allow both DataFrames and Arrays to be passed to
Poverty and Weights routines
"
# ArrayOrFrame = Union{ DataFrame, AbstractArray{<:Number,2}}


text:
using MedicalTransforms
using Test
using TestSetExtensions

using MedicalTransforms: crop_foreground


text:
module TileCoding

mutable struct IHT
    size::Integer
    overfullCount::Int64
    dictionary::Dict{Array{Int64, 1}, Int64}
    IHT(sizeval) = new(sizeval, 0, Dict{Array{Int64, 1}, Int64}())
end

capacity(iht::IHT) = iht.size
count(iht::IHT) = length(iht.dictionary)
fullp(iht::IHT) = length(iht.dictionary) >= count(iht)

function getindex!(iht::IHT, obj::Array{Int64, 1}, readonly=false)
    d = iht.dictionary
    if obj in keys(d)
        return d[obj]::Int64
    elseif readonly
        return nothing
    end
    iht_size = capacity(iht)
    iht_count = count(iht)
    if count(iht) >= capacity(iht)
        if iht.overfullCount==0
            println("IHT full, starting to allow collisions")
        end
        iht.overfullCount += 1
        return hash(obj) % capacity(iht)
    end
    d[obj] = iht_count
    return iht_count
end

hashcoords!(coordinates, m, readonly=false) = nothing
hashcoords!(coordinates, m::IHT, readonly=false) = getindex!(m, coordinates, readonly)
hashcoords!(coordinates, m::Integer, readonly=false) = hash(tuple(coordinates)) % m

function tiles!(ihtORsize, numtilings, floats, ints=[], readonly=false)
    qfloats = [floor(f*numtilings) for f in floats]
    tiles = zeros(Int64, numtilings)
    for tiling = 1:numtilings
        tilingX2 = tiling*2
        coords = [tiling]::Array{Int64, 1}
        b = tiling
        for (q_idx, q) in enumerate(qfloats)
            append!(coords, floor((q + b) / numtilings))
            b += tilingX2
        end
        append!(coords, ints)
        tiles[tiling] = hashcoords!(coords, ihtORsize, readonly)
    end
    return tiles
end

function tileswrap!(ihtORsize, numtilings, floats, wrapwidths, ints=[], readonly=false)
    qfloats = [floor(f*numtilings) for f in floats]
    tiles = zeros(Int64, numtilings)
    for tiling = 1:numtilings
        tilingX2 = tiling*2
        coords = convert(Array{Any}, [tiling])
        b = tiling
        for (q_idx, q) in enumerate(qfloats)
            width = nothing
            if length(wrapwidths) >= q_idx
                width = wrapwidths[q_idx]
            end
            c = floor((q + b%numtilings) / numtilings)
            append!(coords, width==nothing ? c : c%width)
            b += tilingX2
        end
        append!(coords, ints)
        tiles[tiling] = hashcoords!(coords, ihtORsize, readonly)
    end
    return tiles
end

end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
"""
    Rician(ν, σ)

The *Rician distribution* with parameters `ν` and `σ` has probability density function:

```math
f(x; \\nu, \\sigma) = \\frac{x}{\\sigma^2} \\exp\\left( \\frac{-(x^2 + \\nu^2)}{2\\sigma^2} \\right) I_0\\left( \\frac{x\\nu}{\\sigma^2} \\right).
```

If shape and scale parameters `K` and `Ω` are given instead, `ν` and `σ` may be computed from them:

```math
\\sigma = \\sqrt{\\frac{\\Omega}{2(K + 1)}}, \\quad \\nu = \\sigma\\sqrt{2K}
```

```julia
Rician()         # Rician distribution with parameters ν=0 and σ=1
Rician(ν, σ)     # Rician distribution with parameters ν and σ

params(d)        # Get the parameters, i.e. (ν, σ)
shape(d)         # Get the shape parameter K = ν²/2σ²
scale(d)         # Get the scale parameter Ω = ν² + 2σ²
```

External links:

* [Rician distribution on Wikipedia](https://en.wikipedia.org/wiki/Rice_distribution)

"""
struct Rician{T<:Real} <: ContinuousUnivariateDistribution
    ν::T
    σ::T
    Rician{T}(ν, σ) where {T} = new{T}(ν, σ)
end

function Rician(ν::T, σ::T; check_args::Bool=true) where {T<:Real}
    @check_args Rician (ν, ν ≥ zero(ν)) (σ, σ ≥ zero(σ))
    return Rician{T}(ν, σ)
end

Rician() = Rician{Float64}(0.0, 1.0)
Rician(ν::Real, σ::Real; check_args::Bool=true) = Rician(promote(ν, σ)...; check_args=check_args)
Rician(ν::Integer, σ::Integer; check_args::Bool=true) = Rician(float(ν), float(σ); check_args=check_args)

@distr_support Rician 0.0 Inf

#### Conversions

function convert(::Type{Rician{T}}, ν::Real, σ::Real) where T<:Real
    Rician(T(ν), T(σ))
end

function convert(::Type{Rician{T}}, d::Rician{S}) where {T <: Real, S <: Real}
    Rician(T(d.ν), T(d.σ); check_args=false)
end

#### Parameters

shape(d::Rician) = d.ν^2 / (2 * d.σ^2)
scale(d::Rician) = d.ν^2 + 2 * d.σ^2

params(d::Rician) = (d.ν, d.σ)
partype(d::Rician{T}) where {T<:Real} = T

#### Statistics

# helper
_Lhalf(x) = exp(x/2) * ((1-x) * besseli(zero(x), -x/2) - x * besseli(oneunit(x), -x/2))

mean(d::Rician) = d.σ * sqrthalfπ * _Lhalf(-d.ν^2/(2 * d.σ^2))
var(d::Rician) = 2 * d.σ^2 + d.ν^2 - halfπ * d.σ^2 * _Lhalf(-d.ν^2/(2 * d.σ^2))^2

function mode(d::Rician)
    m = mean(d)
    _minimize_gss(x -> -pdf(d, x), zero(m), m)
end

# helper: 1D minimization using Golden-section search
function _minimize_gss(f, a, b; tol=1e-12)
    ϕ = (√5 + 1) / 2
    c = b - (b - a) / ϕ
    d = a + (b - a) / ϕ
    while abs(b - a) > tol
        if f(c) < f(d)
            b = d
        else
            a = c
        end
        c = b - (b - a) / ϕ
        d = a + (b - a) / ϕ
    end
    (b + a) / 2
end

#### PDF/CDF/quantile delegated to NoncentralChisq

function quantile(d::Rician, x::Real)
    ν, σ = params(d)
    return sqrt(quantile(NoncentralChisq(2, (ν / σ)^2), x)) * σ
end

function cquantile(d::Rician, x::Real)
    ν, σ = params(d)
    return sqrt(cquantile(NoncentralChisq(2, (ν / σ)^2), x)) * σ
end

function pdf(d::Rician, x::Real)
    ν, σ = params(d)
    result = 2 * x / σ^2 * pdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2)
    return x < 0 || isinf(x) ? zero(result) : result
end

function logpdf(d::Rician, x::Real)
    ν, σ = params(d)
    result = log(2 * abs(x) / σ^2) + logpdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2)
    return x < 0 || isinf(x) ? oftype(result, -Inf) : result
end

function cdf(d::Rician, x::Real)
    ν, σ = params(d)
    result = cdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2)
    return x < 0 ? zero(result) : result
end

function logcdf(d::Rician, x::Real)
    ν, σ = params(d)
    result = logcdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2)
    return x < 0 ? oftype(result, -Inf) : result
end

#### Sampling

function rand(rng::AbstractRNG, d::Rician)
    x = randn(rng) * d.σ + d.ν
    y = randn(rng) * d.σ
    hypot(x, y)
end

#### Fitting

# implementation based on the Koay inversion technique
function fit(::Type{<:Rician}, x::AbstractArray{T}; tol=1e-12, maxiters=500) where T
    μ₁ = mean(x)
    μ₂ = var(x)
    r = μ₁ / √μ₂
    if r < sqrt(π/(4-π))
        ν = zero(float(T))
        σ = scale(fit(Rayleigh, x))
    else
        ξ(θ) = 2 + θ^2 - π/8 * exp(-θ^2 / 2) * ((2 + θ^2) * besseli(0, θ^2 / 4) + θ^2 * besseli(1, θ^2 / 4))^2
        g(θ) = sqrt(ξ(θ) * (1+r^2) - 2)
        θ = g(1)
        for j in 1:maxiters
            θ⁻ = θ
            θ = g(θ)
            abs(θ - θ⁻) < tol && break
        end
        ξθ = ξ(θ)
        σ = convert(float(T), sqrt(μ₂ / ξθ))
        ν = convert(float(T), sqrt(μ₁^2 + (ξθ - 2) * σ^2))
    end
    Rician(ν, σ)
end

# Not implemented:
#   skewness(d::Rician)
#   kurtosis(d::Rician)
#   entropy(d::Rician)
#   mgf(d::Rician, t::Real)
#   cf(d::Rician, t::Real)


text:
import Base: gcd, lcm

function gcd(a::SimplePolynomial, b::SimplePolynomial)
    #try
    d = integerize(_gcd(a, b))
    return d
    # catch
    #     aa = big(a)
    #     bb = big(b)
    #     dd = _gcd(aa,bb)
    #     d = integerize(dd)
    #     return d
    # end
end

function _gcd(a::SimplePolynomial, b::SimplePolynomial)

    if a == 0 && b == 0
        error("The two arguments cannot both be zero")
    end

    if a == 0
        return b
    end

    if b == 0
        return a
    end

    q, r = divrem(a, b)
    return monic(_gcd(b, monic(r)))
end

function lcm(a::SimplePolynomial, b::SimplePolynomial)
    q, r = divrem(a * b, gcd(a, b))
    return monic(q)
end


text:
using TutorialsForDevelopers
using Documenter
using Literate

repo_src = joinpath(@__DIR__,"..","src")
pages_dir = joinpath(@__DIR__,"src","pages")
Sys.rm(pages_dir;recursive=true,force=true)

pages = ["Introduction"=>"index.md"]

for (i,(file_jl,title)) in enumerate(TutorialsForDevelopers.files)

  file, = splitext(file_jl)
  file_md = file * ".md"
  tutorial_title = "# # Tutorial $i: $title"

  Literate.markdown(
    joinpath(repo_src,file_jl),
    pages_dir; name=file, preprocess = content -> string(tutorial_title,"\n",content))

  push!(pages,"$i $title"=>joinpath("pages",file_md))

end

makedocs(;
    modules=[TutorialsForDevelopers],
    authors="Francesc Verdugo <fverdugo@cimne.upc.edu> and contributors",
    repo="https://github.com/gridap/TutorialsForDevelopers/blob/{commit}{path}#L{line}",
    sitename="TutorialsForDevelopers",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://gridap.github.io/TutorialsForDevelopers",
        assets=String[],
    ),
    pages=pages,
)

deploydocs(;
    repo="github.com/gridap/TutorialsForDevelopers",
)

max_stars_repo_name: gridap/TutorialsForDevelopers

max_issues_repo_name: gridap/TutorialsForDevelopers

max_forks_repo_name: gridap/TutorialsForDevelopers


text:
using Oceananigans
using Oceananigans.Units
using GLMakie
using Printf

grid = RectilinearGrid(size=(64, 64), x=(0, 128), z=(-64, 0), halo=(3, 3),
                       topology=(Periodic, Flat, Bounded))

boundary_conditions = (; b = FieldBoundaryConditions(top = FluxBoundaryCondition(1e-8)))
closure = ScalarDiffusivity(ν=1e-4, κ=1e-4)
coriolis = FPlane(f=1e-4)

model = NonhydrostaticModel(; grid, closure, coriolis, boundary_conditions,
                            advection = UpwindBiasedFifthOrder(),
                            timestepper = :RungeKutta3,
                            tracers = :b,
                            buoyancy = BuoyancyTracer())

bᵢ(x, y, z) = 1e-6 * z
ϵ(x, y, z) = 1e-6 * randn() # noise
set!(model, b=bᵢ, u=ϵ, v=ϵ, w=ϵ)

simulation = Simulation(model, Δt=2minutes, stop_time=12hours)

wizard = TimeStepWizard(cfl=0.8, max_change=1.1, max_Δt=2minutes)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

progress(sim) = @printf("Iter: %d, time: %s, Δt: %s, max|w|: %s\n",
                        iteration(sim), prettytime(sim), prettytime(sim.Δt), maximum(abs, model.velocities.w))

simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

# Turbulent kinetic energy
u, v, w = model.velocities
e = Field((u^2 + v^2 + w^2) / 2)
ξ = Field(∂x(w) + ∂z(u))

simulation.output_writers[:fields] = JLD2OutputWriter(model, (; e, w, ξ, b = model.tracers.b),
                                                      schedule = TimeInterval(1minute),
                                                      prefix = "ocean_convection",
                                                      field_slicer = nothing,
                                                      force = true)

run!(simulation)

# Visualization

data_path = simulation.output_writers[:fields].filepath
e = FieldTimeSeries(data_path, "e")
b = FieldTimeSeries(data_path, "b")
x, y, z = nodes(e)
Nt = length(e.times)

fig = Figure(resolution=(1600, 800))
axe = Axis(fig[1, 1])
axξ = Axis(fig[1, 2])
slider = Slider(fig[3, :], range=1:Nt, startvalue=1)
n = slider.value

eⁿ = @lift dropdims(interior(e[$n]), dims=2)
bⁿ = @lift dropdims(interior(b[$n]), dims=2)
hm = heatmap!(axe, x, z, eⁿ; colormap=:solar, colorrange=(0, maximum(e) / 4))
Colorbar(fig[2, 1], hm, vertical=false, flipaxis=false, label="Turbulent kinetic energy (m² s⁻²)")

hm = heatmap!(axξ, x, z, bⁿ; colormap=:deep, colorrange=(minimum(b), maximum(b)))
Colorbar(fig[2, 2], hm, vertical=false, flipaxis=false, label="Buoyancy (m s⁻²)")

title = @lift "Idealized ocean convection at t = " * prettytime(e.times[$n])
Label(fig[0, :], title)

display(fig)

#record(fig, "idealized_ocean_convection.mp4", 1:Nt, framerate=12) do nn
#    @info "Rendering frame $nn of $Nt...
#    n[] = nn
#end


text:
mutable struct joint_pd_override_t <: LCMType
    timestamp::Int64
    position_ind::Int32
    qi_des::Float64
    qdi_des::Float64
    kp::Float64
    kd::Float64
    weight::Float64
end

@lcmtypesetup(joint_pd_override_t)