text:
#
# Methods for SOM visualisation via PyPlot, MatPlotLib:
#
# Plots are always in a range of (0,0) to (100,100) with the legends
# in the region between (100,0) and (150,100).
#

ρ = 1.0                 # global scaling
fig = nothing           # the figure
█ = nothing             # \blockfull: the subplot into which everything is plotted

"""
    sigmoidScale(num, perc, min)

Return a non-linear scaled colour value. 85%-percentile is mapped to
tanh(1.5). Scaled values 0 < x < min are returned as min.
"""
function sigmoidScale(num, perc, min = 0.1)

    TANH_ANCHOR = 1.5     # perc is mapped to 1.5 on tanh-scale

    # keep precise 0:
    #
    if num == 0
        return 0.0
    end

    # set 85%-percentile to 1.5 for scale colour with tanh
    # modify so that 0->0.0, 1->min, oo->1.0:
    #
    perc = floor(perc) + 1
    colNum = (tanh((num-1) / perc * TANH_ANCHOR) * (1-min)) + min

    # start scale at MIN_VAL
    #
    if colNum < min
        colNum = min
    end

    return colNum
end

function getPaperSizes(paper::Symbol; margin = 1.0)

    if paper == :a4
        size = (8.3, 11.7)
    elseif paper == :a4r
        size = (11.7, 8.3)
    elseif paper == :letter
        size = (8.5, 11.0)
    elseif paper == :letterr
        size = (11.0, 8.5)
    else
        size = (10.0, 10.0)
    end

    size = size .- margin
    return size
end

function adjustLineWidth(som::Som, lwd)

    defaultSize = 10
    maxDim = max(som.xdim, som.ydim)

    return Int(ceil(defaultSize / maxDim * lwd))
end

function drawBackGround(som)

    if som.topol == :hexagonal
        x = (-0.75, som.xdim+0.5)
        y = (-0.75, (som.ydim-1)*0.866+0.75)
        █.set_xlim(x)
        █.set_ylim(y)
    elseif som.topol == :spherical
        minx, miny, minz = minimum(som.grid, dims = 1)
        maxx, maxy, maxz = maximum(som.grid, dims = 1)
        x = (minx - 0.5, maxx+0.5)
        y = (miny - 0.5, maxy+0.5)
        z = (minz - 0.5, maxz+0.5)
        █.set_xlim(x)
        █.set_ylim(y)
        █.set_zlim(z)
    else
        x = (-0.75, som.xdim+0.5)
        y = (-0.75, som.ydim-1+0.75)
        █.set_xlim(x)
        █.set_ylim(y)
    end

    █.axis("off")
    █.set_aspect("equal")
    fig.patch.set_facecolor("white")
end

function drawTitle(title)

    title = PyPlot.title(title)
    title.set_position([0.5, 1.05])
    # family = "sans-serif",
    # size = 15
end

function drawColorBar(som, population, colourMap)

    # set 85%-percentile to 1.5 for scale colour with tanh:
    #
    perc85 = percentile(population, 85)
    maxPop = maximum(population)

    # Generate fake ScalarMappable for colorbar:
    #
    sm = cm.ScalarMappable(cmap=colourMap)
    sm.set_array([])

    # # make new axws for colorbar (for better positioning, punt colourbar
    # # into extra axes):
    # #
    # cbAxes = fig.add_axes([0.9, 0.1, 0.03, 0.8])
    # cbAxes.axis("off")
    # cbar = PyPlot.colorbar(sm, cax = cbAxes)
    # # cbar = PyPlot.colorbar(sm)
    # cbar.set_label("Population of neurons")

    divider = ag1.make_axes_locatable(█)
    cax = divider.append_axes("right", size=0.25, pad=0.05)
    cbar = PyPlot.colorbar(sm, cax=cax)
    #cax.axis("off")

    # tick positions depend on absolute numbers of
    ticks = [1, Int(round(maxPop*0.25)), Int(round(maxPop*0.5)),
             Int(round(maxPop*0.75)), maxPop]
    tickLabels = ticks
    tickPos = sigmoidScale.(ticks, perc85)

    cbar.set_ticks(tickPos)
    cbar.set_ticklabels(ticks)
#    cbar.set_ticklabels([8, 4, 1], update_ticks = true)
    cbar.update_ticks()
end

function drawFreqLegend(colours)

    # create proxy patches for legend:
    #
    classes = sort(collect(keys(colours)))
    proxys = [patches.Patch(color=colours[c], label=string(c)) for c in classes]
    PyPlot.legend(handles = proxys,
                  bbox_to_anchor=(1.05, 0.1), loc=3, borderaxespad=0.,
                  frameon=false)
end

"""
    drawSquare(x, y, siz, lwd, fill::Bool, linCol, fillCol)

Draw a square with diameter size around (x,y).
"""
function drawSquare(x, y, siz, lwd, fill::Bool, linCol, fillCol)

    δ = siz / 2
    points = [x-δ y-δ;
              x-δ y+δ;
              x+δ y+δ;
              x+δ y-δ]
              points = points .* ρ

    poly = patches.Polygon(points, closed = true,
                                   fill = fill, facecolor = fillCol,
                                   linewidth = lwd, edgecolor = linCol)
    █.add_artist(poly)
end

"""
    drawHexagon(x, y, siz, lwd, fill::Bool, linCol, fillCol)

Draw an hexagon with diameter size around (x,y).
"""
function drawHexagon(x, y, siz, lwd, fill::Bool, linCol, fillCol)

    δy = 0.5 * siz * tan(π/6)
    h  = 0.5 * siz / cos(π/6)
    δx = 0.5 * siz
    points = [x-δx y-δy;
              x-δx y+δy;
              x    y+h;
              x+δx y+δy;
              x+δx y-δy;
              x    y-h]
    points = points .* ρ

    poly = patches.Polygon(points, closed = true,
                                   fill = fill, facecolor = fillCol,
                                   linewidth = lwd, edgecolor = linCol)
    █.add_artist(poly)
end

"""
    drawNeurons(som::Som, col, lwd::Number)

Draw all neurons as empty squares or hexagons.
"""
function drawNeurons(som::Som; col = "black", lwd = 1)

    # set shape to topol:
    #
    if som.topol == :hexagonal
        draw = drawHexagon
    else
        draw = drawSquare
    end

    # adjust linewidth:
    #
    lwd = adjustLineWidth(som, lwd)

    for i in 1:som.nCodes
        draw(som.grid[i,1], som.grid[i,2], 1.0, lwd, false, col, "none")
    end
end

"""
    drawDesity(som::Som, population, colormap)

Draw the population of neurons as colors.
"""
function drawDensity(som::Som, population, colourMap)

    # set shape to topol:
    #
    if som.topol == :hexagonal
        draw = drawHexagon
    else
        draw = drawSquare
    end

    # set 85%-percentile to 1.5 for scale colour with tanh:
    #
    perc85 = percentile(population, 85)

    cmap = get_cmap(colourMap)
    for i in 1:som.nCodes
        colNum = sigmoidScale( population[i], perc85)
        col = colNum > 0 ? cmap(colNum) : "white"
        # println("$i: pop: $(som.population[i]), $colNum, col: $col")
        draw(som.grid[i,1], som.grid[i,2], 1.0, 0, true, "none", col)
    end
end

function drawPieChart(x, y, freqs, colours, lwd, d)

    r = d / 2
    classes = names(freqs)
    pieStart = 0.0
    pieEnd = 0.0

    for c in sort(classes)

        pieStart = pieEnd
        pieEnd = pieStart + freqs[c] * 360
        wedge = patches.Wedge((x*ρ,y*ρ), r*ρ, pieStart, pieEnd,
                              fill = true, facecolor = colours[c],
                              linewidth = 0, edgecolor = "none")
        █.add_artist(wedge)
    end

    circle = patches.Circle((x*ρ,y*ρ), radius = r*ρ,
                            fill = false,
                            linewidth = lwd, edgecolor = "0.5")
    █.add_artist(circle)

end

function drawPieCharts(som::Som, frequencies::DataFrame, colours)

    # set 85%-percentile to 1.5 for scale pieChart with tanh:
    #
    perc85 = percentile(frequencies.Population, 85)

    # loop neurons:
    #
    for i in 1:som.nCodes

        if frequencies[i, :Population] != 0

            r = sigmoidScale(frequencies[i, :Population], perc85)
            drawPieChart( som.grid[i,1], som.grid[i,2],
                          frequencies[i,5:end],
                          colours, adjustLineWidth(som, 1), r)
        end
    end
end

"""
    printToDevice(device, fName)

Print to file or display.

Prints the current plot to the device `device`. If fName ends with
`.png, .svg, .pdf` this overrides the setting in device.
"""
function printToDevice(device, fName)

    # select device from filename, rather than from device:
    #
    if occursin(r"\.[a-zA-Z]{3}$", fName)
        println("Writing plot to file $fName")
        PyPlot.savefig(fName)
    elseif device != :display
        ext = string(device)
        fName = fName * "." * ext
        println("Writing plot to file $fName")
        PyPlot.savefig(fName)
    else
        PyPlot.show()
    end
end

function drawPopulation(som::Som, population, title, paper, colourMap, device, fName)

    paperSize = getPaperSizes(paper)
    global fig = PyPlot.figure(figsize=paperSize)
    global █ = fig.add_subplot(1,1,1)

    println("DrawPopulation()")
    drawBackGround(som)
    drawTitle(title)

    drawDensity(som, population, colourMap)
    drawNeurons(som, col = "0.5", lwd = 1)
    drawColorBar(som, population, colourMap)

    printToDevice(device, fName)
end

function drawFrequencies(som::Som, frequencies, title, paper,
                         colours, device, fName)

    paperSize = getPaperSizes(paper)
    global fig = PyPlot.figure(figsize=paperSize)
    global █ = fig.add_subplot(1,1,1)
    drawBackGround(som)
    drawTitle(title)

    drawPieCharts(som, frequencies, colours)
    drawNeurons(som, col = "0.5", lwd = 1)
    drawFreqLegend(colours)

    printToDevice(device, fName)

end

function drawSpherePopulation(som::Som, population, detail, title, paper,
                              colourMap, device, fName)

    paperSize = getPaperSizes(paper)
    global fig = PyPlot.figure(figsize=paperSize)
    global █ = fig.add_subplot(1,1,1, projection="3d")
    drawBackGround(som)
    # █.axis("off")
    # █.set_aspect("equal")
    # fig.patch.set_facecolor("white")
    drawTitle(title)

    drawSphereColour(som, population, detail,colourMap)
    #drawColorBar(som, population, colourMap)

    printToDevice(device, fName)
end

function drawSphereFreqs(som::Som, frequencies, detail, title, paper,
                              colours, device, fName)

    paperSize = getPaperSizes(paper)
    global fig = PyPlot.figure(figsize=paperSize)
    global █ = fig.add_subplot(1,1,1, projection="3d")
    █.axis("off")
    █.set_aspect("equal")
    fig.patch.set_facecolor("white")
    drawTitle(title)

    drawSphereClasses(som, frequencies, detail, colours)
    drawFreqLegend(colours)

    printToDevice(device, fName)
end


text:
"""
License: MIT
Julia 1.02
Dynamic programming formulation to count number of matchings of maximum cardinality
"""

# function from the paper
# INPUT:
# G and E - can be provided by read_Graph(filename)
# maxCard - can be prodived by OptimalMatching(G,E)
function Count_Mat_Paper(G, maxCard,E,r,s,aux)
    # reset table of results
    if aux == 1
        global F_table = Dict()
    end
    # do not repeat computations:
    try
        # if aux == 1 then that entry wont exist
        return F_table[r,maxCard,[if s[v]>0 1 else 0 end for v in sort(collect(keys(G)))]]
    catch
        # if any of the entries maxCard, r or s is negative return 0
        if maxCard<0 || r<0 || minimum(values(s))<0
            F_table[r,maxCard,[if s[v]>0 1 else 0 end for v in sort(collect(keys(G)))]]=0
            # aux =1 means we are in the end of the recurssion so return table
            if aux ==1
                return F_table, 0
            end
            # if not in end of the recurssion return 0 (no need to return table since it is a global variable)
            return 0
        end
        # matchings containing more edges than the ones available
        if maxCard>r
            F_table[r,maxCard,[if s[v]>0 1 else 0 end for v in sort(collect(keys(G)))]]=0
            # aux =1 means we are in the end of the recurssion so return table
            if aux ==1
                return F_table, 0
            end
            # if not in end of the recurssion return 0 (no need to return table since it is global variable)
            return 0
        end
        # if we want a matching of cardinality 0, the only choice is M= empty set, i.e., there only one solution
        if maxCard ==0 && r >=0 && minimum(values(s))>=0
            F_table[r,maxCard,[if s[v]>0 1 else 0 end for v in sort(collect(keys(G)))]]=1
            # aux =1 means we are in the end of the recurssion so return table
            if aux ==1
                return F_table, 1
            end
            # if not in end of the recurssion return 0 (no need to return table since it is global variable)
            return 1
        # if there are no edges available, i.e., r=0, and maxCard>0 then there is no matching under these condition
        elseif maxCard >0 && r ==0 && minimum(values(s))>=0
            F_table[r,maxCard,[if s[v]>0 1 else 0 end for v in sort(collect(keys(G)))]]=0
            # aux =1 means we are in the end of the recurssion so return table
            if aux ==1
                return F_table, 0
            end
            # if not in end of the recurssion return 0 (no need to return table since it is global variable)
            return 0
        # OPT >0 and r>0 and min(s)>=0
        else
            No_r = Count_Mat_Paper(G,maxCard,E,r-1,s,0)
            Yes_r = Count_Mat_Paper(G,maxCard-1,E,r-1,Dict(if v in E[r] v=>s[v]-1 else v=>s[v] end for v in sort(collect(keys(G)))),0)
            F_table[r,maxCard,[if s[v]>0 1 else 0 end for v in sort(collect(keys(G)))]] = No_r+Yes_r
            if aux ==1
                return F_table,No_r+Yes_r
            end
            return No_r+Yes_r
        end
    end
end

# read graphs

function Read_Graph(filename)
    G_tpm = Dict()
    f = open(filename)
    lines =  readlines(f)
    close(f)
    num_V, num_E = split(lines[1])
    num_V = parse(Int,num_V)
    num_E = parse(Int,num_E)
    for i=1:num_E
        v1,v2,w = split(lines[1+i])
        v1 = parse(Int,v1)
        v2 = parse(Int,v2)
        if !(v2 in keys(G_tpm))
            G_tpm[v2] = []
        end
        if !(v1 in keys(G_tpm))
            G_tpm[v1] = [v2]
        else
            append!(G_tpm[v1],v2)
        end
    end
    G = Dict(v=>[] for v in keys(G_tpm))
    E = []
    for v1 in keys(G_tpm)
        for v2 in G_tpm[v1]
            if v1 in G_tpm[v2] && !(v1 in G[v2]) && !(v2 in G[v1])
                append!(G[v1],v2)
                append!(G[v2],v1)
                append!(E,[Set([v1,v2])])
            end
        end
    end
    if length(keys(G))==0
        G = Dict(v=>[] for v in 1:num_V)
    end
    return G,E
end

using CPLEX, JuMP
function OptimalMatching(G,E)
    #m = Model(solver=CplexSolver())
    m = Model(with_optimizer(CPLEX.Optimizer))
    # only one core:
    CPX_PARAM_THREADS = 1
    CPX_PARAM_MIPDISPLAY =0
    CPX_PARAM_TILIM = 3600
    num_E = length(E)
    if num_E==0
        return 0
    end
    @variable(m, x[1:num_E], Bin)
    for v in keys(G)
        @constraint(m,sum(x[j[1]] for j in enumerate(E) if v in j[2])<=1)
    end
    @objective(m, Max, sum(x[j[1]] for j in enumerate(E)))
    optimize!(m)
    println("Objective value: ", objective_value(m))
    #println("Objective value: ", getobjectivevalue(m))
    return objective_value(m)
end

ext: jl

max_stars_repo_name: mxmmargarida/KEG

max_issues_repo_name: mxmmargarida/KEG

max_forks_repo_name: mxmmargarida/KEG


text:
module Planck

using Unitful
using Unitful: AbstractQuantity
using Unitful: h, c0, k

export blackbody

# constants in SI units
const _h = ustrip(u"J*s", h)
const _c0 = ustrip(u"m/s", c0)
const _k = ustrip(u"J/K", k)

using Unitful: 𝐋, 𝐓

"""
    blackbody([OT], x, T)

Evaluate a blackbody[^1] at `x`, which is in meters by default. If `OT` is given, the output will be converted to that type, which is convenient for unit conversions. Temperature is assumed to be Kelvin.

!!! warning "Units"

    If you do not use [Unitful.jl](https://github.com/PainterQubits/Unitful.jl), the defaults will be [SI](https://en.wikipedia.org/wiki/International_System_of_Units)

# Examples

```jldoctest
julia> using Unitful

julia> blackbody(6e-7, 5850)
2.583616647617973e13

julia> blackbody(Float32, 6e-7, 5850)
2.5836166f13

julia> blackbody(u"erg/s/cm^2/nm/sr", 600u"nm", 5850u"K")
2.5836166476179734e7 erg nm⁻¹ cm⁻² s⁻¹ sr⁻¹
```

# References
[^1]: [Blackbody radiation](https://en.wikipedia.org/wiki/Planck%27s_law)
"""
blackbody(OT, x, T) = OT(blackbody(x, T))

function blackbody(ν::AbstractQuantity{V,inv(𝐓)}, T) where V
    2 * h * ν^3 / c0^2 / expm1(h * ν / (k * T))
end

function blackbody(λ::AbstractQuantity{V,𝐋}, T) where V
    2 * h * c0^2 / λ^5 / expm1(h * c0 / (λ * k * T))
end

# SI (meters, Kelvin)
function blackbody(λ, T)
    2 * _h * _c0^2 / λ^5 / expm1(_h * _c0 / (λ * _k * T)) # W / m^3
end

end


text:
###############################################################################
#
#   SparsePoly.jl : Generic sparse univariate polynomials over rings
#
###############################################################################

export SparsePolynomialRing

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

parent_type(::Type{SparsePoly{T}}) where {T <: RingElement} = SparsePolyRing{T}

elem_type(::Type{SparsePolyRing{T}}) where {T <: RingElement} = SparsePoly{T}

function isdomain_type(a::Type{SparsePoly{T}}) where T <: RingElement
   return isdomain_type(T)
end

function isexact_type(a::Type{SparsePoly{T}}) where T <: RingElement
   return isexact_type(T)
end

var(a::SparsePolyRing) = a.S

function gen(a::SparsePolyRing)
   return a([one(base_ring(a))], [UInt(1)])
end

###############################################################################
#
#   Basic manipulation
#
###############################################################################

function coeff(x::SparsePoly, i::Int)
   i < 0 && throw(DomainError())
   return x.coeffs[i + 1]
end

one(R::SparsePolyRing) = R(1)

zero(R::SparsePolyRing) = R(0)

iszero(x::SparsePoly) = length(x) = 0

isone(x::SparsePoly) = x == 1

length(x::SparsePoly) = x.length

lead(x::SparsePoly) = x.length == 0 ? base_ring(x)() : x.coeffs[x.length]

trail(x::SparsePoly) = x.length == 0 ? base_ring(x)() : x.coeffs[1]

function normalise(a::SparsePoly, n::Int)
   while n > 0 && iszero(a.coeffs[n])
      n -= 1
   end
   return n
end

base_ring(a::SparsePoly) = base_ring(parent(a))

base_ring(R::SparsePolyRing{T}) where {T <: RingElement} = R.base_ring::parent_type(T)

parent(a::SparsePoly) = a.parent

function Base.deepcopy_internal(a::SparsePoly{T}, dict::ObjectIdDict) where {T <: RingElement}
   Re = Base.deepcopy_internal(a.exps, dict)
   Rc = Array{T}(a.length)
   for i = 1:a.length
      Rc[i] = Base.deepcopy_internal(a.coeffs[i], dict)
   end
   return parent(a)(Rc, Re)
end

###############################################################################
#
#   String I/O
#
###############################################################################

function show(io::IO, x::SparsePoly)
    len = length(x)
    U = string(var(parent(x)))
    if len == 0
      print(io, base_ring(x)(0))
    else
      for i = 1:len
        c = coeff(x, len - i)
        bracket = needs_parentheses(c)
        if i != 1 && !isnegative(c)
           print(io, "+")
        end
        X = x.exps[len - i + 1]
        if !isone(c) && (c != -1 || show_minus_one(typeof(c)))
          if bracket
            print(io, "(")
          end
          show(io, c)
          if bracket
            print(io, ")")
          end
          if c != 1 && !(c == -1 && !show_minus_one(typeof(c))) && X != UInt(0)
             print(io, "*")
          end
        end
        if c == -1 && !show_minus_one(typeof(c))
          print(io, "-")
        end
        if X == UInt(0)
          if c == 1
             print(io, c)
          elseif c == -1 && !show_minus_one(typeof(c))
             print(io, 1)
          end
        end
        n = reinterpret(Int, X)
        if n != 0
           print(io, U)
        end
        if n > 1
           print(io, "^", n)
        end
      end
    end
end

function show(io::IO, p::SparsePolyRing)
   print(io, "Sparse Univariate Polynomial Ring in ")
   print(io, string(p.S))
   print(io, " over ")
   show(io, base_ring(p))
end

show_minus_one(::Type{SparsePoly{T}}) where {T <: RingElement} = show_minus_one(T)

needs_parentheses(a::SparsePoly{T}) where {T <: RingElement} = length(a) > 1

isnegative(x::SparsePoly) = length(x) <= 1 && isnegative(coeff(x, 0))

###############################################################################
#
#   Arithmetic functions
#
###############################################################################

function -(a::SparsePoly{T}) where {T <: RingElement}
   r = parent(a)()
   fit!(r, length(a))
   for i = 1:length(a)
      r.coeffs[i] = -a.coeffs[i]
      r.exps[i] = a.exps[i]
   end
   r.length = a.length
   return r
end

function +(a::SparsePoly{T}, b::SparsePoly{T}) where {T <: RingElement}
   par = parent(a)
   r = par()
   fit!(r, length(a) + length(b))
   i = 1
   j = 1
   k = 1
   while i <= length(a) && j <= length(b)
      if a.exps[i] < b.exps[j]
         r.coeffs[k] = a.coeffs[i]
         r.exps[k] = a.exps[i]
         i += 1
      elseif a.exps[i] == b.exps[j]
         c = a.coeffs[i] + b.coeffs[j]
         if c != 0
            r.coeffs[k] = c
            r.exps[k] = a.exps[i]
         else
            k -= 1
         end
         i += 1
         j += 1
      else
         r.coeffs[k] = b.coeffs[j]
         r.exps[k] = b.exps[j]
         j += 1
      end
      k += 1
   end
   while i <= length(a)
      r.coeffs[k] = a.coeffs[i]
      r.exps[k] = a.exps[i]
      i += 1
      k += 1
   end
   while j <= length(b)
      r.coeffs[k] = b.coeffs[j]
      r.exps[k] = b.exps[j]
      j += 1
      k += 1
   end
   r.length = k - 1
   return r
end

function -(a::SparsePoly{T}, b::SparsePoly{T}) where {T <: RingElement}
   par = parent(a)
   r = par()
   fit!(r, length(a) + length(b))
   i = 1
   j = 1
   k = 1
   while i <= length(a) && j <= length(b)
      if a.exps[i] < b.exps[j]
         r.coeffs[k] = a.coeffs[i]
         r.exps[k] = a.exps[i]
         i += 1
      elseif a.exps[i] == b.exps[j]
         c = a.coeffs[i] - b.coeffs[j]
         if c != 0
            r.coeffs[k] = c
            r.exps[k] = a.exps[i]
         else
            k -= 1
         end
         i += 1
         j += 1
      else
         r.coeffs[k] = -b.coeffs[j]
         r.exps[k] = b.exps[j]
         j += 1
      end
      k += 1
   end
   while i <= length(a)
      r.coeffs[k] = a.coeffs[i]
      r.exps[k] = a.exps[i]
      i += 1
      k += 1
   end
   while j <= length(b)
      r.coeffs[k] = -b.coeffs[j]
      r.exps[k] = b.exps[j]
      j += 1
      k += 1
   end
   r.length = k - 1
   return r
end

struct heap_sr
   exp::UInt
   n::Int
end

function isless(a::heap_sr, b::heap_sr)
   return a.exp < b.exp
end

function ==(a::heap_sr, b::heap_sr)
   return a.exp == b.exp
end

function heapinsert!(xs::Array{heap_sr, 1}, ys::Array{heap_t, 1}, m::Int, exp::UInt)
   i = n = length(xs) + 1
   @inbounds if i != 1 && exp == xs[1].exp
      ys[m] = heap_t(ys[m].i, ys[m].j, xs[1].n)
      xs[1] = heap_sr(exp, m)
      return
   end
   @inbounds while (j = heapparent(i)) >= 1
      if exp == xs[j].exp
         ys[m] = heap_t(ys[m].i, ys[m].j, xs[j].n)
         xs[j] = heap_sr(exp, m)
         return
      elseif exp < xs[j].exp
         i = j
      else
         break
      end
   end
   push!(xs, heap_sr(exp, 0))
   @inbounds while n > i
      xs[n] = xs[heapparent(n)]
      n >>= 1
   end
   xs[i] = heap_sr(exp, m)
   return
end

function heappop!(xs::Array{heap_sr, 1})
   s = length(xs)
   x = xs[1]
   i = 1
   j = 2
   @inbounds while j < s
      if xs[j].exp >= xs[j + 1].exp
         j += 1
      end
      xs[i] = xs[j]
      i = j
      j *= 2
   end
   exp = xs[s].exp
   j = i >> 1
   @inbounds while i > 1 && exp < xs[j].exp
      xs[i] = xs[j]
      i = j
      j >>= 1
   end
   xs[i] = xs[s]
   pop!(xs)
   return
end

function mul_johnson(a::SparsePoly{T}, b::SparsePoly{T}) where {T <: RingElement}
   par = parent(a)
   R = base_ring(par)
   m = length(a)
   n = length(b)
   if m == 0 || n == 0
      return par()
   end
   H = Array{heap_sr}(0)
   I = Array{heap_t}(0)
   # set up heap
   push!(H, heap_sr(a.exps[1] + b.exps[1], 1))
   push!(I, heap_t(1, 1, 0))
   r_alloc = max(m, n) + n
   Rc = Array{T}(r_alloc)
   Re = Array{UInt}(r_alloc)
   k = 0
   c = R()
   Q = Array{Int}(0)
   while !isempty(H)
      exp = H[1].exp
      k += 1
      if k > r_alloc
         r_alloc *= 2
         resize!(Rc, r_alloc)
         resize!(Re, r_alloc)
      end
      first = true
      @inbounds while !isempty(H) && H[1].exp == exp
         x = H[1]
         heappop!(H)
         v = I[x.n]
         if first
            Rc[k] = a.coeffs[v.i]*b.coeffs[v.j]
            Re[k] = exp
            first = false
         else
            Rc[k] = addmul!(Rc[k], a.coeffs[v.i], b.coeffs[v.j], c)
         end
         if v.j < n || v.j == 1
            push!(Q, x.n)
         end
         while (xn = v.next) != 0
            v = I[xn]
            Rc[k] = addmul!(Rc[k], a.coeffs[v.i], b.coeffs[v.j], c)
            if v.j < n || v.j == 1
               push!(Q, xn)
            end
         end
      end
      @inbounds while !isempty(Q)
         xn = pop!(Q)
         v = I[xn]
         if v.j == 1 && v.i < m
            push!(I, heap_t(v.i + 1, 1, 0))
            Collections.heappush!(H, heap_sr(a.exps[v.i + 1] + b.exps[1], length(I)))
         end
         if v.j < n
            I[xn] = heap_t(v.i, v.j + 1, 0)
            heapinsert!(H, I, xn, a.exps[v.i] + b.exps[v.j + 1]) # either chain or insert v into heap
         end
      end
      if Rc[k] == 0
         k -= 1
      end
   end
   resize!(Rc, k)
   resize!(Re, k)
   return parent(a)(Rc, Re)
end

function *(a::SparsePoly{T}, b::SparsePoly{T}) where {T <: RingElement}
   m = length(a)
   n = length(b)
   if m == 0 || n == 0
      return parent(a)()
   end
   if m < n
      return b*a
   end
   r = a*b.coeffs[1]
   r = shift_left!(r, b.exps[1])
   for i = 2:n
      s = a*b.coeffs[i]
      s = shift_left!(s, b.exps[i])
      r += s
   end
   return r
end

###############################################################################
#
#   Euclidean division
#
###############################################################################

function divrem(a::SparsePoly{T}, b::SparsePoly{T}) where {T <: RingElement}
   par = parent(a)
   R = base_ring(par)
   m = length(a)
   n = length(b)
   maxn = a.exps[m]
   n == 0 && error("Division by zero in divrem_monagan_pearce")
   if a.exps[m] < b.exps[n]
      return par(), a
   end
   H = Array{heap_sr}(0)
   I = Array{heap_t}(0)
   # set up heap
   push!(H, heap_sr(maxn - a.exps[m], 1))
   push!(I, heap_t(0, 1, 0))
   q_alloc = max(m - n, n)
   r_alloc = n
   Qc = Array{T}(q_alloc)
   Qe = Array{UInt}(q_alloc)
   Rc = Array{T}(r_alloc)
   Re = Array{UInt}(r_alloc)
   k = 0
   l = 0
   s = n
   c = R()
   qc = R()
   m1 = -R(1)
   mb = -b.coeffs[n]
   Q = Array{Int}(0)
   reuse = Array{Int}(0)
   while !isempty(H)
      exp = H[1].exp
      k += 1
      if k > q_alloc
         q_alloc *= 2
         resize!(Qc, q_alloc)
         resize!(Qe, q_alloc)
      end
      @inbounds while !isempty(H) && H[1].exp == exp
         x = H[1]
         heappop!(H)
         v = I[x.n]
         if v.i == 0
            qc = addmul!(qc, a.coeffs[m + 1 - v.j], m1, c)
         else
            qc = addmul!(qc, b.coeffs[n + 1 - v.i], Qc[v.j], c)
         end
         if v.i != 0 || v.j < m
            push!(Q, x.n)
         else
            push!(reuse, x.n)
         end
         while (xn = v.next) != 0
            v = I[xn]
            if v.i == 0
               qc = addmul!(qc, a.coeffs[m + 1 - v.j], m1, c)
            else
               qc = addmul!(qc, b.coeffs[n + 1 - v.i], Qc[v.j], c)
            end
            if v.i != 0 || v.j < m
               push!(Q, xn)
            else
               push!(reuse, xn)
            end
         end
      end
      @inbounds while !isempty(Q)
         xn = pop!(Q)
         v = I[xn]
         if v.i == 0
            I[xn] = heap_t(0, v.j + 1, 0)
            heapinsert!(H, I, xn, maxn - a.exps[m - v.j]) # either chain or insert into heap
         elseif v.j < k - 1
            I[xn] = heap_t(v.i, v.j + 1, 0)
            heapinsert!(H, I, xn, maxn - b.exps[n + 1 - v.i] - Qe[v.j + 1]) # either chain or insert into heap
         elseif v.j == k - 1
            s += 1
            push!(reuse, xn)
         end
      end
      if qc == 0
         k -= 1
      else
         d1 = maxn - b.exps[n] >= exp
         texp = maxn - b.exps[n] - exp
         if !d1
            l += 1
            if l >= r_alloc
               r_alloc *= 2
               resize!(Rc, r_alloc)
               resize!(Re, r_alloc)
            end
            Rc[l] = -qc
            Re[l] = maxn - exp
            k -= 1
         else
            tq, tr = divrem(qc, mb)
            if tr != 0
               l += 1
               if l >= r_alloc
                  r_alloc *= 2
                  resize!(Rc, r_alloc)
                  resize!(Re, r_alloc)
               end
               Rc[l] = -tr
               Re[l] = maxn - exp
            end
            if tq != 0
               Qc[k] = tq
               Qe[k] = texp
               for i = 2:s
                  if !isempty(reuse)
                     xn = pop!(reuse)
                     I[xn] = heap_t(i, k, 0)
                     heapinsert!(H, I, xn, maxn - b.exps[n + 1 - i] - Qe[k]) # either chain or insert into heap
                  else
                     push!(I, heap_t(i, k, 0))
                     Collections.heappush!(H, heap_sr(maxn - b.exps[n + 1 - i] - Qe[k], length(I)))
                  end
               end
               s = 1
            else
               k -= 1
            end
         end
      end
      qc = zero!(qc)
   end
   resize!(Qc, k)
   resize!(Qe, k)
   resize!(Rc, l)
   resize!(Re, l)
   reverse!(Qc)
   reverse!(Qe)
   reverse!(Rc)
   reverse!(Re)
   return parent(a)(Qc, Qe), parent(a)(Rc, Re)
end

###############################################################################
#
#   Ad hoc arithmetic functions
#
###############################################################################

function *(a::SparsePoly, n::Union{Integer, Rational, AbstractFloat})
   r = parent(a)()
   fit!(r, length(a))
   j = 1
   for i = 1:length(a)
      c = a.coeffs[i]*n
      if c != 0
         r.coeffs[j] = c
         r.exps[j] = a.exps[i]
         j += 1
      end
   end
   r.length = j - 1
   return r
end

function *(a::SparsePoly{T}, n::T) where {T <: RingElem}
   r = parent(a)()
   fit!(r, length(a))
   j = 1
   for i = 1:length(a)
      c = a.coeffs[i]*n
      if c != 0
         r.coeffs[j] = c
         r.exps[j] = a.exps[i]
         j += 1
      end
   end
   r.length = j - 1
   return r
end

*(n::T, a::SparsePoly{T}) where {T <: RingElem} = a*n

*(n::Union{Integer, Rational, AbstractFloat}, a::SparsePoly) = a*n

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(a::SparsePoly{T}, b::SparsePoly{T}) where {T <: RingElement}
   if length(a) != length(b)
      return false
   end
   for i = 1:length(a)
      if a.exps[i] != b.exps[i] || a.coeffs[i] != b.coeffs[i]
         return false
      end
   end
   return true
end

###############################################################################
#
#   Ad hoc comparison
#
###############################################################################

function ==(a::SparsePoly, b::Union{Integer, Rational, AbstractFloat})
   return length(a) == 0 ? b == 0 : a.length == 1 &&
          a.exps[1] == 0 && a.coeffs[1] == b
end

==(a::Union{Integer, Rational, AbstractFloat}, b::SparsePoly) = b == a

function ==(a::SparsePoly{T}, b::T) where T <: RingElem
   return length(a) == 0 ? b == 0 : a.length == 1 &
          a.exps[1] == 0 && a.coeffs[1] == b
end

==(a::T, b::SparsePoly{T}) where T <: RingElem = b == a

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

function from_exp(R::Ring, a::UInt)
   return R(reinterpret(Int, a))
end

# Implement fps algorithm from "Sparse polynomial powering with heaps" by
# Monagan and Pearce, except that we modify the algorithm to return terms
# in ascending order and we fix some issues in the original algorithm
# http://www.cecm.sfu.ca/CAG/papers/SparsePowering.pdf

function pow_fps(f::SparsePoly{T}, k::Int) where {T <: RingElement}
   par = parent(f)
   R = base_ring(par)
   m = length(f)
   H = Array{heap_sr}(0) # heap
   I = Array{heap_t}(0) # auxilliary data for heap nodes
   # set up output poly coeffs and exponents (corresponds to h in paper)
   r_alloc = k*(m - 1) + 1
   Rc = Array{T}(r_alloc)
   Re = Array{UInt}(r_alloc)
   rnext = 1
   # set up g coeffs and exponents (corresponds to g in paper)
   g_alloc = k*(m - 1) + 1
   gc = Array{T}(g_alloc)
   ge = Array{UInt}(g_alloc)
   gnext = 1
   # set up heap
   gc[1] = f.coeffs[1]^(k-1)
   ge[1] = f.exps[1]*(k - 1)
   Rc[1] = f.coeffs[1]*gc[1]
   Re[1] = f.exps[1]*k
   push!(H, heap_sr(f.exps[2] + ge[1], 1))
   push!(I, heap_t(2, 1, 0))
   Q = Array{Int}(0) # corresponds to Q in paper
   topbit = -1 << (sizeof(Int)*8 - 1)
   mask = ~topbit
   largest = fill(topbit, m) # largest j s.t. (i, j) has been in heap
   largest[2] = 1
   # precompute some values
   fik = Array{T}(m)
   for i = 1:m
      fik[i] = from_exp(R, f.exps[i])*(k - 1)
   end
   kp1f1 = k*from_exp(R, f.exps[1])
   gi = Array{T}(1)
   gi[1] = -from_exp(R, ge[1])
   finalexp = f.exps[m]*(k - 1) + f.exps[1]
   # temporary space
   t1 = R()
   C = R() # corresponds to C in paper
   SS = R() # corresponds to S in paper
   temp = R() # temporary space for addmul
   temp2 = R() # temporary space for add
   # begin algorithm
   while !isempty(H)
      exp = H[1].exp
      gnext += 1
      rnext += 1
      if gnext > g_alloc
         g_alloc *= 2
         resize!(gc, g_alloc)
         resize!(ge, g_alloc)
      end
      if rnext > r_alloc
         r_alloc *= 2
         resize!(Rc, r_alloc)
         resize!(Re, r_alloc)
      end
      first = true
      C = zero!(C)
      SS = zero!(SS)
      while !isempty(H) && H[1].exp == exp
         x = H[1]
         heappop!(H)
         v = I[x.n]
         largest[v.i] |= topbit
         t1 = mul!(t1, f.coeffs[v.i], gc[v.j])
         SS = addeq!(SS, t1)
         if exp <= finalexp
            temp2 = add!(temp2, fik[v.i], gi[v.j])
            C = addmul!(C, temp2, t1, temp)
         end
         if first
            ge[gnext] = exp - f.exps[1]
            first = false
         end
         push!(Q, x.n)
         while (xn = v.next) != 0
            v = I[xn]
            largest[v.i] |= topbit
            t1 = mul!(t1, f.coeffs[v.i], gc[v.j])
            SS = addeq!(SS, t1)
            if exp <= finalexp
               temp2 = add!(temp2, fik[v.i], gi[v.j])
               C = addmul!(C, temp2, t1, temp)
            end
            push!(Q, xn)
         end
      end
      reuse = 0
      while !isempty(Q)
         xn = pop!(Q)
         v = I[xn]
         if v.i < m && largest[v.i + 1] == ((v.j - 1) | topbit)
            I[xn] = heap_t(v.i + 1, v.j, 0)
            heapinsert!(H, I, xn, f.exps[v.i + 1] + ge[v.j]) # either chain or insert v into heap
            largest[v.i + 1] = v.j
         else
            reuse = xn
         end
         if v.j < gnext - 1 && (largest[v.i] & mask) <  v.j + 1
            if reuse != 0
               I[reuse] = heap_t(v.i, v.j + 1, 0)
               heapinsert!(H, I, reuse, f.exps[v.i] + ge[v.j + 1]) # either chain or insert v into heap
               reuse = 0
            else
               push!(I, heap_t(v.i, v.j + 1, 0))
               Collections.heappush!(H, heap_sr(f.exps[v.i] + ge[v.j + 1], length(I)))
            end
            largest[v.i] = v.j + 1
         end
      end
      if C != 0
         temp = divexact(C, from_exp(R, exp) - kp1f1)
         SS = addeq!(SS, temp)
         gc[gnext] = divexact(temp, f.coeffs[1])
         push!(gi, -from_exp(R, ge[gnext]))
         if (largest[2] & topbit) != 0
            push!(I, heap_t(2, gnext, 0))
            Collections.heappush!(H, heap_sr(f.exps[2] + ge[gnext], length(I)))
            largest[2] = gnext
         end
      end
      if SS != 0
         Rc[rnext] = SS
         Re[rnext] = ge[gnext] + f.exps[1]
         SS = R()
      else
         rnext -= 1
      end
      if C == 0
         gnext -= 1
      end
   end
   resize!(Rc, rnext)
   resize!(Re, rnext)
   return parent(f)(Rc, Re)
end

function ^(a::SparsePoly{T}, b::Int) where {T <: RingElement}
   b < 0 && throw(DomainError())
   # special case powers of x for constructing polynomials efficiently
   if length(a) == 0
      return parent(a)()
   elseif length(a) == 1
      return parent(a)([coeff(a, 0)^b], [a.exps[1]*b])
   elseif b == 0
      return parent(a)(1)
   elseif b == 1
      return a
   elseif b == 2
      return a*a
   else
      z = a*a
      for i = 3:b
         z *= a
      end
      return z
   end
end

###############################################################################
#
#   Exact division
#
###############################################################################

function divides_monagan_pearce(a::SparsePoly{T}, b::SparsePoly{T}) where {T <: RingElement}
   par = parent(a)
   R = base_ring(par)
   m = length(a)
   n = length(b)
   n == 0 && error("Division by zero in divides_monagan_pearce")
   if m == 0
      return true, par()
   end
   H = Array{heap_sr}(0)
   I = Array{heap_t}(0)
   # set up heap
   push!(H, heap_sr(a.exps[1], 1))
   push!(I, heap_t(0, 1, 0))
   q_alloc = max(m - n, n)
   Qc = Array{T}(q_alloc)
   Qe = Array{UInt}(q_alloc)
   k = 0
   s = n
   c = R()
   qc = R()
   m1 = -R(1)
   mb = -b.coeffs[1]
   Q = Array{Int}(0)
   reuse = Array{Int}(0)
   while !isempty(H)
      exp = H[1].exp
      k += 1
      if k > q_alloc
         q_alloc *= 2
         resize!(Qc, q_alloc)
         resize!(Qe, q_alloc)
      end
      first = true
      d1 = false
      @inbounds while !isempty(H) && H[1].exp == exp
         x = H[1]
         heappop!(H)
         v = I[x.n]
         if first
            d1 = exp >= b.exps[1]
            Qe[k] = exp - b.exps[1]
            first = false
         end
         if v.i == 0
            qc = addmul!(qc, a.coeffs[v.j], m1, c)
         else
            qc = addmul!(qc, b.coeffs[v.i], Qc[v.j], c)
         end
         if v.i != 0 || v.j < m
            push!(Q, x.n)
         else
            push!(reuse, x.n)
         end
         while (xn = v.next) != 0
            v = I[xn]
            if v.i == 0
               qc = addmul!(qc, a.coeffs[v.j], m1, c)
            else
               qc = addmul!(qc, b.coeffs[v.i], Qc[v.j], c)
            end
            if v.i != 0 || v.j < m
               push!(Q, xn)
            else
               push!(reuse, xn)
            end
         end
      end
      @inbounds while !isempty(Q)
         xn = pop!(Q)
         v = I[xn]
         if v.i == 0
            I[xn] = heap_t(0, v.j + 1, 0)
            heapinsert!(H, I, xn, a.exps[v.j + 1]) # either chain or insert into heap
         elseif v.j < k - 1
            I[xn] = heap_t(v.i, v.j + 1, 0)
            heapinsert!(H, I, xn, b.exps[v.i] + Qe[v.j + 1]) # either chain or insert into heap
         elseif v.j == k - 1
            s += 1
            push!(reuse, xn)
         end
      end
      if qc == 0
         k -= 1
      else
         d2, Qc[k] = divides(qc, mb)
         if !d1 || !d2
             return false, par()
         end
         for i = 2:s
            if !isempty(reuse)
               xn = pop!(reuse)
               I[xn] = heap_t(i, k, 0)
               heapinsert!(H, I, xn, b.exps[i] + Qe[k]) # either chain or insert into heap
            else
               push!(I, heap_t(i, k, 0))
               Collections.heappush!(H, heap_sr(b.exps[i] + Qe[k], length(I)))
            end
         end
         s = 1
      end
      qc = zero!(qc)
   end
   resize!(Qc, k)
   resize!(Qe, k)
   return true, parent(a)(Qc, Qe)
end

function divides(a::SparsePoly{T}, b::SparsePoly{T}) where {T <: RingElement}
   d1 = a.exps[a.length]
   d2 = b.exps[b.length] - b.exps[1]
   q_alloc = b.length
   Qe = Array{UInt}(q_alloc)
   Qc = Array{T}(q_alloc)
   r = a
   k = 0
   while length(r) > 0
      if r.exps[1] + d2 > d1
         return false, parent(a)()
      end
      flag, c = divides(r.coeffs[1], b.coeffs[1])
      if !flag
         return false, parent(a)()
      end
      s = b*c
      if r.exps[1] < b.exps[1]
         return false, parent(a)()
      end
      d = r.exps[1] - b.exps[1]
      k += 1
      if k > q_alloc
         q_alloc *= 2
         resize!(Qe, q_alloc)
         resize!(Qc, q_alloc)
      end
      Qe[k] = d
      Qc[k] = c
      s = shift_left!(s, d)
      r -= s
   end
   resize!(Qe, k)
   resize!(Qc, k)
   return true, parent(a)(Qc, Qe)
end

function divexact(a::SparsePoly{T}, b::SparsePoly{T}) where {T <: RingElement}
   d, q = divides(a, b)
   d == false && error("Not an exact division in divexact")
   return q
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divides(a::SparsePoly{T}, b::T) where {T <: RingElem}
   len = a.length
   Qc = Array{T}(len)
   for i = 1:len
      flag, Qc[i] = divides(a.coeffs[i], b)
      if !flag
         return false, parent(a)()
      end
   end
   return true, parent(a)(Qc, a.exps)
end

function divexact(a::SparsePoly{T}, b::T) where {T <: RingElem}
   len = length(a)
   exps = deepcopy(a.exps)
   coeffs = [divexact(a.coeffs[i], b) for i in 1:len]
   return parent(a)(coeffs, exps)
end

function divexact(a::SparsePoly, b::Union{Integer, Rational, AbstractFloat})
   len = length(a)
   exps = deepcopy(a.exps)
   coeffs = [divexact(a.coeffs[i], b) for i in 1:len]
   return parent(a)(coeffs, exps)
end

###############################################################################
#
#   Pseudodivision
#
###############################################################################

function pseudodivrem(a::SparsePoly{T}, b::SparsePoly{T}) where {T <: RingElement}
   par = parent(a)
   R = base_ring(par)
   m = length(a)
   n = length(b)
   maxn = a.exps[m]
   n == 0 && error("Division by zero in pseudodivrem")
   if a.exps[m] < b.exps[n]
      return par(), a
   end
   H = Array{heap_sr}(0)
   I = Array{heap_t}(0)
   # set up heap
   push!(H, heap_sr(maxn - a.exps[m], 1))
   push!(I, heap_t(0, 1, 0))
   q_alloc = max(m - n + 1, n)
   r_alloc = n
   p_alloc = max(m - n + 1, 1)
   Qc = Array{T}(q_alloc)
   Qe = Array{UInt}(q_alloc)
   Qpow = Array{Int}(q_alloc)
   Rc = Array{T}(r_alloc)
   Re = Array{UInt}(r_alloc)
   Pc = Array{T}(p_alloc)
   Pc[1] = b.coeffs[n]
   p = -1 # current power of lead(b)
   p_max = 1
   k = 0
   l = 0
   s = n
   c = R()
   m1 = -R(1)
   mb = -b.coeffs[n]
   Q = Array{Int}(0)
   reuse = Array{Int}(0)
   while !isempty(H)
      exp = H[1].exp
      qc = R()
      k += 1
      if k > q_alloc
         q_alloc *= 2
         resize!(Qc, q_alloc)
         resize!(Qe, q_alloc)
         resize!(Qpow, q_alloc)
      end
      p += 1
      if p > p_alloc
         p_alloc *= 2
         resize!(Pc, p_alloc)
      end
      if p > p_max
         Pc[p] = Pc[p - 1]*b.coeffs[n]
         p_max = p
      end
      @inbounds while !isempty(H) && H[1].exp == exp
         x = H[1]
         heappop!(H)
         v = I[x.n]
         if v.i == 0
            # cannot use mul! here since Pc[p] is one variable less
            if p > 0
               tc = a.coeffs[m + 1 - v.j]*(-Pc[p])
            else
               tc = -a.coeffs[m + 1 - v.j]
            end
            qc = addeq!(qc, tc)
         else
            if p - Qpow[v.j] - 1 > 0
               tc = Qc[v.j]*Pc[p - Qpow[v.j] - 1]
            else
               tc = Qc[v.j]
            end
            qc = addmul!(qc, b.coeffs[n + 1 - v.i], tc, c)
         end
         if v.i != 0 || v.j < m
            push!(Q, x.n)
         else
            push!(reuse, x.n)
         end
         while (xn = v.next) != 0
            v = I[xn]
            if v.i == 0
               if p > 0
                  tc = a.coeffs[m + 1 - v.j]*(-Pc[p])
               else
                  tc = -a.coeffs[m + 1 - v.j]
               end
               qc = addeq!(qc, tc)
            else
               if p - Qpow[v.j] - 1 > 0
                  tc = Qc[v.j]*Pc[p - Qpow[v.j] - 1]
               else
                  tc = Qc[v.j]
               end
               qc = addmul!(qc, b.coeffs[n + 1 - v.i], tc, c)
            end
            if v.i != 0 || v.j < m
               push!(Q, xn)
            else
               push!(reuse, xn)
            end
         end
      end
      @inbounds while !isempty(Q)
         xn = pop!(Q)
         v = I[xn]
         if v.i == 0
            I[xn] = heap_t(0, v.j + 1, 0)
            heapinsert!(H, I, xn, maxn - a.exps[m - v.j]) # either chain or insert into heap
         elseif v.j < k - 1
            I[xn] = heap_t(v.i, v.j + 1, 0)
            heapinsert!(H, I, xn, maxn - b.exps[n + 1 - v.i] - Qe[v.j + 1]) # either chain or insert into heap
         elseif v.j == k - 1
            s += 1
            push!(reuse, xn)
         end
      end
      if qc == 0
         k -= 1
         p -= 1
      else
         d1 = maxn - b.exps[n] >= exp
         texp = maxn - b.exps[n] - exp
         if !d1
            l += 1
            if l >= r_alloc
               r_alloc *= 2
               resize!(Rc, r_alloc)
               resize!(Re, r_alloc)
            end
            Rc[l] = -qc
            Re[l] = maxn - exp
            k -= 1
            p -= 1
         else
            Qc[k] = -qc
            Qe[k] = texp
            Qpow[k] = p
            for i = 2:s
               if !isempty(reuse)
                  xn = pop!(reuse)
                  I[xn] = heap_t(i, k, 0)
                  heapinsert!(H, I, xn, maxn - b.exps[n + 1 - i] - Qe[k]) # either chain or insert into heap
               else
                  push!(I, heap_t(i, k, 0))
                  Collections.heappush!(H, heap_sr(maxn - b.exps[n + 1 - i] - Qe[k], length(I)))
               end
            end
            s = 1
         end
      end
   end
   resize!(Qc, k)
   resize!(Qe, k)
   for i = 1:k
      if p > Qpow[i]
         Qc[i] = mul!(Qc[i], Qc[i], Pc[p - Qpow[i]])
      else
         Qc[i] = Qc[i]
      end
   end
   resize!(Rc, l)
   resize!(Re, l)
   reverse!(Qc)
   reverse!(Qe)
   reverse!(Rc)
   reverse!(Re)
   if a.exps[m] - b.exps[n] > p
      k1 = reinterpret(Int, a.exps[m] - b.exps[n] - p)
      p1 = b.coeffs[n]^k1
      return parent(a)(Qc, Qe)*p1, parent(a)(Rc, Re)*p1
   else
      return parent(a)(Qc, Qe), parent(a)(Rc, Re)
   end
end

function pseudorem_monagan_pearce(a::SparsePoly{T}, b::SparsePoly{T}) where {T <: RingElement}
   par = parent(a)
   R = base_ring(par)
   m = length(a)
   n = length(b)
   maxn = a.exps[m]
   n == 0 && error("Division by zero in pseudodivrem")
   if a.exps[m] < b.exps[n]
      return a
   end
   H = Array{heap_sr}(0)
   I = Array{heap_t}(0)
   # set up heap
   push!(H, heap_sr(maxn - a.exps[m], 1))
   push!(I, heap_t(0, 1, 0))
   q_alloc = max(m - n + 1, n)
   r_alloc = n
   p_alloc = max(m - n + 1, 1)
   Qc = Array{T}(q_alloc)
   Qe = Array{UInt}(q_alloc)
   Qpow = Array{Int}(q_alloc)
   Rc = Array{T}(r_alloc)
   Re = Array{UInt}(r_alloc)
   Pc = Array{T}(p_alloc)
   Pc[1] = b.coeffs[n]
   p = -1 # current power of lead(b)
   p_max = 1
   k = 0
   l = 0
   s = n
   c = R()
   m1 = -R(1)
   mb = -b.coeffs[n]
   Q = Array{Int}(0)
   reuse = Array{Int}(0)
   while !isempty(H)
      exp = H[1].exp
      qc = R()
      k += 1
      if k > q_alloc
         q_alloc *= 2
         resize!(Qc, q_alloc)
         resize!(Qe, q_alloc)
         resize!(Qpow, q_alloc)
      end
      p += 1
      if p > p_alloc
         p_alloc *= 2
         resize!(Pc, p_alloc)
      end
      if p > p_max
         Pc[p] = Pc[p - 1]*b.coeffs[n]
         p_max = p
      end
      @inbounds while !isempty(H) && H[1].exp == exp
         x = H[1]
         heappop!(H)
         v = I[x.n]
         if v.i == 0
            # cannot use mul! here since Pc[p] is one variable less
            if p > 0
               tc = a.coeffs[m + 1 - v.j]*(-Pc[p])
            else
               tc = -a.coeffs[m + 1 - v.j]
            end
            qc = addeq!(qc, tc)
         else
            if p - Qpow[v.j] - 1 > 0
               tc = Qc[v.j]*Pc[p - Qpow[v.j] - 1]
            else
               tc = Qc[v.j]
            end
            qc = addmul!(qc, b.coeffs[n + 1 - v.i], tc, c)
         end
         if v.i != 0 || v.j < m
            push!(Q, x.n)
         else
            push!(reuse, x.n)
         end
         while (xn = v.next) != 0
            v = I[xn]
            if v.i == 0
               if p > 0
                  tc = a.coeffs[m + 1 - v.j]*(-Pc[p])
               else
                  tc = -a.coeffs[m + 1 - v.j]
               end
               qc = addeq!(qc, tc)
            else
               if p - Qpow[v.j] - 1 > 0
                  tc = Qc[v.j]*Pc[p - Qpow[v.j] - 1]
               else
                  tc = Qc[v.j]
               end
               qc = addmul!(qc, b.coeffs[n + 1 - v.i], tc, c)
            end
            if v.i != 0 || v.j < m
               push!(Q, xn)
            else
               push!(reuse, xn)
            end
         end
      end
      @inbounds while !isempty(Q)
         xn = pop!(Q)
         v = I[xn]
         if v.i == 0
            I[xn] = heap_t(0, v.j + 1, 0)
            heapinsert!(H, I, xn, maxn - a.exps[m - v.j]) # either chain or insert into heap
         elseif v.j < k - 1
            I[xn] = heap_t(v.i, v.j + 1, 0)
            heapinsert!(H, I, xn, maxn - b.exps[n + 1 - v.i] - Qe[v.j + 1]) # either chain or insert into heap
         elseif v.j == k - 1
            s += 1
            push!(reuse, xn)
         end
      end
      if qc == 0
         k -= 1
         p -= 1
      else
         d1 = maxn - b.exps[n] >= exp
         texp = maxn - b.exps[n] - exp
         if !d1
            l += 1
            if l >= r_alloc
               r_alloc *= 2
               resize!(Rc, r_alloc)
               resize!(Re, r_alloc)
            end
            Rc[l] = -qc
            Re[l] = maxn - exp
            k -= 1
            p -= 1
         else
            Qc[k] = -qc
            Qe[k] = texp
            Qpow[k] = p
            for i = 2:s
               if !isempty(reuse)
                  xn = pop!(reuse)
                  I[xn] = heap_t(i, k, 0)
                  heapinsert!(H, I, xn, maxn - b.exps[n + 1 - i] - Qe[k]) # either chain or insert into heap
               else
                  push!(I, heap_t(i, k, 0))
                  Collections.heappush!(H, heap_sr(maxn - b.exps[n + 1 - i] - Qe[k], length(I)))
               end
            end
            s = 1
         end
      end
   end
   resize!(Rc, l)
   resize!(Re, l)
   reverse!(Rc)
   reverse!(Re)
   if a.exps[m] - b.exps[n] > p
      k1 = reinterpret(Int, a.exps[m] - b.exps[n] - p)
      p1 = b.coeffs[n]^k1
      return parent(a)(Rc, Re)*p1
   else
      return parent(a)(Rc, Re)
   end
end

function pseudorem(a::SparsePoly{T}, b::SparsePoly{T}) where {T <: RingElement}
   par = parent(a)
   R = base_ring(par)
   m = length(a)
   n = length(b)
   n == 0 && throw(DivideError())
   if a.exps[a.length] < b.exps[b.length]
      return deepcopy(a)
   end
   k = reinterpret(Int, a.exps[a.length] - b.exps[b.length]) + 1
   l = lead(b)
   while a.length > 0 && a.exps[a.length] >= b.exps[b.length]
      s = lead(a)*b
      s = shift_left!(s, a.exps[a.length] - b.exps[b.length])
      a = a*l - s
      k -= 1
   end
   return a*l^k
end

###############################################################################
#
#   Evaluation
#
###############################################################################

function evaluate(a::SparsePoly{T}, b::S) where {S <: RingElem, T <: RingElement}
   if a.length == 0
      return base_ring(a)()
   end
   r = a.coeffs[a.length]
   for i = 1:a.length - 1
      r *= b^(reinterpret(Int, a.exps[a.length - i + 1] - a.exps[a.length - i]))
      r += a.coeffs[a.length - i]
   end
   if a.exps[1] != 0
      r *= b^(reinterpret(Int, a.exps[1]))
   end
   return r
end

function evaluate(a::SparsePoly{T}, b::Rational{S}) where {S <: Integer, T <: RingElement}
   if a.length == 0
      return base_ring(a)()
   end
   r = a.coeffs[a.length]
   for i = 1:a.length - 1
      r *= b^(reinterpret(Int, a.exps[a.length - i + 1] - a.exps[a.length - i]))
      r += a.coeffs[a.length - i]
   end
   if a.exps[1] != 0
      r *= b^(reinterpret(Int, a.exps[1]))
   end
   return r
end

function evaluate(a::SparsePoly{T}, b::Integer) where {T <: RingElement}
   if a.length == 0
      return base_ring(a)()
   end
   R = base_ring(a)
   r = a.coeffs[a.length]
   for i = 1:a.length - 1
      r *= R(b)^(reinterpret(Int, a.exps[a.length - i + 1] - a.exps[a.length - i]))
      r += a.coeffs[a.length - i]
   end
   if a.exps[1] != 0
      r *= R(b)^(reinterpret(Int, a.exps[1]))
   end
   return r
end

###############################################################################
#
#   GCD, content and primitive part
#
###############################################################################

function gcd(a::SparsePoly{T}, b::SparsePoly{T}, ignore_content::Bool = false) where {T <: RingElement}
   # ensure degree in main variable of a is at least that of b
   if b.exps[b.length] > a.exps[a.length]
      (a, b) = (b, a)
   end
   if b == 0
      return deepcopy(a)
   end
   if b == 1
      return deepcopy(b)
   end
   # compute gcd of contents and divide content out
   if !ignore_content
      c1 = content(a)
      c2 = content(b)
      c = gcd(c1, c2)
      a = divexact(a, c1)
      b = divexact(b, c2)
   end
   # check if we are in the univariate case
   constant_coeffs = true
   for i = 1:a.length
      if !isconstant(a.coeffs[i])
         constant_coeffs = false
         break
      end
   end
   if constant_coeffs
      for i = 1:b.length
         if !isconstant(b.coeffs[i])
            constant_coeffs = false
            break
         end
      end
   end
   # if we are in univariate case, convert to dense, take gcd, convert back
   if constant_coeffs
      # convert polys to univariate dense
      R, x = PolynomialRing(base_ring(base_ring(a)), "\$")
      f = R()
      g = R()
      fit!(f, reinterpret(Int, a.exps[a.length] + 1))
      fit!(g, reinterpret(Int, b.exps[b.length] + 1))
      for i = 1:a.length
         f = setcoeff!(f, reinterpret(Int, a.exps[i]), a.coeffs[i].coeffs[1])
      end
      for i = 1:b.length
         g = setcoeff!(g, reinterpret(Int, b.exps[i]), b.coeffs[i].coeffs[1])
      end
      # take gcd of univariate dense polys
      h = gcd(f, g)
      # convert back to sparse polys
      nonzero = 0
      Ac = Array{T}(0)
      Ae = Array{UInt}(0)
      for i = 1:length(h)
         ci = coeff(h, i - 1)
         if ci != 0
            push!(Ac, base_ring(a)(ci))
            push!(Ae, UInt(i - 1))
         end
      end
      r = parent(a)(Ac, Ae)
      if !ignore_content
         return c*r
      else
         return r
      end
   end
   # compute likely degree of gcd
   deg = 0
   # is the lead/trail term a monomial
   lead_monomial = lead(a).length == 1 || lead(b).length == 1
   trail_monomial = trail(a).length == 1 || trail(b).length == 1
   lead_a = lead(a)
   lead_b = lead(b)
   # psr algorithm
   g = one(base_ring(a))
   h = one(base_ring(a))
   while true
      adeg = reinterpret(Int,  a.exps[a.length])
      bdeg = reinterpret(Int,  b.exps[b.length])
      # optimisation: if degree of b is equal to likely degree of gcd, try
      # exact division, learned from Bernard Parisse
      if bdeg == deg
         flag, q = divides(a, b)
         if flag
            break
         end
      end
      d = reinterpret(Int, adeg - bdeg)
      r = pseudorem(a, b)
      # zero remainder
      if r == 0
         break
      end
      # constant remainder
      if r.length == 1 && r.exps[1] == 0
         b = one(parent(a))
         break
      end
      (a, b) = (b, divexact(r, g*h^d))
      g = lead(a)
      if d > 1
         h = divexact(g^d, h^(d - 1))
      else
         h = h^(1 - d)*g^d
      end
   end
   # sometimes don't care about content, e.g. when computing likely gcd degree
   if !ignore_content
      # remove content from b as cheaply as possible, as per Bernard Parisse
      if lead(b).length != 1 && trail(b).length != 1
         if lead_monomial # lead term monomial, so content contains rest
            d = divexact(lead(b), term_content(lead(b)))
            b = divexact(b, d)
         elseif trail_monomial # trail term is monomial, so ditto
            d = divexact(trail(b), term_content(trail(b)))
            b = divexact(b, d)
         else
            glead = gcd(lead_a, lead_b)
            if glead.length == 1 # gcd of lead coeffs monomial
               d = divexact(lead(b), term_content(lead(b)))
               b = divexact(b, d)
            else # last ditched attempt to find easy content
               h = gcd(lead(b), glead)
               h = divexact(h, term_content(h))
               flag, q = divides(b, h)
               if flag
                  b = q
               end
            end
         end
      end
      # remove any monomial content
      b1 = term_content(b.coeffs[1])
      for i = 2:b.length
         b1 = gcd(b1, term_content(b.coeffs[i]))
         if isone(b1)
            break
         end
      end
      b = divexact(b, b1)
      # remove any stubborn content and put back actual content
      return c*primpart(b)
   else
      return b
   end
end

function content(a::SparsePoly{T}) where {T <: RingElement}
   for i = 1:length(a)
      if a.coeffs[i].length == 1
         z = term_content(a.coeffs[1])
         for j = 2:length(a)
            if isone(z)
               return z
            end
            z = gcd(z, term_content(a.coeffs[j]))
         end
         return z
      end
   end
   z = coeff(a, 0)
   for i = 2:length(a)
      if z == 1
         break
      end
      z = gcd(coeff(a, i - 1), z)
   end
   return z
end

function primpart(a::SparsePoly{T}) where {T <: RingElement}
   d = content(a)
   return divexact(a, d)
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function fit!(a::SparsePoly{T}, n::Int) where {T <: RingElement}
   if length(a.coeffs) < n
      resize!(a.coeffs, n)
      resize!(a.exps, n)
   end
   return nothing
end

function addmul!(a::SparsePoly{T}, b::SparsePoly{T}, c::SparsePoly{T}, d::SparsePoly{T}) where {T <: RingElement}
   t = b*c
   t += a
   a.coeffs = t.coeffs
   a.exps = t.exps
   a.length = t.length
   return a
end

function mul!(a::SparsePoly{T}, b::SparsePoly{T}, c::SparsePoly{T}) where {T <: RingElement}
   t = b*c
   a.coeffs = t.coeffs
   a.exps = t.exps
   a.length = t.length
   return a
end

function addeq!(a::SparsePoly{T}, b::SparsePoly{T}) where {T <: RingElement}
   t = a + b
   a.coeffs = t.coeffs
   a.exps = t.exps
   a.length = t.length
   return a
end

function add!(a::SparsePoly{T}, b::SparsePoly{T}, c::SparsePoly{T}) where {T <: RingElement}
   t = b + c
   a.coeffs = t.coeffs
   a.exps = t.exps
   a.length = t.length
   return a
end

function zero!(a::SparsePoly{T}) where {T <: RingElement}
   a.length = 0
   return a
end

function shift_left!(a::SparsePoly{T}, n::UInt) where {T <: RingElement}
   for i = 1:length(a)
      a.exps[i] += n
   end
   return a
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule(::Type{SparsePoly{T}}, ::Type{SparsePoly{T}}) where T <: RingElement = SparsePoly{T}

function promote_rule(::Type{SparsePoly{T}}, ::Type{U}) where {T <: RingElement, U <: RingElement}
   promote_rule(T, U) == T ? SparsePoly{T} : Union{}
end

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::SparsePolyRing{T} where {T <: RingElement})(b::RingElement)
   return a([base_ring(a)(b)], [UInt(0)])
end

function (a::SparsePolyRing{T})() where {T <: RingElement}
   z = SparsePoly{T}()
   z.parent = a
   return z
end

function (a::SparsePolyRing{T})(b::Union{Integer, Rational, AbstractFloat}) where {T <: RingElement}
   z = SparsePoly{T}(base_ring(a)(b))
   z.parent = a
   return z
end

function (a::SparsePolyRing{T})(b::T) where {T <: RingElement}
   parent(b) != base_ring(a) && error("Unable to coerce to polynomial")
   z = SparsePoly{T}(b)
   z.parent = a
   return z
end

function (a::SparsePolyRing{T})(b::SparsePoly{T}) where {T <: RingElement}
   parent(b) != a && error("Unable to coerce polynomial")
   return b
end

function (a::SparsePolyRing{T})(b::Array{T, 1}, m::Array{UInt, 1}) where {T <: RingElement}
   if length(b) > 0
      parent(b[1]) != base_ring(a) && error("Unable to coerce to polynomial")
   end
   z = SparsePoly{T}(b, m)
   z.parent = a
   return z
end

###############################################################################
#
#   PolynomialRing constructor
#
###############################################################################

doc"""
    SparsePolynomialRing(R::Nemo.Ring, s::String; cached::Bool = true)
> Given a base ring `R` and a string `s` specifying how the generator
> (variable) should be printed, return a tuple `S, x` representing the new
> polynomial ring $T = R[x1, x2, ...]$ and the generator $x$ of the polynomial
> ring. By default the parent object `T` will depend only on `R` and `x` and
> will be cached. Setting the optional argument `cached` to `false` will
> prevent the parent object `T` from being cached.
"""
function SparsePolynomialRing(R::Nemo.Ring, s::String; cached::Bool = true)
   U = Symbol(s)
   T = elem_type(R)

   par = SparsePolyRing{T}(R, U, cached)

   return par, gen(par)
end


text:
module PCM

using CxxWrap

include("lib.jl")
using .Lib

reset() = Lib.resetPMU(Lib.getInstance())
cleanup() = Lib.cleanup(Lib.getInstance())

const MAX_CUSTOM_EVENTS = convert(Int, Lib.getMaxCustomCoreEvents(Lib.getInstance()))

#####
##### General Queries
#####

pcm() = Lib.getInstance()
numcores() = convert(Int, Lib.getNumCores(pcm()))

#####
##### Core Counters
#####

struct EventDescription
    event_number::Int64
    umask_value::Int64
    name::Symbol

    function EventDescription(event::Integer, umask::Integer, name = :unknown_event)
        return new(
            convert(Int64, event),
            convert(Int64, umask),
            name
        )
    end
end

function Base.show(io::IO, E::EventDescription)
    print(
        io,
        "PCM Event \'$(E.name)\'-- event: 0x",
        string(E.event_number; base = 16),
        ", umask: 0x",
        string(E.umask_value; base = 16),
    )
    return nothing
end

mutable struct CoreMonitor
    pcm::Lib.PCMRef
    wrapper::Lib.CounterWrapperAllocated
    cores::Vector{Int64}
    events::Vector{EventDescription}

    # Inner constructore - need to get the reference to the singular PCM instance as well
    # as instantiate a `CounterWrapper` for storing counter values.
    function CoreMonitor(; cores::AbstractVector = 1:numcores())
        pcm = Lib.getInstance()
        wrapper = Lib.CounterWrapper()

        # Do some processing on the on the `cores` argument - make sure all values are
        # within the valid core range, materialize it etc.
        cores = sort(Vector(cores))
        if first(cores) < 0
            throw(ArgumentError("Core numbers must be greater than zero"))
        end

        if last(cores) > numcores()
            throw(ArgumentError("Core numbers must be less than or equal to $(numcores())"))
        end

        # Subtract 1 from everything since we're using this to talk with the C++ code.
        cores .= cores .- 1

        monitor = new(
            pcm,
            wrapper,
            cores,
            EventDescription[],
        )

        return monitor
    end
end

function program(monitor::CoreMonitor, events::Vector{EventDescription})
    original_size = length(events)

    # If more than MAX_CUSTOM_EVENTS events are provided, throw an error.
    if length(events) > MAX_CUSTOM_EVENTS
        err = ArgumentError("Number of events must be less than or equal to 4!")
        throw(err)
    end

    # Fill out vector of event descriptions if not enough are provided.
    # This needs to be of length 4 because PCM is not the best coded thing in the world.
    while length(events) < 4
        push!(events, EventDescription(0, 0))
    end

    # Convert the array of structs into two arrays
    event_numbers = [i.event_number for i in events]
    umask_values = [i.umask_value for i in events]
    resize!(events, original_size)

    # Program the event counters.
    err = Lib.program(monitor.pcm, event_numbers, umask_values)

    if err == Lib.PCM_MSRAccessDenied
        throw(error("Cannot Access MSR"))
    elseif err == Lib.PCM_PMUBusy
        throw(error("PMU is Busy"))
    elseif err == Lib.PCM_UnknownError
        throw(error("Unknow Error!"))
    end

    monitor.events = events

    # Sample twice to initialize the counter arrays
    sample!(monitor)
    sample!(monitor)

    return err
end

# Swap the double-buffered counters and store the current counter state into the `after`
# set of counters.
sample!(monitor::CoreMonitor) = Lib.sample(monitor.wrapper, monitor.pcm)

# Get the values of the counters in `monitor`.
function getcounters(monitor)
    cores = monitor.cores
    values = Lib.aggregate_counters(monitor.wrapper, monitor.pcm, cores, length(cores))

    return collect(values)[1:length(monitor.events)]
end

#####
##### Uncore Counters
#####

# Type for measuring Memory and PMM bandwidth
mutable struct UncoreMemoryMonitor
    pcm::Lib.PCMRef

    # Keep vectors of UncorePowerStates, one for each socket.
    # We will swap these back and forth to get the deltas of the metrics we are measuring.
    before::Vector{Lib.ServerUncorePowerState}
    after::Vector{Lib.ServerUncorePowerState}

    # Some statistics for performance

    # Number of memory controllers
    num_memory_controllers::Int64
    num_memory_channels::Int64
    num_sockets::Int64
end

function UncoreMemoryMonitor(num_memory_controllers, num_memory_channels)
    # Instantiate the PCM struct and program the memory controller counters
    pcm = Lib.getInstance()
    Lib.programServerUncoreMemoryMetrics(pcm, -1, -1, true)

    # Instantiate the dual arrays for sampling and do a quick sample so
    before = Vector{Lib.ServerUncorePowerState}()
    after = Vector{Lib.ServerUncorePowerState}()
    _sample!(pcm, before)
    _sample!(pcm, after)

    obj = UncoreMemoryMonitor(
        pcm,
        before,
        after,
        num_memory_controllers,
        num_memory_channels,
        Lib.getNumSockets(pcm)
    )
end

function sample!(U::UncoreMemoryMonitor)
    # Swap before and after state.
    U.before, U.after = U.after, U.before

    # Fill the new `after` state with the values of the counters.
    _sample!(U.pcm, U.after)
    return nothing
end

function _sample!(m::Lib.PCMRef, counters::Vector{Lib.ServerUncorePowerState})
    empty!(counters)
    for i in 1:Lib.getNumSockets(m)
        push!(counters, Lib.getServerUncorePowerState(m, i-1))
    end
end

function getcounter(U::UncoreMemoryMonitor, counter, socket)
    return getcounter(U.before, U.after, counter, U.num_memory_channels, socket + 1)
end

function getcounter(before, after, counter, nchannels, socket)
    return map(1:nchannels) do i
        x = Lib.getMCCounter(
            i-1,                    # memory channel
            convert(Int, counter),  # counter number
            before[socket],         # before uncore counter state
            after[socket],          # after uncore counter state
        )

        # normalize return value
        return convert(Int, x)
    end
end

end # module


text:
srand(0)

performancetests = []

println("==========================")
println("Running performance tests:")
println("==========================")

for t in performancetests
  tfile = t*".jl"
  println("  * $(tfile) *")
  include(string("performance/",tfile))
end


text:
using Test
using NQCDynamics
using Distributions: Normal

@testset "ElectronicDistribution" begin
    diabatic = NonadiabaticDistributions.Diabatic()

    @testset "SingleState" begin
        @test NonadiabaticDistributions.SingleState(1) === NonadiabaticDistributions.SingleState(1, diabatic)
    end

    @testset "ElectronicPopulation" begin
        dist = NonadiabaticDistributions.ElectronicPopulation([1, 0])
        dist2 = NonadiabaticDistributions.ElectronicPopulation([1, 0], diabatic)
        @test dist.statetype === dist2.statetype
        @test dist.populations == dist2.populations
    end
end

@testset "CombinedDistribution" begin
    electronic = NonadiabaticDistributions.SingleState(1)
    nuclear = DynamicalDistribution(Normal(), Normal(), (3, 2))
    combined = electronic * nuclear
    @test combined isa NonadiabaticDistributions.NonadiabaticDistribution
end


text:
module ApproxInference

const VecI  = AbstractVector # Input  Vector
const VecO  = AbstractVector # Output Vector
const VecB  = AbstractVector # Buffer Vector
const VecIO = AbstractVector # In/Out Vector
const MatI  = AbstractMatrix # Input  Matrix
const MatO  = AbstractMatrix # Output Matrix
const MatB  = AbstractMatrix # Buffer Matrix
const MatIO = AbstractMatrix # In/Out Matrix

import LinearAlgebra: BLAS
import LinearAlgebra: LAPACK

if haskey(ENV, "BLAS_THREAD_NUM")
    BLAS.set_num_threads(parse(Int, ENV["BLAS_THREAD_NUM"]))
else
    BLAS.set_num_threads(4)
end

include("./linalg.jl")
include("./macros.jl")
include("./statistics.jl")

abstract type AbstractObjective <: Function end

macro objective(e::Expr)
    e.head ≡ :struct || error("Objective funtion should be defined as `struct ... end`.")
    @inbounds a = e.args[2]
    if a isa Expr && a.head ≡ :<:
        @inbounds a.args[2] = :AbstractObjective
    else
        @inbounds e.args[2] = Expr(:<:, a, :AbstractObjective)
    end
    return e
end

function optimize!     end
function get_residual! end
function get_jacobian! end

include("./lmpack.jl")
include("./vipack.jl")

optimizer(; nd::Int=1, ny::Int=1, args...) = optimizer(nd, ny; args...)
function optimizer(nd::Int, ny::Int; method::String="lm")
    if method ≡ "lm" || method ≡ "LM" || method ≡ "LevenbergMarquardt" || method ≡ "Levenberg-Marquardt"
        return LevenbergMarquardtOptimizer(nd, ny)
    elseif method ≡ "vi" || method ≡ "VI" || method ≡ "VariationalInference" || method ≡ "Variational-Inference"
        return VariationalInferenceOptimizer(nd, ny)
    end
end

export @objective, optimize!, get_residual!, get_jacobian!, optimizer

end # module


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ f493fdb0-022f-11eb-34a6-4b133c23a71a
using OrdinaryDiffEq

# ╔═╡ c803fe10-0231-11eb-1f84-dbe43fd97eb4
using Plots

# ╔═╡ be8523c0-0239-11eb-0439-592e370c5618
using LaTeXStrings #rotulos de las figuras bonitos

# ╔═╡ 5fe018a0-023b-11eb-1e69-d59b562edc9b
using Symbolics
#alternativamente también se puede usar SymEngine pero hay que definir el jacobiano a mano

# ╔═╡ 1c0442a0-a993-11eb-1b2d-55eb8b868dc7
#algunas definiciones útiles
begin
	using Latexify  # si no lo tenéis hay que instalarlo
    Base.show(io::IO, ::MIME"text/html", x::Symbolics.Num) = print(io, latexify(x))
    Base.show(io::IO, ::MIME"text/html", x::Symbolics.Symbolic) = print(io, latexify(x))
    Base.show(io::IO, ::MIME"text/html", x::Symbolics.Vector{Equation}) = print(io, latexify(x))
    Base.show(io::IO, ::MIME"text/html", x::Symbolics.AbstractArray{Num}) = print(io, latexify(x))

	diff(f,variable) = expand_derivatives(Symbolics.derivative(f,variable))
	subs(expr,dict) = substitute(expr, dict)
	expand(x) = simplify(x;polynorm=true) #En ciertas versiones es expand = true
	solve_for(ecuacion,variable) = Symbolics.solve_for(ecuacion, variable)
	jacobian(f,x) = Symbolics.jacobian(f,x)
end

# ╔═╡ 183a6460-a9ca-11eb-340e-872e08f29a9e
using LinearAlgebra

# ╔═╡ c900ad60-095f-11eb-3d18-9f187abdb9b3
using PlutoUI

# ╔═╡ 5c4af6e0-022e-11eb-152f-2153b9b672c1
md"""# Ejemplos de sistemas a controlar

Primer ejemplo

$$J \ddot \theta + MgLsin(\theta)=u$$

"""

# ╔═╡ 6b1fd530-0a28-11eb-1142-47afd1798c3e
html"<button onclick=present()>Present</button>"

# ╔═╡ 85c03440-0274-11eb-3c81-514259334750
md"## Simulamos y pintamos"

# ╔═╡ e3e9d0a0-022c-11eb-04ba-97dcb9ba908a
md"""# Grado relativo
La idea es "derivar hasta que sale u"
"""

# ╔═╡ 6831eb02-023b-11eb-0842-fbe8f8eb8937
md"Vamos a usar cálculo simbólico"

# ╔═╡ 3ff84c30-0248-11eb-33ba-f99304ee8672
begin
	D=x1=x2=t=u=J=M=g=L=0  #eso es cosa de pluto y las macros (si no no es necesario)
	@variables D x1 x2 t u J M g L  # =symbols("D,x1,x2,t,u,J, M, g, L")
end

# ╔═╡ 6d535e50-a99a-11eb-08c1-4d98301e53b3
begin
	#escribir las ecuaciones simbólicas
	x=[x1;x2]
	dx1=x2
	dx2=-M*g*L/J*sin(x1)+u/J
	dx=[dx1;dx2]
end

# ╔═╡ 5525d1a0-a9c7-11eb-1de0-1700c9c20460
y=x1

# ╔═╡ 65893910-a9c7-11eb-2f9f-ebadf375cdbc
dy=diff(y,x1)*dx1+diff(y,x2)*dx2

# ╔═╡ a0571620-a9c7-11eb-070f-297fa9e406b2
ddy=diff(dy,x1)*dx1+diff(dy,x2)*dx2

# ╔═╡ 39ce8af0-023d-11eb-3178-857433b369ae
md"ya aparece u --> grado relativo 2"

# ╔═╡ 4c478bf0-023d-11eb-0dce-439d408010fa
md"¿Es repetitivo podemos hacerlo en general?
 hagamos una **función**"

# ╔═╡ b7e08f60-023d-11eb-13cb-ad5285e900d8
begin
	function derivada(funcion,variables,derivadas)
		d=0
		for i=1:length(variables)
			d=d+diff(funcion,variables[i])*derivadas[i]
		end
		return d
	end

	function derivada(f,x,dx,n) #multiple distpach #esto es nuevo
		d=f
		for i=1:n
			d=derivada(d,x,dx)
		end
		return d
	end
end

# ╔═╡ 697f24d0-a99a-11eb-26e6-a749fb04cb0f
#usar la función
der=derivada(y,x,dx,2)

# ╔═╡ 9ecc5ff0-022d-11eb-3715-51831818d04f
md"""# Cambio de variables
vamos a cambiar u por v para ello despejamos de
$$y^{(\gamma)}=v$$
"""

# ╔═╡ 65ad8360-a99a-11eb-3764-aded1015107e
begin
	#definir v
	v=0
	@variables v
end

# ╔═╡ 4b753d20-a9c8-11eb-320a-d58ec5e3c431
ecuación=der~v

# ╔═╡ 6e31c180-a9c8-11eb-3893-f55415fa4fad
uc=solve_for(ecuación,u)

# ╔═╡ 561e07e2-0250-11eb-1380-4b5dfa861cd6
md"vamos a ver que lo hemos hecho bien"

# ╔═╡ 64db3f40-a99a-11eb-01df-d3ca4846ca33
#sustituir u en la derivada y ver que sale v
expand(subs(der,u=>uc))

# ╔═╡ de497f50-0250-11eb-028e-ffe56cf1ba45
md"La idea es controlar un nuevo sistema:

$$\dot z_1=z_2$$

$$...$$

$$\dot z_{\gamma}=v$$

$$y=z_1$$
"

# ╔═╡ c2f97c40-a99a-11eb-30e7-359effec0a1d
#define z1 y su derivada z2
z1=y

# ╔═╡ 85e1f1f0-a9c9-11eb-2a97-1944f7fcb34b
z2=derivada(y,x,dx)

# ╔═╡ 0db8c9c0-0252-11eb-199c-2d69997c6769
md"veamos que lo hemos hecho bien y la derivada es la corecta"

# ╔═╡ dd56c0c2-a99a-11eb-0c87-8f8e2c6321a8
expand(subs(derivada(z2,x,dx), u=>uc))

# ╔═╡ 72e4e740-022d-11eb-3a63-31023b9d464f
md"""# Difeomorfismo
A partir de ahora controlaremos las **z** pero *OJO* que nuestro objetivo son las **x**

¿Al controlar z **se controla x**?

-Si se puede despejar x(x) entonces si-->difeomorfismo (al menos local)

Para verlo calculamos el jacobiano y vemos si es singular

$$\begin{pmatrix}

\frac{\partial z_1}{\partial x_1}... & \frac{\partial z_1}{\partial x_n} \\
...\\
\frac{\partial z_n}{\partial x_1}... & \frac{\partial z_n}{\partial x_n}

\end{pmatrix}$$

"""

# ╔═╡ f1e0ee80-a99a-11eb-0c06-85f8d7086aac
#definir el vector z, calcular el jacobiano y derivar
z=[z1,z2]

# ╔═╡ 03abf0e0-a9ca-11eb-339c-5599c951871c
jac=jacobian(z,x)

# ╔═╡ 2d9109e0-a9ca-11eb-08a5-7596eee9e4b2
det(jac)

# ╔═╡ 9ec274e0-022d-11eb-085a-7907b8f1eb93
md"""# Estabilización
El sitema en las zs es lineal
Para estabilizar la dinámica se realimentan negativamente las zs con ganancias
$$v=-d_0z_1-d_1z_2-...-d_{\gamma-1}z_{\gamma}$$

Recordemos que $$\dot z_1=z_2$$, $$\ddot z_1= \dot z_2=z_3$$...
$$z_1^{(\gamma)}=v=-d_0z_1-d_1z_2-...-d_{\gamma-1}z_{\gamma}$$

$$z_1^{(\gamma)} +d_{\gamma-1}z^{\gamma-1}+...+d_1\dot z+ d_0z  =0$$

Que es un sistema lineal con polinomio característico

$$s^{(\gamma)} +d_{\gamma-1}s+...+d_1s+ d_0  =0$$

## Elección de polos

Hay que hacer que el polinomio sea hurwitz
Lo más fácil es constuirlo, si queremos unos polos $p_1$, $p_2$... $p_{\gamma}$ podemos hacer que el polinomio sea

$$(s-p_1)(s-p_2)...(s-p_{\gamma})$$
"""

# ╔═╡ aed8d350-0265-11eb-1223-1dd0651632b6
begin
	#expandir (s-p1)*(s-p2)) y comparar con lo deseado
	p1=p2=-1
	s=0
	@variables s
	expand((s-p1)*(s-p2))
end

# ╔═╡ 0d736c40-0266-11eb-3433-ad75c1617314
md"Entonces $d_0$ y $d_1$ son 1 y 2 ;)"

# ╔═╡ 454272d0-a9cb-11eb-36db-ffb885a95a52
string(uc)

# ╔═╡ f54ff5b2-0266-11eb-1df3-49f3c95729c7
md"Ahora lo ponemos todo junto"

# ╔═╡ 2814f040-0267-11eb-3f4d-2933c0c3528a
function control_1_estabilizacion(x,p,t)
	x1,x2=x
	M,g,L,J=p
	d0=1
	d1=2

	y=x1 #o si lo preferís z1
	dy=x2 #z2
	v=-d0*y-d1*dy
	u=J*v + L*M*g*sin(x1)
	return u
end

# ╔═╡ 3b125190-0269-11eb-2473-33cbb7451779
md"Si queremos podemos sustituir todo **ojo a veces no compensa** (por que queda muy largo y difícil de depurar) en este caso si"

# ╔═╡ 639dce00-0269-11eb-191e-d16a91f3bcf9
u_sustituido=subs(uc,[v=>-2*y-3*dy,L=>1,M=>1,g=>1,J=>1])

# ╔═╡ e0b21ca0-0a2c-11eb-0e50-efa6c355e7c5
md"Pero **ojo** se pierde la interpretación física de lo que hace el control"

# ╔═╡ f162f3f0-022d-11eb-30d9-55fb6e39c8c0
md"""# Seguimiento
Una vez linealizado seguir una referencia es fácil, si definimos errores

$$e_0=y_r-y$$

$$\dot e_0= \dot y_r - \dot y = e_1$$

...

$$e_0^{(\gamma)}=y_r^{(\gamma)} - v$$

Es todo igual pero añadiendo $y_r$ y con el signo contrario (aparece -v en vez de v) es decir:

$$v=d_0e_0+d_1e_1+...+d_{\gamma-1}e_{\gamma-1} + y_r^{(\gamma)}$$

**OJO a la notación** no solo cambian los signos sino que *los índices ahora están "bien"* ya que $y=e_0$ no como antes que era $z_1$
"""

# ╔═╡ a7745f70-026b-11eb-1abe-d1944165d231
function referencia_1_a(t)
	yr=sin(t)
	dyr=cos(t)
	ddyr=-sin(t)
	return [yr,dyr,ddyr]
end

# ╔═╡ a0e6f960-026b-11eb-23c9-410e6d61c4a1
function control_1_seguimiento(x,p,t)
	x1,x2=x
	yr,dyr,ddyr=referencia_1_a(t)

	M,g,L,J=p

	d0=2
	d1=3

	y=x1
	dy=x2

	e0=yr-y #esto es lo nuevo
	e1=dyr-dy
	v=d0*e0 + d1*e1 + ddyr #y aquí también aparece la referencia y cambian los signos
	#v=-d0*y-d1*dy esto es lo de antes
	#v=-d0*z1-d1*z2 o en términos de las zetas

	u=J*(v + L*M*g*sin(x1)/J)
	return u
end

# ╔═╡ 506a20d0-026f-11eb-2bde-f94724bd27f3
begin
	#control_1(x,p,t)=0
	#referencia_1(t)=0

	#control_1=control_1_estabilizacion

	control_1=control_1_seguimiento
	referencia_1=referencia_1_a
end

# ╔═╡ 85c0d350-022e-11eb-0873-c118ffde014c
function derivadas_ejemplo1(x,p,t)
	x1,x2=x
	M,g,L,J=p
	u=control_1(x,p,t)#calculamos el control
	dx1=x2
    dx2=-M*g*L/J*sin(x1)+u/J
	dx=[dx1 dx2] #no hace falta return dx
end

# ╔═╡ 922bccde-022d-11eb-17f2-270aae61f923
md"""# Dinámica cero
El ejemplo anterior funciona muy bien por dos motivos:
1) Al derivar obtengo tantos estados como tenía antes (grado relativo máximo)
2) El cambio de variabes de los estados a las derivadas de y era invertible (difeomorfismo)
Eso significa que controlar $y$ y sus derivadas $y, \dot y ... y^{(\gamma)}$ es lo mismo que controlar $x_1...x_\gamma$.

Pero ¿y si no es así? ¿Qué pasa si u aparece "demasiado pronto"?

"""

# ╔═╡ cc0e02e2-04c9-11eb-2e08-b906aaec3959
md"""# Otro ejemplo

Segundo ejemplo

$$\dot x_1=-ax_1 + e^{2x_2}u$$
$$\dot x_2=2x_1x_2 +sin(x_2) + \frac{u}{2} $$
$$\dot x_3=2x_2$$
$$y=x_3$$
"""

# ╔═╡ 8839fb60-095a-11eb-063e-85fd960ea83a
md"Ya sabemos hacerlo, vamos al grano"

# ╔═╡ 54523da0-0958-11eb-2730-cbef0823f684
begin
	a=x3=0
	@variables a x3
	#copia de las ecuaciones poniendo Bs para no pisar variables anteriores
	dxb1=-x1 + exp(2x2)*u
    dxb2=2x1*x2 +sin(x2) + u/2
	dxb3=2x2
	dxb=[dxb1 dxb2 dxb3]
	#estados y salida
	xb=[x1;x2;x3]
	yb=x3
end

# ╔═╡ 2a98e6e0-a227-11eb-026a-bb5831f344bb
let
	x0 = [0.1 0.1]
	parametros=[1 1 1 1] # una g un poco rara ¿no?
	tspan = (0.0,10.0)
	prob = ODEProblem(derivadas_ejemplo1,x0,tspan,parametros);
	sol = solve(prob,Tsit5());

	# Salida y referencia
	fig1=plot(sol,vars=(0,1), xaxis=L"t",yaxis=L"y(t)",label=L"y(t)")
	#referencia, x1 no se usa para nada pero espera que le pases algún estado
	fig1=plot!(fig1,sol,vars=( (t,x1)->(t,referencia_1(t)[1]) , 0, 1 ), 					label=L"y_r", linestyle=[:dot], linewidth=2)

	#Pinto los estados, si no quiero etiquetasp1=
	fig2=plot(sol, xaxis="t",yaxis=L"x(t)", label=:none)
	#fig2=plot(sol,vars=(0,1), xaxis=L"t",yaxis=L"x(t)",label=L"x_1")
	#fig2=plot!(sol,vars=(0,2), ,label=L"x_2")

	#Señal de control, u no se guarda, hay que recalcularlo
	fig3=plot(sol, vars=((t,x1,x2)->(t,control_1([x1,x2],parametros,t)),0, 1, 2),
		    xaxis=L"t",yaxis=L"u(t)",legend=false)
	l = @layout [a ; b ; c]
	plot(fig1,fig2,fig3, layout=l)
end

# ╔═╡ 0610d2ee-0958-11eb-2490-19e0efcd3996
dyb=derivada(yb,xb,dxb,1)

# ╔═╡ 52eade10-095b-11eb-0717-a9686481d875
ddyb=expand(derivada(yb,xb,dxb,2))

# ╔═╡ 5c70004e-095b-11eb-0068-89e38689515b
md"""El grado relativo es 2, pero hay tres estados... *ups*

El jacobiano es obvio que es singular (2x3)
"""

# ╔═╡ ab5f1de0-0960-11eb-2cf6-414898fd1323
jacobian([yb;dyb],xb)

# ╔═╡ c7d09f80-0a28-11eb-3548-11dbb05bbbbb
md"""Vamos a seguir *como si nada* a ver que pasa...
"""

# ╔═╡ 77f78c80-a997-11eb-0b80-8da3a424a52f
string(solve_for(ddyb~v,u))

# ╔═╡ 706b0500-0a2d-11eb-023e-1de0809fe6cf
referencia_2(t)=0*referencia_1_a(t)

# ╔═╡ c8abde00-095c-11eb-237c-a9fd2bd95b8a
function control_2(x,p,t)
	x1,x2,x3=x

	yr,dyr,ddyr=referencia_2(t) #defindia más abajo

	d0=2
	d1=3

	y=x3
	dy=2*x2

	e0=yr-y
	e1=dyr-dy

	v=d0*e0 + d1*e1 + ddyr

	u= v - 2*sin(x2) - 4*x1*x2
	return u
end

# ╔═╡ 97e9d100-04ca-11eb-07ea-3f5f5addbc45
function derivadas_ejemplo2(x,p,t)
	x1,x2,x3=x
	a=p[1]
	u=control_2(x,p,t);#calculamos el control
	dx1=-a*x1 + exp(2x2)*u
    dx2=2x1*x2 +sin(x2) + u/2
	dx3=2x2
	dx=[dx1 dx2 dx3]
end

# ╔═╡ c8a04540-095c-11eb-0b16-0bbf6ae0b6eb
md"""# ¡Funciona!
Pero si cambiamos a=1 por a=-1 ¿Qué pasaría?"""

# ╔═╡ 062714d0-0a2e-11eb-073f-6bd47b020c54
	constante_a=1

# ╔═╡ 667f8b70-095f-11eb-08a8-0f49a920e3cb
@bind tf Slider(10:100)

# ╔═╡ b117ab50-04cc-11eb-1580-fdafb3c3213b
let
	x0 = [0.1 0.1 0.1]
	parametros=[constante_a]
	tspan = (0.0,tf)
	prob = ODEProblem(derivadas_ejemplo2,x0,tspan,parametros);
	sol = solve(prob,Tsit5());

	# Salida y referencia
	fig1=plot(sol,vars=(0,3), xaxis=L"t",yaxis=L"y(t)",label=L"y(t)")
	fig1=plot!(fig1,sol,vars=( (t,x1)->(t,referencia_2(t)[1]) , 0, 1 ), 					label=L"y_r", linestyle=[:dot], linewidth=2)

	fig2=plot(sol, xaxis="t",yaxis=L"x(t)", label=:none)
	#fig2=plot(sol,vars=(0,1), xaxis=L"t",yaxis=L"x(t)",label=L"x_1")
	#fig2=plot!(sol,vars=(0,2), ,label=L"x_2")

	#Señal de control, u no se guarda, hay que recalcularlo
	fig3=plot(sol, vars=((t,x1,x2, x3)->(t,control_2([x1,x2,x3],parametros,t)),0, 1, 2, 3),
		    xaxis=L"t",yaxis=L"u(t)",legend=false)
	l = @layout [a ; b ; c]
	plot(fig1,fig2,fig3, layout=l)
end

# ╔═╡ 1f8be93e-0962-11eb-14f4-73a59ba80a72
md""" Que raro... para ver que ha pasado veamos la dinámica cero, recordemos las ecuaciones:

$\dot x_1=-ax_1 + e^{2x_2}u$
$\dot x_2=2x_1x_2 +sin(x_2) + \frac{u}{2}$
$\dot x_3=2x_2$
$y=x_3$

Si la salida $y(t)=0$ cero entonces:

$y=x_3=0 \to \dot x_3=0=2x_2 \to x_2=0 \to \dot x_2=0$
$0=0 +sin(0) + \frac{u}{2} \to u=0$

$\dot x_1=-ax_1$

Esto es estable si a es positivo pero *inestable* si a es negativo, eso podría explica por qué x_1 se aleja del origen (dinámica cero inestable)
"""

# ╔═╡ 41656a80-0a7c-11eb-1070-db649ce73728
md"""Para comprobar la teoría ponemos la referenica a cero, cuando la salida se haga cero lo que queda es "la dimámica cero" y la teoría nos dice:

-Dinámica cero asintóticamete estable $\to$ el sistema es **localmente** estable

-Dinámica cero inestable $\to$ el sistema es inestable

**OJO** lo que la teoría no dice es que pueda seguir una trayectoria, en este caso la sigue pero podría ser que no (al seguir una referencia *general* la salida no es cero y si es muy distinta de cero nos alejamos de la zona conocida)
"""

# ╔═╡ Cell order:
# ╟─5c4af6e0-022e-11eb-152f-2153b9b672c1
# ╟─6b1fd530-0a28-11eb-1142-47afd1798c3e
# ╠═85c0d350-022e-11eb-0873-c118ffde014c
# ╠═506a20d0-026f-11eb-2bde-f94724bd27f3
# ╠═f493fdb0-022f-11eb-34a6-4b133c23a71a
# ╠═c803fe10-0231-11eb-1f84-dbe43fd97eb4
# ╠═be8523c0-0239-11eb-0439-592e370c5618
# ╟─85c03440-0274-11eb-3c81-514259334750
# ╟─2a98e6e0-a227-11eb-026a-bb5831f344bb
# ╟─e3e9d0a0-022c-11eb-04ba-97dcb9ba908a
# ╟─6831eb02-023b-11eb-0842-fbe8f8eb8937
# ╠═5fe018a0-023b-11eb-1e69-d59b562edc9b
# ╠═1c0442a0-a993-11eb-1b2d-55eb8b868dc7
# ╠═3ff84c30-0248-11eb-33ba-f99304ee8672
# ╠═6d535e50-a99a-11eb-08c1-4d98301e53b3
# ╠═5525d1a0-a9c7-11eb-1de0-1700c9c20460
# ╠═65893910-a9c7-11eb-2f9f-ebadf375cdbc
# ╠═a0571620-a9c7-11eb-070f-297fa9e406b2
# ╟─39ce8af0-023d-11eb-3178-857433b369ae
# ╟─4c478bf0-023d-11eb-0dce-439d408010fa
# ╟─b7e08f60-023d-11eb-13cb-ad5285e900d8
# ╠═697f24d0-a99a-11eb-26e6-a749fb04cb0f
# ╟─9ecc5ff0-022d-11eb-3715-51831818d04f
# ╠═65ad8360-a99a-11eb-3764-aded1015107e
# ╠═4b753d20-a9c8-11eb-320a-d58ec5e3c431
# ╠═6e31c180-a9c8-11eb-3893-f55415fa4fad
# ╟─561e07e2-0250-11eb-1380-4b5dfa861cd6
# ╠═64db3f40-a99a-11eb-01df-d3ca4846ca33
# ╟─de497f50-0250-11eb-028e-ffe56cf1ba45
# ╠═c2f97c40-a99a-11eb-30e7-359effec0a1d
# ╠═85e1f1f0-a9c9-11eb-2a97-1944f7fcb34b
# ╟─0db8c9c0-0252-11eb-199c-2d69997c6769
# ╠═dd56c0c2-a99a-11eb-0c87-8f8e2c6321a8
# ╟─72e4e740-022d-11eb-3a63-31023b9d464f
# ╠═f1e0ee80-a99a-11eb-0c06-85f8d7086aac
# ╠═03abf0e0-a9ca-11eb-339c-5599c951871c
# ╠═183a6460-a9ca-11eb-340e-872e08f29a9e
# ╠═2d9109e0-a9ca-11eb-08a5-7596eee9e4b2
# ╟─9ec274e0-022d-11eb-085a-7907b8f1eb93
# ╠═aed8d350-0265-11eb-1223-1dd0651632b6
# ╟─0d736c40-0266-11eb-3433-ad75c1617314
# ╠═454272d0-a9cb-11eb-36db-ffb885a95a52
# ╟─f54ff5b2-0266-11eb-1df3-49f3c95729c7
# ╠═2814f040-0267-11eb-3f4d-2933c0c3528a
# ╟─3b125190-0269-11eb-2473-33cbb7451779
# ╠═639dce00-0269-11eb-191e-d16a91f3bcf9
# ╟─e0b21ca0-0a2c-11eb-0e50-efa6c355e7c5
# ╟─f162f3f0-022d-11eb-30d9-55fb6e39c8c0
# ╠═a7745f70-026b-11eb-1abe-d1944165d231
# ╠═a0e6f960-026b-11eb-23c9-410e6d61c4a1
# ╟─922bccde-022d-11eb-17f2-270aae61f923
# ╟─cc0e02e2-04c9-11eb-2e08-b906aaec3959
# ╠═97e9d100-04ca-11eb-07ea-3f5f5addbc45
# ╟─8839fb60-095a-11eb-063e-85fd960ea83a
# ╠═54523da0-0958-11eb-2730-cbef0823f684
# ╠═0610d2ee-0958-11eb-2490-19e0efcd3996
# ╠═52eade10-095b-11eb-0717-a9686481d875
# ╟─5c70004e-095b-11eb-0068-89e38689515b
# ╠═ab5f1de0-0960-11eb-2cf6-414898fd1323
# ╟─c7d09f80-0a28-11eb-3548-11dbb05bbbbb
# ╠═77f78c80-a997-11eb-0b80-8da3a424a52f
# ╠═c8abde00-095c-11eb-237c-a9fd2bd95b8a
# ╠═706b0500-0a2d-11eb-023e-1de0809fe6cf
# ╟─c8a04540-095c-11eb-0b16-0bbf6ae0b6eb
# ╠═062714d0-0a2e-11eb-073f-6bd47b020c54
# ╠═c900ad60-095f-11eb-3d18-9f187abdb9b3
# ╠═667f8b70-095f-11eb-08a8-0f49a920e3cb
# ╟─b117ab50-04cc-11eb-1580-fdafb3c3213b
# ╟─1f8be93e-0962-11eb-14f4-73a59ba80a72
# ╟─41656a80-0a7c-11eb-1070-db649ce73728

max_stars_repo_path: Tema 4.jl

max_stars_repo_name: Dictino/CNL

max_issues_repo_path: Tema 4.jl

max_issues_repo_name: Dictino/CNL

max_forks_repo_path: Tema 4.jl

max_forks_repo_name: Dictino/CNL


text:
using Base.Test

tests = [
	"genxsv",
    "readxsv"]

println("Running tests:")

for my_test in [join([t,".jl"]) for t=tests]
    try
        include(my_test)
        println("\t\033[1m\033[32mPASSED\033[0m: $(my_test)")
    catch e
        anyerrors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(my_test)")
        rethrow(e)
    end
end


text:
include("IOM/src/share/constants.jl")
include("IOM/src/share/PolelikeCoordinate.jl")
include("IOM/src/share/BasicMatrixOperators.jl")
include("IOM/src/share/AdvancedMatrixOperators.jl")

using .PolelikeCoordinate
using NCDatasets
using ArgParse, JSON
using Statistics

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--scheme"
            help = "Scheme: AGA2020, KSC2018, CO2012."
            arg_type = String
            required = true
    end

    return parse_args(s)
end

parsed = parse_commandline()

JSON.print(parsed, 4)

scheme = parsed["scheme"]

Nϵ = 31
Nγ = 31

ϵ_vec  = range(0.5, 1.0, length=Nϵ) / 86400.0
γ_vec = range(0, 0.2, length=Nγ)

Δϵ  = ϵ_vec[2] - ϵ_vec[1]
Δγ  = γ_vec[2] - γ_vec[1]
lat_rng = [-20.0, 20.0]
yr_rng = [1, 20]
σ_w = 10e-5  # 10m/day tolerance
σ_u = 1e-2  # 1cm/s   tolerance
N_layer = 5

output_file = "eps_gam_logpost2D_$(scheme).nc"

logpost_w = zeros(Float64, Nϵ, Nγ)
logpost_u = zeros(Float64, Nϵ, Nγ)
logpost_uw = zeros(Float64, Nϵ, Nγ)

# Load domain and construct operators
domain_file = "CESM_domains/domain.ocn.gx1v6.090206.nc"
gf = PolelikeCoordinate.CurvilinearSphericalGridFile(
    domain_file;
    R  = Re,
    Ω  = Ω,
)

mask_sT = reshape(gf.mask, 1, size(gf.mask)...)
gd_slab = PolelikeCoordinate.genGrid(gf, [0, -50.0])

println("Constructing operators...")
@time amo_slab = AdvancedMatrixOperators(;
    gd     = gd_slab,
    mask_T     = mask_sT,
    deepmask_T = mask_sT,
)

println("Making matrices")

f_sT = 2 * gd_slab.Ω * sin.(gd_slab.ϕ_T)
β_sT = (2 * gd_slab.Ω / gd_slab.R) * cos.(gd_slab.ϕ_T)
f2_sT = f_sT.^2

T_DIVx_T = dropzeros(amo_slab.T_DIVx_U * amo_slab.U_interp_T)
T_DIVy_T = dropzeros(amo_slab.T_DIVy_V * amo_slab.V_interp_T)

files = readdir("hist")

cnt = 0
for f in files

    global cnt

#    if cnt >= 5
#        break
#    end
    m = match(r"paper2021_CTL_POP2\.pop\.h\.(?<yr>\d\d\d\d)-\d\d\.nc", f)

    if m == nothing
        println("Skip $(f)")
        continue
    else
        if ! ( yr_rng[1] <= parse(Int64, m[:yr]) <= yr_rng[2] )
            println("Skip $(f)")
            continue
        end
        println("Loading $(f)")
        cnt += 1
    end

    ds = Dataset("hist/$(f)", "r")

    TAUX   = convert(Array{Float64}, nomissing(ds["TAUX"][:, :, 1], NaN)) / 10.0
    TAUY   = convert(Array{Float64}, nomissing(ds["TAUY"][:, :, 1], NaN)) / 10.0
    u_OGCM = convert(Array{Float64}, nomissing(ds["UVEL"][:, :, 1:N_layer, 1], NaN)) / 100.0
    v_OGCM = convert(Array{Float64}, nomissing(ds["VVEL"][:, :, 1:N_layer, 1], NaN)) / 100.0
    w_OGCM = convert(Array{Float64}, nomissing(ds["WVEL"][:, :, 6, 1], NaN)) / 100.0
    z_w_top = - convert(Array{Float64}, nomissing(ds["z_w_top"][1:N_layer], NaN)) / 100.0
    z_w_bot = - convert(Array{Float64}, nomissing(ds["z_w_bot"][1:N_layer], NaN)) / 100.0
    close(ds)

    valid_idx = (lat_rng[1] .<= gf.yc .<= lat_rng[2]) .& (gf.mask .== 1.0) .& ( isfinite.(w_OGCM) )

    function checkValid(x, idx)
        valid_x = x[idx]
        isfinite_x = isfinite.(valid_x)

        if all(isfinite_x)
            #println("[checkValid] PASS.")
        else
            println("[checkValid] NOT PASS.")
            throw(ErrorException("There are $(length(valid_x) - sum(isfinit_x)) pts that are not finite."))
        end

    end

    checkValid(TAUX, valid_idx)
    checkValid(TAUY, valid_idx)
    checkValid(w_OGCM, valid_idx)

    TAUX = reshape(TAUX, 1, size(TAUX)...)
    TAUY = reshape(TAUY, 1, size(TAUY)...)

    # process
    TAUX_east, TAUY_north = PolelikeCoordinate.project(
        gd_slab,
        TAUX,
        TAUY,
        direction=:Backward,
        grid=:T,
    )

    curlτ_sT = reshape(
        amo_slab.T_CURLx_T * view(TAUY, :) + amo_slab.T_CURLy_T * view(TAUX, :),
        1, gd_slab.Nx, gd_slab.Ny,
    )

    for (i, ϵ) in enumerate(ϵ_vec), (j, γ) in enumerate(γ_vec)

        #println("[$i, $j] Doing ϵ = $(ϵ); γ = $(γ)")

        ϵ_sT  = f_sT * 0 .+ ϵ
        ϵ2_sT  = ϵ_sT.^2

        ϵ2invβ_sT = ϵ2_sT * (gd_slab.R / 2.0 / gd_slab.Ω)
        invD_sT   = (f2_sT + ϵ2_sT).^(-1.0)

        if scheme == "AGA2020"
            VFX_east   = (   f_sT .* TAUY_north                               ) .* invD_sT / ρ_sw
            VFY_north  = ( - f_sT .* TAUX_east + γ * ϵ2invβ_sT .* curlτ_sT    ) .* invD_sT / ρ_sw
        elseif scheme == "CO2012"
            VFX_east   = (   ϵ_sT .* TAUX_east + f_sT .* TAUY_north   ) .* invD_sT / ρ_sw
            VFY_north  = ( - f_sT .* TAUX_east + ϵ_sT .* TAUY_north  ) .* invD_sT / ρ_sw
        elseif scheme == "KSC2018"
            VFX_east   = (   ϵ_sT .* TAUX_east + f_sT .* TAUY_north  ) .* invD_sT / ρ_sw
            VFY_north  = ( - f_sT .* TAUX_east ) .* invD_sT / ρ_sw
        end

        VFX, VFY = PolelikeCoordinate.project(
            gd_slab,
            VFX_east,
            VFY_north,
            direction=:Forward,
            grid=:T,
        )

        w = reshape(( T_DIVx_T * reshape(VFX,:) + T_DIVy_T * reshape(VFY,:) ), gd_slab.Nx, gd_slab.Ny)

        #println("w size = ", size(w))
        #println("w_OGCM size = ", size(w_OGCM))
        Δw = w .- w_OGCM

        tmp = - ( Δw / σ_w ).^2.0 / 2.0
        #println(size(tmp))
        #if ! all( isfinite.(tmp[valid_idx]) )
        #    println(sum(isnan.(tmp[valid_idx])))
        #    throw(ErrorException("Not all is finite"))
        #end

        logpost_w[i, j] += sum(tmp[valid_idx])

        #=
        Dataset("w_check.nc", "c") do ds

            defDim(ds, "Nx", gd_slab.Nx)
            defDim(ds, "Ny", gd_slab.Ny)
            for (varname, vardata, vardim, attrib) in [
                ("w",  w,    ("Nx", "Ny",), Dict()),
                ("w_OGCM",  w_OGCM,    ("Nx", "Ny",), Dict()),
            ]

                var = defVar(ds, varname, Float64, vardim)
                var.attrib["_FillValue"] = 1e20

                for (k, v) in attrib
                    var.attrib[k] = v
                end

                rng = []
                for i in 1:length(vardim)-1
                    push!(rng, Colon())
                end
                push!(rng, 1:size(vardata)[end])
                var[rng...] = vardata

            end

        end
        break
        =#

        # approximate U
        Δz = reshape(z_w_top - z_w_bot, 1, 1, :)
        H  = sum(Δz[1:N_layer])
        VFX_OGCM = reshape(view(sum(Δz .* u_OGCM, dims=3), :, :, 1), 1, gd_slab.Nx, gd_slab.Ny)
        VFY_OGCM = reshape(view(sum(Δz .* v_OGCM, dims=3), :, :, 1), 1, gd_slab.Nx, gd_slab.Ny)

        VFX_east_OGCM, VFY_north_OGCM = PolelikeCoordinate.project(
            gd_slab,
            VFX_OGCM,
            VFY_OGCM,
            direction=:Backward,
            grid=:T,
        )

        Δu  = (VFX_east  - VFX_OGCM) / H
        #Δv  = (VFY_north - VFY_OGCM) / H
        #tmp = - (Δu.^2 + Δv.^2) / (2 * σ_u^2.0)
        tmp = - (Δu.^2) / (2 * σ_u^2.0)
        logpost_u[i, j] += sum( view(tmp, 1, :, :)[valid_idx] )

    end

end

println("Processed $(cnt) files.")

function processLogpost(logpost, Δϵ, Δγ)
    logpost = logpost .- maximum(logpost)
    Δϵ2 = Δϵ * Δγ

    post = exp.(logpost)
    post /= sum(post) * Δϵ2

    return logpost, post
end

logpost_uw = logpost_w + logpost_u

logpost_w, post_w = processLogpost(logpost_w, Δϵ, Δγ)
logpost_u, post_u = processLogpost(logpost_u, Δϵ, Δγ)
logpost_uw, post_uw = processLogpost(logpost_uw, Δϵ, Δγ)

println("Writing output: $(output_file)")
Dataset(output_file, "c") do ds

    defDim(ds, "Neps", Nϵ)
    defDim(ds, "Ngam", Nγ)
    for (varname, vardata, vardim, attrib) in [
        ("logpost_w",  logpost_w,  ("Neps", "Ngam",), Dict()),
        ("post_w",     post_w,     ("Neps", "Ngam"), Dict()),
        ("logpost_u",  logpost_u,  ("Neps", "Ngam",), Dict()),
        ("post_u",     post_u,     ("Neps", "Ngam"), Dict()),
        ("logpost_uw", logpost_uw, ("Neps", "Ngam",), Dict()),
        ("post_uw",    post_uw,    ("Neps", "Ngam"), Dict()),
        ("epsilon",        ϵ_vec,    ("Neps",), Dict()),
        ("gamma",  γ_vec,    ("Ngam",), Dict()),
    ]

        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

# 1. project TAUX and TAUY to true east and north
# 2. compute curl
# 3. compute horizontal velocity
# 4. compute divergence and thus w_50m

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module HSA

using CEnum

include("libhsa_common.jl")
include("libhsa_enums.jl")
include("libhsa_types.jl")
include("libhsa_api.jl")

end


text:
import DandelionWebSockets: getrequestheaders, HTTPHandshakeLogic
import HTTP
import JSON

include("../src/PlutoBindClient.jl")

function get_latest_events()
        r = HTTP.get("https://api.github.com/events", ["Accept" => "application/vnd.github.v3+json"])
        JSON.parse(String(r.body))
end

function send_events(handler::PlutoClient)
        events = get_latest_events()
        @info "Got $(length(events)) events"
        for event in events
                set_bond_value(handler, "bonds/ghevents", event)
                sleep(0.1) # throttle a bit
        end
end

notebook_id, cookie = ARGS

getrequestheaders(h::HTTPHandshakeLogic) = [
        "Sec-WebSocket-Version" => "13",
        "Upgrade" => "websocket",
        "Connection" => "Upgrade",
        "Sec-WebSocket-Key" => h.key,
        "Cookie" => "secret=$cookie"]

client = PlutoClient(notebook_id=notebook_id)
ok = connect(client, "127.0.0.1:1235")
if !ok
        @error "Failed to connect to the notebook"
        exit(1)
end

while true
        send_events(client)
        sleep(.4)
end

take!(client.stop_channel)


text:
@testset "nqueens.jl" begin
    @testset "fill_with_generator!" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        n=8

        generator = SeaPearl.NQueensGenerator(n)

        for _ in 1:3
            SeaPearl.fill_with_generator!(model, generator)

            @test length(keys(model.variables)) == n
            @test length(model.constraints) == 3
            empty!(model)
        end
    end
end


text:
f(x) = x.^2
println("Square it: $(f(2))")

##Plots.jl is extremely slow for the first plot.  Much faster after
using Plots
gr() # Choose a backend

##Alternatively, use Gadfly
#using Gadfly
p = plot(1:7, f(1:7)) # This will plot to the plot pane
display(p) #In Juno and the REPL, you need to display a plot


text:
using  Gadfly

set_default_plot_size(6inch, 3inch)

x= 1:10
s = [-1,-1,1,1,-1,-1,1,1,-1,-1]
pa = plot(x=x, y=x.^2, Geom.hair, Geom.point)
pb = plot(x=s.*(x.^2), y=x, Geom.hair(orientation=:horizontal), Geom.point, color=string.(s), Theme(key_position=:none))
hstack(pa, pb)