text:
import Base: ==, abs, dot, -, +, cross, first, last, isempty, *, size

#########################
# Drawing simple shapes #
#########################

"""Point in R^2"""
immutable Point
    x::Real
    y::Real
end
Point(p::Tuple{Real, Real}) = Point(p[1], p[2])
==(p::Point, q::Point) = (p.x == q.x) && (p.y == q.y)
-(p::Point, q::Point) = Point(p.x-q.x, p.y-q.y)
+(p::Point, q::Point) = Point(p.x+q.x, p.y+q.y)
*(p::Point, a::Number) = Point(a*p.x, a*p.y)
*(a::Number, p::Point) = Point(a*p.x, a*p.y)

dot(p::Point, q::Point) = p.x*q.x + p.y*q.y
abs(p::Point) = sqrt(dot(p, p))
cross(p::Point, q::Point) = p.x*q.y - p.y*q.x

first(p::Point) = p.x
last(p::Point) = p.y

# ----------------------------------------------------------------------------------------

# We are after making closed polygons. These are some primitives which constitute them
abstract Curve

immutable Line <: Curve
    p0::Point
    p1::Point

    function Line(p0, p1)
        @assert !(p0 == p1)
        new(p0, p1)
    end
end
+(l::Line, vec::Point) = Line(l.p0+vec, l.p1+vec)
first(l::Line) = l.p0
last(l::Line) = l.p1

immutable CircleArc <: Curve
    p0::Point
    center::Point
    p1::Point
    # NOTE: order of p0, p1 translates to orientation
    function CircleArc(p0, center, p1)
        @assert !(p0 == center) && !(p1 == center) && !(p0 == p1)
        new(p0, center, p1)
    end
end
+(c::CircleArc, vec::Point) = CircleArc(c.p0+vec, c.center+vec, c.p1+vec)
first(c::CircleArc) = c.p0
last(c::CircleArc) = c.p1

immutable EllipseArc <: Curve
    p0::Point
    center::Point
    p1::Point
    # NOTE: this is a defintion similar to GMSH with the condition that start is at major
    # point; so p0(start), center, p0, p1(end) is how you'd do it in GMSH. Beware that this
    # convenction enforces orientation
    function EllipseArc(p0, center, p1)
        @assert !(p0 == center) && !(p1 == center) && !(p0 == p1)
        @assert abs(p0-center) > abs(p1-center)
        new(p0, center, p1)
    end
end
+(c::EllipseArc, vec::Point) = EllipseArc(c.p0+vec, c.center+vec, c.p1+vec)
first(c::EllipseArc) = c.p0
last(c::EllipseArc) = c.p1

# ----------------------------------------------------------------------------------------

"""Shapes of cells"""
abstract Shape

immutable Circle <: Shape
    center::Point
    radius::Real

    function Circle(center, radius)
        @assert radius > 0
        new(center, radius)
    end
end
+(shape::Circle, vec::Point) = Circle(shape.center+vec, shape.radius)

immutable Ellipse <: Shape
    center::Point
    size_x::Real
    size_y::Real

    function Ellipse(center, x_axis, y_axis)
        @assert x_axis > 0 && y_axis > 0
        new(center, x_axis, y_axis)
    end
end
+(shape::Ellipse, vec::Point) = Ellipse(shape.center+vec, shape.size_x, shape.size_y)

immutable Rectangle <: Shape
    ll::Point
    size_x::Real
    size_y::Real

    function Rectangle(ll, size_x, size_y)
        @assert size_x > 0 && size_y > 0
        new(ll, size_x, size_y)
    end
end
+(shape::Rectangle, vec::Point) = Rectangle(shape.ll+vec, shape.size_x, shape.size_y)

immutable Square <: Shape
    ll::Point
    size::Real

    function Square(ll, size)
        @assert size > 0
        new(ll, size)
    end
end
+(shape::Square, vec::Point) = Square(shape.ll+vec, shape.size)

immutable ClosedPolygon <: Shape
    points::Vector{Point}

    function ClosedPolygon(points)
        @assert length(points) > 2
        @assert first(points) != last(points)
        new(points)
    end
end
+(shape::ClosedPolygon, vec::Point) = ClosedPolygon(map(p -> p+vec, shape.points))

Triangle(v0::Point, v1::Point, v2::Point) = ClosedPolygon([v0, v1, v2])

# This is a equilateral N-sided polygon with given volume
function NGon(N::Int, V::Real, x0::Real=0, y0::Real=0)
    @assert N > 2 && V > 0
    R = sqrt(2*V/N/sin(2*pi/N))
    angles = 2*pi/N*(1:N)
    ClosedPolygon(map(Point, zip(x0 + R*cos(angles), y0 + R*sin(angles))))
end

function NGonR(N::Int, R::Real, x0::Real=0, y0::Real=0)
    @assert N > 2 && R > 0
    angles = 2*pi/N*(1:N)
    ClosedPolygon(map(Point, zip(x0 + R*cos(angles), y0 + R*sin(angles))))
end

function make_path{T<:Curve}(curves::Vector{T})
    curve = first(curves)
    orientation = [1]
    path = Vector{Curve}([curve])

    for index in 2:length(curves)
        next = curves[index]
        linked = false
        # The first curve is orientated and the curves should follow each other
        # So either they agree on orientation
        link = (last(orientation) == 1) ? last(curve) : first(curve)

        if link == first(next)
            push!(orientation, 1)
            linked = true
        end
        # Or they don't
        if link == last(next)
            push!(orientation, -1)
            linked = true
        end
        # But they must be ordered
        @assert linked "Linking $(index): $(last(curve)), $(first(next)), $(last(next)) | $(next)"
        push!(path, next)
        curve = next
    end
    (path, orientation)
end

immutable Loop <: Shape
    curves::Vector{Curve}
    orientation::Vector{Int}

    function Loop(curves)
        @assert !isempty(curves)

        curves, orientation = make_path(curves)
        # Now see if it is closed
        last(orientation) == 1 && @assert first(first(curves)) == last(last(curves))
        last(orientation) == -1 && @assert first(first(curves)) == first(last(curves))

        new(curves, orientation)
    end
end
+(shape::Loop, vec::Point) = Loop(map(p -> p+vec, shape.curves))

# The points should be unique
function is_degenerate(poly::ClosedPolygon)
    points = poly.points
    n = length(poly.points)
    for i in 1:n
        pointi = points[i]
        for j in i+1:n-1
            pointj = points[j]
            pointi == pointj && return true
        end
    end
    false
end

# Checking convexity (not used anywhere, just for fun)
function is_convex(poly::ClosedPolygon)
    points = poly.points
    n = length(points)

    for k in 1:n
        prev = (k == 1) ? n : k-1
        next = (k == n) ? 1 : k+1

        outflow = points[next]-points[k]
        inflow = points[k]-points[prev]

        sign(cross(outflow, inflow)) > 0 && return false
    end
    true
end

# ----------------------------------------------------------------------------------------

"""Smallest rectangle that contains the shape"""
immutable BoundingBox
    ll::Point
    ur::Point

    function BoundingBox(ll, ur)
        @assert first(ll) <= first(ur)
        @assert last(ll) <= last(ur)
        new(ll, ur)
    end
end
size(bbox::BoundingBox, dim::Int) = (dim == 1) ? (bbox.ur.x-bbox.ll.x) : (bbox.ur.y-bbox.ll.y)

BoundingBox(l::Line) = BoundingBox([l.p0, l.p1])
BoundingBox(c::CircleArc) = BoundingBox([c.p0, c.p1])
BoundingBox(c::EllipseArc) = BoundingBox([c.p0, c.p1])

BoundingBox(shape::Circle) = BoundingBox(shape.center - Point(shape.radius, shape.radius),
                                         shape.center + Point(shape.radius, shape.radius))

BoundingBox(shape::Ellipse) = BoundingBox(shape.center - Point(shape.size_x, shape.size_y),
                                          shape.center + Point(shape.size_x, shape.size_y))

BoundingBox(shape::Rectangle) = BoundingBox(shape.ll,
                                            shape.ll + Point(shape.size_x, shape.size_y))

BoundingBox(shape::Square) = BoundingBox(shape.ll,
                                         shape.ll + Point(shape.size, shape.size))

BoundingBox(points::Vector{Point}) = BoundingBox(Point(minimum(map(first, points)),
                                                      minimum(map(last, points))),
                                                Point(maximum(map(first, points)),
                                                      maximum(map(last, points))))

BoundingBox(shape::ClosedPolygon) = BoundingBox(shape.points)

BoundingBox(loop::Loop) = BoundingBox(map(BoundingBox, loop.curves))

BoundingBox{T<:Shape}(shapes::Vector{T}) = BoundingBox(map(BoundingBox, shapes))

BoundingBox(b::BoundingBox) = b
BoundingBox(b::BoundingBox, B::BoundingBox) = BoundingBox([b.ll, b.ur, B.ll, B.ur])

function BoundingBox(boxes::Vector{BoundingBox})
    length(boxes) == 1 && return first(boxes)
    # Otherwise, don't want recursion
    BoundingBox([[b.ll for b in boxes]..., [b.ur for b in boxes]...])
end

# Collision between bounding boxes
function collides(b::BoundingBox, B::BoundingBox, tol=1E-13)
    const ll, ur, LL, UR = b.ll, b.ur, B.ll, B.ur
    !(((ur.x+tol < LL.x) || (UR.x + tol < ll.x)) || ((ur.y+tol < LL.y) || (UR.y + tol < ll.y)))
end

collides{T<:Shape}(b::BoundingBox, shape::T) = collides(b, BoundingBox(shape))
collides{T<:Shape}(shape::T, B::BoundingBox) = collides(BoundingBox(shape), B)
collides{T<:Shape, S<:Shape}(b::T, B::S) = collides(BoundingBox(b), BoundingBox(B))

# ----------------------------------------------------------------------------------------

"""Canvas is where we draw the shapes"""
type Canvas
    shapes::Vector{Shape}
    bbox::BoundingBox
end
Canvas{T<:Shape}(shapes::Vector{T})=Canvas(shapes, BoundingBox(shapes))
Canvas() = Canvas(Vector{Shape}(), BoundingBox(Point(-Inf, -Inf), Point(Inf, Inf)))

collides{T<:Shape}(canvas::Canvas, shape::T) = collides(canvas.bbox, BoundingBox(shape))
isempty(canvas::Canvas) = length(canvas.shapes) == 0

"""Add shape - modifies canvas, recomputed bbox. Only allowed for noncolliding"""
function +(canvas::Canvas, shape::Shape)
    if !isempty(canvas)
        @assert !collides(canvas, shape)
        canvas.bbox = BoundingBox(canvas.bbox, BoundingBox(shape))
    else
        canvas.bbox = BoundingBox(shape)
    end
    push!(canvas.shapes, shape)
    canvas
end

function +(canvas::Canvas, shape::ClosedPolygon)
    if !isempty(canvas)
        @assert !is_degenerate(shape)
        @assert !collides(canvas, shape)
        canvas.bbox = BoundingBox(canvas.bbox, BoundingBox(shape))
    else
        canvas.bbox = BoundingBox(shape)
    end
    push!(canvas.shapes, shape)
    canvas
end

BoundingBox(c::Canvas, C::Canvas) = BoundingBox(c.bbox, C.bbox)

function +(c::Canvas, C::Canvas)
    isempty(C) && return c
    # A copy
    if isempty(c)
        c.shapes = C.shapes
        c.bbox = C.bbox
        return c
    end

    c.bbox = BoundingBox(c, C)
    c.shapes = vcat(c.shapes, C.shapes)
    c
end

"""For later drawing increase the bounding rectangle"""
function set_bbox!(canvas::Canvas, ll::Point, ur::Point)
    @assert ll.x < canvas.bbox.ll.x && ll.y < canvas.bbox.ll.y
    @assert ur.x > canvas.bbox.ur.x && ur.y > canvas.bbox.ur.y
    canvas.bbox = BoundingBox(ll, ur)
end

function set_bbox!(canvas::Canvas, padx::Real, pady::Real)
    @assert padx > 0 && pady > 0
    ll = canvas.bbox.ll - Point(padx, pady)
    ur = canvas.bbox.ur + Point(padx, pady)
    canvas.bbox = BoundingBox(ll, ur)
end

#############################
# Drawing composites/tissue #
#############################

max_stars_repo_name: MiroK/emi-cylinders

max_issues_repo_name: MiroK/emi-cylinders

max_forks_repo_name: MiroK/emi-cylinders


text:
#= PROCESSES
Examples of Levy measures (functions ν_p, ν_n and σ_κ)
Tempered stable, truncated stable and Watanabe class
(compatible with the second versions of the functions above)
=#

using SpecialFunctions, GSL, Base.Threads, PoissonRandom

#= TemperedStable (tempered stable), which includes the models: CGMY and KoBoL

Levy measure:
ν(dx)/dx = c_±|x|^{-1-α_±} e^{-λ_±|x|}

Input:
1. b: drift
2-3. c_p, c_n: nonnegative factors such that c_p + c_n > 0
4-5. α_p, α_n: nonnegative factors in [0,2)
6-7. λ_p, λ_n: positive factors

Output:
Functions (b_ν,sν_p,sν_n,sum_sν_p,sum_sν_n,σ_ν), where
1. b_ν(x) = b-∫_{(-1,1)∖(-x,x)}y ν(dy),
2-5. sν_p, sν_n, sum_sν_p and sum_sν_n are as above and
6. σ_ν(x) = √(∫_{(-x,x)}y^2 ν(dy)).
=#
function TemperedStable_init(b::Number,c_p::Number,c_n::Number,α_p::Number,α_n::Number,λ_p::Number,λ_n::Number)
    b_ν(x) = b - ( c_p*λ_p^(α_p-1)*(sf_gamma_inc(1-α_p, x*λ_p)-sf_gamma_inc(1-α_p, λ_p)) -
        c_n*λ_n^(α_n-1)*(sf_gamma_inc(1-α_n, x*λ_n)-sf_gamma_inc(1-α_n, λ_n)) )
    σ_ν(x) = sqrt(c_p*λ_p^(α_p-2)*gamma(2-α_p) * gamma_inc(2-α_p, x*λ_p, 0)[1] +
        c_n*λ_n^(α_n-2)*gamma(2-α_n) * gamma_inc(2-α_n, x*λ_n, 0)[1])
    (sν_p,sum_sν_p) = tempered_power_sampler(α_p,λ_p,c_p)
    (sν_n,sum_sν_n) = tempered_power_sampler(α_n,λ_n,c_n)

    return (b_ν,sν_p,sν_n,sum_sν_p,sum_sν_n,σ_ν)
end

# (Rosinski's series representation)
# Samples the jumps of the measure y ↦ c y^{-α-1}e^{-λy}1_{y>x}dy
function tempered_power_sampler(α::Number,λ::Number,c::Number)
    function s(x::Number,t::Number)
        n = pois_rand(c*t*x^(-α)/α)
        X = min.( rand(n) .^ (-1/α) .* x , rand(Exponential(1/λ),n) .* rand(n) .^ (1/α))
        return X[X .> x]
    end
    function sum_s(x::Number,t::Number)
        n = pois_rand(c*t*x^(-α)/α)
        ap = 1/α
        an = -ap
        dist = Exponential(1/λ)

        if n > 64
            X = zeros(nthreads())
            Y = zeros(nthreads())
            @threads for i = 1:n
                Y[threadid()] = min(rand() ^ an * x , rand(dist) * rand() ^ ap)
                if Y[threadid()] > x
                    X[threadid()] += Y[threadid()]
                end
            end
            return sum(X)
        else
            X = min.(rand(n) .^ an .* x , rand(dist,n) .* rand(n) .^ ap)
            return sum(X[X .> x])
        end
    end
    function sum_s(x::Number,y::Number,t::Number)
        n = pois_rand(c*t*y^(-α)/α)
        ap = 1/α
        an = -ap
        dist = Exponential(1/λ)

        if n > 64
            X1 = zeros(nthreads())
            X2 = zeros(nthreads())
            Y = zeros(nthreads())
            @threads for i = 1:n
                Y[threadid()] = min(rand() ^ an * y , rand(dist) * rand() ^ ap)
                if Y[threadid()] > x
                    X1[threadid()] += Y[threadid()]
                    X2[threadid()] += Y[threadid()]
                elseif Y[threadid()] > y
                    X2[threadid()] += Y[threadid()]
                end
            end
            return (sum(X1),sum(X2))
        else
            X = min.(rand(n) .^ an .* y , rand(dist,n) .* rand(n) .^ ap)
            return (sum(X[X .> x]), sum(X[X .> y]))
        end
    end
    return (s,sum_s)
end

#= Stable process

Levy measure:
ν(dx)/dx = c_±|x|^{-1-α_±}

Input:
1. b: drift
2-3. c_p, c_n: nonnegative factors such that c_p + c_n > 0
4-5. α_p, α_n: nonnegative factors in [0,2)
6-7. λ_p, λ_n: positive factors

Output:
Functions (b_ν,sν_p,sν_n,sum_sν_p,sum_sν_n,σ_ν), where
1. b_ν(x) = b-∫_{(-1,1)∖(-x,x)}y ν(dy),
2-5. sν_p, sν_n, sum_sν_p and sum_sν_n are as above and
6. σ_ν(x) = √(∫_{(-x,x)}y^2 ν(dy)).
=#
function Stable_init(b::Number,c_p::Number,c_n::Number,α_p::Number,α_n::Number)
    b_ν(x) = b + c_p * (α_p != 1 ? (x^(1 - α_p) - 1) / (1-α_p) : log(x)) + c_p * (α_n != 1 ? (x^(1 - α_n) - 1) / (1-α_n) : log(x))
    σ_ν(x) = sqrt(c_p * x^(2 - α_p) / (2-α_p) + c_n * x^(2 - α_n) / (2-α_n))
    (sν_p,sum_sν_p) = power_sampler(α_p,c_p)
    (sν_n,sum_sν_n) = power_sampler(α_n,c_n)

    return (b_ν,sν_p,sν_n,sum_sν_p,sum_sν_n,σ_ν)
end

# Inversion of tail Levy measure
# Samples the jumps of the measure y↦cy^{-α-1}dy
function power_sampler(α::Number,c::Number)
    function s(x::Number,t::Number)
        return rand(pois_rand(c*t*x^(-α)/α)) .^ (-1/α) .* x
    end
    function sum_s(x::Number,t::Number)
        n = pois_rand(c*t*x^(-α)/α)
        a = -1/α
        if n > 32
            X = zeros(nthreads())
            @threads for i = 1:n
                X[threadid()] += rand() ^ a
            end
            return sum(X) * x
        else
            return sum(rand(n) .^ a) * x
        end
    end
    function sum_s(x::Number,y::Number,t::Number)
        n = pois_rand(c*t*y^(-α)/α)
        sc = x/y
        a = -1/α

        if n > 32
            X1 = zeros(nthreads())
            X2 = zeros(nthreads())
            Y = zeros(nthreads())
            @threads for i = 1:n
                Y[threadid()] = rand() ^ a
                X2[threadid()] += Y[threadid()]
                if Y[threadid()] > sc
                    X1[threadid()] += Y[threadid()]
                end
            end
            return (sum(X1) * x, sum(X2) * x)
        else
            X = rand(n) .^ a
            return (sum(X[X .> sc]) * x, sum(X) * x)
        end
    end
    return (s,sum_s)
end

#= (subclass of) Watanabe (continuous singular)

General Levy measure: for some m∈ℕ, m≥2,
ν = ∑_{n∈ℤ}c_{n,+}δ_{m^n}+c_{n,-}δ_{-m^n}
With ∑_{n∈ℕ}(c_{-n,+}+c_{-n,-})b^(-n)+∑_{n∈ℕ}c_{n,+}+c_{n,-}<∞
and max_{n∈ℤ}(c_+ + c_-)<∞.

Our class:
ν = ∑_{n∈ℕ}c_+ δ_{m^(-n)} + c_- δ_{-m^(-n)}.

Note:
Infinite activity iff min(r_+,r_-)≤1
Blumenthal-Getoor index = max(0,-log(max(r_+,r_-))/log(b)))

Input:
1. b: drift
2. m: integer ≥2
3-4. c_p, c_n: nonnegative numbers

Output:
Functions (b_ν,sν_p,sν_n,sum_sν_p,sum_sν_n,σ_ν), where
1. b_ν(x) = b-∫_{(-1,1)∖(-x,x)}y ν(dy),
2-5. sν_p, sν_n, sum_sν_p and sum_sν_n are as above and
6. σ_ν(x) = √(∫_{(-x,x)}y^2 ν(dy)).
=#
function Watanabe_init(b::Number,m::Integer,c_p::Number,c_n::Number)
    function b_ν(x::Number)
        return b - (c_p - c_n) * m^(-( x >= 1 ? 0. : floor(-log(x)/log(m))+1 )) / (1-1/m)
    end
    function σ_ν(x::Number)
        return sqrt((c_p + c_n) * m^(-2*( x >= 1 ? 0. : floor(-log(x)/log(m))+1 )) / (1-1/m^2))
    end
    (sν_p,sum_sν_p) = wat_sampler(c_p,m)
    (sν_n,sum_sν_n) = wat_sampler(c_n,m)

    return (b_ν,sν_p,sν_n,sum_sν_p,sum_sν_n,σ_ν)
end

# Computes the sum c * |{n∈ℕ:m^n≤1/x}|
function wat(c::Number,m::Integer)
    local f(x) = x >= 1 ? 0. :
        c * ceil(-log(x)/log(m))
    return f
end

# Samples from the density p(k) = 1 / |{n∈ℕ:m^n≤1/x}|
function wat_sampler(c::Number,m::Integer)
    function s(x::Number,t::Number)
        if x >= 1
            return  Float64[]
        else
            return [1/m^ceil(rand() * floor(-log(x)/log(m))) for i = 1:Int(pois_rand(c * t * ceil(-log(x)/log(m))))]
        end
    end
    function sum_s(x::Number,t::Number)
        return sum(s(x,t))
    end
    function sum_s(x::Number,y::Number,t::Number)
        aux = s(y,t)
        return (sum(aux .* (aux .> x)), sum(aux))
    end
    return (s,sum_s)
end

ext: jl


text:
using Images

function filter_gauss_localsigma(image, mask_order::Number=3, c::Number=0.05)

  # Checks if the mask is right
  if mod(mask_order,2)!=1
    error("The input value for the mask order must be an odd number >=3.")
  end

  # Stores the input image's type and convert it into a matrix
  tipo = typeof(image)
  try
    image = convert(Array{Float32}, image)
  catch er
    error(er," --- Problem converting the single band image into an array.")
  end

  # Takes the image's dimensions and creates 'a' (a variable to help in the convolution)
  d = size(image)
  a = floor(Int,mask_order/2)

  # Variables to receive the masks created and the filtered image
  mask = Matrix(mask_order,mask_order)
  img1 = Matrix(d[1]-mask_order+1,d[2]-mask_order+1)

  # the convolution process itself
  for i = (a+1):(d[1]-a)
    for j = (a+1):(d[2]-a)

      ai=i-a;
      ia=i+a;
      aj=j-a;
      ja=j+a;

      σ = (c/std(image[ai:ia,aj:ja], corrected=true))^2

      if σ==0.0 || (1/(2*pi*σ^2))==0.0
        mask[:,:] = 0.0
        mask[a+1,a+1] = 1.0
      else
        for x=-a:a
          for y=-a:a
            mask[x+a+1,y+a+1] = convert(Float32,exp(-(x^2 + y^2)/(2*(σ^2))))
          end
        end
        mask = mask/sum(mask)
      end

      img1[i-a,j-a] = convert(Float32,sum(mask.*image[ai:ia,aj:ja]))

      if img1[i-a,j-a] > 1.0
          img1[i-a,j-a] = 1.0
      end
    end
  end

  # Converts the filtered image to the same input image type
  img1 = convert(tipo, img1)

  # Returns the result image and end of the function
  return(img1)
end


text:
import HMCUtilities:
    VariableConstraint,
    constrain,
    free

struct TestAffineConstraint{N,TM,TMinv,Tb} <: HMCUtilities.VariableConstraint{N,N}
    M::TM
    Minv::TMinv
    b::Tb

    function TestAffineConstraint(M, b)
        n = length(b)
        @assert size(M) == (n, n)
        Minv = inv(M)
        return new{n,typeof(M),typeof(Minv),typeof(b)}(M, Minv, b)
    end
end

HMCUtilities.constrain(c::TestAffineConstraint, y) = c.Minv * (y - c.b)
HMCUtilities.free(c::TestAffineConstraint, x) = c.M * x + c.b

struct TestSphericalConstraint <: HMCUtilities.VariableConstraint{3,2} end

function HMCUtilities.constrain(c::TestSphericalConstraint, y)
    sinθ, cosθ = sincos(y[1])
    sinϕ, cosϕ = sincos(y[2])
    return [sinθ * cosϕ, sinθ * sinϕ, cosθ]
end

function HMCUtilities.free(c::TestSphericalConstraint, x)
    ϕ = atan(x[2], x[1])
    θ = acos(x[3])
    return [θ, ϕ]
end

struct TestSquareConstraint <: HMCUtilities.UnivariateConstraint end

HMCUtilities.constrain(c::TestSquareConstraint, y::Real) = y^2

HMCUtilities.free(c::TestSquareConstraint, x::Real) = sqrt(x)

@testset "defaults" begin
    @testset "one-to-one" begin
        M = [1.0 2.0;
             3.0 4.0]
        b = [5.0, 6.0]

        c = TestAffineConstraint(M, b)
        xs = [zeros(2), ones(2), [1.0, 2.0]]
        ys = [b, [8.0, 13.0], [10.0, 17.0]]
        logπxs = [-dot(x, x) / 2 for x in xs]
        ∇x_logπxs = [-x for x in xs]
        Minv = inv(M)
        logπys = logπxs .- log(2)
        ∇y_logπys = [Minv'x for x in ∇x_logπxs]

        for (x, y, logπx, ∇x_logπx, y_exp, logπy_exp, ∇y_logπy_exp) in zip(xs, ys, logπxs, ∇x_logπxs, ys, logπys, ∇y_logπys)
            test_constraint(
                c,
                x,
                y,
                logπx,
                ∇x_logπx,
                y_exp,
                logπy_exp,
                ∇y_logπy_exp,
                test_type_stability = false
            )
        end
    end

    @testset "non one-to-one" begin
        c = TestSphericalConstraint()
        xs = [[sqrt(3 / 2), sqrt(3 / 2), 1] ./ 2]
        ys = [π .* [1 / 3, 1 / 4]]
        κμ = 10.0 * [1.0, 0.0, 0.0]
        logπxs = [dot(x, κμ) for x in xs]  # vMF
        ∇x_logπxs = [κμ for x in xs]
        logπys = [(5 * sqrt(6) - log(4 / 3)) / 2]
        ∇y_logπys = [[5 / sqrt(2) + 1 / sqrt(3), -5 * sqrt(3 / 2)]]

        for (x, y, logπx, ∇x_logπx, y_exp, logπy_exp, ∇y_logπy_exp) in zip(xs, ys, logπxs, ∇x_logπxs, ys, logπys, ∇y_logπys)
            test_constraint(
                c,
                x,
                y,
                logπx,
                ∇x_logπx,
                y_exp,
                logπy_exp,
                ∇y_logπy_exp,
                test_type_stability = false
            )
        end
    end

    @testset "univariate" begin
        c = TestSquareConstraint()
        x = rand()
        y = -sqrt(x)
        y_exp = abs(y)
        logπx = x^2 / 2
        ∇x_logπx = x
        logπy_exp = logπx + log(2 * abs(y))
        ∇y_logπy_exp = 2 * y^3 + inv(y)

        @testset "scalar" begin
            vtypes = [Float64, Float32]
            test_constraint(
                c,
                x,
                y,
                logπx,
                ∇x_logπx,
                y_exp,
                logπy_exp,
                ∇y_logπy_exp;
                cvtypes=vtypes,
                fvtypes=vtypes,
                test_type_stability = false
            )
        end

        @testset "vector" begin
            vtypes = [Vector{Float64}, Vector{Float32}]
            test_constraint(
                c,
                [x],
                [y],
                logπx,
                [∇x_logπx],
                [y_exp],
                logπy_exp,
                [∇y_logπy_exp];
                cvtypes=vtypes,
                fvtypes=vtypes,
                test_type_stability = false
            )
        end
    end
end

@testset "IdentityConstraint" begin
    @testset "vector n = $n" for n in 1:3
        for m in 1:3
            c = HMCUtilities.IdentityConstraint(n)
            x = randn(n)
            y = copy(x)
            logπx = -dot(x, x) / 2
            ∇x_logπx = -x
            logπy_exp = copy(logπx)
            ∇y_logπy_exp = copy(∇x_logπx)
            test_constraint(c, x, y, logπx, ∇x_logπx, y, logπy_exp, ∇y_logπy_exp)
        end
    end

    @testset "scalar" begin
        c = HMCUtilities.IdentityConstraint(1)
        vtypes = [Float64, Float32]
        x = randn()
        y = copy(x)
        logπx = -dot(x, x) / 2
        ∇x_logπx = -x
        logπy_exp = copy(logπx)
        ∇y_logπy_exp = copy(∇x_logπx)
        test_constraint(
            c,
            x,
            y,
            logπx,
            ∇x_logπx,
            y,
            logπy_exp,
            ∇y_logπy_exp;
            cvtypes=vtypes,
            fvtypes=vtypes
        )
    end
end

@testset "LowerBoundedConstraint" begin
    vtypes = [Float64, Float32]

    @testset "lb=-5, y=0" begin
        c = HMCUtilities.LowerBoundedConstraint(-5.0)
        (y, x, logπx, ∇x_logπx, logπy_exp, ∇y_logπy_exp) = (0, -4, -8, 4, -8, 5)
        test_constraint(
            c,
            x,
            y,
            logπx,
            ∇x_logπx,
            y,
            logπy_exp,
            ∇y_logπy_exp;
            cvtypes=vtypes,
            fvtypes=vtypes
        )

        @testset "clamp" begin
            @test constrain(c, -Inf) ≈ -5 + eps()
        end
    end

    @testset "lb=5, y=-5" begin
        c = HMCUtilities.LowerBoundedConstraint(5.0)
        (y, x, logπx, ∇x_logπx, logπy_exp, ∇y_logπy_exp) = (
            -5, 5.00673794700, -12.5337124350, -5.00673794700, -17.5337124350, 0.966264865075
        )
        test_constraint(
            c,
            x,
            y,
            logπx,
            ∇x_logπx,
            y,
            logπy_exp,
            ∇y_logπy_exp;
            cvtypes=vtypes,
            fvtypes=vtypes
        )

        @testset "clamp" begin
            @test constrain(c, -Inf) ≈ 5 + eps()
        end
    end
end

@testset "UpperBoundedConstraint" begin
    vtypes = [Float64, Float32]

    @testset "ub=-5, y=0" begin
        c = HMCUtilities.UpperBoundedConstraint(-5.0)
        (y, x, logπx, ∇x_logπx, logπy_exp, ∇y_logπy_exp) = (0, -6, -18, 6, -18, -5)
        test_constraint(
            c,
            x,
            y,
            logπx,
            ∇x_logπx,
            y,
            logπy_exp,
            ∇y_logπy_exp;
            cvtypes=vtypes,
            fvtypes=vtypes
        )

        @testset "clamp" begin
            xbound = -5 - eps()
            @test constrain(c, -Inf) ≈ xbound
            @test constrain_with_pushlogpdf(c, -Inf)[1] ≈ xbound
            @test constrain_with_pushlogpdf_grad(c, -Inf)[1] ≈ xbound
        end
    end

    @testset "ub=5, y=-5" begin
        c = HMCUtilities.UpperBoundedConstraint(5.0)
        (y, x, logπx, ∇x_logπx, logπy_exp, ∇y_logπy_exp) = (
            -5, 4.99326205300, -12.4663329650, -4.99326205300, -17.4663329650, 1.03364433507
        )
        test_constraint(
            c,
            x,
            y,
            logπx,
            ∇x_logπx,
            y,
            logπy_exp,
            ∇y_logπy_exp;
            cvtypes=vtypes,
            fvtypes=vtypes
        )

        @testset "clamp" begin
            xbound = 5 - eps()
            @test constrain(c, -Inf) ≈ xbound
            @test constrain_with_pushlogpdf(c, -Inf)[1] ≈ xbound
            @test constrain_with_pushlogpdf_grad(c, -Inf)[1] ≈ xbound
        end
    end
end

@testset "BoundedConstraint" begin
    vtypes = [Float64]

    @testset "lb=0, ub=1" begin
        c = HMCUtilities.BoundedConstraint(0.0, 1.0)
        @testset "y=-20" begin
            (y, x, logπx, ∇x_logπx, logπy_exp, ∇y_logπy_exp) = (
                -20, 2.0611536181902035814e-9, 0.0, 0.0, -20, 1.0
            )
            test_constraint(
                c,
                x,
                y,
                logπx,
                ∇x_logπx,
                y,
                logπy_exp,
                ∇y_logπy_exp;
                cvtypes=vtypes,
                fvtypes=vtypes
            )
        end

        @testset "y=20" begin
            (y, x, logπx, ∇x_logπx, logπy_exp, ∇y_logπy_exp) = (
                20, 0.99999999793884638181, -0.5, -1, -20.5, -1
            )
            test_constraint(
                c,
                x,
                y,
                logπx,
                ∇x_logπx,
                y,
                logπy_exp,
                ∇y_logπy_exp;
                cvtypes=vtypes,
                fvtypes=vtypes
            )
        end

        @testset "clamp" begin
            xbound = 0 + eps()
            @test isapprox(constrain(c, -Inf), xbound, atol=1e-12)
            @test isapprox(constrain_with_pushlogpdf(c, -Inf)[1], xbound, atol=1e-12)
            @test isapprox(constrain_with_pushlogpdf_grad(c, -Inf)[1], xbound, atol=1e-12)

            xbound = 1 - eps()
            @test isapprox(constrain(c, Inf), xbound, atol=1e-12)
            @test isapprox(constrain_with_pushlogpdf(c, Inf)[1], xbound, atol=1e-12)
            @test isapprox(constrain_with_pushlogpdf_grad(c, Inf)[1], xbound, atol=1e-12)
        end
    end

    @testset "lb=-5, ub=5" begin
        c = HMCUtilities.BoundedConstraint(-5.0, 5.0)
        @testset "y=1" begin
            (y, x, logπx, ∇x_logπx, logπy_exp, ∇y_logπy_exp) = (
                1, 2.310585786, -2.669403338, -2.310585786, -1.993341620, -5.005004541
            )
            test_constraint(
                c,
                x,
                y,
                logπx,
                ∇x_logπx,
                y,
                logπy_exp,
                ∇y_logπy_exp;
                cvtypes=vtypes,
                fvtypes=vtypes
            )

        end

        @testset "clamp" begin
            xbound = -5 + eps()
            @test constrain(c, -Inf) ≈ xbound
            @test constrain_with_pushlogpdf(c, -Inf)[1] ≈ xbound
            @test constrain_with_pushlogpdf_grad(c, -Inf)[1] ≈ xbound

            xbound = 5 - eps()
            @test constrain(c, Inf) ≈ xbound
            @test constrain_with_pushlogpdf(c, Inf)[1] ≈ xbound
            @test constrain_with_pushlogpdf_grad(c, Inf)[1] ≈ xbound
        end
    end
end

@testset "TransformConstraint" begin
    (lb, ub) = (-rand(), rand())
    @test_throws AssertionError HMCUtilities.TransformConstraint(ub, lb)

    c = HMCUtilities.TransformConstraint(lb, Inf)
    @test c === HMCUtilities.LowerBoundedConstraint(lb)

    c = HMCUtilities.TransformConstraint(-Inf, ub)
    @test c === HMCUtilities.UpperBoundedConstraint(ub)

    c = HMCUtilities.TransformConstraint(lb, ub)
    @test c === HMCUtilities.BoundedConstraint(lb, ub)

    c = HMCUtilities.TransformConstraint()
    @test c === HMCUtilities.IdentityConstraint(1)
    c = HMCUtilities.TransformConstraint(-Inf, Inf)
    @test c === HMCUtilities.IdentityConstraint(1)
end

@testset "UnitVectorConstraint" begin
    vtypes = [Vector{Float64}, Vector{Float32}]

    @testset "n=$n" for n in 2:4
        c = HMCUtilities.UnitVectorConstraint(n)
        κμ = rand() * normalize(randn(n))

        y = randn(n)
        x = normalize(y)
        y_exp = copy(x)
        logπx = dot(x, κμ)  # vMF
        ∇x_logπx = κμ
        logπy_exp = logπx - dot(y, y) / 2
        ∇y_logπy_exp = (I - x * x') * ∇x_logπx / norm(y) - y

        test_constraint(
            c,
            x,
            y,
            logπx,
            ∇x_logπx,
            y_exp,
            logπy_exp,
            ∇y_logπy_exp;
            cvtypes=vtypes,
            fvtypes=vtypes,
            test_logdetJ_consistency = false
        )

        @testset "clamp" begin
            xbound = normalize(y)
            @test constrain(c, y) ≈ xbound
            @test constrain_with_pushlogpdf(c, y)[1] ≈ xbound
            @test constrain_with_pushlogpdf_grad(c, y)[1] ≈ xbound
        end
    end
end

@testset "UnitVectorScaledConstraint" begin
    vtypes = [Vector{Float64}, Vector{Float32}]

    @testset "n=$n, r=$r" for n in 2:4, r in (1.0, 3.0, 1.0)
        c = HMCUtilities.UnitVectorScaledConstraint(n, r)
        κμ = rand() * normalize(randn(n))

        y = randn(n) .* r
        x = normalize(y)
        y_exp = copy(x) .* r
        logπx = dot(x, κμ)  # vMF
        ∇x_logπx = κμ
        logπy_exp = logπx - dot(y, y) / r^2 / 2
        ∇y_logπy_exp = (I - x * x') * ∇x_logπx / norm(y) - y ./ r^2

        test_constraint(
            c,
            x,
            y,
            logπx,
            ∇x_logπx,
            y_exp,
            logπy_exp,
            ∇y_logπy_exp;
            cvtypes=vtypes,
            fvtypes=vtypes,
            test_logdetJ_consistency = false
        )

        @testset "clamp" begin
            xbound = normalize(y)
            @test constrain(c, y) ≈ xbound
            @test constrain_with_pushlogpdf(c, y)[1] ≈ xbound
            @test constrain_with_pushlogpdf_grad(c, y)[1] ≈ xbound
        end
    end
end

@testset "UnitSimplexConstraint" begin
    vtypes = [Vector{Float64}, Vector{Float32}]

    ys = (
        [-log(2)],
        [-log(2), log(2)],
        [-log(2), 0.0, log(2)]
    )

    xs = (
        [1/3, 2/3],
        [1/5, 8/15, 4/15],
        [1/7, 2/7, 8/21, 4/21]
    )

    Js = (
        [2/9, -2/9],
        [4/25 0.0; -8/75 8/45; -4/75 -8/45],
        [6/49 0.0 0.0;
        -2/49 4/21 0.0;
        -8/147 -8/63 8/63;
        -4/147 -4/63 -8/63]
    )

    ∇y_logdetJs = (
        [1/3],
        [2/5, -1/3],
        [3/7, 0.0, -1/3]
    )

    @testset "n=$n" for n in 2:4
        vtypes = [Vector{Float64}, Vector{Float32}]

        c = HMCUtilities.UnitSimplexConstraint(n)
        alpha = rand(n)

        x = xs[n - 1]
        y = ys[n - 1]
        J = Js[n - 1]
        ∇y_logdetJ = ∇y_logdetJs[n - 1]

        logπx, back = Zygote.pullback(x -> sum((alpha .- 1) .* log.(x)), x)
        ∇x_logπx = back(1.0)[1]

        logdetJ = _logabsdet(J'J)[1] / 2

        logπy_exp = logπx + logdetJ
        ∇y_logπy_exp = J' * ∇x_logπx .+ ∇y_logdetJ

        test_constraint(
            c,
            x,
            y,
            logπx,
            ∇x_logπx,
            y,
            logπy_exp,
            ∇y_logπy_exp;
            cvtypes=vtypes,
            fvtypes=vtypes
        )

        @testset "clamp" begin
            y = [-Inf, randn(n-2)...]
            xbound = constrain(c, y)
            @test isapprox(xbound[1], eps(); atol=1e-16)
            xbound = constrain_with_pushlogpdf(c, y)[1]
            @test isapprox(xbound[1], eps(); atol=1e-16)
            xbound = constrain_with_pushlogpdf_grad(c, y)[1]
            @test isapprox(xbound[1], eps(); atol=1e-16)

            y = [Inf, randn(n-2)...]
            xbound = constrain(c, y)
            @test xbound ≈ [1 - eps(), eps() * ones(n-1)...]
            @test xbound[1] < 1
            @test all(xbound[2:end] .> 0)
            xbound = constrain_with_pushlogpdf(c, y)[1]
            @test xbound ≈ [1 - eps(), eps() * ones(n-1)...]
            @test xbound[1] < 1
            @test all(xbound[2:end] .> 0)
            xbound = constrain_with_pushlogpdf_grad(c, y)[1]
            @test xbound ≈ [1 - eps(), eps() * ones(n-1)...]
            @test xbound[1] < 1
            @test all(xbound[2:end] .> 0)
        end
    end
end

@testset "JointConstraint" begin
    vtypes = [Vector{Float64}, Vector{Float32}]

    cs = [HMCUtilities.IdentityConstraint(3),
          HMCUtilities.LowerBoundedConstraint(0.0),
          HMCUtilities.UpperBoundedConstraint(1.0),
          HMCUtilities.BoundedConstraint(0.0, 1.0),
          HMCUtilities.UnitVectorConstraint(4),
          HMCUtilities.UnitSimplexConstraint(3)]

    y = Float64[]
    x = Float64[]
    y_exp = Float64[]
    logπx = randn()
    ∇x_logπx = Float64[]
    ∇y_logπy_exp = Float64[]
    logdetJ_exp = sum(cs) do (c)
        yi = randn(HMCUtilities.free_dimension(c))
        xi = constrain(c, yi)
        yi_exp = free(c, xi)
        ∇x_logπxi = randn(length(xi))
        _, pushlogpdf_grad = constrain_with_pushlogpdf_grad(c, yi)
        logdetJi, ∇y_logπyi = pushlogpdf_grad(0.0, ∇x_logπxi)
        push!(y, yi...)
        push!(x, xi...)
        push!(y_exp, yi_exp...)
        push!(∇x_logπx, ∇x_logπxi...)
        push!(∇y_logπy_exp, ∇y_logπyi...)
        return logdetJi
    end
    logπy_exp = logπx + logdetJ_exp

    jc = HMCUtilities.JointConstraint(cs...)

    test_constraint(
        jc,
        x,
        y,
        logπx,
        ∇x_logπx,
        y_exp,
        logπy_exp,
        ∇y_logπy_exp;
        cvtypes=vtypes,
        fvtypes=vtypes,
        test_logdetJ_consistency = false
    )

    @testset "IdentityConstraint merged" begin
        cs = [HMCUtilities.IdentityConstraint(3),
              HMCUtilities.IdentityConstraint(2),
              HMCUtilities.UpperBoundedConstraint(1.0),
              HMCUtilities.IdentityConstraint(2),
              HMCUtilities.IdentityConstraint(1)]
        jc = HMCUtilities.JointConstraint(cs...)

        cs2 = [HMCUtilities.IdentityConstraint(5),
               HMCUtilities.UpperBoundedConstraint(1.0),
               HMCUtilities.IdentityConstraint(3)]
        jc2 = HMCUtilities.JointConstraint(cs...)

        @test jc == jc2
    end
end


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export cosine_similarity, cosine_distance, angle_distance

using LinearAlgebra
import LinearAlgebra: normalize, normalize!

function normalize(X::AbstractVector{S}) where S<:AbstractVector
    normalize.(X)
end

function normalize!(X::AbstractVector{S}) where S<:AbstractVector
    for x in X
        normalize!(x)
    end

    X
end

function normalize!(X::AbstractMatrix)
    for x in eachcol(X)
        normalize!(x)
    end

    X
end

"""
    cosine_distance(a, b)::Float64

Computes the cosine distance between two vectors, it expects normalized vectors (see [normalize!](@ref) method).
Please use angle_distance if you are expecting a metric function (cosine_distance is a faster
alternative whenever the triangle inequality is not needed)
"""
function cosine_distance(a, b)::Float64
    1 - dot(a, b)
end

const π_2 = π / 2
"""
    angle_distance(a, b)::Float64

Computes the angle  between two SparseVector objects (sparse vectors).
It supposes that all vectors are normalized (see `normalize!` function)

"""
function angle_distance(a, b)::Float64
    d = dot(a, b)

    if d <= -1.0
        return π
    elseif d >= 1.0
        return 0.0
    elseif d == 0  # turn around for zero vectors, in particular for denominator=0
        return π_2
    else
        return acos(d)
    end
end


text:
struct NoGap <: NoDecision
    value::Float64
end

Base.show(io::IO, ng::NoGap) = print(io, """
    The lower and upper bounds for the root belief were both $(ng.value), so no tree was created.

    Use the default_action solver parameter to specify behavior for this case.
    """)


text:
"""
    AbstractMaximumWeightMaximalMatchingAlgorithm
Abstract type that allows users to pass in their preferred algorithm
"""
abstract type AbstractMaximumWeightMaximalMatchingAlgorithm end

"""
    LPAlgorithm <: AbstractMaximumWeightMaximalMatchingAlgorithm
Forces the maximum_weight_maximal_matching function to use a linear programming formulation.
"""
struct LPAlgorithm <: AbstractMaximumWeightMaximalMatchingAlgorithm end

function maximum_weight_maximal_matching(
    g::Graph,
    w::AbstractMatrix{T},
    algorithm::LPAlgorithm,
    solver = nothing
) where {T<:Real}
    if ! isa(solver, AbstractMathProgSolver)
        error("The keyword argument solver must be an AbstractMathProgSolver, as accepted by JuMP.")
    end

    return maximum_weight_maximal_matching_lp(g, solver, w)
end

"""
    HungarianAlgorithm <: AbstractMaximumWeightMaximalMatchingAlgorithm
Forces the maximum_weight_maximal_matching function to use the Hungarian algorithm.
"""
struct HungarianAlgorithm <: AbstractMaximumWeightMaximalMatchingAlgorithm end

function maximum_weight_maximal_matching(
    g::Graph,
    w::AbstractMatrix{T},
    algorithm::HungarianAlgorithm,
    solver = nothing
) where {T<:Real}
    return maximum_weight_maximal_matching_hungarian(g, w)
end

"""
    maximum_weight_maximal_matching{T<:Real}(g, w::Dict{Edge,T})

Given a bipartite graph `g` and an edge map `w` containing weights associated to edges,
returns a matching with the maximum total weight among the ones containing the
greatest number of edges.

Edges in `g` not present in `w` will not be considered for the matching.

A `cutoff` keyword argument can be given, to reduce computational times
excluding edges with weights lower than the cutoff.

Finally, a specific algorithm can be chosen (`algorithm` keyword argument);
each algorithm has specific dependencies. For instance:

- If `algorithm=HungarianAlgorithm()` (the default), the package Hungarian.jl is used.
  This algorithm is always polynomial in time, with complexity O(n³).
- If `algorithm=LPAlgorithm()`, the package JuMP.jl and one of its supported solvers is required.
  In this case, the algorithm relies on a linear relaxation on of the matching problem, which is
  guaranteed to have integer solution on bipartite graphs. A solver must be provided with
  the `solver` keyword parameter.

The returned object is of type `MatchingResult`.
"""
function maximum_weight_maximal_matching(
        g::Graph,
        w::AbstractMatrix{T};
        cutoff = nothing,
        algorithm::AbstractMaximumWeightMaximalMatchingAlgorithm = HungarianAlgorithm(),
        solver = nothing
    ) where {T<:Real}

    if cutoff != nothing && ! isa(cutoff, Real)
        error("The cutoff value must be of type Real or nothing.")
    end

    if cutoff != nothing
        return maximum_weight_maximal_matching(g, cutoff_weights(w, cutoff), algorithm, solver)
    else
        return maximum_weight_maximal_matching(g, w, algorithm, solver)
    end
end

"""
    cutoff_weights copies the weight matrix with all elements below cutoff set to 0
"""
function cutoff_weights(w::AbstractMatrix{T}, cutoff::R) where {T<:Real, R<:Real}
    wnew = copy(w)
    for j in 1:size(w,2)
        for i in 1:size(w,1)
            if wnew[i,j] < cutoff
                wnew[i,j] = zero(T)
            end
        end
    end
    wnew
end

@deprecate maximum_weight_maximal_matching(g::Graph, solver::AbstractMathProgSolver, w::AbstractMatrix{T}, cutoff::R) where {T<:Real, R<:Real} maximum_weight_maximal_matching(g, w, algorithm=LPAlgorithm(), cutoff=cutoff, solver=solver)


text:
using Random
using Flux
include("base.jl")

GENES = [
    "ConvBlock",
    "DSConvBlock",
    "MaxPool",
    "SkipConnection"
]

MUTATIONS = [
    "NM_InsertConvBlock",
    "NM_IncreaseConvFilters",
    "NM_InsertSkipConnection",
    "ANM_RemoveLayer",
    "ANM_PruneConvLayer",
    "ANM_Conv2DSConv"
]

LAYER_OUTPUT = Dense(256, 1, relu)

########################################
#            LEMONADE GENES            #
########################################

mutable struct MConvBlock <: MBlockGene
    gene_conv::MConv
    gene_bn::MBatchNorm
    gene_relu::MReLU
    in_dims::Tuple
    out_dims::Tuple
end

function MConvBlock(in_dims::Tuple, out_dims::Tuple)
    gene_conv = MConv(in_dims, out_dims)
    gene_bn = MBatchNorm(out_dims)
    gene_relu = MReLU(out_dims)
    MConvBlock(gene_conv, gene_bn, gene_relu, in_dims, out_dims)
end

function get_phenotype(prev_layer::MGene, cur_layer::MConvBlock)
    layers = Array{Any}(nothing, 3)
    layers[1] = get_phenotype(prev_layer, cur_layer.gene_conv)
    layers[2] = get_phenotype(cur_layer.gene_conv, cur_layer.gene_bn)
    layers[3] = get_phenotype(cur_layer.gene_bn, cur_layer.gene_relu)
    return layers
end

mutable struct MDepthwiseConvBlock <: MBlockGene
    gene_conv::MDepthwiseConv
    gene_bn::MBatchNorm
    gene_relu::MReLU
    in_dims::Tuple
    out_dims::Tuple
end

function MDepthwiseConvBlock(in_dims::Tuple, out_dims::Tuple)
    gene_conv = MDepthwiseConv(in_dims, out_dims)
    gene_bn = MBatchNorm(out_dims)
    gene_relu = MReLU(out_dims)
    MDepthwiseConvBlock(gene_conv, gene_bn, gene_relu, in_dims, out_dims)
end

function get_phenotype(prev_layer::MGene, cur_layer::MDepthwiseConvBlock)
    layers = Array{Any}(nothing, 3)
    ### TODO: placeholder?
    if cur_layer.gene_conv.out_dims[3] < prev_layer.out_dims[3] || cur_layer.gene_conv.out_dims[3] % prev_layer.out_dims[3] != 0
        println("changing channels")
        # cur_layer.gene_conv.out_dims[3] = prev_layer.out_dims[3] * 2
        cur_layer.gene_conv = MDepthwiseConv(cur_layer.gene_conv.in_dims, (prev_layer.out_dims[1], prev_layer.out_dims[2], prev_layer.out_dims[3] * 2))
        cur_layer.out_dims = cur_layer.gene_conv.out_dims
        cur_layer.gene_bn = MBatchNorm(cur_layer.gene_conv.out_dims)
        cur_layer.gene_relu = MReLU(cur_layer.gene_bn.out_dims)
    end
    ###
    layers[1] = get_phenotype(prev_layer, cur_layer.gene_conv)
    layers[2] = get_phenotype(cur_layer.gene_conv, cur_layer.gene_bn)
    layers[3] = get_phenotype(cur_layer.gene_bn, cur_layer.gene_relu)
    return layers
end

########################################
#          LEMONADE MUTATIONS          #
########################################

function mutate!(ind::Individual)
    # mutation = sample(MUTATIONS)
    mutation = "NM_InsertConvBlock"

    if mutation == "NM_InsertConvBlock"
        # genotype = ind.genotype

        insert!(ind.genotype, 2, MConv((40, 40, 16), (40, 40, 16)))
        # println(genotype)
        # insert!(ind.model, 2, Conv())
    end

end

max_stars_repo_name: mahmoud-al-najar/neuroevolution

max_issues_repo_name: mahmoud-al-najar/neuroevolution

max_forks_repo_name: mahmoud-al-najar/neuroevolution


text:
module MetaTools

include("parsed_function.jl")
include("emit_function.jl")
include("commutative.jl")

end # module


text:
using ExpressPathToRegex
using Test
using JSON

tests = JSON.parsefile(dirname(@__FILE__) * "/tests.json")

exec = (re::Regex, str) -> begin
  m = match(re, str)
  m==nothing ? nothing : [m.match; m.captures]
end

let test_path = "/user/:id",
    test_param = ExpressPathToRegex.Token(
      "id",
      "/",
      "/",
      false,
      false,
      "[^\\/]+?"
    )

    begin # arguments"
      begin # should accept an array of keys as the second argument"
        keys = []
        re, ret_keys = path_to_regex(test_path, keys, match_end=false)

        @test ret_keys == keys
        @test keys == Any[test_param]
        @test exec(re, "/user/123/show")==["/user/123", "123"]
      end

      begin # should work with keys as null"
        re, ret_keys = path_to_regex(test_path, match_end=false)

        @test ret_keys == [test_param]
        @test exec(re, "/user/123/show")==["/user/123", "123"]
      end
    end

    begin # tokens"
      tokens = ExpressPathToRegex.parse(test_path)

      begin # should expose method to compile tokens to regexp"
        re = ExpressPathToRegex.tokens_to_regex(tokens)

        @test exec(re, "/user/123")==["/user/123", "123"]
      end

      begin # should expose method to compile tokens to a path function"
        fn = ExpressPathToRegex.tokens_to_function(tokens)

        @test fn(Dict("id" => "123")) == "/user/123"
      end
    end
end

begin # rules"
  for test in tests
    path = test[1]
    opts = test[2]
    tokens = map(test[3]) do token
      if !(typeof(token) <: AbstractString)
        token = ExpressPathToRegex.Token(
          typeof(token["name"]) <: Integer ? string(token["name"]+1) : token["name"], # keys start at 1
          token["prefix"],
          token["delimiter"],
          token["optional"],
          token["repeat"],
          token["pattern"]
          )
      end
      token
    end
    matchCases = test[4]
    compileCases = length(test)>=5 ? test[5] : []
    opts_keyargs = []
    if opts != nothing
      for opt in opts
        push!(opts_keyargs, (Symbol(opt[1] == "end" ? "match_end" : opt[1]), opt[2]))
      end
    end

    keys = filter((token) -> !(typeof(token) <: AbstractString), tokens)

    # Parsing and compiling is only supported with string input.
    if typeof(path) <: AbstractString
      begin # should parse"
        @test ExpressPathToRegex.parse(path) == tokens
      end

      begin # compile"
        toPath = ExpressPathToRegex.compile(path)

        for io in compileCases
          input = io[1]
          output = io[2]

          # higher numeric token names aka keys by 1
          if input!=nothing
            parsed_input = Dict()
            for row in input
              try
                parsed_input[string(parse(Int, row[1])+1)]=row[2]
              catch
                parsed_input[row[1]] = row[2]
              end
            end
            input = parsed_input
          end

          if output != nothing
            begin # should compile using string(input)
              @test toPath(input) == output
            end
          else
            begin # should not compile using string(input)
              @test_throws ArgumentError toPath(input)
            end
          end
        end
      end

      begin # match (opts != nothing ? " using " * string(opts) : "")
        for io in matchCases
          input = io[1]
          output = io[2]
          begin "should" * (output != nothing ? " " : " not ") * "match " * string(input)
            re, ret_keys = path_to_regex(path; opts_keyargs...)

            @test ret_keys == keys
            @test exec(re, input) == output
          end
        end
      end
    end
  end
end

begin # compile errors"
  begin # should throw when a required param is undefined"
    toPath = ExpressPathToRegex.compile("/a/:b/c")

    # Expected "b" to be defined"
    @test_throws ArgumentError toPath(nothing)
  end

  begin # should throw when it does not match the pattern"
    toPath = ExpressPathToRegex.compile("/:foo(\\d+)")

    # Expected "foo" to match "\d+"
    @test_throws ArgumentError toPath(Dict("foo"=>"abc"))
  end

  begin # should throw when expecting a repeated value"
    toPath = ExpressPathToRegex.compile("/:foo+")

    # Expected "foo" to not be empty"
    @test_throws ArgumentError toPath(Dict("foo"=>[]))
  end

  begin # should throw when not expecting a repeated value"
    toPath = ExpressPathToRegex.compile("/:foo")

    # Expected "foo" to not repeat"
    @test_throws ArgumentError toPath(Dict("foo"=>[]))
  end

  begin # should throw when repeated value does not match"
    toPath = ExpressPathToRegex.compile("/:foo(\\d+)+")

    # Expected all "foo" to match "\d+"
    @test_throws ArgumentError toPath(Dict("foo"=>[1, 2, 3, 'a']))
  end
end


text:
#=
Copyright (c) 2015, Intel Corporation
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.
=#

module SeqTest
using ParallelAccelerator

#ParallelAccelerator.DomainIR.set_debug_level(4)
#ParallelAccelerator.ParallelIR.set_debug_level(4)
#ParallelAccelerator.CGen.set_debug_level(4)
#ParallelAccelerator.set_debug_level(4)

# testing sequential code taken from quant example
@acc function mmult(u::Array{Float64, 2}, v::Array{Float64, 2})
    m, l0 = size(u)
    l, n = size(v)
    # assert(l0 == l)
    r = Array(Float64, m, n)
    for i = 1:m
        for j = 1:n
            s = convert(Float64, 0)
            for k = 1:l
                s += u[i, k] * v[k, j]
            end
            r[i, j] = s
        end
    end
    return r
end

function test1()
    A = [2.0 1.0; 2.0 2.0]
    B = [4.0 2.0; 1.0 8.0]
    return mmult(A, B)
end

end

using Base.Test

println("Testing sequential code...")
@test SeqTest.test1() ≈ [9.0 12.0; 10.0 20.0]
println("Done testing sequential code.")


text:
function _class_group(c::Vector{BigInt})
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = Qx(c)
  K, a = number_field(f, cached = false)
  OK = lll(maximal_order(K))
  C, mC = class_group(OK)
  return BigInt(order(C))
end

function _class_group_batch(polys::Vector{fmpq_poly})
  res = Dict()
  for i in 1:length(polys)
    f = polys[i]
    c = BigInt[FlintZZ(coeff(f, j)) for j in 0:degree(f)]
    res[i] = @spawn _class_group(c)
  end

  cls = Dict()

  @sync for i in 1:length(polys)
    @async cls[i] = fetch(res[i])
  end
  return cls
end

#mutable struct NFDB
#  path::String
#  metadata::Dict
#  head::Int
#  length::Int
#  mmapped::Bool
#  mmap::Vector{UInt}
#  hash::UInt
#  field_names::Vector{Symbol}
#
#  function NFDB(path::String)
#    metadata = Dict()
#    f = open(path)
#    head = 0
#    while true
#      line = readline(f)
#      if line[1] == '#'
#        head += 1
#        if !(':' in line)
#          continue
#        end
#        i = 2
#        while line[i] != ':'
#          i += 1
#        end
#        key = strip(line[2:i-1])
#        val = strip(line[i+1:end])
#        metadata[key] = val
#      else
#        break
#      end
#    end
#    seekstart(f)
#    length = countlines(f)
#    close(f)
#    z =  new()
#    z.path = path
#    z.metadata = metadata
#    z.head = head
#    z.length = length - head
#    z.mmapped = false
#    return z
#  end
#end

mutable struct NFDBRecord{T}
  data::Dict{Symbol, Any}
  K::AnticNumberField

  function NFDBRecord{T}(data) where {T}
    z = new{T}()
    z.data = data
    return z
  end
end

mutable struct NFDB{T}
  meta::Dict{Symbol, Any}
  fields::Vector{NFDBRecord{T}}
  prop::Vector{Tuple{Symbol, Int}} # 0 not, 1 there, 2 partial

  function NFDB{T}() where {T}
    z = new{T}()
    z.meta = Dict{Symbol, Any}()
    z.fields = NFDBRecord{T}[]
    return z
  end
end

Base.length(D::NFDB) = length(D.fields)

function update_properties!(D::NFDB)
  prop = Tuple{Symbol, Int}[ (x[1], 0) for x in properties_comp ]

  if length(D.fields) == 0
    D.prop = prop
    return prop
  end

  for k in 1:length(prop)
    s = prop[k][1]
    j = 0
    for i in 1:length(D.fields)
      if haskey(D.fields[i], s)
        j += 1
      end
    end
    if j == length(D.fields)
      prop[k] = (s, 1)
    elseif j == 0
      prop[k] = (s, 0)
    else
      prop[k] = (s, 2)
    end
  end
  D.prop = prop

  return D
end

function NFDB(L::Vector{AnticNumberField}; compute = [])
  res = NFDB{1}()
  for K in L
    D = _create_record(K, compute = compute)
    push!(res.fields, D)
  end

  update_properties!(res)
  return res
end

function NFDB(L::Vector{NFDBRecord{1}})
  z = NFDB{1}()
  z.fields = L
  update_properties!(z)
  return z
end

function add_meta!(D::NFDB, p::Pair)
  D.meta[p[1]] = p[2]
  return D
end

function get_meta!(D::NFDB, p::Symbol)
  if haskey(D.meta, p)
    return D.meta[p]
  else
    error("Information $p not found in metadata")
  end
end

function Base.get(D::NFDB, i::Int)
  it = Iterators.drop(eachline(D.path), D.head + i - 1)
  line = IOBuffer(first(it))
end

function _get_NFDBRecord()
end

struct NFDBRecordInfo
  name_tuples
end

const properties_comp = Dict(:id => (Int, x -> UInt(hash(x))),
                             :deg => (Int, x -> degree(x)),
                             :poly => (fmpq_poly, x -> defining_polynomial(x)),
                             :discriminant => (fmpz, x -> discriminant(maximal_order(x))),
                             :signature => (Tuple{Int, Int}, x -> signature(x)),
                             :class_number => (fmpz, x -> order(class_group(maximal_order(x))[1])),
                             :class_group => (Vector{fmpz}, x -> elementary_divisors(class_group(maximal_order(x))[1])),
                             :is_cm => (Bool, x -> is_cm_field(x)[1]),
                             :relative_class_number => (fmpz, x -> begin
                                                          fl, tau = is_cm_field(x)
                                                          @assert fl
                                                          k, mk = fixed_field(x, tau)
                                                          hk = order(class_group(lll(maximal_order(k)))[1])
                                                          hK = order(class_group(maximal_order(x))[1])
                                                          @assert mod(hK, hk) == 0
                                                          return divexact(hK, hk)
                                                        end),
                              :is_normal => (Bool, x -> is_normal(x)),
                              :automorphism_group => (Tuple{Int, Int}, x -> find_small_group(automorphism_group(x)[1])[1]),
                              :regulator => (arb, x -> regulator(maximal_order(x))),
                              :lmfdb_label => (String, x -> ""),
                              :is_abelian => (Bool, x -> is_abelian(automorphism_group(x)[1])),
                              :non_simple => (Vector{fmpq_poly}, x -> non_simple_extension(x)),
                              :galois_group => (Tuple{Int, Int}, x -> error()))

for (k, v) in properties_comp
  @eval ($k)(D::NFDBRecord) = D[Symbol($k)]::($(v[1]))
end

Base.getindex(D::NFDBRecord, s) = getindex(D.data, s)

properties(D::NFDBRecord) = collect(keys(D.data))

function field(D::NFDBRecord; cached = false)
  if isdefined(D, :K)
    return D.K
  else
    f = D[:poly]
    K, a = NumberField(f, "a", cached = false)
    if cached
      D.K = K
    end
    return K
  end
end

function setindex!(D::NFDBRecord, s, k::Symbol)
  if !(k in record_info_v1.name_tuples)
    throw(error("asdsD"))
  end
  if haskey(D.data, k)
    throw(error("Data for $k already exists"))
  end

  if !(s isa properties_comp[k][1])
    throw(error("$s has the wrong type (expected $(properties_comp[k][1]))"))
  end

  D.data[k] = s
  return D
end

function compute!(D::NFDBRecord, s::Symbol)
  if haskey(D.data, s)
    return D.data[s]
  end
  K = field(D)
  d = _get(K, s)
  D.data[s] = d
  return d
end

function compute!(D::NFDB, S::Vector{Symbol})
  for s in S
    compute!(D, s)
  end
end

function compute!(D::NFDB, s::Symbol)
  PB = Pkg.GitTools.MiniProgressBar(header = "Computing $s")
  PB.max = length(D.fields)
  rate = 0.0
  length_eta = 0
  for i in 1:length(D.fields)
    compute!(D[i], s)
    PB.current = i
    Pkg.GitTools.showprogress(stdout, PB)#ETA)
    flush(stdout)
  end
end

# You may add stuff at the end
const record_info_v1 = NFDBRecordInfo([:id,
                                       :poly,
                                       :polredabs,
                                       :deg,
                                       :discriminant,
                                       :signature,
                                       :class_number,
                                       :class_group,
                                       :regulator,
                                       :is_cm,
                                       :relative_class_number,
                                       :is_tamely_ramified,
                                       :is_normal,
                                       :automorphism_group,
                                       :lmfdb_label,
                                       :is_abelian,
                                       :non_simple,
                                       :galois_group])

@assert length(record_info_v1.name_tuples) <= 56

function Base.write(io::IO, D::NFDBRecord{1})
  o = zero(UInt64)
  i = 1
  j = 0
  for d in record_info_v1.name_tuples
    o = _set_flag(o, i, haskey(D.data, d))
    if haskey(D.data, d)
      j += 1
    end
    i += 1
  end
  o = _set_version(o, UInt8(1))

  print(io, repr(o))
  print(io, ",")

  m = 0
  for d in record_info_v1.name_tuples
    if haskey(D.data, d)
      m += 1
      print(io, _stringify(D.data[d]))
      if m < j
        print(io, ",")
      end
    end
  end
  return nothing
end

Base.haskey(D::NFDBRecord, s::Symbol) = haskey(D.data, s)

function Base.get(D::NFDBRecord{1}, s::Symbol)
  if haskey(D, s)
    return D.data[s]
  else
    return missing
  end
end

function _get(K, s)
  if haskey(properties_comp, s)
    return (properties_comp[s][2])(K)
  else
    throw(error("Invalid property :$(s) of number field"))
  end
end

function _create_record(K::AnticNumberField; compute = [])
  f = defining_polynomial(K)
  data = Dict{Symbol, Any}()
  data[:poly] = f
  data[:signature] = signature(f)
  data[:deg] = degree(f)
  for p in compute
    c = _get(K, p)
    data[p] = c
  end
  D = NFDBRecord{1}(data)
  D.K = K
  return D
end

# 1 <= i <= 56
_get_flag(flags::UInt, i::Int; offset = 0) = Bool((flags >> (offset + i - 1)) & 1)

# 1 <= i <= 56
function _set_flag(flags::UInt, i::Int, b::Bool; offset = 0)
  if b
    return flags | 1 << (offset + i - 1)
  else
    return flags & ~(1 << (offset + i - 1))
  end
end

_get_version(flags::UInt) = (flags >> 56) % UInt8

_set_version(flags::UInt, v::UInt8) = flags |= (UInt(v) << 56)

_parse_as(x) = x

_parse_as(::Type{fmpq_poly}) = Vector{fmpq}

_parse_as(::Type{Vector{fmpq_poly}}) = Vector{Vector{fmpq}}

create(::Type{fmpq_poly}, v::Vector{fmpq}) = Hecke.Globals.Qx(v)

create(::Type{Vector{fmpq_poly}}, v::Vector{Vector{fmpq}}) = [ Hecke.Globals.Qx(w) for w in v]

create(T, v) = v

function _stringify(x::fmpq_poly)
  return _stringify([coeff(x, i) for i in 0:degree(x)])
end

function _stringify(x::Vector{T}) where {T}
  s = "["
  for i in 1:(length(x) - 1)
    s = s * _stringify(x[i]) * ","
  end
  s = s * _stringify(x[end]) * "]"
  return s
end

_stringify(x::Tuple{Int, Int}) = string(x)

_stringify(x) = string(x)

_stringify(x::arb) = _string(x)

function get_record(io::IO)
  data = Dict{Symbol, Any}()
  b, v  = _parse(UInt, io)
  ver = _get_version(v)
  @assert ver == 1
  b = Base.read(io, UInt8)
  for (i, s) in enumerate(record_info_v1.name_tuples)
    if !_get_flag(v, i)
      continue
    end
    T = properties_comp[s][1]
    S = _parse_as(T)
    z = _parse(S, io, b)
    b, d = z
    dd = create(T, d)
    data[s] = dd

    if eof(io)
      break
    end
    while !eof(io) && b == UInt8(' ')
      b = Base.read(io, UInt8)
    end
    @assert b == UInt8(',')
    b = Base.read(io, UInt8)
  end
  return NFDBRecord{Int(ver)}(data)
end

function _parse(::Type{Bool}, io, start = Base.read(io, UInt8))
  if start == UInt8('t')
    fl = true
  elseif start == UInt8('f')
    fl = false
  else
    throw(error("Not possible"))
  end
  if eof(io)
    b = UInt8(255)
  else
    b = Base.read(io, UInt8)
  end
  return b, fl
end

# I am too lazy for UInt
function _parse(::Type{UInt}, io, start = Base.read(io, UInt8))
  n = IOBuffer()
  b = start
  while !eof(io) && b == UInt8(' ')
    b = Base.read(io, UInt8)
  end

  @assert b == UInt('0')
  b = Base.read(io, UInt8)
  @assert b == UInt('x')
  b = Base.read(io, UInt8)

  lo = UInt8('0')
  up = UInt8('9')
  lo2 = UInt8('a')
  up2 = UInt8('f')

  while (lo <= b <= up) || (lo2 <= b <= up2)
    Base.write(n, b)
    if eof(io)
      b = UInt8(255)
      break
    end
    b = Base.read(io, UInt8)
  end

  nu = String(take!(n))

  num = parse(UInt, nu, base = 16)

  return b, num
end

function _parse(::Type{Int}, io, start = Base.read(io, UInt8))
  w = UInt8(' ')
  b = start
  while !eof(io) && b == w
    b = Base.read(io, UInt8)
  end
  exp = 10
  res = 0
  neg = false
  if b == UInt8('-')
    neg = true
    b = Base.read(io, UInt8)
  end
  lo = UInt8('0')
  up = UInt8('9')
  while !eof(io) && lo <= b <= up
    res = res * exp + (b - lo)
    b = Base.read(io, UInt8)
  end
  if neg
    res = -res
  end
  return b, res
end

function _parse(::Type{String}, io, start = Base.read(io, UInt8))
  b = start
  res = IOBuffer()
  @assert b == UInt8('"')
  b = Base.read(io, UInt8)
  while b != UInt8('"')
    Base.write(res, b)
    b = Base.read(io, UInt8)
  end

  if eof(io)
    b = UInt8(255)
  else
    b = Base.read(io, UInt8)
  end

  return b, String(take!(res))
end

function _parse(::Type{Vector{T}}, io, start = Base.read(io, UInt8)) where {T}
  res = T[]
  w = UInt8(' ')
  b = start
  while !eof(io) && b == w
    b = Base.read(io, UInt8)
  end
  @assert b == UInt8('[')
  b = Base.read(io, UInt8)
  while !eof(io) && b == w
    b = Base.read(io, UInt8)
  end
  while b != UInt8(']')
    b, q = _parse(T, io, b)
    push!(res, q)
    if b == UInt8(']')
      break
    end
    while !eof(io) && b == w
      b = Base.read(io, UInt8)
    end
    b = Base.read(io, UInt8)
    while !eof(io) && b == w
      b = Base.read(io, UInt8)
    end
  end
  if eof(io)
    b = UInt8(255)
  else
    b = Base.read(io, UInt8)
  end
  return b, res
end

function _parse(::Type{fmpz}, io, start = Base.read(io, UInt8))
  n = IOBuffer()
  b = start
  while !eof(io) && b == UInt8(' ')
    b = Base.read(io, UInt8)
  end
  lo = UInt8('0')
  up = UInt8('9')
  mi = UInt8('-')

  while (lo <= b <= up) || b == mi
    Base.write(n, b)
    if eof(io)
      b = UInt8(255)
      break
    end
    b = Base.read(io, UInt8)
  end

  nu = String(take!(n))

  if length(nu) < 18
    num = fmpz(parse(Int, nu))
  else
    num = parse(fmpz, nu)
  end

  return b, num
end

function _parse(::Type{fmpq}, io, start = Base.read(io, UInt8))
  n = IOBuffer()
  d = IOBuffer()
  read_den = false
  b = start
  lo = UInt8('0')
  up = UInt8('9')
  mi = UInt8('-')
  while ((lo <= b <= up) || b == mi || (!read_den && b == UInt('/')))
    if b == UInt8('/')
      read_den = true
      skip(io, 1)
      @goto read
    end

    if !read_den
      Base.write(n, b)
    else
      Base.write(d, b)
    end

    @label read
    if eof(io)
      b = UInt8(255)
      break
    end
    b = Base.read(io, UInt8)
  end

  nu = String(take!(n))

  if length(nu) < 18
    num = fmpz(parse(Int, nu))
  else
    num = parse(fmpz, nu)
  end

  if !read_den
    return b, fmpq(num)
  else
    de = String(take!(d))
    if length(de) < 18
      den = fmpz(parse(Int, de))
    else
      den = parse(fmpz, de)
    end
    return b, num//den
  end
end

function _parse(::Type{Tuple{Int, Int}}, io, start = Base.read(io, UInt8))
  w = UInt8(' ')
  b = start
  while !eof(io) && b == w
    b = Base.read(io, UInt8)
  end
  @assert b == UInt8('(')
  b, x1 = _parse(Int, io)
  while !eof(io) && b == UInt8(' ')
    b = Base.read(io, UInt8)
  end
  @assert b == UInt8(',')
  b, x2 = _parse(Int, io)
  @assert b == UInt8(')')
  if eof(io)
    return UInt8(255), (x1, x2)
  else
    b = Base.read(io, UInt8)
  end
  return b, (x1, x2)
end

function _parse(::Type{Tuple{S, T}}, io, start = Base.read(io, UInt8)) where {S, T}
  w = UInt8(' ')
  b = start
  while !eof(io) && b == w
    b = Base.read(io, UInt8)
  end
  @assert b == UInt8('(')
  b, x1 = _parse(S, io)
  while !eof(io) && b == UInt8(' ')
    b = Base.read(io, UInt8)
  end
  @assert b == UInt8(',')
  b, x2 = _parse(T, io)
  @assert b == UInt8(')')
  if eof(io)
    return UInt8(255), (x1, x2)
  else
    b = Base.read(io, UInt8)
  end
  return b, (x1, x2)
end

function perm_from_string(s)
  c, p = AbstractAlgebra.Generic.parse_cycles(s)
  if length(c) == 0
    n = 1
  else
    n = maximum(c)
  end
  cdec = AbstractAlgebra.Generic.cycledec(c, p, n)
  return SymmetricGroup(cdec.cptrs[end]-1)(cdec)
end

function _parse(::Type{Perm{Int}}, io, start = Base.read(io, UInt8))
  b = start
  @assert b == UInt8('p')
  b = Base.read(io, UInt8)
  @assert b == UInt8('e')
  b = Base.read(io, UInt8)
  @assert b == UInt8('r')
  b = Base.read(io, UInt8)
  @assert b == UInt8('m')
  b = Base.read(io, UInt8)
  @assert b == UInt8('"')
  b = Base.read(io, UInt8)
  res = IOBuffer()
  while b != UInt8('"')
    Base.write(res, b)
    b = Base.read(io, UInt8)
  end

  if eof(io)
    b = UInt8(255)
  else
    b = Base.read(io, UInt8)
  end

  return b, perm_from_string(String(take!(res)))
end

function _parse(::Type{arb}, io, start = Base.read(io, UInt8))
  n = IOBuffer()
  b = start
  while !eof(io) && b == UInt8(' ')
    b = Base.read(io, UInt8)
  end

  while b != UInt8(']')
    Base.write(n, b)
    if eof(io)
      b = UInt8(255)
      break
    end
    b = Base.read(io, UInt8)
  end
  nu = String(take!(n))

  RR = ArbField(64)

  if eof(io)
    b = UInt8(255)
  else
    b = Base.read(io, UInt8)
  end

  return b, RR(nu)
end

function _string(a::arb)
  s = string(a)
  if s[1] != '['
    s = "[" * s * "]"
  end
  return s
end

function _stringify(a::Bool)
  if a
    return "t"
  else
    return "f"
  end
end

function Base.show(io::IO, D::NFDB)
  k = length(D.meta)
  if k == 0
    println(io, "Number field table no metadata")
  else
    println(io, "Number field table with metadata:")
  end

  for (i, (p, e)) in enumerate(D.meta)
    print(io, p, ": ", e)
    print(io, "\n")
  end

  print(io, "Data present: ")
  for (s, i) in D.prop
    if i == 1
      print(io, s, ", ")
    elseif i == 1
      print(io, "(", s, "), ")
    end
  end
  println(io)

  print(io, "No. fields: ", length(D.fields))
end

function Base.show(io::IO, D::NFDBRecord)
  println(io, "Number field record with data:")
  for (p, e) in D.data
    print(io, p, ": ", e)
    print(io, "\n")
  end
end

function Base.write(f::String, D::NFDB)
  if isfile(f)
    throw(error("File $f already exists. Please move first"))
  end
  open(f, "w") do io
    Base.write(io, D)
  end
end

function Base.write(io::IO, D::NFDB)
  for (p, e) in D.meta
    print(io, "# ", p, ":", e, "\n")
  end
  for i in 1:length(D.fields)
    Base.write(io, D.fields[i])
    print(io, "\n")
  end
end

Base.getindex(D::NFDB, i::Int) = D.fields[i]

Base.eltype(::Type{NFDB{T}}) where {T} = NFDBRecord{T}

function Base.read(file::String, ::Type{NFDB})
  metadata = Dict()
  f = open(file)
  head = 0
  local line
  while !eof(f)
    line = readline(f)
    if line[1] == '#'
      head += 1
      if !(':' in line)
        continue
      end
      i = 2
      while line[i] != ':'
        i += 1
      end
      key = Symbol(strip(line[2:i-1]))
      val = strip(line[i+1:end])
      metadata[key] = val
    else
      break
    end
  end
  DB = NFDB{1}()
  if line[1] != '#'
    D = get_record(IOBuffer(line))
    push!(DB.fields, D)
    while !eof(f)
      line = readline(f)
      D = get_record(IOBuffer(line))
      push!(DB.fields, D)
    end
  end
  close(f)
  DB.meta = metadata

  update_properties!(DB)

  return DB
end

function _latexify(f::fmpq_poly, dollar = true)
  if dollar
    return "\$" * AbstractAlgebra.obj_to_latex_string(f) * "\$"
  else
    return AbstractAlgebra.obj_to_latex_string(f)
  end
end

function _latexify(fs::Vector{fmpq_poly})
  s = IOBuffer()
  print(s, "\$")
  for i in 1:length(fs)
    print(s, _latexify(fs[i], false))
    if i < length(fs)
      print(s, ", ")
    end
  end
  print(s, "\$")
  return String(take!(s))
end

_latexify(a::fmpz) = AbstractAlgebra.obj_to_latex_string(a)

_latexify(a::Int) = AbstractAlgebra.obj_to_latex_string(a)

_latexify(a::Tuple{Int, Int}) = string(a)

function latex_table(d::NFDBRecord; entries = [:deg, :automorphism_group, :poly, :non_simple, :class_number, :discriminant], mpage_for_poly = true)
  io = IOBuffer()
  for i in 1:length(entries)
    if entries[i] == :discriminant
      print(io, _latexify(factor(d[entries[i]]), withnumber = false, approx = false))
    elseif entries[i] == :poly
      if mpage_for_poly
        print(io, "\\begin{minipage}{4cm}\\begin{center}\n\$\n\\begin{aligned}\n")
        print(io, _latexify(d[entries[i]], false))
        print(io, "\n\\end{aligned}\n\$\n\\end{center}\\end{minipage}\n")
      else
        print(io, _latexify(d[entries[i]]))
      end
    else
      print(io, _latexify(d[entries[i]]))
    end
    if i < length(entries)
      print(io, " & ")
    end
  end
  print(io, "\\\\")
  return String(take!(io))
end

function latex_table(D::NFDB; entries = [:deg, :automorphism_group, :poly, :non_simple, :class_number, :discriminant])
  io = IOBuffer()
  for d in D
    print(io, latex_table(d, entries = entries))
    print(io, "\n")
  end
  return String(take!(io))
end

function _latexify(f::Fac{fmpz}; withnumber = true, approx = true, cdot = false, sign = true)
  l = Tuple{fmpz, Int}[(p, e) for (p, e) in f]
  sort!(l, by = x -> x[1])
  if length(l) == 0
    return "\$1\$"
  end
  m = prod([p^e for (p, e) in f]) * unit(f)
  if withnumber
    s = "\$ $m = "
  else
    s = "\$"
  end
  if sign && m < 0
    s = s * "-"
  end
  for i in 1:(length(l) - 1)
    if l[i][2] == 1
      s = s * "$(l[i][1])\\cdot "
    else
      if cdot
        s = s * "$(l[i][1])^{$(l[i][2])} \\cdot "
      else
        s = s * "$(l[i][1])^{$(l[i][2])} "
      end
    end
  end
  if l[end][2] == 1
    s = s * "$(l[end][1])"
  else
    s = s * "$(l[end][1])^{$(l[end][2])}"
  end
  if approx
    s = s * "\\approx " * (m < 0 ? "-" : "") * "10^{$(Int(ceil(log(10, abs(m) * 1.0))))}"
  end
  s = s * "\$"
  return s
end

################################################################################
#
#    Merge
#
################################################################################

function Base.merge!(R::NFDB{1}, D1::NFDB{1})
  sizehint!(R.fields, length(R) + length(D1))
  append!(R.fields, D1.fields)
  update_properties!(R)
  return R
end

function Base.merge(D1::NFDB{1}, D2::NFDB{1})
  D3 = NFDB{1}()
  merge!(D3, D1)
  merge!(D3, D2)
  return D3
end

function Base.merge(D::Vector{NFDB{1}})
  if length(D) == 1
    return D
  end

  R = NFDB{1}()
  for i in 1:length(D)
    merge!(R, D[i])
  end

  return R
end

names32 = [ "C32", "C2.C4^2", "C4*C8", "C2.OD16", "C2^2:C8", "C2^2.D4",
           "C2^3.C4", "OD16.C2", "D4:C4", "Q8:C4", "C4wrC2", "C4:C8", "C8:C4",
           "C2.D8", "C8.C4", "C2*C16", "OD32", "D16", "SD32", "Q32", "C2*C4^2",
           "C2*C2^2:C4", "C2*C4:C4", "(C2*C4):C4", "C4*D4", "C4*Q8",
           "C2^2wrC2", "C2^2:D4", "C2^2:Q8", "(C2^2*C4):C2", "C4.D4", "C4.Q8",
           "C4^2:C2", "C4:D4", "C4:Q8", "C2^2*C8", "C2*OD16", "OD16:C2",
           "C2*D8", "C2*SD16", "C2*Q16", "D8:C2", "C8:C2^2", "SD16:C2",
           "C2^3*C4", "C2^2*D4", "C2^2*Q8", "C2*D4:C2", "Q8:C2^2", "C4.C2^3",
           "C2^5" ]

names64 = [ "C64", "C8^2", "C8:C8", "C2^3:C8", "(C2*C4):C8", "D4:C8", "Q8:C8",
           "C2^2.SD16", "C2.(C2^2.D4)", "C2.(C2^3.C4)", "C2.C4wrC2", "C4.D8",
           "C4.Q16", "C2.Q8:C4", "C2.C8:C4", "C2.(C2.D8)", "C2.(C4*C8)",
           "(C2*OD16).C2", "C2^3.Q8", "(C2*C4).D4", "C2.D4:C4", "OD16.C4",
           "C2^2:C4:C4", "OD16:C4", "(C2*C8):C4", "C4*C16", "C2.OD32",
           "OD32.C2", "C2^2:C16", "C2^3.C8", "D4.C8", "C2wrC4", "(C2^2*C4):C4",
           "C4^2:C4", "C2^2.D4.C2", "(OD16.C2):C2", "OD16.C2.C2", "C2^2.D8",
           "C2.Q32", "D8.C4", "OD32:C2", "(C2*D8).C2", "Q16.C4", "C4:C16",
           "C8.C8", "C16:C4", "C2.D16", "C2.SD32", "C16.C4", "C2*C32", "OD64",
           "D32", "SD64", "Q64", "C4^3", "C2*C2.C4^2", "C4.C4^2", "C4*C2^2:C4",
           "C4*C4:C4", "C2^5.C2", "(C2^3*C4).C2", "C2^2.(C2*D4)", "C4.C4:C4",
           "C2.C4^2.C2", "C2.C4:Q8", "C2^2:(C4:C4)", "C2^3.D4", "C4:C4:C4",
           "C2^4.C2^2", "C4:(C4:C4)", "C4:(C2^2:C4)", "C2.(C4*Q8)", "C2.C4:D4",
           "C2.(C4.D4)", "(C2^2*D4).C2", "(C2^2*Q8).C2", "(C2.C4^2):C2",
           "C2.C4^2:C2", "(C2*C4).Q8", "C2^2.(C2*Q8)", "C2.C2^2:Q8",
           "C2^2.D4:C2", "C2*C4*C8", "C2*C2.OD16", "C4*OD16", "C2.(C2*C4^2)",
           "C2*C2^2:C8", "C2^2:OD16", "(C2^2*C8):C2", "C2*C2^2.D4",
           "(C2*D4):C4", "C2*C2^3.C4", "C2*OD16.C2", "(C2^3.C4):C2",
           "C2*D4:C4", "C2*Q8:C4", "D4:C2:C4", "Q8:(C2*C4)", "(C2*C8):C2^2",
           "(C2*Q8):C4", "C2*C4wrC2", "OD16:C2^2", "C2*C4:C8", "C4:OD16",
           "(C2^2*C8).C2", "C2*C8:C4", "C2*C2.D8", "C8.(C2*C4)", "C8:(C2*C4)",
           "C2*C8.C4", "(C8.C4):C2", "C4.OD16", "C4^2.C4", "C2^3.(C2*C4)",
           "C8*D4", "C2.(C4*D4)", "(C2*D4).C4", "C4*D8", "C4*SD16", "C4*Q16",
           "SD16:C4", "Q16:C4", "D8:C4", "(C4*C8):C2", "OD16:C2:C2", "C8*Q8",
           "(C2*Q8).C4", "(C2*D8):C2", "C2.(C2*SD16)", "C4:C4.C2^2",
           "(C2*Q8):C2^2", "C2^2:Q16", "(C2*Q16):C2", "D4:D4", "C4wrC2:C2",
           "D4.D4", "C2.C2^2wrC2", "C2wrC2^2", "(C2^2.D4):C2", "C4:D8",
           "C4:D4.C2", "(C2*SD16).C2", "C4:Q16", "(C2*SD16):C2", "C4:C8:C2",
           "C2^2:SD16", "C2^2:D8", "(C2.D8):C2", "C4.(C2*D4)", "C2.C8:C2^2",
           "C2^2:Q8.C2", "SD16:C2:C2", "C8:C2^2:C2", "Q8.D4", "D4:Q8", "Q8:Q8",
           "D4:C4.C2", "C2.(C2*Q16)", "D4.Q8", "Q8.Q8", "C2^2:C8:C2",
           "C2^2:D4.C2", "C2.D8:C2", "C2^2:C8.C2", "C2^2.Q16", "C8:C4:C2",
           "C4.SD16", "(C4*C8).C2", "Q8:C4:C2", "C4:Q8:C2", "D4:C4:C2",
           "(C2*C8).C2^2", "C4:SD16", "(C2*D4).C2^2", "(C2*Q16).C2", "C8.D4",
           "C8:D4", "(C2.OD16):C2", "C8:Q8", "C8.Q8", "C2.(C2*D8)", "C8:C4.C2",
           "C2^2*C16", "C2*OD32", "C16.C2^2", "C2*D16", "C2*SD32", "C2*Q32",
           "D16:C2", "C16:C2^2", "Q32:C2", "C2^2*C4^2", "C2^2*C2^2:C4",
           "C2^2*C4:C4", "C2*(C2*C4):C4", "C2*C4*D4", "C2*C4*Q8", "C4*D4:C2",
           "C2^3:(C2*C4)", "C2.(C2^3*C4)", "D4:(C2*C4)", "C2*C2^2wrC2",
           "C2*C2^2:D4", "C2*C2^2:Q8", "C2*(C2^2*C4):C2", "C4^2:C2^2",
           "C2*C4.D4", "C2*C4.Q8", "C2*C4^2:C2", "(C2*C4^2):C2", "C2*C4:D4",
           "C2*C4:Q8", "(C2*C4):D4", "(C2*C4^2).C2", "C2^4:C2^2",
           "(C2^2*D4):C2", "(C2^2*Q8):C2", "C2.(C2^2*D4)", "C2.Q8:C2^2",
           "(C4*D4):C2", "C4:D4:C2", "(C4*Q8):C2", "C4^2:C2:C2", "C2^3:Q8",
           "(C2*C4):Q8", "D4^2", "(C2*D4):C2^2", "C4:(D4:C2)", "C2^2:(D4:C2)",
           "D4*Q8", "Q8:D4", "C2^2wrC2.C2", "(C4.Q8):C2", "C4^2.C2^2",
           "(C4*D4).C2", "(C4.D4):C2", "(C2*Q8).C2^2", "(C4*Q8).C2", "Q8^2",
           "C2^3.C2^3", "(C2*C4).C2^3", "C2^2wrC2:C2", "C2^2:D4:C2",
           "C2^2:Q8:C2", "C2^2.C2^4", "C2^3*C8", "C2^2*OD16", "C2*OD16:C2",
           "(C2*OD16):C2", "C2^2*D8", "C2^2*SD16", "C2^2*Q16", "C2*D8:C2",
           "C2*C8:C2^2", "C2*SD16:C2", "SD16:C2^2", "C8.C2^3", "C4.C2^4",
           "D10.C2^2", "C2^4*C4", "C2^3*D4", "C2^3*Q8", "C2^2*D4:C2",
           "C2*Q8:C2^2", "C2*C4.C2^3", "D4.C2^3", "C2^6" ]

function has_obviously_relative_class_number_not_one(K::AnticNumberField, is_normal::Bool = true, maxdeg::Int = degree(K))
  if is_normal
    subs = subfields_normal(K)
  else
    subs = subfields(K)
  end

  sort!(subs, by = x -> degree(x[1]))

  for (L, mL) in subs
    if degree(L) > min(degree(K) - 1, maxdeg)
      continue
    end
    fl, tau = is_cm_field(L)
    if !fl
      continue
    end

    l = fixed_field(L, tau)[1]

    h = order(class_group(lll(maximal_order(L)))[1])
    hp = order(class_group(lll(maximal_order(l)))[1])
    @assert mod(h, hp) == 0
    hm = divexact(h, hp)
    if hm == 3 || hm > 4
      return true, L
    end
  end
  return false, K
end

################################################################################
#
#  Iterator interface
#
################################################################################

function Base.iterate(D::NFDB, i = 1)
  if i > length(D)
    return nothing
  end
  return D[i], i + 1
end


text:
using DiscreteFiltering
using LinearAlgebra
using Plots
using ProgressLogging
using OrdinaryDiffEq: RK4

## Parameters

# Domain
a = 0.0
b = 2π
domain = PeriodicIntervalDomain(a, b)

# Exact solutions (advection equation)
u(x, t) = sin(x - t) + 3 / 5 * cos(5(x - t)) + 7 / 25 * sin(20(x - 1 - t))
u_int(x, t) = -cos(x - t) + 3 / 25 * sin(5(x - t)) - 7 / 25 / 20 * cos(20(x - 1 - t))

## Time
T = 1.0

# Ode solver tolerances
# tols = (;)
# tols = (; abstol = 1e-6, reltol = 1e-4)
# tols = (; abstol = 1e-7, reltol = 1e-5)
tols = (; abstol = 1e-10, reltol = 1e-8)
degmax = 1000
subspacedim = 1000
λ = 1e-8

# Amplitudes
c = [
    [1],
    [1, 3],
    [1, 3 / 5, 1 / 25],
    rand(5),
]

# Frequencies
ω = [
    [0],
    [2, 3],
    [1, 5, 20],
    (rand(1:10, 5)),
]

# Phase-shifts
ϕ = [
    [π / 2],
    [0, 0],
    [0, 0, 20],
    2π * rand(5),
]

# Number of time steps for fitting C̄
nₜ = 500
t = LinRange(0, T, nₜ)
UU = sum_of_sines.([domain], c, ω, ϕ)
u₀_list = [UU[1] for UU ∈ UU]
U₀_list = [UU[2] for UU ∈ UU]

# Number of mesh points
nrefine = 10
MM = floor.(Int, 10 .^ LinRange(1, 4, nrefine))
NN = MM
# NN = floor.(Int, 5 // 4 .* MM)
# λλ = [0, 0, 0, 0, 1e-1]

# Errors
err = zeros(nrefine)
err_bar = zeros(nrefine)
err_allbar = zeros(nrefine)

## Solve
@time for i = 1:nrefine
    # λ = λλ[i]
    λ = 1e-8
    λ_ridge = 1e-8
    # M = 100
    M = MM[i]
    # N = NN[i]
    N = 1000 # NN[i]
    # M = N

    println("Solving for M = $M and N = $N")
    # Discretization
    x = discretize(domain, M)
    ξ = discretize(domain, N)
    Δx = (b - a) / M

    # Filter
    h₀ = 1.0Δx
    # h₀ = (b - a) / 200
    # h₀ = Δx / 2
    h(x) = h₀ * (1 - 1 / 2 * cos(x))
    # h(x) = 0.49999Δx
    filter = TopHatFilter(h)
    # filter = ConvolutionalFilter(x -> 2h(x), x -> abs(x) ≤ h(x) ? 1/ 2h(x) : zero(x))
    # h(x) = 2h₀
    # σ = h₀/2
    # filter = GaussianFilter(h, σ)

    C̄ = fit_Cbar(domain, filter, u₀_list, U₀_list, M, t; λ, method = :ridge)

    # Equations
    equation = AdvectionEquation(domain, IdentityFilter())
    equation_filtered = AdvectionEquation(domain, filter)

    # Solve discretized problem
    sol = solve(
        equation,
        ξ -> u(ξ, 0.0),
        (0.0, T),
        M,
        M;
        method = "discretizefirst",
        subspacedim,
        tols...,
        degmax,
        λ,
    )

    # Solve filtered-then-discretized problem
    # sol_bar = solve(
    #     equation_filtered,
    #     ξ -> u(ξ, 0.0),
    #     (0.0, T),
    #     M,
    #     M;
    #     method = "filterfirst",
    #     subspacedim,
    #     tols...,
    #     degmax,
    #     λ,
    # )

    # Solve discretized-then-filtered problem
    sol_allbar = solve(
        equation_filtered,
        ξ -> u(ξ, 0.0),
        (0.0, T),
        M,
        N;
        # method = "discretizefirst",
        # method = "discretizefirst-without-R",
        method = "discretizefirst-fit-Cbar",
        # solver = RK4(),
        subspacedim,
        tols...,
        degmax,
        λ,
        λ_ridge,
        C̄,
    )

    # Exact filtered solution
    ū = apply_filter_int(ξ -> u_int(ξ, T), filter, domain)
    # ū = apply_filter(ξ -> u(ξ, T), filter, domain)

    # Relative error
    # u_exact = u.(ξ, T)
    u_exact = u.(x, T)
    ū_exact = ū.(x)
    # err[i] = norm(sol(T) - u_exact) / norm(u_exact)
    err[i] = norm(sol(T) - ū_exact) / norm(ū_exact)
    # err_bar[i] = norm(sol_bar(T) - ū.(ξ)) / norm(ū.(ξ))
    # err_bar[i] = norm(sol_bar(T) - ū.(x)) / norm(ū.(x))
    err_allbar[i] = norm(sol_allbar(T) - ū_exact) / norm(ū_exact)
    # println("λ = $λ, err = $(err_allbar[i])")
    println("λ = $λ, err_bar = $(err_bar[i]), err_allbar = $(err_allbar[i])")
end

## Set PGFPlotsX plotting backend to export to tikz
pgfplotsx()

## In terminal plots
unicodeplots()

## Default backend
gr()

##
p = plot(
    xaxis = :log10,
    yaxis = :log10,
    minorgrid = true,
    # size = (400, 300),
    legend = :bottomleft,
    # xlims = (10, 10^4),
    ylims = (1e-6, 1e-0),
    xticks = 10 .^ (1:4),
    xlabel = "M",
    # title = "N = 1000"
);
plot!(p, MM, err_allbar; label = "Data driven Cbar, HF", marker = :d, color = 2);
# plot!(p, MM, err_allbar_lf; label = "Data driven Cbar, LF", marker = :d, color = 1);
plot!(p, MM, err; label = "Without filter, HF", marker = :c, linestyle = :dash, color = 2);
# plot!(p, MM, err_lf; label = "Without filter, LF", marker = :c, linestyle = :dash, color = 1);
# display(p)
# output = "output/advection/Cbar.tikz"
output = "output/advection/Cbar.pdf"
savefig(p, output)


text:
function ELBO(model::GP, pr_mean, kernel)
    setpr_mean!(model.f, pr_mean)
    setkernel!(model.f, kernel)
    compute_kernel_matrices!(model, true)
    return log_py(model)
end

@traitfn function ELBO(model::TGP, pr_means, kernels) where {TGP <: AbstractGP; IsFull{TGP}}
    setpr_means!(model, pr_means)
    setkernels!(model, kernels)
    compute_kernel_matrices!(model, true)
    return ELBO(model)
end

@traitfn function ELBO(
    model::TGP, pr_means, kernels, Zs
) where {TGP <: AbstractGP; !IsFull{TGP}}
    setpr_means!(model, pr_means)
    setkernels!(model, kernels)
    setZs!(model, Zs)
    compute_kernel_matrices!(model, true)
    return ELBO(model)
end


text:
@testset "legend_merging" begin
    mwe_data = (; t = 1:20, y1 = ones(20), y2 = rand(20))
    mwe_names = ["f_n", "f_d"]
    plt = data(mwe_data) *
        mapping("t", ["y1", "y2"] .=> "y"; color = dims(1) => (i -> mwe_names[i]), marker=dims(1)) *
        (visual(Lines) + visual(Scatter))
    @test_throws ArgumentError draw(plt)
end


text:
using Documenter, LLVM

makedocs(
    modules = [LLVM],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "LLVM.jl",
    pages = [
        "Home"    => "index.md",
        "Manual"  => [
            "man/usage.md",
            "man/troubleshooting.md"
        ],
        "Library" => [
            "lib/api.md",
            "lib/interop.md"
        ]
    ],
    doctest = true
)

deploydocs(
    repo = "github.com/maleadt/LLVM.jl.git"
)


text:
using OctreeBH
using StaticArrays
using PyPlot
#plt.style.use("dark_background")
#plt.style.use("default")

function plot_tree_grid(node::Node{N,T}, ix, iy, cc, ax) where {N,T}
	@assert N>=2
	if isLeaf(node)
    	xmin = node.center[ix] - 0.5*node.length[ix]
    	xmax = node.center[ix] + 0.5*node.length[ix]
    	ymin = node.center[iy] - 0.5*node.length[iy]
    	ymax = node.center[iy] + 0.5*node.length[iy]
    	color=cc
    	ax.plot([xmin,xmin],[ymin,ymax], c=color)
    	ax.plot([xmin,xmax],[ymin,ymin], c=color)
    	ax.plot([xmax,xmax],[ymin,ymax], c=color)
    	ax.plot([xmin,xmax],[ymax,ymax], c=color)
	else
		@inbounds for i in 1:2^N
	        #println("open this node")
			plot_tree_grid(node.child[i], ix, iy, cc, ax)
	    end
	end
end

const BOXSIZE = 1.0

const T = Float64
const N = 2 #spatial dimension
const Npart = 300 #number of particles
const hsml0 = 0.15 #searching radius

using Random
Random.seed!(1114)

#randomly distributing particles
X = [@SVector rand(N) for _ in 1:Npart]

topnode_length = @SVector(ones(N)) * BOXSIZE  #actual length of tree
center = topnode_length .* 0.5

hsml = ones(Npart) .* hsml0
mass = ones(Npart) #particle mass
part = [Data{N,T}(SVector(X[i]), i, hsml[i], mass[i]) for i in eachindex(X)]

#build the tree
@time tree = buildtree(part, center, topnode_length);

@show max_depth = get_max_tree_depth(tree)

#set_max_depth_AMR!(tree, max_depth)

boxsizes = @SVector(ones(N)) * BOXSIZE  #for periodic B.C.

ms = 6
ix, iy = 1, 2

clf()
fig, ax = subplots(1, 1, figsize=(8, 8))
ax.axis([0,boxsizes[ix],0,boxsizes[iy]])
ax.set_aspect(1)
#ax.plot(getindex.(X,ix), getindex.(X,iy), ".", c="tab:blue", ms=ms)
@time plot_tree_grid(tree, ix, iy, "grey", ax)
savefig("balance_none.png")

#max_depth = 9
@time for i in max_depth:-1:1
	println("balancing the tree at level ", i)
	OctreeBH.balance!(tree, i, tree.length[1])

	clf()
	fig, ax = subplots(1, 1, figsize=(8, 8))
	ax.axis([0,boxsizes[ix],0,boxsizes[iy]])
	ax.set_aspect(1)
	@time plot_tree_grid(tree, ix, iy, "grey", ax)
	savefig("balance_" * string(i) * ".png")
end


text:
using JuMP; using NCDatasets; using AxisArrays; using CalliopeJuMP.Util;

function load_energy_balance_constraints(model_dict)
    sets = model_dict["sets"]
    constraint_dict = Dict()

    model_dict["expressions"]["system_balance"] = @expression(
        model_dict["backend_model"],
        [loc_carrier=sets["loc_carriers"],
                       timestep=sets["timesteps"]], 0
    )

    constraint_dict["system_balance_constraint"] = (
        build_constraint(
            ["loc_carriers_system_balance_constraint", "timesteps"],
            "system_balance", model_dict
        )
    )

    if haskey(sets, "loc_techs_finite_resource_supply_plus")
        constraint_dict["resource_availability_supply_plus"] = (
            build_constraint(["loc_techs_finite_resource_supply_plus", "timesteps"],
                "resource_availability_supply_plus", model_dict)
        )
    end

    if haskey(sets, "loc_techs_finite_resource_supply")
        constraint_dict["balance_supply_constraint"] = (
            build_constraint(["loc_techs_finite_resource_supply", "timesteps"],
                "balance_supply", model_dict)
        )
    end

    if haskey(sets, "loc_techs_balance_demand_constraint")
        constraint_dict["balance_demand_constraint"] = (
            build_constraint(["loc_techs_balance_demand_constraint", "timesteps"],
                "balance_demand", model_dict)
        )
    end

    if haskey(sets, "loc_techs_balance_transmission_constraint")
        constraint_dict["balance_transmission_constraint"] = (
            build_constraint(["loc_techs_balance_transmission_constraint", "timesteps"],
                "balance_transmission", model_dict)
        )
    end

    if haskey(sets, "loc_techs_balance_supply_plus_constraint")
        constraint_dict["balance_supply_plus_constraint"] = (
            build_constraint(["loc_techs_balance_supply_plus_constraint", "timesteps"],
                "balance_supply_plus", model_dict)
        )
    end

    if haskey(sets, "loc_techs_balance_storage_constraint")
        constraint_dict["balance_storage_constraint"] = (
            build_constraint(["loc_techs_balance_storage_constraint", "timesteps"],
                "balance_storage", model_dict)
        )
    end

    return constraint_dict
end

function system_balance_constraint_rule(backend_model, set_indices, model_dict)
    """
    System balance ensures that, within each location, the production and
    consumption of each carrier is balanced.
    """
    loc_carrier, timestep = set_indices
    _prod, _con, exporting = get_loc_tech_carriers(model_dict, loc_carrier)
    variables = model_dict["variables"]
    expressions = model_dict["expressions"]
    expressions["system_balance"][loc_carrier, timestep] = (
        sum(variables["carrier_prod"][loc_tech_carrier, timestep]
            for loc_tech_carrier in _prod) +
        sum(variables["carrier_con"][loc_tech_carrier, timestep]
            for loc_tech_carrier in _con)
    )

    return @constraint(backend_model,
        expressions["system_balance"][loc_carrier, timestep] == 0
    )
end

# Supply Balance
function balance_supply_constraint_rule(backend_model, set_indices, model_dict)
    loc_tech, timestep = set_indices
    variables = model_dict["variables"]
    sets = model_dict["sets"]
    loc_tech_carrier = model_dict["parameters"]["lookup_loc_techs"][loc_tech]

    energy_eff = get_param(model_dict, "energy_eff", [loc_tech], timestep)

    if energy_eff == 0
        energy_eff = 0
    else
        energy_eff = 1 / get_param(model_dict, "energy_eff", [loc_tech], timestep)
    end

    @expression(backend_model, prod,
        energy_eff * variables["carrier_prod"][loc_tech_carrier, timestep]
    )

    @expression(backend_model, available_resource,
        get_param(model_dict, "resource", [loc_tech], timestep) *
        get_param(model_dict, "resource_scale", [loc_tech], timestep)
    )

    if loc_tech_is_in(sets, loc_tech, "loc_techs_area")
        available_resource *= variables["resource_area"][loc_tech]
    end

    # creating the constraint just requires putting together the affine expressions
    if get_param(model_dict, "force_resource", [loc_tech], timestep) == true
        return @constraint(backend_model, prod == available_resource)
    else
        return @constraint(backend_model, prod <= available_resource)
    end
end

# Demand Balance
function balance_demand_constraint_rule(backend_model, set_indices, model_dict)
    """
    Limit consumption from demand techs to their required resource
    """
    loc_tech, timestep = set_indices
    variables = model_dict["variables"]
    sets = model_dict["sets"]
    loc_tech_carrier = model_dict["parameters"]["lookup_loc_techs"][loc_tech]

    @expression(backend_model, con,
        variables["carrier_con"][loc_tech_carrier, timestep]
        * get_param(model_dict, "energy_eff", [loc_tech], timestep)
    )

    @expression(backend_model, available_resource,
        get_param(model_dict, "resource", [loc_tech], timestep) *
        get_param(model_dict, "resource_scale", [loc_tech], timestep)
    )

    if loc_tech_is_in(sets, loc_tech, "loc_techs_area")
        available_resource *= variables["resource_area"][loc_tech]
    end
    # creating the constraint just requires putting together the affine expressions
    if get_param(model_dict, "force_resource", [loc_tech], timestep) == true
        return @constraint(backend_model, con == available_resource)
    else
        return @constraint(backend_model, con >= available_resource)
    end
end

# Resource availablity
function resource_availability_supply_plus_constraint_rule(backend_model, set_indices, model_dict)
    """
    Limit production from supply_plus techs to their available resource
    """
    loc_tech, timestep = set_indices
    variables = model_dict["variables"]
    sets = model_dict["sets"]

    @expression(backend_model, available_resource,
        get_param(model_dict, "resource", [loc_tech], timestep) *
        get_param(model_dict, "resource_scale", [loc_tech], timestep) *
        get_param(model_dict, "resource_eff", [loc_tech], timestep)
    )

    if loc_tech_is_in(sets, loc_tech, "loc_techs_area")
        available_resource *= variables["resource_area"][loc_tech]
    end

    if get_param(model_dict, "force_resource", [loc_tech], timestep) == true
        return @constraint(backend_model, variables["resource_con"][loc_tech, timestep] == available_resource)
    else
        return @constraint(backend_model, variables["resource_con"][loc_tech, timestep] <= available_resource)
    end
end

# Transmission Balance
function balance_transmission_constraint_rule(backend_model, set_indices, model_dict)
    """
    Balance carrier production and consumption of transmission technologies
    """
    parameters = model_dict["parameters"]
    variables = model_dict["variables"]
    loc_tech, timestep = set_indices

    loc_tech_carrier = parameters["lookup_loc_techs"][loc_tech]
    remote_loc_tech = parameters["lookup_remotes"][loc_tech]
    remote_loc_tech_carrier = parameters["lookup_loc_techs"][remote_loc_tech]

    return @constraint(backend_model, variables["carrier_prod"][loc_tech_carrier, timestep]
        == -1 * variables["carrier_con"][remote_loc_tech_carrier, timestep]
        * get_param(model_dict, "energy_eff", [loc_tech], timestep))
end

# Supply_plus Balance
function balance_supply_plus_constraint_rule(backend_model, set_indices, model_dict)
    """
    Balance carrier production and resource consumption of supply_plus technologies
    alongside any use of resource storage
    """
    loc_tech, timestep = set_indices
    variables = model_dict["variables"]
    sets = model_dict["sets"]
    parameters = model_dict["parameters"]
    loc_tech_carrier = parameters["lookup_loc_techs"][loc_tech]
    total_eff = (get_param(model_dict, "energy_eff", [loc_tech], timestep) *
        get_param(model_dict, "parasitic_eff", [loc_tech], timestep))

    if total_eff == 0
        _prod = 0
    else
        _prod = variables["carrier_prod"][loc_tech_carrier, timestep] / total_eff
    end

    if !loc_tech_is_in(sets, loc_tech, "loc_techs_store")
        return @constraint(backend_model, variables["resource_con"][loc_tech, timestep] == _prod)
    else
        timestep_index = findin(sets["timesteps"], [timestep])[1]
        if timestep_index == 1
            storage_previous_step = get_param(model_dict, "storage_initial", [loc_tech])
        else
            storage_loss = get_param(model_dict, "storage_loss", [loc_tech])
            previous_timestep = model_dict["sets"]["timesteps"][timestep_index - 1]
            time_resolution = model_dict["parameters"]["timestep_resolution"][previous_timestep]
            storage_previous_step = (
                ((1 - storage_loss) ^ time_resolution) *
                variables["storage"][loc_tech, previous_timestep]
            )
        end
        return @constraint(backend_model,
            variables["storage"][loc_tech, timestep] ==
            storage_previous_step + variables["resource_con"][loc_tech, timestep] - _prod
        )
    end
end

# Storage Balance
function balance_storage_constraint_rule(backend_model, set_indices, model_dict)
    loc_tech, timestep = set_indices
    variables = model_dict["variables"]
    sets = model_dict["sets"]
    parameters = model_dict["parameters"]
    loc_tech_carrier = parameters["lookup_loc_techs"][loc_tech]

    energy_eff = get_param(model_dict, "energy_eff", [loc_tech], timestep)

    if energy_eff == 0
        _prod = 0
    else
        _prod = variables["carrier_prod"][loc_tech_carrier, timestep] / energy_eff
    end

    _con = variables["carrier_con"][loc_tech_carrier, timestep] * energy_eff

    timestep_index = findin(sets["timesteps"], [timestep])[1]
    if timestep_index == 1
        storage_previous_step = get_param(model_dict, "storage_initial", [loc_tech])
    else
        storage_loss = get_param(model_dict, "storage_loss", [loc_tech])
        previous_timestep = model_dict["sets"]["timesteps"][timestep_index - 1]
        time_resolution = model_dict["parameters"]["timestep_resolution"][previous_timestep]
        storage_previous_step = (
            ((1 - storage_loss) ^ time_resolution) *
            variables["storage"][loc_tech, previous_timestep]
        )
    end

    return @constraint(backend_model,
        variables["storage"][loc_tech, timestep] == storage_previous_step - _prod - _con
    )

end

max_stars_repo_name: calliope-project/CalliopeJuMP

max_issues_repo_name: calliope-project/CalliopeJuMP

max_forks_repo_name: calliope-project/CalliopeJuMP