From 5b038dd141b07b8f03a519a89f3abcae0b1e9c97 Mon Sep 17 00:00:00 2001 From: ShantanuKodgirwar Date: Thu, 24 Oct 2024 20:48:19 +0200 Subject: [PATCH] moving the brownian curve implementation + associated refactorings --- geosss/sphere.py | 79 -------------------------- geosss/spherical_curve.py | 79 ++++++++++++++++++++++++++ scripts/curve_Nd.py | 14 ++--- scripts/curve_Nd_parallelized_nruns.py | 9 +-- 4 files changed, 90 insertions(+), 91 deletions(-) diff --git a/geosss/sphere.py b/geosss/sphere.py index d6b4234..9f48094 100644 --- a/geosss/sphere.py +++ b/geosss/sphere.py @@ -128,82 +128,3 @@ def rotate(theta, u=u, v=v, x=x): ) return rotate - - -def brownian_curve(n_points=100, dimension=6, step_size=0.05, seed=1234): - """ - Generate smooth points on a unit d-sphere using Brownian motion. - - Parameters: - n_points: number of points to generate - dimension: dimension of the sphere - step_size: step size for the Brownian motion, increasing it will make the points jump more - seed: random seed - """ - rng = np.random.default_rng(seed) - - # Initialize the first point on the dimension-1 -sphere - points = np.zeros((n_points, dimension)) - points[0] = radial_projection(rng.standard_normal(dimension)) - - # Generate subsequent points via Brownian motion (small random steps) - for i in range(1, n_points): - step = rng.normal(size=dimension) * step_size - new_point = points[i - 1] + step - - # Project back to the unit sphere - points[i] = radial_projection(new_point) - - return points - - -def constrained_brownian_curve(n_points=100, dimension=6, step_size=0.05, seed=1234): - """ - Generate smooth points on a unit d-sphere using constrained Brownian motion - to avoid loops and overlaps. - - Parameters: - n_points: number of points to generate - dimension: dimension of the sphere (must be at least 2) - step_size: step size for the motion - seed: random seed - """ - rng = np.random.default_rng(seed) - - # Initialize the first point on the (dimension - 1)-sphere - points = np.zeros((n_points, dimension)) - points[0] = radial_projection(rng.standard_normal(dimension)) - - # Initialize the direction of motion (tangent vector at the first point) - v = rng.standard_normal(dimension) - v -= np.dot(v, points[0]) * points[0] # Make v orthogonal to points[0] - v /= np.linalg.norm(v) # Normalize the tangent vector - - for i in range(1, n_points): - # Generate a small random perturbation orthogonal to both v and points[i - 1] - random_step = rng.standard_normal(dimension) - - # Make the random perturbation tangent to the sphere at points[i - 1] - random_step -= np.dot(random_step, points[i - 1]) * points[i - 1] - - # Make the random perturbation orthogonal to the current direction v - random_step -= np.dot(random_step, v) * v - - # Normalize the random_step so that it lies on the tangent plane - random_step /= np.linalg.norm(random_step) - - # Scale the random_step with the given step size - random_step *= step_size - - # Update the tangent vector v - v += random_step - - # Ensure the updated tangent vector v remains tangent to the sphere - v -= np.dot(v, points[i - 1]) * points[i - 1] - v /= np.linalg.norm(v) - - # Move the point along v by a small angle (step_size) - angle = step_size - points[i] = points[i - 1] * np.cos(angle) + v * np.sin(angle) - - return points diff --git a/geosss/spherical_curve.py b/geosss/spherical_curve.py index 47d65cf..1eac49e 100644 --- a/geosss/spherical_curve.py +++ b/geosss/spherical_curve.py @@ -100,3 +100,82 @@ def find_nearest(self, point): distances.append(d) points.append(y) return points[np.argmin(distances)] + + +def brownian_curve(n_points=100, dimension=6, step_size=0.05, seed=1234): + """ + Generate smooth points on a unit d-sphere using Brownian motion. + + Parameters: + n_points: number of points to generate + dimension: dimension of the sphere + step_size: step size for the Brownian motion, increasing it will make the points jump more + seed: random seed + """ + rng = np.random.default_rng(seed) + + # Initialize the first point on the dimension-1 -sphere + points = np.zeros((n_points, dimension)) + points[0] = sphere.radial_projection(rng.standard_normal(dimension)) + + # Generate subsequent points via Brownian motion (small random steps) + for i in range(1, n_points): + step = rng.normal(size=dimension) * step_size + new_point = points[i - 1] + step + + # Project back to the unit sphere + points[i] = sphere.radial_projection(new_point) + + return points + + +def constrained_brownian_curve(n_points=100, dimension=6, step_size=0.05, seed=1234): + """ + Generate smooth points on a unit d-sphere using constrained Brownian motion + to avoid loops and overlaps. + + Parameters: + n_points: number of points to generate + dimension: dimension of the sphere (must be at least 2) + step_size: step size for the motion + seed: random seed + """ + rng = np.random.default_rng(seed) + + # Initialize the first point on the (dimension - 1)-sphere + points = np.zeros((n_points, dimension)) + points[0] = sphere.radial_projection(rng.standard_normal(dimension)) + + # Initialize the direction of motion (tangent vector at the first point) + v = rng.standard_normal(dimension) + v -= np.dot(v, points[0]) * points[0] # Make v orthogonal to points[0] + v /= np.linalg.norm(v) # Normalize the tangent vector + + for i in range(1, n_points): + # Generate a small random perturbation orthogonal to both v and points[i - 1] + random_step = rng.standard_normal(dimension) + + # Make the random perturbation tangent to the sphere at points[i - 1] + random_step -= np.dot(random_step, points[i - 1]) * points[i - 1] + + # Make the random perturbation orthogonal to the current direction v + random_step -= np.dot(random_step, v) * v + + # Normalize the random_step so that it lies on the tangent plane + random_step /= np.linalg.norm(random_step) + + # Scale the random_step with the given step size + random_step *= step_size + + # Update the tangent vector v + v += random_step + + # Ensure the updated tangent vector v remains tangent to the sphere + v -= np.dot(v, points[i - 1]) * points[i - 1] + v /= np.linalg.norm(v) + + # Move the point along v by a small angle (step_size) + angle = step_size + points[i] = points[i - 1] * np.cos(angle) + v * np.sin(angle) + + return points diff --git a/scripts/curve_Nd.py b/scripts/curve_Nd.py index 1db9246..ae1cfd3 100644 --- a/scripts/curve_Nd.py +++ b/scripts/curve_Nd.py @@ -10,7 +10,7 @@ import geosss as gs from geosss.distributions import CurvedVonMisesFisher, Distribution -from geosss.spherical_curve import SlerpCurve, SphericalCurve +from geosss.spherical_curve import SlerpCurve, SphericalCurve, brownian_curve mpl.rcParams["mathtext.fontset"] = "cm" # Use Computer Modern font @@ -307,7 +307,7 @@ def argparser(): n_runs = args.n_runs # sampler runs (default: 1), for ess computations `n_runs=10` # optional controls - brownian_curve = True # brownian curve or curve with fixed knots + is_brownian_curve = True # brownian curve or curve with fixed knots reprod_switch = True # seeds samplers for reproducibility show_plots = True # show the plots savefig = True # save the plots @@ -319,7 +319,7 @@ def argparser(): setup_logging(savedir, kappa) # Define curve on the sphere - if not brownian_curve: + if not is_brownian_curve: knots = np.array( [ [-0.25882694, 0.95006168, 0.17433133], @@ -339,18 +339,16 @@ def argparser(): if n_dim > knots.shape[1]: knots = np.pad(knots, ((0, 0), (n_dim - knots.shape[1], 0))) - curve = SlerpCurve(knots) - else: # generates a smooth curve on the sphere with brownian motion - knots = gs.sphere.brownian_curve( + knots = brownian_curve( n_points=10, dimension=n_dim, step_size=0.5, # larger step size will result in more spread out points seed=4562, ) - curve = SlerpCurve(knots) - # curve = SlerpCurve.random_curve(n_knots=100, seed=4562, dimension=n_dim) + + curve = SlerpCurve(knots) # Initialize based on dimensionality initial = ( diff --git a/scripts/curve_Nd_parallelized_nruns.py b/scripts/curve_Nd_parallelized_nruns.py index 452343b..9c9cd49 100644 --- a/scripts/curve_Nd_parallelized_nruns.py +++ b/scripts/curve_Nd_parallelized_nruns.py @@ -8,7 +8,7 @@ import geosss as gs from geosss.distributions import CurvedVonMisesFisher, Distribution -from geosss.spherical_curve import SlerpCurve +from geosss.spherical_curve import SlerpCurve, brownian_curve def setup_logging( @@ -298,12 +298,12 @@ def main(): assert n_dim >= 3, msg # optional controls - brownian_curve = True # brownian curve or curve with fixed knots + is_brownian_curve = True # brownian curve or curve with fixed knots reprod_switch = True # seeds samplers for reproducibility rerun_if_samples_exists = False # rerun even if samples file exists # creating a target as a curve on the sphere - if not brownian_curve: + if not is_brownian_curve: knots = np.array( [ [-0.25882694, 0.95006168, 0.17433133], @@ -325,13 +325,14 @@ def main(): else: # generates a smooth curve on the sphere with brownian motion - knots = gs.sphere.brownian_curve( + knots = brownian_curve( n_points=10, dimension=n_dim, step_size=0.5, # larger step size will result in more spread out points seed=4562, ) + logging.info(f"Target curve: {knots}") # defining the curve on the sphere curve = SlerpCurve(knots)