From d793f4a0e17a99a3215063a8c8c664d9c3c439f7 Mon Sep 17 00:00:00 2001 From: Lukas Winkler Date: Thu, 30 Jun 2022 13:39:30 +0200 Subject: [PATCH 1/2] spectra fixes --- spectra_computation.py | 1 - spectra_plot.py | 5 +++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/spectra_computation.py b/spectra_computation.py index 29f8b80..0e2c9b5 100644 --- a/spectra_computation.py +++ b/spectra_computation.py @@ -80,7 +80,6 @@ def cross_run(waveforms: list, resolutions: list, Lbox: int, time: str): if __name__ == '__main__': input("are you sure you want to run this? This might need a large amount of memory") Lbox = 100 - k0 = 2 * 3.14159265358979323846264338327950 / Lbox resolutions = [128, 256, 512] waveforms = ["DB2", "DB4", "DB8", "shannon"] diff --git a/spectra_plot.py b/spectra_plot.py index 4cf07ed..bdceb93 100644 --- a/spectra_plot.py +++ b/spectra_plot.py @@ -123,9 +123,10 @@ def create_plot(mode): comp_p1 = comp_data["P1"] end_p1 /= comp_p1 - ax_ics.loglog(ics_k, ics_p1, color=colors[j]) - ax_end.loglog(end_k, end_p1, color=colors[j]) + ax_ics.semilogx(ics_k, ics_p1, color=colors[j]) + ax_end.semilogx(end_k, end_p1, color=colors[j]) for ax in [ax_ics, ax_end]: + ax.set_ylim(0.9,1.10) ax.axvline( k0 * resolution, color=colors[j], From 07cb2b8e9d6b77e6e0642e0c719b7e57fec530e7 Mon Sep 17 00:00:00 2001 From: Lukas Winkler Date: Thu, 30 Jun 2022 18:14:20 +0200 Subject: [PATCH 2/2] end of day --- compare_all.py | 28 ++++++++++++++++++---------- compare_halos.py | 19 ++++++++++--------- read_vr_files.py | 8 ++++---- sizes.py | 17 +++++++++++------ spectra_plot.py | 22 +++++++++------------- 5 files changed, 52 insertions(+), 42 deletions(-) diff --git a/compare_all.py b/compare_all.py index 00906c1..76f5271 100644 --- a/compare_all.py +++ b/compare_all.py @@ -1,15 +1,23 @@ from compare_halos import compare_halo_resolutions +something_failed = False for wf in ["DB2", "DB4", "DB8", "shannon"]: for res in [128, 256, 512]: for source_wf in [wf, "shannon"]: - compare_halo_resolutions( - ref_waveform=source_wf, - reference_resolution=128, - comp_waveform=wf, - comparison_resolution=res, - plot=False, - plot3d=False, - velo_halos=True, - single=False - ) + try: + compare_halo_resolutions( + ref_waveform=source_wf, + reference_resolution=128, + comp_waveform=wf, + comparison_resolution=res, + plot=False, + plot3d=False, + velo_halos=True, + single=False + ) + except: + something_failed = True + continue + +if something_failed: + print("\n\nAt least one comparison failed") diff --git a/compare_halos.py b/compare_halos.py index fdc1cc6..015362a 100644 --- a/compare_halos.py +++ b/compare_halos.py @@ -85,14 +85,12 @@ def compare_halo_resolutions( if ref_halo["cNFW"] < 0: print("NEGATIVE") print(ref_halo["cNFW"]) + raise ValueError() if len(halo_particle_ids) < 50: # TODO: decide on a lower size limit (and also apply it to comparison halo?) print(f"halo is too small ({len(halo_particle_ids)}") print("skipping") continue - if ref_halo_mass < 2: - print("skip weird, small mass halo") - continue print("LEN", len(halo_particle_ids), ref_halo.Mass_tot) offset_x, offset_y = ref_halo.X, ref_halo.Y # cumulative_mass_profile(particles_in_ref_halo, ref_halo, ref_meta, plot=plot) @@ -132,6 +130,9 @@ def compare_halo_resolutions( if len(nearby_halos) < 10: print(f"only {len(nearby_halos)} halos, expanding to 50xRvir") nearby_halos = set(df_comp_halo.loc[halo_distances < ref_halo.Rvir * 50].index.to_list()) + if len(nearby_halos) < 10: + print(f"only {len(nearby_halos)} halos, expanding to 150xRvir") + nearby_halos = set(df_comp_halo.loc[halo_distances < ref_halo.Rvir * 150].index.to_list()) if not nearby_halos: raise Exception("no halos are nearby") # TODO @@ -204,9 +205,6 @@ def compare_halo_resolutions( # particles_in_comp_halo: DataFrame = df_comp.loc[df_comp["FOFGroupIDs"] == halo] particle_ids_in_comp_halo = comp_halo_lookup[halo_id] mass_factor_limit = 5 - if comp_halo_masses[halo_id] < 2: - print("skip weird, small mass halo") - continue if not (1 / mass_factor_limit < (comp_halo_masses[halo_id] / ref_halo_mass) < mass_factor_limit): # print("mass not similar, skipping") @@ -216,7 +214,9 @@ def compare_halo_resolutions( halo_size = len(particle_ids_in_comp_halo) # df = particles_in_comp_halo.join(halo_particles, how="inner", rsuffix="ref") shared_particles = particle_ids_in_comp_halo.intersection(halo_particle_ids) - shared_size = len(shared_particles) + union_particles=particle_ids_in_comp_halo.union(halo_particle_ids) + + shared_size = len(shared_particles)/len(union_particles) # print(shared_size) if not shared_size: continue @@ -225,6 +225,7 @@ def compare_halo_resolutions( df = df_comp.loc[list(shared_particles)] if plot: color = f"C{i + 1}" + comp_halo: pd.Series = df_comp_halo.loc[halo_id] ax.scatter(apply_offset_to_list(df["X"], offset_x), apply_offset_to_list(df["Y"], offset_y), s=1, alpha=.3, c=color) @@ -253,7 +254,7 @@ def compare_halo_resolutions( np.array([ref_halo.X, ref_halo.Y, ref_halo.Z]) - np.array([comp_halo.X, comp_halo.Y, comp_halo.Z]) ) / ref_halo.Rvir halo_data["distance"] = distance - halo_data["match"] = best_halo_match / len(halo_particle_ids) + halo_data["match"] = best_halo_match compared_halos.append(halo_data) # exit() if plot: @@ -301,7 +302,7 @@ if __name__ == '__main__': reference_resolution=128, comparison_resolution=256, plot=False, - plot3d=False, + plot3d=True, plot_cic=False, velo_halos=True, single=False, diff --git a/read_vr_files.py b/read_vr_files.py index cf6ca04..5264a14 100644 --- a/read_vr_files.py +++ b/read_vr_files.py @@ -88,8 +88,8 @@ def read_velo_halos(directory: Path,veloname="vroutput"): Returns a dataframe containing all scalar properties of the halos (https://velociraptor-stf.readthedocs.io/en/latest/output.html), """ - group_catalog = h5py.File(directory / f"{veloname}.catalog_groups") - group_properties = h5py.File(directory / f"{veloname}.properties") + group_catalog = h5py.File(directory / f"{veloname}.catalog_groups.0") + group_properties = h5py.File(directory / f"{veloname}.properties.0") scalar_properties = {} for k, v in group_properties.items(): if not isinstance(v, Dataset): @@ -125,11 +125,11 @@ def read_velo_halo_particles( and two dictionaries mapping the halo IDs to sets of particle IDs """ df = read_velo_halos(directory) - particle_catalog = h5py.File(directory / "vroutput.catalog_particles") + particle_catalog = h5py.File(directory / "vroutput.catalog_particles.0") particle_ids = np.asarray(particle_catalog["Particle_IDs"]) particle_catalog_unbound = h5py.File( - directory / "vroutput.catalog_particles.unbound" + directory / "vroutput.catalog_particles.unbound.0" ) particle_ids_unbound = particle_catalog_unbound["Particle_IDs"][:] diff --git a/sizes.py b/sizes.py index 96b03a3..6928398 100644 --- a/sizes.py +++ b/sizes.py @@ -12,19 +12,24 @@ file = Path(argv[1]) df = pd.read_csv(file) -with pd.option_context('display.max_rows', None): - print(df[["ref_npart", "comp_npart", "ref_cNFW", "comp_cNFW"]]) # df = df.iloc fig: Figure = plt.figure() ax: Axes = fig.gca() # hist2d, log? -x_col = "ref_cNFW" -y_col = "comp_cNFW" +matches:pd.DataFrame=df["match"] +# print(matches) +# exit() +print(matches.median()) +print(matches.std()) +exit() -# x_col = "ref_Mvir" -# y_col = "comp_Mvir" +# x_col = "ref_cNFW" +# y_col = "comp_cNFW" +# +x_col = "ref_Mvir" +y_col = "comp_Mvir" min_x = min([min(df[x_col]), min(df[y_col])]) max_x = max([max(df[x_col]), max(df[y_col])]) diff --git a/spectra_plot.py b/spectra_plot.py index bdceb93..e03b9b4 100644 --- a/spectra_plot.py +++ b/spectra_plot.py @@ -29,7 +29,8 @@ columns = [ ] # linestyles = ["solid", "dashed", "dotted"] -colors = ["C1", "C2", "C3", "C4"] +colors=[f"C{i}" for i in range(10)] +# colors = ["C1", "C2", "C3", "C4"] def spectra_data( @@ -96,10 +97,10 @@ def create_plot(mode): verticalalignment="top", transform=ax.transAxes, ) - for res in [128]: + for j, res in enumerate(resolutions[:-1] if mode == "cross" else resolutions): ax.axvline( k0 * res, - color="gray", + color=colors[j], linestyle="dashed", label=f"{res}", ) @@ -126,13 +127,8 @@ def create_plot(mode): ax_ics.semilogx(ics_k, ics_p1, color=colors[j]) ax_end.semilogx(end_k, end_p1, color=colors[j]) for ax in [ax_ics, ax_end]: - ax.set_ylim(0.9,1.10) - ax.axvline( - k0 * resolution, - color=colors[j], - linestyle="dashed", - label=f"{resolution}", - ) + ax.set_ylim(0.9, 1.10) + # fig.suptitle(f"Power Spectra {time}") #Not needed for paper # fig.tight_layout() @@ -145,15 +141,15 @@ def create_plot(mode): ics_k = ics_data["k [Mpc]"] ics_pcross = ics_data["Pcross"] - ax_ics.semilogx(ics_k, ics_pcross, color=colors[j], label=f'{res1} vs {res2}') + ax_ics.semilogx(ics_k, ics_pcross, color=colors[j+3], label=f'{res1} vs {res2}') end_data = spectra_data(waveform, res1, res2, Lbox, 'end') end_k = end_data["k [Mpc]"] end_pcross = end_data["Pcross"] - ax_end.semilogx(end_k, end_pcross, color=colors[j], label=f'{res1} vs {res2}') + ax_end.semilogx(end_k, end_pcross, color=colors[j+3], label=f'{res1} vs {res2}') - ax_end.set_xlim(right=k0 * 256) + ax_end.set_xlim(right=k0 * resolutions[-1]) ax_end.set_ylim(0.8, 1.02) if bottom_row: ax_end.legend()