Compute the ordinary Procrustes sum of squares of two matrices in Julialang. All credits to procOPA {shapes}

proc_OSS.jl

#This is just a quick translation of an ordinary Procrustes sum of squares method, as done in the R function procOPA.
#All credits to Dryden, I. L. (2014). shapes package. R Foundation for Statistical Computing, Vienna, Austria. Contributed package. Version 1-1.10. URL http://www.R-project.org

#The function are, for now, as they are. I'll be working on them in the future...

#Compute a matrix with ones on the main diagonal
#and -1/n elsewhere
function scaled_ones(n)
    eye(n) - (1/n) * ones(n,n)
end

#Scale the matrix A
function scale_matrix(A)
    return scaled_ones(size(A,1)) * A
end

function rsq(A,B)
    tbatab = transpose(B) * A * transpose(A) * B
    return sum(sqrt(abs(eigvals(tbatab))))
end

function proc_reflect(A,B)
    Abar = scale_matrix(A)
    Bbar = scale_matrix(B)
    return rsq(Abar,Bbar) / sum(diag(transpose(Bbar) * Bbar))
end

function proc_rotate_reflect(A,B)
    X = transpose(scale_matrix(A)) * scale_matrix(B)
    Xsvd = svdfact!(X)
    return transpose(Xsvd.Vt) * transpose(Xsvd.U)
end

function proc_OSS(A,B)
    R = proc_rotate_reflect(A,B)
    s = proc_reflect(A,B)
    Ahat = scale_matrix(A)
    Bhat = scale_matrix(B) * R * s
    resids = Ahat - Bhat
    return sum(diag(transpose(resids) * resids))
end

#P.S. this is a clone of an anonymous gist I mistakenly published. I'll clean up stuff as soon as github's fellows allow me to do it.