diff --git a/.travis.yml b/.travis.yml
index 9231118..8df6985 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,9 +5,10 @@ python:
 
 install:
   - pip install -r requirements.txt
+  - pip install -r optional_requirements.txt
 
 script:
-  - pytest
+  - pytest 
 
 notifications:
   email: frederike.duembgen@epfl.ch
diff --git a/CHANGELOG.md b/CHANGELOG.md
index fdd7564..65e0574 100644
--- a/CHANGELOG.md
+++ b/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
diff --git a/CONTRIBUTE.rst b/CONTRIBUTE.rst
new file mode 100644
index 0000000..5dc9a7a
--- /dev/null
+++ b/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.
diff --git a/README.rst b/README.rst
index 3ad708e..b88a991 100644
--- a/README.rst
+++ b/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
 *************
diff --git a/docs/source/conf.py b/docs/source/conf.py
index b1a8a91..7d9a8cd 100644
--- a/docs/source/conf.py
+++ b/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.
diff --git a/optional_requirements.txt b/optional_requirements.txt
new file mode 100644
index 0000000..596bcca
--- /dev/null
+++ b/optional_requirements.txt
@@ -0,0 +1,2 @@
+cvxpy >= 1.0.6
+cvxopt >= 1.1.9
diff --git a/pylocus/__init__.py b/pylocus/__init__.py
index 2fb9426..44b35bf 100644
--- a/pylocus/__init__.py
+++ b/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.
diff --git a/pylocus/algorithms.py b/pylocus/algorithms.py
index 5cdb63a..b07a73c 100644
--- a/pylocus/algorithms.py
+++ b/pylocus/algorithms.py
@@ -58,6 +58,26 @@ 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.
 
@@ -65,7 +85,7 @@ def procrustes(anchors, X, scale=True, print_out=False):
     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
diff --git a/pylocus/basics.py b/pylocus/basics.py
index cdf5e13..d201915 100644
--- a/pylocus/basics.py
+++ b/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])
diff --git a/pylocus/edm_completion.py b/pylocus/edm_completion.py
index ea73812..ebd79a2 100644
--- a/pylocus/edm_completion.py
+++ b/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
 
diff --git a/pylocus/lateration.py b/pylocus/lateration.py
index ec44866..365b620 100644
--- a/pylocus/lateration.py
+++ b/pylocus/lateration.py
@@ -3,10 +3,18 @@
 
 import numpy as np
 from scipy.linalg import eigvals, eigvalsh
-from cvxpy import *
+
+try:
+    from cvxpy import *
+except:
+    print("WARNING from pylocs.lateration module: Failed to load cvxpy. This might lead to errors later on.")
+
 
 from pylocus.basics import assert_print, assert_all_print
 
+class GeometryError(Exception):
+    pass
+
 
 def get_lateration_parameters(all_points, indices, index, edm, W=None):
     """ Get parameters relevant for lateration from full all_points, edm and W.
@@ -42,24 +50,10 @@ def get_lateration_parameters(all_points, indices, index, edm, W=None):
     return anchors, w, r2
 
 
-def SRLS(anchors, w, r2, rescale=False, z=None, print_out=False):
-    '''Squared range least squares (SRLS)
-
-    Algorithm written by A.Beck, P.Stoica in "Approximate and Exact solutions of Source Localization Problems".
-
-    :param anchors: anchor points (Nxd)
-    :param w: weights for the measurements (Nx1)
-    :param r2: squared distances from anchors to point x. (Nx1)
-    :param rescale: Optional parameter. When set to True, the algorithm will
-        also identify if there is a global scaling of the measurements.  Such a
-        situation arise for example when the measurement units of the distance is
-        unknown and different from that of the anchors locations (e.g. anchors are
-        in meters, distance in centimeters).
-    :param z: Optional parameter. Use to fix the z-coordinate of localized point.
-    :param print_out: Optional parameter, prints extra information.
+def solve_GTRS(A, b, D, f):
+    from scipy import optimize
+    from scipy.linalg import sqrtm
 
-    :return: estimated position of point x.
-    '''
     def y_hat(_lambda):
         lhs = ATA + _lambda * D
         assert A.shape[0] == b.shape[0]
@@ -70,25 +64,73 @@ def y_hat(_lambda):
         try:
             return np.linalg.solve(lhs, rhs)
         except:
-            return np.zeros((lhs.shape[1],))
+            # TODO: It was empirically found that it is essential that the default is 
+            # not zero, but a small negative value. It is not clear why this is the case
+            # from a mathematical point of view. 
+            return np.full((lhs.shape[1],), -1e-3)
 
     def phi(_lambda):
         yhat = y_hat(_lambda).reshape((-1, 1))
         sol = np.dot(yhat.T, np.dot(D, yhat)) + 2 * np.dot(f.T, yhat)
         return sol.flatten()
 
-    def phi_prime(_lambda):
-        # TODO: test this.
-        B = np.linalg.inv(ATA + _lambda * D)
-        C = A.T.dot(b) - _lambda*f
-        y_prime = -B.dot(D.dot(B.dot(C)) - f)
-        y = y_hat(_lambda)
-        return 2*y.T.dot(D).dot(y_prime) + 2*f.T.dot(y_prime)
+    ATA = np.dot(A.T, A)
 
-    from scipy import optimize
-    from scipy.linalg import sqrtm
+    try:
+        eig = np.sort(np.real(eigvalsh(a=D, b=ATA)))
+    except:
+        raise GeometryError('Degenerate configuration.')
+    if np.abs(eig[-1]) < 1e-10:
+        lower_bound = -1e3
+    else:
+        lower_bound = - 1.0 / eig[-1] 
+
+    inf = 1e5
+    xtol = 1e-12
+
+    lambda_opt = 0
+    # We will look for the zero of phi between lower_bound and inf. 
+    # Therefore, the two have to be of different signs. 
+    if (phi(lower_bound) > 0) and (phi(inf) < 0): 
+        try: 
+            # brentq is considered the best rootfinding routine. 
+            lambda_opt, r = optimize.brentq(phi, lower_bound, inf, xtol=xtol, full_output=True)
+            assert r.converged
+        except:
+            print('SRLS error: brentq did not converge even though we found an estimate for lower and upper bonud. Setting lambda to 0.')
+    else: 
+        try: 
+            lambda_opt = optimize.newton(phi, lower_bound, maxiter=10000, tol=xtol)
+            assert phi(lambda_opt) < xtol, 'did not find solution of phi(lambda)=0:={}'.format(phi(lambda_opt))
+        except AssertionError:
+            print('SRLS error: newton did not converge. Setting lambda to 0.')
+
+    yhat = y_hat(lambda_opt)
+    return yhat
+
+
+def SRLS(anchors, w, r2, rescale=False, z=None, print_out=False):
+    """ Squared range least squares (SRLS)
+
+    Algorithm written by A.Beck, P.Stoica in "Approximate and Exact solutions of Source Localization Problems".
+
+    :param anchors: anchor points (Nxd)
+    :param w: weights for the measurements (Nx1)
+    :param r2: squared distances from anchors to point x. (Nx1)
+    :param rescale: Optional parameter. When set to True, the algorithm will
+        also identify if there is a global scaling of the measurements.  Such a
+        situation arise for example when the measurement units of the distance is
+        unknown and different from that of the anchors locations (e.g. anchors are
+        in meters, distance in centimeters).
+    :param z: Optional parameter. Use to fix the z-coordinate of localized point.
+    :param print_out: Optional parameter, prints extra information.
+
+    :return: estimated position of point x.
+    """
 
     n, d = anchors.shape
+    if type(r2) == list:
+        r2 = np.array(r2).reshape((-1, 1))
     assert r2.shape[1] == 1 and r2.shape[0] == n, 'r2 has to be of shape Nx1'
     assert w.shape[1] == 1 and w.shape[0] == n, 'w has to be of shape Nx1'
     if z is not None:
@@ -96,119 +138,92 @@ def phi_prime(_lambda):
 
     if rescale and z is not None:
         raise NotImplementedError('Cannot run rescale for fixed z.')
+
     if rescale and n < d + 2:
-        raise ValueError(
-            'A minimum of d + 2 ranges are necessary for rescaled ranging.')
-    elif n < d + 1 and z is None:
-        raise ValueError(
-            'A minimum of d + 1 ranges are necessary for ranging.')
-    elif n < d:
-        raise ValueError(
-            'A minimum of d ranges are necessary for ranging.')
+        raise ValueError('A minimum of d + 2 ranges are necessary for rescaled ranging.')
+    elif z is None and n < d + 1:
+        raise ValueError('A minimum of d + 1 ranges are necessary for ranging.')
+    elif z is not None and n < d:
+        raise ValueError('A minimum of d ranges are necessary for ranging.')
+
+    if rescale: 
+        return SRLS_rescale(anchors, w, r2, print_out)
+
+    if z is not None: 
+        return SRLS_fixed_z(anchors, w, r2, z)
+
+    A = np.c_[-2 * anchors, np.ones((n, 1))]
+    b = r2 - np.power(np.linalg.norm(anchors, axis=1), 2).reshape(r2.shape)
 
     Sigma = np.diagflat(np.power(w, 0.5))
+    A = Sigma.dot(A)
+    b = Sigma.dot(b)
 
-    if rescale:
-        A = np.c_[-2 * anchors, np.ones((n, 1)), -r2]
-    else:
-        if z is None:
-            A = np.c_[-2 * anchors, np.ones((n, 1))]
-        else:
-            A = np.c_[-2 * anchors[:, :2], np.ones((n, 1))]
+    D = np.zeros((d + 1, d + 1))
+    D[:-1, :-1] = np.eye(D.shape[0]-1)
 
-        A = Sigma.dot(A)
+    f = np.c_[np.zeros((1, d)), -0.5].T
 
-    if rescale:
-        b = - np.power(np.linalg.norm(anchors, axis=1), 2).reshape(r2.shape)
-    else:
-        b = r2 - np.power(np.linalg.norm(anchors, axis=1), 2).reshape(r2.shape)
-        if z is not None:
-            b = b + 2 * anchors[:, 2].reshape((-1, 1)) * z - z**2
-        b = Sigma.dot(b)
+    yhat = solve_GTRS(A, b, D, f)
+    return yhat[:d]
 
-    ATA = np.dot(A.T, A)
 
-    if rescale:
-        D = np.zeros((d + 2, d + 2))
-        D[:d, :d] = np.eye(d)
-    else:
-        if z is None:
-            D = np.zeros((d + 1, d + 1))
-        else:
-            D = np.zeros((d, d))
-        D[:-1, :-1] = np.eye(D.shape[0]-1)
-
-    if rescale:
-        f = np.c_[np.zeros((1, d)), -0.5, 0.].T
-    elif z is None:
-        f = np.c_[np.zeros((1, d)), -0.5].T
-    else:
-        f = np.c_[np.zeros((1, 2)), -0.5].T
-
-    eig = np.sort(np.real(eigvalsh(a=D, b=ATA)))
-    if (print_out):
-        print('ATA:', ATA)
-        print('rank:', np.linalg.matrix_rank(A))
-        print('eigvals:', eigvals(ATA))
-        print('condition number:', np.linalg.cond(ATA))
-        print('generalized eigenvalues:', eig)
-
-    #eps = 0.01
-    if eig[-1] > 1e-10:
-        lower_bound = - 1.0 / eig[-1] 
-    else:
-        print('Warning: biggest eigenvalue is zero!')
-        lower_bound = -1e-5
+def SRLS_rescale(anchors, w, r2, print_out=False):
+    """ Helper routine for SRLS. """
+    n, d = anchors.shape
 
-    inf = 1e5
-    xtol = 1e-12
+    A = np.c_[-2 * anchors, np.ones((n, 1)), -r2]
+    b = - np.power(np.linalg.norm(anchors, axis=1), 2).reshape(r2.shape)
 
-    lambda_opt = 0
-    # We will look for the zero of phi between lower_bound and inf. 
-    # Therefore, the two have to be of different signs. 
-    if (phi(lower_bound) > 0) and (phi(inf) < 0): 
-        # brentq is considered the best rootfinding routine. 
-        try: 
-            lambda_opt = optimize.brentq(phi, lower_bound, inf, xtol=xtol)
-        except:
-            print('SRLS error: brentq did not converge even though we found an estimate for lower and upper bonud. Setting lambda to 0.')
-    else: 
-        try: 
-            lambda_opt = optimize.newton(phi, lower_bound, fprime=phi_prime, maxiter=1000, tol=xtol, verbose=True)
-            assert phi(lambda_opt) < xtol, 'did not find solution of phi(lambda)=0:={}'.format(phi(lambda_opt))
-        except:
-            print('SRLS error: newton did not converge. Setting lambda to 0.')
+    Sigma = np.diagflat(np.power(w, 0.5))
+    A = Sigma.dot(A)
+    b = Sigma.dot(b)
 
-    if (print_out):
-        print('phi(lower_bound)', phi(lower_bound))
-        print('phi(inf)', phi(inf))
-        print('phi(lambda_opt)', phi(lambda_opt))
-        pos_definite = ATA + lambda_opt*D
-        eig = np.sort(np.real(eigvals(pos_definite)))
-        print('should be strictly bigger than 0:', eig)
+    D = np.zeros((d + 2, d + 2))
+    D[:d, :d] = np.eye(d)
 
-    # Compute best estimate
-    yhat = y_hat(lambda_opt)
+    f = np.c_[np.zeros((1, d)), -0.5, 0.].T
 
-    if print_out and rescale:
+    yhat = solve_GTRS(A, b, D, f)
+
+    if print_out:
         print('Scaling factor :', yhat[-1])
+    return yhat[:d], yhat[-1]
 
-    if rescale:
-        return yhat[:d], yhat[-1]
-    elif z is None:
-        return yhat[:d]
-    else:
-        return np.r_[yhat[0], yhat[1], z]
+
+def SRLS_fixed_z(anchors, w, r2, z):
+    """ Helper routine for SRLS. """
+    n, d = anchors.shape
+
+    Sigma = np.diagflat(np.power(w, 0.5))
+
+    A = np.c_[-2 * anchors[:, :2], np.ones((n, 1))]
+
+    b = r2 - np.power(np.linalg.norm(anchors, axis=1), 2).reshape(r2.shape)
+    b = b + 2 * anchors[:, 2].reshape((-1, 1)) * z - z**2
+
+    A = Sigma.dot(A)
+    b = Sigma.dot(b)
+
+    ATA = np.dot(A.T, A)
+
+    D = np.zeros((d, d))
+    D[:-1, :-1] = np.eye(D.shape[0]-1)
+
+    f = np.c_[np.zeros((1, 2)), -0.5].T
+
+    yhat = solve_GTRS(A, b, D, f)
+    return np.r_[yhat[0], yhat[1], z]
 
 
 def PozyxLS(anchors, W, r2, print_out=False):
-    ''' Algorithm used by pozyx (https://www.pozyx.io/Documentation/how_does_positioning_work)
+    """ Algorithm used by pozyx (https://www.pozyx.io/Documentation/how_does_positioning_work)
 
     :param anchors: anchor points
     :param r2: squared distances from anchors to point x.
 
     :returns: estimated position of point x.
-    '''
+    """
     N = anchors.shape[0]
     anchors_term = np.sum(np.power(anchors[:-1], 2), axis=1)
     last_term = np.sum(np.power(anchors[-1], 2))
@@ -223,7 +238,7 @@ def RLS(anchors, W, r, print_out=False, grid=None, num_points=10):
 
     Algorithm written by A.Beck, P.Stoica in "Approximate and Exact solutions of Source Localization Problems".
 
-    :param anchors: anchor points
+    :param anchors: anchor point coordinates, N x d
     :param r2: squared distances from anchors to point x.
     :param grid: where to search for solution.  (min, max) where min and max  are 
     lists of d elements, d being the dimension of the setup. If not given, the 
diff --git a/pylocus/mds.py b/pylocus/mds.py
index 54325bb..9108fed 100644
--- a/pylocus/mds.py
+++ b/pylocus/mds.py
@@ -1,7 +1,10 @@
 #!/usr/bin/env python
 # module MDS
 import numpy as np
-from cvxpy import *
+try:
+    from cvxpy import *
+except:
+    print("WARNING from pylocs.mds module: Failed to load cvxpy. This might lead to errors later on.")
 
 from pylocus.basics import eigendecomp
 
@@ -20,7 +23,13 @@ def x_from_eigendecomp(factor, u, dim):
 
 
 def MDS(D, dim, method='simple', theta=False):
-    """ recover points from euclidean distance matrix using classic MDS algorithm. 
+    """ Recover points from euclidean distance matrix using classic MDS algorithm. 
+    
+    :param D: Euclidean distance matrix.
+    :param dim: dimension of points. 
+    :param method: method to use. Can be either "simple", "advanced" or "geometric". 
+    :param theta: set to True if using this function in a 1D-angle-setting. 
+
     """
     N = D.shape[0]
     if method == 'simple':
@@ -168,7 +177,3 @@ def signedMDS(cdm, W=None):
     x_est -= np.min(x_est)
 
     return x_est, A, np.linalg.pinv(A)
-
-
-if __name__ == "__main__":
-    print('nothing happens when running this module.')
diff --git a/requirements.txt b/requirements.txt
index 67489cc..51aa5e0 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,3 +1,2 @@
-cvxpy >= 1.0.6
-cvxopt
-matplotlib
+numpy >= 1.14.3
+matplotlib >= 3.0.2
diff --git a/setup.py b/setup.py
index 310f3d6..4e844c7 100644
--- a/setup.py
+++ b/setup.py
@@ -11,7 +11,7 @@
 setup(
     name='pylocus',
 
-    version='0.0.3',
+    version='0.0.4',
 
     description='Localization Package',
 
@@ -58,7 +58,7 @@
     # your project is installed. For an analysis of "install_requires" vs pip's
     # requirements files see:
     # https://packaging.python.org/en/latest/requirements.html
-    install_requires=['numpy', 'cvxpy'],
+    install_requires=['matplotlib', 'numpy'],
 
     # List additional groups of dependencies here (e.g. development
     # dependencies). You can install these using the following syntax,
diff --git a/test/test_common.py b/test/test_common.py
index 08d4508..66b306e 100644
--- a/test/test_common.py
+++ b/test/test_common.py
@@ -35,7 +35,6 @@ def test_multiple(self):
             for i in range(self.n_it): # seed 381 used to fail.
                 np.random.seed(i)
                 self.test_zero_noise(it=i)
-                self.test_zero_noise_relaxed()
 
         def test_zero_noise(self, it=0):
             print('TestCommon:test_zero_noise')
@@ -46,24 +45,4 @@ def test_zero_noise(self, it=0):
                         points_estimate = self.call_method(method)
                         if points_estimate is None:
                             continue
-                        error = np.linalg.norm(self.pts.points - points_estimate) 
-                        self.assertTrue(error < self.eps, 
-                                        'error (method={}, it={}, N={}, d={}): {} not smaller than {}'.format(method, it, N, d, error, self.eps))
-
-        def test_zero_noise_relaxed(self):
-            print('TestCommon:test_zero_noise_relaxed')
-            success = 0
-            total = 0
-            for N in self.N_relaxed:
-                for d in (2, 3):
-                    self.create_points(N, d)
-                    for method in self.methods: 
-                        points_estimate = self.call_method(method)
-                        if points_estimate is None:
-                            continue
-                        error = np.linalg.norm(self.pts.points - points_estimate) 
-                        if error < self.eps:
-                            success +=1
-                        total +=1
-            rate = success/total*100
-            self.assertTrue(rate >= self.success_rate, 'noiseless success rate below {}: {}'.format(self.success_rate, rate))
+                        np.testing.assert_allclose(self.pts.points, points_estimate, self.eps)
diff --git a/test/test_completion.py b/test/test_completion.py
index 3348752..9d47965 100644
--- a/test/test_completion.py
+++ b/test/test_completion.py
@@ -12,7 +12,7 @@ def setUp(self):
         self.pts = PointSet(N=5, d=3)
         self.pts.set_points(mode='convex')
         W = np.ones(self.pts.edm.shape)
-        indices_missing = [[4,1],[0,3]]
+        indices_missing = ([4,1],[0,3])
         W[indices_missing] = 0.0
         W = np.multiply(W, W.T)
         np.fill_diagonal(W, 0.0)
diff --git a/test/test_srls.py b/test/test_srls.py
index 0c881eb..75dfa91 100644
--- a/test/test_srls.py
+++ b/test/test_srls.py
@@ -3,12 +3,15 @@
 
 import unittest
 import numpy as np
-from .test_common import BaseCommon
+import sys
+sys.path.append('./test/')
+
+from test_common import BaseCommon
 
 from pylocus.point_set import PointSet, create_from_points
 from pylocus.simulation import create_noisy_edm
 from pylocus.algorithms import reconstruct_srls
-from pylocus.lateration import SRLS, get_lateration_parameters
+from pylocus.lateration import SRLS, get_lateration_parameters, GeometryError
 
 class TestSRLS(BaseCommon.TestAlgorithms):
     def setUp(self):
@@ -28,13 +31,13 @@ def create_points(self, N=10, d=2):
     def call_method(self, method=''):
         print('TestSRLS:call_method')
         if method == '' or method == 'normal':
-            return reconstruct_srls(self.pts.edm, self.pts.points, n=self.n,
+            return reconstruct_srls(self.pts.edm, self.pts.points, 
                                     W=np.ones(self.pts.edm.shape))
         elif method == 'rescale':
-            return reconstruct_srls(self.pts.edm, self.pts.points, n=self.n,
+            return reconstruct_srls(self.pts.edm, self.pts.points, 
                                     W=np.ones(self.pts.edm.shape), rescale=True)
         elif method == 'fixed' and self.pts.d == 3:
-            return reconstruct_srls(self.pts.edm, self.pts.points, n=self.n,
+            return reconstruct_srls(self.pts.edm, self.pts.points, 
                                     W=np.ones(self.pts.edm.shape), rescale=False,
                                     z=self.pts.points[0, 2])
 
@@ -100,7 +103,7 @@ def test_srls_rescale(self):
 
         # Normal ranging
         x_srls = SRLS(anchors, w, r2)
-        self.assertTrue(np.allclose(x, x_srls))
+        np.testing.assert_allclose(x, x_srls)
 
         # Rescaled ranging
         x_srls_resc, scale = SRLS(anchors, w, sigma * r2, rescale=True)
@@ -111,11 +114,21 @@ def test_srls_fixed(self):
         print('TestSRLS:test_srls_fixed')
         self.create_points(N=10, d=3)
         zreal = self.pts.points[0, 2]
-        xhat = reconstruct_srls(self.pts.edm, self.pts.points, n=self.n, 
+        xhat = reconstruct_srls(self.pts.edm, self.pts.points,  
                                 W=np.ones(self.pts.edm.shape), rescale=False, 
-                                z=self.pts.points[0, 2])
-        self.assertEqual(xhat[0, 2], zreal)
-        np.testing.assert_allclose(xhat, self.pts.points)
+                                z=zreal)
+        if xhat is not None:
+            np.testing.assert_allclose(xhat[0, 2], zreal)
+            np.testing.assert_allclose(xhat, self.pts.points)
+
+    def test_srls_fail(self):
+        anchors = np.array([[11.881,  3.722,  1.5  ],
+                            [11.881,  14.85,   1.5 ],
+                            [11.881,  7.683,  1.5  ]])
+        w = np.ones((3, 1))
+        distances = [153.32125426, 503.96654466, 234.80741129] 
+        z = 1.37
+        self.assertRaises(GeometryError, SRLS, anchors, w, distances, False, z)
 
     def zero_weights(self, noise=0.1):
         index = np.arange(self.n)
@@ -130,7 +143,7 @@ def zero_weights(self, noise=0.1):
         edm_anchors = np.delete(edm_noisy, indices, axis=0)
         edm_anchors = np.delete(edm_anchors, indices, axis=1)
         missing_anchors = reconstruct_srls(
-            edm_anchors, points_missing.points, n=self.n, W=None,
+            edm_anchors, points_missing.points, W=None,
             print_out=False)
 
         # missing distances
@@ -139,7 +152,7 @@ def zero_weights(self, noise=0.1):
         weights[index, indices] = 0.0
 
         missing_distances = reconstruct_srls(
-            edm_noisy, self.pts.points, n=self.n, W=weights)
+            edm_noisy, self.pts.points, W=weights)
         left_distances = np.delete(range(self.pts.N), indices)
 
         self.assertTrue(np.linalg.norm(
diff --git a/test/test_srls_fail.py b/test/test_srls_fail.py
index c96b429..2a6c997 100644
--- a/test/test_srls_fail.py
+++ b/test/test_srls_fail.py
@@ -46,10 +46,10 @@ def setUp(self):
         self.weights[1:, 0] = weights
 
     def test_strange_case(self):
-        xhat = reconstruct_srls(self.edm, self.all_points, n=1,
+        xhat = reconstruct_srls(self.edm, self.all_points, 
                                 W=self.weights)
         print(xhat[0])
-        xhat = reconstruct_srls(self.edm, self.all_points, n=1,
+        xhat = reconstruct_srls(self.edm, self.all_points, 
                                 W=self.weights, rescale=False,
                                 z=self.all_points[0, 2])
         print(xhat[0])