text:
"""
    densityplot(x, y; kwargs...)

Description
============

Draws a density plot for the given points.

The first vector `x` should contain the horizontal
positions for all the points. The second vector `y` should then
contain the corresponding vertical positions respectively. This
means that the two vectors must be of the same length and
ordering.

Usage
======

    densityplot(x, y; name = "", title = "", xlabel = "", ylabel = "", labels = true, border = :solid, margin = 3, padding = 1, color = :auto, width = 40, height = 15, xlim = (0, 0), ylim = (0, 0), grid = false)

Arguments
==========

- **`x`** : The horizontal position for each point.
  If omitted, the axes of `y` will be used as `x`.

- **`y`** : The vertical position for each point.

- **`name`** : Annotation of the current drawing to be displayed
  on the right

$DOC_PLOT_PARAMS

- **`height`** : Number of character rows that should be used
  for plotting.

- **`xlim`** : Plotting range for the x axis.
  `(0, 0)` stands for automatic.

- **`ylim`** : Plotting range for the y axis.
  `(0, 0)` stands for automatic.

- **`grid`** : If `true`, draws grid-lines at the origin.

Returns
========

A plot object of type `Plot{DensityCanvas}`

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
=========

```julia-repl
julia> densityplot(randn(1000), randn(1000), title = "Density Plot")
                      Density Plot
        ┌────────────────────────────────────────┐
    2.9 │                    ░                   │
        │                ░    ░  ░               │
        │                ░░░  ░ ░ ░      ░       │
        │             ░░░ ░░▒░▒░░░ ▒ ░           │
        │          ░░ ░░ ░▒░░▒▒▒▒░░░░░░░         │
        │         ░░ ░▒░░░░░▓▒▒░▒░▒░▓░▒░░░       │
        │        ░ ░░░░▓▓▒▓▒▒█▒▓▒▒▓▒▒░▒ ░ ░      │
        │           ░ ▒▒▒▒▓▓▒▓▓▓▓▒█▒▒░░░░░░ ░    │
        │     ░   ░░░░░░░▒▒▓░░▓▒█▒▒▓▓▒▒ ░    ░   │
        │           ░ ░░▒▒░░▒░░▒░░░░ ░░  ░ ░     │
        │         ░   ░ ░ ░░░░░░ ░░ ░░ ░         │
        │              ░  ░ ░ ░  ░ ░   ░░░       │
        │                ░░  ░                   │
        │                                        │
   -3.3 │                                        │
        └────────────────────────────────────────┘
        -3.4                                   2.9
```

See also
=========

[`Plot`](@ref), [`scatterplot`](@ref), [`DensityCanvas`](@ref)
"""
function densityplot(
    x::AbstractVector,
    y::AbstractVector;
    color::UserColorType = :auto,
    grid = false,
    name = "",
    kw...,
)
    plot = Plot(x, y, DensityCanvas; grid = grid, kw...)
    scatterplot!(plot, x, y; color = color, name = name)
end

densityplot!(plot::Plot{<:DensityCanvas}, x::AbstractVector, y::AbstractVector; kw...) =
    scatterplot!(plot, x, y; kw...)


text:
todisplay = VERSION <= VersionNumber("1.2") ? true : false

using QHull
using Random
using SparseArrays
using LinearAlgebra
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using CAGD
if todisplay
	using ViewerGL
	GL = ViewerGL
	include("../examples/3d/views.jl")
end

V,CV = Lar.rod()([3,1])
ch = QHull.chull(convert(Lar.Points,V'))
FV = ch.simplices
EV = Lar.simplexFacets(FV)
GL.VIEW([ GL.GLGrid(V,EV,GL.COLORS[1],1) ]);

cyl = Lar.Struct([ Lar.t(0,0,0), (V,FV,EV) ]) #(0,0,-1.5)
cyl = Lar.struct2lar(cyl)

#V,FV = ([1.0 1.0 -0.5 -0.5 -0.5 -0.5; 0.0 0.0 0.86603 0.86603 -0.86603 -0.86603; 0.0 3.0 0.0 3.0 0.0 3.0],
#Array{Int64,1}[[3, 4, 1, 2], [3, 4, 5, 6], [5, 6, 1, 2], [3, 1, 5], [2, 4, 6]])
#EV = [[1,3],[3,5],[5,1], [2,4],[4,6],[6,2], [1,2],[3,4],[5,6]]
#cyl = (V,FV,EV)

tris = Lar.Struct([
    cyl,
    Lar.Struct([Lar.r(pi/2,0,0), cyl ]) ,
    Lar.Struct([Lar.r(0,pi/2,0), cyl ]) ,
    Lar.Struct([Lar.r(0,pi,0), cyl ]) ,
])
V,FV,EV = Lar.struct2lar(tris)

model = CAGD.Model(V)
scs = convert(Lar.ChainOp, Lar.coboundary_1(model.G, FV, EV))
CAGD.addModelCells!(model, 1, EV, signed = true)
CAGD.addModelCells!(model, 2, scs)

if todisplay  displayModel(model)  end

split_model = CAGD.pairwise_decomposition(model)

if todisplay  displayModel(split_model)  end

congr_model = CAGD.mergeModelVertices(split_model, signed_merge=true)

if todisplay  displayModel(congr_model)  end
gift_model = deepcopy(congr_model)
FCcycles, bicon_comps = CAGD.tgw(congr_model, 3, atol = atol)
CAGD.addModelCells!(gift_model, 3, convert(Lar.ChainOp, FC'))

if todisplay  viewExplode(gift_model, complete = false)  end

# Since first face is inner, then retrieving the external is needed

externals = [bc[CAGD.getExternalBoundary(congr_model, 3, FCcycles, bc)] for bc in bicon_comps]
internals = [setdiff(bicon_comps[i], externals[i]) for i = 1 : length(externals)]
assembled_model = deepcopy(congr_model)

FC = FCcycles[:, internals...]
JK = ones(Int8, length(externals))
FC = [FCcycles * SparseArrays.sparse(externals, JK, JK, size(FCcycles, 2), 1) FC]
CAGD.addModelCells!(assembled_model, 3, convert(Lar.ChainOp, FC'))

if todisplay  viewExplode(assembled_model)  end

ext: jl


text:
using Revise
using Random
using POMDPs
using LinearAlgebra
using POMDPModelTools
using POMDPPolicies
using POMDPSimulators
using IterTools
using AutomotiveDrivingModels
using AutomotivePOMDPs: polygon, overlap
using Parameters
using GridInterpolations
using StaticArrays
using AutoUrban
using AutoViz
using ProgressMeter
using BenchmarkTools
using Profile
using ProfileView
using Printf

includet("environment.jl")
includet("mdp_definition.jl")

const NTESTS = 1000

mdp = MergingMDP()
statespace = collect(states(mdp));
actionspace = actions(mdp)
rng = MersenneTwister(1)
test_states = rand(rng, statespace, NTESTS)

function benchmark_transition(mdp::MergingMDP, test_states::Vector{MergingState}, actionspace::UnitRange{Int64})
    for s in test_states
        for a in actionspace
            d = transition(mdp, s, a)
        end
    end
end

println("Starting transition benchmark")

@btime benchmark_transition($mdp, $test_states, $actionspace)

Profile.clear()
@profile benchmark_transition(mdp, test_states, actionspace)

ProfileView.view()


text:
export spConnect

function spConnect(model;pre=[],post=[],expr="",p=0.1)
  var=match(r"^.*(?==)",expr).match
  value=match(r"(?<==).*",expr).match

  if !contains(model.connection,var)
    model.connection*="@parallel for p in workers() @eval $var=spzeros(N,N_tot) end"
  end

  for preTemp in [pre;],postTemp in [post;]
    i=model.groups[postTemp][:localindex]
    j=model.groups[preTemp][:localindex]
    processor=model.groups[postTemp][:processor]
    model.connection*="@spawnat $(processor+1) $var[:,$j]=$value*sprand($(length(i)),$(length(j)),$p)\n"
  end
end # end of function spConnect


text:
using BinDeps

@BinDeps.setup

nettle = library_dependency("nettle", aliases = ["libnettle","libnettle-4-6"])

@windows_only begin
  using WinRPM
  provides(WinRPM.RPM, "libnettle-4-6", nettle, os = :Windows )
end

@osx_only begin
  using Homebrew
  provides( Homebrew.HB, "nettle", nettle, os = :Darwin )
end

provides( AptGet, "libnettle4", nettle )
provides( Yum, "nettle", nettle )

julia_usrdir = normpath(JULIA_HOME*"/../") # This is a stopgap, we need a better builtin solution to get the included libraries
libdirs = String["$(julia_usrdir)/lib"]
includedirs = String["$(julia_usrdir)/include"]
env = {"HOGWEED_LIBS" => "-L$(libdirs[1]) -L$(BinDeps.libdir(nettle)) -lhogweed -lgmp",
       "NETTLE_LIBS" => "-L$(libdirs[1]) -L$(BinDeps.libdir(nettle)) -lnettle -lgmp",
       "LD_LIBRARY_PATH" => join([libdirs[1];BinDeps.libdir(nettle)],":")}

provides( Sources,
          URI("http://www.lysator.liu.se/~nisse/archive/nettle-2.7.1.tar.gz"),
          SHA="bc71ebd43435537d767799e414fce88e521b7278d48c860651216e1fc6555b40",
          nettle )
provides( BuildProcess,
          Autotools(lib_dirs = libdirs,
                    include_dirs = includedirs,
                    env = env,
                    configure_options = ["--disable-openssl", "--libdir=$(BinDeps.libdir(nettle))"]),
          nettle )

@BinDeps.install [:nettle => :nettle]


text:
"""
    parse_kw(ps::ParseState)

Dispatch function for when the parser has reached a keyword.
"""
function parse_kw(ps::ParseState)
    k = kindof(ps.t)
    if k === Tokens.IF
        return @default ps @closer ps :block parse_if(ps)
    elseif k === Tokens.LET
        return @default ps @closer ps :block parse_blockexpr(ps, Let)
    elseif k === Tokens.TRY
        return @default ps @closer ps :block parse_try(ps)
    elseif k === Tokens.FUNCTION
        return @default ps @closer ps :block parse_blockexpr(ps, FunctionDef)
    elseif k === Tokens.MACRO
        return @default ps @closer ps :block parse_blockexpr(ps, Macro)
    elseif k === Tokens.BEGIN
        @static if VERSION < v"1.4"
            return @default ps @closer ps :block parse_blockexpr(ps, Begin)
        else
            if ps.closer.inref
                ret = mKEYWORD(ps)
            else
                return @default ps @closer ps :block parse_blockexpr(ps, Begin)
            end
        end
    elseif k === Tokens.QUOTE
        return @default ps @closer ps :block parse_blockexpr(ps, Quote)
    elseif k === Tokens.FOR
        return @default ps @closer ps :block parse_blockexpr(ps, For)
    elseif k === Tokens.WHILE
        return @default ps @closer ps :block parse_blockexpr(ps, While)
    elseif k === Tokens.BREAK
        return INSTANCE(ps)
    elseif k === Tokens.CONTINUE
        return INSTANCE(ps)
    elseif k === Tokens.IMPORT
        return parse_imports(ps)
    elseif k === Tokens.USING
        return parse_imports(ps)
    elseif k === Tokens.EXPORT
        return parse_export(ps)
    elseif k === Tokens.MODULE
        return @default ps @closer ps :block parse_blockexpr(ps, ModuleH)
    elseif k === Tokens.BAREMODULE
        return @default ps @closer ps :block parse_blockexpr(ps, BareModule)
    elseif k === Tokens.CONST
        return @default ps parse_const(ps)
    elseif k === Tokens.GLOBAL
        return @default ps parse_global(ps)
    elseif k === Tokens.LOCAL
        return @default ps parse_local(ps)
    elseif k === Tokens.RETURN
        return @default ps parse_return(ps)
    elseif k === Tokens.END
        if ps.closer.square
            ret = mKEYWORD(ps)
        else
            ret = mErrorToken(ps, mIDENTIFIER(ps), UnexpectedToken)
        end
        return ret
    elseif k === Tokens.ELSE || k === Tokens.ELSEIF || k === Tokens.CATCH || k === Tokens.FINALLY
        return mErrorToken(ps, mIDENTIFIER(ps), UnexpectedToken)
    elseif k === Tokens.ABSTRACT
        return @default ps parse_abstract(ps)
    elseif k === Tokens.PRIMITIVE
        return @default ps parse_primitive(ps)
    elseif k === Tokens.TYPE
        return mIDENTIFIER(ps)
    elseif k === Tokens.STRUCT
        return @default ps @closer ps :block parse_blockexpr(ps, Struct)
    elseif k === Tokens.MUTABLE
        return @default ps @closer ps :block parse_mutable(ps)
    elseif k === Tokens.OUTER
        return mIDENTIFIER(ps)
    else
        return mErrorToken(ps, Unknown)
    end
end

function parse_const(ps::ParseState)
    kw = mKEYWORD(ps)
    arg = parse_expression(ps)
    if !(is_assignment(unwrapbracket(arg)) || (typof(arg) === Global && is_assignment(unwrapbracket(arg.args[2]))))
        arg = mErrorToken(ps, arg, ExpectedAssignment)
    end
    ret = EXPR(Const, EXPR[kw, arg])
    return ret
end

function parse_global(ps::ParseState)
    kw = mKEYWORD(ps)
    arg = parse_expression(ps)

    return EXPR(Global, EXPR[kw, arg])
end

function parse_local(ps::ParseState)
    kw = mKEYWORD(ps)
    arg = parse_expression(ps)

    return EXPR(Local, EXPR[kw, arg])
end

function parse_return(ps::ParseState)
    kw = mKEYWORD(ps)
    # Note to self: Nothing could be treated as implicit and added
    # during conversion to Expr.
    args = closer(ps) ? NOTHING() : parse_expression(ps)

    return EXPR(Return, EXPR[kw, args])
end

function parse_abstract(ps::ParseState)
    if kindof(ps.nt) === Tokens.TYPE
        kw1 = mKEYWORD(ps)
        kw2 = mKEYWORD(next(ps))
        sig = @closer ps :block parse_expression(ps)
        ret = EXPR(Abstract, EXPR[kw1, kw2, sig, accept_end(ps)])
    else
        ret = mIDENTIFIER(ps)
    end
    return ret
end

function parse_primitive(ps::ParseState)
    if kindof(ps.nt) === Tokens.TYPE
        kw1 = mKEYWORD(ps)
        kw2 = mKEYWORD(next(ps))
        sig = @closer ps :ws @closer ps :wsop parse_expression(ps)
        arg = @closer ps :block parse_expression(ps)
        ret = EXPR(Primitive, EXPR[kw1, kw2, sig, arg, accept_end(ps)])
    else
        ret = mIDENTIFIER(ps)
    end
    return ret
end

function parse_mutable(ps::ParseState)
    if kindof(ps.nt) === Tokens.STRUCT
        kw = mKEYWORD(ps)
        next(ps)
        ret = parse_blockexpr(ps, Mutable)
        pushfirst!(ret, kw)
        update_span!(ret)
    else
        ret = mIDENTIFIER(ps)
    end
    return ret
end

function parse_imports(ps::ParseState)
    kw = mKEYWORD(ps)
    kwt = is_import(kw) ? Import : Using

    arg = parse_dot_mod(ps)

    if !iscomma(ps.nt) && !iscolon(ps.nt)
        ret = EXPR(kwt, vcat(kw, arg))
    elseif iscolon(ps.nt)
        ret = EXPR(kwt, vcat(kw, arg))
        push!(ret, mOPERATOR(next(ps)))

        arg = parse_dot_mod(ps, true)
        append!(ret, arg)
        prevpos = position(ps)
        while iscomma(ps.nt)
            accept_comma(ps, ret)
            arg = parse_dot_mod(ps, true)
            append!(ret, arg)
            prevpos = loop_check(ps, prevpos)
        end
    else
        ret = EXPR(kwt, vcat(kw, arg))
        prevpos = position(ps)
        while iscomma(ps.nt)
            accept_comma(ps, ret)
            arg = parse_dot_mod(ps)
            append!(ret, arg)
            prevpos = loop_check(ps, prevpos)
        end
    end

    return ret
end

function parse_export(ps::ParseState)
    args = EXPR[mKEYWORD(ps)]
    append!(args, parse_dot_mod(ps))

    prevpos = position(ps)
    while iscomma(ps.nt)
        push!(args, mPUNCTUATION(next(ps)))
        arg = parse_dot_mod(ps)[1]
        push!(args, arg)
        prevpos = loop_check(ps, prevpos)
    end

    return EXPR(Export, args)
end

"""
    parse_blockexpr_sig(ps::ParseState, head)

Utility function to parse the signature of a block statement (i.e. any statement preceding
the main body of the block). Returns `nothing` in some cases (e.g. `begin end`)
"""
function parse_blockexpr_sig(ps::ParseState, head)
    if head === Struct || head == Mutable || head === While
        return @closer ps :ws parse_expression(ps)
    elseif head === For
        return parse_iterators(ps)
    elseif head === FunctionDef || head === Macro
        sig = @closer ps :inwhere @closer ps :ws parse_expression(ps)
        if convertsigtotuple(sig)
            sig = EXPR(TupleH, sig.args)
        end
        prevpos = position(ps)
        while kindof(ps.nt) === Tokens.WHERE && kindof(ps.ws) != Tokens.NEWLINE_WS
            sig = @closer ps :inwhere @closer ps :ws parse_operator_where(ps, sig, INSTANCE(next(ps)), false)
            prevpos = loop_check(ps, prevpos)
        end
        return sig
    elseif head === Let
        if isendoflinews(ps.ws)
            return nothing
        else
            arg = @closer ps :comma @closer ps :ws  parse_expression(ps)
            if iscomma(ps.nt) || !(is_wrapped_assignment(arg) || isidentifier(arg))
                arg = EXPR(Block, EXPR[arg])
                prevpos = position(ps)
                while iscomma(ps.nt)
                    accept_comma(ps, arg)
                    startbyte = ps.nt.startbyte
                    nextarg = @closer ps :comma @closer ps :ws parse_expression(ps)
                    push!(arg, nextarg)
                    prevpos = loop_check(ps, prevpos)
                end
            end
            return arg
        end
    elseif head === Do
        sig = EXPR(TupleH, EXPR[])
        prevpos = position(ps)
        @closer ps :comma @closer ps :block while !closer(ps)
            @closer ps :ws a = parse_expression(ps)
            push!(sig, a)
            if kindof(ps.nt) === Tokens.COMMA
                accept_comma(ps, sig)
            elseif @closer ps :ws closer(ps)
                break
            end
            prevpos = loop_check(ps, prevpos)
        end
        return sig
    elseif head === ModuleH || head === BareModule
        return isidentifier(ps.nt) ? mIDENTIFIER(next(ps)) :
            @precedence ps 15 @closer ps :ws parse_expression(ps)
    end
    return nothing
end

function parse_do(ps::ParseState, pre::EXPR)
    ret = parse_blockexpr(next(ps), Do)
    pushfirst!(ret, pre)
    update_span!(ret)
    return ret
end

"""
    parse_blockexpr(ps::ParseState, head)

General function for parsing block expressions comprised of a series of statements
terminated by an `end`.
"""
function parse_blockexpr(ps::ParseState, head)
    kw = mKEYWORD(ps)
    sig = parse_blockexpr_sig(ps, head)
    blockargs = parse_block(ps, EXPR[], (Tokens.END,), docable(head))

    if sig === nothing
        EXPR(head, EXPR[kw, EXPR(Block, blockargs), accept_end(ps)])
    elseif (head === FunctionDef || head === Macro) && is_either_id_op_interp(sig)
        EXPR(head, EXPR[kw, sig, accept_end(ps)])
    else
        EXPR(head, EXPR[kw, sig, EXPR(Block, blockargs), accept_end(ps)])
    end
end

"""
    parse_if(ps, ret, nested=false, puncs=[])

Parse an `if` block.
"""
function parse_if(ps::ParseState, nested=false)
    # Parsing
    kw = mKEYWORD(ps)
    if isendoflinews(ps.ws)
        cond = mErrorToken(ps, MissingConditional)
    else
        cond = @closer ps :ws parse_expression(ps)
    end
    ifblockargs = parse_block(ps, EXPR[], (Tokens.END, Tokens.ELSE, Tokens.ELSEIF))

    if nested
        ret = EXPR(If, EXPR[cond, EXPR(Block, ifblockargs)])
    else
        ret = EXPR(If, EXPR[kw, cond, EXPR(Block, ifblockargs)])
    end

    elseblockargs = EXPR[]
    if kindof(ps.nt) === Tokens.ELSEIF
        push!(ret, mKEYWORD(next(ps)))
        push!(elseblockargs, parse_if(ps, true))
    end
    elsekw = kindof(ps.nt) === Tokens.ELSE
    if kindof(ps.nt) === Tokens.ELSE
        push!(ret, mKEYWORD(next(ps)))
        parse_block(ps, elseblockargs)
    end

    # Construction
    if !(isempty(elseblockargs) && !elsekw)
        push!(ret, EXPR(Block, elseblockargs))
    end
    !nested && accept_end(ps, ret)

    return ret
end

function parse_try(ps::ParseState)
    kw = mKEYWORD(ps)
    ret = EXPR(Try, EXPR[kw])

    tryblockargs = parse_block(ps, EXPR[], (Tokens.END, Tokens.CATCH, Tokens.FINALLY))
    push!(ret, EXPR(Block, tryblockargs))

    #  catch block
    if kindof(ps.nt) === Tokens.CATCH
        next(ps)
        push!(ret, mKEYWORD(ps))
        # catch closing early
        if kindof(ps.nt) === Tokens.FINALLY || kindof(ps.nt) === Tokens.END
            caught = FALSE()
            catchblock = EXPR(Block, EXPR[])
        else
            if isendoflinews(ps.ws)
                caught = FALSE()
            else
                caught = @closer ps :ws parse_expression(ps)
            end

            catchblockargs = parse_block(ps, EXPR[], (Tokens.END, Tokens.FINALLY))
            if !(is_either_id_op_interp(caught) || kindof(caught) === Tokens.FALSE)
                pushfirst!(catchblockargs, caught)
                caught = FALSE()
            end
            catchblock = EXPR(Block, catchblockargs)
        end
    else
        caught = FALSE()
        catchblock = EXPR(Block, EXPR[])
    end
    push!(ret, caught)
    push!(ret, catchblock)

    # finally block
    if kindof(ps.nt) === Tokens.FINALLY
        if isempty(catchblock.args)
            ret.args[4] = setparent!(FALSE(), ret)
        end
        push!(ret, mKEYWORD(next(ps)))
        finallyblockargs = parse_block(ps)
        push!(ret, EXPR(Block, finallyblockargs))
    end

    push!(ret, accept_end(ps))
    return ret
end


text:
################################################################################
#
#          EllCrv/Torsion.jl : Computing torsion points on elliptic curves
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# (C) 2016 Tommy Hofmann
# (C) 2016 Robin Ammon
# (C) 2016 Sofia Brenner
# (C) 2022 Jeroen Hanselman
#
################################################################################

export is_torsion_point, torsion_points, torsion_structure, torsion_bound, pr_torsion_basis,
       division_polynomial, division_polynomial_univariate, torsion_points_division_poly, order,
       torsion_points_lutz_nagell

################################################################################
#
#  Order of a point
#
################################################################################

@doc Markdown.doc"""
    order(P::EllCrvPt{nf_elem}) -> fmpz

Returns the order of the point $P$ or $0$ if the order is infinite.
"""
function order(P::EllCrvPt{T}) where T<:Union{nf_elem, fmpq}
  if T==fmpq
    N = 12
  else
    E = parent(P)
    N = torsion_bound(E, 20)
  end

  Q = P
  for i in 1:N
    if !is_finite(Q)
      return fmpz(i)
    end
    Q = Q + P
  end

  return fmpz(0)
end

################################################################################
#
#  Test if point is torsion point
#
################################################################################

@doc Markdown.doc"""
    is_torsion_point(P::EllCrvPt{fmpq}) -> Bool

Returns whether the point $P$ is a torsion point.
"""
function is_torsion_point(P::EllCrvPt{T}) where T <: Union{nf_elem,fmpq}
  o = order(P)
  return o != 0
end

@doc Markdown.doc"""
    is_torsion_point(P::EllCrvPt{fq_nmod}) -> Bool

Returns whether the point $P$ is a torsion point.
"""
function is_torsion_point(P::EllCrvPt{fq_nmod})
  return true
end

################################################################################
#
#  Torsion points over QQ
#
################################################################################

# via theorem of Lutz-Nagell
@doc Markdown.doc"""
    torsion_points_lutz_nagell(E::EllCrv{fmpq}) -> Vector{EllCrvPt}

Computes the rational torsion points of an elliptic curve using the
>Lutz-Nagell theorem.
"""
function torsion_points_lutz_nagell(F::EllCrv{fmpq})

  if F.short == false
    (G, trafo, ruecktrafo) = short_weierstrass_model(F)
  else
    G = F
  end

  # transform the curve to an equivalent one with integer coefficients, if necessary
  (E, trafo_int, trafo_rat) = integral_model(G)

  res = [infinity(E)]
  d = discriminant(E)

  # Lutz-Nagell: necessary: y = 0 or y^2 divides d

  ycand = collect(squaredivisors(numerator(d))) # candidates for y-coordinate

  push!(ycand,0)

  pcand = Tuple{fmpz, fmpz}[] # candidates for torsion points

  Zx, x = PolynomialRing(FlintZZ, "x")

  _, _, _, a4, a6 = a_invars(E)

  # Lutz-Nagell: coordinates of torsion points need to be in ZZ
  for i = 1:length(ycand)
    # are there corresponding integer x-values?
    xcand = zeros(x^3 + numerator(a4)*x + numerator(a6) - ycand[i]^2)
    if length(xcand) != 0
      for j = 1: length(xcand)
        push!(pcand, (xcand[j], ycand[i])) # add to candidates
      end
    end
  end

  # check if candidates are torsion points
  if length(pcand) != 0
    for i = 1:length(pcand)
      P = E([pcand[i][1],pcand[i][2]])

      # Mazur: Order of torsion point is at most 12
      Q = P
      for j = 2:12
        Q = Q + P
        if is_infinite(Q)
          push!(res, P)
          break
        end
      end
    end
  end

  torsionpoints = res

  #if F.short == false
  #  for i = 1:length(torsionpoints)
  #    torsionpoints[i]= ruecktrafo(trafo_rat(torsionpoints[i]))
  #  end
  #else
  #  for i = 1:length(torsionpoints)
  #    torsionpoints[i] = trafo_rat(torsionpoints[i])
  #  end
  #end
  return torsionpoints
end

# via division polynomials
@doc Markdown.doc"""
    torsion_points_division_poly(E::EllCrv{fmpq}) -> Array{EllCrvPt}

Computes the rational torsion points of a rational elliptic curve $E$ using
division polynomials.
"""
function torsion_points_division_poly(F::EllCrv{fmpq})

  # if necessary, transform curve into short form
  if F.short == false
    (G, trafo, ruecktrafo) = short_weierstrass_model(F)
  else
    G = F
  end

  # transform the curve to an equivalent one with integer coefficients, if necessary
  (E, trafo_int, trafo_rat) = integral_model(G)

  # check whether we already know the torsion points or not
  if isdefined(E, :torsion_points)
    if F.short == false
      F.torsion_points = EllCrvPt{fmpq}[]
      for i = 1:length(E.torsion_points)
        push!(F.torsion_points, ruecktrafo(trafo_rat(E.torsion_points[i])))
      end
    else
      F.torsion_points = EllCrvPt{fmpq}[]
      for i = 1:length(E.torsion_points)
        push!(F.torsion_points, trafo_rat(E.torsion_points[i]))
      end
    end
    return F.torsion_points::Vector{EllCrvPt{fmpq}}
  end

  # curve has integer coefficients
  _, _, _, A, B = map(numerator, a_invars(E))

  torsionpoints = [infinity(E)]

  # points of order 2 (point has order 2 iff y-coordinate is zero)
  # (note: these points are not detected by the division polynomials)

  Zx, x = PolynomialRing(FlintZZ, "x")

  s = zeros(x^3 + A*x + B) # solutions of x^3 + Ax + B = 0
  if length(s) != 0
    for i = 1:length(s)
      P = E([s[i], 0])
      push!(torsionpoints, P)
    end
  end

  # Mazur: order of torsion point is at most 12
  for n = 7 : 12 # points of smaller order will also be found
    Psi = division_polynomial_univariate(E,n)[1]
    z = roots(Psi)

    if length(z) == 0
      continue
    end

    # z contains candidates for x-coordinates of torsion points
    for i = 1:length(z)

      # are there corresponding integer-square y-coordinates?
      if isinteger(z[i])
        zi = numerator(z[i])
        ysquarecand = zi^3 + A*zi + B

        if ysquarecand < 0
          continue
        end

        a = isqrtrem(ysquarecand)

        if a[1] == 0 # then ysquarecand = 0, so have found torsion point
          P = E([zi, ysquarecand])

        # if not already contained in torsionpoints, add P
          if !(P in torsionpoints)
            push!(torsionpoints, P)
          end

          continue
        end

      # now ysquarecand > 0
        if a[2] == 0 # then y is a square of an integer
          P = E([zi, a[1]])
          Q = E([zi, -a[1]])

        # if not already contained in torsionpoints, add P and Q
          if !(P in torsionpoints)
            push!(torsionpoints, P)
            push!(torsionpoints, Q)
          end
        end
      end
    end
  end

  if F.short == false
    for i = 1:length(torsionpoints)
      torsionpoints[i]= ruecktrafo(trafo_rat(torsionpoints[i]))
    end
  else
    for i = 1:length(torsionpoints)
      torsionpoints[i] = trafo_rat(torsionpoints[i])
    end
  end

  return torsionpoints::Vector{EllCrvPt{fmpq}}
end

# function for users
@doc Markdown.doc"""
    torsion_points(E::EllCrv{fmpq}) -> Vector{EllCrvPt{fmpq}}

Returns the rational torsion points of $E$.
"""
function torsion_points(E::EllCrv{fmpq})
  if isdefined(E, :torsion_points)
    return E.torsion_points::Vector{EllCrvPt{fmpq}}
  end

  t = torsion_points_division_poly(E::EllCrv{fmpq})
  E.torsion_points = t
  return t
end

################################################################################
#
#  Torsion structure
#
################################################################################

@doc Markdown.doc"""
    torsion_structure(E::EllCrv{fmpq}) -> (A::Vector{fmpz},
                                           B::Vector{EllCrvPt{fmpq}}

Compute the structure of the rational torsion group of an elliptic curve $E$.
Then `A` is an array with `A = [n]` resp. `A = [n,m]` such that the
rational torsion of $E$ is isomorphic to $\mathbf Z/n\mathbf Z$ resp.
$\mathbf Z/n\mathbf Z \times \mathbf Z/m\mathbf Z$.
And `B` is an array of points with `B = [P]` and $P$ has order $n$ resp.
`B = [P, Q]` and $P$ has order $n$, $Q$ has order $m$.
"""
function torsion_structure(E::EllCrv{fmpq})
  T = torsion_points(E)
  grouporder = length(T)
  orders = fmpz[]

  for i in 1:grouporder
    push!(orders, 0)
  end

 # determine orders of group elements
  for i in 1:grouporder
    j = 1
    while (j < 13) && (orders[i] == 0)
      if (j*T[i]).is_infinite == true
        orders[i] = j
      end
    j = j + 1
    end
  end

  # find generators
  if in(grouporder, orders) == true # is the group cyclic?
    k = something(findfirst(isequal(grouporder), orders), 0)
    return fmpz[grouporder], [T[k]]
  else # group not cyclic
    m = div(grouporder, 2)
    k1 = something(findfirst(isequal(2), orders), 0)
    k2 = something(findlast(isequal(m), orders), 0) # findlast to get different points if m = 2
    points = [T[k1], T[k2]]
  return (fmpz[2, m], points)
  end
end

################################################################################
#
#  Torsion over more general fields
#
################################################################################

@doc Markdown.doc"""
    torsion_bound(E::EllCrv{nf_elem}, n::Int) -> fmpz

Bound the order of the torsion subgroup of $E by considering
the order of the reduction of $E$ modulo $n$ discinct primes
with good reduction
"""
function torsion_bound(E::EllCrv{nf_elem}, n::Int)

  R = ring_of_integers(base_field(E))
  badp = bad_primes(E)
  p = next_prime(1)
  i = 0
  bound = 0
  while i < n
    L = prime_ideals_over(R, p)
    p = next_prime(p)
    for P in L
      if !(P in badp)
        i=i+1
        EP = modp_reduction(E,P)
        bound = gcd(bound, order(EP))
      end
    end
  end
  return(bound)
end

################################################################################
#
#  p^r-Torsion part
#
################################################################################

#Adapted from Sage: ell_generic.py
@doc Markdown.doc"""
    pr_torsion_basis(E::EllCrv{nf_elem}, p::fmpz, r = Int) -> Vector{EllCrvPt}

Compute a basis for E[p^r]. When r is given the algorithm stops searching after
having found p^r points.
"""
function pr_torsion_basis(E::EllCrv, p, r = typemax(Int))

  if !isprime(p)
    throw(DomainError(p, "p should be a prime number"))
  end

  #First we find all the p-torsion points
  p_torsion = division_points(infinity(E), p)
  p_rank = Int(log(fmpz(p), fmpz(length(p_torsion))))

  if p_rank == 0
    return Tuple{elem_type(E), Int}[]
  elseif p_rank == 1
  #If the dimension of the p-torsion is 1.
    P = p_torsion[1]
    if is_infinite(P)
      P=p_torsion[2]
    end
    k = 1
    if r==1
      return [(P,k)]
    end

    #We keep dividing P by p until we have found a generator for the p^r-torsion.
    pts = division_points(P, p)

    while length(pts) > 0
      k += 1
      P = pts[1]
      if r <= k
        return [(P,k)]
      end
      points = division_points(P, p)
    end
    return [(P,k)]
  else  #The p-torsion has rank 2
    P1 = popfirst!(p_torsion)
    while is_infinite(P1)
      P1 = popfirst!(p_torsion)
    end
    P2 = popfirst!(p_torsion)
    #We find a linearly independent basis for the p-torsion.
    while linearly_dependent(P1,P2)
     P2 = popfirst!(p_torsion)
    end

    k = 1
    log_order = 2
    if r<= 2
      return [(P1,1),(P2,1)]
    end

    #We keep dividing P1 and P2 by p until this is no longer possible
    pts1 = division_points(P1, p)
    pts2 = division_points(P2, p)
    while length(pts1) > 0 && length(pts2) > 0
      k += 1
      P1 = pts1[1]
      P2 = pts2[1]
      log_order += 2
      if r<=log_order
        return [(P1,k),(P2,k)]
      end
      pts1 = division_points(P1, p)
      pts2 = division_points(P2, p)
    end
  #Now k is the maximal integer for which P1,P2 are a basis
  #for the p^k torsion

  #But there could still be a p^k torsion point which can be divided
  #further, so we search for a linear combination of P1 and P2 to
  #see if this is the case.

    if length(pts1) > 0
      pts = pts1
    elseif length(pts2) > 0
      P1, P2 = P2, P1
      pts = pts2
    else
      for Q in [P1+a*P2 : a in (1:p-1)]
          # Q runs through P1+a*P2 for a=1,2,...,p-1
        pts = division_points(Q, p)
        if len(pts) > 0
          P1 = Q
          break
        end
      end
    end
    if len(pts)==0
      return [(P1,k),(P2,k)]
    end
  #If we have found a P1 that we can divide further,
  #we continue trying to divide P1 by p. If we fail
  #we replace P1 again by a linear combination P1+a*P2 for
  #some a in [1..p-1]. We continue with this until we can no longer
  #divide any linear combination and we have found the full torsion structure
  #Z/p^n x Z/p^k for some n >k.

    n = k
    while true
      P1 = pts[1]
      n += 1
      log_order += 1
      if r <= log_order
        return [(P1,n),(P2,k)]
      end

      pts = division_points(P1, p)
      if len(pts) == 0
        for Q in [P1+a*P2: a in (1:p-1)]
        # Q runs through P1+a*P2 for a=1,2,...,p-1
          pts = division_points(Q, p)
          if len(pts) > 0
            break
          end
        end
        if len(pts) == 0
          return [(P1,n),(P2,k)]
        end
      end
    end
  end
end

#Returns [m, n] with m >n. This is inconsistent with the way torsion_structure
#returns elements for EllCrv over QQ.
@doc Markdown.doc"""
    torsion_structure(E::EllCrv{nf_elem}) -> (A::Vector{fmpz},
                                           B::Vector{EllCrvPt{nf_elem}}

Compute the structure of the rational torsion group of an elliptic curve $E$.
Then `A` is an array with `A = [n]` resp. `A = [n,m]` such that the
rational torsion of $E$ is isomorphic to $\mathbf Z/n\mathbf Z$ resp.
$\mathbf Z/n\mathbf Z \times \mathbf Z/m\mathbf Z$.
And `B` is an array of points with `B = [P]` and $P$ has order $n$ resp.
`B = [P, Q]` and $P$ has order $n$, $Q$ has order $m$.
"""
function torsion_structure(E::EllCrv{nf_elem})

  T1 = T2 = infinity(E)
  k1 = k2 = 1
  bound = torsion_bound(E, 20)
  for (p,r) in factor(bound)
    ptor = pr_torsion_basis(E, p, r)
    if length(ptor) > 0
      T1 += ptor[1][1]
      k1 *= p^(ptor[1][2])
    end
    if length(ptor) > 1
      T2 += ptor[2][1]
      k2 *= p^(ptor[2][2])
    end
  end

  if k1 == 1
    structure = [fmpz(1)]
    gens = [infinity(E)]
  elseif k2 == 1
    structure = [fmpz(k1)]
    gens = [T1]
  else
    structure = [fmpz(k1), fmpz(k2)]
    gens = [T1, T2]
  end
  return structure, gens
end

################################################################################
#
#  Misc
#
################################################################################

function linearly_dependent(P1::EllCrvPt{T}, P2::EllCrvPt{T}) where T

#if !istorsion_point(P1)||!istorsion_point(P2)
#  throw(DomainError, "Points need to be torsion_points")
#end

  A = P1
  B = P2

  while is_finite(A)
    while is_finite(B)
      if is_infinite(A+B)
        return true
      end
      B+=P2
    end
    A+=P1
    B+=P2
  end
  return false
end

################################################################################
#
#  Division polynomials
#
################################################################################
"""
    division_polynomial_univariate(E::EllCrv, n::Int, [x]) -> Poly

Compute the n-th univariate division polynomial of an elliptic curve defined
over a field k following Mazur and Tate. By default the result is a univariate polynomial over the base ring of `E`.
When `x` is given, the output is evaluated using the given value for `x`.

A triple of objects is returned:
- The n-th division polynomial as a univariate polynomial with a mutiplicity
   of 2 at the non-zero two-torsion points.
- The n-th division polynomial as a univariate polynomial divided by the
  univariate 2-torsion polynomial when n is even.
- The complementary factor, i.e. the first output divided by the second output.
"""
function division_polynomial_univariate(E::EllCrv, n::S, x = PolynomialRing(base_field(E),"x")[2]) where S<:Union{Integer, fmpz}

  R = parent(x)
  if is_short_weierstrass_model(E)
    poly = divpol_g_short(E,n,x)
    if mod(n,2) == 0
      _, _, _, A, B = a_invars(E)
      twotorsfactor = 4*(x^3+A*x+B)
    else
      twotorsfactor = one(R)
    end
  else
    poly = divpol_g(E,n,x)
      if mod(n,2) == 0
        b2, b4, b6 = b_invars(E)
        twotorsfactor = 4*x^3+b2*x^2+2*b4*x+b6
      else
        twotorsfactor = one(R)
      end
  end
  return twotorsfactor*poly, poly, twotorsfactor
end

"""
    division_polynomial(E::EllCrv, n::Int, x, y) -> Poly

Compute the n-th division polynomial of an elliptic curve defined over a field
k following Mazur and Tate. When x and or y are given the output is
automatically evaluated using the given values.
"""
function division_polynomial(E::EllCrv, n::S, x = PolynomialRing(base_field(E),"x")[2], y = PolynomialRing(parent(x),"y")[2]) where S<:Union{Integer, fmpz}
  R = parent(y)
  if is_short_weierstrass_model(E)
     if mod(n,2) == 0
      return 2*y*divpol_g_short(E,n,x)
    else
      return R(divpol_g_short(E,n,x))
    end
  else
    a1, _, a3 = a_invars(E)
    if mod(n,2) == 0
      return (2*y + a1*x + a3)*divpol_g(E,n,x)
    else
      return R(divpol_g(E,n,x))
    end
  end
end

function divpol_g_short(E::EllCrv, n::S, x = PolynomialRing(base_field(E),"x")[2]) where S<:Union{Integer, fmpz}

  Kx = parent(x)
  _, _, _, A, B = a_invars(E)

  B6sqr = (4*x^3+4*A*x+4*B)^2

  if n == 1 || n == 2
    return one(parent(x))
  elseif n == 3
    return 3*x^4 + 6*(A)*x^2 + 12*(B)*x - (A)^2
  elseif n == 4
    return 2*(x^6 + 5*(A)*x^4 + 20*(B)*x^3 - 5*(A)^2*x^2 - 4*(A)*(B)*x - 8*(B)^2 - (A)^3)
  elseif mod(n,2) == 0
    m = div(n,2)
    return (divpol_g_short(E,m,x))*(divpol_g_short(E,m+2, x)*divpol_g_short(E,m-1,x)^2 - divpol_g_short(E,m-2,x)*divpol_g_short(E,m+1,x)^2)
  else m = div(n-1,2)
    m = div(n-1,2)
    part1 = divpol_g_short(E,m+2,x)  * divpol_g_short(E,m,x)^3
    part2 = divpol_g_short(E,m-1,x) * divpol_g_short(E,m+1,x)^3
    if mod(m,2) == 0
      return B6sqr * part1 - part2
    else
      return part1 - B6sqr * part2
    end
  end
end

function divpol_g(E::EllCrv, n::S, x = PolynomialRing(base_field(E),"x")[2]) where S<:Union{Integer, fmpz}

  Kx = parent(x)

  b2, b4, b6, b8 = E.b_invars
  B4 = 6*x^2+b2*x+b4
  B6sqr = (4*x^3+b2*x^2+2*b4*x+b6)^2
  B8 = 3*x^4 + b2*x^3 + 3*b4*x^2 + 3*b6*x + b8

  if n == 1 || n ==2
    return one(Kx)
  elseif n == 3
    return B8
  elseif n == 4
    return -B6sqr+B4*B8
  elseif mod(n,2) == 0
    m = div(n-2,2)
    return divpol_g(E,m+3,x)*divpol_g(E,m+1,x)*divpol_g(E,m,x)^2 - divpol_g(E,m+1,x)*divpol_g(E,m-1,x)*divpol_g(E,m+2,x)^2
  else
    m = div(n-1,2)
    part1 = divpol_g(E,m+2,x)  * divpol_g(E,m,x)^3
    part2 = divpol_g(E,m-1,x) * divpol_g(E,m+1,x)^3
    if mod(m,2) == 0
      return B6sqr * part1 - part2
    else
      return part1 - B6sqr * part2
    end
  end
end


text:
# Physical parameter monitoring.

"""
    check_plasma_characteristics(n, v, T, B)

Display characteristic plasma parameters given density `n`, bulk velocity `v`, temperature
`T`, and magnetic field strength `B` in SI units.
"""
function check_plasma_characteristics(n, v, T, B)
   γ = 5/3    # adiabatic index

   p = n*kB*T # [Pa]

   # Characteristic parameters
   di  = √(mᵢ*ϵ₀/n)*c/qᵢ           # ion inertial length, [m]
   fci = qᵢ * B / (mᵢ * 2π)        # [/s]
   vA  = B / √(μ₀ * n * mᵢ)        # Alfvén speed, [m/s]
   vS  = √(γ * p / (n * mᵢ))       # sonic speed, [m/s]
   vT  = √(kB * T / mᵢ)            # thermal speed, [m/s]
   rᵢ  = mᵢ * vT / (qᵢ * B)        # gyroradius, [m]

   println("--------------------------------------------------")
   println("* Characteristic plasma properties")
   println("Density             : ", rpad(round(n/1e6; digits=2), 9), "amu/cc")
   println("Velocity            : ", rpad(round(v/1e3; digits=2), 9), "km/s")
   println("Pressure            : ", rpad(round(p*1e9; digits=3), 9), "nPa")
   println("Temperature         : ", rpad(round(T; digits=3), 9), "K")
   println("Magnetic field      : ", rpad(round(B*1e9; digits=2), 9), "nT")
   println("Ion inertial length : ", rpad(round(di/1e3; digits=2),9), "km")
   println("Ion gyroradius      : ", rpad(round(rᵢ/1e3; digits=2),9), "km")
   println("Gyrofrequency       : ", rpad(round(fci; digits=2),   9), "Hz")
   println("Alfvén speed        : ", rpad(round(vA/1e3; digits=2),9), "km/s")
   println("Sonic speed         : ", rpad(round(vS/1e3; digits=2),9), "km/s")
   println("--------------------------------------------------")
end


text:
#############################################################################
# trace_logm(X, C) returns tr(C*logm(X)) where X and C are a positive definite
# matrices and C is constant.
#
# trace_logm is concave in X.
#
# This function implements the semidefinite programming approximation given
# in the reference below.  Parameters m and k control the accuracy of this
# approximation: m is the number of quadrature nodes to use and k the number
# of square-roots to take. See reference for more details.
#
# Implementation uses the expression
#   tr(C*logm(X)) = -tr(C*D_{op}(I||X))
# where D_{op} is the operator relative entropy:
#   D_{op}(X||Y) = X^{1/2}*logm(X^{1/2} Y^{-1} X^{1/2})*X^{1/2}
#
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#
#REFERENCE
#   Ported from CVXQUAD which is based on the paper: "Lieb's concavity
#   theorem, matrix geometric means and semidefinite optimization" by Hamza
#   Fawzi and James Saunderson (arXiv:1512.03401)
#############################################################################

struct TraceLogm <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr,}
    size::Tuple{Int, Int}
    C::AbstractMatrix
    m::Integer
    k::Integer

    function TraceLogm(X::AbstractExpr, C::AbstractMatrix, m::Integer, k::Integer)
        children = (X,)
        if size(X) != size(C)
            throw(DimensionMismatch("X and C must be the same size"))
        end
        n = size(X)[1]
        if size(X) != (n, n)
            throw(DimensionMismatch("X and C must be square"))
        end
        if norm(C - C') > 1e-6
            throw(DomainError(C, "C must be Hermitian"))
        end
        if any(eigvals(Hermitian(C)) .< -1e-6)
            throw(DomainError(C, "C must be positive semidefinite"))
        end
        return new(:trace_logm, hash(children), children, (1, 1), C, m, k)
    end
end

function sign(atom::TraceLogm)
    return NoSign()
end

function monotonicity(atom::TraceLogm)
    return (NoMonotonicity(),)
end

function curvature(atom::TraceLogm)
    return ConcaveVexity()
end

function evaluate(atom::TraceLogm)
    X = evaluate(atom.children[1])
    return trace_logm(X, atom.C)
end

const MatrixOrConstant = Union{AbstractMatrix, Constant}

function trace_logm(X::AbstractExpr, C::MatrixOrConstant, m::Integer=3, k::Integer=3)
    #println("trace_logm general case")
    return TraceLogm(X, evaluate(C), m, k)
end

function trace_logm(X::MatrixOrConstant, C::MatrixOrConstant, m::Integer=3, k::Integer=3)
    #println("trace_logm constant X")
    eye = Matrix(1.0*I, size(X))
    return -quantum_relative_entropy(C, X) + quantum_relative_entropy(C, eye);
end

function conic_form!(atom::TraceLogm, unique_conic_forms)
    if !has_conic_form(unique_conic_forms, atom)
        X = atom.children[1]
        C = atom.C
        m = atom.m
        k = atom.k
        eye = Matrix(1.0*I, size(X))

        is_complex = sign(X) == ComplexSign() || sign(Constant(C)) == ComplexSign()
        if is_complex
            τ = ComplexVariable(size(X))
        else
            τ = Variable(size(X))
        end
        conic_form!(τ in RelativeEntropyEpiCone(eye, X, m, k), unique_conic_forms)

        # It's already a real mathematically, but need to make it a real type.
        t = real(-tr(C * τ))
        cache_conic_form!(unique_conic_forms, atom, maximize(t))
    end
    return get_conic_form(unique_conic_forms, atom)
end


text:
"""
    For loading networks with stoichiometry stored in matrices.
    Assumed that the substrate and product stoichiometry matrices
    are stored as numspecies by numrxs matrices, with entry (i,j)
    giving the stoichiometric coefficient of species i within rx j.
"""

"""
    For dense matrices
"""
function loadrxnetwork(ft::MatrixNetwork, networkname::String,
                        rateexprs::AbstractVector,
                        substoich::AbstractMatrix,
                        prodstoich::AbstractMatrix;
                        species=Symbol[],
                        params=Symbol[])

    sz = size(substoich)
    @assert sz == size(prodstoich)
    numspecs = sz[1]
    numrxs = sz[2]

    # create the network
    rn = eval(Meta.parse("@empty_reaction_network $networkname"))

    # create the species if none passed in
    isempty(species) && (species = [Symbol("S",i) for i=1:numspecs])
    foreach(i -> addspecies!(rn, species[i]), 1:numspecs)

    # create the parameters
    foreach(param -> addparam!(rn, param), params)

    # create the reactions
    subs = Vector{Pair{Symbol,Int}}()
    prods = Vector{Pair{Symbol,Int}}()
    for j = 1:numrxs
        empty!(subs)
        empty!(prods)

        # stoich
        for i = 1:numspecs
            scoef = substoich[i,j]
            (scoef > 0) && push!(subs, species[i] => scoef)

            pcoef = prodstoich[i,j]
            (pcoef > 0) && push!(prods, species[i] => pcoef)
        end

        addreaction!(rn, rateexprs[j], Tuple(subs), Tuple(prods))
    end

    ParsedReactionNetwork(rn, nothing)
end

"""
    For sparse matrices
"""
function loadrxnetwork(ft::MatrixNetwork, networkname::String,
                        rateexprs::AbstractVector,
                        substoich::SparseMatrixCSC,
                        prodstoich::SparseMatrixCSC;
                        species=Symbol[],
                        params=Symbol[])

    sz = size(substoich)
    @assert sz == size(prodstoich)
    numspecs = sz[1]
    numrxs = sz[2]

    # create the network
    rn = eval(Meta.parse("@empty_reaction_network $networkname"))

    # create the species if none passed in
    isempty(species) && (species = [Symbol("S",i) for i=1:numspecs])
    foreach(i -> addspecies!(rn, species[i]), 1:numspecs)

    # create the parameters
    foreach(param -> addparam!(rn, param), params)

    # create the reactions
    subs = Vector{Pair{Symbol,Int}}()
    prods = Vector{Pair{Symbol,Int}}()
    srows = rowvals(substoich)
    svals = nonzeros(substoich)
    prows = rowvals(prodstoich)
    pvals = nonzeros(prodstoich)
    for j = 1:numrxs
        empty!(subs)
        empty!(prods)

        for ir in nzrange(substoich, j)
           i     = srows[ir]
           scoef = svals[ir]
           (scoef > 0) && push!(subs, species[i] => scoef)
        end

        for ir in nzrange(prodstoich, j)
            i     = prows[ir]
            pcoef = pvals[ir]
            (pcoef > 0) && push!(prods, species[i] => pcoef)
         end

         addreaction!(rn, rateexprs[j], Tuple(subs), Tuple(prods))
     end

    ParsedReactionNetwork(rn, nothing)
end


text:
"""
    test_project_extras(package::Union{Module, PkgId})
    test_project_extras(packages::Vector{Union{Module, PkgId}})

Check that test target of the root project and test project
(`test/Project.toml`) are consistent.  This is useful for supporting
Julia < 1.2 while recording test-only dependency compatibility in
`test/Project.toml`.
"""
function test_project_extras(packages)
    @testset "$(result.label)" for result in analyze_project_extras(packages)
        @debug result.label result
        @test result ⊜ true
    end
end

function project_toml_path(dir)
    candidates = joinpath.(dir, ["Project.toml", "JuliaProject.toml"])
    i = findfirst(isfile, candidates)
    i === nothing && return candidates[1], false
    return candidates[i], true
end

analyze_project_extras(packages) = map(_analyze_project_extras, aspkgids(packages))

is_julia12_or_later(compat::AbstractString) = is_julia12_or_later(semver_spec(compat))
is_julia12_or_later(compat::VersionSpec) = isempty(compat ∩ semver_spec("1.0 - 1.1"))

function _analyze_project_extras(pkg::PkgId)
    label = string(pkg)

    result = root_project_or_failed_lazytest(pkg)
    result isa LazyTestResult && return result
    root_project_path = result

    pkgpath = dirname(dirname(Base.locate_package(pkg)))
    test_project_path, found = project_toml_path(joinpath(pkgpath, "test"))
    if !found
        return LazyTestResult(label, "test/Project.toml file does not exist.", true)
    end
    root_project = TOML.parsefile(root_project_path)
    test_project = TOML.parsefile(test_project_path)

    # Ignore root project's extras if only supporting julia 1.2 or later.
    # See: # https://julialang.github.io/Pkg.jl/v1/creating-packages/#Test-specific-dependencies-in-Julia-1.2-and-above-1
    julia_version = get(get(root_project, "compat", Dict()), "julia", "1")
    if is_julia12_or_later(julia_version)
        return LazyTestResult(
            label,
            string(
                "Supporting only post-1.2 `julia` ($julia_version); ",
                "ignoring root project.",
            ),
            true,
        )
    end

    # `extras_deps`: test-only dependencies according to /Project.toml
    all_extras_deps = get(root_project, "extras", Dict())
    target = Set{String}(get(get(root_project, "targets", Dict()), "test", []))
    extras_deps = setdiff(
        Set{Pair{String,String}}(p for p in all_extras_deps if first(p) in target),
        Set{Pair{String,String}}(get(root_project, "deps", [])),
    )

    # `test_deps`: test-only dependencies according to /test/Project.toml:
    test_deps = setdiff(
        Set{Pair{String,String}}(get(test_project, "deps", [])),
        Set{Pair{String,String}}(get(root_project, "deps", [])),
        [root_project["name"] => root_project["uuid"]],
    )

    not_in_extras = setdiff(test_deps, extras_deps)
    not_in_test = setdiff(extras_deps, test_deps)
    if isempty(not_in_extras) && isempty(not_in_test)
        return LazyTestResult(
            label,
            """
            Root and test projects are consistent.
            Root project: $root_project_path
            Test project: $test_project_path
            """,
            true,
        )
    else
        msg = sprint() do io
            println(io, "Root and test projects are inconsistent.")
            if !isempty(not_in_extras)
                println(io, "Test dependencies not in root project ($root_project_path):")
                for (name, uuid) in sort!(collect(not_in_extras))
                    println(io, "    $name = \"$uuid\"")
                end
            end
            if !isempty(not_in_test)
                println(io, "Dependencies not in test project ($test_project_path):")
                for (name, uuid) in sort!(collect(not_in_test))
                    println(io, "    $name = \"$uuid\"")
                end
            end
        end
        return LazyTestResult(label, msg, false)
    end
end


text:
const _pm_task_fields = ("type", "measurements_dict",)

const _index_types = Union{Int64,String}
const _name_types  = String

const _accepted_task_types = (
    "Input_Handler",
    "Output_Handler",
    "Processor",
    "Memory",
)

build_taskgraph(c::MapConstructor) = build_taskgraph(c, parse(c))
function build_taskgraph(c::MapConstructor, json_dict::Dict)
    taskgraph = Taskgraph()
    options = json_dict[_options_path_]
    parse_input(taskgraph, json_dict["task_structure"], options)

    # post-processing routines
    for op in taskgraph_ops(c)
        @debug "Running Taskgraph Transform $op"
        taskgraph = op(taskgraph, options)
    end
    return taskgraph
end

taskgraph_ops(::PMConstructor) = (
    transform_task_types,
    compute_edge_metadata,
    apply_link_weights,
)

function parse_input(t::Taskgraph, tasklist, options)
    # Check if packet_links will be included in the netlist.
    use_packet_links = options[:use_packet_links]

    # One iteration to collect all of the task nodes.
    for task in tasklist
        name = task["name"]
        metadata = getkeys(task, _pm_task_fields)
        add_node(t, TaskgraphNode(name, metadata = metadata))
    end

    # Iterate through again, collect and add all edges. Each edge can be
    # uniquely represented as a tuple
    # (link_class, source_task, source_index, dest_task, dest_index).
    #
    # Keep track of these tuples to avoid adding redundant links.
    seen_tuples = Set()
    packet_warn = true
    for task in tasklist
        # leaf_node_dict should be a dict of dicts. Where:
        # - keys are link class ("Circuit_Link", "Packet_Link" etc)
        # - values eventually end up with the link tuple.
        leaf_node_dict = task["leaf_node_dict"]
        for (link_class, link_dict) in leaf_node_dict
            if !use_packet_links && link_class == "Packet_Link"
                continue
            end
            # Deep dictionary walking ...
            for io_decl in values(link_dict), link_def in values(io_decl)
                # Extract the link tuple fields.
                source_task  = type_sanitize(_name_types,link_def["source_task"])
                source_index = type_sanitize(_index_types, link_def["source_index"])
                dest_task    = type_sanitize(_name_types,link_def["dest_task"])
                dest_index   = type_sanitize(_index_types, link_def["dest_index"])
                # Construct and check link tuple
                link_tuple = (link_class, source_task, source_index, dest_task, dest_index)
                in(link_tuple, seen_tuples) && continue
                push!(seen_tuples, link_tuple)

                # Determine whether or not this link should be routed.
                # Right now, only Packet Links are not supposed to be routed.
                # TODO: Fix this.
                route_link = !(link_class == "Packet_Link")

                # Create metadata for this link.
                metadata = Dict{String,Any}(
                      "pm_class"            => link_class,
                      "source_task"         => source_task,
                      "source_index"        => source_index,
                      "dest_task"           => dest_task,
                      "dest_index"          => dest_index,
                      "route_link"          => route_link,
                      "measurements_dict"   => link_def["measurements_dict"],
                     )

                new_edge = TaskgraphEdge(source_task, dest_task, metadata = metadata)
                add_edge(t, new_edge)
            end
        end
    end
    @info """
    Taskgraph Constructed.

    Number of Tasks: $(num_nodes(t))

    Number of Links: $(num_edges(t))
    """
end

#-------------------------------------------------------------------------------
# Project manager taskgraph transforms.
#-------------------------------------------------------------------------------
function compute_edge_metadata(t::Taskgraph, options::Dict)
    for edge in getedges(t)
        edge.metadata["link_class"] = lowercase(edge.metadata["pm_class"])
        # Preserve the destination fifo if the edge is a circuit link and it
        # does not end at an output handler.
        if edge.metadata["pm_class"] == "Circuit_Link"
            neighbor = getsinks(edge)
            @assert length(neighbor) == 1

            neighbor_node = getnode(t, first(neighbor))

            # I came across an edge case when doing so mappings for ASAP2 with
            # the output mux. For Asap3/4, there is only one output port, so
            # setting "preserve_dest" has no affect. But for Asap2, it's nice
            # to not have this set to the router can use any of the output
            # ports.
            if isoutput(neighbor_node)
                edge.metadata["preserve_dest"] = false
            else
                edge.metadata["preserve_dest"] = true
            end
        else
            edge.metadata["preserve_dest"] = false
        end
    end
    return t
end

"""
    transform_task_types(t::Taskgraph, options::Dict)

Transform the Project Manager Task Types to the task type used internally by
the Mapper.
"""
function transform_task_types(t::Taskgraph, options::Dict)
    memory_neighbors = String[]
    for task in getnodes(t)
        # Get the project manager class
        task_type = task.metadata["type"]

        # Brute-force decoding.
        if task_type == "Input_Handler"
            make_input!(task)
        elseif task_type == "Output_Handler"
            make_output!(task)
        elseif task_type == "Processor"
            make_proc!(task)
        elseif task_type == "Memory"
            # join the input and output to determine if 1 port or 2
            neighbors = union(innode_names(t, task.name), outnode_names(t, task.name))
            num_neighbors = length(unique(neighbors))

            if num_neighbors in (1,2)
                make_memory!(task, num_neighbors)
            else
                error("""
                    Expected memory $(task.name) to have 1 or 2 neighbors.
                    Found $num_neighbors.
                    """)
            end
            append!(memory_neighbors, neighbors)
        end
    end

    # Give "memory_processor" attributes to all neighbors of memory processors.
    for name in unique(memory_neighbors)
        task = getnode(t, name)
        if !isproc(task) && !ismemoryproc(task)
            error("Expected neighbors of memories to have type \"processor\".")
        end
        make_memoryproc!(task)
    end
    return t
end

function apply_link_weights(t::Taskgraph, options::Dict)
    # Number of binary digits to round the weights to. Making this number too
    # larger results in some links being ignored entirely, which is pretty bad
    # for routing purposes. Setting it to 3 seems to hit a sweet spot.
    ndigits = 3

    # Set a minimum linke weight to ensure that no link gets assigned a weight
    # of 0. This causes problems for both routing and placement.
    minimum_link_weight = 2.0 ^ (-ndigits)

    # Weight to assign memory links.
    memory_link_weight = 5.0

    # Need to still apply weight in order to get memory modules to work correctly.
    # Just assigns unit weights to each non-memory link and a weight of 5.0
    # to each memory link.
    use_profiled_links = options[:use_profiled_links]

    # Determine the maximium and minimum number of writes over the entire
    # colletion of edges. Do this in a single pass through the input
    # dictionary as this will work correctly even if the "num_writes" field
    # of the measurements dict is not available.
    min_writes = typemax(Int64)
    max_writes = typemin(Int64)
    for edge in getedges(t)
        measurements = edge.metadata["measurements_dict"]
        source_task = getnode(t, first(getsources(edge)))
        if haskey(measurements, "num_writes")
            num_writes = measurements["num_writes"]
            # Keep track of minimum and maximum number of writes.
            min_writes = min(min_writes, num_writes)
            max_writes = max(max_writes, num_writes)
        end
    end

    range = max_writes - min_writes

    @debug """
        Min Writes: $min_writes
        Max Writes: $max_writes
        Minimum link weight: $minimum_link_weight
    """

    # Assign weight to each edge based on the number of writes compared with
    # the average.
    for edge in getedges(t)
        measurements = edge.metadata["measurements_dict"]
        source_task = getnode(t, first(getsources(edge)))
        dest_task   = getnode(t, first(getsinks(edge)))

        # Since memory processors must be located directly next to their
        # respective memories, assign a high cost to memory links.
        if edge.metadata["pm_class"] in ("Memory_Request_Link", "Memory_Response_Link")
                edge.metadata["cost"] = memory_link_weight

        # Default value if per-link weight is not being used, or if for some
        # reason measurements associated with a link are missing.
        elseif !use_profiled_links || !haskey(measurements, "num_writes")
            edge.metadata["cost"] = 1.0

        # By default, scale the weight of a link with the number of writes
        # on the link. Scale between the minimum and maximum number of writes
        # so the weight is between "minimum_link_weight" and 1.0
        else
            num_writes = measurements["num_writes"]
            scaled_cost = (num_writes - min_writes) / range
            cost = max(
                       ceil(scaled_cost, digits = ndigits, base = 2),
                       minimum_link_weight
                      )

            edge.metadata["cost"] = cost
        end
    end
    return t
end


text:
using SossDistributions
using Test

@testset "SossDistributions.jl" begin
    # Write your own tests here.
end


text:
#
# Factors Reduce Tests
#

@testset "Factors Redude" begin

@testset "main" begin
    x = Dimension(:X, [2, 5, 7])
    y = Dimension(:Y, [2, 3])
    z = Dimension(:Z, ['a', 'c'])

    ϕ = Factor([x, y, z], [1, 2, 3, 2, 3, 4, 4, 6, 7, 8, 10, 16])

    df_original = DataFrame(ϕ)

    @test_throws ArgumentError reducedim(*, ϕ, :waldo)
    # make sure it didn't change ϕ
    @test DataFrame(ϕ) == df_original

    @test sum(broadcast(*, ϕ, :Z, 0)).potential == reshape([0.0], ())

    df = DataFrame(X = [2, 5, 7], potential = [123.0, 165.0, 237.0])

    ϕ2 = broadcast(*, ϕ, :Z, [1, 10])
    ϕ2 = sum(ϕ2, [:Y, :Z])

    # ϕ didn't change
    @test DataFrame(ϕ2) != df_original
    @test DataFrame(ϕ2) == df

    df2 = DataFrame(X = [2, 5, 7], potential = [15, 21, 30])
    @test DataFrame(sum(ϕ, [:Y, :Z])) == df2

    @test first(Z(ϕ).potential) == sum(ϕ.potential)
end

@testset "singleton" begin
    ϕ = Factor(1:16, :X=>4, :Z=>4)
    z = Z(ϕ)
    @test ndims(z) == 0
    @test z.potential == reshape([136], ())
    @test names(z) == Symbol[]

    ϕ2 = sum(ϕ, :Z)
    @test ϕ2.potential == [28, 32, 36, 40]
    @test_approx_eq prod(ϕ2).potential 1290240.0

    ϕs = Factor(23)
    @test minimum(ϕs).potential == reshape([23], ())
end

end


text:
@testset "Dataset splitting helpers" begin

    @testset "slice computation" begin
        @test slicesof([11], 2) == [(1, 1, 1, 6), (1, 7, 1, 11)]
        a = [30, 40, 50]
        @test slicesof(a, 1) == [(1, 1, 3, 50)]

        @test slicesof(a, 2) == [(1, 1, 2, 30), (2, 31, 3, 50)]

        @test slicesof(a, 3) == [(1, 1, 2, 10), (2, 11, 3, 10), (3, 11, 3, 50)]

        @test slicesof(a, 4) ==
              [(1, 1, 1, 30), (2, 1, 2, 30), (2, 31, 3, 20), (3, 21, 3, 50)]

        @test slicesof(a, 5) ==
              [(1, 1, 1, 24), (1, 25, 2, 18), (2, 19, 3, 2), (3, 3, 3, 26), (3, 27, 3, 50)]

        @test slicesof(a, 10) == [
            (1, 1, 1, 12),
            (1, 13, 1, 24),
            (1, 25, 2, 6),
            (2, 7, 2, 18),
            (2, 19, 2, 30),
            (2, 31, 3, 2),
            (3, 3, 3, 14),
            (3, 15, 3, 26),
            (3, 27, 3, 38),
            (3, 39, 3, 50),
        ]

    end

    @testset "slice collection" begin
        s = slicesof([4, 4], 3)

        @test vcollectSlice(i -> repeat([i], 4), s[1])[:, 1] == [1, 1, 1]
        @test vcollectSlice(i -> repeat([i], 4), s[2])[:, 1] == [1, 2, 2]
        @test vcollectSlice(i -> repeat([i], 4), s[3])[:, 1] == [2, 2]
    end
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    SingleBridgeOptimizer{BT<:AbstractBridge}(model::MOI.ModelLike)

Return [`MOI.Bridges.AbstractBridgeOptimizer`](@ref) that always bridges any
variables constrained on creation supported by the bridge `BT`.

This is in contrast with the [`MOI.Bridges.LazyBridgeOptimizer`](@ref),
which only bridges the variables constrained on creation if they are supported
by the bridge `BT` and unsupported by `model`.

!!! warning
    Two `SingleBridgeOptimizer`s cannot be used together as both of them assume
    that the underlying model only returns variable indices with nonnegative
    values. Use [`MOI.Bridges.LazyBridgeOptimizer`](@ref) instead.

## Example

```jldoctest var_singlebridgeoptimizer; setup=:(using MathOptInterface; const MOI = MathOptInterface)
julia> struct MyNewBridge{T} <: MOI.Bridges.Variable.AbstractBridge end

julia> bridge = MOI.Bridges.Variable.SingleBridgeOptimizer{MyNewBridge{Float64}}(
           MOI.Utilities.Model{Float64}(),
       )
MOIB.Variable.SingleBridgeOptimizer{MyNewBridge{Float64}, MOIU.Model{Float64}}
with 0 variable bridges
with inner model MOIU.Model{Float64}
```

## Implementation notes

All bridges should simplify the creation of `SingleBridgeOptimizer`s by defining
a constant that wraps the bridge in a `SingleBridgeOptimizer`.
```jldoctest var_singlebridgeoptimizer
julia> const MyNewBridgeModel{T,OT<:MOI.ModelLike} =
           MOI.Bridges.Variable.SingleBridgeOptimizer{MyNewBridge{T},OT};
```
This enables users to create bridged models as follows:
```jldoctest var_singlebridgeoptimizer
julia> MyNewBridgeModel{Float64}(MOI.Utilities.Model{Float64}())
MOIB.Variable.SingleBridgeOptimizer{MyNewBridge{Float64}, MOIU.Model{Float64}}
with 0 variable bridges
with inner model MOIU.Model{Float64}
```
"""
mutable struct SingleBridgeOptimizer{BT<:AbstractBridge,OT<:MOI.ModelLike} <:
               MOI.Bridges.AbstractBridgeOptimizer
    model::OT
    map::Map # index of bridged variable -> variable bridge
    var_to_name::Dict{MOI.VariableIndex,String}
    name_to_var::Union{Dict{String,MOI.VariableIndex},Nothing}
    con_to_name::Dict{MOI.ConstraintIndex,String}
    name_to_con::Union{Dict{String,MOI.ConstraintIndex},Nothing}
end

function SingleBridgeOptimizer{BT}(model::OT) where {BT,OT<:MOI.ModelLike}
    return SingleBridgeOptimizer{BT,OT}(
        model,
        Map(),
        Dict{MOI.VariableIndex,String}(),
        nothing,
        Dict{MOI.ConstraintIndex,String}(),
        nothing,
    )
end

bridges(::MOI.Bridges.AbstractBridgeOptimizer) = EmptyMap()

bridges(bridge::SingleBridgeOptimizer) = bridge.map

MOI.Bridges.supports_constraint_bridges(::SingleBridgeOptimizer) = false

function MOI.Bridges.supports_bridging_constrained_variable(
    ::SingleBridgeOptimizer{BT},
    S::Type{<:MOI.AbstractSet},
) where {BT}
    return supports_constrained_variable(BT, S)
end

function MOI.Bridges.is_variable_bridged(
    ::SingleBridgeOptimizer,
    ::Type{<:MOI.AbstractSet},
)
    return true
end

function MOI.Bridges.is_bridged(
    b::SingleBridgeOptimizer,
    S::Type{<:MOI.AbstractSet},
)
    return MOI.Bridges.supports_bridging_constrained_variable(b, S)
end

function MOI.Bridges.is_bridged(
    ::SingleBridgeOptimizer,
    ::Type{<:MOI.AbstractFunction},
    ::Type{<:MOI.AbstractSet},
)
    return false
end

function MOI.Bridges.is_bridged(
    ::SingleBridgeOptimizer,
    ::Type{<:MOI.AbstractScalarFunction},
)
    return false
end

function MOI.Bridges.bridge_type(
    ::SingleBridgeOptimizer{BT},
    ::Type{<:MOI.AbstractSet},
) where {BT}
    return BT
end

MOI.Bridges.bridging_cost(::SingleBridgeOptimizer, args...) = 1.0

MOI.Bridges.recursive_model(b::SingleBridgeOptimizer) = b.model