I have trouble reusing code (made for Illustris-1) for extracting parts of the merger tree. I identified the problem as code not reading the correct subhalo as root of the merger tree.
The code in question is lightly modified your code in sublink.py, more precisely the code of the function numMergers. That function have index=0 as starting point for reading the tree, so rootID = tree['SubhaloID'][index], in other words we find the correct root in the tree (provided by loadTree with snapshot number and subfindID) through root index which is always 0 in singular tree. From that starting point we can explore the length and breath of the tree.
This works flawlessly when used on the Illustris sublink tree, but completely fails in TNG300 sublink tree, and I have no idea why.
This is what I get for several roots in Illustris (values are root snapshot number, index of the root, rootID and first progenitor ID):
and this is what I get for several roots in TNG300 (the same values):
I load tree with function: tree=il.sublink.loadTree(basePath,99,subID,fields=fields,onlyMPB=False) but root is never on the snapshot 99 (and it's always on snapshot 135 in Illustris, which is generated with the same code), here it's on 43. and 26. snapshot, and mass of that subhalos are not even close to the correct roots which should be on 99. snapshot. It's as if index 0 is not the root. And if it's not, how can I find the root?
Can you help me solve this problem please?
Can you provide a specific example, what subID are you looking at for TNG300-1 snapshot 99?
These are first 36 subhalos of my sample: 11,12,14,15,16,17,18,19,20,21,22,23,25,26,27,29,30,32,34,36,37,39,41,43,11755,11756,11757,11759,11760,11762,11763,11764,11765,11767,11768,11769. The subhalos 11 and 12 are the ones from the picture above. After first 36 subhalos my program doesn't execute because of the (unknown) index error:
IndexError: index 2 is out of bounds for axis 0 with size 2
Not one of the first 36 of the sample have root at index 0.
I'm afraid I cannot debug your program, as I don't know the code. Can we identify an issue in a existing (or simple) piece of code? For now, I can load the tree and count the mergers of subhalo 11 without issue:
In : basePath = 'sims.TNG/TNG300-1/output/'
In : import illustris_python as il
In : tree = il.sublink.loadTree(basePath, 99, 11)
In : tree['SnapNum']
Out: array([99, 98, 97, ..., 15, 46, 95], dtype=int16)
In : tree['SubhaloID']
Out: array([1104150, 1104151, 1104152, ..., 1116695, 1116696, 1116697])
In : il.sublink.numMergers(tree)
Now my colleague tells me that simple using the code:
gives her nonsense, partial tree not starting from snapshot 99. It seems to me like a same problem.
Was there changes in illustris_python package version from Illustris-1?
You should certainly make sure you are using the current version of illustris_python. As you can see above, the tree loaded starts at snapshot 99. So to identify a problem, we need a case which produces an unexpected error, or unexpected behavior.
Uh...this is difficult because literaly when I execute this code:
sys.path = ["/data_4tb/tools"] + sys.path
import illustris_python as il
basePath = "/data_4tb/tng300/TNG300-1/postprocessing/"
fields = ['SubhaloID','NextProgenitorID','MainLeafProgenitorID','FirstProgenitorID','SnapNum','SubfindID']
I get this: [43 42 41 ..., 29 28 27]
I'm kind of lost here. I updated illustris_python library, but as I expected it didn't change anything.
Is it possible that merger tree is just partially downloaded, but checksums show everything ok?
Some days ago I asked you if there is some list where I could check manually if lengths of my tree parts are correct.
Because most of my tree parts are 12.2GB long, but about fifth of the parts are shorter (12GB, 12.1GB, 10.2GB... even 2.1GB for the 124. part) which seems a bit strange to me.
Checksums shows everything ok though.
Yes either the files are bad, or the code differs. Have you modified the illustris_python package?
Also, I do get an error running your exact code above:
TypeError: only integer scalar arrays can be converted to a scalar index
493 for use with the hyperslab selection routines
--> 495 start, stop, step = exp.indices(length)
496 # Now if step > 0, then start and stop are in [0, length];
497 # if step < 0, they are in [-1, length - 1] (Python 2.6b2 and later;
using python 3.8 and h5py 2.10.0. Are you using current versions? To fix it I have to instead set subID = 11 without the list.
subID = 11
You could try anything, if in doubt, on the Lab service, where the data files are the original ones.
I just wanted to update discussion in case anybody else have the same problem.
It was damaged offsets file. As soon as I replaced it, everything started working perfectly.