src/LinSolve.jl

"""
Linear solvers for small problems.
"""
module LinSolve

using ..Metaprogramming

export linsolve,
       TupleMatrix

const TupleMatrix{M,N} = NTuple{M, NTuple{N, Float64}}

"""
`linsolve(AB :: TupleMatrix{M,N}, :: Type{TupleMatrix{M, K}}) :: TupleMatrix{M, K}`

where

`TupleMatrix{M, N} = NTuple{M, NTuple{N, Float64}}`

“Immutable” Gaussian elimination on tuples: solve AX=B for X,
Both A and B are stored in AB. The second type parameter indicates the size of B.
"""
@polly function linsolve₀(AB :: TupleMatrix{M,N}, :: Val{K}) :: TupleMatrix{M, K} where {N,M,K}
    @assert(M == N - K)

    k = 0

    # Convert to row-echelon form
    for h = 1:(M-1)
        # Find pivotable column (has some non-zero entries in rows ≥ h)
        v = 0.0
        î = h
        while k ≤ N-1 && v == 0
            k = k + 1
            v = abs(AB[h][k])
            # Find row ≥ h of maximum absolute value in this column
            for i=(h+1):M
                local v′ = abs(AB[i][k])
                if v′ > v
                    î = i
                    v = v′
                end
            end
        end

        if v > 0
            AB = (
                AB[1:(h-1)]...,
                AB[î],
                (let ĩ = (i==î ? h : i), h̃ = î, f = AB[ĩ][k] / AB[h̃][k]
                    ((0.0 for _ = 1:k)..., (AB[ĩ][j]-AB[h̃][j]*f for j = (k+1):N)...,)
                end for i = (h+1):M)...
            )
        end
    end

    # Solve UAX=UB for X where UA with U presenting the transformations above an
    # upper triangular matrix.
    X = ()
    for i=M:-1:1
        r = .+(AB[i][M+1:end], (-AB[i][j].*X[j-i] for j=i+1:M)...)./AB[i][i]
        X = (r, X...)
    end
    return X
end

@inline function linsolve₀(AB :: TupleMatrix{M,N}) :: NTuple{M,Float64} where {N,M}
    X = linsolve₀(AB, Val(1))
    return ((x for (x,) ∈ X)...,)
end

#@generated function linsolve₁(AB :: TupleMatrix{M,N}, :: Type{TupleMatrix{M, K}}) :: TupleMatrix{M, K} where {N,M,K}
function generate_linsolve(AB :: Symbol, M :: Int, N :: Int, K :: Int)
    @assert(M == N - K)
    # The variables of ABN collect the stepwise stages of the transformed matrix.
    # Initial i-1 entries of ABN[i] are never used, as the previous steps have already
    # finalised the corresponding rows, but are included as “missing” (at compile time)
    # for the sake of indexing clarity. The M-1:th row-wise step finalises the row-echelon
    # form, so ABN has M-1 rows itself.
    step_ABN(step) = ((missing for i=1:step-1)...,
                      (gensym("ABN_$(step)_$(i)") for i ∈ step:M)...,)
    ABN = ((step_ABN(step) for step=1:M-1)...,)
    # UAB “diagonally” refers to ABN to collate the final rows of the transformed matrix UAB.
    # Since the M-1:th row-wise step finalises the row-echelon form, the last row comes already
    # from the M-1:th step.
    UAB = ((ABN[i][i] for i=1:M-1)..., ABN[M-1][M])
    # The variables of X collect the rows of the solution to AX=B.
    X = ((gensym("X_$(i)") for i ∈ 1:M)...,)

    # Convert to row-echelon form. On each step we strip leading zeroes from ABN.
    # In the end UAB should be upper-triangular.
    convert_row(ABNout, h, ABNin) = quote
        # Find pivotable column (has some non-zero entries in rows ≥ h)
        $(ABNout[h]) = $(ABNin[h])
        local v = abs($(ABNout[h])[1])
        # Find row ≥ h of maximum absolute value in this column
        $(sequence_exprs(:(
            let v′ = abs($(ABNin[i])[1])
                if v′ > v
                    $(ABNout[h]), $(ABNin[i]) = $(ABNin[i]), $(ABNout[h])
                    v = v′
                end
            end
        ) for i=(h+1):M))

        $(lift_exprs(ABNout[h+1:M])) = v > 0 ? ( $(lift_exprs( :(
            # Transform
            $(ABNin[i])[2:$N-$h+1] .- $(ABNout[h])[2:$N-$h+1].*( $(ABNin[i])[1] / $(ABNout[h])[1]) 
        ) for i=h+1:M)) ) : ( $(lift_exprs( :(
            # Strip leading zeroes
            $(ABNin[i])[2:$N-$h+1]
        ) for i=h+1:M)) )
    end

    # Solve UAX=UB for X where UA with U presenting the transformations above an
    # upper triangular matrix.
    solve_row(UAB, i) = :(
        $(X[i]) = $(lift_exprs( :(
            +($(UAB[i])[$M-$i+1+$k], $(( :( -$(UAB[i])[$j-$i+1]*$(X[j])[$k] ) for j=i+1:M)...))
        ) for k=1:K )) ./ $(UAB[i])[1]
    )

    return X, quote
        $(lift_exprs(ABN[1][i] for i=1:M)) = $(lift_exprs( :( $AB[$i] ) for i=1:M))
        $(convert_row(ABN[1], 1, ABN[1]))
        $((convert_row(ABN[h], h, ABN[h-1]) for h = 2:(M-1))...)
        $((solve_row(UAB, i) for i=M:-1:1)...)
    end
end

@inline @generated function linsolve₁(AB :: TupleMatrix{M,N}, :: Val{K}) :: TupleMatrix{M, K} where {N,M,K}
    X, solver = generate_linsolve(:AB, M, N, K)
    return quote
        $solver
        return $(lift_exprs( X[i] for i=1:M ))
    end
end

@inline @generated function linsolve₁(AB :: TupleMatrix{M,N}) :: NTuple{M,Float64} where {N,M}
    X, solver = generate_linsolve(:AB, M, N, 1)
    return quote
        $solver
        return $(lift_exprs( :( $(X[i])[1] ) for i=1:M ))
    end
end

const linsolve = linsolve₁

function tuplify(M, N, A, b)
    ((((A[i][j] for j=1:N)..., b[j]) for i=1:M)...,)
end


function compare(; dim=5, n_matrices=10000, n_testvectors=100)
    testmatrices=[]
    testvectors=[]
    while length(testmatrices)<n_matrices
        A=randn(dim, dim)
        push!(testmatrices, A)
    end
    test_vectors=[randn(dim) for _ = 1:n_testvectors]


    function evaluate(fn)
        for A ∈ testmatrices
            for b ∈ testvectors
                fn(A, b)
            end
        end
    end

    function evaluate_and_report(fn, name)
        printstyled("Evaluating $name…\n", color=:cyan)
        @time evaluate(fn)
    end

    evaluate_and_report("tuple-linsolve, ungenerated") do A, b
        linsolve₀(tuplify(A, b))
    end

    evaluate_and_report("tuple-linsolve, generated") do A, b
        linsolve₁(tuplify(A, b))
    end

    evaluate_and_report("backslash") do A, b
        A \ b
    end
end

end # module