jzuhone / pyxsim

Simulating X-ray observations from astrophysical sources.
http://hea-www.cfa.harvard.edu/~jzuhone/pyxsim
Other
20 stars 7 forks source link

make_photons TypeError: 'float' object is not subscriptable #35

Closed benopp99 closed 2 years ago

benopp99 commented 2 years ago

This happens after looping through several halos, and does not happen at the same place in the same halo if I change e.g. the exposure time (otherwise completely repeatable), and almost seems to reach an error after a similar amount of time. Attaching my script too at end.

Traceback (most recent call last): File "test_script_4.py", line 143, in n_photons, n_cells = pyxsim.make_photons(simput_basename + ".my_photons", my_sphere, redshift_project, area, simput_exp_time, thermal_model) File "/srv/conda/envs/notebook/lib/python3.7/site-packages/pyxsim/photon_list.py", line 291, in make_photons chunk_data = source_model(chunk) File "/srv/conda/envs/notebook/lib/python3.7/site-packages/pyxsim/source_models.py", line 308, in call kT_sorted = kT[idxs] TypeError: 'float' object is not subscriptable

Thanks, Ben

Script: import yt import numpy as np import pyxsim import soxs import sys import astropy import h5py

from illustris_python import groupcat as groups

from yt.extensions.astro_analysis.halo_analysis import HaloCatalog

def _metallicity(field, data): if((sim=="IllustrisTNG")): return data["PartType0","GFM_Metallicity"].to("Zsun") if(sim=="SIMBA"): return data["gas","metallicity"].to("Zsun")

def hot_gas(pfilter, data): pfilter1 = data[pfilter.filtered_type, "temperature"] > 2.0e5 pfilter2 = data[pfilter.filtered_type, "temperature"] < 4.5e8 pfilter3 = data["PartType0", "StarFormationRate"] == 0.0

if((sim=="IllustrisTNG")):
    pfilter4 = data["PartType0", "GFM_CoolingRate"] < 0.0
    return ((pfilter1 & pfilter2 & pfilter3) & pfilter4)
else:
    return ((pfilter1 & pfilter2 & pfilter3))

sim = "IllustrisTNG"

sim = "SIMBA"

sim = sys.argv[1]

run = "1P_50"

run = sys.argv[2] snap = "033" redshift_project = 0.01 axis = "z" exp_time_ksec = 3. exp_time = exp_time_ksec1e+03 simput_exp_time = exp_time10. area = 10000.0 sky_center = [0.0, 0.0] # in degrees

snapshotfilename = "Simulations/" + sim + "/" + run + "/snap" + snap + ".hdf5" halocat_filename = "Simulations/" + sim + "/" + run + "/fof_subhalotab" + snap + ".hdf5" output_basename = sim + "." + run + "." + snap

cat_filename = output_basename + ".halodata.cat" cat_file = open(cat_filename,"w") cat_file.write("#ID lg(M200c) xcoord ycoord zcoord\n") cat_file.close()

print(halocat_filename)

ds = yt.load(snapshot_filename) yt.add_particle_filter("hot_gas", function=hot_gas, filtered_type='gas', requires=["temperature"]) ds.add_particle_filter("hot_gas")

halo_file = h5py.File(halocat_filename,"r")

h = halo_file["Header"].attrs["HubbleParam"] redshift = halo_file["Header"].attrs["Redshift"] print("h, redshift= ",h,redshift)

M200c = halo_file["Group"]["Group_M_Crit200"] / h * 1e+10 nhalo = len(M200c)

Grp_Pos = halo_file["Group"]["GroupPos"] print("GrpPos= ", Grp_Pos,Grp_Pos[0],Grp_Pos[1],Grp_Pos[2])

Grp_FirstSub = halo_file["Group"]["GroupFirstSub"]

Sub_Pos = halo_file["Subhalo"]["SubhaloPos"]

Sub_Mass = halo_file["Subhalo"]["SubhaloMass"]

Sub_Len = halo_file["Subhalo"]["SubhaloLen"]

print("SubPos_GrpCens= ", Sub_Pos[Grp_FirstSub[0]],Sub_Pos[Grp_FirstSub[1]],Sub_Pos[Grp_FirstSub[2]])

halo = h5py.File(halocat_filename)

ad = ds.all_data()

print("ds.field_list= ", ds.field_list)

print("type(ds)= ", type(ds))

print("ds.fields.gas.C_metallicity= ", ds.fields.gas.C_metallicity)

print("ds.fields.gas.C_fraction= ", ds.fields.gas.C_fraction)

mass_centre = ad.quantities.center_of_mass()

print ("mass_centre= ", mass_centre)

for i in range(nhalo):

halo_id = i
if((M200c[i]<1e+12) | (M200c[i]>2e+12)):
    continue
else:
    print("lg(M200c)= ",np.log10(M200c[i]))
    print("Sub_Pos[Grp_FirstSub[i]]= ", Sub_Pos[Grp_FirstSub[i]])

cat_file = open(cat_filename,"a")
cat_file.write("%4d %5.2f %9.3f %9.3f %9.3f \n"%(i, np.log10(M200c[i]), Sub_Pos[Grp_FirstSub[i]][0], Sub_Pos[Grp_FirstSub[i]][1], Sub_Pos[Grp_FirstSub[i]][2]))
cat_file.close()

my_sphere = ds.sphere(Sub_Pos[Grp_FirstSub[i]], (1000., "kpc"))

Xray_basename = output_basename + "." + str(int(exp_time/1000)) + "ks.z" + str(redshift_project) + "." + axis + ".grp" + str(halo_id)
simput_basename = output_basename + "." + str(int(simput_exp_time/1000)) + "ks.z" + str(redshift_project) + "." + axis + ".grp" + str(halo_id)

do_simput = True
do_soxs = True

if(do_simput):

    emin = 0.02 # minimum energy in keV
    emax = 12.0 # maximum energy in keV
    nchan = 1200 # should ideally result in spectral resolution that is better than the detector you are using

    kT_min = 0.025 # minimum temperature in keV
    kT_max = 64.0 # maximum temperature in keV
    n_kT = 10000 # number of temperature bins

    print("ds.fields.gas.metallicity= ", ds.fields.gas.metallicity) 

    var_elem = {"C": ("hot_gas","C_fraction"),"N": ("hot_gas","N_fraction"),"O": ("hot_gas","O_fraction"),"Ne": ("hot_gas","Ne_fraction"),"Mg": ("hot_gas","Mg_fraction"),"Si": ("hot_gas","Si_fraction"), "S": ("hot_gas","S_fraction"), "Ca": ("hot_gas","Ca_fraction"),"Fe": ("hot_gas","Fe_fraction")}
    if(sim=="IllustrisTNG"):
        var_elem = {"C": ("hot_gas","C_fraction"),"N": ("hot_gas","N_fraction"),"O": ("hot_gas","O_fraction"),"Ne": ("hot_gas","Ne_fraction"),"Mg": ("hot_gas","Mg_fraction"),"Si": ("hot_gas","Si_fraction"),"Fe": ("hot_gas","Fe_fraction")}

    if(sim=="SIMBA"):
        var_elem = {"C": ("hot_gas","C_metallicity"), "N": ("hot_gas","N_metallicity"), "O": ("hot_gas","O_metallicity"), "Ne": ("hot_gas","Ne_metallicity"), "Mg": ("hot_gas","Mg_metallicity"), "Si": ("hot_gas","Si_metallicity"), "S": ("hot_gas","S_metallicity"), "Ca": ("hot_gas","Ca_metallicity"), "Fe": ("hot_gas","Fe_metallicity")}

    thermal_model = pyxsim.ThermalSourceModel("apec",emin,emax,nchan,kT_min=kT_min,kT_max=kT_max,n_kT=n_kT,max_density=5.0e-25,Zmet=("hot_gas", "metallicity"), var_elem=var_elem,emission_measure_field=("hot_gas","emission_measure"),temperature_field=("hot_gas","temperature"))
    ###thermal_model = pyxsim.ThermalSourceModel("apec",emin,emax,nchan,kT_min=kT_min,kT_max=kT_max,n_kT=n_kT,max_density=5.0e-25,Zmet=("hot_gas", "metallicity"),emission_measure_field=("hot_gas","emission_measure"),temperature_field=("hot_gas","temperature"))

    n_photons, n_cells = pyxsim.make_photons(simput_basename + ".my_photons", my_sphere, redshift_project, area, simput_exp_time, thermal_model)

    nH = 0.02 # 10^20 atoms/cm**2
    abs_model = "tbabs" #pyxsim.TBabsModel(nH)

    n_events = pyxsim.project_photons(simput_basename + ".my_photons", simput_basename + ".my_events", axis, sky_center,absorb_model=abs_model, nH=nH)
    events = pyxsim.EventList(simput_basename + ".my_events.h5")

    events.write_to_simput(simput_basename, overwrite=True)

    simput_file = "%s_simput.fits" % simput_basename

if(do_soxs):

soxs.instrument_simulator(simput_file, Xray_basename + ".hdxi_nobg_evt.fits", exp_time, "lynx_hdxi",sky_center, overwrite=True, foreground=False, ptsrc_bkgnd=False, instr_bkgnd=False)

soxs.write_image(Xray_basename + ".hdxi_nobg_evt.fits", Xray_basename + ".hdxi_nobg_img.fits", emin=0.5, emax=2.0, overwrite=True,reblock=4)

soxs.write_spectrum(Xray_basename + ".hdxi_nobg_evt.fits", Xray_basename + ".hdxi_spec.pha", overwrite=True)

fig,ax = soxs.plot_spectrum(Xray_basename + ".hdxi_spec.pha", xmin=0.2, xmax=7.0)

fig.savefig(Xray_basename + ".spec_hdxi_nobg.png")

    soxs.instrument_simulator(simput_file, Xray_basename + ".axis_nobg_evt.fits", exp_time, "axis",sky_center, overwrite=True, foreground=False, ptsrc_bkgnd=False, instr_bkgnd=False)
    soxs.write_image(Xray_basename + ".axis_nobg_evt.fits", Xray_basename + ".axis_nobg_img.fits", emin=0.5, emax=2.0, overwrite=True,reblock=4)
    soxs.write_spectrum(Xray_basename + ".axis_nobg_evt.fits", Xray_basename + ".axis_spec.pha", overwrite=True)
    fig,ax = soxs.plot_spectrum(Xray_basename + ".axis_spec.pha", xmin=0.2, xmax=7.0)
    fig.savefig(Xray_basename + ".spec_axis_nobg.png")

    soxs.instrument_simulator(simput_file, Xray_basename + ".acisi_nobg_evt.fits", exp_time, "chandra_acisi_cy0",sky_center, overwrite=True, foreground=False, ptsrc_bkgnd=False, instr_bkgnd=False)
    soxs.write_image(Xray_basename + ".acisi_nobg_evt.fits", Xray_basename + ".acisi_nobg_img.fits", emin=0.5, emax=2.0, overwrite=True, reblock=4)
    soxs.write_spectrum(Xray_basename + ".acisi_nobg_evt.fits", Xray_basename + ".acisi_spec.pha", overwrite=True)
    fig,ax = soxs.plot_spectrum(Xray_basename + ".acisi_spec.pha", xmin=0.2, xmax=7.0)
    fig.savefig(Xray_basename + ".spec_acisi_nobg.png")
jzuhone commented 2 years ago

See these lines here:

https://github.com/jzuhone/pyxsim/blob/b9cacebf01da2440db38584e16a79034425cf01c/pyxsim/source_models.py#L297-L304

Can you put a print("orig_ncells = ", orig_ncells) right after line 298? You may have to reinstall from source with pip install -e . to make it work.

It would appear that what you have is a situation where the number of particles in a chunk is 1, which would be very odd but possible, I guess.

benopp99 commented 2 years ago

Weird, I can't do a pip reinstall since I can't find a setup.py in the original pip install. So pip install -e . doesn't work.

ERROR: File "setup.py" not found. Directory cannot be installed in editable mode: /home/jovyan/.local/lib/python3.7/site-packages/pyxsim

This is a rather basic pip question, and I'm missing something (tried also pip install -e /srv/conda/envs/notebook/lib/python3.7/site-packages/pyxsim but this gave me same issue).

-Ben

On Tue, Aug 10, 2021 at 7:50 PM John ZuHone @.***> wrote:

See these lines here:

https://github.com/jzuhone/pyxsim/blob/b9cacebf01da2440db38584e16a79034425cf01c/pyxsim/source_models.py#L297-L304

Can you put a print("orig_ncells = ", orig_ncells) right after line 298? You may have to reinstall from source with pip install -e . to make it work.

It would appear that what you have is a situation where the number of particles in a chunk is 1, which would be very odd but possible, I guess.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/jzuhone/pyxsim/issues/35#issuecomment-896436389, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMEF72PLI6GS2DAQR5FS5CTT4HJNTANCNFSM5B5E5NJQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .

jzuhone commented 2 years ago

You need to download a new copy from source, you can't pip install from the one that pip installed originally.

benopp99 commented 2 years ago

First stoppage with printing orig_ncells. I don't see anything obvious.

Preparing spectrum table : 100%|█████████████████████████████████████████████████████████| 201/201 [00:41<00:00, 4.79it/s]

Processing cells/particles : 0%| | 0/18560 [00:00<?, ?it/s]orig_ncells = 0 Processing cells/particles : 0%| | 0/18560 [00:05<?, ?it/s]orig_ncells = 0 Processing cells/particles : 0%| | 0/18560 [00:10<?, ?it/s]orig_ncells = 0 Processing cells/particles : 0%| | 0/18560 [00:15<?, ?it/s]orig_ncells = 6549 Processing cells/particles : 35%|█████████████████▋ | 6549/18560 [00:49<00:06, 1838.06it/s]orig_ncells = 0 Processing cells/particles : 35%|██████████████████▎ | 6549/18560 [00:59<04:09, 48.11it/s]orig_ncells = 12 Processing cells/particles : 35%|██████████████████▎ | 6558/18560 [01:20<12:04, 16.58it/s]orig_ncells = 228 Processing cells/particles : 36%|██████████████████▋ | 6673/18560 [01:55<25:18, 7.83it/s]orig_ncells = 0 Processing cells/particles : 37%|███████████████████ | 6789/18560 [02:14<27:10, 7.22it/s]orig_ncells = 1 Traceback (most recent call last): File "test_script_4.py", line 143, in n_photons, n_cells = pyxsim.make_photons(simput_basename + ".my_photons", my_sphere, redshift_project, area, simput_exp_time, thermal_model) File "/home/jovyan/.local/lib/python3.7/site-packages/pyxsim-3.0.0-py3.7-linux-x86_64.egg/pyxsim/photon_list.py", line 291, in make_photons chunk_data = source_model(chunk) File "/home/jovyan/.local/lib/python3.7/site-packages/pyxsim-3.0.0-py3.7-linux-x86_64.egg/pyxsim/source_models.py", line309, in call kT_sorted = kT[idxs] TypeError: 'float' object is not subscriptable

On Wed, Aug 11, 2021 at 10:14 AM John ZuHone @.***> wrote:

You need to download a new copy from source, you can't pip install from the one that pip installed originally.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/jzuhone/pyxsim/issues/35#issuecomment-896961399, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMEF72OWZT4V3LM2Y6Q4U7LT4KOYBANCNFSM5B5E5NJQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .

benopp99 commented 2 years ago

Looks like it's a special case when orig_ncells = 1. Baffles me.

Processing cells/particles : 13%|██████▋ | 1849/14063 [00:33<01:36, 126.44it/s]orig_ncells = 0 Processing cells/particles : 14%|███████▏ | 1939/14063 [00:50<04:31, 44.61it/s]orig_ncells = 92 Processing cells/particles : 14%|███████▏ | 1941/14063 [01:14<10:54, 18.51it/s]orig_ncells = 1213 Processing cells/particles : 23%|███████████▉ | 3244/14063 [01:57<05:50, 30.88it/s]orig_ncells = 61 Processing cells/particles : 23%|████████████ | 3249/14063 [02:23<10:03, 17.91it/s]orig_ncells = 2081 Processing cells/particles : 37%|███████████████████ | 5150/14063 [02:57<02:28, 59.92it/s]orig_ncells = 0 orig_ncells = 31 Processing cells/particles : 38%|███████████████████▉ | 5388/14063 [03:34<06:44, 21.43it/s]orig_ncells = 1 Traceback (most recent call last): File "test_script_4.py", line 143, in n_photons, n_cells = pyxsim.make_photons(simput_basename + ".my_photons", my_sphere, redshift_project, area, simput_exp_time, thermal_model) File "/home/jovyan/.local/lib/python3.7/site-packages/pyxsim-3.0.0-py3.7-linux-x86_64.egg/pyxsim/photon_list.py", line 291, in make_photons chunk_data = source_model(chunk) File "/home/jovyan/.local/lib/python3.7/site-packages/pyxsim-3.0.0-py3.7-linux-x86_64.egg/pyxsim/source_models.py", line309, in call kT_sorted = kT[idxs] TypeError: 'float' object is not subscriptable

On Wed, Aug 11, 2021 at 4:08 PM Benjamin Oppenheimer @.***> wrote:

First stoppage with printing orig_ncells. I don't see anything obvious.

Preparing spectrum table : 100%|█████████████████████████████████████████████████████████| 201/201 [00:41<00:00, 4.79it/s]

Processing cells/particles : 0%| | 0/18560 [00:00<?, ?it/s]orig_ncells = 0 Processing cells/particles : 0%| | 0/18560 [00:05<?, ?it/s]orig_ncells = 0 Processing cells/particles : 0%| | 0/18560 [00:10<?, ?it/s]orig_ncells = 0 Processing cells/particles : 0%| | 0/18560 [00:15<?, ?it/s]orig_ncells = 6549 Processing cells/particles : 35%|█████████████████▋ | 6549/18560 [00:49<00:06, 1838.06it/s]orig_ncells = 0 Processing cells/particles : 35%|██████████████████▎ | 6549/18560 [00:59<04:09, 48.11it/s]orig_ncells = 12 Processing cells/particles : 35%|██████████████████▎ | 6558/18560 [01:20<12:04, 16.58it/s]orig_ncells = 228 Processing cells/particles : 36%|██████████████████▋ | 6673/18560 [01:55<25:18, 7.83it/s]orig_ncells = 0 Processing cells/particles : 37%|███████████████████ | 6789/18560 [02:14<27:10, 7.22it/s]orig_ncells = 1 Traceback (most recent call last): File "test_script_4.py", line 143, in n_photons, n_cells = pyxsim.make_photons(simput_basename + ".my_photons", my_sphere, redshift_project, area, simput_exp_time, thermal_model) File "/home/jovyan/.local/lib/python3.7/site-packages/pyxsim-3.0.0-py3.7-linux-x86_64.egg/pyxsim/photon_list.py", line 291, in make_photons chunk_data = source_model(chunk) File "/home/jovyan/.local/lib/python3.7/site-packages/pyxsim-3.0.0-py3.7-linux-x86_64.egg/pyxsim/source_models.py", line309, in call kT_sorted = kT[idxs] TypeError: 'float' object is not subscriptable

On Wed, Aug 11, 2021 at 10:14 AM John ZuHone @.***> wrote:

You need to download a new copy from source, you can't pip install from the one that pip installed originally.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/jzuhone/pyxsim/issues/35#issuecomment-896961399, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMEF72OWZT4V3LM2Y6Q4U7LT4KOYBANCNFSM5B5E5NJQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .

jzuhone commented 2 years ago

Interesting, but it's exactly what I expected. Will fix.

benopp99 commented 2 years ago

Thanks! Have noted that there are cases where orig_ncells = 1 but it does go on.

-Ben

On Wed, Aug 11, 2021 at 6:44 PM John ZuHone @.***> wrote:

Interesting, but it's exactly what I expected. Will fix.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/jzuhone/pyxsim/issues/35#issuecomment-897260153, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMEF72OPUYJDUXPSNZK7C2TT4MKPVANCNFSM5B5E5NJQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .

jzuhone commented 2 years ago

Closed by #36.