text:
using Gtk.ShortNames, GtkReactive, Graphics, Colors

win = Window("Drawing")
c = canvas(UserUnit)       # create a canvas with user-specified coordinates
push!(win, c)

const lines = Signal([])   # the list of lines that we'll draw
const newline = Signal([]) # the in-progress line (will be added to list above)

const drawing = Signal(false)  # this will become true if we're actively dragging

# c.mouse.buttonpress is a `Reactive.Signal` that updates whenever the
# user clicks the mouse inside the canvas. The value of this signal is
# a MouseButton which contains position and other information.

# We're going to define a callback function that runs whenever the
# button is clicked. If we just wanted to print the value of the
# returned button object, we could just say
#     map(println, c.mouse.buttonpress)
# However, here our function is longer than `println`, so
# we're going to use Julia's do-block syntax to define the function:
sigstart = map(c.mouse.buttonpress) do btn
    # This is the beginning of the function body, operating on the argument `btn`
    if btn.button == 1 && btn.modifiers == 0 # is it the left button, and no shift/ctrl/alt keys pressed?
        push!(drawing, true)   # activate dragging
        push!(newline, [btn.position])  # initialize the line with the current position
    end
end

const dummybutton = MouseButton{UserUnit}()
# See the Reactive.jl documentation for `filterwhen`
sigextend = map(filterwhen(drawing, dummybutton, c.mouse.motion)) do btn
    # while dragging, extend `newline` with the most recent point
    push!(newline, push!(value(newline), btn.position))
end

sigend = map(c.mouse.buttonrelease) do btn
    if btn.button == 1
        push!(drawing, false)  # deactivate dragging
        # append our new line to the overall list
        push!(lines, push!(value(lines), value(newline)))
        # For the next click, make sure `newline` starts out empty
        push!(newline, [])
    end
end

# Because `draw` isn't a one-line function, we again use do-block syntax:
redraw = draw(c, lines, newline) do cnvs, lns, newl  # the function body takes 3 arguments
    fill!(cnvs, colorant"white")   # set the background to white
    set_coordinates(cnvs, BoundingBox(0, 1, 0, 1))  # set coordinates to 0..1 along each axis
    ctx = getgc(cnvs)   # gets the "graphics context" object (see Cairo/Gtk)
    for l in lns
        drawline(ctx, l, colorant"blue")  # draw old lines in blue
    end
    drawline(ctx, newl, colorant"red")    # draw new line in red
end

function drawline(ctx, l, color)
    isempty(l) && return
    p = first(l)
    move_to(ctx, p.x, p.y)
    set_source(ctx, color)
    for i = 2:length(l)
        p = l[i]
        line_to(ctx, p.x, p.y)
    end
    stroke(ctx)
end

showall(win)

max_stars_repo_name: ashwani-rathee/Chatbots

max_issues_repo_name: ashwani-rathee/Chatbots

max_forks_repo_name: ashwani-rathee/Chatbots


text:
# Original Java version Copyright 2005 Google Inc.
# Julia translation Copyright 2022 Ted Dunning
#
# 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.

md"""
An S1Interval represents a closed interval on a unit circle (also known as a 1-dimensional
sphere). It is capable of representing the empty interval (containing no points), the full
interval (containing all points), and zero-length intervals (containing a single point).

Points are represented by the angle they make with the positive x-axis in the range [-π, π].
An interval is represented by its lower and upper bounds (both inclusive, since the interval is
closed). The lower bound may be greater than the upper bound, in which case the interval is
"inverted" (i.e. it passes through the point (-1, 0)).

Note that the point (-1, 0) has two valid representations, π and -π. The normalized
representation of this point internally is π, so that endpoints of normal intervals are in the
range (-π, π]. However, we take advantage of the point -π to construct two special intervals:
the full() interval is [-π, π], and the Empty() interval is [π, -π].
"""
struct S1Interval
    lo::Float64
    hi::Float64

    S1Interval(low::Number, high::Number) = S1Interval(low, high, false)
    S1Interval(a::S1Interval) = S1Interval(a.lo, a.hi)
    function S1Interval(lowx::Number, highx::Number, checked::Bool)
        low, high = Float64(lowx), Float64(highx)
        lo, hi = low, high
        if !checked
            if low == -S2_M_PI && high != S2_M_PI
                lo = S2_M_PI
            end
            if high == -S2_M_PI && low != S2_M_PI
                hi = S2_M_PI
            end
        end
        new(lo, hi)
    end
end

md"""The empty interval"""
Base.empty(T::Type{S1Interval}) = S1Interval(S2_M_PI, -S2_M_PI, true)

md"""The full interval"""
full(T::Type{S1Interval}) = S1Interval(-S2_M_PI, S2_M_PI, true)

md"Construct an interval containing a single point"
function fromPoint(::Type{S1Interval}, p::Number)
    if p == -S2_M_PI
        p = S2_M_PI
    end

    return S1Interval(p, p)
end

md"""
Convenience method to construct the minimal interval containing the two given points. This is
equivalent to starting with an empty interval and calling addPoint() twice, but it is more
efficient.
"""
function fromPointPair(::Type{S1Interval}, p1::Number, p2::Number)::S1Interval
    if (p1 == -S2_M_PI)
        p1 = S2_M_PI
    end
    if (p2 == -S2_M_PI)
        p2 = S2_M_PI
    end
    if (positiveDistance(p1, p2) <= S2_M_PI)
        S1Interval(p1, p2)
    else
        S1Interval(p2, p1)
    end
end

md"""
An interval is valid if neither bound exceeds π in absolute value, and the value -π appears
only in the empty() and full() intervals.
"""
Base.isvalid(a::S1Interval) = abs(a.lo) <= S2_M_PI && abs(a.hi) <= S2_M_PI &&
    !(a.lo == -S2_M_PI && a.hi != S2_M_PI) && !(a.hi == -S2_M_PI && a.lo != S2_M_PI)

md"Returns true if the interval contains all points on the unit circle"
isfull(a::S1Interval) = a.hi - a.lo == 2 * S2_M_PI

md"Returns true if the interval is empty, i.e. it contains no points"
Base.isempty(a::S1Interval) = a.lo - a.hi == 2 * S2_M_PI

md"Returns true if lo() > hi(). (This is also true for empty intervals)"
isinverted(a::S1Interval) = a.lo > a.hi

md"""
Returns the midpoint of the interval. For full and empty intervals, the result is arbitrary
"""
function center(a::S1Interval)
    center = 0.5 * (a.lo + a.hi)
    if !isinverted(a)
        return center
    end
    # Return the center in the range (-π, π].
    return (center <= 0) ? (center + S2_M_PI) : (center - S2_M_PI)
end

md"Returns the length of the interval. The length of an empty interval is negative"
function Base.length(a::S1Interval)
    length = a.hi - a.lo
    if length >= 0
      return length
    end
    length += 2 * S2_M_PI
    # Empty intervals have a negative length.
    return (length > 0) ? length : -1
end

md"""
Return the complement of the interior of the interval. An interval and its complement have the
same boundary but do not share any interior values. The complement operator is not a bijection,
since the complement of a singleton interval (containing a single value) is the same as the
complement of an empty interval.
"""
function complement(a::S1Interval)
    if (a.lo == a.hi)
        return full(S1Interval)
    end
    # Handles empty and full.
    return S1Interval(a.hi, a.lo, true)
end

md"""
Return the midpoint of the complement of the interval. For full and empty intervals, the result
is arbitrary. For a singleton interval (containing a single point), the result is its antipodal
point on S1.
"""
function complementcenter(a::S1Interval)
    if a.lo != a.hi
        return center(complement(a))
    else
        # Singleton
        return (a.hi <= 0) ? (a.hi + S2_M_PI) : (a.hi - S2_M_PI)
    end
end

md"Returns true if the interval (which is closed) contains the point `p`"
function Base.contains(a::S1Interval, px::Number)
    p = Float64(px)
    # Works for empty, full, and singleton intervals.
    if p == -S2_M_PI
        p = S2_M_PI
    end
    if (isinverted(a))
        return (p >= a.lo || p <= a.hi) && !isempty(a)
    else
        return a.lo <= p <= a.hi
    end
end

md"Returns true if the interior of the interval contains the point `p`"
function interiorcontains(a::S1Interval, px::Number)
    p = Float64(px)
    # Works for empty, full, and singleton intervals.
    if p == -S2_M_PI
        p = S2_M_PI
    end

    if isinverted(a)
        return p > a.lo || p < a.hi
    else
        return (a.lo < p < a.hi) || isfull(a)
    end
end

md"""
Returns true if the interval `a` contains the interval `b`. Works for empty, full, and
singleton intervals.
"""
function Base.contains(a::S1Interval, b::S1Interval)
    # It might be helpful to compare the structure of these tests to
    # the simpler Contains(double) method above.
    if isinverted(a)
        if isinverted(b)
            return b.lo >= a.lo && b.hi <= a.hi
        end
        return (b.lo >= a.lo || b.hi <= a.hi) && !isempty(a)
    else
        if isinverted(b)
            return isfull(a) || isempty(b)
        end
    end
    return b.lo >= a.lo && b.hi <= a.hi
end

md"""
Returns true if the interior of interval `a` contains the entire interval `b`. Note that
`interiorcontains(x, x)` is true only when `x` is the empty or full interval, and
`interiorcontains(x, S1Interval(p,p))` is equivalent to `InteriorContains(x, p)`.
"""
function interiorcontains(a::S1Interval, b::S1Interval)
    if isinverted(a)
        if !isinverted(b)
            return b.lo > a.lo || b.hi < a.hi
        end
        return (b.lo > a.lo && b.hi < a.hi) || isempty(b)
    else
        if isinverted(b)
            return isfull(a) || isempty(b)
        end
        return (b.lo > a.lo && b.hi < a.hi) || isfull(a)
    end
end

md"""
Returns true if the two intervals contain any points in common. Note that the point +/-π has
two representations, so the intervals [-π,-3] and [2,π] intersect, for example.
"""
function intersects(a::S1Interval, b::S1Interval)
    if (isempty(a) || isempty(b))
        return false
    end
    if isinverted(a)
        # Every non-empty inverted interval contains π
        return isinverted(b) || b.lo <= a.hi || b.hi >= a.lo
    else
      if isinverted(b)
          return b.lo <= a.hi || b.hi >= a.lo
      end
      return b.lo <= a.hi && b.hi >= a.lo
    end
end

md"""
Returns true if the interior of this interval contains any point of the interval {@code y}
(including its boundary). Works for empty, full, and singleton intervals.
"""
function interiorintersects(a::S1Interval, b::S1Interval)
    if isempty(a) || isempty(b) || a.lo == a.hi
        return false
    end
    if isinverted(a)
        return isinverted(b) || b.lo < a.hi || b.hi > a.lo
    else
        if isinverted(b)
            return b.lo < a.hi || b.hi > a.lo
        end
      return (b.lo < a.hi && b.hi > a.lo) || isfull(a)
    end
end

md"""
Return the Hausdorff distance from `a` to interval `b`. For two S1Intervals `a` and `b`,
this distance is defined by `h(x, y) = max_{p in x} min_{q in y} d(p, q)`, where `d(.,.)`
is measured along `S1`.
"""
function getDirectedHausdorffDistance(a::S1Interval, b::S1Interval)
    if contains(b, a)
        # this accounts for empty `a`
        return 0.0
    end
    if isempty(b)
        # max distance in `S1`
        return S2_M_PI
    end

    bComplementCenter = complementcenter(b)
    if contains(a, bComplementCenter)
        return positiveDistance(b.hi, bComplementCenter)
    else
        # The Hausdorff distance is realized by either two hi endpoints or two
        # lo endpoints, whichever is farther apart.
        hiHi = contains(S1Interval(b.hi, bComplementCenter), a.hi) ? positiveDistance(b.hi, a.hi) : 0
        loLo = contains(S1Interval(bComplementCenter, b.lo), a.lo) ? positiveDistance(a.lo, b.lo) : 0
        return max(hiHi, loLo)
    end
end

md"""
Expands the interval by the minimum amount necessary so that it contains the point {@code p}
(an angle in the range [-π, π]).
"""
function addPoint(a::S1Interval, p::Number)
    if p == -S2_M_PI
        p = S2_M_PI
    end

    if contains(a, p)
        return a
    end

    if isempty(a)
        return fromPoint(S1Interval, p)
    else
        # Compute distance from p to each endpoint.
        dlo = positiveDistance(p, a.lo)
        dhi = positiveDistance(a.hi, p)
        if (dlo < dhi)
            return S1Interval(p, a.hi)
        else
            return S1Interval(a.lo, p)
        end
        # Adding a point can never turn a non-full interval into a full one.
    end
end

md"""
Returns the closest point in the interval to the point {@code p}. The interval must be
non-empty.
"""
function clampPoint(a::S1Interval, p::Number)
    if p == -S2_M_PI
        p = S2_M_PI
    end

    if contains(a, p)
        return p
    end

    # Compute distance from p to each endpoint.
    dlo = positiveDistance(p, a.lo)
    dhi = positiveDistance(a.hi, p)
    return (dlo < dhi) ? a.lo : a.hi
end

md"""
Returns a new interval that has been expanded on each side by the distance {@code margin}. If
"margin" is negative, then shrink the interval on each side by "margin" instead. The resulting
interval may be empty or full. Any expansion (positive or negative) of a full interval remains
full, and any expansion of an empty interval remains empty.
"""
function expanded(a::S1Interval, margin::Number)
    if margin >= 0
        if isempty(a)
            return a
        end
        # Check whether this interval will be full after expansion, allowing
        # for a 1-bit rounding error when computing each endpoint.
        if length(a) + 2 * margin + 2 * eps(Float64) >= 2 * S2_M_PI
            return full(S1Interval)
        end
    else
        if isfull(a)
            return a
        end
        # Check whether this interval will be empty after expansion, allowing
        # for a 1-bit rounding error when computing each endpoint.
        if length(a) + 2 * margin - 2 * eps(Float64) <= 0
            return empty(S1Interval)
        end
    end

    lo = rem(a.lo - margin, 2 * S2_M_PI, RoundNearest)
    hi = rem(a.hi + margin, 2 * S2_M_PI, RoundNearest)

    return S1Interval(lo, hi)
end

md"""
Returns the smallest interval that contains interval `a` and the interval `b`.
"""
function Base.union(a::S1Interval, b::S1Interval)
    lo = a.lo
    hi = a.hi

    # The `isfull(b)` case is handled correctly in all cases by the code
    # below, but can follow three separate code paths depending on whether
    # `a` is inverted, is non-inverted but contains π, or neither.
    if !isempty(b)
        if contains(a, b.lo)
            if contains(a, b.hi)
                # Either this interval contains y, or the union of the two
                # intervals is the full interval.
                if (!contains(a, b))
                    return full(S1Interval)
                end
            else
                hi = b.hi
            end
        elseif contains(a, b.hi)
            lo = b.lo
        elseif isempty(a) || contains(b, lo)
            # Interval `a` contains neither endpoint of `b`. This means that either `b`
            # contains all of `a`, or the two intervals are disjoint.
            lo = b.lo
            hi = b.hi
        else
            # Check which pair of endpoints are closer together.
            dlo = positiveDistance(b.hi, lo)
            dhi = positiveDistance(hi, b.lo)
            if dlo < dhi
                lo = b.lo
            else
                hi = b.hi
            end
        end
    end
    return S1Interval(lo, hi)
end

md"""
Returns the smallest interval that contains the intersection of interval `a` with `b`
Note that the region of intersection may consist of two disjoint intervals.
"""
function intersection(a::S1Interval, b::S1Interval)
    lo, hi = a.lo, a.hi
    # The `isfull(b)` case is handled correctly in all cases by the code below, but can follow three
    # separate code paths depending on whether this interval is inverted, is non-inverted but
    # contains π, or neither.

    if isempty(b)
        return empty(S1Interval)
    elseif contains(a, b.lo)
        if contains(a, b.hi)
            # Either `a` contains `b`, or the region of intersection consists of two disjoint
            # subintervals. In either case, we want to set the interval to the shorter of the two
            # original intervals.
            if length(b) < length(a)
                return S1Interval(b.lo, b.hi, true) # isfull(a) code path
            end
        else
            return S1Interval(b.lo, hi, true)
        end
    elseif contains(a, b.hi)
        return S1Interval(lo, b.hi, true)
    else
        # Interval `a` contains neither endpoint of `b`. This means that either `b`
        # contains all of this interval, or the two intervals are disjoint.
        if !contains(b, lo)
            return empty(S1Interval)
        end
    end
    return S1Interval(lo, hi)
end

md"""
Returns true if this interval can be transformed into the interval {@code y} by moving each
endpoint by at most "maxError" (and without the endpoints crossing, which would invert the
interval). Empty and full intervals are considered to start at an arbitrary point on the unit
circle, thus any interval with (length <= 2*maxError) matches the empty interval, and any
interval with (length >= 2*π - 2*maxError) matches the full interval.
"""
function Base.isapprox(a::S1Interval, b::S1Interval; atol=1e-15)
    # Full and empty intervals require special cases because the "endpoints"
    # are considered to be positioned arbitrarily.
    if isempty(a)
        return length(b) <= 2 * atol
    end
    if isempty(b)
        return length(a) <= 2 * atol
    end
    if isfull(a)
        return length(b) >= 2 * (S2_M_PI - atol)
    end
    if isfull(b)
        return length(a) >= 2 * (S2_M_PI - atol)
    end

    # The purpose of the last test below is to verify that moving the endpoints
    # does not invert the interval, e.g. [-1e20, 1e20] vs. [1e20, -1e20].
    return abs(rem(b.lo - a.lo, 2 * S2_M_PI, RoundNearest)) <= atol &&
        abs(rem(b.hi - a.hi, 2 * S2_M_PI, RoundNearest)) <= atol &&
        abs(length(a) - length(b)) <= 2 * atol
end

function isequals(a::S1Interval, b::S1Interval)
    return a.lo == b.lo && a.hi == b.hi
end

md"""
Computes the distance from {@code a} to {@code b} in the range [0, 2*π). This is equivalent to
{@code drem(b - a - π, 2 * π) + π}, except that it is more numerically stable
(it does not lose precision for very small positive distances).
"""
function positiveDistance(a::Number, b::Number)
    d = Float64(b - a);
    if d >= 0
        return d
    else
        # We want to ensure that if b == π and a == (-π + eps),
        # the return result is approximately 2π and not zero.
        return (b + S2_M_PI) - (a - S2_M_PI)
    end
end

ext: jl


text:
module BatteriesBattery

using ..Ahorn, Maple

@mapdef Entity "batteries/battery" Battery(x::Integer, y::Integer, maxCharge::Integer=500, initalCharge::Integer=500, dischargeRate::Integer=80, oneUse::Bool=false, onlyFits::Integer=-1, ignoreBarriers::Bool=false)

const placements = Ahorn.PlacementDict(
    "Battery" => Ahorn.EntityPlacement(
        Battery,
        "point"
    ),
    "Battery (Permanent)" => Ahorn.EntityPlacement(
        Battery,
        "point",
        Dict{String, Any}(
            "initalCharge" => 500,
            "dischargeRate" => 0,
            "oneUse" => true,
        )
    )
)

sprite = "batteries/battery/full0"

function Ahorn.selection(entity::Battery)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y - 8)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Battery, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, -8)

end

max_stars_repo_name: PenguinOwl/Batteries

max_issues_repo_name: PenguinOwl/Batteries

max_forks_repo_name: PenguinOwl/Batteries


text:
using BayesNets
using Discretizers
using JLD

basepath = "./results/Test1/"
paths = [string(basepath,i) for i in 0:4]

logname = "test_log_BN.txt"
logfile = "$basepath$logname"

scorename = "test_BN_p_dists.txt"
#scorefile = "$basepath$i/$scorename"
moveCats = [1,2,3]
numNextMoves = length(moveCats)
useLaneType = true
if contains(basepath,".1")
  useLaneType = false
end
laneTypeIndex = 3
laneTypeEncodedLen = length([0,0,0,0])
numYawBins = 20

max_parents = 7

#truncate the log file
f = open(logfile, "w")
close(f)

function Distributions.fit(::Type{DiscreteBayesNet}, data::DataFrame, params::GreedyHillClimbing, indexOfOut::Int;
    ncategories::Vector{Int} = map!(i->infer_number_of_instantiations(data[i]), Array(Int, ncol(data)), 1:ncol(data)),
    )

    n = ncol(data)
    parent_list = map!(i->Int[], Array(Vector{Int}, n), 1:n)
    datamat = convert(Matrix{Int}, data)'
    score_components = bayesian_score_components(parent_list, ncategories, datamat, params.prior, params.cache)

    while true
        best_diff = 0.0
        best_parent_list = parent_list
        i = indexOfOut
        # 1) add an edge (j->i)
        if length(parent_list[i]) < params.max_n_parents
            for j in deleteat!(collect(1:n), parent_list[i])
                    if adding_edge_preserves_acyclicity(parent_list, j, i)
                        new_parents = sort!(push!(copy(parent_list[i]), j))
                        new_component_score = bayesian_score_component(i, new_parents, ncategories, datamat, params.prior, params.cache)
                    if new_component_score - score_components[i] > best_diff
                        best_diff = new_component_score - score_components[i]
                        best_parent_list = deepcopy(parent_list)
                        best_parent_list[i] = new_parents
                    end
                end
             end
        end

        # 2) remove an edge
        for (idx, j) in enumerate(parent_list[i])

            new_parents = deleteat!(copy(parent_list[i]), idx)
            new_component_score = bayesian_score_component(i, new_parents, ncategories, datamat, params.prior, params.cache)
            if new_component_score - score_components[i] > best_diff
                best_diff = new_component_score - score_components[i]
                best_parent_list = deepcopy(parent_list)
                best_parent_list[i] = new_parents
            end

            # 3) flip an edge
            new_parent_list = deepcopy(parent_list) # TODO: make this more efficient
            deleteat!(new_parent_list[i], idx)

            if adding_edge_preserves_acyclicity(new_parent_list, i, j)
                sort!(push!(new_parent_list[j], i))
                new_diff = bayesian_score_component(i, new_parent_list[i], ncategories, datamat, params.prior, params.cache) - score_components[i]
                new_diff += bayesian_score_component(j, new_parent_list[j], ncategories, datamat, params.prior, params.cache) - score_components[j]
                if new_diff > best_diff
                    best_diff = new_diff
                    best_parent_list = new_parent_list
                end
            end
        end

        if best_diff > 0.0
            parent_list = best_parent_list
            score_components = bayesian_score_components(parent_list, ncategories, datamat, params.prior, params.cache)
        else
            break
        end
    end
    # construct the BayesNet
    cpds = Array(DiscreteCPD, n)
    varnames = names(data)
    for j in 1:n
        name = varnames[j]
        parents = varnames[parent_list[j]]
        cpds[j] = fit(DiscreteCPD, data, name, parents, params.prior,
                      parental_ncategories=ncategories[parent_list[j]],
                      target_ncategories=ncategories[j])
    end
    BayesNet(cpds)
end

open(logfile, "a") do logf
  write(logf,"Done making discretizers\n")
  write(logf, string(now(),"\n"))
end
allFeatures = Dict()
allTargets = Dict()
allFids = Dict()
nFids = 0
println("Loading features/targets for BN")
open(logfile, "a") do logf
  write(logf,"Loading features/targets for BN\n")
  write(logf, string(now(),"\n"))
end
numFeatures = 0
for path in paths
    i = path[end]
    f = open("$path/featureSet")
    allFeatures[i] = readdlm(f)
    allFeatures[i] = allFeatures[i][1:500,:]
    close(f)

    f = open("$path/targetSet")
    allTargets[i] = readdlm(f)
    allTargets[i] = allTargets[i][1:500,:]
    close(f)

    if numFeatures == 0
        numFeatures = length(allFeatures[i][1,:])
    end
    f = open("$path/Fids")
    allFids[i] = readdlm(f)
    #allFids[i] = allFids[i][1:1000]
    close(f)
    nFids += length(allFids[i])
end
moveIndex = numFeatures+1

println("Done loading features/targets for BN")
println("NumFeatures == ", numFeatures)
open(logfile, "a") do logf
  write(logf,"Done loading features/targets for BN\n")
  write(logf, string(now(),"\n"))
end

typeToEdges = Dict(
    "Vy"=>[0.0,0.5,20,40,60],
    "Vx"=>[-1000,-20,-5,-0.5,0.5,5,20,1000],
    "A"=>[-5000,-20,-5,-0.5,0.5,5,20,5000],
    "dist" => [0.0, 5, 25, 100, 500, 10000],
    "hdwy" => [0.0, 5, 25, 100, 250, 10000],)
idToType = Dict(
    4=>"Vy",
    5=>"A",
    6=> "Vx",
    7=> "A",
    8=> "yaw",
    9=> "hdwy",
    10=> "dist",)
#above dict was made when lanetype was integer
indexToID = zeros(Int64, (1, numFeatures+1+laneTypeEncodedLen))   #+1 for when no lane type, shift needs buffer
coreLen = 9 + laneTypeEncodedLen
histLen = 9
numHist = 4
totHist = histLen * numHist
traffLen = 8
for index in 1:numFeatures
    ID = index
    if index <= coreLen
      ID = index
      if useLaneType
        if index >= 4
          ID = index - laneTypeEncodedLen + 1
        end
      end
    elseif index <= coreLen + totHist
      ID = ((index - coreLen) % histLen)
      if ID == 0
        ID = histLen
      end
      ID += 1
    else
      ID = ((index - (coreLen + totHist)) % traffLen)
      if ID == 0
        ID = traffLen
      elseif ID < 3
        ID = traffLen + 1
      end
      ID += 1
    end
    indexToID[index] = ID
end
println("Done making discretizer dictionaries")

#big for loop
for i in keys(allFeatures)
  #truncate score file
  f = open("$basepath$i/$scorename", "w")
  close(f)
  println("Starting feature/target formatting for CV #: $i")
  open(logfile, "a") do logf
    write(logf,"Starting feature/target formatting for CV #: $i \n")
    write(logf, string(now(),"\n"))
  end
  testlines = allFeatures[i]
  testactuals = allTargets[i]

  trainlines = zeros(Float64, (1,numFeatures))
  trainactuals = zeros(Float64, (1,numFeatures))
  first = true
  for j in keys(allFeatures)
    if j != i
      if first == true
        trainlines = allFeatures[j]
        trainactuals = allTargets[j]
        first = false
      else
        trainlines = vcat(trainlines, allFeatures[j])
        trainactuals = vcat(trainactuals, allTargets[j])
      end
    end
  end
  println("Done feature/target formatting for CV #: $i")
  println("Starting to discretize for CV #: $i")
  open(logfile, "a") do logf
    write(logf,"Done feature/target formatting for CV #: $i \n")
    write(logf,"Starting to discretize for CV #: $i \n")
    write(logf, string(now(),"\n"))
  end
  #make discretizers with all lines, outlier problems
  allLines = vcat(trainlines, testlines)
  allactuals = vcat(trainactuals, testactuals)
  discAllLines = zeros(Int64, (length(trainlines[:,1])+length(testlines[:,1]), numFeatures+1))
  discTrainlines = zeros(Int64, (length(trainlines[:,1]), numFeatures+1))
  discTestlines = zeros(Int64, (length(testlines[:,1]), numFeatures+1))
  #DiscreteAlready = [1,2,3,11,12,20,21,29,30]
  for k in 1:numFeatures
    id = indexToID[k]
    if haskey(idToType, id) #if not, already discretized
      F_type = idToType[id]
      if haskey(typeToEdges, F_type)
        Disc = LinearDiscretizer(typeToEdges[F_type])
      else
        Disc = LinearDiscretizer(binedges(DiscretizeUniformWidth(numYawBins),allLines[:,k]))
      end
      DiscretizedAll = encode(Disc, allLines[:,k])
      DiscretizedTrain = encode(Disc, trainlines[:,k])
      DiscretizedTest = encode(Disc, testlines[:,k])
    else
      DiscretizedAll = allLines[:,k]
      DiscretizedTrain = trainlines[:,k]
      DiscretizedTest = testlines[:,k]
    end
    discTrainlines[:,k] = encode(CategoricalDiscretizer(DiscretizedAll), DiscretizedTrain)
    discTestlines[:,k] = encode(CategoricalDiscretizer(DiscretizedAll), DiscretizedTest)
    discAllLines[:,k] = encode(CategoricalDiscretizer(DiscretizedAll), DiscretizedAll)
  end
  moveDiscretizer = CategoricalDiscretizer(trainactuals)
  discTrainlines[:,moveIndex] = encode(moveDiscretizer,trainactuals)
  discTestlines[:,moveIndex] = encode(moveDiscretizer,testactuals)
  discAllLines[:,moveIndex] = encode(moveDiscretizer,allactuals)
  println("Done discretizing for CV #: $i")
  open(logfile, "a") do logf
    write(logf,"Done discretizing for CV #: $i \n")
    write(logf, string(now(),"\n"))
  end
  alldata = convert(DataFrame, discAllLines) #used to find num bins
  traindata = convert(DataFrame, discTrainlines)
  testdata = convert(DataFrame, discTestlines)

  println("Starting to fit BN for CV #: $i")
  open(logfile, "a") do logf
    write(logf,"Starting to fit BN for CV #: $i \n")
    write(logf, string(now(),"\n"))
  end
  max_parents = min(numFeatures, max_parents)
  println("Max parents: ", max_parents)
  params = GreedyHillClimbing(ScoreComponentCache(traindata), max_n_parents=max_parents, prior=UniformPrior())
  num_bins_all = map!(i->infer_number_of_instantiations(alldata[i]), Array(Int, ncol(alldata)), 1:ncol(alldata))
  bnDis = @time fit(DiscreteBayesNet, traindata, params, moveIndex; ncategories=num_bins_all)
  println("Done fitting BN for CV #: $i")
  open(logfile, "a") do logf
    write(logf,"Done fitting BN for CV #: $i \n")
    write(logf, string(now(),"\n"))
    println("Starting to score BN for CV #: $i")
    write(logf,"Starting to score BN for CV #: $i\n")
    write(logf, string(now(),"\n"))
  end
  save(string("$basepath$i","/BN_model.jld"), "bnDis", bnDis)
  score = 0
  numNaN = 0
  nrows = length(testlines[:,1])
  moveCats = unique(collect(testdata[:,moveIndex]))
  numNextMoves = length(moveCats)
  for index in 1:nrows
      actualNextMove = convert(Int64,decode(moveDiscretizer,testdata[index,moveIndex]))
      p_dists = zeros(Float64, (1, numNextMoves))
      for move in moveCats
        featureline = testdata[index,:]
        featureline[moveIndex] = encode(moveDiscretizer, move)
        p_dists[move] = pdf(bnDis, featureline)
      end
      for x in p_dists
        if !(x > 1) && !(x < 1) && !(x == 1)
          p_dists = [1.0/numNextMoves for i in 1:numNextMoves]
          numNaN += 1
          break
        end
      end
      p_dists = p_dists ./ (sum(p_dists) + 0.0000000000000001)
      p_right = p_dists[actualNextMove]
      score += (1.0  - p_right)
      scoref = open("$basepath$i/$scorename", "a")
      writedlm(scoref, reshape(p_dists, (1,numNextMoves)))
      close(scoref)
  end
  println("Score:", score)
  println("Num Predicitons:", nrows)
  println("Num NaN:", numNaN)
  scoref = open("$basepath$i/$scorename", "a")
  #writedlm(scoref, score)
  #writedlm(scoref, numNaN)
  close(scoref)
  println("Done scoring BN for CV #: $i")
  open(logfile, "a") do logf
    write(logf,"Done scoring BN for CV #: $i \n")
    write(logf, string(now(),"\n"))
  end
end #end of for i in allfeatures keys
#end  #for log file

max_stars_repo_name: zhangzhangzf/IntentionPrediction

max_issues_repo_name: zhangzhangzf/IntentionPrediction

max_forks_repo_name: zhangzhangzf/IntentionPrediction


text:
include("../required.jl")

# Calculate nestedness
ρ(B)
nodf(B)

max_stars_repo_name: kiristern/reseaux

max_issues_repo_name: kiristern/reseaux

max_forks_repo_name: kiristern/reseaux


text:
#
#  File:    qcqo1.jl
#
#  Purpose: Demonstrates how to solve small quadratic and quadratically
#           constrained optimization problem using the MOSEK Python API.
##

using Mosek
using Printf
# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:

# Set up and input bounds and linear coefficients
bkc   = [ MSK_BK_LO ]
blc   = [ 1.0 ]
buc   = [ Inf ]

bkx   = [ MSK_BK_LO
          MSK_BK_LO
          MSK_BK_LO ]
blx   = [ 0.0,  0.0, 0.0 ]
bux   = [ Inf,  Inf, Inf ]

c     = [ 0.0, -1.0, 0.0 ]

asub  = [ 1 ,2, 3 ]
aval  = [ 1.0, 1.0, 1.0 ]

numvar = length(bkx)
numcon = length(bkc)

# Create a task
maketask() do task
    # Append 'numcon' empty constraints.
    # The constraints will initially have no bounds.
    appendcons(task,numcon)

    #Append 'numvar' variables.
    # The variables will initially be fixed at zero (x=0).
    appendvars(task,numvar)

    #Optionally add a constant term to the objective.
    putcfix(task,0.0)
    # Set the linear term c_j in the objective.
    putclist(task,[1:numvar;],c)

    # Set the bounds on variable j
    # blx[j] <= x_j <= bux[j]
    putvarboundslice(task,1,numvar+1,bkx,blx,bux)
    # Input column j of A
    putarow(task,1,asub,aval)

    putconbound(task,1,bkc[1],blc[1],buc[1])

    # Set up and input quadratic objective

    qsubi = [ 1,   2,    3,   3   ]
    qsubj = [ 1,   2,    1,   3   ]
    qval  = [ 2.0, 0.2, -1.0, 2.0 ]

    putqobj(task,qsubi,qsubj,qval)

    # The lower triangular part of the Q^0
    # matrix in the first constraint is specified.
    # This corresponds to adding the term
    # - x0^2 - x1^2 - 0.1 x2^2 + 0.2 x0 x2

    qsubi = [  1,    2,    3,   3   ]
    qsubj = [  1,    2,    3,   1   ]
    qval  = [ -2.0, -2.0, -0.2, 0.2 ]

    # put Q^0 in constraint with index 0.

    putqconk(task,1, qsubi,qsubj, qval)

    # Input the objective sense (minimize/maximize)
    putobjsense(task,MSK_OBJECTIVE_SENSE_MINIMIZE)

    # Optimize the task
    optimize(task,"mosek://solve.mosek.com:30080")
    # Print a summary containing information
    # about the solution for debugging purposes
    solutionsummary(task,MSK_STREAM_MSG)
    prosta = getprosta(task,MSK_SOL_ITR)
    solsta = getsolsta(task,MSK_SOL_ITR)

    if solsta == MSK_SOL_STA_OPTIMAL
        # Output a solution
        xx = getxx(task,MSK_SOL_ITR)
        @printf("Optimal solution: %s\n", xx')
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        println("Primal or dual infeasibility.\n")
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        println("Primal or dual infeasibility.\n")
    elseif solsta == MSK_SOL_STA_UNKNOWN
        println("Unknown solution status")
    else
        println("Other solution status")
    end
end


text:
# SignalExtensions.jl

module SignalExtensions

using OffsetArrays

import Base: eltype, getindex, setindex!, eachindex, collect, nzrange

import Base: transpose, ctranspose, conj, reverse

import Base: *

# A collection of useful functions
include("util.jl")
# Definition of abstract Sequence type
include("sequences.jl")
# Definition of z-transform and Fourier transform
include("transforms.jl")
# The various types of extensions
include("extensions.jl")
# Extension types that store a vector
include("extension_sequences.jl")
# Lazy representations of sequence manipulations
include("lazy.jl")
# A number of special sequences
include("special_sequences.jl")

## We list all functions exported by the package below.
## This list is exhaustive and may serve as a reference for the functionality.

# From util.jl
export reverse_indices

# From sequences.jl:
# - types
export Sequence
# - iterator methods
export firstindex, lastindex, each_nonzero_index
# - properties
export iscompact

# From transforms.jl
# - types
export ZTransform, FourierTransform
# - methods
export ztransform, fouriertransform

# From extensions.jl
# - types
export PeriodicExtension, ZeroPadding, ConstantPadding, SymmetricExtension
# - methods
export element, periodic_index

# From extension_sequences.jl
# - types
export ExtensionSequence, PeriodicSequence, CompactSequence, ZeroPaddingSequence,
    ConstantPaddingSequence, SymmetricSequence
# - convenience constructors
export symmetric_extension_wholepoint_even, symmetric_extension_wholepoint_odd,
    symmetric_extension_halfpoint_even, symmetric_extension_halfpoint_odd
# - other methods
export subvector, first_subindex, last_subindex, sublength
export left_parity, right_parity, left_symmetry, right_symmetry

# From lazy.jl
# - types
export ModulatedSequence, DownsampledSequence, UpsampledSequence, ShiftedSequence,
    ReversedSequence, Convolution
# - methods
export downsample, upsample, convolve, shift, modulate, reverse
export sequence_shift, samplefactor

# From special_sequences.jl
export ZeroSequence, DiracSequence

end # module


text:
const deflayout = Plots.EmptyLayout(height=:auto, width=:auto)

"""
    plot_all(data, plottables, [label])

Plot all data for selected plottables.
"""
function plot_all end

plot_all(data, plottables) =  plot_all(data, plottables, "")
plot_all(data, plottables, label) = begin
    plots = []
    plot_all!(plots, data, plottables, label)
    num_plots = length(plots)
    if num_plots > 0
        row_layout = Plots.GridLayout(num_plots, 1)
        return plot(plots..., layout=row_layout)
    else
        return nothing
    end
end

"""
    plot_all!(plots, data, plottables, label)

Recursive plotting methods for (almost) all types.
"""
function plot_all! end

plot_all!(plots, data::Associative, plottables::Associative, label...) =
    for key in keys(plottables)
        plot_all!(plots, data[key], plottables[key], key)
    end
plot_all!(plots, data::Any, plottables::Associative, label...) = begin
    if haskey(plottables, :__data)
        index
        for key in keys(plottables)
            p = plot()
            data_plot!(p, getfield(data, key), plottables[key], key)
            push!(plots, p)
        end
    else
        for key in keys(plottables)
            plot_all!(plots, getfield(data, key), plottables[key], key)
        end
    end
end
plot_all!(plots, data::Tuple, plottables::Tuple, label...) =
    for i in 1:length(data)
        plot_all!(plots, data[i], plottables[i])
    end
plot_all!(plots, data::AxisArray{T,2}, plottables::AbstractArray, label...) where T = begin
    p = plot()
    labels = axes(data)[1]
    for i = 1:length(data[:, 1])
        plotit!(p, data, data[i, :], labels[i])
    end
    push_if_plotted!(plots, p)
end
plot_all!(plots, data::AxisArray{T,3}, plottables::AbstractArray, label...) where T = begin
    p = plot()
    labels1 = axes(data)[1]
    labels2 = axes(data)[2]
    rows = length(data[:, 1, 1])
    cols = length(data[1, :, 1])
    for i = 1:rows
        for j = 1:cols
            if plottables[i, j]
                l = string(labels1[i], ",", labels2[j])
                plotit!(p, data, data[i, j, :], l)
            end
        end
    end
    push_if_plotted!(plots, p)
end
@require DataFrames begin
    plot_all!(plots, data::DataFrame, plottables, label...) = begin
        p = plot()
        for (i, n) in enumerate(names(data))
            if plottables[i]
                plotit!(p, data, data[n], n)
            end
        end
        push_if_plotted!(plots, p)
    end
end
plot_all!(plots, data, plottables::Bool, label...) = begin
    if !plottables return end

    p = plot()
    if length(label) > 0
        l = label[1]
    else
        l = ""
    end
    plotit!(p, data, data, l)
    push_if_plotted!(plots, p)
end

data_plot!(p, data, plottables::Associative, keys) = begin
    for key in keys(plottables)
        data_plot!(p, data, plottables[key], (keys..., key))
    end
end
data_plot!(p, data, plottables::Bool, keys) = begin
    if plottables
        plotit!(p, data, [recursive_getfield(data, keys) for i in length(data)], keys[end])
    end
end

recursive_getfield(data, keys) =
    recursive_getfield(getfield(data, keys[1]), Base.tail(keys))
recursive_getfield(data, keys::Tuple{}) = data

"""
Return the intersect of the array range and the default plotting range.
"""
inbounds_range(data) = intersect(1:length(data), defaults[:plot_range])

"""
If it isn't empty push a plot to an array of plots,
"""
push_if_plotted!(plots, p) = if p.n > 0 push!(plots, p) end


text:
"""
    rvec(n, ρ::Matrix, margins::Vector{<:UnivariateDistribution})

Generate samples for a list of marginal distributions and a correaltion structure.

# Examples
```julia-repl
julia> using Distributions

julia> import LinearAlgebra: diagind

julia> margins = [Normal(3, 1), LogNormal(3, 1), Exponential(3)]

julia> R = fill(0.5, 3, 3); R[diagind(R)] .= 1.0;

julia> rvec(10, R, margins)
#>10×3 Matrix{Float64}:
 3.71109  82.0696   8.14332
 1.23022  38.8599   2.64595
 2.07222   3.76843  1.13465
 2.82434  11.6953   0.891066
 2.37599  10.0552   1.80555
 2.85431  25.4935   3.40865
 3.85253  21.2241   3.67532
 3.70605  59.2439   2.02502
 2.34066   1.89257  0.619948
 3.83507  16.787    0.66837
```
"""
function rvec end
for T in (Float64, Float32, Float16)
    @eval function rvec(n::Int, ρ::Matrix{$T}, margins::Vector{<:UD})
        d = length(margins)
        r,s = size(ρ)

        !(r == s == d) && throw(DimensionMismatch("The number of margins must match the size of the correlation matrix."))
        !iscorrelation(ρ) && throw(ValidCorrelationError())

        Z = SharedMatrix{$T}(_rmvn(n, ρ))
        @inbounds @threads for i in 1:d
            Z[:,i] = normal_to_margin(margins[i], Z[:,i])
        end
        sdata(Z)
    end
    @eval rvec(n::Real, ρ::Matrix{$T}, margins::Vector{<:UD}) = rvec(Int(n), ρ, margins)
end


text:
function suspend(cur_player::AbstractPlayer)
  cur_game = cur_player.game

  cur_game.in_focus = false

  cur_player.clock.push = Nullable{UUIDs.UUID}()
  cur_player.clock.drop = Nullable{UUIDs.UUID}()

  tetris_js(
    cur_player.game.scope,
    JSString("""
      \$(".js-tetris-container").addClass("cs-has-splash");
      \$(".js-splash-screen").removeClass("cs-disappear");

      \$(".js-pause-button").closest("a").removeClass("cs-active");
    """)
  )

  return
end


text:
function f=TBaroWave(x,y,eta,Param)
y=abs(y);
PhiPrim=Param.u0/2*((Param.f0-Param.beta0*Param.y0)...
  *(y-Param.Ly/2-Param.Ly/(2*pi)*sin(2*pi*y/Param.Ly))...
  +Param.beta0/2*(y*y-Param.Ly*y/pi*sin(2*pi*y/Param.Ly)...
  -Param.Ly*Param.Ly/(pi*pi*2)*cos(2*pi/Param.Ly)-...
  -Param.Ly*Param.Ly/3-Param.Ly*Param.Ly/(pi*pi*2)));
TAvg=Param.T0*eta^(Param.Rd*Param.Lapse/Param.Grav);
f=TAvg+PhiPrim/Param.Rd*(2/Param.b^2*(log(eta))^2-1)...
  *exp(-(log(eta)/Param.b)^2);
end


text:
abstract type AbstractCompactLattice{T,N_ARR,N_CDT} <: AbstractLattice{T,N_ARR,N_CDT} end

@doc """
A lattice of points, each an `N_ARR` coordinate, evenly distributed.

Can be constructed using default AbstractLattice constructor
"""
struct CompactLattice{T,N_ARR} <: AbstractCompactLattice{T,N_ARR,N_ARR}
    arr::Array{NTuple{N_ARR,T},N_ARR}
end
difference(lattice::AbstractCompactLattice, edge) = abs_difference(edge)

# TODO: Should be completely uniform in periodic case?
export simpson_weights
function simpson_weights(lattice::CompactLattice{T,1}) where T
    @assert all(size(lattice) .% 2 .== 1)
    w = ones(T, size(lattice)...)
    w[1] = 0.5
    w[end] = 0.5
    return w
end
function simpson_weights(lattice::CompactLattice{T,2}) where T
    # http://mathfaculty.fullerton.edu/mathews/n2003/SimpsonsRule2DMod.html
    @assert all(size(lattice) .% 2 .== 1)
    w = ones(T, size(lattice)...)
    w[2:2:end-1,:] .*= 4.0
    w[3:2:end-2,:] .*= 2.0
    w[:,2:2:end-1] .*= 4.0
    w[:,3:2:end-2] .*= 2.0
    return w
end
Base.getindex(lat::CompactLattice, dx) = getindex(lat.arr, dx)

# Use default lattice plotting


text:
module Tmp

f(n) = 2n

end


text:
module Utils

using Base.Meta: isexpr
using ..Specialized

function spec_lhs_rhs(ex)
    if isexpr(ex, :where)
        lhs, rhs = spec_lhs_rhs(ex.args[1])
        return Expr(:where, lhs, ex.args[2:end]...), rhs
    elseif isexpr(ex, :call)
        f = ex.args[1]
        @assert f isa Symbol
        spec = Expr(:call, :($Specialized.$f), ex.args[2:end]...)
        @debug "@def -> spec_lhs_rhs" ex spec
        rhs = strip_typeassert(spec)
        if isexpr(get(rhs.args, 2, nothing), :parameters)
            rhs.args[2].args .= map(rhs.args[2].args) do x
                if isexpr(x, :kw)
                    Expr(:kw, x.args[1], x.args[1])
                else
                    x
                end
            end
        end
        return (spec, rhs)
    else
        error("Cannot handle expression: ", ex)
    end
end

function strip_typeassert(ex)
    if isexpr(ex, :(::))
        length(ex.args) == 2 ||
            error("argument name is required for `@def`. Got:\n", ex)
        return ex.args[1]
    elseif ex isa Expr
        return Expr(ex.head, map(strip_typeassert, ex.args)...)
    else
        return ex
    end
end

"""
    @def f(...) = ...
    @def function f(...) ... end

Expand function definition for function `f` to:

```julia
@inline \$f(...) = Specialized.\$f(...)
@inline Specialized.\$f(...) = ...
```
"""
macro def(ex)
    @assert ex.head in [:function, :(=)]
    @assert length(ex.args) == 2
    lhs, rhs = spec_lhs_rhs(ex.args[1])
    block = Expr(:block, __source__, rhs)
    quote
        @inline $(ex.args[1]) = $block
        @inline $lhs = $(ex.args[2])
    end |> esc
end

end  # module
using .Utils: @def


text:
abstract type AbstractEstimation end

#contains the estimation problem
mutable struct Estimation <: AbstractEstimation
  Θ::Parameters
  Π::Priors
  F::Posteriors

  G::Portfolio
  TG::TrisectedPortfolio
  P::Portfolio

  nburn::Int
  nsamples::Int
  ncomplete::Int

  outparameters::Matrix{Float64}
  outweights::Matrix{Float64}
end

#creates the initial estimation object
function Estimation(U::Universe; Θ::Parameters = Parameters(),
  Π::Priors = Priors(),
  F::Posteriors = Posteriors(),
  nburn::Int = N_BURN,
  nsamples::Int = N_SAMPLES)::Estimation

  G::Portfolio = Portfolio(U)
  TG::TrisectedPortfolio = TrisectedPortfolio(G)
  P::Portfolio = initializetestportfolio(G)

  ncomplete::Int = 0

  outparameters::Matrix{Float64} = Matrix{Float64}(undef, length(copy2vector(Θ)), nsamples)
  outweights::Matrix{Float64} = Matrix{Float64}(undef, U.M, nsamples)

  return Estimation(Θ, Π, F, G, TG, P, nburn, nsamples, ncomplete, outparameters, outweights)
end

function Estimation(;datapath = DATA_PATH,
  dataname = DATA_NAME,
  datebegin::Date = DATE_BEGIN,
  dateend::Date = DATE_END,
  Θ::Parameters = Parameters(),
  Π::Priors = Priors(),
  F::Posteriors = Posteriors(),
  nburn::Int = N_BURN,
  nsamples::Int = N_SAMPLES)::Estimation

  return Estimation(Universe(datapath=datapath, dataname=dataname,
    datebegin=datebegin, dateend=dateend),
  Θ=Θ, Π=Π, F=F, nburn=nburn, nsamples=nsamples)
end

function mergeestimations!(E::Estimation...)::Estimation
  if length(E) == 1
    return E[1]
  elseif length(E) == 2
    E[1].nsamples += E[2].nsamples
    E[1].outparameters = [E[1].outparameters E[2].outparameters]
    E[1].outweights = [E[1].outweights E[2].outweights]
    return E[1]
  end

  #pivot and recurse
  pivot::Int = length(E) ÷ 2
  Eout::Estimation = mergeestimations!(
    mergeestimations!(E[1:pivot]...),
    mergeestimations!(E[(pivot+1):end]...))

  #print("length(E): $(length(E)) result: $(sum(Eout.outparameters))")
  return Eout
end

#runs a single estimation pass
function estimationpass!(E::Estimation)::Nothing

  redrawtestportfolio!(E.P)
  #E.Θ.SG = var(E.G) #mitigates error propogation
  E.Θ.SGP = cov(E.G, E.P) #where should I do this? try it here

  rotatetrisected!(E.TG)

  E.F.σ²G!(E.Θ, E.Π) #draw from the posteriors
  E.F.ζ²G!(E.Θ, E.Π)
  E.F.ζ²P!(E.Θ, E.Π)
  E.F.SG!(E.Θ, E.Π)
  E.F.SGP!(E.Θ, E.Π)

  #compute the implied weights
  ω::Vector{Float64} = getω(E.Θ, E.TG, E.P)
  rescaletrisected!(E.TG, ω)

  #record the output if we are past the burn-in
  E.ncomplete += 1
  if E.ncomplete > E.nburn
    sampleidx::Int = E.ncomplete - E.nburn

    E.outparameters[:, sampleidx] .= copy2vector(E.Θ)
    E.outweights[:, sampleidx] .= E.G.wₜ
  end

  return nothing
end

function recordestimation(E::Estimation,
  outname = OUT_NAME, outpath::String = OUT_PATH, record::Bool = true,
  timestamp::Bool = false)::Nothing

  if timestamp
    parameterspath::String = "$outpath\\$(outname)_par_$(now()).csv"
    weightspath::String = "$outpath\\$(outname)_w_$(now()).csv"
  else
    parameterspath = "$outpath\\$(outname)_par.csv"
    weightspath = "$outpath\\$(outname)_w.csv"
  end

  dfΘ::DataFrame = DataFrame(pass = collect(1:E.nsamples),
    sigma2G = E.outparameters[1,:], zeta2G = E.outparameters[2,:], zeta2P = E.outparameters[3,:],
    SG = E.outparameters[4,:], SGP = E.outparameters[5,:])
  dfw::DataFrame = DataFrame(E.outweights')

  ((oldn::Symbol, newn::Symbol)->rename!(dfw, oldn => newn)).(names(dfw), setdiff(names(E.G.U.df), [:date]))

  dfw = [DataFrame(pass = collect(1:E.nsamples)) dfw]

  CSV.write(parameterspath, dfΘ, missingstring="NA")
  CSV.write(weightspath, dfw, missingstring="NA")

  return nothing

end

#launches and performes the estimation
function estimate!(E::Estimation; verbose::Bool = true)::Estimation

  npasses::Int = E.nburn + E.nsamples

  verbose && println("Begining estimation....")
  for i::Int ∈ 1:npasses
    estimationpass!(E)
    verbose && (i % 10_000 == 0) && println("Completed $i passes out of $npasses")
  end
  verbose && println("Estimation complete and successful!")

  return E
end

#helper function which starts from the investment universe
function estimate(U::Universe; verbose::Bool = true, nburn::Int = N_BURN,
  nsamples::Int = N_SAMPLES)::Estimation

  E::Estimation = Estimation(U, nburn=nburn, nsamples=nsamples)
  estimate!(E, verbose=verbose)

  return E
end

function estimateminvariance(;datapath = DATA_PATH,
  dataname = DATA_NAME,
  datebegin::Date = DATE_BEGIN,
  dateend::Date = DATE_END,
  nburn::Int = N_BURN,
  nsamples::Int = N_SAMPLES,
  parallel::Bool = (nsamples > length(workers())*nburn), #simple heuristic for deciding if we will parllelize
  recordestimate::Bool=true)::Nothing

  U::Universe = Universe(datapath=datapath, dataname=dataname,
    datebegin=datebegin, dateend=dateend)
  local E::Estimation

  if !parallel
    E = estimate(U, nburn=nburn, nsamples=nsamples)
    recordestimate && recordestimation(E::Estimation)
    return E
  end

  #parllelized MCMC
  local pids = workers()
  local np::Int = length(pids)
  local targetassignmentsize::Int = nsamples ÷ np
  local assignments::Vector{Int} = Vector{Int}(undef, np)

  assignments[1:(end-1)] .= targetassignmentsize #divy up the assignments
  assignments[end] = nsamples - sum(assignments[1:(end-1)])

  nestimations::Int = length(assignments)
  futures::Vector{Future} = Vector{Future}()
  sizehint!(futures, np)

  #deploy the workers
  for i ∈ 1:np
    push!(futures,
      (@spawnat pids[i] estimate(U, nburn=nburn, nsamples=assignments[i])))
  end

  estimations::Vector{Estimation} = (fetch).(futures)
  estimationstest::Vector{Estimation} = deepcopy(estimations)

  E = mergeestimations!(estimations...)
  recordestimate && recordestimation(E)

  return nothing

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
"""
Provides a definition of a GROW method for constructing Koza trees.
"""
type KozaGrowBuilderDefinition <: KozaBuilderDefinition
  min_depth::Int
  max_depth::Int
  prob_terminal::Float

  KozaGrowBuilderDefinition() = new(1, 8, 0.5)
  KozaGrowBuilderDefinition(min::Int, max::Int) = new(min, max, 0.5)
  KozaGrowBuilderDefinition(min::Int, max::Int, p_terminal::Float) =
    new(min, max, p_terminal)
end

"""
Implements the GROW method for constructing Koza trees.
"""
type KozaGrowBuilder{T <: KozaTree} <: KozaBuilder
  min_depth::Int
  max_depth::Int
  prob_terminal::Float

  KozaGrowBuilder(min::Int, max::Int, p::Float) = new(min, max, p)
end

"""
Composes a Koza Grow builder, for a given sub-type of Koza Tree, using a
provided definition.
"""
function compose!(d::KozaGrowBuilderDefinition, t::Type)
  println("BUILDING GROW")
  KozaGrowBuilder{t}(d.min_depth, d.max_depth, d.prob_terminal)
end

grow() = KozaGrowBuilderDefinition()
grow(mn::Int, mx::Int) = KozaGrowBuilderDefinition(mn, mx)
function grow(min::Int, max::Int, p::Float)
  KozaGrowBuilderDefinition(min, max, p)
end
function grow(f::Function)
  def = KozaGrowBuilderDefinition()
  f(def)
  def
end

"""
Constructs a Koza Tree belonging to a given representation using the GROW
method, according to a set of provided specifics.
"""
function build{T}(b::KozaGrowBuilder{T}, r::KozaTreeRepresentation)
  t = T()
  t.root = build(KozaGrowBuilder, r, t, 1, rand(b.min_depth:b.max_depth), b.prob_terminal)
  return t
end

build(::Type{KozaGrowBuilder},
  r::KozaTreeRepresentation,
  p::KozaParent,
  d::Integer,
  md::Integer,
  pb::Float
) =
  if rand() <= pb || d == md
    return sample_terminal(r, p)
  else
    n = sample_non_terminal(r, p)
    for i in 1:arity(n)
      n.children[i] = build(KozaGrowBuilder, r, n, d+1, md, pb)
    end
    return n
  end


text:
using DataFrames
using CSV
using VegaLite

main = CSV.File("/home/graham_s/VirtualWorlds/projects/action_for_children/england/results/mainresults/by_sys_and_year.csv", delim='\t') |> DataFrame

main

function makeplot( main :: DataFrame, measure :: AbstractString )

    v = main |>
           @vlplot(
            title = "Projected Total Young People in Staying Put, 2019-2024",
            width=600,
            height=500,
            x={"year:o"}
           ) +
           @vlplot(
               mark={:area,color=:blue,opacity=0.1},
               title="",
               y={"min_payments_sys_1:q", title=""},
               y2={"max_payments_sys_1:q", title=""}
           )  +
           @vlplot(
               mark={:area,color=:blue,opacity=0.1},
               title="",
               y={"pct_10_payments_sys_1:q",title=""},
               y2={"pct_90_payments_sys_1:q",title=""}
           ) +
           @vlplot(
               mark={:area,color=:blue,opacity=0.3},
               title="",
               y={"pct_25_payments_sys_1:q", title=""},
               y2={"pct_75_payments_sys_1:q", title=""}
           ) +
           @vlplot(
               mark={:line,color=:black},
               y="avg_payments_sys_1:q",
           )+
           @vlplot(
               mark={:area,color=:red,opacity=0.1},
               title="",
               y={"min_payments_sys_2:q", title=""},
               y2={"max_payments_sys_2:q", title=""}
           )  +
           @vlplot(
               mark={:area,color=:red,opacity=0.1},
               title="",
               y={"pct_20_payments_sys_2:q",title=""},
               y2={"pct_90_payments_sys_2:q",title=""}
           ) +
           @vlplot(
               mark={:area,color=:red,opacity=0.3},
               title="",
               y={"pct_25_payments_sys_2:q", title=""},
               y2={"pct_75_payments_sys_2:q", title=""}
           ) +
           @vlplot(
               mark={:line,color=:black},
               y="avg_payments_sys_2:q",
           ) +
           @vlplot(
               mark={:area,color=:green,opacity=0.1},
               title="",
               y={"min_payments_sys_3:q", title=""},
               y2={"max_payments_sys_3:q", title=""}
           )  +
           @vlplot(
               mark={:area,color=:green,opacity=0.1},
               title="",
               y={"pct_30_payments_sys_3:q",title=""},
               y2={"pct_90_payments_sys_3:q",title=""}
           ) +
           @vlplot(
               mark={:area,color=:green,opacity=0.3},
               title="",
               y={"pct_35_payments_sys_3:q", title=""},
               y2={"pct_75_payments_sys_3:q", title=""}
           ) +
           @vlplot(
               mark={:line,color=:black},
               y="avg_payments_sys_3:q",
           )
    v
end

function makeplot2( main :: DataFrame, measure :: AbstractString )

    colours = [:red, :blue, :green, :orange, :purple]
    titles = [
    "current modelled",
    "option 2(a), with Oldham style fee",
    "option 2(a), without fee",
    "option 2(b), with HB from all benefit recipients",
    "option 3 - 2(a) with taper"
    ]
    v = main |>
           @vlplot(
            title = "Projected Annual Costs to Local Authorities of Staying Put, 2019-2024",
            width=600,
            height=500,
            legend= { :title="" },
            x={"year:o"}
           )
    for sys in 1:5
        v +=
           @vlplot(
               mark={:area,color=colours[sys],opacity=0.1},
               title="",
               y={"min_payments_sys_$sys:q", title=""},
               y2={"max_payments_sys_$sys:q", title=""}
           )  +
           @vlplot(
               mark={:area,color=colours[sys],opacity=0.1},
               title="",
               y={"pct_10_payments_sys_$sys:q",title=""},
               y2={"pct_90_payments_sys_$sys:q",title=""}
           ) +
           @vlplot(
               mark={:area,color=colours[sys],opacity=0.1},
               title="",
               y={"pct_25_payments_sys_$sys:q", title=""},
               y2={"pct_75_payments_sys_$sys:q", title=""}
           ) +
           @vlplot(
               mark={:line,color=colours[sys]},
               y={"avg_payments_sys_$sys:q", label=titles[sys]}
           )
       end
    v
end

save( "results/images/main_agg_2.svg", v )
save( "results/images/main_agg_2.png", v )
save( "results/images/main_agg_2.pdf", v )

names( main )

by_la = CSV.File("/home/graham_s/VirtualWorlds/projects/action_for_children/england/results/mainresults/by_la_sys_and_year_merged_with_grant.csv", delim='\t') |> DataFrame

main = CSV.File("/home/graham_s/VirtualWorlds/projects/action_for_children/england/results/mainresults/by_sys_and_year.csv", delim='\t') |> DataFrame


text:
## Logistic regression

module logistic_regression

using Distributions

struct data
    y::Vector
    x::Array{Float64}
    m::Vector
end

## Logistic regression nobs
function nobs(data::data)
    nx = size(data.x)[1]
    ny = length(data.y)
    nm = length(data.m)
    if (nx != ny)
        error("number of rows in of x is not equal to the length of y")
    elseif (nx != nm)
        error("number of rows in of x is not equal to the length of m")
    elseif (ny != nm)
        error("length of y is not equal to the length of m")
    end
    nx
end

function loglik(theta::Vector,
                data::data,
                i::Int64)
    eta = sum(data.x[i, :] .* theta)
    mu = exp.(eta)./(1 .+ exp.(eta))
    ll = data.y[i] .* log.(mu) + (data.m[i] - data.y[i]) .* log.(1 .- mu)
    ll
end

function ef(theta::Vector,
            data::data,
            i::Int64)
    eta = sum(data.x[i, :] .* theta)
    mu = exp.(eta)./(1 .+ exp.(eta))
    data.x[i, :] * (data.y[i] - data.m[i] * mu)
end

function simulate(theta::Vector,
                  x::Matrix{Float64},
                  m::Vector{Int64})
    n = size(x)[2]
    y = rand.(Binomial.(m, cdf.(Logistic(), x * theta)))
    data(y, x, m);
end

function loglik_alt(theta::Vector,
                    data::data,
                    i::Int64)
    eta = sum(data.x[i, :] .* theta)
    mu = exp.(eta)./(1 .+ exp.(eta))
    ll = data.y[i] .* log.(mu) + (data.m[i] - data.y[i]) .* log.(1 .- mu)
    ll - 10 * sum(theta.^2) / length(data.y)
end

function l2_penalty_obj(theta::Vector,
                        data::data)
    - 10 * sum(theta.^2)
end

function l2_penalty_ef(theta::Vector,
                       data::data)
    - 20 * theta
end

end


text:
function optimfunwrapper(x::Vector, g::Vector, site, var)
    g === nothing && (g = zeros(Float64, length(x)))
    return PLsiteAndGrad!(x, g, site,  var)
end

function optimfunwrapper(x::Vector, g::Vector, var)
    g === nothing && (g = zeros(Float64, length(x)))
    return PLsiteAndGradSym!(x, g, var)
end

function ComputeScore(Jmat::Array{Float64,2}, var::PlmVar, min_separation::Int)

    q = var.q
    N = var.N

    JJ=reshape(Jmat[1:end-q,:], q,q,N-1,N)

    Jtemp1=zeros( q,q,Int(N*(N-1)/2))
    Jtemp2=zeros( q,q,Int(N*(N-1)/2))
    l = 1

    for i=1:(N-1)
        for j=(i+1):N
            Jtemp1[:,:,l]=JJ[:,:,j-1,i]; #J_ij as estimated from from g_i.
            Jtemp2[:,:,l]=JJ[:,:,i,j]'; #J_ij as estimated from from g_j.
            l=l+1;
        end
    end

    Jtensor = zeros(q,q,N,N)
    l = 1
    for i = 1:N-1
        for j=i+1:N
            Jtensor[:,:,i,j] = Jtemp1[:,:,l]
            Jtensor[:,:,j,i] = Jtemp2[:,:,l]'
            l += 1
        end
    end

    ASFN = zeros(N,N)
    for i=1:N,j=1:N
        i!=j && (ASFN[i,j] =sum(Jtensor[:,:,i,j].^2))
    end

    J1=zeros(q,q,Int(N*(N-1)/2))
    J2=zeros(q,q,Int(N*(N-1)/2))

    for l=1:Int(N*(N-1)/2)
        J1[:,:,l] = Jtemp1[:,:,l]-repeat(mean(Jtemp1[:,:,l], dims=1),q,1)-repeat(mean(Jtemp1[:,:,l], dims=2),1,q) .+ mean(Jtemp1[:,:,l])
        J2[:,:,l] = Jtemp2[:,:,l]-repeat(mean(Jtemp2[:,:,l], dims=1),q,1)-repeat(mean(Jtemp2[:,:,l], dims=2),1,q) .+ mean(Jtemp2[:,:,l])
    end
    J = 0.5 * ( J1 + J2 )

    htensor = fill(0.0, q,N)
    for i in 1:N
        htensor[:,i] = Jmat[end-q+1:end,i]
    end

    FN = zeros(Float64, N,N)
    l = 1

    for i=1:N-1
        for j=i+1:N
            # norm(A, p::Real=2)
            # If A is a matrix and p=2, then this is equivalent to the Frobenius norm.
            FN[i,j] = norm(J[1:q-1,1:q-1,l], 2)
#            FN[i,j] = vecnorm(J[:,:,l],2)
            FN[j,i] =FN[i,j]
            l+=1
        end
    end
    FN=GaussDCA.correct_APC(FN)
    score = GaussDCA.compute_ranking(FN,min_separation)
    return score, FN, 0.5*(permutedims(Jtensor,[2,1,4,3])+Jtensor), htensor
end

function ReadFasta(filename::AbstractString,max_gap_fraction::Real, theta::Any, remove_dups::Bool)

    Z = GaussDCA.read_fasta_alignment(filename, max_gap_fraction)
    if remove_dups
        Z, _ = GaussDCA.remove_duplicate_seqs(Z)
    end

    N, M = size(Z)
    q = round(Int,maximum(Z))

    q > 32 && error("parameter q=$q is too big (max 31 is allowed)")
    W , Meff = GaussDCA.compute_weights(Z,q,theta)
    W .*= 1.0/Meff
    Zint=round.(Int,Z)
    return W, Zint,N,M,q
end

function sumexp(vec::Array{Float64,1})
    mysum = 0.0
    @inbounds @simd for i=1:length(vec)
        mysum += exp(vec[i])
    end
    return mysum
end

function Base.deepcopy_internal(x::DecVar, d::IdDict)
    haskey(d, x) && return d[x]
    dmask_c = Base.deepcopy_internal(x.dmask, d)
    xc = DecVar(x.fracdec, x.fracmax, dmask_c)
    d[x] = xc
    return xc
end

max_stars_repo_name: YoshitakaMo/PlmDCA

max_issues_repo_name: YoshitakaMo/PlmDCA

max_forks_repo_name: YoshitakaMo/PlmDCA


text:
# See https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array-1
using Revise, BinaryTraits

import Base: size, getindex, setindex!

# -----------------------------------------------------------------------------
# AbstractArray interface
# -----------------------------------------------------------------------------

@trait Dimension prefix Has,No
@implement HasDimension by size()::Tuple

# Bottom is used due to duck typing of the value `v` below
const Bottom = Base.Bottom

@trait LinearIndexing prefix Has,No
@implement HasLinearIndexing by getindex(i::Int)
@implement HasLinearIndexing by setindex!(v::Bottom, i::Int)

const IntVarArg = Vararg{Int, N} where N

@trait CartesianIndexing prefix Has,No
@implement HasCartesianIndexing by getindex(I::IntVarArg)
@implement HasCartesianIndexing by setindex!(v::Bottom, I::IntVarArg)

# -----------------------------------------------------------------------------
# Example: 1-D Int array
# -----------------------------------------------------------------------------
const Array1DInt = Array{Int,1}

@assign Array1DInt with HasDimension
@check(Array1DInt)
#=
julia> @check(Array1DInt)
✅ Array{Int64,1} has implemented:
1. DimensionTrait: HasDimension ⇢ size(::<Type>)::Tuple
=#

@assign Array1DInt with HasLinearIndexing
@check(Array1DInt)
#=
julia> @check(Array1DInt)
✅ Array{Int64,1} has implemented:
1. LinearIndexingTrait: HasLinearIndexing ⇢ getindex(::<Type>, ::Int64)::Any
2. LinearIndexingTrait: HasLinearIndexing ⇢ setindex!(::<Type>, ::Union{}, ::Int64)::Any
3. DimensionTrait: HasDimension ⇢ size(::<Type>)::Tuple
=#

# 1D array is a specialized version of CartesianIndexing.
# Let's verify.
@assign Array1DInt with HasCartesianIndexing
@check(Array1DInt)
#=
julia> @check(Array1DInt)
✅ Array{Int64,1} has implemented:
1. LinearIndexingTrait: HasLinearIndexing ⇢ getindex(::<Type>, ::Int64)::Any
2. LinearIndexingTrait: HasLinearIndexing ⇢ setindex!(::<Type>, ::Union{}, ::Int64)::Any
3. CartesianIndexingTrait: HasCartesianIndexing ⇢ getindex(::<Type>, ::Vararg{Int64,N} where N)::Any
4. CartesianIndexingTrait: HasCartesianIndexing ⇢ setindex!(::<Type>, ::Union{}, ::Vararg{Int64,N} where N)::Any
5. DimensionTrait: HasDimension ⇢ size(::<Type>)::Tuple
=#

# -----------------------------------------------------------------------------
# Example: SquaresVector
# -----------------------------------------------------------------------------

struct SquaresVector <: AbstractArray{Int, 1}
    count::Int
end

Base.size(S::SquaresVector) = (S.count,)
Base.getindex(S::SquaresVector, i::Int) = i*i

@assign SquaresVector with HasDimension,HasLinearIndexing
@check(SquaresVector)
#=
julia> @check(SquaresVector)
✅ SquaresVector has implemented:
1. LinearIndexingTrait: HasLinearIndexing ⇢ getindex(::<Type>, ::Int64)::Any
2. LinearIndexingTrait: HasLinearIndexing ⇢ setindex!(::<Type>, ::Union{}, ::Int64)::Any
3. DimensionTrait: HasDimension ⇢ size(::<Type>)::Tuple
=#


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
#
# This file depends upon sudoku_common.jl

"""
    row_block_permute(puzzle,row_block_a,row_block_b)

Swap entire block-rows within a SolvablePuzzle.
"""
function row_block_permute!(puzzle::Matrix,row_block_a::Integer,row_block_b::Integer)
    rank=get_rank(puzzle)
    if row_block_a < 1 || row_block_a > rank
        throw(DomainError("Invalid row block A"))
    elseif row_block_b < 1 || row_block_b > rank
        throw(DomainError("Invalid row block B"))
    end
    start_row_a = rank*(row_block_a-1)+1
    stop_row_a = rank*(row_block_a-1)+rank
    start_row_b = rank*(row_block_b-1)+1
    stop_row_b = rank*(row_block_b-1)+rank
    rank_squared = rank*rank
    # copy block A from source
    block = puzzle[start_row_a:stop_row_a,:]
    # write block B into A
    puzzle[start_row_a:stop_row_a,:] = puzzle[start_row_b:stop_row_b,:]
    # write block A into B
    puzzle[start_row_b:stop_row_b,:] = block[:,:]
end

"""
    col_block_permute(puzzle,col_block_a,col_block_b)

Swap entire block-columns within a SolvablePuzzle.
"""
function col_block_permute!(puzzle::Matrix,col_block_a::Integer,col_block_b::Integer)
    rank = get_rank(puzzle)
    if col_block_a < 1 || col_block_a > rank
        throw(DomainError("Invalid col block A"))
    elseif col_block_b < 1 || col_block_b > rank
        throw(DomainError("Invalid col block B"))
    end
    start_col_a = rank*(col_block_a-1)+1
    stop_col_a = rank*(col_block_a-1)+rank
    start_col_b = rank*(col_block_b-1)+1
    stop_col_b = rank*(col_block_b-1)+rank
    rank_squared = rank*rank
    # copy block A from source
    block = puzzle[:,start_col_a:stop_col_a]
    # write block B into A
    puzzle[:,start_col_a:stop_col_a] = puzzle[:,start_col_b:stop_col_b]
    # write block A into B
    puzzle[:,start_col_b:stop_col_b] = block[:,:]
end

"""
    row_permute!(puzzle,row_block,subrow_a,subrow_b)

Swap rows within a SolvablePuzzle block column.
"""
function row_permute!(puzzle::Matrix,row_block::Integer,subrow_a::Integer,subrow_b::Integer)
    rank = get_rank(puzzle)
    if row_block < 1 || row_block > rank
        throw(DomainError("Invalid row block"))
    elseif subrow_a < 1 || subrow_a > rank
        throw(DomainError("Invalid subrow A"))
    elseif subrow_b < 1 || subrow_b > rank
        throw(DomainError("Invalid subrow B"))
    end
    start_row = rank*(row_block-1)+1
    stop_row = rank*(row_block-1)+rank
    permutation = [i for i in 1:rank]
    permutation[subrow_a] = subrow_b
    permutation[subrow_b] = subrow_a
    # copy permuted rows from source
    block = puzzle[(start_row-1).+permutation,:]
    # update source with permuted rows
    puzzle[start_row:stop_row,:] = block[:,:]
end

"""
    col_permute!(puzzle,col_block,subcol_a,subcol_b)

Swap columns within a SolvablePuzzle block column.
"""
function col_permute!(puzzle::Matrix,col_block::Integer,subcol_a::Integer,subcol_b::Integer)
    rank=get_rank(puzzle)
    if col_block < 1 || col_block > rank
        throw(DomainError("Invalid col block"))
    elseif subcol_a < 1 || subcol_a > rank
        throw(DomainError("Invalid subcol A"))
    elseif subcol_b < 1 || subcol_b > rank
        throw(DomainError("Invalid subcol B"))
    end
    start_col = rank*(col_block-1)+1
    stop_col = rank*(col_block-1)+rank
    permutation = [i for i in 1:rank]
    permutation[subcol_a] = subcol_b
    permutation[subcol_b] = subcol_a
    # copy permuted rows from source
    block = puzzle[:,(start_col-1).+permutation]
    # update source with permuted rows
    puzzle[:,start_col:stop_col] = block[:,:]
end

"""
    mirror_horizontal!(puzzle)

Swap columns across a SolvablePuzzle.
"""
function mirror_horizontal!(puzzle::Matrix)
    rank=get_rank(puzzle)
    if rank > 1 # blockwise swapping
        half_rank = Integer(floor(rank/2))
        for i = 1:half_rank # swap blocks
            col_block_permute!(puzzle,i,rank-i+1)
            for j = 1:half_rank # swap rows in blocks
                col_permute!(puzzle,i,j,rank-j+1)
                col_permute!(puzzle,rank-i+1,j,rank-j+1)
            end # end subrow loop
        end # end block loop
    end # end rank switch
end

"""
    mirror_vertical!(puzzle)

Swap rows across a SolvablePuzzle.
"""
function mirror_vertical!(puzzle::Matrix)
    rank=get_rank(puzzle)
    # no mirroring for rank 1
    if rank > 1 # blockwise swapping
        half_rank = Integer(floor(rank/2))
        for i = 1:half_rank # swap blocks
            row_block_permute!(puzzle,i,rank-i+1)
            for j = 1:half_rank # swap rows in blocks
                row_permute!(puzzle,i,j,rank-j+1)
                row_permute!(puzzle,rank-i+1,j,rank-j+1)
            end # end subrow loop
        end # end block loop
    end # end rank switch
end

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
# This is part of Quiz.jl. License is MIT.

# 64-bit pod real numbers
N64 = Union{Float64, Int64}

# Unit-ed 64-bit types
Q64 = Quantity{<:N64}

# All 64-bit pod and unit-ed types
A64 = Union{N64, Q64}

export N64, Q64, A64

n64(x::Integer) = Int64(x)
n64(x::Real) = Float64(x)

function val_unit(x::Number)
    u = unit(x)
    u isa Unitful.AffineUnits ? (n64(x.val), u) : (n64(x / u), u)
end


text:
@unix_only run(`make`)
# @windows_only Int == Int32 ? run(`mingw32-make`) : run(`mingw32-make`)
@windows_only run(`mingw32-make`)


text:
# sexp uses helper functions defined in PhyloNetworks

@doc doc"""
    function sexp(net::HybridNetwork)

Export a HybridNework object to the R language
as either `phylo` or `evonet` object (depending on degree of hybridization)
recognized by the R package `ape`.
Used by the `$object` syntax and by `@rput` to use a Julia object in R:
see the examples below. Makes it easy to plot a Julia tree or network
using plotting facilities in R.

code inspired from [Phylo.jl](https://github.com/richardreeve/Phylo.jl/blob/master/src/rcall.jl)

# Examples

```julia-repl
julia> using RCall
julia> using PhyloNetworks
julia> net = readTopology("(((A:.2,(B:.1)#H1:.1::0.9):.1,(C:.11,#H1:.01::0.1):.19):.1,D:.4);");
R> library(ape); # type $ to switch from julia to R
R> $net

Evolutionary network with 1 reticulation

               --- Base tree ---
Phylogenetic tree with 4 tips and 5 internal nodes.

Tip labels:
[1] "A" "B" "C" "D"

Rooted; includes branch lengths.
julia> @rput net # press the delete key to switch from R back to julia
R> net

Evolutionary network with 1 reticulation

--- Base tree ---
Phylogenetic tree with 4 tips and 5 internal nodes.

Tip labels:
[1] "A" "B" "C" "D"

Rooted; includes branch lengths.

R> str(net)
List of 7
 $ edge               : int [1:8, 1:2] 5 5 6 6 7 8 8 9 6 4 ...
 $ reticulation.length: num 0.01
 $ Nnode              : int 5
 $ edge.length        : num [1:8] 0.1 0.4 0.1 0.19 0.11 0.2 0.1 0.1
 $ reticulation       : int [1, 1:2] 7 9
 $ reticulation.gamma : num 0.1
 $ tip.label          : chr [1:4] "A" "B" "C" "D"
 - attr(*, "class")= chr [1:2] "evonet" "phylo"
NULL

R> plot(net)
```
""" #"
function sexp(net::HybridNetwork)
    PhyloNetworks.resetNodeNumbers!(net)
    ntips = length(net.leaf)
    totalnodes = length(net.node)
    Nnode = totalnodes - ntips
    o = sortperm([n.number for n in net.leaf])
    tipLabel = [net.leaf[i].name for i in o]
    edge = PhyloNetworks.majoredgematrix(net)
    phy = Dict{Symbol, Any}() # dictionary exported with regular sexp at the end
    phy[:Nnode] = Nnode
    phy[Symbol("tip.label")] = tipLabel
    phy[:edge] = edge
    edgeLength = PhyloNetworks.majoredgelength(net)
    if any(.!ismissing.(edgeLength))
        phy[Symbol("edge.length")] = edgeLength
    end
    if net.numHybrids > 0
        reticulation = PhyloNetworks.minorreticulationmatrix(net)
        reticulationGamma = PhyloNetworks.minorreticulationgamma(net)
        reticulationLength = PhyloNetworks.minorreticulationlength(net)
        phy[:reticulation] = reticulation
        if any(.!ismissing.(reticulationGamma))
            phy[Symbol("reticulation.gamma")] = reticulationGamma
        end
        if any(.!ismissing.(reticulationLength))
            phy[Symbol("reticulation.length")] = reticulationLength
        end
    end
    sobj = RCall.protect(sexp(phy)) # RObject
    if net.numHybrids == 0
        setclass!(sobj, sexp("phylo"))
    else
        setclass!(sobj, sexp(["evonet", "phylo"]))
    end
    RCall.unprotect(1)
    return(sobj)
end

@doc doc"""
    rexport(net::HybridNetwork; mainTree=false, useEdgeLength=true)

Create an RObject of class `phylo` (and `evonet` depending on the number
of hybridizations) recognized by the `ape` library in R (S3 object). This
RObject can be evaluated using the tools available in the `ape` library in R.
For example, we can visualize the network using `ape`'s `plot` function.

not exported: [`sexp`](@ref) is the best way to go.

# Arguments

- useEdgeLength: if true, export edge lengths from `net`.
- mainTree: if true, minor hybrid edges are omitted.

# Examples

```julia-repl
julia> net = readTopology("(((A,(B)#H1:::0.9),(C,#H1:::0.1)),D);");
julia> phy = rexport(net)
RCall.RObject{RCall.VecSxp}
$Nnode
[1] 5

$edge
     [,1] [,2]
[1,]    5    6
[2,]    5    4
[3,]    6    8
[4,]    6    7
[5,]    7    3
[6,]    8    1
[7,]    8    9
[8,]    9    2

$tip.label
[1] "A" "B" "C" "D"

$reticulation
     [,1] [,2]
[1,]    7    9

$reticulation.gamma
[1] 0.1

attr(,"class")
[1] "evonet" "phylo"

julia> using RCall

julia> R"library(ape)";

julia> phy
RCall.RObject{RCall.VecSxp}

    Evolutionary network with 1 reticulation

               --- Base tree ---
Phylogenetic tree with 4 tips and 5 internal nodes.

Tip labels:
[1] "A" "B" "C" "D"

Rooted; no branch lengths.

R> phy

Evolutionary network with 1 reticulation

               --- Base tree ---
Phylogenetic tree with 4 tips and 5 internal nodes.

Tip labels:
[1] "A" "B" "C" "D"

Rooted; no branch lengths.

R> str(phy)
List of 5
$ Nnode             : int 5
$ edge              : int [1:8, 1:2] 5 5 6 6 7 8 8 9 6 4 ...
$ tip.label         : chr [1:4] "A" "B" "C" "D"
$ reticulation      : int [1, 1:2] 7 9
$ reticulation.gamma: num 0.1
- attr(*, "class")= chr [1:2] "evonet" "phylo"
```
""" #"
function rexport(net::HybridNetwork; mainTree::Bool=false, useEdgeLength::Bool=true)
# worry about R object created within the function not accessible from outside:
# can it be garbage collected?

    if mainTree == true && net.numHybrids > 0
        net = majorTree(net)
    end
    PhyloNetworks.resetNodeNumbers!(net)
    ntips = length(net.leaf)
    totalnodes = length(net.node)
    Nnode = totalnodes - ntips
    o = sortperm([n.number for n in net.leaf])
    tipLabel = [net.leaf[i].name for i in o]
    edge = PhyloNetworks.majoredgematrix(net)
    R"""
    phy = list(Nnode = $Nnode, edge = $edge, tip.label = $tipLabel)
    """
    if useEdgeLength == true
        edgeLength = PhyloNetworks.majoredgelength(net)
        if any(.!ismissing.(edgeLength))
            R"""
            phy[['edge.length']] = $edgeLength
            """
        end
    end
    if net.numHybrids > 0
        reticulation = PhyloNetworks.minorreticulationmatrix(net)
        reticulationGamma = PhyloNetworks.minorreticulationgamma(net)
        R"""
        phy[['reticulation']] = $reticulation
        class(phy) <- c("evonet", "phylo")
        """
        if any(.!ismissing.(reticulationGamma))
            R"phy[['reticulation.gamma']] = $reticulationGamma"
        end
        if useEdgeLength # extract minor edge lengths
            reticulationLength = PhyloNetworks.minorreticulationlength(net)
            if any(.!ismissing.(reticulationLength))
                R"""
                phy[['reticulation.length']] = $reticulationLength
                """
            end
        end
    elseif net.numHybrids == 0
        R"""
        class(phy) = "phylo"
        """
    end
    phy = reval("phy")
    return(phy)
end


text:
using SqState
using Test

@testset "SqState.jl" begin

    include("read.jl")
    include("utils.jl")
    include("polynomial.jl")
    include("wigner.jl")

end


text:
using Documenter
using Arrow

makedocs(;
    modules=[Arrow],
    repo="https://github.com/JuliaData/Arrow.jl/blob/{commit}{path}#L{line}",
    sitename="Arrow.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaData.github.io/Arrow.jl",
        assets=String[],
    ),
    pages = [
        "Home" => "index.md",
        "User Manual" => "manual.md",
        "API Reference" => "reference.md"
    ]
)

deploydocs(;
    repo="github.com/JuliaData/Arrow.jl",
)


text:
export PandaBot
using PandaRobot
include("../dynamics/pandadynamics.jl")

using RigidBodySim
using RigidBodyDynamics

mutable struct PandaBot{T<:AbstractFloat} <: Robot
  # Robot Properties
  total_mass::T
  n_motors::Int
  # J::Matrix{T}
  # Jinv::Matrix{T}
  # n_thrusters::Int
  r::T
  # hard_limit_vel::T
  # hard_limit_accel::T
  # hard_limit_omega::T
  # hard_limit_alpha::T
  btCollisionObject

  # PandaRobot structure
  pan#::Panda <: PandaRobot

  # Robot state
  q # joints # same as RigidBodyDynamics.configuration(state)
  q_dot
  q_ddot
  state # state (includes joints)

  # Frame
  world_frame

  # End effector
  EE_id::Int
  EE_link
  EE_link_frame
  EE_link_point
  EE_link_radius::T
  EE_path
  EE_point3d_inWorldframe_limit_min
  EE_point3d_inWorldframe_limit_max
  # EE_Jp
end
function PandaBot{T}() where T
  # Robot Properties
  total_mass = 0.0
  n_motors = 9
  # s = 0.5*0.305   # each side of cube is 30.5cm
  r = sqrt(3)*0.5*0.305   # inflate to sphere
  # hard_limit_vel = 0.5
  # hard_limit_accel = 0.1
  # hard_limit_omega = 45*π/180
  # hard_limit_alpha = 50*π/180
  # J = 0.1083*Eye(3)
  # Jinv = inv(J)
  btCollisionObject = BT.sphere(SVector{3}(zeros(T,3)), r) # Update this

  # PandaRobot structure
  # Note that it includes the definition of all the links and their relative definition
  pan = Panda()

  # Robot state
  q = zeros(n_motors)
  q_dot = zeros(n_motors)
  q_ddot = zeros(n_motors)
  state = RigidBodyDynamics.MechanismState(pan.mechanism)
    set_configuration!(state,q)
    set_velocity!(state,q_dot)

  # Frame
  world_frame = RigidBodyDynamics.root_frame(pan.mechanism)

  # End effector (hand)
  EE_id = 11
  EE_link = RigidBodyDynamics.bodies(pan.mechanism)[EE_id]
  EE_link_frame = RigidBodyDynamics.default_frame(EE_link)
  EE_link_point = RigidBodyDynamics.Point3D(EE_link_frame,0.0,0.0,0.1)
  EE_link_radius = 0.02
  EE_path = RigidBodyDynamics.path(pan.mechanism, root_body(pan.mechanism),EE_link)
     # Get limits
    q0 = [0.0;0.0;0.0;0.0;0.0;π;0.0;0.0;0.0]
    set_configuration!(state,q0)
    lims_up = RigidBodyDynamics.transform(state, EE_link_point, world_frame)
  EE_point3d_inWorldframe_limit_min, EE_point3d_inWorldframe_limit_max = -maximum(lims_up.v)*ones(3), maximum(lims_up.v)*ones(3)
  # Jacobian (joints to point)
  # EE_Jp = point_jacobian(state, EE_path, RigidBodyDynamics.transform(state, EE_link_point, world_frame))  )

  # new PandaBot instance
  return PandaBot{T}(total_mass,n_motors,r,btCollisionObject,pan,q,q_dot,q_ddot,state,world_frame,
    EE_id, EE_link, EE_link_frame, EE_link_point, EE_link_radius, EE_path,
    EE_point3d_inWorldframe_limit_min, EE_point3d_inWorldframe_limit_max)
  # return PandaBot{T}(total_mass,n_motors,pan,q,state)
end
PandaBot(::Type{T} = Float64; kwargs...) where {T} = PandaBot{T}(; kwargs...)

##################
# Initializations
##################
# Careful, x_init, x_goal are in configuration space!
function init_traj_straightline(TOP::TrajectoryOptimizationProblem{PandaBot{T}, PandaDyn, E}) where {T,E}
  model, x_init, x_goal = TOP.PD.model, TOP.PD.x_init, TOP.PD.x_goal
  x_dim, u_dim, N, tf_guess = model.x_dim, model.u_dim, TOP.N, TOP.tf_guess
  N = TOP.N

  X = hcat(linspace(x_init, x_goal, N)...)
  U = zeros(u_dim, N)
  Trajectory(X, U, tf_guess)
end

##################
# Access Values
##################
function get_EE_point_inWorldFrame(manipulator::PandaBot)
  EE_point_inWorldFrame = RigidBodyDynamics.transform(manipulator.state, manipulator.EE_link_point, manipulator.world_frame)
  return EE_point_inWorldFrame
end
function get_mass_matrix(manipulator::PandaBot)
  mass_matrix(manipulator.state)
end
function get_torques(manipulator::PandaBot)
  inverse_dynamics(manipulator.state, manipulator.q_dot)
end

##################
# Update Values
##################
function update_state_robot(manipulator::PandaBot)
      set_configuration!(manipulator.state, manipulator.q)
      # set_configuration!(mvis, configuration(state))
end
function update_q_robot(manipulator::PandaBot)
      manipulator.q = manipulator.state.q
      # set_configuration!(mvis, configuration(state))
end


text:
module SequenceTests

using Test
using WeaklyHard

##############################
### Generate all sequences ###
##############################

function generate_all_seq(k)
    return [(1 << k) | i for i in 0:2^k-1]
end # function

function generate_all_valid(lambda::T, seqs::Vector{R}, k) where {T <: Constraint, R <: Integer}
    IntType = BigInt
    if k + 1 < 8
        IntType = Int8
    elseif k + 1 < 16
        IntType = Int16
    elseif k + 1 < 32
        IntType = Int32
    elseif k + 1 < 64
        IntType = Int64
    end
    cmp_w = IntType(2^k - 1)
    return Set{IntType}(map(x -> IntType(x) & IntType(cmp_w), filter(x -> is_satisfied(lambda, x), seqs)))
end # function
function generate_all_valid(Lambda::Set{T}, seqs::Vector{R}, k) where {T <: Constraint, R <: Integer}
    IntType = BigInt
    if k + 1 < 8
        IntType = Int8
    elseif k + 1 < 16
        IntType = Int16
    elseif k + 1 < 32
        IntType = Int32
    elseif k + 1 < 64
        IntType = Int64
    end
    cmp_w = IntType(2^k - 1)
    return Set{IntType}(map(x -> IntType(x) & IntType(cmp_w), filter(x -> is_satisfied(Lambda, x), seqs)))
end # function

##########################
### All Sequence tests ###
##########################

k_max = 15

@testset "Test all sequences generated by automaton - Not RowHitConstraint" begin

    for k in 1:k_max

        IntType = BigInt
        if k < 8
            IntType = Int8
        elseif k < 16
            IntType = Int16
        elseif k < 32
            IntType = Int32
        elseif k < 64
            IntType = Int64
        end

        all_seqs    = generate_all_seq(k)

        for x in 0:k
            lambda      = AnyMissConstraint(x, k)
            all_valid   = generate_all_valid(lambda, all_seqs, k)
            g           = build_automaton(lambda)
            all_g_seq   = all_sequences(g, k)
            @test isempty(setdiff(all_valid, all_g_seq))
            @test isempty(setdiff(all_g_seq, all_valid))

            lambda      = AnyHitConstraint(x, k)
            all_valid   = generate_all_valid(lambda, all_seqs, k)
            g           = build_automaton(lambda)
            all_g_seq   = all_sequences(g, k)
            @test isempty(setdiff(all_valid, all_g_seq))
            @test isempty(setdiff(all_g_seq, all_valid))

            if x < k
                lambda      = RowMissConstraint(x)
                all_valid   = generate_all_valid(lambda, all_seqs, k)
                g           = build_automaton(lambda)
                all_g_seq   = all_sequences(g, k)
                @test isempty(setdiff(all_valid, all_g_seq))
                @test isempty(setdiff(all_g_seq, all_valid))
            end
        end # for
    end # for
end # testset

@testset "Test all sequences generated by automaton - Only RowHitConstraint" begin

    for k in 1:k_max

        IntType = BigInt
        if 2*k + 2 < 8
            IntType = Int8
        elseif 2*k + 2 < 16
            IntType = Int16
        elseif 2*k + 2 < 32
            IntType = Int32
        elseif 2*k + 2 < 64
            IntType = Int64
        end

        all_seqs    = generate_all_seq(k)

        for x in 0:ceil(k/2)
            lambda      = RowHitConstraint(x, k)
            all_valid   = generate_all_valid(lambda, all_seqs, k)
            g           = build_automaton(lambda)
            all_g_seq   = all_sequences(g, k)
            @test isempty(setdiff(all_valid, all_g_seq))
            @test isempty(setdiff(all_g_seq, all_valid))
        end # for
    end #for
end # testset

@testset "Test all sequences generated by automaton - Sets of All but RowHitConstraint" begin

    k       = k_max
    n_tests = 100
    n_cstr  = 4
    constraint_names = ["AnyMissConstraint", "AnyHitConstraint", "RowMissConstraint"]

    all_seqs    = generate_all_seq(k)

    for n_c in 2:n_cstr
        for _ in 1:n_tests

            Lambda = Set{Constraint}()
            while length(Lambda) < n_c
                Lambda = Set{Constraint}()
                for _ in 1:n_c
                    lk = rand(1:k)
                    lx = rand(0:lk)
                    lambda_type = rand(constraint_names)
                    l = getfield(WeaklyHard, Symbol(lambda_type))(lx, lk)
                    push!(Lambda, l)
                end # for
            end # while

            all_valid   = generate_all_valid(Lambda, all_seqs, k)
            g           = build_automaton(Lambda)
            minimize_automaton!(g)
            all_g_seq   = all_sequences(g, k)
            @test isempty(setdiff(all_valid, all_g_seq))
            @test isempty(setdiff(all_g_seq, all_valid))
            if !isempty(setdiff(all_valid, all_g_seq)) || !isempty(setdiff(all_g_seq, all_valid))
                println(Lambda)
            end
        end # for
    end # for
end # testset

@testset "Test all sequences generated by automaton - Only Sets of RowHitConstraint" begin

    k       = k_max
    n_tests = 100
    n_cstr  = 4

    all_seqs    = generate_all_seq(k)

    for n_c in 2:n_cstr
        for _ in 1:n_tests

            Lambda = Set{RowHitConstraint}()
            while length(Lambda) < n_c
                Lambda = Set{RowHitConstraint}()
                for _ in 1:n_c
                    lk = rand(1:k)
                    lx = rand(0:lk)
                    l = RowHitConstraint(lx, lk)
                    push!(Lambda, l)
                end # for
            end # while

            all_valid   = generate_all_valid(Lambda, all_seqs, k)
            g           = build_automaton(Lambda)
            minimize_automaton!(g)
            all_g_seq   = all_sequences(g, k)
            @test isempty(setdiff(all_valid, all_g_seq))
            @test isempty(setdiff(all_g_seq, all_valid))
            if !isempty(setdiff(all_valid, all_g_seq)) || !isempty(setdiff(all_g_seq, all_valid))
                println(Lambda)
            end
        end # for
    end # for
end # testset

@testset "Test all sequences generated by automaton - Arbitrary sets" begin

    k       = k_max
    n_tests = 100
    n_cstr  = 4
    constraint_names = ["AnyMissConstraint", "AnyHitConstraint", "RowMissConstraint", "RowHitConstraint"]

    all_seqs    = generate_all_seq(k)

    for n_c in 2:n_cstr
        for _ in 1:n_tests

            Lambda = Set{Constraint}()
            while length(Lambda) < n_c
                Lambda = Set{Constraint}()
                for _ in 1:n_c
                    lk = rand(1:k)
                    lx = rand(0:lk)
                    lambda_type = rand(constraint_names)
                    l = getfield(WeaklyHard, Symbol(lambda_type))(lx, lk)
                    push!(Lambda, l)
                end # for
            end # while

            all_valid   = generate_all_valid(Lambda, all_seqs, k)
            g           = build_automaton(Lambda)
            minimize_automaton!(g)
            all_g_seq   = all_sequences(g, k)
            @test isempty(setdiff(all_valid, all_g_seq))
            @test isempty(setdiff(all_g_seq, all_valid))
            if !isempty(setdiff(all_valid, all_g_seq)) || !isempty(setdiff(all_g_seq, all_valid))
                println(Lambda)
            end
        end # for
    end # for
end # testset

######################
### Sequence tests ###
######################

N = 200

@testset "Random sequence tests - AnyMissConstraint" begin
    n_max = 0
    for k in 1:k_max
        for x in 0:k
            lambda      = AnyMissConstraint(x, k)
            automaton   = build_automaton(lambda)
            for n in k_max:N
                seq         = random_sequence(automaton, n)
                @test is_satisfied(lambda, seq)
            end # for

            n_max       = max(n_max, length(automaton))
        end # for
    end # for
    println("Biggest automaton for AnyMissConstraint: $n_max")
end # testset
@testset "Random sequence tests - AnyHitConstraint" begin
    n_max = 0
    for k in 1:k_max
        for x in 0:k
            lambda      = AnyHitConstraint(x, k)
            automaton   = build_automaton(lambda)
            for n in k_max:N
                seq         = random_sequence(automaton, n)
                @test is_satisfied(lambda, seq)
            end # for

            n_max       = max(n_max, length(automaton))
        end # for
    end # for
    println("Biggest automaton for AnyHitConstraint: $n_max")
end # testset
@testset "Random sequence tests - RowMissConstraint" begin
    n_max = 0
    for k in 1:k_max
        for x in 0:k
            lambda      = RowMissConstraint(x, k)
            automaton   = build_automaton(lambda)
            for n in k_max:N
                seq         = random_sequence(automaton, n)
                @test is_satisfied(lambda, seq)
            end # for

            n_max       = max(n_max, length(automaton))
        end # for
    end # for
    println("Biggest automaton for RowMissConstraint: $n_max")
end # testset
@testset "Random sequence tests - RowHitConstraint" begin
    n_max = 0
    for k in 1:k_max
        for x in 0:k
            lambda      = RowHitConstraint(x, k)
            automaton   = build_automaton(lambda)
            for n in k_max:N
                seq         = random_sequence(automaton, n)
                @test is_satisfied(lambda, seq)
            end # for

            n_max       = max(n_max, length(automaton))
        end # for
    end # for
    println("Biggest automaton for RowHitConstraint: $n_max")
end # testset

end # module


text:
#!/usr/bin/env julia
using Base.Test
using RdRand

function bench_rdrand(N)
    for i=1:N
        rdrand()
    end
end

function bench_rand(N)
    for i=1:N
        rand()
    end
end

#warmup
bench_rdrand(1)
bench_rand(1)

N = 10_000_000
t1 = @elapsed bench_rdrand(N)
t2 = @elapsed bench_rand(N)

@printf("%i rdrand() calls took %.3f seconds\n", N, t1)
@printf("%.3f nanoseconds per call\n", t1/N*1e9)

@printf("%i rand() calls took %.3f seconds\n", N, t2)
@printf("%.3f nanoseconds per call\n", t2/N*1e9)


text:
using Base.Test
using PropertyGraph
using UUID
using Compat

function graphpopulationtest()
	# create a graph
	g = Graph()
	g.tracker = ChangeTracker()

	@test g.tracker.ischanged == false

	nodea = add!(g, Vertex("Node", @compat Dict{String,Any}("name"=>"a", "val"=>50)))
	nodeb = add!(g, Vertex("Node", @compat Dict{String,Any}("name"=>"b", "val"=>51)))
	nodec = add!(g, Vertex("Node", @compat Dict{String,Any}("name"=>"c", "val"=>60)))
	noded = add!(g, Vertex("Node", @compat Dict{String,Any}("name"=>"d", "val"=>70)))

	edgeab = add!(g, Edge("NodeToNode", nodea, nodeb, @compat Dict{String,Any}("val"=>100)))
	edgeac = add!(g, Edge("NodeToNode", nodea, nodec, @compat Dict{String,Any}("val"=>160)))

	@test g.tracker.ischanged == true
	@test haskey(g.tracker.newvertices, nodea.id)
	@test haskey(g.tracker.newvertices, nodeb.id)
	@test haskey(g.tracker.newvertices, nodec.id)
	@test haskey(g.tracker.newvertices, noded.id)
	@test haskey(g.tracker.newedges, edgeab.id)
	@test haskey(g.tracker.newedges, edgeac.id)

	@test length(g.tracker.newvertices) == 4
	@test length(g.tracker.newedges) == 2
	@test length(g.tracker.changededges) == 0
	@test length(g.tracker.changedvertices) == 0
	@test length(g.tracker.changedgraphs) == 0
	@test length(g.tracker.deletededges) == 0
	@test length(g.tracker.deletedvertices) == 0

	setpropertyvalues!(nodea, @compat Dict{String,Any}("val"=>150,"new"=>"new value"))
	setpropertyvalue!(edgeab, "val", 200)

	@test haskey(g.tracker.changedvertices, nodea.id)
	@test haskey(g.tracker.changededges, edgeab.id)
	@test length(g.tracker.newvertices) == 4
	@test length(g.tracker.newedges) == 2
	@test length(g.tracker.changededges) == 1
	@test length(g.tracker.changedvertices) == 1
	@test length(g.tracker.changedgraphs) == 0
	@test length(g.tracker.deletededges) == 0
	@test length(g.tracker.deletedvertices) == 0

	clearchanges(g.tracker)

	@test g.tracker.ischanged == false
	@test length(g.tracker.newvertices) == 0
	@test length(g.tracker.newedges) == 0
	@test length(g.tracker.changededges) == 0
	@test length(g.tracker.changedvertices) == 0
	@test length(g.tracker.changedgraphs) == 0
	@test length(g.tracker.deletededges) == 0
	@test length(g.tracker.deletedvertices) == 0

	remove!(g, noded)

	@test g.tracker.ischanged == true
	@test haskey(g.tracker.deletedvertices, noded.id)
	@test length(g.tracker.newvertices) == 0
	@test length(g.tracker.newedges) == 0
	@test length(g.tracker.changededges) == 0
	@test length(g.tracker.changedvertices) == 0
	@test length(g.tracker.changedgraphs) == 0
	@test length(g.tracker.deletededges) == 0
	@test length(g.tracker.deletedvertices) == 1

	remove!(g, edgeab)

	@test haskey(g.tracker.deletededges, edgeab.id)
	@test length(g.tracker.newvertices) == 0
	@test length(g.tracker.newedges) == 0
	@test length(g.tracker.changededges) == 0
	@test length(g.tracker.changedvertices) == 0
	@test length(g.tracker.changedgraphs) == 0
	@test length(g.tracker.deletededges) == 1
	@test length(g.tracker.deletedvertices) == 1

	remove!(g, nodec) # removing a node will auto-remove its edges

	@test haskey(g.tracker.deletedvertices, nodec.id)
	@test haskey(g.tracker.deletededges, edgeac.id)
	@test length(g.tracker.newvertices) == 0
	@test length(g.tracker.newedges) == 0
	@test length(g.tracker.changededges) == 0
	@test length(g.tracker.changedvertices) == 0
	@test length(g.tracker.changedgraphs) == 0
	@test length(g.tracker.deletededges) == 2
	@test length(g.tracker.deletedvertices) == 2

	clearchanges(g.tracker)
end


text:
using Documenter, ABC

makedocs(
    format = :html,
    sitename = "ABC.jl",
    modules = [ABC],
#    doctest = true,
#    clean = false,
    pages = [
        "index.md", "page1.md",
    ],
    html_prettyurls = !("local" in ARGS)

)

deploydocs(
    repo   = "github.com/eford/ABC.jl.git",
    julia  = "0.6",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
import Ecto.English: pluralize, singularize
using Base.Test

@test "users" == pluralize("User")
@test "user" == singularize("users")


text:
function convertodays(x)
    if ismissing(x)
        return missing
    else
        return Day(x)
    end
end


text:
"""
    abstract type AbstractLine{T, N, S} <: AbstractGeometry{T, N} end

T: eltype
N: Dimension
S: Coordinate System: `:cartesian` or `:cylindrical`
"""
abstract type AbstractLine{T, N, S} <: AbstractGeometry{T, N} end

"""
    struct Line{T, N, S} <: AbstractLine{T, N, S}

<----A----------B---->
"""
struct Line{T, N, S} <: AbstractLine{T, N, S}
    org::AbstractCoordinatePoint{T, N, S}
    dir::AbstractCoordinateVector{T, N, S}
end

function Line(p1::CartesianPoint{T}, p2::CartesianPoint{T})::Line{T, 3, :cartesian} where {T}
    return Line{T, 3, :cartesian}(p1, p2 - p1)
end
function in(pt::CartesianPoint{T}, l::Line{T, 3, :cartesian}; atol::Real = T(0))::Bool where {T}
    norm((pt - l.org) × l.dir) <= atol
end

"""
    struct Ray{T, N, S} <: AbstractLine{T, N, S}

    [A----------B---->
"""
struct Ray{T, N, S} <: AbstractLine{T, N, S}
    org::AbstractCoordinatePoint{T, N, S}
    dir::AbstractCoordinateVector{T, N, S}
end
function Ray(p1::CartesianPoint{T}, p2::CartesianPoint{T})::Ray{T, 3, :cartesian} where {T}
    return Ray{T, 3, :cartesian}(p1, p2 - p1)
end
function in(pt::CartesianPoint{T}, l::Ray{T, 3, :cartesian}; atol::Real = T(0))::Bool where {T}
    d::CartesianVector{T} = pt - l.org
    if norm(d × l.dir) <= atol
        return d ⋅ l.dir >= 0
    end
    return false
end

"""
    struct LineSegment{T, N, S} <: AbstractLine{T, N, S}

    [A----------B]
"""
struct LineSegment{T, N, S} <: AbstractLine{T, N, S}
    org::AbstractCoordinatePoint{T, N, S}
    dir::AbstractCoordinateVector{T, N, S}
end

function LineSegment(p1::CartesianPoint{T}, p2::CartesianPoint{T})::LineSegment{T, 3, :cartesian} where {T}
    return LineSegment{T, 3, :cartesian}(p1, p2 - p1)
end
function in(pt::CartesianPoint{T}, l::LineSegment{T, 3, :cartesian}; atol::T = T(0))::Bool where {T}
    d::CartesianVector{T} = pt - l.org
    if norm(d × l.dir) <= atol
        d2::CartesianVector{T} = d - l.dir
        return (d ⋅ l.dir >= 0) && (d2 ⋅ l.dir <= 0)
    end
    return false
end

function iscolinear(l1::AbstractLine{T, 3, :cartesian}, l2::AbstractLine{T, 3, :cartesian}, atol::Real = T(0))::Bool where {T}
    return norm(l1.dir × l2.dir) <= atol # l1.dir == l2.dir || norm(l1.dir × l2.dir) == 0
end

# Plot Recipes:
@recipe function f(l::AbstractLine{T, 3, :cartesian}) where {T}
    x::Vector{T} = [l.org.x, l.org.x + l.dir.x]
    y::Vector{T} = [l.org.y, l.org.y + l.dir.y]
    z::Vector{T} = [l.org.z, l.org.z + l.dir.z]
    x, y, z
end

@recipe function f(l::AbstractLine{T, 2, :cartesian}) where {T}
    x::Vector{T} = [l.org.x, l.org.x + l.dir.x]
    y::Vector{T} = [l.org.y, l.org.y + l.dir.y]
    x, y
end

@recipe function f(ls::Array{<:AbstractLine{T, 3, :cartesian}, 1}) where {T}
    linecolor --> 1
    for l in ls
        @series begin
            label := ""
            l
        end
    end
    @series begin
        label --> ""
        [], []
    end
end