text:
using Oscar
using Polymake
using InvertedIndices
using Combinatorics
using LinearAlgebra
# StackyFan struct
export StackyFan, makeStackyFan, addStackStructure
export stackyWeights, getRayStack, encode
# Fan
export toric_blowup, makeSimplicial, makeSmooth, starSubdivision
export getConeRank, getDimension, getConeFaces
export coneMultiplicity, getMultiplicity
export findFaceContainingRay, findMinimalCone, interiorPoints, findBarycenter
export rootConstruction, rootConstructionDistinguished, rootConstructionDistinguishedIndices
# Stacky fan
export stackyBlowup, coneRayDecomposition
export findStackyBarycenter, findStackyRayMatrix
# Utility
export slicematrix, rowMinors, convertIncidenceMatrix, convertToIncidence
export coneConvert, getCones, getConesPolymake
# Display
export plot3dSimpCone, plot2dCone, showSimpFan, showSimpStackyFan, coneVectorOrder, plot3dCone
export showFan, showStackyFan

"""
    Structure to store information of a stacky fan - this is a fan together with a dictionary assigning stacky values to each ray.

# Properties:
-> `fan` - the underlying fan, as a polymake object
-> `scalars` - the array of stacky values
-> `stacks` - a dictionary assigning stacky values to each ray.

"""
struct StackyFan
    fan::Polymake.BigObjectAllocated
    stacks::Dict{String, Int64}
    # Constructors for the StackyFan object
    StackyFan(
        fan::Polymake.BigObjectAllocated,
        stacks::Dict{String, Int64}) = new(fan, stacks)
    StackyFan(
        rays::Array{Int64, 2},
        cones::Array{Array{Int64, 1}, 1},
        scalars::Array{Int64, 1}) = makeStackyFan(rays, cones, scalars)
    StackyFan(
        fan::Polymake.BigObjectAllocated,
        scalars::Array{Int64, 1}) = addStackStructure(fan, scalars)
end

## Helper functions

"""
    makeStackyFan(::Array{Int64,2},::Array{Array{Int64,1},1},::Array{Int64,1}))

Function to generate a stacky fan from a matrix representing rays as row vectors, a vector of vectors representing the rays contained in each cone, and a vector of stacky values to be assigned the rays. The second input should be zero-indexed.

# Examples
```jldoctest
julia> makeStackyFan([1 0; 1 1; 1 2],[[0,1],[1,2]],[2,2,2])
[ 2 ,  2 ,  2 ]
```
"""
function makeStackyFan(
    rays::Array{<:Number,2},
    cones::Array{Array{Int64,1},1},
    scalars::Array{Int64,1})

    # Construct a normal fan from the given rays and cones
    fan = fulton.NormalToricVariety(INPUT_RAYS=rays, INPUT_CONES=cones)

    # Construct the dictionary
    stack_rays = mapslices(encode, convert(Array{Int64,2},Array(Polymake.common.primitive(fan.RAYS))), dims=2)
    pairs = map((x,y) -> (x,y), stack_rays, scalars)
    stacks = Dict(pairs)

    return(StackyFan(fan, stacks))
end

"""
    addStackStructure(::Polymake.BigObjectAllocated, ::Array{Int64, 1})

Function to generate a stacky fan from a given fan and a set of scalars.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 1 1; 1 2],INPUT_CONES=[[0,1],[1,2]]);

julia> stackyWeights(addStackStructure(X,[2,2,2]))
[ 2 ,  2 ,  2 ]
```
"""
function addStackStructure(
    fan::Polymake.BigObjectAllocated,
    scalars::Array{Int64, 1})

    # Construct the dictionary
    stack_rays = mapslices(encode, Polymake.common.primitive(fan.RAYS), dims=2)
    pairs = map((x,y) -> (x,y), stack_rays, scalars)
    stacks = Dict(pairs)

    return(StackyFan(fan, stacks))
end

"""
    encode(::Union{Polymake.VectorAllocated{Polymake.Rational},Polymake.VectorAllocated{Polymake.Integer},Vector{Int64}})

    Internal function that converts a vector, representing a ray in the fan,
to a string in order to allow for hashing for the dictionary.

#Examples
```jldoctest
julia> encode([1,0,2,5])
"1,0,2,5"
```
"""
function encode(objects::Union{Polymake.VectorAllocated{Polymake.Rational},Polymake.VectorAllocated{Polymake.Integer},Vector{Int64}})
    return(foldl((x,y) -> string(x, ',', y), objects))
end

"""
    stackyWeights(::StackyFan)

    Returns a list of the stacky weights of the rays of the given stacky fan with the same order as the rays of the fan.

#Examples
```jldoctest
julia> F=makeStackyFan([1 0; 1 1; 1 2; 1 3],[[0,1],[1,2],[2,3]],[1,2,3,4]);

julia> stackyWeights(F)
[ 1 ,  2 ,  3 ,  4 ]
```
"""
function stackyWeights(sf::StackyFan)
    rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(sf.fan.RAYS)))
    #println("Stacky weights ray matrix: $rayMatrix")
    rayList=slicematrix(rayMatrix)
    out=Int64[]
    for ray in rayList
        stack=sf.stacks[encode(ray)]
        push!(out,stack)
    end
    return out
end

## API functions

"""
    getRayStack(::StackyFan, ::Array{Int64, 1})

    Get the scalar associated with a ray in the given stacky fan structure.

# Examples
```jldoctest
julia> F=makeStackyFan([1 0; 1 1; 1 2; 1 3],[[0,1],[1,2],[2,3]],[1,2,3,4]);

julia> getRayStack(F,[1,2])
3
```
"""
function getRayStack(sf::StackyFan, ray::Array{Int64, 1})
    return sf.stacks[encode(ray)]
end

"""
    rootConstruction(::StackyFan, ::Array{Int64, 1})

Given a fan and a set of scalars corresponding to the rays of the fan,
performs a root construction on the fan by multiplying the stack scalars
by the given values.

rootConstruction returns a new StackyFan object, and does not modify the input.

# Examples
```jldoctest

julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2],[1,2,3]]);

julia> SX = StackyFan(X, [2,3,5,7]);

julia> stackyWeights(rootConstruction(SX, [1, 4, 2, 1]))
[ 2 ,  12 ,  10 ,  7 ]
```
"""
function rootConstruction(
    sf::StackyFan,
    scalars::Array{Int64, 1})

    # Multiply the scalars of the fan by the given values
    return StackyFan(sf.fan, stackyWeights(sf) .* scalars)
end

"""
    rootConstructionDistinguishedIndices(::StackyFan, ::Array{Int64, 1}, ::Array{Int64, 1})

    Given a fan, the indices of the distinguished rays in the fan rays (as an incidence matrix), and
a set of scalars corresponding to the rays of the fan, performs a root
construction on the fan by multiplying the stack scalars by the given values.

    rootConstructionDistinguishedIndices returns a new StackyFan object,
and does not modify the input.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2],[1,2,3]]);

julia> SX = StackyFan(X, [2,3,5,7]);

julia> stackyWeights(rootConstructionDistinguishedIndices(SX, [0,1,1,0], [4, 2, 1, 3]))
[ 2 ,  6 ,  5 ,  7 ]
```
"""
function rootConstructionDistinguishedIndices(
    sf::StackyFan,
    distIndices::Array{Int64, 1},
    scalars::Array{Int64, 1})

    numRays = size(sf.fan.RAYS, 1)
    fullScalars = fill(1, numRays)
    for i in 1:numRays
        if distIndices[i]==1 && scalars[i] != 0
            fullScalars[i] = scalars[i]
        end
    end
    # Multiply the scalars of the fan by the given values
    return rootConstruction(sf, fullScalars)
end

"""
    rootConstructionDistinguished(
        ::StackyFan,
        ::Polymake.Matrix{Polymake.Rational},
        ::Array{Int64, 1})

    Given a fan, a set of distinguished rays, and a set of scalars of equal size,
performs a root construction on the fan on the distinguished rays by multiplying
the stack scalars by the given values.

    rootConstructionDistinguished returns a new StackyFan object,
and does not modify the input.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2],[1,2,3]]);

julia> SX = StackyFan(X, [2,3,5,7]);

julia> distinguished = X.RAYS[[2,3],:];

julia> stackyWeights(rootConstructionDistinguished(SX, distinguished, [4, 2]))
[ 2 ,  12 ,  10 ,  7 ]
```
"""
function rootConstructionDistinguished(
    sf::StackyFan,
    rays::Polymake.Matrix{Polymake.Rational},
    scalars::Array{Int64, 1})

    # Check that the rays and scalars are the same size
    #if (size(rays, 1) != length(scalars))
    #    error("Inputs are not of equal size")
    #end

    encoded_rays = mapslices(encode, rays, dims=2)
    # Make a copy of the dictionary
    newStacks = copy(sf.stacks)
    for i in 1:length(encoded_rays)
        ray = encoded_rays[i]
        # Multiply the scalar of the corresponding ray
        newStacks[ray] *= scalars[i]
    end

    # Convert the dictionary to an array of scalars matching the indices
    #newScalars = mapslices(ray -> newStacks[encode(ray)], sf.fan.RAYS, dims=2)
    newScalars = Array{Int64, 1}()
    for i in 1:size(sf.fan.RAYS, 1)
        push!(newScalars, newStacks[encode(sf.fan.RAYS[i,:])])
    end

    return StackyFan(sf.fan, newScalars)
end

"""
    findBarycenter(::Union{AbstractSet,AbstractVector},::Polymake.BigObjectAllocated)

    Takes a normal toric variety X and a set s corresponding to a subset of rays of X, and outputs a polymake vector corresponding to the barycenter of those rays.

# Examples
```jldoctest makeSmoothWithDependencies
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0;1 1; 0 1],INPUT_CONES=[[0,1],[1,2]]);

julia> s=[1,2];

julia> findBarycenter(s,X)
pm::Matrix<pm::Integer>
2 1
```
"""
function findBarycenter(s::Union{AbstractSet,AbstractVector},X::Polymake.BigObjectAllocated)
#     rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(X.RAYS)))
#     rays = rowMinors(rayMatrix, s)
#     dim=size(rays,2)
#     bary=zeros(Int64,dim,1)
#     for i in 1:size(rays,1)
#         bary+=rays[i,:]
#     end
#     bary=Polymake.common.primitive(bary)
#     return vec(bary)
# end

    rays = rowMinors(Array(X.RAYS), s)
    dim=size(rays,2)
    bary=zeros(Polymake.Rational,dim,1)
    for i in 1:size(rays,1)
        bary+=rays[i,:]
    end
    bary=Polymake.common.primitive(transpose(bary))
    return bary
end

"""
    findStackyBarycenter(::Union{AbstractSet,AbstractVector},::StackyFan)

    Takes a stacky fan SX and a set s corresponding to a subset of rays of SX, calculates the 'stacky rays' corresponding to those rays (the rays times their stacky values), and find the barycenter of the stacky rays.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 1 2],INPUT_CONES=[[0,1]]);

julia> F=addStackStructure(X,[2,3]);

julia> findStackyBarycenter([1,2],F)
[ 5 ,  6 ]
```
"""
function findStackyBarycenter(s::Union{AbstractSet,AbstractVector},SX::StackyFan)
    rayMatrix=convert(Array{Int64,2}, Array(Polymake.common.primitive(SX.fan.RAYS)))
    # Multiply the rays by their stacky values
    stackMatrix = diagm(stackyWeights(SX)) * rayMatrix
    rays = rowMinors(stackMatrix, s)
    dim=size(rays,2)
    bary=zeros(Int64,dim,1)
    for i in 1:size(rays,1)
        bary+=rays[i,:]
    end
    return vec(bary)
end

"""
    findStackyRayMatrix(::StackyFan)

    Outputs the ray matrix of the given stacky fan, such that each primitive ray is multiplied by its stacky weights.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 2 3; 1 5],INPUT_CONES=[[0,1],[1,2]]);

julia> F=addStackStructure(X,[1,2,3]);

julia> findStackyRayMatrix(F)
3×2 Matrix{Int64}:
 1   0
 4   6
 3  15
```
"""
function findStackyRayMatrix(sf::StackyFan)
    rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(sf.fan.RAYS)))
    sW=stackyWeights(sf)
    return sW .* rayMatrix
end

"""
    getConeRank(::AbstractMatrix,::AbstractVector)

Takes a matrix and a vector containing indices corresponding to rows of a matrix, and calculates the rank of the matrix consisting only of those rows.

#Examples
```jldoctest makeSmoothWithDependencies
julia> v=[1,2]

julia> M=[0 1; 1 1; 1 0]

julia> getConeRank(v,M)
2
```
"""
function getConeRank(coneRayIndices::AbstractVector, rayMatrix::AbstractMatrix)
    coneRays = rowMinors(rayMatrix,coneRayIndices)
    return rank(Matrix(coneRays))
end

"""
    getDimension(::Polymake.BigObjectAllocated)

Returns the ambient dimension of a normal toric variety.

#Examples
```jldoctest makeSmoothWithDependencies
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 2 0;0 0 1;0 1 0; 1 1 1],INPUT_CONES=[[0,1,2],[0,2,3,4]])

julia> getDimension(X)
3
```
"""
function getDimension(X)
    return size(X.RAYS, 2)
end

"""
    getConeFaces(::Polymake.BigObjectAllocated,::AbstractVector,::AbstractMatrix)

Takes a fan, its ray matrix, and a vector corresponding to one of its cones, and returns a list of maximal strict faces of that cone.

#Examples
```jldoctest makeSmoothWithDependencies
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0; 1 1 0; 1 0 1; 1 1 1],INPUT_CONES=[[0,1,2,3]])

julia> getConeFaces(X,[1,2,3,4],Array(X.RAYS))
[[ 1 ,  2 ], [ 1 ,  3 ], [ 3 ,  4 ], [ 2 ,  4 ]]
```
"""
function getConeFaces(fan::Polymake.BigObjectAllocated,cone::AbstractVector,rayMatrix::AbstractMatrix)
    lattice = fan.HASSE_DIAGRAM
    faces = @Polymake.convert_to Array{Set{Int}} lattice.FACES
    cone_faces=[]
    c = rank(Array(rowMinors(rayMatrix, cone))) - 1
    rank_c_subcone_indices = @Polymake.convert_to Array{Int} Polymake.graph.nodes_of_rank(lattice,c)
    rank_c_subcones = [faces[i + 1] for i in rank_c_subcone_indices]
    for subcone in rank_c_subcones
        new_cone = [i+1 for i in subcone]
        if all((i -> i in cone).(new_cone))
            push!(cone_faces, new_cone)
        end
    end
    return cone_faces
end

"""
    toric_blowup(::Union{AbstractSet,AbstractVector},::Polymake.BigObjectAllocated,::AbstractVector)

    Takes a normal toric variety X, a set s corresponding to a subset of rays of X, and a (optional) polymake vector, v, blow up X at v. If v is not provided, blow up X at the barycenter of s.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 1 2],INPUT_CONES=[[0,1]]);

julia> B=toric_blowup([0,1],X,[1,1]);

julia> B.INPUT_RAYS
pm::Matrix<pm::Rational>
1 0
1 2
1 1

julia> B.INPUT_CONES
pm::IncidenceMatrix<pm::NonSymmetric>
{0 2}
{1 2}
```
"""
function toric_blowup(s, X, v)
    s = [i + 1 for i in s]
    if v==nothing
        v=findBarycenter(s,X)
    end
    if size(v,2)==1
         v=transpose(v)
    end
    coneList = convertIncidenceMatrix(X.MAXIMAL_CONES)
    # Extracting the indices of all the cones in X that contains the set of rays s
    starIndex = findall((t) -> all(((i) -> i in t).(s)), coneList)
    star = [coneList[i] for i in starIndex]
    rayMatrix = X.RAYS

    lattice = X.HASSE_DIAGRAM
    faces = @Polymake.convert_to Array{Set{Int}} lattice.FACES

    # Get all the subcones of X that is contained in one of the cones in star and has the same rank
    clStar = []
    # Iterate over star
    for t in star
        # Get the rank of t
        c = rank(Array(rowMinors(rayMatrix, t))) - 1
        # Get all the subcones of X with rank c
        rank_c_subcone_indices = @Polymake.convert_to Array{Int} Polymake.graph.nodes_of_rank(lattice,c)
        rank_c_subcones = [faces[i + 1] for i in rank_c_subcone_indices]
        # Iterate over rank_c_subcones, and put the cones that is contained in t into clStar
        for cone in rank_c_subcones
            new_cone = [i+1 for i in cone]
            if all((i -> i in t).(new_cone))
                push!(clStar, new_cone)
            end
        end
    end
    clStar = unique(clStar)

    n = size(rayMatrix, 1) + 1
    # Filter out the cones in star from conelist
    coneList = filter(x -> !(x in star), coneList)

    if length(s) == 1
        # If s consists of a single ray, find all the cones in clStar that does not contain s
        newCones = []
        for t in clStar
            if !(s[1] in t)
                push!(newCones, sort(push!(t, s[1])))
            end
        end
        # return newCones plus coneList
        finalCones = [[i - 1 for i in cone] for cone in append!(coneList, newCones)]
        return Polymake.fulton.NormalToricVariety(INPUT_RAYS = Array(X.RAYS), INPUT_CONES = finalCones)
    end
    newCones = []
    for t in clStar
        # Find all the cones in clStar that does not contain at least one ray in s
        # QUESTION: Why seperate this from the one element case? Any won't work with one element list?
        if any(((i) -> !(i in t)).(s))
            push!(newCones, push!(t, n))
        end
    end
    # return newCones plus coneList
    finalRays = vcat((X.RAYS),v)
    finalCones = [[i - 1 for i in cone] for cone in append!(coneList, newCones)]
    return Polymake.fulton.NormalToricVariety(INPUT_RAYS = finalRays, INPUT_CONES = finalCones)
end

"""
    makeSimplicial(::Polymake.BigObjectAllocated)

Takes in a normal toric variety and returns a simplicial toric variety by subdividing (blowing up) the non-simplicial maximal cones.

#Examples
```jldoctest makeSmoothWithDependencies
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2,3]])

julia> X.SIMPLICIAL
false

julia> makeSimplicial(X).SIMPLICIAL
true

julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0 0;0 1 0 0;0 0 1 0;1 -1 1 0; 1 0 -2 0],INPUT_CONES=[[0,1,2,3],[0,4]])

julia> X.SIMPLICIAL
false

julia> makeSimplicial(X).SIMPLICIAL
true
```
"""
function makeSimplicial(X::Polymake.BigObjectAllocated)
    Y = copy(X)
    while (true)
        # If the initial toric variety is simplicial, the program terminates and returns it.
        if Y.SIMPLICIAL==true
            break
        end
        #Maximal cones and ray matrix
        coneList = convertIncidenceMatrix(Y.MAXIMAL_CONES)
        rayMatrix = Y.RAYS
        badCone = nothing
        for i in 1:size(coneList,1)
            cone = coneList[i]
            if (getConeRank(cone, rayMatrix) != size(cone)[1])
                badCone = cone
            end
        end
        if (badCone == nothing)
            # All cones are linearly independent
            break
        else
            # Find the first ray that is contained in more than one orbit
            # and subdivide at that ray, using toricBlowup

            # Get faces (need to replace this)
            edges = getConeFaces(Y,badCone,rayMatrix)
            # Find the first ray that is contained in more than one orbit
            i = 1
            while count(r->(badCone[i] in r), edges) == 1
                i += 1
            end
            # Subdivide at the cone containing just that ray
            badCone=[i-1 for i in badCone]
            Y = toric_blowup(badCone, Y,nothing)
            #Y = toric_blowup([badCone[i]], Y,nothing)
        end
        # Repeat this process until there are no more bad cones
    end
    return Y
end

"""
    makeSmooth(::Polymake.BigObjectAllocated)

Takes in a normal toric variety X and output a new smooth toric variety by iteratively blowing up. In the language of fans, these blowups are achieved by subdividing non-smooth cones.

#Examples
```jldoctest makeSmoothWithDependencies
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[4 -1; 0 1],INPUT_CONES=[[0, 1]])

julia> X.SMOOTH_FAN
false

julia> makeSmooth(X).SMOOTH_FAN
true

julia> makeSmooth(X).INPUT_RAYS
pm::Matrix<pm::Rational>
1 -1/4
0 1
1 0

julia>X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2,3]])

julia>X.SMOOTH_FAN
false

julia>makeSmooth(X).SMOOTH_FAN
true

julia>makeSmooth(X).INPUT_RAYS
pm::Matrix<pm::Rational>
1 0 0
1 1 0
1 0 1
1 1 1
2 1 1

julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;0 1 0;0 0 1;0 -1 -1; -1 0 -1; -2 -1 0],INPUT_CONES=[[0,1,2],[0,1,3],[1,3,4],[1,2,4],[2,4,5],[0,2,5],[0,3,5],[3,4,5]])

julia> X.SMOOTH_FAN
false

julia> makeSmooth(X).SMOOTH_FAN
true
```
"""
function makeSmooth(X::Polymake.BigObjectAllocated)
    Y  = copy(X)
    while(true)
        coneList = convertIncidenceMatrix(Y.MAXIMAL_CONES)
        rayMatrix = Array(Polymake.common.primitive(Y.RAYS))
        k = 1
        # Iterate through the coneList, getting the index of the first cone not smooth
        for coneSet in coneList
            # Getting the number of rays in coneSet
            S=size(coneSet)[1]
            coneRays=rowMinors(rayMatrix,coneSet)
            # Checking whether this cone is smooth
            smoothCheck=Polymake.fan.check_fan_objects(Polymake.polytope.Cone(RAYS=coneRays)).SMOOTH_FAN
            if !smoothCheck
                # If the cone is not simplicial or not smooth, we have found the cone that we need to make smooth
                break
            else
                k+=1
            end
        end
        # At this point, all the cones are smooth. The program terminates.
        if k == size(coneList,1)+1
            break
        end

        # Get the cone that we found to be not smooth
        sigma=coneList[k]
        sigmaRays=slicematrix(rowMinors(rayMatrix,sigma))
        tau=0; tauRays=0; tauCone=0
        # Iterate over the subcones of sigma, finding tau, the smallest one that is not smooth
        for subset in collect(powerset(sigma))
            if size(subset,1) > 1
                S=size(subset)[1]
                subsetRays=rowMinors(rayMatrix,subset)
                subsetCone=Polymake.polytope.Cone(RAYS=subsetRays)
                smoothCheck=Polymake.fan.check_fan_objects(subsetCone).SMOOTH_FAN
                if !smoothCheck
                    tau=subset
                    tauRays=subsetRays
                    tauCone=subsetCone
                    break
                end
            end
        end

        # Getting the Hilbert Basis of tau
        H=slicematrix(Matrix(tauCone.HILBERT_BASIS_GENERATORS[1]))
        rayIndex=0
        # Iterate over the Hilbert Basis, finding the first ray that is not the generator of sigma
        for i in 1:size(H,1)
            if !(H[i] in sigmaRays)
                rayIndex=i
                break
            end
        end
        if rayIndex==0
            # Every Hilbert Basis of tau is a generator of sigma. Make Y simplicial is sufficient to make sigma smooth
            Y=makeSimplicial(Y)
        else
            # blowupRay is not a generator of sigma, blow up tau at blowupRay
            blowupRay=H[rayIndex]
            tau=[i-1 for i in tau]
            Y=toric_blowup(tau,Y,transpose(blowupRay))
        end
    end
    return Y
end

"""
    stackyBlowup(::StackyFan,::Array{Int64,1},::Array{Int64,1})

    Takes a stacky fan sf, a ray excep, and a cone, and subdivides the stacky fan at the given ray. Crucially, the given cone should be the minimal cone containing the exceptional ray. The cone input should be zero-indexed.

#examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 1 2],INPUT_CONES=[[0,1]]);

julia> F=addStackStructure(X,[2,3]);

julia> stackyWeights(stackyBlowup(F,[0,1],[1,1]))
[ 2 ,  1 ,  3 ]
```
"""
function stackyBlowup(sf::StackyFan, cone::Array{Int64,1}, excep::Array{Int64,1})
    # Express the exceptional ray as a scalar multiple of a primitive ray
    # Use this scalar as the stacky weight in the resulting stacky fan
    G=gcd(excep)
    excep=Polymake.common.primitive(excep)

    # Perform toric blowup at the given ray
    blowup = toric_blowup(cone, sf.fan, excep)
    sf.stacks[encode(excep)] = G

    return(StackyFan(blowup, sf.stacks))
end

"""
    getConesPolymake(sf::StackyFan)

    Returns a list of cones of a stacky fan, with the cones represented as polymake objects.

#Examples
```jldoctest
julia> F=makeStackyFan([1 0; 1 1; 1 2; 1 3],[[0,1],[1,2],[2,3]],[1,1,1,1])

julia> getConesPolymake(F)[1].RAYS
pm::Matrix<pm::Rational>
1 0
1 1
```
"""
function getConesPolymake(sf::StackyFan)
    formatted = convertIncidenceMatrix(sf.fan.MAXIMAL_CONES)
    cones = map((x) -> Polymake.polytope.Cone(
        INPUT_RAYS=sf.fan.RAYS[x,:]), formatted)
    return(cones)
end

"""
    slicematrix(::AbstractMatrix{<:Number})

    Take a two-dimensional matrix and output a list of its row vectors.

# Examples
```jldoctest
julia> A=[1 2; 3 4];

julia> slicematrix(A)
[[ 1 ,  2 ], [ 3 ,  4 ]]
```
"""
function slicematrix(A::AbstractMatrix{<:Number})
    return [A[i, :] for i in 1:size(A,1)]
end

"""
    rowMinors(::AbstractMatrix{<:Number},::Union{AbstractSet,AbstractVector})

    Identical to slicematrix, except only returns row vectors indexed by a set S.

# Examples
```jldoctest
julia> A=[1 2 3;4 5 6; 7 8 9];

julia> S=Set([1,3]);

julia> rowMinors(A,S)
2×3 LinearAlgebra.Transpose{Int64,Array{Int64,2}}:
 1  2  3
 7  8  9
```
"""
function rowMinors(A::AbstractMatrix{<:Number},S::Union{AbstractSet,AbstractVector})
    outList=[]
    slices=slicematrix(A)
    for i in 1:size(slices,1)
        if i in S
            append!(outList,[slices[i]])
        end
    end
    return Array(transpose(hcat(outList...)))
end

"""
    convertIncidenceMatrix(::Polymake.IncidenceMatrixAllocated{Polymake.NonSymmetric})

    Takes a Polymake incidence matrix (e.g., the output of X.MAXIMAL_CONES for a toric variety X) and outputs a list of vectors, with each vector recording the indices marked on a given row of the incidence matrix.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0;1 1; 0 1],INPUT_CONES=[[0,1],[1,2]]);

julia> M=X.MAXIMAL_CONES;

julia> convertIncidenceMatrix(M)
[[ 1 ,  2 ], [ 2 ,  3 ]]
```
"""
function convertIncidenceMatrix(A::Polymake.IncidenceMatrixAllocated{Polymake.NonSymmetric})
    A=Array(A)
    dim1=size(A,1)
    dim2=size(A,2)
    out=[]
    for i in 1:dim1
        members=[]
        for j in 1:dim2
            if A[i,j]==true
                append!(members,j)
            end
        end
        append!(out,[members])
    end
    return convert(Array{Array{Int64,1},1}, out)
end

"""
    coneMultiplicity(C::Polymake.BigObjectAllocated)

    Returns the multiplicity of a polyhedral cone (inputted as a Polymake object): here, the multiplicity is defined as the index of the sublattice generated by the rays of the cone, inside the full integer lattice contained in the linear subspace generated by the edges of the cone.

# Examples
```jldoctest

julia> C=Polymake.polytope.Cone(INPUT_RAYS=[1 0; 1 2]);

julia> coneMultiplicity(C)
2
```
"""
function coneMultiplicity(C::Polymake.BigObjectAllocated)
    A=Polymake.common.primitive(C.RAYS)
    M=matrix(ZZ,[fmpz.(y) for y in A])
    SNF=Nemo.snf(M)
    mult=1
    for i in 1:size(SNF,1)
        mult*=SNF[i,i]
    end
    return mult
end

"""
    getMultiplicity(::Array{Int64,1},::Array{Int64,2})

    Same functionality as coneMultiplicity, but calculates the cone rays as a subset of the columns of a ray matrix rather than from a Polymake cone object.

# Examples
```jldoctest
julia> getMultiplicity([1,2],[1 0; 1 2; 1 3])
2
```
"""
function getMultiplicity(cone::Array{Int64,1},rayMatrix::Array{Int64,2})
    A=rowMinors(rayMatrix,cone)
    M=matrix(ZZ,[fmpz.(y) for y in A])
    SNF=Nemo.snf(M)
    mult=1
    for i in 1:size(SNF,1)
        mult*=SNF[i,i]
    end
    return mult
end

"""
    coneConvert(::abstractVector{Int64},::abstractMatrix{Int64})

    Takes a matrix where the columns represent rays, and a list of indices, and forms a Polymake cone object generated by the rays corresponding to those indices.

# Examples
```jldoctest

julia> typeof(coneConvert([1, 2, 4],[1 0 0; 0 1 0; 0 0 1; 1 1 1]))
Polymake.BigObjectAllocated
```
"""
function coneConvert(cone::Array{Int64,1},rayMatrix::Array{Int64,2})
    coneRays=rowMinors(rayMatrix,cone)
    C=Polymake.polytope.Cone(RAYS=coneRays)
    return C
end

"""
    getCones(X::Polymake.BigObjectAllocated)

    Returns all the cones of a fan X as a list of lists, with each interior list containing the indices of the rays generating a given cone.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0; 1 1 0; 1 0 1; 1 1 ],INPUT_CONES=[[0,1,2,3]]);

julia> getCones(X)
[[ 0 ,  1 ,  2 ,  3 ], [ 0 ,  1 ], [ 0 ,  2 ], [ 2 ,  3 ], [ 1 ,  3 ], [ 0 ], [ 1 ], [ 2 ], [ 3 ]]
```
"""
function getCones(X::Polymake.BigObjectAllocated)
    lattice=X.HASSE_DIAGRAM
    faces=@Polymake.convert_to Array{Set{Int}} lattice.FACES
    out=Array{Int64,1}[]
    for i in 2:(size(faces,1)-1)
        newface=Array(@Polymake.convert_to Array{Int} faces[i])
        push!(out,[i+1 for i in newface])
    end
    return out
end

"""
    findFaceContainingRay(::Polymake.BigObjectAllocated,::Array{Int64, 1})

Given a cone and a ray, finds a maximal proper face of the cone containing that ray.

# Examples
```jldoctest
julia> C=Polymake.polytope.Cone(INPUT_RAYS=[1 0 0; 1 0 1; 1 1 0; 1 1 1]);

julia> findFaceContainingRay(C,[1,1,1]).RAYS
pm::Matrix<pm::Rational>
1 0 1
1 1 1
```
"""
function findFaceContainingRay(C::Polymake.BigObjectAllocated,v::Array{Int64, 1})
    rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(C.RAYS)))
    faces=convertIncidenceMatrix(C.RAYS_IN_FACETS)
    if faces==[[]]
        return nothing
    end
    for face in faces
        faceCone=coneConvert(face,rayMatrix)
        if Polymake.polytope.contains(faceCone, v)
            return faceCone
        end
    end
    return nothing
end

"""
    findMinimalCone(::Polymake.BigObjectAllocated,::Array{Int64,1})

Given a fan and a ray, finds the minimal cone fo the fan containing the ray.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0 0; 1 0 1 0; 1 1 0 0; 1 1 1 0; 0 0 0 1],INPUT_CONES=[[0,1,2,3,4]]);

julia> findMinimalCone(X,[1,2,0,0])
[1, 3]
```
"""
function findMinimalCone(X::Polymake.BigObjectAllocated,v::Array{Int64, 1})
    rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(X.RAYS)))
    cones=convertIncidenceMatrix(X.MAXIMAL_CONES)
    #print(cones)
    startCone=nothing
    for cone in cones
        #print(cone)
        polyCone=coneConvert(cone,rayMatrix)
        if Polymake.polytope.contains(polyCone, v)
            startCone=polyCone
        end
    end
    if startCone==nothing
        error("The given ray is not contained in any cone of the fan.")
    end
    currentCone=startCone
    while(true)
        nextCone=findFaceContainingRay(currentCone,v)
        if nextCone==nothing
            break
        end
        currentCone=nextCone
    end
    currentRays=slicematrix(Array(currentCone.RAYS))
    fanRays=slicematrix(rayMatrix)
    indices=findall(x->x in currentRays,fanRays)
    return indices
end

"""
    starSubdivision(::Polymake.BigObjectAllocated,::Array{Int64,1})

    Blows up the input fan at the input ray. A wrapper integrating findMinimalCone and toric_blowup.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0 0; 1 0 1 0; 1 1 0 0; 1 1 1 0; 0 0 0 1],INPUT_CONES=[[0,1,2,3,4]]);

julia> starSubdivision(X,[2,1,0,0]).RAYS
pm::Matrix<pm::Rational>
1 0 0 0
1 0 1 0
0 0 0 1
1 1/2 0 0
1 1 0 0
1 1 1 0
```
"""
function starSubdivision(X::Polymake.BigObjectAllocated, v::Array{Int64, 1})
    minimalCone = findMinimalCone(X, v)
    s = [i - 1 for i in minimalCone]
    v = transpose(v)
    return toric_blowup(s, X, v)
end

"""
    convertToIncidence(v::Array{Int64,1},l::Int64)

Returns a vector of length l, with entries of 1 indexed by v and entries of 0 everywhere else.

# Examples
```jldoctest
julia> convertToIncidence([2,3,5],6)
[ 0 , 1 , 1 , 0 , 1 , 0 ]
```
"""
function convertToIncidence(v::Array{Int64,1},l::Int64)
    out=[]
    for j in 1:l
        if j in v
            append!(out,1)
        else
            append!(out,0)
        end
    end
    return out
end

"""
    interiorPoints(::Polymake.BigObjectAllocated)

    Finds all interior lattice points contained in the fundamental region of a given cone. When multiple interior lattice points lie along the same ray, only the point closest to the origin is returned. Notably,

# Examples
```jldoctest
julia> C=Polymake.polytope.Cone(INPUT_RAYS=[1 2; 2 1]);

julia> interiorPoints(C)
[[ 1 ,  1 ]]
```
"""
function interiorPoints(C::Polymake.BigObjectAllocated)
    rayMatrix=Array(Polymake.common.primitive(C.RAYS))
    l=size(rayMatrix,1)
    dim=size(rayMatrix,2)
    if rank(rayMatrix)<l
        print(rayMatrix)
        error("Input cone is not simplicial.")
    end
    subsets=collect(powerset([1:l;]))
    vertices=[]
    for elt in subsets #vertices of the fundamental region are in correspondence with subsets of the generators of the cone, by summing the generators in a subset to obtain a vertex
        vert=zeros(Polymake.Rational,1,dim)
        for i in 1:l
            if i in elt
                vert+=rayMatrix[[i],:]
            end
        end
        append!(vertices,[vert])
    end
    V=vcat(vertices...)
    VH=hcat(ones(Polymake.Rational,size(V,1)),V)
    P=Polymake.polytope.Polytope(POINTS=VH) #make a Polymake polytope object from the vertices of the fundamental region found in the last step
    if size(P.INTERIOR_LATTICE_POINTS,1)==0
        return nothing
    end
    intPoints=Array(P.INTERIOR_LATTICE_POINTS)[:,2:(dim+1)] #find all the interior lattice points
    validPoints=[]
    #return intPoints
    for i in 1:size(intPoints,1) #throw out all points that are integer multiples of other points
        point=intPoints[i,:]
        if gcd(point)==1
            append!(validPoints,[point])
        end
    end
    return validPoints
end

"""
    coneRayDecomposition(::Array{Int64,1},::Array{Int64,2},::Array{Int64,1},::Array{Int64,1})

    This function takes in a cone (a vector of indices of cone generators in rayMatrix), a ray, and a stacky structure for rayMatrix. It first multiplies all generators of the cone by their stacky values, and then finds an expression for the ray as a sum of these stacky generators. The output is a vector of coefficients of the above representation in terms of the rays in rayMatrix, with zeros as coefficients for all rays not in the given cone.

# Examples
```jldoctest
julia> coneRayDecomposition([1,2,3],[3 5 7; 8 16 9;2 1 3;1 1 1],[2,2,3],[1,1,1,1])
[ 6 ,  5 ,  52 ,  0 ]
```
"""
function coneRayDecomposition(cone,rayMatrix,ray,stack)
    stackMatrix=diagm(stack)*rayMatrix # multiply all rays by stack values
    coneRays=rowMinors(stackMatrix,cone) # find the (stackified) rays in the given cone
    if rank(coneRays)<size(coneRays,1)
        error("The given cone is not simplicial.")
    end
    B=Polymake.common.null_space(hcat(transpose(coneRays),-ray)) # Express the input ray in terms of the stackified cone generators
    N=convert(Array{Int64,1},vec(B))
    if size(N,1)==0
        error("The given ray is not in the span of the cone generators.")
    end
    if N[end]<0 #since the nullspace has arbitrary sign, fix it so the coefficients are all positive
        N*=-1
    end
    pop!(N)
    out=zeros(Int64,size(rayMatrix,1))
    for i in 1:size(N,1)#rewrite the coefficients vector in terms of all the rays in rayMatrix, by padding with zeros when appropriate.
        out[cone[i]]=N[i]
    end
    return out
end

"""
=========VISUALIZATION FUNCTIONALITY==========
"""

using Plots

"""
    plot3dSimpCone(::Array{Array{Int64,1},1})

    Give a list of three 3-dimensional vectors, plots the polygon defined by those vectors and the origin via Plots.jl.

# Examples
```jldoctest
julia> plot3dSimpCone([[1,0,0],[0,1,0],[0,0,1]]);
```
"""
function plot3dSimpCone(L::Array{Array{Int64,1},1})
    l=size(L,1)
    if l==3
        (x1,y1,z1)=L[1]
        (x2,y2,z2)=L[2]
        (x3,y3,z3)=L[3]
        X=[0,x1,0,x2,0,x3,x2,x1,x3]
        Y=[0,y1,0,y2,0,y3,y2,y1,y3]
        Z=[0,z1,0,z2,0,z3,z2,z1,z3]
        plot!(X,Y,Z)
    elseif l==2
        (x1,y1,z1)=L[1]
        (x2,y2,z2)=L[2]
        X=[0,x1,x2,0]
        Y=[0,y1,y2,0]
        Z=[0,z1,z2,0]
        plot!(X,Y,Z)
    elseif l==1
        (x1,y1,z1)=L[1]
        X=[0,x1]
        Y=[0,y1]
        Z=[0,z1]
        plot!(X,Y,Z)
    else
        error("Input cone is empty or not simplicial.")
    end
end

"""
    plot2dCone(::Array{Array{Int64,1},1})

    Given a list of two 2-dimesnional vectors, plots the polygon defined by those vectors and the origin via Plots.jl.

# Examples
```jldoctest
julia> plot2dCone([[1,0],[0,1]]);
```
"""
function plot2dCone(L::Array{Array{Int64,1},1})
    l=size(L,1)
    if l==2
        (x1,y1)=L[1]
        (x2,y2)=L[2]
        S=Shape([(0,0),(x1,y1),(x2,y2)])
        plot!(S)
    elseif l==1
        (x1,y1)=L[1]
        X=[0,x1]
        Y=[0,y1]
        plot!(X,Y)
    else
        error("Input cone is empty or not simplicial.")
    end
end

"""
    showSimpFan(X::Polymake.BigObjectAllocated)

    Plots the input simplicial fan as a collection of polygons or polyhedra defined by the input rays and the origin via Plots.jl.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 2 5],INPUT_CONES=[[0,1]]);

julia> showSimpFan(X);
```
"""
function showSimpFan(X::Polymake.BigObjectAllocated)
    plot()
    rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(X.RAYS)))
    dim=size(rayMatrix,2)
    if !(dim in [2,3])
         println("Fan visualization is only supported in dimensions 2 and 3.")
    end
    maxCones=convertIncidenceMatrix(X.MAXIMAL_CONES)
    for cone in maxCones
        coneRays=slicematrix(rowMinors(rayMatrix,cone))
        if dim==2
            plot2dCone(coneRays)
        elseif dim==3
            plot3dSimpCone(coneRays)
        end
    end
    plot!(legend=false)
end

"""
    showSimpStackyFan(::StackyFan;::Bool=true)

    Extends the functionality of showSimpFan to stacky simplicial fans. If the stackypoints input is set to true, a red dot is shown at the location of each primitive vector multiplied by its stacky weight.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 2 5],INPUT_CONES=[[0,1]]);

julia> F=addStackStructure(X,[1,2]);

julia> showSimpStackyFan(F);
```
"""
function showSimpStackyFan(F::StackyFan;stackypoints::Bool=true)
    X=F.fan
    plot()
    rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(X.RAYS)))
    dim=size(rayMatrix,2)
    if !(dim in [2,3])
         println("Fan visualization is only supported in dimensions 2 and 3.")
    end
    maxCones=convertIncidenceMatrix(X.MAXIMAL_CONES)
    for cone in maxCones
        coneRays=slicematrix(rowMinors(rayMatrix,cone))
        if dim==2
            plot2dCone(coneRays)
        elseif dim==3
            plot3dSimpCone(coneRays)
        end
    end
    if stackypoints==true
        stackyPoints=slicematrix(findStackyRayMatrix(F))
        for point in stackyPoints
            if dim==2
                (x,y)=point
                plot!([x],[y],shape=:circle,color=:red)
            elseif dim==3
                (x,y,z)=point
                plot!([x],[y],[z],shape=:circle,color=:red)
            end
        end
    end
    plot!(legend=false)
end

"""
    coneVectorOrder(::Polymake.BigObjectAllocated)

    Takes a polyhedral cone in 2 or 3 dimensions and outputs a list of its defining rays, arranged in counterclockwise order around the exterior of the cone.

# Examples
```jldoctest
julia> C=Polymake.polytope.Cone(INPUT_RAYS=[1 0 0; 1 0 1; 1 1 0; 1 1 1]);

julia> coneVectorOrder(C)
[[1
, 1, 0], [1, 1, 1], [1, 0, 1], [1, 0, 0]]
```
"""
function coneVectorOrder(C::Polymake.BigObjectAllocated)
    rayList=slicematrix(convert(Array{Int64,2},Array(Polymake.common.primitive(C.RAYS))))
    orderedRayList=Array{Int64,1}[]
    ordering=Array(C.RIF_CYCLIC_NORMAL[1])
    ordering=[i+1 for i in ordering]
    for i in ordering
        push!(orderedRayList,rayList[i])
    end
    return orderedRayList
end

"""
    plot3dCone(::Array{Array{Int64,1},1})

    Plots the 3-dimensional cone defined by the given list of rays. The input rays are assumed to be in counterclockwise or clockwise order around the cone; coneVectorOrder() can be used to obtain this ordering from an arbitrary cone.

# Examples
```jldoctest
julia> plot3dCone([[1,1,0],[1,1,1],[1,0,1],[1,0,0]]);
```
"""
function plot3dCone(L::Array{Array{Int64,1},1})
    l=size(L,1)
    X=Int64[0]
    Y=Int64[0]
    Z=Int64[0]
    for i in 1:l
        push!(X,L[i][1])
        push!(Y,L[i][2])
        push!(Z,L[i][3])
        push!(X,0)
        push!(Y,0)
        push!(Z,0)
    end
    for i in 1:l
        push!(X,L[i][1])
        push!(Y,L[i][2])
        push!(Z,L[i][3])
    end
    push!(X,L[1][1])
    push!(Y,L[1][2])
    push!(Z,L[1][3])
    plot!(X,Y,Z)
end

"""
    showFan(X::Polymake.BigObjectAllocated)

    Plots the input fan as a collection of polygons or polyhedra defined by the input rays and the origin via Plots.jl.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0; 1 0 1; 1 1 0; 1 1 1],INPUT_CONES=[[0,1,2,3]]);

julia> showFan(X);
```
"""
function showFan(X::Polymake.BigObjectAllocated)
    plot()
    rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(X.RAYS)))
    dim=size(rayMatrix,2)
    if !(dim in [2,3])
         println("Fan visualization is only supported in dimensions 2 and 3.")
    end
    maxCones=convertIncidenceMatrix(X.MAXIMAL_CONES)
    for cone in maxCones
        if dim==2
            coneRays=slicematrix(rowMinors(rayMatrix,cone))
            plot2dCone(coneRays)
        elseif dim==3
            C=coneConvert(cone,rayMatrix)
            plot3dCone(coneVectorOrder(C))
        end
    end
    plot!(legend=false)
end

"""
    showStackyFan(::StackyFan;::Bool=true)

    Extends the functionality of showFan to stacky fans. If the stackypoints input is set to true, a red dot is shown at the location of each primitive vector multiplied by its stacky weight.

# Examples
```jldoctest
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0; 1 0 1; 1 1 0; 1 1 1],INPUT_CONES=[[0,1,2,3]]);

julia> F=addStackStructure(X,[1,2,2,5]);

julia> showStackyFan(F);
```
"""
function showStackyFan(F::StackyFan;stackypoints::Bool=true)
    X=F.fan
    plot()
    rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(X.RAYS)))
    dim=size(rayMatrix,2)
    if !(dim in [2,3])
         println("Fan visualization is only supported in dimensions 2 and 3.")
    end
    maxCones=convertIncidenceMatrix(X.MAXIMAL_CONES)
    for cone in maxCones
        if dim==2
            coneRays=slicematrix(rowMinors(rayMatrix,cone))
            plot2dCone(coneRays)
        elseif dim==3
            C=coneConvert(cone,rayMatrix)
            plot3dCone(coneVectorOrder(C))
        end
    end
    if stackypoints==true
        stackyPoints=slicematrix(findStackyRayMatrix(F))
        for point in stackyPoints
            if dim==2
                (x,y)=point
                plot!([x],[y],shape=:circle,color=:red)
            elseif dim==3
                (x,y,z)=point
                plot!([x],[y],[z],shape=:circle,color=:red)
            end
        end
    end
    plot!(legend=false)
end

max_stars_repo_name: nbenabla/destackification

max_issues_repo_name: nbenabla/destackification

max_forks_repo_name: nbenabla/destackification


text:
using Decimals
using Compat.Test

@testset "Normalization" begin

@test Decimal(1, 151100, -4) == Decimal(1, 1511, -2)
@test Decimal(0, 100100, -5) == Decimal(0, 1001, -3)
@test normalize(Decimal(1, 151100, -4)) == Decimal(1, 1511, -2)
@test normalize(Decimal(0, 100100, -5)) == Decimal(0, 1001, -3)
@test parse(Decimal, "3.0") == Decimal(0, 3, 0)
@test parse(Decimal, "3.0") == Decimal(0, 30, -1)
@test parse(Decimal, "3.1400") == Decimal(0, 314, -2)
@test parse(Decimal, "1234") == Decimal(0, 1234, 0)

end

ext: jl


text:
#
# domains
#

export Domain, domain

struct Domain
    handle::nvtxDomainHandle_t

    function Domain(name::String)
        handle = nvtxDomainCreateA(name)
        new(handle)
    end
end

Base.unsafe_convert(::Type{nvtxDomainHandle_t}, dom::Domain) = dom.handle

unsafe_destroy!(dom::Domain) = nvtxDomainDestroy(dom)

function Domain(f::Function, name::String)
    dom = Domain(name)
    f(dom)
    unsafe_destroy!(dom)
end

#
# markers
#

export mark

mark(msg::String) = nvtxMarkA(msg)

#
# ranges
#

export Range, start_range, stop_range, @range

struct Range
    id::nvtxRangeId_t
end

Base.convert(::Type{nvtxRangeId_t}, range::Range) = range.id

"""
    start_range(msg)

Create and start a new range. The range is not automatically stopped, use
[`end_range(::Range)`](@ref) for that.

Use this API if you need overlapping ranges, for scope-based use [`@range`](@ref) instead.
"""
start_range(msg::String) = nvtxRangeStartA(msg)
end_range(r::Range) = nvtxRangeEnd(r)

push_range(msg::String) = nvtxRangePushA(msg)
pop_range() = nvtxRangePop()

"""
    @range "msg" ex
    @range function ... end

Create a new range and execute `ex`. The range is popped automatically afterwards.

See also: [`range`](@ref)
"""
macro range(msg, ex)
    quote
        push_range($(esc(msg)))
        local ret = $(esc(ex))
        pop_range()
        ret
    end
end
macro range(ex)
    def = splitdef(ex)
    def[:body] = quote
        $push_range($(string(def[:name])))
        try
            $(def[:body])
        finally
            $pop_range()
        end
    end
    esc(combinedef(def))
end


text:
using Test
using IBMQClient

@testset "schema" begin
    include("schema.jl")
end

@testset "adapters" begin
    include("adapters.jl")
end


text:
include("clustering.jl")
include("models.jl")

println( models.test_matrix )
println( typeof( models.test_matrix ) )
println( models.test_matrix[CartesianIndex(1, 1)] )
println( clustering.connectivity(models.test_matrix,8.0,true) )
println( clustering.calculate_hierarchy_matrix(models.test_distance_matrix) )

max_stars_repo_name: richardtjornhammar/ingenuous

max_issues_repo_name: richardtjornhammar/ingenuous

max_forks_repo_name: richardtjornhammar/ingenuous


text:
using PkgBenchmark
using LightGraphs

testdatadir = joinpath(dirname(@__FILE__), "..", "test", "testdata")
benchdatadir = joinpath(dirname(@__FILE__), "data")
paramsfile = joinpath(benchdatadir, "params.jld")

println("testdatadir = $testdatadir")
println("paramsfile = $paramsfile")

dg1fn = joinpath(testdatadir, "graph-5k-50k.jgz")

DIGRAPHS = Dict{String,DiGraph}(
    "complete100"   => CompleteDiGraph(100),
    "5000-50000"    => LightGraphs.load(dg1fn)["graph-5000-50000"],
    "path500"       => PathDiGraph(500)
)

GRAPHS = Dict{String,Graph}(
    "complete100"   => CompleteGraph(100),
    "tutte"         => smallgraph(:tutte),
    "path500"       => PathGraph(500),
    "5000-49947"    => Graph(DIGRAPHS["5000-50000"])
)

@benchgroup "LightGraphsBenchmarks" begin
  include("core.jl")
  include("parallel/egonets.jl")
  include("insertions.jl")
  include("edges.jl")
  include("centrality.jl")
  include("connectivity.jl")
  include("max-flow.jl")
  include("traversals.jl")
end


text:
int_rules_4_1_0_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.0.1*(a*sin)^m*(b*trg)^n =#
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := Antiderivative(DeactivateTrig(u, x), x) <-- SimplifyFlag && FunctionOfTrigOfLinearQ(u, x), Antiderivative((~u), ~x) := Antiderivative(DeactivateTrig(u, x), x) <-- FunctionOfTrigOfLinearQ(u, x))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * cos(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => ((~a * sin(~e + ~f * ~x)) ^ (~m + 1) * (~b * cos(~e + ~f * ~x)) ^ (~n + 1)) / (~a * ~b * ~f * (~m + 1)) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (EqQ(~m + ~n + 2, 0) && NeQ(~m, -1))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * cos(~(e') + ~(f') * ~x) ^ ~(n'), ~x) => (1 / (~a * ~f)) * Subst(Antiderivative((~x) ^ ~m * (1 - (~x) ^ 2 / (~a) ^ 2) ^ ((~n - 1) / 2), ~x), ~x, ~a * sin(~e + ~f * ~x)) <-- FreeQ([~a, ~e, ~f, ~m], ~x) && (IntegerQ((~n - 1) / 2) && Not(IntegerQ((~m - 1) / 2) && LtQ(0, ~m, ~n)))
    @apply_utils Antiderivative((~(a') * cos(~(e') + ~(f') * ~x)) ^ ~(m') * sin(~(e') + ~(f') * ~x) ^ ~(n'), ~x) => (-1 / (~a * ~f)) * Subst(Antiderivative((~x) ^ ~m * (1 - (~x) ^ 2 / (~a) ^ 2) ^ ((~n - 1) / 2), ~x), ~x, ~a * cos(~e + ~f * ~x)) <-- FreeQ([~a, ~e, ~f, ~m], ~x) && (IntegerQ((~n - 1) / 2) && Not(IntegerQ((~m - 1) / 2) && (GtQ(~m, 0) && LeQ(~m, ~n))))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * cos(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (-(~a) * (~a * sin(~e + ~f * ~x)) ^ (~m - 1) * (~b * cos(~e + ~f * ~x)) ^ (~n + 1)) / (~b * ~f * (~n + 1)) + (((~a) ^ 2 * (~m - 1)) / ((~b) ^ 2 * (~n + 1))) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ (~m - 2) * (~b * cos(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (GtQ(~m, 1) && (LtQ(~n, -1) && (IntegersQ(2 * ~m, 2 * ~n) || EqQ(~m + ~n, 0))))
    @apply_utils Antiderivative((~(a') * cos(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~a * (~a * cos(~e + ~f * ~x)) ^ (~m - 1) * (~b * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~b * ~f * (~n + 1)) + (((~a) ^ 2 * (~m - 1)) / ((~b) ^ 2 * (~n + 1))) * Antiderivative((~a * cos(~e + ~f * ~x)) ^ (~m - 2) * (~b * sin(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (GtQ(~m, 1) && (LtQ(~n, -1) && (IntegersQ(2 * ~m, 2 * ~n) || EqQ(~m + ~n, 0))))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * cos(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (-(~a) * (~b * cos(~e + ~f * ~x)) ^ (~n + 1) * (~a * sin(~e + ~f * ~x)) ^ (~m - 1)) / (~b * ~f * (~m + ~n)) + (((~a) ^ 2 * (~m - 1)) / (~m + ~n)) * Antiderivative((~b * cos(~e + ~f * ~x)) ^ ~n * (~a * sin(~e + ~f * ~x)) ^ (~m - 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && (GtQ(~m, 1) && (NeQ(~m + ~n, 0) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * cos(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~a * (~b * sin(~e + ~f * ~x)) ^ (~n + 1) * (~a * cos(~e + ~f * ~x)) ^ (~m - 1)) / (~b * ~f * (~m + ~n)) + (((~a) ^ 2 * (~m - 1)) / (~m + ~n)) * Antiderivative((~b * sin(~e + ~f * ~x)) ^ ~n * (~a * cos(~e + ~f * ~x)) ^ (~m - 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && (GtQ(~m, 1) && (NeQ(~m + ~n, 0) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * cos(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~b * cos(~e + ~f * ~x)) ^ (~n + 1) * (~a * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~a * ~b * ~f * (~m + 1)) + ((~m + ~n + 2) / ((~a) ^ 2 * (~m + 1))) * Antiderivative((~b * cos(~e + ~f * ~x)) ^ ~n * (~a * sin(~e + ~f * ~x)) ^ (~m + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && (LtQ(~m, -1) && IntegersQ(2 * ~m, 2 * ~n))
    @apply_utils Antiderivative((~(a') * cos(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (-((~b * sin(~e + ~f * ~x)) ^ (~n + 1)) * (~a * cos(~e + ~f * ~x)) ^ (~m + 1)) / (~a * ~b * ~f * (~m + 1)) + ((~m + ~n + 2) / ((~a) ^ 2 * (~m + 1))) * Antiderivative((~b * sin(~e + ~f * ~x)) ^ ~n * (~a * cos(~e + ~f * ~x)) ^ (~m + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && (LtQ(~m, -1) && IntegersQ(2 * ~m, 2 * ~n))
    @apply_utils Antiderivative(sqrt(~(a') * sin(~(e') + ~(f') * ~x)) * sqrt(~(b') * cos(~(e') + ~(f') * ~x)), ~x) => ((sqrt(~a * sin(~e + ~f * ~x)) * sqrt(~b * cos(~e + ~f * ~x))) / sqrt(sin(2 * ~e + 2 * ~f * ~x))) * Antiderivative(sqrt(sin(2 * ~e + 2 * ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x)
    @apply_utils Antiderivative(1 / (sqrt(~(a') * sin(~(e') + ~(f') * ~x)) * sqrt(~(b') * cos(~(e') + ~(f') * ~x))), ~x) => (sqrt(sin(2 * ~e + 2 * ~f * ~x)) / (sqrt(~a * sin(~e + ~f * ~x)) * sqrt(~b * cos(~e + ~f * ~x)))) * Antiderivative(1 / sqrt(sin(2 * ~e + 2 * ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x)
    #= Antiderivative(((~a')*sin((~e')+(~f')*(~x)))^(~m)*((~b')*cos((~e')+(~f')*(~x)))^(~n),~x) :=   (a*sin(e+f*x))^m*(b*cos(e+f*x))^n/(a*tan(e+f*x))^m*Antiderivative((a*tan(e+f*x) )^m,x) <-- FreeQ([a,b,e,f,m,n],x) && EqQ(m+n,0) =#
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * cos(~(e') + ~(f') * ~x)) ^ ~n, ~x) => With([k = Denominator(~m)], ((k * ~a * ~b) / ~f) * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) / ((~a) ^ 2 + (~b) ^ 2 * (~x) ^ (2k)), ~x), ~x, (~a * sin(~e + ~f * ~x)) ^ (1 / k) / (~b * cos(~e + ~f * ~x)) ^ (1 / k))) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (EqQ(~m + ~n, 0) && (GtQ(~m, 0) && LtQ(~m, 1)))
    @apply_utils Antiderivative((~(a') * cos(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => With([k = Denominator(~m)], ((-k * ~a * ~b) / ~f) * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) / ((~a) ^ 2 + (~b) ^ 2 * (~x) ^ (2k)), ~x), ~x, (~a * cos(~e + ~f * ~x)) ^ (1 / k) / (~b * sin(~e + ~f * ~x)) ^ (1 / k))) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (EqQ(~m + ~n, 0) && (GtQ(~m, 0) && LtQ(~m, 1)))
    #= Antiderivative(((~a')*sin((~e')+(~f')*(~x)))^(~m)*((~b')*cos((~e')+(~f')*(~x)))^(~n),~x) :=   b^(2*IntPart((n-1)/2)+1)*(b*cos(e+f*x))^(2*FracPart((n-1)/2))/(a*f*( cos(e+f*x)^2)^FracPart((n-1)/2))* Subst(Antiderivative(x^m*(1-x^2/a^2)^((n-1)/2),x),x,a*sin(e+f*x)) <--  FreeQ([a,b,e,f,m,n],x) && (RationalQ(n) || Not(RationalQ(m)) &&  (EqQ(b,1) || NeQ(a,1))) =#
    #= Antiderivative(((~a')*cos((~e')+(~f')*(~x)))^(~m)*((~b')*sin((~e')+(~f')*(~x)))^(~n),~x) :=   -b^(2*IntPart((n-1)/2)+1)*(b*sin(e+f*x))^(2*FracPart((n-1)/2))/(a*f* (sin(e+f*x)^2)^FracPart((n-1)/2))* Subst(Antiderivative(x^m*(1-x^2/a^2)^((n-1)/2),x),x,a*cos(e+f*x)) <--  FreeQ([a,b,e,f,m,n],x) =#
    @apply_utils Antiderivative((~(a') * cos(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((-((~b) ^ (2 * IntPart((~n - 1) / 2) + 1)) * (~b * sin(~e + ~f * ~x)) ^ (2 * FracPart((~n - 1) / 2)) * (~a * cos(~e + ~f * ~x)) ^ (~m + 1)) / (~a * ~f * (~m + 1) * (sin(~e + ~f * ~x) ^ 2) ^ FracPart((~n - 1) / 2))) * HypergeometricFunctions._₂F₁((1 + ~m) / 2, (1 - ~n) / 2, (3 + ~m) / 2, cos(~e + ~f * ~x) ^ 2) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && SimplerQ(~n, ~m)
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * cos(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((~b) ^ (2 * IntPart((~n - 1) / 2) + 1) * (~b * cos(~e + ~f * ~x)) ^ (2 * FracPart((~n - 1) / 2)) * (~a * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~a * ~f * (~m + 1) * (cos(~e + ~f * ~x) ^ 2) ^ FracPart((~n - 1) / 2))) * HypergeometricFunctions._₂F₁((1 + ~m) / 2, (1 - ~n) / 2, (3 + ~m) / 2, sin(~e + ~f * ~x) ^ 2) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~b * (~a * sin(~e + ~f * ~x)) ^ (~m + 1) * (~b * sec(~e + ~f * ~x)) ^ (~n - 1)) / (~a * ~f * (~m + 1)) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (EqQ((~m - ~n) + 2, 0) && NeQ(~m, -1))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~a * ~b * (~a * sin(~e + ~f * ~x)) ^ (~m - 1) * (~b * sec(~e + ~f * ~x)) ^ (~n - 1)) / (~f * (~n - 1)) - (((~a) ^ 2 * (~b) ^ 2 * (~m - 1)) / (~n - 1)) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ (~m - 2) * (~b * sec(~e + ~f * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (GtQ(~n, 1) && (GtQ(~m, 1) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~a * sin(~e + ~f * ~x)) ^ (~m + 1) * (~b * sec(~e + ~f * ~x)) ^ (~n + 1)) / (~a * ~b * ~f * (~m - ~n)) - ((~n + 1) / ((~b) ^ 2 * (~m - ~n))) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ ~m * (~b * sec(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (GtQ(~n, 1) && IntegersQ(2 * ~m, 2 * ~n))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~a * sin(~e + ~f * ~x)) ^ (~m + 1) * (~b * sec(~e + ~f * ~x)) ^ (~n + 1)) / (~a * ~b * ~f * (~m + 1)) - ((~n + 1) / ((~a) ^ 2 * (~b) ^ 2 * (~m + 1))) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ (~m + 2) * (~b * sec(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (LtQ(~n, -1) && (LtQ(~m, -1) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~a * sin(~e + ~f * ~x)) ^ (~m + 1) * (~b * sec(~e + ~f * ~x)) ^ (~n + 1)) / (~a * ~b * ~f * (~m - ~n)) - ((~n + 1) / ((~b) ^ 2 * (~m - ~n))) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ ~m * (~b * sec(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (LtQ(~n, -1) && (NeQ(~m - ~n, 0) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (-(~a) * ~b * (~a * sin(~e + ~f * ~x)) ^ (~m - 1) * (~b * sec(~e + ~f * ~x)) ^ (~n - 1)) / (~f * (~m - ~n)) + (((~a) ^ 2 * (~m - 1)) / (~m - ~n)) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ (~m - 2) * (~b * sec(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && (GtQ(~m, 1) && (NeQ(~m - ~n, 0) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~b * (~a * sin(~e + ~f * ~x)) ^ (~m + 1) * (~b * sec(~e + ~f * ~x)) ^ (~n - 1)) / (~a * ~f * (~m + 1)) + (((~m - ~n) + 2) / ((~a) ^ 2 * (~m + 1))) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ (~m + 2) * (~b * sec(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && (LtQ(~m, -1) && IntegersQ(2 * ~m, 2 * ~n))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~b * cos(~e + ~f * ~x)) ^ ~n * (~b * sec(~e + ~f * ~x)) ^ ~n * Antiderivative((~a * sin(~e + ~f * ~x)) ^ ~m / (~b * cos(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (IntegerQ(~m - 1 / 2) && IntegerQ(~n - 1 / 2))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (1 / (~b) ^ 2) * (~b * cos(~e + ~f * ~x)) ^ (~n + 1) * (~b * sec(~e + ~f * ~x)) ^ (~n + 1) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ ~m / (~b * cos(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (Not(IntegerQ(~m)) && (Not(IntegerQ(~n)) && LtQ(~n, 1)))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~b) ^ 2 * (~b * cos(~e + ~f * ~x)) ^ (~n - 1) * (~b * sec(~e + ~f * ~x)) ^ (~n - 1) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ ~m / (~b * cos(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (Not(IntegerQ(~m)) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~a * ~b) ^ IntPart(~n) * (~a * sin(~e + ~f * ~x)) ^ FracPart(~n) * (~b * csc(~e + ~f * ~x)) ^ FracPart(~n) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ (~m - ~n), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (Not(IntegerQ(~m)) && Not(IntegerQ(~n)))
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.0.1 (a sin)^m (b trg)^n.jl


text:
mutable struct MaximinObjectiveCore
    points::Vector{Float64}
    IIF::Matrix{Float64}

    MaximinObjectiveCore() = new(Vector{Float64}(), Matrix{Float64}(undef, 0, 0))
end


text:
excludedim(a, dim)= (a[1:dim-1],a[dim+1:end])

# takes SP activations in 1 timestep
sparseness(aᵢ)= count(aᵢ)/length(aᵢ)

"""
`entropy(aᵢₜ; tdim=2)` given a history of SP activations in time returns the SP's total entropy.
 - `t` is assumed to be the 2nd dimension. If `t` is another dimension, set `tdim`
"""
entropy(aᵢₜ; tdim=2, trials=10)= begin
  # should take a vector in time foreach minicolumn
  entropy_lim(aₜ)= aₜ ≈ 0 ? 0.0 :
    aₜ ≈ 1 ? 0.0 : entropyf(aₜ)
  # aₜᵢ:= (t,i,j,...)
  # @view aₜᵢ[:,i,j,...]
  c= Iterators.product((excludedim(size(aᵢₜ), tdim).|> CartesianIndices)...)
  mapreduce(+, c) do (pre,post)
    entropy_lim(sparseness(@view aᵢₜ[pre,:,post]))
  end / length(c)
end
entropyf(x)= -x*log2(x) -(1-x)*log2(1-x)

entropy(sp::SpatialPooler, data; trials=10)= begin
  sampleSize= floor(Int,size(data,2)/trials)
  sampleIdx= randperm(size(data,2))
  trial(i)= begin
    zset= data[:,sampleIdx[(i-1)*sampleSize+1 : i*sampleSize]]
    mapslices(z-> sp(z)|> vec, zset, dims=1)|> entropy
  end
  map(trial, 1:trials)
end

"""
`addnoise(sdr,k)` flips `k ∈ [0,1]`- fraction of random active bits of the SDR and the same amount of inactive,
adding k-level noise while maintaining sparsity.
"""
addnoise(sdr,k)= begin
  k<0 ? k=0 : ( k>1 ? k=1 : nothing)
  nFlip= min( round(Int,k*count(sdr)), length(sdr) - count(sdr) )
  act= (sdr|> findall|> shuffle)[1:nFlip]
  inact= (.~sdr|> findall|> shuffle)[1:nFlip]
  flip!(x,i)= begin
    @views x[i].= .~x[i]
    x
  end
  @> sdr copy flip!(act) flip!(inact)
end

"""
`(metric, changeCurve) = noiserobustness(sp::SpatialPooler, data; trials=15)` runs repeated experiments on the SP
with increasingly distorted stimuli to gauge the noise robustness.
"""
noiserobustness(sp::SpatialPooler, data; trials=10)= begin
  @unpack szᵢₙ = sp.params
  k_step= 19;

  zⁿ(z,k)= addnoise(z,k)

  # The SP activation is stochastic: sp(z) != sp(z). Estimate how many bits differ on average for the same input
  ref= sp(data[:,1])
  spTieNoiseEstimate= 0 #mapreduce(_ -> count(ref) - count(sp(data[:,1]) .& ref), +, 1:100) / 100 / count(ref)

  k_range= 0:1/k_step:1
  sampleSize= floor(Int,size(data,2)/trials)
  changeWithNoise= zeros(length(k_range),trials)
  sampleIdx= randperm(size(data,2))
  trial(i)= begin
    zset= data[:,sampleIdx[(i-1)*sampleSize+1 : i*sampleSize]]
    # ∫{0..1} f(k)dk
    mapreduce(+, enumerate(k_range)) do (kᵢ,k)
      changeWithNoise[kᵢ,i]= mapreduce(+, 1:sampleSize) do s
        z= zset[:,s]; zⁿ= addnoise(z,k)
        a= sp(z); aⁿ= sp(zⁿ)
        count(a .& aⁿ) / count(a)
      end / sampleSize
      changeWithNoise[kᵢ,i]
    end / (k_step+1)
  end
  result= map(trial, 1:trials)
  (result, mapslices(mean,changeWithNoise,dims=[2]), spTieNoiseEstimate, k_range)
end

"""
`stabilitySP(aᵢ,paᵢ)` given a vector `aᵢ` of SP activations for a set of test inputs at
time j, and the activations on the same set of test inputs at a previous time k `paᵢ`,
shows a metric of how much the SP's response has changed from training.
`aᵢ` dims: {activation, samples}
"""
stabilitySP(aᵢ,paᵢ)=
  mapreduce(+, 1:size(aᵢ,2)) do i
    @views count(aᵢ[:,i] .& paᵢ[:,i]) / count(paᵢ[:,i])
  end / size(aᵢ,2)

max_stars_repo_name: Oblynx/HTMexperiments

max_issues_repo_name: Oblynx/HTMexperiments

max_forks_repo_name: Oblynx/HTMexperiments


text:
import Base: getindex, show

@doc raw"""
    BasisSet

Object holding a set of BasisFunction objects associated with an array of atoms.

# Fields

| Name        | Type                               |   Description |
|:------------|:-----------------------------------|:-----------------------------------------------------------|
|`atoms`      | `Vector{Atom}`                     | An array of `Molecules.Atom` objects  |
|`name`       | `String`                           | String holding the basis set name  |
|`basis`      | `Vector{Vector{BasisFunction}}`    | An Array of arrays with BasisFunction          |
|`natoms`     | `Int32`                            | Number of atoms in the BasisSet |
|`nbas`       | `Int32`                            | Number of basis functions (Note, not equal the number of BasisFunction objects) |
|`nshells`    | `Int32`                            | Number of shells, i.e. BasisFunction objects |
|`lc_atoms`   | `Array{Int32,1}`                   | Integer array maping data to libcint |
|`lc_bas`     | `Array{Int32,1}`                   | Integer array maping data to libcint |
|`lc_env`     | `Array{Float64,1}`                 | Float64 array maping data to libcint |

# Example

Build a basis set from default options
```julia
julia> water = Molecules.parse_string(
"O        1.2091536548      1.7664118189     -0.0171613972
 H        2.1984800075      1.7977100627      0.0121161719
 H        0.9197881882      2.4580185570      0.6297938832"
)
julia> bset = BasisSet("sto-3g", water)
sto-3g Basis Set
Number of shells: 5
Number of basis:  7

O: 1s 2s 1p
H: 1s
H: 1s
```
The BasisSet object can be accessed as two-dimensional array.
```julia
julia> bset[1,2] # Show the second basis for the first atom (O 2s)
S shell with 1 basis built from 3 primitive gaussians

χ₀₀  =    0.8486970052⋅Y₀₀⋅exp(-5.033151319⋅r²)
     +    1.1352008076⋅Y₀₀⋅exp(-1.169596125⋅r²)
     +    0.8567529838⋅Y₀₀⋅exp(-0.38038896⋅r²)
```
You can also create your own crazy mix!
Lets create one S and one P basis functions for H
```julia
julia> h2 = Molecules.parse_string(
   "H 0.0 0.0 0.0
    H 0.0 0.0 0.7"
)
julia> s = BasisFunction(0, [0.5215367271], [0.122])
julia> p = BasisFunction(1, [1.9584045349], [0.727])
```
The basis set is constructed with an array of atoms (Vector{Atom}) and a corresponding array of Vector{BasisFunction}
holding all basis functions for that particular atom. In this example, we consider an unequal treatment for the
two atoms in the H₂ molecule.
```
julia> shells = [
    [s],  # A s function on the first hydrogen
    [s,p] # One s and one p function on the second hydrogen
]
julia> BasisSet(h2, "UnequalHydrogens", shells)
UnequalHydrogens Basis Set
Number of shells: 3
Number of basis:  5

H: 1s
H: 1s 1p
```
"""
struct BasisSet
    name::String
    atoms::Vector{Atom}
    basis::Vector{Vector{BasisFunction}}
    natoms::Cint
    nbas::Cint
    nshells::Cint
    lc_atoms::Array{Cint,1}
    lc_bas::Array{Cint, 1}
    lc_env::Array{Cdouble,1}
end

function BasisSet(name::String, atoms::Vector{T}) where T <: Atom

    basis = Vector{BasisFunction}[]

    # Cache entries so repeated atoms don't need to be looked up all the time
    cache = Dict()
    for i in eachindex(atoms)
        S = symbol(atoms[i])
        if S in keys(cache)
            push!(basis, cache[S])
        else
            push!(basis, read_basisset(name, S))
            cache[S] = basis[end]
        end
    end
    BasisSet(name, atoms, basis)
end

function BasisSet(name::String, atoms::Vector{T}, basis::Vector{Vector{BasisFunction}}) where T <: Atom

    length(atoms) != length(basis) ? throw(ArgumentError("Error combining atoms and basis functions")) : nothing

    ATM_SLOTS = 6
    BAS_SLOTS = 8

    natm = length(atoms)
    nbas = 0
    nshells = 0
    nexps = 0
    nprims = 0

    for i in eachindex(atoms)
        for b in basis[i]
            nshells += 1
            nbas += 2*b.l + 1
            nexps += length(b.exp)
            nprims += length(b.coef)
        end
    end

    lc_atm = zeros(Cint, natm*ATM_SLOTS)
    lc_bas = zeros(Cint, nshells*BAS_SLOTS)
    env = zeros(Cdouble, 3*natm+nexps+nprims)

    # Prepare the lc_atom input
    off = 0
    ib = 0
    for i = eachindex(atoms)
        A = atoms[i]
        # lc_atom has ATM_SLOTS (6) "spaces" for each atom
        # The first one (Z_INDEX) is the atomic number
        lc_atm[1 + ATM_SLOTS*(i-1)] = Cint(A.Z)
        # The second one is the env index address for xyz
        lc_atm[2 + ATM_SLOTS*(i-1)] = off
        env[off+1:off+3] .= A.xyz ./ Molecules.bohr_to_angstrom
        off += 3
        # The remaining 4 slots are zero.

        # Prepare the lc_bas input
        for j = eachindex(basis[i])
            B = basis[i][j]
            Ne = length(B.exp)
            Nc = length(B.coef)
            # lc_bas has BAS_SLOTS for each basis set
            # The first one is the index of the atom starting from 0
            lc_bas[1 + BAS_SLOTS*ib] = i-1
            # The second one is the angular momentum
            lc_bas[2 + BAS_SLOTS*ib] = B.l
            # The third is the number of primitive functions
            lc_bas[3 + BAS_SLOTS*ib] = Nc
            # The fourth is the number of contracted functions
            lc_bas[4 + BAS_SLOTS*ib] = 1
            # The fifth is a κ parameter
            lc_bas[5 + BAS_SLOTS*ib] = 0
            # Sixth is the env index address for exponents
            lc_bas[6 + BAS_SLOTS*ib] = off
            env[off+1:off+Ne] .= B.exp
            off += Ne
            # Seventh is the env index address for contraction coeff
            lc_bas[7 + BAS_SLOTS*ib] = off
            env[off+1:off+Nc] .= B.coef
            off += Nc
            # Eigth, nothing
            ib += 1
        end
    end

    return BasisSet(name, atoms, basis, natm, nbas, nshells, lc_atm, lc_bas, env)
end

"""
    GaussianBasis.gto_norm(n::Signed, a::AbstractFloat)

Function that returns the normalization factor for a gaussian basis function rⁿ⋅exp(-a⋅r²)
"""
function gto_norm(n, a)
   # normalization factor of function rⁿ exp(-ar²)
    s = 2^(2n+3) * factorial(n+1) * (2a)^(n+1.5) / (factorial(2n+2) * √π)
    return √s
end

"""
    GaussianBasis.normalize_basisfunction!(B::BasisFunction)

Modify the BasisFunction object by normalizing each primitive.
"""
function normalize_basisfunction!(B::BasisFunction)
    for i = eachindex(B.coef)
        B.coef[i] *= gto_norm(B.l, B.exp[i])
    end
end

function getindex(B::BasisSet, N::Int)
    return B.basis[N]
end

function getindex(B::BasisSet, i::Int, j::Int)
    return B.basis[i][j]
end

function string_repr(B::BasisFunction)
    # Generate Unicode symbol for sub number
    l_sub = Char(0x2080 + B.l)

    # Unicode for superscript is a bit messier, so gotta use control flow
    l_sup = B.l == 1 ? Char(0x00B9) :
            B.l in [2,3] ? Char(0x00B1 + B.l) :
            Char(0x2070 + B.l)

    nbas = 2*B.l + 1
    mvals = collect(-B.l:B.l)
    nprim = length(B.exp)

    # Reverse Dict(symbol=>num) to get Symbols from B.l
    Lsymbol = Dict(value => key for (key, value) in AMDict)[B.l]
    out = "$(Lsymbol) shell with $nbas basis built from $nprim primitive gaussians\n\n"
    for m in mvals
        # Add sub minus sign (0x208B) if necessary
        m_sub = m < 0 ? Char(0x208B)*Char(0x2080 - m) : Char(0x2080 + m)
        out *= format("{:<4s} = ","χ$(l_sub)$m_sub")
        for i in eachindex(B.coef)

            if i > 1
                out *= B.coef[i] > 0 ? "\n     + " : "\n     - "
            end

            #out *= "$(abs(B.coef[i]))⋅Y$(l_sub)$m_sub"
            out *= format("{:>15.10f}⋅Y$(l_sub)$m_sub", abs(B.coef[i]))

            if B.l != 0
                out *= "⋅r$l_sup"
            end

            out *= "⋅exp(-$(B.exp[i])⋅r²)"
        end
        out *="\n\n"
    end
    return strip(out)
end

function string_repr(B::BasisSet)
    out  =  "$(B.name) Basis Set\n"
    out *= "Number of shells: $(B.nshells)\n"
    out *= "Number of basis:  $(B.nbas)\n\n"

    l_to_symbol = Dict(
        0 => "s",
        1 => "p",
        2 => "d",
        3 => "f",
        4 => "g",
        5 => "h",
        6 => "i",
    )

    for i in eachindex(B.atoms)
        A = B.atoms[i]
        # Count how many times s,p,d appears for numbering
        count = zeros(Int16, 7)
        out *= "$(symbol(A)): "
        for bfs in B.basis[i]
            L = bfs.l
            count[L+1] += 1
            out *= "$(count[L+1])$(l_to_symbol[L]) "
        end
        out *="\n"
    end

    return strip(out)
end

# Pretty printing
function show(io::IO, ::MIME"text/plain", X::T) where T<:Union{BasisFunction, BasisSet}
    print(io, string_repr(X))
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Connectivity{PL,N}

A connectivity list of `N` indices representing a [`Polytope`](@ref)
of type `PL`. Indices are taken from a global vector of [`Point`](@ref).

Connectivity objects are constructed with the [`connect`](@ref) function.
"""
struct Connectivity{PL<:Polytope,N}
  list::NTuple{N,Int}
end

"""
    paramdim(connectivity)

Return the parametric dimension of the `connectivity`.
"""
paramdim(::Type{Connectivity{PL,N}}) where {PL,N} = paramdim(PL)
paramdim(c::Connectivity) = paramdim(typeof(c))

"""
    indices(connectivity)

Return the list of indices of the `connectivity`.
"""
indices(c::Connectivity) = c.list

"""
    connect(list, PL)

Connect a `list` of indices from a global vector of [`Point`](@ref)
into a [`Polytope`](@ref) of type `PL`.

## Example
```
Δ = connect((1,2,3), Triangle)
```
"""
connect(list::Tuple, PL::Type{<:Polytope}) = Connectivity{PL,length(list)}(list)

"""
    materialize(connec, points)

Materialize a face using the `connec` list and a global vector of `points`.
"""
function materialize(connec::Connectivity{PL},
                     points::AbstractVector{P}) where {PL<:Polytope,P<:Point}
  PL(view(points, SVector(connec.list...)))
end

function Base.show(io::IO, c::Connectivity{PL}) where {PL}
  print(io, "$PL$(c.list)")
end


text:
using BitBasis
using Test
using Aqua
Aqua.test_all(BitBasis)

const ⊗ = kron

@testset "test bit literal" begin
    include("bit_str.jl")
end

@testset "test operations" begin
    include("operations.jl")
end

@testset "test reorder" begin
    include("reorder.jl")
end

@testset "test utils" begin
    include("utils.jl")
end

@testset "test iterate control" begin
    include("iterate_control.jl")
end


text:
#=
function edit(style::Style, m::Module, custumization::Dict{Symbol,Any})

end

function edit(style::Style, f::Function, custumization::Dict{Symbol,Any})

end

function edit(style::Style, t::Type, custumization::Dict{Symbol,Any})

end
=#


text:
# This is a Julia script to
# plot the simulated orbits

"Check if the value of the option is the last argument"
function check_last(i)
    if i + 1 == length(ARGS)
        println("The last argument should be the input directory.")
        exit(1)
    end
end

# Define default values for optional arguments
POSTFIX = ""

# Parse the options
for i in eachindex(ARGS)
    # A number of integration iterations used
    if ARGS[i] == "-n"
        check_last(i)
        try
            global N = parse(UInt, ARGS[i+1])
        catch
            println("Couldn't parse the value of the `-n` argument.")
            exit(1)
        end
    end
    # A number of simulation iterations
    if ARGS[i] == "-s"
        check_last(i)
        try
            global S = parse(UInt, ARGS[i+1])
        catch
            println("Couldn't parse the value of the `-h` argument.")
            exit(1)
        end
    end
    # A postfix for the names of output files
    if ARGS[i] == "--postfix"
        check_last(i)
        try
            global POSTFIX = " ($(ARGS[i+1]))"
        catch
            println("Couldn't parse the value of the `--postfix` argument.")
            exit(1)
        end
    end
end

# Check for required arguments (exc)
if !@isdefined(N) || !@isdefined(S) || length(ARGS) == 0
    println("""
        Usage:
        julia --project=. scripts/simulations.jl -n <N> -s <S> <input> [--postfix <POSTFIX>]"""
    )
    exit(1)
end

# Define the input directory
INPUT_DIR = ARGS[end]

println('\n', " "^4, "> Loading the packages...")

using LaTeXStrings
using Plots
using Statistics

# Use the GR backend for plots
gr()

# Change some of the default parameters for plots
default(fontfamily = "Computer Modern", dpi = 300, legend = :topright)

# Define the number of minor ticks for heatmaps
minorticks = 8

# Define the paths to output directories
CURRENT_DIR = @__DIR__
ROOT_DIR = basename(CURRENT_DIR) == "scripts" ? dirname(CURRENT_DIR) : CURRENT_DIR
PLOTS_DIR = joinpath(ROOT_DIR, "plots")
SIMULATIONS_DIR = joinpath(PLOTS_DIR, "simulations")
OUTPUT_DIR = joinpath(SIMULATIONS_DIR, basename(INPUT_DIR))

# Make sure the needed directories exist
mkpath(OUTPUT_DIR)

"Plot a heatmap from the coordinates data"
function plot_heatmap(name, object_dir, x, y, xlabel, ylabel, plane)
    # Plot the orbit
    p = plot(
        x[1:100:(N+1)],
        y[1:100:(N+1)];
        label = "",
        title = name,
        xlabel,
        ylabel,
        size = (400, 400),
        minorticks,
        minorgrid = true
    )

    # Define the ticks step
    h₁ = abs(xticks(p)[1][1][2] - xticks(p)[1][1][1]) / minorticks
    h₂ = abs(yticks(p)[1][1][2] - yticks(p)[1][1][1]) / minorticks

    # Define the ranges of the heatmap matrix
    r₁ = range(
        (xticks(p)[1][1][1] - h₁ * floor((xticks(p)[1][1][1] - xlims(p)[1]) / h₁)) - 5 * h₁,
        (xticks(p)[1][1][end] + h₁ * floor((xlims(p)[2] - xticks(p)[1][1][end]) / h₁)) + 5 * h₁;
        step = h₁
    )
    r₂ = range(
        (yticks(p)[1][1][1] - h₂ * floor((yticks(p)[1][1][1] - ylims(p)[1]) / h₂)) - 5 * h₂,
        (yticks(p)[1][1][end] + h₂ * floor((ylims(p)[2] - yticks(p)[1][1][end]) / h₂)) + 5 * h₂;
        step = h₂
    )

    # Get lengths of the ranges
    l₁ = length(r₁)
    l₂ = length(r₂)

    # Define the heatmap matrix
    m = Matrix{Int}(undef, l₂ - 1, l₁ - 1)
    mₚ = Matrix{Int}(undef, l₂ - 1, l₁ - 1)
    m .= 0

    # For each simulation iteration
    for i = 1:(S+1)
        # Get a copy of matrix on the previous step
        mₚ .= m
        # For each pair of the coordinates
        for j = ((i-1)*N+i):(i*N+i)
            # Check if the orbit is in the limits of the grid
            if (r₁[begin] ≤ x[j] ≤ r₁[end]) && (r₂[begin] ≤ y[j] ≤ r₂[end])
                # Get the indices of the square
                k = 1 + floor(Int, abs(x[j] - r₁[begin]) / h₁)
                l = 1 + floor(Int, abs(y[j] - r₂[begin]) / h₂)

                # Mark the crossed squares, once
                if m[l, k] == mₚ[l, k]
                    m[l, k] += 1
                end
            end
        end
    end

    # Plot the heatmap
    heatmap!(
        p,
        r₁,
        r₂,
        m / (S + 1);
        color = cgrad([RGBA{Float64}(0.0, 0.0, 0.0, 0.0); cgrad(:inferno)[2:end]]),
        xticks = [xticks(p)[1][1][begin] - h₁ * minorticks; xticks(p)[1][1]; xticks(p)[1][1][end] + h₁ * 8],
        yticks = [yticks(p)[1][1][begin] - h₂ * minorticks; yticks(p)[1][1]; yticks(p)[1][1][end] + h₂ * 8]
    )

    # Redraw the orbit
    plot!(
        p,
        x[1:100:(N+1)],
        y[1:100:(N+1)];
        label = "",
        color = palette(:default)[1]
    )

    # Point out the starting position
    scatter!(p, [x[1],], [y[1],]; label = "")

    # Save the figure as PDF and PNG
    savefig(p, joinpath(object_dir, "$(name) (Simulated orbits, $(plane))$(POSTFIX).pdf"))
    savefig(p, joinpath(object_dir, "$(name) (Simulated orbits, $(plane))$(POSTFIX).png"))
end

# For each path in the output directory
for path in readdir(INPUT_DIR; join = true)
    # Check if the path is a directory
    if isdir(path)
        # Define the paths to the binary files
        r_path = joinpath(path, "r.bin")
        z_path = joinpath(path, "z.bin")
        x_path = joinpath(path, "x.bin")
        y_path = joinpath(path, "y.bin")
        apo_path = joinpath(path, "apo.bin")
        peri_path = joinpath(path, "peri.bin")

        # Get the ID of the object
        name = basename(path)

        # Define the object directory
        object_dir = joinpath(OUTPUT_DIR, name)
        mkpath(object_dir)

        # Plot the histogram of apocentric distances if the corresponding data file exists
        if isfile(apo_path) && isfile(peri_path)
            # Prepare the array
            apo = Vector{Float64}(undef, (S + 1))

            # Read the data
            read!(apo_path, apo)

            println(
                " "^4, "> Plotting the histogram of apocentric distances from \"$(name)\"...", '\n',
                " "^4, "  > max: ", maximum(apo), '\n',
                " "^4, "  > min: ", minimum(apo), '\n',
                " "^4, "  > no var.: ", apo[1]
            )

            # Calculate the number of bins
            bins = ceil(Int, (maximum(apo) - minimum(apo)) / (3.49 * std(apo) / S^(1 / 3)) * 1.5)

            # Plot the histogram of apocentric distances
            p = histogram(
                apo;
                bins,
                label = "",
                title = name,
                xlabel = L"r_a \; [\mathrm{kpc}]",
                xminorticks = 10,
                color = "#80cdfd"
            )
            plot!(p, yticks = range(0, ceil(Int, ylims(p)[2]); step = 5))
            vline!(p, [apo[1],]; label = "", lw = 1.5)

            # Save the figure as PDF and PNG
            savefig(p, joinpath(object_dir, "$(name) (Apocentric distances)$(POSTFIX).pdf"))
            savefig(p, joinpath(object_dir, "$(name) (Apocentric distances)$(POSTFIX).png"))

        end

        # Plot the histogram of pericentric distances if the corresponding data file exists
        if isfile(peri_path)
            # Prepare the array
            peri = Vector{Float64}(undef, (S + 1))

            # Read the data
            read!(peri_path, peri)

            # Calculate the number of bins
            bins = ceil(Int, (maximum(peri) - minimum(peri)) / (3.49 * std(peri) / S^(1 / 3)) * 1.5)

            println(
                " "^4, "> Plotting the histogram of pericentric distances from \"$(name)\"...", '\n',
                " "^4, "  > max: ", maximum(peri), '\n',
                " "^4, "  > min: ", minimum(peri), '\n',
                " "^4, "  > no var.: ", peri[1]
            )

            # Plot the histogram of pericentric distances
            p = histogram(
                peri;
                bins,
                label = "",
                title = name,
                xlabel = L"r_p \; [\mathrm{kpc}]",
                xminorticks = 10,
                color = "#80cdfd"
            )
            plot!(p, yticks = range(0, ceil(Int, ylims(p)[2]); step = 5))
            vline!(p, [peri[1],]; label = "", lw = 1.5)

            # Save the figure as PDF and PNG
            savefig(p, joinpath(object_dir, "$(name) (Pericentric distances)$(POSTFIX).pdf"))
            savefig(p, joinpath(object_dir, "$(name) (Pericentric distances)$(POSTFIX).png"))
        end

        # Plot the heatmap of simulated orbits in the RZ plane if the corresponding data files exist
        if isfile(r_path) && isfile(z_path)
            # Prepare the arrays
            r = Vector{Float64}(undef, (N + 1) * (S + 1))
            z = Vector{Float64}(undef, (N + 1) * (S + 1))

            try
                # Read the data
                read!(r_path, r)
                read!(z_path, z)

                println(" "^4, "> Plotting the heatmap of simulated orbits in the RZ plane from \"$(name)\"...")
                plot_heatmap(name, object_dir, r, z, L"R \; [\mathrm{kpc}]", L"Z \; [\mathrm{kpc}]", "RZ")
            catch
            end
        end

        # Plot the heatmap of simulated orbits in the XY plane if the corresponding data files exist
        if isfile(x_path) && isfile(y_path)
            # Prepare the arrays
            x = Vector{Float64}(undef, (N + 1) * (S + 1))
            y = Vector{Float64}(undef, (N + 1) * (S + 1))

            try
                # Read the data
                read!(x_path, x)
                read!(y_path, y)

                println(" "^4, "> Plotting the heatmap of simulated orbits in the XY plane from \"$(name)\"...")
                plot_heatmap(name, object_dir, x, y, L"X \; [\mathrm{kpc}]", L"Y \;\, [\mathrm{kpc}]", "XY")
            catch
            end
        end
    end
end

println()


text:
module Angular
using Compat
using LinearAlgebra
using Statistics
using WignerSymbols: HalfInteger, clebschgordan

# TODO: move this to WignerSymbols
Base.round(::Type{HalfInteger}, x) = HalfInteger(round(Int, 2x), 2)
Base.:*(a::HalfInteger, b::HalfInteger) = (a.numerator * b.numerator) // 4 # converts to Rational

"""
    @halfint(x) -> HalfInteger

Performs a compile-time conversion of a number literal into a `HalfInteger`. This allows
for an easy construction of `HalfInteger` values, from e.g. integers and rationals.

```jldoctest
julia> @halfint(3//2)
3/2

julia> @halfint(3//2) |> typeof
HalfInteger
```

Trying to use `@halfint` on a non-literal value will yield a compile-time error:

```jldoctest
julia> function foo(q)
           @halfint(q)
       end
ERROR: LoadError: ArgumentError: @halfint macro needs a compile-time constant as an argument
```
"""
macro halfint(x)
    value = try
        eval(x)
    catch
        throw(ArgumentError("@halfint macro needs a compile-time constant as an argument"))
    end
    convert(HalfInteger, value)
end

"""
    halfint(x) -> HalfInteger

Converts any valid numeric value to a `HalfInteger`. A shorthand for `convert(HalfInteger, x)`.
"""
halfint(x) = convert(HalfInteger, x)

include("combinatorics.jl")
include("linearalgebra.jl")
include("vectorspaces.jl")
include("compositeoperators.jl")
include("tensorproducts.jl")
include("fockspace.jl")
include("angularmomentum.jl")
include("lscoupling.jl")

end # module


text:
#=
Author:= Shreeda Bhat
License:= MIT
UTF-8 =#

# julia has many logical operator, !, && and ||
#=
! - Logical Not
&& - Logical AND
|| - Logical OR
=#

###### NOT operator########
println("##################Let's start with Not operator#########")

println(!true)# contradicts as True i.e False

println(!false)

########### AND operator##############
println("#######################AND operator###################")

println(true && true)# True

println(true && false)# False

println(false && false)# False

################OR operator############

println("#######################OR operator###################")

println(true || true )# True

println(true || false)# True

println(false || false)# False

#*********************#Comparision operator#**************#
#=
> - Greater Than
< - Lesser Than
== - Equal To
=#
println("#*********************#Comparision operator#**************#")
println(2<56)#True

print(25>23)#True

print(25<2)#False

println(23>25)#False

println(25==25)#True

println(25==50-25)#True

println(25==20)#False

println((10-56>1)&& (10==10)||(10-6>1))#True

#Julia evaluates from left to right for logical expressions
println(true&&false||true) # True

#=

&&  and || are the boolean operators

| and & are the bitwise operators =#


text:
using Test
using SCIP

@testset "direct library calls" begin
    include("direct_library_calls.jl")
end

@testset "managed memory" begin
    include("scip_data.jl")
end

# new type definitions in module (needs top level)
include("conshdlr_support.jl")

@testset "constraint handlers" begin
    include("conshdlr.jl")
end

if SCIP.SCIP_versionnumber() >= v"7"
    include("sepa_support.jl")
    @testset "separators" begin
        include("sepa.jl")
    end
    @testset "cut callbacks" begin
        include("cutcallback.jl")
    end
else
    @info "Separation and callbacks not tested for SCIP versions below 7"
end

const MOI_BASE_EXCLUDED = [
    "Semicontinuous",
    "ScalarAffineFunction_Semiinteger",
    "ScalarQuadraticFunction_Semiinteger",
    "VectorAffineFunction_GeometricMeanCone",
    "Indicator_GreaterThan",
    "Indicator_LessThan",
    "Indicator_ACTIVATE_ON_ZERO", # odd MOI bug?
    "test_constraint_get_ConstraintIndex", # accessing constraint from string name
    "test_model_ListOfConstraintAttributesSet",
    "BoundAlreadySet", # see TODO,
    "ListOfConstraintIndices", # MathOptInterface.ListOfModelAttributesSet
    "ListOfConstraintTypesPresent",
    "ScalarAffineFunction_ConstraintName", # get(::SCIP.Optimizer, ::Type{MathOptInterface.ConstraintIndex}, ::String)
    "UnsupportedAttribute",    # test_model_copy_to_UnsupportedAttribute: MOI.copy_to(model, BadVariableAttributeModel())
    "UnsupportedConstraint",
    "test_model_delete", # MOI.ListOfConstraintTypesPresent
    "duplicate_VariableName", # two identical variable names should error
    "test_modification_delete_variables_in_a_batch",
    "test_modification_set_function_single_variable",
    "test_modification_set_scalaraffine_",
    "test_modification_set_singlevariable_",
    "test_modification_transform_",
    "test_nonlinear_", # None of tests provide expression graphs in the evaluator.
    "FEASIBILITY_SENSE", # TODO
    "ObjectiveFunction_ScalarAffineFunction",
    "test_objective_set_via_modify", # ListOfModelAttributesSet
]

# second order cone
# SCIP needs non-negative variable on right-hand side of second-order cone constraint
append!(MOI_BASE_EXCLUDED,[
    "test_conic_RotatedSecondOrderCone_INFEASIBLE",
    "test_conic_RotatedSecondOrderCone_VectorAffineFunction",
    "test_conic_RotatedSecondOrderCone_VectorOfVariables",
    "test_conic_RotatedSecondOrderCone_out_of_order",
    "test_conic_SecondOrderCone_INFEASIBLE",
    "test_conic_SecondOrderCone_Non",
    "test_conic_SecondOrderCone_VectorAffineFunction",
    "test_conic_SecondOrderCone_VectorOfVariables",
    "test_conic_SecondOrderCone_negative_initial_bound",
    "test_conic_SecondOrderCone_negative_post_bound",
    "test_conic_SecondOrderCone_no_initial_bound",
    "test_conic_SecondOrderCone_out_of_order",
    "test_constraint_PrimalStart_DualStart_SecondOrderCone",
    "test_quadratic_Integer_SecondOrderCone",
])

@testset "MathOptInterface tests (bridged)" begin
    include("MOI_wrapper_bridged.jl")
end

@testset "MathOptInterface tests (bridged & cached)" begin
    include("MOI_wrapper_cached.jl")
end

@testset "MathOptInterface tests (direct)" begin
    include("MOI_wrapper_direct.jl")
end

@testset "MathOptInterface additional tests" begin
    include("MOI_additional.jl")
end

@testset "MathOptInterface nonlinear expressions" begin
    include("MOI_nonlinear_exprs.jl")
end

include("MOI_conshdlr.jl")


text:
# Common graduated values classification
struct EqualIntervals <:  AbstractGraduatedMapClassificator
    k::Int
    EqualIntervals(k::Int = 5) = return new(k);
end

struct Quantiles <: AbstractGraduatedMapClassificator
    k::Int
    Quantiles(k::Int = 5) = return new(k);
end

struct NaturalBreaks <: AbstractGraduatedMapClassificator
    k::Int
    NaturalBreaks(k::Int = 5) = return new(k);
end

struct CustomBreaks <: AbstractGraduatedMapClassificator
    k::Int
    bins::Union{Vector{Int}, Vector{Float64}}
    CustomBreaks(bins::Union{Vector{Int}, Vector{Float64}}) = return new(length(bins) + 1, bins);
end

# Bounds
function mapclassifybounds(mcr::EqualIntervals, x::Vector{T} where T<:Real)::Tuple{Vector{Float64}, Vector{Float64}}
    k = mcr.k
    minx = minimum(x)
    maxx = maximum(x)
    stepinc = (maxx - minx) / k

    lbound = zeros(k)
    ubound = zeros(k)

    for i in 1:k
        lbound[i] = minx + (i - 1) * stepinc
        ubound[i] = minx + i * stepinc
    end

    lbound, ubound
end

function mapclassifybounds(mcr::Quantiles, x::Vector{T} where T<:Real)::Tuple{Vector{Float64}, Vector{Float64}}
    k = mcr.k
    qx = quantile(x, range(1/k, 1, length = k))

    lbound = zeros(k)
    lbound[1] = minimum(x)
    lbound[2:k] = qx[1:k-1]

    ubound = qx

    lbound, ubound
end

function mapclassifybounds(mcr::CustomBreaks, x::Vector{T} where T<:Real)::Tuple{Vector{Float64}, Vector{Float64}}
    bins = mcr.bins

    lbound = vcat(minimum(x), bins)
    ubound = vcat(bins, maximum(x))

    lbound, ubound
end

function mapclassifybounds(mcr::NaturalBreaks, x::Vector{T} where T<:Real)::Tuple{Vector{Float64}, Vector{Float64}}
    k = mcr.k
    n = length(x)
    xs = sort(x)

    (k < n) || throw(ArgumentError("number of classes should be lower than number of observations"))

    # Create matrices and set initial values
    mat1 = zeros(Int64,   n, k) # Lower class limits
    mat2 = zeros(Float64, n, k) # Variance combinations

    mat1[1, :]     .= 1
    mat2[2:end, :] .= Inf   # Set to the most larger possible value

    # Loop through each observation
    v = 0.0
    @inbounds for l in 2:n
        s1 = 0.0 # sum
        s2 = 0.0 # sum of squares
        w = 0 # Number of observations considered

        # Loop through observations till l in reverse order
        for m in reverse(1:l)
            val = xs[m]

            # Increase sum of sqaures and observations
            s1 = s1 + val
            s2 = s2 + val * val
            w = w + 1

            # Variance
            v = s2 - (s1 * s1) / w

            # Update lower bound and variance if newer is lower
            prev = m - 1
            if prev != 0
                for j in 2:k
                    if mat2[l, j] >= (v + mat2[prev, j - 1])
                        mat1[l, j] = m
                        mat2[l, j] = v + mat2[prev, j - 1]
                    end
                end
            end
        end

        # Store lower class limits and variance
        mat1[l, 1] = 1
        mat2[l, 1] = v

    end

    # Get class bounds
    kclass = zeros(k + 1)
    kclass[end] = xs[end]
    kclass[1] = xs[1]

    # Get class bounds
    lbound = zeros(k)
    ubound = zeros(k)

    pivot = n
    for ki in reverse(2:k)
        pivot = mat1[pivot, ki] - 1
        lbound[ki] = xs[pivot]
    end

    lbound[1] = minimum(xs)
    ubound[1:k-1] = lbound[2:k]
    ubound[k] = maximum(xs)

    lbound, ubound
end

# Structure for Graduated  Map Classification
struct GraduatedMapClassification <: AbstractMapClassification
    mcr::AbstractMapClassificator
    k::Int
    group::Vector{Int}
    ngroup::Vector{Int}
    lbound::Vector{Float64}
    ubound::Vector{Float64}
    lower::Symbol
    upper::Symbol
end

function mapclassify(mcr::AbstractGraduatedMapClassificator, x::Vector{T} where T<:Real;
    lower::Symbol = :open, upper::Symbol = :closed)::GraduatedMapClassification

    ((lower == :open) || (lower == :closed )) || throw(ArgumentError("lower bound must be :open or :closed"))
    ((upper == :open) || (upper == :closed )) || throw(ArgumentError("upper bound must be :open or :closed"))
    (lower !== upper) || throw(ArgumentError("lower and upper bounds cannot be both :open or :closed"))

    lbound, ubound = mapclassifybounds(mcr, x)

    k = mcr.k
    group = zeros(Int, length(x))
    ngroup = zeros(Int, k)

    for i in 1:k
        if i == 1
            # First group
            if upper == :open
                subset = (x .>= lbound[i]) .& (x .< ubound[i])
            elseif upper == :closed
                subset = (x .>= lbound[i]) .& (x .<= ubound[i])
            end
        elseif i == k
            # Last group
            if lower == :open
                subset = (x .> lbound[i]) .& (x .<= ubound[i])
            else
                subset = (x .>= lbound[i]) .& (x .<= ubound[i])
            end
        else i < k
            # Intermediate groups
            if lower == :open && upper == :closed
                subset = (x .> lbound[i]) .& (x .<= ubound[i])
            elseif lower == :closed && upper == :open
                subset = (x .>= lbound[i]) .& (x .< ubound[i])
            end
        end

        group[subset] .= i
        ngroup[i] = sum(subset)
    end

    return GraduatedMapClassification(mcr, k, group, ngroup, lbound, ubound, lower, upper)
end

function maplabels(mc::GraduatedMapClassification; digits::Int = 2, sep::String = ", ", counts::Bool = true)::Vector{String}
    k = mc.k
    ngroup = mc.ngroup
    lbound = mc.lbound
    ubound = mc.ubound
    lower = mc.lower
    upper = mc.upper

    labels = Vector{String}(undef, k)

    # Label for each group in the classification
    for i in 1:k
        labeli = ""

        # Lower interval Symbol
        if i == 1
            labeli = labeli * "["
        else
            if lower == :open
                labeli = labeli * "("
            elseif lower == :closed
                labeli = labeli * "["
            end
        end

        # Range
        labeli = labeli * string(round(lbound[i], digits = digits)) *
            sep *
            string(round(ubound[i], digits = digits))

        # Upper interval symbol
        if i == k
            labeli = labeli * "]"
        else
            if upper == :open
                labeli = labeli * ")"
            elseif upper == :closed
                labeli = labeli * "]"
            end
        end

        # Counts
        if counts
            labeli = labeli * " (" * string(ngroup[i]) * ")"
        end

        labels[i] = labeli
    end

    labels
end

bounds(mc::GraduatedMapClassification)::Tuple{Vector{Float64}, Vector{Float64}} = return mc.lbound, mc.ubound;

levels(mc::GraduatedMapClassification)::Vector{String} = return string.(1:mc.k);


text:
#############################################################################
# ECOS.jl
# Wrapper around the ECOS solver https://github.com/ifa-ethz/ecos
# See http://github.com/JuliaOpt/ECOS.jl
#############################################################################
# test/mpb.jl
# Test the MathProgBase.jl interface for the ECOS.jl solver wrapper
#############################################################################

using Base.Test
importall MathProgBase
using ECOS

s = ECOSSolver()

# Problem 1 - all vars in nonneg cone
# min -3x - 2y - 4z
# st    x +  y +  z == 3
#            y +  z == 2
#       x>=0 y>=0 z>=0
# Opt solution = -11
# x = 1, y = 0, z = 2
m = MathProgBase.model(s)
MathProgBase.loadconicproblem!(m,
                    [-3.0, -2.0, -4.0],
                    [ 1.0   1.0   1.0;
                      0.0   1.0   1.0],
                    [ 3.0,  2.0],
                    [(:NonNeg, 1:3)])
MathProgBase.optimize!(m)
@test MathProgBase.status(m) == :Optimal
@test_approx_eq_eps MathProgBase.getobjval(m) -11 1e-6
@test_approx_eq_eps MathProgBase.getsolution(m)[1] 1.0 1e-6
@test_approx_eq_eps MathProgBase.getsolution(m)[2] 0.0 1e-6
@test_approx_eq_eps MathProgBase.getsolution(m)[3] 2.0 1e-6

# Problem 2 - mixed free, nonneg, nonpos, zero, shuffled cones
# min  3x + 2y - 4z + 0s
# st    x           -  s  == -4    (i.e. x >= -4)
#            y            == -3
#       x      +  z       == 12
#       x free
#       y <= 0
#       z >= 0
#       s zero
# Opt solution = -82
# x = -4, y = -3, z = 16, s == 0
m = MathProgBase.model(s)
MathProgBase.loadconicproblem!(m,
                    [ 3.0,  2.0, -4.0,  0.0],
                    [ 1.0   0.0   0.0  -1.0;
                      0.0   1.0   0.0   0.0;
                      1.0   0.0   1.0   0.0],
                    [-4.0, -3.0, 12.0],
                    [(:Free,1), (:NonNeg,3), (:Zero,4), (:NonPos,2)])
MathProgBase.optimize!(m)
@test MathProgBase.status(m) == :Optimal
@test_approx_eq_eps MathProgBase.getobjval(m) -82 1e-6
@test_approx_eq_eps MathProgBase.getsolution(m)[1] -4.0 1e-6
@test_approx_eq_eps MathProgBase.getsolution(m)[2] -3.0 1e-6
@test_approx_eq_eps MathProgBase.getsolution(m)[3] 16.0 1e-6
@test_approx_eq_eps MathProgBase.getsolution(m)[4]  0.0 1e-6

# Problem 3 - SOC
# min 0x - 1y - 1z
#  st  x            == 1
#      x >= ||(y,z)||
m = MathProgBase.model(s)
MathProgBase.loadconicproblem!(m,
                    [ 0.0, -1.0, -1.0],
                    [ 1.0   0.0   0.0],
                    [ 1.0],
                    [(:SOC,1:3)])
MathProgBase.optimize!(m)
@test MathProgBase.status(m) == :Optimal
@test_approx_eq_eps MathProgBase.getobjval(m) -sqrt(2.0) 1e-6
@test_approx_eq_eps MathProgBase.getsolution(m)[1] 1.0 1e-6
@test_approx_eq_eps MathProgBase.getsolution(m)[2] 1.0/sqrt(2.0) 1e-6
@test_approx_eq_eps MathProgBase.getsolution(m)[3] 1.0/sqrt(2.0) 1e-6