theochem / procrustes Goto Github PK

This is a long-term goal and it's mentioned in theochem/Selector#55.

Do you mean having 1D (SMILES), 2D or 3D structures as input for Procrustes, and then Procrustes would compute the similarities between the molecules?
I am thinking of having such functionality as a tools module. What do you think? @PaulWAyers

Turn the code blocks given in tutorials into Python scripts (and then reference them in the documentation), so it can be easily downloaded and executed. It is also important to check that they run properly as the code changes.

Example: http://theochem.github.io/horton/2.1.1/user_postproc_aim.html#python-interface-to-the-partitioning-code

In permutation Procrustes, the error is computed by explicit multiplications like (A x P) or (P1 x A x P2). This is extremely inefficient because one is only rearranging rows/columns of A (an O(n^2) operation) but one is using matrix multiplication (which is O(n^3)) to do this. This can be done efficiently by:

1. For each permutation matrix, construct the corresponding permutation of indices; then constructed the matrix with permuted rows/columns with fancy indexing.
2. For each permutation matrix and the original matrix, construct a sparse matrix; use sparse matrix multiplication. The sparse matrix representation for A should be constructed only once or otherwise you'll give up much of the gains due to overhead.

(a) would be preferred in my view. It amounts to a more efficient way to evaluate the error function for k-opt in its application to Procrustes.

Can you check why the testing is failing for https://github.com/theochem/procrustes/runs/8172507556?check_suite_focus=true? I have tried to rerun the testing and it didn't help. I don't know if it's Scipy or Numpy version related or operating system related. The tests work before we merge the PR.

Thank you. @banrovegrie

@fwmeng88 should we add a tutorial on how to compute the distance between molecular structures/isomers, along the lines of what is done in the project with @rayhe88 and @Nick6618 .

This is a reminder to myself that the documentation, jupyter notebooks, and docstrings should be fixed in corresponding to the changes of function returns.

We should get the https://procrustes.qcdevs.org/ web site working.

This positive semi-definite Procrustes problem is trying to

More details can be found at Woodgate, K. G. (1993, December). A new algorithm for the positive semi-definite Procrustes problem. In Proceedings of 32nd IEEE Conference on Decision and Control (pp. 3596-3601). IEEE..

This will be on the to-do list of our package, but not for this stage.

I am on Python 3.8.0. When I try to install procrustes from pip, I get this error:

ERROR: Could not find a version that satisfies the requirement qc-procrustes (from versions: none)
ERROR: No matching distribution found for qc-procrustes

Also, I can't clone the repo, since I get permission denied:

Cloning into 'procrustes'...
[email protected]: Permission denied (publickey).
fatal: Could not read from remote repository.

Please make sure you have the correct access rights
and the repository exists.

This is a reminder to myself about updating the documentation.

This is on Python 3.9. pip install qc-procrustes results in this error:

ERROR: Could not find a version that satisfies the requirement qc-procrustes (from versions: none)
ERROR: No matching distribution found for qc-procrustes

I tried doing protein alignment with IOData as you have built for the manuscript and the code snippet seems not working.

import numpy as np

from iodata import load_one
from procrustes import rotational

# load PDB
pdb = load_one("2hhb.pdb")

# get coordinates of C_alpha atoms in chains A & C (in angstrom)
chainid = pdb.extra['chainid']
attypes = pdb.atffparams['attypes']
ca_a = pdb.atcoords[(chainid == 'A') & (attypes == 'CA')] * 0.529177
ca_c = pdb.atcoords[(chainid == 'C') & (attypes == 'CA')] * 0.529177

print("RMSD of intial coordinates:")   # output: 39.4
print(np.sqrt(np.mean(np.sum((ca_a - ca_c)**2, axis=1))))

# rotational Procrustes analysis
result = rotational(ca_a, ca_c, translate=True)

# compute transformed (translated & rotated) coordinates of chain A
ca_at = np.dot(result.new_a, result.array_u)

print("RMSD of transformed coordinates:")   # output: 0.23
print(np.sqrt(np.mean(np.sum((ca_at - result.new_b)**2, axis=1))))

I guess the problem is that there is no attribute called chainid after checking the latest source code.
Do you know which part I get things wrong? Thank you.
@FarnazH

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

An interesting paper proposed to solve the rotational Procrustes problem without SVD, J. Geom. Symmetry Phys., 53, 21–53, 2019. We can support this in the next major release.

Now the scaling factor and translation vector are calculated in the utilities but they are not returned in the answer (result = orthogonal(a, b, scale=True, translate=True)). For 3D-graphics calculations, these would be handy to scale and translate any new points found from the distance matrix (I get preliminary coordinates from MDS) back into the real coordinates of b.

Move the functions for kopt_heuristic to a new module called kopt. They are currently in utils module.

During testing with github action, the follow test failed "test_generic_rectangular_translate_scale" for test_generic.py" with m = 829, n = 878 on python 3.7.

The error is
def _raise_linalgerror_svd_nonconvergence(err, flag): raise LinAlgError("SVD did not converge") E numpy.linalg.LinAlgError: SVD did not converge.

For more info on the github action, see https://github.com/Ali-Tehrani/procrustes/runs/2058961187.

However, running this test on my own machine it passes. So it seems the error is hardware related. Google searching this returns a 2019 SciPy issue [1]. I would recommend to decrease the matrix size so that it isn't hardware/MKL/OpenBlas dependent.

[1] - scipy/scipy#10032

I remember that at some point we talked about changing the argument names to array_a and array_b to a and b. I still think it is good to be consistent with numpy and scipy and use shorter argument names.
See: https://github.com/theochem/procrustes/blob/master/procrustes/orthogonal.py#L37
What do you think @PaulWAyers?

I'm happy with putting references centrally. I don't worry so much about keeping them up-to-date, or having long documentation (especially if the references are just at the very end anyway), since most of the references we should be including are just "the classics." But I know that because our documentation is embedded in the code, it can be annoying to have a lot of documentation before the code even starts.

I guess the question, then, is whether it is good to have the specific reference where a reader can get more information about the underlying algorithm/method and its justification. I would tend towards including references in the code I guess (where one might wonder what's going on) but as long as there is a source for the references, it is OK I think.

Originally posted by @PaulWAyers in #68 (comment)

The pinv2 in Scipy has been deprecated since 1.7.0, and it has been removed for 1.9.0 and above, scipy/scipy#16245. Originally, the (Moore-Penrose) pseudo-inverse is computed by SVD including all 'large' singular values (using scipy.linalg.pinv2) and scipy.linalg.pinv obtains the (Moore-Penrose) pseudo-inverse by using least-squares solver. The least-squares implementation is less efficient, but more robust, than the SVD implementation (

So, for newer versions of Scipy >= 1.9.0, it would be a problem to use pinv2. Therefore, I think the easiest solution is to stick to pinv and remove pinv2. Another option is to provide different solutions for different versions of Scipy greater than 1.9.0 or less than 1.9.0. What do you prefer? @PaulWAyers @Ali-Tehrani

Hello,
Since scipy 1.9.0 , the function pinv2 has been subsumed in pinv, is it possible to remove it's use in generic.py

P.S
Attaching link to the release notes https://docs.scipy.org/doc/scipy/release.1.9.0.html?highlight=pinv2

Actually, why make a permutation matrix (mostly zeros and ones) and multiply by A (cost of O(n^3)). You can directly compute A with permuted rows/columns, then compute the error in (permuted) A. This has cost O(n^2) (computing the error) instead of O(n^3). IF the error is reduced, you can return the permutation, or even the permutation matrix (which can be constructed at the very end).

Originally posted by @PaulWAyers in #62 (comment)

The psdp_woodgate (#170) has set up a good starting point for implementing PSDP algorithms. Thanks! @banrovegrie

We are missing a few initial setting up such as scaling, translation, etc. You can add this code blocks like

 # check inputs
    new_a, new_b = setup_input_arrays(
        a,
        b,
        unpad_col,
        unpad_row,
        pad,
        translate,
        scale,
        check_finite,
        weight,
    )
    if new_a.shape != new_b.shape:
        raise ValueError(
            f"Shape of A and B does not match: {new_a.shape} != {new_b.shape} "
            "Check pad, unpad_col, and unpad_row arguments."
        )

And you will need to change the finial output accordingly, such as new_a=new_a and new_b=new_b. Sample codes can be found at

Given that you want the first PR to be merged in a quicker manner, I have merged #170. You can fix this in a new PR together with adding testing codes or you can just create a small PR for this.

This is a reminder for myself that a few things need to be fixed.

• Update the header
• Update contact information using QC-Dev community
• Update the footnote information in conf.py
• Update the equations and make it consistent with our manuscript.
• Update installation guide
• Update quick start example
• Add API rst file for kopt module
• fix examples, following . This related to #26
• Update docstrings if necessary: equations, argument names, and default parameters. This relates to #25.
• Add RST file style checker

In the docs, the transform_mode argument for permutation_2sided should be one of single_undirected, single_directed, double. In the code, it can only be single or double, and then the directed/undirected mode is chosen based on testing the suitability of the matrix.

The Procrustes method can be solved where the transformation is a Toeplitz matrix
http://dx.doi.org/10.1155/2013/696019
I don't know any applications of this, so I consider it low priority. But it does look straightforward.

• Jupyter notebooks for all the tutorial examples
• Avoid using private functions
• Use https://mybinder.org/ to host the jupyter notebooks

Currently, we allow weights for each dimension, and then use a (diagonal) weight matrix. In some applications in quantum mechanics, it would be useful to have a positive semidefinite weight matrix (which would be symmetric in all cases I can imagine). The idea would be to let the user pass a vector (then construct a diagonal weight matrix) or a matrix (then use the matrix).

In the cases I'm thinking about, the weight matrix would be the atomic-orbital overlap matrix or the one-electron reduced density matrix.

I a bit confused with the generalized Procrustes method. I understand that every matrix has a right-hand-side transformation, but does it needs to be an orthogonal transformation (as implemented)? The documentation does not say anything about it and just used T. See: https://github.com/theochem/procrustes/blob/master/procrustes/generalized.py#L87

I noticed that the binder doesn't work.

The weigh argument is added to all flavors of Procrustes but it does not show up in the mathematical formulation of the problem, also the documentation is not clear at all on how it would be used.

When performing generalized Procrustes analysis, it is very likely to have missing values in the configuration matrices. Casper J. Albers and John C. Gower proposed a strategy to deal with such a situation, _Albers, C. J., & Gower, J. C. (2010). A general approach to handling missing values in Procrustes analysis. Advances in Data Analysis and Classification, 4(4), 223-237.. A very nice feature of their algorithm is they can also deal with centring and/or standardisation.

This issue relates to #76.

Make sure that similar to IOData, we can install Procrustes through:

pip install qc-procrustes

We need to add some descriptions for this package.

In chapter 7.5.2 of Procrustes Problems by J. C. Gower and G. G. Dijksterhuis an algorithm for 2-sided orthogonal Procrustes is presented. This generalizes the current treatment, which requires that the matrices being matched be symmetric and square. This is also explained in section 3.3 of Pythagoras Papadimitriou's Ph.D. thesis at the Univ. of Manchester.
papad93.pdf

Right now, the object is to find an orthogonal matrices, Q1 and Q2, such that || Q1 A Q2 - B || is minimized, where A and B are symmetric square matrices and have the same shape (i.e., both are mxm matrices). However, this reference shows that this can work even as long as A and B have the same shape; they need not be square or symmetric.

This could be an embellishment on the existing method (it contains the old method as a special case) or it could be implemented as a separate method. I'd tend to do the latter, and leave the easy case as a routine that can be used where applicable.

@fanwang I couldn't regenerate the protein alignment results you had presented (using https://github.com/theochem/procrustes/blob/master/docs/notebooks/Protein_Structure_Alignment.ipynb). Can you also please take a look at it?
The screenshot below shows the final RMSD I have got (which doesn't match yours) and the error when plotting the transformed coordinates.

Hi, Thanks for the project! I'm interested in solving procrustes using a more robust metric of dispersion than F norm.

I saw a couple papers online that solve e.g. smooth surrogates of l1, but wondering if you have any pointers in particular.

In our implementations, when two input matrices have different shapes (unbalanced), zero-padding is used to make them balanced. An alternative method is proposed in a recent SIAM paper,

LH Zhang, WH Yang, C Shen, J Ying, SIAM J. Matrix Anal. Appl., 41(3), 957–983..

In this paper, they showed that the orthogonal Procrustes can be transformed into an eigenvalue minimization problem, in which an adapted self-consistent field (SCF) iteration is used to solve the problem.

We can leave this for future implementation. Anyone interested to implement the algorithm is also welcome!

Analogous to what was done with k-opt, we should add functionality to optimize a user-specified function over the set of doubly-stochastic matrices. (Doubly stochastic matrices are nonnegative matrices that have row-sum and column-sum equal to 1.)

• scipy.opt could be used for the (local) optimization starting from a user-specified starting point. This is a classic problem for the slsqp algorithm but the trust-constr is also appropriate. Which algorithm could be a user option. Or just choose one; this is intended as a useful utility and should not be confused with a cutting-edge treatment for this problem.
• The default starting point could be constructed as (1) take a random orthogonal matrix; (2) square all the elements).
• as a user option, the function could return the permutation matrix that is closest to the doubly-stochastic matrix that is found.

This sort of problem showed up in @wilhadams work. For an objective function of the form in one-sided permutation Procrustes, it amounts to a (slow, iterative) alternative to the Hungarian algorithm for the assignment problem.

Note that in Procrustes problems of a quite general (tensor-ish) sense, it is usually important to maximize (or minimize the negative of) the cosine/overlap between the optimizee and the target, as direct minimization of the error doubles the degree of the optimization.

The original "algorithm 3" has been peer reviewed (I think; it is published in a conference proceedings and that is never a sure thing).
http://www.scielo.org.co/pdf/rcm/v55n1/0034-7426-rcm-55-01-109.pdf

However, I really like Fanwang's suggestion. The Oviedo paper (algorithm 3) performs similarly to what they call PARATAN and what the above reference calls PARTAN. The above reference seems to have good performance, and has four closely related algorithms which could be implemented simulataneously, switching between the algorithms as specified by the user by a keyword. This is explained at the top of Section 5.

Specifically, the

• projected gradient method
• fast gradient method
• parallel tangents method
• semi-analytic fast gradient method

This gives a nice alternative to the first two algorithms, which focus on a more analytic approach (more like the other methods we have considered). This is a numerical approach, and perhaps might be faster (hard to know until we test).

Originally posted by @PaulWAyers in #189 (comment)

We currently support the solution of the general (unconstrained transformation) one-sided Procrustes problem but not the two-sided unconstrained Procrustes problem. A statement, and solution, of this problem can be found in section 3.2 of the Ph.D. thesis of Pythagoras Papadimitriou (Univ. Manchester). It would be good to support this case also.
papad93.pdf

For example, instead of e_opt, I think error would be a more intuitive argument.
See: https://github.com/theochem/procrustes/blob/master/procrustes/utils.py#L562

To be completed.

From #69 & #77: Add a function variable "lapack_driver" to orthogonal, permutation, rotational, and symmetric functions with docstring lapack_driver : {"gesvd", "gesdd"}, optional used in the singular value decomposition function from SciPy. Only allowed two options, with "gesvd" being less-efficient than "gesdd" but is more robust. Default is "gesvd".