Merge pull request #24 from LCAV/feature/srls_improvement

Algorithm improvements

.travis.yml

@@ -5,9 +5,10 @@ python:
 
 install:
   - pip install -r requirements.txt
+  - pip install -r optional_requirements.txt
 
 script:
-  - pytest
+  - pytest 
 
 notifications:
   email: [email protected]

CHANGELOG.md

@@ -1,6 +1,15 @@
 # Changelog
 
-## Unreleased 
+## Unreleased
+
+## [0.0.4] - 2019-09-12
+### Changed
+
+- Added more detailed documentation. 
+- More flexibility in passing anchor coordinates: Before, one had to add n zero-vectors to the 
+beginning of the coordinates matrix, one for each point to be localized. Now, one can pass
+the anchor coordinates only, and the zero-vectors will be added automatically where needed.
+- Cleaned up requirements.
 
 ## [0.0.3] - 2018-07-16
 ### Changed

CONTRIBUTE.rst

@@ -0,0 +1,24 @@
+How to contribute
+=================
+
+Release new version
+-------------------
+
+1. Change the version number in 
+- setup.py
+- docs/source/conf.py 
+
+2. Add a new section to CHANGELOG.md
+
+3. Run the following commands
+
+  python3 setup.py sdist bdist_wheel
+  python3 -m twine upload --repository-url https://test.pypi.org/legacy/ dist/*
+  python3 -m pip install --index-url https://test.pypi.org/simple/ --no-deps example-pkg-your-username
+  python3 -c "import pylocus.lateration" # test that it worked.
+
+4. If all looks ok, then run
+
+   python3 -m twine upload dist/*
+
+5. If above does not work, make sure that the information in ~/.pypirc is up to date.

README.rst

@@ -1,5 +1,5 @@
-pylocus 
-=======
+Welcome to pylocus 
+==================
 
 .. image:: https://travis-ci.org/LCAV/pylocus.svg?branch=master
     :target: https://travis-ci.org/LCAV/pylocus
@@ -8,7 +8,7 @@ Python Localization Package
 ---------------------------
 
 
-MDS, lateration and other algorithms, used for localization using distances and/or angles.
+This package contains Multidimensional Scaling, lateration, and other algorithms useful for localization using distances and/or angles.
 
 Local install
 *************
@@ -51,7 +51,7 @@ Depending on which parts of the project you are using, you might need to install
 
 .. code-block:: bash
 
-  pip install -r requirements.txt
+  pip install -r optional_requirements.txt
 
 Documentation
 *************

docs/source/conf.py

@@ -94,7 +94,7 @@
 # The short X.Y version.
 version = u'0.0'
 # The full version, including alpha/beta/rc tags.
-release = u'0.0.3'
+release = u'0.0.4'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

optional_requirements.txt

@@ -0,0 +1,2 @@
+cvxpy >= 1.0.6
+cvxopt >= 1.1.9

pylocus/__init__.py

@@ -1,7 +1 @@
-from .settings import *
-from .algorithms import *
-from .basics import *
-from .plots_cti import *
-from .point_set import *
-
-__version__ = "0.0.3"
+### This file exists only to tell python that this folder is a package.

pylocus/algorithms.py

@@ -58,14 +58,34 @@ def classical_mds(D):
     return MDS(D, 1, 'geometric')
 
 
+def complete_points(all_points, N):
+    ''' add zero-rows to top of all_points to reach total of N. 
+    :param all_points: m x d array, where m <= N 
+    :param N: final dimension of all_points.
+
+    :return: number of added lines (n), new array all_points of size Nxd
+    
+    '''
+    m = all_points.shape[0]
+    d = all_points.shape[1]
+    n = 1 # number of points to localize
+    if m < N:
+        n = N-m
+        all_points = np.r_[np.zeros((n, d)), all_points]
+        assert all_points.shape == (N, d)
+    elif m > N:
+        raise ValueError("Cannot have more anchor points than edm entries.")
+    return n, all_points
+
+
 def procrustes(anchors, X, scale=True, print_out=False):
     """ Fit X to anchors by applying optimal translation, rotation and reflection.
 
     Given m >= d anchor nodes (anchors in R^(m x d)), return transformation
     of coordinates X (output of EDM algorithm) optimally matching anchors in least squares sense.
 
     :param anchors: Matrix of shape m x d, where m is number of anchors, d is dimension of setup.
-    :param X: Matrix of shape N x d, where the last m points will be used to find fit with the anchors. 
+    :param X: Matrix of shape N x d, of which the last m rows will be used to find fit with the anchors. 
     :param scale: set to True if the point set should be scaled to match the anchors.
     
     :return: the transformed vector X, the rotation matrix, translation vector, and scaling factor.
@@ -74,6 +94,8 @@ def centralize(X):
         n = X.shape[0]
         ones = np.ones((n, 1))
         return X - np.multiply(1 / n * np.dot(ones.T, X), ones)
+
+    assert X.shape[1] == anchors.shape[1], 'Anchors and X must be of shape (mxd) and (Nxd), respectively.'
     m = anchors.shape[0]
     N, d = X.shape
     assert m >= d, 'Have to give at least d anchor nodes.'
@@ -86,10 +108,11 @@ def centralize(X):
     sigmaxy = 1 / m * np.dot((anchors - muy).T, X_m - mux)
     try:
         U, D, VT = np.linalg.svd(sigmaxy)
-    except np.LinAlgError:
+    except:
         print('strange things are happening...')
         print(sigmaxy)
         print(np.linalg.matrix_rank(sigmaxy))
+        raise
     #this doesn't work and doesn't seem to be necessary! (why?)
     #  S = np.eye(D.shape[0])
     #  if (np.linalg.det(U)*np.linalg.det(VT.T) < 0):
@@ -109,6 +132,7 @@ def centralize(X):
     R = np.dot(U, VT)
     t = muy.T - c * np.dot(R, mux.T)
     X_transformed = (c * np.dot(R, (X - mux).T) + muy.T).T
+    assert np.allclose(X_transformed.shape, X.shape)
     return X_transformed, R, t, c
 
 
@@ -167,13 +191,24 @@ def reconstruct_cdm(dm, absolute_angles, all_points, W=None):
     return Y
 
 
-def reconstruct_mds(edm, all_points, completion='optspace', mask=None, method='geometric', print_out=False, n=1):
+def reconstruct_mds(edm, all_points, completion='optspace', mask=None, method='geometric', print_out=False):
     """ Reconstruct point set using MDS and matrix completion algorithms.
+
+    :param edm: Euclidean distance matrix, NxN. 
+    :param all_points: Mxd vector of anchor coordinates. if M < N, M-N 0-row-vectors will be added to the beginning of the array. If M=N and n is not specified, we assume that the first row corresponds to the point to be localized. 
+    :param completion: can be 'optspace' or 'alternate'. Algorithm to use for matrix completion. See pylocus.edm_completion for details.
+    :param mask: Optional mask of missing entries.  
+    :param method: method to be used for MDS algorithm. See method "MDS" from pylocus.mds module for details. 
+
     """
+
     from .point_set import dm_from_edm
     from .mds import MDS
-    N = all_points.shape[0]
     d = all_points.shape[1]
+    N = edm.shape[0]
+
+    n, all_points = complete_points(all_points, N)
+
     if mask is not None:
         edm_missing = np.multiply(edm, mask)
         if completion == 'optspace':
@@ -191,8 +226,8 @@ def reconstruct_mds(edm, all_points, completion='optspace', mask=None, method='g
             print('{}: relative error:{}'.format(completion, err))
         edm = edm_complete
     Xhat = MDS(edm, d, method, False).T
+    assert (~np.isnan(Xhat)).all()
     Y, R, t, c = procrustes(all_points[n:], Xhat, True)
-    #Y, R, t, c = procrustes(all_points, Xhat, True)
     return Y
 
 
@@ -206,11 +241,19 @@ def reconstruct_sdp(edm, all_points, W=None, print_out=False, lamda=1000, **kwar
     return Xhat, edm_complete
 
 
-def reconstruct_srls(edm, all_points, W=None, print_out=False, n=1, rescale=False,
-                     z=None):
+def reconstruct_srls(edm, all_points, W=None, print_out=False, rescale=False, z=None):
     """ Reconstruct point set using S(quared)R(ange)L(east)S(quares) method.
+
+    See reconstruct_mds for edm and all_points parameters.
+
+    :param W: optional weights matrix, same dimension as edm. 
+    :param rescale, z: optional parameters for SRLS. See SRLS documentation for explanation (module pylocus.lateration)
     """
     from .lateration import SRLS, get_lateration_parameters
+
+    N = edm.shape[0]
+    n, all_points = complete_points(all_points, N)
+
     Y = all_points.copy()
     indices = range(n)
     for index in indices:
@@ -222,7 +265,17 @@ def reconstruct_srls(edm, all_points, W=None, print_out=False, n=1, rescale=Fals
             print('w', w)
             print('r2', r2)
 
-        srls = SRLS(anchors, w, r2, rescale, z, print_out)
+        try:
+            srls = SRLS(anchors, w, r2, rescale, z, print_out)
+            if z is not None:
+                assert srls[2] == z
+
+        except Exception as e:
+            print(e)
+            print("Something went wrong; probably bad geometry. (All anchors in the same plane, two distances are exactly the same, etc.)")
+            print("anchors, w, r2, z:", anchors, w, r2, z)
+            return None
+
         if rescale:
             srls = srls[0]  # second element of output is the scale
         Y[index, :] = srls

pylocus/basics.py

@@ -82,6 +82,10 @@ def eigendecomp(G, d):
     assert (lamda_sorted[
             :d] > -1e-10).all(), "{} not all positive!".format(lamda_sorted[:d])
 
+    # Set the small negative values of 
+    # lamda to zero.
+    lamda_sorted[lamda_sorted < 0] = 0
+
     u = u[:, indices]
     factor = np.empty((N,), dtype=lamda.dtype)
     np.sqrt(lamda_sorted[:d], out=factor[0:d])

pylocus/edm_completion.py

@@ -1,7 +1,11 @@
 #!/usr/bin/env python
 # module EDM_COMPLETION
 import numpy as np
-from cvxpy import *
+
+try:
+    from cvxpy import *
+except:
+    print("WARNING from pylocs.edm_completion module: Failed to load cvxpy. This might lead to errors later on.")
 
 from pylocus.basics import get_edm