treverhines/RBF

Add sequence start argument to min_energy_nodes

iraikov opened this issue · 2 comments

Hello,

Sometimes it is useful to provide a different starting point to the Halton sequence used in min_energy_nodes. The attached patch adds option start to min_energy_nodes and rejection_sampling, when then is passed to the underlying HaltonSequence object. Thanks for the consideration.

diff --git a/rbf/pde/nodes.py b/rbf/pde/nodes.py
index 72816e3..21b212e 100644
--- a/rbf/pde/nodes.py
+++ b/rbf/pde/nodes.py
@@ -428,7 +428,7 @@ def prepare_nodes(nodes, domain,
   return nodes, groups, normals
 
   
-def min_energy_nodes(n, domain, rho=None, build_rtree=False, 
+def min_energy_nodes(n, domain, rho=None, build_rtree=False, start=0,
                      **kwargs):
   '''
   Generates nodes within a two or three dimensional. This first
@@ -536,7 +536,7 @@ def min_energy_nodes(n, domain, rho=None, build_rtree=False,
     def rho(x): 
         return np.ones(x.shape[0])
 
-  nodes = rejection_sampling(n, rho, domain)
+  nodes = rejection_sampling(n, rho, domain, start=start)
   out = prepare_nodes(nodes, domain, rho=rho, **kwargs)
   return out                      
 
diff --git a/rbf/pde/sampling.pyx b/rbf/pde/sampling.pyx
index 82b1a07..f12900e 100644
--- a/rbf/pde/sampling.pyx
+++ b/rbf/pde/sampling.pyx
@@ -32,7 +32,7 @@ cdef double distance(list a, list b):
     return sqrt(out)        
             
 
-def rejection_sampling(size, rho, domain, max_sample_size=1000000):
+def rejection_sampling(size, rho, domain, start=0, max_sample_size=1000000):
     '''
     Returns points within the boundaries defined by `vert` and `smp`
     and with density `rho`. The nodes are generated by rejection
@@ -71,7 +71,7 @@ def rejection_sampling(size, rho, domain, max_sample_size=1000000):
     # and 1
     rng = HaltonSequence(1, prime_index=0)
     # form a Halton sequence to generate the test points
-    pnt_rng = HaltonSequence(domain.dim, prime_index=1)
+    pnt_rng = HaltonSequence(domain.dim, prime_index=1, start=start)
     # initiate array of points
     points = np.zeros((0, domain.dim), dtype=float)
     # node counter

Sounds good. I made this change to master

Thank you!