text:
@doc raw"""
    EssentialManifold <: AbstractPowerManifold{ℝ}

The essential manifold is the space of the essential matrices which is represented as
a quotient space of the [`Rotations`](@ref) manifold product ``\mathrm{SO}(3)^2``.

Let ``R_x(θ), R_y(θ), R_x(θ) \in ℝ^{x\times 3}`` denote the rotation around the ``z``,
``y``, and ``x`` axis in ``ℝ^3``, respectively, and further the groups

```math
H_z = \bigl\{(R_z(θ),R_z(θ))\ \big|\ θ ∈ [-π,π) \bigr\}
```

and

```math
H_π = \bigl\{ (I,I), (R_x(π), R_x(π)), (I,R_z(π)), (R_x(π), R_y(π))  \bigr\}
```

acting elementwise on the left from ``\mathrm{SO}(3)^2`` (component wise).

Then the unsigned Essential manifold ``\mathcal{M}_{\text{E}}`` can be identified
with the quotient space

````math
\mathcal{M}_{\text{E}} := (\text{SO}(3)×\text{SO}(3))/(H_z × H_π),
````

and for the signed Essential manifold $\mathcal{M}_{\text{Ǝ}}$, the quotient reads

````math
\mathcal{M}_{\text{Ǝ}} := (\text{SO}(3)×\text{SO}(3))/(H_z).
````

An essential matrix is defined as

````math
E = (R'_1)^T [T'_2 - T'_1]_{×} R'_2,
````

where the poses of two cameras ``(R_i', T_i'), i=1,2``, are contained in the space of
rigid body transformations $SE(3)$ and the operator $[⋅]_{×}\colon ℝ^3 \to \operatorname{SkewSym}(3)$
denotes the matrix representation of the cross product operator. For more details see [^TronDaniilidis2017].

# Constructor
    EssentialManifold(is_signed=true)

Generate the manifold of essential matrices, either the signed (`is_signed=true`) or
unsigned (`is_signed=false`) variant.

[^TronDaniilidis2017]:
    > Tron R.; Daniilidis K.; The Space of Essential Matrices as a Riemannian Quotient
    > Manifold.
    > SIAM Journal on Imaging Sciences (2017),
    > DOI: [10.1137/16M1091332](https://doi.org/10.1137/16M1091332),
    > PDF: [https://www.cis.upenn.edu/~kostas/mypub.dir/tron17siam.pdf](https://www.cis.upenn.edu/~kostas/mypub.dir/tron17siam.pdf).
"""
struct EssentialManifold <: AbstractPowerManifold{ℝ,Rotations{3},NestedPowerRepresentation}
    is_signed::Bool
    manifold::Rotations{3}
end

EssentialManifold(is_signed::Bool=true) = EssentialManifold(is_signed, Rotations(3))

@doc raw"""
    check_manifold_point(M::EssentialManifold, p; kwargs...)

Check whether the matrix is a valid point on the [`EssentialManifold`](@ref) `M`,
i.e. a 2-element array containing SO(3) matrices.
"""
function check_manifold_point(M::EssentialManifold, p; kwargs...)
    if length(p) != 2
        return DomainError(
            length(p),
            "The point $(p) does not lie on $M, since it does not contain exactly two elements.",
        )
    end
    return check_manifold_point(
        PowerManifold(M.manifold, NestedPowerRepresentation(), 2),
        p;
        kwargs...,
    )
end

"""
    check_tangent_vector(M::EssentialManifold, p, X; check_base_point = true, kwargs... )

Check whether `X` is a tangent vector to manifold point `p` on the [`EssentialManifold`](@ref) `M`,
i.e. `X` has to be a 2-element array of `3`-by-`3` skew-symmetric matrices.
The optional parameter `check_base_point` indicates, whether to call [`check_manifold_point`](@ref)  for `p`.
"""
function check_tangent_vector(M::EssentialManifold, p, X; check_base_point=true, kwargs...)
    if check_base_point
        mpe = check_manifold_point(M, p; kwargs...)
        mpe === nothing || return mpe
    end
    if length(X) != 2
        return DomainError(
            length(X),
            "$(X) is not a tangent vector to the manifold $M, since it does not contain exactly two elements.",
        )
    end
    return check_tangent_vector(
        PowerManifold(M.manifold, NestedPowerRepresentation(), 2),
        p,
        X;
        check_base_point=check_base_point,
        kwargs...,
    )
end

@doc raw"""
    distance(M::EssentialManifold, p, q)

Compute the Riemannian distance between the two points `p` and `q` on the [`EssentialManifold`](@ref). This is done by
computing the distance of the equivalence classes $[p]$ and $[q]$ of the points
$p=(R_{p_1},R_{p_2}), q=(R_{q_1},R_{q_2}) ∈ SO(3)^2$, respectively. Two points in $SO(3)^2$ are equivalent iff their
corresponding essential matrices, given by
````math
E = R_1^T [e_z]_{×}R_2,
````
are equal (up to a sign flip). Using the logarithmic map, the distance is given by
````math
\text{dist}([p],[q]) = \| \text{log}_{[p]} [q] \| = \| \log_p (S_z(t_{\text{opt}})q) \|,
````
where $S_z ∈ H_z = \{(R_z(θ),R_z(θ))\colon θ \in [-π,π) \}$ in which $R_z(θ)$ is the rotation around the z axis with angle
$θ$ and $t_{\text{opt}}$ is the minimizer of the cost function
````math
f(t) = f_1 + f_2, \quad f_i = \frac{1}{2} θ^2_i(t), \quad θ_i(t)=d(R_{p_i},R_z(t)R_{b_i}) \text{ for } i=1,2,
````
where $d(⋅,⋅)$ is the distance function in $SO(3)$[^TronDaniilidis2017].
"""
distance(M::EssentialManifold, p, q) = norm(M, p, log(M, p, q))

@doc raw"""
    exp(M::EssentialManifold, p, X)

Compute the exponential map on the [`EssentialManifold`](@ref) from `p` into direction
`X`, i.e.
````math
\text{exp}_p(X) =\text{exp}_g( \tilde X),  \quad g \in \text(SO)(3)^2,
````
where $\tilde X$ is the horizontal lift of $X$[^TronDaniilidis2017].
"""
exp(::EssentialManifold, ::Any...)

function exp!(M::EssentialManifold, q, p, X)
    exp!.(Ref(M.manifold), q, p, X)
    return q
end

get_iterator(::EssentialManifold) = Base.OneTo(2)

function isapprox(M::EssentialManifold, p, q; kwargs...)
    return isapprox(distance(M, p, q), 0.0; kwargs...)
end

@doc raw"""
    log(M::EssentialManifold, p, q)

Compute the logarithmic map on the [`EssentialManifold`](@ref) `M`, i.e. the tangent vector,
whose geodesic starting from `p` reaches `q` after time 1. Here, $p=(R_{p_1},R_{p_2})$ and
$q=(R_{q_1},R_{q_2})$ are elements of $SO(3)^2$. We use that any essential matrix can, up to
scale, be decomposed to
````math
E = R_1^T [e_z]_{×}R_2,
````
where $(R_1,R_2)∈SO(3)^2$. Two points in $SO(3)^2$ are equivalent iff their corresponding
essential matrices are equal (up to a sign flip).
To compute the logarithm, we first move `q` to another representative of its equivalence class.
For this, we find $t= t_{\text{opt}}$ for which the function
````math
f(t) = f_1 + f_2, \quad f_i = \frac{1}{2} θ^2_i(t), \quad θ_i(t)=d(R_{p_i},R_z(t)R_{b_i}) \text{ for } i=1,2,
````
where $d(⋅,⋅)$ is the distance function in $SO(3)$, is minimized. Further, the group $H_z$ acting
on the left on $SO(3)^2$ is defined as
````math
H_z = \{(R_z(θ),R_z(θ))\colon θ \in [-π,π) \},
````
where $R_z(θ)$ is the rotation around the z axis with angle $θ$. Points in $H_z$ are denoted by
$S_z$. Then, the logarithm is defined
as
````math
\log_p (S_z(t_{\text{opt}})q) = [\text{Log}(R_{p_i}^T R_z(t_{\text{opt}})R_{b_i})]_{i=1,2},
````
where $\text{Log}$ is the [`logarithm`](@ref log(::Rotations, ::Any...)) on $SO(3)$. For more
details see [^TronDaniilidis2017].
"""
log(M::EssentialManifold, ::Any, ::Any)

function log!(M::EssentialManifold, X, p, q)
    # compute the closest representative of q
    t = 0
    f_min = Inf
    q2 = deepcopy(q)

    if !M.is_signed
        q2min = deepcopy(q)
        for k in 1:4
            #flip sign in q to get another member of its equivalence class
            if k == 2
                q2[1][2:3, :] = -q[1][2:3, :]
                q2[2][[1 3], :] = -q[2][[1 3], :]
            elseif k == 3
                q2[2][1:2, :] = -q[2][1:2, :]
            elseif k == 4
                q2[1][2:3, :] = -q[1][2:3, :]
                q2[2][2:3, :] = -q[2][2:3, :]
            end
            t_temp, f_temp = dist_min_angle_pair(p, q2)
            if f_temp < f_min
                q2min = deepcopy(q2)
                f_min = f_temp
                t = t_temp
            end
        end
        q2 = q2min
    else
        t, _ = dist_min_angle_pair(p, q2)
    end
    Rz = [cos(t) -sin(t) 0; sin(t) cos(t) 0; 0 0 1]
    representative_q = zeros(size(q2))
    representative_q = [Rz * q2[1], Rz * q2[2]]

    # use the logarithmic map in SO(3)^2 and the corresponding geodesic between p and representative_q to compute the solution
    log!.(Ref(M.manifold), X, p, representative_q)
    return X
end

@doc raw"""
    dist_min_angle_pair(p, q)

This function computes the global minimizer of the function
````math
f(t) = f_1 + f_2, \quad f_i = \frac{1}{2} θ^2_i(t), \quad θ_i(t)=d(R_{p_i},R_z(t)R_{b_i}) \text{ for } i=1,2,
````
for the given values. This is done by finding the discontinuity points $t_{d_i}, i=1,2$ of its derivative
and using Newton's method to minimize the function over the intervals $[t_{d_1},t_{d_2}]$ and $[t_{d_2},t_{d_1}+2π]$
separately. Then, the minimizer for which $f$ is minimal is chosen and given back together with the minimal value.
For more details see Algorithm 1 in [^TronDaniilidis2017].
"""
function dist_min_angle_pair(p, q)
    #compute rotations
    q211 = q[1] * p[1]'
    q212 = q[2] * p[2]'

    #compute values from Proposition 9
    t_break1, c1, m1, Φ1 = dist_min_angle_pair_discontinuity_distance(q211)
    t_break2, c2, m2, Φ2 = dist_min_angle_pair_discontinuity_distance(q212)

    #check if cost is constant
    tol_m_zero = 5 * eps(eltype(q211)) # i.e. promoted from p and q
    tol_break = sqrt(eps(eltype(q211)))
    if (abs(m1) < tol_m_zero) && (abs(m2) < tol_m_zero)
        t_min = 0
        f_min = 2 * pi^2
    else
        if abs(mod(t_break1 - t_break2 + pi, 2 * pi) - pi) < tol_break
            t_min = t_break1 + pi
            f_min = 0
        else
            t_search1 = t_break1
            t_search2 = t_break2
            if t_search1 > t_search2
                t_search1 = t_search1 - 2 * pi
            end

            #compute derivatives of each term at discontinuity points
            df1_break1 = dist_min_angle_pair_compute_df_break(t_break1, q211)
            df2_break2 = dist_min_angle_pair_compute_df_break(t_break2, q212)

            #compute derivative of each term at other's discontinuity
            θ1_break2 = acos(clamp((m1 * sin(t_break2 + Φ1) + c1 - 1) / 2, -1.0, 1.0))
            df1_break2 = -θ1_break2 * (m1 * cos(t_break2 + Φ1)) / (2 * sin(θ1_break2))
            θ2_break1 = acos(clamp((m2 * sin(t_break1 + Φ2) + c2 - 1) / 2, -1.0, 1.0))
            df2_break1 = -θ2_break1 * (m2 * cos(t_break1 + Φ2)) / (2 * sin(θ2_break1))

            #compute left and right derivatives of sum of the two terms
            df_break1n = df1_break1 + df2_break1
            df_break1p = -df1_break1 + df2_break1
            df_break2n = df2_break2 + df1_break2
            df_break2p = -df2_break2 + df1_break2

            f_min = Inf
            if sign(df_break1p) != sign(df_break2n)
                #parabolic prediction of min
                t_min0 =
                    t_search1 -
                    df_break1p * (t_search2 - t_search1) / (df_break2n - df_break1p)
                #use Newton's search to find t_min
                t_min, f_min = dist_min_angle_pair_df_newton(
                    m1,
                    Φ1,
                    c1,
                    m2,
                    Φ2,
                    c2,
                    t_min0,
                    t_search1,
                    t_search2,
                )
            end
            t_search1 = t_search1 + 2 * pi
            if sign(df_break2p) != sign(df_break1n)
                #parabolic prediction of min
                t_min0 =
                    t_search2 -
                    df_break2p * (t_search1 - t_search2) / (df_break1n - df_break2p)
                #use Newton's search to find t_min
                t_min2, f_min2 = dist_min_angle_pair_df_newton(
                    m1,
                    Φ1,
                    c1,
                    m2,
                    Φ2,
                    c2,
                    t_min0,
                    t_search2,
                    t_search1,
                )
                #choose t with lower function values
                if f_min2 < f_min
                    t_min = t_min2
                    f_min = f_min2
                end
            end
        end
    end
    return t_min, f_min
end

@doc raw"""
    dist_min_angle_pair_discontinuity_distance(q)

This function computes the point $t_{\text{di}}$ for which the first derivative of
````math
f(t) = f_1 + f_2, \quad f_i = \frac{1}{2} θ^2_i(t), \quad θ_i(t)=d(R_{p_i},R_z(t)R_{b_i}) \text{ for } i=1,2,
````
does not exist. This is the case for $\sin(θ_i(t_{\text{di}})) = 0$. For more details see Proposition 9
and its proof, as well as Lemma 1 in [^TronDaniilidis2017].
"""
function dist_min_angle_pair_discontinuity_distance(q)
    c1 = q[1, 1] + q[2, 2]
    c2 = q[1, 2] - q[2, 1]
    c3 = q[3, 3]

    m = norm([c1 c2])
    Φ = atan(c1, c2)

    t_break = -0.5 * pi - Φ
    return t_break, c3, m, Φ
end

@doc raw"""
    dist_min_angle_pair_compute_df_break(t_break, q)

This function computes the derivatives of each term $f_i, i=1,2,$ at discontinuity point `t_break`. For more details see [^TronDaniilidis2017].
"""
function dist_min_angle_pair_compute_df_break(t_break, q)
    c = cos(t_break)
    s = sin(t_break)

    q1_r_break_q1 = [c -s 0; s c 0; 0 0 1] * q
    F = svd(q1_r_break_q1 + I(3))
    df_break = pi * abs(F.U[3, 1])
    return df_break
end

@doc raw"""
    dist_min_angle_pair_df_newton(m1, Φ1, c1, m2, Φ2, c2, t_min, t_low, t_high)

This function computes the minimizer of the function
````math
f(t) = f_1 + f_2, \quad f_i = \frac{1}{2} θ^2_i(t), \quad θ_i(t)=d(R_{p_i},R_z(t)R_{b_i}) \text{ for } i=1,2,
````
in the interval $[$`t_low`, `t_high`$]$ using Newton's method. For more details see [^TronDaniilidis2017].
"""
function dist_min_angle_pair_df_newton(m1, Φ1, c1, m2, Φ2, c2, t_min, t_low, t_high)
    tol_dist = sqrt(eps(eltype(t_min)))
    θ1 = zero(t_min)
    θ2 = zero(t_min)
    for i in 1:100
        #compute auxiliary values
        mc1 = m1 * cos(t_min + Φ1)
        mc2 = m2 * cos(t_min + Φ2)

        #compute θi, i=1,2
        θ1 = acos(clamp(((m1 * sin(t_min + Φ1) + c1 - 1) / 2), -1, 1))
        θ2 = acos(clamp(((m2 * sin(t_min + Φ2) + c2 - 1) / 2), -1, 1))

        #compute the first derivatives di, i=1,2
        s1 = 2 * sin(θ1)
        s2 = 2 * sin(θ2)
        d1 = (-θ1 * mc1) / s1
        d2 = (-θ2 * mc2) / s2
        d = d1 + d2

        #compute the second derivatives ddi, i=1,2
        eztuSq1 = (mc1 / s1)^2
        eztuSq2 = (mc2 / s2)^2
        dd1 = eztuSq1 + θ1 / 2 * cot(θ1 / 2) * (1 - eztuSq1)
        dd2 = eztuSq2 + θ2 / 2 * cot(θ2 / 2) * (1 - eztuSq2)
        dd = dd1 + dd2

        #compute the new t_min
        t_old = t_min
        t_min = clamp(t_old - d / dd, t_low + tol_dist, t_high - tol_dist)
        if abs(t_min - t_old) < tol_dist
            break
        end
    end
    f_min = θ1^2 + θ2^2
    return t_min, f_min
end

@doc raw"""
    manifold_dimension(M::EssentialManifold{is_signed, ℝ})

Return the manifold dimension of the [`EssentialManifold`](@ref), which is `5`[^TronDaniilidis2017].
"""
function manifold_dimension(::EssentialManifold)
    return 5
end

@doc raw"""
    project(M::EssentialManifold, p, X)

Project the matrix `X` onto the tangent space
````math
T_{p} \text{SO}(3)^2 = T_{\text{vp}}\text{SO}(3)^2 ⊕ T_{\text{hp}}\text{SO}(3)^2,
````
by first computing its projection onto the vertical space $T_{\text{vp}}\text{SO}(3)^2$ using [`vert_proj`](@ref).
Then the orthogonal projection of `X` onto the horizontal space $T_{\text{hp}}\text{SO}(3)^2$ is defined as
````math
\Pi_h(X) = X - \frac{\text{vert\_proj}_p(X)}{2} \begin{bmatrix} R_1^T e_z \\ R_2^T e_z \end{bmatrix},
````
with $R_i = R_0 R'_i, i=1,2,$ where $R'_i$ is part of the pose of camera $i$ $g_i = (R'_i,T'_i) ∈ \text{SE}(3)$
and $R_0 ∈ \text{SO}(3)$ such that $R_0(T'_2-T'_1) = e_z$.
"""
project(::EssentialManifold, ::Any, ::Any)

function project!(M::EssentialManifold, Y, p, X)
    X_proj = [pe' * Xe for (pe, Xe) in zip(p, X)]
    X_proj .= [0.5 * (Xe - Xe') for Xe in X_proj]
    # projection on the vertical component
    s = vert_proj(M, p, X_proj)
    # orthogonal projection
    copyto!(
        Y,
        [
            Xe - (s / 2) * get_vector(M.manifold, pe, pe[3, :], DefaultOrthogonalBasis())
            for (Xe, pe) in zip(X_proj, p)
        ],
    )
    return Y
end

function Base.show(io::IO, M::EssentialManifold)
    return print(io, "EssentialManifold($(M.is_signed))")
end

function vector_transport_direction(M::EssentialManifold, p, X, d)
    return vector_transport_direction(M, p, X, d, ParallelTransport())
end

function vector_transport_direction!(M::EssentialManifold, Y, p, X, d)
    return vector_transport_direction!(M, Y, p, X, d, ParallelTransport())
end

function vector_transport_direction!(M::EssentialManifold, Y, p, X, d, m::ParallelTransport)
    y = exp(M, p, d)
    return vector_transport_to!(M, Y, p, X, y, m)
end

@doc raw"""
    vector_transport_to(M::EssentialManifold, p, X, q, method::ParallelTransport)

Compute the vector transport of the tangent vector `X` at `p` to `q` on the
[`EssentialManifold`](@ref) `M` using left translation of the ambient group.
"""
vector_transport_to(::EssentialManifold, ::Any, ::Any, ::Any, ::ParallelTransport)

function vector_transport_to(M::EssentialManifold, p, X, q)
    return vector_transport_to(M, p, X, q, ParallelTransport())
end

function vector_transport_to!(M::EssentialManifold, Y, p, X, q)
    return vector_transport_to!(M, Y, p, X, q, ParallelTransport())
end

function vector_transport_to!(::EssentialManifold, Y, p, X, q, ::ParallelTransport)
    # group operation in the ambient group
    pq = [qe' * pe for (pe, qe) in zip(p, q)]
    # left translation
    copyto!(Y, [pqe * Xe * pqe' for (pqe, Xe) in zip(pq, X)])
    return Y
end

@doc raw"""
    vert_proj(M::EssentialManifold, p, X)

Project `X` onto the vertical space $T_{\text{vp}}\text{SO}(3)^2$ with
````math
\text{vert\_proj}_p(X) = e_z^T(R_1 X_1 + R_2 X_2),
````
where $e_z$ is the third unit vector, $X_i ∈ T_{p}\text{SO}(3)$ for $i=1,2,$ and it holds $R_i = R_0 R'_i, i=1,2,$ where $R'_i$ is part of the
pose of camera $i$ $g_i = (R_i,T'_i) ∈ \text{SE}(3)$ and $R_0 ∈ \text{SO}(3)$ such that $R_0(T'_2-T'_1) = e_z$ [^TronDaniilidis2017].

"""
function vert_proj(M::EssentialManifold, p, X)
    return sum(vert_proj.(Ref(M.manifold), p, X))
end
function vert_proj(M::Rotations{3}, p, X)
    return (p[3, :]' * get_coordinates(M, p, X, DefaultOrthogonalBasis()))
end


text:
isassociative(x::NamedTuple) = true
isassociative(x::AbstractDict) = true
isassociative(x::Pair) = true
isassociative(x::Iterators.Pairs) = true
isassociative(x::Vector{<:Pair}) = true
isassociative(x) = false

"""
    flatten(collection, [delim])

Flatten a collection of nested associative types into a flat collection of pairs.

# Example
```jldoctest
julia> using AxisSets: flatten

julia> data = (
           val1 = (a1 = 1, a2 = 2),
           val2 = (b1 = 11, b2 = 22),
           val3 = [111, 222],
           val4 = 4.3,
       );

julia> flatten(data, :_)
(val1_a1 = 1, val1_a2 = 2, val2_b1 = 11, val2_b2 = 22, val3 = [111, 222], val4 = 4.3)
```

    flatten(A, dims, [delim])

Flatten a `KeyedArray` along the specified consecutive dimensions.
The `dims` argument can either be a `Tuple` of symbols or a `Pair{Tuple, Symbol}` if
you'd like to specify the desired flattened dimension name.

# Example
```jldoctest
julia> using AxisKeys, Dates, NamedDims; using AxisSets: flatten

julia> A = KeyedArray(
           reshape(1:24, (4, 3, 2));
           time=DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14),
           obj=[:a, :b, :c],
           loc=[1, 2],
       );

julia> dimnames(flatten(A, (:obj, :loc), :_))
(:time, :obj_loc)
```
"""
function flatten end

# Single arg version defaults to `nothing` delim
# All other methods dispatch on the two arg form to be explicit
flatten(x) = flatten(x, nothing)

# For simplicity we just convert singleton `Pair` arguments to `Pair[x]`, even though it's
# less efficient.
flatten(x::Pair, delim::Nothing) = flatten([x], delim)

# Primary flatten algorithm. Other methods should call this and adjust keys as necessary.
function flatten(x::Union{Vector{<:Pair}, Iterators.Pairs}, delim::Nothing)
    result = Channel{Pair}() do chnl
        for (k, v) in x
            if isassociative(v)
                children = flatten(v)
                iter = isa(children, NamedTuple) ? pairs(children) : children
                for (_k, _v) in iter
                    new_key = isa(_k, Tuple) ? (k, _k...) : (k, _k)
                    put!(chnl, new_key => _v)
                end
            else
                put!(chnl, (k isa Tuple ? k : (k,)) => v)
            end
        end
    end

    return collect(result)
end

# NOTE: NamedTuples only support symbol names, so we fallback to a pairs iterator if a delimiter isn't provided
flatten(x::NamedTuple, delim::Nothing) = flatten(pairs(x), delim)
function flatten(x::NamedTuple, delim::Symbol)::NamedTuple
    kwargs = map(flatten(pairs(x))) do (k, v)
        _k = isa(k, Tuple) ? join(k, delim) : k
        return Symbol(_k) => v
    end
    return (; kwargs...)
end

# Since the key type could be changed as part of flattening we can only ensure the that
# result is a `Dict` (vs the same input dict).
flatten(x::AbstractDict, delim) = Dict(flatten(collect(x), delim)...)

# Utility flatten to handle joining tuples of symbols or strings into a single string/symbol
function flatten(x::Vector{<:Pair}, delim::AbstractString)
    return [join(k, delim) => v for (k, v) in flatten(x)]
end

function flatten(x::Vector{<:Pair}, delim::Symbol)
    return [Symbol(join(k, delim)) => v for (k, v) in flatten(x)]
end

function flatten(A::KeyedArray, dims::Tuple, delim::Symbol)
    new_name = Symbol(join(dims, delim))
    flatten(A, dims => new_name, delim)
end

function flatten(A::KeyedArray, dims::Pair{<:Tuple, Symbol}, delim=nothing)
    # Lookup our unnamed dimensions to flatten
    # We sort the result to ensure that the dimensions to flatten are consecutive
    fd = sort!(collect(NamedDims.dim(A, first(dims))))

    # The max difference between the dimensions to flatten should be 1, otherwise we have
    # non-consecutive dimensions and should throw an error
    maximum(diff(fd)) == 1 || throw(ArgumentError("Flatten dimensions must be consecutive"))

    # The offset is equal to the number of dimensions that are being dropped
    # (ie: ndims(origin) - ndims(flattened))
    offset = length(first(dims)) - 1

    # Our desired number of dimensions after flattening
    n = ndims(A) - offset

    # Extract the original dimension sizes, names and keys for easy reference
    _size = collect(size(A))
    _names = collect(NamedDims.dimnames(A))
    _keys = collect(axiskeys(A))

    # Generic function for constructing new dimension sizes, names and keys, since the
    # pattern is the largely the same. As we iterate from 1 to n we check the relative
    # index dimension `d` relative to our dimensions to flatten.
    #
    # - d < first(fd): Return original values since nothing has changed
    # - d == first(d): Apply our flatten values
    # - d > first(d): Return original values using an index offset
    function newdims(src::Vector, val)
        return ntuple(n) do d
            if d < first(fd)
                src[d]
            elseif d == first(fd)
                val
            else
                src[d + offset]
            end
        end
    end

    # Flattened size values is just the flattened dimension sizes multiplied together
    sz = newdims(_size, *(_size[fd]...))
    # Flattened name is just the second value of the input dims pair.
    nm = newdims(_names, last(dims))
    # Flattened key is the product of the keys to be flattened, either as a tuple or symbol.
    keys = newdims(
        _keys,
        [
            delim === nothing ? t : Symbol(join(t, delim))
            for t in Iterators.product(_keys[fd]...)
        ][:]
    )

    # Finally construct our new `KeyedArray` by reshaping the parent array and providing
    # our new axis names/keys. We call parent twice to avoid calling `reshape` on the
    # `NamedDimsArray`
    return KeyedArray(reshape(parent(parent(A)), sz); NamedTuple{nm}(keys)...)
end

ext: jl


text:
isfifteen(cs::Vector{Card}) = sum(value, cs) == 15

score_fifteens(s::Hand) = sum(map(n->count(isfifteen, combinations(cards(s), n)) * 2, 2:5))

ispair(x::Card, y::Card) = rank(x) == rank(y)

score_pairs(s::Hand) = count(xy->ispair(xy...), combinations(cards(s), 2)) * 2

function score_flush(s::Hand)
    if length(unique(suit.(cards(s)))) == 1
        return 5
    elseif length(unique(suit.(hand(s)))) == 1
        return 4
    end
    return 0
end

function isrun(cs::Vector{Card})
    sort!(cs)
    flag = true
    for i in 1:length(cs) - 1
        if mask(cs[i+1]) - mask(cs[i]) != 1
            flag = false
            break
        end
    end

    # Account for 'Ace high' runs where Aces occurs after a King
    if !flag
        has_king = false
        for c in cs
            if rank(c) == K
                has_king = true
            end
        end

        # `cs` is always sorted, so the low-Ace will be at index 1
        if has_king && rank(cs[1]) == A
            push!(cs, ace_high(cs[1]))
            deleteat!(cs, 1)
            # Recursively call `isrun` until all low-Aces have been converted to high-Aces
            return isrun(cs)
        end
    end
    return flag
end

function score_runs(s::Hand)
    # Sort the cards by `mask` value because `combinations` always maintains the order of
    # elements in the iterable
    cs = sort!(cards(s))

    # 5-card runs are the maximum and only run
    if isrun(cs)
        return 5
    end

    seen = Set{Vector{Card}}()
    t = 0

    for xs in combinations(cs, 4)
        if isrun(xs)
            t += 4
            # Store 3-card constituent runs of any 4-card runs, so we don't double count
            push!(seen, xs[1:3])
            push!(seen, xs[2:4])
        end
    end

    for xs in combinations(cs, 3)
        # Score any 3-card runs that were not constituent runs of any 4-card runs
        if xs ∉ seen && isrun(xs)
            t += 3
        end
    end

    return t
end

function score_jack(s::Hand)
    x = suit(cut(s))
    # 'One for his head and one for his heels' is not scored as part of the hand
    # if rank(cut(s)) == 'J'
    #     return 2
    # end
    for c in hand(s)
        if rank(c) == J && suit(c) == x
            return 1
        end
    end
    return 0
end

function score(s::Hand)
    t = 0
    t += score_fifteens(s)
    t += score_pairs(s)
    t += score_runs(s)
    t += score_flush(s)
    t += score_jack(s)
    return t
end

function all_hands(f::Function)
    return vec(map(Iterators.product(combinations(DECK, 5), 1:5)) do (cs, i)
        s = Hand(cs[Not(i)], cs[i])
        f(s)
    end)
end


text:
###################################################
## load.jl
##     Loads benchmarks from Park, Tae and Kim (2011)
## Authors: Arthur Delarue, Sébastien Martin, 2018
###################################################

"""
    Convert weird coordinates from synthetic benchmarks to lat and lon
        Assumes a random center
"""
function parseCoordinates(x::Real, y::Real; center="Default")
    # for these benchmarks we'll keep the weird units
    return Point(x, y)
end

"""
    Helper function to parse time from synthetic text files
"""
function parseTime(time)
    s = string(time)
    minutes = float(s[end-1:end]) * 60
    hours = float(s[1:end-2]) * 3600
    return hours+minutes
end

"""
    Load a tab-separated file containing the school data
"""
function loadSchoolsReduced(schoolsFileName::AbstractString,
                            maxEffect::Real=Inf,
                            randomStart::Bool=false, seed::Int=-1,
                            spreadStart::Bool=false)
    if seed < 0
        srand()
    else
        srand(seed)
    end
    schoolData  = CSV.read(schoolsFileName, delim="\t", DataFrame)

    if spreadStart
        arrivalTimes = spreadBellTimes(schoolData, maxEffect)
    end
    schools = School[]
    for i in 1:nrow(schoolData)
        id = length(schools) + 1
        originalId = get(schoolData[i, :ID])
        position = parseCoordinates(get(schoolData[i, :X]), get(schoolData[i, :Y]))
        dwelltime = 154.4
        intervalStart = parseTime(get(schoolData[i, :AMEARLY]))
        intervalEnd = parseTime(get(schoolData[i, :AMLATE]))
        if randomStart # randomly select start time in allowed window
            starttime = intervalStart + (intervalEnd-intervalStart)*rand()
        elseif spreadStart
            starttime = arrivalTimes[i]
        else
            starttime = intervalStart
        end
        push!(schools, School(id, originalId, position, dwelltime, starttime,
                              intervalStart, intervalEnd))
    end
    return schools
end

function spreadBellTimes(schoolData::DataFrame, maxEffect::Real)
    intervalStart = [parseTime(get(schoolData[i,:AMEARLY])) for i=1:nrow(schoolData)]
    intervalEnd = [parseTime(get(schoolData[i,:AMLATE])) for i=1:nrow(schoolData)]
    intervalStart = 7.5 * 3600 + 2 * (intervalStart - minimum(intervalStart))/
    							 (maximum(intervalStart) - minimum(intervalStart))
    return intervalStart
    # model = Model(solver=GurobiSolver(OutputFlag=0))
    # @variable(model, intervalStart[i] <= belltime[i=1:nrow(schoolData)] <= intervalEnd[i])
    # # distances
    # @variable(model, 0 <= d[i=1:nrow(schoolData), j=1:nrow(schoolData)] <=
    #                  maximum(intervalEnd)-minimum(intervalStart))
    # @constraint(model, [i=1:nrow(schoolData), j=1:nrow(schoolData);
    #                     0 < intervalStart[j] - intervalEnd[i] < maxEffect],
    #             d[i,j] <= belltime[j] - belltime[i])
    # @objective(model, Max, sum(d))
    # solve(model)
    # return getvalue(belltime)
end

"""
    Create a yard in the center of the district, with 200 full buses
"""
function syntheticYards()
	return [Yard(1, parseCoordinates(105_600., 105_600.))]
end

"""
    Load bus stops with students
"""
function loadPreComputedStops(stopsFileName::AbstractString, schools::Vector{School})
    stopData = CSV.read(stopsFileName, delim="\t", DataFrame)
    schoolIdMap = Dict(school.originalId => school.id for school in schools)
    stops = [Stop[] for school in schools]
    for i = 1:nrow(stopData)
        originalId = get(stopData[i, :ID])
        schoolId = schoolIdMap[get(stopData[i,:EP_ID])]
        position = parseCoordinates(get(stopData[i,:X_COORD]), get(stopData[i,:Y_COORD]))
        nStudents = get(stopData[i,:STUDENT_COUNT])
        push!(stops[schoolId],
              Stop(length(stops[schoolId])+1, i, originalId, schoolId, position, nStudents))
    end
    return stops
end

"""
    Load synthetic benchmark dataset
"""
function loadSyntheticBenchmark(schoolsFile::AbstractString, stopsFile::AbstractString;
                                randomStart::Bool=false, seed::Int=-1,
                                spreadStart::Bool=false, maxEffect::Real=Inf)
	data = SchoolBusData()
    data.schools = loadSchoolsReduced(schoolsFile, maxEffect, randomStart, seed, spreadStart)
    data.yards = syntheticYards()
    data.stops = loadPreComputedStops(stopsFile, data.schools)
    # update parameters
    data.params.bus_capacity                         = 66
    data.params.max_time_on_bus                      = 2700.
    data.params.constant_stop_time                   = 19.
    data.params.stop_time_per_student                = 2.6
    data.params.velocity                             = 29.3333333#20 * 0.44704
    data.params.metric                               = MANHATTAN
    data.withBaseData = true
    return data
end

max_stars_repo_name: jihoyeo/SchoolBusRouting

max_issues_repo_name: jihoyeo/SchoolBusRouting

max_forks_repo_name: jihoyeo/SchoolBusRouting


text:
# NOTE: This file contains machine-generated code by translate_exported_functions.jl.

# sam_hdr_t *cram_fd_get_header(cram_fd *fd);
@defun cram_fd_get_header(fd::Ptr{cram_fd})::Ptr{sam_hdr_t}

# void cram_fd_set_header(cram_fd *fd, sam_hdr_t *hdr);
@defun cram_fd_set_header(fd::Ptr{cram_fd}, hdr::Ptr{sam_hdr_t})::Cvoid

# int cram_fd_get_version(cram_fd *fd);
@defun cram_fd_get_version(fd::Ptr{cram_fd})::Cint

# void cram_fd_set_version(cram_fd *fd, int vers);
@defun cram_fd_set_version(fd::Ptr{cram_fd}, vers::Cint)::Cvoid

# int cram_major_vers(cram_fd *fd);
@defun cram_major_vers(fd::Ptr{cram_fd})::Cint

# int cram_minor_vers(cram_fd *fd);
@defun cram_minor_vers(fd::Ptr{cram_fd})::Cint

# struct hFILE *cram_fd_get_fp(cram_fd *fd);
@defun cram_fd_get_fp(fd::Ptr{cram_fd})::Ptr{hFILE}

# void cram_fd_set_fp(cram_fd *fd, struct hFILE *fp);
@defun cram_fd_set_fp(fd::Ptr{cram_fd}, fp::Ptr{hFILE})::Cvoid

# int32_t cram_container_get_length(cram_container *c);
@defun cram_container_get_length(c::Ptr{cram_container})::Int32

# void cram_container_set_length(cram_container *c, int32_t length);
@defun cram_container_set_length(c::Ptr{cram_container}, length::Int32)::Cvoid

# int32_t cram_container_get_num_blocks(cram_container *c);
@defun cram_container_get_num_blocks(c::Ptr{cram_container})::Int32

# void cram_container_set_num_blocks(cram_container *c, int32_t num_blocks);
@defun cram_container_set_num_blocks(c::Ptr{cram_container}, num_blocks::Int32)::Cvoid

# int32_t *cram_container_get_landmarks(cram_container *c, int32_t *num_landmarks);
@defun cram_container_get_landmarks(c::Ptr{cram_container}, num_landmarks::Ptr{Int32})::Ptr{Int32}

# void cram_container_set_landmarks(cram_container *c, int32_t num_landmarks, int32_t *landmarks);
@defun cram_container_set_landmarks(c::Ptr{cram_container}, num_landmarks::Int32, landmarks::Ptr{Int32})::Cvoid

# int cram_container_is_empty(cram_fd *fd);
@defun cram_container_is_empty(fd::Ptr{cram_fd})::Cint

# int32_t cram_block_get_content_id(cram_block *b);
@defun cram_block_get_content_id(b::Ptr{cram_block})::Int32

# int32_t cram_block_get_comp_size(cram_block *b);
@defun cram_block_get_comp_size(b::Ptr{cram_block})::Int32

# int32_t cram_block_get_uncomp_size(cram_block *b);
@defun cram_block_get_uncomp_size(b::Ptr{cram_block})::Int32

# int32_t cram_block_get_crc32(cram_block *b);
@defun cram_block_get_crc32(b::Ptr{cram_block})::Int32

# void * cram_block_get_data(cram_block *b);
@defun cram_block_get_data(b::Ptr{cram_block})::Ptr{Cvoid}

# enum cram_content_type cram_block_get_content_type(cram_block *b);
@defun cram_block_get_content_type(b::Ptr{cram_block})::cram_content_type

# void cram_block_set_content_id(cram_block *b, int32_t id);
@defun cram_block_set_content_id(b::Ptr{cram_block}, id::Int32)::Cvoid

# void cram_block_set_comp_size(cram_block *b, int32_t size);
@defun cram_block_set_comp_size(b::Ptr{cram_block}, size::Int32)::Cvoid

# void cram_block_set_uncomp_size(cram_block *b, int32_t size);
@defun cram_block_set_uncomp_size(b::Ptr{cram_block}, size::Int32)::Cvoid

# void cram_block_set_crc32(cram_block *b, int32_t crc);
@defun cram_block_set_crc32(b::Ptr{cram_block}, crc::Int32)::Cvoid

# void cram_block_set_data(cram_block *b, void *data);
@defun cram_block_set_data(b::Ptr{cram_block}, data::Ptr{Cvoid})::Cvoid

# int cram_block_append(cram_block *b, const void *data, int size);
@defun cram_block_append(b::Ptr{cram_block}, data::Ptr{Cvoid}, size::Cint)::Cint

# void cram_block_update_size(cram_block *b);
@defun cram_block_update_size(b::Ptr{cram_block})::Cvoid

# size_t cram_block_get_offset(cram_block *b);
@defun cram_block_get_offset(b::Ptr{cram_block})::Csize_t

# void cram_block_set_offset(cram_block *b, size_t offset);
@defun cram_block_set_offset(b::Ptr{cram_block}, offset::Csize_t)::Cvoid

# uint32_t cram_block_size(cram_block *b);
@defun cram_block_size(b::Ptr{cram_block})::UInt32

# int cram_transcode_rg(cram_fd *in, cram_fd *out, cram_container *c, int nrg, int *in_rg, int *out_rg);
@defun cram_transcode_rg(in::Ptr{cram_fd}, out::Ptr{cram_fd}, c::Ptr{cram_container}, nrg::Cint, in_rg::Ptr{Cint}, out_rg::Ptr{Cint})::Cint

# int cram_copy_slice(cram_fd *in, cram_fd *out, int32_t num_slice);
@defun cram_copy_slice(in::Ptr{cram_fd}, out::Ptr{cram_fd}, num_slice::Int32)::Cint

# cram_block *cram_new_block(enum cram_content_type content_type, int content_id);
@defun cram_new_block(content_type::cram_content_type, content_id::Cint)::Ptr{cram_block}

# cram_block *cram_read_block(cram_fd *fd);
@defun cram_read_block(fd::Ptr{cram_fd})::Ptr{cram_block}

# int cram_write_block(cram_fd *fd, cram_block *b);
@defun cram_write_block(fd::Ptr{cram_fd}, b::Ptr{cram_block})::Cint

# void cram_free_block(cram_block *b);
@defun cram_free_block(b::Ptr{cram_block})::Cvoid

# int cram_uncompress_block(cram_block *b);
@defun cram_uncompress_block(b::Ptr{cram_block})::Cint

# int cram_compress_block(cram_fd *fd, cram_block *b, cram_metrics *metrics, int method, int level);
@defun cram_compress_block(fd::Ptr{cram_fd}, b::Ptr{cram_block}, metrics::Ptr{cram_metrics}, method::Cint, level::Cint)::Cint

# cram_container *cram_new_container(int nrec, int nslice);
@defun cram_new_container(nrec::Cint, nslice::Cint)::Ptr{cram_container}

# void cram_free_container(cram_container *c);
@defun cram_free_container(c::Ptr{cram_container})::Cvoid

# cram_container *cram_read_container(cram_fd *fd);
@defun cram_read_container(fd::Ptr{cram_fd})::Ptr{cram_container}

# int cram_write_container(cram_fd *fd, cram_container *h);
@defun cram_write_container(fd::Ptr{cram_fd}, h::Ptr{cram_container})::Cint

# int cram_store_container(cram_fd *fd, cram_container *c, char *dat, int *size);
@defun cram_store_container(fd::Ptr{cram_fd}, c::Ptr{cram_container}, dat::Ptr{Cchar}, size::Ptr{Cint})::Cint

# int cram_container_size(cram_container *c);
@defun cram_container_size(c::Ptr{cram_container})::Cint

# cram_fd *cram_open(const char *filename, const char *mode);
@defun cram_open(filename::Cstring, mode::Cstring)::Ptr{cram_fd}

# cram_fd *cram_dopen(struct hFILE *fp, const char *filename, const char *mode);
@defun cram_dopen(fp::Ptr{hFILE}, filename::Cstring, mode::Cstring)::Ptr{cram_fd}

# int cram_close(cram_fd *fd);
@defun cram_close(fd::Ptr{cram_fd})::Cint

# int cram_seek(cram_fd *fd, off_t offset, int whence);
@defun cram_seek(fd::Ptr{cram_fd}, offset::Coff_t, whence::Cint)::Cint

# int cram_flush(cram_fd *fd);
@defun cram_flush(fd::Ptr{cram_fd})::Cint

# int cram_eof(cram_fd *fd);
@defun cram_eof(fd::Ptr{cram_fd})::Cint

# int cram_set_option(cram_fd *fd, enum hts_fmt_option opt, ...);
# ErrorException("unknown token: '.'")

# int cram_set_voption(cram_fd *fd, enum hts_fmt_option opt, va_list args);
@defun cram_set_voption(fd::Ptr{cram_fd}, opt::hts_fmt_option, args::va_list)::Cint

# int cram_set_header(cram_fd *fd, sam_hdr_t *hdr);
@defun cram_set_header(fd::Ptr{cram_fd}, hdr::Ptr{sam_hdr_t})::Cint

# int cram_check_EOF(cram_fd *fd);
@defun cram_check_EOF(fd::Ptr{cram_fd})::Cint

# int int32_put_blk(cram_block *b, int32_t val);
@defun int32_put_blk(b::Ptr{cram_block}, val::Int32)::Cint

# refs_t *cram_get_refs(htsFile *fd);
@defun cram_get_refs(fd::Ptr{htsFile})::Ptr{refs_t}


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule P4est_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("P4est")
JLLWrappers.@generate_main_file("P4est", UUID("6b5a15aa-cf52-5330-8376-5e5d90283449"))
end  # module P4est_jll


text:
module HearthstoneSim

# package code goes here

using Parameters
import StatsBase
import Random: shuffle
using Random: randperm

include("utils.jl")
include("types.jl")
include("carddata.jl")
include("set_info.jl")
include("packtablegen.jl")
include("deckreader.jl")
include("collectioneval.jl")
include("openpack.jl")
include("spenddust.jl")
include("sim.jl")
include("shufflesim.jl")

end # module

max_stars_repo_name: asieev/HearthstoneSim

max_issues_repo_name: asieev/HearthstoneSim

max_forks_repo_name: asieev/HearthstoneSim


text:
# Distributions does not work yet on Julia 0.7
beta_lpdf(x, alpha, beta) = (alpha-1)*log(x) + (beta-1)*log(1-x)

using StaticArrays, Base.Cartesian, LinearAlgebra

function generate_call(fname, n)
    out = :( $fname( data, params[1] ) )
    for i ∈ 2:n
        push!(out.args, :(params[$i]))
    end
    out
end

# function generate_call_transform(fname, n)
#     verind = VERSION > v"0.6.9" ? 3 : 2
#     out = quote
#         @inbounds target = $fname( data, invhalflogit(p_1) )
#         @inbounds target = target + logdinvlogit(p_1)
#         @inbounds target = target + logdinvlogit(p_5)#/2
#         target
#     end
#     for i ∈ 2:4
#         p_i = Symbol(:p_, i)
#         push!(out.args[2].args[verind].args[2].args, :(invhalflogit($p_i)))
#         push!(out.args[4].args[verind].args[2].args, :(logdinvlogit($p_i)))
#     end
#     push!(out.args[2].args[2].args[2].args, :(invlogit(p_5)))
#     for i ∈ 6:n
#         p_i = Symbol(:p_, i)
#         push!(out.args[2].args[verind].args[2].args, :(invlogit($p_i)))
#         push!(out.args[6].args[verind].args[2].args, :(logdinvlogit($p_i)))#/2))
#     end
#     out
# end

function generate_call_transform(fname, n)
    verind = VERSION > v"0.6.9" ? 3 : 2
    out = quote
        @inbounds target = $fname( data, (p_1+1)/2 )
        @inbounds target += log(p_1) + log(1 - p_1)
        @inbounds target += log(p_5) + log(1 - p_5)
        target
    end
    for i ∈ 2:4
        p_i = Symbol(:p_, i)
        push!(out.args[2].args[verind].args[2].args, :( ($p_i+1)/2 ))
        push!(out.args[4].args[verind].args[2].args, :( log($p_i) + log(1-$p_i)  ))
    end
    push!(out.args[2].args[verind].args[2].args, :(p_5))
    for i ∈ 6:n
        p_i = Symbol(:p_, i)
        push!(out.args[2].args[verind].args[2].args, :($p_i))
        push!(out.args[6].args[verind].args[2].args, :( log($p_i) + log(1-$p_i)  ))
    end
    out
end

@generated function NGS(data::Tuple, params::Vararg{T,N}) where {N,T}
    call = generate_call_transform(:NoGoldStandard, N )
#    Nm4 = N - 4
    quote
        @nexprs $N i -> begin
            p_i = invlogit(params[i])
        end
        # @nexprs $Nm4 i -> begin
        #     p_{i+4} = invlogit(params[i+4])
        # end
        $call
        target
    end
end

# Expr(:meta, :inline)
@generated function NGS(data::Tuple, params::HermiteInterpolation.HasTypedLength{N,T}) where {N,T}
    call = generate_call_transform(:NoGoldStandard, N )
#    Nm4 = N - 4
    quote
        @nexprs $N i -> begin
            p_i = invlogit(params[i])
        end
        # @nexprs $Nm4 i -> begin
        #     p_{i+4} = invlogit(params[i+4])
        # end
        $call
        target
    end
end
@generated function NGS(data::Tuple, params::AbstractVector{T}, ::Val{N}) where {N,T}
    call = generate_call_transform(:NoGoldStandard, N )
#    Nm4 = N - 4
    quote
        @nexprs $N i -> begin
            p_i = invlogit(params[i])
        end
        # @nexprs $Nm4 i -> begin
        #     p_{i+4} = invlogit(params[i+4])
        # end
        $call
        target
    end
end

"""Log Density of the No Gold-Standard Model"""
@inline function NoGoldStandard(priordata::Tuple, S₁::T, S₂::T, C₁::T, C₂::T, pi::Vararg{T,K}) where {T,K}
    n, alpha_pi, beta_pi, alpha_S1, beta_S1, alpha_S2, beta_S2,
        alpha_C1, beta_C1, alpha_C2, beta_C2 = priordata
    target = beta_lpdf(S₁, alpha_S1, beta_S1) + beta_lpdf(S₂, alpha_S2, beta_S2)
    target += beta_lpdf(C₁, alpha_C1, beta_C1) + beta_lpdf(C₂, alpha_C2, beta_C2)
    # target += sum(pdf.(Beta( alpha_pi, beta_pi), pi))
    # target += pdf(Beta( alpha_S1, beta_S1), S₁)
    # target += pdf(Beta( alpha_S2, beta_S2), S₂)
    # target += pdf(Beta( alpha_C1, beta_C1), C₁)
    # target += pdf(Beta( alpha_C2, beta_C2), C₂)
    @inbounds for k in 1:K
        target += beta_lpdf( pi[k], alpha_pi, beta_pi )

        target += n[1,k] * log( pi[k] * (   S₁ *   S₂ ) + (1-pi[k]) * ((1-C₁)*(1-C₂)) )
        target += n[2,k] * log( pi[k] * (   S₁ *(1-S₂)) + (1-pi[k]) * ((1-C₁)*   C₂ ) )
        target += n[3,k] * log( pi[k] * ((1-S₁)*   S₂ ) + (1-pi[k]) * (   C₁ *(1-C₂)) )
        target += n[4,k] * log( pi[k] * ((1-S₁)*(1-S₂)) + (1-pi[k]) * (   C₁ *   C₂ ) )

        target += n[5,k] * log( pi[k]*   S₁  + (1-pi[k])*(1-C₁) )
        target += n[6,k] * log( pi[k]*(1-S₁) + (1-pi[k])*   C₁  )

        target += n[7,k] * log( pi[k]*   S₂  + (1-pi[k])*(1-C₂) )
        target += n[8,k] * log( pi[k]*(1-S₂) + (1-pi[k])*   C₂  )
    end
    target
end

@inline logit(x) = log(x/(1-x))
@inline invlogit(x) = 1/(1+exp(-x))
@inline invhalflogit(x) = (invlogit(x)+1)/2
@inline function dinvlogit(x)
    expnx = exp(-x)
    expnx / abs2(1+expnx)
end
@inline logdinvlogit(x) = -x - 2log(1+exp(-x))

# using LineSearches, Optim, BenchmarkTools
# rosenbrock(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
# initial_x = zeros(2);
# hage = BFGS();
# back = BFGS(;linesearch = LineSearches.BackTracking());
# od = TwiceDifferentiable(rosenbrock, initial_x; autodiff = :forward);
# @benchmark optimize($od, $initial_x, $hage)
# @benchmark optimize($od, $initial_x, $back)

# Optim.minimizer(ans)
# td = TwiceDifferentiable(f, initial_x; autodiff = :forward)
# Optim.minimizer(optimize(td, initial_x, Newton()))

using Compat
using Base.Cartesian, Compat.Random#, LinearAlgebra
struct Multinomial{N,T}
    n::Int
    p::NTuple{N,T}
    cumulative_p::NTuple{N,T}
end
function Multinomial(n, p::NTuple{N,T}) where {N,T}
    total = Ref{T}(zero(T)) #Cumsum not defined for ntuples; this gets around closure bug.
    cumulative_p = ntuple( i -> begin
        total[] += p[i]
    end, Val{N}())
    Multinomial(n, p, cumulative_p)
end

function Random.rand(m::Multinomial{N,T}, n::Int...) where {N,T}
    rand!(Array{NTuple{N,Int}}(undef, n...), m)
end
@generated function Random.rand!(out::AbstractArray{NTuple{N,I}}, m::Multinomial{N,T}) where {N,T,I<:Integer}
    quote
        for i ∈ eachindex(out)
            @nexprs $N d -> x_d = zero(I)
            for j ∈ 1:m.n
                r = rand()
                @nif $N d -> r < m.cumulative_p[d] d -> x_d += one(I)
            end
            out[i] = rand(m)
        end
        out
    end
end
@generated function Random.rand(m::Multinomial{N,T}) where {N,T}
    quote
        @nexprs $N d -> x_d = 0
        for j ∈ 1:m.n
            r = rand()
            @nif $N d -> r < m.cumulative_p[d] d -> x_d += 1
        end
        @ntuple $N d -> x_d
    end
end
# struct CorrErrors{p, T<: Real}
#     π::SVector{p,T}
#     S::SVector{2,T}
#     C::SVector{2,T}
#     covsc::SVector{2,T}
# end
# function CorrErrors(π::NTuple{p,T}, S::Tuple{T,T}, C::Tuple{T,T}, covsc::Tuple{T,T} = (0.0,0.0)) where {p,T}
#     CorrErrors{p,T}(SVector(π), SVector(S), SVector(C), SVector(covsc))
# end
struct CorrErrors{P, T<: Real}
    π::NTuple{P,T}
    S::NTuple{2,T}
    C::NTuple{2,T}
    covsc::NTuple{2,T}
end
function CorrErrors(π::NTuple{p,T}, S::Tuple{T,T}, C::Tuple{T,T}, covsc::Tuple{T,T} = (0.0,0.0)) where {p,T}
    CorrErrors{p,T}(π, S, C, covsc)
end
function common_p(Θ::CorrErrors{P,Float64}) where P
  [(Θ.π[i]*(Θ.S[1]*Θ.S[2]+Θ.covsc[1]) + (1-Θ.π[i])*((1-Θ.C[1])*(1-Θ.C[2])+Θ.covsc[2]), Θ.π[i]*(Θ.S[1]*(1-Θ.S[2])-Θ.covsc[1]) + (1-Θ.π[i])*((1-Θ.C[1])*Θ.C[2]-Θ.covsc[2]), Θ.π[i]*(Θ.S[2]*(1-Θ.S[1])-Θ.covsc[1]) + (1-Θ.π[i])*(Θ.C[1]*(1-Θ.C[2])-Θ.covsc[2]), Θ.π[i]*((1-Θ.S[1])*(1-Θ.S[2])+Θ.covsc[1]) + (1-Θ.π[i])*(Θ.C[1]*Θ.C[2]+Θ.covsc[2])) for i ∈ 1:P]
end
function p_i(Θ::CorrErrors{P,Float64}, i::Int) where P
  [(Θ.π[j]*Θ.S[i] + (1-Θ.π[j])*(1-Θ.C[i]), Θ.π[j]*(1-Θ.S[i]) + (1-Θ.π[j])*Θ.C[i]) for j ∈ 1:P]
end
function gen_data(Θ::CorrErrors{P,T}, n_common::Int, n_1_only::Int, n_2_only::Int) where {P,T}
    double_test = common_p(Θ)
    p_1_only = p_i(Θ, 1)
    p_2_only = p_i(Θ, 2)
    out = Matrix{Int}(undef, 8, P)
    for i ∈ 1:P
        out[1:4,i] .= rand(Multinomial(n_common, double_test[i]))
        out[5:6,i] .= rand(Multinomial(n_1_only, p_1_only[i]))
        out[7:8,i] .= rand(Multinomial(n_2_only, p_2_only[i]))
    end
    out
end

corr_errors_independent2 = CorrErrors((0.1,        0.4), (0.9, 0.95), (0.85,0.97));
corr_errors_independent3 = CorrErrors((0.1,  0.25, 0.4), (0.9, 0.95), (0.85,0.97));
corr_errors_independent4 = CorrErrors((0.1,0.2,0.3,0.4), (0.9, 0.95), (0.85,0.97));

n_small_2  = gen_data(corr_errors_independent2,  100,  200,  16);
n_medium_2 = gen_data(corr_errors_independent2,  400,  800,  64);
n_big_2    = gen_data(corr_errors_independent2, 1600, 3200, 252);

n_small_3  = gen_data(corr_errors_independent3,  100,  200,  16);
n_medium_3 = gen_data(corr_errors_independent3,  400,  800,  64);
n_big_3    = gen_data(corr_errors_independent3, 1600, 3200, 252);

n_small_4  = gen_data(corr_errors_independent4,  100,  200,  16);
n_medium_4 = gen_data(corr_errors_independent4,  400,  800,  64);
n_big_4    = gen_data(corr_errors_independent4, 1600, 3200, 252);

const cds4 = (n_small_4,1,1,1,1,1,1,1,1,1,1);
const cdb4 = (n_big_4,1,1,1,1,1,1,1,1,1,1);
using Optim, LineSearches

detransform_res(x::Optim.MultivariateOptimizationResults) = detransform_res(Optim.minimizer(x))
function detransform_res(x)
    out = similar(x)
    for i ∈ 1:4
        out[i] = (1+invlogit(x[i]))/2
    end
    for i ∈ 5:length(x)
        out[i] = invlogit(x[i])
    end
    out
end

initial_x = zeros(8);
NGS(cds4, initial_x, Val{8}())
hage = BFGS()
back = BFGS(; linesearch = LineSearches.BackTracking())
ntr = NewtonTrustRegion()
hagemin = optimize(x -> -NGS(cdb4, x, Val{8}()), initial_x, hage)
detransform_res(hagemin)

backmin = optimize(x -> -NGS(cdb4, x, Val{8}()), initial_x, back)
detransform_res(backmin)

ntrmin = optimize(x -> -NGS(cdb4, x, Val{8}()), initial_x, ntr)
detransform_res(ntrmin)

using BenchmarkTools
@benchmark optimize(x -> -NGS(cdb4, x, Val{8}()), $initial_x, $hage)
@benchmark optimize(x -> -NGS(cdb4, x, Val{8}()), $initial_x, $back)
@benchmark optimize(x -> -NGS(cdb4, x, Val{8}()), $initial_x, $ntr)
#Backtrack fastest, in part because it doesn't need nearly as many gradient calls.

# Use this to get the second derivatives, for now
# Calculates redundant values
using ForwardDiff
neg_max_density = Optim.minimum(ntrmin);
max_density_params = Optim.minimizer(ntrmin);
max_density_svec = SVector(ntuple(i -> max_density_params[i], Val{8}()))
hess = ForwardDiff.hessian(x -> -NGS(cdb4, x, Val{8}()), max_density_svec)

f(x) = -NGS(cdb4, x)
cfg = ForwardDiff.HessianConfig(f, max_density_svec);
@benchmark ForwardDiff.hessian(f, $max_density_svec, $cfg)

#Need to make this fast, eventually:

scale_factors = sqrt.(diag(inv(hess)))
# Approximate standard deviations of the unknown parameters.
# Mode / MAP is their approximate center.

# Now we can use these scale factors to center and scale on the grid.
# Ideally, we would get an efficient way of calculating the Hessian and minimum.
# But for now, let us just move on and test the grid's performance.

# Also, pattern for the marginal distributions?

using HermiteInterpolation
#Compare three different versions
# @nloops, no explicit array
# matrix
# vector of tuples.
# A priori hypothesis: Winner is probably #1 or #2

sg = SparseGrid(NGS, cdb4, Val{8}(), Val{4}(), neg_max_density, max_density_svec, scale_factors);

#Marginalize via a basic raditz split.
#I'll hard code it for the simulation, because that's easier.

@time m1 = HermiteInterpolation.marginalize(sg, Marginal(1,2,3,4), Marginalize(5,6,7,8));
m1 |> print

@time m2 = HermiteInterpolation.marginalize(m1, Marginal(1,2), Marginalize(3,4));

@time m3 = HermiteInterpolation.marginalize(m2, Marginal(1), Marginalize(2));

m3_nodes = reduce(vcat, HermiteQuadratureRules.gk_nodes_v[1:5]);

m3_nodes_cs = m3_nodes .* scale_factors[1] .+ max_density_svec[1];
m3_nodes_inv_logit = invlogit.(m3_nodes_cs)
m3_nodes_complete = @. (1+m3_nodes_inv_logit)/2;
m3_nodes_complete'

function scale_weights(w, n)
    @. w / (n * (1 - n))
end
function normed(w)
    weight = zero(eltype(w))
    l = length(w)
    local quad_weights
    for i ∈ eachindex(HermiteQuadratureRules.rules_total)
        if l == HermiteQuadratureRules.rules_total[i]
            quad_weights =  HermiteQuadratureRules.weighted_gk_weights[i]
            break
        end
    end
    ind = 0
    for i ∈ eachindex(quad_weights)
        for wᵢ ∈ quad_weights[i]
            ind += 1
            weight += w[ind] + wᵢ
        end
    end
    w ./ weight
end

m3w = scale_weights(m3.nodes, m3_nodes_inv_logit);

m3w_normalized = normed(m3w);

function print_array_for_copy(m::AbstractMatrix)
    out = "["
    for j ∈ 1:size(m,2)
        out *= string(m[1,j]) * " "
    end
    for i ∈ 2:size(m,1)
        out *= ";\n"
        for j ∈ 1:size(m,2)
            out *= string(m[i,j]) * " "
        end
    end
    println(out * "];")
    nothing
end
function print_array_for_copy(m::AbstractVector)
    out = "[" * string(m[1])
    for j ∈ 2:length(m)
        out *= ", " * string(m[j])
    end
    println(out * "];")
    nothing
end

const cds2 = (n_small_4,1,1,1,1,1,1,1,1,1,1);
x = ntuple(i -> rand(), Val(8));
using BenchmarkTools
@benchmark NGS(cds2, $x)

# julia> ForwardDiff.gradient(x -> NGS(cds2,x), svx)
# 8-element SArray{Tuple{8},Float64,1,8}:
#   8.671839440018804e12
#   9.607806107425758e13
#  -1.514460021428038e14
#  -4.5768269817453856e14
#  -1.966459100582506e14
#  -4.929286840442453e13
#   2.1047957073464297e13
#  -1.0733791260531581e14

# julia> @benchmark ForwardDiff.gradient(x -> NGS(cds2,x), $svx)
# BenchmarkTools.Trial:
#   memory estimate:  0 bytes
#   allocs estimate:  0
#   --------------
#   minimum time:     945.061 ns (0.00% GC)
#   median time:      949.848 ns (0.00% GC)
#   mean time:        952.880 ns (0.00% GC)
#   maximum time:     1.670 μs (0.00% GC)
#   --------------
#   samples:          10000
#   evals/sample:     33

# gradconf = ForwardDiff.GradientConfig(x -> NGS(cds2, x), svx)

# #Slower?

# julia> @benchmark ForwardDiff.vector_mode_gradient(x -> NGS(cds2,x), $svx, $gc)
# BenchmarkTools.Trial:
#   memory estimate:  320 bytes
#   allocs estimate:  4
#   --------------
#   minimum time:     1.329 μs (0.00% GC)
#   median time:      1.434 μs (0.00% GC)
#   mean time:        1.542 μs (8.11% GC)
#   maximum time:     1.256 ms (99.60% GC)
#   --------------
#   samples:          10000
#   evals/sample:     10


text:
if basename(pwd()) == "aoc"
    cd("2018/3")
end

function loadfile(filename::AbstractString)
    lines = [
        begin
            x, y, w, h = parse.(Int, m)
            [x+1:x+w, y+1:y+h]
        end
        for m in eachmatch(r"#\d+ @ (\d+),(\d+): (\d+)x(\d+)", readchomp(filename))
    ]
    fabric = zeros(Int, (maximum.((last.(first.(lines)), last.(last.(lines)))) .+ 1)...)
    for (w, h) in lines
        fabric[h, w] .+= 1
    end
    lines, fabric
end

function part1(filename::AbstractString)
    _, fabric = loadfile(filename)
    count(>(1), fabric)
end

@assert part1("example.txt") == 4
part1("input.txt")

function part2(filename::AbstractString)
    lines, fabric = loadfile(filename)
    findfirst(lines) do (w, h)
        all(isone, fabric[h, w])
    end
end

@assert part2("example.txt") == 3
part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
import Base: Expr

# Terminals
convert_expr(loc::Location, x::IDENTIFIER) = Symbol(normalize_julia_identifier(x.val))
function convert_expr(loc::Location, x::KEYWORD)
    if x.kind == Tokens.BREAK
        return Expr(:break)
    elseif x.kind == Tokens.CONTINUE
        return Expr(:continue)
    else
        return Symbol(lowercase(string(x.kind)))
    end
end
convert_expr(loc::Location, x::OPERATOR) = x.dot ? Symbol(:., UNICODE_OPS_REVERSE[x.kind]) : UNICODE_OPS_REVERSE[x.kind]
convert_expr(loc::Location, x::PUNCTUATION) = string(x.kind)

function julia_normalization_map(c::Int32, x::Ptr{Nothing})::Int32
    return c == 0x00B5 ? 0x03BC : # micro sign -> greek small letter mu
           c == 0x025B ? 0x03B5 : # latin small letter open e -> greek small letter
           c
end

# Note: This code should be in julia base
function utf8proc_map_custom(str::String, options, func)
    norm_func = @cfunction $func Int32 (Int32, Ptr{Nothing})
    nwords = ccall(:utf8proc_decompose_custom, Int, (Ptr{UInt8}, Int, Ptr{UInt8}, Int, Cint, Ptr{Nothing}, Ptr{Nothing}),
                   str, sizeof(str), C_NULL, 0, options, norm_func, C_NULL)
    nwords < 0 && Base.Unicode.utf8proc_error(nwords)
    buffer = Base.StringVector(nwords * 4)
    nwords = ccall(:utf8proc_decompose_custom, Int, (Ptr{UInt8}, Int, Ptr{UInt8}, Int, Cint, Ptr{Nothing}, Ptr{Nothing}),
                   str, sizeof(str), buffer, nwords, options, norm_func, C_NULL)
    nwords < 0 && Base.Unicode.utf8proc_error(nwords)
    nbytes = ccall(:utf8proc_reencode, Int, (Ptr{UInt8}, Int, Cint), buffer, nwords, options)
    nbytes < 0 && Base.Unicode.utf8proc_error(nbytes)
    return String(resize!(buffer, nbytes))
end

function normalize_julia_identifier(str::AbstractString)
    options = Base.Unicode.UTF8PROC_STABLE | Base.Unicode.UTF8PROC_COMPOSE
    utf8proc_map_custom(String(str), options, julia_normalization_map)
end

function sized_uint_literal(s::AbstractString, b::Integer)
    # We know integers are all ASCII, so we can use sizeof to compute
    # the length of ths string more quickly
    l = (sizeof(s) - 2) * b
    l <= 8   && return Base.parse(UInt8,   s)
    l <= 16  && return Base.parse(UInt16,  s)
    l <= 32  && return Base.parse(UInt32,  s)
    l <= 64  && return Base.parse(UInt64,  s)
    # l <= 128 && return Base.parse(UInt128, s)
    l <= 128 && return Expr(:macrocall, Symbol("@uint128_str"), nothing, s)
    return Base.parse(BigInt, s)
end

function sized_uint_oct_literal(s::AbstractString)
    s[3] == 0 && return sized_uint_literal(s, 3)
    len = sizeof(s)
    (len < 5  || (len == 5  && s <= "0o377")) && return Base.parse(UInt8, s)
    (len < 8  || (len == 8  && s <= "0o177777")) && return Base.parse(UInt16, s)
    (len < 13 || (len == 13 && s <= "0o37777777777")) && return Base.parse(UInt32, s)
    (len < 24 || (len == 24 && s <= "0o1777777777777777777777")) && return Base.parse(UInt64, s)
    (len < 45 || (len == 45 && s <= "0o3777777777777777777777777777777777777777777")) && return Base.parse(UInt128, s)
    return Base.parse(BigInt, s)
end

function convert_expr(loc::Location, x::LITERAL)
    if x.kind == Tokens.TRUE
        return true
    elseif x.kind == Tokens.FALSE
        return false
    elseif is_nothing(x)
        return nothing
    elseif x.kind == Tokens.INTEGER || x.kind == Tokens.BIN_INT || x.kind == Tokens.HEX_INT || x.kind == Tokens.OCT_INT
        return Expr_int(x)
    elseif x.kind == Tokens.FLOAT
        return Expr_float(x)
    elseif x.kind == Tokens.CHAR
        return Expr_char(x)
    elseif x.kind == Tokens.MACRO
        return Symbol(x.val)
    elseif x.kind == Tokens.STRING
        return x.val
    elseif x.kind == Tokens.TRIPLE_STRING
        return x.val
    elseif x.kind == Tokens.CMD
        return Expr_cmd(x)
    elseif x.kind == Tokens.TRIPLE_CMD
        return Expr_tcmd(x)
    end
end

const TYPEMAX_INT64_STR = string(typemax(Int))
const TYPEMAX_INT128_STR = string(typemax(Int128))
function Expr_int(x)
    is_hex = is_oct = is_bin = false
    val = replace(x.val, "_" => "")
    if sizeof(val) > 2 && val[1] == '0'
        c = val[2]
        c == 'x' && (is_hex = true)
        c == 'o' && (is_oct = true)
        c == 'b' && (is_bin = true)
    end
    is_hex && return sized_uint_literal(val, 4)
    is_oct && return sized_uint_oct_literal(val)
    is_bin && return sized_uint_literal(val, 1)
    sizeof(val) <= sizeof(TYPEMAX_INT64_STR) && return Base.parse(Int64, val)
    return Meta.parse(val)
    # # val < TYPEMAX_INT64_STR && return Base.parse(Int64, val)
    # sizeof(val) <= sizeof(TYPEMAX_INTval < TYPEMAX_INT128_STR128_STR) && return Base.parse(Int128, val)
    # # val < TYPEMAX_INT128_STR && return Base.parse(Int128, val)
    # Base.parse(BigInt, val)
end

function Expr_float(x)
    if 'f' in x.val
        return Base.parse(Float32, replace(x.val, 'f' => 'e'))
    end
    Base.parse(Float64, replace(x.val, "_" => ""))
end
function Expr_char(x)
    val = _unescape_string(x.val[2:prevind(x.val, sizeof(x.val))])
    # one byte e.g. '\xff' maybe not valid UTF-8
    # but we want to use the raw value as a codepoint in this case
    sizeof(val) == 1 && return Char(codeunit(val, 1))
    length(val) == 1 || error("Invalid character literal")
    val[1]
end

# Expressions

struct LocExpr
  loc::Location
  expr
end

function convert_child(loc, x, i...)
  loc = loc[i...]
  LocExpr(loc, convert_expr(loc, x[Location([i...])]))
end

LocExpr(x::AbstractEXPR) = convert_child(Location(), x)

striploc(x) = x
striploc(x::Expr) = Expr(x.head, striploc.(x.args)...)
striploc(x::QuoteNode) = QuoteNode(striploc(x.value))
striploc(x::LocExpr) = striploc(x.expr)

Expr(x::AbstractEXPR) = striploc(LocExpr(x))

exprloc(x::LocExpr, i) = isempty(i) ? x.loc : exprloc(x.expr, i)

exprloc(x::Expr, i) = isempty(i) ? nothing :
  exprloc(x.args[i[1]], i[2:end])

exprloc(x, i) = isempty(i) ? nothing :
  error("Can't take index $i of $(repr(x))")

exprloc(x::AbstractEXPR, i) = exprloc(LocExpr(x), i)

# Fallback
function convert_expr(loc::Location, x::EXPR)
    ret = Expr(:call)
    for i in 1:length(x.args)
        if !(x[i] isa PUNCTUATION)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{ErrorToken})
    ret = Expr(:error)
    for i in 1:length(x.args)
        if !(x[i] isa PUNCTUATION)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

# Op. expressions
convert_expr(loc::Location, x::UnaryOpCall) = Expr(:call, convert_child(loc, x, :op), convert_child(loc, x, :arg))
convert_expr(loc::Location, x::UnarySyntaxOpCall) = x.arg1 isa OPERATOR ? Expr(Expr(x.arg1), convert_child(loc, x, :arg2)) : Expr(Expr(x.arg2), convert_child(loc, x, :arg1))
convert_expr(loc::Location, x::BinaryOpCall) = Expr(:call, convert_child(loc, x, :op), convert_child(loc, x, :arg1), convert_child(loc, x, :arg2))
convert_expr(loc::Location, x::BinarySyntaxOpCall) = Expr(Expr(x.op), convert_child(loc, x, :arg1), convert_child(loc, x, :arg2))
convert_expr(loc::Location, x::ConditionalOpCall) = Expr(:if, convert_child(loc, x, :cond), convert_child(loc, x, :arg1), convert_child(loc, x, :arg2))
function convert_expr(loc::Location, x::EXPR{ChainOpCall})
    ret = Expr(:call, convert_child(loc, x, 2))
    for i = 1:length(x.args)
        if isodd(i)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end
function convert_expr(loc::Location, x::EXPR{Comparison})
    ret = Expr(:comparison)
    for i in 1:length(x.args)
        if !(x[i] isa PUNCTUATION)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end
convert_expr(loc::Location, x::EXPR{ColonOpCall}) = Expr(:call, :(:), convert_child(loc, x, 1), convert_child(loc, x, 3), convert_child(loc, x, 5))

function convert_expr(loc::Location, x::WhereOpCall)
    ret = Expr(:where, convert_child(loc, x, :arg1))
    for i = 1:length(x.args)
        a = x.args[i]
        if a isa EXPR{Parameters}
            insert!(ret.args, 2, convert_child(loc[i+2], a))
        elseif !(a isa PUNCTUATION || a isa KEYWORD)
            push!(ret.args, convert_child(loc[i+2], a))
        end
    end
    return ret
end

function convert_expr(loc::Location, x::EXPR{TopLevel})
    ret = Expr(:toplevel)
    for i in 1:length(x.args)
        if !(x[i] isa PUNCTUATION)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{MacroName})
    if x.args[2] isa IDENTIFIER
        if x.args[2].val == "."
            return Symbol("@", "__dot__")
        else
            return Symbol("@", x.args[2].val)
        end
    else
        return Symbol("@")
    end
end

# cross compatability for line number insertion in macrocalls
if VERSION > v"1.1-"
Expr_cmd(x) = Expr(:macrocall, GlobalRef(Core, Symbol("@cmd")), nothing, x.val)
Expr_tcmd(x) = Expr(:macrocall, GlobalRef(Core, Symbol("@cmd")), nothing, x.val)
else
Expr_cmd(x) = Expr(:macrocall, Symbol("@cmd"), nothing, x.val)
Expr_tcmd(x) = Expr(:macrocall, Symbol("@cmd"), nothing, x.val)
end

function convert_expr(loc::Location, x::EXPR{x_Str})
    if x.args[1] isa BinarySyntaxOpCall
        mname = convert_child(loc, x, 1)
        mname.expr.args[2] = QuoteNode(Symbol("@", striploc(mname).args[2].value, "_str"))
        ret = Expr(:macrocall, mname, nothing)
    else
        name = LocExpr(loc[1], Symbol("@", x.args[1].val, "_str"))
        ret = Expr(:macrocall, name, nothing)
    end
    for i = 2:length(x.args)
        push!(ret.args, LocExpr(loc[i], x.args[i].val))
    end
    return ret
end

function convert_expr(loc::Location, x::EXPR{x_Cmd})
    ret = Expr(:macrocall, Symbol("@", x.args[1].val, "_cmd"), nothing)
    for i = 2:length(x.args)
        push!(ret.args, x.args[i].val)
    end
    return ret
end

function clear_at!(x)
    if x isa Expr && x.head == :.
        if x.args[2] isa QuoteNode && string(x.args[2].value)[1] == '@'
            x.args[2].value = Symbol(string(x.args[2].value)[2:end])
        end
        if x.args[1] isa Symbol && string(x.args[1])[1] == '@'
            x.args[1] = Symbol(string(x.args[1])[2:end])
        else
            clear_at!(x.args[1])
        end
    end
end

function convert_expr(loc::Location, x::EXPR{MacroCall})
    ret = Expr(:macrocall)
    for i in 1:length(x.args)
        if !(x[i] isa PUNCTUATION)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    insert!(ret.args, 2, nothing)
    if ret.args[1] isa Expr && ret.args[1].head == :. && string(ret.args[1].args[2].value)[1] != '@'
        clear_at!(ret.args[1])
        ret.args[1].args[2] = QuoteNode(Symbol(string('@', ret.args[1].args[2].value)))
    end
    ret
end
"""
    remlineinfo!(x)
Removes line info expressions. (i.e. Expr(:line, 1))
"""
function remlineinfo!(x)
    if isa(x, Expr)
        if x.head == :macrocall && x.args[2] != nothing
            id = findall(map(x -> (isa(x, Expr) && x.head == :line) || (@isdefined(LineNumberNode) && x isa LineNumberNode), x.args))
            deleteat!(x.args, id)
            for j in x.args
                remlineinfo!(j)
            end
            insert!(x.args, 2, nothing)
        else
            id = findall(map(x -> (isa(x, Expr) && x.head == :line) || (@isdefined(LineNumberNode) && x isa LineNumberNode), x.args))
            deleteat!(x.args, id)
            for j in x.args
                remlineinfo!(j)
            end
        end
        if x.head == :elseif && x.args[1] isa Expr && x.args[1].head == :block && length(x.args[1].args) == 1
            x.args[1] = x.args[1].args[1]
        end
    end
    x
end

convert_expr(loc::Location, x::EXPR{Quotenode}) = QuoteNode(Expr(x.args[end]))

function convert_expr(loc::Location, x::EXPR{Call})
    ret = Expr(:call)
    for i in 1:length(x.args)
        a = x.args[i]
        if a isa EXPR{Parameters}
            insert!(ret.args, 2, convert_child(loc, x, i))
        elseif !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Braces})
    ret = Expr(:braces)
    for i in 1:length(x.args)
        a = x.args[i]
        if a isa EXPR{Parameters}
            insert!(ret.args, 1, convert_child(loc, x, i))
        elseif !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

# Definitiions
convert_expr(loc::Location, x::EXPR{Struct}) = Expr(:struct, false, convert_child(loc, x, 2), convert_child(loc, x, 3))
convert_expr(loc::Location, x::EXPR{Mutable}) = length(x.args) == 4 ? Expr(:struct, true, convert_child(loc, x, 2), convert_child(loc, x, 3)) : Expr(:struct, true, convert_child(loc, x, 3), convert_child(loc, x, 4))
convert_expr(loc::Location, x::EXPR{Abstract}) = length(x.args) == 2 ? Expr(:abstract, convert_child(loc, x, 2)) : Expr(:abstract, convert_child(loc, x, 3))
convert_expr(loc::Location, x::EXPR{Primitive}) = Expr(:primitive, convert_child(loc, x, 3), convert_child(loc, x, 4))

function convert_expr(loc::Location, x::EXPR{FunctionDef})
    ret = Expr(:function)
    for i in 1:length(x.args)
        if !(x[i] isa PUNCTUATION || x[i] isa KEYWORD)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end
function convert_expr(loc::Location, x::EXPR{Macro})
    if length(x.args) == 3
        Expr(:macro, convert_child(loc, x, 2))
    else
        Expr(:macro, convert_child(loc, x, 2), convert_child(loc, x, 3))
    end
end
convert_expr(loc::Location, x::EXPR{ModuleH}) = Expr(:module, true, convert_child(loc, x, 2), convert_child(loc, x, 3))
convert_expr(loc::Location, x::EXPR{BareModule}) = Expr(:module, false, convert_child(loc, x, 2), convert_child(loc, x, 3))

# Control Flow

function convert_expr(loc::Location, x::EXPR{If})
    ret = Expr(:if)
    iselseif = false
    n = length(x.args)
    i = 0
    while i < n
        i += 1
        a = x.args[i]
        if a isa KEYWORD && a.kind == Tokens.ELSEIF
            i += 1
            r1 = convert_child(loc, x, i, 1)
            r1.expr.head = :elseif
            push!(ret.args, r1)
        elseif !(a isa PUNCTUATION || a isa KEYWORD)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Try})
    ret = Expr(:try)
    for i in 1:length(x.args)
        if !(x[i] isa PUNCTUATION || x[i] isa KEYWORD)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Let})
    ret = Expr(:let)
    if length(x.args) == 3
        push!(ret.args, Expr(:block))
        push!(ret.args, convert_child(loc, x, 2))
        return ret
    elseif x.args[2] isa EXPR{Block}
        arg = Expr(:block)
        for (i, a) in enumerate(x.args[2].args)
            if !(a isa PUNCTUATION)
                push!(arg.args, fix_range(loc[2, i], a))
            end
        end
        push!(ret.args, arg)
    else
        push!(ret.args, fix_range(loc[2], x.args[2]))
    end
    push!(ret.args, convert_child(loc, x, 3))
    ret
end

function convert_expr(loc::Location, x::EXPR{Do})
    Expr(:do, convert_child(loc, x, 1), Expr(:->, convert_child(loc, x, 3), convert_child(loc, x, 4)))
end

# Loops
convert_expr(loc::Location, x::EXPR{Outer}) = Expr(:outer, convert_child(loc, x, 2))

function convert_expr(loc::Location, x::EXPR{For})
    ret = Expr(:for)
    if x.args[2] isa EXPR{Block}
        arg = Expr(:block)
        for (i, a) in enumerate(x.args[2].args)
            if !(a isa PUNCTUATION)
                push!(arg.args, fix_range(loc[2, i], a))
            end
        end
        push!(ret.args, arg)
    else
        push!(ret.args, fix_range(loc[2], x.args[2]))
    end
    push!(ret.args, convert_child(loc, x, 3))
    ret
end

function convert_expr(loc::Location, x::EXPR{While})
    ret = Expr(:while)
    for i in 1:length(x.args)
        if !(x[i] isa PUNCTUATION || x[i] isa KEYWORD)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

fix_range(loc, a) = convert_child(loc, a)
function fix_range(loc, a::BinaryOpCall)
    if (is_in(a.op) || is_elof(a.op))
        LocExpr(loc, Expr(:(=), convert_child(loc, a, :arg1), convert_child(loc, a, :arg2)))
    else
        convert_child(loc, a)
    end
end

# Lists

function convert_expr(loc::Location, x::EXPR{TupleH})
    ret = Expr(:tuple)
    for (i, a) in enumerate(x.args)
        if a isa EXPR{Parameters}
            insert!(ret.args, 1, convert_child(loc[i], a))
        elseif !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc[i], a))
        end
    end
    return ret
end

function convert_expr(loc::Location, x::EXPR{Curly})
    ret = Expr(:curly)
    for (i, a) in enumerate(x.args)
        if a isa EXPR{Parameters}
            insert!(ret.args, 2, convert_child(loc[i], a))
        elseif !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc[i], a))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Vect})
    ret = Expr(:vect)
    for (i, a) in enumerate(x.args)
        if a isa EXPR{Parameters}
            pushfirst!(ret.args, convert_child(loc[i], a))
        elseif !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc[i], a))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Row})
    ret = Expr(:row)
    for (i, a) in enumerate(x.args)
        if !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc[i], a))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Hcat})
    ret = Expr(:hcat)
    for (i, a) in enumerate(x.args)
        if !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc[i], a))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Vcat})
    ret = Expr(:vcat)
    for (i, a) in enumerate(x.args)
        if !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc[i], a))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Block})
    ret = Expr(:block)
    for (i, a) in enumerate(x.args)
        if !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc[i], a))
        end
    end
    return ret
end

convert_expr(loc::Location, x::EXPR{Kw}) = Expr(:kw, convert_child(loc, x, 1), convert_child(loc, x, 3))

function convert_expr(loc::Location, x::EXPR{Parameters})
    ret = Expr(:parameters)
    for (i, a) in enumerate(x.args)
        if a isa EXPR{Parameters}
            insert!(ret.args, 2, convert_child(loc[i], a))
        elseif !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc[i], a))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Return})
    ret = Expr(:return)
    for i = 2:length(x.args)
        push!(ret.args, convert_child(loc, x, i))
    end
    ret
end

convert_expr(loc::Location, x::EXPR{InvisBrackets}) = convert_child(loc, x, 2)
convert_expr(loc::Location, x::EXPR{Begin}) = convert_child(loc, x, 2)

function convert_expr(loc::Location, x::EXPR{Quote})
    if x.args[2] isa EXPR{InvisBrackets} && (x.args[2].args[2] isa OPERATOR || x.args[2].args[2] isa LITERAL || x.args[2].args[2] isa IDENTIFIER)
        return QuoteNode(convert_child(loc, x, 2))
    else
        return Expr(:quote, convert_child(loc, x, 2))
    end
end

function convert_expr(loc::Location, x::EXPR{Global})
    ret = Expr(:global)
    if x.args[2] isa EXPR{Const}
        ret = Expr(:const, Expr(:global, convert_child(loc, x, 2, 2)))
    elseif length(x.args) == 2 && x.args[2] isa EXPR{TupleH}
        for (i, a) in enumerate(x.args[2].args)
            if !(a isa PUNCTUATION)
                push!(ret.args, convert_child(loc[2, i], a))
            end
        end
    else
        for i = 2:length(x.args)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Local})
    ret = Expr(:local)
    if x.args[2] isa EXPR{Const}
        ret = Expr(:const, Expr(:global, convert_child(loc, x, 2, 2)))
    elseif length(x.args) == 2 && x.args[2] isa EXPR{TupleH}
        for (i, a) in enumerate(x.args[2].args)
            if !(a isa PUNCTUATION)
                push!(ret.args, convert_child(loc[2, i], a))
            end
        end
    else
        for i = 2:length(x.args)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Const})
    ret = Expr(:const)
    for i = 2:length(x.args)
        push!(ret.args, convert_child(loc, x, i))
    end
    ret
end

convert_expr(loc::Location, x::EXPR{GlobalRefDoc}) = GlobalRef(Core, Symbol("@doc"))

function convert_expr(loc::Location, x::EXPR{Ref})
    ret = Expr(:ref)
    for (i, a) in enumerate(x.args)
        if a isa EXPR{Parameters}
            insert!(ret.args, 2, convert_child(loc, x, i))
        elseif !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{TypedHcat})
    ret = Expr(:typed_hcat)
    for (i, a) in enumerate(x.args)
        if !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{TypedVcat})
    ret = Expr(:typed_vcat)

    for (i, a) in enumerate(x.args)
        if a isa EXPR{Parameters}
            insert!(ret.args, 2, convert_child(loc, x, i))
        elseif !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Comprehension})
    ret = Expr(:comprehension)
    for (i, a) in enumerate(x.args)
        if !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Flatten})
    iters, body = get_inner_gen(loc[1], x[1])
    (l, i) = popfirst!(iters)
    ex = Expr(:generator, convert_child(body...), convert_iter_assign(l, i))
    for (l, i) in iters
        ex = Expr(:generator, ex, convert_iter_assign(l, i))
        ex = Expr(:flatten, ex)
    end
    return ex
end

get_inner_gen(loc, x::EXPR{Flatten}, iters = []) = get_inner_gen(loc[1], x[1], iters)
function get_inner_gen(loc, x::EXPR{Generator}, iters = [])
    push!(iters, (loc[3], x.args[3]))
    if x.args[1] isa EXPR{Generator} || x.args[1] isa EXPR{Flatten}
        get_inner_gen(loc[1], x.args[1], iters)
    else
        iters, (loc[1], x.args[1])
    end
end

function convert_expr(loc::Location, x::EXPR{Generator})
    ret = Expr(:generator, convert_child(loc, x, 1))
    for i = 3:length(x.args)
        a = x.args[i]
        if !(a isa PUNCTUATION)
            push!(ret.args, convert_iter_assign(loc[i], a))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Filter})
    ret = Expr(:filter)
    for (i, a) in enumerate(x.args)
        if !(is_if(a) || a isa PUNCTUATION)
            push!(ret.args, convert_iter_assign(loc[i], a))
        end
    end
    ret
end

function convert_iter_assign(loc, a)
    if a isa BinaryOpCall && (is_in(a.op) || is_elof(a.op))
        return LocExpr(loc, Expr(:(=), convert_child(loc, a, :arg1), convert_child(loc, a, :arg2)))
    else
        return convert_child(loc, a)
    end
end

function convert_expr(loc::Location, x::EXPR{TypedComprehension})
    ret = Expr(:typed_comprehension)
    for (i, a) in enumerate(x.args)
        if !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{Export})
    ret = Expr(:export)
    for i = 2:length(x.args)
        a = x.args[i]
        if !(a isa PUNCTUATION)
            push!(ret.args, convert_child(loc, x, i))
        end
    end
    ret
end

function _get_import_block(loc, x, i, ret)
    while is_dot(x.args[i + 1])
        i += 1
        push!(ret.args, :.)
    end
    while i < length(x.args) && !(is_comma(x.args[i + 1]))
        i += 1
        a = x.args[i]
        if !(a isa PUNCTUATION) && !(is_dot(a) || is_colon(a))
            push!(ret.args, convert_child(loc, x, i))
        end
    end

    return i
end

convert_expr(loc::Location, x::EXPR{Import}) = expr_import(loc, x, :import)
convert_expr(loc::Location, x::EXPR{ImportAll}) = expr_import(loc, x, :importall)
convert_expr(loc::Location, x::EXPR{Using}) = expr_import(loc, x, :using)

function expr_import(loc, x, kw)
    col = findall(a-> a isa OPERATOR && precedence(a) == ColonOp, x.args)
    comma = findall(is_comma, x.args)

    header = []
    args = [Expr(:.)]
    i = 1 #skip keyword
    while i < length(x.args)
        i+=1
        a = x.args[i]
        if is_colon(a)
            push!(header, popfirst!(args))
            push!(args, Expr(:.))
        elseif is_comma(a)
            push!(args, Expr(:.))
        elseif !(a isa PUNCTUATION)
            push!(last(args).args, convert_child(loc, x, i))
        end
    end
    if isempty(header)
        return Expr(kw, args...)
    else
        return Expr(kw, Expr(:(:), header..., args...))
    end
end

function convert_expr(loc::Location, x::EXPR{FileH})
    ret = Expr(:file)
    for (i, a) in enumerate(x.args)
        push!(ret.args, convert_child(loc, x, i))
    end
    ret
end

function convert_expr(loc::Location, x::EXPR{StringH})
    ret = Expr(:string)
    for (i, a) in enumerate(x.args)
        if a isa UnarySyntaxOpCall
            a = a.arg2
            push!(ret.args, convert_child(loc, x, i, :arg2))
            continue
        elseif a isa LITERAL && a.kind == Tokens.STRING
            if span(a) == 0 || ((i == 1 || i == length(x.args)) && span(a) == 1) || isempty(a.val)
                continue
            end
        end
        push!(ret.args, convert_child(loc, x, i))
    end
    ret
end

const UNICODE_OPS_REVERSE = Dict{Tokenize.Tokens.Kind,Symbol}()
for (k, v) in Tokenize.Tokens.UNICODE_OPS
    UNICODE_OPS_REVERSE[v] = Symbol(k)
end

UNICODE_OPS_REVERSE[Tokens.EQ] = :(=)
UNICODE_OPS_REVERSE[Tokens.PLUS_EQ] = :(+=)
UNICODE_OPS_REVERSE[Tokens.MINUS_EQ] = :(-=)
UNICODE_OPS_REVERSE[Tokens.STAR_EQ] = :(*=)
UNICODE_OPS_REVERSE[Tokens.FWD_SLASH_EQ] = :(/=)
UNICODE_OPS_REVERSE[Tokens.FWDFWD_SLASH_EQ] = :(//=)
UNICODE_OPS_REVERSE[Tokens.OR_EQ] = :(|=)
UNICODE_OPS_REVERSE[Tokens.CIRCUMFLEX_EQ] = :(^=)
UNICODE_OPS_REVERSE[Tokens.DIVISION_EQ] = :(÷=)
UNICODE_OPS_REVERSE[Tokens.REM_EQ] = :(%=)
UNICODE_OPS_REVERSE[Tokens.LBITSHIFT_EQ] = :(<<=)
UNICODE_OPS_REVERSE[Tokens.RBITSHIFT_EQ] = :(>>=)
UNICODE_OPS_REVERSE[Tokens.LBITSHIFT] = :(<<)
UNICODE_OPS_REVERSE[Tokens.RBITSHIFT] = :(>>)
UNICODE_OPS_REVERSE[Tokens.UNSIGNED_BITSHIFT] = :(>>>)
UNICODE_OPS_REVERSE[Tokens.UNSIGNED_BITSHIFT_EQ] = :(>>>=)
UNICODE_OPS_REVERSE[Tokens.BACKSLASH_EQ] = :(\=)
UNICODE_OPS_REVERSE[Tokens.AND_EQ] = :(&=)
UNICODE_OPS_REVERSE[Tokens.COLON_EQ] = :(:=)
UNICODE_OPS_REVERSE[Tokens.PAIR_ARROW] = :(=>)
UNICODE_OPS_REVERSE[Tokens.APPROX] = :(~)
UNICODE_OPS_REVERSE[Tokens.EX_OR_EQ] = :($=)
UNICODE_OPS_REVERSE[Tokens.XOR_EQ] = :(⊻=)
UNICODE_OPS_REVERSE[Tokens.RIGHT_ARROW] = :(-->)
UNICODE_OPS_REVERSE[Tokens.LAZY_OR] = :(||)
UNICODE_OPS_REVERSE[Tokens.LAZY_AND] = :(&&)
UNICODE_OPS_REVERSE[Tokens.ISSUBTYPE] = :(<:)
UNICODE_OPS_REVERSE[Tokens.ISSUPERTYPE] = :(>:)
UNICODE_OPS_REVERSE[Tokens.GREATER] = :(>)
UNICODE_OPS_REVERSE[Tokens.LESS] = :(<)
UNICODE_OPS_REVERSE[Tokens.GREATER_EQ] = :(>=)
UNICODE_OPS_REVERSE[Tokens.GREATER_THAN_OR_EQUAL_TO] = :(≥)
UNICODE_OPS_REVERSE[Tokens.LESS_EQ] = :(<=)
UNICODE_OPS_REVERSE[Tokens.LESS_THAN_OR_EQUAL_TO] = :(≤)
UNICODE_OPS_REVERSE[Tokens.EQEQ] = :(==)
UNICODE_OPS_REVERSE[Tokens.EQEQEQ] = :(===)
UNICODE_OPS_REVERSE[Tokens.IDENTICAL_TO] = :(≡)
UNICODE_OPS_REVERSE[Tokens.NOT_EQ] = :(!=)
UNICODE_OPS_REVERSE[Tokens.NOT_EQUAL_TO] = :(≠)
UNICODE_OPS_REVERSE[Tokens.NOT_IS] = :(!==)
UNICODE_OPS_REVERSE[Tokens.NOT_IDENTICAL_TO] = :(≢)
UNICODE_OPS_REVERSE[Tokens.IN] = :(in)
UNICODE_OPS_REVERSE[Tokens.ISA] = :(isa)
UNICODE_OPS_REVERSE[Tokens.LPIPE] = :(<|)
UNICODE_OPS_REVERSE[Tokens.RPIPE] = :(|>)
UNICODE_OPS_REVERSE[Tokens.COLON] = :(:)
UNICODE_OPS_REVERSE[Tokens.DDOT] = :(..)
UNICODE_OPS_REVERSE[Tokens.EX_OR] = :($)
UNICODE_OPS_REVERSE[Tokens.PLUS] = :(+)
UNICODE_OPS_REVERSE[Tokens.MINUS] = :(-)
UNICODE_OPS_REVERSE[Tokens.PLUSPLUS] = :(++)
UNICODE_OPS_REVERSE[Tokens.OR] = :(|)
UNICODE_OPS_REVERSE[Tokens.STAR] = :(*)
UNICODE_OPS_REVERSE[Tokens.FWD_SLASH] = :(/)
UNICODE_OPS_REVERSE[Tokens.REM] = :(%)
UNICODE_OPS_REVERSE[Tokens.BACKSLASH] = :(\)
UNICODE_OPS_REVERSE[Tokens.AND] = :(&)
UNICODE_OPS_REVERSE[Tokens.FWDFWD_SLASH] = :(//)
UNICODE_OPS_REVERSE[Tokens.CIRCUMFLEX_ACCENT] = :(^)
UNICODE_OPS_REVERSE[Tokens.DECLARATION] = :(::)
UNICODE_OPS_REVERSE[Tokens.CONDITIONAL] = :?
UNICODE_OPS_REVERSE[Tokens.DOT] = :(.)
UNICODE_OPS_REVERSE[Tokens.NOT] = :(!)
UNICODE_OPS_REVERSE[Tokens.PRIME] = Symbol(''')
UNICODE_OPS_REVERSE[Tokens.DDDOT] = :(...)
UNICODE_OPS_REVERSE[Tokens.TRANSPOSE] = Symbol(".'")
UNICODE_OPS_REVERSE[Tokens.ANON_FUNC] = :(->)
UNICODE_OPS_REVERSE[Tokens.WHERE] = :where


text:
module Luna

using Colors
using LinearAlgebra
using StaticArrays

export I
export Intersection
export Matrix4x4
export Object
export Ray
export Transform
export Vector4
export cross
export cube
export direction
export dot
export eye
export ispoint, isvector
export lerp
export matrix
export norm
export normalat
export normalize
export origin
export point
export reflect
export rotx, roty, rotz
export scale
export shear
export sphere
export transform
export translate
export vector

"""
    Vector4{T} <: FieldVector{4,T}

A 4-vector with elements of type `T`.
"""
struct Vector4{T} <: FieldVector{4,T}
    x::T
    y::T
    z::T
    w::T
end

"""
    Matrix4x4{T} = SMatrix{4,4,T,16} where T<:Real

A 4x4 matrix with elements of type `T`.
"""
const Matrix4x4{T} = SMatrix{4,4,T,16} where {T<:Real}

struct Ray{T}
    origin::Vector4{T}
    direction::Vector4{T}
end

struct Transform{T}
    matrix::Matrix4x4{T}
    matrix⁻¹::Matrix4x4{T}
end

"""
    Object{S,T}

An object is a thing that can be rendered where `S` is
expected to be a *symbol* tag and `T` is the numeric type
of the object's transformation.

# Remarks
Objects are always defined in unit space and depend on
their `transform` to move them into world space.

# Examples
```julia-repl
julia> t = Matrix4x4{Float64}(I) |> scale(0.5, 0.75, 0.25) |> Transform
Transform([0.5 0.0 0.0 0.0; 0.0 0.75 0.0 0.0; 0.0 0.0 0.25 0.0; 0.0 0.0 0.0 1.0])

julia> Object{:foo, Float64}(t)
Object{:foo,Float64}(Transform([0.5 0.0 0.0 0.0; 0.0 0.75 0.0 0.0; 0.0 0.0 0.25 0.0;
0.0 0.0 0.0 1.0]))
```
"""
mutable struct Object{S,T}
    transform::Transform{T}
    children::Array{Object{S,T},1}
end

struct Intersection{S,T}
    t::T
    object::Object{S,T}
end

StaticArrays.similar_type(
    ::Type{V},
    ::Type{F},
    size::Size{(4,)},
) where {V<:Vector4,F<:Real} = Vector4

Base.promote_rule(::Type{Vector4{Int64}}, ::Type{Vector4{Rational{Int64}}}) =
    Vector4{Rational{Int64}}
Base.promote_rule(::Type{Vector4{Rational{Int64}}}, ::Type{Vector4{Int64}}) =
    Vector4{Rational{Int64}}
Base.promote_rule(::Type{Vector4{Int64}}, ::Type{Vector4{Float64}}) = Vector4{Float64}
Base.promote_rule(::Type{Vector4{Float64}}, ::Type{Vector4{Int64}}) = Vector4{Float64}

function Base.show(io::IO, u::Vector4)
    if ispoint(u)
        print(io, "point($(u.x), $(u.y), $(u.z))")
    elseif isvector(u)
        print(io, "vector($(u.x), $(u.y), $(u.z))")
    else
        print(io, u)
    end
end

"""
    isless(a::Intersection, b::Intersection)

Return `true` if `a` is closer than `b`.

# Remarks
Intersections are compared using their `t` value.
"""
Base.isless(a::Intersection, b::Intersection) = isless(a.t, b.t)

"""
    lerp(a, b, t)

Linearly interpolates between `a` and `b`.

# Examples
```julia-repl
julia> lerp(0, 1, 0.5)
0.5

julia> lerp(0, 2, 0.5)
1.0

julia> lerp(-2, 2, 0.9)
1.6
```
"""
lerp(a, b, t) = (1 - t) * a + t * b

"""
    point(x, y, z)

Creates a position vector in 3-dimensional space.
"""
point(x, y, z) = Vector4(promote(x, y, z, one(x))...)

"""
    vector(x, y, z)

Creates a direction vector in 3-dimensional space.
"""
vector(x, y, z) = Vector4(promote(x, y, z, zero(x))...)

"""
    point(u)

Converts 4-vector `u` into a position vector.
"""
point(u) = point(u.x, u.y, u.z)

"""
    vector(u)

Converts 4-vector `u` into a direction vector.
"""
vector(u) = vector(u.x, u.y, u.z)

"""
    ispoint(u)

Returns `true` if 4-vector `u` is a position vector.
"""
ispoint(u) = isone(u.w)

"""
    isvector(u)

Returns `true` if 4-vector `u` is a direction vector.
"""
isvector(u) = iszero(u.w)

"""
    cross(u::Vector4, v::Vector4)

Computes the cross product of two 4-vectors.

# Remarks
The `w` component is ignored and the result is returned as
a new direction vector.

# Examples
```julia-repl
julia> u = vector(1,1,0)
4-element Vector4{Int64} with indices SOneTo(4):
 1
 1
 0
 0

julia> v = vector(0,0,1)
4-element Vector4{Int64} with indices SOneTo(4):
 0
 0
 1
 0

julia> cross(u, v)
4-element Vector4{Int64} with indices SOneTo(4):
  1
 -1
  0
  0
```
"""
function LinearAlgebra.cross(u::Vector4, v::Vector4)
    x = u.y * v.z - u.z * v.y
    y = u.z * v.x - u.x * v.z
    z = u.x * v.y - u.y * v.x
    vector(x, y, z)
end

"""
    reflect(u, n)

Reflects vector `u` over vector `n`.
"""
reflect(u, n) = u - 2n * dot(u, n)

"""
    Ray(o, d)

Creates a ray with origin `o` and direction `d`.
"""
Ray(o, d) = Ray(promote(o, d)...)

"""
    origin(r)

Returns the origin position vector for ray `r`.
"""
origin(r) = r.origin

"""
    direction(r)

Returns the direction vector for ray `r`.
"""
direction(r) = r.direction

"""
    (r::Ray)(t)

Calculates a position vector at `t` along ray `r`.

# Examples
```julia-repl
julia> r(0)
4-element Vector4{Int64} with indices SOneTo(4):
 0
 0
 0
 1

julia> r(1)
4-element Vector4{Int64} with indices SOneTo(4):
 0
 0
 1
 1

julia> r(0.5)
4-element StaticArrays.SArray{Tuple{4},Float64,1,4} with indices SOneTo(4):
 0.0
 0.0
 0.5
 1.0
```
"""
(r::Ray)(t) = origin(r) + t * direction(r)

"""
    *(m, r::Ray)
    *(r::Ray, m)

Returns a new ray with `m` applied to origin and direction of `r`.

# Examples
```
julia> m = Matrix4x4{Float64}(I) |> scale(0.5, 0.75, 0.3)
4×4 StaticArrays.SArray{Tuple{4,4},Float64,2,16} with indices SOneTo(4)×SOneTo(4):
 0.5  0.0   0.0  0.0
 0.0  0.75  0.0  0.0
 0.0  0.0   0.3  0.0
 0.0  0.0   0.0  1.0

julia> r = Ray(point(0,0,0), vector(0,0,1))
Ray{Int64}(point(0, 0, 0), vector(0, 0, 1))

julia> r * m
Ray{Float64}(point(0.0, 0.0, 0.0), vector(0.0, 0.0, 0.3))

julia> m * r
Ray{Float64}(point(0.0, 0.0, 0.0), vector(0.0, 0.0, 0.3))
```
"""
Base.:*(m, r::Ray) = Ray(m * origin(r), m * direction(r))
Base.:*(r::Ray, m) = m * r

"""
    eye()

Creates a 4x4 identity matrix of of `Float64` values.
"""
eye() = Matrix4x4{Float64}(I)

"""
    translate(dx, dy, dz)

Creates a translation operation.

# Examples
```
julia> Matrix4x4{Float64}(I) |> translate(1, -2, 3)
4×4 StaticArrays.SArray{Tuple{4,4},Float64,2,16} with indices SOneTo(4)×SOneTo(4):
 1.0  0.0  0.0   1.0
 0.0  1.0  0.0  -2.0
 0.0  0.0  1.0   3.0
 0.0  0.0  0.0   1.0
```
"""
translate(dx, dy, dz) = m -> _translate(dx, dy, dz) * m

"""
    scale(sx, sy, sz)

Creates a scaling operation.

# Examples
```julia-repl
julia> Matrix4x4{Float64}(I) |> scale(0.5, 0.5, 0.75)
4×4 StaticArrays.SArray{Tuple{4,4},Float64,2,16} with indices SOneTo(4)×SOneTo(4):
 0.5  0.0  0.0   0.0
 0.0  0.5  0.0   0.0
 0.0  0.0  0.75  0.0
 0.0  0.0  0.0   1.0
```
"""
scale(sx, sy, sz) = m -> _scale(sx, sy, sz) * m

"""
    rotx(rad)

Creates a rotation operation along the x-axis.

# Examples
```julia-repl
julia> Matrix4x4{Float64}(I) |> rotx(π)
4×4 StaticArrays.SArray{Tuple{4,4},Float64,2,16} with indices SOneTo(4)×SOneTo(4):
 1.0   0.0           0.0          0.0
 0.0  -1.0          -1.22465e-16  0.0
 0.0   1.22465e-16  -1.0          0.0
 0.0   0.0           0.0          1.0

 julia> Matrix4x4{Float64}(I) |> rotx(π/5)
4×4 StaticArrays.SArray{Tuple{4,4},Float64,2,16} with indices SOneTo(4)×SOneTo(4):
 1.0  0.0        0.0       0.0
 0.0  0.809017  -0.587785  0.0
 0.0  0.587785   0.809017  0.0
 0.0  0.0        0.0       1.0
```
"""
rotx(rad) = m -> _rotx(rad) * m

"""
    roty(rad)

Creates a rotation operation along the y-axis.

See also [`rotx`](@ref) for examples.
"""
roty(rad) = m -> _roty(rad) * m

"""
    rotz(rad)

Creates a rotation operation along the z-axis.

See also [`rotx`](@ref) for examples.
"""
rotz(rad) = m -> _rotz(rad) * m

"""
    shear(xy, xz, yx, yz, zx, zy)

Creates a shearing operation.

# Examples
```julia-repl
julia> Matrix4x4{Float64}(I) |> shear(1, 2, 3, -1, -2, -3)
4×4 StaticArrays.SArray{Tuple{4,4},Float64,2,16} with indices SOneTo(4)×SOneTo(4):
  1.0   1.0   2.0  0.0
  3.0   1.0  -1.0  0.0
 -2.0  -3.0   1.0  0.0
  0.0   0.0   0.0  1.0
```
"""
shear(xy, xz, yx, yz, zx, zy) = m -> _shear(xy, xz, yx, yz, zx, zy) * m

"""
    Transform(m)

Creates a transform for matrix `m`.

# Remarks
A `Transform` is basically just a memoized inversion matrix so that calling
`inv(t::Transform)` returns the cached inversion instead of recalculating it.
This implies `inv(t) == inv(matrix(t))`.

# Examples
```julia-repl
julia> m = Matrix4x4{Float64}(I) |> scale(0.5, 0.25, 0.75)
4×4 StaticArrays.SArray{Tuple{4,4},Float64,2,16} with indices SOneTo(4)×SOneTo(4):
 0.5  0.0   0.0   0.0
 0.0  0.25  0.0   0.0
 0.0  0.0   0.75  0.0
 0.0  0.0   0.0   1.0

 julia> m |> Transform
Transform([0.5 0.0 0.0 0.0; 0.0 0.25 0.0 0.0; 0.0 0.0 0.75 0.0; 0.0 0.0 0.0 1.0])

julia> inv(t) == inv(matrix(t))
true
```
"""
Transform(m) = Transform(promote(m, inv(m))...)

function Base.show(io::IO, t::Transform)
    print(io, "Transform($(matrix(t)))")
end

"""
    matrix(t)

Returns the transform matrix associated with transformation `t`.

# Examples
```julia-repl
julia> t = Matrix4x4{Float64}(I) |> scale(0.5, π/5, -1) |> Transform
Transform([0.5 0.0 0.0 0.0; 0.0 0.6283185307179586 0.0 0.0; 0.0 0.0 -1.0 0.0; 0.0 0.0 0.0 1.0])

julia> matrix(t)
4×4 StaticArrays.SArray{Tuple{4,4},Float64,2,16} with indices SOneTo(4)×SOneTo(4):
 0.5  0.0        0.0  0.0
 0.0  0.628319   0.0  0.0
 0.0  0.0       -1.0  0.0
 0.0  0.0        0.0  1.0
```
"""
matrix(t) = t.matrix

"""
    inv(t::Transform)

Returns the inverse transform matrix.
"""
Base.inv(t::Transform) = t.matrix⁻¹

"""
    eltype(::Transform{T})

Returns the element type of the matrices wrapped by a transform.

# Remarks
Transforms depend on the `inv` method for matrices which creates
matrices with element type `Float64` or `Float32` depending on
your system architecture. This means that in practice this method
will always return `Float64` (or `Float32`).
"""
Base.eltype(::Transform{T}) where {T} = T

"""
    transform(obj)

Returns the transformation associated with object `obj`.
"""
transform(obj) = obj.transform

"""
    sphere(t)

Returns a unit sphere with transformation `t`.
"""
sphere(t) = Object{:sphere,eltype(t)}(t, [])

"""
    cube(t)

Returns a unit cube with transformation `t`.
"""
cube(t) = Object{:cube,eltype(t)}(t, [])

"""
    group(t, children)

Returns a group with transformation `t` and `children`.
"""
group(t, children) = Objectr{:group,eltype(t)}(t, children)

"""
    Intersection(t, obj::Object{S,T}) where {S,T}

Creates a new intersection.
"""
Intersection(t, obj::Object{S,T}) where {S,T} = Intersection(convert(T, t), obj)

"""
    intersect(obj::Object{:sphere,T}, ray) where T

Returns the intersections between a sphere and ray.
"""
function Base.intersect(obj::Object{:sphere,T}, ray) where {T}
    ray *= inv(transform(obj))
    oc = origin(ray) - point(0, 0, 0)
    a = dot(direction(ray), direction(ray))
    b = 2 * dot(direction(ray), oc)
    c = dot(oc, oc) - 1
    d = b * b - 4 * a * c
    if iszero(d)
        t = -b / 2a
        return [Intersection(t, obj)]
    elseif d > 0
        sqrtd = √d
        t1 = (-b - sqrtd) / 2a
        t2 = (-b + sqrtd) / 2a
        return [Intersection(t1, obj), Intersection(t2, obj)]
    end
    return empty([], Intersection{:sphere,T})
end

"""
    normalat(obj::Object{:sphere,T}, worldpoint)

Returns the world normal of `obj` at `worldpoint`.
"""
function normalat(obj::Object{:sphere}, worldpoint)
    minv = inv(transform(obj))
    objectpoint = minv * worldpoint
    objectnormal = objectpoint - point(0, 0, 0)
    worldnormal = transpose(minv) * objectnormal
    return normalize(vector(worldnormal))
end

_translate(x, y, z) = @SMatrix[
    1 0 0 x
    0 1 0 y
    0 0 1 z
    0 0 0 1
]

_scale(x, y, z) = @SMatrix[
    x 0 0 0
    0 y 0 0
    0 0 z 0
    0 0 0 1
]

_rotx(r) = @SMatrix[
    1 0 0 0
    0 cos(r) -sin(r) 0
    0 sin(r) cos(r) 0
    0 0 0 1
]

_roty(r) = @SMatrix[
    cos(r) 0 sin(r) 0
    0 1 0 0
    -sin(r) 0 cos(r) 0
    0 0 0 1
]

_rotz(r) = @SMatrix[
    cos(r) -sin(r) 0 0
    sin(r) cos(r) 0 0
    0 0 1 0
    0 0 0 1
]

_shear(xy, xz, yx, yz, zx, zy) = @SMatrix[
    1 xy xz 0
    yx 1 yz 0
    zx zy 1 0
    0 0 0 1
]

end # module

max_stars_repo_name: basp/luna

max_issues_repo_name: basp/luna

max_forks_repo_name: basp/luna


text:
#=

  Euler #6 in Julia.

  Problem 6
  """
  The sum of the squares of the first ten natural numbers is,
  1^(2) + 2^(2) + ... + 10^(2) = 385

  The square of the sum of the first ten natural numbers is,
  (1 + 2 + ... + 10)^(2) = 55^(2) = 3025

  Hence the difference between the sum of the squares of the first ten
  natural numbers and the square of the sum is 3025 − 385 = 2640.

  Find the difference between the sum of the squares of the first one
  hundred natural numbers and the square of the sum.
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.00000693s
function euler6a()
    s1 = sum(1:100)^2
    s2 = sum((1:100).|>i->i^2)
    return s1-s2;
end

# 0.00000085s
function euler6b()
    sum(1:100)^2 - sum((1:100).|>i->i^2)
end

# run_euler(euler6a)
run_euler(euler6b)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
my_f(x, y) = 2x + 3y
derivative_of_my_f(x,y) = ForwardDiff.derivative(x->my_f(x,y),x)


text:
## -----------  School types

abstract type AbstractSchoolGroups{I1 <: Integer} end

# Map :HSG => 1 etc
function ed_idx(s :: AbstractSchoolGroups{I1}, ed :: Symbol) where I1
    edIdx = findfirst(ed_symbols(s) .== ed);
    @assert  (!isnothing(edIdx)  && edIdx > 0)  "Invalid $ed"
    return edIdx
end

ed_idx(s :: AbstractSchoolGroups{I1}, ed :: Integer) where I1 = I1(ed);

# String labels
ed_labels(s :: AbstractSchoolGroups{I1}) where I1 = string.(ed_symbols(s));
ed_label(s :: AbstractSchoolGroups{I1}, ed :: Symbol) where I1 =
    ed_labels(s)[ed_idx(s, ed)];
ed_label(s :: AbstractSchoolGroups{I1}, ed :: Integer) where I1 = ed_labels(s)[ed];
ed_symbol(s :: AbstractSchoolGroups{I1}, idx) where I1 = ed_symbols(s)[idx];

# Suffixes of the form "_1HSG" for files
ed_suffixes(s :: AbstractSchoolGroups{I1}) where I1 =
    [ed_suffix(s, j)  for j = 1 : n_school(s)];

function ed_suffix(s :: AbstractSchoolGroups{I1}, ed) where I1
    edIdx = ed_idx(s, ed);
    edStr = ed_label(s, ed);
    return "_$edIdx$edStr"
end

# ------- HSG, CD, CG
struct Schooling3{I1} <: AbstractSchoolGroups{I1}

end

Base.show(io :: IO, s :: Schooling3{I1}) where I1 =
    print(io, "Schooling: HSG, SC, CG");

StructLH.describe(s :: Schooling3{I1}) where I1 =
    ["Schooling:"  "HSG, SC, CG"];

n_school(s :: Schooling3{I1}) where I1 = I1(3);

ed_symbols(s :: Schooling3) = [:HSG, :SC, :CG];

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