ellipticity_ryden04 has a shape mismatch bug

[Saransh-cpp]

Describe the Bug

The implementation right now keeps sampling values until it gets all positive values -

glass/glass/shapes.py

Lines 142 to 151 in 2eb3e8b

	if size is None:
	size = np.broadcast(gamma, sigma_gamma, mu, sigma).shape
	# draw gamma and epsilon from truncated normal -- eq.s (10)-(11)
	# first sample unbounded normal, then rejection sample truncation
	eps = rng.normal(mu, sigma, size=size)
	bad = eps > 0
	while np.any(bad):
	eps[bad] = rng.normal(mu, sigma, size=eps[bad].shape)
	bad = eps > 0

Let's say if mu = [-1.85, -2.85], gamma = 0.89, gamma = 0.222, sigma_gamma = 0.056 -

• size = (2,)
• eps.shape = (2,)
• rng.normal(mu, sigma) will have the shape (2,) because numpy broadcasts the parameters internally

If one value in eps is bad (<0) -

• bad = [True, False]
• eps[bad].shape = (1,)

then this will not work -

eps[bad] = rng.normal(mu, sigma, size=eps[bad].shape)

as a (2,) array cannot be resized to (1,).

---

This bug was hidden by the float value is not indexable error when the broadcasting rule was not present.

I am not sure what should be done here. Should we consider using something like truncated normal here?

To Reproduce

from glass import ellipticity_ryden04


while True:
    ellipticity_ryden04([-1.85, -2.85], 0.89, 0.222, 0.056)  # rerun till it errors out

Expected Behaviour

No shape mismatch error.

Actual Behaviour

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
      1 while True:
----> 2     ellipticity_ryden04([-1.85, -2.85], 0.89, 0.222, 0.056)

File ~/Code/UCL/glass/glass/shapes.py:146, in ellipticity_ryden04(mu, sigma, gamma, sigma_gamma, size, rng)
    144 bad = eps > 0
    145 while np.any(bad):
--> 146     eps[bad] = rng.normal(mu, sigma, size=eps[bad].shape)
    147     bad = eps > 0
    148 gam = rng.normal(gamma, sigma_gamma, size=size)

File numpy/random/_generator.pyx:1288, in numpy.random._generator.Generator.normal()

File _common.pyx:621, in numpy.random._common.cont()

File _common.pyx:539, in numpy.random._common.cont_broadcast_2()

File _common.pyx:250, in numpy.random._common.validate_output_shape()

ValueError: Output size (1,) is not compatible with broadcast dimensions of inputs (2,).

Version In Use

main

Additional Context

- Python version: 3.11
- Operating system: macOS

[ntessore]

bad = np.where(eps > 0) ?

I see the problem. We probably need to broadcast_arrays() everything to the same shape, then select on all inputs. I don't see how there is another way in the general case.