text:
module Remez

using Compat

export ratfn_minimax

"""
Evaluate a polynomial.

Arguments:
   coeffs   Array of BigFloats giving the coefficients of the polynomial.
            Starts with the constant term, i.e. coeffs[i] is the
            coefficient of x^(i-1) (because Julia arrays are 1-based).
   x        Point at which to evaluate the polynomial.

Return value is a BigFloat.
"""
function poly_eval(coeffs::Array{BigFloat}, x::BigFloat)
    n = length(coeffs)
    if n == 0
        return BigFloat(0)
    end
    y = coeffs[n]
    for i = n-1:-1:1
        y = coeffs[i] + x*y
    end
    y
end

"""
Evaluate a rational function.

Arguments:
   ncoeffs  Array of BigFloats giving the coefficients of the numerator.
            Starts with the constant term, and 1-based, as above.
   dcoeffs  Array of BigFloats giving the coefficients of the denominator.
            Starts with the constant term, and 1-based, as above.
   x        Point at which to evaluate the function.

Return value is a BigFloat.
"""
function ratfn_eval(ncoeffs::Array{BigFloat}, dcoeffs::Array{BigFloat},
                    x::BigFloat)
    return poly_eval(ncoeffs, x) / poly_eval(dcoeffs, x)
end

# ----------------------------------------------------------------------
# Least-squares fitting of a rational function to a set of (x,y)
# points.
#
# We use this to get an initial starting point for the Remez
# iteration. Therefore, it doesn't really need to be particularly
# accurate; it only needs to be good enough to wiggle back and forth
# across the target function the right number of times (so as to give
# enough error extrema to start optimising from) and not have any
# poles in the target interval.
#
# Least-squares fitting of a _polynomial_ is actually a sensible thing
# to do, and minimises the rms error. Doing the following trick with a
# rational function P/Q is less sensible, because it cannot be made to
# minimise the error function (P/Q-f)^2 that you actually wanted;
# instead it minimises (P-fQ)^2. But that should be good enough to
# have the properties described above.
#
# Some theory: suppose you're trying to choose a set of parameters a_i
# so as to minimise the sum of squares of some error function E_i.
# Basic calculus says, if you do this in one variable, just
# differentiate and solve for zero. In this case, that works fine even
# with multiple variables, because you _partially_ differentiate with
# respect to each a_i, giving a system of equations, and that system
# turns out to be linear so we just solve it as a matrix.
#
# In this case, our parameters are the coefficients of P and Q; to
# avoid underdetermining the system we'll fix Q's constant term at 1,
# so that our error function (as described above) is
#
# E = \sum (p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d)^2
#
# where the sum is over all (x,y) coordinate pairs. Setting dE/dp_j=0
# (for each j) gives an equation of the form
#
# 0 = \sum 2(p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d) x^j
#
# and setting dE/dq_j=0 gives one of the form
#
# 0 = \sum 2(p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d) y x^j
#
# And both of those row types, treated as multivariate linear
# equations in the p,q values, have each coefficient being a value of
# the form \sum x^i, \sum y x^i or \sum y^2 x^i, for various i. (Times
# a factor of 2, but we can throw that away.) So we can go through the
# list of input coordinates summing all of those things, and then we
# have enough information to construct our matrix and solve it
# straight off for the rational function coefficients.

# Arguments:
#    f        The function to be approximated. Maps BigFloat -> BigFloat.
#    xvals    Array of BigFloats, giving the list of x-coordinates at which
#             to evaluate f.
#    n        Degree of the numerator polynomial of the desired rational
#             function.
#    d        Degree of the denominator polynomial of the desired rational
#             function.
#    w        Error-weighting function. Takes two BigFloat arguments x,y
#             and returns a scaling factor for the error at that location.
#             A larger value indicates that the error should be given
#             greater weight in the square sum we try to minimise.
#             If unspecified, defaults to giving everything the same weight.
#
# Return values: a pair of arrays of BigFloats (N,D) giving the
# coefficients of the returned rational function. N has size n+1; D
# has size d+1. Both start with the constant term, i.e. N[i] is the
# coefficient of x^(i-1) (because Julia arrays are 1-based). D[1] will
# be 1.
function ratfn_leastsquares(f::Function, xvals::Array{BigFloat}, n, d,
                            w = (x,y)->BigFloat(1))
    # Accumulate sums of x^i y^j, for j={0,1,2} and a range of x.
    # Again because Julia arrays are 1-based, we'll have sums[i,j]
    # being the sum of x^(i-1) y^(j-1).
    maxpow = max(n,d) * 2 + 1
    sums = zeros(BigFloat, maxpow, 3)
    for x = xvals
        y = f(x)
        weight = w(x,y)
        for i = 1:1:maxpow
            for j = 1:1:3
                sums[i,j] += x^(i-1) * y^(j-1) * weight
            end
        end
    end

    # Build the matrix. We're solving n+d+1 equations in n+d+1
    # unknowns. (We actually have to return n+d+2 coefficients, but
    # one of them is hardwired to 1.)
    matrix = Array{BigFloat}(undef, n+d+1, n+d+1)
    vector = Array{BigFloat}(undef, n+d+1)
    for i = 0:1:n
        # Equation obtained by differentiating with respect to p_i,
        # i.e. the numerator coefficient of x^i.
        row = 1+i
        for j = 0:1:n
            matrix[row, 1+j] = sums[1+i+j, 1]
        end
        for j = 1:1:d
            matrix[row, 1+n+j] = -sums[1+i+j, 2]
        end
        vector[row] = sums[1+i, 2]
    end
    for i = 1:1:d
        # Equation obtained by differentiating with respect to q_i,
        # i.e. the denominator coefficient of x^i.
        row = 1+n+i
        for j = 0:1:n
            matrix[row, 1+j] = sums[1+i+j, 2]
        end
        for j = 1:1:d
            matrix[row, 1+n+j] = -sums[1+i+j, 3]
        end
        vector[row] = sums[1+i, 3]
    end

    # Solve the matrix equation.
    all_coeffs = matrix \ vector

    # And marshal the results into two separate polynomial vectors to
    # return.
    ncoeffs = all_coeffs[1:n+1]
    dcoeffs = vcat([1], all_coeffs[n+2:n+d+1])
    return (ncoeffs, dcoeffs)
end

# ----------------------------------------------------------------------
# Golden-section search to find a maximum of a function.

# Arguments:
#    f        Function to be maximised/minimised. Maps BigFloat -> BigFloat.
#    a,b,c    BigFloats bracketing a maximum of the function.
#
# Expects:
#    a,b,c are in order (either a<=b<=c or c<=b<=a)
#    a != c             (but b can equal one or the other if it wants to)
#    f(a) <= f(b) >= f(c)
#
# Return value is an (x,y) pair of BigFloats giving the extremal input
# and output. (That is, y=f(x).)
function goldensection(f::Function, a::BigFloat, b::BigFloat, c::BigFloat)
    # Decide on a 'good enough' threshold.
    epsbits = precision(BigFloat)
    threshold = abs(c-a) * 2^(-epsbits/2)

    # We'll need the golden ratio phi, of course. Or rather, in this
    # case, we need 1/phi = 0.618...
    one_over_phi = 2 / (1 + sqrt(BigFloat(5)))

    # Flip round the interval endpoints so that the interval [a,b] is
    # at least as large as [b,c]. (Then we can always pick our new
    # point in [a,b] without having to handle lots of special cases.)
    if abs(b-a) < abs(c-a)
        a,  c  = c,  a
    end

    # Evaluate the function at the initial points.
    fa = f(a)
    fb = f(b)
    fc = f(c)

    while abs(c-a) > threshold

        # Check invariants.
        @assert(a <= b <= c || c <= b <= a)
        @assert(fa <= fb >= fc)

        # Subdivide the larger of the intervals [a,b] and [b,c]. We've
        # arranged that this is always [a,b], for simplicity.
        d = a + (b-a) * one_over_phi

        # Now we have an interval looking like this (possibly
        # reversed):
        #
        #    a            d       b            c
        #
        # and we know f(b) is bigger than either f(a) or f(c). We have
        # two cases: either f(d) > f(b), or vice versa. In either
        # case, we can narrow to an interval of 1/phi the size, and
        # still satisfy all our invariants (three ordered points,
        # [a,b] at least the width of [b,c], f(a)<=f(b)>=f(c)).
        fd = f(d)
        if fd > fb
            a,  b,  c  = a,  d,  b
            fa, fb, fc = fa, fd, fb
        else
            a,  b,  c  = c,  b,  d
            fa, fb, fc = fc, fb, fd
        end
    end

    return (b, fb)
end

# ----------------------------------------------------------------------
# Find the extrema of a function within a given interval.

# Arguments:
#    f         The function to be approximated. Maps BigFloat -> BigFloat.
#    grid      A set of points at which to evaluate f. Must be high enough
#              resolution to make extrema obvious.
#
# Returns an array of (x,y) pairs of BigFloats, with each x,y giving
# the extremum location and its value (i.e. y=f(x)).
function find_extrema(f::Function, grid::Array{BigFloat})
    len = length(grid)
    extrema = Tuple{BigFloat, BigFloat}[]
    for i = 1:1:len
        # We have to provide goldensection() with three points
        # bracketing the extremum. If the extremum is at one end of
        # the interval, then the only way we can do that is to set two
        # of the points equal (which goldensection() will cope with).
        prev = max(1, i-1)
        next = min(i+1, len)

        # Find our three pairs of (x,y) coordinates.
        xp, xi, xn = grid[prev], grid[i], grid[next]
        yp, yi, yn = f(xp), f(xi), f(xn)

        # See if they look like an extremum, and if so, ask
        # goldensection() to give a more exact location for it.
        if yp <= yi >= yn
            push!(extrema, goldensection(f, xp, xi, xn))
        elseif yp >= yi <= yn
            x, y = goldensection(x->-f(x), xp, xi, xn)
            push!(extrema, (x, -y))
        end
    end
    return extrema
end

# ----------------------------------------------------------------------
# Winnow a list of a function's extrema to give a subsequence of a
# specified length, with the extrema in the subsequence alternating
# signs, and with the smallest absolute value of an extremum in the
# subsequence as large as possible.
#
# We do this using a dynamic-programming approach. We work along the
# provided array of extrema, and at all times, we track the best set
# of extrema we have so far seen for each possible (length, sign of
# last extremum) pair. Each new extremum is evaluated to see whether
# it can be added to any previously seen best subsequence to make a
# new subsequence that beats the previous record holder in its slot.

# Arguments:
#    extrema   An array of (x,y) pairs of BigFloats giving the input extrema.
#    n         Number of extrema required as output.
#
# Returns a new array of (x,y) pairs which is a subsequence of the
# original sequence. (So, in particular, if the input was sorted by x
# then so will the output be.)
function winnow_extrema(extrema::Array{Tuple{BigFloat,BigFloat}}, n)
    # best[i,j] gives the best sequence so far of length i and with
    # sign j (where signs are coded as 1=positive, 2=negative), in the
    # form of a tuple (cost, actual array of x,y pairs).
    best = fill((BigFloat(0), Tuple{BigFloat,BigFloat}[]), n, 2)

    for (x,y) = extrema
        if y > 0
            sign = 1
        elseif y < 0
            sign = 2
        else
            # A zero-valued extremum cannot possibly contribute to any
            # optimal sequence, so we simply ignore it!
            continue
        end

        for i = 1:1:n
            # See if we can create a new entry for best[i,sign] by
            # appending our current (x,y) to some previous thing.
            if i == 1
                # Special case: we don't store a best zero-length
                # sequence :-)
                candidate = (abs(y), [(x,y)])
            else
                othersign = 3-sign # map 1->2 and 2->1
                oldscore, oldlist = best[i-1, othersign]
                newscore = min(abs(y), oldscore)
                newlist = vcat(oldlist, [(x,y)])
                candidate = (newscore, newlist)
            end
            # If our new candidate improves on the previous value of
            # best[i,sign], then replace it.
            if candidate[1] > best[i,sign][1]
                best[i,sign] = candidate
            end
        end
    end

    # Our ultimate return value has to be either best[n,1] or
    # best[n,2], but it could be either. See which one has the higher
    # score.
    if best[n,1][1] > best[n,2][1]
        ret = best[n,1][2]
    else
        ret = best[n,2][2]
    end
    # Make sure we did actually _find_ a good answer.
    @assert(length(ret) == n)
    return ret
end

# ----------------------------------------------------------------------
# Construct a rational-function approximation with equal and
# alternating weighted deviation at a specific set of x-coordinates.

# Arguments:
#    f         The function to be approximated. Maps BigFloat -> BigFloat.
#    coords    An array of BigFloats giving the x-coordinates. There should
#              be n+d+2 of them.
#    n, d      The degrees of the numerator and denominator of the desired
#              approximation.
#    prev_err  A plausible value for the alternating weighted deviation.
#              (Required to kickstart a binary search in the nonlinear case;
#              see comments below.)
#    w         Error-weighting function. Takes two BigFloat arguments x,y
#              and returns a scaling factor for the error at that location.
#              The returned approximation R should have the minimum possible
#              maximum value of abs((f(x)-R(x)) * w(x,f(x))). Optional
#              parameter, defaulting to the always-return-1 function.
#
# Return values: a pair of arrays of BigFloats (N,D) giving the
# coefficients of the returned rational function. N has size n+1; D
# has size d+1. Both start with the constant term, i.e. N[i] is the
# coefficient of x^(i-1) (because Julia arrays are 1-based). D[1] will
# be 1.
function ratfn_equal_deviation(f::Function, coords::Array{BigFloat},
                               n, d, prev_err::BigFloat,
                               w = (x,y)->BigFloat(1))
    @assert(length(coords) == n+d+2)

    if d == 0
        # Special case: we're after a polynomial. In this case, we
        # have the particularly easy job of just constructing and
        # solving a system of n+2 linear equations, to find the n+1
        # coefficients of the polynomial and also the amount of
        # deviation at the specified coordinates. Each equation is of
        # the form
        #
        #   p_0 x^0 + p_1 x^1 + ... + p_n x^n ± e/w(x) = f(x)
        #
        # in which the p_i and e are the variables, and the powers of
        # x and calls to w and f are the coefficients.

        matrix = Array{BigFloat}(undef, n+2, n+2)
        vector = Array{BigFloat}(undef, n+2)
        currsign = +1
        for i = 1:1:n+2
            x = coords[i]
            for j = 0:1:n
                matrix[i,1+j] = x^j
            end
            y = f(x)
            vector[i] = y
            matrix[i, n+2] = currsign / w(x,y)
            currsign = -currsign
        end

        outvector = matrix \ vector

        ncoeffs = outvector[1:n+1]
        dcoeffs = [BigFloat(1)]
        return ncoeffs, dcoeffs
    else
        # For a nontrivial rational function, the system of equations
        # we need to solve becomes nonlinear, because each equation
        # now takes the form
        #
        #   p_0 x^0 + p_1 x^1 + ... + p_n x^n
        #   --------------------------------- ± e/w(x) = f(x)
        #     x^0 + q_1 x^1 + ... + q_d x^d
        #
        # and multiplying up by the denominator gives you a lot of
        # terms containing e × q_i. So we can't do this the really
        # easy way using a matrix equation as above.
        #
        # Fortunately, this is a fairly easy kind of nonlinear system.
        # The equations all become linear if you switch to treating e
        # as a constant, so a reasonably sensible approach is to pick
        # a candidate value of e, solve all but one of the equations
        # for the remaining unknowns, and then see what the error
        # turns out to be in the final equation. The Chebyshev
        # alternation theorem guarantees that that error in the last
        # equation will be anti-monotonic in the input e, so we can
        # just binary-search until we get the two as close to equal as
        # we need them.

        function try_e(e)
            # Try a given value of e, derive the coefficients of the
            # resulting rational function by setting up equations
            # based on the first n+d+1 of the n+d+2 coordinates, and
            # see what the error turns out to be at the final
            # coordinate.
            matrix = Array{BigFloat}(undef, n+d+1, n+d+1)
            vector = Array{BigFloat}(undef, n+d+1)
            currsign = +1
            for i = 1:1:n+d+1
                x = coords[i]
                y = f(x)
                y_adj = y - currsign * e / w(x,y)
                for j = 0:1:n
                    matrix[i,1+j] = x^j
                end
                for j = 1:1:d
                    matrix[i,1+n+j] = -x^j * y_adj
                end
                vector[i] = y_adj
                currsign = -currsign
            end

            outvector = matrix \ vector

            ncoeffs = outvector[1:n+1]
            dcoeffs = vcat([BigFloat(1)], outvector[n+2:n+d+1])

            x = coords[n+d+2]
            y = f(x)
            last_e = (ratfn_eval(ncoeffs, dcoeffs, x) - y) * w(x,y) * -currsign

            return ncoeffs, dcoeffs, last_e
        end

        epsbits = precision(BigFloat)
        threshold = 2^(-epsbits/2) # convergence threshold

        # Start by trying our previous iteration's error value. This
        # value (e0) will be one end of our binary-search interval,
        # and whatever it caused the last point's error to be, that
        # (e1) will be the other end.
        e0 = prev_err
        nc, dc, e1 = try_e(e0)
        if abs(e1-e0) <= threshold
            # If we're _really_ lucky, we hit the error right on the
            # nose just by doing that!
            return nc, dc
        end
        s = sign(e1-e0)

        # Verify by assertion that trying our other interval endpoint
        # e1 gives a value that's wrong in the other direction.
        # (Otherwise our binary search won't get a sensible answer at
        # all.)
        nc, dc, e2 = try_e(e1)
        @assert(sign(e2-e1) == -s)

        # Now binary-search until our two endpoints narrow enough.
        local emid
        while abs(e1-e0) > threshold
            emid = (e1+e0)/2
            nc, dc, enew = try_e(emid)
            if sign(enew-emid) == s
                e0 = emid
            else
                e1 = emid
            end
        end

        return nc, dc
    end
end

"""
    N,D,E,X = ratfn_minimax(f, interval, n, d, w)

Top-level function to find a minimax rational-function approximation.

Arguments:

 * f         The function to be approximated. Maps BigFloat -> BigFloat.
 * interval  A tuple giving the endpoints of the interval
             (in either order) on which to approximate f.
 * n, d      The degrees of the numerator and denominator of the desired
             approximation.
 * w         Error-weighting function. Takes two BigFloat arguments x,y
             and returns a scaling factor for the error at that location.
             The returned approximation R should have the minimum possible
             maximum value of abs((f(x)-R(x)) * w(x,f(x))). Optional
             parameter, defaulting to the always-return-1 function.

Return values: a tuple (N,D,E,X), where

 * N,D       A pair of arrays of BigFloats giving the coefficients
             of the returned rational function. N has size n+1; D
             has size d+1. Both start with the constant term, i.e.
             N[i] is the coefficient of x^(i-1) (because Julia
             arrays are 1-based). D[1] will be 1.
 * E         The maximum weighted error (BigFloat).
 * X         An array of pairs of BigFloats giving the locations of n+2
             points and the weighted error at each of those points. The
             weighted error values will have alternating signs, which
             means that the Chebyshev alternation theorem guarantees
             that any other function of the same degree must exceed
             the error of this one at at least one of those points.
"""
function ratfn_minimax(f, interval, n, d,
                       w = (x,y)->BigFloat(1))
    # We start off by finding a least-squares approximation. This
    # doesn't need to be perfect, but if we can get it reasonably good
    # then it'll save iterations in the refining stage.
    #
    # Least-squares approximations tend to look nicer in a minimax
    # sense if you evaluate the function at a big pile of Chebyshev
    # nodes rather than uniformly spaced points. These values will
    # also make a good grid to use for the initial search for error
    # extrema, so we'll keep them around for that reason too.

    # Construct the grid.
    lo = BigFloat(minimum(interval))
    hi = BigFloat(maximum(interval))

    local grid
    let
        mid = (hi+lo)/2
        halfwid = (hi-lo)/2
        nnodes = 16 * (n+d+1)
        grid = [ mid - halfwid * cospi(big(i)/big(nnodes)) for i=0:nnodes ]
    end

    # Find the initial least-squares approximation.
    (nc, dc) = ratfn_leastsquares(f, grid, n, d, w)

    # Threshold of convergence. We stop when the relative difference
    # between the min and max (winnowed) error extrema is less than
    # this.
    #
    # This is set to the cube root of machine epsilon on a more or
    # less empirical basis, because the rational-function case will
    # not converge reliably if you set it to only the square root.
    # (Repeatable by using the --test mode.) On the assumption that
    # input and output error in each iteration can be expected to be
    # related by a simple power law (because it'll just be down to how
    # many leading terms of a Taylor series are zero), the cube root
    # was the next thing to try.
    epsbits = precision(BigFloat)
    threshold = 2^(-epsbits/3)

    # Main loop.
    while true
        # Find all the error extrema we can.
        function compute_error(x)
            real_y = f(x)
            approx_y = ratfn_eval(nc, dc, x)
            return (approx_y - real_y) * w(x, real_y)
        end
        extrema = find_extrema(compute_error, grid)

        # Winnow the extrema down to the right number, and ensure they
        # have alternating sign.
        extrema = winnow_extrema(extrema, n+d+2)

        # See if we've finished.
        min_err = minimum([abs(y) for (x,y) = extrema])
        max_err = maximum([abs(y) for (x,y) = extrema])
        variation = (max_err - min_err) / max_err
        if variation < threshold
            return nc, dc, max_err, extrema
        end

        # If not, refine our function by equalising the error at the
        # extrema points, and go round again.
        (nc, dc) = ratfn_equal_deviation(f, map(x->x[1], extrema),
                                         n, d, max_err, w)
    end
end

# ----------------------------------------------------------------------
# Check if a polynomial is well-conditioned for accurate evaluation in
# a given interval by Horner's rule.
#
# This is true if at every step where Horner's rule computes
# (coefficient + x*value_so_far), the constant coefficient you're
# adding on is of larger magnitude than the x*value_so_far operand.
# And this has to be true for every x in the interval.
#
# Arguments:
#    coeffs    The coefficients of the polynomial under test. Starts with
#              the constant term, i.e. coeffs[i] is the coefficient of
#              x^(i-1) (because Julia arrays are 1-based).
#    lo, hi    The bounds of the interval.
#
# Return value: the largest ratio (x*value_so_far / coefficient), at
# any step of evaluation, for any x in the interval. If this is less
# than 1, the polynomial is at least somewhat well-conditioned;
# ideally you want it to be more like 1/8 or 1/16 or so, so that the
# relative rounding error accumulated at each step are reduced by
# several factors of 2 when the next coefficient is added on.

function wellcond(coeffs, lo, hi)
    x = max(abs(lo), abs(hi))
    worst = 0
    so_far = 0
    for i = length(coeffs):-1:1
        coeff = abs(coeffs[i])
        so_far *= x
        if coeff != 0
            thisval = so_far / coeff
            worst = max(worst, thisval)
            so_far += coeff
        end
    end
    return worst
end

end # module


text:
function _bounds(id::I,
                 new_state::DiseaseState,
                 extents::EventExtents{T},
                 obs::EventObservations{T},
                 events::Events{T},
                 network::TransmissionNetwork) where {I <: Integer, T <: EpidemicModel}
  if new_state == State_E
    lowerbounds = [events.infection[id] - extents.exposure[2]]
    upperbounds = [events.infection[id] - extents.exposure[1]]
    path_to = _pathway_to(id, network, depth = 1)
    if length(path_to) > 1
      parent_host = path_to[2]
      push!(lowerbounds, events.infection[parent_host])
      if (T == SEIR) && !isnan(events.removal[parent_host])
        push!(upperbounds, events.removal[parent_host])
      end
    end
  elseif new_state == State_I
    path_from = _pathway_from(id, network, depth = 1)
    if T in [SEIR; SEI]
      lowerbounds = [obs.infection[id] - extents.infection[2]
                     events.exposure[id] + extents.exposure[1]]
      upperbounds = [obs.infection[id] - extents.infection[1]
                     events.exposure[id] + extents.exposure[2]]
      if length(path_from) > 1
        child_hosts = path_from[2:end]
        append!(upperbounds, events.exposure[child_hosts])
      end
    elseif T in [SIR; SI]
      lowerbounds = [obs.infection[id] - extents.infection[2]]
      upperbounds = [obs.infection[id] - extents.infection[1]]
      path_to = _pathway_to(id, network, depth = 1)
      if length(path_from) > 1
        child_hosts = path_from[2:end]
        append!(upperbounds, events.infection[child_hosts])
      end
      if length(path_to) > 1
        parent_host = path_to[2]
        push!(lowerbounds, events.infection[parent_host])
        if (T == SIR) && !isnan(events.removal[parent_host])
          push!(upperbounds, events.removal[parent_host])
        end
      end
    end
  elseif new_state == State_R
    path_from = _pathway_from(id, network, depth = 1)
    lowerbounds = [obs.removal[id] - extents.removal[2]
                   obs.infection[id]]
    upperbounds = [obs.removal[id] - extents.removal[1]]
    if length(path_from) > 1
      child_hosts = path_from[2:end]
      if T == SEIR
        append!(lowerbounds, events.exposure[child_hosts])
      elseif T == SIR
        append!(lowerbounds, events.infection[child_hosts])
      end
    end
  end
  @debug "Transition of $id into $new_state with bounds: \n  [max($(round.(lowerbounds, digits=3))),\n   min($(round.(upperbounds, digits=3)))]"
  lowerbound = maximum(lowerbounds)
  upperbound = minimum(upperbounds)
  if lowerbound >= upperbound
    @warn "Lower >= upper bound for transition of i=$id into $new_state" lower = lowerbounds upper = upperbounds
  end
  return lowerbound, upperbound
end

function _bounds(last_event::Event{T},
                 extents::EventExtents{T},
                 obs::EventObservations{T},
                 events::Events{T},
                 network::TransmissionNetwork) where T <: EpidemicModel
  return _bounds(last_event.individual,
                 last_event.new_state,
                 extents, obs, events, network)
end

function _bounds(id::I,
                 new_state::DiseaseState,
                 extents::EventExtents{T},
                 obs::EventObservations{T},
                 events::Events{T}) where {I <: Integer, T <: EpidemicModel}
  if new_state == State_E
    lowerbounds = [events.infection[id] - extents.exposure[2]]
    upperbounds = [events.infection[id] - extents.exposure[1]]
  elseif new_state == State_I
    if T in [SEIR; SEI]
      lowerbounds = [obs.infection[id] - extents.infection[2]
                     events.exposure[id] + extents.exposure[1]]
      upperbounds = [obs.infection[id] - extents.infection[1]
                     events.exposure[id] + extents.exposure[2]]
    elseif T in [SIR; SI]
      lowerbounds = [obs.infection[id] - extents.infection[2]]
      upperbounds = [obs.infection[id] - extents.infection[1]]
    end
  elseif new_state == State_R
    lowerbounds = [obs.removal[id] - extents.removal[2]
                   obs.infection[id]]
    upperbounds = [obs.removal[id] - extents.removal[1]]
  end
  @debug "Transition of $id into $new_state with bounds: \n  [max($(round.(lowerbounds, digits=3))),\n   min($(round.(upperbounds, digits=3)))]"
  lowerbound = maximum(lowerbounds)
  upperbound = minimum(upperbounds)
  if lowerbound >= upperbound
    @warn "Lower >= upper bound for transition of i=$id into $new_state" lower = lowerbounds upper = upperbounds
  end
  return lowerbound, upperbound
end

function _bounds(last_event::Event{T},
                 extents::EventExtents{T},
                 obs::EventObservations{T},
                 events::Events{T}) where T <: EpidemicModel
  return _bounds(last_event.individual,
                 last_event.new_state,
                 extents, obs, events)
end

ext: jl


text:
@testset "alldifferent" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:10] <= 5, Int)
    @constraint(m, x in CS.AllDifferent())
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = get_constraints_by_type(com, CS.AllDifferentConstraint)[1]

    # doesn't check the length
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 2, 3])
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [2, 2, 3])

    sorted_min = [1, 1, 2, 2, 3]
    sorted_max = [5, 5, 4, 4, 2]
    @test CS.get_alldifferent_extrema(sorted_min, sorted_max, 3) == (1 + 2 + 3, 5 + 4 + 3)

    sorted_min = [1, 3, 3, 3, 4]
    sorted_max = [9, 7, 7, 7, 5]
    @test CS.get_alldifferent_extrema(sorted_min, sorted_max, 3) == (1 + 3 + 4, 9 + 7 + 6)
    @test CS.get_alldifferent_extrema(sorted_min, sorted_max, 5) ==
          (1 + 3 + 4 + 5 + 6, 9 + 7 + 6 + 5 + 4)

    constr_indices = constraint.indices
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[2],
        5,
    )
    @test CS.fix!(com, com.search_space[constr_indices[2]], 5)
    @test !CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        5,
    )

    # need to create a backtrack_vec to reverse pruning
    dummy_backtrack_obj = CS.BacktrackObj(com)
    dummy_backtrack_obj.step_nr = 1
    push!(com.backtrack_vec, dummy_backtrack_obj)
    # reverse previous fix
    CS.reverse_pruning!(com, 1)
    # now setting it to 5 should be feasible
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        5,
    )

    com.c_backtrack_idx = 1

    # feasible and no changes
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices
        @test sort(CS.values(com.search_space[ind])) == -5:5
    end
    @test CS.fix!(com, com.search_space[constr_indices[2]], 5)
    @test CS.rm!(com, com.search_space[constr_indices[1]], 1)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices[3:end]
        @test sort(CS.values(com.search_space[ind])) == -5:4
    end
    @test sort(CS.values(com.search_space[1])) == [-5, -4, -3, -2, -1, 0, 2, 3, 4]
    @test sort(CS.values(com.search_space[2])) == [5]

    # 3 and 4 are taken by indices 3 and 4 so not available at other positions
    @test CS.remove_below!(com, com.search_space[constr_indices[3]], 3)
    @test CS.remove_below!(com, com.search_space[constr_indices[4]], 3)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices[3:4]
        @test sort(CS.values(com.search_space[ind])) == 3:4
    end
    for ind in constr_indices[5:end]
        @test sort(CS.values(com.search_space[ind])) == -5:2
    end
    @test sort(CS.values(com.search_space[1])) == [-5, -4, -3, -2, -1, 0, 2]
    @test sort(CS.values(com.search_space[2])) == [5]

    # we don't need -5
    @test CS.remove_below!(com, com.search_space[constr_indices[1]], -4)
    for ind in constr_indices[5:end]
        @test CS.remove_below!(com, com.search_space[constr_indices[ind]], -4)
    end
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)

    # but we need -4 to have enough values available
    @test CS.remove_below!(com, com.search_space[constr_indices[1]], -3)
    for ind in constr_indices[5:end]
        @test CS.remove_below!(com, com.search_space[constr_indices[ind]], -3)
    end
    @test !CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
end

@testset "all different with gap in variables" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, x[1:4], CS.Integers([-5, -2, 3, 0, 7]))
    @constraint(m, x in CS.AllDifferent())
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = get_constraints_by_type(com, CS.AllDifferentConstraint)[1]
    @test CS.is_constraint_solved(
        constraint,
        constraint.fct,
        constraint.set,
        [-2, 0, 7, -5],
    )
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)

    constr_indices = constraint.indices
    for ind in constr_indices
        @test sort(CS.values(com.search_space[ind])) == [-5, -2, 0, 3, 7]
    end
    @test CS.fix!(com, com.search_space[constr_indices[1]], -5)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices[2:4]
        @test sort(CS.values(com.search_space[ind])) == [-2, 0, 3, 7]
    end
    @test CS.rm!(com, com.search_space[constr_indices[2]], -2)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    @test sort(CS.values(com.search_space[2])) == [0, 3, 7]
    for ind in constr_indices[3:4]
        @test sort(CS.values(com.search_space[ind])) == [-2, 0, 3, 7]
    end

    @test CS.remove_below!(com, com.search_space[constr_indices[3]], 2)
    @test CS.remove_below!(com, com.search_space[constr_indices[4]], 2)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    @test CS.isfixed(com.search_space[2])
    @test CS.value(com.search_space[2]) == 0
    for ind in constr_indices[3:4]
        @test sort(CS.values(com.search_space[ind])) == [3, 7]
    end
end

@testset "all different with huge gap in variables" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 0 <= x[1:4] <= 3, Int)
    @variable(m, 1000 <= y[1:4] <= 3000, Int)
    @constraint(m, [x...,y...] in CS.AllDifferent())
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = get_constraints_by_type(com, CS.AllDifferentConstraint)[1]
    constr_indices = constraint.indices
    @test CS.is_constraint_solved(
        constraint,
        constraint.fct,
        constraint.set,
        [0, 1, 2, 3, 1000, 2000],
    )
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)

    @test CS.fix!(com, com.search_space[constr_indices[5]], 1000)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices[6:8]
        @test sort(CS.values(com.search_space[ind])) == 1001:3000
    end
    @test CS.fix!(com, com.search_space[constr_indices[1]], 0)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices[2:4]
        @test sort(CS.values(com.search_space[ind])) == 1:3
    end
end

@testset "all different is_constraint_violated test" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:10] <= 5, Int)
    @constraint(m, x in CS.AllDifferent())
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = get_constraints_by_type(com, CS.AllDifferentConstraint)[1]

    variables = com.search_space
    @test CS.fix!(com, variables[1], 1; check_feasibility = false)
    @test CS.fix!(com, variables[2], 1; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:10] <= 5, Int)
    @constraint(m, x in CS.AllDifferent())
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = get_constraints_by_type(com, CS.AllDifferentConstraint)[1]

    variables = com.search_space
    @test CS.fix!(com, variables[1], 1; check_feasibility = false)
    @test CS.fix!(com, variables[2], 2; check_feasibility = false)
    @test !CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)
end

@testset "all 8queens solutions" begin
    n = 8
    model = Model(optimizer_with_attributes(CS.Optimizer, "all_optimal_solutions"=>true, "logging"=>[]))

    @variable(model, 1 <= x[1:n] <= n, Int)
    @constraint(model, x in CS.AllDifferent())
    @constraint(model, [x[i] + i for i in 1:n] in CS.AllDifferent())
    @constraint(model, [x[i] - i for i in 1:n] in CS.AllDifferent())

    optimize!(model)

    status = JuMP.termination_status(model)

    @test status == MOI.OPTIMAL
    num_sols = MOI.get(model, MOI.ResultCount())
    @test num_sols == 92
    for sol in 1:num_sols
        x_val = convert.(Integer,JuMP.value.(x; result=sol))
        @test allunique(x_val)
        @test allunique([x_val[i] + i for i in 1:n])
        @test allunique([x_val[i] - i for i in 1:n])
    end
end

@testset "domainerror due to term not discrete " begin
    n = 8
    model = Model(optimizer_with_attributes(CS.Optimizer, "all_optimal_solutions"=>true, "logging"=>[]))

    @variable(model, 1 <= x[1:n] <= n, Int)
    @constraint(model, x in CS.AllDifferent())
    # 1.5 not allowed
    @constraint(model, [1.5*x[i] + i for i in 1:n] in CS.AllDifferent())
    @constraint(model, [x[i] - i for i in 1:n] in CS.AllDifferent())

    @test_throws DomainError optimize!(model)
end


text:
using NBodyIPs, JuLIP, BenchmarkTools, StaticArrays
using Test

using NBodyIPs: BondLengthDesc, transform, transform_d, fcut, fcut_d,
                invariants, invariants_d, invariants_ed
using NBodyIPs.Polys: NBPoly

nbasis = [0, 10, 10, 10, 10]
nbasis = [0, 10, 50, 100, 300]

println("Setting up the test systems ...")
r0 = rnn(:Cu)
TRANSFORM = ExpTransform(3, r0)
at = rattle!(bulk(:Cu, cubic=true) * 2, 0.05)
rcut3 = 3.1 * r0
D3 = BondLengthDesc(TRANSFORM, CosCut(0.66*rcut3, rcut3) )
rcut4 = 2.1 * r0
D4 = BondLengthDesc(TRANSFORM, CosCut(0.66*rcut4, rcut4) )
rcut5 = 1.5 * r0
D5 = BondLengthDesc(TRANSFORM, CosCut(0.66*rcut5, rcut5) )
DD = [nothing, D3, D3, D4, D5]

random_nbody(N, ntup) = (
   (N <= 3) ? NBPoly( [tuple( [rand(1:4, ((N*(N-1))÷2)); 0]... ) for n = 1:ntup],
                     (0.1 .+ rand(ntup))/factorial(N), DD[N] )
           : NBPoly( [tuple( rand(0:N, ((N*(N-1))÷2+1))... ) for n = 1:ntup],
                   (0.1 .+ rand(ntup))/factorial(N), DD[N] ) )

BB = [ nothing; [ [random_nbody(N, 1) for _=1:nbasis[N]] for N = 2:5 ] ]
at = rattle!(bulk(:Cu, cubic=true) * 2, 0.05)

println("Performance Tests")
println("-----------------")
println("to get a rough estimate of the cost per atom take the `separate`")
println("timing and multiply by the factor provided")
for N = 2:5
   B = BB[N]
   println("[$N-body], length(B) == $(length(B)), factor = $(1/(length(at)*length(B)))")
   print(" E separate:" ); @btime ([energy(b, $at) for b in $B])
   print(" E combined:" ); @btime energy( $B, $at )
   print(" F separate:" ); @btime ([forces(b, $at) for b in $B])
   print(" F combined:" ); @btime forces( $B, $at )
   print(" S separate:" ); @btime ([virial(b, $at) for b in $B])
   print(" S combined:" ); @btime virial( $B, $at )
end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
module Radar

include("fmcw.jl")

end


text:
## UV based file operations ##

module FS

const S_IRUSR = 0o400
const S_IWUSR = 0o200
const S_IXUSR = 0o100
const S_IRWXU = 0o700
const S_IRGRP = 0o040
const S_IWGRP = 0o020
const S_IXGRP = 0o010
const S_IRWXG = 0o070
const S_IROTH = 0o004
const S_IWOTH = 0o002
const S_IXOTH = 0o001
const S_IRWXO = 0o007

export File,
       # open,
       # close,
       write,
       unlink,
       rename,
       sendfile,
       JL_O_WRONLY,
       JL_O_RDONLY,
       JL_O_RDWR,
       JL_O_APPEND,
       JL_O_CREAT,
       JL_O_EXCL,
       JL_O_TRUNC,
       JL_O_TEMPORARY,
       JL_O_SHORT_LIVED,
       JL_O_SEQUENTIAL,
       JL_O_RANDOM,
       JL_O_NOCTTY,
       S_IRUSR, S_IWUSR, S_IXUSR, S_IRWXU,
       S_IRGRP, S_IWGRP, S_IXGRP, S_IRWXG,
       S_IROTH, S_IWOTH, S_IXOTH, S_IRWXO

import Base: uvtype, uvhandle, eventloop, fd, position, stat, close, write, read, read!, readbytes, isopen,
            _sizeof_uv_fs, uv_error

include("file_constants.jl")

abstract AbstractFile <: IO

type File <: AbstractFile
    path::String
    open::Bool
    handle::Int32
    File(path::String) = new(path,false,-1)
    File(fd::RawFD) = new("",true,fd.fd)
end

type AsyncFile <: AbstractFile
    path::String
    open::Bool
end

isopen(f::Union(File,AsyncFile)) = f.open

uvtype(::File) = Base.UV_RAW_FD
uvhandle(file::File) = file.handle

_uv_fs_result(req) = ccall(:jl_uv_fs_result,Int32,(Ptr{Void},),req)

function open(f::File,flags::Integer,mode::Integer)
    req = c_malloc(_sizeof_uv_fs)
    ret = ccall(:uv_fs_open,Int32,(Ptr{Void},Ptr{Void},Ptr{Uint8},Int32,Int32,Ptr{Void}),
                eventloop(),req,bytestring(f.path),flags,mode,C_NULL)
    f.handle = _uv_fs_result(req)
    ccall(:uv_fs_req_cleanup,Void,(Ptr{Void},),req)
    c_free(req)
    uv_error("open",ret)
    f.open = true
    f
end
open(f::String,flags,mode) = open(File(f),flags,mode)
open(f::String,flags) = open(f,flags,0)

function close(f::File)
    if !f.open
        error("file is already closed")
    end
    err = ccall(:jl_fs_close, Int32, (Int32,), f.handle)
    uv_error("close",err)
    f.handle = -1
    f.open = false
    f
end

function unlink(p::String)
    err = ccall(:jl_fs_unlink, Int32, (Ptr{Uint8},), bytestring(p))
    uv_error("unlink",err)
end
function unlink(f::File)
    if isempty(f.path)
      error("no path associated with this file")
    end
    if f.open
        close(f)
    end
    unlink(f.path)
    f
end

# For move command
function rename(src::String, dst::String)
    err = ccall(:jl_fs_rename, Int32, (Ptr{Uint8}, Ptr{Uint8}), bytestring(src),
                bytestring(dst))

    # on error, default to cp && rm
    if err < 0
        # Note that those two functions already handle their errors.
        # first copy
        sendfile(src, dst)

        # then rm
        unlink(src)
    end
end

# For copy command
function sendfile(src::String, dst::String)
    flags = JL_O_RDONLY
    src_file = open(src, flags)
    if !src_file.open
        error("Src file is not open")
    end

    flags = JL_O_CREAT | JL_O_RDWR
    mode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP| S_IROTH | S_IWOTH
    dst_file = open(dst, flags, mode)
    if !dst_file.open
        error("Dst file is not open")
    end

    src_stat = stat(src_file)
    err = ccall(:jl_fs_sendfile, Int32, (Int32, Int32, Int64, Csize_t),
                fd(src_file), fd(dst_file), 0, src_stat.size)
    uv_error("sendfile", err)

    if src_file.open
        close(src_file)
    end
    if dst_file.open
        close(dst_file)
    end
end

function write(f::File, buf::Ptr{Uint8}, len::Integer, offset::Integer=-1)
    if !f.open
        error("file is not open")
    end
    err = ccall(:jl_fs_write, Int32, (Int32, Ptr{Uint8}, Csize_t, Csize_t),
                f.handle, buf, len, offset)
    uv_error("write",err)
    len
end

function write(f::File, c::Uint8)
    if !f.open
        error("file is not open")
    end
    err = ccall(:jl_fs_write_byte, Int32, (Int32, Cchar), f.handle, c)
    uv_error("write",err)
    1
end

function write{T}(f::File, a::Array{T})
    if isbits(T)
        write(f,pointer(a),length(a)*sizeof(eltype(a)))
    else
        invoke(write, (IO, Array), f, a)
    end
end

function truncate(f::File, n::Integer)
    req = Base.c_malloc(_sizeof_uv_fs)
    err = ccall(:uv_fs_ftruncate,Int32,(Ptr{Void},Ptr{Void},Int32,Int64,Ptr{Void}),
                eventloop(),req,f.handle,n,C_NULL)
    c_free(req)
    uv_error("ftruncate", err)
    f
end

function read(f::File, ::Type{Uint8})
    if !f.open
        error("file is not open")
    end
    ret = ccall(:jl_fs_read_byte, Int32, (Int32,), f.handle)
    uv_error("read", ret)
    return uint8(ret)
end

function read!{T}(f::File, a::Array{T}, nel=length(a))
    if nel < 0 || nel > length(a)
        throw(BoundsError())
    end
    if isbits(T)
        nb = nel*sizeof(T)
        ret = ccall(:jl_fs_read, Int32, (Int32, Ptr{Void}, Csize_t),
                    f.handle, a, nb)
        uv_error("read",ret)
    else
        invoke(read, (IO, Array), s, a)
    end
    a
end

nb_available(f::File) = filesize(f) - position(f)

function readbytes!(f::File, b::Array{Uint8}, nb=length(b))
    nr = min(nb, nb_available(f))
    if length(b) < nr
        resize!(b, nr)
    end
    read!(f, b, nr)
    return nr
end
readbytes(io::File) = read!(io, Array(Uint8, nb_available(io)))
readbytes(io::File, nb) = read!(io, Array(Uint8, min(nb, nb_available(io))))

function readbytes(f::File)
    a = Array(Uint8, nb_available(f))
    read!(f,a)
    a
end

const SEEK_SET = int32(0)
const SEEK_CUR = int32(1)
const SEEK_END = int32(2)

function position(f::File)
    ret = ccall(:jl_lseek, Coff_t,(Int32,Coff_t,Int32),f.handle,0,SEEK_CUR)
    systemerror("lseek", ret == -1)
    ret
end

fd(f::File) = RawFD(f.handle)
stat(f::File) = stat(fd(f))

end

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
module QSWalk

using Expokit
using SparseArrays
using LinearAlgebra

include("utils.jl")
include("dirac.jl")
include("operator.jl")
include("demoralization.jl")
include("evolution.jl")

"""
QSWalk provides package for Julia programming language which enables
high-performance analysis of quantum stochastic walks. There are two main
advantages of the presented packages over the existing software. First, it can
be use to describe quantum stochastic walks in the local, as well as global
regime. Second, it enables the user to seamlessly utilize parallel computing
capabilities.

"""
QSWalk

end


text:
#__precompile__()

module FirstOrderSolvers

using ProximalOperators
using SparseArrays
using Printf
import ValueHistories
using LinearAlgebra
import LinearAlgebra: mul!, Transpose, dot, norm

export Feasibility

include("cones.jl")
include("types.jl")
include("status.jl")
include("utilities/conjugategradients.jl")
include("utilities/affinepluslinear.jl")
include("problemforms/HSDE/HSDEStatus.jl")
include("problemforms/HSDE/HSDE.jl")
include("problemforms/HSDE/HSDEAffine.jl")

include("problemforms/Feasibility/Feasibility.jl")
include("problemforms/Feasibility/FeasibilityStatus.jl")

include("FOSSolverInterface.jl")  # MathProgBase interface
include("solverwrapper.jl")

include("wrappers/longstep.jl")
include("wrappers/linesearch.jl")
include("solvers/defaults.jl")
include("solvers/solvers.jl")
#include("lapack.jl")
#include("solvermethods.jl")
#include("line_search_methods.jl")

end # module


text:
## macro for argument checking

macro strip_linenos(expr)
    return esc(Base.remove_linenums!(expr))
end

"""
    @check_args(
        D,
        @setup(statements...),
        (arg₁, cond₁, message₁),
        (cond₂, message₂),
        ...,
    )

A convenience macro that generates checks of arguments for a distribution of type `D`.

The macro expects that a boolean variable of name `check_args` is defined and generates
the following Julia code:
```julia
Distributions.check_args(check_args) do
    \$(statements...)
    cond₁ || throw(DomainError(arg₁, \$(string(D, ": ", message₁))))
    cond₂ || throw(ArgumentError(\$(string(D, ": ", message₂))))
    ...
end
```

The `@setup` argument can be elided if no setup code is needed. Moreover, error messages
can be omitted. In this case the message `"the condition \$(cond) is not satisfied."` is
used.
"""
macro check_args(D, setup_or_check, checks...)
    # Extract setup statements
    if Meta.isexpr(setup_or_check, :macrocall) && setup_or_check.args[1] == Symbol("@setup")
        setup_stmts = Any[esc(ex) for ex in setup_or_check.args[3:end]]
    else
        setup_stmts = []
        checks = (setup_or_check, checks...)
    end

    # Generate expressions for each condition
    conds_exprs = map(checks) do check
        if Meta.isexpr(check, :tuple, 3)
            # argument, condition, and message specified
            arg = check.args[1]
            cond = check.args[2]
            message = string(D, ": ", check.args[3])
            return :(($(esc(cond))) || throw(DomainError($(esc(arg)), $message)))
        elseif Meta.isexpr(check, :tuple, 2)
            cond_or_message = check.args[2]
            if cond_or_message isa String
                # only condition and message specified
                cond = check.args[1]
                message = string(D, ": ", cond_or_message)
                return :(($(esc(cond))) || throw(ArgumentError($message)))
            else
                # only argument and condition specified
                arg = check.args[1]
                cond = cond_or_message
                message = string(D, ": the condition ", cond, " is not satisfied.")
                return :(($(esc(cond))) || throw(DomainError($(esc(arg)), $message)))
            end
        else
            # only condition specified
            cond = check
            message = string(D, ": the condition ", cond, " is not satisfied.")
            return :(($(esc(cond))) || throw(ArgumentError($message)))
        end
    end

    return @strip_linenos quote
        Distributions.check_args($(esc(:check_args))) do
            $(__source__)
            $(setup_stmts...)
            $(conds_exprs...)
        end
    end
end

"""
    check_args(f, check::Bool)

Perform check of arguments by calling a function `f`.

If `check` is `false`, the checks are skipped.
"""
function check_args(f::F, check::Bool) where {F}
    check && f()
    nothing
end

ChainRulesCore.@non_differentiable check_args(::Any, ::Bool)

##### Utility functions

isunitvec(v::AbstractVector) = (norm(v) - 1.0) < 1.0e-12

isprobvec(p::AbstractVector{<:Real}) =
    all(x -> x ≥ zero(x), p) && isapprox(sum(p), one(eltype(p)))

# get a type wide enough to represent all a distributions's parameters
# (if the distribution is parametric)
# if the distribution is not parametric, we need this to be a float so that
# inplace pdf calculations, etc. allocate storage correctly
@inline partype(::Distribution) = Float64

# because X == X' keeps failing due to floating point nonsense
function isApproxSymmmetric(a::AbstractMatrix{Float64})
    tmp = true
    for j in 2:size(a, 1)
        for i in 1:(j - 1)
            tmp &= abs(a[i, j] - a[j, i]) < 1e-8
        end
    end
    return tmp
end

"""
    ispossemdef(A, k) -> Bool
Test whether a matrix is positive semi-definite with specified rank `k` by
checking that `k` of its eigenvalues are positive and the rest are zero.
# Examples
```jldoctest; setup = :(using Distributions: ispossemdef)
julia> A = [1 0; 0 0]
2×2 Matrix{Int64}:
 1  0
 0  0
julia> ispossemdef(A, 1)
true
julia> ispossemdef(A, 2)
false
```
"""
function ispossemdef(X::AbstractMatrix, k::Int;
                     atol::Real=0.0,
                     rtol::Real=(minimum(size(X))*eps(real(float(one(eltype(X))))))*iszero(atol))
    _check_rank_range(k, minimum(size(X)))
    ishermitian(X) || return false
    dp, dz, dn = eigsigns(Hermitian(X), atol, rtol)
    return dn == 0 && dp == k
end
function ispossemdef(X::AbstractMatrix;
                     atol::Real=0.0,
                     rtol::Real=(minimum(size(X))*eps(real(float(one(eltype(X))))))*iszero(atol))
    ishermitian(X) || return false
    dp, dz, dn = eigsigns(Hermitian(X), atol, rtol)
    return dn == 0
end

function _check_rank_range(k::Int, n::Int)
    0 <= k <= n || throw(ArgumentError("rank must be between 0 and $(n) (inclusive)"))
    nothing
end

#  return counts of the number of positive, zero, and negative eigenvalues
function eigsigns(X::AbstractMatrix,
                  atol::Real=0.0,
                  rtol::Real=(minimum(size(X))*eps(real(float(one(eltype(X))))))*iszero(atol))
    eigs = eigvals(X)
    eigsigns(eigs, atol, rtol)
end
function eigsigns(eigs::Vector{<: Real}, atol::Real, rtol::Real)
    tol = max(atol, rtol * eigs[end])
    eigsigns(eigs, tol)
end
function eigsigns(eigs::Vector{<: Real}, tol::Real)
    dp = count(x -> tol < x, eigs)        #  number of positive eigenvalues
    dz = count(x -> -tol < x < tol, eigs) #  number of numerically zero eigenvalues
    dn = count(x -> x < -tol, eigs)       #  number of negative eigenvalues
    return dp, dz, dn
end


text:
#######################################################################
#
# An example of creating of a Pareto chart with  XlsxWriter.
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("chart_pareto.xlsx")
	ws = add_worksheet!(wb)
	# Formats used in the workbook.
	bold = add_format!(wb, Dict("bold"=> 1))

	percent_format = add_format!(wb, Dict("num_format"=> "0.0%"))

	# Widen the columns for visibility.
	set_column!(ws, "A:A", 15)
	set_column!(ws, "B:C", 10)

	# Add the worksheet data that the charts will refer to.
	headings = ["Reason", "Number", "Percentage"]

	reasons = [
		"Traffic", "Child care", "Public Transport", "Weather",
		"Overslept", "Emergency",
	]

	numbers  = [60,   40,    20,  15,  10,    5]
	percents = [0.44, 0.667, 0.8, 0.9, 0.967, 1]

	write_row!(ws, "A1", headings, fmt=bold)
	write_column!(ws, "A2", reasons)
	write_column!(ws, "B2", numbers)
	write_column!(ws, "C2", percents, fmt=percent_format)

	# Create a new column chart. This will be the primary chart.
	column_chart = add_chart!(wb, Dict("type"=> "column"))

	# Add a series.
	add_series!(column_chart, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$B\$2:\$B\$7",
	))

	# Add a chart title.
	set_title!(column_chart, Dict("name"=> "Reasons for lateness"))

	# Turn off the chart legend.
	set_legend!(column_chart, Dict("position"=> "none"))

	# Set the title and scale of the Y axes. Note, the secondary axis is set from
	# the primary chart.
	set_y_axis!(column_chart, Dict(
		"name"=> "Respondents (number)",
		"min"=> 0,
		"max"=> 120
	))
	set_y2_axis!(column_chart, Dict("max"=> 1))

	# Create a new line chart. This will be the secondary chart.
	line_chart = add_chart!(wb, Dict("type"=> "line"))

	# Add a series, on the secondary axis.
	add_series!(line_chart, Dict(
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$C\$2:\$C\$7",
		"marker"=>     Dict("type"=> "automatic"),
		"y2_axis"=>    1,
	))

	# Combine the charts.
	combine!(column_chart, line_chart)

	# Insert the chart into the worksheet.
	insert_chart!(ws, "F2", column_chart)

	close(wb)
	isfile("chart_pareto.xlsx")
end

test()


text:
import Reinforce: action, actions, finished, ismdp, maxsteps, reset!, reward, state, step!

###############################################################################
#  Example Usage
###############################################################################

mutable struct FooEnv <: AbstractEnvironment
  s::Int  # state
  r::Int  # reward

  FooEnv() = new(1, 0)
end

state(env::FooEnv) = env.s
reward(env::FooEnv) = env.r
reset!(env::FooEnv) = (env.s = 1; env.r = 0; env)
step!(env::FooEnv, s, a) = (env.s += 1; env.r = -1; (env.r, env.s))
maxsteps(env::FooEnv) = 3
actions(env::FooEnv, s′) = [1, 2, 3]

struct FooPolicy <: AbstractPolicy
end

action(π::FooPolicy, r, s, A) = rand(A)

# Iterating a Episode:
#
# ep = Episode(FooEnv(), FooPolicy())
# for (s, a, r, s′) ∈ ep
#   ...
# end


text:
using Test
using QContinuousMeasurement
using HDF5
using Random

typedict(x) = Dict(string(fn) => getfield(x, fn) for fn ∈ fieldnames(typeof(x)))

# Opens the FileWriter
@testset "FileWriter" begin
    Ntraj = 100

    params = ModelParameters(outpoints=200)

    timevector = collect(get_time(params))

    filename = tempname() # Get a temporary filename

    @info "Temporary file: $filename"
    FI = rand(length(timevector), Ntraj)
    QFI = rand(length(timevector), Ntraj)
    xi2y = rand(length(timevector), Ntraj)
    FIstrong = rand(length(timevector), Ntraj)

    writer = FileWriter(filename, params, Ntraj, ["FI", "QFI", "xi2y", "something", "FIstrong"])

    for i = 1:Ntraj
        tuple = (FI = FI[:, i],
                 QFI = QFI[:, i],
                 xi2y = xi2y[:, i],
                 FIstrong = FIstrong[:, i],
                 extrafield = xi2y[:, i]) # This field should be ignored
        put!(writer, tuple)
    end

    close(writer)

    fid = h5open(filename, "r")

    @testset "datasets" begin
        @test read(fid["t"]) == timevector
        @test read(fid["FI"]) == FI
        @test read(fid["QFI"]) == QFI
        @test read(fid["xi2y"]) == xi2y
        @test iszero(read(fid["something"]))
        @test read(attributes(fid)["stored_traj"]) == Ntraj
        close(fid)
    end
    @testset "attrs" begin
        h5attr = h5readattr(filename, "/")
        pop!(h5attr, "stored_traj")
        @test h5attr == typedict(params)
    end

end


text:
module MyExample

using ForwardDiff

include("extral_file.jl")
export my_f, derivative_of_my_f
end # module


text:
# # clip-06-10.jl

# Load Julia packages (libraries) needed  for the snippets in chapter 0

using StatisticalRethinking, Optim
#gr(size=(600,600));

# ### snippet 3.2

p_grid = range(0, step=0.001, stop=1)
prior = ones(length(p_grid))
likelihood = [pdf(Binomial(9, p), 6) for p in p_grid]
posterior = likelihood .* prior
posterior = posterior / sum(posterior)

# ### snippet 3.3
# Draw 10000 samples from this posterior distribution

N = 10000
samples = sample(p_grid, Weights(posterior), N);

# In StatisticalRethinkingJulia samples will always be stored
# in an MCMCChains.Chains object.

chn = MCMCChains.Chains(reshape(samples, N, 1, 1), ["toss"]);

# Describe the chain

MCMCChains.describe(chn)

# ### snippet 3.6

v = 0.0
for i in 1:length(p_grid)
  global v
  if p_grid[i] < 0.5
    v += posterior[i]
  end
end
v

# ### snippet 3.7

mapreduce(p -> p < 0.5 ? 1 : 0, +, samples) / N   |> display

# ### snippet 3.8

mapreduce(p -> (p > 0.5 && p < 0.75) ? 1 : 0, +, samples) / N   |> display

# ### snippet 3.9

quantile(samples, 0.8)

# ### snippet 3.10

quantile(samples, [0.1, 0.9])

# End of `03/clip-06-10.jl`


text:
using RecipesBase

function xs_ys_barline(walls, ys)
    @argcheck length(walls) == length(ys) + 1
    xs_ret = Float64[]
    ys_ret = Float64[]

    push!(xs_ret, walls[1])
    push!(ys_ret, 0)

    for i in 1:length(walls)-1
        x_left = walls[i]
        x_right = walls[i+1]
        y = ys[i]
        push!(xs_ret, x_left)
        push!(ys_ret, y)
        push!(xs_ret, x_right)
        push!(ys_ret, y)
    end
    push!(xs_ret, walls[end])
    push!(ys_ret, 0)

    xs_ret, ys_ret
end

function spacings(xs)
    xs[2:end] - xs[1:end-1]
end

@recipe function plot(grid::VegasGrid{1})
    # seriestype --> :vline
    # walls = first(grid.boundaries)
    xs = first(grid.boundaries)
    if get(plotattributes, :seriestype, nothing) == :vline
        xs
    else
        ys = 1 ./ spacings(xs)
        xs_ys_barline(xs, ys)
    end
end


text:
import Statistics: mean
using Printf
using DelimitedFiles
import Images
import ImageFiltering
using Base.Threads
using MultivariateStats
using Statistics
#data augmentation
"""
    flip_volume_front(x, y)

Flip a volume so that the longitudinal (x) axis becomes the vertical (z) axis
and adjust the landmakr coordinates accordingly.
"""
function flip_volume_front(x, y)
  newx = zeros(Float32, size(x,3),size(x,2),size(x,1),size(x,4))
  newy = deepcopy(y)
  Threads.@threads for ind in 1:size(y,2)
    for z in 1:size(x, 1)
      newx[:,:,z,ind] = transpose(x[z,:,:,ind])
    end
    for cor in 1:3:size(y,1)
      newy[cor, ind] = y[cor+2, ind]
      newy[cor+1, ind] = y[cor+1, ind]
      newy[cor+2, ind] = y[cor, ind]
    end
  end
  return newx, newy
end

"""
    flip_volume_side(x, y)

Flip a volume so that the lateral (y) axis becomes the vertical (z) axis
and adjust the landmark coordinates accordingly.
"""
function flip_volume_side(x, y)
  newx = zeros(Float32, size(x,2),size(x,3),size(x,2),size(x,4))
  newy = deepcopy(y)
  Threads.@threads for ind in 1:size(y,2)
    for z in 1:size(x, 2)
      newx[:,:,z,ind] = x[:,z,:,ind]
    end
    for cor in 1:3:size(y,1)
      newy[cor, ind] = y[cor, ind]
      newy[cor+1, ind] = y[cor+2, ind]
      newy[cor+2, ind] = y[cor+1, ind]
    end
  end
  return newx, newy
end

"""
    function mirror_vol(x, y)

A form of data augmentation. Mirrors a volume (exchanges x and y axis) and returns
concatenation of original data and mirrored data.
"""
function mirror_vol(x, y)
  copyx = x
  copyy = y
  copyx = cat(x, copyx, dims=length(size(x)))
  copyy = hcat(y, copyy)
  inds = size(x)[length(size(x))]
  for ind in 1:inds
    for slice in 1:size(x)[3]
      if length(size(x)) == 5
        trans = x[:,:,slice, 1, ind]
        copyx[:,:,slice, 1, ind+inds] = trans'
      else
        trans = x[:,:,slice, ind]
        copyx[:,:,slice, ind+inds] = trans'
      end
    end
    for cor in 1:3:size(y)[1]
      copyy[cor, ind+inds] = y[cor+1, ind]
      copyy[cor+1, ind+inds] = y[cor, ind]
    end
  end
  return copyx, copyy
end

"""
  flip_2D(x, y)

Takes 2D images in a 4D tensor and landmark data in 2D tensor and returns
the flipped (clockwise) images, aswell
as the coordinates for the flipped images. only works on square images.
"""
function flip_2D(x, y)
  copyx = zeros(Float32, size(x, 2),size(x, 1),size(x, 3),size(x, 4))
  copyy = deepcopy(y)
  inds = size(y, 2)
  size_x = size(x, 1)
  for i in 1:inds
    for cor in 1:2:size(y, 1)
      copyy[cor+1, i] = size(x,2)/10 - y[cor, i]
      copyy[cor, i] = y[cor+1, i]
    end
    for img in 1:size(x, 3)
      for o in 1:size_x
        copyx[:, size_x-(o-1), img, i] = x[o, :, img, i]
      end
    end
  end
  return copyx, copyy
end

"""
    flip_3D(x, y)

Takes 3D volumes in a 4D tensor and landmark data in 2D tensor and returns
the flipped (clockwise) volumes, aswell
as the coordinates for the flipped volumes.
"""
function flip_3D(x, y)
  copyx = zeros(Float32, size(x, 2),size(x, 1),size(x, 3),size(x, 4))
  copyy = deepcopy(y)
  inds = size(y, 2)
  size_x = size(x, 1)
  Threads.@threads for i in 1:inds
    for cor in 1:3:size(y, 1)
      copyy[cor+1, i] = size(x,2)/10 - y[cor, i]
      copyy[cor, i] = y[cor+1, i]
    end
    for o in 1:size_x
      copyx[:, size_x-(o-1), :, i] = x[o, :, :, i]
    end
  end
  return copyx, copyy
end

# """
#     jitter(x, y)
#
# Jitters around images based on their corresponding 2D landmarks, so that
# all landmarks and at least 8 pixels around them will still be inside
# the image.
# """
# function jitter(x, y)
#   copyx = zeros(Float32, size(x, 2),size(x, 1),size(x, 3),size(x, 4))
#   copyy = y
#   copyx = cat(x, copyx, dims=length(size(x)))
#   copyy = hcat(y, copyy)
#   inds = size(y, 2)
#   for i in 1:inds
#     max_x = 0
#     max_y = 0
#     for cor in 1:2:size(y,1)
#       if y[cor+1, i] > max_y
#         max_y = round(y[cor+1, i], digits=1)
#       end
#       if y[cor, i] > max_x
#         max_x = round(y[cor, i], digits=1)
#       end
#     end
#     max_x = max_x*10 + 8
#     max_y = max_y*10 + 8
#     minval = minimum(x[:,:,1,i])
#     if max_x >= size(x,1)
#       jitterx = 1
#     else
#       jitterx = rand(1:size(x,1)-floor(Int, max_x))
#       for o in 1:jitterx
#         copyx[o, :, 1, inds + i] .= minval
#       end
#     end
#     if max_y >= size(x,2)
#       jittery = 1
#     else
#       jittery = rand(1:size(x,2)-floor(Int, max_y))
#       for o in 1: jittery
#         copyx[:, o, 1, inds + i] .= minval
#       end
#     end
#     for o in jitterx:size(x,1)
#       copyx[o, jittery+1:end,1,inds+i] = x[(o-(jitterx-1)), 1:end-jittery,1,i]
#     end
#     for o in 1:2:size(y,1)
#       copyy[o,inds+i] = y[o,i] + jitterx/10
#       copyy[o+1,inds+i] = y[o+1,i] + jittery/10
#     end
#   end
#   return copyx, copyy
# end

"""
    jitter_2D(imgs, landmarks, padding)

Jitters around 2D images based on their landmarks. The landmarked part of the
images plus #(padding) pixels around them will not be moved out of the
picture. returns jittered images and landmarks.
"""
function jitter_2D(imgs, landmarks, padding)
  coordinates = to_3d_array_2d(landmarks)
  out = deepcopy(imgs)
  Threads.@threads for ind in 1:size(landmarks, 2)
    size_x = size(imgs, 1)
    size_y = size(imgs, 2)
    min_x = max(0, minimum(coordinates[:,1,ind]) * 10 - padding)
    max_x = min(size_x, maximum(coordinates[:,1,ind]) * 10 + padding)
    min_y = max(0, minimum(coordinates[:,2,ind]) * 10 - padding)
    max_y = min(size_y, maximum(coordinates[:,2,ind]) * 10 + padding)
    jitter_x = rand(-floor(Int, min_x):(size_x-floor(Int, max_x)))
    jitter_y = rand(-floor(Int, min_y):(size_y-floor(Int, max_y)))
    for img in 1:size(imgs, 3)
        filler = maximum(imgs[:,:,img,ind])
        out[:,:,img,ind] .= filler
        out[max(1, jitter_x+1):min(size_x, size_x+jitter_x), max(1, jitter_y+1):min(size_y, size_y+jitter_y),
          img, ind] .= imgs[max(1, -jitter_x+1):min(size_x, size_x-jitter_x),
          max(1, -jitter_y+1):min(size_y, size_y-jitter_y), img, ind]
    end
    coordinates[:,1,ind] .= coordinates[:,1,ind] .+ jitter_x/10
    coordinates[:,2,ind] .= coordinates[:,2,ind] .+ jitter_y/10
  end
  lms_out = to_2d_array_2d(coordinates)
  return out, lms_out
end

"""
    jitter_3D(volumes, landmarks, padding)

Jitters around binarized volumes based on their corresponding x/y landmarks, so that
the relevant object (smaller voxel value) will still be fully inside the volume but
moved around randomly inside the volume.
"""
function jitter_3D(volumes, landmarks, padding)
  coordinates = to_3d_array(landmarks)
  out = deepcopy(volumes)
  Threads.@threads for ind in 1:size(landmarks, 2)
    filler = maximum(volumes[:,:,:,ind])
    size_x = size(volumes[:,:,:,ind], 1)
    size_y = size(volumes[:,:,:,ind], 2)
    size_z = size(volumes[:,:,:,ind], 3)
    min_x = max(0, minimum(coordinates[:,1,ind]) * 10 - padding)
    max_x = min(size_x, maximum(coordinates[:,1,ind]) * 10 + padding)
    min_y = max(0, minimum(coordinates[:,2,ind]) * 10 - padding)
    max_y = min(size_y, maximum(coordinates[:,2,ind]) * 10 + padding)
    min_z = max(0, minimum(coordinates[:,3,ind]) * 10 - padding)
    max_z = min(size_z, maximum(coordinates[:,3,ind]) * 10 + padding)
    jitter_x = rand(-floor(Int, min_x):(size_x-floor(Int, max_x)))
    jitter_y = rand(-floor(Int, min_y):(size_y-floor(Int, max_y)))
    jitter_z = rand(-floor(Int, min_z):(size_z-floor(Int, max_z)))
    out[:,:,:,ind] .= filler
    out[max(1, jitter_x+1):min(size_x, size_x+jitter_x), max(1, jitter_y+1):min(size_y, size_y+jitter_y),
      max(1, jitter_z+1):min(size_z, size_z+jitter_z), ind] .= volumes[max(1, -jitter_x+1):min(size_x, size_x-jitter_x),
      max(1, -jitter_y+1):min(size_y, size_y-jitter_y),
        max(1, -jitter_z+1):min(size_z, size_z-jitter_z), ind]
    coordinates[:,1,ind] .= coordinates[:,1,ind] .+ jitter_x/10
    coordinates[:,2,ind] .= coordinates[:,2,ind] .+ jitter_y/10
    coordinates[:,3,ind] .= coordinates[:,3,ind] .+ jitter_z/10
  end
  lms_out = to_2d_array(coordinates)
  return out, lms_out
end

"""
    rotate_2d(img, lms, deg)

Rotates an Image by deg degrees around its center in a counter clockwise direction.
The 2d landmarks (coordinates x individuals) will be adjusted to the rotated image.
Returns rotated image and adjusted landmarks.
"""
function rotate_2d(img, lms, deg)
  deg = deg*pi/180
  rotated = Images.imrotate(img, -deg)
  size_x = size(rotated, 1)
  size_y = size(rotated, 2)
  out = Images.imresize(rotated, (size(img, 1), size(img, 2)))
  n_points = floor(Int, length(lms)/2)
  coords = zeros(2,n_points)
  coords_out = zeros(n_points*2, 1)
  for i in 1:n_points
    coords[1,i] = lms[i*2-1]
    coords[2,i] = lms[i*2]
  end
  coords[1,:] = coords[1,:] .- (size(img,1)/2)/10
  coords[2,:] = coords[2,:] .- (size(img,2)/2)/10
  for i in 1:n_points
    coords_out[i*2-1] = (coords[1,i]*cos(deg) - coords[2,i]*sin(deg) + (size_x/2)/10) * size(img,1)/size_x
    coords_out[i*2] = (coords[1,i]*sin(deg) + coords[2,i]*cos(deg) + (size_y/2)/10) * size(img,2)/size_y
  end
  return out, coords_out
end

"""
    rotate_images(img, lms, deg)

Takes a 3d tensor [res1 x res1 x n] with n images and rotates all of them
in a counterclockwise direction around their center by deg degrees. Will adjust the 2d coordinates
of the landmark array [coords x n]. Returns rotated images and adjusted landmarks.
"""
function rotate_images(imgs, lms, deg)
  n_inds = size(imgs, 4)
  out = deepcopy(imgs)
  lms_out = deepcopy(lms)
  for i in 1:n_inds
    for img in 1:size(imgs, 3)
      filler = maximum(imgs[:,:,img,i])
      if img == 1
        rotated, rotated_lms = rotate_2d(imgs[:,:,img,i], lms[:,i], deg)
        lms_out[:,i:i] .= rotated_lms
        rotated[findall(x->isnan(x)==true, rotated)] .= filler
      else
        rotated_img = Images.imrotate(imgs[:,:,img,i], -deg*pi/180)
        rotated = Images.imresize(rotated_img, (size(imgs, 1), size(imgs, 2)))
        rotated[findall(x->isnan(x)==true, rotated)] .= filler
      end
      out[:,:,img,i] .= rotated

    end
  end
  return out, lms_out
end

"""
    rotate_volumes(vols, lms, deg)

Takes a 4d tensor [res1 x res2 x res3 x n] with n volumes and rotates all of them
around their center by deg degrees. The rotation is along the z-axis, so the z-coordinates
will not be affected. x- and y coordinates of the landmark array [coords x n] will be adjusted.
Returns rotated volumes and adjusted landmarks.
"""
function rotate_volumes(vols, lms, deg)
  n_inds = size(lms, 2)
  n_points = floor(Int, size(lms, 1)/3)
  x_y_only = zeros(n_points*2, n_inds)
  for i in 1:n_points
    x_y_only[i*2-1:i*2,:] .= lms[i*3-2:i*3-1,:]
  end
  x_y_rotated = deepcopy(x_y_only)
  out = deepcopy(vols)
  out_lms = deepcopy(lms)
  for i in 1:n_inds
    filler = maximum(vols[:,:,:,i])
    for img in 1:size(vols, 3)
      if img == 1
        rotated, rotated_lms = rotate_2d(vols[:,:,img,i], x_y_only[:,i], deg)
        rotated[findall(x->isnan(x)==true, rotated)] .= filler
        out[:,:,img,i] .= rotated
        x_y_rotated[:,i] .= rotated_lms[:,1]
      else
        rotated_img = Images.imrotate(vols[:,:,img,i], -deg*pi/180)
        rotated = Images.imresize(rotated_img, (size(vols, 1), size(vols, 2)))
        rotated[findall(x->isnan(x)==true, rotated)] .= filler
        out[:,:,img,i] .= rotated
      end
    end
  end
  for i in 1:n_points
    out_lms[i*3-2:i*3-1,:] .= x_y_rotated[i*2-1:i*2,:]
  end
  return out, out_lms
end

max_stars_repo_name: imi-diener/AutoLM

max_issues_repo_name: imi-diener/AutoLM

max_forks_repo_name: imi-diener/AutoLM


text:
function part1(input)
    x = sort(parse.(Int, split(readchomp(input), ",")))
    return sum(abs.(x .- x[(end + 1) ÷ 2]))
end

part2_cost(n) = abs(n) * (abs(n) + 1) / 2
part2_total_cost(x, target) = sum(part2_cost(xi - target) for xi in x)

function part2(input)
    x = parse.(Int, split(readchomp(input), ","))
    target = sum(x) / length(x)
    return Int(minimum(part2_total_cost.(Ref(x), floor(target):ceil(target))))
end


text:
################################################################################
#
#          EllCrv/EllCrvFin.jl : Elliptic curves over finite fields
#
################################################################################

################################################################################
#
#  Submodules
#
################################################################################

include("EllCrv.jl")
include("../src/Misc/BigIntInterface.jl")
include("../src/Misc/PolyFac.jl")

################################################################################
#
#  Imports
#
################################################################################

using Random

################################################################################
#
#  Exports
#
################################################################################

export order, order_via_schoof, rand_point, crypto_curve

################################################################################
#
#  Random Point Generator
#
################################################################################

function issquare(x::Union{AbstractAlgebra.gfelem{Int32}, AbstractAlgebra.FinFieldElem})
    k = parent(x)
    S, t = AbstractAlgebra.PolynomialRing(k, "t", cached = false)

    f = t^2 - x

    # Distinct Degree Factorization for now.
    fac = factor_ddf(f)

    p = first(keys(fac.fac))

    if fac[p] == 2
        root = -AbstractAlgebra.coeff(p, 0)
        return true, k(root)
    elseif length(fac) == 2
        root = -coeff(p, 0)
        return true, k(root)
    else
        return false, zero(k)
    end
end

function rand_point(E::EllCrv)
    k = E.base_field
    (a4, a6) = a_invars(E)

    if E.short == false
        error("Not implemented for long form.")
    end

    while true

        x = rand(k)
        square = x^3 + a4*x + a6

        a = AbstractAlgebra.issquare(square)
        if a[1] == true
            y = a[2]
            P = E((x, y))
            return P
        end
    end
end

################################################################################
#
#  Elliptic Curve Generator over Finite Fields (Naive implementation at first.)
#  https://pdfs.semanticscholar.org/c60c/46a9ee6f90958f609f464723f6bf92835feb.pdf
#
################################################################################

function crypto_curve(r0::BigInt, k0::BigInt, h0::BigInt)
    if prod(r0, k0, h0) <= 0
        error("r0, k0, h0 must be positive integers.")
    end

    if r0 < 2^BigInt(159) || k0 < 4 || h0 < 200
        error("Invalid input.")
    end

    (d, p, r, k) = find_prime(r0, k0, h0)

    (E, G) = find_curve(d, p, r, k)

    return (d, p, r, k, E, G)
end

function find_prime(r0::BigInt, k0::BigInt, h0::BigInt)
    print("Enter a (rational) prime p: ")
    input = readline()
    p = parse(BigInt, chomp(input))

    m = 50

    if p == 0
        find_prime_delta_fixed(r0, k0, h0)
    elseif is_prime(p, m) == true && floor(log2(p)) == floor(log2(r0*k0))
        find_discriminant(p, r0, k0, h0)
    else
        error("Invalid prime.")
    end
end

function find_prime_delta_fixed(r0::BigInt, k0::BigInt, h0::BigInt)
    print("Enter a discriminant: ")
    input = readline()
    d = parse(BigInt, chomp(input))

    abs_d = abs(d)

    if k0 >= 4
        if mod(d, 8) == 1 || mod(abs_d, 8) == 7
            b = BigInt(floor(log2(4*r0)))
            i0 = 2^b
            i1 = 2^(b + 1) - 1
            interval = collect(BigInt, i0:i1)

            if find_prime_1_mod_8(r0, 4, d, interval)
                return true
            else
                error("Unable to find prime.")
            end
        else
            error("Get discriminant congruent to 1 mod 8,
            or its norm congruent to 7 mod 8.")
        end
    elseif k0 >= 2
        if mod(d, 16) == 8 || 12
            b = BigInt(floor(log2(r0*k0)))
            i0 = 2^b
            i1 = 2^(b + 1) - 1
            interval = collect(BigInt, i0:i1)

            if find_prime_0_mod_4(r0, k0, d, interval)
                return true
            else
                error("Unable to find prime.")
            end
        else
            error("Get discriminant congruent to 8 or 12 modulo 16.")
        end
    else
        if mod(d, 8) == 5
            b = BigInt(floor(log2(r0)))
            i0 = 2^b
            i1 = 2^(b + 1) - 1
            interval = collect(BigInt, i0:i1)

            if find_prime_5_mod_8(r0, 1, d, interval)
                return true
            else
                error("Unable to find prime.")
            end
        else
            error("Get discrimint congruent to 5 mod 8.")
        end
    end
end

function find_prime_0_mod_4(r0::BigInt, k0::BigInt, d::BigInt,
    interval::Array{BigInt, 1})
    # Number of successive tested pairs (t, y).
    # If not successful after T tries, init a new pair.
    const T = 2000

    # Bit of 2 in binary representation of t.
    t_bit_1 = Int

    abs_delta = assign(d)

    ###
end


text:
# This file was generated, do not modify it. # hide
for m in models(matching(X, y))
    println("Model name = ",m.name,", ","Prediction type = ",m.prediction_type,", ","Package name = ",m.package_name);
end

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using TargetedLearning: Parameters

facts("Testing Parameters") do
    context("Regimen") do
        @fact regimen(1) --> StaticRegimen(1.0)
        d = regimen([1, 0, 1])
        @fact d.a --> [1.0, 0.0, 1.0]
        @fact d === regimen(d) --> true
    end

    srand(412)
    n,p = 100, 10
    w = [ones(n) rand(n,p)]
    a = round(rand(n))
    y = round(rand(n))
    gn1 = predict(lreg(w,a), w)

    context("fluccovar") do
        ey1 = Mean(1.0)
        @fact fluccovar(ey1, a) --> a
        @fact fluccovar(ey1, ones(n)) --> ones(n)

        ate = ATE(ones(n), 0.0)
        @fact fluccovar(ate, a) --> 2a-1
        @fact fluccovar(ate, ones(n)) --> ones(n)
        @fact fluccovar(ate, zeros(n)) --> -ones(n)

        d = [ones(div(n,2)); zeros(div(n,2))]
        eyd = Mean(d)
        gnd = ifelse(d.==1, gn1, 1-gn1)
        hAA = ifelse(a.==d, 1./gnd, 0.0)
        @fact fluccovar(eyd, a)./gnd --> hAA
        hA1 = ifelse(1.==d, 1./gnd, 0.0)
        @fact fluccovar(eyd, ones(n)) --> float(d .== 1)
        hA0 = ifelse(0.==d, 1./gnd, 0.0)
        @fact fluccovar(eyd, zeros(n)) --> float(d .== 0)
    end
end


text:
using PlotlyJS, FileIO, JLD2, SeisIO

# This script makes two time_shifted or stretched traces for validation of SeisXcorrelation

#----------------------------------#
dt     = 0.01; # sampling time of waveform
# control time shift
tarrival_1 = 30 #wave arrival time at station 1
tarrival_2 = 30.05 #wave arrival time at station 2
# control relative frequency
f1 = 0.1 # source sine frequency Hz
f2 = f1*1.01

npts = 6000; # length of time series

foname = "./sinStretchTest.jld2"
#----------------------------------#
#Period of sine
T1 = 1/f1
T2 = 1/f2

npts_sine1 = round(Int64, T1/dt)
npts_sine2 = round(Int64, T2/dt)

# make waveform time series
s1   = sin.(LinRange(0, 2*pi, npts_sine1));
s2   = sin.(LinRange(0, 2*pi, npts_sine2));

u1 = zeros(Float64, npts)
u2 = zeros(Float64, npts)

tp1 = round(Int64, tarrival_1/dt)
tp2 = round(Int64, tarrival_2/dt)

u1[tp1:tp1+npts_sine1-1] = s1
u2[tp2:tp2+npts_sine2-1] = s2

tvec  = collect(( 0:npts - 1 ) .* dt);

Δtlag_true = tarrival_2 - tarrival_1 # true arrival time, which should be equal to the time at maximam coherency of cross-correlation function

# package generated signal in SeisChannel and save
t=[1 0; npts 0]
Ch1 = SeisChannel(name="BP.u1u1.fst", id="BP.u1u1.fst", x=u1, fs=1/dt, t=t)
Ch1.misc["dlerror"] = 0
Ch2 = SeisChannel(name="BP.u2u2.snd", id="BP.u2u2.snd", x=u2, fs=1/dt, t=t)
Ch2.misc["dlerror"] = 0
jldopen(foname, "a+") do file
    file["info/DLtimestamplist"] = ["time1"]
    file["info/stationlist"] = ["BP.u1u1.fst", "BP.u2u2.snd"]
    file["info/DL_time_unit"] = npts*dt
    file["info/f1"] = f1
    file["info/f2"] = f2
    file["info/Δtlag_true"] = Δtlag_true
    file["time1/BP.u1u1.fst"] = Ch1
    file["time1/BP.u2u2.snd"] = Ch2
end

#plot curve
function lineplot1()
    tr1 = scatter(;x=tvec, y=u1, mode="lines+markers", name="u1")
    tr2 = scatter(;x=tvec, y=u2, mode="lines+markers", name="u2")
    layout = Layout(
        xaxis=attr(title="u1"),
        yaxis=attr(title="u2"),
        title="True time shift = $Δtlag_true [s]"
    )
    p = plot([tr1, tr2], layout)
    display(p)
    readline()
end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    LinearPath

Traverse an object with `N` elements in order `1:N`.
"""
struct LinearPath <: Path end

traverse(object, ::LinearPath) = 1:nelements(object)


text:
using GR
using LaTeXStrings

f(x) = x*exp(-2x*im)
x = 0:0.01:π
w = f.(x)

plot(real.(w), imag.(w), x, linewidth=5, title=L"f(x) = xe^{-2xi} \forall x \in [0,\pi]", xlabel=L"\Re e\{f(x)\}", ylabel=L"\Im e\{f(x)\}")


text:
# Perfect matching in complete bipartite graphs.

## Components shared by both correlated and uncorrelated arms.

abstract type PerfectBipartiteMatchingSolver end
function build!(solver::PerfectBipartiteMatchingSolver, n_arms::Int) end

abstract type PerfectBipartiteMatching <: CombinatorialInstance{Tuple{Int, Int}} end

function initial_state(instance::PerfectBipartiteMatching)
  n = round(Int, sqrt(instance.n_arms))
  zero_counts = Dict((i, j) => 0 for i in 1:n, j in 1:n)
  zero_rewards = Dict((i, j) => 0.0 for i in 1:n, j in 1:n)
  return State{Tuple{Int, Int}}(0, 0.0, 0.0, zero_counts, zero_rewards, copy(zero_rewards))
end

function is_feasible(instance::PerfectBipartiteMatching, arms::Vector{Tuple{Int, Int}})
  if length(arms) > instance.n_arms
    return false
  end

  # Check whether each node is taken at most once on each side.
  lefts = sort([arm[1] for arm in arms])
  rights = sort([arm[2] for arm in arms])

  # No duplicates on either side.
  return lefts == unique(lefts) && rights == unique(rights)
end

# Any partial matching is feasible.
is_partially_acceptable(instance::PerfectBipartiteMatching, arms::Vector{Tuple{Int, Int}}) = is_feasible(instance, arms)

## Uncorrelated arms.

struct UncorrelatedPerfectBipartiteMatching <: PerfectBipartiteMatching
  # Mandatory properties.
  n_arms::Int
  optimal_average_reward::Float64

  # Probability distributions for the arm rewards.
  reward::Matrix{Distribution}

  # Internal solver.
  solver::PerfectBipartiteMatchingSolver

  function UncorrelatedPerfectBipartiteMatching(reward::Matrix{Distribution}, solver::PerfectBipartiteMatchingSolver)
    n = size(reward, 1)
    if n != size(reward, 2)
      error("Graph is not bipartite complete: reward matrix must be square.")
    end

    # Prepare the solver to be used (if required).
    build!(solver, n)

    avg_reward = Dict{Tuple{Int, Int}, Float64}(i => mean(reward[i...]) for i in all_arm_indices(reward))
    opt_sol = solve_linear(solver, avg_reward)
    opt = sum(avg_reward[i] for i in opt_sol)

    # Done!
    return new(n ^ 2, opt, reward, solver)
  end
end

## Correlated arms.

struct CorrelatedPerfectBipartiteMatching <: PerfectBipartiteMatching
  # Mandatory properties.
  n_arms::Int
  optimal_average_reward::Float64
  all_arm_indices::Vector{Tuple{Int, Int}}

  # Probability distributions for the arm rewards.
  reward::MultivariateDistribution

  # Internal solver.
  solver::PerfectBipartiteMatchingSolver

  """
  Builds a perfect bipartite matching instance where rewards can be correlated.

  * `reward` is the distribution of the rewards, which can be as complex as required.
    It is supposed to return a vector of size \$\\left(\\# arms\\right)^2\$, with
    element \$\\left(\\# arms \\times (i-1)\\right) + j\$ corresponds to choosing the
    matching between \$i\$ and \$j\$.
  * `solver` is a perfect bipartite matching solver. Nothing is specific to correlated
    matching at this point: the only difference is in the bandit algorithm.
  """
  function CorrelatedPerfectBipartiteMatching(reward::MultivariateDistribution, solver::PerfectBipartiteMatchingSolver)
    n = length(reward)
    n_arms = round(Int, sqrt(n))
    all_arm_indices = vec(Tuple{Int, Int}[(i, j) for i in 1:n_arms, j in 1:n_arms])

    # Prepare the solver to be used (if required).
    build!(solver, n_arms)

    avg_reward_vector = mean(reward)
    avg_reward = Dict{Tuple{Int, Int}, Float64}(i => v for (v, i) in zip(avg_reward_vector, all_arm_indices))
    opt_sol = solve_linear(solver, avg_reward)
    opt = sum(avg_reward[i] for i in opt_sol)

    # Done!
    return new(n, opt, all_arm_indices, reward, solver)
  end
end


text:
using StrainRecon
using Base.Test
# include("../src/asqp.jl")
@testset "ASQP" begin

print("\t\tsolve 2D QP without equality constraints...\n\t\t")
n = 2

G = eye(n)
c = 0*[1;1.]

Ae = zeros(0,n)
be = zeros(0)
me = size(Ae,1)

Ai = eye(n)
bi = 0*rand(n)
mi = size(Ai,1)

x0  = -1*rand(n)

xopt,iter,Wc = asqp(G,c, Ae, be, Ai, bi, x0) # start from infeasible x0 .< 0
@test all((Ai*xopt-bi).>=0)

x0 = rand(n) + 1
xopt,iter,Wc = asqp(G,c, Ae, be, Ai, bi, x0) # start from feasible x0 .> 1
@test all((Ai*xopt-bi).>=0)

# x* = [1;0]
Ai = [1 0.0]
G = eye(2);
c = zeros(2,1);
bi = [1.0]
Ae = zeros(0,2)
be = zeros(0)
xopt, iter, Wc = asqp(G,c,Ae,be,Ai,bi) # start from infeasible x0 = [0;0]
@test all(xopt.==[1;0])
@test all(Ai*xopt-bi.>=0)

x0 = [1;2]
xopt, iter, Wc = asqp(G,c,Ae,be,Ai,bi,x0) # start from feasible x0 = [1;2]
@test all(xopt.==[1;0])
@test all(Ai*xopt-bi.>=0)

# x* = [10, 8]
G = [1 -1; -1 2]
c = [-2; -6]
A = [1 1; -1 2; 2 1]
b = [2; 2; 3]
Ae = zeros(0,2)
be = zeros(0)
xopt, iter, Wc = asqp(G,c,Ae,be,A,b) # start from infeasible x0 = [0;0]
@test all((A*xopt-b).>=0)
@test all(xopt .== [10;8])

x0 = [1;2]
xopt, iter, Wc = asqp(G,c,Ae,be,A,b,x0) # start from feasible x0 = [1;2]
@test all((A*xopt-b).>=0)
@test all(xopt .== [10;8])

# x* = [3;5;7], unconstrained example, G is positive definite, global min exists
G = [5 -2 -1; -2 4 3; -1 3 5]
A = zeros(3,3)
c = [2; -35; -47]
b = [0; 0; 0]
Ae = zeros(0,3)
be = zeros(0)

xopt, iter, Wc = asqp(G,c,Ae,be,A,b)
@test all((A*xopt-b).>=0)
@test norm(xopt - [3.0;5.0;7.0])/norm([3.;5;7]) < 1e-14 # xopt = x* but comparison fails

print("\t\tsolve 2D QP with inequality and equality constraints...\n\t\t")

# x* = [0.5;0.5]

n = 2

G = eye(n)
c = 0*[1;1.]

Ae = [1 1]
be = ones(1)

Ai = eye(n)
bi = 0.0*ones(n)

x0  = [-0.5;0.5]

xopt,iter,Wc = asqp(G,c, Ae, be, Ai, bi, x0) # x0 is infeasible for both equality and inequality constraints
@test all(Ai*xopt.>=bi)
@test norm(Ae*xopt-be) < 1e-14

x0  = -1*rand(n)
xopt,iter,Wc = asqp(G,c, Ae, be, Ai, bi, x0) # x0 is infeasible for both equality and inequality constraints
@test all(Ai*xopt.>=bi)
@test norm(Ae*xopt-be) < 1e-14 # xopt not feasible for equality constraints

x0  = rand(n)+1
xopt,iter,Wc = asqp(G,c, Ae, be, Ai, bi, x0) # x0 is feasible for only inequality constraints
@test all(Ai*xopt.>=bi)
@test norm(Ae*xopt-be) < 1e-14 # xopt not feasible for equality constraints

    # x* = [0.25;0.75]
    G = [4 1;1 2];
    Ai = [1 0;0 1];
    bi = [0;0]
    Ae = [1 1];
    c = [1; 1];
    be = [1];

    x0 = 0*c
    xopt, iter, Wc = asqp(G,c,Ae,be,Ai,bi,x0) # x0 is feasible for only inequality constraint
    @test all(Ai*xopt.>=bi)
    @test norm(Ae*xopt-be) < 1e-14

    x0 = [2;-1]
    xopt, iter, Wc = asqp(G,c,Ae,be,Ai,bi,x0) # x0 is feasible for only equality constraint
    @test all(Ai*xopt.>=bi)
    @test norm(Ae*xopt-be) < 1e-14

    x0 = -1*rand(2)
    xopt, iter, Wc = asqp(G,c,Ae,be,Ai,bi,x0) # x0 is infeasible for inequality and equality constraints
    @test all(Ai*xopt.>=bi)
    @test norm(Ae*xopt-be) < 1e-14

    x0 = [1e-3;1-1e-3]
    xopt, iter, Wc = asqp(G,c,Ae,be,Ai,bi,x0) # x0 is feasible for both inequality and equality constraints
    @test all(Ai*xopt.>=bi)
    @test norm(Ae*xopt-be) < 1e-14

end


text:
# =======================================================
#         CompleteBipartiteGraph struct
# =======================================================

"""
    CompleteBipartiteGraph <: AbstractGraph

A structure representing an undirected bipartite graph.

See also: [`LightGraphs.complete_bipartite_graph`](@ref)
"""
struct CompleteBipartiteGraph{T<:Integer} <: LG.AbstractGraph{T}
    m::T
    n::T

    function CompleteBipartiteGraph{T}(m, n) where {T}

        m = convert(T, m)
        n = convert(T, n)
        m >= zero(T) || throw(ArgumentError("m must be >= 0"))
        n >= zero(T) || throw(ArgumentError("n must be >= 0"))
        _, does_overflow = Base.Checked.add_with_overflow(m, n)
        does_overflow && throw(ArgumentError("m + n must not exceed typemax($T)"))

        return new{T}(m, n)
    end
end

function CompleteBipartiteGraph(m::Integer, n::Integer)
    m, n = promote(m, n)
    return CompleteBipartiteGraph{typeof(m)}(m, n)
end

# =======================================================
#         traits
# =======================================================

LG.is_directed(::Type{<:CompleteBipartiteGraph}) = false

# =======================================================
#        vertices
# =======================================================

Base.eltype(::Type{CompleteBipartiteGraph{T}}) where {T} = T
Base.eltype(g::CompleteBipartiteGraph) = eltype(typeof(g))

LG.nv(g::CompleteBipartiteGraph) = g.m + g.n

LG.vertices(g::CompleteBipartiteGraph) = Base.OneTo(nv(g))

LG.has_vertex(g::CompleteBipartiteGraph, v) = v in vertices(g)

# =======================================================
#        edges
# =======================================================

LG.edgetype(::Type{CompleteBipartiteGraph{T}}) where {T} = Edge{T}
LG.edgetype(g::CompleteBipartiteGraph{T}) where {T} = edgetype(typeof(g))

function LG.ne(g::CompleteBipartiteGraph)
    m = Int(g.m)
    n = Int(g.n)

    return m * n
end

function LG.has_edge(g::CompleteBipartiteGraph{T}, u, v) where {T}

    u, v = minmax(T(u), T(v))

    m = g.m
    n = g.n
    return (u ∈ Base.OneTo(m)) & ((v - m) ∈ (Base.OneTo(n)))
end

# ---- edges vector -----------------------------------

LG.edges(g::CompleteBipartiteGraph) = SimpleEdgeVector(g)

function Base.size(edgevec::SimpleEdgeVector{V, G}) where {V, G <: CompleteBipartiteGraph}
    g = edgevec.graph
    return (ne(g), )
end

@inline function Base.getindex(edgevec::SimpleEdgeVector{V, G}, i::Integer) where {V, G <: CompleteBipartiteGraph}

    @boundscheck i ∈ Base.OneTo(length(edgevec)) || throw(BoundsError(edgevec, i))

    g = edgevec.graph
    T = eltype(g)
    m::T = g.m
    n::T = g.n

    div, rem = T.(divrem(i - 1, n))

    u = div + one(T)
    v = rem + one(T) + m
    return Edge{T}(u, v)
end

Base.IndexStyle(::Type{<:SimpleEdgeVector{V, G}}) where {V, G <: CompleteBipartiteGraph} = IndexLinear()

# =======================================================
#        neighbors
# =======================================================

@inline function LG.outneighbors(g::CompleteBipartiteGraph, v::Integer)
    @boundscheck v ∈ vertices(g) || throw(BoundsError(g, v))
    return OutNeighborVector(g, eltype(g)(v))
end

LG.inneighbors(g::CompleteBipartiteGraph, v::Integer) = outneighbors(g, v)

# ---- neighbors vector -----------------------------------

function Base.size(nbs::OutNeighborVector{V, G}) where {V, G <: CompleteBipartiteGraph}

    g = nbs.graph
    T = eltype(g)
    m::T = g.m
    n::T = g.n
    v::T = nbs.vertex

    return (Int(ifelse(v <= m, n, m)), )
end

@inline function Base.getindex(nbs::OutNeighborVector{V, G}, i::Integer) where {V, G <: CompleteBipartiteGraph}

    @boundscheck i ∈ Base.OneTo(length(nbs)) || throw(BoundsError(nbs, i))

    g = nbs.graph
    T = eltype(g)
    m::T = g.m
    n::T = g.n
    v::T = nbs.vertex

    return ifelse(v <= m, T(i) + m, T(i))
end

Base.IndexStyle(::Type{<:OutNeighborVector{V, G}}) where {V, G <: CompleteBipartiteGraph} = IndexLinear()

# =======================================================
#         converting
# =======================================================

Base.convert(::Type{SimpleGraph}, g::CompleteBipartiteGraph{T}) where {T} = complete_bipartite_graph(g.m, g.n)
Base.convert(::Type{SimpleGraph{T}}, g::CompleteBipartiteGraph) where {T} = complete_bipartite_graph(T(g.m), T(g.n))


text:
using Knet,BenchmarkTools

mul100(a,b)=(for i=1:100; c=a*b; end; Knet.@cuda(cudart,cudaDeviceSynchronize,()))

N = (10,100,1000)
T = (Float32,Float64)
for t in T, m in N, k in N, n in N
    info((t,m,k,n))
    a = KnetArray(rand(t, m, k))
    b = KnetArray(rand(t, k, n))
    r = @benchmarkable mul100($a,$b)
    s = min(1000,div(3e9,(100*m*n*sizeof(t))))
    println(run(r, samples=s))
    gpuinfo("before gc")
    a=b=nothing; gc(); sleep(1)
    gpuinfo("after  gc")
    knetgc()
    gpuinfo("after kgc")
end


text:
##Create callbacks for vaccine function based on 1st Dec 2020 start for sim

function start_vac1_30perc!(integrator)
    integrator.p[29] = fully_vaccinated_30perc[end] / (t_end - t_0)
end
function start_vac1_50perc!(integrator)
    integrator.p[29] = fully_vaccinated_50perc[end] / (t_end - t_0)
end
function start_vac1_70perc!(integrator)
    integrator.p[29] = fully_vaccinated_70perc[end] / (t_end - t_0)
end
function start_vac1_30perc_rapid!(integrator)
    integrator.p[29] = fully_vaccinated_30perc[end] / (t_end_rapid - t_0)
end
function start_vac1_50perc_rapid!(integrator)
    integrator.p[29] = fully_vaccinated_50perc[end] / (t_end_rapid - t_0)
end
function start_vac1_70perc_rapid!(integrator)
    integrator.p[29] = fully_vaccinated_70perc[end] / (t_end_rapid - t_0)
end

function start_vac2_30perc!(integrator)
    integrator.p[30] = fully_vaccinated_30perc[end] / (t_end - t_0)
end
function start_vac2_50perc!(integrator)
    integrator.p[30] = fully_vaccinated_50perc[end] / (t_end - t_0)
end
function start_vac2_70perc!(integrator)
    integrator.p[30] = fully_vaccinated_70perc[end] / (t_end - t_0)
end
function start_vac2_30perc_rapid!(integrator)
    integrator.p[30] = fully_vaccinated_30perc[end] / (t_end_rapid - t_0)
end
function start_vac2_50perc_rapid!(integrator)
    integrator.p[30] = fully_vaccinated_50perc[end] / (t_end_rapid - t_0)
end
function start_vac2_70perc_rapid!(integrator)
    integrator.p[30] = fully_vaccinated_70perc[end] / (t_end_rapid - t_0)
end

function stop_vac1!(integrator)
    integrator.p[29] = 0.0
end
function stop_vac2!(integrator)
    integrator.p[30] = 0.0
end

vac_cbs = [CallbackSet(),
    CallbackSet(
        PresetTimeCallback([(Date(2021, 9, 27) - Date(2020, 12, 1)).value - 56 + 14], start_vac1_30perc!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) - Date(2020, 12, 1)).value + 14], start_vac2_30perc!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) + Month(18) - Date(2020, 12, 1)).value - 56 + 14], stop_vac1!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) + Month(18) - Date(2020, 12, 1)).value + 14], stop_vac2!, save_positions = (false, false))
    ),
    CallbackSet(
        PresetTimeCallback([(Date(2021, 9, 27) - Date(2020, 12, 1)).value - 56 + 14], start_vac1_50perc!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) - Date(2020, 12, 1)).value + 14], start_vac2_50perc!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) + Month(18) - Date(2020, 12, 1)).value - 56 + 14], stop_vac1!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) + Month(18) - Date(2020, 12, 1)).value + 14], stop_vac2!, save_positions = (false, false))
    ),
    CallbackSet(
        PresetTimeCallback([(Date(2021, 9, 27) - Date(2020, 12, 1)).value - 56 + 14], start_vac1_70perc!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) - Date(2020, 12, 1)).value + 14], start_vac2_70perc!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) + Month(18) - Date(2020, 12, 1)).value - 56 + 14], stop_vac1!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) + Month(18) - Date(2020, 12, 1)).value + 14], stop_vac2!, save_positions = (false, false))
    ),
    CallbackSet(
        PresetTimeCallback([(Date(2021, 9, 27) - Date(2020, 12, 1)).value - 56 + 14], start_vac1_30perc_rapid!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) - Date(2020, 12, 1)).value + 14], start_vac2_30perc_rapid!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) + Month(18) - Date(2020, 12, 1)).value - 56 + 14], stop_vac1!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) + Month(18) - Date(2020, 12, 1)).value + 14], stop_vac2!, save_positions = (false, false))
    ),
    CallbackSet(
        PresetTimeCallback([(Date(2021, 9, 27) - Date(2020, 12, 1)).value - 56 + 14], start_vac1_50perc_rapid!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) - Date(2020, 12, 1)).value + 14], start_vac2_50perc_rapid!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) + Month(18) - Date(2020, 12, 1)).value - 56 + 14], stop_vac1!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) + Month(18) - Date(2020, 12, 1)).value + 14], stop_vac2!, save_positions = (false, false))
    ),
    CallbackSet(
        PresetTimeCallback([(Date(2021, 9, 27) - Date(2020, 12, 1)).value - 56 + 14], start_vac1_70perc_rapid!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) - Date(2020, 12, 1)).value + 14], start_vac2_70perc_rapid!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) + Month(18) - Date(2020, 12, 1)).value - 56 + 14], stop_vac1!, save_positions = (false, false)),
        PresetTimeCallback([(Date(2021, 9, 27) + Month(18) - Date(2020, 12, 1)).value + 14], stop_vac2!, save_positions = (false, false))
    )
]