using MAT

function readmat(file, names...)
    n = length(names)
    objects = Array{Any}(undef,n)
    f1 = matopen(file)
    for k = 1:n
        objects[k] = read(f1, names[k])
    end
    close(f1)
    return tuple(objects...)
end

"""
    test_small()

Test special cases.
"""
function test_small()
    # dimension 0
    B = spzeros(0, 0)
    blu = basiclu_object(0)
    factorize(blu, B)
    rhs = zeros(0)
    lhs = solve(blu, rhs, 'N')
    lhs = solve(blu, rhs, 'T')
    rhs = spzeros(0)
    lhs = solve(blu, rhs, 'N')
    lhs = solve(blu, rhs, 'T')
    lhs = solve_for_update(blu, rhs, getsol=true)

    # dimension 1
    B = sparse(1.0I, 1, 1)
    blu = basiclu_object(1)
    factorize(blu, B)
    rhs = [2.0]
    lhs = solve(blu, rhs, 'N')
    @assert norm(lhs-rhs, Inf) == 0.0
    lhs = solve(blu, rhs, 'T')
    @assert norm(lhs-rhs, Inf) == 0.0
    rhs = SparseVector{Float64,Int64}(1, [1], [2.0])
    lhs = solve(blu, rhs, 'N')
    @assert norm(lhs-rhs, Inf) == 0.0
    lhs = solve(blu, rhs, 'T')
    @assert norm(lhs-rhs, Inf) == 0.0
    lhs = solve_for_update(blu, rhs, getsol=true)
    @assert norm(lhs-rhs, Inf) == 0.0

    # dimension 1 singular
    B = spzeros(1, 1)
    blu = basiclu_object(1)
    factorize(blu, B, check=false)
    info = getinfo(blu)
    @assert info.nPivot == 0
    # basiclu has inserted a slack column for the zero pivot
    rhs = [2.0]
    lhs = solve(blu, rhs, 'N')
    @assert norm(lhs-rhs, Inf) == 0.0
    lhs = solve(blu, rhs, 'T')
    @assert norm(lhs-rhs, Inf) == 0.0
    rhs = SparseVector{Float64,Int64}(1, [1], [2.0])
    lhs = solve(blu, rhs, 'N')
    @assert norm(lhs-rhs, Inf) == 0.0
    lhs = solve(blu, rhs, 'T')
    @assert norm(lhs-rhs, Inf) == 0.0
    lhs = solve_for_update(blu, rhs, getsol=true)
    @assert norm(lhs-rhs, Inf) == 0.0
end

"""
    test_factorize(testdir::String, trans::Bool=false)

For all `*.mat` files in `testdir` read matrix `B` and factorize it. Monitor
residual of factorization and forward/backward solves.

`trans` specifies whether `B` or its transposed is factorized.
"""
function test_factorize(testdir::String, trans::Bool=false)
    files = readdir(testdir)
    for f in files
        if length(f) < 4 || f[end-3:end] != ".mat"
            continue
        end
        @printf(" %-24s", f)
        B, = readmat(string(testdir, f), "B")
        Bt = sparse(transpose(B))
        if trans
            (B,Bt) = (Bt,B)
        end
        m = size(B,1)
        blu = basiclu_object(m)
        factorize(blu, B)
        L,U,p,q = getfactors(blu)
        res = norm(L*U-B[p,q], Inf)

        # test dense solve
        rhs = ones(m)
        lhs = solve(blu, rhs, 'N')
        res = max(res, norm(B*lhs-rhs, Inf))
        lhs = solve(blu, rhs, 'T')
        res = max(res, norm(Bt*lhs-rhs, Inf))

        # test sparse solve
        rhs = SparseVector{Float64,Int64}(m, [1;m], [1.0;1.0])
        lhs = solve(blu, rhs, 'N')
        res = max(res, norm(B*lhs-rhs, Inf))
        lhs = solve(blu, rhs, 'T')
        res = max(res, norm(Bt*lhs-rhs, Inf))

        # test solve for update
        lhs = solve_for_update(blu, rhs, getsol=true)
        res = max(res, norm(B*lhs-rhs, Inf))
        lhs = solve_for_update(blu, m-1, getsol=true)
        rhs = zeros(m); rhs[m-1] = 1
        res = max(res, norm(Bt*lhs-rhs, Inf))
        @printf("%.2e\n", res)
    end
    nothing
end

"""
    test_update(testdir::String)

For all `*.mat` files in `testdir` read matrix `A` and the pivot sequence.
Starting from the slack basis apply pivot operations. Monitor residuals to
forward/backward solves after each 100 updates.
"""
function test_update(testdir::String)
    files = readdir(testdir)
    for f in files
        if length(f) < 4 || f[end-3:end] != ".mat"
            continue
        end
        @printf(" %-24s", f)
        A,invar,outvar = readmat(string(testdir, f), "A", "invar", "outvar")
        invar = invar[:]
        outvar = outvar[:]
        m,n = size(A)
        A1 = [A sparse(1.0I, m, m)]
        basis = Vector{Int64}(undef,m)
        basis[:] = collect(1:m) .+ n # slack basis
        res,nfactor,nforrest,nperm = test_update(A1, basis, invar, outvar)
        @printf("%.2e      %5d factor, %6d forrest, %6d perm\n",
                res, nfactor, nforrest, nperm)
    end
    nothing
end

function test_update(A::SparseMatrixCSC{Float64,Int64}, basis::Vector{Int64}, invar::Vector{Int64},
                     outvar::Vector{Int64})
    m,n = size(A)
    niter = length(invar)
    map2basis = zeros(Int64, n)
    map2basis[basis] = 1:m
    blu = basiclu_object(m)
    param = getparam(blu)
    #param.relPivotTol = 0.5
    setparam!(blu, param)
    factorize(blu, A[:,basis])
    res = 0.0
    for k = 1:niter
        if invar[k] == outvar[k] # bound flip
            continue
        end
        j = invar[k]
        p = map2basis[outvar[k]]
        @assert p > 0
        lhs = solve_for_update(blu, A[:,j], getsol=true)
        piv = lhs[p]
        solve_for_update(blu, p)
        piverr = update(blu, piv)
        basis[p] = j
        map2basis[j] = p
        map2basis[outvar[k]] = 0
        cost = getinfo(blu, :updateCost)
        if piverr > 1e-10 || cost > 1.0
            factorize(blu, A[:,basis])
        end
        if k%100 == 0
            B = A[:,basis]
            rhs = SparseVector{Float64,Int64}(m, [1;m], [1.0;1.0])
            lhs = solve(blu, rhs, 'N')
            res = max(res, norm(B*lhs-rhs, Inf))
            lhs = solve(blu, rhs, 'T')
            res = max(res, norm(B'*lhs-rhs, Inf))
        end
    end
    info = getinfo(blu)
    return res, info.nFactorize, info.nForrestTotal, info.nUpdateTotal - info.nForrestTotal
end

"""
    test_maxvolume(testdir::String)

For all `*.mat` files in `testdir` read matrix `A` and call
[`maxvolbasis`](@ref) to construct a basis.
"""
function test_maxvolume(testdir::String)
    files = readdir(testdir)
    for f in files
        if length(f) < 4 || f[end-3:end] != ".mat"
            continue
        end
        @printf("%s\n", f)
        A, = readmat(string(testdir, f), "A")
        m, n = size(A)
        basis, obj = basiclu.maxvolbasis(A, maxpass=5)
        AI = [A I]
        B = AI[:, basis]
        # refactorize to check that basis is nonsingular
        setparam!(obj, :relPivotTol, 0.3)
        factorize(obj, B)
        info = basiclu.getinfo(obj)
        nslack = sum(basis.>n)
        @printf("%d slack, condestL = %.2e, condestU = %.2e\n\n",
                nslack, info.condestL, info.condestU)
    end
    nothing
end