text:
export graph_sastre_exp, graph_sastre_poly, graph_sastre_yks_degopt

"""
    (graph,cref)=graph_sastre_exp(k,method=:auto)

Computes a polynomial evaluation approximating the exponential using `k` matrix
multiplications following a `method` given in the reference. The schemes are
embedded into `degop`-format, see [`graph_degopt`](@ref).

The methods are (from the paper referenced below):

  - `:ps_degopt`, Paterson–Stockmeyer method embedded into `degopt`-format.
  - `:y1s`, given by equations (34)-(35)
  - `:z1ps`, given by equations (34)-(35) and (52)
  - `:h2m`, given by equations (34)-(35) and (69)

Not all combinations of `k` and `method` are implemented. Available ones are:

  - `k<3`, `method`=`:ps_degopt`
  - `k=3`, `method`=`:y1s`, as per Table 4 in the reference
  - `k=4`, `method`=`:y1s`
  - `k=6`, `method`=`:h2m`, as per Table 11 in the reference
  - `k=8`, `method`=`:z1ps`, as per Table 7 in the reference

The default option `method=:auto` will choose method according to the value of
`k`, as prescribed above.

**Reference**

1. J. Sastre. "Efficient evaluation of matrix polynomials". Linear Algebra and
   its Applications, 539:229-250, 2018.
   DOI: [10.1016/j.laa.2017.11.010](https://doi.org/10.1016/j.laa.2017.11.010)
"""
function graph_sastre_exp(k, method = :auto)
    if (method == :auto)
        if (k < 3)
            method = :ps_degopt
        elseif (k == 3) || (k == 4)
            method = :y1s
        elseif (k == 6)
            method = :h2m
        elseif (k == 8)
            method = :z1ps
        else
            error("This number of multiplication is not available.")
        end
    end

    if (k < 3) && (method == :ps_degopt)
        if k == 2
            deg = 4
        elseif k == 1
            deg = 2
        else
            deg = 1
        end
        return graph_ps_degopt(1 ./ factorial.(0:deg))

    elseif (k == 3) && (method == :y1s) # Table 4
        e0 = 2.974307204847627
        e2 = 1.225521150112075e-1
        d1 = 8.765009801785554e-1
        d2 = 7.665265321119147e-2
        c3 = 1.992047682223989e-2
        c4 = 4.980119205559973e-3
        f0 = 1.0
        f1 = 1.0
        f2 = 0.5

        c = [c3; c4]
        d = [d1; d2]
        e = [e0; NaN; e2]

        f = [f0; f1; f2]

        s = 2
        p = 0
        return graph_sastre_z1ps_degopt(s, p, c, d, e, f, NaN)

    elseif (k == 4) && (method == :y1s) # Not tabulated?
        e0 = 5.018851944498568
        e = [e0, NaN, 0.03806343180936604, 0.017732587443103232]
        c = [0.002193172316532563, 0.0002741465395665704, 4.569108992776174e-5]
        d = [1.3093238729699403, 0.1955094199013519, 0.01626158346315151]
        f = [1.0, 1.0, 0.5, 0.1168293067115003]

        s = 3
        p = 0
        return graph_sastre_z1ps_degopt(s, p, c, d, e, f, NaN)

    elseif (k == 6) && (method == :h2m) # Table 11 and equation (69)
        e0 = 7.922322450524197
        e2 = 5.317514832355802e-2
        d1 = 2.868706220817633e-1
        d2 = -3.289442879547955e-2
        c3 = 1.052151783051235e-3
        c4 = 4.675683454147702e-4
        f0 = 3.096467971936040
        f1 = 7.723603212944010e-1
        f2 = 1.673139636901279e-1
        c = [c3; c4]
        d = [d1; d2]
        e = [e0; NaN; e2]
        f = [f0; f1; f2]

        e0p = 1.930814505527068
        e2p = 2.771400028062960e-2
        d1p = 3.968985915411500e-1
        d2p = 2.219811707032801e-2
        c3p = 2.688394980266927e-3
        c4p = 4.675683454147702e-4
        f0p = 0.0
        f1p = 2 * 1.614743005681339e-1 #Error in paper? Seems we should use 2*f1'
        f2p = 8.092036376147299e-2
        cp = [c3p; c4p]
        dp = [d1p; d2p]
        ep = [e0p; NaN; e2p]
        fp = [f0p; f1p; f2p]

        s = 2
        graph_sastre_h2m_degopt(s, (c, cp), (d, dp), (e, ep), (f, fp))

    elseif (k == 8) && (method == :z1ps) # Table 7
        c10 = -6.140022498994532e-17
        c9 = -9.210033748491798e-16
        c8 = -1.980157255925737e-14
        c7 = -4.508311519886735e-13
        c6 = -1.023660713518307e-11
        d5 = -1.227011356117036e-10
        d4 = -6.770221628797445e-9
        d3 = -1.502070379373464e-7
        d2 = -3.013961104055248e-6
        d1 = -5.893435534477677e-5
        e5 = -3.294026127901678e-10
        e4 = -2.785084196756015e-9
        e3 = -4.032817333361947e-8
        e2 = -5.100472475630675e-7
        e0 = -1.023463999572971e-3
        f5 = 4.024189993755686e-13
        f4 = 7.556768134694921e-12
        f3 = 1.305311326377090e-10
        f2 = 2.087675698786810e-9
        f1 = 2.505210838544172e-8
        f0 = 2.755731922398589e-7

        c = [c6; c7; c8; c9; c10]
        d = [d1; d2; d3; d4; d5]
        e = [e0; NaN; e2; e3; e4; e5]

        f = [f0; f1; f2; f3; f4; f5]

        a = 1 ./ factorial.(0:9)

        s = 5
        p = 10
        return graph_sastre_z1ps_degopt(s, p, c, d, e, f, a)

    else
        error("Not implemented k=$k and method=$method")
    end
end

"""
    (graph,cref)=graph_sastre_poly(b)

Computes the degree-8 polynomial

    p(z)=b[1]+z*b[2]+z^2*b[3]+...+z^8*b[9]

according to Example 3.1 in the reference.

**Reference**

1. J. Sastre. "Efficient evaluation of matrix polynomials". Linear Algebra and
   its Applications, 539:229-250, 2018. DOI:
   [10.1016/j.laa.2017.11.010](https://doi.org/10.1016/j.laa.2017.11.010)
"""
function graph_sastre_poly(b)
    # Equations (16) - (32)
    if (size(b, 1) != 9)
        error("Not implemented for length(b)=$k")
    end

    b0 = b[1]
    b1 = b[2]
    b2 = b[3]
    b3 = b[4]
    b4 = b[5]
    b5 = b[6]
    b6 = b[7]
    b7 = b[8]
    b8 = b[9]

    c4 = sqrt(b8) # plus minus?
    c3 = b7 / (2 * c4)
    d2_plus_e2 = (b6 - c3^2) / c4
    d1 = (b5 - c3 * d2_plus_e2) / c4

    e2_num_sqrt =
        (d1 - (c3 / c4) * d2_plus_e2)^2 +
        4 * (c3 / c4) * (b3 + (c3^2 / c4) * d1 - (c3 / c4) * b4)

    e2_num = (c3 / c4) * d2_plus_e2 - d1 + sqrt(e2_num_sqrt)  #plus minus
    e2 = e2_num / (2 * c3 / c4)
    d2 = d2_plus_e2 - e2
    f2 = b2
    f1 = b1
    f0 = b0

    e0 = (b3 - d1 * e2) / c3 # Not explicitly documented?

    e = [e0; NaN; e2]
    c = [c3; c4]
    f = [f0; f1; f2]
    d = [d1; d2]

    s = 2
    p = s
    return graph_sastre_z1ps_degopt(s, p, c, d, e, f, NaN)
end

"""
    (graph,cref)=graph_sastre_yks_degopt(k,s,c)

Transforms the polynomial evaluation format given by equations (62)-(65) in the
reference to `degop`-format. The `graph` is a representation of ```y_{k,s}```.
Input `c` is a grouping of the coefficients as given by the representation
(62)-(65). `c` is a `Vector{Vector{Vector}}` of length `k+1`, representing

    [
    [c_i^{(0,1)}, c_i^{(0,2)}]
    [c_i^{(1,1)}, c_i^{(1,2)}, c_i^{(1,3)}, c_i^{(1,4)}, c_i^{(1,5)}, c_i^{(1,6)}]
    ...
    [c_i^{(k,1)}, c_i^{(k,2)}, c_i^{(k,3)}, c_i^{(k,4)}, c_i^{(k,5)}, c_i^{(k,6)}]
    ]

Hence, `c[1]` contains two vectors, the first of length `s`-1 and the second of
length `s`. (Note: In the first vector the constant for I is set to zero) and
`c[2]` up to `c[k+1]` conatins six vectors: Even-numbered vectors have length
`s`+1 Odd-numbered vectors have length j-1, where j is the intex in `c`, e.g.,
`c[2][1]` has one element and `c[3][5]` have two. For example, equations
(57)-(59) are implemented as:

    c = [
    [[c15, c16], [0.0, 0, 0]], # (57)
    [[1.0], [0, c13, c14], [1.0], [c11, 0, c12], [c10], [0.0, 0, 0]], # (58)
    [[0.0, 1], [0, c8, c9], [c7, 1], [0, c6, 0], [c4, c5], [c1, c2, c3]], # (59)
    ]

**Reference**

1. J. Sastre. "Efficient evaluation of matrix polynomials". Linear Algebra and
   its Applications, 539:229-250, 2018. DOI:
   [10.1016/j.laa.2017.11.010](https://doi.org/10.1016/j.laa.2017.11.010)
"""
function graph_sastre_yks_degopt(k, s, c)
    T = eltype(eltype(eltype(c)))
    x = Vector{Tuple{Vector{T},Vector{T}}}()

    for j = 1:s-1
        push!(
            x,
            (vcat(zero(T), one(T), zeros(T, j - 1)), vcat(zeros(T, j), one(T))),
        )
    end

    # y_0s first term
    push!(
        x,
        (vcat(zeros(T, s), one(T)), vcat(zero(T), c[1][1][1:s])), #C_i^{(0,1)}
    )
    if (k == 0)
        z = vcat(c[1][2][1:s+1], one(T)) #C_i^{(0,2)}
        return graph_degopt(x, z)
    else

        # y_js = g_js + sum_i=0^(j-1) c_i^(j,5) y_is + sum_i=0^s c_i^(j,6) x^i
        # Where g_js are the result of the multiplication that is saved. By
        # expanding y_is into similar sums there is a recursion of equivalent
        # adjusted coefficients that needs to be calculated as iterations
        # proceed. Other terms have to adjusted for. Store recursion in vector
        # below.

        # Vector of adjustments of monomial coeffs
        adj_coeffs_mon = [c[1][2][1:s+1]]
        # Vector of adjustments of coeffs for pre-calculated-terms
        adj_coeffs_mult = [zeros(T, 0)]

        # y_1s first term - incorporating second term of y0s
        a = c[2][2][1:s+1] + c[2][1][1] * adj_coeffs_mon[1]
        b = c[2][4][1:s+1] + c[2][3][1] * adj_coeffs_mon[1]
        push!(x, (vcat(a, c[2][1][1]), vcat(b, c[2][3][1])))

        push!(adj_coeffs_mon, c[2][6][1:s+1] + c[2][5][1] * adj_coeffs_mon[1])
        push!(adj_coeffs_mult, [c[2][5][1]])

        for j = 2:k
            # yjs first term - incorporating second and third terms of lower yls
            # for l < j.
            a1 = c[j+1][2][1:s+1]
            b1 = c[j+1][4][1:s+1]
            c1 = c[j+1][6][1:s+1]
            for l = 1:j
                a1 += c[j+1][1][l] * adj_coeffs_mon[l]
                b1 += c[j+1][3][l] * adj_coeffs_mon[l]
                c1 += c[j+1][5][l] * adj_coeffs_mon[l]
            end

            a2 = c[j+1][1][1:j]
            b2 = c[j+1][3][1:j]
            c2 = c[j+1][5][1:j]
            for l = 2:j
                a2[1:l-1] += c[j+1][1][l] * adj_coeffs_mult[l][1:l-1]
                b2[1:l-1] += c[j+1][3][l] * adj_coeffs_mult[l][1:l-1]
                c2[1:l-1] += c[j+1][5][l] * adj_coeffs_mult[l][1:l-1]
            end

            push!(x, (vcat(a1, a2), vcat(b1, b2)))

            push!(adj_coeffs_mon, c1)
            push!(adj_coeffs_mult, c2)
        end

        z = vcat(adj_coeffs_mon[k+1], adj_coeffs_mult[k+1], one(T))
    end

    return graph_degopt(x, z)
end

# Internal use only
# Transforms formula (34)-(35) + (52) to degopt form
# Input are the coefficients given in the paper
# s int
# d=[d1,...ds]
# c=[c_(s+1)...c_(2*s)]  # size = s
# e=[e0;NaN;e2;...e_s]  # size = s+1
# f=[f0;...f_s] # size = s+1
# a=[a0;...a_(p-1)] # size = p   Can be NaN if p=s
# Nof mult: s+1+p/s = s+1+v
function graph_sastre_z1ps_degopt(s, p, c, d, e, f, a)
    T = eltype(c)
    x = Vector{Tuple{Vector{T},Vector{T}}}()

    for j = 1:s-1
        push!(
            x,
            (vcat(zero(T), one(T), zeros(T, j - 1)), vcat(zeros(T, j), one(T))),
        )
    end

    # y0s
    push!(x, (vcat(zeros(T, s), one(T)), vcat(zero(T), c[1:s])))

    # first term y1s
    push!(
        x,
        (
            vcat(zero(T), d[1:s], one(T)),
            vcat(zero(T), zero(T), e[3:(s+1)], one(T)),
        ),
    )

    # y1s
    y1s = vcat(f[1:s+1], e[1], one(T))
    if (p == 0) || (p == s)
        z = y1s
        return graph_degopt(x, z)
    else #Apply PS-scheme as in (52)
        # Assumed to be an integer, according to paper, p = v*s
        v = convert(Int, p / s)

        # Evaluate y1s*x^s
        push!(x, (y1s, vcat(zeros(T, s), one(T), zeros(T, 2))))

        # Evaluate rest of the multiplications with PS-scheme
        for i in reverse(0:v-2)
            lower_idx = (i + 1) * s + 1
            upper_idx = lower_idx + s - 1
            idx = lower_idx:upper_idx
            c = view(a, idx)
            # Coeffs for P_{v-i} + A^s*P_{v+1-i}
            xl = vcat(c, zero(T), zeros(T, v - i), one(T))
            xr = vcat(zeros(T, s), one(T), zeros(T, v + 1 - i)) # Coeffs for A^s
            push!(x, (xl, xr)) # A^s*( P_{v-i} + A^s*P_{v+1-i} )
        end

        z = vcat(a[1:s], zeros(T, 2 + v), one(T))
        return graph_degopt(x, z)
    end
end

# Internal use only
# Transforms formula (34)-(35) + (69) to degopt form
# Input is analogous to graph_sastre_z1ps_degopt() but with tuples of vectors
# # Nof mult: s+4
function graph_sastre_h2m_degopt(s, c, d, e, f)
    T = eltype(eltype(c))
    x = Vector{Tuple{Vector{T},Vector{T}}}()

    for j = 1:s-1
        push!(
            x,
            (vcat(zero(T), one(T), zeros(T, j - 1)), vcat(zeros(T, j), one(T))),
        )
    end

    # y0s
    push!(x, (vcat(zeros(T, s), one(T)), vcat(zero(T), c[1][1:s])))
    # first term y1s
    push!(
        x,
        (
            vcat(zero(T), d[1][1:s], one(T)),
            vcat(zero(T), zero(T), e[1][3:(s+1)], one(T)),
        ),
    )
    # y1s (adjusted for two more multiplications)
    y1s = vcat(f[1][1:s+1], e[1][1], one(T), zeros(T, 2))

    # y0s_p
    push!(
        x,
        (
            vcat(zeros(T, s), one(T), zeros(T, 2)),
            vcat(zero(T), c[2][1:s], zeros(T, 2)),
        ),
    )
    # first term y1s_p
    push!(
        x,
        (
            vcat(zero(T), d[2][1:s], zeros(T, 2), one(T)),
            vcat(zero(T), zero(T), e[2][3:(s+1)], zeros(T, 2), one(T)),
        ),
    )
    # y1s_p (skipping two multiplications in the middle)
    y1s_p = vcat(f[2][1:s+1], zeros(T, 2), e[2][1], one(T))

    push!(x, (vcat(y1s), y1s_p))

    z = vcat(one(T), zeros(s + 4), one(T))
    return graph_degopt(x, z)
end


text:
using LinearAlgebra

# Covers how to consider a square violation
"""
    sqrDistEuclid(pt, projFunc, p = 2)

Provide a point `pt` and a function `pt -> projFunc(pt)` and this function will
return the squared distance to that point. A `p = 2` norm is assumed, but that
can be adjusted if needed.

Primarily used to visualize square penalty constraints.
"""
function sqrDistEuclid(pt, projFunc, p = 2)
    #=
    Choose a given projection method and determine the squared distance
    at this point.
    =#
    return norm(pt - projFunc(pt), p)^2
end

# Covers Projections

# First, Lines from vectors at the origin
@doc raw"""
    projLine(pt, vec)

Project a point `pt` onto a line described by `vec`. `pt` should have the same
dimensionality as `vec`.

```math
Proj_v(p) = v \ \frac{p ⋅ v}{v ⋅ v}
```

See also: [`orthoProjLine`](@ref)
"""
function projLine(pt, vec)
    #=
    This function projects a point pt onto a line described by vec

    proj_v(p) = v * (p'v)/(v'v)
    =#
    return vec * (pt'vec) / (vec'vec)
end

@doc raw"""
    orthoProjLine(pt, vec)

Project a point `pt` onto the orthogonal of a line described  by `vec`. `pt`
should have the same dimensionality as `vec`.

```math
Proj_{v_⟂}(p) = p - Proj_v(p)
```

See also: [`projLine`](@ref)
"""
function orthoProjLine(pt, vec)
    #=
    This function projects a point pt onto the orthogonal of a line described
    by vec

    proj_v(p) = v * (p'v)/(v'v)
    proj_v⟂(p) = pt - v * (p'v)/(v'v)
    =#
    return pt - projLine(pt, vec)
end

# Second, Positive Orthant
@doc raw"""
    projPosOrth(pt)

Project a point `pt` onto the positive orthant. This is most commonly used for
constraints of the form `τ ≥ 0`.

```math
Proj_+(p) = \max.(p, 0)
```

"""
function projPosOrth(pt)
    #=
    This function is for a constraint of the form "τ ≥ 0."
    When τ_i ≥ 0, the variable τ is in the positive orthant, thus, the
    projection is simply τ_i.
    When τ_i < 0, the variable τ is not in the positive orthant, thus, the
    projection is distance to the positive orthant
    =#

    return max.(pt, 0)
end

# Third, Affine Case
@doc raw"""
    projAffineEq(a, b, x)

Project a point `x` onto an n-dimensional hyperplane. Equivalently, this is a
single affine constraint of the form `a'x = b`. Note that `a` and `x` should
have the same dimensionality and `b` is a real number. This function does
NOT handle matrix versions of the form `Ax = b`. To handle cases of this form,
split `A` and `b` by row and pass the same `x`. In particular:

Turn `Ax = b` into `A[1, :]'x = b1`, `A[2, :]'x = b2`, etc.

Important, there is a singularity at `a = 0` corresponding to a hyperplane
independent of `x`! Please avoid this.

```math
Proj_{Ax = b}(x) = x + \frac{ba}{a ⋅ a} - a \ \frac{a ⋅ x}{a ⋅ a}
```

See also [`projAffineIneq`](@ref)
"""
function projAffineEq(a, b, x)
    #=
    Projection for Affine Equality Case (a'x = b)

    We can show that the projection is given by
    proj(x) = x + (b - a'x) * a / a'a
            = x + a * [b] / a'a - a * (a'x) / a'a
    =#
    if a'a == 0
        println("Divide by zero error in projection, returning zero.")
        return 0
    end
    t1 = a * b / a'a
    # println("T1 = $t1")
    t2 = a * (a'x) / a'a
    # println("T2 = $t2")
    return x + t1 - t2
end

"""
    projAffineIneq(a, b, x)

Project a point `x` onto an n-dimensional half space. Equivalently, this is a
single affine constraint of the form `a'x ≤ b`. This function is entirely based
on [`projAffineEq`](@ref).

Please see [`projAffineEq`](@ref) for details. When the constraint is satisfied,
the projection is simply `x` itself.
"""
function projAffineIneq(a, b, x)
    #=
    Projection for Affine Inequality Case (a'x ≤ b)

    We can show that the projection is given by
    projAffineEq(x) for a'x > b
    x               for a'x ≤ b
    =#
    # println("a = $(size(a))")
    # display(a)
    # println("x = $(size(x))")
    if (a'x)[1] > b
        return projAffineEq(a, b, x)
    end
    return x
end

@doc raw"""
    projSecondOrderCone(v, s, filled = true, p = 2)

Project `(v, s)` onto the second order cone described by
```math
C = \{(x, t) ∈ R^{n + 1} \ | \ ||x||_2 ≤ t\}
```
if the cone is filled (i.e. `filled = true`). Otherwise, project onto the
boundary of the second order cone described by
```math
C = \{(x, t) ∈ R^{n + 1} \ | \ ||x||_2 = t\}
```

The parameter `p` would allow you to select a different cone, however, for now
only the second order cone (`p = 2`) is functional.

---
Three projection cases:
1. The point is below the tip. This returns 0.
2. The point is interior to the filled cone. This returns `(v, s)` in the
filled case and projects onto the boundary in the unfilled case.
3. The point is exterior to the cone, but not below the tip. This projects
onto the boundary of the cone.

```math
Proj_{Cone \ Border}(v, s) = \frac{1}{2} (1 + \frac{s}{||v||}) \ [v, ||v||]
```

"""
function projSecondOrderCone(v, s, filled = true, p = 2)
    #=
    Projection for Second-Order Cone (AKA quadratic cone or the Lorentz cone)

    The second-order cone is C = {(x, t) ∈ R^n+1 | ||x||2 ≤ t}. Using the
    2-norm. Projection onto it is given by

    proj(v, s) =
    0                                   for ||v|| ≤ -s  (Below the tip)
    (v, s)                              for ||v|| ≤ s   (In the cone)
    (1/2)(1 + s/||v||)(v, ||v||)        for ||v|| ≥ |s| (Onto the cone)

    Note that (|s| = absolute value of s)
    =#

    if norm(v, p) ≤ -s
        # println("Below Tip")
        return zeros(size([v; s]))
    elseif norm(v, p) ≤ s
        # println("Inside already")
        if filled
            return [v; s]
        end
        if !filled
            return (1/2) * (1 + s / norm(v, p)) * [v; norm(v, p)]
        end
    elseif norm(v, p) ≥ abs(s)
        # println("Outside")
        return (1/2) * (1 + s / norm(v, p)) * [v; norm(v, p)]
    end

    # The above cases should be all inclusive, but as a safety I put an error.
    error("Second Order Cone Conditions ERROR")

end

ext: jl

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
#############################################################################
# sum.jl
# Handles sums of expressions.
# All expressions and atoms are subtpyes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################

import Base.sum
export sum

### Sum Atom
struct SumAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int, Int}

    function SumAtom(x::AbstractExpr)
        children = (x,)
        return new(:sum, hash(children), children, (1, 1))
    end
end

function sign(x::SumAtom)
    return sign(x.children[1])
end

# The monotonicity
function monotonicity(x::SumAtom)
    return (Nondecreasing(),)
end

# If we have h(x) = f o g(x), the chain rule says h''(x) = g'(x)^T f''(g(x))g'(x) + f'(g(x))g''(x);
# this represents the first term
function curvature(x::SumAtom)
    return ConstVexity()
end

function evaluate(x::SumAtom)
    return sum(evaluate(x.children[1]))
end

# Suppose x was of the form
# x = Ay where A was a coefficient. Then sum(x) can also be considered
# sum(A, 1) * y
function conic_form!(x::SumAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, x)
        objective = conic_form!(x.children[1], unique_conic_forms)
        new_obj = copy(objective)
        for var in keys(new_obj)
            re = sum(new_obj[var][1], dims=1)
            im = sum(new_obj[var][2], dims=1)
            new_obj[var] = (re,im)
        end
        cache_conic_form!(unique_conic_forms, x, new_obj)
    end
    return get_conic_form(unique_conic_forms, x)
end

# Dispatch to an internal helper function that handles the dimension argument in
# the same manner as Base, with dims=: denoting a regular sum
sum(x::AbstractExpr; dims=:) = _sum(x, dims)

_sum(x::AbstractExpr, ::Colon) = SumAtom(x)

function _sum(x::AbstractExpr, dimension::Integer)
    if dimension == 1
        return Constant(ones(1, x.size[1]), Positive()) * x
    elseif dimension == 2
        return x * Constant(ones(x.size[2], 1), Positive())
    else
        error("Sum not implemented for dimension $dimension")
    end
end

Base.@deprecate sum(x::AbstractExpr, dim::Int) sum(x, dims=dim)


text:
using RMInsurance
using DataFrames
using Test
using CSV

### TESTS
curr_dir = @__DIR__

println("start test ChainLadder_Test.jl...")
# We use the example from Mack's 1993-paper
# Cumulated claim amounts

## First code cell
cum_triangle = CSV.File(curr_dir * "/ChainLadderMack_Input.csv",
                        header=false, ignoreemptyrows=true) |> DataFrame
cum_triangle = Matrix(cum_triangle)
cum_triangle = convert.(Float64, cum_triangle)

## Next code cell
triangle =cum2claims(cum_triangle)
@test cum_triangle ≈ claims2cum(triangle)

mack = Mack(cum_triangle, cum=true)

# Test: development factor
@test  [3.491, 1.747, 1.457, 1.174, 1.104, 1.086, 1.054,
        1.077,  1.018] ≈ round.(mack.f, digits=3)
# Test: future claims
@test sum(mack.futureclaims) ≈ mack.tot_res
# Test: standard error
@test [0.80, 0.26, 0.19, 0.27, 0.29, 0.26,
       0.22, 0.23, 0.29] ≈ round.(sqrt.(mack.mse) ./ mack.res, digits=2)[2:end]
# Test: total standard error
@test 0.13 ≈ round(√(mack.tot_mse)/mack.tot_res, digits = 2)

println("...end test ChainLadder_Test.jl")


text:
using Base.Test

function print_evalfile(fdir, fname)
    println("EVALUATING ", fname, ":\n*************")
    result = evalfile(joinpath(fdir, fname))
    println("*************\n")
    return result
end

fdir, this_fname = splitdir(@__FILE__)
fdir = joinpath(fdir, "../src")
cd(fdir)

for fname in readdir(".")
    if endswith(fname, ".jl") && fname != this_fname
        print_evalfile(fdir, fname)
    end
end

max_stars_repo_name: inkydragon/JuliaByExample

max_issues_repo_name: inkydragon/JuliaByExample

max_forks_repo_name: inkydragon/JuliaByExample


text:
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Heerden2013")

@time begin

    import Chemostat_Heerden2013
    const ChH = Chemostat_Heerden2013
    const Hd = ChH.HeerdenData

    using UtilsJL
    const UJL = UtilsJL

    using Plots
end

## -------------------------------------------------------------------
fileid = "1.1"
function mysavefig(p, pname; params...)
    fname = UJL.mysavefig(p, string(fileid, "_", pname), Hd.HEERDEN_FIGURES_DIR; params...)
    @info "Plotting" fname
end

## -------------------------------------------------------------------
# D vs X
let
    p = plot(; title = "Heerden",
        xlabel = string("Xv (", Hd.unit(:DCW), ")"),
        ylabel = string("D (", Hd.unit(:D), ")")
    )
    scatter!(p, Hd.val(:D), Hd.val(:DCW);
        label = "", m = 8, color = :black
    )
    mysavefig(p, "X_vs_D")
end

## -------------------------------------------------------------------
# BALANCE
let
    ps = Plots.Plot[]
    for met in [:GLC]
        p = plot(; title = string("Balance: ", met),
            xlabel = "feed",
            ylabel = "exch + drain"
        )

        exps = 1:13

        feed = Hd.cval.(met, exps) .* Hd.val.(:D, exps)
        exch = Hd.uval.(met, exps) .* Hd.val.(:DCW, exps) .|> abs
        drain = Hd.sval.(met, exps) .* Hd.val.(:D, exps) .|> abs

        scatter!(p, feed, exch .+ drain;
            label = "", m = 8, color = :black
        )
        vals = [feed; exch .+ drain] |> sort
        plot!(p, vals, vals;
            label = "", ls = :dash, lw = 3, alpha = 0.6, color = :black
        )
        push!(ps, p)
    end
    mysavefig(ps, "Balances")
end


text:
struct BackPath
    sides::NTuple{2, HorizonSide}
    path::Vector{Int}
    function BackPath(robot, sides::NTuple{2, HorizonSide})
        local path=Int[]
        while !isborder(robot, sides[1]) || !isborder(robot, sides[2])
            for s in side
                push!(path, movements!(robot, s))
            end
        end
        return new(reverse(inverse.(sides)), path)
    end
end

function back(robot, backpath::BackPath)
    i = 1
    for n in backpath.path
        movemets!(robot, backpath.sides[i])
        i = i % 2 + 1
    end

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

back_path = BackPath(r, (Sud, West))

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct DistributionDensity{D<:Distribution{Multivariate,Continuous}} <: DistLikeDensity
    d::D
end

Base.convert(::Type{AbstractDensity}, d::Distribution{Multivariate,Continuous}) =
    DistributionDensity(d)

Base.convert(::Type{DistLikeDensity}, d::Distribution{Multivariate,Continuous}) =
    DistributionDensity(d)

Base.parent(density::DistributionDensity) = density.d

function density_logval(
    density::DistributionDensity,
    params::Any
)
    Distributions.logpdf(density.d, params)
end

param_bounds(density::DistributionDensity) = NoParamBounds(length(density.d))

params_shape(density::DistributionDensity) = valshape(density.d)

Distributions.sampler(density::DistributionDensity) = bat_sampler(parent(density))

Random.Sampler(rng::AbstractRNG, density::DistributionDensity, repetition::Val{1}) = sampler(density)

Statistics.cov(density::DistributionDensity) = cov(density.d)

param_bounds(density::DistributionDensity{<:NamedTupleDist}) = param_bounds(density.d)


text:
using LightGraphs, SimpleWeightedGraphs, LinearAlgebra
"""
    SunFlowerGraph(; N = 400)

SUNFLOWERGRAPH construct a simple weighted sunflower graph with N vertices. Edge weights are the reciprocal of Euclidean distances.

# Input Arguments
- `N::Int64`: default is 400, the number of vertices. Requires N > 26.

# Output Argument
- `G::SimpleWeightedGraph{Int64,Float64}`: a simple weighted graph of the sunflower.
- `L::Matrix{Float64}`: the weighted unnormalized graph Laplacian matrix.
- `X::Matrix{Float64}`: a matrix whose i-th row represent the 2D coordinates of the i-th node.

"""
function SunFlowerGraph(; N = 400)
    c=1.0/N; θ=(sqrt(5.0)-1)*π;
    X = zeros(N,2); for k=1:N X[k,:]=c*(k-1)*[cos((k-1)*θ) sin((k-1)*θ)]; end

    G = SimpleWeightedGraph(N)
    for k = 2:8
        G.weights[1,k] = 1/norm(X[1,:] - X[k,:])
        G.weights[k,1] = 1/norm(X[1,:] - X[k,:])
    end

    for k = 1:N
        if k+8 <= N
            G.weights[k,k+8] = 1/norm(X[k,:] - X[k+8,:])
            G.weights[k+8,k] = 1/norm(X[k,:] - X[k+8,:])
        end
        if k+13 <= N
            G.weights[k,k+13] = 1/norm(X[k,:] - X[k+13,:])
            G.weights[k+13,k] = 1/norm(X[k,:] - X[k+13,:])
        end
    end
    W = weights(G) #weighted adjacency_matrix
    Lw = Diagonal(sum(W, dims = 2)[:]) - W #weighted laplacian_matrix
    return G, Lw, X
end


text:
using Documenter, TweedieDistributions

makedocs(;
    modules=[TweedieDistributions],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/jkbest2/TweedieDistributions.jl/blob/{commit}{path}#L{line}",
    sitename="TweedieDistributions.jl",
    authors="John Best",
    assets=String[],
)

deploydocs(;
    repo="github.com/jkbest2/TweedieDistributions.jl",
)


text:
# using IJulia; notebook(detached=true; dir = "/home/fritzsch/fri/JAC.jl/tutorials")
using IJulia; jupyterlab(detached=true; dir = "/home/fritzsch/fri/JAC.jl/tutorials")


text:
# samplepath.jl

"""
    abstract type AbstractSamplePath <: Any

subtype: SamplePath.

Its subtype implements a continuous samplepath.
"""
abstract type AbstractSamplePath end

"""
    isstrictlyincreasing(x::AbstractVector{T})::Bool where T <: Number
    isstrictlyincreasing(x::AbstractRange{T})::Bool where T <: Number

Internal function, not exported.

Tests whether a vector of numbers is strictly increasing.
"""
function isstrictlyincreasing(x)
  for i in 1:length(x)-1
    x[i+1] <= x[i] && return false
  end
  return true
end

isstrictlyincreasing(x::AbstractRange) = step(x) > 0

"""
    SamplePath(timeinterval::S, samplevalues::T) where
        {S<:AbstractVector{Float64}, T<:AbstractVector{Float64}}

 SamplePath implements a continuous samplepath where the timeinterval is not necessarily
 equally spaced. Sample value samplevalues[k] is the value of the process at
 time timeinterval[k]. timeinterval is a strictly increasing vector. timeinterval and
 samplevalues are of equal length.

 # Performance hint:

 Use AbstractRange objects for SamplePath.timeinterval whenever possible.

# Constructors:

```julia
    SamplePath(timeinterval, samplevalues)
    SamplePath(timeinterval, f::Function) =
        SamplePath(timeinterval, f.(timeinterval))
```

# Examples

```julia
t = 0.0:0.01:2.0
X = SamplePath(t, sinpi)
```

```julia
t = 0.0:0.01:1.0
v = map(x->floor(10*x), t)
X = SamplePath(t, v)
```
"""
struct SamplePath{S, T} <: AbstractSamplePath where {S<:AbstractVector{Float64},
        T<:AbstractVector{Float64}}
  timeinterval::S
  samplevalues::T

  function SamplePath(timeinterval::S, samplevalues::T) where {
      S<:AbstractVector{Float64}, T<:AbstractVector{Float64}}
    length(timeinterval) == length(samplevalues) ||
      throw(DimensionMismatch("Length of timeinterval should be equal to the length of the
        samplevalues vector"))
    isstrictlyincreasing(timeinterval) ||
      throw(ArgumentError("Timeinterval should be increasing"))
    new{S,T}(timeinterval, samplevalues)
  end
end

#constructor
SamplePath(timeinterval, f::Function) = SamplePath(timeinterval, f.(timeinterval))

"""
    Base.step(X::SamplePath{<:AbstractRange})

Returns the step size of SamplePath.timeinterval.

# Example

```julia
t = 0.0:0.01:2.0
X = SamplePath(t, sinpi)
step(X) == step(t)
```
"""
Base.step(X::SamplePath{<:AbstractRange}) = Base.step(X.timeinterval)

"""
    Base.length(X::SamplePath)

Returns the length of the vector timeinterval == length vector samplevalues.

# Examples
```julia
X = SamplePath([0.,1.,2.], [3.,5., -1.])
length(X) # == 3
```

```julia
X = SamplePath(0.0:0.1:2π, sin)
length(X) == length(0.0:0.1:2π)
```
"""
Base.length(X::SamplePath)=length(X.timeinterval)

# For the functions below test need to be written.

"""
    Base.iterate(X::SamplePath, state=1)

Iterate over the sample values.

# Example
```julia
X = SamplePath([1.0, 2.0, 3.0], x->x^2)
for value in X
    println(value)
end
# Or in the reverse
for value in Iterators.Reverse(X)
    println(value)
end
```
"""
Base.iterate(X::SamplePath, state=1) = state > length(X) ? nothing : (X.samplevalues[state], state + 1)
Base.iterate(rX::Base.Iterators.Reverse{SamplePath{S,T}}, state=length(rX)) where {S, T} = state < 1 ? nothing : (rX.itr.samplevalues[state], state - 1)

"""
    Base.eltype(::Type{SamplePath})

Outputs element type of SamplePath, which is Float64.
"""
Base.eltype(::Type{SamplePath{S,T}}) where {S,T} = Float64

"""
    Base.getindex(X::SamplePath, i::Int)

Returns the i-th observation of the samplepath, which is X.samplevalues[i].
"""
Base.getindex(X::SamplePath, i::Int) = X.samplevalues[i]

"""
    Base.firstindex(X::SamplePath)

The first index of the SamplePath is 1.
"""
Base.firstindex(X::SamplePath) = 1

"""
    Base.lastindex(X::SamplePath)

The last index of the SamplePath is length(X).
"""
Base.lastindex(X::SamplePath) = length(X)


text:
struct Executor
    stepdefmatcher::StepDefinitionMatcher
    presenter::RealTimePresenter

    Executor(matcher::StepDefinitionMatcher, presenter::RealTimePresenter = QuietRealTimePresenter()) = new(matcher, presenter)
end

abstract type StepExecutionResult end

struct NoStepDefinitionFound <: StepExecutionResult
    step::Gherkin.ScenarioStep
end
struct NonUniqueMatch <: StepExecutionResult
    locations::Vector{StepDefinitionLocation}
end
struct SuccessfulStepExecution <: StepExecutionResult end
struct StepFailed <: StepExecutionResult
    assertion::String
end
struct UnexpectedStepError <: StepExecutionResult
    ex::Exception
    stack::StackTrace
end
struct SkippedStep <: StepExecutionResult end

issuccess(::SuccessfulStepExecution) = true
issuccess(::StepExecutionResult) = false

struct ScenarioResult
    steps::Vector{StepExecutionResult}
    scenario::Scenario
end

function executescenario(executor::Executor, scenario::Gherkin.Scenario)
    present(executor.presenter, scenario)
    context = StepDefinitionContext()
    steps = Vector{StepExecutionResult}(length(scenario.steps))
    fill!(steps, SkippedStep())
    lastexecutedstep = 0

    # Find a unique step definition for each step and execute it.
    # Present each step, first as it is about to be executed, and then once more with the result.
    for i = 1:length(scenario.steps)
        present(executor.presenter, scenario.steps[i])
        steps[i] = try
            stepdefinition = findstepdefinition(executor.stepdefmatcher, scenario.steps[i])
            try
                Base.invokelatest(stepdefinition.definition, context)
            catch ex
                UnexpectedStepError(ex, catch_stacktrace())
            end
        catch ex
            if ex isa NoMatchingStepDefinition
                NoStepDefinitionFound(scenario.steps[i])
            elseif ex isa NonUniqueStepDefinition
                NonUniqueMatch(ex.locations)
            else
                rethrow(ex)
            end
        end
        present(executor.presenter, scenario.steps[i], steps[i])
        lastexecutedstep = i
        if !issuccess(steps[i])
            break
        end
    end

    # Present any remaining steps as skipped.
    for k = lastexecutedstep + 1:length(scenario.steps)
        present(executor.presenter, scenario.steps[k])
        present(executor.presenter, scenario.steps[k], steps[k])
    end

    ScenarioResult(steps, scenario)
end

function executescenario(executor::Executor, outline::Gherkin.ScenarioOutline)
    scenarios = transformoutline(outline)
    [executescenario(executor, scenario)
     for scenario in scenarios]
end

struct FeatureResult
    feature::Feature
    scenarioresults::Vector{ScenarioResult}
end

function executefeature(executor::Executor, feature::Gherkin.Feature)
    present(executor.presenter, feature)
    scenarioresults = [executescenario(executor, s) for s in feature.scenarios]
    FeatureResult(feature, reduce(vcat, [], scenarioresults))
end


text:
module TestDataFrames

using Base.Test
using DataArrays, DataFrames
using DataFramesMeta

df = DataFrame(A = 1:3, B = [2, 1, 2])
x = [2, 1, 0]

@test  @with(df, :A + 1)   ==  df[:A] + 1
@test  @with(df, :A + :B)  ==  df[:A] + df[:B]
@test  @with(df, :A + x)   ==  df[:A] + x
x = @with df begin
    res = 0.0
    for i in 1:length(:A)
        res += :A[i] * :B[i]
    end
    res
end
@test  x == sum(df[:A] .* df[:B])
@test  @with(df, df[:A .> 1, ^([:B, :A])]) == df[df[:A] .> 1, [:B, :A]]
@test  @with(df, DataFrame(a = :A * 2, b = :A + :B)) == DataFrame(a = df[:A] * 2, b = df[:A] + df[:B])

@test  @ix(df, :A .> 1)           == df[df[:A] .> 1,:]
@test  @ix(df, :B .> 1)           == df[df[:B] .> 1,:]
@test  @ix(df, :A .> x)           == df[df[:A] .> x,:]
@test  @ix(df, :B .> x)           == df[df[:B] .> x,:]
@test  @ix(df, :A .> :B)          == df[df[:A] .> df[:B],:]
@test  @ix(df, :A .> 1, [:B, :A]) == df[df[:A] .> 1, [:B, :A]]

@test  @where(df, :A .> 1)         == df[df[:A] .> 1,:]
@test  @where(df, :B .> 1)         == df[df[:B] .> 1,:]
@test  @where(df, :A .> x)         == df[df[:A] .> x,:]
@test  @where(df, :B .> x)         == df[df[:B] .> x,:]
@test  @where(df, :A .> :B)        == df[df[:A] .> df[:B],:]

@byrow(df, if :A > :B; :A = 0 end)
@test  df == DataFrame(A = [1, 0, 0], B = [2, 1, 2])

df = DataFrame(A = 1:3, B = [2, 1, 2])  # Restore df
y = 0
function f(x)
    @byrow(df, if :A + :B == 3; x += 1 end)
    x
end
y = f(y)
@test  y == 2

end # module


text:
using Serialization

function sleep_random()
    t = ELECTION_TIMEOUT + rand(Int)%23 * CANDIDATE_SLEEP_GRANULARITY
    sleep(t)
end

is_majority(votes) = return (votes < (length(nodes) + 1)/2 ? false : true)

"""
    get_node_info(id)

Given a node id return node info.
"""
function get_node_info(id)
    for i = 1:length(nodes)
        (nodes[i].nodeId == id) && return nodes[i]
    end
    return nothing
end

"""
    byte_array(x)

Serialize Julia type *x* such that it can be written to
a db or file
"""
function byte_array(x)
    iob = IOBuffer()
    serialize(iob, x)
    iob.data
end

struct KVSSerializeException <: Exception
    msg::String
end

"""
    array_to_type{T}(arr, ::Type{T})

Deserialize from Array{UInt8} back to Julia type.
"""
function array_to_type(arr, ::Type{T}) where T
    iob = IOBuffer(arr)
    seek(iob, 0)
    t = deserialize(iob)
    if ! isa(t, T)
        throw(KVSSerializeException("Could not deserialize to type"))
    end
    t
end

function array_to_type(arr)
    iob = IOBuffer(arr)
    seek(iob, 0)
    t = deserialize(iob)
    t
end

const DEFAUT_BASE_DIR = "/tmp"
function basedir()
    b = DEFAUT_BASE_DIR
    try
        b =  ENV["RAVANA_HOME"] * "/.ravana/"
    catch e
        b = DEFAUT_BASE_DIR * "/.ravana/"
    end

    try
        stat(b).inode == 0 && Base.mkdir(b)
    catch e
        error("Cannot access $(b)")
    end
    b
end


text:
# low-level wrappers of the CUSOLVER library

function cusolverDnCreate()
  handle = Ref{cusolverDnHandle_t}()
  cusolverDnCreate(handle)
  return handle[]
end

function cusolverSpCreate()
  handle = Ref{cusolverSpHandle_t}()
  cusolverSpCreate(handle)
  return handle[]
end

function cusolverGetProperty(property::CUDAapi.libraryPropertyType)
  value_ref = Ref{Cint}()
  cusolverGetProperty(property, value_ref)
  value_ref[]
end


text:
@testset "rule types" begin
    @testset "iterating and indexing rules" begin
        _, rule = frule(dummy_identity, 1)
        i = 0
        for r in rule
            @test r === rule
            i += 1
        end
        @test i == 1  # rules only iterate once, yielding themselves
        @test rule[1] == rule
        @test_throws BoundsError rule[2]
    end

    @testset "Rule" begin
        @testset "show" begin
            @test occursin(r"^Rule\(.*foo.*\)$", repr(Rule(function foo() 1 end)))
            @test occursin(r"^Rule\(.*identity.*\)$", repr(Rule(identity)))

            @test occursin(r"^Rule\(.*identity.*\,.*\+.*\)$", repr(Rule(identity, +)))
        end
    end

    @testset "WirtingerRule" begin
        myabs2(x) = abs2(x)

        function ChainRulesCore.frule(::typeof(myabs2), x)
            return abs2(x), AbstractRule(
                typeof(x),
                Rule(Δx -> Δx * x'),
                Rule(Δx -> Δx * x)
            )
        end

        # real input
        x = rand(Float64)
        f, _df = @inferred frule(myabs2, x)
        @test f === x^2

        df = @inferred _df(One())
        @test df === x + x

        Δ = rand(Complex{Int64})
        df = @inferred _df(Δ)
        @test df === Δ * (x + x)

        # complex input
        z = rand(Complex{Float64})
        f, _df = @inferred frule(myabs2, z)
        @test f === abs2(z)

        df = @inferred _df(One())
        @test df === Wirtinger(z', z)

        Δ = rand(Complex{Int64})
        df = @inferred _df(Δ)
        @test df === Wirtinger(Δ * z', Δ * z)
    end
end


text:
module GLAPACK
# missing or inadequate wrappers: tgsen, tgsyl
using LinearAlgebra

const liblapack = Base.liblapack_name

import LinearAlgebra.BLAS.@blasfunc

import LinearAlgebra: BlasFloat, BlasInt, LAPACKException,
    DimensionMismatch, SingularException, PosDefException, chkstride1, checksquare

import LinearAlgebra.LAPACK: chkargsok, chklapackerror, chktrans

for (tgsen, elty) in
    ((:dtgsen_, :Float64),
     (:stgsen_, :Float32))
    @eval begin
        function tgsenx!(select::AbstractVector{BlasInt},
                         S::AbstractMatrix{$elty}, T::AbstractMatrix{$elty},
                         Q::AbstractMatrix{$elty}, Z::AbstractMatrix{$elty};
                         pnorm = :Frobenius)
            chkstride1(select, S, T, Q, Z)
            n, nt, nq, nz = checksquare(S, T, Q, Z)
            if n != nt
                throw(DimensionMismatch("dimensions of S, ($n,$n), and T, ($nt,$nt), must match"))
            end
            if n != nq
                throw(DimensionMismatch("dimensions of S, ($n,$n), and Q, ($nq,$nq), must match"))
            end
            if n != nz
                throw(DimensionMismatch("dimensions of S, ($n,$n), and Z, ($nz,$nz), must match"))
            end
            lds = max(1, stride(S, 2))
            ldt = max(1, stride(T, 2))
            ldq = max(1, stride(Q, 2))
            ldz = max(1, stride(Z, 2))
            m = Ref{BlasInt}(sum(select))
            alphai = similar(T, $elty, n)
            alphar = similar(T, $elty, n)
            beta = similar(T, $elty, n)
            lwork = BlasInt(-1)
            work = Vector{$elty}(undef, 1)
            liwork = BlasInt(-1)
            iwork = Vector{BlasInt}(undef, 1)
            pl = Ref{$elty}()
            pr = Ref{$elty}()
            dif = Vector{$elty}(undef, 2)
            info = Ref{BlasInt}()
            select = convert(Array{BlasInt}, select)
            ijob = BlasInt((pnorm == :Frobenius) ? 4 : 5)
            for i = 1:2  # first call returns lwork as work[1] and liwork as iwork[1]
                ccall((@blasfunc($tgsen), liblapack), Cvoid,
                       (Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ptr{BlasInt},
                        Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                        Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
                        Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                        Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
                        Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt},
                        Ptr{BlasInt}),
                    ijob, 1, 1, select,
                    n, S, lds, T,
                    ldt, alphar, alphai, beta,
                    Q, ldq, Z, ldz,
                    m, pl, pr, dif,
                    work, lwork, iwork, liwork,
                    info)
                chklapackerror(info[])
                if i == 1 # only estimated optimal lwork, liwork
                    lwork  = 2*BlasInt(real(work[1]))
                    resize!(work, lwork)
                    liwork = BlasInt(real(iwork[1]))
                    resize!(iwork, liwork)
                end
            end
            S, T, complex.(alphar, alphai), beta, Q, Z, pl[], pr[], dif
        end
    end
end

for (tgsen, elty, relty) in
    ((:ztgsen_, :ComplexF64, :Float64),
     (:ctgsen_, :ComplexF32, :Float32))
    @eval begin
        function tgsenx!(select::AbstractVector{BlasInt},
                         S::AbstractMatrix{$elty}, T::AbstractMatrix{$elty},
                         Q::AbstractMatrix{$elty}, Z::AbstractMatrix{$elty};
                         pnorm = :Frobenius)
            chkstride1(select, S, T, Q, Z)
            n, nt, nq, nz = checksquare(S, T, Q, Z)
            if n != nt
                throw(DimensionMismatch("dimensions of S, ($n,$n), and T, ($nt,$nt), must match"))
            end
            if n != nq
                throw(DimensionMismatch("dimensions of S, ($n,$n), and Q, ($nq,$nq), must match"))
            end
            if n != nz
                throw(DimensionMismatch("dimensions of S, ($n,$n), and Z, ($nz,$nz), must match"))
            end
            lds = max(1, stride(S, 2))
            ldt = max(1, stride(T, 2))
            ldq = max(1, stride(Q, 2))
            ldz = max(1, stride(Z, 2))
            m = Ref{BlasInt}(sum(select))
            alpha = similar(T, $elty, n)
            beta = similar(T, $elty, n)
            lwork = BlasInt(-1)
            work = Vector{$elty}(undef, 1)
            liwork = BlasInt(-1)
            iwork = Vector{BlasInt}(undef, 1)
            pl = Ref{$relty}()
            pr = Ref{$relty}()
            dif = Vector{$relty}(undef, 2)
            info = Ref{BlasInt}()
            select = convert(Array{BlasInt}, select)
            ijob = BlasInt((pnorm == :Frobenius) ? 4 : 5)
            for i = 1:2  # first call returns lwork as work[1] and liwork as iwork[1]
                ccall((@blasfunc($tgsen), liblapack), Cvoid,
                      (Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ptr{BlasInt},
                       Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                       Ref{BlasInt}, Ptr{$elty}, Ptr{$elty},
                       Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                       Ptr{BlasInt}, Ptr{$relty}, Ptr{$relty}, Ptr{$relty},
                       Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt},
                       Ptr{BlasInt}),
                      ijob, 1, 1, select,
                      n, S, lds, T,
                      ldt, alpha, beta,
                      Q, ldq, Z, ldz,
                      m, pl, pr, dif,
                      work, lwork, iwork, liwork,
                      info)
                chklapackerror(info[])
                if i == 1 # only estimated optimal lwork, liwork
                    lwork  = 2*BlasInt(real(work[1]))
                    resize!(work, lwork)
                    liwork = BlasInt(real(iwork[1]))
                    resize!(iwork, liwork)
                end
            end
            if m[] != sum(select)
                @warn "surprise! tgsen returns m=$(m[]), expected $(sum(select))"
            end
            S, T, alpha, beta, Q, Z, pl[], pr[], dif
        end # func
    end # eval
end # type loop

for (fn, elty, relty) in ((:dtgsyl_, :Float64, :Float64),
                   (:stgsyl_, :Float32, :Float32),
                   (:ztgsyl_, :ComplexF64, :Float64),
                   (:ctgsyl_, :ComplexF32, :Float32))
    @eval begin
        function tgsyl!(transa::AbstractChar,
                        A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty},
                        C::AbstractMatrix{$elty}, D::AbstractMatrix{$elty},
                        E::AbstractMatrix{$elty}, F::AbstractMatrix{$elty},
                        ijob::Int=1)
            # require_one_based_indexing(A, B, C)
            chkstride1(A, B, C)
            m, n = checksquare(A, B)
            lda = max(1, stride(A, 2))
            ldb = max(1, stride(B, 2))
            m1, n1 = size(C)
            if m != m1 || n != n1
                throw(DimensionMismatch("dimensions of A, ($m,$n), and C, ($m1,$n1), must match"))
            end
            ldc = max(1, stride(C, 2))

            # require_one_based_indexing(D, E, F)
            chkstride1(D, E, F)
            m2, n2 = checksquare(D, E)
            if m != m2 || n != n2
                throw(DimensionMismatch("dimensions of (A,B), ($m,$n), and (D,E), ($m2,$n2), must match"))
            end
            ldd = max(1, stride(D, 2))
            lde = max(1, stride(E, 2))
            m1, n1 = size(F)
            if m != m1 || n != n1
                throw(DimensionMismatch("dimensions of A, ($m,$n), and F, ($m1,$n1), must match"))
            end
            ldf = max(1, stride(F, 2))

            scale = Ref{$relty}(1.0)
            dif = Ref{$relty}(0.0)
            info  = Ref{BlasInt}()
            lwork = BlasInt(-1)
            work = Vector{$elty}(undef, 1)
            # only need m+n+2 for complex but m+n+6 for real
            iwork = Vector{BlasInt}(undef, m+n+6)
            for i=1:2 # first call to get workspace
              ccall((@blasfunc($fn), liblapack), Cvoid,
                  (Ref{UInt8}, Ref{BlasInt}, # trans, ijob
                   Ref{BlasInt}, Ref{BlasInt}, # m,n
                   Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                   Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                   Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                   Ptr{$relty}, Ptr{$relty}, # scale, dif
                   Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, # work, lwork, iwork
                   Ptr{BlasInt}),
                    transa, ijob, m, n,
                    A, lda, B, ldb, C, ldc,
                    D, ldd, E, lde, F, ldf,
                    scale, dif, work, lwork, iwork, info)
                chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(real(work[1]))
                    resize!(work, lwork)
                end
            end
            C, F, scale[], dif[]
        end # func
    end # eval
end # type loop

end # module


text:
using LightXML
using PyPlot

lattice_file = ARGS[1]

lfile = parse_file(lattice_file)
graph = LightXML.root(lfile)

xmlvertices = get_elements_by_tagname(graph, "VERTEX")
x = Array{Float64}(length(xmlvertices))
y = Array{Float64}(length(xmlvertices))
for (i,v) in enumerate(xmlvertices)
  cords = map(x->parse(Float64,x), split(content(find_element(v, "COORDINATE")), " "))
  x[i] = cords[1]
  y[i] = cords[2]
end

scatter(x,y, c=1:length(x), cmap="Blues")
axis("equal")
savefig("$lattice_file.pdf")

close("all")

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
module LaborCoercion

# There is a population of mass 1 of identical (coercive) producers and a
# population of mass L < 1 of identical agents
L = Var([0, 1])

# Each producer has a project that yields x > 0 units of a consumption good if
# successful and zero units if unsuccessful
x = Var(Real⁺)

# Each producer is initially randomly matched with a worker with probability L

# A producer with productivity x who is matched with an agent chooses a level of guns
# η > 0 corresponds to the cost of guns
η = Var(Real⁺)
guncost(g) = η * χ(g)

payoff(contract, agent) = agent.ū - contract.g

# and commits to a "contract" spiecfying output-dependent wage-puinshment

struct Contract
  wʰ      # Wage for high output
  wˡ      # Wage for low output
  pʰ      # Punishment for high output
  pʰ      # Punishment for low output
end

# If the agent accepts the contract offer, then she chooses effort level a ∈ [0, 1] at private cost c(a)

function ok()
end

end


text:
module ArrayViewsAPL

using Base.Cartesian

import Base: convert, copy, eltype, getindex, length, ndims, parent, parentdims,
    parentindexes, pointer, setindex!, similar, size, stride, strides

export
    # types
    View,
    # functions
    sliceview,
    subview

typealias NonSliceIndex Union(UnitRange{Int}, StepRange{Int,Int}, Vector{Int})
typealias ViewIndex Union(Int, NonSliceIndex)

# Since there are no multidimensional range objects, we only permit 1d indexes
type View{T,N,P<:AbstractArray,I<:(ViewIndex...)} <: AbstractArray{T,N}
    parent::P
    indexes::I
    dims::NTuple{N,Int}
end

# Simple utilities
eltype{T,N,P,I}(V::View{T,N,P,I}) = T
eltype{T,N,P,I}(::Type{View{T,N,P,I}}) = T
ndims{T,N,P,I}(V::View{T,N,P,I}) = N
ndims{T,N,P,I}(::Type{View{T,N,P,I}}) = N
size(V::View) = V.dims
size(V::View, d::Integer) = d <= ndims(V) ? (@inbounds ret = V.dims[d]; ret) : 1
length(V::View) = prod(V.dims)

similar(V::View, T, dims::Dims) = similar(V.parent, T, dims)
copy(V::View) = copy!(similar(V.parent, size(V)), V)

parent(V::View) = V.parent
parentindexes(V::View) = V.indexes

## View creation

stagedfunction sliceview{T,NP}(A::AbstractArray{T,NP}, I::ViewIndex...)
    N = 0
    sizeexprs = Any[]
    for k = 1:length(I)
        i = I[k]
        if !(i <: Real)
            N += 1
            push!(sizeexprs, :(length(I[$k])))
        end
    end
    dims = :(tuple($(sizeexprs...)))
    :(ArrayViewsAPL.View{$T,$N,$A,$I}(A, I, $dims))
end

# Conventional style (drop trailing singleton dimensions, keep any other singletons)
stagedfunction subview{T,NP}(A::AbstractArray{T,NP}, I::ViewIndex...)
    sizeexprs = Any[]
    Itypes = Any[]
    Iexprs = Any[]
    N = length(I)
    while N > 0 && I[N] <: Real
        N -= 1
    end
    for k = 1:length(I)
        if k <= N
            push!(sizeexprs, :(length(I[$k])))
        end
        if k < N && I[k] <: Real
            push!(Itypes, UnitRange{Int})
            push!(Iexprs, :(int(I[$k]):int(I[$k])))
        else
            push!(Itypes, I[k])
            push!(Iexprs, :(I[$k]))
        end
    end
    dims = :(tuple($(sizeexprs...)))
    Iext = :(tuple($(Iexprs...)))
    It = tuple(Itypes...)
    :(ArrayViewsAPL.View{$T,$N,$A,$It}(A, $Iext, $dims))
end

# Constructing from another View
# This "pops" the old View and creates a more compact one
stagedfunction sliceview{T,NV,PV,IV}(V::View{T,NV,PV,IV}, I::ViewIndex...)
    N = 0
    sizeexprs = Any[]
    indexexprs = Any[]
    Itypes = Any[]
    k = 0
    for j = 1:length(IV)
        if IV[j] <: Real
            push!(indexexprs, :(V.indexes[$j]))
            push!(Itypes, IV[j])
        else
            k += 1
            if k < length(I) || k == NV || j == length(IV)
                if !(I[k] <: Real)
                    N += 1
                    push!(sizeexprs, :(length(I[$k])))
                end
                push!(indexexprs, :(V.indexes[$j][I[$k]]))
                push!(Itypes, rangetype(IV[j], I[k]))
            else
                # We have a linear index that spans more than one dimension of the parent
                N += 1
                push!(sizeexprs, :(length(I[$k])))
                push!(indexexprs, :(ArrayViewsAPL.merge_indexes(V.indexes[$j:end], size(V.parent)[$j:end], I[$k])))
                push!(Itypes, Array{Int, 1})
                break
            end
        end
    end
    for i = k+1:length(I)
        if !(I[i] <: Real)
            N += 1
            push!(sizeexprs, :(length(I[$i])))
        end
        push!(indexexprs, :(I[$i]))
        push!(Itypes, I[i])
    end
    Inew = :(tuple($(indexexprs...)))
    dims = :(tuple($(sizeexprs...)))
    It = tuple(Itypes...)
    :(ArrayViewsAPL.View{$T,$N,$PV,$It}(V.parent, $Inew, $dims))
end

stagedfunction subview{T,NV,PV,IV}(V::View{T,NV,PV,IV}, I::ViewIndex...)
    N = length(I)
    while N > 0 && I[N] <: Real
        N -= 1
    end
    sizeexprs = Any[]
    indexexprs = Any[]
    Itypes = Any[]
    k = 0
    for j = 1:length(IV)
        if IV[j] <: Real
            push!(indexexprs, :(V.indexes[$j]))
            push!(Itypes, IV[j])
        else
            k += 1
            if k <= N
                push!(sizeexprs, :(length(I[$k])))
            end
            if k < N && I[k] <: Real
                # convert scalar to a range
                push!(indexexprs, :(V.indexes[$j][int(I[$k]):int(I[$k])]))
                push!(Itypes, rangetype(IV[j], UnitRange{Int}))
            elseif k < length(I) || j == length(IV)
                # simple indexing
                push!(indexexprs, :(V.indexes[$j][I[$k]]))
                push!(Itypes, rangetype(IV[j], I[k]))
            else
                # We have a linear index that spans more than one dimension of the parent
                push!(indexexprs, :(ArrayViewsAPL.merge_indexes(V.indexes[$j:end], size(V.parent)[$j:end], I[$k])))
                push!(Itypes, Array{Int, 1})
                break
            end
        end
    end
    for i = k+1:length(I)
        if i <= N
            push!(sizeexprs, :(length(I[$i])))
        end
        push!(indexexprs, :(I[$i]))
        push!(Itypes, I[i])
    end
    Inew = :(tuple($(indexexprs...)))
    dims = :(tuple($(sizeexprs...)))
    It = tuple(Itypes...)
    :(ArrayViewsAPL.View{$T,$N,$PV,$It}(V.parent, $Inew, $dims))
end

function rangetype(T1, T2)
    rt = Base.return_types(getindex, (T1, T2))
    length(rt) == 1 || error("Can't infer return type")
    rt[1]
end

subview(A::AbstractArray, I::Union(ViewIndex, Colon)...) = subview(A, ntuple(length(I), i-> isa(I[i], Colon) ? (1:size(A,i)) : I[i])...)
sliceview(A::AbstractArray, I::Union(ViewIndex, Colon)...) = sliceview(A, ntuple(length(I), i-> isa(I[i], Colon) ? (1:size(A,i)) : I[i])...)

## Strides
stagedfunction strides(V::View)
    T,N,P,I = V.parameters
    all(map(x->x<:RangeIndex, I)) || error("strides valid only for RangeIndex indexing")
    strideexprs = Array(Any, N+1)
    strideexprs[1] = 1
    i = 1
    Vdim = 1
    for i = 1:length(I)
        if !(I[i]==Int)
            strideexprs[Vdim+1] = copy(strideexprs[Vdim])
            strideexprs[Vdim] = :(step(V.indexes[$i])*$(strideexprs[Vdim]))
            Vdim += 1
        end
        strideexprs[Vdim] = :(size(V.parent, $i) * $(strideexprs[Vdim]))
    end
    :(tuple($(strideexprs[1:N]...)))
end

stride(V::View, d::Integer) = d <= ndims(V) ? strides(V)[d] : strides(V)[end] * size(V)[end]

## Pointer conversion (for ccall)
function first_index(V::View)
    f = 1
    s = 1
    for i = 1:length(V.indexes)
        f += (first(V.indexes[i])-1)*s
        s *= size(V.parent, i)
    end
    f
end

convert{T,N,P<:Array,I<:(RangeIndex...)}(::Type{Ptr{T}}, V::View{T,N,P,I}) =
    pointer(V.parent) + (first_index(V)-1)*sizeof(T)

convert{T,N,P<:Array,I<:(RangeIndex...)}(::Type{Ptr{Void}}, V::View{T,N,P,I}) =
    convert(Ptr{Void}, convert(Ptr{T}, V))

pointer(V::View, i::Int) = pointer(V, ind2sub(size(V), i))

function pointer{T,N,P<:Array,I<:(RangeIndex...)}(V::View{T,N,P,I}, is::(Int...))
    index = first_index(V)
    strds = strides(V)
    for d = 1:length(is)
        index += (is[d]-1)*strds[d]
    end
    return pointer(V.parent, index)
end

## Convert
convert{T,S,N}(::Type{Array{T,N}}, V::View{S,N}) = copy!(Array(T, size(V)), V)

## Scalar indexing
# Low dimensions: avoid splatting
vars = Expr[]
varsInt = Expr[]
varsOther = Expr[]
vars_toindex = Expr[]
for i = 1:4
    sym = symbol(string("i",i))
    push!(vars, Expr(:quote, sym))
    push!(varsInt, :($sym::Int))
    push!(varsOther, :($sym::Union(Real, AbstractVector)))
    push!(vars_toindex, :(Base.to_index($sym)))
    ex = i == 1 ? quote
         getindex(V::View, $sym::Real) = getindex(V, Base.to_index($sym))
        setindex!(V::View, v, $sym::Real) = setindex!(V, v, Base.to_index($sym))
         getindex(V::View, $sym::AbstractVector{Bool}) = getindex(V, Base.to_index($sym))
        setindex!(V::View, v, $sym::AbstractVector{Bool}) = setindex!(V, v, Base.to_index($sym))
    end : quote
         getindex(V::View, $(varsOther...)) = getindex(V::View, $(vars_toindex...))
        setindex!(V::View, v, $(varsOther...)) = setindex!(V::View, v, $(vars_toindex...))
    end
    @eval begin
        stagedfunction getindex(V::View, $(varsInt...))
            T, N, P, IV = V.parameters
            exhead, ex = index_generate(ndims(P), IV, :V, [$(vars...)])
            quote
                $exhead
                $ex
            end
        end
        stagedfunction setindex!(V::View, v, $(varsInt...))
            T, N, P, IV = V.parameters
            exhead, ex = index_generate(ndims(P), IV, :V, [$(vars...)])
            quote
                $exhead
                $ex = v
            end
        end
        $ex
    end
end
# V[] notation (extracts the first element)
stagedfunction getindex(V::View)
    T, N, P, IV = V.parameters
    Isyms = ones(Int, N)
    exhead, ex = index_generate(ndims(P), IV, :V, Isyms)
    quote
        $exhead
        $ex
    end
end
# Splatting variants
stagedfunction getindex(V::View, I::Int...)
    T, N, P, IV = V.parameters
    Isyms = [:(I[$d]) for d = 1:length(I)]
    exhead, ex = index_generate(ndims(P), IV, :V, Isyms)
    quote
        $exhead
        $ex
    end
end
stagedfunction setindex!(V::View, v, I::Int...)
    T, N, P, IV = V.parameters
    Isyms = [:(I[$d]) for d = 1:length(I)]
    exhead, ex = index_generate(ndims(P), IV, :V, Isyms)
    quote
        $exhead
        $ex = v
    end
end

# Indexing with non-scalars. For now, this returns a copy, but changing that
# is just a matter of deleting the explicit call to copy.
getindex(V::View, I::ViewIndex...) = copy(subview(V, I...))
getindex{T,N}(V::View{T,N}, I::AbstractArray{Bool,N}) = copy(subview(V, find(I)))   # this could be much better optimized
getindex{T,N}(V::View{T,N}, I::Union(Real, AbstractVector)...) = getindex(V, Base.to_index(I)...)

function setindex!{T}(V::View{T,1}, v, I::AbstractArray{Bool,1})
    length(I) == length(V) || throw(DimensionMismatch("logical vector must match array length"))
    setindex!(V, v, Base.to_index(I))
end
function setindex!{T,N}(V::View{T,N}, v, I::AbstractArray{Bool,1})
    length(I) == length(V) || throw(DimensionMismatch("logical vector must match array length"))
    setindex!(V, v, Base.to_index(I))
end
function setindex!{T,N}(V::View{T,N}, v, I::AbstractArray{Bool,N})
    size(I) == size(V) || throw(DimensionMismatch("size of Boolean mask must match array size"))
    _setindex!(V, v, find(I))  # this could be better optimized
end
setindex!{T,N}(V::View{T,N}, v, I::Union(Real,AbstractVector)...) = setindex!(V, v, Base.to_index(I)...)
setindex!(V::View, x, J::Union(Int,AbstractVector)...) = _setindex!(V, x, J...)
stagedfunction _setindex!(V::View, x, J::Union(Real,AbstractVector)...)
    gen_setindex_body(length(J))
end

function gen_setindex_body(N::Int)
    quote
        Base.Cartesian.@nexprs $N d->(J_d = J[d])
        Base.Cartesian.@ncall $N checkbounds V J
        Base.Cartesian.@nexprs $N d->(I_d = Base.to_index(J_d))
        if !isa(x, AbstractArray)
            Base.Cartesian.@nloops $N i d->(1:length(I_d)) d->(@inbounds j_d = Base.unsafe_getindex(I_d, i_d)) begin
                @inbounds (Base.Cartesian.@nref $N V j) = x
            end
        else
            X = x
            Base.Cartesian.@ncall $N Base.setindex_shape_check X I
            k = 1
            Base.Cartesian.@nloops $N i d->(1:length(I_d)) d->(@inbounds j_d = Base.unsafe_getindex(I_d, i_d)) begin
                @inbounds (Base.Cartesian.@nref $N V j) = X[k]
                k += 1
            end
        end
        V
    end
end

# NP is parent dimensionality, Itypes is the tuple typeof(V.indexes)
# NP may not be equal to length(Itypes), because a view of a 2d matrix A
# can be constructed as V = A[5:13] or as V = A[2:4, 1:3, 1].
function index_generate(NP, Itypes, Vsym, Isyms)
    if isempty(Isyms)
        Isyms = Any[1]  # this handles the syntax getindex(V)
    end
    exhead = :nothing
    NV = 0
    for I in Itypes
        NV += !(I == Int)
    end
    if length(Isyms) < NV
        # Linear indexing in the last index
        n = NV - length(Isyms)
        m = length(Isyms)
        strides = [gensym() for i = 1:n]
        indexes = [gensym() for i = 1:n+1]
        resid = gensym()
        linblock = Array(Expr, 2n+2)
        linblock[1] = :($(strides[1]) = size($Vsym, $m))
        for k = 2:n
            m += 1
            linblock[k] = :($(strides[k]) = $(strides[k-1]) * size($Vsym, $m))
        end
        k = n+1
        linblock[k] = :($resid = $(Isyms[end])-1)
        for i = n:-1:1
            k += 1
            linblock[k] = quote
                $(indexes[i+1]), $resid = divrem($resid, $(strides[i]))
                $(indexes[i+1]) += 1
            end
        end
        linblock[end] = :($(indexes[1]) = $resid+1)
        exhead = Expr(:block, linblock...)
        pop!(Isyms)
        append!(Isyms, indexes)
    end
    L = length(Itypes)
    indexexprs = Array(Any, L)
    j = 0
    for i = 1:L
        if Itypes[i] <: Real # && L-i+1 > length(Isyms)-j  # consume Isyms if we're running out of V.indexes
            indexexprs[i] = :($Vsym.indexes[$i])
        else
            j += 1
            indexexprs[i] = :(ArrayViewsAPL.unsafe_getindex($Vsym.indexes[$i], $(Isyms[j])))  # TODO: make Range bounds-checking respect @inbounds
        end
    end
    # Append any extra indexes. Must be trailing 1s or it will cause a BoundsError.
    if L < NP && j < length(Isyms)
        # This view was created as V = A[5:13], so appending them would generate interpretive confusion.
        # Instead, use double-indexing, i.e., A[indexes1...][indexes2...], where indexes2 contains the leftovers.
        return exhead, :($Vsym.parent[$(indexexprs...)][$(Isyms[j+1:end]...)])
    end
    for k = j+1:length(Isyms)
        push!(indexexprs, Isyms[k])
    end
    exhead, :($Vsym.parent[$(indexexprs...)])
end

unsafe_getindex(v::Real, ind::Int) = v
unsafe_getindex(v::Range, ind::Int) = first(v) + (ind-1)*step(v)
unsafe_getindex(v::AbstractArray, ind::Int) = v[ind]
unsafe_getindex(v, ind::Real) = unsafe_getindex(v, Base.to_index(ind))

## Merging indexes
# A view created like V = A[2:3:8, 5:2:17] can later be indexed as V[2:7],
# creating a new 1d view.
# In such cases we have to collapse the 2d space spanned by the ranges.
#
# API:
#    merge_indexes(indexes::NTuple, dims::Dims, linindex)
# where dims encodes the trailing dimensions of the parent array,
# indexes encodes the view's trailing indexes into the parent array,
# and linindex encodes the subset of these elements that we'll select.
#
# The generic algorithm makes use of div to convert elements
# of linindex into a cartesian index into indexes, looks up
# the corresponding cartesian index into the parent, and then uses
# dims to convert back to a linear index into the parent array.
#
# However, a common case is linindex::UnitRange.
# Since div is slow and in(j::Int, linindex::UnitRange) is fast,
# it can be much faster to generate all possibilities and
# then test whether the corresponding linear index is in linindex.
# One exception occurs when only a small subset of the total
# is desired, in which case we fall back to the div-based algorithm.
stagedfunction merge_indexes(indexes::NTuple, dims::Dims, linindex::UnitRange{Int})
    N = length(indexes)
    N > 0 || error("Cannot merge empty indexes")
    quote
        n = length(linindex)
        Base.Cartesian.@nexprs $N d->(I_d = indexes[d])
        L = 1
        Base.Cartesian.@nexprs $N d->(L *= length(I_d))
        if n < 0.1L   # this has not been carefully tuned
            return merge_indexes_div(indexes, dims, linindex)
        end
        Pstride_1 = 1   # parent strides
        Base.Cartesian.@nexprs $(N-1) d->(Pstride_{d+1} = Pstride_d*dims[d])
        Istride_1 = 1   # indexes strides
        Base.Cartesian.@nexprs $(N-1) d->(Istride_{d+1} = Istride_d*length(I_d))
        Base.Cartesian.@nexprs $N d->(counter_d = 1) # counter_0 is a linear index into indexes
        Base.Cartesian.@nexprs $N d->(offset_d = 1)  # offset_0 is a linear index into parent
        k = 0
        index = Array(Int, n)
        Base.Cartesian.@nloops $N i d->(1:length(I_d)) d->(offset_{d-1} = offset_d + (I_d[i_d]-1)*Pstride_d; counter_{d-1} = counter_d + (i_d-1)*Istride_d) begin
            if in(counter_0, linindex)
                index[k+=1] = offset_0
            end
        end
        index
    end
end
merge_indexes(indexes::NTuple, dims::Dims, linindex) = merge_indexes_div(indexes, dims, linindex)

# This could be written as a regular function, but performance
# will be better using Cartesian macros to avoid the heap and
# an extra loop.
stagedfunction merge_indexes_div(indexes::NTuple, dims::Dims, linindex)
    N = length(indexes)
    N > 0 || error("Cannot merge empty indexes")
    Istride_N = symbol("Istride_$N")
    quote
        Base.Cartesian.@nexprs $N d->(I_d = indexes[d])
        Pstride_1 = 1   # parent strides
        Base.Cartesian.@nexprs $(N-1) d->(Pstride_{d+1} = Pstride_d*dims[d])
        Istride_1 = 1   # indexes strides
        Base.Cartesian.@nexprs $(N-1) d->(Istride_{d+1} = Istride_d*length(I_d))
        n = length(linindex)
        L = $(Istride_N) * length(indexes[end])
        index = Array(Int, n)
        for i = 1:n
            k = linindex[i] # k is the indexes-centered linear index
            1 <= k <= L || throw(BoundsError())
            k -= 1
            j = 0  # j will be the new parent-centered linear index
            Base.Cartesian.@nexprs $N d->(d < $N ?
                begin
                    c, k = divrem(k, Istride_{$N-d+1})
                    j += (ArrayViewsAPL.unsafe_getindex(I_{$N-d+1}, c+1)-1)*Pstride_{$N-d+1}
                end : begin
                    j += ArrayViewsAPL.unsafe_getindex(I_1, k+1)
                end)
            index[i] = j
        end
        index
    end
end

## Compatability
# deprecate?
function parentdims(s::View)
    nd = ndims(s)
    dimindex = Array(Int, nd)
    sp = strides(s.parent)
    sv = strides(s)
    j = 1
    for i = 1:ndims(s.parent)
        r = s.indexes[i]
        if j <= nd && (isa(r,Range) ? sp[i]*step(r) : sp[i]) == sv[j]
            dimindex[j] = i
            j += 1
        end
    end
    dimindex
end

end


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

using GLMakie, VideoIO, Makie, Dates, Printf, Colors, ImageTransformations, FFMPEG

mss = [Millisecond(25), Millisecond(125)] # intervals between frames for the 2 different FPS
fs = Int[] # number of frames per FPS
ts = [Millisecond(0)] # the timestamps
for (i, ms) in enumerate(mss)
    n = 1
    while last(ts) < Second(5i)
        push!(ts, last(ts) + ms)
        n += 1
    end
    push!(fs, n)
end
fs[1] -= 1 # a correction for how `mp4fpsmod` works.
dar = 2 # or PAR, the physical geometry of the pixels.
w = length(ts) # width of the image
h = dar*w # scaling the height
frame = Node(1)
point = lift(frame) do fn
    [(fn, dar*fn)]
end
msg = lift(frame, point) do fn, p
    ti = @sprintf "%03i" fn
    t = @sprintf "%-012s" Time(0) + ts[fn]
    x, y = p[]
    """
    Frame #: $ti
    Time: $t
    x: $x
    y: $y"""
end
limits = FRect3D(Vec3f0(1, 2, 0), Vec3f0(w - 1, 2w - 2, 0))
scene = Scene(resolution = (w, h), limits = limits,  scale_plot = false, show_axis = false)#, padding=(0,0), backgroundcolor=:red)
text!(scene, msg, position = (2, w), align = (:left,  :bottom))#, color = :white)
scatter!(scene, point, markersize = 10, marker = :+)#, color = :white)
imgstack = []
imgstack = map(1:w) do i
    frame[] = i
    AbstractPlotting.update_limits!(scene)
    AbstractPlotting.update!(scene)
    imresize(Gray.(GLMakie.scene2image(scene)), h, h)
end
props = [:color_range=>2, :priv_data => ("crf"=>"0","preset"=>"medium")]
y = "cfr.mp4"
encodevideo(y,imgstack,framerate=30,AVCodecContextProperties=props)
f1, f2 = fs
fps1, fps2 = Millisecond(Second(1)) .÷ mss
run(`mp4fpsmod -r $f1:$fps1  -r $f2:$fps2 -o vfr.mp4 cfr.mp4`)
ffmpeg_exe(` -y -i vfr.mp4 -c copy -bsf:v h264_metadata=sample_aspect_ratio=0.5 testvideo.mp4`)

f = VideoIO.openvideo("testvideo.mp4")
pixelaspectratio = VideoIO.aspect_ratio(f)
@assert pixelaspectratio == 0.5

f = VideoIO.openvideo("testvideo.mp4")
ts2 = Millisecond[]
while !eof(f)
    read(f)
    l = gettime(f)
    push!(ts2, Millisecond(round(Int, 1000l)))
    # push!(ts2, Millisecond(1))
end
@assert ts2 == ts

rm("cfr.mp4")
rm("vfr.mp4")