From ffa33e0cb28a712f4cf7f7bf37877bfa23356809 Mon Sep 17 00:00:00 2001 From: Alex Fox Date: Fri, 10 Oct 2025 21:54:41 -0600 Subject: [PATCH 1/5] fixed minor bug --- PyOCN/c_src/flowgrid.c | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/PyOCN/c_src/flowgrid.c b/PyOCN/c_src/flowgrid.c index 3e93213..627935e 100644 --- a/PyOCN/c_src/flowgrid.c +++ b/PyOCN/c_src/flowgrid.c @@ -204,7 +204,10 @@ FlowGrid *fg_copy(FlowGrid *G){ Status fg_destroy(FlowGrid *G){ if (G != NULL){ - if (G->vertices != NULL) free(G->vertices); G->vertices = NULL; + if (G->vertices != NULL){ + free(G->vertices); + G->vertices = NULL; + } free(G); } return SUCCESS; From 713f7f0d69667b02069e82b4f350affa826b76b2 Mon Sep 17 00:00:00 2001 From: Alex Fox Date: Sat, 11 Oct 2025 18:04:12 -0600 Subject: [PATCH 2/5] fixed a bug where temperature was incorrectly computed in single-outlet networks --- PyOCN/c_src/ocn.c | 13 ++++++------- sandbox.py | 11 ++++++++--- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/PyOCN/c_src/ocn.c b/PyOCN/c_src/ocn.c index 34798b3..d906b1c 100644 --- a/PyOCN/c_src/ocn.c +++ b/PyOCN/c_src/ocn.c @@ -22,11 +22,12 @@ * @param temperature The current temperature. * @return true if the new state is accepted, false otherwise. */ -static inline bool simulate_annealing(double energy_new, double energy_old, double inv_temperature, rng_state_t *rng){ - const double delta_energy = energy_new - energy_old; +bool simulate_annealing(double energy_new, double energy_old, double temperature, rng_state_t *rng){ + double accpt_prob = rng_uniformdouble(rng); + double delta_energy = energy_new - energy_old; if (delta_energy <= 0.0) return true; // Always accept improvements - const double p = exp(-delta_energy * inv_temperature); - return rng_uniformdouble(rng) < p; + double p = exp(-delta_energy / temperature); + return accpt_prob < p; } /** @@ -97,8 +98,6 @@ Status ocn_single_erosion_event(FlowGrid *G, double gamma, double temperature){ CartPair dims = G->dims; linidx_t nverts = (linidx_t)dims.row * (linidx_t)dims.col; - double inv_temperature = 1.0 / temperature; - double energy_old, energy_new; energy_old = G->energy; @@ -172,7 +171,7 @@ Status ocn_single_erosion_event(FlowGrid *G, double gamma, double temperature){ update_drained_area(G, -da_inc, a_down_old); // remove drainage from old path update_drained_area(G, da_inc, a_down_new); // add drainage to new path energy_new = ocn_compute_energy(G, gamma); // recompute energy from scratch - if (simulate_annealing(energy_new, energy_old, inv_temperature, &G->rng)){ + if (simulate_annealing(energy_new, energy_old, temperature, &G->rng)){ G->energy = energy_new; return SUCCESS; } diff --git a/sandbox.py b/sandbox.py index 10053f5..86baa45 100644 --- a/sandbox.py +++ b/sandbox.py @@ -1,10 +1,15 @@ -from timeit import timeit +import matplotlib.pyplot as plt import PyOCN as po ocn = po.OCN.from_net_type( net_type="V", - dims=(74, 74), + dims=(100, 100), + wrap=True, random_state=8472, ) -print(timeit(ocn.fit, number=10)/10) \ No newline at end of file +ocn.fit(pbar=True, n_iterations=100**2*100, constant_phase=0.25, cooling_rate=1) + +plt.plot(ocn.history[:, 0], ocn.history[:, 1]) +plt.plot(ocn.history[:, 0], ocn.history[:, 2]) +plt.show() \ No newline at end of file From ee4358be1b674bc1b4bb5d0d45ff4b62c8473b7d Mon Sep 17 00:00:00 2001 From: Alex Fox Date: Sat, 11 Oct 2025 18:28:09 -0600 Subject: [PATCH 3/5] added a test to confirm that convergence occurs with a greedy cooling schedule --- sandbox.py | 19 ++++++++++++++----- tests/test_basic.py | 21 +++++++++++++++++++++ 2 files changed, 35 insertions(+), 5 deletions(-) diff --git a/sandbox.py b/sandbox.py index 86baa45..1f621de 100644 --- a/sandbox.py +++ b/sandbox.py @@ -1,15 +1,24 @@ import matplotlib.pyplot as plt import PyOCN as po +import numpy as np ocn = po.OCN.from_net_type( - net_type="V", - dims=(100, 100), + net_type="E", + dims=(16, 16), wrap=True, - random_state=8472, + random_state=84712, ) -ocn.fit(pbar=True, n_iterations=100**2*100, constant_phase=0.25, cooling_rate=1) +ocn.fit_custom_cooling(lambda t: np.ones_like(t)*1e-7, pbar=True, n_iterations=16*16*100, max_iterations_per_loop=1) -plt.plot(ocn.history[:, 0], ocn.history[:, 1]) + + +print(ocn.history.shape) +print(np.max(np.diff(ocn.history[:, 1]))) +print(np.quantile(np.diff(ocn.history[:, 1]), 0.999)) + +plt.plot(ocn.history[:-1, 0], -np.diff(ocn.history[:, 1])) plt.plot(ocn.history[:, 0], ocn.history[:, 2]) +# plt.xscale("log") +plt.yscale("log") plt.show() \ No newline at end of file diff --git a/tests/test_basic.py b/tests/test_basic.py index 7f498b0..603fbbc 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -252,6 +252,27 @@ def test_error_handling(self): with self.assertRaises((ValueError, TypeError)): po.OCN.from_net_type("V", dims=(5,)) # Should need 2D dims + def test_greedy_optimization(self): + ocn = po.OCN.from_net_type( + net_type="V", + dims=(16, 16), + wrap=True, + random_state=8472, + ) + # use a temperature of 1e-7 since we were bitten before by accidentally using 1/temperature. + # Using temperature=0 could lead to infinities, which might behave unexpectedly. + # this way, if our mass is wrong and blows up somewhere, we know that we will catch it. + ocn.fit_custom_cooling(lambda t: np.ones_like(t)*1e-7, pbar=False, n_iterations=16**2*100, max_iterations_per_loop=1) + self.assertLessEqual(np.quantile(np.diff(ocn.history[:, 1]), 0.999) - 1e-7, 0, "Energy did not decrease monotonically.") + + ocn = po.OCN.from_net_type( + net_type="E", + dims=(16, 16), + wrap=True, + random_state=8472, + ) + ocn.fit_custom_cooling(lambda t: np.ones_like(t)*1e-7, pbar=False, n_iterations=16**2*100, max_iterations_per_loop=1) + self.assertLessEqual(np.quantile(np.diff(ocn.history[:, 1]), 0.999) - 1e-7, 0, "Energy did not decrease monotonically.") if __name__ == "__main__": unittest.main() \ No newline at end of file From 3cee192e9fd5c3bbb2d39f647e0547242ebf2fee Mon Sep 17 00:00:00 2001 From: Alex Fox Date: Sat, 11 Oct 2025 18:29:15 -0600 Subject: [PATCH 4/5] update version --- CHANGELOG.md | 5 ++++- docs/conf.py | 2 +- pyproject.toml | 2 +- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 71f304f..5f08d45 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,4 +16,7 @@ - fixed a bug preventing early exit in some situations during optimization when the convergence criterion is met # 1.3.20251010 -- fixed critical bug preventing export to xarray datasets \ No newline at end of file +- fixed critical bug preventing export to xarray datasets + +# 1.3.20251011 +- fixed a bug where optimization did not work for single-outlet networks \ No newline at end of file diff --git a/docs/conf.py b/docs/conf.py index 1730f3b..0457131 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -14,7 +14,7 @@ project = 'PyOCN' copyright = '2025, Alexander S. Fox' author = 'Alexander S. Fox' -release = '1.3.20251010' +release = '1.3.20251011' # -- General configuration --------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration diff --git a/pyproject.toml b/pyproject.toml index acceae7..2983a9a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "PyOCN" -version = "1.3.20251010" +version = "1.3.20251011" description = "Optimal Channel Networks (OCN) in Python with a C core" readme = "readme.md" requires-python = ">=3.10" From b2183e8ac7903ee1f4e340a024abb765460c003c Mon Sep 17 00:00:00 2001 From: Alex Fox Date: Sat, 11 Oct 2025 18:32:43 -0600 Subject: [PATCH 5/5] replaced static/const/inline for annealing --- PyOCN/c_src/ocn.c | 8 ++++---- sandbox.py | 13 +++++++------ 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/PyOCN/c_src/ocn.c b/PyOCN/c_src/ocn.c index d906b1c..8e0d75b 100644 --- a/PyOCN/c_src/ocn.c +++ b/PyOCN/c_src/ocn.c @@ -22,11 +22,11 @@ * @param temperature The current temperature. * @return true if the new state is accepted, false otherwise. */ -bool simulate_annealing(double energy_new, double energy_old, double temperature, rng_state_t *rng){ - double accpt_prob = rng_uniformdouble(rng); - double delta_energy = energy_new - energy_old; +static inline bool simulate_annealing(double energy_new, double energy_old, double temperature, rng_state_t *rng){ + const double accpt_prob = rng_uniformdouble(rng); + const double delta_energy = energy_new - energy_old; if (delta_energy <= 0.0) return true; // Always accept improvements - double p = exp(-delta_energy / temperature); + const double p = exp(-delta_energy / temperature); return accpt_prob < p; } diff --git a/sandbox.py b/sandbox.py index 1f621de..92bb830 100644 --- a/sandbox.py +++ b/sandbox.py @@ -9,16 +9,17 @@ wrap=True, random_state=84712, ) -ocn.fit_custom_cooling(lambda t: np.ones_like(t)*1e-7, pbar=True, n_iterations=16*16*100, max_iterations_per_loop=1) +ocn.fit_custom_cooling(lambda t: np.ones_like(t)*1e-1, pbar=True, n_iterations=16*16*100, max_iterations_per_loop=1) +ocn.fit(max_iterations_per_loop=1) -print(ocn.history.shape) -print(np.max(np.diff(ocn.history[:, 1]))) -print(np.quantile(np.diff(ocn.history[:, 1]), 0.999)) +# print(ocn.history.shape) +# print(np.max(np.diff(ocn.history[:, 1]))) +# print(np.quantile(np.diff(ocn.history[:, 1]), 0.999)) -plt.plot(ocn.history[:-1, 0], -np.diff(ocn.history[:, 1])) +plt.plot(ocn.history[:, 0], ocn.history[:, 1]) plt.plot(ocn.history[:, 0], ocn.history[:, 2]) # plt.xscale("log") -plt.yscale("log") +# plt.yscale("log") plt.show() \ No newline at end of file