From 7c3749ed786d002917bce2b4c235a67b4e5d5fb4 Mon Sep 17 00:00:00 2001 From: Chris Granade Date: Tue, 7 Mar 2017 17:38:53 +1100 Subject: [PATCH 1/6] Initial implementation of convex distance. --- src/qinfer/__init__.py | 1 + src/qinfer/geometry.py | 93 ++++++++++++++++++++++++++++++++++++++++++ src/qinfer/utils.py | 36 ++++++++++++++++ 3 files changed, 130 insertions(+) create mode 100644 src/qinfer/geometry.py diff --git a/src/qinfer/__init__.py b/src/qinfer/__init__.py index 7413ecb..4e29ab6 100644 --- a/src/qinfer/__init__.py +++ b/src/qinfer/__init__.py @@ -31,6 +31,7 @@ from qinfer._exceptions import * +from qinfer.geometry import * from qinfer.gpu_models import * from qinfer.perf_testing import * from qinfer.expdesign import * diff --git a/src/qinfer/geometry.py b/src/qinfer/geometry.py new file mode 100644 index 0000000..32e0587 --- /dev/null +++ b/src/qinfer/geometry.py @@ -0,0 +1,93 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- +## +# geometry.py: Computational geometry utilities. +## +# © 2016 Chris Ferrie (csferrie@gmail.com) and +# Christopher E. Granade (cgranade@gmail.com) +# +# This file is a part of the Qinfer project. +# Licensed under the AGPL version 3. +## +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Affero General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with this program. If not, see . +## + +## FEATURES ################################################################### + +from __future__ import absolute_import +from __future__ import division + +## EXPORTS ################################################################### + +__all__ = [ + # TODO: + # 'simple_est_rabi' +] + +## IMPORTS ################################################################### + +import numpy as np + +from qinfer.utils import requires_optional_module + +try: + import cvxopt + import cvxopt.solvers +except: + cvxopt = None + +## FUNCTIONS ################################################################# + +@requires_optional_module('cvxopt') +def convex_distance(test_point, convex_points): + r""" + Returns the distance :math:`d(\vec{x}, S) = + \min_{\vec{z} \in \operatorname{conv}(S)} \| \vec{x} - \vec{z} \|^2` from + a point :math:`\vec{x}` to the convex closure of a set of points + :math:`S`. + """ + # TODO: document this much better. + + n_points, dimension = convex_points.shape + + # Start by centering the convex region of interest. + # NB: don't use in-place (-=) here, since we want a copy. + convex_points = convex_points - test_point + + # Construct the quadratic program parameters P and q required by CVXOPT. + # Since we want min x^T P x, with x representing the weight of each point + # in our set of a convex combination, we need to construct the outer + # product + # P_{ik} = sum_j y_ij y_jk + # where y_ij is the jth component of the ith point in our set. + P = cvxopt.matrix(np.dot(convex_points, convex_points.T)) + # We have no linear objective, so we set q = 0. + # Similarly, q is a vector of all ones representing the sum + # over points in the set. + q = cvxopt.matrix(np.zeros((n_points, 1))) + # To enfoce x_i ≥ 0, we need G = -\mathbb{1} and h = 0, + # such that CVXOPT's constraint Gx ≥ h is correct. + G = cvxopt.matrix(-np.eye(n_points)) + h = cvxopt.matrix(0.0, (n_points, 1)) + + # We do, however, have a linear constraint A x = 1. + A = cvxopt.matrix(1.0, (1, n_points)) + b = cvxopt.matrix(1.0) + + # Now we solve with cvxopt. + solution = cvxopt.solvers.qp(P, q, G, h, A, b) + + # Return the obtained primal objective, noting that CVXOPT finds + # min_x 1/2 x^T P X, so that we must multiply by 2. + return 2 * solution['primal objective'] diff --git a/src/qinfer/utils.py b/src/qinfer/utils.py index 8cac929..2cf9da9 100644 --- a/src/qinfer/utils.py +++ b/src/qinfer/utils.py @@ -32,6 +32,8 @@ ## IMPORTS #################################################################### from builtins import range +from importlib import import_module +from distutils.version import LooseVersion import warnings @@ -48,6 +50,40 @@ ## FUNCTIONS ################################################################## +def get_optional_module(module_name, required_version=None): + """ + Attempts to load an optional module, silently returning + ``None`` if the module is missing or is not of the required + version. + """ + try: + module = import_module(module_name) + if ( + required_version is not None and + module.__version__ < LooseVersion(required_version) + ): + return None + except ImportError, AttributeError: + return None + + return module + + +def requires_optional_module(module_name, required_version=None): + def decorator(fn): + if get_optional_module(module_name, required_version) is not None: + return fn + else: + raise ImportError( + "Function {} requires version {} of the optional module " + "{}, which was not found.".format( + fn.__name__, required_version, module_name + ) + ) + + return decorator + +# FIXME: generalize to use above *_optional_module functions. def get_qutip_module(required_version='3.2'): """ Attempts to return the qutip module, but From c3006c32191e91020b86bf1bac5fc46dd9a9fa56 Mon Sep 17 00:00:00 2001 From: Chris Granade Date: Tue, 7 Mar 2017 21:45:14 +1100 Subject: [PATCH 2/6] Improvements to optional module loading. --- src/qinfer/utils.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/qinfer/utils.py b/src/qinfer/utils.py index 2cf9da9..dd1cb5c 100644 --- a/src/qinfer/utils.py +++ b/src/qinfer/utils.py @@ -33,6 +33,7 @@ from builtins import range from importlib import import_module +from functools import wraps from distutils.version import LooseVersion import warnings @@ -63,7 +64,7 @@ def get_optional_module(module_name, required_version=None): module.__version__ < LooseVersion(required_version) ): return None - except ImportError, AttributeError: + except (ImportError, AttributeError): return None return module @@ -74,13 +75,23 @@ def decorator(fn): if get_optional_module(module_name, required_version) is not None: return fn else: - raise ImportError( + msg = ( "Function {} requires version {} of the optional module " "{}, which was not found.".format( fn.__name__, required_version, module_name ) + if required_version is not None else + "Function {} requires the optional module {}, which was " + "not found.".format( + fn.__name__, module_name + ) ) + @wraps(fn) + def dummy_fn(*args, **kwwargs): + raise ImportError(msg) + return dummy_fn + return decorator # FIXME: generalize to use above *_optional_module functions. From 31b23157adffb519a366bbfeac0ee27562c3dbad Mon Sep 17 00:00:00 2001 From: Chris Granade Date: Wed, 8 Mar 2017 10:47:51 +1100 Subject: [PATCH 3/6] Added tests of new geometry module. --- src/qinfer/tests/test_geometry.py | 72 +++++++++++++++++++++++++++++++ src/qinfer/utils.py | 19 +++++--- 2 files changed, 86 insertions(+), 5 deletions(-) create mode 100644 src/qinfer/tests/test_geometry.py diff --git a/src/qinfer/tests/test_geometry.py b/src/qinfer/tests/test_geometry.py new file mode 100644 index 0000000..e171928 --- /dev/null +++ b/src/qinfer/tests/test_geometry.py @@ -0,0 +1,72 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- +## +# test_geometry.py: Checks correctness of computational geometry functionality. +## +# © 2014 Chris Ferrie (csferrie@gmail.com) and +# Christopher E. Granade (cgranade@gmail.com) +# +# This file is a part of the Qinfer project. +# Licensed under the AGPL version 3. +## +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Affero General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with this program. If not, see . +## + +## FEATURES ################################################################### + +from __future__ import absolute_import +from __future__ import division # Ensures that a/b is always a float. + +## IMPORTS #################################################################### + +import numpy as np +from numpy.testing import assert_almost_equal + +from qinfer.tests.base_test import DerandomizedTestCase + +from qinfer.geometry import convex_distance +from qinfer.utils import assert_sigfigs_equal, requires_optional_module + +## TESTS ##################################################################### + +class TestConvexGeometry(DerandomizedTestCase): + + @requires_optional_module("cvxopt", if_missing="skip") + def test_convex_distance_known_cases(self): + """ + Convex geometry: Checks convex_distance for known cases. + """ + assert_almost_equal( + convex_distance(np.ones((3,)), np.eye(3)), + np.linalg.norm(np.ones(3,), ord=1) + ) + + # TODO: add cases. + + @requires_optional_module("cvxopt", if_missing="skip") + def test_convex_distance_interior_point(self, n_points=500, dimension=5): + """ + Convex geometry: Checks that convex_distance is zero for interior pts. + """ + S = np.random.random((n_points, dimension)) + alpha = np.random.random((n_points)) + alpha /= alpha.sum() + + z = np.dot(alpha, S) + + assert_almost_equal( + convex_distance(z, S), + 0 + ) + diff --git a/src/qinfer/utils.py b/src/qinfer/utils.py index dd1cb5c..2628ea4 100644 --- a/src/qinfer/utils.py +++ b/src/qinfer/utils.py @@ -37,6 +37,7 @@ from distutils.version import LooseVersion import warnings +import unittest import numpy as np import numpy.linalg as la @@ -70,7 +71,9 @@ def get_optional_module(module_name, required_version=None): return module -def requires_optional_module(module_name, required_version=None): +def requires_optional_module(module_name, required_version=None, if_missing='raise'): + # TODO: document + def decorator(fn): if get_optional_module(module_name, required_version) is not None: return fn @@ -87,10 +90,15 @@ def decorator(fn): ) ) - @wraps(fn) - def dummy_fn(*args, **kwwargs): - raise ImportError(msg) - return dummy_fn + if if_missing == "raise": + @wraps(fn) + def dummy_fn(*args, **kwwargs): + raise ImportError(msg) + + return dummy_fn + + elif if_missing == "skip": + return unittest.skip(msg) return decorator @@ -316,6 +324,7 @@ def particle_covariance_mtx(weights,locations): return cov +# FIXME: move ellipsoid_volume, mvee, and in_ellipsoid to geometry.py. def ellipsoid_volume(A=None, invA=None): """ From 6dae269b5880e9455e43296a4d8799e19e99691c Mon Sep 17 00:00:00 2001 From: Chris Granade Date: Wed, 8 Mar 2017 10:58:39 +1100 Subject: [PATCH 4/6] Fixed one bug detected by new test cases. --- src/qinfer/geometry.py | 11 +++++++---- src/qinfer/tests/test_geometry.py | 20 ++++++++++++-------- 2 files changed, 19 insertions(+), 12 deletions(-) diff --git a/src/qinfer/geometry.py b/src/qinfer/geometry.py index 32e0587..801a8f5 100644 --- a/src/qinfer/geometry.py +++ b/src/qinfer/geometry.py @@ -31,8 +31,7 @@ ## EXPORTS ################################################################### __all__ = [ - # TODO: - # 'simple_est_rabi' + # TODO ] ## IMPORTS ################################################################### @@ -89,5 +88,9 @@ def convex_distance(test_point, convex_points): solution = cvxopt.solvers.qp(P, q, G, h, A, b) # Return the obtained primal objective, noting that CVXOPT finds - # min_x 1/2 x^T P X, so that we must multiply by 2. - return 2 * solution['primal objective'] + # min_x 1/2 x^T P X, so that we must multiply by 2, then take the + # square root. Taking the square root, we have to be a little + # careful, since finite tolerances can cause the objective to + # be negative, even though we constrained it to be at least zero. + # Thus, we take the max of the objective and zero. + return np.sqrt(2 * max(0, solution['primal objective'])) diff --git a/src/qinfer/tests/test_geometry.py b/src/qinfer/tests/test_geometry.py index e171928..7fc684d 100644 --- a/src/qinfer/tests/test_geometry.py +++ b/src/qinfer/tests/test_geometry.py @@ -36,20 +36,23 @@ from qinfer.tests.base_test import DerandomizedTestCase from qinfer.geometry import convex_distance -from qinfer.utils import assert_sigfigs_equal, requires_optional_module +from qinfer.utils import ( + assert_sigfigs_equal, requires_optional_module, +) +from qinfer.perf_testing import numpy_err_policy ## TESTS ##################################################################### class TestConvexGeometry(DerandomizedTestCase): @requires_optional_module("cvxopt", if_missing="skip") - def test_convex_distance_known_cases(self): + def test_convex_distance_known_cases(self, dimension=3): """ Convex geometry: Checks convex_distance for known cases. """ assert_almost_equal( - convex_distance(np.ones((3,)), np.eye(3)), - np.linalg.norm(np.ones(3,), ord=1) + convex_distance(np.ones((dimension,)), np.eye(dimension)), + np.linalg.norm(((dimension - 1) / dimension) * np.ones(dimension,), ord=2) ) # TODO: add cases. @@ -65,8 +68,9 @@ def test_convex_distance_interior_point(self, n_points=500, dimension=5): z = np.dot(alpha, S) - assert_almost_equal( - convex_distance(z, S), - 0 - ) + with numpy_err_policy(invalid='raise'): + assert_almost_equal( + convex_distance(z, S), + 0 + ) From 1c6be1e9928a7979128b5f78d2d7bbb14d4fbdc0 Mon Sep 17 00:00:00 2001 From: Chris Granade Date: Wed, 8 Mar 2017 15:59:20 +1100 Subject: [PATCH 5/6] =?UTF-8?q?Initial=20version=20of=20=CE=B5-approx=20co?= =?UTF-8?q?nvex=20hull?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/qinfer/geometry.py | 189 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 186 insertions(+), 3 deletions(-) diff --git a/src/qinfer/geometry.py b/src/qinfer/geometry.py index 801a8f5..b9e24b2 100644 --- a/src/qinfer/geometry.py +++ b/src/qinfer/geometry.py @@ -31,12 +31,14 @@ ## EXPORTS ################################################################### __all__ = [ - # TODO + 'convex_distance', + 'approximate_convex_hull' ] ## IMPORTS ################################################################### import numpy as np +from scipy.spatial.distance import cdist from qinfer.utils import requires_optional_module @@ -48,6 +50,8 @@ ## FUNCTIONS ################################################################# +# FIXME: add outgoing citations. + @requires_optional_module('cvxopt') def convex_distance(test_point, convex_points): r""" @@ -56,7 +60,7 @@ def convex_distance(test_point, convex_points): a point :math:`\vec{x}` to the convex closure of a set of points :math:`S`. """ - # TODO: document this much better. + # TODO: document this much better, esp. params/returns. n_points, dimension = convex_points.shape @@ -85,7 +89,9 @@ def convex_distance(test_point, convex_points): b = cvxopt.matrix(1.0) # Now we solve with cvxopt. - solution = cvxopt.solvers.qp(P, q, G, h, A, b) + solution = cvxopt.solvers.qp(P, q, G, h, A, b, + options={'show_progress': False} + ) # Return the obtained primal objective, noting that CVXOPT finds # min_x 1/2 x^T P X, so that we must multiply by 2, then take the @@ -94,3 +100,180 @@ def convex_distance(test_point, convex_points): # be negative, even though we constrained it to be at least zero. # Thus, we take the max of the objective and zero. return np.sqrt(2 * max(0, solution['primal objective'])) + +def rand_argmin(arr, tol=1e-7): + minimum = np.min(arr) + return np.random.choice(np.nonzero(arr <= minimum + tol)[0]) + +def expand_hull_directed_search(partial_hull, candidate_points, interior_thresh=1e-7): + r""" + Uses the Sartipizadeh and Vincent 2016 algorithm to expand + an approximate convex hull by finding a point :math:`\vec{x}` + solving + + .. math:: + + \vec{x} = \operatorname{arg\,min}_{\vec{x} \in S \setminus E} + \max_{\vec{z} \in S \setminus E} + d(\vec{z}, E \cup \{\vec{x}\}) + + where :math:`E` is the set of points comprising the approximate + hull so far, and where :math:`S` is the set of points that we + wish to approximate the convex hull of. That is, this function + finds a point to add to an approximate hull that minimizes + the worst-case distance to any point not currently in the hull. + """ + # TODO: document this much better, esp. params/returns. + + # returns: + # index of point to add to the hull + # current error of the approximate hull + # masked matrix of known convex distances + + # Our strategy will be to build up a matrix D_{i, j} + # defined as: + # + # Dᵢⱼ ≔ d(xᵢ, E ∪ {xⱼ}) + # + # where the set of candidate points S \ E = {xᵢ} is represented + # as the array candidate_points. + # + # NB: our notation is different from VC16, since they reuse + # the letter E. + + n_partial, dimension = partial_hull.shape + n_candidates, dimension_check = candidate_points.shape + assert dimension == dimension_check + + # Start by allocating arrays to represent how far into + # each maximization we've considered for the minimization over + # candidates. + n_distances_considered = np.ones((n_candidates,), dtype=int) + + # We will also need to remember the list of points interior + # to each candidate so that they can be removed later. + interior_points_found = [[] for idx in range(n_candidates)] + + # We then define how to compute each element. + closest_miss = [np.inf] + def distance(idx_max, idx_min): + if idx_max == idx_min: + this_distance = 0 + else: + this_distance = convex_distance( + candidate_points[idx_max], + np.concatenate([ + partial_hull, + candidate_points[idx_min, np.newaxis] + ], axis=0) + ) + if this_distance < closest_miss[0]: + closest_miss[0] = this_distance + + if this_distance < interior_thresh: + interior_points_found[idx_min].append(idx_max) + return this_distance + + # Next, we allocate an array for the worst maxima we've seen for + # each candidate. + worst_distances = np.array([ + distance(0, idx_min) + for idx_min in range(n_candidates) + ]) + + # Using this, we can find our best candidate so far + idx_best_candidate = rand_argmin(worst_distances) + + # ITERATION: Look at each "row" in turn until we get to the bottom. At each step, + # evalaute the next element of the "column" with the best (smallest) + # maximum. + while n_distances_considered[idx_best_candidate] < n_candidates: + + next_distance = distance(n_distances_considered[idx_best_candidate], idx_best_candidate) + n_distances_considered[idx_best_candidate] += 1 + + if next_distance > worst_distances[idx_best_candidate]: + worst_distances[idx_best_candidate] = next_distance + + idx_best_candidate = rand_argmin(worst_distances) + + # COMPLETION: Return the best element we've found so far, along with the value + # at that element, and the indices of interior points for the newly + # expanded hull. + approximation_error = worst_distances[idx_best_candidate] + idxs_interior_points = interior_points_found[idx_best_candidate] + + print(closest_miss) + + return idx_best_candidate, approximation_error, idxs_interior_points + + +def approximate_convex_hull(points, max_n_vertices, desired_approximation_error, + interior_thresh=1e-7, + seed_method='extreme' + ): + + n_points, dimension = points.shape + max_n_vertices = min(n_points - 1, max_n_vertices) + + # FIXME: we should shuffle. + + # Start by picking a point a seed point to grow the hull from. + if seed_method == 'extreme': + idx_dimension = np.random.randint(dimension) + sign = np.random.choice([-1, 1]) + idx_seed_points = [ + (sign * points[:, idx_dimension]).argmin() + ] + elif seed_method == 'centroid': + centroid = np.mean(points, axis=0) + centroid_distances = cdist(points, centroid[None, :])[:, 0] + idx_seed_points = [centroid_distances.argmin(), centroid_distances.argmax()] + + + # We maintain masks to indicate which points are in the approximate + # hull, and which masks are candidates (haven't been eliminated as + # interior). + hull_mask = np.zeros((n_points,), dtype=bool) + hull_mask[idx_seed_points] = True + + candidate_mask = ~hull_mask.copy() + + # Iterate until we hit max_n_vertices or epsilon_desired, expanding + # the hull each time and eliminating the interior points. + while True: + partial_hull = points[hull_mask] + candidate_points = points[candidate_mask] + idxs_partial_hull = np.nonzero(hull_mask)[0] + + idx_best_candidate, approximation_error, idxs_interior_points = \ + expand_hull_directed_search(partial_hull, candidate_points, interior_thresh=interior_thresh) + + # Add the new best candidate to the approximate hull, + # and eliminate any candidate points interior to the new point. + idxs_candidates = np.nonzero(candidate_mask)[0] + hull_mask[idxs_candidates[idx_best_candidate]] = True + candidate_mask[idxs_candidates[idxs_interior_points]] = False + + # Eliminate anything from the hull_mask + # that's interior to the new point. + n_hull_removed = 0 + reduced_hull_mask = hull_mask.copy() + for idx_previous_hull_point in idxs_partial_hull: + reduced_hull_mask[idx_previous_hull_point] = False + if convex_distance(points[idx_previous_hull_point], points[reduced_hull_mask]) <= interior_thresh: + # Remove this as a hull point. + n_hull_removed += 1 + hull_mask[idx_previous_hull_point] = False + else: + # The point wasn't interior, so add it back in. + reduced_hull_mask[idx_previous_hull_point] = True + + print("Eliminated {} interior points.".format(len(idxs_interior_points) - 1 + n_hull_removed)) + + # If we hit either of our stopping criteria, return now. + if ( + approximation_error <= desired_approximation_error or + np.sum(hull_mask) >= max_n_vertices + ): + return np.nonzero(hull_mask)[0], approximation_error From 9aa4394034ce360a54b1a4cba8fb2b4648130dfa Mon Sep 17 00:00:00 2001 From: Chris Granade Date: Mon, 20 Mar 2017 15:06:51 +1100 Subject: [PATCH 6/6] WIP on new heuristics to reduce initial pt cloud size. --- src/qinfer/geometry.py | 38 ++++++++++++++++++++++++++++++++------ 1 file changed, 32 insertions(+), 6 deletions(-) diff --git a/src/qinfer/geometry.py b/src/qinfer/geometry.py index b9e24b2..8ee1265 100644 --- a/src/qinfer/geometry.py +++ b/src/qinfer/geometry.py @@ -40,7 +40,7 @@ import numpy as np from scipy.spatial.distance import cdist -from qinfer.utils import requires_optional_module +from qinfer.utils import requires_optional_module, uniquify try: import cvxopt @@ -220,11 +220,15 @@ def approximate_convex_hull(points, max_n_vertices, desired_approximation_error, # Start by picking a point a seed point to grow the hull from. if seed_method == 'extreme': - idx_dimension = np.random.randint(dimension) - sign = np.random.choice([-1, 1]) - idx_seed_points = [ + # http://www.sciencedirect.com/science/article/pii/S2405896315009866 + # idx_dimension = np.random.randint(dimension) + # sign = np.random.choice([-1, 1]) + idx_seed_points = uniqify([ (sign * points[:, idx_dimension]).argmin() - ] + for idx_dimension in range(dimension) + for sign in [-1, 1] + ]) + elif seed_method == 'centroid': centroid = np.mean(points, axis=0) centroid_distances = cdist(points, centroid[None, :])[:, 0] @@ -239,6 +243,9 @@ def approximate_convex_hull(points, max_n_vertices, desired_approximation_error, candidate_mask = ~hull_mask.copy() + # Before we iterate, remove candidates that are interior to the seed hull. + # TODO + # Iterate until we hit max_n_vertices or epsilon_desired, expanding # the hull each time and eliminating the interior points. while True: @@ -248,7 +255,7 @@ def approximate_convex_hull(points, max_n_vertices, desired_approximation_error, idx_best_candidate, approximation_error, idxs_interior_points = \ expand_hull_directed_search(partial_hull, candidate_points, interior_thresh=interior_thresh) - + # Add the new best candidate to the approximate hull, # and eliminate any candidate points interior to the new point. idxs_candidates = np.nonzero(candidate_mask)[0] @@ -277,3 +284,22 @@ def approximate_convex_hull(points, max_n_vertices, desired_approximation_error, np.sum(hull_mask) >= max_n_vertices ): return np.nonzero(hull_mask)[0], approximation_error + +# FIXME: remove this. +if __name__ == "__main__": + + dims = np.arange(2, 11) + times = np.empty_like(dims, dtype=float) + n_points = 100 + + from qinfer.perf_testing import timing + + for idx_dim, dim in enumerate(dims): + print("======= DIM {} =======".format(dim)) + points = np.random.randn(n_points, dim) + points /= np.linalg.norm(points, axis=0) + points *= np.random.random((n_points,))[:, None] + with timing() as t: + hull_pts, eps = approximate_convex_hull(points, n_points - 1, 0, seed_method='extreme') + times[idx_dim] = t.delta_t + print(t)