text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "TMatrix"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/JuliaRemoteSensing/TMatrix-fortran.git", "bf0d1bf15ff382b562427492895f8f3d0fb6b177")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/TMatrix-fortran/
mkdir build && cd build
meson --cross-file=${MESON_TARGET_TOOLCHAIN}
meson install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms(;experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libtmatrixro", :tmatrix_random_orientation),
    LibraryProduct("libtmatrixfo", :tmatrix_fixed_orientation)
]

# Dependencies that must be installed before this package can be built
dependencies = [Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
##### NEW KEYS #####

struct ShapeKey <: Gadfly.GuideElement
    title::AbstractString
    labels::Vector{String}
    pos::Vector
end
ShapeKey(;title="Shape", labels=[""], pos=Float64[]) = ShapeKey(title, labels, pos)

"""
    Guide.shapekey[(; title="Shape", labels=[""], pos=Float64[])]
    Guide.shapekey(title, labels, pos)

Enable control of the auto-generated shapekey.  Set the key `title` and the item `labels`.
`pos` overrides [Theme(key_position=)](@ref Gadfly) and can be in either
relative (e.g. [0.7w, 0.2h] is the lower right quadrant), absolute (e.g. [0mm,
0mm]), or plot scale (e.g. [0,0]) coordinates.
"""
const shapekey = ShapeKey

function Guide.render(guide::Guide.ShapeKey, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)

    (theme.key_position == :none) && return Gadfly.Guide.PositionedGuide[]
    gpos = guide.pos
    (theme.key_position == :inside) && (gpos == Float64[]) &&  (gpos = [0.7w, 0.25h])

    # Aesthetics for keys: shape_key_title, shape_label (Function), shape_key_shapes (AbstractDict)
    nshapes = length(unique(aes.shape))
    guide_title = (guide.title!="Shape" || aes.shape_key_title==nothing) ? guide.title : aes.shape_key_title
    shape_key_labels = !(guide.labels==[""]) ? guide.labels : aes.shape_label(1:nshapes)

    colors = [nothing]
    if (aes.shape_key_title !=nothing)  && (aes.color_key_title==aes.shape_key_title)
        colors = collect(keys(aes.color_key_colors))
    end

    title_context, title_width = render_key_title2(guide_title, theme)
    ctxs = render_discrete_key(shape_key_labels, title_context, title_width, theme, shapes=1:nshapes, colors=colors)

    position = right_guide_position
    if gpos != Float64[]
        position = over_guide_position
        ctxs = [compose(context(), (context(gpos[1],gpos[2]), ctxs[1]))]
    elseif theme.key_position == :left
        position = left_guide_position
    elseif theme.key_position == :top
        position = top_guide_position
    elseif theme.key_position == :bottom
        position = bottom_guide_position
    end

    return [Guide.PositionedGuide(ctxs, 0, position)]
end

function render_discrete_key(labels::Vector{String}, title_ctx::Context, title_width::Measure, theme::Gadfly.Theme;
    colors=[nothing], aes_color_label=nothing, shapes=[nothing], sizes=[nothing])

    n = max(length(colors), length(shapes), length(sizes))
    shape1 = shapes[1]
    shapes = (shape1==nothing) ? fill(theme.key_swatch_shape, n) : theme.point_shapes[shapes]
    (colors[1]==nothing) && (colors = fill((theme.key_swatch_color==nothing) ? theme.default_color : theme.key_swatch_color, n))
    (sizes[1]==nothing) && (sizes = fill((theme.key_swatch_size==nothing) ? theme.point_size : theme.key_swatch_size, n))

    # only consider layouts with a reasonable number of columns
    maxcols = theme.key_max_columns < 1 ? 1 : theme.key_max_columns
    maxcols = min(n, maxcols)

    extents = text_extents(theme.key_label_font, theme.key_label_font_size, values(labels)...)
    text_widths, text_heights = first.(extents), last.(extents)

    ypad = 1.0mm
    title_height = title_ctx.box.a[2]
    swatch_size = 2*maximum(sizes)
    entry_height = max(swatch_size, maximum(text_heights)) + ypad

    # return a context with a lyout of numcols columns
    function make_layout(numcols)
        colrows = Array{Int}(undef, numcols)
        m = n
        for i in 1:numcols
            colrows[i] = ceil(Int, m/(1+numcols-i))
            m -= colrows[i]
        end

        xpad = 1mm
        colwidths = Array{Measure}(undef, numcols)
        m = 0
        for (i, nrows) in enumerate(colrows)
            if m == n
                colwidths[i] = 0mm
            else
                colwidth = maximum(text_widths[m+1:m+nrows])
                colwidth += entry_height + xpad
                colwidths[i] = colwidth
                m += nrows
            end
        end

        ctxwidth = sum(colwidths)
        ctxheight = entry_height * colrows[1] + title_height

        ctxp = ctxpromise() do draw_context
            yoff = 0.5h - ctxheight/2
            outerctx = context()

            compose!(outerctx, (context(xpad, yoff), title_ctx))

            ctx = context(0, yoff+title_height, ctxwidth, ctxheight-title_height, units=UnitBox(0, 0, 1, colrows[1]))

            m = 0
            xpos = 0w
            for (colwidth, nrows) in zip(colwidths, colrows)

                x = [0.5cy]
                clrs = colors[m+1:m+nrows]
                shps = shapes[m+1:m+nrows]
                szs = sizes[m+1:m+nrows]

                swatches_shapes = [f(x, [y-0.5].*cy, [s]) for (y,(f,s)) in enumerate(zip(shps, szs))]
                sw1 = [(context(), s, fill(c), fillopacity(theme.alphas[1]), stroke(theme.discrete_highlight_color(c)))
                     for (s,c) in zip(swatches_shapes, clrs)]
                swatches = compose!(context(), linewidth(theme.highlight_width), sw1...)

                swatch_labels = compose!(context(),
                text([entry_height+xpad], [0.5:nrows;]*cy, labels[m+1:m+nrows], [hleft], [vcenter]),
                    font(theme.key_label_font),
                    fontsize(theme.key_label_font_size),
                    fill(theme.key_label_color))

                col = compose!(context(xpos, 0, colwidth), swatches, swatch_labels)

                if aes_color_label != nothing
                    classes = [svg_color_class_from_label(aes_color_label([c])[1]) for c in clrs]
                    #class_jscalls = ["data(\"color_class\", \"$(c)\")" for c in classes]
                    compose!(col,
                        svgclass(classes),
                        jscall(["""
                            data(\"color_class\", \"$(c)\")
                            .click(Gadfly.colorkey_swatch_click)
                            """ for c in classes]))
                end
                compose!(ctx, col)

                m += nrows
                xpos += colwidth
            end

            return compose!(outerctx, ctx,
                            # defeat webkit's asinine default drag behavior
                            jscall("drag(function() {}, function() {}, function() {})"),
                            svgclass("guide colorkey"))
        end

        return compose!(context(minwidth=max(title_width, ctxwidth),
                    minheight=ctxheight, units=UnitBox()),
                    ctxp)
    end

    return map(make_layout, 1:maxcols)
end

function render_key_title2(title::AbstractString, theme::Gadfly.Theme)
    title_width, title_height = max_text_extents(theme.key_title_font, theme.key_title_font_size, title)

    if theme.guide_title_position == :left
        title_form = text(0.0w, title_height, title, hleft, vbottom)
    elseif theme.guide_title_position == :center
        title_form = text(0.5w, title_height, title, hcenter, vbottom)
    elseif theme.guide_title_position == :right
        title_form = text(1.0w, title_height, title, hright, vbottom)
    else
        error("$(theme.guide_title_position) is not a valid guide title position")
    end

    title_padding = 1.5mm
    title_context = compose!(
        context(0w, 0h, 1w, title_height + title_padding),
        title_form,
        stroke(nothing),
        font(theme.key_title_font),
        fontsize(theme.key_title_font_size),
        fill(theme.key_title_color))

    return title_context, title_width
end

struct SizeKey <: Gadfly.GuideElement
    title::AbstractString
    labels::Vector{String}
    pos::Vector{Compose.MeasureOrNumber}
    visible::Bool
end
SizeKey(;title="Size", labels=String[], pos=[], visible=true) = SizeKey(title, labels, pos, visible)
SizeKey(v::Nothing) = SizeKey(visible=false)
SizeKey(title::AbstractString, labels::Vector{String}, pos::Vector) = SizeKey(title, labels, pos, true)

"""
    Guide.sizekey[(; title="size", labels=String[], pos=[])]
    Guide.sizekey(title, labels, pos)

Enable control of the sizekey.  Set the key `title` and the item `labels`.
`pos` overrides [Theme(key_position=)](@ref Gadfly) and can be in either
relative (e.g. [0.7w, 0.2h] is the lower right quadrant), absolute (e.g. [0mm,
0mm]), or plot scale (e.g. [0,0]) coordinates. Currently `Guide.sizekey` will only work by adding
`Scale.size_discrete2` to the plot. `Guide.sizekey(nothing)` will hide the key.
"""
const sizekey = SizeKey

function render(guide::SizeKey, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)

    (theme.key_position==:none || !guide.visible || aes.size_key_vals===nothing) && return PositionedGuide[]
    gpos = guide.pos
    (theme.key_position==:inside) && isempty(gpos) &&  (gpos = [0.7w, 0.25h])

    # Aesthetics for keys: size_key_title, size_label (Function), size_key_vals (AbstractDict)
    nsizes = length(unique(aes.size))
    guide_title = (guide.title≠"Size" || aes.size_key_title===nothing) ? guide.title : aes.size_key_title
    sizes = collect(keys(aes.size_key_vals))
    size_key_labels = isempty(guide.labels) ? aes.size_label(sizes) : guide.labels

    colors = [nothing]
    if (aes.size_key_title≠nothing)  && (aes.color_key_title==aes.size_key_title)
        colors = collect(keys(aes.color_key_colors))
    end

    title_context, title_width = render_key_title2(guide_title, theme)
    ctxs = render_discrete_key(size_key_labels, title_context, title_width, theme, sizes=sizes, colors=colors)

    position = right_guide_position
    if !isempty(gpos)
        position = over_guide_position
        ctxs = [compose(context(), (context(gpos[1], gpos[2]), ctxs[1]))]
    elseif theme.key_position == :left
        position = left_guide_position
    elseif theme.key_position == :top
        position = top_guide_position
    elseif theme.key_position == :bottom
        position = bottom_guide_position
    end

    return [PositionedGuide(ctxs, 0, position)]
end

ext: jl


text:
module LARS

using LinearAlgebra, StatsBase

end # module


text:
using MLStyle

_wrap_block(ex::Expr) = ex
_wrap_block(o) = Expr(:block, o)
_simplify_block(ex) = _wrap_block(simplify_ex(ex))

function simplify_ex(ex::Expr)
    @match ex begin
        Expr(:for, Expr(:block, assigns...), body) =>
            foldr(assigns, init=_simplify_block(body)) do ass, last
                Expr(:for, simplify_ex(ass), last)
            end
        Expr(:let, Expr(:block, assigns...), body) =>
            foldr(assigns, init=_simplify_block(body)) do ass, last
                Expr(:let, simplify_ex(ass), last)
            end
        Expr(:(=), lhs && Expr(:call, _...), rhs) =>
            Expr(:function, simplify_ex(lhs), simplify_ex(rhs))
        Expr(:->, lhs && Expr(:tuple, _...), rhs) =>
            Expr(:function, simplify_ex(lhs), simplify_ex(rhs))
        Expr(:(->), lhs, rhs) =>
            Expr(:function, Expr(:tuple, simplify_ex(lhs)), simplify_ex(rhs))
        Expr(:do, Expr(:call, f, args...), lam) =>
            Expr(:call, simplify_ex(f), simplify_ex(lam), map(simplify_ex, args)...)
        Expr(head, args...) =>
            Expr(head, map(simplify_ex, args)...)
    end
end

simplify_ex(s) = s


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

################################################################################
#
#  Independence test
#
################################################################################

# see Robledo, p. 47
@doc Markdown.doc"""
    isindependent(S::Array{EllCrvPt{fmpq}}) -> Bool

Tests whether a given set of points $S$ on a rational elliptic curve
is linearly independent. Returns true if they are independent, otherwise false.
This function may return false results.
"""
function isindependent(P)
  epsilon = 10.0^(-8)
  r = length(P)
  M = Matrix{Float64}(r,r)

  for i in 1:r
    for j in 1:r
      M[i,j] = canonical_height(P[i] + P[j]) - canonical_height(P[i]) - canonical_height(P[j])
    end
  end

  determinante = det(M)

  if abs(determinante) < epsilon
    return false
  else
    return true
  end
end

################################################################################
#
#  Arithmetic geometric mean
#
################################################################################

@doc Markdown.doc"""
    agm(x::Float64, y::Float64, e::Int) -> Float64
  Returns the arithmetic-geometric mean of $x$ and $y$.
"""
function agm(x::Float64, y::Float64, e::Int = 5)
    0 < y && 0 < y && 0 < e || throw(DomainError())
    err = e*eps(x)
    (g, a) = extrema([x, y])
    while err < (a - g)
        ap = a
        a = 0.5*(a + g)
        g = sqrt(ap*g)
    end
    return a
end

################################################################################
#
#  Real period
#
################################################################################

# see Cohen
@doc Markdown.doc"""
    real_period(E::EllCrv{fmpz}) -> Float64
  Returns the real period of an elliptic curve $E$ with integer coefficients.
"""
function real_period(E)
  a1 = numerator(E.coeff[1])
  a2 = numerator(E.coeff[2])
  a3 = numerator(E.coeff[3])
  a4 = numerator(E.coeff[4])
  a6 = numerator(E.coeff[5])

  b2, b4, b6, b8 = get_b_integral(E)

  delta = discriminant(E)
  f(x) = x^3 + (Float64(b2)/4)*x^2 + (Float64(b4)/2)*x + (Float64(b6)/4)
  root = fzeros(f)

  if delta < 0 # only one real root
    e1 = root[1]
    a = 3*e1 + (Float64(b2)/4)
    b = sqrt(3*e1^2 + (Float64(b2)/2)*e1 + (Float64(b4)/2))
    lambda = 2*pi / agm(2*sqrt(b), sqrt(2*b + a))
  else
    root = sort(root)
    e1 = root[1]
    e2 = root[2]
    e3 = root[3]
    w1 = pi / agm(sqrt(e3-e1), sqrt(e3-e2))
    lambda = 2*w1
  end

  return lambda
end

################################################################################
#
#  Logarithmic height of a rational number
#
################################################################################

@doc Markdown.doc"""
    height(x::fmpq) -> Float64

Computes the height of a rational number $x$.
"""
function log_height(x::fmpq)
  a = Float64(numerator(x))
  b = Float64(denominator(x))
  return log(max(abs(a),abs(b)))
end

# every rational point is given by P = (a/c^2, b/c^3), gcd(a,c) = gcd(b,c) = 1. then h(P) = max(|a|, c^2)
# but is it always in this form?
@doc Markdown.doc"""
    naive_height(P::EllCrvPt{fmpq}) -> Float64

Computes the naive height of a point $P$.
"""
function naive_height(P)
  x = P.coordx
  a = Float64(numerator(x))
  c2 = Float64(denominator(x))

  h = log(max(abs(a), abs(c2)))
  return h
end

# p.75 Cremona
@doc Markdown.doc"""
    points_with_bounded_naive_height(E:EllCrv, B::Int) -> Array{EllCrvPt}

Computes all rational points on a curve $E$ with integer coefficients which have naive height <= $B$.
"""
function points_with_bounded_naive_height(E, B)
  a1 = numerator(E.coeff[1])
  a2 = numerator(E.coeff[2])
  a3 = numerator(E.coeff[3])
  a4 = numerator(E.coeff[4])
  a6 = numerator(E.coeff[5])

  # 2-torsion
  f(x) = x^3 + Float64(a2)*x^2 + Float64(a4)*x + Float64(a6)
  torsiontwo = sort(fzeros(f))
  x0 = torsiontwo[1]

  R, z = PolynomialRing(FlintZZ, "z")

  points = []

  # iterate over possible values for c and a
  k = Int(floor(exp(Float64(B)/2)))
  for c in 1:k
    for a in Int(floor(max(c^2*x0, -exp(Float64(B))))):Int(floor(exp(Float64(B))))
      if gcd(a, c) == 1
      # look for possible values for b; they are the zeros of g
        g = z^2 + (a1*c*a + a3*c^3)*z - (a^3 + a2*c^2*a^2 + a4*c^4*a + a6*c^6)
        zeros = zeros(g)

        if length(zeros) != 0
          for b in zeros
            P = E([FlintQQ(a,c^2), FlintQQ(b, c^3)])
            push!(points, P)
          end
        end
      end
    end
  end

  return points
end

@doc Markdown.doc"""
torsion_points_via_height(E::EllCrv{fmpz}) ->  Array{EllCrvPt}
Returns the rational torsion points of a curve $E$ with integer coefficients.
"""
function torsion_points_via_height(E::EllCrv{fmpq})

  if E.short == true
    E = EllipticCurve([0, 0, 0, E.coeff[1], E.coeff[2]])
  end

  jay = j(E)
  hj = log_height(jay) # height of the j-invariant
  jfloat = Float64(numerator(jay))/Float64(denominator(jay))

  delta = numerator(discriminant(E))
  b2, b4, b6, b8 = get_b_integral(E)
  twostar = 2
  if b2 == 0
    twostar = 1
  end

  # mu(E)
  mu = (1/6)*( log(abs(Float64(delta))) + log(max(1, abs(jfloat))) ) + log(max(1, abs(Float64(b2)/12))) + log(twostar)

  B = (1/12)*hj + mu + 1.922

  # all torsion points have naive height <= B, see Cremona, p. 77
  torsion_candidates = points_with_bounded_naive_height(E, B)
  torsion_points = [infinity(E)]

  # check which points of the candidates are torsion points (order of a torsion point is <= 12)
  for P in torsion_candidates
    istorsion = false
    k = 7
    while (istorsion == false) && (k <= 12)
      if (k*P).isinfinite == true
        istorsion = true
      end
      k = k + 1
    end

    if istorsion == true
      push!(torsion_points, P)
    end
  end

  return torsion_points
end

@doc Markdown.doc"""
independent_points_up_to(E::EllCrv{fmpq}, B::Int) -> Array{EllCrvPt}
Returns a maximal set of independent points with naive height <= $B$
"""
function independent_points_up_to(E::EllCrv{fmpq},B::IntegerUnion)

  if E.short == true
      E = EllipticCurve([0, 0, 0, E.coeff[1], E.coeff[2]])
  end

  points = points_with_bounded_naive_height(E,B)
  counter = 1
  M_ind = Matrix{Float64}(0,0)
  M_cand = Matrix{Float64}(1,1)
  points_ind = []

  for p in points
    istorsion = false
    i = 7
    while (i < 13) && (istorsion == false)
      if (i*p).isinfinite == true
        istorsion = true
      end
      i = i + 1
    end

    if istorsion == true
      continue
    end

    push!(points_ind, p)
    for i = 1:length(points_ind)
      M_cand[i, counter] = canonical_height(points_ind[i] + points_ind[counter]) - canonical_height(points_ind[i]) - canonical_height(points_ind[counter])
      M_cand[counter, i] = M_cand[i, counter]
    end

    if abs(det(M_cand)) < 10.0^(-8)
      pop!(points_ind)
    else
      counter = counter + 1
      M_ind = M_cand
      M_cand = Matrix{Float64}(size(M_cand)[1] + 1, size(M_cand)[1] + 1)

      for i = 1:size(M_cand)[1] - 1
        for j = 1:size(M_cand)[1] - 1
          M_cand[i, j] = M_ind[i, j]
        end
      end
    end
  end

  return points_ind
end

################################################################################
#
#  BSD
#
################################################################################

@doc Markdown.doc"""
mod_red(E::EllCrv, B::Int) -> (P::Array{Int}, N::Array{Nemo.fmpz})
input: E elliptic curve given in long form over ZZ
output: arrays P, N, where
  P contains all primes smaller than B (for which E/F_p is non-singular)
  N[i] = #E(F_P[i])
"""
function mod_red(E, B)
    minmodel = EllipticCurve(laska_kraus_connell(E)) # global minimal model for E
    P = primes(B) # primes smaller than B
    N = Nemo.fmpz[]

    for i in 1:length(P)
        p = P[i]
        R = GF(p, cached = false)
        Ep = EllipticCurve([R(numerator(minmodel.coeff[1])), R(numerator(minmodel.coeff[2])), R(numerator(minmodel.coeff[3])), R(numerator(minmodel.coeff[4])), R(numerator(minmodel.coeff[5]))],  false) # reduction of E mod p

        if  discriminant(Ep) != 0 # can only determine group order if curve is non-singular
            ord = order_best(Ep)
            push!(N, ord)
        else
            P[i] = 0
        end
    end

    P = deleteat!(P, findin(P, 0)) # delete all zeros from P

    return P, N
end

@doc Markdown.doc"""
check_weak_bsd(E::EllCrv, B::Int) -> (a::Float64, b::Float64)
checks weak bsd-conjecture for elliptic curve $E$ given in long form over ZZ, positive integer $B$
returns linear regression values for log(log(B)) and sum of log(N_p/p) for p <= B
"""
function check_weak_bsd(E, B)

    (P, N) = mod_red(E, B)
    a = length(P)
    logprod = Float64[]
    loglogB = Float64[]

    # initial value
    push!(logprod, log(Int(N[1]) / P[1]) ) # N is nemo.fmpz, P is Int64
    push!(loglogB, log(log( P[1] + 1 )) )

    for i in 2:(a - 1)
        push!(logprod, log(Int(N[i]) / P[i]) + logprod[i-1] )
        push!(loglogB, log(log( float(P[i] + 1 ))) )
    end

    # last value
    push!(logprod, log(Int(N[a]) / P[a]) + logprod[a-1] )
    push!(loglogB, log(log(B)) )

    a, b = linreg(loglogB, logprod)
    plot(loglogB, logprod, "o")
    plot(loglogB, [a + b*i for i in loglogB])

    return a, b
end


text:
function parse_line(l)
    name, spd, dur, rst = match(r"(\w+)[ \D]+(\d+)[ \D]+(\d+)[ \D]+(\d+).*", l).captures
    (name, parse(Int, spd), parse(Int, dur), parse(Int, rst))
end

const data = parse_line.(readlines())
const time = 2503

function distance(spd, dur, rst, time)
    dst = zeros(Int, time+1)
    run = dur
    slp = 0
    for i in 2:time+1
        if run > 0
            dst[i] = dst[i-1] + spd
            run -= 1
            run == 0 && (slp = rst)
        else
            dst[i] = dst[i-1]
            slp -= 1
            slp == 0 && (run = dur)
        end
    end
    # @show dst
    dst
end

# Part 1
maximum(data) do (_, spd, dur, rst)
    distance(spd, dur, rst, time)[end]
end |> println

all_scores = reduce(hcat, [distance(spd, dur, rst, time) for (_, spd, dur, rst) in data])
best_score = maximum(all_scores; dims=2)

# Part 2
maximum(data) do (_, spd, dur, rst)
    (distance(spd, dur, rst, time) .== best_score) |> sum
end - 1 |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
#=

# 226: Terminal flux calculation via test functions, nD, boundary reaction
([source code](SOURCE_URL))

=#

module Example226_BoundaryIntegral

using VoronoiFVM, GridVisualize,ExtendableGrids

function main(;n=10,Plotter=nothing,verbose=false, unknown_storage=:sparse,dim=2)
    n=[101,21,5]
    X=collect(range(0.0,1,length=n[dim]))
    if dim==1
        grid=simplexgrid(X)
        Γ_where_T_equal_1=[2]
        Γ_where_T_equal_0=[1]
    elseif dim==2
        grid=simplexgrid(X,X)
        Γ_where_T_equal_1=[2]
        Γ_where_T_equal_0=[4]
    elseif dim==3
        grid=simplexgrid(X,X,X)
        Γ_where_T_equal_1=[2]
        Γ_where_T_equal_0=[4]
    end

    function storage(f,u,node)
        f.=u
    end

    function flux(f,u,edge)
	f[1]=u[1,1]-u[1,2]
    end

    function breaction(f,u,node)
        if node.region==Γ_where_T_equal_1[1]
            f[1]= u[1]^2
        end
    end

    physics=VoronoiFVM.Physics(flux=flux,
	                       storage=storage,
	                       breaction=breaction)

    system=VoronoiFVM.System(grid,physics)
    enable_species!(system,1,[1])
    boundary_dirichlet!(system,1,Γ_where_T_equal_0[1],1.0);
    inival=unknowns(system,inival=0.0)

    U=solve(inival,system)

    tf=VoronoiFVM.TestFunctionFactory(system)
    T=testfunction(tf,Γ_where_T_equal_0,Γ_where_T_equal_1)

    scalarplot(grid,U[1,:],Plotter=Plotter,zplane=0.50001)
    I=integrate(system,T,U)
    B=integrate(system,breaction,U; boundary=true)
    isapprox(-I[1], B[Γ_where_T_equal_1[1]],rtol=1.0e-12)
end

function test()
    main(dim=1, unknown_storage=:sparse )  ? true : return false
    main(dim=1, unknown_storage=:dense  )  ? true : return false
    main(dim=2, unknown_storage=:sparse )  ? true : return false
    main(dim=2, unknown_storage=:dense  )  ? true : return false
    main(dim=3, unknown_storage=:sparse )  ? true : return false
    main(dim=3, unknown_storage=:dense  )  ? true : return false
end

end


text:
"""
    BoseHubbardMom1D2C(add::BoseFS2C; ua=1.0, ub=1.0, ta=1.0, tb=1.0, v=1.0, kwargs...)

Implements a one-dimensional Bose Hubbard chain in momentum space with a two-component
Bose gas.

```math
\\hat{H} = \\hat{H}_a + \\hat{H}_b + \\frac{V}{M}\\sum_{kpqr} b^†_{r} a^†_{q} b_p a_k δ_{r+q,p+k}
```

# Arguments

* `add`: the starting address.
* `ua`: the `u` parameter for Hamiltonian a.
* `ub`: the `u` parameter for Hamiltonian b.
* `ta`: the `t` parameter for Hamiltonian a.
* `tb`: the `t` parameter for Hamiltonian b.
* `v`: the inter-species interaction parameter V.
Further keyword arguments are passed on to the constructor of [`HubbardMom1D`](@ref).

# See also

* [`BoseFS2C`](@ref)
* [`BoseHubbardReal1D2C`](@ref)

"""
struct BoseHubbardMom1D2C{T,HA,HB,V} <: TwoComponentHamiltonian{T}
    ha::HA
    hb::HB
end

function BoseHubbardMom1D2C(add::BoseFS2C; ua=1.0, ub=1.0, ta=1.0, tb=1.0, v=1.0, args...)
    ha = HubbardMom1D(add.bsa; u=ua, t=ta, args...)
    hb = HubbardMom1D(add.bsb; u=ub, t=tb, args...)
    T = promote_type(eltype(ha), eltype(hb))
    return BoseHubbardMom1D2C{T,typeof(ha),typeof(hb),v}(ha, hb)
end

function Base.show(io::IO, h::BoseHubbardMom1D2C)
    addr = starting_address(h)
    ua = h.ha.u
    ub = h.hb.u
    ta = h.ha.t
    tb = h.hb.t
    v = h.v
    print(io, "BoseHubbardMom1D2C($addr; ua=$ua, ub=$ub, ta=$ta, tb=$tb, v=$v)")
end

function starting_address(h::BoseHubbardMom1D2C)
    return BoseFS2C(starting_address(h.ha), starting_address(h.hb))
end

LOStructure(::Type{<:BoseHubbardMom1D2C{<:Real}}) = IsHermitian()

Base.getproperty(h::BoseHubbardMom1D2C, s::Symbol) = getproperty(h, Val(s))
Base.getproperty(h::BoseHubbardMom1D2C, ::Val{:ha}) = getfield(h, :ha)
Base.getproperty(h::BoseHubbardMom1D2C, ::Val{:hb}) = getfield(h, :hb)
Base.getproperty(h::BoseHubbardMom1D2C{<:Any,<:Any,<:Any,V}, ::Val{:v}) where {V} = V

function num_offdiagonals(ham::BoseHubbardMom1D2C, add::BoseFS2C)
    M = num_modes(add)
    sa = num_occupied_modes(add.bsa)
    sb = num_occupied_modes(add.bsb)
    return num_offdiagonals(ham.ha, add.bsa) + num_offdiagonals(ham.hb, add.bsb) + sa*(M-1)*sb
    # number of excitations that can be made
end

function diagonal_element(ham::BoseHubbardMom1D2C, add::BoseFS2C)
    M = num_modes(add)
    onrep_a = onr(add.bsa)
    onrep_b = onr(add.bsb)
    interaction2c = Int32(0)
    for p in 1:M
        iszero(onrep_b[p]) && continue
        for k in 1:M
            interaction2c += onrep_a[k]*onrep_b[p] # b†_p b_p a†_k a_k
        end
    end
    return (
        diagonal_element(ham.ha, add.bsa) +
        diagonal_element(ham.hb, add.bsb) +
        ham.v/M*interaction2c
    )
end

function get_offdiagonal(ham::BoseHubbardMom1D2C, add::BoseFS2C, chosen)
    return offdiagonals(ham, add)[chosen]
end

"""
    OffdiagonalsBoseMom1D2C

Specialized [`AbstractOffdiagonals`](@ref) that keep track of number of off-diagonals and
number of occupied sites in both components of the address.
"""
struct OffdiagonalsBoseMom1D2C{
    A<:BoseFS2C,T,H<:TwoComponentHamiltonian{T},O1,O2,M1,M2
} <: AbstractOffdiagonals{A,T}
    hamiltonian::H
    address::A
    length::Int
    offdiags_a::O1
    offdiags_b::O2
    map_a::M1
    map_b::M2
end

function offdiagonals(h::BoseHubbardMom1D2C{T}, a::BoseFS2C) where {T}
    offdiags_a = offdiagonals(h.ha, a.bsa)
    offdiags_b = offdiagonals(h.hb, a.bsb)
    map_a = OccupiedModeMap(a.bsa)
    map_b = OccupiedModeMap(a.bsb)
    occ_a = length(map_a)
    occ_b = length(map_b)
    len = length(offdiags_a) + length(offdiags_b) + occ_a * (num_modes(a) - 1) * occ_b

    return OffdiagonalsBoseMom1D2C(h, a, len, offdiags_a, offdiags_b, map_a, map_b)
end

function Base.getindex(s::OffdiagonalsBoseMom1D2C{A,T}, i)::Tuple{A,T} where {A,T}
    @boundscheck 1 ≤ i ≤ s.length || throw(BoundsError(s, i))
    num_hops_a = length(s.offdiags_a)
    num_hops_b = length(s.offdiags_b)
    if i ≤ num_hops_a
        new_a, matrix_element = s.offdiags_a[i]
        new_add = A(new_a, s.address.bsb)
    elseif i ≤ num_hops_a + num_hops_b
        i -= num_hops_a
        new_b, matrix_element = s.offdiags_b[i]
        new_add = A(s.address.bsa, new_b)
    else
        i -= num_hops_a + num_hops_b
        new_a, new_b, val = momentum_transfer_excitation(
            s.address.bsa, s.address.bsb, i, s.map_a, s.map_b
        )
        new_add = A(new_a, new_b)
        matrix_element = s.hamiltonian.v/num_modes(new_add) * val
    end
    return new_add, matrix_element
end

Base.size(s::OffdiagonalsBoseMom1D2C) = (s.length,)


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

statement_1 = "This is cat "
# This is string
println("$statement_1: This is String")
println("##############################################################")

character_1 = 'h'
println("$character_1:This is Character")

println("##############################################################")

statement_3 = """This is also string dont be shocked"""

println(statement_3)
# character_2 = 'he'
# println("$character_2:This is not a Character Error will be raised so I'm commenting ")
#=
Let's do some string concatination
=#

statement_2 = "Are we cats "

println("#######################  1st-Method  ###########################")
println("$statement_1$statement_2")

println("#######################  2nd-Method(bit different lets do some magic here )  ###########################")

println("$statement_1, Are we ",10," cats")

println("#######################  3rd method(string interpolation)  ###########################")

println("$statement_1"*"$statement_2")


text:
using SpatialGrids
using StaticArrays
using Test

include("spatial_grid_test.jl")


text:
using Test
using V5RPC
using UUIDs

@testset "V5Packet" begin
    d = Vector{UInt8}("Hello")
    i = uuid4()
    p = V5RPC.V5Packet(d, i)
    @test p.REPLY == true
    p.flags = 0x0
    @test p.REPLY == false
    io = IOBuffer()
    write(io, p)
    seekstart(io)
    q = read(io, V5RPC.V5Packet)
    @test p.requestid == q.requestid
    @test p.payload == q.payload
end #testset


text:
#!/usr/bin/env julia

## This rootfs includes enough of a host toolchain to build the LLVM passes.
## Eventually, this image will probably be replaced with the actual builder image,
## as that will have the necessary toolchains as well, but that image is not built yet.

include("rootfs_utils.jl")

# Build debian-based image with the following extra packages:
packages = [
    "bash",
    "build-essential",
    "cmake",
    "curl",
    "gfortran",
    "git",
    "less",
    "libatomic1",
    "m4",
    "perl",
    "pkg-config",
    "python",
    "python3",
    "wget",
]
tarball_path = debootstrap("llvm-passes"; packages)

# Upload it
upload_rootfs_image(tarball_path; tag_name = "v1")

max_stars_repo_name: Luapulu/julia

max_issues_repo_name: Luapulu/julia

max_forks_repo_name: Luapulu/julia


text:
module mmodt6gen4
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: T3block, T6block, T3toT6, T6toT3
using MeshSteward: vconnected, vnewnumbering, compactify
using Test
function test()
    connectivity = T6block(2.0, 0.75*pi, 6, 5, :b)
    locs = connectivity.right.attributes["geom"]
    for i in 1:length(locs)
        r, a = locs[i][1]+2.7, locs[i][2]
        locs[i] = (cos(a)*r, sin(a)*r)
    end
    @test nshapes(connectivity.left) == 60
    connectivity3 = T6toT3(connectivity)
    @test nshapes(connectivity3.left) == nshapes(connectivity.left)
    @test nshapes(connectivity3.right) == 42

    vtkwrite("mmodt6gen4", connectivity3)
    try rm("mmodt6gen4.vtu"); catch end
    true
end
end
using .mmodt6gen4
mmodt6gen4.test()

module mmodt6gen5
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite, summary
using MeshSteward: Q4block, transform, fusevertices, cat, renumbered, withvertices
using MeshSteward: vconnected, vnewnumbering, compactify
using Test
function test()
    conn1 = Q4block(2.0, 0.75*pi, 13, 12)
    conn2 = Q4block(3.0, 0.75*pi, 7, 12)
    transform(conn2, x -> [x[1]+2, x[2]])
    # vtkwrite("mmodt6gen5-1", conn1)
    # vtkwrite("mmodt6gen5-2", conn2)
    locs1 = conn1.right.attributes["geom"]
    locs2 = conn2.right.attributes["geom"]
    tolerance = 1.0e-3
    nlocs1, ni1 = fusevertices(locs1, locs2, tolerance)
    conn1 = withvertices(conn1, nlocs1)
    conn2 = withvertices(conn2, nlocs1)
    conn1 = renumbered(conn1, ni1)
    @test summary(conn1) == "(elements, vertices): elements = 156 x Q4, vertices = 273 x P1 {geom,}"
    @test summary(conn2) == "(elements, vertices): elements = 84 x Q4, vertices = 273 x P1 {geom,}"
    # vtkwrite("mmodt6gen5-1", conn1)
    # vtkwrite("mmodt6gen5-2", conn2)

    connectivity = cat(conn1, conn2)
    locs = connectivity.right.attributes["geom"]
    for i in 1:length(locs)
        r, a = locs[i][1]+2.7, locs[i][2]
        locs[i] = (cos(a)*r, sin(a)*r)
    end

    # vtkwrite("mmodt6gen5", connectivity)
    # try rm("mmodt6gen5.vtu"); catch end
    true
end
end
using .mmodt6gen5
mmodt6gen5.test()

module mmodt6gen6
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite, summary
using MeshSteward: Q4block, transform, fusevertices, cat, renumbered, withvertices
using MeshSteward: vconnected, vnewnumbering, compactify
using Test
function test()
    conn1 = Q4block(2.0, 0.75*pi, 13, 12)
    conn2 = Q4block(3.0, 0.75*pi, 7, 12)
    transform(conn2, x -> [x[1]+2, x[2]])
    # vtkwrite("mmodt6gen6-1", conn1)
    # vtkwrite("mmodt6gen6-2", conn2)
    locs1 = conn1.right.attributes["geom"]
    locs2 = conn2.right.attributes["geom"]
    tolerance = 1.0e-3
    nlocs1, ni1 = fusevertices(locs1, locs2, tolerance)
    conn1 = withvertices(conn1, nlocs1)
    conn2 = withvertices(conn2, nlocs1)
    conn1 = renumbered(conn1, ni1)
    @test summary(conn1) == "(elements, vertices): elements = 156 x Q4, vertices = 273 x P1 {geom,}"
    @test summary(conn2) == "(elements, vertices): elements = 84 x Q4, vertices = 273 x P1 {geom,}"
    # vtkwrite("mmodt6gen6-1", conn1)
    # vtkwrite("mmodt6gen6-2", conn2)

    connectivity = cat(conn1, conn2)
    transform(connectivity, x -> begin
        r, a = x[1]+2.7, x[2]
        [cos(a)*r, sin(a)*r]
    end)
    # locs = connectivity.right.attributes["geom"]
    # for i in 1:length(locs)
    #     r, a = locs[i][1]+2.7, locs[i][2]
    #     locs[i] = (cos(a)*r, sin(a)*r)
    # end

    vtkwrite("mmodt6gen6", connectivity)
    try rm("mmodt6gen6.vtu"); catch end
    true
end
end
using .mmodt6gen6
mmodt6gen6.test()

module mmodt6gen7
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite, summary
using MeshSteward: Q4block, transform, fusevertices, cat, renumbered, withvertices
using MeshSteward: vconnected, vnewnumbering, compactify, mergeirs
using Test
function test()
    conn1 = Q4block(2.0, 0.75*pi, 13, 12)
    conn2 = Q4block(3.0, 0.75*pi, 7, 12)
    transform(conn2, x -> [x[1]+2, x[2]])
    connectivity  = mergeirs(conn1, conn2, 0.001)

    transform(connectivity, x -> begin
        r, a = x[1]+2.7, x[2]
        [cos(a)*r, sin(a)*r]
    end)

    vtkwrite("mmodt6gen7", connectivity)
    try rm("mmodt6gen7.vtu"); catch end
    true
end
end
using .mmodt6gen7
mmodt6gen7.test()

module mt4symrcm1a
using StaticArrays
using MeshCore.Exports
using MeshSteward.Exports
using Test
function test()
    # connectivity = T4block(1.0, 2.0, 3.0, 1, 1, 1, :a)
    connectivity = T4block(1.0, 2.0, 3.0, 7, 9, 13, :a)
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    # @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (1120, 4914)

    connectivity = minimize_profile(connectivity)

    vtkwrite("mt4symrcm1a", connectivity)
    try rm("mt4symrcm1a.vtu"); catch end
    true
end
end
using .mt4symrcm1a
mt4symrcm1a.test()


text:
include("prepareDB.jl")
##
fdata = filter(r ->
    r.Area != "PIR" &&
    r.Area != "ENT" &&
    r.Area != "OT" &&
    r.Areac!=
    # r.Mouse != "FM14",
    r.Valid_frs,
    fulldata)
    # lvA = union(fdata.Area)
    # fdata.Area = categorical(fdata.Area; levels = lvA)
##
function prov_define_bold(df0, signal; av_pres = false, method = :Peak, peak_to = :all, nmax = 4, window = 4:10, nanmath = true)
    fmean, fsem = nanmath ? (safemean, safesem) : (mean, sem)
    if  method == :Max
        #average response over all presentation
        df1 = average_presentation(df0,signal; nanmath = nanmath)
        gd1 = groupby(df1,[:Genotype, :Mouse, :Stim, :Area, :Concentration])
        signal = :Bold
        #find the first n max values in the bold signal
        df2 = combine(gd1, signal => (s -> FMRIOpto.max_bold(s; max = nmax, nanmath = nanmath)) => :Bold)
    elseif method == :Window
        #average all bold values per mouse and concentration in a time window
        df1 = combine(groupby(df0, [:Time,:Concentration,:Area,:Stim,:Mouse,:Genotype]), signal => (s -> mean(skipmissing(s))) => :Bold)
        filter!(r -> r.Time in window, df1)
        gd1 = groupby(df1,[:Genotype, :Mouse, :Stim, :Area, :Concentration])
        df2 = combine(gd1, :Bold => (s -> mean(skipmissing(s))) => :Bold)
    elseif method == :None
        df2 = filter(r -> r.Time in window, df0)
        df2.Bold = df2[:,signal]
    else
        error("unrecognised bold processing method:/n options are :Peak, :Max, :Timewindow")
    end
    return df2
end

function prov_plot_bold(df0; nanmath = true, lim = (-0.4,0.7))
    fmean = nanmath ? safemean : mean
    fsem = nanmath ? safesem : sem
    #calculate mean repsonse over concentrations per area
    df = combine(groupby(df0,[:Area, :Stim, :Concentration, :Genotype]),  :Bold .=> [fmean, fsem] .=> [:mean, :sem])
    df1 = filter(r-> r.Genotype == "sert",df)
    # identify regions with missing concentration info
    res = combine(groupby(df1,[:Area]), :mean => (m -> .!(any(isnan.(m)))) => :to_keep)
    to_keep = filter(r -> r.to_keep, res)[:,:Area]
    # prepare a grid plot with the results for each regions
    df1.Color = [x ? :red : :black for x in df1.Stim]
    # posdict = Dict(v=>p for (v,p) in zip(union(df1.Concentration), collect(1:length(union(df1.Concentration)))))
    # df1.Pos = [get(posdict, x,0) for x in df1.Concentration]
    # sort!(df1,[:Stim,:Pos])
    sort!(df1,[:Stim,:Concentration])
    # fconc && filter!(r -> r.Concentration != 0.0001 ,df1)
    plt = []
        combine(groupby(df1, :Area)) do dd
            pp = @df dd scatter(:Concentration, :mean, yerror = :sem, group = :Stim, color = :Color,
                title = :Area[1], legend = false, ylims = lim)
            @df dd plot!(:Concentration, :mean, color = :Color, group = :Stim)
            Plots.abline!(0,0,color = :black, linestyle = :dash)
            push!(plt,pp)
        end
    return plt, df, to_keep
end
function prov_odour_response(df0; farea = :none)
    # use only control data to determine odour concentration dose dependency
    df1 = filter(r -> !r.Stim, df0)
    #remove lowest concentration
    # fconc && filter!(r -> r.Concentration != 0.0001, df1)
    # filters area with insufficient response, identified from plot_bold function
    farea != :none && filter!(r -> r.Area in farea, df1)
    # identify if the model requires to incorporate time
    if iscolumn(df0,:Time)
        f0 = @formula(Bold ~ 1 + Time + Concentration + (1+Time|Mouse))
        f1 = @formula(Bold ~ 1 + Time + Concentration*Area + (1+Time|Mouse))
    else
        f0 = @formula(Bold ~ 1 + Concentration + (1|Mouse))
        f1 = @formula(Bold ~ 1 + Concentration*Area + (1|Mouse))
    end
    # #replace NaN values used for plotting to missing for model
    # df1.Bold = [isnan(x) ? missing : x for x in df1.Bold]
    # model fit and test
    m0 = fit(MixedModel, f0, df1)
    m1 = fit(MixedModel, f1, df1)
    mp = MixedModels.likelihoodratiotest(m0,m1)
    # identifies area that significantly respond to odour concentration
    res = DataFrame(Area = coefnames(m1), P = m1.pvalues, Z = m1.beta ./ m1.stderror)
    res2 = filter!(r -> occursin("Concentration & Area: ", r.Area), res)
    res2.Area = replace.(res2.Area, "Concentration & Area: " => "")
    to_keep = filter(r -> r.P <= 0.049, res2)[:,:Area]

    return m1, mp, to_keep
end

function prov_stim_interaction(df0, farea)
    # filters area lacking concentration response identified by odour_response function
    df1 = filter(r -> r.Genotype == "sert" && r.Area in farea, df0)
    # fconc && filter!(r -> r.Concentration != 0.0001, df0)
    # identify if the model requires to incorporate time
    if iscolumn(df0,:Time)
        f1 = @formula(Bold ~ 1 + Time + Concentration*Area + Stim + zerocorr(1+Time+Concentration|Mouse) + zerocorr(1+Time+Concentration|Area))
        f2 = @formula(Bold ~ 1 + Time + Concentration + Area + Stim + Concentration & Area + Concentration & Area & Stim + zerocorr(1+Time+Concentration|Mouse) + zerocorr(1+Time+Concentration|Area))
        f3 = @formula(Bold ~ 1 + Time + Concentration*Area*Stim + zerocorr(1+Time+Concentration|Mouse) + zerocorr(1+Time+Concentration|Area))
    else
        f1 = @formula(Bold ~ 1 + Concentration*Area + Stim + zerocorr(1+Concentration|Mouse) + zerocorr(1+Concentration|Area))
        f2 = @formula(Bold ~ 1 + Concentration + Area + Stim + Concentration & Area + Concentration & Area & Stim + zerocorr(1+Concentration|Mouse) + zerocorr(1+Concentration|Area))
        f3 = @formula(Bold ~ 1 + Concentration*Area*Stim + zerocorr(1+Concentration|Mouse) + zerocorr(1+Concentration|Area))
    end
    #replace NaN values used for plotting to missing for model
    # df1.Bold = [isnan(x) ? missing : x for x in df1.Bold]
    m1 = fit(MixedModel, f1,df1)
    m2 = fit(MixedModel, f2,df1)
    m3 = fit(MixedModel, f3,df1)
    p_simple = MixedModels.likelihoodratiotest(m1,m2)
    p_complex = MixedModels.likelihoodratiotest(m2,m3)
    res = model_res(m2)
    filter!(r -> occursin("Stim", r.Area), res)
    sort!(res,:P)
    return df1, m2, m3, p_complex, res
end

function prov_analysis_bold(df0, signal; fconc = true, av_pres = true, method = :Peak, peak_to = :all, nmax = 6, window = 4:12, logarithmic = false, nanmath = true, lim = (-0.4,0.7))
    #adjust bold signal
    bold_df = prov_define_bold(df0,signal; method = method, av_pres = av_pres, peak_to = peak_to, nmax = nmax, window = window, nanmath = nanmath)
    fconc && filter!(r -> r.Concentration != 0.0001, bold_df)
    logarithmic && (bold_df.Concentration = log10.(bold_df.Concentration))
    if method == :None
        gd = groupby(bold_df,[:Genotype, :Mouse, :Stim, :Area, :Concentration])
        pltdf = combine(gd, :Bold => (b -> mean(skipmissing(b))) => :Bold)
    else
        pltdf = transform(bold_df, :Bold => ByRow(b -> (ismissing(b) ? NaN : b)) => :Bold)
    end
    plt, plt_df, to_keep = prov_plot_bold(pltdf; nanmath = true, lim = lim)
    display(plot(plt..., size = (900*1.5,1200*1.5)))#layout = (6,3))
    ## find areas that respond to concentration
    m1, mp, to_keep2 = prov_odour_response(bold_df; farea = to_keep)
    # assess stim concentration interaction
    stim_df, m2, m3, p_complex, res = prov_stim_interaction(bold_df, to_keep2)
    show(m2)
    show(m3)
    show(p_complex)
    return bold_df, plt, plt_df, to_keep, m1, mp, to_keep2, stim_df, m2, m3, p_complex, res
end
##
bold_df, plt, plt_df, to_keep, m1, mp, to_keep2, stim_df, m2, m3, p_complex, res =
    prov_analysis_bold(fdata, :BOLD_90; fconc = true, av_pres = true, method = :Window,
    peak_to = :all, nmax = 6, window = 4:10, logarithmic = true,
    nanmath = true, lim = (-0.4,0.7))
##
union(bold_df.Genotype)
test_df = filter(r -> r.Genotype == "sert" && r.Area in to_keep && r.Area in to_keep2, bold_df)
f1 = @formula(Bold ~ 1 + Concentration * Area + Stim + Stim & Concentration + zerocorr(1+Concentration|Mouse) + zerocorr(1+Concentration|Area))
f2 = @formula(Bold ~ 1 + Concentration * Area * Stim + zerocorr(1+Concentration|Mouse) + zerocorr(1+Concentration|Area))
m1 = fit(MixedModel,f1,test_df)
m2 = fit(MixedModel,f2,test_df)
MixedModels.likelihoodratiotest(m1,m2)
##
f_add = @formula(Bold ~ 1 + Concentration * Area + Stim + Stim & Area + Stim & Concentration + zerocorr(1+Concentration|Mouse) + zerocorr(1+Concentration|Area))
f_mul = @formula(Bold ~ 1 + Concentration * Area + Stim + Stim & Concentration + Stim & Concentration & Area + zerocorr(1+Concentration|Mouse) + zerocorr(1+Concentration|Area))
m_add = fit(MixedModel,f_add,test_df)
m_mul = fit(MixedModel,f_mul,test_df)
##
function AICc(model)
    aic(model) + ((2*(dof(model)^2) + 2*dof(model))/(nobs(model) - dof(model) - 1))
end
function AICc_test(candidate, simpler)
    exp((AICc(candidate) - AICc(simpler))/2)
end

AICc_test(m_add,m_mul)

MixedModels.likelihoodratiotest(m_add,m2)
MixedModels.likelihoodratiotest(m_mul,m2)
##

df00 = fdata
signal = :BOLD_90
fconc = true
av_pres = true
method = :Window
peak_to = :all
nmax = 6
window = 4:10
logarithmic = true

df0 = prov_define_bold(df00, signal; av_pres = av_pres, method = :Window, window = window, nanmath = true)
fconc && filter!(r -> r.Concentration != 0.0001, df0)
logarithmic && (df0.Concentration = log10.(df0.Concentration))
open_html_table(df0)
pltdf = transform(df0, :Bold => ByRow(b -> (ismissing(b) ? NaN : b)) => :Bold)
plt, plt_df, to_keep = prov_plot_bold(pltdf)
plot(plt...)
sum(isnan.(df00.BOLD_90))
sum(ismissing.(df00.BOLD_90))
m1, mp, to_keep2 = prov_odour_response(df0; farea = to_keep)
to_keep2


text:
@testset "Point Arithmetic" begin

    p = Point(1,1)
    @test eltype(p) == Int
    @test_throws MethodError Point(1,1.0)

    p1 = Point(2,2)
    @test p + p1 == Point(3,3)

    p2 = Point(2.0,2.0)
    @test p + p1 == Point(3.0,3.0)

    @test 2*p2 == Point(4.0,4.0)
    @test p2 == Point(2.0,2.0)
    @test p2 != Point(2.0,1.9999999)

    arr = vcat(p2,Point(1.0,1.9999999))
    @test in(p2,arr)

end

@testset "MLine construction from vectors" begin
    x = collect(0:0.1:1)
    y = rand(11)
    m = MLine(x,y)
    @test length(m) == 11
    @test size(m) == (11,)
end

@testset "MLine construction from vector of points" begin
    p1 = Point(1.0,2.3)
    p2 = Point(2.0,2.0)
    p3 = Point(3.0,5.2)
    m2 = MLine([p1;p2;p3])
    @test length(m2) == 3
    @test size(m2) == (3,)
    @test extrema(getx(m2)) == (1.0,3.0)
    @test getx(m2) == [1.0;2.0;3.0]
    @test extrema(gety(m2)) == (2.0,5.2)
    @test eltype(m2) == Point{Float64}
    @test m2[1] == p1
    @test m2[3] == p3
    @test DCEGM.min_x(m2) == 1.0
    @test DCEGM.max_x(m2) == 3.0
    @test DCEGM.min_y(m2) == 2.0
    @test DCEGM.max_y(m2) == 5.2
end

# @testset "Interpolations" begin
#     @testset "interpolate at log(1)" begin
#         x = collect(0:0.1:1)
#         y = log.(1.+x)
#         m = MLine(x,y)
#         i = interp(m,[0.0])
#         @test i[1] == 0.0
#     end
# end

@testset "Modifying MLine methods" begin

    @testset "prepend" begin
        x = collect(0:0.1:1)
        y = log.(1 .+ x)
        L = MLine(x,y)
        prepend!(L,[Point(18.0,-1.1)])
        @test size(L)==(12,)
        @test L[1] == Point(18.0,-1.1)
        @test L[2] == Point(0.0,0.0)
    end

    @testset "append" begin
        x = collect(0:0.1:1)
        y = log.(1 .+ x)
        L = MLine(x,y)
        append!(L,[Point(18.0,-1.1)])
        @test size(L)==(12,)
        @test L[end] == Point(18.0,-1.1)
        @test L[1] == Point(0.0,0.0)
    end

    @testset "delete!" begin
        x = collect(0:0.1:1)
        y = rand(11)
        L = MLine(x,y)
        delete!(L,9)
        @test size(L)==(10,)
        @test L[9] == Point(x[10],y[10])
        @test getx(L) == vcat(x[1:8],x[10:11])
    end
    @testset "insert!" begin
        x = collect(0:0.1:1)
        y = log.(1 .+ x)
        L = MLine(x,y)
        insert!(L,Point(1.1,3.3),9)
        @test size(L)==(12,)
        @test L[9] == Point(1.1,3.3)
        @test L[10] == Point(x[9],y[9])
        @test L[8] == Point(x[8],y[8])
    end
    @testset "splitat" begin
        x = collect(0:0.1:1)
        y = log.(1 .+ x)
        L = MLine(x,y)
        n,o=splitat(L,1)
        @test size(n)==(1,)
        @test size(o)==(11,)

        n,o=splitat(L,2)
        @test size(n)==(2,)
        @test size(o)==(10,)
        @test n[2] == L[2]
        @test o[2] == L[3]
        @test o[1] == L[2]
        @test o[end] == L[end]

        n,o=splitat(L,2,false)
        @test size(n)==(2,)
        @test size(o)==(9,)
        @test n[2] == L[2]
        @test o[2] == L[4]
        @test o[1] == L[3]
    end
    @testset "sortx!" begin
        x = collect(0:0.1:1)
        y = rand(11)
        L = MLine(x,y)
        insert!(L,Point(2.0,2.0),2)
        @test !(issorted(L))
        DCEGM.sortx!(L)
        @test issorted(L)
        @test issorted(getx(L))
    end
    @testset "floor a line's y values" begin
        x = collect(0.0:0.1:1.5)
        L1 = MLine(x,log.(x))
        @test !all(gety(L1) .>= 0.0)
        @test length(L1) == length(x)
        DCEGM.floory!(L1,0.0)
        @test all(gety(L1) .>= 0)
        @test length(L1) == length(x)
    end
end


text:
## completion signal

const completion_signal_base = _packet_offsets[findfirst(x->x==:completion_signal,_packet_names)]

@generated function _completion_signal()
    JuliaContext() do ctx
        T_int8 = LLVM.Int8Type(ctx)
        T_int64 = LLVM.Int64Type(ctx)
        _as = convert(Int, AS.Constant)
        T_ptr_i8 = LLVM.PointerType(T_int8, _as)
        T_ptr_i64 = LLVM.PointerType(T_int64, _as)

        # create function
        llvm_f, _ = create_function(T_int64)
        mod = LLVM.parent(llvm_f)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry", ctx)
            position!(builder, entry)

            # get the kernel dispatch pointer
            intr_typ = LLVM.FunctionType(T_ptr_i8)
            intr = LLVM.Function(mod, "llvm.amdgcn.dispatch.ptr", intr_typ)
            ptr = call!(builder, intr)

            # load the index
            signal_ptr_i8 = inbounds_gep!(builder, ptr, [ConstantInt(completion_signal_base, ctx)])
            signal_ptr = bitcast!(builder, signal_ptr_i8, T_ptr_i64)
            signal = load!(builder, signal_ptr)
            ret!(builder, signal)
        end

        call_function(llvm_f, UInt64)
    end
end

signal_completion(value::Int64) = device_signal_store!(_completion_signal(), value)

## misc. intrinsics
@inline sendmsg(x1, x2=Int32(0)) = ccall("llvm.amdgcn.s.sendmsg", llvmcall, Cvoid, (Int32, Int32), x1, x2)
@inline sendmsghalt(x1, x2=Int32(0)) = ccall("llvm.amdgcn.s.sendmsghalt", llvmcall, Cvoid, (Int32, Int32), x1, x2)
@inline endpgm() = @asmcall("s_endpgm", "", true)


text:
module Spec

using ExprTools: splitdef, combinedef
using Test: @test
using MacroTools: postwalk
using Cassette: Cassette, overdub, prehook, @context

export @test, @validated, @pre_spec, @post_spec

@context ValidationCtx

"""
    @validated code

run `code` in a special context where validation tests are run before and/or after
each registered function is called.

Examples:

    using Spec
    f(x) = √x + 1
    @pre_spec  f(x) = @test x >= 0
    @post_spec f(x) = @test isfinite(__result__)

    julia> @validated f(1)
    2.0

    julia> @validated f(-1)
    Test Failed at REPL[16]:1
      Expression: x >= 0
       Evaluated: -1 >= 0
    ERROR: There was an error during testing

    julia> @validated f(Inf)
    Test Failed at REPL[17]:1
      Expression: isfinite(var"##result#254")
    ERROR: There was an error during testing

    julia> f(Inf)
    Inf

*Any* function encountered during the execution context will have it's pre and/or post
validation tests run; the functions do not need to be at the 'top-level'.

    julia> g(x) = f(x) + 1
    g (generic function with 1 method)

    julia> @validated g(Inf)
    Test Failed at REPL[17]:1
      Expression: isfinite(var"##result#254")
    ERROR: There was an error during testing

    julia> g(Inf)
    Inf

Running code within the `@validated` context *will* impose a performance penalty, but
functions with `@pre_spec`s and/or `@post_spec`s will **not** suffer any performance
penalty outside of a `@validated` context.

    julia> using BenchmarkTools

    julia> h(x) = √x + 2 # never touches f
    h (generic function with 1 method)

    julia> @btime g(($Ref(1))[]);
      1.329 ns (0 allocations: 0 bytes)

    julia> @btime h(($Ref(1))[]);
      1.329 ns (0 allocations: 0 bytes)
"""
macro validated(ex)
    esc(:($Cassette.overdub($ValidationCtx(), () -> $ex)))
end

"""
    @pre_spec fdef

Take a function definition `fdef` and create a spec hook to be run within a
`@validated` context before the specified function gets run.

    using Spec
    f(x) = √x + 1
    @pre_spec  f(x) = @test x >= 0

    julia> @validated f(1)
    2.0

    julia> @validated f(-1)
    Test Failed at REPL[16]:1
      Expression: x >= 0
       Evaluated: -1 >= 0
    ERROR: There was an error during testing
"""
macro pre_spec(fdef)
    d = splitdef(fdef)
    get!(d, :args, [])
    pushfirst!(d[:args], :(::typeof($(d[:name]))))
    pushfirst!(d[:args], :(::$ValidationCtx))
    d[:name] = :($Cassette.prehook)
    (esc ∘ combinedef)(d)
end

"""
    @post_spec fdef

Take a function definition `fdef` and create a spec hook to be run within a
`@validated` context before the specified function gets run. Use `__result__` to
directly reference the result of the function call inside the post_spec.

    using Spec
    f(x) = √x + 1
    @post_spec f(x) = @test isfinite(__result__)

    julia> f(Inf)
    Inf

    julia> @validated f(Inf)
    Test Failed at REPL[17]:1
      Expression: isfinite(var"##result#254")
    ERROR: There was an error during testing
"""
macro post_spec(fdef)
    res = gensym(:result)
    d = splitdef(fdef)
    get!(d, :args, [])
    pushfirst!(d[:args], :(::$ValidationCtx), res, :(::typeof($(d[:name]))))
    d[:name] = :($Cassette.posthook)
    d[:body] = postwalk(d[:body]) do x
        x == :(__result__) ? (res) : x
    end
    (esc ∘ combinedef)(d)
end

end # module


text:
#RSDeltaSigmaPort: Compute power values.
#-------------------------------------------------------------------------------

"""`dbm(v,R=50) = 10*log10(v^2/R*1000)`

The equivalent in dBm of an RMS voltage v
"""
dbm(v; R::Float64=50.0) = 10*log10.(abs.(v .^ 2)/R) .+ 30;

"""`dbp(x) = 10*log10(x)`

The dB equivalent of the power x
"""
dbp(x) = 10*log10.(abs.(x))

"""`dbv(x) = 20*log10(abs(x))`

The dB equivalent of the voltage x
"""
dbv(x) = 20*log10.(abs.(x))

"`y = rms(x; no_dc=false)`"
function rms(x; no_dc::Bool=false)
	if no_dc
		x = x - mean(x)
	end
	return norm(x)/sqrt(length(x))
end

"""`v = undbm(p, z=50) = sqrt(z*10^(p/10-3))`

RMS voltage equivalent to a power p in dBm
"""
undbm(p; z::Float64=50.0) = sqrt.(z*10 .^ (p/10 .- 3))

"""`y = undbp(x)`

Convert x from dB to a power
"""
undbp(x) = 10 .^ (x/10)

"""`y = undbv(x)`

Convert x from dB to a voltage
"""
undbv(x) = 10 .^ (x/20)

#Last line


text:
module DICOMDictionary

using EzXML

function make_dicom_dictionary(outputfile; url = "http://dicom.nema.org/medical/dicom/current/source/docbook/part06/part06.xml", xmlfile = "dicom.xml")
    download(url, xmlfile)
    xml = readxml(xmlfile)
    meta_registry = parse_chapter(xml, "Registry of DICOM File Meta Elements")
    data_registry = parse_chapter(xml, "Registry of DICOM Data Elements")
    open(outputfile, "w+") do io
        write(io, "dcm_dict = Dict(\n")
        write_registry(io, meta_registry)
        write_registry(io, data_registry)
        write(io, ")\n")
    end
    println("DICOM dictionary saved at: $outputfile")
    return nothing
end

function parse_chapter(xml, chaptername)
    # Traverse the XML tree until we reach table
    chapter = find_chapter(xml, chaptername)
    contents = elements(chapter)
    table = find_node(contents, "table")
    tablenodes = elements(table)
    tablebody = find_node(tablenodes, "tbody")
    tablerows = elements(tablebody)
    # Each tablenode represents a DICOM attribute
    registry = []
    for row in tablerows
        entry = parse_tablerow(row)
        if entry.tag !== (0x0000,0x0000) && !isempty(entry.keyword)
            # (0x0000, 0x0000) is failure to parse tag
            # this is because some tags contain 'xx'
            push!(registry, entry)
        end
    end
    return registry
end

function find_chapter(xml, chaptername)
    parents = elements(xml.root)
    chapters = find_nodes(parents, "chapter")
    for chapter in chapters
        children = elements(chapter)
        titles = find_nodes(children, "title")
        for title in titles
            if nodecontent(title) == chaptername
                return chapter
            end
        end
    end
    error("$chaptername not found")
end

find_nodes(nodes, name) = nodes[@. nodename(nodes) == name]
find_node(nodes, name) = only(find_nodes(nodes, name))

function parse_tablerow(nodes::EzXML.Node)
    d = strip.(nodecontent.(elements(nodes)))
    tag = intify(d[1])
    keyword = replace(d[3], "\u200b" => "") # Delete zero-width spaces
    vr = d[4]
    # Some elements can have variable VR, e.g. "US or SS". Pick the first one.
    if length(vr) > 2
        vr = vr[1:2]
    end
    return (; tag = tag, keyword = keyword, vr = vr, vm = d[5])
end

function intify(tag)
    regex = match(r"\(([0-9A-F]{4}),([0-9A-F]{4})\)", tag)
    if isnothing(regex)
        extracted_tag = ["0000", "0000"]
    else
        extracted_tag = regex.captures
    end
    (group, element) = @. parse(UInt16, "0x" * extracted_tag)
    return ((group, element))
end

function write_registry(io, registry)
   for entry in registry
        out = "$(entry.tag) => [:$(entry.keyword), \"$(entry.vr)\", \"$(entry.vm)\"],\n"
        write(io, out)
    end
    return nothing
end

export make_dicom_dictionary

end


text:
export UpdateSampleIterate

"""
    UpdateSampleIterate <: UpdateRule

 Update using by iterating over another recursive rule, while sampling from the input

## Parameters
* `rule::UpdateRule`: rule to use at each iteration (default is Elo() )
* `batch_size::Int`: the size of batches to use (default is 1)
* `n_samples::Int`: the number of random samples to take
* ``

## Not implemented yet
* `batching_rule::` ??? how to select batches ???

```
"""
struct UpdateSampleIterate <: UpdateRule
    rule::UpdateRule
    batch_size::Int
    n_samples::Int
    function UpdateSampleIterate(rule::UpdateRule, batch_size::Int, n_samples::Int)
        @check_args(UpdateSampleIterate, batch_size >= one(batch_size))
        @check_args(UpdateSampleIterate, n_samples >= one(n_samples))
        @check_args(UpdateSampleIterate, update_info(rule)[:state_model] == "recursive")
        new(rule, batch_size, n_samples)
    end
end
UpdateSampleIterate(;rule::UpdateRule=UpdateElo(), batch_size::Int=1, n_samples::Int=1) = UpdateSampleIterate( rule, batch_size, n_samples)

function update_info( rule::UpdateSampleIterate )
    info = Dict(
                :name => "SampleIterate",
                :reference => "none",
                :computation => "sequential",
                :state_model => "recursive",
                :input => update_info( rule.rule )[:input],
                :output => update_info( rule.rule )[:output],
                :model => update_info( rule.rule )[:model],
                :ties => update_info( rule.rule )[:ties],
                :factors => update_info( rule.rule )[:factors],
                :parameters => append!( ["rule(default=Elo)", "batch_size(default=1)"], update_info( rule.rule )[:parameters] ),
                :record => true
                )
    return info
end
update_info( ::Type{UpdateSampleIterate} ) =  error("Need to specify the sub-rule")

function update_ratings( rule::UpdateSampleIterate,
                         input_ratings::RatingsList,
                         input_competitions::DataFrame;
                         record::Union{Missing,RatingsTable}=missing)
    n = size(input_competitions,1)
    m = length( input_ratings.players )
    I = player_indexes( input_ratings.players )
    r = deepcopy(input_ratings)

    # select a batch, and then update using the appropriate rule
    d = input_competitions # just an abbreviation
    for (j,i) = enumerate( 1 : rule.batch_size : rule.n_samples )
        rows = d[ rand(1:n, rule.batch_size), :]
        r = update_ratings(rule.rule, r, rows)
        if !ismissing(record) && size(record,1)>=j
            record[j] = r # record results if required, and space allocated
        elseif !ismissing(record)
            println("Warn: out of space to record ratings: j=$j of $(size(record,1))")
        end
    end

    # output ratings list
    return r
end

function predict_outcome(rule::UpdateSampleIterate,
                         ratingA::Real, ratingB::Real,
                         factorA::Union{Missing,Real}, factorB::Union{Missing,Real})
    return predict_outcome(rule.rule, ratingA, ratingB, factorA, factorB)
end


text:
module TikzGraphs

export plot, Layouts

using Compat

preamble = readstring(joinpath(dirname(@__FILE__), "..", "src", "preamble.tex"))

using TikzPictures
using LightGraphs
using Compat

module Layouts

export Layered, Spring, SimpleNecklace

abstract Layout

type Layered <: Layout
end

type Spring <: Layout
  randomSeed
  Spring(;randomSeed=42) = new(randomSeed)
end

type SimpleNecklace <: Layout
end

end

using .Layouts

plot(g::LightGraphs.SimpleGraph) = plot(g, Layered())
plot{T<:AbstractString}(g::LightGraphs.SimpleGraph, labels::Vector{T}) = plot(g, Layered(), labels)

function plotHelper(g::LightGraphs.SimpleGraph, libraryname::AbstractString, layoutname::AbstractString, options::AbstractString)
  o = IOBuffer()
  println(o, "\\graph [$layoutname, $options] {")
  for e in LightGraphs.edges(g)
    a = src(e)
    b = dst(e)
    println(o, "$a -> $b;")
  end
  # include isolated nodes
  for v in LightGraphs.vertices(g)
    if indegree(g, v) == 0 && outdegree(g, v) == 0
      println(o, "$v;")
    end
  end
  println(o, "};")
  mypreamble = preamble * "\n\\usegdlibrary{$libraryname}"
  TikzPicture(takebuf_string(o), preamble=mypreamble)
end

function plotHelper{T<:AbstractString}(g::LightGraphs.SimpleGraph, libraryname::AbstractString, layoutname::AbstractString, options::AbstractString, labels::Vector{T})
  o = IOBuffer()
  println(o, "\\graph [$layoutname, $options] {")
  for v in LightGraphs.vertices(g)
    println(o, "$v/\"$(labels[v])\",")
  end
  println(o, ";")
  for e in LightGraphs.edges(g)
    a = src(e)
    b = dst(e)
    println(o, "$a -> $b;")
  end
  println(o, "};")
  mypreamble = preamble * "\n\\usegdlibrary{$libraryname}"
  TikzPicture(takebuf_string(o), preamble=mypreamble)
end

function plot(g::LightGraphs.SimpleGraph, p::Layered)
  plotHelper(g, "layered", "layered layout", "")
end

function plot(g::LightGraphs.SimpleGraph, p::Spring)
  options = "random seed = $(p.randomSeed)"
  plotHelper(g, "force", "spring layout", options)
end

function plot(g::LightGraphs.SimpleGraph, p::SimpleNecklace)
  plotHelper(g, "circular", "simple necklace layout", "")
end

function plot{T<:AbstractString}(g::LightGraphs.SimpleGraph, p::Layered, labels::Vector{T})
  plotHelper(g, "layered", "layered layout", "", labels)
end

function plot{T<:AbstractString}(g::LightGraphs.SimpleGraph, p::Spring, labels::Vector{T})
  options = "random seed = $(p.randomSeed)"
  plotHelper(g, "force", "spring layout", options, labels)
end

function plot{T<:AbstractString}(g::LightGraphs.SimpleGraph, p::SimpleNecklace, labels::Vector{T})
  plotHelper(g, "circular", "simple necklace layout", "", labels)
end

end # module


text:
#########
# THIS SECTION JUST STEALS FROM STATSBASE BUT REMOVES TYPE RESTRICTIONS AND ADDS SIMD

# Skewness
# This is Type 1 definition according to Joanes and Gill (1998)
"""
    skew(v, m=mean(v))

Compute the standardized skewness of iterable `v`.
"""
function skew(v, m::T) where {T<:Union{Missing,Real}}
    n = length(v)
    cm2 = 0.0   # empirical 2nd centered moment (variance)
    cm3 = 0.0   # empirical 3rd centered moment
    @simd for v_i in v
        z = v_i - m
        z2 = z * z

        cm2 += z2
        cm3 += z2 * z
    end
    cm3 /= n
    cm2 /= n
    return cm3 / sqrt(cm2 * cm2 * cm2)  # this is much faster than cm2^1.5
end

skew(v) = skew(v,mean(v))

# (excess) Kurtosis
# This is Type 1 definition according to Joanes and Gill (1998)
"""
    kurt(v, m=mean(v))

Compute the excess kurtosis of iterable `v`.
"""
function kurt(v, m::T) where {T<:Union{Missing,Real}}
    n = length(v)
    cm2 = 0.0  # empirical 2nd centered moment (variance)
    cm4 = 0.0  # empirical 4th centered moment
    @simd for v_i in v
        z = v_i - m
        z2 = z * z
        cm2 += z2
        cm4 += z2 * z2
    end
    cm4 /= n
    cm2 /= n
    return (cm4 / (cm2 * cm2)) - 3.0
end

kurt(v) = kurt(v,mean(v))

"""
    autocorrelate(x,n=1;demean=true)

Compute the autocorrelation of lag `n` of `x`. `demean` defaults to true.
"""
function autocorrelate(x,n::Int=1; demean::Bool=true)
    T = typeof(zero(eltype(x)) / 1)
    z::Vector{T} = demean ? x .- mean(x) : x
    zz = dot(z, z)
    ac = @views dot(z[1+n:end],z[1:end-n]) / zz
    return ac
end

"""
    autocorrelate(f,x,n=1;demean=true)

Compute the autocorrelation of lag `n` of `f(x)`. `demean` defaults to true and takes place after `f` is applied.
"""
function autocorrelate(f::Function,x,n::Int=1;demean::Bool=true)
    y = f.(x)
    T = typeof(zero(eltype(y)) / 1)
    z::Vector{T} = demean ? y .- mean(y) : y
    zz = dot(z, z)
    ac = @views dot(z[1+n:end],z[1:end-n]) / zz
    return ac
end

## static kernels are terrible here for long windows, need some branching on n to happen
function varianceratio(x::AbstractVector{T},n::Int) where {T}
    n == 1 && return one(T)
    k = makekernel(sum,-(n-1):0)
    vr = varianceratio(x,k,n)
    return vr
end

function varianceratio(x::AbstractVector{T},k::Kernel,n::Int) where {T}
    v1 = var(x)
    summed = applyrolling(k,x)
    vn = @views var(summed[n:end])
    vr = vn / (n*v1)
    return vr
end

# the below are optimized versions for frequently used horizons
# uses beyond the first will not incur the cost of kernel construction (~250ns)
vr2(x) = varianceratio(x,makekernel(sum,-1:0),2)
vr4(x) = varianceratio(x,makekernel(sum,-3:0),4)
vr6(x) = varianceratio(x,makekernel(sum,-5:0),6)
vr8(x) = varianceratio(x,makekernel(sum,-7:0),8)
vr12(x) = varianceratio(x,makekernel(sum,-11:0),12)

"""
    conditionalcor(f,x,y)

Compute the correlation between `x` and `y` conditional on some event defined by the function `f(x,y)`.
The function `f` must return an object that can index into `x` and `y`.

__**Example**__
```
x = randn(1000)
y = randn(1000)
conditionalcor((xx,yy) -> xx .<= 0,x,y)
conditionalcor((xx,yy) -> (xx .<= 0) .& (yy .<= 0),x,y)
```
"""
function conditionalcor(f,x,y)
    loc = f(x,y)
    out = @views cor(x[loc],y[loc])
    return out
end

downsidecor(x,y;q=0.5) = conditionalcor((xx,yy) -> xx .<= quantile(x,q),x,y)

upsidecor(x,y;q=0.5) = conditionalcor((xx,yy) -> xx .>= quantile(x,q),x,y)

loggrowth(X,n=1; scal=100.0) = scal.*(log.(X) .- lag(log.(X),n))

"""
    correlogram(x,y;leadlags::Int=10)

Create a correlogram from `x` and `y`. By default, uses up to lag 10. Total output is `2*leadlags+1` in size.
If the plot peaks to the left of 0, `x` leads `y`.
"""
function correlogram(x,y;leadlags::Int=10)
    leadlagrange = -leadlags:leadlags
    N = length(leadlagrange)
    cors = zeros(N)
    k = 1
    for i in -leadlags:0
        cors[k] = @views cor(x[1:end+i],y[(-i+1):end])
        k += 1
    end

    for i in 1:leadlags
        cors[k] = @views cor(x[i+1:end],y[1:end-i])
        k += 1
    end
    return cors,leadlagrange
end


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module TestPyFlowGraphs
using Compat
using Test

using Catlab.WiringDiagrams
using SemanticFlowGraphs

using ..IntegrationTest: db, read_raw_graph, read_semantic_graph

# Raw flow graph
################

# Deserialize Python raw flow graph from GraphML.
diagram = read_raw_graph(joinpath("python", "pandas_read_sql"))
@test nboxes(diagram) == 2
b1, b2 = boxes(diagram)
@test isnothing(b1.value.annotation)
@test b2.value.annotation == "python/pandas/read-sql-table"
@test [ p.annotation for p in output_ports(b1) ] ==
  [ "python/sqlalchemy/engine" ]
@test [ p.annotation for p in input_ports(b2)[1:2] ] ==
  [ "python/builtins/str", "python/sqlalchemy/engine" ]
@test [ p.annotation for p in output_ports(b2) ] ==
  [ "python/pandas/data-frame" ]

# Semantic flow graph
#####################

# Read SQL table using pandas and SQLAlchemy.
semantic = read_semantic_graph(joinpath("python", "pandas_read_sql"))
d = WiringDiagram([], concepts(db, ["table"]))
engine = add_box!(d, Box([], concepts(db, ["sql-database"])))
cons = add_box!(d, construct(pair(concept(db, "sql-table-database"),
                                  concept(db, "sql-table-name"))))
read_table = add_box!(d, concept(db, "read-table"))
add_wires!(d, [
  (engine, 1) => (cons, 1),
  (cons, 1) => (read_table, 1),
  (read_table, 1) => (output_id(d), 1),
])
@test semantic == d

# K-means clustering on the Iris dataset using NumPy and SciPy.
semantic = read_semantic_graph(joinpath("python", "scipy_clustering_kmeans"))
kmeans_fit = Hom("fit",
  otimes(concept(db, "k-means"), concept(db, "data")),
  concept(db, "k-means"))
d = WiringDiagram([], concepts(db, ["array","array"]))
file = add_box!(d, construct(pair(concepts(db, ["tabular-file", "filename"])...)))
read_file = add_box!(d, concept(db, "read-tabular-file"))
transform = add_box!(d, Box(concepts(db, ["array"]), concepts(db, ["array"])))
kmeans = add_box!(d, construct(pair(concepts(db,
  ["k-means", "clustering-model-n-clusters"])...)))
fit = add_box!(d, kmeans_fit)
centroids = add_box!(d, concept(db, "k-means-centroids"))
clusters = add_box!(d, concept(db, "clustering-model-clusters"))
add_wires!(d, [
  (file, 1) => (read_file, 1),
  (read_file, 1) => (transform, 1),
  (kmeans, 1) => (fit, 1),
  (transform, 1) => (fit, 2),
  (fit, 1) => (clusters, 1),
  (fit, 1) => (centroids, 1),
  (centroids, 1) => (output_id(d), 1),
  (clusters, 1) => (output_id(d), 2),
])
@test semantic == d

# K-means clustering on the Iris dataset using pandas and scikit-learn.
semantic = read_semantic_graph(joinpath("python", "sklearn_clustering_kmeans"))
d = WiringDiagram([], concepts(db, ["array"]))
file = add_box!(d, construct(pair(concepts(db, ["tabular-file", "filename"])...)))
read_file = add_box!(d, concept(db, "read-tabular-file"))
transform = add_box!(d, Box(concepts(db, ["table"]), concepts(db, ["array"])))
kmeans = add_box!(d, construct(concept(db, "k-means")))
fit = add_box!(d, concept(db, "fit"))
clusters = add_box!(d, concept(db, "clustering-model-clusters"))
add_wires!(d, [
  (file, 1) => (read_file, 1),
  (read_file, 1) => (transform, 1),
  (kmeans, 1) => (fit, 1),
  (transform, 1) => (fit, 2),
  (fit, 1) => (clusters, 1),
  (clusters, 1) => (output_id(d), 1),
])
@test semantic == d

# Compare sklearn clustering models using a cluster similarity metric.
semantic = read_semantic_graph(joinpath("python", "sklearn_clustering_metrics"))
clustering_fit = Hom("fit",
  otimes(concept(db, "clustering-model"), concept(db, "data")),
  concept(db, "clustering-model"))
d = WiringDiagram([], concepts(db, ["array", "k-means", "agglomerative-clustering"]))
make_data = add_box!(d, Box([], concepts(db, ["array", "array"])))
kmeans = add_box!(d, construct(concept(db, "k-means")))
kmeans_fit = add_box!(d, clustering_fit)
kmeans_clusters = add_box!(d, concept(db, "clustering-model-clusters"))
agglom = add_box!(d, construct(concept(db, "agglomerative-clustering")))
agglom_fit = add_box!(d, clustering_fit)
agglom_clusters = add_box!(d, concept(db, "clustering-model-clusters"))
score = add_box!(d, Box(concepts(db, ["array", "array"]), []))
add_wires!(d, [
  (make_data, 1) => (kmeans_fit, 2),
  (make_data, 1) => (agglom_fit, 2),
  (make_data, 2) => (output_id(d), 1),
  (kmeans, 1) => (kmeans_fit, 1),
  (kmeans_fit, 1) => (kmeans_clusters, 1),
  (kmeans_fit, 1) => (output_id(d), 2),
  (kmeans_clusters, 1) => (score, 1),
  (agglom, 1) => (agglom_fit, 1),
  (agglom_fit, 1) => (agglom_clusters, 1),
  (agglom_fit, 1) => (output_id(d), 3),
  (agglom_clusters, 1) => (score, 2),
])
@test semantic == d

# Errors metrics for linear regression using sklearn.
semantic = read_semantic_graph(joinpath("python", "sklearn_regression_metrics"))
d = WiringDiagram([], [])
file = add_box!(d, construct(pair(concepts(db, ["tabular-file", "filename"])...)))
read_file = add_box!(d, concept(db, "read-tabular-file"))
data_x = add_box!(d, Box(concepts(db, ["table"]), concepts(db, ["table"])))
data_y = add_box!(d, Box(concepts(db, ["table"]), concepts(db, ["column"])))
ols = add_box!(d, construct(concept(db, "least-squares")))
fit = add_box!(d, concept(db, "fit-supervised"))
predict = add_box!(d, concept(db, "predict"))
error_l1 = add_box!(d, concept(db, "mean-absolute-error"))
error_l2 = add_box!(d, concept(db, "mean-squared-error"))
add_wires!(d, [
  (file, 1) => (read_file, 1),
  (read_file, 1) => (data_x, 1),
  (read_file, 1) => (data_y, 1),
  (ols, 1) => (fit, 1),
  (data_x, 1) => (fit, 2),
  (data_y, 1) => (fit, 3),
  (fit, 1) => (predict, 1),
  (data_x, 1) => (predict, 2),
  (data_y, 1) => (error_l1, 1),
  (predict, 1) => (error_l1, 2),
  (data_y, 1) => (error_l2, 1),
  (predict, 1) => (error_l2, 2),
])
@test semantic == d

# Linear regression on an R dataset using statsmodels.
semantic = read_semantic_graph(joinpath("python", "statsmodels_regression"))
d = WiringDiagram([], concepts(db, ["linear-regression"]))
r_data = add_box!(d, construct(pair(concept(db, "r-dataset-name"),
                                    concept(db, "r-dataset-package"))))
read_table = add_box!(d, concept(db, "read-table"))
ols = add_box!(d, construct(concept(db, "least-squares")))
eval_formula = add_box!(d, concept(db, "evaluate-formula-supervised"))
fit = add_box!(d, concept(db, "fit-supervised"))
add_wires!(d, [
  (r_data, 1) => (read_table, 1),
  (read_table, 1) => (eval_formula, 2),
  (ols, 1) => (fit, 1),
  (eval_formula, 1) => (fit, 2),
  (eval_formula, 2) => (fit, 3),
  (fit, 1) => (output_id(d), 1),
])
@test semantic == d

end

max_stars_repo_name: epatters/flowgraphs

max_issues_repo_name: epatters/flowgraphs

max_forks_repo_name: epatters/flowgraphs