theochem / procrustes

Python library for finding the optimal transformation(s) that makes two matrices as close as possible to each other.
https://procrustes.qcdevs.org/
GNU General Public License v3.0
109 stars 20 forks source link

Non square data in softassign function - cannot reshape array #120

Closed yg42 closed 2 years ago

yg42 commented 2 years ago

Dear all, I do not know how to get the permutation matrix (like what I would get from the softassign function). I tried the following code, that works in the case of of a square matrix a (d=n), but not for a non square matrix. The error I get is

File "/home/yann/.virtualenvs/work/lib/python3.8/site-packages/procrustes/softassign.py", line 201, in softassign
    c_tensor = array_c.reshape(row_num, row_num, row_num, row_num)

ValueError: cannot reshape array of size 400 into shape (10,10,10,10)
import numpy as np
from procrustes import *

from scipy.stats import ortho_group

# random input 10x2 matrix A
n=10
d=2
a = np.random.rand(n, d)

# random orthogonal 2x2 matrix T
t = ortho_group.rvs(d)

# target matrix B (which is a shifted AxT)
b = np.dot(a, t) + np.random.rand(1, d)

# orthogonal Procrustes analysis with translation
result = orthogonal(a, b, scale=True, translate=True)

# display Procrustes results
# error (expected to be zero)
print(result.error)

res = softassign(a, b)
print(res['t'])

How could I get the permutation matrix directly from the orthogonal function ? Is this a bug ?

thank you very much for your time best regards

yann

PaulWAyers commented 2 years ago

Softassign is a developmental feature (in a package under development itself) so it's far from bulletproof.

The (current) softassign in procrustes is intended to only support one specific use case, as an implementation of the 2-sided permutation Procrustes problem with one transformation. Thus the input matrices must be square, and we need to (and will) update our documentation to indicate that the a and b matrices must be square. To use this, you could pad the matrices with zeros (and then remove the default option to eliminate zero rows/columns from softassign).

At a scientific level, the easiest way to find a permutation matrix from a orthogonal matrix is to find the closest permutation matrix to a given rotation matrix. This is just a one-sided permutation matrix with a = identity and b = target_rotation_matrix.

FarnazH commented 2 years ago

The following code will work with the current implementation, but the API of softassign is changing to become consistent with other procrustes functions:

# add zero columns to make A and B square
n, m = a.shape
square_a = np.hstack((a, np.zeros((n, n - m))))

n, m = b.shape
square_b = np.hstack((b, np.zeros((n, n - m))))

# run softassign without removing the zero columns
res = softassign(square_a, square_b, remove_zero_col=False, remove_zero_row=False)
print("error = ", res.error)
print("permutation matrix = ")
print(res.t)
yg42 commented 2 years ago

Thank you very much for this answer. We will test this solution in our problem.

FanwangM commented 2 years ago

Thank you for your interest in our package. I will close this issue for now. Please feel free to reopen it if needed. @yg42