Rank calculus for JuliaLang/julia#4774

julia.jl

include("rankstability.jl")

info("Julia semantics")

axiomsused = Dict()
@axiom matmul Mat * Mat --> Mat
@axiom mattrans transpose(Mat) --> Mat
@axiom matdiv Mat / Mat --> Mat
@axiom matvec Mat * Vec --> Vec

@axiom vecmat Vec * Mat --> Mat
@axiom vectrans transpose(Vec) --> Mat

x = Vec()
y = Vec()
A = Mat()

@tryout inner(x, y)
@tryout outer(x, y)
@tryout bilinear(x, A, y)
@tryout rayleigh(A, x)
@show x'', x'' == x
@show (A*x)', x'*A', (A*x)' == (x'*A')

println("Axioms used:")
println(join(sort(collect(keys(axiomsused))), "\n"))
println()

mathematica.jl

include("rankstability.jl")
info("Mathematica semantics")

#Mathematica's Dot and Transpose functions

@axiom matmul Mat * Mat --> Mat
@axiom mattrans transpose(Mat) --> Mat
@axiom matdiv Mat / Mat --> Mat
@axiom scaldiv Scal / Scal --> Scal

@axiom vecvec Vec * Vec --> Scal
@axiom matvec Mat * Vec --> Vec
@axiom vecmat Vec * Mat --> Vec
#@axiom vectrans Transpose[{a, b, c}] is an error

x = Vec()
y = Vec()
A = Mat()

@tryout inner(x, y)
@tryout outer(x, y)
@tryout bilinear(x, A, y)
@tryout rayleigh(A, x)
@show x'', x'' == x
@show (A*x)', x'*A', (A*x)' == (x'*A')

println("Axioms used:")
println(join(sort(collect(keys(axiomsused))), "\n"))
println()

matlab.jl

include("rankstability.jl")
info("MATLAB semantics")

@axiom matmul Mat * Mat --> Mat
@axiom mattrans transpose(Mat) --> Mat
@axiom matdiv Mat / Mat --> Mat

x = Mat()
y = Mat()
A = Mat()

@tryout inner(x, y)
@tryout outer(x, y)
@tryout bilinear(x, A, y)
@tryout rayleigh(A, x)
@show x'', x'' == x
@show (A*x)', x'*A', (A*x)' == (x'*A')

println("Axioms used:")
println(join(sort(collect(keys(axiomsused))), "\n"))
println()

numpy.jl

include("rankstability.jl")
info("NumPy semantics")

#In NumPy,
# * is numpy.dot
# ' is (numpy.matrix(...).transpose()
#
# OLD numpy had numpy.matrix which did _not_ compose with numpy.ndarrays, so
# numpy.matrix does not interact with any true vectors.
#
# Quote from PEP 465:
# The result is a strong consensus among both numpy developers and developers
# of downstream packages that numpy.matrix should essentially never be used,
# because of the problems caused by having conflicting duck types for arrays.

@axiom matmul Mat * Mat --> Mat
@axiom mattrans transpose(Mat) --> Mat
@axiom matvec Mat * Vec --> Vec

@axiom vecmat Vec * Mat --> Vec
@axiom vecvec Vec * Vec --> Scal
@axiom vectrans transpose(Vec) --> Vec

@axiom scaldiv Scal / Scal --> Scal

x = Vec()
y = Vec()
A = Mat()

@tryout inner(x, y)
@tryout outer(x, y)
@tryout bilinear(x, A, y)
@tryout rayleigh(A, x)
@show x'', x'' == x
@show (A*x)', x'*A', (A*x)' == (x'*A')

println("Axioms used:")
println(join(sort(collect(keys(axiomsused))), "\n"))
println()

py35.jl

include("rankstability.jl")
info("Python 3.5/NumPy semantics")

#In NumPy,
# * is Python's @ operator for matmul
# ' is (numpy.matrix(...).transpose()

#Python's @ for matmul is discussed in PEP 465
#http://legacy.python.org/dev/peps/pep-0465/
#
#Quote:
#
#1d vector inputs are promoted to 2d by prepending or appending a '1' to the
#shape, the operation is performed, and then the added dimension is removed
#from the output. The 1 is always added on the "outside" of the shape:
#prepended for left arguments, and appended for right arguments. The result is
#that matrix @ vector and vector @ matrix are both legal (assuming compatible
#shapes), and both return 1d vectors; vector @ vector returns a scalar.

@axiom matmul Mat * Mat --> Mat
@axiom matvec Mat * Vec --> Vec
#@axiom matdiv Mat / Mat --> Mat

@axiom vecmat Vec * Mat --> Vec
@axiom vecvec Vec * Vec --> Scal
@axiom vectrans transpose(Vec) --> Vec

@axiom scaldiv Scal / Scal --> Scal

x = Vec()
y = Vec()
A = Mat()

@tryout inner(x, y)
@tryout outer(x, y)
@tryout bilinear(x, A, y)
@tryout rayleigh(A, x)
@show x'', x'' == x
@show (A*x)', x'*A', (A*x)' == (x'*A')

println("Axioms used:")
println(join(sort(collect(keys(axiomsused))), "\n"))
println()

r.jl

include("rankstability.jl")
info("R semantics")

@axiom matvec Mat * Vec --> Mat
@axiom matdiv Mat / Mat --> Mat
@axiom mattrans transpose(Mat) --> Mat

@axiom vecmat Vec * Mat --> Mat
@axiom vectrans transpose(Vec) --> Mat

x = Vec()
y = Vec()
A = Mat()

@tryout inner(x, y)
@tryout outer(x, y)
@tryout bilinear(x, A, y)
@tryout rayleigh(A, x)
@show x'', x'' == x
@show (A*x)', x'*A', (A*x)' == (x'*A')

println("Axioms used:")
println(join(sort(collect(keys(axiomsused))), "\n"))
println()

rankstability.jl

immutable Qty{N} end

typealias Scal Qty{0}
typealias Vec Qty{1}
typealias Mat Qty{2}

######################
# Fundamental axioms #

import Base: *, /, transpose
axiomsused = Dict()

macro axiom(name, expr)
    #expr looks like :(Mat * Mat --> Mat)
    @assert expr.head == :(-->)
    @assert expr.args[1].head == :call

    funcname = expr.args[1].args[1]
    nargs = length(expr.args[1].args) - 1
    axiomname = string(name)
    axiomout = expr.args[end]

    if nargs == 1
        output = quote
            global $funcname
            function $funcname(x::$(expr.args[1].args[2]))
                axiomsused[$axiomname] = get(axiomsused, $axiomname, 0) + 1
                $axiomout()
            end
        end
    elseif nargs == 2
        output = quote
            global $funcname
            function $funcname(x::$(expr.args[1].args[2]), y::$(expr.args[1].args[3]))
                axiomsused[$axiomname] = get(axiomsused, $axiomname, 0) + 1
                $axiomout()
            end
        end
    else
        error("Not implemented")
    end

    argtypes = expr.args[1].args[2:end]
    sigs = join(map(string, argtypes), ", ")
    info("Defining $axiomname: $funcname($sigs) --> $axiomout")
    output
end

@axiom matmul Mat * Mat --> Mat
@assert Mat() * Mat() == Mat()
######################
# Derived operations #

inner(x, y) = x'y
outer(x, y) = x*y'
function bilinear(x, A, y)
    z = (x'*A)*y
    if z == x'*(A*y)
        return z
    else
        error("Not associative")
    end
end
rayleigh(A, x) = x'A*x/x'x

macro tryout(expr)
    output = quote
        try
            println($(esc(expr)), " --> ", typeof($expr))
        catch exc
            println($(esc(expr)), " --> ", exc)#typeof($expr))
            #println(_EXPR, " is a $exc")
        end
    end
    #Now replace _EXPR with the actual expression
    output.args[end].args[1].args[2].args[2] =
    output.args[end].args[end].args[end].args[2] = string(expr)

    output
end