<div>Matt, Stephen,</div><div><br></div>Thanks Matt for all the suggestions. I've made most of the changes, but some I wasn't quite sure what you meant:<div><br></div><div>1) Use booleans for parameters such as use_halos. True and False, not 0 and 1.<br>
- Done</div><div>2) Iterate over the halos:<br>- Done (I assume this is faster because it doesn't try to find which halo every time I call it?)</div><div>3) Index variables should be clearly index variables; "halo" is an<br>
object, but you use it as an index.</div><div>- Done, switched to calling halo_num to avoid ambiguity</div><div>4) cm => com</div><div>- Done<br>5) You set ax to be an array twice, which will double the memory<br>requirements. I don't think ax needs to be an array at all.<br>
- Switched to using list [], not sure why an array of arrays is better than list of arrays, but I take your word for it</div><div>6) Line 109, just modify ax in place. Fancy-indexing in numpy will (I<br>believe) implicitly make a copy when you modify:<br>
<br>In [1]: import numpy as na<br>In [2]: arr = na.random.random(100)<br>In [3]: ai = arr > 0.8<br>In [4]: bb = arr[ai]<br>In [5]: bb -= 1.0<br>In [6]: arr.max()<br>Out[6]: 0.99642788644174751<br>In [7]: bb.max()<br>Out[7]: -0.0035721135582524877<br>
<br></div><div>- Not quite sure was meant by the example</div><div><br>7) Line 113, use range or xrange, and use singular "axis" not "axes"</div><div>- Done<br>8) Line 114, you double the memory by making a new array again.<br>
- Done, switched to list</div><div>9) Line 117, abs is a numpy function.</div><div>- Done, switched to na.abs()<br>10) Line 121, this makes four new arrays, each of which is the length<br>of the number of particles. Construct a temporary array and do<br>
in-place math. You don't need "where" here; use booleans.</div><div>- Old code left from my using copy and paste from the testing phase.</div><div>Most things if used only once I avoid using memory by calling it a name</div>
<div>11) Your calculation of the inertial tensor is nice, but I would again<br>recommend that you avoid creating four temporary arrays when<br>calculating the magnitude of 'ax'.</div><div>- Not changed, since this is the fastest way I can think of. Since I </div>
<div>realized the memory requirement of particles is much less than fields, </div><div>memory should not be the main concern here but speed. But I'm all</div><div>ears if you or Stephen have a faster/efficient way.</div>
<div>12) I would encourage you to investigate the eigenvector code if you<br>suspect it gives LH system coordinates.</div><div>- Not done<br>13) "for iter in range(0, maxiterations)") should just be for iter in<br>
xrange(maxiterations). There is no need to specify 0.</div><div>- Done<br>14) Please use print without parenthesis. Python 3 has transformed<br>print into a function, but the formatting rules are also different.</div><div>
- Done<br>15) Avoid using where!</div><div>- Done<br>16) use .size instead of len() for arrays</div><div>- Done<br>17) Feel free to throw an error. Define a subclass of Exception and<br>then raise an instance of it.</div>
<div>- Will look into how to write an error class in YT..., very useful I think but want to see if the changes I have so far is adequate.<br><br></div><div>updated script:</div><div><div><a href="http://paste.enzotools.org/show/1799/">http://paste.enzotools.org/show/1799/</a></div>
</div><div><br></div><div>So the ones I have trouble with are #6, 10, 11, all involving the phrase "in place". I had to wiki the phrase because I'm not familiar with it, and I think you mean "input is overwritten by output" as per <a href="http://en.wikipedia.org/wiki/In-place_algorithm">http://en.wikipedia.org/wiki/In-place_algorithm</a></div>
<div><br></div><div>When I solve for the 3 different cases of boundary condition, I had to double (1/3 axes * 3 cases = 1 extra) the particle position memory, because I only check 1 axis at a time, so it's not as bad as if I do all 3 axes and 3 cases at once. I can paste the code bit of "cases" into the next line where it is used twice by na.choose(). That should make the code "in place" but it would be a lot less readable. And since "cases" is used twice, having it in memory also makes things easier. I just wanted to clarify why it is left the way it is, but if you still think the in place method is better, I can make the change easily.</div>
<div><br></div><div>Another point which I didn't understand was the difference between using array of array compare to list of arrays. From what I read</div><div><a href="http://stackoverflow.com/questions/176011/python-list-vs-array-when-to-use">http://stackoverflow.com/questions/176011/python-list-vs-array-when-to-use</a></div>
<div>people seem to prefer list over array, and use arrays only when the same data type is absolutely require. Another point is that arrays are more memory efficient, because lists uses extra byte incase we want to append things to it etc, but if we want to do math to it or modify it, it is better use lists. In my case I multiply the rotation matrix to the xys positions, so I should use list because of doing math. Is that the correct understanding?</div>
<div><br></div><div>From</div><div>G.S.</div><div><br><div class="gmail_quote">---------- Forwarded message ----------<br>From: <b class="gmail_sendername">Matthew Turk</b> <span dir="ltr"><<a href="mailto:matthewturk@gmail.com" target="_blank">matthewturk@gmail.com</a>></span><br>
Date: Tue, Sep 6, 2011 at 4:11 AM<br>Subject: Re: [Yt-dev] ellipsoid part 2<br>To: <a href="mailto:yt-dev@lists.spacepope.org" target="_blank">yt-dev@lists.spacepope.org</a><br><br><br>Hi Geoffrey,<br>
<br>
Great to hear from you again on this project! I have several comments<br>
about the code interspersed; they are all relatively minor, but there<br>
are several.<br>
<div><br>
On Mon, Sep 5, 2011 at 10:23 PM, Geoffrey So <<a href="mailto:gsiisg@gmail.com" target="_blank">gsiisg@gmail.com</a>> wrote:<br>
> Hello Stephen,<br>
> This is part 2 of the ellipsoid parameter finding. The method is outlined<br>
> mainly in the Dubinski and Carlberg 1991 paper. But since we are doing this<br>
> in YT and we already know the particle list to begin with, I put in the<br>
> option that the user can specify if they want all the particles in the<br>
> region to be taken into account or only the particles found by the halo<br>
> finder algorithm. This is an option that I think should be amended to the<br>
> first ellipsoid parameter finder.<br>
<br>
</div>Is this the parameter halo_only? I would like to note that<br>
creation_time < 0 is not reliable, and will make your code<br>
Enzo-specific. I recommend for now using the boolean selector IsStar:<br>
<br>
dm_index = dd["IsStar"]<br>
<br>
Additionally, I would recommend not constructing a new array (which<br>
will be slower than allowing the array selections to be left as is)<br>
and to instead set ax (which should be renamed, as 'ax' is<br>
idiomatically used in yt to reference an axis) to a list.<br>
<br>
I would go a bit further and suggest you not use the parameter<br>
halo_only as it is implicit in the halo finding. If the halo finding<br>
subselected on dark matter, then you would want that in the ellipsoid.<br>
If it didn't subselect, then you probably also don't want to<br>
subselect in the ellipsoid finding.<br>
<div><br>
> The script form is linked here:<br>
> <a href="http://paste.enzotools.org/show/1782/" target="_blank">http://paste.enzotools.org/show/1782/</a><br>
> Once I have written out the methodology for myself, I'll work on the YT<br>
> documentation, will probably need some help in this area when it comes to<br>
> it.<br>
<br>
</div>I have a few comments:<br>
<br>
1) Use booleans for parameters such as use_halos. True and False, not 0 and 1.<br>
2) Iterate over the halos:<br>
<br>
for halo in halos:<br>
3) Index variables should be clearly index variables; "halo" is an<br>
object, but you use it as an index.<br>
4) cm => com<br>
5) You set ax to be an array twice, which will double the memory<br>
requirements. I don't think ax needs to be an array at all.<br>
6) Line 109, just modify ax in place. Fancy-indexing in numpy will (I<br>
believe) implicitly make a copy when you modify:<br>
<br>
In [1]: import numpy as na<br>
In [2]: arr = na.random.random(100)<br>
In [3]: ai = arr > 0.8<br>
In [4]: bb = arr[ai]<br>
In [5]: bb -= 1.0<br>
In [6]: arr.max()<br>
Out[6]: 0.99642788644174751<br>
In [7]: bb.max()<br>
Out[7]: -0.0035721135582524877<br>
<br>
7) Line 113, use range or xrange, and use singular "axis" not "axes"<br>
8) Line 114, you double the memory by making a new array again.<br>
9) Line 117, abs is a numpy function.<br>
10) Line 121, this makes four new arrays, each of which is the length<br>
of the number of particles. Construct a temporary array and do<br>
in-place math. You don't need "where" here; use booleans.<br>
11) Your calculation of the inertial tensor is nice, but I would again<br>
recommend that you avoid creating four temporary arrays when<br>
calculating the magnitude of 'ax'.<br>
12) I would encourage you to investigate the eigenvector code if you<br>
suspect it gives LH system coordinates.<br>
13) "for iter in range(0, maxiterations)") should just be for iter in<br>
xrange(maxiterations). There is no need to specify 0.<br>
14) Please use print without parenthesis. Python 3 has transformed<br>
print into a function, but the formatting rules are also different.<br>
15) Avoid using where!<br>
16) use .size instead of len() for arrays<br>
17) Feel free to throw an error. Define a subclass of Exception and<br>
then raise an instance of it.<br>
<br>
One fun thing would also be to plot projections of "Ones" using these<br>
ellipsoids as sources. That would help us see them better, too!<br>
<div><br>
> Doron Lemez has used this method in his paper and helped me tested the<br>
> accuracy of this script. I was wondering do I need to ask him if he wants<br>
> his name in the credits, or just mention his name as a special thanks?<br>
<br>
</div>Which credits do you mean? The CREDITS file is reserved for code<br>
contributions. Feel free in your documentation to mention that he<br>
assisted you, if you feel it is warranted. Please be sure that your<br>
documentation includes the ".. sectionauthor" directive to indicate<br>
that you wrote it.<br>
<br>
Please issue a pull request once you have included this routine into<br>
the halo finding module; we can directly review it then.<br>
<br>
Additionally, for the documentation modifications, I would encourage<br>
you to fork yt-doc and make your changes in source/analysis_modules/<br>
under the halo finding section. Include a few (small file size)<br>
images, too.<br>
<br>
Thanks, GS!<br>
<br>
-Matt<br>
<div><br>
> One<br>
> thing he did in the paper that was especially useful is defined the starting<br>
> sphere to have 2x the virial radius. Without that specification, I run into<br>
> non-convergence sometimes, where after many iterations I end up with no<br>
> particles. And I want to mention again that there's no guarantee that any<br>
> or all the haloes will converge with this method.<br>
> I tried to adhere to the PEP8 formatting, but old habits die hard so I'm<br>
> sure I made mistakes in the formatting, let me know if you see them.<br>
> From<br>
> G.S.<br>
</div>> _______________________________________________<br>
> Yt-dev mailing list<br>
> <a href="mailto:Yt-dev@lists.spacepope.org" target="_blank">Yt-dev@lists.spacepope.org</a><br>
> <a href="http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org" target="_blank">http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org</a><br>
><br>
><br>
_______________________________________________<br>
Yt-dev mailing list<br>
<a href="mailto:Yt-dev@lists.spacepope.org" target="_blank">Yt-dev@lists.spacepope.org</a><br>
<a href="http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org" target="_blank">http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org</a><br>
</div><br>
</div>