Dev

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
+- fixed critical bug preventing export to xarray datasets
+
+# 1.3.20251011
+- fixed a bug where optimization did not work for single-outlet networks

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;

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){
+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
-    const double p = exp(-delta_energy * inv_temperature);
-    return rng_uniformdouble(rng) < p;
+    const 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;
         }

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

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"

sandbox.py

@@ -1,10 +1,25 @@
-from timeit import timeit
 
+import matplotlib.pyplot as plt
 import PyOCN as po
+import numpy as np
 
 ocn = po.OCN.from_net_type(
-    net_type="V",
-    dims=(74, 74),
-    random_state=8472,
+    net_type="E",
+    dims=(16, 16),
+    wrap=True,
+    random_state=84712,
 )
-print(timeit(ocn.fit, number=10)/10)
+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))
+
+plt.plot(ocn.history[:, 0], ocn.history[:, 1])
+plt.plot(ocn.history[:, 0], ocn.history[:, 2])
+# plt.xscale("log")
+# plt.yscale("log")
+plt.show()

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()