Merge branch 'test-hybridIA' into 'development'

putting understanding of hybridIA into code

See merge request damask/DAMASK!667

python/damask/util.py

@@ -416,7 +416,7 @@ def project_equal_area(vector: _np.ndarray,
                     -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
 
 
-def hybrid_IA(dist: _np.ndarray,
+def hybrid_IA(dist: _FloatSequence,
               N: int,
               rng_seed: _NumpyRngSeed = None) -> _np.ndarray:
     """
@@ -425,19 +425,25 @@ def hybrid_IA(dist: _np.ndarray,
     Parameters
     ----------
     dist : numpy.ndarray
-        Distribution to be approximated
+        Distribution to be approximated.
     N : int
         Number of samples to draw.
     rng_seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
         A seed to initialize the BitGenerator. Defaults to None.
         If None, then fresh, unpredictable entropy will be pulled from the OS.
 
+    Returns
+    -------
+    hist : numpy.ndarray, shape (N)
+        Integer approximation of the distribution.
+
     """
-    N_opt_samples,N_inv_samples = (max(_np.count_nonzero(dist),N),0)                                 # random subsampling if too little samples requested
+    N_opt_samples = max(_np.count_nonzero(dist),N)                                                  # random subsampling if too little samples requested
+    N_inv_samples = 0
 
     scale_,scale,inc_factor = (0.0,float(N_opt_samples),1.0)
     while (not _np.isclose(scale, scale_)) and (N_inv_samples != N_opt_samples):
-        repeats = _np.rint(scale*dist).astype(_np.int64)
+        repeats = _np.rint(scale*_np.array(dist)).astype(_np.int64)
         N_inv_samples = _np.sum(repeats)
         scale_,scale,inc_factor = (scale,scale+inc_factor*0.5*(scale - scale_), inc_factor*2.0) \
                                    if N_inv_samples < N_opt_samples else \

python/tests/test_util.py

@@ -43,7 +43,7 @@ class TestUtil:
 
 
     @pytest.mark.parametrize('rv',[stats.rayleigh(),stats.weibull_min(1.2),stats.halfnorm(),stats.pareto(2.62)])
-    def test_hybridIA(self,rv):
+    def test_hybridIA_distribution(self,rv):
         bins = np.linspace(0,10,100000)
         centers = (bins[1:]+bins[:-1])/2
         N_samples = bins.shape[0]-1000
@@ -52,6 +52,21 @@ class TestUtil:
         dist_sampled = np.histogram(centers[selected],bins)[0]/N_samples*np.sum(dist)
         assert np.sqrt(((dist - dist_sampled) ** 2).mean()) < .025 and selected.shape[0]==N_samples
 
+    def test_hybridIA_constant(self):
+       N_bins = np.random.randint(20,400)
+       m = np.random.randint(1,20)
+       N_samples = m * N_bins
+       dist = np.ones(N_bins)*np.random.rand()
+       assert np.all(np.sort(util.hybrid_IA(dist,N_samples))==np.arange(N_samples).astype(int)//m)
+
+    def test_hybridIA_linear(self):
+       N_points = np.random.randint(10,200)
+       m = np.random.randint(1,20)
+       dist = np.arange(N_points)
+       N_samples = m * np.sum(dist)
+       assert np.all(np.bincount(util.hybrid_IA(dist*np.random.rand(),N_samples)) == dist*m)
+
+
     @pytest.mark.parametrize('point,direction,normalize,keepdims,answer',
                              [
                               ([1,0,0],'z',False,True, [1,0,0]),