text:
"""
You can easily define your own plotting recipes with convenience methods:

```
@userplot type GroupHist
    args
end

@recipe function f(gh::GroupHist)
    # set some attributes, add some series, using gh.args as input
end

# now you can plot like:
grouphist(rand(1000,4))
```
"""
macro userplot(expr)
    _userplot(expr)
end

function _userplot(expr::Expr)
    if expr.head != :type
        errror("Must call userplot on a type/immutable expression.  Got: $expr")
    end

    typename = expr.args[2]
    funcname = Symbol(lowercase(string(typename)))
    funcname2 = Symbol(funcname, "!")

    # return a code block with the type definition and convenience plotting methods
    esc(quote
        $expr
        export $funcname, $funcname2
        $funcname(args...; kw...) = plot($typename(args); kw...)
        $funcname2(args...; kw...) = plot!($typename(args); kw...)
    end)
end

function _userplot(sym::Symbol)
    _userplot(:(type $sym
            args
    end))
end

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

const _series_recipe_deps = Dict()

function series_recipe_dependencies(st::Symbol, deps::Symbol...)
    _series_recipe_deps[st] = deps
end

function seriestype_supported(st::Symbol)
    seriestype_supported(backend(), st)
end

# returns :no, :native, or :recipe depending on how it's supported
function seriestype_supported(pkg::AbstractBackend, st::Symbol)
    # is it natively supported
    if is_seriestype_supported(pkg, st)
        return :native
    end

    haskey(_series_recipe_deps, st) || return :no

    supported = true
    for dep in _series_recipe_deps[st]
        if seriestype_supported(pkg, dep) == :no
            supported = false
        end
    end
    supported ? :recipe : :no
end

macro deps(st, args...)
    :(Plots.series_recipe_dependencies($(quot(st)), $(map(quot, args)...)))
end

# get a list of all seriestypes
function all_seriestypes()
    sts = Set{Symbol}(keys(_series_recipe_deps))
    for bsym in backends()
        btype = _backendType[bsym]
        sts = union(sts, Set{Symbol}(supported_seriestypes(btype())))
    end
    sort(collect(sts))
end

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

num_series(x::AMat) = size(x,2)
num_series(x) = 1

RecipesBase.apply_recipe{T}(d::KW, ::Type{T}, plt::Plot) = throw(MethodError("Unmatched plot recipe: $T"))

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

# for seriestype `line`, need to sort by x values
@recipe function f(::Type{Val{:line}}, x, y, z)
    indices = sortperm(x)
    x := x[indices]
    y := y[indices]
    if typeof(z) <: AVec
        z := z[indices]
    end
    seriestype := :path
    ()
end
@deps line path

function hvline_limits(axis::Axis)
    vmin, vmax = axis_limits(axis)
    if vmin >= vmax
        if isfinite(vmin)
            vmax = vmin + 1
        else
            vmin, vmax = 0.0, 1.1
        end
    end
    vmin, vmax
end

@recipe function f(::Type{Val{:hline}}, x, y, z)
    xmin, xmax = hvline_limits(d[:subplot][:xaxis])
    n = length(y)
    newx = repmat(Float64[xmin, xmax, NaN], n)
    newy = vec(Float64[yi for i=1:3,yi=y])
    x := newx
    y := newy
    seriestype := :path
    ()
end
@deps hline path

@recipe function f(::Type{Val{:vline}}, x, y, z)
    ymin, ymax = hvline_limits(d[:subplot][:yaxis])
    n = length(y)
    newx = vec(Float64[yi for i=1:3,yi=y])
    newy = repmat(Float64[ymin, ymax, NaN], n)
    x := newx
    y := newy
    seriestype := :path
    ()
end
@deps vline path

# ---------------------------------------------------------------------------
# steps

function make_steps(x, y, st)
    n = length(x)
    n == 0 && return zeros(0),zeros(0)
    newx, newy = zeros(2n-1), zeros(2n-1)
    for i=1:n
        idx = 2i-1
        newx[idx] = x[i]
        newy[idx] = y[i]
        if i > 1
            newx[idx-1] = x[st == :steppre ? i-1 : i]
            newy[idx-1] = y[st == :steppre ? i   : i-1]
        end
    end
    newx, newy
end

# create a path from steps
@recipe function f(::Type{Val{:steppre}}, x, y, z)
    d[:x], d[:y] = make_steps(x, y, :steppre)
    seriestype := :path

    # create a secondary series for the markers
    if d[:markershape] != :none
        @series begin
            seriestype := :scatter
            x := x
            y := y
            label := ""
            primary := false
            ()
        end
        markershape := :none
    end
    ()
end
@deps steppre path scatter

# create a path from steps
@recipe function f(::Type{Val{:steppost}}, x, y, z)
    d[:x], d[:y] = make_steps(x, y, :steppost)
    seriestype := :path

    # create a secondary series for the markers
    if d[:markershape] != :none
        @series begin
            seriestype := :scatter
            x := x
            y := y
            label := ""
            primary := false
            ()
        end
        markershape := :none
    end
    ()
end
@deps steppost path scatter

# ---------------------------------------------------------------------------
# sticks

# create vertical line segments from fill
@recipe function f(::Type{Val{:sticks}}, x, y, z)
    n = length(x)
    fr = d[:fillrange]
    if fr == nothing
        yaxis = d[:subplot][:yaxis]
        fr = if yaxis[:scale] == :identity
            0.0
        else
            min(axis_limits(yaxis)[1], minimum(y))
        end
    end
    newx, newy = zeros(3n), zeros(3n)
    for i=1:n
        rng = 3i-2:3i
        newx[rng] = [x[i], x[i], NaN]
        newy[rng] = [cycle(fr,i), y[i], NaN]
    end
    x := newx
    y := newy
    fillrange := nothing
    seriestype := :path

    # create a secondary series for the markers
    if d[:markershape] != :none
        @series begin
            seriestype := :scatter
            x := x
            y := y
            label := ""
            primary := false
            ()
        end
        markershape := :none
    end
    ()
end
@deps sticks path scatter

# ---------------------------------------------------------------------------
# bezier curves

# get the value of the curve point at position t
function bezier_value(pts::AVec, t::Real)
    val = 0.0
    n = length(pts)-1
    for (i,p) in enumerate(pts)
        val += p * binomial(n, i-1) * (1-t)^(n-i+1) * t^(i-1)
    end
    val
end

# create segmented bezier curves in place of line segments
@recipe function f(::Type{Val{:curves}}, x, y, z; npoints = 30)
    args = z != nothing ? (x,y,z) : (x,y)
    newx, newy = zeros(0), zeros(0)
    fr = d[:fillrange]
    newfr = fr != nothing ? zeros(0) : nothing
    newz = z != nothing ? zeros(0) : nothing
    # lz = d[:line_z]
    # newlz = lz != nothing ? zeros(0) : nothing

    # for each line segment (point series with no NaNs), convert it into a bezier curve
    # where the points are the control points of the curve
    for rng in iter_segments(args...)
        length(rng) < 2 && continue
        ts = linspace(0, 1, npoints)
        nanappend!(newx, map(t -> bezier_value(cycle(x,rng), t), ts))
        nanappend!(newy, map(t -> bezier_value(cycle(y,rng), t), ts))
        if z != nothing
            nanappend!(newz, map(t -> bezier_value(cycle(z,rng), t), ts))
        end
        if fr != nothing
            nanappend!(newfr, map(t -> bezier_value(cycle(fr,rng), t), ts))
        end
        # if lz != nothing
        #     lzrng = cycle(lz, rng) # the line_z's for this segment
        #     push!(newlz, 0.0)
        #     append!(newlz, map(t -> lzrng[1+floor(Int, t * (length(rng)-1))], ts))
        # end
    end

    x := newx
    y := newy
    if z == nothing
        seriestype := :path
    else
        seriestype := :path3d
        z := newz
    end
    if fr != nothing
        fillrange := newfr
    end
    # if lz != nothing
    #     # line_z := newlz
    #     linecolor := (isa(d[:linecolor], ColorGradient) ? d[:linecolor] : cgrad())
    # end
    # Plots.DD(d)
    ()
end
@deps curves path

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

# create a bar plot as a filled step function
@recipe function f(::Type{Val{:bar}}, x, y, z)
    nx, ny = length(x), length(y)
    axis = d[:subplot][isvertical(d) ? :xaxis : :yaxis]
    cv = [discrete_value!(axis, xi)[1] for xi=x]
    x = if nx == ny
        cv
    elseif nx == ny + 1
        0.5diff(cv) + cv[1:end-1]
    else
        error("bar recipe: x must be same length as y (centers), or one more than y (edges).\n\t\tlength(x)=$(length(x)), length(y)=$(length(y))")
    end

    # compute half-width of bars
    bw = d[:bar_width]
    hw = if bw == nothing
        0.5mean(diff(x))
    else
        Float64[0.5cycle(bw,i) for i=1:length(x)]
    end

    # make fillto a vector... default fills to 0
    fillto = d[:fillrange]
    if fillto == nothing
        fillto = 0
    end

    # create the bar shapes by adding x/y segments
    xseg, yseg = Segments(), Segments()
    for i=1:ny
        center = x[i]
        hwi = cycle(hw,i)
        yi = y[i]
        fi = cycle(fillto,i)
        push!(xseg, center-hwi, center-hwi, center+hwi, center+hwi, center-hwi)
        push!(yseg, yi, fi, fi, yi, yi)
    end

    # widen limits out a bit
    expand_extrema!(axis, widen(extrema(xseg.pts)...))

    # switch back
    if !isvertical(d)
        xseg, yseg = yseg, xseg
    end

    # reset orientation
    orientation := default(:orientation)

    x := xseg.pts
    y := yseg.pts
    seriestype := :shape
    ()
end
@deps bar shape

# ---------------------------------------------------------------------------
# Histograms

# edges from number of bins
function calc_edges(v, bins::Integer)
    vmin, vmax = extrema(v)
    linspace(vmin, vmax, bins+1)
end

# just pass through arrays
calc_edges(v, bins::AVec) = bins

# find the bucket index of this value
function bucket_index(vi, edges)
    for (i,e) in enumerate(edges)
        if vi <= e
            return max(1,i-1)
        end
    end
    return length(edges)-1
end

function my_hist(v, bins; normed = false, weights = nothing)
    edges = calc_edges(v, bins)
    counts = zeros(length(edges)-1)

    # add a weighted count
    for (i,vi) in enumerate(v)
        idx = bucket_index(vi, edges)
        counts[idx] += (weights == nothing ? 1.0 : weights[i])
    end

    # normalize by bar area?
    norm_denom = normed ? sum(diff(edges) .* counts) : 1.0
    if norm_denom == 0
        norm_denom = 1.0
    end

    edges, counts ./ norm_denom
end

@recipe function f(::Type{Val{:histogram}}, x, y, z)
    edges, counts = my_hist(y, d[:bins],
                               normed = d[:normalize],
                               weights = d[:weights])
    x := edges
    y := counts
    seriestype := :bar
    ()
end
@deps histogram bar

# ---------------------------------------------------------------------------
# Histogram 2D

# if tuple, map out bins, otherwise use the same for both
calc_edges_2d(x, y, bins) = calc_edges(x, bins), calc_edges(y, bins)
calc_edges_2d{X,Y}(x, y, bins::Tuple{X,Y}) = calc_edges(x, bins[1]), calc_edges(y, bins[2])

# the 2D version
function my_hist_2d(x, y, bins; normed = false, weights = nothing)
    xedges, yedges = calc_edges_2d(x, y, bins)
    counts = zeros(length(yedges)-1, length(xedges)-1)

    # add a weighted count
    for i=1:length(x)
        r = bucket_index(y[i], yedges)
        c = bucket_index(x[i], xedges)
        counts[r,c] += (weights == nothing ? 1.0 : weights[i])
    end

    # normalize to cubic area of the imaginary surface towers
    norm_denom = normed ? sum((diff(yedges) * diff(xedges)') .* counts) : 1.0
    if norm_denom == 0
        norm_denom = 1.0
    end

    xedges, yedges, counts ./ norm_denom
end

centers(v::AVec) = 0.5 * (v[1:end-1] + v[2:end])

@recipe function f(::Type{Val{:histogram2d}}, x, y, z)
    xedges, yedges, counts = my_hist_2d(x, y, d[:bins],
                                              normed = d[:normalize],
                                              weights = d[:weights])
    for (i,c) in enumerate(counts)
        if c == 0
            counts[i] = NaN
        end
    end
    x := centers(xedges)
    y := centers(yedges)
    z := Surface(counts)
    linewidth := 0
    seriestype := :heatmap
    ()
end
@deps histogram2d heatmap

# ---------------------------------------------------------------------------
# scatter 3d

@recipe function f(::Type{Val{:scatter3d}}, x, y, z)
    seriestype := :path3d
    if d[:markershape] == :none
        markershape := :circle
    end
    linewidth := 0
    linealpha := 0
    ()
end

# note: don't add dependencies because this really isn't a drop-in replacement

# ---------------------------------------------------------------------------
# contourf - filled contours

@recipe function f(::Type{Val{:contourf}}, x, y, z)
    fillrange := true
    seriestype := :contour
    ()
end

# ---------------------------------------------------------------------------
# Error Bars

function error_style!(d::KW)
    d[:seriestype] = :path
    d[:linecolor] = d[:markerstrokecolor]
    d[:linewidth] = d[:markerstrokewidth]
    d[:label] = ""
end

# if we're passed a tuple of vectors, convert to a vector of tuples
function error_zipit(ebar)
    if istuple(ebar)
        collect(zip(ebar...))
    else
        ebar
    end
end

function error_coords(xorig, yorig, ebar)
    # init empty x/y, and zip errors if passed Tuple{Vector,Vector}
    x, y = zeros(0), zeros(0)

    # for each point, create a line segment from the bottom to the top of the errorbar
    for i = 1:max(length(xorig), length(yorig))
        xi = cycle(xorig, i)
        yi = cycle(yorig, i)
        ebi = cycle(ebar, i)
        nanappend!(x, [xi, xi])
        e1, e2 = if istuple(ebi)
            first(ebi), last(ebi)
        elseif isscalar(ebi)
            ebi, ebi
        else
            error("unexpected ebi type $(typeof(ebi)) for errorbar: $ebi")
        end
        nanappend!(y, [yi - e1, yi + e2])
    end
    x, y
end

# we will create a series of path segments, where each point represents one
# side of an errorbar
@recipe function f(::Type{Val{:yerror}}, x, y, z)
    error_style!(d)
    markershape := :hline
    d[:x], d[:y] = error_coords(d[:x], d[:y], error_zipit(d[:yerror]))
    ()
end
@deps yerror path

@recipe function f(::Type{Val{:xerror}}, x, y, z)
    error_style!(d)
    markershape := :vline
    d[:y], d[:x] = error_coords(d[:y], d[:x], error_zipit(d[:xerror]))
    ()
end
@deps xerror path

# TODO: move quiver to PlotRecipes

# ---------------------------------------------------------------------------
# quiver

# function apply_series_recipe(d::KW, ::Type{Val{:quiver}})
function quiver_using_arrows(d::KW)
    d[:label] = ""
    d[:seriestype] = :path
    if !isa(d[:arrow], Arrow)
        d[:arrow] = arrow()
    end

    velocity = error_zipit(d[:quiver])
    xorig, yorig = d[:x], d[:y]

    # for each point, we create an arrow of velocity vi, translated to the x/y coordinates
    x, y = zeros(0), zeros(0)
    for i = 1:max(length(xorig), length(yorig))
        # get the starting position
        xi = cycle(xorig, i)
        yi = cycle(yorig, i)

        # get the velocity
        vi = cycle(velocity, i)
        vx, vy = if istuple(vi)
            first(vi), last(vi)
        elseif isscalar(vi)
            vi, vi
        elseif isa(vi,Function)
            vi(xi, yi)
        else
            error("unexpected vi type $(typeof(vi)) for quiver: $vi")
        end

        # add the points
        nanappend!(x, [xi, xi+vx, NaN])
        nanappend!(y, [yi, yi+vy, NaN])
    end

    d[:x], d[:y] = x, y
    # KW[d]
end

# function apply_series_recipe(d::KW, ::Type{Val{:quiver}})
function quiver_using_hack(d::KW)
    d[:label] = ""
    d[:seriestype] = :shape

    velocity = error_zipit(d[:quiver])
    xorig, yorig = d[:x], d[:y]

    # for each point, we create an arrow of velocity vi, translated to the x/y coordinates
    pts = P2[]
    for i = 1:max(length(xorig), length(yorig))

        # get the starting position
        xi = cycle(xorig, i)
        yi = cycle(yorig, i)
        p = P2(xi, yi)

        # get the velocity
        vi = cycle(velocity, i)
        vx, vy = if istuple(vi)
            first(vi), last(vi)
        elseif isscalar(vi)
            vi, vi
        elseif isa(vi,Function)
            vi(xi, yi)
        else
            error("unexpected vi type $(typeof(vi)) for quiver: $vi")
        end
        v = P2(vx, vy)

        dist = norm(v)
        arrow_h = 0.1dist          # height of arrowhead
        arrow_w = 0.5arrow_h       # halfwidth of arrowhead
        U1 = v ./ dist             # vector of arrowhead height
        U2 = P2(-U1[2], U1[1])     # vector of arrowhead halfwidth
        U1 *= arrow_h
        U2 *= arrow_w

        ppv = p+v
        nanappend!(pts, P2[p, ppv-U1, ppv-U1+U2, ppv, ppv-U1-U2, ppv-U1])
    end

    d[:x], d[:y] = Plots.unzip(pts[2:end])
    # KW[d]
end

# function apply_series_recipe(d::KW, ::Type{Val{:quiver}})
@recipe function f(::Type{Val{:quiver}}, x, y, z)
    if :arrow in supported_attrs()
        quiver_using_arrows(d)
    else
        quiver_using_hack(d)
    end
    ()
end
@deps quiver shape path

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

# TODO: move OHLC to PlotRecipes finance.jl

type OHLC{T<:Real}
  open::T
  high::T
  low::T
  close::T
end
Base.convert(::Type{OHLC}, tup::Tuple) = OHLC(tup...)
# Base.tuple(ohlc::OHLC) = (ohlc.open, ohlc.high, ohlc.low, ohlc.close)

# get one OHLC path
function get_xy(o::OHLC, x, xdiff)
    xl, xm, xr = x-xdiff, x, x+xdiff
    ox = [xl, xm, NaN,
          xm, xm, NaN,
          xm, xr]
    oy = [o.open, o.open, NaN,
          o.low, o.high, NaN,
          o.close, o.close]
    ox, oy
end

# get the joined vector
function get_xy(v::AVec{OHLC}, x = 1:length(v))
    xdiff = 0.3mean(abs(diff(x)))
    x_out, y_out = zeros(0), zeros(0)
    for (i,ohlc) in enumerate(v)
        ox,oy = get_xy(ohlc, x[i], xdiff)
        nanappend!(x_out, ox)
        nanappend!(y_out, oy)
    end
    x_out, y_out
end

# these are for passing in a vector of OHLC objects
# TODO: when I allow `@recipe f(::Type{T}, v::T) = ...` definitions to replace convertToAnyVector,
#       then I should replace these with one definition to convert to a vector of 4-tuples

# to squash ambiguity warnings...
@recipe f(x::AVec{Function}, v::AVec{OHLC}) = error()
@recipe f{R1<:Number,R2<:Number,R3<:Number,R4<:Number}(x::AVec{Function}, v::AVec{Tuple{R1,R2,R3,R4}}) = error()

# this must be OHLC?
@recipe f{R1<:Number,R2<:Number,R3<:Number,R4<:Number}(x::AVec, ohlc::AVec{Tuple{R1,R2,R3,R4}}) = x, OHLC[OHLC(t...) for t in ohlc]

@recipe function f(x::AVec, v::AVec{OHLC})
    seriestype := :path
    get_xy(v, x)
end

@recipe function f(v::AVec{OHLC})
    seriestype := :path
    get_xy(v)
end

# the series recipe, when passed vectors of 4-tuples

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

# TODO: everything below here should be either changed to a
#       series recipe or moved to PlotRecipes

# "Sparsity plot... heatmap of non-zero values of a matrix"
# function spy{T<:Real}(z::AMat{T}; kw...)
#     mat = map(zi->float(zi!=0), z)'
#     xn, yn = size(mat)
#     heatmap(mat; leg=false, yflip=true, aspect_ratio=:equal,
#         xlim=(0.5, xn+0.5), ylim=(0.5, yn+0.5),
#         kw...)
# end

# Only allow matrices through, and make it seriestype :spy so the backend can
# optionally handle it natively.

@userplot Spy

@recipe function f(g::Spy)
    @assert length(g.args) == 1 && typeof(g.args[1]) <: AbstractMatrix
    seriestype := :spy
    mat = g.args[1]
    n,m = size(mat)
    Plots.SliceIt, 1:m, 1:n, Surface(mat)
end

@recipe function f(::Type{Val{:spy}}, x,y,z)
    yflip := true
    aspect_ratio := 1
    rs, cs, zs = findnz(z.surf)
    xlim := extrema(cs)
    ylim := extrema(rs)
    if d[:markershape] == :none
        markershape := :circle
    end
    if d[:markersize] == default(:markersize)
        markersize := 1
    end
    markerstrokewidth := 0
    marker_z := zs
    label := ""
    x := cs
    y := rs
    z := nothing
    seriestype := :scatter
    ()
end

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

"Adds a+bx... straight line over the current plot"
function abline!(plt::Plot, a, b; kw...)
    plot!(plt, [extrema(plt)...], x -> b + a*x; kw...)
end

abline!(args...; kw...) = abline!(current(), args...; kw...)

# -------------------------------------------------
# Dates

@recipe f(::Type{Date}, dt::Date) = (dt -> convert(Int,dt), dt -> string(convert(Date,dt)))
@recipe f(::Type{DateTime}, dt::DateTime) = (dt -> convert(Int,dt), dt -> string(convert(DateTime,dt)))


text:
"""
    struct Relationship

A structure designed to capture the relationship between two individuals
"""
struct Relationship
    common_ancestor::Integer # the vertex id of the common ancestor
    individual_a::Integer # the vertex id of the first individual
    individual_b::Integer # the vertex id of the second individual
    a_to_ancestor::AInteger # the relationship of the first individual to the common ancestor
    b_to_ancestor::AInteger # the relationship of the second individual to the common ancestor
end

function cousinth(generational)
    if generational == 2
        return "1st cousin"
    elseif generational == 3
        return "2nd cousin"
    elseif generational == 4
        return "3rd cousin"
    else
        return "$(generational - 1)th cousin"
    end
end

function grandparentith(removal)
    if removal == 4
        return "2nd"
    elseif removal == 5
        return "3rd"
    else
        return "$(removal - 2)th"
    end
end

"""
    lookup_relationship(AInteger, AInteger)

Given two individuals with a common ancestor where the first parameter is the generational distance
from the first individual to the common ancestor and the second parameter is the generational
distance from the second individual to the common ancestor, this function provides a string
representation of how the two individuals are related.

"""
function lookup_relationship(da::Integer, db::Integer)
    if da > db
        removal = da - db
    else
        removal = db - da
    end
    if da == 0 && db == 0
        return "self"
    elseif da == 0 || db == 0
        # direct ancestor
        if removal == 1
            return "parent"
        elseif removal == 2
            return "grandparent"
        elseif removal == 3
            return "great grandparent"
        else
            return "$(grandparentith(removal)) great grandparent"
        end
    else
        generational = min(da, db)
        if generational == 1
            if removal == 0
                return "sibling"
            elseif removal == 1
                return "aunt/uncle"
            elseif removal == 2
                return "great aunt/uncle"
            else
                return "$(grandparentith(removal + 1)) great aunt/uncle"
            end
        elseif removal == 0
            return cousinth(generational)
        else
            return "$(cousinth(generational)) $removal removed"
        end
    end
end

function lookup_relationship(r::Relationship)
    a = r.a_to_ancestor
    b = r.b_to_ancestor
    if a == red_one(a)
        a = one(a)
    end
    if b == red_one(b)
        b = one(b)
    end
    if a == 0 || b == 0
        return "No Relationship"
    end
    return lookup_relationship(MSB(a), MSB(b))
end

function calculate_relationship(a::Vector{T}, b::Vector{T}) where T <: AInteger
    common_ancestor = 0
    closest_relationship = (T(0), T(0))
    individual_a = 0
    individual_b = 0
    for (index, value) in enumerate(zip(a, b))
        if value[1] == 1 || value[1] == red_one(value[1])
            individual_a = index
        end
        if value[2] == 1 || value[2] == red_one(value[1])
            individual_b = index
        end
        if value[1] != 0 && value[2] != 0
            if value[1] + value[2] < closest_relationship[1] + closest_relationship[2]
                common_ancestor = index
                closest_relationship = value
            end
        end
    end
    return Relationship(common_ancestor, individual_a, individual_b, closest_relationship[1], closest_relationship[2])
end

ext: jl


text:
module RobustShortestPath

# package code goes here

using LightGraphs

include("misc.jl")
include("one.jl")
include("two.jl")

export
	get_shortest_path,
	get_robust_path,
	get_robust_path_two

end # module


text:
using ThArrays
using Tracker: forward, data
using Random
using Test

using ThArrays.TrackerAD: _th, _tr

@testset "Torch and Tracker Mixed Gradient" begin

    Random.seed!(0);

    a = rand(3, 2)
    b = rand(3, 2)

    @testset "Simple Mixed AD" begin
        # all with th(PyTorch Backend)
        f1(x, y) = sum(sin.(_th(x)) + sin.(_th(y)))
        # sin with tr, (+, sum) wtih th
        f2(x, y) = sum(_th(sin.(x)) + _th(sin.(y)))
        # all with tr(Tracker Backend)
        f3(x, y) = sum((sin.(x)) + (sin.(y)))
        # (sin, +) with th, sum with tr
        f4(x, y) = sum(_tr(sin.(_th(x)) + sin.(_th(y))))
        # (sin, +) with tr, sum with th
        f5(x, y) = sum(_th(sin.(x) + sin.(y)))

        y1, back1 = forward(f1, a, b)
        y2, back2 = forward(f2, a, b)
        y3, back3 = forward(f3, a, b)
        y4, back4 = forward(f4, a, b)
        y5, back5 = forward(f5, a, b)

        b1 = data(back1(2))
        b2 = data(back2(2))
        b3 = data(back3(2))
        b4 = data(back4(2))
        b5 = data(back5(2))

        # @show y1, y2, y3, y4, y5
        # @show b1, b2, b3, b4, b5
        @test b1 == b2 == b3 == b4 == b5
    end

end


text:
using ParticleInCell
using LinearAlgebra

@testset "FDTD solver" begin

    dimx, dimy = 2π, 2π
    nx, ny = 128, 128
    dt = 1e-4
    nstep = 8

    mesh = TwoDGrid(dimx, nx, dimy, ny)

    ex = zeros(nx+1, ny+1)
    ey = zeros(nx+1, ny+1)
    bz = zeros(nx+1, ny+1)
    jx = zeros(nx+1, ny+1)
    jy = zeros(nx+1, ny+1)

    fdtd = FDTD(mesh)

    ω = sqrt(2)

    xn = mesh.x
    yn = mesh.y |> transpose
    xc = (mesh.x[1:end-1] .+ mesh.x[2:end]) ./ 2
    yc = (mesh.y[1:end-1] .+ mesh.y[2:end]) ./ 2 |> transpose

    # Ex and Ey are set to zero at t = 0.0
    # Bz is set at  t = dt/2
    t = 0.5dt
    fdtd.bz .= -cos.(xc) .* cos.(yc) .* cos(ω * t)

    sol_ex(x, y) = +cos.(x) .* sin.(y) ./ ω
    sol_ey(x, y) = -sin.(x) .* cos.(y) ./ ω
    sol_bz(x, y) = -cos.(x) .* cos.(y)

    for istep = 1:nstep # Loop over time

        ampere_maxwell!(fdtd, mesh, jx, jy, dt)

        t = t + 0.5dt

        @test maximum(abs.(fdtd.ex .- sol_ex(xc, yn) .* sin(ω * t))) < 1e-6
        @test maximum(abs.(fdtd.ey .- sol_ey(xn, yc) .* sin(ω * t))) < 1e-6

        update_fields!(ex, ey, bz, mesh, fdtd)

        @test maximum(abs.(ex .- sol_ex(xn, yn) .* sin(ω * t))) < 1e-6
        @test maximum(abs.(ey .- sol_ey(xn, yn) .* sin(ω * t))) < 1e-6

        faraday!(fdtd, mesh, dt)

        t = t + 0.5dt

        @test maximum(abs.(fdtd.bz .- sol_bz(xc, yc) .* cos(ω * t))) < 1e-6

        update_fields!(ex, ey, bz, mesh, fdtd)

        @test maximum(abs.(bz .- sol_bz(xn, yn) .* cos(ω * t))) < 1e-3

    end # next time step

end


text: abstract type Extrapolation end


text:
"""
Handling convergence ciretron for iterative algorithms

@author : Spencer Lyon <spencer.lyon@stern.nyu.edu>
@date : 2015-04-13 11:29:12

"""
module IterationManagers

abstract IterationManager
abstract IterationState{T}

export
# types
    IterationManager, IterationState, IterTolManager,  DefaultManager,
    IterManager, DefaultState, ExtraState,

# functions
    finished, update!, managed_iteration,
    verbose, print_now, display_iter, num_iter

for s in ["default", "extra"]
    include("states/$s.jl")
end

for m in ["itertol", "tol", "iter"]
    include("managers/$m.jl")
end

include("api.jl")

end # module


text:
#Make comparison plot with theoretical model
function make_energy_plot(  p::Param ;
                            theory_comparison = :true,
                            gamma_comparison  = :false,
                            gamma_list = nothing,
                            gamma_legend = nothing,
                            mscale = log10 ,
                            oname = "default_name",
                            xlims = nothing,
                            ylims = nothing ,
                            energy_data_every = 10)

    #Welcome information
    println("\n \n")
    println("___________________________________________________________")
    println(" ")
    println(" Making Energy Plot from file: /data" , p.folder  ,p.fname)
    println(" ")

    #Default values
    my_xlims = xlims
    my_ylims = ylims

    if( xlims === nothing )
        my_xlims = [0 , p.tmax]
    end
    if( ylims === nothing )
        my_ylims = [1E-4 , 5.0]
    end

    #Activatecairo
    set_my_cairo_theme()

    #Create figure
    fig = Figure(resolution = (650, 400))
    ax1 = fig[1, 1] = Axis(fig , yscale = mscale)

    if( theory_comparison)

        #Theoretical model
        x , y = get_simple_theoretical_model_data(p0 , my_xlims[2])
        lines!(ax1 , x , y , linewidth = 1 ,  linestyle = :dash , color = "red" , label = "Model")

    end

    if(gamma_comparison)

        for (i,g) in enumerate(gamma_list)

            mcolor = [ "green" , "blue", "red"]

            #Theoretical model
            x = 0:0.1:my_xlims[2]
            y = exp.((-2*g).*x)
            lines!(ax1 , x ,  y, linewidth = 1 ,  linestyle = :dash , color = mcolor[i] , label = gamma_legend[i] )#("γ = " * "$(g).") )

        end
    end

    #Actual data
    x , y = get_energy(p , energy_data_every , p.tmax)
    lines!(ax1 , x , y/y[1] , linewidth = 1 ,  linestyle = :solid , color = "black", label = "Simulation")

    #Configure plot
    ax1.xlabel = "Time"
    ax1.ylabel = "Energy Density [ ϵ(t) / ϵ(0) ]"

    limits!(ax1 , my_xlims[1]  , my_xlims[2] ,   my_ylims[1]  ,   my_ylims[2]   )

    #Legend
    axislegend(ax1 , position = :lb)

    fdata   = "./data"
    outname   = fdata * p.folder *  "$(oname).png"

    #outname = abspath(p.folder, p.fname * ".png")

    save(outname, fig, px_per_unit = 10)

    println(" Image saved to: " , outname )

    println("___________________________________________________________")

end

#Make plot for several parameters
function make_energy_plot_params(   s  , legends ;
                                    gamma_comparison  = :false,
                                    gamma_list = nothing,
                                    gamma_legend = nothing,
                                    mscale = log10 ,
                                    oname = "default_name",
                                    xlims = nothing,
                                    ylims = nothing ,
                                    energy_data_every = 10,
                                    plot_title = "")

    #Welcome information
    println("\n \n")
    println("___________________________________________________________")
    println(" ")
    println(" Making Energy Plot from files: ")

    for p in s
        println("/data" , p.folder  ,p.fname)
    end

    println(" ")

    #Default values
    my_xlims = xlims
    my_ylims = ylims

    if( xlims === nothing )
        my_xlims = [0 , s[1].tmax]
    end

    if( ylims === nothing )
        my_ylims = [1E-5 , 5.0]
    end

    #Activatecairo
    set_my_cairo_theme()

    #Create figure
    fig = Figure(resolution = (650, 400) )
    ax1 = fig[1, 1] = Axis(fig , yscale = mscale , title = plot_title )

    pallete = ["#111111", "#65ADC2", "#233B43", "#E84646", "#C29365", "#362C21", "#316675", "#168E7F", "#109B37"]
    #colormap( "Oranges" , 10)

    for (i , p) in enumerate(s)

        x , y = get_energy(p , energy_data_every , p.tmax)

        lines!(ax1 , x , y./y[1] , color = pallete[4+i] , linewidth = 1, label = legends[i])

    end

    #Gamma comparison
    if(gamma_comparison)

        for (i,g) in enumerate(gamma_list)

            #Theoretical model
            x = 0:0.1:3000
            y = exp.((-2*g).*x)
            lines!(ax1 , x ,  y, linewidth = 1 ,  linestyle = :dash , color = pallete[4+i]  )#("γ = " * "$(g).") )

        end
    end

    #Configure plot
    ax1.xlabel = "Time"
    ax1.ylabel = "Energy Density [ ϵ(t) / ϵ(0) ]"

    limits!(ax1 , my_xlims[1]  , my_xlims[2] ,   my_ylims[1]  ,   my_ylims[2]   )

    #Legend
    axislegend(ax1 ,position = :rt)

    fdata   = "./data/"
    outname   = fdata *  "$(oname).png"
    save(outname, fig, px_per_unit = 10)

    println(" Image saved to: " , outname )

    println("___________________________________________________________")

end


text:
# using ArviZ
using CSV
using DataFrames
using Distributions
using NPZ
using StatsPlots
using Random
using Transits
using NestedSamplers
using StatsBase

rootdir(args...) = joinpath(@__DIR__, args...)
datadir(args...) = rootdir("data", args...)

data = npzread(datadir("generated_data_100ppm.npz"))
const t = data["t"]
const flux = data["flux"]
data = npzread(datadir("generated_data.npz"))
const true_mod = data["flux"]

ground_truth = Dict(
    "period" => 3.5,
    "t0" => 1.3,
    "us" => [0.5, 0.2],
    "r" => 0.03,
    "a" => 10.0,
    "yerr" => 1e-4
)

phases = @. (t - ground_truth["t0"] + ground_truth["period"] / 2) % ground_truth["period"] - ground_truth["period"] / 2
maxphase = 0.1
inds = @. abs(phases) < maxphase

# show data
scatter(phases[inds], flux[inds], yerr=ground_truth["yerr"], c=:black, label="data")
plot!(phases[inds], true_mod[inds], label="truth")

orbit = KeplerianOrbit(ground_truth["a"], ground_truth["period"], 0, ground_truth["t0"], 0)
law = PolynomialLimbDark(ground_truth["us"])
f = @. law(orbit, t, ground_truth["r"])
plot(phases[inds], f[inds], label="transits")

# Define a function to optimize.
function loglike(X)
    logP, t0, r, loga, yerr = @view X[1:5]
    us = @view X[end-1:end]

    orbit = KeplerianOrbit(exp(loga), exp(logP), 0, t0, 0)
    ld = QuadLimbDark(us)
    ll = zero(eltype(X))
    for idx in eachindex(flux)
        dist = Normal(ld(orbit, t[idx], r), yerr)
        ll += logpdf(dist, flux[idx])
    end
    return ll
end

function lightcurve(X)
    logP, t0, r, loga, yerr = @view X[1:5]
    us = @view X[end-1:end]

    orbit = KeplerianOrbit(exp(loga), exp(logP), 0, t0, 0)
    ld = QuadLimbDark(us)
    return @. ld(orbit, t, r)
end

X0 = [log(3.5), 1.3, 0.03, log(10), 1e-4, 0.5, 0.2]

## Distributions
priors = Dict(
    :logP => Normal(log(ground_truth["period"]), 0.1),
    :t0 => Normal(ground_truth["t0"], 10),
    :r => Uniform(0.01, 0.1),
    :loga => Normal(log(ground_truth["a"]), 0.1),
    :yerr => Truncated(Cauchy(0, 0.1), 0, Inf),
)

function kiptform!(X)
    sqrtq1 = sqrt(first(X))
    twoq2 = 2 * last(X)
    X[begin] = sqrtq1 * twoq2
    X[end] = sqrtq1 * (1 - twoq2)
    return X
end

kiptform(X) = kiptform!(copy(X))

function prior_transform(X)
    out = copy(X)
    for (i, k) in enumerate((:logP, :t0, :r, :loga, :yerr))
        out[i] = quantile(priors[k], X[i])
    end
    us = @view out[end-1:end]
    kiptform!(us)
    return out
end

## sampling
model = NestedModel(loglike, prior_transform)
splr = Nested(7, 1000; proposal=Proposals.RWalk(ratio=0.9))

chain, state = sample(model, splr, param_names=["logP", "t0", "r", "loga", "yerr", "us[1]", "us[2]"])

max_stars_repo_name: mileslucas/yeehaw

max_issues_repo_name: mileslucas/yeehaw

max_forks_repo_name: mileslucas/yeehaw


text:
const basemap_Float32 = [0.0f0, 1.0f0, 2.0f0, NaN32]
const basemap_Float64 = [0.0, 1.0, 2.0, NaN]
const dosage_unit_Float64 = 1/(1 << 14)
const dosage_unit_Float32 = Float32(dosage_unit_Float64)

"""
    alt_allele_dosage!(buf, genobuf, p, v; genoldbuf)

Computes unphased biallelic dosage of ALT allele.

- `buf`: stores dosage values.
- `genobuf`: stores genotype values.
- `p`: a `Pgen` object.
- `v`: a `Variant` object.
- `genoldbuf`: most recent non-LD-compressed genotypes.

Returns:
- `buf`
- `genobuf`
- `offset`: end of dosage record on the current variant record track.
"""
function alt_allele_dosage!(buf::AbstractVector{T}, genobuf::AbstractVector{UInt8},
    p::Pgen, v::Variant;
    genoldbuf::Union{Nothing, Vector{UInt8}}=nothing
    ) where T <: AbstractFloat
    n_samples = p.header.n_samples
    @assert length(buf) == n_samples && length(genobuf) == n_samples

    _, variant_record, offset = get_genotypes!(genobuf, p, v; ldbuf=genoldbuf)

    # convert genotype to floating-point values
    basemap = T == Float32 ? basemap_Float32 : basemap_Float64
    @inbounds for i in 1:n_samples
        buf[i] = basemap[genobuf[i] + 1]
    end

    @assert v.record_type & 0x08 == 0 "Multiallelic case unsupported"
    if v.record_type & 0x08 != 0
        # TODO: compute offset for multiallelic genotype track
    end

    # TODO: test this.
    if v.record_type & 0x10 != 0 # skip phase info (aux track 2)
        offset = _phase_offset(p, genobuf, variant_record, offset)
    end

    if v.record_type & 0x60 == 0
        return buf, variant_record, offset
    elseif v.record_type & 0x20 !== 0 && v.record_type & 0x40 == 0
        # track 3 is a difflist
        dl, offset = parse_difflist(variant_record, offset, p.header.bytes_per_sample_id, false)
        offset = _get_difflist_dosage!(buf, p, dl, variant_record, offset)
    elseif v.record_type & 0x20 == 0 && v.record_type & 0x40 != 0
        # track 3 does not exist, 0xffff is missing value in track 4
        dosages = reinterpret(UInt16, @view(variant_record[offset + 1 : offset + 2 * n_samples]))
        dosage_unit = T == Float32 ? dosage_unit_Float32 : dosage_unit_Float64
        for i in 1:n_samples
            d = dosages[i]
            buf[i] = d == 0xffff ? NaN : d * dosage_unit
        end
        offset += 2 * n_samples
    else # track 3 is a bitarray for dosage existence in track 4.
        bv_nbytes = (n_samples + 7) >> 3
        bv = BitsVector(pointer(variant_record, offset + 1), 1, n_samples)
        #BitsVector(@view(variant_record[offset + 1 : offset + bv_nbytes]), 1, n_samples)
        offset += bv_nbytes
        dosage_len = sum(bv)
        dosages = reinterpret(UInt16, @view(variant_record[offset + 1 : offset + 2 * dosage_len]))
        dosage_unit = T == Float32 ? dosage_unit_Float32 : dosage_unit_Float64
        j = 1
        for i in 1:n_samples
            if bv[i] == 0x01
                buf[i] = dosages[j] * dosage_unit
                j += 1
            end
            if j > dosage_len
                break
            end
        end
        offset += 2 * dosage_len
    end
    return buf, genobuf, offset
end

"""
    ref_allele_dosage!(buf, genobuf, p, v; genoldbuf)

Computes unphased biallelic dosage of REF allele.

- `buf`: stores dosage values.
- `genobuf`: stores genotype values.
- `p`: a `Pgen` object.
- `v`: a `Variant` object.
- `genoldbuf`: most recent non-LD-compressed genotypes.

Returns:
- `buf`
- `genobuf`
- `offset`: end of dosage record on the current variant record track.
"""
function ref_allele_dosage!(buf::AbstractVector{T}, genobuf::AbstractVector{UInt8},
    p::Pgen, v::Variant;
    genoldbuf::Union{Nothing, Vector{UInt8}}=nothing
    ) where T <: AbstractFloat
    _, offset = alt_allele_dosage!(buf, genobuf, p, v; genoldbuf=genoldbuf)
    @inbounds for i in 1:p.header.n_samples
        buf[i] = 2.0 - buf[i]
    end
    return buf, genobuf, offset
end

"""
    alt_allele_dosage(p, v)

Computes unphased biallelic dosage of ALT allele.

- `p`: a `Pgen` object.
- `v`: a `Variant` object.

Returns:
- `buf`: stores dosage values.
- `genobuf`: stores genotype values.
- `offset`: end of dosage record on the current variant record track.
"""
function alt_allele_dosage(p::Pgen, v::Variant)
    n_samples = p.header.n_samples
    buf = Vector{Float32}(undef, n_samples)
    genobuf = Vector{UInt8}(undef, n_samples)
    alt_allele_dosage!(buf, genobuf, p, v)
end

"""
    ref_allele_dosage(p, v)

Computes unphased biallelic dosage of REF allele.

- `p`: a `Pgen` object.
- `v`: a `Variant` object.

Returns:
- `buf`: stores dosage values.
- `genobuf`: stores genotype values.
- `offset`: end of dosage record on the current variant record track.
"""
function ref_allele_dosage(p::Pgen, v::Variant)
    n_samples = p.header.n_samples
    buf = Vector{Float32}(undef, n_samples)
    genobuf = Vector{UInt8}(undef, n_samples)
    ref_allele_dosage!(buf, genobuf, p, v)
end

# This function is untested.
function _phase_offset(p::Pgen, g::AbstractVector{UInt8},
    variant_record::AbstractVector{UInt8}, offset::Integer)
    hetero_count = 0
    @inbounds for i in 1:p.header.n_samples
        if g[i] == 0x01
            hetero_count += 1
        end
    end
    phasepresent = variant_record[offset + 1] & 0x01 != 0
    phase_bits = 0
    if phasepresent
        phasepresentbv_nbytes = (hetero_count >> 3) + 1 #((hetero_count + 1) + 7) >> 3
        phase_bits = 0
        for i in 1:phasepresentbv_nbytes
            phase_bits += count_ones(variant_record[offset + i])
        end
        phase_bits -= 1 # remove initial set bit.
        offset += phasepresentbv_nbytes
    else
        phase_bits = hetero_count + 1 # should include phasepresent bit
    end
    phase_nbytes = (phase_bits + 7) >> 3
    offset += phase_nbytes
    return offset
end

# replaces the dosage in `buf` with the values in the `DiffList` `dl`.
function _get_difflist_dosage!(buf::Vector{T}, p::Pgen, dl::DiffList,
    variant_record::AbstractVector{UInt8}, offset::Integer) where T <: AbstractFloat
    ngroups = (dl.len + 63) ÷ 64
    if dl.len == 0
        return offset
    end
    dosage_unit = T == Float32 ? dosage_unit_Float32 : dosage_unit_Float64
    dosages = reinterpret(UInt16, @view(variant_record[offset + 1 : offset + 2 * dl.len]))
    for gid in 1:ngroups
        parse_difflist_sampleids!(p.difflist_cache, p.difflist_cache_incr, dl, gid)
        for (idx, sampleid) in enumerate(p.difflist_cache)
            totalidx = 64 * (gid - 1) + idx
            if totalidx > dl.len
                break
            end
            buf[sampleid] = dosages[totalidx] * dosage_unit
        end
    end
    return offset + 2 * dl.len
end


text:
# Recurrent inference machines for seismic imaging
# Same notation as Putzky and Welling paper (2017)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using JUDI.TimeModeling, JUDI4Flux, LinearAlgebra, InvertibleNetworks
using PyPlot, Random, JLD, Flux
import Flux.Optimise.update!

# Training data
D = load("../data/overthrust_images_train.jld")
ntrain = length(D["m"])

# Crop images and models to 400 x 120 (must be evenly dividable by 4)
for j=1:ntrain
    D["m"][j] = D["m"][j][1:400, 1:120]
    D["m0"][j] = D["m0"][j][1:400, 1:120]
    D["dm"][j] = D["dm"][j][1:400, 1:120]
end

# Use one sample to set up operators
m = D["m"][1]
m0 = D["m0"][1]
dm = vec(D["dm"][1])

# Set up model structure
n = size(m0)
d = (25., 25.)
o = (0., 0.)

# Setup info and model structure
nsrc = 1
model0 = Model(n, d, o, m0)

# Set up receiver geometry
nxrec = 301
xrec = range(400f0, stop=9600f0, length=nxrec)
yrec = 0f0
zrec = range(250f0, stop=250f0, length=nxrec)

# receiver sampling and recording time
time = 2000f0   # receiver recording time [ms]
dt = 1f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dt, t=time, nsrc=nsrc)

# Set up source geometry (cell array with source locations for each shot)
num_simsource = 21
xsrc = range(500f0, stop=9500f0, length=num_simsource)
ysrc = range(0f0, stop=0f0, length=num_simsource)
zsrc = range(20f0, stop=20f0, length=num_simsource)

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dt, t=time)

# setup wavelet
f0 = 0.015f0     # MHz
wavelet = zeros(Float32, srcGeometry.nt[1], num_simsource)
for j=1:num_simsource
    wavelet[:, j] = ricker_wavelet(time, dt, f0) * randn(1)[1]/sqrt(num_simsource)
end
q = judiVector(srcGeometry, wavelet)
q_sim = zeros(Float32, size(wavelet))

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model0)
info = Info(prod(n), nsrc, ntComp)

####################################################################################################

# Return data as julia array
opt = Options(return_array=true, dt_comp=1f0)

# Setup operators
F0 = judiModeling(info, model0; options=opt)
Pr = judiProjection(info, recGeometry)
Ps = judiProjection(info, srcGeometry)

# Modeling operator and Jacobian
F0 = Pr*F0*Ps'  # initial model
J = judiJacobian(F0, q)

####################################################################################################

# Objective function
function loss(L, J, d, η)

    # Initiliaze w/ zeros
    nx, ny = J.model.n
    η0 = zeros(Float32, nx, ny, 1, 1)
    s0 = zeros(Float32, nx, ny, L.L[1].C.k - 1, 1)

    # Forward pass
    η_, s_ = L.forward(η0, s0, d, J)

    # Residual and function value
    Δη = η_ - η
    f = .5f0*norm(Δη)^2

    # Backward pass (set gradients)
    L.backward(Δη, 0f0.*s0, η_, s_, d, J)

    return f
end

# Dimensions
n_in = 32
n_hidden = 64
batchsize = 1
maxiter = 2
Ψ(η) = identity(η)

# Unrolled loop
L = NetworkLoop(n[1], n[2], n_in, n_hidden, batchsize, maxiter, Ψ; type="HINT")

# Get network parameters and overwrite w/ saved values
iter_start = 1
P_curr = get_params(L)
# P_save = load(join(["network_params_iteration_", string(iter_start), ".jld"]))["P"]
# for j=1:length(P_curr)
#     P_curr[j].data = P_save[j].data
# end

# Optimization parameters
opt = Flux.ADAM(1f-3)
train_iter = 1000
indices = randperm(ntrain)

# Training loop
for j=iter_start+1:train_iter
    # Draw image + velocity from training data
    i = indices[j]
    print("Source no: ", i, "\n")
    η = D["dm"][i]
    m0 = D["m0"][i]

    # Draw random source
    for k=1:num_simsource
        q_sim[:,k] = ricker_wavelet(time, dt, f0) * randn(Float32, 1)[1]/sqrt(1f0*num_simsource)
    end

    # Generate observed data on the fly
    J.model.m = m0
    J.source[1] = q_sim
    d = J*vec(η)

    # Evaluate objective and gradients
    @time f = loss(L, J, d, η)
    print("Iteration: ", j, "; f(x) = ", f, "\n")

    # Update weights
    P = get_params(L)
    for p in P
        update!(opt, p.data, p.grad)
    end
    clear_grad!(L)

    # Save intermediate results
    if mod(j, 100) == 0
        P = get_params(L)
        save(join(["network_params_iteration_", string(j), ".jld"]), "P", P, "opt", opt)
    end
end

####################################################################################################
# Evaluate trained network for new data

# Step 1: Draw new m0, q_sim and d
i = randperm(ntrain)[1]
η = D["dm"][i]
m0 = D["m0"][i]
for k=1:num_simsource
    q_sim[:,k] = ricker_wavelet(time, dt, f0) * randn(Float32, 1)[1]/sqrt(1f0*num_simsource)
end

# Set up Jacobian
J.model.m = m0
J.source[1] = q_sim
d = J*vec(η)    # Observed data

# Step 2: Compute predicted image
η0 = zeros(Float32, n[1], n[2], 1, 1)   # zero initial guess
s0 = zeros(Float32, n[1], n[2], n_in-1, 1)
η_ = L.forward(η0, s0, d, J)[1]

# Comparison to RTM
rtm = J'*d

# Plot
η_ = η_ / norm(η_, 2)
η = η / norm(η, 2)
rtm = rtm / norm(rtm, 2)

figure(figsize=(6, 7))
subplot(3,1,1); imshow(reshape(η, model0.n)', cmap="gray", vmin=-2e-2, vmax=2e-2); title("True image")
subplot(3,1,2); imshow(reshape(rtm, model0.n)', cmap="gray", vmin=-2e-2, vmax=2e-2); title("RTM")
subplot(3,1,3); imshow(reshape(η_, model0.n)', cmap="gray", vmin=-2e-2, vmax=2e-2); title("i-RIM")

max_stars_repo_name: slimgroup/LoopUnrolledSeismicImaging

max_issues_repo_name: slimgroup/LoopUnrolledSeismicImaging

max_forks_repo_name: slimgroup/LoopUnrolledSeismicImaging


text:
mutable struct Population
    ngen # number of generations
    psze # population size
    nslc # number of selected solutions in population
    nelt # number of elites to keep for the next generation
    ffun # fitness function
    fscl # fitness scale, set to -1 for minimisation
    vars # variables to pass to fitness function
    fpop # vector of fitness over population
    bstf # best fitness so far
    nevl # number of evaluation
end

function Population(ngen, psze, nslc, nelt, ffun, fscl, vars)
    psze < nslc && throw(ArgumentError(
        "Population size `psze` must be greater than number of selected individuals `nslc`."))
    nslc < nelt && throw(ArgumentError(
        "Number of selected individuals `nslc` must be greater than number of elites `nelt`."))
    fpop = Vector{Float64}()
    bstf = -Inf
    Population(ngen, psze, nslc, nelt, ffun, fscl, vars, fpop, bstf, 0)
end

Population(f, v::AbstractVariable; ngen, psze, nslc, nelt, fscl = 1) =
    Population(ngen, psze, nslc, nelt, f, fscl, VariableSet(v))

Population(f, v::Tuple; ngen, psze, nslc, nelt, fscl = 1) =
    Population(ngen, psze, nslc, nelt, f, fscl, VariableSet(v...))

Population(f, v::VariableSet; ngen, psze, nslc, nelt, fscl = 1) =
    Population(ngen, psze, nslc, nelt, f, fscl, v)

function update_population(p)
    i = sortperm(p.fpop, rev = true)
    b = i[1]
    s = i[1:p.nslc]
    e = i[1:p.nelt]
    for v ∈ p.vars
        update_variable!(v, b, s, e)
    end
    if p.fpop[b] > p.bstf
        p.bstf = p.fpop[b]
    end
    p.fpop = p.fpop[e]
    p
end

# add option: cb before or after update!(o)
# 1 seule callback
function optimise!(p; cbbs = (p) -> (), cbas = (p) -> (), mt = false)
    isempty(p.fpop) || throw(ArgumentError("A population can not be optimised more than once."))
    for _ ∈ 1:p.ngen
        psze = isempty(p.fpop) ? p.psze : p.psze - p.nelt
        evls = []
        if mt
            for _ ∈ 1:Threads.nthreads()
                push!(evls, [])
            end
            Threads.@threads for _ ∈ 1:psze
                x = [rand(v) for v ∈ p.vars]
                f = p.fscl * p.ffun(x...)
                push!(evls[Threads.threadid()], (x, f))
            end
            evls = vcat(evls...)
        else
            for _ ∈ 1:psze
                x = [rand(v) for v ∈ p.vars]
                f = p.fscl * p.ffun(x...)
                push!(evls, (x, f))
            end
        end
        for e ∈ evls
            x, f = e
            for (v, s) ∈ zip(p.vars, x)
                push!(v.vpop, s)
            end
            push!(p.fpop, f)
            p.nevl += 1
        end
        cbbs(p)
        update_population(p)
        cbas(p)
        # if termination criteria break
    end
    p
end

function optimum(p::Population)
    if length(p.vars) == 1
        return (p.vars.bstv, p.bstf / p.fscl)
    end
    nm = names(p.vars)
    bv = tuple([v.bstv for v ∈ p.vars]...)
    (NamedTuple{nm}(bv), p.bstf / p.fscl)
end

# models

# current state of population : fitness values and variable values (summary)
# (dictionnary : fitness => fpop, name = var)
# :fitness : mean(fpop ./ fscl) std(fpop ./ fscl)
# :variables : :mean :std :model
# model : DiscreteVariable v.modl * v.wdth (arithmetic on didtributions are allowed)
#         DiscreteBoundedVariable v.modl * v.wdth + v.lbnd
function state(p)
    model(v::AbstractVariable) = v.modl
    model(v::DiscreteVariable) = v.modl * v.wdth
    model(v::DiscreteBoundedVariable) = v.modl * v.wdth + v.lbnd
    nm = []
    st = []
    for (n, v) ∈ zip(names(p.vars), p.vars)
        push!(nm, n)
        #push!(st, v.vpop)
        m = model(v)
        push!(st, (mean = mean(m), std = std(m), model = m))
    end
    push!(nm, :fitness)
    fpop = p.fpop ./ p.fscl
    push!(st, (mean = mean(fpop), std = std(fpop)))
    NamedTuple{tuple(nm...)}(tuple(st...))
end

neval(p) = p.nevl


text:
export XorEnv, fitness, new_gen!

using NeuroEvolution

"Random BERL environment: returns random fitness"
mutable struct XorEnv <: BERLenv
    name::String
    memory::Array # Array of storable data
    X::Array{Array{Float64}}
    y::Array{Array{Float64}}
end

environments["xor"]=XorEnv

function XorEnv(cfg::Dict)
    cfg["env"]="XOR"
    cfg["n_in"]=2
    cfg["n_out"]=1
    X, y = xor_dataset(2, 100)
    XorEnv("XOR", [], X, y)
end

"Returns the fitness of one Cambrian Individual in this environment"
function fitness(indiv::Cambrian.Individual, env::XorEnv)
    y_pred = []
    for x in env.X
        push!(y_pred, process(indiv, x))
    end
    [sum(log_fitness.(env.y, y_pred)) / length(env.X)]
end

"Sets up the env for the next generation"
function new_gen!(env::XorEnv, cfg::Dict=nothing)
    env.X, env.y = xor_dataset(2, 100) # Generate new data
end

## Utils
function xor(a::Float64, b::Float64)
    if a + b == 1
        1.
    else
        0.
    end
end

"Array to array XOR"
function xor(a::Array{Float64})
    v = a[1]
    for i in 2:length(a)
        v = xor(v, a[i])
    end
    [v]
end

function rand_bin(len::Int64)
    1.0 * rand(0:1, len)
end

function xor_dataset(len::Int64, n_records::Int64)
    X::Array{Array{Float64}}=[]
    y::Array{Array{Float64}}=[]
    for i in 1:n_records
        l = rand_bin(len)
        push!(X, l)
        push!(y, xor(l))
    end
    X, y
end

function log_loss(y_true::Float64, y_pred::Float64)
    y_pred = maximum([minimum([y_pred, 1-10^-15]), 10^-15])
    if y_true == 1.0
        -log(y_pred)
    else
        -log(1-y_pred)
    end
end

function log_loss(y_true::Int64, y_pred::Int64)
    log_loss(float(y_true), float(y_pred))
end

function log_loss(y_true::Array, y_pred::Array)
    log_loss(y_true[1], y_pred[1])
end

function log_fitness(y_true, y_pred)
    max_loss = log_loss(y_true, 1 .- y_true)
    max_loss - log_loss(y_true, y_pred)
end


text:
using Test, ArgParse
using SignalingDimension

try print_test_results catch err
    include("../utilities.jl")
    using .ScriptUtilities: print_test_results
end

# parsing command line arguments
arg_table = ArgParseSettings()
@add_arg_table! arg_table begin
    "-Y"
        help = "max number of outputs to consider"
        arg_type = Int64
        default = 45
    "--dir"
        help = "directory to print results .txt to"
        arg_type = String
        default = "./script/facet_verifications/results/"
    "--stdout"
        help = "Print output to STDOUT. If not flagged, output is printed to `.txt` file."
        action = :store_true
end
args = parse_args(ARGS, arg_table)

# verification script: verifies all ambiguous guessing facets with Y ∈ [4, Y_max]
function verify_ambiguous_guessing_facet(Y_max :: Int64)
    @testset "Testing all maximum likelihood facets with `Y ∈ [4, $Y_max]` or less." begin
        @testset "scanning over all `d ∈ [2, $Y - 2]`" for Y in 4:Y_max
            @testset "Verifying Y = $Y and d = $d" for d in 2:Y-2
                vertices = aff_ind_ambiguous_guessing_vertices(Y, d)
                facet = ambiguous_guessing_facet(Y, d)

                @test verify_facet(d, facet, vertices)
            end
        end
    end
end

# outputting result to file or stdout
print_test_results(
    verify_ambiguous_guessing_facet,
    params=[args["Y"]],
    stdout=args["stdout"],
    dir=args["dir"]
)


text:
module DataUtil

export DataParam, make_dataset, centering, train_val_split, KFold

using Random, Statistics, Distributions
using LinearAlgebra
using InvertedIndices
using Parameters

@with_kw struct DataParam{A,B,C,D,E,F}
    N::A
    β₀::B
    𝛃::C
    𝚺::D = I
    𝛍ₓ::E = zeros(length(𝛃)); @assert length(𝛍ₓ) == length(𝛃)
    𝚺ₓ::F = I
end

# make_dataset
function make_dataset(β₀, 𝛃, N)
    d = length(𝛃)
    𝐗 = randn(N, d)
    𝐲 = β₀ .+ 𝐗 * 𝛃 .+ randn(N)
    return (𝐗, 𝐲)
end

function make_dataset(param::DataParam)
    @unpack N, β₀, 𝛃, 𝚺, 𝛍ₓ, 𝚺ₓ = param
    𝐗 = rand(MvNormal(𝛍ₓ, 𝚺ₓ), N) |> transpose
    𝐲 = rand(MvNormal(β₀ .+ 𝐗 * 𝛃, 𝚺))
    return (𝐗, 𝐲)
end

# centering
function centering(𝐗::AbstractMatrix)
    X̄ = mean(𝐗, dims=1)
    return 𝐗 .- X̄
end

function centering(𝐲::AbstractVector)
    ȳ = mean(𝐲)
    return 𝐲 .- ȳ
end

centering(𝐗, 𝐲) = centering(𝐗), centering(𝐲)

# train_val_split
function train_val_split(𝐗::AbstractMatrix, val_begin, val_end; do_centering=true)
    𝐗_train, 𝐗_val = @views 𝐗[Not(val_begin:val_end), :], 𝐗[val_begin:val_end, :]
    return do_centering ? centering(𝐗_train, 𝐗_val) : (𝐗_train, 𝐗_val)
end

function train_val_split(𝐲::AbstractVector, val_begin, val_end; do_centering=true)
    𝐲_train, 𝐲_val = @views 𝐲[Not(val_begin:val_end)], 𝐲[val_begin:val_end]
    return do_centering ? centering(𝐲_train, 𝐲_val) : (𝐲_train, 𝐲_val)
end

function train_val_split(𝐗::AbstractMatrix; ratio=0.5, do_centering=true)
    N = size(𝐗, 1)
    train_size = round(Int, N * ratio)
    return train_val_split(𝐗, train_size + 1, N; do_centering=do_centering)
end

function train_val_split(𝐲::AbstractVector; ratio=0.5, do_centering=true)
    N = length(𝐲)
    train_size = round(Int, N * ratio)
    return train_val_split(𝐲, train_size + 1, N; do_centering=do_centering)
end

# KFold
struct KFold{S,T <: AbstractArray}
    𝑘::S # 𝑘-fold
    parts::T # partition index list (the size is 𝑘+1, parts[begin]=1, parts[end]=N)
end

KFold(N::Integer, 𝑘::Integer) = KFold(𝑘, [min(1 + i * ceil(Int, N / 𝑘), N + 1) for i = 0:𝑘])

end # module DataUtil


text:
module OptMethods
using TargetFunction

function Fletcher_Reeves(

    start_point, fun::TargetFunction.ObjFun, iteration=10, alpha=0.1)

    x = start_point
    track = x

    k = 1
    while k < iteration

        if k < 2

            g = fun.g_x(x)
            beta = 0.0
            s = -g
        else

            g = fun.g_x(x)
            m_old = (g_old'*g_old)[1]
            m_new = (g'*g)[1]
            beta = m_new / m_old
        end

        s = -g + beta * s
        alpha_k = _backtracking_line_search(fun, x, s)
        g_old = fun.g_x(x)
        x = x + alpha_k * s

        k += 1
        track = [track x]
    end
    return track
end
function _backtracking_line_search(fun::TargetFunction.ObjFun, x, s, c=1e-4)

    alpha = 1.0
    p = 0.6
    c = 0.0001
    diff = 1

    f_alpha = fun.about_alpha(x, s)
    g_alpha = fun.about_alpha_prime(x, s)
    f_alpha_0 = f_alpha(0.0)
    g_alpha_0 = g_alpha(0.0)

    i = 0

    while (diff > 0.0 && i < 1000)

        f_x = fun.f_x(x)
        diff = f_alpha(alpha) - f_alpha_0 - c * alpha * g_alpha_0
        diff = diff[1]
        alpha *= p
        i += 1
    end

    return alpha
end

end

max_stars_repo_name: wyli/gradient

max_issues_repo_name: wyli/gradient

max_forks_repo_name: wyli/gradient


text:
using Aqua
using ProbabilisticCircuits
using Test

@testset "Aqua tests" begin
    Aqua.test_all(ProbabilisticCircuits,
                    ambiguities = false)
end


text:
using Pkg
Pkg.activate(".")

using Reproduce

function main()

    as = ArgParseSettings()
    @add_arg_table as begin
        "config"
        arg_type=String
        "--path"
        arg_type=String
        default=""
        "--numworkers"
        arg_type=Int
        default=4
        "--numjobs"
        action=:store_true
    end
    parsed = parse_args(as)

    experiment = Reproduce.parse_experiment_from_config(parsed["config"], parsed["path"])

    pre_experiment(experiment)
    ret = job(experiment; num_workers=parsed["numworkers"])
    post_experiment(experiment, ret)

end

main()


text:
function compute_c_n(u_n::Array{T,1}) where {T <: Real}
# Transform the u_n coefficients to c_n, which are coefficients
# of the basis in which the P(G_n) functions are computed.
n = length(u_n)
c_n = zeros(T,n+3)
a_n = zeros(T,n+1)
# Pre-allocated memory - just need to zero terms:
a_n[1] = one(T)  # Add in the first constant coefficient term
for i=1:n
  # Compute the contribution to a_n*\mu^n
  bcoeff = one(T)
  for j=0:i
#    a_n[j+1] -= u_n[i]*binomial(i,j)*(-1)^j
    a_n[j+1] -= u_n[i]*bcoeff*(-1)^j
    bcoeff *= (i-j)/(j+1)
#    println("i: ",i," j: ",j," a_i: ",a_n[j+1])
  end
end
# Now, compute the c_n coefficients:
for j=n:-1:2
  c_n[j+1] = a_n[j+1]/(j+2)+c_n[j+3]
end
c_n[2] = a_n[2]+3*c_n[4]
c_n[1] = a_n[1]+2*c_n[3]
return c_n[1:n+1]
end

function compute_c_n_grad(u_n::Array{T,1}) where {T <: Real}
# Transform the u_n coefficients to c_n, which are coefficients
# of the basis in which the P(G_n) functions are computed.
# Compute the derivatives of the flux with respect to the u coefficients.
n = length(u_n)
# We define c_n with two extra elements which are zero:
c_n = zeros(T,n+3)
a_n = zeros(T,n+1)
dadu = zeros(T,n+1,n)
dcdu = zeros(T,n+3,n)
a_n[1] = one(T)  # Add in the first constant coefficient term
for i=1:n
  # Compute the contribution to a_n*\mu^n
  bcoeff = one(T)
  for j=0:i
#    a_n[j+1] -= u_n[i]*binomial(i,j)*(-1)^j
    a_n[j+1] -= u_n[i]*bcoeff*(-1)^j
#    dadu[j+1,i] -= binomial(i,j)*(-1)^j
    dadu[j+1,i] -= bcoeff*(-1)^j
    bcoeff *= (i-j)/(j+1)
#    println("i: ",i," j: ",j," a_i: ",a_n[j+1])
  end
end
# Now, compute the c_n coefficients and propagate derivatives:
for j=n:-1:2
  c_n[j+1] = a_n[j+1]/(j+2)+c_n[j+3]
  for i=1:n
    dcdu[j+1,i] = dadu[j+1,i]/(j+2) + dcdu[j+3,i]
  end
end
c_n[2] = a_n[2]+3*c_n[4]
for i=1:n
  dcdu[2,i] = dadu[2,i] + 3*dcdu[4,i]
end
c_n[1] = a_n[1]+2*c_n[3]
for i=1:n
  dcdu[1,i] = dadu[1,i] + 2*dcdu[3,i]
end
return c_n[1:n+1],dcdu[1:n+1,:]
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
using Test
using ABM4OSN
using LightGraphs

n = ABM4OSN.create_network(100, 5)

@test nv(n) == 100
@test is_directed(n)
@test eltype(n) == Int64


text:
include("test_CG.jl")
include("test_FFTOp.jl")
include("test_MWNI.jl")
include("test_POCS.jl")
using SeisProcessing


text:
"""
```
cols = im2col(img, blocksize, [blocktype or strides])
```

creates a copy of all image blocks of size `blocksize` and arranges them as
the columns of output array `cols`. The third parameter can either be a tuple
of `strides` or a `blocktype` string (`"sliding"`, `"distinct"`); the default
value is `"sliding"`, which corresponds to a stride of 1 in all dimensions.

This function is similar to Matlab's `im2col`, but supports N-D images.
See http://www.mathworks.com/help/images/ref/im2col.html
"""

function im2col{T,N}(A::AbstractArray{T,N}, sz::NTuple{N,Int}, blocktype::AbstractString="sliding")
  if blocktype == "sliding"
    return im2col(A, sz, tuple(ones(Int,N)...))
  elseif blocktype == "distinct"
    return im2col(A, sz, sz)
  else
    throw(ArgumentError("Invalid 'blocktype' (must be 'sliding' or 'distinct')."))
  end
end

function im2col{T,N}(A::AbstractArray{T,N}, sz::NTuple{N,Int}, stride::NTuple{N,Int})
  check_im2col_col2im(size(A),sz,stride)
  ncols::Int = mapreduce( d -> 1 + div(size(A,d)-sz[d], stride[d]), *, 1:N)
  im2col!(Array(T,prod(sz),ncols), A, sz, stride)
end

"""
```
im2col!(cols, img, blocksize, [strides])
```

stores all image blocks of size `blocksize` in the pre-allocated output `cols`.

See also: `im2col`.
"""

@generated function im2col!{T,N}(cols::AbstractArray{T,2}, A::AbstractArray{T,N}, sz::NTuple{N,Int}, stride::NTuple{N,Int}=tuple(ones(Int,N)...))
  quote
    check_im2col_col2im(size(A),sz,stride,size(cols))
    c = 0
    @inbounds begin
      @nloops $N i d->1:stride[d]:(size(A,d)-sz[d]+1) begin
        r = 0; c += 1
        @nloops $N j d->0:(sz[d]-1) begin
          r += 1
          cols[r,c] = @nref $N A d->i_d+j_d
        end
      end
    end
    cols
  end
end

function check_im2col_col2im{N}(szA::NTuple{N,Int}, szB::NTuple{N,Int}, stride::NTuple{N,Int}, szcols::Union{Void,NTuple{2,Int}}=nothing)
  all(d -> stride[d] >= 1, 1:N)                    || throw(DimensionMismatch("All strides must be >= 1."))
  all(d -> szA[d] >= szB[d] >= 1, 1:N)             || throw(DimensionMismatch("Window size too small or large."))
  all(d -> mod(szA[d]-szB[d],stride[d]) == 0, 1:N) || throw(DimensionMismatch("Image and window sizes not compatible with strides."))
  if szcols != nothing
    ncols = mapreduce( d -> 1 + div(szA[d]-szB[d], stride[d]), *, 1:N)
    szcols == (prod(szB),ncols) || throw(DimensionMismatch("Image and window sizes not compatible with 'cols' matrix."))
  end
  nothing
end

"""
```
A = col2im(cols, blocksize, sizeA, [strides, accu_op, accu_init, norm_op] or [blocktype, overlap])
```

uses each column of the matrix `cols` as a block of size `blocksize` to
form the output image `A` of size `sizeA`.

There are two sets of optional parameters:

- `strides`:   tuple of block strides (default is a stride of 1 in all dimensions)
- `accu_op`:   function to accumulate values if blocks overlap (default: +)
- `accu_init`: initial value for all pixels of output `A` (default: 0)
- `norm_op`:   function to normalize each pixel with its count of overlapping blocks (default: /)

or

- `blocktype`: string to denote strides, can be `"sliding"` (default: all strides = 1)
               or `"distinct"` (strides = blocksize)
- `overlap`:   string to denote overlap policy (`accu_op`,`accu_init`,`norm_op`),
               can be `"sum"` (default: +,0,nothing) or `"average"` (+,0,/)

This function is similar to Matlab's `col2im`, but supports N-D images and overlapping blocks.
See http://www.mathworks.com/help/images/ref/col2im.html
"""

function col2im{T,N}(cols::AbstractArray{T,2}, sz::NTuple{N,Int}, szA::NTuple{N,Int}, blocktype::AbstractString="sliding", overlap::AbstractString="average")
  if blocktype == "sliding"
    overlap in ("sum","average") || throw(ArgumentError("Invalid 'overlap' (must be 'sum' or 'average')."))
    return col2im(cols, sz, szA, tuple(ones(Int,N)...), +, zero(T), overlap == "average" ? (/) : nothing)
  elseif blocktype == "distinct"
    return col2im(cols, sz, szA, sz)
  else
    throw(ArgumentError("Invalid 'blocktype' (must be 'sliding' or 'distinct')."))
  end
end

function col2im{T,N}(cols::AbstractArray{T,2}, sz::NTuple{N,Int}, szA::NTuple{N,Int}, stride::NTuple{N,Int}, accu_op::Function=+, accu_init=zero(T), norm_op::Union{Void,Function}=/)
  check_im2col_col2im(szA,sz,stride,size(cols))
  do_normalize = any(d -> stride[d] < sz[d], 1:N) && norm_op != nothing
  !do_normalize || isa(norm_op(cols[1],1),T) ||
    (U = typeof(norm_op(cols[1],1)); throw(ArgumentError("'norm_op' causes type change from $T to $U: try changing the type of 'cols' to $U beforehand.")))
  A = col2im!(fill!(Array{T}(szA),convert(T,accu_init)), cols, sz, stride, accu_op)
  do_normalize && broadcast!(norm_op, A, A, overlap_counts(szA,sz,stride))
  A
end

overlap_counts{N}(szA::NTuple{N,Int}, sz::NTuple{N,Int}, stride::NTuple{N,Int}) =
  col2im( im2col(ones(Int,szA),sz,stride), sz, szA, stride, +, 0, nothing)

"""
```
col2im!(A, cols, blocksize, [strides, accu_op])
```

rearranges the columns of `cols` as blocks of size `blocksize` in the pre-allocated output `A`.

See also: `col2im`.
"""

col2im!{T,N}(A::AbstractArray{T,N}, cols::AbstractArray{T,2}, sz::NTuple{N,Int}, stride::NTuple{N,Int}=tuple(ones(Int,N)...), op::Function=+) =
  col2im!(A,cols,sz,stride,Type{symbol(op)})

@generated function col2im!{T,N,op}(A::AbstractArray{T,N}, cols::AbstractArray{T,2}, sz::NTuple{N,Int}, stride::NTuple{N,Int}, ::Type{Type{op}})
  quote
    check_im2col_col2im(size(A),sz,stride,size(cols))
    c = 0
    @inbounds begin
      @nloops $N i d->1:stride[d]:(size(A,d)-sz[d]+1) begin
        r = 0; c += 1
        @nloops $N j d->0:(sz[d]-1) begin
          r += 1
          (@nref $N A d->i_d+j_d) = $op((@nref $N A d->i_d+j_d), cols[r,c])
        end
      end
    end
    A
  end
end


text:
#=
Copyright 2021 BlackRock, Inc.

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.
=#

using PiecewiseQuadratics
using Documenter
using Test

@testset "Doctests" begin
    DocMeta.setdocmeta!(PiecewiseQuadratics, :DocTestSetup, :(using PiecewiseQuadratics);
                        recursive=true)
    doctest(PiecewiseQuadratics)
end

@testset "Intervals" begin
    include("test_intervals.jl")
end
@testset "Bounded Quadratics" begin
    include("test_bounded_quadratics.jl")
end
@testset "Piecewise Quadratics" begin
    include("test_piecewise_quadratics.jl")
end
@testset "Convex envelope" begin
    include("test_envelope.jl")
end


text:
function train!(model, X, y, opt=Adam())
    # ps = Flux.params(model)
    return loss(X, y) = ELBO(model, X, y)
    # Flux.train!(loss,ps,(X,y),opt)
end


text:
module BlockHankelMatrices

using BlockArrays

include("typedef.jl")

export BlockHankel

end


text:
using StippleKatex
using Test

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


text:
using MPRebinSpectra
using Test

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


text:
##################################################
function symbols(x::AbstractString; kwargs...)
    out = sympy.symbols(x; kwargs...)
end
symbols(x::Symbol; kwargs...) = symbols(string(x); kwargs...)
symbols(xs::T...; kwargs...) where {T <: SymbolicObject} = xs

"""
    @vars x y z

Define symbolic values, possibly with assumptions

Examples:
```
@vars x y
@vars a b real=true
```

"""
macro vars(x...)
    q = Expr(:block)
    as = []    # running list of assumptions to be applied
    ss = []    # running list of symbols created
    for s in reverse(x)
        if isa(s, Expr)    # either an assumption or a named variable
            if s.head == :(=)
                s.head = :kw
                push!(as, s)
            elseif s.head == :(=>)
                push!(ss, s.args[1])
                push!(q.args, Expr(:(=), esc(s.args[1]), Expr(:call, :symbols, s.args[2], map(esc,as)...)))
            end
        elseif isa(s, Symbol)   # raw symbol to be created
            push!(ss, s)
            push!(q.args, Expr(:(=), esc(s), Expr(:call, :symbols, Expr(:quote, s), map(esc,as)...)))
        else
            throw(AssertionError("@syms expected a list of symbols and assumptions"))
        end
    end
    push!(q.args, Expr(:tuple, map(esc,reverse(ss))...)) # return all of the symbols we created
    q
end

# old alias for @vars; should deprecate...
macro syms(x...)
    q = Expr(:block)
    as = []    # running list of assumptions to be applied
    ss = []    # running list of symbols created
    for s in reverse(x)
        if isa(s, Expr)    # either an assumption or a named variable
            if s.head == :(=)
                s.head = :kw
                push!(as, s)
            elseif s.head == :(=>)
                push!(ss, s.args[1])
                push!(q.args, Expr(:(=), esc(s.args[1]), Expr(:call, :symbols, s.args[2], map(esc,as)...)))
            end
        elseif isa(s, Symbol)   # raw symbol to be created
            push!(ss, s)
            push!(q.args, Expr(:(=), esc(s), Expr(:call, :symbols, Expr(:quote, s), map(esc,as)...)))
        else
            throw(AssertionError("@syms expected a list of symbols and assumptions"))
        end
    end
    push!(q.args, Expr(:tuple, map(esc,reverse(ss))...)) # return all of the symbols we created
    q
end

## avoid PyObject conversion as possible
Sym(x::T) where {T <: Number} = sympy.sympify(x)
Sym(x::Rational{T}) where {T} = Sym(numerator(x))/Sym(denominator(x))
function Sym(x::Complex{Bool})
    !x.re && x.im && return IM
    !x.re && !x.im && return zero(Sym)
    x.re && !x.im && return Sym(1)
    x.re && x.im && return Sym(1) + IM
end
Sym(xs::Symbol...) = Tuple(Sym.((string(x) for x in xs)))
Sym(x::AbstractString) = sympy.symbols(x)
Sym(s::SymbolicObject) = s
Sym(x::Irrational{T}) where {T} = convert(Sym, x)

convert(::Type{Sym}, s::AbstractString) = Sym(s)

sympify(s, args...; kwargs...) = sympy.sympify(s, args...; kwargs...)

SymMatrix(s::SymMatrix) = s
SymMatrix(s::Sym) = sympy.ImmutableMatrix([s])
SymMatrix(A::Matrix) = sympy.ImmutableMatrix([A[i,:] for i in 1:size(A)[1]])


text:
using Random

sym(x) = Symbol(x)

function populate_dependency!(doc, comp)
    empty!(comp.deps)
    for x in list_encapsulation(doc, comp)
        push!(comp.deps, sym(x["component"]))
    end
end

function add_component!(doc, name, node, populate=true)
    comp = Component(name, node, Set{Symbol}())
    populate && populate_dependency!(doc, comp)
    push!(doc.comps, comp)
    comp
end

function add_connection!(doc, c1, c2, vars)
    conn = Connection(c1, c2, vars)
    push!(doc.conns, conn)
    conn
end

"""
    load_cellml loads a CellML file and populates a Document structures with
    the components and connections.

    If resolve = true, load_cellml also resolve the imported components.
"""
function load_cellml(path; resolve=true)
    xml = readxml(path)
    doc = Document(path, [xml], Component[], Connection[])

    for c in list_components(xml)
        add_component!(doc, sym(c["name"]), c)
    end

    for conn in list_connections(xml)
        c1, c2 = sym.(components_of(get_connection_component(conn)))
        vars = [make_var(c1,v1) => make_var(c2,v2)
                for (v1,v2) in variables_of.(list_connection_variables(conn))]
        add_connection!(doc, c1, c2, vars)
    end

    resolve && resolve_imports!(doc)

    return doc
end

"""
    find_component finds a component by name
"""
function find_component(comps::Array{Component}, name)
    i = findfirst(x -> x.name == Symbol(name), comps)
    return i == nothing ? nothing : comps[i]
end

find_component(doc::Document, name) = find_component(doc.comps, name)

# implicit_name(sym) = Symbol("$(string(sym))_$(randstring(6))")
implicit_name(sym) = Symbol("$(string(sym))")

"""
    resolve_imports recursivelly resolves the imported components of doc.
"""
function resolve_imports!(doc::Document)
    for ϵ in list_imports(doc)  # εισαγωγή == import
        href = ϵ["xlink:href"]
        @info "importing $href"
        path = joinpath(splitdir(doc.path)[1], href)
        child = load_cellml(path)
        append!(doc.xmls, child.xmls)

        L = Dict{Symbol, Symbol}()      # L keeps the import names of components
        D = Set{Symbol}()               # D is the set of imported dependencies
        C = Set{Symbol}()

        # first, import explicitly imported component
        for σ in list_import_components(ϵ)  # συστατικό == component
            old_name = sym(σ["component_ref"])
            new_name = sym(σ["name"])
            old_comp = find_component(child, old_name)

            add_component!(doc, new_name, old_comp.node)

            D = D ∪ old_comp.deps
            L[old_name] = new_name
            push!(C, old_name)
        end

        # second, import implicitly imported component, as defined by closure
        #for c in find_closure(child, collect(keys(L)))
        for old_name in D
            if old_name ∉ C
                old_comp = find_component(child, old_name)
                new_name = implicit_name(old_name) # generates a name with a random suffix to prevent name collision
                add_component!(doc, new_name, old_comp.node, false)
                L[old_name] = new_name
            end
        end

        # third, we also need to import connections from the child CellML file
        for conn in connections(child)
            c1, c2 = components(conn)
            if haskey(L,c1) && haskey(L,c2)
                d1, d2 = L[c1], L[c2]
                vars = [Var(d1,v1.var) => Var(d2,v2.var) for (v1,v2) in conn.vars]
                add_connection!(doc, d1, d2, vars)
            end
        end
    end
end

"""
    find_closure finds the transitive closure of a list of componenets (l) minus
    the list itself, i.e., it returns the list of components in doc which are
    reachable through a chain of connections starting from any component in l.
"""
function find_closure(doc::Document, l)
    n = copy(l)

    done = false
    while !done
        done = true
        for (c1,c2) in components.(connections(doc))
            if c1 ∈ n && c2 ∉ n
                push!(n, c2)
                done = false
            elseif c1 ∉ n && c2 ∈ n
                push!(n, c1)
                done = false
            end
        end
    end

    [x for x in n if x ∉ l]
end

"""
Given a path to a directory containing multiple CellML files.

Returns the top level files.

Used as a helper function with the CellML Physiome Model repositories.
"""
function list_top_cellml_files(dir)
    files = filter(f->endswith(f, ".cellml"), readdir(dir))

    if length(files) == 1
        return joinpath.(dir, files)
    end

    imported = Set{String}()
    for f in files
        xml = readxml(joinpath(dir, f))
        for n in CellMLToolkit.list_imports(xml)
            push!(imported, n["xlink:href"])
        end
    end

    if !isempty(imported)
        printstyled("imported files are $imported\n"; color=:yellow)
    end

    joinpath.(dir, [f for f in files if f ∉ imported])
end


text:
export
    TranslateFilter

mutable struct TranslateFilter <: Ranger.AbstractNode
    base::Nodes.NodeData
    transform::Nodes.TransformProperties{Float64}

    # Collection of nodes.
    children::Array{Ranger.AbstractNode,1}

    # This node's parent translation component
    translate_component::Math.AffineTransform{Float64}

    function TranslateFilter(world::Ranger.World, name::String, parent::Ranger.AbstractNode)
        o = new()

        o.base = Nodes.NodeData(name, parent, world)
        o.transform = Nodes.TransformProperties{Float64}()
        o.children = Array{Ranger.AbstractNode,1}[]
        o.translate_component = Math.AffineTransform{Float64}()

        o
    end
end

# Filters are special in that they overload the visit() method
# Because this is a translate filter we "filter out" everything
# but translation component from the immediate parent.
function Nodes.visit(node::TranslateFilter, context::Rendering.RenderContext, interpolation::Float64)
    if !Nodes.is_visible(node)
        return;
    end

    Rendering.save!(context)

    children = Nodes.get_children(node)
    if children ≠ nothing
        for child in children
            Rendering.save!(context)

            if Nodes.has_parent(node)
                parent = node.base.parent
                inverse = parent.transform.inverse

                # This removes the immediate parent's transform effects
                Rendering.apply!(context, inverse)

                # Re-introduce only the parent's translation component by
                # excluding the other components
                Nodes.calc_filtered_transform!(parent.transform, false, true, true, node.translate_component)

                Rendering.apply!(context, node.translate_component)
            else
                println("TranslateFilter::visit: ", node, " has NO parent")
                return;
            end

            # Now visit the child with the modified context
            Nodes.visit(child, context, interpolation)

            Rendering.restore!(context)
        end
    end

    Rendering.restore!(context)
end

# --------------------------------------------------------
# Grouping
# --------------------------------------------------------
function Nodes.get_children(node::TranslateFilter)
    node.children
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
mutable struct MapProposeState{T}
    choices::DynamicChoiceMap
    weight::Float64
    retvals::Vector{T}
end

function process_new!(gen_fn::Map{T,U}, args::Tuple, key::Int,
                      state::MapProposeState{T}) where {T,U}
    local subtrace::U
    kernel_args = get_args_for_key(args, key)
    (submap, weight, retval) = propose(gen_fn.kernel, kernel_args)
    set_submap!(state.choices, key, submap)
    state.weight += weight
    state.retvals[key] = retval
end

function propose(gen_fn::Map{T,U}, args::Tuple) where {T,U}
    len = length(args[1])
    choices = choicemap()
    state = MapProposeState{T}(choices, 0., Vector{T}(undef,len))
    for key=1:len
        process_new!(gen_fn, args, key, state)
    end
    (state.choices, state.weight, PersistentVector{T}(state.retvals))
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

const curmod = current_module()
const curmod_name = fullname(curmod)
const curmod_prefix = "$(["$m." for m in curmod_name]...)"

using Base.Test

@test_throws MethodError convert(Enum, 1.0)

@test_throws ArgumentError eval(:(@enum Foo))

@enum Fruit apple orange kiwi
@test typeof(Fruit) == DataType
@test isbits(Fruit)
@test typeof(apple) <: Fruit <: Enum
@test Int(apple) == 0
@test Int(orange) == 1
@test Int(kiwi) == 2
@test Fruit(0) == apple
@test Fruit(1) == orange
@test Fruit(2) == kiwi
@test_throws ArgumentError Fruit(3)
@test_throws ArgumentError Fruit(-1)
@test Fruit(0x00) == apple
@test Fruit(big(0)) == apple
@test_throws MethodError Fruit(0.0)
@test typemin(Fruit) == apple
@test typemax(Fruit) == kiwi
@test convert(Fruit,0) == apple
@test convert(Fruit,1) == orange
@test convert(Fruit,2) == kiwi
@test_throws ArgumentError convert(Fruit,3)
@test_throws ArgumentError convert(Fruit,-1)
@test convert(UInt8,apple) === 0x00
@test convert(UInt16,orange) === 0x0001
@test convert(UInt128,kiwi) === 0x00000000000000000000000000000002
@test typeof(convert(BigInt,apple)) <: BigInt
@test convert(BigInt,apple) == 0
@test convert(Bool,apple) == false
@test convert(Bool,orange) == true
@test_throws InexactError convert(Bool,kiwi)
@test instances(Fruit) == (apple, orange, kiwi)

f(x::Fruit) = "hey, I'm a Fruit"
@test f(apple) == "hey, I'm a Fruit"

d = Dict(apple=>"apple",orange=>"orange",kiwi=>"kiwi")
@test d[apple] == "apple"
@test d[orange] == "orange"
@test d[kiwi] == "kiwi"
vals = [apple,orange,kiwi]
for (i,enum) in enumerate(instances(Fruit))
    @test enum == vals[i]
end

@enum(QualityofFrenchFood, ReallyGood)
@test length(instances(QualityofFrenchFood)) == 1
@test typeof(ReallyGood) <: QualityofFrenchFood <: Enum
@test Int(ReallyGood) == 0

@enum Binary _zero=0 _one=1 _two=10 _three=11
@test Int(_zero) === 0
@test Int(_one) === 1
@test Int(_two) === 10
@test Int(_three) === 11
@enum Negative _neg1=-1 _neg2=-2
@test Int(_neg1) === -1
@test Int(_neg2) === -2
@test_throws InexactError convert(UInt8, _neg1)
@enum Negative2 _neg5=-5 _neg4 _neg3
@test Int(_neg5) === -5
@test Int(_neg4) === -4
@test Int(_neg3) === -3

@test_throws ArgumentError eval(:(@enum Test1 _zerofp=0.0))
@test_throws ArgumentError eval(:(@enum Test11 _zerofp2=0.5))
@enum Test111 _zerobi=BigInt(1)
@test Integer(_zerobi) == 1

# can't use non-identifiers as enum members
@test_throws ArgumentError eval(:(@enum(Test2, ?)))
@test_throws ArgumentError eval(:(@enum Test22 1=2))

# other Integer types of enum members
# TODO - not yet supported
#=
@enum Test3 _one_Test3=0x01 _two_Test3=0x02 _three_Test3=0x03
@test typeof(_one_Test3.val) <: UInt8
@test _one_Test3.val === 0x01
@test length(instances(Test3)) == 3

@enum Test4 _one_Test4=0x01 _two_Test4=0x0002 _three_Test4=0x03
@test _one_Test4.val === 0x0001
@test _two_Test4.val === 0x0002
@test _three_Test4.val === 0x0003
@test typeof(_one_Test4.val) <: UInt16

@enum Test5 _one_Test5=0x01 _two_Test5=0x00000002 _three_Test5=0x00000003
@test _one_Test5.val === 0x00000001
@test _two_Test5.val === 0x00000002
@test _three_Test5.val === 0x00000003
@test typeof(_one_Test5.val) <: UInt32

@enum Test6 _one_Test6=0x00000000000000000000000000000001 _two_Test6=0x00000000000000000000000000000002
@test _one_Test6.val === 0x00000000000000000000000000000001
@test _two_Test6.val === 0x00000000000000000000000000000002
@test typeof(_one_Test6.val) <: UInt128

@enum Test7 _zero_Test7=0b0 _one_Test7=0b1 _two_Test7=0b10
@test _zero_Test7.val === 0x00
@test _one_Test7.val === 0x01
@test _two_Test7.val === 0x02
@test typeof(_zero_Test7.val) <: UInt8

@test_throws ArgumentError eval(:(@enum Test8 _zero="zero"))
@test_throws ArgumentError eval(:(@enum Test9 _zero='0'))

@enum Test8 _zero_Test8=zero(Int64)
@test typeof(_zero_Test8.val) <: Int64
@test _zero_Test8.val === Int64(0)

@enum Test9 _zero_Test9 _one_Test9=0x01 _two_Test9
@test typeof(_zero_Test9.val) <: Int
@test _zero_Test9.val === 0
@test typeof(_one_Test9.val) <: Int
@test _one_Test9.val === 1
@test typeof(_two_Test9.val) <: Int
@test _two_Test9.val === 2

@enum Test10 _zero_Test10=0x00 _one_Test10 _two_Test10
@test typeof(_zero_Test10.val) <: UInt8
@test _zero_Test10.val === 0x00
@test typeof(_one_Test10.val) <: UInt8
@test _one_Test10.val === 0x01
@test typeof(_two_Test10.val) <: UInt8
@test _two_Test10.val === 0x02
=#

# test macro handles keyword arguments
@enum(Test11, _zero_Test11=2,
              _one_Test11,
              _two_Test11=5,
              _three_Test11)

@test Int(_zero_Test11) == 2
@test Int(_one_Test11) == 3
@test Int(_two_Test11) == 5
@test Int(_three_Test11) == 6

# don't allow enum value to overflow
@test_throws ArgumentError @eval(@enum EnumOvf x=typemax(Int32) y)

# test for unique Enum values
@test_throws ArgumentError eval(:(@enum(Test14, _zero_Test14, _one_Test14, _two_Test14=0)))

@test repr(apple) == "apple::$(string(Fruit)) = 0"
@test string(apple) == "apple"

@test reprmime("text/plain", Fruit) == "Enum $(string(Fruit)):\napple = 0\norange = 1\nkiwi = 2"
@test reprmime("text/plain", orange) == "orange::$(curmod_prefix)Fruit = 1"

@enum LogLevel DEBUG INFO WARN ERROR CRITICAL
@test DEBUG < CRITICAL

# serialization
let b = IOBuffer()
    serialize(b, apple)
    seekstart(b)
    @test deserialize(b) === apple
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
"""
    function compute_IoM!(obj)

Computes the integrals of motion of `obj.ψ` and saves them in respective fields of `obj`.
`obj` can be a `::Sim` or `::Calc` object. The fields are:

```
obj.KE  # Array containing the kinetic energy K(x)
obj.PE  # Array containing the potential energy V(x)
obj.E   # Array containing the energy H(x)
obj.N   # Array containing the norm N(x)
obj.P   # Array containing the momentum P(x)
obj.dE  # Array containing the energy error δE(x)
obj.dN  # Array containing the norm error δN(x)
obj.dP  # Array containing the momentum error δP(x)
```

Result can be plotted using `plot(obj, :IoM)`
"""
function compute_IoM!(obj)
    @info "Computing integrals of motions"
    ψ² = abs2.(obj.ψ)
    ψ̃² = abs2.(obj.ψ̃)
    # The norm is preserved extremely well so if the W.F. is normalized we do not worry about dividing the IoM by the norm to improve performance
    #if normalize
        #obj.N .= sum(ψ², dims=1)[:]./obj.box.Nₜ
        #obj.PE .= -0.5*sum(ψ².^2,dims=1)[:]./(obj.N*obj.box.Nₜ)
        #obj.KE .= 0.5*sum((obj.box.ω.^2 .* ψ̃²),dims=1)[:]./obj.N
        #obj.P .= 2*imag.(sum(im * obj.box.ω.* ψ̃²,dims=1)[:]./obj.N)
        #obj.P .= 2*sum(obj.box.ω.* ψ̃²,dims=1)[:]./obj.N
        #obj.E .= obj.KE .+ obj.PE
    #else
    # Why does the stuff from NumericalIntegration.jl not work?
    #obj.N .= integrate(obj.box.t, ψ², SimpsonEvenFast())./obj.T
    #obj.PE .= -0.5*integrate(obj.box.t, (ψ²).^2, SimpsonEvenFast())./obj.T
    obj.N .= sum(ψ², dims=1)[:]./obj.box.Nₜ
    obj.PE .= -0.5*sum(ψ².^2,dims=1)[:]./obj.box.Nₜ
    obj.KE .= 0.5*sum((obj.box.ω.^2 .* ψ̃²),dims=1)[:]
    #obj.P .= 2*imag.(sum(im * obj.box.ω.* ψ̃²,dims=1)[:])
    obj.P .= 2*sum(obj.box.ω.* ψ̃²,dims=1)[:]
    obj.E .= obj.KE .+ obj.PE
    #end

    @info "Integrals of motion computed."
end

"""
    function params(; kwargs...)

Computes parameters `λ, T, Ω` given either of `λ, T, Ω, a` and possibly `m` for dnoidal seed. Used only for `exp` and `dn` seeds.
"""
function params(;m=0.0, kwargs...)
    if length(kwargs) != 1
        throw(ArgumentError("You have either specified too few or too many parameters. You must specify one and only one of the following options: λ, Ω, T, a."))
    end
    param = Dict(kwargs)
    if m ==0.0
        if :a in keys(param)
            λ = im * sqrt(2 * param[:a])
            T = π/sqrt(1 - imag(λ)^2)
            Ω = 2π/T
            @info "Passed a = $(param[:a]), computed λ = $λ, T = $T and Ω = $Ω"
        elseif :λ  in keys(param)
            λ = param[:λ]
            T = π/sqrt(1 - imag(λ)^2)
            Ω = 2π/T
            @info "Passed λ=$λ, computed T = $T and Ω = $Ω"
        elseif :Ω in keys(param)
            λ = im * sqrt((1 - (param[:Ω] / 2)^2))
            T = π/sqrt(1 - imag(λ)^2)
            Ω = param[:Ω]
            @info "Passed Ω=$Ω, computed λ = $λ and T = $T"
        elseif :T in keys(param)
            λ = im * sqrt((1 - ((2*π/param[:T]) / 2)^2))
            T = param[:T]
            Ω = 2π/T
            @info "Passed T = $T, computed λ = $λ and Ω = $Ω"
        end
    elseif m > 0 && m <= 1
        if :a in keys(param)
            @error "Passing a not yet supported in m != 0 mode. Please pass λ"
        elseif :λ  in keys(param)
            λ = param[:λ]
            Ω = real(2*sqrt(1 + (λ - m/4/λ)^2))
            T = 2π/Ω
            @info "Passed λ=$λ, computed T = $T and Ω = $Ω"
        elseif :Ω in keys(param)
            @error "Passing Ω not yet supported in m != 0 mode. Please pass λ"
        elseif :T in keys(param)
            @error "Passing T not yet supported in m != 0 mode. Please pass λ"
        end
    elseif m > 1 || m < 0
        @error "Wrong range for m. m must be between 0 and 1"
    end

    return λ, T, Ω
end #compute_parameters

###########################################################################
# ψ₀
###########################################################################
"""
    function ψ₀_periodic(coeff::Array, box::Box, Ω; phase=0)

Computes an initial wavefunction for the `box::Box`, with fundamental frequency `Ω`
and coefficients ``A_{1\\ldots n}`` = `coeff` and an overall phase ``e^{i \\phi t}`` where ``\\phi`` = `phase`. i.e. ``\\psi_0`` is of the form:

``
\\psi(x=0, t) = e^{i \\phi t} (A_0 + 2 \\sum_{1}^{n} A_m \\cos(m \\Omega t))
``

where:

``
A_0 = \\sqrt{1 - 2 \\sum_{m=1}^n |A_m|^2}
``
"""
function ψ₀_periodic(coeff::Array, box::Box, Ω; phase=0)
    @info "Initializing periodic ψ₀"
    for (n, An) in enumerate(coeff)
        if abs(An) >= 1
            @error "The absolute value of the coefficient A($(n+1)) = $(An) is greater than 1. psi_0 needs to be normalizable."
        end #if
    end #for

    @info "Computing A₀ to preserve normalization."
    A0 = sqrt(1 - 2 * sum([abs(An) .^ 2 for An in coeff]))
    @info "Computed A₀ = $A0"
    #A0 = 1
    if phase != 0
        str = "ψ₀ = exp(i $phase t) ($A0 + "
    else
        str = "ψ₀ = $A0 + "
    end

    ψ₀ = A0 * ones(box.Nₜ)
    for (n, An) in enumerate(coeff)
        ψ₀ += 2 * An * cos.(n * Ω * box.t)
        str = string(str, "2 × $An × cos($n × $(Ω) t)")
    end #for
    # Multiply by the overall phase
    ψ₀ = exp.(im * phase * box.t) .* ψ₀

    if phase != 0
        str = string(str, ")")
    end

    @info str

    return ψ₀, A0
end #psi0_periodic

"""
    ψ₀_DT(λ, tₛ, xₛ, X₀, box; seed="exp", f = Dict(:α=> 0.0, :γ => 0.0, :δ=>0.0))

Computes an initial condition using the DT characterized by `λ, xₛ, tₛ` and `seed` at x = `X₀` for the extended equation characterized by a dictionary of parameters `f` (currently only `α`) is supported. `box::Box` is the simulation box.

See also: [`Calc`](@ref), [`Sim`](@ref)
"""
function ψ₀_DT(λ, tₛ, xₛ, X₀, box; seed="exp", f = Dict(:α=> 0.0, :γ => 0.0, :δ=>0.0))
    xᵣ = X₀=>X₀+1e-5
    T = abs(box.t[1]*2)
    Nₜ = box.Nₜ
    box_dt = Box(xᵣ, T, Nₓ=1, Nₜ = box.Nₜ)

    calc = Calc(λ, tₛ, xₛ, seed, box_dt, f = f)
    solve!(calc)

    return calc.ψ[:, 1]
end

"""
    λ_maximal(λ₁, N; m = 0)

Computed a maximal intensity set of `λ` of order `N` given `λ₁` and possibly `m` for dnoidal bakcground.

See also: [`λ_given_m`](@ref)
"""
function λ_maximal(λ₁, N; m = 0)
    ν₁ = imag(λ₁)
    #ν_min = sqrt(1 - 1/N^2)
    mp = 1 - m
    C_N = sqrt((N^2-1)*mp*(N^2 - mp))
    H_N = N^2*(mp + 1) - 2*mp
    ν_min = sqrt(2*C_N + H_N)/(2*N)
    if ν₁ <= ν_min
        throw(ArgumentError("λ = $λ₁ not big enough for N = $N, need at least λ = $ν_min im"))
    end
    n = (1:N)
    #λ = sqrt.(n.^2 .* (ν₁^2 - 1) .+ 1)*im
    G_n = (m^2 .* n.^2) .+ 8*(m-2).*(n.^2 .- 1)*ν₁.^2 .+ 16 .* n.^2  .* ν₁^4
    λ = sqrt.(G_n .+ sqrt.(G_n.^2 .- 64*m^2*ν₁^4))./(4*sqrt(2)*ν₁)*im
end

"""
    λ_given_m(m; q = 2)

Computed a `λ` that is matched to the dnoidal background given by `m` with an integer `q`. See paper for more details.

See also: [`λ_maximal`](@ref)
"""
function λ_given_m(m; q = 2)
    F = π/(2*q*Elliptic.K(m))
    λ = 0.5 * sqrt(2 - 2*F^2 - m + 2*sqrt((F^2-1)*(F^2-1+m)))*im
end

"""
    λ_given_f(f, ν)

Computes `λ = v + i ν` such that `v` results in breather to soliton conversion for the extended NLS characterized by the dictionary of parameters `f`.
"""
function λ_given_f(f, ν)
    α = f[:α]
    γ = f[:γ]
    δ = f[:δ]
    @info "Computing real part of eigenvalue for breather to soliton conversion"
    @info "Got α = $α, γ = $γ, δ = $δ"
    poly = reverse([64*δ, -24*γ, -8*(α+2*δ+8*δ*ν^2), γ*(4 + 8*ν^2)+1])
    r = roots(poly)

    if imag(r[1]) == 0
        @info "Computed v = $(real(r[1]))"
        return real(r[1]) + ν*im
    elseif imag(r[2]) == 0
        @info "Computed v = $(real(r[2]))"
        return real(r[2]) + ν*im
    elseif imag(r[3]) == 0
        @info "Computed v = $(real(r[3]))"
        return real(r[3]) + ν*im
    end
end

"""
    PHF(calc::Calc)

Computed the peak height formula for any `calc::Calc` objecti irrespective of seed.
"""
function PHF(calc::Calc)
    s = 2*sum(imag.(calc.λ))
    if calc.seed == "exp" || calc.seed == "dn"
        ψ₀₀ = 1
    elseif calc.seed == "cn"
        ψ₀₀ = sqrt(calc.m)
    elseif calc.seed == "0"
        ψ₀₀ = 0
    end
    peak = ψ₀₀ + s
end


text:
# This file was generated, do not modify it. # hide
m1 = machine(model1, X, y)
m2 = machine(model2, X, y)


text:
using GigaSOM
using Test

@test hello("Julia") == "Hello, Julia"


text:
"""
    Peters{N,TF,TV<:SVector{N,TF},TA<:SMatrix{N,N,TF}}

Two-dimensional aerodynamic model based on Peters' finite state model with `N` state
variables, inputs ``u, \\omega, \\dot{v}, \\dot{\\omega}`` and parameters ``a, b, a_0,
\\alpha_0, c_{d,0}, c_{m,0}``
"""
struct Peters{N,TF,TV<:SVector{N,TF},TA<:SMatrix{N,N,TF}}
    A::TA
    b::TV
    c::TV
end

"""
    Peters{N, TF=Float64}()

Initialize a model with type [`Peters`](@ref) with `N` state variables
"""
Peters()

Peters{N}() where N = Peters{N,Float64}()

function Peters{N,TF}() where {N,TF}

    A, b, c = peters_constants(N, TF)

    return Peters(A, b, c)
end

# --- Submodel Creation --- #

function Submodel(model::Peters{N}) where N

    Abar = model.A
    bbar = model.b
    cbar = model.c

    # residual function
    fresid = (dx, x, y, p, t) -> peters_residual(dx, x, y, p, t; Abar, cbar)

    # number of state, input, and parameters (use Val(N) to use inferrable dimensions)
    nx = Val(N)
    ny = Val(4)
    np = Val(6)

    # jacobian definitions
    ratejac = Invariant(Abar)
    statejac = Linear((dx, x, y, p, t) -> peters_state_jacobian(dx, x, y, p, t; cbar))
    inputjac = Nonlinear((dx, x, y, p, t) -> peters_input_jacobian(dx, x, y, p, t; Abar, cbar))
    paramjac = Nonlinear((dx, x, y, p, t) -> peters_parameter_jacobian(dx, x, y, p, t; Abar, cbar))
    tgrad = Zeros()

    # convenience functions for setting states, inputs, and parameters
    setstate = peters_set_states!
    setinput = peters_set_inputs!
    setparam = peters_set_parameters!

    # convenience functions for separating states, inputs, and parameters
    sepstate = peters_separate_states
    sepinput = peters_separate_inputs
    sepparam = peters_separate_parameters

    # model definition
    return Submodel{false}(fresid, nx, ny, np;
        ratejac = ratejac,
        statejac = statejac,
        inputjac = inputjac,
        paramjac = paramjac,
        tgrad = tgrad,
        setstate = setstate,
        setinput = setinput,
        setparam = setparam,
        sepstate = sepstate,
        sepinput = sepinput,
        sepparam = sepparam,
    )
end

# --- Internal Methods --- #

# function for defining model constants
function peters_constants(N, TF)

    b = zeros(TF, N)
    for n = 1:N-1
        b[n] = (-1)^(n-1)*factorial(big(N + n - 1))/factorial(big(N - n - 1))*
            1/factorial(big(n))^2
    end
    b[N] = (-1)^(N-1)

    c = zeros(TF, N)
    for n = 1:N
        c[n] = 2/n
    end

    d = zeros(TF, N)
    d[1] = 1/2

    D = zeros(TF, N, N)
    for m in 1:N-1
        n = m + 1
        D[n, m] = 1/(2*n)
    end
    for m in 2:N
        n = m - 1
        D[n, m] = -1/(2*n)
    end

    A = D + d*b' + c*d' + 1/2*c*b'

    return SMatrix{N,N}(A), SVector{N}(b), SVector{N}(c)
end

# residual function
function peters_residual(dx, x, y, p, t; Abar, cbar)

    dλ = similar_type(cbar, eltype(dx))(dx)
    λ = similar_type(cbar, eltype(x))(x)
    u, ω, vdot, ωdot = y
    a, b, a0, α0, cd0, cm0 = p

    return Abar*dλ - cbar*(vdot + u*ω + (b/2-a*b)*ωdot) + u/b*λ
end

# rate jacobian function
peters_rate_jacobian(; Abar) = Abar

# state jacobian function
function peters_state_jacobian(dx, x, y, p, t; cbar)

    u, ω, vdot, ωdot = y
    a, b, a0, α0, cd0, cm0 = p

    return u/b*Diagonal(ones(similar_type(cbar)))
end

# input jacobian function
function peters_input_jacobian(dx, x, y, p, t; Abar, cbar)

    λ = similar_type(cbar, eltype(x))(x)
    u, ω, vdot, ωdot = y
    a, b, a0, α0, cd0, cm0 = p

    dλ_u = -cbar*ω + λ/b
    dλ_ω = -u*cbar
    dλ_vdot = -cbar
    dλ_ωdot = -(b/2-a*b)*cbar

    return hcat(dλ_u, dλ_ω, dλ_vdot, dλ_ωdot)
end

# parameter jacobian function
function peters_parameter_jacobian(dx, x, y, p, t; Abar, cbar)

    λ = similar_type(cbar, eltype(x))(x)
    u, ω, vdot, ωdot = y
    a, b, a0, α0, cd0, cm0 = p

    dλ_a = b*ωdot*cbar
    dλ_b = -cbar*(1/2-a)*ωdot - u/b^2*λ
    dλ_a0 = zero(cbar)
    dλ_α0 = zero(cbar)
    dλ_cd0 = zero(cbar)
    dλ_cm0 = zero(cbar)

    return hcat(dλ_a, dλ_b, dλ_a0, dλ_α0, dλ_cd0, dλ_cm0)
end

# convenience function for defining this model's state vector
peters_set_states!(x; lambda) = x .= lambda

# convenience function for defining this model's input vector
function peters_set_inputs!(y; u, omega, vdot, omegadot)
    y[1] = u
    y[2] = omega
    y[3] = vdot
    y[4] = omegadot
    return y
end

# convenience function for defining this model's parameter vector
function peters_set_parameters!(p; a, b, a0, alpha0, cd0, cm0)
    p[1] = a
    p[2] = b
    p[3] = a0
    p[4] = alpha0
    p[5] = cd0
    p[6] = cm0
    return p
end

# convenience function for separating this model's state vector
peters_separate_states = (x) -> (lambda=x,)

# convenience function for separating this model's input vector
peters_separate_inputs = (y) -> (u=y[1], omega=y[2], vdot=y[3], omegadot=y[4])

# convenience function for separating this model's parameter vector
peters_separate_parameters = (p) -> (a=p[1], b=p[2], a0=p[3], alpha0=p[4], cd0=p[5], cm0=p[6])

# aerodynamic loads per unit span
function peters_loads(a, b, ρ, c, a0, α0, cd0, cm0, bbar, u, v, ω, vdot, ωdot, λ)
    # circulatory load factor
    tmp1 = a0*ρ*u*b
    # non-circulatory load factor
    tmp2 = pi*ρ*b^3
    # constant based on geometry
    d = b/2 - a*b
    # induced flow velocity
    λ0 = 1/2 * bbar'*λ
    # Velocity Magnitude (squared)
    V2 = u^2 + v^2
    # Mach Number (squared)
    M2 = V2/c^2
    # Prandtl-Glauert correction factor
    beta = sqrt(1 - min(0.99, M2))
    # normal force at the reference point
    N = tmp1*(v + d*ω - λ0 - u*α0) + tmp2*(vdot/b + u/b*ω - a*ωdot)
    # axial force at the reference point
    A = -a0*ρ*b*(v + d*ω - λ0 - u*α0)^2
    # moment at reference point
    M = -tmp2*(vdot/2 + u*ω + b*(1/8 - a/2)*ωdot) + 2*ρ*b^2*u^2*cm0 + (b/2 + a*b)*N
    # apply compressibility correction
    N = N / beta
    A = A / beta
    M = M / beta
    # add skin friction drag
    A += ρ*b*u^2*cd0
    # return loads
    return SVector(N, A, M)
end