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

Spiral profile #1277

wants to merge 2 commits into from

Conversation

esheldon
Copy link
Contributor

@esheldon esheldon commented Mar 6, 2024

Adds a new profile Spiral that creates a logarithmic spiral object, represented by point sources.

The points are created in three dimensions in a disk, with possible inclination along the line of sight and rotations in the plane of the sky. The number of spiral arms is also configurable.

Some examples

two spiral arms

Points in 3d as well as a rendered image

Edge on

More arms

Full scenes

rgb-1021620208
rgb-119650250

The spiral is a lograthmic spiral in 3d that is projected onto 2d.  The
profile is represented by a number of point sources.  The number
of arms, inclination, rotation and hlr/disk height are configurable,
as well as the pitch angle of the arms.

Unit tests are also added in tests/test_spiral.py
@cwwalter
Copy link
Member

cwwalter commented Mar 6, 2024

Cool! Would you still add a bulge component separately?

@esheldon
Copy link
Contributor Author

esheldon commented Mar 6, 2024

Yes, that can be added separately.

In the color images above, the galaxies have up to three parts

  • spiral using the Spiral class
  • knots, using the exact same Spiral profile parameters but with fewer points, more flux per point, and bluer color
  • bulge as a DeVauc profile

@esheldon
Copy link
Contributor Author

esheldon commented Mar 6, 2024

Here you can see the knots clearly
rgb-979104074

@rmjarvis
Copy link
Member

rmjarvis commented Mar 6, 2024

This is quite cool, Erin! The images look great.

@esheldon esheldon requested a review from rmjarvis March 6, 2024 22:26
@esheldon
Copy link
Contributor Author

esheldon commented Mar 6, 2024

Some things for discussion

  • Mike mentioned that nowadays we should be able to replace the _sbp property with something better.
  • The use of points in 3d to represent the profile makes the code simple, but a lot of points are needed to make it look right when the spiral is well resolved. This can make it slow to render, even with photon shooting. Is there some way to hack up an analytic form?

@rmandelb
Copy link
Member

rmandelb commented Mar 7, 2024

Thanks Erin - this is super cool. I was going to ask about computational expense and then you commented on it yourself... it feels to me like there must be some trick one could use, relying on the fact that each point is fundamentally the same exact light profile after PSF convolution, modulo sub-pixel shifts. (I assume we are still making the approximation that there is a single well-defined PSF for the entire Spiral object... i.e. no PSF spatial variation across the Spiral object and no color gradients within the Spiral)

Copy link
Member

@rmjarvis rmjarvis left a comment

Choose a reason for hiding this comment

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

Thanks Erin. I have some requests in detail, but the overall code looks very good. This will be a nice addition!


self._orig_rng = rng.duplicate()
self._rng = rng
self._np_rng = np.random.default_rng(self._rng.raw())
Copy link
Member

Choose a reason for hiding this comment

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

Better would be self._np.rng = self._rng.np.
cf. http://galsim-developers.github.io/GalSim/_build/html/deviate.html#galsim.BaseDeviate.as_numpy_generator
Or just use self._rng.np whenever you want np.random functions.

from .position import PositionD
from coord import degrees, radians

GOLDEN_SPIRAL_PITCH = 17.03239
Copy link
Member

Choose a reason for hiding this comment

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

This deserves a comment giving the source.

# = arctan(ln(phi) / (pi/2))
# cf. https://en.wikipedia.org/wiki/Golden_spiral

Copy link
Member

Choose a reason for hiding this comment

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

I'd also prefer this be a class value in SpiralParticles, rather than a module level constant.


GOLDEN_SPIRAL_PITCH = 17.03239

PIXEL_SCALE = 0.2
Copy link
Member

Choose a reason for hiding this comment

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

I think you don't use this. Which is good, since this shouldn't be hard-coded.


@doc_inherit
def magnify(self, scale):
return self.expand(scale)
Copy link
Member

Choose a reason for hiding this comment

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

This is wrong. It should be expand(mu**0.5). (Where mu is what we normally call the arg here.) Probably should just omit this override though.

def withFlux(self, flux):
kw = self._get_constructor_kw()
kw['flux'] = flux
return Spiral(**kw)
Copy link
Member

Choose a reason for hiding this comment

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

Not sure if this one is necessary or helpful. Rescaling the flux value isn't slow in k-space. (Probably isn't helpful for RandomKnots either.)


def _sample_exponential(rng, r0, n):
return rng.exponential(
scale=r0 * 2.375,
Copy link
Member

Choose a reason for hiding this comment

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

What is 2.375?

np.testing.assert_almost_equal(sp.centroid.y, pts[:, 1].mean())

sp.calculateHLR()
sp.calculate_e1e2()
Copy link
Member

Choose a reason for hiding this comment

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

Should probably assert something about the return values of these two, rather than just that they run without error.

conv2 = galsim.Convolve(sp2, psf)
im1 = conv1.drawImage()
im2 = conv2.drawImage()
assert im1 == im2
Copy link
Member

Choose a reason for hiding this comment

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

I guess this is asserting that xvalue_accuracy and kvalue_accuracy are not used by the Spiral class anywhere. Seems more useful to test that something that does matter causes the image to be different.

# Check that image is not sensitive to use of rng by other objects.
rng3 = galsim.BaseDeviate(seed)
sp3 = galsim.Spiral(npoints, half_light_radius=hlr, rng=rng3)
rng3.discard(523)
Copy link
Member

Choose a reason for hiding this comment

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

I'm not sure this is an interesting test. I think the rng is only used by the constructor. Once the points are made, even if the rng wasn't duplicated, I think this would still work. Perhaps a better test that you are sequestering the rng would be that copying it still produces the same image.

rng3.discard(523)
sp4 = sp3.withFlux(sp3.flux)
sp5 = sp3.copy()
sp6 = sp3.shift(0,0)
...

These should all be equivalent to sp3, and I think exercise the step of duplicating the rng.


check_pickle(sp, irreprable=True)
check_pickle(conv1, irreprable=True)
check_pickle(conv1, lambda x: x.drawImage(scale=1), irreprable=True)
Copy link
Member

Choose a reason for hiding this comment

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

I'd like these to all work with irrerable=False.

@rmjarvis rmjarvis added this to the v2.6 milestone Jun 6, 2024
@rmandelb
Copy link
Member

Hi @esheldon - thanks have been quiet on this PR for a while, so I wanted to check in and let you know that we've been discussing tagging GalSim v2.6 in about a month from now. Do you want this feature in the next tagged version, and if so, do you anticipate that timeline working for you to address the review?

@esheldon
Copy link
Contributor Author

I started this as a fun project with my son while on vacation. I haven't had any time since to work on it.

If I don't get to Mike's comments in time, please don't wait for me.

@rmandelb
Copy link
Member

Ah, the vacation project; I can relate to that phenomenon. Thanks for letting us know... we won't wait on this, but if you do find time for it, we can definitely include it.

@rmandelb rmandelb removed this from the v2.6 milestone Jul 17, 2024
@rmjarvis rmjarvis added the feature request Request for a new feature in GalSim label Jul 24, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
feature request Request for a new feature in GalSim
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants