<div dir="ltr"><div><div><div><div><div><div><div><div>Hi Nathan,<br><br></div>You're right, the artifacts I saw were related to pixelization. When I do as you suggest:<br><br></div><div>ad = ds.all_data()<br></div><div>my_cut_region = ad.cut_region('(obj["temperatu<wbr>re"] > 1e4) & (obj["temperature"] < 1e6)')<br></div><div><div>slice = ds.slice(axis, coort, data_source=my_cut_region)</div><div>temp = slice['temperature']</div>print (temp.min(), temp.max())<br><br>Here's the result:<br>(10002.1923828 K, 999992.375 K)<br><br></div>So, the temperature values are all within the cut_region. Thanks for the suggestion! <br><br></div>As for the morton indices, I was surprised that they were all identical as well. After playing around with it, I discovered that I there's a difference in how straight Python handles types versus Cython: in particular, the code in my Python function pos_to_morton() used to restrict the input number to a valid range resulted in overflows consistently (hence the identical values). The line in question is:<br><br>ii[i] = np.min([(1 << ORDER_MAX)-1, np.max([0, ii[i]])])<br><br></div>I used this line instead of i64clip() from fp_utils.pyx. If I use this:<br><br>ii[i] &= np.uint64((1 << ORDER_MAX) - 1)<br><br></div>I get unique morton indices for the entire domain with the yt algorithm. So, the implementation in yt is probably fine, assuming that Cython handles the variable types appropriately in i64clip(). Although, It might make sense to move to Baert's algorithm in future for super high resolution simulations, since he uses 21 bits, vs. 20 bits (ORDER_MAX in yt/utilities/lib/geometry_lib.pyx). <br><br></div>My apologies for the wild goose chase(s), and thanks so much for the help!<br><br></div>Cheers,<br></div>Jason<br><div><div><div><div><div><div><br></div></div></div></div></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Oct 6, 2016 at 11:26 AM, Nathan Goldbaum <span dir="ltr"><<a href="mailto:nathan12343@gmail.com" target="_blank">nathan12343@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi Jason,<div><br></div><div>I'm having some trouble digesting your last message. Unfortunately it's hard for me to say exactly what's happening because I can't run the tests you've described locally. It would help me to discuss the issues you're having if you can share a test script that I can run locally. You can share it via the yt pastebin service to share it.</div><div><br></div><div>If you can't share the dataset you're using, you can use one of the publicly available test datasets on <a href="http://yt-project.org/data" target="_blank">yt-project.org/data</a>.</div><div><br></div><div>If you can share the dataset you're using, you might find the yt curldrop to be useful for sharing it: <a href="https://docs.hub.yt/services.html#curldrop" target="_blank">https://docs.hub.yt/servic<wbr>es.html#curldrop</a></div><div><br></div><div>Just to speculate, I bet the artifacts you're seeing are due to the pixelization rather than the data selection. One way to verify would be to create a slice data object that uses the cut_region you've defined as a data source and then query the temperature field and get the extrema of that data. Something like:</div><div><br></div><div>slice = ds.slice(axis, coort, data_source=my_cut_region)</div><div>temp = slice['temperature']</div><div>print (temp.min(), temp.max())</div><div><br></div><div>If it is due to the pixelization, it might be tricky to fix, but you should be able to at least set the colorbar range in your plots manually so the issue is less visible.</div><div><br></div><div>I'm surprised that the morton indices you're calculating are all identical given the fact that the CutRegionSelector seems to be mostly working, since the CutRegionSelector depends on all cells having unique Morton indices. Keep in mind that position_to_morton was not written with the intention of being part of the public API, so it's not really meant to be used in the way you're using it, I think.</div><span class="HOEnZb"><font color="#888888"><div><br></div><div>-Nathan</div></font></span><div><div class="h5"><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Oct 6, 2016 at 10:12 AM, Jason Galyardt <span dir="ltr"><<a href="mailto:jason.galyardt@gmail.com" target="_blank">jason.galyardt@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div><div><div><div><div><div><div>Hi Nathan,<br><br></div>The solution you came up with for <a href="https://bitbucket.org/yt_analysis/yt/issues/941/cutregionselector-doesnt-work-with-some" target="_blank">issue 941</a> that uses the morton indices + a tolerance factor seems to mostly work for my FLASH simulations for both single conditionals, as well as the union of 2 conditionals. The attached slice plot in temperature was created from one of my FLASH simulations with the following conditionals: 1e4 K < T < 1e6 K. It looks like there may be some small residual issues: there are a few cells here and there that seem to have T ~= 10**-9 K, which is odd, since T > 1e4 K was one of the conditionals. <br><br></div>In a related but separate issue, I came across a problem with the position_to_morton() function in geometry_utils.pyx. I looked at the morton indices for my FLASH simulation: I get identical values for all cells. I think this is basically a loss of precision, though I've not investigated fully. Since the position_to_morton() function uses magic numbers from a <a href="http://stackoverflow.com/questions/1024754/how-to-compute-a-3d-morton-number-interleave-the-bits-of-3-ints" target="_blank">StackOverflow post</a>, I thought I'd search for an alternative. I found one on <a href="http://www.forceflow.be/2013/10/07/morton-encodingdecoding-through-bit-interleaving-implementations/" target="_blank">Jeroen Baert's blog</a>. This magic bit implementation spreads the bits out differently and restricts the input to the first 21 bits in each dimension (so the morton index for each cell is a 63 bit number). I wrote a very pedestrian python script to test both of these implementations against my FLASH simulation:<br><br><a href="http://paste.yt-project.org/show/6844/" target="_blank">http://paste.yt-project.org/sh<wbr>ow/6844/</a><br><br>I used two different values of the cell width (delta below and in the script above): one uses the same width as is used in yt/utilities/lib/geometry_util<wbr>s.pyx:<br><br>delta[i] = (dle[i] - dre[i]) / (1 << ORDER_MAX)<br><br></div>where dle is the domain left edge array, and dre is the domain right edge array, and ORDER_MAX = 20. The other cell width is the minimum cell width in each dimension. I found that using the yt implementation of cell width (the ORDER_MAX one), all the morton indices are identical for the yt implementation of position_to_morton() used in yt. Whereas, the yt implementation of position_to_morton() works fine with the cell width set to the minimum found in each dimension (i.e. all morton indices are unique).<br><br></div>Jeroen Baert's morton index generation algorithm generates unique morton indices for both implementations of the cell width array. The output of this script when run on my FLASH simulations is below. <br><br></div>Using Jeroen Baert's morton index generation algorithm doesn't seem to affect the issue of the cut_region artifacts, as far as I can tell. That is to say, the artifacts are still there when you don't include the tolerance factor code. I've tested this by inserting Jeroen Baert's algorithm into geometry_utils.pyx from Nathan's latest changeset (5cf30893afdf), rebuilding, and then generating the same slice plot as attached. The plots are identical to the eye (down to the small regions of T ~ 10**-9 K). <br><br>So, I'm not sure what this all means because I don't have a clear picture of how the morton indices are used in yt. I do find it odd that the hash used in the cut_region() code still works when all the morton indices are identical though. This is probably not terribly helpful; any ideas on how we can shed more light on either of these issues? Anything else I can do?<br><br></div><div><br></div>Cheers,<br></div>Jason<br><br>------<br></div>Jason Galyardt<br></div>University of Georgia<br><br><div><div><br><div><div><div><div><div><div><div><div><div>Here's the output of the script, as promised:<br>------------------------------<wbr>----------------------------<br><br>yt : [INFO ] 2016-10-06 10:05:05,130 Particle file found: smithCloud_hdf5_part_0115<br>yt : [INFO ] 2016-10-06 10:05:05,137 integer runtime parameter checkpointfilenumber overwrites a simulation scalar of the same name<br>yt : [INFO ] 2016-10-06 10:05:05,137 integer runtime parameter particlefilenumber overwrites a simulation scalar of the same name<br>yt : [INFO ] 2016-10-06 10:05:05,138 integer runtime parameter plotfilenumber overwrites a simulation scalar of the same name<br>yt : [INFO ] 2016-10-06 10:05:05,151 Parameters: current_time = 1.81358087048e+15<br>yt : [INFO ] 2016-10-06 10:05:05,151 Parameters: domain_dimensions = [16 16 96]<br>yt : [INFO ] 2016-10-06 10:05:05,151 Parameters: domain_left_edge = [ -6.17136050e+21 -6.17136050e+21 -6.78849660e+22]<br>yt : [INFO ] 2016-10-06 10:05:05,151 Parameters: domain_right_edge = [ 6.17136050e+21 6.17136050e+21 6.17136050e+21]<br>yt : [INFO ] 2016-10-06 10:05:05,152 Parameters: cosmological_simulation = 0.0<br>level # grids # cells # cells^3<br>------------------------------<wbr>----------------<br> 0 48 24576 30<br> 1 288 147456 53<br> 2 1856 950272 99<br> 3 10776 5517312 177<br>------------------------------<wbr>----------------<br> 12968 6639616<br><br><br>t = 1.81358087e+15 = 1.81358087e+15 s = 5.74689099e+07 years<br><br>Smallest Cell:<br> Width: 3.125e-05 Mpc<br> Width: 3.125e+01 pc<br> Width: 6.446e+06 AU<br> Width: 9.643e+19 cm<br>len(ytMorton1): 5812736<br><br>Using delta = [min(dx), min(dy), min(dz)]:<br>yt Morton: number of identical values: 0<br>Baert Morton: number of identical values: 0<br><br>Using delta[i] = (dle[i] - dre[i]) / (1 << ORDER_MAX):<br>yt Morton: number of identical values: 5812735<br>Baert Morton: number of identical values: 0<br><br><br></div><div><br></div></div></div></div></div></div></div></div></div></div></div></div><div class="m_-7828266142145116872m_1265035036860908473HOEnZb"><div class="m_-7828266142145116872m_1265035036860908473h5"><div class="gmail_extra"><br><div class="gmail_quote">On Sat, Oct 1, 2016 at 8:05 PM, Nathan Goldbaum <span dir="ltr"><<a href="mailto:nathan12343@gmail.com" target="_blank">nathan12343@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi Jason,<div><br></div><div>I'm CCing yt-dev here so that this project discussion happens in public. I hope that's ok.</div><div><br></div><div>Yes, it would be great if you could take a stab at this. Unfortunately the main blocker here is choosing a good, fast, periodic KDTree that can be used easily from within cython. Having this available would be great for many other things besides the CutRegionSelection. We in fact already have a copy of scipy's C KDTree implementation vendored inside yt that we'd like to get rid of. I believe there are plans to either add a dependency on or vendor a fast cython kdtree that Meagan Lang has been working on:</div><div><br></div><div><a href="https://bitbucket.org/langmm/cykdtree" target="_blank">https://bitbucket.org/langmm/c<wbr>ykdtree</a></div><div><br></div><div>Earlier today I took a stab at implementing a naive O(N^2) search in CutRegionSelector.select_cell. Unfortunately it's very, very slow to do this so we're going to need some sort of acceleration in order to realistically accomplish this, either via .a tree approach or by using a hash search like we're currently doing via python's builtin set object.</div><div><br></div><div>An alternate way around this that maintains the O(1) speed of the current hashmap-based approach would be to compute unique integer indexes for each cell inside the cut_region and save those instead of just saving the floating point position of each cell. When I originally wrote the CutRegionSelector there wasn't anything like that in yt, but since then we've added a way to compute morton indices for each cell, which should do the trick.</div><div><br></div><div>In fact, I went ahead and took a stab at this after coming up with that idea while writing this e-mail:</div><div><br></div><div><a href="https://bitbucket.org/yt_analysis/yt/pull-requests/2405" target="_blank">https://bitbucket.org/yt_analy<wbr>sis/yt/pull-requests/2405</a></div><div><br></div><div>And unfortunately it still doesn't work perfectly - I still see artifacts (although they're less severe) in the image produced by the test script in the issue.</div><div><br></div><div>I think my attempt in PR 2405 might be going on the right track - you could try to poke around with it, or test it with your data.</div><div><br></div><div>Sorry to not have a satisfactory answer here. I was very excited about initially being able to implement cut_region data sources for SlicePlots, but the issue with FLASH data has been an annoying thorn in my side for a while, it would be great if we could figure out a way to fix it once and for all.</div><span class="m_-7828266142145116872m_1265035036860908473m_-8138335024614259689HOEnZb"><font color="#888888"><div><br></div><div>-Nathan</div></font></span><div><div class="m_-7828266142145116872m_1265035036860908473m_-8138335024614259689h5"><div><br><div class="gmail_extra"><div class="gmail_quote">On Sat, Oct 1, 2016 at 4:27 PM, Jason Galyardt <span dir="ltr"><<a href="mailto:jason.galyardt@gmail.com" target="_blank">jason.galyardt@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div><div>Hi Nathan,<br><br></div>Unfortunately, contours won't be enough for me -- I need the actual data region to do analysis on the gas contained within (i.e. more than visualization). I'm not very familiar with the yt codebase, but I'd like to see whether it's worth my time to work on this issue. I'm not very familiar with the yt codebase or Cython, but my programming skills are decent; how much work do you think this would take? Could I use an existing KDTree library, such as the one in SciPy or <a href="https://github.com/jtsiomb/kdtree" target="_blank">this C implementation</a><code>? <br><br></code></div><code>Thanks,<br></code></div><code>Jason<br></code><div><div><code><br></code></div></div></div><div class="m_-7828266142145116872m_1265035036860908473m_-8138335024614259689m_4337933984561522302gmail-HOEnZb"><div class="m_-7828266142145116872m_1265035036860908473m_-8138335024614259689m_4337933984561522302gmail-h5"><div class="gmail_extra"><br><div class="gmail_quote">On Sat, Oct 1, 2016 at 3:44 PM, Nathan Goldbaum <span dir="ltr"><<a href="mailto:nathan12343@gmail.com" target="_blank">nathan12343@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Hi Jason,<div><br></div><div>Unfortunately no. I think I might be able to come back to this issue soonish, we're going to be introducing a new fast periodic KDTree which could be exploited here.</div><div><br></div><div>For now I think you're<span></span> going to need to make this plot in a different way, without using a cut_region. Perhaps by overplotting contours?</div><span class="m_-7828266142145116872m_1265035036860908473m_-8138335024614259689m_4337933984561522302gmail-m_-2056896066610442086HOEnZb"><font color="#888888"><div><br></div></font></span><div><span class="m_-7828266142145116872m_1265035036860908473m_-8138335024614259689m_4337933984561522302gmail-m_-2056896066610442086HOEnZb"><font color="#888888">Nathan</font></span><div><div class="m_-7828266142145116872m_1265035036860908473m_-8138335024614259689m_4337933984561522302gmail-m_-2056896066610442086h5"><br><br>On Saturday, October 1, 2016, Jason Galyardt <<a href="mailto:jason.galyardt@gmail.com" target="_blank">jason.galyardt@gmail.com</a>> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Hi Nathan,<br><br></div><div>Is there a previous version of yt in which compound cut regions work?<br><br></div><div>Thanks,<br></div><div>Jason<br></div><div><div><div><br>-----<br><br><div dir="ltr">Hi Jason,<div><br></div><div>This is a known issue:</div><div><br></div><div><a href="https://bitbucket.org/yt_analysis/yt/issues/941/cutregionselector-doesnt-work-with-some" target="_blank">https://bitbucket.org/yt_analy<wbr>sis/yt/issues/941/cutregionsel<wbr>ector-doesnt-work-with-some</a><br></div><div><br></div><div>Unfortunately
fixing it will require a modification to the cut_region implementation
to do a better job of searching for cells that belong to the region for
chained selectors.</div><div><br></div><div>-Nathan</div></div><br>On Fri, Sep 30, 2016 at 1:25 PM, Jason Galyardt <span dir="ltr"><<a>jason.galyardt@gmail.com</a>></span> wrote:<br><div><div><div><div><div>Dear all,<br><br></div>I
would like to generate a compound cut region for a FLASH simulation,
but the output I'm getting is odd. I've pasted my simple script here:<br><br><a href="http://paste.yt-project.org/show/6839/" target="_blank">http://paste.yt-project.org/sh<wbr>ow/6839/</a><br><br></div>It
loads a simulation file, creates a compound cut_region (1e4 K < T
< 1e6 K), and applies it to a slice plot. The resulting image is
banded, as you can see by the attached plots. Am I missing something, or
is this a problem related to yt?<br><br></div>FYI, I get the same
behavior when I use a projection plot or only one of the conditionals
(e.g. T < 1e6 K). I've tried both the stock 3.3.1 version, as well as
the self-contained development version that I installed earlier this
afternoon using install_script.sh. Let me know if any other info is
useful.<br><br></div>Cheers,<br></div>Jason<br><br></div></div></div></div>
</blockquote></div></div></div>
</blockquote></div><br></div>
</div></div></blockquote></div><br></div></div></div></div></div>
</blockquote></div><br></div>
</div></div></blockquote></div><br></div></div></div></div>
</blockquote></div><br></div>