possible bugs in analytical mass function code and halo profiler code
I think I've identified a bug in the analytical mass function code, and a possible bug or improvement in the halo profiling code: 1) In yt/analysis_modules/halo_mass_function/halo_mass_function.py line 327, shouldn't dn_M_z be multiplied by h^3 instead of h^4? 2) In yt/analysis_modules/halo_profiler/halo_filters.py, I've found that the way the virial radius is determined often doesn't work very well for me. a) I commonly find that more than the first few bins in the profile have overDensity <= virial_overdensity, even though there's nothing wrong with the halo profile. It's just a transient, perhaps caused by too small bin sizes. Just a little bit farther from the center the profile becomes reasonable again, with overDensity > virial_overdensity and decreasing outwards. In these cases, this condition (line 100) if (overDensity[1] <= virial_overdensity): index = 0 results in this halo getting an incorrect virial radius and getting tossed out as "not virialized". b) Related to the previous point, I've found that determining the virial radius (the bin where overDensity drops below virial_overdensity) from the inside-out, as in (line 105) for q in (na.arange(len(overDensity)-2))+2: if (overDensity[q] < virial_overdensity): index = q - 1 break doesn't work as well as an outside-in approach. I've gotten much better results, i.e. fewer rejected halos that were perfectly fine and better virial radii, by doing this: if (overDensity[-1] >= virial_overdensity): index = -2 else: for q in (na.arange(len(overDensity),0,-1)-1): if (overDensity[q] < virial_overdensity) and (overDensity[q-1] > virial_overdensity): index = q - 1 break 3) Lastly, one question about the radial profiles. I've noticed that the 'CellVolume' field is rarely equal to 4.0/3.0 * pi * ('RadiusMpc' * 3.09e24)**3. Typically 'CellVolume' is about 35 - 40% larger. Maybe this is because the grid cells partly "overshoot" the edge of the sphere? But I'm surprised that this is a 35% effect... This potentially affects the determination of the virial radius. Thanks, Mike -- ********************************************************************* * * * Dr. Michael Kuhlen Theoretical Astrophysics Center * * email: mqk@astro.berkeley.edu UC Berkeley * * cell phone: (831) 588-1468 601 Campbell Hall * * skype username: mikekuhlen Berkeley, CA 94720 * * * *********************************************************************
Hi Mike, I'll leave the rest for Britton, Stephen, etc, but I can address your final point.
3) Lastly, one question about the radial profiles. I've noticed that the 'CellVolume' field is rarely equal to 4.0/3.0 * pi * ('RadiusMpc' * 3.09e24)**3. Typically 'CellVolume' is about 35 - 40% larger. Maybe this is because the grid cells partly "overshoot" the edge of the sphere? But I'm surprised that this is a 35% effect... This potentially affects the determination of the virial radius.
CellVolume is always calculated by summing up the values of the cells that are included, which are determined by the AMRSphere object. There are a couple potential weaknesses here, but they will be consistent with the cells that have been included. This is done in two routines: yt/data_objects/object_finding_mixin.py line 162: ObjectFindingMixin.find_sphere_grids yt/data_objects/data_containers.py line 2936: AMRSphereBase._get_cut_mask (also see _is_fully_enclosed just above) This uses the RadiusCode field, which is defined on line 746 of yt/data_objects/universal_fields.py. The is_fully_contained (which you can verify is not causing the problem inserting "return False" at the top of that routine and checking again) routing determines whether or not the grid requires additional masking. The masking is done by comparing RadiusCode to the radius of the sphere. RadiusCode comes from the center of the cells. If you catch too many root grid cells, for instance, this may result in an overestimate. It's possible that this could be 30-40%, if it's integrated over the entire sphere, catching cells here and there on the exterior that contribute far too much. (You could, at best, be adding on an entire ~3/4 of a cell, and effectively extending the radius by sqrt(3) dx.) It's worth investigating if this is indeed where the issue is coming from. That aside, the selection method is identical everywhere, so I do believe that taking the sum of CellVolume is the correct method, rather than using 4/3 pi r^3. -Matt
Thanks,
Mike
-- ********************************************************************* * * * Dr. Michael Kuhlen Theoretical Astrophysics Center * * email: mqk@astro.berkeley.edu UC Berkeley * * cell phone: (831) 588-1468 601 Campbell Hall * * skype username: mikekuhlen Berkeley, CA 94720 * * * ********************************************************************* _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
Hi Mike, I've found that the issue of overdensity being lower in the innermost bins can be resolved with a recentering of the sphere before profiling. I believe that the current default is to use the location of maximum dark matter density as the center of the sphere. Stephen has done a lot of work on this and has written some functionality to recenter profiles on the location of maximum baryon density, for example. I'm not sure what the right answer is for how profiles should be centered. As Stephen has found, it gets messier when a halo is undergoing a significant merger. All that aside, I think I agree with calculating the virial properties from the outside in, rather than the inside out. I support that change being made to the code. If you've already got it implemented, feel free to push it. Nice work! Britton On Tue, May 24, 2011 at 3:26 PM, Matthew Turk <matthewturk@gmail.com> wrote:
Hi Mike,
I'll leave the rest for Britton, Stephen, etc, but I can address your final point.
3) Lastly, one question about the radial profiles. I've noticed that the 'CellVolume' field is rarely equal to 4.0/3.0 * pi * ('RadiusMpc' * 3.09e24)**3. Typically 'CellVolume' is about 35 - 40% larger. Maybe this is because the grid cells partly "overshoot" the edge of the sphere? But I'm surprised that this is a 35% effect... This potentially affects the determination of the virial radius.
CellVolume is always calculated by summing up the values of the cells that are included, which are determined by the AMRSphere object. There are a couple potential weaknesses here, but they will be consistent with the cells that have been included. This is done in two routines:
yt/data_objects/object_finding_mixin.py line 162: ObjectFindingMixin.find_sphere_grids yt/data_objects/data_containers.py line 2936: AMRSphereBase._get_cut_mask (also see _is_fully_enclosed just above)
This uses the RadiusCode field, which is defined on line 746 of yt/data_objects/universal_fields.py.
The is_fully_contained (which you can verify is not causing the problem inserting "return False" at the top of that routine and checking again) routing determines whether or not the grid requires additional masking. The masking is done by comparing RadiusCode to the radius of the sphere. RadiusCode comes from the center of the cells. If you catch too many root grid cells, for instance, this may result in an overestimate. It's possible that this could be 30-40%, if it's integrated over the entire sphere, catching cells here and there on the exterior that contribute far too much. (You could, at best, be adding on an entire ~3/4 of a cell, and effectively extending the radius by sqrt(3) dx.) It's worth investigating if this is indeed where the issue is coming from.
That aside, the selection method is identical everywhere, so I do believe that taking the sum of CellVolume is the correct method, rather than using 4/3 pi r^3.
-Matt
Thanks,
Mike
-- ********************************************************************* * * * Dr. Michael Kuhlen Theoretical Astrophysics Center * * email: mqk@astro.berkeley.edu UC Berkeley * * cell phone: (831) 588-1468 601 Campbell Hall * * skype username: mikekuhlen Berkeley, CA 94720 * * * ********************************************************************* _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
Ah, that's a good point about the re-centering. I see that in the yt 'stable' branch there's a 'use_density_center' option, although the docs say "This is generally not needed." But I'm using the 'yt' branch, and there 'use_density_center' has been replaced with a 'recenter' option. It looks like I have to pass a function that performs the halo re-centering. What does this function have to look like? Is there an example anywhere?
All that aside, I think I agree with calculating the virial properties from the outside in, rather than the inside out. I support that change being made to the code. If you've already got it implemented, feel free to push it.
I just tried pushing this change, but apparently I don't have commit permissions for yt. I think Matt set that up for me a while ago with the same password as for Enzo, but that doesn't work for me now. Cheers, Mike On Tue, May 24, 2011 at 1:20 PM, Britton Smith <brittonsmith@gmail.com> wrote:
Hi Mike,
I've found that the issue of overdensity being lower in the innermost bins can be resolved with a recentering of the sphere before profiling. I believe that the current default is to use the location of maximum dark matter density as the center of the sphere. Stephen has done a lot of work on this and has written some functionality to recenter profiles on the location of maximum baryon density, for example. I'm not sure what the right answer is for how profiles should be centered. As Stephen has found, it gets messier when a halo is undergoing a significant merger.
All that aside, I think I agree with calculating the virial properties from the outside in, rather than the inside out. I support that change being made to the code. If you've already got it implemented, feel free to push it.
Nice work!
Britton
On Tue, May 24, 2011 at 3:26 PM, Matthew Turk <matthewturk@gmail.com> wrote:
Hi Mike,
I'll leave the rest for Britton, Stephen, etc, but I can address your final point.
3) Lastly, one question about the radial profiles. I've noticed that the 'CellVolume' field is rarely equal to 4.0/3.0 * pi * ('RadiusMpc' * 3.09e24)**3. Typically 'CellVolume' is about 35 - 40% larger. Maybe this is because the grid cells partly "overshoot" the edge of the sphere? But I'm surprised that this is a 35% effect... This potentially affects the determination of the virial radius.
CellVolume is always calculated by summing up the values of the cells that are included, which are determined by the AMRSphere object. There are a couple potential weaknesses here, but they will be consistent with the cells that have been included. This is done in two routines:
yt/data_objects/object_finding_mixin.py line 162: ObjectFindingMixin.find_sphere_grids yt/data_objects/data_containers.py line 2936: AMRSphereBase._get_cut_mask (also see _is_fully_enclosed just above)
This uses the RadiusCode field, which is defined on line 746 of yt/data_objects/universal_fields.py.
The is_fully_contained (which you can verify is not causing the problem inserting "return False" at the top of that routine and checking again) routing determines whether or not the grid requires additional masking. The masking is done by comparing RadiusCode to the radius of the sphere. RadiusCode comes from the center of the cells. If you catch too many root grid cells, for instance, this may result in an overestimate. It's possible that this could be 30-40%, if it's integrated over the entire sphere, catching cells here and there on the exterior that contribute far too much. (You could, at best, be adding on an entire ~3/4 of a cell, and effectively extending the radius by sqrt(3) dx.) It's worth investigating if this is indeed where the issue is coming from.
That aside, the selection method is identical everywhere, so I do believe that taking the sum of CellVolume is the correct method, rather than using 4/3 pi r^3.
-Matt
Thanks,
Mike
-- ********************************************************************* * * * Dr. Michael Kuhlen Theoretical Astrophysics Center * * email: mqk@astro.berkeley.edu UC Berkeley * * cell phone: (831) 588-1468 601 Campbell Hall * * skype username: mikekuhlen Berkeley, CA 94720 * * * ********************************************************************* _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
-- ********************************************************************* * * * Dr. Michael Kuhlen Theoretical Astrophysics Center * * email: mqk@astro.berkeley.edu UC Berkeley * * cell phone: (831) 588-1468 601 Campbell Hall * * skype username: mikekuhlen Berkeley, CA 94720 * * * *********************************************************************
Hi Mike,
Ah, that's a good point about the re-centering. I see that in the yt 'stable' branch there's a 'use_density_center' option, although the docs say "This is generally not needed." But I'm using the 'yt' branch, and there 'use_density_center' has been replaced with a 'recenter' option. It looks like I have to pass a function that performs the halo re-centering. What does this function have to look like? Is there an example anywhere?
Yes, the new recentering stuff is slated to be fully released with the 2.2 release. I've added docs to the yt-docs mercurial repository, but bitbucket is having issues right now, so I can't point you to them, and I'm assuming you don't have a copy locally already. So I've pasted the relevant stuff here: http://paste.enzotools.org/show/1667/ Let me know if that's not helpful! -- Stephen Skory s@skory.us http://stephenskory.com/ 510.621.3687 (google voice)
The yt-docs are nominally available in the directory $YT_DEST/src/yt-supplemental/yt-docs , which gets updated when you do "yt instinfo -u", but that repository may not be up to date with the yt-doc repo on BB because of how subrepos work in hg. Bitbucket is working for me, by the way: http://hg.enzotools.org/yt-doc Mike, if you email me with your BB username, I'll add you to the push ACL. -Matt On Tue, May 24, 2011 at 5:26 PM, Stephen Skory <s@skory.us> wrote:
Hi Mike,
Ah, that's a good point about the re-centering. I see that in the yt 'stable' branch there's a 'use_density_center' option, although the docs say "This is generally not needed." But I'm using the 'yt' branch, and there 'use_density_center' has been replaced with a 'recenter' option. It looks like I have to pass a function that performs the halo re-centering. What does this function have to look like? Is there an example anywhere?
Yes, the new recentering stuff is slated to be fully released with the 2.2 release. I've added docs to the yt-docs mercurial repository, but bitbucket is having issues right now, so I can't point you to them, and I'm assuming you don't have a copy locally already. So I've pasted the relevant stuff here:
http://paste.enzotools.org/show/1667/
Let me know if that's not helpful!
-- Stephen Skory s@skory.us http://stephenskory.com/ 510.621.3687 (google voice) _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
Hi Mike,
1) In yt/analysis_modules/halo_mass_function/halo_mass_function.py line 327, shouldn't dn_M_z be multiplied by h^3 instead of h^4?
For easy reference here's the function we're talking about: https://bitbucket.org/yt_analysis/yt/src/f8d1ca7f701c/yt/analysis_modules/ha... Thanks for pointing this out, and now that I take a closer look at it I'm slightly confused, too, and I may need a pinch hitter to weigh in on this (AKA passing the buck). The history of that code is I took some C code written by Brian O'Shea and I translated it into Python. The C version required re-compiling for any changes, linked to the GNU scientific library (which is ubiquitous, of course), and involved a multi-step process to get the end product. Converting it to Python eliminated these issues, and including it in yt made it more convenient. This section is basically verbatim from Brian (code and comments), which is clear from the extraneous end-of-expression semicolons I was too lazy to remove. You're probably getting the h^3 power based on the comments a few lines above #327. But as I look at those comments, it looks like to me that a) the power should actually be h^2 (= h^2 * h * h^-1) and that b) line 327 should actually be a divide to remove the powers of h. Brian, if you have a moment, can you dig into your memory vault and take a look at what Mike has pointed out? Or are we just confused by the comments which are misleading? Thanks! -- Stephen Skory s@skory.us http://stephenskory.com/ 510.621.3687 (google voice)
Hi Stephen I think it should be h^3, not h^2, because as the comment above says, "rho0 has units of h^2, dsigmadm has units of h, and massarray has units of h^-1". The h^-1 factors from massarray actually cancel, since you multiply by (self.massarray[i+1] - self.massarray[i]) / self.massarray[i]. And it should be multiplication, not division, because (for example) [rho0] = (Msun/h) / (Mpc/h)**3 = Msun/Mpc**3 * h^2, so you have to multiply rho0 by h^2 to go from (Msun/h)/(Mpc/h)**3 to Msun/Mpc**3. Years ago I wrote my own code to calculate P&S, S&T, Jenkins, etc. mass functions, and I did not use units scaled by h. When I change the yt code to multiplication by h^3 I find agreement between the yt mass functions and my own. Cheers, Mike On Tue, May 24, 2011 at 1:20 PM, Stephen Skory <s@skory.us> wrote:
Hi Mike,
1) In yt/analysis_modules/halo_mass_function/halo_mass_function.py line 327, shouldn't dn_M_z be multiplied by h^3 instead of h^4?
For easy reference here's the function we're talking about: https://bitbucket.org/yt_analysis/yt/src/f8d1ca7f701c/yt/analysis_modules/ha...
Thanks for pointing this out, and now that I take a closer look at it I'm slightly confused, too, and I may need a pinch hitter to weigh in on this (AKA passing the buck). The history of that code is I took some C code written by Brian O'Shea and I translated it into Python. The C version required re-compiling for any changes, linked to the GNU scientific library (which is ubiquitous, of course), and involved a multi-step process to get the end product. Converting it to Python eliminated these issues, and including it in yt made it more convenient. This section is basically verbatim from Brian (code and comments), which is clear from the extraneous end-of-expression semicolons I was too lazy to remove.
You're probably getting the h^3 power based on the comments a few lines above #327. But as I look at those comments, it looks like to me that a) the power should actually be h^2 (= h^2 * h * h^-1) and that b) line 327 should actually be a divide to remove the powers of h.
Brian, if you have a moment, can you dig into your memory vault and take a look at what Mike has pointed out? Or are we just confused by the comments which are misleading? Thanks!
-- Stephen Skory s@skory.us http://stephenskory.com/ 510.621.3687 (google voice) _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
-- ********************************************************************* * * * Dr. Michael Kuhlen Theoretical Astrophysics Center * * email: mqk@astro.berkeley.edu UC Berkeley * * cell phone: (831) 588-1468 601 Campbell Hall * * skype username: mikekuhlen Berkeley, CA 94720 * * * *********************************************************************
Hi all, To follow up: I concur that it should be h^3, not h^4. I think this is a typo in the original code. Thanks for catching this! --Brian On Tue, May 24, 2011 at 4:52 PM, Michael Kuhlen <mqk@astro.berkeley.edu>wrote:
Hi Stephen
I think it should be h^3, not h^2, because as the comment above says, "rho0 has units of h^2, dsigmadm has units of h, and massarray has units of h^-1". The h^-1 factors from massarray actually cancel, since you multiply by (self.massarray[i+1] - self.massarray[i]) / self.massarray[i].
And it should be multiplication, not division, because (for example) [rho0] = (Msun/h) / (Mpc/h)**3 = Msun/Mpc**3 * h^2, so you have to multiply rho0 by h^2 to go from (Msun/h)/(Mpc/h)**3 to Msun/Mpc**3.
Years ago I wrote my own code to calculate P&S, S&T, Jenkins, etc. mass functions, and I did not use units scaled by h. When I change the yt code to multiplication by h^3 I find agreement between the yt mass functions and my own.
Cheers,
Mike
On Tue, May 24, 2011 at 1:20 PM, Stephen Skory <s@skory.us> wrote:
Hi Mike,
1) In yt/analysis_modules/halo_mass_function/halo_mass_function.py line 327, shouldn't dn_M_z be multiplied by h^3 instead of h^4?
For easy reference here's the function we're talking about:
https://bitbucket.org/yt_analysis/yt/src/f8d1ca7f701c/yt/analysis_modules/ha...
Thanks for pointing this out, and now that I take a closer look at it I'm slightly confused, too, and I may need a pinch hitter to weigh in on this (AKA passing the buck). The history of that code is I took some C code written by Brian O'Shea and I translated it into Python. The C version required re-compiling for any changes, linked to the GNU scientific library (which is ubiquitous, of course), and involved a multi-step process to get the end product. Converting it to Python eliminated these issues, and including it in yt made it more convenient. This section is basically verbatim from Brian (code and comments), which is clear from the extraneous end-of-expression semicolons I was too lazy to remove.
You're probably getting the h^3 power based on the comments a few lines above #327. But as I look at those comments, it looks like to me that a) the power should actually be h^2 (= h^2 * h * h^-1) and that b) line 327 should actually be a divide to remove the powers of h.
Brian, if you have a moment, can you dig into your memory vault and take a look at what Mike has pointed out? Or are we just confused by the comments which are misleading? Thanks!
-- Stephen Skory s@skory.us http://stephenskory.com/ 510.621.3687 (google voice) _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
-- ********************************************************************* * * * Dr. Michael Kuhlen Theoretical Astrophysics Center * * email: mqk@astro.berkeley.edu UC Berkeley * * cell phone: (831) 588-1468 601 Campbell Hall * * skype username: mikekuhlen Berkeley, CA 94720 * * * ********************************************************************* _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
Brian and Mike,
To follow up: I concur that it should be h^3, not h^4. I think this is a typo in the original code. Thanks for catching this!
I guess that settles that. Thanks Brian and Mike! I've just pushed a change to the halo mass function. -- Stephen Skory s@skory.us http://stephenskory.com/ 510.621.3687 (google voice)
participants (5)
-
Brian O'Shea -
Britton Smith -
Matthew Turk -
Michael Kuhlen -
Stephen Skory