text:
"""
Abstract supertype for layerwise training specification.
May be specifications for a normal RBM layer (see `TrainLayer`) or
multiple combined specifications for a partitioned layer
(see `TrainPartitionedLayer`).
"""
abstract type AbstractTrainLayer end

mutable struct TrainLayer <: AbstractTrainLayer
   epochs::Int
   usedefaultepochs::Bool
   learningrate::Float64
   learningrates::Vector{Float64}
   usedefaultlearningrate::Bool
   sdlearningrate::Float64
   sdlearningrates::Vector{Float64}
   monitoring::Function
   rbmtype::DataType
   nhidden::Int
   nvisible::Int
   batchsize::Int
   usedefaultbatchsize::Bool
   pcd::Bool
   cdsteps::Int
   startrbm::AbstractRBM
   optimizer::AbstractOptimizer
   optimizers::Vector{AbstractOptimizer}
   fisherscoring::Int
   inversedispersion::Array{Float64,1}
   dispersionquantile::Float64
   estimatedispersion::String
   lambda::Float64
   zeroinflation::Bool
end

"""
Specify parameters for training one RBM-layer in a DBM.

# Optional keyword arguments:
* The optional keyword arguments `rbmtype`, `nhidden`, `epochs`,
  `learningrate`/`learningrates`, `sdlearningrate`/`sdlearningrates`,
  `batchsize`, `pcd`, `cdsteps`, `startrbm` and `optimizer`/`optimizers`
  are passed to `fitrbm`. For a detailed description, see there.
  If a negative value is specified for `learningrate` or `epochs`, this indicates
  that a corresponding default value should be used
  (parameter defined by call to `stackrbms`).
* `monitoring`: also like in `fitrbm`, but may take a `DataDict` as third argument
   (see function `stackrbms` and its argument `monitoringdata`).
* `nvisible`: Number of visible units in the RBM. Only relevant for partitioning.
   This parameter is derived as much as possible by `stackrbms`.
   For `MultimodalDBM`s with a partitioned first layer, it is necessary to specify
   the number of visible nodes for all but at most one partition in the input layer.
"""
function TrainLayer(;
      epochs::Int = -1,
      learningrate::Float64 = -Inf,
      learningrates::Vector{Float64} = Vector{Float64}(),
      sdlearningrate::Float64 = 0.0,
      sdlearningrates::Vector{Float64} = Vector{Float64}(),
      monitoring = nomonitoring,
      rbmtype::DataType = BernoulliRBM,
      nhidden::Int = -1,
      nvisible::Int = -1,
      batchsize::Int = -1,
      pcd::Bool = true,
      cdsteps::Int = 1,
      startrbm::AbstractRBM = NoRBM(),
      optimizer::AbstractOptimizer = NoOptimizer(),
      optimizers::Vector{<:AbstractOptimizer} = Vector{AbstractOptimizer}(),
      fisherscoring::Int = 0,
      inversedispersion::Array{Float64,1} = (ones(1500) .* 1.0),   #TODO: Adapt to proper number of dimensions
      dispersionquantile::Float64 = 0.5,
      estimatedispersion::String = "gene",
      lambda::Float64 = 10.0,
      zeroinflation::Bool = false)

   usedefaultepochs = (epochs < 0)
   usedefaultlearningrate = (learningrate < 0 && isempty(learningrates))
   usedefaultbatchsize = (batchsize <= 0)
   TrainLayer(
         (usedefaultepochs ? 10 : epochs),
         usedefaultepochs,
         (usedefaultlearningrate ? 0.005 : learningrate),
         learningrates,
         usedefaultlearningrate,
         sdlearningrate, sdlearningrates,
         monitoring, rbmtype, nhidden, nvisible,
         (usedefaultbatchsize ? 1 : batchsize),
         usedefaultbatchsize,
         pcd, cdsteps,
         startrbm,
         optimizer, optimizers, fisherscoring, inversedispersion, dispersionquantile, estimatedispersion, lambda, zeroinflation)
end

"""
Encapsulates a vector of `TrainLayer` objects for training a partitioned layer.
"""
struct TrainPartitionedLayer <: AbstractTrainLayer
   parts::Vector{TrainLayer}
end

const AbstractTrainLayers = Vector{<:AbstractTrainLayer}

"""
A dictionary containing names of data sets as keys and the data sets (matrices
with samples in rows) as values.
"""
const DataDict = Dict{String,Array{Float64,2}}

"""
Converts a vector to a vector of the most specific type that all
elements share as common supertype.
"""
function converttomostspecifictype(v::Vector)
   Vector{mostspecifictype(v)}(v)
end

"""
    mostspecifictype(v)
Returns the most specific supertype for all elements in the vector `v`.
"""
function mostspecifictype(v::Vector)
   mapreduce(typeof, typejoin, v)
end

"""
    propagateforward(rbm, datadict, factor)
Returns a new `DataDict` containing the same labels as the given `datadict` but
as mapped values it contains the hidden potential in the `rbm` of the original
datasets. The factor is applied for calculating the hidden potential and is 1.0
by default.
"""
function propagateforward(rbm::AbstractRBM, datadict::DataDict, factor::Float64 = 1.0)
   DataDict(p[1] => hiddenpotential(rbm, p[2], factor) for p = pairs(datadict))
end

function partitioneddata(datadict::DataDict, visrange::UnitRange{Int})
   DataDict(p[1] => p[2][:, visrange] for p = pairs(datadict))
end

function setdefaultsforunspecified!(trainlayer::TrainLayer,
         learningrate::Float64, epochs::Int, batchsize::Int,
         optimizer::AbstractOptimizer)

   if trainlayer.usedefaultlearningrate
      trainlayer.learningrate = learningrate
   end
   if trainlayer.usedefaultepochs
      trainlayer.epochs = epochs
   end
   if trainlayer.usedefaultbatchsize
      trainlayer.batchsize = batchsize
   end
   if trainlayer.optimizer === NoOptimizer()
      trainlayer.optimizer = optimizer
   end

   if isempty(trainlayer.learningrates)
      # learning rates for each epoch have neither been explicitly specified
      # --> set same value for each epoch
      trainlayer.learningrates =
            fill(trainlayer.learningrate, trainlayer.epochs)
   end
   if isempty(trainlayer.sdlearningrates)
      # same for sdlearningrates
      trainlayer.sdlearningrates =
            fill(trainlayer.sdlearningrate, trainlayer.epochs)
   end
   trainlayer
end

function setdefaultsforunspecified!(trainpartitionedlayer::TrainPartitionedLayer,
      learningrate::Float64, epochs::Int, batchsize::Int,
      optimizer::AbstractOptimizer)

   setdefaultsforunspecified!(trainpartitionedlayer.parts,
         learningrate, epochs, batchsize, optimizer)
end

function setdefaultsforunspecified!(trainlayers::Vector{T},
      learningrate::Float64, epochs::Int, batchsize::Int,
      optimizer::AbstractOptimizer
      ) where {T <: AbstractTrainLayer}

   for trainlayer in trainlayers
      setdefaultsforunspecified!(trainlayer,
            learningrate, epochs, batchsize, optimizer)
   end
   trainlayers
end

"""
    stackrbms(x; ...)
Performs greedy layerwise training for Deep Belief Networks or greedy layerwise
pretraining for Deep Boltzmann Machines and returns the trained model.

# Optional keyword arguments (ordered by importance):
* `predbm`: boolean indicating that the greedy layerwise training is
   pre-training for a DBM.
   If its value is false (default), a DBN is trained.
* `nhiddens`: vector containing the number of nodes of the i'th hidden layer in
   the i'th entry
* `epochs`: number of training epochs
* `learningrate`: learningrate, default 0.005
* `batchsize`: size of minibatches, defaults to 1
* `trainlayers`: a vector of `TrainLayer` objects. With this argument it is possible
   to specify the training parameters for each layer/RBM individually.
   If the number of training epochs and the learning rate are not specified
   explicitly for a layer, the values of `epochs` and `learningrate` are used.
   For more information see help of `TrainLayer`.
* `monitoringdata`: a data dictionary (see type `DataDict`)
   The data is propagated forward through the
   network to monitor higher levels.
   If a non-empty dictionary is given, the monitoring functions in the
   `trainlayers`-arguments must accept a `DataDict` as third argument.
* `samplehidden`: boolean indicating that consequent layers are to be trained
   with sampled values instead of the deterministic potential,
   which is the default.
"""
function stackrbms(x::Array{Float64,2};
      nhiddens::Vector{Int} = Vector{Int}(),
      epochs::Int = 10,
      predbm::Bool = false,
      samplehidden::Bool = false,
      learningrate::Float64 = 0.005,
      batchsize::Int = 1,
      optimizer::AbstractOptimizer = NoOptimizer(),
      trainlayers::AbstractTrainLayers = Vector{TrainLayer}(),
      monitoringdata::DataDict = DataDict())

   stackrbms_checkmonitoringdata(x, monitoringdata)

   trainlayers = stackrbms_preparetrainlayers(trainlayers, x, epochs,
         learningrate, nhiddens, batchsize, optimizer)

   nrbms = length(trainlayers)
   dbmn = Vector{AbstractRBM}(undef, nrbms)

   upfactor = downfactor = 1.0
   if predbm
      upfactor = 2.0
   end

   dbmn[1] = stackrbms_trainlayer(x, trainlayers[1];
         monitoringdata = monitoringdata,
         upfactor = upfactor, downfactor = downfactor)

   hiddenval = x
   for i = 2:nrbms
      hiddenval = hiddenpotential(dbmn[i-1], hiddenval, upfactor)

      if samplehidden
         bernoulli!(hiddenval)
      end

      if !isempty(monitoringdata)
         monitoringdata = propagateforward(dbmn[i-1], monitoringdata, upfactor)
      end

      if predbm
         upfactor = downfactor = 2.0
         if i == nrbms
            upfactor = 1.0
         end
      else
         upfactor = downfactor = 1.0
      end

      dbmn[i] = stackrbms_trainlayer(hiddenval, trainlayers[i];
            monitoringdata = monitoringdata,
            upfactor = upfactor, downfactor = downfactor)
   end

   dbmn = converttomostspecifictype(dbmn)
   dbmn
end

function stackrbms_checkmonitoringdata(x::Matrix{Float64}, monitoringdata::DataDict)
   if !isempty(monitoringdata)
      if any(map(m -> size(m, 2), values(monitoringdata)) .!= size(x, 2))
         error("Matrices in the dictionary `monitoringdata` must have the same number of columns as the input matrix.")
      end
   end
end

""" Prepares the layerwise training specifications for `stackrbms` """
function stackrbms_preparetrainlayers(
      trainlayers::AbstractTrainLayers,
      x::Matrix{Float64},
      epochs::Int,
      learningrate::Float64,
      nhiddens::Vector{Int},
      batchsize::Int,
      optimizer::AbstractOptimizer)

   if isempty(trainlayers)
      # construct default "trainlayers"
      if isempty(nhiddens)
         nhiddens = fill(size(x,2), 2) # default value for nhiddens
      end
      trainlayers = map(n -> TrainLayer(nhidden = n), nhiddens)
      setdefaultsforunspecified!(trainlayers,
            learningrate, epochs, batchsize, optimizer)
      return trainlayers
   end

   # We are here, Argument "trainlayers" has been specified
   # --> check for correct specification and derive/set some parameters

   # avoid modification of given object
   trainlayers = deepcopy(trainlayers)

   if !isempty(nhiddens)
      warn("Argument `nhiddens` not used.")
   end

   setdefaultsforunspecified!(trainlayers,
         learningrate, epochs, batchsize, optimizer)

   function derive_nvisibles!(layer::AbstractTrainLayer,
         prevlayer::AbstractTrainLayer)
      # do nothing
   end

   function derive_nvisibles!(layer::TrainPartitionedLayer,
         prevlayer::TrainPartitionedLayer)

      if length(layer.parts) == length(prevlayer.parts)
         for j in eachindex(layer.parts)
            if layer.parts[j].nvisible < 0
               layer.parts[j].nvisible = prevlayer.parts[j].nhidden
            end
         end
      end
      if any(map(t -> t.nvisible < 0, layer.parts))
         error("Parameter `nvisible` could not be derived.")
      end
   end

   # derive number of visible nodes for higher layers
   for i = 2:length(trainlayers)
      derive_nvisibles!(trainlayers[i], trainlayers[i-1])
   end

   # check number of visible nodes for first layer
   if isa(trainlayers[1], TrainPartitionedLayer)
      partswithoutnvisible = findall(t -> t.nvisible <= 0, trainlayers[1].parts)
      if length(partswithoutnvisible) > 1
         error("The number of visible nodes is unspecified for more " *
               "than one segment in the first layer. ")
      end

      nvariables = size(x, 2)
      nvisiblestotal = mapreduce(t -> max(t.nvisible, 0), +, trainlayers[1].parts)

      if length(partswithoutnvisible) == 1
         # if only one is unspecified, derive it
         partwithunspecifiednvisible = trainlayers[1].parts[partswithoutnvisible[1]]
         partwithunspecifiednvisible.nvisible = nvariables - nvisiblestotal
      elseif nvisiblestotal != nvariables
         error("Number of visible nodes for first layer " *
               "($nvisiblestotal) is not of same length as " *
               "number of columns in `x` ($(size(x,2))).")
      end
   end

   trainlayers
end

""" Trains a layer without partitioning for `stackrbms`. """
function stackrbms_trainlayer(x::Matrix{Float64},
      trainlayer::TrainLayer;
      monitoringdata::DataDict = DataDict(),
      upfactor::Float64 = 1.0, downfactor::Float64 = 1.0)

   if isempty(monitoringdata)
      monitoring = trainlayer.monitoring
   elseif trainlayer.monitoring != nomonitoring
      monitoring = (rbm, epoch) ->
            trainlayer.monitoring(rbm, epoch, monitoringdata)
   end

   fitrbm(x;
         upfactor = upfactor, downfactor = downfactor,
         nhidden = trainlayer.nhidden,
         epochs = trainlayer.epochs,
         rbmtype = trainlayer.rbmtype,
         learningrates = trainlayer.learningrates,
         batchsize = trainlayer.batchsize,
         cdsteps = trainlayer.cdsteps,
         pcd = trainlayer.pcd,
         sdlearningrate = trainlayer.sdlearningrate,
         sdlearningrates = trainlayer.sdlearningrates,
         monitoring = monitoring,
         startrbm = trainlayer.startrbm,
         optimizer = trainlayer.optimizer,
         optimizers = trainlayer.optimizers,
         fisherscoring = trainlayer.fisherscoring,
         inversedispersion = trainlayer.inversedispersion,
         dispersionquantile = trainlayer.dispersionquantile,
         estimatedispersion = trainlayer.estimatedispersion,
         zeroinflation = trainlayer.zeroinflation,
         lambda = trainlayer.lambda)
end

""" Trains a partitioned layer for `stackrbms`. """
function stackrbms_trainlayer(x::Matrix{Float64},
      trainpartitionedlayer::TrainPartitionedLayer;
      monitoringdata::DataDict = DataDict(),
      upfactor::Float64 = 1.0, downfactor::Float64 = 1.0)

   visranges = ranges(map(t -> t.nvisible, trainpartitionedlayer.parts))

   # prepare the arguments before the for-loop for each call
   # to avoid unnecessary copying
   trainingargs = map(
         i -> (
               x[:, visranges[i]],
               trainpartitionedlayer.parts[i],
               partitioneddata(monitoringdata, visranges[i])
         ),
         eachindex(trainpartitionedlayer.parts))

   rbms = @distributed (vcat) for arg in trainingargs
      stackrbms_trainlayer(arg[1], arg[2];
            monitoringdata = arg[3],
            upfactor = upfactor, downfactor = downfactor)
   end

   commontype = mostspecifictype(rbms)
   PartitionedRBM{commontype}(Vector{commontype}(rbms))
end


text:
function d1hat(n::Int64, dx::Float64 )
    operator = zeros(n,n) # initialize
    for k in 1:n
        try operator[k, k+1] = 1. catch; Nothing end
        try operator[k, k-1] = -1. catch; Nothing end
    end

    return operator./(2*dx)
end

function d2hat(n::Int64, dx::Float64 )
    operator = zeros(n,n)
    for k in 1:n
        operator[k,k] = -2.
        try operator[k, k+1] = 1. catch; Nothing end
        try operator[k, k-1] = 1. catch; Nothing end
    end

    return operator./(dx^2)
end

ext: jl


text:
@testset "inset.jl" begin
    @testset "InSet()" begin
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        x = SeaPearl.IntVar(0, 3, "x", trailer)

        constraint = SeaPearl.InSet(x, a, trailer)

        @test constraint in a.onDomainChange
        @test constraint in x.onDomainChange
        @test constraint.active.value
    end

    @testset "propagate!(::InSet)" begin
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        x = SeaPearl.IntVar(0, 3, "x", trailer)
        constraint = SeaPearl.InSet(x, a, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test length(x.domain) == 2
        @test !(0 in x.domain)
        @test !(1 in x.domain)
        @test 2 in x.domain
        @test constraint.active.value
        @test prunedDomains == SeaPearl.CPModification("x" => [0, 1])

        SeaPearl.assign!(x, 3)

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test SeaPearl.required_values(a.domain) == Set{Int}([3])
        @test SeaPearl.possible_not_required_values(a.domain) == Set{Int}([2, 4, 5])
        @test !constraint.active.value

        # Infeasibility
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        x = SeaPearl.IntVar(0, 1, "x", trailer)
        constraint = SeaPearl.InSet(x, a, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test !SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        #Infeasibility with excluding
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        x = SeaPearl.IntVar(2, 3, "x", trailer)
        constraint = SeaPearl.InSet(x, a, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        SeaPearl.exclude!(a.domain, 2)
        SeaPearl.exclude!(a.domain, 3)

        @test !SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        # Deactivation
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        x = SeaPearl.IntVar(2, 3, "x", trailer)
        constraint = SeaPearl.InSet(x, a, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        SeaPearl.require!(a.domain, 2)
        SeaPearl.require!(a.domain, 3)

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
        @test !constraint.active.value
    end
end


text:
export rational_reconstruction, farey_lift, div, valence, leading_coefficient,
       trailing_coefficient, constant_coefficient, factor_mod_pk,
       factor_mod_pk_init, hensel_lift, rres, rresx,
       coefficients

function PolynomialRing(R::Ring; cached::Bool = false)
  return PolynomialRing(R, "x", cached = cached)
end

function PolynomialRing(R::FlintRationalField, a::Symbol; cached::Bool = true)
  Qx = FmpqPolyRing(R, a, cached)
  return Qx, gen(Qx)
end

function PolynomialRing(R::FlintIntegerRing, a::Symbol; cached::Bool = true)
  Zx = FmpzPolyRing(a, cached)
  return Zx, gen(Zx)
end

function fmpz(a::Generic.Res{Nemo.fmpz})
  return a.data
end

function fmpz(a::Nemo.nmod)
  return fmpz(a.data)
end

function lift(::FlintIntegerRing, a::Generic.Res{Nemo.fmpz})
  return a.data
end

function (::FlintIntegerRing)(a::Generic.Res{Nemo.fmpz})
  return a.data
end

function lift(::FlintIntegerRing, a::Nemo.nmod)
  return fmpz(a.data)
end

function (::FlintIntegerRing)(a::Nemo.nmod)
  return fmpz(a.data)
end

function div(f::PolyElem, g::PolyElem)
  q, r = divrem(f,g)
  return q
end

function rem(f::PolyElem, g::PolyElem)
  return mod(f, g)
end

function ismonic(a::PolyElem)
  return isone(lead(a))
end

@doc Markdown.doc"""
  valence(f::PolyElem) -> RingElem

 The last non-zero coefficient of f
"""
function valence(f::PolyElem)
  for i=0:degree(f)
    c = coeff(f, i)
    if !iszero(c)
      return c
    end
  end
  return c
end

@doc Markdown.doc"""
  leading_coefficient(f::PolyElem) -> RingElem

 The last leading coefficient of f.
"""
leading_coefficient(f::PolyElem) = lead(f)

@doc Markdown.doc"""
  trailing_coefficient(f::PolyElem) -> RingElem
  constant_coefficient(f::PolyElem) -> RingElem

 The constant coefficient of f.
"""
function trailing_coefficient(f::PolyElem)
  if iszero(f)
    return base_ring(f)(0)
  end
  return coeff(f, 0)
end

@doc Markdown.doc"""
    induce_rational_reconstruction(a::fmpz_poly, M::fmpz) -> fmpq_poly
Apply {{{rational_reconstruction}}} to each coefficient of $a$, resulting
in either a fail (return (false, s.th.)) or (true, g) for some rational
polynomial $g$ s.th. $g \equiv a \bmod M$.
"""
function induce_rational_reconstruction(a::fmpz_poly, M::fmpz)
  b = PolynomialRing(FlintQQ, parent(a).S, cached = false)[1]()
  for i=0:degree(a)
    fl, x,y = rational_reconstruction(coeff(a, i), M)
    if fl
      setcoeff!(b, i, x//y)
    else
      return false, b
    end
  end
  return true, b
end

const constant_coefficient = trailing_coefficient

function resultant(f::fmpz_poly, g::fmpz_poly, d::fmpz, nb::Int)
  z = fmpz()
  ccall((:fmpz_poly_resultant_modular_div, :libflint), Nothing,
     (Ref{fmpz}, Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz}, Int),
     z, f, g, d, nb)
  return z
end

################################################################################
#
#  Random polynomial
#
################################################################################

@doc Markdown.doc"""
    Base.rand(Rt::PolyRing{T}, n::Int) where T <: ResElem{S} where S <: Union{fmpz, Integer} -> PolElem{T}
Find a random polynomial of degree(n)
"""
function Base.rand(Rt::PolyRing{T}, n::Int) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  f = Rt()
  R = base_ring(Rt)
  for i=0:n
    setcoeff!(f, i, rand(R))
  end
  return f
end

################################################################################
#
# Hensel
#
################################################################################
#cannot use Ref here, has to be Ptr as the underlying stuff is
#not Julia allocated (but flint)
mutable struct fmpz_poly_raw  ## fmpz_poly without parent like in c
  coeffs::Ptr{Nothing}
  alloc::Int
  length::Int

  function fmpz_poly_raw()
    error("should not get called")
    z = new()
    ccall((:fmpz_poly_init, :libflint), Nothing, (Ref{fmpz_poly},), z)
    finalizer(_fmpz_poly_raw_clear_fn, z)
    return z
  end
end

function _fmpz_poly_raw_clear_fn(a::fmpz_poly)
  ccall((:fmpz_poly_clear, :libflint), Nothing, (Ref{fmpz_poly},), a)
end

mutable struct fmpz_poly_factor
  c::Int   # really an fmpz  - but there is no fmpz_raw to be flint compatible
  poly::Ptr{fmpz_poly_raw}
  exp::Ptr{Int}
  _num::Int
  _alloc::Int

  function fmpz_poly_factor()
    z = new()
    ccall((:fmpz_poly_factor_init, :libflint), Nothing,
            (Ref{fmpz_poly_factor}, ), z)
    finalizer(_fmpz_poly_factor_clear_fn, z)
    return z
  end
end

function _fmpz_poly_factor_clear_fn(a::fmpz_poly_factor)
  ccall((:fmpz_poly_factor_clear, :libflint), Nothing,
          (Ref{fmpz_poly_factor}, ), a)
end

function factor_to_dict(a::fmpz_poly_factor)
  res = Dict{fmpz_poly,Int}()
  Zx,x = PolynomialRing(FlintZZ, "x", cached = false)
  for i in 1:a._num
    f = Zx()
    ccall((:fmpz_poly_set, :libflint), Nothing, (Ref{fmpz_poly}, Ref{fmpz_poly_raw}), f, a.poly+(i-1)*sizeof(fmpz_poly_raw))
    res[f] = unsafe_load(a.exp, i)
  end
  return res
end

function show(io::IO, a::fmpz_poly_factor)
  ccall((:fmpz_poly_factor_print, :libflint), Nothing, (Ref{fmpz_poly_factor}, ), a)
end

mutable struct HenselCtx
  f::fmpz_poly
  p::UInt

  LF :: fmpz_poly_factor
  link::Ptr{Int}
  v::Ptr{fmpz_poly_raw}
  w::Ptr{fmpz_poly_raw}
  N::UInt
  prev::UInt
  r::Int  #for the cleanup
  lf:: Nemo.nmod_poly_factor

  function HenselCtx(f::fmpz_poly, p::fmpz)
    a = new()
    a.f = f
    a.p = UInt(p)
    Zx,x = PolynomialRing(FlintZZ, "x", cached=false)
    Rx,x = PolynomialRing(GF(UInt(p), cached=false), "x", cached=false)
    a.lf = Nemo.nmod_poly_factor(UInt(p))
    ccall((:nmod_poly_factor, :libflint), UInt,
          (Ref{Nemo.nmod_poly_factor}, Ref{gfp_poly}), (a.lf), Rx(f))
    r = a.lf.num
    a.r = r
    a.LF = fmpz_poly_factor()
    @assert r > 1  #flint restriction
    a.v = ccall((:flint_malloc, :libflint), Ptr{fmpz_poly_raw}, (Int, ), (2*r-2)*sizeof(fmpz_poly_raw))
    a.w = ccall((:flint_malloc, :libflint), Ptr{fmpz_poly_raw}, (Int, ), (2*r-2)*sizeof(fmpz_poly_raw))
    for i=1:(2*r-2)
      ccall((:fmpz_poly_init, :libflint), Nothing, (Ptr{fmpz_poly_raw}, ), a.v+(i-1)*sizeof(fmpz_poly_raw))
      ccall((:fmpz_poly_init, :libflint), Nothing, (Ptr{fmpz_poly_raw}, ), a.w+(i-1)*sizeof(fmpz_poly_raw))
    end
    a.link = ccall((:flint_calloc, :libflint), Ptr{Int}, (Int, Int), 2*r-2, sizeof(Int))
    a.N = 0
    a.prev = 0
    finalizer(HenselCtx_free, a)
    return a
  end

  function free_fmpz_poly_array(p::Ref{fmpz_poly_raw}, r::Int)
    for i=1:(2*r-2)
      ccall((:fmpz_poly_clear, :libflint), Nothing, (Ref{fmpz_poly_raw}, ), p+(i-1)*sizeof(fmpz_poly_raw))
    end
    ccall((:flint_free, :libflint), Nothing, (Ref{fmpz_poly_raw}, ), p)
  end
  function free_int_array(a::Ref{Int})
    ccall((:flint_free, :libflint), Nothing, (Ref{Int}, ), a)
  end
  function HenselCtx_free(a::HenselCtx)
    free_fmpz_poly_array(a.v, a.r)
    free_fmpz_poly_array(a.w, a.r)
    free_int_array(a.link)
  end
end

function show(io::IO, a::HenselCtx)
  println("factorisation of $(a.f) modulo $(a.p)^$(a.N)")
  if a.N > 0
    d = factor_to_dict(a.LF)
    println("currently: $d")
  end
end

function start_lift(a::HenselCtx, N::Int)
  a.prev = ccall((:_fmpz_poly_hensel_start_lift, :libflint), UInt,
       (Ref{fmpz_poly_factor}, Ref{Int}, Ref{fmpz_poly_raw}, Ref{fmpz_poly_raw}, Ref{fmpz_poly}, Ref{Nemo.nmod_poly_factor}, Int),
       a.LF, a.link, a.v, a.w, a.f, a.lf, N)
  a.N = N
end

function continue_lift(a::HenselCtx, N::Int)
  a.prev = ccall((:_fmpz_poly_hensel_continue_lift, :libflint), Int,
       (Ref{fmpz_poly_factor}, Ref{Int}, Ref{fmpz_poly_raw}, Ref{fmpz_poly_raw}, Ref{fmpz_poly}, UInt, UInt, Int, Ref{fmpz}),
       a.LF, a.link, a.v, a.w, a.f, a.prev, a.N, N, fmpz(a.p))
  a.N = N
end

@doc Markdown.doc"""
  factor_mod_pk(f::fmpz_poly, p::Int, k::Int) -> Dict{fmpz_poly, Int}

 For f that is square-free modulo p, return the factorisation modulo p^k.
"""
function factor_mod_pk(f::fmpz_poly, p::Int, k::Int)
  H = HenselCtx(f, fmpz(p))
  start_lift(H, k)
  return factor_to_dict(H.LF)
end

@doc Markdown.doc"""
  factor_mod_pk_init(f::fmpz_poly, p::Int) -> HenselCtx

 For f that is square-free modulo p, return a structure that allows to compute
 the factorisation modulo p^k for any k
"""
function factor_mod_pk_init(f::fmpz_poly, p::Int)
  H = HenselCtx(f, fmpz(p))
  return H
end

@doc Markdown.doc"""
  factor_mod_pk(H::HenselCtx, k::Int) -> RingElem

 Using the result of factor_mod_pk_init, return a factorisation modulo p^k
"""
function factor_mod_pk(H::HenselCtx, k::Int)
  @assert k>= H.N
  if H.N == 0
    start_lift(H, k)
  else
    continue_lift(H, k)
  end
  return factor_to_dict(H.LF)
end

#I think, experimentally, that p = Q^i, p1 = Q^j and j<= i is the condition to make it tick.
function hensel_lift!(G::fmpz_poly, H::fmpz_poly, A::fmpz_poly, B::fmpz_poly, f::fmpz_poly, g::fmpz_poly, h::fmpz_poly, a::fmpz_poly, b::fmpz_poly, p::fmpz, p1::fmpz)
  ccall((:fmpz_poly_hensel_lift, :libflint), Nothing, (Ref{fmpz_poly}, Ref{fmpz_poly},  Ref{fmpz_poly},  Ref{fmpz_poly},  Ref{fmpz_poly},  Ref{fmpz_poly},  Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz_poly}, Ref{fmpz}, Ref{fmpz}), G, H, A, B, f, g, h, a, b, p, p1)
end

@doc Markdown.doc"""
  hensel_lift(f::fmpz_poly, g::fmpz_poly, h::fmpz_poly, p::fmpz, k::Int) -> (fmpz_poly, fmpz_poly)

 Given f = gh modulo p for g, h coprime modulo p, compute G, H s.th. f = GH mod p^k and
 G = g mod p, H = h mod p.
"""
function hensel_lift(f::fmpz_poly, g::fmpz_poly, h::fmpz_poly, p::fmpz, k::Int)
  Rx, x = PolynomialRing(ResidueField(FlintZZ, p, cached=false), cached=false)
  fl, a, b = gcdx(Rx(g), Rx(h))
  @assert isone(fl)
  @assert k>= 2
  ## if one of the cofactors is zero, this crashes.
  ## this can only happen if one of the factors is one. In this case, the other
  ## is essentially f and f would be a legal answer. Probably reduced mod p^k
  ## with all non-negative coefficients
  ## for now:
  @assert !iszero(a) && !iszero(b)
  a = lift(parent(g), a)
  b = lift(parent(g), b)
  G = parent(g)()
  H = parent(g)()
  A = parent(g)()
  B = parent(g)()
  g = deepcopy(g)
  h = deepcopy(h)

  # the idea is to have a good chain of approximations, ie.
  # to reach p^10, one should do p, p^2, p^3, p^5, p^10
  # rather than p, p^2, p^4, p^8, p^10
  # the chain has the same length, but smaller entries.
  l = Int[k]
  while k>1
    k = div(k+1, 2)
    push!(l, k)
  end
  ll = Int[]
  for i=length(l)-1:-1:1
    push!(ll, l[i] - l[i+1])
  end
  P = p
  for i in ll
    p1 = p^i
    hensel_lift!(G, H, A, B, f, g, h, a, b, P, p1)
    G, g = g, G
    H, h = h, H
    A, a = a, A
    B, b = b, B
    P *= p1
  end
  return g, h
end

@doc Markdown.doc"""
  hensel_lift(f::fmpz_poly, g::fmpz_poly, p::fmpz, k::Int) -> (fmpz_poly, fmpz_poly)

 Given f and g such that g is a divisor of f mod p and g and f/g are coprime, compute a hensel lift of g modulo p^k.
"""
function hensel_lift(f::fmpz_poly, g::fmpz_poly, p::fmpz, k::Int)
  Rx, x = PolynomialRing(ResidueField(FlintZZ, p, cached=false), cached=false)
  h = lift(parent(f), div(Rx(f), Rx(g)))
  return hensel_lift(f, g, h, p, k)[1]
end

function fmpq_poly_to_nmod_poly_raw!(r::nmod_poly, a::fmpq_poly)
  ccall((:fmpq_poly_get_nmod_poly, :libflint), Nothing, (Ref{nmod_poly}, Ref{fmpq_poly}), r, a)
end

function fmpq_poly_to_gfp_poly_raw!(r::gfp_poly, a::fmpq_poly)
  ccall((:fmpq_poly_get_nmod_poly, :libflint), Nothing, (Ref{gfp_poly}, Ref{fmpq_poly}), r, a)
end

function fmpq_poly_to_fmpz_mod_poly_raw!(r::fmpz_mod_poly, a::fmpq_poly, t1::fmpz_poly = fmpz_poly(), t2::fmpz = fmpz())
  ccall((:fmpq_poly_get_numerator, :libflint), Nothing, (Ref{fmpz_poly}, Ref{fmpq_poly}), t1, a)
  ccall((:fmpz_mod_poly_set_fmpz_poly, :libflint), Nothing, (Ref{fmpz_mod_poly}, Ref{fmpz_poly}), r, t1)
  ccall((:fmpq_poly_get_denominator, :libflint), Nothing, (Ref{fmpz}, Ref{fmpq_poly}), t2, a)
  if !isone(t2)
    res = ccall((:fmpz_invmod, :libflint), Cint, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), t2, t2, base_ring(r).modulus)
    @assert res != 0
    ccall((:fmpz_mod_poly_scalar_mul_fmpz, :libflint), Nothing, (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz}), r, r, t2)
  end
end

function fmpq_poly_to_gfp_fmpz_poly_raw!(r::gfp_fmpz_poly, a::fmpq_poly, t1::fmpz_poly = fmpz_poly(), t2::fmpz = fmpz())
  ccall((:fmpq_poly_get_numerator, :libflint), Nothing, (Ref{fmpz_poly}, Ref{fmpq_poly}), t1, a)
  ccall((:fmpz_mod_poly_set_fmpz_poly, :libflint), Nothing, (Ref{gfp_fmpz_poly}, Ref{fmpz_poly}), r, t1)
  ccall((:fmpq_poly_get_denominator, :libflint), Nothing, (Ref{fmpz}, Ref{fmpq_poly}), t2, a)
  if !isone(t2)
    res = ccall((:fmpz_invmod, :libflint), Cint, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), t2, t2, base_ring(r).modulus)
    @assert res != 0
    ccall((:fmpz_mod_poly_scalar_mul_fmpz, :libflint), Nothing, (Ref{gfp_fmpz_poly}, Ref{gfp_fmpz_poly}, Ref{fmpz}), r, r, t2)
  end
end

function fmpq_poly_to_nmod_poly(Rx::Nemo.NmodPolyRing, f::fmpq_poly)
  g = Rx()
  fmpq_poly_to_nmod_poly_raw!(g, f)
  return g
end

function fmpq_poly_to_gfp_poly(Rx::Nemo.GFPPolyRing, f::fmpq_poly)
  g = Rx()
  fmpq_poly_to_gfp_poly_raw!(g, f)
  return g
end

function fmpz_poly_to_nmod_poly_raw!(r::nmod_poly, a::fmpz_poly)
  ccall((:fmpz_poly_get_nmod_poly, :libflint), Nothing,
                  (Ref{nmod_poly}, Ref{fmpz_poly}), r, a)

end

function fmpz_poly_to_gfp_poly_raw!(r::gfp_poly, a::fmpz_poly)
  ccall((:fmpz_poly_get_nmod_poly, :libflint), Nothing,
                  (Ref{gfp_poly}, Ref{fmpz_poly}), r, a)

end

function fmpz_poly_to_nmod_poly(Rx::Nemo.NmodPolyRing, f::fmpz_poly)
  g = Rx()
  fmpz_poly_to_nmod_poly_raw!(g, f)
  return g
end

function fmpq_poly_to_fmpz_mod_poly(Rx::Nemo.FmpzModPolyRing, f::fmpq_poly)
  g = Rx()
  fmpq_poly_to_fmpz_mod_poly_raw!(g, f)
  return g
end

function fmpq_poly_to_gfp_fmpz_poly(Rx::Nemo.GFPFmpzPolyRing, f::fmpq_poly)
  g = Rx()
  fmpq_poly_to_gfp_fmpz_poly_raw!(g, f)
  return g
end

function fmpz_poly_to_fmpz_mod_poly_raw!(r::fmpz_mod_poly, a::fmpz_poly)
  ccall((:fmpz_poly_get_fmpz_mod_poly, :libflint), Nothing,
                  (Ref{fmpz_mod_poly}, Ref{fmpz_poly}), r, a)

end

function fmpz_poly_to_fmpz_mod_poly(Rx::Nemo.FmpzModPolyRing, f::fmpz_poly)
  g = Rx()
  fmpz_poly_to_fmpz_mod_poly_raw!(g, f)
  return g
end

################################################################################
#
#  Deflate and inflate
#
################################################################################

@doc Markdown.doc"""
    inflate(f::PolyElem, n::Int64) -> PolyElem
Given a polynomial $f$ in $x$, return $f(x^n)$, ie. multiply
all exponents by $n$.
"""
function inflate(x::PolyElem, n::Int64)
  y = parent(x)()
  for i=0:degree(x)
    setcoeff!(y, n*i, coeff(x, i))
  end
  return y
end

@doc Markdown.doc"""
    deflate(f::PolyElem, n::Int64) -> PolyElem
Given a polynomial $f$ in $x^n$, write it as a polynomial in $x$, ie. divide
all exponents by $n$.
"""
function deflate(x::PolyElem, n::Int64)
  y = parent(x)()
  for i=0:div(degree(x), n)
    setcoeff!(y, i, coeff(x, n*i))
  end
  return y
end

@doc Markdown.doc"""
    deflate(x::PolyElem) -> PolyElem
Deflate the polynomial $f$ maximally, ie. find the largest $n$ s.th.
$f$ can be deflated by $n$, ie. $f$ is actually a polynomial in $x^n$.
"""
function deflate(x::PolyElem)
  g = 0
  for i=0:degree(x)
    if coeff(x, i) != 0
      g = gcd(g, i)
      if g==1
        return x, 1
      end
    end
  end
  return deflate(x, g), g
end

################################################################################
#
#  change_base_ring
#
################################################################################

@doc Markdown.doc"""
    change_ring(f::PolyElem, Rx::PolyRing) -> PolyElem{elem_type(base_ring(Rx))}
Given a polynomial $f$ in $S[x]$, returns a polynomial in $Rx$ whose coefficient
are the coefficient of $f$ coerced to the base ring of $Rx$
"""
function change_ring(f::PolyElem, Rx::PolyRing)
  R = base_ring(Rx)
  coeffs = Vector{elem_type(R)}(undef, degree(f)+1)
  for i = 1:length(coeffs)
    coeffs[i] = R(coeff(f, i-1))
  end
  return Rx(coeffs)
end

@doc Markdown.doc"""
    change_base_ring(f::PolyElem, R::Ring) -> PolyElem{elem_type(R)}
Given a polynomial $f$ over a ring $S$, returns a polynomial whose coefficient
are the coefficient of $f$ coerced to $R$
"""
function change_base_ring(f::PolyElem, R::Ring)
  Rx, x = PolynomialRing(R, "_\$", cached = false)
  coeffs = Vector{elem_type(R)}(undef, degree(f)+1)
  for i = 1:length(coeffs)
    coeffs[i] = R(coeff(f, i-1))
  end
  return Rx(coeffs)
end

function change_base_ring(f::PolyElem, g::T) where T <: Union{Function, Map}
  h0 = g(coeff(f, 0))
  R = parent(h0)
  Rt, t = PolynomialRing(R, cached = false)
  coeffs = Vector{typeof(h0)}(undef, degree(f)+1)
  coeffs[0] = h0
  for i=1:degree(f)
    coeffs[i] = g(coeff(f, i))
  end
  return Rt(coeffs)
end

################################################################################
#
#  Reduced resultant
#
################################################################################

@doc Markdown.doc"""
    rres(f::fmpz_poly, g::fmpz_poly) -> fmpz
The reduced resultant of $f$ and $g$,
that is a generator for the ideal $(f, g) \cap Z$
"""
function rres(f::fmpz_poly, g::fmpz_poly)
  return rres_bez(f,g)
end

function rres_hnf(f::fmpz_poly, g::fmpz_poly)
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  s = Nemo.Generic.sylvester_matrix(f, g)
  h = hnf(s)
  return abs(h[nrows(h), ncols(h)])
end

function rres_bez(f::fmpz_poly, g::fmpz_poly)
  Nemo.check_parent(f, g)
  Qx = PolynomialRing(FlintQQ, "x", cached = false)[1]
  f1 = Qx(f)
  g1 = Qx(g)
  d, q, w = gcdx(f1, g1)
  if iszero(q) || iszero(w)
    if isconstant(f) || isconstant(g)
      if isconstant(f) && isconstant(g)
        return gcd(coeff(f, 0), coeff(g, 0))
      end
      if isconstant(f)
        if !isone(gcd(lead(g), coeff(f, 0)))
          cg = content(g - coeff(g, 0))
          ann = divexact(coeff(f, 0), gcd(coeff(f, 0), cg))
          return gcd(coeff(f, 0), ann*coeff(g, 0))
        else
          return coeff(f, 0)
        end
      end
      if !isone(gcd(lead(f), coeff(g, 0)))
        cf = content(f - coeff(f, 0))
        ann = divexact(coeff(g, 0), gcd(coeff(g, 0), cf))
        return gcd(coeff(g, 0), ann*coeff(f, 0))
      else
        return coeff(g, 0)
      end
    end
    return fmpz(0)
  end
  return lcm(denominator(q), denominator(w))
end

@doc Markdown.doc"""
    rresx(f::fmpz_poly, g::fmpz_poly) -> r, u, v
The reduced resultant, ie. a generator for the intersect
of the ideal generated by $f$ and $g$ with the integers.
As well as polynomials $u$ and $v$ s.th. $r = uf+vg$,
$\deg u < \deg g$ and $\deg v < \deg f$.
"""
function rresx(f::fmpz_poly, g::fmpz_poly)
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  Qx = PolynomialRing(FlintQQ, "x", cached = false)[1]
  g, q, w = gcdx(Qx(f), Qx(g))
  l = lcm(denominator(q), denominator(w))
  Zx = parent(f)
  return l, Zx(l*q), Zx(l*w)
end

struct PolyCoeffs{T <: RingElem}
  f::T
end

function coefficients(f::PolyElem)
  return PolyCoeffs(f)
end

function Base.iterate(PC::PolyCoeffs{<:PolyElem}, st::Int = -1)
  st += 1
  if st > degree(PC.f)
    return nothing
  else
    return coeff(PC.f, st), st
  end
end

Base.IteratorEltype(M::PolyElem) = Base.HasEltype()
Base.eltype(M::PolyElem{T}) where {T} = T

Base.IteratorSize(M::PolyCoeffs{<:PolyElem}) = Base.HasLength()
Base.length(M::PolyCoeffs{<:PolyElem}) = degree(M.f)+1

function Base.lastindex(a::PolyCoeffs{<:PolyElem})
  return degree(a.f)
end

function Base.getindex(a::PolyCoeffs{<:PolyElem}, i::Int)
  return coeff(a.f, i)
end

################################################################################
#
#  fmpq_poly with denominator 1 to fmpz_poly
#
################################################################################

function (a::FmpzPolyRing)(b::fmpq_poly)
  (!isone(denominator(b))) && error("Denominator has to be 1")
  z = a()
  ccall((:fmpq_poly_get_numerator, :libflint), Nothing,
              (Ref{fmpz_poly}, Ref{fmpq_poly}), z, b)
  return z
end

################################################################################
#
#  Power sums
#
################################################################################

@doc Markdown.doc"""
    polynomial_to_power_sums(f::PolyElem{T}, n::Int=degree(f)) -> Array{T, 1}
Uses Newton (or Newton-Girard) formulas to compute the first $n$
power sums from the coefficients of $f$.
"""
function polynomial_to_power_sums(f::PolyElem{T}, n::Int=degree(f)) where T <: FieldElem
  d = degree(f)
  R = base_ring(f)
  S = PowerSeriesRing(R, n+1, "gen(S)")[1]
  #careful: converting to power series and derivative do not commute
  #I also don't quite get this: I thought this was just the log,
  #but it isn't
  A = S([coeff(reverse(derivative(f)), i) for i=0:d-1], d, n+1, 0)
  B = S([coeff(reverse(f), i) for i=0:d], d+1, n+1, 0)
  L = A*inv(B)
  s = T[coeff(L, i) for i=1:n]
  return s
end

#plain vanilla recursion
function polynomial_to_power_sums(f::PolyElem{T}, n::Int=degree(f)) where T
  if n == 0
    return elem_type(base_ring(f))[]
  end
  d = degree(f)
  R = base_ring(f)
  E = T[(-1)^i*coeff(f, d-i) for i=0:min(d, n)] #should be the elementary symm.
  while length(E) <= n
    push!(E, R(0))
  end
  P = T[]

  push!(P, E[1+1])
  for k=2:n
    push!(P, (-1)^(k-1)*k*E[k+1] + sum((-1)^(k-1+i)*E[k-i+1]*P[i] for i=1:k-1))
  end
  return P
end

@doc Markdown.doc"""
    power_sums_to_polynomial(P::Array{T, 1}) -> PolyElem{T}
Uses the Newton (or Newton-Girard) identities to obtain the polynomial
coefficients (the elementary symmetric functions) from the power sums.
"""
function power_sums_to_polynomial(P::Array{T, 1}) where T <: FieldElem
  d = length(P)
  R = parent(P[1])
  S = PowerSeriesRing(R, d, "gen(S)")[1]
  s = S(P, length(P), d, 0)

  r = -integral(s)
  r1 = exp(r)
  #=
  if false
    r = S(T[R(1), -P[1]], 2, 2, 0)
    la = [d+1]
    while la[end]>1
      push!(la, div(la[end]+1, 2))
    end
    n = length(la)-1
    while n > 0
      set_prec!(r, la[n])
      rr = derivative(r)*inv(r)
      md = -(rr+s)
      m = S([R(1)], 1, la[n], 0)+integral(md)
      r *= m
      n -= 1
    end
    println("new exp $r")
  end
  =#
  @assert iszero(valuation(r1))
  Rx, x = PolynomialRing(R, cached = false)
  return Rx([Nemo.polcoeff(r1, d-i) for i=0:d])
end

function power_sums_to_polynomial(P::Array{T, 1}) where T
  E = T[1]
  R = parent(P[1])
  last_non_zero = 0
  for k=1:length(P)
    push!(E, divexact(sum((-1)^(i-1)*E[k-i+1]*P[i] for i=1:k), R(k)))
    if E[end] != 0
      last_non_zero = k
    end
  end
  E = E[1:last_non_zero+1]
  d = length(E) #the length of the resulting poly...
  for i=1:div(d, 2)
    E[i], E[d-i+1] = (-1)^(d-i)*E[d-i+1], (-1)^(i-1)*E[i]
  end
  return PolynomialRing(R, cached = false)[1](E)
end

##############################################################
# all of this should be in Nemo/AbstractAlgebra
#
#TODO:
# expand systematically for all finite fields
# and for fmpz/fmpq poly
# for fun: ispower(a::nf_elem)
#

function factor(f::PolyElem, R::Field)
  f1 = change_base_ring(f, R)
  return factor(f1)
end

function factor(f::fmpq_poly, R::T) where T <: Union{Nemo.FqNmodFiniteField, Nemo.GaloisField}
  Rt, t = PolynomialRing(R, "t", cached=false)
  return factor(Rt(f))
end

function roots(f::PolyElem, R::Field)
  f1 = change_base_ring(f, R)
  return roots(f1)
end

function roots(f::fmpq_poly, R::T) where T <: Union{Nemo.FqNmodFiniteField, Nemo.GaloisField}
  Rt, t = PolynomialRing(R, "t", cached=false)
  fp = PolynomialRing(FlintZZ, cached = false)[1](f*denominator(f))
  fpp = Rt(fp)
  return roots(fpp)
end

function roots(f::gfp_poly, K::FqNmodFiniteField)
  @assert characteristic(K) == characteristic(base_ring(f))
  Kx = PolynomialRing(K, cached = false)[1]
  ff = Kx()
  for i=0:degree(f)
    setcoeff!(ff, i, lift(coeff(f, i)))
  end
  return roots(ff)
end

function ispower(a::Nemo.fq_nmod, m::Int)
  s = size(parent(a))
  if gcd(s-1, m) == 1
    return true, a^invmod(m, s-1)
  end
  St, t = PolynomialRing(parent(a), "t", cached=false)
  f = t^m-a
  rt = roots(f)
  if length(rt) > 0
    return true, rt[1]
  else
    return false, a
  end
end

function roots(f::T) where T <: Union{fq_nmod_poly, fq_poly} # should be in Nemo and
                                    # made available for all finite fields I guess.
  q = size(base_ring(f))
  x = gen(parent(f))
  if degree(f) < q
    x = powmod(x, q, f)-x
  else
    x = x^Int(q)-x
  end
  f = gcd(f, x)
  l = factor(f).fac
  return elem_type(base_ring(f))[-trailing_coefficient(x) for x = keys(l) if degree(x)==1]
end

# generic fall back
# ...
function roots(f::PolyElem)
  lf = factor(f)
  return elem_type(base_ring(f))[-trailing_coefficient(x) for x= keys(lf.fac) if degree(x)==1]
end

function ispower(a::RingElem, n::Int)
  if isone(a) || iszero(a)
    return true, a
  end
  if isone(-a) && isodd(n)
    return true, a
  end
  R = parent(a)
  Rt, t = PolynomialRing(R, "t", cached = false)
  r = roots(t^n-a)
  if length(r) == 0
    return false, a
  else
    return true, r[1]
  end
end

function root(a::RingElem, n::Int)
  fl, b = ispower(a, n)
  fl || error("element does not have a $n-th root")
  return b
end

function setcoeff!(z::fq_nmod_poly, n::Int, x::fmpz)
   ccall((:fq_nmod_poly_set_coeff_fmpz, :libflint), Nothing,
         (Ref{fq_nmod_poly}, Int, Ref{fmpz}, Ref{FqNmodFiniteField}),
         z, n, x, base_ring(parent(z)))
     return z
end

###############################################################################
#
#  Sturm sequence
#
###############################################################################

#See Wikipedia as a reference
function _divide_by_content(f::fmpz_poly)

  p = primpart(f)
  if sign(lead(f))== sign(lead(p))
    return p
  else
    return -p
  end
end

function sturm_sequence(f::fmpz_poly)

  g = f
  h = _divide_by_content(derivative(g))
  seq = fmpz_poly[g,h]
  while true
    r = _divide_by_content(pseudorem(g,h))
    if r != 0
      push!(seq, -r)
      g, h = h, -r
    else
      break
    end
  end
  return seq

end

function _number_changes(a::Array{Int,1})

  nc = 0
  filter!(x -> x != 0, a)
  for i = 2:length(a)
    if sign(a[i]) != sign(a[i-1])
      nc += 1
    end
  end
  return nc

end

function number_positive_roots(f::fmpz_poly)

  s = sturm_sequence(f)
  evinf = Int[sign(coeff(x, degree(x))) for x in s]
  ev0 = Int[sign(coeff(x,0)) for x in s]
  return _number_changes(ev0)-_number_changes(evinf)

end

function number_positive_roots(f::PolyElem{nf_elem}, P::InfPlc)
  fsq = squarefree_factorization(f)
  p = 0
  for (g, e) in fsq
    p = p + _number_positive_roots_sqf(g, P) * e
  end
  return p
end

function _number_positive_roots_sqf(f::PolyElem{nf_elem}, P::InfPlc; start_prec::Int = 32)
  prec = start_prec
  while true
    coeffs = Vector{acb}(undef, length(f))
    c = evaluate(coeff(f, 0), P, prec)
    coeffs[1] = c
    C = parent(c)
    Cx, x = PolynomialRing(C, "x", cached = false)
    for i in 1:degree(f)
      coeffs[i + 1] = evaluate(coeff(f, i), P, prec)
    end
    g = Cx(coeffs)
    rts = real.(Hecke.roots(g))
    if any(contains_zero, rts)
      prec = 2 * prec
    else
      return count(ispositive, rts)
    end
  end
end

################################################################################
#
#  Squarefree factorization in characteristic 0
#
################################################################################

# This is Musser's algorithm
function squarefree_factorization(f::PolyElem)
  @assert iszero(characteristic(base_ring(f)))
  c = lead(f)
  f = divexact(f, c)
  di = gcd(f, derivative(f))
  ei = divexact(f, di)
  i = 1
  res = Dict{typeof(f), Int}()
  while !isconstant(ei)
    eii = gcd(di, ei)
    dii = divexact(di, eii)
    if degree(eii) != degree(ei)
      res[divexact(ei, eii)] = i
    end
    i = i +1
    di = dii
    ei = eii
  end
  return Fac(parent(f)(c), res)
end

################################################################################
#
#  Squarefree factorization for fmpq_poly
#
################################################################################

@doc Markdown.doc"""
    factor_squarefree(x::fmpq_poly)
Returns the squarefree factorization of $x$.
"""
function factor_squarefree(x::fmpq_poly)
   res, z = _factor_squarefree(x)
   return Fac(parent(x)(z), res)
end

function _factor_squarefree(x::fmpq_poly)
   res = Dict{fmpq_poly, Int}()
   y = fmpz_poly()
   ccall((:fmpq_poly_get_numerator, :libflint), Nothing,
         (Ref{fmpz_poly}, Ref{fmpq_poly}), y, x)
   fac = Nemo.fmpz_poly_factor()
   ccall((:fmpz_poly_factor_squarefree, :libflint), Nothing,
              (Ref{Nemo.fmpz_poly_factor}, Ref{fmpz_poly}), fac, y)
   z = fmpz()
   ccall((:fmpz_poly_factor_get_fmpz, :libflint), Nothing,
            (Ref{fmpz}, Ref{Nemo.fmpz_poly_factor}), z, fac)
   f = fmpz_poly()
   for i in 1:fac.num
      ccall((:fmpz_poly_factor_get_fmpz_poly, :libflint), Nothing,
            (Ref{fmpz_poly}, Ref{Nemo.fmpz_poly_factor}, Int), f, fac, i - 1)
      e = unsafe_load(fac.exp, i)
      res[parent(x)(f)] = e
   end
   return res, fmpq(z, denominator(x))

end

function charpoly_mod(M::Generic.Mat{nf_elem}; integral::Bool = false, normal::Bool = false, proof::Bool = true)
  K = base_ring(M)
  p = p_start
  Kt, t = PolynomialRing(K, cached = false)
  f = Kt()
  f_last = f
  d = fmpz(1)
  stable = 5
  max_stable = 5
  while true
    p = next_prime(p)
    me = modular_init(K, p)
    if normal && length(me.fld) < degree(K)
      continue
    end
    t = Hecke.modular_proj(M, me)
    if !isdefined(me, :fldx)
      me.fldx = [PolynomialRing(x, "_x", cached=false)[1] for x = me.fld]
      me.Kx = Kt
    end

    fp = map(i-> charpoly(me.fldx[i], t[i]), 1:length(t))
    gp= Hecke.modular_lift(fp, me)
    if iszero(f)
      f = gp
    else
      f, d = induce_crt(f, d, gp, fmpz(p), true)
      if integral
        fl = true
        gg = f
      else
        fl, gg = induce_rational_reconstruction(f, d)
      end

      if fl && gg == f_last
        stable -= 1
        if stable <= 0
          break
        end
      else
        stable = max_stable
      end
      f_last = gg
    end
  end
  if !proof
    return f_last
  end
  error("Proof not implemented")
end
#=
function cyclic_generators(A::MatElem{T}) where {T <: FieldElem}
  b = matrix(base_ring(A), 0, nrows(A), [])
  g = matrix(base_ring(A), 0, nrows(A), [])
  while nrows(b) < nrows(A)
    if nrows(g) == 0
      g = zero_matrix(base_ring(A), 1, nrows(A))
      g[1,1] = 1
    else
      i = findfirst(j-> b[j,j] == 0, 1:nrows(b))
      if i == nothing
        i = nrows(b)+1
      end
      g = vcat(g, zero_matrix(base_ring(A), 1, nrows(A)))
      g[nrows(g), i] = 1
    end
    b = extend_cyclic_subspace(A::MatElem{T}, b::MatElem{T}, g)
    if nrows(b) == nrows(A)
      return g
    end
  end
end

function extend_cyclic_subspace(A::MatElem{T}, b::MatElem{T}, g) where {T <: FieldElem}
  while true
    g = vcat(g, g*A)
    cleanvect(b, g) #currently does only single rows...
    i = findfirst(i->iszero_row(g, i), 1:nrows(g))
    if i != nothing
      b = vcat(b, view(g, 1:i-1, 1:ncols(g)))
      rk, b = rref!(b)
      @assert rk == nrows(b)
      return b
    end
    A *= A
  end
end
function cyclic_subspace(A::MatElem{T}, b::MatElem{T}) where {T <: FieldElem}
  b = deepcopy!(b)
  rk, b = rref!(b)
  bv = view(b, 1:rk, 1:ncols(b))
  if rk == 0 || rk == ncols(b)
    return bv
  end
  while true
    b2 = bv*A
    b = vcat(bv, b2)
    rk_new, b = rref!(b)
    if rk_new == rk
      return bv
    end
    rk= rk_new
    bv = view(b, 1:rk, 1:ncols(b))
    if rk == ncols(b)
      return bv
    end
    A *= A
  end
end
=#
#=
  plan for proof:
    if f is irreducible (or at least square-free), then there are
      (many) primes p s.th. f is square-free mod p
    then that means there are vectors b s.th. the
    space <M^i b | i> = everyhting, at least mod p, so in general.
    Now f(M)b = 0 implies f(M) = 0.

    if f is known to be integral, then one can use arb to compute the
      complex version and use this to derive bounds...

    There are also bounds on the coefficients which are sometimes tight
=#

#=
function roots(f::fmpz_poly)
  g = gcd(f, derivative(f))
  p = p_start
  while true
    p = next_prime(p)
    R = GF(p)
    Rx,x = PolynomialRing(R, cached = false)
    gp = Rx(g)
    if !issquarefree(gp)
      continue
    end
    #TODO: try a few primes to find best (with fewest roots)
    rp = roots(gp)
    if length(rp) == 0
      return []
    end

  end

end

=#


text:
# Similar function to quadprog routine in MATLAB. Here f(x) = (1/2)*x'*H*x + f'*x,
# where H is a matrix, f a vector and we have equality as well as
# boundary restrictions such as Ax <= b, (Aeq)*x = beq, lb <= x <= ub.

# This routine uses the Ipopt Solver, speeding up the optimization.
# but one can use other solvers, just change the usesolver as input of the
# function

using JuMP, Ipopt

function quadprog(H,f; A=[],b=[],Aeq = [],beq = [],lb = [],ub = [],
  usesolver = IpoptSolver(), starter = zeros(size(H,2)))

  r,s = size(H)

  m = Model(solver = usesolver)

  #Initialize the point, you can specify the initial point as you want
  #using the starter input on the function.

  if starter == zeros(size(H,2))
    @variable(m,x[1:s])
  else
    @variable(m, x[i=1:s], start = starter[i])
  end

  #Creates inequality constraint, if A or b are different from the empty array

  if !(isempty(A)) || !(isempty(b))
    @constraintref loopConstraint1[1:length(b)]
    for i = 1:length(b)
      loopConstraint1[i] = @constraint(m,dot(A[i,:],x) <= b[i])
    end
  end

  #Creates equality constraint if it exists.

  if !(isempty(Aeq)) || !(isempty(beq))
    @constraintref loopConstraint2[1:length(beq)]
    for i = 1:length(beq)
      loopConstraint2[i] = @constraint(m, dot(Aeq[i,:],x) == beq[i])
    end
  end

  #Creates the lower bound constraint if it exists

  if !(isempty(lb))
    @constraintref loopConstraint3[1:length(lb)]
    for i = 1:length(lb)
      loopConstraint3[i] = @constraint(m,x[i] >= lb[i])
    end
  end

  #Creates the upper bound constraint if it exists

  if !(isempty(ub))
    @constraintref loopConstraint3[1:length(ub)]
    for i = 1:length(ub)
      loopConstraint4[i] = @constraint(m,x[i] <= ub[i])
    end
  end

  #Specify the objective function, just as in Matlab

  @objective(m, Min, 0.5*dot(H*x,x) + dot(f,x))

  #print the problem on the REPL

  print(m)

  #solve the minimization problem

  solve(m)

  #return the best solution given by the solver

  return getvalue(x)

end

max_stars_repo_name: chronchi/Quadprog

max_issues_repo_name: chronchi/Quadprog

max_forks_repo_name: chronchi/Quadprog


text:
using Base: @deprecate

@deprecate ci(args...) confint(args...)

@deprecate MultinomialLRT MultinomialLRTest
@deprecate OneSampleHotellingT2 OneSampleHotellingT2Test
@deprecate EqualCovHotellingT2 EqualCovHotellingT2Test
@deprecate UnequalCovHotellingT2 UnequalCovHotellingT2Test

@deprecate BartlettsTest BartlettTest

@deprecate confint(x::HypothesisTest, alpha::Real; kwargs...) confint(x; level=1-alpha, kwargs...)


text:
using BenchmarkTools
using SolidStateDetectors
using Plots
using Unitful

using SolidStateDetectors: rb_even, rb_odd

### Inner loop of SOR
T = Float32
sim_car = Simulation{T}(SSD_examples[:CGD]);

calculate_electric_potential!(
    sim_car, depletion_handling = false,
    convergence_limit = 1e-4,
    max_tick_distance = 1u"mm",
    refinement_limits = [0.2, 0.1],
    # refinement_limits = missing,
); plot(sim_car.electric_potential, y = 0)

# plot(
#     plot(sim_car.point_types, y = 0),
#     plot(sim_car.electric_potential, y = 0),
# )

pcs_car = SolidStateDetectors.PotentialCalculationSetup(
    sim_car.detector,
    sim_car.electric_potential.grid,
    # Grid(sim_car, max_tick_distance = 1u"mm"),
    SolidStateDetectors.material_properties[SolidStateDetectors.materials["vacuum"]],
    sim_car.electric_potential.data
);

nthreads = 1 #Base.Threads.nthreads()
is_weighting_potential = Val{false}()
depletion_handling = Val{false}()
only2d = Val{false}()
update_even_points = Val{true}()
update_uneven_points = Val{false}()

idx3 = 2
even_points = true
update_even_points = Val(even_points)
is_weighting_potential = Val(false)
rb_tar_idx, rb_src_idx = even_points ? (rb_even, rb_odd) : (rb_odd, rb_even)
depletion_handling = Val{false}()
idx3iseven = Val(iseven(idx3))
idx3istwo = Val(idx3 == 2)

# calculate_weights_for_innerloop! & innerloop!
depletion_handling = Val{false}()

idx3_is_even = Val(iseven(idx3))

iz = idx3
inz = iz - 1

pwwzr        = pcs_car.geom_weights[3][1, inz]
pwwzl        = pcs_car.geom_weights[3][2, inz]
pwΔmpz       = pcs_car.geom_weights[3][3, inz]
Δz_ext_inv_r = pcs_car.geom_weights[3][4, inz + 1]
Δz_ext_inv_l = pcs_car.geom_weights[3][4, inz]

line_weights = Array{T, 2}(undef, size(pcs_car.potential, 1) - 2, 6)

iy = 2
iny = iy - 1

pwwyr  = pcs_car.geom_weights[2][1, iny]
pwwyl  = pcs_car.geom_weights[2][2, iny]
pwΔmpy = pcs_car.geom_weights[2][3, iny]
pwΔmpy_pwΔmpz = pwΔmpy * pwΔmpz
Δy_ext_inv_r_pwΔmpz  = pcs_car.geom_weights[2][4, iny + 1] * pwΔmpz
Δy_ext_inv_l_pwΔmpz  = pcs_car.geom_weights[2][4, iny]     * pwΔmpz
Δz_ext_inv_r_pwΔmpy = Δz_ext_inv_r * pwΔmpy
Δz_ext_inv_l_pwΔmpy = Δz_ext_inv_l * pwΔmpy

pwwyr_pwwzr = pwwyr * pwwzr
pwwyr_pwwzl = pwwyr * pwwzl
pwwyl_pwwzr = pwwyl * pwwzr
pwwyl_pwwzl = pwwyl * pwwzl

SolidStateDetectors.middleloop!(idx3, rb_tar_idx, rb_src_idx, pcs_car,
                        update_even_points, depletion_handling, is_weighting_potential, only2d, idx3iseven, idx3istwo)

@code_warntype SolidStateDetectors.middleloop!(idx3, rb_tar_idx, rb_src_idx, pcs_car,
                        update_even_points, depletion_handling, is_weighting_potential, only2d, idx3iseven, idx3istwo)

@btime SolidStateDetectors.middleloop!($idx3, $rb_tar_idx, $rb_src_idx, $pcs_car,
                        $update_even_points, $depletion_handling, $is_weighting_potential, $only2d, $idx3iseven, $idx3istwo)

SolidStateDetectors.calculate_weights_for_innerloop!(line_weights, pcs_car, iy, iny, iz, inz,
            update_even_points, idx3_is_even,
            pwwzr, pwwzl, pwwyr, pwwyl,
            pwwyr_pwwzr, pwwyr_pwwzl, pwwyl_pwwzr, pwwyl_pwwzl,
            pwΔmpy_pwΔmpz,
            Δz_ext_inv_r_pwΔmpy, Δz_ext_inv_l_pwΔmpy,
            Δy_ext_inv_r_pwΔmpz, Δy_ext_inv_l_pwΔmpz)

@btime SolidStateDetectors.calculate_weights_for_innerloop!($line_weights, $pcs_car, $iy, $iny, $iz, $inz,
            $update_even_points, $idx3_is_even,
            $pwwzr, $pwwzl, $pwwyr, $pwwyl,
            $pwwyr_pwwzr, $pwwyr_pwwzl, $pwwyl_pwwzr, $pwwyl_pwwzl,
            $pwΔmpy_pwΔmpz,
            $Δz_ext_inv_r_pwΔmpy, $Δz_ext_inv_l_pwΔmpy,
            $Δy_ext_inv_r_pwΔmpz, $Δy_ext_inv_l_pwΔmpz)

SolidStateDetectors.innerloop!(line_weights, pcs_car, iy, iny, iz, inz, rb_tar_idx, rb_src_idx,
            update_even_points, idx3_is_even,
            depletion_handling, is_weighting_potential, only2d)

@code_warntype SolidStateDetectors.innerloop!(line_weights, pcs_car, iy, iny, iz, inz, rb_tar_idx, rb_src_idx,
            update_even_points, idx3_is_even,
            depletion_handling, is_weighting_potential, only2d)

@btime SolidStateDetectors.innerloop!($line_weights, $pcs_car, $iy, $iny, $iz, $inz, $rb_tar_idx, $rb_src_idx,
            $update_even_points, $idx3_is_even,
            $depletion_handling, $is_weighting_potential, $only2d)

@code_llvm SolidStateDetectors.innerloop!(line_weights, pcs_car, iy, iny, iz, inz, rb_tar_idx, rb_src_idx,
            update_even_points, idx3_is_even,
            depletion_handling, is_weighting_potential, only2d)

using KernelAbstractions
using CUDAKernels, SolidStateDetectors, Unitful
using CUDAKernels.CUDA: CuArray

T = Float32
sim_car = Simulation{T}(SSD_examples[:CGD]);

calculate_electric_potential!(
    sim_car, depletion_handling = false,
    device_array_type = CuArray,
    convergence_limit = 1e-4,
    max_tick_distance = 1u"mm",
    refinement_limits = [0.2, 0.1],
    max_n_iterations = 2000,
    # refinement_limits = missing,
)

begin
    @info "Grid size: $(size(sim.electric_potential.data))"
    # SolidStateDetectors.outerloop!(pcs_car, nthreads, update_even_points, depletion_handling, is_weighting_potential, only2d)
    for dp in (false, true)
        depletion_handling = Val{dp}()
        @info "Depletion handling: $dp"
        for nt in (1, 2, 4, 8, 16, 32, 64)
            @info "N Threads: $nt"
            @btime SolidStateDetectors.outerloop!($pcs_car, $nt, $update_even_points, $depletion_handling, $is_weighting_potential, $only2d)
        end
    end
end


text:
using Documenter, openBF

makedocs(
    format = :html,
    sitename = "openBF",
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Guide" => "man/guide.md",
            "man/config.md",
            "man/examples.md"
        ]
    ])


text:
"""
    LayerCache()
Holds image matrices of layer frames in case [`show_layer_frame`](@ref) is called.

# Arguments
- `frames::UnitRange`:The frame range for which layer's frames are to be viewed
- `frame_counter::Int`: internal counter to keep count of the layer's frame being placed
- `layer_frames::Union{UnitRange,Int}`: The frame/frames of the layer to be viewed
- `position::Array{Point}`: The position of the layer at each frame that is to be viewed
- `settings_cache::Array{LayerSetting}` : Tha layer settings of each frame of the layer to be viewed
- `matrix_cache::Array`: a list that holds the image matrices of the layer frames to be viewed
"""
mutable struct LayerCache
    frames::UnitRange
    frame_counter::Int
    layer_frames::Union{UnitRange,Int}
    # todo while dropping support for julia 1.5
    # change type to Base.ReinterpretArray{ARGB32,2,UInt32,Matrix{UInt32},false}
    position::Array{Point}
    settings_cache::Array{LayerSetting}
    matrix_cache::Array

    LayerCache() = new(0:0, 0, 0, [], [], [])
end


text:
# constant used for computation
EPSSCALE = 10

# Computes the Wasserstein distance between a univariate discrete distribution
# and a univariate distribution with an implemented function 'cdf'.
#
# Args:
#   Q - univariate SteinDiscrete distribution; support of Q must be sorted
#   targetcdf - CDF of the (other) univariate distribution P
#
# Returns:
#   Tuple of Wasserstein distance and upper bound on numerical integration
#   error.
function wasserstein1d(Q::SteinDiscrete,
                       targetcdf::Function;
                       lower=-Inf,
                       upper=Inf)
    X = Q.support;
    if size(X,2) != 1
        error("function only defined by univariate distributions")
    end
    if !issorted(X)
        error("Support of Q must be sorted")
    end
    q = Q.weights;
    n = length(q);

    # Compute \int_t |Q(t) - P(t)| dt where Q(t) is the cdf of Q and P(t)
    # is the cdf of P
    # First compute contribution of the interval (lower, x_1)
    # We need to be watchful for machine precision issues
    (integral, max_error) = quadgk(t->thresholdpos(targetcdf(t), EPSSCALE*eps(0.0)), lower, X[1])
    # Cumulative Q weight
    weight = 0;
    for i in 1:(n-1)
        # Compute integral contributed by each (x_i, x_{i+1}] interval
        weight = weight + q[i];
        (estimate, error) = quadgk(t->thresholdpos(abs(targetcdf(t) - weight), EPSSCALE*eps(weight)), X[i], X[i+1]);
        integral = integral + estimate;
        max_error = max_error + error;
    end
    # Add contribution of (x_n, upper)
    (estimate, error) = quadgk(t->thresholdpos(1-targetcdf(t), EPSSCALE*eps(-1.0)), X[n], upper)
    integral = integral + estimate;
    max_error = max_error + error;
    (integral, max_error)
end

# Computes the Wasserstein distance between a univariate discrete distribution
# (represented by points and weights) and a univariate distribution with an
# implemented function 'cdf'.
#
# Args:
#   points - n x 1 array of support points
#   weights - n x 1 array of real-valued weights associated with support points
#   target - univariate SteinDistribution implementing the function 'cdf'
#
# Returns:
#   Tuple of Wasserstein distance and upper bound on numerical integration
#   error.
function wasserstein1d(points::Array{Float64};
                       weights=fill(1/size(points,1), size(points,1)),
                       targetcdf=nothing)
    isa(targetcdf, Function) ||
        error("targetcdf must be the cdf of the continuous distribution.")
    wasserstein1d(SteinDiscrete(points, weights), targetcdf)
end

# This is a helper function that used to avoid weird
# machine precision issues when integrating the cdf.
function thresholdpos(x::Float64, threshold::Float64)
    (x <= threshold) ? 0.0 : x
end

# This function computes the earth movers distance between two
# discrete distributions.
#
# It solves a large LP in order to compute the EMD.
function wassersteindiscrete{T<:Number}(;
    xpoints::AbstractArray{T,2}=[],
    xweights::Array{Float64,1}=fill(1.0/size(xpoints,1), size(xpoints,1)),
    ypoints::AbstractArray{T,2}=[],
    yweights::Array{Float64,1}=fill(1.0/size(ypoints,1), size(ypoints,1)),
    solver=nothing)
    # make sure samples are same dimensionality
    @assert size(xpoints, 2) == size(ypoints, 2)
    # make sure solver is defined
    isa(solver, AbstractMathProgSolver) ||
        error("Must specify solver of type AbstractMathProgSolver")
    # setup discrete distributions
    xdist = SteinDiscrete(xpoints, xweights)
    nx = size(xdist.support, 1)
    ydist = SteinDiscrete(ypoints, yweights)
    ny = size(ydist.support, 1)
    # prepare the manhattan extension graph
    complete_support = [xdist.support; ydist.support]
    graph_support, edgeset = makel1spanner(complete_support)
    # compute the supply for the network problem
    ntotal = size(graph_support,1)
    qx = zeros(ntotal)
    qx[1:nx] = xdist.weights
    qy = zeros(ntotal)
    qy[(nx+1):(nx+ny)] = ydist.weights
    # now solve the min-cost max-flow problem
    (emd, flow, status) = mincostflow(graph_support,
                                      edgeset,
                                      qx,
                                      qy,
                                      solver)
    return (emd, size(graph_support,1), length(edgeset), status)
end

# This function computes the earth movers distance between two discrete
# univariate distributions.
# This is simply int |F1(z) - F2(z)| dz.
function wasserstein1d{T<:Number}(;
    xpoints::AbstractArray{T,1}=[],
    xweights::Array{Float64,1}=fill(1.0/length(xpoints), length(xpoints)),
    ypoints::AbstractArray{T,1}=[],
    yweights::Array{Float64,1}=fill(1.0/size(ypoints,1), size(ypoints,1)))
    # setup discrete distributions
    xdist = SteinDiscrete(xpoints, xweights)
    ydist = SteinDiscrete(ypoints, yweights)
    # now grab the goodies
    xsupp = xdist.support
    qx = xdist.weights
    ysupp = ydist.support
    qy = ydist.weights
    # initialize the indices we walk up the vectors
    ix = 0; iy = 0
    wasserstein = 0.0
    Fx = 0.0
    Fy = 0.0
    currval = -Inf
    if xsupp[1] <= ysupp[1]
        ix = 1
        Fx = qx[1]
        currval = xsupp[1]
    end
    if xsupp[1] >= ysupp[1]
        iy = 1
        Fy = qy[1]
        currval = ysupp[1]
    end
    while ix != length(xsupp) || iy != length(ysupp)
        # handle edge case when we've cycled through the distribution
        if ix == length(xsupp)
            iy += 1
            wasserstein += (ysupp[iy] - currval) * (1.0 - Fy)
            Fy += qy[iy]
            currval = ysupp[iy]
            continue
        end
        if iy == length(ysupp)
            ix += 1
            wasserstein += (xsupp[ix] - currval) * (1.0 - Fx)
            Fx += qx[ix]
            currval = xsupp[ix]
            continue
        end
        # now do the look ahead
        nextval = min(xsupp[ix+1], ysupp[iy+1])
        wasserstein += (nextval - currval) * abs(Fx - Fy)
        if xsupp[ix+1] <= ysupp[iy+1]
            ix += 1
            Fx += qx[ix]
            currval = xsupp[ix]
        else
            iy +=1
            Fy += qy[iy]
            currval = ysupp[iy]
        end
    end
    (wasserstein, ix + iy, ix + iy - 1, :Optimal)
end

# This function approximates the wasserstein metric by using two facts:
# 1. We can compute the EMD for two different discrete distributions.
# 2. The empirical wasserstein metric is root{n}-consistent for the
# true wasserstein metric.
function approxwasserstein{T<:Number}(;
    points::AbstractArray{T,2}=[],
    weights::Array{Float64,1}=fill(1.0/size(points,1), size(points,1)),
    targetsamplegen=nothing,
    replicates::Int=30,
    alpha::Float64=0.95,
    solver=nothing)
    # make sure targetsamplegen is given
    isa(targetsamplegen, Function) ||
        error("Must supply targetsamplegen as a function to generate random samples.")
    # compute a few estimates of wasserstein metric
    wassersteinestimates = Array{Float64}(replicates)
    for ii in 1:replicates
        sample = targetsamplegen()
        (emd, numnodes, numedges, status) =
            wassersteindiscrete(xpoints=points,
                                xweights=weights,
                                ypoints=sample,
                                solver=solver)
        wassersteinestimates[ii] = emd
    end
    # compute CI
    noise = Distributions.TDist(replicates - 1)
    wassmean = mean(wassersteinestimates)
    wassstd = std(wassersteinestimates) / sqrt(replicates)
    zalpha = abs(Distributions.quantile(noise, (1.0 - alpha) / 2))
    return (wassmean - zalpha * wassstd, wassmean + zalpha * wassstd)
end


text:
import Base: convert

mutable struct MessageIterator
    data::Vector{UInt8}
    messages::OrderedDict{Int64, OrderedDict{Int64, String}}
    tag::Vector{Char}
    value::Vector{Char}
    state::Bool
    message_id::Int64
    function MessageIterator(data::Vector{UInt8})
        return new(data, OrderedDict{Int64, OrderedDict{Int64, String}}(), Char[], Char[], true, 0)
    end
end

length(fix::MessageIterator) = length(fix.data)

function start(fix::MessageIterator)
    return 1
end

function next(fix::MessageIterator, idx::Int64)::Tuple{MessageIterator, Int64}
    x = fix.data[idx]
    if x == 0x01
        int_tag = parse(Int64, String(fix.tag))
        str_value = String(fix.value)
        !haskey(fix.messages, fix.message_id) && (fix.messages[fix.message_id] = OrderedDict{Int64, String}())
        fix.messages[fix.message_id][int_tag] = str_value
        resize!(fix.tag, 0)
        resize!(fix.value, 0)
        if int_tag == 10
            fix.message_id += 1
        end
        fix.state = true
    elseif x == 0x3d
        fix.state = !fix.state
    else
        if fix.state
            push!(fix.tag, Char(x))
        else
            push!(fix.value, Char(x))
        end
    end

    return (fix, idx + 1)
end

function done(fix::MessageIterator, idx::Int64)
    return length(fix.data) > idx
end

function parse(data::Vector{UInt8})::OrderedDict{Int64, OrderedDict{Int64, String}}
    f = MessageIterator(data)
    for k in eachindex(data)
        next(f, k)
    end
    return f.messages
end

function convert(tags::Dict{Int64, String}, msg::OrderedDict{Int64, String})::OrderedDict{String, String}
    out = OrderedDict{String, String}()
    for (k, v) in msg
        out[tags[k]] = v
    end
    return out
end


text:
# Autogenerated wrapper script for Readline_jll for i686-linux-musl
export libhistory, libreadline

using Ncurses_jll
JLLWrappers.@generate_wrapper_header("Readline")
JLLWrappers.@declare_library_product(libhistory, "libhistory.so.8")
JLLWrappers.@declare_library_product(libreadline, "libreadline.so.8")
function __init__()
    JLLWrappers.@generate_init_header(Ncurses_jll)
    JLLWrappers.@init_library_product(
        libhistory,
        "lib/libhistory.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "1001.grid-illumination.jl" begin
    @test grid_illumination(5, [(0,0),(4,4)], [(1,1),(1,1)]) == [1, 1]
    @test grid_illumination(5, [(0,0),(4,4)], [(1,1),(1,0)]) == [1, 0]
    @test grid_illumination(5, [(0,0),(0,4)], [(0,4),(0,1),(1,4)]) == [1, 1, 0]
end


text:
function load_dream_challenge_data()
	directory = joinpath(PROJECT_ROOT, "data", "dream_challenge")
	log_message("importing dream challenge data from $directory")
	experiment = create_experiment("dream_challenge")
	for f in readdir(directory)
		if !endswith(f,".txt") || endswith(f, "README.txt")
			continue
		end
		data_type_string = replace(f, "DREAM7_DrugSensitivity1_" => "")
		data_type_string = replace(data_type_string, ".txt" => "")
		key_type = Gene
		if data_type_string == "Exomeseq"
			data_type = ExomeSeq
		elseif data_type_string == "GeneExpression"
			data_type = GeneExpression
		elseif data_type_string == "Methylation"
			data_type = Methylation
		elseif data_type_string == "RNAseq_expressed_calls"
			data_type = RNASeqCall
		elseif data_type_string == "RNAseq_quantification"
			data_type = RNASeq
		elseif data_type_string == "RPPA"
			key_type = Protein
			data_type = RPPA
		elseif data_type_string == "SNP6_gene_level"
			data_type = CNV
		elseif data_type_string == "Drug_Response_Training" || data_type_string == "test_data"
			# responses
		else
			log_message("couldn't determine data type $data_type_string")
			continue
		end
		log_message("importing file $f")
		# make nullable
		df = CSV.read(joinpath(directory, f), DataFrame, delim='\t',rows_for_type_detect=1200, missingstring="NA", strings=:raw)
		if data_type_string == "Drug_Response_Training"
			df = CSV.read(joinpath(directory, f), DataFrame, delim='\t',rows_for_type_detect=1200, missingstring="NA", strings=:raw)
		end
		log_message("creating genes or proteins...")
		if data_type_string == "Drug_Response_Training" || data_type_string == "test_data"
			drugs = names(df)[2:end]
			# println(typeof(df[:CellLine]))
			cell_lines = df[:CellLine]
			cl_objects = map(cl_id -> get_cell_line!(experiment, cl_id, "Breast cancer"), cell_lines)
			for d in drugs
				vals = df[d]
				o = Outcome(string(d), "IC50")
				add_results!(o, cl_objects[.!ismissing.(df[d])], Float64.(vals[.!ismissing.(df[d])]))
				drug = get!(experiment.drugs, string(d), Drug(string(d)))
				add_outcome(experiment, drug, o)
			end
			if data_type_string == "test_data"
				[cl.in_test_set = true for cl in cl_objects]
			end
			continue
		elseif data_type_string == "test_data"
			drugs
			continue
		elseif data_type_string == "RPPA"
			i = 1
			all_proteins = df[[:Antibody_ID, :FullyValidated]]
			n = size(all_proteins, 1)
			for p_row in eachrow(all_proteins)
				is_fully_validated = p_row[:FullyValidated] == "Yes" ? true : false
				add_protein(experiment, p_row[:Antibody_ID], fully_validated = is_fully_validated)
			end
			i += 1
			if i%5000 == 0 log_message("processed $i of $n proteins...") end
		else
			i = 1
			n = size(df, 1)
			df_names = names(df)
			for data_row in eachrow(df)
				gene_object = nothing
				has_entrez_id = false
				has_hgnc_id = false
				has_ensembl_id = false
				entrez_id = -1
				ensembl_id = ""
				hgnc_id = ""

				# check for different id types and try to retrieve gene object
				if in(:EntrezID, df_names)
					has_entrez_id = true
					entrez_id = data_row[:EntrezID]
					try
						gene_object = get_gene(experiment, entrez_id)
					catch end
				end

				if in(:Ensembl_ID, df_names)
					has_ensembl_id = true
					ensembl_id = data_row[:Ensembl_ID]
					if gene_object == nothing
						try
							gene_object = get_gene(experiment, ensembl_id, "ensembl_id")
						catch end
					else
						add_gene_id!(experiment, gene_object, ensembl_id, "ensembl_id")
					end
				end

				if in(:HGNC_ID, df_names)
					has_hgnc_id = true
					hgnc_id = data_row[:HGNC_ID]
					if gene_object == nothing
						try
							gene_object = get_gene(experiment, hgnc_id, "hgnc_id")
						catch end
					else
						add_gene_id!(experiment, gene_object, hgnc_id, "hgnc_id")
					end
				end

				if has_entrez_id
					if gene_object == nothing
						#no gene found for any id type
						gene_object = add_gene(experiment, entrez_id)
						if has_ensembl_id add_gene_id!(experiment, gene_object, ensembl_id, "ensembl_id") end
						if has_hgnc_id add_gene_id!(experiment, gene_object, hgnc_id, "hgnc_id") end
					else
						# check if object had no entrez id before
						if gene_object.entrez_id == -1
							add_gene_id!(experiment, gene_object, entrez_id, "entrez_id")
						else
							# found gene object has another entrez id (which is unique) so create a new gene
							gene_object = add_gene(experiment, entrez_id)
							if has_ensembl_id add_gene_id!(experiment, gene_object, ensembl_id, "ensembl_id") end
							if has_hgnc_id add_gene_id!(experiment, gene_object, hgnc_id, "hgnc_id") end
						end
					end
				end

				if has_ensembl_id
					if gene_object == nothing
						# if still no gene object, no entrez id is present
						gene_object = add_gene(experiment, ensembl_id, "ensembl_id")
						if has_hgnc_id add_gene_id!(experiment, gene_object, hgnc_id, "hgnc_id") end
					else
						# ensembl ids also are uniqe
						if gene_object.ensembl_id == ""
							add_gene_id!(experiment, gene_object, ensembl_id, "ensembl_id")
						else
							gene_object = add_gene(experiment, ensembl_id, "ensembl_id")
							if has_entrez_id add_gene_id!(experiment, gene_object, entrez_id, "entrez_id") end
							if has_hgnc_id add_gene_id!(experiment, gene_object, hgnc_id, "hgnc_id") end
						end
					end
				end

				if has_hgnc_id
					if gene_object == nothing
						# if still no gene object, hgnc id is the only id present
						gene_object = add_gene(experiment, hgnc_id, "hgnc_id")
					else
						if gene_object.hgnc_id != hgnc_id add_gene_id!(experiment, gene_object, hgnc_id, "hgnc_id") end
					end
				end

				if in(Symbol("CancerGene?"), df_names)
					is_cancer_gene = data_row[Symbol("CancerGene?")] == "No" ? false : true
					gene_object.cancer_gene = is_cancer_gene
				end
				i+=1
				if i%5000 == 0 log_message("processed $i of $n genes...") end
			end
		end
		log_message("importing cell line data for $data_type")
		try
			for cl in get_cell_line_names_from_data_frame(data_type, df)
				cl_obj = get_cell_line!(experiment, cl, "Breast cancer")
				data_view = get_dataview!(cl_obj, data_type, DataView{key_type, data_type}(cl))
				populate_data_view!(data_view, df, experiment)
			end
			add_view!(experiment, data_type)
		catch exc
			rethrow(exc)
		end
	end
	log_message("number of genes now: Entrez=$(length(experiment.genes)), HGNC=$(length(experiment.genes_by_hgnc)) Ensembl=$(length(experiment.genes_by_ensembl))")
	# experiment
end

get_cell_line_names_from_data_frame(::Type{ExomeSeq}, df::DataFrame) = map(string, map(v->uppercase(v), unique(df[:CellLine])))
get_cell_line_names_from_data_frame(::Type{GeneExpression}, df::DataFrame) = extract_column_names(df, 2)
get_cell_line_names_from_data_frame(::Type{Methylation}, df::DataFrame) = extract_column_names(df, 5)
get_cell_line_names_from_data_frame(::Type{<:Union{RNASeq, RNASeqCall, RPPA, CNV}}, df::DataFrame) = extract_column_names(df, 3)

extract_column_names(df::DataFrame, first_used_column::Int64) = map(uppercase, map(string, names(df)[first_used_column:end]))

function populate_data_view!(data_view::DataView{Gene,ExomeSeq}, df::DataFrame, experiment::Experiment)
	cl_df = DataFrames.filter(r -> r[:CellLine] == data_view.cell_line_id, df)
	for data_row in eachrow(cl_df)
		gene = get_gene(experiment, data_row[:HGNC_ID])
		num_cosmic = data_row[Symbol("#Cosmic")]
		is_cancer_gene = data_row[Symbol("CancerGene?")] == "No" ? false : true
		variant_effect = data_row[:Type]
		protein_change = ismissing(data_row[:Summary]) ? "None" : data_row[:Summary]
		nucleotid_change = data_row[:AltBase]
		variant_confidence = Float64(data_row[:Confidence])
		norm_zygosity = data_row[Symbol("Zygosity(norm)")] == 0.0 ? "unk" : data_row[Symbol("Zygosity(norm)")]
		norm_reference_count = data_row[Symbol("RefCount(norm)")] == "." ? 0 : data_row[Symbol("RefCount(norm)")]
		norm_variant_count = data_row[Symbol("altCount(norm)")] == "." ? 0 : data_row[Symbol("altCount(norm)")]
		tumor_zygosity = data_row[Symbol("Zygosity(tumor)")]
		tumor_reference_count = Float64(data_row[Symbol("RefCount(tumor)")])
		tumor_variant_count = Float64(data_row[Symbol("AltCount(tumor)")])
		reference_mismatch_avg = data_row[Symbol("Avg#Mismatch(ref)")]
		variant_mismatch_avg = data_row[Symbol("Avg#Mismatch(alt)")]
		reference_mismatch_sum = data_row[Symbol("MismatchQualitySum(ref)")]
		variant_mismatch_sum = data_row[Symbol("MismatchQualitySum(alt)")]
		reference_dist3effective_avg = data_row[Symbol("DistanceEffective3'end(ref)")]
		variant_dist3effective_avg = data_row[Symbol("DistanceEffective3'end(alt)")]
		details = data_row[:Details]
		# data_summary = reference_mismatch_avg + variant_mismatch_avg + reference_mismatch_sum + variant_mismatch_sum + reference_dist3effective_avg + variant_dist3effective_avg

		#create some exome data view, create a key type for exome sequencing as a combination
		#of gene and protein change
		# exome_data = ExomeSeq(num_cosmic, variant_effect, protein_change,
		# 		nucleotid_change, variant_confidence,
		# 		norm_zygosity, norm_reference_count, norm_variant_count,
		# 		tumor_zygosity, tumor_reference_count, tumor_variant_count,
		# 		reference_mismatch_avg, variant_mismatch_avg,
		# 		reference_mismatch_sum, variant_mismatch_sum,
		# 		reference_dist3effective_avg, variant_dist3effective_avg,
		# 		details, data_summary)
		exome_data = ExomeSeq(protein_change, reference_mismatch_sum=reference_mismatch_sum, reference_mismatch_avg=reference_mismatch_avg,
								reference_dist3effective_avg=reference_dist3effective_avg, variant_mismatch_sum=variant_mismatch_sum,
								variant_mismatch_avg=variant_mismatch_avg, variant_dist3effective_avg=variant_dist3effective_avg,
								num_cosmic=num_cosmic, variant_effect=variant_effect, nucleotid_change=nucleotid_change, variant_confidence=variant_confidence,
								norm_zygosity=norm_zygosity, norm_reference_count=norm_reference_count, norm_variant_count=norm_variant_count,
								tumor_zygosity=tumor_zygosity, tumor_reference_count=tumor_reference_count, tumor_variant_count=tumor_variant_count,
								details=details)
		add_measurement!(data_view, gene, exome_data)
	end
end

function populate_data_view!(data_view::DataView{Gene,GeneExpression}, df::DataFrame, experiment::Experiment)
	cl_df = df[[:HGNC_ID, Symbol(data_view.cell_line_id)]]
	for data_row in eachrow(cl_df)
		if ismissing(data_row[2]) continue end
		gene = get_gene(experiment, data_row[1])
		gene_expression = GeneExpression(data_row[2])
		add_measurement!(data_view, gene, gene_expression)
	end
end

function populate_data_view!(data_view::DataView{Gene,Methylation}, df::DataFrame, experiment::Experiment)
	cl_df = df[[:HGNC_ID, Symbol(data_view.cell_line_id), :CGct1, :Cct1, :Illumina_ID]]
	for data_row in eachrow(cl_df)
		gene = get_gene(experiment, data_row[1])
		methylation = Methylation(data_row[2], data_row[3], data_row[4], illumina_id = data_row[5])
		add_measurement!(data_view, gene, methylation)
	end
end

function populate_data_view!(data_view::DataView{Gene,RNASeq}, df::DataFrame, experiment::Experiment)
	cl_df = df[[:Ensembl_ID, Symbol(data_view.cell_line_id)]]
	for data_row in eachrow(cl_df)
		gene = get_gene(experiment, data_row[1], "ensembl_id")
		add_measurement!(data_view, gene, RNASeq(data_row[2]))
	end
end

function populate_data_view!(data_view::DataView{Gene,RNASeqCall}, df::DataFrame, experiment::Experiment)
	cl_df = df[[:Ensembl_ID, Symbol(data_view.cell_line_id)]]
	for data_row in eachrow(cl_df)
		gene = get_gene(experiment, data_row[1], "ensembl_id")
		value = Bool(data_row[2])
		add_measurement!(data_view, gene, RNASeqCall(value))
	end
end

function populate_data_view!(data_view::DataView{Protein,RPPA}, df::DataFrame, experiment::Experiment)
	cl_df = df[[:Antibody_ID, Symbol(data_view.cell_line_id)]]
	for data_row in eachrow(cl_df)
		rppa = RPPA(data_row[2])
		protein = get_protein(experiment, data_row[1])
		add_measurement!(data_view, protein, rppa)
	end
end

function populate_data_view!(data_view::DataView{Gene,CNV}, df::DataFrame, experiment::Experiment)
	cl_df = df[[:EntrezID, Symbol(data_view.cell_line_id)]]
	for data_row in eachrow(cl_df)
		if ismissing(data_row[2]) continue end
		gene_level_cnv = CNV(data_row[2])
		gene = get_gene(experiment, data_row[1])
		add_measurement!(data_view, gene, gene_level_cnv)
	end
end


text:
"""
```julia
    SoftMaxLikelihood()
```
Multiclass likelihood with [Softmax transformation](https://en.wikipedia.org/wiki/Softmax_function):

```math
p(y=i|{fₖ}) = exp(fᵢ)/ ∑ₖexp(fₖ)
```

There is no possible augmentation for this likelihood
"""
struct SoftMaxLikelihood{T<:Real} <: MultiClassLikelihood{T}
    Y::AbstractVector{BitVector} #Mapping from instances to classes
    class_mapping::AbstractVector{Any} # Classes labels mapping
    ind_mapping::Dict{Any,Int} # Mapping from label to index
    y_class::AbstractVector{Int64} #GP Index for each sample
    θ::AbstractVector{AbstractVector{T}} # Variational parameter of Polya-Gamma distribution
    function SoftMaxLikelihood{T}() where {T<:Real}
        new{T}()
    end
    function SoftMaxLikelihood{T}(Y::AbstractVector{<:BitVector},
    class_mapping::AbstractVector, ind_mapping::Dict{<:Any,<:Int},y_class::AbstractVector{<:Int}) where {T<:Real}
        new{T}(Y,class_mapping,ind_mapping,y_class)
    end
end

function SoftMaxLikelihood()
    SoftMaxLikelihood{Float64}()
end

function pdf(l::SoftMaxLikelihood,f::AbstractVector)
    StatsFuns.softmax(f)
end

function pdf(l::SoftMaxLikelihood,y::Int,f::AbstractVector{<:Real})
    StatsFuns.softmax(f)[y]
end

function Base.show(io::IO,model::SoftMaxLikelihood{T}) where T
    print(io,"Softmax likelihood")
end

function init_likelihood(likelihood::SoftMaxLikelihood{T},nLatent::Integer,nSamplesUsed::Integer) where T
    if inference isa GibbsSampling
        θ = [abs.(rand(T,nSamplesUsed))*2 for i in 1:nLatent]
        LogisticSoftMaxLikelihood{T}(likelihood.Y,likelihood.class_mapping,likelihood.ind_mapping,likelihood.y_class,θ)
    else
        return likelihood
    end
end

function sample_local!(model::VGP{T,<:SoftMaxLikelihood,<:GibbsSampling}) where {T}
    pg = PolyaGammaDist()
    model.likelihood.θ .= broadcast((y::BitVector,γ::AbstractVector{<:Real},μ::AbstractVector{<:Real},i::Int64)->draw.([pg],1.0,μ-logsumexp()),model.likelihood.Y,model.likelihood.γ,model.μ,1:model.nLatent)
    return nothing #TODO FINISH AT SOME POINT
end

function grad_samples(model::AbstractGP{T,<:SoftMaxLikelihood},samples::AbstractMatrix{T},index::Integer) where {T}
    class = model.likelihood.y_class[index]
    grad_μ = zeros(T,model.nLatent)
    grad_Σ = zeros(T,model.nLatent)
    nSamples = size(samples,1)
    samples .= mapslices(StatsFuns.softmax,samples,dims=2)
    t = 0.0
    @inbounds for i in 1:nSamples
        s = samples[i,class]::T
        @views g_μ = grad_softmax(samples[i,:],class)/s
        grad_μ += g_μ
        @views h = diaghessian_softmax(samples[i,:],class)/s
        grad_Σ += h - abs2.(g_μ)
    end
    for k in 1:model.nLatent
        model.inference.∇μE[k][index] = grad_μ[k]/nSamples
        model.inference.∇ΣE[k][index] = 0.5.*grad_Σ[k]/nSamples
    end
end

function log_like_samples(model::AbstractGP{T,<:SoftMaxLikelihood},samples::AbstractMatrix,index::Integer) where {T}
    class = model.likelihood.y_class[index]
    nSamples = size(samples,1)
    loglike = mapslices(logsumexp,samples,dims=2)/nSamples
end

function grad_softmax(s::AbstractVector{<:Real},i::Integer)
    return (δ.(i,eachindex(s))-s)*s[i]
end

function diaghessian_softmax(s::AbstractVector{<:Real},i::Integer)
    return s[i]*(abs2.(δ.(i,eachindex(s))-s) - s.*(1.0.-s))
end

function hessian_softmax(s::AbstractVector{<:Real},i::Integer)
    m = length(s)
    hessian = zeros(m,m)
    for j in 1:m
        for k in 1:m
            hessian[j,k] = s[i]*((δ(i,k)-s[k])*(δ(i,j)-s[j])-s[j]*(δ(j,k)-s[k]))
        end
    end
    return hessian
end


text:
# -------------- #
# latex printing #
# -------------- #
const latex_subs = let
    repl_latex = REPL.REPLCompletions.latex_symbols
    unicode = Dict([(v, k) for (k, v) in repl_latex])
    tex_name = Dict([(strip(k, '\\'), k) for (k, v) in repl_latex])
    poppers = [k for k in keys(tex_name) if length(k)==1]
    for p in poppers
        pop!(tex_name, p)
    end
    merge(unicode, tex_name)
end

function _latex!(io::IO, v::Symbol, n::Union{Nothing,Integer}=nothing)
    # split the symbol into parts
    v_string = string(v)
    v_string = strip(v_string, '_')  # undo Dolang normalization
    front_ix = something(findfirst(_x -> _x == '_', v_string), 0)
    if front_ix == 0
        # no subscripts or superscripts
        front = v_string
        back = ""
    else
        front = v_string[1:front_ix-1]
        back = v_string[front_ix:end]
    end

    mods = []

    # figure how which of bar, star, dot we have
    while true
        m = match(r"(?<front>\w+)(?<mod>bar|star|dot)", front)
        # if we didn't match, we are done
        m === nothing && break
        front = m["front"]
        push!(mods, m["mod"])
    end

    # now we can get the variable name
    var = get(latex_subs, front, front)

    # add mods
    "dot" in mods && print(io, "\\dot{")
    "bar" in mods && print(io, "\\overline{")

    # now add variable name and close dot/bar braces
    print(io, var)
    "dot" in mods && print(io, "}")
    "bar" in mods && print(io, "}")

    # now work on subscripts
    _subs = collect(m.match for m in eachmatch(r"(?<!_)_([a-zA-Z0-9]+)", back))
    subscripts = strip.(_subs, ['_'])
    has_subscript = !isempty(subscripts)

    # now print subscripts
    (has_subscript || isa(n, Integer)) && print(io, "_{")
    if has_subscript
        for sub in subscripts
            _latex!(io, sub)
            print(io, ",")
        end
    end

    if isa(n, Integer)
        t = n > 0 ? print(io, "t+", n) :
        n < 0 ? print(io, "t-", -n) :
        print(io, "t")
    else
        # we have one too many commas here, delete one
        # io.size -= 1; io.ptr -= 1
    end

    (has_subscript || isa(n, Integer)) && print(io, "}")

    # now work on superscripts
    _sups = collect(m.match for m in eachmatch(r"__([a-zA-Z0-9]+)", back))
    superscripts = strip.(_sups, ['_'])
    has_superscript = !isempty(superscripts)
    (has_superscript || "star" in mods) &&  print(io, "^{")

    if "star" in mods
        print(io, "\\star")
        has_superscript && print(io, ",")
    end

    if has_superscript
        for sup in superscripts[1:end-1]
            _latex!(io, sup)
            print(io, ",")
        end
        _latex!(io, superscripts[end])
    end
    (has_superscript || "star" in mods) &&  print(io, "}")
end

_latex!(io::IO, s::AbstractString) = _latex!(io, Meta.parse(s))

_latex!(io::IO, x::Tuple{Symbol,T}) where {T<:Integer} = _latex!(io, x[1], x[2])

_latex!(io::IO, x::Number) = print(io, string(x))

function _latex!(io::IO, ex::Expr)
    if is_time_shift(ex)
        _latex!(io, ex.args[1], ex.args[2])
        return
    end

    if ex.head == :(=)
        _latex!(io, :($(ex.args[1]) == $(ex.args[2]) ))
        return
    end

    if ex.head == :block && length(ex.args) == 1  # quoted expression
        _latex!(io, ex.args[1])
        return
    end

    if ex.head == :call
        if ex.args[1] == :(==)
            _latex!(io, ex.args[2])
            print(io, " = ")
            _latex!(io, ex.args[3])
            return
        end
        f = ex.args[1]
        if f in ARITH_SYMBOLS
            if length(ex.args) == 2  # sometimes a unary negation
                print(io, f)
                _latex!(io, ex.args[2])
            elseif length(ex.args) == 3

                if f == :^
                    print(io, "\\left(")
                    _latex!(io, ex.args[2])
                    print(io, "\\right)")
                    print(io, f)
                    print(io, "{")
                    _latex!(io, ex.args[3])
                    print(io, "}")
                elseif f == :/
                    print(io, "\\frac{")
                    _latex!(io, ex.args[2])
                    print(io, "}{")
                    _latex!(io, ex.args[3])
                    print(io, "}")
                elseif f == :*
                    _latex!(io, ex.args[2])
                    print(io, " ")
                    _latex!(io, ex.args[3])
                elseif f == :-
                    print(io, "\\left(")
                    _latex!(io, ex.args[2])
                    print(io, f)
                    _latex!(io, ex.args[3])
                    print(io, "\\right)")
                else
                    _latex!(io, ex.args[2])
                    print(io, f)
                    _latex!(io, ex.args[3])
                end

            elseif length(ex.args) > 3  # chained + or *
                for arg in ex.args[2:end-1]
                    _latex!(io, arg)
                    print(io, f)
                end
                _latex!(io, ex.args[end])
            end
            return
        end

        # not arith
        print(io, "\\text{$f}\\left(")
        for arg in ex.args[2:end]
            _latex!(io, arg)
        end
        print(io, "\\right)")
        return
    end

    error("not sure what to do here")
end

function latex(ex::Union{Expr,Symbol})
    io = IOBuffer()
    _latex!(io, ex)
    String(take!(io))
end

latex(s::String) = latex(Meta.parse(s))

max_issues_repo_name: EconForge/Dolang


text:
# interfacing with other packages

using LinearAlgebra
using LinearAlgebra: BlasComplex, BlasFloat, BlasReal

function mv_wrapper(transa::SparseChar, alpha::Number, A::CuSparseMatrix{T}, X::DenseCuVector{T},
                    beta::Number, Y::CuVector{T}) where {T}
    mv!(transa, alpha, A, X, beta, Y, 'O')
end

function mm_wrapper(transa::SparseChar, transb::SparseChar, alpha::Number,
                    A::CuSparseMatrix{T}, B::CuMatrix{T}, beta::Number, C::CuMatrix{T}) where {T}
    if version() <= v"10.3.1"
        # Generic mm! doesn't support transposed B on CUDA10
        return mm2!(transa, transb, alpha, A, B, beta, C, 'O')
    end
    mm!(transa, transb, alpha, A, B, beta, C, 'O')
end

tag_wrappers = ((identity, identity),
                (T -> :(HermOrSym{T, <:$T}), A -> :(parent($A))))
op_wrappers = (
    (identity, 'N', identity),
    (T -> :(Transpose{<:T, <:$T}), 'T', A -> :(parent($A))),
    (T -> :(Adjoint{<:T, <:$T}), 'C', A -> :(parent($A)))
)
for (taga, untaga) in tag_wrappers, (wrapa, transa, unwrapa) in op_wrappers
    TypeA = wrapa(taga(:(CuSparseMatrix{T})))

    @eval begin
        function LinearAlgebra.mul!(C::CuVector{T}, A::$TypeA, B::DenseCuVector{T},
                                    alpha::Number, beta::Number) where {T <: BlasFloat}
            mv_wrapper($transa, alpha, $(untaga(unwrapa(:A))), B, beta, C)
        end
    end

    for (tagb, untagb) in tag_wrappers, (wrapb, transb, unwrapb) in op_wrappers
        TypeB = wrapb(tagb(:(DenseCuMatrix{T})))

        isadjoint(expr) = expr.head == :curly && expr.args[1] == :Adjoint
        if isadjoint(TypeA) || isadjoint(TypeB)
            # CUSPARSE defines adjoints only for complex inputs. For real inputs we run
            # adjoints as tranposes so that we can still support the whole API surface of
            # LinearAlgebra.
            @eval begin
                function LinearAlgebra.mul!(C::CuMatrix{T}, A::$TypeA, B::$TypeB,
                                            alpha::Number, beta::Number) where {T <: BlasComplex}
                    mm_wrapper($transa, $transb, alpha, $(untaga(unwrapa(:A))),
                               $(untagb(unwrapb(:B))), beta, C)
                end
            end

            transa_real = transa == 'C' ? 'T' : transa
            transb_real = transb == 'C' ? 'T' : transb
            @eval begin
                function LinearAlgebra.mul!(C::CuMatrix{T}, A::$TypeA, B::$TypeB,
                                            alpha::Number, beta::Number) where {T <: BlasReal}
                    mm_wrapper($transa_real, $transb_real, alpha, $(untaga(unwrapa(:A))),
                               $(untagb(unwrapb(:B))), beta, C)
                end
            end
        else
            @eval begin
                function LinearAlgebra.mul!(C::CuMatrix{T}, A::$TypeA, B::$TypeB,
                                            alpha::Number, beta::Number) where {T <: BlasFloat}
                    mm_wrapper($transa, $transb, alpha, $(untaga(unwrapa(:A))),
                               $(untagb(unwrapb(:B))), beta, C)
                end
            end
        end
    end
end

Base.:(\)(A::UpperTriangular{T, <:AbstractCuSparseMatrix{T}},
          B::DenseCuMatrix{T}) where {T<:BlasFloat} =
    sm('N',A,B,'O')
Base.:(\)(A::LowerTriangular{T, <:AbstractCuSparseMatrix{T}},
          B::DenseCuMatrix{T}) where {T<:BlasFloat} =
    sm('N',A,B,'O')
Base.:(\)(A::UpperTriangular{<:Any, <:Transpose{T, <:AbstractCuSparseMatrix{T}}},
          B::DenseCuMatrix{T}) where {T<:BlasFloat} =
    sm('T',parent(A),B,'O')
Base.:(\)(A::LowerTriangular{<:Any, <:Transpose{T, <:AbstractCuSparseMatrix{T}}},
          B::DenseCuMatrix{T}) where {T<:BlasFloat} =
    sm('T',parent(A),B,'O')
Base.:(\)(A::UpperTriangular{<:Any, <:Adjoint{T, <:AbstractCuSparseMatrix{T}}},
          B::DenseCuMatrix{T}) where {T<:BlasFloat} =
    sm('C',parent(A),B,'O')
Base.:(\)(A::LowerTriangular{<:Any, <:Adjoint{T, <:AbstractCuSparseMatrix{T}}},
          B::DenseCuMatrix{T}) where {T<:BlasFloat} =
    sm('C',parent(A),B,'O')

# TODO: some metaprogramming to reduce the amount of definitions here

Base.:(\)(A::UpperTriangular{T, <:AbstractCuSparseMatrix{T}},
          B::DenseCuVector{T}) where {T<:BlasFloat} =
    sv2('N', 'U', one(T), parent(A), B,'O')
Base.:(\)(A::LowerTriangular{T, <:AbstractCuSparseMatrix{T}},
          B::DenseCuVector{T}) where {T<:BlasFloat} =
    sv2('N', 'L', one(T), parent(A), B,'O')
Base.:(\)(A::UpperTriangular{<:Any, <:Transpose{T, <:AbstractCuSparseMatrix{T}}},
          B::DenseCuVector{T}) where {T<:BlasFloat} =
    sv2('T', 'L', one(T), parent(parent(A)), B, 'O')
Base.:(\)(A::LowerTriangular{<:Any, <:Transpose{T, <:AbstractCuSparseMatrix{T}}},
          B::DenseCuVector{T}) where {T<:BlasFloat} =
    sv2('T', 'U', one(T), parent(parent(A)), B, 'O')
Base.:(\)(A::UpperTriangular{<:Any, <:Adjoint{T, <:AbstractCuSparseMatrix{T}}},
          B::DenseCuVector{T}) where {T<:BlasFloat} =
    sv2('C', 'L', one(T), parent(parent(A)), B, 'O')
Base.:(\)(A::LowerTriangular{<:Any, <:Adjoint{T, <:AbstractCuSparseMatrix{T}}},
          B::DenseCuVector{T}) where {T<:BlasFloat} =
    sv2('C', 'U', one(T), parent(parent(A)), B, 'O')

Base.:(\)(A::UnitUpperTriangular{T, <:AbstractCuSparseMatrix{T}},
          B::DenseCuVector{T}) where {T<:BlasFloat} =
    sv2('N', 'U', one(T), parent(A), B, 'O', unit_diag=true)
Base.:(\)(A::UnitLowerTriangular{T, <:AbstractCuSparseMatrix{T}},
          B::DenseCuVector{T}) where {T<:BlasFloat} =
    sv2('N', 'L', one(T), parent(A), B, 'O', unit_diag=true)
Base.:(\)(A::UnitUpperTriangular{<:Any, <:Transpose{T, <:AbstractCuSparseMatrix{T}}},
          B::DenseCuVector{T}) where {T<:BlasFloat} =
    sv2('T', 'L', one(T), parent(parent(A)), B, 'O', unit_diag=true)
Base.:(\)(A::UnitLowerTriangular{<:Any, <:Transpose{T, <:AbstractCuSparseMatrix{T}}},
          B::DenseCuVector{T}) where {T<:BlasFloat} =
    sv2('T', 'U', one(T), parent(parent(A)), B, 'O', unit_diag=true)
Base.:(\)(A::UnitUpperTriangular{<:Any, <:Adjoint{T, <:AbstractCuSparseMatrix{T}}},
          B::DenseCuVector{T}) where {T<:BlasFloat} =
    sv2('C', 'L', one(T), parent(parent(A)), B, 'O', unit_diag=true)
Base.:(\)(A::UnitLowerTriangular{<:Any, <:Adjoint{T, <:AbstractCuSparseMatrix{T}}},
          B::DenseCuVector{T}) where {T<:BlasFloat} =
    sv2('C', 'U', one(T), parent(parent(A)), B, 'O', unit_diag=true)

Base.:(+)(A::Union{CuSparseMatrixCSR,CuSparseMatrixCSC},
          B::Union{CuSparseMatrixCSR,CuSparseMatrixCSC}) = geam(A,B,'O','O','O')
Base.:(-)(A::Union{CuSparseMatrixCSR,CuSparseMatrixCSC},
          B::Union{CuSparseMatrixCSR,CuSparseMatrixCSC}) = geam(A,-one(eltype(A)),B,'O','O','O')