[Yt-dev] ellipsoid part 2

[Matthew Turk]

Hi Geoffrey,

Great to hear from you again on this project!  I have several comments
about the code interspersed; they are all relatively minor, but there
are several.

On Mon, Sep 5, 2011 at 10:23 PM, Geoffrey So <gsiisg at gmail.com> wrote:
> Hello Stephen,
> This is part 2 of the ellipsoid parameter finding.  The method is outlined
> mainly in the Dubinski and Carlberg 1991 paper.  But since we are doing this
> in YT and we already know the particle list to begin with, I put in the
> option that the user can specify if they want all the particles in the
> region to be taken into account or only the particles found by the halo
> finder algorithm. This is an option that I think should be amended to the
> first ellipsoid parameter finder.

Is this the parameter halo_only?  I would like to note that
creation_time < 0 is not reliable, and will make your code
Enzo-specific.  I recommend for now using the boolean selector IsStar:

dm_index = dd["IsStar"]

Additionally, I would recommend not constructing a new array (which
will be slower than allowing the array selections to be left as is)
and to instead set ax (which should be renamed, as 'ax' is
idiomatically used in yt to reference an axis) to a list.

I would go a bit further and suggest you not use the parameter
halo_only as it is implicit in the halo finding.  If the halo finding
subselected on dark matter, then you would want that in the ellipsoid.
 If it didn't subselect, then you probably also don't want to
subselect in the ellipsoid finding.

> The script form is linked here:
> http://paste.enzotools.org/show/1782/
> Once I have written out the methodology for myself, I'll work on the YT
> documentation, will probably need some help in this area when it comes to
> it.

I have a few comments:

1) Use booleans for parameters such as use_halos.  True and False, not 0 and 1.
2) Iterate over the halos:

for halo in halos:
3) Index variables should be clearly index variables; "halo" is an
object, but you use it as an index.
4) cm => com
5) You set ax to be an array twice, which will double the memory
requirements.  I don't think ax needs to be an array at all.
6) Line 109, just modify ax in place.  Fancy-indexing in numpy will (I
believe) implicitly make a copy when you modify:

In [1]: import numpy as na
In [2]: arr = na.random.random(100)
In [3]: ai = arr > 0.8
In [4]: bb = arr[ai]
In [5]: bb -= 1.0
In [6]: arr.max()
Out[6]: 0.99642788644174751
In [7]: bb.max()
Out[7]: -0.0035721135582524877

7) Line 113, use range or xrange, and use singular "axis" not "axes"
8) Line 114, you double the memory by making a new array again.
9) Line 117, abs is a numpy function.
10) Line 121, this makes four new arrays, each of which is the length
of the number of particles.  Construct a temporary array and do
in-place math.  You don't need "where" here; use booleans.
11) Your calculation of the inertial tensor is nice, but I would again
recommend that you avoid creating four temporary arrays when
calculating the magnitude of 'ax'.
12) I would encourage you to investigate the eigenvector code if you
suspect it gives LH system coordinates.
13) "for iter in range(0, maxiterations)") should just be for iter in
xrange(maxiterations).  There is no need to specify 0.
14) Please use print without parenthesis.  Python 3 has transformed
print into a function, but the formatting rules are also different.
15) Avoid using where!
16) use .size instead of len() for arrays
17) Feel free to throw an error.  Define a subclass of Exception and
then raise an instance of it.

One fun thing would also be to plot projections of "Ones" using these
ellipsoids as sources.  That would help us see them better, too!

> Doron Lemez has used this method in his paper and helped me tested the
> accuracy of this script.  I was wondering do I need to ask him if he wants
> his name in the credits, or just mention his name as a special thanks?

Which credits do you mean?  The CREDITS file is reserved for code
contributions.  Feel free in your documentation to mention that he
assisted you, if you feel it is warranted.  Please be sure that your
documentation includes the ".. sectionauthor" directive to indicate
that you wrote it.

Please issue a pull request once you have included this routine into
the halo finding module; we can directly review it then.

Additionally, for the documentation modifications, I would encourage
you to fork yt-doc and make your changes in source/analysis_modules/
under the halo finding section.  Include a few (small file size)
images, too.

Thanks, GS!

-Matt

>  One
> thing he did in the paper that was especially useful is defined the starting
> sphere to have 2x the virial radius.  Without that specification, I run into
> non-convergence sometimes, where after many iterations I end up with no
> particles.  And I want to mention again that there's no guarantee that any
> or all the haloes will converge with this method.
> I tried to adhere to the PEP8 formatting, but old habits die hard so I'm
> sure I made mistakes in the formatting, let me know if you see them.
> From
> G.S.
> _______________________________________________
> Yt-dev mailing list
> Yt-dev at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
>
>