Indexerror within make_photon function (reopened)

[Rongjun-ANU]

Sorry, the previous one was closed by accident

Hi, I meet an index error when using make_photon function in pyxsim

I have loaded a yt object called ds_plt1555000; the shape is 1281281024 I used ad_plt1555000 = ds_plt1555000.all_data()

Here is the issue:

sp_plt1555000 = ds_plt1555000.sphere("c", (10, "kpc"))

Generate the photon list
n_photons, n_cells = pyxsim.make_photons("plt1555000_photon_list",
data_source=sp_plt1555000,
redshift=redshift, area=area,
exp_time=exp_time, source_model=source_model)

pyxsim : [INFO ] 2023-09-26 12:19:14,484 Cosmology: h = 0.71, omega_matter = 0.27, omega_lambda = 0.73
pyxsim : [INFO ] 2023-09-26 12:19:14,486 Using emission measure field '('gas', 'emission_measure')'.
pyxsim : [INFO ] 2023-09-26 12:19:14,486 Using temperature field '('gas', 'temperature')'.
Preparing spectrum table : 100%
138/138 [00:37<00:00, 81.02it/s]
Processing cells/particles : 100%
16777216/16777216 [03:33<00:00, 92569.44it/s]

---

IndexError Traceback (most recent call last)
/home/rongjun/VSCAA/ASTR8010/try/Dash.ipynb Cell 59 line 1
7 region = ds_plt1555000.region(domain_center, domain_left_edge, domain_right_edge)
8 # ds_p = ds_plt1555000.force_periodicity()
9
10 # Generate the photon list
---> 11 n_photons, n_cells = pyxsim.make_photons("plt1555000_photon_list",
12 data_source=sp_plt1555000,
13 redshift=redshift, area=area,
14 exp_time=exp_time, source_model=source_model)

File ~/.local/lib/python3.8/site-packages/pyxsim/photon_list.py:416, in make_photons(photon_prefix, data_source, redshift, area, exp_time, source_model, point_sources, parameters, center, dist, cosmology, velocity_fields, bulk_velocity, observer, fields_to_keep)
413 d["energy"].resize((p_size,))
415 for i, ax in enumerate("xyz"):
--> 416 pos = chunk[p_fields[i]][idxs].to_value("kpc")
417 # Fix photon coordinates for regions crossing a periodic boundary
418 if ds.periodicity[i]:

File ~/.local/lib/python3.8/site-packages/unyt/array.py:1722, in unyt_array.getitem(self, item)
1721 def getitem(self, item):
-> 1722 ret = super(unyt_array, self).getitem(item)
1723 if getattr(ret, "shape", None) == ():
1724 ret = unyt_quantity(ret, bypass_validation=True, name=self.name)

IndexError: index 1189003 is out of bounds for axis 0 with size 1048576

/home/rongjun/.local/lib/python3.8/site-packages/unyt/array.py(1722)getitem()
1720
1721 def getitem(self, item):
-> 1722 ret = super(unyt_array, self).getitem(item)
1723 if getattr(ret, "shape", None) == ():
1724 ret = unyt_quantity(ret, bypass_validation=True, name=self.name)

---

In pdb, I navigate to "--> 416 pos = chunk[p_fields[i]][idxs].to_value("kpc")"
I can see that chunk.shape' is (16777216,). That is good because 128*128*1024=16777216. But chunk[p_fields[i]]' is actually only (1048576,)
And `p_fields[i]' is ('index', 'x') as i=0 here
So I don't know why the size of ('index', 'x') field reduces from 16777216 to 1048576 in make_photon function.

I also see that before using make_photon, sp_plt1555000[('index', 'x')].shape is still (16777217,), but after make_photon, it reduces to (1048576,)

I also found that, idxs will change once I rerun the make_photon function. So every time I use make_photon, it generates a different list of indexes and there is always some index greater than 1048576, and this returns indexerror.

Finally, Here is the link to my script: https://github.com/Rongjun-ANU/Indexerror/blob/main/Dash-Code.ipynb

[jzuhone]

Hi @Rongjun-ANU ,

So I looked at your notebook, and there's a lot going on in there. In particular, I'm noticing that you appear to be redefining fields in a number of places that should exist already, like ("gas", "mass"), and that the other issue is that you are creating a covering grid (which returns 3D NumPy arrays for fields) and re-defining field definitions for the dataset in terms of them.

Field definitions are just functions, so they need to be able to pass the data in whatever shape in comes in and do the required manipulations, preserving the shape of the arrays. Is there a particular reason you have to convert the data to 3D before using pyXSIM? Unless I am misunderstanding the use case, this should not be necessary in most circumstances.

[Rongjun-ANU]

Hi @jzuhone

1. Actually, fields like ("gas", "mass") are not included in my data initially. The on-disk fields of my data only include

   [('boxlib', 'gasDensity'), ('boxlib', 'gasEnergy'), ('boxlib', 'gasInternalEnergy'), ('boxlib', 'scalar_0'), ('boxlib', 'scalar_1'), ('boxlib', 'scalar_2'), ('boxlib', 'x-GasMomentum'), ('boxlib', 'y-GasMomentum'), ('boxlib', 'z-GasMomentum')]

   So that's the reason why I need to define this fields to satisfy the requirement of function make_source_fields and make_photon.

2. I use the covering grid as I want to have a 3d visualization of some fields. But even I rerun the code without covering gird (I.e., 1d array everywhere), it still raises indexxerror.
3. Also, I noticed that every time I rerun the code and see "index xxx is out of bounds for axis 0 with size 1048576", this index value "xxx" changes every time.

[Rongjun-ANU]

Here are the fields included in my data initially:

[('boxlib', 'cell_volume'),
 ('boxlib', 'dx'),
 ('boxlib', 'dy'),
 ('boxlib', 'dz'),
 ('boxlib', 'gasDensity'),
 ('boxlib', 'gasEnergy'),
 ('boxlib', 'gasInternalEnergy'),
 ('boxlib', 'path_element_x'),
 ('boxlib', 'path_element_y'),
 ('boxlib', 'path_element_z'),
 ('boxlib', 'scalar_0'),
 ('boxlib', 'scalar_1'),
 ('boxlib', 'scalar_2'),
 ('boxlib', 'volume'),
 ('boxlib', 'x'),
 ('boxlib', 'x-GasMomentum'),
 ('boxlib', 'y'),
 ('boxlib', 'y-GasMomentum'),
 ('boxlib', 'z'),
 ('boxlib', 'z-GasMomentum'),
 ('gas', 'cell_volume'),
 ('gas', 'dx'),
 ('gas', 'dy'),
 ('gas', 'dz'),
 ('gas', 'path_element_x'),
 ('gas', 'path_element_y'),
 ('gas', 'path_element_z'),
 ('gas', 'volume'),
 ('gas', 'x'),
 ('gas', 'y'),
 ('gas', 'z'),
 ('index', 'cell_volume'),
 ('index', 'cylindrical_radius'),
 ('index', 'cylindrical_theta'),
 ('index', 'cylindrical_z'),
 ('index', 'dx'),
 ('index', 'dy'),
 ('index', 'dz'),
 ('index', 'grid_indices'),
 ('index', 'grid_level'),
 ('index', 'morton_index'),
 ('index', 'ones'),
 ('index', 'ones_over_dx'),
 ('index', 'path_element_x'),
 ('index', 'path_element_y'),
 ('index', 'path_element_z'),
 ('index', 'radius'),
 ('index', 'spherical_phi'),
 ('index', 'spherical_radius'),
 ('index', 'spherical_theta'),
 ('index', 'virial_radius_fraction'),
 ('index', 'volume'),
 ('index', 'x'),
 ('index', 'y'),
 ('index', 'z'),
 ('index', 'zeros')]

[jzuhone]

@Rongjun-ANU can you post this dataset somewhere so I can play with it?

[Rongjun-ANU]

sure, https://drive.google.com/drive/folders/1Xf7nJ_cHu9e8h8kqCt9i09yKrFT1xjp4?usp=drive_link

[jzuhone]

@Rongjun-ANU just want to double-check internalEnergy is indeed in erg and not erg/cm**3?

[Rongjun-ANU]

Hi @jzuhone Sorry for the confusion. The original fields provided in my data are purely values without units. However, when I go to check the place where I compute the temp (temperature), I think the unit of internalEnergy should be erg/g.

Also, I think I have partially solved the indexerror.

I extract some fields (density, mass, temperature, emission_measure) from original object. Then set a new yt object in a new notebook with these fields and still keep the same properties (dim:1281281024, cell size dx=2.3578125e+19, 'cm', i.e. roughly 8pc, etc.). And finally, it works!

The reason why I say I just partially solve it is that I still cannot figure out why the previous codes are not working. However, I suspect the index issue arises when I try to add fields that are supposed to exist. For example, we do have a density field called ('boxlib', 'gasDensity'), but you know it is dimensionless. Also, some functions in pyxsim require a density field named ('gas', 'density'), so I have to extract ('boxlib', 'gasDensity') and times it with the unit (g/cm^3) and add this new yt array as a derived field ('gas', 'density') in my previous ipynb. And probably, such action may cause indexerror.

[jzuhone]

@Rongjun-ANU gasInternalEnergy cannot be erg/g because its value is too small for this system--it's more likely to be erg or erg/cm**3.

As I noted before, I would not redefine fields for this dataset using the covering grid extraction, because yt is not expecting that and it breaks a lot of things. I created a notebook that shows how to map your fields to the yt fields:

https://gist.github.com/jzuhone/26e984707be5052ec2c0409f810a91e6

Note that you get all of the fields you need to run pyXSIM in this way. pyXSIM accepts sphere, region, and disk objects from yt, but not covering grids, at least when making photons. You can use covering grids when making X-ray fields, however.

Let me know if you can make a photon list in this way, and let me know if you have any further questions.

[Rongjun-ANU]

Hi @jzuhone , thank you. I can make a photon list now.