diff --git a/README.rst b/README.rst index 9f004c6..3ad708e 100644 --- a/README.rst +++ b/README.rst @@ -25,6 +25,13 @@ To perform a local install on your computer, run from this folder level (where s pip install -e . This installs the package using symbolic links, avoiding the need for a reinstall whenever the source code is changed. +If you use conda, then + +.. code-block:: bash + + conda develop . + +does the same trick. Install ******* @@ -49,21 +56,6 @@ Depending on which parts of the project you are using, you might need to install Documentation ************* -You can import the package and its submodules as usual, for example: - -.. code-block:: python - - from pylocus import algorithms - - Xhat = algorithms.reconstruct_mds(...) - -or - -.. code-block:: python - - from pylocus.algorithms import * - - Xhat = reconstruct_mds(...) - +See the tutorials folder for some exmaple scripts on how to use this package. More scripts will be added soon. -Documentation of all available functionalities can be found on ReadTheDocs: http://pylocus.readthedocs.org/en/latest/ +This is a constantly growing package and documentation is work-in-progress. The current version can be found on ReadTheDocs: http://pylocus.readthedocs.org/en/latest/ diff --git a/pylocus/algorithms.py b/pylocus/algorithms.py index 5a80b4d..01ae1cb 100644 --- a/pylocus/algorithms.py +++ b/pylocus/algorithms.py @@ -2,36 +2,53 @@ # module ALGORITHMS import numpy as np - -def execute_method(method, noisy_edm=None, real_points=None, W=None, z=None, **kwargs): +IMPLEMENTED_METHODS = ['MDS', + 'MDSoptspace', + 'MDSalternate', + 'SDR', + 'ACD', + 'dwMDS', + 'SRLS'] + + +def execute_method(method, measured_matrix=None, all_points=None, W=None, **kwargs): + if method not in IMPLEMENTED_METHODS: + raise NotImplementedError( + 'method {} is not implemented.'.format(method)) if method == 'MDS': xhat = reconstruct_mds( - noisy_edm, real_points=real_points, method='geometric') + measured_matrix, all_points=all_points, method='geometric') if method == 'MDSoptspace': - xhat = reconstruct_mds(noisy_edm, real_points=real_points, + xhat = reconstruct_mds(measured_matrix, all_points=all_points, method='geometric', mask=W, completion='optspace', print_out=False) if method == 'MDSalternate': - xhat = reconstruct_mds(noisy_edm, real_points=real_points, + xhat = reconstruct_mds(measured_matrix, all_points=all_points, method='geometric', mask=W, completion='alternate', print_out=False) if method == 'SDR': - x_SDRold, EDMbest = reconstruct_sdp( - noisy_edm, W=W, real_points=real_points) + x_SDRold, __ = reconstruct_sdp( + measured_matrix, W=W, all_points=all_points) + # TODO # Added to avoid strange "too large to be a matrix" error - N, d = real_points.shape + N, d = all_points.shape xhat = np.zeros((N, d)) xhat[:, :] = x_SDRold if method == 'ACD': X0 = kwargs.get('X0', None) - xhat, costs = reconstruct_acd(noisy_edm, W=W, X0=X0) + if X0 is None: + raise NameError('Need to provide X0 for method ACD.') + xhat, costs = reconstruct_acd(measured_matrix, W=W, X0=X0) if method == 'dwMDS': X0 = kwargs.pop('X0', None) - xhat, costs = reconstruct_dwmds(noisy_edm, W=W, X0=X0, **kwargs) + if X0 is None: + raise NameError('Need to provide X0 for method dwMDS.') + xhat, costs = reconstruct_dwmds(measured_matrix, W=W, X0=X0, **kwargs) if method == 'SRLS': - n = kwargs.get('n', None) + n = kwargs.get('n', 1) + z = kwargs.get('z', None) rescale = kwargs.get('rescale', False) - xhat = reconstruct_srls(noisy_edm, real_points, + xhat = reconstruct_srls(measured_matrix, all_points, n=n, W=W, rescale=rescale, z=z) return xhat @@ -95,16 +112,16 @@ def centralize(X): return X_transformed, R, t, c -def reconstruct_emds(edm, Om, real_points, method=None, **kwargs): +def reconstruct_emds(edm, Om, all_points, method=None, **kwargs): """ Reconstruct point set using E(dge)-MDS. """ from .point_set import dm_from_edm - N = real_points.shape[0] - d = real_points.shape[1] + N = all_points.shape[0] + d = all_points.shape[1] dm = dm_from_edm(edm) if method is None: from .mds import superMDS - Xhat, __ = superMDS(real_points[0, :], N, d, Om=Om, dm=dm) + Xhat, __ = superMDS(all_points[0, :], N, d, Om=Om, dm=dm) else: C = kwargs.get('C', None) b = kwargs.get('b', None) @@ -115,24 +132,24 @@ def reconstruct_emds(edm, Om, real_points, method=None, **kwargs): if method == 'iterative': from .mds import iterativeEMDS Xhat, __ = iterativeEMDS( - real_points[0, :], N, d, KE=KE_noisy, C=C, b=b) + all_points[0, :], N, d, KE=KE_noisy, C=C, b=b) elif method == 'relaxed': from .mds import relaxedEMDS Xhat, __ = relaxedEMDS( - real_points[0, :], N, d, KE=KE_noisy, C=C, b=b) + all_points[0, :], N, d, KE=KE_noisy, C=C, b=b) else: raise NameError('Undefined method', method) - Y, R, t, c = procrustes(real_points, Xhat, scale=False) + Y, R, t, c = procrustes(all_points, Xhat, scale=False) return Y -def reconstruct_cdm(dm, absolute_angles, real_points, W=None): +def reconstruct_cdm(dm, absolute_angles, all_points, W=None): """ Reconstruct point set from angle and distance measurements, using coordinate difference matrices. """ from pylocus.point_set import dmi_from_V, sdm_from_dmi, get_V from pylocus.mds import signedMDS - N = real_points.shape[0] + N = all_points.shape[0] V = get_V(absolute_angles, dm) @@ -146,17 +163,17 @@ def reconstruct_cdm(dm, absolute_angles, real_points, W=None): points_y = signedMDS(sdmy, W) Xhat = np.c_[points_x, points_y] - Y, R, t, c = procrustes(real_points, Xhat, scale=False) + Y, R, t, c = procrustes(all_points, Xhat, scale=False) return Y -def reconstruct_mds(edm, real_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, n=1): """ Reconstruct point set using MDS and matrix completion algorithms. """ from .point_set import dm_from_edm from .mds import MDS - N = real_points.shape[0] - d = real_points.shape[1] + N = all_points.shape[0] + d = all_points.shape[1] if mask is not None: edm_missing = np.multiply(edm, mask) if completion == 'optspace': @@ -174,30 +191,30 @@ def reconstruct_mds(edm, real_points, completion='optspace', mask=None, method=' print('{}: relative error:{}'.format(completion, err)) edm = edm_complete Xhat = MDS(edm, d, method, False).T - Y, R, t, c = procrustes(real_points[n:], Xhat, True) - #Y, R, t, c = procrustes(real_points, Xhat, True) + Y, R, t, c = procrustes(all_points[n:], Xhat, True) + #Y, R, t, c = procrustes(all_points, Xhat, True) return Y -def reconstruct_sdp(edm, real_points, W=None, print_out=False, lamda=1000, **kwargs): +def reconstruct_sdp(edm, all_points, W=None, print_out=False, lamda=1000, **kwargs): """ Reconstruct point set using semi-definite rank relaxation. """ from .edm_completion import semidefinite_relaxation edm_complete = semidefinite_relaxation( edm, lamda=lamda, W=W, print_out=print_out, **kwargs) - Xhat = reconstruct_mds(edm_complete, real_points, method='geometric') + Xhat = reconstruct_mds(edm_complete, all_points, method='geometric') return Xhat, edm_complete -def reconstruct_srls(edm, real_points, W=None, print_out=False, n=1, rescale=False, +def reconstruct_srls(edm, all_points, W=None, print_out=False, n=1, rescale=False, z=None): """ Reconstruct point set using S(quared)R(ange)L(east)S(quares) method. """ from .lateration import SRLS, get_lateration_parameters - Y = real_points.copy() + Y = all_points.copy() indices = range(n) for index in indices: - anchors, w, r2 = get_lateration_parameters(real_points, indices, index, + anchors, w, r2 = get_lateration_parameters(all_points, indices, index, edm, W) if print_out: print('SRLS parameters:') @@ -216,6 +233,9 @@ def reconstruct_acd(edm, X0, W=None, print_out=False, tol=1e-10, sweeps=10): """ Reconstruct point set using alternating coordinate descent. :param X0: Nxd matrix of starting points. + :param tol: Stopping criterion: if the stepsize in all coordinate directions + is less than tol for 2 consecutive sweeps, we stop. + :param sweep: Maximum number of sweeps. One sweep goes through all coordintaes and points once. """ def get_unique_delta(delta, i, coord, print_out=False): @@ -306,8 +326,13 @@ def loop_point(i): def reconstruct_dwmds(edm, X0, W=None, n=None, r=None, X_bar=None, print_out=False, tol=1e-10, sweeps=100): """ Reconstruct point set using d(istributed)w(eighted) MDS. + Refer to paper "Distributed Weighted-Multidimensional Scaling for Node Localization in Sensor Networks" for + implementation details (doi.org/10.1145/1138127.1138129) + :param X0: Nxd matrix of starting points. :param n: Number of points of unknown position. The first n points in X0 and edm are considered unknown. + :param tol: Stopping criterion: when the cost is below this level, we stop. + :param sweeps: Maximum number of sweeps. """ from .basics import get_edm from .distributed_mds import get_b, get_Si diff --git a/pylocus/algorithms.py.orig b/pylocus/algorithms.py.orig new file mode 100644 index 0000000..ee35009 --- /dev/null +++ b/pylocus/algorithms.py.orig @@ -0,0 +1,388 @@ +#!/usr/bin/env python +# module ALGORITHMS +import numpy as np + +<<<<<<< HEAD + +def execute_method(method, noisy_edm=None, real_points=None, W=None, z=None, **kwargs): +======= +IMPLEMENTED_METHODS = ['MDS', + 'MDSoptspace', + 'MDSalternate', + 'SDR', + 'ACD', + 'dwMDS', + 'SRLS'] + + +def execute_method(method, measured_matrix=None, all_points=None, W=None, **kwargs): + if method not in IMPLEMENTED_METHODS: + raise NotImplementedError( + 'method {} is not implemented.'.format(method)) +>>>>>>> 5790607c7fc2d3683148ed8d79d54e3b658262c0 + if method == 'MDS': + xhat = reconstruct_mds( + measured_matrix, all_points=all_points, method='geometric') + if method == 'MDSoptspace': + xhat = reconstruct_mds(measured_matrix, all_points=all_points, + method='geometric', mask=W, + completion='optspace', print_out=False) + if method == 'MDSalternate': + xhat = reconstruct_mds(measured_matrix, all_points=all_points, + method='geometric', mask=W, + completion='alternate', print_out=False) + if method == 'SDR': + x_SDRold, __ = reconstruct_sdp( + measured_matrix, W=W, all_points=all_points) + # TODO + # Added to avoid strange "too large to be a matrix" error + N, d = all_points.shape + xhat = np.zeros((N, d)) + xhat[:, :] = x_SDRold + if method == 'ACD': + X0 = kwargs.get('X0', None) + if X0 is None: + raise NameError('Need to provide X0 for method ACD.') + xhat, costs = reconstruct_acd(measured_matrix, W=W, X0=X0) + if method == 'dwMDS': + X0 = kwargs.pop('X0', None) + if X0 is None: + raise NameError('Need to provide X0 for method dwMDS.') + xhat, costs = reconstruct_dwmds(measured_matrix, W=W, X0=X0, **kwargs) + if method == 'SRLS': + n = kwargs.get('n', 1) + rescale = kwargs.get('rescale', False) +<<<<<<< HEAD + xhat = reconstruct_srls(noisy_edm, real_points, + n=n, W=W, rescale=rescale, z=z) +======= + xhat = reconstruct_srls(measured_matrix, all_points, + n=n, W=W, rescale=rescale) +>>>>>>> 5790607c7fc2d3683148ed8d79d54e3b658262c0 + return xhat + + +def classical_mds(D): + from .mds import MDS + return MDS(D, 1, 'geometric') + + +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 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. + """ + def centralize(X): + n = X.shape[0] + ones = np.ones((n, 1)) + return X - np.multiply(1 / n * np.dot(ones.T, X), ones) + m = anchors.shape[0] + N, d = X.shape + assert m >= d, 'Have to give at least d anchor nodes.' + X_m = X[N - m:, :] + ones = np.ones((m, 1)) + + mux = 1 / m * np.dot(ones.T, X_m) + muy = 1 / m * np.dot(ones.T, anchors) + sigmax = 1 / m * np.linalg.norm(X_m - mux)**2 + sigmaxy = 1 / m * np.dot((anchors - muy).T, X_m - mux) + try: + U, D, VT = np.linalg.svd(sigmaxy) + except np.LinAlgError: + print('strange things are happening...') + print(sigmaxy) + print(np.linalg.matrix_rank(sigmaxy)) + #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): + # print('switching') + # S[-1,-1] = -1.0 + # else: + # print('not switching') + # c = np.trace(np.dot(np.diag(D),S))/sigmax + # R = np.dot(U, np.dot(S,VT)) + if (scale): + c = np.trace(np.diag(D)) / sigmax + else: + c = np.trace(np.diag(D)) / sigmax + if (print_out): + print('Optimal scale would be: {}. Setting it to 1 now.'.format(c)) + c = 1.0 + 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 + return X_transformed, R, t, c + + +def reconstruct_emds(edm, Om, all_points, method=None, **kwargs): + """ Reconstruct point set using E(dge)-MDS. + """ + from .point_set import dm_from_edm + N = all_points.shape[0] + d = all_points.shape[1] + dm = dm_from_edm(edm) + if method is None: + from .mds import superMDS + Xhat, __ = superMDS(all_points[0, :], N, d, Om=Om, dm=dm) + else: + C = kwargs.get('C', None) + b = kwargs.get('b', None) + if C is None or b is None: + raise NameError( + 'Need constraints C and b for reconstruct_emds in iterative mode.') + KE_noisy = np.multiply(np.outer(dm, dm), Om) + if method == 'iterative': + from .mds import iterativeEMDS + Xhat, __ = iterativeEMDS( + all_points[0, :], N, d, KE=KE_noisy, C=C, b=b) + elif method == 'relaxed': + from .mds import relaxedEMDS + Xhat, __ = relaxedEMDS( + all_points[0, :], N, d, KE=KE_noisy, C=C, b=b) + else: + raise NameError('Undefined method', method) + Y, R, t, c = procrustes(all_points, Xhat, scale=False) + return Y + + +def reconstruct_cdm(dm, absolute_angles, all_points, W=None): + """ Reconstruct point set from angle and distance measurements, using coordinate difference matrices. + """ + from pylocus.point_set import dmi_from_V, sdm_from_dmi, get_V + from pylocus.mds import signedMDS + + N = all_points.shape[0] + + V = get_V(absolute_angles, dm) + + dmx = dmi_from_V(V, 0) + dmy = dmi_from_V(V, 1) + + sdmx = sdm_from_dmi(dmx, N) + sdmy = sdm_from_dmi(dmy, N) + + points_x = signedMDS(sdmx, W) + points_y = signedMDS(sdmy, W) + + Xhat = np.c_[points_x, points_y] + Y, R, t, c = procrustes(all_points, Xhat, scale=False) + return Y + + +def reconstruct_mds(edm, all_points, completion='optspace', mask=None, method='geometric', print_out=False, n=1): + """ Reconstruct point set using MDS and matrix completion algorithms. + """ + from .point_set import dm_from_edm + from .mds import MDS + N = all_points.shape[0] + d = all_points.shape[1] + if mask is not None: + edm_missing = np.multiply(edm, mask) + if completion == 'optspace': + from .edm_completion import optspace + edm_complete = optspace(edm_missing, d + 2) + elif completion == 'alternate': + from .edm_completion import rank_alternation + edm_complete, errs = rank_alternation( + edm_missing, d + 2, print_out=False, edm_true=edm) + else: + raise NameError('Unknown completion method {}'.format(completion)) + if (print_out): + err = np.linalg.norm(edm_complete - edm)**2 / \ + np.linalg.norm(edm)**2 + print('{}: relative error:{}'.format(completion, err)) + edm = edm_complete + Xhat = MDS(edm, d, method, False).T + Y, R, t, c = procrustes(all_points[n:], Xhat, True) + #Y, R, t, c = procrustes(all_points, Xhat, True) + return Y + + +def reconstruct_sdp(edm, all_points, W=None, print_out=False, lamda=1000, **kwargs): + """ Reconstruct point set using semi-definite rank relaxation. + """ + from .edm_completion import semidefinite_relaxation + edm_complete = semidefinite_relaxation( + edm, lamda=lamda, W=W, print_out=print_out, **kwargs) + Xhat = reconstruct_mds(edm_complete, all_points, method='geometric') + return Xhat, edm_complete + + +def reconstruct_srls(edm, all_points, W=None, print_out=False, n=1, rescale=False, + z=None): + """ Reconstruct point set using S(quared)R(ange)L(east)S(quares) method. + """ + from .lateration import SRLS, get_lateration_parameters + Y = all_points.copy() + indices = range(n) + for index in indices: + anchors, w, r2 = get_lateration_parameters(all_points, indices, index, + edm, W) + if print_out: + print('SRLS parameters:') + print('anchors', anchors) + print('w', w) + print('r2', r2) + + srls = SRLS(anchors, w, r2, rescale, z, print_out) + if rescale: + srls = srls[0] # second element of output is the scale + Y[index, :] = srls + return Y + + +def reconstruct_acd(edm, X0, W=None, print_out=False, tol=1e-10, sweeps=10): + """ Reconstruct point set using alternating coordinate descent. + + :param X0: Nxd matrix of starting points. + :param tol: Stopping criterion: if the stepsize in all coordinate directions + is less than tol for 2 consecutive sweeps, we stop. + :param sweep: Maximum number of sweeps. One sweep goes through all coordintaes and points once. + """ + + def get_unique_delta(delta, i, coord, print_out=False): + delta_unique = delta.copy() + ftest = [] + for de in delta_unique: + X_ktest = X_k.copy() + X_ktest[i, coord] += de + ftest.append(f(X_ktest, edm, W)) + # choose delta_unique with lowest cost. + if print_out: + print(ftest) + print(delta_unique) + delta_unique = delta_unique[ftest == min(ftest)] + if print_out: + print(delta_unique) + if len(delta_unique) > 1: + # if multiple delta_uniques give the same cost, choose the biggest step. + delta_unique = max(delta_unique) + return delta_unique + + def sweep(): + for i in range(N): + loop_point(i) + + def loop_coordinate(coord, i): + delta = get_step_size(i, coord, X_k, edm, W) + if len(delta) > 1: + delta_unique = get_unique_delta(delta, i, coord) + elif len(delta) == 1: + delta_unique = delta[0] + else: + print('Warning: did not find delta!', delta) + delta_unique = 0.0 + try: + X_k[i, coord] += delta_unique + except: + get_unique_delta(delta, i, coord, print_out=True) + raise + + cost_this = f(X_k, edm, W) + costs.append(cost_this) + + if delta_unique <= tol: + coordinates_converged[i, coord] += 1 + else: + if print_out: + print('======= coordinate {} of {} not yet converged'.format(coord, i)) + coordinates_converged[i, coord] = 0 + + def loop_point(i): + coord_counter = 0 + while not (coordinates_converged[i] >= 2).all(): + if print_out: + print('==== point {}'.format(i)) + coord_counter += 1 + for coord in range(d): + if print_out: + print('======= coord {}'.format(coord)) + loop_coordinate(coord, i) + if coord_counter > coord_n_it: + break + from .distributed_mds import get_step_size, f + + N, d = X0.shape + if W is None: + W = np.ones((N, N)) - np.eye(N) + + X_k = X0.copy() + + costs = [] + coordinates_converged = np.zeros(X_k.shape) + coord_n_it = 3 + for sweep_counter in range(sweeps): + if print_out: + print('= sweep', sweep_counter) + sweep() + if (coordinates_converged >= 2).all(): + if (print_out): + print('acd: all coordinates converged after {} sweeps.'.format( + sweep_counter)) + return X_k, costs + if (print_out): + print('acd: did not converge after {} sweeps'.format(sweep_counter + 1)) + return X_k, costs + + +def reconstruct_dwmds(edm, X0, W=None, n=None, r=None, X_bar=None, print_out=False, tol=1e-10, sweeps=100): + """ Reconstruct point set using d(istributed)w(eighted) MDS. + + Refer to paper "Distributed Weighted-Multidimensional Scaling for Node Localization in Sensor Networks" for + implementation details (doi.org/10.1145/1138127.1138129) + + :param X0: Nxd matrix of starting points. + :param n: Number of points of unknown position. The first n points in X0 and edm are considered unknown. + :param tol: Stopping criterion: when the cost is below this level, we stop. + :param sweeps: Maximum number of sweeps. + """ + from .basics import get_edm + from .distributed_mds import get_b, get_Si + + N, d = X0.shape + if r is None and n is None: + raise ValueError('either r or n have to be given.') + elif n is None: + n = r.shape[0] + + if W is None: + W = np.ones((N, N)) - np.eye(N) + + X_k = X0.copy() + + costs = [] + # don't have to ignore i=j, because W[i,i] is zero. + a = np.sum(W[:n, :n], axis=1).flatten() + 2 * \ + np.sum(W[:n, n:], axis=1).flatten() + if r is not None: + a += r.flatten() + for k in range(sweeps): + S = 0 + for i in range(n): + edm_estimated = get_edm(X_k) + bi = get_b(i, edm_estimated, W, edm, n) + if r is not None and X_bar is not None: + X_k[i] = 1 / a[i] * (r[i] * X_bar[i, :] + + X_k.T.dot(bi).flatten()) + Si = get_Si(i, edm_estimated, edm, W, n, r, X_bar[i], X_k[i]) + else: + X_k[i] = 1 / a[i] * X_k.T.dot(bi).flatten() + Si = get_Si(i, edm_estimated, edm, W, n) + S += Si + costs.append(S) + if k > 1 and abs(costs[-1] - costs[-2]) < tol: + if (print_out): + print('dwMDS: converged after', k) + break + return X_k, costs + + +if __name__ == "__main__": + print('nothing happens when running this module.') diff --git a/pylocus/basics.py b/pylocus/basics.py index 9a5cc5b..cdf5e13 100644 --- a/pylocus/basics.py +++ b/pylocus/basics.py @@ -104,7 +104,7 @@ def divide_where_nonzero(divide_this, by_this): zero_mask = (by_this == 0) if zero_mask.size: result[zero_mask] = 0.0 - nonzero_mask = mask # creates a view + nonzero_mask = zero_mask # creates a view # overwrites memory of zero_mask np.logical_not(zero_mask, out=nonzero_mask) result[nonzero_mask] = divide_this[nonzero_mask] / \ diff --git a/pylocus/edm_completion.py b/pylocus/edm_completion.py index 956b71f..ea73812 100644 --- a/pylocus/edm_completion.py +++ b/pylocus/edm_completion.py @@ -12,13 +12,11 @@ def optspace(edm_missing, rank, niter=500, tol=1e-6, print_out=False): Uses OptSpace algorithm to complete and denoise EDM. The problem being solved is X,S,Y = argmin_(X,S,Y) || W ° (D - XSY') ||_F^2 - Args: - edm_missing: EDM with 0 where no measurement was taken - rank: expected rank of complete EDM - niter, tol: see opt_space module for description. + :param edm_missing: EDM with 0 where no measurement was taken. + :param rank: expected rank of complete EDM. + :param niter, tol: see opt_space module for description. - Returns: - Completed matrix. + :return: Completed matrix. """ from .opt_space import opt_space N = edm_missing.shape[0] @@ -30,19 +28,17 @@ def optspace(edm_missing, rank, niter=500, tol=1e-6, print_out=False): def rank_alternation(edm_missing, rank, niter=50, print_out=False, edm_true=None): - """Complete missing EDM entries using rank alternation. + """Complete and denoise EDM using rank alternation. Iteratively impose rank and strucutre to complete marix entries - Args: - edm_missing: EDM with 0 where no measurement was taken - rank: expected rank of complete EDM - niter: maximum number of iterations - edm: if given, the relative EDM error is tracked. + :param edm_missing: EDM with 0 where no measurement was taken. + :param rank: expected rank of complete EDM. + :param niter: maximum number of iterations. + :param edm: if given, the relative EDM error is tracked. - Returns: - Completed matrix and array of errors (empty if no true edm is given). - The matrix is of the correct structure, but might not have the right measured entries. + :return: Completed matrix and array of errors (empty if no true edm is given). + The matrix is of the correct structure, but might not have the right measured entries. """ from pylocus.basics import low_rank_approximation @@ -69,6 +65,22 @@ def rank_alternation(edm_missing, rank, niter=50, print_out=False, edm_true=None def semidefinite_relaxation(edm_missing, lamda, W=None, print_out=False, **kwargs): + """Complete and denoise EDM using semidefinite relaxation. + + Returns solution to the relaxation of the following problem: + D = argmin || W * (D - edm_missing) || + s.t. D is EDM + where edm_missing is measured matrix, W is a weight matrix, and * is pointwise multiplication. + + Refer to paper "Euclidean Distance Matrices - Essential Theory, Algorithms and Applications", + Algorithm 5, for details. (https://www.doi.org/%2010.1109/MSP.2015.2398954) + + :param edm_missing: EDM with 0 where no measurement was taken. + :param lamda: Regularization parameter. + :param W: Optional mask. If no mask is given, a binary mask is created based on missing elements of edm_missing. + If mask is given + :param kwargs: more options passed to the solver. See cvxpy documentation for all options. + """ from .algorithms import reconstruct_mds def kappa(gram): @@ -93,15 +105,14 @@ def kappa_cvx(gram, n): V = np.c_[-np.ones((n - 1, 1)) / np.sqrt(n), np.eye(n - 1) - np.ones((n - 1, n - 1)) / (n + np.sqrt(n))].T - # Creates a n-1 by n-1 positive semidefinite variable. H = Variable((n - 1, n - 1), PSD=True) G = V * H * V.T # * is overloaded - print(G.shape) edm_optimize = kappa_cvx(G, n) if method == 'maximize': obj = Maximize(trace(H) - lamda * norm(multiply(W, (edm_optimize - edm_missing)), p=1)) + # TODO: add a reference to paper where "minimize" is used instead of maximize. elif method == 'minimize': obj = Minimize(trace(H) + lamda * norm(multiply(W, (edm_optimize - edm_missing)), p=1)) @@ -131,12 +142,36 @@ def kappa_cvx(gram, n): def completion_acd(edm, X0, W=None, tol=1e-6, sweeps=3): + """ Complete an denoise EDM using alternating decent. + + The idea here is to simply run reconstruct_acd for a few iterations, + yieding a position estimate, which can in turn be used + to get a completed and denoised edm. + + :param edm: noisy matrix (NxN) + :param X0: starting points (Nxd) + :param W: optional weight matrix. + :param tol: Stopping criterion of iterative algorithm. + :param sweeps: Maximum number of sweeps. + """ from .algorithms import reconstruct_acd Xhat, costs = reconstruct_acd(edm, X0, W, tol=tol, sweeps=sweeps) return get_edm(Xhat) def completion_dwmds(edm, X0, W=None, tol=1e-10, sweeps=100): + """ Complete an denoise EDM using dwMDS. + + The idea here is to simply run reconstruct_dwmds for a few iterations, + yieding a position estimate, which can in turn be used + to get a completed and denoised edm. + + :param edm: noisy matrix (NxN) + :param X0: starting points (Nxd) + :param W: optional weight matrix. + :param tol: Stopping criterion of iterative algorithm. + :param sweeps: Maximum number of sweeps. + """ from .algorithms import reconstruct_dwmds Xhat, costs = reconstruct_dwmds(edm, X0, W, n=1, tol=tol, sweeps=sweeps) return get_edm(Xhat) diff --git a/pylocus/edm_completion.py.orig b/pylocus/edm_completion.py.orig new file mode 100644 index 0000000..fce3950 --- /dev/null +++ b/pylocus/edm_completion.py.orig @@ -0,0 +1,185 @@ +#!/usr/bin/env python +# module EDM_COMPLETION +import numpy as np +from cvxpy import * + +from pylocus.basics import get_edm + + +def optspace(edm_missing, rank, niter=500, tol=1e-6, print_out=False): + """Complete and denoise EDM using OptSpace algorithm. + + Uses OptSpace algorithm to complete and denoise EDM. The problem being solved is + X,S,Y = argmin_(X,S,Y) || W ° (D - XSY') ||_F^2 + + :param edm_missing: EDM with 0 where no measurement was taken. + :param rank: expected rank of complete EDM. + :param niter, tol: see opt_space module for description. + + :return: Completed matrix. + """ + from .opt_space import opt_space + N = edm_missing.shape[0] + X, S, Y, __ = opt_space(edm_missing, r=rank, niter=niter, + tol=tol, print_out=print_out) + edm_complete = X.dot(S.dot(Y.T)) + edm_complete[range(N), range(N)] = 0.0 + return edm_complete + + +def rank_alternation(edm_missing, rank, niter=50, print_out=False, edm_true=None): + """Complete and denoise EDM using rank alternation. + + Iteratively impose rank and strucutre to complete marix entries + + :param edm_missing: EDM with 0 where no measurement was taken. + :param rank: expected rank of complete EDM. + :param niter: maximum number of iterations. + :param edm: if given, the relative EDM error is tracked. + + :return: Completed matrix and array of errors (empty if no true edm is given). + The matrix is of the correct structure, but might not have the right measured entries. + + """ + from pylocus.basics import low_rank_approximation + errs = [] + N = edm_missing.shape[0] + edm_complete = edm_missing.copy() + edm_complete[edm_complete == 0] = np.mean(edm_complete[edm_complete > 0]) + for i in range(niter): + # impose matrix rank + edm_complete = low_rank_approximation(edm_complete, rank) + + # impose known entries + edm_complete[edm_missing > 0] = edm_missing[edm_missing > 0] + + # impose matrix structure + edm_complete[range(N), range(N)] = 0.0 + edm_complete[edm_complete < 0] = 0.0 + edm_complete = 0.5 * (edm_complete + edm_complete.T) + + if edm_true is not None: + err = np.linalg.norm(edm_complete - edm_true) + errs.append(err) + return edm_complete, errs + + +def semidefinite_relaxation(edm_missing, lamda, W=None, print_out=False, **kwargs): + """Complete and denoise EDM using semidefinite relaxation. + + Returns solution to the relaxation of the following problem: + D = argmin || W * (D - edm_missing) || + s.t. D is EDM + where edm_missing is measured matrix, W is a weight matrix, and * is pointwise multiplication. + + Refer to paper "Euclidean Distance Matrices - Essential Theory, Algorithms and Applications", + Algorithm 5, for details. (https://www.doi.org/%2010.1109/MSP.2015.2398954) + + :param edm_missing: EDM with 0 where no measurement was taken. + :param lamda: Regularization parameter. + :param W: Optional mask. If no mask is given, a binary mask is created based on missing elements of edm_missing. + If mask is given + :param kwargs: more options passed to the solver. See cvxpy documentation for all options. + """ + from .algorithms import reconstruct_mds + + def kappa(gram): + n = len(gram) + e = np.ones(n) + return np.outer(np.diag(gram), e) + np.outer(e, np.diag(gram).T) - 2 * gram + + def kappa_cvx(gram, n): + e = np.ones((n, 1)) + return reshape(diag(gram), (n, 1)) * e.T + e * reshape(diag(gram), (1, n)) - 2 * gram + + method = kwargs.pop('method', 'maximize') + options = {'solver': 'CVXOPT'} + options.update(kwargs) + + if W is None: + W = (edm_missing > 0) + else: + W[edm_missing == 0] = 0.0 + + n = edm_missing.shape[0] + V = np.c_[-np.ones((n - 1, 1)) / np.sqrt(n), np.eye(n - 1) - + np.ones((n - 1, n - 1)) / (n + np.sqrt(n))].T + + H = Variable((n - 1, n - 1), PSD=True) + G = V * H * V.T # * is overloaded + edm_optimize = kappa_cvx(G, n) + + if method == 'maximize': + obj = Maximize(trace(H) - lamda * +<<<<<<< HEAD + norm(multiply(W, (edm_optimize - edm_missing)), p=1)) +======= + norm(multiply(W, (edm_optimize - edm_missing)))) + # TODO: add a reference to paper where "minimize" is used instead of maximize. +>>>>>>> 5790607c7fc2d3683148ed8d79d54e3b658262c0 + elif method == 'minimize': + obj = Minimize(trace(H) + lamda * + norm(multiply(W, (edm_optimize - edm_missing)), p=1)) + + prob = Problem(obj) + + total = prob.solve(**options) + if print_out: + print('total cost:', total) + print('SDP status:', prob.status) + + if H.value is not None: + Gbest = V.dot(H.value).dot(V.T) + if print_out: + print('eigenvalues:', np.sum(np.linalg.eigvals(Gbest)[2:])) + edm_complete = kappa(Gbest) + else: + edm_complete = edm_missing + + if (print_out): + if H.value is not None: + print('trace of H:', np.trace(H.value)) + print('other cost:', lamda * + norm(multiply(W, (edm_complete - edm_missing)), p=1).value) + + return np.array(edm_complete) + + +def completion_acd(edm, X0, W=None, tol=1e-6, sweeps=3): + """ Complete an denoise EDM using alternating decent. + + The idea here is to simply run reconstruct_acd for a few iterations, + yieding a position estimate, which can in turn be used + to get a completed and denoised edm. + + :param edm: noisy matrix (NxN) + :param X0: starting points (Nxd) + :param W: optional weight matrix. + :param tol: Stopping criterion of iterative algorithm. + :param sweeps: Maximum number of sweeps. + """ + from .algorithms import reconstruct_acd + Xhat, costs = reconstruct_acd(edm, X0, W, tol=tol, sweeps=sweeps) + return get_edm(Xhat) + + +def completion_dwmds(edm, X0, W=None, tol=1e-10, sweeps=100): + """ Complete an denoise EDM using dwMDS. + + The idea here is to simply run reconstruct_dwmds for a few iterations, + yieding a position estimate, which can in turn be used + to get a completed and denoised edm. + + :param edm: noisy matrix (NxN) + :param X0: starting points (Nxd) + :param W: optional weight matrix. + :param tol: Stopping criterion of iterative algorithm. + :param sweeps: Maximum number of sweeps. + """ + from .algorithms import reconstruct_dwmds + Xhat, costs = reconstruct_dwmds(edm, X0, W, n=1, tol=tol, sweeps=sweeps) + return get_edm(Xhat) + + +if __name__ == "__main__": + print('nothing happens when running this module.') diff --git a/pylocus/lateration.py b/pylocus/lateration.py index 8fbb9d9..ec44866 100644 --- a/pylocus/lateration.py +++ b/pylocus/lateration.py @@ -8,14 +8,14 @@ from pylocus.basics import assert_print, assert_all_print -def get_lateration_parameters(real_points, indices, index, edm, W=None): - """ Get parameters relevant for lateration from full real_points, edm and W. +def get_lateration_parameters(all_points, indices, index, edm, W=None): + """ Get parameters relevant for lateration from full all_points, edm and W. """ if W is None: W = np.ones(edm.shape) # delete points that are not considered anchors - anchors = np.delete(real_points, indices, axis=0) + anchors = np.delete(all_points, indices, axis=0) r2 = np.delete(edm[index, :], indices) w = np.delete(W[index, :], indices) @@ -167,7 +167,6 @@ def phi_prime(_lambda): # 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): - print('lower bound:', lower_bound) # brentq is considered the best rootfinding routine. try: lambda_opt = optimize.brentq(phi, lower_bound, inf, xtol=xtol) diff --git a/pylocus/simulation.py b/pylocus/simulation.py index f5d67d0..ff9f891 100644 --- a/pylocus/simulation.py +++ b/pylocus/simulation.py @@ -3,28 +3,28 @@ import numpy as np -def weights_one(N, noise, noise_edm): +def weights_one(N): weights = np.ones((N, N)) weights[range(N), range(N)] = 0.0 return weights -def weights_zero(N, noise, noise_edm): - weights = weights_one(N, noise, noise_edm) +def weights_zero(N): + weights = weights_one(N) weights[0, 1] = 0.0 weights[1, 0] = 0.0 return weights def weights_linear(N, noise, noise_edm): - weights = weights_one(N, noise, noise_edm) + weights = weights_one(N) weights[0, 1] = noise_edm / noise weights[1, 0] = weights[0, 1] return weights def weights_quadratic(N, noise, noise_edm): - weights = weights_one(N, noise, noise_edm) + weights = weights_one(N) weights[0, 1] = (noise_edm / noise)**2 weights[1, 0] = weights[0, 1] return weights @@ -99,5 +99,19 @@ def create_mask(N, method='all', nmissing=0): return weights +def create_weights(N, method='one', noise=0.0, noise_edm=None): + if method == 'one': + return weights_one(N) + elif method == 'zero': + return weights_zero(N) + elif method == 'linear': + return weights_linear(N, noise, noise_edm) + elif method == 'quadratic': + return weights_quadratic(N, noise, noise_edm) + else: + raise NameError("Unknown method", method) + + + if __name__ == "__main__": print('nothing happens when running this module.') diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..bb9b37f --- /dev/null +++ b/setup.cfg @@ -0,0 +1,208 @@ +[yapf] +################################################################# +## IMPORTANT PARAMETERS ## +################################################################# + +based_on_style = pep8 + +# The number of spaces required before a trailing comment. +spaces_before_comment = 2 + +# Set to True to prefer splitting before 'and' or 'or' rather than +# after. +split_before_logical_operator=True + +# Use the Tab character for indentation. +use_tabs=False + +# The column limit. +column_limit=100 + +# Indent width used for line continuations. +continuation_indent_width=4 + + +################################################################# +## NOT SO IMPORTANT PARAMETERS ## +################################################################# + +# Align closing bracket with visual indentation. +#align_closing_bracket_with_visual_indent=True + +# Allow dictionary keys to exist on multiple lines. For example: +# +# x = { +# ('this is the first element of a tuple', +# 'this is the second element of a tuple'): +# value, +# } +allow_multiline_dictionary_keys=True + +# Allow lambdas to be formatted on more than one line. +allow_multiline_lambdas=True + +# Allow splits before the dictionary value. +allow_split_before_dict_value=True + +# Insert a blank line before a class-level docstring. +blank_line_before_class_docstring=False + +# Insert a blank line before a 'def' or 'class' immediately nested +# within another 'def' or 'class'. For example: +# +# class Foo: +# # <------ this blank line +# def method(): +# ... +blank_line_before_nested_class_or_def=False + +# Do not split consecutive brackets. Only relevant when +# dedent_closing_brackets is set. For example: +# +# call_func_that_takes_a_dict( +# { +# 'key1': 'value1', +# 'key2': 'value2', +# } +# ) +# +# would reformat to: +# +# call_func_that_takes_a_dict({ +# 'key1': 'value1', +# 'key2': 'value2', +# }) +#coalesce_brackets=True + +# Put closing brackets on a separate line, dedented, if the bracketed +# expression can't fit in a single line. Applies to all kinds of brackets, +# including function definitions and calls. For example: +# +# config = { +# 'key1': 'value1', +# 'key2': 'value2', +# } # <--- this bracket is dedented and on a separate line +# +# time_series = self.remote_client.query_entity_counters( +# entity='dev3246.region1', +# key='dns.query_latency_tcp', +# transform=Transformation.AVERAGE(window=timedelta(seconds=60)), +# start_ts=now()-timedelta(days=3), +# end_ts=now(), +# ) # <--- this bracket is dedented and on a separate line +#dedent_closing_brackets=True + +# Place each dictionary entry onto its own line. +each_dict_entry_on_separate_line=True + +# The regex for an i18n comment. The presence of this comment stops +# reformatting of that line, because the comments are required to be +# next to the string they translate. +i18n_comment= + +# The i18n function call names. The presence of this function stops +# reformattting on that line, because the string it has cannot be moved +# away from the i18n comment. +i18n_function_call= + +# Indent the dictionary value if it cannot fit on the same line as the +# dictionary key. For example: +# +# config = { +# 'key1': +# 'value1', +# 'key2': value1 + +# value2, +# } +indent_dictionary_value=False + +# The number of columns to use for indentation. +indent_width=4 + +# Join short lines into one line. E.g., single line 'if' statements. +join_multiple_lines=True + +# Do not include spaces around selected binary operators. For example: +# +# 1 + 2 * 3 - 4 / 5 +# +# will be formatted as follows when configured with a value "*,/": +# +# 1 + 2*3 - 4/5 +# +no_spaces_around_selected_binary_operators=set() + +# Use spaces around default or named assigns. +spaces_around_default_or_named_assign=False + +# Use spaces around the power operator. +spaces_around_power_operator=False + +# Insert a space between the ending comma and closing bracket of a list, +# etc. +space_between_ending_comma_and_closing_bracket=True + +# Split before arguments if the argument list is terminated by a +# comma. +split_arguments_when_comma_terminated=False + +# Set to True to prefer splitting before '&', '|' or '^' rather than +# after. +split_before_bitwise_operator=True + +# Split before a dictionary or set generator (comp_for). For example, note +# the split before the 'for': +# +# foo = { +# variable: 'Hello world, have a nice day!' +# for variable in bar if variable != 42 +# } +split_before_dict_set_generator=True + +# Split after the opening paren which surrounds an expression if it doesn't +# fit on a single line. +split_before_expression_after_opening_paren=False + +# If an argument / parameter list is going to be split, then split before +# the first argument. +split_before_first_argument=False + +# Split named assignments onto individual lines. +split_before_named_assigns=True + +# The penalty for splitting right after the opening bracket. +split_penalty_after_opening_bracket=300 + +# The penalty for splitting the line after a unary operator. +split_penalty_after_unary_operator=10000 + +# The penalty for splitting right before an if expression. +split_penalty_before_if_expr=0 + +# The penalty of splitting the line around the '&', '|', and '^' +# operators. +split_penalty_bitwise_operator=300 + +# The penalty for characters over the column limit. +split_penalty_excess_character=4500 + +# The penalty incurred by adding a line split to the unwrapped line. The +# more line splits added the higher the penalty. +split_penalty_for_added_line_split=30 + +# The penalty of splitting a list of "import as" names. For example: +# +# from a_very_long_or_indented_module_name_yada_yad import (long_argument_1, +# long_argument_2, +# long_argument_3) +# +# would reformat to something like: +# +# from a_very_long_or_indented_module_name_yada_yad import ( +# long_argument_1, long_argument_2, long_argument_3) +split_penalty_import_names=0 + +# The penalty of splitting the line around the 'and' and 'or' +# operators. +split_penalty_logical_operator=300 + diff --git a/test/test_emds.py b/test/test_emds.py index 4b438bc..9d47657 100644 --- a/test/test_emds.py +++ b/test/test_emds.py @@ -29,10 +29,10 @@ def call_method(self, method=''): print('TestEMDS:call_method with', method) if method == '': return reconstruct_emds(self.pts.edm, Om=self.pts.Om, - real_points=self.pts.points) + all_points=self.pts.points) else: return reconstruct_emds(self.pts.edm, Om=self.pts.Om, - real_points=self.pts.points, C=self.C, b=self.b, method=method) + all_points=self.pts.points, C=self.C, b=self.b, method=method) if __name__ == "__main__": diff --git a/test/test_sdr.py b/test/test_sdr.py index 313f22f..53d719f 100644 --- a/test/test_sdr.py +++ b/test/test_sdr.py @@ -26,7 +26,7 @@ def create_points(self, N=10, d=3): def call_method(self, method=''): print('TestSDP:call_method') - Xhat, edm = reconstruct_sdp(self.pts.edm, real_points=self.pts.points, + Xhat, edm = reconstruct_sdp(self.pts.edm, all_points=self.pts.points, solver='CVXOPT', method='maximize') return Xhat @@ -43,7 +43,7 @@ def test_parameters(self): print('testing options', options, eps) self.eps = eps points_estimate, __ = reconstruct_sdp( - self.pts.edm, real_points=self.pts.points, method='maximize', **options) + self.pts.edm, all_points=self.pts.points, method='maximize', **options) error = np.linalg.norm(self.pts.points - points_estimate) self.assertTrue(error < self.eps, 'with options {} \nerror: {} not smaller than {}'.format( options, error, self.eps)) diff --git a/test/test_srls_fail.py b/test/test_srls_fail.py index 1c20dc7..c96b429 100644 --- a/test/test_srls_fail.py +++ b/test/test_srls_fail.py @@ -19,7 +19,7 @@ class TestSRLSFail(unittest.TestCase): """ trying to reproduce a fail of SRLS. """ def setUp(self): # from error logs: - self.real_points = np.array([[0.89, 2.87, 1.22], + self.all_points = np.array([[0.89, 2.87, 1.22], [2.12000000001, 0.0, 1.5], [5.83999999999, 2.5499998, 1.37999999], [5.83999999999, 4.6399997, 1.57000001], @@ -28,7 +28,7 @@ def setUp(self): [2.879999999999, 3.33999999999, 0.75], [5.759999999998, 3.46, 1.0], [3.54, 6.8700001, 1.139999]]) - self.anchors = self.real_points[1:, :] + self.anchors = self.all_points[1:, :] self.xhat = [0.8947106, 2.87644499, 1.22189597] distances = [1.1530589508606381, 1.6797066823110891, 1.1658051261713582, 0.42021624764667165, 1.1196561208517539, 0.44373006779135082, 1.3051970883541162, 2.6628083049963136] @@ -46,12 +46,12 @@ def setUp(self): self.weights[1:, 0] = weights def test_strange_case(self): - xhat = reconstruct_srls(self.edm, self.real_points, n=1, + xhat = reconstruct_srls(self.edm, self.all_points, n=1, W=self.weights) print(xhat[0]) - xhat = reconstruct_srls(self.edm, self.real_points, n=1, + xhat = reconstruct_srls(self.edm, self.all_points, n=1, W=self.weights, rescale=False, - z=self.real_points[0, 2]) + z=self.all_points[0, 2]) print(xhat[0]) if __name__ == "__main__": diff --git a/tutorials/README.md b/tutorials/README.md new file mode 100644 index 0000000..0370405 --- /dev/null +++ b/tutorials/README.md @@ -0,0 +1,6 @@ +# Tutorials for pylocus + +This folder contains some example scripts for how to use pylocus. Before running these scripts make sure +you have pylocus installed. + +- example_srls.py: example of how to use srls for localizing one point given noisy distance measurements. diff --git a/tutorials/example_srls.py b/tutorials/example_srls.py new file mode 100644 index 0000000..16ed860 --- /dev/null +++ b/tutorials/example_srls.py @@ -0,0 +1,28 @@ +#! /usr/bin/env python3 +# -*- coding: utf-8 -*- + +import numpy as np + +from pylocus.point_set import PointSet +from pylocus.algorithms import reconstruct_srls +from pylocus.simulation import create_noisy_edm, create_mask, create_weights +from pylocus.basics import mse, rmse + +points = PointSet(N=5, d=2) +points.set_points('random') +print("point to localize:", points.points[0, :]) +print("anchors:", points.points[1:, :]) + +std = 0.1 +edm_noisy = create_noisy_edm(points.edm, noise=std) + +mask = create_mask(points.N , method='none') +weights = create_weights(points.N, method='one') +weights = np.multiply(mask, weights) + +points_estimated = reconstruct_srls(edm_noisy, points.points, W=weights) +error = mse(points_estimated[0, :], points.points[0, :]) + +print("estimated point: {}, original point: {}, mse: {:2.2e}".format(points_estimated[0, :], + points.points[0, :], + error))