Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Spiral profile #1277

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
get coverage to 100%, do some cleanup
esheldon committed Mar 6, 2024
commit 11f1d577b68eafa46f95008933130d7339d06961
85 changes: 29 additions & 56 deletions galsim/spiral.py
Original file line number Diff line number Diff line change
@@ -259,14 +259,6 @@ def input_half_light_radius(self):
"""
return self._half_light_radius

def _get_particles(self):
if not hasattr(self, '_part'):
raise RuntimeError(
'Cannot get particles, Spiral was constructed from '
'points'
)
return self._part

@property
def narms(self):
"""
@@ -655,16 +647,24 @@ def __init__(
rel_height=0.1,
inclination=None,
rotation=None,
shift=None,
):

self.rng = rng
self.narms = int(narms)
if self.narms <= 0:
raise ValueError(f'number of arms should be > 0, got {self.narms}')

if self.narms <= 1:
raise GalSimValueError(
f'number of arms should be > 0, got {self.narms}',
self.narms,
)

self.hlr = float(hlr)

if self.hlr <= 0:
raise ValueError(f'hlr must be > 0, got {self.hlr}')
raise GalSimValueError(
f'hlr must be > 0, got {self.hlr}',
self.hlr,
)

self.r0 = self.hlr / 1.67835
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could use Exponential._hlr_factor rather than a bare number.


@@ -684,12 +684,10 @@ def __init__(
self.pitch = pitch
self.inclination = inclination
self.rotation = rotation
self.shift = shift

_check_angle(self.pitch, 'pitch')
_check_angle(self.inclination, 'inclination')
_check_angle(self.rotation, 'rotation')
_check_sequence(self.shift, 'shift', 2)

def sample(self, n):
"""
@@ -702,33 +700,30 @@ def sample(self, n):

theta = np.log(r / self.r0) / self.pitch.tan()

if self.angle_fuzz is not None:
theta += theta * self.rng.normal(
scale=self.angle_fuzz,
size=theta.size,
)
theta += theta * self.rng.normal(
scale=self.angle_fuzz,
size=theta.size,
)

# split points into arms
if self.narms > 1:
angle = 2 * np.pi / self.narms
frac = int(r.size / self.narms)
start = frac
for i in range(self.narms-1):
start = (i+1) * frac
end = (i+2) * frac
theta[start:end] += (i+1) * angle
angle = 2 * np.pi / self.narms
frac = int(r.size / self.narms)
start = frac
for i in range(self.narms-1):
start = (i+1) * frac
end = (i+2) * frac
theta[start:end] += (i+1) * angle

xyz = np.zeros((n, 3))
xyz[:, 0] = r * np.cos(theta)
xyz[:, 1] = r * np.sin(theta)

if self.xy_fuzz is not None:
xyz[:, 0] += self.rng.normal(
scale=self.xy_fuzz * self.hlr, size=r.size,
)
xyz[:, 1] += self.rng.normal(
scale=self.xy_fuzz * self.hlr, size=r.size,
)
xyz[:, 0] += self.rng.normal(
scale=self.xy_fuzz * self.hlr, size=r.size,
)
xyz[:, 1] += self.rng.normal(
scale=self.xy_fuzz * self.hlr, size=r.size,
)

xyz[:, 2] = self._get_z(r.size)

@@ -744,10 +739,6 @@ def sample(self, n):
xyz[:, 0], xyz[:, 1], self.rotation,
)

if self.shift is not None:
xyz[:, 0] += self.shift[0]
xyz[:, 1] += self.shift[1]

return xyz
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it useful to return 3 coords? I think once you've finished with the rotations, you just use x, y. Maybe this should just return xyz[:,:2].


def _get_z(self, n):
@@ -792,21 +783,3 @@ def _check_angle(angle, name):
raise GalSimValueError(
f'{name} should be a galsim.Angle, got {tangle}', angle
)


def _check_sequence(data, name, nel):
if data is None:
return

try:
ndata = len(data)
if ndata != nel:
raise GalSimValueError(
f'expected {nel} element sequence for {name}, got {ndata}',
data
)
except TypeError:
raise GalSimValueError(
f'expected sequence for {name}, got {type(data)}',
data
)
89 changes: 83 additions & 6 deletions tests/test_spiral.py
Original file line number Diff line number Diff line change
@@ -136,6 +136,18 @@ def test_spiral_valid_input(use_config):
kw['rng'] = rng
sp = galsim.Spiral(**kw)

kw2 = kw.copy()
pitch = 20 * galsim.degrees
kw2['pitch'] = pitch
tmp_sp = galsim.Spiral(**kw2)
assert tmp_sp.pitch == pitch

del kw2['inclination']
galsim.Spiral(**kw2)

del kw2['rotation']
galsim.Spiral(**kw2)

pts = sp.points
nobj = pts.shape[0]
assert nobj == kw['npoints'], (
@@ -228,6 +240,53 @@ def test_spiral_invalid_inputs():
with pytest.raises(GalSimRangeError):
galsim.Spiral(**this_kw)

with pytest.raises(RuntimeError):
# need npoints= or points=
galsim.Spiral(half_light_radius=2)

with pytest.raises(GalSimValueError):
galsim.Spiral(points=np.zeros((10, 4)))

with pytest.raises(GalSimValueError):
galsim.Spiral(points='blah')

with pytest.raises(TypeError):
galsim.Spiral(npoints=100, half_light_radius=1, rng='blah')

with pytest.raises(TypeError):
galsim.Spiral(npoints=100, half_light_radius=1).shear(
galsim.Shear(g1=0.5),
g2='blah',
)

with pytest.raises(TypeError):
galsim.Spiral(npoints=100, half_light_radius=1).shear(0.1)

with pytest.raises(TypeError):
galsim.Spiral(npoints=100, half_light_radius=1).shear(0.1, 0.2, 0.3)

with pytest.raises(TypeError):
galsim.Spiral(npoints=100, half_light_radius=1).shear()

with pytest.raises(TypeError):
galsim.Spiral(npoints=100, half_light_radius=1).rotate(0)

nprng = np.random.default_rng(88)
with pytest.raises(GalSimValueError):
galsim.spiral.SpiralParticles(hlr=1, narms=-3, rng=nprng)

with pytest.raises(GalSimValueError):
galsim.spiral.SpiralParticles(hlr=-1, rng=nprng)

with pytest.raises(GalSimValueError):
galsim.spiral.SpiralParticles(hlr=1, rng=nprng, pitch=0)

with pytest.raises(GalSimValueError):
galsim.spiral.SpiralParticles(hlr=1, inclination=0, rng=nprng)

with pytest.raises(GalSimValueError):
galsim.spiral.SpiralParticles(hlr=1, rotation=0, rng=nprng)


def test_transforms():

@@ -237,6 +296,7 @@ def test_transforms():
kw['rng'] = rng

sp = galsim.Spiral(**kw)
spp = galsim.Spiral(points=sp.points)
x, y = sp.points[:, 0], sp.points[:, 1]

# flux
@@ -254,11 +314,17 @@ def test_transforms():
tsp = sp.dilate(dilate)
assert tsp.input_half_light_radius == kw['half_light_radius'] * dilate

tsp = spp.dilate(dilate)
assert tsp.input_half_light_radius is None

# expand and magnify
tsp = sp.expand(dilate)
assert tsp.input_half_light_radius == kw['half_light_radius'] * dilate
assert tsp.flux == kw['flux'] / dilate**2

tsp = spp.expand(dilate)
assert tsp.input_half_light_radius is None

tsp = sp.magnify(dilate)
assert tsp.input_half_light_radius == kw['half_light_radius'] * dilate
assert tsp.flux == kw['flux'] / dilate**2
@@ -283,6 +349,9 @@ def test_transforms():
(tsp.points[:, 1] == yp)
)

tsp = spp.transform(dudx, dudy, dvdx, dvdy)
assert tsp.input_half_light_radius is None

# rotate
rot = 20 * galsim.degrees
tsp = sp.rotate(rot)
@@ -302,6 +371,7 @@ def test_transforms():
# shear
sh = galsim.Shear(g1=0.1, g2=-0.05)
tsp = sp.shear(sh)
tspkw = sp.shear(g1=sh.g1, g2=sh.g2)

# we already tested transform above, use it again
mat = sh.getMatrix()
@@ -311,11 +381,13 @@ def test_transforms():
dvdx=mat[1, 0],
dvdy=mat[1, 1],
)
assert np.all(
(tsp.points[:, 0] == tsp2.points[:, 0])
&
(tsp.points[:, 1] == tsp2.points[:, 1])
)

for _tsp in (tsp, tspkw):
assert np.all(
(_tsp.points[:, 0] == tsp2.points[:, 0])
&
(_tsp.points[:, 1] == tsp2.points[:, 1])
)


def test_spiral_hlr():
@@ -432,4 +504,9 @@ def test_spiral_sed():


if __name__ == '__main__':
test_spiral_valid_input(True)
testfns = [
v for k, v in vars().items()
if k[:5] == 'test_' and callable(v)
]
for testfn in testfns:
testfn()